Code that worked now fails

[thk686]

I have this bit of code:

// Copyright Timothy H. Keitt 2015 // See license for odeintr package

#include <Rcpp.h> // [[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<stepper_type>(atol, rtol);

typedef std::vector vec_type; static std::vector<vec_type> 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).

[ThomasThelen]

Which version of odeint were you using when it was working, and which version did you upgrade to that caused this behavior?

[thk686]

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.