tests/testthat/test-nlmixrGrad.R

In nlmixr2est: Nonlinear Mixed Effects Models in Population PK/PD, Estimation Routines

test_that("a Gill83 search that ends without an accepted interval restores theta", {
  # x[2] is flat to roundoff, so its search runs out of iterations ("Constant
  # Grad").  It is searched first and df[1] depends on x[2], so an x[2] left at
  # its last probe shows up as a wrong derivative for x[1].
  f <- function(x) (x[1] - 1)^2 + 1 + 1e-9 * x[2]^2 + 1e4 * (x[1] - 1) * x[2]
  g <- nlmixr2Gill83(f, c(1, 2))
  expect_equal(as.character(g$info), c("Good", "Constant Grad"))
  expect_equal(g$df[1], 2e4, tolerance = 1e-7)
  # with the last parameter left out, the base objective was never evaluated
  g <- nlmixr2Gill83(f, c(1, 2), which = c(TRUE, FALSE))
  expect_equal(as.character(g$info), c("Good", "Not Assessed"))
  expect_equal(g$f[1], 1)
  expect_equal(g$df[1], 2e4, tolerance = 1e-7)
})

test_that("nlmixr2Gill83() and nlmixr2Hess() never modify the caller's vector", {
  p <- c(1, 2)
  p0 <- p + 0
  h <- nlmixr2Hess(p, function(x) (x[1] - 1)^2 + 1)
  expect_identical(p, p0)
  expect_equal(h, matrix(c(2, 0, 0, 0), 2), tolerance = 1e-6)
  # an error part-way through a search used to leave the vector at the probe
  f <- function(x) if (x[1] > 1.0001) stop("outside") else sum(x^2)
  expect_error(nlmixr2Gill83(f, p), "outside")
  expect_identical(p, p0)
  # every call gets its own vector: one the objective keeps is not changed by
  # the probes after it, so only the base evaluation was taken at p
  acc <- new.env(parent = emptyenv())
  acc$x <- list()
  f2 <- function(x) {
    acc$x[[length(acc$x) + 1L]] <- x
    sum((x - c(1, 2))^2) + 1
  }
  p2 <- c(1, 2)
  h <- nlmixr2Hess(p2, f2)
  expect_identical(p2, p0)
  expect_equal(h, diag(2, 2), tolerance = 1e-6)
  expect_identical(sum(vapply(acc$x, identical, logical(1), p0)), 1L)
})

test_that("nlmixr2Gill83() and nlmixr2Hess() use the gill* arguments", {
  g <- nlmixr2Gill83(
    sin,
    1,
    gillRtol = 1e-6,
    gillK = 4L,
    gillStep = 3,
    gillFtol = 1e-3
  )
  expect_equal(g$gillRtol, 1e-6)
  expect_equal(g$gillK, 4L)
  expect_equal(g$gillStep, 3)
  expect_equal(g$gillFtol, 1e-3)
  # an accepted interval is 2*sqrt(|f|*gillRtol/|f''|)
  expect_equal(as.character(g$info), "Good")
  expect_equal(g$hf, 2 * sqrt(abs(sin(1)) * 1e-6 / abs(g$df2)))
  expect_false(identical(nlmixr2Hess(1, sin, gillRtol = 1e-4), nlmixr2Hess(1, sin)))
  expect_error(nlmixr2Gill83(sin, 1, gillK = -1L), "gillK")
})

test_that("gillK = 0 determines no step size", {
  # the search would start from gillStep * 2 * (1 + |x|) * sqrt(epsA / (1 + |f|)),
  # epsA = |f| * gillRtol; that is what is reported, unsearched
  .f <- sin(1)
  .epsA <- abs(.f) * sqrt(.Machine$double.eps)
  .h0 <- 2 * (2 * (1 + 1) * sqrt(.epsA / (1 + abs(.f))))
  g <- nlmixr2Gill83(sin, 1, gillK = 0L)
  expect_identical(as.character(g$info), "Not Assessed")
  expect_identical(g$hf, .h0)
  expect_identical(g$hphi, .h0)
  expect_true(is.na(g$df) && is.na(g$df2) && is.na(g$err))
  # a function flat in x is not searched either
  g <- nlmixr2Gill83(function(x) 1 + 1e-9 * x^2, 2, gillK = 0L)
  expect_identical(as.character(g$info), "Not Assessed")
  # nlmixr2Hess differences with that interval: |x| * rEpsC + aEpsC = .h0
  .h <- 1 * (.h0 * 0.5) + .h0 * 0.5
  .H <- (-sin(1 + 2 * .h) + 16 * sin(1 + .h) - 30 * .f + 16 * sin(1 - .h) - sin(1 - 2 * .h)) /
    (12 * .h * .h)
  expect_equal(nlmixr2Hess(1, sin, gillK = 0L), matrix(.H))
  # and so does the first nlmixr2GradFun gradient, a forward difference
  gf <- nlmixr2GradFun(sin, gillK = 0L, print = 0)
  gf$eval(1)
  expect_equal(gf$grad(1), (sin(1 + .h) - .f) / .h)
})

test_that("nlmixr2GradFun() gradients leave the point alone", {
  gf <- nlmixr2GradFun(
    function(x) if (x[1] > 1) NA_real_ else sum(x^2) + 1,
    print = 0
  )
  x <- c(1, 1)
  gf$eval(x)
  gf$grad(x) # the first gradient is the Gill search
  # x[1]'s forward leg is NA, so it takes the backward difference
  g1 <- gf$grad(x)
  g2 <- gf$grad(x)
  expect_identical(x, c(1, 1))
  expect_identical(g1, g2)
  expect_equal(g1[2], 2, tolerance = 1e-3)
})

test_that("nlmixr2GradFun() prints wrapped rows, underlining the last continuation", {
  withr::local_envvar(RSTUDIO = NA)
  nmSuppressMsg() # a fit run under suppressMessages() leaves the C-level print silenced
  gf <- nlmixr2GradFun(
    function(x) sum((x - c(1, 2, 3))^2 * c(1, 10, 100)),
    thetaNames = c("t1", "t2", "t3"),
    useColor = TRUE,
    printNcol = 2
  )
  expect_identical(
    capture.output(invisible(gf$eval(c(0.5, 0.5, 0.5)))),
    c(
      "|-----+---------------+-----------+-----------|",
      "|    #| Objective Fun |        t1 |        t2 |",
      "\033[4m|.....................|        t3 |...........|\033[0m",
      "|    1|     647.75000 |    0.5000 |    0.5000 |",
      "|.....................|    0.5000 |...........|\033[0m"
    )
  )
})

test_that("nlmixr2Hess() differences the parameters `which` leaves out of the search (issue 1140)", {
  # a quadratic: the 5-point stencil is exact at any step
  f <- function(x) (x[1] - 1)^2 + 3 * (x[2] - 2)^2 + x[1] * x[2] + 5
  h <- matrix(c(2, 1, 1, 6), 2)
  p <- c(0.5, 1)
  expect_equal(nlmixr2Hess(p, f), h, tolerance = 1e-6)
  # a parameter left out of the search is still differenced, about f(p)
  expect_equal(nlmixr2Hess(p, f, which = c(FALSE, TRUE)), h, tolerance = 1e-6)
  expect_equal(nlmixr2Hess(p, f, which = c(TRUE, FALSE)), h, tolerance = 1e-6)
  expect_equal(nlmixr2Hess(p, f, which = c(FALSE, FALSE)), h, tolerance = 1e-6)
  # the base objective is reported on every row, the left-out ones too
  expect_identical(nlmixr2Gill83(f, p, which = c(FALSE, TRUE))$f, rep(f(p), 2))
  expect_identical(nlmixr2Gill83(f, p, which = c(FALSE, FALSE))$f, rep(NA_real_, 2))
  # a left-out parameter is differenced with the interval gillK = 0 gives it
  f2 <- function(x) exp(x[1]) + sin(x[2]) + x[1] * x[2]
  expect_identical(
    nlmixr2Hess(p, f2, which = c(FALSE, TRUE))[1, 1],
    nlmixr2Hess(p, f2, gillK = 0L)[1, 1]
  )
  expect_identical(
    nlmixr2Hess(p, f2, which = c(FALSE, TRUE))[2, 2],
    nlmixr2Hess(p, f2)[2, 2]
  )
})

test_that("the objective sees the caller's names and attributes", {
  f <- function(x) unname((x["a"] - 1)^2 + x["b"]^2)
  gf <- nlmixr2GradFun(f, print = 0)
  x <- c(a = 2, b = 1)
  expect_equal(gf$eval(x), 2)
  # the first gradient is the Gill search
  expect_equal(unname(gf$grad(x)), c(2, 2), tolerance = 1e-3)
  expect_equal(unname(gf$grad(x)), c(2, 2), tolerance = 1e-3)
  h <- nlmixr2Hess(c(a = 2, b = 1), function(x) unname(x["a"]^2 + 3 * x["b"]^2))
  expect_equal(unname(h), diag(c(2, 6)), tolerance = 1e-3)
  fw <- function(x) sum(x^2) * attr(x, "weight")
  gw <- nlmixr2GradFun(fw, print = 0)
  xw <- structure(c(1, 2), weight = 3)
  expect_equal(gw$eval(xw), 15)
  expect_equal(gw$grad(xw), c(6, 12), tolerance = 1e-3)
  expect_equal(nlmixr2Hess(xw, fw), diag(c(6, 6)), tolerance = 1e-3)
})

test_that("nlmixr2GradFun() takes a Gill step factor of at least 1, as foceiControl() does", {
  f <- function(x) sum(x^2)
  for (.bad in list(0.5, 0, Inf, NA_real_, c(2, 3))) {
    expect_error(nlmixr2GradFun(f, gillStep = .bad), "gillStep")
  }
  expect_error(nlmixr2GradFun(f, gillStep = 0.5), "Element 1 is not >= 1", fixed = TRUE)
})

test_that("the gillFtol fallback of a high-error search keeps a step, not a slope (#1194)", {
  # slope ~0 at 0.5 fails the 50% error test; every trial has |phif| < gillFtol
  f <- function(x) 1 + (x - 0.5)^2 + (x - 0.5)^4
  g0 <- nlmixr2Gill83(f, 0.5)
  expect_identical(as.character(g0$info), "High Grad Error")
  g <- nlmixr2Gill83(f, 0.5, gillFtol = 10)
  expect_identical(as.character(g$info), "High Grad Error")
  # the step was a central-difference slope (about 0) before the fix
  expect_gt(g$hphi, 0)
  expect_identical(g$hphi, g$hf)
  expect_false(g$hf == g0$hf)
  # slope, curvature and error all belong to the trial at that step
  .h <- g$hf
  .fp <- f(0.5 + .h)
  .fn <- f(0.5 - .h)
  .f <- f(0.5)
  expect_equal(g$df, (.fp - .f) / .h, tolerance = 1e-12)
  expect_equal(g$df2, (.fp - 2 * .f + .fn) / (.h * .h), tolerance = 1e-12)
  expect_equal(g$err, .h * abs(g$df2) / 2 + 2 * abs(.f) * g$gillRtol / .h, tolerance = 1e-12)
})

test_that("the gillFtol fallback skips a step rejected for roundoff (#1194)", {
  # the first shrink of h overshoots into roundoff (Ch > 0.1), so the search
  # backs off to h0; the shrunk step used to be the fallback
  f <- function(x) 1 + 1e4 * x^2
  g0 <- nlmixr2Gill83(f, 0, gillStep = 50)
  g <- nlmixr2Gill83(f, 0, gillStep = 50, gillFtol = 1e10)
  expect_identical(as.character(g$info), "High Grad Error")
  expect_identical(g$hf, g0$hphi)
  expect_identical(g$hphi, g0$hphi)
})

test_that("the gillFtol fallback skips the growing steps rejected for roundoff (#1194)", {
  # h grows from roundoff (Ch > 0.1); only those rejected steps have
  # |phif| < gillFtol, so there is no fallback trial
  f <- function(x) 1 + 0.01 * x^2
  g0 <- nlmixr2Gill83(f, 0)
  expect_identical(as.character(g0$info), "High Grad Error")
  g <- nlmixr2Gill83(f, 0, gillFtol = 0.01 * g0$hphi * 0.6)
  expect_identical(g$hf, g0$hf)
  expect_identical(g$hphi, g0$hphi)
  expect_identical(g$err, g0$err)
})

test_that("a constant-gradient Gill83 search reports its roundoff error (#1194)", {
  g <- nlmixr2Gill83(function(x) 1 + 1e-9 * x^2, 2)
  expect_identical(as.character(g$info), "Constant Grad")
  expect_equal(g$err, 2 * g$gillRtol / g$hf)
})

test_that("an accurate Gill83 derivative is accepted (#1194)", {
  # a negative slope always failed 'err <= 0.5 * df'
  g <- nlmixr2Gill83(function(x) 1 - x + x^2, 0)
  expect_identical(as.character(g$info), "Good")
  expect_equal(g$df, -1, tolerance = 1e-3)
  # an interval accepted after shrinking the step compared df with a central
  # difference that was never computed (0)
  g <- nlmixr2Gill83(function(x) 1 + x + 500 * x^2, 0)
  expect_identical(as.character(g$info), "Good")
  expect_equal(g$df, 1, tolerance = 1e-2)
})

test_that("a Gill83 search that backs off a shrunk step keeps that step's curvature (#1197)", {
  # h0 is accurate (Ch < 0.001) but the curvature collapses at h0/10, so the
  # search rejects the shrunk step for roundoff and backs off to h0
  x0 <- 0.1
  f <- function(x) 1 + x^2 + 1e7 * (x - x0)^4
  g <- nlmixr2Gill83(f, x0, gillStep = 10)
  expect_identical(as.character(g$info), "Good")
  h <- g$hphi
  phi <- (f(x0 + h) - 2 * f(x0) + f(x0 - h)) / h^2
  expect_equal(g$df2, phi, tolerance = 1e-6)
  expect_equal(g$hf, 2 * sqrt(abs(f(x0)) * g$gillRtol / abs(phi)), tolerance = 1e-6)
  # the rejected step is roundoff-dominated: its curvature fails the condition test
  hk <- h / 10
  phiK <- (f(x0 + hk) - 2 * f(x0) + f(x0 - hk)) / hk^2
  expect_gt(4 * abs(f(x0)) * g$gillRtol / (hk^2 * abs(phiK)), 0.1)
})