diff --git a/NAMESPACE b/NAMESPACE index 98d051b15..373cc51ed 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -284,6 +284,7 @@ export(.rxSEeqUsr) export(.rxSens) export(.rxToSE) export(.rxTransInfo) +export(.rxTransform) export(.rxWithOptions) export(.rxWithSink) export(.rxWithSinkBoth) @@ -331,6 +332,7 @@ export(assertVariableExists) export(assertVariableName) export(assertVariableNew) export(binomProbs) +export(boxCoxInv) export(cvPost) export(dfWishart) export(erf) @@ -584,6 +586,8 @@ export(warnRxBounded) export(write.template.server) export(write.template.ui) export(xlab) +export(yeoJohnson) +export(yeoJohnsonInv) export(ylab) export(zeroRe) import(data.table) diff --git a/R/utils.R b/R/utils.R index d1de668ed..6f7d9661c 100644 --- a/R/utils.R +++ b/R/utils.R @@ -386,6 +386,67 @@ gammapInv <- function(a, p) { gammapInva <- function(x, p) { .Call(`_gammapInva`, x, p, PACKAGE = "rxode2") } +#' rxode2 general transformation function +#' +#' @param x value that will be transformed +#' @param lambda lambda value for the transformation +#' @param transform transformation to use (can be integer or string +#' matching supported transformations) +#' @param low lower bound for the transformation +#' @param high upper bound for the transformation +#' @param inverse boolean if the inverse transformation should be performed +#' @return transformed value +#' @export +#' @author Matthew L. Fidler +#' @keywords internal +#' @examples +#' +#' logit(0.25) +#' +#' .rxTransform(0.25, transform="logit") +#' +#' expit(-1.09) +#' +#' .rxTransform(-1.09, transform="logit", inverse=TRUE) +#' +.rxTransform <- function(x, lambda=1.0, + low = 0.0, high = 1.0, + transform=c("boxCox", "yeoJohnson", "untransformed", + "lnorm", "logit", "logit + yeoJohnson", + "probit", "probit + yeoJohnson", + "logit + boxCox", "probit + boxCox"), + inverse=FALSE) { + if (is.integer(transform)) { + } else { + transform <- factor(match.arg(transform), + levels=c("boxCox", "yeoJohnson", "untransformed", + "lnorm", "logit", "logit + yeoJohnson", + "probit", "probit + yeoJohnson", "logit + boxCox", + "probit + boxCox")) + transform <- as.integer(transform)-1L + } + if (length(lambda) > 1 || + length(low) > 1 || + length(high) > 1 || + length(transform) > 1 || + length(inverse) > 1) { + .df <- data.frame(x = x, lambda = lambda, low = low, high = high, + transform=transform, inverse=inverse) + vapply(1:nrow(.df), + function(i) { + .powerD(.df$x[i], .df$lambda[i], .df$low[i], .df$high[i], + .df$transform[i], .df$inverse[i]) + }, numeric(1), USE.NAMES = FALSE) + } else { + checkmate::assertNumeric(x, any.missing = FALSE) + checkmate::assertNumeric(lambda, any.missing = FALSE) + checkmate::assertNumeric(low, any.missing = FALSE) + checkmate::assertNumeric(high, any.missing = FALSE) + checkmate::assertInteger(transform, any.missing = FALSE) + checkmate::assertLogical(inverse, any.missing = FALSE) + .Call(`_rxode2_powerD`, x, low, high, lambda, transform, inverse) + } +} #' logit and inverse logit (expit) functions #' @@ -435,22 +496,27 @@ gammapInva <- function(x, p) { #' logitNormInfo(logit(0.25), sd = 0.1) #' #' logitNormInfo(logit(1, 0, 10), sd = 1, low = 0, high = 10) +#' #' @export logit <- function(x, low = 0, high = 1) { - .Call(`_logit`, x, low, high, PACKAGE = "rxode2") + .rxTransform(x, 1.0, low, high, 4L, FALSE) } #' @rdname logit #' @export expit <- function(alpha, low = 0, high = 1) { - .Call(`_expit`, alpha, low, high, PACKAGE = "rxode2") + .rxTransform(alpha, 1.0, low, high, 4L, TRUE) } #' @rdname logit #' @export logitNormInfo <- function(mean = 0, sd = 1, low = 0, high = 1, abs.tol = 1e-6, ...) { - .fM1 <- function(x) .Call(`_expit`, x, low, high, PACKAGE = "rxode2") * dnorm(x, mean = mean, sd = sd) + .fM1 <- function(x) { + expit(x, low, high) * dnorm(x, mean = mean, sd = sd) + } .m <- integrate(.fM1, -Inf, Inf, abs.tol = abs.tol, ...)$value - .fV <- function(x) (.Call(`_expit`, x, low, high, PACKAGE = "rxode2") - .m)^2 * dnorm(x, mean = mean, sd = sd) + .fV <- function(x){ + (expit(x, low, high) - .m)^2 * dnorm(x, mean = mean, sd = sd) + } .v <- integrate(.fV, -Inf, Inf, abs.tol = abs.tol, ...)$value c(mean = .m, var = .v, cv = sqrt(.v) / .m) } @@ -470,26 +536,60 @@ logitNormInfo <- function(mean = 0, sd = 1, low = 0, high = 1, abs.tol = 1e-6, . #' probitNormInfo(probit(1, 0, 10), sd = 1, low = 0, high = 10) #' @export probit <- function(x, low = 0, high = 1) { - .Call(`_probit`, x, low, high, PACKAGE = "rxode2") + .rxTransform(x, 1.0, low, high, 6L, FALSE) } #' @rdname probit #' @export probitInv <- function(x, low = 0, high = 1) { - .Call(`_probitInv`, x, low, high, PACKAGE = "rxode2") + .rxTransform(x, 1.0, low, high, 6L, TRUE) } #' @rdname logit #' @export probitNormInfo <- function(mean = 0, sd = 1, low = 0, high = 1, abs.tol = 1e-6, ...) { - .fM1 <- function(x) .Call(`_probitInv`, x, low, high, PACKAGE = "rxode2") * dnorm(x, mean = mean, sd = sd) + .fM1 <- function(x) probitInv(x, low, high) * dnorm(x, mean = mean, sd = sd) .m <- integrate(.fM1, -Inf, Inf, abs.tol = abs.tol, ...)$value - .fV <- function(x) (.Call(`_probitInv`, x, low, high, PACKAGE = "rxode2") - .m)^2 * dnorm(x, mean = mean, sd = sd) + .fV <- function(x) (probitInv(x, low, high) - .m)^2 * dnorm(x, mean = mean, sd = sd) .v <- integrate(.fV, -Inf, Inf, abs.tol = abs.tol, ...)$value c(mean = .m, var = .v, cv = sqrt(.v) / .m) } +#' boxCox/yeoJohnson and inverse boxCox/yeoJohnson functions +#' +#' @param x input value(s) to transform +#' @param lambda lambda value for the transformation +#' @return values from boxCox and boxCoxInv +#' @export +#' +#' boxCox(10, 0.5) +#' +#' boxCoxInv(4.32, 0.5) +#' +#' yeoJohson(10, 0.5) +#' +#' yeoJohnsonInv(4.32, 0.5) +#' +boxCox <- function(x, lambda = 1.0) { + checkmate::assertNumeric(x, lower=0.0, any.missing=FALSE) + .rxTransform(x, lambda, low=0.0, high=1.0, 0L, FALSE) +} +#' @rdname boxCox +#' @export +boxCoxInv <- function(x, lambda = 1.0) { + .rxTransform(x, lambda, low=0.0, high=1.0, 0L, TRUE) +} +#' @rdname boxCox +#' @export +yeoJohnson <- function(x, lambda = 1.0) { + .rxTransform(x, lambda, low=0.0, high=1.0, 1L, FALSE) +} +#' @rdname boxCox +#' @export +yeoJohnsonInv <- function(x, lambda = 1.0) { + .rxTransform(x, lambda, low=0.0, high=1.0, 1L, TRUE) +} #' Get/Set the number of threads that rxode2 uses #' #' @param threads NULL (default) rereads environment variables. 0 diff --git a/man/boxCox.Rd b/man/boxCox.Rd new file mode 100644 index 000000000..9fe20f297 --- /dev/null +++ b/man/boxCox.Rd @@ -0,0 +1,28 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/utils.R +\name{boxCox} +\alias{boxCox} +\alias{boxCoxInv} +\alias{yeoJohnson} +\alias{yeoJohnsonInv} +\title{boxCox/yeoJohnson and inverse boxCox/yeoJohnson functions} +\usage{ +boxCox(x, lambda = 1) + +boxCoxInv(x, lambda = 1) + +yeoJohnson(x, lambda = 1) + +yeoJohnsonInv(x, lambda = 1) +} +\arguments{ +\item{x}{input value(s) to transform} + +\item{lambda}{lambda value for the transformation} +} +\value{ +values from boxCox and boxCoxInv +} +\description{ +boxCox/yeoJohnson and inverse boxCox/yeoJohnson functions +} diff --git a/man/dot-rxTransform.Rd b/man/dot-rxTransform.Rd new file mode 100644 index 000000000..6dd0db20b --- /dev/null +++ b/man/dot-rxTransform.Rd @@ -0,0 +1,52 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/utils.R +\name{.rxTransform} +\alias{.rxTransform} +\title{rxode2 general transformation function} +\usage{ +.rxTransform( + x, + lambda = 1, + low = 0, + high = 1, + transform = c("boxCox", "yeoJohnson", "untransformed", "lnorm", "logit", + "logit + yeoJohnson", "probit", "probit + yeoJohnson", "logit + boxCox", + "probit + boxCox"), + inverse = FALSE +) +} +\arguments{ +\item{x}{value that will be transformed} + +\item{lambda}{lambda value for the transformation} + +\item{low}{lower bound for the transformation} + +\item{high}{upper bound for the transformation} + +\item{transform}{transformation to use (can be integer or string +matching supported transformations)} + +\item{inverse}{boolean if the inverse transformation should be performed} +} +\value{ +transformed value +} +\description{ +rxode2 general transformation function +} +\examples{ + +logit(0.25) + +.rxTransform(0.25, transform="logit") + +expit(-1.09) + +.rxTransform(-1.09, transform="logit", inverse=TRUE) + +} +\author{ +Matthew L. Fidler +} +\keyword{internal} diff --git a/man/logit.Rd b/man/logit.Rd index 580f392cf..a90eb22ad 100644 --- a/man/logit.Rd +++ b/man/logit.Rd @@ -65,4 +65,5 @@ expit(-1.09) logitNormInfo(logit(0.25), sd = 0.1) logitNormInfo(logit(1, 0, 10), sd = 1, low = 0, high = 10) + } diff --git a/man/rxode2.Rd b/man/rxode2.Rd index 3574082d4..1f11b5643 100644 --- a/man/rxode2.Rd +++ b/man/rxode2.Rd @@ -215,6 +215,148 @@ suggest the package emphasis on pharmacometrics modeling, including pharmacokinetics (PK), pharmacodynamics (PD), disease progression, drug-disease modeling, etc. +The ODE-based model specification may be coded inside four places: +\itemize{ +\item Inside a \code{rxode2({})} block statements: +}\if{html}{\out{ + +}} + + +\if{html}{\out{
}}\preformatted{library(rxode2) +mod <- rxode2(\{ + # simple assignment + C2 <- centr/V2 + + # time-derivative assignment + d/dt(centr) <- F*KA*depot - CL*C2 - Q*C2 + Q*C3; +\}) +}\if{html}{\out{
}} + +\if{html}{\out{
}}\preformatted{## using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’ +}\if{html}{\out{
}} +\itemize{ +\item Inside a \code{rxode2("")} string statement: +}\if{html}{\out{ + +}} + + +\if{html}{\out{
}}\preformatted{mod <- rxode2(" + # simple assignment + C2 <- centr/V2 + + # time-derivative assignment + d/dt(centr) <- F*KA*depot - CL*C2 - Q*C2 + Q*C3; +") +}\if{html}{\out{
}} + +\if{html}{\out{
}}\preformatted{## using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’ +}\if{html}{\out{
}} +\itemize{ +\item In a file name to be loaded by rxode2: +}\if{html}{\out{ + +}} + + +\if{html}{\out{
}}\preformatted{writeLines(" + # simple assignment + C2 <- centr/V2 + + # time-derivative assignment + d/dt(centr) <- F*KA*depot - CL*C2 - Q*C2 + Q*C3; +", "modelFile.rxode2") +mod <- rxode2(filename='modelFile.rxode2') +unlink("modelFile.rxode2") +}\if{html}{\out{
}} +\itemize{ +\item In a model function which can be parsed by \code{rxode2}: +}\if{html}{\out{ + +}} + + +\if{html}{\out{
}}\preformatted{mod <- function() \{ + model(\{ + # simple assignment + C2 <- centr/V2 + + # time-derivative assignment + d/dt(centr) <- F*KA*depot - CL*C2 - Q*C2 + Q*C3; + \}) +\} + +mod <- rxode2(mod) # or simply mod() if the model is at the end of the function + +# These model functions often have residual components and initial +# (`ini(\{\})`) conditions attached as well. For example the +# theophylline model can be written as: + +one.compartment <- function() \{ + ini(\{ + tka <- 0.45 # Log Ka + tcl <- 1 # Log Cl + tv <- 3.45 # Log V + eta.ka ~ 0.6 + eta.cl ~ 0.3 + eta.v ~ 0.1 + add.sd <- 0.7 + \}) + model(\{ + ka <- exp(tka + eta.ka) + cl <- exp(tcl + eta.cl) + v <- exp(tv + eta.v) + d/dt(depot) = -ka * depot + d/dt(center) = ka * depot - cl / v * center + cp = center / v + cp ~ add(add.sd) + \}) +\} + +# after parsing the model +mod <- one.compartment() +}\if{html}{\out{
}} + +For the block statement, character string or text file an internal +\code{rxode2} compilation manager translates the ODE system into C, compiles +it and loads it into the R session. The call to \code{rxode2} produces an +object of class \code{rxode2} which consists of a list-like structure +(environment) with various member functions. + +For the last type of model (a model function), a call to \code{rxode2} +creates a parsed \code{rxode2} ui that can be translated to the \code{rxode2} +compilation model. + +\if{html}{\out{
}}\preformatted{mod$simulationModel +}\if{html}{\out{
}} + +\if{html}{\out{
}}\preformatted{## using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’ + +## rxode2 3.0.0 model named rx_ee5e0e3c742f5a87d21a1eef8d862958 model (ready). +## x$state: depot, center +## x$stateExtra: cp +## x$params: tka, tcl, tv, add.sd, eta.ka, eta.cl, eta.v, rxerr.cp +## x$lhs: ka, cl, v, cp, ipredSim, sim +}\if{html}{\out{
}} + +\if{html}{\out{
}}\preformatted{# or +mod$simulationIniModel +}\if{html}{\out{
}} + +\if{html}{\out{
}}\preformatted{## using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’ + +## rxode2 3.0.0 model named rx_67efed5f8cd3901f46cb7f70686d23cd model (ready). +## x$state: depot, center +## x$stateExtra: cp +## x$params: tka, tcl, tv, add.sd, eta.ka, eta.cl, eta.v, rxerr.cp +## x$lhs: ka, cl, v, cp, ipredSim, sim +}\if{html}{\out{
}} + +This is the same type of function required for \code{nlmixr2} estimation and +can be extended and modified by model piping. For this reason will be +focused on in the documentation. + This basic model specification consists of one or more statements optionally terminated by semi-colons \verb{;} and optional comments (comments are delimited by \verb{#} and an end-of-line). diff --git a/src/init.c b/src/init.c index dbadca008..8227ae7fe 100644 --- a/src/init.c +++ b/src/init.c @@ -109,8 +109,6 @@ SEXP _gammaqInva(SEXP, SEXP); double expit(double, double, double); double logit(double, double, double); -SEXP _expit(SEXP, SEXP, SEXP); -SEXP _logit(SEXP, SEXP, SEXP); SEXP _linCmtParse(SEXP vars, SEXP inStr, SEXP verbose); SEXP _rxode2_linCmtGen(SEXP linCmt, SEXP vars, SEXP linCmtSens, SEXP verbose); @@ -261,8 +259,6 @@ void sortIds(rx_solve* rx, int ini); void handleTlast(double *time, rx_solving_options_ind *ind); -SEXP _probit(SEXP xS, SEXP lowS, SEXP highS); -SEXP _probitInv(SEXP xS, SEXP lowS, SEXP highS); SEXP _rxode2_rxrandnV(SEXP, SEXP); SEXP _rxode2_rxErf(SEXP); @@ -539,8 +535,11 @@ SEXP _rxode2_rxode2Ptr(void) { return ret; } +SEXP _rxode2_powerD(SEXP, SEXP, SEXP, SEXP, SEXP, + SEXP); void R_init_rxode2(DllInfo *info){ R_CallMethodDef callMethods[] = { + {"_rxode2_powerD", (DL_FUNC) &_rxode2_powerD, 6}, {"_rxode2_rxode2Ptr", (DL_FUNC) &_rxode2_rxode2Ptr, 0}, {"_rxode2_iniDparserPtr", (DL_FUNC) &_rxode2_iniDparserPtr, 1}, {"_iniPreciseSumsPtr", (DL_FUNC) &iniPreciseSumsPtr, 1}, @@ -682,8 +681,6 @@ void R_init_rxode2(DllInfo *info){ {"_gammaqInv", (DL_FUNC) _gammaqInv, 2}, {"_gammapInva", (DL_FUNC) _gammapInv, 2}, {"_gammaqInva", (DL_FUNC) _gammaqInv, 2}, - {"_expit", (DL_FUNC) _expit, 3}, - {"_logit", (DL_FUNC) _logit, 3}, {"_calcDerived", (DL_FUNC) _calcDerived, 4}, {"_linCmtParse", (DL_FUNC) _linCmtParse, 3}, {"_rxode2_linCmtGen", (DL_FUNC) _rxode2_linCmtGen, 4}, @@ -692,8 +689,6 @@ void R_init_rxode2(DllInfo *info){ {"getRxThreads_R", (DL_FUNC) getRxThreads_R, 1}, {"setRxthreads", (DL_FUNC) setRxthreads, 3}, {"_rxHasOpenMp", (DL_FUNC) _rxHasOpenMp, 0}, - {"_probit", (DL_FUNC) _probit, 3}, - {"_probitInv", (DL_FUNC) _probitInv, 3}, {"_rxode2_isNullZero", (DL_FUNC) _rxode2_isNullZero, 1}, {"_rxode2_invWR1d", (DL_FUNC) _rxode2_invWR1d, 3}, {"_rxode2_rxSimThetaOmega", (DL_FUNC) _rxode2_rxSimThetaOmega, 28}, diff --git a/src/utilc.c b/src/utilc.c index 3b4a2c0a6..80851462b 100644 --- a/src/utilc.c +++ b/src/utilc.c @@ -596,117 +596,17 @@ double probitInv(double alpha, double low, double high) { return _powerDi(alpha, 1.0, 6, low, high); } -SEXP _probit(SEXP xS, SEXP lowS, SEXP highS) { +SEXP _rxode2_powerD(SEXP xS, SEXP lowS, SEXP highS, SEXP lambdaS, SEXP yjS, SEXP inverseS) { int typex = TYPEOF(xS); int typelow = TYPEOF(lowS); int typehigh = TYPEOF(highS); - double low, high; - if (Rf_length(lowS) != 1){ - Rf_errorcall(R_NilValue, _("'low' must be a numeric of length 1")); - } - if (Rf_length(highS) != 1){ - Rf_errorcall(R_NilValue, _("'high' must be a numeric of length 1")); - } - if (typelow == INTSXP){ - low = (double)(INTEGER(lowS)[0]); - } else if (typelow == REALSXP) { - low = REAL(lowS)[0]; - } else { - Rf_errorcall(R_NilValue, _("'low' must be a numeric of length 1")); - } - if (typehigh == INTSXP){ - high = (double)(INTEGER(highS)[0]); - } else if (typehigh == REALSXP) { - high = REAL(highS)[0]; - } else { - Rf_errorcall(R_NilValue, _("'high' must be a numeric of length 1")); - } - if (high <= low) { - Rf_errorcall(R_NilValue, _("'high' must be greater than 'low'")); - } - int lenx = Rf_length(xS); - double *xD = NULL; - int *xI = NULL; - int isD=0; - if (typex == REALSXP){ - isD=1; - xD = REAL(xS); - } else if (typex == INTSXP){ - xI = INTEGER(xS); - } - SEXP ret = PROTECT(Rf_allocVector(REALSXP, lenx)); - double *retD = REAL(ret); - if (isD){ - for (int i = lenx; i--;){ - retD[i] = probit(xD[i], low, high); - } - } else { - for (int i = lenx; i--;){ - retD[i] = probit((double)(xI[i]), low, high); - } - } - UNPROTECT(1); - return ret; -} - -SEXP _probitInv(SEXP xS, SEXP lowS, SEXP highS) { - int typex = TYPEOF(xS); - int typelow = TYPEOF(lowS); - int typehigh = TYPEOF(highS); - double low, high; - if (Rf_length(lowS) != 1){ - Rf_errorcall(R_NilValue, _("'low' must be a numeric of length 1")); - } - if (Rf_length(highS) != 1){ - Rf_errorcall(R_NilValue, _("'high' must be a numeric of length 1")); - } - if (typelow == INTSXP){ - low = (double)(INTEGER(lowS)[0]); - } else if (typelow == REALSXP) { - low = REAL(lowS)[0]; - } else { - Rf_errorcall(R_NilValue, _("'low' must be a numeric of length 1")); - } - if (typehigh == INTSXP){ - high = (double)(INTEGER(highS)[0]); - } else if (typehigh == REALSXP) { - high = REAL(highS)[0]; - } else { - Rf_errorcall(R_NilValue, _("'high' must be a numeric of length 1")); - } - if (high <= low) { - Rf_errorcall(R_NilValue, _("'high' must be greater than 'low'")); - } - int lenx = Rf_length(xS); - double *xD = NULL; - int *xI = NULL; - int isD=0; - if (typex == REALSXP){ - isD=1; - xD = REAL(xS); - } else if (typex == INTSXP){ - xI = INTEGER(xS); - } - SEXP ret = PROTECT(Rf_allocVector(REALSXP, lenx)); - double *retD = REAL(ret); - if (isD){ - for (int i = lenx; i--;){ - retD[i] = probitInv(xD[i], low, high); - } - } else { - for (int i = lenx; i--;){ - retD[i] = probitInv((double)(xI[i]), low, high); - } + int typelambda = TYPEOF(lambdaS); + int inverse = INTEGER(inverseS)[0]; + int yj = INTEGER(yjS)[0]; + double low, high, lambda; + if (Rf_length(lambdaS) != 1){ + Rf_errorcall(R_NilValue, _("'lambda' must be a numeric of length 1")); } - UNPROTECT(1); - return ret; -} - -SEXP _logit(SEXP xS, SEXP lowS, SEXP highS) { - int typex = TYPEOF(xS); - int typelow = TYPEOF(lowS); - int typehigh = TYPEOF(highS); - double low, high; if (Rf_length(lowS) != 1){ Rf_errorcall(R_NilValue, _("'low' must be a numeric of length 1")); } @@ -730,59 +630,12 @@ SEXP _logit(SEXP xS, SEXP lowS, SEXP highS) { if (high <= low) { Rf_errorcall(R_NilValue, _("'high' must be greater than 'low'")); } - int lenx = Rf_length(xS); - double *xD = NULL; - int *xI = NULL; - int isD=0; - if (typex == REALSXP){ - isD=1; - xD = REAL(xS); - } else if (typex == INTSXP){ - xI = INTEGER(xS); - } - SEXP ret = PROTECT(Rf_allocVector(REALSXP, lenx)); - double *retD = REAL(ret); - if (isD){ - for (int i = lenx; i--;){ - retD[i] = logit(xD[i], low, high); - } - } else { - for (int i = lenx; i--;){ - retD[i] = logit((double)(xI[i]), low, high); - } - } - UNPROTECT(1); - return ret; -} - - -SEXP _expit(SEXP xS, SEXP lowS, SEXP highS) { - int typex = TYPEOF(xS); - int typelow = TYPEOF(lowS); - int typehigh = TYPEOF(highS); - double low, high; - if (Rf_length(lowS) != 1){ - Rf_errorcall(R_NilValue, _("'low' must be a numeric of length 1")); - } - if (Rf_length(highS) != 1){ - Rf_errorcall(R_NilValue, _("'high' must be a numeric of length 1")); - } - if (typelow == INTSXP){ - low = (double)(INTEGER(lowS)[0]); - } else if (typelow == REALSXP) { - low = REAL(lowS)[0]; + if (typelambda == INTSXP) { + lambda = (double)(INTEGER(lambdaS)[0]); + } else if (typelambda == REALSXP) { + lambda = REAL(lambdaS)[0]; } else { - Rf_errorcall(R_NilValue, _("'low' must be a numeric of length 1")); - } - if (typehigh == INTSXP){ - high = (double)(INTEGER(highS)[0]); - } else if (typehigh == REALSXP) { - high = REAL(highS)[0]; - } else { - Rf_errorcall(R_NilValue, _("'high' must be a numeric of length 1")); - } - if (high <= low) { - Rf_errorcall(R_NilValue, _("'high' must be greater than 'low'")); + Rf_errorcall(R_NilValue, _("'lambda' must be a numeric of length 1")); } int lenx = Rf_length(xS); double *xD = NULL; @@ -796,13 +649,25 @@ SEXP _expit(SEXP xS, SEXP lowS, SEXP highS) { } SEXP ret = PROTECT(Rf_allocVector(REALSXP, lenx)); double *retD = REAL(ret); - if (isD){ - for (int i = lenx; i--;){ - retD[i] = expit(xD[i], low, high); + if (inverse) { + if (isD) { + for (int i = lenx; i--;){ + retD[i] = _powerDi(xD[i], lambda, yj, low, high); + } + } else { + for (int i = lenx; i--;){ + retD[i] = _powerDi((double)(xI[i]), lambda, yj, low, high); + } } } else { - for (int i = lenx; i--;){ - retD[i] = expit((double)(xI[i]), low, high); + if (isD) { + for (int i = lenx; i--;){ + retD[i] = _powerD(xD[i], lambda, yj, low, high); + } + } else { + for (int i = lenx; i--;){ + retD[i] = _powerD((double)(xI[i]), lambda, yj, low, high); + } } } UNPROTECT(1);