marray.lisp 8.33 KB
Newer Older
1
;; A "marray" is an array in both GSL and CL
2
;; Liam Healy 2008-04-06 21:23:41EDT
Liam Healy's avatar
Liam Healy committed
3
;; Time-stamp: <2009-12-06 19:50:09EST 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 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
18
19
20
21
22
23
24
25
26
27
28
29
30
31
;;; We don't allocate or free the C array data, because that is
;;; handled by foreign-array.  We can use the GSL functions
;;; *_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)
  (declare (ignore initargs)) ; required by &rest arg for foreign-array?
  ;; 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
111
112
113
114
115
	 *array-element-types*)))

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

137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
;;; 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))>

152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
(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))))

168
169
170
171
172
173
(defun quote-numeric-list (s)
  (if (and (listp s) (numberp (first s)))
      `',s
      s))

(defvar *row-separator* '^)
174
(export *row-separator*)
175
176
177

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

178
179
180
181
182
183
(set-dispatch-macro-character
 #\# #\m
 (lambda (stream subchar arg)
   (declare (ignore subchar))
   (read-char stream)
   (let ((list (read-delimited-list #\) stream)))
184
185
186
187
188
189
190
191
     `(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)))))))
192

193
;;;;****************************************************************************
194
;;;; Copy to and from bare mpointers 
195
196
;;;;****************************************************************************

197
198
199
200
201
202
(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."))

(defmethod copy-making-destination ((pointer #.+foreign-pointer-class+))
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
  (foreign-pointer-method
   pointer
   ;; Default assumption when destination isn't given in #'copy is
   ;; that this should make a vector-double-float.
   (copy-to-destination pointer 'vector-double-float)))

(defmethod copy-to-destination
    ((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
219
     (c-array:lookup-type class-name *class-element-type*)))))
220
221
222
223
224
225

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