ROOT   Reference Guide
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
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 with 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)
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());
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
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();
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 - NCalls = ", fcn.NumOfCalls());
269
271
272 print.Debug("After Gradient - NCalls = ", fcn.NumOfCalls());
273
274
275 // move Error updator after Gradient since the Value is cached inside
276
277 MinimumError e = ErrorUpdator().Update(s0, p, gc, lambda);
278
279 edm = Estimator().Estimate(g, s0.Error());
280
281 print.Debug("Updated new point:", "\n FVAL ", p.Fval(), "\n Parameter", p.Vec(), "\n Gradient", g.Vec(),
282 "\n InvHessian", e.InvHessian(), "\n Hessian", e.Hessian(), "\n Edm", edm);
283
284 if (edm < 0.) {
285 print.Warn("Matrix not pos.def., Edm < 0");
286
287 MnPosDef psdf;
288 s0 = psdf(s0, prec);
289 edm = Estimator().Estimate(g, s0.Error());
290 if (edm < 0.) {
291 result.push_back(s0);
292 if (TraceIter())
293 TraceIteration(result.size() - 1, result.back());
294 return FunctionMinimum(seed, result, fcn.Up());
295 }
296 }
297
298 // check lambda according to step
299 if (p.Fval() < s0.Fval())
300 // fcn is decreasing along the step
301 lambda *= 0.1;
302 else {
303 lambda *= 10;
304 // if close to minimum stop to avoid oscillations around minimum value
305 if (edm < 0.1)
306 break;
307 }
308
309 print.Debug("finish iteration -", result.size(), "lambda =", lambda, "f1 =", p.Fval(), "f0 =", s0.Fval(),
310 "num of calls =", fcn.NumOfCalls(), "edm =", edm);
311
312 result.push_back(MinimumState(p, e, g, edm, fcn.NumOfCalls()));
313 if (TraceIter())
314 TraceIteration(result.size() - 1, result.back());
315
316 print.Info(MnPrint::Oneline(result.back(), result.size()));
317
318 edm *= (1. + 3. * e.Dcovar());
319
320 } while (edm > edmval && fcn.NumOfCalls() < maxfcn);
321
322 if (fcn.NumOfCalls() >= maxfcn) {
323 print.Warn("Call limit exceeded");
324
326 }
327
328 if (edm > edmval) {
329 if (edm < std::fabs(prec.Eps2() * result.back().Fval())) {
330 print.Warn("Machine accuracy limits further improvement");
331
332 return FunctionMinimum(seed, result, fcn.Up());
333 } else if (edm < 10 * edmval) {
334 return FunctionMinimum(seed, result, fcn.Up());
335 } else {
336
337 print.Warn("No convergence; Edm", edm, "is above tolerance", 10 * edmval);
338
340 }
341 }
342 // std::cout<<"result.back().Error().Dcovar()= "<<result.back().Error().Dcovar()<<std::endl;
343
344 print.Debug("Exiting successfully", "Ncalls", fcn.NumOfCalls(), "FCN", result.back().Fval(), "Edm", edm, "Requested",
345 edmval);
346
347 return FunctionMinimum(seed, result, fcn.Up());
348}
349
350} // namespace Minuit2
351
352} // namespace ROOT
#define g(i)
Definition RSha256.hxx:105
#define s0(x)
Definition RSha256.hxx:90
#define e(i)
Definition RSha256.hxx:103
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void gc
FunctionMinimum Minimum(const MnFcn &fMnFcn, const GradientCalculator &fGradienCalculator, const MinimumSeed &fMinimumSeed, const MnStrategy &fMnStrategy, unsigned int maxfcn, double edmval) const override
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) ...
interface class for gradient calculators
unsigned int size() const
Definition LAVector.h:231
void TraceIteration(int iter, const MinimumState &state) const
MinimumError keeps the inv.
const MnAlgebraicSymMatrix & InvHessian() const
const MnAlgebraicVector & Vec() const
const MinimumParameters & Parameters() const
Definition MinimumSeed.h:29
const MnMachinePrecision & Precision() const
Definition MinimumSeed.h:33
const MinimumState & State() const
Definition MinimumSeed.h:28
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
const MinimumError & Error() const
const MnAlgebraicVector & Vec() 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:40
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:147
void Info(const Ts &... args)
Definition MnPrint.h:141
void Warn(const Ts &... args)
Definition MnPrint.h:135
API class for defining four levels of strategies: low (0), medium (1), high (2), very high (>=3); act...
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...