Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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
67 MnAlgebraicVector s(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
#define e(i)
Definition RSha256.hxx:103
Class describing a symmetric matrix of size n.
Definition LASymMatrix.h:45
unsigned int Nrow() const
unsigned int size() const
unsigned int size() const
Definition LAVector.h:227
MinimumError keeps the inv.
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
const MinimumError & Error() const
const MinimumParameters & Parameters() const
const FunctionGradient & Gradient() const
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:138
void Warn(const Ts &... args)
Definition MnPrint.h:126
LAVector eigenvalues(const LASymMatrix &mat)
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...