27#include "Math/Minimizer1D.h"
59 print.
Debug(
"doing line search along given step and gdel = s^t*g = ",
gdel);
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();
130#ifdef TRY_OPPOSIT_DIR
138#ifdef TRY_OPPOSITE_DIR
147 print.
Trace(
"slam larger than max value - set to",
slamax);
155 print.
Trace(
"slam smaller than",
slamin,
"return");
179#ifdef TRY_OPPOSITE_DIR
211 print.
Trace(
"after initial 2-point iter:",
'\n',
" x0, x1, x2:",
p0.X(),
p1.X(),
slam,
'\n',
" f0, f1, f2:",
p0.Y(),
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()) {
268 print.
Trace(
"f3", f3,
"f3-p(2-0).Y()", f3 -
p2.Y(), f3 -
p1.Y(), f3 -
p0.Y());
270 if (f3 >
p0.Y() && f3 >
p1.Y() && f3 >
p2.Y()) {
271 print.
Trace(
"f3 worse than all three previous");
284 print.
Debug(
"Exhausted max number of iterations ",
maxiter,
" return");
290 if (
p0.Y() >
p1.Y() &&
p0.Y() >
p2.Y())
292 else if (
p1.Y() >
p0.Y() &&
p1.Y() >
p2.Y())
310 print.
Debug(
"f1, f2 =",
p0.Y(),
p1.Y(),
'\n',
"x1, x2 =",
p0.X(),
p1.X(),
'\n',
"x, f =",
xvmin,
fvmin);
323 MnPrint print(
"MnLineSearch::CubicSearch");
325 print.Trace(
"gdel",
gdel,
"g2del",
g2del,
"step", step);
335 for (
unsigned int i = 0; i < step.
size(); i++) {
338 double ratio = std::fabs(
st.Vec()(i) / step(i));
348 double f0 =
st.Fval();
349 double f1 =
fcn(
st.Vec() + step);
352 print.Debug(
"f0",
f0,
"f1",
f1);
390 print.Debug(
"Vec:\n ",
bVec);
394 print.Warn(
"Inversion failed - return");
430 print.Debug(
"try with slam 2",
slam2,
"f2", f2);
453 print.Debug(
"try with slam 1",
slam1,
"f3", f3);
480 print.Debug(
"iter",
niter,
"test approx deriv ad second deriv at",
slam,
"fp", fp);
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);
489 print.Debug(
"deriv", df,
df2);
492 if (std::fabs(df) <
prec.Eps() && std::fabs(
df2) <
prec.Eps()) {
497 }
else if (std::fabs(
df2) <= 0) {
523 ProjectedFcn(
const MnFcn &
fcn,
const MinimumParameters &
pa,
const MnAlgebraicVector &step)
524 : fFcn(
fcn), fPar(
pa), fStep(step)
531 double DoEval(
double x)
const {
return fFcn(fPar.Vec() +
x * fStep); }
534 const MinimumParameters &fPar;
538MnPoint MnLineSearch::BrentSearch(
const MnFcn &
fcn,
const MinimumParameters &
st,
const MnAlgebraicVector &step,
539 double gdel,
double g2del,
const MnMachinePrecision &
prec)
const
541 MnPrint print(
"MnLineSearch::BrentSearch");
543 print.Debug(
"gdel",
gdel,
"g2del",
g2del);
545 print.Debug([&](std::ostream &os) {
546 for (
unsigned int i = 0; i < step.size(); ++i) {
548 os <<
"step(i) " << step(i) <<
'\n';
549 std::cout <<
"par(i) " <<
st.Vec()(i) <<
'\n';
559 double f0 =
st.Fval();
560 double f1 =
fcn(
st.Vec() + step);
565 print.Debug(
"f0",
f0,
"f1",
f1);
608 print.Debug(
"Vec:\n ",
bVec);
612 print.Warn(
"Inversion failed - return");
645 print.Debug(
"slam1",
slam1,
"slam2",
slam2,
"f(slam1)",
fs1,
"f(slam2)",
fs2);
675 print.Debug(
"f(0)", func(0.),
"f1", func(1.0),
"f(3)", func(3.0),
"f(5)", func(5.0));
691 print.Debug(
"iter", iter,
"a",
a,
"b",
b,
"x0", x0,
"fa", fa,
"fb", fb,
"f0",
f0);
692 if (fa <=
f0 || fb <=
f0) {
705 if (std::fabs((fa - f2) / (
a -
x2)) >
toler) {
717 x0 = x0 + dir * delta;
720 if (std::fabs((
f0 - fa) / (x0 -
a)) <
toler) {
728 print.Info(
"A: Done a reset - scan in opposite direction!");
747 if (std::fabs((fb - f2) / (
x2 -
b)) >
toler) {
759 x0 = x0 + dir * delta;
762 if (std::fabs((
f0 - fb) / (x0 -
b)) <
toler) {
770 print.Info(
"B: Done a reset - scan in opposite direction!");
773 if (x0 >
a && x0 <
b)
787 print.Debug(
"new x0", x0,
"f0",
f0);
799 if (
f0 > fa ||
f0 > fb)
804 print.Info(
"1D minimization using",
minType);
806 ROOT::Math::Minimizer1D min(
minType);
808 min.SetFunction(func, x0,
a,
b);
811 MnPrint::info(
"result of GSL 1D minimization:",
ret,
"niter", min.Iterations(),
"xmin", min.XMinimum(),
"fmin",
814 return {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.
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,...
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)
void Trace(const Ts &... args)
Type
Enumeration with One Dimensional Minimizer Algorithms.
int iterate(rng_state_t *X)
LAVector MnAlgebraicVector
Namespace for new ROOT classes and functions.