marray.lisp 10.4 KB
Newer Older
1
;; A "marray" is an array in both GSL and CL
2
;; Liam Healy 2008-04-06 21:23:41EDT
3
;; Time-stamp: <2010-01-20 15:31:57EST marray.lisp>
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
;;
;; Copyright 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
;;
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program.  If not, see <http://www.gnu.org/licenses/>.
20
21
22
23

(in-package :gsl)

;;;;****************************************************************************
Liam Healy's avatar
Liam Healy committed
24
;;;; The class marray and its construction
25
26
;;;;****************************************************************************

Liam Healy's avatar
Liam Healy committed
27
(export 'marray)
28
(defclass marray (mobject c-array:foreign-array)
29
30
  ((block-pointer :initarg :block-pointer :reader block-pointer)
   (total-size :reader size))
31
32
  (:documentation
   "A superclass for arrays represented in GSL and CL."))
33

Liam Healy's avatar
Liam Healy committed
34
;;; We don't allocate or free the C array data, because that is
35
;;; handled by c-array:foreign-array.  We can use the GSL functions
Liam Healy's avatar
Liam Healy committed
36
37
38
39
40
41
;;; *_alloc_from_block because they will allocate only the structure,
;;; but we must use CFFI to allocate the block structure; otherwise
;;; GSL would try to allocate the C array.  We do not use any of the
;;; GSL free functions because they would free the C array data.

(defmethod initialize-instance :after ((object marray) &rest initargs)
42
  (declare (ignore initargs)) ; required by &rest arg for c-array:foreign-array?
Liam Healy's avatar
Liam Healy committed
43
44
45
46
47
48
  ;; Need to check that all allocations succeeded.
  ;; Don't do anything if mpointer has been assigned (shouldn't happen)
  ;; or this is a permutation or combination (they have their own methods).
  (unless (and (slot-boundp object 'mpointer) (mpointer object))
    (let ((blockptr (cffi:foreign-alloc 'gsl-block-c)))
      (setf (cffi:foreign-slot-value blockptr 'gsl-block-c 'size)
49
	    (size object)
50
51
	    (slot-value object 'block-pointer)
	    blockptr)
Liam Healy's avatar
Liam Healy committed
52
53
      (let ((array-struct (alloc-from-block object blockptr)))
	(setf (slot-value object 'mpointer) array-struct)
54
	#-native (set-struct-array-pointer object)
Liam Healy's avatar
Liam Healy committed
55
56
57
58
	(tg:finalize object
		     (lambda ()
		       (cffi:foreign-free blockptr)
		       (cffi:foreign-free array-struct)))))))
59

60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
(defun set-struct-array-pointer (object)
  "Set the pointer in the 'data slot of the foreign structs to be the
   current c-pointer value.  In non-native implementations this need
   be called only once when the marray is made.  In native implementations
   it is called whenever mpointer is requested because of the possibility
   that a GC moved the pointer."
  (setf (cffi:foreign-slot-value (block-pointer object) 'gsl-block-c 'data)
	(c-pointer object)
	;; alloc-from-block automatically copies over the data pointer
	;; from the block to the vector/matrix; we must do that manually here
	(cffi:foreign-slot-value
	 (slot-value object 'mpointer)
	 (if (typep object 'matrix) 'gsl-matrix-c 'gsl-vector-c)
	 'data)
	(c-pointer object)))

#+native
(defmethod mpointer ((object marray))
  "Compute the c-pointer of the array and place it in the GSL struct
   because the stored version is untrustworthy unless it was computed
   within the same native-pointer-protect form as this mpointer
   extraction."
  (set-struct-array-pointer object)
Liam Healy's avatar
Liam Healy committed
83
  (call-next-method))
84

85
(defmethod make-load-form ((object marray) &optional env)
86
  (declare (ignore env))
87
  `(make-marray
88
    ',(element-type object)
89
    :dimensions ',(dimensions object)
90
91
    :initial-contents
    ',(loop for elt across
Liam Healy's avatar
Liam Healy committed
92
93
94
	   (subseq (c-array:original-array object)
		   (c-array:offset object)
		   (+ (c-array:offset object) (size object)))
95
96
	   collect elt)))

97
98
99
100
101
102
103
104
105
106
107
;;;;****************************************************************************
;;;; Definition of specific data classes
;;;;****************************************************************************

(defparameter *class-element-type* nil
  "The mapping between the class name and the CL element type.")

(defun data-class-name (category element-type)
  "The class name from the type of element."
  ;; e.g. (data-class-name 'vector '(unsigned-byte 8))
  ;; -> VECTOR-UNSIGNED-BYTE-8
108
  (if (member category '(vector matrix))
109
110
111
      (intern
       (format nil "~a-~a" category (c-array:cl-single element-type :gsl))
       :gsl)
112
      category))
113
114
115
116
117
118
119
120
121
122
123
124

(defun data-defclass (category superclass)
  "Define all the subclasses based on the known element types."
  (cons 'progn
	(mapcan
	 (lambda (element-type-cl)
	   (let* ((class-name (data-class-name category element-type-cl)))
	     ;; Define the class
	     `((defclass ,class-name
		   (,superclass)
		 ((element-type :initform ',element-type-cl
				:allocation :class)))
125
126
	       ;; Push mapping onto *class-element-type*
	       (pushnew ',(cons class-name element-type-cl)
127
128
			*class-element-type* :test #'equal)
	       (export ',class-name))))
129
	 c-array:*array-element-types*)))
130
131
132
133
134

;;;;****************************************************************************
;;;; Make data from either the dimensions provided or from the initial values
;;;;****************************************************************************

135
136
(export 'make-marray)
(defun make-marray
137
    (class-or-element-type &rest keys &key dimensions initial-contents cl-array
138
     &allow-other-keys)
139
  "Make a GSLL array with the given element type,
140
   :dimensions, and :initial-contents, :initial-element or :cl-array.
141
142
143
144
145
146
147
148
149
150
151
152
153
154
   If the :cl-array is supplied, it should be a CL array generated
   with #'make-ffa."
  (apply #'make-instance
	 (if (subtypep class-or-element-type 'marray)
	     class-or-element-type
	     (data-class-name
	      (if
	       (or
		(and dimensions (listp dimensions) (eql (length dimensions) 2))
		(and initial-contents (listp (first initial-contents)))
		(and cl-array (eql (length (array-dimensions cl-array)) 2)))
	       'matrix 'vector)
	      class-or-element-type))
	 keys))
155

156
157
(eval-when (:compile-toplevel :load-toplevel :execute)
  (pushnew 'marray grid:*grid-data-superclasses*))
Liam Healy's avatar
Liam Healy committed
158
159
160

(defmethod grid:make-grid-data
    ((type (eql 'marray)) dimensions rest-spec
161
162
     &key initial-element)
  (make-marray (grid:spec-scalar-p rest-spec)
Liam Healy's avatar
Liam Healy committed
163
164
165
	       :dimensions dimensions
	       :initial-element initial-element))

166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
;;; The reader macro #m will read a list of arguments, evaluating the
;;; contents, and construct a marray from it.  If the symbol ^ occurs
;;; in the list, it indicates the object being made is a matrix, and
;;; this symbol indicates the end of each row.  For compatibility
;;; reasons, a list of numeric values may be given for each row also.

;;; Examples:
;;; #m(1 2 (exp 1))
;;; #<VECTOR-DOUBLE-FLOAT #(1.0d0 2.0d0 2.7182817459106445d0)>
;;; #m(1 2 3 ^ pi (cos (/ pi 4)) -12)
;;; #<MATRIX-DOUBLE-FLOAT #2A((1.0d0 2.0d0 3.0d0)
;;;                       (3.141592653589793d0 0.7071067811865476d0 -12.0d0))>
;;; #2m(2 1 3 -1)
;;; #<VECTOR-COMPLEX-DOUBLE-FLOAT #(#C(2.0d0 1.0d0) #C(3.0d0 -1.0d0))>

181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
(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))))

197
198
199
200
201
202
(defun quote-numeric-list (s)
  (if (and (listp s) (numberp (first s)))
      `',s
      s))

(defvar *row-separator* '^)
203
(export *row-separator*)
204
205
206

(defun conslist (l) (cons 'list l))

207
208
209
210
211
212
(set-dispatch-macro-character
 #\# #\m
 (lambda (stream subchar arg)
   (declare (ignore subchar))
   (read-char stream)
   (let ((list (read-delimited-list #\) stream)))
213
214
215
216
217
218
219
220
     `(make-marray
       ',(hashm-numeric-code arg)
       :initial-contents
       ,(if (find *row-separator* list)
	    (conslist
	     (mapcar 'conslist
		     (cl-utilities:split-sequence *row-separator* list)))
	    `(list ,@(mapcar 'quote-numeric-list list)))))))
221

222
;;;;****************************************************************************
223
;;;; Copy to and from bare mpointers 
224
225
;;;;****************************************************************************

226
227
228
229
230
(defgeneric contents-from-pointer (pointer struct-type &optional element-type)
  (:documentation
   "Create a contents list from the GSL object of type struct-type
    referenced by pointer."))

231
(defmethod c-array:copy-making-destination ((pointer #.+foreign-pointer-class+))
232
233
234
235
  (foreign-pointer-method
   pointer
   ;; Default assumption when destination isn't given in #'copy is
   ;; that this should make a vector-double-float.
Liam Healy's avatar
Liam Healy committed
236
   (c-array:copy-to-destination pointer 'vector-double-float)))
237

Liam Healy's avatar
Liam Healy committed
238
(defmethod c-array:copy-to-destination
239
240
241
242
243
244
245
246
247
    ((pointer #.+foreign-pointer-class+) (class-name symbol))
  (foreign-pointer-method
   pointer
   (make-marray
    class-name
    :initial-contents
    (contents-from-pointer
     pointer
     (if (subtypep class-name 'matrix) 'gsl-matrix-c 'gsl-vector-c)
Liam Healy's avatar
Liam Healy committed
248
     (c-array:lookup-type class-name *class-element-type*)))))
249
250
251
252
253
254

;; Some functions in solve-minimize-fit return a pointer to a GSL
;; vector of double-floats.  This function turns that into a
;; foreign-friendly array.  There is no choice but to copy over the
;; data even on native implementations; because GSL is doing the
;; mallocing, the data are not CL-accessible.
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275


;;;;****************************************************************************
;;;; Convenient defaulting of defmfun arguments 
;;;;****************************************************************************

(defun make-marray-or-default
    (default dimensions
     &optional dont-make-if-nil (element-type 'double-float) initial-element)
  "If default is T or dont-make-if-nil, make a marray and
   returned with the dimensions and element-type.  If default is a
   marray, it is returned after being synced for non-native."
  (if (member default (list t dont-make-if-nil))
      (if initial-element
	  (make-marray
	   element-type :dimensions dimensions :initial-element initial-element)
	  (make-marray element-type :dimensions dimensions))
      (progn #-native
	     (when (typep default 'marray)
	       (c-array:copy-cl-to-c default))
	     default)))