6 (in-package 'lisp-stat
)
8 (export '(newtonmax nelmeadmax
))
10 (import '(ls-basics::new-minfo-internals ls-basics
::minfo-maximize
))
13 ;;;; Mode Info Prototype
16 (defproto minfo-proto
'(internals))
18 #+xlisp
(send minfo-proto
:add-method
:isnew
#'|minfo-isnew|
)
19 #+xlisp
(send minfo-proto
:add-method
:maximize
#'|minfo-maximize|
)
20 #+xlisp
(send minfo-proto
:add-method
:loglaplace
#'|minfo-loglap|
)
22 (defmeth minfo-proto
:isnew
(&rest args
)
23 (setf (slot-value 'internals
) (apply #'new-minfo-internals args
)))
25 (defmeth minfo-proto
:maximize
(&rest args
)
26 (apply #'minfo-maximize
(slot-value 'internals
) args
))
28 (defmeth minfo-proto
:x
() (aref (slot-value 'internals
) 3))
29 (defmeth minfo-proto
:scale
() (aref (slot-value 'internals
) 4))
30 (defmeth minfo-proto
:derivstep
() (aref (aref (slot-value 'internals
) 9) 1))
31 (defmeth minfo-proto
:tilt
() (aref (aref (slot-value 'internals
) 9) 6))
33 (defmeth minfo-proto
:f
(&optional
(val nil set
))
35 (send self
:set-no-vals-supplied
)
36 (setf (aref (slot-value 'internals
) 0) val
))
37 (aref (slot-value 'internals
) 0))
39 (defmeth minfo-proto
:set-no-vals-supplied
()
40 (setf (aref (aref (slot-value 'internals
) 8) 6) 0))
42 (defmeth minfo-proto
:exptilt
(&optional
(val nil set
))
44 (let ((old (send self
:exptilt
)))
45 (setf (aref (aref (slot-value 'internals
) 8) 7) (if val
1 0))
46 (if (and (not (or (and old val
) (and (not old
) (not val
))))
47 (/= (send self
:tilt
) 0.0))
48 (send self
:set-no-vals-supplied
))))
49 (= 1 (aref (aref (slot-value 'internals
) 8) 7)))
51 (defmeth minfo-proto
:newtilt
(&optional
(val nil set
))
53 (setf (aref (aref (slot-value 'internals
) 9) 7) (float val
))
54 (if (/= (send self
:tilt
) 0.0) (send self
:set-no-vals-supplied
)))
55 (aref (aref (slot-value 'internals
) 9) 7))
57 (defmeth minfo-proto
:gfuns
(&optional
(val nil set
))
59 (if (or (not (consp val
))
60 (not (every #'functionp val
)))
61 (error "not all functions"))
62 (setf (aref (slot-value 'internals
) 1) val
)
63 (setf (aref (aref (slot-value 'internals
) 8) 1) (length val
))
64 (setf (aref (slot-value 'internals
) 10) (repeat 1.0 (length val
)))
65 (if (/= (send self
:tilt
) 0.0) (send self
:set-no-vals-supplied
)))
66 (aref (slot-value 'internals
) 1))
68 (defmeth minfo-proto
:cfuns
(&optional
(val nil set
))
70 (if (or (not (consp val
))
71 (not (every #'functionp val
)))
72 (error "not all functions"))
73 (setf (aref (slot-value 'internals
) 2) val
)
74 (setf (aref (aref (slot-value 'internals
) 8) 2) (length val
))
75 (setf (aref (slot-value 'internals
) 7) (repeat 0.0 (length val
)))
76 (setf (aref (slot-value 'internals
) 11) (repeat 0.0 (length val
)))
77 (send self
:set-no-vals-supplied
))
78 (aref (slot-value 'internals
) 2))
80 (defmeth minfo-proto
:ctarget
(&optional
(val nil set
))
82 (if (/= (length val
) (length (send self
:ctarget
)))
83 (error "bad target length"))
84 (setf (aref (slot-value 'internals
) 7) val
))
85 (aref (slot-value 'internals
) 7))
87 (defmeth minfo-proto
:fvals
()
88 (let* ((fv (aref (slot-value 'internals
) 5))
89 (n (length (send self
:x
)))
91 (grad (select fv
(iseq 1 n
)))
92 (hess (matrix (list n n
) (select fv
(iseq (+ 1 n
) (+ n
(* n n
)))))))
93 (list val grad hess
)))
95 (defmeth minfo-proto
:copy
()
96 (let ((obj (make-object minfo-proto
))
97 (internals (copy-seq (slot-value 'internals
))))
98 (dotimes (i (length internals
))
99 (let ((x (aref internals i
)))
101 (setf (aref internals i
) (copy-seq x
)))))
102 (send obj
:add-slot
'internals internals
)
105 (defmeth minfo-proto
:derivscale
()
106 (let* ((step (^ machine-epsilon
(/ 1 6)))
107 (hess (numhess (send self
:f
) (send self
:x
) (send self
:scale
) step
))
108 (scale (pmax (abs (send self
:x
)) (sqrt (abs (/ (diagonal hess
)))))))
109 (setf hess
(numhess (send self
:f
) (send self
:x
) scale step
))
110 (setf scale
(pmax (abs (send self
:x
)) (sqrt (abs (/ (diagonal hess
))))))
111 (setf (aref (slot-value 'internals
) 4) scale
)
112 (setf (aref (aref (slot-value 'internals
) 9) 1) step
)))
114 (defmeth minfo-proto
:verbose
(&optional
(val nil set
))
116 (setf (aref (aref (slot-value 'internals
) 8) 5)
117 (cond ((integerp val
) val
)
120 (aref (aref (slot-value 'internals
) 8) 5))
122 (defmeth minfo-proto
:backtrack
(&optional
(val nil set
))
123 (if set
(setf (aref (aref (slot-value 'internals
) 8) 4) (if val
1 0)))
124 (aref (aref (slot-value 'internals
) 8) 4))
126 (defmeth minfo-proto
:maxiter
(&optional
(val nil set
))
127 (if set
(setf (aref (aref (slot-value 'internals
) 8) 3)
128 (if (integerp val
) val -
1)))
129 (aref (aref (slot-value 'internals
) 8) 3))
131 (defmeth minfo-proto
:tiltscale
(&optional
(val nil set
))
133 (if (/= (length val
) (length (send self
:gfuns
)))
134 (error "wrong size tilt scale sequence"))
135 (setf (aref (slot-value 'internals
) 10) val
))
136 (aref (slot-value 'internals
) 10))
140 ;;;; Newton's Method with Backtracking
144 (defun newtonmax (f start
&key
150 "Args:(f start &key scale derivstep (verbose 1) return-derivs)
151 Maximizes F starting from START using Newton's method with backtracking.
152 If RETURN-DERIVS is NIL returns location of maximum; otherwise returns
153 list of location, unction value, gradient and hessian at maximum.
154 SCALE should be a list of the typical magnitudes of the parameters.
155 DERIVSTEP is used in numerical derivatives and VERBOSE controls printing
156 of iteration information. COUNT-LIMIT limits the number of iterations"
157 (let ((verbose (if verbose
(if (integerp verbose
) verbose
1) 0))
158 (minfo (send minfo-proto
:new f start
159 :scale scale
:derivstep derivstep
)))
160 (send minfo
:maxiter count-limit
)
161 (send minfo
:derivscale
)
162 (send minfo
:maximize verbose
)
164 (cons (send minfo
:x
) (- (send minfo
:fvals
)))
169 ;;;; Nelder-Mead Simplex Method
173 (defun nelmeadmax (f start
&key
175 (epsilon (sqrt machine-epsilon
))
182 "Args: (f start &key (size 1) (epsilon (sqrt machine-epsilon))
183 (count-limit 500) (verbose t) alpha beta gamma delta)
184 Maximizes F using the Nelder-Mead simplex method. START can be a
185 starting simplex - a list of N+1 points, with N=dimension of problem,
186 or a single point. If start is a single point you should give the
187 size of the initial simplex as SIZE, a sequence of length N. Default is
188 all 1's. EPSILON is the convergence tolerance. ALPHA-DELTA can be used to
189 control the behavior of simplex algorithm."
190 (let ((s (send simplex-proto
:new f start size
)))
191 (do ((best (send s
:best-point
) (send s
:best-point
))
192 (count 0 (+ count
1))
194 ((or (< (send s
:relative-range
) epsilon
) (>= count count-limit
))
195 (if (and verbose
(>= count count-limit
))
196 (format t
"Iteration limit exceeded.~%"))
197 (send s
:point-location
(send s
:best-point
)))
198 (setf next
(send s
:extrapolate-from-worst
(- alpha
)))
199 (if (send s
:is-worse best next
)
200 (setf next
(send s
:extrapolate-from-worst gamma
))
201 (when (send s
:is-worse next
(send s
:second-worst-point
))
202 (setf next
(send s
:extrapolate-from-worst beta
))
203 (if (send s
:is-worse next
(send s
:worst-point
))
204 (send s
:shrink-to-best delta
))))
206 (format t
"Value = ~10g~%"
207 (send s
:point-value
(send s
:best-point
)))))))
211 ;;; Simplex Prototype
214 (defproto simplex-proto
'(f simplex
))
220 (defmeth simplex-proto
:make-point
(x)
221 (let ((f (send self
:f
)))
223 (let ((val (funcall f x
)))
224 (cons (if (consp val
) (car val
) val
) x
))
227 (defmeth simplex-proto
:point-value
(x) (car x
))
229 (defmeth simplex-proto
:point-location
(x) (cdr x
))
231 (defmeth simplex-proto
:is-worse
(x y
)
232 (< (send self
:point-value x
) (send self
:point-value y
)))
235 ;;; Making New Simplices
238 (defmeth simplex-proto
:isnew
(f start
&optional size
)
239 (send self
:simplex start size
)
243 ;;; Slot Accessors and Mutators
246 (defmeth simplex-proto
:simplex
(&optional new size
)
249 (if (and (consp new
) (sequencep (car new
)))
250 (if (/= (length new
) (+ 1 (length (car new
))))
251 (error "bad simplex data")
253 (let* ((n (length new
))
254 (size (if size size
(repeat 1 n
)))
255 ; (pts (- (* 2 (uniform-rand (repeat n (+ n 1)))) 1)))
256 (diag (* 2 size
(- (random (repeat 2 n
)) .5)))
257 (pts (cons (repeat 0 n
)
258 (mapcar #'(lambda (x) (coerce x
'list
))
259 (column-list (diagonal diag
))))))
260 (mapcar #'(lambda (x) (+ (* size x
) new
)) pts
)))))
261 (setf (slot-value 'simplex
)
262 (mapcar #'(lambda (x) (send self
:make-point x
)) simplex
))
263 (send self
:sort-simplex
)))
264 (slot-value 'simplex
))
266 (defmeth simplex-proto
:f
(&optional f
)
268 (setf (slot-value 'f
) f
)
270 (mapcar #'(lambda (x) (send self
:point-location x
))
271 (send self
:simplex
))))
272 (send self
:simplex simplex
)))
275 (defmeth simplex-proto
:sort-simplex
()
277 (setf (slot-value 'simplex
)
278 (sort (slot-value 'simplex
)
279 #'(lambda (x y
) (send self
:is-worse x y
))))))
282 ;;; Other Methods Using List Representation of SImplex
285 (defmeth simplex-proto
:best-point
() (car (last (send self
:simplex
))))
286 (defmeth simplex-proto
:worst-point
() (first (send self
:simplex
)))
287 (defmeth simplex-proto
:second-worst-point
() (second (send self
:simplex
)))
288 (defmeth simplex-proto
:replace-point
(new old
)
289 (let* ((simplex (send self
:simplex
))
290 (n (position old simplex
)))
292 (setf (nth n simplex
) new
)
293 (send self
:sort-simplex
))))
294 (defmeth simplex-proto
:mean-opposite-face
(x)
295 (let ((face (mapcar #'(lambda (x) (send self
:point-location x
))
296 (remove x
(send self
:simplex
)))))
297 (/ (apply #'+ face
) (length face
))))
300 ;;; Iteration Step Methods
303 (defmeth simplex-proto
:extrapolate-from-worst
(fac)
304 (let* ((worst (send self
:worst-point
))
305 (wloc (send self
:point-location worst
))
306 (delta (- (send self
:mean-opposite-face worst
) wloc
))
307 (new (send self
:make-point
(+ wloc
(* (- 1 fac
) delta
)))))
308 (if (send self
:is-worse worst new
) (send self
:replace-point new worst
))
311 (defmeth simplex-proto
:shrink-to-best
(fac)
312 (let* ((best (send self
:best-point
))
313 (bloc (send self
:point-location best
)))
314 (dolist (x (copy-list (send self
:simplex
)))
315 (if (not (eq x best
))
316 (send self
:replace-point
317 (send self
:make-point
320 (- (send self
:point-location x
) bloc
))))
323 (defmeth simplex-proto
:relative-range
()
324 (let ((best (send self
:point-value
(send self
:best-point
)))
325 (worst (send self
:point-value
(send self
:worst-point
))))
326 (* 2 (/ (abs (- best worst
)) (+ 1 (abs best
) (abs worst
))))))