#### Summary:
When I call `integrate_ode_rk45` from the `model` block, I can …use literals for the mandatory real and integer data arguments `x_r` and `x_i`. So for instance, I can write `integrate_ode_rk45(ode, y0, t0, Ts, theta, {0.0}, {0})`. If I try to do this in an auxiliary function defined in the `functions` block, compilation of the stan model fails with a `CompileError`.
#### Description:
I find it convenient to use an auxiliary function (defined in the functions block) to integrate ODEs for my stan model (for instance in combination with the `map_rect` function). However, I found that one has to be careful with the real and integer data arguments of e.g. `integrate_ode_rk45`. If I use such an auxiliary function, these have to be defined in the `data` or `transformed data` blocks, and I can't pass dummy literal values (e.g. `{0.0}` instead of `x_r`).
The Stan "pre-compiler" does not give an error message. Instead the C compiler fails (without a clear description of the error).
#### Reproducible Steps:
As an example, I've included the following model. This version causes the `CompileError`.
If I comment out the call to `bad_auxiliary_function` in the `model` block, and use the `good_auxiliary_function` instead, all works fine. I can also call `bad_auxiliary_function` from the `generated quantities` block without any errors.
```
// ode_integrator_bug.stan
functions {
real[] ode(real t, real[] y, real[] theta, real[] x_r, int[] x_i) {
real dy[2];
dy[1] = theta[1] * y[2];
dy[2] = -y[1];
return dy;
}
real[] bad_auxiliary_function(data real[] ts, real omega) {
real theta[1]; real y0[2]; real result[size(ts), 2];
theta = {omega};
y0 = {0.0, 1.0};
result = integrate_ode_rk45(ode, y0, 0.0, ts, theta, {0.0}, {0});
return result[:, 1];
}
real[] good_auxiliary_function(data real[] ts, real omega, data real[] x_r, int[] x_i) {
real theta[1]; real y0[2]; real result[size(ts), 2];
theta = {omega};
y0 = {0.0, 1.0};
result = integrate_ode_rk45(ode, y0, 0.0, ts, theta, x_r, x_i);
return result[:, 1];
}
}
data {
int N;
real Ts[N]; // time points
real Xs[N]; // observations
}
transformed data {
// mandatory data for ODE integrator
real x_r[0];
int x_i[0];
}
parameters {
real<lower=0> omega;
real<lower=0> sigma;
}
model {
real xhats[N]; real theta[1]; real y0[2]; real result[N, 2];
// integrate ODEs in model block...
theta = {omega};
y0 = {0.0, 1.0};
result = integrate_ode_rk45(ode, y0, 0.0, Ts, theta, {0.0}, {0}); // no problem here...
// likelihood
//Xs ~ normal(result[:, 1], sigma);
// alternatively: use an auxiliary function to call the integrator
//xhats = good_auxiliary_function(Ts, omega, x_r, x_i); // works fine...
xhats = bad_auxiliary_function(Ts, omega); // compilation fails
// likelihood
Xs ~ normal(xhats, sigma);
// prior
omega ~ normal(1, 0.1);
sigma ~ normal(0, 1);
}
generated quantities {
real xhats[N];
xhats = bad_auxiliary_function(Ts, omega); // no problem here
}
```
For easy reproducibility, here's my python script to compile and run the model.
```
#!/usr/bin/env python
# coding: utf-8
import pystan
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
## compile the stan model
with open("ode_integrator_bug.stan", 'r') as f:
stan_model = f.read()
sm = pystan.StanModel(model_code=stan_model, verbose=True)
## create some data
N = 20
Ts = np.linspace(1e-2, 3*np.pi, N)
Xs = np.sin(Ts) + scipy.stats.norm.rvs(loc=0, scale=0.1, size=N)
## run the model
samples = sm.sampling(data={"N" : N, "Ts" : Ts, "Xs" : Xs}, iter=1000, chains=1)
chain = samples.extract(permuted=True)
xhats = chain["xhats"]
## make a figure
fig, ax = plt.subplots(1, 1, figsize=(7, 5))
ax.plot(Ts, xhats.T, color='gray', zorder=1)
ax.scatter(Ts, Xs, color='k', zorder=2, label='data')
fig.savefig("oscillations.pdf", bbox_inches='tight')
```
#### Current Output:
This is the full error reported by python:
```
Traceback (most recent call last):
File "/usr/lib/python3.6/distutils/unixccompiler.py", line 118, in _compile
extra_postargs)
File "/usr/lib/python3.6/distutils/ccompiler.py", line 909, in spawn
spawn(cmd, dry_run=self.dry_run)
File "/usr/lib/python3.6/distutils/spawn.py", line 36, in spawn
_spawn_posix(cmd, search_path, dry_run=dry_run)
File "/usr/lib/python3.6/distutils/spawn.py", line 159, in _spawn_posix
% (cmd, exit_status))
distutils.errors.DistutilsExecError: command 'x86_64-linux-gnu-gcc' failed with exit status 1
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File "stan-ode-bug.py", line 12, in <module>
sm = pystan.StanModel(model_code=stan_model, verbose=True)
File "/usr/local/lib/python3.6/dist-packages/pystan/model.py", line 349, in __init__
build_extension.run()
File "/usr/lib/python3.6/distutils/command/build_ext.py", line 339, in run
self.build_extensions()
File "/usr/lib/python3.6/distutils/command/build_ext.py", line 448, in build_extensions
self._build_extensions_serial()
File "/usr/lib/python3.6/distutils/command/build_ext.py", line 473, in _build_extensions_serial
self.build_extension(ext)
File "/usr/lib/python3.6/distutils/command/build_ext.py", line 533, in build_extension
depends=ext.depends)
File "/usr/lib/python3.6/distutils/ccompiler.py", line 574, in compile
self._compile(obj, src, ext, cc_args, extra_postargs, pp_opts)
File "/usr/lib/python3.6/distutils/unixccompiler.py", line 120, in _compile
raise CompileError(msg)
distutils.errors.CompileError: command 'x86_64-linux-gnu-gcc' failed with exit status 1
```
The same problem occurs when I use **cmdstan** to `make` the model. This error might contain more useful information:
```
$ make ~/stan_ode_bug/ode_integrator_bug
--- Translating Stan model to C++ code ---
bin/stanc /home/chvandorp/stan_ode_bug/ode_integrator_bug.stan --o=/home/chvandorp/stan_ode_bug/ode_integrator_bug.hpp
Model name=ode_integrator_bug_model
Input file=/home/chvandorp/stan_ode_bug/ode_integrator_bug.stan
Output file=/home/chvandorp/stan_ode_bug/ode_integrator_bug.hpp
--- Linking C++ model ---
g++ -Wall -I . -isystem stan/lib/stan_math/lib/eigen_3.3.3 -isystem stan/lib/stan_math/lib/boost_1.66.0 -isystem stan/lib/stan_math/lib/sundials_3.1.0/include -std=c++1y -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE -DBOOST_DISABLE_ASSERTS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -Wno-unused-function -Wno-uninitialized -I src -isystem stan/src -isystem stan/lib/stan_math/ -DFUSION_MAX_VECTOR_SIZE=12 -Wno-unused-local-typedefs -DEIGEN_NO_DEBUG -DNO_FPRINTF_OUTPUT -pipe src/cmdstan/main.cpp -O3 -o /home/chvandorp/stan_ode_bug/ode_integrator_bug -include /home/chvandorp/stan_ode_bug/ode_integrator_bug.hpp stan/lib/stan_math/lib/sundials_3.1.0/lib/libsundials_nvecserial.a stan/lib/stan_math/lib/sundials_3.1.0/lib/libsundials_cvodes.a stan/lib/stan_math/lib/sundials_3.1.0/lib/libsundials_idas.a
In file included from <command-line>:0:0:
/home/chvandorp/stan_ode_bug/ode_integrator_bug.hpp: In instantiation of ‘std::vector<typename boost::math::tools::promote_args<RT1, RT2>::type> ode_integrator_bug_model_namespace::bad_auxiliary_function(const std::vector<T>&, const T1__&, std::ostream*) [with T0__ = double; T1__ = stan::math::var; typename boost::math::tools::promote_args<RT1, RT2>::type = stan::math::var; std::ostream = std::basic_ostream<char>]’:
/home/chvandorp/stan_ode_bug/ode_integrator_bug.hpp:435:61: required from ‘T__ ode_integrator_bug_model_namespace::ode_integrator_bug_model::log_prob(std::vector<T_l>&, std::vector<int>&, std::ostream*) const [with bool propto__ = true; bool jacobian__ = true; T__ = stan::math::var; std::ostream = std::basic_ostream<char>]’
stan/src/stan/model/log_prob_grad.hpp:43:13: required from ‘double stan::model::log_prob_grad(const M&, std::vector<double>&, std::vector<int>&, std::vector<double>&, std::ostream*) [with bool propto = true; bool jacobian_adjust_transform = true; M = ode_integrator_bug_model_namespace::ode_integrator_bug_model; std::ostream = std::basic_ostream<char>]’
stan/src/stan/services/util/initialize.hpp:167:18: required from ‘std::vector<double> stan::services::util::initialize(Model&, stan::io::var_context&, RNG&, double, bool, stan::callbacks::logger&, stan::callbacks::writer&) [with bool Jacobian = true; Model = ode_integrator_bug_model_namespace::ode_integrator_bug_model; RNG = boost::random::additive_combine_engine<boost::random::linear_congruential_engine<unsigned int, 40014, 0, 2147483563>, boost::random::linear_congruential_engine<unsigned int, 40692, 0, 2147483399> >]’
stan/src/stan/services/diagnose/diagnose.hpp:54:29: required from ‘int stan::services::diagnose::diagnose(Model&, stan::io::var_context&, unsigned int, unsigned int, double, double, double, stan::callbacks::interrupt&, stan::callbacks::logger&, stan::callbacks::writer&, stan::callbacks::writer&) [with Model = ode_integrator_bug_model_namespace::ode_integrator_bug_model]’
src/cmdstan/command.hpp:149:57: required from ‘int cmdstan::command(int, const char**) [with Model = ode_integrator_bug_model_namespace::ode_integrator_bug_model]’
src/cmdstan/main.cpp:8:50: required from here
/home/chvandorp/stan_ode_bug/ode_integrator_bug.hpp:119:54: error: no matching function for call to ‘integrate_ode_rk45(ode_integrator_bug_model_namespace::ode_functor__, std::vector<stan::math::var>&, double, const std::vector<double>&, std::vector<stan::math::var>&, std::vector<stan::math::var>, std::vector<int>, std::ostream*&)’
stan::math::assign(result, integrate_ode_rk45(ode_functor__(), y0, 0.0, ts, theta, static_cast<std::vector<local_scalar_t__> >(stan::math::array_builder<local_scalar_t__ >().add(0.0).array()), static_cast<std::vector<int> >(stan::math::array_builder<int >().add(0).array()), pstream__));
~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
In file included from stan/lib/stan_math/stan/math/prim/arr.hpp:44:0,
from stan/lib/stan_math/stan/math/prim/mat.hpp:325,
from stan/lib/stan_math/stan/math/rev/mat.hpp:12,
from stan/lib/stan_math/stan/math.hpp:4,
from stan/src/stan/model/model_header.hpp:4,
from /home/chvandorp/stan_ode_bug/ode_integrator_bug.hpp:3,
from <command-line>:0:
stan/lib/stan_math/stan/math/prim/arr/functor/integrate_ode_rk45.hpp:70:1: note: candidate: template<class F, class T1, class T2> std::vector<std::vector<typename stan::return_type<T_shared_param, T_job_param>::type> > stan::math::integrate_ode_rk45(const F&, const std::vector<T>&, double, const std::vector<double>&, const std::vector<T_l>&, const std::vector<double>&, const std::vector<int>&, std::ostream*, double, double, int)
integrate_ode_rk45(const F& f, const std::vector<T1>& y0, double t0,
^~~~~~~~~~~~~~~~~~
stan/lib/stan_math/stan/math/prim/arr/functor/integrate_ode_rk45.hpp:70:1: note: template argument deduction/substitution failed:
In file included from <command-line>:0:0:
/home/chvandorp/stan_ode_bug/ode_integrator_bug.hpp:119:54: note: cannot convert ‘stan::math::array_builder<T>::array() [with T = stan::math::var]()’ (type ‘std::vector<stan::math::var>’) to type ‘const std::vector<double>&’
stan::math::assign(result, integrate_ode_rk45(ode_functor__(), y0, 0.0, ts, theta, static_cast<std::vector<local_scalar_t__> >(stan::math::array_builder<local_scalar_t__ >().add(0.0).array()), static_cast<std::vector<int> >(stan::math::array_builder<int >().add(0).array()), pstream__));
~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
make/models:14: recipe for target '/home/chvandorp/stan_ode_bug/ode_integrator_bug' failed
make: *** [/home/chvandorp/stan_ode_bug/ode_integrator_bug] Error 1
```
#### Expected Output:
The model works fine if I comment out the call to `bad_auxiliary_function` in the model block and use the `good_auxiliary_function` instead.
#### PyStan Version:
2.18.1.0
#### Python Version:
Python 3.6.7
Cython 0.29.4
#### Operating System:
Ubuntu 18.04.2 LTS
GCC 8.2.0