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-ffi-int
  11        :lisp-stat-object-system
  12        :lisp-stat-types
  13        :lisp-stat-sequence
  14        :lisp-stat-compound-data
  15        :lisp-stat-math
  16        :lisp-stat-basics
  17        :lisp-stat-matrix
  18        :lisp-stat-linalg-data
  19        :lisp-stat-linalg)
  20  (:shadowing-import-from :lisp-stat-object-system
  21                          slot-value call-method call-next-method)
  22  (:export
  23      ;; derivatives
  24      numgrad numhess
  25
  26      ;; optimization
  27      newtonmax nelmeadmax))
  28
  29
  30
  31 ;;; FIXME:AJR: There is a need to figure out the proper symbols to
  32 ;;; export.  more importantly should there be any specialty package
  33 z;;; that are exported for maximization?
  34
  35 (in-package :lisp-stat-optimize)
  36
  37 (defvar *maximize-callback-function* nil
  38   "Used in generic optimization to determine function name -- symbol or string?")
  39
  40 (defvar *maximize-callback-arg* nil
  41   "args to function to maximize")
  42
  43
  44 ;;;
  45 ;;; CFFI support using library for optimization work.
  46 ;;;
  47
  48 ;; There is a problem with this particular approach, in terms of
  49 ;; circular dependencies.  We can not have this out-of-object call
  50 ;; into optimize, at least not from here.
  51 (cffi:defcallback ccl-maximize-callback :void ((n :int)
  52                                                (px :pointer)
  53                                                (pfval :pointer)
  54                                                (pgrad :pointer)
  55                                                (phess :pointer)
  56                                                (pderivs :pointer))
  57    (lisp-stat-optimize::maximize-callback n px pfval pgrad phess pderivs))
  58
  59 (cffi:defcfun ("register_maximize_callback" register-maximize-callback)
  60     :void (x :pointer))
  61 (register-maximize-callback (cffi:callback ccl-maximize-callback))
  62
  63 (cffi:defcfun ("ccl_numgrad_front" ccl-numgrad-front)
  64     :int (x :int) (y :pointer) (z :pointer) (u :double) (v :pointer))
  65 (defun numgrad-front (x y z u v)
  66   (ccl-numgrad-front x y z (float u 1d0) v))
  67
  68 (cffi:defcfun ("ccl_numhess_front" ccl-numhess-front)
  69     :int (x :int) (y :pointer) (z :pointer) (u :pointer) (v :pointer) (w :double) (a :pointer))
  70 (defun numhess-front (x y z u v w a)
  71   (ccl-numhess-front x y z u v (float w 1d0) a))
  72
  73 (cffi:defcfun ("ccl_minfo_maximize" ccl-minfo-maximize)
  74     :int (x :pointer) (y :pointer) (z :pointer) (u :pointer) (v :pointer) (w :int))
  75 (defun base-minfo-maximize (x y z u v w)
  76   (ccl-minfo-maximize x y z u v w))
  77
  78
  79
  80 ;;;;
  81 ;;;; minfo basics (internal??)
  82 ;;;;
  83
  84 (defun init-minfo-ipar-values (n ipars &key
  85                                (TRUE 1)
  86                                (FALSE 0)
  87                                (k 0)
  88                                (m 0)
  89                                (itnlimit -1)
  90                                (backtrack TRUE)
  91                                (verbose 0)
  92                                (vals_suppl FALSE)
  93                                (exptilt TRUE)
  94                                (count 0)
  95                                (termcode 0))
  96   "Initialize ipars (iteration parameters) by destructive modification."
  97   (setf (aref ipars 0) n)
  98   (setf (aref ipars 1) m)
  99   (setf (aref ipars 2) k)
 100   (setf (aref ipars 3) itnlimit)
 101   (setf (aref ipars 4) backtrack)
 102   (setf (aref ipars 5) verbose)
 103   (setf (aref ipars 6) vals_suppl)
 104   (setf (aref ipars 7) exptilt)
 105   (setf (aref ipars 8) count)
 106   (setf (aref ipars 9) termcode))
 107
 108 (defun init-minfo-dpar-values (h dpars &key
 109                                (typf 1.0)
 110                                (gradtol -1.0)
 111                                (steptol -1.0)
 112                                (maxstep -1.0)
 113                                (dflt 0.0)
 114                                (tilt 0.0)
 115                                (newtilt 0.0)
 116                                (hessadd 0.0))
 117   "Initialize dpars (derivative parameters) by destructive modification."
 118   (setf (aref dpars 0) typf)
 119   (setf (aref dpars 1) h)
 120   (setf (aref dpars 2) gradtol)
 121   (setf (aref dpars 3) steptol)
 122   (setf (aref dpars 4) maxstep)
 123   (setf (aref dpars 5) dflt)
 124   (setf (aref dpars 6) tilt)
 125   (setf (aref dpars 7) newtilt)
 126   (setf (aref dpars 8) hessadd))
 127
 128 (defun init-minfo-internals (n h internals)
 129   (let ((ipars (aref internals 8))
 130         (dpars (aref internals 9)))
 131     (init-minfo-ipar-values n ipars)
 132     (init-minfo-dpar-values h dpars)))
 133
 134 (defun new-minfo-internals (f x &key scale ((:derivstep h) -1.0))
 135   (check-sequence x)
 136   (check-real x)
 137   (check-one-real h)
 138   (let ((n (length x)))
 139     (when scale
 140           (check-sequence scale)
 141           (check-real scale)
 142           (if (/= n (length scale)) (error "scale and x not the same length")))
 143     (let ((internals (make-array 12)))
 144       (setf (aref internals 0) f)
 145       (setf (aref internals 3) (if (consp x) (copy-list x) (coerce x 'list)))
 146       (setf (aref internals 4)
 147             (if scale (copy-seq scale) (make-array n :initial-element 1.0)))
 148       (setf (aref internals 5) (make-list (+ 1 n (* n n))))
 149       (setf (aref internals 8) (make-array 10))
 150       (setf (aref internals 9) (make-array 9))
 151       (init-minfo-internals n h internals)
 152       internals)))
 153
 154 (defun minfo-maximize (internals &optional verbose)
 155   "This function does what?"
 156   (let* ((f (aref internals 0))
 157          (x (aref internals 3))
 158          (fvals (aref internals 5))
 159          (n (length x))
 160          (v (if verbose (if (integerp verbose) verbose 1) -1)))
 161     (setf (aref internals 3) (copy-list x))
 162     (setf (aref internals 5) (copy-list fvals))
 163     (let ((*maximize-callback-function* f)
 164           (*maximize-callback-arg* (make-list n)))
 165       (let* ((x (aref internals 3))
 166              (scale (aref internals 4))
 167              (fvals (aref internals 5))
 168              (ip (aref internals 8))
 169              (dp (aref internals 9))
 170              (px (la-data-to-vector x mode-re))
 171              (pscale (la-data-to-vector scale mode-re))
 172              (pfvals (la-vector (length fvals) mode-re))
 173              (pip (la-data-to-vector ip mode-in))
 174              (pdp (la-data-to-vector dp mode-re)))
 175         (unwind-protect
 176             (progn
 177               (base-minfo-maximize px pfvals pscale pip pdp v)) ;; access to C
 178           (la-vector-to-data px n mode-re x)
 179           (la-vector-to-data pfvals (+ 1 n (* n n)) mode-re fvals)
 180           (la-vector-to-data pip (length ip) mode-in ip)
 181           (la-vector-to-data pdp (length dp) mode-re dp))
 182         (get-buf)))))
 183
 184
 185
 186 ;;;;
 187 ;;;; Mode Info Prototype
 188 ;;;;
 189
 190 (defproto minfo-proto '(internals))
 191
 192 #+xlisp (send minfo-proto :add-method :isnew #'|minfo-isnew|)
 193 #+xlisp (send minfo-proto :add-method :maximize #'|minfo-maximize|)
 194 #+xlisp (send minfo-proto :add-method :loglaplace #'|minfo-loglap|)
 195 #-xlisp
 196 (defmeth minfo-proto :isnew (&rest args)
 197   (setf (slot-value 'internals) (apply #'new-minfo-internals args)))
 198 #-xlisp
 199 (defmeth minfo-proto :maximize (&rest args)
 200   (apply #'minfo-maximize (slot-value 'internals) args))
 201
 202 (defmeth minfo-proto :x () (aref (slot-value 'internals) 3))
 203 (defmeth minfo-proto :scale () (aref (slot-value 'internals) 4))
 204 (defmeth minfo-proto :derivstep () (aref (aref (slot-value 'internals) 9) 1))
 205 (defmeth minfo-proto :tilt () (aref (aref (slot-value 'internals) 9) 6))
 206
 207 (defmeth minfo-proto :f (&optional (val nil set))
 208   (when set
 209         (send self :set-no-vals-supplied)
 210         (setf (aref (slot-value 'internals) 0) val))
 211   (aref (slot-value 'internals) 0))
 212
 213 (defmeth minfo-proto :set-no-vals-supplied ()
 214   (setf (aref (aref (slot-value 'internals) 8) 6) 0))
 215
 216 (defmeth minfo-proto :exptilt (&optional (val nil set))
 217   (if set
 218       (let ((old (send self :exptilt)))
 219         (setf (aref (aref (slot-value 'internals) 8) 7) (if val 1 0))
 220         (if (and (not (or (and old val) (and (not old) (not val))))
 221                  (/= (send self :tilt) 0.0))
 222             (send self :set-no-vals-supplied))))
 223   (= 1 (aref (aref (slot-value 'internals) 8) 7)))
 224
 225 (defmeth minfo-proto :newtilt (&optional (val nil set))
 226   (when set
 227         (setf (aref (aref (slot-value 'internals) 9) 7) (float val))
 228         (if (/= (send self :tilt) 0.0) (send self :set-no-vals-supplied)))
 229   (aref (aref (slot-value 'internals) 9) 7))
 230
 231 (defmeth minfo-proto :gfuns (&optional (val nil set))
 232   (when set
 233         (if (or (not (consp val))
 234                 (not (every #'functionp val)))
 235             (error "not all functions"))
 236         (setf (aref (slot-value 'internals) 1) val)
 237         (setf (aref (aref (slot-value 'internals) 8) 1) (length val))
 238         (setf (aref (slot-value 'internals) 10) (repeat 1.0 (length val)))
 239         (if (/= (send self :tilt) 0.0) (send self :set-no-vals-supplied)))
 240   (aref (slot-value 'internals) 1))
 241
 242 (defmeth minfo-proto :cfuns (&optional (val nil set))
 243   (when set
 244         (if (or (not (consp val))
 245                 (not (every #'functionp val)))
 246             (error "not all functions"))
 247         (setf (aref (slot-value 'internals) 2) val)
 248         (setf (aref (aref (slot-value 'internals) 8) 2) (length val))
 249         (setf (aref (slot-value 'internals) 7) (repeat 0.0 (length val)))
 250         (setf (aref (slot-value 'internals) 11) (repeat 0.0 (length val)))
 251         (send self :set-no-vals-supplied))
 252   (aref (slot-value 'internals) 2))
 253
 254 (defmeth minfo-proto :ctarget (&optional (val nil set))
 255   (when set
 256         (if (/= (length val) (length (send self :ctarget)))
 257             (error "bad target length"))
 258         (setf (aref (slot-value 'internals) 7) val))
 259   (aref (slot-value 'internals) 7))
 260
 261 (defmeth minfo-proto :fvals ()
 262   (let* ((fv (aref (slot-value 'internals) 5))
 263          (n (length (send self :x)))
 264          (val (select fv 0))
 265          (grad (select fv (iseq 1 n)))
 266          (hess (matrix (list n n) (select fv (iseq (+ 1 n) (+ n (* n n)))))))
 267     (list val grad hess)))
 268
 269 (defmeth minfo-proto :copy ()
 270   (let ((obj (make-object minfo-proto))
 271         (internals (copy-seq (slot-value 'internals))))
 272     (dotimes (i (length internals))
 273              (let ((x (aref internals i)))
 274                (if (sequencep x)
 275                    (setf (aref internals i) (copy-seq x)))))
 276     (send obj :add-slot 'internals internals)
 277     obj))
 278
 279 (defmeth minfo-proto :derivscale ()
 280   (let* ((step (^ machine-epsilon (/ 1 6)))
 281          (hess (numhess (send self :f) (send self :x) (send self :scale) step))
 282          (scale (pmax (abs (send self :x)) (sqrt (abs (/ (diagonal hess)))))))
 283     (setf hess (numhess (send self :f) (send self :x) scale step))
 284     (setf scale (pmax (abs (send self :x)) (sqrt (abs (/ (diagonal hess))))))
 285     (setf (aref (slot-value 'internals) 4) scale)
 286     (setf (aref (aref (slot-value 'internals) 9) 1) step)))
 287
 288 (defmeth minfo-proto :verbose (&optional (val nil set))
 289   (when set
 290         (setf (aref (aref (slot-value 'internals) 8) 5)
 291               (cond ((integerp val) val)
 292                     ((null val) 0)
 293                     (t 1))))
 294   (aref (aref (slot-value 'internals) 8) 5))
 295
 296 (defmeth minfo-proto :backtrack (&optional (val nil set))
 297   (if set (setf (aref (aref (slot-value 'internals) 8) 4) (if val 1 0)))
 298   (aref (aref (slot-value 'internals) 8) 4))
 299
 300 (defmeth minfo-proto :maxiter (&optional (val nil set))
 301     (if set (setf (aref (aref (slot-value 'internals) 8) 3)
 302                   (if (integerp val) val -1)))
 303     (aref (aref (slot-value 'internals) 8) 3))
 304
 305 (defmeth minfo-proto :tiltscale (&optional (val nil set))
 306   (when set
 307         (if (/= (length val) (length (send self :gfuns)))
 308             (error "wrong size tilt scale sequence"))
 309         (setf (aref (slot-value 'internals) 10) val))
 310   (aref (slot-value 'internals) 10))
 311
 312 ;;;;
 313 ;;;;
 314 ;;;; Newton's Method with Backtracking
 315 ;;;;
 316 ;;;;
 317
 318 (defun newtonmax (f start &key
 319                     scale
 320                     (derivstep -1.0)
 321                     (count-limit -1)
 322                     (verbose 1)
 323                     return-derivs)
 324 "Args:(f start &key scale derivstep (verbose 1) return-derivs)
 325 Maximizes F starting from START using Newton's method with backtracking.
 326 If RETURN-DERIVS is NIL returns location of maximum; otherwise returns
 327 list of location, unction value, gradient and hessian at maximum.
 328 SCALE should be a list of the typical magnitudes of the parameters.
 329 DERIVSTEP is used in numerical derivatives and VERBOSE controls printing
 330 of iteration information. COUNT-LIMIT limits the number of iterations"
 331   (let ((verbose (if verbose (if (integerp verbose) verbose 1) 0))
 332         (minfo (send minfo-proto :new f start
 333                      :scale scale :derivstep derivstep)))
 334     (send minfo :maxiter count-limit)
 335     (send minfo :derivscale)
 336     (send minfo :maximize verbose)
 337     (if return-derivs
 338         (cons (send minfo :x) (- (send minfo :fvals)))
 339         (send minfo :x))))
 340
 341 ;;;;
 342 ;;;;
 343 ;;;; Nelder-Mead Simplex Method
 344 ;;;;
 345 ;;;;
 346
 347 (defun nelmeadmax (f start &key
 348                      (size 1)
 349                      (epsilon (sqrt machine-epsilon))
 350                      (count-limit 500)
 351                      (verbose t)
 352                      (alpha 1.0)
 353                      (beta 0.5)
 354                      (gamma 2.0)
 355                      (delta 0.5))
 356 "Args: (f start &key (size 1) (epsilon (sqrt machine-epsilon))
 357           (count-limit 500) (verbose t) alpha beta gamma delta)
 358 Maximizes F using the Nelder-Mead simplex method. START can be a
 359 starting simplex - a list of N+1 points, with N=dimension of problem,
 360 or a single point. If start is a single point you should give the
 361 size of the initial simplex as SIZE, a sequence of length N. Default is
 362 all 1's. EPSILON is the convergence tolerance. ALPHA-DELTA can be used to
 363 control the behavior of simplex algorithm."
 364     (let ((s (send simplex-proto :new f start size)))
 365       (do ((best (send s :best-point) (send s :best-point))
 366            (count 0 (+ count 1))
 367            next)
 368           ((or (< (send s :relative-range) epsilon) (>= count count-limit))
 369            (if (and verbose (>= count count-limit))
 370                (format t "Iteration limit exceeded.~%"))
 371            (send s :point-location (send s :best-point)))
 372           (setf next (send s :extrapolate-from-worst (- alpha)))
 373           (if (send s :is-worse best next)
 374               (setf next (send s :extrapolate-from-worst gamma))
 375               (when (send s :is-worse next (send s :second-worst-point))
 376                     (setf next (send s :extrapolate-from-worst beta))
 377                     (if (send s :is-worse next (send s :worst-point))
 378                         (send s :shrink-to-best delta))))
 379           (if verbose
 380               (format t "Value = ~10g~%"
 381                       (send s :point-value (send s :best-point)))))))
 382
 383
 384 ;;;
 385 ;;; Simplex Prototype
 386 ;;;
 387
 388 (defproto simplex-proto '(f simplex))
 389
 390 ;;;
 391 ;;; Simplex Points
 392 ;;;
 393
 394 (defmeth simplex-proto :make-point (x)
 395   (let ((f (send self :f)))
 396     (if f
 397         (let ((val (funcall f x)))
 398           (cons (if (consp val) (car val) val) x))
 399         (cons nil x))))
 400
 401 (defmeth simplex-proto :point-value (x) (car x))
 402
 403 (defmeth simplex-proto :point-location (x) (cdr x))
 404
 405 (defmeth simplex-proto :is-worse (x y)
 406   (< (send self :point-value x) (send self :point-value y)))
 407
 408 ;;;
 409 ;;; Making New Simplices
 410 ;;;
 411
 412 (defmeth simplex-proto :isnew (f start &optional size)
 413   (send self :simplex start size)
 414   (send self :f f))
 415
 416 ;;;
 417 ;;; Slot Accessors and Mutators
 418 ;;;
 419
 420 (defmeth simplex-proto :simplex (&optional new size)
 421   (if new
 422       (let ((simplex
 423              (if (and (consp new) (sequencep (car new)))
 424                  (if (/= (length new) (+ 1 (length (car new))))
 425                      (error "bad simplex data")
 426                      (copy-list new))
 427                  (let* ((n (length new))
 428                         (size (if size size (repeat 1 n)))
 429                      ;   (pts (- (* 2 (uniform-rand (repeat n (+ n 1)))) 1)))
 430                         (diag (* 2 size (- (random (repeat 2 n)) .5)))
 431                         (pts (cons (repeat 0 n)
 432                                    (mapcar #'(lambda (x) (coerce x 'list))
 433                                            (column-list (diagonal diag))))))
 434                    (mapcar #'(lambda (x) (reduce #'+ (list (* size x) new))) pts)))))
 435         (setf (slot-value 'simplex)
 436               (mapcar #'(lambda (x) (send self :make-point x)) simplex))
 437         (send self :sort-simplex)))
 438   (slot-value 'simplex))
 439
 440 (defmeth simplex-proto :f (&optional f)
 441   (when f
 442         (setf (slot-value 'f) f)
 443         (let ((simplex
 444                (mapcar #'(lambda (x) (send self :point-location x))
 445                        (send self :simplex))))
 446           (send self :simplex simplex)))
 447   (slot-value 'f))
 448
 449 (defmeth simplex-proto :sort-simplex ()
 450   (if (send self :f)
 451       (setf (slot-value 'simplex)
 452             (sort (slot-value 'simplex)
 453                   #'(lambda (x y) (send self :is-worse x y))))))
 454
 455 ;;;
 456 ;;; Other Methods Using List Representation of SImplex
 457 ;;;
 458
 459 (defmeth simplex-proto :best-point () (car (last (send self :simplex))))
 460 (defmeth simplex-proto :worst-point () (first (send self :simplex)))
 461 (defmeth simplex-proto :second-worst-point () (second (send self :simplex)))
 462 (defmeth simplex-proto :replace-point (new old)
 463   (let* ((simplex (send self :simplex))
 464          (n (position old simplex)))
 465     (when n
 466           (setf (nth n simplex) new)
 467           (send self :sort-simplex))))
 468 (defmeth simplex-proto :mean-opposite-face (x)
 469   (let ((face (mapcar #'(lambda (x) (send self :point-location x))
 470                       (remove x (send self :simplex)))))
 471     (/ (reduce #'+ face) (length face))))
 472
 473 ;;;
 474 ;;; Iteration Step Methods
 475 ;;;
 476
 477 (defmeth simplex-proto :extrapolate-from-worst (fac)
 478   (let* ((worst (send self :worst-point))
 479          (wloc (send self :point-location worst))
 480          (delta (- (send self :mean-opposite-face worst) wloc))
 481          (new (send self :make-point (+ wloc (* (- 1 fac) delta)))))
 482     (if (send self :is-worse worst new) (send self :replace-point new worst))
 483     new))
 484
 485 (defmeth simplex-proto :shrink-to-best (fac)
 486   (let* ((best (send self :best-point))
 487          (bloc (send self :point-location best)))
 488     (dolist (x (copy-list (send self :simplex)))
 489             (if (not (eq x best))
 490                 (send self :replace-point
 491                       (send self :make-point
 492                             (+ bloc
 493                                (* fac
 494                                   (- (send self :point-location x) bloc))))
 495                       x)))))
 496
 497 (defmeth simplex-proto :relative-range ()
 498   (let ((best (send self :point-value (send self :best-point)))
 499         (worst (send self :point-value (send self :worst-point))))
 500     (* 2 (/ (abs (- best worst)) (+ 1 (abs best) (abs worst))))))
 501
 502
 503
 504
 505 ;;;;
 506 ;;;; Maximization and Numerical Derivatives
 507 ;;;;
 508
 509
 510 (defun data2double (n data ptr)
 511   (declare (fixnum n))
 512   (let* ((seq (compound-data-seq data))
 513          (elem (make-next-element seq)))
 514     (if (/= (length seq) n) (error "bad data size"))
 515     (dotimes (i n)
 516       (declare (fixnum i))
 517       (la-put-double ptr i (get-next-element elem i)))))
 518
 519 (defun maximize-callback (n px pfval pgrad phess pderivs)
 520   (la-vector-to-data px n mode-re *maximize-callback-arg*)
 521   (let* ((val (funcall *maximize-callback-function* *maximize-callback-arg*))
 522          (derivs (if (consp val) (- (length val) 1) 0)))
 523     (la-put-integer pderivs 0 derivs)
 524     (la-put-double pfval 0 (if (consp val) (first val) val))
 525     (if (<= 1 derivs) (data2double n (second val) pgrad))
 526     (if (<= 2 derivs) (data2double (* n n) (third val) phess))))
 527
 528 (defun numgrad (f x &optional scale (h -1.0))
 529 "Args: (f x &optional scale derivstep)
 530 Computes the numerical gradient of F at X."
 531   (check-sequence x)
 532   (check-real x)
 533   (when scale
 534         (check-sequence scale)
 535         (check-real scale))
 536   (check-one-real h)
 537   (let* ((n (length x))
 538          (result (make-list n)))
 539     (if (and scale (/= n (length scale)))
 540         (error "scale not the same length as x"))
 541     (let ((*maximize-callback-function* f)
 542           (*maximize-callback-arg* (make-list n)))
 543       (let ((px (la-data-to-vector x mode-re))
 544             (pgrad (la-vector n mode-re))
 545             (pscale (la-data-to-vector
 546                      (if scale scale (make-list n :initial-element 1.0))
 547                      mode-re)))
 548         (unwind-protect
 549             (progn
 550               (numgrad-front n px pgrad h pscale)
 551               (la-vector-to-data pgrad n mode-re result))
 552           (la-free-vector px)
 553           (la-free-vector pgrad)
 554           (la-free-vector pscale))))
 555       result))
 556
 557 (defun numhess (f x &optional scale (h -1.0) all)
 558 "Args: (f x &optional scale derivstep)
 559 Computes the numerical Hessian matrix of F at X."
 560   (check-sequence x)
 561   (check-real x)
 562   (when scale
 563         (check-sequence scale)
 564         (check-real scale))
 565   (check-one-real h)
 566   (let* ((n (length x))
 567          (result (if all
 568                      (list nil (make-list n) (make-array (list n n)))
 569                      (make-array (list n n)))))
 570     (if (and scale (/= n (length scale)))
 571         (error "scale not the same length as x"))
 572     (let ((*maximize-callback-function* f)
 573           (*maximize-callback-arg* (make-list n)))
 574       (let ((hess-data (compound-data-seq (if all (third result) result)))
 575             (px (la-data-to-vector x mode-re))
 576             (pf (la-vector 1 mode-re))
 577             (pgrad (la-vector n mode-re))
 578             (phess (la-vector (* n n) mode-re))
 579             (pscale (la-data-to-vector
 580                      (if scale scale (make-list n :initial-element 1.0))
 581                      mode-re)))
 582         (unwind-protect
 583             (progn
 584               (numhess-front n px pf pgrad phess h pscale)
 585               (when all
 586                     (setf (first result) (la-get-double pf 0))
 587                     (la-vector-to-data pgrad n mode-re (second result)))
 588               (la-vector-to-data phess (* n n) mode-re hess-data))
 589           (la-free-vector pf)
 590           (la-free-vector px)
 591           (la-free-vector pgrad)
 592           (la-free-vector phess)
 593           (la-free-vector pscale))))
 594     result))
 595