From 3a33bed96de064672dd3c6d019d1ce3740528a0e Mon Sep 17 00:00:00 2001
From: WongDWai <dwong338@uwo.ca>
Date: Fri, 21 Aug 2026 10:42:41 -0400
Subject: [PATCH 1/8] working on sim.arg crash at high rho. Tried to add
 carrier.id to distinguish individuals; recombination parent lineages inherit
 source pathogen id, coalescent inherits a child's id. do.infection counts
 unique carrier.id for rhyper().

---
 R/Pathogen.R                 | 17 +++++++++++--
 R/simARG.R                   |  8 +++++++
 R/simInnerTree.R             | 46 ++++++++++++++++++++++++++++++------
 tests/testthat/test_simARG.R | 40 +++++++++++++++++++++++++++++++
 4 files changed, 102 insertions(+), 9 deletions(-)

diff --git a/R/Pathogen.R b/R/Pathogen.R
index b81f3b7..3562bf4 100644
--- a/R/Pathogen.R
+++ b/R/Pathogen.R
@@ -15,18 +15,28 @@
 #'        Pathogen objects or there will be a circular reference problem!)
 #' @param breakpoint:  integer, genomic position of recombination breakpoint;
 #'        NA for non-recombinant lineages
+#' @param carrier.id:  identifier for the physical host individual this
+#'        lineage's genome currently occupies. Distinct from the Pathogen's
+#'        own name/lineage identity: recombination splits one lineage into
+#'        two Pathogen objects that are still the same physical genome, so
+#'        both inherit the same carrier.id. New (non-recombinant) Pathogens
+#'        get a fresh carrier.id. Used to count physical individuals rather
+#'        than tracked lineages (see InnerTree registry / .do.infection),
+#'        so recombination-driven lineage growth doesn't inflate apparent
+#'        population size.
 #' @export
 Pathogen <- R6Class(
   "Pathogen",
   public = list(
     initialize = function(name=NA, start.time=NA, end.time=NA, parents=list(),
-                          children=list(), breakpoint=NA) {
+                          children=list(), breakpoint=NA, carrier.id=NA) {
       private$name <- name
       private$start.time <- start.time
       private$end.time <- end.time
       private$parents <- parents
       private$children <- children
       private$breakpoint <- breakpoint
+      private$carrier.id <- if (is.na(carrier.id)) name else carrier.id
     },
     
     print = function() {
@@ -41,6 +51,8 @@ Pathogen <- R6Class(
     # immutable attributes
     get.name = function() { private$name },
     get.end.time = function() { private$end.time },
+    get.carrier.id = function() { private$carrier.id },
+    set.carrier.id = function(id) { private$carrier.id <- id },
     
     # mutables
     get.start.time = function() { private$start.time },
@@ -78,6 +90,7 @@ Pathogen <- R6Class(
     end.time = NULL,
     parents = NULL,
     children = NULL,
-    breakpoint = NULL
+    breakpoint = NULL,
+    carrier.id = NULL
   )
 )
diff --git a/R/simARG.R b/R/simARG.R
index 251df9e..2c3c5c8 100644
--- a/R/simARG.R
+++ b/R/simARG.R
@@ -205,6 +205,14 @@ sim.arg <- function(outer, rho = 1e-4, seq.length = 9000L) {
   parent.left  <- inner$new.pathogen(time)
   parent.right <- inner$new.pathogen(time)
 
+  # both new lineages are still the SAME physical genome as `pathogen` --
+  # recombination splits which genomic interval each lineage tracks, not
+  # which individual carries it. Propagate carrier.id so downstream code
+  # (e.g. the superinfection bottleneck draw) counts physical individuals
+  # rather than tracked lineages.
+  parent.left$set.carrier.id(pathogen$get.carrier.id())
+  parent.right$set.carrier.id(pathogen$get.carrier.id())
+
   # record parent-child relationships (recombination has two parents)
   parent.left$add.child(pathogen)
   parent.right$add.child(pathogen)
diff --git a/R/simInnerTree.R b/R/simInnerTree.R
index e5c9cd7..bfd5b44 100644
--- a/R/simInnerTree.R
+++ b/R/simInnerTree.R
@@ -194,7 +194,15 @@ sim.inner.tree <- function(outer) {
   
   if (is.infected[[e$from.comp]]) {
     # superinfection
-    count <- recipient$count.pathogens()
+    # count PHYSICAL INDIVIDUALS (unique carrier.id), not tracked lineages --
+    # recombination can split one individual's genome into two Pathogen
+    # objects sharing the same carrier.id, which would otherwise inflate
+    # this count past the real population size and break the rhyper() draw
+    # below (see carrier.id design note on Pathogen)
+    paths.in.recipient <- recipient$get.pathogens()
+    carrier.ids <- sapply(paths.in.recipient, function(p) p$get.carrier.id())
+    unique.carriers <- unique(carrier.ids)
+    count <- length(unique.carriers)
     
     # determine bottleneck and population sizes
     expr <- inner$get.model()$get.bottleneck.size(e$to.comp)
@@ -202,13 +210,30 @@ sim.inner.tree <- function(outer) {
     expr <- inner$get.model()$get.pop.size(e$to.comp)
     p.size <- eval(parse(text=expr), envir=envir)
     
-    n.transfer <- rhyper(1, count, p.size-count, b.size)
+    # count is now the number of genuinely distinct physical individuals
+    # (carrier.id-deduped), but it can still legitimately exceed p.size if
+    # more distinct lineages have converged on this host than its nominal
+    # population size anticipates. rhyper() can't accept a negative second
+    # argument, so fall back to sampling from the nominal population size
+    # in that case rather than crashing.
+    count.for.draw <- min(count, p.size)
+    n.transfer <- rhyper(1, count.for.draw, p.size-count.for.draw, b.size)
     if (n.transfer > 0) {
-      for (i in 1:n.transfer) {
-        path <- recipient$remove.pathogen(1)
-        source$add.pathogen(path)
-        event$pathogen1 <- path$get.name()
-        inner$add.event(event)
+      # transfer whole individuals: every Pathogen sharing a selected
+      # carrier.id moves together, since they represent one physical genome
+      selected.carriers <- sample(unique.carriers, n.transfer)
+      for (cid in selected.carriers) {
+        moving.names <- sapply(paths.in.recipient[carrier.ids == cid],
+                                function(p) p$get.name())
+        for (pname in moving.names) {
+          current <- recipient$get.pathogens()
+          idx <- which(sapply(current, function(p) p$get.name()) == pname)
+          if (length(idx) != 1) next
+          path <- recipient$remove.pathogen(idx)
+          source$add.pathogen(path)
+          event$pathogen1 <- path$get.name()
+          inner$add.event(event)
+        }
       }
       if (!source.is.active) { active$add.host(source) }
       
@@ -321,6 +346,13 @@ sim.inner.tree <- function(outer) {
     p2$set.start.time(time)
     
     anc <- inner$new.pathogen(time)  # sets end.time
+    # anc is not a new physical individual -- it's the merged ancestor of
+    # p1 and p2, so it must inherit an existing carrier.id rather than
+    # default to a fresh one, or every coalescent event would spuriously
+    # inflate the carrier count. p1/p2 may originally have had different
+    # carrier.ids (two lineages that happen to share an ancestor without
+    # being the same individual); arbitrarily keep p1's.
+    anc$set.carrier.id(p1$get.carrier.id())
     host$add.pathogen(anc)
     
     # assign ancestral/descendant relations
diff --git a/tests/testthat/test_simARG.R b/tests/testthat/test_simARG.R
index c049d5b..fe31744 100644
--- a/tests/testthat/test_simARG.R
+++ b/tests/testthat/test_simARG.R
@@ -207,3 +207,43 @@ test_that("resolve.arg produces genuinely divergent topology (hand-built positiv
   expect_true(is.monophyletic(phy2, c("A","C")))
   expect_false(is.monophyletic(phy2, c("A","B")))
 })
+test_that("sim.arg does not crash at high rho (carrier.id regression)", {
+  settings <- read_yaml("test_Superinfection.yaml")
+  settings$Parameters$sigma <- 0.05
+  mod <- Model$new(settings)
+  set.seed(33)
+  dyn <- tryCatch(sim.dynamics(mod, max.attempts=10), error=function(e) NULL)
+  if (is.null(dyn)) skip("could not build dynamics for this seed")
+  outer <- tryCatch(
+    withCallingHandlers(sim.outer.tree(dyn), warning=function(w) invokeRestart("muffleWarning")),
+    error=function(e) NULL)
+  if (is.null(outer)) skip("could not build outer tree for this seed")
+
+  # rho=5 previously crashed with "missing value where TRUE/FALSE needed"
+  # (rhyper() given a negative population argument once recombination-
+  # created lineages inflated count past the host's nominal population
+  # size). carrier.id propagation in .do.recombination and .do.coalescent
+  # fixes this by counting physical individuals, not tracked lineages.
+  arg <- tryCatch(sim.arg(outer, rho=5, seq.length=9000), error=function(e) e)
+  expect_false(inherits(arg, "error"))
+  expect_true(!is.null(arg$breakpoints))
+})
+
+test_that("carrier.id: recombination children inherit parent's carrier.id", {
+  parent <- Pathogen$new(name="P1", end.time=1)
+  child1 <- Pathogen$new(name="P2", end.time=0.5)
+  child2 <- Pathogen$new(name="P3", end.time=0.5)
+  child1$set.carrier.id(parent$get.carrier.id())
+  child2$set.carrier.id(parent$get.carrier.id())
+  expect_equal(child1$get.carrier.id(), parent$get.carrier.id())
+  expect_equal(child2$get.carrier.id(), parent$get.carrier.id())
+})
+
+test_that("carrier.id: coalescent ancestor inherits a child's carrier.id, not a fresh one", {
+  p1 <- Pathogen$new(name="P1", end.time=1)
+  p2 <- Pathogen$new(name="P2", end.time=1)
+  anc <- Pathogen$new(name="P3", end.time=2)
+  anc$set.carrier.id(p1$get.carrier.id())
+  expect_equal(anc$get.carrier.id(), p1$get.carrier.id())
+  expect_false(anc$get.carrier.id() == anc$get.name())
+})

From f2db7879f24a016df9c6f18d9b1b1f8c3a723a82 Mon Sep 17 00:00:00 2001
From: WongDWai <dwong338@uwo.ca>
Date: Fri, 21 Aug 2026 11:03:57 -0400
Subject: [PATCH 2/8] Revert "working on sim.arg crash at high rho. Tried to
 add carrier.id to distinguish individuals; recombination parent lineages
 inherit source pathogen id, coalescent inherits a child's id. do.infection
 counts unique carrier.id for rhyper()."

This reverts commit 3a33bed96de064672dd3c6d019d1ce3740528a0e.
---
 R/Pathogen.R                 | 17 ++-----------
 R/simARG.R                   |  8 -------
 R/simInnerTree.R             | 46 ++++++------------------------------
 tests/testthat/test_simARG.R | 40 -------------------------------
 4 files changed, 9 insertions(+), 102 deletions(-)

diff --git a/R/Pathogen.R b/R/Pathogen.R
index 3562bf4..b81f3b7 100644
--- a/R/Pathogen.R
+++ b/R/Pathogen.R
@@ -15,28 +15,18 @@
 #'        Pathogen objects or there will be a circular reference problem!)
 #' @param breakpoint:  integer, genomic position of recombination breakpoint;
 #'        NA for non-recombinant lineages
-#' @param carrier.id:  identifier for the physical host individual this
-#'        lineage's genome currently occupies. Distinct from the Pathogen's
-#'        own name/lineage identity: recombination splits one lineage into
-#'        two Pathogen objects that are still the same physical genome, so
-#'        both inherit the same carrier.id. New (non-recombinant) Pathogens
-#'        get a fresh carrier.id. Used to count physical individuals rather
-#'        than tracked lineages (see InnerTree registry / .do.infection),
-#'        so recombination-driven lineage growth doesn't inflate apparent
-#'        population size.
 #' @export
 Pathogen <- R6Class(
   "Pathogen",
   public = list(
     initialize = function(name=NA, start.time=NA, end.time=NA, parents=list(),
-                          children=list(), breakpoint=NA, carrier.id=NA) {
+                          children=list(), breakpoint=NA) {
       private$name <- name
       private$start.time <- start.time
       private$end.time <- end.time
       private$parents <- parents
       private$children <- children
       private$breakpoint <- breakpoint
-      private$carrier.id <- if (is.na(carrier.id)) name else carrier.id
     },
     
     print = function() {
@@ -51,8 +41,6 @@ Pathogen <- R6Class(
     # immutable attributes
     get.name = function() { private$name },
     get.end.time = function() { private$end.time },
-    get.carrier.id = function() { private$carrier.id },
-    set.carrier.id = function(id) { private$carrier.id <- id },
     
     # mutables
     get.start.time = function() { private$start.time },
@@ -90,7 +78,6 @@ Pathogen <- R6Class(
     end.time = NULL,
     parents = NULL,
     children = NULL,
-    breakpoint = NULL,
-    carrier.id = NULL
+    breakpoint = NULL
   )
 )
diff --git a/R/simARG.R b/R/simARG.R
index 2c3c5c8..251df9e 100644
--- a/R/simARG.R
+++ b/R/simARG.R
@@ -205,14 +205,6 @@ sim.arg <- function(outer, rho = 1e-4, seq.length = 9000L) {
   parent.left  <- inner$new.pathogen(time)
   parent.right <- inner$new.pathogen(time)
 
-  # both new lineages are still the SAME physical genome as `pathogen` --
-  # recombination splits which genomic interval each lineage tracks, not
-  # which individual carries it. Propagate carrier.id so downstream code
-  # (e.g. the superinfection bottleneck draw) counts physical individuals
-  # rather than tracked lineages.
-  parent.left$set.carrier.id(pathogen$get.carrier.id())
-  parent.right$set.carrier.id(pathogen$get.carrier.id())
-
   # record parent-child relationships (recombination has two parents)
   parent.left$add.child(pathogen)
   parent.right$add.child(pathogen)
diff --git a/R/simInnerTree.R b/R/simInnerTree.R
index bfd5b44..e5c9cd7 100644
--- a/R/simInnerTree.R
+++ b/R/simInnerTree.R
@@ -194,15 +194,7 @@ sim.inner.tree <- function(outer) {
   
   if (is.infected[[e$from.comp]]) {
     # superinfection
-    # count PHYSICAL INDIVIDUALS (unique carrier.id), not tracked lineages --
-    # recombination can split one individual's genome into two Pathogen
-    # objects sharing the same carrier.id, which would otherwise inflate
-    # this count past the real population size and break the rhyper() draw
-    # below (see carrier.id design note on Pathogen)
-    paths.in.recipient <- recipient$get.pathogens()
-    carrier.ids <- sapply(paths.in.recipient, function(p) p$get.carrier.id())
-    unique.carriers <- unique(carrier.ids)
-    count <- length(unique.carriers)
+    count <- recipient$count.pathogens()
     
     # determine bottleneck and population sizes
     expr <- inner$get.model()$get.bottleneck.size(e$to.comp)
@@ -210,30 +202,13 @@ sim.inner.tree <- function(outer) {
     expr <- inner$get.model()$get.pop.size(e$to.comp)
     p.size <- eval(parse(text=expr), envir=envir)
     
-    # count is now the number of genuinely distinct physical individuals
-    # (carrier.id-deduped), but it can still legitimately exceed p.size if
-    # more distinct lineages have converged on this host than its nominal
-    # population size anticipates. rhyper() can't accept a negative second
-    # argument, so fall back to sampling from the nominal population size
-    # in that case rather than crashing.
-    count.for.draw <- min(count, p.size)
-    n.transfer <- rhyper(1, count.for.draw, p.size-count.for.draw, b.size)
+    n.transfer <- rhyper(1, count, p.size-count, b.size)
     if (n.transfer > 0) {
-      # transfer whole individuals: every Pathogen sharing a selected
-      # carrier.id moves together, since they represent one physical genome
-      selected.carriers <- sample(unique.carriers, n.transfer)
-      for (cid in selected.carriers) {
-        moving.names <- sapply(paths.in.recipient[carrier.ids == cid],
-                                function(p) p$get.name())
-        for (pname in moving.names) {
-          current <- recipient$get.pathogens()
-          idx <- which(sapply(current, function(p) p$get.name()) == pname)
-          if (length(idx) != 1) next
-          path <- recipient$remove.pathogen(idx)
-          source$add.pathogen(path)
-          event$pathogen1 <- path$get.name()
-          inner$add.event(event)
-        }
+      for (i in 1:n.transfer) {
+        path <- recipient$remove.pathogen(1)
+        source$add.pathogen(path)
+        event$pathogen1 <- path$get.name()
+        inner$add.event(event)
       }
       if (!source.is.active) { active$add.host(source) }
       
@@ -346,13 +321,6 @@ sim.inner.tree <- function(outer) {
     p2$set.start.time(time)
     
     anc <- inner$new.pathogen(time)  # sets end.time
-    # anc is not a new physical individual -- it's the merged ancestor of
-    # p1 and p2, so it must inherit an existing carrier.id rather than
-    # default to a fresh one, or every coalescent event would spuriously
-    # inflate the carrier count. p1/p2 may originally have had different
-    # carrier.ids (two lineages that happen to share an ancestor without
-    # being the same individual); arbitrarily keep p1's.
-    anc$set.carrier.id(p1$get.carrier.id())
     host$add.pathogen(anc)
     
     # assign ancestral/descendant relations
diff --git a/tests/testthat/test_simARG.R b/tests/testthat/test_simARG.R
index fe31744..c049d5b 100644
--- a/tests/testthat/test_simARG.R
+++ b/tests/testthat/test_simARG.R
@@ -207,43 +207,3 @@ test_that("resolve.arg produces genuinely divergent topology (hand-built positiv
   expect_true(is.monophyletic(phy2, c("A","C")))
   expect_false(is.monophyletic(phy2, c("A","B")))
 })
-test_that("sim.arg does not crash at high rho (carrier.id regression)", {
-  settings <- read_yaml("test_Superinfection.yaml")
-  settings$Parameters$sigma <- 0.05
-  mod <- Model$new(settings)
-  set.seed(33)
-  dyn <- tryCatch(sim.dynamics(mod, max.attempts=10), error=function(e) NULL)
-  if (is.null(dyn)) skip("could not build dynamics for this seed")
-  outer <- tryCatch(
-    withCallingHandlers(sim.outer.tree(dyn), warning=function(w) invokeRestart("muffleWarning")),
-    error=function(e) NULL)
-  if (is.null(outer)) skip("could not build outer tree for this seed")
-
-  # rho=5 previously crashed with "missing value where TRUE/FALSE needed"
-  # (rhyper() given a negative population argument once recombination-
-  # created lineages inflated count past the host's nominal population
-  # size). carrier.id propagation in .do.recombination and .do.coalescent
-  # fixes this by counting physical individuals, not tracked lineages.
-  arg <- tryCatch(sim.arg(outer, rho=5, seq.length=9000), error=function(e) e)
-  expect_false(inherits(arg, "error"))
-  expect_true(!is.null(arg$breakpoints))
-})
-
-test_that("carrier.id: recombination children inherit parent's carrier.id", {
-  parent <- Pathogen$new(name="P1", end.time=1)
-  child1 <- Pathogen$new(name="P2", end.time=0.5)
-  child2 <- Pathogen$new(name="P3", end.time=0.5)
-  child1$set.carrier.id(parent$get.carrier.id())
-  child2$set.carrier.id(parent$get.carrier.id())
-  expect_equal(child1$get.carrier.id(), parent$get.carrier.id())
-  expect_equal(child2$get.carrier.id(), parent$get.carrier.id())
-})
-
-test_that("carrier.id: coalescent ancestor inherits a child's carrier.id, not a fresh one", {
-  p1 <- Pathogen$new(name="P1", end.time=1)
-  p2 <- Pathogen$new(name="P2", end.time=1)
-  anc <- Pathogen$new(name="P3", end.time=2)
-  anc$set.carrier.id(p1$get.carrier.id())
-  expect_equal(anc$get.carrier.id(), p1$get.carrier.id())
-  expect_false(anc$get.carrier.id() == anc$get.name())
-})

From a641953b2adaf6654d702b7b0260be9c8fac805e Mon Sep 17 00:00:00 2001
From: WongDWai <dwong338@uwo.ca>
Date: Fri, 21 Aug 2026 11:07:33 -0400
Subject: [PATCH 3/8] Replace opaque rhyper() crash with explicit diagnostic
 error. Reverted earlier carrier.id because I don't think it was modelled in a
 biologically accurate way.

---
 R/simInnerTree.R | 23 +++++++++++++++++++++++
 1 file changed, 23 insertions(+)

diff --git a/R/simInnerTree.R b/R/simInnerTree.R
index e5c9cd7..8e00e10 100644
--- a/R/simInnerTree.R
+++ b/R/simInnerTree.R
@@ -202,6 +202,29 @@ sim.inner.tree <- function(outer) {
     expr <- inner$get.model()$get.pop.size(e$to.comp)
     p.size <- eval(parse(text=expr), envir=envir)
     
+    if (count > p.size) {
+      # count (n.active.lineages) is tracked ancestral segment/lineage
+      # objects, not necessarily distinct physical genomes -- once
+      # recombination is active, lineage count can legitimately exceed
+      # the population's census size, and the correct relationship
+      # between the two is not yet resolved (see issue tracker: rhyper()
+      # here implicitly assumes count occupies count distinct slots
+      # among p.size exchangeable individuals, which breaks once one
+      # genome can carry several ancestral segments). Fail loudly with
+      # a clear diagnostic rather than silently capping or crashing on
+      # an opaque NA, until the bottleneck/occupancy model is revisited.
+      stop(sprintf(
+        "Tracked lineage count (%d) exceeds pathogen population size (%d) ",
+        count, p.size),
+        "in host ", recipient$get.name(), " at time ", e$time, ". ",
+        "This means more ancestral lineages/segments are being tracked ",
+        "than the model's nominal population size anticipates (likely ",
+        "from recombination). The bottleneck sampling here assumes ",
+        "count occupies count distinct slots among p.size individuals, ",
+        "which is not valid once lineage count and physical individual ",
+        "count can diverge -- needs a proper occupancy/carrier model, ",
+        "not a silent cap.")
+    }
     n.transfer <- rhyper(1, count, p.size-count, b.size)
     if (n.transfer > 0) {
       for (i in 1:n.transfer) {

From da3dd6ae8d953c20b817aa906c0a208957397340 Mon Sep 17 00:00:00 2001
From: WongDWai <dwong338@uwo.ca>
Date: Fri, 28 Aug 2026 07:27:04 -0400
Subject: [PATCH 4/8] Fix recombination lineage growth. sample a recombination
 parent from a fixed pool of p.size lineages (a proportion of which are
 currently active), only activating a previously-inactive one when chosen.
 Total pool size never changes.

---
 R/Host.R                     | 33 ++++++++++++++-
 R/Model.R                    | 29 ++++++++++++-
 R/Pathogen.R                 |  5 ++-
 R/simARG.R                   | 75 ++++++++++++++++++++++++++++++----
 R/simInnerTree.R             | 21 ++++++++++
 tests/testthat/test_simARG.R | 79 ++++++++++++++++++++++++++++++++++++
 6 files changed, 230 insertions(+), 12 deletions(-)

diff --git a/R/Host.R b/R/Host.R
index b2bda7e..9057645 100644
--- a/R/Host.R
+++ b/R/Host.R
@@ -82,6 +82,35 @@ Host <- R6Class(
       private$pathogens[[length(private$pathogens)+1]] <- new.pathogen
     },
     
+    # --- fixed-size lineage pool (for recombination parent sampling) ---
+    # Fixes exponential lineage growth under recombination: instead of
+    # always creating a brand-new ancestral lineage, recombination should
+    # sample a parent from a FIXED pool of n (= population size) possible
+    # lineages, only "activating" a previously-inactive one when chosen.
+    # Total pool size never changes; only how many slots are active does.
+    init.pool = function(n) {
+      if (is.null(private$pool.size)) {
+        private$pool.size <- n
+        private$lineage.pool <- vector("list", n)
+      }
+    },
+    get.pool.size = function() { private$pool.size },
+    is.pool.initialized = function() { !is.null(private$pool.size) },
+    is.slot.active = function(slot.id) { !is.null(private$lineage.pool[[slot.id]]) },
+    get.slot.occupant = function(slot.id) { private$lineage.pool[[slot.id]] },
+    activate.slot = function(slot.id, pathogen) {
+      private$lineage.pool[[slot.id]] <- pathogen
+    },
+    deactivate.slot = function(slot.id) {
+      private$lineage.pool[slot.id] <- list(NULL)  # preserves pool length
+    },
+    get.active.slots = function() {
+      which(!sapply(private$lineage.pool, is.null))
+    },
+    get.inactive.slots = function() {
+      which(sapply(private$lineage.pool, is.null))
+    },
+    
     remove.pathogen = function(idx) {
       path <- private$pathogens[[idx]]
       private$pathogens[[idx]] <- NULL
@@ -136,7 +165,9 @@ Host <- R6Class(
     sampling.time = NULL,
     sampling.comp = NULL,
     unsampled = NULL,
-    pathogens = NULL
+    pathogens = NULL,
+    lineage.pool = NULL,
+    pool.size = NULL
   )
 )
 
diff --git a/R/Model.R b/R/Model.R
index 7ef2b1d..44dc439 100644
--- a/R/Model.R
+++ b/R/Model.R
@@ -316,7 +316,21 @@ Model <- R6Class(
             stop(src, "`coalescent.rate` should be a number or R expression,",
                  "not an associative array.")
           }
-          private$check.expression(params$coalescent.rate, env)
+          # coalescent.rate is allowed to reference `k` (current active
+          # lineage count) and `p.size` (compartment population size) at
+          # runtime -- these are simulation-state variables that don't
+          # exist at model-construction time, so validate against a copy
+          # of `env` with placeholder values instead of the real `env`.
+          # k=2 matches the minimum lineage count coalescent.rate is ever
+          # evaluated for; p.size uses the compartment's real declared
+          # pop.size when available, falling back to a placeholder.
+          coal.env <- list2env(as.list(env), parent = parent.env(env))
+          assign("k", 2, envir = coal.env)
+          p.size.placeholder <- tryCatch(
+            eval(parse(text = params$pop.size), envir = env),
+            error = function(e) 1)
+          assign("p.size", p.size.placeholder, envir = coal.env)
+          private$check.expression(params$coalescent.rate, coal.env)
           private$coalescent.rates[[src]] <- params$coalescent.rate
         }
         
@@ -328,6 +342,19 @@ Model <- R6Class(
           }
           private$check.expression(params$pop.size, env)
           private$pop.sizes[[src]] <- params$pop.size
+        } else if (!is.null(params$coalescent.rate) &&
+                   params$coalescent.rate != "Inf") {
+          # coalescence is enabled for this compartment (finite rate) but
+          # pop.size was never explicitly set, so it's silently using the
+          # unconfigured default (100) baked into this package -- this
+          # default has no connection to the compartment's own `size:`
+          # field or any other yaml value, and coalescent.rate expressions
+          # that reference population size (or were chosen assuming a
+          # particular population size) may be silently mismatched with it.
+          warning(src, ": coalescent.rate is set but pop.size is not -- ",
+                  "using unconfigured default pop.size=100. If ",
+                  "coalescent.rate was chosen with a particular population ",
+                  "size in mind, set pop.size explicitly in the yaml.")
         }
       }  # end loop
       
diff --git a/R/Pathogen.R b/R/Pathogen.R
index b81f3b7..67ebcde 100644
--- a/R/Pathogen.R
+++ b/R/Pathogen.R
@@ -41,6 +41,8 @@ Pathogen <- R6Class(
     # immutable attributes
     get.name = function() { private$name },
     get.end.time = function() { private$end.time },
+    get.slot.id = function() { private$slot.id },
+    set.slot.id = function(id) { private$slot.id <- id },
     
     # mutables
     get.start.time = function() { private$start.time },
@@ -78,6 +80,7 @@ Pathogen <- R6Class(
     end.time = NULL,
     parents = NULL,
     children = NULL,
-    breakpoint = NULL
+    breakpoint = NULL,
+    slot.id = NA
   )
 )
diff --git a/R/simARG.R b/R/simARG.R
index 251df9e..406a8b3 100644
--- a/R/simARG.R
+++ b/R/simARG.R
@@ -63,8 +63,11 @@ sim.arg <- function(outer, rho = 1e-4, seq.length = 9000L) {
         if (ev$type == "coalescent") {
           .do.coalescent(ev$host, inner, event.time, envir = env)
         } else {
+          host.obj <- active$get.host.by.name(ev$host)
+          p.size.expr <- mod$get.pop.size(host.obj$get.compartment())
+          p.size.val <- eval(parse(text = p.size.expr), envir = env)
           bp <- .do.recombination(ev$host, ev$pathogen, inner, event.time,
-                                  seq.length = seq.length)
+                                  seq.length = seq.length, p.size = p.size.val)
           breakpoints[[bp$child]] <- bp$position
         }
       }
@@ -133,6 +136,19 @@ sim.arg <- function(outer, rho = 1e-4, seq.length = 9000L) {
 
     # coalescence rate for this host (requires 2+ lineages)
     if (k >= 2) {
+      # expose current lineage count (k) and the compartment's nominal
+      # population size (p.size) to the coalescent.rate expression, so a
+      # yaml can write a population-size-scaled rate (e.g. "2/p.size")
+      # instead of only a flat constant. Needed so coalescence can
+      # naturally pull lineage count back toward p.size as it grows,
+      # rather than assuming it independently -- see rhyper() bottleneck
+      # crash / carrier.id discussion.
+      assign("k", k, envir = envir)
+      size.expr <- mod$get.pop.size(comp)
+      p.size.local <- tryCatch(eval(parse(text = size.expr), envir = envir),
+                                error = function(e) NA)
+      assign("p.size", p.size.local, envir = envir)
+
       expr <- mod$get.coalescent.rate(comp)
       rate <- eval(parse(text = expr), envir = envir)
       if (rate > 0) {
@@ -190,20 +206,58 @@ sim.arg <- function(outer, rho = 1e-4, seq.length = 9000L) {
 #' @keywords internal
 #' @noRd
 .do.recombination <- function(host.name, pathogen, inner, time,
-                              seq.length = 9000L) {
+                              seq.length = 9000L, p.size = NULL) {
   active <- inner$get.active()
   host   <- active$get.host.by.name(host.name)
 
+  # lazily initialize this host's FIXED lineage pool (Art's fix: sample
+  # recombination parents from a fixed pool of p.size lineages, only
+  # activating a previously-inactive one when chosen, instead of always
+  # creating a brand-new lineage de novo. Total pool size never changes.)
+  if (!host$is.pool.initialized()) {
+    if (is.null(p.size)) {
+      stop(".do.recombination: p.size must be provided to initialize ",
+           "host's lineage pool on first use")
+    }
+    host$init.pool(p.size)
+  }
+
+  # ensure the recombining pathogen already occupies a pool slot -- if
+  # this is the first pool-tracked event involving it, assign one now
+  if (is.na(pathogen$get.slot.id())) {
+    free <- host$get.inactive.slots()
+    slot <- if (length(free) > 0) free[1] else 1
+    pathogen$set.slot.id(slot)
+    host$activate.slot(slot, pathogen)
+  }
+  own.slot <- pathogen$get.slot.id()
+
   # sample breakpoint uniformly across genome
   breakpoint <- sample.int(seq.length - 1L, 1L)
   pathogen$set.breakpoint(breakpoint)
-
-  # end the current lineage at this recombination event
   pathogen$set.start.time(time)
 
-  # create two parental lineages — left and right of breakpoint
-  parent.left  <- inner$new.pathogen(time)
-  parent.right <- inner$new.pathogen(time)
+  # LEFT parent: continues in the SAME slot as the child
+  parent.left <- inner$new.pathogen(time)
+  parent.left$set.slot.id(own.slot)
+  host$activate.slot(own.slot, parent.left)
+
+  # RIGHT parent: sample a slot from the fixed pool (excluding own slot)
+  pool.size <- host$get.pool.size()
+  other.slots <- setdiff(seq_len(pool.size), own.slot)
+
+  if (length(other.slots) == 0) {
+    parent.right <- parent.left
+  } else {
+    sampled.slot <- if (length(other.slots) == 1) other.slots else sample(other.slots, 1)
+    if (host$is.slot.active(sampled.slot)) {
+      parent.right <- host$get.slot.occupant(sampled.slot)
+    } else {
+      parent.right <- inner$new.pathogen(time)
+      parent.right$set.slot.id(sampled.slot)
+      host$activate.slot(sampled.slot, parent.right)
+    }
+  }
 
   # record parent-child relationships (recombination has two parents)
   parent.left$add.child(pathogen)
@@ -216,9 +270,12 @@ sim.arg <- function(outer, rho = 1e-4, seq.length = 9000L) {
   idx <- which(sapply(paths, function(p) p$get.name()) == pathogen$get.name())
   if (length(idx) == 1) host$remove.pathogen(idx)
   host$add.pathogen(parent.left)
-  host$add.pathogen(parent.right)
+  already.present <- any(sapply(host$get.pathogens(), function(p) {
+    p$get.name() == parent.right$get.name()
+  }))
+  if (!already.present) host$add.pathogen(parent.right)
 
-  # log the recombination event (breakpoint not stored in log — fixed schema)
+  # log the recombination event (breakpoint not stored in log -- fixed schema)
   event <- list(
     time = time, event = "recombination",
     from.comp = host$get.compartment(), to.comp = NA,
diff --git a/R/simInnerTree.R b/R/simInnerTree.R
index 8e00e10..b8a574d 100644
--- a/R/simInnerTree.R
+++ b/R/simInnerTree.R
@@ -344,6 +344,27 @@ sim.inner.tree <- function(outer) {
     p2$set.start.time(time)
     
     anc <- inner$new.pathogen(time)  # sets end.time
+
+    # reconcile the lineage pool: the ancestor represents the same
+    # physical individual as whichever of p1/p2 already occupied a pool
+    # slot (from a prior recombination event). If both occupied slots,
+    # keep one for the ancestor and free the other -- two active
+    # lineages coalescing means one fewer active individual going
+    # forward. If neither occupied a slot, this coalescence never
+    # touched the pool, so the ancestor stays unassigned too.
+    p1.slot <- p1$get.slot.id()
+    p2.slot <- p2$get.slot.id()
+    if (!is.na(p1.slot)) {
+      anc$set.slot.id(p1.slot)
+      host$activate.slot(p1.slot, anc)
+      if (!is.na(p2.slot) && p2.slot != p1.slot) {
+        host$deactivate.slot(p2.slot)
+      }
+    } else if (!is.na(p2.slot)) {
+      anc$set.slot.id(p2.slot)
+      host$activate.slot(p2.slot, anc)
+    }
+
     host$add.pathogen(anc)
     
     # assign ancestral/descendant relations
diff --git a/tests/testthat/test_simARG.R b/tests/testthat/test_simARG.R
index c049d5b..d4d9ab4 100644
--- a/tests/testthat/test_simARG.R
+++ b/tests/testthat/test_simARG.R
@@ -207,3 +207,82 @@ test_that("resolve.arg produces genuinely divergent topology (hand-built positiv
   expect_true(is.monophyletic(phy2, c("A","C")))
   expect_false(is.monophyletic(phy2, c("A","B")))
 })
+test_that("Host lineage pool: basic activate/deactivate/query", {
+  h <- Host$new(name="H1", compartment="I")
+  expect_false(h$is.pool.initialized())
+  h$init.pool(5)
+  expect_true(h$is.pool.initialized())
+  expect_equal(h$get.pool.size(), 5)
+  expect_equal(length(h$get.inactive.slots()), 5)
+  expect_equal(length(h$get.active.slots()), 0)
+
+  p <- Pathogen$new(name="P1", end.time=1)
+  h$activate.slot(2, p)
+  expect_true(h$is.slot.active(2))
+  expect_equal(h$get.slot.occupant(2)$get.name(), "P1")
+  expect_equal(length(h$get.active.slots()), 1)
+  expect_equal(length(h$get.inactive.slots()), 4)
+
+  h$deactivate.slot(2)
+  expect_false(h$is.slot.active(2))
+  expect_equal(length(h$get.active.slots()), 0)
+  # pool size must not shrink after deactivation
+  expect_equal(h$get.pool.size(), 5)
+})
+
+test_that("Host lineage pool: init.pool is idempotent", {
+  h <- Host$new(name="H1", compartment="I")
+  h$init.pool(10)
+  h$init.pool(999)  # should be ignored, pool already initialized
+  expect_equal(h$get.pool.size(), 10)
+})
+
+test_that("sim.arg does not exceed pop.size at high rho (fixed lineage pool)", {
+  settings <- read_yaml("test_Superinfection.yaml")
+  settings$Parameters$sigma <- 0.05
+  mod <- Model$new(settings)
+  set.seed(33)
+  dyn <- tryCatch(sim.dynamics(mod, max.attempts=10), error=function(e) NULL)
+  if (is.null(dyn)) skip("could not build dynamics for this seed")
+  outer <- tryCatch(
+    withCallingHandlers(sim.outer.tree(dyn), warning=function(w) invokeRestart("muffleWarning")),
+    error=function(e) NULL)
+  if (is.null(outer)) skip("could not build outer tree for this seed")
+
+  # rho=5,7,10 previously exceeded pop.size and crashed/errored before
+  # the fixed lineage pool (Art's design: sample recombination parents
+  # from a fixed pool of p.size lineages instead of always creating a
+  # new lineage de novo). Now they should all succeed.
+  for (rho in c(5, 7, 10)) {
+    arg <- tryCatch(sim.arg(outer, rho=rho, seq.length=9000), error=function(e) e)
+    expect_false(inherits(arg, "error"),
+                 info=paste("rho =", rho, "should not error with fixed lineage pool"))
+  }
+})
+
+test_that("resolve.arg produces valid trees on top of the fixed lineage pool", {
+  settings <- read_yaml("test_Superinfection.yaml")
+  settings$Parameters$sigma <- 0.05
+  mod <- Model$new(settings)
+  set.seed(33)
+  dyn <- tryCatch(sim.dynamics(mod, max.attempts=10), error=function(e) NULL)
+  if (is.null(dyn)) skip("could not build dynamics for this seed")
+  outer <- tryCatch(
+    withCallingHandlers(sim.outer.tree(dyn), warning=function(w) invokeRestart("muffleWarning")),
+    error=function(e) NULL)
+  if (is.null(outer)) skip("could not build outer tree for this seed")
+
+  n.sampled <- outer$get.sampled()$count.type()
+  arg <- sim.arg(outer, rho=2, seq.length=9000)
+  res <- resolve.arg(arg, seq.length=9000)
+
+  valid <- sapply(res$local.trees, function(lt) {
+    phy <- lt$phylo
+    inherits(phy, "phylo") &&
+      length(phy$tip.label) == n.sampled &&
+      sum(duplicated(phy$tip.label)) == 0 &&
+      !any(is.na(phy$edge.length)) &&
+      !any(phy$edge.length < 0, na.rm=TRUE)
+  })
+  expect_true(all(valid))
+})

From 95e9036d89758926f04dbe31f8d6f640547a2f41 Mon Sep 17 00:00:00 2001
From: WongDWai <dwong338@uwo.ca>
Date: Fri, 28 Aug 2026 09:28:05 -0400
Subject: [PATCH 5/8] Remove orphaned k/p.size exposure on coalescent.rate.

---
 R/Model.R  | 16 +---------------
 R/simARG.R | 12 ------------
 2 files changed, 1 insertion(+), 27 deletions(-)

diff --git a/R/Model.R b/R/Model.R
index 44dc439..a73646b 100644
--- a/R/Model.R
+++ b/R/Model.R
@@ -316,21 +316,7 @@ Model <- R6Class(
             stop(src, "`coalescent.rate` should be a number or R expression,",
                  "not an associative array.")
           }
-          # coalescent.rate is allowed to reference `k` (current active
-          # lineage count) and `p.size` (compartment population size) at
-          # runtime -- these are simulation-state variables that don't
-          # exist at model-construction time, so validate against a copy
-          # of `env` with placeholder values instead of the real `env`.
-          # k=2 matches the minimum lineage count coalescent.rate is ever
-          # evaluated for; p.size uses the compartment's real declared
-          # pop.size when available, falling back to a placeholder.
-          coal.env <- list2env(as.list(env), parent = parent.env(env))
-          assign("k", 2, envir = coal.env)
-          p.size.placeholder <- tryCatch(
-            eval(parse(text = params$pop.size), envir = env),
-            error = function(e) 1)
-          assign("p.size", p.size.placeholder, envir = coal.env)
-          private$check.expression(params$coalescent.rate, coal.env)
+          private$check.expression(params$coalescent.rate, env)
           private$coalescent.rates[[src]] <- params$coalescent.rate
         }
         
diff --git a/R/simARG.R b/R/simARG.R
index 406a8b3..02fa2df 100644
--- a/R/simARG.R
+++ b/R/simARG.R
@@ -136,18 +136,6 @@ sim.arg <- function(outer, rho = 1e-4, seq.length = 9000L) {
 
     # coalescence rate for this host (requires 2+ lineages)
     if (k >= 2) {
-      # expose current lineage count (k) and the compartment's nominal
-      # population size (p.size) to the coalescent.rate expression, so a
-      # yaml can write a population-size-scaled rate (e.g. "2/p.size")
-      # instead of only a flat constant. Needed so coalescence can
-      # naturally pull lineage count back toward p.size as it grows,
-      # rather than assuming it independently -- see rhyper() bottleneck
-      # crash / carrier.id discussion.
-      assign("k", k, envir = envir)
-      size.expr <- mod$get.pop.size(comp)
-      p.size.local <- tryCatch(eval(parse(text = size.expr), envir = envir),
-                                error = function(e) NA)
-      assign("p.size", p.size.local, envir = envir)
 
       expr <- mod$get.coalescent.rate(comp)
       rate <- eval(parse(text = expr), envir = envir)

From 5dea260713eaba1579ebdfe03cc0a292f5b25a40 Mon Sep 17 00:00:00 2001
From: WongDWai <dwong338@uwo.ca>
Date: Fri, 28 Aug 2026 10:47:07 -0400
Subject: [PATCH 6/8] adding explicit pop.size in yamls.

---
 tests/testthat/test_Superinfection.yaml | 1 +
 tests/testthat/test_simInnerTree.R      | 1 +
 tests/testthat/test_superinfection.R    | 2 +-
 3 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/tests/testthat/test_Superinfection.yaml b/tests/testthat/test_Superinfection.yaml
index fc96ce3..7e329a7 100644
--- a/tests/testthat/test_Superinfection.yaml
+++ b/tests/testthat/test_Superinfection.yaml
@@ -24,6 +24,7 @@ Compartments:
     size: 2
     bottleneck.size: 1
     coalescent.rate: 0.01
+    pop.size: 100    # was silently defaulting to this -- confirm intended within-host pop size with Art
   I_samp:
     infected: true
     size: 0
diff --git a/tests/testthat/test_simInnerTree.R b/tests/testthat/test_simInnerTree.R
index e6c147e..22e6880 100644
--- a/tests/testthat/test_simInnerTree.R
+++ b/tests/testthat/test_simInnerTree.R
@@ -3,6 +3,7 @@ require(twt)
 # generate test fixtures
 settings <- yaml.load_file("test_SIR.yaml")
 settings$Compartments$I$coalescent.rate <- 1.0
+settings$Compartments$I$pop.size <- 100  # was silently defaulting to this
 mod <- Model$new(settings)
 set.seed(276)
 dynamics <- sim.dynamics(mod)
diff --git a/tests/testthat/test_superinfection.R b/tests/testthat/test_superinfection.R
index 4d79ceb..6662525 100644
--- a/tests/testthat/test_superinfection.R
+++ b/tests/testthat/test_superinfection.R
@@ -393,7 +393,7 @@ test_that("model: compartment I flagged infected=TRUE", {
   expect_true(mod$get.infected("I"))
 })
 
-test_that("model YAML: compartment I has pop.size=2 and bottleneck.size=1", {
+test_that("model YAML: compartment I has size=2 and bottleneck.size=1", {
   try_model(SUPERINF_INNER_PATH)
   cfg <- read_yaml(SUPERINF_INNER_PATH)
   expect_equal(cfg$Compartments$I$size, 2)

From cbb6e89f66ce9c3ff5b23b4735a7503775d20621 Mon Sep 17 00:00:00 2001
From: WongDWai <80217512+WongDWai@users.noreply.github.com>
Date: Fri, 28 Aug 2026 12:58:06 -0400
Subject: [PATCH 7/8] Update simARG.R

---
 R/simARG.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/simARG.R b/R/simARG.R
index 02fa2df..b7feb4a 100644
--- a/R/simARG.R
+++ b/R/simARG.R
@@ -198,7 +198,7 @@ sim.arg <- function(outer, rho = 1e-4, seq.length = 9000L) {
   active <- inner$get.active()
   host   <- active$get.host.by.name(host.name)
 
-  # lazily initialize this host's FIXED lineage pool (Art's fix: sample
+  # initialize this host's FIXED lineage pool (sample
   # recombination parents from a fixed pool of p.size lineages, only
   # activating a previously-inactive one when chosen, instead of always
   # creating a brand-new lineage de novo. Total pool size never changes.)

From d7d73b8422ce4fae49c1a7a3bf5267135d77df26 Mon Sep 17 00:00:00 2001
From: WongDWai <80217512+WongDWai@users.noreply.github.com>
Date: Fri, 28 Aug 2026 12:59:12 -0400
Subject: [PATCH 8/8] Update test_Superinfection.yaml

---
 tests/testthat/test_Superinfection.yaml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/testthat/test_Superinfection.yaml b/tests/testthat/test_Superinfection.yaml
index 7e329a7..1f2db63 100644
--- a/tests/testthat/test_Superinfection.yaml
+++ b/tests/testthat/test_Superinfection.yaml
@@ -24,7 +24,7 @@ Compartments:
     size: 2
     bottleneck.size: 1
     coalescent.rate: 0.01
-    pop.size: 100    # was silently defaulting to this -- confirm intended within-host pop size with Art
+    pop.size: 100    # was silently defaulting to this
   I_samp:
     infected: true
     size: 0