lsst::meas::algorithms::shapelet::FittedPsf Class Reference

#include <FittedPsf.h>

Public Member Functions
	FittedPsf (const ConfigFile &params)
int	getPsfOrder () const
int	getPsfSize () const
int	getFitOrder () const
int	getFitSize () const
int	getNpca () const
double	getSigma () const
void	setSigma (double sigma)
double	getXMin () const
double	getXMax () const
double	getYMin () const
double	getYMax () const
const Bounds &	getBounds () const
void	calculate (const std::vector< Position > &pos, const std::vector< BVec > &psf, const std::vector< double > &nu, std::vector< long > &flags)
void	interpolate (Position pos, BVec &b) const
double	interpolateSingleElement (Position pos, int i) const
FittedPsfAtXY	operator() (Position pos) const
BVec	getMean () const

Private Member Functions
void	interpolateVector (Position pos, DVectorView b) const

Private Attributes
const ConfigFile &	_params
int	_psfOrder
double	_sigma
int	_fitOrder
int	_fitSize
int	_nPca
Bounds	_bounds
std::auto_ptr< DVector >	_avePsf
std::auto_ptr< DMatrix >	_mV_transpose
std::auto_ptr< DMatrix >	_f

Friends
class	FittedPsfAtXY

Detailed Description

Definition at line 20 of file FittedPsf.h.

Constructor & Destructor Documentation

lsst::meas::algorithms::shapelet::FittedPsf::FittedPsf

(

const ConfigFile &

params

)

Definition at line 309 of file FittedPsf.cc.

                                                 : 
        _params(params), _psfOrder(_params.read<int>("psf_order")),
        _fitOrder(_params.read<int>("fitpsf_order")),
        _fitSize((_fitOrder+1)*(_fitOrder+2)/2)
    { }

Member Function Documentation

void lsst::meas::algorithms::shapelet::FittedPsf::calculate	(	const std::vector< Position > &	pos,
		const std::vector< BVec > &	psf,
		const std::vector< double > &	nu,
		std::vector< long > &	flags
	)

Definition at line 83 of file FittedPsf.cc.

    {
        const int nStars = pos.size();
        const int psfSize = (_psfOrder+1)*(_psfOrder+2)/2;
        const double nSigmaClip = _params.read("fitpsf_nsigma_outlier",3);

        // This is an empirical fit to the chisq level that corresponds to 
        // a 3-sigma outlier for more than 1 dimension.
        // I calculated the values up to n=100 and for n>30, they form a pretty
        // good approximation to a straight line.
        // This is almost certainly wrong for nSigmaClip != 3, so if we start 
        // choosing other values for nSigmaClip, it might be worth doing this 
        // right.
        // That means calculating the 1-d critical value for the given nSigma.
        // e.g. nSigma = 3 -> alpha = P(chisq > 9) = 0.0027.
        // Then calculate the critical value of chisq for that alpha with 
        // the full degrees of freedom = psfSize
        const double chisqLevel = 0.14*psfSize + 2.13;

        const double outlierThresh = nSigmaClip * nSigmaClip * chisqLevel;
        dbg<<"outlierThresh = "<<outlierThresh<<std::endl;

        int nGoodPsf, nOutliers, dof;
        double chisq;
        do {

            // Calculate the average psf vector
            _avePsf.reset(new DVector(psfSize));
            Assert(psfSize == int(_avePsf->size()));
            _avePsf->setZero();
            nGoodPsf = 0;
            for(int n=0;n<nStars;++n) if ( flags[n]==0 ) {
                Assert(psf[n].getSigma() == _sigma);
                *_avePsf += psf[n].vec();
                ++nGoodPsf;
            }
            if (nGoodPsf == 0) {
                dbg<<"ngoodpsf = 0 in FittedPsf::calculate\n";
                throw ProcessingException("No good stars found for interpolation.");
            }
            *_avePsf /= double(nGoodPsf);

            // Rotate the vectors into their eigen directions.
            // The matrix V is stored to let us get back to the original basis.
            DMatrix mM(nGoodPsf,psfSize);
            DDiagMatrix inverseSigma(nGoodPsf);
            int i=0;
            for(int n=0;n<nStars;++n) if ( flags[n]==0 ) {
                Assert(int(psf[n].size()) == psfSize);
                Assert(i < nGoodPsf);
                mM.row(i) = psf[n].vec() - *_avePsf;
                inverseSigma(i) = nu[n];
                _bounds += pos[n];
                ++i;
            }
            Assert(i == nGoodPsf);
            xdbg<<"bounds = "<<_bounds<<std::endl;
            mM = inverseSigma EIGEN_asDiag() * mM;

            int nPcaTot = std::min(nGoodPsf,psfSize);
            DDiagMatrix mS(nPcaTot);
#ifdef USE_TMV
            DMatrixView mU = mM.colRange(0,nPcaTot);
            _mV.reset(new tmv::Matrix<double,tmv::RowMajor>(nPcaTot,psfSize));
            if (nGoodPsf > psfSize) {
                SV_Decompose(mU.view(),mS.view(),_mV->view(),true);
            } else {
                *_mV = mM;
                SV_Decompose(_mV->transpose(),mS.view(),mU.transpose());
            }
            xdbg<<"In FittedPSF: SVD S = "<<mS.diag()<<std::endl;
#else
            DMatrix mU(mM.TMV_colsize(),nPcaTot);
            _mV_transpose.reset(new DMatrix(psfSize,nPcaTot));
            if (nGoodPsf > psfSize) {
                Eigen::JacobiSVD<DMatrix> svd(
                    TMV_colRange(mM,0,nPcaTot), Eigen::ComputeThinU | Eigen::ComputeThinV
                );
                mU = svd.matrixU();
                mS = svd.singularValues();
                *_mV_transpose = svd.matrixV().transpose();
            } else {
                Eigen::JacobiSVD<DMatrix> svd(
                    mM, Eigen::ComputeThinU | Eigen::ComputeThinV
                );
                mU = svd.matrixU();
                mS = svd.singularValues();
                *_mV_transpose = svd.matrixV().transpose();
            }
            xdbg<<"In FittedPSF: SVD S = "<<EIGEN_Transpose(mS)<<std::endl;
#endif
            if (_params.keyExists("fitpsf_npca")) {
                _nPca = _params["fitpsf_npca"];
                dbg<<"npca = "<<_nPca<<" from parameter file\n";
            } else {
                double thresh = mS(0);
                if (_params.keyExists("fitpsf_pca_thresh")) 
                    thresh *= double(_params["fitpsf_pca_thresh"]);
                else thresh *= std::numeric_limits<double>::epsilon();
                dbg<<"thresh = "<<thresh<<std::endl;
                for(_nPca=1;_nPca<int(mM.TMV_rowsize());++_nPca) {
                    if (mS(_nPca) < thresh) break;
                }
                dbg<<"npca = "<<_nPca<<std::endl;
            }
#ifdef USE_TMV
            mU.colRange(0,_nPca) *= mS.subDiagMatrix(0,_nPca);
#else
            TMV_colRange(mU,0,_nPca) *= mS.TMV_subVector(0,_nPca).asDiagonal();
#endif
            xdbg<<"After U *= S\n";
            // U S = M(orig) * Vt

            while (nGoodPsf <= _fitSize && _fitSize > 1) {
                --_fitOrder;
                _fitSize = (_fitOrder+1)*(_fitOrder+2)/2;
                dbg<<"Too few good stars... reducing order of fit to "<<
                    _fitOrder<<std::endl;
            }
            DMatrix mP(nGoodPsf,_fitSize);
            mP.setZero();
            i=0;
            for(int n=0;n<nStars;++n) if ( flags[n]==0 ) {
                xdbg<<"n = "<<n<<" / "<<nStars<<std::endl;
#ifdef USE_TMV
                setPRow(_fitOrder,pos[n],_bounds,mP.row(i));
#else
                DVector mProwi(mP.TMV_rowsize());
                setPRow(_fitOrder,pos[n],_bounds,mProwi);
                mP.row(i) = mProwi.transpose();
#endif
                ++i;
            }
            Assert(i == nGoodPsf);
            mP = inverseSigma EIGEN_asDiag() * mP;
            xdbg<<"after mP = sigma * mP\n";

#ifdef USE_TMV
            _f.reset(new DMatrix(TMV_colRange(mU,0,_nPca)/mP));
#else
            _f.reset(new DMatrix(_fitSize,_nPca));
            *_f = mP.colPivHouseholderQr().solve(TMV_colRange(mU,0,_nPca));
#endif
            xdbg<<"Done making FittedPSF\n";

            //
            // Remove outliers from the fit using the empirical covariance matrix
            // of the data with respect to the fitted values.
            //
            xdbg<<"Checking for outliers:\n";

            // Calculate the covariance matrix
            DMatrix cov(psfSize,psfSize);
            cov.setZero();
            for(int n=0;n<nStars;++n) if ( flags[n]==0 ) {
                const BVec& data = psf[n];
                DVector fit(psfSize);
                interpolateVector(pos[n],TMV_vview(fit));
                DVector diff = data.vec() - fit;
#ifdef USE_TMV
                cov += diff ^ diff;
#else
                cov += diff * diff.transpose();
#endif
            }

            chisq = (mP * *_f - TMV_colRange(mU,0,_nPca)).TMV_normSq();
            dbg<<"chisq calculation #1 = "<<chisq<<std::endl;
            dof = nGoodPsf - _fitSize;

            if (dof > 0) { 
                cov /= double(dof);
            }
#ifdef USE_TMV
            cov.divideUsing(tmv::SV);
            cov.saveDiv();
            cov.setDiv();
            dbg<<"cov S = "<<cov.svd().getS().diag()<<std::endl;
#else
            Eigen::JacobiSVD<DMatrix> cov_svd(cov, Eigen::ComputeThinU | Eigen::ComputeThinV);
            dbg<<"cov S = "<<EIGEN_Transpose(cov_svd.singularValues())<<std::endl;

#endif

            // Clip out 3 sigma outliers:
            nOutliers = 0;
            chisq = 0;
            for(int n=0;n<nStars;++n) if ( flags[n]==0 ) {
                const BVec& data = psf[n];
                DVector fit(psfSize);
                interpolateVector(pos[n],TMV_vview(fit));
                DVector diff = data.vec() - fit;
#ifdef USE_TMV
                double dev = diff * cov.inverse() * diff;
#else
                DVector temp(psfSize);
                temp = cov_svd.solve(diff);
                double dev = (diff.transpose() * temp)(0,0);
#endif
                chisq += dev;
                if (dev > outlierThresh) {
                    xdbg<<"n = "<<n<<" is an outlier.\n";
                    xdbg<<"data = "<<data.vec()<<std::endl;
                    xdbg<<"fit = "<<fit<<std::endl;
                    xdbg<<"diff = "<<diff<<std::endl;
#ifdef USE_TMV
                    xdbg<<"diff/cov = "<<diff/cov<<std::endl;
#else
                    xdbg<<"diff/cov = "<<temp<<std::endl;
#endif
                    xdbg<<"dev = "<<dev<<std::endl;
                    ++nOutliers;
                    flags[n] |= PSF_INTERP_OUTLIER;
                }
            }
            dbg<<"ngoodpsf = "<<nGoodPsf<<std::endl;
            dbg<<"nOutliers = "<<nOutliers<<std::endl;
            dbg<<"chisq calculation #2 = "<<chisq<<std::endl;

        } while (nOutliers > 0);
    }

const Bounds& lsst::meas::algorithms::shapelet::FittedPsf::getBounds ( ) const

inline

Definition at line 45 of file FittedPsf.h.

{ return _bounds; }

int lsst::meas::algorithms::shapelet::FittedPsf::getFitOrder ( ) const

inline

Definition at line 35 of file FittedPsf.h.

{ return _fitOrder; }

int lsst::meas::algorithms::shapelet::FittedPsf::getFitSize ( ) const

inline

Definition at line 36 of file FittedPsf.h.

{ return _fitSize; }

BVec lsst::meas::algorithms::shapelet::FittedPsf::getMean

(

)

const

int lsst::meas::algorithms::shapelet::FittedPsf::getNpca ( ) const

inline

Definition at line 37 of file FittedPsf.h.

{ return _nPca; }

int lsst::meas::algorithms::shapelet::FittedPsf::getPsfOrder ( ) const

inline

Definition at line 33 of file FittedPsf.h.

{ return _psfOrder; }

int lsst::meas::algorithms::shapelet::FittedPsf::getPsfSize ( ) const

inline

Definition at line 34 of file FittedPsf.h.

{ return (_psfOrder+1)*(_psfOrder+2)/2; }

double lsst::meas::algorithms::shapelet::FittedPsf::getSigma ( ) const

inline

Definition at line 38 of file FittedPsf.h.

{ return _sigma; }

double lsst::meas::algorithms::shapelet::FittedPsf::getXMax ( ) const

inline

Definition at line 42 of file FittedPsf.h.

{ return _bounds.getXMax(); }

double lsst::meas::algorithms::shapelet::FittedPsf::getXMin ( ) const

inline

Definition at line 41 of file FittedPsf.h.

{ return _bounds.getXMin(); }

double lsst::meas::algorithms::shapelet::FittedPsf::getYMax ( ) const

inline

Definition at line 44 of file FittedPsf.h.

{ return _bounds.getYMax(); }

double lsst::meas::algorithms::shapelet::FittedPsf::getYMin ( ) const

inline

Definition at line 43 of file FittedPsf.h.

{ return _bounds.getYMin(); }

void lsst::meas::algorithms::shapelet::FittedPsf::interpolate ( Position  pos, BVec &  b  ) const

inline

Definition at line 65 of file FittedPsf.h.

        {
            Assert(_avePsf.get());
            //Assert(_mV.get());
            Assert(b.getOrder() == _psfOrder);
            Assert(static_cast<int>(b.size()) == _avePsf->size());
            b.setSigma(_sigma);
            interpolateVector(pos,TMV_vview(b.vec()));
        }

double lsst::meas::algorithms::shapelet::FittedPsf::interpolateSingleElement	(	Position	pos,
		int	i
	)		const

Definition at line 315 of file FittedPsf.cc.

    {
        DVector P(_fitSize);
#ifdef USE_TMV
        setPRow(_fitOrder,pos,_bounds,P.view());
        DVector b1 = P * (*_f);
        double bi = b1 * _mV->col(i,0,_nPca);
#else
        setPRow(_fitOrder,pos,_bounds,P);
        DVector b1 = _f->transpose() * P;
        double bi = EIGEN_ToScalar(
            EIGEN_Transpose(b1) * _mV_transpose->TMV_rowpart(i,0,_nPca));
#endif
        bi += (*_avePsf)(i);
        return bi;
    }

void lsst::meas::algorithms::shapelet::FittedPsf::interpolateVector ( Position  pos, DVectorView  b  ) const

private

Definition at line 332 of file FittedPsf.cc.

    {
        DVector P(_fitSize);
#ifdef USE_TMV
        setPRow(_fitOrder,pos,_bounds,P.view());
        DVector b1 = P * (*_f);
        b = b1 * _mV->rowRange(0,_nPca);
#else
        setPRow(_fitOrder,pos,_bounds,P);
        DVector b1 = _f->transpose() * P;
        b = TMV_colRange(*_mV_transpose,0,_nPca) * b1;
#endif
        b += *_avePsf;
    }

FittedPsfAtXY lsst::meas::algorithms::shapelet::FittedPsf::operator() ( Position  pos ) const

inline

Definition at line 130 of file FittedPsf.h.

    { return FittedPsfAtXY(*this,pos); }

void lsst::meas::algorithms::shapelet::FittedPsf::setSigma ( double  sigma )

inline

Definition at line 39 of file FittedPsf.h.

{ _sigma = sigma; }

Friends And Related Function Documentation

friend class FittedPsfAtXY

friend

Definition at line 82 of file FittedPsf.h.

Member Data Documentation

std::auto_ptr<DVector> lsst::meas::algorithms::shapelet::FittedPsf::_avePsf

private

Definition at line 99 of file FittedPsf.h.

Bounds lsst::meas::algorithms::shapelet::FittedPsf::_bounds

private

Definition at line 98 of file FittedPsf.h.

std::auto_ptr<DMatrix> lsst::meas::algorithms::shapelet::FittedPsf::_f

private

Definition at line 105 of file FittedPsf.h.

int lsst::meas::algorithms::shapelet::FittedPsf::_fitOrder

private

Definition at line 95 of file FittedPsf.h.

int lsst::meas::algorithms::shapelet::FittedPsf::_fitSize

private

Definition at line 96 of file FittedPsf.h.

std::auto_ptr<DMatrix> lsst::meas::algorithms::shapelet::FittedPsf::_mV_transpose

private

Definition at line 103 of file FittedPsf.h.

int lsst::meas::algorithms::shapelet::FittedPsf::_nPca

private

Definition at line 97 of file FittedPsf.h.

const ConfigFile& lsst::meas::algorithms::shapelet::FittedPsf::_params

private

Definition at line 91 of file FittedPsf.h.

int lsst::meas::algorithms::shapelet::FittedPsf::_psfOrder

private

Definition at line 93 of file FittedPsf.h.

double lsst::meas::algorithms::shapelet::FittedPsf::_sigma

private

Definition at line 94 of file FittedPsf.h.

---

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

• /home/lsstsw/stack/Linux64/meas_algorithms/11.0-2-gb8b8ce7/include/lsst/meas/algorithms/shapelet/FittedPsf.h
• /home/lsstsw/stack/Linux64/meas_algorithms/11.0-2-gb8b8ce7/src/shapelet/FittedPsf.cc