R/hjk.R In Rmpfr: R MPFR - Multiple Precision Floating-Point Reliable

Documented in hjkMpfr

```####
####  h o o k e j e e v e s . R  Hooke-Jeeves Minimization Algorithm
####

## From: John C Nash <[email protected]>
## To: Martin Maechler <[email protected]>, Hans Werner Borchers
## 	<[email protected]>
## Subject: Re: Rmpfr for optimization? Minor success.
## Date: Tue, 5 Jun 2012 12:37:19 -0400

## Changing to hjk routine was a bit easier to deal with. I found main changes
## were to wrap output with as.numeric() to allow cat() to function.

## I'll get no prizes for tidy code, but it is running an n=2 Chebyquad
## minimization, and seems to be working on an n=6. This may be a good way to
## encourage the sale of cpu power.

## Best, JN

hjkMpfr <- function(par, fn, control = list(), ...) {
## Following fails when par is mpfr number JN120605
## if (!is.numeric(par))
##    stop("Argument 'par' must be a numeric vector.", call. = FALSE)
n <- length(par)
if (n == 1)
stop("For univariate functions use some different method.", call. = FALSE)

##-- Control list handling ----------
cntrl <- list(tol      = 1.e-06,
maxfeval = Inf,   # set to Inf if no limit wanted
maximize = FALSE, # set to TRUE  for maximization
target   = Inf,   # set to Inf for no restriction
info     = FALSE) # for printing interim information
nmsCo <- match.arg(names(control), choices = names(cntrl), several.ok = TRUE)
if (!is.null(names(control))) cntrl[nmsCo] <- control

tol      <- cntrl\$tol;
maxfeval <- cntrl\$maxfeval
maximize <- cntrl\$maximize
target   <- cntrl\$target
info <- cntrl\$info

scale <- if (maximize) -1 else 1
fun <- match.fun(fn)
f <- function(x) scale * fun(x, ...)

##-- Setting steps and stepsize -----
nsteps <- floor(log2(1/tol))	# number of steps
steps  <- 2^c(-(0:(nsteps-1)))      # decreasing step size
dir <- diag(1, n, n)                # orthogonal directions

x <- par                            # start point
fx <- f(x)                          # smallest value so far
fcount <- 1                         # counts number of function calls

if (info) cat(sprintf("step  nofc %-12s | %20s\n",
"fmin", "xpar"))

##-- Start the main loop ------------
ns <- 0
while (ns < nsteps && fcount < maxfeval && abs(fx) < target) {
ns <- ns + 1
hjs    <- .hjsearch(x, f, steps[ns], dir, fcount, maxfeval, target)
x      <- hjs\$x
fx     <- hjs\$fx
## found  <- hjs\$found
fcount <- fcount + hjs\$finc

if (info)
cat(sprintf("%4d %5d %-12.7g | %-20.15g %-20.15g%s\n",
ns, fcount, as.numeric(fx/scale),
as.numeric(x[1]), as.numeric(x[2]),
if(n > 2)" ...."))
}

conv <-
if (fcount > maxfeval) {
warning("Function evaluation limit exceeded -- may not converge.")
FALSE
} else if (abs(fx) > target) {
warning("Function exceeds min/max value -- may not converge.")
FALSE
} else
TRUE
fx <- fx / scale                    # undo scaling
list(par = x, value = fx,
convergence = conv, feval = fcount, niter = ns)
}

##  Search with a single scale -----------------------------
.hjsearch <- function(xb, f, h, dir, fcount, maxfeval, target) {
xc <- x <- xb
finc <- 0
hje  <- .hjexplore(xb, xc, f, h, dir)
x    <- hje\$x
fx   <- hje\$fx
found <- hje\$found
finc <- finc + hje\$numf

## Pattern move
while (found) {
d  <- x-xb
xb <- x
xc <- x+d
fb <- fx
hje  <- .hjexplore(xb, xc, f, h, dir, fb)
x    <- hje\$x
fx   <- hje\$fx
found <- hje\$found
finc <- finc + hje\$numf

if (!found) {  # pattern move failed
hje  <- .hjexplore(xb, xb, f, h, dir, fb)
x    <- hje\$x
fx   <- hje\$fx
found <- hje\$found
finc <- finc + hje\$numf
}
if (fcount + finc > maxfeval || abs(fx) > target) break
}

list(x = x, fx = fx, found=found, finc=finc)
}

##  Exploratory move ---------------------------------------
.hjexplore <- function(xb, xc, f, h, dir, fbold) {
n <- length(xb)
x <- xb

if (missing(fbold)) {
fb <- f(x)
numf <- 1
} else {
fb <- fbold
numf <- 0
}

fx <- fb
xt <- xc
found <- FALSE # do we find a better point ?
dirh <- h * dir
fbold <- fx
for (k in sample.int(n, n)) {       # resample orthogonal directions
p <- xt + (d. <- dirh[, k])
fp <- f(p)
numf <- numf + 1

if (fp >= fb) {
p <- xt - d.
fp <- f(p)
numf <- numf + 1
}
if (fp < fb) {
found <- TRUE
xt <- p
fb <- fp
}
}
if(found) {
x <- xt
fx <- fb
}
list(x = x, fx = fx, found=found, numf = numf)
}
```

Try the Rmpfr package in your browser

Any scripts or data that you put into this service are public.

Rmpfr documentation built on July 31, 2018, 1 a.m.