share/minpack/minpack-interface.lisp

   1 ;;; -*- Mode: lisp -*-
   2
   3 ;;; Simple Maxima interface to minpack routines
   4
   5 (in-package #:maxima)
   6
   7 (defmvar $debug_minpack nil
   8   "Set to true to enable debugging messages from minpack routines")
   9
  10 (defmfun $minpack_lsquares (fcns vars init-x
  11                                  &key
  12                                  (jacobian t)
  13                                  (tolerance #.(sqrt double-float-epsilon)))
  14   "Minimize the sum of the squares of m functions in n unknowns (n <= m)
  15
  16    VARS    list of the variables
  17    INIT-X  initial guess
  18    FCNS    list of the m functions
  19
  20    Optional keyword args (key = val)
  21    TOLERANCE  tolerance in solution
  22    JACOBIAN   If true, maxima computes the Jacobian directly from FCNS.
  23               If false, the Jacobian is internally computed using a
  24               forward-difference approximation.
  25               Otherwise, it is a function returning the Jacobian"
  26   (unless (and (listp fcns) (eq (caar fcns) 'mlist))
  27     (merror "~M: ~M is not a list of functions" %%pretty-fname fcns))
  28   (unless (and (listp vars) (eq (caar vars) 'mlist))
  29     (merror "~M: ~M is not a list of variables" %%pretty-fname vars))
  30   (unless (and (listp init-x) (eq (caar init-x) 'mlist))
  31     (merror "~M: ~M is not a list of initial values" %%pretty-fname init-x))
  32   (setf tolerance ($float tolerance))
  33   (unless (and (realp tolerance) (plusp tolerance))
  34     (merror "~M: tolerance must be a non-negative real number, not: ~M"
  35             %%pretty-fname tolerance))
  36
  37   (let* ((n (length (cdr vars)))
  38          (m (length (cdr fcns)))
  39          (x (make-array n :element-type 'double-float
  40                         :initial-contents (mapcar #'(lambda (z)
  41                                                       ($float z))
  42                                                   (cdr init-x))))
  43          (fvec (make-array m :element-type 'double-float))
  44          (fjac (make-array (* m n) :element-type 'double-float))
  45          (ldfjac m)
  46          (info 0)
  47          (ipvt (make-array n :element-type 'f2cl-lib:integer4))
  48          (fv (coerce-float-fun fcns vars))
  49          (fj (cond ((eq jacobian t)
  50                     ;; T means compute it ourselves
  51                     (meval `(($jacobian) ,fcns ,vars)))
  52                    (jacobian
  53                     ;; Use the specified Jacobian
  54                     )
  55                    (t
  56                     ;; No jacobian at all
  57                     nil))))
  58     ;; Massage the Jacobian into a function
  59     (when jacobian
  60       (setf fj (coerce-float-fun fj vars)))
  61     (cond
  62       (jacobian
  63        ;; Jacobian given (or computed by maxima), so use lmder1
  64        (let* ((lwa (+ m (* 5 n)))
  65               (wa (make-array lwa :element-type 'double-float)))
  66          (labels ((fcn-and-jacobian (m n x fvec fjac ldfjac iflag)
  67                     (declare (type f2cl-lib:integer4 m n ldfjac iflag)
  68                              (type (cl:array double-float (*)) x fvec fjac))
  69                     (ecase iflag
  70                       (1
  71                        ;; Compute function at point x, placing result
  72                        ;; in fvec (subseq needed because sometimes
  73                        ;; we're called with vector that is longer than
  74                        ;; we want.  Perfectly valid Fortran, though.)
  75                        (let ((val (apply 'funcall fv (subseq (coerce x 'list) 0 n))))
  76                          (unless (consp val)
  77                            (merror "Unable to evaluate function at the point ~M"
  78                                    (list* '(mlist) (subseq (coerce x 'list) 0 n))))
  79                          (when $debug_minpack
  80                            (format t "f(~{~A~^, ~}) =~%[~@<~{~A~^, ~:_~}~:>]~%"
  81                                    (coerce x 'list)
  82                                    (cdr val)))
  83                          (replace fvec (mapcar #'(lambda (z)
  84                                                    (cl:float z 1d0))
  85                                                (cdr val)))))
  86                       (2
  87                        ;; Compute Jacobian at point x, placing result in fjac
  88                        (let ((j (apply 'funcall fj (subseq (coerce x 'list) 0 n))))
  89                          (unless (consp j)
  90                            (merror "Unable to evaluate Jacobian at the point ~M"
  91                                    (list* '(mlist) (subseq (coerce x 'list) 0 n))))
  92                          ;; Extract out elements of Jacobian and place into
  93                          ;; fjac, in column-major order.
  94                          (let ((row-index 0))
  95                            (dolist (row (cdr j))
  96                              (let ((col-index 0))
  97                                (dolist (col (cdr row))
  98                                  (setf (aref fjac (+ row-index (* ldfjac col-index)))
  99                                        (cl:float col 1d0))
 100                                  (incf col-index)))
 101                              (incf row-index))))))
 102                     (values m n nil nil nil ldfjac iflag)))
 103            (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5 ldfjac var-6 info)
 104                (minpack:lmder1 #'fcn-and-jacobian
 105                                m
 106                                n
 107                                x
 108                                fvec
 109                                fjac
 110                                ldfjac
 111                                tolerance
 112                                info
 113                                ipvt
 114                                wa
 115                                lwa)
 116              (declare (ignore ldfjac var-0 var-1 var-2 var-3 var-4 var-5 var-6))
 117
 118              ;; Return a list of the solution and the info flag
 119              (list '(mlist)
 120                    (list* '(mlist) (coerce x 'list))
 121                    (minpack:enorm m fvec)
 122                    info)))))
 123       (t
 124        ;; No Jacobian given so we need to use differences to compute
 125        ;; a numerical Jacobian.  Use lmdif1.
 126        (let* ((lwa (+ m (* 5 n) (* m n)))
 127               (wa (make-array lwa :element-type 'double-float)))
 128          (labels ((fval (m n x fvec iflag)
 129                     (declare (type f2cl-lib:integer4 m n ldfjac iflag)
 130                              (type (cl:array double-float (*)) x fvec fjac))
 131                     ;; Compute function at point x, placing result in fvec
 132                     (let ((val (apply 'funcall fv (subseq (coerce x 'list) 0 n))))
 133                       (unless (consp val)
 134                         (merror "Unable to evaluate function at the point ~M"
 135                                 (list* '(mlist) (subseq (coerce x 'list) 0 n))))
 136                       (when $debug_minpack
 137                         (format t "f(~{~A~^, ~}) =~%[~@<~{~A~^, ~:_~}~:>]~%"
 138                                 (coerce x 'list)
 139                                 (cdr val)))
 140                       (replace fvec (mapcar #'(lambda (z)
 141                                                 (cl:float z 1d0))
 142                                             (cdr val))))
 143                     (values m n nil nil iflag)))
 144            (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5 info)
 145                (minpack:lmdif1 #'fval
 146                                m
 147                                n
 148                                x
 149                                fvec
 150                                tolerance
 151                                info
 152                                ipvt
 153                                wa
 154                                lwa)
 155              (declare (ignore var-0 var-1 var-2 var-3 var-4 var-5))
 156
 157              ;; Return a list of the solution and the info flag
 158              (list '(mlist)
 159                    (list* '(mlist) (coerce x 'list))
 160                    (minpack:enorm m fvec)
 161                    info))))))))
 162
 163 (defmfun $minpack_solve (fcns vars init-x
 164                               &key
 165                               (jacobian t)
 166                               (tolerance #.(sqrt double-float-epsilon)))
 167   "Solve the system of n equations in n unknowns
 168
 169    VARS    list of the n variables
 170    INIT-X  initial guess
 171    FCNS    list of the n functions
 172
 173    Optional keyword args (key = val)
 174    TOLERANCE  tolerance in solution
 175    JACOBIAN   If true, maxima computes the Jacobian directly from FCNS.
 176               If false, the Jacobian is internally computed using a
 177               forward-difference approximation.
 178               Otherwise, it is a function returning the Jacobian"
 179   (unless (and (listp fcns) (eq (caar fcns) 'mlist))
 180     (merror "~M: ~M is not a list of functions"
 181             %%pretty-fname fcns))
 182   (unless (and (listp vars) (eq (caar vars) 'mlist))
 183     (merror "~M:  ~M is not a list of variables"
 184             %%pretty-fname vars))
 185   (unless (and (listp init-x) (eq (caar init-x) 'mlist))
 186     (merror "~M: ~M is not a list of initial values"
 187             %%pretty-fname init-x))
 188   (unless (and (realp tolerance) (plusp tolerance))
 189     (merror "~M: tolerance must be a non-negative real number, not: ~M"
 190             %%pretty-fname tolerance))
 191
 192   (let* ((n (length (cdr vars)))
 193          (x (make-array n :element-type 'double-float
 194                         :initial-contents (mapcar #'(lambda (z)
 195                                                       ($float z))
 196                                                   (cdr init-x))))
 197          (fvec (make-array n :element-type 'double-float))
 198          (fjac (make-array (* n n) :element-type 'double-float))
 199          (ldfjac n)
 200          (info 0)
 201          (fv (coerce-float-fun fcns vars))
 202          (fj (cond ((eq jacobian t)
 203                     ;; T means compute it ourselves
 204                     (mfuncall '$jacobian fcns vars))
 205                    (jacobian
 206                     ;; Use the specified Jacobian
 207                     )
 208                    (t
 209                     ;; No jacobian at all
 210                     nil))))
 211     ;; Massage the Jacobian into a function
 212     (when jacobian
 213       (setf fj (coerce-float-fun fj vars)))
 214     (cond
 215       (jacobian
 216        ;; Jacobian given (or computed by maxima), so use lmder1
 217        (let* ((lwa (/ (* n (+ n 13)) 2))
 218               (wa (make-array lwa :element-type 'double-float)))
 219          (labels ((fcn-and-jacobian (n x fvec fjac ldfjac iflag)
 220                     (declare (type f2cl-lib:integer4 n ldfjac iflag)
 221                              (type (cl:array double-float (*)) x fvec fjac))
 222                     (ecase iflag
 223                       (1
 224                        ;; Compute function at point x, placing result
 225                        ;; in fvec (subseq needed because sometimes
 226                        ;; we're called with vector that is longer than
 227                        ;; we want.  Perfectly valid Fortran, though.)
 228                        (let ((val (apply 'funcall fv (subseq (coerce x 'list) 0 n))))
 229                          (replace fvec (mapcar #'(lambda (z)
 230                                                    (cl:float z 1d0))
 231                                                (cdr val)))))
 232                       (2
 233                        ;; Compute Jacobian at point x, placing result in fjac
 234                        (let ((j (apply 'funcall fj (subseq (coerce x 'list) 0 n))))
 235                          ;; Extract out elements of Jacobian and place into
 236                          ;; fjac, in column-major order.
 237                          (let ((row-index 0))
 238                            (dolist (row (cdr j))
 239                              (let ((col-index 0))
 240                                (dolist (col (cdr row))
 241                                  (setf (aref fjac (+ row-index (* ldfjac col-index)))
 242                                        (cl:float col 1d0))
 243                                  (incf col-index)))
 244                              (incf row-index))))))
 245                     (values n nil nil nil ldfjac iflag)))
 246            (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 ldfjac var-6 info)
 247                (minpack:hybrj1 #'fcn-and-jacobian
 248                                n
 249                                x
 250                                fvec
 251                                fjac
 252                                ldfjac
 253                                tolerance
 254                                info
 255                                wa
 256                                lwa)
 257              (declare (ignore ldfjac var-0 var-1 var-2 var-3 var-4 var-6))
 258              ;; Return a list of the solution and the info flag
 259              (list '(mlist)
 260                    (list* '(mlist) (coerce x 'list))
 261                    (minpack:enorm n fvec)
 262                    info)))))
 263       (t
 264        ;; No Jacobian given so we need to use differences to compute
 265        ;; a numerical Jacobian.  Use lmdif1.
 266        (let* ((lwa (/ (* n (+ (* 3 n) 13)) 2))
 267               (wa (make-array lwa :element-type 'double-float)))
 268          (labels ((fval (n x fvec iflag)
 269                     (declare (type f2cl-lib:integer4 n ldfjac iflag)
 270                              (type (cl:array double-float (*)) x fvec fjac))
 271                     ;; Compute function at point x, placing result in fvec
 272                     (let ((val (apply 'funcall fv (subseq (coerce x 'list) 0 n))))
 273                       (replace fvec (mapcar #'(lambda (z)
 274                                                 (cl:float z 1d0))
 275                                             (cdr val))))
 276                     (values n nil nil iflag)))
 277            (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 info)
 278                (minpack:hybrd1 #'fval
 279                                n
 280                                x
 281                                fvec
 282                                tolerance
 283                                info
 284                                wa
 285                                lwa)
 286              (declare (ignore var-0 var-1 var-2 var-3 var-4))
 287
 288              ;; Return a list of the solution and the info flag
 289              (list '(mlist)
 290                    (list* '(mlist) (coerce x 'list))
 291                    (minpack:enorm n fvec)
 292                    info))))))))