Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
MnHesse.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/MnHesse.h"
11 #include "Minuit2/MnUserParameterState.h"
12 #include "Minuit2/MnUserFcn.h"
13 #include "Minuit2/FCNBase.h"
14 #include "Minuit2/MnPosDef.h"
15 #include "Minuit2/HessianGradientCalculator.h"
16 #include "Minuit2/Numerical2PGradientCalculator.h"
17 #include "Minuit2/InitialGradientCalculator.h"
18 #include "Minuit2/MinimumState.h"
19 #include "Minuit2/VariableMetricEDMEstimator.h"
20 #include "Minuit2/FunctionMinimum.h"
21 
22 //#define DEBUG
23 
24 #if defined(DEBUG) || defined(WARNINGMSG)
25 #include "Minuit2/MnPrint.h"
26 #endif
27 #if defined(DEBUG) && !defined(WARNINGMSG)
28 #define WARNINGMSG
29 #endif
30 
31 #include "Minuit2/MPIProcess.h"
32 
33 namespace ROOT {
34 
35  namespace Minuit2 {
36 
37 
38 MnUserParameterState MnHesse::operator()(const FCNBase& fcn, const std::vector<double>& par, const std::vector<double>& err, unsigned int maxcalls) const {
39  // interface from vector of params and errors
40  return (*this)(fcn, MnUserParameterState(par, err), maxcalls);
41 }
42 
43 MnUserParameterState MnHesse::operator()(const FCNBase& fcn, const std::vector<double>& par, unsigned int nrow, const std::vector<double>& cov, unsigned int maxcalls) const {
44  // interface from vector of params and covariance
45  return (*this)(fcn, MnUserParameterState(par, cov, nrow), maxcalls);
46 }
47 
48 MnUserParameterState MnHesse::operator()(const FCNBase& fcn, const std::vector<double>& par, const MnUserCovariance& cov, unsigned int maxcalls) const {
49  // interface from vector of params and covariance
50  return (*this)(fcn, MnUserParameterState(par, cov), maxcalls);
51 }
52 
53 MnUserParameterState MnHesse::operator()(const FCNBase& fcn, const MnUserParameters& par, unsigned int maxcalls) const {
54  // interface from MnUserParameters
55  return (*this)(fcn, MnUserParameterState(par), maxcalls);
56 }
57 
58 MnUserParameterState MnHesse::operator()(const FCNBase& fcn, const MnUserParameters& par, const MnUserCovariance& cov, unsigned int maxcalls) const {
59  // interface from MnUserParameters and MnUserCovariance
60  return (*this)(fcn, MnUserParameterState(par, cov), maxcalls);
61 }
62 
63 MnUserParameterState MnHesse::operator()(const FCNBase& fcn, const MnUserParameterState& state, unsigned int maxcalls) const {
64  // interface from MnUserParameterState
65  // create a new Minimum state and use that interface
66  unsigned int n = state.VariableParameters();
67  MnUserFcn mfcn(fcn, state.Trafo(),state.NFcn());
68  MnAlgebraicVector x(n);
69  for(unsigned int i = 0; i < n; i++) x(i) = state.IntParameters()[i];
70  double amin = mfcn(x);
71  Numerical2PGradientCalculator gc(mfcn, state.Trafo(), fStrategy);
72  MinimumParameters par(x, amin);
73  FunctionGradient gra = gc(par);
74  MinimumState tmp = (*this)(mfcn, MinimumState(par, MinimumError(MnAlgebraicSymMatrix(n), 1.), gra, state.Edm(), state.NFcn()), state.Trafo(), maxcalls);
75 
76  return MnUserParameterState(tmp, fcn.Up(), state.Trafo());
77 }
78 
79 void MnHesse::operator()(const FCNBase& fcn, FunctionMinimum& min, unsigned int maxcalls) const {
80  // interface from FunctionMinimum to be used after minimization
81  // use last state from the minimization without the need to re-create a new state
82  // do not reset function calls and keep updating them
83  MnUserFcn mfcn(fcn, min.UserState().Trafo(),min.NFcn());
84  MinimumState st = (*this)( mfcn, min.State(), min.UserState().Trafo(), maxcalls);
85  min.Add(st);
86 }
87 
88 MinimumState MnHesse::operator()(const MnFcn& mfcn, const MinimumState& st, const MnUserTransformation& trafo, unsigned int maxcalls) const {
89  // internal interface from MinimumState and MnUserTransformation
90  // Function who does the real Hessian calculations
91 
92  const MnMachinePrecision& prec = trafo.Precision();
93  // make sure starting at the right place
94  double amin = mfcn(st.Vec());
95  double aimsag = sqrt(prec.Eps2())*(fabs(amin)+mfcn.Up());
96 
97  // diagonal Elements first
98 
99  unsigned int n = st.Parameters().Vec().size();
100  if(maxcalls == 0) maxcalls = 200 + 100*n + 5*n*n;
101 
102  MnAlgebraicSymMatrix vhmat(n);
103  MnAlgebraicVector g2 = st.Gradient().G2();
104  MnAlgebraicVector gst = st.Gradient().Gstep();
105  MnAlgebraicVector grd = st.Gradient().Grad();
106  MnAlgebraicVector dirin = st.Gradient().Gstep();
107  MnAlgebraicVector yy(n);
108 
109 
110  // case gradient is not numeric (could be analytical or from FumiliGradientCalculator)
111 
112  if(st.Gradient().IsAnalytical() ) {
113  Numerical2PGradientCalculator igc(mfcn, trafo, fStrategy);
114  FunctionGradient tmp = igc(st.Parameters());
115  gst = tmp.Gstep();
116  dirin = tmp.Gstep();
117  g2 = tmp.G2();
118  }
119 
120  MnAlgebraicVector x = st.Parameters().Vec();
121 
122 #ifdef DEBUG
123  std::cout << "\nMnHesse " << std::endl;
124  std::cout << " x " << x << std::endl;
125  std::cout << " amin " << amin << " " << st.Fval() << std::endl;
126  std::cout << " grd " << grd << std::endl;
127  std::cout << " gst " << gst << std::endl;
128  std::cout << " g2 " << g2 << std::endl;
129  std::cout << " Gradient is analytical " << st.Gradient().IsAnalytical() << std::endl;
130 #endif
131 
132 
133  for(unsigned int i = 0; i < n; i++) {
134 
135  double xtf = x(i);
136  double dmin = 8.*prec.Eps2()*(fabs(xtf) + prec.Eps2());
137  double d = fabs(gst(i));
138  if(d < dmin) d = dmin;
139 
140 #ifdef DEBUG
141  std::cout << "\nDerivative parameter " << i << " d = " << d << " dmin = " << dmin << std::endl;
142 #endif
143 
144 
145  for(unsigned int icyc = 0; icyc < Ncycles(); icyc++) {
146  double sag = 0.;
147  double fs1 = 0.;
148  double fs2 = 0.;
149  for(unsigned int multpy = 0; multpy < 5; multpy++) {
150  x(i) = xtf + d;
151  fs1 = mfcn(x);
152  x(i) = xtf - d;
153  fs2 = mfcn(x);
154  x(i) = xtf;
155  sag = 0.5*(fs1+fs2-2.*amin);
156 
157 #ifdef DEBUG
158  std::cout << "cycle " << icyc << " mul " << multpy << "\t sag = " << sag << " d = " << d << std::endl;
159 #endif
160  // Now as F77 Minuit - check taht sag is not zero
161  if (sag != 0) goto L30; // break
162  if(trafo.Parameter(i).HasLimits()) {
163  if(d > 0.5) goto L26;
164  d *= 10.;
165  if(d > 0.5) d = 0.51;
166  continue;
167  }
168  d *= 10.;
169  }
170 
171 L26:
172 #ifdef WARNINGMSG
173 
174  // get parameter name for i
175  // (need separate scope for avoiding compl error when declaring name)
176  {
177  const char * name = trafo.Name( trafo.ExtOfInt(i));
178  MN_INFO_VAL2("MnHesse: 2nd derivative zero for Parameter ", name);
179  MN_INFO_MSG("MnHesse fails and will return diagonal matrix ");
180  }
181 #endif
182 
183  for(unsigned int j = 0; j < n; j++) {
184  double tmp = g2(j) < prec.Eps2() ? 1. : 1./g2(j);
185  vhmat(j,j) = tmp < prec.Eps2() ? 1. : tmp;
186  }
187 
188  return MinimumState(st.Parameters(), MinimumError(vhmat, MinimumError::MnHesseFailed()), st.Gradient(), st.Edm(), mfcn.NumOfCalls());
189 
190 L30:
191  double g2bfor = g2(i);
192  g2(i) = 2.*sag/(d*d);
193  grd(i) = (fs1-fs2)/(2.*d);
194  gst(i) = d;
195  dirin(i) = d;
196  yy(i) = fs1;
197  double dlast = d;
198  d = sqrt(2.*aimsag/fabs(g2(i)));
199  if(trafo.Parameter(i).HasLimits()) d = std::min(0.5, d);
200  if(d < dmin) d = dmin;
201 
202 #ifdef DEBUG
203  std::cout << "\t g1 = " << grd(i) << " g2 = " << g2(i) << " step = " << gst(i) << " d = " << d
204  << " diffd = " << fabs(d-dlast)/d << " diffg2 = " << fabs(g2(i)-g2bfor)/g2(i) << std::endl;
205 #endif
206 
207 
208  // see if converged
209  if(fabs((d-dlast)/d) < Tolerstp()) break;
210  if(fabs((g2(i)-g2bfor)/g2(i)) < TolerG2()) break;
211  d = std::min(d, 10.*dlast);
212  d = std::max(d, 0.1*dlast);
213  }
214  vhmat(i,i) = g2(i);
215  if(mfcn.NumOfCalls() > maxcalls) {
216 
217 #ifdef WARNINGMSG
218  //std::cout<<"maxcalls " << maxcalls << " " << mfcn.NumOfCalls() << " " << st.NFcn() << std::endl;
219  MN_INFO_MSG("MnHesse: maximum number of allowed function calls exhausted.");
220  MN_INFO_MSG("MnHesse fails and will return diagonal matrix ");
221 #endif
222 
223  for(unsigned int j = 0; j < n; j++) {
224  double tmp = g2(j) < prec.Eps2() ? 1. : 1./g2(j);
225  vhmat(j,j) = tmp < prec.Eps2() ? 1. : tmp;
226  }
227 
228  return MinimumState(st.Parameters(), MinimumError(vhmat, MinimumError::MnHesseFailed()), st.Gradient(), st.Edm(), mfcn.NumOfCalls());
229  }
230 
231  }
232 
233 #ifdef DEBUG
234  std::cout << "\n Second derivatives " << g2 << std::endl;
235 #endif
236 
237  if(fStrategy.Strategy() > 0) {
238  // refine first derivative
239  HessianGradientCalculator hgc(mfcn, trafo, fStrategy);
240  FunctionGradient gr = hgc(st.Parameters(), FunctionGradient(grd, g2, gst));
241  // update gradient and step values
242  grd = gr.Grad();
243  gst = gr.Gstep();
244  }
245 
246  //off-diagonal Elements
247  // initial starting values
248  if (n > 0) {
249  MPIProcess mpiprocOffDiagonal(n*(n-1)/2,0);
250  unsigned int startParIndexOffDiagonal = mpiprocOffDiagonal.StartElementIndex();
251  unsigned int endParIndexOffDiagonal = mpiprocOffDiagonal.EndElementIndex();
252 
253  unsigned int offsetVect = 0;
254  for (unsigned int in = 0; in<startParIndexOffDiagonal; in++)
255  if ((in+offsetVect)%(n-1)==0) offsetVect += (in+offsetVect)/(n-1);
256 
257  for (unsigned int in = startParIndexOffDiagonal;
258  in<endParIndexOffDiagonal; in++) {
259 
260  int i = (in+offsetVect)/(n-1);
261  if ((in+offsetVect)%(n-1)==0) offsetVect += i;
262  int j = (in+offsetVect)%(n-1)+1;
263 
264  if ((i+1)==j || in==startParIndexOffDiagonal)
265  x(i) += dirin(i);
266 
267  x(j) += dirin(j);
268 
269  double fs1 = mfcn(x);
270  double elem = (fs1 + amin - yy(i) - yy(j))/(dirin(i)*dirin(j));
271  vhmat(i,j) = elem;
272 
273  x(j) -= dirin(j);
274 
275  if (j%(n-1)==0 || in==endParIndexOffDiagonal-1)
276  x(i) -= dirin(i);
277 
278  }
279 
280  mpiprocOffDiagonal.SyncSymMatrixOffDiagonal(vhmat);
281  }
282 
283  //verify if matrix pos-def (still 2nd derivative)
284 
285 #ifdef DEBUG
286  std::cout << "Original error matrix " << vhmat << std::endl;
287 #endif
288 
289  MinimumError tmpErr = MnPosDef()(MinimumError(vhmat,1.), prec);
290 
291 #ifdef DEBUG
292  std::cout << "Original error matrix " << vhmat << std::endl;
293 #endif
294 
295  vhmat = tmpErr.InvHessian();
296 
297 #ifdef DEBUG
298  std::cout << "PosDef error matrix " << vhmat << std::endl;
299 #endif
300 
301 
302  int ifail = Invert(vhmat);
303  if(ifail != 0) {
304 
305 #ifdef WARNINGMSG
306  MN_INFO_MSG("MnHesse: matrix inversion fails!");
307  MN_INFO_MSG("MnHesse fails and will return diagonal matrix.");
308 #endif
309 
310  MnAlgebraicSymMatrix tmpsym(vhmat.Nrow());
311  for(unsigned int j = 0; j < n; j++) {
312  double tmp = g2(j) < prec.Eps2() ? 1. : 1./g2(j);
313  tmpsym(j,j) = tmp < prec.Eps2() ? 1. : tmp;
314  }
315 
316  return MinimumState(st.Parameters(), MinimumError(tmpsym, MinimumError::MnInvertFailed()), st.Gradient(), st.Edm(), mfcn.NumOfCalls());
317  }
318 
319  FunctionGradient gr(grd, g2, gst);
320  VariableMetricEDMEstimator estim;
321 
322  // if matrix is made pos def returns anyway edm
323  if(tmpErr.IsMadePosDef()) {
324  MinimumError err(vhmat, MinimumError::MnMadePosDef() );
325  double edm = estim.Estimate(gr, err);
326 #ifdef WARNINGMSG
327  MN_INFO_MSG("MnHesse: matrix was forced pos. def. ");
328 #endif
329  return MinimumState(st.Parameters(), err, gr, edm, mfcn.NumOfCalls());
330  }
331 
332  //calculate edm for good errors
333  MinimumError err(vhmat, 0.);
334  double edm = estim.Estimate(gr, err);
335 
336 #ifdef DEBUG
337  std::cout << "\nHesse is ACCURATE. New state from MnHesse " << std::endl;
338  std::cout << "Gradient " << grd << std::endl;
339  std::cout << "Second Deriv " << g2 << std::endl;
340  std::cout << "Gradient step " << gst << std::endl;
341  std::cout << "Error " << vhmat << std::endl;
342  std::cout << "edm " << edm << std::endl;
343 #endif
344 
345 
346  return MinimumState(st.Parameters(), err, gr, edm, mfcn.NumOfCalls());
347 }
348 
349 /*
350  MinimumError MnHesse::Hessian(const MnFcn& mfcn, const MinimumState& st, const MnUserTransformation& trafo) const {
351 
352  const MnMachinePrecision& prec = trafo.Precision();
353  // make sure starting at the right place
354  double amin = mfcn(st.Vec());
355  // if(fabs(amin - st.Fval()) > prec.Eps2()) std::cout<<"function Value differs from amin by "<<amin - st.Fval()<<std::endl;
356 
357  double aimsag = sqrt(prec.Eps2())*(fabs(amin)+mfcn.Up());
358 
359  // diagonal Elements first
360 
361  unsigned int n = st.Parameters().Vec().size();
362  MnAlgebraicSymMatrix vhmat(n);
363  MnAlgebraicVector g2 = st.Gradient().G2();
364  MnAlgebraicVector gst = st.Gradient().Gstep();
365  MnAlgebraicVector grd = st.Gradient().Grad();
366  MnAlgebraicVector dirin = st.Gradient().Gstep();
367  MnAlgebraicVector yy(n);
368  MnAlgebraicVector x = st.Parameters().Vec();
369 
370  for(unsigned int i = 0; i < n; i++) {
371 
372  double xtf = x(i);
373  double dmin = 8.*prec.Eps2()*fabs(xtf);
374  double d = fabs(gst(i));
375  if(d < dmin) d = dmin;
376  for(int icyc = 0; icyc < Ncycles(); icyc++) {
377  double sag = 0.;
378  double fs1 = 0.;
379  double fs2 = 0.;
380  for(int multpy = 0; multpy < 5; multpy++) {
381  x(i) = xtf + d;
382  fs1 = mfcn(x);
383  x(i) = xtf - d;
384  fs2 = mfcn(x);
385  x(i) = xtf;
386  sag = 0.5*(fs1+fs2-2.*amin);
387  if(sag > prec.Eps2()) break;
388  if(trafo.Parameter(i).HasLimits()) {
389  if(d > 0.5) {
390  std::cout<<"second derivative zero for Parameter "<<i<<std::endl;
391  std::cout<<"return diagonal matrix "<<std::endl;
392  for(unsigned int j = 0; j < n; j++) {
393  vhmat(j,j) = (g2(j) < prec.Eps2() ? 1. : 1./g2(j));
394  return MinimumError(vhmat, 1., false);
395  }
396  }
397  d *= 10.;
398  if(d > 0.5) d = 0.51;
399  continue;
400  }
401  d *= 10.;
402  }
403  if(sag < prec.Eps2()) {
404  std::cout<<"MnHesse: internal loop exhausted, return diagonal matrix."<<std::endl;
405  for(unsigned int i = 0; i < n; i++)
406  vhmat(i,i) = (g2(i) < prec.Eps2() ? 1. : 1./g2(i));
407  return MinimumError(vhmat, 1., false);
408  }
409  double g2bfor = g2(i);
410  g2(i) = 2.*sag/(d*d);
411  grd(i) = (fs1-fs2)/(2.*d);
412  gst(i) = d;
413  dirin(i) = d;
414  yy(i) = fs1;
415  double dlast = d;
416  d = sqrt(2.*aimsag/fabs(g2(i)));
417  if(trafo.Parameter(i).HasLimits()) d = std::min(0.5, d);
418  if(d < dmin) d = dmin;
419 
420  // see if converged
421  if(fabs((d-dlast)/d) < Tolerstp()) break;
422  if(fabs((g2(i)-g2bfor)/g2(i)) < TolerG2()) break;
423  d = std::min(d, 10.*dlast);
424  d = std::max(d, 0.1*dlast);
425  }
426  vhmat(i,i) = g2(i);
427  }
428 
429  //off-diagonal Elements
430  for(unsigned int i = 0; i < n; i++) {
431  x(i) += dirin(i);
432  for(unsigned int j = i+1; j < n; j++) {
433  x(j) += dirin(j);
434  double fs1 = mfcn(x);
435  double elem = (fs1 + amin - yy(i) - yy(j))/(dirin(i)*dirin(j));
436  vhmat(i,j) = elem;
437  x(j) -= dirin(j);
438  }
439  x(i) -= dirin(i);
440  }
441 
442  return MinimumError(vhmat, 0.);
443  }
444  */
445 
446  } // namespace Minuit2
447 
448 } // namespace ROOT