Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
A matrix package for common lisp building on work by Mark Hoemmen, Evan Monroig, Tamas Papp and Rif.
Common Lisp Other
Tree: 7a943ea3ff

Fetching latest commit…

Cannot retrieve the latest commit at this time

Failed to load latest commit information.
bugs
src
.gitignore
ChangeLog
README.org
lisp-matrix.asd
lispMatrixTests.R

README.org

-*- mode: org -*-

Introduction

Lisp-Matrix is a reasonably modern, flexible system for numeric algebra in Lisp.

Rif, Tomas Papp, and Mark H contributed various critical pieces (code, design, code, testing, code, benchmarks, code, and more code).

I’m just a dwarf standing on the shoulders of giants.

Design

Basic approach:

Storage can be done at the lisp or external level. Selection is user-specific, since you might have various people doing one or either, depending oupon the nature of the computations that are required. For exmpel large LAPACK-ish computatiosn are required to be done externally whereas stuff which is more lisp-cerntreic culd be done internally.

We avoid deep copying whenever possible. current scablability of storage is not nearly as great as the storgeage required to do interesting things with biomedical data or economic data. Whoops. This means tht we spend time on strucutres where copying is optional, done when needed to break the connectino from the initial data storage.

We, of course, means “them”.

Requirements:

Common Lisp Systems

Works with: SBCL Others: not clear (though CLISP is targetted)

Packages:

From Rif:

org.middleangle.cl-blapack

org.middleangle.foreign-numeric-vector

From Tomas:

ffa

array-operations

From Mark:

lisp-matrix integration (also found in ffa and cl-blapack).

In general:

cffi (depends on babel, alexandria )

cl-utilities

iterate

metabang-bind

asdf-system-connections

lift (depends on trivial-timeout)

These dependencies need to be better worked out.

Others (AJR?)

xarray (generic array-like accessors)

Documentation

Need to autogenerate approach for documenting what we can do with this. Until then, simple reference:

(make-matrix)  instantiates any supported matrix type:
	         = lisp/foriegn
		 = single/double/complex-single/complex-double/integer

(mref mat x y) get/set

(bind2 mat1 mat2 :by [:row|:column] )

(diagonal mat)

(m* mat1 mat2) => selection of the correct ZYYmm type (gemm for general mat mult)
(m+ mat1 mat2)
(m- mat1 mat3)
(axpy a mat1 mat2) => (scalar * matrix) + matrix

Usage

Demo (broken things)

;;; Precursor systems
(in-package :cl-user)
;; (asdf:oos 'asdf:compile-op 'ffa :force t)
;; (asdf:oos 'asdf:compile-op 'array-operations :force t)

;; (asdf:oos 'asdf:compile-op 'org.middleangle.foreign-numeric-vector :force t)
;; (asdf:oos 'asdf:compile-op 'org.middleangle.cl-blapack :force t) ;  :force t

;;; The main thing...
;; (delete-package 'lisp-matrix) ;; fails, but we need to cleanup a bit more.

;; (asdf:oos 'asdf:compile-op 'lisp-matrix :force t)
;; (asdf:oos 'asdf:compile-op 'lisp-matrix)
;; (asdf:oos 'asdf:load-op 'lisp-matrix)

;; (asdf:oos 'asdf:compile-op 'cffi :force t)

(in-package :lisp-matrix-unittests)
;; Tests = 69, Failures = 0, Errors = 12 ;; 26.2.2009
(run-tests :suite 'lisp-matrix-ut)
(describe (run-tests :suite 'lisp-matrix-ut))
;; or simply...
(run-lisp-matrix-tests)
(describe  (run-lisp-matrix-tests))

;; failures:

;; Note that when unit tests fail in m*- tests, it seems to do with a
;; "macro vs defun" problem, related to compile-time vs. run-time
;; evaluation that I (tony) am not quite understanding, causing a
;; possible increase in the number of errors beyond the number
;; reported above.
;;
;; The current two errors are:  
;; * foreign arrays with integer values are not supported.
;; * mixed CL-BLAPACK calls are not yet supported (lisp/foreign stored
;;   matrix-like calls).
;; I'm sure there will be more.

(in-package :lisp-matrix-user)

;; (lisp-matrix-unittests:run-lisp-matrix-tests)
;; (describe (lisp-matrix-unittests:run-lisp-matrix-tests))

(describe 
 (lift::run-test
  :test-case  'lisp-matrix-unittests::strided-matrix-column-access
  :suite 'lisp-matrix-ut-vectors))


;; Here is what we need to fix, based on the above:
;; #  creation of foreign-array matrices which are integer valued
;;    fails.


;; Just a reminder:
;; (typep -1 '(integer 0 *))   ;=> nil
;; (typep  2 '(integer 0 *))   ;=> T
;; (typep  3 '(integer -1 2))  ;=> nil
;; (typep  2 '(integer -1 2))  ;=> T

;;; FIXME FOLLOWING ERRORS: MIGRATE INTO UNITTESTS...

(progn  ;;#FIXME: writing out R matrices -- as strings and via RCLG

  (defparameter *x-temp*
    (make-matrix 4 5
		 :implementation :lisp-array
		 :element-type 'double-float
		 :initial-contents #2A((11d0 12d0 13d0 14d0 15d0)
				       (21d0 22d0 23d0 24d0 25d0)
				       (31d0 32d0 33d0 34d0 35d0)
				       (41d0 42d0 43d0 44d0 45d0))))

  ;; bad:  (min (values (list 4d0 2d0 3d0 5d0 3d0)))
  (reduce #'min (list 4d0 2d0 3d0 5d0 3d0))
  (reduce #'min (list 2d0 4d0 3d0 5d0 3d0))
  (reduce #'min (list 4d0 3d0 5d0 3d0 2d0))

  (reduce #'(lambda (x y) (concatenate 'string x y))
	  "test"
	  " "
	  (list "a2" " s3 " "asdf")
	  "end.")

  (defun lispmatrix2r (m &key (rvarname "my.mat"))
    "Write out a string that can be used to read in the matrix into R.
Used for creating verfication scripts and test cases."
    (check-type m matrix-like)
    (apply 
     #'concatenate 'string
     (format nil "~%~s <- matrix ( data = c(" rvarname)
     (let ((result (list)))
		    (dotimes (i (matrix-dimension m 0))
		      (dotimes (j (matrix-dimension m 1))
			(cons (format nil "~d," (mref m i j)) result)))
		    (reverse result))
     (list  (format nil "), nrows=~d, ncols=~d, by.row=TRUE)"
	     (matrix-dimension m 0)
	     (matrix-dimension m 1)))))

  (lispmatrix2R *x-temp*)


  (let ((result (make-array (list 3 5) :element-type 'string)))
    (dotimes (i 3)
      (dotimes (j 5)
	(format t "~s ~s ~%" i j)
	(setf (aref result i j) (format t "(~d ~d)," i j))))
    (reverse result))

  )


#+nil 
(progn   ;; QR decomp

  (let* ((state1 (make-random-state))
	 (state2 (make-random-state state1)))
    (m= (rand 2 3 :state state1)
	(rand 2 3 :state state2)))

  ;;; Problems here...
  (geqrf (make-matrix 2 2 :initial-contents #2A(( 1d0 2d0 ) (2d0 1d0))))
  (geqrf (make-matrix 2 2 :initial-contents '(( 1d0 2d0 ) (2d0 1d0))))
  ;;  (make-vector 2 :type :column :initial-contents '((1d0)(1d0))))

  )


#+nil
(progn ;; FIXME: R's apply across array indicies

  ;; Thought 1 (currently not planned for implementation)
  ;; consider using affi as a general iterator/walker generator.
  ;; So, R has a notion of apply, sapply, tapply, lapply -- what we
  ;; should do is something like
  ;;
  ;;     (map-matrix with-fn this-matrix
  ;;                 :by iterator
  ;;                 :result-type 'list)
  ;;
  ;; silly or for later:        :computation-type [:parallel|:serial]
  ;;
  ;; or similar, where :result-type is something that can be coerced to
  ;; from a sequence, and computation-type might drive whether there are
  ;; dependencies or not.   (this last is probably too premature).

  ;; The basic idea is to use vector functions (taking a vector, and
  ;; returning a object) and use them to provide an object that can be
  ;; part of a list (or generally, a sequence of homogeneous objects).

  ;; Reviewing Tamas Papp's affi package provides one approach to this
  ;; challenge.  He suggests that an obvious approach would be to
  ;; break up the 2 actions needed for selection consist of describing
  ;; the mapping from array to structure, and then walking the
  ;; structure to extract (for copy or use).  For our needs, we need a
  ;; means of doing this to partition the space, and then
  ;; post-partition, deciding which partitions need to be considered
  ;; for further processing, and which ones get discarded.

  ;; So to clarify how this might work: 
  ;; 1. we need a function which takes a matrix and creates a list of
  ;; matrix-like or vector-like elements.
  ;; 2. we have functions which operate in general on matrix-like or
  ;; vector-like objects.
  ;; 3. we use mapcar or similar to create the results.  
  ;; 3a. multi-value return could be used to create multiple lists of
  ;; vector-like or matrix-like objects, for example to get a complex
  ;; computation using inner-products.   So for instance:
  ;;   list1: v1a v2a v3a
  ;;   list2: m1  m2  m3
  ;;   list3: v1b v2b v3b
  ;; and we compute
  ;;   (list-of (IP v#a m1 v#b )) 
  ;; via
  ;;   (mapcar #'IP (list-of-vector-matrix-vector M))

  ;; We would need such an "extractor" to make things work out right.  
  #+nil(mapcar #'function-on-matrix (make-list-of-matrices original-matrix)) 


  (list->vector-like (list 1d0 2d0 3d0) :orientation :row)

  (make-vector 3 :type :column 
	       :initial-contents
	       (mapcar #'(lambda (x) (list (coerce x 'double-float)))
		       (list 1d0 2d0 3d0)))

  (make-vector 3 :type :row 
	       :initial-contents
	       (list (mapcar  #'(lambda (x) (coerce x 'double-float))
			      (list 1d0 2d0 3d0))))

  ;; The following approach would be required to do a proper map-back.
  #+nil(list->vector-like (map 'list #'function-of-2-args (list1) (list2)) :type :row) ; or :column
  ;; this would take a list and create an appropriate vector-like of
  ;; the appropriate type.

  ;; Thought 2, the current immediate approach:
  ;; What we currently do is break it out into components.

  (defparameter *m1-app* (ones 2 3))
  (let ((col-list (list-of-columns *m1-app*)))
    (dotimes (i (length col-list))
	  (princ (v= (nth i col-list)
		      (ones 2 1)))))

  (list-of-columns *m1-app*)
  (list-of-rows *m1-app*)
  
  (mapcar #'princ (list-of-columns *m1-app*))

  (format nil "R-Apply approach"))


#+nil
(progn
  ;; Studies in Class inheritance

  (subtypep 'LA-SIMPLE-VECTOR-DOUBLE 'VECTOR-LIKE)
  (subtypep 'LA-SLICE-VECVIEW-DOUBLE 'VECTOR-LIKE)
  (subtypep 'LA-SIMPLE-VECTOR-DOUBLE 'LA-SLICE-VECVIEW-DOUBLE)
  (subtypep  'LA-SLICE-VECVIEW-DOUBLE 'LA-SIMPLE-VECTOR-DOUBLE)

  (subtypep 'FA-SIMPLE-VECTOR-DOUBLE 'MATRIX-LIKE)

  ;;; weird!
  (m- (make-vector 2 :initial-contents '((1d0 1d0)))
      (make-vector 2 :initial-contents '((1d0 1d0))))

  (let ((*default-implementation* :foreign-array))
    (m- (make-vector 2 :initial-contents '((1d0 1d0)))
	(make-vector 2 :initial-contents '((1d0 1d0)))))

  (let ((*default-implementation* :lisp-array))
    (m- (make-vector 2 :initial-contents '((1d0 1d0)))
	(make-vector 2 :initial-contents '((1d0 1d0)))))

  (m- (make-vector 2
		   :implementation :lisp-array
		   :initial-contents '((1d0 1d0)))
      (make-vector 2
		   :implementation :foreign-array
		   :initial-contents '((1d0 1d0))))

  (typep  (first *lm-result*) 'vector-like)
  (typep  (first *lm-result*) 'matrix-like)
  (typep  (second *lm-result*) 'vector-like)
  (typep  (second *lm-result*) 'matrix-like)
  (typep *x-temp* 'vector-like)
  (typep *x-temp* 'matrix-like) ; => T ,rest of this paragraph are false.

  (m- *x-temp* *x-temp*))

Demo (working things)

Demos for Lisp Matrix (encoded within progn’s)

  1. instantiating matrices and vectors
  2. inversion using BLAS/LAPACK
(in-package :lisp-matrix-user)

(progn ;; data object instantiation

  (defparameter *m01*
    (make-matrix
     6 5
     :initial-contents '((11d0 12d0 13d0 14d0 15d0)
			 (21d0 22d0 23d0 24d0 25d0)
			 (31d0 32d0 33d0 34d0 35d0)
			 (41d0 42d0 43d0 44d0 45d0)
			 (51d0 52d0 53d0 54d0 55d0)
			 (61d0 62d0 63d0 64d0 65d0)))
    "6x5 matrix with entries representing row+1,col+1 values, for
     test purposes.")

  (documentation  '*m01* 'variable)

  (defparameter *m1-ex*  (make-matrix 2 5
			   :implementation :lisp-array  ;; :foreign-array
			   :element-type 'double-float)
    "quick variable initialized to zeros")
    
  (defparameter *m2-la-int*
    (make-matrix 2 5
		 :implementation :lisp-array  ;; :foreign-array
		 :element-type 'integer ; 'double-float
		 ;; :initial-contents (list 1 2 3 4 5 6 7 8 9 10)
		 :initial-contents #2A((1 2 3 4 5)
				       (6 7 8 9 10)))
    "placeholder 2")

  ;; Currently we can make a foriegn matrix of doubles, but not a
  ;; foriegn matrix of integers.
  (defparameter *m2-fa*
    (make-matrix
     2 5
     :implementation :foreign-array 
     :element-type 'double-float
     :initial-contents #2A(( 1d0 2d0 3d0 4d0 5d0)
			   ( 6d0 7d0 8d0 9d0 10d0)))
    "placeholder 2")

  (defparameter *m2-la*
    (make-matrix
     2 5
     :implementation :lisp-array 
     :element-type 'double-float
     :initial-contents #2A(( 1d0 2d0 3d0 4d0 5d0)
			   ( 6d0 7d0 8d0 9d0 10d0)))
    "placeholder 2")


  (defparameter *m3-fa*
    (make-matrix
     2 2
     :implementation :foreign-array 
     :element-type 'double-float
     :initial-contents #2A(( 1d0 2d0 )
			   ( 6d0 7d0 )))
    "placeholder 2")

  (defparameter *m3-la*
    (make-matrix
     2 2
     :implementation :lisp-array 
     :element-type 'double-float
     :initial-contents #2A(( 1d0 2d0 )
			   ( 6d0 7d0 )))
    "placeholder 2")

    
  (defparameter *m01b*
    (strides *m01* :nrows 2 :ncols 3
	     :row-stride 2
	     :row-offset 1 :col-offset 1))
  
  (defparameter *m01c* 
    (window *m01*
	    :nrows 2 :ncols 3
	    :row-offset 2 :col-offset 1))
					; EVAL BELOW TO SETUP DATA


  ;; data for lls estimation
  (defparameter *xv*
    (make-vector
     8
     :type :row ;; default, not usually needed!
     :initial-contents '((1d0 3d0 2d0 4d0 3d0 5d0 4d0 6d0))))

  ;; col vector
  (defparameter *xv2*
    (make-vector
     8
     :type :column
     :initial-contents '((1d0)
			 (3d0)
			 (2d0)
			 (4d0)
			 (3d0)
			 (5d0)
			 (4d0)
			 (6d0))))

  (v= *xv* *xv2*) ; => T
  (m= *xv* *xv2*) ; => nil

  (defparameter *xv+1*
    (make-matrix
     8 2
     :initial-contents '((1d0 1d0)
			 (1d0 3d0)
			 (1d0 2d0)
			 (1d0 4d0)
			 (1d0 3d0)
			 (1d0 5d0)
			 (1d0 4d0)
			 (1d0 6d0))))

  (defparameter *xv+1a*
    (make-matrix
     8 2
     :initial-contents #2A((1d0 1d0)
			   (1d0 3d0)
			   (1d0 2d0)
			   (1d0 4d0)
			   (1d0 3d0)
			   (1d0 5d0)
			   (1d0 4d0)
			   (1d0 6d0))))

  (defparameter *xv+1b*
    (bind2
     (ones 8 1)
     (make-matrix
      8 1
      :initial-contents '((1d0)
			  (3d0)
			  (2d0)
			  (4d0)
			  (3d0)
			  (5d0)
			  (4d0)
			  (6d0)))
     :by :column))

  (m= *xv+1a* *xv+1b*) ; => T

  (defparameter *xm*
    (make-matrix
     2 8
     :initial-contents '((1d0 3d0 2d0 4d0 3d0 5d0 4d0 6d0)
			 (1d0 2d0 3d0 4d0 5d0 6d0 7d0 8d0))))

  (defparameter *y*
    (make-vector
     8
     :type :row
     :initial-contents '((1d0 2d0 3d0 4d0 5d0 6d0 7d0 8d0))))

  (defparameter *y2*
    (make-vector
     8
     :type :column
     :initial-contents '((1d0)
			 (2d0)
			 (3d0)
			 (4d0)
			 (5d0)
			 (6d0)
			 (7d0)
			 (8d0))))
  (transpose *y2*)




  (format nil "Data set up"))

#+nil
(progn 
  ;; Tests for square matrices...
  (trap2mat (rand 3 3))

  (trap2mat (make-matrix 3 3
			 :initial-contents #2A((1d0 2d0 3d0)
					       (4d0 5d0 6d0)
					       (7d0 8d0 9d0))))
  (trap2mat (make-matrix 3 3
			 :initial-contents #2A((1d0 2d0 3d0)
					       (4d0 5d0 6d0)
					       (7d0 8d0 9d0)))
	    :type :lower)
  (trap2mat (make-matrix 3 3
			 :initial-contents #2A((1d0 2d0 3d0)
					       (4d0 5d0 6d0)
					       (7d0 8d0 9d0)))
	    :type :upper)

  ;; need to write unit tests for square and rect matrices.
  )


#+nil
(progn
  ;; factorization and inversion via LAPACK

  ;; LU
  (let ((test-eye (eye 7 7)))
    (m* test-eye (minv-lu test-eye)))

  ;; Cholesky
  (let ((myrand (rand 4 4)))
    (princ myrand)
    (princ (matrix-like-symmetric-p (m* (transpose myrand) myrand)))
    (princ (m*  (m* (transpose myrand) myrand)
		(minv-cholesky  (m* (transpose myrand) myrand))))))


(progn  
  ;; Using xGEQRF routines for supporting linear regression.

  ;; Question: Need to incorporate the xGEQRF routines, to support
  ;; linear regression work?

  ;; LAPACK suggests to use the xGELSY driver (GE general matrix, LS
  ;; least squares, need to lookup Y intent (used to be an X alg, see
  ;; release notes).

  (let ((a (rand 10 5)))
    (geqrf a)))

TODO.lisp : things that don’t work but should lm-demo.lisp : things that might work but should

Demo (getting started)

(in-package :cl-user)
(asdf:oos 'asdf:load-op :lisp-matrix)

Demo (more working things)

;;; This file illustrates some common actions in the course of working
;;; with matrices using lisp-matrix.  It is important to note that
;;; there are better ways to do this, that this are to help introduce
;;; usage, not describe best practices for using this system.

;;; = Precursor systems
;;  (asdf:oos 'asdf:compile-op 'ffa :force t)
;;  (asdf:oos 'asdf:compile-op 'org.middleangle.foreign-numeric-vector :force t)
;;  (asdf:oos 'asdf:compile-op 'org.middleangle.cl-blapack :force t)

;;; = The maing thing...
;; (asdf:oos 'asdf:compile-op 'lisp-matrix :force t)
;; (asdf:oos 'asdf:compile-op 'lisp-matrix)

;;; And the only thing that ought to be required;
(asdf:oos 'asdf:load-op 'lisp-matrix)

;;; Check status of the installation...

(in-package :lisp-matrix-unittests)
(run-lisp-matrix-tests)

;; if the above describes errors, here is how we figure out what bug
;; report to write...

(describe  (run-lisp-matrix-tests))

;;; Now we can use it, either by importing the symbols into the
;;; current package by:

;; (use-package :lisp-matrix)

;;; or by trying it out in the -user package, before implementing for
;;; production usage.

(in-package :lisp-matrix-user)

;; (lisp-matrix-unittests:run-lisp-matrix-tests)
;; (describe (lisp-matrix-unittests:run-lisp-matrix-tests))

;;; We wrap these up into a progn for simple overall evaluation, but
;;; stepping through them is fine as well.

(progn 
  
  ;; make some matrices
  (defparameter *m1* (make-matrix 2 5
			:implementation :lisp-array  ;; :foreign-array
			:element-type 'double-float)
    "placeholder 1")
  
  ;; works, as it should.  Indexing is zero-based, so we get the first
  ;; element by...
  (mref *m1* 0 0)
  (mref *m1* 1 3)
  (setf (mref *m1* 1 3) 1.2d0)
  *m1*


  ;; increase complexity

  (defparameter *m2* (make-matrix 2 5
			:implementation :lisp-array  ;; :foreign-array
			:element-type 'integer ; 'double-float
			;; :initial-contents (list 1 2 3 4 5 6 7 8 9 10)
			:initial-contents #2A(( 1 2 3 4 5)
					      ( 6 7 8 9 10)))
    "placeholder 2")

  (defparameter *m2a*
    (make-matrix 2 5
		 :implementation :lisp-array  ;; :foreign-array
		 :element-type 'integer ; 'double-float
		 :initial-contents '((1 2 3 4 5)
				     (6 7 8 9 10)))
    "placeholder...")

  ;; Currently we can make a foriegn matrix of doubles, but not a
  ;; foreign matrix of integers.  If we are working with smaller
  ;; matrices and are not doing a great deal of matrix algebra, then
  ;; we probably prefer :lisp-array rather than :foreign-array.
  (defvar *m2b*
    (make-matrix 2 5
		 :implementation :foreign-array 
		 :element-type 'double-float
		 :initial-contents #2A(( 1d0 2d0 3d0 4d0 5d0)
				       ( 6d0 7d0 8d0 9d0 10d0)))
    "placeholder 2")
  *m2b*

  (mref *m2b* 0 2) ;; => 3
  *m2b*
  (transpose *m2b*)

  ;; simple subsetting is simple
  (m= (row *m2b* 0)
      (col (transpose *m2b*) 0)) ; => nil, orientation
  (v= (row *m2b* 0)
      (col (transpose *m2b*) 0)) ; => T, no orientation worries

  (m= (col *m2b* 0)
      (row (transpose *m2b*) 0))
  (v= (col *m2b* 0)
      (row (transpose *m2b*) 0))


  (defvar *m3*
    (make-matrix 6 5 :initial-contents '((1d0 2d0 3d0 4d0 5d0)
					 (6d0  7d0  8d0  9d0  10d0)
					 (11d0 12d0 13d0 14d0 15d0)
					 (16d0 17d0 18d0 19d0 20d0)
					 (21d0 22d0 23d0 24d0 25d0)
					 (26d0 27d0 28d0 29d0 30d0)))
    "placeholder 3")

  (row *m3* 2)
  (col *m3* 1)


  (= (mref *m3* 0 1)
     (mref (transpose *m3*) 1 0))

  (=  (mref *m3* 2 2)
      (mref (transpose *m3*) 2 2))

  *m3*
  (transpose *m3*)

  ;;; Now we play with striding and slicing subsets.  These work well
  ;;; for simple subsetting which can be done by counting/enumeration
  ;;; on some form of regular scale.

  ;;; In addition, equality is somewhat important for numerical
  ;;; issues.  Right.  Anyway, for matrices it is mostly clear what to
  ;;; do, but for vectors, which are inheriting from matrices, we have
  ;;; 2 issues.  The first is the obvious, the numerical values, and
  ;;; the second is not quite obvious, which is the metadata
  ;;; surrounding the difference between an MxN and NxM matrix.  For
  ;;; the first, think about v= and for the second, m= is the right
  ;;; function.

  (defvar *m4* (strides *m3* :nrows 2 :row-stride 2)
    "yet another placeholder.")
  *m4*
  (m= (row *m4* 0)
      (make-matrix 1 5 :initial-contents '((1d0 2d0 3d0 4d0 5d0))))
  (m= (row *m4* 1)
      (make-matrix 1 5 :initial-contents '((11d0 12d0 13d0 14d0 15d0))))
  ;; note the redoing for the columns -- different!
  (m= (col *m4* 0)
      (make-matrix 2 1 :initial-contents '((1d0) (11d0))))
  (m= (col *m4* 1)
      (make-matrix 2 1 :initial-contents '((2d0) (12d0))))

  (v= (row *m4* 0) (col (transpose *m4*) 0))
  (v= (col *m4* 0) (row (transpose *m4*) 0))

  *m4*
  (row *m4* 0)
  (col *m4* 4)


  (let* ((*default-element-type* '(complex double-float))
	 (m1 (axpy #C(1.0d0 0.0d0)
		   (ones 2 2)
		   (scal #C(1.5d0 0.0d0)
			 (ones 2 2))))
	 (m2 (scal #C(2.5d0 0.0d0) (ones 2 2)))
	 (m3 (axpy #C(-1.0d0 0.0d0)
		   (ones 2 2)
		   (scal #C(1.5d0 0.0d0) (ones 2 2))))
	 (m4 (scal #C(0.5d0 0.0d0) (ones 2 2))))
    (format t "~A ~A ~%"
	    (m= m1 m2)
	    (m= m3 m4)))

  (m+ (row m3 1) (row m3 2))
  (m- (row m3 1) (row m3 2))

  )



;;; EXAMPLES TO DEMONSTRATE


;;; consider the following matrix:
;;; n1= 11 12 13
;;;     21 22 23
(defparameter *n1*
  (make-matrix 2 3
	       :implementation :lisp-array
	       :element-type 'double-float
	       :initial-contents #2A ((11d0 12d0 13d0)
				      (21d0 22d0 23d0))))
*n1*
;;; then storage in row-major orientation would be a sequence
;;;     11 12 13 21 22 23
;;; while in column-major orientation it would be
;;;     11 21 12 22 13 23 
;;; At this point, consider the following.  Suppose we have a matview
;;; with dims 1x3, row/col offset 1,0:
;;; n2= 21 22 23
(defparameter *n2*
  (window *n1*
	  :nrows 1 :ncols 3
	  :row-offset 1 :col-offset 0))
*n2*
;;; or alternatively dims 2x2, row/col offset 0,1:
;;; n3= 12 13
;;;     22 23
(defparameter *n3*
  (window *n1*
	  :nrows 2 :ncols 2
	  :row-offset 0 :col-offset 1))
*n3*
;;;
;;; for the first, we see that, by orientation, we have the following:
;;;     .. .. .. 21 22 23   (row-major)
;;;     .. 21 .. 22 .. 23   (column-major)
;;; 
;;; so we see that for
;;; row-major:    index=3 (ncols), stride=1
;;; column-major: index=1 (ncols), stride=2 (nrows)
;;; 
;;; for the second, by orientation, we have:
;;;     .. 12 13 .. 22 23  (row-major)
;;;     .. 12 22 .. 13 23  (column-major)
;;; 
;;; so we see that for
;;; row-major:    index=1 (ncols), stride=2 (ncols)
;;; column-major: index=1,(nrows), stride=3 (nrows)
;;; 
;;; Consider a more complex matrix:
;;; 
;;; o1= 11 12 13 14 15
;;;     21 22 23 24 25
;;;     31 32 33 34 35
;;;     41 42 43 44 45
(defparameter *o1*
  (make-matrix 4 5
	       :implementation :lisp-array
	       :element-type 'double-float
	       :initial-contents #2A ((11d0 12d0 13d0 14d0 15d0)
				      (21d0 22d0 23d0 24d0 25d0)
				      (31d0 32d0 33d0 34d0 35d0)
				      (41d0 42d0 43d0 44d0 45d0))))
*o1*
;;; row-major:
;;;    o1= 11 12 13 14 15 21 22 23 24 25 31 32 33 34 35 41 42 43 44 45
;;; col-major: 
;;;    o1= 11 21 31 41 12 22 32 42 13 23 33 43 14 24 34 44 15 25 35 45
;;;
;;;
;;; Then a matview, dims 3, offset 2,1 :
;;;
;;; o2= 32 33 34
;;;     42 43 44
(defparameter *o2*
  (window *o1*
	  :nrows 2 :ncols 3
	  :row-offset 2 :col-offset 1))
*o2*
;;;
;;; and a strided matview, indexed, could be (offset 2,3; row-stride 2)
;;;
;;; o3= 23 24 25
;;;     43 44 45
(defparameter *o3*
  (strides *o1*
	   :nrows 2 :ncols 3
	   :row-offset 1 :col-offset 2
	   :row-stride 2 :col-stride 1))
*o3*
;;; and for where this sits in the original matrix...
;;;
;;; and now to pull out the rows and columns via slicing on a strided
;;; matrix, we have the following approaches, for the zero-th column:
;;;     23
;;;     43
(slice *o3* :offset 0 :stride 1 :nelts (nrows *o3*) :type :column)
(parent *o3*)
;;; and for the 2nd column (3rd, since we are zero counting).
;;;     25
;;;     45
(slice *o3* :offset 4 :stride 1 :nelts (nrows *o3*) :type :column)
;;; and for the 1st row (2nd, again zero-counting):
;;;     43 44 45
(slice *o3* :offset 1 :stride 2 :nelts (ncols *o3*) :type :row)
;;; 
(orientation *o3*)

;; convert between foriegn-array and lisp-array.

;; operate ()

;; do some blas/lapack

;; output

;; Windowing -- simple, works!
(m= (col *c* 0)
    (make-matrix 3 1 :initial-contents '((16d0) (21d0) (26d0))))
(m= (col *c* 1) 
    (make-matrix 3 1 :initial-contents '((17d0) (22d0) (27d0))))
(m= (col *c* 2)
    (make-matrix 3 1 :initial-contents '((18d0) (23d0) (28d0))))
(m= (col *c* 3)
    (make-matrix 3 1 :initial-contents '((19d0) (24d0) (29d0))))
(m= (col *c* 4)
    (make-matrix 3 1 :initial-contents '((20d0) (25d0) (30d0))))

(m= (col *d* 0)
    (make-matrix 3 1 :initial-contents '((18d0) (23d0) (28d0))))
(m= (col *d* 1) 
    (make-matrix 3 1 :initial-contents '((19d0) (24d0) (29d0))))

;; do we want this as part of the API? Currently fails.
;; (m= (col *c* 4)
;;     (col *c* 4)
;;     (make-matrix 3 1 :initial-contents '((20d0) (25d0) (30d0))))


;;;;;;;;


;; strided matrix col access
m01b
(orientation m01b)
(unit-strides-p m01b) ;; false, it's explicitly strided
(parent m01b)
(orientation  (parent m01b))
(unit-strides-p (parent m01b)) ;; true, it's the original...

;; Windowed matrix
(orientation m01c)
(row m01c 0) ; Y
(row m01c 1) ; Y
(col m01c 0) ; Y
(col m01c 1) ; Y
(col m01c 2) ; Y

;; slice matrix access to rows
(row m01b 0) ; Y
(row m01b 1) ; Y
(orientation m01b) (offset m01b)
(row-offset m01b) (col-offset m01b)
(col m01b 0) ; N
(col m01b 1) ; N...
(col m01b 2)
(col m01b 3)

(slice m01b :offset 0 :stride 2 :nelts (ncols m01b) :type :row)
(slice (parent m01b) ; equiv on parent
       :offset 1
       :stride 2
       :nelts (ncols m01b)
       :type :row)
;; 
(slice m01b :offset 1 :stride 2 :nelts (ncols m01b) :type :row)
(slice (parent m01b) ; equiv on parent
       :offset 1
       :stride 2
       :nelts (ncols m01b)
       :type :row)

;; slice matrix access to columns
(slice m01b :offset 0 :stride 1 :nelts (nrows m01b) :type :column)
(col m01b 0)
(slice m01b :offset 2 :stride 1 :nelts (nrows m01b) :type :column)
(col m01b 1)
(slice m01b :offset 4 :stride 1 :nelts (nrows m01b) :type :column)
(col m01b 2)
(slice m01b :offset 6 :stride 1 :nelts (nrows m01b) :type :column)
(col m01b 3)
(offset m01b)
(row-stride m01b) ; => 2
(col-stride m01b) ; => 1

  (m= (col m01b 0)
      (make-matrix 2 1 :initial-contents '((11d0) (31d0))))
  (m= (col m01b 1)
      (make-matrix 2 1 :initial-contents '((12d0) (32d0))))
  (m= (col m01b 2)
      (make-matrix 2 1 :initial-contents '((13d0) (33d0))))
  (m= (col m01b 3)
      (make-matrix 2 1 :initial-contents '((14d0) (34d0))))
  (m= (col m01b 4)
      (make-matrix 2 1 :initial-contents '((15d0) (35d0))))
  (row m01b 0)
  (row m01b 1)
  (col m01b 0)
  (col m01b 1)

  
  ;; FIXME: there are bugs in slicing/striding with transposed
  ;; matrices. 

  ;; the following are correct, but..
  (row m01 0)
  (row m01 1)
  (row m01 2)
  (row m01 3)

  (col m01 0)
  (col m01 1)
  (col m01 2)
  (col m01 3)

  m01
  (transpose m01)
  (row (transpose m01) 0)
  (row (transpose m01) 1) ; wrong: grab bad column, AND by 1 (pushed up)
  (row (transpose m01) 2) ; ditto, wrong by 2
  (row (transpose m01) 3) ; etc...wrong by 3

  (row (transpose m01) 0)
  (transpose (row (transpose m01) 0))

  m01
  (transpose m01)
  (col (transpose m01) 0)
  (col (transpose m01) 1) ; last rather than first
  (col (transpose m01) 2) ;
  (col (transpose m01) 3) ; ditto above


  (v= (row m01 0)
      (col (transpose m01) 0)) ;; works
  
  (m= (row m01 0)
      (col (transpose m01) 0)) ;; fails, since dims unequal
  
  m01
  (transpose m01)
  ;; given the above...
  ;; FIXME: Big Barf!
  (v= (row m01 1)
      (col (transpose m01) 1) ) ;; fails badly.  Real badly.
  
  (v= (col m01 1)
      (row (transpose m01) 1) ) ;; fails, but closer...
  
  (col m01 1)
  (col (transpose m01) 1) ;; this is the problem, indexing issue...
  
  
  ;; and the same problem.
  m3 
  (transpose m3)
  (v= (col m3 1) (row (transpose m3) 1))
  (v= (row m3 1) (col (transpose m3) 1))
	  
  ;; Striding and Slicing issues:
  ;; Strides provide matrix sections; slicing provides vector'd sections.

  ;; STRIDING
  m01
  (strides m01 :nrows 2 :row-stride 2)  ;; view just rows 1 and 3 from m01
  (strides m01 :nrows 3) ;; first 3 rows
  (strides m01 :ncols 3 :col-stride 2) ;; cols 1, 3 ,5
  (strides m01 :ncols 2) ;; first 2 cols
  m01

  ;; SLICING
  m01
  (slice m01 :offset 5 :stride  2 :nelts 3 :type :row)
  ;; col 2 
  (slice m01 :offset 5 :stride  2 :nelts 3 :type :row)


  (slice (transpose m01) :offset 5 :stride  2 :nelts 3 :type :row)
  (slice m01
	 :offset 5
	 :stride  2
	 :nelts 3
	 :type :row)
  (slice (transpose m01) :offset 5 :stride  2 :nelts 3 :type :row)

  ;; slicing isn't affected by transposition -- doesn't affect the
  ;; counting.  Would have suggested that column-major or row-major.
  ;; Should this be the case?  (need to migrate to unit-tests).

  (v=  (slice m01 :offset 5 :stride  2 :nelts 3 :type :row)
       (slice (transpose m01) :offset 5 :stride  2 :nelts 3 :type :row))
  (v=  (slice m01 :offset 5 :stride  2 :nelts 3 :type :row)
       (slice (transpose m01) :offset 5 :stride  2 :nelts 3 :type :column))
  ;; and note the above -- vector equality doesn't depend on orientation...

  (slice m01 :offset 1 :stride  2 :nelts 3 :type :column)
  (slice m01 :offset 1 :stride  0 :nelts 3 :type :column)
  ;; :type   : provides the form to provide output for
  ;; :offset : number of observations (in "col/row major"
  ;;           matrix-dependent order) to skip over before starting
  ;;           extraction
  ;; :stride : 0 = repeat same value; 1, as ordered, 2 every other, 
  ;;           etc... 


  ;; Alternative approach for slicing, via Tamas's AFFI package:
  (defparameter *my-idx* (affi:make-affi '(5 6))) ; -> generator
  (affi:calculate-index *my-idx* #(1 1)) ; -> 7 



  ;; FIXME: need to get the foriegn-friendly arrays package involved
  ;; to create integer matrices.  Or do we just throw an error that
  ;; says to use lisp-arrays?
  (make-matrix 2 5
	       :implementation :foreign-array 
	       :element-type 'integer 
	       :initial-contents #2A(( 1 2 3 4 5)
				     ( 6 7 8 9 10)))


  ;; FIXME -- indexing with mref not checked against dims, doesn't
  ;; barf correctly.  (now is checked, but badly/poorly -- this FIXME
  ;; is about better optimization, NOT about it failing to work, which
  ;; was the original problem).
  m01
  (assert-valid-matrix-index m01 1 8)
  (assert-valid-matrix-index m01 8 1)
  (mref m01 1 8) ; good -- we throw an error... but
  (mref m01 8 1) ; BAD! barfs, not protecting against first index...
  (setf (mref m01 7 7) 1.2d0)
  m01
  
  
  ;; FIXME: the following has no applicable method -- only for
  ;; doubles, not integers.  
  (m* m2 (transpose m2))
  ;; but we can multiple doubles, but...
  (m* m01 (transpose m01))






(progn 
  (defparameter *a*
    (make-matrix 6 5 :initial-contents '((1d0 2d0 3d0 4d0 5d0)
					 (6d0  7d0  8d0  9d0  10d0)
					 (11d0 12d0 13d0 14d0 15d0)
					 (16d0 17d0 18d0 19d0 20d0)
					 (21d0 22d0 23d0 24d0 25d0)
					 (26d0 27d0 28d0 29d0 30d0))))
  (defparameter *b* (strides *a* :nrows 3 :row-stride 2))
  (defparameter *b1* (strides *a* :nrows 2 :ncols 3 :row-stride 2 :col-stride 1))
  (defparameter *c* (window *a* :nrows 3 :row-offset 3))
  (defparameter *d* (window *a* :nrows 3 :ncols 2 :row-offset 3 :col-offset 2))
  (format nil "Data initialized"))

(orientation *b*)

;; Striding
(typep *b* 'lisp-matrix::strided-matview)
(typep *b* 'lisp-matrix::window-matview)
(typep *b* 'strided-matview)
(typep *b* 'window-matview)

(parent *b*)
(offset *b*) (offset *a*)
(row-offset *a*) (col-offset *a*)
(row-offset *b*) (col-offset *b*)
(row-offset *c*) (row-offset *c*)
(col-stride *b*)  (row-stride *b*) (nrows (parent *b*))

(equal  (data *a*)
	(data *b*))
;; col 0 =  1  3  5 indicies; currently getting  1 13 25  (+ 12, not + 2)
;; col 1 =  7  9 11 indicies
;;
(m= (princ  (col *b* 0))
    (princ  (make-matrix 3 1 :initial-contents '((1d0) (11d0) (21d0)))))
(m= (col *b* 1) 
    (make-matrix 3 1 :initial-contents '((2d0) (12d0) (22d0))))
(m= (col *b* 2)
    (make-matrix 3 1 :initial-contents '((3d0) (13d0) (23d0))))
(m= (col *b* 3)
    (make-matrix 3 1 :initial-contents '((4d0) (14d0) (24d0))))
(m= (col *b* 4)
    (make-matrix 3 1 :initial-contents '((5d0) (15d0) (25d0))))

Tasks [1/3]

[#B] Migrate DITZ issues into this file.

  • State “TODO” from “” [2010-06-07 Mon 16:53]

[#B] Refactor src into “lisp-matrix”, “support”, etc.

  • State “TODO” from “” [2010-06-07 Mon 16:39]

[#B] Lisp-matrix in own package

  • State “DONE” from “CURR” [2010-06-07 Mon 16:39]
    Finished a while back.
  • State “CURR” from “TODO” [2010-06-07 Mon 16:39]
  • State “TODO” from “” [2010-06-07 Mon 16:39]

supports a lisp-matrix-user playground.

Disserata

say what.

Something went wrong with that request. Please try again.