two-species model partially coded for filtering

NAMESPACE

@@ -66,6 +66,7 @@ export(sir_pomp)
 export(sirs_pomp)
 export(stew)
 export(treeplot)
+export(twospecies_pomp)
 export(viewport)
 export(yaml)
 import(ggplot2)

R/moran.R

@@ -8,11 +8,11 @@
 ##' @aliases Moran
 ##' @include getinfo.R
 ##' @family Genealogy processes
-##' @param mu event rate
-##' @param psi sampling rate.
+##' @aliases Moran
+##' @param mu per capita event rate
+##' @param psi per capita sampling rate
 ##' @param n population size
-##' @param time final time
-##' @param t0 initial time
+##' @inheritParams sir
 ##' @return \code{runMoran} and \code{continueMoran} return objects of class \sQuote{gpsim} with \sQuote{model} attribute \dQuote{Moran}.
 ##' @references
 ##' \Moran1958
@@ -21,8 +21,8 @@ NULL
 ##' @rdname moran
 ##' @export
 runMoran <- function (
-  time,  t0 = 0, n = 100,
-  mu = 1, psi = 1
+  time, t0 = 0,
+  mu = 1, psi = 1, n = 100
 ) {
   params <- c(mu=mu,psi=psi)
   ivps <- c(n=n)

R/seir.R

@@ -3,16 +3,15 @@
 ##' The population is structured by infection progression.
 ##'
 ##' @name seir
-##' @aliases SEIR
 ##' @family Genealogy processes
-##' @include sir.R
+##' @aliases SEIR
 ##' @param Beta transmission rate
 ##' @param sigma progression rate
 ##' @param gamma recovery rate
 ##' @param psi per capita sampling rate
 ##' @param omega rate of waning of immunity
 ##' @param S0,E0,I0,R0 initial sizes of S, E, I, R compartments, respectively.
-##' @inheritParams runSIR
+##' @inheritParams sir
 ##' @return \code{runSEIR} and \code{continueSEIR} return objects of class \sQuote{gpsim} with \sQuote{model} attribute \dQuote{SEIR}.
 ##' @references
 ##' \King2024

R/twospecies.R

@@ -26,10 +26,29 @@ NULL
 ##' @export
 runTwoSpecies <- function (
   time, t0 = 0,
-  Beta11 = 4, Beta12 = 0, Beta21 = 0, Beta22 = 4, gamma1 = 1, gamma2 = 1, psi1 = 1, psi2 = 0, omega1 = 0, omega2 = 0, b1 = 0, b2 = 0, d1 = 0, d2 = 0, iota1 = 0, iota2 = 0, S1_0 = 100, S2_0 = 100, I1_0 = 0, I2_0 = 10, R1_0 = 0, R2_0 = 0
+  Beta11 = 0.5, Beta12 = 0.5, Beta21 = 0.06, Beta22 = 1.5,
+  gamma1 = 0.5, gamma2 = 0.5, psi1 = 2, psi2 = 0,
+  omega1 = 0, omega2 = 0,
+  b1 = 0.2, b2 = 0.5,
+  d1 = 0.2, d2 = 0.5,
+  iota1 = 0, iota2 = 0.001,
+  S1_0 = 1000, S2_0 = 300,
+  I1_0 = 0, I2_0 = 0,
+  R1_0 = 0, R2_0 = 700
 ) {
-  params <- c(Beta11=Beta11,Beta12=Beta12,Beta21=Beta21,Beta22=Beta22,gamma1=gamma1,gamma2=gamma2,psi1=psi1,psi2=psi2,omega1=omega1,omega2=omega2,b1=b1,b2=b2,d1=d1,d2=d2,iota1=iota1,iota2=iota2)
-  ivps <- c(S1_0=S1_0,S2_0=S2_0,I1_0=I1_0,I2_0=I2_0,R1_0=R1_0,R2_0=R2_0)
+  params <- c(
+    Beta11=Beta11,Beta12=Beta12,Beta21=Beta21,Beta22=Beta22,
+    gamma1=gamma1,gamma2=gamma2,
+    psi1=psi1,psi2=psi2,
+    omega1=omega1,omega2=omega2,
+    b1=b1,b2=b2,d1=d1,d2=d2,
+    iota1=iota1,iota2=iota2
+  )
+  ivps <- c(
+    S1_0=S1_0,S2_0=S2_0,
+    I1_0=I1_0,I2_0=I2_0,
+    R1_0=R1_0,R2_0=R2_0
+  )
   x <- .Call(P_makeTwoSpecies,params,ivps,t0)
   .Call(P_runTwoSpecies,x,time) |>
     structure(model="TwoSpecies",class=c("gpsim","gpgen"))
@@ -48,3 +67,58 @@ continueTwoSpecies <- function (
   .Call(P_runTwoSpecies,x,time) |>
     structure(model="TwoSpecies",class=c("gpsim","gpgen"))
 }
+
+##' @name twospecies_pomp
+##' @rdname twospecies
+##' @include seir.R
+##' @param x genealogy in \pkg{phylopomp} format.
+##' @return
+##' \code{twospecies_pomp} returns a \sQuote{pomp} object.
+##' @details
+##' \code{twospecies_pomp} constructs a \sQuote{pomp} object containing a given set of data and a TwoSpecies model.
+##' @importFrom pomp pomp onestep
+##' @export
+twospecies_pomp <- function (
+  x,
+  Beta11, Beta12, Beta21, Beta22,
+  gamma1, gamma2, psi1, psi2, omega1, omega2,
+  b1, b2, d1, d2, iota1, iota2,
+  S1_0, S2_0, I1_0, I2_0, R1_0, R2_0
+)
+{
+  x |> gendat() -> gi
+  ic <- as.integer(c(S1_0,S2_0,I1_0,I2_0,R1_0,R2_0))
+  names(ic) <- c("S1_0","S2_0","I1_0","I2_0","R1_0","R2_0")
+  if (any(ic < 0))
+    pStop(paste(sQuote(names(ic)),collapse=","),
+      " must be nonnegative integers.")
+  pomp(
+    data=NULL,
+    t0=gi$nodetime[1L],
+    times=gi$nodetime[-1L],
+    params=c(
+      Beta11=Beta11,Beta12=Beta12,Beta21=Beta21,Beta22=Beta22,
+      gamma1=gamma1,gamma2=gamma2,psi1=psi1,psi2=psi2,
+      omega1=omega1,omega2=omega2,
+      b1=b1,b2=b2,d1=d1,d2=d2,
+      iota1=iota1,iota2=iota2,
+      S1_0=S1_0,S2_0=S2_0,I1_0=I1_0,I2_0=I2_0,R1_0=R1_0,R2_0=R2_0
+    ),
+    userdata=gi,
+    nstatevars=12L+gi$nsample,
+    rinit="twospecies_rinit",
+    rprocess=onestep("twospecies_gill"),
+    dmeasure="twospecies_dmeas",
+    statenames=c(
+      "S1","S2","I1","I2","R1","R2","N1","N2",
+      "ll","node","ellE","ellI","color"
+    ),
+    paramnames=c(
+      "Beta11","Beta12","Beta21","Beta22",
+      "gamma1","gamma2","psi1","psi2","omega1","omega2",
+      "b1","b2","d1","d2","iota1","iota2",
+      "S1_0","S2_0","I1_0","I2_0","R1_0","R2_0"
+    ),
+    PACKAGE="phylopomp"
+  )
+}