# R/nlmixrGrad.R In nlmixr: Nonlinear Mixed Effects Models in Population PK/PD

#### Documented in nlmixrGill83nlmixrGradFunnlmixrHess

##' Get the optimal forward difference interval by Gill83 method
##'
##' @inheritParams base::do.call
##' @param which Which parameters to calculate the forward difference
##'     and optimal forward difference interval
##' @inheritParams foceiControl
##' @author Matthew Fidler
##'
##' @return
##'
##' A data frame with the following columns:
##' \itemize{
##' \item{info}{Gradient evaluation/forward difference information}
##' \item{hf}{Forward difference final estimate}
##' \item{df}{Derivative estimate}
##' \item{df2}{2nd Derivative Estimate}
##' \item{err}{Error of the final estimate derivative}
##' \item{aEps}{Absolute difference for forward numerical differences}
##' \item{rEps}{Relative Difference for backward numerical differences}
##' \item{aEpsC}{Absolute difference for central numerical differences}
##' \item{rEpsC}{Relative difference for central numerical differences}
##' }
##'
##' The \code{info} returns one of the following:
##' \itemize{
##' \item{Not Assessed}{Gradient wasn't assessed}
##' \item{Good}{Success in Estimating optimal forward difference interval}
##' \item{High Grad Error}{Large error; Derivative estimate error \code{fTol} or more of the derivative}
##' \item{Constant Grad}{Function constant or nearly constant for this parameter}
##' \item{Odd/Linear Grad}{Function odd or nearly linear, df = K, df2 ~ 0}
##' \item{Grad changes quickly}{df2 increases rapidly as h decreases}
##' }
##' @examples
##'
##' ## These are taken from the numDeriv::grad examples to show how
##' ## simple gradients are assessed with nlmixrGill83
##'
##' nlmixrGill83(sin, pi)
##'
##' nlmixrGill83(sin, (0:10)*2*pi/10)
##'
##' func0 <- function(x){ sum(sin(x))  }
##' nlmixrGill83(func0 , (0:10)*2*pi/10)
##'
##' func1 <- function(x){ sin(10*x) - exp(-x) }
##' curve(func1,from=0,to=5)
##'
##' x <- 2.04
##' numd1 <- nlmixrGill83(func1, x)
##' exact <- 10*cos(10*x) + exp(-x)
##' c(numd1$df, exact, (numd1$df - exact)/exact)
##'
##' x <- c(1:10)
##' numd1 <- nlmixrGill83(func1, x)
##' exact <- 10*cos(10*x) + exp(-x)
##' cbind(numd1=numd1$df, exact, err=(numd1$df - exact)/exact)
##'
##' sc2.f <- function(x){
##'   n <- length(x)
##'    sum((1:n) * (exp(x) - x)) / n
##' }
##'
##' sc2.g <- function(x){
##'   n <- length(x)
##'   (1:n) * (exp(x) - 1) / n
##' }
##'
##' x0 <- rnorm(100)
##' exact <- sc2.g(x0)
##'
##' g <- nlmixrGill83(sc2.f, x0)
##'
##' max(abs(exact - g$df)/(1 + abs(exact))) ##' ##' @export nlmixrGill83 <- function(what, args, envir = parent.frame(), which, gillRtol = sqrt(.Machine$double.eps), gillK = 10L, gillStep = 2, gillFtol = 0) {
if (missing(which)) {
which <- rep(TRUE, length(args))
}
return(nlmixrGill83_(what, args, envir, which,
gillRtol = sqrt(.Machine$double.eps), gillK = 10L, gillStep = 2, gillFtol = 0 )) } .nlmixrGradInfo <- new.env(parent = emptyenv()) ##' Create a gradient function based on gill numerical differences ##' ##' @param thetaNames Names for the theta parameters ##' @inheritParams nlmixrGill83 ##' @inheritParams foceiControl ##' @param theta for the internal functions theta is the parameter ##' values ##' @param md5 the md5 identifier for the internal gradient function ##' information. ##' ##' @return A list with eval, grad, hist and unscaled ##' functions. This is an internal module used with dynmodel ##' ##' @keywords internal ##' ##' @examples ##' ##' func0 <- function(x){ sum(sin(x)) } ##' ##' ## This will printout every interation or when print=X ##' gf <- nlmixrGradFun(func0) ##' ##' ## x ##' x <- (0:10)*2*pi/10; ##' gf$eval(x)
##' gf$grad(x) ##' ##' ## x2 ##' x2 <- x+0.1 ##' gf$eval(x2)
##' gf$grad(x2) ##' ##' ## Gives the parameter history as a data frame ##' gf$hist()
##'
##' @export
nlmixrGradFun <- function(what, envir = parent.frame(), which, thetaNames,
gillRtol = sqrt(.Machine\$double.eps), gillK = 10L, gillStep = 2, gillFtol = 0,
useColor = crayon::has_color(),
printNcol = floor((getOption("width") - 23) / 12),
print = 1) {
.md5 <- digest::digest(list(what, gillRtol, gillK, gillStep, gillFtol))
.nlmixrGradInfo[["isRstudio"]] <- (Sys.getenv("RSTUDIO") == "1")
if (!missing(which)) {
.nlmixrGradInfo[[paste0(.md5, ".w")]] <- which
}
if (!missing(thetaNames)) {
}
.nlmixrGradInfo[[paste0(.md5, ".n")]] <- 0L
.nlmixrGradInfo[[paste0(.md5, ".f")]] <- what
.nlmixrGradInfo[[paste0(.md5, ".e")]] <- envir
.nlmixrGradInfo[[paste0(.md5, ".rtol")]] <- gillRtol
.nlmixrGradInfo[[paste0(.md5, ".k")]] <- gillK
.nlmixrGradInfo[[paste0(.md5, ".s")]] <- gillStep
.nlmixrGradInfo[[paste0(.md5, ".ftol")]] <- gillFtol
.eval <- eval(parse(text = paste0("function(theta){
nlmixrEval_(theta, \"", .md5, "\");
}")))
.grad <- eval(parse(text = paste0("function(theta){
nlmixrGrad_(theta, \"", .md5, "\");
}")))
.hist <- eval(parse(text = paste0("function(){
nlmixrParHist_(md5=\"", .md5, "\");
}")))
.unscaled <- eval(parse(text = paste0("function(theta){
nlmixrUnscaled_(theta,md5=\"", .md5, "\");
}")))
return(list(eval = .eval, grad = .grad, hist = .hist, unscaled = .unscaled))
}

##' Calculate Hessian
##'
##' Unlike stats::optimHess which assumes the gradient is accurate,
##' nlmixrHess does not make as strong an assumption that the gradient
##' is accurate but takes more function evaluations to calculate the
##' Hessian.  In addition, this procedures optimizes the forward
##' difference interval by \code{\link{nlmixrGill83}}
##'
##' If you have an analytical gradient function, you should use
##' stats::optimHess
##'
##' @inheritParams stats::optimHess
##' @param ... Extra arguments sent to \code{\link{nlmixrGill83}}
##' @inheritParams base::do.call
##' @author Matthew Fidler
##' @return Hessian matrix based on Gill83
##' @export
##' @references
##'
##' \url{https://v8doc.sas.com/sashtml/ormp/chap5/sect28.htm}
##'
##' @examples
##'  func0 <- function(x){ sum(sin(x))  }
##'  x <- (0:10)*2*pi/10
##'  nlmixrHess(x, func0)
##'
##' fr <- function(x) {   ## Rosenbrock Banana function
##'     x1 <- x[1]
##'     x2 <- x[2]
##'     100 * (x2 - x1 * x1)^2 + (1 - x1)^2
##' }
##' grr <- function(x) { ## Gradient of 'fr'
##'     x1 <- x[1]
##'     x2 <- x[2]
##'     c(-400 * x1 * (x2 - x1 * x1) - 2 * (1 - x1),
##'        200 *      (x2 - x1 * x1))
##' }
##'
##' h1 <- optimHess(c(1.2,1.2), fr, grr)
##'
##' h2 <- optimHess(c(1.2,1.2), fr)
##'
##' ## in this case h3 is closer to h1 where the gradient is known
##'
##' h3 <- nlmixrHess(c(1.2,1.2), fr)
nlmixrHess <- function(par, fn, ..., envir = parent.frame()) {
.gill <- nlmixrGill83(fn, par, envir = envir, ...)
return(nlmixrHess_(par, fn, envir, .gill))
}


## Try the nlmixr package in your browser

Any scripts or data that you put into this service are public.

nlmixr documentation built on March 27, 2022, 5:05 p.m.