LSSTApplications  17.0+11,17.0+34,17.0+56,17.0+57,17.0+59,17.0+7,17.0-1-g377950a+33,17.0.1-1-g114240f+2,17.0.1-1-g4d4fbc4+28,17.0.1-1-g55520dc+49,17.0.1-1-g5f4ed7e+52,17.0.1-1-g6dd7d69+17,17.0.1-1-g8de6c91+11,17.0.1-1-gb9095d2+7,17.0.1-1-ge9fec5e+5,17.0.1-1-gf4e0155+55,17.0.1-1-gfc65f5f+50,17.0.1-1-gfc6fb1f+20,17.0.1-10-g87f9f3f+1,17.0.1-11-ge9de802+16,17.0.1-16-ga14f7d5c+4,17.0.1-17-gc79d625+1,17.0.1-17-gdae4c4a+8,17.0.1-2-g26618f5+29,17.0.1-2-g54f2ebc+9,17.0.1-2-gf403422+1,17.0.1-20-g2ca2f74+6,17.0.1-23-gf3eadeb7+1,17.0.1-3-g7e86b59+39,17.0.1-3-gb5ca14a,17.0.1-3-gd08d533+40,17.0.1-30-g596af8797,17.0.1-4-g59d126d+4,17.0.1-4-gc69c472+5,17.0.1-6-g5afd9b9+4,17.0.1-7-g35889ee+1,17.0.1-7-gc7c8782+18,17.0.1-9-gc4bbfb2+3,w.2019.22
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/log/Log.h"
32 
35 
36 #include "ndarray.h"
37 #include "ndarray/eigen.h"
38 
39 #define DEBUG_MATRIX 0
40 #define DEBUG_MATRIX2 0
41 
42 namespace afwDet = lsst::afw::detection;
43 namespace afwMath = lsst::afw::math;
44 namespace afwGeom = lsst::afw::geom;
45 namespace afwImage = lsst::afw::image;
47 
48 namespace lsst {
49 namespace ip {
50 namespace diffim {
51 
52  /* Unique identifier for solution */
54 
56  Eigen::MatrixXd mMat,
57  Eigen::VectorXd bVec,
58  bool fitForBackground
59  ) :
60  _id(++_SolutionId),
61  _mMat(mMat),
62  _bVec(bVec),
63  _aVec(),
64  _solvedBy(NONE),
65  _fitForBackground(fitForBackground)
66  {};
67 
69  bool fitForBackground
70  ) :
71  _id(++_SolutionId),
72  _mMat(),
73  _bVec(),
74  _aVec(),
75  _solvedBy(NONE),
76  _fitForBackground(fitForBackground)
77  {};
78 
80  _id(++_SolutionId),
81  _mMat(),
82  _bVec(),
83  _aVec(),
84  _solvedBy(NONE),
85  _fitForBackground(true)
86  {};
87 
89  solve(_mMat, _bVec);
90  }
91 
93  return getConditionNumber(_mMat, conditionType);
94  }
95 
96  double KernelSolution::getConditionNumber(Eigen::MatrixXd const& mMat,
97  ConditionNumberType conditionType) {
98  switch (conditionType) {
99  case EIGENVALUE:
100  {
101  Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eVecValues(mMat);
102  Eigen::VectorXd eValues = eVecValues.eigenvalues();
103  double eMax = eValues.maxCoeff();
104  double eMin = eValues.minCoeff();
105  LOGL_DEBUG("TRACE3.ip.diffim.KernelSolution.getConditionNumber",
106  "EIGENVALUE eMax / eMin = %.3e", eMax / eMin);
107  return (eMax / eMin);
108  break;
109  }
110  case SVD:
111  {
112  Eigen::VectorXd sValues = mMat.jacobiSvd().singularValues();
113  double sMax = sValues.maxCoeff();
114  double sMin = sValues.minCoeff();
115  LOGL_DEBUG("TRACE3.ip.diffim.KernelSolution.getConditionNumber",
116  "SVD eMax / eMin = %.3e", sMax / sMin);
117  return (sMax / sMin);
118  break;
119  }
120  default:
121  {
123  "Undefined ConditionNumberType : only EIGENVALUE, SVD allowed.");
124  break;
125  }
126  }
127  }
128 
129  void KernelSolution::solve(Eigen::MatrixXd const& mMat,
130  Eigen::VectorXd const& bVec) {
131 
132  if (DEBUG_MATRIX) {
133  std::cout << "M " << std::endl;
134  std::cout << mMat << std::endl;
135  std::cout << "B " << std::endl;
136  std::cout << bVec << std::endl;
137  }
138 
139  Eigen::VectorXd aVec = Eigen::VectorXd::Zero(bVec.size());
140 
141  boost::timer t;
142  t.restart();
143 
144  LOGL_DEBUG("TRACE2.ip.diffim.KernelSolution.solve",
145  "Solving for kernel");
146  _solvedBy = LU;
147  Eigen::FullPivLU<Eigen::MatrixXd> lu(mMat);
148  if (lu.isInvertible()) {
149  aVec = lu.solve(bVec);
150  } else {
151  LOGL_DEBUG("TRACE3.ip.diffim.KernelSolution.solve",
152  "Unable to determine kernel via LU");
153  /* LAST RESORT */
154  try {
155 
157  Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eVecValues(mMat);
158  Eigen::MatrixXd const& rMat = eVecValues.eigenvectors();
159  Eigen::VectorXd eValues = eVecValues.eigenvalues();
160 
161  for (int i = 0; i != eValues.rows(); ++i) {
162  if (eValues(i) != 0.0) {
163  eValues(i) = 1.0/eValues(i);
164  }
165  }
166 
167  aVec = rMat * eValues.asDiagonal() * rMat.transpose() * bVec;
168  } catch (pexExcept::Exception& e) {
169 
170  _solvedBy = NONE;
171  LOGL_DEBUG("TRACE3.ip.diffim.KernelSolution.solve",
172  "Unable to determine kernel via eigen-values");
173 
174  throw LSST_EXCEPT(pexExcept::Exception, "Unable to determine kernel solution");
175  }
176  }
177 
178  double time = t.elapsed();
179  LOGL_DEBUG("TRACE3.ip.diffim.KernelSolution.solve",
180  "Compute time for matrix math : %.2f s", time);
181 
182  if (DEBUG_MATRIX) {
183  std::cout << "A " << std::endl;
184  std::cout << aVec << std::endl;
185  }
186 
187  _aVec = aVec;
188  }
189 
190  /*******************************************************************************************************/
191 
192  template <typename InputT>
194  lsst::afw::math::KernelList const& basisList,
195  bool fitForBackground
196  )
197  :
198  KernelSolution(fitForBackground),
199  _cMat(),
200  _iVec(),
201  _ivVec(),
202  _kernel(),
203  _background(0.0),
204  _kSum(0.0)
205  {
206  std::vector<double> kValues(basisList.size());
208  new afwMath::LinearCombinationKernel(basisList, kValues)
209  );
210  };
211 
212  template <typename InputT>
215  throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return solution");
216  }
217  return _kernel;
218  }
219 
220  template <typename InputT>
223  throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return image");
224  }
227  );
228  (void)_kernel->computeImage(*image, false);
229  return image;
230  }
231 
232  template <typename InputT>
235  throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return background");
236  }
237  return _background;
238  }
239 
240  template <typename InputT>
243  throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return ksum");
244  }
245  return _kSum;
246  }
247 
248  template <typename InputT>
252  throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return solution");
253  }
254 
256  }
257 
258  template <typename InputT>
260  lsst::afw::image::Image<InputT> const &templateImage,
261  lsst::afw::image::Image<InputT> const &scienceImage,
263  ) {
264 
265  afwMath::Statistics varStats = afwMath::makeStatistics(varianceEstimate, afwMath::MIN);
266  if (varStats.getValue(afwMath::MIN) < 0.0) {
268  "Error: variance less than 0.0");
269  }
270  if (varStats.getValue(afwMath::MIN) == 0.0) {
272  "Error: variance equals 0.0, cannot inverse variance weight");
273  }
274 
275  lsst::afw::math::KernelList basisList =
277 
278  unsigned int const nKernelParameters = basisList.size();
279  unsigned int const nBackgroundParameters = _fitForBackground ? 1 : 0;
280  unsigned int const nParameters = nKernelParameters + nBackgroundParameters;
281 
282  std::vector<std::shared_ptr<afwMath::Kernel> >::const_iterator kiter = basisList.begin();
283 
284  /* Ignore buffers around edge of convolved images :
285  *
286  * If the kernel has width 5, it has center pixel 2. The first good pixel
287  * is the (5-2)=3rd pixel, which is array index 2, and ends up being the
288  * index of the central pixel.
289  *
290  * You also have a buffer of unusable pixels on the other side, numbered
291  * width-center-1. The last good usable pixel is N-width+center+1.
292  *
293  * Example : the kernel is width = 5, center = 2
294  *
295  * ---|---|-c-|---|---|
296  *
297  * the image is width = N
298  * convolve this with the kernel, and you get
299  *
300  * |-x-|-x-|-g-|---|---| ... |---|---|-g-|-x-|-x-|
301  *
302  * g = first/last good pixel
303  * x = bad
304  *
305  * the first good pixel is the array index that has the value "center", 2
306  * the last good pixel has array index N-(5-2)+1
307  * eg. if N = 100, you want to use up to index 97
308  * 100-3+1 = 98, and the loops use i < 98, meaning the last
309  * index you address is 97.
310  */
311 
312  /* NOTE - we are accessing particular elements of Eigen arrays using
313  these coordinates, therefore they need to be in LOCAL coordinates.
314  This was written before ndarray unification.
315  */
316  afwGeom::Box2I goodBBox = (*kiter)->shrinkBBox(templateImage.getBBox(afwImage::LOCAL));
317  unsigned int const startCol = goodBBox.getMinX();
318  unsigned int const startRow = goodBBox.getMinY();
319  // endCol/Row is one past the index of the last good col/row
320  unsigned int endCol = goodBBox.getMaxX() + 1;
321  unsigned int endRow = goodBBox.getMaxY() + 1;
322 
323  boost::timer t;
324  t.restart();
325 
326  /* Eigen representation of input images; only the pixels that are unconvolved in cimage below */
327  Eigen::MatrixXd eigenTemplate = imageToEigenMatrix(templateImage).block(startRow,
328  startCol,
329  endRow-startRow,
330  endCol-startCol);
331  Eigen::MatrixXd eigenScience = imageToEigenMatrix(scienceImage).block(startRow,
332  startCol,
333  endRow-startRow,
334  endCol-startCol);
335  Eigen::MatrixXd eigeniVariance = imageToEigenMatrix(varianceEstimate).block(
336  startRow, startCol, endRow-startRow, endCol-startCol
337  ).array().inverse().matrix();
338 
339  /* Resize into 1-D for later usage */
340  eigenTemplate.resize(eigenTemplate.rows()*eigenTemplate.cols(), 1);
341  eigenScience.resize(eigenScience.rows()*eigenScience.cols(), 1);
342  eigeniVariance.resize(eigeniVariance.rows()*eigeniVariance.cols(), 1);
343 
344  /* Holds image convolved with basis function */
345  afwImage::Image<PixelT> cimage(templateImage.getDimensions());
346 
347  /* Holds eigen representation of image convolved with all basis functions */
348  std::vector<Eigen::MatrixXd> convolvedEigenList(nKernelParameters);
349 
350  /* Iterators over convolved image list and basis list */
351  typename std::vector<Eigen::MatrixXd>::iterator eiter = convolvedEigenList.begin();
352  /* Create C_i in the formalism of Alard & Lupton */
353  for (kiter = basisList.begin(); kiter != basisList.end(); ++kiter, ++eiter) {
354  afwMath::convolve(cimage, templateImage, **kiter, false); /* cimage stores convolved image */
355 
356  Eigen::MatrixXd cMat = imageToEigenMatrix(cimage).block(startRow,
357  startCol,
358  endRow-startRow,
359  endCol-startCol);
360  cMat.resize(cMat.size(), 1);
361  *eiter = cMat;
362 
363  }
364 
365  double time = t.elapsed();
366  LOGL_DEBUG("TRACE3.ip.diffim.StaticKernelSolution.build",
367  "Total compute time to do basis convolutions : %.2f s", time);
368  t.restart();
369 
370  /*
371  Load matrix with all values from convolvedEigenList : all images
372  (eigeniVariance, convolvedEigenList) must be the same size
373  */
374  Eigen::MatrixXd cMat(eigenTemplate.col(0).size(), nParameters);
375  typename std::vector<Eigen::MatrixXd>::iterator eiterj = convolvedEigenList.begin();
376  typename std::vector<Eigen::MatrixXd>::iterator eiterE = convolvedEigenList.end();
377  for (unsigned int kidxj = 0; eiterj != eiterE; eiterj++, kidxj++) {
378  cMat.col(kidxj) = eiterj->col(0);
379  }
380  /* Treat the last "image" as all 1's to do the background calculation. */
381  if (_fitForBackground)
382  cMat.col(nParameters-1).fill(1.);
383 
384  _cMat = cMat;
385  _ivVec = eigeniVariance.col(0);
386  _iVec = eigenScience.col(0);
387 
388  /* Make these outside of solve() so I can check condition number */
389  _mMat = _cMat.transpose() * (_ivVec.asDiagonal() * _cMat);
390  _bVec = _cMat.transpose() * (_ivVec.asDiagonal() * _iVec);
391  }
392 
393  template <typename InputT>
395  LOGL_DEBUG("TRACE3.ip.diffim.StaticKernelSolution.solve",
396  "mMat is %d x %d; bVec is %d; cMat is %d x %d; vVec is %d; iVec is %d",
397  _mMat.rows(), _mMat.cols(), _bVec.size(),
398  _cMat.rows(), _cMat.cols(), _ivVec.size(), _iVec.size());
399 
400  if (DEBUG_MATRIX) {
401  std::cout << "C" << std::endl;
402  std::cout << _cMat << std::endl;
403  std::cout << "iV" << std::endl;
404  std::cout << _ivVec << std::endl;
405  std::cout << "I" << std::endl;
406  std::cout << _iVec << std::endl;
407  }
408 
409  try {
411  } catch (pexExcept::Exception &e) {
412  LSST_EXCEPT_ADD(e, "Unable to solve static kernel matrix");
413  throw e;
414  }
415  /* Turn matrices into _kernel and _background */
416  _setKernel();
417  }
418 
419  template <typename InputT>
422  throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot make solution");
423  }
424 
425  unsigned int const nParameters = _aVec.size();
426  unsigned int const nBackgroundParameters = _fitForBackground ? 1 : 0;
427  unsigned int const nKernelParameters =
429  if (nParameters != (nKernelParameters + nBackgroundParameters))
430  throw LSST_EXCEPT(pexExcept::Exception, "Mismatched sizes in kernel solution");
431 
432  /* Fill in the kernel results */
433  std::vector<double> kValues(nKernelParameters);
434  for (unsigned int idx = 0; idx < nKernelParameters; idx++) {
435  if (std::isnan(_aVec(idx))) {
437  str(boost::format("Unable to determine kernel solution %d (nan)") % idx));
438  }
439  kValues[idx] = _aVec(idx);
440  }
441  _kernel->setKernelParameters(kValues);
442 
444  new ImageT(_kernel->getDimensions())
445  );
446  _kSum = _kernel->computeImage(*image, false);
447 
448  if (_fitForBackground) {
449  if (std::isnan(_aVec(nParameters-1))) {
451  str(boost::format("Unable to determine background solution %d (nan)") %
452  (nParameters-1)));
453  }
454  _background = _aVec(nParameters-1);
455  }
456  }
457 
458 
459  template <typename InputT>
461  throw LSST_EXCEPT(pexExcept::Exception, "Uncertainty calculation not supported");
462 
463  /* Estimate of parameter uncertainties comes from the inverse of the
464  * covariance matrix (noise spectrum).
465  * N.R. 15.4.8 to 15.4.15
466  *
467  * Since this is a linear problem no need to use Fisher matrix
468  * N.R. 15.5.8
469  *
470  * Although I might be able to take advantage of the solution above.
471  * Since this now works and is not the rate limiting step, keep as-is for DC3a.
472  *
473  * Use Cholesky decomposition again.
474  * Cholkesy:
475  * Cov = L L^t
476  * Cov^(-1) = (L L^t)^(-1)
477  * = (L^T)^-1 L^(-1)
478  *
479  * Code would be:
480  *
481  * Eigen::MatrixXd cov = _mMat.transpose() * _mMat;
482  * Eigen::LLT<Eigen::MatrixXd> llt = cov.llt();
483  * Eigen::MatrixXd error2 = llt.matrixL().transpose().inverse()*llt.matrixL().inverse();
484  */
485  }
486 
487  /*******************************************************************************************************/
488 
489  template <typename InputT>
491  lsst::afw::math::KernelList const& basisList,
492  bool fitForBackground
493  ) :
494  StaticKernelSolution<InputT>(basisList, fitForBackground)
495  {};
496 
497  template <typename InputT>
499  lsst::afw::image::Image<InputT> const &templateImage,
500  lsst::afw::image::Image<InputT> const &scienceImage,
503  ) {
504 
505  afwMath::Statistics varStats = afwMath::makeStatistics(varianceEstimate, afwMath::MIN);
506  if (varStats.getValue(afwMath::MIN) < 0.0) {
508  "Error: variance less than 0.0");
509  }
510  if (varStats.getValue(afwMath::MIN) == 0.0) {
512  "Error: variance equals 0.0, cannot inverse variance weight");
513  }
514 
515  /* Full footprint of all input images */
517  new afwDet::Footprint(std::make_shared<afwGeom::SpanSet>(templateImage.getBBox())));
518 
519  afwMath::KernelList basisList =
521  std::vector<std::shared_ptr<afwMath::Kernel> >::const_iterator kiter = basisList.begin();
522 
523  /* Only BAD pixels marked in this mask */
524  afwImage::MaskPixel bitMask =
529 
530  /* Create a Footprint that contains all the masked pixels set above */
532  afwDet::FootprintSet maskFpSet(pixelMask, threshold, true);
533 
534  /* And spread it by the kernel half width */
535  int growPix = (*kiter)->getCtr().getX();
536  afwDet::FootprintSet maskedFpSetGrown(maskFpSet, growPix, true);
537 
538 #if 0
539  for (typename afwDet::FootprintSet::FootprintList::iterator
540  ptr = maskedFpSetGrown.getFootprints()->begin(),
541  end = maskedFpSetGrown.getFootprints()->end();
542  ptr != end;
543  ++ptr) {
544 
545  afwDet::setMaskFromFootprint(finalMask,
546  (**ptr),
548  }
549 #endif
550 
551  afwImage::Mask<afwImage::MaskPixel> finalMask(pixelMask.getDimensions());
552  for (auto const & foot : *(maskedFpSetGrown.getFootprints())) {
553  foot->getSpans()->setMask(finalMask, afwImage::Mask<afwImage::MaskPixel>::getPlaneBitMask("BAD"));
554  }
555  pixelMask.writeFits("pixelmask.fits");
556  finalMask.writeFits("finalmask.fits");
557 
558 
559  ndarray::Array<int, 1, 1> maskArray =
560  ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
561  fullFp->getSpans()->flatten(maskArray, finalMask.getArray(), templateImage.getXY0());
562  auto maskEigen = ndarray::asEigenMatrix(maskArray);
563 
564  ndarray::Array<InputT, 1, 1> arrayTemplate =
565  ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
566  fullFp->getSpans()->flatten(arrayTemplate, templateImage.getArray(), templateImage.getXY0());
567  auto eigenTemplate0 = ndarray::asEigenMatrix(arrayTemplate);
568 
569  ndarray::Array<InputT, 1, 1> arrayScience =
570  ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
571  fullFp->getSpans()->flatten(arrayScience, scienceImage.getArray(), scienceImage.getXY0());
572  auto eigenScience0 = ndarray::asEigenMatrix(arrayScience);
573 
574  ndarray::Array<afwImage::VariancePixel, 1, 1> arrayVariance =
575  ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
576  fullFp->getSpans()->flatten(arrayVariance, varianceEstimate.getArray(), varianceEstimate.getXY0());
577  auto eigenVariance0 = ndarray::asEigenMatrix(arrayVariance);
578 
579  int nGood = 0;
580  for (int i = 0; i < maskEigen.size(); i++) {
581  if (maskEigen(i) == 0.0)
582  nGood += 1;
583  }
584 
585  Eigen::VectorXd eigenTemplate(nGood);
586  Eigen::VectorXd eigenScience(nGood);
587  Eigen::VectorXd eigenVariance(nGood);
588  int nUsed = 0;
589  for (int i = 0; i < maskEigen.size(); i++) {
590  if (maskEigen(i) == 0.0) {
591  eigenTemplate(nUsed) = eigenTemplate0(i);
592  eigenScience(nUsed) = eigenScience0(i);
593  eigenVariance(nUsed) = eigenVariance0(i);
594  nUsed += 1;
595  }
596  }
597 
598 
599  boost::timer t;
600  t.restart();
601 
602  unsigned int const nKernelParameters = basisList.size();
603  unsigned int const nBackgroundParameters = this->_fitForBackground ? 1 : 0;
604  unsigned int const nParameters = nKernelParameters + nBackgroundParameters;
605 
606  /* Holds image convolved with basis function */
607  afwImage::Image<InputT> cimage(templateImage.getDimensions());
608 
609  /* Holds eigen representation of image convolved with all basis functions */
610  std::vector<Eigen::VectorXd> convolvedEigenList(nKernelParameters);
611 
612  /* Iterators over convolved image list and basis list */
613  typename std::vector<Eigen::VectorXd>::iterator eiter = convolvedEigenList.begin();
614 
615  /* Create C_i in the formalism of Alard & Lupton */
616  for (kiter = basisList.begin(); kiter != basisList.end(); ++kiter, ++eiter) {
617  afwMath::convolve(cimage, templateImage, **kiter, false); /* cimage stores convolved image */
618 
619  ndarray::Array<InputT, 1, 1> arrayC =
620  ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
621  fullFp->getSpans()->flatten(arrayC, cimage.getArray(), cimage.getXY0());
622  auto eigenC0 = ndarray::asEigenMatrix(arrayC);
623 
624  Eigen::VectorXd eigenC(nGood);
625  int nUsed = 0;
626  for (int i = 0; i < maskEigen.size(); i++) {
627  if (maskEigen(i) == 0.0) {
628  eigenC(nUsed) = eigenC0(i);
629  nUsed += 1;
630  }
631  }
632 
633  *eiter = eigenC;
634  }
635  double time = t.elapsed();
636  LOGL_DEBUG("TRACE3.ip.diffim.StaticKernelSolution.buildWithMask",
637  "Total compute time to do basis convolutions : %.2f s", time);
638  t.restart();
639 
640  /* Load matrix with all convolved images */
641  Eigen::MatrixXd cMat(eigenTemplate.size(), nParameters);
642  typename std::vector<Eigen::VectorXd>::iterator eiterj = convolvedEigenList.begin();
643  typename std::vector<Eigen::VectorXd>::iterator eiterE = convolvedEigenList.end();
644  for (unsigned int kidxj = 0; eiterj != eiterE; eiterj++, kidxj++) {
645  cMat.block(0, kidxj, eigenTemplate.size(), 1) =
646  Eigen::MatrixXd(*eiterj).block(0, 0, eigenTemplate.size(), 1);
647  }
648  /* Treat the last "image" as all 1's to do the background calculation. */
649  if (this->_fitForBackground)
650  cMat.col(nParameters-1).fill(1.);
651 
652  this->_cMat = cMat;
653  this->_ivVec = eigenVariance.array().inverse().matrix();
654  this->_iVec = eigenScience;
655 
656  /* Make these outside of solve() so I can check condition number */
657  this->_mMat = this->_cMat.transpose() * this->_ivVec.asDiagonal() * (this->_cMat);
658  this->_bVec = this->_cMat.transpose() * this->_ivVec.asDiagonal() * (this->_iVec);
659  }
660 
661 
662  template <typename InputT>
664  lsst::afw::image::Image<InputT> const &templateImage,
665  lsst::afw::image::Image<InputT> const &scienceImage,
668  ) {
669 
670  afwMath::Statistics varStats = afwMath::makeStatistics(varianceEstimate, afwMath::MIN);
671  if (varStats.getValue(afwMath::MIN) < 0.0) {
673  "Error: variance less than 0.0");
674  }
675  if (varStats.getValue(afwMath::MIN) == 0.0) {
677  "Error: variance equals 0.0, cannot inverse variance weight");
678  }
679 
680  lsst::afw::math::KernelList basisList =
682  std::vector<std::shared_ptr<afwMath::Kernel> >::const_iterator kiter = basisList.begin();
683 
684  /* Only BAD pixels marked in this mask */
686  afwImage::MaskPixel bitMask =
690  sMask &= bitMask;
691  /* TBD: Need to figure out a way to spread this; currently done in Python */
692 
693  unsigned int const nKernelParameters = basisList.size();
694  unsigned int const nBackgroundParameters = this->_fitForBackground ? 1 : 0;
695  unsigned int const nParameters = nKernelParameters + nBackgroundParameters;
696 
697  /* NOTE - we are accessing particular elements of Eigen arrays using
698  these coordinates, therefore they need to be in LOCAL coordinates.
699  This was written before ndarray unification.
700  */
701  /* Ignore known EDGE pixels for speed */
702  afwGeom::Box2I shrunkLocalBBox = (*kiter)->shrinkBBox(templateImage.getBBox(afwImage::LOCAL));
703  LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
704  "Limits of good pixels after convolution: %d,%d -> %d,%d (local)",
705  shrunkLocalBBox.getMinX(), shrunkLocalBBox.getMinY(),
706  shrunkLocalBBox.getMaxX(), shrunkLocalBBox.getMaxY());
707 
708  /* Subimages are addressed in raw pixel coordinates */
709  unsigned int startCol = shrunkLocalBBox.getMinX();
710  unsigned int startRow = shrunkLocalBBox.getMinY();
711  unsigned int endCol = shrunkLocalBBox.getMaxX();
712  unsigned int endRow = shrunkLocalBBox.getMaxY();
713  /* Needed for Eigen block slicing operation */
714  endCol += 1;
715  endRow += 1;
716 
717  boost::timer t;
718  t.restart();
719 
720  /* Eigen representation of input images; only the pixels that are unconvolved in cimage below */
721  Eigen::MatrixXi eMask = maskToEigenMatrix(sMask).block(startRow,
722  startCol,
723  endRow-startRow,
724  endCol-startCol);
725 
726  Eigen::MatrixXd eigenTemplate = imageToEigenMatrix(templateImage).block(startRow,
727  startCol,
728  endRow-startRow,
729  endCol-startCol);
730  Eigen::MatrixXd eigenScience = imageToEigenMatrix(scienceImage).block(startRow,
731  startCol,
732  endRow-startRow,
733  endCol-startCol);
734  Eigen::MatrixXd eigeniVariance = imageToEigenMatrix(varianceEstimate).block(
735  startRow, startCol, endRow-startRow, endCol-startCol
736  ).array().inverse().matrix();
737 
738  /* Resize into 1-D for later usage */
739  eMask.resize(eMask.rows()*eMask.cols(), 1);
740  eigenTemplate.resize(eigenTemplate.rows()*eigenTemplate.cols(), 1);
741  eigenScience.resize(eigenScience.rows()*eigenScience.cols(), 1);
742  eigeniVariance.resize(eigeniVariance.rows()*eigeniVariance.cols(), 1);
743 
744  /* Do the masking, slowly for now... */
745  Eigen::MatrixXd maskedEigenTemplate(eigenTemplate.rows(), 1);
746  Eigen::MatrixXd maskedEigenScience(eigenScience.rows(), 1);
747  Eigen::MatrixXd maskedEigeniVariance(eigeniVariance.rows(), 1);
748  int nGood = 0;
749  for (int i = 0; i < eMask.rows(); i++) {
750  if (eMask(i, 0) == 0) {
751  maskedEigenTemplate(nGood, 0) = eigenTemplate(i, 0);
752  maskedEigenScience(nGood, 0) = eigenScience(i, 0);
753  maskedEigeniVariance(nGood, 0) = eigeniVariance(i, 0);
754  nGood += 1;
755  }
756  }
757  /* Can't resize() since all values are lost; use blocks */
758  eigenTemplate = maskedEigenTemplate.block(0, 0, nGood, 1);
759  eigenScience = maskedEigenScience.block(0, 0, nGood, 1);
760  eigeniVariance = maskedEigeniVariance.block(0, 0, nGood, 1);
761 
762 
763  /* Holds image convolved with basis function */
764  afwImage::Image<InputT> cimage(templateImage.getDimensions());
765 
766  /* Holds eigen representation of image convolved with all basis functions */
767  std::vector<Eigen::MatrixXd> convolvedEigenList(nKernelParameters);
768 
769  /* Iterators over convolved image list and basis list */
770  typename std::vector<Eigen::MatrixXd>::iterator eiter = convolvedEigenList.begin();
771  /* Create C_i in the formalism of Alard & Lupton */
772  for (kiter = basisList.begin(); kiter != basisList.end(); ++kiter, ++eiter) {
773  afwMath::convolve(cimage, templateImage, **kiter, false); /* cimage stores convolved image */
774 
775  Eigen::MatrixXd cMat = imageToEigenMatrix(cimage).block(startRow,
776  startCol,
777  endRow-startRow,
778  endCol-startCol);
779  cMat.resize(cMat.size(), 1);
780 
781  /* Do masking */
782  Eigen::MatrixXd maskedcMat(cMat.rows(), 1);
783  int nGood = 0;
784  for (int i = 0; i < eMask.rows(); i++) {
785  if (eMask(i, 0) == 0) {
786  maskedcMat(nGood, 0) = cMat(i, 0);
787  nGood += 1;
788  }
789  }
790  cMat = maskedcMat.block(0, 0, nGood, 1);
791  *eiter = cMat;
792  }
793 
794  double time = t.elapsed();
795  LOGL_DEBUG("TRACE3.ip.diffim.StaticKernelSolution.build",
796  "Total compute time to do basis convolutions : %.2f s", time);
797  t.restart();
798 
799  /*
800  Load matrix with all values from convolvedEigenList : all images
801  (eigeniVariance, convolvedEigenList) must be the same size
802  */
803  Eigen::MatrixXd cMat(eigenTemplate.col(0).size(), nParameters);
804  typename std::vector<Eigen::MatrixXd>::iterator eiterj = convolvedEigenList.begin();
805  typename std::vector<Eigen::MatrixXd>::iterator eiterE = convolvedEigenList.end();
806  for (unsigned int kidxj = 0; eiterj != eiterE; eiterj++, kidxj++) {
807  cMat.col(kidxj) = eiterj->col(0);
808  }
809  /* Treat the last "image" as all 1's to do the background calculation. */
810  if (this->_fitForBackground)
811  cMat.col(nParameters-1).fill(1.);
812 
813  this->_cMat = cMat;
814  this->_ivVec = eigeniVariance.col(0);
815  this->_iVec = eigenScience.col(0);
816 
817  /* Make these outside of solve() so I can check condition number */
818  this->_mMat = this->_cMat.transpose() * this->_ivVec.asDiagonal() * this->_cMat;
819  this->_bVec = this->_cMat.transpose() * this->_ivVec.asDiagonal() * this->_iVec;
820 
821  }
822 
823  /* NOTE - this was written before the ndarray unification. I am rewriting
824  buildSingleMask to use the ndarray stuff */
825  template <typename InputT>
827  lsst::afw::image::Image<InputT> const &templateImage,
828  lsst::afw::image::Image<InputT> const &scienceImage,
830  lsst::afw::geom::Box2I maskBox
831  ) {
832 
833  afwMath::Statistics varStats = afwMath::makeStatistics(varianceEstimate, afwMath::MIN);
834  if (varStats.getValue(afwMath::MIN) < 0.0) {
836  "Error: variance less than 0.0");
837  }
838  if (varStats.getValue(afwMath::MIN) == 0.0) {
840  "Error: variance equals 0.0, cannot inverse variance weight");
841  }
842 
843  lsst::afw::math::KernelList basisList =
845 
846  unsigned int const nKernelParameters = basisList.size();
847  unsigned int const nBackgroundParameters = this->_fitForBackground ? 1 : 0;
848  unsigned int const nParameters = nKernelParameters + nBackgroundParameters;
849 
850  std::vector<std::shared_ptr<afwMath::Kernel> >::const_iterator kiter = basisList.begin();
851 
852  /*
853  NOTE : If we are using these views in Afw's Image space, we need to
854  make sure and compensate for the XY0 of the image:
855 
856  afwGeom::Box2I fullBBox = templateImage.getBBox();
857  int maskStartCol = maskBox.getMinX();
858  int maskStartRow = maskBox.getMinY();
859  int maskEndCol = maskBox.getMaxX();
860  int maskEndRow = maskBox.getMaxY();
861 
862 
863  If we are going to be doing the slicing in Eigen matrices derived from
864  the images, ignore the XY0.
865 
866  afwGeom::Box2I fullBBox = templateImage.getBBox(afwImage::LOCAL);
867 
868  int maskStartCol = maskBox.getMinX() - templateImage.getX0();
869  int maskStartRow = maskBox.getMinY() - templateImage.getY0();
870  int maskEndCol = maskBox.getMaxX() - templateImage.getX0();
871  int maskEndRow = maskBox.getMaxY() - templateImage.getY0();
872 
873  */
874 
875 
876  afwGeom::Box2I shrunkBBox = (*kiter)->shrinkBBox(templateImage.getBBox());
877 
878  LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
879  "Limits of good pixels after convolution: %d,%d -> %d,%d",
880  shrunkBBox.getMinX(), shrunkBBox.getMinY(),
881  shrunkBBox.getMaxX(), shrunkBBox.getMaxY());
882 
883  unsigned int const startCol = shrunkBBox.getMinX();
884  unsigned int const startRow = shrunkBBox.getMinY();
885  unsigned int const endCol = shrunkBBox.getMaxX();
886  unsigned int const endRow = shrunkBBox.getMaxY();
887 
888  /* NOTE: no endCol/endRow += 1 for block slicing, since we are doing the
889  slicing using Afw, not Eigen
890 
891  Eigen arrays have index 0,0 in the upper right, while LSST images
892  have 0,0 in the lower left. The y-axis is flipped. When translating
893  Images to Eigen matrices in ipDiffim::imageToEigenMatrix this is
894  accounted for. However, we still need to be aware of this fact if
895  addressing subregions of an Eigen matrix. This is why the slicing is
896  done in Afw, its cleaner.
897 
898  Please see examples/maskedKernel.cc for elaboration on some of the
899  tests done to make sure this indexing gets done right.
900 
901  */
902 
903 
904  /* Inner limits; of mask box */
905  int maskStartCol = maskBox.getMinX();
906  int maskStartRow = maskBox.getMinY();
907  int maskEndCol = maskBox.getMaxX();
908  int maskEndRow = maskBox.getMaxY();
909 
910  /*
911 
912  |---------------------------|
913  | Kernel Boundary |
914  | |---------------------| |
915  | | Top | |
916  | |......_________......| |
917  | | | | | |
918  | | L | Box | R | |
919  | | | | | |
920  | |......---------......| |
921  | | Bottom | |
922  | |---------------------| |
923  | |
924  |---------------------------|
925 
926  4 regions we want to extract from the pixels: top bottom left right
927 
928  */
929  afwGeom::Box2I tBox = afwGeom::Box2I(afwGeom::Point2I(startCol, maskEndRow + 1),
930  afwGeom::Point2I(endCol, endRow));
931 
932  afwGeom::Box2I bBox = afwGeom::Box2I(afwGeom::Point2I(startCol, startRow),
933  afwGeom::Point2I(endCol, maskStartRow - 1));
934 
935  afwGeom::Box2I lBox = afwGeom::Box2I(afwGeom::Point2I(startCol, maskStartRow),
936  afwGeom::Point2I(maskStartCol - 1, maskEndRow));
937 
938  afwGeom::Box2I rBox = afwGeom::Box2I(afwGeom::Point2I(maskEndCol + 1, maskStartRow),
939  afwGeom::Point2I(endCol, maskEndRow));
940 
941  LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
942  "Upper good pixel region: %d,%d -> %d,%d",
943  tBox.getMinX(), tBox.getMinY(), tBox.getMaxX(), tBox.getMaxY());
944  LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
945  "Bottom good pixel region: %d,%d -> %d,%d",
946  bBox.getMinX(), bBox.getMinY(), bBox.getMaxX(), bBox.getMaxY());
947  LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
948  "Left good pixel region: %d,%d -> %d,%d",
949  lBox.getMinX(), lBox.getMinY(), lBox.getMaxX(), lBox.getMaxY());
950  LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
951  "Right good pixel region: %d,%d -> %d,%d",
952  rBox.getMinX(), rBox.getMinY(), rBox.getMaxX(), rBox.getMaxY());
953 
955  boxArray.push_back(tBox);
956  boxArray.push_back(bBox);
957  boxArray.push_back(lBox);
958  boxArray.push_back(rBox);
959 
960  int totalSize = tBox.getWidth() * tBox.getHeight();
961  totalSize += bBox.getWidth() * bBox.getHeight();
962  totalSize += lBox.getWidth() * lBox.getHeight();
963  totalSize += rBox.getWidth() * rBox.getHeight();
964 
965  Eigen::MatrixXd eigenTemplate(totalSize, 1);
966  Eigen::MatrixXd eigenScience(totalSize, 1);
967  Eigen::MatrixXd eigeniVariance(totalSize, 1);
968  eigenTemplate.setZero();
969  eigenScience.setZero();
970  eigeniVariance.setZero();
971 
972  boost::timer t;
973  t.restart();
974 
975  int nTerms = 0;
976  typename std::vector<afwGeom::Box2I>::iterator biter = boxArray.begin();
977  for (; biter != boxArray.end(); ++biter) {
978  int area = (*biter).getWidth() * (*biter).getHeight();
979 
980  afwImage::Image<InputT> siTemplate(templateImage, *biter);
981  afwImage::Image<InputT> siScience(scienceImage, *biter);
982  afwImage::Image<InputT> sVarEstimate(varianceEstimate, *biter);
983 
984  Eigen::MatrixXd eTemplate = imageToEigenMatrix(siTemplate);
985  Eigen::MatrixXd eScience = imageToEigenMatrix(siScience);
986  Eigen::MatrixXd eiVarEstimate = imageToEigenMatrix(sVarEstimate).array().inverse().matrix();
987 
988  eTemplate.resize(area, 1);
989  eScience.resize(area, 1);
990  eiVarEstimate.resize(area, 1);
991 
992  eigenTemplate.block(nTerms, 0, area, 1) = eTemplate.block(0, 0, area, 1);
993  eigenScience.block(nTerms, 0, area, 1) = eScience.block(0, 0, area, 1);
994  eigeniVariance.block(nTerms, 0, area, 1) = eiVarEstimate.block(0, 0, area, 1);
995 
996  nTerms += area;
997  }
998 
999  afwImage::Image<InputT> cimage(templateImage.getDimensions());
1000 
1001  std::vector<Eigen::MatrixXd> convolvedEigenList(nKernelParameters);
1002  typename std::vector<Eigen::MatrixXd>::iterator eiter = convolvedEigenList.begin();
1003  /* Create C_i in the formalism of Alard & Lupton */
1004  for (kiter = basisList.begin(); kiter != basisList.end(); ++kiter, ++eiter) {
1005  afwMath::convolve(cimage, templateImage, **kiter, false); /* cimage stores convolved image */
1006  Eigen::MatrixXd cMat(totalSize, 1);
1007  cMat.setZero();
1008 
1009  int nTerms = 0;
1010  typename std::vector<afwGeom::Box2I>::iterator biter = boxArray.begin();
1011  for (; biter != boxArray.end(); ++biter) {
1012  int area = (*biter).getWidth() * (*biter).getHeight();
1013 
1014  afwImage::Image<InputT> csubimage(cimage, *biter);
1015  Eigen::MatrixXd esubimage = imageToEigenMatrix(csubimage);
1016  esubimage.resize(area, 1);
1017  cMat.block(nTerms, 0, area, 1) = esubimage.block(0, 0, area, 1);
1018 
1019  nTerms += area;
1020  }
1021 
1022  *eiter = cMat;
1023 
1024  }
1025 
1026  double time = t.elapsed();
1027  LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
1028  "Total compute time to do basis convolutions : %.2f s", time);
1029  t.restart();
1030 
1031  /*
1032  Load matrix with all values from convolvedEigenList : all images
1033  (eigeniVariance, convolvedEigenList) must be the same size
1034  */
1035  Eigen::MatrixXd cMat(eigenTemplate.col(0).size(), nParameters);
1036  typename std::vector<Eigen::MatrixXd>::iterator eiterj = convolvedEigenList.begin();
1037  typename std::vector<Eigen::MatrixXd>::iterator eiterE = convolvedEigenList.end();
1038  for (unsigned int kidxj = 0; eiterj != eiterE; eiterj++, kidxj++) {
1039  cMat.col(kidxj) = eiterj->col(0);
1040  }
1041  /* Treat the last "image" as all 1's to do the background calculation. */
1042  if (this->_fitForBackground)
1043  cMat.col(nParameters-1).fill(1.);
1044 
1045  this->_cMat = cMat;
1046  this->_ivVec = eigeniVariance.col(0);
1047  this->_iVec = eigenScience.col(0);
1048 
1049  /* Make these outside of solve() so I can check condition number */
1050  this->_mMat = this->_cMat.transpose() * this->_ivVec.asDiagonal() * this->_cMat;
1051  this->_bVec = this->_cMat.transpose() * this->_ivVec.asDiagonal() * this->_iVec;
1052  }
1053  /*******************************************************************************************************/
1054 
1055 
1056  template <typename InputT>
1058  lsst::afw::math::KernelList const& basisList,
1059  bool fitForBackground,
1060  Eigen::MatrixXd const& hMat,
1062  )
1063  :
1064  StaticKernelSolution<InputT>(basisList, fitForBackground),
1065  _hMat(hMat),
1066  _policy(policy)
1067  {};
1068 
1069  template <typename InputT>
1071  Eigen::MatrixXd vMat = this->_cMat.jacobiSvd().matrixV();
1072  Eigen::MatrixXd vMatvMatT = vMat * vMat.transpose();
1073 
1074  /* Find pseudo inverse of mMat, which may be ill conditioned */
1075  Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eVecValues(this->_mMat);
1076  Eigen::MatrixXd const& rMat = eVecValues.eigenvectors();
1077  Eigen::VectorXd eValues = eVecValues.eigenvalues();
1078  double eMax = eValues.maxCoeff();
1079  for (int i = 0; i != eValues.rows(); ++i) {
1080  if (eValues(i) == 0.0) {
1081  eValues(i) = 0.0;
1082  }
1083  else if ((eMax / eValues(i)) > maxCond) {
1084  LOGL_DEBUG("TRACE3.ip.diffim.RegularizedKernelSolution.estimateRisk",
1085  "Truncating eValue %d; %.5e / %.5e = %.5e vs. %.5e",
1086  i, eMax, eValues(i), eMax / eValues(i), maxCond);
1087  eValues(i) = 0.0;
1088  }
1089  else {
1090  eValues(i) = 1.0 / eValues(i);
1091  }
1092  }
1093  Eigen::MatrixXd mInv = rMat * eValues.asDiagonal() * rMat.transpose();
1094 
1095  std::vector<double> lambdas = _createLambdaSteps();
1096  std::vector<double> risks;
1097  for (unsigned int i = 0; i < lambdas.size(); i++) {
1098  double l = lambdas[i];
1099  Eigen::MatrixXd mLambda = this->_mMat + l * _hMat;
1100 
1101  try {
1102  KernelSolution::solve(mLambda, this->_bVec);
1103  } catch (pexExcept::Exception &e) {
1104  LSST_EXCEPT_ADD(e, "Unable to solve regularized kernel matrix");
1105  throw e;
1106  }
1107  Eigen::VectorXd term1 = (this->_aVec.transpose() * vMatvMatT * this->_aVec);
1108  if (term1.size() != 1)
1109  throw LSST_EXCEPT(pexExcept::Exception, "Matrix size mismatch");
1110 
1111  double term2a = (vMatvMatT * mLambda.inverse()).trace();
1112 
1113  Eigen::VectorXd term2b = (this->_aVec.transpose() * (mInv * this->_bVec));
1114  if (term2b.size() != 1)
1115  throw LSST_EXCEPT(pexExcept::Exception, "Matrix size mismatch");
1116 
1117  double risk = term1(0) + 2 * (term2a - term2b(0));
1118  LOGL_DEBUG("TRACE4.ip.diffim.RegularizedKernelSolution.estimateRisk",
1119  "Lambda = %.3f, Risk = %.5e",
1120  l, risk);
1121  LOGL_DEBUG("TRACE5.ip.diffim.RegularizedKernelSolution.estimateRisk",
1122  "%.5e + 2 * (%.5e - %.5e)",
1123  term1(0), term2a, term2b(0));
1124  risks.push_back(risk);
1125  }
1126  std::vector<double>::iterator it = min_element(risks.begin(), risks.end());
1127  int index = distance(risks.begin(), it);
1128  LOGL_DEBUG("TRACE3.ip.diffim.RegularizedKernelSolution.estimateRisk",
1129  "Minimum Risk = %.3e at lambda = %.3e", risks[index], lambdas[index]);
1130 
1131  return lambdas[index];
1132 
1133  }
1134 
1135  template <typename InputT>
1136  Eigen::MatrixXd RegularizedKernelSolution<InputT>::getM(bool includeHmat) {
1137  if (includeHmat == true) {
1138  return this->_mMat + _lambda * _hMat;
1139  }
1140  else {
1141  return this->_mMat;
1142  }
1143  }
1144 
1145  template <typename InputT>
1147 
1148  LOGL_DEBUG("TRACE3.ip.diffim.RegularizedKernelSolution.solve",
1149  "cMat is %d x %d; vVec is %d; iVec is %d; hMat is %d x %d",
1150  this->_cMat.rows(), this->_cMat.cols(), this->_ivVec.size(),
1151  this->_iVec.size(), _hMat.rows(), _hMat.cols());
1152 
1153  if (DEBUG_MATRIX2) {
1154  std::cout << "ID: " << (this->_id) << std::endl;
1155  std::cout << "C:" << std::endl;
1156  std::cout << this->_cMat << std::endl;
1157  std::cout << "Sigma^{-1}:" << std::endl;
1158  std::cout << Eigen::MatrixXd(this->_ivVec.asDiagonal()) << std::endl;
1159  std::cout << "Y:" << std::endl;
1160  std::cout << this->_iVec << std::endl;
1161  std::cout << "H:" << std::endl;
1162  std::cout << _hMat << std::endl;
1163  }
1164 
1165 
1166  this->_mMat = this->_cMat.transpose() * this->_ivVec.asDiagonal() * this->_cMat;
1167  this->_bVec = this->_cMat.transpose() * this->_ivVec.asDiagonal() * this->_iVec;
1168 
1169 
1170  /* See N.R. 18.5
1171 
1172  Matrix equation to solve is Y = C a + N
1173  Y = vectorized version of I (I = image to not convolve)
1174  C_i = K_i (x) R (R = image to convolve)
1175  a = kernel coefficients
1176  N = noise
1177 
1178  If we reweight everything by the inverse square root of the noise
1179  covariance, we get a linear model with the identity matrix for
1180  the noise. The problem can then be solved using least squares,
1181  with the normal equations
1182 
1183  C^T Y = C^T C a
1184 
1185  or
1186 
1187  b = M a
1188 
1189  with
1190 
1191  b = C^T Y
1192  M = C^T C
1193  a = (C^T C)^{-1} C^T Y
1194 
1195 
1196  We can regularize the least square problem
1197 
1198  Y = C a + lambda a^T H a (+ N, which can be ignored)
1199 
1200  or the normal equations
1201 
1202  (C^T C + lambda H) a = C^T Y
1203 
1204 
1205  The solution to the regularization of the least squares problem is
1206 
1207  a = (C^T C + lambda H)^{-1} C^T Y
1208 
1209  The approximation to Y is
1210 
1211  C a = C (C^T C + lambda H)^{-1} C^T Y
1212 
1213  with smoothing matrix
1214 
1215  S = C (C^T C + lambda H)^{-1} C^T
1216 
1217  */
1218 
1219  std::string lambdaType = _policy.getString("lambdaType");
1220  if (lambdaType == "absolute") {
1221  _lambda = _policy.getDouble("lambdaValue");
1222  }
1223  else if (lambdaType == "relative") {
1224  _lambda = this->_mMat.trace() / this->_hMat.trace();
1225  _lambda *= _policy.getDouble("lambdaScaling");
1226  }
1227  else if (lambdaType == "minimizeBiasedRisk") {
1228  double tol = _policy.getDouble("maxConditionNumber");
1229  _lambda = estimateRisk(tol);
1230  }
1231  else if (lambdaType == "minimizeUnbiasedRisk") {
1233  }
1234  else {
1235  throw LSST_EXCEPT(pexExcept::Exception, "lambdaType in Policy not recognized");
1236  }
1237 
1238  LOGL_DEBUG("TRACE3.ip.diffim.RegularizedKernelSolution.solve",
1239  "Applying kernel regularization with lambda = %.2e", _lambda);
1240 
1241 
1242  try {
1243  KernelSolution::solve(this->_mMat + _lambda * _hMat, this->_bVec);
1244  } catch (pexExcept::Exception &e) {
1245  LSST_EXCEPT_ADD(e, "Unable to solve static kernel matrix");
1246  throw e;
1247  }
1248  /* Turn matrices into _kernel and _background */
1250  }
1251 
1252  template <typename InputT>
1254  std::vector<double> lambdas;
1255 
1256  std::string lambdaStepType = _policy.getString("lambdaStepType");
1257  if (lambdaStepType == "linear") {
1258  double lambdaLinMin = _policy.getDouble("lambdaLinMin");
1259  double lambdaLinMax = _policy.getDouble("lambdaLinMax");
1260  double lambdaLinStep = _policy.getDouble("lambdaLinStep");
1261  for (double l = lambdaLinMin; l <= lambdaLinMax; l += lambdaLinStep) {
1262  lambdas.push_back(l);
1263  }
1264  }
1265  else if (lambdaStepType == "log") {
1266  double lambdaLogMin = _policy.getDouble("lambdaLogMin");
1267  double lambdaLogMax = _policy.getDouble("lambdaLogMax");
1268  double lambdaLogStep = _policy.getDouble("lambdaLogStep");
1269  for (double l = lambdaLogMin; l <= lambdaLogMax; l += lambdaLogStep) {
1270  lambdas.push_back(pow(10, l));
1271  }
1272  }
1273  else {
1274  throw LSST_EXCEPT(pexExcept::Exception, "lambdaStepType in Policy not recognized");
1275  }
1276  return lambdas;
1277  }
1278 
1279  /*******************************************************************************************************/
1280 
1282  lsst::afw::math::KernelList const& basisList,
1283  lsst::afw::math::Kernel::SpatialFunctionPtr spatialKernelFunction,
1286  ) :
1287  KernelSolution(),
1288  _spatialKernelFunction(spatialKernelFunction),
1289  _constantFirstTerm(false),
1290  _kernel(),
1291  _background(background),
1292  _kSum(0.0),
1293  _policy(policy),
1294  _nbases(0),
1295  _nkt(0),
1296  _nbt(0),
1297  _nt(0) {
1298 
1299  bool isAlardLupton = _policy.getString("kernelBasisSet") == "alard-lupton";
1300  bool usePca = _policy.getBool("usePcaForSpatialKernel");
1301  if (isAlardLupton || usePca) {
1302  _constantFirstTerm = true;
1303  }
1304  this->_fitForBackground = _policy.getBool("fitForBackground");
1305 
1306  _nbases = basisList.size();
1307  _nkt = _spatialKernelFunction->getParameters().size();
1308  _nbt = _fitForBackground ? _background->getParameters().size() : 0;
1309  _nt = 0;
1310  if (_constantFirstTerm) {
1311  _nt = (_nbases - 1) * _nkt + 1 + _nbt;
1312  } else {
1313  _nt = _nbases * _nkt + _nbt;
1314  }
1315 
1316  Eigen::MatrixXd mMat(_nt, _nt);
1317  Eigen::VectorXd bVec(_nt);
1318  mMat.setZero();
1319  bVec.setZero();
1320 
1321  this->_mMat = mMat;
1322  this->_bVec = bVec;
1323 
1325  new afwMath::LinearCombinationKernel(basisList, *_spatialKernelFunction)
1326  );
1327 
1328  LOGL_DEBUG("TRACE3.ip.diffim.SpatialKernelSolution",
1329  "Initializing with size %d %d %d and constant first term = %s",
1330  _nkt, _nbt, _nt,
1331  _constantFirstTerm ? "true" : "false");
1332 
1333  }
1334 
1335  void SpatialKernelSolution::addConstraint(float xCenter, float yCenter,
1336  Eigen::MatrixXd const& qMat,
1337  Eigen::VectorXd const& wVec) {
1338 
1339  LOGL_DEBUG("TRACE5.ip.diffim.SpatialKernelSolution.addConstraint",
1340  "Adding candidate at %f, %f", xCenter, yCenter);
1341 
1342  /* Calculate P matrices */
1343  /* Pure kernel terms */
1344  Eigen::VectorXd pK(_nkt);
1345  std::vector<double> paramsK = _spatialKernelFunction->getParameters();
1346  for (int idx = 0; idx < _nkt; idx++) { paramsK[idx] = 0.0; }
1347  for (int idx = 0; idx < _nkt; idx++) {
1348  paramsK[idx] = 1.0;
1349  _spatialKernelFunction->setParameters(paramsK);
1350  pK(idx) = (*_spatialKernelFunction)(xCenter, yCenter); /* Assume things don't vary over stamp */
1351  paramsK[idx] = 0.0;
1352  }
1353  Eigen::MatrixXd pKpKt = (pK * pK.transpose());
1354 
1355  Eigen::VectorXd pB;
1356  Eigen::MatrixXd pBpBt;
1357  Eigen::MatrixXd pKpBt;
1358  if (_fitForBackground) {
1359  pB = Eigen::VectorXd(_nbt);
1360 
1361  /* Pure background terms */
1362  std::vector<double> paramsB = _background->getParameters();
1363  for (int idx = 0; idx < _nbt; idx++) { paramsB[idx] = 0.0; }
1364  for (int idx = 0; idx < _nbt; idx++) {
1365  paramsB[idx] = 1.0;
1366  _background->setParameters(paramsB);
1367  pB(idx) = (*_background)(xCenter, yCenter); /* Assume things don't vary over stamp */
1368  paramsB[idx] = 0.0;
1369  }
1370  pBpBt = (pB * pB.transpose());
1371 
1372  /* Cross terms */
1373  pKpBt = (pK * pB.transpose());
1374  }
1375 
1376  if (DEBUG_MATRIX) {
1377  std::cout << "Spatial weights" << std::endl;
1378  std::cout << "pKpKt " << pKpKt << std::endl;
1379  if (_fitForBackground) {
1380  std::cout << "pBpBt " << pBpBt << std::endl;
1381  std::cout << "pKpBt " << pKpBt << std::endl;
1382  }
1383  }
1384 
1385  if (DEBUG_MATRIX) {
1386  std::cout << "Spatial matrix inputs" << std::endl;
1387  std::cout << "M " << qMat << std::endl;
1388  std::cout << "B " << wVec << std::endl;
1389  }
1390 
1391  /* first index to start the spatial blocks; default=0 for non-constant first term */
1392  int m0 = 0;
1393  /* how many rows/cols to adjust the matrices/vectors; default=0 for non-constant first term */
1394  int dm = 0;
1395  /* where to start the background terms; this is always true */
1396  int mb = _nt - _nbt;
1397 
1398  if (_constantFirstTerm) {
1399  m0 = 1; /* we need to manually fill in the first (non-spatial) terms below */
1400  dm = _nkt-1; /* need to shift terms due to lack of spatial variation in first term */
1401 
1402  _mMat(0, 0) += qMat(0,0);
1403  for(int m2 = 1; m2 < _nbases; m2++) {
1404  _mMat.block(0, m2*_nkt-dm, 1, _nkt) += qMat(0,m2) * pK.transpose();
1405  }
1406  _bVec(0) += wVec(0);
1407 
1408  if (_fitForBackground) {
1409  _mMat.block(0, mb, 1, _nbt) += qMat(0,_nbases) * pB.transpose();
1410  }
1411  }
1412 
1413  /* Fill in the spatial blocks */
1414  for(int m1 = m0; m1 < _nbases; m1++) {
1415  /* Diagonal kernel-kernel term; only use upper triangular part of pKpKt */
1416  _mMat.block(m1*_nkt-dm, m1*_nkt-dm, _nkt, _nkt) +=
1417  (pKpKt * qMat(m1,m1)).triangularView<Eigen::Upper>();
1418 
1419  /* Kernel-kernel terms */
1420  for(int m2 = m1+1; m2 < _nbases; m2++) {
1421  _mMat.block(m1*_nkt-dm, m2*_nkt-dm, _nkt, _nkt) += qMat(m1,m2) * pKpKt;
1422  }
1423 
1424  if (_fitForBackground) {
1425  /* Kernel cross terms with background */
1426  _mMat.block(m1*_nkt-dm, mb, _nkt, _nbt) += qMat(m1,_nbases) * pKpBt;
1427  }
1428 
1429  /* B vector */
1430  _bVec.segment(m1*_nkt-dm, _nkt) += wVec(m1) * pK;
1431  }
1432 
1433  if (_fitForBackground) {
1434  /* Background-background terms only */
1435  _mMat.block(mb, mb, _nbt, _nbt) +=
1436  (pBpBt * qMat(_nbases,_nbases)).triangularView<Eigen::Upper>();
1437  _bVec.segment(mb, _nbt) += wVec(_nbases) * pB;
1438  }
1439 
1440  if (DEBUG_MATRIX) {
1441  std::cout << "Spatial matrix outputs" << std::endl;
1442  std::cout << "mMat " << _mMat << std::endl;
1443  std::cout << "bVec " << _bVec << std::endl;
1444  }
1445 
1446  }
1447 
1450  throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return image");
1451  }
1453  new afwImage::Image<afwMath::Kernel::Pixel>(_kernel->getDimensions())
1454  );
1455  (void)_kernel->computeImage(*image, false, pos[0], pos[1]);
1456  return image;
1457  }
1458 
1460  /* Fill in the other half of mMat */
1461  for (int i = 0; i < _nt; i++) {
1462  for (int j = i+1; j < _nt; j++) {
1463  _mMat(j,i) = _mMat(i,j);
1464  }
1465  }
1466 
1467  try {
1469  } catch (pexExcept::Exception &e) {
1470  LSST_EXCEPT_ADD(e, "Unable to solve spatial kernel matrix");
1471  throw e;
1472  }
1473  /* Turn matrices into _kernel and _background */
1474  _setKernel();
1475  }
1476 
1480  throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return solution");
1481  }
1482 
1483  return std::make_pair(_kernel, _background);
1484  }
1485 
1486  void SpatialKernelSolution::_setKernel() {
1487  /* Report on condition number */
1488  double cNumber = this->getConditionNumber(EIGENVALUE);
1489 
1490  if (_nkt == 1) {
1491  /* Not spatially varying; this fork is a specialization for convolution speed--up */
1492 
1493  /* Set the basis coefficients */
1494  std::vector<double> kCoeffs(_nbases);
1495  for (int i = 0; i < _nbases; i++) {
1496  if (std::isnan(_aVec(i))) {
1497  throw LSST_EXCEPT(
1500  "I. Unable to determine spatial kernel solution %d (nan). Condition number = %.3e") % i % cNumber));
1501  }
1502  kCoeffs[i] = _aVec(i);
1503  }
1504  lsst::afw::math::KernelList basisList =
1506  _kernel.reset(
1507  new afwMath::LinearCombinationKernel(basisList, kCoeffs)
1508  );
1509  }
1510  else {
1511 
1512  /* Set the kernel coefficients */
1514  kCoeffs.reserve(_nbases);
1515  for (int i = 0, idx = 0; i < _nbases; i++) {
1516  kCoeffs.push_back(std::vector<double>(_nkt));
1517 
1518  /* Deal with the possibility the first term doesn't vary spatially */
1519  if ((i == 0) && (_constantFirstTerm)) {
1520  if (std::isnan(_aVec(idx))) {
1521  throw LSST_EXCEPT(
1524  "II. Unable to determine spatial kernel solution %d (nan). Condition number = %.3e") % idx % cNumber));
1525  }
1526  kCoeffs[i][0] = _aVec(idx++);
1527  }
1528  else {
1529  for (int j = 0; j < _nkt; j++) {
1530  if (std::isnan(_aVec(idx))) {
1531  throw LSST_EXCEPT(
1534  "III. Unable to determine spatial kernel solution %d (nan). Condition number = %.3e") % idx % cNumber));
1535  }
1536  kCoeffs[i][j] = _aVec(idx++);
1537  }
1538  }
1539  }
1540  _kernel->setSpatialParameters(kCoeffs);
1541  }
1542 
1543  /* Set kernel Sum */
1544  std::shared_ptr<ImageT> image (new ImageT(_kernel->getDimensions()));
1545  _kSum = _kernel->computeImage(*image, false);
1546 
1547  /* Set the background coefficients */
1548  std::vector<double> bgCoeffs(_fitForBackground ? _nbt : 1);
1549  if (_fitForBackground) {
1550  for (int i = 0; i < _nbt; i++) {
1551  int idx = _nt - _nbt + i;
1552  if (std::isnan(_aVec(idx))) {
1555  "Unable to determine spatial background solution %d (nan)") % (idx)));
1556  }
1557  bgCoeffs[i] = _aVec(idx);
1558  }
1559  }
1560  else {
1561  bgCoeffs[0] = 0.;
1562  }
1563  _background->setParameters(bgCoeffs);
1564  }
1565 
1566 /***********************************************************************************************************/
1567 //
1568 // Explicit instantiations
1569 //
1570  typedef float InputT;
1571 
1572  template class StaticKernelSolution<InputT>;
1573  template class MaskedKernelSolution<InputT>;
1574  template class RegularizedKernelSolution<InputT>;
1575 
1576 }}} // end of namespace lsst::ip::diffim
uint64_t * ptr
Definition: RangeSet.cc:88
int getHeight() const noexcept
Definition: Box.h:175
Eigen::MatrixXi maskToEigenMatrix(lsst::afw::image::Mask< lsst::afw::image::MaskPixel > const &mask)
estimate sample minimum
Definition: Statistics.h:76
a container for holding hierarchical configuration data in memory.
Definition: Policy.h:169
virtual std::pair< std::shared_ptr< lsst::afw::math::Kernel >, double > getSolutionPair()
Eigen::VectorXd _bVec
Derived least squares B vector.
T endl(T... args)
SpatialKernelSolution(lsst::afw::math::KernelList const &basisList, lsst::afw::math::Kernel::SpatialFunctionPtr spatialKernelFunction, lsst::afw::math::Kernel::SpatialFunctionPtr background, lsst::pex::policy::Policy policy)
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:117
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.
const std::string getString(const std::string &name) const
return a string value associated with the given name .
Definition: Policy.h:631
#define LOGL_DEBUG(logger, message...)
Log a debug-level message using a varargs/printf style interface.
Definition: Log.h:489
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:280
STL class.
LSST DM logging module built on log4cxx.
RegularizedKernelSolution(lsst::afw::math::KernelList const &basisList, bool fitForBackground, Eigen::MatrixXd const &hMat, lsst::pex::policy::Policy policy)
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:463
Represent a 2-dimensional array of bitmask pixels.
Definition: Mask.h:78
lsst::afw::image::Image< lsst::afw::math::Kernel::Pixel > ImageT
int getMaxY() const noexcept
Definition: Box.h:149
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:741
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:168
T make_pair(T... args)
int getWidth() const noexcept
Definition: Box.h:174
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:383
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...
std::shared_ptr< lsst::afw::image::Image< lsst::afw::math::Kernel::Pixel > > makeKernelImage(lsst::afw::geom::Point2D const &pos)
int getMaxX() const noexcept
Definition: Box.h:148
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
bool getBool(const std::string &name) const
return a boolean value associated with the given name.
Definition: Policy.h:589
Eigen::MatrixXd _mMat
Derived least squares M matrix.
#define DEBUG_MATRIX
double getDouble(const std::string &name) const
return a double value associated with the given name.
Definition: Policy.h:617
int getMinX() const noexcept
Definition: Box.h:144
lsst::geom::Point2I getXY0() const
Return the image&#39;s origin.
Definition: ImageBase.h:341
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
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)
afw::table::Key< afw::table::Array< ImagePixelT > > image
Reports invalid arguments.
Definition: Runtime.h:66
T isnan(T... args)
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::afw::geom::Box2I maskBox)
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:275
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:374
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects...
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:54
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:59
Image Subtraction helper functions.
int end
int getMinY() const noexcept
Definition: Box.h:145
T reserve(T... args)