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::RegularizedKernelSolution< InputT > Class Template Reference

#include <KernelSolution.h>

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

Public Types

typedef std::shared_ptr< RegularizedKernelSolution< 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

 RegularizedKernelSolution (lsst::afw::math::KernelList const &basisList, bool fitForBackground, Eigen::MatrixXd const &hMat, lsst::daf::base::PropertySet const &ps)
 
virtual ~RegularizedKernelSolution ()
 
void solve ()
 
double getLambda ()
 
double estimateRisk (double maxCond)
 
Eigen::MatrixXd getM (bool includeHmat=true)
 
virtual void solve (Eigen::MatrixXd const &mMat, Eigen::VectorXd const &bVec)
 
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 ()
 
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::RegularizedKernelSolution< InputT >

Definition at line 148 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 150 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

◆ RegularizedKernelSolution()

template<typename InputT >
lsst::ip::diffim::RegularizedKernelSolution< InputT >::RegularizedKernelSolution ( lsst::afw::math::KernelList const &  basisList,
bool  fitForBackground,
Eigen::MatrixXd const &  hMat,
lsst::daf::base::PropertySet const &  ps 
)

Definition at line 1067 of file KernelSolution.cc.

1073  :
1074  StaticKernelSolution<InputT>(basisList, fitForBackground),
1075  _hMat(hMat),
1076  _ps(ps.deepCopy())
1077  {};

◆ ~RegularizedKernelSolution()

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

Definition at line 157 of file KernelSolution.h.

157 {};

Member Function Documentation

◆ _setKernel()

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

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
double _background
Derived differential background estimate.
std::shared_ptr< lsst::afw::math::Kernel > _kernel
Derived single-object convolution kernel.
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
protectedinherited

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 
)
virtualinherited

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.
Eigen::VectorXd _iVec
Vectorized I.
Eigen::VectorXd _ivVec
Inverse variance.
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)

◆ estimateRisk()

template<typename InputT >
double lsst::ip::diffim::RegularizedKernelSolution< InputT >::estimateRisk ( double  maxCond)

Definition at line 1080 of file KernelSolution.cc.

1080  {
1081  Eigen::MatrixXd vMat = this->_cMat.jacobiSvd().matrixV();
1082  Eigen::MatrixXd vMatvMatT = vMat * vMat.transpose();
1083 
1084  /* Find pseudo inverse of mMat, which may be ill conditioned */
1085  Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eVecValues(this->_mMat);
1086  Eigen::MatrixXd const& rMat = eVecValues.eigenvectors();
1087  Eigen::VectorXd eValues = eVecValues.eigenvalues();
1088  double eMax = eValues.maxCoeff();
1089  for (int i = 0; i != eValues.rows(); ++i) {
1090  if (eValues(i) == 0.0) {
1091  eValues(i) = 0.0;
1092  }
1093  else if ((eMax / eValues(i)) > maxCond) {
1094  LOGL_DEBUG("TRACE3.ip.diffim.RegularizedKernelSolution.estimateRisk",
1095  "Truncating eValue %d; %.5e / %.5e = %.5e vs. %.5e",
1096  i, eMax, eValues(i), eMax / eValues(i), maxCond);
1097  eValues(i) = 0.0;
1098  }
1099  else {
1100  eValues(i) = 1.0 / eValues(i);
1101  }
1102  }
1103  Eigen::MatrixXd mInv = rMat * eValues.asDiagonal() * rMat.transpose();
1104 
1105  std::vector<double> lambdas = _createLambdaSteps();
1106  std::vector<double> risks;
1107  for (unsigned int i = 0; i < lambdas.size(); i++) {
1108  double l = lambdas[i];
1109  Eigen::MatrixXd mLambda = this->_mMat + l * _hMat;
1110 
1111  try {
1112  KernelSolution::solve(mLambda, this->_bVec);
1113  } catch (pexExcept::Exception &e) {
1114  LSST_EXCEPT_ADD(e, "Unable to solve regularized kernel matrix");
1115  throw e;
1116  }
1117  Eigen::VectorXd term1 = (this->_aVec.transpose() * vMatvMatT * this->_aVec);
1118  if (term1.size() != 1)
1119  throw LSST_EXCEPT(pexExcept::Exception, "Matrix size mismatch");
1120 
1121  double term2a = (vMatvMatT * mLambda.inverse()).trace();
1122 
1123  Eigen::VectorXd term2b = (this->_aVec.transpose() * (mInv * this->_bVec));
1124  if (term2b.size() != 1)
1125  throw LSST_EXCEPT(pexExcept::Exception, "Matrix size mismatch");
1126 
1127  double risk = term1(0) + 2 * (term2a - term2b(0));
1128  LOGL_DEBUG("TRACE4.ip.diffim.RegularizedKernelSolution.estimateRisk",
1129  "Lambda = %.3f, Risk = %.5e",
1130  l, risk);
1131  LOGL_DEBUG("TRACE5.ip.diffim.RegularizedKernelSolution.estimateRisk",
1132  "%.5e + 2 * (%.5e - %.5e)",
1133  term1(0), term2a, term2b(0));
1134  risks.push_back(risk);
1135  }
1136  std::vector<double>::iterator it = min_element(risks.begin(), risks.end());
1137  int index = distance(risks.begin(), it);
1138  LOGL_DEBUG("TRACE3.ip.diffim.RegularizedKernelSolution.estimateRisk",
1139  "Minimum Risk = %.3e at lambda = %.3e", risks[index], lambdas[index]);
1140 
1141  return lambdas[index];
1142 
1143  }
#define LSST_EXCEPT_ADD(e, m)
Add the current location and a message to an existing exception before rethrowing it.
Definition: Exception.h:54
T min_element(T... args)
T push_back(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
virtualinherited

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
virtualinherited

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
virtualinherited

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  }

◆ getLambda()

template<typename InputT >
double lsst::ip::diffim::RegularizedKernelSolution< InputT >::getLambda ( )
inline

Definition at line 159 of file KernelSolution.h.

159 {return _lambda;}

◆ getM() [1/2]

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

Definition at line 64 of file KernelSolution.h.

64 {return _mMat;}

◆ getM() [2/2]

template<typename InputT >
Eigen::MatrixXd lsst::ip::diffim::RegularizedKernelSolution< InputT >::getM ( bool  includeHmat = true)

Definition at line 1146 of file KernelSolution.cc.

1146  {
1147  if (includeHmat == true) {
1148  return this->_mMat + _lambda * _hMat;
1149  }
1150  else {
1151  return this->_mMat;
1152  }
1153  }

◆ getSolutionPair()

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

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()

template<typename InputT >
std::shared_ptr< lsst::afw::image::Image< lsst::afw::math::Kernel::Pixel > > lsst::ip::diffim::StaticKernelSolution< InputT >::makeKernelImage
virtualinherited

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::RegularizedKernelSolution< InputT >::solve
virtual

Reimplemented from lsst::ip::diffim::StaticKernelSolution< InputT >.

Definition at line 1156 of file KernelSolution.cc.

1156  {
1157 
1158  LOGL_DEBUG("TRACE3.ip.diffim.RegularizedKernelSolution.solve",
1159  "cMat is %d x %d; vVec is %d; iVec is %d; hMat is %d x %d",
1160  this->_cMat.rows(), this->_cMat.cols(), this->_ivVec.size(),
1161  this->_iVec.size(), _hMat.rows(), _hMat.cols());
1162 
1163  if (DEBUG_MATRIX2) {
1164  std::cout << "ID: " << (this->_id) << std::endl;
1165  std::cout << "C:" << std::endl;
1166  std::cout << this->_cMat << std::endl;
1167  std::cout << "Sigma^{-1}:" << std::endl;
1168  std::cout << Eigen::MatrixXd(this->_ivVec.asDiagonal()) << std::endl;
1169  std::cout << "Y:" << std::endl;
1170  std::cout << this->_iVec << std::endl;
1171  std::cout << "H:" << std::endl;
1172  std::cout << _hMat << std::endl;
1173  }
1174 
1175 
1176  this->_mMat = this->_cMat.transpose() * this->_ivVec.asDiagonal() * this->_cMat;
1177  this->_bVec = this->_cMat.transpose() * this->_ivVec.asDiagonal() * this->_iVec;
1178 
1179 
1180  /* See N.R. 18.5
1181 
1182  Matrix equation to solve is Y = C a + N
1183  Y = vectorized version of I (I = image to not convolve)
1184  C_i = K_i (x) R (R = image to convolve)
1185  a = kernel coefficients
1186  N = noise
1187 
1188  If we reweight everything by the inverse square root of the noise
1189  covariance, we get a linear model with the identity matrix for
1190  the noise. The problem can then be solved using least squares,
1191  with the normal equations
1192 
1193  C^T Y = C^T C a
1194 
1195  or
1196 
1197  b = M a
1198 
1199  with
1200 
1201  b = C^T Y
1202  M = C^T C
1203  a = (C^T C)^{-1} C^T Y
1204 
1205 
1206  We can regularize the least square problem
1207 
1208  Y = C a + lambda a^T H a (+ N, which can be ignored)
1209 
1210  or the normal equations
1211 
1212  (C^T C + lambda H) a = C^T Y
1213 
1214 
1215  The solution to the regularization of the least squares problem is
1216 
1217  a = (C^T C + lambda H)^{-1} C^T Y
1218 
1219  The approximation to Y is
1220 
1221  C a = C (C^T C + lambda H)^{-1} C^T Y
1222 
1223  with smoothing matrix
1224 
1225  S = C (C^T C + lambda H)^{-1} C^T
1226 
1227  */
1228 
1229  std::string lambdaType = _ps->getAsString("lambdaType");
1230  if (lambdaType == "absolute") {
1231  _lambda = _ps->getAsDouble("lambdaValue");
1232  }
1233  else if (lambdaType == "relative") {
1234  _lambda = this->_mMat.trace() / this->_hMat.trace();
1235  _lambda *= _ps->getAsDouble("lambdaScaling");
1236  }
1237  else if (lambdaType == "minimizeBiasedRisk") {
1238  double tol = _ps->getAsDouble("maxConditionNumber");
1239  _lambda = estimateRisk(tol);
1240  }
1241  else if (lambdaType == "minimizeUnbiasedRisk") {
1243  }
1244  else {
1245  throw LSST_EXCEPT(pexExcept::Exception, "lambdaType in PropertySet not recognized");
1246  }
1247 
1248  LOGL_DEBUG("TRACE3.ip.diffim.RegularizedKernelSolution.solve",
1249  "Applying kernel regularization with lambda = %.2e", _lambda);
1250 
1251 
1252  try {
1253  KernelSolution::solve(this->_mMat + _lambda * _hMat, this->_bVec);
1254  } catch (pexExcept::Exception &e) {
1255  LSST_EXCEPT_ADD(e, "Unable to solve static kernel matrix");
1256  throw e;
1257  }
1258  /* Turn matrices into _kernel and _background */
1260  }
void _setKernel()
Set kernel after solution.
#define DEBUG_MATRIX2

◆ 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  }
#define DEBUG_MATRIX

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
protectedinherited

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
protectedinherited

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
protectedinherited

Vectorized I.

Definition at line 106 of file KernelSolution.h.

◆ _ivVec

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

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
protectedinherited

Derived single-object convolution kernel.

Definition at line 109 of file KernelSolution.h.

◆ _kSum

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

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: