Logo ROOT  
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"
20#include "Minuit2/LaSum.h"
21#include <limits>
22
23#if defined(DEBUG) || defined(WARNINGMSG)
24#include "Minuit2/MnPrint.h"
25#endif
26
27
28namespace ROOT {
29
30namespace Minuit2 {
31
32
33
34double 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
#define MN_INFO_MSG(str)
Definition: MnPrint.h:110
#define s0(x)
Definition: RSha256.hxx:90
#define h(i)
Definition: RSha256.hxx:106
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
interface class for gradient calculators
Class describing a symmetric matrix of size n.
Definition: LASymMatrix.h:51
unsigned int size() const
Definition: LAVector.h:198
MinimumError keeps the inv.
Definition: MinimumError.h:26
const MnAlgebraicVector & Vec() const
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
Definition: MinimumState.h:29
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
int Invert(LASymMatrix &)
Definition: LaInverse.cxx:22
double sum_of_elements(const LASymMatrix &)
VSD Structures.
Definition: StringConv.hxx:21