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"
11 #include "Minuit2/MinimumState.h"
13 #include "Minuit2/MnPrint.h"
14 
15 #include <algorithm>
16 
17 namespace ROOT {
18 
19 namespace Minuit2 {
20 
21 LAVector eigenvalues(const LASymMatrix &);
22 
24 {
25  // interface from minimum state
26  MinimumError err = (*this)(st.Error(), prec);
27  return MinimumState(st.Parameters(), err, st.Gradient(), st.Edm(), st.NFcn());
28 }
29 
31 {
32  MnPrint print("MnPosDef");
33 
34  // make error matrix positive defined returning a new corrected minimum error state
35 
36  MnAlgebraicSymMatrix err(e.InvHessian());
37  if (err.size() == 1 && err(0, 0) < prec.Eps()) {
38  err(0, 0) = 1.;
40  }
41  if (err.size() == 1 && err(0, 0) > prec.Eps()) {
42  return e;
43  }
44  // std::cout<<"MnPosDef init matrix= "<<err<<std::endl;
45 
46  double epspdf = std::max(1.e-6, prec.Eps2());
47  double dgmin = err(0, 0);
48 
49  for (unsigned int i = 0; i < err.Nrow(); i++) {
50  if (err(i, i) <= 0 /* prec.Eps2() */)
51  print.Warn("non-positive diagonal element in covariance matrix[", i, "] =", err(i, i));
52 
53  if (err(i, i) < dgmin)
54  dgmin = err(i, i);
55  }
56  double dg = 0.;
57  if (dgmin <= 0) {
58  // dg = 1. + epspdf - dgmin;
59  dg = 0.5 + epspdf - dgmin;
60  // dg = 0.5*(1. + epspdf - dgmin);
61 
62  print.Warn("Added to diagonal of Error matrix a value", dg);
63 
64  // std::cout << "Error matrix " << err << std::endl;
65  }
66 
67  MnAlgebraicVector s(err.Nrow());
68  MnAlgebraicSymMatrix p(err.Nrow());
69  for (unsigned int i = 0; i < err.Nrow(); i++) {
70  err(i, i) += dg;
71  if (err(i, i) < 0.)
72  err(i, i) = 1.;
73  s(i) = 1. / std::sqrt(err(i, i));
74  for (unsigned int j = 0; j <= i; j++) {
75  p(i, j) = err(i, j) * s(i) * s(j);
76  }
77  }
78 
79  // std::cout<<"MnPosDef p: "<<p<<std::endl;
81  double pmin = eval(0);
82  double pmax = eval(eval.size() - 1);
83  // std::cout<<"pmin= "<<pmin<<" pmax= "<<pmax<<std::endl;
84  pmax = std::max(std::fabs(pmax), 1.);
85  if (pmin > epspdf * pmax)
86  return MinimumError(err, e.Dcovar());
87 
88  double padd = 0.001 * pmax - pmin;
89 
90  for (unsigned int i = 0; i < err.Nrow(); i++)
91  err(i, i) *= (1. + padd);
92 
93  print.Debug([&](std::ostream &os) {
94  os << "Eigenvalues:";
95  for (unsigned i = 0; i < err.Nrow(); ++i)
96  os << "\n " << eval(i);
97  });
98 
99  // std::cout<<"MnPosDef final matrix: "<<err<<std::endl;
100 
101  print.Warn("Matrix forced pos-def by adding to diagonal", padd);
102 
104 }
105 
106 } // namespace Minuit2
107 
108 } // namespace ROOT
ROOT::Minuit2::MnMachinePrecision::Eps2
double Eps2() const
eps2 returns 2*sqrt(eps)
Definition: MnMachinePrecision.h:41
ROOT::Minuit2::MinimumState::Edm
double Edm() const
Definition: MinimumState.h:57
ROOT::Minuit2::BasicMinimumError::MnMadePosDef
Definition: BasicMinimumError.h:35
e
#define e(i)
Definition: RSha256.hxx:103
ROOT::Minuit2::MinimumState
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
Definition: MinimumState.h:29
MnMachinePrecision.h
ROOT::Minuit2::LAVector
Definition: LAVector.h:32
ROOT::Minuit2::MnPrint::Debug
void Debug(const Ts &... args)
Definition: MnPrint.h:138
ROOT::Minuit2::eigenvalues
LAVector eigenvalues(const LASymMatrix &mat)
Definition: LaEigenValues.cxx:19
ROOT::Minuit2::MinimumState::Gradient
const FunctionGradient & Gradient() const
Definition: MinimumState.h:55
ROOT::Minuit2::LASymMatrix
Class describing a symmetric matrix of size n.
Definition: LASymMatrix.h:45
ROOT::Minuit2::LASymMatrix::size
unsigned int size() const
Definition: LASymMatrix.h:272
TGeant4Unit::s
static constexpr double s
Definition: TGeant4SystemOfUnits.h:162
ROOT::Minuit2::MinimumState::NFcn
int NFcn() const
Definition: MinimumState.h:58
ROOT::Minuit2::MnPosDef::operator()
MinimumState operator()(const MinimumState &, const MnMachinePrecision &) const
Definition: MnPosDef.cxx:23
ROOT::Minuit2::MnMachinePrecision::Eps
double Eps() const
eps returns the smallest possible number so that 1.+eps > 1.
Definition: MnMachinePrecision.h:38
ROOT::Minuit2::MnPrint::Warn
void Warn(const Ts &... args)
Definition: MnPrint.h:126
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::MinimumState::Error
const MinimumError & Error() const
Definition: MinimumState.h:54
ROOT::Minuit2::MnMachinePrecision
Sets the relative floating point (double) arithmetic precision.
Definition: MnMachinePrecision.h:32
sqrt
double sqrt(double)
MnPosDef.h
ROOT::Minuit2::LASymMatrix::Nrow
unsigned int Nrow() const
Definition: LASymMatrix.h:274
ROOT::Minuit2::LAVector::size
unsigned int size() const
Definition: LAVector.h:227
ROOT::Minuit2::MinimumError
MinimumError keeps the inv.
Definition: MinimumError.h:26
ROOT::Minuit2::MinimumState::Parameters
const MinimumParameters & Parameters() const
Definition: MinimumState.h:50
MnPrint.h
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