29#include "Math/Minimizer1D.h"
61 print.
Debug(
"gdel", gdel,
"step", step);
63 double overal = 1000.;
64 double undral = -100.;
73 for (
unsigned int i = 0; i < step.
size(); i++) {
76 double ratio = std::fabs(st.
Vec()(i) / step(i));
82 if (std::fabs(slamin) < prec.
Eps())
84 slamin *= prec.
Eps2();
86 double f0 = st.
Fval();
87 double f1 = fcn(st.
Vec() + step);
89 double fvmin = st.
Fval();
96 double toler8 = toler;
97 double slamax = slambg;
111 print.
Debug(
"flast", flast,
"f0", f0,
"flast-f0", flast - f0,
"slam", slam);
119 double denom = 2. * (flast - f0 - gdel * slam) / (slam * slam);
121 print.
Debug(
"denom", denom);
123 slam = -gdel / denom;
128 print.
Debug(
"new slam", slam);
130#ifdef TRY_OPPOSIT_DIR
132 bool slamIsNeg =
false;
137 print.
Debug(
"slam is negative - set to", slamax);
138#ifdef TRY_OPPOSITE_DIR
139 slamNeg = 2 * slam - 1;
141 print.
Debug(
"slam is negative - compare values between", slamNeg,
"and", slamax);
148 print.
Debug(
"slam larger than mac value - set to", slamax);
152 print.
Debug(
"slam too small - set to", toler8);
157 print.
Debug(
"slam smaller than", slamin,
"return");
163 if (std::fabs(slam - 1.) < toler8 && p1.
Y() < p0.
Y()) {
169 if (std::fabs(slam - 1.) < toler8)
177 f2 = fcn(st.
Vec() + slam * step);
181#ifdef TRY_OPPOSITE_DIR
184 double f2alt = fcn(st.
Vec() + slamNeg * step);
197 if (std::fabs(p0.
Y() - fvmin) < std::fabs(fvmin) * prec.
Eps()) {
201 toler8 = toler * slam;
202 overal = slam - toler8;
207 }
while (
iterate && niter < maxiter);
208 if (niter >= maxiter) {
213 print.
Debug(
"after initial 2-point iter:",
'\n',
" x0, x1, x2:", p0.
X(), p1.
X(), slam,
'\n',
" f0, f1, f2:", p0.
Y(),
220 slamax = std::max(slamax, alpha * std::fabs(xvmin));
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();
228 slam = xvmin + slamax;
230 slam = xvmin - slamax;
231 print.
Debug(
"xvmin", xvmin,
"slopem", slopem,
"slam", slam);
235 if (slam > xvmin + slamax)
236 slam = xvmin + slamax;
237 if (slam < xvmin - slamax)
238 slam = xvmin - slamax;
248 print.
Debug(
"slam", slam,
"undral", undral,
"overal", overal);
253 print.
Debug(
"iterate on f3- slam", niter,
"slam", slam,
"xvmin", xvmin);
256 double toler9 = std::max(toler8, std::fabs(toler8 * slam));
258 if (std::fabs(p0.
X() - slam) < toler9 || std::fabs(p1.X() - slam) < toler9 ||
259 std::fabs(p2.
X() - slam) < toler9) {
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());
272 if (f3 > p0.
Y() && f3 > p1.Y() && f3 > p2.
Y()) {
273 print.
Debug(
"f3 worse than all three previous");
275 overal = std::min(overal, slam - toler8);
277 undral = std::max(undral, slam + toler8);
278 slam = 0.5 * (slam + xvmin);
279 print.
Debug(
"new slam", slam);
283 }
while (
iterate && niter < maxiter);
284 if (niter >= maxiter) {
291 if (p0.
Y() > p1.Y() && p0.
Y() > p2.
Y())
293 else if (p1.Y() > p0.
Y() && p1.Y() > p2.
Y())
297 print.
Debug(
"f3", f3,
"fvmin", fvmin,
"xvmin", xvmin);
303 overal = std::min(overal, slam - toler8);
305 undral = std::max(undral, slam + toler8);
309 }
while (niter < maxiter);
311 print.
Debug(
"f1, f2 =", p0.
Y(), p1.
Y(),
'\n',
"x1, x2 =", p0.
X(), p1.
X(),
'\n',
"x, f =", xvmin, fvmin);
324 MnPrint print(
"MnLineSearch::CubicSearch");
326 print.Debug(
"gdel", gdel,
"g2del", g2del,
"step", step);
329 double overal = 100.;
330 double undral = -100.;
336 for (
unsigned int i = 0; i < step.
size(); i++) {
339 double ratio = std::fabs(st.
Vec()(i) / step(i));
345 if (std::fabs(slamin) < prec.
Eps())
347 slamin *= prec.
Eps2();
349 double f0 = st.
Fval();
350 double f1 = fcn(st.
Vec() + step);
351 double fvmin = st.
Fval();
353 print.Debug(
"f0", f0,
"f1",
f1);
359 double toler8 = toler;
360 double slamax = slambg;
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;
391 print.Debug(
"Vec:\n ", bVec);
394 if (!cubicMatrix.
Invert()) {
395 print.Warn(
"Inversion failed - return");
396 return MnParabolaPoint(xvmin, fvmin);
399 cubicCoeff = cubicMatrix * bVec;
400 print.Debug(
"Cubic:\n ", cubicCoeff);
402 double ddd = cubicCoeff(1) * cubicCoeff(1) - 4 * cubicCoeff(0) * cubicCoeff(2);
403 double slam1, slam2 = 0;
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));
412 slam1 = -gdel / g2del;
416 print.Debug(
"slam1", slam1,
"slam2", slam2);
426 if (std::fabs(slam2) < toler)
427 slam2 = (slam2 >= 0) ? slamax : -slamax;
429 double f2 = fcn(st.
Vec() + slam2 * step);
431 print.Debug(
"try with slam 2", slam2,
"f2", f2);
449 if (std::fabs(slam1) < toler)
450 slam1 = (slam1 >= 0) ? -slamax : slamax;
452 double f3 = fcn(st.
Vec() + slam1 * step);
454 print.Debug(
"try with slam 1", slam1,
"f3", f3);
473 bool iterate2 =
false;
481 print.Debug(
"iter", niter,
"test approx deriv ad second deriv at", slam,
"fp", fp);
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);
490 print.Debug(
"deriv", df, df2);
493 if (std::fabs(df) < prec.
Eps() && std::fabs(df2) < prec.
Eps()) {
495 slam = (slam >= 0) ? -slamax : slamax;
497 fp = fcn(st.
Vec() + slam * step);
498 }
else if (std::fabs(df2) <= 0) {
501 return MnParabolaPoint(slam, fp);
512 return MnParabolaPoint(slam, fp);
515 }
while (niter < maxiter);
517 return MnParabolaPoint(xvmin, fvmin);
524 ProjectedFcn(
const MnFcn &fcn,
const MinimumParameters &pa,
const MnAlgebraicVector &step)
525 : fFcn(fcn), fPar(pa), fStep(step)
532 double DoEval(
double x)
const {
return fFcn(fPar.Vec() +
x * fStep); }
535 const MinimumParameters &fPar;
539MnParabolaPoint MnLineSearch::BrentSearch(
const MnFcn &fcn,
const MinimumParameters &st,
const MnAlgebraicVector &step,
540 double gdel,
double g2del,
const MnMachinePrecision &prec)
const
542 MnPrint print(
"MnLineSearch::BrentSearch");
544 print.Debug(
"gdel", gdel,
"g2del", g2del);
546 print.Debug([&](std::ostream &os) {
547 for (
unsigned int i = 0; i < step.size(); ++i) {
549 os <<
"step(i) " << step(i) <<
'\n';
550 std::cout <<
"par(i) " << st.Vec()(i) <<
'\n';
556 ProjectedFcn func(fcn, st, step);
560 double f0 = st.Fval();
561 double f1 = fcn(st.Vec() + step);
564 double fvmin = st.Fval();
566 print.Debug(
"f0", f0,
"f1",
f1);
577 double undral = -1000;
578 double overal = 1000;
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;
609 print.Debug(
"Vec:\n ", bVec);
612 if (!cubicMatrix.
Invert()) {
613 print.Warn(
"Inversion failed - return");
614 return MnParabolaPoint(xvmin, fvmin);
617 cubicCoeff = cubicMatrix * bVec;
618 print.Debug(
"Cubic:\n ", cubicCoeff);
620 double ddd = cubicCoeff(1) * cubicCoeff(1) - 4 * cubicCoeff(0) * cubicCoeff(2);
621 double slam1, slam2 = 0;
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));
630 slam1 = -gdel / g2del;
644 double fs1 = func(slam1);
645 double fs2 = func(slam2);
646 print.Debug(
"slam1", slam1,
"slam2", slam2,
"f(slam1)", fs1,
"f(slam2)", fs2);
664 double a = x0 - astart;
665 double b = x0 + astart;
669 double relTol = 0.05;
676 print.Debug(
"f(0)", func(0.),
"f1", func(1.0),
"f(3)", func(3.0),
"f(5)", func(5.0));
683 double delta = delta0;
685 bool scanmin =
false;
692 print.Debug(
"iter", iter,
"a",
a,
"b",
b,
"x0", x0,
"fa", fa,
"fb", fb,
"f0", f0);
693 if (fa <= f0 || fb <= f0) {
706 if (std::fabs((fa - f2) / (
a -
x2)) > toler) {
718 x0 = x0 + dir * delta;
721 if (std::fabs((f0 - fa) / (x0 -
a)) < toler) {
722 delta = 10 * delta0 / (10. * (nreset + 1));
729 print.Info(
"A: Done a reset - scan in opposite direction!");
732 if (x0 < b && x0 >
a)
748 if (std::fabs((fb - f2) / (
x2 -
b)) > toler) {
760 x0 = x0 + dir * delta;
763 if (std::fabs((f0 - fb) / (x0 -
b)) < toler) {
764 delta = 10 * delta0 / (10. * (nreset + 1));
771 print.Info(
"B: Done a reset - scan in opposite direction!");
774 if (x0 >
a && x0 <
b)
788 print.Debug(
"new x0", x0,
"f0", f0);
798 }
while (
iterate && iter < maxiter);
800 if (f0 > fa || f0 > fb)
803 return MnParabolaPoint(xvmin, fvmin);
805 print.Info(
"1D minimization using", minType);
807 ROOT::Math::Minimizer1D min(minType);
809 min.SetFunction(func, x0,
a,
b);
810 int ret = min.Minimize(maxiter, absTol, relTol);
812 MnPrint::info(
"result of GSL 1D minimization:", ret,
"niter", min.Iterations(),
"xmin", min.XMinimum(),
"fmin",
815 return MnParabolaPoint(min.XMinimum(), min.FValMinimum());
static const double x2[5]
static const double x1[5]
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
SMatrix: a generic fixed size D1 x D2 Matrix class.
bool Invert()
Invert a square Matrix ( this method changes the current matrix).
SVector: a generic fixed size Vector class.
unsigned int size() const
const MnAlgebraicVector & Vec() const
Wrapper class to FCNBase interface used internally by Minuit.
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,...
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)
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.
double B() const
Accessor to the coefficient of the linear term.
double C() const
Accessor to the coefficient of the constant term.
double Min() const
Calculates the x coordinate of the Minimum of the parabola.
double A() const
Accessor to the coefficient of the quadratic term.
void Debug(const Ts &... args)
Type
Enumeration with One Dimensional Minimizer Algorithms.
int iterate(rng_state_t *X)
LAVector MnAlgebraicVector
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...