Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse files

doc cleanup

  • Loading branch information...
commit 092e5b07185167a4436917f161986dbfa1648328 1 parent 1e9b51e
@blindglobe authored
Showing with 976 additions and 964 deletions.
  1. +975 −960 README.org
  2. +1 −4 src/package.lisp
View
1,935 README.org
@@ -5,7 +5,7 @@
Lisp-Matrix intends to be a reasonably modern, flexible system for
numeric algebra in Lisp.
- Rif, Tomas Papp, and Mark H contributed various critical pieces
+ Rif, Tamas Papp, and Mark H contributed various critical pieces
(code, design, code, testing, code, benchmarks, code, and more
code).
@@ -69,40 +69,57 @@
** Packages:
*** From Rif:
+
**** org.middleangle.cl-blapack
+
+ This covers the overall BLAS and LAPACK integration.
+
**** org.middleangle.foreign-numeric-vector
+ This provides a data store for C/FORTRAN code that can be
+ accessed by common-lisp.
+
*** From Tomas:
+
**** ffa
+
**** array-operations
-From Mark:
+*** From Mark:
**** lisp-matrix integration (also found in ffa and cl-blapack).
In general:
+*** From Others
+
+ These are dependencies that arise from the original authors.
+
**** cffi (depends on babel, alexandria )
**** cl-utilities
**** iterate
**** metabang-bind
**** asdf-system-connections
-**** lift (depends on trivial-timeout)
+**** lift
+
+ (depends on trivial-timeout)
-These dependencies need to be better worked out.
+ Am not sure that this is still the case (trivial-timeout dependency)
*** Others (AJR?)
**** lift (unit testing)
**** xarray (generic array-like accessors)
+ This is not so much a dependency as an enhancer.
* Documentation
-#+srcname: loadit
-#+begin_src lisp
+#+NAME: loadit
+#+BEGIN_SRC lisp
(asdf:oos 'asdf:load-op :lisp-matrix)
-#+end_src
+ ;; (ql:quickload :lisp-matrix)
+#+END_SRC
Need to autogenerate approach for documenting what we can do with
this. Until then, simple reference.
@@ -112,10 +129,11 @@ These dependencies need to be better worked out.
- single/double/complex-single/complex-double/integer
- (TODO: need to consider normal or mmap'd structures as well)
by:
-#+srcname:
-#+begin_src lisp
+
+#+NAME: matrixExData
+#+BEGIN_SRC lisp
(make-matrix )
-#+end_src
+#+END_SRC
right now, we are being numerical analysts, and only allow for a
single modality, i.e. lisp-integer, foriegn-doubleFloat, etc.
@@ -124,19 +142,22 @@ These dependencies need to be better worked out.
typed matrices/arrays.
Referencing elements is done using the xarray system, so that needs
- to be a dependency of this.
+ to be a dependency of this. (one can use the native system, but it
+ would be so much better to have a uniform table-access and
+ manipulation API, xarray or grid or affi or...)
-#+srcname:
+#+name: NewMetaAPI
#+begin_src lisp
- (xref mat
- x
- y
+ (xref mat x y
:return-as 'matrix) ; for a single mat[x,y] value
(xref mat
(rows x1 x2 x3)
(columns y1 y2 y3)) ; for a 3x3 matrix restricted
; to the appropriate rows and
- ; columns
+ ; columns. A better approach
+ ; might be to use a
+ ; cross-product API, or
+ ; serial-list-then-row-or-column-major-fill-to-spec-d-format
(xref mat
(except-for-rows x1 x2 x3)
@@ -153,1119 +174,1113 @@ These dependencies need to be better worked out.
#+end_src
-
-
+#+NAME: OldNativeAPI
#+begin_src lisp
-
-
(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
#+end_src
+
* Usage
-** Demo (broken things)
+** Demo (working things)
+ Demos for Lisp Matrix (encoded within progn's)
+
+ 1. instantiating matrices and vectors
+ 2. inversion using BLAS/LAPACK
+
#+begin_src common-lisp
-;;; Precursor systems
-(in-package :cl-user)
-;; (asdf:oos 'asdf:compile-op 'ffa :force t)
-;; (asdf:oos 'asdf:compile-op 'array-operations :force t)
+ (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)))
+#+end_src
-;; (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
+TODO.lisp : things that don't work but should
+lm-demo.lisp : things that might work but should
+** Demo (getting started)
+#+begin_src common-lisp
+ (in-package :cl-user)
+ (asdf:oos 'asdf:load-op :lisp-matrix)
+#+end_src
-;;; The main thing...
-;; (delete-package 'lisp-matrix) ;; fails, but we need to cleanup a bit more.
+** Demo (more working things)
+#+begin_src common-lisp
+;;; 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)
-;; (asdf:oos 'asdf:load-op 'lisp-matrix)
-;; (asdf:oos 'asdf:compile-op 'cffi :force t)
+;;; 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)
-;; 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)
+
+;; if the above describes errors, here is how we figure out what bug
+;; report to write...
+
(describe (run-lisp-matrix-tests))
-;; failures:
+;;; Now we can use it, either by importing the symbols into the
+;;; current package by:
-;; 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.
+;; (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))
-(describe
- (lift::run-test
- :test-case 'lisp-matrix-unittests::strided-matrix-column-access
- :suite 'lisp-matrix-ut-vectors))
-
+;;; We wrap these up into a progn for simple overall evaluation, but
+;;; stepping through them is fine as well.
-;; Here is what we need to fix, based on the above:
-;; # creation of foreign-array matrices which are integer valued
-;; fails.
+(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*
-;; 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
+ ;; increase complexity
-;;; FIXME FOLLOWING ERRORS: MIGRATE INTO UNITTESTS...
+ (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")
-(progn ;;#FIXME: writing out R matrices -- as strings and via RCLG
+ (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...")
- (defparameter *x-temp*
- (make-matrix 4 5
- :implementation :lisp-array
+ ;; 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((11d0 12d0 13d0 14d0 15d0)
- (21d0 22d0 23d0 24d0 25d0)
- (31d0 32d0 33d0 34d0 35d0)
- (41d0 42d0 43d0 44d0 45d0))))
+ :initial-contents #2A(( 1d0 2d0 3d0 4d0 5d0)
+ ( 6d0 7d0 8d0 9d0 10d0)))
+ "placeholder 2")
+ *m2b*
- ;; 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))
+ (mref *m2b* 0 2) ;; => 3
+ *m2b*
+ (transpose *m2b*)
- (reduce #'(lambda (x y) (concatenate 'string x y))
- "test"
- " "
- (list "a2" " s3 " "asdf")
- "end.")
+ ;; 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
- (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)))))
+ (m= (col *m2b* 0)
+ (row (transpose *m2b*) 0))
+ (v= (col *m2b* 0)
+ (row (transpose *m2b*) 0))
- (lispmatrix2R *x-temp*)
+ (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")
- (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))
+ (row *m3* 2)
+ (col *m3* 1)
- )
+ (= (mref *m3* 0 1)
+ (mref (transpose *m3*) 1 0))
-#+nil
-(progn ;; QR decomp
+ (= (mref *m3* 2 2)
+ (mref (transpose *m3*) 2 2))
- (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))))
-
- )
+ *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.
-#+nil
-(progn ;; FIXME: R's apply across array indicies
+ ;;; 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.
- ;; 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).
+ (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))))
- ;; 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).
+ (v= (row *m4* 0) (col (transpose *m4*) 0))
+ (v= (col *m4* 0) (row (transpose *m4*) 0))
- ;; 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.
+ *m4*
+ (row *m4* 0)
+ (col *m4* 4)
- ;; 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))
+ (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))
- (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.
+;;; EXAMPLES TO DEMONSTRATE
- ;; 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)))))
+;;; 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*)
- (list-of-columns *m1-app*)
- (list-of-rows *m1-app*)
-
- (mapcar #'princ (list-of-columns *m1-app*))
+;; convert between foriegn-array and lisp-array.
- (format nil "R-Apply approach"))
+;; operate ()
+;; do some blas/lapack
-#+nil
-(progn
- ;; Studies in Class inheritance
+;; output
- (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)
+;; 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))))
- (subtypep 'FA-SIMPLE-VECTOR-DOUBLE 'MATRIX-LIKE)
+(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))))
- ;;; weird!
- (m- (make-vector 2 :initial-contents '((1d0 1d0)))
- (make-vector 2 :initial-contents '((1d0 1d0))))
+;; 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))))
- (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.
+;; 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...
- (m- *x-temp* *x-temp*))
+;; Windowed matrix
+(orientation m01c)
+(row m01c 0) ; Y
+(row m01c 1) ; Y
+(col m01c 0) ; Y
+(col m01c 1) ; Y
+(col m01c 2) ; Y
-#+end_src
+;; 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)
-** Demo (working things)
- Demos for Lisp Matrix (encoded within progn's)
+(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)
- 1. instantiating matrices and vectors
- 2. inversion using BLAS/LAPACK
-
-#+begin_src common-lisp
+;; 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
-(in-package :lisp-matrix-user)
+ (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)
-(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")
+
+ ;; FIXME: there are bugs in slicing/striding with transposed
+ ;; matrices.
- ;; 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")
+ ;; the following are correct, but..
+ (row m01 0)
+ (row m01 1)
+ (row m01 2)
+ (row m01 3)
- (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")
+ (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
- (defparameter *m3-fa*
- (make-matrix
- 2 2
- :implementation :foreign-array
- :element-type 'double-float
- :initial-contents #2A(( 1d0 2d0 )
- ( 6d0 7d0 )))
- "placeholder 2")
+ (row (transpose m01) 0)
+ (transpose (row (transpose m01) 0))
- (defparameter *m3-la*
- (make-matrix
- 2 2
- :implementation :lisp-array
- :element-type 'double-float
- :initial-contents #2A(( 1d0 2d0 )
- ( 6d0 7d0 )))
- "placeholder 2")
+ 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
-
- (defparameter *m01b*
- (strides *m01* :nrows 2 :ncols 3
- :row-stride 2
- :row-offset 1 :col-offset 1))
+
+ (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...
+
- (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"))
+ ;; 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.
-#+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.
- )
+ ;; 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)
-#+nil
-(progn
- ;; factorization and inversion via LAPACK
- ;; LU
- (let ((test-eye (eye 7 7)))
- (m* test-eye (minv-lu test-eye)))
+ (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)
- ;; 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))))))
+ ;; 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...
-(progn
- ;; Using xGEQRF routines for supporting linear regression.
+ (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...
- ;; 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).
+ ;; 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
- (let ((a (rand 10 5)))
- (geqrf a)))
-#+end_src
-TODO.lisp : things that don't work but should
-lm-demo.lisp : things that might work but should
-** Demo (getting started)
-#+begin_src common-lisp
- (in-package :cl-user)
- (asdf:oos 'asdf:load-op :lisp-matrix)
-#+end_src
-** Demo (more working things)
-#+begin_src common-lisp
-;;; 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.
+ ;; 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)))
-;;; = 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)
+ ;; 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))
-;;; 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*
+ (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)))
- "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)
-
+ (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"))
- (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)))
+(orientation *b*)
- (m+ (row m3 1) (row m3 2))
- (m- (row m3 1) (row m3 2))
+;; 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))))
+#+end_src
+** Demo (broken things)
+#+begin_src common-lisp
+;;; Precursor systems
+(in-package :cl-user)
+;; (asdf:oos 'asdf:compile-op 'ffa :force t)
+;; (asdf:oos 'asdf:compile-op 'array-operations :force t)
-;;; EXAMPLES TO DEMONSTRATE
+;; (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.
-;;; 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*)
+;; (asdf:oos 'asdf:compile-op 'lisp-matrix :force t)
+;; (asdf:oos 'asdf:compile-op 'lisp-matrix)
+;; (asdf:oos 'asdf:load-op 'lisp-matrix)
-;; convert between foriegn-array and lisp-array.
+;; (asdf:oos 'asdf:compile-op 'cffi :force t)
-;; operate ()
+(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))
-;; do some blas/lapack
+;; failures:
-;; output
+;; 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.
-;; 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))))
+(in-package :lisp-matrix-user)
-(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))))
+;; (lisp-matrix-unittests:run-lisp-matrix-tests)
+;; (describe (lisp-matrix-unittests:run-lisp-matrix-tests))
-;; 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))))
+(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.
-;; 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...
+;; 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
-;; Windowed matrix
-(orientation m01c)
-(row m01c 0) ; Y
-(row m01c 1) ; Y
-(col m01c 0) ; Y
-(col m01c 1) ; Y
-(col m01c 2) ; Y
+;;; FIXME FOLLOWING ERRORS: MIGRATE INTO UNITTESTS...
-;; 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)
+(progn ;;#FIXME: writing out R matrices -- as strings and via RCLG
-(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)
+ (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))))
-;; 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
+ ;; 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))
- (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)
+ (reduce #'(lambda (x y) (concatenate 'string x y))
+ "test"
+ " "
+ (list "a2" " s3 " "asdf")
+ "end.")
-
- ;; FIXME: there are bugs in slicing/striding with transposed
- ;; matrices.
+ (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)))))
- ;; the following are correct, but..
- (row m01 0)
- (row m01 1)
- (row m01 2)
- (row m01 3)
+ (lispmatrix2R *x-temp*)
- (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
+ (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))
- (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
+#+nil
+(progn ;; QR decomp
- (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.
+ (let* ((state1 (make-random-state))
+ (state2 (make-random-state state1)))
+ (m= (rand 2 3 :state state1)
+ (rand 2 3 :state state2)))
- ;; 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
+ ;;; 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))))
- ;; 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)
+#+nil
+(progn ;; FIXME: R's apply across array indicies
- ;; 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).
+ ;; 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).
- (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...
+ ;; 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).
- (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...
+ ;; 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))
- ;; 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
+ ;; 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)
- ;; 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)))
+ (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))))
- ;; 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
-
+ ;; 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*)
- ;; 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))
+ (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))))
-(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"))
+ (let ((*default-implementation* :foreign-array))
+ (m- (make-vector 2 :initial-contents '((1d0 1d0)))
+ (make-vector 2 :initial-contents '((1d0 1d0)))))
-(orientation *b*)
+ (let ((*default-implementation* :lisp-array))
+ (m- (make-vector 2 :initial-contents '((1d0 1d0)))
+ (make-vector 2 :initial-contents '((1d0 1d0)))))
-;; Striding
-(typep *b* 'lisp-matrix::strided-matview)
-(typep *b* 'lisp-matrix::window-matview)
-(typep *b* 'strided-matview)
-(typep *b* 'window-matview)
+ (m- (make-vector 2
+ :implementation :lisp-array
+ :initial-contents '((1d0 1d0)))
+ (make-vector 2
+ :implementation :foreign-array
+ :initial-contents '((1d0 1d0))))
-(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*))
+ (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.
-(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))))
+ (m- *x-temp* *x-temp*))
#+end_src
+
* Tasks [1/3]
** TODO [#B] Migrate DITZ issues into this file.
- State "TODO" from "" [2010-06-07 Mon 16:53]
View
5 src/package.lisp
@@ -24,7 +24,6 @@
:org.middleangle.foreign-numeric-vector
:org.middleangle.cl-blapack
:ffa)
- (:import-from :fnv) ;; do we really need this? We are using it!
(:export
;; base classes (need we export more?) and informtion
@@ -135,15 +134,13 @@
potrf potri potrs minv-cholesky msolve-cholesky ;; cholesky, symm matrices
geqrf ;; qr
-
matrix-like-symmetric-p
cross-product ;; outer-product
;; data transforms
list->vector-like vector-like->list
- trap2mat
- ))
+ trap2mat ))
(defpackage :lisp-matrix-user
Please sign in to comment.
Something went wrong with that request. Please try again.