Made zero tolerance adaptive for higher orders of diff. and accuracy

Author: Fifis

R/gradient.R

@@ -36,7 +36,7 @@
 checkDimensions <- function(FUN, x, f0 = NULL, func = NULL,
                             elementwise = NA, vectorised = NA, multivalued = NA,
                             deriv.order = 1, acc.order = 2, side = 0, h = NULL,
-                            zero.tol = sqrt(.Machine$double.eps), cores = 1, preschedule = TRUE, cl = NULL, ...) {
+                            zero.tol = NULL, cores = 1, preschedule = TRUE, cl = NULL, ...) {
   if (missing(FUN)) {
     if (is.function(func)) {
       FUN <- func
@@ -273,8 +273,8 @@ generateGrid <- function(x, h, stencils, elementwise, vectorised) {
 #'   \code{"DV"}, or \code{"SW"} to be used in [gradstep()]). The default value is
 #'   described in \code{?GenD}.
 #' @param zero.tol Small positive integer: if \code{abs(x) >= zero.tol}, then, the automatically
-#'   guessed step size is relative (\code{x} multiplied by the step), unless an auto-selection
-#'   procedure is requested; otherwise, it is absolute.
+#'   guessed step size is relative or near-relative (\code{x} multiplied by the step),
+#'   unless an auto-selection procedure is requested; otherwise, it is absolute.
 #' @param f0 Optional numeric: if provided, used to determine the vectorisation type
 #'   to save time. If FUN(x) must be evaluated (e.g. second derivatives), saves one evaluation.
 #' @param h0 Numeric scalar of vector: initial step size for automatic search with
@@ -368,7 +368,7 @@ generateGrid <- function(x, h, stencils, elementwise, vectorised) {
 #' system.time(GenD(fslow, x, side = 1, acc.order = 2))
 GenD <- function(FUN, x, elementwise = NA, vectorised = NA, multivalued = NA,
                  deriv.order = 1L, side = 0, acc.order = 2L, stencil = NULL,
-                 h = NULL, zero.tol = sqrt(.Machine$double.eps),  h0 = NULL, control = list(),
+                 h = NULL, zero.tol = NULL,  h0 = NULL, control = list(),
                  f0 = NULL, cores = 1, preschedule = TRUE, cl = NULL,
                  func = NULL, method = NULL, method.args = list(), ...) {
   if (is.function(x) && !is.function(FUN)) {
@@ -629,7 +629,7 @@ GenD <- function(FUN, x, elementwise = NA, vectorised = NA, multivalued = NA,
 #' Grad(LReLU, seq(-1, 1, 0.1))
 Grad <- function(FUN, x, elementwise = NA, vectorised = NA, multivalued = NA,
                  deriv.order = 1L, side = 0, acc.order = 2, stencil = NULL,
-                 h = NULL, zero.tol = sqrt(.Machine$double.eps), h0 = NULL, control = list(),
+                 h = NULL, zero.tol = NULL, h0 = NULL, control = list(),
                  f0 = NULL, cores = 1, preschedule = TRUE, cl = NULL,
                  func = NULL, method = NULL, method.args = list(), ...) {
   if (is.function(x) && !is.function(FUN)) {
@@ -714,7 +714,7 @@ Grad <- function(FUN, x, elementwise = NA, vectorised = NA, multivalued = NA,
 #' @export
 Jacobian <- function(FUN, x, elementwise = NA, vectorised = NA, multivalued = NA,
                      deriv.order = 1L, side = 0, acc.order = 2, stencil = NULL,
-                     h = NULL, zero.tol = sqrt(.Machine$double.eps), h0 = NULL,
+                     h = NULL, zero.tol = NULL, h0 = NULL,
                      control = list(), f0 = NULL,
                      cores = 1, preschedule = TRUE, cl = NULL,
                      func = NULL, method = NULL, method.args = list(), ...) {

R/step-select.R

@@ -14,27 +14,32 @@
 #' # The step-selection function is piecewise linear in log-coordinates
 #' plot(-12:4, stepx(10^(-12:4)), log = "y", type = "l")
 #' stepx(10^(-10:2), deriv.order = 2, acc.order = 4)
-stepx <- function(x, deriv.order = 1, acc.order = 2, zero.tol = sqrt(.Machine$double.eps)) {
+stepx <- function(x, deriv.order = 1, acc.order = 2, zero.tol = NULL) {
   x <- abs(x)
   n <- length(x)
-  i1 <- x < zero.tol
-  i3 <- x > 1
-  i2 <- (!i1) & (!i3)
-  ret <- rep(zero.tol, length(x))
 
   if (length(deriv.order) == 1) deriv.order <- rep(deriv.order, n)
   if (length(acc.order) == 1) acc.order <- rep(acc.order, n)
   if (length(deriv.order) != n) stop("The argument 'deriv.order' must have length 1 or length(x).")
   if (length(acc.order) != n) stop("The argument 'acc.order' must have length 1 or length(x).")
-
   ad <- deriv.order + acc.order
+
+  if (is.null(zero.tol)) zero.tol <- .Machine$double.eps^(1/(ad-0.5))
+  if (length(zero.tol) == 1) zero.tol <- rep(zero.tol, n)  # When it is passed from outside
+  i1 <- x < zero.tol
+  i3 <- x > 1
+  i2 <- (!i1) & (!i3)
+  ret <- zero.tol  # Initialised with i1 case everywhere
+
   if (any(i3)) ret[i3] <- x[i3] * .Machine$double.eps^(1/ad[i3])
   if (any(i2)) {  # Exponential interpolation via log(y) ~ a + b*log(x)
-    # f(sqrt(macheps)) = macheps^0.5, f(1) = macheps^(1/(a+d))
-    leps <- log(.Machine$double.eps)
-    a <- leps / ad[i2]
-    b <- 1 - a / leps * 2
-    ret[i2] <- exp(a + b * log(x[i2]))
+    # f(zero.tol) = macheps^(1/(a+d-0.5)), f(1) = macheps^(1/(a+d))
+    leps  <- log(.Machine$double.eps)
+    b0    <- leps / ad[i2]  # Passes through (x=1, y = eps^(1/ad))
+    lz    <- log(zero.tol[i2])
+    y0log <- leps / (ad[i2] - 0.5)  # log y at x = zero.tol
+    b1    <- (y0log - b0) / lz  # Slope in log-log space
+    ret[i2] <- exp(b0 + b1 * log(x[i2]))
   }
   return(ret)
 }

man/GenD.Rd

@@ -27,7 +27,7 @@ test_that("Jacobian works for scalar functions", {
   }
   p0 <- structure(c(0.1, 0.2, 0.3), names = c("A", "B", "C"))
   h1.jac <- Jacobian(h1, p0)
-  h2.jac <- pnd::Jacobian(h2, p0)
+  h2.jac <- Jacobian(h2, p0)
   expect_true(inherits(h1.jac, "GenD"))
   expect_true(inherits(h2.jac, "GenD"))
   expect_identical(dim(h1.jac), c(1L, 3L))

man/Grad.Rd

@@ -20,9 +20,10 @@ test_that("step size is valid for various values of x", {
   h0 <- stepx(10^(-16:3))
   h <- attr(GenD(sin, x = 10^(-16:3)), "step.size")
   expect_identical(h0, h)
-  expect_true(all(diff(h) >= 0))
-  expect_identical(h[1], h[2])  # Constant step for small x
-  expect_identical(h[length(h)], 1000*.Machine$double.eps^(1/3))
+  expect_true(all(diff(h0) >= 0))
+  expect_identical(h0[1], h0[2])  # Constant step for small x
+  # The following depends on the defaults in stepx!
+  expect_identical(h0[length(h0)], 1000*.Machine$double.eps^(1/3))
 
   expect_true(all(diff(stepx(10^(-16:3), deriv.order = 2, acc.order = 4)) >= 0))
 })