Logo ROOT  
Reference Guide
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 
11 #include "Minuit2/MinimumState.h"
12 #include "Minuit2/LaSum.h"
13 #include "Minuit2/LaProd.h"
14 #include "Minuit2/MnPrint.h"
15 
16 namespace ROOT {
17 
18 namespace Minuit2 {
19 
20 double inner_product(const LAVector &, const LAVector &);
21 double similarity(const LAVector &, const LASymMatrix &);
22 double sum_of_elements(const LASymMatrix &);
23 
24 MinimumError
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 /*
76 MinimumError 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
ROOT::Minuit2::DavidonErrorUpdator::Update
virtual MinimumError Update(const MinimumState &, const MinimumParameters &, const FunctionGradient &) const
Definition: DavidonErrorUpdator.cxx:25
v0
@ v0
Definition: rootcling_impl.cxx:3665
ROOT::Minuit2::MinimumState
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
Definition: MinimumState.h:29
ROOT::Minuit2::LAVector
Definition: LAVector.h:32
ROOT::Minuit2::MnPrint::Debug
void Debug(const Ts &... args)
Definition: MnPrint.h:138
ROOT::Minuit2::LASymMatrix
Class describing a symmetric matrix of size n.
Definition: LASymMatrix.h:45
ROOT::Minuit2::MinimumParameters::Vec
const MnAlgebraicVector & Vec() const
Definition: MinimumParameters.h:38
ROOT::Minuit2::FunctionGradient::Vec
const MnAlgebraicVector & Vec() const
Definition: FunctionGradient.h:41
ROOT::Minuit2::FunctionGradient
Definition: FunctionGradient.h:21
LaProd.h
ROOT::Minuit2::MnAlgebraicVector
LAVector MnAlgebraicVector
Definition: MnMatrix.h:28
ROOT::Minuit2::MnPrint::Warn
void Warn(const Ts &... args)
Definition: MnPrint.h:126
s0
#define s0(x)
Definition: RSha256.hxx:90
ROOT::Minuit2::MinimumParameters
Definition: MinimumParameters.h:21
ROOT::Minuit2::inner_product
double inner_product(const LAVector &, const LAVector &)
Definition: LaInnerProduct.cxx:18
ROOT::Minuit2::sum_of_elements
double sum_of_elements(const LASymMatrix &)
Definition: LaSumOfElements.cxx:26
ROOT::Minuit2::similarity
double similarity(const LAVector &, const LASymMatrix &)
Definition: LaVtMVSimilarity.cxx:20
ROOT::Minuit2::MinimumError
MinimumError keeps the inv.
Definition: MinimumError.h:26
ROOT::Minuit2::Outer_product
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;.
Definition: LaOuterProduct.h:26
LaSum.h
DavidonErrorUpdator.h
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