BDI-pathogens · Jun 9, 2021
diff --git a/‎DESCRIPTION
+2 b/‎DESCRIPTION
+2
diff --git a/‎Makefile
+1-1 b/‎Makefile
+1-1
diff --git a/‎NAMESPACE
+2 b/‎NAMESPACE
+2
diff --git a/‎R/Model.R
+159-42 b/‎R/Model.R
+159-42
diff --git a/‎R/Network.R
+1-1 b/‎R/Network.R
+1-1
diff --git a/‎R/Parameters.R
+1-1 b/‎R/Parameters.R
+1-1
diff --git a/‎R/Strain.R
+79 b/‎R/Strain.R
+79
@@ -59,7 +59,9 @@ Collate:
   Network.R
   Model.R
   Simulation.R
+  Strain.R
   VaccineSchedule.R
+  Vaccine.R
 Suggests:
     testthat
 RoxygenNote: 7.1.1
 
@@ -47,7 +47,7 @@ CFLAGS = -g -Wall -fmessage-length=0 -I$(INC) $(shell gsl-config --cflags) -O0
 LDFLAGS = -L$(LIB) $(shell gsl-config --libs)
 
 # Swig's input
-SWIG_INPUT = src/disease.h src/ward.h src/nurse.h src/network_utils.i src/input.h src/individual.h src/hospital.h src/params.h src/structure.h src/constant.h src/doctor.h src/utilities.h src/model_utils.i src/covid19.i src/list.h src/network.h src/model.h src/interventions.h src/params_utils.i src/demographics.h src/strain.h
+SWIG_INPUT = src/disease.h src/ward.h src/nurse.h src/network_utils.i src/vaccine_utils.i src/strain_utils.i src/input.h src/individual.h src/hospital.h src/params.h src/structure.h src/constant.h src/doctor.h src/utilities.h src/model_utils.i src/covid19.i src/list.h src/network.h src/model.h src/interventions.h src/params_utils.i src/demographics.h src/strain.h
 
 # Swig's output
 SWIG_OUTPUT_PY = src/covid19_wrap.o src/covid19_wrap.c src/covid19.py src/_covid19.cpython-37m-darwin.so src/build src/covid19.egg-info
 
@@ -14,8 +14,10 @@ export("Network")
 export("Environment")
 export("Agent")
 export("Simulation")
+export("Strain")
 export("COVID19IBM")
 export("VaccineSchedule")
+export("Vaccine")
 
 # wrapper function for R6 calsses
 export("Parameters.new")
 
@@ -23,6 +23,7 @@ SWIG_calculate_R_instanteous <- calculate_R_instanteous
 SWIG_seed_infect_by_idx <- seed_infect_by_idx
 SWIG_add_new_strain <- add_new_strain
 SWIG_destroy_model <- destroy_model
+SWIG_set_cross_immunity_probability <- set_cross_immunity_probability
 
 
 #' R6Class Model
@@ -225,8 +226,13 @@ Model <- R6Class( classname = 'Model', cloneable = FALSE,
       enum <- get_base_param_from_enum(param)
       if (!is.null(enum)) {
         # multi-value parameter (C array)
-        getter <- get(paste0("get_model_param_", enum$base_name))
-        result <- getter( self$c_model, enum$index )
+        getter_str = paste0("get_model_param_", enum$base_name)
+        if (exists(getter_str)) {
+          getter <- get(paste0("get_model_param_", enum$base_name))
+          result <- getter( self$c_model, enum$index )
+        } else {
+          result = private$params_object$get_param(param)
+        }
       } else {
         # single-value parameter
         getter_str <- paste0("get_model_param_", param)
@@ -465,26 +471,81 @@ Model <- R6Class( classname = 'Model', cloneable = FALSE,
     #' @param strain_idx The idx of the strain the person is infected with
     #' @param network_id The network ID.
     #' @return \code{TRUE} on success, \code{FALSE} otherwise.
-    seed_infect_by_idx = function(ID, strain_idx = 1, network_id = -1 )
+    seed_infect_by_idx = function(ID, strain_idx = 0, strain = NULL, network_id = -1 )
     {
       n_total <- private$c_params$n_total
-
       if (ID < 0 || ID >= n_total) {
         stop("ID out of range (0<=ID<n_total)")
       }
-      
+
+      if( !is.null( strain ) )
+      {
+        if (!is.R6(strain) || !('Strain' %in% class(strain)))
+          stop("argument strain must be an object of type Strain")
+
+        strain_idx = strain$idx()
+      }
+
+      n_strains = self$c_model$n_initialised_strains;
+      if( strain_idx < 0  || strain_idx >= n_strains )
+        stop( "strain_idx out of range (0 <= strain_idx < self$c_model$n_initialized_strains)" )
+
       res <- SWIG_seed_infect_by_idx(self$c_model, ID, strain_idx, network_id)
       return(as.logical(res))
     },
 
     #' @description Adds a new strain (variant)
     #' Wrapper for C API \code{add_new_strain}.
     #' @param transmission_multiplier The relative transmission rate of the strain
-    #' @return \code{strain_idx} The index assigned to the new strain
-    add_new_strain = function(ID, transmission_multiplier = 1 )
+    #' @param hospitalised_fraction the fraction of symptomatic (not mild) who progress to hospital [default: None is no change)]
+    #' @return \code{Strain} A Strain object representing this strain
+    add_new_strain = function( transmission_multiplier = 1, hospitalised_fraction = NA )
+    {
+
+      max_n_strains = self$get_param( "max_n_strains" )
+      n_strains = self$c_model$n_initialised_strains;
+
+      if( n_strains == max_n_strains )
+        stop( "cannot add any more strains - increase the parameter max_n_strains at the initialisation of the model" )
+
+      if( is.na( hospitalised_fraction ) )
+      {
+        hospitalised_fraction = c()
+        for( idx in 1:length( AgeGroupEnum ) )
+          hospitalised_fraction[ idx ] = self$get_param(
+            sprintf( "hospitalised_fraction%s", names( AgeGroupEnum[idx])) )
+      }
+
+      c_model_ptr <- private$c_model_ptr()
+      strain_idx<-.Call('R_add_new_strain',c_model_ptr,transmission_multiplier,
+                        hospitalised_fraction, PACKAGE='OpenABMCovid19');
+
+      return( Strain$new( self, strain_idx ) )
+    },
+
+    #' @description Set the cross_immunity matrix
+    #' Wrapper for C API \code{set_cross_immunity_probability}.
+    #' @param cross_immunity the cross immunity matrix
+    #' @return
+
+    set_cross_immunity_matrix = function( cross_immunity )
     {
-      strain_idx <- add_new_strain(self$c_model, transmission_multiplier )
-      return(strain_idx)
+        max_n_strains = self$get_param("max_n_strains")
+
+        if( !is.matrix( cross_immunity ) )
+          stop( "cross_immunity is a matrix ")
+
+        if( max( nrow( cross_immunity ), ncol( cross_immunity ) ) > max_n_strains )
+          stop( "cross_immunnity matrix is of maximum size max_n_strains")
+
+        for( i in 1:nrow( cross_immunity) )
+          for( j in 1:ncol( cross_immunity) )
+          {
+            if( cross_immunity[i,j] < 0 || cross_immunity[i,j] > 1 )
+              stop( "cross_immunity must be between 0 and 1")
+
+            SWIG_set_cross_immunity_probability( self$c_model, i-1, j-1, cross_immunity[i,j] )
+          }
     },
 
     #' @description Get the list of network IDs
@@ -535,43 +596,90 @@ Model <- R6Class( classname = 'Model', cloneable = FALSE,
       return(as.data.frame(tmp))
     },
 
-    #' @description Vaccinate an individual.
-    #' Wrapper for C API \code{intervention_vaccinate_by_idx}.
-    #' @param ID The ID of the individual (must be \code{0 <= ID <= n_total}).
-    #' @param vaccine_type The type of vaccine, see \code{\link{VACCINE_TYPES}}.
-    #' @param efficacy Probability that the person is successfully vaccinated
+    #' @description Add a new vaccine.
+    #' Wrapper for C API \code{add_vaccine}.
+    #' @param full_efficacy Probability that the person is successfully vaccinated
     #'   (must be \code{0 <= efficacy <= 1}).
+    #' @param symptoms_efficacy Probability that the person is successfully vaccinated
+    #'   against getting symptoms (must be \code{0 <= efficacy <= 1}).
+    #' @param sever_efficacy Probability that the person is successfully vaccinated
+    #'   against getting severer symptoms (must be \code{0 <= efficacy <= 1}).
     #' @param time_to_protect Delay before it takes effect (in days).
     #' @param vaccine_protection_period The duration of the vaccine before it
     #'   wanes.
+    #' @return Vaccine object
+    add_vaccine = function(
+      full_efficacy     = 1.0,
+      symptoms_efficacy = 1.0,
+      severe_efficacy   = 1.0,
+      time_to_protect   = 14,
+      vaccine_protection_period = 1000
+    )
+    {
+      if( time_to_protect < 1 )
+        stop( "vaccine must take at least one day to take effect" )
+
+      if( vaccine_protection_period <= time_to_protect )
+        stop( "vaccine must protect for longer than it takes to by effective" )
+
+      n_strains = self$get_param( "max_n_strains" )
+
+      if( length( full_efficacy ) == 1 )
+        full_efficacy = rep( full_efficacy,  n_strains )
+
+      if( length( full_efficacy ) != n_strains )
+        stop( "full_efficacy must be a float or a list of length max_n_strains" )
+
+      if( length( symptoms_efficacy ) == 1 )
+        symptoms_efficacy = rep( symptoms_efficacy,  n_strains )
+
+      if( length( symptoms_efficacy ) != n_strains )
+        stop( "symptoms_efficacy must be a float or a list of length max_n_strains" )
+
+      if( length( severe_efficacy ) == 1 )
+        severe_efficacy = rep( severe_efficacy,  n_strains )
+
+      if( length( severe_efficacy ) != n_strains )
+        stop( "severe_efficacy must be a float or a list of length max_n_strains" )
+
+      for( idx in 1:n_strains )
+      {
+        if( full_efficacy[ idx ] < 0  | full_efficacy[ idx ] > 1 )
+          stop( "full_efficacy must be between 0 and 1" )
+
+        if( symptoms_efficacy[ idx ] < 0  | symptoms_efficacy[ idx ] > 1 )
+          stop( "symptoms_efficacy must be between 0 and 1" )
+
+        if( severe_efficacy[ idx ] < 0  | severe_efficacy[ idx ] > 1 )
+          stop( "severe_efficacy must be between 0 and 1" )
+      }
+
+      c_model_ptr <- private$c_model_ptr()
+      idx<-.Call('R_add_vaccine', c_model_ptr, full_efficacy, symptoms_efficacy,
+                 severe_efficacy, time_to_protect, vaccine_protection_period,
+                 PACKAGE='OpenABMCovid19');
+
+      return( Vaccine$new( self, idx ) )
+    },
+
+    #' @description Vaccinate an individual.
+    #' Wrapper for C API \code{intervention_vaccinate_by_idx}.
+    #' @param ID The ID of the individual (must be \code{0 <= ID <= n_total}).
+    #' @param vaccine The of vaccine object to be given
     #' @return Logical value, \code{TRUE} if vaccinated \code{FALSE} otherwise.
-    vaccinate_individual = function(
-      ID,
-      vaccine_type = 0,
-      efficacy = 1.0,
-      time_to_protect = 14,
-      vaccine_protection_period = 1000 )
+    vaccinate_individual = function( ID, vaccine )
     {
       n_total <- private$c_params$n_total
 
       if (ID < 0 || ID >= n_total) {
         stop("ID out of range (0<=ID<n_total)")
       }
-      if (efficacy < 0 || efficacy > 1) {
-        stop("efficacy must be between 0 and 1")
-      }
-      if (time_to_protect < 1) {
-        stop("vaccine must take at least one day to take effect")
-      }
-      if (vaccine_protection_period <= time_to_protect) {
-        stop("vaccine must protect for longer than it takes to by effective")
-      }
-      if (!is.numeric(vaccine_type) || ! vaccine_type %in% VACCINE_TYPES) {
-        stop("vaccine type must be listed in VaccineTypesEnum")
-      }
 
-      res <- intervention_vaccinate_by_idx(self$c_model, ID, vaccine_type,
-        efficacy, time_to_protect, vaccine_protection_period)
+      if (!is.R6(vaccine) || !('Vaccine' %in% class(vaccine)))
+        stop("argument vaccine must be an object of type Vaccine")
+
+      res <- intervention_vaccinate_by_idx(self$c_model, ID, vaccine$c_vaccine )
+
       return(as.logical(res))
     },
 
@@ -584,14 +692,23 @@ Model <- R6Class( classname = 'Model', cloneable = FALSE,
       if (!is.R6(schedule) || !('VaccineSchedule' %in% class(schedule))) {
         stop("argument VaccineSchedule must be an object of type VaccineSchedule")
       }
-      return(as.logical(intervention_vaccinate_age_group(
-        self$c_model,
-        schedule$fraction_to_vaccinate,
-        schedule$vaccine_type,
-        schedule$efficacy,
-        schedule$time_to_protect,
-        schedule$vaccine_protection_period,
-        schedule$total_vaccinated)))
+
+      vaccine = schedule$vaccine
+      if (!is.R6(vaccine) || !('Vaccine' %in% class(vaccine))) {
+        stop("argument schedule$vaccine must be an object of type Vaccine")
+      }
+
+      c_model_ptr    = private$c_model_ptr()
+      c_vaccine_ptr  = vaccine$c_vaccine@ref
+      total_vaccinated = schedule$total_vaccinated
+      fraction_to_vaccinate = schedule$fraction_to_vaccinate
+
+      new_vaccinated <-.Call( "R_intervention_vaccinate_age_group",
+            c_model_ptr, fraction_to_vaccinate, c_vaccine_ptr )
+
+      schedule$total_vaccinated = total_vaccinated + new_vaccinated
+
+      return( sum( new_vaccinated ) )
     },
 
     #' @description Gets information about all individuals. Wrapper for
 
@@ -18,7 +18,7 @@ Network <- R6Class( classname = 'Network', cloneable = FALSE,
 
     #' the C network R pointer object
     c_network_ptr = function() {
-      return( private$c_network()@ref )
+      return( self$c_network@ref )
     },
 
     #' check the C network still exists
 
@@ -106,7 +106,7 @@ Parameters <- R6Class( classname = 'Parameters', cloneable = FALSE,
 
     #' the C params R pointer object
     c_params_ptr = function() {
-      return( private$.c_params()@ref )
+      return( self$c_params@ref )
     },
 
     #' check the C params still exists
 
@@ -0,0 +1,79 @@
+#' R6Class Strain
+#'
+#' @description
+#' Strain object has information about each new strain
+#'
+#' @examples
+#' # Add new strain with increased transmissibility
+#' strain = model$add_new_strain( transmission_multiplier = 1.3 )
+#'
+Strain <- R6Class( classname = 'Strain', cloneable = FALSE,
+
+   private = list(
+     #' the strain ID
+     id        = NULL,
+
+     #' .c_strain External pointer, reference to \code{strain} C struct.
+     .c_strain = NULL,
+
+     #' the C strain R pointer object
+     c_strain_ptr = function() {
+       return( self$c_strain@ref )
+     },
+
+     #' check the C strain still exists
+     c_strain_valid = function() {
+       return( !is_null_xptr( private$.c_strain@ref ))
+     }
+   ),
+
+   active = list(
+     #' @field c_strain the C strain R pointer object (SWIG wrapped)
+     c_strain = function( val = NULL )
+     {
+       if( is.null( val ) )
+       {
+         if( private$c_strain_valid() )
+           return( private$.c_strain )
+         stop( "c_strain is no longer valid - create a new strain")
+       }
+       else
+         stop( "cannot set c_strain" )
+     }
+   ),
+
+   public = list(
+
+      #' @param model R6 Model object
+      #' @param stain_id The strain ID.
+      initialize = function( model, strain_id )
+      {
+         if( !is.R6(model) || !inherits( model, "Model"))
+            stop( "model must be a R6 class of type Model")
+
+         private$.c_strain <- get_strain_by_id( model$c_model, strain_id )
+         private$id         <- strain_id
+     },
+
+      #' @description Wrapper for C API \code{strain$idx()}.
+      #' @return the index of the strain
+      idx = function() {
+         return(strain_idx( self$c_strain ))
+      },
+
+      #' @description Wrapper for C API \code{strain$transmission_multiplier()}.
+      #' @return the transmission_multiplier of the strain
+      transmission_multiplier = function() {
+         return(strain_transmission_multiplier( self$c_strain ))
+      },
+
+      #' @description Wrapper for C API \code{strain$hospitalised_fraction()}.
+      #' @return the hospitalised fraction for the strain
+      hospitalised_fraction = function()
+      {
+         c_strain_ptr = private$c_strain_ptr();
+         return( .Call('R_strain_hospitalised_fraction',c_strain_ptr,PACKAGE='OpenABMCovid19') )
+      }
+   )
+
+)