Compilation error on using the Differential-Algebraic equation (DAE) solver on time that depends on the parameters

Modeling
1

Summary

I’m trying to solve a system of differential equations using the dae solver which is available on Stan, but I encounter a compilation issue when I try to evaluate a system at a time that depends on the parameters.

The Issue

To illustrate what I mean, I made use of the dae on a system of the form
S_2' + S_1 = 0
S_2 - S_1' = 0

Because Stan requires me to have a model and datasets, I’ve implemented a linear regression model with dummy data.

The model file looks like the following:

functions {
  vector residual(real time, vector state, vector deriv) {
    vector[2] res;
    res[1] = deriv[2] + state[1];
    res[2] = state[2] - deriv[1];

    return res;
  }
}

data {
  int N1;
  array[N1] real x;
  array[N1] real yexp;
  array[N1] real error;
}

transformed data {
}

parameters {
  real m;
  real b;
}

transformed parameters {
  array[N1] real y;
  for (i in 1:N1) {
     y[i] = m*x[i] + b;
   }

  array[N1] vector[2] S;
  S = dae(residual, [1, 1]', [1, -1]', 0.0, {1, 2, 3});

  print(S);
}

model {
  m ~ normal(0, 10);
  b ~ normal(0, 10);

  yexp ~ normal(y, error);
}

generated quantities {
}

The key to the issue here is the function call. In this case the model works as expected, and it print the correct solution of the previous system.

However, if we replace the time, which in the previous case is given by the array {1,2,3} by something that depends on the parameters, e.g.: {m, 2*m, 3*m}, then the program gives an error at compile time.
Using cmdstanpy to compile the program reveals what I believe to be the source of the error, which is repeated across the tracelog multiple times:

cannot convert 'std::vector<stan::math::var_value<double>, 
std::allocator<stan::math::var_value<double> > >
(std::initializer_list<stan::math::var_value<double> 
>{((const stan::math::var_value<double>*)
(& const stan::math::var_value<double> 
[3]{m, stan::math::operator*<int>(2, m), stan::math::operator*<int>(3, m)})), 3},
std::allocator<stan::math::var_value<double> >())' 
(type 'std::vector<stan::math::var_value<double>, 
std::allocator<stan::math::var_value<double> > >') 
to type 'const std::vector<double, std::allocator<double> >&'

(new lines were added on the previous code block to improve readability)

The full tracelog when compiling the previous model with the modification indicated above using cmdstanpy is available in the following file:
log.txt (29.3 KB)

It seems to me that the dae was coded in such a way as that the time which is provided to it is forced to be a constant, and therefore cannot vary as the parameters vary during the sampling.

My use case

Obviously the previous model is simply for debugging purposes, but I do in fact have a system where I would like to compute the solution of a differential algebraic system of equations for a specific value of “time”, where that time depends on the parameters.

In my case I have a system which looks like
(E^2 - 2 \lambda) e^{\lambda/E^2} = \Omega_m (1+z)^3 + \Omega_r (1+z)^4 ,

where \Omega_m and \Omega_r are parameters, and \lambda is derived from the previous two parameters.
Here z plays the role of time, given that it is an independent variable, and I’m trying to obtain the value of E for each corresponding z.

Additionally I am also interesting in knowing the values of the integral
\int_0^z 1/E(x) dx ,

which thanks to the user nhuurre in a previous post, can be coded in using the following residual function

vector residual(real z, vector state, vector deriv, real lambda, real Omega_m) {
    real E = state[1];
    real lhs = (E^2 - 2*lambda)*exp(lambda/E^2);   // left hand side of the equation for E
    real rhs = Omega_m*(1+z)^3 + Omega_r*(1+z)^4;  // right hand side of the equation for E

    // compute the residuals
    vector[2] res;
    res[1] = rhs - lhs;       // E satisfies the algebraic relation, meaning that res[1] should be 0
    res[2] = 1/E - deriv[2];  // the derivative of S₂ is in fact 1/E, this should also be 0

    return res;
  }

and can be solved by called the dae with

S = dae(residual, [1, 0.0]', [deriv_0, 1]', 0.0, redshift, lambda, Omega_m);

where the initial conditions where placed by hand, expect for deriv_0 which was computed previously, and redshift is an array of values of z which were declared in the data block, meaning that they are constant.

The problem now arises because I wish to compute the value of this time for a single value of “time”, which I refer to as z_*, that is a function (only) of the parameters, just like the MWE before.

How will I be able to use Stan to solve this numerical problem?

Further unrelated issue

To make matters a little worse, I’m also interesting in computing the value of yet another integral, which looks like
\int_{z_*}^\infty \frac{1}{E(x)*f(x)} dx

where f(x) is a function that depends on both the “time” and the parameters.

My idea was to integrate the previous integrand numerically, which required that on each step of the integrator a call to the dae to obtain the value of E would have to be performed, which in my mind could increase significantly the code run-time, and one of the things which I enjoyed when I used the dae was how fast it was.

I’m not quite sure if there is another approach or not, and hence I’m mentioning this issue here.

References

Stan Functions Reference - 11.3 (DAE Solver)

2

I don’t think there’s a simple way to solve this. What you’ve encountered should be a bug. The current DAE solver implementation doesn’t support time to be a parameter.

EDIT: on second thought IIUC you’d like to integrate u(z)\equiv 1/E(z). Would converting the equation to an ODE work? As it is can be turned into the canonical form of ODE:

h(u)e^{g(u)}=f(z) \quad \longrightarrow \quad u'=\frac{f'}{h_ue^g+he^gg_u}\equiv F(u, z, ...).

though one must be careful about possible singularities where the denominator \rightarrow 0.

3

A possible workaround: consider the identity

\int_{0}^{z}f\left(x\right)\mathrm{d}x=\int_{0}^{1}zf\left(zx\right)\mathrm{d}x

Rescale the time parameter and the derivative.

// integrate x from 0 to 1 <--> integrate z from 0 to z_star
vector residual(real x, vector state, vector deriv, real lambda, real Omega_m, real z_star) {
    real z = x*z_star;
    real E = state[1];
    real lhs = (E^2 - 2*lambda)*exp(lambda/E^2);   // left hand side of the equation for E
    real rhs = Omega_m*(1+z)^3 + Omega_r*(1+z)^4;  // right hand side of the equation for E

    // compute the residuals
    vector[2] res;
    res[1] = rhs - lhs;       // E satisfies the algebraic relation, meaning that res[1] should be 0
    res[2] = z_star/E - deriv[2];  // the derivative of S₂ is in fact 1/E, this should also be 0

    return res;
}
4

I also thought so, but there might be a good reason as to why this is required to be a constant.
Regardless, it would be nice if the solver was more dynamic and would accept time that depended on the parameters.

Thank you for the suggestion!

Just to ensure we are on the same page, you are suggesting that I write an ODE of that form, use the ODE solver in Stan to compute for u(z) \equiv 1/E(z) and then integrate u(z) it with the integral solver which is also available in Stan, correct?

Because if so that would also allow me to easily integrate another expression that I have which is of the form
\int_\infty^z \frac{1}{E(x)f(x)} dx
where the function f(x) also depends on the parameters.

If I got the idea correct, it seems to be that there is a disadvantage to it, which is the fact that if I make two separate integrals, one for \int_0^z 1/E(x) dx and another for \int_\infty^z \frac{1}{E(x)f(x)} dx it will perform redundant computations on 1/E(z) obtained from the ODE solver, which might significantly slow down the code.
Is there a way to cache the results of a function?

5

Yes that is also a good idea, nice and simple!
yizhang also suggested a different approach, I will see which one works the best for my use case.
Thanks!

6

@yizhang - for the compilation error, is the times parameter of the dae function intended to be data only? If so we should promote that to a stanc3 compile time error.

7

Yes. The c++ impl considers time only as data.

Edit: I see it’s been fixed. Thanks!

1 Like
8

Is this meant to be the standard behavior or are you planning to add the ability to use time dependent parameters to your roadmap? (or the third option, is this something which is still an open question?).

9

This is the fix of the stan transpiler so the error will be caught and use will get meaningful msg saying the time must be data in the DAE solver. We can add supporting time as param on the roadmap but currently I lack time for it.

10

I’m not sure whether I should make a new thread, given that the initial problem was “solved” (in the sense that it is not expected to use time that depend on the parameters as I was trying to do).
Please let me know if you think that I should make the previous comment by yizhang as a solution as well.

Anyways, I’m trying the approach you mentioned earlier to turn this problem into a canonical ODE, but the code is giving me an unknown error (literally) and I have no idea on what could be wrong (I’ve double checked my math and my code).

The model

First of all, I’ve omitted minor things during the previous description of my model, and I shall make the model clear here.

In my use case, the parameters in my system are
\theta = (h, \Omega_b, \Omega_c, \Omega_r)

and a quantity \lambda which can be derived from the parameters which is given by
\lambda = \frac{1}{2} + W_0\left( -\frac{\Omega_b + \Omega_c + \Omega_r}{2e^{1/2}} \right)

where W_0 is the 0-th branch of the Lambert function.

As for the data, the following two quantities will be used as datasets:
R = \sqrt{\Omega_b + \Omega_c} \int_0^{z_*} u(z) dx
w_b = \Omega_b*h^2

where u(z) \equiv 1/E(z), as you defined previously and z_* depends only on the parameters (the expression to compute its value from the parameters is fairly complicated and does not bring any meaningful information).

The experiment gave us that
R = 1.710 \pm 0.019
w_b = 0.0224 \pm 0.0001

Using the equations you have so kindly provided I obtained:
u' = \frac{e^{-\lambda u^2}}{\lambda u (u^{-2} - 2 \lambda) - u^{-3}} \left[\frac{3}{2}(\Omega_b + \Omega_c)(1+z)^2 + 2 \Omega_r(1+z)^3) \right]

To compute the vale of R I solve the previous differential equation using one of Stan’s ODE solver, to obtain u(z), and integrate it using one of Stan’s integration routines, to finally obtain \int_0^{z_*} u(z) dx.

Stan model file

/* block for user defined functions */
functions {
  // auxiliary function to solve for u(z) ≡ 1/E(z) using ODE solver
  vector ode(real z, vector y, array[] real faketheta) {
    real Omega_b = faketheta[2];
    real Omega_c = faketheta[3];
    real Omega_r = faketheta[4];
    real lambda = faketheta[5];
    real u = y[1];

    real f1 = exp(-lambda*u^2) / (lambda*u*(u^(-2)-2*lambda) - u^(-3));
    real f2 = 1.5*(Omega_b + Omega_c)*(1+z)^2 + 2*Omega_r*(1+z)^3;

    vector[1] output;
    output[1] = f1*f2;

    return output;
  }

  // returns u(z) ≡ 1/E(z) to integrate
  real u(real z, real xc, array[] real faketheta, array[] real x_r, array[] int x_i) {
    // initial conditions
    real z0 = 0;
    vector[1] init;
    init[1] = 1;

    // obtain u(z) using ODE solver
    vector[1] vec = ode_rk45(ode, init, z0, {z}, faketheta)[1];
    real sol = vec[1];

    return sol;
  }
}

/* block for datasets used */
data {
  real Rexp;
  real Rexp_err;
  real wbexp;
  real wbexp_err;
}

/* block for transformed data */
transformed data {
  // create null data values to give to integrate_1d because they are both required arguments
  array[0] real x_r;
  array[0] int x_i;
}

/* block for model parameters */
parameters {
  real<lower=0> h;
  real<lower=0> Omega_b;
  real<lower=0> Omega_c;
  real Omega_r;
}

/* block for transformed parameters */
transformed parameters {
  // compute λ
  real lambda = 0.5 + lambert_w0( -(Omega_b + Omega_c + Omega_r)/(2*exp(0.5)) );

  // parameter array (labeled as fake because it includes λ, which is not a true parameter)
  array[5] real faketheta = {h, Omega_b, Omega_c, Omega_r, lambda};

  // compute z_star
  real g1 = 0.0783*(Omega_b*h^2)^(-0.238)/(1 + 39.5*(Omega_b*h^2)^(-0.763));
  real g2 = 0.560/(1 + 21.1*(Omega_b*h^2)^1.81);
  real z_star = 1048 * (1 + 0.00124*(Omega_b*h^2)^(-0.738)) * (1 + g1*((Omega_b + Omega_c)*h^2)^g2);

  // get integral of u(z) ≡ 1/E(z) from 0 to z_star
  real intofu = integrate_1d(u, 0, z_star, faketheta, x_r, x_i);

  // wb
  real wb = Omega_b*h^2;

  // R
  real R = sqrt(Omega_b + Omega_c) * intofu;
}

/* block for model definition */
model {
  // priors
  h ~ normal(0.7, 10);
  Omega_b ~ normal(0.05, 10);
  Omega_c ~ normal(0.27, 10);
  Omega_r ~ normal(0, 10);

  // likelihood
  Rexp ~ normal(R, Rexp_err);
  wbexp ~ normal(wb, wbexp_err);
}

/* block for generated quantites */
generated quantities {
}

The error

Compiling and running the previous model using CmdStanPy (with show_console=True, a single chain and sensible initial values) gives the following output:

15:28:43 - cmdstanpy - INFO - compiling stan file /home/joseferreira/tmp/model.stan to exe file /home/joseferreira/tmp/model
15:29:06 - cmdstanpy - INFO - compiled model executable: /home/joseferreira/tmp/model
15:29:06 - cmdstanpy - INFO - Chain [1] start processing
Chain [1] method = sample (Default)
Chain [1] sample
Chain [1] num_samples = 1000 (Default)
Chain [1] num_warmup = 1000 (Default)
Chain [1] save_warmup = 0 (Default)
Chain [1] thin = 1 (Default)
Chain [1] adapt
Chain [1] engaged = 1 (Default)
Chain [1] gamma = 0.050000000000000003 (Default)
Chain [1] delta = 0.80000000000000004 (Default)
Chain [1] kappa = 0.75 (Default)
Chain [1] t0 = 10 (Default)
Chain [1] init_buffer = 75 (Default)
Chain [1] term_buffer = 50 (Default)
Chain [1] window = 25 (Default)
Chain [1] algorithm = hmc (Default)
Chain [1] hmc
Chain [1] engine = nuts (Default)
Chain [1] nuts
Chain [1] max_depth = 10 (Default)
Chain [1] metric = diag_e (Default)
Chain [1] metric_file =  (Default)
Chain [1] stepsize = 1 (Default)
Chain [1] stepsize_jitter = 0 (Default)
Chain [1] num_chains = 1 (Default)
Chain [1] id = 1 (Default)
Chain [1] data
Chain [1] file = /tmp/tmppcx6vxnf/2kwnz5li.json
Chain [1] init = /tmp/tmppcx6vxnf/zyatktdh.json
Chain [1] random
Chain [1] seed = 6437
Chain [1] output
Chain [1] file = /tmp/tmppcx6vxnf/modelobz90mtx/model-20221110152906.csv
Chain [1] diagnostic_file =  (Default)
Chain [1] refresh = 100 (Default)
Chain [1] sig_figs = -1 (Default)
Chain [1] profile_file = profile.csv (Default)
Chain [1] num_threads = 1 (Default)
Chain [1] 
Chain [1] 
15:29:06 - cmdstanpy - INFO - Chain [1] done processing
15:29:06 - cmdstanpy - ERROR - Chain [1] error: terminated by signal 11 Unknown error -11
Traceback (most recent call last):
  File "/home/joseferreira/tmp/mwe-cmdstanpy.py", line 18, in <module>
    fit = model.sample(data=data, show_console=True, chains=1, inits=inits)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/joseferreira/.micromamba/envs/cmdstanpy/lib/python3.11/site-packages/cmdstanpy/model.py", line 1188, in sample
    raise RuntimeError(msg)
RuntimeError: Error during sampling:

Command and output files:
RunSet: chains=1, chain_ids=[1], num_processes=1
 cmd (chain 1):
	['/home/joseferreira/tmp/model', 'id=1', 'random', 'seed=6437', 'data', 'file=/tmp/tmppcx6vxnf/2kwnz5li.json', 'init=/tmp/tmppcx6vxnf/zyatktdh.json', 'output', 'file=/tmp/tmppcx6vxnf/modelobz90mtx/model-20221110152906.csv', 'method=sample', 'algorithm=hmc', 'adapt', 'engaged=1']
 retcodes=[-11]
 per-chain output files (showing chain 1 only):
 csv_file:
	/tmp/tmppcx6vxnf/modelobz90mtx/model-20221110152906.csv
 console_msgs (if any):
	/tmp/tmppcx6vxnf/modelobz90mtx/model-20221110152906_0-stdout.txt
Consider re-running with show_console=True if the above output is unclear!

The important part is that Stan spits out an unknown error (terminated by signal 11 Unknown error -11) and exits, providing me no further information on what caused it.

Additionally, it speaks of two different files (/tmp/tmppcx6vxnf/modelobz90mtx/model-20221110152906.csv and /tmp/tmppcx6vxnf/modelobz90mtx/model-20221110152906_0-stdout.txt) which are nowhere to be found.
Although I suspect these two files will not bring any meaningful information

My main guess was that these initial values are in a situation where the denominator goes to zero, but that would trigger a warning, right?

My second guess was that integrating something which is being outputted by an ODE solver might give some sort of issues, but that shouldn’t be the case, right? Worst case scenario, it would simply take ages.

Besides the previous two rather wild guesses, I’m stuck.
Could this be due to Stan itself?

System

Everything took place in a new virtual environment, with packages installed from conda-forge
Python version: 3.11
CmdStanPy version: 1.0.8
CmdStan version: 2.30.1
OS: Ubuntu 18.04 (x64).
Compiler: gcc version 7.5.0 (Ubuntu 7.5.0-3ubuntu1~18.04) - however, when I installed CmdStanPy, it also installed a newer version of gcc (10.4.0), which does not seem to be available in my PATH. This is the one which shows up when I write gcc -v, but it is the one which is being used by CmdStan?

11

That model does seem to cause a segmentation fault, which is definitely worthy of a bug report. Can you share the data/inits you used to avoid any spurious issues?

12

Compiling with full debugging I got


Program received signal SIGSEGV, Segmentation fault.
boost::numeric::odeint::detail::for_each8<__gnu_cxx::__normal_iterator<double*, std::vector<double, std::allocator<double> > >, __gnu_cxx::__normal_iterator<double const*, std::vector<double, std::allocator<double> > >, __gnu_cxx::__normal_iterator<double const*, std::vector<double, std::allocator<double> > >, __gnu_cxx::__normal_iterator<double const*, std::vector<double, std::allocator<double> > >, __gnu_cxx::__normal_iterator<double const*, std::vector<double, std::allocator<double> > >, __gnu_cxx::__normal_iterator<double const*, std::vector<double, std::allocator<double> > >, __gnu_cxx::__normal_iterator<double const*, std::vector<double, std::allocator<double> > >, __gnu_cxx::__normal_iterator<double const*, std::vector<double, std::allocator<double> > >, boost::numeric::odeint::default_operations::scale_sum7<double, double, double, double, double, double, double> > (op=..., first8=0, first7=<error reading variable: Cannot access memory at address 0x8>, first6=<error reading variable: Cannot access memory at address 0x8>, first5=<error reading variable: Cannot access memory at address 0x8>, first4=<error reading variable: Cannot access memory at address 0x8>, first3=0, first2=0, last1=..., first1=0) at stan/lib/stan_math/lib/boost_1.78.0/boost/numeric/odeint/algebra/detail/for_each.hpp:90
90                  op( *first1++ , *first2++ , *first3++ , *first4++ , *first5++ , *first6++ , *first7++ , *first8++ );
(gdb) bt
#0  boost::numeric::odeint::detail::for_each8<__gnu_cxx::__normal_iterator<double*, std::vector<double, std::allocator<double> > >, __gnu_cxx::__normal_iterator<double const*, std::vector<double, std::allocator<double> > >, __gnu_cxx::__normal_iterator<double const*, std::vector<double, std::allocator<double> > >, __gnu_cxx::__normal_iterator<double const*, std::vector<double, std::allocator<double> > >, __gnu_cxx::__normal_iterator<double const*, std::vector<double, std::allocator<double> > >, __gnu_cxx::__normal_iterator<double const*, std::vector<double, std::allocator<double> > >, __gnu_cxx::__normal_iterator<double const*, std::vector<double, std::allocator<double> > >, __gnu_cxx::__normal_iterator<double const*, std::vector<double, std::allocator<double> > >, boost::numeric::odeint::default_operations::scale_sum7<double, double, double, double, double, double, double> > (op=..., first8=0, first7=<error reading variable: Cannot access memory at address 0x8>, first6=<error reading variable: Cannot access memory at address 0x8>, first5=<error reading variable: Cannot access memory at address 0x8>, first4=<error reading variable: Cannot access memory at address 0x8>, first3=0, first2=0, last1=..., first1=0) at stan/lib/stan_math/lib/boost_1.78.0/boost/numeric/odeint/algebra/detail/for_each.hpp:90
#1  boost::numeric::odeint::range_algebra::for_each8<std::vector<double, std::allocator<double> >, std::vector<double, std::allocator<double> > const, std::vector<double, std::allocator<double> > const, std::vector<double, std::allocator<double> > const, std::vector<double, std::allocator<double> > const, std::vector<double, std::allocator<double> > const, std::vector<double, std::allocator<double> > const, std::vector<double, std::allocator<double> > const, boost::numeric::odeint::default_operations::scale_sum7<double, double, double, double, double, double, double> > (op=..., s8=..., s7=std::vector of length 0, capacity 0, s6=std::vector of length 0, capacity 0, s5=std::vector of length 0, capacity 0, s4=std::vector of length 0, capacity 0, s3=..., s2=..., s1=std::vector of length 1, capacity 1 = {...}) at stan/lib/stan_math/lib/boost_1.78.0/boost/numeric/odeint/algebra/range_algebra.hpp:83
#2  boost::numeric::odeint::runge_kutta_dopri5<std::vector<double, std::allocator<double> >, double, std::vector<double, std::allocator<double> >, double, boost::numeric::odeint::range_algebra, boost::numeric::odeint::default_operations, boost::numeric::odeint::initially_resizer>::calc_state<std::vector<double, std::allocator<double> >, std::vector<double, std::allocator<double> >, std::vector<double, std::allocator<double> >, std::vector<double, std::allocator<double> >, std::vector<double, std::allocator<double> > > (t_new=0, deriv_new=..., t_old=0, deriv_old=..., x_old=..., x=std::vector of length 1, capacity 1 = {...}, t=6.4432420282839148e-273, this=0x7fffffffa290) at stan/lib/stan_math/lib/boost_1.78.0/boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp:267
#3  boost::numeric::odeint::dense_output_runge_kutta<boost::numeric::odeint::controlled_runge_kutta<boost::numeric::odeint::runge_kutta_dopri5<std::vector<double, std::allocator<double> >, double, std::vector<double, std::allocator<double> >, double, boost::numeric::odeint::range_algebra, boost::numeric::odeint::default_operations, boost::numeric::odeint::initially_resizer>, boost::numeric::odeint::default_error_checker<double, boost::numeric::odeint::range_algebra, boost::numeric::odeint::default_operations>, boost::numeric::odeint::default_step_adjuster<double, double>, boost::numeric::odeint::initially_resizer, boost::numeric::odeint::explicit_error_stepper_fsal_tag>, boost::numeric::odeint::explicit_controlled_stepper_fsal_tag>::calc_state<std::vector<double, std::allocator<double> > > (x=std::vector of length 1, capacity 1 = {...}, t=6.4432420282839148e-273, this=0x7fffffffa290) at stan/lib/stan_math/lib/boost_1.78.0/boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp:354
#4  boost::numeric::odeint::checked_stepper<boost::numeric::odeint::dense_output_runge_kutta<boost::numeric::odeint::controlled_runge_kutta<boost::numeric::odeint::runge_kutta_dopri5<std::vector<double, std::allocator<double> >, double, std::vector<double, std::allocator<double> >, double, boost::numeric::odeint::range_algebra, boost::numeric::odeint::default_operations, boost::numeric::odeint::initially_resizer>, boost::numeric::odeint::default_error_checker<double, boost::numeric::odeint::range_algebra, boost::numeric::odeint::default_operations>, boost::numeric::odeint::default_step_adjuster<double, double>, boost::numeric::odeint::initially_resizer, boost::numeric::odeint::explicit_error_stepper_fsal_tag>, boost::numeric::odeint::explicit_controlled_stepper_fsal_tag>, boost::numeric::odeint::max_step_checker, boost::numeric::odeint::dense_output_stepper_tag>::calc_state<std::vector<double, std::allocator<double> > > (x=std::vector of length 1, capacity 1 = {...}, t=6.4432420282839148e-273, this=<synthetic pointer>) at stan/lib/stan_math/lib/boost_1.78.0/boost/numeric/odeint/integrate/check_adapter.hpp:164
#5  boost::numeric::odeint::detail::integrate_times<boost::numeric::odeint::checked_stepper<boost::numeric::odeint::dense_output_runge_kutta<boost::numeric::odeint::controlled_runge_kutta<boost::numeric::odeint::runge_kutta_dopri5<std::vector<double, std::allocator<double> >, double, std::vector<double, std::allocator<double> >, double, boost::numeric::odeint::range_algebra, boost::numeric::odeint::default_operations, boost::numeric::odeint::initially_resizer>, boost::numeric::odeint::default_error_checker<double, boost::numeric::odeint::range_algebra, boost::numeric::odeint::default_operations>, boost::numeric::odeint::default_step_adjuster<double, double>, boost::numeric::odeint::initially_resizer, boost::numeric::odeint::explicit_error_stepper_fsal_tag>, boost::numeric::odeint::explicit_controlled_stepper_fsal_tag>, boost::numeric::odeint::max_step_checker, boost::numeric::odeint::dense_output_stepper_tag>, std::reference_wrapper<stan::math::coupled_ode_system<crash_model_namespace::ode_variadic2_functor__, double, std::vector<double, std::allocator<double> > const&> >, std::vector<double, std::allocator<double> >, __gnu_cxx::__normal_iterator<double*, std::vector<double, std::allocator<double> > >, double, boost::numeric::odeint::checked_observer<stan::math::ode_rk45_tol_impl<crash_model_namespace::ode_variadic2_functor__, Eigen::Matrix<double, -1, 1, 0, -1, 1>, double, double, std::vector<double, std::allocator<double> >, (void*)0>(char const*, crash_model_namespace::ode_variadic2_functor__ const&, Eigen::Matrix<double, -1, 1, 0, -1, 1> const&, double, std::vector<double, std::allocator<double> > const&, double, double, long, std::ostream*, std::vector<double, std::allocator<double> > const&)::{lambda(std::vector<double, std::allocator<double> > const&, double)#3}, boost::numeric::odeint::max_step_checker> >(boost::numeric::odeint::checked_stepper<boost::numeric::odeint::dense_output_runge_kutta<boost::numeric::odeint::controlled_runge_kutta<boost::numeric::odeint::runge_kutta_dopri5<std::vector<double, std::allocator<double> >, double, std::vector<double, std::allocator<double> >, double, boost::numeric::odeint::range_algebra, boost::numeric::odeint::default_operations, boost::numeric::odeint::initially_resizer>, boost::numeric::odeint::default_error_checker<double, boost::numeric::odeint::range_algebra, boost::numeric::odeint::default_operations>, boost::numeric::odeint::default_step_adjuster<double, double>, boost::numeric::odeint::initially_resizer, boost::numeric::odeint::explicit_error_stepper_fsal_tag>, boost::numeric::odeint::explicit_controlled_stepper_fsal_tag>, boost::numeric::odeint::max_step_checker, boost::numeric::odeint::dense_output_stepper_tag>, std::reference_wrapper<stan::math::coupled_ode_system<crash_model_namespace::ode_variadic2_functor__, double, std::vector<double, std::allocator<double> > const&> >, std::vector<double, std::allocator<double> >&, __gnu_cxx::__normal_iterator<double*, std::vector<double, std::allocator<double> > >, __gnu_cxx::__normal_iterator<double*, std::vector<double, std::allocator<double> > >, double, boost::numeric::odeint::checked_observer<stan::math::ode_rk45_tol_impl<crash_model_namespace::ode_variadic2_functor__, Eigen::Matrix<double, -1, 1, 0, -1, 1>, double, double, std::vector<double, std::allocator<double> >, (void*)0>(char const*, crash_model_namespace::ode_variadic2_functor__ const&, Eigen::Matrix<double, -1, 1, 0, -1, 1> const&, double, std::vector<double, std::allocator<double> > const&, double, double, long, std::ostream*, std::vector<double, std::allocator<double> > const&)::{lambda(std::vector<double, std::allocator<double> > const&, double)#3}, boost::numeric::odeint::max_step_checker>, boost::numeric::odeint::dense_output_stepper_tag) (observer=..., dt=0.10000000000000001, end_time=0, start_time=6.4432420282839148e-273, start_state=std::vector of length 1, capacity 1 = {...}, system=..., stepper=...) at stan/lib/stan_math/lib/boost_1.78.0/boost/numeric/odeint/integrate/detail/integrate_times.hpp:149
#6  boost::numeric::odeint::integrate_times<boost::numeric::odeint::dense_output_runge_kutta<boost::numeric::odeint::controlled_runge_kutta<boost::numeric::odeint::runge_kutta_dopri5<std::vector<double, std::allocator<double> >, double, std::vector<double, std::allocator<double> >, double, boost::numeric::odeint::range_algebra, boost::numeric::odeint::default_operations, boost::numeric::odeint::initially_resizer>, boost::numeric::odeint::default_error_checker<double, boost::numeric::odeint::range_algebra, boost::numeric::odeint::default_operations>, boost::numeric::odeint::default_step_adjuster<double, double>, boost::numeric::odeint::initially_resizer, boost::numeric::odeint::explicit_error_stepper_fsal_tag>, boost::numeric::odeint::explicit_controlled_stepper_fsal_tag>, std::reference_wrapper<stan::math::coupled_ode_system<crash_model_namespace::ode_variadic2_functor__, double, std::vector<double, std::allocator<double> > const&> >, std::vector<double, std::allocator<double> >, __gnu_cxx::__normal_iterator<double*, std::vector<double, std::allocator<double> > >, double, stan::math::ode_rk45_tol_impl<crash_model_namespace::ode_variadic2_functor__, Eigen::Matrix<double, -1, 1, 0, -1, 1>, double, double, std::vector<double, std::allocator<double> >, (void*)0>(char const*, crash_model_namespace::ode_variadic2_functor__ const&, Eigen::Matrix<double, -1, 1, 0, -1, 1> const&, double, std::vector<double, std::allocator<double> > const&, double, double, long, std::ostream*, std::vector<double, std::allocator<double> > const&)::{lambda(std::vector<double, std::allocator<double> > const&, double)#3}, boost::numeric::odeint::max_step_checker>(boost::numeric::odeint::dense_output_runge_kutta<boost::numeric::odeint::controlled_runge_kutta<boost::numeric::odeint::runge_kutta_dopri5<std::vector<double, std::allocator<double> >, double, std::vector<double, std::allocator<double> >, double, boost::numeric::odeint::range_algebra, boost::numeric::odeint::default_operations, boost::numeric::odeint::initially_resizer>, boost::numeric::odeint::default_error_checker<double, boost::numeric::odeint::range_algebra, boost::numeric::odeint::default_operations>, boost::numeric::odeint::default_step_adjuster<double, double>, boost::numeric::odeint::initially_resizer, boost::numeric::odeint::explicit_error_stepper_fsal_tag>, boost::numeric::odeint::explicit_controlled_stepper_fsal_tag>, std::reference_wrapper<stan::math::coupled_ode_system<crash_model_namespace::ode_variadic2_functor__, double, std::vector<double, std::allocator<double> > const&> >, std::vector<double, std::allocator<double> >&, __gnu_cxx::__normal_iterator<double*, std::vector<double, std::allocator<double> > >, __gnu_cxx::__normal_iterator<double*, std::vector<double, std::allocator<double> > >, double, stan::math::ode_rk45_tol_impl<crash_model_namespace::ode_variadic2_functor__, Eigen::Matrix<double, -1, 1, 0, -1, 1>, double, double, std::vector<double, std::allocator<double> >, (void*)0>(char const*, crash_model_namespace::ode_variadic2_functor__ const&, Eigen::Matrix<double, -1, 1, 0, -1, 1> const&, double, std::vector<double, std::allocator<double> > const&, double, double, long, std::ostream*, std::vector<double, std::allocator<double> > const&)::{lambda(std::vector<double, std::allocator<double> > const&, double)#3}, boost::numeric::odeint::max_step_checker) (checker=..., observer=..., dt=0.10000000000000001, times_end=..., times_start=..., start_state=std::vector of length 1, capacity 1 = {...}, system=..., stepper=...) at stan/lib/stan_math/lib/boost_1.78.0/boost/numeric/odeint/integrate/integrate_times.hpp:53
#7  stan::math::ode_rk45_tol_impl<crash_model_namespace::ode_variadic2_functor__, Eigen::Matrix<double, -1, 1, 0, -1, 1>, double, double, std::vector<double, std::allocator<double> >, (void*)0> (function_name=function_name@entry=0x5555556948d9 "ode_rk45", f=..., y0_arg=Eigen::Matrix<double,1,1,ColMajor> (data ptr: 0x5555557203f0) = {...}, t0=<optimized out>, t0@entry=0, ts=std::vector of length 1, capacity 1 = {...}, relative_tolerance=<optimized out>, relative_tolerance@entry=9.9999999999999995e-07, absolute_tolerance=<optimized out>, absolute_tolerance@entry=9.9999999999999995e-07, max_num_steps=<optimized out>, max_num_steps@entry=1000000, msgs=0x7fffffffac30) at stan/lib/stan_math/stan/math/prim/functor/ode_rk45.hpp:145
#8  0x000055555558cbec in stan::math::ode_rk45<crash_model_namespace::ode_variadic2_functor__, Eigen::Matrix<double, -1, 1, 0, -1, 1>, double, double, std::vector<double, std::allocator<double> >, (void*)0> (msgs=0x7fffffffac30, ts=std::vector of length 1, capacity 1 = {...}, t0=0, y0=Eigen::Matrix<double,1,1,ColMajor> (data ptr: 0x5555557203f0) = {...}, f=...) at stan/lib/stan_math/stan/math/prim/functor/ode_rk45.hpp:254
#9  crash_model_namespace::u<double, double, double, double, (void*)0> (z=<optimized out>, faketheta=std::vector of length 5, capacity 5 = {...}, pstream__=0x7fffffffac30, x_i=..., x_r=..., xc=<optimized out>) at ../../ml/stanc3/crash.hpp:171
#10 0x00005555555a6c61 in crash_model_namespace::u_functor__::operator()<double, double, double, double, (void*)0> (pstream__=<optimized out>, x_i=..., x_r=..., faketheta=..., xc=<synthetic pointer>: <optimized out>, z=@0x7fffffffa5f0: 6.4432420282839148e-273, this=<optimized out>) at ../../ml/stanc3/crash.hpp:205
#11 integrate_1d_adapter<crash_model_namespace::u_functor__>::operator()<double, double, double> (x_i=..., x_r=..., theta=..., msgs=<optimized out>, xc=<synthetic pointer>: <optimized out>, x=@0x7fffffffa5f0: 6.4432420282839148e-273, this=<optimized out>) at stan/lib/stan_math/stan/math/prim/functor/integrate_1d_adapter.hpp:24
#12 stan::math::integrate_1d_impl<integrate_1d_adapter<crash_model_namespace::u_functor__>, std::vector<double, std::allocator<double> >, std::vector<double, std::allocator<double> >, std::vector<int, std::allocator<int> >, (void*)0>(integrate_1d_adapter<crash_model_namespace::u_functor__> const&, double, double, double, std::ostream*, std::vector<double, std::allocator<double> > const&, std::vector<double, std::allocator<double> > const&, std::vector<int, std::allocator<int> > const&)::{lambda(auto:1 const&, auto:2 const&)#1}::operator()<double, double>(double const&, double const&) const (xc=<synthetic pointer>: <optimized out>, x=@0x7fffffffa5f0: 6.4432420282839148e-273, __closure=<optimized out>) at stan/lib/stan_math/stan/math/prim/functor/integrate_1d.hpp:187
#13 stan::math::integrate<stan::math::integrate_1d_impl<integrate_1d_adapter<crash_model_namespace::u_functor__>, std::vector<double, std::allocator<double> >, std::vector<double, std::allocator<double> >, std::vector<int, std::allocator<int> >, (void*)0>(integrate_1d_adapter<crash_model_namespace::u_functor__> const&, double, double, double, std::ostream*, std::vector<double, std::allocator<double> > const&, std::vector<double, std::allocator<double> > const&, std::vector<int, std::allocator<int> > const&)::{lambda(auto:1 const&, auto:2 const&)#1}>(stan::math::integrate_1d_impl<integrate_1d_adapter<crash_model_namespace::u_functor__>, std::vector<double, std::allocator<double> >, std::vector<double, std::allocator<double> >, std::vector<int, std::allocator<int> >, (void*)0>(integrate_1d_adapter<crash_model_namespace::u_functor__> const&, double, double, double, std::ostream*, std::vector<double, std::allocator<double> > const&, std::vector<double, std::allocator<double> > const&, std::vector<int, std::allocator<int> > const&)::{lambda(auto:1 const&, auto:2 const&)#1} const&, double, double, double)::{lambda(double, double)#4}::operator()(double, double) const (__closure=<optimized out>, __closure=<optimized out>, xc=<optimized out>, x=<optimized out>) at stan/lib/stan_math/stan/math/prim/functor/integrate_1d.hpp:141
#14 _ZZN5boost4math10quadrature9tanh_sinhIdNS0_8policies6policyINS3_14default_policyES5_S5_S5_S5_S5_S5_S5_S5_S5_S5_S5_S5_EEE9integrateIZN4stan4math9integrateIZNSA_17integrate_1d_implI20integrate_1d_adapterIN21crash_model_namespace11u_functor__EEJSt6vectorIdSaIdEESJ_SH_IiSaIiEEELPv0EEEdRKT_dddPSoDpRKT0_EUlSP_RKT0_E_EEdSP_dddEUlddE2_EEKDTclcl7declvalISN_EEcl7declvalIdEEcl7declvalIdEEEESN_dddPdS12_PmENKUlddE_clEdd (z=<optimized out>, zc=<optimized out>, __closure=<optimized out>) at stan/lib/stan_math/lib/boost_1.78.0/boost/math/quadrature/tanh_sinh.hpp:244
#15 _ZZN5boost4math10quadrature9tanh_sinhIdNS0_8policies6policyINS3_14default_policyES5_S5_S5_S5_S5_S5_S5_S5_S5_S5_S5_S5_EEE9integrateIZN4stan4math9integrateIZNSA_17integrate_1d_implI20integrate_1d_adapterIN21crash_model_namespace11u_functor__EEJSt6vectorIdSaIdEESJ_SH_IiSaIiEEELPv0EEEdRKT_dddPSoDpRKT0_EUlSP_RKT0_E_EEdSP_dddEUlddE2_EEKDTclcl7declvalISN_EEcl7declvalIdEEcl7declvalIdEEEESN_dddPdS12_PmENKUlddE_clEdd (__closure=<optimized out>, z=<optimized out>, zc=<optimized out>) at stan/lib/stan_math/lib/boost_1.78.0/boost/math/quadrature/tanh_sinh.hpp:241
#16 0x00005555555a6e1b in _ZNK5boost4math10quadrature6detail16tanh_sinh_detailIdNS0_8policies6policyINS4_14default_policyES6_S6_S6_S6_S6_S6_S6_S6_S6_S6_S6_S6_EEE9integrateIZNS1_9tanh_sinhIdS7_E9integrateIZN4stan4math9integrateIZNSE_17integrate_1d_implI20integrate_1d_adapterIN21crash_model_namespace11u_functor__EEJSt6vectorIdSaIdEESN_SL_IiSaIiEEELPv0EEEdRKT_dddPSoDpRKT0_EUlST_RKT0_E_EEdST_dddEUlddE2_EEKDTclcl7declvalISR_EEcl7declvalIdEEcl7declvalIdEEEESR_dddPdS16_PmEUlddE_EES14_SR_S16_S16_PKcdddS17_ (this=0x55555571fa70, f=..., error=error@entry=0x7fffffffa7a8, L1=L1@entry=0x7fffffffa7b0, function=0x55555569498c "tanh_sinh<%1%>::integrate", left_min_complement=left_min_complement@entry=8.9002954340288055e-308, right_min_complement=right_min_complement@entry=8.9002954340288055e-308, tolerance=tolerance@entry=1.4901161193847656e-08, levels=levels@entry=0x7fffffffa7b8) at stan/lib/stan_math/lib/boost_1.78.0/boost/math/quadrature/detail/tanh_sinh_detail.hpp:238
#17 0x000055555558d883 in _ZN5boost4math10quadrature9tanh_sinhIdNS0_8policies6policyINS3_14default_policyES5_S5_S5_S5_S5_S5_S5_S5_S5_S5_S5_S5_EEE9integrateIZN4stan4math9integrateIZNSA_17integrate_1d_implI20integrate_1d_adapterIN21crash_model_namespace11u_functor__EEJSt6vectorIdSaIdEESJ_SH_IiSaIiEEELPv0EEEdRKT_dddPSoDpRKT0_EUlSP_RKT0_E_EEdSP_dddEUlddE2_EEKDTclcl7declvalISN_EEcl7declvalIdEEcl7declvalIdEEEESN_dddPdS12_Pm (levels=0x7fffffffa7b8, L1=0x7fffffffa7b0, error=0x7fffffffa7a8, tolerance=1.4901161193847656e-08, b=<optimized out>, a=<optimized out>, f=..., this=0x7fffffffa7e0) at stan/lib/stan_math/lib/boost_1.78.0/boost/math/quadrature/tanh_sinh.hpp:251
#18 stan::math::integrate<stan::math::integrate_1d_impl<integrate_1d_adapter<crash_model_namespace::u_functor__>, std::vector<double, std::allocator<double> >, std::vector<double, std::allocator<double> >, std::vector<int, std::allocator<int> >, (void*)0>(integrate_1d_adapter<crash_model_namespace::u_functor__> const&, double, double, double, std::ostream*, std::vector<double, std::allocator<double> > const&, std::vector<double, std::allocator<double> > const&, std::vector<int, std::allocator<int> > const&)::{lambda(auto:1 const&, auto:2 const&)#1}>(stan::math::integrate_1d_impl<integrate_1d_adapter<crash_model_namespace::u_functor__>, std::vector<double, std::allocator<double> >, std::vector<double, std::allocator<double> >, std::vector<int, std::allocator<int> >, (void*)0>(integrate_1d_adapter<crash_model_namespace::u_functor__> const&, double, double, double, std::ostream*, std::vector<double, std::allocator<double> > const&, std::vector<double, std::allocator<double> > const&, std::vector<int, std::allocator<int> > const&)::{lambda(auto:1 const&, auto:2 const&)#1} const&, double, double, double) (relative_tolerance=1.4901161193847656e-08, b=<optimized out>, a=0, f=...) at stan/lib/stan_math/stan/math/prim/functor/integrate_1d.hpp:152
#19 stan::math::integrate_1d_impl<integrate_1d_adapter<crash_model_namespace::u_functor__>, std::vector<double, std::allocator<double> >, std::vector<double, std::allocator<double> >, std::vector<int, std::allocator<int> >, (void*)0> (f=..., b=<optimized out>, msgs=<optimized out>, relative_tolerance=1.4901161193847656e-08, a=0) at stan/lib/stan_math/stan/math/prim/functor/integrate_1d.hpp:186
#20 0x00005555555d3566 in stan::math::integrate_1d<crash_model_namespace::u_functor__> (relative_tolerance=1.4901161193847656e-08, msgs=0x7fffffffac30, x_i=std::vector of length 1, capacity 1 = {...}, x_r=std::vector of length 1, capacity 1 = {...}, theta=std::vector of length 5, capacity 5 = {...}, b=<optimized out>, a=0, f=...) at stan/lib/stan_math/stan/math/prim/functor/integrate_1d.hpp:245
#21 crash_model_namespace::crash_model::log_prob_impl<false, true, std::vector<double, std::allocator<double> >, std::vector<int, std::allocator<int> >, (void*)0, (void*)0> (this=0x555555720d20, params_r__=..., params_i__=..., pstream__=0x7fffffffac30) at ../../ml/stanc3/crash.hpp:381
#22 0x000055555564b801 in stan::model::model_base::log_prob<false, true, double> (msgs=0x7fffffffac30, params_i=std::vector of length 0, capacity 0, params_r=std::vector of length 4, capacity 4 = {...}, this=0x555555720d20) at stan/src/stan/model/model_base.hpp:549
#23 stan::services::util::initialize<true, stan::model::model_base, stan::io::var_context, boost::random::additive_combine_engine<boost::random::linear_congruential_engine<unsigned int, 40014u, 0u, 2147483563u>, boost::random::linear_congruential_engine<unsigned int, 40692u, 0u, 2147483399u> > > (model=..., init=..., rng=..., init_radius=init_radius@entry=2, print_timing=print_timing@entry=true, logger=..., init_writer=...) at stan/src/stan/services/util/initialize.hpp:129
#24 0x000055555566a691 in stan::services::sample::hmc_nuts_diag_e_adapt<stan::model::model_base> (model=..., init=..., init_inv_metric=..., random_seed=random_seed@entry=1770809485, chain=chain@entry=1, init_radius=init_radius@entry=2, num_warmup=num_warmup@entry=1000, num_samples=1000, num_thin=1, save_warmup=false, refresh=100, stepsize=stepsize@entry=1, stepsize_jitter=stepsize_jitter@entry=0, max_depth=10, delta=delta@entry=0.80000000000000004, gamma=gamma@entry=0.050000000000000003, kappa=kappa@entry=0.75, t0=t0@entry=10, init_buffer=75, term_buffer=50, window=25, interrupt=..., logger=..., init_writer=..., sample_writer=warning: RTTI symbol not found for class 'stan::callbacks::unique_stream_writer<std::ostream>'
..., diagnostic_writer=warning: RTTI symbol not found for class 'stan::callbacks::unique_stream_writer<std::ostream>'
...) at stan/src/stan/services/sample/hmc_nuts_diag_e_adapt.hpp:73
#25 0x000055555566d315 in stan::services::sample::hmc_nuts_diag_e_adapt<stan::model::model_base> (diagnostic_writer=warning: RTTI symbol not found for class 'stan::callbacks::unique_stream_writer<std::ostream>'
..., sample_writer=warning: RTTI symbol not found for class 'stan::callbacks::unique_stream_writer<std::ostream>'
..., init_writer=..., logger=..., interrupt=..., window=25, term_buffer=50, init_buffer=75, t0=10, kappa=0.75, gamma=0.050000000000000003, delta=0.80000000000000004, max_depth=10, stepsize_jitter=0, stepsize=1, refresh=100, save_warmup=false, num_thin=1, num_samples=1000, num_warmup=1000, init_radius=2, chain=1, random_seed=1770809485, init=..., model=...) at stan/src/stan/services/sample/hmc_nuts_diag_e_adapt.hpp:151
#26 stan::services::sample::hmc_nuts_diag_e_adapt<stan::model::model_base, std::shared_ptr<stan::io::var_context>, stan::callbacks::writer, stan::callbacks::unique_stream_writer<std::ostream>, stan::callbacks::unique_stream_writer<std::ostream> > (model=..., num_chains=num_chains@entry=1, init=std::vector of length 1, capacity 1 = {...}, random_seed=random_seed@entry=1770809485, init_chain_id=init_chain_id@entry=1, init_radius=init_radius@entry=2, num_warmup=num_warmup@entry=1000, num_samples=1000, num_thin=1, save_warmup=false, refresh=100, stepsize=stepsize@entry=1, stepsize_jitter=stepsize_jitter@entry=0, max_depth=10, delta=delta@entry=0.80000000000000004, gamma=gamma@entry=0.050000000000000003, kappa=kappa@entry=0.75, t0=t0@entry=10, init_buffer=75, term_buffer=50, window=25, interrupt=..., logger=..., init_writer=std::vector of length 1, capacity 1 = {...}, sample_writer=warning: RTTI symbol not found for class 'stan::callbacks::unique_stream_writer<std::ostream>'
std::vector of length 1, capacity 1 = {...}, diagnostic_writer=warning: RTTI symbol not found for class 'stan::callbacks::unique_stream_writer<std::ostream>'
std::vector of length 1, capacity 1 = {...}) at stan/src/stan/services/sample/hmc_nuts_diag_e_adapt.hpp:340
#27 0x00005555555fb04b in cmdstan::command (argc=<optimized out>, argv=<optimized out>) at src/cmdstan/command.hpp:991
#28 0x00005555555881aa in main (argc=<optimized out>, argv=<optimized out>) at src/cmdstan/main.cpp:6

Which suggested to me that the ODE call is accessing memory it shouldn’t somehow. My guess is that you have arrays of different sizes which need to be the same somewhere, I’m still poking around

13

Yes, of course, here’s how I am running this model:

from cmdstanpy import CmdStanModel
import numpy as np

program = "model.stan"

model = CmdStanModel(stan_file=program)

data = {"Rexp": 1.710,
        "Rexp_err": 0.019,
        "wbexp": 0.0224,
        "wbexp_err": 0.0001}

inits = {"h": 0.7,
         "Omega_b": 0.05,
         "Omega_c": 0.25,
         "Omega_r": 0.0001}

fit = model.sample(data=data, show_console=True, chains=1, inits=inits)

Is it possible to change the compilation flags that are sent to gcc to give further information on performing illegal memory access? (I’m very far away from being an expert on the matter, but I do remember there being flags on gcc to assist us with that, during runtime of course).

14

I’m trying to go to a lower level and attempt to debug this model.
If you think this doesn’t belong here, please do let me know if I should open a new thread or even a bug report at the Stan Github repository.

I’ve cloned the cmdstan repository, compiled the model and ran the MCMC with the the previously mentioned dataset.
The result was identical to yours and the ones I had obtained before:

$ ./tmp/model sample data file=tmp/data.json
method = sample (Default)
  sample
    num_samples = 1000 (Default)
    num_warmup = 1000 (Default)
    save_warmup = 0 (Default)
    thin = 1 (Default)
    adapt
      engaged = 1 (Default)
      gamma = 0.050000000000000003 (Default)
      delta = 0.80000000000000004 (Default)
      kappa = 0.75 (Default)
      t0 = 10 (Default)
      init_buffer = 75 (Default)
      term_buffer = 50 (Default)
      window = 25 (Default)
    algorithm = hmc (Default)
      hmc
        engine = nuts (Default)
          nuts
            max_depth = 10 (Default)
        metric = diag_e (Default)
        metric_file =  (Default)
        stepsize = 1 (Default)
        stepsize_jitter = 0 (Default)
    num_chains = 1 (Default)
id = 1 (Default)
data
  file = tmp/data.json
init = 2 (Default)
random
  seed = 2024983781 (Default)
output
  file = output.csv (Default)
  diagnostic_file =  (Default)
  refresh = 100 (Default)
  sig_figs = -1 (Default)
  profile_file = profile.csv (Default)
num_threads = 1 (Default)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: Error in function boost::math::lambert_w0<double>: Expected z >= -e^-1 (-0.367879...) but got -1.2928604982088452. (in 'tmp/model.stan', line 64, column 2 to column 80)
Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: Error in function boost::math::lambert_w0<double>: Expected z >= -e^-1 (-0.367879...) but got -1.2164425942816695. (in 'tmp/model.stan', line 64, column 2 to column 80)
[1]    556839 segmentation fault (core dumped)  ./tmp/model sample data file=tmp/data.json

I had the idea of trying to catch undefined behavior by making use of the -fsanitize flags (I’m using the Section 3.11 of the gcc manual as a reference).
However, I’m not sure where to add said flags in the makefile provided at the top level directory.
Any hints on where to place these flags?

Also, is it possible that it’s my model that has a bug which is not being caught by stanc?
I’ve double check the model, with a focus on the indexing, but I might have missed something.
Is it possible to catch that in advance? (I think it’s more likely that it is my fault rather than Stan’s fault).

15

If you make a file called local in the sub-folder of CmdStan called make/, that file will be included. There is a local.example as a template. I believe all you need to do is add CXXFLAGS += -fsanitize

Yes, most indexing bugs will show up at runtime, rather than compile time. We try to have nice error messages rather than seg faults though, so this will likely result in a bug report/issue on github eventually!

If you’re curious, there is an idea called Rice’s theorem which states “all non-trivial semantic properties of programs are undecidable”. In short, even simple questions like “will this program index out of bounds”, while possibly simple in some cases, are essentially impossible to answer in general. This is why stanc doesn’t really try.

16

It seems like it might be a bug specifically in ode_rk45, using ckcr does not seg fault and adams or bdf raise a different error

17

Hurray, it runs!
It seems to be taking too long and is rejecting quite a few samples.
Obviously that’s just my model probably being misspecified, I will have to look into it myself.

As a follow up on your previous post, injecting the -fsanitize flag in local.example and compiling the model with ode_rk45 didn’t do anything besides printing the segmentation fault.
Maybe I’m not doing it correctly,
Regardless, it seems that you already found the source of the bug to be somewhere in ode_rk45, I don’t think I would be able to do much more than that.

Should I open up a bug report on the Stan’s math library Github page?
If so, should I try to do a smaller MWE of the segmentation fault?

So that’s how the make utility knows that make stanc is meant to compile the stanc binary, cool!
Out of curiosity, if you don’t mind asking, what file is called when I write make <my-model-without-the-.stan>? Is it the top-level makefile?

Also, thanks for the theorem that you have sent me, it will definitely be future read material as soon as I have the time to do so!

1 Like
18

local.example is ignored, you need a file just called local (no extension)

This would be great, yes.

Yes. The top level makefile includes a bunch of others, like the make/local file.

The actual rule which processes make mymodel (where mymodel.stan is a file that exists) is in the file make/program:

This uses a “wildcard target” (%). The way to read this is that a make target with no extension (or .exe on windows, handled by the $(EXE) variable) relies on %.hpp, the same name but ending in .hpp. Right above that, you can see another rule which tells make that %.hpp depends on %.stan

So, if you say make mymodel, make will first go “Okay, great, I need to make mymodel.hpp”. It will then say “Oh, that means I need mymodel.stan”

The makefiles are difficult to read because of makes’ weird syntax and the amount of work we need to do to make them work on all platforms, but the basic idea is fundamentally nothing crazy

19

Thank you for you explanation, once again, this time on makefiles!

Compiling the program with -fsanitize=undefined and -fsanitize=address reveals the following output:

method = sample (Default)
  sample
    num_samples = 1000 (Default)
    num_warmup = 1000 (Default)
    save_warmup = 0 (Default)
    thin = 1 (Default)
    adapt
      engaged = 1 (Default)
      gamma = 0.050000000000000003 (Default)
      delta = 0.80000000000000004 (Default)
      kappa = 0.75 (Default)
      t0 = 10 (Default)
      init_buffer = 75 (Default)
      term_buffer = 50 (Default)
      window = 25 (Default)
    algorithm = hmc (Default)
      hmc
        engine = nuts (Default)
          nuts
            max_depth = 10 (Default)
        metric = diag_e (Default)
        metric_file =  (Default)
        stepsize = 1 (Default)
        stepsize_jitter = 0 (Default)
    num_chains = 1 (Default)
id = 1 (Default)
data
  file = tmp/data.json
init = 2 (Default)
random
  seed = 2126629878 (Default)
output
  file = output.csv (Default)
  diagnostic_file =  (Default)
  refresh = 100 (Default)
  sig_figs = -1 (Default)
  profile_file = profile.csv (Default)
num_threads = 1 (Default)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: Error in function boost::math::lambert_w0<double>: Expected z >= -e^-1 (-0.367879...) but got -1.5587280086931916. (in 'tmp/model.stan', line 64, column 2 to column 80)
Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: Error in function boost::math::lambert_w0<double>: Expected z >= -e^-1 (-0.367879...) but got -1.2707744865573449. (in 'tmp/model.stan', line 64, column 2 to column 80)
Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: Error in function boost::math::lambert_w0<double>: Expected z >= -e^-1 (-0.367879...) but got -1.0979054103223893. (in 'tmp/model.stan', line 64, column 2 to column 80)
Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: Error in function boost::math::lambert_w0<double>: Expected z >= -e^-1 (-0.367879...) but got -2.6601930641890781. (in 'tmp/model.stan', line 64, co

However, has I was making a bug report this in the Stan’s math Github repository, I thought about doing a MWE, just for the sake of making it easy for others to understand the problem.
To my surprise, the MWE compiles and runs properly!
Here’s the MWE in Stan:

functions {
  vector ode(real t, vector y, array[] real theta) {
    real Omega_m = theta[1];
    real Omega_r = theta[2];
    real lambda = theta[3];
    real u = y[1];

    real f1 = exp(-lambda*u^2) / (lambda*u*(u^(-2)-2*lambda) - u^(-3));
    real f2 = 1.5*(Omega_m)*(1+t)^2 + 2*Omega_r*(1+t)^3;

    vector[1] uderiv;
    uderiv[1] = f1*f2;

    return uderiv;
  }
}

data {
  array[10] real t;
  array[10] real uobs;
  array[10] real error;
}

transformed data {
}

parameters {
  real Omega_m;
  real Omega_r;
  real lambda;
}

transformed parameters {
  // parameter array
  array[3] real theta = {Omega_m, Omega_r, lambda};

  // initial conditions
  real t0 = 0;
  vector[1] init;
  init[1] = 1;

  // obtain u(t) using ODE solver
  array[10] real u;
  array[10] vector[1] sol = ode_rk45(ode, init, t0, t, theta);
  for (i in 1:10) {
    u[i] = sol[i][1];
  }
}

model {
  // priors
  Omega_m ~ normal(0.3, 10);
  Omega_r ~ normal(0, 10);
  lambda ~ normal(1, 10);

  // likelihood
  u ~ normal(uobs, error);
}

generated quantities {
}

The changes I made were defining a new parameter \Omega_m which is the sum of two other parameters \Omega_b + \Omega_c, I got rid of everything that wasn’t related to the ODE, which meant getting rid of the parameter h and the integral I was computing before.

And here’s how it’s being run:

# imports
from random import gauss, uniform
import numpy as np
#import arviz as az
#import pandas
import stan

# model file
with open("tmp/mwe.stan", "r") as file:
    program = file.read()

# data file
data = {"t": np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10]),
        "uobs": np.array([0.58827315, 0.34228862, 0.22566855, 0.16229165, 0.1237196, 0.09827532, 0.08047538, 0.06745748, 0.05760086, 0.04992747]),
        "error": np.array([0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001])}

# configure number of chains/samples/warmup
chains = 4
samples = 1000
warmup = 500

# initial conditions for each chain
init = []
for i in range(0, chains):
    init.append({})
    init[i]["Omega_m"] = float(0.3)
    init[i]["Omega_r"] = float(0.0001)
    init[i]["lambda"] = float(0.4)

# run!
posterior = stan.build(program, data=data)
fit = posterior.sample(num_chains=chains, num_samples=samples, num_warmup=warmup, init=init, save_warmup=True)

# do stuff after the program runs
from IPython import embed
embed()

I created the dummy dataset by solving the differential equation with scipy.integrate.odeint and the parameter values of \Omega_m = 0.3, \Omega_r = 0.001 and \lambda = 0.4.
I should mention that I’m able to recover the values of the parameters I have used in creating the mock dataset, meaning that the ODE itself should be well-coded.

So, it seems that this quest isn’t over… there is something in my model which is causing, somehow, a reference to a null pointer…
Could it be that the integral was forcing the ODE solver to go somewhere it didn’t like and, by changing to a different ODE solver, it would go over that point and is somehow safe?
This doesn’t feel right, because if it were something like that, then most likely an infinity would should up or something like that, and Stan would reject that sample, instead of crashing with a segmentation fault due to a reference to a null pointer.

20

Okay I think my idea was somehow correct.
By adding an integration routine of top of the previous model then the program crashes when running it with ode_rk45 but not with ode_ckrk. So it seems that the combo integrate_1d + ode_rk45 is the cause of this null pointer being referenced.
However, even though it runs with ode_ckrk, it takes forever to actually get something out of it, as a lot of samples seem to be rejected due to the error estimate of the integral exceeding the given relative tolerance.

Given that this does seem like a bug and not me somehow performing illegal operations on Stan which are then caught further down at runtime, I’ve decided to open a bug report in Stan’s math library Github page, with the following url: #2848.

1 Like