29#include "Math/Minimizer1D.h"
73 for (
unsigned int i = 0; i < step.
size(); i++) {
76 double ratio = std::fabs(
st.Vec()(i) / step(i));
86 double f0 =
st.Fval();
87 double f1 =
fcn(
st.Vec() + step);
130#ifdef TRY_OPPOSIT_DIR
138#ifdef TRY_OPPOSITE_DIR
148 print.
Debug(
"slam larger than mac value - set to",
slamax);
157 print.
Debug(
"slam smaller than",
slamin,
"return");
181#ifdef TRY_OPPOSITE_DIR
213 print.
Debug(
"after initial 2-point iter:",
'\n',
" x0, x1, x2:",
p0.X(),
p1.X(),
slam,
'\n',
" f0, f1, f2:",
p0.Y(),
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()) {
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");
291 if (
p0.Y() >
p1.Y() &&
p0.Y() >
p2.Y())
293 else if (
p1.Y() >
p0.Y() &&
p1.Y() >
p2.Y())
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);
336 for (
unsigned int i = 0; i < step.
size(); i++) {
339 double ratio = std::fabs(
st.Vec()(i) / step(i));
349 double f0 =
st.Fval();
350 double f1 =
fcn(
st.Vec() + step);
353 print.Debug(
"f0",
f0,
"f1",
f1);
391 print.Debug(
"Vec:\n ",
bVec);
395 print.Warn(
"Inversion failed - return");
431 print.Debug(
"try with slam 2",
slam2,
"f2", f2);
454 print.Debug(
"try with slam 1",
slam1,
"f3", f3);
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()) {
498 }
else if (std::fabs(
df2) <= 0) {
501 return MnParabolaPoint(
slam, fp);
512 return MnParabolaPoint(
slam, fp);
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';
560 double f0 =
st.Fval();
561 double f1 =
fcn(
st.Vec() + step);
566 print.Debug(
"f0",
f0,
"f1",
f1);
609 print.Debug(
"Vec:\n ",
bVec);
613 print.Warn(
"Inversion failed - return");
646 print.Debug(
"slam1",
slam1,
"slam2",
slam2,
"f(slam1)",
fs1,
"f(slam2)",
fs2);
676 print.Debug(
"f(0)", func(0.),
"f1", func(1.0),
"f(3)", func(3.0),
"f(5)", func(5.0));
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) {
729 print.Info(
"A: Done a reset - scan in opposite direction!");
748 if (std::fabs((fb - f2) / (
x2 -
b)) >
toler) {
760 x0 = x0 + dir * delta;
763 if (std::fabs((
f0 - fb) / (x0 -
b)) <
toler) {
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);
800 if (
f0 > fa ||
f0 > fb)
805 print.Info(
"1D minimization using",
minType);
807 ROOT::Math::Minimizer1D min(
minType);
809 min.SetFunction(func, x0,
a,
b);
812 MnPrint::info(
"result of GSL 1D minimization:",
ret,
"niter", min.Iterations(),
"xmin", min.XMinimum(),
"fmin",
815 return MnParabolaPoint(min.XMinimum(), min.FValMinimum());
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
Option_t Option_t TPoint TPoint const char x2
Option_t Option_t TPoint TPoint const char x1
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.
SVector: a generic fixed size Vector class.
unsigned int size() 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.
This class defines a parabola of the form a*x*x + b*x + c.
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...