source: Sophya/trunk/SophyaLib/NTools/simplex.cc@ 2650

Last change on this file since 2650 was 2650, checked in by ansari, 21 years ago

Ajout classe MinZSimplex (minimisation simplex, fichier simplex.h .cc) et MAJ Makefile - Reza 09/02/2005

File size: 11.2 KB
Line 
1#include "sopnamsp.h"
2#include "simplex.h"
3#include "ntuple.h"
4#include <math.h>
5
6#include "timing.h"
7
8//---------------------------------------------------------------
9//------------------- Classe MinZFunction -------------------
10//---------------------------------------------------------------
11// Interface de classe de function multivariable pour le SimplexMinmizer
12
13MinZFunction::MinZFunction(unsigned int nvar)
14 : mNVar(nvar)
15{
16}
17
18MinZFunction::~MinZFunction()
19{
20}
21
22//---------------------------------------------------------------
23//------------------- Classe MinZFuncXi2 --------------------
24//---------------------------------------------------------------
25MinZFuncXi2::MinZFuncXi2(GeneralXi2* gxi2, GeneralFitData* gd)
26 : mGXi2(gxi2) , mGData(gd), MinZFunction(gxi2->NPar())
27{
28}
29
30MinZFuncXi2::~MinZFuncXi2()
31{
32}
33
34double MinZFuncXi2::Value(double const xp[])
35{
36 int ndataused;
37 return mGXi2->Value(*mGData, const_cast<double *>(xp), ndataused);
38}
39
40//---------------------------------------------------------------
41//------------------- Classe MinZTestFunc -------------------
42//---------------------------------------------------------------
43class MinZTestFunc : public MinZFunction {
44public:
45 MinZTestFunc(int sel);
46 virtual double Value(double const xp[]);
47 string ToString();
48 Vector OptParms();
49protected:
50 static int ISelToNvar(int isel);
51 int mSel;
52};
53
54int MinZTestFunc::ISelToNvar(int isel)
55{
56 if (isel == 0) return 1;
57 if (isel == 1) return 1;
58 else if (isel == 2) return 1;
59 else if (isel == 3) return 2;
60 else if (isel == 4) return 3;
61 else return 1;
62}
63
64MinZTestFunc::MinZTestFunc(int sel)
65 : MinZFunction(ISelToNvar(sel))
66{
67 if ((sel < 0) || (sel > 4)) sel = 0;
68 mSel = sel;
69}
70
71string MinZTestFunc::ToString()
72{
73 string rs;
74 if (mSel == 0) {
75 rs = "-x+(x-2)^2";
76 }
77 else if (mSel == 1) {
78 rs = "0.1*x^2-3exp(-(x-2)^2)-5*exp(-0.5*(x+3)^2)";
79 }
80 else if (mSel == 2) {
81 rs = "0.1*x^2-3exp(-(x-2)^2)+5*exp(-0.5*(x+3)^2)";
82 }
83 else if (mSel == 3) {
84 rs = "1.3*(x-50.35)^2+25*(y+3.14)^2";
85 }
86 else if (mSel == 4) {
87 rs = "(x-2.2)^2+2.*(y+3.6)^2+3.*(z-1.1)^2";
88 }
89 else rs = "????";
90 return rs;
91}
92
93Vector MinZTestFunc::OptParms()
94{
95 Vector xx;
96 if (mSel == 0) {
97 Vector rv(1);
98 rv = 2.5;
99 return rv;
100 }
101 else if (mSel == 1) {
102 Vector rv(1);
103 rv = -2.883;
104 return rv;
105 }
106 else if (mSel == 2) {
107 Vector rv(1);
108 rv = 1.812;
109 return rv;
110 }
111 else if (mSel == 3) {
112 Vector rv(2);
113 rv(0) = 50.35;
114 rv(1) = -3.14;
115 return rv;
116 }
117 else if (mSel == 4) {
118 Vector rv(3);
119 rv(0) = 2.2;
120 rv(1) = -3.6;
121 rv(2) = 1.1;
122 return rv;
123 }
124 else xx = 0.;
125 return xx ;
126}
127
128
129double MinZTestFunc::Value(double const xp[])
130{
131 double retval = 0;
132 if (mSel == 0) {
133 double x = xp[0];
134 retval = -x+(x-2.)*(x-2.);
135 }
136 else if ((mSel == 1) || (mSel == 2)) {
137 double x = xp[0];
138 retval = 0.1*x*x;
139 x = xp[0]-2.;
140 x = x*x;
141 retval -= 3*exp(-x);
142 x = xp[0]+3.;
143 x = 0.5*x*x;
144 if (mSel == 1) retval -= 5*exp(-x);
145 else retval += 5*exp(-x);
146 }
147 else if (mSel == 3) {
148 double x = xp[0]-50.35;
149 double y = xp[1]+3.14;
150 retval = 1.3*x*x+25.*y*y;
151 }
152 else if (mSel == 4) {
153 double x = xp[0]-2.2;
154 double y = xp[1]+3.6;
155 double z = xp[2]-1.1;
156 retval = x*x+2.*y*y+3.*z*z;
157 }
158 else retval = 0.;
159 return retval;
160}
161
162//---------------------------------------------------------------
163//------------------- Classe MinZSimplex --------------------
164//---------------------------------------------------------------
165string __Vec2Str4MinZ_AutoTest(Vector& xx)
166{
167 string rs;
168 char buff[32];
169 for(int i=0; i<xx.Size(); i++) {
170 sprintf(buff," %g " , xx(i));
171 rs += buff;
172 }
173 return rs;
174}
175
176int MinZSimplex::AutoTest(int tsel, int prtlev)
177{
178 int rc = 0;
179 cout << " --- MinZSimplex::AutoTest() --- TSel= " << tsel << " PrtLev=" << prtlev << endl;
180 for(int i=0; i<5; i++) {
181 if ((tsel >= 0) && (tsel != i)) continue;
182 cout << " ======= Test avec ISel= " << i;
183 Vector xx;
184 MinZTestFunc mzf(i);
185 cout << " - Func= " << mzf.ToString() << endl;
186 Vector rv = mzf.OptParms();
187 xx = rv;
188 for(int j=0; j<2; j++) {
189 double vi = 50.*(j-0.5);
190 for(int k=0; k<2; k++) {
191 double vs = (k == 0) ? 1. : 10. ;
192 cout << "--[" << j << "," << k
193 << "] Initialisation avec IniPoint= " << vi << " IniStep= " << vs << endl;
194 MinZSimplex simplex(&mzf);
195 xx = vi;
196 simplex.SetInitialPoint(xx);
197 xx = vs;
198 simplex.SetInitialStep(xx);
199 simplex.SetPrtLevel(prtlev);
200 int rcs = simplex.Minimize(xx);
201 Vector diff = rv-xx;
202 double d2 = diff.Norm2();
203 cout << " Rc(simplex.Minimize() = " << rc << " NIter= "
204 << simplex.NbIter() << " ===> Distance^2= " << d2
205 << "\nConverged to " << __Vec2Str4MinZ_AutoTest(xx)
206 << " Best Value= " << __Vec2Str4MinZ_AutoTest(rv)
207 << " Diff = " << __Vec2Str4MinZ_AutoTest(diff) << endl;
208 if ((rcs > 5) || (d2 > 0.5)) rc ++;
209 }
210 }
211 }
212 cout << " --- MinZSimplex::AutoTest() --- Rc=" << rc << " -- END ----- " << endl;
213 return rc;
214}
215
216MinZSimplex::MinZSimplex(MinZFunction *mzf)
217 : mZF(mzf) , mPoint0(mZF->NVar()) , mStep0(mZF->NVar())
218{
219 SetMaxIter();
220 SetControls();
221 Vector xx(NDim());
222 xx = 0.;
223 SetInitialPoint(xx);
224 xx = 1.0;
225 SetInitialStep(xx);
226 SetStopTolerance();
227 mIter = -1;
228 mStop = -1;
229 SetPrtLevel();
230}
231
232MinZSimplex::~MinZSimplex()
233{
234}
235
236int MinZSimplex::Minimize(Vector& fpoint)
237{
238 // vector< TVector<r_8> > splx;
239 Vector splx[100];
240 Vector Y(NDim()+1);
241 // On calcule le simplex initial
242 // N = NDim, N+1 points (pp) ds l'espace a N dimensions
243 // Point0, Point0 + Step0(i) e_i
244 Vector pp,ppc;
245 pp = mPoint0;
246 //ppc = pp;
247 //splx.push_back(ppc);
248 splx[0] = pp;
249 int i,j,k;
250 for(i=0; i<NDim(); i++) {
251 Vector pps;
252 pps = mPoint0;
253 pps(i) += mStep0(i);
254 //splx.push_back(pps);
255 splx[i+1] = pps;
256 }
257 int mpts = NDim()+1;
258 // calcul des valeurs de la fonction sur les sommets
259 for(i=0; i<mpts; i++)
260 Y(i) = Value(splx[i]);
261
262 int iter = 0;
263 mIter = iter;
264 mStop = 0;
265
266 int nbugrtol2 = 0;
267 bool stop = false, stop0=false;
268 int rc = 0;
269 int ilo, ihi, inhi;
270 int move = 0;
271 char* smov[6] = { "None", "Reflection", "ReflecExpand", "ContractHigh", "ContractLow", "ExpandHigh" };
272 int movcnt[6] = {0,0,0,0,0,0};
273
274 int nrep1=0, nrep2=0;
275 FindMinMax12(Y, ilo, ihi, inhi);
276 double yhilast = Y(ihi);
277 yhilast += fabs(yhilast);
278
279 while (!stop) { //
280 FindMinMax12(Y, ilo, ihi, inhi);
281 double ymean = (fabs(Y(ihi))+fabs(Y(ilo)));
282 if (ymean < mTol0) { stop0 = true; ymean = mTol0; }
283 double rtol1 = 2.*fabs(Y(ihi)-Y(ilo))/ymean;
284 double ym2 = (fabs(yhilast)+fabs(Y(ihi)));
285 if (ym2 < mTol0) ym2 = mTol0;
286 double rtol2 = 2.*(yhilast-Y(ihi))/ym2;
287 yhilast = Y(ihi);
288 if (rtol2 < 0.) {
289 if (move != 40) {
290 cout << " !!!! MinZSimplex::Minimize() BUG RTol2< 0. --> Chs " << endl;
291 nbugrtol2++;
292 }
293 else nrep2 = 0;
294 rtol2 = -rtol2;
295 }
296 if (PrtLevel() > 1)
297 cout << "--MinZSimplex::Minimize() - Iter=" << iter
298 << " Move= " << move << " (" << smov[move/10] << ")" << endl;
299 if (PrtLevel() > 2)
300 cout << "..ILO=" << ilo << " IHI=" << ihi << " INHI=" << inhi
301 << " Y(ILO)=" << Y(ilo) << " Y(IHI)=" << Y(ihi) << "\n"
302 << "...YMean_Abs=" << ymean << " RTOL1=" << rtol1 << " RTOL2=" << rtol2 << endl;
303 if (PrtLevel() > 3) {
304 for(i=0; i<mpts; i++) {
305 cout << "....Simplex[" << i << "]= ";
306 for(j=0; j<NDim(); j++) cout << splx[i](j) << " , ";
307 cout << " Y=Value= " << Y(i) << endl;
308 }
309 }
310 if (rtol1 < mTol1) nrep1++;
311 else nrep1 = 0;
312 if (rtol2 < mTol2) nrep2++;
313 else nrep2 = 0;
314
315 if (stop0) { mStop = 1; rc = 0; stop = true; break; }
316 if (nrep1 > mRep1) { mStop = 2; rc = 0; stop = true; break; }
317 if (nrep2 > mRep2) { mStop = 3; rc = 0; stop = true; break; }
318 if (iter > MaxIter() ) { mStop = 0, rc = iter; break; }
319 iter++;
320 if (iter > 0) movcnt[move/10]++;
321
322 // Next iteration, on modifie le simplex
323 // Calcul du centre de gravite su simplex, hors le point le + haut
324 Vector pbar(NDim());
325 pbar = 0.;
326 for(i=0; i<mpts; i++) {
327 if (i == ihi) continue;
328 pbar += splx[i];
329 }
330 pbar /= (double)NDim();
331 // On calcule le sommet oppose a point IHI (le + haut)
332 Vector pr, prr;
333 double YPR, YPRR;
334 pr = (1.+Alpha())*pbar-Alpha()*splx[ihi];
335 YPR = Value(pr);
336 if (YPR < Y(ilo)) { // Amelioaration par rapport au meilleur point,
337 // on va plus loin d'un facteur gamma
338 prr = Gamma()*pr+(1.-Gamma())*pbar;
339 YPRR = Value(prr);
340 if (YPRR < Y(ilo)) { // On remplace le IHI par YPRR
341 splx[ihi] = prr;
342 Y(ihi) = YPRR;
343 move = 20;
344 }
345 else { // sinon, on remplace par YPR
346 splx[ihi] = pr;
347 Y(ihi) = YPR;
348 move = 10;
349 }
350 }
351 else { // Moins bon que le meilleur point ..
352 if (YPR > Y(inhi)) { // Plus mauvais que le second plus haut (INHI)
353 if (YPR < Y(ihi)) { // Mais meilleur que le plus haut (IHI)
354 splx[ihi] = pr; // On remplace donc le plus haut
355 Y(ihi) = YPR;
356 move = 11;
357 }
358 else { // Plus mauvais que le plus mauvais IHI
359 // on tente avec un point intermediaire
360 prr = Beta()*splx[ihi]+(1.-Beta())*pbar;
361 YPRR = Value(prr);
362 if (YPRR < Y(ihi)) { // Le point intermediaire ameliore les choses
363 splx[ihi] = prr; // On remplace donc le point le + haut
364 Y(ihi) = YPRR;
365 move = 30;
366 }
367 else {
368 // On tente aussi de rester du meme cote, mais aller plus loin
369 prr = Gamma2()*splx[ihi]+(1.-Gamma2())*pbar;
370 YPRR = Value(prr);
371 if (YPRR < Y(ihi)) { // Le point intermediaire ameliore les choses
372 splx[ihi] = prr; // On remplace donc le point le + haut
373 Y(ihi) = YPRR;
374 move = 50;
375 }
376 else {
377 // Rien n'y fait, on contracte autour du meilleur point
378 for(i=0; i<mpts; i++) {
379 if (i == ilo) continue;
380 splx[i] = Beta2()*splx[i]+(1.-Beta())*splx[ilo];
381 Y(i) = Value(splx[i]);
382 move = 40;
383 }
384 }
385 }
386 }
387 }
388 else { // Meilleur que le IHI et le INHI
389 splx[ihi] = pr; // On remplace le plus haut
390 Y(ihi) = YPR;
391 move = 12;
392 }
393 }
394 } // Fin de la boucle while principale
395
396 fpoint = splx[ilo];
397 mIter = iter;
398
399 if (PrtLevel() > 0) {
400 string sr;
401 StopReason(sr);
402 cout << "-----MinZSimplex::Minimize()/Ended - NIter=" << iter
403 << " Moves[0..5]= " << movcnt[0] << "," << movcnt[1] << ","
404 << movcnt[2] << "," << movcnt[3] << ","
405 << movcnt[4] << "," << movcnt[5]
406 << "\n..MinZSimplex Stop=" << StopReason() << " -> " << sr << endl;
407
408 if (nbugrtol2 > 0) cout << "MinZSimplex::Minimize()/Warning - nbugrtol2= " << nbugrtol2 << endl;
409 }
410 return rc;
411}
412
413int MinZSimplex::StopReason(string& s)
414{
415 char* sr[5] = { "NoConverg, MaxIterReached", "OK, fm<Tol0", "OK, Df/f<Tol1",
416 "OK, [Df/f max]Iter<Tol2" "Error - Wrong StopReason" };
417 int stop = mStop;
418 if ((stop < 0) || (stop > 3)) stop = 4;
419 s = sr[stop];
420 return mStop;
421}
422
423int MinZSimplex::FindMinMax12(Vector& fval, int& ilo, int& ihi, int& inhi)
424{
425 ilo = 0;
426 if (fval(0) > fval(1)) { ihi = 0; inhi = 1; }
427 else { ihi = 1; inhi = 0; }
428
429 for(int k=0; k<fval.Size(); k++) {
430 if (fval(k) < fval(ilo)) ilo = k;
431 if (fval(k) > fval(ihi)) {
432 inhi = ihi;
433 ihi = k;
434 }
435 else if (fval(k) > fval(inhi)) {
436 if (k != ihi) inhi = k; // ce test n'est peut-etre pas necessaire ???
437 }
438 }
439 return ilo;
440}
Note: See TracBrowser for help on using the repository browser.