Logo ROOT   6.10/09
Reference Guide
RMinimizer.cxx
Go to the documentation of this file.
1 
2 #include "TRInterface.h"
3 #include "Math/RMinimizer.h"
4 #include "Math/IFunction.h"
5 #include <TVectorD.h>
6 #include "Math/BasicMinimizer.h"
7 
8 namespace ROOT {
9  namespace Math{
10 
11  /// function wrapper for the function to be minimized
13  /// function wrapper for the gradient of the function to be minimized
15  /// integer for the number of function calls
16  int gNCalls = 0;
17 
18  ///function to return the function values at point x
19  double minfunction(const std::vector<double> & x){
20  gNCalls++;
21  //return (*gFunction)(x.GetMatrixArray());
22  return (*gFunction)(x.data());
23  }
24  ///function to return the gradient values at point y
26  unsigned int size = y.GetNoElements();
27  const double * yy = y.GetMatrixArray();
28  double z[size];
29  gGradFunction->Gradient(yy,z);
30  TVectorD zz(size,z);
31  return zz;
32  }
33 
34  /*Default constructor with option for the method of minimization, can be any of the following:
35  *
36  *"Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN", "Brent" (Brent only for 1D minimization)
37  */
39  fMethod=method;
40  if (fMethod.empty() || fMethod=="Migrad") fMethod="BFGS";
41  }
42 
43  ///returns number of function calls
44  unsigned int RMinimizer::NCalls() const { return gNCalls; }
45 
46  ///function for finding the minimum
48 
49  //Set the functions
52 
53  gNCalls = 0;
54 
55  //pass functions and variables to R
57 
58  r["minfunction"] = ROOT::R::TRFunctionExport(minfunction);
59  r["mingradfunction"] = ROOT::R::TRFunctionExport(mingradfunction);
60  r["method"] = fMethod.c_str();
61  std::vector<double> stepSizes(StepSizes(), StepSizes()+NDim());
62  std::vector<double> values(X(), X()+NDim());
63  r["ndim"] = NDim();
64  int ndim = NDim();
65  r["stepsizes"] = stepSizes;
66  r["initialparams"] = values;
67 
68  //check if optimx is available
69  bool optimxloaded = FALSE;
70  r["optimxloaded"] = optimxloaded;
71  r.Execute("optimxloaded<-library(optimx,logical.return=TRUE)");
72  //int ibool = r.ParseEval("optimxloaded").ToScalar<Int_t>();
73  int ibool = r.Eval("optimxloaded");
74  if (ibool==1) optimxloaded=kTRUE;
75 
76  //string for the command to be processed in R
77  TString cmd;
78 
79  //optimx is available and loaded
80  if (optimxloaded==kTRUE) {
81  if (!gGradFunction) {
82  // not using gradient function
83  cmd = TString::Format("result <- optimx( initialparams, minfunction,method='%s',control = list(ndeps=stepsizes,maxit=%d,trace=%d,abstol=%e),hessian=TRUE)",fMethod.c_str(),MaxIterations(),PrintLevel(),Tolerance());
84  }
85  else {
86  // using user provided gradient
87  cmd = TString::Format("result <- optimx( initialparams, minfunction,mingradfunction, method='%s', control = list(ndeps=stepsizes,maxit=%d,trace=%d,abstol=%e),hessian=TRUE)",fMethod.c_str(),MaxIterations(),PrintLevel(),Tolerance());
88 
89  }
90  }
91 
92  //optimx is not available
93  else {
94  if (!gGradFunction) {
95  // not using gradient function
96  cmd = TString::Format("result <- optim( initialparams, minfunction,method='%s',control = list(ndeps=stepsizes,maxit=%d,trace=%d,abstol=%e),hessian=TRUE)",fMethod.c_str(),MaxIterations(),PrintLevel(),Tolerance());
97  }
98  else {
99  // using user provided gradient
100  cmd = TString::Format("result <- optim( initialparams, minfunction,mingradfunction, method='%s', control = list(ndeps=stepsizes,maxit=%d,trace=%d,abstol=%e),hessian=TRUE)",fMethod.c_str(),MaxIterations(),PrintLevel(),Tolerance());
101  }
102  }
103  //execute the minimization in R
104  std::cout << "Calling R with command " << cmd << std::endl;
105  r.Execute(cmd.Data());
106 
107  //results with optimx
108  if (optimxloaded){
109  //get result from R
110  r.Execute("par<-coef(result)");
111  //get hessian matrix (in list form)
112  r.Execute("hess<-attr(result,\"details\")[,\"nhatend\"]");
113  //convert hess to a matrix
114  r.Execute("hess<-sapply(hess,function(x) x)");
115  //convert to square matrix
116  r.Execute("hess<-matrix(hess,c(ndim,ndim))");
117  //find covariant matrix from inverse of hess
118  r.Execute("cov<-solve(hess)");
119  //get errors from the sqrt of the diagonal of cov
120  r.Execute("errors<-sqrt(abs(diag(cov)))");
121  }
122 
123  //results with optim
124  else {
125  r.Execute("par<-result$par");
126  r.Execute("hess<-result$hessian");
127  r.Execute("cov<-solve(hess)");
128  r.Execute("errors<-sqrt(abs(diag(cov)))");
129  }
130 
131  //return the minimum to ROOT
132  //TVectorD vector = gR->ParseEval("par").ToVector<Double_t>();
133  std::vector<double> vectorPar = r["par"];
134 
135  //get errors and matrices from R
136  // ROOT::R::TRObjectProxy p = gR->ParseEval("cov");
137  // TMatrixD cm = p.ToMatrix<Double_t>();
138  TMatrixD cm = r["cov"];
139  // p = gR->ParseEval("errors");
140  // TVectorD err = p.ToVector<Double_t>();
141  std::vector<double> err = r["errors"];
142  // p = gR->ParseEval("hess");
143  // TMatrixD hm = p.ToMatrix<Double_t>();
144  TMatrixD hm = r["hess"];
145 
146  //set covariant and Hessian matrices and error vector
147  fCovMatrix.ResizeTo(ndim,ndim);
148  fHessMatrix.ResizeTo(ndim,ndim);
149  //fErrors.ResizeTo(ndim);
150  fCovMatrix = cm;
151  fErrors = err;
152  fHessMatrix = hm;
153 
154  //get values and show minimum
155  const double *min=vectorPar.data();
156  SetFinalValues(min);
157  SetMinValue((*gFunction)(min));
158  std::cout<<"Value at minimum ="<<MinValue()<<std::endl;
159 
160  return kTRUE;
161  }
162 #ifdef LATER
163  //Returns the ith jth component of the covarient matrix
164  double RMinimizer::CovMatrix(unsigned int i, unsigned int j) const {
165  unsigned int ndim = NDim();
166  if (fCovMatrix==0) return 0;
167  if (i > ndim || j > ndim) return 0;
168  return fCovMatrix[i][j];
169  }
170  // //Returns the full parameter error vector
171  // TVectorD RMinimizer::RErrors() const {
172  // return fErrors;
173  // }
174  //Returns the ith jth component of the Hessian matrix
175  double RMinimizer::HessMatrix(unsigned int i, unsigned int j) const {
176  unsigned int ndim = NDim();
177  if (fHessMatrix==0) return 0;
178  if (i > ndim || j > ndim) return 0;
179  return fHessMatrix[i][j];
180  }
181 #endif
182  } // end namespace MATH
183 }
Interface (abstract class) for multi-dimensional functions providing a gradient calculation.
Definition: IFunction.h:330
const ROOT::Math::IMultiGradFunction * gGradFunction
function wrapper for the gradient of the function to be minimized
Definition: RMinimizer.cxx:14
virtual const double * X() const
return pointer to X values at the minimum
void Execute(const TString &code)
Method to eval R code.
Definition: TRInterface.cxx:82
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
RMinimizer(Option_t *method)
Default constructor.
Definition: RMinimizer.cxx:38
const char Option_t
Definition: RtypesCore.h:62
TMatrixD fCovMatrix
covariant matrix
Definition: RMinimizer.h:37
std::string fMethod
minimizer method to be used, must be of a type listed in R optim or optimx descriptions ...
Definition: RMinimizer.h:33
Basic string class.
Definition: TString.h:129
std::vector< double > fErrors
vector of parameter errors
Definition: RMinimizer.h:36
virtual void Gradient(const double *x, double *grad) const
Evaluate all the vector of function derivatives (gradient) at a point x.
Definition: IFunction.h:350
virtual double MinValue() const
return minimum function value
virtual TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t=-1)
Set size of the matrix to nrows x ncols New dynamic elements are created, the overlapping part of the...
Definition: TMatrixT.cxx:1210
int PrintLevel() const
minimizer configuration parameters
Definition: Minimizer.h:411
virtual bool Minimize()
Function to find the minimum.
Definition: RMinimizer.cxx:47
virtual const double * StepSizes() const
accessor methods
virtual unsigned int NCalls() const
Returns the number of function calls.
Definition: RMinimizer.cxx:44
int gNCalls
integer for the number of function calls
Definition: RMinimizer.cxx:16
double HessMatrix(unsigned int i, unsigned int j) const
Returns the ith jth component of the Hessian matrix.
Double_t x[n]
Definition: legend1.C:17
static TString Format(const char *fmt,...)
Static method which formats a string using a printf style format descriptor and return a TString...
Definition: TString.cxx:2345
void SetMinValue(double val)
unsigned int MaxIterations() const
max iterations
Definition: Minimizer.h:417
virtual unsigned int NDim() const
number of dimensions
const ROOT::Math::IMultiGenFunction * ObjFunction() const
return pointer to used objective function
double Tolerance() const
absolute tolerance
Definition: Minimizer.h:420
virtual double CovMatrix(unsigned int ivar, unsigned int jvar) const
return covariance matrices element for variables ivar,jvar if the variable is fixed the return value ...
Definition: RMinimizer.h:66
Element * GetMatrixArray()
Definition: TVectorT.h:78
Documentation for the abstract class IBaseFunctionMultiDim.
Definition: IFunction.h:62
const ROOT::Math::IMultiGenFunction * gFunction
function wrapper for the function to be minimized
Definition: RMinimizer.cxx:12
void SetFinalValues(const double *x)
TRandom2 r(17)
Int_t GetNoElements() const
Definition: TVectorT.h:76
TMatrixD fHessMatrix
Hessian matrix.
Definition: RMinimizer.h:38
const ROOT::Math::IMultiGradFunction * GradObjFunction() const
return pointer to used gradient object function (NULL if gradient is not supported) ...
Double_t y[n]
Definition: legend1.C:17
static TRInterface & Instance()
static method to get an TRInterface instance reference
TVectorD mingradfunction(TVectorD y)
function to return the gradient values at point y
Definition: RMinimizer.cxx:25
Namespace for new Math classes and functions.
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
#define FALSE
double minfunction(const std::vector< double > &x)
function to return the function values at point x
Definition: RMinimizer.cxx:19
Int_t Eval(const TString &code, TRObject &ans)
Method to eval R code and you get the result in a reference to TRObject.
Definition: TRInterface.cxx:63
const Bool_t kTRUE
Definition: RtypesCore.h:91
This is a class to pass functions from ROOT to R
const char * Data() const
Definition: TString.h:347