LSSTApplications  17.0+124,17.0+14,17.0+73,18.0.0+37,18.0.0+80,18.0.0-4-g68ffd23+4,18.1.0-1-g0001055+12,18.1.0-1-g03d53ef+5,18.1.0-1-g1349e88+55,18.1.0-1-g2505f39+44,18.1.0-1-g5315e5e+4,18.1.0-1-g5e4b7ea+14,18.1.0-1-g7e8fceb+4,18.1.0-1-g85f8cd4+48,18.1.0-1-g8ff0b9f+4,18.1.0-1-ga2c679d+1,18.1.0-1-gd55f500+35,18.1.0-10-gb58edde+2,18.1.0-11-g0997b02+4,18.1.0-13-gfe4edf0b+12,18.1.0-14-g259bd21+21,18.1.0-19-gdb69f3f+2,18.1.0-2-g5f9922c+24,18.1.0-2-gd3b74e5+11,18.1.0-2-gfbf3545+32,18.1.0-26-g728bddb4+5,18.1.0-27-g6ff7ca9+2,18.1.0-3-g52aa583+25,18.1.0-3-g8ea57af+9,18.1.0-3-gb69f684+42,18.1.0-3-gfcaddf3+6,18.1.0-32-gd8786685a,18.1.0-4-gf3f9b77+6,18.1.0-5-g1dd662b+2,18.1.0-5-g6dbcb01+41,18.1.0-6-gae77429+3,18.1.0-7-g9d75d83+9,18.1.0-7-gae09a6d+30,18.1.0-9-gc381ef5+4,w.2019.45
LSSTDataManagementBasePackage
KernelSolution.cc
Go to the documentation of this file.
1 // -*- lsst-c++ -*-
11 #include <iterator>
12 #include <cmath>
13 #include <algorithm>
14 #include <limits>
15 
16 #include <memory>
17 #include "boost/timer.hpp"
18 
19 #include "Eigen/Core"
20 #include "Eigen/Cholesky"
21 #include "Eigen/QR"
22 #include "Eigen/LU"
23 #include "Eigen/Eigenvalues"
24 #include "Eigen/SVD"
25 
26 #include "lsst/afw/math.h"
27 #include "lsst/afw/geom.h"
28 #include "lsst/afw/image.h"
29 #include "lsst/afw/detection.h"
30 #include "lsst/geom.h"
31 #include "lsst/log/Log.h"
33 
36 
37 #include "ndarray.h"
38 #include "ndarray/eigen.h"
39 
40 #define DEBUG_MATRIX 0
41 #define DEBUG_MATRIX2 0
42 
43 namespace geom = lsst::geom;
44 namespace afwDet = lsst::afw::detection;
45 namespace afwMath = lsst::afw::math;
46 namespace afwGeom = lsst::afw::geom;
47 namespace afwImage = lsst::afw::image;
49 
50 namespace lsst {
51 namespace ip {
52 namespace diffim {
53 
54  /* Unique identifier for solution */
56 
58  Eigen::MatrixXd mMat,
59  Eigen::VectorXd bVec,
60  bool fitForBackground
61  ) :
62  _id(++_SolutionId),
63  _mMat(mMat),
64  _bVec(bVec),
65  _aVec(),
66  _solvedBy(NONE),
67  _fitForBackground(fitForBackground)
68  {};
69 
71  bool fitForBackground
72  ) :
73  _id(++_SolutionId),
74  _mMat(),
75  _bVec(),
76  _aVec(),
77  _solvedBy(NONE),
78  _fitForBackground(fitForBackground)
79  {};
80 
82  _id(++_SolutionId),
83  _mMat(),
84  _bVec(),
85  _aVec(),
86  _solvedBy(NONE),
87  _fitForBackground(true)
88  {};
89 
91  solve(_mMat, _bVec);
92  }
93 
95  return getConditionNumber(_mMat, conditionType);
96  }
97 
98  double KernelSolution::getConditionNumber(Eigen::MatrixXd const& mMat,
99  ConditionNumberType conditionType) {
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  }
130 
131  void KernelSolution::solve(Eigen::MatrixXd const& mMat,
132  Eigen::VectorXd const& bVec) {
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  }
191 
192  /*******************************************************************************************************/
193 
194  template <typename InputT>
196  lsst::afw::math::KernelList const& basisList,
197  bool fitForBackground
198  )
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  };
213 
214  template <typename InputT>
217  throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return solution");
218  }
219  return _kernel;
220  }
221 
222  template <typename InputT>
225  throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return image");
226  }
229  );
230  (void)_kernel->computeImage(*image, false);
231  return image;
232  }
233 
234  template <typename InputT>
237  throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return background");
238  }
239  return _background;
240  }
241 
242  template <typename InputT>
245  throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return ksum");
246  }
247  return _kSum;
248  }
249 
250  template <typename InputT>
254  throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return solution");
255  }
256 
258  }
259 
260  template <typename InputT>
262  lsst::afw::image::Image<InputT> const &templateImage,
263  lsst::afw::image::Image<InputT> const &scienceImage,
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 =
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  for (kiter = basisList.begin(); kiter != basisList.end(); ++kiter, ++eiter) {
356  afwMath::convolve(cimage, templateImage, **kiter, false); /* cimage stores convolved image */
357 
358  Eigen::MatrixXd cMat = imageToEigenMatrix(cimage).block(startRow,
359  startCol,
360  endRow-startRow,
361  endCol-startCol);
362  cMat.resize(cMat.size(), 1);
363  *eiter = cMat;
364 
365  }
366 
367  double time = t.elapsed();
368  LOGL_DEBUG("TRACE3.ip.diffim.StaticKernelSolution.build",
369  "Total compute time to do basis convolutions : %.2f s", time);
370  t.restart();
371 
372  /*
373  Load matrix with all values from convolvedEigenList : all images
374  (eigeniVariance, convolvedEigenList) must be the same size
375  */
376  Eigen::MatrixXd cMat(eigenTemplate.col(0).size(), nParameters);
377  typename std::vector<Eigen::MatrixXd>::iterator eiterj = convolvedEigenList.begin();
378  typename std::vector<Eigen::MatrixXd>::iterator eiterE = convolvedEigenList.end();
379  for (unsigned int kidxj = 0; eiterj != eiterE; eiterj++, kidxj++) {
380  cMat.col(kidxj) = eiterj->col(0);
381  }
382  /* Treat the last "image" as all 1's to do the background calculation. */
383  if (_fitForBackground)
384  cMat.col(nParameters-1).fill(1.);
385 
386  _cMat = cMat;
387  _ivVec = eigeniVariance.col(0);
388  _iVec = eigenScience.col(0);
389 
390  /* Make these outside of solve() so I can check condition number */
391  _mMat = _cMat.transpose() * (_ivVec.asDiagonal() * _cMat);
392  _bVec = _cMat.transpose() * (_ivVec.asDiagonal() * _iVec);
393  }
394 
395  template <typename InputT>
397  LOGL_DEBUG("TRACE3.ip.diffim.StaticKernelSolution.solve",
398  "mMat is %d x %d; bVec is %d; cMat is %d x %d; vVec is %d; iVec is %d",
399  _mMat.rows(), _mMat.cols(), _bVec.size(),
400  _cMat.rows(), _cMat.cols(), _ivVec.size(), _iVec.size());
401 
402  if (DEBUG_MATRIX) {
403  std::cout << "C" << std::endl;
404  std::cout << _cMat << std::endl;
405  std::cout << "iV" << std::endl;
406  std::cout << _ivVec << std::endl;
407  std::cout << "I" << std::endl;
408  std::cout << _iVec << std::endl;
409  }
410 
411  try {
413  } catch (pexExcept::Exception &e) {
414  LSST_EXCEPT_ADD(e, "Unable to solve static kernel matrix");
415  throw e;
416  }
417  /* Turn matrices into _kernel and _background */
418  _setKernel();
419  }
420 
421  template <typename InputT>
424  throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot make solution");
425  }
426 
427  unsigned int const nParameters = _aVec.size();
428  unsigned int const nBackgroundParameters = _fitForBackground ? 1 : 0;
429  unsigned int const nKernelParameters =
431  if (nParameters != (nKernelParameters + nBackgroundParameters))
432  throw LSST_EXCEPT(pexExcept::Exception, "Mismatched sizes in kernel solution");
433 
434  /* Fill in the kernel results */
435  std::vector<double> kValues(nKernelParameters);
436  for (unsigned int idx = 0; idx < nKernelParameters; idx++) {
437  if (std::isnan(_aVec(idx))) {
439  str(boost::format("Unable to determine kernel solution %d (nan)") % idx));
440  }
441  kValues[idx] = _aVec(idx);
442  }
443  _kernel->setKernelParameters(kValues);
444 
446  new ImageT(_kernel->getDimensions())
447  );
448  _kSum = _kernel->computeImage(*image, false);
449 
450  if (_fitForBackground) {
451  if (std::isnan(_aVec(nParameters-1))) {
453  str(boost::format("Unable to determine background solution %d (nan)") %
454  (nParameters-1)));
455  }
456  _background = _aVec(nParameters-1);
457  }
458  }
459 
460 
461  template <typename InputT>
463  throw LSST_EXCEPT(pexExcept::Exception, "Uncertainty calculation not supported");
464 
465  /* Estimate of parameter uncertainties comes from the inverse of the
466  * covariance matrix (noise spectrum).
467  * N.R. 15.4.8 to 15.4.15
468  *
469  * Since this is a linear problem no need to use Fisher matrix
470  * N.R. 15.5.8
471  *
472  * Although I might be able to take advantage of the solution above.
473  * Since this now works and is not the rate limiting step, keep as-is for DC3a.
474  *
475  * Use Cholesky decomposition again.
476  * Cholkesy:
477  * Cov = L L^t
478  * Cov^(-1) = (L L^t)^(-1)
479  * = (L^T)^-1 L^(-1)
480  *
481  * Code would be:
482  *
483  * Eigen::MatrixXd cov = _mMat.transpose() * _mMat;
484  * Eigen::LLT<Eigen::MatrixXd> llt = cov.llt();
485  * Eigen::MatrixXd error2 = llt.matrixL().transpose().inverse()*llt.matrixL().inverse();
486  */
487  }
488 
489  /*******************************************************************************************************/
490 
491  template <typename InputT>
493  lsst::afw::math::KernelList const& basisList,
494  bool fitForBackground
495  ) :
496  StaticKernelSolution<InputT>(basisList, fitForBackground)
497  {};
498 
499  template <typename InputT>
501  lsst::afw::image::Image<InputT> const &templateImage,
502  lsst::afw::image::Image<InputT> const &scienceImage,
505  ) {
506 
507  afwMath::Statistics varStats = afwMath::makeStatistics(varianceEstimate, afwMath::MIN);
508  if (varStats.getValue(afwMath::MIN) < 0.0) {
510  "Error: variance less than 0.0");
511  }
512  if (varStats.getValue(afwMath::MIN) == 0.0) {
514  "Error: variance equals 0.0, cannot inverse variance weight");
515  }
516 
517  /* Full footprint of all input images */
519  new afwDet::Footprint(std::make_shared<afwGeom::SpanSet>(templateImage.getBBox())));
520 
521  afwMath::KernelList basisList =
523  std::vector<std::shared_ptr<afwMath::Kernel> >::const_iterator kiter = basisList.begin();
524 
525  /* Only BAD pixels marked in this mask */
526  afwImage::MaskPixel bitMask =
531 
532  /* Create a Footprint that contains all the masked pixels set above */
534  afwDet::FootprintSet maskFpSet(pixelMask, threshold, true);
535 
536  /* And spread it by the kernel half width */
537  int growPix = (*kiter)->getCtr().getX();
538  afwDet::FootprintSet maskedFpSetGrown(maskFpSet, growPix, true);
539 
540 #if 0
541  for (typename afwDet::FootprintSet::FootprintList::iterator
542  ptr = maskedFpSetGrown.getFootprints()->begin(),
543  end = maskedFpSetGrown.getFootprints()->end();
544  ptr != end;
545  ++ptr) {
546 
547  afwDet::setMaskFromFootprint(finalMask,
548  (**ptr),
550  }
551 #endif
552 
553  afwImage::Mask<afwImage::MaskPixel> finalMask(pixelMask.getDimensions());
554  for (auto const & foot : *(maskedFpSetGrown.getFootprints())) {
555  foot->getSpans()->setMask(finalMask, afwImage::Mask<afwImage::MaskPixel>::getPlaneBitMask("BAD"));
556  }
557  pixelMask.writeFits("pixelmask.fits");
558  finalMask.writeFits("finalmask.fits");
559 
560 
561  ndarray::Array<int, 1, 1> maskArray =
562  ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
563  fullFp->getSpans()->flatten(maskArray, finalMask.getArray(), templateImage.getXY0());
564  auto maskEigen = ndarray::asEigenMatrix(maskArray);
565 
566  ndarray::Array<InputT, 1, 1> arrayTemplate =
567  ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
568  fullFp->getSpans()->flatten(arrayTemplate, templateImage.getArray(), templateImage.getXY0());
569  auto eigenTemplate0 = ndarray::asEigenMatrix(arrayTemplate);
570 
571  ndarray::Array<InputT, 1, 1> arrayScience =
572  ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
573  fullFp->getSpans()->flatten(arrayScience, scienceImage.getArray(), scienceImage.getXY0());
574  auto eigenScience0 = ndarray::asEigenMatrix(arrayScience);
575 
576  ndarray::Array<afwImage::VariancePixel, 1, 1> arrayVariance =
577  ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
578  fullFp->getSpans()->flatten(arrayVariance, varianceEstimate.getArray(), varianceEstimate.getXY0());
579  auto eigenVariance0 = ndarray::asEigenMatrix(arrayVariance);
580 
581  int nGood = 0;
582  for (int i = 0; i < maskEigen.size(); i++) {
583  if (maskEigen(i) == 0.0)
584  nGood += 1;
585  }
586 
587  Eigen::VectorXd eigenTemplate(nGood);
588  Eigen::VectorXd eigenScience(nGood);
589  Eigen::VectorXd eigenVariance(nGood);
590  int nUsed = 0;
591  for (int i = 0; i < maskEigen.size(); i++) {
592  if (maskEigen(i) == 0.0) {
593  eigenTemplate(nUsed) = eigenTemplate0(i);
594  eigenScience(nUsed) = eigenScience0(i);
595  eigenVariance(nUsed) = eigenVariance0(i);
596  nUsed += 1;
597  }
598  }
599 
600 
601  boost::timer t;
602  t.restart();
603 
604  unsigned int const nKernelParameters = basisList.size();
605  unsigned int const nBackgroundParameters = this->_fitForBackground ? 1 : 0;
606  unsigned int const nParameters = nKernelParameters + nBackgroundParameters;
607 
608  /* Holds image convolved with basis function */
609  afwImage::Image<InputT> cimage(templateImage.getDimensions());
610 
611  /* Holds eigen representation of image convolved with all basis functions */
612  std::vector<Eigen::VectorXd> convolvedEigenList(nKernelParameters);
613 
614  /* Iterators over convolved image list and basis list */
615  typename std::vector<Eigen::VectorXd>::iterator eiter = convolvedEigenList.begin();
616 
617  /* Create C_i in the formalism of Alard & Lupton */
618  for (kiter = basisList.begin(); kiter != basisList.end(); ++kiter, ++eiter) {
619  afwMath::convolve(cimage, templateImage, **kiter, false); /* cimage stores convolved image */
620 
621  ndarray::Array<InputT, 1, 1> arrayC =
622  ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
623  fullFp->getSpans()->flatten(arrayC, cimage.getArray(), cimage.getXY0());
624  auto eigenC0 = ndarray::asEigenMatrix(arrayC);
625 
626  Eigen::VectorXd eigenC(nGood);
627  int nUsed = 0;
628  for (int i = 0; i < maskEigen.size(); i++) {
629  if (maskEigen(i) == 0.0) {
630  eigenC(nUsed) = eigenC0(i);
631  nUsed += 1;
632  }
633  }
634 
635  *eiter = eigenC;
636  }
637  double time = t.elapsed();
638  LOGL_DEBUG("TRACE3.ip.diffim.StaticKernelSolution.buildWithMask",
639  "Total compute time to do basis convolutions : %.2f s", time);
640  t.restart();
641 
642  /* Load matrix with all convolved images */
643  Eigen::MatrixXd cMat(eigenTemplate.size(), nParameters);
644  typename std::vector<Eigen::VectorXd>::iterator eiterj = convolvedEigenList.begin();
645  typename std::vector<Eigen::VectorXd>::iterator eiterE = convolvedEigenList.end();
646  for (unsigned int kidxj = 0; eiterj != eiterE; eiterj++, kidxj++) {
647  cMat.block(0, kidxj, eigenTemplate.size(), 1) =
648  Eigen::MatrixXd(*eiterj).block(0, 0, eigenTemplate.size(), 1);
649  }
650  /* Treat the last "image" as all 1's to do the background calculation. */
651  if (this->_fitForBackground)
652  cMat.col(nParameters-1).fill(1.);
653 
654  this->_cMat = cMat;
655  this->_ivVec = eigenVariance.array().inverse().matrix();
656  this->_iVec = eigenScience;
657 
658  /* Make these outside of solve() so I can check condition number */
659  this->_mMat = this->_cMat.transpose() * this->_ivVec.asDiagonal() * (this->_cMat);
660  this->_bVec = this->_cMat.transpose() * this->_ivVec.asDiagonal() * (this->_iVec);
661  }
662 
663 
664  template <typename InputT>
666  lsst::afw::image::Image<InputT> const &templateImage,
667  lsst::afw::image::Image<InputT> const &scienceImage,
670  ) {
671 
672  afwMath::Statistics varStats = afwMath::makeStatistics(varianceEstimate, afwMath::MIN);
673  if (varStats.getValue(afwMath::MIN) < 0.0) {
675  "Error: variance less than 0.0");
676  }
677  if (varStats.getValue(afwMath::MIN) == 0.0) {
679  "Error: variance equals 0.0, cannot inverse variance weight");
680  }
681 
682  lsst::afw::math::KernelList basisList =
684  std::vector<std::shared_ptr<afwMath::Kernel> >::const_iterator kiter = basisList.begin();
685 
686  /* Only BAD pixels marked in this mask */
688  afwImage::MaskPixel bitMask =
692  sMask &= bitMask;
693  /* TBD: Need to figure out a way to spread this; currently done in Python */
694 
695  unsigned int const nKernelParameters = basisList.size();
696  unsigned int const nBackgroundParameters = this->_fitForBackground ? 1 : 0;
697  unsigned int const nParameters = nKernelParameters + nBackgroundParameters;
698 
699  /* NOTE - we are accessing particular elements of Eigen arrays using
700  these coordinates, therefore they need to be in LOCAL coordinates.
701  This was written before ndarray unification.
702  */
703  /* Ignore known EDGE pixels for speed */
704  geom::Box2I shrunkLocalBBox = (*kiter)->shrinkBBox(templateImage.getBBox(afwImage::LOCAL));
705  LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
706  "Limits of good pixels after convolution: %d,%d -> %d,%d (local)",
707  shrunkLocalBBox.getMinX(), shrunkLocalBBox.getMinY(),
708  shrunkLocalBBox.getMaxX(), shrunkLocalBBox.getMaxY());
709 
710  /* Subimages are addressed in raw pixel coordinates */
711  unsigned int startCol = shrunkLocalBBox.getMinX();
712  unsigned int startRow = shrunkLocalBBox.getMinY();
713  unsigned int endCol = shrunkLocalBBox.getMaxX();
714  unsigned int endRow = shrunkLocalBBox.getMaxY();
715  /* Needed for Eigen block slicing operation */
716  endCol += 1;
717  endRow += 1;
718 
719  boost::timer t;
720  t.restart();
721 
722  /* Eigen representation of input images; only the pixels that are unconvolved in cimage below */
723  Eigen::MatrixXi eMask = maskToEigenMatrix(sMask).block(startRow,
724  startCol,
725  endRow-startRow,
726  endCol-startCol);
727 
728  Eigen::MatrixXd eigenTemplate = imageToEigenMatrix(templateImage).block(startRow,
729  startCol,
730  endRow-startRow,
731  endCol-startCol);
732  Eigen::MatrixXd eigenScience = imageToEigenMatrix(scienceImage).block(startRow,
733  startCol,
734  endRow-startRow,
735  endCol-startCol);
736  Eigen::MatrixXd eigeniVariance = imageToEigenMatrix(varianceEstimate).block(
737  startRow, startCol, endRow-startRow, endCol-startCol
738  ).array().inverse().matrix();
739 
740  /* Resize into 1-D for later usage */
741  eMask.resize(eMask.rows()*eMask.cols(), 1);
742  eigenTemplate.resize(eigenTemplate.rows()*eigenTemplate.cols(), 1);
743  eigenScience.resize(eigenScience.rows()*eigenScience.cols(), 1);
744  eigeniVariance.resize(eigeniVariance.rows()*eigeniVariance.cols(), 1);
745 
746  /* Do the masking, slowly for now... */
747  Eigen::MatrixXd maskedEigenTemplate(eigenTemplate.rows(), 1);
748  Eigen::MatrixXd maskedEigenScience(eigenScience.rows(), 1);
749  Eigen::MatrixXd maskedEigeniVariance(eigeniVariance.rows(), 1);
750  int nGood = 0;
751  for (int i = 0; i < eMask.rows(); i++) {
752  if (eMask(i, 0) == 0) {
753  maskedEigenTemplate(nGood, 0) = eigenTemplate(i, 0);
754  maskedEigenScience(nGood, 0) = eigenScience(i, 0);
755  maskedEigeniVariance(nGood, 0) = eigeniVariance(i, 0);
756  nGood += 1;
757  }
758  }
759  /* Can't resize() since all values are lost; use blocks */
760  eigenTemplate = maskedEigenTemplate.block(0, 0, nGood, 1);
761  eigenScience = maskedEigenScience.block(0, 0, nGood, 1);
762  eigeniVariance = maskedEigeniVariance.block(0, 0, nGood, 1);
763 
764 
765  /* Holds image convolved with basis function */
766  afwImage::Image<InputT> cimage(templateImage.getDimensions());
767 
768  /* Holds eigen representation of image convolved with all basis functions */
769  std::vector<Eigen::MatrixXd> convolvedEigenList(nKernelParameters);
770 
771  /* Iterators over convolved image list and basis list */
772  typename std::vector<Eigen::MatrixXd>::iterator eiter = convolvedEigenList.begin();
773  /* Create C_i in the formalism of Alard & Lupton */
774  for (kiter = basisList.begin(); kiter != basisList.end(); ++kiter, ++eiter) {
775  afwMath::convolve(cimage, templateImage, **kiter, false); /* cimage stores convolved image */
776 
777  Eigen::MatrixXd cMat = imageToEigenMatrix(cimage).block(startRow,
778  startCol,
779  endRow-startRow,
780  endCol-startCol);
781  cMat.resize(cMat.size(), 1);
782 
783  /* Do masking */
784  Eigen::MatrixXd maskedcMat(cMat.rows(), 1);
785  int nGood = 0;
786  for (int i = 0; i < eMask.rows(); i++) {
787  if (eMask(i, 0) == 0) {
788  maskedcMat(nGood, 0) = cMat(i, 0);
789  nGood += 1;
790  }
791  }
792  cMat = maskedcMat.block(0, 0, nGood, 1);
793  *eiter = cMat;
794  }
795 
796  double time = t.elapsed();
797  LOGL_DEBUG("TRACE3.ip.diffim.StaticKernelSolution.build",
798  "Total compute time to do basis convolutions : %.2f s", time);
799  t.restart();
800 
801  /*
802  Load matrix with all values from convolvedEigenList : all images
803  (eigeniVariance, convolvedEigenList) must be the same size
804  */
805  Eigen::MatrixXd cMat(eigenTemplate.col(0).size(), nParameters);
806  typename std::vector<Eigen::MatrixXd>::iterator eiterj = convolvedEigenList.begin();
807  typename std::vector<Eigen::MatrixXd>::iterator eiterE = convolvedEigenList.end();
808  for (unsigned int kidxj = 0; eiterj != eiterE; eiterj++, kidxj++) {
809  cMat.col(kidxj) = eiterj->col(0);
810  }
811  /* Treat the last "image" as all 1's to do the background calculation. */
812  if (this->_fitForBackground)
813  cMat.col(nParameters-1).fill(1.);
814 
815  this->_cMat = cMat;
816  this->_ivVec = eigeniVariance.col(0);
817  this->_iVec = eigenScience.col(0);
818 
819  /* Make these outside of solve() so I can check condition number */
820  this->_mMat = this->_cMat.transpose() * this->_ivVec.asDiagonal() * this->_cMat;
821  this->_bVec = this->_cMat.transpose() * this->_ivVec.asDiagonal() * this->_iVec;
822 
823  }
824 
825  /* NOTE - this was written before the ndarray unification. I am rewriting
826  buildSingleMask to use the ndarray stuff */
827  template <typename InputT>
829  lsst::afw::image::Image<InputT> const &templateImage,
830  lsst::afw::image::Image<InputT> const &scienceImage,
832  lsst::geom::Box2I maskBox
833  ) {
834 
835  afwMath::Statistics varStats = afwMath::makeStatistics(varianceEstimate, afwMath::MIN);
836  if (varStats.getValue(afwMath::MIN) < 0.0) {
838  "Error: variance less than 0.0");
839  }
840  if (varStats.getValue(afwMath::MIN) == 0.0) {
842  "Error: variance equals 0.0, cannot inverse variance weight");
843  }
844 
845  lsst::afw::math::KernelList basisList =
847 
848  unsigned int const nKernelParameters = basisList.size();
849  unsigned int const nBackgroundParameters = this->_fitForBackground ? 1 : 0;
850  unsigned int const nParameters = nKernelParameters + nBackgroundParameters;
851 
852  std::vector<std::shared_ptr<afwMath::Kernel> >::const_iterator kiter = basisList.begin();
853 
854  /*
855  NOTE : If we are using these views in Afw's Image space, we need to
856  make sure and compensate for the XY0 of the image:
857 
858  geom::Box2I fullBBox = templateImage.getBBox();
859  int maskStartCol = maskBox.getMinX();
860  int maskStartRow = maskBox.getMinY();
861  int maskEndCol = maskBox.getMaxX();
862  int maskEndRow = maskBox.getMaxY();
863 
864 
865  If we are going to be doing the slicing in Eigen matrices derived from
866  the images, ignore the XY0.
867 
868  geom::Box2I fullBBox = templateImage.getBBox(afwImage::LOCAL);
869 
870  int maskStartCol = maskBox.getMinX() - templateImage.getX0();
871  int maskStartRow = maskBox.getMinY() - templateImage.getY0();
872  int maskEndCol = maskBox.getMaxX() - templateImage.getX0();
873  int maskEndRow = maskBox.getMaxY() - templateImage.getY0();
874 
875  */
876 
877 
878  geom::Box2I shrunkBBox = (*kiter)->shrinkBBox(templateImage.getBBox());
879 
880  LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
881  "Limits of good pixels after convolution: %d,%d -> %d,%d",
882  shrunkBBox.getMinX(), shrunkBBox.getMinY(),
883  shrunkBBox.getMaxX(), shrunkBBox.getMaxY());
884 
885  unsigned int const startCol = shrunkBBox.getMinX();
886  unsigned int const startRow = shrunkBBox.getMinY();
887  unsigned int const endCol = shrunkBBox.getMaxX();
888  unsigned int const endRow = shrunkBBox.getMaxY();
889 
890  /* NOTE: no endCol/endRow += 1 for block slicing, since we are doing the
891  slicing using Afw, not Eigen
892 
893  Eigen arrays have index 0,0 in the upper right, while LSST images
894  have 0,0 in the lower left. The y-axis is flipped. When translating
895  Images to Eigen matrices in ipDiffim::imageToEigenMatrix this is
896  accounted for. However, we still need to be aware of this fact if
897  addressing subregions of an Eigen matrix. This is why the slicing is
898  done in Afw, its cleaner.
899 
900  Please see examples/maskedKernel.cc for elaboration on some of the
901  tests done to make sure this indexing gets done right.
902 
903  */
904 
905 
906  /* Inner limits; of mask box */
907  int maskStartCol = maskBox.getMinX();
908  int maskStartRow = maskBox.getMinY();
909  int maskEndCol = maskBox.getMaxX();
910  int maskEndRow = maskBox.getMaxY();
911 
912  /*
913 
914  |---------------------------|
915  | Kernel Boundary |
916  | |---------------------| |
917  | | Top | |
918  | |......_________......| |
919  | | | | | |
920  | | L | Box | R | |
921  | | | | | |
922  | |......---------......| |
923  | | Bottom | |
924  | |---------------------| |
925  | |
926  |---------------------------|
927 
928  4 regions we want to extract from the pixels: top bottom left right
929 
930  */
931  geom::Box2I tBox = geom::Box2I(geom::Point2I(startCol, maskEndRow + 1),
932  geom::Point2I(endCol, endRow));
933 
934  geom::Box2I bBox = geom::Box2I(geom::Point2I(startCol, startRow),
935  geom::Point2I(endCol, maskStartRow - 1));
936 
937  geom::Box2I lBox = geom::Box2I(geom::Point2I(startCol, maskStartRow),
938  geom::Point2I(maskStartCol - 1, maskEndRow));
939 
940  geom::Box2I rBox = geom::Box2I(geom::Point2I(maskEndCol + 1, maskStartRow),
941  geom::Point2I(endCol, maskEndRow));
942 
943  LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
944  "Upper good pixel region: %d,%d -> %d,%d",
945  tBox.getMinX(), tBox.getMinY(), tBox.getMaxX(), tBox.getMaxY());
946  LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
947  "Bottom good pixel region: %d,%d -> %d,%d",
948  bBox.getMinX(), bBox.getMinY(), bBox.getMaxX(), bBox.getMaxY());
949  LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
950  "Left good pixel region: %d,%d -> %d,%d",
951  lBox.getMinX(), lBox.getMinY(), lBox.getMaxX(), lBox.getMaxY());
952  LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
953  "Right good pixel region: %d,%d -> %d,%d",
954  rBox.getMinX(), rBox.getMinY(), rBox.getMaxX(), rBox.getMaxY());
955 
956  std::vector<geom::Box2I> boxArray;
957  boxArray.push_back(tBox);
958  boxArray.push_back(bBox);
959  boxArray.push_back(lBox);
960  boxArray.push_back(rBox);
961 
962  int totalSize = tBox.getWidth() * tBox.getHeight();
963  totalSize += bBox.getWidth() * bBox.getHeight();
964  totalSize += lBox.getWidth() * lBox.getHeight();
965  totalSize += rBox.getWidth() * rBox.getHeight();
966 
967  Eigen::MatrixXd eigenTemplate(totalSize, 1);
968  Eigen::MatrixXd eigenScience(totalSize, 1);
969  Eigen::MatrixXd eigeniVariance(totalSize, 1);
970  eigenTemplate.setZero();
971  eigenScience.setZero();
972  eigeniVariance.setZero();
973 
974  boost::timer t;
975  t.restart();
976 
977  int nTerms = 0;
978  typename std::vector<geom::Box2I>::iterator biter = boxArray.begin();
979  for (; biter != boxArray.end(); ++biter) {
980  int area = (*biter).getWidth() * (*biter).getHeight();
981 
982  afwImage::Image<InputT> siTemplate(templateImage, *biter);
983  afwImage::Image<InputT> siScience(scienceImage, *biter);
984  afwImage::Image<InputT> sVarEstimate(varianceEstimate, *biter);
985 
986  Eigen::MatrixXd eTemplate = imageToEigenMatrix(siTemplate);
987  Eigen::MatrixXd eScience = imageToEigenMatrix(siScience);
988  Eigen::MatrixXd eiVarEstimate = imageToEigenMatrix(sVarEstimate).array().inverse().matrix();
989 
990  eTemplate.resize(area, 1);
991  eScience.resize(area, 1);
992  eiVarEstimate.resize(area, 1);
993 
994  eigenTemplate.block(nTerms, 0, area, 1) = eTemplate.block(0, 0, area, 1);
995  eigenScience.block(nTerms, 0, area, 1) = eScience.block(0, 0, area, 1);
996  eigeniVariance.block(nTerms, 0, area, 1) = eiVarEstimate.block(0, 0, area, 1);
997 
998  nTerms += area;
999  }
1000 
1001  afwImage::Image<InputT> cimage(templateImage.getDimensions());
1002 
1003  std::vector<Eigen::MatrixXd> convolvedEigenList(nKernelParameters);
1004  typename std::vector<Eigen::MatrixXd>::iterator eiter = convolvedEigenList.begin();
1005  /* Create C_i in the formalism of Alard & Lupton */
1006  for (kiter = basisList.begin(); kiter != basisList.end(); ++kiter, ++eiter) {
1007  afwMath::convolve(cimage, templateImage, **kiter, false); /* cimage stores convolved image */
1008  Eigen::MatrixXd cMat(totalSize, 1);
1009  cMat.setZero();
1010 
1011  int nTerms = 0;
1012  typename std::vector<geom::Box2I>::iterator biter = boxArray.begin();
1013  for (; biter != boxArray.end(); ++biter) {
1014  int area = (*biter).getWidth() * (*biter).getHeight();
1015 
1016  afwImage::Image<InputT> csubimage(cimage, *biter);
1017  Eigen::MatrixXd esubimage = imageToEigenMatrix(csubimage);
1018  esubimage.resize(area, 1);
1019  cMat.block(nTerms, 0, area, 1) = esubimage.block(0, 0, area, 1);
1020 
1021  nTerms += area;
1022  }
1023 
1024  *eiter = cMat;
1025 
1026  }
1027 
1028  double time = t.elapsed();
1029  LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
1030  "Total compute time to do basis convolutions : %.2f s", time);
1031  t.restart();
1032 
1033  /*
1034  Load matrix with all values from convolvedEigenList : all images
1035  (eigeniVariance, convolvedEigenList) must be the same size
1036  */
1037  Eigen::MatrixXd cMat(eigenTemplate.col(0).size(), nParameters);
1038  typename std::vector<Eigen::MatrixXd>::iterator eiterj = convolvedEigenList.begin();
1039  typename std::vector<Eigen::MatrixXd>::iterator eiterE = convolvedEigenList.end();
1040  for (unsigned int kidxj = 0; eiterj != eiterE; eiterj++, kidxj++) {
1041  cMat.col(kidxj) = eiterj->col(0);
1042  }
1043  /* Treat the last "image" as all 1's to do the background calculation. */
1044  if (this->_fitForBackground)
1045  cMat.col(nParameters-1).fill(1.);
1046 
1047  this->_cMat = cMat;
1048  this->_ivVec = eigeniVariance.col(0);
1049  this->_iVec = eigenScience.col(0);
1050 
1051  /* Make these outside of solve() so I can check condition number */
1052  this->_mMat = this->_cMat.transpose() * this->_ivVec.asDiagonal() * this->_cMat;
1053  this->_bVec = this->_cMat.transpose() * this->_ivVec.asDiagonal() * this->_iVec;
1054  }
1055  /*******************************************************************************************************/
1056 
1057 
1058  template <typename InputT>
1060  lsst::afw::math::KernelList const& basisList,
1061  bool fitForBackground,
1062  Eigen::MatrixXd const& hMat,
1064  )
1065  :
1066  StaticKernelSolution<InputT>(basisList, fitForBackground),
1067  _hMat(hMat),
1068  _ps(ps.deepCopy())
1069  {};
1070 
1071  template <typename InputT>
1073  Eigen::MatrixXd vMat = this->_cMat.jacobiSvd().matrixV();
1074  Eigen::MatrixXd vMatvMatT = vMat * vMat.transpose();
1075 
1076  /* Find pseudo inverse of mMat, which may be ill conditioned */
1077  Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eVecValues(this->_mMat);
1078  Eigen::MatrixXd const& rMat = eVecValues.eigenvectors();
1079  Eigen::VectorXd eValues = eVecValues.eigenvalues();
1080  double eMax = eValues.maxCoeff();
1081  for (int i = 0; i != eValues.rows(); ++i) {
1082  if (eValues(i) == 0.0) {
1083  eValues(i) = 0.0;
1084  }
1085  else if ((eMax / eValues(i)) > maxCond) {
1086  LOGL_DEBUG("TRACE3.ip.diffim.RegularizedKernelSolution.estimateRisk",
1087  "Truncating eValue %d; %.5e / %.5e = %.5e vs. %.5e",
1088  i, eMax, eValues(i), eMax / eValues(i), maxCond);
1089  eValues(i) = 0.0;
1090  }
1091  else {
1092  eValues(i) = 1.0 / eValues(i);
1093  }
1094  }
1095  Eigen::MatrixXd mInv = rMat * eValues.asDiagonal() * rMat.transpose();
1096 
1097  std::vector<double> lambdas = _createLambdaSteps();
1098  std::vector<double> risks;
1099  for (unsigned int i = 0; i < lambdas.size(); i++) {
1100  double l = lambdas[i];
1101  Eigen::MatrixXd mLambda = this->_mMat + l * _hMat;
1102 
1103  try {
1104  KernelSolution::solve(mLambda, this->_bVec);
1105  } catch (pexExcept::Exception &e) {
1106  LSST_EXCEPT_ADD(e, "Unable to solve regularized kernel matrix");
1107  throw e;
1108  }
1109  Eigen::VectorXd term1 = (this->_aVec.transpose() * vMatvMatT * this->_aVec);
1110  if (term1.size() != 1)
1111  throw LSST_EXCEPT(pexExcept::Exception, "Matrix size mismatch");
1112 
1113  double term2a = (vMatvMatT * mLambda.inverse()).trace();
1114 
1115  Eigen::VectorXd term2b = (this->_aVec.transpose() * (mInv * this->_bVec));
1116  if (term2b.size() != 1)
1117  throw LSST_EXCEPT(pexExcept::Exception, "Matrix size mismatch");
1118 
1119  double risk = term1(0) + 2 * (term2a - term2b(0));
1120  LOGL_DEBUG("TRACE4.ip.diffim.RegularizedKernelSolution.estimateRisk",
1121  "Lambda = %.3f, Risk = %.5e",
1122  l, risk);
1123  LOGL_DEBUG("TRACE5.ip.diffim.RegularizedKernelSolution.estimateRisk",
1124  "%.5e + 2 * (%.5e - %.5e)",
1125  term1(0), term2a, term2b(0));
1126  risks.push_back(risk);
1127  }
1128  std::vector<double>::iterator it = min_element(risks.begin(), risks.end());
1129  int index = distance(risks.begin(), it);
1130  LOGL_DEBUG("TRACE3.ip.diffim.RegularizedKernelSolution.estimateRisk",
1131  "Minimum Risk = %.3e at lambda = %.3e", risks[index], lambdas[index]);
1132 
1133  return lambdas[index];
1134 
1135  }
1136 
1137  template <typename InputT>
1138  Eigen::MatrixXd RegularizedKernelSolution<InputT>::getM(bool includeHmat) {
1139  if (includeHmat == true) {
1140  return this->_mMat + _lambda * _hMat;
1141  }
1142  else {
1143  return this->_mMat;
1144  }
1145  }
1146 
1147  template <typename InputT>
1149 
1150  LOGL_DEBUG("TRACE3.ip.diffim.RegularizedKernelSolution.solve",
1151  "cMat is %d x %d; vVec is %d; iVec is %d; hMat is %d x %d",
1152  this->_cMat.rows(), this->_cMat.cols(), this->_ivVec.size(),
1153  this->_iVec.size(), _hMat.rows(), _hMat.cols());
1154 
1155  if (DEBUG_MATRIX2) {
1156  std::cout << "ID: " << (this->_id) << std::endl;
1157  std::cout << "C:" << std::endl;
1158  std::cout << this->_cMat << std::endl;
1159  std::cout << "Sigma^{-1}:" << std::endl;
1160  std::cout << Eigen::MatrixXd(this->_ivVec.asDiagonal()) << std::endl;
1161  std::cout << "Y:" << std::endl;
1162  std::cout << this->_iVec << std::endl;
1163  std::cout << "H:" << std::endl;
1164  std::cout << _hMat << std::endl;
1165  }
1166 
1167 
1168  this->_mMat = this->_cMat.transpose() * this->_ivVec.asDiagonal() * this->_cMat;
1169  this->_bVec = this->_cMat.transpose() * this->_ivVec.asDiagonal() * this->_iVec;
1170 
1171 
1172  /* See N.R. 18.5
1173 
1174  Matrix equation to solve is Y = C a + N
1175  Y = vectorized version of I (I = image to not convolve)
1176  C_i = K_i (x) R (R = image to convolve)
1177  a = kernel coefficients
1178  N = noise
1179 
1180  If we reweight everything by the inverse square root of the noise
1181  covariance, we get a linear model with the identity matrix for
1182  the noise. The problem can then be solved using least squares,
1183  with the normal equations
1184 
1185  C^T Y = C^T C a
1186 
1187  or
1188 
1189  b = M a
1190 
1191  with
1192 
1193  b = C^T Y
1194  M = C^T C
1195  a = (C^T C)^{-1} C^T Y
1196 
1197 
1198  We can regularize the least square problem
1199 
1200  Y = C a + lambda a^T H a (+ N, which can be ignored)
1201 
1202  or the normal equations
1203 
1204  (C^T C + lambda H) a = C^T Y
1205 
1206 
1207  The solution to the regularization of the least squares problem is
1208 
1209  a = (C^T C + lambda H)^{-1} C^T Y
1210 
1211  The approximation to Y is
1212 
1213  C a = C (C^T C + lambda H)^{-1} C^T Y
1214 
1215  with smoothing matrix
1216 
1217  S = C (C^T C + lambda H)^{-1} C^T
1218 
1219  */
1220 
1221  std::string lambdaType = _ps->getAsString("lambdaType");
1222  if (lambdaType == "absolute") {
1223  _lambda = _ps->getAsDouble("lambdaValue");
1224  }
1225  else if (lambdaType == "relative") {
1226  _lambda = this->_mMat.trace() / this->_hMat.trace();
1227  _lambda *= _ps->getAsDouble("lambdaScaling");
1228  }
1229  else if (lambdaType == "minimizeBiasedRisk") {
1230  double tol = _ps->getAsDouble("maxConditionNumber");
1231  _lambda = estimateRisk(tol);
1232  }
1233  else if (lambdaType == "minimizeUnbiasedRisk") {
1235  }
1236  else {
1237  throw LSST_EXCEPT(pexExcept::Exception, "lambdaType in PropertySet not recognized");
1238  }
1239 
1240  LOGL_DEBUG("TRACE3.ip.diffim.RegularizedKernelSolution.solve",
1241  "Applying kernel regularization with lambda = %.2e", _lambda);
1242 
1243 
1244  try {
1245  KernelSolution::solve(this->_mMat + _lambda * _hMat, this->_bVec);
1246  } catch (pexExcept::Exception &e) {
1247  LSST_EXCEPT_ADD(e, "Unable to solve static kernel matrix");
1248  throw e;
1249  }
1250  /* Turn matrices into _kernel and _background */
1252  }
1253 
1254  template <typename InputT>
1256  std::vector<double> lambdas;
1257 
1258  std::string lambdaStepType = _ps->getAsString("lambdaStepType");
1259  if (lambdaStepType == "linear") {
1260  double lambdaLinMin = _ps->getAsDouble("lambdaLinMin");
1261  double lambdaLinMax = _ps->getAsDouble("lambdaLinMax");
1262  double lambdaLinStep = _ps->getAsDouble("lambdaLinStep");
1263  for (double l = lambdaLinMin; l <= lambdaLinMax; l += lambdaLinStep) {
1264  lambdas.push_back(l);
1265  }
1266  }
1267  else if (lambdaStepType == "log") {
1268  double lambdaLogMin = _ps->getAsDouble("lambdaLogMin");
1269  double lambdaLogMax = _ps->getAsDouble("lambdaLogMax");
1270  double lambdaLogStep = _ps->getAsDouble("lambdaLogStep");
1271  for (double l = lambdaLogMin; l <= lambdaLogMax; l += lambdaLogStep) {
1272  lambdas.push_back(pow(10, l));
1273  }
1274  }
1275  else {
1276  throw LSST_EXCEPT(pexExcept::Exception, "lambdaStepType in PropertySet not recognized");
1277  }
1278  return lambdas;
1279  }
1280 
1281  /*******************************************************************************************************/
1282 
1284  lsst::afw::math::KernelList const& basisList,
1285  lsst::afw::math::Kernel::SpatialFunctionPtr spatialKernelFunction,
1288  ) :
1289  KernelSolution(),
1290  _spatialKernelFunction(spatialKernelFunction),
1291  _constantFirstTerm(false),
1292  _kernel(),
1293  _background(background),
1294  _kSum(0.0),
1295  _ps(ps.deepCopy()),
1296  _nbases(0),
1297  _nkt(0),
1298  _nbt(0),
1299  _nt(0) {
1300 
1301  bool isAlardLupton = _ps->getAsString("kernelBasisSet") == "alard-lupton";
1302  bool usePca = _ps->getAsBool("usePcaForSpatialKernel");
1303  if (isAlardLupton || usePca) {
1304  _constantFirstTerm = true;
1305  }
1306  this->_fitForBackground = _ps->getAsBool("fitForBackground");
1307 
1308  _nbases = basisList.size();
1309  _nkt = _spatialKernelFunction->getParameters().size();
1310  _nbt = _fitForBackground ? _background->getParameters().size() : 0;
1311  _nt = 0;
1312  if (_constantFirstTerm) {
1313  _nt = (_nbases - 1) * _nkt + 1 + _nbt;
1314  } else {
1315  _nt = _nbases * _nkt + _nbt;
1316  }
1317 
1318  Eigen::MatrixXd mMat(_nt, _nt);
1319  Eigen::VectorXd bVec(_nt);
1320  mMat.setZero();
1321  bVec.setZero();
1322 
1323  this->_mMat = mMat;
1324  this->_bVec = bVec;
1325 
1327  new afwMath::LinearCombinationKernel(basisList, *_spatialKernelFunction)
1328  );
1329 
1330  LOGL_DEBUG("TRACE3.ip.diffim.SpatialKernelSolution",
1331  "Initializing with size %d %d %d and constant first term = %s",
1332  _nkt, _nbt, _nt,
1333  _constantFirstTerm ? "true" : "false");
1334 
1335  }
1336 
1337  void SpatialKernelSolution::addConstraint(float xCenter, float yCenter,
1338  Eigen::MatrixXd const& qMat,
1339  Eigen::VectorXd const& wVec) {
1340 
1341  LOGL_DEBUG("TRACE5.ip.diffim.SpatialKernelSolution.addConstraint",
1342  "Adding candidate at %f, %f", xCenter, yCenter);
1343 
1344  /* Calculate P matrices */
1345  /* Pure kernel terms */
1346  Eigen::VectorXd pK(_nkt);
1347  std::vector<double> paramsK = _spatialKernelFunction->getParameters();
1348  for (int idx = 0; idx < _nkt; idx++) { paramsK[idx] = 0.0; }
1349  for (int idx = 0; idx < _nkt; idx++) {
1350  paramsK[idx] = 1.0;
1351  _spatialKernelFunction->setParameters(paramsK);
1352  pK(idx) = (*_spatialKernelFunction)(xCenter, yCenter); /* Assume things don't vary over stamp */
1353  paramsK[idx] = 0.0;
1354  }
1355  Eigen::MatrixXd pKpKt = (pK * pK.transpose());
1356 
1357  Eigen::VectorXd pB;
1358  Eigen::MatrixXd pBpBt;
1359  Eigen::MatrixXd pKpBt;
1360  if (_fitForBackground) {
1361  pB = Eigen::VectorXd(_nbt);
1362 
1363  /* Pure background terms */
1364  std::vector<double> paramsB = _background->getParameters();
1365  for (int idx = 0; idx < _nbt; idx++) { paramsB[idx] = 0.0; }
1366  for (int idx = 0; idx < _nbt; idx++) {
1367  paramsB[idx] = 1.0;
1368  _background->setParameters(paramsB);
1369  pB(idx) = (*_background)(xCenter, yCenter); /* Assume things don't vary over stamp */
1370  paramsB[idx] = 0.0;
1371  }
1372  pBpBt = (pB * pB.transpose());
1373 
1374  /* Cross terms */
1375  pKpBt = (pK * pB.transpose());
1376  }
1377 
1378  if (DEBUG_MATRIX) {
1379  std::cout << "Spatial weights" << std::endl;
1380  std::cout << "pKpKt " << pKpKt << std::endl;
1381  if (_fitForBackground) {
1382  std::cout << "pBpBt " << pBpBt << std::endl;
1383  std::cout << "pKpBt " << pKpBt << std::endl;
1384  }
1385  }
1386 
1387  if (DEBUG_MATRIX) {
1388  std::cout << "Spatial matrix inputs" << std::endl;
1389  std::cout << "M " << qMat << std::endl;
1390  std::cout << "B " << wVec << std::endl;
1391  }
1392 
1393  /* first index to start the spatial blocks; default=0 for non-constant first term */
1394  int m0 = 0;
1395  /* how many rows/cols to adjust the matrices/vectors; default=0 for non-constant first term */
1396  int dm = 0;
1397  /* where to start the background terms; this is always true */
1398  int mb = _nt - _nbt;
1399 
1400  if (_constantFirstTerm) {
1401  m0 = 1; /* we need to manually fill in the first (non-spatial) terms below */
1402  dm = _nkt-1; /* need to shift terms due to lack of spatial variation in first term */
1403 
1404  _mMat(0, 0) += qMat(0,0);
1405  for(int m2 = 1; m2 < _nbases; m2++) {
1406  _mMat.block(0, m2*_nkt-dm, 1, _nkt) += qMat(0,m2) * pK.transpose();
1407  }
1408  _bVec(0) += wVec(0);
1409 
1410  if (_fitForBackground) {
1411  _mMat.block(0, mb, 1, _nbt) += qMat(0,_nbases) * pB.transpose();
1412  }
1413  }
1414 
1415  /* Fill in the spatial blocks */
1416  for(int m1 = m0; m1 < _nbases; m1++) {
1417  /* Diagonal kernel-kernel term; only use upper triangular part of pKpKt */
1418  _mMat.block(m1*_nkt-dm, m1*_nkt-dm, _nkt, _nkt) +=
1419  (pKpKt * qMat(m1,m1)).triangularView<Eigen::Upper>();
1420 
1421  /* Kernel-kernel terms */
1422  for(int m2 = m1+1; m2 < _nbases; m2++) {
1423  _mMat.block(m1*_nkt-dm, m2*_nkt-dm, _nkt, _nkt) += qMat(m1,m2) * pKpKt;
1424  }
1425 
1426  if (_fitForBackground) {
1427  /* Kernel cross terms with background */
1428  _mMat.block(m1*_nkt-dm, mb, _nkt, _nbt) += qMat(m1,_nbases) * pKpBt;
1429  }
1430 
1431  /* B vector */
1432  _bVec.segment(m1*_nkt-dm, _nkt) += wVec(m1) * pK;
1433  }
1434 
1435  if (_fitForBackground) {
1436  /* Background-background terms only */
1437  _mMat.block(mb, mb, _nbt, _nbt) +=
1438  (pBpBt * qMat(_nbases,_nbases)).triangularView<Eigen::Upper>();
1439  _bVec.segment(mb, _nbt) += wVec(_nbases) * pB;
1440  }
1441 
1442  if (DEBUG_MATRIX) {
1443  std::cout << "Spatial matrix outputs" << std::endl;
1444  std::cout << "mMat " << _mMat << std::endl;
1445  std::cout << "bVec " << _bVec << std::endl;
1446  }
1447 
1448  }
1449 
1452  throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return image");
1453  }
1455  new afwImage::Image<afwMath::Kernel::Pixel>(_kernel->getDimensions())
1456  );
1457  (void)_kernel->computeImage(*image, false, pos[0], pos[1]);
1458  return image;
1459  }
1460 
1462  /* Fill in the other half of mMat */
1463  for (int i = 0; i < _nt; i++) {
1464  for (int j = i+1; j < _nt; j++) {
1465  _mMat(j,i) = _mMat(i,j);
1466  }
1467  }
1468 
1469  try {
1471  } catch (pexExcept::Exception &e) {
1472  LSST_EXCEPT_ADD(e, "Unable to solve spatial kernel matrix");
1473  throw e;
1474  }
1475  /* Turn matrices into _kernel and _background */
1476  _setKernel();
1477  }
1478 
1482  throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return solution");
1483  }
1484 
1485  return std::make_pair(_kernel, _background);
1486  }
1487 
1488  void SpatialKernelSolution::_setKernel() {
1489  /* Report on condition number */
1490  double cNumber = this->getConditionNumber(EIGENVALUE);
1491 
1492  if (_nkt == 1) {
1493  /* Not spatially varying; this fork is a specialization for convolution speed--up */
1494 
1495  /* Set the basis coefficients */
1496  std::vector<double> kCoeffs(_nbases);
1497  for (int i = 0; i < _nbases; i++) {
1498  if (std::isnan(_aVec(i))) {
1499  throw LSST_EXCEPT(
1501  str(boost::format(
1502  "I. Unable to determine spatial kernel solution %d (nan). Condition number = %.3e") % i % cNumber));
1503  }
1504  kCoeffs[i] = _aVec(i);
1505  }
1506  lsst::afw::math::KernelList basisList =
1508  _kernel.reset(
1509  new afwMath::LinearCombinationKernel(basisList, kCoeffs)
1510  );
1511  }
1512  else {
1513 
1514  /* Set the kernel coefficients */
1516  kCoeffs.reserve(_nbases);
1517  for (int i = 0, idx = 0; i < _nbases; i++) {
1518  kCoeffs.push_back(std::vector<double>(_nkt));
1519 
1520  /* Deal with the possibility the first term doesn't vary spatially */
1521  if ((i == 0) && (_constantFirstTerm)) {
1522  if (std::isnan(_aVec(idx))) {
1523  throw LSST_EXCEPT(
1525  str(boost::format(
1526  "II. Unable to determine spatial kernel solution %d (nan). Condition number = %.3e") % idx % cNumber));
1527  }
1528  kCoeffs[i][0] = _aVec(idx++);
1529  }
1530  else {
1531  for (int j = 0; j < _nkt; j++) {
1532  if (std::isnan(_aVec(idx))) {
1533  throw LSST_EXCEPT(
1535  str(boost::format(
1536  "III. Unable to determine spatial kernel solution %d (nan). Condition number = %.3e") % idx % cNumber));
1537  }
1538  kCoeffs[i][j] = _aVec(idx++);
1539  }
1540  }
1541  }
1542  _kernel->setSpatialParameters(kCoeffs);
1543  }
1544 
1545  /* Set kernel Sum */
1546  std::shared_ptr<ImageT> image (new ImageT(_kernel->getDimensions()));
1547  _kSum = _kernel->computeImage(*image, false);
1548 
1549  /* Set the background coefficients */
1550  std::vector<double> bgCoeffs(_fitForBackground ? _nbt : 1);
1551  if (_fitForBackground) {
1552  for (int i = 0; i < _nbt; i++) {
1553  int idx = _nt - _nbt + i;
1554  if (std::isnan(_aVec(idx))) {
1556  str(boost::format(
1557  "Unable to determine spatial background solution %d (nan)") % (idx)));
1558  }
1559  bgCoeffs[i] = _aVec(idx);
1560  }
1561  }
1562  else {
1563  bgCoeffs[0] = 0.;
1564  }
1565  _background->setParameters(bgCoeffs);
1566  }
1567 
1568 /***********************************************************************************************************/
1569 //
1570 // Explicit instantiations
1571 //
1572  typedef float InputT;
1573 
1574  template class StaticKernelSolution<InputT>;
1575  template class MaskedKernelSolution<InputT>;
1576  template class RegularizedKernelSolution<InputT>;
1577 
1578 }}} // end of namespace lsst::ip::diffim
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
uint64_t * ptr
Definition: RangeSet.cc:88
RegularizedKernelSolution(lsst::afw::math::KernelList const &basisList, bool fitForBackground, Eigen::MatrixXd const &hMat, lsst::daf::base::PropertySet const &ps)
afw::table::Key< afw::table::Array< ImagePixelT > > image
int getHeight() const noexcept
Definition: Box.h:188
Eigen::MatrixXi maskToEigenMatrix(lsst::afw::image::Mask< lsst::afw::image::MaskPixel > const &mask)
estimate sample minimum
Definition: Statistics.h:76
virtual std::pair< std::shared_ptr< lsst::afw::math::Kernel >, double > getSolutionPair()
Eigen::VectorXd _bVec
Derived least squares B vector.
T endl(T... args)
virtual std::shared_ptr< lsst::afw::image::Image< lsst::afw::math::Kernel::Pixel > > makeKernelImage()
std::shared_ptr< lsst::afw::math::Function2< double > > SpatialFunctionPtr
Definition: Kernel.h:114
Use (pixels & (given mask))
Definition: Threshold.h:48
A Threshold is used to pass a threshold value to detection algorithms.
Definition: Threshold.h:43
T end(T... args)
KernelSolvedBy _solvedBy
Type of algorithm used to make solution.
Provides consistent interface for LSST exceptions.
Definition: Exception.h:107
StaticKernelSolution(lsst::afw::math::KernelList const &basisList, bool fitForBackground)
#define DEBUG_MATRIX2
double _background
Derived differential background estimate.
Eigen::VectorXd _iVec
Vectorized I.
#define LOGL_DEBUG(logger, message...)
Log a debug-level message using a varargs/printf style interface.
Definition: Log.h:504
Eigen::VectorXd _ivVec
Inverse variance.
A class to evaluate image statistics.
Definition: Statistics.h:215
std::pair< std::shared_ptr< lsst::afw::math::LinearCombinationKernel >, lsst::afw::math::Kernel::SpatialFunctionPtr > getSolutionPair()
Declaration of classes to store the solution for convolution kernels.
iterator end() const
Return an STL compliant iterator to the end of the image.
Definition: Image.cc:268
STL class.
LSST DM logging module built on log4cxx.
T push_back(T... args)
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)
void _setKernel()
Set kernel after solution.
virtual std::shared_ptr< lsst::afw::math::Kernel > getKernel()
A base class for image defects.
static int _SolutionId
Unique identifier for solution.
void _setKernelUncertainty()
Not implemented.
bool _fitForBackground
Background terms included in fit.
lsst::geom::Box2I getBBox(ImageOrigin origin=PARENT) const
Definition: ImageBase.h:462
Represent a 2-dimensional array of bitmask pixels.
Definition: Mask.h:77
lsst::afw::image::Image< lsst::afw::math::Kernel::Pixel > ImageT
int getMaxY() const noexcept
Definition: Box.h:162
virtual double getConditionNumber(ConditionNumberType conditionType)
int _id
Unique ID for object.
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:354
A kernel that is a linear combination of fixed basis kernels.
Definition: Kernel.h:738
T make_pair(T... args)
int getWidth() const noexcept
Definition: Box.h:187
Eigen::MatrixXd imageToEigenMatrix(lsst::afw::image::Image< PixelT > const &img)
Turns a 2-d Image into a 2-d Eigen Matrix.
T reset(T... args)
T dynamic_pointer_cast(T... args)
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)
static MaskPixelT getPlaneBitMask(const std::vector< std::string > &names)
Return the bitmask corresponding to a vector of plane names OR&#39;d together.
Definition: Mask.cc:379
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...
SpatialKernelSolution(lsst::afw::math::KernelList const &basisList, lsst::afw::math::Kernel::SpatialFunctionPtr spatialKernelFunction, lsst::afw::math::Kernel::SpatialFunctionPtr background, lsst::daf::base::PropertySet const &ps)
int getMaxX() const noexcept
Definition: Box.h:161
Class to describe the properties of a detected object from an image.
Definition: Footprint.h:62
double getValue(Property const prop=NOTHING) const
Return the value of the desired property (if specified in the constructor)
Definition: Statistics.cc:1056
Eigen::MatrixXd _mMat
Derived least squares M matrix.
#define DEBUG_MATRIX
int getMinX() const noexcept
Definition: Box.h:157
lsst::geom::Point2I getXY0() const
Return the image&#39;s origin.
Definition: ImageBase.h:340
T size(T... args)
double _kSum
Derived kernel sum.
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
A set of Footprints, associated with a MaskedImage.
Definition: FootprintSet.h:53
std::shared_ptr< lsst::afw::image::Image< lsst::afw::math::Kernel::Pixel > > makeKernelImage(lsst::geom::Point2D const &pos)
MaskedKernelSolution(lsst::afw::math::KernelList const &basisList, bool fitForBackground)
void addConstraint(float xCenter, float yCenter, Eigen::MatrixXd const &qMat, Eigen::VectorXd const &wVec)
T begin(T... args)
Class for storing generic metadata.
Definition: PropertySet.h:67
Reports invalid arguments.
Definition: Runtime.h:66
T isnan(T... args)
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)
iterator begin() const
Return an STL compliant iterator to the start of the image.
Definition: Image.cc:263
Eigen::MatrixXd const & getM()
Eigen::VectorXd _aVec
Derived least squares solution matrix.
std::shared_ptr< lsst::afw::math::Kernel > _kernel
Derived single-object convolution kernel.
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.
lsst::geom::Extent2I getDimensions() const
Return the image&#39;s size; useful for passing to constructors.
Definition: ImageBase.h:373
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects...
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)
std::shared_ptr< FootprintList > getFootprints()
: Return the Footprints of detected objects
Definition: FootprintSet.h:156
#define LSST_EXCEPT_ADD(e, m)
Add the current location and a message to an existing exception before rethrowing it...
Definition: Exception.h:54
An integer coordinate rectangle.
Definition: Box.h:55
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:58
Image Subtraction helper functions.
int end
int getMinY() const noexcept
Definition: Box.h:158
T reserve(T... args)