Using external C++ function in the model block

Modeling
1

I have an external C++ function. The header file looks like:

#ifndef tacppnew_hpp
#define tacppnew_hpp

#include <boost/math/distributions/normal.hpp>
#include <stan/math.hpp>
#include <vector>
#include <iostream>
#include <limits>
#include <sstream>



template <typename T_a__, typename T_b__,typename T_theta__>
typename boost::math::tools::promote_args <T_a__, T_b__, T_theta__>::type
test3( const double& a, const double& b, const double beta, const double mu0, const double mu1, const double gamma, const double delta, const double sigma, const double M0, std::ostream& msgs) {
    
     using stan::math::var;
    
   
   
    auto const f2= [&](const double x, const double xc, const std::vector<var> theta,  std::vector<double> x_r, std::vector<int> x_i,  std::ostream& msgs) {
        
        return   ( 0.5 * erfc(-(M0-(theta[1]+theta[2]*exp(- theta[3]* pow((7.0+log(x)),theta[4]))))/(theta[5]*sqrt(2.0))) );
                 
    };
    

    std::vector<double> theta = {beta, mu0, mu1, gamma, delta, sigma};

    var Q = stan::math::integrate_1d(f2, a, b, theta,{},{}, msgs);
     return Q;
    
}
#endif /* tacppnew_hpp */

I want to use this external C++ function in the model block. For this, I declared it in the function block. I can’t understand how to use it in the model block. The current model block looks like:

model{
//priors

K~ exponential(1);
cp~exponential(1);
p~normal(1.12,0.16) T[1,];

//likelihood
real[N] partcA;
for (i==1){
partcA[i] =  partC( t[1], t[2], beta, mu0, mu1, gamma,
delta,  sigma,  M0) 
}
for(i in 2: (N-1)){
partcA[i] += partC( t[i]-t[1:(i-1)], t[i+1]-t[1:(i-1)], beta, mu0, mu1, gamma,
delta,  sigma,  M0) 
}
for (i==N){
partcA[i] += partC( t[i]-t[1:(i-1)], T-t[1:(i-1)], beta, mu0, mu1, gamma,
delta,  sigma,  M0) 
}
real final = partA ( N,   K,  cp,  p, t,    rt,  M0,  T) -  partB- sum(partcA)
}

Am I doing it correctly? Please help.

2

You have to declare (but not define) the function with the same arguments (excluding pstream__) in the functions block of the Stan program. Then call it the same way you would call it if it were defined in the Stan language. There are examples in https://cran.r-project.org/web/packages/rstan/vignettes/external.html .

3

Hi,

thank you. I looked at your provided link. I did not define, I only declared it in the functions block. The example in the provided link does not have anything in the model block.But in my case, I have to do it. I want to know, shall I use it just as in C++ ? Please note, I have tested my external code and it parses successfully.

4

Well, you just call it in the model block like any other function. There was an example of that at

5

Hi,

Thank you for this. I have to do an integration inside my log likelihood function. Therefore I am using the external C++ function. But it seems that whenever I am using that in the “model block”, it shows me error about loading shared object. But if I keep the model block empty, it does not show any error. What is the solution for this?

6

It either means that you haven’t defined one of the specializations that it needs or you are running into the issue that we haven’t figured out where you are calling an externally defined C++ function from within another Stan function, as in the example I linked to.

7

I guess in the example (you gave me) the problem is that the model does not works if it has a function that calls another function. But in my case, I am not doing that. I am directly calling the external C++ function.

8

Then it is probably missing a necessary template specialization

9

The “integrate_1d” in Stan math library is using four templates: F, T_a, T_b, T_theta…

template<typename F , typename T_a , typename T_b , typename T_theta >
var stan::math::integrate_1d	(	const F & 	f,
const T_a & 	a,
const T_b & 	b,
const std::vector< T_theta > & 	theta,
const std::vector< double > & 	x_r,
const std::vector< int > & 	x_i,
std::ostream & 	msgs,
const double 	tolerance = std::sqrt(std::numeric_limits<double>::epsilon()) 
)

I found that I did not use template F. Is it a possible reason for this error?

10

Maybe. But I would try with

auto stan::math::integrate_1d(...);

rather than hard-coding the output type.

11

Could you please explain more? I am using the Stan::math::integrate_1d inside my C++ code which is as below:


template <typename T_a__, typename T_b__,typename T_theta__>
typename boost::math::tools::promote_args <T_a__, T_b__, T_theta__>::type
partC ( const double& a, const double& b, const double beta, const double mu0, const double mu1, const double gamma, const double delta, const double sigma, const double M0, std::ostream& msgs) {
    
     using stan::math::var;
    
   
   
    auto const f2= [&](const double x, const double xc, const std::vector<var> theta,  std::vector<double> x_r, std::vector<int> x_i,  std::ostream& msgs) {
        
        return   ( 0.5 * erfc(-(M0-(theta[1]+theta[2]*exp(- theta[3]* pow((7.0+log(x)),theta[4]))))/(theta[5]*sqrt(2.0))) );
                 
    };
    

    std::vector<double> theta = {beta, mu0, mu1, gamma, delta, sigma};

    var Q = stan::math::integrate_1d(f2, a, b, theta,{},{}, msgs);
     return Q;
    
}
#endif /* tacppnew_hpp */

Are you suggesting to do it that way? Many thanks in advance

12

I was suggesting

auto Q = stan::math::integrate_1d(f2, a, b, theta,{},{}, msgs);
1 Like
13

How you do that?
calling the external C++ Function in rstan: function or model blocks.