ROOT  6.07/01
Reference Guide
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
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
const char * Name(unsigned int) const
void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
Definition: Ifit.C:26
static void PrintFcn(std::ostream &os, double value, bool endline=true)
Definition: MnPrint.cxx:51
static int Level()
Definition: MnPrint.cxx:47
#define assert(cond)
Definition: unittest.h:542
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
determines the relative floating point arithmetic precision.
const MnUserTransformation & Trafo() const
#define MN_INFO_MSG(str)
Definition: MnPrint.h:110
std::pair< FunctionGradient, MnAlgebraicVector > DeltaGradient(const MinimumParameters &, const FunctionGradient &) const
Double_t x[n]
Definition: legend1.C:17
class performing the numerical gradient calculation
API class for calculating the numerical covariance matrix (== 2x Inverse Hessian == 2x Inverse 2nd de...
Definition: MnHesse.h:40
virtual MinimumSeed operator()(const MnFcn &, const GradientCalculator &, const MnUserParameterState &, const MnStrategy &) 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
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)
double Eps2() const
eps2 returns 2*sqrt(eps)
#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...
double Estimate(const FunctionGradient &, const MinimumError &) const
const MnMachinePrecision & Precision() const
unsigned int ExtOfInt(unsigned int internal) const
const std::vector< double > & IntParameters() const
Class to calculate an initial estimate of the gradient.
unsigned int Strategy() const
Definition: MnStrategy.h:39
bool HasNegativeG2(const FunctionGradient &, const MnMachinePrecision &) const
In case that one of the components of the second derivative g2 calculated by the numerical Gradient c...
const MnAlgebraicVector & Grad() const
MinimumError keeps the inv.
Definition: MinimumError.h:26
unsigned int NumOfCalls() const
Definition: MnFcn.h:43
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
Definition: MinimumState.h:29
int printLevel
const MnUserCovariance & IntCovariance() const
API class for defining three levels of strategies: low (0), medium (1), high (>=2); acts on: Migrad (...
Definition: MnStrategy.h:27
const MnAlgebraicVector & G2() const
const Int_t n
Definition: legend1.C:16
interface class for gradient calculators