forked from root-project/root
-
Notifications
You must be signed in to change notification settings - Fork 0
/
TRFunctionExport.h
185 lines (149 loc) · 6.4 KB
/
TRFunctionExport.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
// @(#)root/r:$Id$
// Author: Omar Zapata 16/06/2013
/*************************************************************************
* Copyright (C) 2013-2015, Omar Andres Zapata Mesa *
* All rights reserved. *
* *
* For the licensing terms see $ROOTSYS/LICENSE. *
* For the list of contributors see $ROOTSYS/README/CREDITS. *
*************************************************************************/
#ifndef ROOT_R_TRFunctionExport
#define ROOT_R_TRFunctionExport
#ifndef ROOT_R_TRInternalFunction
#include<TRInternalFunction.h>
#endif
namespace ROOT {
namespace R {
/**
\class TRFunctionExport
This is a class to pass functions from ROOT to R
<center><h2>TRFunctionExport class</h2></center>
<p>
The TRFunctionExport class lets you pass ROOT's functions to R's environment<br>
</p>
<p>
The next example was based in <br>
<a href="http://root.cern.ch/root/html/tutorials/fit/NumericalMinimization.C.html">
http://root.cern.ch/root/html/tutorials/fit/NumericalMinimization.C.html
</a><br>
<a href="http://stat.ethz.ch/R-manual/R-devel/library/stats/html/optim.html">
http://stat.ethz.ch/R-manual/R-devel/library/stats/html/optim.html</a><br>
</p>
Let \f$ f(x,y)=(x-1)^{2} + 100(y-x^{2})^{2} \f$ , which is called the Rosenbrock
function.
It's known that this function has a minimum when \f$ y = x^{2}\f$ , and \f$ x = 1.\f$
Let's get the minimum using R's optim package through ROOTR's interface.
In the code this function was called "Double_t RosenBrock(const TVectorD xx )", because for
optim, the input in your function definition must be a single vector.
The Gradient is formed by
\f$ \frac{\partial f}{\partial x} = -400x(y - x^{2}) - 2(1 - x) \f$
\f$ \frac{\partial f}{\partial y} = 200(y - x^{2}); \f$
The "TVectorD RosenBrockGrad(const TVectorD xx )" function
must have a single vector as the argument a it will return a single vetor.
\code{.cpp}
#include<TRInterface.h>
//in the next function the pointer *double must be changed by TVectorD, because the pointer has no
//sense in R's environment.
Double_t RosenBrock(const TVectorD xx )
{
const Double_t x = xx[0];
const Double_t y = xx[1];
const Double_t tmp1 = y-x*x;
const Double_t tmp2 = 1-x;
return 100*tmp1*tmp1+tmp2*tmp2;
}
TVectorD RosenBrockGrad(const TVectorD xx )
{
const Double_t x = xx[0];
const Double_t y = xx[1];
TVectorD grad(2);
grad[0]=-400 * x * (y - x * x) - 2 * (1 - x);
grad[1]=200 * (y - x * x);
return grad;
}
void Minimization()
{
ROOT::R::TRInterface &r=ROOT::R::TRInterface::Instance();
//passing RosenBrock function to R
r["RosenBrock"]<<ROOT::R::TRFunctionExport(RosenBrock);
//passing RosenBrockGrad function to R
r["RosenBrockGrad"]<<ROOT::R::TRFunctionExport(RosenBrockGrad);
//the option "method" could be "Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN","Brent"
//the option "control" lets you put some constraints like:
//"maxit" The maximum number of iterations
//"abstol" The absolute convergence tolerance.
//"reltol" Relative convergence tolerance.
r<<"result <- optim( c(0.01,0.01), RosenBrock,method='BFGS',control = list(maxit = 1000000) )";
//Getting results from R
TVectorD min=r.Eval("result$par");
std::cout.precision(8);
//printing results
std::cout<<"-----------------------------------------"<<std::endl;
std::cout<<"Minimum x="<<min[0]<<" y="<<min[1]<<std::endl;
std::cout<<"Value at minimum ="<<RosenBrock(min)<<std::endl;
//using the gradient
r<<"optimHess(result$par, RosenBrock, RosenBrockGrad)";
r<<"hresult <- optim(c(-1.2,1), RosenBrock, NULL, method = 'BFGS', hessian = TRUE)";
//getting the minimum calculated with the gradient
TVectorD hmin=r.Eval("hresult$par");
//printing results
std::cout<<"-----------------------------------------"<<std::endl;
std::cout<<"Minimization with the Gradient"<<endl;
std::cout<<"Minimum x="<<hmin[0]<<" y="<<hmin[1]<<std::endl;
std::cout<<"Value at minimum ="<<RosenBrock(hmin)<<std::endl;
}
\endcode
Output
\code
Processing Minimization.C...
-----------------------------------------
Minimum x=0.99980006 y=0.99960016
Value at minimum =3.9974288e-08
-----------------------------------------
Minimization with the Gradient
Minimum x=0.99980443 y=0.99960838
Value at minimum =3.8273828e-08
\endcode
<h2>Users Guide </h2>
<a href="http://oproject.org/tiki-index.php?page=ROOT+R+Users+Guide"> http://oproject.org/tiki-index.php?page=ROOT+R+Users+Guide</a><br>
<a href="https://root.cern.ch/drupal/content/how-use-r-root-root-r-interface"> https://root.cern.ch/drupal/content/how-use-r-root-root-r-interface</a>
@ingroup R
*/
class TRInterface;
class TRFunctionExport: public TObject {
friend class TRInterface;
friend SEXP Rcpp::wrap<TRFunctionExport>(const TRFunctionExport &f);
protected:
TRInternalFunction *f; //Internar Function to export
public:
/**
Default TRFunctionExport constructor
*/
TRFunctionExport();
/**
TRFunctionExport copy constructor
\param fun other TRFunctionExport
*/
TRFunctionExport(const TRFunctionExport &fun);
/**
TRFunctionExport template constructor that supports a lot of function's prototypes
\param fun supported function to be wrapped by Rcpp
*/
template<class T> TRFunctionExport(T fun)
{
f = new TRInternalFunction(fun);
}
/**
function to assign function to export,
template method that supports a lot of function's prototypes
\param fun supported function to be wrapped by Rcpp
*/
template<class T> void SetFunction(T fun)
{
f = new TRInternalFunction(fun);
}
ClassDef(TRFunctionExport, 0) //
};
}
}
#endif