Commit 83cd8040 authored by Liam Healy's avatar Liam Healy
Browse files

GSL array making with make-array* and #m

GSL arrays are garbage-collected objects, and so can be made like any
other lisp object; they do not need to be make inside a wrapper form.
#m is a macro which makes them; it expands to a #'make-array* call.
Examples and usage ported to this new form.  Compiled only; not tested
yet.  No effort to port to non-native implementation yet; tests as
expanded in the tests/ directory are not ported yet.
parent 5d155e1c
;; Basis splines.
;; Liam Healy 2008-02-18 14:43:20EST basis-splines.lisp
;; Time-stamp: <2008-08-23 15:35:59EDT basis-splines.lisp>
;; Time-stamp: <2008-11-30 23:43:58EST basis-splines.lisp>
;; $Id$
(in-package :gsl)
......@@ -75,13 +75,13 @@
(bw (basis-spline order nbreak))
(mw (fit-workspace ndata ncoeffs))
(rng (random-number-generator *mt19937* 0))
(B (vector-double-float ncoeffs))
(c (vector-double-float ncoeffs))
(cov (matrix-double-float (list ncoeffs ncoeffs)))
(w (vector-double-float ndata))
(x (vector-double-float ndata))
(y (vector-double-float ndata))
(Xmatrix (matrix-double-float (list ndata ncoeffs)))
(B (make-array* 'double-float :dimensions ncoeffs))
(c (make-array* 'double-float :dimensions ncoeffs))
(cov (make-array* 'double-float :dimensions (list ncoeffs ncoeffs)))
(w (make-array* 'double-float :dimensions ndata))
(x (make-array* 'double-float :dimensions ndata))
(y (make-array* 'double-float :dimensions ndata))
(Xmatrix (make-array* 'double-float :dimensions (list ndata ncoeffs)))
(sigma 0.1d0))
;; The data to be fitted.
(dotimes (i ndata)
;; Tests of array functions
;; Liam Healy 2008-10-20 22:41:48EDT array-tests.lisp
;; Time-stamp: <2008-11-16 10:28:21EST array-tests.lisp>
;; Time-stamp: <2008-11-30 22:51:15EST array-tests.lisp>
;; $Id: $
;;; Generate each file with #'write-test-to-file, e.g.
......@@ -24,8 +24,8 @@
matrix-set-all-m+ :no-complex
(letm ((m1 (array-default '(3 3) t))
(m2 (array-default '(3 3))))
(let ((m1 (array-default '(3 3) t))
(m2 (array-default '(3 3))))
(set-all m1 (scalar-default))
(cl-array (m+ m1 m2))))
;; Combinations
;; Liam Healy, Sun Mar 26 2006 - 11:51
;; Time-stamp: <2008-11-29 15:09:07EST combination.lisp>
;; Time-stamp: <2008-11-30 19:17:47EST combination.lisp>
;; $Id$
(in-package :gsl)
......@@ -29,11 +29,9 @@
(arglist-only combination "A combination." n k sync-array-on-exit)
(defun make-data-combination (nk)
(export 'make-combination)
(defun make-combination (nk)
"Make the combination object with the data array."
;; The replaces the function of make-data-from-array and
;; make-data-from-dimensions because the class needs different
;; initialization than other kinds of data.
(let ((k (second nk)))
;; "Data" is bulk arrayed data, like vectors, matrices, permutations,
;; combinations, or histograms.
;; Liam Healy 2008-04-06 21:23:41EDT
;; Time-stamp: <2008-11-29 15:08:23EST data.lisp>
;; Time-stamp: <2008-11-30 19:27:55EST data.lisp>
;; $Id$
(in-package :gsl)
......@@ -80,7 +80,6 @@
(defparameter *class-element-type* nil
;;; This is needed for make-data-from-dimensions
"The mapping between the class name and the CL element type.")
(defun data-class-name (category element-type)
......@@ -104,6 +103,7 @@
((element-type :initform ',element-type-cl
:allocation :class)))
;; Define the letm-expansion method
;;#+(or) ; REMOVE
(defmethod letm-expansion
(symbol (type (eql ',class-name)) args body)
(expand-data symbol type (first args) (second args) body))
......@@ -111,6 +111,7 @@
(pushnew ',(cons class-name element-type-cl)
*class-element-type* :test #'equal)
;; Set up dummy function to provide arglist and documentation
;;#+(or) ; REMOVE
,(format nil "A ~(~a~) of element type ~(~a~)."
......@@ -123,29 +124,47 @@
;;;; Make data from either the dimensions provided or from the initial values
(defun make-data-from-array (class array)
"Make the data object from the CL array make with make-array*."
(make-instance class
:cl-array array
:mpointer nil ; this will be set by :before method below.
#-native :c-pointer #-native nil ; this will be set by defmfun
:dimensions (array-dimensions array)
:total-size (array-total-size array)))
(defun make-data-from-dimensions (type dimensions)
(make-array* dimensions (lookup-type type *class-element-type*))))
(defun make-data (type array-or-dimensions)
(if (eq type 'combination)
(make-data-combination array-or-dimensions)
(if (typep array-or-dimensions 'array)
(make-data-from-array type array-or-dimensions)
(make-data-from-dimensions type array-or-dimensions))))
(export 'make-array*)
(defun make-array* (element-type &rest keys)
"Make a GSLL array with the given element type,
:dimensions, :initial-contents and/or :initial-element."
(let ((ffa (apply #'make-ffa element-type keys)))
(if (eql (array-rank ffa) 2) 'matrix 'vector)
:cl-array ffa
:mpointer nil
#-native :c-pointer #-native nil ; this will be set by defmfun
:dimensions (array-dimensions ffa)
:total-size (array-total-size ffa))))
(defun hashm-numeric-code (n)
"Get the appropriate element type for the numeric code n"
(case n
((nil 1) 'double-float)
(2 '(complex double-float))
(3 'single-float)
(4 '(complex single-float))
(7 '(signed-byte 8))
(8 '(unsigned-byte 8))
(15 '(signed-byte 16))
(16 '(unsigned-byte 16))
(31 '(signed-byte 32))
(32 '(unsigned-byte 32))
(63 '(signed-byte 64))
(64 '(unsigned-byte 64))))
#\# #\m
(lambda (stream subchar arg)
(declare (ignore subchar))
(read-char stream)
(let ((list (read-delimited-list #\) stream)))
`(make-array* ',(hashm-numeric-code arg) :initial-contents ',list))))
;;;; Expand letm bindings
;;;; Expand letm bindings TO BE OBSOLETE
(defun component-type (eltype)
......@@ -199,6 +218,7 @@
"Create an array (vector or matrix) from the literal contents."
`(abody ,@(mapcar (lambda (e) `(quote ,e)) type-contents)))
;;#+(or) ; REMOVE
(defun expand-data (symbol type init-or-spec sync-exit body)
"Expand the form for a gsl-data object. The symbol is bound
within the body to the object. The argument
;; Use the foreign-friendly arrays package.
;; Liam Healy 2008-03-22 15:40:08EDT
;; Time-stamp: <2008-11-29 11:16:03EST foreign-friendly.lisp>
;; Time-stamp: <2008-11-30 18:44:50EST foreign-friendly.lisp>
;; $Id$
;;; Foreign-friendly arrays (original implementation by Tamas Papp)
......@@ -14,32 +14,58 @@
;;;; Make arrays for possible foreign use
(export '(make-array*))
(defun make-array*
(dimensions element-type
&key (initial-element nil initial-element-p)
(initial-contents nil initial-contents-p))
(defun make-ffa
&key dimensions (initial-element nil initial-element-p)
(initial-contents nil initial-contents-p)
"Make an array of one or two dimensions for possible use in foreign code.
Syntax is similar to make-array, but note that element-type
is mandatory and limited to certain types."
(assert (member element-type *array-element-types* :test 'equal)
"Specified element-type must be one of *array-element-types*.")
(if initial-contents-p
(if (and (subtypep element-type 'complex)
(= (length initial-contents)
(* 2 (if (listp dimensions) (apply #'* dimensions) dimensions))))
(make-ffa dimensions element-type
(loop for (re im) on initial-contents by #'cddr
collect (complex re im)))
(make-ffa dimensions element-type :initial-contents initial-contents))
(if initial-element-p
(make-ffa dimensions element-type :initial-element initial-element)
(make-ffa dimensions element-type))))
(defun make-ffa (dimensions element-type &key
(let* ((initial-matrix (listp (first initial-contents)))
(and initial-contents-p
(subtypep element-type 'complex)
(if initial-matrix
(caar initial-contents)
(first initial-contents))
(unless dimensions
(if initial-contents-p
(setf dimensions
(if initial-matrix
;; matrix
(list (length initial-contents)
(ceiling (length (first initial-contents))
(if complex-initial-real 2 1)))
;; vector
(ceiling (length initial-contents)
(if complex-initial-real 2 1))))
(error "dimensions must be specified if contents are not")))
(when (and initial-matrix initial-contents-p)
;; flatten matrix spec
(setf initial-contents (mapcar #'append initial-contents)))
(when complex-initial-real
(setf initial-contents (list-complex-from-real initial-contents)))
(apply 'make-ffa-1d
dimensions element-type
(when initial-contents-p
(list :initial-contents initial-contents))
(when initial-element-p
(list :initial-element initial-element))))))
(defun list-complex-from-real (list)
"Collect a list of complex numbers from real components."
(assert (evenp (length list)) (list)
"Cannot collect reals into complex from an odd number of reals.")
(loop for (re im) on list by #'cddr collect (complex re im)))
(defun make-ffa-1d (dimensions element-type &key
(initial-element 0 initial-element-p)
(initial-contents nil initial-contents-p))
"Make an array that is either one-dimensional or displaced to a
......@@ -60,6 +86,9 @@
;; contents given, copy or coerce
(when (listp (first initial-contents))
;; flatten list for matrices
(setf initial-contents (apply 'append initial-contents)))
(assert (= (length initial-contents) length))
(if (typep initial-contents (list 'vector element-type))
(copy-seq initial-contents)
;; Permutations
;; Liam Healy, Sun Mar 26 2006 - 11:51
;; Time-stamp: <2008-10-25 18:52:27EDT permutation.lisp>
;; Time-stamp: <2008-11-30 23:32:48EST permutation.lisp>
;; $Id$
(in-package :gsl)
......@@ -287,7 +287,7 @@
(swap-elements perm-1 1 3)
(cl-array perm-1))
(letm ((perm-1 (permutation 4 t)) ;permute-vector
(intvec (vector-signed-byte-32 (a 11 22 33 44))))
(intvec #31m(11 22 33 44)))
(set-identity perm-1)
(swap-elements perm-1 1 3)
(swap-elements perm-1 0 2)
......@@ -84,6 +84,7 @@ will be visible; however, the function should not be used outside a
<div class="content">
<p>FIX: #m, #nm, <code>make-array*</code>
GSLL has many functions that work on vectors (one-dimensional arrays)
and matrices (two-dimensional arrays). GSLL supports all array element
......@@ -204,7 +205,7 @@ gamma-randist, negative-binomial.
<address><a href="">Liam Healy</a></address>
<!-- hhmts start -->
Time-stamp: <2008-11-15 21:53:03EST documentation.html>
Time-stamp: <2008-11-30 19:28:51EST documentation.html>
<!-- hhmts end -->
;; Eigenvectors and eigenvalues
;; Liam Healy, Sun May 21 2006 - 19:52
;; Time-stamp: <2008-10-25 11:34:13EDT eigensystems.lisp>
;; Time-stamp: <2008-11-30 23:37:31EST eigensystems.lisp>
;; $Id$
(in-package :gsl)
......@@ -158,13 +158,12 @@
(defun eigenvalue-eigenvectors-example ()
(letm ((evecs (matrix-double-float '(3 3)))
(evals (vector-double-float 3))
(letm ((evecs (make-array* 'double-float :dimensions '(3 3)))
(evals (make-array* 'double-float :dimensions 3))
(ws (eigen-symmv 3))
(mat (matrix-double-float
(a (20.0d0 -10.0d0 0.0d0)
(-10.0d0 30.0d0 0.0d0)
(0.0d0 0.0d0 40.0d0)))))
(mat #m((20.0d0 -10.0d0 0.0d0)
(-10.0d0 30.0d0 0.0d0)
(0.0d0 0.0d0 40.0d0))))
(eigenvalues-eigenvectors mat evals evecs ws)
(values (cl-array evals) (cl-array evecs))))
;; Definition of GSLL system
;; Liam Healy
;; Time-stamp: <2008-11-29 11:18:52EST gsll.asd>
;; Time-stamp: <2008-11-30 23:17:52EST gsll.asd>
;; $Id$
(asdf:defsystem "gsll"
......@@ -94,7 +94,7 @@
;; Skip fft for now, I'm not sure how it works in C
(:file "numerical-integration" :depends-on (init))
(:module random
:depends-on (init)
:depends-on (init data)
((:file "rng-types")
(:file "generators" :depends-on (rng-types))
;; Define examples.
;; Liam Healy 2008-09-07 21:00:48EDT generate-tests.lisp
;; Time-stamp: <2008-11-16 12:13:10EST generate-examples.lisp>
;; Time-stamp: <2008-11-30 23:02:06EST generate-examples.lisp>
;; $Id: $
;;; Define examples that can be displayed by users with the function
......@@ -102,27 +102,28 @@
;;;; Generate forms for all array types
(defun array-default (spec &optional no-init sync-on-exit)
(defun array-default (spec &optional no-init)
"Make an array of the current type and intialize from the pool."
(declare (special default-element-type starting-element))
(let ((matrixp (listp spec)))
`(,(data-class-name (if matrixp 'matrix 'vector) default-element-type)
,(if no-init
`',spec ; no initial values, just dimension
(cons ; macro to make initial-contents
(if matrixp
(loop repeat (first spec)
(if no-init
;; No initial values, just dimension
`(make-array* ',default-element-type :dimensions ',spec)
;; Initial contents
',(if matrixp
(loop repeat (first spec)
default-element-type (second spec) starting-element)
(incf starting-element (second spec)))
default-element-type (second spec) starting-element)
(incf starting-element (second spec)))
default-element-type spec starting-element)
(incf starting-element spec)))))
default-element-type spec starting-element)
(incf starting-element spec)))))))
(defun scalar-default (&optional float-type)
"Make a scalar of the current type from the pool. For complex
......@@ -162,7 +163,7 @@
(setf starting-element 0)
'(array-default vector-default scalar-default)))))
'(array-default scalar-default)))))
(defmacro generate-all-array-tests (name element-types test)
`(save-test ,name ,@(generate-all-array-tests-body element-types test)))
;; Example spline
;; Liam Healy, Sat Nov 10 2007 - 21:18
;; Time-stamp: <2008-03-09 19:21:47EDT spline-example.lisp>
;; Time-stamp: <2008-11-30 23:45:40EST spline-example.lisp>
;; $Id$
(in-package :gsl)
......@@ -20,8 +20,8 @@
(multiple-value-bind (xarr yarr)
(letm ((acc (acceleration))
(cxarr (vector-double-float xarr))
(cyarr (vector-double-float yarr))
(cxarr (make-array* 'double-float :dimensions xarr))
(cyarr (make-array* 'double-float :dimensions yarr))
(spline (spline *cubic-spline-interpolation* cxarr cyarr)))
(loop for xi from (aref xarr 0) below (aref xarr 9) by 0.01d0
collect (list xi (evaluate-spline spline xi acc))))))
;; Exponential of a matrix
;; Liam Healy 2008-08-10 17:25:35EDT exponential.lisp
;; Time-stamp: <2008-08-10 18:27:04EDT exponential.lisp>
;; Time-stamp: <2008-11-30 23:35:44EST exponential.lisp>
;; $Id: $
(in-package :gsl)
......@@ -25,8 +25,8 @@
;;; A matrix of the form [[0, 1], [-1, 0]]
;;; when exponentiated gives [[cos x, sin x], [-sin x, cos x]]
(letm ((mat (matrix-double-float (a (0.0d0 1.0d0) (-1.0d0 0.0d0))))
(exp (matrix-double-float '(2 2))))
(letm ((mat #m((0.0d0 1.0d0) (-1.0d0 0.0d0)))
(exp (make-array* 'double-float :dimensions '(2 2))))
(cl-array (matrix-exponential mat exp)))
#2A((0.5403023058681384d0 0.841470984807895d0)
(-0.841470984807895d0 0.5403023058681384d0))
;; LU decomposition
;; Liam Healy, Thu Apr 27 2006 - 12:42
;; Time-stamp: <2008-11-29 15:30:45EST lu.lisp>
;; Time-stamp: <2008-11-30 23:44:35EST lu.lisp>
;; $Id$
(in-package :gsl)
......@@ -136,10 +136,10 @@
(export 'invert-matrix)
(defun invert-matrix (mat)
"Invert the matrix."
(letm ((mmat (matrix-double-float mat))
(letm ((mmat mat)
(dim (array-dimension mat 0))
(per (permutation dim))
(inv (matrix-double-float (list dim dim))))
(inv (make-array* 'double-float :dimensions (list dim dim))))
(LU-decomposition mmat per)
(lu-invert mmat per inv)
(cl-array inv)))
;; Monte Carlo Integration
;; Liam Healy Sat Feb 3 2007 - 17:42
;; Time-stamp: <2008-10-25 11:44:35EDT monte-carlo.lisp>
;; Time-stamp: <2008-11-30 23:04:24EST monte-carlo.lisp>
;; $Id$
(in-package :gsl)
......@@ -294,22 +294,22 @@
(defun random-walk-plain-example (&optional (nsamples 500000))
(letm ((ws (monte-carlo-plain 3))
(lower (vector-double-float (a 0.0d0 0.0d0 0.0d0)))
(upper (vector-double-float (a* pi pi pi)))
(lower #m(0.0d0 0.0d0 0.0d0))
(upper (make-array* 'double-float :initial-contents (list pi pi pi)))
(rng (random-number-generator *mt19937* 0)))
(monte-carlo-integrate-plain monte-carlo-g lower upper nsamples rng ws)))
(defun random-walk-miser-example (&optional (nsamples 500000))
(letm ((ws (monte-carlo-miser 3))
(lower (vector-double-float (a 0.0d0 0.0d0 0.0d0)))
(upper (vector-double-float (a* pi pi pi)))
(lower #m(0.0d0 0.0d0 0.0d0))
(upper (make-array* 'double-float :initial-contents (list pi pi pi)))
(rng (random-number-generator *mt19937* 0)))
(monte-carlo-integrate-miser monte-carlo-g lower upper nsamples rng ws)))
(defun random-walk-vegas-example (&optional (nsamples 500000))
(letm ((ws (monte-carlo-vegas 3))
(lower (vector-double-float (a 0.0d0 0.0d0 0.0d0)))
(upper (vector-double-float (a* pi pi pi)))
(lower #m(0.0d0 0.0d0 0.0d0))
(upper (make-array* 'double-float :initial-contents (list pi pi pi)))
(rng (random-number-generator *mt19937* 0)))
(monte-carlo-integrate-vegas monte-carlo-g lower upper nsamples rng ws)))
;; Polynomials
;; Liam Healy, Tue Mar 21 2006 - 18:33
;; Time-stamp: <2008-10-25 11:47:01EDT polynomial.lisp>
;; Time-stamp: <2008-11-30 23:38:39EST polynomial.lisp>
;; $Id$
(in-package :gsl)
......@@ -146,11 +146,11 @@
w of the appropriate size. The n-1 roots are returned in
the packed complex array z of length 2(n-1), alternating
real and imaginary parts."
(let ((len (length coefficients)))
(letm ((coef (vector-double-float coefficients))
(answer (vector-double-float (* 2 (1- len))))
(let ((len (total-size coefficients)))
;; Should this be making a complex array?
(letm ((answer (make-array* 'double-float :dimensions (* 2 (1- len))))
(ws (complex-workspace len)))
(values-list (polynomial-solve-ws coef ws answer)))))
(values-list (polynomial-solve-ws coefficients ws answer)))))
(defmfun polynomial-solve-ws (coefficients workspace answer-pd)
......@@ -172,16 +172,16 @@
(save-test polynomial
(letm ((xa (vector-double-float (a 0.0d0 1.0d0 2.0d0 3.0d0)))
(ya (vector-double-float (a 2.5d0 7.2d0 32.7d0 91.0d0)))
(dd (vector-double-float 4)))
(letm ((xa #m(0.0d0 1.0d0 2.0d0 3.0d0))
(ya (#m(2.5d0 7.2d0 32.7d0 91.0d0) ))
(dd (make-array* 'double-float :dimensions 4)))
(divided-difference dd xa ya)
(polynomial-eval-divided-difference dd xa 0.0d0)
(polynomial-eval-divided-difference dd xa 1.0d0)
(polynomial-eval-divided-difference dd xa 2.0d0)
(polynomial-eval-divided-difference dd xa 3.0d0)))
(letm ((vec (vector-double-float (a 1.0d0 2.0d0 3.0d0))))
(letm ((vec #m(1.0d0 2.0d0 3.0d0)))
(polynomial-eval vec -1.0d0))
(solve-quadratic 1.0d0 0.0d0 1.0d0)
(solve-quadratic 1.0d0 -2.0d0 1.0d0)
......@@ -189,5 +189,5 @@
(solve-cubic -6.0d0 -13.0d0 42.0d0)
(solve-cubic-complex -1.0d0 1.0d0 -1.0d0)
;; Example from GSL manual
(polynomial-solve (a double-float -1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0)))
(polynomial-solve #m(-1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0)))
;; Dirichlet distribution
;; Liam Healy, Sun Oct 29 2006
;; Time-stamp: <2008-10-25 12:13:42EDT dirichlet.lisp>
;; Time-stamp: <2008-11-30 23:17:06EST dirichlet.lisp>
;; $Id$
(in-package :gsl)
......@@ -57,15 +57,15 @@
;;; Examples and unit test
(save-test dirichlet
(letm ((rng (random-number-generator *mt19937* 0))
(alpha (vector-double-float (a 1.0d0 2.0d0 3.0d0 4.0d0)))
(theta (vector-double-float 4)))
(alpha #m(1.0d0 2.0d0 3.0d0 4.0d0))
(theta (make-array* 'double-float :dimensions 4)))
(dirichlet rng alpha theta)
(cl-array theta))
(letm ((alpha (vector-double-float (a 1.0d0 2.0d0 3.0d0 4.0d0)))
(theta (vector-double-float (a 0.1d0 0.3d0 0.4d0 0.2d0))))
(letm ((alpha #m(1.0d0 2.0d0 3.0d0 4.0d0))
(theta #m(1.0d0 2.0d0 3.0d0 4.0d0)))
(dirichlet-pdf alpha theta))
(letm ((alpha (vector-double-float (a 1.0d0 2.0d0 3.0d0 4.0d0)))
(theta (vector-double-float (a 0.1d0 0.3d0 0.4d0 0.2d0))))
(letm ((alpha #m(1.0d0 2.0d0 3.0d0 4.0d0))
(theta #m(1.0d0 2.0d0 3.0d0 4.0d0)))
(dirichlet-log-pdf alpha theta)))
;; Discrete random variables
;; Liam Healy, Sat Nov 11 2006 - 21:51
;; Time-stamp: <2008-10-25 12:14:54EDT discrete.lisp>
;; Time-stamp: <2008-11-30 23:13:24EST discrete.lisp>
;; $Id$
(in-package :gsl)
......@@ -53,15 +53,12 @@
(save-test discrete
;; Must have two letms because the vector value is not set until
;; the body, but the discrete-random needs that set value.
(letm ((probabilities (vector-double-float (a 0.25d0 0.5d0 0.25d0))))
(letm ((probabilities #m(0.25d0 0.5d0 0.25d0)))
(letm ((table (discrete-random probabilities))
(rng (random-number-generator *mt19937* 0)))
(loop for i from 0 to 10