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