From 741e3ba5d8daf7e6ff8d2440e77a85b9e6cc83c2 Mon Sep 17 00:00:00 2001 From: Paul Lietar Date: Fri, 17 May 2024 16:52:28 +0100 Subject: [PATCH] Revert "Rewrite the exponential decay process in C++. (#285)" This reverts commit b3376d33cda29cc5e4d7217fefad4b367f9f9167. --- R/RcppExports.R | 4 -- R/processes.R | 3 +- src/RcppExports.cpp | 13 ---- src/processes.cpp | 117 -------------------------------- tests/testthat/test-processes.R | 23 ------- 5 files changed, 1 insertion(+), 159 deletions(-) delete mode 100644 src/processes.cpp delete mode 100644 tests/testthat/test-processes.R diff --git a/R/RcppExports.R b/R/RcppExports.R index ad658f2c..b1ae9854 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -45,10 +45,6 @@ rainfall <- function(t, g0, g, h, floor) { .Call(`_malariasimulation_rainfall`, t, g0, g, h, floor) } -exponential_process_cpp <- function(variable, rate) { - .Call(`_malariasimulation_exponential_process_cpp`, variable, rate) -} - solver_get_states <- function(solver) { .Call(`_malariasimulation_solver_get_states`, solver) } diff --git a/R/processes.R b/R/processes.R index 12f8c3db..a3396d5f 100644 --- a/R/processes.R +++ b/R/processes.R @@ -285,9 +285,8 @@ create_processes <- function( #' @param rate the exponential rate #' @noRd create_exponential_decay_process <- function(variable, rate) { - stopifnot(inherits(variable, "DoubleVariable")) decay_rate <- exp(-1/rate) - exponential_process_cpp(variable$.variable, decay_rate) + function(timestep) variable$queue_update(variable$get_values() * decay_rate) } #' @title Create and initialise lagged_infectivity object diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index e2a6b9d6..9eed5816 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -179,18 +179,6 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } -// exponential_process_cpp -Rcpp::XPtr exponential_process_cpp(Rcpp::XPtr variable, const double rate); -RcppExport SEXP _malariasimulation_exponential_process_cpp(SEXP variableSEXP, SEXP rateSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< Rcpp::XPtr >::type variable(variableSEXP); - Rcpp::traits::input_parameter< const double >::type rate(rateSEXP); - rcpp_result_gen = Rcpp::wrap(exponential_process_cpp(variable, rate)); - return rcpp_result_gen; -END_RCPP -} // solver_get_states std::vector solver_get_states(Rcpp::XPtr solver); RcppExport SEXP _malariasimulation_solver_get_states(SEXP solverSEXP) { @@ -362,7 +350,6 @@ static const R_CallMethodDef CallEntries[] = { {"_malariasimulation_carrying_capacity", (DL_FUNC) &_malariasimulation_carrying_capacity, 8}, {"_malariasimulation_eggs_laid", (DL_FUNC) &_malariasimulation_eggs_laid, 3}, {"_malariasimulation_rainfall", (DL_FUNC) &_malariasimulation_rainfall, 5}, - {"_malariasimulation_exponential_process_cpp", (DL_FUNC) &_malariasimulation_exponential_process_cpp, 2}, {"_malariasimulation_solver_get_states", (DL_FUNC) &_malariasimulation_solver_get_states, 1}, {"_malariasimulation_solver_set_states", (DL_FUNC) &_malariasimulation_solver_set_states, 3}, {"_malariasimulation_solver_step", (DL_FUNC) &_malariasimulation_solver_step, 1}, diff --git a/src/processes.cpp b/src/processes.cpp deleted file mode 100644 index b5db4bad..00000000 --- a/src/processes.cpp +++ /dev/null @@ -1,117 +0,0 @@ -#include -#include - -/** - * An iterator adaptor which yields the same values as the underlying iterator, - * but scaled by a pre-determined factor. - * - * This is used by the exponential_process below to scale an std::vector by a - * constant. - * - * There are two straightforward ways of performing the operation. The first is - * to create an empty vector, use `reserve(N)` to pre-allocate the vector and - * then call `push_back` with each new value. The second way would be to create - * a zero-initialised vector of size N and then use `operator[]` to fill in the - * values. - * - * Unfortunately both approaches have significant overhead. In the former, the - * use of `push_back` requires repeated checks as to whether the vector needs - * growing, despite the prior reserve call. These calls inhibits optimizations - * such as unrolling and auto-vectorization of the loop. The latter approach - * requires an initial memset when zero-initializing the vector, even though the - * vector then gets overwritten entirely. Sadly gcc fails to optimize out either - * of those. Ideally we want a way to create a pre-sized but uninitialised - * std::vector we can write to ourselves, but there is no API in the standard - * library to do this. All existing workarounds end up with an std::vector with - * non-default item type or allocators. - * - * There is however a way out! std::vector has a constructor which accepts a - * pair of iterators and fills the vector with values from the iterators. Using - * `std::distance` on the iterator pair it can even pre-allocate the vector to - * the right size. No zero-initialisation or no capacity checks, just one - * allocation and a straightforward easily optimizable loop. All we need is an - * iterator yielding the right values, hence `scale_iterator`. In C++20 we would - * probably be able to use the new ranges library as our iterators. - * - * How much does this matter? On microbenchmarks, for small and medium sized - * vector (<= 1M doubles), this version is about 30% faster than the - * zero-initialising implementation and 60% faster than the one which uses - * push_back. For larger vector sizes the difference is less pronounced, - * possibly because caches become saturated. At the time of writing, on a - * real-word run of malariasimulation with a population size of 1M the overall - * speedup is about 2-3%. - * - * https://wolchok.org/posts/cxx-trap-1-constant-size-vector/ - * https://codingnest.com/the-little-things-the-missing-performance-in-std-vector/ - * https://lemire.me/blog/2012/06/20/do-not-waste-time-with-stl-vectors/ - */ -template -struct scale_iterator { - using iterator_category = std::forward_iterator_tag; - using difference_type = typename std::iterator_traits::difference_type; - using value_type = typename std::iterator_traits::value_type; - using pointer = typename std::iterator_traits::pointer; - - // We skirt the rules a bit by returning a prvalue from `operator*`, even - // though the C++17 (and prior) standard says forward iterators are supposed - // to return a reference type (ie. a glvalue). Because the scaling is - // applied on the fly, there is no glvalue we could return a reference to. - // - // An input iterator would be allowed to return a prvalue, but the - // std::vector constructor wouldn't be able to figure out the length ahead - // of time if we were an input iterator. - // - // C++20 actually introduces parallel definitions of input and forward - // iterators, which relax this requirement, so under that classification our - // implementation in correct. - // - // In practice though, this does not really matter. We only use this - // iterator in one specific context, and the vector constructor doesn't do - // anything elaborate that we would be upsetting. - using reference = value_type; - - scale_iterator(underlying_iterator it, value_type factor) : it(it), factor(factor) {} - reference operator*() { - return factor * (*it); - } - bool operator==(const scale_iterator& other) { - return it == other.it; - } - bool operator!=(const scale_iterator& other) { - return it != other.it; - } - scale_iterator& operator++() { - it++; - return *this; - } - scale_iterator operator++(int) { - return scale_iterator(it++, factor); - } - - private: - underlying_iterator it; - value_type factor; -}; - -template -scale_iterator make_scale_iterator(T&& it, typename std::iterator_traits::value_type scale) { - return scale_iterator(std::forward(it), scale); -} - -//[[Rcpp::export]] -Rcpp::XPtr exponential_process_cpp( - Rcpp::XPtr variable, - const double rate -){ - return Rcpp::XPtr( - new process_t([=](size_t t){ - const std::vector& values = variable->get_values(); - std::vector new_values( - make_scale_iterator(values.cbegin(), rate), - make_scale_iterator(values.cend(), rate)); - - variable->queue_update(std::move(new_values), std::vector()); - }), - true - ); -} diff --git a/tests/testthat/test-processes.R b/tests/testthat/test-processes.R deleted file mode 100644 index 29f0d6b0..00000000 --- a/tests/testthat/test-processes.R +++ /dev/null @@ -1,23 +0,0 @@ -test_that("exponential_decay_process works as expected", { - # This rate gives a halving at every timestep - rate <- -1 / log(0.5) - - v <- individual::DoubleVariable$new(c(0,0.5,1,2,4,10)) - p <- create_exponential_decay_process(v, rate) - - individual:::execute_any_process(p, 1) - v$.update() - - expect_equal(v$get_values(), c(0, 0.25, 0.5, 1, 2, 5)) - - individual:::execute_any_process(p, 2) - v$.update() - - expect_equal(v$get_values(), c(0, 0.125, 0.25, 0.5, 1, 2.5)) -}) - -test_that("exponential_decay_process fails on IntegerVariable", { - rate <- -1 / log(0.5) - v <- individual::IntegerVariable$new(c(0,1,2,3)) - expect_error(create_exponential_decay_process(v, rate)) -})