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