# Signal imputation in ODE's

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,
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

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?

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