Merge pull request #152 from mrc-ide/mrc-6351

Support events that happen on times that are exactly hit

Author: richfitz

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: dust2
 Title: Next Generation dust
-Version: 0.3.20
+Version: 0.3.21
 Authors@R: c(person("Rich", "FitzJohn", role = c("aut", "cre"),
                     email = "[email protected]"),
              person("Imperial College of Science, Technology and Medicine",

inst/include/dust2/continuous/solver.hpp

@@ -329,9 +329,17 @@ class solver {
   real_type apply_events(real_type t0, real_type h, const real_type* y,
                          const events_type<real_type>& events,
                          ode::internals<real_type>& internals) {
-    size_t idx_first = events.size();
     real_type t1 = t0 + h;
-    real_type sign = 0;
+
+    // It might be worth saving this storage space in the solver, but
+    // I doubt it matters in pratice.  We need to save all the found
+    // events and their signs, even though probably only one will be
+    // found, because of the possibility that multiple events could
+    // happen at the same time (e.g., two time-scheduled events or two
+    // functions that happen to hit their roots at the same time).
+    std::vector<bool> found(events.size());
+    std::vector<real_type> sign(events.size());
+    bool found_any = false;
 
     for (size_t idx_event = 0; idx_event < events.size(); ++idx_event) {
       const auto& e = events[idx_event];
@@ -350,20 +358,44 @@ class solver {
         // interpolation is expected to be quite fast and accurate.
         constexpr real_type eps = 1e-6;
         constexpr size_t steps = 100;
-        auto root = lostturnip::find_result<real_type>(fn, t0, t1, eps, steps);
-        idx_first = idx_event;
-        t1 = root.x;
-        sign = f_t0 < 0 ? 1 : -1;
+        t1 = lostturnip::find_result<real_type>(fn, t0, t1, eps, steps).x;
+        // Currently untested - in the case where we have two roots
+        // that would have been crossed in this time window, the one
+        // we are currently considering happens first so pre-empts the
+        // previously found events.
+        if (found_any) {
+          std::fill(found.begin(), found.end(), false);
+        }
+      } else if (!(f_t1 == 0 && f_t0 != 0)) {
+        // Consider the case where jump to a root *exactly* at t1;
+        // this happens in coincident roots and with roots that are
+        // based in time, and which we arrange for the solver to stop
+        // at (e.g., while using simulate()).
+        //
+        // This test is the inverse of this though, because in the
+        // case where we *don't* get an exact root we should skip the
+        // bookkeeping below and try the next event.
+        continue;
       }
-      if (idx_first < events.size()) {
-        internals.last.interpolate(t1, y_next_.data());
-        events[idx_first].action(t1, sign, y_next_.data());
-        // We need to modify the history here so that search will find
-        // the right point.
-        internals.last.t1 = t1;
-        internals.events.push_back({t1, idx_first, sign});
+      sign[idx_event] = f_t0 < 0 ? 1 : -1;
+      found[idx_event] = true;
+      found_any = true;
+    }
+
+    // If we found at least one event, then reset the solver state
+    // back to the point of the event and apply all the events in
+    // turn.
+    if (found_any) {
+      internals.last.interpolate(t1, y_next_.data());
+      for (size_t idx_event = 0; idx_event < events.size(); ++idx_event) {
+        if (found[idx_event]) {
+          events[idx_event].action(t1, sign[idx_event], y_next_.data());
+          internals.events.push_back({t1, idx_event, sign[idx_event]});
+        }
       }
+      internals.last.t1 = t1;
     }
+
     return t1;
   }
 

tests/testthat/examples/event-time.cpp

@@ -0,0 +1,73 @@
+#include <dust2/common.hpp>
+
+// [[dust2::class(change)]]
+// [[dust2::time_type(continuous)]]
+// [[dust2::parameter(r, rank = 0)]]
+// [[dust2::parameter(n, rank = 0)]]
+// [[dust2::parameter(t_change, rank = 1)]]
+// [[dust2::parameter(delta, rank = 1)]]
+class change {
+public:
+  change() = delete;
+
+  using real_type = double;
+
+  struct shared_state {
+    real_type r;
+    std::vector<real_type> t_change;
+    std::vector<real_type> delta;
+  };
+
+  struct internal_state {};
+
+  using rng_state_type = monty::random::generator<real_type>;
+
+  static dust2::packing packing_state(const shared_state& shared) {
+    return dust2::packing{{"y", {}}};
+  }
+
+  static void initial(real_type time,
+                      const shared_state& shared,
+                      internal_state& internal,
+                      rng_state_type& rng_state,
+                      real_type * state) {
+    state[0] = 0;
+  }
+
+  static void rhs(real_type time,
+                  const real_type * state,
+                  const shared_state& shared,
+                  internal_state& internal,
+                  real_type * state_deriv) {
+    state_deriv[0] = shared.r;
+  }
+
+  static shared_state build_shared(cpp11::list pars) {
+    const real_type r1 = dust2::r::read_real(pars, "r");
+    const size_t n = dust2::r::read_int(pars, "n");
+    std::vector<real_type> t_change(n);
+    std::vector<real_type> delta(n);
+    dust2::r::read_real_vector(pars, n, t_change.data(), "t_change", true);
+    dust2::r::read_real_vector(pars, n, delta.data(), "delta", false);
+    return shared_state{r1, t_change, delta};
+  }
+
+  static void update_shared(cpp11::list pars, shared_state& shared) {
+  }
+
+  static auto events(const shared_state& shared, internal_state& internal) {
+    const auto n = shared.t_change.size();
+    dust2::ode::events_type<real_type> events;
+    events.reserve(n);
+    for (size_t i = 0; i < n; ++i) {
+      auto test = [&, i](const double t, const real_type* y) {
+        return t - shared.t_change[i];
+      };
+      auto action = [&, i](const double t, const double sign, double* y) {
+        y[0] += shared.delta[i];
+      };
+      events.push_back(dust2::ode::event<real_type>({}, test, action));
+    }
+    return events;
+  }
+};

tests/testthat/test-zzz-events.R

@@ -28,3 +28,52 @@ test_that("can run system with roots and events", {
   ## Overall solution:
   expect_equal(y[1, ], cmp$y, tolerance = 1e-6)
 })
+
+
+test_that("can run events with events in time", {
+  gen <- dust_compile("examples/event-time.cpp", quiet = FALSE, debug = TRUE)
+
+  # A mix of times that we will hit exactly and bracket
+  pars <- list(r = 0.2, n = 3, t_change = c(2, 5.1234, 7), delta = rnorm(3))
+  control <- dust_ode_control(
+    debug_record_step_times = TRUE,
+    save_history = TRUE
+  )
+  sys <- dust_system_create(gen, pars, ode_control = control)
+
+  t <- seq(0, 10, length.out = 101)
+  y <- dust_system_simulate(sys, t)
+
+  info <- dust_system_internals(sys, include_history = TRUE)
+  expect_equal(
+    info$events[[1]],
+    data_frame(time = pars$t_change, index = 1:3, sign = 1)
+  )
+  expect_equal(
+    drop(y),
+    t * 0.2 + c(0, cumsum(pars$delta))[findInterval(t, pars$t_change) + 1])
+})
+
+
+test_that("can cope with coincident events", {
+  gen <- dust_compile("examples/event-time.cpp", quiet = FALSE, debug = TRUE)
+
+  pars <- list(r = 0.2, n = 3, t_change = c(2, 2, 3), delta = c(1, 3, 5))
+  control <- dust_ode_control(
+    debug_record_step_times = TRUE,
+    save_history = TRUE
+  )
+  sys <- dust_system_create(gen, pars, ode_control = control)
+
+  t <- seq(0, 10, length.out = 101)
+  y <- dust_system_simulate(sys, t)
+
+  info <- dust_system_internals(sys, include_history = TRUE)
+  expect_equal(
+    info$events[[1]],
+    data_frame(time = pars$t_change, index = 1:3, sign = 1)
+  )
+  expect_equal(
+    drop(y),
+    t * 0.2 + c(0, cumsum(pars$delta))[findInterval(t, pars$t_change) + 1])
+})