vector.lisp 17.5 KB
Newer Older
1
2
;; Vectors
;; Liam Healy, Sun Mar 26 2006 - 11:51
3
;; Time-stamp: <2008-02-10 11:06:54EST vector.lisp>
4
;; $Id: $
5
6
7

(in-package :gsl)

8
9
10
11
12
13
14
15
16
17
;;; Vectors are specified in a letm binding with
;;;  (vector-double size-or-initial &optional zero)
;;;  (vector-single size-or-initial &optional zero)
;;;  (vector-fixnum size-or-initial &optional zero)
;;;  (vector-complex size-or-initial &optional zero)
;;; where size-or-initial is a positive integer indicating the
;;; size, and zero indicates that all elements should be set to zero,
;;; or, size-or-initial is a sequence to which the vector
;;; should be initially set.

18
19
20
;;; No mechanism for C stream input/output yet.
;;; Generalize check-gsl-status to optionally signal errors, use here?
;;; Functions like write-binary etc. as a single function, selecting the C fn with typecase?
21
22
;;; #'subvector, #'subvector-stride cause crash, see notes 2006-03-30
;;; #'vector-complex-real, #'vector-complex-imag need structure definition
23

24
;;; GSL bug?:  no gsl_vector_complex_add, etc.
25

26
;;; Need to build real vector out of view pointer.
27

28
29
30
;;; Need to #'cl-invalidate when (setf gsl-aref) called, see Mon Nov 26 2007.


31
;;;;****************************************************************************
32
;;;; Vector structure and CL object
33
34
35
;;;;****************************************************************************

#|
36
37
38
;;FDL
The size is simply the number of vector elements.  The range of
valid indices runs from 0 to size-1.  The stride is the
39
step-size from one element to the next in physical memory, measured in
40
units of the appropriate datatype.  The pointer data gives the
41
location of the first element of the vector in memory.  The pointer
42
block stores the location of the memory block in which the vector
43
elements are located (if any).  If the vector owns this block then the
44
owner field is set to one and the block will be deallocated when the
45
vector is freed.  If the vector points to a block owned by another
46
object then the owner field is zero and any underlying block will not be
47
48
49
50
deallocated with the vector.
|#

;;; GSL-vector definition
liam's avatar
liam committed
51
(cffi:defcstruct gsl-vector-c
52
53
54
55
56
57
  (size :size)
  (stride :size)
  (data :pointer)
  (block :pointer)
  (owner :int))

58
59
(defclass gsl-vector (gsl-data) ())

60
;;; Allocation, freeing, reading and writing
61
62
63
64
(defdata "vector" vector-double double-float gsl-vector)
(defdata "vector_float" vector-single single-float gsl-vector)
(defdata "vector_int" vector-fixnum fixnum gsl-vector)
(defdata "vector_complex" vector-complex complex gsl-vector)
65

66
67
;;; (defdata "vector_uint" vector-unsigned-fixnum (unsigned-byte 64) gsl-vector)
;;; doesn't work 2007-03-04.
68
69
70
71
;;; (deftype unsigned-fixnum () `(integer 0 ,most-positive-fixnum))
;;; (defdata "vector_uint" vector-unsigned-fixnum unsigned-fixnum gsl-vector)
;;; doesn't work either 2008-02-10.

72
73
74
75

(defmacro defun-gsl-vdsfc (&rest args)
  "A defun-gsl for vectors of double, single, fixnum, and complex."
  (defun-gsl-all
76
      '(double single fixnum complex)
77
78
79
80
81
82
83
84
      '(:double :float :int gsl-complex)
    "vector"
    'gsl-vector
    args))

(defmacro defun-gsl-vdsf (&rest args)
  "A defun-gsl for vectors of double, single, and fixnum."
  (defun-gsl-all
85
      '(double single fixnum)
86
87
88
89
      '(:double :float :int)
    "vector"
    'gsl-vector
    args))
90

91
(defmethod gsl-array ((object gsl-vector))
92
  (cffi:foreign-slot-value (pointer object) 'gsl-vector-c 'data))
93

94
95
96
97
98
99
100
101
102
103
(defun make-data-from-pointer (pointer &optional (class 'gsl-vector-double) size)
  "Given a C pointer to a GSL data type, make the CL object."
  (make-instance
   class
   :pointer pointer
   :storage-size
   (or size (cffi:foreign-slot-value pointer 'gsl-vector-c 'size))))

(export 'vector-data)
(defun vector-data (pointer)
104
105
106
107
  "A pointer to the GSL array with the data contents, from the
   sruct pointer."
  (cffi:foreign-slot-value pointer 'gsl-vector-c 'data))

108
;;;;****************************************************************************
109
;;;; Getting values
110
111
;;;;****************************************************************************

112
(defun-gsl-vdsfc gsl-aref ((vector gsl-vector) &rest indices)
113
114
    "gsl_vector_get"
  (((pointer vector) :pointer) ((first indices) :size))
115
  :c-return :c-base-type
116
117
  :documentation			; FDL
  "The ith element of the vector.")
118

119
120
121
122
123
124
125
(defun-gsl vref (pointer index)
  "gsl_vector_get"
  ((pointer :pointer) (index :size))
  :c-return :double
  :index nil
  :documentation "An element of the vector of doubles, computed from the pointer.")

126
127
128
(defun-gsl gsl-vector-ptr (vector i)
  "gsl_vector_ptr" (((pointer vector) :pointer) (i :size))
  :c-return :pointer
129
130
  :documentation			; FDL
  "The ith element of the vector as a pointer.")
131

132
;;;;****************************************************************************
133
;;;; Setting values
134
135
;;;;****************************************************************************

136
(defun-gsl-vdsfc (setf gsl-aref) (value (vector gsl-vector) &rest indices)
137
  "gsl_vector_set"
138
  (((pointer vector) :pointer) ((first indices) :size) (value :c-base-type))
139
  :c-return :void
140
141
  :documentation			; FDL
  "Set an element of the vector.")
142
143
144
145
146
147

(defun-gsl (setf vref) (value pointer index)
  "gsl_vector_set"
  ((pointer :pointer) (index :size) (value :double))
  :c-return :void
  :index nil
148
149
  :documentation			; FDL
 "Set an element of the vector of doubles, using its pointer.")
150

151
(defun-gsl-vdsfc set-all ((object gsl-vector) value)
152
  "gsl_vector_set_all"
153
  (((pointer object) :pointer) (value :c-base-type))
154
  :c-return :void)
155

156
(defun-gsl-vdsfc set-zero ((object gsl-vector))
157
158
  "gsl_vector_set_zero" (((pointer object) :pointer))
  :c-return :void)
159

160
(defun-gsl-vdsfc set-basis ((vector gsl-vector) index)
161
  "gsl_vector_set_basis" (((pointer vector) gsl-vector-c) (index :size))
162
163
164
  :invalidate (vector)
  :documentation			; FDL
  "Set the index element to 1, and the rest to 0.")
165

166
167
168
169
170
;;;;****************************************************************************
;;;; Views
;;;;****************************************************************************

(cffi:defcstruct gsl-vector-view
171
  (vector gsl-vector-c))
172

173
174
175
;;; broken
;;; (letm ((vec (vector-double  #(-3.21d0 1.0d0 12.8d0)))) (subvector vec 1 2))

176
177
178
(export '(subvector subvector-stride))

(defgeneric subvector (vector offset size)
179
  (:documentation			; FDL
180
   "Return a vector view of a subvector of another vector
181
182
  v.  The start of the new vector is offset by offset elements
  from the start of the original vector.  The new vector has size
183
  elements."))
184

185
(defgeneric subvector-stride (vector offset stride size)
186
187
188
189
190
191
192
  (:documentation			; FDL
   "A vector view of a subvector of another vector
  v with an additional stride argument. The subvector is formed in
  the same way as for #'subvector but the new vector has
  n elements with a step-size of stride from one element to
  the next in the original vector.  Mathematically, the i-th element
  of the new vector v' is given by
193
  v'(i) = v->data[(offset + i*stride)*v->stride]
194
  where the index i runs from 0 to n-1.
195
  Note that subvector views give direct access to the underlying elements
196
197
198
199
200
201
202
203
204
205
206
207
  of the original vector."))

(defun-gsl-vdsfc subvector ((vector gsl-vector) offset size)
    "gsl_vector_subvector"
  (((pointer vector) gsl-vector-c) (offset :size) (size :size))
  :c-return :pointer)

(defun-gsl-vdsfc subvector-stride ((vector gsl-vector) offset stride size)
  "gsl_vector_subvector_with_stride"
  (((pointer vector) gsl-vector-c)
   (offset :size) (stride :size) (size :size))
  :c-return :pointer)
208
209
210

;;; These require that the gsl-vector-complex structure be defined.
#|
211
212
213
(defun-gsl vector-complex-real (vector)
  "gsl_vector_complex_real" ((vector gsl-vector-complex))
  :c-return :pointer
214
  :documentation
215
  "A vector view of the real parts of the complex vector v.")
216

217
218
219
(defun-gsl vector-complex-imag (vector)
  "gsl_vector_complex_imag"((vector gsl-vector-complex))
  :c-return :pointer
220
  :documentation
221
  "A vector view of the imaginary parts of the complex vector v.")
222
223
|#

224
225
226
(defun-gsl vector-array (base size)
  "gsl_vector_view_array" ((base :pointer) (size :size))
  :c-return :pointer
227
  :documentation			; FDL
228
  "A vector view of an array.  The start of the new
229
  vector is given by base and has n elements.")
230

231
(defun-gsl vector-array-stride (base stride size)
232
  "gsl_vector_view_array_with_stride"
233
234
  ((base :pointer) (stride :size) (size :size))
  :c-return :pointer
235
  :documentation			; FDL
236
  "A vector view of an array with stride.  The start of the new
237
  vector is given by base.")
238

239
240
241
242
;;;;****************************************************************************
;;;; Copying
;;;;****************************************************************************

243
(defun-gsl-vdsfc copy ((destination gsl-vector) (source gsl-vector))
244
  "gsl_vector_memcpy"
245
246
  (((pointer destination) gsl-vector-c) ((pointer source) gsl-vector-c))
  :invalidate (destination)
247
248
249
  :documentation			; FDL
  "Copy the elements of the vector source into the
   vector destination.  The two vectors must have the same length.")
250

251
(defun-gsl-vdsfc swap ((v gsl-vector) (w gsl-vector))
252
253
  "gsl_vector_swap" (((pointer v) gsl-vector-c) ((pointer w) gsl-vector-c))
  :invalidate (v w)
254
255
  :documentation			; FDL
  "Exchange the elements of the vectors v and w
256
257
   by copying.  The two vectors must have the same length.")

258
259
260
261
;;;;****************************************************************************
;;;; Exchanging elements
;;;;****************************************************************************

262
(export '(swap-elements vector-reverse))
263

264
(defgeneric swap-elements (vec i j)
265
266
  (:documentation			; FDL
   "Exchange the i-th and j-th elements of the vector vec in-place."))
267
268
269

(defgeneric vector-reverse (vec)
  (:documentation
270
   "Exchange the i-th and j-th elements of the vector vec in-place."))
271

272
(defun-gsl-vdsfc swap-elements ((vec gsl-vector) i j)
273
  "gsl_vector_swap_elements" (((pointer vec) gsl-vector-c) (i :size) (j :size))
274
275
276
277
  :after ((when (listp (cl-invalid vec))
	    (push (list i) (cl-invalid vec))
	    (push (list j) (cl-invalid vec))))
  :return (vec))
278

279
(defun-gsl-vdsfc vector-reverse ((vec gsl-vector))
280
281
  "gsl_vector_reverse" (((pointer vec) gsl-vector-c))
  :invalidate (vec)
282
283
  :documentation			; FDL
  "Reverse the order of the elements of the vector vec.")
284

285
286
287
288
;;;;****************************************************************************
;;;; Arithmetic operations
;;;;****************************************************************************

289
(defun-gsl-vdsf gsl+ ((a gsl-vector) (b gsl-vector))
290
291
  "gsl_vector_add" (((pointer a) gsl-vector-c) ((pointer b) gsl-vector-c))
  :invalidate (a)
292
293
294
  :documentation			; FDL
  "Add the elements of vector b to the elements of
  vector a, a'_i = a_i + b_i. The two vectors must have the
295
  same length.")
296

297
(defun-gsl-vdsf gsl- ((a gsl-vector) (b gsl-vector))
298
299
  "gsl_vector_sub" (((pointer a) gsl-vector-c) ((pointer b) gsl-vector-c))
  :invalidate (a)
300
301
302
  :documentation			; FDL
  "Subtract the elements of vector b from the elements of
   vector a, a'_i = a_i - b_i.  The two vectors must have the same length.")
303

304
(defun-gsl-vdsf gsl* ((a gsl-vector) (b gsl-vector))
305
306
  "gsl_vector_mul" (((pointer a) gsl-vector-c) ((pointer b) gsl-vector-c))
  :invalidate (a)
307
308
309
  :documentation			; FDL
  "Multiply the elements of vector a by the elements of
  vector b, a'_i = a_i * b_i. The two vectors must have the same length.")
310

311
(defun-gsl-vdsf gsl/ ((a gsl-vector) (b gsl-vector))
312
313
  "gsl_vector_div" (((pointer a) gsl-vector-c) ((pointer b) gsl-vector-c))
  :invalidate (a)
314
315
316
  :documentation			; FDL
  "Divide the elements of vector a by the elements of
  vector b, a'_i = a_i / b_i. The two vectors must have the same length.")
317

318
(defun-gsl-vdsf gsl*c ((a gsl-vector) x)
319
320
  "gsl_vector_scale" (((pointer a) gsl-vector-c) (x :double))
  :invalidate (a)
321
322
  :documentation			; FDL
  "Multiply the elements of vector a by the constant factor x, a'_i = x a_i.")
323

324
(defun-gsl-vdsf gsl+c ((a gsl-vector) x)
325
326
  "gsl_vector_add_constant" (((pointer a) gsl-vector-c) (x :double))
  :invalidate (a)
327
328
  :documentation			; FDL
  "Add the constant value x to the elements of the vector a, a'_i = a_i + x.")
329

330
331
332
333
;;;;****************************************************************************
;;;; Maximum and minimum elements
;;;;****************************************************************************

334
(defun-gsl-vdsf gsl-max ((v gsl-vector))
335
  "gsl_vector_max" (((pointer v) gsl-vector-c))
336
337
338
  :c-return :c-base-type
  :documentation			; FDL
  "The maximum value in the vector v.")
liam's avatar
liam committed
339

340
(defun-gsl-vdsf gsl-min ((v gsl-vector))
341
  "gsl_vector_min" (((pointer v) gsl-vector-c))
342
343
344
  :c-return :c-base-type
  :documentation			; FDL
  "The minimum value in the vector v.")
liam's avatar
liam committed
345

346
347
348
(defun-gsl-vdsf gsl-minmax ((v gsl-vector))
  "gsl_vector_minmax"
  (((pointer v) gsl-vector-c) (min :c-base-type) (max :c-base-type))
349
350
351
  :c-return :void
  :documentation			; FDL
  "The minimum and maximum values in the vector v.")
liam's avatar
liam committed
352

353
(defun-gsl-vdsf gsl-max-index ((v gsl-vector))
354
  "gsl_vector_max_index" (((pointer v) gsl-vector-c))
355
356
357
  :c-return :size
  :documentation			; FDL
  "The index of the maximum value in the vector v.
liam's avatar
liam committed
358
   When there are several equal minimum elements then the lowest index is
359
   returned.")
liam's avatar
liam committed
360

361
(defun-gsl-vdsf gsl-min-index ((v gsl-vector))
362
  "gsl_vector_min_index" (((pointer v) gsl-vector-c))
363
364
365
366
  :c-return :size
  :documentation			; FDL
  "The index of the minimum value in the vector v.  When there are several
  equal minimum elements then the lowest index is returned.")
liam's avatar
liam committed
367

368
(defun-gsl-vdsf gsl-minmax-index ((v gsl-vector))
liam's avatar
liam committed
369
  "gsl_vector_minmax_index"
370
  (((pointer v) gsl-vector-c) (imin :size) (imax :size))
371
372
373
  :c-return :void
  :documentation			; FDL
  "The indices of the minimum and maximum values in the vector v.
liam's avatar
liam committed
374
  When there are several equal minimum elements then the lowest index is
375
  returned.")
liam's avatar
liam committed
376

377
378
379
380
;;;;****************************************************************************
;;;; Properties
;;;;****************************************************************************

381
(defun-gsl-vdsfc gsl-zerop ((v gsl-vector))
382
  "gsl_vector_isnull" (((pointer v) gsl-vector-c))
383
384
385
  :c-return :boolean
  :documentation			; FDL
  "All elements of vector v are zero.")
liam's avatar
liam committed
386

387
;;;;****************************************************************************
388
;;;; Examples and unit tests
389
390
391
392
393
;;;;****************************************************************************

(lisp-unit:define-test vector-fixnum
  (lisp-unit:assert-eql			;(setf gsl-aref), gsl-aref
   77
394
395
396
   (letm ((intvec (vector-fixnum 4)))
     (setf (gsl-aref intvec 1) 77)
     (gsl-aref intvec 1)))
397
398
  (lisp-unit:assert-equalp		;(setf data)
   #(4 6 8 2)
399
400
401
   (letm ((intvec (vector-fixnum 4)))
     (setf (data intvec) #(4 6 8 2))
     (data intvec)))
402
403
  (lisp-unit:assert-equalp		;set-zero
   #(0 0 0 0)
404
405
406
   (letm ((intvec (vector-fixnum 4)))
     (set-zero intvec)
     (data intvec)))
407
408
  (lisp-unit:assert-equalp		;set-all
   #(44 44 44 44)
409
410
411
   (letm ((intvec (vector-fixnum 4)))
     (set-all intvec 44)
     (data intvec)))
412
413
  (lisp-unit:assert-equalp		;set-basis
   #(0 1 0 0)
414
415
416
   (letm ((intvec (vector-fixnum 4)))
     (set-basis intvec 1)
     (data intvec)))
417
418
  (lisp-unit:assert-equalp		;vector-reverse
   #(4 3 2 1)
419
420
421
   (letm ((intvec (vector-fixnum #(1 2 3 4))))
     (vector-reverse intvec)
     (data intvec)))
422
423
  (lisp-unit:assert-eql			;gsl-min
   -12
424
425
   (letm ((intvec (vector-fixnum #(-1 -12 8 3))))
     (gsl-min intvec)))
426
427
  (lisp-unit:assert-eql			;gsl-max
   8
428
429
   (letm ((intvec (vector-fixnum #(-1 -12 8 3))))
     (gsl-max intvec)))
430
431
  (lisp-unit:assert-equal		;gsl-minmax
   '(-12 8)
432
433
   (letm ((intvec (vector-fixnum #(-1 -12 8 3))))
     (multiple-value-list (gsl-minmax intvec))))
434
435
  (lisp-unit:assert-eql			;gsl-min-index
   1
436
437
   (letm ((intvec (vector-fixnum #(-1 -12 8 3))))
     (gsl-min-index intvec)))
438
439
  (lisp-unit:assert-eql			;gsl-max-index
   2
440
441
   (letm ((intvec (vector-fixnum #(-1 -12 8 3))))
     (gsl-max-index intvec)))
442
443
  (lisp-unit:assert-equal		;gsl-minmax-index
   '(1 2)
444
445
   (letm ((intvec (vector-fixnum #(-1 -12 8 3))))
     (multiple-value-list (gsl-minmax-index intvec))))
446
  (lisp-unit:assert-equalp		;copy
447
   #(1 2 3 4)
448
449
450
451
   (letm ((intvec1 (vector-fixnum #(1 2 3 4)))
	  (intvec2 (vector-fixnum 4)))
     (copy intvec2 intvec1)
     (data intvec2)))
452
  (lisp-unit:assert-equalp		;swap
453
   #(5 6 7 8 1 2 3 4)
454
455
456
457
   (letm ((intvec1 (vector-fixnum #(1 2 3 4)))
	  (intvec2 (vector-fixnum #(5 6 7 8))))
     (swap intvec2 intvec1)
     (concatenate 'vector (data intvec1) (data intvec2))))
458
  (lisp-unit:assert-equalp		;swap-elements
459
   #(1 4 3 2)
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
   (letm ((intvec (vector-fixnum #(1 2 3 4))))
     (swap-elements intvec 1 3)
     (data intvec))))

(lisp-unit:define-test vector-double
  (lisp-unit:assert-equal
   '("-0.321000000000d+01" "0.100000000000d+01" "0.128000000000d+02")
   (lisp-unit:fp-sequence 
    (letm ((vec (vector-double 3)))
      (setf (gsl-aref vec 0) -3.21d0
	    (gsl-aref vec 1) 1.0d0
	    (gsl-aref vec 2) 12.8d0
	    (cl-invalid vec) t)
      (data vec)))
   (lisp-unit:assert-equal
    '("-0.321000000000d+01" "0.100000000000d+01" "0.128000000000d+02")
    (lisp-unit:fp-sequence 
     (letm ((vec (vector-double 3)))
       (setf (data vec) #(-3.21d0 1.0d0 12.8d0))
       (data vec))))
   (lisp-unit:assert-equal
    '("-0.321000000000d+01" "0.100000000000d+01" "0.128000000000d+02")
    (lisp-unit:fp-sequence 
     (letm ((vec (vector-double #(-3.21d0 1.0d0 12.8d0))))
       (data vec))))
   (lisp-unit:assert-equal
    '("0.000000000000d+01" "0.100000000000d+01" "0.000000000000d+01"
      "0.000000000000d+01" "0.000000000000d+01")
    (lisp-unit:fp-sequence 
     (letm ((base (vector-double 5)))
       (set-basis base 1)
       (data base))))
   (lisp-unit:assert-equal
    '("-0.321000000000d+01" "0.100000000000d+01" "0.128000000000d+02")
    (lisp-unit:fp-sequence 
     (letm ((vec1 (vector-double #(-3.21d0 1.0d0 12.8d0)))
	    (vec2 (vector-double 3)))
       (copy vec2 vec1)
       (data vec2))))))