Merge pull request #411 from mrc-ide/mrc-4347

mrc-4347: cope with missing data

Author: r-ash

.gitattributes

@@ -0,0 +1,2 @@
+R/cpp11.R linguist-generated=true
+src/cpp11.cpp linguist-generated=true

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: dust
 Title: Iterate Multiple Realisations of Stochastic Models
-Version: 0.14.7
+Version: 0.14.8
 Authors@R: c(person("Rich", "FitzJohn", role = c("aut", "cre"),
                     email = "[email protected]"),
              person("Alex", "Hill", role = "aut"),

R/cpp11.R

@@ -1,4 +1,4 @@
-## Generated by dust (version 0.14.4) - do not edit
+## Generated by dust (version 0.14.8) - do not edit
 logistic <- R6::R6Class(
   "dust",
   cloneable = FALSE,

R/dust.R

@@ -59,8 +59,11 @@ class sir {
 
   real_type compare_data(const real_type * state, const data_type& data,
                          rng_state_type& rng_state) {
-    const real_type incidence_modelled = state[4];
     const real_type incidence_observed = data.incidence;
+    if (std::isnan(data.incidence)) {
+      return 0;
+    }
+    const real_type incidence_modelled = state[4];
     const real_type lambda = incidence_modelled +
       dust::random::exponential(rng_state, shared->exp_noise);
     return dust::density::poisson(incidence_observed, lambda, true);

inst/examples/sir.cpp

@@ -59,10 +59,6 @@ filter(T * obj,
     obj->run(time);
     obj->compare_data(weights, d->second);
 
-    // TODO: we should cope better with the case where all weights
-    // are 0; I think that is the behaviour in the model (or rather
-    // the case where there is no data and so we do not resample)
-    //
     // TODO: we should cope better with the case where one filter
     // has become impossible but others continue, but that's hard!
     auto wi = weights.begin();

inst/include/dust/filter.hpp

@@ -8,6 +8,7 @@
 #include <vector>
 
 #include "dust/random/random.hpp"
+#include "dust/utils.hpp"
 
 namespace dust {
 namespace filter {
@@ -17,23 +18,30 @@ void resample_weight(typename std::vector<real_type>::const_iterator w,
                      size_t n, real_type u, size_t offset,
                      typename std::vector<size_t>::iterator idx) {
   const real_type tot = std::accumulate(w, w + n, static_cast<real_type>(0));
-  real_type ww = 0.0, uu0 = tot * u / n, du = tot / n;
-  size_t j = offset;
-  const size_t end = n + offset;
-  for (size_t i = 0; i < n; ++i) {
-    // We could accumulate uu by adding du at each iteration but that
-    // suffers roundoff error here with floats.
-    const real_type uu = uu0 + i * du;
-    // The second clause (i.e., j - offset < n) should never be hit
-    // but prevents any invalid read if we have pathalogical 'u' that
-    // is within floating point eps of 1
-    while (ww < uu && j < end) {
-      ww += *w;
-      ++w;
-      ++j;
+  if (tot == 0 && dust::utils::all_zero<real_type>(w, w + n)) {
+    for (size_t i = 0; i < n; ++i) {
+      *idx = offset + i;
+      ++idx;
+    }
+  } else {
+    real_type ww = 0.0, uu0 = tot * u / n, du = tot / n;
+    size_t j = offset;
+    const size_t end = n + offset;
+    for (size_t i = 0; i < n; ++i) {
+      // We could accumulate uu by adding du at each iteration but that
+      // suffers roundoff error here with floats.
+      const real_type uu = uu0 + i * du;
+      // The second clause (i.e., j - offset < n) should never be hit
+      // but prevents any invalid read if we have pathalogical 'u' that
+      // is within floating point eps of 1
+      while (ww < uu && j < end) {
+        ww += *w;
+        ++w;
+        ++j;
+      }
+      *idx = j == 0 ? 0 : j - 1;
+      ++idx;
     }
-    *idx = j == 0 ? 0 : j - 1;
-    ++idx;
   }
 }
 

inst/include/dust/filter_tools.hpp

@@ -301,12 +301,16 @@ void find_intervals(const real_type * cum_weights,
     real_type start_val = par_idx > 0 ? cum_weights[par_idx * n_particles_each - 1] : 0;
     real_type normalising_constant =
       cum_weights[(par_idx + 1) * n_particles_each - 1] - start_val;
-    real_type u_particle = normalising_constant /
-                        static_cast<real_type>(n_particles_each) *
-                        (u[par_idx] + i % n_particles_each);
-    index[i] = binary_interval_search(
-      cum_weights + par_idx * n_particles_each,
-      n_particles_each, u_particle, start_val) + par_idx * n_particles_each;
+    if (normalising_constant == 0) {
+      index[i] = i;
+    } else {
+      real_type u_particle = normalising_constant /
+                          static_cast<real_type>(n_particles_each) *
+                          (u[par_idx] + i % n_particles_each);
+      index[i] = binary_interval_search(
+        cum_weights + par_idx * n_particles_each,
+        n_particles_each, u_particle, start_val) + par_idx * n_particles_each;
+    }
 #ifdef __NVCC__
   }
 #else

inst/include/dust/gpu/kernels.hpp

@@ -4,8 +4,8 @@
 
 #define DUST_VERSION_MAJOR 0
 #define DUST_VERSION_MINOR 14
-#define DUST_VERSION_PATCH 7
-#define DUST_VERSION_STRING "0.14.7"
-#define DUST_VERSION_CODE 1407
+#define DUST_VERSION_PATCH 8
+#define DUST_VERSION_STRING "0.14.8"
+#define DUST_VERSION_CODE 1408
 
 #endif

inst/include/dust/random/version.hpp

@@ -107,6 +107,11 @@ T clamp(T x, T min, T max) {
   return std::max(std::min(x, max), min);
 }
 
+template <typename real_type, typename It>
+bool all_zero(It begin, It end) {
+  return std::all_of(begin, end, [](real_type x) { return x == 0; });
+}
+
 }
 }