vuniroot.R
vuniroot <-
function (f, interval, ..., lower=pmin(interval[,1], interval[,2]),
upper=pmax(interval[,1], interval[,2]),
f.lower = f(lower, ...), f.upper = f(upper, ...),
extendInt = c("no", "yes", "downX", "upX"),
check.conv = FALSE, tol = .Machine$double.eps^0.25,
maxiter = 1000, trace = 0, n=NULL)
{
if (!missing(interval) && ncol(interval) != 2L)
stop("'interval' must be a matrix with two columns")
if (all(!is.numeric(lower) | !is.numeric(upper) | lower >= upper))
stop("lower < upper is not fulfilled")
if (is.null(n) && length(lower) == 1 && length(upper) == 1 && length(f.lower) > 1)
n <- length(f.lower)
## warning("Bounds have length 1: possibly replicate for the bounds to have the correct length")
if (!is.null(n) && length(lower) == 1 && length(upper) == 1) {
lower = rep(lower,n)
upper = rep(upper,n)
f.lower = f(lower, ...)
f.upper = f(upper, ...)
}
if (any(is.na(f.lower)))
stop("f.lower = f(lower) is NA at least once")
if (any(is.na(f.upper)))
stop("f.upper = f(upper) is NA at least once")
Sig <- switch(match.arg(extendInt), yes = NULL, downX = -1,
no = 0, upX = 1, stop("invalid 'extendInt'; please report"))
truncate <- function(x) pmax.int(pmin(x, .Machine$double.xmax),
-.Machine$double.xmax)
NAF <- function(x) ifelse(is.na(x),FALSE,x)
f.low. <- truncate(f.lower)
f.upp. <- truncate(f.upper)
fun <- function(x) f(x, ...)
doX <- (is.null(Sig) && any(f.low. * f.upp. > 0, na.rm=TRUE) ||
(is.numeric(Sig) &&
(any(Sig * f.low. > 0, na.rm=TRUE) || any(Sig * f.upp. < 0, na.rm=TRUE))))
if (doX) {
if (trace)
cat(sprintf("search in [%g,%g]%s", lower, upper,
if (trace >= 2)
"\n"
else " ... "))
Delta <- function(l,u) 0.01 * cbind(pmax(1e-04, abs(l)), pmax(1e-04, abs(u)))
it <- 0L
delta <- Delta(lower, upper)
if (is.null(Sig)) {
while (any(i <- NAF(f.lower * f.upper > 0)) &&
(any((iFl <- is.finite(lower)) | (iFu <- is.finite(upper))))) {
if ((it <- it + 1L) > maxiter)
stop(gettextf("no sign change found in %d iterations",
it - 1), domain = NA)
if (any(j <- iFl & i)) {
ol <- lower[j]
of <- f.lower[j]
lower[j] <- lower[j] - delta[j,1]
f.lower[j] <- fun(lower)[j]
if (any(k <- is.na(f.lower[j]))) {
lower[j][k] <- ol[k]
f.lower[j][k] <- of[k]
delta[j,1][k] <- delta[j,1][k]/4
}
}
if (any(j <- iFu & i)) {
ol <- upper[j]
of <- f.upper[j]
upper[j] <- upper[j] + delta[j,2]
f.upper[j] <- fun(upper)[j]
if (any(k <- is.na(f.upper[j]))) {
upper[j][k] <- ol[k]
f.upper[j][k] <- of[k]
delta[j,2][k] <- delta[j,2][k]/4
}
}
if (trace >= 2)
cat(sprintf(" .. modified lower,upper: (%15g,%15g)\n",
lower, upper))
delta[i & (iFl | iFu)] <- 2 * delta[i & (iFl | iFu)]
}
}
else {
while (any(i <- NAF(Sig * f.lower > 0))) {
if ((it <- it + 1L) > maxiter)
stop(gettextf("no sign change found in %d iterations",
it - 1), domain = NA)
lower[i] <- lower[i]-delta[i,1]
f.lower[i] <- fun(lower)[i]
if (trace >= 2)
cat(sprintf(" .. modified lower: %g\n", lower))
delta[i] <- 2 * delta[i]
}
while (any(i <- NAF(Sig * f.upper < 0))) {
if ((it <- it + 1L) > maxiter)
stop(gettextf("no sign change found in %d iterations",
it - 1), domain = NA)
upper[i] <- upper[i] + delta[i,2]
f.upper[i] <- f(upper, ...)[i]
if (trace >= 2)
cat(sprintf(" .. modified upper: %g\n", upper))
delta[i] <- 2 * delta[i]
}
}
if (trace && trace < 2)
cat(sprintf("extended to [%g, %g] in %d steps\n",
lower, upper, it))
}
if (any(!NAF(as.vector(f.lower * f.upper <= 0))))
stop(if (doX)
"did not succeed extending the interval endpoints for f(lower) * f(upper) <= 0"
else "f() values at end points not of opposite sign")
if (check.conv) {
val <- tryCatch(.Call("vunirootRcpp", fun, lower, upper, f.lower, f.upper, maxiter, tol, PACKAGE="rstpm2"),
warning = function(w) w)
if (inherits(val, "warning"))
stop("convergence problem in zero finding: ", conditionMessage(val))
}
else {
val <- .Call("vunirootRcpp", fun, lower, upper, f.lower, f.upper, as.integer(maxiter), tol, PACKAGE="rstpm2")
}
iter <- as.integer(val[[2L]])
if (any(iter < 0)) {
(if (check.conv)
stop
else warning)(sprintf(ngettext(maxiter, "_NOT_ converged in %d iteration",
"_NOT_ converged in %d iterations"), maxiter), domain = NA)
iter[iter<0] <- maxiter
}
if (doX)
iter <- iter + it
else it <- NA_integer_
list(root = val[[1L]], f.root = f(val[[1L]], ...), iter = iter,
init.it = it, estim.prec = val[[3L]])
}
voptimize <- function (f, interval, ...,
lower=pmin(interval[,1], interval[,2]),
upper=pmax(interval[,1], interval[,2]),
maximum = FALSE, tol = .Machine$double.eps^0.25)
{
if (maximum) {
val <- .Call("voptimizeRcpp", function(arg) -f(arg, ...),
lower, upper, tol, PACKAGE="rstpm2")
list(maximum = val, objective = f(val, ...))
}
else {
val <- .Call("voptimizeRcpp", function(arg) f(arg, ...),
lower, upper, tol, PACKAGE="rstpm2")
list(minimum = val, objective = f(val, ...))
}
}
voptimise <- voptimize
vintegrate = function(f, lower, upper, ..., subdivisions = 100L, rel.tol=.Machine$double.eps^0.25,
abs.tol=rel.tol) {
f <- match.fun(f)
ff <- function(x) f(x, ...)
ny = max(length(lower), length(upper))
if (ny>1 && length(lower)==1) lower = rep(lower,ny)
if (ny>1 && length(upper)==1) upper = rep(upper,ny)
limit <- as.integer(subdivisions)
if (limit < 1L || (abs.tol <= 0 && rel.tol < max(50 * .Machine$double.eps,
5e-29)))
stop("invalid parameter values")
if (any(is.na(lower) | is.na(upper)))
stop("a limit is NA or NaN")
if (!(all(is.finite(lower)) || all(!is.finite(lower))))
stop("lower bounds are a mix of finite and infinite values")
if (!(all(is.finite(upper)) || all(!is.finite(upper))))
stop("upper bounds are a mix of finite and infinite values")
if (all(is.finite(lower)) && all(is.finite(upper))) {
res = vdqagsRcpp(ff, lower, upper, rel.tol, abs.tol, limit, ny)
res$value = drop(res$value)
res$abs.err = drop(res$abs.err)
res$abserr = NULL
} else {
if (all(is.finite(lower))) {
inf <- 1L
bound <- lower
}
else if (all(is.finite(upper))) {
inf <- -1L
bound <- upper
}
else {
inf <- 2L
bound <- rep(0,ny) # ignored:)
}
res = vdqagiRcpp(ff, bound, inf, rel.tol, abs.tol, limit, ny)
res$value = drop(res$value)
res$abserr = drop(res$abserr)
}
res
}
