Logo ROOT  
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"
19#include "Minuit2/MnMatrix.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
45namespace 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);
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
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
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
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_INFO_MSG(str)
Definition: MnPrint.h:110
#define MN_INFO_VAL(x)
Definition: MnPrint.h:116
#define MN_ERROR_MSG(str)
Definition: MnPrint.h:113
const MnAlgebraicVector & Gstep() const
const MnAlgebraicVector & Grad() const
const MnAlgebraicVector & G2() const
interface class for gradient calculators
HessianGradientCalculator: class to calculate Gradient for Hessian.
std::pair< FunctionGradient, MnAlgebraicVector > DeltaGradient(const MinimumParameters &, const FunctionGradient &) const
Class to calculate an initial estimate of the gradient.
Class describing a symmetric matrix of size n.
Definition: LASymMatrix.h:51
MinimumError keeps the inv.
Definition: MinimumError.h:26
MinimumSeed contains the starting values for the minimization produced by the SeedGenerator.
Definition: MinimumSeed.h:31
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
unsigned int NumOfCalls() const
Definition: MnFcn.h:43
API class for calculating the numerical covariance matrix (== 2x Inverse Hessian == 2x Inverse 2nd de...
Definition: MnHesse.h:40
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
static void PrintFcn(std::ostream &os, double value, bool endline=true)
Definition: MnPrint.cxx:51
static int Level()
Definition: MnPrint.cxx:47
virtual MinimumSeed operator()(const MnFcn &, const GradientCalculator &, const MnUserParameterState &, const MnStrategy &) const
API class for defining three levels of strategies: low (0), medium (1), high (>=2); acts on: Migrad (...
Definition: MnStrategy.h:27
unsigned int Strategy() const
Definition: MnStrategy.h:39
class which holds the external user and/or internal Minuit representation of the parameters and error...
const MnMachinePrecision & Precision() const
const std::vector< double > & IntParameters() const
const MnUserTransformation & Trafo() const
const MnUserCovariance & IntCovariance() const
unsigned int ExtOfInt(unsigned int internal) const
const char * Name(unsigned int) const
In case that one of the components of the second derivative g2 calculated by the numerical Gradient c...
bool HasNegativeG2(const FunctionGradient &, const MnMachinePrecision &) const
class performing the numerical gradient calculation
double Estimate(const FunctionGradient &, const MinimumError &) const
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
VSD Structures.
Definition: StringConv.hxx:21