Call GNU Scientific Library from stan

Hi all,

I’m writing a function in a stan model in which I want call a intergation precedure in GSL library. Does anyone have experiences with this? Could you please advise me to write this function.
Here is my code
functions {
real integrate(real a, real b, real alpha );
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>

real f (real x, void * params) {
real alpha = *(real ) params;
real f = log(alpha
x) / sqrt(x);
return f;
}

gsl_integration_workspace * w
= gsl_integration_workspace_alloc (1000);
real result;
real error;

gsl_function F;
F.function = &f;
F.params = α
gsl_integration_qng (&F, a, b, 0, 1e-7, 1000,w, &result, &error);
return result;
}

Thanks and best regards,
Nhat

It looks like your are trying to calculate an integral. For that you can abuse the ODE integrators (described on the forum) or you can use the develop version of stan-math where a 1D integrator is defined such that you can use it through the C++ extension facility which is described as part of rstan and cmdstan.

Anything else is probably difficult as you need the gradients of whatever you call (and non-templated C++ code suited for plugging into stan-math is non-trivial).

In general, the Stan Math Library cannot include GNU Scientific Library code because the latter is licensed under the GPL. In particular Stan programs, users are free to utilize GNU Scientific Library code via the mechanism for including external code. But in your case, you are probably better off including the integration functionality from Boost until we can wrap it ourselves in Stan Math.

Hi all,
I found on the github of yixuan ( https://github.com/yixuan/RcppNumerical) a function to compute the integral using Gauss Kronrod method but the code was written in Rcpp.
Here is the code:

// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(RcppNumerical)]]
#include <RcppNumerical.h>
using namespace Numer;
class function: public Func
{
private:
double nu;
public:
function(double nu_) : nu(nu_) {}

double operator()(const double& x) const
{
return exp(x*nu);
}
};

// [[Rcpp::export]]
double integrate_GK(double ll, double ul,double nu)
{
function f(nu);
double err_est;
int err_code;
const double res = integrate(f, ll, ul, err_est, err_code);
return res;
}
/*** R
integrate_GK(0,1,1)
*/

library(rstan)
mc <-

functions {

real integrate_GK(real ll, real ul,real nu);

}
model {

} // use the fib() function somehow

stan_model(model_code = mc, model_name = “external”, allow_undefined = TRUE,
includes = paste0(’\n#include “’,
file.path(getwd(), ‘integrate_GK.cpp’), '”\n’))

I want to use the function integrate_GK in the stan model. But I still dont know how to do it. I try to save the code in .cpp and using the instruction https://cran.r-project.org/web/packages/rstan/vignettes/external.html but I cannot sucess to make the model compile. The error are “Error in .Call(“CPP_stanc280”, model_code, model_cppname, allow_undefined, : “CPP_stanc280” not resolved from current namespace (rstan)”.

My main question is: can we use a defined function in Rcpp and call it a stan model?
Does anyone know how handle this issue? Could you please help me to solve it!

Thank you very much!
Nhat

Your rstan does not seem to be installed properly.

Thanks bgoodri!

You mean that if I install my rstan properly. It is feasible that I can call a Rcpp defined function in a stan model?

Nhat

That was already feasible. You just need to install rstan properly to avoid this error:

Thanks Ben!

I will fix my rstan.

Best,
Nhat

Hi all,

I tried a simple example
mc5 <-

functions {
vector timesTwo(vector x);
}

model_tmp4<-stan_model(model_code = mc5, model_name = “externa4”, allow_undefined = TRUE,
includes = paste0(’\n#include “’,
file.path(getwd(), ‘timesTwo.hpp’), '”\n’))

here is the timesTwo.hpp file

#include <Rcpp.h>
using namespace Rcpp;
NumericVector timesTwo(NumericVector x);
// [[Rcpp::export]]
NumericVector timesTwo(NumericVector x) {
return x * 2;
}

I have error:
Error in compileCode(f, code, language = language, verbose = verbose) :
Compilation ERROR, function(s)/method(s) not created! In file included from C:/Users/Nhat Le Thanh Hoang/Documents/R/win-library/3.3/BH/include/boost/config.hpp:39:0,
from C:/Users/Nhat Le Thanh Hoang/Documents/R/win-library/3.3/BH/include/boost/math/tools/config.hpp:13,
from C:/Users/Nhat Le Thanh Hoang/Documents/R/win-library/3.3/StanHeaders/include/stan/math/rev/core/var.hpp:7,
from C:/Users/Nhat Le Thanh Hoang/Documents/R/win-library/3.3/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
from C:/Users/Nhat Le Thanh Hoang/Documents/R/win-library/3.3/StanHeaders/include/stan/math/rev/core.hpp:12,
from C:/Users/Nhat Le Thanh Hoang/Documents/R/win-library/3.3/StanHeaders/include/stan/math/rev/mat.hpp:4,
from C:/Users/Nhat Le Thanh Hoang/Documents/R/win-library/3.3/StanHeaders/include/stan/math.hpp:4,
from C:/Users/Nhat Le Thanh Hoang/Documents/R/win-
In addition: Warning message:
running command ‘C:/PROGRA~1/R/R-33~1.3/bin/x64/R CMD SHLIB file1f4849c31000.cpp 2> file1f4849c31000.cpp.err.txt’ had status 1

Does anyone know why Rcpp function doesn’t work on the stan model?

Thanks!
Nhat

It is supposed to be

model_tmp4 <- stan_model(model_code = mc5, model_name = “externa4”, allow_undefined = TRUE,
                includes = paste0('\n#include "', file.path(getwd(), "TimesTwo.hpp"), '"\n'))

but I doubt it is going to work with that NumericVector typedef.

Dear Ben,

I tried another example without NumericVector typedef.
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]

double timesTwo(double x) {
return x * 2;
}

I have the same error as before. But if I removed the first three lines
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]

the code worked. I think the main issue is these three lines. Do you have any idea to fix this?

Many thanks!
Nhat

I imagine it is the using namespace Rcpp; that clashes with something. If you have to use Rcpp functions, you should prefix them with Rcpp::.

Hi Ben,

I re-install rtstan and rtools. When I tested Rtools with the code Sys.getenv(‘PATH’)
the output is:
“C:\Program Files\R\R-3.5.0\bin\x64;c:\Rtools\bin;c:\Rtools\mingw_32\bin;c:\Rtools\bin;c:\Rtools\gcc-4.6.3\bin;c:\Rtools\bin;c:\Rtools\bin;c:\Rtools\gcc-4.6.3\bin;c:\Rtools\bin;c:\Rtools\bin;C:\Program Files\Microsoft MPI\Bin\;c:\Rtools\bin;C:\WINDOWS\system32;C:\WINDOWS;C:\WINDOWS\System32\Wbem;C:\WINDOWS\System32\WindowsPowerShell\v1.0\;C:\Program Files\PuTTY\;C:\Strawberry\c\bin;C:\Strawberry\perl\site\bin;C:\Strawberry\perl\bin;C:\Users\Nhat Le Thanh Hoang\AppData\Local\Microsoft\WindowsApps;”

It begins with C:\Program Files\R\R-3.5.0\bin\x64; not with c:\Rtools\bin;c:\Rtools\mingw_32…
Does it mean I didn’t install rtools properly? Hence rstan is not correct as well?

Thanks,
Nhat

I think it is OK for the PATH to begin with the directory to R, but you have some duplicative entries. Not sure if that makes any difference.

Dear Ben,
I tried to make a file.hpp with the library Rcpp. I never success to make the code run. Do you have any example that you included the user-defined function like timesTwo to the stan model?

Many thanks!
Nhat

There are two examples in
https://cran.r-project.org/web/packages/rstan/vignettes/external.html

Hi Ben,
I follow your instruction 1D numerical integration (again) and https://www.boost.org/doc/libs/1_67_0/libs/math/doc/html/quadrature.html to use boost library to write a function in .cpp with Rcpp to compute a integral of exp(x) function.
Here is the code.

#include <Rcpp.h>
#include <boost/math/quadrature/gauss_kronrod.hpp>
using namespace Rcpp;

// [[Rcpp::plugins(cpp14)]]
// [[Rcpp::depends(BH)]]

// [[Rcpp::export]]
double integral_GK(double ll, double ul,double nu) {
boost::math::quadrature::gauss_kronrod<double, 15> integrator;
auto f1 = [&](double x) {
return std::exp(x);
};
double error;
double ll_new=ll * nu;
double ul_new=ul * nu;
double Q = integrator.integrate(f1, ll_new, ul_new, 0, 0, &error)/nu;
return Q;
}

The code works well. What I’m trying is to translate this code to .hpp code as your instruction in https://cran.r-project.org/web/packages/rstan/vignettes/external.html so that I can inlcude in stan model, but I’m not successful to make the code run.

Here is the .hpp code that I tried:
template <typename T0__, typename T1__,typename T2__>
typename boost::math::tools::promote_args<T0__, T1__,T2__>::type
integral_GK(const T0__& ll, const T1__& ul,const T2__& nu, std::ostream* pstream__) {
boost::math::quadrature::gauss_kronrod<double, 15> integrator;
auto f1 = [&](double x) {
return std::exp(x);
};
double error;
double ll_new=ll * nu;
double ul_new=ul * nu;
double Q = integrator.integrate(f1, ll_new, ul_new, 0, 0, &error)/nu;
return Q;
}

and I go to the model_header.hpp and I included the header
#include <boost/math/quadrature/gauss_kronrod.hpp>

Here is a part of the error:
In file included from C:/Users/nhatlth/Documents/R/win-library/3.3/BH/include/boost/config.hpp:39:0,
from C:/Users/nhatlth/Documents/R/win-library/3.3/BH/include/boost/math/tools/config.hpp:13,
from C:/Users/nhatlth/Documents/R/win-library/3.3/StanHeaders/include/stan/math/rev/core/var.hpp:7,
from C:/Users/nhatlth/Documents/R/win-library/3.3/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
from C:/Users/nhatlth/Documents/R/win-library/3.3/StanHeaders/include/stan/math/rev/core.hpp:12,
from C:/Users/nhatlth/Documents/R/win-library/3.3/StanHeaders/include/stan/math/rev/mat.hpp:4,
from C:/Users/nhatlth/Documents/R/win-library/3.3/StanHeaders/include/stan/math.hpp:4,
from C:/Users/nhatlth/Documents/R/win-library/3.3/StanHeaders/include/src/stan/model/model_header.hpp:4,
from file29e417e718b4.cpp:8:
C:/Users/nhatlth/Documents/R/win-library/3.3/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: “BOOST_NO_CXX11_RVALUE_REFERENCES” redefined

define BOOST_NO_CXX11_RVALUE_REFERENCES

^
:0:0: note: this is the location of the previous definition
In file included from C:/Users/nhatlth/Documents/R/win-library/3.3/BH/include/boost/math/quadrature/gauss_kronrod.hpp:18:0,
from C:/Users/nhatlth/Documents/R/win-library/3.3/StanHeaders/include/src/stan/model/model_header.hpp:20,
from file29e417e718b4.cpp:8:
C:/Users/nhatlth/Documents/R/win-library/3.3/BH/include/boost/math/special_functions/legendre_stieltjes.hpp: In member function ‘std::vector boost::math::legendre_stieltjes::zeros() const’:
C:/Users/nhatlth/Documents/R/win-library/3.3/BH/include/boost/math/special_functions/legendre_stieltjes.hpp:205:18: error: ‘g’ does not name a type
auto g = [&](Real t) { return this->operator()(t); };
^
C:/Users/nhatlth/Documents/R/win-library/3.3/BH/include/boost/math/special_functions/legendre_stieltjes.hpp:206:18: error: ‘p’ does not name a type
auto p = boost::math::tools::bisect(g, lower_bound, upper_bound, tol);
^
C:/Users/nhatlth/Documents/R/win-library/3.3/BH/include/boost/math/special_functions/legendre_stieltjes.hpp:208:31: error: ‘p’ was not declared in this scope
Real x_nk_guess = p.first + (p.second - p.first)*half();
^
C:/Users/nhatlth/Documents/R/win-library/3.3/BH/include/boost/math/special_functions/legendre_stieltjes.hpp:211:18: error: ‘f’ does not name a type
auto f = [&] (Real x) { Real Pn = this->operator()(x);
^…

Can you please help me to translate this code into a code.hpp so that I can run in stan model. I know that I ask too much from you. But it seems that I’m stuck here since my C++ knowledge is very limited.

Thank you very much in advance!
Nhat

Judging from error messages like this:

error: ‘g’ does not name a type

where g was declared with auto, I suspect that part of the problem is that you are trying to compile C++11 code with a compiler that is expecting older C++98 syntax (which, AFAIK, uses auto in a very different and mostly useless way that’s a holdover from legacy C).

What compiler are you using, and what flags are you passing to it?

Hi James,

What I used is rstan and the code was in C++. I checked that in my computer there are many the C++ complier: MV C++ 2005 restributable to MV C++ 2015 Restributable and DEV C++ as well. I don’t know what flag means.
The thing is if I run the code in Rcpp format, it works well. So I think, the computer understands the syntax and know all the library. But if I include the code in the stan model then for some reason they cannot link to the correct library.

In the intsruction https://cran.r-project.org/web/packages/rstan/vignettes/external.html. They say that we can call a function in Boost Library Math whose header are pulled by the Stan Math Library by using the prefix boost::math.

The same thing that I want to use is to call the class gauss_kronrod<double, 15> from the boost::math::quadrature library by defining a new class integrator as boost::math::quadrature::gauss_kronrod<double, 15> integrator;

I think at this step the stan don’t know where the library boost::math::quadrature comes from. My question is what should I do to be able to add the line using namespace boost::math::quadrature;
in the program source:
"Program source:
1:
2: // includes from the plugin
3:
4:
5: // user includes
6: #define STAN__SERVICES__COMMAND_HPP// Code generated by Stan version 2.17.0
7:
8: #include <stan/model/model_header.hpp>
9:
10: namespace model30c83e0371b6_external_namespace {
11:
12: using std::istream;
13: using std::string;
14: using std::stringstream;
15: using std::vector;
16: using stan::io::dump;
17: using stan::math::lgamma;
18: using stan::model::prob_grad;
19: using namespace stan::math;
20:
21: typedef Eigen::Matrix<double,Eigen::Dynamic,1> vector_d;
22: typedef Eigen::Matrix<double,1,Eigen::Dynamic> row_vector_d;
23: typedef Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> matrix_d;
24:
25: static int current_statement_begin__;
26:
27: stan::io::program_reader prog_reader__() {
28: stan::io::program_reader reader;
29: reader.add_event(0, 0, “start”, “model30c83e0371b6_external”);
30: reader.add_event(9, 9, “end”, “model30c83e0371b6_external”);
31: return reader;
32: }
33:
34: template <typename T0__, typename T1__, typename T2__>
35: typename boost::math::tools::promote_args<T0__, T1__, T2__>::type
36: integral_GK(const T0__& ll,
37: const T1__& ul,
38: const T2__& nu, std::ostream* pstream__);
39:
40:
41: #include “C:/Users/nhatlth/Documents/exp_GK.hpp”
42: class model30c83e0371b6_external : public prob_grad {
43: private:
44: public:
45: model30c83e0371b6_external(stan::io::var_context& context__,
46: std::ostream* pstream__ = 0)
47: : prob_grad(0) {
48: ctor_body(context__, 0, pstream__);
49: }
50:
51: model30c83e0371b6_external(stan::io::var_context& context__,
52: unsigned int random_seed__,
53: std::ostream* pstream__ = 0)
54: : prob_grad(0) {
55: ctor_body(context__, random_seed__, pstream__);
56: }

"

Thanks,
Nhat

A flag is a command-line option. For example, if you were compiling C++11 code with an older g++ compiler, you might need to compile it like this:

g++ -std=c++11 -c my_code.cpp

not like this:

g++ -c my_code.cpp

Here, -std=c++11 is a flag.

What that probably really means is that the code that you used with Rcpp was compiled with a different compiler than the one that is giving you the error.

In any case, you really should be showing us exactly what steps you took that led to the error. Chances are the error is due to some step that you haven’t yet mentioned.