diff --git a/NEWS.md b/NEWS.md index 3febe2b7e..dc9722cc0 100644 --- a/NEWS.md +++ b/NEWS.md @@ -2,8 +2,9 @@ - Rework the `tad()` and related functions so they use the same interface as compartments (this way they do not depend on the order - of compartments); See #815. Also allow dummy variables to ignore - state requirements (for parsing mu-referenced equations). + of compartments); See #815. For mu-referencing, Also allow dummy + variables to ignore state requirements (ie `podo(depot)` in a single + line will not error when parsing mu-referenced equations). - Add `getRxNpars` to api. This allows the development version of `babelmixr2` to better check what model is loaded and unload/reload @@ -12,6 +13,12 @@ - Add `rxUdfUiControl()` to rxode2 user function to get control information from something like `nlmixr2` +- Bug fix for tracking time after dose when dosing to 2 compartments + occur at the exact same time (#804, #819) + +- Change `transit()` model so that it uses `tad0()`, `podo0()` and + related functions for a bit more stable simulation and estimation + # rxode2 3.0.2 - Bug fix for `api`, the censoring function pointer has been updated diff --git a/R/d.R b/R/d.R index 8b62a6d1e..1ce974b25 100644 --- a/R/d.R +++ b/R/d.R @@ -546,6 +546,8 @@ .rxD$podo <- list(function(a) { return("0") }) +.rxD$podo0 <- .rxD$podo +.rxD$dose0 <- .rxD$dose .rxD$tlast <- list(function(a) { return("0") @@ -553,6 +555,10 @@ .rxD$tfirst <- list(function(a) { return("0") }) + +.rxD$tlast0 <- .rxD$tlast +.rxD$tfirst0 <- .rxD$tfirst + .rxD$first <- list(function(a) { return("0") }) diff --git a/R/parseFuns.R b/R/parseFuns.R index aded6523c..b94e31252 100644 --- a/R/parseFuns.R +++ b/R/parseFuns.R @@ -12,34 +12,35 @@ "fprec", "fround", "ftrunc", "transit", "gammaq", "gammapDer", "gammapInv", "gammapInva", "gammaqInv", "gammaqInva", "lowergamma", "uppergamma", "max", "min", "logit", "expit", "probit", "probitInv", -"tlast", "tfirst", "lag", "lead", "dose", "podo", "dabs", "dabs2", -"abs1", "dabs1", "erfinv", "abs0", "dosenum", "first", "last", -"diff", "is.nan", "is.na", "is.finite", "is.infinite", "llikPois", -"llikPoisDlambda", "llikBinom", "llikBinomDprob", "llikNbinom", -"llikNbinomDprob", "llikNbinomMu", "llikNbinomMuDmu", "llikBeta", -"llikBetaDshape1", "llikBetaDshape2", "llikT", "llikTDdf", "llikTDmean", -"llikTDsd", "llikChisq", "llikChisqDdf", "llikExp", "llikExpDrate", -"llikF", "llikFDdf1", "llikFDdf2", "llikGeom", "llikGeomDprob", -"llikUnif", "llikUnifDalpha", "llikUnifDbeta", "llikWeibull", -"llikWeibullDshape", "llikWeibullDscale", "llikGamma", "llikGammaDshape", -"llikGammaDrate", "llikCauchy", "llikCauchyDlocation", "llikCauchyDscale", -"llikNorm", "llikNormDmean", "llikNormDsd", "llikXPois", "llikXPoisDlambda", -"llikXBinom", "llikXBinomDprob", "llikXNbinomMu", "llikXNbinomMuDmu", -"llikXNbinom", "llikXNbinomDprob", "llikXBeta", "llikXBetaDshape1", -"llikXBetaDshape2", "llikXT", "llikXTDdf", "llikXTDmean", "llikXTDsd", -"llikXChisq", "llikXChisqDdf", "llikXExp", "llikXExpDrate", "llikXF", -"llikXFDdf1", "llikXFDdf2", "llikXGeom", "llikXGeomDprob", "llikXUnif", -"llikXUnifDalpha", "llikXUnifDbeta", "llikXWeibull", "llikXWeibullDshape", -"llikXWeibullDscale", "llikXGamma", "llikXGammaDshape", "llikXGammaDrate", -"llikXCauchy", "llikXCauchyDlocation", "llikXCauchyDscale", "llikXNorm", -"llikXNormDmean", "llikXNormDsd", "ReLU", "dReLU", "GELU", "dGELU", -"d2GELU", "d3GELU", "d4GELU", "ELU", "dELU", "d2ELU", "d2aELU", -"dELUa", "d2ELUa", "softplus", "dsoftplus", "d2softplus", "d3softplus", -"d4softplus", "SELU", "dSELU", "lReLU", "dlReLU", "PReLU", "dPReLU", -"d2PReLU", "dPReLUa", "dPReLUa1", "Swish", "dSwish", "linCmt", -"rnorm", "rxnorm", "rxbinom", "rbinom", "rxcauchy", "rcauchy", -"rchisq", "rxchisq", "rexp", "rxexp", "rbeta", "rxbeta", "rgeom", -"rxgeom", "rxpois", "rpois", "rxt", "rt") +"tlast", "tlast0", "tfirst", "tfirst0", "lag", "lead", "dose", +"podo", "dose0", "podo0", "dabs", "dabs2", "abs1", "dabs1", "erfinv", +"abs0", "dosenum", "first", "last", "diff", "is.nan", "is.na", +"is.finite", "is.infinite", "llikPois", "llikPoisDlambda", "llikBinom", +"llikBinomDprob", "llikNbinom", "llikNbinomDprob", "llikNbinomMu", +"llikNbinomMuDmu", "llikBeta", "llikBetaDshape1", "llikBetaDshape2", +"llikT", "llikTDdf", "llikTDmean", "llikTDsd", "llikChisq", "llikChisqDdf", +"llikExp", "llikExpDrate", "llikF", "llikFDdf1", "llikFDdf2", +"llikGeom", "llikGeomDprob", "llikUnif", "llikUnifDalpha", "llikUnifDbeta", +"llikWeibull", "llikWeibullDshape", "llikWeibullDscale", "llikGamma", +"llikGammaDshape", "llikGammaDrate", "llikCauchy", "llikCauchyDlocation", +"llikCauchyDscale", "llikNorm", "llikNormDmean", "llikNormDsd", +"llikXPois", "llikXPoisDlambda", "llikXBinom", "llikXBinomDprob", +"llikXNbinomMu", "llikXNbinomMuDmu", "llikXNbinom", "llikXNbinomDprob", +"llikXBeta", "llikXBetaDshape1", "llikXBetaDshape2", "llikXT", +"llikXTDdf", "llikXTDmean", "llikXTDsd", "llikXChisq", "llikXChisqDdf", +"llikXExp", "llikXExpDrate", "llikXF", "llikXFDdf1", "llikXFDdf2", +"llikXGeom", "llikXGeomDprob", "llikXUnif", "llikXUnifDalpha", +"llikXUnifDbeta", "llikXWeibull", "llikXWeibullDshape", "llikXWeibullDscale", +"llikXGamma", "llikXGammaDshape", "llikXGammaDrate", "llikXCauchy", +"llikXCauchyDlocation", "llikXCauchyDscale", "llikXNorm", "llikXNormDmean", +"llikXNormDsd", "ReLU", "dReLU", "GELU", "dGELU", "d2GELU", "d3GELU", +"d4GELU", "ELU", "dELU", "d2ELU", "d2aELU", "dELUa", "d2ELUa", +"softplus", "dsoftplus", "d2softplus", "d3softplus", "d4softplus", +"SELU", "dSELU", "lReLU", "dlReLU", "PReLU", "dPReLU", "d2PReLU", +"dPReLUa", "dPReLUa1", "Swish", "dSwish", "linCmt", "rnorm", +"rxnorm", "rxbinom", "rbinom", "rxcauchy", "rcauchy", "rchisq", +"rxchisq", "rexp", "rxexp", "rbeta", "rxbeta", "rgeom", "rxgeom", +"rxpois", "rpois", "rxt", "rt") .parseEnv$.parseNum <- c(lgamma = 1, abs = 1, acos = 1, acosh = 1, asin = 1, asinh = 1, atan = 1, atan2 = 2, atanh = 1, beta = 2, cos = 1, cosh = 1, erf = 1, erfc = 1, exp = 1, gamma = 1, linCmtA = 20, linCmtC = 20, @@ -52,28 +53,29 @@ fprec = 2, fround = 2, ftrunc = 2, transit = NA, gammaq = 2, gammapDer = 2, gammapInv = 2, gammapInva = 2, gammaqInv = 2, gammaqInva = 2, lowergamma = 2, uppergamma = 2, max = NA, min = NA, logit = NA, expit = NA, probit = NA, probitInv = NA, tlast = NA, -tfirst = NA, lag = NA, lead = NA, dose = NA, podo = NA, dabs = 1, -dabs2 = 1, abs1 = 1, dabs1 = 1, erfinv = 1, abs0 = 1, dosenum = 0, -first = 1, last = 1, diff = 1, is.nan = 1, is.na = 1, is.finite = 1, -is.infinite = 1, llikPois = 2, llikPoisDlambda = 2, llikBinom = 3, -llikBinomDprob = 3, llikNbinom = 3, llikNbinomDprob = 3, llikNbinomMu = 3, -llikNbinomMuDmu = 3, llikBeta = 3, llikBetaDshape1 = 3, llikBetaDshape2 = 3, -llikT = 4, llikTDdf = 4, llikTDmean = 4, llikTDsd = 4, llikChisq = 2, -llikChisqDdf = 2, llikExp = 2, llikExpDrate = 2, llikF = 3, llikFDdf1 = 3, -llikFDdf2 = 3, llikGeom = 2, llikGeomDprob = 2, llikUnif = 3, -llikUnifDalpha = 3, llikUnifDbeta = 3, llikWeibull = 3, llikWeibullDshape = 3, -llikWeibullDscale = 3, llikGamma = 3, llikGammaDshape = 3, llikGammaDrate = 3, -llikCauchy = 3, llikCauchyDlocation = 3, llikCauchyDscale = 3, -llikNorm = 3, llikNormDmean = 3, llikNormDsd = 3, llikXPois = 3, -llikXPoisDlambda = 3, llikXBinom = 4, llikXBinomDprob = 4, llikXNbinomMu = 4, -llikXNbinomMuDmu = 4, llikXNbinom = 4, llikXNbinomDprob = 4, -llikXBeta = 4, llikXBetaDshape1 = 4, llikXBetaDshape2 = 4, llikXT = 5, -llikXTDdf = 5, llikXTDmean = 5, llikXTDsd = 5, llikXChisq = 3, -llikXChisqDdf = 3, llikXExp = 3, llikXExpDrate = 3, llikXF = 4, -llikXFDdf1 = 4, llikXFDdf2 = 4, llikXGeom = 3, llikXGeomDprob = 3, -llikXUnif = 4, llikXUnifDalpha = 4, llikXUnifDbeta = 4, llikXWeibull = 4, -llikXWeibullDshape = 4, llikXWeibullDscale = 4, llikXGamma = 4, -llikXGammaDshape = 4, llikXGammaDrate = 4, llikXCauchy = 4, llikXCauchyDlocation = 4, +tlast0 = NA, tfirst = NA, tfirst0 = NA, lag = NA, lead = NA, +dose = NA, podo = NA, dose0 = NA, podo0 = NA, dabs = 1, dabs2 = 1, +abs1 = 1, dabs1 = 1, erfinv = 1, abs0 = 1, dosenum = 0, first = 1, +last = 1, diff = 1, is.nan = 1, is.na = 1, is.finite = 1, is.infinite = 1, +llikPois = 2, llikPoisDlambda = 2, llikBinom = 3, llikBinomDprob = 3, +llikNbinom = 3, llikNbinomDprob = 3, llikNbinomMu = 3, llikNbinomMuDmu = 3, +llikBeta = 3, llikBetaDshape1 = 3, llikBetaDshape2 = 3, llikT = 4, +llikTDdf = 4, llikTDmean = 4, llikTDsd = 4, llikChisq = 2, llikChisqDdf = 2, +llikExp = 2, llikExpDrate = 2, llikF = 3, llikFDdf1 = 3, llikFDdf2 = 3, +llikGeom = 2, llikGeomDprob = 2, llikUnif = 3, llikUnifDalpha = 3, +llikUnifDbeta = 3, llikWeibull = 3, llikWeibullDshape = 3, llikWeibullDscale = 3, +llikGamma = 3, llikGammaDshape = 3, llikGammaDrate = 3, llikCauchy = 3, +llikCauchyDlocation = 3, llikCauchyDscale = 3, llikNorm = 3, +llikNormDmean = 3, llikNormDsd = 3, llikXPois = 3, llikXPoisDlambda = 3, +llikXBinom = 4, llikXBinomDprob = 4, llikXNbinomMu = 4, llikXNbinomMuDmu = 4, +llikXNbinom = 4, llikXNbinomDprob = 4, llikXBeta = 4, llikXBetaDshape1 = 4, +llikXBetaDshape2 = 4, llikXT = 5, llikXTDdf = 5, llikXTDmean = 5, +llikXTDsd = 5, llikXChisq = 3, llikXChisqDdf = 3, llikXExp = 3, +llikXExpDrate = 3, llikXF = 4, llikXFDdf1 = 4, llikXFDdf2 = 4, +llikXGeom = 3, llikXGeomDprob = 3, llikXUnif = 4, llikXUnifDalpha = 4, +llikXUnifDbeta = 4, llikXWeibull = 4, llikXWeibullDshape = 4, +llikXWeibullDscale = 4, llikXGamma = 4, llikXGammaDshape = 4, +llikXGammaDrate = 4, llikXCauchy = 4, llikXCauchyDlocation = 4, llikXCauchyDscale = 4, llikXNorm = 4, llikXNormDmean = 4, llikXNormDsd = 4, ReLU = 1, dReLU = 1, GELU = 1, dGELU = 1, d2GELU = 1, d3GELU = 1, d4GELU = 1, ELU = 2, dELU = 2, d2ELU = 2, d2aELU = 2, dELUa = 2, diff --git a/R/symengine.R b/R/symengine.R index 07bf6babc..e8331aec2 100644 --- a/R/symengine.R +++ b/R/symengine.R @@ -142,11 +142,15 @@ regIfOrElse <- rex::rex(or(regIf, regElse)) "probit" = NA, "probitInv" = NA, "tlast" = NA, + "tlast0" = NA, "tfirst" = NA, + "tfirst0" = NA, "lag" = NA, "lead" = NA, "dose" =NA, "podo" =NA, + "dose0" =NA, + "podo0" =NA, "dabs" = 1, "dabs2" = 1, "abs1" = 1, @@ -662,7 +666,7 @@ rxToSE <- function(x, envir = NULL, progress = FALSE, } .rxToSEDualVarFunction <- c("tlast", "tlast0", "tad", "tad0", "tafd", "tafd0", - "dose", "podo") + "dose", "podo", "dose0", "podo0") #' Change rxode2 syntax to symengine syntax for symbols and numbers #' @@ -1090,6 +1094,9 @@ rxToSE <- function(x, envir = NULL, progress = FALSE, if (identical(x[[1]], quote(`podo`))) { return(paste0("podo(", .rxLastAssignedDdt, ")")) } + if (identical(x[[1]], quote(`podo0`))) { + return(paste0("podo0(", .rxLastAssignedDdt, ")")) + } } else if (.len == 2L) { if (length(x[[2]]) != 1) { stop(as.character(x[[1]]), "() must be used with a state", call. = FALSE) @@ -1284,11 +1291,11 @@ rxToSE <- function(x, envir = NULL, progress = FALSE, .bio <- .rxToSE(x[[4]], envir = envir) if (isEnv) envir$..curCall <- .lastCall return(paste0( - "exp(log((", .bio, ")*(podo(", .rxLastAssignedDdt, ")))+log(", + "exp(log((", .bio, ")*(podo0(", .rxLastAssignedDdt, ")))+log(", .n, " + 1)-log(", .mtt, ")+(", .n, ")*((log(", .n, "+1)-log(", .mtt, - "))+log(t-tlast(", .rxLastAssignedDdt, ")))-((", .n, "+1)/(", .mtt, - "))*(t-tlast(", .rxLastAssignedDdt, "))-lgamma(1+", .n, "))" + "))+log(t-tlast0(", .rxLastAssignedDdt, ")))-((", .n, "+1)/(", .mtt, + "))*(t-tlast0(", .rxLastAssignedDdt, "))-lgamma(1+", .n, "))" )) } else if (length(x) == 3) { if (isEnv) { @@ -1298,7 +1305,7 @@ rxToSE <- function(x, envir = NULL, progress = FALSE, .n <- .rxToSE(x[[2]], envir = envir) .mtt <- .rxToSE(x[[3]], envir = envir) if (isEnv) envir$..curCall <- .lastCall - return(paste0("exp(log(podo(", .rxLastAssignedDdt, "))+(log(", .n, "+1)-log(", .mtt, "))+(", .n, ")*((log(", .n, "+1)-log(", .mtt, "))+ log(t-tlast(", .rxLastAssignedDdt, ")))-((", .n, " + 1)/(", .mtt, "))*(t-tlast(",.rxLastAssignedDdt, "))-lgamma(1+", .n, "))")) + return(paste0("exp(log(podo0(", .rxLastAssignedDdt, "))+(log(", .n, "+1)-log(", .mtt, "))+(", .n, ")*((log(", .n, "+1)-log(", .mtt, "))+ log(t-tlast0(", .rxLastAssignedDdt, ")))-((", .n, " + 1)/(", .mtt, "))*(t-tlast0(",.rxLastAssignedDdt, "))-lgamma(1+", .n, "))")) } else { stop("'transit' can only take 2-3 arguments", call. = FALSE) } @@ -1351,7 +1358,7 @@ rxToSE <- function(x, envir = NULL, progress = FALSE, } else if (identical(x[[1]], quote(`tad`))) { return(.rxToSETad(x, envir = envir, progress = progress, isEnv=isEnv)) } else if (identical(x[[1]], quote(`tad0`))) { - return(.rxToSETad(x, envir = envir, progress = progress, isEnv=isEnv)) + return(.rxToSETad0(x, envir = envir, progress = progress, isEnv=isEnv)) } else if (identical(x[[1]], quote(`lag`)) || identical(x[[1]], quote(`lead`))) { return(.rxToSELagOrLead(x, envir = envir, progress = progress, isEnv=isEnv)) @@ -1361,10 +1368,12 @@ rxToSE <- function(x, envir = NULL, progress = FALSE, return(.rxToSETlastOrTafd0(x, envir = envir, progress = progress, isEnv=isEnv)) } else if (identical(x[[1]], quote(`tlast`)) || identical(x[[1]], quote(`tfirst`)) || - identical(x[[1]], quote(`last0`)) || - identical(x[[1]], quote(`first0`)) || + identical(x[[1]], quote(`tlast0`)) || + identical(x[[1]], quote(`tfirst0`)) || identical(x[[1]], quote(`dose`)) || - identical(x[[1]], quote(`podo`))) { + identical(x[[1]], quote(`podo`)) || + identical(x[[1]], quote(`dose0`)) || + identical(x[[1]], quote(`podo0`))) { return(.rxToSETlastOrTfirst(x, envir = envir, progress = progress, isEnv=isEnv)) } else if (identical(x[[1]], quote(`psigamma`))) { return(.rxToSEPsigamma(x, envir = envir, progress = progress, isEnv=isEnv)) @@ -2375,7 +2384,8 @@ rxFromSE <- function(x, unknownDerivatives = c("forward", "central", "error"), ")" ) return(.ret) - } else if (any(paste(.ret0[[1]]) == c("tlast", "tfirst", "dose", "podo"))) { + } else if (any(paste(.ret0[[1]]) == c("tlast", "tfirst", "dose", "podo", + "tlast0", "first0", "dose0", "podo0"))) { if (length(.ret0) == 1L) { return(paste0(.ret0[[1]], "()")) } else if (length(.ret0) == 2L) { diff --git a/data/rxReservedKeywords.rda b/data/rxReservedKeywords.rda index 1ec5bb6b3..e05da3bbc 100644 Binary files a/data/rxReservedKeywords.rda and b/data/rxReservedKeywords.rda differ diff --git a/data/rxResidualError.rda b/data/rxResidualError.rda index 5a104fee5..1f9f7e2d7 100644 Binary files a/data/rxResidualError.rda and b/data/rxResidualError.rda differ diff --git a/data/rxSyntaxFunctions.rda b/data/rxSyntaxFunctions.rda index c88115476..f96214a93 100644 Binary files a/data/rxSyntaxFunctions.rda and b/data/rxSyntaxFunctions.rda differ diff --git a/inst/include/rxode2_model_shared.h b/inst/include/rxode2_model_shared.h index 073000d55..0f3d1e3bd 100644 --- a/inst/include/rxode2_model_shared.h +++ b/inst/include/rxode2_model_shared.h @@ -219,9 +219,14 @@ static inline double Rx_pow_di_(double a, double b, rx_solve *rx) { #define _logitInv1(x) expit(x, 0.0, 1.0) #define _logitInv2(x, y) expit(x, y, 1.0) #define _podo0() (_solveData->subjects[_cSub].curDose) +#define _podo00() (ISNA(_solveData->subjects[_cSub].curDose) ? 0 : _solveData->subjects[_cSub].curDose) #define _podo1(x) (_solveData->subjects[_cSub].curDoseS[x]) +#define _podo01(x) (ISNA(_solveData->subjects[_cSub].curDoseS[x]) ? 0 : _solveData->subjects[_cSub].curDoseS[x]) + #define _dose0() (_solveData->subjects[_cSub].curDose) #define _dose1(x) (_solveData->subjects[_cSub].curDoseS[x]) +#define _dose00() (ISNA(_solveData->subjects[_cSub].curDose) ? 0 : _solveData->subjects[_cSub].curDose) +#define _dose01(x) (ISNA(_solveData->subjects[_cSub].curDoseS[x]) ? 0 : _solveData->subjects[_cSub].curDoseS[x]) #define _tad0() (t-_solveData->subjects[_cSub].tlast) #define _tad1(x) (t-_solveData->subjects[_cSub].tlastS[x]) #define _tad00() (ISNA(_solveData->subjects[_cSub].tlast)? 0 : (t- _solveData->subjects[_cSub].tlast)) diff --git a/inst/include/rxode2parseHandleEvid.h b/inst/include/rxode2parseHandleEvid.h index 165b402bb..c99ae4e18 100644 --- a/inst/include/rxode2parseHandleEvid.h +++ b/inst/include/rxode2parseHandleEvid.h @@ -248,9 +248,10 @@ static inline void handleTlastInline(double *time, rx_solving_options_ind *ind) } else { evid = getEvid(ind, ind->ix[ind->idx]); } - if (op->neq + op->extraCmt != 0 && ind->tlast != _time && + if (op->neq + op->extraCmt != 0 && isDose(evid) && - ind->cmt < op->neq + op->extraCmt) { + ind->cmt < op->neq + op->extraCmt && + ind->tlastS[ind->cmt] != _time) { double curDose = getDoseIndex(ind, ind->idx), tinf = NA_REAL; if (handleTlastInlineUpateDosingInformation(ind, &curDose, &tinf) == 0) return; ind->dosenum++; diff --git a/src/genModelVars.h b/src/genModelVars.h index 60efbbc4f..26f76e958 100644 --- a/src/genModelVars.h +++ b/src/genModelVars.h @@ -250,6 +250,12 @@ static inline int sortStateVectorsErrHandle(int prop, int pass, int i) { if ((prop & propDose) != 0) { sAppend(&sbt, "'dose(%s)', ", buf); } + if ((prop & propPodo0) != 0) { + sAppend(&sbt, "'podo0(%s)', ", buf); + } + if ((prop & propDose0) != 0) { + sAppend(&sbt, "'dose0(%s)', ", buf); + } // Take off trailing "', sbt.o -= 2; sbt.s[sbt.o] = 0; diff --git a/src/parseFuns.h b/src/parseFuns.h index c3b36afcd..fc136758d 100644 --- a/src/parseFuns.h +++ b/src/parseFuns.h @@ -95,6 +95,9 @@ typedef struct transFunctions { int isDose; int isPodo; + int isDose0; + int isPodo0; + int isTfirst; int isTfirst0; @@ -140,6 +143,10 @@ static inline void transFunctionsIni(transFunctions *tf) { tf->isDose = 0; tf->isPodo = 0; + tf->isDose0 = 0; + tf->isPodo0 = 0; + + tf->isTfirst = 0; tf->isTfirst0 = 0; diff --git a/src/parseFunsDosing.h b/src/parseFunsDosing.h index 2ad5cfbe8..bb45fb6ff 100644 --- a/src/parseFunsDosing.h +++ b/src/parseFunsDosing.h @@ -36,7 +36,10 @@ static inline int isFunctionTadType(transFunctions *tf) { (tf->isTfirst0 = !strcmp("tfirst0", tf->v)) || (tf->isDose = !strcmp("dose", tf->v)) || - (tf->isPodo = !strcmp("podo", tf->v)); + (tf->isPodo = !strcmp("podo", tf->v)) || + + (tf->isPodo0 = !strcmp("podo0", tf->v)) || + (tf->isDose0 = !strcmp("dose0", tf->v)); } static inline int handleFunctionTadEmptyCcode(transFunctions *tf,char *v2) { if (allSpaces(v2)){ @@ -101,6 +104,10 @@ static inline int handleFunctionTadSingleStateCcode(transFunctions *tf,char *v2) tb.dprop[tb.id] += propDose; } else if (tf->isPodo && (tb.dprop[tb.id] & propPodo) == 0) { tb.dprop[tb.id] += propPodo; + } else if (tf->isPodo0 && (tb.dprop[tb.id] & propPodo0) == 0) { + tb.dprop[tb.id] += propPodo0; + } else if (tf->isDose0 && (tb.dprop[tb.id] & propDose0) == 0) { + tb.dprop[tb.id] += propDose0; } return 1; } diff --git a/src/tran.h b/src/tran.h index 38819f336..07585345b 100644 --- a/src/tran.h +++ b/src/tran.h @@ -389,4 +389,7 @@ char *getLine (char *src, int line, int *lloc); #define propPodo 8192 #define propDose 16384 +#define propPodo0 32768 +#define propDose0 65536 + #endif // __TRAN_H__ diff --git a/tests/testthat/test-dsl.R b/tests/testthat/test-dsl.R index 4ad417168..51bb9c536 100644 --- a/tests/testthat/test-dsl.R +++ b/tests/testthat/test-dsl.R @@ -210,21 +210,22 @@ rxTest({ expect_equal( rxToSE(transit(n, mtt, bio)), - "exp(log((bio)*(podo()))+log(n + 1)-log(mtt)+(n)*((log(n+1)-log(mtt))+log(t-tlast()))-((n+1)/(mtt))*(t-tlast())-lgamma(1+n))") + "exp(log((bio)*(podo0()))+log(n + 1)-log(mtt)+(n)*((log(n+1)-log(mtt))+log(t-tlast0()))-((n+1)/(mtt))*(t-tlast0())-lgamma(1+n))") expect_equal( rxToSE(transit(n, mtt)), - "exp(log(podo())+(log(n+1)-log(mtt))+(n)*((log(n+1)-log(mtt))+ log(t-tlast()))-((n + 1)/(mtt))*(t-tlast())-lgamma(1+n))") + "exp(log(podo0())+(log(n+1)-log(mtt))+(n)*((log(n+1)-log(mtt))+ log(t-tlast0()))-((n + 1)/(mtt))*(t-tlast0())-lgamma(1+n))") tmp <- rxode("d/dt(depot) <- transit(n, mtt, bio)-ka*depot\nd/dt(center)=ka*depot-kel*center") + tmp2 <- rxS(tmp) tmp3 <- tmp2$rx__d_dt_depot__ - expect_equal(rxFromSE(tmp3), "-ka*depot+exp(n*(-log(mtt)+log1p(n)+log(t-tlast(depot)))-(t-tlast(depot))*(1+n)/mtt-log(mtt)+log(bio*podo(depot))+log1p(n)-lgamma1p(n))") + expect_equal(rxFromSE(tmp3), "-ka*depot+exp(n*(-log(mtt)+log1p(n)+log(t-tlast0(depot)))-(1+n)*(t-tlast0(depot))/mtt-log(mtt)+log(bio*podo0(depot))+log1p(n)-lgamma1p(n))") tmp <- rxode("d/dt(depot) <- transit(n, mtt) - ka*depot\nd/dt(center)=ka*depot-kel*center") tmp2 <- rxS(tmp) tmp3 <- tmp2$rx__d_dt_depot__ - expect_equal(rxFromSE(tmp3), "-ka*depot+exp(n*(-log(mtt)+log1p(n)+log(t-tlast(depot)))-(t-tlast(depot))*(1+n)/mtt-log(mtt)+log1p(n)+log(podo(depot))-lgamma1p(n))") + expect_equal(rxFromSE(tmp3), "-ka*depot+exp(n*(-log(mtt)+log1p(n)+log(t-tlast0(depot)))-(1+n)*(t-tlast0(depot))/mtt-log(mtt)+log1p(n)+log(podo0(depot))-lgamma1p(n))") }) @@ -519,6 +520,13 @@ rxTest({ expect_error(rxToSE("tad(matt+f)")) }) + test_that("tad0()", { + expect_equal(rxToSE("tad0()"), "(t-tlast0())") + expect_equal(rxToSE("tad0(matt)"), "(t-tlast0(matt))") + expect_error(rxToSE("tad0(matt,f)")) + expect_error(rxToSE("tad0(matt+f)")) + }) + test_that("tafd()", { expect_equal(rxToSE("tafd()"), "(t-tfirst())") expect_equal(rxToSE("tafd(matt)"), "(t-tfirst(matt))") diff --git a/tests/testthat/test-transit.R b/tests/testthat/test-transit.R index dafca22ad..7f132a77c 100644 --- a/tests/testthat/test-transit.R +++ b/tests/testthat/test-transit.R @@ -204,4 +204,36 @@ transit = matt + fun ") expect_s3_class(mod, "rxode2") }) + + test_that("transit compartment works well with dual absorption (#804, #819)", { + + mod <- function() { + ini({ + ## Table 3 from Savic 2007 + cl <- 17.2 # (L/hr) + vc <- 45.1 # L + ka <- 0.38 # 1/hr + mtt <- 1.37 # hr + f2 <-0.5 # Fraction of 1st Order portion + n <- 20.101 + }) + model({ + k <- cl/vc + bio <- 1-f2 + d/dt(depot1) <- transit(n,mtt,bio)-ka*depot1 + d/dt(depot2) <- -ka*depot2 + f(depot2) <-f2 + d/dt(cen) <- ka*depot1 + ka*depot2-k*cen + }) + } + + ev1 <- et(0, 7, length.out=200) %>% + et(amt=20, cmt='depot1', evid=7) %>% + et(amt=20, cmt='depot2', evid=1) + + case1 <- rxSolve(mod, ev1) + + expect_true(max(case1$depot1) > 7.7) + + }) })