lsst::ip::diffim::StaticKernelSolution< InputT > Class Template Reference

#include <KernelSolution.h>

Inheritance diagram for lsst::ip::diffim::StaticKernelSolution< InputT >:

Public Types
typedef boost::shared_ptr < StaticKernelSolution< InputT > >	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
	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 lsst::afw::math::Kernel::Ptr	getKernel ()
virtual lsst::afw::image::Image < lsst::afw::math::Kernel::Pixel > ::Ptr	makeKernelImage ()
virtual double	getBackground ()
virtual double	getKsum ()
virtual std::pair < boost::shared_ptr < lsst::afw::math::Kernel > , double >	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

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

Protected Attributes
boost::shared_ptr < Eigen::MatrixXd >	_cMat
	K_i x R. More...
boost::shared_ptr < Eigen::VectorXd >	_iVec
	Vectorized I. More...
boost::shared_ptr < Eigen::VectorXd >	_ivVec
	Inverse variance. More...
lsst::afw::math::Kernel::Ptr	_kernel
	Derived single-object convolution kernel. More...
double	_background
	Derived differential background estimate. More...
double	_kSum
	Derived kernel sum. More...
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...

Additional Inherited Members
Static Protected Attributes inherited from lsst::ip::diffim::KernelSolution
static int	_SolutionId = 0
	Unique identifier for solution. More...

Detailed Description

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

Definition at line 82 of file KernelSolution.h.

Member Typedef Documentation

template<typename InputT >

typedef boost::shared_ptr<StaticKernelSolution<InputT> > lsst::ip::diffim::StaticKernelSolution< InputT >::Ptr

Definition at line 84 of file KernelSolution.h.

Constructor & Destructor Documentation

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.

        :
        KernelSolution(fitForBackground),
        _cMat(),
        _iVec(),
        _ivVec(),
        _kernel(),
        _background(0.0),
        _kSum(0.0)
    {
        std::vector<double> kValues(basisList.size());
        _kernel = boost::shared_ptr<afwMath::Kernel>( 
            new afwMath::LinearCombinationKernel(basisList, kValues) 
            );
    };

template<typename InputT >

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

inlinevirtual

Definition at line 88 of file KernelSolution.h.

{};

Member Function Documentation

template<typename InputT >

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

protected

Set kernel after solution.

Definition at line 433 of file KernelSolution.cc.

                                                  {
        if (_solvedBy == KernelSolution::NONE) {
            throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot make solution");
        }

        unsigned int const nParameters           = _aVec->size();
        unsigned int const nBackgroundParameters = _fitForBackground ? 1 : 0;
        unsigned int const nKernelParameters     = 
            boost::dynamic_pointer_cast<afwMath::LinearCombinationKernel>(_kernel)->getKernelList().size();
        if (nParameters != (nKernelParameters + nBackgroundParameters)) 
            throw LSST_EXCEPT(pexExcept::Exception, "Mismatched sizes in kernel solution");

        /* Fill in the kernel results */
        std::vector<double> kValues(nKernelParameters);
        for (unsigned int idx = 0; idx < nKernelParameters; idx++) {
            if (std::isnan((*_aVec)(idx))) {
                throw LSST_EXCEPT(pexExcept::Exception, 
                                  str(boost::format("Unable to determine kernel solution %d (nan)") % idx));
            }
            kValues[idx] = (*_aVec)(idx);
        }
        _kernel->setKernelParameters(kValues);

        ImageT::Ptr image (
            new ImageT(_kernel->getDimensions())
            );
        _kSum  = _kernel->computeImage(*image, false);              
        
        if (_fitForBackground) {
            if (std::isnan((*_aVec)(nParameters-1))) {
                throw LSST_EXCEPT(pexExcept::Exception, 
                                  str(boost::format("Unable to determine background solution %d (nan)") % 
                                      (nParameters-1)));
            }
            _background = (*_aVec)(nParameters-1);
        }
    }        

template<typename InputT >

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

protected

Not implemented.

Definition at line 473 of file KernelSolution.cc.

                                                             {
        throw LSST_EXCEPT(pexExcept::Exception, "Uncertainty calculation not supported");

        /* Estimate of parameter uncertainties comes from the inverse of the
         * covariance matrix (noise spectrum).  
         * N.R. 15.4.8 to 15.4.15
         * 
         * Since this is a linear problem no need to use Fisher matrix
         * N.R. 15.5.8
         *
         * Although I might be able to take advantage of the solution above.
         * Since this now works and is not the rate limiting step, keep as-is for DC3a.
         *
         * Use Cholesky decomposition again.
         * Cholkesy:
         * Cov       =  L L^t
         * Cov^(-1)  = (L L^t)^(-1)
         *           = (L^T)^-1 L^(-1)
         * 
         * Code would be:
         * 
         * Eigen::MatrixXd             cov    = (*_mMat).transpose() * (*_mMat);
         * Eigen::LLT<Eigen::MatrixXd> llt    = cov.llt();
         * Eigen::MatrixXd             error2 = llt.matrixL().transpose().inverse()*llt.matrixL().inverse();
         */
    }

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.

          {

        afwMath::Statistics varStats = afwMath::makeStatistics(varianceEstimate, afwMath::MIN);
        if (varStats.getValue(afwMath::MIN) < 0.0) {
            throw LSST_EXCEPT(pexExcept::Exception, 
                              "Error: variance less than 0.0");
        }
        if (varStats.getValue(afwMath::MIN) == 0.0) {
            throw LSST_EXCEPT(pexExcept::Exception, 
                              "Error: variance equals 0.0, cannot inverse variance weight");
        }

        lsst::afw::math::KernelList basisList = 
            boost::dynamic_pointer_cast<afwMath::LinearCombinationKernel>(_kernel)->getKernelList();
        
        unsigned int const nKernelParameters     = basisList.size();
        unsigned int const nBackgroundParameters = _fitForBackground ? 1 : 0;
        unsigned int const nParameters           = nKernelParameters + nBackgroundParameters;

        std::vector<boost::shared_ptr<afwMath::Kernel> >::const_iterator kiter = basisList.begin();
        
        /* Ignore buffers around edge of convolved images :
         * 
         * If the kernel has width 5, it has center pixel 2.  The first good pixel
         * is the (5-2)=3rd pixel, which is array index 2, and ends up being the
         * index of the central pixel.
         * 
         * You also have a buffer of unusable pixels on the other side, numbered
         * width-center-1.  The last good usable pixel is N-width+center+1.
         * 
         * Example : the kernel is width = 5, center = 2
         * 
         * ---|---|-c-|---|---|
         * 
         * the image is width = N
         * convolve this with the kernel, and you get
         * 
         * |-x-|-x-|-g-|---|---| ... |---|---|-g-|-x-|-x-|
         * 
         * g = first/last good pixel
         * x = bad
         * 
         * the first good pixel is the array index that has the value "center", 2
         * the last good pixel has array index N-(5-2)+1
         * eg. if N = 100, you want to use up to index 97
         * 100-3+1 = 98, and the loops use i < 98, meaning the last
         * index you address is 97.
         */

        /* NOTE - we are accessing particular elements of Eigen arrays using
           these coordinates, therefore they need to be in LOCAL coordinates.
           This was written before ndarray unification.
        */
        afwGeom::Box2I goodBBox = (*kiter)->shrinkBBox(templateImage.getBBox(afwImage::LOCAL));
        unsigned int const startCol = goodBBox.getMinX();
        unsigned int const startRow = goodBBox.getMinY();
        // endCol/Row is one past the index of the last good col/row
        unsigned int endCol = goodBBox.getMaxX() + 1;
        unsigned int endRow = goodBBox.getMaxY() + 1;

        boost::timer t;
        t.restart();
        
        /* Eigen representation of input images; only the pixels that are unconvolved in cimage below */
        Eigen::MatrixXd eigenTemplate = imageToEigenMatrix(templateImage).block(startRow, 
                                                                                startCol, 
                                                                                endRow-startRow, 
                                                                                endCol-startCol);
        Eigen::MatrixXd eigenScience = imageToEigenMatrix(scienceImage).block(startRow, 
                                                                              startCol, 
                                                                              endRow-startRow, 
                                                                              endCol-startCol);
        Eigen::MatrixXd eigeniVariance = imageToEigenMatrix(varianceEstimate).block(
            startRow, startCol, endRow-startRow, endCol-startCol
        ).array().inverse().matrix();

        /* Resize into 1-D for later usage */
        eigenTemplate.resize(eigenTemplate.rows()*eigenTemplate.cols(), 1);
        eigenScience.resize(eigenScience.rows()*eigenScience.cols(), 1);
        eigeniVariance.resize(eigeniVariance.rows()*eigeniVariance.cols(), 1);
        
        /* Holds image convolved with basis function */
        afwImage::Image<PixelT> cimage(templateImage.getDimensions());
        
        /* Holds eigen representation of image convolved with all basis functions */
        std::vector<boost::shared_ptr<Eigen::MatrixXd> > convolvedEigenList(nKernelParameters);
        
        /* Iterators over convolved image list and basis list */
        typename std::vector<boost::shared_ptr<Eigen::MatrixXd> >::iterator eiter = 
            convolvedEigenList.begin();
        /* Create C_i in the formalism of Alard & Lupton */
        for (kiter = basisList.begin(); kiter != basisList.end(); ++kiter, ++eiter) {
            afwMath::convolve(cimage, templateImage, **kiter, false); /* cimage stores convolved image */

            boost::shared_ptr<Eigen::MatrixXd> cMat (
                new Eigen::MatrixXd(imageToEigenMatrix(cimage).block(startRow, 
                                                                     startCol, 
                                                                     endRow-startRow, 
                                                                     endCol-startCol))
                );
            cMat->resize(cMat->rows()*cMat->cols(), 1);
            *eiter = cMat;

        } 

        double time = t.elapsed();
        pexLog::TTrace<5>("lsst.ip.diffim.StaticKernelSolution.build", 
                          "Total compute time to do basis convolutions : %.2f s", time);
        t.restart();
        
        /* 
           Load matrix with all values from convolvedEigenList : all images
           (eigeniVariance, convolvedEigenList) must be the same size
        */
        Eigen::MatrixXd cMat(eigenTemplate.col(0).size(), nParameters);
        typename std::vector<boost::shared_ptr<Eigen::MatrixXd> >::iterator eiterj = 
            convolvedEigenList.begin();
        typename std::vector<boost::shared_ptr<Eigen::MatrixXd> >::iterator eiterE = 
            convolvedEigenList.end();
        for (unsigned int kidxj = 0; eiterj != eiterE; eiterj++, kidxj++) {
            cMat.col(kidxj) = (*eiterj)->col(0);
        }
        /* Treat the last "image" as all 1's to do the background calculation. */
        if (_fitForBackground)
            cMat.col(nParameters-1).fill(1.);

        _cMat.reset(new Eigen::MatrixXd(cMat));
        _ivVec.reset(new Eigen::VectorXd(eigeniVariance.col(0)));
        _iVec.reset(new Eigen::VectorXd(eigenScience.col(0)));

        /* Make these outside of solve() so I can check condition number */
        _mMat.reset(new Eigen::MatrixXd((*_cMat).transpose() * ((*_ivVec).asDiagonal() * (*_cMat))));
        _bVec.reset(new Eigen::VectorXd((*_cMat).transpose() * ((*_ivVec).asDiagonal() * (*_iVec))));
    }

template<typename InputT >

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

virtual

Definition at line 235 of file KernelSolution.cc.

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

template<typename InputT >

lsst::afw::math::Kernel::Ptr lsst::ip::diffim::StaticKernelSolution< InputT >::getKernel ( )

virtual

Definition at line 215 of file KernelSolution.cc.

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

template<typename InputT >

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

virtual

Definition at line 243 of file KernelSolution.cc.

                                                 {
        if (_solvedBy == KernelSolution::NONE) {
            throw LSST_EXCEPT(pexExcept::Exception, "Kernel not solved; cannot return ksum");
        }
        return _kSum;
    }

template<typename InputT >

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

virtual

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

template<typename InputT >

lsst::afw::image::Image< lsst::afw::math::Kernel::Pixel >::Ptr lsst::ip::diffim::StaticKernelSolution< InputT >::makeKernelImage ( )

virtual

Definition at line 223 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);              
        return image;
    }

template<typename InputT >

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

virtual

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

Reimplemented in lsst::ip::diffim::RegularizedKernelSolution< InputT >.

Definition at line 401 of file KernelSolution.cc.

                                             {
        pexLog::TTrace<5>("lsst.ip.diffim.StaticKernelSolution.solve", 
                          "mMat is %d x %d; bVec is %d; cMat is %d x %d; vVec is %d; iVec is %d", 
                          (*_mMat).rows(), (*_mMat).cols(), (*_bVec).size(),
                          (*_cMat).rows(), (*_cMat).cols(), (*_ivVec).size(), (*_iVec).size());

        /* If I put this here I can't check for condition number before solving */
        /*
        _mMat.reset(new Eigen::MatrixXd((*_cMat).transpose() * ((*_ivVec).asDiagonal() * (*_cMat))));
        _bVec.reset(new Eigen::VectorXd((*_cMat).transpose() * ((*_ivVec).asDiagonal() * (*_iVec))));
        */

        if (DEBUG_MATRIX) {
            std::cout << "C" << std::endl;
            std::cout << (*_cMat) << std::endl;
            std::cout << "iV" << std::endl;
            std::cout << (*_ivVec) << std::endl;
            std::cout << "I" << std::endl;
            std::cout << (*_iVec) << std::endl;
        }

        try {
            KernelSolution::solve();
        } catch (pexExcept::Exception &e) {
            LSST_EXCEPT_ADD(e, "Unable to solve static kernel matrix");
            throw e;
        }
        /* Turn matrices into _kernel and _background */
        _setKernel();
    }

Member Data Documentation

template<typename InputT >

double lsst::ip::diffim::StaticKernelSolution< InputT >::_background

protected

Derived differential background estimate.

Definition at line 109 of file KernelSolution.h.

template<typename InputT >

boost::shared_ptr<Eigen::MatrixXd> lsst::ip::diffim::StaticKernelSolution< InputT >::_cMat

protected

K_i x R.

Definition at line 104 of file KernelSolution.h.

template<typename InputT >

boost::shared_ptr<Eigen::VectorXd> lsst::ip::diffim::StaticKernelSolution< InputT >::_iVec

protected

Vectorized I.

Definition at line 105 of file KernelSolution.h.

template<typename InputT >

boost::shared_ptr<Eigen::VectorXd> lsst::ip::diffim::StaticKernelSolution< InputT >::_ivVec

protected

Inverse variance.

Definition at line 106 of file KernelSolution.h.

template<typename InputT >

lsst::afw::math::Kernel::Ptr lsst::ip::diffim::StaticKernelSolution< InputT >::_kernel

protected

Derived single-object convolution kernel.

Definition at line 108 of file KernelSolution.h.

template<typename InputT >

double lsst::ip::diffim::StaticKernelSolution< InputT >::_kSum

protected

Derived kernel sum.

Definition at line 110 of file KernelSolution.h.

---

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

• /home/lsstsw/stack/Linux64/ip_diffim/11.0+1/include/lsst/ip/diffim/KernelSolution.h
• /home/lsstsw/stack/Linux64/ip_diffim/11.0+1/src/KernelSolution.cc