src/lsp/maximize.lsp

   1 ;;;; maxmize.lsp -- maximization and derivative code for XLISP-STAT
   2 ;;;; XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney
   3 ;;;; Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz
   4 ;;;; You may give out copies of this software; for conditions see the
   5 ;;;; file COPYING included with this distribution.
   6
   7 (provide "maximize")
   8
   9 (in-package "XLISP")
  10
  11 (export '(numgrad numhess newtonmax nelmeadmax))
  12
  13 ;;**** allow typf to be specified?
  14
  15 (defparameter *default-derivative-step* (^ machine-epsilon (/ 1.0 6.0)))
  16
  17 (defmacro make-fundata () `(make-array 7))
  18 (defmacro fundata-f (fd) `(aref ,fd 0))
  19 (defmacro fundata-arg (fd) `(aref ,fd 1))
  20 (defmacro fundata-typx (fd) `(aref ,fd 2))
  21 (defmacro fundata-fsum (fd) `(aref ,fd 3))
  22 (defmacro fundata-n (fd) `(aref ,fd 4))
  23 (defmacro fundata-changesign (fd) `(aref ,fd 5))
  24 (defmacro fundata-h (fd) `(aref ,fd 6))
  25
  26 (defmacro get-function-value (f x arg cs)
  27   (let ((xsym (gensym "X"))
  28         (asym (gensym "ARG")))
  29     `(let ((,xsym ,x)
  30            (,asym ,arg))
  31        (replace ,asym ,xsym)
  32        (let* ((v (funcall ,f ,asym))
  33               (fv (if (consp v) (first v) v)))
  34          (if ,cs (- fv) fv)))))
  35
  36 (defun change-array-sign (x)
  37   (let ((x (compound-data-seq x)))
  38     (if (typep x '(vector float))
  39         (blas-dscal (length x) -1.0 x 0 1)
  40         (dotimes (i (length x)) (setf (aref x i) (- (aref x i)))))))
  41
  42
  43 ;;;;
  44 ;;;; Functions Evaluation Routines
  45 ;;;;
  46
  47 ;;;
  48 ;;; All Hessian evaluations by numerical derivatives assume the gradient is
  49 ;;; evaluated first at the same location. The results are cached away.
  50 ;;;
  51
  52 (defun make-function-data (f n changesign h typx)
  53   (let ((data (make-fundata)))
  54     (setf (fundata-f data) f)
  55     (setf (fundata-arg data) (make-list n))
  56     (setf (fundata-typx data)
  57           (make-array n :element-type 'float :initial-element 1.0))
  58     (setf (fundata-fsum data) (make-array n :element-type 'float))
  59     (setf (fundata-n data) n)
  60     (setf (fundata-changesign data) changesign)
  61     (setf (fundata-h data) (if (and h (> h 0.0)) h *default-derivative-step*))
  62     (when typx
  63           (let ((newtypx (fundata-typx data)))
  64             (dotimes (i n)
  65               (let ((txi (elt typx i)))
  66                 (setf (aref newtypx i) (if (> txi 0.0) txi 1.0))))))
  67     data))
  68
  69 (defun eval-function (func x grad hess)
  70   (let* ((f (fundata-f func))
  71          (arg (fundata-arg func))
  72          (typx (fundata-typx func))
  73          (fsum (fundata-fsum func))
  74          (n (fundata-n func))
  75          (cs (fundata-changesign func))
  76          (h (fundata-h func)))
  77     (replace arg x)
  78     (let* ((result (funcall f arg))
  79            (val (if (consp result) (car result) result)))
  80       (when result
  81             (let ((fderivs (if (consp result) (- (length result) 1) 0)))
  82               (when cs (setf val (- val)))
  83               (when grad
  84                     (cond
  85                      ((> fderivs 0)
  86                       (replace grad (car (cdr result)))
  87                       (when cs (change-array-sign grad)))
  88                      (t
  89                       (cmpgrad n x grad fsum func h typx))))
  90               (when hess
  91                     (cond
  92                      ((= fderivs 2)
  93                       (replace (compound-data-seq hess)
  94                                (compound-data-seq (car (cdr (cdr result)))))
  95                       (when cs (change-array-sign hess)))
  96                      (t
  97                       ;;**** fix this for hess from grad?
  98                       (when (and (= fderivs 1) (null grad))
  99                             ;; kludge to get fsum for analytic gradients
 100                             (cmpgrad n x fsum fsum func h typx))
 101                       (cmphess n x hess val fsum func h typx))))))
 102       val)))
 103
 104 (defun cmpgrad (n x grad fsum ffun h typx)
 105   (let ((f (fundata-f ffun))
 106         (arg (fundata-arg ffun))
 107         (cs (fundata-changesign ffun)))
 108     (dotimes (i n)
 109       (let ((old-xi (aref x i))
 110             (hi (if typx (* (aref typx i) h) h))
 111             (f1 0.0)
 112             (f2 0.0))
 113         (setf (aref x i) (+ old-xi hi))
 114         (setf f1 (get-function-value f x arg cs))
 115         (setf (aref x i) (- old-xi hi))
 116         (setf f2 (get-function-value f x arg cs))
 117         (setf (aref x i) old-xi)
 118         (setf (aref grad i) (/ (- f1 f2) (* 2.0 hi)))
 119         (setf (aref fsum i) (+ f1 f2))))))
 120
 121 (defun cmphess (n x hess f fsum ffun h typx)
 122   (let ((fun (fundata-f ffun))
 123         (arg (fundata-arg ffun))
 124         (cs (fundata-changesign ffun)))
 125     (dotimes (i n)
 126       (let ((hi (if typx (* (aref typx i) h) h))
 127             (old-xi (aref x i)))
 128         (setf (aref hess i i) (/ (- (aref fsum i) (* 2 f)) (* hi hi)))
 129         (dotimes (j i)
 130           (let ((hj (if typx (* (aref typx j) h) h))
 131                 (old-xj (aref x j))
 132                 (f1 0.0)
 133                 (f2 0.0))
 134             (setf (aref x i) (+ old-xi hi))
 135             (setf (aref x j) (+ old-xj hj))
 136             (setf f1 (get-function-value fun x arg cs))
 137             (setf (aref x i) (- old-xi hi))
 138             (setf (aref x j) (- old-xj hj))
 139             (setf f2 (get-function-value fun x arg cs))
 140             (setf (aref x i) old-xi)
 141             (setf (aref x j) old-xj)
 142             (let ((hij (/ (- (+ (* 2 f) f1 f2) (aref fsum i) (aref fsum j))
 143                           (* 2.0 hi hj))))
 144               (setf (aref hess i j) hij)
 145               (setf (aref hess j i) hij))))))))
 146
 147 (defun numgrad (f x &optional scale h all cs)
 148   (let* ((n (length x))
 149          (fd (make-function-data f n cs h scale))
 150          (xv (make-array n :initial-contents x))
 151          (grad (make-array n))
 152          (fv (eval-function fd xv grad nil)))
 153     (coerce grad 'list)))
 154
 155 (defun numhess (f x &optional scale h all cs)
 156   (let* ((n (length x))
 157          (fd (make-function-data f n cs h scale))
 158          (xv (make-array n :initial-contents x))
 159          (grad (make-array n))
 160          (hess (make-array (list n n)))
 161          (fv (eval-function fd xv grad hess)))
 162     (if all (list fv (coerce grad 'list) hess) hess)))
 163
 164
 165 ;;;;
 166 ;;;; Nonlinear optimization modules adapted from Dennis and Schnabel,
 167 ;;;; "Numerical Methods for Unconstrained Optimization and Nonlinear
 168 ;;;; Equations."
 169 ;;;;
 170
 171
 172 ;;;;
 173 ;;;; Mode Internals Data
 174 ;;;;
 175
 176 (defmacro make-internals () `(make-array 7))
 177 (defmacro internals-f (i) `(aref ,i, 0))
 178 (defmacro internals-x (i) `(aref ,i, 1))
 179 (defmacro internals-scale (i) `(aref ,i, 2))
 180 (defmacro internals-delf (i) `(aref ,i, 3))
 181 (defmacro internals-hessf (i) `(aref ,i, 4))
 182 (defmacro internals-ipars (i) `(aref ,i, 5))
 183 (defmacro internals-dpars (i) `(aref ,i, 6))
 184
 185 (defmacro make-ipars () `(make-array 11))
 186 (defmacro ipars-n (ipars) `(aref ,ipars 0))
 187 (defmacro ipars-maxiter (i) `(aref ,i 1))
 188 (defmacro ipars-backtrack (i) `(aref ,i 2))
 189 (defmacro ipars-verbose (i) `(aref ,i 3))
 190 (defmacro ipars-vals-suppl (i) `(aref ,i 4))
 191 (defmacro ipars-count (i) `(aref ,i 5))
 192 (defmacro ipars-termcode (i) `(aref ,i 6))
 193 (defmacro ipars-change-sign (i) `(aref ,i 7))
 194 (defmacro ipars-maxtaken (i) `(aref ,i 8))
 195 (defmacro ipars-consecmax (i) `(aref ,i 9))
 196 (defmacro ipars-retcode (i) `(aref ,i 10))
 197
 198 (defmacro make-dpars () `(make-array 8))
 199 (defmacro dpars-typf(d) `(aref ,d 0))
 200 (defmacro dpars-derivstep (d) `(aref ,d 1))
 201 (defmacro dpars-gradtol (d) `(aref ,d 2))
 202 (defmacro dpars-steptol (d) `(aref ,d 3))
 203 (defmacro dpars-maxstep (d) `(aref ,d 4))
 204 (defmacro dpars-hessadd (d) `(aref ,d 5))
 205 (defmacro dpars-f (d) `(aref ,d 6))
 206 (defmacro dpars-new-f (d) `(aref ,d 7))
 207
 208 (defun newinternals (f x scale h)
 209   (let ((n (length x))
 210         (ip (make-ipars))
 211         (dp (make-dpars))
 212         (internals (make-internals)))
 213     (setf (ipars-n ip) n)
 214     (setf (ipars-maxiter ip ) -1)
 215     (setf (ipars-backtrack ip) t)
 216     (setf (ipars-verbose ip) 0)
 217     (setf (ipars-vals-suppl ip) nil)
 218     (setf (ipars-count ip) 0)
 219     (setf (ipars-termcode ip) 0)
 220     (setf (ipars-change-sign ip) t)
 221
 222     (setf (dpars-typf dp) 1.0)
 223     (setf (dpars-derivstep dp) h)
 224     (setf (dpars-gradtol dp) -1.0)
 225     (setf (dpars-steptol dp) -1.0)
 226     (setf (dpars-maxstep dp) -1.0)
 227     (setf (dpars-hessadd dp) 0.0)
 228
 229     (setf (internals-f internals) f)
 230     (setf (internals-x internals)
 231           (make-array n :element-type 'float :initial-contents x))
 232     (setf (internals-scale internals)
 233           (if scale
 234               (make-array n :element-type 'float :initial-contents scale)
 235               (make-array n :element-type 'float :initial-element 1.0)))
 236     (setf (internals-ipars internals) ip)
 237     (setf (internals-dpars internals) dp)
 238     (setf (internals-delf internals) (make-array n :element-type 'float))
 239     (setf (internals-hessf internals)
 240           (make-array (list n n) :element-type 'float))
 241     internals))
 242
 243
 244 ;;;;
 245 ;;;; Definitions and Globals
 246 ;;;;
 247
 248 (defparameter *init-grad-frac* 0.001)
 249 (defparameter *consec-max-limit* 5)
 250 (defparameter *alpha* 0.0001)
 251 (defparameter *max-step-factor* 1000)
 252 (defparameter *gradtol-power* (/ 1.0 3.0))
 253 (defparameter *steptol-power* (/ 2.0 3.0))
 254 (defparameter *itnlimit* 100)
 255 (defparameter *verbose* 0)
 256 (defparameter *use-search* t)
 257
 258 (defparameter *termcodes*
 259   (vector "not yet terminated"
 260           "gradient size is less than gradient tolerance"
 261           "step size is less than step tolerance"
 262           "no satisfactory step found in backtracking"
 263           "iteration limit exceeded"
 264           "maximum size step taken 5 iterations in a row"))
 265
 266
 267 ;;;;
 268 ;;;; Iteration Data
 269 ;;;;
 270
 271 (defmacro make-iteration () `(make-array 12))
 272 (defmacro iteration-ffun (i) `(aref ,i 0))
 273 (defmacro iteration-x (i) `(aref ,i 1))
 274 (defmacro iteration-new-x (i) `(aref ,i 2))
 275 (defmacro iteration-sx (i) `(aref ,i 3))
 276 (defmacro iteration-delf (i) `(aref ,i 4))
 277 (defmacro iteration-new-delf (i) `(aref ,i 5))
 278 (defmacro iteration-qnstep (i) `(aref ,i 6))
 279 (defmacro iteration-L (i) `(aref ,i 7))
 280 (defmacro iteration-hessf (i) `(aref ,i 8))
 281 (defmacro iteration-work (i) `(aref ,i 9))
 282 (defmacro iteration-ip (i) `(aref ,i 10))
 283 (defmacro iteration-dp (i) `(aref ,i 11))
 284
 285 (defmacro iteration-n (i) `(ipars-n (iteration-ip ,i)))
 286 (defmacro iteration-itnlimit (i) `(ipars-maxiter (iteration-ip ,i)))
 287 (defmacro iteration-backtrack (i) `(ipars-backtrack (iteration-ip ,i)))
 288 (defmacro iteration-verbose (i) `(ipars-verbose (iteration-ip ,i)))
 289 (defmacro iteration-vals-suppl (i) `(ipars-vals-suppl (iteration-ip ,i)))
 290 (defmacro iteration-count (i) `(ipars-count (iteration-ip ,i)))
 291 (defmacro iteration-termcode (i) `(ipars-termcode (iteration-ip ,i)))
 292 (defmacro iteration-maxtaken (i) `(ipars-maxtaken (iteration-ip ,i)))
 293 (defmacro iteration-consecmax (i) `(ipars-consecmax (iteration-ip ,i)))
 294 (defmacro iteration-retcode (i) `(ipars-retcode (iteration-ip ,i)))
 295 (defmacro iteration-change-sign (i) `(ipars-change-sign (iteration-ip ,i)))
 296 (defmacro iteration-typf (i) `(dpars-typf (iteration-dp ,i)))
 297 (defmacro iteration-gradtol (i) `(dpars-gradtol (iteration-dp ,i)))
 298 (defmacro iteration-steptol (i) `(dpars-steptol (iteration-dp ,i)))
 299 (defmacro iteration-maxstep (i) `(dpars-maxstep (iteration-dp ,i)))
 300 (defmacro iteration-hessadd (i) `(dpars-hessadd (iteration-dp ,i)))
 301 (defmacro iteration-f (i) `(dpars-f (iteration-dp ,i)))
 302 (defmacro iteration-new-f (i) `(dpars-new-f (iteration-dp ,i)))
 303
 304
 305 ;;;;
 306 ;;;; Stopping Criteria
 307 ;;;;
 308
 309 (defun gradsize (iter new)
 310   (let ((f (iteration-f iter))
 311         (typf (iteration-typf iter))
 312         (x (iteration-x iter))
 313         (sx (iteration-sx iter))
 314         (delf (if new (iteration-new-delf iter) (iteration-delf iter))))
 315     (uncmin-gradsize delf x sx (max (abs f) typf))))
 316
 317 (defun incrsize (iter)
 318   (let ((n (iteration-n iter))
 319         (new-x (iteration-new-x iter))
 320         (x (iteration-x iter))
 321         (sx (iteration-sx iter))
 322         (work (iteration-work iter)))
 323     (blas-dcopy n x 0 1 work 0 1)
 324     (blas-daxpy n -1.0 new-x 0 1 work 0 1)
 325     (uncmin-maxrelsize work x sx)))
 326
 327 (defun stoptest0 (iter)
 328   (setf (iteration-consecmax iter) 0)
 329   (setf (iteration-termcode iter)
 330         (if (<= (gradsize iter nil)
 331                 (* *init-grad-frac* (iteration-gradtol iter)))
 332             1
 333             0))
 334   (iteration-termcode iter))
 335
 336 (defun stoptest (iter)
 337   (let ((retcode (iteration-retcode iter))
 338         (gradtol (iteration-gradtol iter))
 339         (steptol (iteration-steptol iter))
 340         (count (iteration-count iter))
 341         (itnlimit (iteration-itnlimit iter))
 342         (maxtaken (iteration-maxtaken iter))
 343         (consecmax (iteration-consecmax iter))
 344         (termcode 0))
 345     (cond
 346      ((= retcode 1) (setf termcode 3))
 347      ((<= (gradsize iter t) gradtol) (setf termcode 1))
 348      ((<= (incrsize iter) steptol) (setf termcode 2))
 349      ((>= count itnlimit) (setf termcode 4))
 350      (maxtaken (incf consecmax)
 351                (if (>= consecmax *consec-max-limit*) (setf termcode 5)))
 352      (t (setf consecmax 0)))
 353
 354     (setf (iteration-consecmax iter) consecmax)
 355     (setf (iteration-termcode iter) termcode)
 356     termcode))
 357
 358
 359 ;;;;
 360 ;;;; Function and Derivative Evaluation
 361 ;;;;
 362
 363 (defun eval-funval (iter)
 364   (setf (iteration-f iter)
 365         (eval-function (iteration-ffun iter) (iteration-x iter) nil nil)))
 366
 367 (defun eval-next-funval (iter)
 368   (setf (iteration-new-f iter)
 369         (eval-function (iteration-ffun iter) (iteration-new-x iter) nil nil)))
 370
 371 (defun eval-gradient (iter)
 372   (eval-function (iteration-ffun iter)
 373                  (iteration-x iter)
 374                  (iteration-delf iter)
 375                  nil))
 376
 377 (defun eval-next-gradient (iter)
 378   (eval-function (iteration-ffun iter)
 379                  (iteration-new-x iter)
 380                  (iteration-new-delf iter)
 381                  nil))
 382
 383 (defun eval-hessian (iter)
 384   (eval-function (iteration-ffun iter)
 385                  (iteration-x iter)
 386                  nil
 387                  (iteration-hessf iter)))
 388
 389
 390 ;;;;
 391 ;;;; Backtracking Line Search
 392 ;;;;
 393 ;;;; Modified to quit after *itnlimit* steps in case regular termination
 394 ;;;; doesn't work (because of NaN's).
 395
 396 (defun linesearch (iter)
 397   (let ((n (iteration-n iter))
 398         (x (iteration-x iter))
 399         (new-x (iteration-new-x iter))
 400         (qnstep (iteration-qnstep iter)))
 401     (cond
 402      ((not (iteration-backtrack iter))
 403       (setf (iteration-maxtaken iter) nil)
 404       (blas-dcopy n x 0 1 new-x 0 1)
 405       (blas-daxpy n 1.0 qnstep 0 1 new-x 0 1)
 406       (eval-next-funval iter)
 407       (setf (iteration-retcode iter) 0))
 408      (t
 409       (let ((maxstep (iteration-maxstep iter))
 410             (delf (iteration-delf iter))
 411             (sx (iteration-sx iter))
 412             (work (iteration-work iter)))
 413         (setf (iteration-maxtaken iter) nil)
 414         (setf (iteration-retcode iter) 2)
 415
 416         (dotimes (i n)
 417           (setf (aref work i) (* (aref sx i) (aref qnstep i))))
 418         (let ((newtlen (blas-dnrm2 n work 0 1)))
 419           (when (> newtlen maxstep)
 420                 (dotimes (i n)
 421                   (setf (aref qnstep i)
 422                         (* (aref qnstep i) (/ maxstep newtlen))))
 423                 (setf newtlen maxstep))
 424           (let* ((initslope (blas-ddot n delf 0 1 qnstep 0 1))
 425                  (rellength (uncmin-maxrelsize qnstep x sx))
 426                  (minlambda (if (= rellength 0.0)
 427                                 2.0
 428                                 (/ (iteration-steptol iter) rellength)))
 429                  (lambda 1.0)
 430                  (lambdaprev 1.0)
 431                  (fprev 1.0)
 432                  (count 0))
 433             (loop
 434              (incf count)
 435              (when (<= (iteration-retcode iter) 1) (return))
 436              (blas-dcopy n x 0 1 new-x 0 1)
 437              (blas-daxpy n lambda qnstep 0 1 new-x 0 1)
 438              (eval-next-funval iter)
 439              (cond
 440               ((<= (iteration-new-f iter)
 441                    (+ (iteration-f iter) (* *alpha* lambda initslope)))
 442                (setf (iteration-retcode iter) 0)
 443                (when (and (= lambda 1.0) (> newtlen (* 0.99 maxstep)))
 444                      (setf (iteration-maxtaken iter) t)))
 445               ((or (< lambda minlambda) (>= count *itnlimit*))
 446                (setf (iteration-retcode iter) 1)
 447                (setf (iteration-new-f iter) (iteration-f iter))
 448                (blas-dcopy n x 0 1 new-x 0 1))
 449               (t
 450                (setf lambdatemp
 451                      (uncmin-linesearch (iteration-f iter) initslope
 452                                         lambda (iteration-new-f iter)
 453                                         lambdaprev fprev))
 454                (setf lambdaprev lambda)
 455                (setf fprev (iteration-new-f iter))
 456                (setf lambda (max (* 0.1 lambda) lambdatemp))
 457                (when (> (iteration-verbose iter) 0)
 458                      (format t "Backtracking: lambda = ~,6g~%"
 459                              lambda))))))))))))
 460
 461
 462 ;;;;
 463 ;;;; Status Printing Functions
 464 ;;;;
 465
 466 (defun print-header (iter)
 467   (when (< 0 (iteration-verbose iter))
 468         (format t "Iteration ~d.~%" (iteration-count iter))))
 469
 470 (defun print-status (iter)
 471   (let ((n (iteration-n iter))
 472         (f (iteration-f iter))
 473         (x (iteration-x iter))
 474         (delf (iteration-delf iter))
 475         (hessf (iteration-hessf iter))
 476         (verbose (iteration-verbose iter))
 477         (change-sign (iteration-change-sign iter))
 478         (tcode (iteration-termcode iter)))
 479     (when (< 0 verbose)
 480           (format t "Criterion value = ~,6g~%" (if change-sign (- f) f))
 481           (when (< 1 verbose)
 482                 (format t "Location = <")
 483                 (dotimes (i n)
 484                   (format t (if (< i (- n 1)) "~,6g " "~,6g>~%") (aref x i)))
 485                 (when (< 2 verbose)
 486                       (format t "Gradient = <")
 487                       (dotimes (i n)
 488                         (let ((di (aref delf i)))
 489                           (format t
 490                                   (if (< i (- n 1)) "~,6g " "~,6g>~%")
 491                                   (if change-sign (- di) di))))
 492                       (when (< 3 verbose)
 493                             (format t "Hessian:~%")
 494                             (dotimes (i n)
 495                               (dotimes (j n)
 496                                 (let ((hij (aref hess (+ (* i n) j))))
 497                                   (format t
 498                                           (if (< j (- n 1)) "~,6g " "~,6g~%")
 499                                           (if change-sign (- hij) hij))))))))
 500           (when (/= tcode 0)
 501                 (format t
 502                         "Reason for termination: ~a.~%"
 503                         (aref *termcodes* tcode))))))
 504
 505
 506 ;;;;
 507 ;;;; Iteration Driver
 508 ;;;;
 509
 510 (defun findqnstep (iter)
 511   (let ((sx (iteration-sx iter))
 512         (hessf (iteration-hessf iter))
 513         (L (iteration-L iter))
 514         (delf (iteration-delf iter))
 515         (qnstep (iteration-qnstep iter)))
 516     (setf (iteration-hessadd iter) (uncmin-modelhess sx hessf L))
 517     (uncmin-cholsolve delf L  qnstep)))
 518
 519 (defun iterupdate (iter)
 520   (let ((n (iteration-n iter)))
 521     (setf (iteration-f iter) (iteration-new-f iter))
 522     (blas-dcopy n (iteration-new-x iter) 0 1 (iteration-x iter) 0 1)
 523     (blas-dcopy n (iteration-new-delf iter) 0 1 (iteration-delf iter) 0 1)))
 524
 525
 526 ;;;;
 527 ;;;; External Interface Routines
 528 ;;;;
 529
 530 (defun mindriver (iter trfun)
 531   (setf (iteration-consecmax iter) 0)
 532   (setf (iteration-count iter) 0)
 533   (setf (iteration-termcode iter) 0)
 534   (unless (iteration-vals-suppl iter)
 535           (eval-funval iter)
 536           (eval-gradient iter)
 537           (eval-hessian iter))
 538
 539   (stoptest0 iter)
 540   (print-header iter)
 541   (print-status iter)
 542   (loop
 543    (when (/= 0 (iteration-termcode iter)) (return))
 544    (incf (iteration-count iter))
 545    (print-header iter)
 546    (findqnstep iter)
 547    (linesearch iter)
 548    (eval-next-gradient iter)
 549    (stoptest iter)
 550    (iterupdate iter)
 551    (eval-hessian iter)
 552    (print-status iter)
 553    (when trfun (funcall trfun))))
 554
 555 (defun minresultstring (code)
 556   (if (<= 0 code (- (length *termcodes*) 1))
 557       (aref *termcodes* code)
 558       "unknown return code"))
 559
 560 (defun minsetup (ffun x typx ip dp delf hessf)
 561   (let ((iter (make-iteration))
 562         (n (ipars-n ip)))
 563     (setf (iteration-ip iter) ip)
 564     (setf (iteration-dp iter) dp)
 565
 566     (setf (iteration-ffun iter) ffun)
 567     (setf (iteration-x iter) x)
 568     (setf (iteration-new-x iter) (make-array n :element-type 'float))
 569     (setf (iteration-sx iter)
 570           (make-array n :element-type 'float :initial-element 1.0))
 571     (setf (iteration-delf iter) delf)
 572     (setf (iteration-new-delf iter) (make-array n :element-type 'float))
 573     (setf (iteration-qnstep iter) (make-array n :element-type 'float))
 574
 575     (when typx
 576           (let ((sx (iteration-sx iter)))
 577             (dotimes (i n)
 578               (let ((txi (aref typx i)))
 579                 (if (< 0.0 txi)
 580                     (setf (aref sx i) (/ 1.0 txi)))))))
 581
 582     (setf (iteration-L iter) (make-array (list n n) :element-type 'float))
 583     (setf (iteration-hessf iter) hessf)
 584     (setf (iteration-work iter) (make-array n :element-type 'float))
 585
 586     (when (< (iteration-verbose iter) 0)
 587           (setf (iteration-verbose iter) *verbose*))
 588     (when (< (iteration-itnlimit iter) 0)
 589           (setf (iteration-itnlimit iter) *itnlimit*))
 590
 591     (when (<= (iteration-typf iter) 0.0)
 592           (setf (iteration-typf iter) 1.0))
 593     (when (<= (iteration-gradtol iter) 0.0)
 594           (setf (iteration-gradtol iter) (^ machine-epsilon *gradtol-power*)))
 595     (when (<= (iteration-steptol iter) 0.0)
 596           (setf (iteration-steptol iter) (^ machine-epsilon *steptol-power*)))
 597
 598     (when (<= (iteration-maxstep iter) 0.0)
 599           (let ((dx (iteration-x iter))
 600                 (sx (iteration-sx iter))
 601                 (w (iteration-work iter)))
 602             (dotimes (i n)
 603               (setf (aref w i) (* (aref dx i) (aref sx i))))
 604             (let ((nx0 (blas-dnrm2 n w 0 1))
 605                   (nsx (blas-dnrm2 n sx 0 1)))
 606               (setf (iteration-maxstep iter)
 607                     (* *max-step-factor* (max nx0 nsx))))))
 608     iter))
 609 ;;;;
 610 ;;;; Mode Info Accessors
 611 ;;;;
 612
 613 (defmacro mode-info-f () `(aref (slot-value 'internals) 0))
 614 (defmacro mode-info-x () `(aref (slot-value 'internals) 1))
 615 (defmacro mode-info-scale () `(aref (slot-value 'internals) 2))
 616 (defmacro mode-info-delf () `(aref (slot-value 'internals) 3))
 617 (defmacro mode-info-hessf () `(aref (slot-value 'internals) 4))
 618 (defmacro mode-info-ipars () `(aref (slot-value 'internals) 5))
 619 (defmacro mode-info-dpars () `(aref (slot-value 'internals) 6))
 620
 621 (defmacro ipars-n (ipars) `(aref ,ipars 0))
 622 (defmacro ipars-maxiter (ipars) `(aref ,ipars 1))
 623 (defmacro ipars-backtrack (ipars) `(aref ,ipars 2))
 624 (defmacro ipars-verbose (ipars) `(aref ,ipars 3))
 625
 626 (defmacro dpars-derivstep (dpars) `(aref ,dpars 1))
 627 (defmacro dpars-f (dpars) `(aref ,dpars 6))
 628
 629
 630 ;;;;
 631 ;;;; Mode Info Prototype
 632 ;;;;
 633
 634 (defproto minfo-proto '(internals))
 635
 636 (defmeth minfo-proto :isnew (f x &key scale (derivstep -1.0))
 637   (setf (slot-value 'internals) (newinternals f x scale derivstep))
 638   nil)
 639
 640 (defmeth minfo-proto :maximize (&optional verbose trfun)
 641   (let* ((internals (slot-value 'internals))
 642          (f (internals-f internals))
 643          (x (internals-x internals))
 644          (scale (internals-scale internals))
 645          (ip (internals-ipars internals))
 646          (dp (internals-dpars internals))
 647          (delf (internals-delf internals))
 648          (hessf (internals-hessf internals))
 649          (n (ipars-n ip))
 650          (func (make-function-data f n t (dpars-derivstep dp) scale))
 651          (iter (minsetup func x scale ip dp delf hessf)))
 652     (when (and verbose (integerp verbose))
 653           (setf (ipars-verbose ip) verbose))
 654     (when (> (ipars-verbose ip) 0)
 655           (format t "maximizing...~%"))
 656     (mindriver iter trfun)
 657     (setf (ipars-vals-suppl ip) t)
 658     (minresultstring (ipars-termcode ip))))
 659
 660 (defmeth minfo-proto :x () (coerce (mode-info-x) 'list))
 661 (defmeth minfo-proto :scale () (coerce (mode-info-scale) 'list))
 662 (defmeth minfo-proto :derivstep () (dpars-derivstep (mode-info-dpars)))
 663
 664 (defmeth minfo-proto :f (&optional (val nil set))
 665   (when set
 666         (send self :set-no-vals-supplied)
 667         (setf (mode-info-f) val))
 668   (mode-info-f))
 669
 670 (defmeth minfo-proto :fvals ()
 671   (let* ((n (ipars-n (mode-info-ipars)))
 672          (val (dpars-f (mode-info-dpars)))
 673          (grad (coerce (mode-info-delf) 'list))
 674          (hess (make-array (list n n))))
 675     (replace (compound-data-seq hess)
 676              (compound-data-seq (mode-info-hessf)))
 677     (list val grad hess)))
 678
 679 (defmeth minfo-proto :copy ()
 680   (let ((obj (make-object minfo-proto))
 681         (internals (copy-seq (slot-value 'internals))))
 682     (dotimes (i (length internals))
 683              (let ((x (aref internals i)))
 684                (if (sequencep x)
 685                    (setf (aref internals i) (copy-seq x)))))
 686     (send obj :add-slot 'internals internals)
 687     obj))
 688
 689 (defmeth minfo-proto :derivscale ()
 690   (let* ((x (send self :x))
 691          (step (^ machine-epsilon (/ 1 6)))
 692          (hess (numhess (send self :f) x (send self :scale) step))
 693          (scale (pmax (abs x) (sqrt (abs (/ (diagonal hess)))))))
 694     (setf hess (numhess (send self :f) x scale step))
 695     (setf scale (pmax (abs x) (sqrt (abs (/ (diagonal hess))))))
 696     (setf (mode-info-scale) (coerce scale '(vector float)))
 697     (setf (dpars-derivstep (mode-info-dpars)) step)))
 698
 699 (defmeth minfo-proto :verbose (&optional (val nil set))
 700   (when set
 701         (setf (ipars-verbose (mode-info-ipars))
 702               (cond ((integerp val) val)
 703                     ((null val) 0)
 704                     (t 1))))
 705   (ipars-verbose (mode-info-ipars)))
 706
 707 (defmeth minfo-proto :backtrack (&optional (val nil set))
 708   (if set (setf (ipars-backtrack (mode-info-ipars)) (if val 1 0)))
 709   (ipars-backtrack (mode-info-ipars)))
 710
 711 (defmeth minfo-proto :maxiter (&optional (val nil set))
 712     (if set (setf (ipars-maxiter (mode-info-ipars))
 713                   (if (integerp val) val -1)))
 714     (ipars-maxiter (mode-info-ipars)))
 715
 716
 717 ;;;;
 718 ;;;;
 719 ;;;; Newton's Method with Backtracking
 720 ;;;;
 721 ;;;;
 722
 723 (defun newtonmax (f start &key
 724                     scale
 725                     (derivstep -1.0)
 726                     (count-limit -1)
 727                     (verbose 1)
 728                     return-derivs
 729                     trace)
 730 "Args:(f start &key scale derivstep (verbose 1) return-derivs)
 731 Maximizes F starting from START using Newton's method with backtracking.
 732 If RETURN-DERIVS is NIL returns location of maximum; otherwise returns
 733 list of location, unction value, gradient and hessian at maximum.
 734 SCALE should be a list of the typical magnitudes of the parameters.
 735 DERIVSTEP is used in numerical derivatives and VERBOSE controls printing
 736 of iteration information. COUNT-LIMIT limits the number of iterations"
 737   (let ((verbose (if verbose (if (integerp verbose) verbose 1) 0))
 738         (minfo (send minfo-proto :new f start
 739                      :scale scale :derivstep derivstep)))
 740     (send minfo :maxiter count-limit)
 741     (send minfo :derivscale)
 742     (send minfo :maximize
 743           verbose
 744           (if trace #'(lambda () (funcall trace minfo))))
 745     (if return-derivs
 746         (cons (send minfo :x) (- (send minfo :fvals)))
 747         (send minfo :x))))
 748
 749 ;;;;
 750 ;;;;
 751 ;;;; Nelder-Mead Simplex Method
 752 ;;;;
 753 ;;;;
 754
 755 (defun nelmeadmax (f start &key
 756                      (size 1)
 757                      (epsilon (sqrt machine-epsilon))
 758                      (count-limit 500)
 759                      (verbose t)
 760                      (alpha 1.0)
 761                      (beta 0.5)
 762                      (gamma 2.0)
 763                      (delta 0.5))
 764 "Args: (f start &key (size 1) (epsilon (sqrt machine-epsilon))
 765           (count-limit 500) (verbose t) alpha beta gamma delta)
 766 Maximizes F using the Nelder-Mead simplex method. START can be a
 767 starting simplex - a list of N+1 points, with N=dimension of problem,
 768 or a single point. If start is a single point you should give the
 769 size of the initial simplex as SIZE, a sequence of length N. Default is
 770 all 1's. EPSILON is the convergence tolerance. ALPHA-DELTA can be used to
 771 control the behavior of simplex algorithm."
 772     (let ((s (send simplex-proto :new f start size)))
 773       (do ((best (send s :best-point) (send s :best-point))
 774            (count 0 (+ count 1))
 775            next)
 776           ((or (< (send s :relative-range) epsilon) (>= count count-limit))
 777            (if (and verbose (>= count count-limit))
 778                (format t "Iteration limit exceeded.~%"))
 779            (send s :point-location (send s :best-point)))
 780           (setf next (send s :extrapolate-from-worst (- alpha)))
 781           (if (send s :is-worse best next)
 782               (setf next (send s :extrapolate-from-worst gamma))
 783               (when (send s :is-worse next (send s :second-worst-point))
 784                     (setf next (send s :extrapolate-from-worst beta))
 785                     (if (send s :is-worse next (send s :worst-point))
 786                         (send s :shrink-to-best delta))))
 787           (if verbose
 788               (format t "Value = ~g~%"
 789                       (send s :point-value (send s :best-point)))))))
 790
 791
 792 ;;;
 793 ;;; Simplex Prototype
 794 ;;;
 795
 796 (defproto simplex-proto '(f simplex))
 797
 798 ;;;
 799 ;;; Simplex Points
 800 ;;;
 801
 802 (defmeth simplex-proto :make-point (x)
 803   (let ((f (send self :f)))
 804     (if f
 805         (let ((val (funcall f x)))
 806           (cons (if (consp val) (car val) val) x))
 807         (cons nil x))))
 808
 809 (defmeth simplex-proto :point-value (x) (car x))
 810
 811 (defmeth simplex-proto :point-location (x) (cdr x))
 812
 813 (defmeth simplex-proto :is-worse (x y)
 814   (< (send self :point-value x) (send self :point-value y)))
 815
 816 ;;;
 817 ;;; Making New Simplices
 818 ;;;
 819
 820 (defmeth simplex-proto :isnew (f start &optional size)
 821   (send self :simplex start size)
 822   (send self :f f))
 823
 824 ;;;
 825 ;;; Slot Accessors and Mutators
 826 ;;;
 827
 828 (defmeth simplex-proto :simplex (&optional new size)
 829   (if new
 830       (let ((simplex
 831              (if (and (consp new) (sequencep (car new)))
 832                  (if (/= (length new) (+ 1 (length (car new))))
 833                      (error "bad simplex data")
 834                      (copy-list new))
 835                  (let* ((n (length new))
 836                         (size (if size size (repeat 1 n)))
 837                      ;   (pts (- (* 2 (uniform-rand (repeat n (+ n 1)))) 1)))
 838                         (diag (* 2 size (- (random (repeat 2 n)) .5)))
 839                         (pts (cons (repeat 0 n)
 840                                    (mapcar #'(lambda (x) (coerce x 'list))
 841                                            (column-list (diagonal diag))))))
 842                    (mapcar #'(lambda (x) (+ (* size x) new)) pts)))))
 843         (setf (slot-value 'simplex)
 844               (mapcar #'(lambda (x) (send self :make-point x)) simplex))
 845         (send self :sort-simplex)))
 846   (slot-value 'simplex))
 847
 848 (defmeth simplex-proto :f (&optional f)
 849   (when f
 850         (setf (slot-value 'f) f)
 851         (let ((simplex
 852                (mapcar #'(lambda (x) (send self :point-location x))
 853                        (send self :simplex))))
 854           (send self :simplex simplex)))
 855   (slot-value 'f))
 856
 857 (defmeth simplex-proto :sort-simplex ()
 858   (if (send self :f)
 859       (setf (slot-value 'simplex)
 860             (sort (slot-value 'simplex)
 861                   #'(lambda (x y) (send self :is-worse x y))))))
 862
 863 ;;;
 864 ;;; Other Methods Using List Representation of SImplex
 865 ;;;
 866
 867 (defmeth simplex-proto :best-point () (car (last (send self :simplex))))
 868 (defmeth simplex-proto :worst-point () (first (send self :simplex)))
 869 (defmeth simplex-proto :second-worst-point () (second (send self :simplex)))
 870 (defmeth simplex-proto :replace-point (new old)
 871   (let* ((simplex (send self :simplex))
 872          (n (position old simplex)))
 873     (when n
 874           (setf (nth n simplex) new)
 875           (send self :sort-simplex))))
 876 (defmeth simplex-proto :mean-opposite-face (x)
 877   (let ((face (mapcar #'(lambda (x) (send self :point-location x))
 878                       (remove x (send self :simplex)))))
 879     (/ (apply #'+ face) (length face))))
 880
 881 ;;;
 882 ;;; Iteration Step Methods
 883 ;;;
 884
 885 (defmeth simplex-proto :extrapolate-from-worst (fac)
 886   (let* ((worst (send self :worst-point))
 887          (wloc (send self :point-location worst))
 888          (delta (- (send self :mean-opposite-face worst) wloc))
 889          (new (send self :make-point (+ wloc (* (- 1 fac) delta)))))
 890     (if (send self :is-worse worst new) (send self :replace-point new worst))
 891     new))
 892
 893 (defmeth simplex-proto :shrink-to-best (fac)
 894   (let* ((best (send self :best-point))
 895          (bloc (send self :point-location best)))
 896     (dolist (x (copy-list (send self :simplex)))
 897             (if (not (eq x best))
 898                 (send self :replace-point
 899                       (send self :make-point
 900                             (+ bloc
 901                                (* fac
 902                                   (- (send self :point-location x) bloc))))
 903                       x)))))
 904
 905 (defmeth simplex-proto :relative-range ()
 906   (let ((best (send self :point-value (send self :best-point)))
 907         (worst (send self :point-value (send self :worst-point))))
 908     (* 2 (/ (abs (- best worst)) (+ 1 (abs best) (abs worst))))))