Merge pull request #222 from stan-dev/bmw/case-studies-updates

Update ODE conversion case study

Author: jgabry

knitr/convert-odes/convert_odes.Rmd

@@ -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;
   }
@@ -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:
@@ -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, ...)
 ```
 
@@ -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
@@ -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];
@@ -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}
@@ -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);
   }
 }
@@ -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);
 }
 ```
@@ -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);
 }
@@ -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
@@ -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;
 
@@ -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);
   }
 }
@@ -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);
 }
 ```
@@ -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);
 }