Include a C++ function in a Stan code in Julia

Algorithms
1

Dear all,

I am new to Stan and Julia.

I need to include a C++ function in a stan model but I don’t know how I can do it.

The function is remainder and then I include the associated library but I get an error:

SystemError: opening file "<cmath>": No such file or directory

Here a draft:

stancode = "functions {
   real[] P_growth(real t,       // time
                   real[] x,      // state
                   real[] theta, // parameters
                   real[] x_r,   // environmental data
                   int[] x_i){
    real gamma  = theta[1];
    real lambda = theta[2];

    real growth = gamma*x[1];

    #include <cmath>
    if (remainder(t,24) < 14) {
        real loss   = lambda*x[1]*x[1];
    } else {
        real loss   = 0;
    }
    
    
    return {growth - loss};
  }
}
data {
    int<lower = 0> N;           // num obs
    real t_obs[N];              // obs times
    real<lower = 0> y[N];       // observed variable at measurement times
    real sigma;
}
parameters {
    real<lower=0,upper=1> theta[2];      // parameters
    real<lower=0> x0[1];
}
transformed parameters {
    real x[N,1] = integrate_ode_rk45(P_growth, x0, -1, t_obs, theta,
                                      rep_array(0.0, 0), rep_array(0, 0),
                                      1e-6, 1e-5, 1e3);
}
model {
    //theta[1] ~ normal(0.1, 2);
    //theta[2] ~ normal(0.1, 2);
    x0       ~ normal(1.0, 10);
    y[1:N]   ~ normal(x[1:N,1], sigma); // obs
}";

I tried a bit regarding this link: https://cran.r-project.org/web/packages/rstan/vignettes/external.html

and created a function for the remainder:

// define the remainder
functions {
            int remainder(real n);
            int remainder(real n) {
            return remainder(n,24);
    }
}

But I don’t really know where to place it in the Stan prog. and I tried several attempts that didn’t work.

Thanks a lot for your time and help!

2

To the best of my knowledge, that is a feature that is unique to Rstan.

3

(And PyStan).

I think you can manually edit the .hpp file created by the stanc. And then add needed commandline flags for include folders.

1 Like
4

Thanks a lot for the info. Do you mean the file in the cmdstan folder?

But there is not an easy easy way to call classical maths functions in the stancode? Maybe I completely misunderstand my initial error:

A returning function was expected but an undeclared identifier 'remainder' was supplied.

Do I need to declare remainder before to use it?

5

The shortest version of the answer is that you can’t use your own C++ code in Stan via the Julia interface. That feature is only available in Rstan and Pystan. In that case you need to save the C++ file as a separate file (eg remainder.hpp), put the signature in the functions block, and call Rstan as


mc <- 
'
functions { real remainder(real x); }
transformed data { real remaindee_pi = remainder(pi()); }
'
stan_model(model_code = mc, model_name = "external", allow_undefined = TRUE,
           includes = paste0('\n#include "', 
                             file.path(getwd(), 'remainder.hpp'), '"\n'))

If you are very good with C++ and you know how the C++ Stan library works, you can do what Ari suggested and manually edit the C++ code that Stan produces. But this is a really advanced move.

Tl;dr: You can’t do what you want to do through the Julia interface.

Edit: copy and paste error

1 Like
6

Thanks a lot Daniel. It seems a bit too much where I feel there is clearly easiest way to do it.

I tried t%24

But the issue now is that t is a real and not a integer.

7

From the doc for remainder(numer, denom) can be written in stan as


remainder = numer - round(numer/denom) * denom

Ref: http://www.cplusplus.com/reference/cmath/remainder/

2 Likes
8

I thought I can’t use round too inside Julia/C++ but actually yes it is working well like that! Many thanks, @anon75146577 Daniel, really helpful. ;)

9

It’s a Stan function too!
https://mc-stan.org/docs/2_23/functions-reference/step-functions.html

2 Likes
10

Yeah, forgot to remind about this!