Slightly awkward syntax when passing data arguments from function to function


#1

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;
  }
}

#2

I think the array_builder error has been fixed in develop.


#3

I cannot confirm that.

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

compiles to

stan::math::assign(y_hat, 
  integrate_ode_bdf(
    pk_rhs_functor__(), 
    stan::model::rvalue(x_r,
      stan::model::cons_list(
        stan::model::index_min_max(17, 17), 
        stan::model::nil_index_list()), "x_r"), 
    get_base1(x_r,18,"x_r",1), 
    stan::model::rvalue(x_r, 
      stan::model::cons_list(
        stan::model::index_min_max(9, 16), 
        stan::model::nil_index_list()), "x_r"), 
    theta, 
    static_cast<std::vector<local_scalar_t__> >(
      stan::math::array_builder<local_scalar_t__ >().add(
        get_base1(x_r,19,"x_r",1)).array()), 
    x_i, pstream__));

The problematic part isn’t array builder itself. The problem is that local_scalar_t__ is stan::math::var and the built array therefore becomes a std::vector<stan::math::var> while integrate_ode_bdf expects a std::vector<double>. This happens with the develop branch I pulled today.

The relevant part of the error message is this:

src/cmdstan/main.cpp:8:50:   required from here
./../tnfalpha/mdls/gt-cpd/oral_1cmt_mm_clearance_re.hpp:170:52: 
error: no matching function for call to 
‘integrate_ode_bdf(
  oral_1cmt_mm_clearance_re_model_namespace::pk_rhs_functor__, 
  stan::model::rvalue_return<std::vector<double>, 
  stan::model::cons_index_list<stan::model::index_min_max, 
  stan::model::nil_index_list> >::type, const double&, 
  stan::model::rvalue_return<std::vector<double>, 
  stan::model::cons_index_list<stan::model::index_min_max, 
  stan::model::nil_index_list> >::type, std::vector<stan::math::var>&, 
  std::vector<stan::math::var>, const std::vector<int>&, std::ostream*&)’

Here it is also visible that x_r[17:17] get’s nicely translated to a std::vector<double>, but { x_r[19] } does not. But maybe that is also simply expected behaviour in terms of how stanc translates the code.