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
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) ...
Fumili gradient calculator using external gradient provided by FCN Note that the computed Hessian and...
const MnAlgebraicSymMatrix & GetHessian() const
Class describing a symmetric matrix of size n.
Definition: LASymMatrix.h:45
unsigned int size() const
Definition: LAVector.h:231
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:135
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.