Logo ROOT   6.10/09
Reference Guide
MnSeedGenerator.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 
11 #include "Minuit2/MinimumSeed.h"
12 #include "Minuit2/MnFcn.h"
18 #include "Minuit2/MinimumError.h"
19 #include "Minuit2/MnMatrix.h"
22 #include "Minuit2/MnLineSearch.h"
24 #include "Minuit2/MinimumState.h"
26 #include "Minuit2/MnStrategy.h"
27 #include "Minuit2/MnHesse.h"
33 
34 //#define DEBUG
35 
36 //#if defined(DEBUG) || defined(WARNINGMSG)
37 #include "Minuit2/MnPrint.h"
38 //#endif
39 
40 
41 
42 #include <math.h>
43 
44 
45 namespace ROOT {
46 
47  namespace Minuit2 {
48 
49 
51 
52 
53  // find seed (initial minimization point) using the calculated gradient
54  unsigned int n = st.VariableParameters();
55  const MnMachinePrecision& prec = st.Precision();
56 
57 #ifdef DEBUG
58  std::cout << "MnSeedGenerator: operator() - var par = " << n << " mnfcn pointer " << &fcn << std::endl;
59 #endif
60 
61  int printLevel = MnPrint::Level();
62 
63  // initial starting values
65  for(unsigned int i = 0; i < n; i++) x(i) = st.IntParameters()[i];
66  double fcnmin = fcn(x);
67 
68  if (printLevel > 1) {
69  std::cout << "MnSeedGenerator: for initial parameters FCN = ";
70  MnPrint::PrintFcn(std::cout,fcnmin);
71  }
72 
73  MinimumParameters pa(x, fcnmin);
74  FunctionGradient dgrad = gc(pa);
75  MnAlgebraicSymMatrix mat(n);
76  double dcovar = 1.;
77  if(st.HasCovariance()) {
78  for(unsigned int i = 0; i < n; i++)
79  for(unsigned int j = i; j < n; j++) mat(i,j) = st.IntCovariance()(i,j);
80  dcovar = 0.;
81  } else {
82  for(unsigned int i = 0; i < n; i++)
83  mat(i,i) = (fabs(dgrad.G2()(i)) > prec.Eps2() ? 1./dgrad.G2()(i) : 1.);
84  }
85  MinimumError err(mat, dcovar);
86 
87  double edm = VariableMetricEDMEstimator().Estimate(dgrad, err);
88  MinimumState state(pa, err, dgrad, edm, fcn.NumOfCalls());
89 
90  if (printLevel >1) {
91  MnPrint::PrintState(std::cout, state, "MnSeedGenerator: Initial state: ");
92  }
93 
95  if(ng2ls.HasNegativeG2(dgrad, prec)) {
96 #ifdef DEBUG
97  std::cout << "MnSeedGenerator: Negative G2 Found: " << std::endl;
98  std::cout << x << std::endl;
99  std::cout << dgrad.Grad() << std::endl;
100  std::cout << dgrad.G2() << std::endl;
101 #endif
102  state = ng2ls(fcn, state, gc, prec);
103 
104  if (printLevel >1) {
105  MnPrint::PrintState(std::cout, state, "MnSeedGenerator: Negative G2 found - new state: ");
106  }
107 
108  }
109 
110 
111  if(stra.Strategy() == 2 && !st.HasCovariance()) {
112  //calculate full 2nd derivative
113 #ifdef DEBUG
114  std::cout << "MnSeedGenerator: calling MnHesse " << std::endl;
115 #endif
116  MinimumState tmp = MnHesse(stra)(fcn, state, st.Trafo());
117 
118  if (printLevel >1) {
119  MnPrint::PrintState(std::cout, tmp, "MnSeedGenerator: run Hesse - new state: ");
120  }
121 
122  return MinimumSeed(tmp, st.Trafo());
123  }
124 
125  return MinimumSeed(state, st.Trafo());
126 }
127 
128 
130  // find seed (initial point for minimization) using analytical gradient
131  unsigned int n = st.VariableParameters();
132  const MnMachinePrecision& prec = st.Precision();
133 
134  // initial starting values
135  MnAlgebraicVector x(n);
136  for(unsigned int i = 0; i < n; i++) x(i) = st.IntParameters()[i];
137  double fcnmin = fcn(x);
138  MinimumParameters pa(x, fcnmin);
139 
140  InitialGradientCalculator igc(fcn, st.Trafo(), stra);
141  FunctionGradient tmp = igc(pa);
142  FunctionGradient grd = gc(pa);
143  FunctionGradient dgrad(grd.Grad(), tmp.G2(), tmp.Gstep());
144 
145  if(gc.CheckGradient()) {
146  bool good = true;
147  HessianGradientCalculator hgc(fcn, st.Trafo(), MnStrategy(2));
148  std::pair<FunctionGradient, MnAlgebraicVector> hgrd = hgc.DeltaGradient(pa, dgrad);
149  for(unsigned int i = 0; i < n; i++) {
150  if(fabs(hgrd.first.Grad()(i) - grd.Grad()(i)) > hgrd.second(i)) {
151 #ifdef WARNINGMSG
152  MN_INFO_MSG("MnSeedGenerator:gradient discrepancy of external Parameter too large");
153  int externalParameterIndex = st.Trafo().ExtOfInt(i);
154  const char * parameter_name = st.Trafo().Name(externalParameterIndex);
155  MN_INFO_VAL(parameter_name);
156  MN_INFO_VAL(externalParameterIndex);
157  MN_INFO_VAL2("internal",i);
158 #endif
159  good = false;
160  }
161  }
162  if(!good) {
163 #ifdef WARNINGMSG
164  MN_ERROR_MSG("Minuit does not accept user specified Gradient. To force acceptance, override 'virtual bool CheckGradient() const' of FCNGradientBase.h in the derived class.");
165 #endif
166  assert(good);
167  }
168  }
169 
170  MnAlgebraicSymMatrix mat(n);
171  double dcovar = 1.;
172  if(st.HasCovariance()) {
173  for(unsigned int i = 0; i < n; i++)
174  for(unsigned int j = i; j < n; j++) mat(i,j) = st.IntCovariance()(i,j);
175  dcovar = 0.;
176  } else {
177  for(unsigned int i = 0; i < n; i++)
178  mat(i,i) = (fabs(dgrad.G2()(i)) > prec.Eps2() ? 1./dgrad.G2()(i) : 1.);
179  }
180  MinimumError err(mat, dcovar);
181  double edm = VariableMetricEDMEstimator().Estimate(dgrad, err);
182  MinimumState state(pa, err, dgrad, edm, fcn.NumOfCalls());
183 
184  NegativeG2LineSearch ng2ls;
185  if(ng2ls.HasNegativeG2(dgrad, prec)) {
186  Numerical2PGradientCalculator ngc(fcn, st.Trafo(), stra);
187  state = ng2ls(fcn, state, ngc, prec);
188  }
189 
190  if(stra.Strategy() == 2 && !st.HasCovariance()) {
191  //calculate full 2nd derivative
192  MinimumState tmpState = MnHesse(stra)(fcn, state, st.Trafo());
193  return MinimumSeed(tmpState, st.Trafo());
194  }
195 
196  return MinimumSeed(state, st.Trafo());
197 }
198 
199  } // namespace Minuit2
200 
201 } // namespace ROOT
#define MN_INFO_VAL2(loc, x)
Definition: MnPrint.h:130
#define MN_ERROR_MSG(str)
Definition: MnPrint.h:113
double Estimate(const FunctionGradient &, const MinimumError &) const
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
static void PrintFcn(std::ostream &os, double value, bool endline=true)
Definition: MnPrint.cxx:51
static int Level()
Definition: MnPrint.cxx:47
bool HasNegativeG2(const FunctionGradient &, const MnMachinePrecision &) const
MinimumSeed contains the starting values for the minimization produced by the SeedGenerator.
Definition: MinimumSeed.h:31
Class describing a symmetric matrix of size n.
Definition: LASymMatrix.h:51
unsigned int Strategy() const
Definition: MnStrategy.h:39
unsigned int ExtOfInt(unsigned int internal) const
const std::vector< double > & IntParameters() const
determines the relative floating point arithmetic precision.
unsigned int NumOfCalls() const
Definition: MnFcn.h:43
#define MN_INFO_MSG(str)
Definition: MnPrint.h:110
Double_t x[n]
Definition: legend1.C:17
class performing the numerical gradient calculation
const MnMachinePrecision & Precision() const
API class for calculating the numerical covariance matrix (== 2x Inverse Hessian == 2x Inverse 2nd de...
Definition: MnHesse.h:40
const MnAlgebraicVector & G2() const
const char * Name(unsigned int) const
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
const MnUserTransformation & Trafo() const
HessianGradientCalculator: class to calculate Gradient for Hessian.
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
#define MN_INFO_VAL(x)
Definition: MnPrint.h:116
class which holds the external user and/or internal Minuit representation of the parameters and error...
Class to calculate an initial estimate of the gradient.
In case that one of the components of the second derivative g2 calculated by the numerical Gradient c...
virtual MinimumSeed operator()(const MnFcn &, const GradientCalculator &, const MnUserParameterState &, const MnStrategy &) const
const MnAlgebraicVector & Grad() const
double Eps2() const
eps2 returns 2*sqrt(eps)
const MnUserCovariance & IntCovariance() const
MinimumError keeps the inv.
Definition: MinimumError.h:26
std::pair< FunctionGradient, MnAlgebraicVector > DeltaGradient(const MinimumParameters &, const FunctionGradient &) const
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
Definition: MinimumState.h:29
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