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

Last change on this file since 3846 was 3572, checked in by cmv, 17 years ago

char* -> const char* pour regler les problemes de deprecated string const... + comparaison unsigned signed + suppression EVOL_PLANCK rz+cmv 07/02/2009
File size: 7.5 KB
Line 
1// This may look like C code, but it is really -*- C++ -*-
2#ifndef GENERALFIT_SEEN
3#define GENERALFIT_SEEN
4
5#include "pexceptions.h"
6#include "tvector.h"
7#include "generaldata.h"
8#include "generalfunc.h"
9
10namespace SOPHYA {
11
12//================================================================
13// GeneralXi2
14//================================================================
15
16//! Classe de Xi2 a plusieurs parametres
17class GeneralXi2 {
18public:
19 GeneralXi2(unsigned int nPar);
20 virtual ~GeneralXi2();
21
22 /*!
23 Valeur du Xi2 a definir par l'utilisateur (fct virtuelle pure)
24 a partir des donnees de `data'. l'utilisateur doit egalement
25 retourner le nombre de points de mesure utilises dans le calcul
26 du Xi2 (`ndataused').
27 */
28 virtual double Value(GeneralFitData& data, double* parm, int& ndataused)=0;
29 virtual double Derivee(GeneralFitData& data, int i, double* parm);
30 virtual double Derivee2(GeneralFitData& data, int i,int j, double* parm);
31
32 void SetDeltaParm(int numPar, double delta=0.);
33 void SetDeltaParm(double const* dparam);
34
35 //! Retourne le nombre de parametres Ai.
36 inline int NPar() const {return mNPar;}
37
38protected:
39 const int mNPar; //!< nombre de parametres
40
41 double *deltaParm;
42};
43
44//================================================================
45// GENERALFIT
46//================================================================
47
48//! Classe de fit d'une GeneralFunction sur une GeneralFitData
49class GeneralFit {
50public:
51 GeneralFit(GeneralFunction* f);
52 GeneralFit(GeneralXi2* f);
53 ~GeneralFit();
54
55 void WriteStep(const char *filename = NULL);
56 void SetDebug(int level = 0);
57 void SetMaxStep(int n = 100);
58 void SetLambda_Fac(double fac = 10.);
59 void SetStopChi2(double s = 0.01);
60 void SetEps(double ep = 1.e-8);
61 void SetEps(int n,double ep = 1.e-8);
62 void SetStopMx(int nstopmx = 3, double stopchi2 = -1.);
63 void SetStopLent(int nstoplent = 3);
64 void SetFunction(GeneralFunction*);
65 void SetFuncXi2(GeneralXi2*);
66 void SetData(GeneralFitData*);
67 void SetParam(int n,double value, double step
68 ,double min=1., double max=-1.);
69 void SetParam(int n,string const&
70 ,double value, double step
71 ,double min=1., double max=-1.);
72 void SetParam(int n,double value);
73 void SetStep(int n,double value);
74 void SetMinStepDeriv(int i,double val = 0.);
75 void SetMinStepDeriv(double val = 0.);
76 void SetBound(int n,double min,double max);
77 void SetBound(int n);
78 void SetUnBound(int n);
79 void SetUnBound();
80 void SetFix(int n,double v);
81 void SetFix(int n);
82 void SetFree(int n);
83 void SetFree();
84
85 double GetParm(int);
86 TVector<r_8> GetParm();
87 double GetParmErr(int);
88 double GetCoVar(int,int);
89 double GetStep(int n);
90 double GetMax(int n);
91 double GetMin(int n);
92 //! Retourne le Chi2
93 inline double GetChi2() const {return Chi2;};
94 //! Retourne le Chi2 reduit
95 inline double GetChi2Red() const {
96 if(mNddl<=0) return (double) mNddl;
97 return Chi2/(double) mNddl;
98 };
99 //! Retourne la precision de convergence pour le parametre i.
100 inline double GetEps(int i) const {return Eps(i);};
101 //! Retourne le nombre de degres de liberte
102 inline int GetNddl() const {return mNddl;};
103 //! Retourne le nombre d'iterations
104 inline int GetNStep() const {return nStep;};
105 //! Retourne le nombre de variables
106 inline int GetNVar() const {return mNVar;};
107 //! Retourne le nombre de parametres
108 inline int GetNPar() const {return mNPar;};
109 //! Retourne le nombre de parametres libres
110 inline int GetNFree() const {return mNParFree;};
111 //! Retourne le nombre de parametres bornes
112 inline int GetNBound() const {return mNParBound;};
113 //! Retourne le nstop de convergence
114 inline int GetNStop() const {return nStop;};
115 //! Retourne le nstop de convergence lente.
116 inline int GetNStopLent() const {return nStopLent;};
117 //! Retourne le pointeur sur la GeneralFunction utilisee.
118 inline GeneralFunction* GetFunction() const {return mFunction;};
119 //! Retourne le pointeur sur la GeneralFitData utilisee.
120 inline GeneralFitData* GetGData() const {return mData;};
121
122 void PrintStatus();
123 void PrintFit();
124 void PrintParm(int n);
125 void PrintParm();
126
127 int Fit();
128 double ReCalChi2(int& nddl, double* par = NULL);
129 GeneralFitData DataResidus(bool clean=true);
130 GeneralFitData DataFunction(bool clean=true);
131 void PrintFitErr(int rc);
132
133protected:
134 int mNtry; //!< numero d'appel de la routine de fit.
135 int mNVar; //!< nombre de variables f(x,y,z,...)
136 int mNPar; //!< nombre de parametres
137 int mNParFree; //!< nombre de parametres libres
138 int mNParBound; //!< nombre de parametres bornes
139 GeneralFunction* mFunction;
140 GeneralXi2* mFuncXi2;
141 GeneralFitData* mData;
142
143 TVector<r_8> Param;
144 TVector<r_8> errParam;
145 TVector<r_8> stepParam;
146 TVector<r_8> minParam;
147 TVector<r_8> maxParam;
148 TVector<r_8> minStepDeriv;
149 TVector<r_8> Eps;
150 unsigned short int* fixParam;
151 unsigned short int* boundParam;
152 string* nameParam;
153
154 double Lambda_Fac;
155 double stopChi2;
156 int maxStep;
157 int nStopMx;
158 double stopChi2SMx;
159 int nStopLent;
160 int debugLevel;
161 FILE *FileStep;
162
163 TMatrix<r_8> ATGA;
164 TVector<r_8> BETA;
165 TMatrix<r_8> ATGA_Try;
166 TVector<r_8> BETA_Try;
167 TVector<r_8> C;
168 TVector<r_8> D;
169
170 double Chi2;
171 int mNddl;
172 int nStep;
173 int nStop, nStopL;
174 double Lambda;
175
176 // Fonctions privees
177 void write_in_step(double ci2,TVector<r_8>& par);
178 void General_Init(void);
179 void TryFunc(TVector<r_8>& par,TVector<r_8>& par_tr);
180 void TryXi2(TVector<r_8>& par,TVector<r_8>& par_tr);
181 void CheckSanity();
182 void Set_Bound_C_D(int i);
183 void Set_Bound_C_D();
184 double p_vers_tr(int i,double p);
185 TVector<r_8> p_vers_tr(TVector<r_8> const& p);
186 void p_vers_tr(TVector<r_8> const& p,TVector<r_8>& tr);
187 double tr_vers_p(int i,double tr);
188 TVector<r_8> tr_vers_p(TVector<r_8> const& tr);
189 void tr_vers_p(TVector<r_8> const& tr,TVector<r_8>& p);
190 double c_dp_vers_dtr(int i,double tr);
191 TVector<r_8> dp_vers_dtr(TVector<r_8> const& dp,TVector<r_8> const& tr);
192 void dp_vers_dtr(TVector<r_8> const& dp,TVector<r_8> const& tr,TVector<r_8>& dtr);
193 double c_dtr_vers_dp(int i,double tr);
194 TVector<r_8> dtr_vers_dp(TVector<r_8> const& dtr,TVector<r_8> const& tr);
195 int put_in_limits_for_deriv(TVector<r_8> const& p,TVector<r_8>& dp,double dist=0.66);
196 inline void dtr_vers_dp(TVector<r_8> const& dtr,TVector<r_8> const& tr,TVector<r_8>& dp)
197 { for(int i=0;i<mNPar;i++)
198 { if( fixParam[i] ) continue;
199 if( ! boundParam[i] ) dp(i) = dtr(i);
200 else dp(i) = D(i)/(1.+tr(i)*tr(i)) * dtr(i); }
201 };
202};
203
204} // Fin du namespace
205
206#endif
Note: See TracBrowser for help on using the repository browser.