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 ;;; FIXME:AJR: There is a need to figure out the proper symbols to
  17 ;;; export.  more importantly should there be any specialty packages
  18 ;;; that are exported for maximization?
  19
  20 (in-package :lisp-stat-optimize)
  21
  22 ;;;;
  23 ;;;; minfo basics (internal??)
  24 ;;;;
  25
  26 (defun init-minfo-ipar-values (n ipars)
  27   (let* ((TRUE 1)
  28          (FALSE 0)
  29          (k 0)
  30          (m 0)
  31          (itnlimit -1)
  32          (backtrack TRUE)
  33          (verbose 0)
  34          (vals_suppl FALSE)
  35          (exptilt TRUE)
  36          (count 0)
  37          (termcode 0))
  38     (setf (aref ipars 0) n)
  39     (setf (aref ipars 1) m)
  40     (setf (aref ipars 2) k)
  41     (setf (aref ipars 3) itnlimit)
  42     (setf (aref ipars 4) backtrack)
  43     (setf (aref ipars 5) verbose)
  44     (setf (aref ipars 6) vals_suppl)
  45     (setf (aref ipars 7) exptilt)
  46     (setf (aref ipars 8) count)
  47     (setf (aref ipars 9) termcode)))
  48
  49 (defun init-minfo-dpar-values (h dpars)
  50   (let ((typf 1.0)
  51         (gradtol -1.0)
  52         (steptol -1.0)
  53         (maxstep -1.0)
  54         (dflt 0.0)
  55         (tilt 0.0)
  56         (newtilt 0.0)
  57         (hessadd 0.0))
  58     (setf (aref dpars 0) typf)
  59     (setf (aref dpars 1) h)
  60     (setf (aref dpars 2) gradtol)
  61     (setf (aref dpars 3) steptol)
  62     (setf (aref dpars 4) maxstep)
  63     (setf (aref dpars 5) dflt)
  64     (setf (aref dpars 6) tilt)
  65     (setf (aref dpars 7) newtilt)
  66     (setf (aref dpars 8) hessadd)))
  67
  68 (defun init-minfo-internals (n h internals)
  69   (let ((ipars (aref internals 8))
  70         (dpars (aref internals 9)))
  71     (init-minfo-ipar-values n ipars)
  72     (init-minfo-dpar-values h dpars)))
  73
  74 (defun new-minfo-internals (f x &key scale ((:derivstep h) -1.0))
  75   (check-sequence x)
  76   (check-real x)
  77   (check-one-real h)
  78   (let ((n (length x)))
  79     (when scale
  80           (check-sequence scale)
  81           (check-real scale)
  82           (if (/= n (length scale)) (error "scale and x not the same length")))
  83     (let ((internals (make-array 12)))
  84       (setf (aref internals 0) f)
  85       (setf (aref internals 3) (if (consp x) (copy-list x) (coerce x 'list)))
  86       (setf (aref internals 4)
  87             (if scale (copy-seq scale) (make-array n :initial-element 1.0)))
  88       (setf (aref internals 5) (make-list (+ 1 n (* n n))))
  89       (setf (aref internals 8) (make-array 10))
  90       (setf (aref internals 9) (make-array 9))
  91       (init-minfo-internals n h internals)
  92       internals)))
  93
  94 (defun minfo-maximize (internals &optional verbose)
  95   "This function does what?"
  96   (let* ((f (aref internals 0))
  97          (x (aref internals 3))
  98          (fvals (aref internals 5))
  99          (n (length x))
 100          (v (if verbose (if (integerp verbose) verbose 1) -1)))
 101     (setf (aref internals 3) (copy-list x))
 102     (setf (aref internals 5) (copy-list fvals))
 103     (let ((*maximize-callback-function* f)
 104           (*maximize-callback-arg* (make-list n)))
 105       (let* ((x (aref internals 3))
 106              (scale (aref internals 4))
 107              (fvals (aref internals 5))
 108              (ip (aref internals 8))
 109              (dp (aref internals 9))
 110              (px (la-data-to-vector x mode-re))
 111              (pscale (la-data-to-vector scale mode-re))
 112              (pfvals (la-vector (length fvals) mode-re))
 113              (pip (la-data-to-vector ip mode-in))
 114              (pdp (la-data-to-vector dp mode-re)))
 115         (unwind-protect
 116             (progn
 117               (base-minfo-maximize px pfvals pscale pip pdp v)) ;; access to C
 118           (la-vector-to-data px n mode-re x)
 119           (la-vector-to-data pfvals (+ 1 n (* n n)) mode-re fvals)
 120           (la-vector-to-data pip (length ip) mode-in ip)
 121           (la-vector-to-data pdp (length dp) mode-re dp))
 122         (get-buf)))))
 123
 124
 125
 126 ;;;;
 127 ;;;; Mode Info Prototype
 128 ;;;;
 129
 130 (defproto minfo-proto '(internals))
 131
 132 #+xlisp (send minfo-proto :add-method :isnew #'|minfo-isnew|)
 133 #+xlisp (send minfo-proto :add-method :maximize #'|minfo-maximize|)
 134 #+xlisp (send minfo-proto :add-method :loglaplace #'|minfo-loglap|)
 135 #-xlisp
 136 (defmeth minfo-proto :isnew (&rest args)
 137   (setf (slot-value 'internals) (apply #'new-minfo-internals args)))
 138 #-xlisp
 139 (defmeth minfo-proto :maximize (&rest args)
 140   (apply #'minfo-maximize (slot-value 'internals) args))
 141
 142 (defmeth minfo-proto :x () (aref (slot-value 'internals) 3))
 143 (defmeth minfo-proto :scale () (aref (slot-value 'internals) 4))
 144 (defmeth minfo-proto :derivstep () (aref (aref (slot-value 'internals) 9) 1))
 145 (defmeth minfo-proto :tilt () (aref (aref (slot-value 'internals) 9) 6))
 146
 147 (defmeth minfo-proto :f (&optional (val nil set))
 148   (when set
 149         (send self :set-no-vals-supplied)
 150         (setf (aref (slot-value 'internals) 0) val))
 151   (aref (slot-value 'internals) 0))
 152
 153 (defmeth minfo-proto :set-no-vals-supplied ()
 154   (setf (aref (aref (slot-value 'internals) 8) 6) 0))
 155
 156 (defmeth minfo-proto :exptilt (&optional (val nil set))
 157   (if set
 158       (let ((old (send self :exptilt)))
 159         (setf (aref (aref (slot-value 'internals) 8) 7) (if val 1 0))
 160         (if (and (not (or (and old val) (and (not old) (not val))))
 161                  (/= (send self :tilt) 0.0))
 162             (send self :set-no-vals-supplied))))
 163   (= 1 (aref (aref (slot-value 'internals) 8) 7)))
 164
 165 (defmeth minfo-proto :newtilt (&optional (val nil set))
 166   (when set
 167         (setf (aref (aref (slot-value 'internals) 9) 7) (float val))
 168         (if (/= (send self :tilt) 0.0) (send self :set-no-vals-supplied)))
 169   (aref (aref (slot-value 'internals) 9) 7))
 170
 171 (defmeth minfo-proto :gfuns (&optional (val nil set))
 172   (when set
 173         (if (or (not (consp val))
 174                 (not (every #'functionp val)))
 175             (error "not all functions"))
 176         (setf (aref (slot-value 'internals) 1) val)
 177         (setf (aref (aref (slot-value 'internals) 8) 1) (length val))
 178         (setf (aref (slot-value 'internals) 10) (repeat 1.0 (length val)))
 179         (if (/= (send self :tilt) 0.0) (send self :set-no-vals-supplied)))
 180   (aref (slot-value 'internals) 1))
 181
 182 (defmeth minfo-proto :cfuns (&optional (val nil set))
 183   (when set
 184         (if (or (not (consp val))
 185                 (not (every #'functionp val)))
 186             (error "not all functions"))
 187         (setf (aref (slot-value 'internals) 2) val)
 188         (setf (aref (aref (slot-value 'internals) 8) 2) (length val))
 189         (setf (aref (slot-value 'internals) 7) (repeat 0.0 (length val)))
 190         (setf (aref (slot-value 'internals) 11) (repeat 0.0 (length val)))
 191         (send self :set-no-vals-supplied))
 192   (aref (slot-value 'internals) 2))
 193
 194 (defmeth minfo-proto :ctarget (&optional (val nil set))
 195   (when set
 196         (if (/= (length val) (length (send self :ctarget)))
 197             (error "bad target length"))
 198         (setf (aref (slot-value 'internals) 7) val))
 199   (aref (slot-value 'internals) 7))
 200
 201 (defmeth minfo-proto :fvals ()
 202   (let* ((fv (aref (slot-value 'internals) 5))
 203          (n (length (send self :x)))
 204          (val (select fv 0))
 205          (grad (select fv (iseq 1 n)))
 206          (hess (matrix (list n n) (select fv (iseq (+ 1 n) (+ n (* n n)))))))
 207     (list val grad hess)))
 208
 209 (defmeth minfo-proto :copy ()
 210   (let ((obj (make-object minfo-proto))
 211         (internals (copy-seq (slot-value 'internals))))
 212     (dotimes (i (length internals))
 213              (let ((x (aref internals i)))
 214                (if (sequencep x)
 215                    (setf (aref internals i) (copy-seq x)))))
 216     (send obj :add-slot 'internals internals)
 217     obj))
 218
 219 (defmeth minfo-proto :derivscale ()
 220   (let* ((step (^ machine-epsilon (/ 1 6)))
 221          (hess (numhess (send self :f) (send self :x) (send self :scale) step))
 222          (scale (pmax (abs (send self :x)) (sqrt (abs (/ (diagonal hess)))))))
 223     (setf hess (numhess (send self :f) (send self :x) scale step))
 224     (setf scale (pmax (abs (send self :x)) (sqrt (abs (/ (diagonal hess))))))
 225     (setf (aref (slot-value 'internals) 4) scale)
 226     (setf (aref (aref (slot-value 'internals) 9) 1) step)))
 227
 228 (defmeth minfo-proto :verbose (&optional (val nil set))
 229   (when set
 230         (setf (aref (aref (slot-value 'internals) 8) 5)
 231               (cond ((integerp val) val)
 232                     ((null val) 0)
 233                     (t 1))))
 234   (aref (aref (slot-value 'internals) 8) 5))
 235
 236 (defmeth minfo-proto :backtrack (&optional (val nil set))
 237   (if set (setf (aref (aref (slot-value 'internals) 8) 4) (if val 1 0)))
 238   (aref (aref (slot-value 'internals) 8) 4))
 239
 240 (defmeth minfo-proto :maxiter (&optional (val nil set))
 241     (if set (setf (aref (aref (slot-value 'internals) 8) 3)
 242                   (if (integerp val) val -1)))
 243     (aref (aref (slot-value 'internals) 8) 3))
 244
 245 (defmeth minfo-proto :tiltscale (&optional (val nil set))
 246   (when set
 247         (if (/= (length val) (length (send self :gfuns)))
 248             (error "wrong size tilt scale sequence"))
 249         (setf (aref (slot-value 'internals) 10) val))
 250   (aref (slot-value 'internals) 10))
 251
 252 ;;;;
 253 ;;;;
 254 ;;;; Newton's Method with Backtracking
 255 ;;;;
 256 ;;;;
 257
 258 (defun newtonmax (f start &key
 259                     scale
 260                     (derivstep -1.0)
 261                     (count-limit -1)
 262                     (verbose 1)
 263                     return-derivs)
 264 "Args:(f start &key scale derivstep (verbose 1) return-derivs)
 265 Maximizes F starting from START using Newton's method with backtracking.
 266 If RETURN-DERIVS is NIL returns location of maximum; otherwise returns
 267 list of location, unction value, gradient and hessian at maximum.
 268 SCALE should be a list of the typical magnitudes of the parameters.
 269 DERIVSTEP is used in numerical derivatives and VERBOSE controls printing
 270 of iteration information. COUNT-LIMIT limits the number of iterations"
 271   (let ((verbose (if verbose (if (integerp verbose) verbose 1) 0))
 272         (minfo (send minfo-proto :new f start
 273                      :scale scale :derivstep derivstep)))
 274     (send minfo :maxiter count-limit)
 275     (send minfo :derivscale)
 276     (send minfo :maximize verbose)
 277     (if return-derivs
 278         (cons (send minfo :x) (- (send minfo :fvals)))
 279         (send minfo :x))))
 280
 281 ;;;;
 282 ;;;;
 283 ;;;; Nelder-Mead Simplex Method
 284 ;;;;
 285 ;;;;
 286
 287 (defun nelmeadmax (f start &key
 288                      (size 1)
 289                      (epsilon (sqrt machine-epsilon))
 290                      (count-limit 500)
 291                      (verbose t)
 292                      (alpha 1.0)
 293                      (beta 0.5)
 294                      (gamma 2.0)
 295                      (delta 0.5))
 296 "Args: (f start &key (size 1) (epsilon (sqrt machine-epsilon))
 297           (count-limit 500) (verbose t) alpha beta gamma delta)
 298 Maximizes F using the Nelder-Mead simplex method. START can be a
 299 starting simplex - a list of N+1 points, with N=dimension of problem,
 300 or a single point. If start is a single point you should give the
 301 size of the initial simplex as SIZE, a sequence of length N. Default is
 302 all 1's. EPSILON is the convergence tolerance. ALPHA-DELTA can be used to
 303 control the behavior of simplex algorithm."
 304     (let ((s (send simplex-proto :new f start size)))
 305       (do ((best (send s :best-point) (send s :best-point))
 306            (count 0 (+ count 1))
 307            next)
 308           ((or (< (send s :relative-range) epsilon) (>= count count-limit))
 309            (if (and verbose (>= count count-limit))
 310                (format t "Iteration limit exceeded.~%"))
 311            (send s :point-location (send s :best-point)))
 312           (setf next (send s :extrapolate-from-worst (- alpha)))
 313           (if (send s :is-worse best next)
 314               (setf next (send s :extrapolate-from-worst gamma))
 315               (when (send s :is-worse next (send s :second-worst-point))
 316                     (setf next (send s :extrapolate-from-worst beta))
 317                     (if (send s :is-worse next (send s :worst-point))
 318                         (send s :shrink-to-best delta))))
 319           (if verbose
 320               (format t "Value = ~10g~%"
 321                       (send s :point-value (send s :best-point)))))))
 322
 323
 324 ;;;
 325 ;;; Simplex Prototype
 326 ;;;
 327
 328 (defproto simplex-proto '(f simplex))
 329
 330 ;;;
 331 ;;; Simplex Points
 332 ;;;
 333
 334 (defmeth simplex-proto :make-point (x)
 335   (let ((f (send self :f)))
 336     (if f
 337         (let ((val (funcall f x)))
 338           (cons (if (consp val) (car val) val) x))
 339         (cons nil x))))
 340
 341 (defmeth simplex-proto :point-value (x) (car x))
 342
 343 (defmeth simplex-proto :point-location (x) (cdr x))
 344
 345 (defmeth simplex-proto :is-worse (x y)
 346   (< (send self :point-value x) (send self :point-value y)))
 347
 348 ;;;
 349 ;;; Making New Simplices
 350 ;;;
 351
 352 (defmeth simplex-proto :isnew (f start &optional size)
 353   (send self :simplex start size)
 354   (send self :f f))
 355
 356 ;;;
 357 ;;; Slot Accessors and Mutators
 358 ;;;
 359
 360 (defmeth simplex-proto :simplex (&optional new size)
 361   (if new
 362       (let ((simplex
 363              (if (and (consp new) (sequencep (car new)))
 364                  (if (/= (length new) (+ 1 (length (car new))))
 365                      (error "bad simplex data")
 366                      (copy-list new))
 367                  (let* ((n (length new))
 368                         (size (if size size (repeat 1 n)))
 369                      ;   (pts (- (* 2 (uniform-rand (repeat n (+ n 1)))) 1)))
 370                         (diag (* 2 size (- (random (repeat 2 n)) .5)))
 371                         (pts (cons (repeat 0 n)
 372                                    (mapcar #'(lambda (x) (coerce x 'list))
 373                                            (column-list (diagonal diag))))))
 374                    (mapcar #'(lambda (x) (+ (* size x) new)) pts)))))
 375         (setf (slot-value 'simplex)
 376               (mapcar #'(lambda (x) (send self :make-point x)) simplex))
 377         (send self :sort-simplex)))
 378   (slot-value 'simplex))
 379
 380 (defmeth simplex-proto :f (&optional f)
 381   (when f
 382         (setf (slot-value 'f) f)
 383         (let ((simplex
 384                (mapcar #'(lambda (x) (send self :point-location x))
 385                        (send self :simplex))))
 386           (send self :simplex simplex)))
 387   (slot-value 'f))
 388
 389 (defmeth simplex-proto :sort-simplex ()
 390   (if (send self :f)
 391       (setf (slot-value 'simplex)
 392             (sort (slot-value 'simplex)
 393                   #'(lambda (x y) (send self :is-worse x y))))))
 394
 395 ;;;
 396 ;;; Other Methods Using List Representation of SImplex
 397 ;;;
 398
 399 (defmeth simplex-proto :best-point () (car (last (send self :simplex))))
 400 (defmeth simplex-proto :worst-point () (first (send self :simplex)))
 401 (defmeth simplex-proto :second-worst-point () (second (send self :simplex)))
 402 (defmeth simplex-proto :replace-point (new old)
 403   (let* ((simplex (send self :simplex))
 404          (n (position old simplex)))
 405     (when n
 406           (setf (nth n simplex) new)
 407           (send self :sort-simplex))))
 408 (defmeth simplex-proto :mean-opposite-face (x)
 409   (let ((face (mapcar #'(lambda (x) (send self :point-location x))
 410                       (remove x (send self :simplex)))))
 411     (/ (reduce #'+ face) (length face))))
 412
 413 ;;;
 414 ;;; Iteration Step Methods
 415 ;;;
 416
 417 (defmeth simplex-proto :extrapolate-from-worst (fac)
 418   (let* ((worst (send self :worst-point))
 419          (wloc (send self :point-location worst))
 420          (delta (- (send self :mean-opposite-face worst) wloc))
 421          (new (send self :make-point (+ wloc (* (- 1 fac) delta)))))
 422     (if (send self :is-worse worst new) (send self :replace-point new worst))
 423     new))
 424
 425 (defmeth simplex-proto :shrink-to-best (fac)
 426   (let* ((best (send self :best-point))
 427          (bloc (send self :point-location best)))
 428     (dolist (x (copy-list (send self :simplex)))
 429             (if (not (eq x best))
 430                 (send self :replace-point
 431                       (send self :make-point
 432                             (+ bloc
 433                                (* fac
 434                                   (- (send self :point-location x) bloc))))
 435                       x)))))
 436
 437 (defmeth simplex-proto :relative-range ()
 438   (let ((best (send self :point-value (send self :best-point)))
 439         (worst (send self :point-value (send self :worst-point))))
 440     (* 2 (/ (abs (- best worst)) (+ 1 (abs best) (abs worst))))))
 441
 442
 443