rework initial condition specification for some models;

change argument names for `curtail()`

Author: kingaa

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: phylopomp
 Type: Package
 Title: Phylodynamic Inference for POMP Models
-Version: 0.17.0.0
-Date: 2026-03-30
+Version: 0.17.1.0
+Date: 2026-04-01
 Authors@R: c(person(given=c("Aaron","A."),family="King",role=c("aut","cre"),email="kingaa@umich.edu",comment=c(ORCID="0000-0001-6159-3207")),
 	  person(given=c("Qianying"),family="Lin",role=c("aut"),comment=c(ORCID="0000-0001-8620-9910"))
 	  )

NAMESPACE

@@ -126,6 +126,7 @@ importFrom(pomp,bake)
 importFrom(pomp,covariate_table)
 importFrom(pomp,freeze)
 importFrom(pomp,onestep)
+importFrom(pomp,parameter_trans)
 importFrom(pomp,pomp)
 importFrom(pomp,stew)
 importFrom(scales,alpha)

R/curtail.R

@@ -5,11 +5,11 @@
 ##' @name curtail
 ##' @include getinfo.R
 ##' @inheritParams getInfo
-##' @param troot new root time for curtailed genealogy
-##' @return A curtailed genealogy object.
+##' @param t0,time numeric scalars; determine the time interval for curtailed genealogy
+##' @return a genealogy object: the curtailed version of the input genealogy.
 ##' @example examples/curtail.R
 ##' @rdname curtail
 ##' @export
-curtail <- function (object, time = NA, troot = NA) {
-  .Call(P_curtail,geneal(object),time,troot)
+curtail <- function (object, time = NA, t0 = NA) {
+  .Call(P_curtail,geneal(object),time,t0)
 }

R/lbdp_pomp.R

@@ -8,9 +8,8 @@
 lbdp_pomp <- function (x, lambda, mu, psi, chi = 0, n0 = 1)
 {
   x |> gendat() -> gi
-  n0 <- round(n0)
-  if (n0 < 0)
-    pStop(sQuote("n0")," must be a nonnegative integer.")
+  if (round(n0) < 0)
+    pStop(sQuote("n0")," must be nonnegative.")
   pomp(
     data=NULL,
     t0=gi$nodetime[1L],

R/seir.R

@@ -10,7 +10,6 @@
 ##' @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 sir
 ##' @return \code{runSEIR} and \code{continueSEIR} return objects of class \sQuote{gpsim} with \sQuote{model} attribute \dQuote{SEIR}.
 ##' @references
@@ -22,13 +21,18 @@ NULL
 ##' @export
 runSEIR <- function (
   time, t0 = 0,
-  Beta = 4, sigma = 1, gamma = 1, psi = 1, omega = 0, S0 = 100, E0 = 5, I0 = 5, R0 = 0
+  Beta = 4, sigma = 1, gamma = 1, psi = 1, omega = 0,
+  S0 = 0.9, E0 = 0.05, I0 = 0.05, R0 = 0, pop = 100
 ) {
   params <- c(Beta=Beta,sigma=sigma,gamma=gamma,psi=psi,omega=omega)
-  ivps <- c(S0=S0,E0=E0,I0=I0,R0=R0)
+  ivps <- c(S0,E0,I0,R0)
+  ivps <- structure(
+    as.integer(round(pop*ivps/sum(ivps))),
+    names=c("S0","E0","I0","R0")
+  )
   if (any(ivps < 0))
     pStop(paste(sQuote(names(ivps)),collapse=","),
-      " must be nonnegative integers.")
+      " must be nonnegative.")
   x <- .Call(P_makeSEIR,params,ivps,t0)
   .Call(P_runSEIR,x,time) |>
     structure(model="SEIR",class=c("gpsim","gpgen"))

R/seir_pomp.R

@@ -2,6 +2,9 @@
 ##' @rdname seir
 ##' @include lbdp.R sir.R seir.R
 ##' @param x genealogy in \pkg{phylopomp} format.
+##' @param S0,E0,I0,R0 initial conditions;
+##' non-negative numbers that specify the relative occupancies of the compartments at the inital time.
+##' @param pop host population size
 ##' @return
 ##' \code{seirs_pomp} returns a \sQuote{pomp} object.
 ##' @details
@@ -11,22 +14,21 @@
 seirs_pomp <- function (
   x,
   Beta, sigma, gamma, psi, omega = 0,
-  S0, E0, I0, R0
+  S0, E0, I0, R0, pop
 )
 {
   x |> gendat() -> gi
-  ic <- as.integer(c(S0,E0,I0,R0))
-  names(ic) <- c("S0","E0","I0","R0")
-  if (any(ic < 0))
-    pStop(paste(sQuote(names(ic)),collapse=","),
-      " must be nonnegative integers.")
+  ivps <- structure(c(S0,E0,I0,R0),names=c("S0","E0","I0","R0"))
+  if (any(ivps < 0))
+    pStop(paste(sQuote(names(ivps)),collapse=","),
+      " must be nonnegative.")
   pomp(
     data=NULL,
     t0=gi$nodetime[1L],
     times=gi$nodetime[-1L],
     params=c(
       Beta=Beta,sigma=sigma,gamma=gamma,psi=psi,omega=omega,
-      ic,N=sum(ic)
+      ivps,pop=pop
     ),
     userdata=gi,
     nstatevars=8L + gi$nsample,
@@ -39,7 +41,7 @@ seirs_pomp <- function (
     ),
     paramnames=c(
       "Beta","sigma","gamma","psi","omega",
-      "S0","E0","I0","R0","N"
+      "S0","E0","I0","R0","pop"
     ),
     PACKAGE="phylopomp"
   )

R/sir.R

@@ -10,7 +10,6 @@
 ##' @param gamma recovery rate.
 ##' @param psi sampling rate.
 ##' @param omega immunity waning rate
-##' @param S0,I0,R0 initial sizes of susceptible, infected, and recovered populations, respectively.
 ##' @param time final time
 ##' @param t0 initial time
 ##' @return \code{runSIR} and \code{continueSIR} return objects of class \sQuote{gpsim} with \sQuote{model} attribute \dQuote{SIR}.
@@ -26,13 +25,17 @@ NULL
 runSIR <- function (
   time,  t0 = 0,
   Beta = 2, gamma = 1, psi = 1, omega = 0,
-  S0 = 100, I0 = 2, R0 = 0
+  S0 = 100, I0 = 2, R0 = 0, pop = 102
 ) {
   params <- c(Beta=Beta,gamma=gamma,psi=psi,omega=omega)
-  ivps <- c(S0=S0,I0=I0,R0=R0)
+  ivps <- c(S0,I0,R0)
+  ivps <- structure(
+    as.integer(round(pop*ivps/sum(ivps))),
+    names=c("S0","I0","R0")
+  )
   if (any(ivps < 0))
     pStop(paste(sQuote(names(ivps)),collapse=","),
-      " must be nonnegative integers.")
+      " must be nonnegative.")
   x <- .Call(P_makeSIR,params,ivps,t0)
   .Call(P_runSIR,x,time) |>
     structure(model="SIR",class=c("gpsim","gpgen"))

R/sir_pomp.R

@@ -2,34 +2,36 @@
 ##' @rdname sir
 ##' @include sir.R
 ##' @param x genealogy in \pkg{phylopomp} format (i.e., an object that inherits from \sQuote{gpgen}).
+##' @param S0,I0,R0 initial conditions;
+##' non-negative numbers that specify the relative occupancies of the compartments at the inital time.
+##' @param pop host population size
 ##' @details
 ##' \code{sir_pomp} constructs a \sQuote{pomp} object containing a given set of data and a SIR model.
 ##' @return
 ##' \code{sir_pomp} and \code{sirs_pomp} return \sQuote{pomp} objects.
 ##' @importFrom pomp pomp onestep
 ##' @export
 sir_pomp <- function (
-  x, Beta, gamma, psi, omega = 0, S0, I0, R0
+  x, Beta, gamma, psi, omega = 0, S0, I0, R0, pop
 ) {
   x |> gendat() -> gi
-  ic <- as.integer(c(S0,I0,R0))
-  names(ic) <- c("S0","I0","R0")
-  if (any(ic < 0))
-    pStop(paste(sQuote(names(ic)),collapse=","),
-      " must be nonnegative integers.")
+  ivps <- structure(c(S0,I0,R0),names=c("S0","I0","R0"))
+  if (any(ivps < 0))
+    pStop(paste(sQuote(names(ivps)),collapse=","),
+      " must be nonnegative.")
   pomp(
     data=NULL,
     t0=gi$nodetime[1L],
     times=gi$nodetime[-1L],
-    params=c(Beta=Beta,gamma=gamma,psi=psi,omega=omega,ic,N=sum(ic)),
+    params=c(Beta=Beta,gamma=gamma,psi=psi,omega=omega,ivps,pop=pop),
     userdata=gi,
     rinit="sirs_rinit",
     rprocess=onestep("sirs_gill"),
     dmeasure="sirs_dmeas",
     statenames=c("S","I","R","ll","ell","node"),
     paramnames=c(
       "Beta","gamma","psi","omega",
-      "S0","I0","R0","N"
+      names(ivps),"pop"
     ),
     PACKAGE="phylopomp"
   )

R/strains.R

@@ -9,9 +9,6 @@
 ##' @param Beta1,Beta2,Beta3 transmission rate for strains 1, 2, 3, respectively
 ##' @param gamma recovery rate
 ##' @param psi1,psi2,psi3 sampling rates
-##' @param S0 initial size of susceptible population
-##' @param I1_0,I2_0,I3_0 initial numbers of strain-specific infections
-##' @param R0 initial size of immune population
 ##' @inheritParams sir
 ##' @return \code{runStrains} and \code{continueStrains} return objects of class \sQuote{gpsim} with \sQuote{model} attribute \dQuote{Strains}.
 ##'
@@ -24,18 +21,22 @@ runStrains <- function (
   Beta1 = 5, Beta2 = 5, Beta3 = 5,
   gamma = 1,
   psi1 = 1, psi2 = 0, psi3 = 0,
-  S0 = 10000, I1_0 = 10, I2_0 = 10, I3_0 = 10, R0 = 0
+  S0 = 1, I1_0 = 0.01, I2_0 = 0.01, I3_0 = 0.01, R0 = 0,
+  pop=10000
 ) {
   params <- c(
     Beta1=Beta1,Beta2=Beta2,Beta3=Beta3,
     gamma=gamma,
     psi1=psi1,psi2=psi2,psi3=psi3
   )
-  ivps <- c(
-    S0=S0,
-    I1_0=I1_0,I2_0=I2_0,I3_0=I3_0,
-    R0=R0
+  ivps <- c(S0,I1_0,I2_0,I3_0,R0)
+  ivps <- structure(
+    as.integer(round(pop*ivps/sum(ivps))),
+    names=c("S0","I1_0","I2_0","I3_0","R0")
   )
+  if (any(ivps < 0))
+    pStop(paste(sQuote(names(ivps)),collapse=","),
+      " must be nonnegative.")
   x <- .Call(P_makeStrains,params,ivps,t0)
   .Call(P_runStrains,x,time) |>
     structure(model="Strains",class=c("gpsim","gpgen"))

R/strains_pomp.R

@@ -2,24 +2,28 @@
 ##' @rdname strains
 ##' @include strains.R
 ##' @param x genealogy in \pkg{phylopomp} format (i.e., an object that inherits from \sQuote{gpgen}).
+##' @param S0,I1_0,I2_0,I3_0,R0 initial conditions; non-negative numbers that specify the relative occupancies of the compartments at the inital time.
+##' @param pop host population size
 ##' @details
 ##' \code{strains_pomp} constructs a \sQuote{pomp} object containing a given set of data and the Strains model.
 ##' @return
 ##' \code{strains_pomp} returns a \sQuote{pomp} object.
-##' @importFrom pomp pomp onestep
+##' @importFrom pomp pomp onestep parameter_trans
 ##' @export
 strains_pomp <- function (
   x,
   Beta1, Beta2, Beta3, gamma,
-  psi1, psi2, psi3,
+  psi1, psi2, psi3, pop,
   S0, I1_0, I2_0, I3_0, R0
 ) {
   x |> gendat(obscure=FALSE) -> gi
-  ic <- as.integer(c(S0,I1_0,I2_0,I3_0,R0))
-  names(ic) <- c("S0","I1_0","I2_0","I3_0","R0")
-  if (any(ic < 0))
-    pStop(paste(sQuote(names(ic)),collapse=","),
-      " must be nonnegative integers.")
+  ivps <- structure(
+    c(S0,I1_0,I2_0,I3_0,R0),
+    names=c("S0","I1_0","I2_0","I3_0","R0")
+  )
+  if (any(ivps < 0))
+    pStop(paste(sQuote(names(ivps)),collapse=","),
+      " must be nonnegative.")
   pomp(
     data=NULL,
     t0=gi$nodetime[1L],
@@ -28,20 +32,27 @@ strains_pomp <- function (
       Beta1=Beta1,Beta2=Beta2,Beta3=Beta3,
       gamma=gamma,
       psi1=psi1,psi2=psi2,psi3=psi3,
-      ic,N=sum(ic)
+      ivps,pop=pop
     ),
     userdata=gi,
     rinit="strains_rinit",
     rprocess=onestep("strains_gill"),
     dmeasure="strains_dmeas",
+    partrans=parameter_trans(
+      log=c(
+        "Beta1","Beta2","Beta3","gamma",
+        "psi1","psi2","psi3"
+      ),
+      barycentric=c(names(ivps))
+    ),
     statenames=c(
       "S","I1","I2","I3","R",
       "ll","ell1","ell2","ell3","node"
     ),
     paramnames=c(
       "Beta1","Beta2","Beta3","gamma",
       "psi1","psi2","psi3",
-      names(ic),"N"
+      names(ivps),"pop"
     ),
     PACKAGE="phylopomp"
   )