Public Member Functions | Protected Member Functions | Protected Attributes | List of all members
ana::GradientDescent Class Reference

#include "/cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N20-11-30/CAFAna/Fit/GradientDescent.h"

Inheritance diagram for ana::GradientDescent:

Public Member Functions

 GradientDescent (const ROOT::Minuit2::FCNGradientBase &func, const ROOT::Minuit2::MnUserParameters &pars)
 
virtual ROOT::Minuit2::FunctionMinimum operator() (unsigned int maxfcn, double tolerance) override
 
virtual const ROOT::Minuit2::ModularFunctionMinimizer & Minimizer () const override
 
virtual ROOT::Minuit2::ModularFunctionMinimizer & Minimizer () override
 

Protected Member Functions

ROOT::Minuit2::FunctionMinimum Package (const std::vector< double > &pt, double chi, int ncalls) const
 
double Magnitude (const std::vector< double > &xs) const
 
void MakeUnit (std::vector< double > &xs) const
 

Protected Attributes

const ROOT::Minuit2::FCNGradientBase & fFunc
 
const ROOT::Minuit2::MnUserParameters & fPars
 

Detailed Description

A minimalistic gradient descent fitter to complement MINUIT's more elaborate offerings

Definition at line 9 of file GradientDescent.h.

Constructor & Destructor Documentation

ana::GradientDescent::GradientDescent ( const ROOT::Minuit2::FCNGradientBase &  func,
const ROOT::Minuit2::MnUserParameters &  pars 
)

Definition at line 17 of file GradientDescent.cxx.

References operator()().

19  : ROOT::Minuit2::MnApplication(func, gState, gStrat),
20  fFunc(func), fPars(pars)
21  {
22  }
const ROOT::Minuit2::MnUserParameters & fPars
ROOT::Minuit2::MnUserParameterState gState
ROOT::Minuit2::MnStrategy gStrat
double func(double x, double y)
const ROOT::Minuit2::FCNGradientBase & fFunc
std::string pars("Th23Dmsq32")

Member Function Documentation

double ana::GradientDescent::Magnitude ( const std::vector< double > &  xs) const
protected

Definition at line 130 of file GradientDescent.cxx.

References runNovaSAM::ret, std::sqrt(), and submit_syst::x.

Referenced by MakeUnit(), Minimizer(), and operator()().

131  {
132  double ret = 0;
133  for(double x: xs) ret += x*x;
134  return sqrt(ret);
135  }
T sqrt(T number)
Definition: d0nt_math.hpp:156
void ana::GradientDescent::MakeUnit ( std::vector< double > &  xs) const
protected

Definition at line 138 of file GradientDescent.cxx.

References Magnitude(), and submit_syst::x.

Referenced by Minimizer(), and operator()().

139  {
140  const double mag = Magnitude(xs);
141  for(double& x: xs) x /= mag;
142  }
double Magnitude(const std::vector< double > &xs) const
virtual const ROOT::Minuit2::ModularFunctionMinimizer& ana::GradientDescent::Minimizer ( ) const
inlineoverridevirtual

Definition at line 18 of file GradientDescent.h.

19  {
20  abort();
21  }
virtual ROOT::Minuit2::ModularFunctionMinimizer& ana::GradientDescent::Minimizer ( )
inlineoverridevirtual

Definition at line 23 of file GradientDescent.h.

References Magnitude(), MakeUnit(), Package(), and gen_flatrecord::pt.

24  {
25  abort();
26  }
ROOT::Minuit2::FunctionMinimum ana::GradientDescent::operator() ( unsigned int  maxfcn,
double  tolerance 
)
overridevirtual

Definition at line 26 of file GradientDescent.cxx.

References e, fFunc, fPars, stan::math::grad(), MECModelEnuComparisons::i, Magnitude(), MakeUnit(), Package(), gen_flatrecord::pt, and fillBadChanDBTables::step.

Referenced by GradientDescent().

27  {
28  // Initialize at the input parameters
29  std::vector<double> pt = fPars.Params();
30  const unsigned int N = pt.size();
31 
32  // Use the user errors as a crude way to set the initial step size
33  double step = Magnitude(fPars.Errors());
34 
35  // Usually we reuse the chisq from the last iteration, so need to
36  // initialize it once out here.
37  double chi = fFunc(pt);
38  int ncalls = 1;
39 
40  while(true){
41  // std::cout << "New direction, chisq = " << chi << std::endl;
42 
43  // Evaluate gradient to figure out a direction to step in
44  std::vector<double> grad = fFunc.Gradient(pt);
45  ++ncalls;
46  // Make gradient into a unit vector but keep the magnitude around
47  const double gradmag = Magnitude(grad);
48  MakeUnit(grad);
49 
50  // Take step in that direction
51  std::vector<double> trialpt = pt;
52  for(unsigned int i = 0; i < N; ++i) trialpt[i] -= grad[i]*step;
53 
54  // And evaluate the chisq there
55  const double trialchi = fFunc(trialpt);
56  ++ncalls;
57 
58  // std::cout << "step = " << step << " -> chisq " << trialchi << std::endl;
59 
60  // Estimate the second derivative from the two points and one
61  // gradient
62  const double d2 = (trialchi+gradmag*step-chi)/(step*step);
63  // We expect the function to be convex
64  if(d2 > 0){
65  // std::cout << " d2 = " << d2;
66  // If so, we can compute a better step size, one that would place us
67  // right at the minimum if the function only had quadratic terms.
68  step = gradmag/(2*d2);
69  // std::cout << " new step = " << step << std::endl;
70  }
71 
72  while(true){
73  // Keep trying steps until we find one that reduces the chisq
74  std::vector<double> newpt = pt;
75  for(unsigned int i = 0; i < N; ++i) newpt[i] -= grad[i]*step;
76 
77  const double newchi = fFunc(newpt);
78  ++ncalls;
79  // std::cout << " step = " << step << " -> chisq = " << newchi << std::endl;
80 
81  if(newchi > chi){
82  // If the chisq went up instead, try again with smaller step
83  step /= 2;
84  if(step < 1e-8){
85  // std::cout << "Step too small!" << std::endl;
86  // Maybe it's just because we're already in the minimum. If we
87  // can't improve even with very tiny steps, call it a day.
88  return Package(pt, chi, ncalls);
89  }
90  continue;
91  }
92  else{
93  if(chi-newchi < 1e-5){
94  // std::cout << "Small chisq step" << std::endl;
95  // If the improvement we did get is really tiny we're also likely
96  // already at the minimum
97  return Package(newpt, newchi, ncalls);
98  }
99 
100  // In all other cases (ie we took a reasonable step and found a
101  // reasonably better chisq) we want to update our state, take another
102  // look at the gradient vector to figure out which direction to go
103  // next, and preserve our step size, which is some kind of good
104  // estimate of the scale of the function.
105  pt = newpt;
106  chi = newchi;
107  break;
108  }
109  } // end line search
110  } // end while (searching for better pt)
111  }
const ROOT::Minuit2::MnUserParameters & fPars
double Magnitude(const std::vector< double > &xs) const
static void grad(vari *vi)
Definition: grad.hpp:30
ROOT::Minuit2::FunctionMinimum Package(const std::vector< double > &pt, double chi, int ncalls) const
const ROOT::Minuit2::FCNGradientBase & fFunc
void MakeUnit(std::vector< double > &xs) const
Float_t e
Definition: plot.C:35
ROOT::Minuit2::FunctionMinimum ana::GradientDescent::Package ( const std::vector< double > &  pt,
double  chi,
int  ncalls 
) const
protected

Definition at line 115 of file GradientDescent.cxx.

References MECModelEnuComparisons::i, seed, and makeDatasetsPage::state.

Referenced by Minimizer(), and operator()().

116  {
117  // In practice we only use the chisq and the parameter values, so don't be
118  // too careful about setting this all up.
119  const unsigned int N = pt.size();
120  ROOT::Minuit2::MnAlgebraicVector vec(N);
121  for(unsigned int i = 0; i < N; ++i) vec(i) = pt[i];
122  ROOT::Minuit2::MinimumParameters params(vec, chi);
123  ROOT::Minuit2::MinimumState state(params, 0, ncalls);
124  ROOT::Minuit2::MnUserTransformation trans(pt, std::vector<double>(N));
125  ROOT::Minuit2::MinimumSeed seed(state, trans);
126  return ROOT::Minuit2::FunctionMinimum(seed, {state}, chi);
127  }
unsigned int seed
Definition: runWimpSim.h:102
Eigen::VectorXd vec

Member Data Documentation

const ROOT::Minuit2::FCNGradientBase& ana::GradientDescent::fFunc
protected

Definition at line 35 of file GradientDescent.h.

Referenced by operator()().

const ROOT::Minuit2::MnUserParameters& ana::GradientDescent::fPars
protected

Definition at line 36 of file GradientDescent.h.

Referenced by operator()().


The documentation for this class was generated from the following files: