merged from ansiClib
[CommonLispStat.git] / maximize.lsp
blobf4b56beb0626fdfb917736539f8045830f60aade
1 (provide "maximize")
3 #+:CLtL2
4 (in-package lisp-stat)
5 #-:CLtL2
6 (in-package 'lisp-stat)
8 (export '(newtonmax nelmeadmax))
10 (import '(ls-basics::new-minfo-internals ls-basics::minfo-maximize))
12 ;;;;
13 ;;;; Mode Info Prototype
14 ;;;;
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|)
21 #-xlisp
22 (defmeth minfo-proto :isnew (&rest args)
23 (setf (slot-value 'internals) (apply #'new-minfo-internals args)))
24 #-xlisp
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))
34 (when 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))
43 (if 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))
52 (when 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))
58 (when 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))
69 (when 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))
81 (when 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)))
90 (val (select fv 0))
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)))
100 (if (sequencep x)
101 (setf (aref internals i) (copy-seq x)))))
102 (send obj :add-slot 'internals internals)
103 obj))
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))
115 (when set
116 (setf (aref (aref (slot-value 'internals) 8) 5)
117 (cond ((integerp val) val)
118 ((null val) 0)
119 (t 1))))
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))
132 (when 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))
138 ;;;;
139 ;;;;
140 ;;;; Newton's Method with Backtracking
141 ;;;;
142 ;;;;
144 (defun newtonmax (f start &key
145 scale
146 (derivstep -1.0)
147 (count-limit -1)
148 (verbose 1)
149 return-derivs)
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)
163 (if return-derivs
164 (cons (send minfo :x) (- (send minfo :fvals)))
165 (send minfo :x))))
167 ;;;;
168 ;;;;
169 ;;;; Nelder-Mead Simplex Method
170 ;;;;
171 ;;;;
173 (defun nelmeadmax (f start &key
174 (size 1)
175 (epsilon (sqrt machine-epsilon))
176 (count-limit 500)
177 (verbose t)
178 (alpha 1.0)
179 (beta 0.5)
180 (gamma 2.0)
181 (delta 0.5))
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))
193 next)
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))))
205 (if verbose
206 (format t "Value = ~10g~%"
207 (send s :point-value (send s :best-point)))))))
211 ;;; Simplex Prototype
214 (defproto simplex-proto '(f simplex))
217 ;;; Simplex Points
220 (defmeth simplex-proto :make-point (x)
221 (let ((f (send self :f)))
222 (if f
223 (let ((val (funcall f x)))
224 (cons (if (consp val) (car val) val) x))
225 (cons nil 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)
240 (send self :f f))
243 ;;; Slot Accessors and Mutators
246 (defmeth simplex-proto :simplex (&optional new size)
247 (if new
248 (let ((simplex
249 (if (and (consp new) (sequencep (car new)))
250 (if (/= (length new) (+ 1 (length (car new))))
251 (error "bad simplex data")
252 (copy-list new))
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)
267 (when f
268 (setf (slot-value 'f) f)
269 (let ((simplex
270 (mapcar #'(lambda (x) (send self :point-location x))
271 (send self :simplex))))
272 (send self :simplex simplex)))
273 (slot-value 'f))
275 (defmeth simplex-proto :sort-simplex ()
276 (if (send self :f)
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)))
291 (when n
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))
309 new))
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
318 (+ bloc
319 (* fac
320 (- (send self :point-location x) bloc))))
321 x)))))
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))))))