ROOT  6.06/09
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"
11 #include "Minuit2/MinimumState.h"
13 
14 #if defined(DEBUG) || defined(WARNINGMSG)
15 #include "Minuit2/MnPrint.h"
16 #endif
17 
18 #include <algorithm>
19 
20 
21 namespace ROOT {
22 
23  namespace Minuit2 {
24 
25 
26 LAVector 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 
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 
69  MnAlgebraicVector s(err.Nrow());
70  MnAlgebraicSymMatrix p(err.Nrow());
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
Namespace for new ROOT classes and functions.
Definition: ROOT.py:1
unsigned int size() const
Definition: LAVector.h:198
Class describing a symmetric matrix of size n.
Definition: LASymMatrix.h:51
MinimumState operator()(const MinimumState &, const MnMachinePrecision &) const
Definition: MnPosDef.cxx:29
determines the relative floating point arithmetic precision.
const FunctionGradient & Gradient() const
Definition: MinimumState.h:63
const MinimumError & Error() const
Definition: MinimumState.h:62
double sqrt(double)
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
double Eps2() const
eps2 returns 2*sqrt(eps)
const MnAlgebraicSymMatrix & InvHessian() const
Definition: MinimumError.h:60
const MinimumParameters & Parameters() const
Definition: MinimumState.h:58
double Eps() const
eps returns the smallest possible number so that 1.+eps > 1.
static Vc_ALWAYS_INLINE int_v max(const int_v &x, const int_v &y)
Definition: vector.h:440
MinimumError keeps the inv.
Definition: MinimumError.h:26
LAVector eigenvalues(const LASymMatrix &mat)
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
Definition: MinimumState.h:29