Logo ROOT   6.10/09
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 
10 #include "Minuit2/SimplexBuilder.h"
12 #include "Minuit2/MnFcn.h"
13 #include "Minuit2/MinimumSeed.h"
15 #include "Minuit2/MinimumState.h"
16 
17 #if defined(DEBUG) || defined(WARNINGMSG)
18 #include "Minuit2/MnPrint.h"
19 #endif
20 
21 
22 namespace ROOT {
23 
24  namespace Minuit2 {
25 
26 
27 //#define DEBUG 1
28 class GradientCalculator;
29 class MnStrategy;
30 FunctionMinimum 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();
40  MnAlgebraicVector x = seed.Parameters().Vec();
41  MnAlgebraicVector step = 10.*seed.Gradient().Gstep();
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
const FunctionGradient & Gradient() const
Definition: MinimumSeed.h:49
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
const MnAlgebraicVector & Vec() const
MinimumSeed contains the starting values for the minimization produced by the SeedGenerator.
Definition: MinimumSeed.h:31
double Up() const
Definition: MnFcn.cxx:35
const MnMachinePrecision & Precision() const
Definition: MinimumSeed.h:51
determines the relative floating point arithmetic precision.
double beta(double x, double y)
Calculates the beta function.
unsigned int NumOfCalls() const
Definition: MnFcn.h:43
#define MN_INFO_MSG(str)
Definition: MnPrint.h:110
double sqrt(double)
Double_t x[n]
Definition: legend1.C:17
const MinimumParameters & Parameters() const
Definition: MinimumSeed.h:47
class holding the full result of the minimization; both internal and external (MnUserParameterState) ...
static void PrintState(std::ostream &os, const MinimumState &state, const char *msg, int iter=-1)
Definition: MnPrint.cxx:58
Wrapper class to FCNBase interface used internally by Minuit.
Definition: MnFcn.h:33
double gamma(double x)
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
void Update(double, const MnAlgebraicVector &)
unsigned int size() const
Definition: LAVector.h:198
const MnAlgebraicVector & Gstep() const
class describing the simplex set of points (f(x), x ) which evolve during the minimization iteration ...
const std::vector< std::pair< double, MnAlgebraicVector > > & Simplex() const
virtual FunctionMinimum Minimum(const MnFcn &, const GradientCalculator &, const MinimumSeed &, const MnStrategy &, unsigned int, double) const
double Eps2() const
eps2 returns 2*sqrt(eps)
void TraceIteration(int iter, const MinimumState &state) const
MnAlgebraicVector Dirin() const
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
Definition: MinimumState.h:29
Definition: first.py:1
API class for defining three levels of strategies: low (0), medium (1), high (>=2); acts on: Migrad (...
Definition: MnStrategy.h:27
const Int_t n
Definition: legend1.C:16
interface class for gradient calculators