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 |
|
---|
9 | namespace SOPHYA {
|
---|
10 |
|
---|
11 | //================================================================
|
---|
12 | // GeneralFunction
|
---|
13 | //================================================================
|
---|
14 |
|
---|
15 | //! Classe de fonctions parametrees a plusieurs variables
|
---|
16 | class GeneralFunction {
|
---|
17 | public:
|
---|
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 |
|
---|
35 | protected:
|
---|
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
|
---|
48 | class GeneralFunc : public GeneralFunction {
|
---|
49 | public:
|
---|
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 |
|
---|
59 | protected:
|
---|
60 | double (*tmpFun) (double const*, double const*);
|
---|
61 | double (*tmpFunDer)(double const*, double const*, double*);
|
---|
62 | };
|
---|
63 |
|
---|
64 | //================================================================
|
---|
65 | // GeneralXi2
|
---|
66 | //================================================================
|
---|
67 |
|
---|
68 | class GeneralFitData;
|
---|
69 |
|
---|
70 | //! Classe de Xi2 a plusieurs parametres
|
---|
71 | class GeneralXi2 {
|
---|
72 | public:
|
---|
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 |
|
---|
92 | protected:
|
---|
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
|
---|
103 | class GeneralFit {
|
---|
104 | public:
|
---|
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 | Vector 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 |
|
---|
187 | protected:
|
---|
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 | Vector Param;
|
---|
198 | Vector errParam;
|
---|
199 | Vector stepParam;
|
---|
200 | Vector minParam;
|
---|
201 | Vector maxParam;
|
---|
202 | Vector minStepDeriv;
|
---|
203 | Vector 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 | Matrix ATGA;
|
---|
218 | Vector BETA;
|
---|
219 | Matrix ATGA_Try;
|
---|
220 | Vector BETA_Try;
|
---|
221 | Vector C;
|
---|
222 | Vector 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,Vector& par);
|
---|
232 | void General_Init(void);
|
---|
233 | void TryFunc(Vector& par,Vector& par_tr);
|
---|
234 | void TryXi2(Vector& par,Vector& 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 | Vector p_vers_tr(Vector const& p);
|
---|
240 | void p_vers_tr(Vector const& p,Vector& tr);
|
---|
241 | double tr_vers_p(int i,double tr);
|
---|
242 | Vector tr_vers_p(Vector const& tr);
|
---|
243 | void tr_vers_p(Vector const& tr,Vector& p);
|
---|
244 | double c_dp_vers_dtr(int i,double tr);
|
---|
245 | Vector dp_vers_dtr(Vector const& dp,Vector const& tr);
|
---|
246 | void dp_vers_dtr(Vector const& dp,Vector const& tr,Vector& dtr);
|
---|
247 | double c_dtr_vers_dp(int i,double tr);
|
---|
248 | Vector dtr_vers_dp(Vector const& dtr,Vector const& tr);
|
---|
249 | int put_in_limits_for_deriv(Vector const& p,Vector& dp,double dist=0.66);
|
---|
250 | inline void dtr_vers_dp(Vector const& dtr,Vector const& tr,Vector& 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
|
---|