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[1];      // dimension of ODEs
  dydt[1] = - theta[1] * y[1];                   // 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[1]; //one initial value
  real theta[1]; //one parameter next to A0
  real y[1,1]; //ODE solution
  y0[1]=A0; //initial values
  theta[1]=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[1] dydt;
    dydt[1] = -kr[1] * y[1];
    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>[1] kr;    // rate constant to be estimated
  vector[1]            y0;    // initial condition to be estimated
}
model {
  vector[1] 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[1] ~ normal(0.2, 0.02);                      // prior for rate constant
  y0 ~ normal(100,10);                            // prior on initial condition
  y0_obs ~ normal(y0[1], 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[1],   # 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[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[1,1]  fo_ode2(t,A0,k1) {
        real<lower=0> k1;
        vector[1] y0;        //one initial value
        y0[1]=A0;            //initial value
        vector [1] 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[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) {
        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),
                  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!