Setup a root_control object
setup_root_control(FUN, roots_name, ...)a function
a character string identifying the object containing the
arguments passed to FUN
a root_control object
# Setup the default
setup_root_control(start = c(3, 5, 6))
#> An object of class "root_control"
#> Slot ".object_name":
#> [1] "root"
#>
#> Slot ".FUN":
#> function (f, start, maxiter = 100, rtol = 1e-06, atol = 1e-08,
#> ctol = 1e-08, useFortran = TRUE, positive = FALSE, jacfunc = NULL,
#> jactype = "fullint", verbose = FALSE, bandup = 1, banddown = 1,
#> parms = NULL, ...)
#> {
#> initfunc <- NULL
#> if (is.list(f)) {
#> if (!is.null(jacfunc) & "jacfunc" %in% names(f))
#> stop("If 'f' is a list that contains jacfunc, argument 'jacfunc' should be NULL")
#> jacfunc <- f$jacfunc
#> initfunc <- f$initfunc
#> f <- f$func
#> }
#> N <- length(start)
#> if (!is.numeric(start))
#> stop("start conditions should be numeric")
#> if (!is.numeric(maxiter))
#> stop("`maxiter' must be numeric")
#> if (as.integer(maxiter) < 1)
#> stop("maxiter must be >=1")
#> if (!is.numeric(rtol))
#> stop("`rtol' must be numeric")
#> if (!is.numeric(atol))
#> stop("`atol' must be numeric")
#> if (!is.numeric(ctol))
#> stop("`ctol' must be numeric")
#> if (length(atol) > 1 && length(atol) != N)
#> stop("`atol' must either be a scalar, or as long as `start'")
#> if (length(rtol) > 1 && length(rtol) != N)
#> stop("`rtol' must either be a scalar, or as long as `y'")
#> if (length(ctol) > 1)
#> stop("`ctol' must be a scalar")
#> if (useFortran) {
#> if (!is.compiled(f) & is.null(parms)) {
#> Fun1 <- function(time = 0, x, parms = NULL) list(f(x,
#> ...))
#> Fun <- Fun1
#> }
#> else if (!is.compiled(f)) {
#> Fun2 <- function(time = 0, x, parms) list(f(x, parms,
#> ...))
#> Fun <- Fun2
#> }
#> else {
#> Fun <- f
#> f <- function(x, ...) Fun(n = length(start), t = 0,
#> x, f = rep(0, length(start)), 1, 1)$f
#> }
#> JacFunc <- jacfunc
#> if (!is.null(jacfunc))
#> if (!is.compiled(JacFunc) & is.null(parms))
#> JacFunc <- function(time = 0, x, parms = parms) jacfunc(x,
#> ...)
#> else if (!is.compiled(JacFunc))
#> JacFunc <- function(time = 0, x, parms = parms) jacfunc(x,
#> parms, ...)
#> else JacFunc <- jacfunc
#> method <- "stode"
#> if (jactype == "sparse") {
#> method <- "stodes"
#> if (!is.null(jacfunc))
#> stop("jacfunc can not be used when jactype='sparse'")
#> x <- stodes(y = start, time = 0, func = Fun, atol = atol,
#> positive = positive, rtol = rtol, ctol = ctol,
#> maxiter = maxiter, verbose = verbose, parms = parms,
#> initfunc = initfunc)
#> }
#> else x <- steady(y = start, time = 0, func = Fun, atol = atol,
#> positive = positive, rtol = rtol, ctol = ctol, maxiter = maxiter,
#> method = method, jacfunc = JacFunc, jactype = jactype,
#> verbose = verbose, parms = parms, initfunc = initfunc,
#> bandup = bandup, banddown = banddown)
#> precis <- attr(x, "precis")
#> attributes(x) <- NULL
#> x <- unlist(x)
#> if (is.null(parms))
#> reffx <- f(x, ...)
#> else reffx <- f(x, parms, ...)
#> i <- length(precis)
#> }
#> else {
#> if (is.compiled(f))
#> stop("cannot combine compiled code with R-implemented solver")
#> precis <- NULL
#> x <- start
#> jacob <- matrix(nrow = N, ncol = N, data = 0)
#> if (is.null(parms))
#> reffx <- f(x, ...)
#> else reffx <- f(x, parms, ...)
#> if (length(reffx) != N)
#> stop("'f', function must return as many function values as elements in start")
#> for (i in 1:maxiter) {
#> refx <- x
#> pp <- mean(abs(reffx))
#> precis <- c(precis, pp)
#> ewt <- rtol * abs(x) + atol
#> if (max(abs(reffx/ewt)) < 1)
#> break
#> delt <- perturb(x)
#> for (j in 1:N) {
#> x[j] <- x[j] + delt[j]
#> if (is.null(parms))
#> fx <- f(x, ...)
#> else fx <- f(x, parms, ...)
#> jacob[, j] <- (fx - reffx)/delt[j]
#> x[j] <- refx[j]
#> }
#> relchange <- as.numeric(solve(jacob, -1 * reffx))
#> if (max(abs(relchange)) < ctol)
#> break
#> x <- x + relchange
#> if (is.null(parms))
#> reffx <- f(x, ...)
#> else reffx <- f(x, parms, ...)
#> }
#> }
#> names(x) <- names(start)
#> return(list(root = x, f.root = reffx, iter = i, estim.precis = precis[length(precis)]))
#> }
#> <bytecode: 0x7fea475c8ef8>
#> <environment: namespace:rootSolve>
#>
#> Slot ".options":
#> $start
#> [1] 3 5 6
#>
#>
# Also setup the default
setup_root_control(FUN = rootSolve::multiroot,
start = c(3, 5, 6))
#> An object of class "root_control"
#> Slot ".object_name":
#> [1] "root"
#>
#> Slot ".FUN":
#> function (f, start, maxiter = 100, rtol = 1e-06, atol = 1e-08,
#> ctol = 1e-08, useFortran = TRUE, positive = FALSE, jacfunc = NULL,
#> jactype = "fullint", verbose = FALSE, bandup = 1, banddown = 1,
#> parms = NULL, ...)
#> {
#> initfunc <- NULL
#> if (is.list(f)) {
#> if (!is.null(jacfunc) & "jacfunc" %in% names(f))
#> stop("If 'f' is a list that contains jacfunc, argument 'jacfunc' should be NULL")
#> jacfunc <- f$jacfunc
#> initfunc <- f$initfunc
#> f <- f$func
#> }
#> N <- length(start)
#> if (!is.numeric(start))
#> stop("start conditions should be numeric")
#> if (!is.numeric(maxiter))
#> stop("`maxiter' must be numeric")
#> if (as.integer(maxiter) < 1)
#> stop("maxiter must be >=1")
#> if (!is.numeric(rtol))
#> stop("`rtol' must be numeric")
#> if (!is.numeric(atol))
#> stop("`atol' must be numeric")
#> if (!is.numeric(ctol))
#> stop("`ctol' must be numeric")
#> if (length(atol) > 1 && length(atol) != N)
#> stop("`atol' must either be a scalar, or as long as `start'")
#> if (length(rtol) > 1 && length(rtol) != N)
#> stop("`rtol' must either be a scalar, or as long as `y'")
#> if (length(ctol) > 1)
#> stop("`ctol' must be a scalar")
#> if (useFortran) {
#> if (!is.compiled(f) & is.null(parms)) {
#> Fun1 <- function(time = 0, x, parms = NULL) list(f(x,
#> ...))
#> Fun <- Fun1
#> }
#> else if (!is.compiled(f)) {
#> Fun2 <- function(time = 0, x, parms) list(f(x, parms,
#> ...))
#> Fun <- Fun2
#> }
#> else {
#> Fun <- f
#> f <- function(x, ...) Fun(n = length(start), t = 0,
#> x, f = rep(0, length(start)), 1, 1)$f
#> }
#> JacFunc <- jacfunc
#> if (!is.null(jacfunc))
#> if (!is.compiled(JacFunc) & is.null(parms))
#> JacFunc <- function(time = 0, x, parms = parms) jacfunc(x,
#> ...)
#> else if (!is.compiled(JacFunc))
#> JacFunc <- function(time = 0, x, parms = parms) jacfunc(x,
#> parms, ...)
#> else JacFunc <- jacfunc
#> method <- "stode"
#> if (jactype == "sparse") {
#> method <- "stodes"
#> if (!is.null(jacfunc))
#> stop("jacfunc can not be used when jactype='sparse'")
#> x <- stodes(y = start, time = 0, func = Fun, atol = atol,
#> positive = positive, rtol = rtol, ctol = ctol,
#> maxiter = maxiter, verbose = verbose, parms = parms,
#> initfunc = initfunc)
#> }
#> else x <- steady(y = start, time = 0, func = Fun, atol = atol,
#> positive = positive, rtol = rtol, ctol = ctol, maxiter = maxiter,
#> method = method, jacfunc = JacFunc, jactype = jactype,
#> verbose = verbose, parms = parms, initfunc = initfunc,
#> bandup = bandup, banddown = banddown)
#> precis <- attr(x, "precis")
#> attributes(x) <- NULL
#> x <- unlist(x)
#> if (is.null(parms))
#> reffx <- f(x, ...)
#> else reffx <- f(x, parms, ...)
#> i <- length(precis)
#> }
#> else {
#> if (is.compiled(f))
#> stop("cannot combine compiled code with R-implemented solver")
#> precis <- NULL
#> x <- start
#> jacob <- matrix(nrow = N, ncol = N, data = 0)
#> if (is.null(parms))
#> reffx <- f(x, ...)
#> else reffx <- f(x, parms, ...)
#> if (length(reffx) != N)
#> stop("'f', function must return as many function values as elements in start")
#> for (i in 1:maxiter) {
#> refx <- x
#> pp <- mean(abs(reffx))
#> precis <- c(precis, pp)
#> ewt <- rtol * abs(x) + atol
#> if (max(abs(reffx/ewt)) < 1)
#> break
#> delt <- perturb(x)
#> for (j in 1:N) {
#> x[j] <- x[j] + delt[j]
#> if (is.null(parms))
#> fx <- f(x, ...)
#> else fx <- f(x, parms, ...)
#> jacob[, j] <- (fx - reffx)/delt[j]
#> x[j] <- refx[j]
#> }
#> relchange <- as.numeric(solve(jacob, -1 * reffx))
#> if (max(abs(relchange)) < ctol)
#> break
#> x <- x + relchange
#> if (is.null(parms))
#> reffx <- f(x, ...)
#> else reffx <- f(x, parms, ...)
#> }
#> }
#> names(x) <- names(start)
#> return(list(root = x, f.root = reffx, iter = i, estim.precis = precis[length(precis)]))
#> }
#> <bytecode: 0x7fea12964ae0>
#> <environment: namespace:rootSolve>
#>
#> Slot ".options":
#> $start
#> [1] 3 5 6
#>
#>
# Or use uniroot()
setup_root_control(FUN = stats::uniroot,
interval = c(0, 1))
#> An object of class "root_control"
#> Slot ".object_name":
#> [1] "root"
#>
#> Slot ".FUN":
#> function (f, interval, ..., lower = min(interval), upper = max(interval),
#> 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)
#> {
#> if (!missing(interval) && length(interval) != 2L)
#> stop("'interval' must be a vector of length 2")
#> if (!is.numeric(lower) || !is.numeric(upper) || lower >=
#> upper)
#> stop("lower < upper is not fulfilled")
#> if (is.na(f.lower))
#> stop("f.lower = f(lower) is NA")
#> if (is.na(f.upper))
#> stop("f.upper = f(upper) is NA")
#> 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)
#> f.low. <- truncate(f.lower)
#> f.upp. <- truncate(f.upper)
#> doX <- (is.null(Sig) && f.low. * f.upp. > 0 || is.numeric(Sig) &&
#> (Sig * f.low. > 0 || Sig * f.upp. < 0))
#> if (doX) {
#> if (trace)
#> cat(sprintf("search in [%g,%g]%s", lower, upper,
#> if (trace >= 2)
#> "\n"
#> else " ... "))
#> Delta <- function(u) 0.01 * pmax(1e-04, abs(u))
#> it <- 0L
#> if (is.null(Sig)) {
#> delta <- Delta(c(lower, upper))
#> while (isTRUE(f.lower * f.upper > 0) && any(iF <- is.finite(c(lower,
#> upper)))) {
#> if ((it <- it + 1L) > maxiter)
#> stop(gettextf("no sign change found in %d iterations",
#> it - 1), domain = NA)
#> if (iF[1]) {
#> ol <- lower
#> of <- f.lower
#> if (is.na(f.lower <- f(lower <- lower - delta[1],
#> ...))) {
#> lower <- ol
#> f.lower <- of
#> delta[1] <- delta[1]/4
#> }
#> }
#> if (iF[2]) {
#> ol <- upper
#> of <- f.upper
#> if (is.na(f.upper <- f(upper <- upper + delta[2],
#> ...))) {
#> upper <- ol
#> f.upper <- of
#> delta[2] <- delta[2]/4
#> }
#> }
#> if (trace >= 2)
#> cat(sprintf(" .. modified lower,upper: (%15g,%15g)\n",
#> lower, upper))
#> delta <- 2 * delta
#> }
#> }
#> else {
#> delta <- Delta(lower)
#> while (isTRUE(Sig * f.lower > 0)) {
#> if ((it <- it + 1L) > maxiter)
#> stop(gettextf("no sign change found in %d iterations",
#> it - 1), domain = NA)
#> f.lower <- f(lower <- lower - delta, ...)
#> if (trace >= 2)
#> cat(sprintf(" .. modified lower: %g\n", lower))
#> delta <- 2 * delta
#> }
#> delta <- Delta(upper)
#> while (isTRUE(Sig * f.upper < 0)) {
#> if ((it <- it + 1L) > maxiter)
#> stop(gettextf("no sign change found in %d iterations",
#> it - 1), domain = NA)
#> f.upper <- f(upper <- upper + delta, ...)
#> if (trace >= 2)
#> cat(sprintf(" .. modified upper: %g\n", upper))
#> delta <- 2 * delta
#> }
#> }
#> if (trace && trace < 2)
#> cat(sprintf("extended to [%g, %g] in %d steps\n",
#> lower, upper, it))
#> }
#> if (!isTRUE(as.vector(sign(f.lower) * sign(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(.External2(C_zeroin2, function(arg) f(arg,
#> ...), lower, upper, f.lower, f.upper, tol, as.integer(maxiter)),
#> warning = function(w) w)
#> if (inherits(val, "warning"))
#> stop("convergence problem in zero finding: ", conditionMessage(val))
#> }
#> else {
#> val <- .External2(C_zeroin2, function(arg) f(arg, ...),
#> lower, upper, f.lower, f.upper, tol, as.integer(maxiter))
#> }
#> iter <- as.integer(val[2L])
#> if (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 <- 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])
#> }
#> <bytecode: 0x7fea129efb48>
#> <environment: namespace:stats>
#>
#> Slot ".options":
#> $interval
#> [1] 0 1
#>
#>