Logo ROOT  

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