Merge pull request #177 from mrc-ide/dev

Dev

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: individual
 Title: Framework for Specifying and Simulating Individual Based Models
-Version: 0.1.10
+Version: 0.1.11
 Authors@R: c(
   person(
     given = "Giovanni",
@@ -23,6 +23,12 @@ Authors@R: c(
     comment = c(ORCID = "0000-0003-3001-4959"),
     email = '[email protected]'
   ),
+  person(
+    given = "Paul",
+    family = "Liétar",
+    role = c('aut'),
+    email = '[email protected]'
+  ),
   person(
     given = "Imperial College of Science, Technology and Medicine",
     family = "",

NEWS.md

@@ -1,3 +1,8 @@
+# individual 0.1.11
+
+  * Optimised rendering memory usage and speed
+  * Added benchmarks for rendering
+
 # individual 0.1.10
 
   * Vector-based updates are consolidated into one place

R/RcppExports.R

@@ -437,6 +437,18 @@ integer_ragged_variable_queue_shrink_bitset <- function(variable, index) {
     invisible(.Call(`_individual_integer_ragged_variable_queue_shrink_bitset`, variable, index))
 }
 
+create_render_vector <- function(data) {
+    .Call(`_individual_create_render_vector`, data)
+}
+
+render_vector_update <- function(v, index, value) {
+    invisible(.Call(`_individual_render_vector_update`, v, index, value))
+}
+
+render_vector_data <- function(v) {
+    .Call(`_individual_render_vector_data`, v)
+}
+
 execute_process <- function(process, timestep) {
     invisible(.Call(`_individual_execute_process`, process, timestep))
 }
@@ -455,5 +467,5 @@ variable_resize <- function(variable) {
 
 # Register entry points for exported C++ functions
 methods::setLoadAction(function(ns) {
-    .Call('_individual_RcppExport_registerCCallable', PACKAGE = 'individual')
+    .Call(`_individual_RcppExport_registerCCallable`)
 })

R/render.R

@@ -15,8 +15,8 @@ Render <- R6Class(
     #' renderers.
     #' @param timesteps number of timesteps in the simulation.
     initialize = function(timesteps) {
-      private$.timesteps = timesteps
-      private$.vectors[['timestep']] <- seq_len(timesteps)
+      private$.timesteps <- timesteps
+      private$.vectors[['timestep']] <- create_render_vector(seq_len(timesteps))
     },
 
     #' @description
@@ -28,7 +28,7 @@ Render <- R6Class(
       if (name == 'timestep') {
         stop("Cannot set default value for variable 'timestep'")
       }
-      private$.vectors[[name]] = rep(value, private$.timesteps)
+      private$.vectors[[name]] <- create_render_vector(rep(value, private$.timesteps))
     },
 
     #' @description
@@ -41,15 +41,15 @@ Render <- R6Class(
         stop("Please don't name your variable 'timestep'")
       }
       if (!(name %in% names(private$.vectors))) {
-        private$.vectors[[name]] = rep(NA, private$.timesteps)
+        private$.vectors[[name]] <- create_render_vector(rep(NA_real_, private$.timesteps))
       }
-      private$.vectors[[name]][[timestep]] = value
+      render_vector_update(private$.vectors[[name]], timestep, value)
     },
 
     #' @description
     #' Return the render as a \code{\link[base]{data.frame}}.
     to_dataframe = function() {
-      data.frame(private$.vectors)
+      data.frame(lapply(private$.vectors, render_vector_data))
     }
   )
 )

inst/include/RenderVector.h

@@ -0,0 +1,37 @@
+/*
+ *  RenderVector.h
+ *
+ *  Created on: 21 Dec 2023
+ *      Author: pl2113
+ */
+
+#ifndef INST_INCLUDE_RENDER_VECTOR_H_
+#define INST_INCLUDE_RENDER_VECTOR_H_
+
+#include <Rcpp.h>
+
+/**
+ * A thin wrapper around a std::vector<double>, used to provide by-reference
+ * semantics and guaranteed in-place mutation in the Render class.
+ *
+ */
+struct RenderVector {
+    RenderVector(std::vector<double> data) : _data(std::move(data)) { }
+
+    void update(size_t index, double value) {
+        // index is R-style 1-indexed, rather than C's 0-indexing.
+        if (index < 1 || index > _data.size()) {
+            Rcpp::stop("index out-of-bounds");
+        }
+        _data[index - 1] = value;
+    }
+
+    const std::vector<double>& data() const {
+        return _data;
+    }
+
+    private:
+        std::vector<double> _data;
+};
+
+#endif /* INST_INCLUDE_RENDER_VECTOR_H_ */

inst/include/individual_types.h

@@ -15,5 +15,6 @@
 #include "RaggedInteger.h"
 #include "RaggedDouble.h"
 #include "Event.h"
+#include "RenderVector.h"
 
 #endif /* INDIVIDUAL_TYPES_H_ */

src/RcppExports.cpp

@@ -366,7 +366,7 @@ RcppExport SEXP _individual_dummy() {
     if (rcpp_isError_gen) {
         SEXP rcpp_msgSEXP_gen = Rf_asChar(rcpp_result_gen);
         UNPROTECT(1);
-        Rf_error(CHAR(rcpp_msgSEXP_gen));
+        Rf_error("%s", CHAR(rcpp_msgSEXP_gen));
     }
     UNPROTECT(1);
     return rcpp_result_gen;
@@ -1303,6 +1303,40 @@ BEGIN_RCPP
     return R_NilValue;
 END_RCPP
 }
+// create_render_vector
+Rcpp::XPtr<RenderVector> create_render_vector(std::vector<double> data);
+RcppExport SEXP _individual_create_render_vector(SEXP dataSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< std::vector<double> >::type data(dataSEXP);
+    rcpp_result_gen = Rcpp::wrap(create_render_vector(data));
+    return rcpp_result_gen;
+END_RCPP
+}
+// render_vector_update
+void render_vector_update(Rcpp::XPtr<RenderVector> v, size_t index, double value);
+RcppExport SEXP _individual_render_vector_update(SEXP vSEXP, SEXP indexSEXP, SEXP valueSEXP) {
+BEGIN_RCPP
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< Rcpp::XPtr<RenderVector> >::type v(vSEXP);
+    Rcpp::traits::input_parameter< size_t >::type index(indexSEXP);
+    Rcpp::traits::input_parameter< double >::type value(valueSEXP);
+    render_vector_update(v, index, value);
+    return R_NilValue;
+END_RCPP
+}
+// render_vector_data
+std::vector<double> render_vector_data(Rcpp::XPtr<RenderVector> v);
+RcppExport SEXP _individual_render_vector_data(SEXP vSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< Rcpp::XPtr<RenderVector> >::type v(vSEXP);
+    rcpp_result_gen = Rcpp::wrap(render_vector_data(v));
+    return rcpp_result_gen;
+END_RCPP
+}
 // execute_process
 void execute_process(Rcpp::XPtr<process_t> process, size_t timestep);
 RcppExport SEXP _individual_execute_process(SEXP processSEXP, SEXP timestepSEXP) {
@@ -1475,6 +1509,9 @@ static const R_CallMethodDef CallEntries[] = {
     {"_individual_integer_ragged_variable_queue_extend", (DL_FUNC) &_individual_integer_ragged_variable_queue_extend, 2},
     {"_individual_integer_ragged_variable_queue_shrink", (DL_FUNC) &_individual_integer_ragged_variable_queue_shrink, 2},
     {"_individual_integer_ragged_variable_queue_shrink_bitset", (DL_FUNC) &_individual_integer_ragged_variable_queue_shrink_bitset, 2},
+    {"_individual_create_render_vector", (DL_FUNC) &_individual_create_render_vector, 1},
+    {"_individual_render_vector_update", (DL_FUNC) &_individual_render_vector_update, 3},
+    {"_individual_render_vector_data", (DL_FUNC) &_individual_render_vector_data, 1},
     {"_individual_execute_process", (DL_FUNC) &_individual_execute_process, 2},
     {"_individual_variable_get_size", (DL_FUNC) &_individual_variable_get_size, 1},
     {"_individual_variable_update", (DL_FUNC) &_individual_variable_update, 1},

src/render_vector.cpp

@@ -0,0 +1,26 @@
+/*
+ * render_vector.cpp
+ *
+ *  Created on: 21 Dec 2023
+ *      Author: pl2113
+ */
+
+
+#include "../inst/include/RenderVector.h"
+#include <Rcpp.h>
+
+
+//[[Rcpp::export]]
+Rcpp::XPtr<RenderVector> create_render_vector(std::vector<double> data) {
+    return Rcpp::XPtr<RenderVector>(new RenderVector(std::move(data)), true);
+}
+
+//[[Rcpp::export]]
+void render_vector_update(Rcpp::XPtr<RenderVector> v, size_t index, double value) {
+    v->update(index, value);
+}
+
+//[[Rcpp::export]]
+std::vector<double> render_vector_data(Rcpp::XPtr<RenderVector> v) {
+    return v->data();
+}

tests/performance/bench-render.R

@@ -0,0 +1,61 @@
+#
+# bench-render.R
+#
+# Created on: 22 Dec 2023
+#   Author: pl2113
+#
+
+library(individual)
+library(bench)
+library(ggplot2)
+library(scales)
+
+source("./tests/performance/utils.R")
+
+render_single <- bench::press(
+  timesteps=floor(10^seq(3,6,0.25)),
+  {
+    render <- Render$new(timesteps)
+    bench::mark(
+      min_iterations = 50,
+      check = FALSE,
+      render={
+        # Use timesteps/2 to write in the middle of the array
+        render$render("data", 0.5, timesteps/2)
+      })
+  })
+
+render_single %>%
+  simplify_bench_output() %>%
+  ggplot() +
+  aes(x = timesteps, y = as.numeric(time), color=expression, fill=expression, group=as.factor(timesteps):expression) +
+  geom_violin(position=position_dodge(width=0.02), alpha=0.3) +
+  labs(y="time", fill="expression", color="expression") +
+  scale_x_continuous(trans='log10', n.breaks=6, labels = label_comma()) +
+  scale_y_continuous(trans='log10', n.breaks=6, labels = function(x) format(bench::as_bench_time(x))) +
+  ggtitle("Render single timestep benchmark")
+
+render_all <- bench::press(
+  timesteps=floor(10^seq(3,5,0.25)),
+  {
+    data <- runif(timesteps)
+    bench::mark(
+      min_iterations = 5,
+      check = FALSE, 
+      filter_gc = FALSE,
+      render_all={
+        render <- Render$new(timesteps)
+        mapply(function(x, i) render$render("data", x, i), data, seq_along(data))
+      })
+  })
+
+render_all %>%
+  simplify_bench_output(filter_gc=FALSE) %>%
+  ggplot() +
+  aes(x = timesteps, y = as.numeric(time), color=expression, fill=expression, group=as.factor(timesteps):expression) +
+  geom_violin(position=position_dodge(width=0.01), alpha=0.3) +
+  labs(y="time", fill="expression", color="expression") +
+  scale_x_continuous(trans='log10', n.breaks=6, labels = label_comma()) +
+  scale_y_continuous(trans='log10', n.breaks=6, labels = function(x) format(bench::as_bench_time(x))) +
+  ggtitle("Render all timesteps benchmark")
+

tests/performance/utils.R

@@ -48,7 +48,7 @@ create_random_index_bitset <- function(size, limit) {
 #' @description Unnest output to generate histograms or density plots, and remove
 #' all runs where any level of garbage collection was executed.
 #' @param out output of [bench::press] function
-simplify_bench_output <- function(out) {
+simplify_bench_output <- function(out, filter_gc=TRUE) {
   x <- lapply(X = seq_len(nrow(out)), FUN = function(i) {
     # get gc level (if run) as factor
     gc <- rep("none", times = nrow(out$gc[[i]]))
@@ -66,7 +66,10 @@ simplify_bench_output <- function(out) {
     return(out_i)
   })
   out_format <- do.call(what = rbind, args = x)
-  out_format <- out_format[out_format$gc == "none", ]
+  if (filter_gc)
+  {
+    out_format <- out_format[out_format$gc == "none", ]
+  }
   out_format$expression <- as.factor(attr(out_format$expression, "description"))
   return(out_format)
 }

tests/testthat/test-render.R

@@ -51,3 +51,9 @@ test_that("Render default works", {
   rendered <- render$to_dataframe()
   expect_mapequal(true_render, rendered)
 })
+
+test_that("Out of range timestep errors", {
+  render <- Render$new(3)
+  expect_error(render$render('S', 10, 0), "index out-of-bounds")
+  expect_error(render$render('S', 10, 4), "index out-of-bounds")
+})