Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
SimplexBuilder.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
12#include "Minuit2/MnFcn.h"
13#include "Minuit2/MinimumSeed.h"
16#include "Minuit2/MnPrint.h"
17
18namespace ROOT {
19
20namespace Minuit2 {
21
22class GradientCalculator;
23class MnStrategy;
25 const MnStrategy &, unsigned int maxfcn, double minedm) const
26{
27 // find the minimum using the Simplex method of Nelder and Mead (does not use function gradient)
28 // method to find initial simplex is slightly different than in the original Fortran
29 // Minuit since has not been proofed that one to be better
30
32
33 // synchronize print levels
34 MnPrint print("SimplexBuilder", PrintLevel());
35
36 print.Debug("Running with maxfcn", maxfcn, "minedm", minedm);
37
38 const MnMachinePrecision &prec = seed.Precision();
39 MnAlgebraicVector x = seed.Parameters().Vec();
40 MnAlgebraicVector step = 10. * seed.Gradient().Gstep();
41
42 const unsigned int n = x.size();
43
44 // If there are no parameters, we just return a minimum state
45 if (n == 0) {
47 mfcn.NumOfCalls());
48 return FunctionMinimum(seed, {&st, 1}, mfcn.Up());
49 }
50
51 const double wg = 1. / double(n);
52 const double alpha = 1.;
53 const double beta = 0.5;
54 const double gamma = 2.;
55 const double rhomin = 4.;
56 const double rhomax = 8.;
57 const double rho1 = 1. + alpha;
58 // double rho2 = rho1 + alpha*gamma;
59 // change proposed by david sachs (fnal)
60 const double rho2 = 1. + alpha * gamma;
61
62 std::vector<std::pair<double, MnAlgebraicVector>> simpl;
63 simpl.reserve(n + 1);
64 simpl.emplace_back(seed.Fval(), x);
65
66 unsigned int jl = 0;
67 unsigned int jh = 0;
68 double amin = seed.Fval();
69 double aming = seed.Fval();
70
71 for (unsigned int i = 0; i < n; i++) {
72 double dmin = 8. * prec.Eps2() * (std::abs(x(i)) + prec.Eps2());
73 if (step(i) < dmin)
74 step(i) = dmin;
75 x(i) += step(i);
76 double tmp = mfcnCaller(x);
77 if (tmp < amin) {
78 amin = tmp;
79 jl = i + 1;
80 }
81 if (tmp > aming) {
82 aming = tmp;
83 jh = i + 1;
84 }
85 simpl.emplace_back(tmp, x);
86 x(i) -= step(i);
87 }
88 SimplexParameters simplex(simpl, jh, jl);
89
90 print.Debug([&](std::ostream &os) {
91 os << "Initial parameters - min " << jl << " " << amin << " max " << jh << " " << aming << '\n';
92 for (unsigned int i = 0; i < simplex.Simplex().size(); ++i)
93 os << " i = " << i << " x = " << simplex(i).second << " fval(x) = " << simplex(i).first << '\n';
94 });
95
96 double edmPrev = simplex.Edm();
97 int niterations = 0;
98
99 do {
100 jl = simplex.Jl();
101 jh = simplex.Jh();
102 amin = simplex(jl).first;
103 edmPrev = simplex.Edm();
104
105 print.Debug("iteration: edm =", simplex.Edm(), '\n', "--> Min Param is", jl, "pmin", simplex(jl).second,
106 "f(pmin)", amin, '\n', "--> Max param is", jh, simplex(jh).first);
107
108 // std::cout << "ALL SIMPLEX PARAMETERS: "<< std::endl;
109 // for (unsigned int i = 0; i < simplex.Simplex().size(); ++i)
110 // std::cout << " i = " << i << " x = " << simplex(i).second << " fval(x) = " << simplex(i).first <<
111 // std::endl;
112
113 // trace the iterations (need to create a MinimumState although errors and gradient are not existing)
114 if (TraceIter())
115 TraceIteration(niterations, MinimumState(MinimumParameters(simplex(jl).second, simplex(jl).first),
116 simplex.Edm(), mfcn.NumOfCalls()));
117 print.Info(MnPrint::Oneline(simplex(jl).first, simplex.Edm(), mfcn.NumOfCalls(), niterations));
118 niterations++;
119
121 for (unsigned int i = 0; i < n + 1; i++) {
122 if (i == jh)
123 continue;
124 pbar += (wg * simplex(i).second);
125 }
126
127 MnAlgebraicVector pstar = (1. + alpha) * pbar - alpha * simplex(jh).second;
128 const double ystar = mfcnCaller(pstar);
129
130 print.Debug("pbar", pbar, "pstar", pstar, "f(pstar)", ystar);
131
132 if (ystar > amin) {
133 if (ystar < simplex(jh).first) {
134 simplex.Update(ystar, pstar);
135 if (jh != simplex.Jh())
136 continue;
137 }
138 MnAlgebraicVector pstst = beta * simplex(jh).second + (1. - beta) * pbar;
139 double ystst = mfcnCaller(pstst);
140
141 print.Debug("Reduced simplex pstst", pstst, "f(pstst)", ystst);
142
143 if (ystst > simplex(jh).first)
144 break;
145 simplex.Update(ystst, pstst);
146 continue;
147 }
148
149 MnAlgebraicVector pstst = gamma * pstar + (1. - gamma) * pbar;
150 double ystst = mfcnCaller(pstst);
151
152 print.Debug("pstst", pstst, "f(pstst)", ystst);
153
154 const double y1 = (ystar - simplex(jh).first) * rho2;
155 const double y2 = (ystst - simplex(jh).first) * rho1;
156 double rho = 0.5 * (rho2 * y1 - rho1 * y2) / (y1 - y2);
157 if (rho < rhomin) {
158 if (ystst < simplex(jl).first)
159 simplex.Update(ystst, pstst);
160 else
161 simplex.Update(ystar, pstar);
162 continue;
163 }
164 if (rho > rhomax)
165 rho = rhomax;
166 MnAlgebraicVector prho = rho * pbar + (1. - rho) * simplex(jh).second;
167 double yrho = mfcnCaller(prho);
168
169 print.Debug("prho", prho, "f(prho)", yrho);
170
171 if (yrho < simplex(jl).first && yrho < ystst) {
172 simplex.Update(yrho, prho);
173 continue;
174 }
175 if (ystst < simplex(jl).first) {
176 simplex.Update(ystst, pstst);
177 continue;
178 }
179 if (yrho > simplex(jl).first) {
180 if (ystst < simplex(jl).first)
181 simplex.Update(ystst, pstst);
182 else
183 simplex.Update(ystar, pstar);
184 continue;
185 }
186 if (ystar > simplex(jh).first) {
187 pstst = beta * simplex(jh).second + (1. - beta) * pbar;
189 if (ystst > simplex(jh).first)
190 break;
191 simplex.Update(ystst, pstst);
192 }
193
194 print.Debug("End loop : Edm", simplex.Edm(), "pstst", pstst, "f(pstst)", ystst);
195
196 } while ((simplex.Edm() > minedm || edmPrev > minedm) && mfcn.NumOfCalls() < maxfcn);
197
198 jl = simplex.Jl();
199 jh = simplex.Jh();
200 amin = simplex(jl).first;
201
203 for (unsigned int i = 0; i < n + 1; i++) {
204 if (i == jh)
205 continue;
206 pbar += (wg * simplex(i).second);
207 }
208 double ybar = mfcnCaller(pbar);
209 if (ybar < amin)
210 simplex.Update(ybar, pbar);
211 else {
212 pbar = simplex(jl).second;
213 ybar = simplex(jl).first;
214 }
215
216 MnAlgebraicVector dirin = simplex.Dirin();
217 // Scale to sigmas on parameters werr^2 = dirin^2 * (up/edm)
218 dirin *= std::sqrt(mfcn.Up() / simplex.Edm());
219
220 print.Debug("End simplex edm =", simplex.Edm(), "pbar =", pbar, "f(p) =", ybar);
221
222 MinimumState st(MinimumParameters(pbar, dirin, ybar), simplex.Edm(), mfcn.NumOfCalls());
223
224 print.Info("Final iteration", MnPrint::Oneline(st));
225
226 if (TraceIter())
228
229 if (mfcn.NumOfCalls() > maxfcn) {
230 print.Warn("Simplex did not converge, #fcn calls exhausted");
231
233 }
234 if (simplex.Edm() > minedm) {
235 print.Warn("Simplex did not converge, edm > minedm");
236
237 return FunctionMinimum(seed, {&st, 1}, mfcn.Up(), FunctionMinimum::MnAboveMaxEdm);
238 }
239
240 return FunctionMinimum(seed, {&st, 1}, mfcn.Up());
241}
242
243} // namespace Minuit2
244
245} // namespace ROOT
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
Option_t Option_t TPoint TPoint const char y2
Option_t Option_t TPoint TPoint const char y1
class holding the full result of the minimization; both internal and external (MnUserParameterState) ...
interface class for gradient calculators
void TraceIteration(int iter, const MinimumState &state) const
const FunctionGradient & Gradient() const
Definition MinimumSeed.h:32
const MinimumParameters & Parameters() const
Definition MinimumSeed.h:30
const MnMachinePrecision & Precision() const
Definition MinimumSeed.h:34
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
Wrapper class to FCNBase interface used internally by Minuit.
Definition MnFcn.h:34
Sets the relative floating point (double) arithmetic precision.
void Debug(const Ts &... args)
Definition MnPrint.h:135
void Info(const Ts &... args)
Definition MnPrint.h:129
void Warn(const Ts &... args)
Definition MnPrint.h:123
API class for defining four levels of strategies: low (0), medium (1), high (2), very high (>=3); act...
Definition MnStrategy.h:27
FunctionMinimum Minimum(const MnFcn &, const GradientCalculator &, const MinimumSeed &, const MnStrategy &, unsigned int, double) const override
class describing the simplex set of points (f(x), x ) which evolve during the minimization iteration ...
MnAlgebraicVector Dirin() const
void Update(double, const MnAlgebraicVector &)
const std::vector< std::pair< double, MnAlgebraicVector > > & Simplex() const
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...