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 
10 #include "Minuit2/MnLineSearch.h"
11 #include "Minuit2/MnFcn.h"
14 #include "Minuit2/MnParabola.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 
33 namespace ROOT {
34 
35 namespace 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 falg to control the debug
44 */
45 
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 overal = 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  overal = slam - toler8;
203  slamax = overal;
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 > overal)
242  slam = overal;
243  } else {
244  if (slam < undral)
245  slam = undral;
246  }
247 
248  print.Debug("slam", slam, "undral", undral, "overal", overal);
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  overal = std::min(overal, 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  overal = std::min(overal, 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 
321 MnParabolaPoint 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 ot large values
329  double overal = 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 > overal)
423  slam2 = overal;
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 > overal)
447  slam1 = overal;
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
521 class ProjectedFcn : public ROOT::Math::IGenFunction {
522 
523 public:
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 
531 private:
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 
539 MnParabolaPoint 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 overal = 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 > overal)
637  slam1 = overal;
638 
639  if (slam2 < undral)
640  slam2 = undral;
641  if (slam2 > overal)
642  slam2 = overal;
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
ROOT::Minuit2::MnParabolaPoint
A point of a parabola.
Definition: MnParabolaPoint.h:36
ROOT::Minuit2::MnParabola::A
double A() const
Accessor to the coefficient of the quadratic term.
Definition: MnParabola.h:127
ROOT::Minuit2::MnMachinePrecision::Eps2
double Eps2() const
eps2 returns 2*sqrt(eps)
Definition: MnMachinePrecision.h:41
ROOT::Minuit2::MnFcn
Wrapper class to FCNBase interface used internally by Minuit.
Definition: MnFcn.h:30
IFunction.h
ROOT::Math::Minim1D::Type
Type
Enumeration with One Dimensional Minimizer Algorithms.
Definition: GSLMinimizer1D.h:56
MnMachinePrecision.h
ROOT::Minuit2::LAVector
Definition: LAVector.h:32
ROOT::Minuit2::MnPrint::Debug
void Debug(const Ts &... args)
Definition: MnPrint.h:138
ROOT::Math::SVector
SVector: a generic fixed size Vector class.
Definition: SVector.h:75
ROOT::Minuit2::MinimumParameters::Vec
const MnAlgebraicVector & Vec() const
Definition: MinimumParameters.h:38
ROOT::Minuit2::MnParabola::C
double C() const
Accessor to the coefficient of the constant term.
Definition: MnParabola.h:147
ROOT::Minuit2::MnParabola::Min
double Min() const
Calculates the x coordinate of the Minimum of the parabola.
Definition: MnParabola.h:107
ROOT::Math::SMatrix::Invert
bool Invert()
Invert a square Matrix ( this method changes the current matrix).
Definition: SMatrix.icc:412
x
Double_t x[n]
Definition: legend1.C:17
SMatrix.h
b
#define b(i)
Definition: RSha256.hxx:100
MnFcn.h
ROOT::Minuit2::MinimumParameters::Fval
double Fval() const
Definition: MinimumParameters.h:40
x1
static const double x1[5]
Definition: RooGaussKronrodIntegrator1D.cxx:346
ROOT::Minuit2::MnMachinePrecision::Eps
double Eps() const
eps returns the smallest possible number so that 1.+eps > 1.
Definition: MnMachinePrecision.h:38
ROOT::Minuit2::MnAlgebraicVector
LAVector MnAlgebraicVector
Definition: MnMatrix.h:28
SVector.h
MnParabolaPoint.h
ROOT::Math::fabs
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
Definition: UnaryOperators.h:131
ROOT::Minuit2::MinimumParameters
Definition: MinimumParameters.h:21
ROOT::Math::SMatrix
SMatrix: a generic fixed size D1 x D2 Matrix class.
Definition: SMatrix.h:101
h
#define h(i)
Definition: RSha256.hxx:106
ROOT::Minuit2::MnParabolaFactory
Definition: MnParabolaFactory.h:20
a
auto * a
Definition: textangle.C:12
ROOT::Minuit2::MnMachinePrecision
Sets the relative floating point (double) arithmetic precision.
Definition: MnMachinePrecision.h:32
sqrt
double sqrt(double)
MnParabolaFactory.h
ROOT::Minuit2::MnParabola::B
double B() const
Accessor to the coefficient of the linear term.
Definition: MnParabola.h:137
ROOT::Minuit2::MnParabolaPoint::Y
double Y() const
Accessor to the y (second) coordinate.
Definition: MnParabolaPoint.h:70
ROOT::Math::IBaseFunctionOneDim
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
Definition: IFunction.h:135
f1
TF1 * f1
Definition: legend1.C:11
ROOT::Minuit2::LAVector::size
unsigned int size() const
Definition: LAVector.h:227
ROOT::Minuit2::MnParabola
This class defines a parabola of the form a*x*x + b*x + c.
Definition: MnParabola.h:30
MnParabola.h
MnLineSearch.h
x2
static const double x2[5]
Definition: RooGaussKronrodIntegrator1D.cxx:364
MinimumParameters.h
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
LaSum.h
ROOT::Minuit2::MnParabolaPoint::X
double X() const
Accessor to the x (first) coordinate.
Definition: MnParabolaPoint.h:60
MnPrint.h
iterate
int iterate(rng_state_t *X)
Definition: mixmax.icc:34
ROOT
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
Definition: EExecutionPolicy.hxx:4
ROOT::Minuit2::MnPrint
Definition: MnPrint.h:73