Open thk686 opened 7 years ago
Which version of odeint were you using when it was working, and which version did you upgrade to that caused this behavior?
Unfortunately, I don't track that very carefully. I get Boost from: https://github.com/eddelbuettel/bh
I will see if I can think of a way to get the Boost version when building my package.
I have this bit of code:
// Copyright Timothy H. Keitt 2015 // See license for odeintr package
include
// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::depends(BH)]]
include "boost/numeric/odeint.hpp"
namespace odeint = boost::numeric::odeint;
include "boost/numeric/ublas/vector.hpp"
include "boost/numeric/ublas/matrix.hpp"
namespace ublas = boost::numeric::ublas;
;
namespace odeintr { typedef ublas::vector sys_vec;
typedef ublas::matrix sys_mat;
static const std::size_t N = 3;
typedef sys_vec state_type;
static state_type state(N);
typedef odeint::rosenbrock4 stepper_type;
static const double atol = 1e-06, rtol = 1e-06;
static auto stepper = odeint::make_dense_output(atol, rtol);
typedef std::vector vec_type;
static std::vector rec_x(N);
static vec_type rec_t;
const double alpha = 0.04, beta = 10000, gamma = 3e+07;;
struct stiff_system { void operator()(const sys_vec &x, sys_vec &dxdt , double t) {
dxdt[0] = -alpha x[0] + beta x[1] x[2]; dxdt[1] = alpha x[0] - beta x[1] x[2] - gamma x[1] x[1]; dxdt[2] = gamma x[1] x[1]; ; } };
struct stiff_system_jacobi { void operator()(const sys_vec &x, sys_mat &J, const double t, sys_vec &dfdt) { J(0, 0) = -alpha; J(0, 1) = beta x[2]; J(0, 2) = beta x[1]; J(1, 0) = alpha; J(1, 1) = -(beta x[2] + (gamma x[1] + gamma x[1])); J(1, 2) = -(beta x[1]); J(2, 0) = 0; J(2, 1) = gamma x[1] + gamma x[1]; J(2, 2) = 0; dfdt[0] = 0.0; dfdt[1] = 0.0; dfdt[2] = 0.0;; } };
auto sys = std::make_pair(stiff_system(), stiff_system_jacobi());
static void obs(const state_type x, const double t) { for (int i = 0; i != N; ++i) rec_x[i].push_back(x[i]); rec_t.push_back(t); }
}; // namespace odeintr
static void reserve(odeintr::vec_type::size_type n) { odeintr::rec_t.reserve(n); for (auto &i : odeintr::rec_x) i.reserve(n); }
// [[Rcpp::export]] Rcpp::List robertson_get_output() { Rcpp::List out; out("Time") = Rcpp::wrap(odeintr::rec_t); for (int i = 0; i != odeintr::N; ++i) { auto cnam = std::string("X") + std::to_string(i + 1); out(cnam) = Rcpp::wrap(odeintr::rec_x[i]); } out.attr("class") = "data.frame"; int rows_out = odeintr::rec_t.size(); auto rn = Rcpp::IntegerVector::create(NA_INTEGER, -rows_out); out.attr("row.names") = rn; return out; };
// [[Rcpp::export]] void robertson_set_state(Rcpp::NumericVector new_state) { if (new_state.size() != odeintr::N) Rcpp::stop("Invalid initial state"); std::copy(new_state.begin(), new_state.end(), odeintr::state.begin()); }
// [[Rcpp::export]] std::vector
robertson_get_state()
{
return std::vector(odeintr::state.begin(), odeintr::state.end());
}
// [[Rcpp::export]] void robertson_reset_observer() { for (auto &i : odeintr::rec_x) i.resize(0); odeintr::rec_t.resize(0);
}
// [[Rcpp::export]] Rcpp::List robertson_adap(Rcpp::NumericVector init, double duration, double step_size = 1.0, double start = 0.0) { robertson_set_state(init); robertson_reset_observer(); reserve(duration / step_size); odeint::integrate_adaptive(odeintr::stepper, odeintr::sys, odeintr::state, start, start + duration, step_size, odeintr::obs); return robertson_get_output(); }
// [[Rcpp::export]] Rcpp::List robertson_at(Rcpp::NumericVector init, std::vector times,
double step_size = 1.0,
double start = 0.0)
{
robertson_set_state(init);
robertson_reset_observer(); reserve(times.size());
odeint::integrate_const(odeintr::stepper, odeintr::sys, odeintr::state,
start, times[0], step_size);
odeint::integrate_times(odeintr::stepper, odeintr::sys, odeintr::state,
times.begin(), times.end(), step_size, odeintr::obs);
return robertson_get_output();
}
// [[Rcpp::export]] Rcpp::List robertson_continue_at(std::vector times, double step_size = 1.0)
{
double start = odeintr::rec_t.back();
robertson_reset_observer(); reserve(odeintr::rec_t.size() + times.size());
odeint::integrate_const(odeintr::stepper, odeintr::sys, odeintr::state,
start, times[0], step_size);
odeint::integrate_times(odeintr::stepper, odeintr::sys, odeintr::state,
times.begin(), times.end(), step_size, odeintr::obs);
return robertson_get_output();
}
// [[Rcpp::export]] Rcpp::List robertson(Rcpp::NumericVector init, double duration, double step_size = 1.0, double start = 0.0) { robertson_set_state(init); robertson_reset_observer(); reserve(duration / step_size); odeint::integrate_const(odeintr::stepper, odeintr::sys, odeintr::state, start, start + duration, step_size, odeintr::obs); return robertson_get_output(); }
// [[Rcpp::export]] std::vector
robertson_no_record(Rcpp::NumericVector init,
double duration,
double step_size = 1.0,
double start = 0.0)
{
robertson_set_state(init);
odeint::integrate_adaptive(odeintr::stepper, odeintr::sys, odeintr::state,
start, start + duration, step_size);
return robertson_get_state();
}
;
This code was generated as an example in the odeintr package as follows:
Robertson chemical kinetics problem
Robertson = ' dxdt[0] = -alpha x[0] + beta x[1] x[2]; dxdt[1] = alpha x[0] - beta x[1] x[2] - gamma x[1] x[1]; dxdt[2] = gamma x[1] x[1]; ' # Robertson pars = c(alpha = 0.04, beta = 1e4, gamma = 3e7) init.cond = c(1, 0, 0) compile_implicit("robertson", Robertson, pars, TRUE) at = 10 ^ seq(-5, 5, len = 400) x = robertson_at(init.cond, at) # this never terminates
Something changed causing the integration to never complete (or perhaps run 1000's of times more slowly).