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

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

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; }

lsst::ip::diffim::KernelSolution::_id

int _id

Unique ID for object.

Definition KernelSolution.h:72

◆ 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;}

lsst::ip::diffim::KernelSolution::_aVec

Eigen::VectorXd _aVec

Derived least squares solution matrix.

Definition KernelSolution.h:75

◆ 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::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::cpu_timer t;
 
        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");
                        }
                }
 
        t.stop();
        double time = 1e-9 * t.elapsed().wall;
        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;
    }

◆ solve() [2/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();
    }

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/lsst-scipipe-10.0.0/Linux64/ip_diffim/g208c678f98+627fe8cd4e/include/lsst/ip/diffim/KernelSolution.h
/j/snowflake/release/lsstsw/stack/lsst-scipipe-10.0.0/Linux64/ip_diffim/g208c678f98+627fe8cd4e/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