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"
18 #include "Minuit2/MinimumError.h"
19 #include "Minuit2/MinimumState.h"
20 #include "Minuit2/LaSum.h"
21 #include "Minuit2/MnPrint.h"
22 
23 #include <limits>
24 
25 namespace ROOT {
26 
27 namespace Minuit2 {
28 
29 double sum_of_elements(const LASymMatrix &);
30 
31 MinimumError
33 {
34  // dummy methods to suppress unused variable warnings
35  // this member function should never be called within
36  // the Fumili method...
37  s0.Fval();
38  p1.Fval();
39  g1.IsValid();
40  return MinimumError(2);
41 }
42 
44  const GradientCalculator &gc, double lambda) const
45 {
46  // calculate the error matrix using approximation of Fumili
47  // use only first derivatives (see tutorial par. 5.1,5.2)
48  // The Fumili Hessian is provided by the FumiliGRadientCalculator class
49  // we apply also the Marquard lambda factor to increase weight of diagonal term
50  // as suggester in Numerical Receipt for Marquard method
51 
52  // need to downcast to FumiliGradientCalculator
53  FumiliGradientCalculator *fgc = dynamic_cast<FumiliGradientCalculator *>(const_cast<GradientCalculator *>(&gc));
54  assert(fgc != 0);
55 
56  MnPrint print("FumiliErrorUpdator");
57 
58  // get Hessian from Gradient calculator
59 
61 
62  int nvar = p1.Vec().size();
63 
64  // apply Marquard lambda factor
65  double eps = 8 * std::numeric_limits<double>::min();
66  for (int j = 0; j < nvar; j++) {
67  h(j, j) *= (1. + lambda);
68  // if h(j,j) is zero what to do ?
69  if (std::fabs(h(j, j)) < eps) { // should use DBL_MIN
70  // put a cut off to avoid zero on diagonals
71  if (lambda > 1)
72  h(j, j) = lambda * eps;
73  else
74  h(j, j) = eps;
75  }
76  }
77 
78  int ifail = Invert(h);
79  if (ifail != 0) {
80  print.Warn("inversion fails; return diagonal matrix");
81 
82  for (unsigned int i = 0; i < h.Nrow(); i++) {
83  h(i, i) = 1. / h(i, i);
84  }
85  }
86 
87  const MnAlgebraicSymMatrix &v0 = s0.Error().InvHessian();
88 
89  // calculate by how much did the covariance matrix changed
90  // (if it changed a lot since the last step, probably we
91  // are not yet near the Minimum)
92  double dcov = 0.5 * (s0.Error().Dcovar() + sum_of_elements(h - v0) / sum_of_elements(h));
93 
94  return MinimumError(h, dcov);
95 }
96 
97 } // namespace Minuit2
98 
99 } // namespace ROOT
ROOT::Minuit2::FunctionGradient::IsValid
bool IsValid() const
Definition: FunctionGradient.h:42
ROOT::Minuit2::Invert
int Invert(LASymMatrix &)
Definition: LaInverse.cxx:21
v0
@ v0
Definition: rootcling_impl.cxx:3665
ROOT::Minuit2::MinimumState
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
Definition: MinimumState.h:29
MnStrategy.h
ROOT::Minuit2::LASymMatrix
Class describing a symmetric matrix of size n.
Definition: LASymMatrix.h:45
ROOT::Minuit2::MinimumParameters::Vec
const MnAlgebraicVector & Vec() const
Definition: MinimumParameters.h:38
MnMatrix.h
MnUserParameterState.h
ROOT::Minuit2::FunctionGradient
Definition: FunctionGradient.h:21
MnFcn.h
ROOT::Minuit2::MinimumParameters::Fval
double Fval() const
Definition: MinimumParameters.h:40
ROOT::Minuit2::FumiliGradientCalculator::Hessian
const MnAlgebraicSymMatrix & Hessian() const
Definition: FumiliGradientCalculator.h:39
ROOT::Minuit2::MnPrint::Warn
void Warn(const Ts &... args)
Definition: MnPrint.h:126
s0
#define s0(x)
Definition: RSha256.hxx:90
ROOT::Math::fabs
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
Definition: UnaryOperators.h:131
ROOT::Minuit2::MinimumParameters
Definition: MinimumParameters.h:21
FunctionGradient.h
FumiliErrorUpdator.h
h
#define h(i)
Definition: RSha256.hxx:106
ROOT::Minuit2::FumiliErrorUpdator::Update
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) ...
Definition: FumiliErrorUpdator.cxx:43
ROOT::Minuit2::FumiliGradientCalculator
Definition: FumiliGradientCalculator.h:23
FumiliGradientCalculator.h
ROOT::Minuit2::sum_of_elements
double sum_of_elements(const LASymMatrix &)
Definition: LaSumOfElements.cxx:26
ROOT::Minuit2::LAVector::size
unsigned int size() const
Definition: LAVector.h:227
ROOT::Minuit2::MinimumError
MinimumError keeps the inv.
Definition: MinimumError.h:26
MinimumError.h
MinimumParameters.h
LaSum.h
MnPrint.h
ROOT::Minuit2::GradientCalculator
interface class for gradient calculators
Definition: GradientCalculator.h:23
ROOT
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
Definition: EExecutionPolicy.hxx:4
MinimumState.h
ROOT::Minuit2::MnPrint
Definition: MnPrint.h:73