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

Last change on this file since 1516 was 938, checked in by ansari, 26 years ago

Vector,Matrix -> TVector,TMatrix<r_8> cmv 14/4/00
File size: 9.3 KB
RevLine 
[220]1// This may look like C code, but it is really -*- C++ -*-
2#ifndef GENERALFIT_SEEN
3#define GENERALFIT_SEEN
4
[244]5#include "pexceptions.h"
[514]6#include "tvector.h"
[220]7#include "generaldata.h"
8
[552]9namespace SOPHYA {
[307]10
[220]11//================================================================
12// GeneralFunction
13//================================================================
14
[926]15//! Classe de fonctions parametrees a plusieurs variables
[307]16class GeneralFunction {
[220]17public:
18 GeneralFunction(unsigned int nVar, unsigned int nPar);
19 virtual ~GeneralFunction();
20
[914]21 //! Valeur de la fonction a definir par l'utilisateur (fct virtuelle pure)
[220]22 virtual double Value(double const xp[], double const* parm)=0;
[926]23 //! Valeur de la fonction derivee selon les parametres pouvant etre redefinie
[220]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
[914]30 //! Retourne le nombre de variables Xi
[220]31 inline int NVar() const {return mNVar;}
[914]32 //! Retourne le nombre de parametres Ai
[220]33 inline int NPar() const {return mNPar;}
34
35protected:
[914]36 const int mNVar; //!< nombre de variables f(x,y,z,...)
37 const int mNPar; //!< nombre de parametres
[220]38
39 double *deltaParm;
40 double *tmpParm;
41};
42
43//================================================================
44// GeneralFunc
45//================================================================
46
[926]47//! Classe de fonctions parametrees a plusieurs variables type C
[220]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
[926]70//! Classe de Xi2 a plusieurs parametres
[307]71class GeneralXi2 {
[220]72public:
73 GeneralXi2(unsigned int nPar);
74 virtual ~GeneralXi2();
75
[914]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 */
[220]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
[914]89 //! Retourne le nombre de parametres Ai.
[220]90 inline int NPar() const {return mNPar;}
91
92protected:
[914]93 const int mNPar; //!< nombre de parametres
[220]94
95 double *deltaParm;
96};
97
98//================================================================
99// GENERALFIT
100//================================================================
101
[914]102//! Classe de fit d'une GeneralFunction sur une GeneralFitData
[307]103class GeneralFit {
[220]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);
[938]140 TVector<r_8> GetParm();
[220]141 double GetParmErr(int);
142 double GetCoVar(int,int);
143 double GetStep(int n);
144 double GetMax(int n);
145 double GetMin(int n);
[914]146 //! Retourne le Chi2
[220]147 inline double GetChi2() const {return Chi2;};
[914]148 //! Retourne le Chi2 reduit
[220]149 inline double GetChi2Red() const {
150 if(mNddl<=0) return (double) mNddl;
151 return Chi2/(double) mNddl;
152 };
[914]153 //! Retourne la precision de convergence pour le parametre i.
[220]154 inline double GetEps(int i) const {return Eps(i);};
[914]155 //! Retourne le nombre de degres de liberte
[220]156 inline int GetNddl() const {return mNddl;};
[914]157 //! Retourne le nombre d'iterations
[220]158 inline int GetNStep() const {return nStep;};
[914]159 //! Retourne le nombre de variables
[220]160 inline int GetNVar() const {return mNVar;};
[914]161 //! Retourne le nombre de parametres
[220]162 inline int GetNPar() const {return mNPar;};
[914]163 //! Retourne le nombre de parametres libres
[220]164 inline int GetNFree() const {return mNParFree;};
[914]165 //! Retourne le nombre de parametres bornes
[220]166 inline int GetNBound() const {return mNParBound;};
[914]167 //! Retourne le nstop de convergence
[220]168 inline int GetNStop() const {return nStop;};
[914]169 //! Retourne le nstop de convergence lente.
[220]170 inline int GetNStopLent() const {return nStopLent;};
[914]171 //! Retourne le pointeur sur la GeneralFunction utilisee.
[220]172 inline GeneralFunction* GetFunction() const {return mFunction;};
[914]173 //! Retourne le pointeur sur la GeneralFitData utilisee.
[220]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);
[307]183 GeneralFitData DataResidus(bool clean=true);
184 GeneralFitData DataFunction(bool clean=true);
[220]185 void PrintFitErr(int rc);
186
187protected:
[914]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
[220]193 GeneralFunction* mFunction;
194 GeneralXi2* mFuncXi2;
195 GeneralFitData* mData;
196
[938]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;
[220]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
[938]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;
[220]223
224 double Chi2;
225 int mNddl;
226 int nStep;
227 int nStop, nStopL;
228 double Lambda;
229
230 // Fonctions privees
[938]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)
[220]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
[307]258} // Fin du namespace
259
[220]260#endif
Note: See TracBrowser for help on using the repository browser.