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

Last change on this file since 3002 was 2808, checked in by ansari, 20 years ago

MAJ documentation - Reza 14/6/2005
File size: 14.5 KB
RevLine 
[2650]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
[2808]12/*!
13 \class SOPHYA::MinZFunction
14 \ingroup NTools
15 Interface definition for a function object f(x[]) for which MinZSimplex can
16 search the minimum.
17 The pure virtual method Value() should be implemented by the derived classes.
18*/
[2650]19
20MinZFunction::MinZFunction(unsigned int nvar)
21 : mNVar(nvar)
22{
23}
24
25MinZFunction::~MinZFunction()
26{
27}
28
29//---------------------------------------------------------------
30//------------------- Classe MinZFuncXi2 --------------------
31//---------------------------------------------------------------
[2808]32/*!
33 \class SOPHYA::MinZXi2
34 \ingroup NTools
35 Implements the MinZFunction interface using a xi2 calculator
36 \sa GeneralXi2 GeneralFitData
37*/
[2650]38MinZFuncXi2::MinZFuncXi2(GeneralXi2* gxi2, GeneralFitData* gd)
39 : mGXi2(gxi2) , mGData(gd), MinZFunction(gxi2->NPar())
40{
41}
42
43MinZFuncXi2::~MinZFuncXi2()
44{
45}
46
47double MinZFuncXi2::Value(double const xp[])
48{
49 int ndataused;
50 return mGXi2->Value(*mGData, const_cast<double *>(xp), ndataused);
51}
52
53//---------------------------------------------------------------
54//------------------- Classe MinZTestFunc -------------------
55//---------------------------------------------------------------
56class MinZTestFunc : public MinZFunction {
57public:
58 MinZTestFunc(int sel);
59 virtual double Value(double const xp[]);
60 string ToString();
61 Vector OptParms();
62protected:
63 static int ISelToNvar(int isel);
64 int mSel;
65};
66
67int MinZTestFunc::ISelToNvar(int isel)
68{
69 if (isel == 0) return 1;
70 if (isel == 1) return 1;
71 else if (isel == 2) return 1;
72 else if (isel == 3) return 2;
73 else if (isel == 4) return 3;
74 else return 1;
75}
76
77MinZTestFunc::MinZTestFunc(int sel)
78 : MinZFunction(ISelToNvar(sel))
79{
80 if ((sel < 0) || (sel > 4)) sel = 0;
81 mSel = sel;
82}
83
84string MinZTestFunc::ToString()
85{
86 string rs;
87 if (mSel == 0) {
88 rs = "-x+(x-2)^2";
89 }
90 else if (mSel == 1) {
91 rs = "0.1*x^2-3exp(-(x-2)^2)-5*exp(-0.5*(x+3)^2)";
92 }
93 else if (mSel == 2) {
94 rs = "0.1*x^2-3exp(-(x-2)^2)+5*exp(-0.5*(x+3)^2)";
95 }
96 else if (mSel == 3) {
97 rs = "1.3*(x-50.35)^2+25*(y+3.14)^2";
98 }
99 else if (mSel == 4) {
100 rs = "(x-2.2)^2+2.*(y+3.6)^2+3.*(z-1.1)^2";
101 }
102 else rs = "????";
103 return rs;
104}
105
106Vector MinZTestFunc::OptParms()
107{
108 Vector xx;
109 if (mSel == 0) {
110 Vector rv(1);
111 rv = 2.5;
112 return rv;
113 }
114 else if (mSel == 1) {
115 Vector rv(1);
116 rv = -2.883;
117 return rv;
118 }
119 else if (mSel == 2) {
120 Vector rv(1);
121 rv = 1.812;
122 return rv;
123 }
124 else if (mSel == 3) {
125 Vector rv(2);
126 rv(0) = 50.35;
127 rv(1) = -3.14;
128 return rv;
129 }
130 else if (mSel == 4) {
131 Vector rv(3);
132 rv(0) = 2.2;
133 rv(1) = -3.6;
134 rv(2) = 1.1;
135 return rv;
136 }
137 else xx = 0.;
138 return xx ;
139}
140
141
142double MinZTestFunc::Value(double const xp[])
143{
144 double retval = 0;
145 if (mSel == 0) {
146 double x = xp[0];
147 retval = -x+(x-2.)*(x-2.);
148 }
149 else if ((mSel == 1) || (mSel == 2)) {
150 double x = xp[0];
151 retval = 0.1*x*x;
152 x = xp[0]-2.;
153 x = x*x;
154 retval -= 3*exp(-x);
155 x = xp[0]+3.;
156 x = 0.5*x*x;
157 if (mSel == 1) retval -= 5*exp(-x);
158 else retval += 5*exp(-x);
159 }
160 else if (mSel == 3) {
161 double x = xp[0]-50.35;
162 double y = xp[1]+3.14;
163 retval = 1.3*x*x+25.*y*y;
164 }
165 else if (mSel == 4) {
166 double x = xp[0]-2.2;
167 double y = xp[1]+3.6;
168 double z = xp[2]-1.1;
169 retval = x*x+2.*y*y+3.*z*z;
170 }
171 else retval = 0.;
172 return retval;
173}
174
175//---------------------------------------------------------------
176//------------------- Classe MinZSimplex --------------------
177//---------------------------------------------------------------
178string __Vec2Str4MinZ_AutoTest(Vector& xx)
179{
180 string rs;
181 char buff[32];
182 for(int i=0; i<xx.Size(); i++) {
183 sprintf(buff," %g " , xx(i));
184 rs += buff;
185 }
186 return rs;
187}
188
[2808]189/*!
190 \class SOPHYA::MinZSimplex
191 \ingroup NTools
192 This class implements non linear minimization (optimization)
193 in a multidimensional space following the \b Simplex method.
194 A \b Simplex is a geometrical figure made of N+1 points in a
195 N-dimensional space. (triangle in a plane, tetrahedron in 3-d space).
196 The minimization method implemented in this class is based on the
197 algorithm described in "Numerical Recipes, Chapter X".
198
199 The algorithm has been slightly enhanced :
200 - More complex convergence / stop test
201 - A new transformation of the simplex has been included (ExpandHigh)
202
203 For each step, on of the following geometrical transform is performed
204 on the Simplex figure:
205 - Reflection : reflection away from the high point (expansion by factor Alpha)
206 - ReflecExpand : reflection way from the high point and expansion by factor Beta2
207 - ContractHigh : Contraction along the high point (factor Beta)
208 - ContractLow : Contraction toward the low point (factor Beta2)
209 - ExpandHigh : Expansion along the high point
210
211 \sa GeneralFit
212
213 The following sample code shows a usage example:
214 \code
215 include "simplex.h"
216 ...
217 // Define our function to be minimized:
218 class MySFunc : public MinZFunction {
219 public:
220 MySFunc() : MinZFunction(2) {}
221 virtual double Value(double const xp[])
222 { return (xp[0]*xp[0]+2*xp[1]*xp[1]); }
223 };
224
225 ...
226
227 MySFunc mysf;
228 MinZSimplex simplex(&mysf);
229 // Guess the center and step for constructing the initial simplex
230 Vector x0(2); x0 = 1.;
231 Vector step(2); step = 2.;
232 simplex.SetInitialPoint(x0);
233 simplex.SetInitialStep(step);
234 Vector oparm(2);
235 int rc = simplex.Minimize(oparm);
236 if (rc != 0) {
237 string srt;
238 int sr = simplex.StopReason(srt);
239 cout << " Convergence Pb, StopReason= " << sr << " : " << srt << endl;
240 }
241 else {
242 cout << " Converged: NStep= " << simplex.NbIter()
243 << " OParm= " << oparm << endl;
244 }
245 \endcode
246*/
247
248/*!
249 \brief Auto test function
250 \param tsel : select autotest (0,1,2,3,4) , tsel<0 -> all
251 \param prtlev : printlevel
252*/
[2650]253int MinZSimplex::AutoTest(int tsel, int prtlev)
254{
255 int rc = 0;
256 cout << " --- MinZSimplex::AutoTest() --- TSel= " << tsel << " PrtLev=" << prtlev << endl;
257 for(int i=0; i<5; i++) {
258 if ((tsel >= 0) && (tsel != i)) continue;
259 cout << " ======= Test avec ISel= " << i;
260 Vector xx;
261 MinZTestFunc mzf(i);
262 cout << " - Func= " << mzf.ToString() << endl;
263 Vector rv = mzf.OptParms();
264 xx = rv;
265 for(int j=0; j<2; j++) {
266 double vi = 50.*(j-0.5);
267 for(int k=0; k<2; k++) {
268 double vs = (k == 0) ? 1. : 10. ;
269 cout << "--[" << j << "," << k
270 << "] Initialisation avec IniPoint= " << vi << " IniStep= " << vs << endl;
271 MinZSimplex simplex(&mzf);
272 xx = vi;
273 simplex.SetInitialPoint(xx);
274 xx = vs;
275 simplex.SetInitialStep(xx);
276 simplex.SetPrtLevel(prtlev);
277 int rcs = simplex.Minimize(xx);
278 Vector diff = rv-xx;
279 double d2 = diff.Norm2();
280 cout << " Rc(simplex.Minimize() = " << rc << " NIter= "
281 << simplex.NbIter() << " ===> Distance^2= " << d2
282 << "\nConverged to " << __Vec2Str4MinZ_AutoTest(xx)
283 << " Best Value= " << __Vec2Str4MinZ_AutoTest(rv)
284 << " Diff = " << __Vec2Str4MinZ_AutoTest(diff) << endl;
285 if ((rcs > 5) || (d2 > 0.5)) rc ++;
286 }
287 }
288 }
289 cout << " --- MinZSimplex::AutoTest() --- Rc=" << rc << " -- END ----- " << endl;
290 return rc;
291}
292
[2808]293//! Constructor from pointer to MinZFunction object
[2650]294MinZSimplex::MinZSimplex(MinZFunction *mzf)
295 : mZF(mzf) , mPoint0(mZF->NVar()) , mStep0(mZF->NVar())
296{
297 SetMaxIter();
298 SetControls();
299 Vector xx(NDim());
300 xx = 0.;
301 SetInitialPoint(xx);
302 xx = 1.0;
303 SetInitialStep(xx);
304 SetStopTolerance();
305 mIter = -1;
306 mStop = -1;
307 SetPrtLevel();
308}
309
310MinZSimplex::~MinZSimplex()
311{
312}
313
[2808]314//! Perform the minimization
315/*!
316 Return 0 if success
317 \param fpoint : vector containing the optimal point
318
319 Convergence test :
320 \verbatim
321 On minimise f(x) f=mZF->Value() ,
322 f_max = max(f) sur simplex , f_min = min(f) sur simplex
323 fm = (abs(f_max)+abs(f_min))
324 [Delta f] = abs(f_max-f_min)
325 [Delta f/f]simplex = 2.*Delta f / fm
326 fm2 = (abs(f_max)+abs(f_max(iter-1)))
327 [Delta f_max/f_max]iter = [f_max(iter-1)-f_max]/fm2
328 Test d'arret :
329 fm < mTol0 OU
330 [Delta f/f]simplex < mTol1 mRep1 fois de suite OU
331 [Delta f_max/f_max]iter < mTol2 mRep2 fois de suite
332*/
[2650]333int MinZSimplex::Minimize(Vector& fpoint)
334{
335 // vector< TVector<r_8> > splx;
336 Vector splx[100];
337 Vector Y(NDim()+1);
338 // On calcule le simplex initial
339 // N = NDim, N+1 points (pp) ds l'espace a N dimensions
340 // Point0, Point0 + Step0(i) e_i
341 Vector pp,ppc;
342 pp = mPoint0;
343 //ppc = pp;
344 //splx.push_back(ppc);
345 splx[0] = pp;
346 int i,j,k;
347 for(i=0; i<NDim(); i++) {
348 Vector pps;
349 pps = mPoint0;
350 pps(i) += mStep0(i);
351 //splx.push_back(pps);
352 splx[i+1] = pps;
353 }
354 int mpts = NDim()+1;
355 // calcul des valeurs de la fonction sur les sommets
356 for(i=0; i<mpts; i++)
357 Y(i) = Value(splx[i]);
358
359 int iter = 0;
360 mIter = iter;
361 mStop = 0;
362
363 int nbugrtol2 = 0;
364 bool stop = false, stop0=false;
365 int rc = 0;
366 int ilo, ihi, inhi;
367 int move = 0;
368 char* smov[6] = { "None", "Reflection", "ReflecExpand", "ContractHigh", "ContractLow", "ExpandHigh" };
369 int movcnt[6] = {0,0,0,0,0,0};
370
371 int nrep1=0, nrep2=0;
372 FindMinMax12(Y, ilo, ihi, inhi);
373 double yhilast = Y(ihi);
374 yhilast += fabs(yhilast);
375
376 while (!stop) { //
377 FindMinMax12(Y, ilo, ihi, inhi);
378 double ymean = (fabs(Y(ihi))+fabs(Y(ilo)));
379 if (ymean < mTol0) { stop0 = true; ymean = mTol0; }
380 double rtol1 = 2.*fabs(Y(ihi)-Y(ilo))/ymean;
381 double ym2 = (fabs(yhilast)+fabs(Y(ihi)));
382 if (ym2 < mTol0) ym2 = mTol0;
383 double rtol2 = 2.*(yhilast-Y(ihi))/ym2;
384 yhilast = Y(ihi);
385 if (rtol2 < 0.) {
386 if (move != 40) {
387 cout << " !!!! MinZSimplex::Minimize() BUG RTol2< 0. --> Chs " << endl;
388 nbugrtol2++;
389 }
390 else nrep2 = 0;
391 rtol2 = -rtol2;
392 }
393 if (PrtLevel() > 1)
394 cout << "--MinZSimplex::Minimize() - Iter=" << iter
395 << " Move= " << move << " (" << smov[move/10] << ")" << endl;
396 if (PrtLevel() > 2)
397 cout << "..ILO=" << ilo << " IHI=" << ihi << " INHI=" << inhi
398 << " Y(ILO)=" << Y(ilo) << " Y(IHI)=" << Y(ihi) << "\n"
399 << "...YMean_Abs=" << ymean << " RTOL1=" << rtol1 << " RTOL2=" << rtol2 << endl;
400 if (PrtLevel() > 3) {
401 for(i=0; i<mpts; i++) {
402 cout << "....Simplex[" << i << "]= ";
403 for(j=0; j<NDim(); j++) cout << splx[i](j) << " , ";
404 cout << " Y=Value= " << Y(i) << endl;
405 }
406 }
407 if (rtol1 < mTol1) nrep1++;
408 else nrep1 = 0;
409 if (rtol2 < mTol2) nrep2++;
410 else nrep2 = 0;
411
412 if (stop0) { mStop = 1; rc = 0; stop = true; break; }
413 if (nrep1 > mRep1) { mStop = 2; rc = 0; stop = true; break; }
414 if (nrep2 > mRep2) { mStop = 3; rc = 0; stop = true; break; }
415 if (iter > MaxIter() ) { mStop = 0, rc = iter; break; }
416 iter++;
417 if (iter > 0) movcnt[move/10]++;
418
419 // Next iteration, on modifie le simplex
420 // Calcul du centre de gravite su simplex, hors le point le + haut
421 Vector pbar(NDim());
422 pbar = 0.;
423 for(i=0; i<mpts; i++) {
424 if (i == ihi) continue;
425 pbar += splx[i];
426 }
427 pbar /= (double)NDim();
428 // On calcule le sommet oppose a point IHI (le + haut)
429 Vector pr, prr;
430 double YPR, YPRR;
431 pr = (1.+Alpha())*pbar-Alpha()*splx[ihi];
432 YPR = Value(pr);
433 if (YPR < Y(ilo)) { // Amelioaration par rapport au meilleur point,
434 // on va plus loin d'un facteur gamma
435 prr = Gamma()*pr+(1.-Gamma())*pbar;
436 YPRR = Value(prr);
437 if (YPRR < Y(ilo)) { // On remplace le IHI par YPRR
438 splx[ihi] = prr;
439 Y(ihi) = YPRR;
440 move = 20;
441 }
442 else { // sinon, on remplace par YPR
443 splx[ihi] = pr;
444 Y(ihi) = YPR;
445 move = 10;
446 }
447 }
448 else { // Moins bon que le meilleur point ..
449 if (YPR > Y(inhi)) { // Plus mauvais que le second plus haut (INHI)
450 if (YPR < Y(ihi)) { // Mais meilleur que le plus haut (IHI)
451 splx[ihi] = pr; // On remplace donc le plus haut
452 Y(ihi) = YPR;
453 move = 11;
454 }
455 else { // Plus mauvais que le plus mauvais IHI
456 // on tente avec un point intermediaire
457 prr = Beta()*splx[ihi]+(1.-Beta())*pbar;
458 YPRR = Value(prr);
459 if (YPRR < Y(ihi)) { // Le point intermediaire ameliore les choses
460 splx[ihi] = prr; // On remplace donc le point le + haut
461 Y(ihi) = YPRR;
462 move = 30;
463 }
464 else {
465 // On tente aussi de rester du meme cote, mais aller plus loin
466 prr = Gamma2()*splx[ihi]+(1.-Gamma2())*pbar;
467 YPRR = Value(prr);
468 if (YPRR < Y(ihi)) { // Le point intermediaire ameliore les choses
469 splx[ihi] = prr; // On remplace donc le point le + haut
470 Y(ihi) = YPRR;
471 move = 50;
472 }
473 else {
474 // Rien n'y fait, on contracte autour du meilleur point
475 for(i=0; i<mpts; i++) {
476 if (i == ilo) continue;
477 splx[i] = Beta2()*splx[i]+(1.-Beta())*splx[ilo];
478 Y(i) = Value(splx[i]);
479 move = 40;
480 }
481 }
482 }
483 }
484 }
485 else { // Meilleur que le IHI et le INHI
486 splx[ihi] = pr; // On remplace le plus haut
487 Y(ihi) = YPR;
488 move = 12;
489 }
490 }
491 } // Fin de la boucle while principale
492
493 fpoint = splx[ilo];
494 mIter = iter;
495
496 if (PrtLevel() > 0) {
497 string sr;
498 StopReason(sr);
499 cout << "-----MinZSimplex::Minimize()/Ended - NIter=" << iter
500 << " Moves[0..5]= " << movcnt[0] << "," << movcnt[1] << ","
501 << movcnt[2] << "," << movcnt[3] << ","
502 << movcnt[4] << "," << movcnt[5]
503 << "\n..MinZSimplex Stop=" << StopReason() << " -> " << sr << endl;
504
505 if (nbugrtol2 > 0) cout << "MinZSimplex::Minimize()/Warning - nbugrtol2= " << nbugrtol2 << endl;
506 }
507 return rc;
508}
509
[2808]510//! Return the stop reason and fills the corresponding string description
[2650]511int MinZSimplex::StopReason(string& s)
512{
513 char* sr[5] = { "NoConverg, MaxIterReached", "OK, fm<Tol0", "OK, Df/f<Tol1",
514 "OK, [Df/f max]Iter<Tol2" "Error - Wrong StopReason" };
515 int stop = mStop;
516 if ((stop < 0) || (stop > 3)) stop = 4;
517 s = sr[stop];
518 return mStop;
519}
520
521int MinZSimplex::FindMinMax12(Vector& fval, int& ilo, int& ihi, int& inhi)
522{
523 ilo = 0;
524 if (fval(0) > fval(1)) { ihi = 0; inhi = 1; }
525 else { ihi = 1; inhi = 0; }
526
527 for(int k=0; k<fval.Size(); k++) {
528 if (fval(k) < fval(ilo)) ilo = k;
529 if (fval(k) > fval(ihi)) {
530 inhi = ihi;
531 ihi = k;
532 }
533 else if (fval(k) > fval(inhi)) {
534 if (k != ihi) inhi = k; // ce test n'est peut-etre pas necessaire ???
535 }
536 }
537 return ilo;
538}
Note: See TracBrowser for help on using the repository browser.