Possible bug in expose_functions(..) with new ODE interface

I’m simulating a system of differential equations (ODEs). Initially I used integrate_ode_rk45 as tutorials suggested. I then save obtained posterior and do a simulation for new conditions.

According to the deprecation error I changed ODE function from integrate_ode_rk45 to ode_rk45. However, I can’t run expose_functions(..) with the new ODE interface, which is necesary to simulate with new conditions.

An example to reproduce the error is below. (ODEs are taken from here):

library(brms)
library(tidyverse)
library(here)
library(cmdstanr)
library(rstan)

rstan_options(auto_write=TRUE)
options(mc.cores=parallel::detectCores())

# simulated data for exponential decay: c=A0*exp(-k1*time)
set.seed(1)
time <- c(0,1,2,3,4,5,7,9,11)
y_exp <- 100*exp(-0.2*time)  # c0 =100, kr=0.2, exact simulated data
df_sim <- data.frame(time, y_exp)
df_sim <- df_sim %>% mutate(y_obs=y_exp+rnorm(length(time),0,5)) %>% dplyr::select(-y_exp)# random normal error added to simulated data with mean = 0 and sd = 5
N_t <- length(df_sim$time) # number of datapoints

fo_model2 <- "
  vector fo2(real t,                  //time
             vector y,                // the state
             real k1){                // the parameter next to A0

             vector[1] dydt;          // dimension of ODE
             dydt[1] = - k1 * y[1];   // ODE
             return dydt;             // returns a 1-dimensional vector
  }

// function called from brms: integration of ODEs
  vector fo_ode2(data vector time, vector vA0, vector vk1, int N) {
        vector[1] y0 = rep_vector(vA0[1], 1);        //one initial value
        //int N = size(time);
        vector[1] y_hat_ode[N] = ode_rk45(fo2, y0 ,0.0,to_array_1d(time),vk1[1]);
        vector[N] y_hat;
        for(i in 1:N) y_hat[i] = y_hat_ode[i,1];
  return(y_hat);
}
"

df_sim_ode <- df_sim[-1,] #remove the first datapoint at t0

fo_formula <- bf(y_obs~fo_ode2(time,A0,k1, 8),
                 A0~1,
                 k1~1,
                 nl=TRUE, loop=FALSE)

fo_priors <- c(prior(normal(100,10),nlpar=A0),
               prior(normal(0.2,0.1), nlpar=k1),
               prior(cauchy(0,10), class=sigma)
)

fo_result2 <- brm(
  data=df_sim_ode,
  family = gaussian,
  formula=fo_formula,
  prior=fo_priors, 
  inits=0,
  iter=4000,
  chains=4, 
  cores=4, 
  stanvars = stanvar(scode=fo_model2, block="functions"), 
  backend="cmdstanr", 
  file="fo_result2")

# Everything above runs correctly
# This function causes the error
expose_functions(fo_result2, vectorize = FALSE)

Produces this error:

SYNTAX ERROR, MESSAGE(S) FROM PARSER:
Variable "fo2" does not exist.
 error in 'modelf4cf9c8e4d57_User_defined_functions' at line 18, column 45
  -------------------------------------------------
    16:         vector[1] y0 = rep_vector(vA0[1], 1);        //one initial value
    17:         //int N = size(time);
    18:         vector[1] y_hat_ode[N] = ode_rk45(fo2, y0 ,0.0,to_array_1d(time),vk1[1]);
                                                    ^
    19:         vector[N] y_hat;
  -------------------------------------------------

Error in stanc(model_code = mc, model_name = "User-defined functions",  : 
  failed to parse Stan model 'User-defined functions' due to the above error.

I’m using:
R version 4.1.0 (2021-05-18)
Platform: x86_64-pc-linux-gnu (64-bit)
brms_2.18.0
cmdstanr_0.5.2

I’m wondering if anyone knows if it’s a bug or am I doing something wrong?

After some time I’ve realized that the error appears due to the old version of RStan. I managed to find a workaround without updating RStan using a tutorial on expose_cmdstanr_functions() Exposing Stan user-defined functions using CmdStanR and Rcpp (rok-cesnovar.github.io)