maximize.lsp

   1 ;;; -*- mode: lisp -*-
   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.
   5
   6 (defpackage :lisp-stat-optimize
   7  (:use :common-lisp
   8        :lisp-stat-object-system
   9        :lisp-stat-basics)
  10  (:shadowing-import-from :lisp-stat-object-system
  11                          slot-value call-method call-next-method)
  12  (:export newtonmax nelmeadmax))
  13
  14
  15 ;;;;
  16 ;;;; Mode Info Prototype
  17 ;;;;
  18
  19 (defproto minfo-proto '(internals))
  20
  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|)
  24 #-xlisp
  25 (defmeth minfo-proto :isnew (&rest args)
  26   (setf (slot-value 'internals) (apply #'new-minfo-internals args)))
  27 #-xlisp
  28 (defmeth minfo-proto :maximize (&rest args)
  29   (apply #'minfo-maximize (slot-value 'internals) args))
  30
  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))
  35
  36 (defmeth minfo-proto :f (&optional (val nil set))
  37   (when set
  38         (send self :set-no-vals-supplied)
  39         (setf (aref (slot-value 'internals) 0) val))
  40   (aref (slot-value 'internals) 0))
  41
  42 (defmeth minfo-proto :set-no-vals-supplied ()
  43   (setf (aref (aref (slot-value 'internals) 8) 6) 0))
  44
  45 (defmeth minfo-proto :exptilt (&optional (val nil set))
  46   (if 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)))
  53
  54 (defmeth minfo-proto :newtilt (&optional (val nil set))
  55   (when 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))
  59
  60 (defmeth minfo-proto :gfuns (&optional (val nil set))
  61   (when 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))
  70
  71 (defmeth minfo-proto :cfuns (&optional (val nil set))
  72   (when 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))
  82
  83 (defmeth minfo-proto :ctarget (&optional (val nil set))
  84   (when 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))
  89
  90 (defmeth minfo-proto :fvals ()
  91   (let* ((fv (aref (slot-value 'internals) 5))
  92          (n (length (send self :x)))
  93          (val (select fv 0))
  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)))
  97
  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)))
 103                (if (sequencep x)
 104                    (setf (aref internals i) (copy-seq x)))))
 105     (send obj :add-slot 'internals internals)
 106     obj))
 107
 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)))
 116
 117 (defmeth minfo-proto :verbose (&optional (val nil set))
 118   (when set
 119         (setf (aref (aref (slot-value 'internals) 8) 5)
 120               (cond ((integerp val) val)
 121                     ((null val) 0)
 122                     (t 1))))
 123   (aref (aref (slot-value 'internals) 8) 5))
 124
 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))
 128
 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))
 133
 134 (defmeth minfo-proto :tiltscale (&optional (val nil set))
 135   (when 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))
 140
 141 ;;;;
 142 ;;;;
 143 ;;;; Newton's Method with Backtracking
 144 ;;;;
 145 ;;;;
 146
 147 (defun newtonmax (f start &key
 148                     scale
 149                     (derivstep -1.0)
 150                     (count-limit -1)
 151                     (verbose 1)
 152                     return-derivs)
 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)
 166     (if return-derivs
 167         (cons (send minfo :x) (- (send minfo :fvals)))
 168         (send minfo :x))))
 169
 170 ;;;;
 171 ;;;;
 172 ;;;; Nelder-Mead Simplex Method
 173 ;;;;
 174 ;;;;
 175
 176 (defun nelmeadmax (f start &key
 177                      (size 1)
 178                      (epsilon (sqrt machine-epsilon))
 179                      (count-limit 500)
 180                      (verbose t)
 181                      (alpha 1.0)
 182                      (beta 0.5)
 183                      (gamma 2.0)
 184                      (delta 0.5))
 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))
 196            next)
 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))))
 208           (if verbose
 209               (format t "Value = ~10g~%"
 210                       (send s :point-value (send s :best-point)))))))
 211
 212
 213 ;;;
 214 ;;; Simplex Prototype
 215 ;;;
 216
 217 (defproto simplex-proto '(f simplex))
 218
 219 ;;;
 220 ;;; Simplex Points
 221 ;;;
 222
 223 (defmeth simplex-proto :make-point (x)
 224   (let ((f (send self :f)))
 225     (if f
 226         (let ((val (funcall f x)))
 227           (cons (if (consp val) (car val) val) x))
 228         (cons nil x))))
 229
 230 (defmeth simplex-proto :point-value (x) (car x))
 231
 232 (defmeth simplex-proto :point-location (x) (cdr x))
 233
 234 (defmeth simplex-proto :is-worse (x y)
 235   (< (send self :point-value x) (send self :point-value y)))
 236
 237 ;;;
 238 ;;; Making New Simplices
 239 ;;;
 240
 241 (defmeth simplex-proto :isnew (f start &optional size)
 242   (send self :simplex start size)
 243   (send self :f f))
 244
 245 ;;;
 246 ;;; Slot Accessors and Mutators
 247 ;;;
 248
 249 (defmeth simplex-proto :simplex (&optional new size)
 250   (if new
 251       (let ((simplex
 252              (if (and (consp new) (sequencep (car new)))
 253                  (if (/= (length new) (+ 1 (length (car new))))
 254                      (error "bad simplex data")
 255                      (copy-list new))
 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))
 268
 269 (defmeth simplex-proto :f (&optional f)
 270   (when f
 271         (setf (slot-value 'f) f)
 272         (let ((simplex
 273                (mapcar #'(lambda (x) (send self :point-location x))
 274                        (send self :simplex))))
 275           (send self :simplex simplex)))
 276   (slot-value 'f))
 277
 278 (defmeth simplex-proto :sort-simplex ()
 279   (if (send self :f)
 280       (setf (slot-value 'simplex)
 281             (sort (slot-value 'simplex)
 282                   #'(lambda (x y) (send self :is-worse x y))))))
 283
 284 ;;;
 285 ;;; Other Methods Using List Representation of SImplex
 286 ;;;
 287
 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)))
 294     (when n
 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))))
 301
 302 ;;;
 303 ;;; Iteration Step Methods
 304 ;;;
 305
 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))
 312     new))
 313
 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
 321                             (+ bloc
 322                                (* fac
 323                                   (- (send self :point-location x) bloc))))
 324                       x)))))
 325
 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))))))
 330
 331
 332