# New cmdstan ode_solvers and brms

Short summary of the problem
I would like to use brms in combination with the new ode_solvers in Stan via cmdstanr. Ode_solvers can be called from brms is possible with the old solvers as nicely demonstrated by Markus Gesmann (Fitting multivariate ODE models with brms | mages' blog). I can use the new ode_solvers without problem directly in Stan using cmdstanr but when I try to couple this to brms, I fail and I can not figure out what I am doing wrong. To demonstrate my problem I simulated a simple problem: an exponential decay (c=A0exp(-k1time) with two parameters, initial concentration A0 and rate constant k1. This has of course an analytical solution, so I know what should come out (once I get this going I will apply it to more complex problems without an analytical solution; and of course I could do it without brms but I just would like to know what I am doing wrong so that I better understand what is going on, and also because of the nice interface of brms). I am struggling with a syntax error and I cannot get my head around it. It probably has to do with declaration of vectors and real dimensions. I tried many combinations without result. I am just hoping that someone can see what I am doing wrong. The following code shows first the simulated data, a working code that couples brms to the old ode solver, a working stan code that uses the new ode_solvers, and then the attempt to couple brms to the new ode_solvers that does not work. Any suggestion will be much appreciated, many thanks in advance!

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

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

# the next lines of code show the coupling of brms to the old integrators, works without a problem
fo_model <- "
real[] ode_fo(real t, //time
real [] y,         // the rates
real [] theta,     // the parameters
real [] x_r,       // data constant (not used)
int[] x_i){        // data constant (not used)
real dydt;      // dimension of ODEs
dydt = - theta * y;                   // first ODE
return dydt;        // returns a 3-dimensional array
}

// this is the function call from brms, integration of ODEs:
real fo_ode(real t, real A0, real k1) {
real y0; //one initial value
real theta; //one parameter next to A0
real y[1,1]; //ODE solution
y0=A0; //initial values
theta=k1;

y=integrate_ode_rk45(ode_fo, y0,0,rep_array(t,1),theta,rep_array(0.0,0), rep_array(1,1), 0.00001,0.00001,100);
// Return relevant values
return(y[1,1]);
}

"
# t=0 not allowed, will be estimated

df_sim_ode <- df_sim[-1,]

fo_formula <- bf(y_obs~fo_ode(time,A0,k1),
A0~1,
k1~1,
nl=TRUE)

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

fo_result1 <- 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_model, block="functions"), file="fo_result1")

# the following code uses the new integrators directly in Stan, works also fine

write("
functions {                  // input to the ode_integrator:
vector fo(real t,          // time
vector y,        // state
vector kr        // parameters (only one in this case, c0 is automatically estimated)
) {
vector dydt;
dydt = -kr * y;
return dydt;
}
}
data {
int<lower = 1>   N;
real             t0;
real<lower = t0> ts[N];      // each element must be greater than t0
real             y0_obs;     // the observed initial value
real             y_obs[N];   //the observed (measured) values excluding the initial value
}
parameters {
real<lower = 0>      sigma; // experimental sigma
vector<lower = 0> kr;    // rate constant to be estimated
vector            y0;    // initial condition to be estimated
}
model {
vector mu[N] = ode_bdf(fo, y0, t0, ts, kr);  // call to the ode_integrator, result is put in mu
sigma ~ cauchy(0, 5);                           // prior for sigma
kr ~ normal(0.2, 0.02);                      // prior for rate constant
y0 ~ normal(100,10);                            // prior on initial condition
y0_obs ~ normal(y0, sigma);
for (t in 1:N) {
y_obs[t] ~ normal(mu[t], sigma);
}
}

","stan_fo.stan")

# compile the stan model:
fo_model1 <- cmdstan_model("stan_fo.stan")
# the data list for Stan:
fo_data <- list(
N = N_t - 1,                # the first datapoint should not be at time zero because that point is provided by t0 and y0
t0 = 0,
ts = df_sim\$time[-1],       # omitting the first time datapoint at time zero because that needs to be given for y0 (next line)
y0_obs = df_sim\$y_obs,   # the first datapoint is the initial value y0
y_obs = df_sim\$y_obs[-1]    # the measured data without the starting point
)

# sampling from the compiled model:
fo_fit <- fo_model1\$sample(
data = fo_data,
seed = 42L,
refresh = 0,
parallel_chains=4
)

# So far, so good but when coupled to brms things go wrong:

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

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

// function called from brms: integration of ODEs
vector[1,1]  fo_ode2(t,A0,k1) {
real<lower=0> k1;
vector y0;        //one initial value
y0=A0;            //initial value
vector  y_hat[N]=ode_rk45(fo2, y0,0.0,time,k1);
return(y[1,1]);
}

"

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

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

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")

# The software complains about syntax errors
``````
• Operating System: Windows 10
• brms Version: 2.15.0
• cmdstanr version 0.4.0
• R version 4.0.3

I don’t know if that changes anything, but could you try if it works in the latest brms release 2.16.1? I am asking because version 2.16 contains a lot of changes that perhaps also affect your problem.

Thanks Paul, I installed the update, which is a good move anyway, but it did not help. I don’t think it is a brms problem, it is just me not knowing how to call the solvers in the right way. It does not help that I get this cryptic error message “Error: An error occured during compilation! See the message above for more information.” But unfortunately there is no above message, so it is not clear what I am doing wrong.

try out `make_stancode(..., parse = TRUE)` instead of `brm(...)` to see if then a more informative error is printed.

Unfortunately not, same cryptic message. For the time being, I will have to work with pure Stan, I guess, though I find it intriguing that brms does work nicely with the old ode integrators…

Nope… brms is just fine. You just need to get used as to how the brms stan source generation works. The make_stancode thing is the way to go… although I don’t use the parse option as it is easier first to dump that into a file and check it out. Here is what samples for me just fine:

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

vector dydt;          // dimension of ODE
dydt = - k1 * y;   // 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) {
vector y0 = rep_vector(vA0, 1);        //one initial value
int N = size(time);
vector y_hat_ode[N] = ode_rk45(fo2, y0 ,0.0,to_array_1d(time),vk1);
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),
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")

``````
1 Like

Many thanks wds15! That works with me as well. I was indeed struggling with how to couple brms to stan. Still many things to learn, such help from the Stan community is fantastic! Much appreciated!