My current progress is on this stanc3 branch.

The parser accepts closures and creates corresponding callable objects. The object has an API for accessing the captured vars but it doesn’t do anything yet. I didn’t test anything but as long as you capture only data it hopefully works.

Now `integrate_ode_bdf`

always invokes the variadic `ode_bdf_tol`

and here’s a math branch where `ode_bdf`

uses the closure object API.

So, this compiles

```
data {
int<lower = 0> N; // number of measurement times
real ts[N]; // measurement times > 0
vector[2] y_init; // initial measured populations
real<lower = 0> y[N, 2]; // measured populations
}
transformed data {
real x_r[0];
int x_i[0];
}
parameters {
real<lower = 0> alpha;
real<lower = 0> beta;
real<lower = 0> gamma;
real<lower = 0> delta;
vector<lower = 0>[2] z_init;// initial population
real<lower = 0> sigma[2]; // measurement errors
}
transformed parameters {
functions
vector dz_dt(real t, vector z) {
real du_dt = (alpha - beta * z[2]) * z[1];
real dv_dt = (-gamma + delta * z[1]) * z[2];
return [du_dt, dv_dt]';
}
vector[2] z[N] = integrate_ode_bdf(dz_dt, z_init, 0.0, ts);
}
model {
alpha ~ normal(1, 0.5);
gamma ~ normal(1, 0.5);
beta ~ normal(0.05, 0.05);
delta ~ normal(0.05, 0.05);
sigma ~ lognormal(-1, 1);
z_init ~ lognormal(log(10), 1);
for (k in 1:2) {
y_init[k] ~ lognormal(log(z_init[k]), sigma[k]);
y[ , k] ~ lognormal(log(z[, k]), sigma[k]);
}
}
```

The generated C++ class is

```
template <typename captured_t__>
class closure0__ {
const captured_t__ alpha;
const captured_t__ beta;
const captured_t__ delta;
const captured_t__ gamma;
public:
using captured_scalar_t__ = captured_t__;
const int num_vars__;
closure0__(const captured_t__ alpha__, const captured_t__ beta__,
const captured_t__ delta__, const captured_t__ gamma__)
: alpha(alpha__), beta(beta__), delta(delta__),
gamma(gamma__), num_vars__(stan::is_var<captured_t__>::value ? 1 + 1 + 1
+ 1 : 0) { }
template <typename T0__, typename T1__>
Eigen::Matrix<typename boost::math::tools::promote_args<captured_t__,
typename boost::math::tools::promote_args<T0__,
T1__>::type>::type, -1, 1>
operator()(const T0__& t, const Eigen::Matrix<T1__, -1, 1>&
z, std::ostream* pstream__) const {
return closure0_impl__(t, z, alpha, beta, delta, gamma, pstream__);
}
class ValueOf__ {
const double alpha;
const double beta;
const double delta;
const double gamma;
public:
ValueOf__(const closure0__& init) : alpha(value_of(init.alpha)),
beta(value_of(init.beta)), delta(value_of(init.delta)),
gamma(value_of(init.gamma)) { }
template <typename T0__, typename T1__>
Eigen::Matrix<typename boost::math::tools::promote_args<T0__,
T1__>::type, -1, 1>
operator()(const T0__& t, const Eigen::Matrix<T1__, -1, 1>&
z, std::ostream* pstream__) const {
return closure0_impl__(t, z, alpha, beta, delta, gamma, pstream__);
}
};
class DeepCopy__ {
const captured_t__ alpha;
const captured_t__ beta;
const captured_t__ delta;
const captured_t__ gamma;
public:
DeepCopy__(const closure0__& init) : alpha(deep_copy_vars(init.alpha)),
beta(deep_copy_vars(init.beta)), delta(deep_copy_vars(init.delta)),
gamma(deep_copy_vars(init.gamma)) { }
template <typename T0__, typename T1__>
Eigen::Matrix<typename boost::math::tools::promote_args<captured_t__,
typename boost::math::tools::promote_args<T0__,
T1__>::type>::type, -1, 1>
operator()(const T0__& t, const Eigen::Matrix<T1__, -1, 1>&
z, std::ostream* pstream__) const {
return closure0_impl__(t, z, alpha, beta, delta, gamma, pstream__);
}
void accumulate_adjoints(double*) const { /* alpha *//* beta */
/* delta *//* gamma */}
};
void set_zero_adjoints() const { /* alpha *//* beta *//* delta */
/* gamma */ }
void accumulate_adjoints(double*) const { /* alpha *//* beta *//* delta */
/* gamma */}
void save_varis(vari**) const { /* alpha *//* beta *//* delta */
/* gamma */}
};
```

Filling in the missing parts shouldn’t be too hard.

One annoying thing I noticed is that when `ode_bdf`

calls the function it puts the output stream as the third argument but most other callers put the stream in the last place. Not a problem for the above code where the function takes only three arguments but in general we need to do something about the inconsistency.