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#include "Minuit2/MnPrint.h"
14
15#include <algorithm>
16
17namespace ROOT {
18
19namespace Minuit2 {
20
21LAVector 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
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
102
104}
105
106} // namespace Minuit2
107
108} // namespace ROOT
#define e(i)
Definition: RSha256.hxx:103
winID h TVirtualViewer3D TVirtualGLPainter p
Class describing a symmetric matrix of size n.
Definition: LASymMatrix.h:45
unsigned int Nrow() const
Definition: LASymMatrix.h:273
unsigned int size() const
Definition: LASymMatrix.h:271
unsigned int size() const
Definition: LAVector.h:231
MinimumError keeps the inv.
Definition: MinimumError.h:28
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
Definition: MinimumState.h:27
const MinimumError & Error() const
Definition: MinimumState.h:61
const MinimumParameters & Parameters() const
Definition: MinimumState.h:57
Definition: MinimumState.h:62
Sets the relative floating point (double) 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:23
void Debug(const Ts &... args)
Definition: MnPrint.h:147
void Warn(const Ts &... args)
Definition: MnPrint.h:135
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
LAVector eigenvalues(const LASymMatrix &mat)
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.
static constexpr double s