source: Sophya/trunk/SophyaLib/NTools/generalfit.h@ 3853

// This may look like C code, but it is really -*- C++ -*-
#ifndef GENERALFIT_SEEN
#define GENERALFIT_SEEN

#include "pexceptions.h"
#include "tvector.h"
#include "generaldata.h"
#include "generalfunc.h"

namespace SOPHYA {

//================================================================
// GeneralXi2
//================================================================

//! Classe de Xi2 a plusieurs parametres
class GeneralXi2 {
public:
  GeneralXi2(unsigned int nPar);
  virtual ~GeneralXi2();

  /*!
  Valeur du Xi2 a definir par l'utilisateur (fct virtuelle pure)
  a partir des donnees de `data'. l'utilisateur doit egalement
  retourner le nombre de points de mesure utilises dans le calcul
  du Xi2 (`ndataused').
  */
  virtual double Value(GeneralFitData& data, double* parm, int& ndataused)=0;
  virtual double Derivee(GeneralFitData& data, int i, double* parm);
  virtual double Derivee2(GeneralFitData& data, int i,int j, double* parm);

  void SetDeltaParm(int numPar, double delta=0.);
  void SetDeltaParm(double const* dparam);

  //! Retourne le nombre de parametres Ai.
  inline int     NPar() const {return mNPar;}

protected:
  const int mNPar;  //!< nombre de parametres

  double *deltaParm;
};

//================================================================
// GENERALFIT
//================================================================

//! Classe de fit d'une GeneralFunction sur une GeneralFitData
class GeneralFit {
public:
  GeneralFit(GeneralFunction* f);
  GeneralFit(GeneralXi2* f);
  ~GeneralFit();

  void            WriteStep(const char *filename = NULL);
  void            SetDebug(int level = 0);
  void            SetMaxStep(int n = 100);
  void            SetLambda_Fac(double fac = 10.);
  void            SetStopChi2(double s = 0.01);
  void            SetEps(double ep = 1.e-8);
  void            SetEps(int n,double ep = 1.e-8);
  void            SetStopMx(int nstopmx = 3, double stopchi2 = -1.);
  void            SetStopLent(int nstoplent = 3);
  void            SetFunction(GeneralFunction*);
  void            SetFuncXi2(GeneralXi2*);
  void            SetData(GeneralFitData*);
  void            SetParam(int n,double value, double step
                          ,double min=1., double max=-1.);
  void            SetParam(int n,string const&
                          ,double value, double step
                          ,double min=1., double max=-1.);
  void            SetParam(int n,double value);
  void            SetStep(int n,double value);
  void            SetMinStepDeriv(int i,double val = 0.);
  void            SetMinStepDeriv(double val = 0.);
  void            SetBound(int n,double min,double max);
  void            SetBound(int n);
  void            SetUnBound(int n);
  void            SetUnBound();
  void            SetFix(int n,double v);
  void            SetFix(int n);
  void            SetFree(int n);
  void            SetFree();

  double          GetParm(int);
  TVector<r_8>    GetParm();
  double          GetParmErr(int);
  double          GetCoVar(int,int);
  double          GetStep(int n);
  double          GetMax(int n);
  double          GetMin(int n);
  //! Retourne le Chi2
  inline double   GetChi2()    const {return Chi2;};
  //! Retourne le Chi2 reduit
  inline double   GetChi2Red() const {
                        if(mNddl<=0) return (double) mNddl;
                        return Chi2/(double) mNddl;
                                     };
  //! Retourne la precision de convergence pour le parametre i.
  inline double   GetEps(int i) const {return Eps(i);};
  //! Retourne le nombre de degres de liberte
  inline int      GetNddl()    const {return mNddl;};
  //! Retourne le nombre d'iterations
  inline int      GetNStep()   const {return nStep;};
  //! Retourne le nombre de variables
  inline int      GetNVar()    const {return mNVar;};
  //! Retourne le nombre de parametres
  inline int      GetNPar()    const {return mNPar;};
  //! Retourne le nombre de parametres libres
  inline int      GetNFree()   const {return mNParFree;};
  //! Retourne le nombre de parametres bornes
  inline int      GetNBound()  const {return mNParBound;};
  //! Retourne le nstop de convergence
  inline int      GetNStop()   const {return nStop;};
  //! Retourne le nstop de convergence lente.
  inline int      GetNStopLent()   const {return nStopLent;};
  //! Retourne le pointeur sur la GeneralFunction utilisee.
  inline GeneralFunction*  GetFunction() const {return mFunction;};
  //! Retourne le pointeur sur la GeneralFitData utilisee.
  inline GeneralFitData*   GetGData() const {return mData;};

  void            PrintStatus();
  void            PrintFit();
  void            PrintParm(int n);
  void            PrintParm();

  int             Fit();
  double          ReCalChi2(int& nddl, double* par = NULL);
  GeneralFitData  DataResidus(bool clean=true);
  GeneralFitData  DataFunction(bool clean=true);
  void            PrintFitErr(int rc);

protected:
  int             mNtry;       //!< numero d'appel de la routine de fit.
  int             mNVar;       //!< nombre de variables f(x,y,z,...)
  int             mNPar;       //!< nombre de parametres
  int             mNParFree;   //!< nombre de parametres libres
  int             mNParBound;  //!< nombre de parametres bornes
  GeneralFunction*  mFunction;
  GeneralXi2*       mFuncXi2;
  GeneralFitData*   mData;

  TVector<r_8>    Param;
  TVector<r_8>    errParam;
  TVector<r_8>    stepParam;
  TVector<r_8>    minParam;
  TVector<r_8>    maxParam;
  TVector<r_8>    minStepDeriv;
  TVector<r_8>    Eps;
  unsigned short int* fixParam;
  unsigned short int* boundParam;
  string*             nameParam;
  
  double          Lambda_Fac;
  double          stopChi2;
  int             maxStep;
  int             nStopMx;
  double          stopChi2SMx;
  int             nStopLent;
  int             debugLevel;
  FILE            *FileStep;
  
  TMatrix<r_8>    ATGA;
  TVector<r_8>    BETA;
  TMatrix<r_8>    ATGA_Try;
  TVector<r_8>    BETA_Try;
  TVector<r_8>    C;
  TVector<r_8>    D;
  
  double          Chi2;
  int             mNddl;
  int             nStep;
  int             nStop, nStopL;
  double          Lambda;
  
  // Fonctions privees
  void         write_in_step(double ci2,TVector<r_8>& par);
  void         General_Init(void);
  void         TryFunc(TVector<r_8>& par,TVector<r_8>& par_tr);
  void         TryXi2(TVector<r_8>& par,TVector<r_8>& par_tr);
  void         CheckSanity();
  void         Set_Bound_C_D(int i);
  void         Set_Bound_C_D();
  double       p_vers_tr(int i,double p);
  TVector<r_8> p_vers_tr(TVector<r_8> const& p);
  void         p_vers_tr(TVector<r_8> const& p,TVector<r_8>& tr);
  double       tr_vers_p(int i,double tr);
  TVector<r_8> tr_vers_p(TVector<r_8> const& tr);
  void         tr_vers_p(TVector<r_8> const& tr,TVector<r_8>& p);
  double       c_dp_vers_dtr(int i,double tr);
  TVector<r_8> dp_vers_dtr(TVector<r_8> const& dp,TVector<r_8> const& tr);
  void         dp_vers_dtr(TVector<r_8> const& dp,TVector<r_8> const& tr,TVector<r_8>& dtr);
  double       c_dtr_vers_dp(int i,double tr);
  TVector<r_8> dtr_vers_dp(TVector<r_8> const& dtr,TVector<r_8> const& tr);
  int          put_in_limits_for_deriv(TVector<r_8> const& p,TVector<r_8>& dp,double dist=0.66);
  inline void dtr_vers_dp(TVector<r_8> const& dtr,TVector<r_8> const& tr,TVector<r_8>& dp)
              {  for(int i=0;i<mNPar;i++)
                   { if( fixParam[i] ) continue;
                     if( ! boundParam[i] ) dp(i) = dtr(i);
                       else dp(i) = D(i)/(1.+tr(i)*tr(i)) * dtr(i); }
              };
};

} // Fin du namespace

#endif