LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
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/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
43namespace geom = lsst::geom;
45namespace afwMath = lsst::afw::math;
46namespace afwGeom = lsst::afw::geom;
47namespace afwImage = lsst::afw::image;
49
50namespace lsst {
51namespace ip {
52namespace 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
93
95 return getConditionNumber(_mMat, conditionType);
96 }
97
98 double KernelSolution::getConditionNumber(Eigen::MatrixXd const& mMat,
99 ConditionNumberType conditionType) {
100 switch (conditionType) {
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:
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::cpu_timer t;
144
145 LOGL_DEBUG("TRACE2.ip.diffim.KernelSolution.solve",
146 "Solving for kernel");
147 _solvedBy = LU;
148 Eigen::FullPivLU<Eigen::MatrixXd> lu(mMat);
149 if (lu.isInvertible()) {
150 aVec = lu.solve(bVec);
151 } else {
152 LOGL_DEBUG("TRACE3.ip.diffim.KernelSolution.solve",
153 "Unable to determine kernel via LU");
154 /* LAST RESORT */
155 try {
156
158 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eVecValues(mMat);
159 Eigen::MatrixXd const& rMat = eVecValues.eigenvectors();
160 Eigen::VectorXd eValues = eVecValues.eigenvalues();
161
162 for (int i = 0; i != eValues.rows(); ++i) {
163 if (eValues(i) != 0.0) {
164 eValues(i) = 1.0/eValues(i);
165 }
166 }
167
168 aVec = rMat * eValues.asDiagonal() * rMat.transpose() * bVec;
169 } catch (pexExcept::Exception& e) {
170
171 _solvedBy = NONE;
172 LOGL_DEBUG("TRACE3.ip.diffim.KernelSolution.solve",
173 "Unable to determine kernel via eigen-values");
174
175 throw LSST_EXCEPT(pexExcept::Exception, "Unable to determine kernel solution");
176 }
177 }
178
179 t.stop();
180 double time = 1e-9 * t.elapsed().wall;
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>
216 if (_solvedBy == KernelSolution::NONE) {
217 throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return solution");
218 }
219 return _kernel;
220 }
221
222 template <typename InputT>
224 if (_solvedBy == KernelSolution::NONE) {
225 throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return image");
226 }
228 new afwImage::Image<afwMath::Kernel::Pixel>(_kernel->getDimensions())
229 );
230 (void)_kernel->computeImage(*image, false);
231 return image;
232 }
233
234 template <typename InputT>
236 if (_solvedBy == KernelSolution::NONE) {
237 throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return background");
238 }
239 return _background;
240 }
241
242 template <typename InputT>
244 if (_solvedBy == KernelSolution::NONE) {
245 throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return ksum");
246 }
247 return _kSum;
248 }
249
250 template <typename InputT>
253 if (_solvedBy == KernelSolution::NONE) {
254 throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return solution");
255 }
256
257 return std::make_pair(_kernel, _background);
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
278 std::dynamic_pointer_cast<afwMath::LinearCombinationKernel>(_kernel)->getKernelList();
279
280 unsigned int const nKernelParameters = basisList.size();
281 unsigned int const nBackgroundParameters = _fitForBackground ? 1 : 0;
282 unsigned int const nParameters = nKernelParameters + nBackgroundParameters;
283
284 std::vector<std::shared_ptr<afwMath::Kernel> >::const_iterator kiter = basisList.begin();
285
286 /* Ignore buffers around edge of convolved images :
287 *
288 * If the kernel has width 5, it has center pixel 2. The first good pixel
289 * is the (5-2)=3rd pixel, which is array index 2, and ends up being the
290 * index of the central pixel.
291 *
292 * You also have a buffer of unusable pixels on the other side, numbered
293 * width-center-1. The last good usable pixel is N-width+center+1.
294 *
295 * Example : the kernel is width = 5, center = 2
296 *
297 * ---|---|-c-|---|---|
298 *
299 * the image is width = N
300 * convolve this with the kernel, and you get
301 *
302 * |-x-|-x-|-g-|---|---| ... |---|---|-g-|-x-|-x-|
303 *
304 * g = first/last good pixel
305 * x = bad
306 *
307 * the first good pixel is the array index that has the value "center", 2
308 * the last good pixel has array index N-(5-2)+1
309 * eg. if N = 100, you want to use up to index 97
310 * 100-3+1 = 98, and the loops use i < 98, meaning the last
311 * index you address is 97.
312 */
313
314 /* NOTE - we are accessing particular elements of Eigen arrays using
315 these coordinates, therefore they need to be in LOCAL coordinates.
316 This was written before ndarray unification.
317 */
318 geom::Box2I goodBBox = (*kiter)->shrinkBBox(templateImage.getBBox(afwImage::LOCAL));
319 unsigned int const startCol = goodBBox.getMinX();
320 unsigned int const startRow = goodBBox.getMinY();
321 // endCol/Row is one past the index of the last good col/row
322 unsigned int endCol = goodBBox.getMaxX() + 1;
323 unsigned int endRow = goodBBox.getMaxY() + 1;
324
325 boost::timer::cpu_timer t;
326
327 /* Eigen representation of input images; only the pixels that are unconvolved in cimage below */
328 Eigen::MatrixXd eigenTemplate = imageToEigenMatrix(templateImage).block(startRow,
329 startCol,
330 endRow-startRow,
331 endCol-startCol);
332 Eigen::MatrixXd eigenScience = imageToEigenMatrix(scienceImage).block(startRow,
333 startCol,
334 endRow-startRow,
335 endCol-startCol);
336 Eigen::MatrixXd eigeniVariance = imageToEigenMatrix(varianceEstimate).block(
337 startRow, startCol, endRow-startRow, endCol-startCol
338 ).array().inverse().matrix();
339
340 /* Resize into 1-D for later usage */
341 eigenTemplate.resize(eigenTemplate.rows()*eigenTemplate.cols(), 1);
342 eigenScience.resize(eigenScience.rows()*eigenScience.cols(), 1);
343 eigeniVariance.resize(eigeniVariance.rows()*eigeniVariance.cols(), 1);
344
345 /* Holds image convolved with basis function */
346 afwImage::Image<PixelT> cimage(templateImage.getDimensions());
347
348 /* Holds eigen representation of image convolved with all basis functions */
349 std::vector<Eigen::MatrixXd> convolvedEigenList(nKernelParameters);
350
351 /* Iterators over convolved image list and basis list */
352 typename std::vector<Eigen::MatrixXd>::iterator eiter = convolvedEigenList.begin();
353 /* Create C_i in the formalism of Alard & Lupton */
354 afwMath::ConvolutionControl convolutionControl;
355 convolutionControl.setDoNormalize(false);
356 for (kiter = basisList.begin(); kiter != basisList.end(); ++kiter, ++eiter) {
357 afwMath::convolve(cimage, templateImage, **kiter, convolutionControl); /* cimage stores convolved image */
358
359 Eigen::MatrixXd cMat = imageToEigenMatrix(cimage).block(startRow,
360 startCol,
361 endRow-startRow,
362 endCol-startCol);
363 cMat.resize(cMat.size(), 1);
364 *eiter = cMat;
365
366 }
367
368 t.stop();
369 double time = 1e-9 * t.elapsed().wall;
370 LOGL_DEBUG("TRACE3.ip.diffim.StaticKernelSolution.build",
371 "Total compute time to do basis convolutions : %.2f s", time);
372
373 /*
374 Load matrix with all values from convolvedEigenList : all images
375 (eigeniVariance, convolvedEigenList) must be the same size
376 */
377 Eigen::MatrixXd cMat(eigenTemplate.col(0).size(), nParameters);
378 typename std::vector<Eigen::MatrixXd>::iterator eiterj = convolvedEigenList.begin();
379 typename std::vector<Eigen::MatrixXd>::iterator eiterE = convolvedEigenList.end();
380 for (unsigned int kidxj = 0; eiterj != eiterE; eiterj++, kidxj++) {
381 cMat.col(kidxj) = eiterj->col(0);
382 }
383 /* Treat the last "image" as all 1's to do the background calculation. */
384 if (_fitForBackground)
385 cMat.col(nParameters-1).fill(1.);
386
387 _cMat = cMat;
388 _ivVec = eigeniVariance.col(0);
389 _iVec = eigenScience.col(0);
390
391 /* Make these outside of solve() so I can check condition number */
392 _mMat = _cMat.transpose() * (_ivVec.asDiagonal() * _cMat);
393 _bVec = _cMat.transpose() * (_ivVec.asDiagonal() * _iVec);
394 }
395
396 template <typename InputT>
398 LOGL_DEBUG("TRACE3.ip.diffim.StaticKernelSolution.solve",
399 "mMat is %d x %d; bVec is %d; cMat is %d x %d; vVec is %d; iVec is %d",
400 _mMat.rows(), _mMat.cols(), _bVec.size(),
401 _cMat.rows(), _cMat.cols(), _ivVec.size(), _iVec.size());
402
403 if (DEBUG_MATRIX) {
404 std::cout << "C" << std::endl;
405 std::cout << _cMat << std::endl;
406 std::cout << "iV" << std::endl;
407 std::cout << _ivVec << std::endl;
408 std::cout << "I" << std::endl;
409 std::cout << _iVec << std::endl;
410 }
411
412 try {
414 } catch (pexExcept::Exception &e) {
415 LSST_EXCEPT_ADD(e, "Unable to solve static kernel matrix");
416 throw e;
417 }
418 /* Turn matrices into _kernel and _background */
419 _setKernel();
420 }
421
422 template <typename InputT>
424 if (_solvedBy == KernelSolution::NONE) {
425 throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot make solution");
426 }
427
428 unsigned int const nParameters = _aVec.size();
429 unsigned int const nBackgroundParameters = _fitForBackground ? 1 : 0;
430 unsigned int const nKernelParameters =
431 std::dynamic_pointer_cast<afwMath::LinearCombinationKernel>(_kernel)->getKernelList().size();
432 if (nParameters != (nKernelParameters + nBackgroundParameters))
433 throw LSST_EXCEPT(pexExcept::Exception, "Mismatched sizes in kernel solution");
434
435 /* Fill in the kernel results */
436 std::vector<double> kValues(nKernelParameters);
437 for (unsigned int idx = 0; idx < nKernelParameters; idx++) {
438 if (std::isnan(_aVec(idx))) {
440 str(boost::format("Unable to determine kernel solution %d (nan)") % idx));
441 }
442 kValues[idx] = _aVec(idx);
443 }
444 _kernel->setKernelParameters(kValues);
445
447 new ImageT(_kernel->getDimensions())
448 );
449 _kSum = _kernel->computeImage(*image, false);
450
451 if (_fitForBackground) {
452 if (std::isnan(_aVec(nParameters-1))) {
454 str(boost::format("Unable to determine background solution %d (nan)") %
455 (nParameters-1)));
456 }
457 _background = _aVec(nParameters-1);
458 }
459 }
460
461
462 template <typename InputT>
464 throw LSST_EXCEPT(pexExcept::Exception, "Uncertainty calculation not supported");
465
466 /* Estimate of parameter uncertainties comes from the inverse of the
467 * covariance matrix (noise spectrum).
468 * N.R. 15.4.8 to 15.4.15
469 *
470 * Since this is a linear problem no need to use Fisher matrix
471 * N.R. 15.5.8
472 *
473 * Although I might be able to take advantage of the solution above.
474 * Since this now works and is not the rate limiting step, keep as-is for DC3a.
475 *
476 * Use Cholesky decomposition again.
477 * Cholkesy:
478 * Cov = L L^t
479 * Cov^(-1) = (L L^t)^(-1)
480 * = (L^T)^-1 L^(-1)
481 *
482 * Code would be:
483 *
484 * Eigen::MatrixXd cov = _mMat.transpose() * _mMat;
485 * Eigen::LLT<Eigen::MatrixXd> llt = cov.llt();
486 * Eigen::MatrixXd error2 = llt.matrixL().transpose().inverse()*llt.matrixL().inverse();
487 */
488 }
489
490 /*******************************************************************************************************/
491
492 template <typename InputT>
494 lsst::afw::math::KernelList const& basisList,
495 bool fitForBackground
496 ) :
497 StaticKernelSolution<InputT>(basisList, fitForBackground)
498 {};
499
500 template <typename InputT>
502 lsst::afw::image::Image<InputT> const &templateImage,
503 lsst::afw::image::Image<InputT> const &scienceImage,
506 ) {
507
508 afwMath::Statistics varStats = afwMath::makeStatistics(varianceEstimate, afwMath::MIN);
509 if (varStats.getValue(afwMath::MIN) < 0.0) {
511 "Error: variance less than 0.0");
512 }
513 if (varStats.getValue(afwMath::MIN) == 0.0) {
515 "Error: variance equals 0.0, cannot inverse variance weight");
516 }
517
518 /* Full footprint of all input images */
520 new afwDet::Footprint(std::make_shared<afwGeom::SpanSet>(templateImage.getBBox())));
521
522 afwMath::KernelList basisList =
523 std::dynamic_pointer_cast<afwMath::LinearCombinationKernel>(this->_kernel)->getKernelList();
524 std::vector<std::shared_ptr<afwMath::Kernel> >::const_iterator kiter = basisList.begin();
525
526 /* Only BAD pixels marked in this mask */
527 afwImage::MaskPixel bitMask =
532
533 /* Create a Footprint that contains all the masked pixels set above */
534 afwDet::Threshold threshold = afwDet::Threshold(bitMask, afwDet::Threshold::BITMASK, true);
535 afwDet::FootprintSet maskFpSet(pixelMask, threshold, true);
536
537 /* And spread it by the kernel half width */
538 int growPix = (*kiter)->getCtr().getX();
539 afwDet::FootprintSet maskedFpSetGrown(maskFpSet, growPix, true);
540
541#if 0
542 for (typename afwDet::FootprintSet::FootprintList::iterator
543 ptr = maskedFpSetGrown.getFootprints()->begin(),
544 end = maskedFpSetGrown.getFootprints()->end();
545 ptr != end;
546 ++ptr) {
547
548 afwDet::setMaskFromFootprint(finalMask,
549 (**ptr),
551 }
552#endif
553
555 for (auto const & foot : *(maskedFpSetGrown.getFootprints())) {
556 foot->getSpans()->setMask(finalMask, afwImage::Mask<afwImage::MaskPixel>::getPlaneBitMask("BAD"));
557 }
558 pixelMask.writeFits("pixelmask.fits");
559 finalMask.writeFits("finalmask.fits");
560
561
562 ndarray::Array<int, 1, 1> maskArray =
563 ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
564 fullFp->getSpans()->flatten(maskArray, finalMask.getArray(), templateImage.getXY0());
565 auto maskEigen = ndarray::asEigenMatrix(maskArray);
566
567 ndarray::Array<InputT, 1, 1> arrayTemplate =
568 ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
569 fullFp->getSpans()->flatten(arrayTemplate, templateImage.getArray(), templateImage.getXY0());
570 auto eigenTemplate0 = ndarray::asEigenMatrix(arrayTemplate);
571
572 ndarray::Array<InputT, 1, 1> arrayScience =
573 ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
574 fullFp->getSpans()->flatten(arrayScience, scienceImage.getArray(), scienceImage.getXY0());
575 auto eigenScience0 = ndarray::asEigenMatrix(arrayScience);
576
577 ndarray::Array<afwImage::VariancePixel, 1, 1> arrayVariance =
578 ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
579 fullFp->getSpans()->flatten(arrayVariance, varianceEstimate.getArray(), varianceEstimate.getXY0());
580 auto eigenVariance0 = ndarray::asEigenMatrix(arrayVariance);
581
582 int nGood = 0;
583 for (int i = 0; i < maskEigen.size(); i++) {
584 if (maskEigen(i) == 0.0)
585 nGood += 1;
586 }
587
588 Eigen::VectorXd eigenTemplate(nGood);
589 Eigen::VectorXd eigenScience(nGood);
590 Eigen::VectorXd eigenVariance(nGood);
591 int nUsed = 0;
592 for (int i = 0; i < maskEigen.size(); i++) {
593 if (maskEigen(i) == 0.0) {
594 eigenTemplate(nUsed) = eigenTemplate0(i);
595 eigenScience(nUsed) = eigenScience0(i);
596 eigenVariance(nUsed) = eigenVariance0(i);
597 nUsed += 1;
598 }
599 }
600
601
602 boost::timer::cpu_timer t;
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 afwMath::ConvolutionControl convolutionControl;
619 convolutionControl.setDoNormalize(false);
620 for (kiter = basisList.begin(); kiter != basisList.end(); ++kiter, ++eiter) {
621 afwMath::convolve(cimage, templateImage, **kiter, convolutionControl); /* cimage stores convolved image */
622
623 ndarray::Array<InputT, 1, 1> arrayC =
624 ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
625 fullFp->getSpans()->flatten(arrayC, cimage.getArray(), cimage.getXY0());
626 auto eigenC0 = ndarray::asEigenMatrix(arrayC);
627
628 Eigen::VectorXd eigenC(nGood);
629 int nUsed = 0;
630 for (int i = 0; i < maskEigen.size(); i++) {
631 if (maskEigen(i) == 0.0) {
632 eigenC(nUsed) = eigenC0(i);
633 nUsed += 1;
634 }
635 }
636
637 *eiter = eigenC;
638 }
639 t.stop();
640 double time = 1e-9 * t.elapsed().wall;
641 LOGL_DEBUG("TRACE3.ip.diffim.StaticKernelSolution.buildWithMask",
642 "Total compute time to do basis convolutions : %.2f s", time);
643
644 /* Load matrix with all convolved images */
645 Eigen::MatrixXd cMat(eigenTemplate.size(), nParameters);
646 typename std::vector<Eigen::VectorXd>::iterator eiterj = convolvedEigenList.begin();
647 typename std::vector<Eigen::VectorXd>::iterator eiterE = convolvedEigenList.end();
648 for (unsigned int kidxj = 0; eiterj != eiterE; eiterj++, kidxj++) {
649 cMat.block(0, kidxj, eigenTemplate.size(), 1) =
650 Eigen::MatrixXd(*eiterj).block(0, 0, eigenTemplate.size(), 1);
651 }
652 /* Treat the last "image" as all 1's to do the background calculation. */
653 if (this->_fitForBackground)
654 cMat.col(nParameters-1).fill(1.);
655
656 this->_cMat = cMat;
657 this->_ivVec = eigenVariance.array().inverse().matrix();
658 this->_iVec = eigenScience;
659
660 /* Make these outside of solve() so I can check condition number */
661 this->_mMat = this->_cMat.transpose() * this->_ivVec.asDiagonal() * (this->_cMat);
662 this->_bVec = this->_cMat.transpose() * this->_ivVec.asDiagonal() * (this->_iVec);
663 }
664
665
666 template <typename InputT>
668 lsst::afw::image::Image<InputT> const &templateImage,
669 lsst::afw::image::Image<InputT> const &scienceImage,
672 ) {
673
674 afwMath::Statistics varStats = afwMath::makeStatistics(varianceEstimate, afwMath::MIN);
675 if (varStats.getValue(afwMath::MIN) < 0.0) {
677 "Error: variance less than 0.0");
678 }
679 if (varStats.getValue(afwMath::MIN) == 0.0) {
681 "Error: variance equals 0.0, cannot inverse variance weight");
682 }
683
685 std::dynamic_pointer_cast<afwMath::LinearCombinationKernel>(this->_kernel)->getKernelList();
686 std::vector<std::shared_ptr<afwMath::Kernel> >::const_iterator kiter = basisList.begin();
687
688 /* Only BAD pixels marked in this mask */
690 afwImage::MaskPixel bitMask =
694 sMask &= bitMask;
695 /* TBD: Need to figure out a way to spread this; currently done in Python */
696
697 unsigned int const nKernelParameters = basisList.size();
698 unsigned int const nBackgroundParameters = this->_fitForBackground ? 1 : 0;
699 unsigned int const nParameters = nKernelParameters + nBackgroundParameters;
700
701 /* NOTE - we are accessing particular elements of Eigen arrays using
702 these coordinates, therefore they need to be in LOCAL coordinates.
703 This was written before ndarray unification.
704 */
705 /* Ignore known EDGE pixels for speed */
706 geom::Box2I shrunkLocalBBox = (*kiter)->shrinkBBox(templateImage.getBBox(afwImage::LOCAL));
707 LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
708 "Limits of good pixels after convolution: %d,%d -> %d,%d (local)",
709 shrunkLocalBBox.getMinX(), shrunkLocalBBox.getMinY(),
710 shrunkLocalBBox.getMaxX(), shrunkLocalBBox.getMaxY());
711
712 /* Subimages are addressed in raw pixel coordinates */
713 unsigned int startCol = shrunkLocalBBox.getMinX();
714 unsigned int startRow = shrunkLocalBBox.getMinY();
715 unsigned int endCol = shrunkLocalBBox.getMaxX();
716 unsigned int endRow = shrunkLocalBBox.getMaxY();
717 /* Needed for Eigen block slicing operation */
718 endCol += 1;
719 endRow += 1;
720
721 boost::timer::cpu_timer t;
722
723 /* Eigen representation of input images; only the pixels that are unconvolved in cimage below */
724 Eigen::MatrixXi eMask = maskToEigenMatrix(sMask).block(startRow,
725 startCol,
726 endRow-startRow,
727 endCol-startCol);
728
729 Eigen::MatrixXd eigenTemplate = imageToEigenMatrix(templateImage).block(startRow,
730 startCol,
731 endRow-startRow,
732 endCol-startCol);
733 Eigen::MatrixXd eigenScience = imageToEigenMatrix(scienceImage).block(startRow,
734 startCol,
735 endRow-startRow,
736 endCol-startCol);
737 Eigen::MatrixXd eigeniVariance = imageToEigenMatrix(varianceEstimate).block(
738 startRow, startCol, endRow-startRow, endCol-startCol
739 ).array().inverse().matrix();
740
741 /* Resize into 1-D for later usage */
742 eMask.resize(eMask.rows()*eMask.cols(), 1);
743 eigenTemplate.resize(eigenTemplate.rows()*eigenTemplate.cols(), 1);
744 eigenScience.resize(eigenScience.rows()*eigenScience.cols(), 1);
745 eigeniVariance.resize(eigeniVariance.rows()*eigeniVariance.cols(), 1);
746
747 /* Do the masking, slowly for now... */
748 Eigen::MatrixXd maskedEigenTemplate(eigenTemplate.rows(), 1);
749 Eigen::MatrixXd maskedEigenScience(eigenScience.rows(), 1);
750 Eigen::MatrixXd maskedEigeniVariance(eigeniVariance.rows(), 1);
751 int nGood = 0;
752 for (int i = 0; i < eMask.rows(); i++) {
753 if (eMask(i, 0) == 0) {
754 maskedEigenTemplate(nGood, 0) = eigenTemplate(i, 0);
755 maskedEigenScience(nGood, 0) = eigenScience(i, 0);
756 maskedEigeniVariance(nGood, 0) = eigeniVariance(i, 0);
757 nGood += 1;
758 }
759 }
760 /* Can't resize() since all values are lost; use blocks */
761 eigenTemplate = maskedEigenTemplate.block(0, 0, nGood, 1);
762 eigenScience = maskedEigenScience.block(0, 0, nGood, 1);
763 eigeniVariance = maskedEigeniVariance.block(0, 0, nGood, 1);
764
765
766 /* Holds image convolved with basis function */
767 afwImage::Image<InputT> cimage(templateImage.getDimensions());
768
769 /* Holds eigen representation of image convolved with all basis functions */
770 std::vector<Eigen::MatrixXd> convolvedEigenList(nKernelParameters);
771
772 /* Iterators over convolved image list and basis list */
773 typename std::vector<Eigen::MatrixXd>::iterator eiter = convolvedEigenList.begin();
774 /* Create C_i in the formalism of Alard & Lupton */
775 afwMath::ConvolutionControl convolutionControl;
776 convolutionControl.setDoNormalize(false);
777 for (kiter = basisList.begin(); kiter != basisList.end(); ++kiter, ++eiter) {
778 afwMath::convolve(cimage, templateImage, **kiter, convolutionControl); /* cimage stores convolved image */
779
780 Eigen::MatrixXd cMat = imageToEigenMatrix(cimage).block(startRow,
781 startCol,
782 endRow-startRow,
783 endCol-startCol);
784 cMat.resize(cMat.size(), 1);
785
786 /* Do masking */
787 Eigen::MatrixXd maskedcMat(cMat.rows(), 1);
788 int nGood = 0;
789 for (int i = 0; i < eMask.rows(); i++) {
790 if (eMask(i, 0) == 0) {
791 maskedcMat(nGood, 0) = cMat(i, 0);
792 nGood += 1;
793 }
794 }
795 cMat = maskedcMat.block(0, 0, nGood, 1);
796 *eiter = cMat;
797 }
798
799 t.stop();
800 double time = 1e-9 * t.elapsed().wall;
801 LOGL_DEBUG("TRACE3.ip.diffim.StaticKernelSolution.build",
802 "Total compute time to do basis convolutions : %.2f s", time);
803
804 /*
805 Load matrix with all values from convolvedEigenList : all images
806 (eigeniVariance, convolvedEigenList) must be the same size
807 */
808 Eigen::MatrixXd cMat(eigenTemplate.col(0).size(), nParameters);
809 typename std::vector<Eigen::MatrixXd>::iterator eiterj = convolvedEigenList.begin();
810 typename std::vector<Eigen::MatrixXd>::iterator eiterE = convolvedEigenList.end();
811 for (unsigned int kidxj = 0; eiterj != eiterE; eiterj++, kidxj++) {
812 cMat.col(kidxj) = eiterj->col(0);
813 }
814 /* Treat the last "image" as all 1's to do the background calculation. */
815 if (this->_fitForBackground)
816 cMat.col(nParameters-1).fill(1.);
817
818 this->_cMat = cMat;
819 this->_ivVec = eigeniVariance.col(0);
820 this->_iVec = eigenScience.col(0);
821
822 /* Make these outside of solve() so I can check condition number */
823 this->_mMat = this->_cMat.transpose() * this->_ivVec.asDiagonal() * this->_cMat;
824 this->_bVec = this->_cMat.transpose() * this->_ivVec.asDiagonal() * this->_iVec;
825
826 }
827
828 /* NOTE - this was written before the ndarray unification. I am rewriting
829 buildSingleMask to use the ndarray stuff */
830 template <typename InputT>
832 lsst::afw::image::Image<InputT> const &templateImage,
833 lsst::afw::image::Image<InputT> const &scienceImage,
835 lsst::geom::Box2I maskBox
836 ) {
837
838 afwMath::Statistics varStats = afwMath::makeStatistics(varianceEstimate, afwMath::MIN);
839 if (varStats.getValue(afwMath::MIN) < 0.0) {
841 "Error: variance less than 0.0");
842 }
843 if (varStats.getValue(afwMath::MIN) == 0.0) {
845 "Error: variance equals 0.0, cannot inverse variance weight");
846 }
847
849 std::dynamic_pointer_cast<afwMath::LinearCombinationKernel>(this->_kernel)->getKernelList();
850
851 unsigned int const nKernelParameters = basisList.size();
852 unsigned int const nBackgroundParameters = this->_fitForBackground ? 1 : 0;
853 unsigned int const nParameters = nKernelParameters + nBackgroundParameters;
854
855 std::vector<std::shared_ptr<afwMath::Kernel> >::const_iterator kiter = basisList.begin();
856
857 /*
858 NOTE : If we are using these views in Afw's Image space, we need to
859 make sure and compensate for the XY0 of the image:
860
861 geom::Box2I fullBBox = templateImage.getBBox();
862 int maskStartCol = maskBox.getMinX();
863 int maskStartRow = maskBox.getMinY();
864 int maskEndCol = maskBox.getMaxX();
865 int maskEndRow = maskBox.getMaxY();
866
867
868 If we are going to be doing the slicing in Eigen matrices derived from
869 the images, ignore the XY0.
870
871 geom::Box2I fullBBox = templateImage.getBBox(afwImage::LOCAL);
872
873 int maskStartCol = maskBox.getMinX() - templateImage.getX0();
874 int maskStartRow = maskBox.getMinY() - templateImage.getY0();
875 int maskEndCol = maskBox.getMaxX() - templateImage.getX0();
876 int maskEndRow = maskBox.getMaxY() - templateImage.getY0();
877
878 */
879
880
881 geom::Box2I shrunkBBox = (*kiter)->shrinkBBox(templateImage.getBBox());
882
883 LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
884 "Limits of good pixels after convolution: %d,%d -> %d,%d",
885 shrunkBBox.getMinX(), shrunkBBox.getMinY(),
886 shrunkBBox.getMaxX(), shrunkBBox.getMaxY());
887
888 unsigned int const startCol = shrunkBBox.getMinX();
889 unsigned int const startRow = shrunkBBox.getMinY();
890 unsigned int const endCol = shrunkBBox.getMaxX();
891 unsigned int const endRow = shrunkBBox.getMaxY();
892
893 /* NOTE: no endCol/endRow += 1 for block slicing, since we are doing the
894 slicing using Afw, not Eigen
895
896 Eigen arrays have index 0,0 in the upper right, while LSST images
897 have 0,0 in the lower left. The y-axis is flipped. When translating
898 Images to Eigen matrices in ipDiffim::imageToEigenMatrix this is
899 accounted for. However, we still need to be aware of this fact if
900 addressing subregions of an Eigen matrix. This is why the slicing is
901 done in Afw, its cleaner.
902 */
903
904
905 /* Inner limits; of mask box */
906 int maskStartCol = maskBox.getMinX();
907 int maskStartRow = maskBox.getMinY();
908 int maskEndCol = maskBox.getMaxX();
909 int maskEndRow = maskBox.getMaxY();
910
911 /*
912
913 |---------------------------|
914 | Kernel Boundary |
915 | |---------------------| |
916 | | Top | |
917 | |......_________......| |
918 | | | | | |
919 | | L | Box | R | |
920 | | | | | |
921 | |......---------......| |
922 | | Bottom | |
923 | |---------------------| |
924 | |
925 |---------------------------|
926
927 4 regions we want to extract from the pixels: top bottom left right
928
929 */
930 geom::Box2I tBox = geom::Box2I(geom::Point2I(startCol, maskEndRow + 1),
931 geom::Point2I(endCol, endRow));
932
933 geom::Box2I bBox = geom::Box2I(geom::Point2I(startCol, startRow),
934 geom::Point2I(endCol, maskStartRow - 1));
935
936 geom::Box2I lBox = geom::Box2I(geom::Point2I(startCol, maskStartRow),
937 geom::Point2I(maskStartCol - 1, maskEndRow));
938
939 geom::Box2I rBox = geom::Box2I(geom::Point2I(maskEndCol + 1, maskStartRow),
940 geom::Point2I(endCol, maskEndRow));
941
942 LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
943 "Upper good pixel region: %d,%d -> %d,%d",
944 tBox.getMinX(), tBox.getMinY(), tBox.getMaxX(), tBox.getMaxY());
945 LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
946 "Bottom good pixel region: %d,%d -> %d,%d",
947 bBox.getMinX(), bBox.getMinY(), bBox.getMaxX(), bBox.getMaxY());
948 LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
949 "Left good pixel region: %d,%d -> %d,%d",
950 lBox.getMinX(), lBox.getMinY(), lBox.getMaxX(), lBox.getMaxY());
951 LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
952 "Right good pixel region: %d,%d -> %d,%d",
953 rBox.getMinX(), rBox.getMinY(), rBox.getMaxX(), rBox.getMaxY());
954
956 boxArray.push_back(tBox);
957 boxArray.push_back(bBox);
958 boxArray.push_back(lBox);
959 boxArray.push_back(rBox);
960
961 int totalSize = tBox.getWidth() * tBox.getHeight();
962 totalSize += bBox.getWidth() * bBox.getHeight();
963 totalSize += lBox.getWidth() * lBox.getHeight();
964 totalSize += rBox.getWidth() * rBox.getHeight();
965
966 Eigen::MatrixXd eigenTemplate(totalSize, 1);
967 Eigen::MatrixXd eigenScience(totalSize, 1);
968 Eigen::MatrixXd eigeniVariance(totalSize, 1);
969 eigenTemplate.setZero();
970 eigenScience.setZero();
971 eigeniVariance.setZero();
972
973 boost::timer::cpu_timer t;
974
975 int nTerms = 0;
976 typename std::vector<geom::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 afwMath::ConvolutionControl convolutionControl;
1005 convolutionControl.setDoNormalize(false);
1006 for (kiter = basisList.begin(); kiter != basisList.end(); ++kiter, ++eiter) {
1007 afwMath::convolve(cimage, templateImage, **kiter, convolutionControl); /* 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 t.stop();
1029 double time = 1e-9 * t.elapsed().wall;
1030 LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
1031 "Total compute time to do basis convolutions : %.2f s", time);
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") {
1234 _lambda = estimateRisk(std::numeric_limits<double>::max());
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 ) :
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
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
1480 afwMath::Kernel::SpatialFunctionPtr> SpatialKernelSolution::getSolutionPair() {
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 =
1507 std::dynamic_pointer_cast<afwMath::LinearCombinationKernel>(_kernel)->getKernelList();
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 */
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
int end
#define LSST_EXCEPT_ADD(e, m)
Add the current location and a message to an existing exception before rethrowing it.
Definition Exception.h:54
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
afw::table::Key< afw::table::Array< ImagePixelT > > image
Image Subtraction helper functions.
Declaration of classes to store the solution for convolution kernels.
LSST DM logging module built on log4cxx.
#define LOGL_DEBUG(logger, message...)
Log a debug-level message using a varargs/printf style interface.
Definition Log.h:515
uint64_t * ptr
Definition RangeSet.cc:95
T begin(T... args)
Class to describe the properties of a detected object from an image.
Definition Footprint.h:63
A set of Footprints, associated with a MaskedImage.
std::shared_ptr< FootprintList > getFootprints()
: Return the Footprints of detected objects
A Threshold is used to pass a threshold value to detection algorithms.
Definition Threshold.h:43
lsst::geom::Box2I getBBox(ImageOrigin origin=PARENT) const
Definition ImageBase.h:445
lsst::geom::Extent2I getDimensions() const
Return the image's size; useful for passing to constructors.
Definition ImageBase.h:356
lsst::geom::Point2I getXY0() const
Return the image's origin.
Definition ImageBase.h:323
A class to represent a 2-dimensional array of pixels.
Definition Image.h:51
Represent a 2-dimensional array of bitmask pixels.
Definition Mask.h:77
void writeFits(std::string const &fileName, daf::base::PropertySet const *metadata=nullptr, std::string const &mode="w") const
Write a mask to a regular FITS file.
Parameters to control convolution.
void setDoNormalize(bool doNormalize)
void setParameters(std::vector< double > const &params)
Set all function parameters.
Definition Function.h:156
std::vector< double > const & getParameters() const noexcept
Return all function parameters.
Definition Function.h:129
lsst::geom::Extent2I const getDimensions() const
Return the Kernel's dimensions (width, height)
Definition Kernel.h:212
void setSpatialParameters(const std::vector< std::vector< double > > params)
Set the parameters of all spatial functions.
Definition Kernel.cc:110
double computeImage(lsst::afw::image::Image< Pixel > &image, bool doNormalize, double x=0.0, double y=0.0) const
Compute an image (pixellized representation of the kernel) in place.
Definition Kernel.cc:76
A kernel that is a linear combination of fixed basis kernels.
Definition Kernel.h:704
A class to evaluate image statistics.
Definition Statistics.h:222
double getValue(Property const prop=NOTHING) const
Return the value of the desired property (if specified in the constructor)
Class for storing generic metadata.
Definition PropertySet.h:66
An integer coordinate rectangle.
Definition Box.h:55
int getMinY() const noexcept
Definition Box.h:158
int getHeight() const noexcept
Definition Box.h:188
int getMinX() const noexcept
Definition Box.h:157
int getWidth() const noexcept
Definition Box.h:187
int getMaxX() const noexcept
Definition Box.h:161
int getMaxY() const noexcept
Definition Box.h:162
bool _fitForBackground
Background terms included in fit.
virtual double getConditionNumber(ConditionNumberType conditionType)
Eigen::VectorXd _bVec
Derived least squares B vector.
Eigen::MatrixXd const & getM()
Eigen::VectorXd _aVec
Derived least squares solution matrix.
KernelSolvedBy _solvedBy
Type of algorithm used to make solution.
Eigen::MatrixXd _mMat
Derived least squares M matrix.
static int _SolutionId
Unique identifier for solution.
lsst::afw::image::Image< lsst::afw::math::Kernel::Pixel > ImageT
virtual void buildWithMask(lsst::afw::image::Image< InputT > const &templateImage, lsst::afw::image::Image< InputT > const &scienceImage, lsst::afw::image::Image< lsst::afw::image::VariancePixel > const &varianceEstimate, lsst::afw::image::Mask< lsst::afw::image::MaskPixel > const &pixelMask)
virtual void buildSingleMaskOrig(lsst::afw::image::Image< InputT > const &templateImage, lsst::afw::image::Image< InputT > const &scienceImage, lsst::afw::image::Image< lsst::afw::image::VariancePixel > const &varianceEstimate, lsst::geom::Box2I maskBox)
MaskedKernelSolution(lsst::afw::math::KernelList const &basisList, bool fitForBackground)
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)
RegularizedKernelSolution(lsst::afw::math::KernelList const &basisList, bool fitForBackground, Eigen::MatrixXd const &hMat, lsst::daf::base::PropertySet const &ps)
std::pair< std::shared_ptr< lsst::afw::math::LinearCombinationKernel >, lsst::afw::math::Kernel::SpatialFunctionPtr > getSolutionPair()
std::shared_ptr< lsst::afw::image::Image< lsst::afw::math::Kernel::Pixel > > makeKernelImage(lsst::geom::Point2D const &pos)
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)
void addConstraint(float xCenter, float yCenter, Eigen::MatrixXd const &qMat, Eigen::VectorXd const &wVec)
virtual std::pair< std::shared_ptr< lsst::afw::math::Kernel >, double > getSolutionPair()
void _setKernel()
Set kernel after solution.
virtual std::shared_ptr< lsst::afw::math::Kernel > getKernel()
std::shared_ptr< lsst::afw::math::Kernel > _kernel
Derived single-object convolution kernel.
void _setKernelUncertainty()
Not implemented.
virtual std::shared_ptr< lsst::afw::image::Image< lsst::afw::math::Kernel::Pixel > > makeKernelImage()
StaticKernelSolution(lsst::afw::math::KernelList const &basisList, bool fitForBackground)
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)
Provides consistent interface for LSST exceptions.
Definition Exception.h:107
Reports invalid arguments.
Definition Runtime.h:66
T end(T... args)
T endl(T... args)
T isnan(T... args)
T make_pair(T... args)
Statistics makeStatistics(lsst::afw::image::Image< Pixel > const &img, lsst::afw::image::Mask< image::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl=StatisticsControl())
Handle a watered-down front-end to the constructor (no variance)
Definition Statistics.h:361
void convolve(OutImageT &convolvedImage, InImageT const &inImage, KernelT const &kernel, ConvolutionControl const &convolutionControl=ConvolutionControl())
Convolve an Image or MaskedImage with a Kernel, setting pixels of an existing output image.
@ MIN
estimate sample minimum
Definition Statistics.h:66
Eigen::MatrixXd imageToEigenMatrix(lsst::afw::image::Image< PixelT > const &img)
Turns a 2-d Image into a 2-d Eigen Matrix.
Eigen::MatrixXi maskToEigenMatrix(lsst::afw::image::Mask< lsst::afw::image::MaskPixel > const &mask)
T push_back(T... args)
T reserve(T... args)
T reset(T... args)
T size(T... args)
#define DEBUG_MATRIX
#define DEBUG_MATRIX2