ROOT   Reference Guide
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 orginal Fortran
29 // Minuit since has not been proofed that one to be better
30
31 // synchronize print levels
32 MnPrint print("SimplexBuilder", PrintLevel());
33
34 print.Debug("Running with maxfcn", maxfcn, "minedm", minedm);
35
36 const MnMachinePrecision &prec = seed.Precision();
38 MnAlgebraicVector step = 10. * seed.Gradient().Gstep();
39
40 unsigned int n = x.size();
41 double wg = 1. / double(n);
42 double alpha = 1., beta = 0.5, gamma = 2., rhomin = 4., rhomax = 8.;
43 double rho1 = 1. + alpha;
44 // double rho2 = rho1 + alpha*gamma;
45 // change proposed by david sachs (fnal)
46 double rho2 = 1. + alpha * gamma;
47
48 std::vector<std::pair<double, MnAlgebraicVector>> simpl;
49 simpl.reserve(n + 1);
50 simpl.push_back(std::pair<double, MnAlgebraicVector>(seed.Fval(), x));
51
52 unsigned int jl = 0, jh = 0;
53 double amin = seed.Fval(), aming = seed.Fval();
54
55 for (unsigned int i = 0; i < n; i++) {
56 double dmin = 8. * prec.Eps2() * (std::fabs(x(i)) + prec.Eps2());
57 if (step(i) < dmin)
58 step(i) = dmin;
59 x(i) += step(i);
60 double tmp = mfcn(x);
61 if (tmp < amin) {
62 amin = tmp;
63 jl = i + 1;
64 }
65 if (tmp > aming) {
66 aming = tmp;
67 jh = i + 1;
68 }
69 simpl.push_back(std::pair<double, MnAlgebraicVector>(tmp, x));
70 x(i) -= step(i);
71 }
72 SimplexParameters simplex(simpl, jh, jl);
73
74 print.Debug([&](std::ostream &os) {
75 os << "Initial parameters - min " << jl << " " << amin << " max " << jh << " " << aming << '\n';
76 for (unsigned int i = 0; i < simplex.Simplex().size(); ++i)
77 os << " i = " << i << " x = " << simplex(i).second << " fval(x) = " << simplex(i).first << '\n';
78 });
79
80 double edmPrev = simplex.Edm();
81 int niterations = 0;
82 do {
83 jl = simplex.Jl();
84 jh = simplex.Jh();
85 amin = simplex(jl).first;
86 edmPrev = simplex.Edm();
87
88 print.Debug("iteration: edm =", simplex.Edm(), '\n', "--> Min Param is", jl, "pmin", simplex(jl).second,
89 "f(pmin)", amin, '\n', "--> Max param is", jh, simplex(jh).first);
90
91 // std::cout << "ALL SIMPLEX PARAMETERS: "<< std::endl;
92 // for (unsigned int i = 0; i < simplex.Simplex().size(); ++i)
93 // std::cout << " i = " << i << " x = " << simplex(i).second << " fval(x) = " << simplex(i).first <<
94 // std::endl;
95
96 // trace the iterations (need to create a MinimunState although errors and gradient are not existing)
97 if (TraceIter())
98 TraceIteration(niterations, MinimumState(MinimumParameters(simplex(jl).second, simplex(jl).first),
99 simplex.Edm(), mfcn.NumOfCalls()));
100 print.Info(MnPrint::Oneline(simplex(jl).first, simplex.Edm(), mfcn.NumOfCalls(), niterations));
101 niterations++;
102
103 MnAlgebraicVector pbar(n);
104 for (unsigned int i = 0; i < n + 1; i++) {
105 if (i == jh)
106 continue;
107 pbar += (wg * simplex(i).second);
108 }
109
110 MnAlgebraicVector pstar = (1. + alpha) * pbar - alpha * simplex(jh).second;
111 double ystar = mfcn(pstar);
112
113 print.Debug("pbar", pbar, "pstar", pstar, "f(pstar)", ystar);
114
115 if (ystar > amin) {
116 if (ystar < simplex(jh).first) {
117 simplex.Update(ystar, pstar);
118 if (jh != simplex.Jh())
119 continue;
120 }
121 MnAlgebraicVector pstst = beta * simplex(jh).second + (1. - beta) * pbar;
122 double ystst = mfcn(pstst);
123
124 print.Debug("Reduced simplex pstst", pstst, "f(pstst)", ystst);
125
126 if (ystst > simplex(jh).first)
127 break;
128 simplex.Update(ystst, pstst);
129 continue;
130 }
131
132 MnAlgebraicVector pstst = gamma * pstar + (1. - gamma) * pbar;
133 double ystst = mfcn(pstst);
134
135 print.Debug("pstst", pstst, "f(pstst)", ystst);
136
137 double y1 = (ystar - simplex(jh).first) * rho2;
138 double y2 = (ystst - simplex(jh).first) * rho1;
139 double rho = 0.5 * (rho2 * y1 - rho1 * y2) / (y1 - y2);
140 if (rho < rhomin) {
141 if (ystst < simplex(jl).first)
142 simplex.Update(ystst, pstst);
143 else
144 simplex.Update(ystar, pstar);
145 continue;
146 }
147 if (rho > rhomax)
148 rho = rhomax;
149 MnAlgebraicVector prho = rho * pbar + (1. - rho) * simplex(jh).second;
150 double yrho = mfcn(prho);
151
152 print.Debug("prho", prho, "f(prho)", yrho);
153
154 if (yrho < simplex(jl).first && yrho < ystst) {
155 simplex.Update(yrho, prho);
156 continue;
157 }
158 if (ystst < simplex(jl).first) {
159 simplex.Update(ystst, pstst);
160 continue;
161 }
162 if (yrho > simplex(jl).first) {
163 if (ystst < simplex(jl).first)
164 simplex.Update(ystst, pstst);
165 else
166 simplex.Update(ystar, pstar);
167 continue;
168 }
169 if (ystar > simplex(jh).first) {
170 pstst = beta * simplex(jh).second + (1. - beta) * pbar;
171 ystst = mfcn(pstst);
172 if (ystst > simplex(jh).first)
173 break;
174 simplex.Update(ystst, pstst);
175 }
176
177 print.Debug("End loop : Edm", simplex.Edm(), "pstst", pstst, "f(pstst)", ystst);
178
179 } while ((simplex.Edm() > minedm || edmPrev > minedm) && mfcn.NumOfCalls() < maxfcn);
180
181 jl = simplex.Jl();
182 jh = simplex.Jh();
183 amin = simplex(jl).first;
184
185 MnAlgebraicVector pbar(n);
186 for (unsigned int i = 0; i < n + 1; i++) {
187 if (i == jh)
188 continue;
189 pbar += (wg * simplex(i).second);
190 }
191 double ybar = mfcn(pbar);
192 if (ybar < amin)
193 simplex.Update(ybar, pbar);
194 else {
195 pbar = simplex(jl).second;
196 ybar = simplex(jl).first;
197 }
198
199 MnAlgebraicVector dirin = simplex.Dirin();
200 // Scale to sigmas on parameters werr^2 = dirin^2 * (up/edm)
201 dirin *= std::sqrt(mfcn.Up() / simplex.Edm());
202
203 print.Debug("End simplex edm =", simplex.Edm(), "pbar =", pbar, "f(p) =", ybar);
204
205 MinimumState st(MinimumParameters(pbar, dirin, ybar), simplex.Edm(), mfcn.NumOfCalls());
206
207 print.Info("Final iteration", MnPrint::Oneline(st));
208
209 if (TraceIter())
210 TraceIteration(niterations, st);
211
212 if (mfcn.NumOfCalls() > maxfcn) {
213 print.Warn("Simplex did not converge, #fcn calls exhausted");
214
215 return FunctionMinimum(seed, std::vector<MinimumState>(1, st), mfcn.Up(), FunctionMinimum::MnReachedCallLimit);
216 }
217 if (simplex.Edm() > minedm) {
218 print.Warn("Simplex did not converge, edm > minedm");
219
220 return FunctionMinimum(seed, std::vector<MinimumState>(1, st), mfcn.Up(), FunctionMinimum::MnAboveMaxEdm);
221 }
222
223 return FunctionMinimum(seed, std::vector<MinimumState>(1, st), mfcn.Up());
224}
225
226} // namespace Minuit2
227
228} // namespace ROOT
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
Option_t Option_t TPoint TPoint const char y2
Option_t Option_t TPoint TPoint const char y1
const MnAlgebraicVector & Gstep() const
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 MnAlgebraicVector & Vec() const
const FunctionGradient & Gradient() const
Definition: MinimumSeed.h:31
const MinimumParameters & Parameters() const
Definition: MinimumSeed.h:29
const MnMachinePrecision & Precision() const
Definition: MinimumSeed.h:33
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
Definition: MinimumState.h:27
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
Sets the relative floating point (double) arithmetic precision.
double Eps2() const
eps2 returns 2*sqrt(eps)
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
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 beta(double x, double y)
Calculates the beta function.
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
double gamma(double x)
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.
static constexpr double second
Definition: first.py:1