source: Sophya/trunk/SophyaPI/PIext/piafitting.cc@ 375

Last change on this file since 375 was 368, checked in by ercodmgr, 26 years ago

modifs Fenetre de fit et gestion repertoire /func Reza 9/8/99
File size: 25.1 KB
RevLine 
[361]1#include <stdio.h>
2#include <stdlib.h>
3#include <ctype.h>
4
5#include <iostream.h>
6
7#include <typeinfo>
8
9#include "piafitting.h"
10#include "strutil.h"
11
12#include "generalfit.h"
13
14#include "fct1dfit.h"
15#include "fct2dfit.h"
16
17#include "matrix.h"
18#include "cvector.h"
19#include "ntuple.h"
20#include "cimage.h"
21
22#include "histos.h"
23#include "histos2.h"
24#include "ntuple.h"
25#include "hisprof.h"
26
27#include "nobjmgr.h"
28#include "pistdimgapp.h"
29
30//// ---------- Classe de fenetre de fit interactive ------------
[368]31#define FINPARMAX 10 // Nb maxi de parametres de fit
32
[361]33class PIAFitterWind : public PIWindow {
34public :
35 PIAFitterWind(PIStdImgApp* par, PIAFitter* fiter);
36 virtual ~PIAFitterWind();
37 virtual void Show();
38 virtual void Process(PIMessage msg, PIMsgHandler* sender, void* data=NULL);
39 inline void SetObjName(string& nom) { oname = nom; }
[368]40 void SetNbParms(int npar);
[361]41protected :
42 PIStdImgApp* dap;
43 PIAFitter* fitter;
[368]44 PILabel* labn;
[361]45 PIOptMenu* pom[2];
46 PIButton* but[2];
47 PICheckBox* ckb[2];
[368]48 PILabel* lab[2];
49 PILabel* labp[FINPARMAX]; // Maximum FINPARMAX parametres
50 PIText* txtp[FINPARMAX]; // " " "
51 PIText* txts[FINPARMAX]; // " " "
[361]52 string oname;
53};
54
55///////////////////// Fit 1D et 2D //////////////////////////
56
57PIAFitterWind* fwind = NULL;
58
59/* --Methode-- */
60PIAFitter::PIAFitter(PIACmd *piac, PIStdImgApp* app)
61{
62string kw, usage;
63
64kw = "fit";
65usage = "Fitting function to DataObjects (Histo, Histo2D, Vector, ...)";
66usage += "\n Usage: fit nomobj func [Options]";
67usage += "\n [p:p1,...,pn s:s1,...,sn m:m1,...,mn M:M1,...,Mn o:... o:...]";
68piac->RegisterCommand(kw, usage, this, "Fitting");
69
70kw = "fitw";
71usage = "Interactive Window Fit";
72usage += "\n Usage: fit nomobj ";
73piac->RegisterCommand(kw, usage, this, "Fitting");
74
75// Je cree une seule copie de PIAFitterWind
76// Attention, ca ne marchera pas si on detruit this
77if (fwind == NULL) fwind = new PIAFitterWind(app, this);
78}
79
80/* --Methode-- */
81PIAFitter::~PIAFitter()
82{
83}
84
85int PIAFitter::Execute(string& kw, vector<string>& tokens)
86{
87// Fit 1D sur objets 1D. Egalement Fit 2D sur objets 2D.
88if (kw == "fit") {
89 if (tokens.size() < 2) {
90 cout <<"Usage:fit nomobj func \n"
91 <<" [p:p1,...,pn s:s1,...,sn m:m1,...,mn M:M1,...,Mn o:... o:...]\n";
92 return(0);
93 }
94 string p=""; string s=""; string m=""; string M=""; string O="";
95 if (tokens.size()>2)
96 for(int ip=2;ip<tokens.size();ip++) {
97 if(tokens[ip].length()<=2) continue;
98 const char *c = tokens[ip].c_str();
99 if(c[1]!=':') continue;
100 if(c[0]=='p') p=c+2;
101 else if(c[0]=='s') s=c+2;
102 else if(c[0]=='m') m=c+2;
103 else if(c[0]=='M') M=c+2;
104 else if(c[0]=='o') {O += ","; O += c+2;}
105 }
106 Fit12D(tokens[0],tokens[1],p,s,m,M,O);
107}
108else if (kw == "fitw") {
[368]109 if (tokens.size() < 1) {
110 cout << "Usage: fitw nomobj" << endl;
111 return(0);
112 }
[361]113 cout << " *DBG* Affichage PIAFitterWind ** " << endl;
[368]114 fwind->SetObjName(tokens[0]);
[361]115 fwind->Show();
116 }
117
118return(0);
119}
120
121/* --Methode-- cmv 13/10/98 */
122void PIAFitter::Fit12D(string& nom, string& func,
123 string par,string step,string min,string max,
124 string opt)
125//| --------------- Fit d'objets a 1 et 2 dimensions ---------------
126//| nom : nom de l'objet qui peut etre:
127//| fit-1D: Vector,Histo1D,HProf ou GeneraFitData(1D)
128//| fit-2D: Matrix,Histo2D,Image<T> ou GeneraFitData(2D)
129//| func : pnn : fit polynome degre nn avec classe Poly (lineaire) 1D ou 2D
130//| : Pnn : fit polynome degre nn avec GeneralFit (non-lineaire) 1D ou 2D
131//| : gnn : fit gaussienne (hauteur) + polynome de degre nn 1D
132//| : g : fit gaussienne (hauteur) 1D
133//| : enn : fit exponentielle + polynome de degre nn 1D
134//| : e : fit exponentielle 1D
135//| : Gnn : fit gaussienne (volume) + polynome de degre nn 1D
136//| : G : fit gaussienne (volume) 1D
137//| : : fit gaussienne+fond (volume) 2D
138//| : Gi : fit gaussienne+fond integree (volume) 2D
139//| : d : fit DL de gaussienne+fond (volume) 2D
140//| : di : fit DL de gaussienne+fond integree (volume) 2D
141//| : D : fit DL de gaussienne+fond avec coeff variable p6 (volume) 2D
142//| : Di : fit DL de gaussienne+fond integree avec coeff variable p6 (volume) 2D
143//| : M : fit Moffat+fond (expos=p6) (volume) 2D
144//| : Mi : fit Moffat+fond integree (expos=p6) (volume) 2D
145//| par : p1,...,pn : valeur d'initialisation des parametres (def=0)
146//| step : s1,...,sn : valeur des steps de depart (def=1)
147//| min : m1,...,mn : valeur des minima (def=1)
148//| max : M1,...,Mn : valeur des maxima (def=-1) (max<=min : pas de limite)
149//| opt : options "Eaa.b,eaa.b,f,r,caa.b,Xaa.b"
150//| f : generation d'un Objet identique contenant la fonction fittee
151//| r : generation d'un Objet identique contenant les residus
152//| Xaa.b : aa.b valeur du DXi2 d'arret (def=1.e-3)
153//| Naa : aa nombre maximum d'iterations (def=100)
154//| la.b : niveau "a.b" de print: a=niveau de print Fit1/2D
155//| b=niveau de debug GeneralFit
156//| Ii1/i2 numeros des bins X de l'histos utilises pour le fit [i1,i2]
157//|2D Jj1/j2 numeros des bins Y de l'histos utilises pour le fit [j1,j2]
158//| - L'erreur est celle associee a l'objet (si elle existe),
159//| elle est mise a 1 sinon, sauf si E... ou e... est precise:
160//| Eaa.b : si |val|>=1 erreur = aa.b*sqrt(|val|)
161//| si |val|<1 erreur = aa.b
162//| si aa.b <=0 alors aa.b=1.0
163//| E seul est equivalent a E1.0
164//| eaa.b : erreur = aa.b
165//| si aa.b <=0 alors aa.b=1.0
166//| e seul est equivalent a e1.0
167//| xaa.b : demande de centrage: on fit x-aa.b au lieu de x)
168//| x : demande de centrage: on fit x-xc au lieu de x
169//| avec xc=abscisse du milieu de l'histogramme
170//| Actif pour exp+poly 1D, poly 1D
171//| pour gauss+poly 1D, xc est le centre de la gaussienne.
172//|2D yaa.b et y : idem "xaa.b et x" mais pour y
173{
174
175NamedObjMgr omg;
176AnyDataObj* obj = omg.GetObj(nom);
177if (obj == NULL) {
178 cout<<"PIAFitter::Fit12D() Error , Pas d'objet de nom "<<nom<<endl;
179 return;
180}
181if(func.length()<=0)
182 {cout<<"PIAFitter::Fit12D() Donnez un nom de fonction a fitter."<<endl;
183 return;}
184string ctyp = typeid(*obj).name();
185
186int ndim = 0, nbinx=0, nbiny=0, ndata = 0;
187Vector* v = NULL; Histo* h = NULL;
188Matrix* m = NULL; Histo2D* h2 = NULL; RzImage* im = NULL;
189GeneralFitData* g = NULL;
190
191 // 1D
192if (typeid(*obj) == typeid(Vector)) {
193 ndim = 1;
194 v = (Vector*) obj; nbinx = v->NElts(); nbiny = 1;
195 }
196else if ( (typeid(*obj) == typeid(HProf)) || (typeid(*obj) == typeid(Histo)) ) {
197 ndim = 1;
198 h = (Histo*) obj; nbinx = h->NBins(); nbiny = 1;
199 }
200else if (typeid(*obj) == typeid(Matrix)) {
201 ndim = 2;
202 m = (Matrix*) obj; nbinx = m->NCol(); nbiny = m->NRows();
203 }
204else if (typeid(*obj) == typeid(Histo2D)) {
205 ndim = 2;
206 h2 = (Histo2D*) obj; nbinx = h2->NBinX(); nbiny = h2->NBinY();
207 }
208else if (typeid(*obj) == typeid(GeneralFitData)) {
209 g = (GeneralFitData*) obj; nbinx = g->NData(); nbiny = 1;
210 if( g->NVar()==1) ndim = 1;
211 else if(g->NVar()==2) ndim = 2;
212 else {
213 cout<<"GeneralFitData ne peut avoir que 1 ou 2 variables d'abscisse: "
214 <<((GeneralFitData*) obj)->NVar()<<endl; return; }
215 }
216else if (dynamic_cast<RzImage*>(obj)) {
217 ndim = 2;
218 im = (RzImage*) obj; nbinx = im->XSize(); nbiny = im->YSize();
219 }
220else {
221 cout<<"PIAFitter::Fit12D() Error , Objet n'est pas un "
222 <<"Histo1D/HProf/Vector/Histo2D/Image/Matrix/GeneralFitData "<<ctyp<<endl;
223 return;
224 }
225
226ndata = nbinx*nbiny;
227if(ndata<=0)
228 {cout<<"L'objet a "<<nbinx<<","<<nbiny<<" bins ("<<ndata<<")"<<endl; return;}
229
230// Decodage des options et des parametres, mise en forme
231Vector Par(1); Vector Step(1); Vector Min(1); Vector Max(1); DFOPTIONS O;
232DecodeFitsOptions(par,step,min,max,opt,Par,Step,Min,Max,O);
233O.i1 = (O.i1<0||O.i1>=nbinx)? 0: O.i1;
234O.i2 = (O.i2<0||O.i2>=nbinx||O.i2<O.i1)? nbinx-1: O.i2;
235if(ndim>=2) {
236 O.j1 = (O.j1<0||O.j1>=nbiny)? 0: O.j1;
237 O.j2 = (O.j2<0||O.j2>=nbiny||O.j2<O.j1)? nbiny-1: O.j2;
238} else O.j2 = O.j1 = 0;
239if(O.polcx==2) {
240 if(v||m) O.xc = (O.i2-O.i1+1)/2.;
241 else if(h) O.xc = (h->XMin()+h->XMax())/2.;
242 else if(h2) O.xc = (h2->XMin()+h2->XMax())/2.;
243 else if(g) {double mini,maxi; g->GetMinMax(2,mini,maxi); O.xc=(mini+maxi)/2.;}
244 else if(im) {O.xc = im->XOrg() * im->XPxSize()*(O.i2-O.i1+1)/2.;}
245}
246if(O.polcy==2 && ndim>=2) {
247 if(m) O.yc = (O.j2-O.j1+1)/2.;
248 if(h2) O.yc = (h2->YMin()+h2->YMax())/2.;
249 if(g) {double mini,maxi; g->GetMinMax(12,mini,maxi); O.yc=(mini+maxi)/2.;}
250 if(im) {O.yc = im->YOrg() * im->YPxSize()*(O.j2-O.j1+1)/2.;}
251}
252if(O.lp>0)
253 cout<<"Fit["<<nbinx<<","<<nbiny<<"] ("<<ndata<<") dim="<<ndim<<":"
254 <<" Int=["<<O.i1<<","<<O.i2<<"],["<<O.j1<<","<<O.j2<<"]"<<endl
255 <<" Cent="<<O.polcx<<","<<O.polcy<<","<<O.xc<<"+x"<<","<<O.yc<<"+y"
256 <<" TypE="<<O.err_e<<","<<O.err_E
257 <<" StpX2="<<O.stc2<<" Nstep="<<O.nstep
258 <<" lp,lpg="<<O.lp<<","<<O.lpg<<endl;
259
260///////////////////////////////////
261// Remplissage de GeneralFitData //
262///////////////////////////////////
263GeneralFitData mydata(ndim,ndata,0);
264{for(int i=O.i1;i<=O.i2;i++) for(int j=O.j1;j<=O.j2;j++) {
265 double x,y,f,e;
266
267 if(v)
268 {x= (double) i; f=(*v)(i); e=1.;}
269 else if(h)
270 {x=h->BinCenter(i); f=(*h)(i); e=(h->HasErrors())?h->Error(i):1.;}
271 else if(m)
272 {x=(double) i; y=(double) j; f=(*m)(j,i); e=1.;}
273 else if(h2)
274 {float xf,yf; h2->BinCenter(i,j,xf,yf); x=(double)xf; y=(double)yf;
275 f=(*h2)(i,j); e=(h2->HasErrors())?h2->Error(i,j):1.;}
276 else if(im)
277 {x=im->XOrg()+(i+0.5)*im->XPxSize(); y=im->YOrg()+(j+0.5)*im->YPxSize();
278 f=im->DValue(i,j); e=1.;}
279 else if(g&&ndim==1) {x= g->X(i); f=g->Val(i); e=g->EVal(i);}
280 else if(g&&ndim==2) {x= g->X(i); y= g->Y(i); f=g->Val(i); e=g->EVal(i);}
281 else x=y=f=e=0.;
282
283 // Gestion des erreurs a utiliser
284 if(O.err_e>0.) e=O.err_e;
285 else if(O.err_E>0.) {e=(f<-1.||f>1.)?O.err_E*sqrt(fabs(f)):O.err_E;}
286
287 // Remplissage de generalfit
288 if(func[0]=='p') {x -= O.xc; if(ndim>=2) y -= O.yc;}
289 if(ndim==1) mydata.AddData1(x,f,e);
290 else if(ndim==2) mydata.AddData2(x,y,f,e);
291}}
292if(mydata.NData()<=0)
293 {cout<<"Pas de donnees dans GeneralFitData: "<<mydata.NData()<<endl;
294 return;}
295if(O.lpg>1) {
296 mydata.PrintStatus();
297 mydata.PrintData(0);
298 mydata.PrintData(mydata.NData()-1);
299}
300
301////////////////////////////////////////////
302// Identification de la fonction a fitter //
303////////////////////////////////////////////
304GeneralFunction* myfunc = NULL;
305if(func[0]=='p' && ndim==1) {
306 // Fit de polynome sans passer par les GeneralFit
307 int degre = 0;
308 if(func.length()>1) sscanf(func.c_str()+1,"%d",&degre);
309 cout<<"Fit (lineaire) 1D polynome de degre "<<degre<<endl;
310 Poly p1(0);
311 double c2rl = mydata.PolFit(0,p1,degre);
312 cout<<"C2r_lineaire = "<<c2rl<<endl;
313 if(O.lp>0) cout<<p1<<endl;
314 return;
315
316} else if(func[0]=='P' && ndim==1) {
317 // Fit de polynome
318 int degre = 0;
319 if(func.length()>1) sscanf(func.c_str()+1,"%d",&degre);
320 cout<<"Fit polynome 1D de degre "<<degre<<endl;
321 Polyn1D* myf = new Polyn1D(degre,O.xc);
322 myfunc = myf;
323
324} else if(func[0]=='e' && ndim==1) {
325 // Fit d'exponentielle
326 int degre =-1;
327 if(func.length()>1) sscanf(func.c_str()+1,"%d",&degre);
328 cout<<"Fit d'exponentielle+polynome 1D de degre "<<degre<<endl;
329 Exp1DPol* myf;
330 if(degre>=0) myf = new Exp1DPol((unsigned int)degre,O.xc);
331 else myf = new Exp1DPol(O.xc);
332 myfunc = myf;
333
334} else if(func[0]=='g' && ndim==1) {
335 // Fit de gaussienne en hauteur
336 int degre =-1;
337 if(func.length()>1) sscanf(func.c_str()+1,"%d",&degre);
338 cout<<"Fit de Gaussienne_en_hauteur+polynome 1D de degre "<<degre<<endl;
339 Gauss1DPol* myf;
340 if(degre>=0) myf = new Gauss1DPol((unsigned int)degre,((O.polcx)?true:false));
341 else { bool bfg = (O.polcx)?true:false; myf = new Gauss1DPol(bfg); }
342 myfunc = myf;
343
344} else if(func[0]=='G' && ndim==1) {
345 // Fit de gaussienne en volume
346 int degre =-1;
347 if(func.length()>1) sscanf(func.c_str()+1,"%d",&degre);
348 cout<<"Fit de Gaussienne_en_volume+polynome 1D de degre "<<degre<<endl;
349 GaussN1DPol* myf;
350 if(degre>=0) myf = new GaussN1DPol((unsigned int)degre,((O.polcx)?true:false));
351 else { bool bfg = (O.polcx)?true:false; myf = new GaussN1DPol(bfg); }
352 myfunc = myf;
353
354} else if(func[0]=='p' && ndim==2) {
355 // Fit de polynome 2D sans passer par les GeneralFit
356 int degre = 0;
357 if(func.length()>1) sscanf(func.c_str()+1,"%d",&degre);
358 cout<<"Fit (lineaire) polynome 2D de degre "<<degre<<endl;
359 Poly2 p2(0);
360 double c2rl = mydata.PolFit(0,1,p2,degre);
361 cout<<"C2r_lineaire = "<<c2rl<<endl;
362 if(O.lp>0) cout<<p2<<endl;
363 return;
364
365} else if(func[0]=='P' && ndim==2) {
366 // Fit de polynome 2D
367 int degre = 0;
368 if(func.length()>1) sscanf(func.c_str()+1,"%d",&degre);
369 cout<<"Fit polynome 2D de degre "<<degre<<endl;
370 Polyn2D* myf = new Polyn2D(degre,O.xc,O.yc);
371 myfunc = myf;
372
373} else if(func[0]=='G' && ndim==2) {
374 // Fit de gaussienne+fond en volume
375 int integ = 0;
376 if(func.length()>1) if(func[1]=='i') integ=1;
377 cout<<"Fit de Gaussienne+Fond 2D integ="<<integ<<endl;
378 if(integ) {GauRhInt2D* myf = new GauRhInt2D; myfunc = myf;}
379 else {GauRho2D* myf = new GauRho2D; myfunc = myf;}
380
381} else if(func[0]=='d' && ndim==2) {
382 // Fit de DL gaussienne+fond en volume
383 int integ = 0;
384 if(func.length()>1) if(func[1]=='i') integ=1;
385 cout<<"Fit de DL de Gaussienne+Fond 2D integ="<<integ<<endl;
386 if(integ) {GdlRhInt2D* myf = new GdlRhInt2D; myfunc = myf;}
387 else {GdlRho2D* myf = new GdlRho2D; myfunc = myf;}
388
389} else if(func[0]=='D' && ndim==2) {
390 // Fit de DL gaussienne+fond avec coeff variable p6 en volume
391 int integ = 0;
392 if(func.length()>1) if(func[1]=='i') integ=1;
393 cout<<"Fit de DL de Gaussienne+Fond avec coeff variable (p6) 2D integ="<<integ<<endl;
394 if(integ) {Gdl1RhInt2D* myf = new Gdl1RhInt2D; myfunc = myf;}
395 else {Gdl1Rho2D* myf = new Gdl1Rho2D; myfunc = myf;}
396
397} else if(func[0]=='M' && ndim==2) {
398 // Fit de Moffat+fond (volume)
399 int integ = 0;
400 if(func.length()>1) if(func[1]=='i') integ=1;
401 cout<<"Fit de Moffat+Fond (expos=p6) 2D integ="<<integ<<endl;
402 if(integ) {MofRhInt2D* myf = new MofRhInt2D; myfunc = myf;}
403 else {MofRho2D* myf = new MofRho2D; myfunc = myf;}
404
405} else {
406 cout<<"Fonction "<<func<<" inconnue pour la dim "<<ndim<<endl;
407 return;
408}
409
410/////////////////////////
411// Fit avec generalfit //
412/////////////////////////
413if(myfunc->NPar()>Par.NElts())
414 {cout<<"Trop de parametres: "<<myfunc->NPar()<<">"<<Par.NElts()<<endl;
415 if(myfunc) delete myfunc; return;}
416GeneralFit myfit(myfunc);
417myfit.SetDebug(O.lpg);
418myfit.SetData(&mydata);
419myfit.SetStopChi2(O.stc2);
420myfit.SetMaxStep(O.nstep);
421{for(int i=0;i<myfunc->NPar();i++) {
422 char str[10];
423 sprintf(str,"P%d",i);
424 myfit.SetParam(i,str,Par(i),Step(i),Min(i),Max(i));
425}}
426if(O.lp>1) myfit.PrintFit();
427double c2r = -1.;
428int rcfit = (double) myfit.Fit();
429if(O.lp>0) myfit.PrintFit();
430if(rcfit>0) {
431 c2r = myfit.GetChi2Red();
432 cout<<"C2r_Reduit = "<<c2r<<" nstep="<<myfit.GetNStep()<<" rc="<<rcfit<<endl;
433 Vector ParFit(myfunc->NPar());
434 for(int i=0;i<myfunc->NPar();i++) ParFit(i)=myfit.GetParm(i);
435} else {
436 cout<<"echec Fit, rc = "<<rcfit<<" nstep="<<myfit.GetNStep()<<endl;
437 myfit.PrintFitErr(rcfit);
438}
439
440// Mise a disposition des resultats
441if(rcfit>=0 && myfunc && (O.okres>0||O.okfun>0)) {
442 string nomres = nom + "res";
443 string nomfun = nom + "fun";
444 if(v) {
445 if(O.okres) {Vector* ob = v->FitResidus(myfit); if(ob) omg.AddObj(ob,nomres);}
446 if(O.okfun) {Vector* ob = v->FitFunction(myfit); if(ob) omg.AddObj(ob,nomfun);}
447 } else if(h) {
448 if(O.okres) {Histo* ob = h->FitResidus(myfit); if(ob) omg.AddObj(ob,nomres);}
449 if(O.okfun) {Histo* ob = h->FitFunction(myfit); if(ob) omg.AddObj(ob,nomfun);}
450 } else if(m) {
451 if(O.okres) {Matrix* ob = m->FitResidus(myfit); if(ob) omg.AddObj(ob,nomres);}
452 if(O.okfun) {Matrix* ob = m->FitFunction(myfit); if(ob) omg.AddObj(ob,nomfun);}
453 } else if(h2) {
454 if(O.okres) {Histo2D* ob = h2->FitResidus(myfit); if(ob) omg.AddObj(ob,nomres);}
455 if(O.okfun) {Histo2D* ob = h2->FitFunction(myfit); if(ob) omg.AddObj(ob,nomfun);}
456 } else if(im) {
457 if(O.okres) {RzImage* ob = im->FitResidus(myfit); if(ob) omg.AddObj(ob,nomres);}
458 if(O.okfun) {RzImage* ob = im->FitFunction(myfit); if(ob) omg.AddObj(ob,nomfun);}
459 } else if(g) {
460 if(O.okres) {GeneralFitData* ob = g->FitResidus(myfit); if(ob) omg.AddObj(ob,nomres);}
461 if(O.okfun) {GeneralFitData* ob = g->FitFunction(myfit); if(ob) omg.AddObj(ob,nomfun);}
462 }
463}
464
465// Nettoyage
466if(myfunc) delete myfunc;
467return;
468}
469
470
471/* --Function static propre aux routines de fit 1D et 2D-- cmv 13/10/98 */
472void PIAFitter::DecodeFitsOptions(string par,string step,string min,string max,string opt
473 ,Vector& Par,Vector& Step,Vector& Min,Vector& Max,DFOPTIONS& O)
474//| Pour decoder les "string" et remplir les vecteurs du fit (cf commentaires dans Fit1D)
475{
476// set des vecteurs et decodage des string correspondantes
477int NParMax = 100;
478Par.Realloc(NParMax); Step.Realloc(NParMax);
479Min.Realloc(NParMax); Max.Realloc(NParMax);
480{
481 Vector* v=NULL; string* s=NULL;
482 {for(int i=0;i<NParMax;i++) {Par(i)=0.; Step(i)=1.; Min(i)=1.; Max(i)=-1.;}}
483 for(int j=0;j<4;j++) {
484 if(j==0) {v=&Par; s=&par;}
485 else if(j==1) {v=&Step; s=&step;}
486 else if(j==2) {v=&Min; s=&min;}
487 else if(j==3) {v=&Max; s=&max;}
488 if(s->length()>0) *s += ",";
489 for(int i=0;i<NParMax;i++) {
490 if(s->length()<=0) break;
491 sscanf(s->c_str(),"%lf",&(*v)(i));
492 size_t p = s->find_first_of(',') + 1;
493 if(p>=s->length()) *s = ""; else *s = s->substr(p);
494 }
495 }
496}
497
498// Decodage de options de opt
499O.okres = O.okfun = 0;
500O.polcx = O.polcy = 0;
501O.xc = O.yc = 0.;
502O.stc2 = 1.e-3;
503O.nstep = 100;
504O.err_e = O.err_E = -1.;
505O.lp = 1; O.lpg = 0;
506O.i1 = O.j1 = O.i2 = O.j2 = -1;
507
508if(opt.length()<=0) return;
509opt = "," + opt + ",";
510
511if(strstr(opt.c_str(),",r,")) O.okres = 1; // residus
512if(strstr(opt.c_str(),",f,")) O.okfun = 1; // fonction fittee
513if(strstr(opt.c_str(),",x")) { // Demande de centrage (fit X=x-xc)
514 O.polcx = 2; // Le centrage est calcule automatiquement
515 size_t p = opt.find(",x");
516 size_t q = opt.find_first_of(',',p+1);
517 string dum = opt.substr(p,q-p);
518 if(dum.length()>2) {
519 sscanf(dum.c_str(),",x%lf",&O.xc);
520 O.polcx = 1; // Le centrage est fixe par la valeur lue
521 }
522}
523if(strstr(opt.c_str(),",y")) { // Demande de centrage (fit Y=y-yc)
524 O.polcy = 2; // Le centrage est calcule automatiquement
525 size_t p = opt.find(",y");
526 size_t q = opt.find_first_of(',',p+1);
527 string dum = opt.substr(p,q-p);
528 if(dum.length()>2) {
529 sscanf(dum.c_str(),",y%lf",&O.yc);
530 O.polcy = 1; // Le centrage est fixe par la valeur lue
531 }
532}
533if(strstr(opt.c_str(),",E")) { // Erreurs imposees a "sqrt(val)" ou "aa.b*sqrt(val)"
534 size_t p = opt.find(",E");
535 size_t q = opt.find_first_of(',',p+1);
536 string dum = opt.substr(p,q-p);
537 if(dum.length()>2) sscanf(dum.c_str(),",E%lf",&O.err_E);
538 if(O.err_E<=0.) O.err_E = 1.;
539 O.err_e=-1.;
540}
541if(strstr(opt.c_str(),",e")) { // Erreurs imposees a "1" ou "aa.b"
542 size_t p = opt.find(",e");
543 size_t q = opt.find_first_of(',',p+1);
544 string dum = opt.substr(p,q-p);
545 if(dum.length()>2) sscanf(dum.c_str(),",e%lf",&O.err_e);
546 if(O.err_e<=0.) O.err_e = 1.;
547 O.err_E=-1.;
548}
549if(strstr(opt.c_str(),",X")) { // Valeur du StopChi2
550 size_t p = opt.find(",X");
551 size_t q = opt.find_first_of(',',p+1);
552 string dum = opt.substr(p,q-p);
553 if(dum.length()>2) sscanf(dum.c_str(),",X%lf",&O.stc2);
554 if(O.stc2<=0.) O.stc2 = 1.e-3;
555}
556if(strstr(opt.c_str(),",N")) { // Nombre maximum d'iterations
557 size_t p = opt.find(",N");
558 size_t q = opt.find_first_of(',',p+1);
559 string dum = opt.substr(p,q-p);
560 if(dum.length()>2) sscanf(dum.c_str(),",N%d",&O.nstep);
561 if(O.nstep<2) O.nstep = 100;
562}
563if(strstr(opt.c_str(),",l")) { // niveau de print
564 size_t p = opt.find(",l");
565 size_t q = opt.find_first_of(',',p+1);
566 string dum = opt.substr(p,q-p);
567 float ab;
568 if(dum.length()>2) sscanf(dum.c_str(),",l%f",&ab);
569 if(ab<0) ab = 0.;
570 O.lp = (int) ab; O.lpg = int(10.*(ab-(float)O.lp+0.01));
571}
572if(strstr(opt.c_str(),",I")) { // intervalle de fit selon X
573 size_t p = opt.find(",I");
574 size_t q = opt.find_first_of(',',p+1);
575 string dum = opt.substr(p,q-p);
576 if(dum.length()>2) sscanf(dum.c_str(),",I%d/%d",&O.i1,&O.i2);
577}
578if(strstr(opt.c_str(),",J")) { // intervalle de fit selon Y
579 size_t p = opt.find(",J");
580 size_t q = opt.find_first_of(',',p+1);
581 string dum = opt.substr(p,q-p);
582 if(dum.length()>2) sscanf(dum.c_str(),",J%d/%d",&O.j1,&O.j2);
583}
584return;
585}
586
587
588// ---------------- Fenetre de fit interactive -------
589
590
591/* --Methode-- */
592PIAFitterWind::PIAFitterWind(PIStdImgApp* par, PIAFitter* fiter)
593 : PIWindow((PIMsgHandler*)par, "PIAFitter", PIWK_normal, 240, 240, 150, 150)
594{
595dap = par;
596fitter = fiter;
597int bsx, bsy, spx, spy;
598
599// On definit la taille a partir de la taille par defaut des composantes
600// PIApplicationPrefCompSize(bsx, bsy);
601par->PrefCompSz(bsx, bsy);
602spx = bsx/10;
603spy = bsy/4;
604
[368]605int wszx = 3*spx+4*bsx;
[361]606int wszy = 5*bsy+6*spy;
607SetSize(wszx, wszy);
608int cpx = spx;
609int cpy = spy;
610
[368]611cpx = spx+bsx;
612labn = new PILabel(this, "ObjectName", 2.0*bsx, bsy, cpx, cpy);
613labn->SetBinding(PIBK_elastic, PIBK_fixed, PIBK_elastic, PIBK_free);
[361]614
615cpy += spy+bsy;
[368]616cpx = spx;
617ckb[0] = new PICheckBox(this,"Gen.Func", 1001, 2*bsx, bsy, cpx, cpy);
618cpx += bsx*2+spx;
619ckb[1] = new PICheckBox(this,"Gen.Resid", 1002, 2*bsx, bsy, cpx, cpy);
620ckb[0]->SetBinding(PIBK_elastic, PIBK_fixed, PIBK_elastic, PIBK_free);
621ckb[1]->SetBinding(PIBK_elastic, PIBK_fixed, PIBK_elastic, PIBK_free);
[361]622
623cpy += spy+bsy;
[368]624cpx = spx;
625pom[0] = new PIOptMenu(this, "FitFunc", 2*bsx, bsy, cpx, cpy);
[361]626pom[0]->AppendItem("Poly0", 100);
627pom[0]->AppendItem("Poly1", 101);
628pom[0]->AppendItem("Poly2", 102);
629pom[0]->AppendItem("Gaussienne", 103);
630pom[0]->AppendItem("Gauss+Poly0",104);
631pom[0]->AppendItem("Gauss+Poly1",105);
632pom[0]->AppendItem("Gauss+Poly2",106);
[368]633pom[0]->SetBinding(PIBK_elastic, PIBK_fixed, PIBK_elastic, PIBK_free);
[361]634
[368]635cpx += bsx*2+spx;
636pom[1] = new PIOptMenu(this, "NPar", 2*bsx, bsy, cpx, cpy);
637pom[1]->AppendItem("Un NP=1", 201);
638pom[1]->AppendItem("Deux NP=2",202);
639pom[1]->AppendItem("Trois NP=3",203);
640pom[1]->AppendItem("Cinq NP=5",205);
641pom[1]->AppendItem("Huit NP=8",208);
642pom[1]->AppendItem("Dix NP=10",210);
643pom[1]->SetBinding(PIBK_elastic, PIBK_fixed, PIBK_elastic, PIBK_free);
[361]644
645cpy += spy+bsy;
[368]646cpx = 0.5*bsx;
647but[0] = new PIButton(this, "DoFit", 555, bsx*1.5, bsy, cpx, cpy);
648cpx += bsx*1.5+3*spx;
649but[1] = new PIButton(this, "Dismiss", 777, bsx*1.5, bsy, cpx, cpy);
650but[0]->SetBinding(PIBK_elastic, PIBK_fixed, PIBK_elastic, PIBK_free);
651but[1]->SetBinding(PIBK_elastic, PIBK_fixed, PIBK_elastic, PIBK_free);
[361]652
[368]653// Pour les parametres
654cpy += bsy+spy;
655cpx = 0.6*bsx+1.5*spx;
656lab[0] = new PILabel(this, "Init.Val", 1.7*bsx, bsy, cpx, cpy);
657cpx += 1.7*bsx+0.5*spx;
658lab[1] = new PILabel(this, "Step", 1.7*bsx, bsy, cpx, cpy);
659lab[0]->SetBinding(PIBK_elastic, PIBK_fixed, PIBK_elastic, PIBK_free);
660lab[1]->SetBinding(PIBK_elastic, PIBK_fixed, PIBK_elastic, PIBK_free);
661
662for(int i=0; i<FINPARMAX; i++) {
663 char buff[32];
664 cpy += bsy+spy;
665 cpx = spx;
666 sprintf(buff, "P%d",i);
667 labp[i] = new PILabel(this, buff, 0.6*bsx, bsy, cpx, cpy);
668 labp[i]->SetBinding(PIBK_elastic, PIBK_fixed, PIBK_elastic, PIBK_free);
669 sprintf(buff, "Val%d",i);
670 cpx += 0.6*bsx+0.5*spx;
671 txtp[i] = new PIText(this, buff, 1.7*bsx, bsy, cpx, cpy);
672 txtp[i]->SetBinding(PIBK_elastic, PIBK_fixed, PIBK_elastic, PIBK_free);
673 sprintf(buff, "Step%d",i);
674 cpx += 1.7*bsx+0.5*spx;
675 txts[i] = new PIText(this, buff, 1.7*bsx, bsy, cpx, cpy);
676 txts[i]->SetBinding(PIBK_elastic, PIBK_fixed, PIBK_elastic, PIBK_free);
677 }
[361]678}
[368]679
[361]680/* --Methode-- */
681PIAFitterWind::~PIAFitterWind()
682{
[368]683int i;
684for(i=0; i<2; i++) {
[361]685 delete pom[i];
686 delete ckb[i];
687 delete but[i];
688 }
[368]689delete labn;
690delete lab[0];
691delete lab[1];
692for(i=0; i< FINPARMAX; i++) {
693 delete labp[i];
694 delete txtp[i];
695 delete txts[i];
696 }
[361]697}
[368]698
[361]699/* --Methode-- */
[368]700void PIAFitterWind::SetNbParms(int npar)
701{
702if (npar < 1) npar = 1;
703if (npar > FINPARMAX) npar = FINPARMAX;
704
705int bsx, bsy, spx, spy;
706
707dap->PrefCompSz(bsx, bsy);
708spx = bsx/10;
709spy = bsy/4;
710
711int wszx = 3*spx+4*bsx;
712int wszy = (5+npar)*bsy+(6+npar)*spy;
713SetSize(wszx, wszy);
714
715int i;
716for(i=0; i<npar; i++) {
717 labp[i]->Manage();
718 txtp[i]->Manage();
719 txts[i]->Manage();
720 }
721for(i=npar; i<FINPARMAX; i++) {
722 labp[i]->UnManage();
723 txtp[i]->UnManage();
724 txts[i]->UnManage();
725 }
726
727}
728
729/* --Methode-- */
[361]730void PIAFitterWind::Show()
731{
732// Si on veut initialiser des trucs au moment ou apparait la fenetre
[368]733labn->SetLabel(oname);
734SetNbParms(3);
735pom[1]->SetValue(203);
[361]736PIWindow::Show();
737}
738
739/* --Methode-- */
740void PIAFitterWind::Process(PIMessage msg, PIMsgHandler* sender, void* data)
741{
742
743char *mf[7] = {"p0","p1","p2","g","g0","g1","g2"};
744msg = UserMsg(msg);
[368]745
746if ((msg > 200) && (msg <= 200+FINPARMAX)) { // On change le nombre de parametre
747 SetNbParms(msg-200);
748 return;
749 }
750
[361]751if (msg == 777) { Hide(); return; } // On cache la fenetre
752
753if (msg == 555) { // On fait le fit
754 vector<string> args;
755 string cmd = oname;
756 string a;
757 args.push_back(oname);
758 int imf = pom[0]->GetValue()-100;
759 if ((imf < 0) || (imf > 6)) imf = 0;
760 a = mf[imf];
761 cmd += " " + a + " ";
762 args.push_back(a);
763 cmd += a;
764 args.push_back(a);
765 if (ckb[0]->GetState() || ckb[1]->GetState()) {
766 a = "o:";
767 if (ckb[0]->GetState()) a += "f";
768 if (ckb[1]->GetState()) a += "r";
769 cmd += " " + a;
770 args.push_back(a);
771 }
772 cout << "FitComm= fit " << cmd << endl;
773 string kw = "fit";
774 fitter->Execute(kw, args);
775 }
776
777}
Note: See TracBrowser for help on using the repository browser.