Segmentation fault on making use of ode_rk45 with integrate_1d
#2848
Comments
|
Thanks for reporting this. Seems to be an ugly bug! FYI: ode_ckrk uses the CVODES utilities while ode_rk45 is based on odeint. Both should work, of course. |
|
I would like to state that using Could it be that this model is problematic to solve numerically? How would I go about and test that? Regardless of how troublesome the model, it shouldn't give a segmentation fault, right? |
100% |
|
|
|
Looks like I forgot this one...just got the example to run and can confirm that the segfault happens on my system as well. |
|
Alright! After a bit of debugging this I am pretty sure that the problem is in RK45 integrator of odeint. The segfault does even happen if you run this thing without any The problem is caused by the integration for integrate_1d starting at 0 AND the ode integrator is also being asked to start integrating at 0. That will trigger the integrator to evaluate time-points for the function which are super close to 0. So close that the odeint integrator is being put off the rails (which id should NOT). So this slightly modified Stan program runs just fine: functions {
// return u'(t) to use in ODE solver
vector dudt(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;
}
// return u(t) to use in integration routine
real ut(real t, real xc, array[] real theta, array[] real x_r, array[] int x_i) {
// initial conditions
real t0 = 0;
vector[1] init;
init[1] = 1;
// obtain u(t) using ODE solver
real sol = ode_rk45(dudt, init, t0, {t}, theta)[1][1];
return sol;
}
}
data {
array[10] real t;
array[10] real intofuobs;
array[10] real error;
}
transformed data {
// create empty x_r and x_i arrays to provide to integrate_1d given that they are required arguments
array[0] real x_r;
array[0] int x_i;
}
parameters {
real Omega_m;
real Omega_r;
real lambda;
}
transformed parameters {
// parameter array
array[3] real theta = {Omega_m, Omega_r, lambda};
// value of intofu given by the current parameters
array[10] real intofu;
for (i in 1:10) {
// NOTE: integrate from almost 0 here...since odeint will always start at 0 we do not run into issues.
intofu[i] = integrate_1d(ut, 1E-2, t[i], theta, x_r, x_i);
}
}
model {
// priors
Omega_m ~ normal(0.3, 10);
Omega_r ~ normal(0, 10);
lambda ~ normal(1, 10);
// likelihood
intofu ~ normal(intofuobs, error);
}
generated quantities {
}As a next step I would try to compile a simple example which throws over odeint using a small example. I'd guess that the odeint author will fix this quite likely. Not sure if we need this issue here open any longer. At least I am confident that this is nothing wrong in Stan-math here. |
|
Interesting, can I know where is Stan getting the |
|
I solved a much simpler ODE and triggered the same problem. You should turn things into a first oder oder, of course! |
|
ode_rk45 is based on the odeint library part of boost. |
|
Many thanks! |
|
I just saw that there has not been any activity on the odeint GitHub for ~4 years. Not a good sign. So let's keep this open here as we may have to deal with the issue from the Stan-math side. Still... I am very much relieved to see that there is not any severe issues on Stan-math side of things. It always worried me, but we are good here. |
Description
When solving an ODE using
$\left.u^{\prime}=\frac{e^{-\lambda u^2}}{\lambda u\left(u^{-2}-2 \lambda\right)-u^{-3}}\left[\frac{3}{2}\Omega_m(1+z)^2+2 \Omega_r(1+z)^3\right)\right]$ ,
ode_rk45which takes the formwhere$\Omega_m$ , $\Omega_r$ and $\lambda$ are parameters, and then integrating $u(t)$ using
integrate_1d, i.e. computingwhere the vector of observations is the set of values$(t, X, \sigma)$ , yields a segmentation fault.
However, if I make use of
ode_ckrkinstead, it does not give a segmentation fault, and the program runs as expected.Even though it runs as expected it takes a long time due to the fact that the sampler rejects a lot of steps due to the error estimate in the integral being great than the given relative tolerance.
Using
ode_rk45successfullyHere's a Stan model file that makes use of
ode_rk45to solve the previous ODE successfully:Compiling and running this model using
stancreveals no issues whatsoever, and I am able to recover the parameters I have used to generate a mock dataset (given in data.csv).Adding the integration routine to the previous model
If I now add a new function, which returns the value of$u(t)$ , i.e. it solves the ODE for a specific time, so that I can then integrate it using
integrate_1din the transformed parameters block.Doing so the model now looks as follows:
I compiled this model with two additional flags
CXXFLAGS += -fsanitize=undefined -fsanitize=addressand when running this model with a new mock dataset (given in data-int.csv) the following error was caught during runtime:Replacing
ode_rk45withode_rkckBy replacing
ode_rk45withode_rkckthe model runs successfully, making use of the previously mentioned dataset.However, as mentioned before, the model takes a long time to run and a lot of samples are rejected due to the the error estimate in the integral being great than the given relative tolerance.
System information:
Additional information
This bug report was the result of a thread I posted over at the Stan's forums, where the user WardBrian was able to replicate this issue and helped me a lot with the debugging (link).
The text was updated successfully, but these errors were encountered: