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

Last change on this file since 2870 was 938, checked in by ansari, 25 years ago

Vector,Matrix -> TVector,TMatrix<r_8> cmv 14/4/00
File size: 9.3 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
9namespace SOPHYA {
10
11//================================================================
12// GeneralFunction
13//================================================================
14
15//! Classe de fonctions parametrees a plusieurs variables
16class GeneralFunction {
17public:
18 GeneralFunction(unsigned int nVar, unsigned int nPar);
19 virtual ~GeneralFunction();
20
21 //! Valeur de la fonction a definir par l'utilisateur (fct virtuelle pure)
22 virtual double Value(double const xp[], double const* parm)=0;
23 //! Valeur de la fonction derivee selon les parametres pouvant etre redefinie
24 virtual double Val_Der(double const xp[], double const* parm
25 , double* DgDpar);
26
27 void SetDeltaParm(int numPar, double delta=0.);
28 void SetDeltaParm(double const* dparam);
29
30 //! Retourne le nombre de variables Xi
31 inline int NVar() const {return mNVar;}
32 //! Retourne le nombre de parametres Ai
33 inline int NPar() const {return mNPar;}
34
35protected:
36 const int mNVar; //!< nombre de variables f(x,y,z,...)
37 const int mNPar; //!< nombre de parametres
38
39 double *deltaParm;
40 double *tmpParm;
41};
42
43//================================================================
44// GeneralFunc
45//================================================================
46
47//! Classe de fonctions parametrees a plusieurs variables type C
48class GeneralFunc : public GeneralFunction {
49public:
50 GeneralFunc(unsigned int nvar, unsigned int npar
51 ,double (*fun)(double const*,double const*)
52 ,double (*funder)(double const*, double const*, double*)=NULL);
53 virtual ~GeneralFunc();
54
55 virtual double Value(double const xp[], double const* Par);
56 virtual double Val_Der(double const xp[],double const* parm
57 , double* DgDpar);
58
59protected:
60 double (*tmpFun) (double const*, double const*);
61 double (*tmpFunDer)(double const*, double const*, double*);
62};
63
64//================================================================
65// GeneralXi2
66//================================================================
67
68class GeneralFitData;
69
70//! Classe de Xi2 a plusieurs parametres
71class GeneralXi2 {
72public:
73 GeneralXi2(unsigned int nPar);
74 virtual ~GeneralXi2();
75
76 /*!
77 Valeur du Xi2 a definir par l'utilisateur (fct virtuelle pure)
78 a partir des donnees de `data'. l'utilisateur doit egalement
79 retourner le nombre de points de mesure utilises dans le calcul
80 du Xi2 (`ndataused').
81 */
82 virtual double Value(GeneralFitData& data, double* parm, int& ndataused)=0;
83 virtual double Derivee(GeneralFitData& data, int i, double* parm);
84 virtual double Derivee2(GeneralFitData& data, int i,int j, double* parm);
85
86 void SetDeltaParm(int numPar, double delta=0.);
87 void SetDeltaParm(double const* dparam);
88
89 //! Retourne le nombre de parametres Ai.
90 inline int NPar() const {return mNPar;}
91
92protected:
93 const int mNPar; //!< nombre de parametres
94
95 double *deltaParm;
96};
97
98//================================================================
99// GENERALFIT
100//================================================================
101
102//! Classe de fit d'une GeneralFunction sur une GeneralFitData
103class GeneralFit {
104public:
105 GeneralFit(GeneralFunction* f);
106 GeneralFit(GeneralXi2* f);
107 ~GeneralFit();
108
109 void WriteStep(char *filename = NULL);
110 void SetDebug(int level = 0);
111 void SetMaxStep(int n = 100);
112 void SetLambda_Fac(double fac = 10.);
113 void SetStopChi2(double s = 0.01);
114 void SetEps(double ep = 1.e-8);
115 void SetEps(int n,double ep = 1.e-8);
116 void SetStopMx(int nstopmx = 3, double stopchi2 = -1.);
117 void SetStopLent(int nstoplent = 3);
118 void SetFunction(GeneralFunction*);
119 void SetFuncXi2(GeneralXi2*);
120 void SetData(GeneralFitData*);
121 void SetParam(int n,double value, double step
122 ,double min=1., double max=-1.);
123 void SetParam(int n,string const&
124 ,double value, double step
125 ,double min=1., double max=-1.);
126 void SetParam(int n,double value);
127 void SetStep(int n,double value);
128 void SetMinStepDeriv(int i,double val = 0.);
129 void SetMinStepDeriv(double val = 0.);
130 void SetBound(int n,double min,double max);
131 void SetBound(int n);
132 void SetUnBound(int n);
133 void SetUnBound();
134 void SetFix(int n,double v);
135 void SetFix(int n);
136 void SetFree(int n);
137 void SetFree();
138
139 double GetParm(int);
140 TVector<r_8> GetParm();
141 double GetParmErr(int);
142 double GetCoVar(int,int);
143 double GetStep(int n);
144 double GetMax(int n);
145 double GetMin(int n);
146 //! Retourne le Chi2
147 inline double GetChi2() const {return Chi2;};
148 //! Retourne le Chi2 reduit
149 inline double GetChi2Red() const {
150 if(mNddl<=0) return (double) mNddl;
151 return Chi2/(double) mNddl;
152 };
153 //! Retourne la precision de convergence pour le parametre i.
154 inline double GetEps(int i) const {return Eps(i);};
155 //! Retourne le nombre de degres de liberte
156 inline int GetNddl() const {return mNddl;};
157 //! Retourne le nombre d'iterations
158 inline int GetNStep() const {return nStep;};
159 //! Retourne le nombre de variables
160 inline int GetNVar() const {return mNVar;};
161 //! Retourne le nombre de parametres
162 inline int GetNPar() const {return mNPar;};
163 //! Retourne le nombre de parametres libres
164 inline int GetNFree() const {return mNParFree;};
165 //! Retourne le nombre de parametres bornes
166 inline int GetNBound() const {return mNParBound;};
167 //! Retourne le nstop de convergence
168 inline int GetNStop() const {return nStop;};
169 //! Retourne le nstop de convergence lente.
170 inline int GetNStopLent() const {return nStopLent;};
171 //! Retourne le pointeur sur la GeneralFunction utilisee.
172 inline GeneralFunction* GetFunction() const {return mFunction;};
173 //! Retourne le pointeur sur la GeneralFitData utilisee.
174 inline GeneralFitData* GetGData() const {return mData;};
175
176 void PrintStatus();
177 void PrintFit();
178 void PrintParm(int n);
179 void PrintParm();
180
181 int Fit();
182 double ReCalChi2(int& nddl, double* par = NULL);
183 GeneralFitData DataResidus(bool clean=true);
184 GeneralFitData DataFunction(bool clean=true);
185 void PrintFitErr(int rc);
186
187protected:
188 int mNtry; //!< numero d'appel de la routine de fit.
189 int mNVar; //!< nombre de variables f(x,y,z,...)
190 int mNPar; //!< nombre de parametres
191 int mNParFree; //!< nombre de parametres libres
192 int mNParBound; //!< nombre de parametres bornes
193 GeneralFunction* mFunction;
194 GeneralXi2* mFuncXi2;
195 GeneralFitData* mData;
196
197 TVector<r_8> Param;
198 TVector<r_8> errParam;
199 TVector<r_8> stepParam;
200 TVector<r_8> minParam;
201 TVector<r_8> maxParam;
202 TVector<r_8> minStepDeriv;
203 TVector<r_8> Eps;
204 unsigned short int* fixParam;
205 unsigned short int* boundParam;
206 string* nameParam;
207
208 double Lambda_Fac;
209 double stopChi2;
210 int maxStep;
211 int nStopMx;
212 double stopChi2SMx;
213 int nStopLent;
214 int debugLevel;
215 FILE *FileStep;
216
217 TMatrix<r_8> ATGA;
218 TVector<r_8> BETA;
219 TMatrix<r_8> ATGA_Try;
220 TVector<r_8> BETA_Try;
221 TVector<r_8> C;
222 TVector<r_8> D;
223
224 double Chi2;
225 int mNddl;
226 int nStep;
227 int nStop, nStopL;
228 double Lambda;
229
230 // Fonctions privees
231 void write_in_step(double ci2,TVector<r_8>& par);
232 void General_Init(void);
233 void TryFunc(TVector<r_8>& par,TVector<r_8>& par_tr);
234 void TryXi2(TVector<r_8>& par,TVector<r_8>& par_tr);
235 void CheckSanity();
236 void Set_Bound_C_D(int i);
237 void Set_Bound_C_D();
238 double p_vers_tr(int i,double p);
239 TVector<r_8> p_vers_tr(TVector<r_8> const& p);
240 void p_vers_tr(TVector<r_8> const& p,TVector<r_8>& tr);
241 double tr_vers_p(int i,double tr);
242 TVector<r_8> tr_vers_p(TVector<r_8> const& tr);
243 void tr_vers_p(TVector<r_8> const& tr,TVector<r_8>& p);
244 double c_dp_vers_dtr(int i,double tr);
245 TVector<r_8> dp_vers_dtr(TVector<r_8> const& dp,TVector<r_8> const& tr);
246 void dp_vers_dtr(TVector<r_8> const& dp,TVector<r_8> const& tr,TVector<r_8>& dtr);
247 double c_dtr_vers_dp(int i,double tr);
248 TVector<r_8> dtr_vers_dp(TVector<r_8> const& dtr,TVector<r_8> const& tr);
249 int put_in_limits_for_deriv(TVector<r_8> const& p,TVector<r_8>& dp,double dist=0.66);
250 inline void dtr_vers_dp(TVector<r_8> const& dtr,TVector<r_8> const& tr,TVector<r_8>& dp)
251 { for(int i=0;i<mNPar;i++)
252 { if( fixParam[i] ) continue;
253 if( ! boundParam[i] ) dp(i) = dtr(i);
254 else dp(i) = D(i)/(1.+tr(i)*tr(i)) * dtr(i); }
255 };
256};
257
258} // Fin du namespace
259
260#endif
Note: See TracBrowser for help on using the repository browser.