Logo ROOT  
Reference Guide
HessianGradientCalculator.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/MnFcn.h"
17 #include "Minuit2/MnStrategy.h"
18 #include "Minuit2/MnPrint.h"
19 #include "Minuit2/MPIProcess.h"
20 
21 #include <cmath>
22 #include <cassert>
23 
24 namespace ROOT {
25 
26 namespace Minuit2 {
27 
29 {
30  // use initial gradient as starting point
32  FunctionGradient gra = gc(par);
33 
34  return (*this)(par, gra);
35 }
36 
38 operator()(const MinimumParameters &par, const FunctionGradient &Gradient) const
39 {
40  // interface of the base class. Use DeltaGradient for op.
41  std::pair<FunctionGradient, MnAlgebraicVector> mypair = DeltaGradient(par, Gradient);
42 
43  return mypair.first;
44 }
45 
47 {
48  // return the precision
49  return fTransformation.Precision();
50 }
51 
53 {
54  // return number of calculation cycles (defined in strategy)
56 }
57 
59 {
60  // return tolerance on step size (defined in strategy)
62 }
63 
65 {
66  // return gradient tolerance (defines in strategy)
67  return Strategy().GradientTolerance();
68 }
69 
70 std::pair<FunctionGradient, MnAlgebraicVector>
72 {
73  // calculate gradient for Hessian
74  assert(par.IsValid());
75 
76  MnPrint print("HessianGradientCalculator");
77 
78  MnAlgebraicVector x = par.Vec();
79  MnAlgebraicVector grd = Gradient.Grad();
80  const MnAlgebraicVector &g2 = Gradient.G2();
81  // const MnAlgebraicVector& gstep = Gradient.Gstep();
82  // update also gradient step sizes
83  MnAlgebraicVector gstep = Gradient.Gstep();
84 
85  double fcnmin = par.Fval();
86  // std::cout<<"fval: "<<fcnmin<<std::endl;
87 
88  double dfmin = 4. * Precision().Eps2() * (fabs(fcnmin) + Fcn().Up());
89 
90  unsigned int n = x.size();
91  MnAlgebraicVector dgrd(n);
92 
93  MPIProcess mpiproc(n, 0);
94  // initial starting values
95  unsigned int startElementIndex = mpiproc.StartElementIndex();
96  unsigned int endElementIndex = mpiproc.EndElementIndex();
97 
98  for (unsigned int i = startElementIndex; i < endElementIndex; i++) {
99  double xtf = x(i);
100  double dmin = 4. * Precision().Eps2() * (xtf + Precision().Eps2());
101  double epspri = Precision().Eps2() + fabs(grd(i) * Precision().Eps2());
102  double optstp = sqrt(dfmin / (fabs(g2(i)) + epspri));
103  double d = 0.2 * fabs(gstep(i));
104  if (d > optstp)
105  d = optstp;
106  if (d < dmin)
107  d = dmin;
108  double chgold = 10000.;
109  double dgmin = 0.;
110  double grdold = 0.;
111  double grdnew = 0.;
112  for (unsigned int j = 0; j < Ncycle(); j++) {
113  x(i) = xtf + d;
114  double fs1 = Fcn()(x);
115  x(i) = xtf - d;
116  double fs2 = Fcn()(x);
117  x(i) = xtf;
118  // double sag = 0.5*(fs1+fs2-2.*fcnmin);
119  // LM: should I calculate also here second derivatives ???
120 
121  grdold = grd(i);
122  grdnew = (fs1 - fs2) / (2. * d);
123  dgmin = Precision().Eps() * (std::fabs(fs1) + std::fabs(fs2)) / d;
124  // if(fabs(grdnew) < Precision().Eps()) break;
125  if (grdnew == 0)
126  break;
127  double change = fabs((grdold - grdnew) / grdnew);
128  if (change > chgold && j > 1)
129  break;
130  chgold = change;
131  grd(i) = grdnew;
132  // LM : update also the step sizes
133  gstep(i) = d;
134 
135  if (change < 0.05)
136  break;
137  if (fabs(grdold - grdnew) < dgmin)
138  break;
139  if (d < dmin)
140  break;
141  d *= 0.2;
142  }
143 
144  dgrd(i) = std::max(dgmin, std::fabs(grdold - grdnew));
145 
146  print.Debug("HGC Param :", i, "\t new g1 =", grd(i), "gstep =", d, "dgrd =", dgrd(i));
147  }
148 
149  mpiproc.SyncVector(grd);
150  mpiproc.SyncVector(gstep);
151  mpiproc.SyncVector(dgrd);
152 
153  return std::pair<FunctionGradient, MnAlgebraicVector>(FunctionGradient(grd, g2, gstep), dgrd);
154 }
155 
156 } // namespace Minuit2
157 
158 } // namespace ROOT
ROOT::Minuit2::MnMachinePrecision::Eps2
double Eps2() const
eps2 returns 2*sqrt(eps)
Definition: MnMachinePrecision.h:41
n
const Int_t n
Definition: legend1.C:16
ROOT::Minuit2::MPIProcess
Definition: MPIProcess.h:42
MnMachinePrecision.h
ROOT::Minuit2::LAVector
Definition: LAVector.h:32
ROOT::Minuit2::MnPrint::Debug
void Debug(const Ts &... args)
Definition: MnPrint.h:138
MnStrategy.h
MPIProcess.h
ROOT::Minuit2::MnUserTransformation::Precision
const MnMachinePrecision & Precision() const
forwarded interface
Definition: MnUserTransformation.h:117
ROOT::Minuit2::MPIProcess::StartElementIndex
unsigned int StartElementIndex() const
Definition: MPIProcess.h:55
MnUserTransformation.h
ROOT::Minuit2::HessianGradientCalculator::Precision
const MnMachinePrecision & Precision() const
Definition: HessianGradientCalculator.cxx:46
ROOT::Minuit2::MinimumParameters::Vec
const MnAlgebraicVector & Vec() const
Definition: MinimumParameters.h:39
ROOT::Minuit2::MnStrategy::GradientTolerance
double GradientTolerance() const
Definition: MnStrategy.h:42
x
Double_t x[n]
Definition: legend1.C:17
ROOT::Minuit2::HessianGradientCalculator::Strategy
const MnStrategy & Strategy() const
Definition: HessianGradientCalculator.h:50
ROOT::Minuit2::InitialGradientCalculator
Class to calculate an initial estimate of the gradient.
Definition: InitialGradientCalculator.h:27
ROOT::Minuit2::FunctionGradient
Definition: FunctionGradient.h:21
ROOT::Minuit2::MnStrategy::GradientStepTolerance
double GradientStepTolerance() const
Definition: MnStrategy.h:41
ROOT::Minuit2::HessianGradientCalculator::fTransformation
const MnUserTransformation & fTransformation
Definition: HessianGradientCalculator.h:58
MnFcn.h
ROOT::Minuit2::FunctionGradient::G2
const MnAlgebraicVector & G2() const
Definition: FunctionGradient.h:45
ROOT::Minuit2::MinimumParameters::Fval
double Fval() const
Definition: MinimumParameters.h:41
ROOT::Minuit2::MnMachinePrecision::Eps
double Eps() const
eps returns the smallest possible number so that 1.+eps > 1.
Definition: MnMachinePrecision.h:38
ROOT::Minuit2::HessianGradientCalculator::operator()
virtual FunctionGradient operator()(const MinimumParameters &) const
Definition: HessianGradientCalculator.cxx:28
ROOT::Minuit2::HessianGradientCalculator::StepTolerance
double StepTolerance() const
Definition: HessianGradientCalculator.cxx:58
ROOT::Math::fabs
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
Definition: UnaryOperators.h:131
ROOT::Minuit2::MinimumParameters
Definition: MinimumParameters.h:19
FunctionGradient.h
ROOT::Minuit2::MnFcn::Up
double Up() const
Definition: MnFcn.cxx:39
ROOT::Minuit2::MPIProcess::SyncVector
bool SyncVector(ROOT::Minuit2::MnAlgebraicVector &mnvector)
Definition: MPIProcess.cxx:153
ROOT::Minuit2::HessianGradientCalculator::Fcn
const MnFcn & Fcn() const
Definition: HessianGradientCalculator.h:47
ROOT::Minuit2::MnMachinePrecision
Sets the relative floating point (double) arithmetic precision.
Definition: MnMachinePrecision.h:32
ROOT::Minuit2::HessianGradientCalculator::fStrategy
const MnStrategy & fStrategy
Definition: HessianGradientCalculator.h:59
sqrt
double sqrt(double)
ROOT::Minuit2::HessianGradientCalculator::fFcn
const MnFcn & fFcn
Definition: HessianGradientCalculator.h:57
ROOT::Minuit2::FunctionGradient::Grad
const MnAlgebraicVector & Grad() const
Definition: FunctionGradient.h:40
ROOT::Minuit2::MinimumParameters::IsValid
bool IsValid() const
Definition: MinimumParameters.h:42
ROOT::Minuit2::FunctionGradient::Gstep
const MnAlgebraicVector & Gstep() const
Definition: FunctionGradient.h:46
HessianGradientCalculator.h
InitialGradientCalculator.h
ROOT::Minuit2::HessianGradientCalculator::Ncycle
unsigned int Ncycle() const
Definition: HessianGradientCalculator.cxx:52
ROOT::Minuit2::MnStrategy::HessianGradientNCycles
unsigned int HessianGradientNCycles() const
Definition: MnStrategy.h:47
d
#define d(i)
Definition: RSha256.hxx:102
ROOT::Minuit2::HessianGradientCalculator::DeltaGradient
std::pair< FunctionGradient, MnAlgebraicVector > DeltaGradient(const MinimumParameters &, const FunctionGradient &) const
Definition: HessianGradientCalculator.cxx:71
ROOT::Minuit2::MPIProcess::EndElementIndex
unsigned int EndElementIndex() const
Definition: MPIProcess.h:61
MinimumParameters.h
ROOT::Minuit2::HessianGradientCalculator::GradTolerance
double GradTolerance() const
Definition: HessianGradientCalculator.cxx:64
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
ROOT::Minuit2::MnPrint
Definition: MnPrint.h:73