Signal imputation in ODE's

Modeling
1

Dear Stan-Community,

I am working on my first ODE translation from R (using package deSolve) to Stan as a Stan newbie. So, far, the model works with a constant input, but generally I want to input a signal as a time-series into the model. Here is the working example in R first (with some truncations for better readability):

Mader <- function(t, y, parms, input){
  with(as.list(c(parms, y)), {

    # input
    watt <- input(t) / mm
    
    # ...
    
    #return
    list(c(dgp, dvo2, dlam, dlab), 
         atp = atp,  # ...
)
  })
}

# ...

signal <- crossing(
    watt = seq(150, 350, 50),
    dt = 1,
    tstage = 5,
    tphase = seq(0, tstage*60, dt)
  ) %>%
  mutate(
    times = cumsum(dt) - dt,
    tmin = times/60
    )

times <- signal$times

sigimp <- approxfun(signal$times, signal$watt, rule = 2)

## Solve model
out <- ode(
  y = yini, 
  times = times,
  func = Mader, 
  parms = parms,
  method = "rk4",
  input = sigimp)

I guess I adopted everything but the sigimp-function. Here is my attempt for implementation in Stan:

functions {
  real sigimp(real t_out, real[] ts, real[] watt) {
    int left = 1;
    int right = 1;
    real w = 1.0;
    int N_in = size(ts);
    int N_out = size(t_out);
    real watt_out;
    
    while (ts[right] < t_out) {
      right = right + 1;
    }
    while (ts[left + 1] < t_out) {
      left = left + 1;
    }
    w = (ts[right] - t_out) / (ts[right] - ts[left]);
    watt_out = w * watt[left] + (1 - w) * watt[right];
  
    return watt_out;
}
  
  real[] mader(real t, real[] y, real[] parms, real[] x_r, int[] x_i) {
    
    // variables
    int T = x_i[1];
    real ts[T] = x_r[1];
    real watt[T] = x_r[2];
    real w = sigimp(t, ts, watt);
    
  // [...]
  }
}
data {
  int<lower = 0> T;
  real t0;
  real<lower = t0> ts[T];
  real y0[5];
  real parms[23];
  real watt[T];
}

transformed data {
  real x_r[1] = ts;
  real x_r[2] = watt;
  int x_i[1] = T;
}
parameters {
  real<lower = 0, upper = 1> p;
}
model {
}
generated quantities {
  real y[T, 5] = integrate_ode_rk45(mader, y0, t0, ts, parms, x_r, x_i);
}

When compiling the code within R using cmdstanR I get the following error message:

Semantic error in line 26, column 4 to column 24:
   -------------------------------------------------
    24:      // variables
    25:      int T = x_i[1];
    26:      real ts[T] = x_r[1];
             ^
    27:      real watt[T] = x_r[2];
    28:      real w = sigimp(t, ts, watt);
   -------------------------------------------------

My questions:

  1. Do you have any generic examples how to write and implement such a function in Stan as in the provided example with the deSolve package?

  2. If my attempt is close before working, any hints to resolve remaining issues?

Any feedback/criticism regarding the code itself or tips to improve are welcome! =)

Regards
Rob

2

Hi there,

I now updated the code to the newer stan syntax and found a workaround. The real and integer data does not need to be packed in arrays to be passed to the newer ´ode_rk45´ function. However, I did not manage to implement an external interpolation function in the manual ode function, so I wrote it directly in the model function as a workaround. Here is the code:


functions {
  
  vector mader(real t, vector y, array[] real parms, array[] real watt, array[] real ts, int T) {
    
    //signal input
    int id = 2;
    real weight; 
    real w;
    
    while (ts[id] < t) { id = id + 1; }
    if(id > T) { id = T; }

    weight = (ts[id] - t) / (ts[id] - ts[id-1]);
    w = weight * watt[id-1] + (1 - weight) * watt[id];
    
// ...
    
    vector[5] dydt = [datp, dgp, dvo2, dlam, dlab]';
    
    return dydt;
  }
}
data {
  int<lower = 0> T;
  real t0;
  array[T] real<lower = t0> ts;
  array[T] real watt;
  vector[5] y0;
  array[23] real parms;
}
transformed data {
}
parameters {
  real<lower = 0, upper = 1> p;
}
model {
}
generated quantities {
  array[T] vector[5] y = ode_rk45(mader, y0, t0, ts, parms, watt,ts, T);
}

If you have any other suggestions to improve the code, please let me know!

The next problem I have is, that there are some parameters calculated I want to export besides the states from the ODE. An example on thow to do it in R with the deSolve package is given in the first post. Do you have any suggestions? Until now I have tried to not call the integrate function iteratively over the time vector. Are there any solutions to this similar to the provided example in R?

3

I think what you describe is what I have put on the forum under the name of a forcing function for an ode.

4

Thanks for your answer, very appreciated! =)

I have seen your example here beforehand and those by @charlesm93 and tried to mimic that function in the first post using linear interpolation. I guess I was close but I lack of experience in stan. The code in my second post works but without using a separate function.

There are some other hurdles now, but I guess I better open separate issues for them…