lsst::meas::algorithms::shapelet::Function2D Class Reference

#include <Function2D.h>

Inheritance diagram for lsst::meas::algorithms::shapelet::Function2D:

Public Member Functions
	Function2D ()
	Function2D (int xo, int yo)
	Function2D (int xo, int yo, const DMatrix &c)
	Function2D (const Function2D &rhs)
virtual	~Function2D ()
virtual void	write (std::ostream &fout) const =0
virtual std::auto_ptr< Function2D >	copy () const =0
virtual void	addLinear (double a, double b, double c)=0
virtual void	linearPreTransform (double a, double b, double c, double d, double e, double f)=0
virtual std::auto_ptr< Function2D >	dFdX () const =0
virtual std::auto_ptr< Function2D >	dFdY () const =0
virtual std::auto_ptr< Function2D >	conj () const
virtual void	operator*= (double scale)
virtual double	operator() (double x, double y) const
double	operator() (const Position &p) const
virtual void	setTo (double value)
bool	isNonZero () const
int	getXOrder () const
int	getYOrder () const
const DMatrix &	getCoeffs () const
virtual void	simpleFit (int order, const std::vector< Position > &pos, const std::vector< double > &v, const std::vector< bool > &shouldUse, const std::vector< double > *sigList=0, double *chisqOut=0, int *dofOut=0, DMatrix *cov=0)
virtual void	outlierFit (int order, double nsig, const std::vector< Position > &pos, const std::vector< double > &v, std::vector< bool > *shouldUse, const std::vector< double > *sigList=0, double *chisqout=0, int *dofout=0, DMatrix *cov=0)
virtual void	orderFit (int maxOrder, double equivProb, const std::vector< Position > &pos, const std::vector< double > &v, const std::vector< bool > &shouldUse, const std::vector< double > *sigList=0, double *chisqout=0, int *dofout=0, DMatrix *cov=0)
virtual void	linearTransform (double a, double b, double c, const Function2D &f, const Function2D &g)
virtual void	linearTransform (double a, double b, const Function2D &f)
virtual void	operator+= (const Function2D &rhs)=0
virtual void	setFunction (int _xOrder, int _yOrder, const DVector &fVect)=0

Static Public Member Functions
static std::auto_ptr< Function2D >	read (std::istream &fin)

Protected Member Functions
virtual DVector	definePX (int order, double x) const =0
virtual DVector	definePY (int order, double y) const =0

Protected Attributes
int	_xOrder
int	_yOrder
std::auto_ptr< DMatrix >	_coeffs

Private Member Functions
void	doSimpleFit (int xOrder, int yOrder, const std::vector< Position > &pos, const std::vector< double > &v, const std::vector< bool > &shouldUse, DVector *f, const std::vector< double > *sigList=0, int *dof=0, DVector *diff=0, DMatrix *cov=0)

Detailed Description

Definition at line 31 of file Function2D.h.

Constructor & Destructor Documentation

lsst::meas::algorithms::shapelet::Function2D::Function2D ( )

inline

Definition at line 35 of file Function2D.h.

                     : _xOrder(0), _yOrder(0), _coeffs(new DMatrix(1,1)) 
        { (*_coeffs)(0,0) = 0.0; }

lsst::meas::algorithms::shapelet::Function2D::Function2D ( int  xo, int  yo  )

inline

Definition at line 38 of file Function2D.h.

                                   : 
            _xOrder(xo), _yOrder(yo), _coeffs(new DMatrix(xo+1,yo+1)) 
        { _coeffs->setZero(); }

lsst::meas::algorithms::shapelet::Function2D::Function2D ( int  xo, int  yo, const DMatrix &  c  )

inline

Definition at line 42 of file Function2D.h.

                                                     :
            _xOrder(xo),_yOrder(yo),_coeffs(new DMatrix(c)) {}

lsst::meas::algorithms::shapelet::Function2D::Function2D ( const Function2D &  rhs )

inline

Definition at line 45 of file Function2D.h.

                                          :
            _xOrder(rhs._xOrder), _yOrder(rhs._yOrder),
            _coeffs(new DMatrix(*rhs._coeffs)) {}

virtual lsst::meas::algorithms::shapelet::Function2D::~Function2D ( )

inlinevirtual

Definition at line 49 of file Function2D.h.

{}

Member Function Documentation

virtual void lsst::meas::algorithms::shapelet::Function2D::addLinear ( double  a, double  b, double  c  )

pure virtual

Implemented in lsst::meas::algorithms::shapelet::Polynomial2D, lsst::meas::algorithms::shapelet::Constant2D, and lsst::meas::algorithms::shapelet::Legendre2D.

std::auto_ptr< Function2D > lsst::meas::algorithms::shapelet::Function2D::conj ( ) const

virtual

Definition at line 315 of file Function2D.cc.

    {
        std::auto_ptr<Function2D> temp = copy();
        TMV_conjugateSelf(*(temp->_coeffs));
        return temp;
    }

virtual std::auto_ptr<Function2D> lsst::meas::algorithms::shapelet::Function2D::copy ( ) const

pure virtual

Implemented in lsst::meas::algorithms::shapelet::Polynomial2D, lsst::meas::algorithms::shapelet::Constant2D, and lsst::meas::algorithms::shapelet::Legendre2D.

virtual DVector lsst::meas::algorithms::shapelet::Function2D::definePX ( int  order, double  x  ) const

protectedpure virtual

Implemented in lsst::meas::algorithms::shapelet::Polynomial2D, lsst::meas::algorithms::shapelet::Constant2D, and lsst::meas::algorithms::shapelet::Legendre2D.

virtual DVector lsst::meas::algorithms::shapelet::Function2D::definePY ( int  order, double  y  ) const

protectedpure virtual

Implemented in lsst::meas::algorithms::shapelet::Polynomial2D, lsst::meas::algorithms::shapelet::Constant2D, and lsst::meas::algorithms::shapelet::Legendre2D.

virtual std::auto_ptr<Function2D> lsst::meas::algorithms::shapelet::Function2D::dFdX ( ) const

pure virtual

Implemented in lsst::meas::algorithms::shapelet::Polynomial2D, lsst::meas::algorithms::shapelet::Constant2D, and lsst::meas::algorithms::shapelet::Legendre2D.

virtual std::auto_ptr<Function2D> lsst::meas::algorithms::shapelet::Function2D::dFdY ( ) const

pure virtual

Implemented in lsst::meas::algorithms::shapelet::Polynomial2D, lsst::meas::algorithms::shapelet::Constant2D, and lsst::meas::algorithms::shapelet::Legendre2D.

void lsst::meas::algorithms::shapelet::Function2D::doSimpleFit ( int  xOrder, int  yOrder, const std::vector< Position > &  pos, const std::vector< double > &  v, const std::vector< bool > &  shouldUse, DVector *  f, const std::vector< double > *  sigList = 0, int *  dof = 0, DVector *  diff = 0, DMatrix *  cov = 0  )

private

Definition at line 404 of file Function2D.cc.

    {
        // f(x,y) = Sum_pq k_pq px_p py_q
        //        = P . K (where each is a vector in pq)
        // chisq = Sum_n ((v_n - f(x,y))/s_n)^2
        // minchisq => 0 = Sum_n (v_n - f(x,y)) P_pq/s_n^2
        // => [Sum_n P_pq P/s_n^2].K = [Sum_n v_n P_pq/s_n^2]
        // Using ^ to represent an outer product, this can be written:
        //
        //    [Sum_n (P/s_n)^(P/s_n)] . K = [Sum_n (v_n/s_n) (P/s_n)]
        //
        // Or if P' is a matrix with n index for rows, and pq index for columns,
        // with each element P'(n,pq) = P_n,pq/s_n
        // and v' is a vector with v'(n) = v_n/s_n
        //
        // Then, this can be written:
        //
        // P' K = v'
        //
        // The solution to this equation, then, gives our answer for K.
        //
        Assert(pos.size() == vals.size());
        Assert(shouldUse.size() == vals.size());
        Assert((sigList==0) || (sigList->size() == vals.size()));
        xdbg<<"Start SimpleFit: size = "<<pos.size()<<std::endl;
        xdbg<<"order = "<<xOrder<<','<<yOrder<<std::endl;

        int highOrder = std::max(xOrder,yOrder);
        int size = fitSize(xOrder,yOrder);
        xdbg<<"size = "<<size<<std::endl;

        const int nVals = vals.size();

        Assert(int(f->size()) == size);
        Assert(!diff || int(diff->size()) == nVals);
        Assert(!cov || 
               (int(cov->TMV_colsize()) == size && int(cov->TMV_rowsize()) == size));

        int nUse = std::count(shouldUse.begin(),shouldUse.end(),true);
        xdbg<<"nuse = "<<nUse<<std::endl;

        DMatrix P(nUse,size);
        P.setZero();
        DVector V(nUse);

        int ii=0;
        for(int i=0;i<nVals;++i) if (shouldUse[i]) {
            if (sigList) {
                Assert((*sigList)[i] > 0.);
                V(ii) = vals[i]/(*sigList)[i];
            } else {
                V(ii) = vals[i];
            }

            DVector px = definePX(xOrder,pos[i].getX());
            DVector py = definePY(yOrder,pos[i].getY());
            int pq=0;
            for(int pplusq=0;pplusq<=highOrder;++pplusq) { 
                for(int p=std::min(pplusq,xOrder),q=pplusq-p;
                    q<=std::min(pplusq,yOrder);--p,++q,++pq) {
                    Assert(p<int(px.size()));
                    Assert(q<int(py.size()));
                    Assert(pq<int(P.TMV_rowsize()));
                    P(ii,pq) = px[p]*py[q];
                }
            }
            Assert(pq == int(P.TMV_rowsize()));
            if (sigList) P.row(ii) /= (*sigList)[i];
            ++ii;
        }
        Assert(ii==nUse);

        xdbg<<"Done make V,P\n";
        xdbg<<"V = "<<EIGEN_Transpose(V)<<std::endl;
        //P.divideUsing(tmv::QR);
        //P.saveDiv();
        //*f = V/P;
        TMV_QR(P);
        TMV_QR_Solve(P,*f,V);
        xdbg<<"*f = "<<EIGEN_Transpose(*f)<<std::endl;

        if (diff) {
            int k=0;
            for(int i=0;i<nVals;++i) {
                if (shouldUse[i]) {
                    (*diff)(i) = 
                        V(k) - EIGEN_ToScalar(P.row(k) * (*f));
                    ++k;
                } else {
                    (*diff)(i) = 0.;
                }
            }
        }
        if (dof) {
            *dof = P.TMV_colsize() - P.TMV_rowsize();
            if (*dof < 0) *dof = 0;
        }
        if (cov) {
            //P.makeInverseATA(*cov);
            TMV_QR_InverseATA(P,*cov);
        }
        xdbg<<"Done simple fit\n";
    }

const DMatrix& lsst::meas::algorithms::shapelet::Function2D::getCoeffs ( ) const

inline

Definition at line 96 of file Function2D.h.

{ return *_coeffs; }

int lsst::meas::algorithms::shapelet::Function2D::getXOrder ( ) const

inline

Definition at line 92 of file Function2D.h.

{ return _xOrder; }

int lsst::meas::algorithms::shapelet::Function2D::getYOrder ( ) const

inline

Definition at line 94 of file Function2D.h.

{ return _yOrder; }

bool lsst::meas::algorithms::shapelet::Function2D::isNonZero ( ) const

inline

Definition at line 89 of file Function2D.h.

        { return (*_coeffs)(0,0) != 0.0 || _xOrder!=0 || _yOrder!=0; }

virtual void lsst::meas::algorithms::shapelet::Function2D::linearPreTransform ( double  a, double  b, double  c, double  d, double  e, double  f  )

pure virtual

Implemented in lsst::meas::algorithms::shapelet::Polynomial2D, lsst::meas::algorithms::shapelet::Constant2D, and lsst::meas::algorithms::shapelet::Legendre2D.

void lsst::meas::algorithms::shapelet::Function2D::linearTransform ( double  a, double  b, double  c, const Function2D &  f, const Function2D &  g  )

virtual

Definition at line 355 of file Function2D.cc.

    {
        if(dynamic_cast<Constant2D*>(this)) {
            Assert(dynamic_cast<const Constant2D*>(&f));
            Assert(dynamic_cast<const Constant2D*>(&g));
        }
        if(dynamic_cast<Polynomial2D*>(this)) {
            Assert(dynamic_cast<const Polynomial2D*>(&f));
            Assert(dynamic_cast<const Polynomial2D*>(&g));
        }
#ifdef LEGENDRE2D_H
        if(dynamic_cast<Legendre2D*>(this)) {
            Assert(dynamic_cast<const Legendre2D*>(&f));
            Assert(dynamic_cast<const Legendre2D*>(&g));
            Assert(dynamic_cast<const Legendre2D*>(&f)->getBounds() ==
                   dynamic_cast<const Legendre2D*>(&g)->getBounds());
        }
#endif
#ifdef CHEBY2D_H
        if(dynamic_cast<Cheby2D*>(this)) {
            Assert(dynamic_cast<const Cheby2D*>(&f));
            Assert(dynamic_cast<const Cheby2D*>(&g));
            Assert(dynamic_cast<const Cheby2D*>(&f)->getBounds() ==
                   dynamic_cast<const Cheby2D*>(&g)->getBounds());
        }
#endif
        Assert(f.getXOrder() == g.getXOrder());
        Assert(f.getYOrder() == g.getYOrder());

        if (_xOrder != f.getXOrder() || _yOrder != f.getYOrder()) {
            _xOrder = f.getXOrder();
            _yOrder = f.getYOrder();
            _coeffs.reset(new DMatrix(_xOrder+1,_yOrder+1));
            _coeffs->setZero();
        } else _coeffs->setZero();
        for(int i=0;i<=_xOrder;++i) for(int j=0;j<=_yOrder;++j) {
            (*_coeffs)(i,j) = a + b*f.getCoeffs()(i,j) + c*g.getCoeffs()(i,j);
        }
    }

virtual void lsst::meas::algorithms::shapelet::Function2D::linearTransform ( double  a, double  b, const Function2D &  f  )

inlinevirtual

Definition at line 132 of file Function2D.h.

        { linearTransform(a,b,0.,f,f); }

double lsst::meas::algorithms::shapelet::Function2D::operator() ( double  x, double  y  ) const

virtual

Reimplemented in lsst::meas::algorithms::shapelet::Constant2D.

Definition at line 322 of file Function2D.cc.

    {
        DVector px = definePX(_xOrder,x);
        DVector py = definePY(_yOrder,y);
        double result = EIGEN_ToScalar(EIGEN_Transpose(px) * (*_coeffs) * py);
        return result;
    }

double lsst::meas::algorithms::shapelet::Function2D::operator() ( const Position &  p ) const

inline

Definition at line 77 of file Function2D.h.

        { return operator()(p.getX(),p.getY()); }

virtual void lsst::meas::algorithms::shapelet::Function2D::operator*= ( double  scale )

inlinevirtual

Definition at line 72 of file Function2D.h.

        { *_coeffs *= scale; }

virtual void lsst::meas::algorithms::shapelet::Function2D::operator+= ( const Function2D &  rhs )

pure virtual

Implemented in lsst::meas::algorithms::shapelet::Polynomial2D, lsst::meas::algorithms::shapelet::Constant2D, and lsst::meas::algorithms::shapelet::Legendre2D.

void lsst::meas::algorithms::shapelet::Function2D::orderFit ( int  maxOrder, double  equivProb, const std::vector< Position > &  pos, const std::vector< double > &  v, const std::vector< bool > &  shouldUse, const std::vector< double > *  sigList = 0, double *  chisqout = 0, int *  dofout = 0, DMatrix *  cov = 0  )

virtual

Definition at line 602 of file Function2D.cc.

    {
        xdbg<<"Start OrderFit\n";
        DVector fVectMax(fitSize(maxOrder,maxOrder));
        DVector diff(vals.size());
        int dofmax;
        doSimpleFit(maxOrder,maxOrder,pos,vals,shouldUse,
                    &fVectMax,sig,&dofmax,&diff,cov);
        double chisqmax = diff.TMV_normSq();
        xdbg<<"chisq,dof(n="<<maxOrder<<") = "<<chisqmax<<','<<dofmax<<std::endl;
        int tryOrder;
        double chisq=-1.;
        int dof=-1;
        std::auto_ptr<DVector > fVect(0);
        for(tryOrder=0;tryOrder<maxOrder;++tryOrder) {
            fVect.reset(new DVector(fitSize(tryOrder,tryOrder)));
            doSimpleFit(tryOrder,tryOrder,pos,vals,shouldUse,
                        fVect.get(),sig,&dof,&diff,cov);
            chisq = diff.TMV_normSq();
            xdbg<<"chisq,dof(n="<<tryOrder<<") = "<<chisq<<','<<dof<<"....  ";
            if (Equivalent(chisqmax,chisq,dofmax,dof,equivProb)) {
                xdbg<<"equiv\n";
                break;
            }
            xdbg<<"not equiv\n";
        }
        if (tryOrder == maxOrder) {
            setFunction(tryOrder,tryOrder,fVectMax);
            if(chisqOut) *chisqOut = chisqmax;
            if(dofOut) *dofOut = dofmax;
        } else {
            Assert(fVect.get());
            setFunction(tryOrder,tryOrder,*fVect);
            if(chisqOut) *chisqOut = chisq;
            if(dofOut) *dofOut = dof;
        }
    }

void lsst::meas::algorithms::shapelet::Function2D::outlierFit ( int  order, double  nsig, const std::vector< Position > &  pos, const std::vector< double > &  v, std::vector< bool > *  shouldUse, const std::vector< double > *  sigList = 0, double *  chisqout = 0, int *  dofout = 0, DMatrix *  cov = 0  )

virtual

Definition at line 534 of file Function2D.cc.

    {
        xdbg<<"start outlier fit\n";
        const int nVals = vals.size();
        bool isDone=false;
        DVector fVect(fitSize(order,order));
        int dof;
        double chisq=0.;
        double nSigSq = nSig*nSig;
        while (!isDone) {
            DVector diff(nVals);
            xdbg<<"before dosimple\n";
            doSimpleFit(order,order,pos,vals,*shouldUse,&fVect,sig,&dof,&diff,cov);
            xdbg<<"after dosimple\n";
            // Caclulate chisq, keeping the vector diffsq for later when
            // looking for outliers
            chisq = diff.TMV_normSq();
            // If sigmas are given, then chisq should be 1, since diff is
            // already normalized by sig_i.  But if not, then this effectively
            // assumes that all the given errors are off by some uniform factor.

            // Look for outliers, setting isDone = false if any are found
            isDone = true;
            if (dof <= 0) break;
            double thresh=nSigSq*chisq/dof;
            xdbg<<"chisq = "<<chisq<<", thresh = "<<thresh<<std::endl;
            for(int i=0;i<nVals;++i) if( (*shouldUse)[i]) { 
                xdbg<<"pos ="<<pos[i]<<", v = "<<vals[i]<<
                    " diff = "<<diff[i]<<std::endl;
                if (absSq(diff[i]) > thresh) {
                    isDone = false; 
                    (*shouldUse)[i] = false;
                    xdbg<<i<<" ";
                }
            }
            if (!isDone) xdbg<<" are outliers\n";
        }
        setFunction(order,order,fVect);
        if (chisqOut) *chisqOut = chisq;
        if (dofOut) *dofOut = dof;
        xdbg<<"done outlier fit\n";
    }

std::auto_ptr< Function2D > lsst::meas::algorithms::shapelet::Function2D::read ( std::istream &  fin )

static

Definition at line 330 of file Function2D.cc.

    {
        char fc,tc;

        fin >> fc >> tc;
        if (tc != 'D' && tc != 'C') fin.putback(tc);
        std::auto_ptr<Function2D> ret;
        switch(fc) {
          case 'C' : ret.reset(new Constant2D(fin));
                     break;
          case 'P' : ret.reset(new Polynomial2D(fin));
                     break;
#ifdef LEGENDRE2D_H
          case 'L' : ret.reset(new Legendre2D(fin));
                     break;
#endif
#ifdef CHEBY2D_H
          case 'X' : ret.reset(new Cheby2D(fin));
                     break;
#endif
          default: throw std::runtime_error("invalid type");
        }
        return ret;
    }

virtual void lsst::meas::algorithms::shapelet::Function2D::setFunction ( int  _xOrder, int  _yOrder, const DVector &  fVect  )

pure virtual

Implemented in lsst::meas::algorithms::shapelet::Polynomial2D, lsst::meas::algorithms::shapelet::Constant2D, and lsst::meas::algorithms::shapelet::Legendre2D.

virtual void lsst::meas::algorithms::shapelet::Function2D::setTo ( double  value )

inlinevirtual

Definition at line 80 of file Function2D.h.

        {
            if (_xOrder || _yOrder) {
                _xOrder = 0; _yOrder = 0; 
                _coeffs.reset(new DMatrix(1,1));
            }
            (*_coeffs)(0,0) = value;
        }

void lsst::meas::algorithms::shapelet::Function2D::simpleFit ( int  order, const std::vector< Position > &  pos, const std::vector< double > &  v, const std::vector< bool > &  shouldUse, const std::vector< double > *  sigList = 0, double *  chisqOut = 0, int *  dofOut = 0, DMatrix *  cov = 0  )

virtual

Definition at line 513 of file Function2D.cc.

    {
        DVector fVect(fitSize(order,order));
        if (chisqOut) {
            DVector diff(vals.size());
            doSimpleFit(order,order,pos,vals,shouldUse,&fVect,sig,dofOut,&diff,cov);
            *chisqOut = diff.TMV_normSq();
        } else {
            doSimpleFit(order,order,pos,vals,shouldUse,&fVect,sig);
        }
        setFunction(order,order,fVect);
    }

virtual void lsst::meas::algorithms::shapelet::Function2D::write ( std::ostream &  fout ) const

pure virtual

Implemented in lsst::meas::algorithms::shapelet::Polynomial2D, lsst::meas::algorithms::shapelet::Constant2D, and lsst::meas::algorithms::shapelet::Legendre2D.

Member Data Documentation

std::auto_ptr<DMatrix> lsst::meas::algorithms::shapelet::Function2D::_coeffs

protected

Definition at line 146 of file Function2D.h.

int lsst::meas::algorithms::shapelet::Function2D::_xOrder

protected

Definition at line 145 of file Function2D.h.

int lsst::meas::algorithms::shapelet::Function2D::_yOrder

protected

Definition at line 145 of file Function2D.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/Function2D.h
• /home/lsstsw/stack/Linux64/meas_algorithms/11.0-2-gb8b8ce7/src/shapelet/Function2D.cc