KernelSolution.h

Go to the documentation of this file.

// -*- lsst-c++ -*-
#ifndef LSST_IP_DIFFIM_KERNELSOLUTION_H
#define LSST_IP_DIFFIM_KERNELSOLUTION_H

#include <memory>
#include "Eigen/Core"

#include "lsst/afw/math.h"
#include "lsst/afw/geom.h"
#include "lsst/afw/image.h"
#include "lsst/pex/policy.h"

namespace lsst { 
namespace ip { 
namespace diffim {

    /* 
     * @brief Method used to solve for M and B
     */

    class KernelSolution {
    public:
        typedef std::shared_ptr<KernelSolution> Ptr;
        typedef lsst::afw::math::Kernel::Pixel PixelT;
        typedef lsst::afw::image::Image<lsst::afw::math::Kernel::Pixel> ImageT;

        enum KernelSolvedBy {
            NONE          = 0,
            CHOLESKY_LDLT = 1,
            CHOLESKY_LLT  = 2,
            LU            = 3,
            EIGENVECTOR   = 4
        };

        enum ConditionNumberType {
            EIGENVALUE = 0,
            SVD        = 1
        };

        explicit KernelSolution(Eigen::MatrixXd mMat,
                                Eigen::VectorXd bVec,
                                bool fitForBackground);
        explicit KernelSolution(bool fitForBackground);
        explicit KernelSolution();

        virtual ~KernelSolution() {};
        virtual void solve();
        virtual void solve(Eigen::MatrixXd const& mMat, 
                           Eigen::VectorXd const& bVec);
        KernelSolvedBy getSolvedBy() {return _solvedBy;} 
        virtual double getConditionNumber(ConditionNumberType conditionType);
        virtual double getConditionNumber(Eigen::MatrixXd const& mMat, ConditionNumberType conditionType);

        inline Eigen::MatrixXd const& getM() {return _mMat;}
        inline Eigen::VectorXd const& getB() {return _bVec;}
        void printM() {std::cout << _mMat << std::endl;}
        void printB() {std::cout << _bVec << std::endl;}
        void printA() {std::cout << _aVec << std::endl;}
        inline int getId() const { return _id; }

    protected:
        int _id;                                                
        Eigen::MatrixXd _mMat;               
        Eigen::VectorXd _bVec;               
        Eigen::VectorXd _aVec;               
        KernelSolvedBy _solvedBy;                               
        bool _fitForBackground;                                 
        static int _SolutionId;                                 

    };

    template <typename InputT>
    class StaticKernelSolution : public KernelSolution {
    public:
        typedef std::shared_ptr<StaticKernelSolution<InputT> > Ptr;

        StaticKernelSolution(lsst::afw::math::KernelList const& basisList,
                             bool fitForBackground);
        virtual ~StaticKernelSolution() {};

        /* Overrides KernelSolution */
        void solve();

        /* Used by RegularizedKernelSolution */
        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();

    protected:
        Eigen::MatrixXd _cMat;               
        Eigen::VectorXd _iVec;               
        Eigen::VectorXd _ivVec;              

        std::shared_ptr<lsst::afw::math::Kernel> _kernel;                   
        double _background;                                     
        double _kSum;                                           

        void _setKernel();                                      
        void _setKernelUncertainty();                           
    };


    template <typename InputT>
    class MaskedKernelSolution : public StaticKernelSolution<InputT> {
    public:
        typedef std::shared_ptr<MaskedKernelSolution<InputT> > Ptr;

        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::afw::geom::Box2I maskBox);
    };



    template <typename InputT>
    class RegularizedKernelSolution : public StaticKernelSolution<InputT> {
    public:
        typedef std::shared_ptr<RegularizedKernelSolution<InputT> > Ptr;

        RegularizedKernelSolution(lsst::afw::math::KernelList const& basisList,
                                  bool fitForBackground,
                                  Eigen::MatrixXd const& hMat,
                                  lsst::pex::policy::Policy policy
                                  );
        virtual ~RegularizedKernelSolution() {};
        void solve();
        double getLambda() {return _lambda;}
        double estimateRisk(double maxCond);

        /* Include additive term (_lambda * _hMat) in M matrix? */
        Eigen::MatrixXd getM(bool includeHmat = true);

    private:
        Eigen::MatrixXd const _hMat;               
        double _lambda;                                         
        lsst::pex::policy::Policy _policy;

        std::vector<double> _createLambdaSteps();
    };


    class SpatialKernelSolution : public KernelSolution {
    public:
        typedef std::shared_ptr<SpatialKernelSolution> Ptr;

        /* Creates a polynomial SpatialFunction */
        SpatialKernelSolution(lsst::afw::math::KernelList const& basisList,
                              lsst::afw::math::Kernel::SpatialFunctionPtr spatialKernelFunction,
                              lsst::afw::math::Kernel::SpatialFunctionPtr background,
                              lsst::pex::policy::Policy policy
            );

        virtual ~SpatialKernelSolution() {};
        
        void addConstraint(float xCenter, float yCenter,
                           Eigen::MatrixXd const& qMat,
                           Eigen::VectorXd const& wVec);

        void solve();
        std::shared_ptr<lsst::afw::image::Image<lsst::afw::math::Kernel::Pixel>> makeKernelImage(lsst::afw::geom::Point2D const& pos);
        std::pair<std::shared_ptr<lsst::afw::math::LinearCombinationKernel>,
                  lsst::afw::math::Kernel::SpatialFunctionPtr> getSolutionPair();

    private:
        lsst::afw::math::Kernel::SpatialFunctionPtr _spatialKernelFunction; 
        bool _constantFirstTerm;                                            

        std::shared_ptr<lsst::afw::math::LinearCombinationKernel> _kernel;   
        lsst::afw::math::Kernel::SpatialFunctionPtr _background; 
        double _kSum;                                            

        lsst::pex::policy::Policy _policy;                       
        int _nbases;                                             
        int _nkt;                                                
        int _nbt;                                                
        int _nt;                                                 

        void _setKernel();                                       
        void _setKernelUncertainty();                            
    };

}}} // end of namespace lsst::ip::diffim

#endif