smoothmodelfitter.h

// -*- C++ -*-
// 
// 
#ifndef SMOOTHMODELFITTER_H
#define SMOOTHMODELFITTER_H


#include <time.h>
#include <math.h>

#include "ccdimage.h"
#include "measuredstar.h"

class Associations;
class Vect;


// 2d polynomial
class Poly2 
{
public:
  Poly2(const double XMin, const double YMin, 
    const double XMax, const double YMax,
    const int Deg);
  
  void      Monomials(const double &X, const double &Y, Vect &M) const;
  
  bool      Compute(const double X, const double Y, Vect const& coeffs, double& val) const;
  
  //  void      SetCoeffs(const Vect &Coeffs);
  
  
  unsigned  NTerms() const { return nterms_;}
  
  bool      Read(std::istream &S);
  void      Write(std::ostream &S) const;
  
  void      DecreaseDegree();
  
private:
  double ax_, bx_;
  double ay_, by_;
  int    deg_;
  unsigned nterms_;
  Vect coeffs_;
};


class SmoothModelFitter
{
public:
  SmoothModelFitter(Associations& assoc, unsigned int deg);
  
  ~SmoothModelFitter() { }
  
  template<class Op>
  void    LoopOnMeasurements(Op& op) const;
  
  bool    Minimize(unsigned int maxiter);
  
  void    RemoveOutliers(double nsig);
  
  double  GetChi2() const;
  
  int     GetNMeas() const;
  
  int     GetNDoF() const;
  
  int     GetNPar() const { return npar_; }
  
  static bool Value(MeasuredStar& m, Vect const& pars, double& value, int deg);
  
  void    DumpNTuples(std::string const& output_dir) const;
  
  void    DumpParNTuple(std::string const& ntname) const;
  
  void    DumpMeasuredStarNTuple(std::string const& ntname) const;
  
  void    DumpFittedStarNTuple(std::string const& ntname) const;
  
private:
  unsigned int     deg_;
  Associations&    assoc_;
  CcdImageList&    imstd::list_;
  FittedStarList&  fittedstarstd::list_;
  int              npar_;
  
  Poly2            corrs_;
  Vect             calibpars_;
  
  mutable double   chi2_;
  mutable int      nmeas_;
  mutable int      ndof_;
  mutable bool     cleanChi2_;
  
  void update_star_fluxes_();
  void update_calib_pars_();
};



template<class Op>
void SmoothModelFitter::LoopOnMeasurements(Op& op) const
{
  CcdImageIterator imI;
  MeasuredStarIterator stI;
  
  for(imI=imstd::list_.begin();imI!=imstd::list_.end();imI++)
    {
      CcdImage& ccdImage = **imI;
      MeasuredStarList& catalog = ccdImage.CatalogForFit();
      for(stI=catalog.begin();stI!=catalog.end();stI++)
    {
      MeasuredStar& m = **stI;
      if( !m.IsValid() ) continue;
      FittedStar*   f = m.GetFittedStar();
      if(!f) continue;
      op(m);
    }
    }
  
  op.Finalize();
}




#endif