Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update ODE conversion case study #222

Merged
merged 1 commit into from
Sep 12, 2023
Merged

Update ODE conversion case study #222

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
69 changes: 33 additions & 36 deletions knitr/convert-odes/convert_odes.Rmd
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,17 +36,18 @@ system function -- in the new interface none of this packing and unpacking
is necessary.

The new interface also uses `vector` variables for the state rather than
`real[]` variables.
`array[] real` variables.

As an example, in the old solver interface the system function for the ODE
$y' = -\alpha y$ would be written:

```{stan, output.var = "", eval = FALSE}
```stan
functions {
real[] rhs(real t, real[] y, real[] theta, real[] x_r, int[] x_i) {
array[] real rhs(real t, array[] real y, array[] real theta,
array[] real x_r, array[] int x_i) {
real alpha = theta[1];

real yp[1] = { -alpha * y[1] };
array[1] real yp = {-alpha * y[1]};

return yp;
}
Expand Down Expand Up @@ -100,7 +101,7 @@ The new `ode_bdf` solver interface is (the interfaces for `ode_adams` and
`ode_rk45` are the same):

```{stan, output.var = "", eval = FALSE}
vector[] ode_bdf(F f, vector y0, real t0, real[] times, ...)
vector[] ode_bdf(F f, vector y0, real t0, array[] real times, ...)
```

The arguments are:
Expand Down Expand Up @@ -144,7 +145,7 @@ The `ode_bdf_tol` interface is (the interfaces for `ode_rk45_tol`
and `ode_adams_tol` are the same):

```{stan, output.var = "", eval = FALSE}
vector[] ode_bdf_tol(F f, vector y0, real t0, real[] times,
array[] vector ode_bdf_tol(F f, vector y0, real t0, array[] real times,
real rel_tol, real abs_tol, int max_num_steps, ...)
```

Expand Down Expand Up @@ -193,7 +194,7 @@ The two models here come from the Stan
### ODE System Function

In the old SIR system function, `beta`, `kappa`, `gamma`, `xi`, and `delta`,
are packed into the `real[] theta` argument. `kappa` isn't actually a model
are packed into the `array[] real theta` argument. `kappa` isn't actually a model
parameter so it is not clear why it is packed in with the other parameters,
but it is. Promoting `kappa` to a parameter causes there to be more states
in the extended ODE sensitivity system (used to get gradients of the ODE
Expand All @@ -208,14 +209,11 @@ functions {
// theta[3] = gamma, recovery rate
// theta[4] = xi, bacteria production rate
// theta[5] = delta, bacteria removal rate
real[] simple_SIR(real t,
real[] y,
real[] theta,
real[] x_r,
int[] x_i) {
real dydt[4];

dydt[1] = - theta[1] * y[4] / (y[4] + theta[2]) * y[1];
array[] real simple_SIR(real t, array[] real y, array[] real theta,
array[] real x_r, array[] int x_i) {
array[4] real dydt;

dydt[1] = -theta[1] * y[4] / (y[4] + theta[2]) * y[1];
dydt[2] = theta[1] * y[4] / (y[4] + theta[2]) * y[1] - theta[3] * y[2];
dydt[3] = theta[3] * y[2];
dydt[4] = theta[4] * y[2] - theta[5] * y[4];
Expand All @@ -230,7 +228,7 @@ rewritten to explicitly name all the parameters. No separation of data
and parameters is necessary either -- the solver will not add more
equations for arguments that are defined in the `data` and
`transformed data` blocks. The state variables in the new model are also
represented by `vector` variables instead of `real[]` variables.
represented by `vector` variables instead of `array[]` variables.
The new ODE system function is:

```{stan, output.var = "", eval = FALSE}
Expand All @@ -256,14 +254,14 @@ functions {

### Calling the ODE Solver

In the old ODE interface, the parameters are all packed into a `real[]` array
In the old ODE interface, the parameters are all packed into a `array[] real`
before calling the ODE solver:

```{stan, output.var = "", eval = FALSE}
transformed parameters {
real<lower=0> y[N_t, 4];
array[N_t, 4] real<lower=0> y;
{
real theta[5] = {beta, kappa, gamma, xi, delta};
array[5] real theta = {beta, kappa, gamma, xi, delta};
y = integrate_ode_rk45(simple_SIR, y0, t0, t, theta, x_r, x_i);
}
}
Expand All @@ -272,11 +270,11 @@ transformed parameters {
In the new ODE interface each of the arguments is appended on to the ODE
solver function call. The RK45 ODE solver with default tolerances is called
`ode_rk45`. Because the states are handled as `vector` variables, the solver
output is an array of vectors (`vector[]`).
output is an array of vectors (`array[] vector`).

```{stan, output.var = "", eval = FALSE}
transformed parameters {
vector<lower=0>[4] y[N_t] = ode_rk45(simple_SIR, y0, t0, t,
array[N_t] vector<lower=0>[4] y = ode_rk45(simple_SIR, y0, t0, t,
beta, kappa, gamma, xi, delta);
}
```
Expand All @@ -285,7 +283,7 @@ The `ode_rk45_tol` function can be used to specify tolerances manually:

```{stan, output.var = "", eval = FALSE}
transformed parameters {
vector<lower=0>[4] y[N_t] = ode_rk45_tol(simple_SIR, y0, t0, t,
array[N_t] vector<lower=0>[4] = ode_rk45_tol(simple_SIR, y0, t0, t,
1e-6, 1e-6, 1000,
beta, kappa, gamma, xi, delta);
}
Expand All @@ -308,12 +306,10 @@ In the old system function, parameters are manually unpacked from `theta` and

```{stan, output.var = "", eval = FALSE}
functions {
real[] one_comp_mm_elim_abs(real t,
real[] y,
real[] theta,
real[] x_r,
int[] x_i) {
real dydt[1];
array[] real one_comp_mm_elim_abs(real t, array[] real y,
array[] real theta, array[] real x_r,
array[] int x_i) {
array[1] real dydt;
real k_a = theta[1]; // Dosing rate in 1/day
real K_m = theta[2]; // Michaelis-Menten constant in mg/L
real V_m = theta[3]; // Maximum elimination rate in 1/day
Expand All @@ -322,8 +318,9 @@ functions {
real dose = 0;
real elim = (V_m / V) * y[1] / (K_m + y[1]);

if (t > 0)
dose = exp(- k_a * t) * D * k_a / V;
if (t > 0) {
dose = exp(-k_a * t) * D * k_a / V;
}

dydt[1] = dose - elim;

Expand Down Expand Up @@ -365,14 +362,14 @@ the `x_i` argument is required even though it isn't used:

```{stan, output.var = "", eval = FALSE}
transformed data {
real x_r[2] = {D, V};
int x_i[0];
array[2] real x_r = {D, V};
array[2] int x_i;
}
...
transformed parameters {
real C[N_t, 1];
array[N_t, 1] real C;
{
real theta[3] = {k_a, K_m, V_m};
array[3] real theta = {k_a, K_m, V_m};
C = integrate_ode_bdf(one_comp_mm_elim_abs, C0, t0, times, theta, x_r, x_i);
}
}
Expand All @@ -383,7 +380,7 @@ In the new interface the arguments are simply passed at the end of the

```{stan, output.var = "", eval = FALSE}
transformed parameters {
vector[1] mu_C[N_t] = ode_bdf(one_comp_mm_elim_abs, C0, t0, times,
array[N_t] vector[1] mu_C = ode_bdf(one_comp_mm_elim_abs, C0, t0, times,
k_a, K_m, V_m, D, V);
}
```
Expand All @@ -392,7 +389,7 @@ The `ode_bdf_tol` function can be used to specify tolerances manually:

```{stan, output.var = "", eval = FALSE}
transformed parameters {
vector[1] mu_C[N_t] = ode_bdf_tol(one_comp_mm_elim_abs, C0, t0, times,
array[N_t] vector[1] mu_C = ode_bdf_tol(one_comp_mm_elim_abs, C0, t0, times,
1e-8, 1e-8, 1000,
k_a, K_m, V_m, D, V);
}
Expand Down
Loading