From 8a2f7b9953889fa7a267c949d773a92733b50363 Mon Sep 17 00:00:00 2001
From: pratikunterwegs <pratikgupte16@gmail.com>
Date: Thu, 25 Jan 2024 12:41:20 +0000
Subject: [PATCH 01/18] Add `population_change` Cpp struct

---
 inst/include/population.h | 33 +++++++++++++++++++++++++++++----
 1 file changed, 29 insertions(+), 4 deletions(-)

diff --git a/inst/include/population.h b/inst/include/population.h
index e7c45280..56870b90 100644
--- a/inst/include/population.h
+++ b/inst/include/population.h
@@ -3,13 +3,11 @@
 #ifndef INST_INCLUDE_POPULATION_H_
 #define INST_INCLUDE_POPULATION_H_
 
-// [[Rcpp::plugins(cpp14)]]
-// [[Rcpp::depends(BH)]]
-// [[Rcpp::depends(RcppEigen)]]
-
 // clang-format off
 #include <Rcpp.h>
 #include <RcppEigen.h>
+
+#include <vector>
 // clang-format on
 
 // add to namespace population
@@ -41,6 +39,33 @@ inline Rcpp::NumericMatrix get_initial_conditions(
   return initial_conditions;
 }
 
+/// @brief Hold parameters for changes to the model population.
+/// This is special functionality for the diphtheria model applied to
+/// humanitarian camps.
+struct population_change {
+  const std::vector<double> times;
+  const std::vector<Rcpp::NumericVector>
+      value;  // note this is a proportion of the total
+  const int n_demo_groups;
+
+  /// @brief Constructor for the population change struct
+  /// @param times The times at which the population changes. Changes apply to
+  /// the SUSCEPTIBLES compartment only; the assumption is that the camp is the
+  /// locus of the outbreak.
+  /// @param value The proportional change in the population size.
+  /// Each demographic group increases by the same proportion.
+  population_change(const std::vector<double> &times,
+                    const std::vector<Rcpp::NumericVector> &value)
+      : times(times), value(value), n_demo_groups(value[0].size()) {}
+
+  // member function for population change at time t
+  /// @brief Calculate the population change at time t
+  /// @param t The current timestep
+  /// @return An Eigen Array of population size changes.
+  const Eigen::ArrayXd get_population_change(const double &t);
+};
+
+
 }  // namespace population
 
 #endif  // INST_INCLUDE_POPULATION_H_

From a9eeda106901fc160c787fa06cb9bc3963160507 Mon Sep 17 00:00:00 2001
From: pratikunterwegs <pratikgupte16@gmail.com>
Date: Thu, 25 Jan 2024 12:41:51 +0000
Subject: [PATCH 02/18] Struct member fn to calc current pop change

---
 inst/include/population.h | 20 ++++++++++++++++++++
 1 file changed, 20 insertions(+)

diff --git a/inst/include/population.h b/inst/include/population.h
index 56870b90..20f5b3d4 100644
--- a/inst/include/population.h
+++ b/inst/include/population.h
@@ -65,6 +65,26 @@ struct population_change {
   const Eigen::ArrayXd get_population_change(const double &t);
 };
 
+/// @brief Calculate the population change at time t
+/// @param t The current timestep
+/// @return An Eigen Array of population size changes, which defaults to zeros
+/// if no population change is scheduled
+inline const Eigen::ArrayXd population_change::get_population_change(
+    const double &t) {
+  // empty value for timepoints of no population change
+  Eigen::ArrayXd pop_change(n_demo_groups);
+  pop_change.fill(0.0);
+
+  // crudely check for a match between current time and scheduled pop change
+  for (size_t i = 0; i < times.size(); i++) {
+    if (t > times[i] && t < (times[i] + 1.0)) {
+      pop_change = Rcpp::as<Eigen::ArrayXd>(value[i]);
+      return pop_change;
+    }
+  }
+
+  return pop_change;
+}
 
 }  // namespace population
 

From 9d4b0ff330659833794cae78e2454987ecbb2282 Mon Sep 17 00:00:00 2001
From: pratikunterwegs <pratikgupte16@gmail.com>
Date: Thu, 25 Jan 2024 12:44:28 +0000
Subject: [PATCH 03/18] Add pop change mechanic to diphtheria model

---
 inst/include/epidemic_diphtheria.h | 16 +++++++++++++---
 1 file changed, 13 insertions(+), 3 deletions(-)

diff --git a/inst/include/epidemic_diphtheria.h b/inst/include/epidemic_diphtheria.h
index 3c6e3506..c47a3669 100644
--- a/inst/include/epidemic_diphtheria.h
+++ b/inst/include/epidemic_diphtheria.h
@@ -30,6 +30,9 @@ struct epidemic_diphtheria {
   const std::unordered_map<std::string, intervention::rate_intervention>
       interventions;
   const Rcpp::List time_dependence;
+  /// NOTE: not marked const as this interferes with member function
+  /// TODO: diagnose issue and achieve functioning marking as const
+  population::population_change pop_change;
 
   /// @brief Constructor for the diphtheria epidemic struct
   /// @param model_params An unordered map of string-double pairs, with the
@@ -53,11 +56,13 @@ struct epidemic_diphtheria {
       const std::unordered_map<std::string, double>& model_params,
       const std::unordered_map<std::string, intervention::rate_intervention>&
           interventions,
-      const Rcpp::List& time_dependence)
+      const Rcpp::List& time_dependence,
+      const population::population_change& pop_change)
       : model_params(model_params),
         model_params_temp(model_params),
         interventions(interventions),
-        time_dependence(time_dependence) {}
+        time_dependence(time_dependence),
+        pop_change(pop_change) {}
 
   /// @brief Operator for the diphtheria model
   /// @param x The initial state of the population - rows represent age groups
@@ -85,6 +90,11 @@ struct epidemic_diphtheria {
     // 0|1|2|3|4
     // S|E|I|H|R
 
+    // calculate any population changes
+    const Eigen::ArrayXd current_pop_change =
+        pop_change.get_population_change(t);
+
+    // calculate total infections
     const double total_infections = (x.col(2).array()).sum();
 
     // compartmental transitions without accounting for stratified contacts
@@ -107,7 +117,7 @@ struct epidemic_diphtheria {
     // β: transmissibility; σ: infectiousness rate; γ: recovery rate
     // ν: reporting rate; τ1: 1 / time to hospitalisation;
     // τ2: 1 / time to discharge from hospital; η: prop. hospitalised
-    dxdt.col(0) = -sToE;        // -β*S*I
+    dxdt.col(0) = -sToE + current_pop_change;  // -β*S*I + pop movements
     dxdt.col(1) = sToE - eToI;  // β*S*I - σ*E
     dxdt.col(2) = eToI - iToR;  // σ*E - γ*I
     dxdt.col(3) = iToH - hToR;  // τ1*η*ν*I - τ2*H

From 5cf793fd5e814bf503b8deb73133cebb05027aaf Mon Sep 17 00:00:00 2001
From: pratikunterwegs <pratikgupte16@gmail.com>
Date: Thu, 25 Jan 2024 14:05:57 +0000
Subject: [PATCH 04/18] Population change Cpp struct uses Rcpp objs

---
 inst/include/population.h | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/inst/include/population.h b/inst/include/population.h
index 20f5b3d4..9971f96c 100644
--- a/inst/include/population.h
+++ b/inst/include/population.h
@@ -43,20 +43,20 @@ inline Rcpp::NumericMatrix get_initial_conditions(
 /// This is special functionality for the diphtheria model applied to
 /// humanitarian camps.
 struct population_change {
-  const std::vector<double> times;
-  const std::vector<Rcpp::NumericVector>
-      value;  // note this is a proportion of the total
+  const Rcpp::NumericVector times;
+  const Rcpp::List value;  // note these are actually absolute values
   const int n_demo_groups;
 
   /// @brief Constructor for the population change struct
   /// @param times The times at which the population changes. Changes apply to
   /// the SUSCEPTIBLES compartment only; the assumption is that the camp is the
   /// locus of the outbreak.
-  /// @param value The proportional change in the population size.
-  /// Each demographic group increases by the same proportion.
-  population_change(const std::vector<double> &times,
-                    const std::vector<Rcpp::NumericVector> &value)
-      : times(times), value(value), n_demo_groups(value[0].size()) {}
+  /// @param value The ABSOLUTE change in the population size; an Rcpp List
+  /// of Rcpp NumericVectors, each of the same length as the number of
+  /// demographic groups.
+  population_change(const Rcpp::NumericVector &times, const Rcpp::List &value,
+                    const int &n_demo_groups)
+      : times(times), value(value), n_demo_groups(n_demo_groups) {}
 
   // member function for population change at time t
   /// @brief Calculate the population change at time t

From 5f43382015173ca09932d0cd7202534df5e34225 Mon Sep 17 00:00:00 2001
From: pratikunterwegs <pratikgupte16@gmail.com>
Date: Thu, 25 Jan 2024 14:10:54 +0000
Subject: [PATCH 05/18] Diphtheria model src file + RcppExports updated

---
 inst/include/epidemic_diphtheria.h |  8 ++++----
 src/RcppExports.cpp                | 10 ++++++----
 src/epidemic_diphtheria.cpp        | 12 +++++++++++-
 3 files changed, 21 insertions(+), 9 deletions(-)

diff --git a/inst/include/epidemic_diphtheria.h b/inst/include/epidemic_diphtheria.h
index c47a3669..2768dc94 100644
--- a/inst/include/epidemic_diphtheria.h
+++ b/inst/include/epidemic_diphtheria.h
@@ -118,10 +118,10 @@ struct epidemic_diphtheria {
     // ν: reporting rate; τ1: 1 / time to hospitalisation;
     // τ2: 1 / time to discharge from hospital; η: prop. hospitalised
     dxdt.col(0) = -sToE + current_pop_change;  // -β*S*I + pop movements
-    dxdt.col(1) = sToE - eToI;  // β*S*I - σ*E
-    dxdt.col(2) = eToI - iToR;  // σ*E - γ*I
-    dxdt.col(3) = iToH - hToR;  // τ1*η*ν*I - τ2*H
-    dxdt.col(4) = iToR + hToR;  // γ*I + τ2*H
+    dxdt.col(1) = sToE - eToI;                 // β*S*I - σ*E
+    dxdt.col(2) = eToI - iToR;                 // σ*E - γ*I
+    dxdt.col(3) = iToH - hToR;                 // τ1*η*ν*I - τ2*H
+    dxdt.col(4) = iToR + hToR;                 // γ*I + τ2*H
   }
 };
 
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index f8795272..774b2372 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -38,8 +38,8 @@ BEGIN_RCPP
 END_RCPP
 }
 // model_diphtheria_cpp_internal
-Rcpp::List model_diphtheria_cpp_internal(const Eigen::MatrixXd& initial_state, const double& transmissibility, const double& infectiousness_rate, const double& recovery_rate, const double& reporting_rate, const double& prop_hosp, const double& hosp_entry_rate, const double& hosp_exit_rate, const Rcpp::List& rate_interventions, const Rcpp::List& time_dependence, const double& time_end, const double& increment);
-RcppExport SEXP _epidemics_model_diphtheria_cpp_internal(SEXP initial_stateSEXP, SEXP transmissibilitySEXP, SEXP infectiousness_rateSEXP, SEXP recovery_rateSEXP, SEXP reporting_rateSEXP, SEXP prop_hospSEXP, SEXP hosp_entry_rateSEXP, SEXP hosp_exit_rateSEXP, SEXP rate_interventionsSEXP, SEXP time_dependenceSEXP, SEXP time_endSEXP, SEXP incrementSEXP) {
+Rcpp::List model_diphtheria_cpp_internal(const Eigen::MatrixXd& initial_state, const double& transmissibility, const double& infectiousness_rate, const double& recovery_rate, const double& reporting_rate, const double& prop_hosp, const double& hosp_entry_rate, const double& hosp_exit_rate, const Rcpp::List& rate_interventions, const Rcpp::List& time_dependence, const Rcpp::NumericVector& pop_change_times, const Rcpp::List& pop_change_values, const double& time_end, const double& increment);
+RcppExport SEXP _epidemics_model_diphtheria_cpp_internal(SEXP initial_stateSEXP, SEXP transmissibilitySEXP, SEXP infectiousness_rateSEXP, SEXP recovery_rateSEXP, SEXP reporting_rateSEXP, SEXP prop_hospSEXP, SEXP hosp_entry_rateSEXP, SEXP hosp_exit_rateSEXP, SEXP rate_interventionsSEXP, SEXP time_dependenceSEXP, SEXP pop_change_timesSEXP, SEXP pop_change_valuesSEXP, SEXP time_endSEXP, SEXP incrementSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
@@ -53,9 +53,11 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< const double& >::type hosp_exit_rate(hosp_exit_rateSEXP);
     Rcpp::traits::input_parameter< const Rcpp::List& >::type rate_interventions(rate_interventionsSEXP);
     Rcpp::traits::input_parameter< const Rcpp::List& >::type time_dependence(time_dependenceSEXP);
+    Rcpp::traits::input_parameter< const Rcpp::NumericVector& >::type pop_change_times(pop_change_timesSEXP);
+    Rcpp::traits::input_parameter< const Rcpp::List& >::type pop_change_values(pop_change_valuesSEXP);
     Rcpp::traits::input_parameter< const double& >::type time_end(time_endSEXP);
     Rcpp::traits::input_parameter< const double& >::type increment(incrementSEXP);
-    rcpp_result_gen = Rcpp::wrap(model_diphtheria_cpp_internal(initial_state, transmissibility, infectiousness_rate, recovery_rate, reporting_rate, prop_hosp, hosp_entry_rate, hosp_exit_rate, rate_interventions, time_dependence, time_end, increment));
+    rcpp_result_gen = Rcpp::wrap(model_diphtheria_cpp_internal(initial_state, transmissibility, infectiousness_rate, recovery_rate, reporting_rate, prop_hosp, hosp_entry_rate, hosp_exit_rate, rate_interventions, time_dependence, pop_change_times, pop_change_values, time_end, increment));
     return rcpp_result_gen;
 END_RCPP
 }
@@ -92,7 +94,7 @@ END_RCPP
 
 static const R_CallMethodDef CallEntries[] = {
     {"_epidemics_model_default_cpp_internal", (DL_FUNC) &_epidemics_model_default_cpp_internal, 15},
-    {"_epidemics_model_diphtheria_cpp_internal", (DL_FUNC) &_epidemics_model_diphtheria_cpp_internal, 12},
+    {"_epidemics_model_diphtheria_cpp_internal", (DL_FUNC) &_epidemics_model_diphtheria_cpp_internal, 14},
     {"_epidemics_model_vacamole_cpp_internal", (DL_FUNC) &_epidemics_model_vacamole_cpp_internal, 20},
     {NULL, NULL, 0}
 };
diff --git a/src/epidemic_diphtheria.cpp b/src/epidemic_diphtheria.cpp
index 788318e1..fd993d7d 100644
--- a/src/epidemic_diphtheria.cpp
+++ b/src/epidemic_diphtheria.cpp
@@ -32,6 +32,10 @@
 //' @param rate_interventions A named list of `<rate_intervention>` objects.
 //' @param time_dependence A named list of functions for parameter time
 //' dependence.
+//' @param pop_change_times A numeric vector of times and which the population
+//' of susceptibles changes.
+//' @param pop_change_values An Rcpp List of numeric vectors giving the value of
+//' changes to each demographic group at each change in population.
 //' @param time_end The end time of the simulation.
 //' @param increment The time increment of the simulation.
 //' @return A two element list, where the first element is a list of matrices
@@ -46,6 +50,8 @@ Rcpp::List model_diphtheria_cpp_internal(
     const double &reporting_rate, const double &prop_hosp,
     const double &hosp_entry_rate, const double &hosp_exit_rate,
     const Rcpp::List &rate_interventions, const Rcpp::List &time_dependence,
+    const Rcpp::NumericVector &pop_change_times,
+    const Rcpp::List &pop_change_values,
     const double &time_end = 100.0,  // double required by boost solver
     const double &increment = 1.0) {
   // initial conditions from input
@@ -66,9 +72,13 @@ Rcpp::List model_diphtheria_cpp_internal(
       rate_interventions_cpp =
           intervention::rate_intervention_cpp(rate_interventions);
 
+  // prepare population changes if any
+  const population::population_change pop_change(
+      pop_change_times, pop_change_values, initial_state.rows());
+
   // create a diphtheria epidemic with parameters
   epidemics::epidemic_diphtheria this_model(
-      model_params, rate_interventions_cpp, time_dependence);
+      model_params, rate_interventions_cpp, time_dependence, pop_change);
 
   // prepare storage containers for the observer
   std::vector<odetools::state_type> x_vec;  // is a vector of MatrixXd

From 46421ffde9294cb4fc31848c64d603675661f154 Mon Sep 17 00:00:00 2001
From: pratikunterwegs <pratikgupte16@gmail.com>
Date: Thu, 25 Jan 2024 14:11:38 +0000
Subject: [PATCH 06/18] Diphtheria args checker handles pop_changes

---
 R/check_args_diphtheria.R | 46 +++++++++++++++++++++++++++++++++++++--
 1 file changed, 44 insertions(+), 2 deletions(-)

diff --git a/R/check_args_diphtheria.R b/R/check_args_diphtheria.R
index 2b535d78..fcf62875 100644
--- a/R/check_args_diphtheria.R
+++ b/R/check_args_diphtheria.R
@@ -26,6 +26,8 @@
 #' 'exposed' and 'exposed_vaccinated' to the 'infectious' and
 #' 'infectious_vaccinated' compartments;
 #'
+#' - `recovery_rate`: a single number for the recovery rate from the infection;
+#'
 #'  - `reporting_rate`: a single number for the proportion of infectious cases
 #' reported;
 #'
@@ -35,7 +37,14 @@
 #' - `hosp_entry_rate`, `hosp_exit_rate`: two numbers representing the rate of
 #' entry and exit from the 'hospitalised' compartment;
 #'
-#' - `recovery_rate`: a single number for the recovery rate from the infection;
+#' - `rate_interventions`: an Rcpp List giving the interventions on model
+#' parameters;
+#'
+#' - `time_dependence`: an Rcpp List giving the time-dependent effects on model
+#' parameters in the form of R functions;
+#'
+#' - `pop_change_times` and `pop_change_values`: the times and values of changes
+#' in the population of susceptibles;
 #'
 #'  - `time_end`, `increment`: two numbers for the time at which to end the
 #' simulation, and the value by which the simulation time
@@ -101,11 +110,43 @@
     )
   }
 
+  # handle population change mechanic, check for correct demography groups
+  # check for only times being NULL and expect values also NULL
+  if (is.null(mod_args[["pop_change_times"]])) {
+    mod_args[["pop_change_times"]] <- 0
+    mod_args[["pop_change_values"]] <- list(
+      rep(
+        0, length(get_parameter(
+          mod_args[["population"]], "demography_vector"
+        ))
+      )
+    )
+  } else {
+    checkmate::assert_numeric(
+      mod_args[["pop_change_times"]],
+      lower = 0, finite = TRUE,
+      min.len = 1
+    )
+    checkmate::assert_list(
+      mod_args[["pop_change_values"]],
+      any.missing = FALSE,
+      len = length(mod_args[["pop_change_times"]])
+    )
+    invisible(
+      lapply(
+        mod_args[["pop_change_values"]], checkmate::assert_numeric,
+        len = length(get_parameter(
+          mod_args[["population"]], "demography_vector"
+        ))
+      )
+    )
+  }
+
   # return arguments invisibly
   invisible(mod_args)
 }
 
-#' @title Prepare arguments for the Vacamole epidemic function
+#' @title Prepare arguments for the diphtheria model
 #' @name check_prepare_diphtheria_args
 #' @rdname check_prepare_diphtheria_args
 #' @keywords internal
@@ -117,6 +158,7 @@
 
   # modify initial state by proportion vaccinated
   # NOTE: this assumes that the first column is 'susceptible'
+  # NOTE: there is no explicit vaccinated compartment
   initial_state[, 1] <- initial_state[, 1] * (1 - mod_args[["prop_vaccinated"]])
 
   # return mod args without population and prop_vaccinated

From e694bf8170eb579ba860ecac7c334c805a9d8835 Mon Sep 17 00:00:00 2001
From: pratikunterwegs <pratikgupte16@gmail.com>
Date: Thu, 25 Jan 2024 14:17:27 +0000
Subject: [PATCH 07/18] Update diphtheria model R frontend and docs

---
 R/RcppExports.R                      |  8 +++++--
 R/model_diphtheria.R                 | 32 +++++++++++++++++++++++++---
 man/check_prepare_diphtheria_args.Rd |  8 ++++++-
 man/dot-model_diphtheria_cpp.Rd      |  8 +++++++
 man/model_diphtheria.Rd              | 30 +++++++++++++++++++++-----
 5 files changed, 75 insertions(+), 11 deletions(-)

diff --git a/R/RcppExports.R b/R/RcppExports.R
index 1d1ad629..825e8cb9 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -62,6 +62,10 @@
 #' @param rate_interventions A named list of `<rate_intervention>` objects.
 #' @param time_dependence A named list of functions for parameter time
 #' dependence.
+#' @param pop_change_times A numeric vector of times and which the population
+#' of susceptibles changes.
+#' @param pop_change_values An Rcpp List of numeric vectors giving the value of
+#' changes to each demographic group at each change in population.
 #' @param time_end The end time of the simulation.
 #' @param increment The time increment of the simulation.
 #' @return A two element list, where the first element is a list of matrices
@@ -69,8 +73,8 @@
 #' as specified in the initial conditions matrix.
 #' The second list element is a vector of timesteps.
 #' @keywords internal
-.model_diphtheria_cpp <- function(initial_state, transmissibility, infectiousness_rate, recovery_rate, reporting_rate, prop_hosp, hosp_entry_rate, hosp_exit_rate, rate_interventions, time_dependence, time_end = 100.0, increment = 1.0) {
-    .Call(`_epidemics_model_diphtheria_cpp_internal`, initial_state, transmissibility, infectiousness_rate, recovery_rate, reporting_rate, prop_hosp, hosp_entry_rate, hosp_exit_rate, rate_interventions, time_dependence, time_end, increment)
+.model_diphtheria_cpp <- function(initial_state, transmissibility, infectiousness_rate, recovery_rate, reporting_rate, prop_hosp, hosp_entry_rate, hosp_exit_rate, rate_interventions, time_dependence, pop_change_times, pop_change_values, time_end = 100.0, increment = 1.0) {
+    .Call(`_epidemics_model_diphtheria_cpp_internal`, initial_state, transmissibility, infectiousness_rate, recovery_rate, reporting_rate, prop_hosp, hosp_entry_rate, hosp_exit_rate, rate_interventions, time_dependence, pop_change_times, pop_change_values, time_end, increment)
 }
 
 #' @title Run the RIVM Vacamole model
diff --git a/R/model_diphtheria.R b/R/model_diphtheria.R
index 08b1d106..80aa9596 100644
--- a/R/model_diphtheria.R
+++ b/R/model_diphtheria.R
@@ -9,8 +9,10 @@
 #' The model is based on Finger et al. (2019) and is intended to be used in the
 #' context of internally displaced people (IDP) or refugee camps.
 #' This model accommodates age or demographic structure and allows for a
-#' proportion of each demographic group to be vaccinated and to not contribute
-#' to the outbreak.
+#' proportion of each demographic group to be vaccinated at the start of the
+#' outbreak, and thus to not contribute to the outbreak.
+#' The model also allows for changes to the initial population size, to model
+#' influxes or evacuations from camps.
 #'
 #' @param population An object of the `population` class, which holds a
 #' population contact matrix, a demography vector, and the initial conditions
@@ -52,10 +54,16 @@
 #' @param time_end The maximum number of timesteps over which to run the model.
 #' Taken as days, with a default value of 100 days.
 #' @param increment The size of the time increment. Taken as days, with a
+#' @param population_change A two-element list, with elements named `"time"` and
+#' `"values"`, giving the times of population changes, and the corresponding
+#' changes in the population of each demographic group at those times.
+#' `"time"` must be a numeric vector, while `"values"` must be a list of the
+#' length of `"time"`, with each element a numeric vector of the same length as
+#' the number of demographic groups in `population`.
 #' default value of 1 day.
 #' @details
 #'
-#' ## R and Rcpp implementations
+#' ## Rcpp implementations
 #'
 #' `model_diphtheria_cpp()` is a wrapper function for [.model_diphtheria_cpp()],
 #' an internal C++ function that uses Boost _odeint_ solvers for an SEIHR model.
@@ -86,6 +94,15 @@
 #'
 #' - Recovery rate (\eqn{\gamma}, `recovery_rate`): 0.333, assuming an
 #' infectious period following symptoms, of 3 days.
+#' 
+#' ## Modelling population changes
+#' 
+#' This model allows changes to the number of susceptibles in each demographic
+#' group, to represent influxes or evacuations from the camp as would be
+#' expected in humanitarian relief situations.
+#' Users can specify the times and changes (to each demographic group) of
+#' changes using the `population_changes` argument, to examine the effect on
+#' outbreak dynamics.
 #'
 #' @references
 #' Finger, F., Funk, S., White, K., Siddiqui, M. R., Edmunds, W. J., &
@@ -108,6 +125,7 @@ model_diphtheria_cpp <- function(population,
                                  ),
                                  intervention = NULL,
                                  time_dependence = NULL,
+                                 population_change = NULL,
                                  time_end = 100,
                                  increment = 1) {
   # check class on required inputs
@@ -152,6 +170,12 @@ model_diphtheria_cpp <- function(population,
     )
   )
 
+  # check population change and create
+  checkmate::assert_list(
+    population_change,
+    null.ok = TRUE, len = 2L, types = c("numeric", "list")
+  )
+
   # check the time end and increment
   # restrict increment to lower limit of 1e-6
   checkmate::assert_number(time_end, lower = 0, finite = TRUE)
@@ -169,6 +193,8 @@ model_diphtheria_cpp <- function(population,
     prop_vaccinated = prop_vaccinated,
     intervention = intervention,
     time_dependence = time_dependence,
+    pop_change_times = population_change[["time"]],
+    pop_change_values = population_change[["values"]],
     time_end = time_end, increment = increment
   )
 
diff --git a/man/check_prepare_diphtheria_args.Rd b/man/check_prepare_diphtheria_args.Rd
index 59d17b70..192a2915 100644
--- a/man/check_prepare_diphtheria_args.Rd
+++ b/man/check_prepare_diphtheria_args.Rd
@@ -27,13 +27,19 @@ vaccinated individuals, respectively;
 \item \code{infectiousness_rate}: a single number for the transition rate from the
 'exposed' and 'exposed_vaccinated' to the 'infectious' and
 'infectious_vaccinated' compartments;
+\item \code{recovery_rate}: a single number for the recovery rate from the infection;
 \item \code{reporting_rate}: a single number for the proportion of infectious cases
 reported;
 \item \code{prop_hosp}: a single number for the proportion of reported cases that
 need hospitalisation;
 \item \code{hosp_entry_rate}, \code{hosp_exit_rate}: two numbers representing the rate of
 entry and exit from the 'hospitalised' compartment;
-\item \code{recovery_rate}: a single number for the recovery rate from the infection;
+\item \code{rate_interventions}: an Rcpp List giving the interventions on model
+parameters;
+\item \code{time_dependence}: an Rcpp List giving the time-dependent effects on model
+parameters in the form of R functions;
+\item \code{pop_change_times} and \code{pop_change_values}: the times and values of changes
+in the population of susceptibles;
 \item \code{time_end}, \code{increment}: two numbers for the time at which to end the
 simulation, and the value by which the simulation time
 is incremented.
diff --git a/man/dot-model_diphtheria_cpp.Rd b/man/dot-model_diphtheria_cpp.Rd
index 2c9b603c..a412d493 100644
--- a/man/dot-model_diphtheria_cpp.Rd
+++ b/man/dot-model_diphtheria_cpp.Rd
@@ -15,6 +15,8 @@
   hosp_exit_rate,
   rate_interventions,
   time_dependence,
+  pop_change_times,
+  pop_change_values,
   time_end = 100,
   increment = 1
 )
@@ -44,6 +46,12 @@ hospital, represented as 1 / time to discharge \eqn{\tau_2}.}
 \item{time_dependence}{A named list of functions for parameter time
 dependence.}
 
+\item{pop_change_times}{A numeric vector of times and which the population
+of susceptibles changes.}
+
+\item{pop_change_values}{An Rcpp List of numeric vectors giving the value of
+changes to each demographic group at each change in population.}
+
 \item{time_end}{The end time of the simulation.}
 
 \item{increment}{The time increment of the simulation.}
diff --git a/man/model_diphtheria.Rd b/man/model_diphtheria.Rd
index 58513dad..902bd0c4 100644
--- a/man/model_diphtheria.Rd
+++ b/man/model_diphtheria.Rd
@@ -17,6 +17,7 @@ model_diphtheria_cpp(
   prop_vaccinated = 0 * get_parameter(population, "demography_vector"),
   intervention = NULL,
   time_dependence = NULL,
+  population_change = NULL,
   time_end = 100,
   increment = 1
 )
@@ -70,11 +71,18 @@ that is dependent on \code{time} (\code{x} represents a model parameter).
 See \strong{Details} for more information, as well as the vignette on time-
 dependence \code{vignette("time_dependence", package = "epidemics")}.}
 
+\item{population_change}{A two-element list, with elements named \code{"time"} and
+\code{"values"}, giving the times of population changes, and the corresponding
+changes in the population of each demographic group at those times.
+\code{"time"} must be a numeric vector, while \code{"values"} must be a list of the
+length of \code{"time"}, with each element a numeric vector of the same length as
+the number of demographic groups in \code{population}.
+default value of 1 day.}
+
 \item{time_end}{The maximum number of timesteps over which to run the model.
 Taken as days, with a default value of 100 days.}
 
-\item{increment}{The size of the time increment. Taken as days, with a
-default value of 1 day.}
+\item{increment}{The size of the time increment. Taken as days, with a}
 }
 \value{
 A \verb{<data.frame>} with the columns "time", "compartment", "age_group",
@@ -87,11 +95,13 @@ compartmental ordinary differential equation model with the compartments
 The model is based on Finger et al. (2019) and is intended to be used in the
 context of internally displaced people (IDP) or refugee camps.
 This model accommodates age or demographic structure and allows for a
-proportion of each demographic group to be vaccinated and to not contribute
-to the outbreak.
+proportion of each demographic group to be vaccinated at the start of the
+outbreak, and thus to not contribute to the outbreak.
+The model also allows for changes to the initial population size, to model
+influxes or evacuations from camps.
 }
 \details{
-\subsection{R and Rcpp implementations}{
+\subsection{Rcpp implementations}{
 
 \code{model_diphtheria_cpp()} is a wrapper function for \code{\link[=.model_diphtheria_cpp]{.model_diphtheria_cpp()}},
 an internal C++ function that uses Boost \emph{odeint} solvers for an SEIHR model.
@@ -119,6 +129,16 @@ individuals are discharged from hospital after 5 days.
 infectious period following symptoms, of 3 days.
 }
 }
+
+\subsection{Modelling population changes}{
+
+This model allows changes to the number of susceptibles in each demographic
+group, to represent influxes or evacuations from the camp as would be
+expected in humanitarian relief situations.
+Users can specify the times and changes (to each demographic group) of
+changes using the \code{population_changes} argument, to examine the effect on
+outbreak dynamics.
+}
 }
 \references{
 Finger, F., Funk, S., White, K., Siddiqui, M. R., Edmunds, W. J., &

From 62199ff60eae89ffa3c5da18df107fea955ca1fd Mon Sep 17 00:00:00 2001
From: pratikunterwegs <pratikgupte16@gmail.com>
Date: Thu, 25 Jan 2024 14:35:10 +0000
Subject: [PATCH 08/18] Fix formatting

---
 R/model_diphtheria.R | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/R/model_diphtheria.R b/R/model_diphtheria.R
index 80aa9596..26da4c1a 100644
--- a/R/model_diphtheria.R
+++ b/R/model_diphtheria.R
@@ -94,9 +94,9 @@
 #'
 #' - Recovery rate (\eqn{\gamma}, `recovery_rate`): 0.333, assuming an
 #' infectious period following symptoms, of 3 days.
-#' 
+#'
 #' ## Modelling population changes
-#' 
+#'
 #' This model allows changes to the number of susceptibles in each demographic
 #' group, to represent influxes or evacuations from the camp as would be
 #' expected in humanitarian relief situations.

From d334140768c2b703825b8a0d05e8529fc907fb4f Mon Sep 17 00:00:00 2001
From: pratikunterwegs <pratikgupte16@gmail.com>
Date: Thu, 25 Jan 2024 16:59:01 +0000
Subject: [PATCH 09/18] Correct calculation of infectious/hospitalised

---
 inst/include/epidemic_diphtheria.h | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/inst/include/epidemic_diphtheria.h b/inst/include/epidemic_diphtheria.h
index 2768dc94..99c96a26 100644
--- a/inst/include/epidemic_diphtheria.h
+++ b/inst/include/epidemic_diphtheria.h
@@ -109,7 +109,10 @@ struct epidemic_diphtheria {
     Eigen::ArrayXd iToH =
         model_params_temp["prop_hosp"] * model_params_temp["reporting_rate"] *
         model_params_temp["hosp_entry_rate"] * x.col(2).array();
-    Eigen::ArrayXd iToR = model_params_temp["recovery_rate"] * x.col(2).array();
+
+    // recoveries are from the infectious/ed who are NOT hospitalised
+    Eigen::ArrayXd iToR =
+        model_params_temp["recovery_rate"] * (x.col(2).array() - iToH);
     Eigen::ArrayXd hToR =
         model_params_temp["hosp_exit_rate"] * x.col(3).array();
 
@@ -119,7 +122,7 @@ struct epidemic_diphtheria {
     // τ2: 1 / time to discharge from hospital; η: prop. hospitalised
     dxdt.col(0) = -sToE + current_pop_change;  // -β*S*I + pop movements
     dxdt.col(1) = sToE - eToI;                 // β*S*I - σ*E
-    dxdt.col(2) = eToI - iToR;                 // σ*E - γ*I
+    dxdt.col(2) = eToI - iToR - iToH;          // σ*E - γ*I - τ1*η*ν*I
     dxdt.col(3) = iToH - hToR;                 // τ1*η*ν*I - τ2*H
     dxdt.col(4) = iToR + hToR;                 // γ*I + τ2*H
   }

From f79973f9f1601e4380f720f9b5db9a68dbb44289 Mon Sep 17 00:00:00 2001
From: pratikunterwegs <pratikgupte16@gmail.com>
Date: Fri, 26 Jan 2024 09:01:04 +0000
Subject: [PATCH 10/18] Correct pop change time comparison

---
 inst/include/population.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/inst/include/population.h b/inst/include/population.h
index 9971f96c..5a03ccb4 100644
--- a/inst/include/population.h
+++ b/inst/include/population.h
@@ -77,7 +77,7 @@ inline const Eigen::ArrayXd population_change::get_population_change(
 
   // crudely check for a match between current time and scheduled pop change
   for (size_t i = 0; i < times.size(); i++) {
-    if (t > times[i] && t < (times[i] + 1.0)) {
+    if (t >= times[i] && t < (times[i] + 1.0)) {
       pop_change = Rcpp::as<Eigen::ArrayXd>(value[i]);
       return pop_change;
     }

From d6effe482e9d6fcc54126b5ac5fbeaecec7b7b49 Mon Sep 17 00:00:00 2001
From: pratikunterwegs <pratikgupte16@gmail.com>
Date: Fri, 26 Jan 2024 09:55:08 +0000
Subject: [PATCH 11/18] Test population change for diphtheria model

---
 tests/testthat/test-model_diphtheria.R | 77 ++++++++++++++++++++++++++
 1 file changed, 77 insertions(+)

diff --git a/tests/testthat/test-model_diphtheria.R b/tests/testthat/test-model_diphtheria.R
index 6f781c82..b4e7b314 100644
--- a/tests/testthat/test-model_diphtheria.R
+++ b/tests/testthat/test-model_diphtheria.R
@@ -67,3 +67,80 @@ test_that("Diptheria model, basic expectations", {
   # NOTE: no checks for final state equal to demography vector as
   # vaccinated individuals are removed from model
 })
+
+# Expectations for the diphtheria model with changed population sizes
+test_that("Diphtheria model with population size changes", {
+  # model runs with population change functionality
+  # create population change data
+  p <- list(
+    time = 70,
+    values = list(
+      c(1e4, 2e5, 1e5)
+    )
+  )
+
+  expect_no_condition(
+    model_diphtheria_cpp(
+      population = camp_pop,
+      population_change = p
+    )
+  )
+
+  # NOTE: expected final population size is larger than the initial
+  # but identical to the original + added population
+  data <- model_diphtheria_cpp(
+    population = camp_pop,
+    prop_hosp = 0.08,
+    population_change = p,
+    time_end = 200
+  )
+
+  last_value <- aggregate(
+    value ~ demography_group,
+    data = data[data$time == max(data$time), ], FUN = "sum"
+  )
+  last_value
+
+  expect_identical(
+    last_value$value,
+    camp_pop$demography_vector + p$values[[1]],
+    tolerance = 1e-6
+  )
+
+  ## Multiple population changes including decreases
+  p <- list(
+    time = c(70, 80),
+    values = list(
+      c(1e4, 2e5, 1e5), # influx to camp
+      c(-9e3, -1.5e5, -0.5e5) # evacuation from camp
+    )
+  )
+
+  expect_no_condition(
+    model_diphtheria_cpp(
+      population = camp_pop,
+      population_change = p
+    )
+  )
+
+  # NOTE: expected final population size is larger than the initial
+  # but identical to the original + net added population
+  data <- model_diphtheria_cpp(
+    population = camp_pop,
+    prop_hosp = 0.08,
+    population_change = p,
+    time_end = 200
+  )
+
+  last_value <- aggregate(
+    value ~ demography_group,
+    data = data[data$time == max(data$time), ], FUN = "sum"
+  )
+  last_value
+
+  expect_identical(
+    last_value$value,
+    camp_pop$demography_vector + Reduce(x = p$values, f = `+`),
+    tolerance = 1e-6
+  )
+})

From 3074972cf4eb679782489241963b8a11c439bb0f Mon Sep 17 00:00:00 2001
From: pratikunterwegs <pratikgupte16@gmail.com>
Date: Fri, 26 Jan 2024 09:55:25 +0000
Subject: [PATCH 12/18] Misc tests and checks for diphtheria model

---
 R/model_diphtheria.R                   | 6 ++++++
 tests/testthat/test-model_diphtheria.R | 4 ++--
 2 files changed, 8 insertions(+), 2 deletions(-)

diff --git a/R/model_diphtheria.R b/R/model_diphtheria.R
index 26da4c1a..fe234464 100644
--- a/R/model_diphtheria.R
+++ b/R/model_diphtheria.R
@@ -175,6 +175,12 @@ model_diphtheria_cpp <- function(population,
     population_change,
     null.ok = TRUE, len = 2L, types = c("numeric", "list")
   )
+  if (!is.null(population_change)) {
+    checkmate::assert_names(
+      names(population_change),
+      identical.to = c("time", "values")
+    )
+  }
 
   # check the time end and increment
   # restrict increment to lower limit of 1e-6
diff --git a/tests/testthat/test-model_diphtheria.R b/tests/testthat/test-model_diphtheria.R
index b4e7b314..19578469 100644
--- a/tests/testthat/test-model_diphtheria.R
+++ b/tests/testthat/test-model_diphtheria.R
@@ -52,11 +52,11 @@ test_that("Diptheria model, basic expectations", {
   # check for identical numbers of individuals at start and end
   # NOTE: high tolerance because hospitalised compartment is not directly
   # linked to infectious compartment per Finger et al. model structure.
-  # leads to more individuals at final state than initial state
+  # leads to different individuals at final state than initial state
   expect_identical(
     sum(output[output$time == min(output$time), ]$value),
     sum(output[output$time == max(output$time), ]$value),
-    tolerance = 100
+    tolerance = 1e-6
   )
   # check that all age groups in the simulation are the same
   # size as the demography vector

From b08b5244145e23d8d3d6eb0def2cfe1f79d48e61 Mon Sep 17 00:00:00 2001
From: pratikunterwegs <pratikgupte16@gmail.com>
Date: Fri, 26 Jan 2024 13:52:59 +0000
Subject: [PATCH 13/18] Improve diphtheria checks and tests

---
 R/check_args_diphtheria.R              | 18 ++++--
 tests/testthat/test-model_diphtheria.R | 88 ++++++++++++++++++++++++++
 2 files changed, 102 insertions(+), 4 deletions(-)

diff --git a/R/check_args_diphtheria.R b/R/check_args_diphtheria.R
index fcf62875..75143460 100644
--- a/R/check_args_diphtheria.R
+++ b/R/check_args_diphtheria.R
@@ -134,10 +134,20 @@
     )
     invisible(
       lapply(
-        mod_args[["pop_change_values"]], checkmate::assert_numeric,
-        len = length(get_parameter(
-          mod_args[["population"]], "demography_vector"
-        ))
+        mod_args[["pop_change_values"]],
+        FUN = function(x) {
+          stopifnot(
+            "`population_change` `values` must be same length as demography" =
+              checkmate::test_numeric(
+                x,
+                len = length(
+                  get_parameter(
+                    mod_args[["population"]], "demography_vector"
+                  )
+                )
+              )
+          )
+        }
       )
     )
   }
diff --git a/tests/testthat/test-model_diphtheria.R b/tests/testthat/test-model_diphtheria.R
index 19578469..176db9ff 100644
--- a/tests/testthat/test-model_diphtheria.R
+++ b/tests/testthat/test-model_diphtheria.R
@@ -144,3 +144,91 @@ test_that("Diphtheria model with population size changes", {
     tolerance = 1e-6
   )
 })
+
+# Test for poor specification of the population change mechanic
+test_that("Diphtheria model handles population_change errors", {
+  # badly specified pop change
+  p <- list(
+    time = "some time",
+    values = list(
+      c(1, 2, 3)
+    )
+  )
+  expect_error(
+    model_diphtheria_cpp(
+      camp_pop,
+      population_change = p
+    ),
+    regexp = "May only contain the following types: \\{numeric,list\\}"
+  )
+
+  # wrong name of first element
+  p <- list(
+    timesteps = 10,
+    values = list(
+      c(1, 2, 3)
+    )
+  )
+  expect_error(
+    model_diphtheria_cpp(
+      camp_pop,
+      population_change = p
+    ),
+    regexp = "(Names must be identical)*(time)*(values)"
+  )
+
+  # wrong length of values - must be same length as demography groups
+  p <- list(
+    time = 10,
+    values = list(
+      c(1, 2)
+    )
+  )
+  expect_error(
+    model_diphtheria_cpp(
+      camp_pop,
+      population_change = p
+    ),
+    regexp = "`population_change` `values` must be same length as demography"
+  )
+})
+
+# Tests for expected failures/input checks on interventions and time-dependence
+test_that("Diphtheria model input checks", {
+  # expect no contacts interventions allowed
+  expect_error(
+    model_diphtheria_cpp(
+      camp_pop,
+      intervention = list(
+        contacts = no_contacts_intervention(camp_pop) # a dummy intervention
+      )
+    ),
+    regexp = "May only contain the following types: \\{rate_intervention\\}"
+  )
+
+  # expect failure on rate interventions targeting disallowed parameters
+  expect_error(
+    model_diphtheria_cpp(
+      camp_pop,
+      intervention = list(
+        beta = intervention(
+          type = "rate",
+          time_begin = 10, time_end = 20,
+          reduction = 0.1
+        )
+      )
+    ),
+    regexp = "(must be a subset of)*(but has additional elements)"
+  )
+
+  # expect failure on time dependence targeting disallowed parameters
+  expect_error(
+    model_diphtheria_cpp(
+      camp_pop,
+      time_dependence = list(
+        beta = function(time, x) x
+      )
+    ),
+    regexp = "(must be a subset of)*(but has additional elements)"
+  )
+})

From 1a94499b2bbe1b7159596b10276a8ff37bdcca59 Mon Sep 17 00:00:00 2001
From: pratikunterwegs <pratikgupte16@gmail.com>
Date: Fri, 26 Jan 2024 13:53:37 +0000
Subject: [PATCH 14/18] Add diphtheria model vignette, fixes #157

---
 _pkgdown.yml             |   1 +
 vignettes/diphtheria.Rmd | 199 +++++++++++++++++++++++++++++++++++++++
 2 files changed, 200 insertions(+)
 create mode 100644 vignettes/diphtheria.Rmd

diff --git a/_pkgdown.yml b/_pkgdown.yml
index ab7c5c8f..29d96ef6 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -28,6 +28,7 @@ articles:
   contents:
   - vacamole
   - ebola_model
+  - diphtheria
   - benchmarking
 
 reference:
diff --git a/vignettes/diphtheria.Rmd b/vignettes/diphtheria.Rmd
new file mode 100644
index 00000000..d0622d1a
--- /dev/null
+++ b/vignettes/diphtheria.Rmd
@@ -0,0 +1,199 @@
+---
+title: "Modelling a diphtheria outbreak in a humanitarian camp setting"
+output:
+  bookdown::html_vignette2:
+    fig_caption: yes
+    code_folding: show
+pkgdown:
+  as_is: true
+bibliography: references.json
+link-citations: true
+vignette: >
+  %\VignetteIndexEntry{Modelling a diphtheria outbreak in a humanitarian camp setting}
+  %\VignetteEncoding{UTF-8}
+  %\VignetteEngine{knitr::rmarkdown}
+editor_options: 
+  chunk_output_type: console
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>",
+  message = FALSE,
+  warning = FALSE,
+  fig.width = 5,
+  fig.height = 4,
+  dpi = 300
+)
+```
+
+::: {.alert .alert-warning}
+**New to _epidemics_?** It may help to read the ["Get started"](epidemics.html) vignette first!
+:::
+
+This vignette shows how to model an outbreak of diphtheria (or other diseases) within the setting of a humanitarian aid camp, where the camp population may fluctuate due to external factors, including influxes from the affected area and evacuations to other camps or to other areas.
+In such situations, implementing large-scale public health measures such as non-pharmaceutical interventions on social contacts or mass vaccination may be impractical.
+
+_epidemics_ provides a simple SEIHR compartmental model based on @finger2019, in which it is possible to vary the population of each demographic group throughout the model's simulation time.
+This model only tracks infections by demographic groups, and does not include demographic variation in contacts, as this is assumed to be less important in the camp setting.
+The model also does not allow interventions on social contacts, and does not include a 'vaccinated' compartment.
+
+However, the model does allow for seasonlity in model parameters and interventions on model parameters.
+Similarly, the model allows for a proportion of the initial camp population to be considered vaccinated and thus immune from infection.
+
+```{r setup}
+library(epidemics)
+library(dplyr)
+library(ggplot2)
+```
+
+## Modelling an outbreak with pre-existing immunity
+
+We create a population object corresponding to the Kutupalong camp in Bangladesh, rounded to the nearest 100, as described in Additional file 1 provided with @finger2019.
+This population has three age groups, < 5 years, 5 -- 14 years, $\qeq$ 15 years.
+We assume that only one individual is infectious in each age group.
+
+```{r camp-population}
+# three age groups with five compartments SEIHR
+n_age_groups <- 3
+demography_vector <- c(83000, 108200, 224600)
+initial_conditions <- matrix(0, nrow = n_age_groups, ncol = 5)
+
+# 1 individual in each group is infectious
+initial_conditions[, 1] <- demography_vector - 1
+initial_conditions[, 3] <- rep(1, n_age_groups)
+
+# camp social contacts are assumed to be uniform
+camp_pop <- population(
+  contact_matrix = matrix(1, nrow = n_age_groups, ncol = n_age_groups),
+  demography_vector = demography_vector,
+  initial_conditions = initial_conditions / demography_vector
+)
+
+camp_pop
+```
+
+We assume, following @finger2019, that 20% of the 5 -- 14 year-olds are vaccinated against and immune to diphtheria, but that coverage is much lower among other age groups.
+
+```{r}
+# 20% of 5-14 year-olds are vaccinated
+prop_vaccinated <- c(0.05, 0.2, 0.05)
+```
+
+We run the model with its default parameters, assuming that:
+- diphtheria has an $\R_0$ of 4.0 and a mean infectious period of 4.5 days, giving a transmissibility ($\beta$) of about 0.889;
+- diphtheria has a pre-infectious or incubation period of 3 days, giving an infectiousness rate ($\sigma$) of about 0.33; and
+- the recovery rate of diphtheria is about 0.33.
+
+We also stick to default assumptions regarding: case reporting (that about 3% of cases are reported), the proportion of reported needing hospitalisation (1%), the time taken by cases needing hospitalisation to seek and be admitted to hospital (5 days, giving a hospitalisation rate of 0.2), and the time spent in hospital (5 days, giving a hospitalisation recovery rate of 0.2).
+
+We assume there are no interventions or seasonal effects, and run the model and plot the outcomes.
+
+```{r diphtheria-basic}
+data <- model_diphtheria_cpp(
+  population = camp_pop,
+  prop_vaccinated = prop_vaccinated
+)
+```
+
+```{r class.source = 'fold-hide', fig.cap="Model results from a single run showing the number of individuals infectious with diphtheria over 100 days of the outbreak."}
+filter(data, compartment == "infectious") |>
+  ggplot() +
+  geom_line(
+    aes(time, value, colour = demography_group)
+  ) +
+  scale_y_continuous(
+    labels = scales::comma
+  ) +
+  scale_colour_brewer(
+    palette = "Dark2",
+    name = "Age group",
+    labels = c("<5", "5-15", ">15")
+  ) +
+  expand_limits(
+    x = c(0, 101)
+  ) +
+  theme_bw() +
+  theme(
+    legend.position = "top"
+  ) +
+  labs(
+    x = "Simulation time (days)",
+    y = "Individuals infectious"
+  )
+```
+
+## Modelling an outbreak with changing population sizes
+
+We now model the same outbreak, but with fluctuation in susceptible individuals towards the end of the outbreak, to illustrate the effect this could have on outbreak dynamics.
+In this example, we do not assume any prior immunity.
+
+We prepare a population change schedule as a named list giving the times of each change, and the corresponding changes to each demographic group.
+
+Note that the model assumes that these changes apply _only to the susceptible compartment_, as we assume that the wider population entering the camp is not yet affected by diphtheria, and that infected or hospitalised individuals do not leave the camp.
+
+```{r pop_change}
+# susceptibles increase by about 12%, 92%, and 89% of initial sizes
+pop_change <- list(
+  time = 70,
+  values = list(
+    c(1e4, 1e5, 2e5)
+  )
+)
+```
+
+Here, the number of individuals increases by about 75% overall, which is well within reported values in real world scenarios [@finger2019].
+
+```{r}
+data <- model_diphtheria_cpp(
+  population = camp_pop,
+  population_change = pop_change
+)
+```
+
+```{r class.source = 'fold-hide', fig.cap="Model results from a single run showing the number of individuals infectious with diphtheria over 100 days of the outbreak, with an increase in the camp population size."}
+filter(data, compartment == "infectious") |>
+  ggplot() +
+  geom_line(
+    aes(time, value, colour = demography_group)
+  ) +
+  geom_vline(
+    xintercept = pop_change$time,
+    colour = "red",
+    linetype = "dashed",
+    linewidth = 0.2
+  ) +
+  annotate(
+    geom = "text",
+    x = pop_change$time,
+    y = 20e3,
+    label = "Population increase",
+    angle = 90,
+    vjust = "inward",
+    colour = "red"
+  ) +
+  scale_y_continuous(
+    labels = scales::comma
+  ) +
+  scale_colour_brewer(
+    palette = "Dark2",
+    name = "Age group",
+    labels = c("<5", "5-15", ">15")
+  ) +
+  expand_limits(
+    x = c(0, 101)
+  ) +
+  theme_bw() +
+  theme(
+    legend.position = "top"
+  ) +
+  labs(
+    x = "Simulation time (days)",
+    y = "Individuals infectious"
+  )
+```
+
+This example shows how an increase in the number of susceptibles in the population can lead to a rise in the number of infections and extend the duration of an outbreak.
+
+## References

From ac6ff571ec6515bd2c9ac6c4084fd8f03d7167e5 Mon Sep 17 00:00:00 2001
From: pratikunterwegs <pratikgupte16@gmail.com>
Date: Fri, 26 Jan 2024 13:53:51 +0000
Subject: [PATCH 15/18] Fix mistake in ebola vignette text

---
 vignettes/ebola_model.Rmd | 1 -
 1 file changed, 1 deletion(-)

diff --git a/vignettes/ebola_model.Rmd b/vignettes/ebola_model.Rmd
index 5a59520e..eae690ea 100644
--- a/vignettes/ebola_model.Rmd
+++ b/vignettes/ebola_model.Rmd
@@ -41,7 +41,6 @@ The transitions between compartments are based on an Erlang model developed by @
 The model does not include demographic variation in contacts, as ebola spreads primarily among contacts caring for infected individuals, making age or demographic structure less important.
 
 This model can currently accommodate interventions on model rates only (as there are no demographic groups for contacts interventions).
-The only rate in the model that can be influenced by an intervention is the transmission rate $\beta$.
 
 ```{r setup}
 # some initial setup to load necessary packages

From 218aecbaf2c8dae2c50003d1efa8f153a34e5780 Mon Sep 17 00:00:00 2001
From: pratikunterwegs <pratikgupte16@gmail.com>
Date: Fri, 26 Jan 2024 14:09:21 +0000
Subject: [PATCH 16/18] Update WORDLIST

---
 inst/WORDLIST            | 3 +++
 vignettes/diphtheria.Rmd | 4 ++--
 2 files changed, 5 insertions(+), 2 deletions(-)

diff --git a/inst/WORDLIST b/inst/WORDLIST
index ec21b8db..438f4acd 100644
--- a/inst/WORDLIST
+++ b/inst/WORDLIST
@@ -30,6 +30,7 @@ Katriona
 Klinkenberg
 Korthals
 Krzywinski
+Kutupalong
 LMIC
 Leeuwen
 Leung
@@ -80,6 +81,7 @@ epiparameter
 epirecipes
 etc
 finalsize
+geq
 ggplot
 gh
 github
@@ -88,6 +90,7 @@ kylieainslie
 md
 odeint
 odin
+olds
 org
 packagename
 pkg
diff --git a/vignettes/diphtheria.Rmd b/vignettes/diphtheria.Rmd
index d0622d1a..775e09f7 100644
--- a/vignettes/diphtheria.Rmd
+++ b/vignettes/diphtheria.Rmd
@@ -39,7 +39,7 @@ _epidemics_ provides a simple SEIHR compartmental model based on @finger2019, in
 This model only tracks infections by demographic groups, and does not include demographic variation in contacts, as this is assumed to be less important in the camp setting.
 The model also does not allow interventions on social contacts, and does not include a 'vaccinated' compartment.
 
-However, the model does allow for seasonlity in model parameters and interventions on model parameters.
+However, the model does allow for seasonality in model parameters and interventions on model parameters.
 Similarly, the model allows for a proportion of the initial camp population to be considered vaccinated and thus immune from infection.
 
 ```{r setup}
@@ -51,7 +51,7 @@ library(ggplot2)
 ## Modelling an outbreak with pre-existing immunity
 
 We create a population object corresponding to the Kutupalong camp in Bangladesh, rounded to the nearest 100, as described in Additional file 1 provided with @finger2019.
-This population has three age groups, < 5 years, 5 -- 14 years, $\qeq$ 15 years.
+This population has three age groups, < 5 years, 5 -- 14 years, $\geq$ 15 years.
 We assume that only one individual is infectious in each age group.
 
 ```{r camp-population}

From c179c936842e6395f2a52254c5cc8bc365a11d91 Mon Sep 17 00:00:00 2001
From: Pratik Gupte <pratik.gupte@lshtm.ac.uk>
Date: Mon, 5 Feb 2024 13:00:16 +0000
Subject: [PATCH 17/18] Edits to vignettes from code review

Co-authored-by: Adam Kucharski <adam.kucharski@lshtm.ac.uk>
---
 vignettes/diphtheria.Rmd | 32 ++++++++++++++++----------------
 1 file changed, 16 insertions(+), 16 deletions(-)

diff --git a/vignettes/diphtheria.Rmd b/vignettes/diphtheria.Rmd
index 775e09f7..400bd327 100644
--- a/vignettes/diphtheria.Rmd
+++ b/vignettes/diphtheria.Rmd
@@ -32,12 +32,12 @@ knitr::opts_chunk$set(
 **New to _epidemics_?** It may help to read the ["Get started"](epidemics.html) vignette first!
 :::
 
-This vignette shows how to model an outbreak of diphtheria (or other diseases) within the setting of a humanitarian aid camp, where the camp population may fluctuate due to external factors, including influxes from the affected area and evacuations to other camps or to other areas.
-In such situations, implementing large-scale public health measures such as non-pharmaceutical interventions on social contacts or mass vaccination may be impractical.
+This vignette shows how to model an outbreak of diphtheria (or a similar acute directly-transmitted infectious disease) within the setting of a humanitarian aid camp, where the camp population may fluctuate due to external factors, including influxes from crisis affected areas or evacuations to other camps or areas.
+In such situations, implementing large-scale public health measures such contact tracing and quarantine, or introducing reactive mass vaccination may be challenging.
 
-_epidemics_ provides a simple SEIHR compartmental model based on @finger2019, in which it is possible to vary the population of each demographic group throughout the model's simulation time.
-This model only tracks infections by demographic groups, and does not include demographic variation in contacts, as this is assumed to be less important in the camp setting.
-The model also does not allow interventions on social contacts, and does not include a 'vaccinated' compartment.
+_epidemics_ provides a simple SEIHR compartmental model based on @finger2019, in which it is possible to vary the population of each demographic group throughout the model's simulation time and explore the resulting epidemic dynamics.
+This baseline model only tracks infections by demographic groups, and does not include variation in contacts between demographic groups (e.g. by age or occupation), as contacts are likely to be less clearly stratified in a camp setting.
+The baseline model also does not allow interventions that target social contacts (such as social distancing or quarantine), and does not include a 'vaccinated' compartment. It is therefore suited to analysis of a rapidly spreading infection prior to the introduction of any reactive vaccine campaign.
 
 However, the model does allow for seasonality in model parameters and interventions on model parameters.
 Similarly, the model allows for a proportion of the initial camp population to be considered vaccinated and thus immune from infection.
@@ -50,7 +50,7 @@ library(ggplot2)
 
 ## Modelling an outbreak with pre-existing immunity
 
-We create a population object corresponding to the Kutupalong camp in Bangladesh, rounded to the nearest 100, as described in Additional file 1 provided with @finger2019.
+We create a population object corresponding to the Kutupalong camp in Cox's Bazar, Bangladesh in 2017-18, rounded to the nearest 100, as described in Additional file 1 provided with @finger2019.
 This population has three age groups, < 5 years, 5 -- 14 years, $\geq$ 15 years.
 We assume that only one individual is infectious in each age group.
 
@@ -64,7 +64,7 @@ initial_conditions <- matrix(0, nrow = n_age_groups, ncol = 5)
 initial_conditions[, 1] <- demography_vector - 1
 initial_conditions[, 3] <- rep(1, n_age_groups)
 
-# camp social contacts are assumed to be uniform
+# camp social contact rates are assumed to be uniform within and between age groups
 camp_pop <- population(
   contact_matrix = matrix(1, nrow = n_age_groups, ncol = n_age_groups),
   demography_vector = demography_vector,
@@ -74,7 +74,7 @@ camp_pop <- population(
 camp_pop
 ```
 
-We assume, following @finger2019, that 20% of the 5 -- 14 year-olds are vaccinated against and immune to diphtheria, but that coverage is much lower among other age groups.
+We assume, following @finger2019, that 20% of the 5 -- 14 year-olds were vaccinated against (and immune to) diphtheria prior to the start of the outbreak, but that coverage is much lower among other age groups.
 
 ```{r}
 # 20% of 5-14 year-olds are vaccinated
@@ -82,13 +82,13 @@ prop_vaccinated <- c(0.05, 0.2, 0.05)
 ```
 
 We run the model with its default parameters, assuming that:
-- diphtheria has an $\R_0$ of 4.0 and a mean infectious period of 4.5 days, giving a transmissibility ($\beta$) of about 0.889;
+- diphtheria has an $R_0$ of 4.0 and a mean infectious period of 4.5 days, giving a transmissibility ($\beta$) of about 0.889;
 - diphtheria has a pre-infectious or incubation period of 3 days, giving an infectiousness rate ($\sigma$) of about 0.33; and
 - the recovery rate of diphtheria is about 0.33.
 
-We also stick to default assumptions regarding: case reporting (that about 3% of cases are reported), the proportion of reported needing hospitalisation (1%), the time taken by cases needing hospitalisation to seek and be admitted to hospital (5 days, giving a hospitalisation rate of 0.2), and the time spent in hospital (5 days, giving a hospitalisation recovery rate of 0.2).
+We also make several assumptions regarding clinical progression and reporting. Specifically: case reporting (that about 3% of infections are reported as cases), the proportion of reported cases needing hospitalisation (1%), the time taken by cases needing hospitalisation to seek and be admitted to hospital (5 days, giving a daily hospitalisation rate of 0.2 among cases), and time spent in hospital (5 days, giving a daily hospitalisation recovery rate of 0.2).
 
-We assume there are no interventions or seasonal effects, and run the model and plot the outcomes.
+Finally, we assume there are no interventions or seasonal effects that affect the dynamics of transmission during the outbreak. We then run the model and plot the outcomes.
 
 ```{r diphtheria-basic}
 data <- model_diphtheria_cpp(
@@ -126,12 +126,12 @@ filter(data, compartment == "infectious") |>
 
 ## Modelling an outbreak with changing population sizes
 
-We now model the same outbreak, but with fluctuation in susceptible individuals towards the end of the outbreak, to illustrate the effect this could have on outbreak dynamics.
-In this example, we do not assume any prior immunity.
+We now model the same outbreak, but an increase in susceptible individuals towards the end of the outbreak, to illustrate the effect that an influx of non-immune individuals could could have on outbreak dynamics.
+In this example, we do not assume any prior immunity among new arrivals into the population.
 
 We prepare a population change schedule as a named list giving the times of each change, and the corresponding changes to each demographic group.
 
-Note that the model assumes that these changes apply _only to the susceptible compartment_, as we assume that the wider population entering the camp is not yet affected by diphtheria, and that infected or hospitalised individuals do not leave the camp.
+Note that the model assumes that these changes apply _only to the susceptible compartment_, as we assume that the wider population entering the camp is not yet affected by diphtheria (i.e. no infected or recovered arrival), and that already infected or hospitalised individuals do not leave the camp.
 
 ```{r pop_change}
 # susceptibles increase by about 12%, 92%, and 89% of initial sizes
@@ -143,7 +143,7 @@ pop_change <- list(
 )
 ```
 
-Here, the number of individuals increases by about 75% overall, which is well within reported values in real world scenarios [@finger2019].
+Here, the population size of the camp increases by about 75% overall, which is similar to reported values in the Kutupalong camp scenario [@finger2019].
 
 ```{r}
 data <- model_diphtheria_cpp(
@@ -194,6 +194,6 @@ filter(data, compartment == "infectious") |>
   )
 ```
 
-This example shows how an increase in the number of susceptibles in the population can lead to a rise in the number of infections and extend the duration of an outbreak.
+This example shows how an increase in the number of susceptibles in the population can lead to a rise in the transmission potential of the infection and therefore extend the duration of an outbreak.
 
 ## References

From dc227b69713ba9d893863db3b51e927ee2f34999 Mon Sep 17 00:00:00 2001
From: pratikunterwegs <pratikgupte16@gmail.com>
Date: Mon, 5 Feb 2024 15:17:32 +0000
Subject: [PATCH 18/18] Update vignette and wordlist

---
 inst/WORDLIST            |  1 +
 vignettes/diphtheria.Rmd | 36 +++++++++++++++++++++++++++++-------
 2 files changed, 30 insertions(+), 7 deletions(-)

diff --git a/inst/WORDLIST b/inst/WORDLIST
index 438f4acd..5283ee33 100644
--- a/inst/WORDLIST
+++ b/inst/WORDLIST
@@ -4,6 +4,7 @@ Altman
 Analogs
 BH
 BMC
+Bazar
 Bjørnstad
 CMD
 COVID
diff --git a/vignettes/diphtheria.Rmd b/vignettes/diphtheria.Rmd
index 400bd327..5f538d66 100644
--- a/vignettes/diphtheria.Rmd
+++ b/vignettes/diphtheria.Rmd
@@ -24,7 +24,7 @@ knitr::opts_chunk$set(
   warning = FALSE,
   fig.width = 5,
   fig.height = 4,
-  dpi = 300
+  dpi = 150
 )
 ```
 
@@ -35,6 +35,10 @@ knitr::opts_chunk$set(
 This vignette shows how to model an outbreak of diphtheria (or a similar acute directly-transmitted infectious disease) within the setting of a humanitarian aid camp, where the camp population may fluctuate due to external factors, including influxes from crisis affected areas or evacuations to other camps or areas.
 In such situations, implementing large-scale public health measures such contact tracing and quarantine, or introducing reactive mass vaccination may be challenging.
 
+::: {.alert .alert-info}
+The [vignette on modelling a vaccination campaign](modelling_vaccination.html) shows how to model the introduction of a mass vaccination campaign for a fixed, stable population.
+::: 
+
 _epidemics_ provides a simple SEIHR compartmental model based on @finger2019, in which it is possible to vary the population of each demographic group throughout the model's simulation time and explore the resulting epidemic dynamics.
 This baseline model only tracks infections by demographic groups, and does not include variation in contacts between demographic groups (e.g. by age or occupation), as contacts are likely to be less clearly stratified in a camp setting.
 The baseline model also does not allow interventions that target social contacts (such as social distancing or quarantine), and does not include a 'vaccinated' compartment. It is therefore suited to analysis of a rapidly spreading infection prior to the introduction of any reactive vaccine campaign.
@@ -64,7 +68,8 @@ initial_conditions <- matrix(0, nrow = n_age_groups, ncol = 5)
 initial_conditions[, 1] <- demography_vector - 1
 initial_conditions[, 3] <- rep(1, n_age_groups)
 
-# camp social contact rates are assumed to be uniform within and between age groups
+# camp social contact rates are assumed to be uniform within and between
+# age groups
 camp_pop <- population(
   contact_matrix = matrix(1, nrow = n_age_groups, ncol = n_age_groups),
   demography_vector = demography_vector,
@@ -82,6 +87,7 @@ prop_vaccinated <- c(0.05, 0.2, 0.05)
 ```
 
 We run the model with its default parameters, assuming that:
+
 - diphtheria has an $R_0$ of 4.0 and a mean infectious period of 4.5 days, giving a transmissibility ($\beta$) of about 0.889;
 - diphtheria has a pre-infectious or incubation period of 3 days, giving an infectiousness rate ($\sigma$) of about 0.33; and
 - the recovery rate of diphtheria is about 0.33.
@@ -150,12 +156,24 @@ data <- model_diphtheria_cpp(
   population = camp_pop,
   population_change = pop_change
 )
+
+# summarise population change in susceptibles
+data_pop_size <- filter(data, compartment == "susceptible") |>
+  group_by(time) |>
+  summarise(
+    total_susceptibles = sum(value)
+  )
 ```
 
-```{r class.source = 'fold-hide', fig.cap="Model results from a single run showing the number of individuals infectious with diphtheria over 100 days of the outbreak, with an increase in the camp population size."}
-filter(data, compartment == "infectious") |>
-  ggplot() +
+```{r class.source = 'fold-hide', fig.cap="Model results from a single run showing the number of individuals infectious with diphtheria over 100 days of the outbreak, with an increase in the camp population size. Shaded blue region shows the number of individuals susceptible to infection (right-hand side Y axis)."}
+ggplot() +
+  geom_area(
+    data = data_pop_size,
+    aes(time, total_susceptibles / 10),
+    fill = "steelblue", alpha = 0.5
+  ) +
   geom_line(
+    data = filter(data, compartment == "infectious"),
     aes(time, value, colour = demography_group)
   ) +
   geom_vline(
@@ -174,11 +192,15 @@ filter(data, compartment == "infectious") |>
     colour = "red"
   ) +
   scale_y_continuous(
-    labels = scales::comma
+    labels = scales::comma,
+    sec.axis = dup_axis(
+      trans = function(x) x * 10,
+      name = "Individuals susceptible"
+    )
   ) +
   scale_colour_brewer(
     palette = "Dark2",
-    name = "Age group",
+    name = "Age group (individuals infectious)",
     labels = c("<5", "5-15", ">15")
   ) +
   expand_limits(