Logo ROOT  
Reference Guide
MnLineSearch.cxx
Go to the documentation of this file.
1// @(#)root/minuit2:$Id$
2// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
3
4/**********************************************************************
5 * *
6 * Copyright (c) 2005 LCG ROOT Math team, CERN/PH-SFT *
7 * *
8 **********************************************************************/
9
11#include "Minuit2/MnFcn.h"
14#include "Minuit2/MnParabola.h"
17#include "Minuit2/LaSum.h"
18
19#include <iostream>
20#include "Minuit2/MnPrint.h"
21#include <algorithm>
22
23#ifdef USE_OTHER_LS
24
25#include "Math/SMatrix.h"
26#include "Math/SVector.h"
27
28#include "Math/IFunction.h"
29#include "Math/Minimizer1D.h"
30
31#endif
32
33//#define DEBUG
34
35namespace ROOT {
36
37 namespace Minuit2 {
38
39
40/** Perform a line search from position defined by the vector st
41 along the direction step, where the length of vector step
42 gives the expected position of Minimum.
43 fcn is Value of function at the starting position ,
44 gdel (if non-zero) is df/dx along step at st.
45 Return a parabola point containing Minimum x position and y (function Value)
46 - add a falg to control the debug
47*/
48
49MnParabolaPoint MnLineSearch::operator()(const MnFcn& fcn, const MinimumParameters& st, const MnAlgebraicVector& step, double gdel, const MnMachinePrecision& prec, bool debug) const {
50
51 //*-*-*-*-*-*-*-*-*-*Perform a line search from position st along step *-*-*-*-*-*-*-*
52 //*-* =========================================
53 //*-* SLAMBG and ALPHA control the maximum individual steps allowed.
54 //*-* The first step is always =1. The max length of second step is SLAMBG.
55 //*-* The max size of subsequent steps is the maximum previous successful
56 //*-* step multiplied by ALPHA + the size of most recent successful step,
57 //*-* but cannot be smaller than SLAMBG.
58 //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
59
60#ifdef DEBUG
61 debug = true;
62#endif
63
64 if (debug) {
65 std::cout<<"gdel= "<<gdel<<std::endl;
66 std::cout<<"step= "<<step<<std::endl;
67 }
68
69 double overal = 1000.;
70 double undral = -100.;
71 double toler = 0.05;
72 double slamin = 0.;
73 double slambg = 5.;
74 double alpha = 2.;
75 int maxiter = 12;
76 // start as in Fortran from 1 and count all the time we evaluate the function
77 int niter = 1;
78
79 for(unsigned int i = 0; i < step.size(); i++) {
80 if(step(i) == 0 ) continue;
81 double ratio = fabs(st.Vec()(i)/step(i));
82 if( slamin == 0) slamin = ratio;
83 if(ratio < slamin) slamin = ratio;
84 }
85 if(fabs(slamin) < prec.Eps()) slamin = prec.Eps();
86 slamin *= prec.Eps2();
87
88 double f0 = st.Fval();
89 double f1 = fcn(st.Vec()+step);
90 niter++;
91 double fvmin = st.Fval();
92 double xvmin = 0.;
93
94 if(f1 < f0) {
95 fvmin = f1;
96 xvmin = 1.;
97 }
98 double toler8 = toler;
99 double slamax = slambg;
100 double flast = f1;
101 double slam = 1.;
102
103 bool iterate = false;
104 MnParabolaPoint p0(0., f0);
105 MnParabolaPoint p1(slam, flast);
106 double f2 = 0.;
107 // quadratic interpolation using the two points p0,p1 and the slope at p0
108 do {
109 // cut toler8 as function goes up
110 iterate = false;
111 //MnParabola pb = MnParabolaFactory()(p0, gdel, p1);
112
113 if (debug) {
114 // std::cout<<"pb.Min() = "<<pb.Min()<<std::endl;
115 std::cout<<"flast, f0= "<<flast<<", "<<f0<<std::endl;
116 std::cout<<"flast-f0= "<<flast-f0<<std::endl;
117 std::cout<<"slam= "<<slam<<std::endl;
118 }
119 // double df = flast-f0;
120 // if(fabs(df) < prec.Eps2()) {
121 // if(flast-f0 < 0.) df = -prec.Eps2();
122 // else df = prec.Eps2();
123 // }
124 // std::cout<<"df= "<<df<<std::endl;
125 // double denom = 2.*(df-gdel*slam)/(slam*slam);
126 double denom = 2.*(flast-f0-gdel*slam)/(slam*slam);
127 if (debug) std::cout<<"denom= "<<denom<<std::endl;
128 if(denom != 0) {
129 slam = - gdel/denom;
130 } else {
131 denom = -0.1*gdel;
132 slam = 1.;
133 }
134 if (debug) std::cout<<"new slam= "<<slam<<std::endl;
135
136
137#ifdef TRY_OPPOSIT_DIR
138 // if is less than zero indicates maximum position. Use then slamax or x = x0 - 2slam + 1
139 bool slamIsNeg = false;
140 double slamNeg = 0;
141#endif
142
143 if(slam < 0.) {
144 if (debug) std::cout<<"slam is negative- set to " << slamax << std::endl;
145#ifdef TRY_OPPOSITE_DIR
146 slamNeg = 2*slam -1;
147 slamIsNeg = true;
148 if (debug) std::cout<<"slam is negative- compare values between "<< slamNeg << " and " << slamax << std::endl;
149#endif
150 slam = slamax;
151 }
152 // std::cout<<"slam= "<<slam<<std::endl;
153 if(slam > slamax) {
154 slam = slamax;
155 if (debug) std::cout << "slam larger than mac value - set to " << slamax << std::endl;
156 }
157
158 if(slam < toler8) {
159 if (debug) std::cout << "slam too small - set to " << toler8 << std::endl;
160 slam = toler8;
161 }
162 // std::cout<<"slam= "<<slam<<std::endl;
163 if(slam < slamin) {
164 if (debug) std::cout << "slam smaller than " << slamin << " return " << std::endl;
165 // std::cout<<"f1, f2= "<<p0.Y()<<", "<<p1.Y()<<std::endl;
166 // std::cout<<"x1, x2= "<<p0.X()<<", "<<p1.X()<<std::endl;
167 // std::cout<<"x, f= "<<xvmin<<", "<<fvmin<<std::endl;
168 return MnParabolaPoint(xvmin, fvmin);
169 }
170 if(fabs(slam - 1.) < toler8 && p1.Y() < p0.Y()) {
171 // std::cout<<"f1, f2= "<<p0.Y()<<", "<<p1.Y()<<std::endl;
172 // std::cout<<"x1, x2= "<<p0.X()<<", "<<p1.X()<<std::endl;
173 // std::cout<<"x, f= "<<xvmin<<", "<<fvmin<<std::endl;
174 return MnParabolaPoint(xvmin, fvmin);
175 }
176 if(fabs(slam - 1.) < toler8) slam = 1. + toler8;
177
178 // if(fabs(gdel) < prec.Eps2() && fabs(denom) < prec.Eps2())
179 // slam = 1000.;
180 // MnAlgebraicVector tmp = step;
181 // tmp *= slam;
182 // f2 = fcn(st.Vec()+tmp);
183 f2 = fcn(st.Vec() + slam*step);
184
185 niter++; // do as in Minuit (count all func evalu)
186
187#ifdef TRY_OPPOSITE_DIR
188 if (slamIsNeg) {
189 // try alternative in the opposite direction
190 double f2alt = fcn(st.Vec() + slamNeg*step);
191 if (f2alt < f2) {
192 slam = slamNeg;
193 f2 = f2alt;
194 undral += slam;
195 }
196 }
197#endif
198 if(f2 < fvmin) {
199 fvmin = f2;
200 xvmin = slam;
201 }
202 // LM : correct a bug using precision
203 if (fabs( p0.Y() - fvmin) < fabs(fvmin)*prec.Eps() ) {
204 // if(p0.Y()-prec.Eps() < fvmin && fvmin < p0.Y()+prec.Eps()) {
205 iterate = true;
206 flast = f2;
207 toler8 = toler*slam;
208 overal = slam - toler8;
209 slamax = overal;
210 p1 = MnParabolaPoint(slam, flast);
211 //niter++;
212 }
213 } while(iterate && niter < maxiter);
214 if(niter >= maxiter) {
215 // exhausted max number of iterations
216 return MnParabolaPoint(xvmin, fvmin);
217 }
218
219 if (debug){
220 std::cout<<"after initial 2-point iter: "<<std::endl;
221 std::cout<<"x0, x1, x2= "<<p0.X()<<", "<<p1.X()<<", "<<slam<<std::endl;
222 std::cout<<"f0, f1, f2= "<<p0.Y()<<", "<<p1.Y()<<", "<<f2<<std::endl;
223 }
224
225 MnParabolaPoint p2(slam, f2);
226
227 // do now the quadratic interpolation with 3 points
228 do {
229 slamax = std::max(slamax, alpha*fabs(xvmin));
230 MnParabola pb = MnParabolaFactory()(p0, p1, p2);
231 if (debug) {
232 std::cout << "\nLS Iteration " << niter << std::endl;
233 std::cout<<"x0, x1, x2= "<<p0.X()<<", "<<p1.X()<<", "<<p2.X() <<std::endl;
234 std::cout<<"f0, f1, f2= "<<p0.Y()<<", "<<p1.Y()<<", "<<p2.Y() <<std::endl;
235 std::cout << "slamax = " << slamax << std::endl;
236 std::cout<<"p2-p0,p1: "<<p2.Y() - p0.Y()<<", "<<p2.Y() - p1.Y()<<std::endl;
237 std::cout<<"a, b, c= "<<pb.A()<<" "<<pb.B()<<" "<<pb.C()<<std::endl;
238 }
239 if(pb.A() < prec.Eps2()) {
240 double slopem = 2.*pb.A()*xvmin + pb.B();
241 if(slopem < 0.) slam = xvmin + slamax;
242 else slam = xvmin - slamax;
243 if (debug) std::cout << "xvmin = " << xvmin << " slopem " << slopem << " slam " << slam << std::endl;
244 } else {
245 slam = pb.Min();
246 // std::cout<<"pb.Min() slam= "<<slam<<std::endl;
247 if(slam > xvmin + slamax) slam = xvmin + slamax;
248 if(slam < xvmin - slamax) slam = xvmin - slamax;
249 }
250 if(slam > 0.) {
251 if(slam > overal) slam = overal;
252 } else {
253 if(slam < undral) slam = undral;
254 }
255
256 if (debug) {
257 std::cout<<" slam= "<<slam<< " undral " << undral << " overal " << overal << std::endl;
258 }
259
260 double f3 = 0.;
261 do {
262
263 if (debug) std::cout << " iterate on f3- slam " << niter << " slam = " << slam << " xvmin " << xvmin << std::endl;
264
265 iterate = false;
266 double toler9 = std::max(toler8, fabs(toler8*slam));
267 // min. of parabola at one point
268 if(fabs(p0.X() - slam) < toler9 ||
269 fabs(p1.X() - slam) < toler9 ||
270 fabs(p2.X() - slam) < toler9) {
271 // std::cout<<"f1, f2, f3= "<<p0.Y()<<", "<<p1.Y()<<", "<<p2.Y()<<std::endl;
272 // std::cout<<"x1, x2, x3= "<<p0.X()<<", "<<p1.X()<<", "<<p2.X()<<std::endl;
273 // std::cout<<"x, f= "<<xvmin<<", "<<fvmin<<std::endl;
274 return MnParabolaPoint(xvmin, fvmin);
275 }
276
277 // take the step
278 // MnAlgebraicVector tmp = step;
279 // tmp *= slam;
280 f3 = fcn(st.Vec() + slam*step);
281 if (debug) {
282 std::cout<<"f3= "<<f3<<std::endl;
283 std::cout<<"f3-p(2-0).Y()= "<<f3-p2.Y()<<" "<<f3-p1.Y()<<" "<<f3-p0.Y()<<std::endl;
284 }
285 // if latest point worse than all three previous, cut step
286 if(f3 > p0.Y() && f3 > p1.Y() && f3 > p2.Y()) {
287 if (debug) {
288 std::cout<<"f3 worse than all three previous"<<std::endl;
289 }
290 if(slam > xvmin) overal = std::min(overal, slam-toler8);
291 if(slam < xvmin) undral = std::max(undral, slam+toler8);
292 slam = 0.5*(slam + xvmin);
293 if (debug) {
294 std::cout<<"new slam= "<<slam<<std::endl;
295 }
296 iterate = true;
297 niter++;
298 }
299 } while(iterate && niter < maxiter);
300 if(niter >= maxiter) {
301 // exhausted max number of iterations
302 return MnParabolaPoint(xvmin, fvmin);
303 }
304
305 // find worst previous point out of three and replace
306 MnParabolaPoint p3(slam, f3);
307 if(p0.Y() > p1.Y() && p0.Y() > p2.Y()) p0 = p3;
308 else if(p1.Y() > p0.Y() && p1.Y() > p2.Y()) p1 = p3;
309 else p2 = p3;
310 if (debug) std::cout << " f3 " << f3 << " fvmin " << fvmin << " xvmin " << xvmin << std::endl;
311 if(f3 < fvmin) {
312 fvmin = f3;
313 xvmin = slam;
314 } else {
315 if(slam > xvmin) overal = std::min(overal, slam-toler8);
316 if(slam < xvmin) undral = std::max(undral, slam+toler8);
317 }
318
319 niter++;
320 } while(niter < maxiter);
321
322 if (debug) {
323 std::cout<<"f1, f2= "<<p0.Y()<<", "<<p1.Y()<<std::endl;
324 std::cout<<"x1, x2= "<<p0.X()<<", "<<p1.X()<<std::endl;
325 std::cout<<"x, f= "<<xvmin<<", "<<fvmin<<std::endl;
326 }
327 return MnParabolaPoint(xvmin, fvmin);
328
329}
330
331#ifdef USE_OTHER_LS
332
333/** Perform a line search using a cubic interpolation using x0, x1 , df/dx(x0) and d2/dx(x0) (second derivative)
334 This is used at the beginning when the second derivative is known to be negative
335*/
336
337MnParabolaPoint MnLineSearch::CubicSearch(const MnFcn& fcn, const MinimumParameters& st, const MnAlgebraicVector& step, double gdel, double g2del, const MnMachinePrecision& prec, bool debug) const {
338
339
340
341 if (debug) {
342 std::cout<<"gdel= "<<gdel<<std::endl;
343 std::cout<<"g2del= "<<g2del<<std::endl;
344 std::cout<<"step= "<<step<<std::endl;
345 }
346
347 // change ot large values
348 double overal = 100.;
349 double undral = -100.;
350 double toler = 0.05;
351 double slamin = 0.;
352 double slambg = 5.;
353 double alpha = 2.;
354
355 for(unsigned int i = 0; i < step.size(); i++) {
356 if(step(i) == 0 ) continue;
357 double ratio = fabs(st.Vec()(i)/step(i));
358 if( slamin == 0) slamin = ratio;
359 if(ratio < slamin) slamin = ratio;
360 }
361 if(fabs(slamin) < prec.Eps()) slamin = prec.Eps();
362 slamin *= prec.Eps2();
363
364 double f0 = st.Fval();
365 double f1 = fcn(st.Vec()+step);
366 double fvmin = st.Fval();
367 double xvmin = 0.;
368 if (debug) std::cout << "f0 " << f0 << " f1 " << f1 << std::endl;
369
370 if(f1 < f0) {
371 fvmin = f1;
372 xvmin = 1.;
373 }
374 double toler8 = toler;
375 double slamax = slambg;
376 double flast = f1;
377 double slam = 1.;
378
379// MnParabolaPoint p0(0., f0);
380// MnParabolaPoint p1(slam, flast);
381
383 ROOT::Math::SVector<double, 3> cubicCoeff; // cubic coefficients to be found
384 ROOT::Math::SVector<double, 3> bVec; // cubic coefficients to be found
385 double x0 = 0;
386
387 // cubic interpolation using the two points p0,p1 and the slope at p0 and the second derivative at p0
388
389 // cut toler8 as function goes up
390 double x1 = slam;
391 cubicMatrix(0,0) = (x0*x0*x0 - x1*x1*x1)/3.;
392 cubicMatrix(0,1) = (x0*x0 - x1*x1)/2.;
393 cubicMatrix(0,2) = (x0 - x1);
394 cubicMatrix(1,0) = x0*x0;
395 cubicMatrix(1,1) = x0;
396 cubicMatrix(1,2) = 1;
397 cubicMatrix(2,0) = 2.*x0;
398 cubicMatrix(2,1) = 1;
399 cubicMatrix(2,2) = 0;
400
401 bVec(0) = f0-f1;
402 bVec(1) = gdel;
403 bVec(2) = g2del;
404
405 //if (debug) std::cout << "Matrix:\n " << cubicMatrix << std::endl;
406 if (debug) std::cout << "Vec:\n " << bVec << std::endl;
407
408 // find the minimum need to invert the matrix
409 if (!cubicMatrix.Invert() ) {
410 std::cout << "Inversion failed - return " << std::endl;
411 return MnParabolaPoint(xvmin, fvmin);
412 }
413
414 cubicCoeff = cubicMatrix * bVec;
415 if (debug) std::cout << "Cubic:\n " << cubicCoeff << std::endl;
416
417
418 double ddd = cubicCoeff(1)*cubicCoeff(1)-4*cubicCoeff(0)*cubicCoeff(2); // b**2 - 4ac
419 double slam1, slam2 = 0;
420 // slam1 should be minimum and slam2 the maximum
421 if (cubicCoeff(0) > 0) {
422 slam1 = (- cubicCoeff(1) - std::sqrt(ddd) ) / ( 2. * cubicCoeff(0) );
423 slam2 = (- cubicCoeff(1) + std::sqrt(ddd) ) / ( 2. * cubicCoeff(0) );
424 }
425 else if (cubicCoeff(0) < 0) {
426 slam1 = (- cubicCoeff(1) + std::sqrt(ddd) ) / ( 2. * cubicCoeff(0) );
427 slam2 = (- cubicCoeff(1) - std::sqrt(ddd) ) / ( 2. * cubicCoeff(0) );
428 }
429 else { // case == 0 (-b/c)
430 slam1 = - gdel/g2del;
431 slam2 = slam1;
432 }
433
434 std::cout << "slam1 & slam 2 " << slam1 << " " << slam2 << std::endl;
435
436 // this should be the minimum otherwise inversion failed and I should do something else
437
438 if (slam2 < undral) slam2 = undral;
439 if (slam2 > overal) slam2 = overal;
440
441
442 // I am stack somewhere - take a large step
443 if (std::fabs(slam2) < toler) slam2 = ( slam2 >=0 ) ? slamax : -slamax;
444
445
446 std::cout << "try with slam 2 " << slam2 << std::endl;
447
448 double f2 = fcn(st.Vec()+slam2*step);
449
450 std::cout << "try with slam 2 " << slam2 << " f2 " << f2 << std::endl;
451
452 double fp;
453 // use this as new minimum
454 //bool noImpr = false;
455 if (f2 < fvmin) {
456 slam = slam2;
457 xvmin = slam;
458 fvmin = f2;
459 fp = fvmin;
460 } else {
461 // try with slam2 if it is better
462
463 if (slam1 < undral) slam1 = undral;
464 if (slam1 > overal) slam1 = overal;
465
466
467 if (std::fabs(slam1) < toler) slam1 = ( slam1 >=0 ) ? -slamax : slamax;
468
469 double f3 = fcn(st.Vec()+slam1*step);
470
471 std::cout << "try with slam 1 " << slam1 << " f3 " << f3 << std::endl;
472
473 if (f3 < fvmin) {
474 slam = slam1;
475 fp = fvmin;
476 xvmin = slam;
477 fvmin = f3;
478 } else {
479 // case both f2 and f3 did not produce a better result
480 if (f2 < f3) {
481 slam = slam1;
482 fp = f2;
483 } else {
484 slam = slam2;
485 fp = f3;
486 }
487 }
488 }
489
490 bool iterate2 = false;
491 int niter = 0;
492
493 int maxiter = 10;
494
495 do {
496 iterate2 = false;
497
498 std::cout << "\n iter" << niter << " test approx deriv ad second deriv at " << slam << " fp " << fp << std::endl;
499
500 // estimate grad and second derivative at new point taking a step of 10-3
501 double h = 0.001*slam;
502 double fh = fcn(st.Vec()+(slam+h)*step);
503 double fl = fcn(st.Vec()+(slam-h)*step);
504 double df = (fh-fl)/(2.*h);
505 double df2 = (fh+fl-2.*fp )/(h*h);
506
507 std::cout << "deriv: " << df << " , " << df2 << std::endl;
508
509 // if I am in a point of still negative derivative
510 if ( std::fabs(df ) < prec.Eps() && std::fabs( df2 ) < prec.Eps() ) {
511 // try in opposite direction with an opposite value
512 slam = ( slam >=0 ) ? -slamax : slamax;
513 slamax *= 10;
514 fp = fcn(st.Vec()+slam*step);
515 }
516 else if ( std::fabs(df2 ) <= 0) { // gradient is significative different than zero then redo a cubic interpolation
517 // from new point
518
519 return MnParabolaPoint(slam, fp); // should redo a cubic interpol. ??
520// niter ++;
521// if (niter > maxiter) break;
522
523// MinimumParameters pa = MinimumParameters(st.Vec() + stepNew, fp);
524// gdel = stepNew(i)
525// MnParabolaPoint pp = CubicSearch(fcn, st, stepNew, df, df2
526
527
528 }
529
530 else
531 return MnParabolaPoint(slam, fp);
532
533 niter++;
534 } while (niter < maxiter);
535
536 return MnParabolaPoint(xvmin, fvmin);
537
538
539}
540
541
542// class describing Fcn function in one dimension
543 class ProjectedFcn : public ROOT::Math::IGenFunction {
544
545 public:
546
547 ProjectedFcn(const MnFcn & fcn, const MinimumParameters & pa, const MnAlgebraicVector & step) :
548 fFcn(fcn),
549 fPar(pa),
550 fStep(step)
551 {}
552
553
554 ROOT::Math::IGenFunction * Clone() const { return new ProjectedFcn(*this); }
555
556 private:
557
558 double DoEval(double x) const {
559 return fFcn(fPar.Vec() + x*fStep);
560 }
561
562 const MnFcn & fFcn;
563 const MinimumParameters & fPar;
564 const MnAlgebraicVector & fStep;
565 };
566
567
568MnParabolaPoint MnLineSearch::BrentSearch(const MnFcn& fcn, const MinimumParameters& st, const MnAlgebraicVector& step, double gdel, double g2del, const MnMachinePrecision& prec, bool debug) const {
569
570 if (debug) {
571 std::cout<<"gdel= "<<gdel<<std::endl;
572 std::cout<<"g2del= "<<g2del<<std::endl;
573 for (unsigned int i = 0; i < step.size(); ++i) {
574 if (step(i) != 0) {
575 std::cout<<"step(i)= "<<step(i)<<std::endl;
576 std::cout<<"par(i) " <<st.Vec()(i) <<std::endl;
577 break;
578 }
579 }
580 }
581
582 ProjectedFcn func(fcn, st, step);
583
584
585// do first a cubic interpolation
586
587 double f0 = st.Fval();
588 double f1 = fcn(st.Vec()+step);
589 f0 = func(0.0);
590 f1 = func(1.);
591 double fvmin = st.Fval();
592 double xvmin = 0.;
593 if (debug) std::cout << "f0 " << f0 << " f1 " << f1 << std::endl;
594
595 if(f1 < f0) {
596 fvmin = f1;
597 xvmin = 1.;
598 }
599// double toler8 = toler;
600// double slamax = slambg;
601// double flast = f1;
602 double slam = 1.;
603
604 double undral = -1000;
605 double overal = 1000;
606
607 double x0 = 0;
608
609// MnParabolaPoint p0(0., f0);
610// MnParabolaPoint p1(slam, flast);
611#ifdef USE_CUBIC
612
614 ROOT::Math::SVector<double, 3> cubicCoeff; // cubic coefficients to be found
615 ROOT::Math::SVector<double, 3> bVec; // cubic coefficients to be found
616
617 // cubic interpolation using the two points p0,p1 and the slope at p0 and the second derivative at p0
618
619 // cut toler8 as function goes up
620 double x1 = slam;
621 cubicMatrix(0,0) = (x0*x0*x0 - x1*x1*x1)/3.;
622 cubicMatrix(0,1) = (x0*x0 - x1*x1)/2.;
623 cubicMatrix(0,2) = (x0 - x1);
624 cubicMatrix(1,0) = x0*x0;
625 cubicMatrix(1,1) = x0;
626 cubicMatrix(1,2) = 1;
627 cubicMatrix(2,0) = 2.*x0;
628 cubicMatrix(2,1) = 1;
629 cubicMatrix(2,2) = 0;
630
631 bVec(0) = f0-f1;
632 bVec(1) = gdel;
633 bVec(2) = g2del;
634
635 //if (debug) std::cout << "Matrix:\n " << cubicMatrix << std::endl;
636 if (debug) std::cout << "Vec:\n " << bVec << std::endl;
637
638 // find the minimum need to invert the matrix
639 if (!cubicMatrix.Invert() ) {
640 std::cout << "Inversion failed - return " << std::endl;
641 return MnParabolaPoint(xvmin, fvmin);
642 }
643
644 cubicCoeff = cubicMatrix * bVec;
645 if (debug) std::cout << "Cubic:\n " << cubicCoeff << std::endl;
646
647
648 double ddd = cubicCoeff(1)*cubicCoeff(1)-4*cubicCoeff(0)*cubicCoeff(2); // b**2 - 4ac
649 double slam1, slam2 = 0;
650 // slam1 should be minimum and slam2 the maximum
651 if (cubicCoeff(0) > 0) {
652 slam1 = (- cubicCoeff(1) - std::sqrt(ddd) ) / ( 2. * cubicCoeff(0) );
653 slam2 = (- cubicCoeff(1) + std::sqrt(ddd) ) / ( 2. * cubicCoeff(0) );
654 }
655 else if (cubicCoeff(0) < 0) {
656 slam1 = (- cubicCoeff(1) + std::sqrt(ddd) ) / ( 2. * cubicCoeff(0) );
657 slam2 = (- cubicCoeff(1) - std::sqrt(ddd) ) / ( 2. * cubicCoeff(0) );
658 }
659 else { // case == 0 (-b/c)
660 slam1 = - gdel/g2del;
661 slam2 = slam1;
662 }
663
664 if (slam1 < undral) slam1 = undral;
665 if (slam1 > overal) slam1 = overal;
666
667 if (slam2 < undral) slam2 = undral;
668 if (slam2 > overal) slam2 = overal;
669
670
671 double fs1 = func(slam1);
672 double fs2 = func(slam2);
673 std::cout << "slam1 & slam 2 " << slam1 << " " << slam2 << std::endl;
674 std::cout << "f(slam1) & f(slam 2) " << fs1 << " " << fs2 << std::endl;
675
676 if (fs1 < fs2) {
677 x0 = slam1;
678 f0 = fs1;
679 }
680 else {
681 x0 = slam2;
682 f0 = fs2;
683 }
684
685#else
686 x0 = xvmin;
687 f0 = fvmin;
688#endif
689
690
691 double astart = 100;
692
693
694 // do a Brent search in a large interval
695 double a = x0 -astart;
696 double b = x0 + astart;
697 //double x0 = 1;
698 int maxiter = 20;
699 double absTol = 0.5;
700 double relTol = 0.05;
701
702 ROOT::Math::Minim1D::Type minType = ROOT::Math::Minim1D::BRENT;
703
704 bool iterate = false;
705 int iter = 0;
706
707 std::cout << "f(0)" << func(0.) << std::endl;
708 std::cout << "f(1)" << func(1.0) << std::endl;
709 std::cout << "f(3)" << func(3.0) << std::endl;
710 std::cout << "f(5)" << func(5.0) << std::endl;
711
712 double fa = func(a);
713 double fb = func(b);
714 //double f0 = func(x0);
715 double toler = 0.01;
716 double delta0 = 5;
717 double delta = delta0;
718 bool enlarge = true;
719 bool scanmin = false;
720 double x2, f2 = 0;
721 double dir = 1;
722 int nreset = 0;
723
724 do {
725
726 std::cout << " iter " << iter << " a , b , x0 " << a << " " << b << " x0 " << x0 << std::endl;
727 std::cout << " fa , fb , f0 " << fa << " " << fb << " f0 " << f0 << std::endl;
728 if (fa <= f0 || fb <= f0) {
729 scanmin = false;
730 if (fa < fb) {
731 if (fa < fvmin) {
732 fvmin = fa;
733 xvmin = a;
734 }
735// double f2 = fa;
736// double x2 = a;
737 if (enlarge) {
738 x2 = a - 1000; // move lower
739 f2 = func(x2);
740 }
741 if ( std::fabs ( (fa -f2 )/(a-x2) ) > toler ) { // significant change in f continue to enlarge interval
742 x0 = a;
743 f0 = fa;
744 a = x2;
745 fa = f2;
746 enlarge = true;
747 } else {
748 // move direction of [a
749 // values change a little, start from central point try with x0 = x0 - delta
750 if (nreset == 0) dir = -1;
751 enlarge = false;
752 x0 = x0 + dir * delta;
753 f0 = func(x0);
754 // if reached limit try opposite direction , direction b]
755 if ( std::fabs (( f0 -fa )/(x0 -a) ) < toler ) {
756 delta = 10 * delta0/(10.*(nreset+1)); // decrease the delta if still not change observed
757 a = x0;
758 fa = f0;
759 x0 = delta;
760 f0 = func(x0);
761 dir *= -1;
762 nreset++;
763 std::cout << " A: Done a reset - scan in opposite direction ! " << std::endl;
764 }
765 delta *= 2; // double delta at next iteration
766 if (x0 < b && x0 > a )
767 scanmin = true;
768 else {
769 x0 = 0;
770 f0 = st.Fval();
771 }
772 }
773 }
774 else { // fb < fa
775 if (fb < fvmin) {
776 fvmin = fb;
777 xvmin = b;
778 }
779 if (enlarge) {
780 x2 = b + 1000; // move upper
781 f2 = func(x2);
782 }
783 if ( std::fabs ( (fb -f2 )/(x2-b) ) > toler ) { // significant change in f continue to enlarge interval
784 f0 = fb;
785 x0 = b;
786 b = x2; //
787 fb = f2;
788 enlarge = true;
789 } else {
790 // here move direction b
791 // values change a little, try with x0 = fa + delta
792 if (nreset == 0) dir = 1;
793 enlarge = false;
794 x0 = x0 + dir * delta;
795 f0 = func(x0);
796 // if reached limit try other side
797 if ( std::fabs ( ( f0 -fb )/(x0-b) ) < toler ) {
798 delta = 10 * delta0/(10.*(nreset+1)); // decrease the delta if still not change observed
799 b = x0;
800 fb = f0;
801 x0 = -delta;
802 f0 = func(x0);
803 dir *= -1;
804 nreset++;
805 std::cout << " B: Done a reset - scan in opposite direction ! " << std::endl;
806 }
807 delta *= 2; // double delta at next iteration
808 if (x0 > a && x0 < b)
809 scanmin = true;
810 else {
811 x0 = 0;
812 f0 = st.Fval();
813 }
814
815 }
816
817 }
818
819 if ( f0 < fvmin) {
820 x0 = xvmin;
821 fvmin = f0;
822 }
823
824 std::cout << " new x0 " << x0 << " " << f0 << std::endl;
825
826 // use golden section
827 iterate = scanmin || enlarge;
828
829 } else { // x0 is a min of [a,b]
830 iterate = false;
831 }
832
833 iter++;
834 } while (iterate && iter < maxiter );
835
836 if ( f0 > fa || f0 > fb)
837 // skip minim 1d try Minuit LS
838 // return (*this) (fcn, st, step, gdel, prec, debug);
839 return MnParabolaPoint(xvmin, fvmin);
840
841
842
843
844 std::cout << "1D minimization using " << minType << std::endl;
845
846 ROOT::Math::Minimizer1D min(minType);
847
848 min.SetFunction(func,x0, a, b);
849 int ret = min.Minimize(maxiter, absTol, relTol);
850
851 std::cout << "result of GSL 1D minimization: " << ret << " niter " << min.Iterations() << std::endl;
852 std::cout << " xmin = " << min.XMinimum() << " fmin " << min.FValMinimum() << std::endl;
853
854 return MnParabolaPoint(min.XMinimum(), min.FValMinimum());
855
856}
857
858#endif
859
860
861
862} // namespace Minuit2
863
864} // namespace ROOT
#define b(i)
Definition: RSha256.hxx:100
#define h(i)
Definition: RSha256.hxx:106
static const double x2[5]
static const double x1[5]
double sqrt(double)
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
Definition: IFunction.h:135
SMatrix: a generic fixed size D1 x D2 Matrix class.
Definition: SMatrix.h:124
bool Invert()
Invert a square Matrix ( this method changes the current matrix).
Definition: SMatrix.icc:412
SVector: a generic fixed size Vector class.
Definition: SVector.h:75
unsigned int size() const
Definition: LAVector.h:198
const MnAlgebraicVector & Vec() const
Wrapper class to FCNBase interface used internally by Minuit.
Definition: MnFcn.h:33
MnParabolaPoint operator()(const MnFcn &, const MinimumParameters &, const MnAlgebraicVector &, double, const MnMachinePrecision &, bool debug=false) const
Perform a line search from position defined by the vector st along the direction step,...
Sets the relative floating point (double) arithmetic precision.
double Eps() const
eps returns the smallest possible number so that 1.+eps > 1.
double Eps2() const
eps2 returns 2*sqrt(eps)
A point of a parabola.
double Y() const
Accessor to the y (second) coordinate.
double X() const
Accessor to the x (first) coordinate.
This class defines a parabola of the form a*x*x + b*x + c.
Definition: MnParabola.h:31
double B() const
Accessor to the coefficient of the linear term.
Definition: MnParabola.h:149
double C() const
Accessor to the coefficient of the constant term.
Definition: MnParabola.h:160
double Min() const
Calculates the x coordinate of the Minimum of the parabola.
Definition: MnParabola.h:116
double A() const
Accessor to the coefficient of the quadratic term.
Definition: MnParabola.h:138
Type
Enumeration with One Dimensional Minimizer Algorithms.
Double_t x[n]
Definition: legend1.C:17
TF1 * f1
Definition: legend1.C:11
int iterate(rng_state_t *X)
Definition: mixmax.icc:34
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
LAVector MnAlgebraicVector
Definition: MnMatrix.h:42
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
Definition: StringConv.hxx:21
auto * a
Definition: textangle.C:12