1####
2####  h o o k e j e e v e s . R  Hooke-Jeeves Minimization Algorithm
3####
4
5## From: John C Nash <nashjc@uottawa.ca>
6## To: Martin Maechler <maechler@stat.math.ethz.ch>, Hans Werner Borchers
7## 	<hwborchers@googlemail.com>
8## Subject: Re: Rmpfr for optimization? Minor success.
9## Date: Tue, 5 Jun 2012 12:37:19 -0400
10
11## Changing to hjk routine was a bit easier to deal with. I found main changes
12## were to wrap output with as.numeric() to allow cat() to function.
13
14## I'll get no prizes for tidy code, but it is running an n=2 Chebyquad
15## minimization, and seems to be working on an n=6. This may be a good way to
16## encourage the sale of cpu power.
17
18## Best, JN
19
20hjkMpfr <- function(par, fn, control = list(), ...) {
21    ## Following fails when par is mpfr number JN120605
22    ## if (!is.numeric(par))
23    ##    stop("Argument 'par' must be a numeric vector.", call. = FALSE)
24    n <- length(par)
25    if (n == 1)
26        stop("For univariate functions use some different method.", call. = FALSE)
27
28    ##-- Control list handling ----------
29    cntrl <- list(tol      = 1.e-06,
30                  maxfeval = Inf,   # set to Inf if no limit wanted
31                  maximize = FALSE, # set to TRUE  for maximization
32                  target   = Inf,   # set to Inf for no restriction
33                  info     = FALSE) # for printing interim information
34    nmsCo <- match.arg(names(control), choices = names(cntrl), several.ok = TRUE)
35    if (!is.null(names(control))) cntrl[nmsCo] <- control
36
37    tol      <- cntrl$tol;
38    maxfeval <- cntrl$maxfeval
39    maximize <- cntrl$maximize
40    target   <- cntrl$target
41    info <- cntrl$info
42
43    scale <- if (maximize) -1 else 1
44    fun <- match.fun(fn)
45    f <- function(x) scale * fun(x, ...)
46
47    ##-- Setting steps and stepsize -----
48    nsteps <- floor(log2(1/tol))	# number of steps
49    steps  <- 2^c(-(0:(nsteps-1)))      # decreasing step size
50    dir <- diag(1, n, n)                # orthogonal directions
51
52    x <- par                            # start point
53    fx <- f(x)                          # smallest value so far
54    fcount <- 1                         # counts number of function calls
55
56    if (info) cat(sprintf("step  nofc %-12s | %20s\n",
57                          "fmin", "xpar"))
58
59    ##-- Start the main loop ------------
60    ns <- 0
61    while (ns < nsteps && fcount < maxfeval && abs(fx) < target) {
62        ns <- ns + 1
63        hjs    <- .hjsearch(x, f, steps[ns], dir, fcount, maxfeval, target)
64        x      <- hjs$x
65        fx     <- hjs$fx
66        ## found  <- hjs$found
67        fcount <- fcount + hjs$finc
68
69	if (info)
70	    cat(sprintf("%4d %5d %-12.7g | %-20.15g %-20.15g%s\n",
71			ns, fcount, as.numeric(fx/scale),
72			as.numeric(x[1]), as.numeric(x[2]),
73			if(n > 2)" ...."))
74    }
75
76    conv <-
77	if (fcount > maxfeval) {
78	    warning("Function evaluation limit exceeded -- may not converge.")
79	    FALSE
80	} else if (abs(fx) > target) {
81	    warning("Function exceeds min/max value -- may not converge.")
82	    FALSE
83	} else
84	    TRUE
85    fx <- fx / scale                    # undo scaling
86    list(par = x, value = fx,
87         convergence = conv, feval = fcount, niter = ns)
88}
89
90##  Search with a single scale -----------------------------
91.hjsearch <- function(xb, f, h, dir, fcount, maxfeval, target) {
92    xc <- x <- xb
93    finc <- 0
94    hje  <- .hjexplore(xb, xc, f, h, dir)
95    x    <- hje$x
96    fx   <- hje$fx
97    found <- hje$found
98    finc <- finc + hje$numf
99
100    ## Pattern move
101    while (found) {
102        d  <- x-xb
103        xb <- x
104        xc <- x+d
105        fb <- fx
106        hje  <- .hjexplore(xb, xc, f, h, dir, fb)
107        x    <- hje$x
108        fx   <- hje$fx
109        found <- hje$found
110        finc <- finc + hje$numf
111
112        if (!found) {  # pattern move failed
113           hje  <- .hjexplore(xb, xb, f, h, dir, fb)
114           x    <- hje$x
115           fx   <- hje$fx
116           found <- hje$found
117           finc <- finc + hje$numf
118        }
119        if (fcount + finc > maxfeval || abs(fx) > target) break
120    }
121
122    list(x = x, fx = fx, found=found, finc=finc)
123}
124
125##  Exploratory move ---------------------------------------
126.hjexplore <- function(xb, xc, f, h, dir, fbold) {
127    n <- length(xb)
128    x <- xb
129
130    if (missing(fbold)) {
131        fb <- f(x)
132        numf <- 1
133    } else {
134        fb <- fbold
135        numf <- 0
136    }
137
138    fx <- fb
139    xt <- xc
140    found <- FALSE # do we find a better point ?
141    dirh <- h * dir
142    fbold <- fx
143    for (k in sample.int(n, n)) {       # resample orthogonal directions
144        p <- xt + (d. <- dirh[, k])
145        fp <- f(p)
146        numf <- numf + 1
147
148        if (fp >= fb) {
149            p <- xt - d.
150            fp <- f(p)
151            numf <- numf + 1
152        }
153        if (fp < fb) {
154            found <- TRUE
155            xt <- p
156            fb <- fp
157        }
158    }
159    if(found) {
160        x <- xt
161        fx <- fb
162    }
163    list(x = x, fx = fx, found=found, numf = numf)
164}
165