regIf <- rex::rex(start, any_spaces, "if", any_spaces, "(", capture(anything), ")", any_spaces, "{", any_spaces, end)
regElse <- rex::rex(start, any_spaces, "else", any_spaces, "{", any_spaces, end)
regEnd <- rex::rex(start, any_spaces, "}", any_spaces, end)
regIfOrElse <- rex::rex(or(regIf, regElse))
## first.arg second.arg and type of function
## RxODE->symengine
.rxSEsingle <- list(
"gammafn" = c("gamma(", ")", "gamma"),
"lgammafn" = c("lgamma(", ")", "lgamma"),
"lgamma" = c("lgamma(", ")", "lgamma"),
"loggamma" = c("lgamma(", ")", "lgamma"),
"digamma" = c("polygamma(0,", ")", "psigamma"),
"trigamma" = c("polygamma(1,", ")", "psigamma"),
"tetragamma" = c("polygamma(2,", ")", "psigamma"),
"pentagamma" = c("polygamma(3,", ")", "psigamma"),
"cospi" = c("cos(pi*(", "))", "cos"),
"sinpi" = c("sin(pi*(", "))", "sin"),
"tanpi" = c("tan(pi*(", "))", "tan"),
"log1p" = c("log(1+", ")", "log"),
"expm1" = c("(exp(", ")-1)", "exp"),
"factorial" = c("gamma(", "+1)", "gamma"),
"lfactorial" = c("lgamma(", "+1)", "lgamma"),
"lgamma1p" = c("lgamma(", "+1)", "lgamma"),
"expm1" = c("(exp(", ")-1)", "exp"),
"log10" = c("log(", ")/log(10)", "log"),
"log2" = c("log(", ")/log(2)", "log"),
"log1pexp" = c("log(1+exp(", "))", "log1pexp"),
"!" = c("rxNot(", ")", ""),
"phi" = c("0.5*(1+erf((", ")/sqrt(2)))"),
"pnorm" = c("0.5*(1+erf((", ")/sqrt(2)))"),
"normcdf" = c("0.5*(1+erf((", ")/sqrt(2)))"),
"qnorm" = c("sqrt(2)*erfinv(2*(", ")-1)"),
"fabs" = c("abs0(", ")")
)
.SEsingle <- list(
"rxNot" = c("(!(", "))"),
"loggamma" = c("lgamma(", ")")
)
.rxSEdouble <- list(
"pow" = c("(", ")^(", ")"),
"R_pow" = c("(", ")^(", ")"),
"R_pow_di" = c("(", ")^(", ")"),
"Rx_pow_di" = c("(", ")^(", ")"),
"Rx_pow" = c("(", ")^(", ")"),
"lbeta" = c("log(beta(", ",", "))"),
"==" = c("rxEq(", ",", ")"),
"!=" = c("rxNeq(", ",", ")"),
">=" = c("rxGeq(", ",", ")"),
"<=" = c("rxLeq(", ",", ")"),
"<" = c("rxLt(", ",", ")"),
">" = c("rxGt(", ",", ")"),
"&&" = c("rxAnd(", ",", ")"),
"||" = c("rxOr(", ",", ")"),
"&" = c("rxAnd(", ",", ")"),
"|" = c("rxOr(", ",", ")")
)
.SEdouble <- list(
"lbeta" = c("lbeta(", ",", ")"),
"rxEq" = c("(", "==", ")"),
"rxNeq" = c("(", "!=", ")"),
"rxGeq" = c("(", ">=", ")"),
"rxLeq" = c("(", "<=", ")"),
"rxGt" = c("(", ">", ")"),
"rxLt" = c("(", "<", ")"),
"rxAnd" = c("(", "&&", ")"),
"rxOr" = c("(", "||", ")")
)
## atan2
.rxSEeq <- 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,
"linCmtB" = 21,
"log" = 1,
"polygamma" = 2,
"rxTBS" = 5,
"rxTBSi" = 5,
"rxTBSd" = 5,
"rxTBSd2" = 5,
"sin" = 1,
"sinh" = 1,
"sqrt" = 1,
"tan" = 1,
"tanh" = 1,
"gammap" = 2,
## C's math.h library
"floor" = 1,
"round" = 1,
"ceil" = 1,
"trunc" = 1,
## Special R functions
"bessel_i" = 3,
"bessel_j" = 2,
"bessel_k" = 3,
"bessel_y" = 2,
"logspace_add" = 2,
"logspace_sub" = 2,
"fmax2" = 2,
"fmin2" = 2,
"sign" = 1,
"fsign" = 2,
"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,
"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
)
.rxOnly <- c(
## Now random number generators
"rnorm" = NA,
"rxnorm" = NA,
"rxbinom" = 2,
"rbinom" = 2,
"rxcauchy" = NA,
"rcauchy" = NA,
"rchisq" = 1,
"rxchisq" = 1,
"rexp" = 1,
"rxexp" = 1,
"rbeta" = 2,
"rxbeta" = 2,
"rgeom" = 1,
"rxgeom" = 1,
"rxpois" = 1,
"rpois" = 1,
"rxt" = 1,
"rt" = 1
)
.rxSEeqUsr <- NULL
.rxCcode <- NULL
.symengineFs <- new.env(parent = emptyenv())
.extraCnow <- ""
.extraC <- function(extraC = NULL) {
if (!is.null(extraC)) {
if (file.exists(extraC)) {
.ret <- sprintf("#include \"%s\"\n", extraC)
} else {
.ret <- paste(extraC, collapse = "\n")
}
} else {
.ret <- ""
}
if (length(.rxCcode) > 0L) {
.ret <- sprintf("%s\n%s\n", .ret, paste(.rxCcode, collapse = "\n"))
}
assignInMyNamespace(".extraCnow", .ret)
return(invisible())
}
#' Add user function to RxODE
#'
#' This adds a user function to RxODE that can be called. If needed,
#' these functions can be differentiated by numerical differences or
#' by adding the derivatives to RxODE's internal derivative table
#' with [rxD()]
#'
#' @param name This gives the name of the user function
#' @param args This gives the arguments of the user function
#' @param cCode This is the C-code for the new function
#' @return nothing
#' @author Matthew L. Fidler
#' @examples
#' \donttest{
#' ## Right now RxODE is not aware of the function f
#' ## Therefore it cannot translate it to symengine or
#' ## Compile a model with it.
#'
#' try(RxODE("a=fun(a,b,c)"))
#'
#' ## Note for this approach to work, it cannot interfere with C
#' ## function names or reserved RxODE specical terms. Therefore
#' ## f(x) would not work since f is an alias for bioaviability.
#'
#' fun <- "
#' double fun(double a, double b, double c) {
#' return a*a+b*a+c;
#' }
#' " ## C-code for function
#'
#' rxFun("fun", c("a", "b", "c"), fun) ## Added function
#'
#' ## Now RxODE knows how to translate this function to symengine
#'
#' rxToSE("fun(a,b,c)")
#'
#' ## And will take a central difference when calculating derivatives
#'
#' rxFromSE("Derivative(fun(a,b,c),a)")
#'
#' ## Of course, you could specify the derivative table manually
#' rxD("fun", list(
#' function(a, b, c) {
#' paste0("2*", a, "+", b)
#' },
#' function(a, b, c) {
#' return(a)
#' },
#' function(a, b, c) {
#' return("0.0")
#' }
#' ))
#'
#' rxFromSE("Derivative(fun(a,b,c),a)")
#'
#' # You can also remove the functions by `rxRmFun`
#'
#' rxRmFun("fun")
#' }
#' @export
rxFun <- function(name, args, cCode) {
if (!is.character(name) || length(name) != 1L) {
stop("name argument must be a length-one character vector", call. = FALSE)
}
if (missing(cCode)) stop("a new function requires a C function so it can be used in RxODE", call. = FALSE)
if (any(name == names(.rxSEeqUsr))) {
stop("already defined user function '", name, "', remove it fist ('rxRmFun')",
call. = FALSE
)
}
suppressWarnings(rxRmFun(name))
assignInMyNamespace(".rxSEeqUsr", c(.rxSEeqUsr, setNames(length(args), name)))
assignInMyNamespace(".rxCcode", c(.rxCcode, setNames(cCode, name)))
assign(name, symengine::Function(name), envir = .symengineFs)
return(invisible())
}
#' @rdname rxFun
#' @export
rxRmFun <- function(name) {
if (!is.character(name) || length(name) != 1L) {
stop("name argument must be a length-one character vector",
call. = FALSE
)
}
if (!any(name == names(.rxSEeqUsr))) {
warning("no user function '", name, "' to remove", call. = FALSE)
}
.w <- which(name == names(.rxSEeqUsr))
if (length(.w) == 1L) assignInMyNamespace(".rxSEeqUsr", .rxSEeqUsr[-.w])
.w <- which(name == names(.rxCcode))
if (length(.w) == 1L) assignInMyNamespace(".rxCcode", .rxCcode[-.w])
if (exists(name, envir = .rxD)) rm(list = name, envir = .rxD)
if (exists(name, envir = .symengineFs)) rm(list = name, envir = .symengineFs)
return(invisible())
}
.SE1p <- c(
"loggamma" = "lgamma1p",
"lgamma" = "lgamma1p",
"log" = "log1p"
)
.SE1m <- c(
"cos" = "cospi",
"sin" = "sinpi",
"tan" = "tanpi"
)
## "rxTBS", "rxTBSd"
.rxSEcnt <- c(
"M_E" = "E",
"M_PI" = "pi",
"M_PI_2" = "pi/2",
"M_PI_4" = "pi/4",
"M_1_PI" = "1/pi",
"M_2_PI" = "2/pi",
"M_2PI" = "2*pi",
"M_SQRT_PI" = "sqrt(pi)",
"M_2_SQRTPI" = "2/sqrt(pi)",
"M_1_SQRT_2PI" = "1/sqrt(2*pi)",
"M_SQRT_2" = "sqrt(2)",
"M_SQRT_3" = "sqrt(3)",
"M_SQRT_32" = "sqrt(32)",
"M_SQRT_2dPI" = "sqrt(2/pi)",
"M_LN_SQRT_PI" = "log(sqrt(pi))",
"M_LN_SQRT_2PI" = "log(sqrt(2*pi))",
"M_LN_SQRT_PId2" = "log(sqrt(pi/2))",
"M_SQRT2" = "sqrt(2)",
"M_SQRT3" = "sqrt(3)",
"M_SQRT32" = "sqrt(32)",
"M_LOG10_2" = "log(2)/log(10)",
"M_LOG2E" = "1/log(2)",
"M_LOG10E" = "1/log(10)",
"M_LN2" = "log(2)",
"M_LN10" = "log(10)"
)
## "rxTBS", "rxTBSd"
.rxSEreserved <- list(
"e" = 2.718281828459045090796,
"E" = 2.718281828459045090796,
"EulerGamma" = 0.57721566490153286060651209008240243104215933593992,
"Catalan" = 0.915965594177219015054603514932384110774,
"GoldenRatio" = 2.118033988749894902526,
"I" = 1i
)
## diff(rxTBS(a,lambda,yj),a)
.rxD <- new.env(parent = emptyenv())
## This environment is a derivative table;
## For example:
## Derivative(f(a,b,c), a) = fa()
## Derivative(f(a,b,c), b) = fb()
## Derivative(f(a,b,c), c) = fc()
## Then
##
## .rxD$f <- list(fa(a,b,c), fb(a,b,c), fc(a,b,c))
##
## fa translates the arguments to the derivative with respect to a
## fb translates the arguments to the derivative with respect to b
##
## If any of the list is NULL then RxODE won't know how to take a
## derivative with respect to the argument.
##
## If the list is shorter than the length of the arguments then the
## argument then the derivative of arguments that are not specified
## cannot be taken.
.rxD$atan2 <- list(
function(y, x) {
return(paste0("(", x, ")/((", x, ")^2+(", y, ")^2)"))
},
function(y, x) {
return(paste0("-(", y, ")/((", x, ")^2+(", y, ")^2)"))
}
)
.rxD$erfinv <- list(
function(x) {
## http://specialfunctionswiki.org/index.php/Derivative_of_inverse_error_function
return(paste0("sqrt(pi)/2*exp((erfinv(", x, "))^2)"))
}
)
.rxD$abs0 <- list(function(x) {
return(paste0("dabs(", x, ")"))
})
.rxD$abs <- list(function(x) {
return(paste0("dabs(", x, ")"))
})
.rxD$abs1 <- list(function(x) {
return(paste0("dabs1(", x, ")"))
})
.rxD$dabs1 <- list(function(x) {
return("0")
})
.rxD$dabs <- list(function(x) {
return(paste0("dabs2(", x, ")"))
})
.rxD$dabs2 <- list(function(x) {
return("0")
})
.rxD$rxTBS <- list(function(a, lambda, yj, hi, low) {
paste0("rxTBSd(", a, ",", lambda, ",", yj, ",", hi, ",", low, ")")
})
.rxD$rxTBSd <- list(function(a, lambda, yj, hi, low) {
paste0("rxTBSd2(", a, ",", lambda, ",", yj, ",", hi, ",", low, ")")
})
.rxD$..k <- 10
.rxD$..tol <- 1e-4
## Approx a==b by
## (1-tanh(k*(a-b))^2)
.rxD$rxEq <- list(
function(a, b) {
.ab <- paste0("(", a, "-", b, ")")
return(paste0(
"(", -2 * .rxD$..k, "*tanh(", .rxD$..k, "*", .ab, ")+",
2 * .rxD$..k, "*tanh(", .rxD$..k, "*", .ab, ")^3)"
))
}, function(a, b) {
.ab <- paste0("(", a, "-", b, ")")
return(paste0(
"(", 2 * .rxD$..k, "*tanh(", .rxD$..k, "*", .ab, ")-",
2 * .rxD$..k, "*tanh(", .rxD$..k, "*", .ab, ")^3)"
))
}
)
.rxD$rxGeq <- list(
function(a, b) {
.delta <- atanh(2 * .rxD$..tol - 1)
## approx is (1/2+1/2*tanh(k*(a-b)-delta))
.ab <- paste0("(", a, "-", b, ")")
## (1/2)*k + (-1/2)*k*tanh(-delta + k*(a - b))^2
return(paste0(
"(", .rxD$..k / 2, "-", .rxD$..k / 2, "*tanh(", -.delta,
"+", .rxD$..k, "*", .ab, ")^2)"
))
}, function(a, b) {
.delta <- atanh(2 * .rxD$..tol - 1)
## approx is (1/2+1/2*tanh(k*(a-b)-delta))
.ab <- paste0("(", a, "-", b, ")")
## (1/2)*k + (-1/2)*k*tanh(-delta + k*(a - b))^2
return(paste0(
"(", -.rxD$..k / 2, "+", .rxD$..k / 2, "*tanh(", -.delta,
"+", .rxD$..k, "*", .ab, ")^2)"
))
}
)
.rxD$rxLeq <- list(
function(a, b) {
.delta <- atanh(2 * .rxD$..tol - 1)
## approx is (1/2-1/2*tanh(k*(a-b)+delta))
.ab <- paste0("(", a, "-", b, ")")
return(paste0(
"(", -.rxD$..k / 2, "+", .rxD$..k / 2, "*tanh(", .delta,
"+", .rxD$..k, "*", .ab, ")^2)"
))
}, function(a, b) {
.delta <- atanh(2 * .rxD$..tol - 1)
## approx is (1/2-1/2*tanh(k*(a-b)+delta))
.ab <- paste0("(", a, "-", b, ")")
return(paste0(
"(", .rxD$..k / 2, "-", .rxD$..k / 2, "*tanh(", .delta,
"+", .rxD$..k, "*", .ab, ")^2)"
))
}
)
.rxD$rxLt <- list(
function(a, b) {
## Approx is 1/2-1/2*tanh(k*(a-b)-delta)
.delta <- atanh(2 * .rxD$..tol - 1)
.ab <- paste0("(", a, "-", b, ")")
## (-1/2)*k + (1/2)*k*tanh(-delta + k*(a - b))^2
return(paste0(
"(", -.rxD$..k / 2, "+", .rxD$..k / 2, "*tanh(", -.delta,
"+", .rxD$..k, "*", .ab, ")^2)"
))
},
function(a, b) {
## Approx is 1/2-1/2*tanh(k*(a-b)-delta)
.delta <- atanh(2 * .rxD$..tol - 1)
.ab <- paste0("(", a, "-", b, ")")
## (-1/2)*k + (1/2)*k*tanh(-delta + k*(a - b))^2
return(paste0(
"(", .rxD$..k / 2, "-", .rxD$..k / 2, "*tanh(", -.delta,
"+", .rxD$..k, "*", .ab, ")^2)"
))
}
)
.rxD$rxGt <- list(
function(a, b) {
## delta <- atanh(2*tol-1);
## 1/2+1/2*tanh(k*(a-b)+delta)
.delta <- atanh(2 * .rxD$..tol - 1)
.ab <- paste0("(", a, "-", b, ")")
## (1/2)*k + (-1/2)*k*tanh(delta + k*(a - b))^2
return(paste0(
"(", .rxD$..k / 2, "-", .rxD$..k / 2, "*tanh(", .delta,
"+", .rxD$..k, "*", .ab, ")^2)"
))
},
function(a, b) {
## delta <- atanh(2*tol-1);
## 1/2+1/2*tanh(k*(a-b)+delta)
.delta <- atanh(2 * .rxD$..tol - 1)
.ab <- paste0("(", a, "-", b, ")")
## (-1/2)*k + (1/2)*k*tanh(delta + k*(a - b))^2
return(paste0(
"(", -.rxD$..k / 2, "+", .rxD$..k / 2, "*tanh(", .delta,
"+", .rxD$..k, "*", .ab, ")^2)"
))
}
)
.rxD$rxAnd <- list(
function(a, b) {
## a*b
return(b)
}, function(a, b) {
## a*b
return(a)
}
)
.rxD$rxOr <- list(
function(a, b) {
## Using DeMorgan's Theorem
## a+b = 1-(1-a)*(1-b)
return(paste0("(1-(", b, "))"))
}, function(a, b) {
return(paste0("(1-(", a, "))"))
}
)
.rxD$rxNot <- list(
function(a) {
## 1 - a
return("(-1)")
}
)
.rxD$tlast <- list(function(a) {
return("0")
})
.rxD$tfirst <- list(function(a) {
return("0")
})
.rxD$first <- list(function(a) {
return("0")
})
.rxD$last <- list(function(a) {
return("0")
})
.rxD$diff <- list(function(a) {
return("0")
})
.rxD$is.nan <- list(function(a) {
return("0")
})
.rxD$is.na <- list(function(a) {
return("0")
})
.rxD$is.finite <- list(function(a) {
return("0")
})
.rxD$is.infinite <- list(function(a) {
return("0")
})
## Approx a>=b by
## 1/2-1/2*tanh(k*x+delta)=1-tol
## 1/2-1+tol=1/2*tanh(k*x+delta)
## atanh(2*tol-1)= delta
## 1/2-1/2*tanh(k*(a-b)+delta)
.linCmtBgen <- function(i) {
.fun <- function(...) {}
body(.fun) <- bquote({
.args <- unlist(list(...))
if (.args[6] != "0") stop("cannot take a second derivative", call. = FALSE)
.args[6] <- .(i)
return(paste0("linCmtB(", paste(.args, collapse = ","), ")"))
})
return(.fun)
}
.rxD$linCmtB <- c(
list(
function(...) {
stop("bad 'linCmtB' derivative", call. = FALSE)
},
function(...) {
stop("bad 'linCmtB' derivative", call. = FALSE)
},
function(...) {
stop("bad 'linCmtB' derivative", call. = FALSE)
},
function(...) {
stop("bad 'linCmtB' derivative", call. = FALSE)
},
function(...) {
stop("bad 'linCmtB' derivative", call. = FALSE)
},
function(...) {
stop("bad 'linCmtB' derivative", call. = FALSE)
}
),
lapply(1:15, .linCmtBgen)
)
#' Add to RxODE's derivative tables
#'
#' @param name Function Name
#' @param derivatives A list of functions. Each function takes the
#' same number of arguments as the original function. The first
#' function will construct the derivative with respect to the
#' first argument; The second function will construct the
#' derivitive with respect to the second argument, and so on.
#' @return nothing
#' @author Matthew Fidler
#' @export
#' @examples
#' ## Add an arbitrary list of derivative functions
#' ## In this case the fun(x,y) is assumed to be 0.5*x^2+0.5*y^2
#'
#' rxD("fun", list(
#' function(x, y) {
#' return(x)
#' },
#' function(x, y) {
#' return(y)
#' }
#' ))
rxD <- function(name, derivatives) {
if (!inherits(derivatives, "list") || length(derivatives) == 0L) {
stop("derivatives must be a list of functions with at least 1 element", call. = FALSE)
}
if (!all(sapply(derivatives, function(x) (inherits(x, "function") || is.null(x))))) {
stop("derivatives must be a list of functions with at least 1 element", call. = FALSE)
}
if (exists(name, envir = .rxD)) {
warning(sprintf(gettext("replacing defined derivatives for '%s'"), name), call. = FALSE)
}
assign(name, derivatives, envir = .rxD)
return(invisible())
}
.promoteLinB <- FALSE
#' RxODE to symengine environment
#'
#' @param x expression
#'
#' @param envir default is `NULL`; Environment to put symengine
#' variables in.
#'
#' @param progress shows progress bar if true.
#'
#' @param promoteLinSens Promote solved linear compartment systems to
#' sensitivity-based solutions.
#'
#' @param unknownDerivatives When handling derivatives from unknown
#' functions, the translator will translate into different types
#' of numeric derivatives. The currently supported methods are:
#'
#' - `forward` for forward differences
#' - `central` for central differences
#' - `error` for throwing an error for unknown derivatives
#' @return An rxode symengine environment
#' @author Matthew L. Fidler
#' @export
rxToSE <- function(x, envir = NULL, progress = FALSE,
promoteLinSens = TRUE) {
assignInMyNamespace(".promoteLinB", promoteLinSens)
if (is(substitute(x), "character")) {
force(x)
} else if (is(substitute(x), "{")) {
x <- deparse1(substitute(x))
if (x[1] == "{") {
x <- x[-1]
x <- x[-length(x)]
}
x <- paste(x, collapse = "\n")
} else {
.xc <- as.character(substitute(x))
x <- substitute(x)
if (length(.xc) == 1) {
.found <- FALSE
.frames <- seq(1, sys.nframe())
.frames <- .frames[.frames != 0]
for (.f in .frames) {
.env <- parent.frame(.f)
if (exists(.xc, envir = .env)) {
.val2 <- try(get(.xc, envir = .env), silent = TRUE)
if (inherits(.val2, "character")) {
.val2 <- eval(parse(text = paste0("quote({", .val2, "})")))
return(.rxToSE(.val2, envir, progress))
} else if (inherits(.val2, "numeric") || inherits(.val2, "integer")) {
return(sprintf("%s", .val2))
}
}
}
}
return(.rxToSE(x, envir, progress))
}
return(.rxToSE(eval(parse(text = paste0("quote({", x, "})"))), envir, progress))
}
.rxChrToSym <- function(x) {
str2lang(paste0(
"rxQ__",
gsub(
" ", "_rxSpace_",
gsub("[.]", "_rxDoT_", x)
),
"__rxQ"
))
}
.rxRepRxQ <- function(x) {
.nchr <- nchar(x)
if (.nchr > 10) {
if (substr(x, 1, 5) == "rxQ__") {
return(deparse1(gsub(
"_rxSpace_", " ",
gsub(
"_rxDoT_", ".",
substr(x, 6, .nchr - 5)
)
)))
}
}
return(x)
}
#' @rdname rxToSE
#' @export
.rxToSE <- function(x, envir = NULL, progress = FALSE) {
rxReq("symengine")
.cnst <- names(.rxSEreserved)
.isEnv <- inherits(envir, "rxS") || inherits(envir, "environment")
if (is.name(x) || is.atomic(x)) {
if (is.character(x)) {
.ret <- .rxChrToSym(x)
if (.isEnv) {
.ret2 <- as.character(.ret)
assign(.ret2, symengine::Symbol(.ret2), envir = envir)
}
return(.ret)
} else {
.ret <- as.character(x)
if (any(.ret == .cnst)) {
.ret <- paste0("rx_SymPy_Res_", .ret)
if (.isEnv && is.name(x)) {
if (substr(x, 1, 1) != ".") {
if (!exists(.ret, envir = envir)) {
assign(.ret, symengine::Symbol(.ret), envir = envir)
}
}
}
return(.ret)
} else {
.ret0 <- .rxSEcnt[.ret]
if (is.na(.ret0)) {
if (.isEnv && is.name(x)) {
## message(.ret)
if (substr(x, 1, 1) != ".") {
if (!exists(.ret, envir = envir) && is.name(x)) {
assign(.ret, symengine::Symbol(.ret), envir = envir)
}
}
}
return(.ret)
}
return(setNames(.ret0, NULL))
}
}
} else if (is.call(x)) {
if (identical(x[[1]], quote(`(`))) {
return(paste0("(", .rxToSE(x[[2]], envir = envir), ")"))
} else if (identical(x[[1]], quote(`{`))) {
.x2 <- x[-1]
if (progress) {
rxProgress(length(.x2))
on.exit({
rxProgressAbort()
})
.ret <- paste(lapply(.x2, function(x) {
rxTick()
.rxToSE(x, envir = envir)
}), collapse = "\n")
rxProgressStop()
} else {
.ret <- paste(lapply(.x2, .rxToSE, envir = envir),
collapse = "\n"
)
}
## Assign and evaluate deferred items.
if (.isEnv) {
for (.var in names(envir$..ddt..)) {
.expr <- envir$..ddt..[[.var]]
.expr <- eval(parse(text = .expr))
assign(.var, .expr, envir = envir)
.rx <- paste0(
rxFromSE(.var), "=",
rxFromSE(.expr)
)
assign("..ddt", c(envir$..ddt, .rx),
envir = envir
)
}
for (.var in names(envir$..sens0..)) {
.expr <- envir$..sens0..[[.var]]
.expr <- eval(parse(text = .expr))
assign(.var, .expr, envir = envir)
.rx <- paste0(
rxFromSE(.var), "=",
rxFromSE(.expr)
)
assign("..sens0", c(envir$..sens0, .rx),
envir = envir
)
}
for (.var in names(envir$..jac0..)) {
.expr <- envir$..jac0..[[.var]]
.expr <- eval(parse(text = .expr))
assign(.var, .expr, envir = envir)
.rx <- paste0(
rxFromSE(.var), "=",
rxFromSE(.expr)
)
assign("..jac0", c(envir$..jac0, .rx),
envir = envir
)
}
}
return(.ret)
} else if (identical(x[[1]], quote(`*`)) ||
identical(x[[1]], quote(`^`)) ||
identical(x[[1]], quote(`+`)) ||
identical(x[[1]], quote(`-`)) ||
identical(x[[1]], quote(`/`))) {
if (length(x) == 3) {
if (identical(x[[1]], quote(`/`))) {
.x2 <- x[[2]]
.x3 <- x[[3]]
## df(%s)/dy(%s)
if (identical(.x2, quote(`d`)) &&
identical(.x3[[1]], quote(`dt`))) {
if (length(.x3[[2]]) == 1) {
.state <- as.character(.x3[[2]]) # .rxToSE(.x3[[2]], envir = envir)
} else {
.state <- .rxToSE(.x3[[2]], envir = envir)
}
return(paste0("rx__d_dt_", .state, "__"))
} else {
if (length(.x2) == 2 && length(.x3) == 2) {
if (identical(.x2[[1]], quote(`df`)) &&
identical(.x3[[1]], quote(`dy`))) {
if (length(.x2[[2]]) == 1) {
.state <- as.character(.x2[[2]])
} else {
.state <- .rxToSE(.x2[[2]], envir = envir)
}
if (length(.x3[[2]]) == 1) {
.var <- as.character(.x3[[2]])
} else {
.var <- .rxToSE(.x3[[2]], envir = envir)
}
return(paste0(
"rx__df_", .state,
"_dy_", .var, "__"
))
}
}
.ret <- paste0(
.rxToSE(.x2, envir = envir),
as.character(x[[1]]),
.rxToSE(.x3, envir = envir)
)
}
} else {
.ret <- paste0(
.rxToSE(x[[2]], envir = envir),
as.character(x[[1]]),
.rxToSE(x[[3]], envir = envir)
)
}
return(.ret)
} else {
## Unary Operators
return(paste(
as.character(x[[1]]),
.rxToSE(x[[2]], envir = envir)
))
}
} else if (identical(x[[1]], quote(`=`)) ||
identical(x[[1]], quote(`<-`)) ||
identical(x[[1]], quote(`~`))) {
.var <- .rxToSE(x[[2]], envir = envir)
.isNum <- FALSE
if (.isEnv) {
if (length(x[[2]]) == 2) {
if (any(as.character(x[[2]][[1]]) == c("alag", "lag", "F", "f", "rate", "dur"))) {
envir$..eventVars <- unique(c(.var, envir$..eventVars))
}
if (as.character(x[[2]][[1]]) == "mtime") {
envir$..mtimeVars <- unique(c(.var, envir$..mtimeVars))
}
}
}
if (inherits(x[[3]], "numeric") || inherits(x[[3]], "integer")) {
.isNum <- TRUE
if (.isEnv) {
if (envir$..doConst) {
assign(.var, x[[3]], envir = envir)
}
}
.expr <- x[[3]]
}
if (.isEnv) {
.expr <- paste0(
"with(envir,",
.rxToSE(x[[3]],
envir = envir
), ")"
)
if (regexpr(rex::rex(or(
.regRate,
.regDur,
.regLag,
.regF,
regIni0,
regDDt
)), .var) != -1) {
if (regexpr(
rex::rex(or(regSens, regSensEtaTheta)),
.var
) != -1) {
.lst <- get("..sens0..", envir = envir)
.lst[[.var]] <- .expr
assign("..sens0..", .lst, envir = envir)
} else {
.lst <- get("..ddt..", envir = envir)
.lst[[.var]] <- .expr
assign("..ddt..", .lst, envir = envir)
}
} else if (regexpr(rex::rex(or(
regDfDy,
regDfDyTh
)), .var) != -1) {
.lst <- get("..jac0..", envir = envir)
.lst[[.var]] <- .expr
assign("..jac0..", .lst, envir = envir)
} else if (!identical(x[[1]], quote(`~`))) {
.expr <- try(eval(parse(text = .expr)), silent = TRUE)
.isNum <- (inherits(.expr, "numeric") || inherits(.expr, "integer"))
if ((.isNum && envir$..doConst) ||
(!.isNum)) {
assign(.var, .expr, envir = envir)
}
.name <- rxFromSE(.var)
.rx <- paste0(
.name, "=",
rxFromSE(.expr)
)
if (regexpr("^(nlmixr|rx)_", .var) == -1) {
if (.isNum) {
names(.rx) <- .name
assign("..lhs0", c(envir$..lhs0, .rx),
envir = envir
)
} else {
if (any(names(envir$..lhs0) == .name)) {
.tmp <- envir$..lhs0
.tmp <- .tmp[names(.tmp) != .name]
assign("..lhs0", .tmp, envir = envir)
}
assign("..lhs", c(envir$..lhs, .rx),
envir = envir
)
}
}
} else {
.expr <- eval(parse(text = .expr))
if (envir$..doConst || !is.numeric(.expr)) {
assign(.var, .expr, envir = envir)
.rx <- paste0(
rxFromSE(.var), "=",
rxFromSE(.expr)
)
}
}
}
} else if (identical(x[[1]], quote(`[`))) {
.type <- toupper(as.character(x[[2]]))
if (any(.type == c("THETA", "ETA"))) {
if (is.numeric(x[[3]])) {
.num <- x[[3]]
if (round(.num) == .num) {
if (.num > 0) {
if (.isEnv) {
if (.type == "THETA") {
if (exists("..maxTheta", envir = envir)) {
.m <- get("..maxTheta", envir = envir)
} else {
.m <- 0
}
.m <- max(.m, .num)
.funs <- envir$..curCall
.funs <- .funs[.funs != ""]
if (length(.funs) > 0) {
.funs <- paste(.funs, collapse = ".")
envir$..extraTheta[[.funs]] <-
unique(c(
envir$..extraTheta[[.funs]],
.num
))
}
assign("..maxTheta", .m, envir = envir)
} else {
if (exists("..maxEta", envir = envir)) {
.m <- get("..maxEta", envir = envir)
} else {
.m <- 0
}
.m <- max(.m, .num)
.funs <- envir$..curCall
.funs <- .funs[.funs != ""]
if (length(.funs) > 0) {
.funs <- paste(.funs, collapse = ".")
envir$..extraEta[[.funs]] <-
unique(c(
envir$..extraEta[[.funs]],
.num
))
}
assign("..maxEta", .m, envir = envir)
}
}
return(paste0(.type, "_", .num, "_"))
} else {
stop("only 'THETA[#]' or 'ETA[#]' are supported", call. = FALSE)
}
} else {
stop("only 'THETA[#]' or 'ETA[#]' are supported", call. = FALSE)
}
} else {
stop("only 'THETA[#]' or 'ETA[#]' are supported", call. = FALSE)
}
} else {
stop("only 'THETA[#]' or 'ETA[#]' are supported", call. = FALSE)
}
} else if (identical(x[[1]], quote(`tad`))) {
.len <- length(x)
if (.len == 1L) {
} else if (.len == 2L) {
if (length(x[[2]]) != 1) {
stop(as.character(x[[1]]), "() must be used with a state", call. = FALSE)
}
return(paste0("(t-tlast(", as.character(x[[2]]), "))"))
} else {
stop(as.character(x[[1]]), "() can have 0-1 arguments", call. = FALSE)
}
return(paste0("(t-tlast())"))
} else if (identical(x[[1]], quote(`lag`)) ||
identical(x[[1]], quote(`lead`))) {
.len <- length(x)
.fun <- as.character(x[[1]])
if (.len == 1L) {
stop(.fun, "() takes 1-2 arguments")
} else if (.len == 2L) {
if (length(x[[2]]) != 1) {
stop(.fun, "() must be used with a variable", call. = FALSE)
}
return(paste0(.fun, "(", as.character(x[[2]]), ")"))
} else if (.len == 3L) {
if (length(x[[2]]) != 1) {
stop(.fun, "() must be used with a variable", call. = FALSE)
}
if (length(x[[3]]) != 1) {
stop(.fun, "(", as.character(x[[2]]), ", #) must have an integer for the number of lagged doses", call. = FALSE)
}
if (regexpr(rex::rex(maybe(one_of("-", "+")), regDecimalint), as.character(x[[3]]), perl = TRUE) == -1) {
stop(.fun, "(", as.character(x[[2]]), ", #) must have an integer for the number of lagged doses", call. = FALSE)
}
return(paste0(.fun, "(", as.character(x[[2]]), ", ", as.character(x[[3]]), ")"))
} else {
stop(as.character(x[[1]]), "() can have 0-1 arguments", call. = FALSE)
}
return(paste0("(t-tfirst())"))
} else if (identical(x[[1]], quote(`tafd`))) {
.len <- length(x)
if (.len == 1L) {
} else if (.len == 2L) {
if (length(x[[2]]) != 1) {
stop(as.character(x[[1]]), "() must be used with a state", call. = FALSE)
}
return(paste0("(t-tfirst(", as.character(x[[2]]), "))"))
} else {
stop(as.character(x[[1]]), "() can have 0-1 arguments", call. = FALSE)
}
return(paste0("(t-tfirst())"))
} else if (identical(x[[1]], quote(`tlast`)) ||
identical(x[[1]], quote(`tfirst`))) {
.len <- length(x)
if (.len == 1L) {
} else if (.len == 2L) {
if (length(x[[2]]) != 1) {
stop(as.character(x[[1]]), "() must be used with a state", call. = FALSE)
}
return(paste0(as.character(x[[1]]), "(", as.character(x[[2]]), ")"))
} else {
stop(as.character(x[[1]]), "() can have 0-1 arguments", call. = FALSE)
}
return(paste0(as.character(x[[1]]), "()"))
} else if (identical(x[[1]], quote(`psigamma`))) {
if (length(x == 3)) {
if (.isEnv) {
.lastCall <- envir$..curCall
envir$..curCall <- c(envir$..curCall, "psigamma")
}
.a <- .rxToSE(x[[2]], envir = envir)
.b <- .rxToSE(x[[3]], envir = envir)
if (.isEnv) envir$..curCall <- .lastCall
return(paste0("polygamma(", .b, ",", .a, ")"))
} else {
stop("'psigamma' takes 2 arguments", call. = FALSE)
}
} else if (identical(x[[1]], quote(`log1pmx`))) {
if (length(x == 2)) {
if (.isEnv) {
.lastCall <- envir$..curCall
envir$..curCall <- c(envir$..curCall, "log")
}
.a <- .rxToSE(x[[2]], envir = envir)
if (.isEnv) envir$..curCall <- .lastCall
return(paste0("(log(1+", .a, ")-(", .a, "))"))
} else {
stop("'log1pmx' only takes 1 argument", call. = FALSE)
}
} else if (identical(x[[1]], quote(`choose`))) {
if (length(x) == 3) {
if (.isEnv) {
.lastCall <- envir$..curCall
envir$..curCall <- c(envir$..curCall, "gamma")
}
.n <- .rxToSE(x[[2]], envir = envir)
.k <- .rxToSE(x[[3]], envir = envir)
if (.isEnv) envir$..curCall <- .lastCall
return(paste0(
"gamma(", .n, "+1)/(gamma(",
.k, "+1)*gamma(", .n, "-(", .k, ")+1))"
))
} else {
stop("'choose' takes 2 arguments", call. = FALSE)
}
} else if (identical(x[[1]], quote(`lchoose`))) {
if (length(x) == 3) {
if (.isEnv) {
.lastCall <- envir$..curCall
envir$..curCall <- c(envir$..curCall, "lgamma")
}
.n <- .rxToSE(x[[2]], envir = envir)
.k <- .rxToSE(x[[3]], envir = envir)
if (.isEnv) envir$..curCall <- .lastCall
return(paste0("(lgamma(", .n, "+1)-lgamma(", .k, "+1)-lgamma(", .n, "-(", .k, ")+1))"))
} else {
stop("'lchoose' takes 2 arguments", call. = FALSE)
}
} else if ((identical(x[[1]], quote(`pnorm`))) |
(identical(x[[1]], quote(`normcdf`))) |
(identical(x[[1]], quote(`phi`)))) {
if (length(x) == 4) {
## pnorm(q, mean, sd)
if (.isEnv) {
.lastCall <- envir$..curCall
envir$..currCall <- c(envir$..curCall, "erf")
}
.q <- .rxToSE(x[[2]], envir = envir)
.mean <- .rxToSE(x[[3]], envir = envir)
.sd <- .rxToSE(x[[4]], envir = envir)
return(paste0("0.5*(1+erf((((", .q, ")-(", .mean, "))/(", .sd, "))/sqrt(2)))"))
} else if (length(x) == 3) {
## pnorm(q, mean)
if (.isEnv) {
.lastCall <- envir$..curCall
envir$..currCall <- c(envir$..curCall, "erf")
}
.q <- .rxToSE(x[[2]], envir = envir)
.mean <- .rxToSE(x[[3]], envir = envir)
return(paste0("0.5*(1+erf((((", .q, ")-(", .mean, ")))/sqrt(2)))"))
} else if (length(x) == 2) {
## pnorm(q)
if (.isEnv) {
.lastCall <- envir$..curCall
envir$..currCall <- c(envir$..curCall, "erf")
}
.q <- .rxToSE(x[[2]], envir = envir)
return(paste0("0.5*(1+erf((", .q, ")/sqrt(2)))"))
} else {
stop("'pnorm' can only take 1-3 arguments", call. = FALSE)
}
} else if (identical(x[[1]], quote(`transit`))) {
if (length(x) == 4) {
## transit(n, mtt, bio)
if (.isEnv) {
.lastCall <- envir$..curCall
envir$..curCall <- c(envir$..curCall, "lgamma")
}
.n <- .rxToSE(x[[2]], envir = envir)
if (.isEnv) {
envir$..curCall <- .lastCall
.lastCall <- envir$..curCall
envir$..curCall <- c(envir$..curCall, "log")
}
.mtt <- .rxToSE(x[[3]], envir = envir)
.bio <- .rxToSE(x[[4]], envir = envir)
if (.isEnv) envir$..curCall <- .lastCall
return(paste0(
"exp(log((", .bio, ")*(podo))+log(",
.n, " + 1)-log(", .mtt, ")+(", .n,
")*((log(", .n, "+1)-log(", .mtt,
"))+log(t))-((", .n, "+1)/(", .mtt,
"))*(t)-lgamma(1+", .n, "))"
))
} else if (length(x) == 3) {
if (.isEnv) {
.lastCall <- envir$..curCall
envir$..curCall <- c(envir$..curCall, "lgamma")
}
.n <- .rxToSE(x[[2]], envir = envir)
.mtt <- .rxToSE(x[[3]], envir = envir)
if (.isEnv) envir$..curCall <- .lastCall
return(paste0("exp(log(podo)+(log(", .n, "+1)-log(", .mtt, "))+(", .n, ")*((log(", .n, "+1)-log(", .mtt, "))+ log(t))-((", .n, " + 1)/(", .mtt, "))*(t)-lgamma(1+", .n, "))"))
} else {
stop("'transit' can only take 2-3 arguments", call. = FALSE)
}
} else {
if (length(x[[1]]) == 1) {
.x1 <- as.character(x[[1]])
.xc <- .rxSEsingle[[.x1]]
if (!is.null(.xc)) {
if (length(x) == 2) {
if (.isEnv) {
.lastCall <- envir$..curCall
envir$..curCall <- c(envir$..curCall, .xc[3])
}
.ret <- paste0(
.xc[1], .rxToSE(x[[2]], envir = envir),
.xc[2]
)
if (.isEnv) envir$..curCall <- .lastCall
return(.ret)
} else {
stop(sprintf("'%s' only acceps 1 argument", .x1), call. = FALSE)
}
}
.xc <- .rxSEdouble[[.x1]]
if (!is.null(.xc)) {
.x1 <- as.character(x[[1]])
if (length(x) == 3) {
.ret <- paste0(
.xc[1], .rxToSE(x[[2]], envir = envir),
.xc[2],
.rxToSE(x[[3]], envir = envir),
.xc[3]
)
return(.ret)
} else {
stop(sprintf("'%s' only acceps 2 arguments", .x1), call. = FALSE)
}
}
}
if (.isEnv) {
.lastCall <- envir$..curCall
envir$..curCall <- c(
envir$..curCall,
as.character(x[[1]])
)
}
.ret0 <- c(list(as.character(x[[1]])), lapply(x[-1], .rxToSE, envir = envir))
if (.isEnv) envir$..curCall <- .lastCall
.SEeq <- c(.rxSEeq, .rxSEeqUsr)
.curName <- paste(.ret0[[1]])
.nargs <- .SEeq[.curName]
if (.promoteLinB && .curName == "linCmtA") {
.ret0 <- c(
list("linCmtB"),
.ret0[2:6],
list("0"),
.ret0[-c(1, 2:6)]
)
.nargs <- .nargs + 1
}
if (!is.na(.nargs)) {
if (.nargs == length(.ret0) - 1) {
.ret <- paste0(.ret0[[1]], "(")
.ret0 <- .ret0[-1]
.ret <- paste0(.ret, paste(unlist(.ret0), collapse = ","), ")")
if (.ret == "exp(1)") {
return("E")
}
return(.ret)
} else {
stop(sprintf(
gettext("'%s' takes %s arguments (has %s)"),
paste(.ret0[[1]]),
.nargs, length(.ret0) - 1
), call. = FALSE)
}
} else {
.fun <- paste(.ret0[[1]])
.ret0 <- .ret0[-1]
if (length(.ret0) == 1L) {
if (any(.fun == c("alag", "lag"))) {
return(paste0("rx_lag_", .ret0[[1]], "_"))
} else if (any(.fun == c("F", "f"))) {
return(paste0("rx_f_", .ret0[[1]], "_"))
} else if (any(.fun == c("rate", "dur"))) {
return(paste0("rx_", .fun, "_", .ret0[[1]], "_"))
}
}
.ret <- paste0("(", paste(unlist(.ret0), collapse = ","), ")")
if (.ret == "(0)") {
return(paste0("rx_", .fun, "_ini_0__"))
} else if (any(.fun == c("cmt", "dvid"))) {
return("")
} else if (any(.fun == c("max", "min"))) {
.ret <- paste0(.fun, "(", paste(unlist(.ret0), collapse = ","), ")")
} else if (.fun == "sum") {
.ret <- paste0("(", paste(paste0("(", unlist(.ret0), ")"), collapse = "+"), ")")
} else if (.fun == "prod") {
.ret <- paste0("(", paste(paste0("(", unlist(.ret0), ")"), collapse = "*"), ")")
} else if (.fun == "probitInv") {
## erf <- function(x) 2 * pnorm(x * sqrt(2)) - 1 (probitInv=pnorm)
if (length(.ret0) == 1) {
.ret <- paste0("0.5*(1+erf((", unlist(.ret0)[1], ")/sqrt(2)))")
} else if (length(.ret0) == 2) {
.ret0 <- unlist(.ret0)
.p <- paste0("0.5*(1+erf((", .ret0[1], ")/sqrt(2)))")
## return (high-low)*p+low;
.ret <- paste0(
"(1.0-(", .ret0[2], "))*(", .p,
")+(", .ret0[2], ")"
)
} else if (length(.ret0) == 3) {
.ret0 <- unlist(.ret0)
.p <- paste0("0.5*(1+erf((", .ret0[1], ")/sqrt(2)))")
.ret <- paste0(
"((", .ret0[3], ")-(", .ret0[2], "))*(", .p,
")+(", .ret0[2], ")"
)
} else {
stop("'probitInv' requires 1-3 arguments",
call. = FALSE
)
}
} else if (.fun == "probit") {
## erfinv <- function (x) qnorm((1 + x)/2)/sqrt(2) (probit=qnorm )
if (length(.ret0) == 1) {
.ret <- paste0("sqrt(2)*erfinv(2*(", unlist(.ret0), ")-1)")
} else if (length(.ret0) == 2) {
.ret0 <- unlist(.ret0)
.p <- paste0(
"((", .ret0[1], ")-(", .ret0[2], "))/(1.0-",
"(", .ret0[2], "))"
)
.ret <- paste0("sqrt(2)*erfinv(2*(", .p, ")-1)")
} else if (length(.ret0) == 3) {
.ret0 <- unlist(.ret0)
## (x-low)/(high-low)
.p <- paste0(
"((", .ret0[1], ")-(", .ret0[2],
"))/((", .ret0[3], ")-(", .ret0[2], "))"
)
.ret <- paste0("sqrt(2)*erfinv(2*(", .p, ")-1)")
} else {
stop("'probit' requires 1-3 arguments",
call. = FALSE
)
}
} else if (.fun == "logit") {
if (length(.ret0) == 1) {
.ret <- paste0("-log(1/(", unlist(.ret0), ")-1)")
} else if (length(.ret0) == 2) {
.ret0 <- unlist(.ret0)
.p <- paste0(
"((", .ret0[1], ")-(", .ret0[2], "))/(1.0-",
"(", .ret0[2], "))"
)
.ret <- paste0("-log(1/(", .p, ")-1)")
} else if (length(.ret0) == 3) {
.ret0 <- unlist(.ret0)
## (x-low)/(high-low)
.p <- paste0(
"((", .ret0[1], ")-(", .ret0[2],
"))/((", .ret0[3], ")-(", .ret0[2], "))"
)
.ret <- paste0("-log(1/(", .p, ")-1)")
} else {
stop("'logit' requires 1-3 arguments",
call. = FALSE
)
}
} else if (any(.fun == c("expit", "invLogit", "logitInv"))) {
if (length(.ret0) == 1) {
.ret <- paste0("1/(1+exp(-(", unlist(.ret0)[1], ")))")
} else if (length(.ret0) == 2) {
.ret0 <- unlist(.ret0)
.p <- paste0("1/(1+exp(-(", .ret0[1], ")))")
## return (high-low)*p+low;
.ret <- paste0(
"(1.0-(", .ret0[2], "))*(", .p,
")+(", .ret0[2], ")"
)
} else if (length(.ret0) == 3) {
.ret0 <- unlist(.ret0)
.p <- paste0("1/(1+exp(-(", .ret0[1], ")))")
.ret <- paste0(
"((", .ret0[3], ")-(", .ret0[2], "))*(", .p,
")+(", .ret0[2], ")"
)
} else {
stop("'expit' requires 1-3 arguments",
call. = FALSE
)
}
} else {
stop(sprintf(gettext("function '%s' or its derivatives are not supported in RxODE"), .fun),
call. = FALSE
)
}
}
}
} else {
stop("unsupported expression", call. = FALSE)
}
}
.rxUnXi <- function(x) {
gsub("_xi_([[:digit:]]*)", "rx_xi_\\1", x)
}
## 0L = error
## 1L = forward
## 2L = central
.rxFromNumDer <- 0L
.rxDelta <- (.Machine$double.eps)^(1 / 3)
.rxD$gammap <- list(
NULL,
function(a, z) {
paste0("gammapDer(", a, ",", z, ")")
}
)
#' @rdname rxToSE
#' @export
rxFromSE <- function(x, unknownDerivatives = c("forward", "central", "error")) {
rxReq("symengine")
.unknown <- c("central" = 2L, "forward" = 1L, "error" = 0L)
assignInMyNamespace(".rxFromNumDer", .unknown[match.arg(unknownDerivatives)])
if (is(substitute(x), "character")) {
return(.rxFromSE(eval(parse(text = paste0("quote({", .rxUnXi(x), "})")))))
} else if (is(substitute(x), "{")) {
x <- deparse1(substitute(x))
if (x[1] == "{") {
x <- x[-1]
x <- x[-length(x)]
}
x <- .rxUnXi(paste(x, collapse = "\n"))
.ret <- .rxFromSE(eval(parse(text = paste0("quote({", x, "})"))))
return(.ret)
} else {
.xc <- as.character(substitute(x))
x <- substitute(x)
if (length(.xc) == 1) {
.found <- FALSE
.frames <- seq(1, sys.nframe())
.frames <- .frames[.frames != 0]
for (.f in .frames) {
.env <- parent.frame(.f)
if (exists(.xc, envir = .env)) {
.val2 <- try(get(.xc, envir = .env), silent = TRUE)
if (inherits(.val2, "character")) {
.val2 <- eval(parse(text = paste0("quote({", .rxUnXi(.val2), "})")))
.ret <- .rxFromSE(.val2)
return(.ret)
} else if (inherits(.val2, "numeric") || inherits(.val2, "integer")) {
return(sprintf("%s", .val2))
} else {
if (!is.null(attr(class(.val2), "package"))) {
if (attr(class(.val2), "package") == "symengine") {
.val2 <- eval(parse(text = paste0("quote({", .rxUnXi(as.character(.val2)), "})")))
.ret <- .rxFromSE(.val2)
return(.ret)
}
}
}
}
}
}
x <- eval(parse(text = paste("quote(", .rxUnXi(paste(deparse1(x), collapse = " ")), ")")))
.ret <- .rxFromSE(x)
return(.ret)
}
x <- .rxUnXi(x)
.ret <- .rxFromSE(eval(parse(text = paste0("quote({", x, "})"))))
return(.ret)
}
.stripP <- function(x) {
if (is.call(x)) {
if (length(x) == 1) {
return(x)
} else if (identical(x[[1]], quote(`(`))) {
return(.stripP(x[[2]]))
} else {
return(x)
}
} else {
return(x)
}
}
.rxM1rmF <- function(x) {
.env <- new.env(parent = emptyenv())
.env$found <- FALSE
.f <- function(x, envir) {
if (is.call(x)) {
if (identical(x[[1]], quote(`*`))) {
if (length(x) == 3) {
.x2 <- as.character(x[[2]])
.x3 <- as.character(x[[3]])
if (length(.x3) == 1) {
if (.x3 == "pi") {
envir$found <- TRUE
return(.rxFromSE(.stripP(x[[2]])))
}
}
if (length(.x2) == 1) {
if (.x2 == "pi") {
envir$found <- TRUE
return(.rxFromSE(.stripP(x[[3]])))
}
}
return(paste0(
.f(x[[2]], envir), "*",
.f(x[[3]], envir)
))
} else {
return(.rxFromSE(x))
}
} else if (identical(x[[1]], quote(`/`))) {
.x2 <- as.character(x[[2]])
.x3 <- as.character(x[[3]])
if (length(.x2) == 1) {
if (.x2 == "pi") {
envir$found <- TRUE
x[[2]] <- 1
return(.rxFromSE(x))
}
return(paste0(
.f(x[[2]], envir), "/",
.f(x[[3]], envir)
))
}
} else {
return(.rxFromSE(x))
}
} else {
return(.rxFromSE(x))
}
}
.ret <- .f(x, .env)
return(list(.ret, .env$found))
}
.rxP1rmF <- function(x) {
.env <- new.env(parent = emptyenv())
.env$found <- FALSE
.f <- function(x, envir) {
if (is.call(x)) {
if (identical(x[[1]], quote(`+`))) {
if (length(x) == 3) {
if (inherits(x[[3]], "numeric")) {
if (x[[3]] == 1) {
envir$found <- TRUE
return(.rxFromSE(.stripP(x[[2]])))
} else {
return(paste0(
.f(x[[2]], envir),
"+", as.character(x[[3]])
))
}
}
if (inherits(x[[2]], "numeric")) {
if (x[[2]] == 1) {
envir$found <- TRUE
return(.rxFromSE(.stripP(x[[3]])))
} else {
return(paste0(
as.character(x[[2]]), "+",
.f(x[[3]], envir)
))
}
}
return(paste0(
.f(x[[2]], envir), "+",
.f(x[[3]], envir)
))
} else {
return(.rxFromSE(x))
}
} else if (identical(x[[1]], quote(`-`))) {
if (length(x) == 3) {
if (inherits(x[[2]], "numeric")) {
if (x[[2]] == 1) {
x <- x[-2]
envir$found <- TRUE
return(.rxFromSE(x))
} else {
return(.rxFromSE(x))
}
} else {
return(.rxFromSE(x))
}
} else {
return(.rxFromSE(x))
}
} else {
return(.rxFromSE(x))
}
} else {
return(.rxFromSE(x))
}
}
.ret <- .f(x, .env)
return(list(.ret, .env$found))
}
.rxFromSEnum <- function(x) {
.ret <- as.character(x)
.retl <- nchar(.ret)
if (.retl > 5) {
.op <- options()
options(digits = 22)
on.exit(options(.op))
.rx <- names(.rxSEcnt)
.val <- setNames(.rxSEcnt, NULL)
E <- exp(1)
for (.i in seq_along(.rx)) {
.tmp <- try(eval(parse(text = paste0("substr(paste(", .val[.i], "), 1, .retl)"))), silent = TRUE)
if (!inherits(.tmp, "try-error")) {
if (.ret == .tmp) {
return(.rx[.i])
}
}
}
}
return(.ret)
}
#' @export
#' @rdname rxToSE
.rxFromSE <- function(x) {
rxReq("symengine")
.cnst <- setNames(
names(.rxSEreserved),
paste0("rx_SymPy_Res_", names(.rxSEreserved))
)
if (is.name(x) || is.atomic(x)) {
.ret <- .rxFromSEnum(x)
.ret <- .rxRepRxQ(.ret)
.ret0 <- .cnst[.ret]
if (!is.na(.ret0)) {
return(.ret0)
}
.ret0 <- .rxSEreserved[[.ret]]
if (!is.null(.ret0)) {
if (is.character(.ret0)) {
return(.ret0)
} else if (is.numeric(.ret0)) {
return(sprintf("%.16f", .ret0))
}
}
if (.ret == "E") {
return("M_E")
}
.ret <- sub(
.regRate, "rate(\\1)",
sub(
.regDur, "dur(\\1)",
sub(
.regLag, "alag(\\1)",
sub(
.regF, "f(\\1)",
sub(
regIni0, "\\1(0)",
sub(
regDfDy, "df(\\1)/dy(\\2)",
sub(
regDfDyTh, "df(\\1)/dy(\\2[\\3])",
sub(
regDDt, "d/dt(\\1)",
sub(
rex::rex(start, regThEt, end),
"\\1[\\2]", .ret
)
)
)
)
)
)
)
)
)
.ret <- sub("[(]rx_SymPy_Res_", "(", .ret)
return(.ret)
} else if (is.call(x)) {
if (identical(x[[1]], quote(`(`))) {
return(paste0("(", .rxFromSE(x[[2]]), ")"))
} else if (identical(x[[1]], quote(`{`))) {
.x2 <- x[-1]
return(paste(lapply(.x2, function(x) {
.ret <- .rxFromSE(x)
return(.ret)
}),
collapse = "\n"
))
} else if (identical(x[[1]], quote(`*`)) ||
identical(x[[1]], quote(`^`)) ||
identical(x[[1]], quote(`+`)) ||
identical(x[[1]], quote(`-`)) ||
identical(x[[1]], quote(`/`))) {
if (length(x) == 3) {
.x1 <- as.character(x[[1]])
.x2 <- x[[2]]
.x2 <- .rxFromSE(.x2)
.x3 <- x[[3]]
.x3 <- .rxFromSE(.x3)
.x3v <- try(eval(parse(text = .x3)), silent = TRUE)
if (inherits(.x3v, "numeric")) {
.x3 <- .rxFromSEnum(.x3v)
}
if (.x1 == "^" && .x3 == "1") {
return(.x2)
}
if (.x1 == "^" && .x3 == "-1") {
return(paste0("(1/(", .x2, "))"))
}
if (.x1 == "^" && is.numeric(x[[3]])) {
if (round(x[[3]]) == x[[3]]) {
return(paste0("Rx_pow_di(", .x2, ",", .x3, ")"))
}
if (.x3 == 0.5) {
if (any(.x2 == c("pi", "M_PI"))) {
return("M_SQRT_PI")
} else if (any(.x2 == c("M_2_PI", "(M_2_PI)"))) {
return("M_SQRT_2dPI")
} else {
return(paste0("sqrt(", .x2, ")"))
}
} else {
return(paste0("Rx_pow(", .x2, ",", .x3, ")"))
}
}
.ret <- paste0(.x2, .x1, .x3)
## FIXME parsing to figure out if *2 or *0.5 *0.4 is in
## expression
if (any(.ret == c("pi*2", "2*pi", "M_PI*2", "2*M_PI"))) {
return("M_2PI")
}
if (any(.ret == c("pi/2", "pi*0.5", "0.5*pi", "M_PI/2", "M_PI*0.5", "0.5*M_PI"))) {
return("M_PI_2")
}
if (any(.ret == c("pi/4", "pi*0.25", "0.25*pi", "M_PI/4", "M_PI*0.25", "0.25*M_PI"))) {
return("M_PI_4")
}
if (any(.ret == c("1/pi", "1/M_PI"))) {
return("M_1_PI")
}
if (any(.ret == c("2/pi", "2/M_PI"))) {
return("M_2_PI")
}
if (any(.ret == c(
"(M_2_PI)^0.5", "(M_2_PI)^(1/2)",
"M_2_PI^0.5", "M_2_PI^(1/2)",
"sqrt((M_2_PI))"
))) {
return("M_SQRT_2dPI")
}
if (any(.ret == c(
"(pi)^0.5", "(pi)^(1/2)",
"pi^0.5", "pi^(1/2)",
"(M_PI)^0.5", "(M_PI)^(1/2)",
"M_PI^0.5", "M_PI^(1/2)"
))) {
return("M_SQRT_PI")
}
if (.ret == "log(2)/log(10)") {
return("M_LOG10_2")
}
if (.ret == "1/log(10)") {
return("M_LOG10E")
}
if (.ret == "1/log(2)") {
return("M_LOG2E")
}
if (any(.ret == c("2/M_SQRT_PI", "2/(M_SQRT_PI)"))) {
return("M_2_SQRTPI")
}
if (any(.ret == c(
"1/sqrt(M_2PI)",
"1/(sqrt((M_2PI)))",
"1/(M_2PI^0.5)", "1/(M_2PI^(1/2))",
"1/((M_2PI)^0.5)", "1/((M_2PI)^(1/2))"
))) {
return("M_1_SQRT_2PI")
}
## if (.x1 == "^") {
## return(paste0("Rx_pow(", .x2, ",", .x3, ")"))
## }
return(.ret)
} else {
## Unary Operators
.ret <- paste0(
as.character(x[[1]]),
.rxFromSE(x[[2]])
)
return(.ret)
}
} else if (identical(x[[1]], quote(`=`)) ||
identical(x[[1]], quote(`<-`)) ||
identical(x[[1]], quote(`~`))) {
.var <- .rxFromSE(x[[2]])
.val <- .rxFromSE(x[[3]])
} else if (identical(x[[1]], quote(`[`))) {
if (any(as.character(x[[2]]) == c("THETA", "ETA"))) {
return(paste0(x[[2]], "[", x[[3]], "]"))
}
stop("[...] expressions not supported",
call. = FALSE
)
} else if (identical(x[[1]], quote(`lag`)) ||
identical(x[[1]], quote(`lead`))) {
.a <- .rxFromSE(x[[2]])
.fun <- as.character(x[[1]])
if (length(x) == 3) {
return(paste0(.fun, "(", .a, ",", .rxFromSE(x[[3]]), ")"))
} else {
return(paste0(.fun, "(", .a, ")"))
}
} else if (identical(x[[1]], quote(`polygamma`))) {
if (length(x == 3)) {
.a <- .rxFromSE(x[[2]])
.b <- .rxFromSE(x[[3]])
if (.a == "0") {
return(paste0("digamma(", .b, ")"))
} else if (.a == "1") {
return(paste0("trigamma(", .b, ")"))
} else if (.a == "2") {
return(paste0("tetragamma(", .b, ")"))
} else if (.a == "3") {
return(paste0("pentagamma(", .b, ")"))
} else {
return(paste0(
"psigamma(", .b, ",",
.a, ")"
))
}
} else {
stop("'polygamma' takes 2 arguments",
call. = FALSE
)
}
} else {
if (length(x[[1]]) == 1) {
.x1 <- as.character(x[[1]])
.xc <- .SEsingle[[.x1]]
if (!is.null(.xc)) {
if (length(x) == 2) {
.x2 <- x[[2]]
if (length(.x2) != 1) {
if (identical(.x2[[1]], quote(`+`))) {
.tmp0 <- .SE1p[.x1]
if (!is.na(.tmp0)) {
.ret <- .rxP1rmF(.x2)
if (.ret[[2]]) {
.r1 <- .ret[[1]]
return(paste0(
.tmp0, "(",
.r1,
")"
))
}
}
return(paste0(.xc[1], .x2[[1]], .xc[2]))
}
}
return(paste0(.xc[1], .rxFromSE(x[[2]]), .xc[2]))
} else {
stop(sprintf("'%s' only acceps 1 argument", .x1),
call. = FALSE
)
}
}
.xc <- .SEdouble[[.x1]]
if (!is.null(.xc)) {
if (length(x) == 3) {
.x1 <- .rxFromSE(x[[1]])
return(paste0(
.xc[1], .rxFromSE(x[[2]]), .xc[2],
.rxFromSE(x[[3]]),
.xc[3]
))
} else {
stop(sprintf("'%s' only acceps 2 arguments", .x1),
call. = FALSE
)
}
}
}
if (length(x) == 2) {
if (identical(x[[1]], quote(`log`))) {
if (length(x[[2]]) == 3) {
if (identical(x[[2]][[1]], quote(`beta`))) {
.tmp <- x[[2]]
.tmp[[1]] <- quote(`lbeta`)
return(.rxFromSE(.tmp))
}
}
}
}
.ret0 <- lapply(lapply(x, .stripP), .rxFromSE)
.SEeq <- c(.rxSEeq, .rxSEeqUsr)
.nargs <- .SEeq[paste(.ret0[[1]])]
if (!is.na(.nargs)) {
if (.nargs == length(.ret0) - 1) {
.x1 <- as.character(.ret0[[1]])
.tmp0 <- .x1
if (.nargs == 1) {
.tmp0 <- .SE1p[.x1]
.x2 <- x[[2]]
if (!is.na(.tmp0)) {
.ret <- .rxP1rmF(.x2)
if (.ret[[2]]) {
if (.tmp0 == "log1p") {
.tmp <- eval(parse(text = paste0("quote(", .ret[[1]], ")")))
if (length(.tmp) > 1) {
if (identical(.tmp[[1]], quote(`exp`))) {
.tmp <- .tmp[[-1]]
.tmp0 <- "log1pexp"
.ret[[1]] <- .rxFromSE(.tmp)
}
}
}
return(paste0(
.tmp0, "(",
.ret[[1]],
")"
))
} else {
.ret <- paste0(
.x1, "(",
.ret[[1]],
")"
)
if (.ret == "log(2)") {
return("M_LN2")
}
if (.ret == "log(10)") {
return("M_LN10")
}
if (.ret == "log(M_SQRT_PI)") {
return("M_LN_SQRT_PI")
}
if (any(.ret == c(
"log(sqrt((M_PI_2)))",
"log(sqrt(M_PI_2))",
"log((M_PI_2)^(1/2))",
"log((M_PI_2)^0.5)",
"log(M_PI_2^(1/2))",
"log(M_PI_2^0.5)"
))) {
return("M_LN_SQRT_PId2")
}
if (any(.ret == c(
"log(sqrt((M_2PI)))",
"log(sqrt(M_2PI))",
"log((M_2PI)^0.5)",
"log((M_2PI)^(1/2))",
"log(M_2PI^0.5)",
"log(M_2PI^(1/2))"
))) {
return("M_LN_SQRT_2PI")
}
return(.ret)
}
}
.tmp0 <- .SE1m[.x1]
if (!is.na(.tmp0)) {
.ret <- .rxM1rmF(.x2)
if (.ret[[2]]) {
return(paste0(
.tmp0, "(",
.ret[[1]],
")"
))
} else {
return(paste0(
.x1, "(",
.ret[[1]],
")"
))
}
}
.tmp0 <- .x1
}
.ret <- paste0(.tmp0, "(")
.ret0 <- .ret0[-1]
.ret <- paste0(
.ret, paste(unlist(.ret0), collapse = ","),
")"
)
if (.ret == "exp(1)") {
return("M_E")
}
if (.ret == "sin(pi)") {
return("0")
}
if (.ret == "cos(pi)") {
return("1")
}
if (.ret == "tan(pi)") {
return("0")
}
if (.ret == "sqrt(3)") {
return("M_SQRT_3")
}
if (.ret == "sqrt(2)") {
return("M_SQRT_2")
}
if (.ret == "sqrt(32)") {
return("M_SQRT_32")
}
if (.ret == "sqrt(pi)") {
return("M_SQRT_PI")
}
if (any(.ret == c("sqrt(M_2_PI)", "sqrt((M_2_PI))"))) {
return("M_SQRT_2dPI")
}
return(.ret)
} else {
stop(sprintf(
gettext("'%s' takes %s arguments"),
paste(.ret0[[1]]),
.nargs
))
}
} else if (identical(x[[1]], quote(`Derivative`))) {
if (length(x) == 3) {
.fun <- as.character(x[[2]])
.var <- .rxFromSE(x[[3]])
if (length(.fun) == 1) {
if (.fun == "abs0") {
return(paste0("abs(", .var, ")"))
}
}
.args <- .fun[-1]
.args <- lapply(.args, .rxFromSE)
.with <- which(.var == .args)
.errD <- function(force = FALSE) {
if (!force && .rxFromNumDer != 0L) {
## Can calculate forward or central
## difference instead.
## Warn
if (.rxFromNumDer == 1L) {
## Forward
.a1 <- .args
.fn <- .fun[1]
.a1[.with] <- paste0("(", .a1[.with], ")+", .rxDelta)
.a2 <- .args
return(paste0(
"(", .fn, "(", paste0(.a1, collapse = ","), ")-",
.fn, "(", paste0(.a2, collapse = ","), "))/", .rxDelta
))
} else if (.rxFromNumDer == 2L) {
## Central
.a1 <- .args
.fn <- .fun[1]
.a1[.with] <- paste0(.a1[.with], "-", (0.5 * .rxDelta))
.a2 <- .args
.a2[.with] <- paste0(.a2[.with], "+", (0.5 * .rxDelta))
return(paste0(
"(", .fn, "(", paste0(.a1, collapse = ","), ")-",
.fn, "(", paste0(.a2, collapse = ","), "))/", .rxDelta
))
} else {
stop("only forward and central differences are supported", call. = FALSE)
}
} else {
stop(sprintf(gettext("cannot figure out the '%s' derivative with respect to '%s'"), .fun[1], .var[1]))
}
}
if (length(.with) != 1) {
.errD(force = TRUE)
}
if (any(.fun[1] == c("lead", "lag"))) {
return("0")
}
if (exists(.fun[1], envir = .rxD)) {
.funLst <- get(.fun[1], envir = .rxD)
if (length(.funLst) < .with) {
return(.errD())
}
.derFun <- .funLst[[.with]]
if (is.null(.derFun)) {
return(.errD())
}
.ret <- try(do.call(.derFun, as.list(.args)), silent = TRUE)
if (inherits(.ret, "try-error")) {
warning(
"an error occurred looking up the derivative for '",
.fun[1], "' using numerical differences instead"
)
return(.errD())
} else {
return(.ret)
}
} else {
if (.rxFromNumDer == 0L) {
stop(sprintf(gettext("RxODE/symengine does not know how to take a derivative of '%s'"), .fun[1]),
call. = FALSE
)
} else {
return(.errD())
}
}
} else {
stop("'Derivative' conversion only takes one function and one argument",
call. = FALSE
)
}
} else if (identical(x[[1]], quote(`Subs`))) {
.fun <- eval(parse(text = paste0("quote(", .rxFromSE(x[[2]]), ")")))
.what <- .stripP(x[[3]])
.with <- .stripP(x[[4]])
.subs <- function(x) {
if (identical(x, .what)) {
return(.with)
} else if (is.call(x)) {
as.call(lapply(x, .subs))
} else if (is.pairlist(x)) {
as.pairlist(lapply(x, .subs))
} else {
return(x)
}
}
.ret <- .subs(.fun)
return(.rxFromSE(.ret))
} else if (any(paste(.ret0[[1]]) == c("max", "min"))) {
.x1 <- as.character(.ret0[[1]])
.ret <- paste0(.x1, "(")
.ret0 <- .ret0[-1]
.ret <- paste0(
.ret, paste(unlist(.ret0), collapse = ","),
")"
)
return(.ret)
} else if (any(paste(.ret0[[1]]) == c("tlast", "tfirst"))) {
if (length(.ret0) == 1L) {
return(paste0(.ret0[[1]], "()"))
} else if (length(.ret0) == 2L) {
if (length(.ret0[[2]]) == 1L) {
return(paste0(.ret0[[1]], "(", .ret0[[2]], ")"))
}
}
stop(paste0(.ret0[[1]], "() takes 0-1 arguments"))
} else {
stop(sprintf(gettext("'%s' not supported in symengine->RxODE"), paste(.ret0[[1]])),
call. = FALSE
)
}
}
} else {
stop("unsupported expression", call. = FALSE)
}
}
.rxFunction <- function(name) {
.f <- function(...) {
1
}
body(.f) <- bquote(return(symengine::FunctionSymbol(.(name), unlist(list(...)))))
return(.f)
}
#' Load a model into a symengine environment
#'
#' @param x RxODE object
#' @param doConst Load constants into the environment as well.
#' @inheritParams rxToSE
#' @return RxODE/symengine environment
#' @author Matthew Fidler
#' @export
rxS <- function(x, doConst = TRUE, promoteLinSens = FALSE) {
rxReq("symengine")
.cnst <- names(.rxSEreserved)
.env <- new.env(parent = loadNamespace("symengine"))
.env$..mv <- rxModelVars(x)
.env$..jac0 <- NULL
.env$..jac0.. <- list()
.env$..ddt <- NULL
.env$..ddt.. <- list()
.env$..sens0 <- NULL
.env$..sens0.. <- list()
.env$..lhs <- NULL
.env$..lhs0 <- NULL
.env$..doConst <- doConst
for (.f in c(
ls(.rxD), "linCmtA", "linCmtB", "rxEq", "rxNeq", "rxGeq", "rxLeq", "rxLt",
"rxGt", "rxAnd", "rxOr", "rxNot", "rxTBS", "rxTBSd", "rxTBSd2", "lag", "lead"
)) {
assign(.f, .rxFunction(.f), envir = .env)
}
for (.v in seq_along(.rxSEreserved)) {
assign(names(.rxSEreserved)[.v], .rxSEreserved[[.v]], envir = .env)
}
.env$..s0 <- symengine::S("0")
.env$..extraTheta <- list()
.env$..extraEta <- list()
.env$..curCall <- character(0)
.env$..eventVars <- NULL
.env$..mtimeVars <- NULL
.env$polygamma <- function(a, b) {
## symengine::subs(symengine::subs(..polygamma, ..a, a), ..b, b)
symengine::psigamma(b, a)
}
.pars <- c(
rxParams(x), rxState(x),
"podo", "t", "time", "tlast", "rx1c", "rx__PTR__"
)
## default lambda/yj values
.env$rx_lambda_ <- symengine::S("1")
.env$rx_yj_ <- symengine::S("2")
.env$rx_low_ <- symengine::S("0")
.env$rx_hi_ <- symengine::S("1")
if (!is.null(.rxSEeqUsr)) {
sapply(names(.rxSEeqUsr), function(x) {
assign(.rxFunction(x), x, envir = .env)
})
}
## EulerGamma=0.57721566490153286060651209008240243104215933593992
## S("I")
## S("pi")
## S("E")
## S("") # EulerGamma
## S("Catalan") = 0.915965594177219015054603514932384110774
## S("GoldenRatio") = 1+sqrt(5)/2
## S("inf")
## S("nan")
sapply(.pars, function(x) {
if (any(.cnst == x)) {
.tmp <- paste0("rx_SymPy_Res_", x)
assign(.tmp, symengine::Symbol(.tmp), envir = .env)
} else {
.tmp <- rxToSE(x)
assign(.tmp, symengine::Symbol(.tmp), envir = .env)
assign(x, symengine::Symbol(x), envir = .env)
}
})
assignInMyNamespace(".promoteLinB", promoteLinSens)
.expr <- eval(parse(text = paste0("quote({", rxNorm(x), "})")))
.ret <- .rxToSE(.expr, .env)
class(.env) <- "rxS"
return(.env)
}
symengineC <- new.env(parent = emptyenv())
symengineC$"**" <- .dslToPow
symengineC$"^" <- .dslToPow
symengineC$S <- function(x) {
sprintf("%s", x)
}
for (f in names(.rxSEeq)) {
symengineC[[f]] <- functionOp(f)
}
symengineC$"(" <- unaryOp("(", ")")
for (op in c("+", "-", "*")) {
symengineC[[op]] <- binaryOp(paste0(" ", op, " "))
}
symengineC[["/"]] <- function(e1, e2) {
sprintf("%s /( (%s == 0) ? %s : %s)", e1, e2, .Machine$double.eps, e2)
}
unknownCsymengine <- function(op) {
force(op)
function(...) {
stop(sprintf("RxODE doesn't support '%s' translation for 'Omega' translation", op),
call. = FALSE
)
}
}
symengineCEnv <- function(expr) {
## Known functions
calls <- allCalls(expr)
callList <- setNames(lapply(calls, unknownCsymengine), calls)
callEnv <- list2env(callList)
rxSymPyFEnv <- cloneEnv(symengineC, callEnv)
names <- allNames(expr)
## Replace time with t.
n1 <- names
n2 <- names
n2 <- gsub(rex::rex("t", capture(numbers)), "REAL(theta)[\\1]", n2)
n2 <- gsub(rex::rex("pi"), "M_PI", n2)
n2 <- gsub(rex::rex("rx_SymPy_Res_"), "", n2)
n2 <- gsub("None", "NA_REAL", n2)
w <- n2[n2 == "t"]
symbol.list <- setNames(as.list(n2), n1)
symbol.env <- list2env(symbol.list, parent = symengineC)
return(symbol.env)
}
seC <- function(x) {
expr <- eval(parse(text = sprintf("quote(%s)", as.character(x))))
.ret <- eval(expr, symengineCEnv(expr))
}
## nocov end
sympyTransit4 <- function(t, n, mtt, bio, podo = "podo", tlast = "tlast") {
ktr <- paste0("((", n, " + 1)/(", mtt, "))")
lktr <- paste0("(log((", n, ") + 1) - log(", mtt, "))")
tc <- paste0("((", t, ")-(", tlast, "))")
paste0(
"exp(log((", bio, ") * (", podo, ")) + ", lktr, " + (",
n, ") * ", "(", lktr, " + log(", t, ")) - ",
ktr, " * (", t, ") - log(gamma(1 + (", n, "))))"
)
}
allNames <- function(x) {
if (is.atomic(x)) {
character()
} else if (is.name(x)) {
as.character(x)
} else if (is.call(x) || is.pairlist(x)) {
children <- lapply(x[-1], allNames)
unique(unlist(children))
} else {
stop("do not know how to handle type ", typeof(x),
call. = FALSE
)
}
}
allCalls <- function(x) {
if (is.atomic(x) || is.name(x)) {
character()
} else if (is.call(x)) {
fname <- as.character(x[[1]])
children <- lapply(x[-1], allCalls)
unique(c(fname, unlist(children)))
} else if (is.pairlist(x)) {
unique(unlist(lapply(x[-1], allCalls), use.names = FALSE))
} else {
stop("do not know how to handle type ", typeof(x), call. = FALSE)
}
}
cloneEnv <- function(env, parent = parent.env(env)) {
list2env(as.list(env), parent = parent)
}
.exists2 <- function(x, where) {
.nc <- try(nchar(x) < 1000, silent = TRUE)
if (inherits(.nc, "try-error")) .nc <- FALSE
if (rxIs(.nc, "logical")) .nc <- FALSE
if (.nc) {
return(exists(x, where))
} else {
return(FALSE)
}
}
## Start error function DSL
rxErrEnvF <- new.env(parent = emptyenv())
for (op in c(
"+", "-", "*", "/", "^", "**",
"!=", "==", "&", "&&", "|", "||"
)) {
op2 <- op
if (op == "**") {
op2 <- "^"
}
rxErrEnvF[[op]] <- binaryOp(paste0(" ", op2, " "))
}
for (op in c("=", "~", "<-")) {
rxErrEnvF[[op]] <- binaryOp(" = ")
}
rxErrEnvF$"{" <- function(...) {
return(sprintf("{\n%s\n}", paste(unlist(list(...)), collapse = "\n")))
}
rxErrEnvF$"(" <- unaryOp("(", ")")
rxErrEnvF$"[" <- function(name, val) {
n <- toupper(name)
err <- gettext("RxODE only supports THETA[#] and ETA[#] numbers")
if (any(n == c("THETA", "ETA")) && is.numeric(val)) {
if (round(val) == val && val > 0) {
if (n == "THETA" && as.numeric(val) <= length(rxErrEnv.init)) {
return(sprintf("THETA[%s]", val))
} else {
return(sprintf("%s[%s]", n, val))
}
} else {
stop(err)
}
} else {
stop(err)
}
}
rxErrEnvF$"if" <- function(lg, tr, fl) {
if (missing(fl)) {
return(sprintf("if (%s) %s", lg, tr))
} else {
return(sprintf("if (%s) %s else %s", lg, tr, fl))
}
}
rxErrEnv.theta <- 1
rxErrEnv.diag.est <- NULL
rxErrEnv.ret <- "rx_r_"
rxErrEnv.init <- NULL
rxErrEnv.lambda <- NULL
rxErrEnv.yj <- NULL
rxErrEnv.combined <- "^2"
rxErrEnv.hasAdd <- FALSE
rxErrEnv.hi <- "1"
rxErrEnv.low <- "0"
.rxErrEnvInit <- function() {
assignInMyNamespace("rxErrEnv.hasAdd", FALSE)
assignInMyNamespace("rxErrEnv.theta", 1)
assignInMyNamespace("rxErrEnv.diag.est", NULL)
assignInMyNamespace("rxErrEnv.ret", "rx_r_")
assignInMyNamespace("rxErrEnv.init", NULL)
assignInMyNamespace("rxErrEnv.lambda", NULL)
assignInMyNamespace("rxErrEnv.yj", NULL)
assignInMyNamespace("rxErrEnv.combined", "^2")
assignInMyNamespace("rxErrEnv.hi", "1")
assignInMyNamespace("rxErrEnv.low", "0")
}
rxErrEnvF$lnorm <- function(est) {
if (rxErrEnv.ret != "rx_r_") {
stop("'lnorm' can only be in an error function", call. = FALSE)
}
if (!is.null(rxErrEnv.lambda)) {
if (rxErrEnv.lambda != "0" && rxErrEnv.yj != "0") {
stop("'lnorm' cannot be used with other data transformations", call. = FALSE)
}
}
if (is.na(est)) {
assignInMyNamespace("rxErrEnv.lambda", "0")
assignInMyNamespace("rxErrEnv.yj", "0")
return("")
} else {
estN <- suppressWarnings(as.numeric(est))
assignInMyNamespace("rxErrEnv.hasAdd", TRUE)
if (is.na(estN)) {
ret <- (sprintf("(%s)%s", est, rxErrEnv.combined))
assignInMyNamespace("rxErrEnv.lambda", "0")
assignInMyNamespace("rxErrEnv.yj", "0")
} else {
theta <- sprintf("THETA[%s]", rxErrEnv.theta)
est <- estN
theta.est <- theta
ret <- (sprintf("(%s)%s", theta.est, rxErrEnv.combined))
tmp <- rxErrEnv.diag.est
tmp[sprintf("THETA[%s]", rxErrEnv.theta)] <- as.numeric(est)
assignInMyNamespace("rxErrEnv.diag.est", tmp)
assignInMyNamespace("rxErrEnv.theta", rxErrEnv.theta + 1)
assignInMyNamespace("rxErrEnv.lambda", "0")
assignInMyNamespace("rxErrEnv.yj", "0")
}
}
return(ret)
}
rxErrEnvF$dlnorm <- rxErrEnvF$lnorm
rxErrEnvF$logn <- rxErrEnvF$lnorm
rxErrEnvF$logitNorm <- function(est, low = "0", hi = "1") {
if (rxErrEnv.ret != "rx_r_") {
stop("'logitNorm' can only be in an error function", call. = FALSE)
}
if (!is.null(rxErrEnv.lambda)) {
if (rxErrEnv.yj != "1") {
if (rxErrEnv.yj != "4" & rxErrEnv.yj != "5") {
stop("'logitNorm' cannot be used with other data transformations", call. = FALSE)
}
}
}
if (is.na(est)) {
if (is.null(rxErrEnv.yj)) {
assignInMyNamespace("rxErrEnv.yj", "4")
} else if (rxErrEnv.yj == "1") {
assignInMyNamespace("rxErrEnv.yj", "5")
} else {
assignInMyNamespace("rxErrEnv.yj", "4")
}
assignInMyNamespace("rxErrEnv.hi", hi)
assignInMyNamespace("rxErrEnv.low", low)
return("")
} else {
estN <- suppressWarnings(as.numeric(est))
assignInMyNamespace("rxErrEnv.hasAdd", TRUE)
if (is.na(estN)) {
ret <- (sprintf("(%s)%s", est, rxErrEnv.combined))
assignInMyNamespace("rxErrEnv.lambda", "0")
if (is.null(rxErrEnv.yj)) {
assignInMyNamespace("rxErrEnv.yj", "4")
} else if (rxErrEnv.yj == "1") {
assignInMyNamespace("rxErrEnv.yj", "5")
} else {
assignInMyNamespace("rxErrEnv.yj", "4")
}
assignInMyNamespace("rxErrEnv.hi", hi)
assignInMyNamespace("rxErrEnv.low", low)
} else {
theta <- sprintf("THETA[%s]", rxErrEnv.theta)
est <- estN
theta.est <- theta
ret <- (sprintf("(%s)%s", theta.est, rxErrEnv.combined))
tmp <- rxErrEnv.diag.est
tmp[sprintf("THETA[%s]", rxErrEnv.theta)] <- as.numeric(est)
assignInMyNamespace("rxErrEnv.diag.est", tmp)
assignInMyNamespace("rxErrEnv.theta", rxErrEnv.theta + 1)
assignInMyNamespace("rxErrEnv.lambda", "0")
if (is.null(rxErrEnv.yj)) {
assignInMyNamespace("rxErrEnv.yj", "4")
} else if (rxErrEnv.yj == "1") {
assignInMyNamespace("rxErrEnv.yj", "5")
} else {
assignInMyNamespace("rxErrEnv.yj", "4")
}
assignInMyNamespace("rxErrEnv.hi", hi)
assignInMyNamespace("rxErrEnv.low", low)
}
}
return(ret)
}
rxErrEnvF$probitNorm <- function(est, low = "0", hi = "1") {
if (rxErrEnv.ret != "rx_r_") {
stop("'probitNorm' can only be in an error function", call. = FALSE)
}
if (!is.null(rxErrEnv.lambda)) {
if (rxErrEnv.yj != "1") {
if (rxErrEnv.yj != "6" & rxErrEnv.yj != "7") {
print(rxErrEnv.yj)
stop("'probitNorm' cannot be used with other data transformations", call. = FALSE)
}
}
}
if (is.na(est)) {
if (is.null(rxErrEnv.yj)) {
assignInMyNamespace("rxErrEnv.yj", "6")
} else if (rxErrEnv.yj == "1") {
assignInMyNamespace("rxErrEnv.yj", "7")
} else {
assignInMyNamespace("rxErrEnv.yj", "6")
}
assignInMyNamespace("rxErrEnv.hi", hi)
assignInMyNamespace("rxErrEnv.low", low)
return("")
} else {
estN <- suppressWarnings(as.numeric(est))
assignInMyNamespace("rxErrEnv.hasAdd", TRUE)
if (is.na(estN)) {
ret <- (sprintf("(%s)%s", est, rxErrEnv.combined))
assignInMyNamespace("rxErrEnv.lambda", "0")
if (is.null(rxErrEnv.yj)) {
assignInMyNamespace("rxErrEnv.yj", "6")
} else if (rxErrEnv.yj == "1") {
assignInMyNamespace("rxErrEnv.yj", "7")
} else {
assignInMyNamespace("rxErrEnv.yj", "6")
}
assignInMyNamespace("rxErrEnv.hi", hi)
assignInMyNamespace("rxErrEnv.low", low)
} else {
theta <- sprintf("THETA[%s]", rxErrEnv.theta)
est <- estN
theta.est <- theta
ret <- (sprintf("(%s)%s", theta.est, rxErrEnv.combined))
tmp <- rxErrEnv.diag.est
tmp[sprintf("THETA[%s]", rxErrEnv.theta)] <- as.numeric(est)
assignInMyNamespace("rxErrEnv.diag.est", tmp)
assignInMyNamespace("rxErrEnv.theta", rxErrEnv.theta + 1)
assignInMyNamespace("rxErrEnv.lambda", "0")
if (is.null(rxErrEnv.yj)) {
assignInMyNamespace("rxErrEnv.yj", "6")
} else if (rxErrEnv.yj == "1") {
assignInMyNamespace("rxErrEnv.yj", "7")
} else {
assignInMyNamespace("rxErrEnv.yj", "6")
}
assignInMyNamespace("rxErrEnv.hi", hi)
assignInMyNamespace("rxErrEnv.low", low)
}
}
return(ret)
}
rxErrEnvF$tbs <- function(lambda) {
if (rxErrEnv.ret != "rx_r_") {
stop("'boxCox' can only be in an error function", call. = FALSE)
}
if (!is.null(rxErrEnv.lambda)) {
if (rxErrEnv.yj != "0" & rxErrEnv.lambda != "0" & rxErrEnv.lambda != "1") {
stop("'boxCox' cannot be used with other data transformations", call. = FALSE)
}
}
estN <- suppressWarnings(as.numeric(lambda))
if (is.na(estN)) {
assignInMyNamespace("rxErrEnv.lambda", lambda)
assignInMyNamespace("rxErrEnv.yj", "0")
} else {
tmp <- rxErrEnv.diag.est
tmp[sprintf("THETA[%s]", rxErrEnv.theta)] <- estN
assignInMyNamespace("rxErrEnv.lambda", sprintf("THETA[%s]", rxErrEnv.theta))
assignInMyNamespace("rxErrEnv.diag.est", tmp)
assignInMyNamespace("rxErrEnv.theta", rxErrEnv.theta + 1)
assignInMyNamespace("rxErrEnv.yj", "0")
}
return("0")
}
rxErrEnvF$boxCox <- rxErrEnvF$tbs
rxErrEnvF$tbsYj <- function(lambda) {
if (rxErrEnv.ret != "rx_r_") {
stop("'yeoJohnson' can only be in an error function", call. = FALSE)
}
if (!is.null(rxErrEnv.lambda)) {
if ((rxErrEnv.yj != "1" & rxErrEnv.yj != "4" & rxErrEnv.yj != "6")) {
stop("'yeoJohnson' cannot be used with other data transformations", call. = FALSE)
}
}
estN <- suppressWarnings(as.numeric(lambda))
if (is.na(estN)) {
assignInMyNamespace("rxErrEnv.lambda", lambda)
if (is.null(rxErrEnv.yj)) {
assignInMyNamespace("rxErrEnv.yj", "1")
} else if (rxErrEnv.yj == "4") {
assignInMyNamespace("rxErrEnv.yj", "5")
} else if (rxErrEnv.yj == "6") {
assignInMyNamespace("rxErrEnv.yj", "7")
} else {
assignInMyNamespace("rxErrEnv.yj", "1")
}
} else {
tmp <- rxErrEnv.diag.est
tmp[sprintf("THETA[%s]", rxErrEnv.theta)] <- estN
assignInMyNamespace("rxErrEnv.lambda", sprintf("THETA[%s]", rxErrEnv.theta))
assignInMyNamespace("rxErrEnv.diag.est", tmp)
assignInMyNamespace("rxErrEnv.theta", rxErrEnv.theta + 1)
if (is.null(rxErrEnv.yj)) {
assignInMyNamespace("rxErrEnv.yj", "1")
} else if (rxErrEnv.yj == "4") {
assignInMyNamespace("rxErrEnv.yj", "5")
} else if (rxErrEnv.yj == "6") {
assignInMyNamespace("rxErrEnv.yj", "7")
} else {
assignInMyNamespace("rxErrEnv.yj", "1")
}
}
return("0")
}
rxErrEnvF$yeoJohnson <- rxErrEnvF$tbsYj
rxErrEnvF$add <- function(est) {
if (rxErrEnv.ret != "rx_r_") {
stop("'add' can only be in an error function", call. = FALSE)
}
assignInMyNamespace("rxErrEnv.hasAdd", TRUE)
estN <- suppressWarnings(as.numeric(est))
if (is.na(estN)) {
ret <- (sprintf("(%s)%s", est, rxErrEnv.combined))
} else {
theta <- sprintf("THETA[%s]", rxErrEnv.theta)
est <- estN
theta.est <- theta
ret <- (sprintf("(%s)%s", theta.est, rxErrEnv.combined))
tmp <- rxErrEnv.diag.est
tmp[sprintf("THETA[%s]", rxErrEnv.theta)] <- as.numeric(est)
assignInMyNamespace("rxErrEnv.diag.est", tmp)
assignInMyNamespace("rxErrEnv.theta", rxErrEnv.theta + 1)
}
return(ret)
}
rxErrEnvF$norm <- rxErrEnvF$add
rxErrEnvF$dnorm <- rxErrEnvF$add
rxErrEnvF$"for" <- function(...) {
stop("'for' is not supported", call. = FALSE)
}
rxErrEnvF$`return` <- function(est) {
if (rxErrEnv.ret == "") {
stop("The PK function should not return anything", call. = FALSE)
}
.extra <- ""
force(est)
if (rxErrEnv.ret == "rx_r_") {
.hi <- rxErrEnv.hi
.low <- rxErrEnv.low
if (is.null(rxErrEnv.lambda)) {
.lambda <- "1"
} else {
.lambda <- rxErrEnv.lambda
}
if (is.null(rxErrEnv.yj)) {
.yj <- "0"
} else {
.yj <- rxErrEnv.yj
}
if (.yj == "0" & .lambda == "1") {
.yj <- "2"
.lambda <- "1"
}
if (.yj == "0" & .lambda == "0") {
.yj <- "3"
.lambda <- "0"
}
.extra <- sprintf("rx_yj_~%s;\nrx_lambda_~%s;\nrx_hi_~%s\nrx_low_~%s\n", .yj, .lambda, .hi, .low)
}
return(sprintf("%s%s=%s;", .extra, rxErrEnv.ret, est))
}
rxErrEnvF$`|` <- binaryOp(" | ")
rxErrEnvF$`||` <- binaryOp(" || ")
rxErrEnvF$`&&` <- binaryOp(" && ")
rxErrEnvF$`<=` <- binaryOp(" <= ")
rxErrEnvF$`>=` <- binaryOp(" >= ")
rxErrEnvF$`==` <- binaryOp(" == ")
rxErrEnvF$prop <- function(est) {
if (rxErrEnv.ret != "rx_r_") {
stop("'prop' can only be in an error function")
}
estN <- suppressWarnings(as.numeric(est))
if (is.na(estN)) {
ret <- (sprintf("(rx_pred_f_)%s * (%s)%s", rxErrEnv.combined, est, rxErrEnv.combined))
} else {
est <- estN
ret <- ""
theta <- sprintf("THETA[%s]", rxErrEnv.theta)
theta.est <- theta
ret <- (sprintf("(rx_pred_f_)%s*(%s)%s", rxErrEnv.combined, theta.est, rxErrEnv.combined))
tmp <- rxErrEnv.diag.est
tmp[sprintf("THETA[%s]", rxErrEnv.theta)] <- as.numeric(est)
assignInMyNamespace("rxErrEnv.diag.est", tmp)
assignInMyNamespace("rxErrEnv.theta", rxErrEnv.theta + 1)
}
return(ret)
}
rxErrEnvF$propT <- function(est) {
if (rxErrEnv.ret != "rx_r_") {
stop("'propT' can only be in an error function")
}
estN <- suppressWarnings(as.numeric(est))
if (is.na(estN)) {
ret <- (sprintf("(rx_pred_)%s * (%s)%s", rxErrEnv.combined, est, rxErrEnv.combined))
} else {
est <- estN
ret <- ""
theta <- sprintf("THETA[%s]", rxErrEnv.theta)
theta.est <- theta
ret <- (sprintf("(rx_pred_)%s*(%s)%s", rxErrEnv.combined, theta.est, rxErrEnv.combined))
tmp <- rxErrEnv.diag.est
tmp[sprintf("THETA[%s]", rxErrEnv.theta)] <- as.numeric(est)
assignInMyNamespace("rxErrEnv.diag.est", tmp)
assignInMyNamespace("rxErrEnv.theta", rxErrEnv.theta + 1)
}
return(ret)
}
rxErrEnvF$pow <- function(est, pow) {
if (rxErrEnv.ret != "rx_r_") {
stop("'pow' can only be in an error function", call. = FALSE)
}
estN <- suppressWarnings(as.numeric(est))
if (is.na(estN)) {
ret <- (sprintf(
"(rx_pred_f_)^(%s%s) * (%s)%s", ifelse(rxErrEnv.combined == "^2", "2*", ""),
pow, est, rxErrEnv.combined
))
} else {
est <- estN
ret <- ""
theta <- sprintf("THETA[%s]", rxErrEnv.theta)
theta2 <- sprintf("THETA[%s]", rxErrEnv.theta + 1)
theta.est <- theta
ret <- (sprintf("(rx_pred_f_)^(%s%s) * (%s)%s", ifelse(rxErrEnv.combined == "^2", "2*", ""), theta2, theta.est, rxErrEnv.combined))
tmp <- rxErrEnv.diag.est
tmp[sprintf("THETA[%s]", rxErrEnv.theta)] <- as.numeric(est)
tmp[sprintf("THETA[%s]", rxErrEnv.theta + 1)] <- as.numeric(pow)
assignInMyNamespace("rxErrEnv.diag.est", tmp)
assignInMyNamespace("rxErrEnv.theta", rxErrEnv.theta + 2)
}
return(ret)
}
rxErrEnvF$powT <- function(est, pow) {
if (rxErrEnv.ret != "rx_r_") {
stop("'powT' can only be in an error function", call. = FALSE)
}
estN <- suppressWarnings(as.numeric(est))
if (is.na(estN)) {
ret <- (sprintf(
"(rx_pred_)^(%s%s) * (%s)%s", ifelse(rxErrEnv.combined == "^2", "2*", ""),
pow, est, rxErrEnv.combined
))
} else {
est <- estN
ret <- ""
theta <- sprintf("THETA[%s]", rxErrEnv.theta)
theta2 <- sprintf("THETA[%s]", rxErrEnv.theta + 1)
theta.est <- theta
ret <- (sprintf("(rx_pred_)^(%s%s) * (%s)%s", ifelse(rxErrEnv.combined == "^2", "2*", ""), theta2, theta.est, rxErrEnv.combined))
tmp <- rxErrEnv.diag.est
tmp[sprintf("THETA[%s]", rxErrEnv.theta)] <- as.numeric(est)
tmp[sprintf("THETA[%s]", rxErrEnv.theta + 1)] <- as.numeric(pow)
assignInMyNamespace("rxErrEnv.diag.est", tmp)
assignInMyNamespace("rxErrEnv.theta", rxErrEnv.theta + 2)
}
return(ret)
}
.convStr <- function(x) {
if (is.atomic(x) || is.name(x)) {
if (is.character(x)) {
.rxChrToSym(x)
} else {
return(x)
}
} else if (is.call(x) || is.pairlist(x)) {
return(as.call(lapply(x, .convStr)))
} else {
stop("do not know how to handle type ", typeof(x),
call. = FALSE
)
}
}
rxErrEnv <- function(expr) {
calls <- allCalls(expr)
callList <- setNames(lapply(calls, functionOp), calls)
callEnv <- list2env(callList)
## Known functions
rxErrFEnv <- cloneEnv(rxErrEnvF, callEnv)
## Symbols
names <- allNames(expr)
n1 <- names
n2 <- names
n2[n2 == "time"] <- "t"
if (any(n2 == "err")) {
stop("use 'return' for errors", call. = FALSE)
}
if (any(n2 == "error")) {
stop("use 'return' for errors", call. = FALSE)
}
if (any(n2 == "rx_r")) {
stop("use 'return' for errors", call. = FALSE)
}
## n2[n2 == "err"] <- "rx_r_";
## n2[n2 == "error"] <- "rx_r_";
## n2[n2 == "f"] <- "rx_pred_f_";
symbol.list <- setNames(as.list(n2), n1)
symbol.env <- list2env(symbol.list, parent = rxErrFEnv)
return(symbol.env)
}
#' Parse PK function for inclusion in RxODE
#'
#' @param x PK function
#' @inheritParams rxParseErr
#' @return RxODE transformed text.
#' @author Matthew L. Fidler
#' @keywords internal
#' @export
rxParsePk <- function(x, init = NULL) {
return(rxParseErr(x, init = init, ret = ""))
}
#' Prepare Pred function for inclusion in RxODE
#'
#' @param x pred function
#' @inheritParams rxParseErr
#' @return RxODE transformed text.
#' @author Matthew L. Fidler
#' @keywords internal
#' @export
rxParsePred <- function(x, init = NULL, err = NULL, addProp = c("combined2", "combined1")) {
if (is.null(err)) {
return(rxParseErr(x, ret = "rx_pred_", init = init, addProp = addProp))
} else {
.ini <- attr(err, "ini")
.errs <- rxExpandIfElse(rxGetModel(err))
.prd <- rxParseErr(x, ret = "rx_pred_", init = init)
.prd <- rxExpandIfElse(rxGetModel(.prd))
if (length(.prd) > 1) {
if (length(.prd) == length(.errs)) {
.prd <- .prd[names(.errs)]
if (any(is.na(.prd))) {
stop("errors and predictions need to have the same conditions ('if'/'then' statements)",
call. = FALSE
)
}
} else if (length(.errs) != 1) {
stop("do not know how to handle this error/pred combination",
call. = FALSE
)
}
}
.ret <- sapply(seq(1, max(length(.errs), length(.prd))), function(en) {
.e <- .errs[min(length(.errs), en)]
.p <- .prd[min(length(.prd), en)]
.reg <- rex::rex("rx_pred_", any_spaces, "=", any_spaces, capture(except_any_of(";\n")), any_of(";\n"))
if (regexpr(rex::rex("rx_yj_~2;\nrx_lambda_~1;\n"), .e) != -1) {
.p <- gsub(.reg, "rx_pred_f_~\\1;\nrx_pred_ = \\1", .p)
} else if (regexpr(rex::rex("rx_yj_~3;\nrx_lambda_~0;\n"), .e) != -1) {
.p <- gsub(.reg, "rx_pred_f_~\\1;\nrx_pred_ = log(\\1)", .p)
} else {
.p <- gsub(.reg, "rx_pred_f_~\\1;\nrx_pred_ = rxTBS(\\1, rx_lambda_, rx_yj_, rx_low_, rx_hi_)", .p)
}
return(gsub("rx_r_", sprintf("%s\nrx_r_", .p), .e))
})
if (length(.ret) == 1L) {
.ret <- setNames(.ret, NULL)
attr(.ret, "ini") <- .ini
return(.ret)
} else {
if (length(.errs) >= length(.prd)) {
.n <- names(.errs)
} else {
.n <- names(.prd)
}
.ord <- order(sapply(.n, nchar))
.n <- .n[.ord]
.ret <- .ret[.ord]
.ret <- paste(sapply(seq_along(.n), function(x) {
if (x > 1) {
if (.n[x] == sprintf("(!%s)", .n[x - 1])) {
return(sprintf("else {\n%s\n}", .ret[x]))
}
}
return(sprintf("if %s {\n%s\n}", .n[x], .ret[x]))
}), collapse = "\n")
attr(.ret, "ini") <- .ini
return(.ret)
}
}
}
#' Prepare Error function for inclusion in RxODE
#'
#' @param x error function
#' @param baseTheta Base theta to start numbering add(.) and prop(.) from.
#' @param ret Internal return type. Should not be changed by the user...
#' @param init Initialization vector
#' @return RxODE transformed text
#' @keywords internal
#' @author Matthew L. Fidler
#' @export
rxParseErr <- function(x, baseTheta, ret = "rx_r_", init = NULL,
addProp = c("combined2", "combined1")) {
addProp <- match.arg(addProp)
if (!missing(baseTheta)) {
assignInMyNamespace("rxErrEnv.theta", baseTheta)
}
if (!missing(ret)) {
assignInMyNamespace("rxErrEnv.ret", ret)
}
if (!missing(init)) {
assignInMyNamespace("rxErrEnv.init", init)
}
if (!missing(init)) {
assignInMyNamespace("rxErrEnv.init", init)
}
if (addProp == "combined2") {
assignInMyNamespace("rxErrEnv.combined", "^2")
} else {
assignInMyNamespace("rxErrEnv.combined", "")
}
if (is(x, "function")) {
x <- rxAddReturn(x, ret != "")
}
if (is(substitute(x), "character")) {
ret <- eval(parse(text = sprintf("RxODE:::rxParseErr(quote({%s}),addProp=\"%s\")", x, addProp)))
ret <- substring(ret, 3, nchar(ret) - 2)
assignInMyNamespace("rxErrEnv.diag.est", NULL)
assignInMyNamespace("rxErrEnv.theta", 1)
assignInMyNamespace("rxErrEnv.ret", "rx_r_")
assignInMyNamespace("rxErrEnv.init", NULL)
return(ret)
} else if (is(substitute(x), "name")) {
ret <- eval(parse(text = sprintf("RxODE:::rxParseErr(%s, addProp=\"%s\")", deparse1(x), addProp)))
assignInMyNamespace("rxErrEnv.diag.est", NULL)
assignInMyNamespace("rxErrEnv.theta", 1)
assignInMyNamespace("rxErrEnv.ret", "rx_r_")
assignInMyNamespace("rxErrEnv.init", NULL)
return(ret)
} else {
ret <- NULL
if (is(x, "character")) {
ret <- eval(parse(text = sprintf("RxODE:::rxParseErr(quote({%s}))", paste(x, collapse = "\n"))))
ret <- substring(ret, 3, nchar(ret) - 2)
} else {
x <- .convStr(x)
ret <- eval(x, rxErrEnv(x))
}
attr(ret, "ini") <- rxErrEnv.diag.est
assignInMyNamespace("rxErrEnv.diag.est", NULL)
assignInMyNamespace("rxErrEnv.theta", 1)
assignInMyNamespace("rxErrEnv.ret", "rx_r_")
assignInMyNamespace("rxErrEnv.init", NULL)
return(ret)
}
}
rxSimpleExprP <- function(x) {
if (is.name(x) || is.atomic(x)) {
return(TRUE)
} else {
return(FALSE)
}
}
#' This function splits a function based on + or - terms
#'
#' It uses the parser and does not disturb terms within other
#' functions. For example:
#'
#' a*exp(b+c)+d*log(e-f)-g*f
#'
#' would return
#'
#' c("a * exp(b + c)", "d * log(e - f)", "- g * f")
#'
#' @param x Quoted R expression for splitting
#' @param level Internal level of parsing
#' @param mult boolean to split based on * and / expressions instead.
#' By default this is turned off.
#' @return character vector of the split expressions
#' @author Matthew L. Fidler
#' @export
#' @keywords internal
rxSplitPlusQ <- function(x, level = 0, mult = FALSE) {
if (is(x, "character") && level == 0) {
return(eval(parse(text = sprintf("rxSplitPlusQ(quote(%s))", x))))
}
if (is.name(x) || is.atomic(x)) {
if (level == 0) {
return(paste(deparse1(x), collapse = ""))
} else {
return(character())
}
} else if (is.call(x)) { # Call recurse_call recursively
if ((mult && ((identical(x[[1]], quote(`*`)) ||
identical(x[[1]], quote(`/`))) && level == 0)) ||
(!mult && ((identical(x[[1]], quote(`+`)) ||
identical(x[[1]], quote(`-`))) && level == 0))) {
if (length(x) == 3) {
if (identical(x[[1]], quote(`+`))) {
one <- paste(deparse1(x[[3]]), collapse = "")
} else if (!mult) {
one <- paste("-", paste(deparse1(x[[3]]), collapse = ""))
} else if (identical(x[[1]], quote(`*`))) {
one <- paste(deparse1(x[[3]]), collapse = "")
} else if (mult) {
one <- paste("1/", paste(deparse1(x[[3]]), collapse = ""))
}
tmp <- rxSplitPlusQ(x[[2]], level = 0, mult = mult)
if (length(tmp) > 0) {
return(c(tmp, one))
} else {
tmp <- paste(deparse1(x[[2]]), collapse = "")
return(c(tmp, one))
}
} else {
## Unary + or -
if (identical(x[[1]], quote(`+`))) {
one <- paste(deparse1(x[[2]]), collapse = "")
} else {
one <- paste("-", paste(deparse1(x[[2]]), collapse = ""))
}
return(one)
}
} else {
tmp <- unlist(lapply(x, rxSplitPlusQ, level = 1, mult = mult))
if (level == 0) {
if (length(tmp) == 0) {
tmp <- paste(deparse1(x), collapse = "")
}
}
return(tmp)
}
} else if (is.pairlist(x)) {
## Call recurse_call recursively
tmp <- unlist(lapply(x, rxSplitPlusQ, level = level, mult = mult))
if (level == 0) {
if (length(tmp) == 0) {
tmp <- paste(deparse1(x), collapse = "")
}
}
return(tmp)
} else { # User supplied incorrect input
stop("Don't know how to handle type '", typeof(x), "'.",
call. = FALSE
)
}
}
.rxSupportedFunsExtra <- FALSE
.rxSupportedFuns <- function(extra = .rxSupportedFunsExtra) {
.ret <- c(
names(.rxSEsingle), names(.rxSEdouble), names(.rxSEeq),
"linCmt", names(.rxOnly), ls(.symengineFs)
)
if (extra) {
.ret <- c(.ret, c(
"rxEq", "rxNeq", "rxGeq", "rxLeq", "rxLt",
"rxGt", "rxAnd", "rxOr", "rxNot", "dabs", "dabs2", "abs0",
"dabs1", "abs1"
))
}
.ret
}
#' Get list of supported functions
#'
#' @return list of supported functions in RxODE
#' @examples
#' rxSupportedFuns()
#' @export
rxSupportedFuns <- function() {
.rxSupportedFuns(FALSE)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.