lsst::ip::diffim::SpatialKernelSolution Class Reference

#include <KernelSolution.h>

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

Public Types
typedef boost::shared_ptr < SpatialKernelSolution >	Ptr
Public Types inherited from lsst::ip::diffim::KernelSolution
enum	KernelSolvedBy {   NONE = 0, CHOLESKY_LDLT = 1, CHOLESKY_LLT = 2, LU = 3,   EIGENVECTOR = 4 }
enum	ConditionNumberType { EIGENVALUE = 0, SVD = 1 }
typedef boost::shared_ptr < KernelSolution >	Ptr
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::pex::policy::Policy policy)
virtual	~SpatialKernelSolution ()
void	addConstraint (float xCenter, float yCenter, boost::shared_ptr< Eigen::MatrixXd > qMat, boost::shared_ptr< Eigen::VectorXd > wVec)
void	solve ()
lsst::afw::image::Image < lsst::afw::math::Kernel::Pixel > ::Ptr	makeKernelImage (lsst::afw::geom::Point2D const &pos)
std::pair < lsst::afw::math::LinearCombinationKernel::Ptr, lsst::afw::math::Kernel::SpatialFunctionPtr >	getSolutionPair ()
Public Member Functions inherited from lsst::ip::diffim::KernelSolution
	KernelSolution (boost::shared_ptr< Eigen::MatrixXd > mMat, boost::shared_ptr< Eigen::VectorXd > bVec, bool fitForBackground)
	KernelSolution (bool fitForBackground)
	KernelSolution ()
virtual	~KernelSolution ()
virtual void	solve (Eigen::MatrixXd mMat, Eigen::VectorXd bVec)
KernelSolvedBy	getSolvedBy ()
virtual double	getConditionNumber (ConditionNumberType conditionType)
virtual double	getConditionNumber (Eigen::MatrixXd mMat, ConditionNumberType conditionType)
boost::shared_ptr < Eigen::MatrixXd >	getM ()
boost::shared_ptr < Eigen::VectorXd >	getB ()
void	printM ()
void	printB ()
void	printA ()
int	getId () const

Private Member Functions
void	_setKernel ()
	Set kernel after solution. More...
void	_setKernelUncertainty ()
	Not implemented. More...

Private Attributes
lsst::afw::math::Kernel::SpatialFunctionPtr	_spatialKernelFunction
	Spatial function for Kernel. More...
bool	_constantFirstTerm
	Is the first term constant. More...
lsst::afw::math::LinearCombinationKernel::Ptr	_kernel
	Spatial convolution kernel. More...
lsst::afw::math::Kernel::SpatialFunctionPtr	_background
	Spatial background model. More...
double	_kSum
	Derived kernel sum. More...
lsst::pex::policy::Policy	_policy
	Policy to control processing. More...
int	_nbases
	Number of basis functions. More...
int	_nkt
	Number of kernel terms. More...
int	_nbt
	Number of background terms. More...
int	_nt
	Total number of terms. More...

Additional Inherited Members
Protected Attributes inherited from lsst::ip::diffim::KernelSolution
int	_id
	Unique ID for object. More...
boost::shared_ptr < Eigen::MatrixXd >	_mMat
	Derived least squares M matrix. More...
boost::shared_ptr < Eigen::VectorXd >	_bVec
	Derived least squares B vector. More...
boost::shared_ptr < 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 inherited from lsst::ip::diffim::KernelSolution
static int	_SolutionId = 0
	Unique identifier for solution. More...

Detailed Description

Definition at line 173 of file KernelSolution.h.

Member Typedef Documentation

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

Definition at line 175 of file KernelSolution.h.

Constructor & Destructor Documentation

lsst::ip::diffim::SpatialKernelSolution::SpatialKernelSolution	(	lsst::afw::math::KernelList const &	basisList,
		lsst::afw::math::Kernel::SpatialFunctionPtr	spatialKernelFunction,
		lsst::afw::math::Kernel::SpatialFunctionPtr	background,
		lsst::pex::policy::Policy	policy
	)

Definition at line 1333 of file KernelSolution.cc.

          :
        KernelSolution(),
        _spatialKernelFunction(spatialKernelFunction),
        _constantFirstTerm(false),
        _kernel(),
        _background(background),
        _kSum(0.0),
        _policy(policy),
        _nbases(0),
        _nkt(0),
        _nbt(0),
        _nt(0) {

        bool isAlardLupton    = _policy.getString("kernelBasisSet") == "alard-lupton";
        bool usePca           = _policy.getBool("usePcaForSpatialKernel");
        if (isAlardLupton || usePca) {
            _constantFirstTerm = true;
        }
        this->_fitForBackground = _policy.getBool("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;
        }
        
        boost::shared_ptr<Eigen::MatrixXd> mMat (new Eigen::MatrixXd(_nt, _nt));
        boost::shared_ptr<Eigen::VectorXd> bVec (new Eigen::VectorXd(_nt));
        (*mMat).setZero();
        (*bVec).setZero();

        this->_mMat = mMat;
        this->_bVec = bVec;

        _kernel = afwMath::LinearCombinationKernel::Ptr(
            new afwMath::LinearCombinationKernel(basisList, *_spatialKernelFunction)
            );

        pexLog::TTrace<5>("lsst.ip.diffim.SpatialKernelSolution", 
                          "Initializing with size %d %d %d and constant first term = %s",
                          _nkt, _nbt, _nt,
                          _constantFirstTerm ? "true" : "false");
        
    }

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

inlinevirtual

Definition at line 184 of file KernelSolution.h.

{};

Member Function Documentation

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

private

Set kernel after solution.

Definition at line 1538 of file KernelSolution.cc.

                                           {
        /* Report on condition number */
        double cNumber = this->getConditionNumber(EIGENVALUE);
        
        if (_nkt == 1) {
            /* Not spatially varying; this fork is a specialization for convolution speed--up */
            
            /* Set the basis coefficients */
            std::vector<double> kCoeffs(_nbases);
            for (int i = 0; i < _nbases; i++) {
                if (std::isnan((*_aVec)(i))) {
                    throw LSST_EXCEPT(
                        pexExcept::Exception, 
                        str(boost::format(
                                "I. Unable to determine spatial kernel solution %d (nan).  Condition number = %.3e") % i % cNumber));
                }
                kCoeffs[i] = (*_aVec)(i);
            }
            lsst::afw::math::KernelList basisList = 
                boost::dynamic_pointer_cast<afwMath::LinearCombinationKernel>(_kernel)->getKernelList();
            _kernel.reset(
                new afwMath::LinearCombinationKernel(basisList, kCoeffs)
                );
        }
        else {
            
            /* Set the kernel coefficients */
            std::vector<std::vector<double> > kCoeffs;
            kCoeffs.reserve(_nbases);
            for (int i = 0, idx = 0; i < _nbases; i++) {
                kCoeffs.push_back(std::vector<double>(_nkt));

                /* Deal with the possibility the first term doesn't vary spatially */
                if ((i == 0) && (_constantFirstTerm)) {
                    if (std::isnan((*_aVec)(idx))) {
                        throw LSST_EXCEPT(
                            pexExcept::Exception, 
                            str(boost::format(
                                    "II. Unable to determine spatial kernel solution %d (nan).  Condition number = %.3e") % idx % cNumber));
                    }
                    kCoeffs[i][0] = (*_aVec)(idx++);
                }
                else {
                    for (int j = 0; j < _nkt; j++) {
                        if (std::isnan((*_aVec)(idx))) {
                            throw LSST_EXCEPT(
                                pexExcept::Exception, 
                                str(boost::format(
                                        "III. Unable to determine spatial kernel solution %d (nan).  Condition number = %.3e") % idx % cNumber));
                        }
                        kCoeffs[i][j] = (*_aVec)(idx++);
                    }
                }
            }
            _kernel->setSpatialParameters(kCoeffs);
        }

        /* Set kernel Sum */
        ImageT::Ptr image (new ImageT(_kernel->getDimensions()));
        _kSum  = _kernel->computeImage(*image, false);              

        /* Set the background coefficients */
        std::vector<double> bgCoeffs(_fitForBackground ? _nbt : 1);
        if (_fitForBackground) {
            for (int i = 0; i < _nbt; i++) {
                int idx = _nt - _nbt + i;
                if (std::isnan((*_aVec)(idx))) {
                    throw LSST_EXCEPT(pexExcept::Exception, 
                                      str(boost::format(
                           "Unable to determine spatial background solution %d (nan)") % (idx)));
                }
                bgCoeffs[i] = (*_aVec)(idx);
            }
        }
        else {
            bgCoeffs[0] = 0.;
        }
        _background->setParameters(bgCoeffs);
    }

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

private

Not implemented.

void lsst::ip::diffim::SpatialKernelSolution::addConstraint	(	float	xCenter,
		float	yCenter,
		boost::shared_ptr< Eigen::MatrixXd >	qMat,
		boost::shared_ptr< Eigen::VectorXd >	wVec
	)

Definition at line 1387 of file KernelSolution.cc.

                                                                                   {
        
        pexLog::TTrace<8>("lsst.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;
        }

    }

std::pair< afwMath::LinearCombinationKernel::Ptr, afwMath::Kernel::SpatialFunctionPtr > lsst::ip::diffim::SpatialKernelSolution::getSolutionPair

(

)

Definition at line 1530 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);
    }

lsst::afw::image::Image< lsst::afw::math::Kernel::Pixel >::Ptr lsst::ip::diffim::SpatialKernelSolution::makeKernelImage

(

lsst::afw::geom::Point2D const &

pos

)

Definition at line 1500 of file KernelSolution.cc.

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

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

virtual

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

Definition at line 1511 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

lsst::afw::math::Kernel::SpatialFunctionPtr lsst::ip::diffim::SpatialKernelSolution::_background

private

Spatial background model.

Definition at line 200 of file KernelSolution.h.

bool lsst::ip::diffim::SpatialKernelSolution::_constantFirstTerm

private

Is the first term constant.

Definition at line 197 of file KernelSolution.h.

lsst::afw::math::LinearCombinationKernel::Ptr lsst::ip::diffim::SpatialKernelSolution::_kernel

private

Spatial convolution kernel.

Definition at line 199 of file KernelSolution.h.

double lsst::ip::diffim::SpatialKernelSolution::_kSum

private

Derived kernel sum.

Definition at line 201 of file KernelSolution.h.

int lsst::ip::diffim::SpatialKernelSolution::_nbases

private

Number of basis functions.

Definition at line 204 of file KernelSolution.h.

int lsst::ip::diffim::SpatialKernelSolution::_nbt

private

Number of background terms.

Definition at line 206 of file KernelSolution.h.

int lsst::ip::diffim::SpatialKernelSolution::_nkt

private

Number of kernel terms.

Definition at line 205 of file KernelSolution.h.

int lsst::ip::diffim::SpatialKernelSolution::_nt

private

Total number of terms.

Definition at line 207 of file KernelSolution.h.

lsst::pex::policy::Policy lsst::ip::diffim::SpatialKernelSolution::_policy

private

Policy to control processing.

Definition at line 203 of file KernelSolution.h.

lsst::afw::math::Kernel::SpatialFunctionPtr lsst::ip::diffim::SpatialKernelSolution::_spatialKernelFunction

private

Spatial function for Kernel.

Definition at line 196 of file KernelSolution.h.

---

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

• /lsst/home/lsstsw/stack/Linux64/ip_diffim/master-gcec1da163f+63/include/lsst/ip/diffim/KernelSolution.h
• /lsst/home/lsstsw/stack/Linux64/ip_diffim/master-gcec1da163f+63/src/KernelSolution.cc