LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
Public Types | Public Member Functions | Protected Member Functions | Protected Attributes | Static Protected Attributes | List of all members
lsst::ip::diffim::StaticKernelSolution< InputT > Class Template Reference

#include <KernelSolution.h>

Inheritance diagram for lsst::ip::diffim::StaticKernelSolution< InputT >:
lsst::ip::diffim::KernelSolution lsst::ip::diffim::MaskedKernelSolution< InputT > lsst::ip::diffim::RegularizedKernelSolution< InputT >

Public Types

typedef std::shared_ptr< StaticKernelSolution< InputT > > Ptr
 
enum  KernelSolvedBy {
  NONE = 0 , CHOLESKY_LDLT = 1 , CHOLESKY_LLT = 2 , LU = 3 ,
  EIGENVECTOR = 4
}
 
enum  ConditionNumberType { EIGENVALUE = 0 , SVD = 1 }
 
typedef lsst::afw::math::Kernel::Pixel PixelT
 
typedef lsst::afw::image::Image< lsst::afw::math::Kernel::PixelImageT
 

Public Member Functions

 StaticKernelSolution (lsst::afw::math::KernelList const &basisList, bool fitForBackground)
 
virtual ~StaticKernelSolution ()
 
void solve ()
 
virtual void build (lsst::afw::image::Image< InputT > const &templateImage, lsst::afw::image::Image< InputT > const &scienceImage, lsst::afw::image::Image< lsst::afw::image::VariancePixel > const &varianceEstimate)
 
virtual std::shared_ptr< lsst::afw::math::KernelgetKernel ()
 
virtual std::shared_ptr< lsst::afw::image::Image< lsst::afw::math::Kernel::Pixel > > makeKernelImage ()
 
virtual double getBackground ()
 
virtual double getKsum ()
 
virtual std::pair< std::shared_ptr< lsst::afw::math::Kernel >, double > getSolutionPair ()
 
virtual void solve (Eigen::MatrixXd const &mMat, Eigen::VectorXd const &bVec)
 
KernelSolvedBy getSolvedBy ()
 
virtual double getConditionNumber (ConditionNumberType conditionType)
 
virtual double getConditionNumber (Eigen::MatrixXd const &mMat, ConditionNumberType conditionType)
 
Eigen::MatrixXd const & getM ()
 
Eigen::VectorXd const & getB ()
 
void printM ()
 
void printB ()
 
void printA ()
 
int getId () const
 

Protected Member Functions

void _setKernel ()
 Set kernel after solution. More...
 
void _setKernelUncertainty ()
 Not implemented. More...
 

Protected Attributes

Eigen::MatrixXd _cMat
 K_i x R. More...
 
Eigen::VectorXd _iVec
 Vectorized I. More...
 
Eigen::VectorXd _ivVec
 Inverse variance. More...
 
std::shared_ptr< lsst::afw::math::Kernel_kernel
 Derived single-object convolution kernel. More...
 
double _background
 Derived differential background estimate. More...
 
double _kSum
 Derived kernel sum. More...
 
int _id
 Unique ID for object. More...
 
Eigen::MatrixXd _mMat
 Derived least squares M matrix. More...
 
Eigen::VectorXd _bVec
 Derived least squares B vector. More...
 
Eigen::VectorXd _aVec
 Derived least squares solution matrix. More...
 
KernelSolvedBy _solvedBy
 Type of algorithm used to make solution. More...
 
bool _fitForBackground
 Background terms included in fit. More...
 

Static Protected Attributes

static int _SolutionId = 0
 Unique identifier for solution. More...
 

Detailed Description

template<typename InputT>
class lsst::ip::diffim::StaticKernelSolution< InputT >

Definition at line 83 of file KernelSolution.h.

Member Typedef Documentation

◆ ImageT

Definition at line 35 of file KernelSolution.h.

◆ PixelT

Definition at line 34 of file KernelSolution.h.

◆ Ptr

Definition at line 85 of file KernelSolution.h.

Member Enumeration Documentation

◆ ConditionNumberType

Enumerator
EIGENVALUE 
SVD 

Definition at line 45 of file KernelSolution.h.

◆ KernelSolvedBy

Enumerator
NONE 
CHOLESKY_LDLT 
CHOLESKY_LLT 
LU 
EIGENVECTOR 

Definition at line 37 of file KernelSolution.h.

Constructor & Destructor Documentation

◆ StaticKernelSolution()

template<typename InputT >
lsst::ip::diffim::StaticKernelSolution< InputT >::StaticKernelSolution ( lsst::afw::math::KernelList const &  basisList,
bool  fitForBackground 
)

Definition at line 195 of file KernelSolution.cc.

199  :
200  KernelSolution(fitForBackground),
201  _cMat(),
202  _iVec(),
203  _ivVec(),
204  _kernel(),
205  _background(0.0),
206  _kSum(0.0)
207  {
208  std::vector<double> kValues(basisList.size());
210  new afwMath::LinearCombinationKernel(basisList, kValues)
211  );
212  };
A kernel that is a linear combination of fixed basis kernels.
Definition: Kernel.h:704
double _background
Derived differential background estimate.
Eigen::VectorXd _iVec
Vectorized I.
std::shared_ptr< lsst::afw::math::Kernel > _kernel
Derived single-object convolution kernel.
Eigen::VectorXd _ivVec
Inverse variance.

◆ ~StaticKernelSolution()

template<typename InputT >
virtual lsst::ip::diffim::StaticKernelSolution< InputT >::~StaticKernelSolution ( )
inlinevirtual

Definition at line 89 of file KernelSolution.h.

89 {};

Member Function Documentation

◆ _setKernel()

template<typename InputT >
void lsst::ip::diffim::StaticKernelSolution< InputT >::_setKernel
protected

Set kernel after solution.

Definition at line 424 of file KernelSolution.cc.

424  {
426  throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot make solution");
427  }
428 
429  unsigned int const nParameters = _aVec.size();
430  unsigned int const nBackgroundParameters = _fitForBackground ? 1 : 0;
431  unsigned int const nKernelParameters =
432  std::dynamic_pointer_cast<afwMath::LinearCombinationKernel>(_kernel)->getKernelList().size();
433  if (nParameters != (nKernelParameters + nBackgroundParameters))
434  throw LSST_EXCEPT(pexExcept::Exception, "Mismatched sizes in kernel solution");
435 
436  /* Fill in the kernel results */
437  std::vector<double> kValues(nKernelParameters);
438  for (unsigned int idx = 0; idx < nKernelParameters; idx++) {
439  if (std::isnan(_aVec(idx))) {
441  str(boost::format("Unable to determine kernel solution %d (nan)") % idx));
442  }
443  kValues[idx] = _aVec(idx);
444  }
445  _kernel->setKernelParameters(kValues);
446 
449  );
450  _kSum = _kernel->computeImage(*image, false);
451 
452  if (_fitForBackground) {
453  if (std::isnan(_aVec(nParameters-1))) {
455  str(boost::format("Unable to determine background solution %d (nan)") %
456  (nParameters-1)));
457  }
458  _background = _aVec(nParameters-1);
459  }
460  }
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
afw::table::Key< afw::table::Array< ImagePixelT > > image
lsst::geom::Extent2I const getDimensions() const
Return the Kernel's dimensions (width, height)
Definition: Kernel.h:212
void setKernelParameters(std::vector< double > const &params)
Set the kernel parameters of a spatially invariant kernel.
Definition: Kernel.h:341
double computeImage(lsst::afw::image::Image< Pixel > &image, bool doNormalize, double x=0.0, double y=0.0) const
Compute an image (pixellized representation of the kernel) in place.
Definition: Kernel.cc:76
bool _fitForBackground
Background terms included in fit.
Eigen::VectorXd _aVec
Derived least squares solution matrix.
KernelSolvedBy _solvedBy
Type of algorithm used to make solution.
lsst::afw::image::Image< lsst::afw::math::Kernel::Pixel > ImageT
Provides consistent interface for LSST exceptions.
Definition: Exception.h:107
T isnan(T... args)
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174

◆ _setKernelUncertainty()

template<typename InputT >
void lsst::ip::diffim::StaticKernelSolution< InputT >::_setKernelUncertainty
protected

Not implemented.

Definition at line 464 of file KernelSolution.cc.

464  {
465  throw LSST_EXCEPT(pexExcept::Exception, "Uncertainty calculation not supported");
466 
467  /* Estimate of parameter uncertainties comes from the inverse of the
468  * covariance matrix (noise spectrum).
469  * N.R. 15.4.8 to 15.4.15
470  *
471  * Since this is a linear problem no need to use Fisher matrix
472  * N.R. 15.5.8
473  *
474  * Although I might be able to take advantage of the solution above.
475  * Since this now works and is not the rate limiting step, keep as-is for DC3a.
476  *
477  * Use Cholesky decomposition again.
478  * Cholkesy:
479  * Cov = L L^t
480  * Cov^(-1) = (L L^t)^(-1)
481  * = (L^T)^-1 L^(-1)
482  *
483  * Code would be:
484  *
485  * Eigen::MatrixXd cov = _mMat.transpose() * _mMat;
486  * Eigen::LLT<Eigen::MatrixXd> llt = cov.llt();
487  * Eigen::MatrixXd error2 = llt.matrixL().transpose().inverse()*llt.matrixL().inverse();
488  */
489  }

◆ build()

template<typename InputT >
void lsst::ip::diffim::StaticKernelSolution< InputT >::build ( lsst::afw::image::Image< InputT > const &  templateImage,
lsst::afw::image::Image< InputT > const &  scienceImage,
lsst::afw::image::Image< lsst::afw::image::VariancePixel > const &  varianceEstimate 
)
virtual

Definition at line 261 of file KernelSolution.cc.

265  {
266 
267  afwMath::Statistics varStats = afwMath::makeStatistics(varianceEstimate, afwMath::MIN);
268  if (varStats.getValue(afwMath::MIN) < 0.0) {
270  "Error: variance less than 0.0");
271  }
272  if (varStats.getValue(afwMath::MIN) == 0.0) {
274  "Error: variance equals 0.0, cannot inverse variance weight");
275  }
276 
277  lsst::afw::math::KernelList basisList =
278  std::dynamic_pointer_cast<afwMath::LinearCombinationKernel>(_kernel)->getKernelList();
279 
280  unsigned int const nKernelParameters = basisList.size();
281  unsigned int const nBackgroundParameters = _fitForBackground ? 1 : 0;
282  unsigned int const nParameters = nKernelParameters + nBackgroundParameters;
283 
284  std::vector<std::shared_ptr<afwMath::Kernel> >::const_iterator kiter = basisList.begin();
285 
286  /* Ignore buffers around edge of convolved images :
287  *
288  * If the kernel has width 5, it has center pixel 2. The first good pixel
289  * is the (5-2)=3rd pixel, which is array index 2, and ends up being the
290  * index of the central pixel.
291  *
292  * You also have a buffer of unusable pixels on the other side, numbered
293  * width-center-1. The last good usable pixel is N-width+center+1.
294  *
295  * Example : the kernel is width = 5, center = 2
296  *
297  * ---|---|-c-|---|---|
298  *
299  * the image is width = N
300  * convolve this with the kernel, and you get
301  *
302  * |-x-|-x-|-g-|---|---| ... |---|---|-g-|-x-|-x-|
303  *
304  * g = first/last good pixel
305  * x = bad
306  *
307  * the first good pixel is the array index that has the value "center", 2
308  * the last good pixel has array index N-(5-2)+1
309  * eg. if N = 100, you want to use up to index 97
310  * 100-3+1 = 98, and the loops use i < 98, meaning the last
311  * index you address is 97.
312  */
313 
314  /* NOTE - we are accessing particular elements of Eigen arrays using
315  these coordinates, therefore they need to be in LOCAL coordinates.
316  This was written before ndarray unification.
317  */
318  geom::Box2I goodBBox = (*kiter)->shrinkBBox(templateImage.getBBox(afwImage::LOCAL));
319  unsigned int const startCol = goodBBox.getMinX();
320  unsigned int const startRow = goodBBox.getMinY();
321  // endCol/Row is one past the index of the last good col/row
322  unsigned int endCol = goodBBox.getMaxX() + 1;
323  unsigned int endRow = goodBBox.getMaxY() + 1;
324 
325  boost::timer t;
326  t.restart();
327 
328  /* Eigen representation of input images; only the pixels that are unconvolved in cimage below */
329  Eigen::MatrixXd eigenTemplate = imageToEigenMatrix(templateImage).block(startRow,
330  startCol,
331  endRow-startRow,
332  endCol-startCol);
333  Eigen::MatrixXd eigenScience = imageToEigenMatrix(scienceImage).block(startRow,
334  startCol,
335  endRow-startRow,
336  endCol-startCol);
337  Eigen::MatrixXd eigeniVariance = imageToEigenMatrix(varianceEstimate).block(
338  startRow, startCol, endRow-startRow, endCol-startCol
339  ).array().inverse().matrix();
340 
341  /* Resize into 1-D for later usage */
342  eigenTemplate.resize(eigenTemplate.rows()*eigenTemplate.cols(), 1);
343  eigenScience.resize(eigenScience.rows()*eigenScience.cols(), 1);
344  eigeniVariance.resize(eigeniVariance.rows()*eigeniVariance.cols(), 1);
345 
346  /* Holds image convolved with basis function */
347  afwImage::Image<PixelT> cimage(templateImage.getDimensions());
348 
349  /* Holds eigen representation of image convolved with all basis functions */
350  std::vector<Eigen::MatrixXd> convolvedEigenList(nKernelParameters);
351 
352  /* Iterators over convolved image list and basis list */
353  typename std::vector<Eigen::MatrixXd>::iterator eiter = convolvedEigenList.begin();
354  /* Create C_i in the formalism of Alard & Lupton */
355  afwMath::ConvolutionControl convolutionControl;
356  convolutionControl.setDoNormalize(false);
357  for (kiter = basisList.begin(); kiter != basisList.end(); ++kiter, ++eiter) {
358  afwMath::convolve(cimage, templateImage, **kiter, convolutionControl); /* cimage stores convolved image */
359 
360  Eigen::MatrixXd cMat = imageToEigenMatrix(cimage).block(startRow,
361  startCol,
362  endRow-startRow,
363  endCol-startCol);
364  cMat.resize(cMat.size(), 1);
365  *eiter = cMat;
366 
367  }
368 
369  double time = t.elapsed();
370  LOGL_DEBUG("TRACE3.ip.diffim.StaticKernelSolution.build",
371  "Total compute time to do basis convolutions : %.2f s", time);
372  t.restart();
373 
374  /*
375  Load matrix with all values from convolvedEigenList : all images
376  (eigeniVariance, convolvedEigenList) must be the same size
377  */
378  Eigen::MatrixXd cMat(eigenTemplate.col(0).size(), nParameters);
379  typename std::vector<Eigen::MatrixXd>::iterator eiterj = convolvedEigenList.begin();
380  typename std::vector<Eigen::MatrixXd>::iterator eiterE = convolvedEigenList.end();
381  for (unsigned int kidxj = 0; eiterj != eiterE; eiterj++, kidxj++) {
382  cMat.col(kidxj) = eiterj->col(0);
383  }
384  /* Treat the last "image" as all 1's to do the background calculation. */
385  if (_fitForBackground)
386  cMat.col(nParameters-1).fill(1.);
387 
388  _cMat = cMat;
389  _ivVec = eigeniVariance.col(0);
390  _iVec = eigenScience.col(0);
391 
392  /* Make these outside of solve() so I can check condition number */
393  _mMat = _cMat.transpose() * (_ivVec.asDiagonal() * _cMat);
394  _bVec = _cMat.transpose() * (_ivVec.asDiagonal() * _iVec);
395  }
#define LOGL_DEBUG(logger, message...)
Log a debug-level message using a varargs/printf style interface.
Definition: Log.h:515
T begin(T... args)
lsst::geom::Box2I getBBox(ImageOrigin origin=PARENT) const
Definition: ImageBase.h:445
lsst::geom::Extent2I getDimensions() const
Return the image's size; useful for passing to constructors.
Definition: ImageBase.h:356
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:51
Parameters to control convolution.
Definition: ConvolveImage.h:50
void setDoNormalize(bool doNormalize)
Definition: ConvolveImage.h:66
A class to evaluate image statistics.
Definition: Statistics.h:220
double getValue(Property const prop=NOTHING) const
Return the value of the desired property (if specified in the constructor)
Definition: Statistics.cc:1047
An integer coordinate rectangle.
Definition: Box.h:55
int getMinY() const noexcept
Definition: Box.h:158
int getMinX() const noexcept
Definition: Box.h:157
int getMaxX() const noexcept
Definition: Box.h:161
int getMaxY() const noexcept
Definition: Box.h:162
Eigen::VectorXd _bVec
Derived least squares B vector.
Eigen::MatrixXd _mMat
Derived least squares M matrix.
T end(T... args)
Statistics makeStatistics(lsst::afw::image::Image< Pixel > const &img, lsst::afw::image::Mask< image::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl=StatisticsControl())
Handle a watered-down front-end to the constructor (no variance)
Definition: Statistics.h:359
void convolve(OutImageT &convolvedImage, InImageT const &inImage, KernelT const &kernel, ConvolutionControl const &convolutionControl=ConvolutionControl())
Convolve an Image or MaskedImage with a Kernel, setting pixels of an existing output image.
@ MIN
estimate sample minimum
Definition: Statistics.h:75
Eigen::MatrixXd imageToEigenMatrix(lsst::afw::image::Image< PixelT > const &img)
Turns a 2-d Image into a 2-d Eigen Matrix.
T size(T... args)
T time(T... args)

◆ getB()

Eigen::VectorXd const& lsst::ip::diffim::KernelSolution::getB ( )
inlineinherited

Definition at line 65 of file KernelSolution.h.

65 {return _bVec;}

◆ getBackground()

template<typename InputT >
double lsst::ip::diffim::StaticKernelSolution< InputT >::getBackground
virtual

Definition at line 235 of file KernelSolution.cc.

235  {
237  throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return background");
238  }
239  return _background;
240  }

◆ getConditionNumber() [1/2]

double lsst::ip::diffim::KernelSolution::getConditionNumber ( ConditionNumberType  conditionType)
virtualinherited

Definition at line 94 of file KernelSolution.cc.

94  {
95  return getConditionNumber(_mMat, conditionType);
96  }
virtual double getConditionNumber(ConditionNumberType conditionType)

◆ getConditionNumber() [2/2]

double lsst::ip::diffim::KernelSolution::getConditionNumber ( Eigen::MatrixXd const &  mMat,
ConditionNumberType  conditionType 
)
virtualinherited

Definition at line 98 of file KernelSolution.cc.

99  {
100  switch (conditionType) {
101  case EIGENVALUE:
102  {
103  Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eVecValues(mMat);
104  Eigen::VectorXd eValues = eVecValues.eigenvalues();
105  double eMax = eValues.maxCoeff();
106  double eMin = eValues.minCoeff();
107  LOGL_DEBUG("TRACE3.ip.diffim.KernelSolution.getConditionNumber",
108  "EIGENVALUE eMax / eMin = %.3e", eMax / eMin);
109  return (eMax / eMin);
110  break;
111  }
112  case SVD:
113  {
114  Eigen::VectorXd sValues = mMat.jacobiSvd().singularValues();
115  double sMax = sValues.maxCoeff();
116  double sMin = sValues.minCoeff();
117  LOGL_DEBUG("TRACE3.ip.diffim.KernelSolution.getConditionNumber",
118  "SVD eMax / eMin = %.3e", sMax / sMin);
119  return (sMax / sMin);
120  break;
121  }
122  default:
123  {
125  "Undefined ConditionNumberType : only EIGENVALUE, SVD allowed.");
126  break;
127  }
128  }
129  }
Reports invalid arguments.
Definition: Runtime.h:66

◆ getId()

int lsst::ip::diffim::KernelSolution::getId ( ) const
inlineinherited

Definition at line 69 of file KernelSolution.h.

69 { return _id; }
int _id
Unique ID for object.

◆ getKernel()

template<typename InputT >
std::shared_ptr< lsst::afw::math::Kernel > lsst::ip::diffim::StaticKernelSolution< InputT >::getKernel
virtual

Definition at line 215 of file KernelSolution.cc.

215  {
217  throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return solution");
218  }
219  return _kernel;
220  }

◆ getKsum()

template<typename InputT >
double lsst::ip::diffim::StaticKernelSolution< InputT >::getKsum
virtual

Definition at line 243 of file KernelSolution.cc.

243  {
245  throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return ksum");
246  }
247  return _kSum;
248  }

◆ getM()

Eigen::MatrixXd const& lsst::ip::diffim::KernelSolution::getM ( )
inlineinherited

Definition at line 64 of file KernelSolution.h.

64 {return _mMat;}

◆ getSolutionPair()

template<typename InputT >
std::pair< std::shared_ptr< lsst::afw::math::Kernel >, double > lsst::ip::diffim::StaticKernelSolution< InputT >::getSolutionPair
virtual

Definition at line 252 of file KernelSolution.cc.

252  {
254  throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return solution");
255  }
256 
258  }
T make_pair(T... args)

◆ getSolvedBy()

KernelSolvedBy lsst::ip::diffim::KernelSolution::getSolvedBy ( )
inlineinherited

Definition at line 60 of file KernelSolution.h.

60 {return _solvedBy;}

◆ makeKernelImage()

Definition at line 223 of file KernelSolution.cc.

223  {
225  throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return image");
226  }
229  );
230  (void)_kernel->computeImage(*image, false);
231  return image;
232  }

◆ printA()

void lsst::ip::diffim::KernelSolution::printA ( )
inlineinherited

Definition at line 68 of file KernelSolution.h.

68 {std::cout << _aVec << std::endl;}
T endl(T... args)

◆ printB()

void lsst::ip::diffim::KernelSolution::printB ( )
inlineinherited

Definition at line 67 of file KernelSolution.h.

67 {std::cout << _bVec << std::endl;}

◆ printM()

void lsst::ip::diffim::KernelSolution::printM ( )
inlineinherited

Definition at line 66 of file KernelSolution.h.

66 {std::cout << _mMat << std::endl;}

◆ solve() [1/2]

template<typename InputT >
void lsst::ip::diffim::StaticKernelSolution< InputT >::solve
virtual

Reimplemented from lsst::ip::diffim::KernelSolution.

Reimplemented in lsst::ip::diffim::RegularizedKernelSolution< InputT >.

Definition at line 398 of file KernelSolution.cc.

398  {
399  LOGL_DEBUG("TRACE3.ip.diffim.StaticKernelSolution.solve",
400  "mMat is %d x %d; bVec is %d; cMat is %d x %d; vVec is %d; iVec is %d",
401  _mMat.rows(), _mMat.cols(), _bVec.size(),
402  _cMat.rows(), _cMat.cols(), _ivVec.size(), _iVec.size());
403 
404  if (DEBUG_MATRIX) {
405  std::cout << "C" << std::endl;
406  std::cout << _cMat << std::endl;
407  std::cout << "iV" << std::endl;
408  std::cout << _ivVec << std::endl;
409  std::cout << "I" << std::endl;
410  std::cout << _iVec << std::endl;
411  }
412 
413  try {
415  } catch (pexExcept::Exception &e) {
416  LSST_EXCEPT_ADD(e, "Unable to solve static kernel matrix");
417  throw e;
418  }
419  /* Turn matrices into _kernel and _background */
420  _setKernel();
421  }
#define LSST_EXCEPT_ADD(e, m)
Add the current location and a message to an existing exception before rethrowing it.
Definition: Exception.h:54
void _setKernel()
Set kernel after solution.
#define DEBUG_MATRIX

◆ solve() [2/2]

void lsst::ip::diffim::KernelSolution::solve ( Eigen::MatrixXd const &  mMat,
Eigen::VectorXd const &  bVec 
)
virtualinherited

Definition at line 131 of file KernelSolution.cc.

132  {
133 
134  if (DEBUG_MATRIX) {
135  std::cout << "M " << std::endl;
136  std::cout << mMat << std::endl;
137  std::cout << "B " << std::endl;
138  std::cout << bVec << std::endl;
139  }
140 
141  Eigen::VectorXd aVec = Eigen::VectorXd::Zero(bVec.size());
142 
143  boost::timer t;
144  t.restart();
145 
146  LOGL_DEBUG("TRACE2.ip.diffim.KernelSolution.solve",
147  "Solving for kernel");
148  _solvedBy = LU;
149  Eigen::FullPivLU<Eigen::MatrixXd> lu(mMat);
150  if (lu.isInvertible()) {
151  aVec = lu.solve(bVec);
152  } else {
153  LOGL_DEBUG("TRACE3.ip.diffim.KernelSolution.solve",
154  "Unable to determine kernel via LU");
155  /* LAST RESORT */
156  try {
157 
159  Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eVecValues(mMat);
160  Eigen::MatrixXd const& rMat = eVecValues.eigenvectors();
161  Eigen::VectorXd eValues = eVecValues.eigenvalues();
162 
163  for (int i = 0; i != eValues.rows(); ++i) {
164  if (eValues(i) != 0.0) {
165  eValues(i) = 1.0/eValues(i);
166  }
167  }
168 
169  aVec = rMat * eValues.asDiagonal() * rMat.transpose() * bVec;
170  } catch (pexExcept::Exception& e) {
171 
172  _solvedBy = NONE;
173  LOGL_DEBUG("TRACE3.ip.diffim.KernelSolution.solve",
174  "Unable to determine kernel via eigen-values");
175 
176  throw LSST_EXCEPT(pexExcept::Exception, "Unable to determine kernel solution");
177  }
178  }
179 
180  double time = t.elapsed();
181  LOGL_DEBUG("TRACE3.ip.diffim.KernelSolution.solve",
182  "Compute time for matrix math : %.2f s", time);
183 
184  if (DEBUG_MATRIX) {
185  std::cout << "A " << std::endl;
186  std::cout << aVec << std::endl;
187  }
188 
189  _aVec = aVec;
190  }

Member Data Documentation

◆ _aVec

Eigen::VectorXd lsst::ip::diffim::KernelSolution::_aVec
protectedinherited

Derived least squares solution matrix.

Definition at line 75 of file KernelSolution.h.

◆ _background

template<typename InputT >
double lsst::ip::diffim::StaticKernelSolution< InputT >::_background
protected

Derived differential background estimate.

Definition at line 110 of file KernelSolution.h.

◆ _bVec

Eigen::VectorXd lsst::ip::diffim::KernelSolution::_bVec
protectedinherited

Derived least squares B vector.

Definition at line 74 of file KernelSolution.h.

◆ _cMat

template<typename InputT >
Eigen::MatrixXd lsst::ip::diffim::StaticKernelSolution< InputT >::_cMat
protected

K_i x R.

Definition at line 105 of file KernelSolution.h.

◆ _fitForBackground

bool lsst::ip::diffim::KernelSolution::_fitForBackground
protectedinherited

Background terms included in fit.

Definition at line 77 of file KernelSolution.h.

◆ _id

int lsst::ip::diffim::KernelSolution::_id
protectedinherited

Unique ID for object.

Definition at line 72 of file KernelSolution.h.

◆ _iVec

template<typename InputT >
Eigen::VectorXd lsst::ip::diffim::StaticKernelSolution< InputT >::_iVec
protected

Vectorized I.

Definition at line 106 of file KernelSolution.h.

◆ _ivVec

template<typename InputT >
Eigen::VectorXd lsst::ip::diffim::StaticKernelSolution< InputT >::_ivVec
protected

Inverse variance.

Definition at line 107 of file KernelSolution.h.

◆ _kernel

template<typename InputT >
std::shared_ptr<lsst::afw::math::Kernel> lsst::ip::diffim::StaticKernelSolution< InputT >::_kernel
protected

Derived single-object convolution kernel.

Definition at line 109 of file KernelSolution.h.

◆ _kSum

template<typename InputT >
double lsst::ip::diffim::StaticKernelSolution< InputT >::_kSum
protected

Derived kernel sum.

Definition at line 111 of file KernelSolution.h.

◆ _mMat

Eigen::MatrixXd lsst::ip::diffim::KernelSolution::_mMat
protectedinherited

Derived least squares M matrix.

Definition at line 73 of file KernelSolution.h.

◆ _SolutionId

int lsst::ip::diffim::KernelSolution::_SolutionId = 0
staticprotectedinherited

Unique identifier for solution.

Definition at line 78 of file KernelSolution.h.

◆ _solvedBy

KernelSolvedBy lsst::ip::diffim::KernelSolution::_solvedBy
protectedinherited

Type of algorithm used to make solution.

Definition at line 76 of file KernelSolution.h.


The documentation for this class was generated from the following files: