optimize.lisp

   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 (in-package :cl-user)
   7
   8 (defpackage :lisp-stat-optimize
   9  (:use :common-lisp
  10        :lisp-stat-object-system
  11        :lisp-stat-basics)
  12  (:shadowing-import-from :lisp-stat-object-system
  13                          slot-value call-method call-next-method)
  14  (:export newtonmax nelmeadmax))
  15
  16 (in-package :lisp-stat-optimize)
  17
  18 ;;;;
  19 ;;;; Mode Info Prototype
  20 ;;;;
  21
  22 (defproto minfo-proto '(internals))
  23
  24 #+xlisp (send minfo-proto :add-method :isnew #'|minfo-isnew|)
  25 #+xlisp (send minfo-proto :add-method :maximize #'|minfo-maximize|)
  26 #+xlisp (send minfo-proto :add-method :loglaplace #'|minfo-loglap|)
  27 #-xlisp
  28 (defmeth minfo-proto :isnew (&rest args)
  29   (setf (slot-value 'internals) (apply #'new-minfo-internals args)))
  30 #-xlisp
  31 (defmeth minfo-proto :maximize (&rest args)
  32   (apply #'minfo-maximize (slot-value 'internals) args))
  33
  34 (defmeth minfo-proto :x () (aref (slot-value 'internals) 3))
  35 (defmeth minfo-proto :scale () (aref (slot-value 'internals) 4))
  36 (defmeth minfo-proto :derivstep () (aref (aref (slot-value 'internals) 9) 1))
  37 (defmeth minfo-proto :tilt () (aref (aref (slot-value 'internals) 9) 6))
  38
  39 (defmeth minfo-proto :f (&optional (val nil set))
  40   (when set
  41         (send self :set-no-vals-supplied)
  42         (setf (aref (slot-value 'internals) 0) val))
  43   (aref (slot-value 'internals) 0))
  44
  45 (defmeth minfo-proto :set-no-vals-supplied ()
  46   (setf (aref (aref (slot-value 'internals) 8) 6) 0))
  47
  48 (defmeth minfo-proto :exptilt (&optional (val nil set))
  49   (if set
  50       (let ((old (send self :exptilt)))
  51         (setf (aref (aref (slot-value 'internals) 8) 7) (if val 1 0))
  52         (if (and (not (or (and old val) (and (not old) (not val))))
  53                  (/= (send self :tilt) 0.0))
  54             (send self :set-no-vals-supplied))))
  55   (= 1 (aref (aref (slot-value 'internals) 8) 7)))
  56
  57 (defmeth minfo-proto :newtilt (&optional (val nil set))
  58   (when set
  59         (setf (aref (aref (slot-value 'internals) 9) 7) (float val))
  60         (if (/= (send self :tilt) 0.0) (send self :set-no-vals-supplied)))
  61   (aref (aref (slot-value 'internals) 9) 7))
  62
  63 (defmeth minfo-proto :gfuns (&optional (val nil set))
  64   (when set
  65         (if (or (not (consp val))
  66                 (not (every #'functionp val)))
  67             (error "not all functions"))
  68         (setf (aref (slot-value 'internals) 1) val)
  69         (setf (aref (aref (slot-value 'internals) 8) 1) (length val))
  70         (setf (aref (slot-value 'internals) 10) (repeat 1.0 (length val)))
  71         (if (/= (send self :tilt) 0.0) (send self :set-no-vals-supplied)))
  72   (aref (slot-value 'internals) 1))
  73
  74 (defmeth minfo-proto :cfuns (&optional (val nil set))
  75   (when set
  76         (if (or (not (consp val))
  77                 (not (every #'functionp val)))
  78             (error "not all functions"))
  79         (setf (aref (slot-value 'internals) 2) val)
  80         (setf (aref (aref (slot-value 'internals) 8) 2) (length val))
  81         (setf (aref (slot-value 'internals) 7) (repeat 0.0 (length val)))
  82         (setf (aref (slot-value 'internals) 11) (repeat 0.0 (length val)))
  83         (send self :set-no-vals-supplied))
  84   (aref (slot-value 'internals) 2))
  85
  86 (defmeth minfo-proto :ctarget (&optional (val nil set))
  87   (when set
  88         (if (/= (length val) (length (send self :ctarget)))
  89             (error "bad target length"))
  90         (setf (aref (slot-value 'internals) 7) val))
  91   (aref (slot-value 'internals) 7))
  92
  93 (defmeth minfo-proto :fvals ()
  94   (let* ((fv (aref (slot-value 'internals) 5))
  95          (n (length (send self :x)))
  96          (val (select fv 0))
  97          (grad (select fv (iseq 1 n)))
  98          (hess (matrix (list n n) (select fv (iseq (+ 1 n) (+ n (* n n)))))))
  99     (list val grad hess)))
 100
 101 (defmeth minfo-proto :copy ()
 102   (let ((obj (make-object minfo-proto))
 103         (internals (copy-seq (slot-value 'internals))))
 104     (dotimes (i (length internals))
 105              (let ((x (aref internals i)))
 106                (if (sequencep x)
 107                    (setf (aref internals i) (copy-seq x)))))
 108     (send obj :add-slot 'internals internals)
 109     obj))
 110
 111 (defmeth minfo-proto :derivscale ()
 112   (let* ((step (^ machine-epsilon (/ 1 6)))
 113          (hess (numhess (send self :f) (send self :x) (send self :scale) step))
 114          (scale (pmax (abs (send self :x)) (sqrt (abs (/ (diagonal hess)))))))
 115     (setf hess (numhess (send self :f) (send self :x) scale step))
 116     (setf scale (pmax (abs (send self :x)) (sqrt (abs (/ (diagonal hess))))))
 117     (setf (aref (slot-value 'internals) 4) scale)
 118     (setf (aref (aref (slot-value 'internals) 9) 1) step)))
 119
 120 (defmeth minfo-proto :verbose (&optional (val nil set))
 121   (when set
 122         (setf (aref (aref (slot-value 'internals) 8) 5)
 123               (cond ((integerp val) val)
 124                     ((null val) 0)
 125                     (t 1))))
 126   (aref (aref (slot-value 'internals) 8) 5))
 127
 128 (defmeth minfo-proto :backtrack (&optional (val nil set))
 129   (if set (setf (aref (aref (slot-value 'internals) 8) 4) (if val 1 0)))
 130   (aref (aref (slot-value 'internals) 8) 4))
 131
 132 (defmeth minfo-proto :maxiter (&optional (val nil set))
 133     (if set (setf (aref (aref (slot-value 'internals) 8) 3)
 134                   (if (integerp val) val -1)))
 135     (aref (aref (slot-value 'internals) 8) 3))
 136
 137 (defmeth minfo-proto :tiltscale (&optional (val nil set))
 138   (when set
 139         (if (/= (length val) (length (send self :gfuns)))
 140             (error "wrong size tilt scale sequence"))
 141         (setf (aref (slot-value 'internals) 10) val))
 142   (aref (slot-value 'internals) 10))
 143
 144 ;;;;
 145 ;;;;
 146 ;;;; Newton's Method with Backtracking
 147 ;;;;
 148 ;;;;
 149
 150 (defun newtonmax (f start &key
 151                     scale
 152                     (derivstep -1.0)
 153                     (count-limit -1)
 154                     (verbose 1)
 155                     return-derivs)
 156 "Args:(f start &key scale derivstep (verbose 1) return-derivs)
 157 Maximizes F starting from START using Newton's method with backtracking.
 158 If RETURN-DERIVS is NIL returns location of maximum; otherwise returns
 159 list of location, unction value, gradient and hessian at maximum.
 160 SCALE should be a list of the typical magnitudes of the parameters.
 161 DERIVSTEP is used in numerical derivatives and VERBOSE controls printing
 162 of iteration information. COUNT-LIMIT limits the number of iterations"
 163   (let ((verbose (if verbose (if (integerp verbose) verbose 1) 0))
 164         (minfo (send minfo-proto :new f start
 165                      :scale scale :derivstep derivstep)))
 166     (send minfo :maxiter count-limit)
 167     (send minfo :derivscale)
 168     (send minfo :maximize verbose)
 169     (if return-derivs
 170         (cons (send minfo :x) (- (send minfo :fvals)))
 171         (send minfo :x))))
 172
 173 ;;;;
 174 ;;;;
 175 ;;;; Nelder-Mead Simplex Method
 176 ;;;;
 177 ;;;;
 178
 179 (defun nelmeadmax (f start &key
 180                      (size 1)
 181                      (epsilon (sqrt machine-epsilon))
 182                      (count-limit 500)
 183                      (verbose t)
 184                      (alpha 1.0)
 185                      (beta 0.5)
 186                      (gamma 2.0)
 187                      (delta 0.5))
 188 "Args: (f start &key (size 1) (epsilon (sqrt machine-epsilon))
 189           (count-limit 500) (verbose t) alpha beta gamma delta)
 190 Maximizes F using the Nelder-Mead simplex method. START can be a
 191 starting simplex - a list of N+1 points, with N=dimension of problem,
 192 or a single point. If start is a single point you should give the
 193 size of the initial simplex as SIZE, a sequence of length N. Default is
 194 all 1's. EPSILON is the convergence tolerance. ALPHA-DELTA can be used to
 195 control the behavior of simplex algorithm."
 196     (let ((s (send simplex-proto :new f start size)))
 197       (do ((best (send s :best-point) (send s :best-point))
 198            (count 0 (+ count 1))
 199            next)
 200           ((or (< (send s :relative-range) epsilon) (>= count count-limit))
 201            (if (and verbose (>= count count-limit))
 202                (format t "Iteration limit exceeded.~%"))
 203            (send s :point-location (send s :best-point)))
 204           (setf next (send s :extrapolate-from-worst (- alpha)))
 205           (if (send s :is-worse best next)
 206               (setf next (send s :extrapolate-from-worst gamma))
 207               (when (send s :is-worse next (send s :second-worst-point))
 208                     (setf next (send s :extrapolate-from-worst beta))
 209                     (if (send s :is-worse next (send s :worst-point))
 210                         (send s :shrink-to-best delta))))
 211           (if verbose
 212               (format t "Value = ~10g~%"
 213                       (send s :point-value (send s :best-point)))))))
 214
 215
 216 ;;;
 217 ;;; Simplex Prototype
 218 ;;;
 219
 220 (defproto simplex-proto '(f simplex))
 221
 222 ;;;
 223 ;;; Simplex Points
 224 ;;;
 225
 226 (defmeth simplex-proto :make-point (x)
 227   (let ((f (send self :f)))
 228     (if f
 229         (let ((val (funcall f x)))
 230           (cons (if (consp val) (car val) val) x))
 231         (cons nil x))))
 232
 233 (defmeth simplex-proto :point-value (x) (car x))
 234
 235 (defmeth simplex-proto :point-location (x) (cdr x))
 236
 237 (defmeth simplex-proto :is-worse (x y)
 238   (< (send self :point-value x) (send self :point-value y)))
 239
 240 ;;;
 241 ;;; Making New Simplices
 242 ;;;
 243
 244 (defmeth simplex-proto :isnew (f start &optional size)
 245   (send self :simplex start size)
 246   (send self :f f))
 247
 248 ;;;
 249 ;;; Slot Accessors and Mutators
 250 ;;;
 251
 252 (defmeth simplex-proto :simplex (&optional new size)
 253   (if new
 254       (let ((simplex
 255              (if (and (consp new) (sequencep (car new)))
 256                  (if (/= (length new) (+ 1 (length (car new))))
 257                      (error "bad simplex data")
 258                      (copy-list new))
 259                  (let* ((n (length new))
 260                         (size (if size size (repeat 1 n)))
 261                      ;   (pts (- (* 2 (uniform-rand (repeat n (+ n 1)))) 1)))
 262                         (diag (* 2 size (- (random (repeat 2 n)) .5)))
 263                         (pts (cons (repeat 0 n)
 264                                    (mapcar #'(lambda (x) (coerce x 'list))
 265                                            (column-list (diagonal diag))))))
 266                    (mapcar #'(lambda (x) (+ (* size x) new)) pts)))))
 267         (setf (slot-value 'simplex)
 268               (mapcar #'(lambda (x) (send self :make-point x)) simplex))
 269         (send self :sort-simplex)))
 270   (slot-value 'simplex))
 271
 272 (defmeth simplex-proto :f (&optional f)
 273   (when f
 274         (setf (slot-value 'f) f)
 275         (let ((simplex
 276                (mapcar #'(lambda (x) (send self :point-location x))
 277                        (send self :simplex))))
 278           (send self :simplex simplex)))
 279   (slot-value 'f))
 280
 281 (defmeth simplex-proto :sort-simplex ()
 282   (if (send self :f)
 283       (setf (slot-value 'simplex)
 284             (sort (slot-value 'simplex)
 285                   #'(lambda (x y) (send self :is-worse x y))))))
 286
 287 ;;;
 288 ;;; Other Methods Using List Representation of SImplex
 289 ;;;
 290
 291 (defmeth simplex-proto :best-point () (car (last (send self :simplex))))
 292 (defmeth simplex-proto :worst-point () (first (send self :simplex)))
 293 (defmeth simplex-proto :second-worst-point () (second (send self :simplex)))
 294 (defmeth simplex-proto :replace-point (new old)
 295   (let* ((simplex (send self :simplex))
 296          (n (position old simplex)))
 297     (when n
 298           (setf (nth n simplex) new)
 299           (send self :sort-simplex))))
 300 (defmeth simplex-proto :mean-opposite-face (x)
 301   (let ((face (mapcar #'(lambda (x) (send self :point-location x))
 302                       (remove x (send self :simplex)))))
 303     (/ (apply #'+ face) (length face))))
 304
 305 ;;;
 306 ;;; Iteration Step Methods
 307 ;;;
 308
 309 (defmeth simplex-proto :extrapolate-from-worst (fac)
 310   (let* ((worst (send self :worst-point))
 311          (wloc (send self :point-location worst))
 312          (delta (- (send self :mean-opposite-face worst) wloc))
 313          (new (send self :make-point (+ wloc (* (- 1 fac) delta)))))
 314     (if (send self :is-worse worst new) (send self :replace-point new worst))
 315     new))
 316
 317 (defmeth simplex-proto :shrink-to-best (fac)
 318   (let* ((best (send self :best-point))
 319          (bloc (send self :point-location best)))
 320     (dolist (x (copy-list (send self :simplex)))
 321             (if (not (eq x best))
 322                 (send self :replace-point
 323                       (send self :make-point
 324                             (+ bloc
 325                                (* fac
 326                                   (- (send self :point-location x) bloc))))
 327                       x)))))
 328
 329 (defmeth simplex-proto :relative-range ()
 330   (let ((best (send self :point-value (send self :best-point)))
 331         (worst (send self :point-value (send self :worst-point))))
 332     (* 2 (/ (abs (- best worst)) (+ 1 (abs best) (abs worst))))))
 333
 334
 335