Logo ROOT   6.18/05
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
13//#include "Minuit2/Numerical2PGradientCalculator.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
38namespace ROOT {
39
40 namespace Minuit2 {
41
42
43
44
45double inner_product(const LAVector&, const LAVector&);
46
47FunctionMinimum 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
197FunctionMinimum 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
#define MN_INFO_MSG(str)
Definition: MnPrint.h:110
#define MN_INFO_VAL(x)
Definition: MnPrint.h:116
#define g(i)
Definition: RSha256.hxx:105
#define s0(x)
Definition: RSha256.hxx:90
#define e(i)
Definition: RSha256.hxx:103
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.
const FumiliErrorUpdator & ErrorUpdator() const
Accessor to the Error updator of the builder.
const VariableMetricEDMEstimator & Estimator() const
Accessor to the EDM (expected vertical distance to the Minimum) estimator.
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) ...
const MnAlgebraicVector & Vec() const
class holding the full result of the minimization; both internal and external (MnUserParameterState) ...
void Add(const MinimumState &state)
const std::vector< ROOT::Minuit2::MinimumState > & States() const
const MinimumError & Error() const
const MinimumState & State() const
const MinimumSeed & Seed() const
interface class for gradient calculators
unsigned int size() const
Definition: LAVector.h:198
void TraceIteration(int iter, const MinimumState &state) const
MinimumError keeps the inv.
Definition: MinimumError.h:26
const MnAlgebraicSymMatrix & InvHessian() const
Definition: MinimumError.h:60
const MnAlgebraicVector & Vec() const
MinimumSeed contains the starting values for the minimization produced by the SeedGenerator.
Definition: MinimumSeed.h:31
const MnUserTransformation & Trafo() const
Definition: MinimumSeed.h:50
const MinimumParameters & Parameters() const
Definition: MinimumSeed.h:47
const MnMachinePrecision & Precision() const
Definition: MinimumSeed.h:51
const MinimumState & State() const
Definition: MinimumSeed.h:46
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
Definition: MinimumState.h:29
const MinimumError & Error() const
Definition: MinimumState.h:62
const MnAlgebraicVector & Vec() const
Definition: MinimumState.h:59
const FunctionGradient & Gradient() const
Definition: MinimumState.h:63
Wrapper class to FCNBase interface used internally by Minuit.
Definition: MnFcn.h:33
double Up() const
Definition: MnFcn.cxx:35
unsigned int NumOfCalls() const
Definition: MnFcn.h:43
API class for calculating the numerical covariance matrix (== 2x Inverse Hessian == 2x Inverse 2nd de...
Definition: MnHesse.h:40
Implements a 1-dimensional minimization along a given direction (i.e.
Definition: MnLineSearch.h:47
determines the relative floating point arithmetic precision.
double Eps() const
eps returns the smallest possible number so that 1.+eps > 1.
double Eps2() const
eps2 returns 2*sqrt(eps)
A point of a parabola.
double Y() const
Accessor to the y (second) coordinate.
double X() const
Accessor to the x (first) coordinate.
Force the covariance matrix to be positive defined by adding extra terms in the diagonal.
Definition: MnPosDef.h:26
static void PrintState(std::ostream &os, const MinimumState &state, const char *msg, int iter=-1)
Definition: MnPrint.cxx:58
static int Level()
Definition: MnPrint.cxx:47
API class for defining three levels of strategies: low (0), medium (1), high (>=2); acts on: Migrad (...
Definition: MnStrategy.h:27
double Estimate(const FunctionGradient &, const MinimumError &) const
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
double inner_product(const LAVector &, const LAVector &)
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21