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

#include <KernelSolution.h>

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

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

Public Member Functions
	MaskedKernelSolution (lsst::afw::math::KernelList const &basisList, bool fitForBackground)
virtual	~MaskedKernelSolution ()
virtual void	buildOrig (lsst::afw::image::Image< InputT > const &templateImage, lsst::afw::image::Image< InputT > const &scienceImage, lsst::afw::image::Image< lsst::afw::image::VariancePixel > const &varianceEstimate, lsst::afw::image::Mask< lsst::afw::image::MaskPixel > pixelMask)
virtual void	buildWithMask (lsst::afw::image::Image< InputT > const &templateImage, lsst::afw::image::Image< InputT > const &scienceImage, lsst::afw::image::Image< lsst::afw::image::VariancePixel > const &varianceEstimate, lsst::afw::image::Mask< lsst::afw::image::MaskPixel > const &pixelMask)
virtual void	buildSingleMaskOrig (lsst::afw::image::Image< InputT > const &templateImage, lsst::afw::image::Image< InputT > const &scienceImage, lsst::afw::image::Image< lsst::afw::image::VariancePixel > const &varianceEstimate, lsst::geom::Box2I maskBox)
void	solve ()
virtual void	solve (Eigen::MatrixXd const &mMat, Eigen::VectorXd const &bVec)
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::Kernel >	getKernel ()
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 ()
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::MaskedKernelSolution< InputT >

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

template<typename InputT >

typedef std::shared_ptr<MaskedKernelSolution<InputT> > lsst::ip::diffim::MaskedKernelSolution< InputT >::Ptr

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

MaskedKernelSolution()

template<typename InputT >

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

Definition at line 493 of file KernelSolution.cc.

          :
        StaticKernelSolution<InputT>(basisList, fitForBackground)
    {};

~MaskedKernelSolution()

template<typename InputT >

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

inlinevirtual

Definition at line 125 of file KernelSolution.h.

{};

Member Function Documentation

_setKernel()

template<typename InputT >

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

protectedinherited

Set kernel after solution.

Definition at line 423 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     =
            std::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);
 
        std::shared_ptr<ImageT> 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);
        }
    }

_setKernelUncertainty()

template<typename InputT >

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

protectedinherited

Not implemented.

Definition at line 463 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();
         */
    }

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 )

virtualinherited

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 =
            std::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<std::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.
        */
        geom::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::cpu_timer t;
 
        /* 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<Eigen::MatrixXd> convolvedEigenList(nKernelParameters);
 
        /* Iterators over convolved image list and basis list */
        typename std::vector<Eigen::MatrixXd>::iterator eiter = convolvedEigenList.begin();
        /* Create C_i in the formalism of Alard & Lupton */
        afwMath::ConvolutionControl convolutionControl;
        convolutionControl.setDoNormalize(false);
        for (kiter = basisList.begin(); kiter != basisList.end(); ++kiter, ++eiter) {
            afwMath::convolve(cimage, templateImage, **kiter, convolutionControl); /* cimage stores convolved image */
 
            Eigen::MatrixXd cMat = imageToEigenMatrix(cimage).block(startRow,
                                                                    startCol,
                                                                    endRow-startRow,
                                                                    endCol-startCol);
            cMat.resize(cMat.size(), 1);
            *eiter = cMat;
 
        }
 
        t.stop();
        double time = 1e-9 * t.elapsed().wall;
        LOGL_DEBUG("TRACE3.ip.diffim.StaticKernelSolution.build",
                   "Total compute time to do basis convolutions : %.2f s", time);
 
        /*
           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<Eigen::MatrixXd>::iterator eiterj = convolvedEigenList.begin();
        typename std::vector<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 = cMat;
        _ivVec = eigeniVariance.col(0);
        _iVec = eigenScience.col(0);
 
        /* Make these outside of solve() so I can check condition number */
        _mMat = _cMat.transpose() * (_ivVec.asDiagonal() * _cMat);
        _bVec = _cMat.transpose() * (_ivVec.asDiagonal() * _iVec);
    }

buildOrig()

template<typename InputT >

void lsst::ip::diffim::MaskedKernelSolution< InputT >::buildOrig ( lsst::afw::image::Image< InputT > const & templateImage, lsst::afw::image::Image< InputT > const & scienceImage, lsst::afw::image::Image< lsst::afw::image::VariancePixel > const & varianceEstimate, lsst::afw::image::Mask< lsst::afw::image::MaskPixel > pixelMask )

virtual

Definition at line 667 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 =
            std::dynamic_pointer_cast<afwMath::LinearCombinationKernel>(this->_kernel)->getKernelList();
        std::vector<std::shared_ptr<afwMath::Kernel> >::const_iterator kiter = basisList.begin();
 
        /* Only BAD pixels marked in this mask */
        lsst::afw::image::Mask<lsst::afw::image::MaskPixel> sMask(pixelMask, true);
        afwImage::MaskPixel bitMask =
            (afwImage::Mask<afwImage::MaskPixel>::getPlaneBitMask("BAD") |
             afwImage::Mask<afwImage::MaskPixel>::getPlaneBitMask("SAT") |
             afwImage::Mask<afwImage::MaskPixel>::getPlaneBitMask("EDGE"));
        sMask &= bitMask;
        /* TBD: Need to figure out a way to spread this; currently done in Python */
 
        unsigned int const nKernelParameters     = basisList.size();
        unsigned int const nBackgroundParameters = this->_fitForBackground ? 1 : 0;
        unsigned int const nParameters           = nKernelParameters + nBackgroundParameters;
 
        /* 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.
        */
        /* Ignore known EDGE pixels for speed */
        geom::Box2I shrunkLocalBBox = (*kiter)->shrinkBBox(templateImage.getBBox(afwImage::LOCAL));
        LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
                   "Limits of good pixels after convolution: %d,%d -> %d,%d (local)",
                   shrunkLocalBBox.getMinX(), shrunkLocalBBox.getMinY(),
                   shrunkLocalBBox.getMaxX(), shrunkLocalBBox.getMaxY());
 
        /* Subimages are addressed in raw pixel coordinates */
        unsigned int startCol = shrunkLocalBBox.getMinX();
        unsigned int startRow = shrunkLocalBBox.getMinY();
        unsigned int endCol   = shrunkLocalBBox.getMaxX();
        unsigned int endRow   = shrunkLocalBBox.getMaxY();
        /* Needed for Eigen block slicing operation */
        endCol += 1;
        endRow += 1;
 
        boost::timer::cpu_timer t;
 
        /* Eigen representation of input images; only the pixels that are unconvolved in cimage below */
        Eigen::MatrixXi eMask = maskToEigenMatrix(sMask).block(startRow,
                                                               startCol,
                                                               endRow-startRow,
                                                               endCol-startCol);
 
        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 */
        eMask.resize(eMask.rows()*eMask.cols(), 1);
        eigenTemplate.resize(eigenTemplate.rows()*eigenTemplate.cols(), 1);
        eigenScience.resize(eigenScience.rows()*eigenScience.cols(), 1);
        eigeniVariance.resize(eigeniVariance.rows()*eigeniVariance.cols(), 1);
 
        /* Do the masking, slowly for now... */
        Eigen::MatrixXd maskedEigenTemplate(eigenTemplate.rows(), 1);
        Eigen::MatrixXd maskedEigenScience(eigenScience.rows(), 1);
        Eigen::MatrixXd maskedEigeniVariance(eigeniVariance.rows(), 1);
        int nGood = 0;
        for (int i = 0; i < eMask.rows(); i++) {
            if (eMask(i, 0) == 0) {
                maskedEigenTemplate(nGood, 0) = eigenTemplate(i, 0);
                maskedEigenScience(nGood, 0) = eigenScience(i, 0);
                maskedEigeniVariance(nGood, 0) = eigeniVariance(i, 0);
                nGood += 1;
            }
        }
        /* Can't resize() since all values are lost; use blocks */
        eigenTemplate    = maskedEigenTemplate.block(0, 0, nGood, 1);
        eigenScience     = maskedEigenScience.block(0, 0, nGood, 1);
        eigeniVariance   = maskedEigeniVariance.block(0, 0, nGood, 1);
 
 
        /* Holds image convolved with basis function */
        afwImage::Image<InputT> cimage(templateImage.getDimensions());
 
        /* Holds eigen representation of image convolved with all basis functions */
        std::vector<Eigen::MatrixXd> convolvedEigenList(nKernelParameters);
 
        /* Iterators over convolved image list and basis list */
        typename std::vector<Eigen::MatrixXd>::iterator eiter = convolvedEigenList.begin();
        /* Create C_i in the formalism of Alard & Lupton */
        afwMath::ConvolutionControl convolutionControl;
        convolutionControl.setDoNormalize(false);
       for (kiter = basisList.begin(); kiter != basisList.end(); ++kiter, ++eiter) {
            afwMath::convolve(cimage, templateImage, **kiter, convolutionControl); /* cimage stores convolved image */
 
            Eigen::MatrixXd cMat = imageToEigenMatrix(cimage).block(startRow,
                                                                    startCol,
                                                                    endRow-startRow,
                                                                    endCol-startCol);
            cMat.resize(cMat.size(), 1);
 
            /* Do masking */
            Eigen::MatrixXd maskedcMat(cMat.rows(), 1);
            int nGood = 0;
            for (int i = 0; i < eMask.rows(); i++) {
                if (eMask(i, 0) == 0) {
                    maskedcMat(nGood, 0) = cMat(i, 0);
                    nGood += 1;
                }
            }
            cMat = maskedcMat.block(0, 0, nGood, 1);
            *eiter = cMat;
        }
 
        t.stop();
        double time = 1e-9 * t.elapsed().wall;
        LOGL_DEBUG("TRACE3.ip.diffim.StaticKernelSolution.build",
                   "Total compute time to do basis convolutions : %.2f s", time);
 
        /*
           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<Eigen::MatrixXd>::iterator eiterj = convolvedEigenList.begin();
        typename std::vector<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 (this->_fitForBackground)
            cMat.col(nParameters-1).fill(1.);
 
        this->_cMat = cMat;
        this->_ivVec = eigeniVariance.col(0);
        this->_iVec = eigenScience.col(0);
 
        /* Make these outside of solve() so I can check condition number */
        this->_mMat = this->_cMat.transpose() * this->_ivVec.asDiagonal() * this->_cMat;
        this->_bVec = this->_cMat.transpose() * this->_ivVec.asDiagonal() * this->_iVec;
 
    }

buildSingleMaskOrig()

template<typename InputT >

void lsst::ip::diffim::MaskedKernelSolution< InputT >::buildSingleMaskOrig ( lsst::afw::image::Image< InputT > const & templateImage, lsst::afw::image::Image< InputT > const & scienceImage, lsst::afw::image::Image< lsst::afw::image::VariancePixel > const & varianceEstimate, lsst::geom::Box2I maskBox )

virtual

Definition at line 831 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 =
            std::dynamic_pointer_cast<afwMath::LinearCombinationKernel>(this->_kernel)->getKernelList();
 
        unsigned int const nKernelParameters     = basisList.size();
        unsigned int const nBackgroundParameters = this->_fitForBackground ? 1 : 0;
        unsigned int const nParameters           = nKernelParameters + nBackgroundParameters;
 
        std::vector<std::shared_ptr<afwMath::Kernel> >::const_iterator kiter = basisList.begin();
 
        /*
           NOTE : If we are using these views in Afw's Image space, we need to
           make sure and compensate for the XY0 of the image:
 
           geom::Box2I fullBBox = templateImage.getBBox();
           int maskStartCol = maskBox.getMinX();
           int maskStartRow = maskBox.getMinY();
           int maskEndCol   = maskBox.getMaxX();
           int maskEndRow   = maskBox.getMaxY();
 
 
          If we are going to be doing the slicing in Eigen matrices derived from
          the images, ignore the XY0.
 
           geom::Box2I fullBBox = templateImage.getBBox(afwImage::LOCAL);
 
           int maskStartCol = maskBox.getMinX() - templateImage.getX0();
           int maskStartRow = maskBox.getMinY() - templateImage.getY0();
           int maskEndCol   = maskBox.getMaxX() - templateImage.getX0();
           int maskEndRow   = maskBox.getMaxY() - templateImage.getY0();
 
        */
 
 
        geom::Box2I shrunkBBox = (*kiter)->shrinkBBox(templateImage.getBBox());
 
        LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
                   "Limits of good pixels after convolution: %d,%d -> %d,%d",
                   shrunkBBox.getMinX(), shrunkBBox.getMinY(),
                   shrunkBBox.getMaxX(), shrunkBBox.getMaxY());
 
        unsigned int const startCol = shrunkBBox.getMinX();
        unsigned int const startRow = shrunkBBox.getMinY();
        unsigned int const endCol   = shrunkBBox.getMaxX();
        unsigned int const endRow   = shrunkBBox.getMaxY();
 
        /* NOTE: no endCol/endRow += 1 for block slicing, since we are doing the
           slicing using Afw, not Eigen
 
           Eigen arrays have index 0,0 in the upper right, while LSST images
           have 0,0 in the lower left.  The y-axis is flipped.  When translating
           Images to Eigen matrices in ipDiffim::imageToEigenMatrix this is
           accounted for.  However, we still need to be aware of this fact if
           addressing subregions of an Eigen matrix.  This is why the slicing is
           done in Afw, its cleaner.
        */
 
 
        /* Inner limits; of mask box */
        int maskStartCol = maskBox.getMinX();
        int maskStartRow = maskBox.getMinY();
        int maskEndCol   = maskBox.getMaxX();
        int maskEndRow   = maskBox.getMaxY();
 
        /*
 
        |---------------------------|
        |      Kernel Boundary      |
        |  |---------------------|  |
        |  |         Top         |  |
        |  |......_________......|  |
        |  |      |       |      |  |
        |  |  L   |  Box  |  R   |  |
        |  |      |       |      |  |
        |  |......---------......|  |
        |  |        Bottom       |  |
        |  |---------------------|  |
        |                           |
        |---------------------------|
 
        4 regions we want to extract from the pixels: top bottom left right
 
        */
        geom::Box2I tBox = geom::Box2I(geom::Point2I(startCol, maskEndRow + 1),
                                       geom::Point2I(endCol, endRow));
 
        geom::Box2I bBox = geom::Box2I(geom::Point2I(startCol, startRow),
                                       geom::Point2I(endCol, maskStartRow - 1));
 
        geom::Box2I lBox = geom::Box2I(geom::Point2I(startCol, maskStartRow),
                                       geom::Point2I(maskStartCol - 1, maskEndRow));
 
        geom::Box2I rBox = geom::Box2I(geom::Point2I(maskEndCol + 1, maskStartRow),
                                       geom::Point2I(endCol, maskEndRow));
 
        LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
                   "Upper good pixel region: %d,%d -> %d,%d",
                   tBox.getMinX(), tBox.getMinY(), tBox.getMaxX(), tBox.getMaxY());
        LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
                   "Bottom good pixel region: %d,%d -> %d,%d",
                   bBox.getMinX(), bBox.getMinY(), bBox.getMaxX(), bBox.getMaxY());
        LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
                   "Left good pixel region: %d,%d -> %d,%d",
                   lBox.getMinX(), lBox.getMinY(), lBox.getMaxX(), lBox.getMaxY());
        LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
                   "Right good pixel region: %d,%d -> %d,%d",
                   rBox.getMinX(), rBox.getMinY(), rBox.getMaxX(), rBox.getMaxY());
 
        std::vector<geom::Box2I> boxArray;
        boxArray.push_back(tBox);
        boxArray.push_back(bBox);
        boxArray.push_back(lBox);
        boxArray.push_back(rBox);
 
        int totalSize = tBox.getWidth() * tBox.getHeight();
        totalSize    += bBox.getWidth() * bBox.getHeight();
        totalSize    += lBox.getWidth() * lBox.getHeight();
        totalSize    += rBox.getWidth() * rBox.getHeight();
 
        Eigen::MatrixXd eigenTemplate(totalSize, 1);
        Eigen::MatrixXd eigenScience(totalSize, 1);
        Eigen::MatrixXd eigeniVariance(totalSize, 1);
        eigenTemplate.setZero();
        eigenScience.setZero();
        eigeniVariance.setZero();
 
        boost::timer::cpu_timer t;
 
        int nTerms = 0;
        typename std::vector<geom::Box2I>::iterator biter = boxArray.begin();
        for (; biter != boxArray.end(); ++biter) {
            int area = (*biter).getWidth() * (*biter).getHeight();
 
            afwImage::Image<InputT> siTemplate(templateImage, *biter);
            afwImage::Image<InputT> siScience(scienceImage, *biter);
            afwImage::Image<InputT> sVarEstimate(varianceEstimate, *biter);
 
            Eigen::MatrixXd eTemplate = imageToEigenMatrix(siTemplate);
            Eigen::MatrixXd eScience = imageToEigenMatrix(siScience);
            Eigen::MatrixXd eiVarEstimate = imageToEigenMatrix(sVarEstimate).array().inverse().matrix();
 
            eTemplate.resize(area, 1);
            eScience.resize(area, 1);
            eiVarEstimate.resize(area, 1);
 
            eigenTemplate.block(nTerms, 0, area, 1) = eTemplate.block(0, 0, area, 1);
            eigenScience.block(nTerms, 0, area, 1) = eScience.block(0, 0, area, 1);
            eigeniVariance.block(nTerms, 0, area, 1) = eiVarEstimate.block(0, 0, area, 1);
 
            nTerms += area;
        }
 
        afwImage::Image<InputT> cimage(templateImage.getDimensions());
 
        std::vector<Eigen::MatrixXd> convolvedEigenList(nKernelParameters);
        typename std::vector<Eigen::MatrixXd>::iterator eiter = convolvedEigenList.begin();
        /* Create C_i in the formalism of Alard & Lupton */
        afwMath::ConvolutionControl convolutionControl;
        convolutionControl.setDoNormalize(false);
        for (kiter = basisList.begin(); kiter != basisList.end(); ++kiter, ++eiter) {
            afwMath::convolve(cimage, templateImage, **kiter, convolutionControl); /* cimage stores convolved image */
            Eigen::MatrixXd cMat(totalSize, 1);
            cMat.setZero();
 
            int nTerms = 0;
            typename std::vector<geom::Box2I>::iterator biter = boxArray.begin();
            for (; biter != boxArray.end(); ++biter) {
                int area = (*biter).getWidth() * (*biter).getHeight();
 
                afwImage::Image<InputT> csubimage(cimage, *biter);
                Eigen::MatrixXd esubimage = imageToEigenMatrix(csubimage);
                esubimage.resize(area, 1);
                cMat.block(nTerms, 0, area, 1) = esubimage.block(0, 0, area, 1);
 
                nTerms += area;
            }
 
            *eiter = cMat;
 
        }
 
        t.stop();
        double time = 1e-9 * t.elapsed().wall;
        LOGL_DEBUG("TRACE3.ip.diffim.MaskedKernelSolution.build",
                   "Total compute time to do basis convolutions : %.2f s", time);
 
        /*
           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<Eigen::MatrixXd>::iterator eiterj = convolvedEigenList.begin();
        typename std::vector<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 (this->_fitForBackground)
            cMat.col(nParameters-1).fill(1.);
 
        this->_cMat = cMat;
        this->_ivVec = eigeniVariance.col(0);
        this->_iVec = eigenScience.col(0);
 
        /* Make these outside of solve() so I can check condition number */
        this->_mMat = this->_cMat.transpose() * this->_ivVec.asDiagonal() * this->_cMat;
        this->_bVec = this->_cMat.transpose() * this->_ivVec.asDiagonal() * this->_iVec;
    }

buildWithMask()

template<typename InputT >

void lsst::ip::diffim::MaskedKernelSolution< InputT >::buildWithMask ( lsst::afw::image::Image< InputT > const & templateImage, lsst::afw::image::Image< InputT > const & scienceImage, lsst::afw::image::Image< lsst::afw::image::VariancePixel > const & varianceEstimate, lsst::afw::image::Mask< lsst::afw::image::MaskPixel > const & pixelMask )

virtual

Definition at line 501 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");
        }
 
        /* Full footprint of all input images */
        std::shared_ptr<afwDet::Footprint> fullFp(
            new afwDet::Footprint(std::make_shared<afwGeom::SpanSet>(templateImage.getBBox())));
 
        afwMath::KernelList basisList =
            std::dynamic_pointer_cast<afwMath::LinearCombinationKernel>(this->_kernel)->getKernelList();
        std::vector<std::shared_ptr<afwMath::Kernel> >::const_iterator kiter = basisList.begin();
 
        /* Only BAD pixels marked in this mask */
        afwImage::MaskPixel bitMask =
            (afwImage::Mask<afwImage::MaskPixel>::getPlaneBitMask("BAD") |
             afwImage::Mask<afwImage::MaskPixel>::getPlaneBitMask("SAT") |
             afwImage::Mask<afwImage::MaskPixel>::getPlaneBitMask("NO_DATA") |
             afwImage::Mask<afwImage::MaskPixel>::getPlaneBitMask("EDGE"));
 
        /* Create a Footprint that contains all the masked pixels set above */
        afwDet::Threshold threshold = afwDet::Threshold(bitMask, afwDet::Threshold::BITMASK, true);
        afwDet::FootprintSet maskFpSet(pixelMask, threshold, true);
 
        /* And spread it by the kernel half width */
        int growPix = (*kiter)->getCtr().getX();
        afwDet::FootprintSet maskedFpSetGrown(maskFpSet, growPix, true);
 
#if 0
        for (typename afwDet::FootprintSet::FootprintList::iterator
                 ptr = maskedFpSetGrown.getFootprints()->begin(),
                 end = maskedFpSetGrown.getFootprints()->end();
             ptr != end;
             ++ptr) {
 
            afwDet::setMaskFromFootprint(finalMask,
                                         (**ptr),
                                         afwImage::Mask<afwImage::MaskPixel>::getPlaneBitMask("BAD"));
        }
#endif
 
        afwImage::Mask<afwImage::MaskPixel> finalMask(pixelMask.getDimensions());
        for (auto const & foot : *(maskedFpSetGrown.getFootprints())) {
            foot->getSpans()->setMask(finalMask, afwImage::Mask<afwImage::MaskPixel>::getPlaneBitMask("BAD"));
        }
        pixelMask.writeFits("pixelmask.fits");
        finalMask.writeFits("finalmask.fits");
 
 
        ndarray::Array<int, 1, 1> maskArray =
            ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
        fullFp->getSpans()->flatten(maskArray, finalMask.getArray(), templateImage.getXY0());
        auto maskEigen = ndarray::asEigenMatrix(maskArray);
 
        ndarray::Array<InputT, 1, 1> arrayTemplate =
            ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
        fullFp->getSpans()->flatten(arrayTemplate, templateImage.getArray(), templateImage.getXY0());
        auto eigenTemplate0 = ndarray::asEigenMatrix(arrayTemplate);
 
        ndarray::Array<InputT, 1, 1> arrayScience =
            ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
        fullFp->getSpans()->flatten(arrayScience, scienceImage.getArray(), scienceImage.getXY0());
        auto eigenScience0 = ndarray::asEigenMatrix(arrayScience);
 
        ndarray::Array<afwImage::VariancePixel, 1, 1> arrayVariance =
            ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
        fullFp->getSpans()->flatten(arrayVariance, varianceEstimate.getArray(), varianceEstimate.getXY0());
        auto eigenVariance0 = ndarray::asEigenMatrix(arrayVariance);
 
        int nGood = 0;
        for (int i = 0; i < maskEigen.size(); i++) {
            if (maskEigen(i) == 0.0)
                nGood += 1;
        }
 
        Eigen::VectorXd eigenTemplate(nGood);
        Eigen::VectorXd eigenScience(nGood);
        Eigen::VectorXd eigenVariance(nGood);
        int nUsed = 0;
        for (int i = 0; i < maskEigen.size(); i++) {
            if (maskEigen(i) == 0.0) {
                eigenTemplate(nUsed) = eigenTemplate0(i);
                eigenScience(nUsed) = eigenScience0(i);
                eigenVariance(nUsed) = eigenVariance0(i);
                nUsed += 1;
            }
        }
 
 
        boost::timer::cpu_timer t;
 
        unsigned int const nKernelParameters     = basisList.size();
        unsigned int const nBackgroundParameters = this->_fitForBackground ? 1 : 0;
        unsigned int const nParameters           = nKernelParameters + nBackgroundParameters;
 
        /* Holds image convolved with basis function */
        afwImage::Image<InputT> cimage(templateImage.getDimensions());
 
        /* Holds eigen representation of image convolved with all basis functions */
        std::vector<Eigen::VectorXd> convolvedEigenList(nKernelParameters);
 
        /* Iterators over convolved image list and basis list */
        typename std::vector<Eigen::VectorXd>::iterator eiter =  convolvedEigenList.begin();
 
        /* Create C_i in the formalism of Alard & Lupton */
        afwMath::ConvolutionControl convolutionControl;
        convolutionControl.setDoNormalize(false);
        for (kiter = basisList.begin(); kiter != basisList.end(); ++kiter, ++eiter) {
            afwMath::convolve(cimage, templateImage, **kiter, convolutionControl); /* cimage stores convolved image */
 
            ndarray::Array<InputT, 1, 1> arrayC =
                ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
            fullFp->getSpans()->flatten(arrayC, cimage.getArray(), cimage.getXY0());
            auto eigenC0 = ndarray::asEigenMatrix(arrayC);
 
            Eigen::VectorXd eigenC(nGood);
            int nUsed = 0;
            for (int i = 0; i < maskEigen.size(); i++) {
                if (maskEigen(i) == 0.0) {
                    eigenC(nUsed) = eigenC0(i);
                    nUsed += 1;
                }
            }
 
            *eiter = eigenC;
        }
        t.stop();
        double time = 1e-9 * t.elapsed().wall;
        LOGL_DEBUG("TRACE3.ip.diffim.StaticKernelSolution.buildWithMask",
                   "Total compute time to do basis convolutions : %.2f s", time);
 
        /* Load matrix with all convolved images */
        Eigen::MatrixXd cMat(eigenTemplate.size(), nParameters);
        typename std::vector<Eigen::VectorXd>::iterator eiterj = convolvedEigenList.begin();
        typename std::vector<Eigen::VectorXd>::iterator eiterE = convolvedEigenList.end();
        for (unsigned int kidxj = 0; eiterj != eiterE; eiterj++, kidxj++) {
            cMat.block(0, kidxj, eigenTemplate.size(), 1) =
                Eigen::MatrixXd(*eiterj).block(0, 0, eigenTemplate.size(), 1);
        }
        /* Treat the last "image" as all 1's to do the background calculation. */
        if (this->_fitForBackground)
            cMat.col(nParameters-1).fill(1.);
 
        this->_cMat = cMat;
        this->_ivVec = eigenVariance.array().inverse().matrix();
        this->_iVec = eigenScience;
 
        /* Make these outside of solve() so I can check condition number */
        this->_mMat = this->_cMat.transpose() * this->_ivVec.asDiagonal() * (this->_cMat);
        this->_bVec = this->_cMat.transpose() * this->_ivVec.asDiagonal() * (this->_iVec);
    }

getB()

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

inlineinherited

Definition at line 65 of file KernelSolution.h.

{return _bVec;}

getBackground()

template<typename InputT >

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

virtualinherited

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

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.

{ return _id; }

getKernel()

template<typename InputT >

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

virtualinherited

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

getKsum()

template<typename InputT >

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

virtualinherited

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

getM()

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

inlineinherited

Definition at line 64 of file KernelSolution.h.

{return _mMat;}

getSolutionPair()

template<typename InputT >

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

virtualinherited

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

getSolvedBy()

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

inlineinherited

Definition at line 60 of file KernelSolution.h.

{return _solvedBy;}

makeKernelImage()

template<typename InputT >

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

virtualinherited

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

printA()

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

inlineinherited

Definition at line 68 of file KernelSolution.h.

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

printB()

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

inlineinherited

Definition at line 67 of file KernelSolution.h.

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

printM()

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

inlineinherited

Definition at line 66 of file KernelSolution.h.

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

solve() [1/2]

template<typename InputT >

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

virtualinherited

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

Definition at line 397 of file KernelSolution.cc.

                                             {
        LOGL_DEBUG("TRACE3.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 (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();
    }

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

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

protectedinherited

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

protectedinherited

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

protectedinherited

Vectorized I.

Definition at line 106 of file KernelSolution.h.

_ivVec

template<typename InputT >

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

protectedinherited

Inverse variance.

Definition at line 107 of file KernelSolution.h.

_kernel

template<typename InputT >

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

protectedinherited

Derived single-object convolution kernel.

Definition at line 109 of file KernelSolution.h.

_kSum

template<typename InputT >

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

protectedinherited

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:

• /j/snowflake/release/lsstsw/stack/lsst-scipipe-9.0.0/Linux64/ip_diffim/g1a2382251a+b4475c5878/include/lsst/ip/diffim/KernelSolution.h
• /j/snowflake/release/lsstsw/stack/lsst-scipipe-9.0.0/Linux64/ip_diffim/g1a2382251a+b4475c5878/src/KernelSolution.cc