1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
|
# hjn.R -- R implementation of J Nash BASIC HJG.BAS 20160705
hjn <- function(par, fn, lower=-Inf, upper=Inf, bdmsk=NULL, control=list(trace=0), ...){
n <- length(par) # number of parameters
if (! is.null(control$maximize) && control$maximize)
stop("Do NOT try to maximize with hjn()")
if (is.null(control$trace)) control$trace <- 0 # just in case
if (is.null(control$stepsize)) {
stepsize <- 1 # initial step size (could put in control())
} else { stepsize <- control$stepsize }
# Want stepsize positive or bounds get messed up
if (is.null(control$stepredn)) {
stepredn <- .1 # defined step reduction (again in control()??)
} else { stepredn <- control$stepredn }
if (is.null(control$maxfeval)) control$maxfeval<-2000*n
if (is.null(control$eps)) control$epsl<-1e-07
steptol <- control$eps
# Hooke and Jeeves with bounds and masks
if (length(upper) == 1) upper <- rep(upper, n)
if (length(lower) == 1) lower <- rep(lower, n)
if (is.null(bdmsk)) {
bdmsk <- rep(1,n)
idx <- 1:n
} else { idx <- which(bdmsk != 0) } # define masks
if (any(lower >= upper)){
warning("hjn: lower >= upper for some parameters -- set masks")
bdmsk[which(lower >= upper)] <- 0
idx <- which(bdmsk != 0)
}
if (control$trace > 0) {
cat("hjn:bdmsk:")
print(bdmsk)
# cat("idx:")
# print(idx)
}
nac <- length(idx)
offset = 100. # get from control() -- used for equality check
if (any(par < lower) || any(par > upper)) stop("hjn: initial parameters out of bounds")
pbase <- par # base parameter set (fold is function value)
f <- fn(par, ...) # fn at base point
fmin <- fold <- f # "best" function so far
pbest <- par # Not really needed
fcount <- 1 # count function evaluations, compare with maxfeval
# cat(fcount, " f=",fold," at ")
# print(par)
# tmp <- readline("cont.")
keepgoing <- TRUE
ccode <- 1 # start assuming won't get to solution before feval limit
while (keepgoing) {
# exploratory search -- fold stays same for whole set of axes
if (control$trace > 0) cat("Exploratory move - stepsize = ",stepsize,"\n")
if (control$trace > 1) {
cat("p-start:")
print(par)
}
for (jj in 1:nac){ # possibly could do this better in R
# use unmasked parameters
j <- idx[jj]
ptmp <- par[j]
doneg <- TRUE # assume we will do negative step
if (ptmp + offset < upper[j] + offset) { # Not on upper bound so do pos step
par[j] <- min(ptmp+stepsize, upper[j])
if ((par[j] + offset) != (ptmp + offset)) {
fcount <- fcount + 1
f <- fn(par, ...)
# cat(fcount, " f=",f," at ")
# print(par)
if (f < fmin) {
fmin <- f
pbest <- par
# cat("*")
doneg <- FALSE # only case where we don't do neg
resetpar <- FALSE
}
# tmp <- readline("cont>")
}
} # end not on upper bound
if (fcount >= control$maxfeval) break
if (doneg) {
resetpar <- TRUE # going to reset the paramter unless we improve
if ((ptmp + offset) > (lower[j] + offset)) { # can do negative step
par[j] <- max((ptmp - stepsize), lower[j])
if ((par[j] + offset) != (ptmp + offset)) {
fcount <- fcount + 1
f <- fn(par, ...)
# cat(fcount, " f=",f," at ")
# print(par)
if (f < fmin) {
fmin <- f
pbest <- par
# cat("*")
resetpar <- FALSE # don't reset parameter
}
# tmp <- readline("cont<")
}
} # neg step possible
} # end doneg
if (resetpar) { par[j] <- ptmp }
} # end loop on axes
if (fcount >= control$maxfeval) {
ccode <- 1
if (control$trace > 0) cat("Function count limit exceeded\n")
ans <- list(par=pbest, value=fmin, counts=c(fcount, NA), convergence=ccode)
return(ans)
}
if (control$trace > 0) {
cat("axial search with stepsize =",stepsize," fn value = ",fmin," after ",fcount," maxfeval =", control$maxfeval,"\n")
}
if (fmin < fold) { # success -- do pattern move (control$trace > 0) cat("Pattern move \n")
if (control$trace > 1) {
cat("PM from:")
print(par)
cat("pbest:")
print(pbest)
}
for (jj in 1:nac) { # Note par is best set of parameters
j <- idx[jj]
ptmp <- 2.0*par[j] - pbase[j]
if (ptmp > upper[j]) ptmp <- upper[j]
if (ptmp < lower[j]) ptmp <- lower[j]
pbase[j] <- par[j]
par[j] <- ptmp
}
fold <- fmin
if (control$trace > 1) {
cat("PM to:")
print(par)
}
# Addition to HJ -- test new base
# fcount <- fcount + 1
# f <- fn(par, ...)
# cat(fcount, " f=",f," at ")
# print(par)
# tmp <- readline("PM point")
# if (f < fmin) {
# if (control$trace > 0) {cat("Use PM point as new base\n")}
# pbest <- pbase <- par
# }
} else { # no success in Axial Seart, so reduce stepsize
if (fcount >= control$maxfeval) {
ccode <- 1
if (control$trace > 0) cat("Function count limit exceeded\n")
ans <- list(par=pbest, value=fmin, counts=c(fcount, NA), convergence=ccode)
return(ans)
}
# first check if point changed
samepoint <- identical((par + offset),(pbase + offset))
if (samepoint) {
stepsize <- stepsize*stepredn
if (control$trace > 1) cat("Reducing step to ",stepsize,"\n")
if (stepsize <= steptol) keepgoing <- FALSE
ccode <- 0 # successful convergence
} else { # return to old base point
if (control$trace > 1) {
cat("return to base at:")
print(pbase)
}
par <- pbase
}
}
if (fcount >= control$maxfeval) {
ccode <- 1
if (control$trace > 0) cat("Function count limit exceeded\n")
ans <- list(par=pbest, value=fmin, counts=c(fcount, NA), convergence=ccode)
return(ans)
}
} # end keepgoing loop
if ( control$trace > 1 ) {
if (identical(pbest, pbase)) {cat("pbase = pbest\n") }
else { cat("BAD!: pbase != pbest\n") }
}
ans <- list(par=pbest, value=fmin, counts=c(fcount, NA), convergence=ccode)
}
|
