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
17#if defined(DEBUG) || defined(WARNINGMSG)
18#include "Minuit2/MnPrint.h"
19#endif
20
21
22namespace ROOT {
23
24 namespace Minuit2 {
25
26
27//#define DEBUG 1
29class MnStrategy;
30FunctionMinimum SimplexBuilder::Minimum(const MnFcn& mfcn, const GradientCalculator&, const MinimumSeed& seed, const MnStrategy&, unsigned int maxfcn, double minedm) const {
31 // find the minimum using the Simplex method of Nelder and Mead (does not use function gradient)
32 // method to find initial simplex is slightly different than in the orginal Fortran
33 // Minuit since has not been proofed that one to be better
34
35#ifdef DEBUG
36 std::cout << "Running Simplex with maxfcn = " << maxfcn << " minedm = " << minedm << std::endl;
37#endif
38
39 const MnMachinePrecision& prec = seed.Precision();
42
43 unsigned int n = x.size();
44 double wg = 1./double(n);
45 double alpha = 1., beta = 0.5, gamma = 2., rhomin = 4., rhomax = 8.;
46 double rho1 = 1. + alpha;
47 //double rho2 = rho1 + alpha*gamma;
48 //change proposed by david sachs (fnal)
49 double rho2 = 1. + alpha*gamma;
50
51
52 std::vector<std::pair<double, MnAlgebraicVector> > simpl; simpl.reserve(n+1);
53 simpl.push_back(std::pair<double, MnAlgebraicVector>(seed.Fval(), x));
54
55 unsigned int jl = 0, jh = 0;
56 double amin = seed.Fval(), aming = seed.Fval();
57
58 for(unsigned int i = 0; i < n; i++) {
59 double dmin = 8.*prec.Eps2()*(fabs(x(i)) + prec.Eps2());
60 if(step(i) < dmin) step(i) = dmin;
61 x(i) += step(i);
62 double tmp = mfcn(x);
63 if(tmp < amin) {
64 amin = tmp;
65 jl = i+1;
66 }
67 if(tmp > aming) {
68 aming = tmp;
69 jh = i+1;
70 }
71 simpl.push_back(std::pair<double, MnAlgebraicVector>(tmp, x));
72 x(i) -= step(i);
73 }
74 SimplexParameters simplex(simpl, jh, jl);
75
76#ifdef DEBUG
77 std::cout << "simplex initial parameters - min " << jl << " " << amin << " max " << jh << " " << aming << std::endl;
78 for (unsigned int i = 0; i < simplex.Simplex().size(); ++i)
79 std::cout << " i = " << i << " x = " << simplex(i).second << " fval(x) = " << simplex(i).first << std::endl;
80#endif
81
82 double edmPrev = simplex.Edm();
83 int niterations = 0;
84 do {
85 jl = simplex.Jl();
86 jh = simplex.Jh();
87 amin = simplex(jl).first;
88 edmPrev = simplex.Edm();
89
90#ifdef DEBUG
91 std::cout << "\n\nsimplex iteration: edm = " << simplex.Edm()
92 << "\n--> Min Param is " << jl << " pmin " << simplex(jl).second << " f(pmin) " << amin
93 << "\n--> Max param is " << jh << " " << simplex(jh).first << std::endl;
94
95 // std::cout << "ALL SIMPLEX PARAMETERS: "<< std::endl;
96 // for (unsigned int i = 0; i < simplex.Simplex().size(); ++i)
97 // std::cout << " i = " << i << " x = " << simplex(i).second << " fval(x) = " << simplex(i).first << std::endl;
98#endif
99
100 // trace the iterations (need to create a MinimunState although errors and gradient are not existing)
101 if (TraceIter() ) TraceIteration(niterations, MinimumState(MinimumParameters(simplex(jl).second,simplex(jl).first), simplex.Edm(), mfcn.NumOfCalls()) );
102 if (PrintLevel() > 1) MnPrint::PrintState(std::cout,simplex(jl).first, simplex.Edm(),mfcn.NumOfCalls(),"Simplex: Iteration # ", niterations);
103 niterations++;
104
105
106 MnAlgebraicVector pbar(n);
107 for(unsigned int i = 0; i < n+1; i++) {
108 if(i == jh) continue;
109 pbar += (wg*simplex(i).second);
110 }
111
112 MnAlgebraicVector pstar = (1. + alpha)*pbar - alpha*simplex(jh).second;
113 double ystar = mfcn(pstar);
114
115#ifdef DEBUG
116 std::cout << " pbar = " << pbar << std::endl;
117 std::cout << " pstar = " << pstar << " f(pstar) = " << ystar << std::endl;
118#endif
119
120 if(ystar > amin) {
121 if(ystar < simplex(jh).first) {
122 simplex.Update(ystar, pstar);
123 if(jh != simplex.Jh()) continue;
124 }
125 MnAlgebraicVector pstst = beta*simplex(jh).second + (1. - beta)*pbar;
126 double ystst = mfcn(pstst);
127#ifdef DEBUG
128 std::cout << "Reduced simplex pstst = " << pstst << " f(pstst) = " << ystst << std::endl;
129#endif
130 if(ystst > simplex(jh).first) break;
131 simplex.Update(ystst, pstst);
132 continue;
133 }
134
135 MnAlgebraicVector pstst = gamma*pstar + (1. - gamma)*pbar;
136 double ystst = mfcn(pstst);
137#ifdef DEBUG
138 std::cout << " pstst = " << pstst << " f(pstst) = " << ystst << std::endl;
139#endif
140
141 double y1 = (ystar - simplex(jh).first)*rho2;
142 double y2 = (ystst - simplex(jh).first)*rho1;
143 double rho = 0.5*(rho2*y1 - rho1*y2)/(y1 - y2);
144 if(rho < rhomin) {
145 if(ystst < simplex(jl).first) simplex.Update(ystst, pstst);
146 else simplex.Update(ystar, pstar);
147 continue;
148 }
149 if(rho > rhomax) rho = rhomax;
150 MnAlgebraicVector prho = rho*pbar + (1. - rho)*simplex(jh).second;
151 double yrho = mfcn(prho);
152#ifdef DEBUG
153 std::cout << " prho = " << prho << " f(prho) = " << yrho << std::endl;
154#endif
155 if(yrho < simplex(jl).first && yrho < ystst) {
156 simplex.Update(yrho, prho);
157 continue;
158 }
159 if(ystst < simplex(jl).first) {
160 simplex.Update(ystst, pstst);
161 continue;
162 }
163 if(yrho > simplex(jl).first) {
164 if(ystst < simplex(jl).first) simplex.Update(ystst, pstst);
165 else simplex.Update(ystar, pstar);
166 continue;
167 }
168 if(ystar > simplex(jh).first) {
169 pstst = beta*simplex(jh).second + (1. - beta)*pbar;
170 ystst = mfcn(pstst);
171 if(ystst > simplex(jh).first) break;
172 simplex.Update(ystst, pstst);
173 }
174#ifdef DEBUG
175 std::cout << "End loop : edm = " << simplex.Edm() << " pstst = " << pstst << " f(pstst) = " << ystst << std::endl;
176#endif
177 } while( (simplex.Edm() > minedm || edmPrev > minedm ) && mfcn.NumOfCalls() < maxfcn);
178
179 jl = simplex.Jl();
180 jh = simplex.Jh();
181 amin = simplex(jl).first;
182
183 MnAlgebraicVector pbar(n);
184 for(unsigned int i = 0; i < n+1; i++) {
185 if(i == jh) continue;
186 pbar += (wg*simplex(i).second);
187 }
188 double ybar = mfcn(pbar);
189 if(ybar < amin) simplex.Update(ybar, pbar);
190 else {
191 pbar = simplex(jl).second;
192 ybar = simplex(jl).first;
193 }
194
195 MnAlgebraicVector dirin = simplex.Dirin();
196 // Scale to sigmas on parameters werr^2 = dirin^2 * (up/edm)
197 dirin *= sqrt(mfcn.Up()/simplex.Edm());
198
199#ifdef DEBUG
200 std::cout << "End simplex " << simplex.Edm() << " pbar = " << pbar << " f(p) = " << ybar << std::endl;
201#endif
202
203
204 MinimumState st(MinimumParameters(pbar, dirin, ybar), simplex.Edm(), mfcn.NumOfCalls());
205
206 if (PrintLevel() > 1)
207 MnPrint::PrintState(std::cout,st,"Simplex: Final iteration");
208 if (TraceIter() ) TraceIteration(niterations, st);
209
210 if(mfcn.NumOfCalls() > maxfcn) {
211#ifdef WARNINGMSG
212 MN_INFO_MSG("Simplex did not converge, #fcn calls exhausted.");
213#endif
214 return FunctionMinimum(seed, std::vector<MinimumState>(1, st), mfcn.Up(), FunctionMinimum::MnReachedCallLimit());
215 }
216 if(simplex.Edm() > minedm) {
217#ifdef WARNINGMSG
218 MN_INFO_MSG("Simplex did not converge, edm > minedm.");
219#endif
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
#define MN_INFO_MSG(str)
Definition: MnPrint.h:110
double sqrt(double)
const MnAlgebraicVector & Gstep() const
class holding the full result of the minimization; both internal and external (MnUserParameterState) ...
void TraceIteration(int iter, const MinimumState &state) const
const MnAlgebraicVector & Vec() const
MinimumSeed contains the starting values for the minimization produced by the SeedGenerator.
Definition: MinimumSeed.h:31
Definition: MinimumSeed.h:49
const MinimumParameters & Parameters() const
Definition: MinimumSeed.h:47
const MnMachinePrecision & Precision() const
Definition: MinimumSeed.h:51
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
Definition: MinimumState.h:29
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
determines the relative floating point arithmetic precision.
double Eps2() const
eps2 returns 2*sqrt(eps)
static void PrintState(std::ostream &os, const MinimumState &state, const char *msg, int iter=-1)
Definition: MnPrint.cxx:58
API class for defining three levels of strategies: low (0), medium (1), high (>=2); acts on: Migrad (...
Definition: MnStrategy.h:27
virtual FunctionMinimum Minimum(const MnFcn &, const GradientCalculator &, const MinimumSeed &, const MnStrategy &, unsigned int, double) const
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< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
VSD Structures.
Definition: StringConv.hxx:21
static constexpr double second
Definition: first.py:1