Design propoal: fold and static events

Developers
1

A bunch of us working on PKPD (@billg, @cbdavis33 , @charlesm93 ) have been thinking about add the features from the Torsten fork to Stan. I came to the following design. The application in mind is with ODE solution and the de-facto PKPD data input from NONMEM’s events specification, but the idea is general for other applications with the same recursive structure. The NONMEM model basically specifies a set of longitudinal events that an ODE system will go through, such as resetting initial condition, modifying the RHS, postprocessing the solution at the time, etc.

Consider a function f that takes argument a, x, y,\dots, and returns the same type as a, we want to support the calculation of (the intermediate as well as final results)

f(f(f(a, x_1, y_1, \dots), x_2, y_2, \dots), \dots)

This amounts to applying a hypothetical “walk” (or “fold”, “reduce” in some langs with the ambiguity whether the intermediate results are saved) function to the recursion, i.e.

walk(f, init, array1, array2, array3, ...)

that is equivalent to

f(f(init, array1[1], array2[1], array3[1], ...), array1[2], array2[2], array3[2],...)...

The arrays contain the same type of elements as x_i, y_i, etc, respectively.

The walk function can also supports arguments that remain the same in every recursion, e.g.

walk(f, init, array1, array2, theta1, theta2)

This is equivalent to

f(f(init, array1[1], array2[1], theta1, theta2), array1[2], array2[2], theta1, theta2)...

The classical example is cumsum

int n;
vector[n] array_to_sum_up;
vector[n] cumsum;
cumsum = walk(add, 0, array_to_sum_up)

To come back to the ODE events, consider a 1-dimentional ODE with initial condition x(t=t_0)=x_0 and right-hand-side f

//t_x = [t, x]
  functions {
    vector solve_f(vector t_x, real t_next, real theta) {
      real t = t_x[1];
      real x = t_x[2];
      vector sol[1], x1[1], ts[1];
      x1[1] = x;
      ts[1] = t_next
      sol = ode_rk45(f, x1, t, ts, theta);
      return [t_next, sol[1]]';
    }
  }
  //... create an array of repeated theta or allow single arg in "reduce"
  array[ns] real ts;
  array[ns] vector[2] sol;
  sol = walk(solve_f, [t0, x0]', ts, theta)

Actual PKPD models would require more sophisticated solve_f to address different kinds of events (i.e. array1, array2, … contain information that require if-else treatment) but it doesn’t change the structure of folding a recursion.

So I want to get a sanity check here on the idea. Happy to hear your thoughts!

1 Like
2

I like it, plus if your function has some good vibes we can add some sugar and have it moonwalk().