Rcpp implementation of final_size (#36)

* Changes resulting from merge commits from old feature branches

* WIP Rcpp implementation fails on matrix division

* Fix asserts for r0 and prop_initial_infected

* Minor style change to final_size example

* Add Rcpp Imports and Linking

* Add fn f1 from final_size.R

* WIP implementing Rcpp version of final_size; prep transmission matrix

* WIP skeleton fn f2 for final_size

* WIP implement fn f2 from final_size.R

* WIP correct fn f2 returns (-beta2 + (diag/x))

* WIP complete Rcpp code for final_size as final_size_cpp

* WIP update fn f1, f2 comments

* WIP using Rcpp plugin for C++11

* WIP include iostream for messages

* WIP add first basic test for final_size_cpp

* WIP initial documentation final_size_cpp

* WIP procedural files Rcpp exports

* Correct fn f1, fix operation order

* WIP minor updates to internal variables

* WIP prevent roxygen2 updates to NAMESPACE

* WIP basic check for final_size_cpp

* WIP add test for polymod data for final_size_cpp

* Add file for correct NAMESAPCE generation

* WIP correct namespace file

* WIP add RcppEigen import to DESCRIPTION

* WIP remove internal funs from src file

* WIP internal functions to header file

* WIP include header in src file

* Add catch testing framework for Cpp internal funs

* Add basic test for funs f1, f2

* WIP update DESCRIPTION for catch testing

* Update RcppExports with catch testing

* WIP remove catch default tests

* WIP add RcppEigen dependency to header file

* Correct header extension name

* Allow prop_suscept to be single value or vector

* WIP test Rcpp fun not base R

* Add tests for r0 effect and error for prop_suscep

* Update Rcpp exports

* Update documentation for single or multiple prop_suscep

* Add test for prop_initial_infected

* Style files with styler

* Correct handling of single prop_suscep

* Remove non-portable flags

* Style files with clang-format

* WIP add include guards and header includes

* Remove include guard from Cpp file

* Procedural changes for Rcpp exports

* Remove extra message

* Remove general export statement

Co-authored-by: Hugo Gruson <[email protected]>
Co-authored-by: Thibaut Jombart <[email protected]>

DESCRIPTION

@@ -38,10 +38,17 @@ Encoding: UTF-8
 Suggests: 
     covr,
     socialmixr,
-    testthat (>= 3.0.0)
+    testthat (>= 3.0.0),
+    xml2
 Config/testthat/edition: 3
 License: MIT + file LICENSE
 Remotes:
     epiforecasts/socialmixr
-Imports:
-    checkmate
+Imports: 
+    checkmate,
+    Rcpp,
+    RcppEigen
+LinkingTo: 
+    Rcpp,
+    RcppEigen,
+    testthat

NAMESPACE

@@ -1,3 +1,8 @@
 # Generated by roxygen2: do not edit by hand
 
 export(final_size)
+export(final_size_cpp)
+exportPattern("^[[:alpha:]]+")
+import(RcppEigen)
+importFrom(Rcpp,evalCpp)
+useDynLib(finalsize)

R/RcppExports.R

@@ -0,0 +1,28 @@
+# Generated by using Rcpp::compileAttributes() -> do not edit by hand
+# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+#' @title Calculate final epidemic size with RcppEigen backend
+#'
+#' @description This function calculates final epidemic size using SIR model
+#' for a heterogeneously mixing population.
+#'
+#' @param r0  Basic reproduction number.
+#' @param contact_matrix  Social contact matrix. Entry $mm_{ij}$ gives average
+#' number of contacts in group $i$ reported by participants in group j
+#' @param demography_vector  Demography vector. Entry $pp_{i}$ gives proportion
+#' of total population in group $i$ (model will normalise if needed)
+#' @param prop_initial_infected Proportion of all age groups that is initially
+#' infected. May be a single number, or a vector of proportions infected.
+#' If a vector, must be the same length as the demography vector, otherwise the
+#' vector will be recycled. Default value is 0.001 for all groups.
+#' @param prop_suscep  Proportion of each group susceptible. May be a single
+#' numeric value or a numeric vector of the same length as the demography
+#' vector.
+#' @param iterations Number of solver iterations
+#'
+#' @keywords epidemic model
+#' @export
+final_size_cpp <- function(r0, contact_matrix, demography_vector, prop_initial_infected, prop_suscep, iterations = 30L) {
+    .Call('_finalsize_final_size_cpp', PACKAGE = 'finalsize', r0, contact_matrix, demography_vector, prop_initial_infected, prop_suscep, iterations)
+}
+

R/catch-routine-registration.R

@@ -0,0 +1,6 @@
+# This dummy function definition is included with the package to ensure that
+# 'tools::package_native_routine_registration_skeleton()' generates the required
+# registration info for the 'run_testthat_tests' symbol.
+(function() {
+  .Call("run_testthat_tests", FALSE, PACKAGE = "finalsize")
+})

R/finalsize-package.R

@@ -0,0 +1,4 @@
+#' @useDynLib finalsize
+#' @import RcppEigen
+#' @importFrom Rcpp evalCpp
+NULL

src/.gitignore

@@ -0,0 +1,3 @@
+*.o
+*.so
+*.dll

src/Makevars

@@ -0,0 +1,2 @@
+PKG_CXXFLAGS = -I../inst/include
+USE_CXX1X = true

src/Makevars.win

@@ -0,0 +1,2 @@
+PKG_CXXFLAGS = -w -I../inst/include
+USE_CXX1X = true

src/RcppExports.cpp

@@ -0,0 +1,56 @@
+// Generated by using Rcpp::compileAttributes() -> do not edit by hand
+// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+#include <Rcpp.h>
+#include <RcppEigen.h>
+
+using namespace Rcpp;
+
+#ifdef RCPP_USE_GLOBAL_ROSTREAM
+Rcpp::Rostream<true>& Rcpp::Rcout = Rcpp::Rcpp_cout_get();
+Rcpp::Rostream<false>& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get();
+#endif
+
+// final_size_cpp
+Eigen::VectorXd final_size_cpp(const double& r0,
+                               const Eigen::MatrixXd& contact_matrix,
+                               const Eigen::VectorXd& demography_vector,
+                               const Eigen::VectorXd& prop_initial_infected,
+                               Eigen::VectorXd prop_suscep,
+                               const int iterations);
+RcppExport SEXP _finalsize_final_size_cpp(SEXP r0SEXP, SEXP contact_matrixSEXP,
+                                          SEXP demography_vectorSEXP,
+                                          SEXP prop_initial_infectedSEXP,
+                                          SEXP prop_suscepSEXP,
+                                          SEXP iterationsSEXP) {
+  BEGIN_RCPP
+  Rcpp::RObject rcpp_result_gen;
+  Rcpp::RNGScope rcpp_rngScope_gen;
+  Rcpp::traits::input_parameter<const double&>::type r0(r0SEXP);
+  Rcpp::traits::input_parameter<const Eigen::MatrixXd&>::type contact_matrix(
+      contact_matrixSEXP);
+  Rcpp::traits::input_parameter<const Eigen::VectorXd&>::type demography_vector(
+      demography_vectorSEXP);
+  Rcpp::traits::input_parameter<const Eigen::VectorXd&>::type
+      prop_initial_infected(prop_initial_infectedSEXP);
+  Rcpp::traits::input_parameter<Eigen::VectorXd>::type prop_suscep(
+      prop_suscepSEXP);
+  Rcpp::traits::input_parameter<const int>::type iterations(iterationsSEXP);
+  rcpp_result_gen = Rcpp::wrap(
+      final_size_cpp(r0, contact_matrix, demography_vector,
+                     prop_initial_infected, prop_suscep, iterations));
+  return rcpp_result_gen;
+  END_RCPP
+}
+
+RcppExport SEXP run_testthat_tests(SEXP);
+
+static const R_CallMethodDef CallEntries[] = {
+    {"_finalsize_final_size_cpp", (DL_FUNC)&_finalsize_final_size_cpp, 6},
+    {"run_testthat_tests", (DL_FUNC)&run_testthat_tests, 1},
+    {NULL, NULL, 0}};
+
+RcppExport void R_init_finalsize(DllInfo* dll) {
+  R_registerRoutines(dll, NULL, CallEntries, NULL, NULL);
+  R_useDynamicSymbols(dll, FALSE);
+}

src/finalsize.cpp

@@ -0,0 +1,121 @@
+#include "finalsize.h"
+
+#include <Rcpp.h>
+#include <RcppEigen.h>
+
+#include <iostream>
+
+// [[Rcpp::plugins(cpp11)]]
+// via the depends attribute we tell Rcpp to create hooks for
+// RcppEigen so that the build process will know what to do
+//
+// [[Rcpp::depends(RcppEigen)]]
+
+//' @title Calculate final epidemic size with RcppEigen backend
+//'
+//' @description This function calculates final epidemic size using SIR model
+//' for a heterogeneously mixing population.
+//'
+//' @param r0  Basic reproduction number.
+//' @param contact_matrix  Social contact matrix. Entry $mm_{ij}$ gives average
+//' number of contacts in group $i$ reported by participants in group j
+//' @param demography_vector  Demography vector. Entry $pp_{i}$ gives proportion
+//' of total population in group $i$ (model will normalise if needed)
+//' @param prop_initial_infected Proportion of all age groups that is initially
+//' infected. May be a single number, or a vector of proportions infected.
+//' If a vector, must be the same length as the demography vector, otherwise the
+//' vector will be recycled. Default value is 0.001 for all groups.
+//' @param prop_suscep  Proportion of each group susceptible. May be a single
+//' numeric value or a numeric vector of the same length as the demography
+//' vector.
+//' @param iterations Number of solver iterations
+//'
+//' @keywords epidemic model
+//' @export
+// [[Rcpp::export]]
+Eigen::VectorXd final_size_cpp(const double &r0,
+                               const Eigen::MatrixXd &contact_matrix,
+                               const Eigen::VectorXd &demography_vector,
+                               const Eigen::VectorXd &prop_initial_infected,
+                               Eigen::VectorXd prop_suscep,
+                               const int iterations = 30) {
+  if (prop_suscep.size() == 1) {
+    const double prop_suscep_ = prop_suscep[0];
+    prop_suscep.resize(demography_vector.size(), 1);
+    prop_suscep.fill(prop_suscep_);
+  } else if (prop_suscep.size() != demography_vector.size()) {
+    Rcpp::stop("Error: prop_suscep must be same size as demography vector");
+  }
+
+  // scale demography vector
+  Eigen::VectorXd pp0 = demography_vector / (demography_vector.sum());
+
+  // largest real eigenvalue of the contact matrix
+  Eigen::EigenSolver<Eigen::MatrixXd> es(contact_matrix, false);
+  Eigen::MatrixXd eig_vals = es.eigenvalues().real();
+  double eig_val_max = eig_vals.maxCoeff();
+  // scale the next generation matrix for max eigenvalue = r0
+  Eigen::MatrixXd mm0 = r0 * (contact_matrix / eig_val_max);
+
+  // define transmission matrix A = mm_{ij} * pp_{j} / pp_{i}
+  Eigen::MatrixXd beta1 = mm0.array();
+  // rowwise division by age group proportion
+  // this is vectorised in final_size.R
+  // with demography vector recycled over columns; i.e., age 1 * row 1 etc.
+  for (size_t i = 0; i < beta1.rows(); i++)
+    beta1.row(i) = beta1.row(i) * prop_suscep[i] / pp0[i];
+
+  Eigen::MatrixXd beta2 = beta1.transpose();
+  // rowwise multiplication with corresponding age group proportion
+  for (size_t i = 0; i < beta2.rows(); i++)
+    beta2.row(i) = beta2.row(i) * pp0[i];
+
+  beta2 = (beta2.array().transpose()).matrix();  // could be more concise?
+
+  // Newton solver for final size equation A(1-x) = -log(x)
+  // get number of age groups
+  const size_t v_size = pp0.size();
+
+  // prepare timesteps as iterations matrix
+  // fill with -1.0
+  Eigen::MatrixXd iterate_output;
+  iterate_output.resize(static_cast<size_t>(iterations), v_size);
+  iterate_output.fill(-1.0);
+
+  // prepare initial infections vector
+  Eigen::VectorXd x0;
+  if (prop_initial_infected.size() == 1) {
+    x0 = prop_initial_infected[0] * pp0;
+  } else {
+    if (prop_initial_infected.size() != pp0.size()) {
+      Rcpp::stop(
+          "Error: prop_initial_infection must be same size as demography "
+          "vector\n");
+    }
+    Rcpp::Rcout << "multiple prop_initial_infected\n";
+
+    x0 = prop_initial_infected.array() * pp0.array();
+  }
+
+  // iterate over timesteps
+  iterate_output.row(0) = x0;  // initial proportion infected
+
+  // holding matrix (vec) for dx, the change in infection proportions
+  Eigen::MatrixXd f1_m, f2_m, dx;
+  // update the size of the current outbreak
+  // take the 1st (0) column of dx, as all cols per row are same
+  // iterate_output.row(i) = iterate_output.row(i - 1) + dx.col(0).array();
+  // starting at second row
+  for (size_t i = 1; i < static_cast<size_t>(iterations); i++) {
+    f2_m = f2(beta2, iterate_output.row(i - 1), v_size);
+    f1_m = -f1(beta2, iterate_output.row(i - 1));
+    dx = f2_m.partialPivLu().solve(f1_m);
+    // crude iteration-and-column wise addition
+    for (size_t j = 0; j < iterate_output.cols(); j++) {
+      iterate_output(i, j) = iterate_output(i - 1, j) + dx.col(0)[j];
+    }
+  }
+
+  // return 1 - (vector of final proportions of each age group infected)
+  return (1.0 - (iterate_output.row(iterate_output.rows() - 1)).array());
+}

src/finalsize.h

@@ -0,0 +1,42 @@
+#ifndef FINALSIZE_H
+#define FINALSIZE_H
+
+#include <Rcpp.h>
+#include <RcppEigen.h>
+
+// [[Rcpp::depends(RcppEigen)]]
+
+/// function f1 defined in final_size.R
+// f1 <- function(beta2, x) {
+//     beta2 %*% (1 - x) + log(x)
+//   }
+inline Eigen::MatrixXd f1(Eigen::MatrixXd beta2, const Eigen::VectorXd x) {
+  Eigen::MatrixXd x_ = (beta2 * ((Eigen::VectorXd::Ones(x.size()) - x))) +
+                       ((x.array().log()).matrix());
+
+  return x_;
+}
+
+/// function f2 defined in final_size.R
+// f2 <- function(beta2, x, size) {
+//     -beta2 + diag(size) / x
+//   }
+inline Eigen::MatrixXd f2(Eigen::MatrixXd beta2, const Eigen::VectorXd x,
+                          const size_t &size) {
+  // make diagonal matrix of dims [size, size]
+  // size should be the number of age groups
+  Eigen::VectorXd v = Eigen::VectorXd::Ones(size);
+  Eigen::MatrixXd diag_size;
+  diag_size.resize(size, size);
+  diag_size.fill(0.0);
+  diag_size.diagonal() = v;
+
+  // divide diagonal matrix of 1s by x
+  // x is the current proportion infected per age group
+  for (size_t i = 0; i < diag_size.rows(); i++)
+    diag_size.row(i) = diag_size.row(i) / x[i];
+
+  return -beta2 + diag_size;
+}
+
+#endif

src/test-internal_funs.cpp

@@ -0,0 +1,49 @@
+/*
+ * This file uses the Catch unit testing library, alongside
+ * testthat's simple bindings, to test a C++ function.
+ *
+ * For your own packages, ensure that your test files are
+ * placed within the `src/` folder, and that you include
+ * `LinkingTo: testthat` within your DESCRIPTION file.
+ */
+
+// All test files should include the <testthat.h>
+// header file.
+#include <Rcpp.h>
+#include <RcppEigen.h>
+#include <testthat.h>
+
+#include "finalsize.h"
+
+// check that function f1 returns output of expected dimensions
+context("Fn f1 return is expected length") {
+  test_that("Function f1 returns a 1 row matrix") {
+    // set up test inputs
+    int n_groups = 3;
+    Eigen::MatrixXd beta2;
+    beta2.resize(n_groups, n_groups);
+    beta2.fill(0.2);  // fill with an arbitrary number
+
+    Eigen::Vector3d x(0.2, 0.2, 0.2);  // same length as n_groups
+    CATCH_REQUIRE(f1(beta2, x).size() == n_groups);
+    CATCH_REQUIRE(f1(beta2, x).rows() == n_groups);
+    CATCH_REQUIRE(f1(beta2, x).cols() == 1);
+  }
+}
+
+// check that function f2 returns output of expected dimensions
+context("Fn f2 return is expected size") {
+  test_that("Function f1 returns a 1 row matrix") {
+    // set up test inputs
+    int n_groups = 3;
+    Eigen::MatrixXd beta2;
+    beta2.resize(n_groups, n_groups);
+    beta2.fill(0.2);  // fill with an arbitrary number
+
+    Eigen::Vector3d x(0.2, 0.2, 0.2);  // same length as n_groups
+
+    CATCH_REQUIRE(f2(beta2, x, n_groups).size() == beta2.size());
+    CATCH_REQUIRE(f2(beta2, x, n_groups).rows() == n_groups);
+    CATCH_REQUIRE(f2(beta2, x, n_groups).cols() == n_groups);
+  }
+}