sample states and within-host model

kingaa · Nov 15, 2024 · 5cdfc85 · 5cdfc85
1 parent f8f5dba
commit 5cdfc85
Show file tree

Hide file tree

Showing 49 changed files with 1,706 additions and 360 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: phylopomp
 Type: Package
 Title: Phylodynamic Inference for POMP Models
-Version: 0.14.8.0
-Date: 2024-11-15
+Version: 0.14.7.5
+Date: 2024-11-13
 Authors@R: c(person(given=c("Aaron","A."),family="King",role=c("aut","cre"),email="[email protected]",comment=c(ORCID="0000-0001-6159-3207")),
 	  person(given=c("Qianying"),family="Lin",role=c("aut"),comment=c(ORCID="0000-0001-8620-9910"))
 	  )
@@ -26,6 +26,7 @@ Collate:
     'diagram.R'
     'foreach.R'
     'gendat.R'
+    'get_states.R'
     'lbdp.R'
     'lbdp_exact.R'
     'lbdp_pomp.R'

diff --git a/R/geneal.R b/R/geneal.R
@@ -18,6 +18,7 @@ geneal <- function (object) {
     modelSI2R = .Call(P_genealSI2R,object),
     modelSIIR = .Call(P_genealSIIR,object),
     modelSIR = .Call(P_genealSIR,object),
+    modelTIMVA = .Call(P_genealTIMVA,object),
     modelTwoSpecies = .Call(P_genealTwoSpecies,object),
     model = structure(object,class=c("gpgen")),
     pStop("unrecognized model ",sQuote(attr(object,"model")))

diff --git a/R/get_states.R b/R/get_states.R
@@ -0,0 +1,29 @@
+##' Get population state history
+##'
+##' Retrieves the history of population states (compartment sizes over time)
+##' from a Markov genealogy process simulation
+##'
+##' @name get_states
+##' @param object A \sQuote{gpsim} object
+##' @return A \code{\link[tibble]{tibble}} containing the population state history
+##' @export
+get_states <- function(object) {
+  if (!inherits(object, "gpsim")) {
+    stop("object must be a gpsim object")
+  }
+
+  switch(
+    paste0("model", as.character(attr(object, "model"))),
+    modelLBDP = .Call(P_get_states_LBDP,object),
+    modelMoran = .Call(P_get_states_Moran,object),
+    modelS2I2R2 = .Call(P_get_states_S2I2R2,object),
+    modelSEIR = .Call(P_get_states_SEIR,object),
+    modelSI2R = .Call(P_get_states_SI2R,object),
+    modelSIIR = .Call(P_get_states_SIIR,object),
+    modelSIR = .Call(P_get_states_SIR,object),
+    modelTIMVA = .Call(P_get_states_TIMVA,object),
+    modelTwoSpecies = .Call(P_get_states_TwoSpecies,object),
+    stop("unrecognized model ", sQuote(attr(object, "model")))
+  ) |>
+    tibble::as_tibble()
+}
diff --git a/R/simulate.R b/R/simulate.R
@@ -35,6 +35,7 @@ simulate.default <- function (object, ...) {
       "- SI2R: Two-deme model of superspreading\n",
       "- SIIR: Two-strain SIR model\n",
       "- SIR: Classical susceptible-infected-recovered model\n",
+      "- TIMVA: Within-host viral dynamic model with ALT\n",
       "- TwoSpecies: Two-host infection model with waning, immigration, demography, and spillover. Hosts are culled upon sampling with a given probability.\n",
       "- SIRS: synonymous with SIR\n",
       "- SEIRS: synonymous with SEIR\n",
@@ -62,6 +63,7 @@ simulate.character <- function (object, time, ...) {
     modelSI2R = runSI2R(time=time,...),
     modelSIIR = runSIIR(time=time,...),
     modelSIR = runSIR(time=time,...),
+    modelTIMVA = runTIMVA(time=time,...),
     modelTwoSpecies = runTwoSpecies(time=time,...),
     modelSIRS = runSIR(time=time,...),
     modelSEIRS = runSEIR(time=time,...),
@@ -89,6 +91,7 @@ simulate.gpsim <- function (object, time, ...) {
     modelSI2R = continueSI2R(object,time=time,...),
     modelSIIR = continueSIIR(object,time=time,...),
     modelSIR = continueSIR(object,time=time,...),
+    modelTIMVA = continueTIMVA(object,time=time,...),
     modelTwoSpecies = continueTwoSpecies(object,time=time,...),
     model = pStop_("no model attribute detected."),
     pStop_("unrecognized model ",sQuote(model),".")

diff --git a/R/timva.R b/R/timva.R
@@ -0,0 +1,58 @@
+##' Within-host viral dynamic model with ALT
+##'
+##' Within-host viral dynamics of four types of cells and the released ALT.
+##'
+##' @name timva
+##' @family Genealogy processes
+##' @aliases TIMVA
+##' @param lambda target cell replenishment rate
+##' @param d target cell death rate
+##' @param Beta transmission rate
+##' @param delta infected cell death rate
+##' @param p free virus particle production rate per infected cell
+##' @param a immune cell activation rate per infected cell
+##' @param kappa infected cell killing rate
+##' @param c free virus particle clearance rate
+##' @param mu immune cell death rate
+##' @param omega proportion of ALT release per dead infected cell
+##' @param s ALT replenishment rate
+##' @param sigma ALT decay rate
+##' @param psi free virus particle sampling rate
+##' @param f probability of getting sampled
+##' @param dt time step for state recording
+##' @param T0 initial size of target cells
+##' @param I0 initial size of infected cells
+##' @param M0 initial size of immune cells
+##' @param V0 initial size of free virus particles
+##' @param A0 initial value of ALT
+##' @inheritParams sir
+##' @return \code{runTIMVA} and \code{continueTIMVA} return objects of class \sQuote{gpsim} with \sQuote{model} attribute \dQuote{TIMVA}.
+##'
+NULL
+
+##' @rdname timva
+##' @export
+runTIMVA <- function (
+  time, t0 = 0,
+  lambda = 10, d = 0, Beta = 1e-4, delta = 0.01, p = 10, a = 0.1, kappa = 0.005, c = 3, mu = 0.05, omega = 0.1, s = 3, sigma = 0.3, psi = 0.1, f = 0.1, dt = 0.1, T0 = 1000, I0 = 0, M0 = 0, V0 = 1, A0 = 0
+) {
+  params <- c(lambda=lambda,d=d,Beta=Beta,delta=delta,p=p,a=a,kappa=kappa,c=c,mu=mu,omega=omega,s=s,sigma=sigma,psi=psi,f=f,dt=dt)
+  ivps <- c(T0=T0,I0=I0,M0=M0,V0=V0,A0=A0)
+  x <- .Call(P_makeTIMVA,params,ivps,t0)
+  .Call(P_runTIMVA,x,time) |>
+    structure(model="TIMVA",class=c("gpsim","gpgen"))
+}
+
+##' @rdname timva
+##' @export
+continueTIMVA <- function (
+  object, time,
+  lambda = NA, d = NA, Beta = NA, delta = NA, p = NA, a = NA, kappa = NA, c = NA, mu = NA, omega = NA, s = NA, sigma = NA, psi = NA, f = NA, dt = NA
+) {
+  params <- c(
+    lambda=lambda,d=d,Beta=Beta,delta=delta,p=p,a=a,kappa=kappa,c=c,mu=mu,omega=omega,s=s,sigma=sigma,psi=psi,f=f,dt=dt
+  )
+  x <- .Call(P_reviveTIMVA,object,params)
+  .Call(P_runTIMVA,x,time) |>
+    structure(model="TIMVA",class=c("gpsim","gpgen"))
+}
diff --git a/R/yaml.R b/R/yaml.R
@@ -23,6 +23,7 @@ yaml <- function (object) {
     modelSI2R = .Call(P_yamlSI2R,object),
     modelSIIR = .Call(P_yamlSIIR,object),
     modelSIR = .Call(P_yamlSIR,object),
+    modelTIMVA = .Call(P_yamlTIMVA,object),
     modelTwoSpecies = .Call(P_yamlTwoSpecies,object),
     model = .Call(P_yaml,object),
     pStop("unrecognized model ",sQuote(attr(object,"model")))

diff --git a/src/decls.h b/src/decls.h
@@ -1,18 +1,6 @@
 /* src/init.c */
 extern void R_init_phylopomp(DllInfo *);
-/* src/lbdp_pomp.c */
-extern void lbdp_rinit(double *, const double *, double, const int *, const int *, const int *, const double *);
-extern void lbdp_gill(double *, const double *, const int *, const int *, const int *, const double *, double, double);
-extern void lbdp_dmeas(double *, const double *, const double *, const double *, int, const int *, const int *, const int *, const int *, const double *, double);
-/* src/seirs_pomp.c */
-extern void seirs_rinit(double *, const double *, double, const int *, const int *, const int *, const double *);
-extern void seirs_gill(double *, const double *, const int *, const int *, const int *, const double *, double, double);
-extern void seirs_dmeas(double *, const double *, const double *, const double *, int, const int *, const int *, const int *, const int *, const double *, double);
 /* src/sirs_pomp.c */
 extern void sirs_rinit(double *, const double *, double, const int *, const int *, const int *, const double *);
 extern void sirs_gill(double *, const double *, const int *, const int *, const int *, const double *, double, double);
 extern void sirs_dmeas(double *, const double *, const double *, const double *, int, const int *, const int *, const int *, const int *, const double *, double);
-/* src/twospecies_pomp.c */
-extern void twospecies_rinit(double *, const double *, double, const int *, const int *, const int *, const double *);
-extern void twospecies_gill(double *, const double *, const int *, const int *, const int *, const double *, double, double);
-extern void twospecies_dmeas(double *, const double *, const double *, const double *, int, const int *, const int *, const int *, const int *, const double *, double);
diff --git a/src/generics.h b/src/generics.h
@@ -63,7 +63,7 @@ SEXP lineage_count (const TYPE& G) {
   return G.lineage_count();
 }
 
-//! data-frame format 
+//! data-frame format
 template <class TYPE>
 SEXP gendat (const TYPE& G) {
   return G.gendat();
@@ -146,6 +146,11 @@ SEXP genealogy (SEXP State) {
     return yaml<TYPE>(State);                           \
   }                                                     \
 
+#define GET_STATES_FN(X,TYPE) SEXP get_states_ ## X (SEXP State) {   \
+    const TYPE &X = State;                                           \
+    return X.format_states();                                        \
+  }                                                                  \
+
 #define GENERICS(X,TYPE)                        \
   extern "C" {                                  \
                                                 \

diff --git a/src/init.c b/src/init.c
@@ -15,14 +15,15 @@ SEXP gendat (SEXP);
 // for each model, there must be
 // one DECLARATIONS line and one METHODS line.
 
-DECLARATIONS(LBDP);
-DECLARATIONS(Moran);
-DECLARATIONS(S2I2R2);
-DECLARATIONS(SEIR);
-DECLARATIONS(SI2R);
-DECLARATIONS(SIIR);
-DECLARATIONS(SIR);
-DECLARATIONS(TwoSpecies);
+DECLARATIONS(LBDP)
+DECLARATIONS(Moran)
+DECLARATIONS(S2I2R2)
+DECLARATIONS(SEIR)
+DECLARATIONS(SI2R)
+DECLARATIONS(SIIR)
+DECLARATIONS(SIR)
+DECLARATIONS(TIMVA)
+DECLARATIONS(TwoSpecies)
 
 static const R_CallMethodDef callMethods[] = {
   METHODS(LBDP),
@@ -32,6 +33,7 @@ static const R_CallMethodDef callMethods[] = {
   METHODS(SI2R),
   METHODS(SIIR),
   METHODS(SIR),
+  METHODS(TIMVA),
   METHODS(TwoSpecies),
   {"parse_newick", (DL_FUNC) &parse_newick, 3},
   {"curtail", (DL_FUNC) &curtail, 2},
@@ -49,8 +51,8 @@ void R_init_phylopomp (DllInfo *info) {
   // Register routines
   R_registerRoutines(info,NULL,callMethods,NULL,extMethods);
   R_useDynamicSymbols(info,TRUE);
-  // R_useDynamicSymbols(info,FALSE);
-  // R_forceSymbols(info,TRUE);
+  //  R_useDynamicSymbols(info,FALSE);
+  //  R_forceSymbols(info,TRUE);
   get_userdata = (get_userdata_t*) R_GetCCallable("pomp","get_userdata");
   get_userdata_double = (get_userdata_double_t*) R_GetCCallable("pomp","get_userdata_double");
   get_userdata_int = (get_userdata_int_t*) R_GetCCallable("pomp","get_userdata_int");

diff --git a/src/init.h b/src/init.h
@@ -4,18 +4,20 @@
 #include <R_ext/Rdynload.h>
 #include "internal.h"
 
-#define DECLARATIONS(X)                                 \
-  SEXP make ## X (SEXP Params, SEXP IVPs, SEXP T0);     \
-  SEXP revive ## X (SEXP State, SEXP Params);           \
-  SEXP run ## X (SEXP State, SEXP Times);               \
-  SEXP geneal ## X (SEXP State);                        \
-  SEXP yaml ## X (SEXP State)
+#define DECLARATIONS(X)							\
+  SEXP make ## X (SEXP Params, SEXP IVPs, SEXP T0);			\
+  SEXP revive ## X (SEXP State, SEXP Params);				\
+  SEXP run ## X (SEXP State, SEXP Times);				\
+  SEXP geneal ## X (SEXP State);					\
+  SEXP yaml ## X (SEXP State);						\
+  SEXP get_states_ ## X (SEXP State);        \
 
 #define METHODS(X)                              \
   {"make" #X, (DL_FUNC) &make ## X, 3},         \
   {"revive" #X, (DL_FUNC) &revive ## X, 2},     \
   {"run" #X, (DL_FUNC) &run ## X, 2},           \
-  {"yaml" #X, (DL_FUNC) &yaml ## X, 1},         \
-  {"geneal" #X, (DL_FUNC) &geneal ## X, 1}      \
+  {"yaml" #X, (DL_FUNC) &yaml ## X, 1},		\
+  {"geneal" #X, (DL_FUNC) &geneal ## X, 1},	\
+  {"get_states_" #X, (DL_FUNC) &get_states_ ## X, 1} \
 
 #endif
diff --git a/src/lbdp.cc b/src/lbdp.cc
@@ -16,6 +16,7 @@ typedef struct {
   double lambda;
   double mu;
   double psi;
+  double dt;
   int n0;
 } lbdp_parameters_t;
 
@@ -29,6 +30,7 @@ std::string lbdp_proc_t::yaml (std::string tab) const {
     + YAML_PARAM(lambda)
     + YAML_PARAM(mu)
     + YAML_PARAM(psi)
+    + YAML_PARAM(dt)
     + YAML_PARAM(n0);
   std::string s = tab + "state:\n"
     + YAML_STATE(n);
@@ -41,6 +43,7 @@ void lbdp_proc_t::update_params (double *p, int n) {
   PARAM_SET(lambda);
   PARAM_SET(mu);
   PARAM_SET(psi);
+  PARAM_SET(dt);
   if (m != n) err("wrong number of parameters!");
 }
 
@@ -65,25 +68,63 @@ double lbdp_proc_t::event_rates (double *rate, int n) const {
 template<>
 void lbdp_genealogy_t::rinit (void) {
   state.n = params.n0;
-  graft(deme,params.n0);
+graft(deme,params.n0);
 }
 
 template<>
 void lbdp_genealogy_t::jump (int event) {
   switch (event) {
   case 0:
-    state.n += 1; birth();
-    break;
-  case 1:
-    state.n -= 1; death();
-    break;
-  case 2:
-    sample();
-    break;
+      state.n += 1; birth();
+      break;
+    case 1:
+      state.n -= 1; death();
+      break;
+    case 2:
+      sample();
+      break;
   default:                      // #nocov
     assert(0);                  // #nocov
     break;                      // #nocov
   }
 }
 
+template<>
+size_t lbdp_proc_t::n_integer_elements() const {
+  return 1;  // Number of integer state variables
+}
+
+template<>
+size_t lbdp_proc_t::n_double_elements() const {
+  return 0;  // Number of double state variables
+}
+
+static const char* LBDP_int_names[] = {"n"};
+static const char* LBDP_dbl_names[] = {""};
+
+template<>
+const char** lbdp_proc_t::integer_names() const {
+  return LBDP_int_names;
+}
+
+template<>
+const char** lbdp_proc_t::double_names() const {
+  return LBDP_dbl_names;
+}
+
+template<>
+void lbdp_proc_t::get_state_elements(size_t i, double *time, int *intg, double *dbl) const {
+  *time = time_history[i];
+  const lbdp_state_t& s = state_history[i];
+    intg[0] = s.n;
+
+}
+
+extern "C" {
+  SEXP get_states_LBDP (SEXP State) {
+    lbdp_genealogy_t x(State);
+    return x.get_states();
+  }
+}
+
 GENERICS(LBDP,lbdp_genealogy_t)