Logo ROOT  
Reference Guide
FumiliBuilder.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/FumiliBuilder.h"
13 //#include "Minuit2/Numerical2PGradientCalculator.h"
14 #include "Minuit2/MinimumState.h"
15 #include "Minuit2/MinimumError.h"
18 #include "Minuit2/MnLineSearch.h"
19 #include "Minuit2/MinimumSeed.h"
20 #include "Minuit2/MnFcn.h"
22 #include "Minuit2/MnPosDef.h"
24 #include "Minuit2/LaSum.h"
25 #include "Minuit2/LaProd.h"
26 #include "Minuit2/MnStrategy.h"
27 #include "Minuit2/MnHesse.h"
28 #include "Minuit2/MnPrint.h"
29 
30 namespace ROOT {
31 
32 namespace Minuit2 {
33 
34 double inner_product(const LAVector &, const LAVector &);
35 
37  const MnStrategy &strategy, unsigned int maxfcn, double edmval) const
38 {
39  // top level function to find minimum from a given initial seed
40  // iterate on a minimum search in case of first attempt is not successful
41 
42  MnPrint print("FumiliBuilder", PrintLevel());
43 
44  edmval *= 0.0001;
45  // edmval *= 0.1; // use small factor for Fumili
46 
47  print.Debug("Convergence when edm <", edmval);
48 
49  if (seed.Parameters().Vec().size() == 0) {
50  return FunctionMinimum(seed, fcn.Up());
51  }
52 
53  // double edm = Estimator().Estimate(seed.Gradient(), seed.Error());
54  double edm = seed.State().Edm();
55 
56  FunctionMinimum min(seed, fcn.Up());
57 
58  if (edm < 0.) {
59  print.Warn("Initial matrix not pos.def.");
60  return min;
61  }
62 
63  std::vector<MinimumState> result;
64  // result.reserve(1);
65  result.reserve(8);
66 
67  result.push_back(seed.State());
68 
69  print.Info("Start iterating until Edm is <", edmval, '\n', "Initial state", MnPrint::Oneline(seed.State()));
70 
71  if (TraceIter())
72  TraceIteration(result.size() - 1, result.back());
73 
74  // do actual iterations
75 
76  // try first with a maxfxn = 50% of maxfcn
77  // Fumili in principle needs much less iterations
78  int maxfcn_eff = int(0.5 * maxfcn);
79  int ipass = 0;
80  double edmprev = 1;
81 
82  do {
83 
84  min = Minimum(fcn, gc, seed, result, maxfcn_eff, edmval);
85 
86  // second time check for validity of function Minimum
87  if (ipass > 0) {
88  if (!min.IsValid()) {
89 
90  print.Warn("FunctionMinimum is invalid");
91 
92  return min;
93  }
94  }
95 
96  // resulting edm of minimization
97  edm = result.back().Edm();
98 
99  print.Debug("Approximate Edm", edm, "npass", ipass);
100 
101  // call always Hesse (error matrix from Fumili is never accurate since is approximate)
102 
103  print.Debug("FumiliBuilder will verify convergence and Error matrix; "
104  "dcov",
105  min.Error().Dcovar());
106 
107  // // recalculate the gradient using the numerical gradient calculator
108  // Numerical2PGradientCalculator ngc(fcn, min.Seed().Trafo(), strategy);
109  // FunctionGradient ng = ngc( min.State().Parameters() );
110  // MinimumState tmp( min.State().Parameters(), min.State().Error(), ng, min.State().Edm(),
111  // min.State().NFcn() );
112 
113  MinimumState st = MnHesse(strategy)(fcn, min.State(), min.Seed().Trafo(), maxfcn);
114  result.push_back(st);
115 
116  print.Info("After Hessian");
117  if (TraceIter())
118  TraceIteration(result.size() - 1, result.back());
119 
120  // check edm
121  edm = st.Edm();
122 
123  print.Debug("Edm", edm, "State", st);
124 
125  // break the loop if edm is NOT getting smaller
126  if (ipass > 0 && edm >= edmprev) {
127  print.Warn("Stop iterations, no improvements after Hesse; current Edm", edm, "previous value", edmprev);
128  break;
129  }
130  if (edm > edmval) {
131  print.Debug("Tolerance not sufficient, continue minimization; Edm", edm, "Requested", edmval);
132  } else {
133  // Case when edm < edmval after Heasse but min is flagged eith a bad edm:
134  // make then a new Function minimum since now edm is ok
135  if (min.IsAboveMaxEdm()) {
136  min = FunctionMinimum(min.Seed(), min.States(), min.Up());
137  break;
138  }
139  }
140 
141  // end loop on iterations
142  // ? need a maximum here (or max of function calls is enough ? )
143  // continnue iteration (re-calculate function Minimum if edm IS NOT sufficient)
144  // no need to check that hesse calculation is done (if isnot done edm is OK anyway)
145  // count the pass to exit second time when function Minimum is invalid
146  // increase by 20% maxfcn for doing some more tests
147  if (ipass == 0)
148  maxfcn_eff = maxfcn;
149  ipass++;
150  edmprev = edm;
151 
152  } while (edm > edmval);
153 
154  // Add latest state (Hessian calculation)
155  min.Add(result.back());
156 
157  return min;
158 }
159 
161  std::vector<MinimumState> &result, unsigned int maxfcn, double edmval) const
162 {
163 
164  // function performing the minimum searches using the FUMILI algorithm
165  // after the modification when I iterate on this functions, so it can be called many times,
166  // the seed is used here only to get precision and construct the returned FunctionMinimum object
167 
168  /*
169  Three options were possible:
170 
171  1) create two parallel and completely separate hierarchies, in which case
172  the FumiliMinimizer would NOT inherit from ModularFunctionMinimizer,
173  FumiliBuilder would not inherit from MinimumBuilder etc
174 
175  2) Use the inheritance (base classes of ModularFunctionMinimizer,
176  MinimumBuilder etc), but recreate the member functions Minimize() and
177  Minimum() respectively (naming them for example minimize2() and
178  minimum2()) so that they can take FumiliFCNBase as Parameter instead FCNBase
179  (otherwise one wouldn't be able to call the Fumili-specific methods).
180 
181  3) Cast in the daughter classes derived from ModularFunctionMinimizer,
182  MinimumBuilder.
183 
184  The first two would mean to duplicate all the functionality already existent,
185  which is a very bad practice and Error-prone. The third one is the most
186  elegant and effective solution, where the only constraint is that the user
187  must know that they have to pass a subclass of FumiliFCNBase to the FumiliMinimizer
188  and not just a subclass of FCNBase.
189  BTW, the first two solutions would have meant to recreate also a parallel
190  structure for MnFcn...
191  **/
192  // const FumiliFCNBase* tmpfcn = dynamic_cast<const FumiliFCNBase*>(&(fcn.Fcn()));
193 
194  const MnMachinePrecision &prec = seed.Precision();
195 
196  const MinimumState &initialState = result.back();
197 
198  double edm = initialState.Edm();
199 
200  MnPrint print("FumiliBuilder");
201 
202  print.Debug("Initial State:", "\n Parameter", initialState.Vec(), "\n Gradient", initialState.Gradient().Vec(),
203  "\n Inv Hessian", initialState.Error().InvHessian(), "\n edm", initialState.Edm(), "\n maxfcn",
204  maxfcn, "\n tolerance", edmval);
205 
206  // iterate until edm is small enough or max # of iterations reached
207  edm *= (1. + 3. * initialState.Error().Dcovar());
208  MnAlgebraicVector step(initialState.Gradient().Vec().size());
209 
210  // initial lambda Value
211  double lambda = 0.001;
212 
213  do {
214 
215  // const MinimumState& s0 = result.back();
216  MinimumState s0 = result.back();
217 
218  step = -1. * s0.Error().InvHessian() * s0.Gradient().Vec();
219 
220  print.Debug("Iteration -", result.size(), "\n Fval", s0.Fval(), "numOfCall", fcn.NumOfCalls(),
221  "\n Internal Parameter values", s0.Vec(), "\n Newton step", step);
222 
223  double gdel = inner_product(step, s0.Gradient().Grad());
224  if (gdel > 0.) {
225  print.Warn("Matrix not pos.def, gdel =", gdel, " > 0");
226 
227  MnPosDef psdf;
228  s0 = psdf(s0, prec);
229  step = -1. * s0.Error().InvHessian() * s0.Gradient().Vec();
230  gdel = inner_product(step, s0.Gradient().Grad());
231 
232  print.Warn("After correction, gdel =", gdel);
233 
234  if (gdel > 0.) {
235  result.push_back(s0);
236  return FunctionMinimum(seed, result, fcn.Up());
237  }
238  }
239 
240  // MnParabolaPoint pp = lsearch(fcn, s0.Parameters(), step, gdel, prec);
241 
242  // if(std::fabs(pp.Y() - s0.Fval()) < prec.Eps()) {
243  // std::cout<<"FumiliBuilder: no improvement"<<std::endl;
244  // break; //no improvement
245  // }
246 
247  // MinimumParameters p(s0.Vec() + pp.X()*step, pp.Y());
248 
249  // if taking a full step
250 
251  // take a full step
252 
253  MinimumParameters p(s0.Vec() + step, fcn(s0.Vec() + step));
254 
255  // check that taking the full step does not deteriorate minimum
256  // in that case do a line search
257  if (p.Fval() >= s0.Fval()) {
258  MnLineSearch lsearch;
259  MnParabolaPoint pp = lsearch(fcn, s0.Parameters(), step, gdel, prec);
260 
261  if (std::fabs(pp.Y() - s0.Fval()) < prec.Eps()) {
262  // std::cout<<"FumiliBuilder: no improvement"<<std::endl;
263  break; // no improvement
264  }
265  p = MinimumParameters(s0.Vec() + pp.X() * step, pp.Y());
266  }
267 
268  print.Debug("Before Gradient", fcn.NumOfCalls());
269 
270  FunctionGradient g = gc(p, s0.Gradient());
271 
272  print.Debug("After Gradient", fcn.NumOfCalls());
273 
274  // FunctionGradient g = gc(s0.Parameters(), s0.Gradient());
275 
276  // move Error updator after Gradient since the Value is cached inside
277 
278  MinimumError e = ErrorUpdator().Update(s0, p, gc, lambda);
279 
280  edm = Estimator().Estimate(g, s0.Error());
281 
282  print.Debug("Updated new point:", "\n FVAL ", p.Fval(), "\n Parameter", p.Vec(), "\n Gradient", g.Vec(),
283  "\n InvHessian", e.Matrix(), "\n Hessian", e.Hessian(), "\n Edm", edm);
284 
285  if (edm < 0.) {
286  print.Warn("Matrix not pos.def., Edm < 0");
287 
288  MnPosDef psdf;
289  s0 = psdf(s0, prec);
290  edm = Estimator().Estimate(g, s0.Error());
291  if (edm < 0.) {
292  result.push_back(s0);
293  if (TraceIter())
294  TraceIteration(result.size() - 1, result.back());
295  return FunctionMinimum(seed, result, fcn.Up());
296  }
297  }
298 
299  // check lambda according to step
300  if (p.Fval() < s0.Fval())
301  // fcn is decreasing along the step
302  lambda *= 0.1;
303  else {
304  lambda *= 10;
305  // if close to minimum stop to avoid oscillations around minimum value
306  if (edm < 0.1)
307  break;
308  }
309 
310  print.Debug("finish iteration -", result.size(), "lambda =", lambda, "f1 =", p.Fval(), "f0 =", s0.Fval(),
311  "num of calls =", fcn.NumOfCalls(), "edm =", edm);
312 
313  result.push_back(MinimumState(p, e, g, edm, fcn.NumOfCalls()));
314  if (TraceIter())
315  TraceIteration(result.size() - 1, result.back());
316 
317  print.Info(MnPrint::Oneline(result.back(), result.size()));
318 
319  edm *= (1. + 3. * e.Dcovar());
320 
321  } while (edm > edmval && fcn.NumOfCalls() < maxfcn);
322 
323  if (fcn.NumOfCalls() >= maxfcn) {
324  print.Warn("Call limit exceeded");
325 
326  return FunctionMinimum(seed, result, fcn.Up(), FunctionMinimum::MnReachedCallLimit());
327  }
328 
329  if (edm > edmval) {
330  if (edm < std::fabs(prec.Eps2() * result.back().Fval())) {
331  print.Warn("Machine accuracy limits further improvement");
332 
333  return FunctionMinimum(seed, result, fcn.Up());
334  } else if (edm < 10 * edmval) {
335  return FunctionMinimum(seed, result, fcn.Up());
336  } else {
337 
338  print.Warn("No convergence; Edm", edm, "is above tolerance", 10 * edmval);
339 
340  return FunctionMinimum(seed, result, fcn.Up(), FunctionMinimum::MnAboveMaxEdm());
341  }
342  }
343  // std::cout<<"result.back().Error().Dcovar()= "<<result.back().Error().Dcovar()<<std::endl;
344 
345  print.Debug("Exiting successfully", "Ncalls", fcn.NumOfCalls(), "FCN", result.back().Fval(), "Edm", edm, "Requested",
346  edmval);
347 
348  return FunctionMinimum(seed, result, fcn.Up());
349 }
350 
351 } // namespace Minuit2
352 
353 } // namespace ROOT
ROOT::Minuit2::MnParabolaPoint
A point of a parabola.
Definition: MnParabolaPoint.h:36
ROOT::Minuit2::MnMachinePrecision::Eps2
double Eps2() const
eps2 returns 2*sqrt(eps)
Definition: MnMachinePrecision.h:41
ROOT::Minuit2::MinimumBuilder::PrintLevel
int PrintLevel() const
Definition: MinimumBuilder.h:38
ROOT::Minuit2::MinimumState::Edm
double Edm() const
Definition: MinimumState.h:57
ROOT::Minuit2::MnFcn
Wrapper class to FCNBase interface used internally by Minuit.
Definition: MnFcn.h:30
e
#define e(i)
Definition: RSha256.hxx:103
GradientCalculator.h
ROOT::Minuit2::MnFcn::NumOfCalls
unsigned int NumOfCalls() const
Definition: MnFcn.h:39
ROOT::Minuit2::FunctionMinimum::Up
double Up() const
Definition: FunctionMinimum.h:84
ROOT::Minuit2::MinimumState
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
Definition: MinimumState.h:29
MnMachinePrecision.h
ROOT::Minuit2::LAVector
Definition: LAVector.h:32
ROOT::Minuit2::MnPrint::Debug
void Debug(const Ts &... args)
Definition: MnPrint.h:138
ROOT::Minuit2::MinimumState::Gradient
const FunctionGradient & Gradient() const
Definition: MinimumState.h:55
MnStrategy.h
ROOT::Minuit2::FunctionMinimum::IsValid
bool IsValid() const
Definition: FunctionMinimum.h:85
ROOT::Minuit2::MnPosDef
Force the covariance matrix to be positive defined by adding extra terms in the diagonal.
Definition: MnPosDef.h:25
ROOT::Minuit2::FunctionMinimum::Seed
const MinimumSeed & Seed() const
Definition: FunctionMinimum.h:67
FunctionMinimum.h
FumiliStandardMaximumLikelihoodFCN.h
MinimumSeed.h
ROOT::Minuit2::MinimumParameters::Vec
const MnAlgebraicVector & Vec() const
Definition: MinimumParameters.h:38
ROOT::Minuit2::FumiliBuilder::Minimum
virtual FunctionMinimum Minimum(const MnFcn &fMnFcn, const GradientCalculator &fGradienCalculator, const MinimumSeed &fMinimumSeed, const MnStrategy &fMnStrategy, unsigned int maxfcn, double edmval) const
Class the member function calculating the Minimum and verifies the result depending on the strategy.
Definition: FumiliBuilder.cxx:36
ROOT::Minuit2::MnPrint::Oneline
Definition: MnPrint.h:85
ROOT::Minuit2::FunctionGradient::Vec
const MnAlgebraicVector & Vec() const
Definition: FunctionGradient.h:41
ROOT::Minuit2::FumiliBuilder::ErrorUpdator
const FumiliErrorUpdator & ErrorUpdator() const
Accessor to the Error updator of the builder.
Definition: FumiliBuilder.h:129
ROOT::Minuit2::FunctionGradient
Definition: FunctionGradient.h:21
MnFcn.h
ROOT::Minuit2::MinimumParameters::Fval
double Fval() const
Definition: MinimumParameters.h:40
LaProd.h
ROOT::Minuit2::MnMachinePrecision::Eps
double Eps() const
eps returns the smallest possible number so that 1.+eps > 1.
Definition: MnMachinePrecision.h:38
ROOT::Minuit2::FunctionMinimum::Error
const MinimumError & Error() const
Definition: FunctionMinimum.h:78
ROOT::Minuit2::MnPrint::Warn
void Warn(const Ts &... args)
Definition: MnPrint.h:126
s0
#define s0(x)
Definition: RSha256.hxx:90
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
FunctionGradient.h
ROOT::Minuit2::MinimumState::Error
const MinimumError & Error() const
Definition: MinimumState.h:54
ROOT::Minuit2::MnFcn::Up
double Up() const
Definition: MnFcn.cxx:39
ROOT::Minuit2::inner_product
double inner_product(const LAVector &, const LAVector &)
Definition: LaInnerProduct.cxx:18
ROOT::Minuit2::MinimumBuilder::TraceIteration
void TraceIteration(int iter, const MinimumState &state) const
Definition: MinimumBuilder.h:49
ROOT::Minuit2::FunctionMinimum::Add
void Add(const MinimumState &state)
Definition: FunctionMinimum.h:62
ROOT::Minuit2::FumiliErrorUpdator::Update
virtual MinimumError Update(const MinimumState &fMinimumState, const MinimumParameters &fMinimumParameters, const GradientCalculator &fGradientCalculator, double lambda) const
Member function that calculates the Error matrix (or the Hessian matrix containing the (approximate) ...
Definition: FumiliErrorUpdator.cxx:43
FumiliBuilder.h
ROOT::Minuit2::MinimumError::Dcovar
double Dcovar() const
Definition: MinimumError.h:65
ROOT::Minuit2::MnMachinePrecision
Sets the relative floating point (double) arithmetic precision.
Definition: MnMachinePrecision.h:32
ROOT::Minuit2::MinimumSeed::Precision
const MnMachinePrecision & Precision() const
Definition: MinimumSeed.h:44
ROOT::Minuit2::VariableMetricEDMEstimator::Estimate
double Estimate(const FunctionGradient &, const MinimumError &) const
Definition: VariableMetricEDMEstimator.cxx:20
MnPosDef.h
MnHesse.h
ROOT::Minuit2::MnLineSearch
Implements a 1-dimensional minimization along a given direction (i.e.
Definition: MnLineSearch.h:40
ROOT::Minuit2::MinimumState::Vec
const MnAlgebraicVector & Vec() const
Definition: MinimumState.h:51
ROOT::Minuit2::MnParabolaPoint::Y
double Y() const
Accessor to the y (second) coordinate.
Definition: MnParabolaPoint.h:70
ROOT::Minuit2::FunctionMinimum
class holding the full result of the minimization; both internal and external (MnUserParameterState) ...
Definition: FunctionMinimum.h:33
ROOT::Minuit2::LAVector::size
unsigned int size() const
Definition: LAVector.h:227
ROOT::Minuit2::BasicFunctionMinimum::MnAboveMaxEdm
Definition: BasicFunctionMinimum.h:38
ROOT::Minuit2::MinimumError
MinimumError keeps the inv.
Definition: MinimumError.h:26
ROOT::Minuit2::FumiliBuilder::Estimator
const VariableMetricEDMEstimator & Estimator() const
Accessor to the EDM (expected vertical distance to the Minimum) estimator.
Definition: FumiliBuilder.h:119
MinimumError.h
ROOT::Minuit2::MnStrategy
API class for defining three levels of strategies: low (0), medium (1), high (>=2); acts on: Migrad (...
Definition: MnStrategy.h:27
MnLineSearch.h
ROOT::Minuit2::MinimumSeed::Parameters
const MinimumParameters & Parameters() const
Definition: MinimumSeed.h:40
ROOT::Minuit2::MinimumBuilder::TraceIter
bool TraceIter() const
Definition: MinimumBuilder.h:40
ROOT::Minuit2::MinimumError::InvHessian
const MnAlgebraicSymMatrix & InvHessian() const
Definition: MinimumError.h:61
ROOT::Minuit2::FunctionMinimum::IsAboveMaxEdm
bool IsAboveMaxEdm() const
Definition: FunctionMinimum.h:93
ROOT::Minuit2::FunctionMinimum::State
const MinimumState & State() const
Definition: FunctionMinimum.h:76
LaSum.h
ROOT::Minuit2::MnParabolaPoint::X
double X() const
Accessor to the x (first) coordinate.
Definition: MnParabolaPoint.h:60
MnPrint.h
ROOT::Minuit2::GradientCalculator
interface class for gradient calculators
Definition: GradientCalculator.h:23
ROOT::Minuit2::FunctionMinimum::States
const std::vector< ROOT::Minuit2::MinimumState > & States() const
Definition: FunctionMinimum.h:68
ROOT::Minuit2::BasicFunctionMinimum::MnReachedCallLimit
Definition: BasicFunctionMinimum.h:36
ROOT::Minuit2::MnHesse
API class for calculating the numerical covariance matrix (== 2x Inverse Hessian == 2x Inverse 2nd de...
Definition: MnHesse.h:39
ROOT
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
Definition: EExecutionPolicy.hxx:4
MinimumState.h
ROOT::Minuit2::MinimumSeed
MinimumSeed contains the starting values for the minimization produced by the SeedGenerator.
Definition: MinimumSeed.h:31
ROOT::Minuit2::MnPrint::Info
void Info(const Ts &... args)
Definition: MnPrint.h:132
ROOT::Minuit2::MinimumSeed::State
const MinimumState & State() const
Definition: MinimumSeed.h:39
ROOT::Minuit2::MnPrint
Definition: MnPrint.h:73
ROOT::Minuit2::MinimumSeed::Trafo
const MnUserTransformation & Trafo() const
Definition: MinimumSeed.h:43
int
g
#define g(i)
Definition: RSha256.hxx:105