August 2006 - common-math-cvs@common-lisp.net

[common-math-cvs] r3 - trunk/common-math/numerics/linear-algebra
by mantoniotti＠common-lisp.net 18 Aug '06

18 Aug '06

Author: mantoniotti Date: Fri Aug 18 11:54:06 2006 New Revision: 3 Added: trunk/common-math/numerics/linear-algebra/determinants.lisp trunk/common-math/numerics/linear-algebra/inversion.lisp Modified: trunk/common-math/numerics/linear-algebra/linear-algebra-pkg.lisp trunk/common-math/numerics/linear-algebra/linear-algebra.system trunk/common-math/numerics/linear-algebra/lu-decomposition.lisp trunk/common-math/numerics/linear-algebra/matrix.lisp Log: Rearranged a few files, added two more and fixed a small number of issues here and there. Added: trunk/common-math/numerics/linear-algebra/determinants.lisp ============================================================================== --- (empty file) +++ trunk/common-math/numerics/linear-algebra/determinants.lisp Fri Aug 18 11:54:06 2006 @@ -0,0 +1,44 @@ +;;;; -*- Mode: Lisp -*- + +;;;; determinants.lisp -- + + +(in-package "CL.MATH.LINEAR-ALGEBRA") + +;;;=========================================================================== +;;; Determinants + +;;; det + +(defgeneric det (m)) + +(defgeneric det/double-float (m)) + + +(defmethod det ((a array)) + (assert (square-matrix-p a)) + (multiple-value-bind (lu permutations parity) + (lu-decompose a) + (declare (ignore permutations)) + (dotimes (i (array-dimension lu 0) parity) + (setf parity (* parity (aref lu i i)))))) + +(defmethod det ((a matrix)) + (det (matrix-data a))) + + +(defmethod det/double-float ((a array)) + (assert (square-matrix-p a)) + (multiple-value-bind (lu permutations parity) + (lu-decompose a) + (declare (type (array double-float (cl:* cl:*)) lu) + (ignore permutations) + (type double-float parity)) + (dotimes (i (array-dimension lu 0) parity) + (setf parity (* parity (aref lu i i)))))) + +(defmethod det/double-float ((a matrix)) + (det/double-float (matrix-data a))) + + +;;; end of file -- determinant.lisp -- Added: trunk/common-math/numerics/linear-algebra/inversion.lisp ============================================================================== --- (empty file) +++ trunk/common-math/numerics/linear-algebra/inversion.lisp Fri Aug 18 11:54:06 2006 @@ -0,0 +1,29 @@ +;;;; -*- Mode: Lisp -*- + +;;;; inversion.lisp -- + + +(in-package "CL.MATH.LINEAR-ALGEBRA") + +;;;=========================================================================== +;;; Inversion + +;;; inverse + +(defgeneric inverse (m)) + +(defgeneric inverse/double-float (m)) + + + +(defmethod inverse ((a matrix)) + (inverse (matrix-data a))) + + + + +(defmethod inverse/double-float ((a matrix)) + (inverse/double-float (matrix-data a))) + + +;;; end of file -- inversion.lisp -- Modified: trunk/common-math/numerics/linear-algebra/linear-algebra-pkg.lisp ============================================================================== --- trunk/common-math/numerics/linear-algebra/linear-algebra-pkg.lisp (original) +++ trunk/common-math/numerics/linear-algebra/linear-algebra-pkg.lisp Fri Aug 18 11:54:06 2006 @@ -13,6 +13,8 @@ (:export ".*" "./" "EXPT*") (:export ".*%2" ".*%1" "./%2" "./%1") + + (:export "DET" "SOLVE") ) Modified: trunk/common-math/numerics/linear-algebra/linear-algebra.system ============================================================================== --- trunk/common-math/numerics/linear-algebra/linear-algebra.system (original) +++ trunk/common-math/numerics/linear-algebra/linear-algebra.system Fri Aug 18 11:54:06 2006 @@ -8,6 +8,9 @@ "common-math-extra" "vector" "matrix" - "lu-decomposition")) + (:file "lu-decomposition" :depends-on ("matrix")) + (:file "determinants" :depends-on ("lu-decomposition")) + (:file "inversion" :depends-on ("lu-decomposition")) + )) ;;;; end of file -- linear-algebra.system -- Modified: trunk/common-math/numerics/linear-algebra/lu-decomposition.lisp ============================================================================== --- trunk/common-math/numerics/linear-algebra/lu-decomposition.lisp (original) +++ trunk/common-math/numerics/linear-algebra/lu-decomposition.lisp Fri Aug 18 11:54:06 2006 @@ -14,13 +14,6 @@ (defparameter *epsilon* 1.0e-20) ; Change this to CL constant. -(defgeneric square-matrix-p (m) - (:method ((m array)) - (and (matrix-array-p m) - (cl:= (array-dimension m 0) (array-dimension m 1)))) - (:method ((m matrix)) - (cl:= (matrix-row-number m) (matrix-column-number m)))) - ;;; Different base type. @@ -283,33 +276,11 @@ (solve/double-float (matrix-data a) b result)) -;;;--------------------------------------------------------------------------- -;;; det - -(defmethod det ((a array)) - (assert (square-matrix-p a)) - (multiple-value-bind (lu permutations parity) - (lu-decompose a) - (declare (ignore permutations)) - (dotimes (i (array-dimension lu 0) parity) - (setf parity (* parity (aref lu i i)))))) - -(defmethod det ((a matrix)) - (det (matrix-data a))) - - -(defmethod det/double-float ((a array)) - (assert (square-matrix-p a)) - (multiple-value-bind (lu permutations parity) - (lu-decompose a) - (declare (type (array double-float (cl:* cl:*)) lu) - (ignore permutations) - (type double-float parity)) - (dotimes (i (array-dimension lu 0) parity) - (setf parity (* parity (aref lu i i)))))) +;;;=========================================================================== +;;; Determinants -(defmethod det/double-float ((a matrix)) - (det/double-float (matrix-data a))) +;;;=========================================================================== +;;; Inverse ;;;=========================================================================== Modified: trunk/common-math/numerics/linear-algebra/matrix.lisp ============================================================================== --- trunk/common-math/numerics/linear-algebra/matrix.lisp (original) +++ trunk/common-math/numerics/linear-algebra/matrix.lisp Fri Aug 18 11:54:06 2006 @@ -146,6 +146,14 @@ (shape-equal-p m1 m2)) + +(defgeneric square-matrix-p (m) + (:method ((m array)) + (cl:= (matrix-row-number m) (matrix-column-number m))) + (:method ((m matrix)) + (cl:= (matrix-row-number m) (matrix-column-number m)))) + + ;;;--------------------------------------------------------------------------- ;;; Creation.

1 0

[common-math-cvs] r2 - in trunk/common-math: . numerics/linear-algebra
by mantoniotti＠common-lisp.net 17 Aug '06

17 Aug '06

Author: mantoniotti Date: Thu Aug 17 12:39:36 2006 New Revision: 2 Modified: trunk/common-math/common-math.lisp trunk/common-math/numerics/linear-algebra/matrix.lisp Log: Added a few missing operations (.*%2 ./%2 .%\\2) and fixed a couple of loose things. Modified: trunk/common-math/common-math.lisp ============================================================================== --- trunk/common-math/common-math.lisp (original) +++ trunk/common-math/common-math.lisp Thu Aug 17 12:39:36 2006 @@ -52,6 +52,8 @@ ;;; purposes and setting up a separate readtable seems a little bit ;;; too much for the time being. +;;; <%2 + (defgeneric <%2 (x y) (:documentation "The <%2 generic function. The binary LESS generic function which is specialized for various @@ -90,6 +92,12 @@ ) +(defmethod <%2 ((x symbol) (y symbol)) + (string< x y)) + + +;;; >%2 + (defgeneric >%2 (x y) (:method ((x number) (y number)) (cl:> x y)) @@ -122,14 +130,13 @@ ) -(defmethod <%2 ((x symbol) (y symbol)) - (string< x y)) - (defmethod >%2 ((x symbol) (y symbol)) (string> y x)) +;;; <=%2 + (defgeneric <=%2 (x y) (:method ((x number) (y number)) (cl:<= x y)) @@ -162,6 +169,8 @@ ) +;;; >=%2 + (defgeneric >=%2 (x y) (:method ((x number) (y number)) (cl:>= x y)) @@ -194,6 +203,8 @@ ) +;;; +%2 + (defgeneric +%2 (x y &optional r) (:method ((x number) (y number) &optional r) (declare (ignore r)) @@ -237,6 +248,8 @@ ) +;;; *%2 + (defgeneric *%2 (x y &optional r) (:method ((x number) (y number) &optional r) (declare (ignore r)) @@ -289,12 +302,16 @@ ) +;;; -%2 + (defgeneric -%2 (x y &optional r) (:method ((x t) (y t) &optional r) (declare (ignore r)) (+%2 x (-%1 y)))) +;;; /%2 + (defgeneric /%2 (x y &optional r) (:method ((x number) (y number) &optional r) (declare (ignore r)) @@ -305,6 +322,8 @@ (*%2 x (/%1 y)))) +;;; +%1 + (defgeneric +%1 (x &optional r) (:method ((x number) &optional r) (declare (ignore r)) @@ -320,6 +339,9 @@ ) + +;;; *%1 + (defgeneric *%1 (x &optional r) (:method ((x number) &optional r) (declare (ignore r)) @@ -335,6 +357,9 @@ ) + +;;; -%1 + (defgeneric -%1 (x &optional r) (:method ((x number) &optional r) (declare (ignore r)) @@ -350,6 +375,9 @@ ) + +;;; /%1 + (defgeneric /%1 (x &optional r) (:method ((x number) &optional r) (declare (ignore r)) @@ -365,6 +393,8 @@ ) +;;; =2% + (defgeneric =%2 (x y) (:method ((x number) (y number)) (cl:= x y)) @@ -396,11 +426,17 @@ ) +;;; =%1 + (defgeneric =%1 (x) (:method ((x t)) T) ) +;;;--------------------------------------------------------------------------- +;;; Other operations + + (defgeneric gcd%2 (x y) (:method ((x integer) (y integer)) (cl:gcd x y))) Modified: trunk/common-math/numerics/linear-algebra/matrix.lisp ============================================================================== --- trunk/common-math/numerics/linear-algebra/matrix.lisp (original) +++ trunk/common-math/numerics/linear-algebra/matrix.lisp Thu Aug 17 12:39:36 2006 @@ -773,20 +773,10 @@ (defmethod *%2 ((y array) (x number) &optional (r (copy-matrix y) r-supplied-p)) + (declare (ignore r-supplied-p)) (*%2 x y r)) -;;; The next one breaks the return type convention. -#| Defined in 'vector.lisp'. -(defmethod *%2 ((x vector) (y vector) &optional r) - (declare (ignore r)) - (assert (conforming-*-dimensions-p x y nil)) - (let ((result 0)) - (dotimes (i (length x) result) - (setf result (+%2 result (* (aref x i) (aref y i))))))) -|# - - ;;;--------------------------------------------------------------------------- ;;; Division. ;;; Only the simple form of division is implemented here. @@ -810,6 +800,8 @@ ;;;--------------------------------------------------------------------------- ;;; Element-wise operations. +;;; .*%2 + (defmethod .*%2 ((x number) (y matrix) &optional (r (copy-matrix y) r-supplied-p)) (when r-supplied-p (assert (shape-equal-p y r))) @@ -850,6 +842,187 @@ (*%2 (row-major-aref x i) (row-major-aref y i))))) +(defmethod .*%2 ((x matrix) (y matrix) &optional (r (copy-matrix y) r-supplied-p)) + (when r-supplied-p (assert (shape-equal-p y r))) + (assert (shape-equal-p x y)) + + (with-slots ((x-data data)) x + (with-slots ((y-data data)) y + (with-slots ((result data)) r + + (dotimes (i (array-total-size result) r) + (setf (row-major-aref result i) + (*%2 (row-major-aref x-data i) (row-major-aref y-data i))))) + ))) + + +(defmethod .*%2 ((x matrix) (y array) &optional (r (copy-matrix y) r-supplied-p)) + (when r-supplied-p (assert (shape-equal-p y r))) + (assert (shape-equal-p x y)) + (.*%2 (matrix-data x) y r)) + + +(defmethod .*%2 ((x array) (y matrix) &optional (r (copy-matrix y) r-supplied-p)) + (when r-supplied-p (assert (shape-equal-p y r))) + (assert (shape-equal-p x y)) + + (with-slots ((y-data data)) y + (with-slots ((result data)) r + (dotimes (i (array-total-size result) r) + (setf (row-major-aref result i) + (*%2 (row-major-aref x i) (row-major-aref y-data i)))) + ))) + + + +;;; ./%2 + +(defmethod ./%2 ((x number) (y matrix) + &optional (r (copy-matrix y) r-supplied-p)) + (when r-supplied-p (assert (shape-equal-p y r))) + (with-slots (data) y + (with-slots ((result data)) r + (dotimes (i (array-total-size data) r) + (setf (row-major-aref result i) (/%2 x (row-major-aref data i)))) + ))) + + +(defmethod ./%2 ((y matrix) (x number) &optional (r (copy-matrix y))) + (.*%2 y (cl:/ x) r)) + + +(defmethod ./%2 ((x number) (y array) + &optional (r (copy-matrix y) r-supplied-p)) + (assert (matrix-array-p y)) + (when r-supplied-p (assert (shape-equal-p y r))) + (let* ((data y) + (result r) + ) + (dotimes (i (array-total-size data) r) + (setf (row-major-aref result i) (/%2 x (row-major-aref data i)))) + )) + + +(defmethod ./%2 ((y array) (x number) &optional (r (copy-matrix y))) + (.*%2 y (cl:/ x) r)) + + + +(defmethod ./%2 ((x array) (y array) &optional (r (copy-matrix y) r-supplied-p)) + (when r-supplied-p (assert (shape-equal-p y r))) + (assert (shape-equal-p x y)) + + (dotimes (i (array-total-size y) r) + (setf (row-major-aref r i) + (/%2 (row-major-aref x i) (row-major-aref y i))))) + + +(defmethod ./%2 ((x matrix) (y matrix) &optional (r (copy-matrix y) r-supplied-p)) + (when r-supplied-p (assert (shape-equal-p y r))) + (assert (shape-equal-p x y)) + + (with-slots ((x-data data)) x + (with-slots ((y-data data)) y + (with-slots ((result data)) r + + (dotimes (i (array-total-size result) r) + (setf (row-major-aref result i) + (/%2 (row-major-aref x-data i) (row-major-aref y-data i))))) + ))) + + +(defmethod ./%2 ((x matrix) (y array) &optional (r (copy-matrix y) r-supplied-p)) + (when r-supplied-p (assert (shape-equal-p y r))) + (assert (shape-equal-p x y)) + (./%2 (matrix-data x) y r)) + + +(defmethod ./%2 ((x array) (y matrix) &optional (r (copy-matrix y) r-supplied-p)) + (when r-supplied-p (assert (shape-equal-p y r))) + (assert (shape-equal-p x y)) + + (with-slots ((y-data data)) y + (with-slots ((result data)) r + (dotimes (i (array-total-size result) r) + (setf (row-major-aref result i) + (/%2 (row-major-aref x i) (row-major-aref y-data i)))) + ))) + + +;;; .\\%2 + +(defmethod .\\%2 ((x number) (y matrix) + &optional (r (copy-matrix y) r-supplied-p)) + (when r-supplied-p (assert (shape-equal-p y r))) + (with-slots (data) y + (with-slots ((result data)) r + (dotimes (i (array-total-size data) r) + (setf (row-major-aref result i) (/%2 (row-major-aref data i) x))) + ))) + + +(defmethod .\\%2 ((y matrix) (x number) &optional (r (copy-matrix y))) + (./%2 x y r)) + + +(defmethod .\\%2 ((x number) (y array) + &optional (r (copy-matrix y) r-supplied-p)) + (assert (matrix-array-p y)) + (when r-supplied-p (assert (shape-equal-p y r))) + (let* ((data y) + (result r) + ) + (dotimes (i (array-total-size data) r) + (setf (row-major-aref result i) (/%2 (row-major-aref data i) x))) + )) + + +(defmethod .\\%2 ((y array) (x number) &optional (r (copy-matrix y))) + (./%2 x y r)) + + + +(defmethod .\\%2 ((x array) (y array) &optional (r (copy-matrix y) r-supplied-p)) + (when r-supplied-p (assert (shape-equal-p y r))) + (assert (shape-equal-p x y)) + + (dotimes (i (array-total-size y) r) + (setf (row-major-aref r i) + (/%2 (row-major-aref y i) (row-major-aref x i))))) + + +(defmethod .\\%2 ((x matrix) (y matrix) &optional (r (copy-matrix y) r-supplied-p)) + (when r-supplied-p (assert (shape-equal-p y r))) + (assert (shape-equal-p x y)) + + (with-slots ((x-data data)) x + (with-slots ((y-data data)) y + (with-slots ((result data)) r + + (dotimes (i (array-total-size result) r) + (setf (row-major-aref result i) + (/%2 (row-major-aref y-data i) (row-major-aref x-data i))))) + ))) + + +(defmethod .\\%2 ((x matrix) (y array) &optional (r (copy-matrix y) r-supplied-p)) + (when r-supplied-p (assert (shape-equal-p y r))) + (assert (shape-equal-p x y)) + (.\\%2 (matrix-data x) y r)) + + +(defmethod .\\%2 ((x array) (y matrix) &optional (r (copy-matrix y) r-supplied-p)) + (when r-supplied-p (assert (shape-equal-p y r))) + (assert (shape-equal-p x y)) + + (with-slots ((y-data data)) y + (with-slots ((result data)) r + (dotimes (i (array-total-size result) r) + (setf (row-major-aref result i) + (/%2 (row-major-aref y-data i) (row-major-aref x i)))) + ))) + + ;;;--------------------------------------------------------------------------- ;;; Transpose

1 0

[common-math-cvs] r1 - in trunk: . common-math common-math/numerics common-math/numerics/linear-algebra
by mantoniotti＠common-lisp.net 16 Aug '06

16 Aug '06

Author: mantoniotti Date: Wed Aug 16 17:07:41 2006 New Revision: 1 Added: trunk/ trunk/common-math/ trunk/common-math/common-math-pkg.lisp trunk/common-math/common-math-user-pkg.lisp trunk/common-math/common-math.lisp trunk/common-math/common-math.system trunk/common-math/numerics/ trunk/common-math/numerics/linear-algebra/ trunk/common-math/numerics/linear-algebra/common-math-extra.lisp trunk/common-math/numerics/linear-algebra/linear-algebra-pkg.lisp trunk/common-math/numerics/linear-algebra/linear-algebra.system trunk/common-math/numerics/linear-algebra/lu-decomposition.lisp trunk/common-math/numerics/linear-algebra/matrix-pkg.lisp trunk/common-math/numerics/linear-algebra/matrix.lisp trunk/common-math/numerics/linear-algebra/vector-matrix-conditions.lisp trunk/common-math/numerics/linear-algebra/vector-matrix.lisp trunk/common-math/numerics/linear-algebra/vector.lisp trunk/common-math/numerics/numerics.system Log: Initial import. Added: trunk/common-math/common-math-pkg.lisp ============================================================================== --- (empty file) +++ trunk/common-math/common-math-pkg.lisp Wed Aug 16 17:07:41 2006 @@ -0,0 +1,47 @@ +;;;; -*- Mode: Lisp -*- + +;;;; common-math-pkg.lisp -- + +(defpackage "COMMON-MATH" (:use "COMMON-LISP") + (:documentation "The COMMON MATH Package. +This package contains `top level' definitions of a `generic math' +package which can be used as a drop in replacement for the standard +CL mathematics routines. Note however, that this will of course +come at a price. + +The COMMON-MATH package addresses the problem of having a generic and +as complete as possible mathematics interface layer over a number of +different algorithms and representations. In this sense it aspires +to bring to CL the flexibility of generic math operations, while +preserving the low level hooks that could be used for (relative) efficiency.") + + (:nicknames "CL.MATH" "CL.MATHEMATICS") + + (:shadow "=" "+" "-" "*" "/" ">" "<" ">=" "<=") + (:shadow "ZEROP" "EXPT" "GCD") + + (:export "=" "+" "-" "*" "/" ">" "<" ">=" "<=") + (:export "ZEROP" "EXPT" "GCD") + + (:export "+POSITIVE-INFINITY+" "+NEGATIVE-INFINITY+") + + (:export "INCREASING" "DECREASING") + + (:export + "+%1" "+%2" + "*%1" "*%2" + "-%1" "-%2" + "/%1" "/%2" + "=%1" "=%2" + "<%2" + ">%2" + "<=%2" + ">=%2" + ) + + ;; Conditions + (:export "UNDEFINED-OPERATION") + + ) + +;;;; end of file -- common-math-pkg.lisp -- Added: trunk/common-math/common-math-user-pkg.lisp ============================================================================== --- (empty file) +++ trunk/common-math/common-math-user-pkg.lisp Wed Aug 16 17:07:41 2006 @@ -0,0 +1,12 @@ +;;;; -*- Mode: Lisp -*- + +;;;; common-math-user-pkg.lisp -- +;;;; The COMMON-MATH-USER package. + +(defpackage "COMMON-MATH-USER" (:use "COMMON-MATH" "CL" #| "LINALG" |# ) + (:shadowing-import-from "COMMON-MATH" "+" "-" "*" "/") + (:shadowing-import-from "COMMON-MATH" "<" ">" "<=" ">=" "=") + (:shadowing-import-from "COMMON-MATH" "ZEROP" "EXPT" "GCD") + ) + +;;;; end of file -- common-math-user-pkg.lisp -- \ No newline at end of file Added: trunk/common-math/common-math.lisp ============================================================================== --- (empty file) +++ trunk/common-math/common-math.lisp Wed Aug 16 17:07:41 2006 @@ -0,0 +1,603 @@ +;;;; -*- Mode: Lisp -*- + +;;;; common-math.lisp -- + +(in-package "CL.MATHEMATICS") + +;;; undefined-operation -- + +(define-condition undefined-operation (error) + ((operator :reader undefined-operation-operator + :initarg :operator) + (arguments :reader undefined-operation-arguments + :initarg :arguments) + ) + (:documentation "The Undefined Operation Condition.") + (:report (lambda (uoc stream) + (format stream "Undefined operation ~S called with ~S." + (undefined-operation-operator uoc) + (undefined-operation-arguments uoc)))) + (:default-initargs :arguments () :operator nil)) + + +;;;--------------------------------------------------------------------------- +;;; Constants and parameters. + +(defconstant +negative-infinity+ '+negative-infinity+ + "The symbolic +NEGATIVE-INFINITY+ constant. +A representation of negative infinity.") +(defconstant +positive-infinity+ '+positive-infinity+ + "The symbolic +POSITIVE-INFINITY+ constant. +A representation of positive infinity.") + +(defparameter *ignore-comparison-errors-p* nil) + + +;;; (defconstant +negative-infinity+ most-negative-fixnum) +;;; (defconstant +positive-infinity+ most-positive-fixnum) + + +;;;--------------------------------------------------------------------------- +;;; Generic functions. +;;; Note that the generic functions corresponding to traditional math +;;; operators are named using a special notation to indicate their +;;; arity. Such 'operator' generic function are named as +;;; +;;; <op>%<arity> +;;; +;;; The #\% character in the name is used to suggest "low level" as +;;; usual. The Prolog-like #\/ could have been used instead but the +;;; #\% has more tradition in CL circles. The #\@ would have made +;;; more sense, but some CL implementations already use it for other +;;; purposes and setting up a separate readtable seems a little bit +;;; too much for the time being. + +(defgeneric <%2 (x y) + (:documentation "The <%2 generic function. +The binary LESS generic function which is specialized for various +combinations of argument types. + +At a minimum there is a guarantee that the method on two NUMBER +arguments dispatching on CL:< is defined.") + (:method ((x number) (y number)) + (cl:< x y)) + + (:method ((x number) (y (eql +negative-infinity+))) + nil) + + (:method ((y (eql +negative-infinity+)) (x number)) + T) + + (:method ((y (eql +negative-infinity+)) (x (eql +negative-infinity+))) + (unless *ignore-comparison-errors-p* + (error 'undefined-operation :operator 'lessp%2))) + + (:method ((y (eql +negative-infinity+)) (x (eql +positive-infinity+))) + T) + + (:method ((y (eql +positive-infinity+)) (x (eql +negative-infinity+))) + nil) + + (:method ((x number) (y (eql +positive-infinity+))) + T) + + (:method ((y (eql +positive-infinity+)) (x number)) + nil) + + (:method ((y (eql +positive-infinity+)) (x (eql +positive-infinity+))) + (unless *ignore-comparison-errors-p* + (error 'undefined-operation :operator 'lessp%2))) + ) + + +(defgeneric >%2 (x y) + (:method ((x number) (y number)) + (cl:> x y)) + + (:method ((x number) (y (eql +negative-infinity+))) + T) + + (:method ((y (eql +negative-infinity+)) (x number)) + nil) + + (:method ((y (eql +negative-infinity+)) (x (eql +negative-infinity+))) + (unless *ignore-comparison-errors-p* + (error 'undefined-operation :operator 'lessp%2))) + + (:method ((y (eql +negative-infinity+)) (x (eql +positive-infinity+))) + nil) + + (:method ((y (eql +positive-infinity+)) (x (eql +negative-infinity+))) + T) + + (:method ((x number) (y (eql +positive-infinity+))) + nil) + + (:method ((y (eql +positive-infinity+)) (x number)) + T) + + (:method ((y (eql +positive-infinity+)) (x (eql +positive-infinity+))) + (unless *ignore-comparison-errors-p* + (error 'undefined-operation :operator 'lessp%2))) + ) + + +(defmethod <%2 ((x symbol) (y symbol)) + (string< x y)) + + +(defmethod >%2 ((x symbol) (y symbol)) + (string> y x)) + + +(defgeneric <=%2 (x y) + (:method ((x number) (y number)) + (cl:<= x y)) + + (:method ((x number) (y (eql +negative-infinity+))) + nil) + + (:method ((y (eql +negative-infinity+)) (x number)) + T) + + (:method ((y (eql +negative-infinity+)) (x (eql +negative-infinity+))) + (unless *ignore-comparison-errors-p* + (error 'undefined-operation :operator 'lessp%2))) + + (:method ((y (eql +negative-infinity+)) (x (eql +positive-infinity+))) + T) + + (:method ((y (eql +positive-infinity+)) (x (eql +negative-infinity+))) + nil) + + (:method ((x number) (y (eql +positive-infinity+))) + T) + + (:method ((y (eql +positive-infinity+)) (x number)) + nil) + + (:method ((y (eql +positive-infinity+)) (x (eql +positive-infinity+))) + (unless *ignore-comparison-errors-p* + (error 'undefined-operation :operator 'lessp%2))) + ) + + +(defgeneric >=%2 (x y) + (:method ((x number) (y number)) + (cl:>= x y)) + + (:method ((x number) (y (eql +negative-infinity+))) + T) + + (:method ((y (eql +negative-infinity+)) (x number)) + nil) + + (:method ((y (eql +negative-infinity+)) (x (eql +negative-infinity+))) + (unless *ignore-comparison-errors-p* + (error 'undefined-operation :operator 'lessp%2))) + + (:method ((y (eql +negative-infinity+)) (x (eql +positive-infinity+))) + nil) + + (:method ((y (eql +positive-infinity+)) (x (eql +negative-infinity+))) + T) + + (:method ((x number) (y (eql +positive-infinity+))) + nil) + + (:method ((y (eql +positive-infinity+)) (x number)) + T) + + (:method ((y (eql +positive-infinity+)) (x (eql +positive-infinity+))) + (unless *ignore-comparison-errors-p* + (error 'undefined-operation :operator 'lessp%2))) + ) + + +(defgeneric +%2 (x y &optional r) + (:method ((x number) (y number) &optional r) + (declare (ignore r)) + (cl:+ x y)) + + (:method ((x number) (y (eql +negative-infinity+)) &optional r) + (declare (ignore r)) + +negative-infinity+) + + (:method ((y (eql +negative-infinity+)) (x number) &optional r) + (declare (ignore r)) + +negative-infinity+) + + (:method ((y (eql +negative-infinity+)) (x (eql +negative-infinity+)) + &optional r) + (declare (ignore r)) + +negative-infinity+) + + (:method ((y (eql +negative-infinity+)) (x (eql +positive-infinity+)) + &optional r) + (declare (ignore r)) + (error 'undefined-operation :operator '+%2)) + + (:method ((x number) (y (eql +positive-infinity+)) &optional r) + (declare (ignore r)) + +positive-infinity+) + + (:method ((y (eql +positive-infinity+)) (x number) &optional r) + (declare (ignore r)) + +positive-infinity+) + + (:method ((y (eql +positive-infinity+)) (x (eql +positive-infinity+)) + &optional r) + (declare (ignore r)) + +negative-infinity+) + + (:method ((y (eql +positive-infinity+)) (x (eql +negative-infinity+)) + &optional r) + (declare (ignore r)) + (error 'undefined-operation :operator '+%2)) + ) + + +(defgeneric *%2 (x y &optional r) + (:method ((x number) (y number) &optional r) + (declare (ignore r)) + (cl:* x y)) + + (:method ((x number) (y (eql +negative-infinity+)) &optional r) + (declare (ignore r)) + (cond ((zerop x) (error 'undefined-operation :operator '*%2)) + ((plusp x) +negative-infinity+) + (t +positive-infinity+))) + + (:method ((y (eql +negative-infinity+)) (x number) &optional r) + (declare (ignore r)) + (cond ((zerop x) (error 'undefined-operation :operator '*%2)) + ((plusp x) +negative-infinity+) + (t +positive-infinity+))) + + + (:method ((y (eql +negative-infinity+)) (x (eql +negative-infinity+)) + &optional r) + (declare (ignore r)) + +positive-infinity+) + + (:method ((y (eql +negative-infinity+)) (x (eql +positive-infinity+)) + &optional r) + (declare (ignore r)) + +negative-infinity+) + + (:method ((x number) (y (eql +positive-infinity+)) &optional r) + (declare (ignore r)) + (cond ((zerop x) (error 'undefined-operation :operator '*%2)) + ((plusp x) +positive-infinity+) + (t +negative-infinity+))) + + (:method ((y (eql +positive-infinity+)) (x number) &optional r) + (declare (ignore r)) + (cond ((zerop x) (error 'undefined-operation :operator '*%2)) + ((plusp x) +positive-infinity+) + (t +negative-infinity+))) + + (:method ((y (eql +positive-infinity+)) (x (eql +positive-infinity+)) + &optional r) + (declare (ignore r)) + +positive-infinity+) + + (:method ((y (eql +positive-infinity+)) (x (eql +negative-infinity+)) + &optional r) + (declare (ignore r)) + +negative-infinity+) + ) + + +(defgeneric -%2 (x y &optional r) + (:method ((x t) (y t) &optional r) + (declare (ignore r)) + (+%2 x (-%1 y)))) + + +(defgeneric /%2 (x y &optional r) + (:method ((x number) (y number) &optional r) + (declare (ignore r)) + (cl:/ x y)) + + (:method ((x t) (y t) &optional r) + (declare (ignore r)) + (*%2 x (/%1 y)))) + + +(defgeneric +%1 (x &optional r) + (:method ((x number) &optional r) + (declare (ignore r)) + (cl:+ x)) + + (:method ((x (eql +positive-infinity+)) &optional r) + (declare (ignore r)) + +positive-infinity+) + + (:method ((x (eql +negative-infinity+)) &optional r) + (declare (ignore r)) + +negative-infinity+) + + ) + +(defgeneric *%1 (x &optional r) + (:method ((x number) &optional r) + (declare (ignore r)) + (cl:* x)) + + (:method ((x (eql +positive-infinity+)) &optional r) + (declare (ignore r)) + +positive-infinity+) + + (:method ((x (eql +negative-infinity+)) &optional r) + (declare (ignore r)) + +negative-infinity+) + + ) + +(defgeneric -%1 (x &optional r) + (:method ((x number) &optional r) + (declare (ignore r)) + (cl:- x)) + + (:method ((x (eql +positive-infinity+)) &optional r) + (declare (ignore r)) + +negative-infinity+) + + (:method ((x (eql +negative-infinity+)) &optional r) + (declare (ignore r)) + +positive-infinity+) + + ) + +(defgeneric /%1 (x &optional r) + (:method ((x number) &optional r) + (declare (ignore r)) + (cl:/ x)) + + (:method ((x (eql +positive-infinity+)) &optional r) + (declare (ignore r)) + 0) + + (:method ((x (eql +negative-infinity+)) &optional r) + (declare (ignore r)) + 0) + ) + + +(defgeneric =%2 (x y) + (:method ((x number) (y number)) (cl:= x y)) + + (:method ((x number) (y (eql +negative-infinity+))) + nil) + + (:method ((y (eql +negative-infinity+)) (x number)) + nil) + + (:method ((y (eql +negative-infinity+)) (x (eql +negative-infinity+))) + (unless *ignore-comparison-errors-p* + (error 'undefined-operation :operator '=%2))) + + (:method ((y (eql +negative-infinity+)) (x (eql +positive-infinity+))) + nil) + + (:method ((y (eql +positive-infinity+)) (x (eql +negative-infinity+))) + nil) + + (:method ((x number) (y (eql +positive-infinity+))) + nil) + + (:method ((y (eql +positive-infinity+)) (x number)) + nil) + + (:method ((y (eql +positive-infinity+)) (x (eql +positive-infinity+))) + (unless *ignore-comparison-errors-p* + (error 'undefined-operation :operator '=%2))) + ) + + +(defgeneric =%1 (x) + (:method ((x t)) T) + ) + + +(defgeneric gcd%2 (x y) + (:method ((x integer) (y integer)) + (cl:gcd x y))) + + +;;;--------------------------------------------------------------------------- +;;; Redefined Common Lisp operators. + +(defun + (&rest args) + (if (null args) + (cl:+) + (let ((n-args (list-length args))) + (cond ((cl:= n-args 1) + (+%1 (first args))) + ((cl:= n-args 2) + (+%2 (first args) (second args))) + (t (+%2 (first args) (apply #'+ (rest args)))) + )))) + + +(defun - (arg1 &rest args) + (if (null args) + (-%1 arg1) + (let ((n-args (list-length args))) + (if (cl:= n-args 1) + (-%2 arg1 (first args)) + ;; The loop is done in order not to assume anything about + ;; the nature of summation. + (loop for arg in args + for minuend =(-%2 arg1 arg) then (-%2 minuend arg) + finally (return minuend)))) + )) + + +(defun * (&rest args) + (if (null args) + (cl:*) + (let ((n-args (list-length args))) + (cond ((cl:= n-args 1) + (*%1 (first args))) + ((cl:= n-args 2) + (*%2 (first args) (second args))) + (t (*%2 (first args) (apply #'* (rest args)))) + )))) + + +(defun / (arg1 &rest args) + (if (null args) + (/%1 arg1) + (let ((n-args (list-length args))) + (if (cl:= n-args 1) + (/%2 arg1 (first args)) + ;; The loop is done in order not to assume anything about + ;; the nature of divisions. + (loop for arg in args + for dividend = (/%2 arg1 arg) then (/%2 dividend arg) + finally (return dividend)))) + )) + + +(defun = (arg1 &rest args) + (if (null args) + (=%1 arg1) + (let ((n-args (list-length args))) + (if (cl:= n-args 1) + (=%2 arg1 (first args)) + ;; The loop is done in order not to assume anything about + ;; the nature of the equality test. + (loop for arg in args + always (=%2 arg1 arg)))) + )) + + +(defun < (arg1 &rest args) + (if (null args) + ;; (<%1 arg1) + t + (let ((n-args (list-length args))) + (if (cl:= n-args 1) + (<%2 arg1 (first args)) + ;; The loop is done in order not to assume anything about + ;; the nature of the equality test. + (loop for a1 = arg1 then a2 + for a2 in args + always (<%2 a1 a2)))) + )) + + +(defun > (arg1 &rest args) + (if (null args) + t + (let ((n-args (list-length args))) + (if (cl:= n-args 1) + (>%2 arg1 (first args)) + ;; The loop is done in order not to assume anything about + ;; the nature of the equality test. + (loop for a1 = arg1 then a2 + for a2 in args + always (>%2 a1 a2)))) + )) + + +(defun <= (arg1 &rest args) + (if (null args) + t + (let ((n-args (list-length args))) + (if (cl:= n-args 1) + (<=%2 arg1 (first args)) + ;; The loop is done in order not to assume anything about + ;; the nature of the equality test. + (loop for a1 = arg1 then a2 + for a2 in args + always (<=%2 a1 a2)))) + )) + + +(defun >= (arg1 &rest args) + (if (null args) + t + (let ((n-args (list-length args))) + (if (cl:= n-args 1) + (>=%2 arg1 (first args)) + ;; The loop is done in order not to assume anything about + ;; the nature of the equality test. + (loop for a1 = arg1 then a2 + for a2 in args + always (>=%2 a1 a2)))) + )) + + +(defun gcd (&rest args) + (if (null args) + (cl:gcd) + (let ((n-args (list-length args))) + (cond ((cl:= n-args 1) + (gcd%2 (first args) (first args))) + ((cl:= n-args 1) + (gcd%2 (first args) (second args))) + (t + (apply #'gcd + (gcd%2 (first args) (second args)) + (cddr args))) + )))) + + +(defgeneric expt (base power) + (:method ((base number) (power number)) + (cl:expt base power))) + + +(defgeneric zerop (n) + (:method ((n number)) + (cl:zerop n)) + + (:method ((x (eql +positive-infinity+))) + nil) + + (:method ((x (eql +negative-infinity+))) + nil) + ) + + + +;;;--------------------------------------------------------------------------- +;;; Other functions. + + +(defun increasing (x) + (declare (type sequence x)) + (loop for i from 0 below (1- (length x)) + for a1 = (elt x 0) then a2 + for a2 = (elt x (1+ i)) + always (<=%2 a1 a2))) + + +(defun decreasing (x) + (declare (type sequence x)) + (loop for i from 0 below (1- (length x)) + for a1 = (elt x 0) then a2 + for a2 = (elt x (1+ i)) + always (>=%2 a1 a2))) + + + +(defun strictly-increasing (x) + (declare (type sequence x)) + (loop for i from 0 below (1- (length x)) + for a1 = (elt x 0) then a2 + for a2 = (elt x (1+ i)) + always (<%2 a1 a2))) + + +(defun strictly-decreasing (x) + (declare (type sequence x)) + (loop for i from 0 below (1- (length x)) + for a1 = (elt x 0) then a2 + for a2 = (elt x (1+ i)) + always (>%2 a1 a2))) + +;;; end of file -- common-math.lisp -- Added: trunk/common-math/common-math.system ============================================================================== --- (empty file) +++ trunk/common-math/common-math.system Wed Aug 16 17:07:41 2006 @@ -0,0 +1,24 @@ +;;;; -*- Mode: Lisp -*- + +;;;; common-math.system -- + +(eval-when (:load-toplevel :execute) + (mk:add-registry-location (make-pathname :name nil + :type nil + :defaults *load-truename*)) + ) + + +(mk:defsystem "common-math" + :components ((:file "common-math-pkg") + (:file "common-math" + :depends-on ("common-math-pkg")) + (:system "numerics") + ;; (:system "polynomials") + ;; (:system "formulas") + (:file "common-math-user-pkg" + :depends-on ("numerics")) + ) + ) + +;;;; end of file -- common-math.system -- Added: trunk/common-math/numerics/linear-algebra/common-math-extra.lisp ============================================================================== --- (empty file) +++ trunk/common-math/numerics/linear-algebra/common-math-extra.lisp Wed Aug 16 17:07:41 2006 @@ -0,0 +1,98 @@ +;;;; -*- Mode: Lisp -*- + +;;;; common-math-extra.lisp -- + +(in-package "CL.MATH.LINEAR-ALGEBRA") + +;;;--------------------------------------------------------------------------- +;;; More functions. Essentially the top level element-wise array operations. + +(defgeneric .*%2 (x y &optional r) + (:method ((x number) (y number) &optional r) + (declare (ignore r)) + (*%2 x y)) + + (:method ((x (eql +negative-infinity+)) (y number) &optional r) + (declare (ignore r)) + (*%2 x y)) + + (:method ((x (eql +positive-infinity+)) (y number) &optional r) + (declare (ignore r)) + (*%2 x y)) + + + (:method ((y number) (x (eql +negative-infinity+)) &optional r) + (declare (ignore r)) + (*%2 x y)) + + + (:method ((y number) (x (eql +positive-infinity+)) &optional r) + (declare (ignore r)) + (*%2 x y)) + + + (:method ((y (eql +positive-infinity+)) (x (eql +positive-infinity+)) &optional r) + (declare (ignore r)) + +positive-infinity+) + + (:method ((y (eql +negative-infinity+)) (x (eql +negative-infinity+)) &optional r) + (declare (ignore r)) + +positive-infinity+) + + (:method ((y (eql +positive-infinity+)) (x (eql +negative-infinity+)) &optional r) + (declare (ignore r)) + +negative-infinity+) + + (:method ((x (eql +negative-infinity+)) (y (eql +positive-infinity+)) &optional r) + (declare (ignore r)) + +negative-infinity+) + ) + + +(defgeneric .*%1 (x &optional r) + (:method ((x number) &optional r) + (declare (ignore r)) + (*%1 x)) + + (:method ((x (eql +positive-infinity+)) &optional r) + (declare (ignore r)) + +positive-infinity+) + + (:method ((x (eql +negative-infinity+)) &optional r) + (declare (ignore r)) + +negative-infinity+) + ) + + +(defgeneric ./%2 (x y &optional r)) + +(defgeneric ./%1 (x &optional r)) + + +;;;--------------------------------------------------------------------------- +;;; The generic top-level functions. + +(defun .* (&rest args) + (if (null args) + (cl:*) + (let ((n-args (list-length args))) + (cond ((cl:= n-args 1) (.*%1 (first args))) + ((cl:= n-args 2) (.*%2 (first args) (second args))) + (t (.*%2 (first args) (apply #'.* (rest args)))) + )))) + + +(defun ./ (arg1 &rest args) + (if (null args) + (./%1 arg1) + (let ((n-args (list-length args))) + (if (cl:= n-args 1) + (./%2 arg1 (first args)) + ;; The loop is done in order not to assume anything about + ;; the nature of divisions. + (loop for arg in args + for dividend = (./%2 arg1 arg) then (./%2 dividend arg) + finally (return dividend)))) + )) + +;;;; end of file -- common-math-extra.lisp -- Added: trunk/common-math/numerics/linear-algebra/linear-algebra-pkg.lisp ============================================================================== --- (empty file) +++ trunk/common-math/numerics/linear-algebra/linear-algebra-pkg.lisp Wed Aug 16 17:07:41 2006 @@ -0,0 +1,19 @@ +;;; -*- Mode: Lisp -*- + +(defpackage "CL.MATH.LINEAR-ALGEBRA" (:use "CL.MATH" "COMMON-LISP") + (:nicknames "CL-MATH-LINALG" "LINALG") + + (:shadow cl:zerop cl:expt cl:gcd) + + (:shadowing-import-from "CL.MATH" "+" "-" "*" "/" "=" ">" "<" ">=" "<=") + + (:export "=" "+" "-" "*" "/" ">" "<" ">=" "<=") + (:export "ZEROP" "EXPT" "GCD") + + (:export ".*" "./" "EXPT*") + + (:export ".*%2" ".*%1" "./%2" "./%1") + + ) + +;;; end of file -- linear-algebra-pkg.lisp -- Added: trunk/common-math/numerics/linear-algebra/linear-algebra.system ============================================================================== --- (empty file) +++ trunk/common-math/numerics/linear-algebra/linear-algebra.system Wed Aug 16 17:07:41 2006 @@ -0,0 +1,13 @@ +;;;; -*- Mode: Lisp -*- + +;;;; linear-algebra.system -- + +(mk:defsystem "linear-algebra" + :components ("linear-algebra-pkg" + "vector-matrix-conditions" + "common-math-extra" + "vector" + "matrix" + "lu-decomposition")) + +;;;; end of file -- linear-algebra.system -- Added: trunk/common-math/numerics/linear-algebra/lu-decomposition.lisp ============================================================================== --- (empty file) +++ trunk/common-math/numerics/linear-algebra/lu-decomposition.lisp Wed Aug 16 17:07:41 2006 @@ -0,0 +1,326 @@ +;;;; -*- Mode: Lisp -*- + +;;;; lu-decomposition.lisp -- + +(in-package "CL.MATH.LINEAR-ALGEBRA") + +(define-condition singular-matrix-in-operation (arithmetic-error) + () + (:report (lambda (cnd stream) + (format stream "Singular matrix passed to operation ~S." + (arithmetic-error-operation cnd))))) + + +(defparameter *epsilon* 1.0e-20) ; Change this to CL constant. + + +(defgeneric square-matrix-p (m) + (:method ((m array)) + (and (matrix-array-p m) + (cl:= (array-dimension m 0) (array-dimension m 1)))) + (:method ((m matrix)) + (cl:= (matrix-row-number m) (matrix-column-number m)))) + + + +;;; Different base type. +;;; Probably a macrolet is in order. + +(defun lu-decompose (m &optional + (lu (make-array (array-dimensions m) + :element-type (array-element-type m) + :initial-element (coerce 0 (array-element-type m)))) + (permutations (make-array (array-dimension m 0) + :initial-element 0 + :element-type 'fixnum))) + (declare (type (array number (cl:* cl:*)) m lu) + (type (vector fixnum) permutations)) + + (assert (square-matrix-p m)) + (assert (square-matrix-p lu)) + (assert (cl:= (length permutations) (array-dimension m 0))) + + (let* ((n (array-dimension m 0)) + (row-scaling-vector (make-array n + :initial-element (coerce 0 (array-element-type m)) + :element-type (array-element-type m))) + (parity 1) + (imax 0) + ) + (declare (dynamic-extent row-scaling-vector) + (type fixnum imax n) + (type (member -1 1) parity) + (type (simple-array number (cl:*)) row-scaling-vector)) + + ;; Copy M into LU and then just operate the in-place algorithm on LU. + ;; If M is EQ to LU, that means that the argumen was supplied and + ;; that the user actually wants a destructive operation. Copying + ;; can be avoided in this case. + (unless (eq m lu) + (dotimes (i n) + (dotimes (j n) + (setf (aref lu i j) (aref m i j))))) + + + ;; Loop over rows to get the implicit scaling information. + ;; Meanwhile check that the matrix is non-singular. + (dotimes (i n) + (loop for k of-type fixnum from 0 below n + maximize (abs (aref lu i k)) into big ;; of-type (number) ; not CL. + finally + (when (cl:zerop big) + (error 'singular-matrix-in-operation + :operation 'lu-decompose + :operands (list m))) + (setf (aref row-scaling-vector i) (/ big)))) + + ;; Loop over the columns of Crout's method. + (dotimes (j n) + (block operate-on-column + ;; (format t "Column ~D IMAX = ~D~%" j imax) + (dotimes (i j) + (let ((sum (aref lu i j))) + (declare (type number sum)) + (dotimes (k i) + (decf sum (* (aref lu i k) (aref lu k j)))) + (setf (aref lu i j) sum))) + + ;; Initialize the search for the largest pivot. + (loop with big of-type number = 0.0 + for i from j below n + do (let ((sum (aref lu i j))) + (declare (type number sum)) + (dotimes (k j) + (decf sum (* (aref lu i k) (aref lu k j)))) + (setf (aref lu i j) sum) + (let ((temp (* (aref row-scaling-vector i) (abs sum)))) + (declare (type number temp)) + (when (>= temp big) + (setf big temp + imax i))) + )) + + ;; (format t "Now IMAX = ~D~%" imax) + (when (/= imax j) + ;; We need to interchange rows. + (dotimes (k n) + (rotatef (aref lu imax k) (aref lu j k))) + (setf parity (- parity)) + (setf (aref row-scaling-vector imax) + (aref row-scaling-vector j))) + (setf (aref permutations j) imax) + ;; (format t "Now J = ~D IMAX = ~D~2%" j imax) + (when (cl:zerop (aref lu j j)) + (setf (aref lu j j) *epsilon*)) + + (when (/= j (1- n)) + (loop with pivot of-type number = (/ (aref lu j j)) + for i from (1+ j) below n + do (setf (aref lu i j) (* pivot (aref lu i j))))) + + )) + (values lu permutations parity) + )) + + +(defun lu-decompose/double-float (m &optional + (lu (make-array (array-dimensions m) + :element-type 'double-float + :initial-element 0.0d0)) + (permutations (make-array (array-dimension m 0) + :initial-element 0 + :element-type 'fixnum))) + (declare (type (array double-float (cl:* cl:*)) m lu) + (type (vector fixnum) permutations)) + + (assert (square-matrix-p m)) + (assert (square-matrix-p lu)) + (assert (cl:= (length permutations) (array-dimension m 0))) + + (let* ((n (array-dimension m 0)) + (row-scaling-vector (make-array n + :initial-element (coerce 0 (array-element-type m)) + :element-type (array-element-type m))) + (parity 1.0d0) + (imax 0) + ) + (declare (dynamic-extent row-scaling-vector) + (type fixnum n imax) + (type (double-float -1.01d0 1.01d0) parity) + (type (simple-array double-float (cl:*)) row-scaling-vector)) + + ;; Copy M into LU and then just operate the in-place algorithm on LU. + ;; If M is EQ to LU, that means that the argumen was supplied and + ;; that the user actually wants a destructive operation. Copying + ;; can be avoided in this case. + (unless (eq m lu) + (dotimes (i n) + (dotimes (j n) + (setf (aref lu i j) (aref m i j))))) + + + ;; Loop over rows to get the implicit scaling information. + ;; Meanwhile check that the matrix is non-singular. + (dotimes (i n) + (loop for k of-type fixnum from 0 below n + maximize (abs (aref lu i k)) into big of-type double-float + finally + (when (cl:zerop big) + (error 'singular-matrix-in-operation + :operation 'lu-decompose + :operands (list m))) + (setf (aref row-scaling-vector i) (/ big)))) + + ;; Loop over the columns of Crout's method. + (dotimes (j n) + (block operate-on-column + ;; (format t "Column ~D IMAX = ~D~%" j imax) + (dotimes (i j) + (let ((sum (aref lu i j))) + (declare (type double-float sum)) + (dotimes (k i) + (decf sum (* (aref lu i k) (aref lu k j)))) + (setf (aref lu i j) sum))) + + ;; Initialize the search for the largest pivot. + (loop with big of-type double-float = 0.0d0 + for i of-type fixnum from j below n + do (let ((sum (aref lu i j))) + (declare (type double-float sum)) + (dotimes (k j) + (decf sum (* (aref lu i k) (aref lu k j)))) + (setf (aref lu i j) sum) + (let ((temp (* (aref row-scaling-vector i) (abs sum)))) + (declare (type double-float temp)) + (when (>= temp big) + (setf big temp + imax i))) + )) + + ;; (format t "Now IMAX = ~D~%" imax) + (when (/= imax j) + ;; We need to interchange rows. + (dotimes (k n) + (rotatef (aref lu imax k) (aref lu j k))) + (setf parity (- parity)) + (setf (aref row-scaling-vector imax) + (aref row-scaling-vector j))) + (setf (aref permutations j) imax) + ;; (format t "Now J = ~D IMAX = ~D~2%" j imax) + (when (cl:zerop (aref lu j j)) + (setf (aref lu j j) *epsilon*)) + + (when (/= j (1- n)) + (loop with pivot of-type double-float = (/ (aref lu j j)) + for i from (1+ j) below n + do (setf (aref lu i j) (* pivot (aref lu i j))))) + + )) + (values lu permutations parity) + )) + + +(defun split-lu-matrix (lu) + (let ((l (make-array (array-dimensions lu) :initial-element 0)) + (u (make-array (array-dimensions lu) :initial-element 0)) + (n (array-dimension lu 0)) + ) + (dotimes (i n) + (dotimes (j i) + (setf (aref l i j) (aref lu i j))) + (setf (aref l i i) 1)) + + (dotimes (i n) + (loop for j from i below n + do (setf (aref u i j) (aref lu i j)))) + + (values l u))) + + +;;;--------------------------------------------------------------------------- +;;; lu-back-substitution + +(defun lu-back-substitution (lu permutations b) + (let ((n (array-dimension lu 0))) + (loop with ii of-type fixnum = 0 + for i of-type fixnum from 0 below n + for ip of-type number = (aref permutations i) + for sum of-type number = (aref b ip) + do (setf (aref b ip) (aref b i)) + do (cond ((not (cl:zerop ii)) + (loop for j of-type fixnum from ii below (1- i) + do (decf sum (* (aref lu i j) (aref b i))))) + ((not (cl:zerop sum)) + (setf (aref b i) sum))) + ) + (loop for i of-type fixnum from (1- n) downto 0 + for sum of-type number = (aref b i) + do (loop for j of-type fixnum from (1+ i) below n + do (decf sum (* (aref lu i j) (aref b i)))) + do (setf (aref b i) (/ sum (aref lu i i)))) + ) + b) + + +;;;--------------------------------------------------------------------------- +;;; solve + +(defmethod solve ((a array) (b vector) &optional (result (copy-seq b))) + (multiple-value-bind (lu permutations) + (lu-decompose a) + (lu-back-substitution lu permutations result))) + +(defmethod solve ((a matrix) (b vector) &optional (result (copy-seq b))) + (solve (matrix-data a) b result)) + + +(defmethod solve/double-float ((a array) (b vector) &optional (result (copy-seq b))) + (multiple-value-bind (lu permutations) + (lu-decompose/double-float a) + (lu-back-substitution lu permutations result))) + +(defmethod solve/double-float ((a matrix) (b vector) &optional (result (copy-seq b))) + (solve/double-float (matrix-data a) b result)) + + +;;;--------------------------------------------------------------------------- +;;; det + +(defmethod det ((a array)) + (assert (square-matrix-p a)) + (multiple-value-bind (lu permutations parity) + (lu-decompose a) + (declare (ignore permutations)) + (dotimes (i (array-dimension lu 0) parity) + (setf parity (* parity (aref lu i i)))))) + +(defmethod det ((a matrix)) + (det (matrix-data a))) + + +(defmethod det/double-float ((a array)) + (assert (square-matrix-p a)) + (multiple-value-bind (lu permutations parity) + (lu-decompose a) + (declare (type (array double-float (cl:* cl:*)) lu) + (ignore permutations) + (type double-float parity)) + (dotimes (i (array-dimension lu 0) parity) + (setf parity (* parity (aref lu i i)))))) + +(defmethod det/double-float ((a matrix)) + (det/double-float (matrix-data a))) + + +;;;=========================================================================== +;;; Test + +(defvar A #2A((1 2 3) + (2 -1 1) + (3 4 -1))) + +(defvar LU-A #2A((1 2 3) + (2 -5 -5) + (3 0.4 -8))) + +;;; end of file -- lu-decomposition.lisp -- Added: trunk/common-math/numerics/linear-algebra/matrix-pkg.lisp ============================================================================== --- (empty file) +++ trunk/common-math/numerics/linear-algebra/matrix-pkg.lisp Wed Aug 16 17:07:41 2006 @@ -0,0 +1,19 @@ +;;; -*- Mode: Lisp -*- + +(defpackage "CL.MATH.POLYNOMIALS" (:use "CL.MATH" "COMMON-LISP") + (:nicknames "POLY") + + (:shadow cl:zerop cl:expt cl:gcd) + + (:shadowing-import-from "CL.MATH" "+" "-" "*" "/" "=" ">" "<" ">=" "<=") + + (:export "=" "+" "-" "*" "/" ">" "<" ">=" "<=") + (:export "ZEROP" "EXPT" "GCD") + + (:export "MAKE-POLYNOMIAL" + "POLYNOMIAL" + "UNIVARIATE-POLYNOMIAL" + ) + ) + +;;; end of file -- polymonials-pkg.lisp -- Added: trunk/common-math/numerics/linear-algebra/matrix.lisp ============================================================================== --- (empty file) +++ trunk/common-math/numerics/linear-algebra/matrix.lisp Wed Aug 16 17:07:41 2006 @@ -0,0 +1,950 @@ +;;;; -*- Mode: Lisp -*- + +;;;; matrix.lisp -- + +(in-package "CL.MATH.LINEAR-ALGEBRA") + +(deftype matrix-array (&optional (type 'number) (m 'cl:*) (n 'cl:*)) + `(array ,type (,m ,n))) + + +(defparameter *default-matrix-element-type* 'double-float) + +(defclass matrix () + ((data :reader matrix-data :initarg :data)) + (:default-initargs :data (make-array '(0 0) :element-type 'number))) + +(defgeneric matrix-p (x + #|&optional + (element-type *default-matrix-element-type*) + |# + ) + (:method ((x matrix)) t) + (:method ((x t)) nil) + (:method ((x vector)) t) ; Note that this is technically correct, but + ; it may cause problems. + (:method ((x array)) (cl:<= (array-rank x) 2)) ; Same here. + ) + + +(defvar *check-element-type* nil) + + +;;;--------------------------------------------------------------------------- +;;; Protocol. + +(defgeneric matrix-array-p (x &key check-element-type) + (:method ((x array) &key (check-element-type *check-element-type*)) + (and (cl:<= (array-rank x) 2) + (or (not check-element-type) + (nth-value 0 (subtypep (array-element-type x) 'number))))) + + (:method ((x t) &key check-element-type) + (declare (ignore check-element-type)) + nil)) + + + +(defgeneric matrix-element-type (m)) + +(defgeneric matrix-row-number (m) + (:method ((x matrix)) (array-dimension (matrix-data x) 0)) + (:method ((x vector)) 1) + (:method ((x array)) + (assert (matrix-array-p x)) + (array-dimension x 0))) + + +(defgeneric matrix-column-number (m) + (:method ((x matrix)) (array-dimension (matrix-data x) 1)) + (:method ((x vector)) (length x)) + (:method ((x array)) + (assert (matrix-array-p x)) + (array-dimension x 1))) + + +(defgeneric matrix-dimensions (m) + (:method ((x matrix)) (array-dimensions (matrix-data x))) + (:method ((x array)) + (assert (matrix-array-p x)) + (array-dimensions x))) + + +(defgeneric matrix-dimension (m axis) + (:method ((m matrix) axis) + (declare (type fixnum axis)) + (array-dimension (matrix-data m) axis)) + (:method ((m array) axis) + (declare (type fixnum axis)) + (array-dimension m axis))) + + +(defmethod matrix-dimension :before ((m array) axis) + (assert (matrix-array-p m))) + + + +(defmethod matrix-element-type ((m matrix)) (array-element-type (matrix-data m))) + +(defmethod matrix-element-type ((m array)) (array-element-type m)) + +(defmethod matrix-element-type :before ((m array)) + (assert (matrix-array-p m))) + + +(defgeneric shape-equal-p (m1 m2) + (:documentation "Returns non-NIL when M1 and M2 have the same `shape'.")) + + +(defmethod shape-equal-p ((m1 matrix) (m2 matrix)) + (with-slots ((m1d data)) + m1 + (with-slots ((m2d data)) + m2 + (and (cl:= (array-dimension m1d 0) + (array-dimension m2d 0)) + (cl:= (array-dimension m1d 1) + (array-dimension m2d 1)))))) + + +(defmethod shape-equal-p ((m1 array) (m2 matrix)) + (assert (matrix-array-p m1 :check-element-type nil)) + (let ((m1d m1)) + (with-slots ((m2d data)) + m2 + (and (cl:= (array-dimension m1d 0) + (array-dimension m2d 0)) + (cl:= (array-dimension m1d 1) + (array-dimension m2d 1)))))) + + +(defmethod shape-equal-p ((m2 matrix) (m1 array)) + (shape-equal-p m1 m2)) + + +(defmethod shape-equal-p ((m1 array) (m2 array)) + (assert (matrix-array-p m1 :check-element-type nil)) + (assert (matrix-array-p m2 :check-element-type nil)) + (let ((m1d m1) + (m2d m2) + ) + (and (cl:= (array-dimension m1d 0) + (array-dimension m2d 0)) + (cl:= (array-dimension m1d 1) + (array-dimension m2d 1))))) + + +(defmethod shape-equal-p ((m1 vector) (m2 vector)) + (cl:= (cl:length m1) (cl:length m2))) + + +(defmethod shape-equal-p ((m1 column) (m2 vector)) + (cl:= 1 (vector-length m1) (cl:length m2))) + + +(defmethod shape-equal-p ((m2 vector) (m1 column)) + (shape-equal-p m1 m2)) + + +;;;--------------------------------------------------------------------------- +;;; Creation. + +(defun make-matrix (n m + &rest array-keys + &key + (element-type *default-matrix-element-type*) + (initial-element (coerce 0 element-type) ie-supplied-p) + (initial-contents () ic-supplied-p) + &allow-other-keys) + (when (and ie-supplied-p ic-supplied-p) + (error "The MAKE-MATRIX :INITIAL-ELEMENT and :INITIAL-CONTENTS keywords may not be simultaneously supplied.")) + + (let ((data (cond ((and ic-supplied-p + (matrix-array-p initial-contents)) + initial-contents) + (ic-supplied-p + (apply #'make-array (list n m) + (append array-keys + (list :element-type element-type)))) + (t + (apply #'make-array (list n m) + (append array-keys + (list :element-type element-type + :initial-element initial-element))) + ))) + ) + (make-instance 'matrix :data data))) + + +(defun make-identity-matrix (n &rest array-keys &key &allow-other-keys) + (when (getf array-keys :initial-contents) + (warn "MAKE-IDENTITY-MATRIX passed an :INITIAL-CONTENTS argument.") + (remf array-keys :initial-contents)) + (when (getf array-keys :initial-element) + (warn "MAKE-IDENTITY-MATRIX passed an :INITIAL-ELEMENT argument.") + (remf array-keys :initial-element)) + (let* ((id (apply #'make-matrix n n array-keys)) + (data (matrix-data id)) + (one (coerce 1 (array-element-type data))) + ) + (dotimes (i n id) + (setf (aref data i i) one)))) + + +(defun eye (n &rest array-keys &key &allow-other-keys) + (apply #'make-identity-matrix n array-keys)) + + +(defun make-zero-matrix (n m + &rest array-keys + &key (element-type *default-matrix-element-type*) + &allow-other-keys) + (when (getf array-keys :initial-element) + (warn "MAKE-ZERO-MATRIX passed and :INITIAL-ELEMENT argument.") + (remf array-keys :initial-element)) + (apply #'make-matrix n m (append array-keys + (list :initial-element (coerce 0 element-type) + :element-type element-type)))) + + +(defun zeros (n m &rest array-keys &key &allow-other-keys) + (apply #'make-zero-matrix n m array-keys)) + + +;;;--------------------------------------------------------------------------- +;;; copy-matrix -- +;;; Uses the code from K Pitman and B. Margolin appeared on CLL a long +;;; long time ago. +;;; Note that at least on LW this code has some of the usual problems +;;; with DOUBLE-FLOATS. +#| +(defun copy-array (array) + (adjust-array (make-array (array-dimensions array) + :displaced-to array + :element-type (array-element-type array)) + (array-dimensions array) + :displaced-to nil)) +|# + + +(defmethod copy-matrix ((m matrix)) + (let* ((array (matrix-data m)) + (result-data + (adjust-array (make-array (array-dimensions array) + :displaced-to array + :element-type (array-element-type array)) + (array-dimensions array) + :displaced-to nil)) + ) + (make-instance 'matrix :data result-data))) + + +(defmethod copy-matrix ((m array)) + (check-type m matrix-array) + (adjust-array (make-array (array-dimensions m) + :displaced-to m + :element-type (array-element-type m)) + (array-dimensions m) + :displaced-to nil)) + + +;;;--------------------------------------------------------------------------- +;;; Equality. + +(defmethod =%2 ((x array) (y array)) + (and (matrix-array-p x) + (matrix-array-p y) + ;; (cl:= (matrix-row-number x) (matrix-row-number y)) + ;; (cl:= (matrix-column-number x) (matrix-column-number y)) + (cl:= (array-dimension x 0) (array-dimension y 0)) + (cl:= (array-dimension x 1) (array-dimension y 1)) + (loop for i from 0 below (array-total-size x) + always (=%2 (row-major-aref x i) (row-major-aref y i))))) + + +(defmethod =%2 ((x matrix) (y matrix)) + (or (eq x y) + (=%2 (matrix-data x) (matrix-data y)))) + + +(defmethod =%2 ((x matrix) (y array)) + (=%2 (matrix-data x) y)) + + +(defmethod =%2 ((y array) (x matrix)) + (=%2 (matrix-data x) y)) + + +;;;--------------------------------------------------------------------------- +;;; zerop + +(defmethod zerop ((x array)) + (and (matrix-array-p x) + (loop for i from 0 below (array-total-size x) + always (zerop (row-major-aref x i))))) + + +(defmethod zerop ((x matrix)) + (zerop (matrix-data x))) + + +;;;--------------------------------------------------------------------------- +;;; Addition. + +(defmethod +%2 :before ((x matrix) (y matrix) + &optional + (r nil r-supplied-p)) + (assert (shape-equal-p x y)) + (when r-supplied-p + (assert (matrix-p x)) + (assert (shape-equal-p x r)))) + + +(defmethod +%2 :before ((x matrix) (y array) + &optional + (r nil r-supplied-p)) + (assert (shape-equal-p x y)) + (when r-supplied-p + (assert (matrix-p x)) + (assert (shape-equal-p x r)))) + + +(defmethod +%2 :before ((y array) (x matrix) + &optional + (r nil r-supplied-p)) + (assert (shape-equal-p x y)) + (when r-supplied-p + (assert (matrix-p x)) + (assert (shape-equal-p x r)))) + + +(defmethod +%2 ((x number) (y matrix) &optional (r (copy-matrix y) r-supplied-p)) + (when r-supplied-p (assert (shape-equal-p y r))) + (with-slots (data) y + (with-slots ((result data)) r + (dotimes (i (array-total-size data) r) + (setf (row-major-aref result i) (+%2 x (row-major-aref data i)))) + ))) + + +(defmethod +%2 ((y matrix) (x number) &optional (r (copy-matrix y) r-supplied-p)) + (when r-supplied-p (assert (shape-equal-p y r))) + (with-slots (data) y + (with-slots ((result data)) r + (dotimes (i (array-total-size data) r) + (setf (row-major-aref result i) (+%2 x (row-major-aref data i)))) + ))) + + + +#+with-slots-slow +(defmethod +%2 ((x matrix) (y matrix) &optional (r (copy-matrix y) r-supplied-p)) + (with-slots ((m1 data)) x + (with-slots ((m2 data)) y + (with-slots ((result data)) r + (dotimes (i (array-total-size m1)) + (setf (row-major-aref result i) + (cl:+ (row-major-aref m1 i) (row-major-aref m2 i)))) + ))) + r + ) + + +(defmethod +%2 ((x matrix) (y matrix) &optional (r (copy-matrix y) r-supplied-p)) + (declare (ignore r-supplied-p)) + (let ((m1 (matrix-data x)) + (m2 (matrix-data y)) + (result (matrix-data r)) + ) + (dotimes (i (array-total-size m1) r) + (setf (row-major-aref result i) + (cl:+ (row-major-aref m1 i) (row-major-aref m2 i)))) + )) + + +#+with-slots-slow +(defmethod +%2 ((x array) (y matrix) &optional (r (copy-matrix y) r-supplied-p)) + (let ((m1 x)) + (with-slots ((m2 data)) y + (with-slots ((result data)) r + (dotimes (i (array-total-size m1)) + (setf (row-major-aref result i) + (cl:+ (row-major-aref m1 i) (row-major-aref m2 i)))) + ))) + r) + + + +(defmethod +%2 ((x array) (y matrix) &optional (r (copy-matrix y) r-supplied-p)) + (declare (ignore r-supplied-p)) + (let ((m1 x) + (m2 (matrix-data y)) + (result (matrix-data r)) + ) + (dotimes (i (array-total-size m1) r) + (setf (row-major-aref result i) + (cl:+ (row-major-aref m1 i) (row-major-aref m2 i)))) + )) + + +#+with-slots-slow +(defmethod +%2 ((y matrix) (x array) &optional (r (copy-matrix y) r-supplied-p)) + (let ((m1 x)) + (with-slots ((m2 data)) y + (with-slots ((result data)) r + (dotimes (i (array-total-size m1)) + (setf (row-major-aref result i) + (cl:+ (row-major-aref m1 i) (row-major-aref m2 i)))) + ))) + r) + + + +(defmethod +%2 ((y matrix) (x array) &optional (r (copy-matrix y) r-supplied-p)) + (declare (ignore r-supplied-p)) + (let ((m1 x) + (m2 (matrix-data y)) + (result (matrix-data r)) + ) + (dotimes (i (array-total-size m1) r) + (setf (row-major-aref result i) + (cl:+ (row-major-aref m1 i) (row-major-aref m2 i)))) + )) + + +;;;--------------------------------------------------------------------------- +;;; Subtraction. + + +(defmethod -%2 :before ((x matrix) (y matrix) + &optional + (r nil r-supplied-p)) + (assert (shape-equal-p x y)) + (when r-supplied-p + (assert (matrix-p x)) + (assert (shape-equal-p x r)))) + + +(defmethod -%2 :before ((x matrix) (y array) + &optional + (r nil r-supplied-p)) + (assert (shape-equal-p x y)) + (when r-supplied-p + (assert (matrix-p x)) + (assert (shape-equal-p x r)))) + + +(defmethod -%2 :before ((y array) (x matrix) + &optional + (r nil r-supplied-p)) + (assert (shape-equal-p x y)) + (when r-supplied-p + (assert (matrix-p x)) + (assert (shape-equal-p x r)))) + + +(defmethod -%2 ((x matrix) (y matrix) &optional (r (copy-matrix y) r-supplied-p)) + (declare (ignore r-supplied-p)) + (with-slots ((m1 data)) x + (with-slots ((m2 data)) y + (with-slots ((result data)) r + (dotimes (i (array-total-size m1)) + (setf (row-major-aref result i) + (cl:- (row-major-aref m1 i) (row-major-aref m2 i)))) + ))) + r + ) + +(defmethod -%2 ((x array) (y matrix) &optional (r (copy-matrix y) r-supplied-p)) + (declare (ignore r-supplied-p)) + (let ((m1 x)) + (with-slots ((m2 data)) y + (with-slots ((result data)) r + (dotimes (i (array-total-size m1)) + (setf (row-major-aref result i) + (cl:- (row-major-aref m1 i) (row-major-aref m2 i)))) + ))) + r) + +(defmethod -%2 ((y matrix) (x array) &optional (r (copy-matrix y) r-supplied-p)) + (declare (ignore r-supplied-p)) + (with-slots ((m1 data)) y + (let ((m2 x)) + (with-slots ((result data)) r + (dotimes (i (array-total-size m1)) + (setf (row-major-aref result i) + (cl:- (row-major-aref m1 i) (row-major-aref m2 i)))) + ))) + r) + + +;;;--------------------------------------------------------------------------- +;;; Multiplication. + +(defgeneric make-*-result-matrix (x y &optional element-type) + (:documentation + "Creates a matrix serve as result of a multiplication.") + + (:method ((x matrix) (y matrix) + &optional (element-type *default-matrix-element-type*)) + (make-matrix (matrix-row-number x) + (matrix-column-number y) + :element-type element-type)) + + (:method ((x array) (y matrix) + &optional (element-type *default-matrix-element-type*)) + (make-matrix (matrix-row-number x) + (matrix-column-number y) + :element-type element-type)) + + (:method ((x matrix) (y array) + &optional (element-type *default-matrix-element-type*)) + (make-matrix (matrix-row-number x) + (matrix-column-number y) + :element-type element-type)) + + (:method ((x vector) (y matrix) + &optional (element-type *default-matrix-element-type*)) + (make-matrix 1 + (matrix-column-number y) + :element-type element-type)) + + (:method ((x matrix) (y vector) + &optional (element-type *default-matrix-element-type*)) + (make-matrix (matrix-row-number x) + 1 + :element-type element-type))) + + + +(defgeneric conforming-*-dimensions-p (x y r) + (:documentation + "Checks whether the dimensions of R are correct to hold XxY.") + + (:method ((x t) (y t) (r t)) + (error "Result object is not a matrix: ~S." r)) + + ;; 1. Matrix results. + + ;; 1.1. Number methods. + + (:method ((x number) (y matrix) (r matrix)) + (shape-equal-p y r)) + + (:method ((x number) (y array) (r matrix)) + (shape-equal-p y r)) + + (:method ((y matrix) (x number) (r matrix)) + (conforming-*-dimensions-p x y r)) + + (:method ((y array) (x number) (r matrix)) + (conforming-*-dimensions-p x y r)) + + ;; 1.2. Matrix and Array arguments. + + (:method ((x array) (y array) (r matrix)) + (and (matrix-array-p x) + (matrix-array-p y) + (cl:= (matrix-row-number x) (matrix-column-number y)) + (cl:= (matrix-row-number x) (matrix-row-number r)) + (cl:= (matrix-column-number y) (matrix-column-number r)))) + + (:method ((x vector) (y array) (r matrix)) + (and (matrix-array-p y) + (cl:= (matrix-row-number x) (matrix-column-number y)) + (cl:= (matrix-row-number x) (matrix-row-number r)) + (cl:= (matrix-column-number y) (matrix-column-number r)))) + + + (:method ((x array) (y vector) (r matrix)) + (and (matrix-array-p x) + (cl:= (matrix-row-number x) (matrix-column-number y)) + (cl:= (matrix-row-number x) (matrix-row-number r)) + (cl:= (matrix-column-number y) (matrix-column-number r)))) + + (:method ((x matrix) (y matrix) (r matrix)) + (conforming-*-dimensions-p (matrix-data x) (matrix-data y) r)) + + (:method ((x array) (y matrix) (r matrix)) + (conforming-*-dimensions-p x (matrix-data y) r)) + + (:method ((x matrix) (y array) (r matrix)) + (conforming-*-dimensions-p (matrix-data x) y r)) + + + ;; 2. Array results. + + ;; 2.1. Number methods. + + (:method ((x number) (y matrix) (r array)) + (shape-equal-p y r)) + + (:method ((x number) (y array) (r array)) + (shape-equal-p y r)) + + (:method ((y matrix) (x number) (r array)) + (conforming-*-dimensions-p x y r)) + + (:method ((y array) (x number) (r array)) + (conforming-*-dimensions-p x y r)) + + ;; 2.2. Matrix and Array arguments. + + (:method ((x array) (y array) (r array)) + (and (matrix-array-p x) + (matrix-array-p y) + (matrix-array-p r) + (cl:= (matrix-row-number x) (matrix-column-number y)) + (cl:= (matrix-row-number x) (matrix-row-number r)) + (cl:= (matrix-column-number y) (matrix-column-number r)))) + + (:method ((x vector) (y array) (r array)) + (and (matrix-array-p y) + (cl:= (matrix-row-number x) (matrix-column-number y)) + (cl:= (matrix-row-number x) (matrix-row-number r)) + (cl:= (matrix-column-number y) (matrix-column-number r)))) + + (:method ((x array) (y vector) (r array)) + (and (matrix-array-p x) + (cl:= (matrix-row-number x) (matrix-column-number y)) + (cl:= (matrix-row-number x) (matrix-row-number r)) + (cl:= (matrix-column-number y) (matrix-column-number r)))) + + (:method ((x matrix) (y matrix) (r array)) + (conforming-*-dimensions-p (matrix-data x) (matrix-data y) r)) + + (:method ((x array) (y matrix) (r array)) + (conforming-*-dimensions-p x (matrix-data y) r)) + + (:method ((x matrix) (y array) (r array)) + (conforming-*-dimensions-p (matrix-data x) y r)) + + ;; 3. Other methods. + + (:method ((x vector) (y vector) (r t)) + ;; Maybe a warning here is needed. + (cl:= (length x) (length y))) + ) + + +;;; *%2 Methods. + +(defmethod *%2 ((x number) (y matrix) &optional (r (copy-matrix y) r-supplied-p)) + (when r-supplied-p (assert (conforming-*-dimensions-p x y r))) + (with-slots ((m data)) y + (with-slots ((result data)) r + (dotimes (i (array-total-size m) r) + (setf (row-major-aref result i) (*%2 x (row-major-aref m i)))) + ))) + + +(defmethod *%2 ((y matrix) (x number) &optional (r (copy-matrix y) r-supplied-p)) + (declare (ignore r-supplied-p)) + (*%2 x y r)) + + +(defmethod *%2 ((x matrix) (y matrix) + &optional (r (make-*-result-matrix x y) r-supplied-p)) + (when r-supplied-p (assert (conforming-*-dimensions-p x y r))) + ;; X is a NxK matrix + ;; Y is a KxM matrix + ;; R is a NxM matrix + (let* ((m1 (matrix-data x)) + (m2 (matrix-data y)) + (result (matrix-data r)) + (n (array-dimension m1 0)) + (k (array-dimension m1 1)) + (m (array-dimension m2 0)) + (z (coerce 0 (array-element-type result))) + ) + (declare (type fixnum n k m)) + (dotimes (i n r) + (dotimes (j m) + (let ((temp z)) + (dotimes (l k) + (setf (aref result i j) + (setf temp (+%2 temp (*%2 (aref m1 i l) (aref m2 l j))))) + )) + )))) + + +(defmethod *%2 ((x array) (y matrix) + &optional (r (make-*-result-matrix x y) r-supplied-p)) + (when r-supplied-p (assert (conforming-*-dimensions-p x y r))) + ;; X is a NxK matrix + ;; Y is a KxM matrix + ;; R is a NxM matrix + (let* ((m1 x) + (m2 (matrix-data y)) + (result (matrix-data r)) + (n (array-dimension m1 0)) + (k (array-dimension m1 1)) + (m (array-dimension m2 0)) + (z (coerce 0 (array-element-type result))) + ) + (declare (type fixnum n k m)) + (dotimes (i n r) + (dotimes (j m) + (let ((temp z)) + (dotimes (l k) + (setf (aref result i j) + (incf temp (* (aref m1 i l) (aref m2 l j)))) + )) + )))) + + +(defmethod *%2 ((x matrix) (y array) + &optional (r (make-*-result-matrix x y) r-supplied-p)) + (when r-supplied-p (assert (conforming-*-dimensions-p x y r))) + ;; X is a NxK matrix + ;; Y is a KxM matrix + ;; R is a NxM matrix + (let* ((m1 (matrix-data x)) + (m2 y) + (result (matrix-data r)) + (n (array-dimension m1 0)) + (k (array-dimension m1 1)) + (m (array-dimension m2 1)) + (z (coerce 0 (array-element-type result))) + ) + (declare (type fixnum n k m)) + (dotimes (i n r) + (dotimes (j m) + (let ((temp z)) + (dotimes (l k) + (setf (aref result i j) + (incf temp (* (aref m1 i l) (aref m2 l j)))) + )) + )))) + + +(defmethod *%2 ((x vector) (y matrix) + &optional (r (make-*-result-matrix x y) r-supplied-p)) + (when r-supplied-p (assert (conforming-*-dimensions-p x y r))) + (let* ((m1 x) + (m2 (matrix-data y)) + (result (matrix-data r)) + (n 1) + (k (length m1)) + (m (array-dimension m2 1)) + (z (coerce 0 (array-element-type result))) + ) + (declare (type fixnum k m) (type (integer 1 1) n)) + (dotimes (i n r) + (dotimes (j m) + (let ((temp z)) + (dotimes (l k) + (setf (aref result i j) + (incf temp (* (aref m1 l) (aref m2 l j)))) + )) + )))) + + +(defmethod *%2 ((x matrix) (y vector) + &optional (r (make-*-result-matrix x y) r-supplied-p)) + (when r-supplied-p (assert (conforming-*-dimensions-p x y r))) + (let* ((m1 y ) + (m2 (matrix-data x)) + (result (matrix-data r)) + (n (array-dimension m1 0)) + (k (length y)) + (m k) + (z (coerce 0 (array-element-type result))) + ) + (declare (type fixnum n k m)) + (dotimes (i n r) + (dotimes (j m) + (let ((temp z)) + (dotimes (l k) + (setf (aref result i j) + (incf temp (* (aref m1 l) (aref m2 l j)))) + )) + )))) + + +(defmethod *%2 ((x number) (y array) &optional (r (copy-matrix y) r-supplied-p)) + (when r-supplied-p (assert (conforming-*-dimensions-p x y r))) + (let ((m y) + (result r) + ) + (dotimes (i (array-total-size m) r) + (setf (row-major-aref result i) (*%2 x (row-major-aref m i)))) + )) + + +(defmethod *%2 ((y array) (x number) &optional (r (copy-matrix y) r-supplied-p)) + (*%2 x y r)) + + +;;; The next one breaks the return type convention. +#| Defined in 'vector.lisp'. +(defmethod *%2 ((x vector) (y vector) &optional r) + (declare (ignore r)) + (assert (conforming-*-dimensions-p x y nil)) + (let ((result 0)) + (dotimes (i (length x) result) + (setf result (+%2 result (* (aref x i) (aref y i))))))) +|# + + +;;;--------------------------------------------------------------------------- +;;; Division. +;;; Only the simple form of division is implemented here. +;;; The general form requires the inverse, which requires extra machinery + +(defmethod /%2 ((x number) (y matrix) &optional (r (copy-matrix y) r-supplied-p)) + (declare (ignore r-supplied-p)) + (*%2 (/ x) y r)) + + +(defmethod /%2 ((y matrix) (x number) &optional (r (copy-matrix y) r-supplied-p)) + (declare (ignore r-supplied-p)) + (*%2 (/ x) y r)) + + +;;;--------------------------------------------------------------------------- +;;; Exponentiation. +;;; ... missing ... + + +;;;--------------------------------------------------------------------------- +;;; Element-wise operations. + +(defmethod .*%2 ((x number) (y matrix) + &optional (r (copy-matrix y) r-supplied-p)) + (when r-supplied-p (assert (shape-equal-p y r))) + (with-slots (data) y + (with-slots ((result data)) r + (dotimes (i (array-total-size data) r) + (setf (row-major-aref result i) (*%2 x (row-major-aref data i)))) + ))) + + +(defmethod .*%2 ((y matrix) (x number) &optional (r (copy-matrix y))) + (.*%2 x y r)) + + +(defmethod .*%2 ((x number) (y array) + &optional (r (copy-matrix y) r-supplied-p)) + (assert (matrix-array-p y)) + (when r-supplied-p (assert (shape-equal-p y r))) + (let* ((data y) + (result r) + ) + (dotimes (i (array-total-size data) r) + (setf (row-major-aref result i) (*%2 x (row-major-aref data i)))) + )) + + +(defmethod .*%2 ((y array) (x number) &optional (r (copy-matrix y))) + (.*%2 x y r)) + + + +(defmethod .*%2 ((x array) (y array) &optional (r (copy-matrix y) r-supplied-p)) + (when r-supplied-p (assert (shape-equal-p y r))) + (assert (shape-equal-p x y)) + + (dotimes (i (array-total-size y) r) + (setf (row-major-aref r i) + (*%2 (row-major-aref x i) (row-major-aref y i))))) + + +;;;--------------------------------------------------------------------------- +;;; Transpose + +(defgeneric conforming-transpose-dimensions-p (x r) + (:documentation + "Checks whether the dimensions of R are correct to hold X'.") + + (:method ((x t) (r t)) + (error "Result object is not a matrix: ~S." r)) + + ;; 1. Matrix results. + + (:method ((x array) (r matrix)) + (and (matrix-array-p x) + (cl:= (matrix-row-number x) (matrix-column-number r)) + (cl:= (matrix-column-number x) (matrix-row-number r)))) + + (:method ((x matrix) (r matrix)) + (and (cl:= (matrix-row-number x) (matrix-column-number r)) + (cl:= (matrix-column-number x) (matrix-row-number r)))) + + + ;; 2. Array results. + + (:method ((x array) (r array)) + (and (matrix-array-p x) + (matrix-array-p r) + (cl:= (matrix-row-number x) (matrix-column-number r)) + (cl:= (matrix-column-number x) (matrix-row-number r)))) + + (:method ((x matrix) (r array)) + (and (matrix-array-p r) + (cl:= (matrix-row-number x) (matrix-column-number r)) + (cl:= (matrix-column-number x) (matrix-row-number r)))) + ) + + + +(defmethod transpose ((m matrix) + &optional + (r (make-matrix (matrix-column-number m) + (matrix-row-number m) + :element-type + (matrix-element-type m)) + r-supplied-p)) + (when r-supplied-p (assert (conforming-transpose-dimensions-p m r))) + (let ((rows (matrix-row-number m)) + (cols (matrix-column-number m)) + (m (matrix-data m)) + (result (matrix-data r)) + ) + (dotimes (i rows r) + (dotimes (j cols) + (setf (aref result j i) (aref m i j)))))) + + +(defmethod transpose ((m array) + &optional + (r (make-array (array-dimensions m) + :element-type + (array-element-type m) + ) + r-supplied-p)) + (when r-supplied-p (assert (conforming-transpose-dimensions-p m r))) + (let ((rows (matrix-row-number m)) + (cols (matrix-column-number m)) + ) + (dotimes (i rows r) + (dotimes (j cols) + (setf (aref r j i) (aref m i j)))))) + + +(defmethod transpose :before ((m array) &optional r) + (declare (ignore r)) + (assert (matrix-array-p m))) + + +;;;--------------------------------------------------------------------------- +;;; Utilities. + +#+simple +(defmethod print-object ((m matrix) stream) + (print-unreadable-object (m stream :identity t :type t) + (format stream "~S" (matrix-data m)))) + + +(defmethod print-object ((m matrix) stream) + (let ((d (matrix-data m))) + (format stream "#M(~S ~S ~S)" + (array-element-type d) + (array-dimensions d) + d))) + +;;; Missing: +;;; Define #M read macro. + + +;;;; end of file -- matrix.lisp -- Added: trunk/common-math/numerics/linear-algebra/vector-matrix-conditions.lisp ============================================================================== --- (empty file) +++ trunk/common-math/numerics/linear-algebra/vector-matrix-conditions.lisp Wed Aug 16 17:07:41 2006 @@ -0,0 +1,27 @@ +;;;; -*- Mode: Lisp -*- + +;;;; vector-matrix-conditions.lisp -- + +(in-package "CL.MATH.LINEAR-ALGEBRA") + +;;;--------------------------------------------------------------------------- +;;; Conditions relative to the vectors and matrix subsystem.. + +(define-condition incompatible-shapes (error) + ((a :accessor incompatible-shape-a + :initarg :a) + (b :accessor incompatible-shape-b + :initarg :b) + ) + (:documentation "The Incompatible Shapes Error. +This error is signalled by functions that get arguments -- mainly +matrices and/or vectors -- whose 'shape' is not compatible with +respect to the semantics of the operation.") + + (:report (lambda (isc stream) + (declare (ignore isc)) + (format stream "The shapes passed have incompatible dimensions."))) + ) + + +;;;; end of file -- vector-matrix-conditions.lisp -- Added: trunk/common-math/numerics/linear-algebra/vector-matrix.lisp ============================================================================== --- (empty file) +++ trunk/common-math/numerics/linear-algebra/vector-matrix.lisp Wed Aug 16 17:07:41 2006 @@ -0,0 +1,613 @@ +;;;; -*- Mode: Lisp -*- + +;;;; vector-matrix.lisp -- +;;;; Protocols and functions common to and dependent on vector and matrix +;;;; definitions. + +(in-package "CL.MATH.LINEAR-ALGEBRA") + +;;;;=========================================================================== +;;;; Protocol. + +;;; join and stack +;;; Note: no 'result' machinery in place yet. + +(defgeneric join%2 (d1 d2)) + +(defgeneric stack%2 (d1 d2)) + + +;;; slice + +(defgeneric slice (e shape-designator)) + + +;;; range + +(defstruct (range (:constructor %range (start &optional end step))) + (start 0 :type (or fixnum (member * cl:*))) + (end nil :type (or null fixnum (member * cl:*))) + (step 1 :type fixnum)) + + +(defstruct (decreasing-range (:include range))) + +(defstruct (increasing-range (:include range))) + + +(defun range (start &optional end (step 1)) + (if end + (if (<= start end) + (make-increasing-range :start start :end end :step step) + (make-decreasing-range :start start :end end :step step)) + (if (plusp step) + (make-increasing-range :start start :end end :step step) + (make-decreasing-range :start start :end end :step step)) + )) + + + +;;; ref + +(defgeneric ref (e index &rest indexes)) + +(defgeneric ref%2 (e i j)) + + +;;;--------------------------------------------------------------------------- +;;; Implementation. + +;;; join%2 -- + +(defmethod join%2 ((d1 vector) (d2 vector)) + (concatenate 'vector d1 d2)) + + +(defmethod join%2 ((d1 string) (d2 string)) + (concatenate 'string d1 d2)) + + +(defmethod join%2 ((d1 list) (d2 list)) + (concatenate 'list d1 d2)) + + +(defmethod join%2 ((d1 vector) (d2 sequence)) + (concatenate 'vector d1 d2)) + + +(defmethod join%2 ((d1 sequence) (d2 vector)) + (concatenate (type-of d1) d1 d2)) + + +(defmethod join%2 ((d1 column) (d2 column)) + (let ((initial-contents + (map 'list #'vector (column-vector d1) (column-vector d2)))) + (declare (dynamic-extent initial-contents)) + (make-array (list (length (column-vector d1)) 2) + :initial-contents initial-contents))) + + +(defmethod join%2 ((d1 vector) (d2 column)) + (unless (cl:= (length (column-vector d2)) 1) + (error "Trying to JOIN a vector (1x~D) to a column (~Dx1)." + (length d1) + (length (column-vector d2)))) + (concatenate 'vector d1 (column-vector d2))) + + +(defmethod join%2 ((d1 column) (d2 vector)) + (unless (cl:= (length (column-vector d1)) 1) + (error "Trying to JOIN a column (~Dx1) to a vector (1x~D)." + (length (column-vector d1)) + (length d2))) + (concatenate 'vector (column-vector d1) d2)) + + +(defmethod join%2 ((d1 column) (d2 array)) + (let ((nr (array-dimension d2 0))) + (assert (cl:= (length (column-vector d1)) nr)) + (let* ((nc (1+ (array-dimension d2 1))) + (result (make-array (list nr nc) + :initial-element (aref d2 0 0) + :element-type (array-element-type d2))) + (col (column-vector d1)) + ) + (dotimes (i nr result) + (setf (aref result i 0) (aref col i)) + (loop for j from 1 below nc + do (setf (aref result i j) (aref d2 i (1- j))))) + ))) + + +(defmethod join%2 ((d1 array) (d2 column)) + (let ((nr (array-dimension d1 0))) + (assert (cl:= (length (column-vector d2)) nr)) + (let* ((nc (1+ (array-dimension d1 1))) + (nc-1 (1- nc)) + (result (make-array (list nr nc) + :initial-element (aref d1 0 0) + :element-type (array-element-type d1))) + (col (column-vector d2)) + ) + + (dotimes (i nr result) + (setf (aref result i nc-1) (aref col i)) + (dotimes (j nc-1) + (setf (aref result i j) (aref d1 i j)))) + ))) + + +(defmethod join%2 ((d1 array) (d2 array)) + (let ((nr1 (array-dimension d1 0)) + (nr2 (array-dimension d2 0)) + ) + (assert (cl:= nr1 nr2)) + (let* ((nc1 (array-dimension d1 1)) + (nc2 (array-dimension d2 1)) + (result (make-array (list nr1 (+ nc1 nc2)) + :initial-element (aref d1 0 0) + :element-type (array-element-type d1))) + ) + + (dotimes (i nr1) + (dotimes (j nc1) + (setf (aref result i j) (aref d1 i j)))) + + (dotimes (i nr2 result) + (dotimes (j nc2) + (setf (aref result i (+ j nc1)) (aref d2 i j)))) + ))) + + +(defmethod join%2 ((d1 matrix) (d2 matrix)) + (let ((nr1 (matrix-row-number d1)) + (nr2 (matrix-row-number d2)) + ) + (assert (cl:= nr1 nr2)) + (let* ((nc1 (matrix-column-number d1)) + (nc2 (matrix-column-number d2)) + (data1 (matrix-data d1)) + (data2 (matrix-data d2)) + (result (make-array (list nr1 (+ nc1 nc2)) + :initial-element (aref data1 0 0) + :element-type (array-element-type data1))) + ) + + (dotimes (i nr1) + (dotimes (j nc1) + (setf (aref result i j) (aref data1 i j)))) + + (dotimes (i nr2 (make-instance 'matrix :data result)) + (dotimes (j nc2) + (setf (aref result i (+ j nc1)) (aref data2 i j)))) + ))) + + +(defmethod join%2 ((d1 array) (d2 matrix)) + (let ((nr1 (matrix-row-number d1)) + (nr2 (matrix-row-number d2)) + ) + (assert (cl:= nr1 nr2)) + (let* ((nc1 (matrix-column-number d1)) + (nc2 (matrix-column-number d2)) + (data1 d1) + (data2 (matrix-data d2)) + (result (make-array (list nr1 (+ nc1 nc2)) + :initial-element (aref data1 0 0) + :element-type (array-element-type data1))) + ) + + (dotimes (i nr1) + (dotimes (j nc1) + (setf (aref result i j) (aref data1 i j)))) + + (dotimes (i nr2 (make-instance 'matrix :data result)) + (dotimes (j nc2) + (setf (aref result i (+ j nc1)) (aref data2 i j)))) + ))) + + +(defmethod join%2 ((d1 matrix) (d2 array)) + (let ((nr1 (matrix-row-number d1)) + (nr2 (matrix-row-number d2)) + ) + (assert (cl:= nr1 nr2)) + (let* ((nc1 (matrix-column-number d1)) + (nc2 (matrix-column-number d2)) + (data1 (matrix-data d1)) + (data2 d2) + (result (make-array (list nr1 (+ nc1 nc2)) + :initial-element (aref data1 0 0) + :element-type (array-element-type data1))) + ) + + (dotimes (i nr1) + (dotimes (j nc1) + (setf (aref result i j) (aref data1 i j)))) + + (dotimes (i nr2 (make-instance 'matrix :data result)) + (dotimes (j nc2) + (setf (aref result i (+ j nc1)) (aref data2 i j)))) + ))) + + +(defmethod join%2 ((d1 vector) (d2 matrix)) + (let ((result (join%2 d1 (matrix-data d2)))) + (make-instance 'matrix :data result))) + + +(defmethod join%2 ((d2 matrix) (d1 vector)) + (let ((result (join%2 (matrix-data d2) d1))) + (make-instance 'matrix :data result))) + + +(defmethod join%2 ((d1 column) (d2 matrix)) + (let ((result (join%2 d1 (matrix-data d2)))) + (make-instance 'matrix :data result))) + + +(defmethod join%2 ((d2 matrix) (d1 column)) + (let ((result (join%2 (matrix-data d2) d1))) + (make-instance 'matrix :data result))) + + +;;; stack%2 -- + +(defmethod stack%2 ((d1 column) (d2 column)) + (column (concatenate 'vector (column-vector d1) (column-vector d2)))) + + +(defmethod stack%2 ((d1 vector) (d2 array)) + (let ((nc (array-dimension d2 1))) + (assert (cl:= (length d1) nc)) + (let* ((nr (1+ (array-dimension d2 0))) + (result (make-array (list nr nc) + :initial-element (aref d2 0 0) + :element-type (array-element-type d2))) + ) + (dotimes (i nc) + (setf (aref result 0 i) (aref d1 i))) + (loop for i from 1 below nr + do (dotimes (j nc) + (setf (aref result i j) (aref d2 (1- i) j)))) + result))) + + +(defmethod stack%2 ((d1 sequence) (d2 array)) + (let ((nc (array-dimension d2 1))) + (assert (cl:= (length d1) nc)) + (let* ((nr (1+ (array-dimension d2 0))) + (result (make-array (list nr nc) + :initial-element (aref d2 0 0) + :element-type (array-element-type d2))) + ) + (dotimes (i nc) + (setf (aref result 0 i) (elt d1 i))) + (loop for i from 1 below nr + do (dotimes (j nc) + (setf (aref result i j) (aref d2 (1- i) j)))) + result))) + + +(defmethod stack%2 ((d1 array) (d2 vector)) + (let ((nc (array-dimension d1 1))) + (assert (cl:= (length d2) nc)) + (let* ((nr (1+ (array-dimension d1 0))) + (nr-1 (1- nr)) + (result (make-array (list nr nc) + :initial-element (aref d1 0 0) + :element-type (array-element-type d1))) + ) + + (dotimes (i nr-1) + (dotimes (j nc) + (setf (aref result i j) (aref d1 i j)))) + (dotimes (i nc) + (setf (aref result nr-1 i) (aref d2 i))) + result))) + + +(defmethod stack%2 ((d1 array) (d2 sequence)) + (let ((nc (array-dimension d1 1))) + (assert (cl:= (length d2) nc)) + (let* ((nr (1+ (array-dimension d1 0))) + (nr-1 (1- nr)) + (result (make-array (list nr nc) + :initial-element (aref d1 0 0) + :element-type (array-element-type d1))) + ) + + (dotimes (i nr-1) + (dotimes (j nc) + (setf (aref result i j) (aref d1 i j)))) + (dotimes (i nc) + (setf (aref result nr-1 i) (elt d2 i))) + result))) + + +(defmethod stack%2 ((d1 array) (d2 array)) + (let ((nc1 (array-dimension d1 1)) + (nc2 (array-dimension d2 1)) + ) + (assert (cl:= nc1 nc2)) + (let* ((nr1 (array-dimension d1 0)) + (nr2 (array-dimension d2 0)) + (result (make-array (list (+ nr1 nr2) nc1) + :initial-element (aref d1 0 0) + :element-type (array-element-type d1))) + ) + + (dotimes (i nr1) + (dotimes (j nc1) + (setf (aref result i j) (aref d1 i j)))) + + (dotimes (i nr2 result) + (dotimes (j nc2) + (setf (aref result (+ i nr1) j) (aref d2 i j)))) + ))) + + +(defmethod stack%2 ((d1 matrix) (d2 matrix)) + (let ((nc1 (matrix-column-number d1)) + (nc2 (matrix-column-number d2)) + ) + (assert (cl:= nc1 nc2)) + (let* ((nr1 (matrix-row-number d1)) + (nr2 (matrix-row-number d2)) + (data1 (matrix-data d1)) + (data2 (matrix-data d2)) + (result (make-array (list (+ nr1 nr2) nc1) + :initial-element (aref data1 0 0) + :element-type (array-element-type data1))) + ) + + (dotimes (i nr1) + (dotimes (j nc1) + (setf (aref result i j) (aref data1 i j)))) + + (dotimes (i nr2 (make-instance 'matrix :data result)) + (dotimes (j nc2) + (setf (aref result (+ i nr1) j) (aref data2 i j)))) + ))) + + +(defmethod stack%2 ((d1 array) (d2 matrix)) + (let ((nc1 (matrix-column-number d1)) + (nc2 (matrix-column-number d2)) + ) + (assert (cl:= nc1 nc2)) + (let* ((nr1 (matrix-row-number d1)) + (nr2 (matrix-row-number d2)) + (data1 d1) + (data2 (matrix-data d2)) + (result (make-array (list (+ nr1 nr2) nc1) + :initial-element (aref data1 0 0) + :element-type (array-element-type data1))) + ) + + (dotimes (i nr1) + (dotimes (j nc1) + (setf (aref result i j) (aref data1 i j)))) + + (dotimes (i nr2 (make-instance 'matrix :data result)) + (dotimes (j nc2) + (setf (aref result (+ i nr1) j) (aref data2 i j)))) + ))) + + + +(defmethod stack%2 ((d1 matrix) (d2 array)) + (let ((nc1 (matrix-column-number d1)) + (nc2 (matrix-column-number d2)) + ) + (assert (cl:= nc1 nc2)) + (let* ((nr1 (matrix-row-number d1)) + (nr2 (matrix-row-number d2)) + (data1 (matrix-data d1)) + (data2 d2) + (result (make-array (list (+ nr1 nr2) nc1) + :initial-element (aref data1 0 0) + :element-type (array-element-type data1))) + ) + + (dotimes (i nr1) + (dotimes (j nc1) + (setf (aref result i j) (aref data1 i j)))) + + (dotimes (i nr2 (make-instance 'matrix :data result)) + (dotimes (j nc2) + (setf (aref result (+ i nr1) j) (aref data2 i j)))) + ))) + + +(defmethod stack%2 ((d1 vector) (d2 matrix)) + (let ((result (stack%2 d1 (matrix-data d2)))) + (make-instance 'matrix :data result))) + + +(defmethod stack%2 ((d2 matrix) (d1 vector)) + (let ((result (stack%2 (matrix-data d2) d1))) + (make-instance 'matrix :data result))) + + +(defmethod stack%2 ((d1 column) (d2 matrix)) + (let ((result (stack%2 d1 (matrix-data d2)))) + (make-instance 'matrix :data result))) + + +(defmethod stack%2 ((d2 matrix) (d1 column)) + (let ((result (stack%2 (matrix-data d2) d1))) + (make-instance 'matrix :data result))) + + +;;; slice -- +#| +(defmethod slice ((v vector) (x fixnum)) + (aref v x)) + + +(defmethod slice ((v vector) (sd list)) + (assert (<= (list-length sd) 2)) + (destructuring-bind (start &optional end) + sd + (let ((start (cond ((or (eq start '*) (eq start 'cl:*)) 0) + ((typep start `(mod ,array-total-size-limit)) start) + (t (error "Unrecognized start ~S in slice designator." + start)))) + (end (cond ((or (null end) (eq end '*) (eq end 'cl:*)) nil) + ((typep end `(mod ,array-total-size-limit)) end) + (t (error "Unrecognized end ~S in slice designator." + end)))) + ) + (subseq v start end)))) + + +(defmethod slice ((v column) (sd list)) + (column (slice (column-vector v) sd))) + + +(defmethod slice ((v vector) (sd (eql '*))) + (copy-vector v)) + + +(defmethod slice ((v vector) (sd (eql 'cl:*))) + (copy-vector v)) + + +(defmethod slice ((v column) (sd (eql '*))) + (copy-vector v)) + + +(defmethod slice ((v column) (sd (eql 'cl:*))) + (copy-vector v)) + + +(defmethod slice (())) + |# + + + +;;; ref -- + +(defmethod ref ((v vector) (i fixnum) &rest indexes) + (when indexes + (error "REF of a vector requires only one index but ~D were supplied." + (1+ (list-length indexes)))) + + + (aref v i)) + + +(defmethod ref ((v vector) (i increasing-range) &rest indexes) + (when indexes + (error "REF of a vector requires only one index or range but a range and ~D indexes were supplied." + (list-length indexes))) + + (loop for vi from (range-start i) + below (or (range-end i) (length v)) + by (range-step i) + collect (aref v vi) into vs + finally (return (coerce vs 'vector)))) + + +(defmethod ref ((v vector) (i decreasing-range) &rest indexes) + (when indexes + (error "REF of a vector requires only one index or range but a range and ~D indexes were supplied." + (list-length indexes))) + + (loop for vi from (or (range-end i) (1- (length v))) + downto (range-start i) + by (abs (range-step i)) + collect (aref v vi) into vs + finally (return (coerce vs 'vector)))) + + +(defmethod ref ((v vector) (i vector) &rest indexes) + (when indexes + (error "REF of a vector requires only one index or subset but a subset and ~D indexes were supplied." + (list-length indexes))) + + (loop for x across i + collect (aref v x) into vs + finally (return (coerce vs 'vector)))) + + +(defmethod ref ((v vector) (i list) &rest indexes) + (when indexes + (error "REF of a vector requires only one index or subset but a subset and ~D indexes were supplied." + (list-length indexes))) + + (loop for x in i + collect (aref v x) into vs + finally (return (coerce vs 'vector)))) + + +(defmethod ref ((v column) (i fixnum) &rest indexes) + (apply #'ref (column-vector v) i indexes)) + + +(defmethod ref ((v column) (i range) &rest indexes) + (column (apply #'ref (column-vector v) i indexes))) + + +(defmethod ref ((v column) (i vector) &rest indexes) + (column (apply #'ref (column-vector v) i indexes))) + + +(defmethod ref ((v column) (i list) &rest indexes) + (column (apply #'ref (column-vector v) i indexes))) + + +(defmethod ref ((m matrix) (i t) &rest indexes) + (when (cl:/= 1 (list-length indexes)) + (error "........")) + + (ref%2 m i (first indexes))) + + + +(defmethod ref%2 ((m matrix) (i fixnum) (j fixnum)) + (aref (matrix-data m) i j)) + + +(defmethod ref%2 ((m matrix) (i fixnum) (j vector)) + (loop with md = (matrix-data m) + for x across j + collecting (aref md i x) into result + finally return (coerce result 'vector))) + + +(defmethod ref%2 ((m matrix) (i vector) (j fixnum)) + (loop with md = (matrix-data m) + for x across i + collecting (aref md x j) into result + finally return (coerce result 'vector))) + + +(defmethod ref%2 ((m matrix) (i vector) (j vector)) + (loop with md = (matrix-data m) + for x across i + collecting (loop for y across j + collecting (aref md x y)) + into result + finally + (return (make-instance 'matrix + :data (make-array (list (length i) + (length j)) + :initial-contents result))))) + + + + + +;;;;=========================================================================== +;;;; Generic Interface Functions. + +(defun join (&rest elements) + (cl:reduce #'join%2 elements)) + + +(defun stack (&rest elements) + (cl:reduce #'stack%2 elements)) + +;;;; end of file -- vector-matrix.lisp -- Added: trunk/common-math/numerics/linear-algebra/vector.lisp ============================================================================== --- (empty file) +++ trunk/common-math/numerics/linear-algebra/vector.lisp Wed Aug 16 17:07:41 2006 @@ -0,0 +1,461 @@ +;;;; -*- Mode: Lisp -*- + +;;;; vector.lisp -- +;;;; There is no wrapper for vectors, but we do have a class (struct) +;;;; for 'column' vectors. +;;;; Apart from that there are some specialized generic functions on +;;;; vectors. + +(in-package "CL.MATH.LINEAR-ALGEBRA") + +(defparameter *default-vector-element-type* 'double-float) + + +(deftype column-array (&optional (type 'cl:number) (length 'cl:*)) + `(array ,type (,length 1))) + + +(defstruct (column (:constructor column (vector))) + "The Column Structure. +A representation for `column' vectors." + (vector #() :type (or list vector))) + + +(defmethod print-object ((c column) (s stream)) + (format s "(~S ~S)" 'column (column-vector c))) + + +;;;--------------------------------------------------------------------------- +;;; Protocol. + +(defgeneric vector-p (v) + (:method ((v vector)) t) + (:method ((v column)) t) + (:method ((v t)) nil) + ) + + +(defgeneric vector-element-type (v) + (:method ((v vector)) (array-element-type v)) + (:method ((v column)) (array-element-type (column-vector v))) + ) + + +(defgeneric vector-length (m) + (:method ((x vector)) (length x)) + (:method ((x column)) (length (column-vector x))) + ) + + +;;;--------------------------------------------------------------------------- +;;; Creation. + +(defun make-vector (n + &rest array-keys + &key + (element-type *default-vector-element-type*) + (column-p nil) + &allow-other-keys + ) + (remf array-keys :column-p) + (let ((new-v (apply #'make-array n + (append array-keys + (list :initial-element (coerce 0 element-type) + :element-type + element-type))))) + (if column-p + (column new-v) + new-v))) + + +(defun make-zero-vector (n + &rest array-keys + &key + (element-type *default-matrix-element-type*) + &allow-other-keys) + (when (getf array-keys :initial-element) + (warn "MAKE-ZERO-MATRIX passed and :INITIAL-ELEMENT argument.") + (remf array-keys :initial-element)) + (apply #'make-vector + n + (append array-keys + (list :initial-element (coerce 0 element-type) + :element-type element-type)))) + + +;;;--------------------------------------------------------------------------- +;;; copy-vector +;;; Uses the code from K Pitman and B. Margolin appeared on CLL a long +;;; long time ago. +;;; Note that at least on LW this code has some of the usual problems +;;; with DOUBLE-FLOATS. +#| +(defun copy-array (array) + (adjust-array (make-array (array-dimensions array) + :displaced-to array + :element-type (array-element-type array)) + (array-dimensions array) + :displaced-to nil)) +|# + +(defgeneric copy-vector (v) + (:documentation "Returns a (shallow) copy of vector V. +The effect are the same as CL:COPY-SEQ for CL vectors while columns +are treated accordingly.")) + + +(defmethod copy-vector ((v vector)) + "COPY-VECTOR method for VECTORs." + (cl:copy-seq v)) + + +(defmethod copy-vector ((v column)) + "COPY-VECTOR method for COLUMN vectors." + (column (cl:copy-seq (column-vector v)))) + + +(defmethod copy-vector ((v array)) + (check-type v column-array) + (adjust-array (make-array (array-dimensions v) + :displaced-to v + :element-type (array-element-type v)) + (array-dimensions v) + :displaced-to nil)) + + + +#| +(defmethod copy-vector ((v vector) + &optional + (result + (adjust-array + (make-array (length v) + :displaced-to v + :element-type (array-element-type v)) + (array-dimensions v) + :displaced-to nil) + result-supplied-p)) + "COPY-VECTOR method for VECTORs." + (when result-supplied-p + (assert (vector-p result)) + (assert (cl:= (length result) (vector-length v))) + (map-into result #'identity v)) + result) + + +(defmethod copy-vector ((c column) + &optional + (result + (let ((v (column-vector c))) + (column + (adjust-array + (make-array (length v) + :displaced-to v + :element-type (array-element-type v)) + (array-dimensions v) + :displaced-to nil))) + result-supplied-p)) + "COPY-VECTOR method for COLUMN vectors." + (when result-supplied-p + (assert (column-p result)) + (assert (cl:= (vector-length result) (vector-length c))) + (map-into (column-vector result) #'identity (column-vector c))) + result) +|# + + + +;;;--------------------------------------------------------------------------- +;;; Equality. + +(defmethod =%2 ((x vector) (y vector)) + (and (cl:= (length x) (length y)) + (every #'=%2 x y))) + + +(defmethod =%2 ((x vector) (y column)) + (and (or (= 1 (length x) (vector-length y)) ; Only cases when this may happen. + (= 0 (length x) (vector-length y))) + (=%2 x (column-vector y)))) + + +(defmethod =%2 ((y column) (x vector)) + (=%2 x y)) + + +;;;--------------------------------------------------------------------------- +;;; zerop + +(defmethod zerop ((x vector)) + (cl:every #'zerop x)) + + +(defmethod zerop ((x column)) + (cl:every #'zerop (column-vector x))) + + +;;;--------------------------------------------------------------------------- +;;; Addition. + +(defmethod +%2 :before ((x vector) (y vector) + &optional + (r nil r-supplied-p)) + (unless (eq y r) + (let ((xl (length x))) + (assert (cl:= xl (length y))) + (when r-supplied-p + (assert (vector-p r)) + (assert (cl:= xl (vector-length r))) + )))) + + +(defmethod +%2 :before ((x vector) (y column) + &optional + (r nil r-supplied-p)) + (declare (ignore r r-supplied-p)) + (assert (or (cl:= 1 (length x) (vector-length y)) + (cl:= 0 (length x) (vector-length y)))) + ) + + +(defmethod +%2 :before ((y column) (x vector) + &optional + (r nil r-supplied-p)) + (declare (ignore r r-supplied-p)) + (assert (or (cl:= 1 (length x) (vector-length y)) + (cl:= 0 (length x) (vector-length y)))) + ) + + +(defmethod +%2 :before ((x number) (y vector) + &optional + (r nil r-supplied-p)) + (unless (eq y r) + (when r-supplied-p + (assert (vector-p r)) + (assert (cl:= (length y) (vector-length r))) + ))) + + +(defmethod +%2 ((x number) (y vector) &optional (r (copy-vector y))) + (dotimes (i (length y) r) + (setf (aref r i) (+%2 x (aref y i))))) + + +(defmethod +%2 ((y vector) (x number) &optional (r (copy-vector y))) + (+%2 x y r)) + + +(defmethod +%2 ((x vector) (y vector) &optional (r (copy-vector y))) + (dotimes (i (length x) r) + (setf (aref r i) (+%2 (aref x i) (aref y i))))) + + +(defmethod +%2 ((x vector) (y column) + &optional (r (copy-vector (column-vector y)))) + ;; This is valid only for vectors and columns of length 1 (or 0). + (+%2 x (column-vector y) r)) + + +(defmethod +%2 ((x column) (y vector) &optional (r (copy-vector y))) + ;; This is valid only for vectors and columns of length 1 (or 0). + (+%2 (column-vector x) y r)) + + +;;;--------------------------------------------------------------------------- +;;; Subtraction. + +(defmethod -%2 :before ((x vector) (y vector) + &optional + (r nil r-supplied-p)) + (let ((xl (length x))) + (assert (cl:= xl (length y))) + (when r-supplied-p + (assert (vector-p r)) + (assert (cl:= xl (vector-length r))) + ))) + + +(defmethod -%2 :before ((x vector) (y column) + &optional + (r nil r-supplied-p)) + (declare (ignore r r-supplied-p)) + (assert (or (cl:= 1 (length x) (vector-length y)) + (cl:= 0 (length x) (vector-length y)))) + ) + + +(defmethod -%2 :before ((y column) (x vector) + &optional + (r nil r-supplied-p)) + (declare (ignore r r-supplied-p)) + (assert (or (cl:= 1 (length x) (vector-length y)) + (cl:= 0 (length x) (vector-length y)))) + ) + + +(defmethod -%2 :before ((x number) (y vector) + &optional + (r nil r-supplied-p)) + (unless (eq y r) + (when r-supplied-p + (assert (vector-p r)) + (assert (cl:= (length y) (vector-length r))) + ))) + + +(defmethod -%2 ((x vector) (y vector) &optional (r (copy-vector y))) + (dotimes (i (length x) r) + (setf (aref r i) (-%2 (aref x i) (aref y i))))) + + + +(defmethod -%2 ((x vector) (y column) + &optional (r (copy-vector (column-vector y)))) + ;; This is valid only for vectors and columns of length 1 (or 0). + (-%2 x (column-vector y) r)) + + +(defmethod -%2 ((x column) (y vector) &optional (r (copy-vector y))) + ;; This is valid only for vectors and columns of length 1 (or 0). + (-%2 (column-vector x) r r)) + + +;;;--------------------------------------------------------------------------- +;;; Multiplication. +;;; Note that left and right multiplications are different! +;;; They follow the matrix multiplication rules. +;;; Here we treat only the vector x column multiplication as the +;;; column x vector has a matrix as a result. + +(defmethod *%2 :before ((x vector) (y column) &optional r) + (declare (ignore r)) + (assert (cl:= (cl:length x) (vector-length y)))) + + +(defmethod *%2 ((x vector) (y column) &optional r) + (declare (ignore r)) + (let ((result (coerce 0 (array-element-type x))) + (y (column-vector y)) + ) + (dotimes (i (cl:length x) result) + (setf result (+%2 result (*%2 (aref x i) (aref y i))))) + )) + +;;; The next one breaks the correct semantics passing of arguments, +;;; but it is provided as a convenience. + +(defmethod *%2 ((x vector) (y vector) &optional r) + (declare (ignore r)) + (assert (cl:= (cl:length x) (cl:length y))) + (let ((result (coerce 0 (array-element-type x)))) + (dotimes (i (cl:length x) result) + (setf result (+%2 result (*%2 (aref x i) (aref y i))))))) + + +;;; Scalar multiplication. + +(defmethod *%2 :before ((x number) (v vector) &optional (r nil r-supplied-p)) + (when r-supplied-p + (assert (cl:vectorp r)) + (assert (cl:= (cl:length v) (vector-length r))))) + + +(defmethod *%2 :before ((x number) (v column) &optional (r nil r-supplied-p)) + (when r-supplied-p + (assert (column-p r)) + (assert (cl:= (vector-length v) (vector-length r))))) + + +(defmethod *%2 ((x number) (v vector) &optional (r (copy-vector v))) + (map-into r (lambda (e) (*%2 x e)) v)) + + +(defmethod *%2 ((v vector) (x number) &optional (r (copy-vector v))) + (map-into r (lambda (e) (*%2 x e)) v)) + + +(defmethod *%2 ((x number) (v column) &optional (r (copy-vector v))) + (map-into (column-vector r) (lambda (e) (*%2 x e)) v) + r) + + +(defmethod *%2 ((v column) (x number) &optional (r (copy-vector v))) + (*%2 x v r)) + + +;;;--------------------------------------------------------------------------- +;;; Division. +;;; Only the simple form of division is implemented here. +;;; The general form requires the inverse, which requires extra machinery + +#| +(defmethod /%2 ((x number) (y matrix) &optional (r (copy-matrix y) r-supplied-p)) + (declare (ignore r-supplied-p)) + (*%2 (/ x) y r)) + + +(defmethod /%2 ((y matrix) (x number) &optional (r (copy-matrix y) r-supplied-p)) + (declare (ignore r-supplied-p)) + (*%2 (/ x) y r)) +|# + +;;;--------------------------------------------------------------------------- +;;; Exponentiation. +;;; All errors here ftb. + + +;;;--------------------------------------------------------------------------- +;;; Transpose + +(defmethod transpose ((v vector) &optional (r (column (copy-vector v)) r-supplied-p)) + (if r-supplied-p + (map-into (column-vector r) #'identity v)) + r) + + +(defmethod transpose ((v column) + &optional + (r (copy-vector (column-vector v)) r-supplied-p)) + (if r-supplied-p + (map-into r #'identity (column-vector v)) + r)) + + +(defmethod transpose :before ((v vector) &optional (r nil r-supplied-p)) + (when r-supplied-p + (assert (column-p r)) + (assert (cl:= (cl:length v) (vector-length r))))) + + +(defmethod transpose :before ((v column) &optional (r nil r-supplied-p)) + (when r-supplied-p + (assert (cl:vectorp r)) + (assert (cl:= (cl:length r) (vector-length v))))) + + +;;;--------------------------------------------------------------------------- +;;; Other functions. + +(defun iota (n &optional (start 0) (step 1 step-supplied-p) + &aux (initial-contents ())) + (cond ((< n start) + (when step-supplied-p (assert (minusp step))) + (loop for x from start above n by step + collect x into xs + finally (setf initial-contents xs))) + ((> n start) + (when step-supplied-p (assert (plusp step))) + (loop for x from start below n by step + collect x into xs + finally (setf initial-contents xs))) + (t ())) + (make-array n :initial-contents initial-contents)) + + + +;;;--------------------------------------------------------------------------- +;;; Utilities. + + +;;;; end of file -- matrix.lisp -- Added: trunk/common-math/numerics/numerics.system ============================================================================== --- (empty file) +++ trunk/common-math/numerics/numerics.system Wed Aug 16 17:07:41 2006 @@ -0,0 +1,18 @@ +;;;; -*- Mode: Lisp -*- + +;;;; numerics.system -- + +(mk:defsystem "numerics" + :components ((:system "linear-algebra") + ) + ) + + +(eval-when (:load-toplevel :execute) + (mk:add-registry-location (make-pathname :name nil + :type nil + :defaults *load-truename*)) + ) + +;;;; end of file -- numerics.system -- +

1 0