Logo ROOT   6.10/09
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 
15 //#define DEBUG
16 
17 #if defined(DEBUG) || defined(WARNINGMSG)
18 #include "Minuit2/MnPrint.h"
19 #endif
20 
21 
22 
23 namespace ROOT {
24 
25  namespace Minuit2 {
26 
27 
28 double inner_product(const LAVector&, const LAVector&);
29 double similarity(const LAVector&, const LASymMatrix&);
30 double sum_of_elements(const LASymMatrix&);
31 
33  const MinimumParameters& p1,
34  const FunctionGradient& g1) const {
35 
36  // update of the covarianze matrix (Davidon formula, see Tutorial, par. 4.8 pag 26)
37  // in case of delgam > gvg (PHI > 1) use rank one formula
38  // see par 4.10 pag 30
39 
40  const MnAlgebraicSymMatrix& v0 = s0.Error().InvHessian();
41  MnAlgebraicVector dx = p1.Vec() - s0.Vec();
42  MnAlgebraicVector dg = g1.Vec() - s0.Gradient().Vec();
43 
44  double delgam = inner_product(dx, dg);
45  double gvg = similarity(dg, v0);
46 
47 
48 #ifdef DEBUG
49  std::cout << "dx = " << dx << std::endl;
50  std::cout << "dg = " << dg << std::endl;
51  std::cout<<"delgam= "<<delgam<<" gvg= "<<gvg<<std::endl;
52 #endif
53 
54  if (delgam == 0 ) {
55 #ifdef WARNINGMSG
56  MN_INFO_MSG("DavidonErrorUpdator: delgam = 0 : cannot update - return same matrix ");
57 #endif
58  return s0.Error();
59  }
60 #ifdef WARNINGMSG
61  if (delgam < 0) MN_INFO_MSG("DavidonErrorUpdator: delgam < 0 : first derivatives increasing along search line");
62 #endif
63 
64  if (gvg <= 0 ) {
65  // since v0 is pos def this gvg can be only = 0 if dg = 0 - should never be here
66 #ifdef WARNINGMSG
67  MN_INFO_MSG("DavidonErrorUpdator: gvg <= 0 : cannot update - return same matrix ");
68 #endif
69  return s0.Error();
70  }
71 
72 
73  MnAlgebraicVector vg = v0*dg;
74 
75  MnAlgebraicSymMatrix vUpd = Outer_product(dx)/delgam - Outer_product(vg)/gvg;
76 
77  if(delgam > gvg) {
78  // use rank 1 formula
79  vUpd += gvg*Outer_product(MnAlgebraicVector(dx/delgam - vg/gvg));
80  }
81 
82  double sum_upd = sum_of_elements(vUpd);
83  vUpd += v0;
84 
85  double dcov = 0.5*(s0.Error().Dcovar() + sum_upd/sum_of_elements(vUpd));
86 
87  return MinimumError(vUpd, dcov);
88 }
89 
90 /*
91 MinimumError DavidonErrorUpdator::Update(const MinimumState& s0,
92  const MinimumParameters& p1,
93  const FunctionGradient& g1) const {
94 
95  const MnAlgebraicSymMatrix& v0 = s0.Error().InvHessian();
96  MnAlgebraicVector dx = p1.Vec() - s0.Vec();
97  MnAlgebraicVector dg = g1.Vec() - s0.Gradient().Vec();
98 
99  double delgam = inner_product(dx, dg);
100  double gvg = similarity(dg, v0);
101 
102 // std::cout<<"delgam= "<<delgam<<" gvg= "<<gvg<<std::endl;
103  MnAlgebraicVector vg = v0*dg;
104 // MnAlgebraicSymMatrix vUpd(v0.Nrow());
105 
106 // MnAlgebraicSymMatrix dd = ( 1./delgam )*outer_product(dx);
107 // dd *= ( 1./delgam );
108 // MnAlgebraicSymMatrix VggV = ( 1./gvg )*outer_product(vg);
109 // VggV *= ( 1./gvg );
110 // vUpd = dd - VggV;
111 // MnAlgebraicSymMatrix vUpd = ( 1./delgam )*outer_product(dx) - ( 1./gvg )*outer_product(vg);
112  MnAlgebraicSymMatrix vUpd = Outer_product(dx)/delgam - Outer_product(vg)/gvg;
113 
114  if(delgam > gvg) {
115 // dx *= ( 1./delgam );
116 // vg *= ( 1./gvg );
117 // MnAlgebraicVector flnu = dx - vg;
118 // MnAlgebraicSymMatrix tmp = Outer_product(flnu);
119 // tmp *= gvg;
120 // vUpd = vUpd + tmp;
121  vUpd += gvg*outer_product(dx/delgam - vg/gvg);
122  }
123 
124 //
125 // MnAlgebraicSymMatrix dd = Outer_product(dx);
126 // dd *= ( 1./delgam );
127 // MnAlgebraicSymMatrix VggV = Outer_product(vg);
128 // VggV *= ( 1./gvg );
129 // vUpd = dd - VggV;
130 //
131 //
132 // double phi = delgam/(delgam - gvg);
133 
134 // MnAlgebraicSymMatrix vUpd(v0.Nrow());
135 // if(phi < 0) {
136 // // rank-2 Update
137 // MnAlgebraicSymMatrix dd = Outer_product(dx);
138 // dd *= ( 1./delgam );
139 // MnAlgebraicSymMatrix VggV = Outer_product(vg);
140 // VggV *= ( 1./gvg );
141 // vUpd = dd - VggV;
142 // }
143 // if(phi > 1) {
144 // // rank-1 Update
145 // MnAlgebraicVector tmp = dx - vg;
146 // vUpd = Outer_product(tmp);
147 // vUpd *= ( 1./(delgam - gvg) );
148 // }
149 //
150 
151 //
152 // if(delgam > gvg) {
153 // // rank-1 Update
154 // MnAlgebraicVector tmp = dx - vg;
155 // vUpd = Outer_product(tmp);
156 // vUpd *= ( 1./(delgam - gvg) );
157 // } else {
158 // // rank-2 Update
159 // MnAlgebraicSymMatrix dd = Outer_product(dx);
160 // dd *= ( 1./delgam );
161 // MnAlgebraicSymMatrix VggV = Outer_productn(vg);
162 // VggV *= ( 1./gvg );
163 // vUpd = dd - VggV;
164 // }
165 //
166 
167  double sum_upd = sum_of_elements(vUpd);
168  vUpd += v0;
169 
170 // MnAlgebraicSymMatrix V1 = v0 + vUpd;
171 
172  double dcov =
173  0.5*(s0.Error().Dcovar() + sum_upd/sum_of_elements(vUpd));
174 
175  return MinimumError(vUpd, dcov);
176 }
177 */
178 
179  } // namespace Minuit2
180 
181 } // namespace ROOT
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
const MnAlgebraicVector & Vec() const
double sum_of_elements(const LASymMatrix &)
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;.
Class describing a symmetric matrix of size n.
Definition: LASymMatrix.h:51
const MnAlgebraicVector & Vec() const
double inner_product(const LAVector &, const LAVector &)
#define MN_INFO_MSG(str)
Definition: MnPrint.h:110
const FunctionGradient & Gradient() const
Definition: MinimumState.h:63
LAVector MnAlgebraicVector
Definition: MnMatrix.h:42
static double p1(double t, double a, double b)
const MinimumError & Error() const
Definition: MinimumState.h:62
virtual MinimumError Update(const MinimumState &, const MinimumParameters &, const FunctionGradient &) const
double similarity(const LAVector &, const LASymMatrix &)
MinimumError keeps the inv.
Definition: MinimumError.h:26
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
Definition: MinimumState.h:29
const MnAlgebraicVector & Vec() const
Definition: MinimumState.h:59
const MnAlgebraicSymMatrix & InvHessian() const
Definition: MinimumError.h:60