share/linearalgebra/linalg-utilities.lisp

   1 ;;  Copyright 2005 by Barton Willis
   2
   3 ;;  This is free software; you can redistribute it and/or
   4 ;;  modify it under the terms of the GNU General Public License,
   5 ;;  http://www.gnu.org/copyleft/gpl.html.
   6
   7 ;; This software has NO WARRANTY, not even the implied warranty of
   8 ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
   9
  10 ($put '$linalgutilities 2 '$version)
  11
  12 (defun inform (level pck msg &rest arg)
  13   (if (member level (member ($get '$infolevel pck) `($debug $verbose $silent)))
  14       (apply 'mtell `(,msg ,@arg))))
  15
  16 (defun $require_nonempty_matrix (m pos fun)
  17   (if (not (and ($matrixp m) (> ($length m) 0) (> ($length ($first m)) 0)))
  18       (merror "The ~:M argument of the function ~:M must be a nonempty matrix" pos fun)))
  19
  20 ;; Why both some and every? Because we want blockmatrixp(matrix()) --> false.
  21
  22 (defun $blockmatrixp (m)
  23   (and ($matrixp m) ($some '$matrixp m) ($every '$matrixp m)))
  24
  25 (defun $require_matrix (m pos fun)
  26   (if (not ($matrixp m))
  27       (merror "The ~:M argument of the function ~:M must be a matrix" pos fun)))
  28
  29 (defun $require_unblockedmatrix (m pos fun)
  30   (if (or (not ($matrixp m)) ($blockmatrixp m))
  31       (merror "The ~:M argument of the function ~:M must be an unblocked matrix" pos fun)))
  32
  33 (defun $require_square_matrix (m pos fun)
  34   (if (not (and ($matrixp m) (= ($length m) ($length ($first m)))))
  35       (merror "The ~:M argument of the function ~:M must be a square matrix" pos fun)))
  36
  37 (defun array-elem (m i j)
  38   (nth j (nth i m)))
  39
  40 (defun $require_symmetric_matrix (m pos fun)
  41   (if (not ($matrixp m)) (merror "The ~:M argument to ~:M must be a matrix" pos fun))
  42   (let ((n ($matrix_size m)))
  43     (if (not (= ($first n) ($second n)))
  44         (merror "The ~:M argument to ~:M must be a square matrix" pos fun))
  45     (if (not ($zeromatrixp (sub m ($transpose m))))
  46         (merror "The ~:M argument to ~:M must be a symmetric matrix" pos fun)))
  47   '$done)
  48
  49 (defun $require_real_symmetric_matrix (m pos fun)
  50   (if (not ($matrixp m)) (merror "The ~:M argument to ~:M must be a matrix" pos fun))
  51   (let ((n ($matrix_size m)))
  52     (if (not (= ($first n) ($second n)))
  53         (merror "The ~:M argument to ~:M must be a square matrix" pos fun))
  54     (if (and ($zeromatrixp (sub m ($transpose m))) ($zeromatrixp (sub m (take '($conjugate) m)))) '$done
  55       (merror "The ~:M argument to ~:M must be a real symmetric matrix" pos fun))))
  56
  57 (defun $require_selfadjoint_matrix (m fun pos)
  58   (if (not ($matrixp m)) (merror "The ~:M argument to ~:M must be a matrix" pos fun))
  59   (let ((n ($matrix_size m)))
  60     (if (not (= ($first n) ($second n)))
  61         (merror "The ~:M argument to ~:M must be a square matrix" pos fun))
  62     (if (not ($zeromatrixp (sub m ($ctranspose m))))
  63         (merror "The ~:M argument to ~:M must be a selfadjoint (hermitian) matrix" pos fun)))
  64   '$done)
  65
  66 ;; matrix() is a 0 x 0 matrix, and matrix([]) is a 1 x 0 matrix.
  67 ;; There is no representation for a 0 x 1 matrix. Currently,
  68 ;; transpose(matrix([])) => matrix(). And that's a bug.
  69
  70 (defun $matrix_size(m)
  71   ($require_matrix m '$first '$matrix_size)
  72   `((mlist) ,($length ($args m)) ,(if ($emptyp ($args m)) 0 ($length ($first ($args m))))))
  73
  74 (defun $require_list (lst pos fun)
  75   (if (not ($listp lst))
  76       (merror "The ~:M argument of the function ~:M must be a list" pos fun)))
  77
  78 (defun $require_posinteger(i pos fun)
  79   (if (not (and (integerp i) (> i 0)))
  80       (merror "The ~:M argument of the function ~:M must be a positive integer" pos fun)))
  81
  82 (defun matrix-map (f mat)
  83   (setq mat (mapcar 'cdr (cdr mat)))
  84   (setq mat (mapcar #'(lambda (s) (mapcar f s)) mat))
  85   (setq mat (mapcar #'(lambda (s) (cons `(mlist simp) s)) mat))
  86   `(($matrix simp) ,@mat))
  87
  88 ;; Map the lisp function fn over the matrix m. This function is block matrix friendly.
  89
  90 (defun full-matrix-map (fn m)
  91   (if (or ($listp m) ($matrixp m))
  92       (cons (car m) (mapcar #'(lambda (s) (full-matrix-map fn s)) (cdr m)))
  93     (funcall fn m)))
  94
  95 (defun $zerofor (mat &optional (fld-name '$generalring))
  96   (let* ((fld ($require_ring fld-name '$second '$zerofor))
  97         (add-id (funcall (mring-mring-to-maxima fld) (funcall (mring-add-id fld)))))
  98     (zerofor mat add-id)))
  99
 100 (defun zerofor (mat zero)
 101   (if (or ($matrixp mat) ($listp mat))
 102       (cons (car mat) (mapcar #'(lambda (s) (zerofor s zero)) (cdr mat)))
 103     zero))
 104
 105 ;; Return an identity matrix that has the same shape as the matrix
 106 ;; mat. The first argument 'mat' should be a square Maxima matrix or a
 107 ;; non-matrix. When 'mat' is a matrix, each entry of 'mat' can be a
 108 ;; square matrix -- thus 'mat' can be a blocked Maxima matrix. The
 109 ;; matrix can be blocked to any (finite) depth.
 110
 111 (defun $identfor (mat &optional (fld-name '$generalring))
 112   (let* ((fld ($require_ring fld-name '$second '$zerofor))
 113          (add-id (funcall (mring-mring-to-maxima fld) (funcall (mring-add-id fld))))
 114          (mult-id (funcall (mring-mring-to-maxima fld) (funcall (mring-mult-id fld)))))
 115     (if ($matrixp mat) (identfor mat add-id mult-id) mult-id)))
 116
 117 (defun identfor (mat zero one)
 118   (let ((i) (acc) (j) (new-mat))
 119     (setq mat (rest mat))
 120     (setq i 0)
 121     (dolist (row mat)
 122       (setq row (rest row))
 123       (setq acc nil)
 124       (setq j 0)
 125       (dolist (aij row)
 126         (push (cond (($matrixp aij)
 127                      (if (= i j) (identfor aij zero one) (zerofor aij zero)))
 128                     ((= i j) one)
 129                     (t zero)) acc)
 130         (incf j))
 131       (incf i)
 132       (push '(mlist) acc)
 133       (push acc new-mat))
 134     (push '($matrix) new-mat)))
 135
 136 (defun $ctranspose (m)
 137   (mfuncall '$transpose (full-matrix-map #'(lambda (s) (simplifya `(($conjugate) ,s) nil)) m)))
 138
 139 (defun $zeromatrixp (m)
 140   (if (or ($matrixp m) ($listp m)) (every '$zeromatrixp (cdr m))
 141     (eq '$zero (csign ($rectform m)))))
 142
 143 (defun array-to-row-list (mat &optional (fn 'identity))
 144   (let ((acc) (r (array-dimensions mat)) (row) (c))
 145     (setq c (second r))
 146     (setq r (first r))
 147     (dotimes (i r)
 148       (setq row nil)
 149       (dotimes (j c)
 150         (push (funcall fn (aref mat i j)) row))
 151       (setq row (reverse row))
 152       (push row acc))
 153     (reverse acc)))
 154
 155 (defun array-to-maxima-matrix (mat &optional (fn 'identity))
 156   (cons '($matrix) (mapcar #'(lambda (s) (cons '(mlist) s)) (array-to-row-list mat fn))))
 157
 158 (defun array-to-maxima-list (ar &optional (fn 'identity))
 159   (cons '(mlist) (mapcar fn (coerce ar 'list))))
 160
 161 (defun maxima-to-array (mat &optional (fn 'identity) typ)
 162   (let ((r ($matrix_size mat)) (c))
 163     (setq c ($second r))
 164     (setq r ($first r))
 165     (setq mat (mapcar #'cdr (cdr mat)))
 166     (setq mat (mapcar #'(lambda (s) (mapcar #'(lambda (w) (funcall fn w)) s)) mat))
 167     (if typ (make-array (list r c) :element-type typ :initial-contents mat)
 168       (make-array (list r c) :initial-contents mat))))