Logo ROOT   6.10/09
Reference Guide
FumiliErrorUpdator.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/MnFcn.h"
12 #include "Minuit2/MnStrategy.h"
17 #include "Minuit2/MnMatrix.h"
18 #include "Minuit2/MinimumError.h"
19 #include "Minuit2/MinimumState.h"
20 #include "Minuit2/LaSum.h"
21 #include <limits>
22 
23 #if defined(DEBUG) || defined(WARNINGMSG)
24 #include "Minuit2/MnPrint.h"
25 #endif
26 
27 
28 namespace ROOT {
29 
30 namespace Minuit2 {
31 
32 
33 
34 double sum_of_elements(const LASymMatrix&);
35 
36 
38  const MinimumParameters& p1,
39  const FunctionGradient& g1) const {
40  // dummy methods to suppress unused variable warnings
41  // this member function should never be called within
42  // the Fumili method...
43  s0.Fval();
44  p1.Fval();
45  g1.IsValid();
46  return MinimumError(2);
47 }
48 
49 
51  const MinimumParameters& p1,
52  const GradientCalculator& gc ,
53  double lambda) const {
54  // calculate the error matrix using approximation of Fumili
55  // use only first derivatives (see tutorial par. 5.1,5.2)
56  // The Fumili Hessian is provided by the FumiliGRadientCalculator class
57  // we apply also the Marquard lambda factor to increase weight of diagonal term
58  // as suggester in Numerical Receipt for Marquard method
59 
60  // need to downcast to FumiliGradientCalculator
61  FumiliGradientCalculator * fgc = dynamic_cast< FumiliGradientCalculator *>( const_cast<GradientCalculator *>(&gc) );
62  assert(fgc != 0);
63 
64 
65  // get Hessian from Gradient calculator
66 
68 
69  int nvar = p1.Vec().size();
70 
71  // apply Marquard lambda factor
72  double eps = 8*std::numeric_limits<double>::min();
73  for (int j = 0; j < nvar; j++) {
74  h(j,j) *= (1. + lambda);
75  // if h(j,j) is zero what to do ?
76  if ( fabs( h(j,j) ) < eps ) { // should use DBL_MIN
77  // put a cut off to avoid zero on diagonals
78  if ( lambda > 1)
79  h(j,j) = lambda*eps;
80  else
81  h(j,j) = eps;
82  }
83  }
84 
85 
86 
87  int ifail = Invert(h);
88  if(ifail != 0) {
89 #ifdef WARNINGMSG
90  MN_INFO_MSG("FumiliErrorUpdator inversion fails; return diagonal matrix.");
91 #endif
92  for(unsigned int i = 0; i < h.Nrow(); i++) {
93  h(i,i) = 1./h(i,i);
94  }
95  }
96 
97 
98  const MnAlgebraicSymMatrix& v0 = s0.Error().InvHessian();
99 
100  // calculate by how much did the covariance matrix changed
101  // (if it changed a lot since the last step, probably we
102  // are not yet near the Minimum)
103  double dcov = 0.5*(s0.Error().Dcovar() + sum_of_elements(h-v0)/sum_of_elements(h));
104 
105 
106 
107  return MinimumError(h, dcov);
108 
109 }
110 
111 
112 } // namespace Minuit2
113 
114 } // namespace ROOT
int Invert(LASymMatrix &)
Definition: LaInverse.cxx:22
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
const MnAlgebraicVector & Vec() const
double sum_of_elements(const LASymMatrix &)
TH1 * h
Definition: legend2.C:5
Class describing a symmetric matrix of size n.
Definition: LASymMatrix.h:51
#define MN_INFO_MSG(str)
Definition: MnPrint.h:110
virtual MinimumError Update(const MinimumState &fMinimumState, const MinimumParameters &fMinimumParameters, const GradientCalculator &fGradientCalculator, double lambda) const
Member function that calculates the Error matrix (or the Hessian matrix containing the (approximate) ...
const MnAlgebraicSymMatrix & Hessian() const
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
unsigned int Nrow() const
Definition: LASymMatrix.h:239
unsigned int size() const
Definition: LAVector.h:198
static double p1(double t, double a, double b)
const MinimumError & Error() const
Definition: MinimumState.h:62
MinimumError keeps the inv.
Definition: MinimumError.h:26
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
Definition: MinimumState.h:29
interface class for gradient calculators
const MnAlgebraicSymMatrix & InvHessian() const
Definition: MinimumError.h:60