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 "Minuit2/MnPrint.h"
22
23#include <limits>
24
25namespace ROOT {
26
27namespace Minuit2 {
28
29double sum_of_elements(const LASymMatrix &);
30
31MinimumError
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
#define s0(x)
Definition: RSha256.hxx:90
#define h(i)
Definition: RSha256.hxx:106
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void gc
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:45
unsigned int size() const
Definition: LAVector.h:227
MinimumError keeps the inv.
Definition: MinimumError.h:28
const MnAlgebraicVector & Vec() const
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
Definition: MinimumState.h:27
void Warn(const Ts &... args)
Definition: MnPrint.h:126
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
int Invert(LASymMatrix &)
Definition: LaInverse.cxx:21
double sum_of_elements(const LASymMatrix &)
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.