Logo ROOT  
Reference Guide
MnPosDef.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
10#include "Minuit2/MnPosDef.h"
13
14#if defined(DEBUG) || defined(WARNINGMSG)
15#include "Minuit2/MnPrint.h"
16#endif
17
18#include <algorithm>
19
20
21namespace ROOT {
22
23 namespace Minuit2 {
24
25
26LAVector eigenvalues(const LASymMatrix&);
27
28
30 // interface from minimum state
31 MinimumError err = (*this)(st.Error(), prec);
32 return MinimumState(st.Parameters(), err, st.Gradient(), st.Edm(), st.NFcn());
33}
34
36 // make error matrix positive defined returning a new corrected minimum error state
37
38 MnAlgebraicSymMatrix err(e.InvHessian());
39 if(err.size() == 1 && err(0,0) < prec.Eps()) {
40 err(0,0) = 1.;
42 }
43 if(err.size() == 1 && err(0,0) > prec.Eps()) {
44 return e;
45 }
46 // std::cout<<"MnPosDef init matrix= "<<err<<std::endl;
47
48 double epspdf = std::max(1.e-6, prec.Eps2());
49 double dgmin = err(0,0);
50
51 for(unsigned int i = 0; i < err.Nrow(); i++) {
52#ifdef WARNINGMSG
53 if(err(i,i) <= 0 /* prec.Eps2() */ )
54 MN_INFO_VAL2("negative or zero diagonal element in covariance matrix",i);
55#endif
56 if(err(i,i) < dgmin) dgmin = err(i,i);
57 }
58 double dg = 0.;
59 if(dgmin <= 0) {
60 //dg = 1. + epspdf - dgmin;
61 dg = 0.5 + epspdf - dgmin;
62 // dg = 0.5*(1. + epspdf - dgmin);
63#ifdef WARNINGMSG
64 MN_INFO_VAL2("added to diagonal of Error matrix a value",dg);
65#endif
66 //std::cout << "Error matrix " << err << std::endl;
67 }
68
71 for(unsigned int i = 0; i < err.Nrow(); i++) {
72 err(i,i) += dg;
73 if(err(i,i) < 0.) err(i,i) = 1.;
74 s(i) = 1./sqrt(err(i,i));
75 for(unsigned int j = 0; j <= i; j++) {
76 p(i,j) = err(i,j)*s(i)*s(j);
77 }
78 }
79
80 //std::cout<<"MnPosDef p: "<<p<<std::endl;
82 double pmin = eval(0);
83 double pmax = eval(eval.size() - 1);
84 //std::cout<<"pmin= "<<pmin<<" pmax= "<<pmax<<std::endl;
85 pmax = std::max(fabs(pmax), 1.);
86 if(pmin > epspdf*pmax) return MinimumError(err, e.Dcovar());
87
88 double padd = 0.001*pmax - pmin;
89#ifdef DEBUG
90 std::cout<<"eigenvalues: "<<std::endl;
91#endif
92 for(unsigned int i = 0; i < err.Nrow(); i++) {
93 err(i,i) *= (1. + padd);
94#ifdef DEBUG
95 std::cout<<eval(i)<<std::endl;
96#endif
97 }
98 // std::cout<<"MnPosDef final matrix: "<<err<<std::endl;
99#ifdef WARNINGMSG
100 MN_INFO_VAL2("matrix forced pos-def by adding to diagonal",padd);
101#endif
103}
104
105 } // namespace Minuit2
106
107} // namespace ROOT
#define MN_INFO_VAL2(loc, x)
Definition: MnPrint.h:130
#define e(i)
Definition: RSha256.hxx:103
double sqrt(double)
Class describing a symmetric matrix of size n.
Definition: LASymMatrix.h:51
unsigned int Nrow() const
Definition: LASymMatrix.h:239
unsigned int size() const
Definition: LASymMatrix.h:237
unsigned int size() const
Definition: LAVector.h:198
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
const MinimumError & Error() const
Definition: MinimumState.h:62
const MinimumParameters & Parameters() const
Definition: MinimumState.h:58
const FunctionGradient & Gradient() const
Definition: MinimumState.h:63
determines the relative floating point arithmetic precision.
double Eps() const
eps returns the smallest possible number so that 1.+eps > 1.
double Eps2() const
eps2 returns 2*sqrt(eps)
MinimumState operator()(const MinimumState &, const MnMachinePrecision &) const
Definition: MnPosDef.cxx:29
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
LAVector eigenvalues(const LASymMatrix &mat)
VSD Structures.
Definition: StringConv.hxx:21
static constexpr double s