calling integrate_ode_rk45 from function causes compilation error #471

chvandorp · 2019-05-08T20:10:15Z

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

The text was updated successfully, but these errors were encountered:

wds15 · 2019-05-08T20:40:45Z

You should not see a compiler error, rather a stan compiler error message. The issue is that temporaries created inside a function are automatically converted to var, but the x_r and x_i must have their base type being double and int.

I am not sure if this is still being fixed for the current parser, but certainly for stanc3. Thanks for the bug report.

syclik · 2019-05-09T14:34:56Z

@wds15, were you able to reproduce the bug? It'd help if you stated that it was confirmed to be a bug.

Do you know if it's a language bug or a Math bug? (meaning the C++ instantiation of the function should be valid in Math, but it currently doesn't handle it?)

wds15 · 2019-05-09T14:36:48Z

I haven't confirmed this, but these things looks familiar to me. It's a parser problem which is created due to over-propagation of just about anything to var in our functions (even temporaries which are clearly created from constants).

rok-cesnovar · 2020-02-24T14:01:23Z

The model still produces a c++ error with stanc3 and develop cmdstan:

stan/src/stan/model/indexing/lvalue.hpp:51:5: error: no match for ‘operator=’ (operand types are ‘std::vector<stan::math::var>’ and ‘std::vector<double>’)
   x = std::forward<U>(y);
     ^
In file included from /usr/include/c++/5/vector:69:0,
                 from /usr/include/c++/5/bits/random.h:34,
                 from /usr/include/c++/5/random:49,
                 from /usr/include/c++/5/bits/stl_algo.h:66,
                 from /usr/include/c++/5/algorithm:62,
                 from stan/lib/stan_math/lib/eigen_3.3.3/Eigen/Core:269,
                 from stan/lib/stan_math/lib/eigen_3.3.3/Eigen/Dense:1,
                 from stan/lib/stan_math/stan/math/prim/fun/Eigen.hpp:13,
                 from stan/lib/stan_math/stan/math/rev.hpp:4,
                 from stan/lib/stan_math/stan/math.hpp:148,
                 from stan/src/stan/model/model_header.hpp:4,
                 from examples/bernoulli/testing.hpp:3:
/usr/include/c++/5/bits/vector.tcc:167:5: note: candidate: std::vector<_Tp, _Alloc>& std::vector<_Tp, _Alloc>::operator=(const std::vector<_Tp, _Alloc>&) [with _Tp = stan::math::var; _Alloc = std::allocator<stan::math::var>]
     vector<_Tp, _Alloc>::
     ^
/usr/include/c++/5/bits/vector.tcc:167:5: note:   no known conversion for argument 1 from ‘std::vector<double>’ to ‘const std::vector<stan::math::var>&’
In file included from /usr/include/c++/5/vector:64:0,
                 from /usr/include/c++/5/bits/random.h:34,
                 from /usr/include/c++/5/random:49,
                 from /usr/include/c++/5/bits/stl_algo.h:66,
                 from /usr/include/c++/5/algorithm:62,
                 from stan/lib/stan_math/lib/eigen_3.3.3/Eigen/Core:269,
                 from stan/lib/stan_math/lib/eigen_3.3.3/Eigen/Dense:1,
                 from stan/lib/stan_math/stan/math/prim/fun/Eigen.hpp:13,
                 from stan/lib/stan_math/stan/math/rev.hpp:4,
                 from stan/lib/stan_math/stan/math.hpp:148,
                 from stan/src/stan/model/model_header.hpp:4,
                 from examples/bernoulli/testing.hpp:3:
/usr/include/c++/5/bits/stl_vector.h:448:7: note: candidate: std::vector<_Tp, _Alloc>& std::vector<_Tp, _Alloc>::operator=(std::vector<_Tp, _Alloc>&&) [with _Tp = stan::math::var; _Alloc = std::allocator<stan::math::var>]
       operator=(vector&& __x) noexcept(_Alloc_traits::_S_nothrow_move())
       ^
/usr/include/c++/5/bits/stl_vector.h:448:7: note:   no known conversion for argument 1 from ‘std::vector<double>’ to ‘std::vector<stan::math::var>&&’
/usr/include/c++/5/bits/stl_vector.h:470:7: note: candidate: std::vector<_Tp, _Alloc>& std::vector<_Tp, _Alloc>::operator=(std::initializer_list<_Tp>) [with _Tp = stan::math::var; _Alloc = std::allocator<stan::math::var>]
       operator=(initializer_list<value_type> __l)
       ^
/usr/include/c++/5/bits/stl_vector.h:470:7: note:   no known conversion for argument 1 from ‘std::vector<double>’ to ‘std::initializer_list<stan::math::var>’

nhuurre · 2020-02-24T20:41:18Z

The original issue has been fixed.
The problem @rok-cesnovar reports is not related to integrate_ode. Here's a small failing example

parameters {
  real x;
}
model {
  real y0[2] = {0.0, 1.0}; // double rhs, stan::math::var lhs
  y0 ~ normal(x, 1);
}

This has been broken since stan-dev/stan#2885.

rok-cesnovar · 2020-02-24T20:44:58Z

Thanks @nhuurre

SteveBronder · 2020-02-25T23:06:36Z

I have a branch here in stan that fixes the lvalue code for the above, however I think this may also be a compiler issue

The generated code with the branch above still throws an error (full hpp here) where the important bit that fails is

examples/blah2.hpp:435:12: error: no match for ‘operator=’ (operand types are 
‘std::vector<stan::math::var>’ and ‘std::vector<double>’)
         y0 = stan::math::array_builder<double>().add(0.0).add(1.0).array();
         ~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

from line 154:

        std::vector<local_scalar_t__> y0;
        y0 = std::vector<local_scalar_t__>(2, 0);

        current_statement__ = 2;
        y0 = stan::math::array_builder<double>().add(0.0).add(1.0).array();

So I think the solution here is to either

Do array_builder<local_scalar_t__> inside of the model block
Use pass types when declaring y0 so it's actually a std::vector<double>

SteveBronder · 2020-02-25T23:50:04Z

Or
3. Use assign for these

rok-cesnovar · 2020-06-03T13:07:09Z

This works now on develop. Closing.

rok-cesnovar transferred this issue from stan-dev/stan Feb 24, 2020

rok-cesnovar mentioned this issue Feb 24, 2020

Allow rvalue and assign to accept eigen expressions stan-dev/stan#2885

Merged

3 tasks

SteveBronder mentioned this issue Feb 25, 2020

Add SFINAE for vector types to assign stan-dev/stan#2889

Closed

3 tasks

rok-cesnovar closed this as completed Jun 3, 2020

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

calling integrate_ode_rk45 from function causes compilation error #471

calling integrate_ode_rk45 from function causes compilation error #471

chvandorp commented May 8, 2019

wds15 commented May 8, 2019

syclik commented May 9, 2019

wds15 commented May 9, 2019

rok-cesnovar commented Feb 24, 2020

nhuurre commented Feb 24, 2020

rok-cesnovar commented Feb 24, 2020

SteveBronder commented Feb 25, 2020

SteveBronder commented Feb 25, 2020

rok-cesnovar commented Jun 3, 2020

calling integrate_ode_rk45 from function causes compilation error #471

calling integrate_ode_rk45 from function causes compilation error #471

Comments

chvandorp commented May 8, 2019

Summary:

Description:

Reproducible Steps:

Current Output:

Expected Output:

PyStan Version:

Python Version:

Operating System:

wds15 commented May 8, 2019

syclik commented May 9, 2019

wds15 commented May 9, 2019

rok-cesnovar commented Feb 24, 2020

nhuurre commented Feb 24, 2020

rok-cesnovar commented Feb 24, 2020

SteveBronder commented Feb 25, 2020

SteveBronder commented Feb 25, 2020

rok-cesnovar commented Jun 3, 2020