marray.lisp 8.65 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: <2009-12-22 23:17:22EST marray.lisp>
4
5
6
7

(in-package :gsl)

;;;;****************************************************************************
Liam Healy's avatar
Liam Healy committed
8
;;;; The class marray and its construction
9
10
;;;;****************************************************************************

11
(defclass marray (mobject c-array:foreign-array)
12
13
  ((block-pointer :initarg :block-pointer :reader block-pointer)
   (total-size :reader size))
14
15
  (:documentation
   "A superclass for arrays represented in GSL and CL."))
16

Liam Healy's avatar
Liam Healy committed
17
;;; We don't allocate or free the C array data, because that is
18
;;; handled by c-array:foreign-array.  We can use the GSL functions
Liam Healy's avatar
Liam Healy committed
19
20
21
22
23
24
;;; *_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)
25
  (declare (ignore initargs)) ; required by &rest arg for c-array:foreign-array?
Liam Healy's avatar
Liam Healy committed
26
27
28
29
30
31
  ;; 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)
32
	    (size object)
33
34
	    (slot-value object 'block-pointer)
	    blockptr)
Liam Healy's avatar
Liam Healy committed
35
36
      (let ((array-struct (alloc-from-block object blockptr)))
	(setf (slot-value object 'mpointer) array-struct)
37
	#-native (set-struct-array-pointer object)
Liam Healy's avatar
Liam Healy committed
38
39
40
41
	(tg:finalize object
		     (lambda ()
		       (cffi:foreign-free blockptr)
		       (cffi:foreign-free array-struct)))))))
42

43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
(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
66
  (call-next-method))
67

68
(defmethod make-load-form ((object marray) &optional env)
69
  (declare (ignore env))
70
  `(make-marray
71
72
73
74
75
    ',(element-type object)
    :initial-contents
    ',(loop for elt across
	   (subseq (original-array object)
		   (offset object)
76
		   (+ (offset object) (size object)))
77
78
	   collect elt)))

79
80
81
82
83
84
85
86
87
88
89
;;;;****************************************************************************
;;;; 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
90
  (if (member category '(vector matrix))
Liam Healy's avatar
Liam Healy committed
91
      (intern (format nil "~a-~a" category (c-array:cl-single element-type))
92
93
	      :gsl)
      category))
94
95
96
97
98
99
100
101
102
103
104
105

(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)))
106
107
	       ;; Push mapping onto *class-element-type*
	       (pushnew ',(cons class-name element-type-cl)
108
109
			*class-element-type* :test #'equal)
	       (export ',class-name))))
110
	 c-array:*array-element-types*)))
111
112
113
114
115

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

116
117
(export 'make-marray)
(defun make-marray
118
    (class-or-element-type &rest keys &key dimensions initial-contents cl-array
119
     &allow-other-keys)
120
  "Make a GSLL array with the given element type,
121
   :dimensions, and :initial-contents, :initial-element or :cl-array.
122
123
124
125
126
127
128
129
130
131
132
133
134
135
   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))
136

Liam Healy's avatar
Liam Healy committed
137
138
139
140
(pushnew 'marray grid:*grid-data-superclasses*)

(defmethod grid:make-grid-data
    ((type (eql 'marray)) dimensions rest-spec
141
142
     &key initial-element)
  (make-marray (grid:spec-scalar-p rest-spec)
Liam Healy's avatar
Liam Healy committed
143
144
145
	       :dimensions dimensions
	       :initial-element initial-element))

146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
;;; 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))>

161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
(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))))

177
178
179
180
181
182
(defun quote-numeric-list (s)
  (if (and (listp s) (numberp (first s)))
      `',s
      s))

(defvar *row-separator* '^)
183
(export *row-separator*)
184
185
186

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

187
188
189
190
191
192
(set-dispatch-macro-character
 #\# #\m
 (lambda (stream subchar arg)
   (declare (ignore subchar))
   (read-char stream)
   (let ((list (read-delimited-list #\) stream)))
193
194
195
196
197
198
199
200
     `(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)))))))
201

202
;;;;****************************************************************************
203
;;;; Copy to and from bare mpointers 
204
205
;;;;****************************************************************************

206
207
208
209
210
(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."))

211
(defmethod c-array:copy-making-destination ((pointer #.+foreign-pointer-class+))
212
213
214
215
  (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
216
   (c-array:copy-to-destination pointer 'vector-double-float)))
217

Liam Healy's avatar
Liam Healy committed
218
(defmethod c-array:copy-to-destination
219
220
221
222
223
224
225
226
227
    ((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
228
     (c-array:lookup-type class-name *class-element-type*)))))
229
230
231
232
233
234

;; 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.