#include <KernelSolution.h>

Inheritance diagram for lsst::ip::diffim::SpatialKernelSolution:

Public Types
typedef std::shared_ptr< SpatialKernelSolution >	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::Pixel >	ImageT

Public Member Functions
	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)

virtual	~SpatialKernelSolution ()

void	addConstraint (float xCenter, float yCenter, Eigen::MatrixXd const &qMat, Eigen::VectorXd const &wVec)

void	solve ()

std::shared_ptr< lsst::afw::image::Image< lsst::afw::math::Kernel::Pixel > >	makeKernelImage (lsst::geom::Point2D const &pos)

std::pair< std::shared_ptr< lsst::afw::math::LinearCombinationKernel >, lsst::afw::math::Kernel::SpatialFunctionPtr >	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 Attributes
int	_id
	Unique ID for object. More...

Eigen::MatrixXd	_mMat
	Derived least squares M matrix. More...

Eigen::VectorXd	_bVec
	Derived least squares B vector. More...

Eigen::VectorXd	_aVec
	Derived least squares solution matrix. More...

KernelSolvedBy	_solvedBy
	Type of algorithm used to make solution. More...

bool	_fitForBackground
	Background terms included in fit. More...

Static Protected Attributes
static int	_SolutionId = 0
	Unique identifier for solution. More...

Detailed Description

Definition at line 174 of file KernelSolution.h.

Member Typedef Documentation

◆ ImageT

typedef lsst::afw::image::Image<lsst::afw::math::Kernel::Pixel> lsst::ip::diffim::KernelSolution::ImageT

inherited

Definition at line 35 of file KernelSolution.h.

◆ PixelT

typedef lsst::afw::math::Kernel::Pixel lsst::ip::diffim::KernelSolution::PixelT

inherited

Definition at line 34 of file KernelSolution.h.

◆ Ptr

typedef std::shared_ptr<SpatialKernelSolution> lsst::ip::diffim::SpatialKernelSolution::Ptr

Definition at line 176 of file KernelSolution.h.

Member Enumeration Documentation

◆ ConditionNumberType

enum lsst::ip::diffim::KernelSolution::ConditionNumberType

inherited

Enumerator
EIGENVALUE
SVD

Definition at line 45 of file KernelSolution.h.

                                  {
             EIGENVALUE = 0,
             SVD        = 1
         };

◆ KernelSolvedBy

enum lsst::ip::diffim::KernelSolution::KernelSolvedBy

inherited

Enumerator
NONE
CHOLESKY_LDLT
CHOLESKY_LLT
LU
EIGENVECTOR

Definition at line 37 of file KernelSolution.h.

                             {
             NONE          = 0,
             CHOLESKY_LDLT = 1,
             CHOLESKY_LLT  = 2,
             LU            = 3,
             EIGENVECTOR   = 4
         };

Constructor & Destructor Documentation

◆ SpatialKernelSolution()

lsst::ip::diffim::SpatialKernelSolution::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
	)

Definition at line 1283 of file KernelSolution.cc.

           :
         KernelSolution(),
         _spatialKernelFunction(spatialKernelFunction),
         _constantFirstTerm(false),
         _kernel(),
         _background(background),
         _kSum(0.0),
         _ps(ps.deepCopy()),
         _nbases(0),
         _nkt(0),
         _nbt(0),
         _nt(0) {
  
         bool isAlardLupton    = _ps->getAsString("kernelBasisSet") == "alard-lupton";
         bool usePca           = _ps->getAsBool("usePcaForSpatialKernel");
         if (isAlardLupton || usePca) {
             _constantFirstTerm = true;
         }
         this->_fitForBackground = _ps->getAsBool("fitForBackground");
  
         _nbases = basisList.size();
         _nkt = _spatialKernelFunction->getParameters().size();
         _nbt = _fitForBackground ? _background->getParameters().size() : 0;
         _nt  = 0;
         if (_constantFirstTerm) {
             _nt = (_nbases - 1) * _nkt + 1 + _nbt;
         } else {
             _nt = _nbases * _nkt + _nbt;
         }
  
         Eigen::MatrixXd mMat(_nt, _nt);
         Eigen::VectorXd bVec(_nt);
         mMat.setZero();
         bVec.setZero();
  
         this->_mMat = mMat;
         this->_bVec = bVec;
  
         _kernel = std::shared_ptr<afwMath::LinearCombinationKernel>(
             new afwMath::LinearCombinationKernel(basisList, *_spatialKernelFunction)
             );
  
         LOGL_DEBUG("TRACE3.ip.diffim.SpatialKernelSolution",
                    "Initializing with size %d %d %d and constant first term = %s",
                    _nkt, _nbt, _nt,
                    _constantFirstTerm ? "true" : "false");
  
     }

◆ ~SpatialKernelSolution()

virtual lsst::ip::diffim::SpatialKernelSolution::~SpatialKernelSolution ( )

inlinevirtual

Definition at line 185 of file KernelSolution.h.

185 {};

Member Function Documentation

◆ addConstraint()

void lsst::ip::diffim::SpatialKernelSolution::addConstraint	(	float	xCenter,
		float	yCenter,
		Eigen::MatrixXd const &	qMat,
		Eigen::VectorXd const &	wVec
	)

Definition at line 1337 of file KernelSolution.cc.

                                                                          {
  
         LOGL_DEBUG("TRACE5.ip.diffim.SpatialKernelSolution.addConstraint",
                    "Adding candidate at %f, %f", xCenter, yCenter);
  
         /* Calculate P matrices */
         /* Pure kernel terms */
         Eigen::VectorXd pK(_nkt);
         std::vector<double> paramsK = _spatialKernelFunction->getParameters();
         for (int idx = 0; idx < _nkt; idx++) { paramsK[idx] = 0.0; }
         for (int idx = 0; idx < _nkt; idx++) {
             paramsK[idx] = 1.0;
             _spatialKernelFunction->setParameters(paramsK);
             pK(idx) = (*_spatialKernelFunction)(xCenter, yCenter); /* Assume things don't vary over stamp */
             paramsK[idx] = 0.0;
         }
         Eigen::MatrixXd pKpKt = (pK * pK.transpose());
  
         Eigen::VectorXd pB;
         Eigen::MatrixXd pBpBt;
         Eigen::MatrixXd pKpBt;
         if (_fitForBackground) {
             pB = Eigen::VectorXd(_nbt);
  
             /* Pure background terms */
             std::vector<double> paramsB = _background->getParameters();
             for (int idx = 0; idx < _nbt; idx++) { paramsB[idx] = 0.0; }
             for (int idx = 0; idx < _nbt; idx++) {
                 paramsB[idx] = 1.0;
                 _background->setParameters(paramsB);
                 pB(idx) = (*_background)(xCenter, yCenter);       /* Assume things don't vary over stamp */
                 paramsB[idx] = 0.0;
             }
             pBpBt = (pB * pB.transpose());
  
             /* Cross terms */
             pKpBt = (pK * pB.transpose());
         }
  
         if (DEBUG_MATRIX) {
             std::cout << "Spatial weights" << std::endl;
             std::cout << "pKpKt " << pKpKt << std::endl;
             if (_fitForBackground) {
                 std::cout << "pBpBt " << pBpBt << std::endl;
                 std::cout << "pKpBt " << pKpBt << std::endl;
             }
         }
  
         if (DEBUG_MATRIX) {
             std::cout << "Spatial matrix inputs" << std::endl;
             std::cout << "M " << qMat << std::endl;
             std::cout << "B " << wVec << std::endl;
         }
  
         /* first index to start the spatial blocks; default=0 for non-constant first term */
         int m0 = 0;
         /* how many rows/cols to adjust the matrices/vectors; default=0 for non-constant first term */
         int dm = 0;
         /* where to start the background terms; this is always true */
         int mb = _nt - _nbt;
  
         if (_constantFirstTerm) {
             m0 = 1;       /* we need to manually fill in the first (non-spatial) terms below */
             dm = _nkt-1;  /* need to shift terms due to lack of spatial variation in first term */
  
             _mMat(0, 0) += qMat(0,0);
             for(int m2 = 1; m2 < _nbases; m2++)  {
                 _mMat.block(0, m2*_nkt-dm, 1, _nkt) += qMat(0,m2) * pK.transpose();
             }
             _bVec(0) += wVec(0);
  
             if (_fitForBackground) {
                 _mMat.block(0, mb, 1, _nbt) += qMat(0,_nbases) * pB.transpose();
             }
         }
  
         /* Fill in the spatial blocks */
         for(int m1 = m0; m1 < _nbases; m1++)  {
             /* Diagonal kernel-kernel term; only use upper triangular part of pKpKt */
             _mMat.block(m1*_nkt-dm, m1*_nkt-dm, _nkt, _nkt) +=
                 (pKpKt * qMat(m1,m1)).triangularView<Eigen::Upper>();
  
             /* Kernel-kernel terms */
             for(int m2 = m1+1; m2 < _nbases; m2++)  {
                 _mMat.block(m1*_nkt-dm, m2*_nkt-dm, _nkt, _nkt) += qMat(m1,m2) * pKpKt;
             }
  
             if (_fitForBackground) {
                 /* Kernel cross terms with background */
                 _mMat.block(m1*_nkt-dm, mb, _nkt, _nbt) += qMat(m1,_nbases) * pKpBt;
             }
  
             /* B vector */
             _bVec.segment(m1*_nkt-dm, _nkt) += wVec(m1) * pK;
         }
  
         if (_fitForBackground) {
             /* Background-background terms only */
             _mMat.block(mb, mb, _nbt, _nbt) +=
                 (pBpBt * qMat(_nbases,_nbases)).triangularView<Eigen::Upper>();
             _bVec.segment(mb, _nbt)         += wVec(_nbases) * pB;
         }
  
         if (DEBUG_MATRIX) {
             std::cout << "Spatial matrix outputs" << std::endl;
             std::cout << "mMat " << _mMat << std::endl;
             std::cout << "bVec " << _bVec << std::endl;
         }
  
     }

◆ getB()

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

inlineinherited

Definition at line 65 of file KernelSolution.h.

65 {return _bVec;}

◆ getConditionNumber() [1/2]

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

virtualinherited

Definition at line 94 of file KernelSolution.cc.

                                                                                {
         return getConditionNumber(_mMat, 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.

                                                                                  {
         switch (conditionType) {
         case EIGENVALUE:
             {
             Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eVecValues(mMat);
             Eigen::VectorXd eValues = eVecValues.eigenvalues();
             double eMax = eValues.maxCoeff();
             double eMin = eValues.minCoeff();
             LOGL_DEBUG("TRACE3.ip.diffim.KernelSolution.getConditionNumber",
                        "EIGENVALUE eMax / eMin = %.3e", eMax / eMin);
             return (eMax / eMin);
             break;
             }
         case SVD:
             {
             Eigen::VectorXd sValues = mMat.jacobiSvd().singularValues();
             double sMax = sValues.maxCoeff();
             double sMin = sValues.minCoeff();
             LOGL_DEBUG("TRACE3.ip.diffim.KernelSolution.getConditionNumber",
                        "SVD eMax / eMin = %.3e", sMax / sMin);
             return (sMax / sMin);
             break;
             }
         default:
             {
             throw LSST_EXCEPT(pexExcept::InvalidParameterError,
                               "Undefined ConditionNumberType : only EIGENVALUE, SVD allowed.");
             break;
             }
         }
     }

◆ getId()

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

inlineinherited

Definition at line 69 of file KernelSolution.h.

69 { return _id; }

◆ getM()

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

inlineinherited

Definition at line 64 of file KernelSolution.h.

64 {return _mMat;}

◆ getSolutionPair()

std::pair< std::shared_ptr< afwMath::LinearCombinationKernel >, afwMath::Kernel::SpatialFunctionPtr > lsst::ip::diffim::SpatialKernelSolution::getSolutionPair ( )

Definition at line 1480 of file KernelSolution.cc.

                                                                                     {
         if (_solvedBy == KernelSolution::NONE) {
             throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return solution");
         }
  
         return std::make_pair(_kernel, _background);
     }

◆ getSolvedBy()

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

inlineinherited

Definition at line 60 of file KernelSolution.h.

60 {return _solvedBy;}

◆ makeKernelImage()

std::shared_ptr< lsst::afw::image::Image< lsst::afw::math::Kernel::Pixel > > lsst::ip::diffim::SpatialKernelSolution::makeKernelImage ( lsst::geom::Point2D const & pos )

Definition at line 1450 of file KernelSolution.cc.

                                                                                                                           {
         if (_solvedBy == KernelSolution::NONE) {
             throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return image");
         }
         std::shared_ptr<afwImage::Image<afwMath::Kernel::Pixel>> image(
             new afwImage::Image<afwMath::Kernel::Pixel>(_kernel->getDimensions())
             );
         (void)_kernel->computeImage(*image, false, pos[0], pos[1]);
         return image;
     }

◆ printA()

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

inlineinherited

Definition at line 68 of file KernelSolution.h.

68 {std::cout << _aVec << std::endl;}

◆ printB()

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

inlineinherited

Definition at line 67 of file KernelSolution.h.

67 {std::cout << _bVec << std::endl;}

◆ printM()

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

inlineinherited

Definition at line 66 of file KernelSolution.h.

66 {std::cout << _mMat << std::endl;}

◆ solve() [1/2]

void lsst::ip::diffim::SpatialKernelSolution::solve ( )

virtual

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

Definition at line 1461 of file KernelSolution.cc.

                                       {
         /* Fill in the other half of mMat */
         for (int i = 0; i < _nt; i++) {
             for (int j = i+1; j < _nt; j++) {
                 _mMat(j,i) = _mMat(i,j);
             }
         }
  
         try {
             KernelSolution::solve();
         } catch (pexExcept::Exception &e) {
             LSST_EXCEPT_ADD(e, "Unable to solve spatial kernel matrix");
             throw e;
         }
         /* Turn matrices into _kernel and _background */
         _setKernel();
     }

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

                                                           {
  
         if (DEBUG_MATRIX) {
             std::cout << "M " << std::endl;
             std::cout << mMat << std::endl;
             std::cout << "B " << std::endl;
             std::cout << bVec << std::endl;
         }
  
         Eigen::VectorXd aVec = Eigen::VectorXd::Zero(bVec.size());
  
         boost::timer t;
         t.restart();
  
         LOGL_DEBUG("TRACE2.ip.diffim.KernelSolution.solve",
                    "Solving for kernel");
                 _solvedBy = LU;
                 Eigen::FullPivLU<Eigen::MatrixXd> lu(mMat);
                 if (lu.isInvertible()) {
                         aVec = lu.solve(bVec);
                 } else {
                         LOGL_DEBUG("TRACE3.ip.diffim.KernelSolution.solve",
                                    "Unable to determine kernel via LU");
                         /* LAST RESORT */
                         try {
  
                                 _solvedBy = EIGENVECTOR;
                                 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eVecValues(mMat);
                                 Eigen::MatrixXd const& rMat = eVecValues.eigenvectors();
                                 Eigen::VectorXd eValues = eVecValues.eigenvalues();
  
                                 for (int i = 0; i != eValues.rows(); ++i) {
                                         if (eValues(i) != 0.0) {
                                                 eValues(i) = 1.0/eValues(i);
                                         }
                                 }
  
                                 aVec = rMat * eValues.asDiagonal() * rMat.transpose() * bVec;
                         } catch (pexExcept::Exception& e) {
  
                                 _solvedBy = NONE;
                                 LOGL_DEBUG("TRACE3.ip.diffim.KernelSolution.solve",
                                            "Unable to determine kernel via eigen-values");
  
                                 throw LSST_EXCEPT(pexExcept::Exception, "Unable to determine kernel solution");
                         }
                 }
  
         double time = t.elapsed();
         LOGL_DEBUG("TRACE3.ip.diffim.KernelSolution.solve",
                    "Compute time for matrix math : %.2f s", time);
  
         if (DEBUG_MATRIX) {
                   std::cout << "A " << std::endl;
             std::cout << aVec << std::endl;
         }
  
         _aVec = aVec;
     }

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.

◆ _bVec

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

protectedinherited

Derived least squares B vector.

Definition at line 74 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.

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

/j/snowflake/release/lsstsw/stack/1a1d771/Linux64/ip_diffim/20.0.0/include/lsst/ip/diffim/KernelSolution.h
/j/snowflake/release/lsstsw/stack/1a1d771/Linux64/ip_diffim/20.0.0/src/KernelSolution.cc

Public Types

Public Member Functions

Protected Attributes

Static Protected Attributes

Detailed Description

Member Typedef Documentation

◆ ImageT

◆ PixelT

◆ Ptr

Member Enumeration Documentation

◆ ConditionNumberType

◆ KernelSolvedBy

Constructor & Destructor Documentation

◆ SpatialKernelSolution()

◆ ~SpatialKernelSolution()

Member Function Documentation

◆ addConstraint()

◆ getB()

◆ getConditionNumber() [1/2]

◆ getConditionNumber() [2/2]

◆ getId()

◆ getM()

◆ getSolutionPair()

◆ getSolvedBy()

◆ makeKernelImage()

◆ printA()

◆ printB()

◆ printM()

◆ solve() [1/2]

◆ solve() [2/2]

Member Data Documentation

◆ _aVec

◆ _bVec

◆ _fitForBackground

◆ _id

◆ _mMat

◆ _SolutionId

◆ _solvedBy