Develop

README

@@ -1,6 +1,50 @@
+Linear Algebra in Common Lisp
+
+Copyright (c) 2011-2014, Odonata Research LLC
+
+Permission is hereby granted, free  of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction,  including without limitation the rights
+to use, copy, modify,  merge,  publish,  distribute,  sublicense, and/or sell
+copies of the  Software,  and  to  permit  persons  to  whom  the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and  this  permission  notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED  "AS IS",  WITHOUT  WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT  NOT  LIMITED  TO  THE  WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE  AND  NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT  HOLDERS  BE  LIABLE  FOR  ANY  CLAIM,  DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
+Copyright (c) 2021-2023, Ten Factor Growth, LLC
+Permission is hereby granted, free  of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction,  including without limitation the rights
+to use, copy, modify,  merge,  publish,  distribute,  sublicense, and/or sell
+copies of the  Software,  and  to  permit  persons  to  whom  the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and  this  permission  notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED  "AS IS",  WITHOUT  WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT  NOT  LIMITED  TO  THE  WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE  AND  NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT  HOLDERS  BE  LIABLE  FOR  ANY  CLAIM,  DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
+
 linear-algebra
 
 linear-algebra is a Common Lisp library of numeric linear algebra routines.
 It approximates the functionality present in the Basic Linear Algebra
 Subroutines (BLAS) Fortran library. The objective of this project is to
 generate linear algebra routines in native ANSI Common Lisp.
+
+Use within quicklisp with (ql:quickload :linear-algebra)

lisp/data-vector.lisp

@@ -36,6 +36,10 @@
   (:documentation
    "A data vector."))
 
+(defmethod print-object ((vector data-vector) stream)
+  (print-unreadable-object (vector stream :type t :identity t)
+    (format stream "~a" (contents vector))))
+
 (defclass row-vector (data-vector)
   ()
   (:documentation
@@ -200,6 +204,9 @@ applying the function to each element of the vectors."
          function
          vectors))
 
+(defmethod reduce-vector ((function function) (vector data-vector) initial-value)
+  (reduce function (contents vector) :initial-value initial-value))
+
 ;;; Data vector transformations
 
 (defmethod apply-rotation :before
@@ -260,6 +267,10 @@ applying the function to each element of the vectors."
   "Return the p-norm of the vector."
   (norm-vector (contents vector) measure))
 
+(defmethod min-vector ((vector data-vector))
+  (loop for element across (contents vector) minimize (abs element)))
+
+
 (defmethod transpose ((vector column-vector))
   "Return a row vector."
   (make-instance
@@ -434,3 +445,101 @@ applying the function to each element of the vectors."
     ((vector1 row-vector) (vector2 column-vector) &optional scalar)
   "Return the dot product of vector1 and vector2."
   (inner-product-vector (contents vector1) (contents vector2) scalar))
+
+(defmethod product :before
+  ((vector1 data-vector) (vector2 data-vector) &optional scalar)
+  "Verify that the dimensions are equal."
+  (declare (ignore scalar))
+  (unless (= (vector-length vector1) (vector-length vector2))
+    (error "VECTOR1 and VECTOR2 are not of equal length.")))
+
+(defmethod product
+    ((vector1 data-vector) (vector2 data-vector) &optional scalar)
+  "Return the dot product of vector1 and vector2."
+  (inner-product-vector (contents vector1) (contents vector2) scalar))
+
+(defmethod vec-equal :before
+  ((vector1 data-vector) (vector2 data-vector))
+  "Verify that the dimensions are equal."
+  (declare (ignore scalar))
+  (unless (= (vector-length vector1) (vector-length vector2))
+    (error "VECTOR1 and VECTOR2 are not of equal length.")))
+
+(defmethod vec-equal
+    ((vector1 data-vector) (vector2 data-vector))
+  (equalp (contents vector1) (contents vector2)))
+
+
+(defmethod elem-divide :before
+  ((vector1 data-vector) (vector2 data-vector))
+  "Verify that the dimensions are equal."
+  (declare (ignore scalar))
+  (unless (= (vector-length vector1) (vector-length vector2))
+    (error "VECTOR1 and VECTOR2 are not of equal length.")))
+
+(defmethod elem-divide
+    ((vector1 data-vector) (vector2 data-vector))
+  (make-instance
+   (common-class-of vector1 vector2)
+   :contents
+   (element-divide-vector
+    (contents vector1) (contents vector2))
+   ))
+
+(defmethod elem-multiply :before
+  ((vector1 data-vector) (vector2 data-vector))
+  "Verify that the dimensions are equal."
+  (declare (ignore scalar))
+  (unless (= (vector-length vector1) (vector-length vector2))
+    (error "VECTOR1 and VECTOR2 are not of equal length.")))
+
+(defmethod elem-multiply
+    ((vector1 data-vector) (vector2 data-vector))
+  (make-instance
+   (common-class-of vector1 vector2)
+   :contents
+   (element-multiply-vector
+    (contents vector1) (contents vector2))
+   ))
+
+(defmethod elem-greater
+    ((vector1 data-vector) (vector2 data-vector))
+  (element-greater-vector
+   (contents vector1) (contents vector2)))
+
+(defmethod vec-every ((vector data-vector) predicate)
+  (every predicate (contents vector)))
+
+
+(defmethod distance :before
+  ((vector1 data-vector) (vector2 data-vector) &optional (measure 2))
+  "Verify that the dimensions are equal."
+  (declare (ignore scalar))
+  (unless (= (vector-length vector1) (vector-length vector2))
+    (error "VECTOR1 and VECTOR2 are not of equal length.")))
+
+(defmethod distance
+    ((vector1 data-vector) (vector2 data-vector) &optional (measure 2))
+  (norm (subtract vector1 vector2) measure))
+
+(defmethod outer-product ((vector1 data-vector)
+                          (vector2 data-vector))
+  (let* ((len1 (vector-length vector1))
+         (len2 (vector-length vector2))
+         (mat-type (cond
+                     ((equal vector1 vector2) 'symmetric-matrix)
+                     ((= len1 len2) 'square-matrix)
+                     (t 'dense-matrix)))
+         (mat (make-instance mat-type
+                             :dimensions (list len1 len2)
+                             :initial-element 0
+                             :element-type (vector-element-type vector1)))
+         )
+    (dotimes (i len1)
+      (dotimes (j len2)
+        (setf
+         (mref mat i j)
+         (* (vref vector1 i) (vref vector2 j)))
+        ))
+    mat
+  ))

lisp/dense-matrix.lisp

@@ -34,6 +34,10 @@
   (:documentation
    "Dense matrix object."))
 
+(defmethod print-object ((matrix dense-matrix) stream)
+  (print-unreadable-object (matrix stream :type t :identity t)
+    (format stream "~a" (contents matrix))))
+
 (defmethod initialize-instance :after
   ((self dense-matrix) &rest initargs
    &key dimensions element-type initial-element initial-contents)
@@ -110,6 +114,10 @@
   "Return true if object is a dense matrix."
   (typep object 'dense-matrix))
 
+(defmethod mat-equal
+    ((matrix1 dense-matrix) (matrix2 dense-matrix))
+  (equalp (contents matrix1) (contents matrix2)))
+
 (defmethod matrix-in-bounds-p
     ((matrix dense-matrix) (row integer) (column integer))
   "Return true if row and column do not exceed the dimensions of matrix."
@@ -490,3 +498,44 @@ matrix with a column vector."
        :contents
        (gauss-invert (contents matrix)))
       (error "The number of rows does not equal columns.")))
+
+(defmethod add-diagonal ((value number) (matrix dense-matrix))
+  (let ((result (copy-matrix matrix))
+	(num    (apply #'min (matrix-dimensions matrix)))
+	)
+    (dotimes (index num result)
+      (setf
+       (mref result index index)
+       (+ (mref matrix index index) value)))
+    result))
+
+(defmethod add-diagonal :before
+    ((vec data-vector) (matrix dense-matrix))
+  (unless (compatible-dimensions-p :solve matrix vector)
+    (error "Matrix~A is incompatible with column vector(~D)."
+           (matrix-dimensions matrix) (vector-length vector))))
+
+(defmethod add-diagonal ((vec data-vector) (matrix dense-matrix))
+  (let ((result (copy-matrix matrix))
+	(num    (min (matrix-dimensions matrix)))
+	)
+    (dotimes (index num result)
+      (setf
+       (mref result index index)
+       (+ (mref matrix index index)
+	  (vref vec index))))
+    result))
+
+(defmethod matrix-diagonal ((matrix dense-matrix))
+  (let* ((num (apply #'min (matrix-dimensions matrix)))
+         (vec (make-instance 'data-vector :size num
+                                          :initial-element 0
+                                          :element-type (matrix-element-type matrix)))
+         )
+    (dotimes (index num )
+      (setf
+       (vref vec    index)
+       (mref matrix index index))
+      )
+    vec
+  ))

lisp/interface/fundamental-ops.lisp

@@ -32,6 +32,14 @@
   (:documentation
    "Return the norm according to measure."))
 
+(defgeneric min-vector (vector)
+  (:documentation
+   "Return the minimum absolute value"))
+
+(defgeneric distance (vector1 vector2 &optional measure)
+  (:documentation
+   "Return the distance between vector1 and vector2 using the measure"))
+
 (defgeneric transpose (vector-or-matrix)
   (:documentation
    "Transpose the vector or matrix."))
@@ -92,3 +100,8 @@
 (defgeneric ninvert (matrix)
   (:documentation
    "Return the invert of the matrix with in-place decomposition."))
+
+(defgeneric add-diagonal (scalar-or-vector matrix)
+  (:documentation
+   "Add a scalar or vector to every diagonal of the matrix"))
+

lisp/interface/matrix.lisp

@@ -99,6 +99,10 @@ MATRIX."))
   (:documentation
    "Set the submatrix of the matrix."))
 
+(defgeneric matrix-trace (matrix)
+  (:documentation
+   "Return a column vector that is the trace of matrix"))
+
 (defgeneric replace-matrix
     (matrix1 matrix2 &key
      start-row1 end-row1
@@ -120,3 +124,11 @@ MATRIX."))
         (values start-row start-column end-row end-column)
         (error "The matrix range (~D:~D,~D:~D) is invalid."
                start-row start-column end-row end-column))))
+
+(defgeneric mat-equal (matrix1 matrix2)
+  (:documentation
+   "Return if the two arrays elments are the same"))
+
+(defgeneric matrix-diagonal (matrix)
+  (:documentation
+   "Return the diagonal of the matrix as a vector"))

lisp/interface/vector.lisp

@@ -89,6 +89,10 @@
    "Destructively modifies the result vector with the result of
 applying the function to each element of the vectors."))
 
+(defgeneric reduce-vector (function vector initial-value)
+  (:documentation
+   "Reduce a vector into a scalar using a function"))
+
 (defmacro dovector ((element vector &optional result) &body body)
   "Iterate over vector returning result."
   (let ((pos (gensym "POS-"))
@@ -108,3 +112,28 @@ applying the function to each element of the vectors."))
 (defgeneric napply-rotation (vector1 vector2 cc ss)
   (:documentation
    "Return the plane rotations of vector1 and vector2 by cc and ss."))
+
+(defgeneric vec-equal (vector1 vector2)
+  (:documentation
+   "Test if the contents of vector1 and vector2 are the same"))
+
+(defgeneric elem-divide (vector1 vector2)
+  (:documentation
+   "Return the element by element division of two vectors"))
+
+(defgeneric elem-multiply (vector1 vector2)
+  (:documentation
+   "Return the element by element multiplation of two vectors"))
+
+(defgeneric elem-greater (vector1 vector2)
+  (:documentation
+   "Return if vector1 is greater than vector 2 on every element"))
+
+(defgeneric vec-every (vector predicate)
+  (:documentation
+   "Return true if every element in the vector satisfies
+  predicate. Return nil on the first element that does not"))
+
+(defgeneric outer-product (vec1 vec2)
+  (:documentation
+   "Outer product of vec1 and vec2 returns an nxm matrix"))

lisp/kernel/binary-operations.lisp

@@ -75,6 +75,16 @@ the operation."))
   "Return the scaled operation."
   (lambda (n1 n2) (- (* scalar1 n1) (* scalar2 n2))))
 
+(defmethod binary-op
+    ((op (eql #'/)) )
+  "Return the scaled operation."
+  (lambda (n1 n2) (/ n1 n2)))
+
+(defmethod binary-op
+    ((op (eql #'*)) )
+  "Return the scaled operation."
+  (lambda (n1 n2) (* n1 n2)))
+
 ;;; Binary vector operations
 
 (defun %vector<-vector1-op-vector2 (operation vector1 vector2)
@@ -124,6 +134,27 @@ addition."
    (scaled-binary-op #'- scalar1 scalar2)
    vector1 vector2))
 
+(defun element-divide-vector (vector1 vector2)
+  "Vector pointwise division."
+  (%vector<-vector1-op-vector2
+   (binary-op #'/)
+   vector1 vector2))
+
+(defun element-multiply-vector (vector1 vector2)
+  "Vector pointwise multiplication."
+  (%vector<-vector1-op-vector2
+   (binary-op #'*)
+   vector1 vector2))
+
+(defun element-greater-vector (vector1 vector2)
+  (loop
+   for element1 across vector1
+   and element2 across vector2
+   when (<= element1 element2)
+     return nil
+   finally (return t)))
+
+
 (defun nsubtract-vector (vector1 vector2 scalar1 scalar2)
   "Destructive vector binary subtraction."
   (%vector1<-vector1-op-vector2
@@ -137,7 +168,7 @@ addition."
    and element2 across vector2
    sum (* element1 element2) into result
    finally
-   (return (if scalar (* scalar result) result))))
+      (return (if scalar (* scalar result) result))))
 
 ;;; Binary array/vector operations