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

#include <KernelSolution.h>

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

Public Types

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

 MaskedKernelSolution (lsst::afw::math::KernelList const &basisList, bool fitForBackground)
 
virtual ~MaskedKernelSolution ()
 
virtual void buildOrig (lsst::afw::image::Image< InputT > const &templateImage, lsst::afw::image::Image< InputT > const &scienceImage, lsst::afw::image::Image< lsst::afw::image::VariancePixel > const &varianceEstimate, lsst::afw::image::Mask< lsst::afw::image::MaskPixel > pixelMask)
 
virtual void buildWithMask (lsst::afw::image::Image< InputT > const &templateImage, lsst::afw::image::Image< InputT > const &scienceImage, lsst::afw::image::Image< lsst::afw::image::VariancePixel > const &varianceEstimate, lsst::afw::image::Mask< lsst::afw::image::MaskPixel > const &pixelMask)
 
virtual void buildSingleMaskOrig (lsst::afw::image::Image< InputT > const &templateImage, lsst::afw::image::Image< InputT > const &scienceImage, lsst::afw::image::Image< lsst::afw::image::VariancePixel > const &varianceEstimate, lsst::geom::Box2I maskBox)
 
void solve ()
 
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::MaskedKernelSolution< InputT >

Definition at line 119 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 121 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

◆ MaskedKernelSolution()

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

Definition at line 494 of file KernelSolution.cc.

497  :
498  StaticKernelSolution<InputT>(basisList, fitForBackground)
499  {};

◆ ~MaskedKernelSolution()

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

Definition at line 125 of file KernelSolution.h.

125 {};

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)

◆ buildOrig()

template<typename InputT >
void lsst::ip::diffim::MaskedKernelSolution< InputT >::buildOrig ( lsst::afw::image::Image< InputT > const &  templateImage,
lsst::afw::image::Image< InputT > const &  scienceImage,
lsst::afw::image::Image< lsst::afw::image::VariancePixel > const &  varianceEstimate,
lsst::afw::image::Mask< lsst::afw::image::MaskPixel pixelMask 
)
virtual

Definition at line 669 of file KernelSolution.cc.

674  {
675 
676  afwMath::Statistics varStats = afwMath::makeStatistics(varianceEstimate, afwMath::MIN);
677  if (varStats.getValue(afwMath::MIN) < 0.0) {
679  "Error: variance less than 0.0");
680  }
681  if (varStats.getValue(afwMath::MIN) == 0.0) {
683  "Error: variance equals 0.0, cannot inverse variance weight");
684  }
685 
686  lsst::afw::math::KernelList basisList =
687  std::dynamic_pointer_cast<afwMath::LinearCombinationKernel>(this->_kernel)->getKernelList();
688  std::vector<std::shared_ptr<afwMath::Kernel> >::const_iterator kiter = basisList.begin();
689 
690  /* Only BAD pixels marked in this mask */
692  afwImage::MaskPixel bitMask =
696  sMask &= bitMask;
697  /* TBD: Need to figure out a way to spread this; currently done in Python */
698 
699  unsigned int const nKernelParameters = basisList.size();
700  unsigned int const nBackgroundParameters = this->_fitForBackground ? 1 : 0;
701  unsigned int const nParameters = nKernelParameters + nBackgroundParameters;
702 
703  /* NOTE - we are accessing particular elements of Eigen arrays using
704  these coordinates, therefore they need to be in LOCAL coordinates.
705  This was written before ndarray unification.
706  */
707  /* Ignore known EDGE pixels for speed */
708  geom::Box2I shrunkLocalBBox = (*kiter)->shrinkBBox(templateImage.getBBox(afwImage::LOCAL));
709  LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
710  "Limits of good pixels after convolution: %d,%d -> %d,%d (local)",
711  shrunkLocalBBox.getMinX(), shrunkLocalBBox.getMinY(),
712  shrunkLocalBBox.getMaxX(), shrunkLocalBBox.getMaxY());
713 
714  /* Subimages are addressed in raw pixel coordinates */
715  unsigned int startCol = shrunkLocalBBox.getMinX();
716  unsigned int startRow = shrunkLocalBBox.getMinY();
717  unsigned int endCol = shrunkLocalBBox.getMaxX();
718  unsigned int endRow = shrunkLocalBBox.getMaxY();
719  /* Needed for Eigen block slicing operation */
720  endCol += 1;
721  endRow += 1;
722 
723  boost::timer t;
724  t.restart();
725 
726  /* Eigen representation of input images; only the pixels that are unconvolved in cimage below */
727  Eigen::MatrixXi eMask = maskToEigenMatrix(sMask).block(startRow,
728  startCol,
729  endRow-startRow,
730  endCol-startCol);
731 
732  Eigen::MatrixXd eigenTemplate = imageToEigenMatrix(templateImage).block(startRow,
733  startCol,
734  endRow-startRow,
735  endCol-startCol);
736  Eigen::MatrixXd eigenScience = imageToEigenMatrix(scienceImage).block(startRow,
737  startCol,
738  endRow-startRow,
739  endCol-startCol);
740  Eigen::MatrixXd eigeniVariance = imageToEigenMatrix(varianceEstimate).block(
741  startRow, startCol, endRow-startRow, endCol-startCol
742  ).array().inverse().matrix();
743 
744  /* Resize into 1-D for later usage */
745  eMask.resize(eMask.rows()*eMask.cols(), 1);
746  eigenTemplate.resize(eigenTemplate.rows()*eigenTemplate.cols(), 1);
747  eigenScience.resize(eigenScience.rows()*eigenScience.cols(), 1);
748  eigeniVariance.resize(eigeniVariance.rows()*eigeniVariance.cols(), 1);
749 
750  /* Do the masking, slowly for now... */
751  Eigen::MatrixXd maskedEigenTemplate(eigenTemplate.rows(), 1);
752  Eigen::MatrixXd maskedEigenScience(eigenScience.rows(), 1);
753  Eigen::MatrixXd maskedEigeniVariance(eigeniVariance.rows(), 1);
754  int nGood = 0;
755  for (int i = 0; i < eMask.rows(); i++) {
756  if (eMask(i, 0) == 0) {
757  maskedEigenTemplate(nGood, 0) = eigenTemplate(i, 0);
758  maskedEigenScience(nGood, 0) = eigenScience(i, 0);
759  maskedEigeniVariance(nGood, 0) = eigeniVariance(i, 0);
760  nGood += 1;
761  }
762  }
763  /* Can't resize() since all values are lost; use blocks */
764  eigenTemplate = maskedEigenTemplate.block(0, 0, nGood, 1);
765  eigenScience = maskedEigenScience.block(0, 0, nGood, 1);
766  eigeniVariance = maskedEigeniVariance.block(0, 0, nGood, 1);
767 
768 
769  /* Holds image convolved with basis function */
770  afwImage::Image<InputT> cimage(templateImage.getDimensions());
771 
772  /* Holds eigen representation of image convolved with all basis functions */
773  std::vector<Eigen::MatrixXd> convolvedEigenList(nKernelParameters);
774 
775  /* Iterators over convolved image list and basis list */
776  typename std::vector<Eigen::MatrixXd>::iterator eiter = convolvedEigenList.begin();
777  /* Create C_i in the formalism of Alard & Lupton */
778  afwMath::ConvolutionControl convolutionControl;
779  convolutionControl.setDoNormalize(false);
780  for (kiter = basisList.begin(); kiter != basisList.end(); ++kiter, ++eiter) {
781  afwMath::convolve(cimage, templateImage, **kiter, convolutionControl); /* cimage stores convolved image */
782 
783  Eigen::MatrixXd cMat = imageToEigenMatrix(cimage).block(startRow,
784  startCol,
785  endRow-startRow,
786  endCol-startCol);
787  cMat.resize(cMat.size(), 1);
788 
789  /* Do masking */
790  Eigen::MatrixXd maskedcMat(cMat.rows(), 1);
791  int nGood = 0;
792  for (int i = 0; i < eMask.rows(); i++) {
793  if (eMask(i, 0) == 0) {
794  maskedcMat(nGood, 0) = cMat(i, 0);
795  nGood += 1;
796  }
797  }
798  cMat = maskedcMat.block(0, 0, nGood, 1);
799  *eiter = cMat;
800  }
801 
802  double time = t.elapsed();
803  LOGL_DEBUG("TRACE3.ip.diffim.StaticKernelSolution.build",
804  "Total compute time to do basis convolutions : %.2f s", time);
805  t.restart();
806 
807  /*
808  Load matrix with all values from convolvedEigenList : all images
809  (eigeniVariance, convolvedEigenList) must be the same size
810  */
811  Eigen::MatrixXd cMat(eigenTemplate.col(0).size(), nParameters);
812  typename std::vector<Eigen::MatrixXd>::iterator eiterj = convolvedEigenList.begin();
813  typename std::vector<Eigen::MatrixXd>::iterator eiterE = convolvedEigenList.end();
814  for (unsigned int kidxj = 0; eiterj != eiterE; eiterj++, kidxj++) {
815  cMat.col(kidxj) = eiterj->col(0);
816  }
817  /* Treat the last "image" as all 1's to do the background calculation. */
818  if (this->_fitForBackground)
819  cMat.col(nParameters-1).fill(1.);
820 
821  this->_cMat = cMat;
822  this->_ivVec = eigeniVariance.col(0);
823  this->_iVec = eigenScience.col(0);
824 
825  /* Make these outside of solve() so I can check condition number */
826  this->_mMat = this->_cMat.transpose() * this->_ivVec.asDiagonal() * this->_cMat;
827  this->_bVec = this->_cMat.transpose() * this->_ivVec.asDiagonal() * this->_iVec;
828 
829  }
Represent a 2-dimensional array of bitmask pixels.
Definition: Mask.h:77
static MaskPixelT getPlaneBitMask(const std::vector< std::string > &names)
Return the bitmask corresponding to a vector of plane names OR'd together.
Definition: Mask.cc:372
Eigen::MatrixXi maskToEigenMatrix(lsst::afw::image::Mask< lsst::afw::image::MaskPixel > const &mask)

◆ buildSingleMaskOrig()

template<typename InputT >
void lsst::ip::diffim::MaskedKernelSolution< InputT >::buildSingleMaskOrig ( lsst::afw::image::Image< InputT > const &  templateImage,
lsst::afw::image::Image< InputT > const &  scienceImage,
lsst::afw::image::Image< lsst::afw::image::VariancePixel > const &  varianceEstimate,
lsst::geom::Box2I  maskBox 
)
virtual

Definition at line 834 of file KernelSolution.cc.

839  {
840 
841  afwMath::Statistics varStats = afwMath::makeStatistics(varianceEstimate, afwMath::MIN);
842  if (varStats.getValue(afwMath::MIN) < 0.0) {
844  "Error: variance less than 0.0");
845  }
846  if (varStats.getValue(afwMath::MIN) == 0.0) {
848  "Error: variance equals 0.0, cannot inverse variance weight");
849  }
850 
851  lsst::afw::math::KernelList basisList =
852  std::dynamic_pointer_cast<afwMath::LinearCombinationKernel>(this->_kernel)->getKernelList();
853 
854  unsigned int const nKernelParameters = basisList.size();
855  unsigned int const nBackgroundParameters = this->_fitForBackground ? 1 : 0;
856  unsigned int const nParameters = nKernelParameters + nBackgroundParameters;
857 
858  std::vector<std::shared_ptr<afwMath::Kernel> >::const_iterator kiter = basisList.begin();
859 
860  /*
861  NOTE : If we are using these views in Afw's Image space, we need to
862  make sure and compensate for the XY0 of the image:
863 
864  geom::Box2I fullBBox = templateImage.getBBox();
865  int maskStartCol = maskBox.getMinX();
866  int maskStartRow = maskBox.getMinY();
867  int maskEndCol = maskBox.getMaxX();
868  int maskEndRow = maskBox.getMaxY();
869 
870 
871  If we are going to be doing the slicing in Eigen matrices derived from
872  the images, ignore the XY0.
873 
874  geom::Box2I fullBBox = templateImage.getBBox(afwImage::LOCAL);
875 
876  int maskStartCol = maskBox.getMinX() - templateImage.getX0();
877  int maskStartRow = maskBox.getMinY() - templateImage.getY0();
878  int maskEndCol = maskBox.getMaxX() - templateImage.getX0();
879  int maskEndRow = maskBox.getMaxY() - templateImage.getY0();
880 
881  */
882 
883 
884  geom::Box2I shrunkBBox = (*kiter)->shrinkBBox(templateImage.getBBox());
885 
886  LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
887  "Limits of good pixels after convolution: %d,%d -> %d,%d",
888  shrunkBBox.getMinX(), shrunkBBox.getMinY(),
889  shrunkBBox.getMaxX(), shrunkBBox.getMaxY());
890 
891  unsigned int const startCol = shrunkBBox.getMinX();
892  unsigned int const startRow = shrunkBBox.getMinY();
893  unsigned int const endCol = shrunkBBox.getMaxX();
894  unsigned int const endRow = shrunkBBox.getMaxY();
895 
896  /* NOTE: no endCol/endRow += 1 for block slicing, since we are doing the
897  slicing using Afw, not Eigen
898 
899  Eigen arrays have index 0,0 in the upper right, while LSST images
900  have 0,0 in the lower left. The y-axis is flipped. When translating
901  Images to Eigen matrices in ipDiffim::imageToEigenMatrix this is
902  accounted for. However, we still need to be aware of this fact if
903  addressing subregions of an Eigen matrix. This is why the slicing is
904  done in Afw, its cleaner.
905 
906  Please see examples/maskedKernel.cc for elaboration on some of the
907  tests done to make sure this indexing gets done right.
908 
909  */
910 
911 
912  /* Inner limits; of mask box */
913  int maskStartCol = maskBox.getMinX();
914  int maskStartRow = maskBox.getMinY();
915  int maskEndCol = maskBox.getMaxX();
916  int maskEndRow = maskBox.getMaxY();
917 
918  /*
919 
920  |---------------------------|
921  | Kernel Boundary |
922  | |---------------------| |
923  | | Top | |
924  | |......_________......| |
925  | | | | | |
926  | | L | Box | R | |
927  | | | | | |
928  | |......---------......| |
929  | | Bottom | |
930  | |---------------------| |
931  | |
932  |---------------------------|
933 
934  4 regions we want to extract from the pixels: top bottom left right
935 
936  */
937  geom::Box2I tBox = geom::Box2I(geom::Point2I(startCol, maskEndRow + 1),
938  geom::Point2I(endCol, endRow));
939 
940  geom::Box2I bBox = geom::Box2I(geom::Point2I(startCol, startRow),
941  geom::Point2I(endCol, maskStartRow - 1));
942 
943  geom::Box2I lBox = geom::Box2I(geom::Point2I(startCol, maskStartRow),
944  geom::Point2I(maskStartCol - 1, maskEndRow));
945 
946  geom::Box2I rBox = geom::Box2I(geom::Point2I(maskEndCol + 1, maskStartRow),
947  geom::Point2I(endCol, maskEndRow));
948 
949  LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
950  "Upper good pixel region: %d,%d -> %d,%d",
951  tBox.getMinX(), tBox.getMinY(), tBox.getMaxX(), tBox.getMaxY());
952  LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
953  "Bottom good pixel region: %d,%d -> %d,%d",
954  bBox.getMinX(), bBox.getMinY(), bBox.getMaxX(), bBox.getMaxY());
955  LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
956  "Left good pixel region: %d,%d -> %d,%d",
957  lBox.getMinX(), lBox.getMinY(), lBox.getMaxX(), lBox.getMaxY());
958  LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
959  "Right good pixel region: %d,%d -> %d,%d",
960  rBox.getMinX(), rBox.getMinY(), rBox.getMaxX(), rBox.getMaxY());
961 
962  std::vector<geom::Box2I> boxArray;
963  boxArray.push_back(tBox);
964  boxArray.push_back(bBox);
965  boxArray.push_back(lBox);
966  boxArray.push_back(rBox);
967 
968  int totalSize = tBox.getWidth() * tBox.getHeight();
969  totalSize += bBox.getWidth() * bBox.getHeight();
970  totalSize += lBox.getWidth() * lBox.getHeight();
971  totalSize += rBox.getWidth() * rBox.getHeight();
972 
973  Eigen::MatrixXd eigenTemplate(totalSize, 1);
974  Eigen::MatrixXd eigenScience(totalSize, 1);
975  Eigen::MatrixXd eigeniVariance(totalSize, 1);
976  eigenTemplate.setZero();
977  eigenScience.setZero();
978  eigeniVariance.setZero();
979 
980  boost::timer t;
981  t.restart();
982 
983  int nTerms = 0;
984  typename std::vector<geom::Box2I>::iterator biter = boxArray.begin();
985  for (; biter != boxArray.end(); ++biter) {
986  int area = (*biter).getWidth() * (*biter).getHeight();
987 
988  afwImage::Image<InputT> siTemplate(templateImage, *biter);
989  afwImage::Image<InputT> siScience(scienceImage, *biter);
990  afwImage::Image<InputT> sVarEstimate(varianceEstimate, *biter);
991 
992  Eigen::MatrixXd eTemplate = imageToEigenMatrix(siTemplate);
993  Eigen::MatrixXd eScience = imageToEigenMatrix(siScience);
994  Eigen::MatrixXd eiVarEstimate = imageToEigenMatrix(sVarEstimate).array().inverse().matrix();
995 
996  eTemplate.resize(area, 1);
997  eScience.resize(area, 1);
998  eiVarEstimate.resize(area, 1);
999 
1000  eigenTemplate.block(nTerms, 0, area, 1) = eTemplate.block(0, 0, area, 1);
1001  eigenScience.block(nTerms, 0, area, 1) = eScience.block(0, 0, area, 1);
1002  eigeniVariance.block(nTerms, 0, area, 1) = eiVarEstimate.block(0, 0, area, 1);
1003 
1004  nTerms += area;
1005  }
1006 
1007  afwImage::Image<InputT> cimage(templateImage.getDimensions());
1008 
1009  std::vector<Eigen::MatrixXd> convolvedEigenList(nKernelParameters);
1010  typename std::vector<Eigen::MatrixXd>::iterator eiter = convolvedEigenList.begin();
1011  /* Create C_i in the formalism of Alard & Lupton */
1012  afwMath::ConvolutionControl convolutionControl;
1013  convolutionControl.setDoNormalize(false);
1014  for (kiter = basisList.begin(); kiter != basisList.end(); ++kiter, ++eiter) {
1015  afwMath::convolve(cimage, templateImage, **kiter, convolutionControl); /* cimage stores convolved image */
1016  Eigen::MatrixXd cMat(totalSize, 1);
1017  cMat.setZero();
1018 
1019  int nTerms = 0;
1020  typename std::vector<geom::Box2I>::iterator biter = boxArray.begin();
1021  for (; biter != boxArray.end(); ++biter) {
1022  int area = (*biter).getWidth() * (*biter).getHeight();
1023 
1024  afwImage::Image<InputT> csubimage(cimage, *biter);
1025  Eigen::MatrixXd esubimage = imageToEigenMatrix(csubimage);
1026  esubimage.resize(area, 1);
1027  cMat.block(nTerms, 0, area, 1) = esubimage.block(0, 0, area, 1);
1028 
1029  nTerms += area;
1030  }
1031 
1032  *eiter = cMat;
1033 
1034  }
1035 
1036  double time = t.elapsed();
1037  LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
1038  "Total compute time to do basis convolutions : %.2f s", time);
1039  t.restart();
1040 
1041  /*
1042  Load matrix with all values from convolvedEigenList : all images
1043  (eigeniVariance, convolvedEigenList) must be the same size
1044  */
1045  Eigen::MatrixXd cMat(eigenTemplate.col(0).size(), nParameters);
1046  typename std::vector<Eigen::MatrixXd>::iterator eiterj = convolvedEigenList.begin();
1047  typename std::vector<Eigen::MatrixXd>::iterator eiterE = convolvedEigenList.end();
1048  for (unsigned int kidxj = 0; eiterj != eiterE; eiterj++, kidxj++) {
1049  cMat.col(kidxj) = eiterj->col(0);
1050  }
1051  /* Treat the last "image" as all 1's to do the background calculation. */
1052  if (this->_fitForBackground)
1053  cMat.col(nParameters-1).fill(1.);
1054 
1055  this->_cMat = cMat;
1056  this->_ivVec = eigeniVariance.col(0);
1057  this->_iVec = eigenScience.col(0);
1058 
1059  /* Make these outside of solve() so I can check condition number */
1060  this->_mMat = this->_cMat.transpose() * this->_ivVec.asDiagonal() * this->_cMat;
1061  this->_bVec = this->_cMat.transpose() * this->_ivVec.asDiagonal() * this->_iVec;
1062  }
int getHeight() const noexcept
Definition: Box.h:188
int getWidth() const noexcept
Definition: Box.h:187
T push_back(T... args)

◆ buildWithMask()

template<typename InputT >
void lsst::ip::diffim::MaskedKernelSolution< InputT >::buildWithMask ( lsst::afw::image::Image< InputT > const &  templateImage,
lsst::afw::image::Image< InputT > const &  scienceImage,
lsst::afw::image::Image< lsst::afw::image::VariancePixel > const &  varianceEstimate,
lsst::afw::image::Mask< lsst::afw::image::MaskPixel > const &  pixelMask 
)
virtual

Definition at line 502 of file KernelSolution.cc.

507  {
508 
509  afwMath::Statistics varStats = afwMath::makeStatistics(varianceEstimate, afwMath::MIN);
510  if (varStats.getValue(afwMath::MIN) < 0.0) {
512  "Error: variance less than 0.0");
513  }
514  if (varStats.getValue(afwMath::MIN) == 0.0) {
516  "Error: variance equals 0.0, cannot inverse variance weight");
517  }
518 
519  /* Full footprint of all input images */
521  new afwDet::Footprint(std::make_shared<afwGeom::SpanSet>(templateImage.getBBox())));
522 
523  afwMath::KernelList basisList =
524  std::dynamic_pointer_cast<afwMath::LinearCombinationKernel>(this->_kernel)->getKernelList();
525  std::vector<std::shared_ptr<afwMath::Kernel> >::const_iterator kiter = basisList.begin();
526 
527  /* Only BAD pixels marked in this mask */
528  afwImage::MaskPixel bitMask =
533 
534  /* Create a Footprint that contains all the masked pixels set above */
536  afwDet::FootprintSet maskFpSet(pixelMask, threshold, true);
537 
538  /* And spread it by the kernel half width */
539  int growPix = (*kiter)->getCtr().getX();
540  afwDet::FootprintSet maskedFpSetGrown(maskFpSet, growPix, true);
541 
542 #if 0
543  for (typename afwDet::FootprintSet::FootprintList::iterator
544  ptr = maskedFpSetGrown.getFootprints()->begin(),
545  end = maskedFpSetGrown.getFootprints()->end();
546  ptr != end;
547  ++ptr) {
548 
549  afwDet::setMaskFromFootprint(finalMask,
550  (**ptr),
552  }
553 #endif
554 
555  afwImage::Mask<afwImage::MaskPixel> finalMask(pixelMask.getDimensions());
556  for (auto const & foot : *(maskedFpSetGrown.getFootprints())) {
557  foot->getSpans()->setMask(finalMask, afwImage::Mask<afwImage::MaskPixel>::getPlaneBitMask("BAD"));
558  }
559  pixelMask.writeFits("pixelmask.fits");
560  finalMask.writeFits("finalmask.fits");
561 
562 
563  ndarray::Array<int, 1, 1> maskArray =
564  ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
565  fullFp->getSpans()->flatten(maskArray, finalMask.getArray(), templateImage.getXY0());
566  auto maskEigen = ndarray::asEigenMatrix(maskArray);
567 
568  ndarray::Array<InputT, 1, 1> arrayTemplate =
569  ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
570  fullFp->getSpans()->flatten(arrayTemplate, templateImage.getArray(), templateImage.getXY0());
571  auto eigenTemplate0 = ndarray::asEigenMatrix(arrayTemplate);
572 
573  ndarray::Array<InputT, 1, 1> arrayScience =
574  ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
575  fullFp->getSpans()->flatten(arrayScience, scienceImage.getArray(), scienceImage.getXY0());
576  auto eigenScience0 = ndarray::asEigenMatrix(arrayScience);
577 
578  ndarray::Array<afwImage::VariancePixel, 1, 1> arrayVariance =
579  ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
580  fullFp->getSpans()->flatten(arrayVariance, varianceEstimate.getArray(), varianceEstimate.getXY0());
581  auto eigenVariance0 = ndarray::asEigenMatrix(arrayVariance);
582 
583  int nGood = 0;
584  for (int i = 0; i < maskEigen.size(); i++) {
585  if (maskEigen(i) == 0.0)
586  nGood += 1;
587  }
588 
589  Eigen::VectorXd eigenTemplate(nGood);
590  Eigen::VectorXd eigenScience(nGood);
591  Eigen::VectorXd eigenVariance(nGood);
592  int nUsed = 0;
593  for (int i = 0; i < maskEigen.size(); i++) {
594  if (maskEigen(i) == 0.0) {
595  eigenTemplate(nUsed) = eigenTemplate0(i);
596  eigenScience(nUsed) = eigenScience0(i);
597  eigenVariance(nUsed) = eigenVariance0(i);
598  nUsed += 1;
599  }
600  }
601 
602 
603  boost::timer t;
604  t.restart();
605 
606  unsigned int const nKernelParameters = basisList.size();
607  unsigned int const nBackgroundParameters = this->_fitForBackground ? 1 : 0;
608  unsigned int const nParameters = nKernelParameters + nBackgroundParameters;
609 
610  /* Holds image convolved with basis function */
611  afwImage::Image<InputT> cimage(templateImage.getDimensions());
612 
613  /* Holds eigen representation of image convolved with all basis functions */
614  std::vector<Eigen::VectorXd> convolvedEigenList(nKernelParameters);
615 
616  /* Iterators over convolved image list and basis list */
617  typename std::vector<Eigen::VectorXd>::iterator eiter = convolvedEigenList.begin();
618 
619  /* Create C_i in the formalism of Alard & Lupton */
620  afwMath::ConvolutionControl convolutionControl;
621  convolutionControl.setDoNormalize(false);
622  for (kiter = basisList.begin(); kiter != basisList.end(); ++kiter, ++eiter) {
623  afwMath::convolve(cimage, templateImage, **kiter, convolutionControl); /* cimage stores convolved image */
624 
625  ndarray::Array<InputT, 1, 1> arrayC =
626  ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
627  fullFp->getSpans()->flatten(arrayC, cimage.getArray(), cimage.getXY0());
628  auto eigenC0 = ndarray::asEigenMatrix(arrayC);
629 
630  Eigen::VectorXd eigenC(nGood);
631  int nUsed = 0;
632  for (int i = 0; i < maskEigen.size(); i++) {
633  if (maskEigen(i) == 0.0) {
634  eigenC(nUsed) = eigenC0(i);
635  nUsed += 1;
636  }
637  }
638 
639  *eiter = eigenC;
640  }
641  double time = t.elapsed();
642  LOGL_DEBUG("TRACE3.ip.diffim.StaticKernelSolution.buildWithMask",
643  "Total compute time to do basis convolutions : %.2f s", time);
644  t.restart();
645 
646  /* Load matrix with all convolved images */
647  Eigen::MatrixXd cMat(eigenTemplate.size(), nParameters);
648  typename std::vector<Eigen::VectorXd>::iterator eiterj = convolvedEigenList.begin();
649  typename std::vector<Eigen::VectorXd>::iterator eiterE = convolvedEigenList.end();
650  for (unsigned int kidxj = 0; eiterj != eiterE; eiterj++, kidxj++) {
651  cMat.block(0, kidxj, eigenTemplate.size(), 1) =
652  Eigen::MatrixXd(*eiterj).block(0, 0, eigenTemplate.size(), 1);
653  }
654  /* Treat the last "image" as all 1's to do the background calculation. */
655  if (this->_fitForBackground)
656  cMat.col(nParameters-1).fill(1.);
657 
658  this->_cMat = cMat;
659  this->_ivVec = eigenVariance.array().inverse().matrix();
660  this->_iVec = eigenScience;
661 
662  /* Make these outside of solve() so I can check condition number */
663  this->_mMat = this->_cMat.transpose() * this->_ivVec.asDiagonal() * (this->_cMat);
664  this->_bVec = this->_cMat.transpose() * this->_ivVec.asDiagonal() * (this->_iVec);
665  }
int end
uint64_t * ptr
Definition: RangeSet.cc:88
Class to describe the properties of a detected object from an image.
Definition: Footprint.h:63
A set of Footprints, associated with a MaskedImage.
Definition: FootprintSet.h:53
A Threshold is used to pass a threshold value to detection algorithms.
Definition: Threshold.h:43
@ BITMASK
Use (pixels & (given mask))
Definition: Threshold.h:48
lsst::geom::Point2I getXY0() const
Return the image's origin.
Definition: ImageBase.h:323
void writeFits(std::string const &fileName, std::shared_ptr< lsst::daf::base::PropertySet const > metadata=std::shared_ptr< lsst::daf::base::PropertySet >(), std::string const &mode="w") const
Write a mask to a regular FITS file.

◆ 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  }

◆ 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
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::StaticKernelSolution< InputT >::solve
virtualinherited

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
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: