R/rootwrap.R

In rootoned: Roots of one-dimensional functions in R-only code

rootwrap <- function(fn=fn, gr=NULL, ri=NULL, method="uniroot", ftrace=TRUE, ...){
  ## A wrapper for a variety of rootfinders
  ## We will try to document here
  # fn = the function for which root is to be found
  # ri = interval (2 element vector) of lower and upper bound to root
  #    OR  if second element is NA, then a guess to start the rootfinding
  # method = character name of the rootfinder (list below)
  # ftrace = logical TRUE if method is to be traced
  # ... = additional arguments to fn or FnTrace (e.g., ifn)
  ## ?? do we want to do a preliminary search? Esp. for GUESS methods
  
  ## OUTPUT
  #  root
  #  froot
  #  rtol
  #  niter (may be same as fncount)
  #  fncount
  #  method
  # ?? NOT INCLUDED grcount -- only newt1d uses gr now

  dots <- list(...)
  groot<-list(ifn=0, igr=0, ftrace=ftrace, fn=fn, gr=gr, label=method, dots=dots)
  # groot<-list(ifn=0, igr=0, ftrace=ftrace, fn=fn, gr=gr, label=method, ...)
  envroot <- list2env(groot)
  #    print(envroot)
  #    envroot
  #    print(envroot$dots)
  
  
FnTrace <- function(x) { 
  # Substitute function to call when rootfinding
  # Evaluate fn(x, ...)
#    cat("FnTrace, envroot$dots =")
#    print(envroot$dots)
    val <- do.call(envroot$fn, c(x=x, envroot$dots))
    envroot$ifn <- envroot$ifn + 1 # probably more efficient ways
    if (envroot$ftrace) {
       cat("f(",x,")=",val," after ",envroot$ifn," ",envroot$label,"\n")
    }
    val
}

grTrace <- function(x) { 
  # Substitute function to call when rootfinding
  # Evaluate fn(x, ...)
  val <- do.call(envroot$gr, c(x=x, envroot$dots))
  envroot$igr <- envroot$igr + 1 # probably more efficient ways
  if (envroot$ftrace) {
    cat("gr(",x,")=",val," after ",envroot$igr," ",envroot$label,"\n")
  }
  val
}

#    cat("test function value at ri[1]=",ri[1],"\n")
#    cat(FnTrace(ri[1]))
  ## Preliminary checks
  # print(ri)
  if (length(ri) != 2) {stop("ri - the interval in which root to be found, must be 2 elements")}
  if (is.na(ri[2])) {
      fguess = ri[1]
  } else if (ri[2] <= ri[1]){stop("Lower bound must be > upper bound")}

  
## Beginning of methods
  
  if (method == "uniroot"){
    tst <- uniroot(FnTrace, interval=ri) # may still need dotargs for tolerances etc.
    tst$froot <- tst$f.root
    tst$f.root <- NULL
    tst$rtol <- tst$estim.prec
    tst$estim.prec <- NULL
    tst$init.it <- NULL
    #   return(tst)
  }
  
  if (method == "root1d"){
    tst <- root1d(f=FnTrace, interval=ri, trace=ftrace)
    tst$iter <- tst$fcount
    tst$fcount <- NULL
  }
  
  if (method == "zeroin"){
    tst <- zeroin(f=FnTrace, ri, trace=ftrace)
  }
  
  if (method == "newt1d"){
    if (is.null(gr)) stop("newt1d in rootwrap MUST have gradient")
    fguess<-ri[1]
    cat("fguess = ", fguess,"\n")
    tst <- newt1d(fn=FnTrace, gr=grTrace, x0=fguess, trace=ftrace)
    tst$iter <- tst$itn
    tst$itn <- NULL
    tst$rtol <- NA
    #   return(tst)
  }
  
  if (method == "newton"){
    fguess<-ri[1]
    tst <- newton(FnTrace, fguess) # newton has dot args in pracma
    tst$froot <- tst$f.root
    tst$f.root <- NULL
    tst$rtol <- tst$estim.prec
    tst$estim.prec <- NULL
    tst$iter <- tst$niter
    tst$niter <- NULL
  }
  
  if (method == "bisect"){
    tst <- bisect(f=FnTrace, ri[1], ri[2]) # No dotargs in pracma
    tst$froot <- tst$f.root
    tst$f.root <- NULL
    tst$rtol <- tst$estim.prec
    tst$estim.prec <- NULL
    #   return(tst)
  }
  
  if (method == "secant"){
    fguess <- ri[1]
    tst <- secant(fun=FnTrace, fguess, (fguess+0.01*(abs(fguess)+1))) # has dotargs in pracma
    tst$froot <- tst$f.root
    tst$f.root <- NULL
    tst$rtol <- tst$estim.prec
    tst$estim.prec <- NULL
    #   return(tst)
  }
  
  if (method == "regulaFalsi"){
    tst <- regulaFalsi(f=FnTrace, ri[1], ri[2]) # no dotargs in pracma
    tst$froot <- tst$f.root
    tst$f.root <- NULL
    tst$rtol <- tst$estim.prec
    tst$estim.prec <- NULL
    tst$iter <- tst$niter
    tst$niter <- NULL
    #   return(tst)
  } 

   if (method == "muller"){
    tst <- muller(f=FnTrace, ri[1], ri[2], 0.5*(ri[1]+ri[2])) # no dotargs in pracma
    tst$iter <- tst$niter
    tst$niter <- NULL
    tst$froot <- tst$fval
    tst$fval <- NULL
    tst$rtol <- tst$reltol
    tst$reltol <- NULL
  #   return(tst)
  } 
  
  
  if (method == "brent"){
    tst <- brent(f=FnTrace, ri[1], ri[2]) # no dotargs in pracma
    tst$froot <- tst$f.root
    tst$f.root <- NULL
    tst$iter <- tst$f.calls
    tst$f.calls <- NULL
    tst$rtol <- tst$estim.prec
    tst$estim.prec <- NULL
  } 
  
  
## End methods
  tst$fncount <- envroot$ifn # save value
  tst$method <- method # save method
#  rm(envroot) # clear the local environment
  tst

}