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

// 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"

namespace SOPHYA {

//================================================================
// GeneralFunction
//================================================================


/*!
  Classe de fonctions parametrees a plusieurs variables:
  \f$ F[x1,x2,x3,...:a1,a2,a3,...] \f$
*/
class GeneralFunction {
public:
  GeneralFunction(unsigned int nVar, unsigned int nPar);
  virtual ~GeneralFunction();

  //! Valeur de la fonction a definir par l'utilisateur (fct virtuelle pure)
  virtual double Value(double const xp[], double const* parm)=0;
  virtual double Val_Der(double const xp[], double const* parm
                        , double* DgDpar);

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

  //! Retourne le nombre de variables Xi
  inline int     NVar() const {return mNVar;}
  //! Retourne le nombre de parametres Ai
  inline int     NPar() const {return mNPar;}

protected:
  const int mNVar;  //!< nombre de variables f(x,y,z,...)
  const int mNPar;  //!< nombre de parametres

  double *deltaParm;
  double *tmpParm;
};

//================================================================
// GeneralFunc
//================================================================

/*!
  Classe de fonctions parametrees a plusieurs variables
  derivant de ``GeneralFunction''. Permet de definir
  une fonction a fiter sans passer par une classe derivee
  en utilisant l'ecriture courante du C. La fonction
  retournant les derivees par rapport aux parametres du fit
  peut etre egalement fournie (optionnel).
*/
class GeneralFunc : public GeneralFunction {
public:
  GeneralFunc(unsigned int nvar, unsigned int npar
        ,double (*fun)(double const*,double const*) 
        ,double (*funder)(double const*, double const*, double*)=NULL);
  virtual ~GeneralFunc();

  virtual double Value(double const xp[], double const* Par);
  virtual double Val_Der(double const xp[],double const* parm
                        , double* DgDpar);

protected:
  double (*tmpFun)   (double const*, double const*);
  double (*tmpFunDer)(double const*, double const*, double*);
};

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

class GeneralFitData;

/*!
  Classe de Xi2 a plusieurs parametres :
  \f$ Xi2[a1,a2,a3,...] \f$
*/
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(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);
  Vector          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;

  Vector          Param;
  Vector          errParam;
  Vector          stepParam;
  Vector          minParam;
  Vector          maxParam;
  Vector          minStepDeriv;
  Vector          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;
  
  Matrix          ATGA;
  Vector          BETA;
  Matrix          ATGA_Try;
  Vector          BETA_Try;
  Vector          C;
  Vector          D;
  
  double          Chi2;
  int             mNddl;
  int             nStep;
  int             nStop, nStopL;
  double          Lambda;
  
  // Fonctions privees
  void       write_in_step(double ci2,Vector& par);
  void       General_Init(void);
  void       TryFunc(Vector& par,Vector& par_tr);
  void       TryXi2(Vector& par,Vector& par_tr);
  void       CheckSanity();
  void       Set_Bound_C_D(int i);
  void       Set_Bound_C_D();
  double     p_vers_tr(int i,double p);
  Vector     p_vers_tr(Vector const& p);
  void       p_vers_tr(Vector const& p,Vector& tr);
  double     tr_vers_p(int i,double tr);
  Vector     tr_vers_p(Vector const& tr);
  void       tr_vers_p(Vector const& tr,Vector& p);
  double     c_dp_vers_dtr(int i,double tr);
  Vector       dp_vers_dtr(Vector const& dp,Vector const& tr);
  void         dp_vers_dtr(Vector const& dp,Vector const& tr,Vector& dtr);
  double     c_dtr_vers_dp(int i,double tr);
  Vector       dtr_vers_dp(Vector const& dtr,Vector const& tr);
  int        put_in_limits_for_deriv(Vector const& p,Vector& dp,double dist=0.66);
  inline void dtr_vers_dp(Vector const& dtr,Vector const& tr,Vector& 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