ROOT   Reference Guide
Searching...
No Matches
DavidonErrorUpdator.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
12#include "Minuit2/LaSum.h"
13#include "Minuit2/LaProd.h"
14#include "Minuit2/MnPrint.h"
15
16namespace ROOT {
17
18namespace Minuit2 {
19
20double inner_product(const LAVector &, const LAVector &);
21double similarity(const LAVector &, const LASymMatrix &);
22double sum_of_elements(const LASymMatrix &);
23
24MinimumError
26{
27
28 // update of the covarianze matrix (Davidon formula, see Tutorial, par. 4.8 pag 26)
29 // in case of delgam > gvg (PHI > 1) use rank one formula
30 // see par 4.10 pag 30
31
32 MnPrint print("DavidonErrorUpdator");
33
34 const MnAlgebraicSymMatrix &v0 = s0.Error().InvHessian();
35 MnAlgebraicVector dx = p1.Vec() - s0.Vec();
36 MnAlgebraicVector dg = g1.Vec() - s0.Gradient().Vec();
37
38 double delgam = inner_product(dx, dg);
39 double gvg = similarity(dg, v0);
40
41 print.Debug("\ndx", dx, "\ndg", dg, "\ndelgam", delgam, "gvg", gvg);
42
43 if (delgam == 0) {
44 print.Warn("delgam = 0 : cannot update - return same matrix");
45 return s0.Error();
46 }
47
48 if (delgam < 0) {
49 print.Warn("delgam < 0 : first derivatives increasing along search line");
50 }
51
52 if (gvg <= 0) {
53 // since v0 is pos def this gvg can be only = 0 if dg = 0 - should never be here
54 print.Warn("gvg <= 0 : cannot update - return same matrix");
55 return s0.Error();
56 }
57
58 MnAlgebraicVector vg = v0 * dg;
59
60 MnAlgebraicSymMatrix vUpd = Outer_product(dx) / delgam - Outer_product(vg) / gvg;
61
62 if (delgam > gvg) {
63 // use rank 1 formula
64 vUpd += gvg * Outer_product(MnAlgebraicVector(dx / delgam - vg / gvg));
65 }
66
67 double sum_upd = sum_of_elements(vUpd);
68 vUpd += v0;
69
70 double dcov = 0.5 * (s0.Error().Dcovar() + sum_upd / sum_of_elements(vUpd));
71
72 return MinimumError(vUpd, dcov);
73}
74
75/*
76MinimumError DavidonErrorUpdator::Update(const MinimumState& s0,
77 const MinimumParameters& p1,
78 const FunctionGradient& g1) const {
79
80 const MnAlgebraicSymMatrix& v0 = s0.Error().InvHessian();
81 MnAlgebraicVector dx = p1.Vec() - s0.Vec();
82 MnAlgebraicVector dg = g1.Vec() - s0.Gradient().Vec();
83
84 double delgam = inner_product(dx, dg);
85 double gvg = similarity(dg, v0);
86
87// std::cout<<"delgam= "<<delgam<<" gvg= "<<gvg<<std::endl;
88 MnAlgebraicVector vg = v0*dg;
89// MnAlgebraicSymMatrix vUpd(v0.Nrow());
90
91// MnAlgebraicSymMatrix dd = ( 1./delgam )*outer_product(dx);
92// dd *= ( 1./delgam );
93// MnAlgebraicSymMatrix VggV = ( 1./gvg )*outer_product(vg);
94// VggV *= ( 1./gvg );
95// vUpd = dd - VggV;
96// MnAlgebraicSymMatrix vUpd = ( 1./delgam )*outer_product(dx) - ( 1./gvg )*outer_product(vg);
97 MnAlgebraicSymMatrix vUpd = Outer_product(dx)/delgam - Outer_product(vg)/gvg;
98
99 if(delgam > gvg) {
100// dx *= ( 1./delgam );
101// vg *= ( 1./gvg );
102// MnAlgebraicVector flnu = dx - vg;
103// MnAlgebraicSymMatrix tmp = Outer_product(flnu);
104// tmp *= gvg;
105// vUpd = vUpd + tmp;
106 vUpd += gvg*outer_product(dx/delgam - vg/gvg);
107 }
108
109//
110// MnAlgebraicSymMatrix dd = Outer_product(dx);
111// dd *= ( 1./delgam );
112// MnAlgebraicSymMatrix VggV = Outer_product(vg);
113// VggV *= ( 1./gvg );
114// vUpd = dd - VggV;
115//
116//
117// double phi = delgam/(delgam - gvg);
118
119// MnAlgebraicSymMatrix vUpd(v0.Nrow());
120// if(phi < 0) {
121// // rank-2 Update
122// MnAlgebraicSymMatrix dd = Outer_product(dx);
123// dd *= ( 1./delgam );
124// MnAlgebraicSymMatrix VggV = Outer_product(vg);
125// VggV *= ( 1./gvg );
126// vUpd = dd - VggV;
127// }
128// if(phi > 1) {
129// // rank-1 Update
130// MnAlgebraicVector tmp = dx - vg;
131// vUpd = Outer_product(tmp);
132// vUpd *= ( 1./(delgam - gvg) );
133// }
134//
135
136//
137// if(delgam > gvg) {
138// // rank-1 Update
139// MnAlgebraicVector tmp = dx - vg;
140// vUpd = Outer_product(tmp);
141// vUpd *= ( 1./(delgam - gvg) );
142// } else {
143// // rank-2 Update
144// MnAlgebraicSymMatrix dd = Outer_product(dx);
145// dd *= ( 1./delgam );
146// MnAlgebraicSymMatrix VggV = Outer_productn(vg);
147// VggV *= ( 1./gvg );
148// vUpd = dd - VggV;
149// }
150//
151
152 double sum_upd = sum_of_elements(vUpd);
153 vUpd += v0;
154
155// MnAlgebraicSymMatrix V1 = v0 + vUpd;
156
157 double dcov =
158 0.5*(s0.Error().Dcovar() + sum_upd/sum_of_elements(vUpd));
159
160 return MinimumError(vUpd, dcov);
161}
162*/
163
164} // namespace Minuit2
165
166} // namespace ROOT
#define s0(x)
Definition RSha256.hxx:90
virtual MinimumError Update(const MinimumState &, const MinimumParameters &, const FunctionGradient &) const
const MnAlgebraicVector & Vec() const
Class describing a symmetric matrix of size n.
Definition LASymMatrix.h:45
MinimumError keeps the inv.
const MnAlgebraicVector & Vec() const
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
void Debug(const Ts &... args)
Definition MnPrint.h:138
void Warn(const Ts &... args)
Definition MnPrint.h:126
ABObj< sym, VectorOuterProduct< ABObj< vec, LAVector, double >, double >, double > Outer_product(const ABObj< vec, LAVector, double > &obj)
LAPACK Algebra function specialize the Outer_product function for LAVector;.
LAVector MnAlgebraicVector
Definition MnMatrix.h:28
double sum_of_elements(const LASymMatrix &)
double similarity(const LAVector &, const LASymMatrix &)
double inner_product(const LAVector &, const LAVector &)
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...