nlmixr2 · mattfidler · Sep 18, 2023 · Sep 18, 2023
diff --git a/NAMESPACE b/NAMESPACE
@@ -86,6 +86,7 @@ S3method(print,rxSolveSimType)
 S3method(print,rxSymInvCholEnv)
 S3method(print,rxUi)
 S3method(print,rxode2)
+S3method(rbind,rxSolve)
 S3method(rxChain2,EventTable)
 S3method(rxChain2,default)
 S3method(rxCompile,character)

diff --git a/R/rbind.R b/R/rbind.R
@@ -0,0 +1,106 @@
+.rbind2rxSove <- function(v, v2) {
+  .tmp <- v$add.dosing # make sure info is calculated
+  .tmp <- v2$add.dosing
+  .cls <- class(v)
+  .env1 <- attr(class(v), ".rxode2.env")
+  .env2 <- attr(class(v2), ".rxode2.env")
+  if (length(.env1$.check.names) != length(.env2$.check.names) &&
+        !all(.env1$.check.names == .env2$.check.names)) {
+    stop("cannot rbind these 2 rxSolve objects", call.=FALSE)
+  }
+  ## if (is.null(.env1$.params.single) || is.null(.env2$.params.single)) {
+  ##   stop("cannot rbind single solve environments", call.=FALSE)
+  ## }
+  .cloneEnv <- new.env(parent=emptyenv())
+  for (.v in ls(.env1, all=TRUE)) {
+    assign(.v, get(.v, envir=.env1), envir=.cloneEnv)
+  }
+  .nStud1 <- .cloneEnv$.args$nStud
+  .nStud2 <- .env2$.args$nStud
+  if (.nStud1 > 1 && .nStud2 > 1) {
+    .args <- .cloneEnv$.args
+    .args$nStud <- .nStud1 + .nStud2
+    .cloneEnv$.args <- .args
+    .v <- v
+    class(.v) <- "data.frame"
+    .v2 <- v2
+    .v2$sim.id <- .nStud1 + .v2$sim.id
+    class(.v2) <- "data.frame"
+    .v <- rbind(.v, .v2)
+    .cloneEnv$.nsub <- .cloneEnv$.nsub + .env2$.nsub
+    .cloneEnv$.et <- NULL
+    .cloneEnv$.args.params <- NULL
+    .cloneEnv$.args.inits <- NULL
+    .cloneEnv$.args.object <- rxModelVars(.cloneEnv$.args.object)
+    .cloneEnv$.args.par0 <- NULL
+    .cloneEnv$.args.params <- NULL
+    .cloneEnv$.check.nrow <- .cloneEnv$.check.nrow + .env2$.check.nrow
+    .cloneEnv$nobs <- .cloneEnv$nobs + .env2$nobs
+    .cloneEnv$.nsim <- .cloneEnv$.nsim + .env2$.nsim
+    .cloneEnv$.nsub <- NULL
+    .cloneEnv$.dadt.counter <- 0L
+    .cloneEnv$.init.dat <- setNames(rep(NA_real_, length(.cloneEnv$.init.dat)),
+                                    names(.cloneEnv$.init.dat))
+    .cloneEnv$.jac.counter <- 0L
+    .cloneEnv$.nsub <- NA_integer_
+    .cloneEnv$.par.pos <- NULL
+    .cloneEnv$.par.pos.ini <- NULL
+    .pd <- .env2$.params.dat
+    .pd$sim.id <- .nStud1 + .pd$sim.id
+
+    .cloneEnv$.params.dat <- rbind(.cloneEnv$.params.dat, .pd)
+    .cloneEnv$.params.single <- NULL
+    .cloneEnv$.par.pos <- NULL
+    .cloneEnv$.par.pos.ini <- NULL
+    .cloneEnv$.et <- NULL
+    if (inherits(.v$id, "factor")) {
+      .cloneEnv$.idLevels <- levels(.v$id)
+    }
+    .cloneEnv$.jac.counter <- .cloneEnv$.slvr.counter <- 0L
+    .cloneEnv$.nsim <- .cloneEnv$.nsim + .env2$.nsim
+    .cloneEnv$.real.update <- FALSE
+    .cloneEnv$.sigma <- NULL
+    .fun <- function(...) {
+      stop("functions don't work on rbound rxSolve", call.=FALSE)
+    }
+    .cloneEnv$.replace.sampling <- .fun
+    .cloneEnv$add.dosing <- .fun
+    .cloneEnv$add.sampling <- .fun
+    .cloneEnv$clear.dosing <- .fun
+    .cloneEnv$clear.sampling <- .fun
+    .cloneEnv$get.dosing <- .fun
+    .cloneEnv$get.EventTable <- .fun
+    .cloneEnv$get.nobs <- .fun
+    .cloneEnv$get.obs.rec <- .fun
+    .cloneEnv$get.sampling <- .fun
+    .cloneEnv$get.units <- .fun
+    .cloneEnv$import <- .fun
+    .cloneEnv$counts.EventTable <- NULL
+    .cloneEnv$get.units <- NULL
+    .cloneEnv$units <- c(dosing="NA", time="NA")
+    .cloneEnv$dll <- NULL
+    attr(.cls, ".rxode2.env") <- .cloneEnv
+    class(.v) <- .cls
+    return(.v)
+  }
+}
+
+#' @export
+rbind.rxSolve <- function(..., deparse.level = 1) {
+  .lst <- list(...)
+  if (length(.lst) >= 2) {
+    .ret <- .rbind2rxSove(.lst[[1]], .lst[[2]])
+    if (length(.lst) == 2) {
+      return(.ret)
+    }
+    return(do.call(rbind.rxSolve,
+                   c(list(.ret),
+                     lapply(seq_along(.lst)[-(1:2)],
+                            function(i){
+                              .lst[[i]]
+                            }))))
+  }
+  if (length(.lst) == 1) return(.lst[[1]])
+  stop("called rbind.rxSolve() with no arguments",
+       call.=FALSE)
+}
diff --git a/R/rxsolve.R b/R/rxsolve.R
@@ -1,4 +1,4 @@
- #' Options, Solving & Simulation of an ODE/solved system
+#' Options, Solving & Simulation of an ODE/solved system
 #'
 #' This uses rxode2 family of objects, file, or model specification to
 #' solve a ODE system.  There are many options for a solved rxode2
@@ -1212,8 +1212,8 @@ rxSolve.function <- function(object, params = NULL, events = NULL, inits = NULL,
     .rx <- object$simulationModel
   }
   list(list(object=.rx, params = params, events = events, inits = inits),
-                       .rxControl,
-                       list(theta = theta, eta = eta))
+       .rxControl,
+       list(theta = theta, eta = eta))
 }
 
 #' @rdname rxSolve
@@ -1238,7 +1238,7 @@ rxSolve.rxUi <- function(object, params = NULL, events = NULL, inits = NULL, ...
   if (is.null(.lst$omega) && is.null(.lst$sigma)) {
     .pred <- TRUE
     if (!.hasIpred && any(rxModelVars(.lst[[1]])$lhs == "ipredSim")) {
-      .lst$drop <- c(.lst$drop, "ipredSim")      
+      .lst$drop <- c(.lst$drop, "ipredSim")
     }
   }
   .ret <- do.call("rxSolve.default", .lst)
@@ -1301,6 +1301,9 @@ rxSolve.nlmixr2FitCore <- rxSolve.nlmixr2FitData
 #' @export
 rxSolve.default <- function(object, params = NULL, events = NULL, inits = NULL, ...,
                             theta = NULL, eta = NULL) {
+  if (inherits(object, "rxSolve") && is.null(.v$env$.ev)) {
+    stop("cannot solve from more than one rxSolve object bound together", call.=FALSE)
+  }
   on.exit({
     .clearPipe()
     .asFunctionEnv$rx <- NULL
@@ -1310,13 +1313,13 @@ rxSolve.default <- function(object, params = NULL, events = NULL, inits = NULL,
   if (rxIs(object, "rxEt")) {
     if (!is.null(events)) {
       stop("events can be pipeline or solving arguments not both",
-        call. = FALSE
-      )
+           call. = FALSE
+           )
     }
     if (is.null(rxode2et::.pipeRx(NA))) {
       stop("need an rxode2 compiled model as the start of the pipeline",
-        call. = FALSE
-      )
+           call. = FALSE
+           )
     } else {
       events <- object
       object <- rxode2et::.pipeRx(NA)
@@ -1328,16 +1331,16 @@ rxSolve.default <- function(object, params = NULL, events = NULL, inits = NULL,
     }
     if (is.null(rxode2et::.pipeRx(NA))) {
       stop("need an rxode2 compiled model as the start of the pipeline",
-        call. = FALSE
-      )
+           call. = FALSE
+           )
     } else {
       .rxParams <- object
       object <- rxode2et::.pipeRx(NA)
     }
     if (is.null(rxode2et::.pipeEvents(NA))) {
       stop("need an rxode2 events as a part of the pipeline",
-        call. = FALSE
-      )
+           call. = FALSE
+           )
     } else {
       events <- rxode2et::.pipeEvents(NA)
       rxode2et::.pipeEvents(NULL)
@@ -1347,29 +1350,29 @@ rxSolve.default <- function(object, params = NULL, events = NULL, inits = NULL,
     events <- rxode2et::.pipeEvents(NA)
   } else if (!is.null(rxode2et::.pipeEvents(NA)) && !is.null(events)) {
     stop("'events' in pipeline AND in solving arguments, please provide just one",
-      call. = FALSE
-    )
+         call. = FALSE
+         )
   } else if (!is.null(rxode2et::.pipeEvents(NA)) && !is.null(params) &&
-    rxIs(params, "event.data.frame")) {
+               rxIs(params, "event.data.frame")) {
     stop("'events' in pipeline AND in solving arguments, please provide just one",
-      call. = FALSE
-    )
+         call. = FALSE
+         )
   }
 
   if (!is.null(rxode2et::.pipeParams(NA)) && is.null(params)) {
     params <- rxode2et::.pipeParams(NA)
   } else if (!is.null(rxode2et::.pipeParams(NA)) && !is.null(params)) {
     stop("'params' in pipeline AND in solving arguments, please provide just one",
-      call. = FALSE
-    )
+         call. = FALSE
+         )
   }
 
   if (!is.null(rxode2et::.pipeInits(NA)) && is.null(inits)) {
     inits <- rxode2et::.pipeInits(NA)
   } else if (!is.null(rxode2et::.pipeInits(NA)) && !is.null(inits)) {
     stop("'inits' in pipeline AND in solving arguments, please provide just one",
-      call. = FALSE
-    )
+         call. = FALSE
+         )
   }
 
   if (.applyParams) {
@@ -1380,13 +1383,13 @@ rxSolve.default <- function(object, params = NULL, events = NULL, inits = NULL,
   .xtra <- list(...)
   if (any(duplicated(names(.xtra)))) {
     stop("duplicate arguments do not make sense",
-      call. = FALSE
-    )
+         call. = FALSE
+         )
   }
   if (any(names(.xtra) == "covs")) {
     stop("covariates can no longer be specified by 'covs'\n  include them in the event dataset\n\nindividual covariates: Can be specified by a 'iCov' dataset\n each each individual covariate has a value\n\ntime varying covariates: modify input event data-frame or\n  'eventTable' to include covariates(https://tinyurl.com/y52wfc2y)\n\nEach approach needs the covariates named to match the variable in the model",
-      call. = FALSE
-    )
+         call. = FALSE
+         )
   }
   .nms <- names(as.list(match.call())[-1])
   .lst <- list(...)
@@ -1400,9 +1403,9 @@ rxSolve.default <- function(object, params = NULL, events = NULL, inits = NULL,
     .both <- intersect(.mv$params, .ctl$keep)
     if (length(.both) > 0) {
       .keep <- .ctl$keep[!(.ctl$keep %in% .both)]
-       if (length(.keep) == 0L) {
-          .keep <- NULL
-       }
+      if (length(.keep) == 0L) {
+        .keep <- NULL
+      }
       .w <- which(names(.ctl) == "keep")
       .ctl[[.w]] <- .keep
       .ctl <- do.call(rxControl,
@@ -1536,8 +1539,8 @@ rxSolve.default <- function(object, params = NULL, events = NULL, inits = NULL,
       events <- c(.theta, .eta)
     } else {
       stop("cannot specify 'params' and 'theta'/'eta' at the same time",
-        call. = FALSE
-      )
+           call. = FALSE
+           )
     }
   }
   if (!is.null(.ctl$iCov)) {

diff --git a/man/meanProbs.Rd b/man/meanProbs.Rd
diff --git a/tests/testthat/test-rbind.R b/tests/testthat/test-rbind.R
@@ -0,0 +1,60 @@
+test_that("rbind", {
+
+  mod2 <- function() {
+    ini({
+      KA <- 2.94E-01
+      TCL <- 1.86E+01
+      V2 <- 4.02E+01
+      Q <- 1.05E+01
+      V3 <- 2.97E+02
+      Kin <- 1
+      Kout <- 1
+      EC50 <- 200
+      eta.Cl ~ 0.2
+      eff.sd <- sqrt(0.05)
+      c2.sd <- sqrt(0.05)
+    })
+    model({
+      CL <- TCL * exp(eta.Cl)
+      d/dt(depot) <- -KA * depot
+      d/dt(centr) <- KA * depot - CL * (centr / V2) - Q * (centr / V2) + Q * (peri / V3)
+      d/dt(peri) <- Q * (centr / V2) - Q * (peri / V3)
+      C2 <- centr/V2
+      d/dt(eff) <- Kin - Kout * (1 - C2 / (EC50 + C2)) * eff
+      eff(0) <- 1000
+      eff ~ add(eff.sd)
+      C2 ~ prop(c2.sd)
+    })
+  }
+
+  ev <- eventTable() %>%
+    add.dosing(dose = 10000, nbr.doses = 10, dosing.interval = 12, dosing.to = 2) %>%
+    add.dosing(dose = 20000, nbr.doses = 5, start.time = 120, dosing.interval = 24, dosing.to = 2) %>%
+    add.sampling(0:240) %>%
+    et(id=1:2) %>%
+    as.data.frame()
+
+  ev$dvid <- ifelse(ev$evid == 0, 1, NA_real_)
+
+  ev2 <- ev[which(ev$dvid == 1), ]
+  ev2$dvid <- 2
+
+  ev <- rbind(ev, ev2) %>%
+    dplyr::arrange(id, time)
+
+  rownames(ev) <- NULL
+
+  pk1 <-
+    suppressWarnings(rxSolve(mod2, nStud = 4, ev, cores = 2))
+
+  pk2 <-
+    suppressWarnings(rxSolve(mod2, nStud = 4, ev, cores = 2))
+
+  b1 <- rbind.rxSolve(pk1, pk2)
+
+  b2 <- rbind.rxSolve(pk1, pk2, pk1)
+
+  expect_true(inherits(b1, "rxSolve"))
+  expect_true(inherits(b2, "rxSolve"))
+  expect_error(rbind.rxSolve())
+})