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.
8 (defpackage :lisp-stat-optimize
10 :lisp-stat-object-system
12 (:shadowing-import-from
:lisp-stat-object-system
13 slot-value call-method call-next-method
)
14 (:export newtonmax nelmeadmax
))
16 ;;; FIXME:AJR: There is a need to figure out the proper symbols to
17 ;;; export. more importantly should there be any specialty packages
18 ;;; that are exported for maximization?
20 (in-package :lisp-stat-optimize
)
23 ;;;; minfo basics (internal??)
26 (defun init-minfo-ipar-values (n ipars
)
38 (setf (aref ipars
0) n
)
39 (setf (aref ipars
1) m
)
40 (setf (aref ipars
2) k
)
41 (setf (aref ipars
3) itnlimit
)
42 (setf (aref ipars
4) backtrack
)
43 (setf (aref ipars
5) verbose
)
44 (setf (aref ipars
6) vals_suppl
)
45 (setf (aref ipars
7) exptilt
)
46 (setf (aref ipars
8) count
)
47 (setf (aref ipars
9) termcode
)))
49 (defun init-minfo-dpar-values (h dpars
)
58 (setf (aref dpars
0) typf
)
59 (setf (aref dpars
1) h
)
60 (setf (aref dpars
2) gradtol
)
61 (setf (aref dpars
3) steptol
)
62 (setf (aref dpars
4) maxstep
)
63 (setf (aref dpars
5) dflt
)
64 (setf (aref dpars
6) tilt
)
65 (setf (aref dpars
7) newtilt
)
66 (setf (aref dpars
8) hessadd
)))
68 (defun init-minfo-internals (n h internals
)
69 (let ((ipars (aref internals
8))
70 (dpars (aref internals
9)))
71 (init-minfo-ipar-values n ipars
)
72 (init-minfo-dpar-values h dpars
)))
74 (defun new-minfo-internals (f x
&key scale
((:derivstep h
) -
1.0))
80 (check-sequence scale
)
82 (if (/= n
(length scale
)) (error "scale and x not the same length")))
83 (let ((internals (make-array 12)))
84 (setf (aref internals
0) f
)
85 (setf (aref internals
3) (if (consp x
) (copy-list x
) (coerce x
'list
)))
86 (setf (aref internals
4)
87 (if scale
(copy-seq scale
) (make-array n
:initial-element
1.0)))
88 (setf (aref internals
5) (make-list (+ 1 n
(* n n
))))
89 (setf (aref internals
8) (make-array 10))
90 (setf (aref internals
9) (make-array 9))
91 (init-minfo-internals n h internals
)
94 (defun minfo-maximize (internals &optional verbose
)
95 "This function does what?"
96 (let* ((f (aref internals
0))
97 (x (aref internals
3))
98 (fvals (aref internals
5))
100 (v (if verbose
(if (integerp verbose
) verbose
1) -
1)))
101 (setf (aref internals
3) (copy-list x
))
102 (setf (aref internals
5) (copy-list fvals
))
103 (let ((*maximize-callback-function
* f
)
104 (*maximize-callback-arg
* (make-list n
)))
105 (let* ((x (aref internals
3))
106 (scale (aref internals
4))
107 (fvals (aref internals
5))
108 (ip (aref internals
8))
109 (dp (aref internals
9))
110 (px (la-data-to-vector x mode-re
))
111 (pscale (la-data-to-vector scale mode-re
))
112 (pfvals (la-vector (length fvals
) mode-re
))
113 (pip (la-data-to-vector ip mode-in
))
114 (pdp (la-data-to-vector dp mode-re
)))
117 (base-minfo-maximize px pfvals pscale pip pdp v
)) ;; access to C
118 (la-vector-to-data px n mode-re x
)
119 (la-vector-to-data pfvals
(+ 1 n
(* n n
)) mode-re fvals
)
120 (la-vector-to-data pip
(length ip
) mode-in ip
)
121 (la-vector-to-data pdp
(length dp
) mode-re dp
))
127 ;;;; Mode Info Prototype
130 (defproto minfo-proto
'(internals))
132 #+xlisp
(send minfo-proto
:add-method
:isnew
#'|minfo-isnew|
)
133 #+xlisp
(send minfo-proto
:add-method
:maximize
#'|minfo-maximize|
)
134 #+xlisp
(send minfo-proto
:add-method
:loglaplace
#'|minfo-loglap|
)
136 (defmeth minfo-proto
:isnew
(&rest args
)
137 (setf (slot-value 'internals
) (apply #'new-minfo-internals args
)))
139 (defmeth minfo-proto
:maximize
(&rest args
)
140 (apply #'minfo-maximize
(slot-value 'internals
) args
))
142 (defmeth minfo-proto
:x
() (aref (slot-value 'internals
) 3))
143 (defmeth minfo-proto
:scale
() (aref (slot-value 'internals
) 4))
144 (defmeth minfo-proto
:derivstep
() (aref (aref (slot-value 'internals
) 9) 1))
145 (defmeth minfo-proto
:tilt
() (aref (aref (slot-value 'internals
) 9) 6))
147 (defmeth minfo-proto
:f
(&optional
(val nil set
))
149 (send self
:set-no-vals-supplied
)
150 (setf (aref (slot-value 'internals
) 0) val
))
151 (aref (slot-value 'internals
) 0))
153 (defmeth minfo-proto
:set-no-vals-supplied
()
154 (setf (aref (aref (slot-value 'internals
) 8) 6) 0))
156 (defmeth minfo-proto
:exptilt
(&optional
(val nil set
))
158 (let ((old (send self
:exptilt
)))
159 (setf (aref (aref (slot-value 'internals
) 8) 7) (if val
1 0))
160 (if (and (not (or (and old val
) (and (not old
) (not val
))))
161 (/= (send self
:tilt
) 0.0))
162 (send self
:set-no-vals-supplied
))))
163 (= 1 (aref (aref (slot-value 'internals
) 8) 7)))
165 (defmeth minfo-proto
:newtilt
(&optional
(val nil set
))
167 (setf (aref (aref (slot-value 'internals
) 9) 7) (float val
))
168 (if (/= (send self
:tilt
) 0.0) (send self
:set-no-vals-supplied
)))
169 (aref (aref (slot-value 'internals
) 9) 7))
171 (defmeth minfo-proto
:gfuns
(&optional
(val nil set
))
173 (if (or (not (consp val
))
174 (not (every #'functionp val
)))
175 (error "not all functions"))
176 (setf (aref (slot-value 'internals
) 1) val
)
177 (setf (aref (aref (slot-value 'internals
) 8) 1) (length val
))
178 (setf (aref (slot-value 'internals
) 10) (repeat 1.0 (length val
)))
179 (if (/= (send self
:tilt
) 0.0) (send self
:set-no-vals-supplied
)))
180 (aref (slot-value 'internals
) 1))
182 (defmeth minfo-proto
:cfuns
(&optional
(val nil set
))
184 (if (or (not (consp val
))
185 (not (every #'functionp val
)))
186 (error "not all functions"))
187 (setf (aref (slot-value 'internals
) 2) val
)
188 (setf (aref (aref (slot-value 'internals
) 8) 2) (length val
))
189 (setf (aref (slot-value 'internals
) 7) (repeat 0.0 (length val
)))
190 (setf (aref (slot-value 'internals
) 11) (repeat 0.0 (length val
)))
191 (send self
:set-no-vals-supplied
))
192 (aref (slot-value 'internals
) 2))
194 (defmeth minfo-proto
:ctarget
(&optional
(val nil set
))
196 (if (/= (length val
) (length (send self
:ctarget
)))
197 (error "bad target length"))
198 (setf (aref (slot-value 'internals
) 7) val
))
199 (aref (slot-value 'internals
) 7))
201 (defmeth minfo-proto
:fvals
()
202 (let* ((fv (aref (slot-value 'internals
) 5))
203 (n (length (send self
:x
)))
205 (grad (select fv
(iseq 1 n
)))
206 (hess (matrix (list n n
) (select fv
(iseq (+ 1 n
) (+ n
(* n n
)))))))
207 (list val grad hess
)))
209 (defmeth minfo-proto
:copy
()
210 (let ((obj (make-object minfo-proto
))
211 (internals (copy-seq (slot-value 'internals
))))
212 (dotimes (i (length internals
))
213 (let ((x (aref internals i
)))
215 (setf (aref internals i
) (copy-seq x
)))))
216 (send obj
:add-slot
'internals internals
)
219 (defmeth minfo-proto
:derivscale
()
220 (let* ((step (^ machine-epsilon
(/ 1 6)))
221 (hess (numhess (send self
:f
) (send self
:x
) (send self
:scale
) step
))
222 (scale (pmax (abs (send self
:x
)) (sqrt (abs (/ (diagonal hess
)))))))
223 (setf hess
(numhess (send self
:f
) (send self
:x
) scale step
))
224 (setf scale
(pmax (abs (send self
:x
)) (sqrt (abs (/ (diagonal hess
))))))
225 (setf (aref (slot-value 'internals
) 4) scale
)
226 (setf (aref (aref (slot-value 'internals
) 9) 1) step
)))
228 (defmeth minfo-proto
:verbose
(&optional
(val nil set
))
230 (setf (aref (aref (slot-value 'internals
) 8) 5)
231 (cond ((integerp val
) val
)
234 (aref (aref (slot-value 'internals
) 8) 5))
236 (defmeth minfo-proto
:backtrack
(&optional
(val nil set
))
237 (if set
(setf (aref (aref (slot-value 'internals
) 8) 4) (if val
1 0)))
238 (aref (aref (slot-value 'internals
) 8) 4))
240 (defmeth minfo-proto
:maxiter
(&optional
(val nil set
))
241 (if set
(setf (aref (aref (slot-value 'internals
) 8) 3)
242 (if (integerp val
) val -
1)))
243 (aref (aref (slot-value 'internals
) 8) 3))
245 (defmeth minfo-proto
:tiltscale
(&optional
(val nil set
))
247 (if (/= (length val
) (length (send self
:gfuns
)))
248 (error "wrong size tilt scale sequence"))
249 (setf (aref (slot-value 'internals
) 10) val
))
250 (aref (slot-value 'internals
) 10))
254 ;;;; Newton's Method with Backtracking
258 (defun newtonmax (f start
&key
264 "Args:(f start &key scale derivstep (verbose 1) return-derivs)
265 Maximizes F starting from START using Newton's method with backtracking.
266 If RETURN-DERIVS is NIL returns location of maximum; otherwise returns
267 list of location, unction value, gradient and hessian at maximum.
268 SCALE should be a list of the typical magnitudes of the parameters.
269 DERIVSTEP is used in numerical derivatives and VERBOSE controls printing
270 of iteration information. COUNT-LIMIT limits the number of iterations"
271 (let ((verbose (if verbose
(if (integerp verbose
) verbose
1) 0))
272 (minfo (send minfo-proto
:new f start
273 :scale scale
:derivstep derivstep
)))
274 (send minfo
:maxiter count-limit
)
275 (send minfo
:derivscale
)
276 (send minfo
:maximize verbose
)
278 (cons (send minfo
:x
) (- (send minfo
:fvals
)))
283 ;;;; Nelder-Mead Simplex Method
287 (defun nelmeadmax (f start
&key
289 (epsilon (sqrt machine-epsilon
))
296 "Args: (f start &key (size 1) (epsilon (sqrt machine-epsilon))
297 (count-limit 500) (verbose t) alpha beta gamma delta)
298 Maximizes F using the Nelder-Mead simplex method. START can be a
299 starting simplex - a list of N+1 points, with N=dimension of problem,
300 or a single point. If start is a single point you should give the
301 size of the initial simplex as SIZE, a sequence of length N. Default is
302 all 1's. EPSILON is the convergence tolerance. ALPHA-DELTA can be used to
303 control the behavior of simplex algorithm."
304 (let ((s (send simplex-proto
:new f start size
)))
305 (do ((best (send s
:best-point
) (send s
:best-point
))
306 (count 0 (+ count
1))
308 ((or (< (send s
:relative-range
) epsilon
) (>= count count-limit
))
309 (if (and verbose
(>= count count-limit
))
310 (format t
"Iteration limit exceeded.~%"))
311 (send s
:point-location
(send s
:best-point
)))
312 (setf next
(send s
:extrapolate-from-worst
(- alpha
)))
313 (if (send s
:is-worse best next
)
314 (setf next
(send s
:extrapolate-from-worst gamma
))
315 (when (send s
:is-worse next
(send s
:second-worst-point
))
316 (setf next
(send s
:extrapolate-from-worst beta
))
317 (if (send s
:is-worse next
(send s
:worst-point
))
318 (send s
:shrink-to-best delta
))))
320 (format t
"Value = ~10g~%"
321 (send s
:point-value
(send s
:best-point
)))))))
325 ;;; Simplex Prototype
328 (defproto simplex-proto
'(f simplex
))
334 (defmeth simplex-proto
:make-point
(x)
335 (let ((f (send self
:f
)))
337 (let ((val (funcall f x
)))
338 (cons (if (consp val
) (car val
) val
) x
))
341 (defmeth simplex-proto
:point-value
(x) (car x
))
343 (defmeth simplex-proto
:point-location
(x) (cdr x
))
345 (defmeth simplex-proto
:is-worse
(x y
)
346 (< (send self
:point-value x
) (send self
:point-value y
)))
349 ;;; Making New Simplices
352 (defmeth simplex-proto
:isnew
(f start
&optional size
)
353 (send self
:simplex start size
)
357 ;;; Slot Accessors and Mutators
360 (defmeth simplex-proto
:simplex
(&optional new size
)
363 (if (and (consp new
) (sequencep (car new
)))
364 (if (/= (length new
) (+ 1 (length (car new
))))
365 (error "bad simplex data")
367 (let* ((n (length new
))
368 (size (if size size
(repeat 1 n
)))
369 ; (pts (- (* 2 (uniform-rand (repeat n (+ n 1)))) 1)))
370 (diag (* 2 size
(- (random (repeat 2 n
)) .5)))
371 (pts (cons (repeat 0 n
)
372 (mapcar #'(lambda (x) (coerce x
'list
))
373 (column-list (diagonal diag
))))))
374 (mapcar #'(lambda (x) (+ (* size x
) new
)) pts
)))))
375 (setf (slot-value 'simplex
)
376 (mapcar #'(lambda (x) (send self
:make-point x
)) simplex
))
377 (send self
:sort-simplex
)))
378 (slot-value 'simplex
))
380 (defmeth simplex-proto
:f
(&optional f
)
382 (setf (slot-value 'f
) f
)
384 (mapcar #'(lambda (x) (send self
:point-location x
))
385 (send self
:simplex
))))
386 (send self
:simplex simplex
)))
389 (defmeth simplex-proto
:sort-simplex
()
391 (setf (slot-value 'simplex
)
392 (sort (slot-value 'simplex
)
393 #'(lambda (x y
) (send self
:is-worse x y
))))))
396 ;;; Other Methods Using List Representation of SImplex
399 (defmeth simplex-proto
:best-point
() (car (last (send self
:simplex
))))
400 (defmeth simplex-proto
:worst-point
() (first (send self
:simplex
)))
401 (defmeth simplex-proto
:second-worst-point
() (second (send self
:simplex
)))
402 (defmeth simplex-proto
:replace-point
(new old
)
403 (let* ((simplex (send self
:simplex
))
404 (n (position old simplex
)))
406 (setf (nth n simplex
) new
)
407 (send self
:sort-simplex
))))
408 (defmeth simplex-proto
:mean-opposite-face
(x)
409 (let ((face (mapcar #'(lambda (x) (send self
:point-location x
))
410 (remove x
(send self
:simplex
)))))
411 (/ (reduce #'+ face
) (length face
))))
414 ;;; Iteration Step Methods
417 (defmeth simplex-proto
:extrapolate-from-worst
(fac)
418 (let* ((worst (send self
:worst-point
))
419 (wloc (send self
:point-location worst
))
420 (delta (- (send self
:mean-opposite-face worst
) wloc
))
421 (new (send self
:make-point
(+ wloc
(* (- 1 fac
) delta
)))))
422 (if (send self
:is-worse worst new
) (send self
:replace-point new worst
))
425 (defmeth simplex-proto
:shrink-to-best
(fac)
426 (let* ((best (send self
:best-point
))
427 (bloc (send self
:point-location best
)))
428 (dolist (x (copy-list (send self
:simplex
)))
429 (if (not (eq x best
))
430 (send self
:replace-point
431 (send self
:make-point
434 (- (send self
:point-location x
) bloc
))))
437 (defmeth simplex-proto
:relative-range
()
438 (let ((best (send self
:point-value
(send self
:best-point
)))
439 (worst (send self
:point-value
(send self
:worst-point
))))
440 (* 2 (/ (abs (- best worst
)) (+ 1 (abs best
) (abs worst
))))))