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