Hi,

I recently learned about Stan’s `map_rect`

and I’m seriously excited about it. Since I work a lot with ODE models, I think I might just be the target group for this feature.

Below, there’s an example of how I would like to use an ODE integrator inside a function (`unit_logLik`

) passed to `map_rect`

. The data argument `x_r`

contains the data, observation times, initial value, initial time and the dose data used in the ODE.

When I wanted to pass the data to `integrate_ode_bdf`

I came upon slightly awkward behaviour. Since the initial value and the parameter I wanted to pass on to my ODE right hand side are only single real values, I first tried passing in

`integrate_ode_bdf(pk_rhs, x_r[17], x_r[18], x_r[9:16], theta, x_r[19], x_i)`

That however, doesn’t work of course since `x_r[17]`

and `x_r[19]`

evaluate to doubles and not to vectors, as expected. My next try was

`integrate_ode_bdf(pk_rhs, { x_r[17] }, x_r[18], x_r[9:16], theta, { x_r[19] }, x_i)`

which works fine for `stanc`

but failes during compilation. As far as my limited knowledge of C++ could take me, the problem is that the Stan code `{ x_r[19] }`

gets compiled to

`stan::math::array_builder<local_scalar_t__>(...)`

where `local_scalar_t__`

is `stan::math::var`

, therefore returning a `std::vector<stan::math::var>`

instead of a `std::vector<double>`

as expected by `integrate_ode_bdf`

. This wasn’t a problem for the initial value, since those can be parameters, I assume. The only solution that worked in the end was to replace `{ x_r[17] }`

and `{ x_r[19] }`

by `x_r[17:17]`

and `x_r[19:19]`

.

It’s not a big deal to do this, but it feels slightly awkward that the first way doesn’t work (even though I understand from the C++ code why it cannot). I’m not sure if `stanc`

could be capable of capturing a mistake like this and warning the user about it? Or I guess it should at least be found in the manual. Which it maybe is? I didn’t find anything…

```
functions {
real[] pk_rhs(real t,
real[] y,
real[] theta,
real[] x_r,
int[] x_i) {
real ka = theta[1];
real f = theta[2];
real Vmax = theta[3];
real Km = theta[4];
real V = theta[5];
real dose = x_r[1];
return { (ka * f * dose * exp(-ka * (t + 2))
- Vmax * y[1] / (Km + y[1])) / V };
}
vector unit_logLik(vector shared,
vector unit_spec,
real[] x_r,
int[] x_i) {
real theta[5] = to_array_1d(shared[1:5]);
real sigma = shared[6];
int n_obs = x_i[1];
real y_hat[n_obs, 1] =
integrate_ode_bdf(pk_rhs, x_r[17:17],
x_r[18], x_r[9:16],
theta, x_r[19:19], x_i);
vector[n_obs] lp;
for (i in 1:n_obs) {
lp[i] = normal_lpdf(x_r[i] | log(y_hat[i,1]), sigma);
}
return lp;
}
}
```