Slightly awkward syntax when passing data arguments from function to function

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


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;

I think the array_builder error has been fixed in develop.

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::model::index_min_max(17, 17), 
        stan::model::nil_index_list()), "x_r"), 
        stan::model::index_min_max(9, 16), 
        stan::model::nil_index_list()), "x_r"), 
    static_cast<std::vector<local_scalar_t__> >(
      stan::math::array_builder<local_scalar_t__ >().add(
    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
error: no matching function for call to 
  stan::model::nil_index_list> >::type, const double&, 
  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.

1 Like

Stanc should not translate to C++ and fail to parse.

The model you posted builds for me as is, and also when I replace x_r[19:19] with { x_r[19] }.

   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] }, x_i);

Do you have an example that translates to C++ then fails to compile under the current develop branch? If so, it’d be super helpful if you could share it, because I can’t reproduce the bug you’re reporting.

Hi Bob!
Yes, I do have an example. I’ve uploaded it here. One version compiles perfectly (model.stan) while the other (model-fail.stan) shows the behaviour I described above. The model gets translated to C++ but fails to compile because of the type of local_scalar_t__.

I compiled this with the develop branch checked out today. I’m working on Ubuntu 18.04 LTS and g++ -v outputs

Using built-in specs.
Target: x86_64-linux-gnu
Configured with: ../src/configure -v --with-pkgversion='Ubuntu 7.3.0-16ubuntu3' --with-bugurl=file:///usr/share/doc/gcc-7/README.Bugs --enable-languages=c,ada,c++,go,brig,d,fortran,objc,obj-c++ --prefix=/usr --with-gcc-major-version-only --with-as=/usr/bin/x86_64-linux-gnu-as --with-ld=/usr/bin/x86_64-linux-gnu-ld --program-suffix=-7 --program-prefix=x86_64-linux-gnu- --enable-shared --enable-linker-build-id --libexecdir=/usr/lib --without-included-gettext --enable-threads=posix --libdir=/usr/lib --enable-nls --with-sysroot=/ --enable-clocale=gnu --enable-libstdcxx-debug --enable-libstdcxx-time=yes --with-default-libstdcxx-abi=new --enable-gnu-unique-object --disable-vtable-verify --enable-libmpx --enable-plugin --enable-default-pie --with-system-zlib --with-target-system-zlib --enable-objc-gc=auto --enable-multiarch --disable-werror --with-arch-32=i686 --with-abi=m64 --with-multilib-list=m32,m64,mx32 --enable-multilib --with-tune=generic --enable-offload-targets=nvptx-none --without-cuda-driver --enable-checking=release --build=x86_64-linux-gnu --host=x86_64-linux-gnu --target=x86_64-linux-gnu
Thread model: posix
gcc version 7.3.0 (Ubuntu 7.3.0-16ubuntu3) 

Tell me if I can provide anything more.

[EDIT]: Maybe worth adding that I compile this with cmdstan and thread support, i.e. I’m adding

CC = g++

to make/local.

All repositories (cmdstan, stan and math) are checked out on their develop branches.

Thanks. I validated that’s indeed a bug in the parser and have opened an issue. Hopefully we’ll get this fixed before the next release. See the Stan GitHub issue #2565.

1 Like

Great! Glad I could help.