diff --git a/NEWS.md b/NEWS.md index d28b21673..6e8c9248f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -41,6 +41,10 @@ the same way as `liblsoda`, so the random number provided will be the same with different solving methods. +- The arguments saved in the `rxSolve` for items like `thetaMat` will + be the reduced matrices used in solving, not the full matrices (this + will likely not break very many items) + ## Possible breaking changes (though unlikely) - `iCov` is no longer merged to the event dataset. This makes solving diff --git a/R/rxsolve.R b/R/rxsolve.R index 12d9ac725..00db8f555 100644 --- a/R/rxsolve.R +++ b/R/rxsolve.R @@ -1764,16 +1764,9 @@ rxSolve.default <- function(object, params = NULL, events = NULL, inits = NULL, .envReset$unload <- FALSE # take care of too many DLLs or not provided simulation errors .names <- NULL - if (inherits(.ctl$thetaMat, "matrix")) { - .mv <- rxModelVars(object) - .col <- colnames(.ctl$thetaMat) - .w <- .col %in% .mv$params - .ignore <- .col[!.w] - if (length(.ignore)>0) { - .minfo(paste0("thetaMat has too many items, ignored: '", paste(.ignore, collapse="', '"), "'")) - } - .names <- c(.names, .col[.w]) - } + + .extraNames <- character(0) + if (inherits(.ctl$omega, "matrix")) { .mv <- rxModelVars(object) .col <- colnames(.ctl$omega) @@ -1782,8 +1775,14 @@ rxSolve.default <- function(object, params = NULL, events = NULL, inits = NULL, if (length(.ignore)>0) { .minfo(paste0("omega has too many items, ignored: '", paste(.ignore, collapse="', '"), "'")) } + .ctl$omega <-.ctl$omega[.w, .w, drop=FALSE] + if (dim(.ctl$omega)[1] == 0) { + .ctl$omega <- NULL + .ctl <- do.call(rxControl, .ctl) + } .names <- c(.names, .col[.w]) } else if ( inherits(.ctl$omega, "character")) { + .extraNames <- c(.extraNames, .ctl$omega) .mv <- rxModelVars(object) .col <- .ctl$omega .w <- .col %in% .mv$params @@ -1801,8 +1800,14 @@ rxSolve.default <- function(object, params = NULL, events = NULL, inits = NULL, if (length(.ignore)>0) { .minfo(paste0("sigma has too many items, ignored: '", paste(.ignore, collapse="', '"), "'")) } + .ctl$sigma <-.ctl$sigma[.w, .w, drop=FALSE] + if (dim(.ctl$sigma)[1] == 0) { + .ctl$sigma <- NULL + .ctl <- do.call(rxControl, .ctl) + } .names <- c(.names, .col[.w]) } else if ( inherits(.ctl$sigma, "character")) { + .extraNames <- c(.extraNames, .ctl$sigma) .mv <- rxModelVars(object) .col <- .ctl$sigma .w <- .col %in% .mv$params @@ -1812,6 +1817,36 @@ rxSolve.default <- function(object, params = NULL, events = NULL, inits = NULL, } .names <- c(.names, .col[.w]) } + + if (inherits(.ctl$thetaMat, "matrix")) { + .mv <- rxModelVars(object) + .col <- colnames(.ctl$thetaMat) + .w <- .col %in% c(.mv$params, .extraNames) + .ignore <- .col[!.w] + if (length(.ignore)>0) { + .minfo(paste0("thetaMat has too many items, ignored: '", paste(.ignore, collapse="', '"), "'")) + } + .ctl$thetaMat <-.ctl$thetaMat[.w, .w, drop=FALSE] + if (dim(.ctl$thetaMat)[1] == 0) { + .ctl$thetaMat <- NULL + .ctl <- do.call(rxControl, .ctl) + } + .names <- c(.names, .col[.w]) + + # now look for zero diagonals + .col <- colnames(.ctl$thetaMat) + .d <- diag(.ctl$thetaMat) + .w <- which(.d == 0) + if (length(.w) > 0) { + .minfo(paste0("thetaMat has zero diagonal items, ignored: '", paste(.col[.w], collapse="', '"), "'")) + .ctl$thetaMat <-.ctl$thetaMat[-.w, -.w, drop=FALSE] + if (dim(.ctl$thetaMat)[1] == 0) { + .ctl$thetaMat <- NULL + .ctl <- do.call(rxControl, .ctl) + } + .names <- c(.names, .col[-.w]) + } + } rxSetCovariateNamesForPiping(NULL) if (length(.ctl$.zeros) > 0) { if (rxIs(params, "rx.event")) { diff --git a/tests/testthat/test-par-solve.R b/tests/testthat/test-par-solve.R index 867a7d81b..b160610fa 100644 --- a/tests/testthat/test-par-solve.R +++ b/tests/testthat/test-par-solve.R @@ -1,4 +1,5 @@ rxTest({ + expectSimWithout <- function(x) { expect_null(x$thetaMat) expect_null(x$omegaList) @@ -27,8 +28,8 @@ rxTest({ # context(sprintf("Test Parallel/Multi-subject Solve (%s)", meth)) mod <- rxode2({ - d / dt(intestine) <- -a * intestine - d / dt(blood) <- a * intestine - b * blood + d/dt(intestine) <- -a * intestine + d/dt(blood) <- a * intestine - b * blood }) et <- eventTable(time.units = "days") @@ -118,10 +119,10 @@ rxTest({ C2 <- prod(centr, 1 / V2) C3 ~ prod(peri, 1 / V3) CL ~ prod(TCL, exp(eta.Cl)) - d / dt(depot) ~ prod(-KA, depot) - d / dt(centr) ~ sum(prod(KA, depot), -prod(CL, C2), -prod(Q, C2), prod(Q, C3)) - d / dt(peri) ~ sum(prod(Q, C2), -prod(Q, C3)) - d / dt(eff) <- sum(Kin, -prod(Kout, sum(1, -prod(C2, 1 / sum(EC50, C2))), eff)) + d/dt(depot) ~ prod(-KA, depot) + d/dt(centr) ~ sum(prod(KA, depot), -prod(CL, C2), -prod(Q, C2), prod(Q, C3)) + d/dt(peri) ~ sum(prod(Q, C2), -prod(Q, C3)) + d/dt(eff) <- sum(Kin, -prod(Kout, sum(1, -prod(C2, 1 / sum(EC50, C2))), eff)) e1 <- err1 e2 <- err2 resp <- sum(eff, e1) @@ -159,10 +160,10 @@ rxTest({ C2 <- centr / V2 C3 ~ peri / V3 CL ~ TCL * exp(eta.Cl) - d / dt(depot) ~ -KA * depot - d / dt(centr) ~ KA * depot - CL * C2 - Q * C2 + Q * C3 - d / dt(peri) ~ Q * C2 - Q * C3 - d / dt(eff) <- Kin - Kout * (1 - C2 / (EC50 + C2)) * eff + d/dt(depot) ~ -KA * depot + d/dt(centr) ~ KA * depot - CL * C2 - Q * C2 + Q * C3 + d/dt(peri) ~ Q * C2 - Q * C3 + d/dt(eff) <- Kin - Kout * (1 - C2 / (EC50 + C2)) * eff eff(0) <- 1000 e1 <- err1 e2 <- err2 diff --git a/tests/testthat/test-rxode-issue-349.R b/tests/testthat/test-rxode-issue-349.R index f23f3be0c..22a636ca1 100644 --- a/tests/testthat/test-rxode-issue-349.R +++ b/tests/testthat/test-rxode-issue-349.R @@ -1,14 +1,15 @@ rxTest({ test_that("Zero variances; RxODE#299", { + mod <- rxode2({ eff(0) <- 1 C2 <- centr / V2 C3 <- peri / V3 CL <- TCl * exp(eta.Cl) ## This is coded as a variable in the model - d / dt(depot) <- -KA * depot - d / dt(centr) <- KA * depot - CL * C2 - Q * C2 + Q * C3 - d / dt(peri) <- Q * C2 - Q * C3 - d / dt(eff) <- Kin - Kout * (1 - C2 / (EC50 + C2)) * eff + d/dt(depot) <- -KA * depot + d/dt(centr) <- KA * depot - CL * C2 - Q * C2 + Q * C3 + d/dt(peri) <- Q * C2 - Q * C3 + d/dt(eff) <- Kin - Kout * (1 - C2 / (EC50 + C2)) * eff e <- eff + eff.err cp <- centr * (1 + cp.err) })