fix: superposition converges instead of hanging when tlast < tau with auc.type "AUClast"

NEWS.md

@@ -6,6 +6,19 @@ the dosing including dose amount and route.
 
 # Development version
 
+* Bug fix: `superposition()` no longer loops forever when steady-state is
+  requested (the default `n.tau=Inf`) and concentrations are extrapolated as
+  zero after `tlast` (for example, `auc.type="AUClast"` when `tlast < tau`).
+  Steady-state is now assessed on the concentrations that can become nonzero,
+  a warning is given when zero concentrations remain in the steady-state
+  profile, and an informative error (instead of an infinite loop) is raised if
+  steady-state cannot be reached within 10000 dosing intervals (#580).
+  Warnings, messages, and printed output from `superposition()` on a
+  `PKNCAconc` object are now raised by the parent process.  Each subject runs
+  through `parallel::mclapply()`, which forks on every platform except
+  Windows, and a condition raised in a forked worker never reaches the
+  caller, so the new warning was previously invisible except on Windows.
+
 * Corrected registry description strings that did not match the
   implementation (#582).  `adj.r.squared.factor` is described as selecting
   the regression with the most points among those within the tolerance of

R/superposition.R

@@ -35,6 +35,12 @@
 #' times: 0 (zero), `dose.times`, `time modulo tau` (shifting `time` for each
 #' dose time as well), `additional.times`, and `tau`.
 #'
+#' When concentrations are extrapolated as zero after `tlast` (for example,
+#' `auc.type = "AUClast"` when `tlast < tau`), some concentrations in the
+#' profile can never become nonzero.  Steady-state (`n.tau = Inf`) is then
+#' assessed on the nonzero concentrations, and a warning is given that zero
+#' concentrations remain in the steady-state profile.
+#'
 #' @seealso [interp.extrap.conc()]
 #' @family Superposition
 #' @export
@@ -48,19 +54,38 @@ superposition.PKNCAconc <- function(conc, ...) {
   # Split the data by grouping and extract just the concentration and
   # time columns
   nested_data <- prepare_PKNCAconc(conc)
+  # parallel::mclapply() forks on every platform except Windows, and a
+  # condition raised in a forked worker never reaches the parent, so collect
+  # each worker's conditions with the result and re-emit them below.
+  # purrr::quietly() also captures messages and printed output, so those are
+  # re-emitted as well; otherwise adding a message() to the calculation later
+  # would silently discard it on every platform.
+  quiet_superposition <- purrr::quietly(superposition.numeric)
   tmp_results <-
     parallel::mclapply(
       X=seq_len(nrow(nested_data)),
       FUN=function(idx) {
-        superposition.numeric(
+        quiet_superposition(
           conc=nested_data$data_conc[[idx]]$conc,
           time=nested_data$data_conc[[idx]]$time,
           ...
         )
       }
     )
+  for (current_warning in unlist(lapply(tmp_results, FUN="[[", "warnings"))) {
+    warning(current_warning, call.=FALSE)
+  }
+  # `messages` keep their trailing newline; `output` has it stripped.
+  for (current_message in unlist(lapply(tmp_results, FUN="[[", "messages"))) {
+    message(current_message, appendLF=FALSE)
+  }
+  current_output <- unlist(lapply(tmp_results, FUN="[[", "output"))
+  current_output <- current_output[nzchar(current_output)]
+  if (length(current_output) > 0) {
+    writeLines(current_output)
+  }
   # Replace the concentration data with the new results
-  nested_data$data_conc <- tmp_results
+  nested_data$data_conc <- lapply(tmp_results, FUN="[[", "result")
   tidyr::unnest(nested_data, cols="data_conc")
 }
 
@@ -212,10 +237,22 @@ superposition.numeric <- function(conc, time, dose.input = NULL,
     # Do the math! (Finally)
     current.tol <- steady.state.tol + 1
     tau.count <- 0
+    # Generous backstop so that a case that cannot reach steady-state errors
+    # instead of looping forever (normal convergence takes far fewer
+    # intervals; the default steady.state.tol needs at most ~1000 even for
+    # negligible accumulation decay).
+    max.tau.count <- 10000L
     # Stop either for reaching steady-state or for reaching the requested number of doses
     while (tau.count < n.tau &
            !is.na(current.tol) &
            current.tol >= steady.state.tol) {
+      if (is.infinite(n.tau) && tau.count >= max.tau.count) {
+        stop(
+          "Superposition did not reach steady-state within ", max.tau.count,
+          " dosing intervals (steady.state.tol=", steady.state.tol,
+          ").  Check the lambda.z, tau, and steady.state.tol inputs, or use a finite n.tau."
+        )
+      }
       prev.conc <- ret$conc
       # Perform the dosing for a single dosing interval.
       for (i in seq_along(dose.times)) {
@@ -238,14 +275,29 @@ superposition.numeric <- function(conc, time, dose.input = NULL,
           )
       }
       tau.count <- tau.count + 1
-      if (any(ret$conc %in% 0)) {
-        # prevent division by 0.  Since not all concentrations are 0,
-        # all values will eventually be nonzero.
-        current.tol <- steady.state.tol + 1
-      } else {
+      zero.conc <- ret$conc %in% 0
+      if (!any(zero.conc)) {
         current.tol <- max(1-(prev.conc/ret$conc))
+      } else if ((tau*(tau.count - 1) - max(dose.times)) > tlast) {
+        # The dosing interval just added drew all of its concentrations from
+        # times after tlast, where extrapolation may be zero (e.g. with
+        # auc.type="AUClast").  A concentration that is still zero can
+        # therefore never become nonzero (a structural zero), so assess
+        # steady-state on the nonzero concentrations only.
+        current.tol <- max(1-(prev.conc[!zero.conc]/ret$conc[!zero.conc]))
+      } else {
+        # A zero concentration may still become nonzero from a later dosing
+        # interval; prevent division by 0 and do not yet assess steady-state.
+        current.tol <- steady.state.tol + 1
       }
     }
+    if (!is.na(current.tol) && current.tol < steady.state.tol && any(ret$conc %in% 0)) {
+      warning(
+        "Zero concentrations remain in the steady-state superposition profile.  ",
+        "They come from concentrations extrapolated as zero after tlast (",
+        signif(tlast, 6), ") with auc.type=", dQuote(auc.type, q = FALSE), "."
+      )
+    }
   }
   ret
 }

tests/testthat/test-superpostion.R

@@ -514,3 +514,154 @@ test_that("PKNCAconc superposition", {
     )
   )
 })
+
+test_that("superposition reaches steady-state despite structural zero concentrations (issue 580)", {
+  # Theoph subjects 6 and 10 have tlast < tau=24, so with auc.type="AUClast"
+  # the concentrations after tlast extrapolate as exactly zero and stay zero at
+  # steady-state.  This formerly looped forever with n.tau=Inf because the
+  # zero-concentration guard reset the steady-state tolerance every interval.
+  d_theoph_6 <- datasets::Theoph[datasets::Theoph$Subject == 6, ]
+  expect_warning(
+    v_ss <-
+      superposition(
+        conc=d_theoph_6$conc, time=d_theoph_6$Time, tau=24, auc.type="AUClast"
+      ),
+    regexp='Zero concentrations remain in the steady-state superposition profile.  They come from concentrations extrapolated as zero after tlast (23.85) with auc.type="AUClast".',
+    fixed=TRUE
+  )
+  # No accumulation is possible because every concentration one or more
+  # dosing intervals after dosing is zero, so the steady-state profile is
+  # exactly the single-dose profile with structural zeros at times 0 and tau.
+  expect_identical(
+    v_ss,
+    data.frame(conc=c(d_theoph_6$conc, 0), time=c(d_theoph_6$Time, 24))
+  )
+  # A finite n.tau with the same input gave the same concentrations before the
+  # fix (all added dosing intervals contribute zero); that is preserved, and
+  # partial accumulation with a finite n.tau does not warn.
+  expect_silent(
+    v_1 <-
+      superposition(
+        conc=d_theoph_6$conc, time=d_theoph_6$Time, tau=24, n.tau=1, auc.type="AUClast"
+      )
+  )
+  expect_identical(v_ss, v_1)
+
+  # The default auc.type="AUCinf" on the same subject is unaffected by the fix:
+  # all concentrations become positive and steady-state is reached silently
+  # (values are pinned from before the fix).
+  expect_silent(
+    v_aucinf <- superposition(conc=d_theoph_6$conc, time=d_theoph_6$Time, tau=24)
+  )
+  expect_equal(
+    v_aucinf,
+    data.frame(
+      conc=c(1.03361737415698, 2.29940375348928, 4.06230162309106,
+             7.37435349523853, 7.18488330467035, 6.28550707796221,
+             5.60636744789876, 4.57905606419521, 3.92005369868033,
+             3.13726988893256, 1.04734393070531, 1.03364150451656),
+      time=c(0, 0.27, 0.58, 1.15, 2.03, 3.57, 5, 7, 9.22, 12.1, 23.85, 24)
+    )
+  )
+
+  # AUClast with tlast > tau converges without a warning: the time-zero
+  # concentration is zero after the first dosing interval, becomes nonzero on
+  # the second, and no zeros remain at steady-state.
+  d_theoph_1 <- datasets::Theoph[datasets::Theoph$Subject == 1, ]
+  expect_silent(
+    v_s1 <-
+      superposition(
+        conc=d_theoph_1$conc, time=d_theoph_1$Time, tau=24, auc.type="AUClast",
+        check.blq=FALSE
+      )
+  )
+  expect_true(all(v_s1$conc > 0))
+
+  # The PKNCAconc method (the form reported in issue 580) returns with one
+  # warning per affected subject (check.blq=FALSE because subject 10 has a
+  # nonzero concentration at time 0)
+  d_theoph_6_10 <-
+    datasets::Theoph[datasets::Theoph$Subject %in% c(6, 10), ]
+  conc_obj <- PKNCAconc(conc~Time|Subject, data=d_theoph_6_10)
+  w_obj <-
+    capture_warnings(
+      v_obj <- superposition(conc_obj, tau=24, auc.type="AUClast", check.blq=FALSE)
+    )
+  expect_length(w_obj, 2)
+  expect_match(
+    w_obj,
+    regexp="Zero concentrations remain in the steady-state superposition profile",
+    fixed=TRUE,
+    all=TRUE
+  )
+  expect_true(all(is.finite(v_obj$conc)))
+})
+
+test_that("superposition errors instead of hanging when steady-state cannot be reached", {
+  # A negligible lambda.z makes each added dosing interval contribute ~clast,
+  # so the relative change shrinks only like 1/n.tau and steady.state.tol=1e-8
+  # cannot be reached within the iteration cap.
+  expect_error(
+    superposition(
+      conc=c(0, 1, 0.5), time=0:2, tau=1,
+      lambda.z=1e-8, clast.pred=0.5, tlast=2,
+      steady.state.tol=1e-8
+    ),
+    regexp="Superposition did not reach steady-state within 10000 dosing intervals",
+    fixed=TRUE
+  )
+})
+
+test_that("superposition warns about structural zeros for a single profile (issue 580)", {
+  # superposition.PKNCAconc() runs each subject through parallel::mclapply(),
+  # which forks everywhere except Windows, and a warning raised in a forked
+  # worker never reaches the parent; the method collects them with
+  # purrr::quietly() and re-emits them.  The PKNCAconc test above cannot show
+  # that on Windows, where mclapply() falls back to lapply(), so pin the
+  # underlying warning here where it is raised.
+  d_theoph_6 <- datasets::Theoph[datasets::Theoph$Subject == 6, ]
+  expect_warning(
+    v_6 <-
+      superposition(
+        conc=d_theoph_6$conc, time=d_theoph_6$Time,
+        tau=24, auc.type="AUClast", check.blq=FALSE
+      ),
+    regexp="Zero concentrations remain in the steady-state superposition profile",
+    fixed=TRUE
+  )
+  expect_equal(nrow(v_6), 12)
+  expect_true(all(is.finite(v_6$conc)))
+})
+
+test_that("superposition.PKNCAconc re-emits everything its workers produce (issue 580)", {
+  # parallel::mclapply() forks everywhere except Windows, and conditions raised
+  # in a forked worker never reach the parent.  purrr::quietly() collects the
+  # warnings, messages, and printed output so the method can re-emit them.
+  # Mock the calculation so all three are produced regardless of the data.
+  d_theoph_6_10 <-
+    datasets::Theoph[datasets::Theoph$Subject %in% c(6, 10), ]
+  conc_obj <- PKNCAconc(conc~Time|Subject, data=d_theoph_6_10)
+  local_mocked_bindings(
+    superposition.numeric=function(conc, time, ...) {
+      warning("mocked warning")
+      message("mocked message")
+      cat("mocked output", fill=TRUE)
+      data.frame(conc=c(1, 2), time=c(0, 1))
+    }
+  )
+  printed <- NULL
+  warned <-
+    capture_warnings(
+      messaged <-
+        capture_messages(
+          # Assign the result so capture.output() does not also see it printed.
+          printed <- capture.output(v_obj <- superposition(conc_obj, tau=24))
+        )
+    )
+  # One of each per subject, and messages are not given an extra blank line.
+  expect_length(warned, 2)
+  expect_equal(unique(warned), "mocked warning")
+  expect_length(messaged, 2)
+  expect_equal(unique(trimws(messaged)), "mocked message")
+  expect_equal(unique(printed), "mocked output")
+})