Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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#include "Minuit2/MnPrint.h"
29
30namespace ROOT {
31
32namespace Minuit2 {
33
34double 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
#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:227
void TraceIteration(int iter, const MinimumState &state) const
MinimumError keeps the inv.
const MnAlgebraicSymMatrix & InvHessian() const
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:43
const MinimumParameters & Parameters() const
Definition MinimumSeed.h:40
const MnMachinePrecision & Precision() const
Definition MinimumSeed.h:44
const MinimumState & State() const
Definition MinimumSeed.h:39
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
const MinimumError & Error() const
const MnAlgebraicVector & Vec() const
const FunctionGradient & Gradient() const
Wrapper class to FCNBase interface used internally by Minuit.
Definition MnFcn.h:30
double Up() const
Definition MnFcn.cxx:39
unsigned int NumOfCalls() const
Definition MnFcn.h:39
API class for calculating the numerical covariance matrix (== 2x Inverse Hessian == 2x Inverse 2nd de...
Definition MnHesse.h:39
Implements a 1-dimensional minimization along a given direction (i.e.
Sets the relative floating point (double) arithmetic precision.
double Eps() const
eps returns the smallest possible number so that 1.+eps > 1.
double Eps2() const
eps2 returns 2*sqrt(eps)
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:25
void Debug(const Ts &... args)
Definition MnPrint.h:138
void Info(const Ts &... args)
Definition MnPrint.h:132
void Warn(const Ts &... args)
Definition MnPrint.h:126
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
double inner_product(const LAVector &, const LAVector &)
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...