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
10 #include "Minuit2/SimplexBuilder.h"
12 #include "Minuit2/MnFcn.h"
13 #include "Minuit2/MinimumSeed.h"
15 #include "Minuit2/MinimumState.h"
16 #include "Minuit2/MnPrint.h"
17
18 namespace ROOT {
19
20 namespace Minuit2 {
21
23 class 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();
37  MnAlgebraicVector x = seed.Parameters().Vec();
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
ROOT::Math::Cephes::gamma
double gamma(double x)
Definition: SpecFuncCephes.cxx:339
ROOT::Minuit2::MnMachinePrecision::Eps2
double Eps2() const
eps2 returns 2*sqrt(eps)
Definition: MnMachinePrecision.h:41
ROOT::Minuit2::MinimumBuilder::PrintLevel
int PrintLevel() const
Definition: MinimumBuilder.h:38
n
const Int_t n
Definition: legend1.C:16
first
Definition: first.py:1
ROOT::Minuit2::MnFcn
Wrapper class to FCNBase interface used internally by Minuit.
Definition: MnFcn.h:30
ROOT::Minuit2::MnFcn::NumOfCalls
unsigned int NumOfCalls() const
Definition: MnFcn.h:39
ROOT::Minuit2::MinimumState
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
Definition: MinimumState.h:27
ROOT::Minuit2::LAVector
Definition: LAVector.h:32
ROOT::Minuit2::MnPrint::Debug
void Debug(const Ts &... args)
Definition: MnPrint.h:138
FunctionMinimum.h
MinimumSeed.h
ROOT::Minuit2::MinimumParameters::Vec
const MnAlgebraicVector & Vec() const
Definition: MinimumParameters.h:39
SimplexBuilder.h
ROOT::Minuit2::SimplexParameters::Update
void Update(double, const MnAlgebraicVector &)
Definition: SimplexParameters.cxx:16
ROOT::Minuit2::MnPrint::Oneline
Definition: MnPrint.h:85
x
Double_t x[n]
Definition: legend1.C:17
ROOT::Math::beta
double beta(double x, double y)
Calculates the beta function.
Definition: SpecFuncMathCore.cxx:111
Definition: MinimumSeed.h:31
MnFcn.h
ROOT::Minuit2::MnPrint::Warn
void Warn(const Ts &... args)
Definition: MnPrint.h:126
ROOT::Math::fabs
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
Definition: UnaryOperators.h:131
ROOT::Minuit2::MnFcn::Up
double Up() const
Definition: MnFcn.cxx:39
size
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
ROOT::Minuit2::FunctionMinimum::MnAboveMaxEdm
@ MnAboveMaxEdm
Definition: FunctionMinimum.h:42
ROOT::Minuit2::MinimumBuilder::TraceIteration
void TraceIteration(int iter, const MinimumState &state) const
Definition: MinimumBuilder.h:49
ROOT::Minuit2::SimplexParameters::Dirin
MnAlgebraicVector Dirin() const
Definition: SimplexParameters.cxx:33
ROOT::Minuit2::MnMachinePrecision
Sets the relative floating point (double) arithmetic precision.
Definition: MnMachinePrecision.h:32
ROOT::Minuit2::MinimumSeed::Precision
const MnMachinePrecision & Precision() const
Definition: MinimumSeed.h:33
double
double
Definition: Converters.cxx:939
sqrt
double sqrt(double)
const MnAlgebraicVector & Gstep() const
ROOT::Minuit2::SimplexParameters::Simplex
const std::vector< std::pair< double, MnAlgebraicVector > > & Simplex() const
Definition: SimplexParameters.h:41
ROOT::Minuit2::FunctionMinimum
class holding the full result of the minimization; both internal and external (MnUserParameterState) ...
Definition: FunctionMinimum.h:36
ROOT::Minuit2::MnStrategy
API class for defining three levels of strategies: low (0), medium (1), high (>=2); acts on: Migrad (...
Definition: MnStrategy.h:27
ROOT::Minuit2::FunctionMinimum::MnReachedCallLimit
@ MnReachedCallLimit
Definition: FunctionMinimum.h:41
ROOT::Minuit2::SimplexParameters::Jh
unsigned int Jh() const
Definition: SimplexParameters.h:49
ROOT::Minuit2::MinimumSeed::Parameters
const MinimumParameters & Parameters() const
Definition: MinimumSeed.h:29
ROOT::Minuit2::MinimumBuilder::TraceIter
bool TraceIter() const
Definition: MinimumBuilder.h:40
ROOT::Minuit2::SimplexParameters::Edm
double Edm() const
Definition: SimplexParameters.h:51
ROOT::Minuit2::SimplexParameters::Jl
unsigned int Jl() const
Definition: SimplexParameters.h:50
MnPrint.h
TGeant4Unit::second
static constexpr double second
Definition: TGeant4SystemOfUnits.h:151
ROOT::Minuit2::SimplexParameters
class describing the simplex set of points (f(x), x ) which evolve during the minimization iteration ...
Definition: SimplexParameters.h:29
ROOT::Minuit2::SimplexBuilder::Minimum
virtual FunctionMinimum Minimum(const MnFcn &, const GradientCalculator &, const MinimumSeed &, const MnStrategy &, unsigned int, double) const
Definition: SimplexBuilder.cxx:24
ROOT::Minuit2::MinimumSeed::Fval
double Fval() const
Definition: MinimumSeed.h:34
ROOT
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
Definition: EExecutionPolicy.hxx:4
SimplexParameters.h
MinimumState.h
ROOT::Minuit2::MinimumSeed
Definition: MinimumSeed.h:23
ROOT::Minuit2::MnPrint::Info
void Info(const Ts &... args)
Definition: MnPrint.h:132
ROOT::Minuit2::MnPrint
Definition: MnPrint.h:73