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

#include <KernelSolution.h>

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

Public Types
typedef std::shared_ptr< RegularizedKernelSolution< 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
	RegularizedKernelSolution (lsst::afw::math::KernelList const &basisList, bool fitForBackground, Eigen::MatrixXd const &hMat, lsst::daf::base::PropertySet const &ps)
virtual	~RegularizedKernelSolution ()
void	solve ()
double	getLambda ()
double	estimateRisk (double maxCond)
Eigen::MatrixXd	getM (bool includeHmat=true)
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. More...
void	_setKernelUncertainty ()
	Not implemented. More...

Protected Attributes
Eigen::MatrixXd	_cMat
	K_i x R. More...
Eigen::VectorXd	_iVec
	Vectorized I. More...
Eigen::VectorXd	_ivVec
	Inverse variance. More...
std::shared_ptr< lsst::afw::math::Kernel >	_kernel
	Derived single-object convolution kernel. More...
double	_background
	Derived differential background estimate. More...
double	_kSum
	Derived kernel sum. More...
int	_id
	Unique ID for object. More...
Eigen::MatrixXd	_mMat
	Derived least squares M matrix. More...
Eigen::VectorXd	_bVec
	Derived least squares B vector. More...
Eigen::VectorXd	_aVec
	Derived least squares solution matrix. More...
KernelSolvedBy	_solvedBy
	Type of algorithm used to make solution. More...
bool	_fitForBackground
	Background terms included in fit. More...

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

Detailed Description

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

Definition at line 148 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<RegularizedKernelSolution<InputT> > lsst::ip::diffim::RegularizedKernelSolution< InputT >::Ptr

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

RegularizedKernelSolution()

template<typename InputT >

lsst::ip::diffim::RegularizedKernelSolution< InputT >::RegularizedKernelSolution	(	lsst::afw::math::KernelList const &	basisList,
		bool	fitForBackground,
		Eigen::MatrixXd const &	hMat,
		lsst::daf::base::PropertySet const &	ps
	)

Definition at line 1059 of file KernelSolution.cc.

        :
        StaticKernelSolution<InputT>(basisList, fitForBackground),
        _hMat(hMat),
        _ps(ps.deepCopy())
    {};

~RegularizedKernelSolution()

template<typename InputT >

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

inlinevirtual

Definition at line 157 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 422 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 462 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 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<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 */
        for (kiter = basisList.begin(); kiter != basisList.end(); ++kiter, ++eiter) {
            afwMath::convolve(cimage, templateImage, **kiter, false); /* cimage stores convolved image */
 
            Eigen::MatrixXd cMat = imageToEigenMatrix(cimage).block(startRow,
                                                                    startCol,
                                                                    endRow-startRow,
                                                                    endCol-startCol);
            cMat.resize(cMat.size(), 1);
            *eiter = cMat;
 
        }
 
        double time = t.elapsed();
        LOGL_DEBUG("TRACE3.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<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);
    }

estimateRisk()

template<typename InputT >

double lsst::ip::diffim::RegularizedKernelSolution< InputT >::estimateRisk

(

double

maxCond

)

Definition at line 1072 of file KernelSolution.cc.

                                                                         {
        Eigen::MatrixXd vMat      = this->_cMat.jacobiSvd().matrixV();
        Eigen::MatrixXd vMatvMatT = vMat * vMat.transpose();
 
        /* Find pseudo inverse of mMat, which may be ill conditioned */
        Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eVecValues(this->_mMat);
        Eigen::MatrixXd const& rMat = eVecValues.eigenvectors();
        Eigen::VectorXd eValues = eVecValues.eigenvalues();
        double eMax = eValues.maxCoeff();
        for (int i = 0; i != eValues.rows(); ++i) {
            if (eValues(i) == 0.0) {
                eValues(i) = 0.0;
            }
            else if ((eMax / eValues(i)) > maxCond) {
                LOGL_DEBUG("TRACE3.ip.diffim.RegularizedKernelSolution.estimateRisk",
                           "Truncating eValue %d; %.5e / %.5e = %.5e vs. %.5e",
                           i, eMax, eValues(i), eMax / eValues(i), maxCond);
                eValues(i) = 0.0;
            }
            else {
                eValues(i) = 1.0 / eValues(i);
            }
        }
        Eigen::MatrixXd mInv    = rMat * eValues.asDiagonal() * rMat.transpose();
 
        std::vector<double> lambdas = _createLambdaSteps();
        std::vector<double> risks;
        for (unsigned int i = 0; i < lambdas.size(); i++) {
            double l = lambdas[i];
            Eigen::MatrixXd mLambda = this->_mMat + l * _hMat;
 
            try {
                KernelSolution::solve(mLambda, this->_bVec);
            } catch (pexExcept::Exception &e) {
                LSST_EXCEPT_ADD(e, "Unable to solve regularized kernel matrix");
                throw e;
            }
            Eigen::VectorXd term1 = (this->_aVec.transpose() * vMatvMatT * this->_aVec);
            if (term1.size() != 1)
                throw LSST_EXCEPT(pexExcept::Exception, "Matrix size mismatch");
 
            double term2a = (vMatvMatT * mLambda.inverse()).trace();
 
            Eigen::VectorXd term2b = (this->_aVec.transpose() * (mInv * this->_bVec));
            if (term2b.size() != 1)
                throw LSST_EXCEPT(pexExcept::Exception, "Matrix size mismatch");
 
            double risk   = term1(0) + 2 * (term2a - term2b(0));
            LOGL_DEBUG("TRACE4.ip.diffim.RegularizedKernelSolution.estimateRisk",
                       "Lambda = %.3f, Risk = %.5e",
                       l, risk);
            LOGL_DEBUG("TRACE5.ip.diffim.RegularizedKernelSolution.estimateRisk",
                       "%.5e + 2 * (%.5e - %.5e)",
                       term1(0), term2a, term2b(0));
            risks.push_back(risk);
        }
        std::vector<double>::iterator it = min_element(risks.begin(), risks.end());
        int index = distance(risks.begin(), it);
        LOGL_DEBUG("TRACE3.ip.diffim.RegularizedKernelSolution.estimateRisk",
                   "Minimum Risk = %.3e at lambda = %.3e", risks[index], lambdas[index]);
 
        return lambdas[index];
 
    }

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

getLambda()

template<typename InputT >

double lsst::ip::diffim::RegularizedKernelSolution< InputT >::getLambda ( )

inline

Definition at line 159 of file KernelSolution.h.

{return _lambda;}

getM() [1/2]

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

inlineinherited

Definition at line 64 of file KernelSolution.h.

{return _mMat;}

getM() [2/2]

template<typename InputT >

Eigen::MatrixXd lsst::ip::diffim::RegularizedKernelSolution< InputT >::getM

(

bool

includeHmat = true

)

Definition at line 1138 of file KernelSolution.cc.

                                                                          {
        if (includeHmat == true) {
            return this->_mMat + _lambda * _hMat;
        }
        else {
            return this->_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::RegularizedKernelSolution< InputT >::solve

virtual

Reimplemented from lsst::ip::diffim::StaticKernelSolution< InputT >.

Definition at line 1148 of file KernelSolution.cc.

                                                  {
 
        LOGL_DEBUG("TRACE3.ip.diffim.RegularizedKernelSolution.solve",
                   "cMat is %d x %d; vVec is %d; iVec is %d; hMat is %d x %d",
                   this->_cMat.rows(), this->_cMat.cols(), this->_ivVec.size(),
                   this->_iVec.size(), _hMat.rows(), _hMat.cols());
 
        if (DEBUG_MATRIX2) {
            std::cout << "ID: " << (this->_id) << std::endl;
            std::cout << "C:" << std::endl;
            std::cout << this->_cMat << std::endl;
            std::cout << "Sigma^{-1}:" << std::endl;
            std::cout << Eigen::MatrixXd(this->_ivVec.asDiagonal()) << std::endl;
            std::cout << "Y:" << std::endl;
            std::cout << this->_iVec << std::endl;
            std::cout << "H:" << std::endl;
            std::cout << _hMat << std::endl;
        }
 
 
        this->_mMat = this->_cMat.transpose() * this->_ivVec.asDiagonal() * this->_cMat;
        this->_bVec = this->_cMat.transpose() * this->_ivVec.asDiagonal() * this->_iVec;
 
 
        /* See N.R. 18.5
 
           Matrix equation to solve is Y = C a + N
           Y   = vectorized version of I (I = image to not convolve)
           C_i = K_i (x) R (R = image to convolve)
           a   = kernel coefficients
           N   = noise
 
           If we reweight everything by the inverse square root of the noise
           covariance, we get a linear model with the identity matrix for
           the noise.  The problem can then be solved using least squares,
           with the normal equations
 
              C^T Y = C^T C a
 
           or
 
              b = M a
 
           with
 
              b = C^T Y
              M = C^T C
              a = (C^T C)^{-1} C^T Y
 
 
           We can regularize the least square problem
 
              Y = C a + lambda a^T H a       (+ N, which can be ignored)
 
           or the normal equations
 
              (C^T C + lambda H) a = C^T Y
 
 
           The solution to the regularization of the least squares problem is
 
              a = (C^T C + lambda H)^{-1} C^T Y
 
           The approximation to Y is
 
              C a = C (C^T C + lambda H)^{-1} C^T Y
 
           with smoothing matrix
 
              S = C (C^T C + lambda H)^{-1} C^T
 
        */
 
        std::string lambdaType = _ps->getAsString("lambdaType");
        if (lambdaType == "absolute") {
            _lambda = _ps->getAsDouble("lambdaValue");
        }
        else if (lambdaType ==  "relative") {
            _lambda  = this->_mMat.trace() / this->_hMat.trace();
            _lambda *= _ps->getAsDouble("lambdaScaling");
        }
        else if (lambdaType ==  "minimizeBiasedRisk") {
            double tol = _ps->getAsDouble("maxConditionNumber");
            _lambda = estimateRisk(tol);
        }
        else if (lambdaType ==  "minimizeUnbiasedRisk") {
            _lambda = estimateRisk(std::numeric_limits<double>::max());
        }
        else {
            throw LSST_EXCEPT(pexExcept::Exception, "lambdaType in PropertySet not recognized");
        }
 
        LOGL_DEBUG("TRACE3.ip.diffim.RegularizedKernelSolution.solve",
                   "Applying kernel regularization with lambda = %.2e", _lambda);
 
 
        try {
            KernelSolution::solve(this->_mMat + _lambda * _hMat, this->_bVec);
        } catch (pexExcept::Exception &e) {
            LSST_EXCEPT_ADD(e, "Unable to solve static kernel matrix");
            throw e;
        }
        /* Turn matrices into _kernel and _background */
        StaticKernelSolution<InputT>::_setKernel();
    }

solve() [2/2]

void lsst::ip::diffim::KernelSolution::solve ( Eigen::MatrixXd const &  mMat, Eigen::VectorXd const &  bVec  )

virtualinherited

Definition at line 131 of file KernelSolution.cc.

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

Member Data Documentation

_aVec

Eigen::VectorXd lsst::ip::diffim::KernelSolution::_aVec

protectedinherited

Derived least squares solution matrix.

Definition at line 75 of file KernelSolution.h.

_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/1a1d771/Linux64/ip_diffim/20.0.0/include/lsst/ip/diffim/KernelSolution.h
• /j/snowflake/release/lsstsw/stack/1a1d771/Linux64/ip_diffim/20.0.0/src/KernelSolution.cc