Initial post because DAE didn’t accept time that depends on the parameters;
This is intended behavior and an exception is now being caught by stanc instead of throwing an error at the level of GCC;
Was provided with a different solution to my problem by turning this problem into solving an ODE instead;
Finding a segmentation fault when using integrate_1d on top of ode_rk45, which have lead me to open the bug report #2848 on Stan’s math Github repository.
The discussion regarding the problems with ode_rk45 or ode_ckrk should be continued in the Github issue report.
However, if you have anything else that could help me bypass this issue and implement my model (where said model is described in post #10 with dataset given in post #13 in this thread) then by all means reply here! Any help will be greatly appreciated!
Given that the other route failed I am now trying to use this approach instead, where I simply re-scale the time and pass z_* as an additional argument to the function.
For the sake of simplicity, I’m just trying to solve the first equation (the algebraic one), i.e.:
where \Omega_m and \Omega_r are parameters, \lambda is computed from the previous parameters using a Lambert function and z plays the role of time.
I’m trying to solve the equation and obtain E(z).
My Stan model file looks like the following:
/* block for user defined functions */
functions {
// system of differential equations to solve
vector dae_system(real x, vector state, vector deriv, array[] real theta, real lambda) {
real Omega_m = theta[1];
real Omega_r = theta[2];
// rescale by a constant instead of a variable which depends on the parameters to keep it simple...
real z=x*1048;
// lhs and rhs of the modified 1st friedmann eq
real lhs = (state[1]^2 - 2*lambda)*exp(lambda/state[1]^2);
real rhs = (Omega_m)*(1+z)^3 + Omega_r*(1+z)^4;
// compute the residuals
vector[1] res;
res[1] = lhs - rhs;
return res;
}
}
/* block for datasets used */
data {
real Omega_m_exp;
real Omega_m_exp_err;
real Omega_r_exp;
real Omega_r_exp_err;
}
/* block for transformed data */
transformed data {
}
/* block for model parameters */
parameters {
real Omega_m;
real Omega_r;
}
/* block for transformed parameters */
transformed parameters {
// debug
print("(Omega_m, Omega_r) = ", Omega_m, ", ", Omega_r);
// lambda
real lambda = 0.5 + lambert_w0( -(Omega_m + Omega_r)/(2*exp(0.5)) );
// parameter array
array[2] real theta = {Omega_m, Omega_r};
// initial conditions to provide to the DAE
real E_0 = 1;
real derivE_0 = 1/(2*exp(lambda)) * 3*(Omega_m)/(1 - lambda + 2*lambda^2);
// call the differential algebraic equation solver
array[1] vector[1] S;
S = dae(dae_system, [E_0]', [derivE_0]', 0, {1}, theta, lambda);
// dummie variables
real Omega_mdummie = 2*Omega_m;
real Omega_rdummie = 2*Omega_r;
}
/* block for model definition */
model {
// priors
Omega_m ~ normal(0.30, 10);
Omega_r ~ normal(0, 10);
// dummie likelihoods
Omega_mdummie ~ normal(Omega_m_exp, Omega_m_exp_err);
Omega_rdummie ~ normal(Omega_r_exp, Omega_r_exp_err);
}
/* block for generated quantites */
generated quantities {
}
Only that now, for some reason, this program seemingly freezes and stops outputting anything.
Only that in reality it stays running in the background eating an increasing amount of RAM to the point that the kernel has to kill it due to lack of memory available!
Here’s the relevant output of dmesg:
[qui nov 17 11:56:52 2022] Out of memory: Killed process 879715 (dae) total-vm:25180500kB, anon-rss:6074124kB, file-rss:64kB, shmem-rss:0kB, UID:1000 pgtables:19240kB oom_score_adj:0
If I am to compile the program with -fsanitize=undefined and -fsanitize=address I find that there is a memory related error:
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/dae/data-dae-reduced.json
init = 2 (Default)
random
seed = 2359428734 (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)
(Omega_m, Omega_r) = 0.946979, -1.42703
(Omega_m, Omega_r) = 0.946979, -1.42703
Gradient evaluation took 0.020271 seconds
1000 transitions using 10 leapfrog steps per transition would take 202.71 seconds.
Adjust your expectations accordingly!
(Omega_m, Omega_r) = 0.946979, -1.42703
(Omega_m, Omega_r) = -1.29396e+06, 2.85427e+10
Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Exception: Error in function boost::math::lambert_w0<double>: Expected z >= -e^-1 (-0.367879...) but got -8655616723.7307339. (in 'tmp/dae/dae.stan', line 47, column 2 to column 70)
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
(Omega_m, Omega_r) = 0.946979, -1.42703
(Omega_m, Omega_r) = -1.29396e+06, 2.85427e+10
Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Exception: Error in function boost::math::lambert_w0<double>: Expected z >= -e^-1 (-0.367879...) but got -8655616723.3131599. (in 'tmp/dae/dae.stan', line 47, column 2 to column 70)
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
(Omega_m, Omega_r) = 0.946979, -1.42703
(Omega_m, Omega_r) = -323489, 7.13567e+09
Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Exception: Error in function boost::math::lambert_w0<double>: Expected z >= -e^-1 (-0.367879...) but got -2163904180.732583. (in 'tmp/dae/dae.stan', line 47, column 2 to column 70)
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
(Omega_m, Omega_r) = 0.946979, -1.42703
(Omega_m, Omega_r) = -80871.9, 1.78392e+09
Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Exception: Error in function boost::math::lambert_w0<double>: Expected z >= -e^-1 (-0.367879...) but got -540976044.91572511. (in 'tmp/dae/dae.stan', line 47, column 2 to column 70)
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
(Omega_m, Omega_r) = 0.946979, -1.42703
(Omega_m, Omega_r) = -20217, 4.4598e+08
Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Exception: Error in function boost::math::lambert_w0<double>: Expected z >= -e^-1 (-0.367879...) but got -135244011.12626737. (in 'tmp/dae/dae.stan', line 47, column 2 to column 70)
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
(Omega_m, Omega_r) = 0.946979, -1.42703
(Omega_m, Omega_r) = -5053.55, 1.11495e+08
Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Exception: Error in function boost::math::lambert_w0<double>: Expected z >= -e^-1 (-0.367879...) but got -33811002.709765173. (in 'tmp/dae/dae.stan', line 47, column 2 to column 70)
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
(Omega_m, Omega_r) = 0.946979, -1.42703
(Omega_m, Omega_r) = -1262.69, 2.78737e+07
Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Exception: Error in function boost::math::lambert_w0<double>: Expected z >= -e^-1 (-0.367879...) but got -8452750.5516611058. (in 'tmp/dae/dae.stan', line 47, column 2 to column 70)
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
(Omega_m, Omega_r) = 0.946979, -1.42703
(Omega_m, Omega_r) = -314.935, 6.96843e+06
Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Exception: Error in function boost::math::lambert_w0<double>: Expected z >= -e^-1 (-0.367879...) but got -2113187.539645826. (in 'tmp/dae/dae.stan', line 47, column 2 to column 70)
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
(Omega_m, Omega_r) = 0.946979, -1.42703
(Omega_m, Omega_r) = -78.0296, 1.74211e+06
Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Exception: Error in function boost::math::lambert_w0<double>: Expected z >= -e^-1 (-0.367879...) but got -528296.76895409555. (in 'tmp/dae/dae.stan', line 47, column 2 to column 70)
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
(Omega_m, Omega_r) = 0.946979, -1.42703
(Omega_m, Omega_r) = -18.793, 435526
Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Exception: Error in function boost::math::lambert_w0<double>: Expected z >= -e^-1 (-0.367879...) but got -132074.08537642061. (in 'tmp/dae/dae.stan', line 47, column 2 to column 70)
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
(Omega_m, Omega_r) = 0.946979, -1.42703
(Omega_m, Omega_r) = -3.98622, 108880
Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Exception: Error in function boost::math::lambert_w0<double>: Expected z >= -e^-1 (-0.367879...) but got -33018.413037123115. (in 'tmp/dae/dae.stan', line 47, column 2 to column 70)
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
(Omega_m, Omega_r) = 0.946979, -1.42703
(Omega_m, Omega_r) = -0.285567, 27219
Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Exception: Error in function boost::math::lambert_w0<double>: Expected z >= -e^-1 (-0.367879...) but got -8254.4938174598556. (in 'tmp/dae/dae.stan', line 47, column 2 to column 70)
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
(Omega_m, Omega_r) = 0.946979, -1.42703
(Omega_m, Omega_r) = 0.638404, 6803.68
Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Exception: Error in function boost::math::lambert_w0<double>: Expected z >= -e^-1 (-0.367879...) but got -2063.5143226668979. (in 'tmp/dae/dae.stan', line 47, column 2 to column 70)
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
(Omega_m, Omega_r) = 0.946979, -1.42703
(Omega_m, Omega_r) = 0.869905, 1699.85
Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Exception: Error in function boost::math::lambert_w0<double>: Expected z >= -e^-1 (-0.367879...) but got -515.76931708716177. (in 'tmp/dae/dae.stan', line 47, column 2 to column 70)
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
(Omega_m, Omega_r) = 0.946979, -1.42703
(Omega_m, Omega_r) = 0.927653, 423.892
Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Exception: Error in function boost::math::lambert_w0<double>: Expected z >= -e^-1 (-0.367879...) but got -128.83316256465974. (in 'tmp/dae/dae.stan', line 47, column 2 to column 70)
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
(Omega_m, Omega_r) = 0.946979, -1.42703
(Omega_m, Omega_r) = 0.942153, 104.903
Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Exception: Error in function boost::math::lambert_w0<double>: Expected z >= -e^-1 (-0.367879...) but got -32.099118220219594. (in 'tmp/dae/dae.stan', line 47, column 2 to column 70)
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
(Omega_m, Omega_r) = 0.946979, -1.42703
(Omega_m, Omega_r) = 0.945798, 25.1554
Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Exception: Error in function boost::math::lambert_w0<double>: Expected z >= -e^-1 (-0.367879...) but got -7.9155991215617183. (in 'tmp/dae/dae.stan', line 47, column 2 to column 70)
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
(Omega_m, Omega_r) = 0.946979, -1.42703
(Omega_m, Omega_r) = 0.946714, 5.21855
Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Exception: Error in function boost::math::lambert_w0<double>: Expected z >= -e^-1 (-0.367879...) but got -1.8697112286928406. (in 'tmp/dae/dae.stan', line 47, column 2 to column 70)
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
(Omega_m, Omega_r) = 0.946979, -1.42703
(Omega_m, Omega_r) = 0.946911, 0.234364
Iteration: 1 / 2000 [ 0%] (Warmup)
(Omega_m, Omega_r) = 0.946979, -1.42703
(Omega_m, Omega_r) = 0.946904, 0.234379
=================================================================
==880793==ERROR: AddressSanitizer: attempting free on address which was not malloc()-ed: 0x603000d6d210 in thread T0
#0 0x7f78dd6fb672 in __interceptor_free /usr/src/debug/gcc/libsanitizer/asan/asan_malloc_linux.cpp:52
#1 0x5568a6bb7269 in stan::mcmc::base_nuts<stan::model::model_base, stan::mcmc::diag_e_metric, stan::mcmc::expl_leapfrog, 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> > >::transition(stan::mcmc::sample&, stan::callbacks::logger&) (/home/undercover/downloads/cmdstan-2.30.1/tmp/dae/dae+0x662269)
#2 0x5568a6bb7a85 in stan::mcmc::adapt_diag_e_nuts<stan::model::model_base, 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> > >::transition(stan::mcmc::sample&, stan::callbacks::logger&) (/home/undercover/downloads/cmdstan-2.30.1/tmp/dae/dae+0x662a85)
#3 0x5568a6b16427 in void stan::services::util::generate_transitions<stan::model::model_base, 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> > >(stan::mcmc::base_mcmc&, int, int, int, int, int, bool, bool, stan::services::util::mcmc_writer&, stan::mcmc::sample&, stan::model::model_base&, 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> >&, stan::callbacks::interrupt&, stan::callbacks::logger&, unsigned long, unsigned long) [clone .constprop.1] (/home/undercover/downloads/cmdstan-2.30.1/tmp/dae/dae+0x5c1427)
#4 0x5568a6b91fea in int stan::services::sample::hmc_nuts_diag_e_adapt<stan::model::model_base>(stan::model::model_base&, stan::io::var_context const&, stan::io::var_context const&, unsigned int, unsigned int, double, int, int, int, bool, int, double, double, int, double, double, double, double, unsigned int, unsigned int, unsigned int, stan::callbacks::interrupt&, stan::callbacks::logger&, stan::callbacks::writer&, stan::callbacks::writer&, stan::callbacks::writer&) (/home/undercover/downloads/cmdstan-2.30.1/tmp/dae/dae+0x63cfea)
#5 0x5568a6b94344 in int 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> >(stan::model::model_base&, unsigned long, std::vector<std::shared_ptr<stan::io::var_context>, std::allocator<std::shared_ptr<stan::io::var_context> > > const&, unsigned int, unsigned int, double, int, int, int, bool, int, double, double, int, double, double, double, double, unsigned int, unsigned int, unsigned int, stan::callbacks::interrupt&, stan::callbacks::logger&, std::vector<stan::callbacks::writer, std::allocator<stan::callbacks::writer> >&, std::vector<stan::callbacks::unique_stream_writer<std::ostream>, std::allocator<stan::callbacks::unique_stream_writer<std::ostream> > >&, std::vector<stan::callbacks::unique_stream_writer<std::ostream>, std::allocator<stan::callbacks::unique_stream_writer<std::ostream> > >&) (/home/undercover/downloads/cmdstan-2.30.1/tmp/dae/dae+0x63f344)
#6 0x5568a6b32c70 in cmdstan::command(int, char const**) (/home/undercover/downloads/cmdstan-2.30.1/tmp/dae/dae+0x5ddc70)
#7 0x5568a68d5a45 in main (/home/undercover/downloads/cmdstan-2.30.1/tmp/dae/dae+0x380a45)
#8 0x7f78dcabb28f (/usr/lib/libc.so.6+0x2328f)
#9 0x7f78dcabb349 in __libc_start_main (/usr/lib/libc.so.6+0x23349)
#10 0x5568a68d6214 in _start ../sysdeps/x86_64/start.S:115
0x603000d6d210 is located 16 bytes inside of 32-byte region [0x603000d6d200,0x603000d6d220)
allocated by thread T0 here:
#0 0x7f78dd6fca89 in __interceptor_malloc /usr/src/debug/gcc/libsanitizer/asan/asan_malloc_linux.cpp:69
#1 0x5568a6949ff1 in Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >::PlainObjectBase<Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> > >(Eigen::DenseBase<Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> > > const&) (/home/undercover/downloads/cmdstan-2.30.1/tmp/dae/dae+0x3f4ff1)
SUMMARY: AddressSanitizer: bad-free /usr/src/debug/gcc/libsanitizer/asan/asan_malloc_linux.cpp:52 in __interceptor_free
==880793==ABORTING
You can easily see is that this re-scaling does exactly the same as if we were to call de DAE on the time interval from 0 to 1024, instead of calling the DAE from 0 to 1 and re-scale the time inside the function.
It’s interesting that this bug also happens in both situations (with and without re-scaling and adjusting the maximum time accordingly), meaning that this is something to do with the system at hand, and not the DAE (or the re-scaling) itself!
Only that for the several values that I have played around in Mathematica, I always get a proper solution to this system, and I’m fairly confident that this equation is well behaved…
Besides, if there were numerical errors, there shouldn’t be a memory error…