maximize.lsp

   1 (provide "maximize")
   2
   3 #+:CLtL2
   4 (in-package lisp-stat)
   5 #-:CLtL2
   6 (in-package 'lisp-stat)
   7
   8 (export '(newtonmax nelmeadmax))
   9
  10 (import '(ls-basics::new-minfo-internals ls-basics::minfo-maximize))
  11
  12 ;;;;
  13 ;;;; Mode Info Prototype
  14 ;;;;
  15
  16 (defproto minfo-proto '(internals))
  17
  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))
  27
  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))
  32
  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))
  38
  39 (defmeth minfo-proto :set-no-vals-supplied ()
  40   (setf (aref (aref (slot-value 'internals) 8) 6) 0))
  41
  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)))
  50
  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))
  56
  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))
  67
  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))
  79
  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))
  86
  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)))
  94
  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))
 104
 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)))
 113
 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))
 121
 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))
 125
 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))
 130
 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))
 137
 138 ;;;;
 139 ;;;;
 140 ;;;; Newton's Method with Backtracking
 141 ;;;;
 142 ;;;;
 143
 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))))
 166
 167 ;;;;
 168 ;;;;
 169 ;;;; Nelder-Mead Simplex Method
 170 ;;;;
 171 ;;;;
 172
 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)))))))
 208
 209
 210 ;;;
 211 ;;; Simplex Prototype
 212 ;;;
 213
 214 (defproto simplex-proto '(f simplex))
 215
 216 ;;;
 217 ;;; Simplex Points
 218 ;;;
 219
 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))))
 226
 227 (defmeth simplex-proto :point-value (x) (car x))
 228
 229 (defmeth simplex-proto :point-location (x) (cdr x))
 230
 231 (defmeth simplex-proto :is-worse (x y)
 232   (< (send self :point-value x) (send self :point-value y)))
 233
 234 ;;;
 235 ;;; Making New Simplices
 236 ;;;
 237
 238 (defmeth simplex-proto :isnew (f start &optional size)
 239   (send self :simplex start size)
 240   (send self :f f))
 241
 242 ;;;
 243 ;;; Slot Accessors and Mutators
 244 ;;;
 245
 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))
 265
 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))
 274
 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))))))
 280
 281 ;;;
 282 ;;; Other Methods Using List Representation of SImplex
 283 ;;;
 284
 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))))
 298
 299 ;;;
 300 ;;; Iteration Step Methods
 301 ;;;
 302
 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))
 310
 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)))))
 322
 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))))))