LSST Applications g0265f82a02+c6dfa2ddaf,g1162b98a3f+b2075782a9,g2079a07aa2+1b2e822518,g2bbee38e9b+c6dfa2ddaf,g337abbeb29+c6dfa2ddaf,g3ddfee87b4+a60788ef87,g50ff169b8f+2eb0e556e8,g52b1c1532d+90ebb246c7,g555ede804d+a60788ef87,g591dd9f2cf+ba8caea58f,g5ec818987f+864ee9cddb,g858d7b2824+9ee1ab4172,g876c692160+a40945ebb7,g8a8a8dda67+90ebb246c7,g8cdfe0ae6a+4fd9e222a8,g99cad8db69+5e309b7bc6,g9ddcbc5298+a1346535a5,ga1e77700b3+df8f93165b,ga8c6da7877+aa12a14d27,gae46bcf261+c6dfa2ddaf,gb0e22166c9+8634eb87fb,gb3f2274832+d0da15e3be,gba4ed39666+1ac82b564f,gbb8dafda3b+5dfd9c994b,gbeb006f7da+97157f9740,gc28159a63d+c6dfa2ddaf,gc86a011abf+9ee1ab4172,gcf0d15dbbd+a60788ef87,gdaeeff99f8+1cafcb7cd4,gdc0c513512+9ee1ab4172,ge79ae78c31+c6dfa2ddaf,geb67518f79+ba1859f325,geb961e4c1e+f9439d1e6f,gee10cc3b42+90ebb246c7,gf1cff7945b+9ee1ab4172,w.2024.12
LSST Data Management Base Package
Loading...
Searching...
No Matches
Public Types | Public Member Functions | Protected Member Functions | Protected Attributes | Static Protected Attributes | List of all members
lsst::ip::diffim::StaticKernelSolution< InputT > Class Template Reference

#include <KernelSolution.h>

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

Public Types

typedef std::shared_ptr< StaticKernelSolution< InputT > > Ptr
 
enum  KernelSolvedBy {
  NONE = 0 , CHOLESKY_LDLT = 1 , CHOLESKY_LLT = 2 , LU = 3 ,
  EIGENVECTOR = 4
}
 
enum  ConditionNumberType { EIGENVALUE = 0 , SVD = 1 }
 
typedef lsst::afw::math::Kernel::Pixel PixelT
 
typedef lsst::afw::image::Image< lsst::afw::math::Kernel::PixelImageT
 

Public Member Functions

 StaticKernelSolution (lsst::afw::math::KernelList const &basisList, bool fitForBackground)
 
virtual ~StaticKernelSolution ()
 
void solve ()
 
virtual void build (lsst::afw::image::Image< InputT > const &templateImage, lsst::afw::image::Image< InputT > const &scienceImage, lsst::afw::image::Image< lsst::afw::image::VariancePixel > const &varianceEstimate)
 
virtual std::shared_ptr< lsst::afw::math::KernelgetKernel ()
 
virtual std::shared_ptr< lsst::afw::image::Image< lsst::afw::math::Kernel::Pixel > > makeKernelImage ()
 
virtual double getBackground ()
 
virtual double getKsum ()
 
virtual std::pair< std::shared_ptr< lsst::afw::math::Kernel >, double > getSolutionPair ()
 
virtual void solve (Eigen::MatrixXd const &mMat, Eigen::VectorXd const &bVec)
 
KernelSolvedBy getSolvedBy ()
 
virtual double getConditionNumber (ConditionNumberType conditionType)
 
virtual double getConditionNumber (Eigen::MatrixXd const &mMat, ConditionNumberType conditionType)
 
Eigen::MatrixXd const & getM ()
 
Eigen::VectorXd const & getB ()
 
void printM ()
 
void printB ()
 
void printA ()
 
int getId () const
 

Protected Member Functions

void _setKernel ()
 Set kernel after solution.
 
void _setKernelUncertainty ()
 Not implemented.
 

Protected Attributes

Eigen::MatrixXd _cMat
 K_i x R.
 
Eigen::VectorXd _iVec
 Vectorized I.
 
Eigen::VectorXd _ivVec
 Inverse variance.
 
std::shared_ptr< lsst::afw::math::Kernel_kernel
 Derived single-object convolution kernel.
 
double _background
 Derived differential background estimate.
 
double _kSum
 Derived kernel sum.
 
int _id
 Unique ID for object.
 
Eigen::MatrixXd _mMat
 Derived least squares M matrix.
 
Eigen::VectorXd _bVec
 Derived least squares B vector.
 
Eigen::VectorXd _aVec
 Derived least squares solution matrix.
 
KernelSolvedBy _solvedBy
 Type of algorithm used to make solution.
 
bool _fitForBackground
 Background terms included in fit.
 

Static Protected Attributes

static int _SolutionId = 0
 Unique identifier for solution.
 

Detailed Description

template<typename InputT>
class lsst::ip::diffim::StaticKernelSolution< InputT >

Definition at line 83 of file KernelSolution.h.

Member Typedef Documentation

◆ ImageT

Definition at line 35 of file KernelSolution.h.

◆ PixelT

Definition at line 34 of file KernelSolution.h.

◆ Ptr

Definition at line 85 of file KernelSolution.h.

Member Enumeration Documentation

◆ ConditionNumberType

Enumerator
EIGENVALUE 
SVD 

Definition at line 45 of file KernelSolution.h.

◆ KernelSolvedBy

Enumerator
NONE 
CHOLESKY_LDLT 
CHOLESKY_LLT 
LU 
EIGENVECTOR 

Definition at line 37 of file KernelSolution.h.

Constructor & Destructor Documentation

◆ StaticKernelSolution()

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

Definition at line 195 of file KernelSolution.cc.

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 };
A kernel that is a linear combination of fixed basis kernels.
Definition Kernel.h:704
double _background
Derived differential background estimate.
Eigen::VectorXd _iVec
Vectorized I.
std::shared_ptr< lsst::afw::math::Kernel > _kernel
Derived single-object convolution kernel.
Eigen::VectorXd _ivVec
Inverse variance.

◆ ~StaticKernelSolution()

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

Definition at line 89 of file KernelSolution.h.

89{};

Member Function Documentation

◆ _setKernel()

template<typename InputT >
void lsst::ip::diffim::StaticKernelSolution< InputT >::_setKernel ( )
protected

Set kernel after solution.

Definition at line 423 of file KernelSolution.cc.

423 {
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 }
445
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 }
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
afw::table::Key< afw::table::Array< ImagePixelT > > image
lsst::geom::Extent2I const getDimensions() const
Return the Kernel's dimensions (width, height)
Definition Kernel.h:212
void setKernelParameters(std::vector< double > const &params)
Set the kernel parameters of a spatially invariant kernel.
Definition Kernel.h:341
double computeImage(lsst::afw::image::Image< Pixel > &image, bool doNormalize, double x=0.0, double y=0.0) const
Compute an image (pixellized representation of the kernel) in place.
Definition Kernel.cc:76
bool _fitForBackground
Background terms included in fit.
Eigen::VectorXd _aVec
Derived least squares solution matrix.
KernelSolvedBy _solvedBy
Type of algorithm used to make solution.
lsst::afw::image::Image< lsst::afw::math::Kernel::Pixel > ImageT
Provides consistent interface for LSST exceptions.
Definition Exception.h:107
T isnan(T... args)

◆ _setKernelUncertainty()

template<typename InputT >
void lsst::ip::diffim::StaticKernelSolution< InputT >::_setKernelUncertainty ( )
protected

Not implemented.

Definition at line 463 of file KernelSolution.cc.

463 {
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 }

◆ build()

template<typename InputT >
void lsst::ip::diffim::StaticKernelSolution< InputT >::build ( lsst::afw::image::Image< InputT > const & templateImage,
lsst::afw::image::Image< InputT > const & scienceImage,
lsst::afw::image::Image< lsst::afw::image::VariancePixel > const & varianceEstimate )
virtual

Definition at line 261 of file KernelSolution.cc.

265 {
266
267 afwMath::Statistics varStats = afwMath::makeStatistics(varianceEstimate, afwMath::MIN);
268 if (varStats.getValue(afwMath::MIN) < 0.0) {
270 "Error: variance less than 0.0");
271 }
272 if (varStats.getValue(afwMath::MIN) == 0.0) {
274 "Error: variance equals 0.0, cannot inverse variance weight");
275 }
276
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. */
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 }
#define LOGL_DEBUG(logger, message...)
Log a debug-level message using a varargs/printf style interface.
Definition Log.h:515
T begin(T... args)
lsst::geom::Box2I getBBox(ImageOrigin origin=PARENT) const
Definition ImageBase.h:445
lsst::geom::Extent2I getDimensions() const
Return the image's size; useful for passing to constructors.
Definition ImageBase.h:356
A class to represent a 2-dimensional array of pixels.
Definition Image.h:51
Parameters to control convolution.
void setDoNormalize(bool doNormalize)
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)
An integer coordinate rectangle.
Definition Box.h:55
int getMinY() const noexcept
Definition Box.h:158
int getMinX() const noexcept
Definition Box.h:157
int getMaxX() const noexcept
Definition Box.h:161
int getMaxY() const noexcept
Definition Box.h:162
Eigen::VectorXd _bVec
Derived least squares B vector.
Eigen::MatrixXd _mMat
Derived least squares M matrix.
T end(T... args)
Statistics makeStatistics(lsst::afw::image::Image< Pixel > const &img, lsst::afw::image::Mask< image::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl=StatisticsControl())
Handle a watered-down front-end to the constructor (no variance)
Definition Statistics.h: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.
T size(T... args)

◆ getB()

Eigen::VectorXd const & lsst::ip::diffim::KernelSolution::getB ( )
inlineinherited

Definition at line 65 of file KernelSolution.h.

65{return _bVec;}

◆ getBackground()

template<typename InputT >
double lsst::ip::diffim::StaticKernelSolution< InputT >::getBackground ( )
virtual

Definition at line 235 of file KernelSolution.cc.

235 {
237 throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return background");
238 }
239 return _background;
240 }

◆ getConditionNumber() [1/2]

double lsst::ip::diffim::KernelSolution::getConditionNumber ( ConditionNumberType conditionType)
virtualinherited

Definition at line 94 of file KernelSolution.cc.

94 {
95 return getConditionNumber(_mMat, conditionType);
96 }
virtual double getConditionNumber(ConditionNumberType conditionType)

◆ getConditionNumber() [2/2]

double lsst::ip::diffim::KernelSolution::getConditionNumber ( Eigen::MatrixXd const & mMat,
ConditionNumberType conditionType )
virtualinherited

Definition at line 98 of file KernelSolution.cc.

99 {
100 switch (conditionType) {
101 case EIGENVALUE:
102 {
103 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eVecValues(mMat);
104 Eigen::VectorXd eValues = eVecValues.eigenvalues();
105 double eMax = eValues.maxCoeff();
106 double eMin = eValues.minCoeff();
107 LOGL_DEBUG("TRACE3.ip.diffim.KernelSolution.getConditionNumber",
108 "EIGENVALUE eMax / eMin = %.3e", eMax / eMin);
109 return (eMax / eMin);
110 break;
111 }
112 case SVD:
113 {
114 Eigen::VectorXd sValues = mMat.jacobiSvd().singularValues();
115 double sMax = sValues.maxCoeff();
116 double sMin = sValues.minCoeff();
117 LOGL_DEBUG("TRACE3.ip.diffim.KernelSolution.getConditionNumber",
118 "SVD eMax / eMin = %.3e", sMax / sMin);
119 return (sMax / sMin);
120 break;
121 }
122 default:
123 {
125 "Undefined ConditionNumberType : only EIGENVALUE, SVD allowed.");
126 break;
127 }
128 }
129 }
Reports invalid arguments.
Definition Runtime.h:66

◆ getId()

int lsst::ip::diffim::KernelSolution::getId ( ) const
inlineinherited

Definition at line 69 of file KernelSolution.h.

69{ return _id; }
int _id
Unique ID for object.

◆ getKernel()

Definition at line 215 of file KernelSolution.cc.

215 {
217 throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return solution");
218 }
219 return _kernel;
220 }

◆ getKsum()

template<typename InputT >
double lsst::ip::diffim::StaticKernelSolution< InputT >::getKsum ( )
virtual

Definition at line 243 of file KernelSolution.cc.

243 {
245 throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return ksum");
246 }
247 return _kSum;
248 }

◆ getM()

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

Definition at line 64 of file KernelSolution.h.

64{return _mMat;}

◆ getSolutionPair()

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

Definition at line 252 of file KernelSolution.cc.

252 {
254 throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return solution");
255 }
256
258 }
T make_pair(T... args)

◆ getSolvedBy()

KernelSolvedBy lsst::ip::diffim::KernelSolution::getSolvedBy ( )
inlineinherited

Definition at line 60 of file KernelSolution.h.

60{return _solvedBy;}

◆ makeKernelImage()

Definition at line 223 of file KernelSolution.cc.

223 {
225 throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return image");
226 }
229 );
230 (void)_kernel->computeImage(*image, false);
231 return image;
232 }

◆ printA()

void lsst::ip::diffim::KernelSolution::printA ( )
inlineinherited

Definition at line 68 of file KernelSolution.h.

T endl(T... args)

◆ printB()

void lsst::ip::diffim::KernelSolution::printB ( )
inlineinherited

Definition at line 67 of file KernelSolution.h.

◆ printM()

void lsst::ip::diffim::KernelSolution::printM ( )
inlineinherited

Definition at line 66 of file KernelSolution.h.

◆ solve() [1/2]

template<typename InputT >
void lsst::ip::diffim::StaticKernelSolution< InputT >::solve ( )
virtual

Reimplemented from lsst::ip::diffim::KernelSolution.

Definition at line 397 of file KernelSolution.cc.

397 {
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;
406 std::cout << "iV" << std::endl;
408 std::cout << "I" << 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 }
#define LSST_EXCEPT_ADD(e, m)
Add the current location and a message to an existing exception before rethrowing it.
Definition Exception.h:54
void _setKernel()
Set kernel after solution.
#define DEBUG_MATRIX

◆ solve() [2/2]

void lsst::ip::diffim::KernelSolution::solve ( Eigen::MatrixXd const & mMat,
Eigen::VectorXd const & bVec )
virtualinherited

Definition at line 131 of file KernelSolution.cc.

132 {
133
134 if (DEBUG_MATRIX) {
135 std::cout << "M " << std::endl;
136 std::cout << mMat << std::endl;
137 std::cout << "B " << std::endl;
138 std::cout << bVec << std::endl;
139 }
140
141 Eigen::VectorXd aVec = Eigen::VectorXd::Zero(bVec.size());
142
143 boost::timer::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 }

Member Data Documentation

◆ _aVec

Eigen::VectorXd lsst::ip::diffim::KernelSolution::_aVec
protectedinherited

Derived least squares solution matrix.

Definition at line 75 of file KernelSolution.h.

◆ _background

template<typename InputT >
double lsst::ip::diffim::StaticKernelSolution< InputT >::_background
protected

Derived differential background estimate.

Definition at line 110 of file KernelSolution.h.

◆ _bVec

Eigen::VectorXd lsst::ip::diffim::KernelSolution::_bVec
protectedinherited

Derived least squares B vector.

Definition at line 74 of file KernelSolution.h.

◆ _cMat

template<typename InputT >
Eigen::MatrixXd lsst::ip::diffim::StaticKernelSolution< InputT >::_cMat
protected

K_i x R.

Definition at line 105 of file KernelSolution.h.

◆ _fitForBackground

bool lsst::ip::diffim::KernelSolution::_fitForBackground
protectedinherited

Background terms included in fit.

Definition at line 77 of file KernelSolution.h.

◆ _id

int lsst::ip::diffim::KernelSolution::_id
protectedinherited

Unique ID for object.

Definition at line 72 of file KernelSolution.h.

◆ _iVec

template<typename InputT >
Eigen::VectorXd lsst::ip::diffim::StaticKernelSolution< InputT >::_iVec
protected

Vectorized I.

Definition at line 106 of file KernelSolution.h.

◆ _ivVec

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

Inverse variance.

Definition at line 107 of file KernelSolution.h.

◆ _kernel

Derived single-object convolution kernel.

Definition at line 109 of file KernelSolution.h.

◆ _kSum

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

Derived kernel sum.

Definition at line 111 of file KernelSolution.h.

◆ _mMat

Eigen::MatrixXd lsst::ip::diffim::KernelSolution::_mMat
protectedinherited

Derived least squares M matrix.

Definition at line 73 of file KernelSolution.h.

◆ _SolutionId

int lsst::ip::diffim::KernelSolution::_SolutionId = 0
staticprotectedinherited

Unique identifier for solution.

Definition at line 78 of file KernelSolution.h.

◆ _solvedBy

KernelSolvedBy lsst::ip::diffim::KernelSolution::_solvedBy
protectedinherited

Type of algorithm used to make solution.

Definition at line 76 of file KernelSolution.h.


The documentation for this class was generated from the following files: