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
12 :lisp-stat-object-system
14 :lisp-stat-compound-data
19 :lisp-stat-linalg-data
21 (:shadowing-import-from
:lisp-stat-object-system
22 slot-value call-method call-next-method
)
23 (:shadowing-import-from
:lisp-stat-math
24 expt
+ -
* / ** mod rem abs
1+ 1- log exp sqrt sin cos tan
25 asin acos atan sinh cosh tanh asinh acosh atanh float random
26 truncate floor ceiling round minusp zerop plusp evenp oddp
27 < <= = /= >= > complex conjugate realpart imagpart phase
28 min max logand logior logxor lognot ffloor fceiling
29 ftruncate fround signum cis
)
35 newtonmax nelmeadmax
))
37 ;; matrix is in statistics, but should that be a predecessor?
39 ;;; FIXME:AJR: There is a need to figure out the proper symbols to
40 ;;; export. more importantly should there be any specialty package
41 ;;; that are exported for maximization?
43 (in-package :lisp-stat-optimize
)
46 (defctype size-t
:unsigned-long
)
48 (defctype size-t
:unsigned-int
)
50 (defvar *maximize-callback-function
* nil
51 "Used in generic optimization to determine function name -- symbol or string?")
53 (defvar *maximize-callback-arg
* nil
54 "args to function to maximize")
58 ;;; CFFI support using library for optimization work.
61 ;; There is a problem with this particular approach, in terms of
62 ;; circular dependencies. We can not have this out-of-object call
63 ;; into optimize, at least not from here.
64 (cffi:defcallback ccl-maximize-callback
:void
((n :int
)
70 (lisp-stat-optimize::maximize-callback n px pfval pgrad phess pderivs
))
72 (cffi:defcfun
("register_maximize_callback" register-maximize-callback
)
74 (register-maximize-callback (cffi:callback ccl-maximize-callback
))
76 (cffi:defcfun
("ccl_numgrad_front" ccl-numgrad-front
)
77 :void
(x size-t
) (y :pointer
) (z :pointer
) (u :double
) (v :pointer
))
78 (defun numgrad-front (x y z u v
)
79 (ccl-numgrad-front x y z
(float u
1d0
) v
))
81 (cffi:defcfun
("ccl_numhess_front" ccl-numhess-front
)
82 :void
(x size-t
) (y :pointer
) (z :pointer
) (u :pointer
) (v :pointer
) (w :double
) (a :pointer
))
83 (defun numhess-front (x y z u v w a
)
84 (ccl-numhess-front x y z u v
(float w
1d0
) a
))
86 (cffi:defcfun
("ccl_minfo_maximize" ccl-minfo-maximize
)
87 :void
(x :pointer
) (y :pointer
) (z :pointer
) (u :pointer
) (v :pointer
) (w :int
))
88 (defun base-minfo-maximize (x y z u v w
)
89 (ccl-minfo-maximize x y z u v w
))
94 ;;;; minfo basics (internal??)
97 (defun init-minfo-ipar-values (n ipars
&key
109 "Initialize ipars (iteration parameters) by destructive modification."
110 (setf (aref ipars
0) n
)
111 (setf (aref ipars
1) m
)
112 (setf (aref ipars
2) k
)
113 (setf (aref ipars
3) itnlimit
)
114 (setf (aref ipars
4) backtrack
)
115 (setf (aref ipars
5) verbose
)
116 (setf (aref ipars
6) vals_suppl
)
117 (setf (aref ipars
7) exptilt
)
118 (setf (aref ipars
8) count
)
119 (setf (aref ipars
9) termcode
))
121 (defun init-minfo-dpar-values (h dpars
&key
130 "Initialize dpars (derivative parameters) by destructive modification."
131 (setf (aref dpars
0) typf
)
132 (setf (aref dpars
1) h
)
133 (setf (aref dpars
2) gradtol
)
134 (setf (aref dpars
3) steptol
)
135 (setf (aref dpars
4) maxstep
)
136 (setf (aref dpars
5) dflt
)
137 (setf (aref dpars
6) tilt
)
138 (setf (aref dpars
7) newtilt
)
139 (setf (aref dpars
8) hessadd
))
141 (defun init-minfo-internals (n h internals
)
142 (let ((ipars (aref internals
8))
143 (dpars (aref internals
9)))
144 (init-minfo-ipar-values n ipars
)
145 (init-minfo-dpar-values h dpars
)))
147 (defun new-minfo-internals (f x
&key scale
((:derivstep h
) -
1.0))
151 (let ((n (length x
)))
153 (check-sequence scale
)
155 (if (/= n
(length scale
)) (error "scale and x not the same length")))
156 (let ((internals (make-array 12)))
157 (setf (aref internals
0) f
)
158 (setf (aref internals
3) (if (consp x
) (copy-list x
) (coerce x
'list
)))
159 (setf (aref internals
4)
160 (if scale
(copy-seq scale
) (make-array n
:initial-element
1.0)))
161 (setf (aref internals
5) (make-list (+ 1 n
(* n n
))))
162 (setf (aref internals
8) (make-array 10))
163 (setf (aref internals
9) (make-array 9))
164 (init-minfo-internals n h internals
)
167 (defun minfo-maximize (internals &optional verbose
)
168 "This function does what?"
169 (let* ((f (aref internals
0))
170 (x (aref internals
3))
171 (fvals (aref internals
5))
173 (v (if verbose
(if (integerp verbose
) verbose
1) -
1)))
174 (setf (aref internals
3) (copy-list x
))
175 (setf (aref internals
5) (copy-list fvals
))
176 (let ((*maximize-callback-function
* f
)
177 (*maximize-callback-arg
* (make-list n
)))
178 (let* ((x (aref internals
3))
179 (scale (aref internals
4))
180 (fvals (aref internals
5))
181 (ip (aref internals
8))
182 (dp (aref internals
9))
183 (px (la-data-to-vector x
+mode-re
+))
184 (pscale (la-data-to-vector scale
+mode-re
+))
185 (pfvals (la-vector (length fvals
) +mode-re
+))
186 (pip (la-data-to-vector ip
+mode-in
+))
187 (pdp (la-data-to-vector dp
+mode-re
+)))
190 (base-minfo-maximize px pfvals pscale pip pdp v
)) ;; access to C
191 (la-vector-to-data px n
+mode-re
+ x
)
192 (la-vector-to-data pfvals
(+ 1 n
(* n n
)) +mode-re
+ fvals
)
193 (la-vector-to-data pip
(length ip
) +mode-in
+ ip
)
194 (la-vector-to-data pdp
(length dp
) +mode-re
+ dp
))
200 ;;;; Mode Info Prototype
204 (defproto minfo-proto
'(internals))
206 #+xlisp
(send minfo-proto
:add-method
:isnew
#'|minfo-isnew|
)
207 #+xlisp
(send minfo-proto
:add-method
:maximize
#'|minfo-maximize|
)
208 #+xlisp
(send minfo-proto
:add-method
:loglaplace
#'|minfo-loglap|
)
210 (defmeth minfo-proto
:isnew
(&rest args
)
211 (setf (slot-value 'internals
) (apply #'new-minfo-internals args
)))
213 (defmeth minfo-proto
:maximize
(&rest args
)
214 (apply #'minfo-maximize
(slot-value 'internals
) args
))
216 (defmeth minfo-proto
:x
() (aref (slot-value 'internals
) 3))
217 (defmeth minfo-proto
:scale
() (aref (slot-value 'internals
) 4))
218 (defmeth minfo-proto
:derivstep
() (aref (aref (slot-value 'internals
) 9) 1))
219 (defmeth minfo-proto
:tilt
() (aref (aref (slot-value 'internals
) 9) 6))
221 (defmeth minfo-proto
:f
(&optional
(val nil set
))
223 (send self
:set-no-vals-supplied
)
224 (setf (aref (slot-value 'internals
) 0) val
))
225 (aref (slot-value 'internals
) 0))
227 (defmeth minfo-proto
:set-no-vals-supplied
()
228 (setf (aref (aref (slot-value 'internals
) 8) 6) 0))
230 (defmeth minfo-proto
:exptilt
(&optional
(val nil set
))
232 (let ((old (send self
:exptilt
)))
233 (setf (aref (aref (slot-value 'internals
) 8) 7) (if val
1 0))
234 (if (and (not (or (and old val
) (and (not old
) (not val
))))
235 (/= (send self
:tilt
) 0.0))
236 (send self
:set-no-vals-supplied
))))
237 (= 1 (aref (aref (slot-value 'internals
) 8) 7)))
239 (defmeth minfo-proto
:newtilt
(&optional
(val nil set
))
241 (setf (aref (aref (slot-value 'internals
) 9) 7) (float val
))
242 (if (/= (send self
:tilt
) 0.0) (send self
:set-no-vals-supplied
)))
243 (aref (aref (slot-value 'internals
) 9) 7))
245 (defmeth minfo-proto
:gfuns
(&optional
(val nil set
))
247 (if (or (not (consp val
))
248 (not (every #'functionp val
)))
249 (error "not all functions"))
250 (setf (aref (slot-value 'internals
) 1) val
)
251 (setf (aref (aref (slot-value 'internals
) 8) 1) (length val
))
252 (setf (aref (slot-value 'internals
) 10) (repeat 1.0 (length val
)))
253 (if (/= (send self
:tilt
) 0.0) (send self
:set-no-vals-supplied
)))
254 (aref (slot-value 'internals
) 1))
256 (defmeth minfo-proto
:cfuns
(&optional
(val nil set
))
258 (if (or (not (consp val
))
259 (not (every #'functionp val
)))
260 (error "not all functions"))
261 (setf (aref (slot-value 'internals
) 2) val
)
262 (setf (aref (aref (slot-value 'internals
) 8) 2) (length val
))
263 (setf (aref (slot-value 'internals
) 7) (repeat 0.0 (length val
)))
264 (setf (aref (slot-value 'internals
) 11) (repeat 0.0 (length val
)))
265 (send self
:set-no-vals-supplied
))
266 (aref (slot-value 'internals
) 2))
268 (defmeth minfo-proto
:ctarget
(&optional
(val nil set
))
270 (if (/= (length val
) (length (send self
:ctarget
)))
271 (error "bad target length"))
272 (setf (aref (slot-value 'internals
) 7) val
))
273 (aref (slot-value 'internals
) 7))
275 (defmeth minfo-proto
:fvals
()
276 (let* ((fv (aref (slot-value 'internals
) 5))
277 (n (length (send self
:x
)))
279 (grad (select fv
(iseq 1 n
)))
280 (hess (matrix (list n n
) (select fv
(iseq (+ 1 n
) (+ n
(* n n
)))))))
281 (list val grad hess
)))
283 (defmeth minfo-proto
:copy
()
286 Make a copy of an minfo instance."
287 (let ((obj (make-object minfo-proto
))
288 (internals (copy-seq (slot-value 'internals
))))
289 (dotimes (i (length internals
))
290 (let ((x (aref internals i
)))
291 (if (typep x
'sequence
)
292 (setf (aref internals i
) (copy-seq x
)))))
293 (send obj
:add-slot
'internals internals
)
296 (defmeth minfo-proto
:derivscale
()
297 (let* ((step (^ machine-epsilon
(/ 1 6)))
298 (hess (numhess (send self
:f
) (send self
:x
) (send self
:scale
) step
))
299 (scale (pmax (abs (send self
:x
)) (sqrt (abs (/ (diagonal hess
)))))))
300 (setf hess
(numhess (send self
:f
) (send self
:x
) scale step
))
301 (setf scale
(pmax (abs (send self
:x
)) (sqrt (abs (/ (diagonal hess
))))))
302 (setf (aref (slot-value 'internals
) 4) scale
)
303 (setf (aref (aref (slot-value 'internals
) 9) 1) step
)))
305 (defmeth minfo-proto
:verbose
(&optional
(val nil set
))
307 (setf (aref (aref (slot-value 'internals
) 8) 5)
308 (cond ((integerp val
) val
)
311 (aref (aref (slot-value 'internals
) 8) 5))
313 (defmeth minfo-proto
:backtrack
(&optional
(val nil set
))
314 (if set
(setf (aref (aref (slot-value 'internals
) 8) 4) (if val
1 0)))
315 (aref (aref (slot-value 'internals
) 8) 4))
317 (defmeth minfo-proto
:maxiter
(&optional
(val nil set
))
318 (if set
(setf (aref (aref (slot-value 'internals
) 8) 3)
319 (if (integerp val
) val -
1)))
320 (aref (aref (slot-value 'internals
) 8) 3))
322 (defmeth minfo-proto
:tiltscale
(&optional
(val nil set
))
324 (if (/= (length val
) (length (send self
:gfuns
)))
325 (error "wrong size tilt scale sequence"))
326 (setf (aref (slot-value 'internals
) 10) val
))
327 (aref (slot-value 'internals
) 10))
331 ;;;; Newton's Method with Backtracking
335 (defun newtonmax (f start
&key
341 "Args:(f start &key scale derivstep (verbose 1) return-derivs)
342 Maximizes F starting from START using Newton's method with backtracking.
343 If RETURN-DERIVS is NIL returns location of maximum; otherwise returns
344 list of location, unction value, gradient and hessian at maximum.
345 SCALE should be a list of the typical magnitudes of the parameters.
346 DERIVSTEP is used in numerical derivatives and VERBOSE controls printing
347 of iteration information. COUNT-LIMIT limits the number of iterations"
348 (let ((verbose (if verbose
(if (integerp verbose
) verbose
1) 0))
349 (minfo (send minfo-proto
:new f start
350 :scale scale
:derivstep derivstep
)))
351 (send minfo
:maxiter count-limit
)
352 (send minfo
:derivscale
)
353 (send minfo
:maximize verbose
)
355 (cons (send minfo
:x
) (- (send minfo
:fvals
)))
359 ;;; Nelder-Mead Simplex Method
362 ;;; Simplex Prototype
364 (defvar simplex-proto
)
365 (defproto simplex-proto
'(f simplex
))
367 (defun nelmeadmax (f start
&key
369 (epsilon (sqrt machine-epsilon
))
376 "Args: (f start &key (size 1) (epsilon (sqrt machine-epsilon))
377 (count-limit 500) (verbose t) alpha beta gamma delta)
378 Maximizes F using the Nelder-Mead simplex method. START can be a
379 starting simplex - a list of N+1 points, with N=dimension of problem,
380 or a single point. If start is a single point you should give the
381 size of the initial simplex as SIZE, a sequence of length N. Default is
382 all 1's. EPSILON is the convergence tolerance. ALPHA-DELTA can be used to
383 control the behavior of simplex algorithm."
384 (let ((s (send simplex-proto
:new f start size
)))
385 (do ((best (send s
:best-point
) (send s
:best-point
))
386 (count 0 (+ count
1))
388 ((or (< (send s
:relative-range
) epsilon
) (>= count count-limit
))
389 (if (and verbose
(>= count count-limit
))
390 (format t
"Iteration limit exceeded.~%"))
391 (send s
:point-location
(send s
:best-point
)))
392 (setf next
(send s
:extrapolate-from-worst
(- alpha
)))
393 (if (send s
:is-worse best next
)
394 (setf next
(send s
:extrapolate-from-worst gamma
))
395 (when (send s
:is-worse next
(send s
:second-worst-point
))
396 (setf next
(send s
:extrapolate-from-worst beta
))
397 (if (send s
:is-worse next
(send s
:worst-point
))
398 (send s
:shrink-to-best delta
))))
400 (format t
"Value = ~10g~%"
401 (send s
:point-value
(send s
:best-point
)))))))
409 (defmeth simplex-proto
:make-point
(x)
410 (let ((f (send self
:f
)))
412 (let ((val (funcall f x
)))
413 (cons (if (consp val
) (car val
) val
) x
))
416 (defmeth simplex-proto
:point-value
(x) (car x
))
418 (defmeth simplex-proto
:point-location
(x) (cdr x
))
420 (defmeth simplex-proto
:is-worse
(x y
)
421 (< (send self
:point-value x
) (send self
:point-value y
)))
424 ;;; Making New Simplices
427 (defmeth simplex-proto
:isnew
(f start
&optional size
)
428 (send self
:simplex start size
)
432 ;;; Slot Accessors and Mutators
435 (defmeth simplex-proto
:simplex
(&optional new size
)
438 (if (and (consp new
) (typep (car new
) 'sequence
))
439 (if (/= (length new
) (+ 1 (length (car new
))))
440 (error "bad simplex data")
442 (let* ((n (length new
))
443 (size (if size size
(repeat 1 n
)))
444 ; (pts (- (* 2 (uniform-rand (repeat n (+ n 1)))) 1)))
445 (diag (* 2 size
(- (random (repeat 2 n
)) .5)))
446 (pts (cons (repeat 0 n
)
447 (mapcar #'(lambda (x) (coerce x
'list
))
448 (column-list (diagonal diag
))))))
449 (mapcar #'(lambda (x) (reduce #'+ (list (* size x
) new
))) pts
)))))
450 (setf (slot-value 'simplex
)
451 (mapcar #'(lambda (x) (send self
:make-point x
)) simplex
))
452 (send self
:sort-simplex
)))
453 (slot-value 'simplex
))
455 (defmeth simplex-proto
:f
(&optional f
)
457 (setf (slot-value 'f
) f
)
459 (mapcar #'(lambda (x) (send self
:point-location x
))
460 (send self
:simplex
))))
461 (send self
:simplex simplex
)))
464 (defmeth simplex-proto
:sort-simplex
()
466 (setf (slot-value 'simplex
)
467 (sort (slot-value 'simplex
)
468 #'(lambda (x y
) (send self
:is-worse x y
))))))
471 ;;; Other Methods Using List Representation of SImplex
474 (defmeth simplex-proto
:best-point
() (car (last (send self
:simplex
))))
475 (defmeth simplex-proto
:worst-point
() (first (send self
:simplex
)))
476 (defmeth simplex-proto
:second-worst-point
() (second (send self
:simplex
)))
477 (defmeth simplex-proto
:replace-point
(new old
)
478 (let* ((simplex (send self
:simplex
))
479 (n (position old simplex
)))
481 (setf (nth n simplex
) new
)
482 (send self
:sort-simplex
))))
483 (defmeth simplex-proto
:mean-opposite-face
(x)
484 (let ((face (mapcar #'(lambda (x) (send self
:point-location x
))
485 (remove x
(send self
:simplex
)))))
486 (/ (reduce #'+ face
) (length face
))))
489 ;;; Iteration Step Methods
492 (defmeth simplex-proto
:extrapolate-from-worst
(fac)
493 (let* ((worst (send self
:worst-point
))
494 (wloc (send self
:point-location worst
))
495 (delta (- (send self
:mean-opposite-face worst
) wloc
))
496 (new (send self
:make-point
(+ wloc
(* (- 1 fac
) delta
)))))
497 (if (send self
:is-worse worst new
) (send self
:replace-point new worst
))
500 (defmeth simplex-proto
:shrink-to-best
(fac)
501 (let* ((best (send self
:best-point
))
502 (bloc (send self
:point-location best
)))
503 (dolist (x (copy-list (send self
:simplex
)))
504 (if (not (eq x best
))
505 (send self
:replace-point
506 (send self
:make-point
509 (- (send self
:point-location x
) bloc
))))
512 (defmeth simplex-proto
:relative-range
()
513 (let ((best (send self
:point-value
(send self
:best-point
)))
514 (worst (send self
:point-value
(send self
:worst-point
))))
515 (* 2 (/ (abs (- best worst
)) (+ 1 (abs best
) (abs worst
))))))
521 ;;;; Maximization and Numerical Derivatives
525 (defun data2double (n data ptr
)
527 (let* ((seq (compound-data-seq data
))
528 (elem (make-next-element seq
)))
529 (if (/= (length seq
) n
) (error "bad data size"))
532 (la-put-double ptr i
(get-next-element elem i
)))))
534 (defun maximize-callback (n px pfval pgrad phess pderivs
)
535 (la-vector-to-data px n
+mode-re
+ *maximize-callback-arg
*)
536 (let* ((val (funcall *maximize-callback-function
* *maximize-callback-arg
*))
537 (derivs (if (consp val
) (- (length val
) 1) 0)))
538 (la-put-integer pderivs
0 derivs
)
539 (la-put-double pfval
0 (if (consp val
) (first val
) val
))
540 (if (<= 1 derivs
) (data2double n
(second val
) pgrad
))
541 (if (<= 2 derivs
) (data2double (* n n
) (third val
) phess
))))
543 (defun numgrad (f x
&optional scale
(h -
1.0))
544 "Args: (f x &optional scale derivstep)
545 Computes the numerical gradient of F at X."
549 (check-sequence scale
)
552 (let* ((n (length x
))
553 (result (make-list n
)))
554 (if (and scale
(/= n
(length scale
)))
555 (error "scale not the same length as x"))
556 (let ((*maximize-callback-function
* f
)
557 (*maximize-callback-arg
* (make-list n
)))
558 (let ((px (la-data-to-vector x
+mode-re
+))
559 (pgrad (la-vector n
+mode-re
+))
560 (pscale (la-data-to-vector
561 (if scale scale
(make-list n
:initial-element
1.0))
565 (numgrad-front n px pgrad h pscale
)
566 (la-vector-to-data pgrad n
+mode-re
+ result
))
568 (la-free-vector pgrad
)
569 (la-free-vector pscale
))))
572 (defun numhess (f x
&optional scale
(h -
1.0) all
)
573 "Args: (f x &optional scale derivstep)
574 Computes the numerical Hessian matrix of F at X."
578 (check-sequence scale
)
581 (let* ((n (length x
))
583 (list nil
(make-list n
) (make-array (list n n
)))
584 (make-array (list n n
)))))
585 (if (and scale
(/= n
(length scale
)))
586 (error "scale not the same length as x"))
587 (let ((*maximize-callback-function
* f
)
588 (*maximize-callback-arg
* (make-list n
)))
589 (let ((hess-data (compound-data-seq (if all
(third result
) result
)))
590 (px (la-data-to-vector x
+mode-re
+))
591 (pf (la-vector 1 +mode-re
+))
592 (pgrad (la-vector n
+mode-re
+))
593 (phess (la-vector (* n n
) +mode-re
+))
594 (pscale (la-data-to-vector
595 (if scale scale
(make-list n
:initial-element
1.0))
599 (numhess-front n px pf pgrad phess h pscale
)
601 (setf (first result
) (la-get-double pf
0))
602 (la-vector-to-data pgrad n
+mode-re
+ (second result
)))
603 (la-vector-to-data phess
(* n n
) +mode-re
+ hess-data
))
606 (la-free-vector pgrad
)
607 (la-free-vector phess
)
608 (la-free-vector pscale
))))