2 ;;; Copyright (c) 2005--2007, by A.J. Rossini <blindglobe@gmail.com>
3 ;;; See COPYRIGHT file for any additional restrictions (BSD license).
4 ;;; Since 1991, ANSI was finally finished. Edited for ANSI Common Lisp.
6 (defpackage :lisp-stat-optimize
8 :lisp-stat-object-system
10 (:shadowing-import-from
:lisp-stat-object-system
11 slot-value call-method call-next-method
)
12 (:export newtonmax nelmeadmax
))
16 ;;;; Mode Info Prototype
19 (defproto minfo-proto
'(internals))
21 #+xlisp
(send minfo-proto
:add-method
:isnew
#'|minfo-isnew|
)
22 #+xlisp
(send minfo-proto
:add-method
:maximize
#'|minfo-maximize|
)
23 #+xlisp
(send minfo-proto
:add-method
:loglaplace
#'|minfo-loglap|
)
25 (defmeth minfo-proto
:isnew
(&rest args
)
26 (setf (slot-value 'internals
) (apply #'new-minfo-internals args
)))
28 (defmeth minfo-proto
:maximize
(&rest args
)
29 (apply #'minfo-maximize
(slot-value 'internals
) args
))
31 (defmeth minfo-proto
:x
() (aref (slot-value 'internals
) 3))
32 (defmeth minfo-proto
:scale
() (aref (slot-value 'internals
) 4))
33 (defmeth minfo-proto
:derivstep
() (aref (aref (slot-value 'internals
) 9) 1))
34 (defmeth minfo-proto
:tilt
() (aref (aref (slot-value 'internals
) 9) 6))
36 (defmeth minfo-proto
:f
(&optional
(val nil set
))
38 (send self
:set-no-vals-supplied
)
39 (setf (aref (slot-value 'internals
) 0) val
))
40 (aref (slot-value 'internals
) 0))
42 (defmeth minfo-proto
:set-no-vals-supplied
()
43 (setf (aref (aref (slot-value 'internals
) 8) 6) 0))
45 (defmeth minfo-proto
:exptilt
(&optional
(val nil set
))
47 (let ((old (send self
:exptilt
)))
48 (setf (aref (aref (slot-value 'internals
) 8) 7) (if val
1 0))
49 (if (and (not (or (and old val
) (and (not old
) (not val
))))
50 (/= (send self
:tilt
) 0.0))
51 (send self
:set-no-vals-supplied
))))
52 (= 1 (aref (aref (slot-value 'internals
) 8) 7)))
54 (defmeth minfo-proto
:newtilt
(&optional
(val nil set
))
56 (setf (aref (aref (slot-value 'internals
) 9) 7) (float val
))
57 (if (/= (send self
:tilt
) 0.0) (send self
:set-no-vals-supplied
)))
58 (aref (aref (slot-value 'internals
) 9) 7))
60 (defmeth minfo-proto
:gfuns
(&optional
(val nil set
))
62 (if (or (not (consp val
))
63 (not (every #'functionp val
)))
64 (error "not all functions"))
65 (setf (aref (slot-value 'internals
) 1) val
)
66 (setf (aref (aref (slot-value 'internals
) 8) 1) (length val
))
67 (setf (aref (slot-value 'internals
) 10) (repeat 1.0 (length val
)))
68 (if (/= (send self
:tilt
) 0.0) (send self
:set-no-vals-supplied
)))
69 (aref (slot-value 'internals
) 1))
71 (defmeth minfo-proto
:cfuns
(&optional
(val nil set
))
73 (if (or (not (consp val
))
74 (not (every #'functionp val
)))
75 (error "not all functions"))
76 (setf (aref (slot-value 'internals
) 2) val
)
77 (setf (aref (aref (slot-value 'internals
) 8) 2) (length val
))
78 (setf (aref (slot-value 'internals
) 7) (repeat 0.0 (length val
)))
79 (setf (aref (slot-value 'internals
) 11) (repeat 0.0 (length val
)))
80 (send self
:set-no-vals-supplied
))
81 (aref (slot-value 'internals
) 2))
83 (defmeth minfo-proto
:ctarget
(&optional
(val nil set
))
85 (if (/= (length val
) (length (send self
:ctarget
)))
86 (error "bad target length"))
87 (setf (aref (slot-value 'internals
) 7) val
))
88 (aref (slot-value 'internals
) 7))
90 (defmeth minfo-proto
:fvals
()
91 (let* ((fv (aref (slot-value 'internals
) 5))
92 (n (length (send self
:x
)))
94 (grad (select fv
(iseq 1 n
)))
95 (hess (matrix (list n n
) (select fv
(iseq (+ 1 n
) (+ n
(* n n
)))))))
96 (list val grad hess
)))
98 (defmeth minfo-proto
:copy
()
99 (let ((obj (make-object minfo-proto
))
100 (internals (copy-seq (slot-value 'internals
))))
101 (dotimes (i (length internals
))
102 (let ((x (aref internals i
)))
104 (setf (aref internals i
) (copy-seq x
)))))
105 (send obj
:add-slot
'internals internals
)
108 (defmeth minfo-proto
:derivscale
()
109 (let* ((step (^ machine-epsilon
(/ 1 6)))
110 (hess (numhess (send self
:f
) (send self
:x
) (send self
:scale
) step
))
111 (scale (pmax (abs (send self
:x
)) (sqrt (abs (/ (diagonal hess
)))))))
112 (setf hess
(numhess (send self
:f
) (send self
:x
) scale step
))
113 (setf scale
(pmax (abs (send self
:x
)) (sqrt (abs (/ (diagonal hess
))))))
114 (setf (aref (slot-value 'internals
) 4) scale
)
115 (setf (aref (aref (slot-value 'internals
) 9) 1) step
)))
117 (defmeth minfo-proto
:verbose
(&optional
(val nil set
))
119 (setf (aref (aref (slot-value 'internals
) 8) 5)
120 (cond ((integerp val
) val
)
123 (aref (aref (slot-value 'internals
) 8) 5))
125 (defmeth minfo-proto
:backtrack
(&optional
(val nil set
))
126 (if set
(setf (aref (aref (slot-value 'internals
) 8) 4) (if val
1 0)))
127 (aref (aref (slot-value 'internals
) 8) 4))
129 (defmeth minfo-proto
:maxiter
(&optional
(val nil set
))
130 (if set
(setf (aref (aref (slot-value 'internals
) 8) 3)
131 (if (integerp val
) val -
1)))
132 (aref (aref (slot-value 'internals
) 8) 3))
134 (defmeth minfo-proto
:tiltscale
(&optional
(val nil set
))
136 (if (/= (length val
) (length (send self
:gfuns
)))
137 (error "wrong size tilt scale sequence"))
138 (setf (aref (slot-value 'internals
) 10) val
))
139 (aref (slot-value 'internals
) 10))
143 ;;;; Newton's Method with Backtracking
147 (defun newtonmax (f start
&key
153 "Args:(f start &key scale derivstep (verbose 1) return-derivs)
154 Maximizes F starting from START using Newton's method with backtracking.
155 If RETURN-DERIVS is NIL returns location of maximum; otherwise returns
156 list of location, unction value, gradient and hessian at maximum.
157 SCALE should be a list of the typical magnitudes of the parameters.
158 DERIVSTEP is used in numerical derivatives and VERBOSE controls printing
159 of iteration information. COUNT-LIMIT limits the number of iterations"
160 (let ((verbose (if verbose
(if (integerp verbose
) verbose
1) 0))
161 (minfo (send minfo-proto
:new f start
162 :scale scale
:derivstep derivstep
)))
163 (send minfo
:maxiter count-limit
)
164 (send minfo
:derivscale
)
165 (send minfo
:maximize verbose
)
167 (cons (send minfo
:x
) (- (send minfo
:fvals
)))
172 ;;;; Nelder-Mead Simplex Method
176 (defun nelmeadmax (f start
&key
178 (epsilon (sqrt machine-epsilon
))
185 "Args: (f start &key (size 1) (epsilon (sqrt machine-epsilon))
186 (count-limit 500) (verbose t) alpha beta gamma delta)
187 Maximizes F using the Nelder-Mead simplex method. START can be a
188 starting simplex - a list of N+1 points, with N=dimension of problem,
189 or a single point. If start is a single point you should give the
190 size of the initial simplex as SIZE, a sequence of length N. Default is
191 all 1's. EPSILON is the convergence tolerance. ALPHA-DELTA can be used to
192 control the behavior of simplex algorithm."
193 (let ((s (send simplex-proto
:new f start size
)))
194 (do ((best (send s
:best-point
) (send s
:best-point
))
195 (count 0 (+ count
1))
197 ((or (< (send s
:relative-range
) epsilon
) (>= count count-limit
))
198 (if (and verbose
(>= count count-limit
))
199 (format t
"Iteration limit exceeded.~%"))
200 (send s
:point-location
(send s
:best-point
)))
201 (setf next
(send s
:extrapolate-from-worst
(- alpha
)))
202 (if (send s
:is-worse best next
)
203 (setf next
(send s
:extrapolate-from-worst gamma
))
204 (when (send s
:is-worse next
(send s
:second-worst-point
))
205 (setf next
(send s
:extrapolate-from-worst beta
))
206 (if (send s
:is-worse next
(send s
:worst-point
))
207 (send s
:shrink-to-best delta
))))
209 (format t
"Value = ~10g~%"
210 (send s
:point-value
(send s
:best-point
)))))))
214 ;;; Simplex Prototype
217 (defproto simplex-proto
'(f simplex
))
223 (defmeth simplex-proto
:make-point
(x)
224 (let ((f (send self
:f
)))
226 (let ((val (funcall f x
)))
227 (cons (if (consp val
) (car val
) val
) x
))
230 (defmeth simplex-proto
:point-value
(x) (car x
))
232 (defmeth simplex-proto
:point-location
(x) (cdr x
))
234 (defmeth simplex-proto
:is-worse
(x y
)
235 (< (send self
:point-value x
) (send self
:point-value y
)))
238 ;;; Making New Simplices
241 (defmeth simplex-proto
:isnew
(f start
&optional size
)
242 (send self
:simplex start size
)
246 ;;; Slot Accessors and Mutators
249 (defmeth simplex-proto
:simplex
(&optional new size
)
252 (if (and (consp new
) (sequencep (car new
)))
253 (if (/= (length new
) (+ 1 (length (car new
))))
254 (error "bad simplex data")
256 (let* ((n (length new
))
257 (size (if size size
(repeat 1 n
)))
258 ; (pts (- (* 2 (uniform-rand (repeat n (+ n 1)))) 1)))
259 (diag (* 2 size
(- (random (repeat 2 n
)) .5)))
260 (pts (cons (repeat 0 n
)
261 (mapcar #'(lambda (x) (coerce x
'list
))
262 (column-list (diagonal diag
))))))
263 (mapcar #'(lambda (x) (+ (* size x
) new
)) pts
)))))
264 (setf (slot-value 'simplex
)
265 (mapcar #'(lambda (x) (send self
:make-point x
)) simplex
))
266 (send self
:sort-simplex
)))
267 (slot-value 'simplex
))
269 (defmeth simplex-proto
:f
(&optional f
)
271 (setf (slot-value 'f
) f
)
273 (mapcar #'(lambda (x) (send self
:point-location x
))
274 (send self
:simplex
))))
275 (send self
:simplex simplex
)))
278 (defmeth simplex-proto
:sort-simplex
()
280 (setf (slot-value 'simplex
)
281 (sort (slot-value 'simplex
)
282 #'(lambda (x y
) (send self
:is-worse x y
))))))
285 ;;; Other Methods Using List Representation of SImplex
288 (defmeth simplex-proto
:best-point
() (car (last (send self
:simplex
))))
289 (defmeth simplex-proto
:worst-point
() (first (send self
:simplex
)))
290 (defmeth simplex-proto
:second-worst-point
() (second (send self
:simplex
)))
291 (defmeth simplex-proto
:replace-point
(new old
)
292 (let* ((simplex (send self
:simplex
))
293 (n (position old simplex
)))
295 (setf (nth n simplex
) new
)
296 (send self
:sort-simplex
))))
297 (defmeth simplex-proto
:mean-opposite-face
(x)
298 (let ((face (mapcar #'(lambda (x) (send self
:point-location x
))
299 (remove x
(send self
:simplex
)))))
300 (/ (apply #'+ face
) (length face
))))
303 ;;; Iteration Step Methods
306 (defmeth simplex-proto
:extrapolate-from-worst
(fac)
307 (let* ((worst (send self
:worst-point
))
308 (wloc (send self
:point-location worst
))
309 (delta (- (send self
:mean-opposite-face worst
) wloc
))
310 (new (send self
:make-point
(+ wloc
(* (- 1 fac
) delta
)))))
311 (if (send self
:is-worse worst new
) (send self
:replace-point new worst
))
314 (defmeth simplex-proto
:shrink-to-best
(fac)
315 (let* ((best (send self
:best-point
))
316 (bloc (send self
:point-location best
)))
317 (dolist (x (copy-list (send self
:simplex
)))
318 (if (not (eq x best
))
319 (send self
:replace-point
320 (send self
:make-point
323 (- (send self
:point-location x
) bloc
))))
326 (defmeth simplex-proto
:relative-range
()
327 (let ((best (send self
:point-value
(send self
:best-point
)))
328 (worst (send self
:point-value
(send self
:worst-point
))))
329 (* 2 (/ (abs (- best worst
)) (+ 1 (abs best
) (abs worst
))))))