blob: 180ac91d72040df478895ecbd750674b61c5002f (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
|
(define-module (extent-sampling wrap)
#:use-module (srfi srfi-1)
#:use-module (srfi srfi-26)
#:use-module (ice-9 match)
#:use-module (system foreign))
(include "load-libs.scm")
(define (make-gsl-alloc allocator-name allocator-args freer-name)
(lambda args
(let ((obj
(apply (pointer->procedure
'* (dynamic-func allocator-name lib-gsl) allocator-args)
args)))
(set-pointer-finalizer! obj (dynamic-func freer-name lib-gsl))
obj)))
(define-public pi
((pointer->procedure double
(dynamic-func "pi" lib-extentsampling)
(list))))
(define-public volume-of-ball
(pointer->procedure double
(dynamic-func "volume_of_ball" lib-extentsampling)
(list int)))
(define-public surface-area-of-ball
(pointer->procedure double
(dynamic-func "surface_area_of_ball" lib-extentsampling)
(list int)))
(define-public angle-between-vectors
(pointer->procedure double
(dynamic-func "angle_between_vectors" lib-extentsampling)
(list '* '*)))
(define-public planar-angle->solid-angle
(pointer->procedure double
(dynamic-func "planar_angle_to_solid_angle" lib-extentsampling)
(list double int)))
;; Random state
(define %gsl-random-state
((make-gsl-alloc "gsl_rng_alloc" (list '*) "gsl_rng_free")
(dereference-pointer
(dynamic-pointer "gsl_rng_default" lib-gsl))))
(define-public (set-gsl-random-state! seed)
((pointer->procedure void
(dynamic-func "gsl_rng_set" lib-gsl)
(list '* unsigned-long))
%gsl-random-state seed))
(define-public random-flat
(cut (pointer->procedure double
(dynamic-func "gsl_ran_flat" lib-gsl)
(list '* double double))
%gsl-random-state <> <>))
;; Vector functions
(define vector-alloc
(make-gsl-alloc "gsl_vector_alloc" (list int) "gsl_vector_free"))
(define-public (vector-size vector)
(match (parse-c-struct vector
(list size_t size_t '* '* int))
((size _ _ _ _) size)))
(define-public (vector-scale! vector k)
"Scale VECTOR by a factor K."
((pointer->procedure int
(dynamic-func "gsl_vector_scale" lib-gsl)
(list '* double))
vector k)
vector)
(define-public (vector-norm2 vector)
((pointer->procedure double
(dynamic-func "gsl_blas_dnrm2" lib-gsl)
(list '*))
vector))
(define-public (vector-add! u v)
"Return the vector sum of vectors U and V. The sum is stored in U."
((pointer->procedure int
(dynamic-func "gsl_vector_add" lib-gsl)
(list '* '*))
u v)
u)
(define-public (basis n dimension)
"Return the basis vector in the Nth canonical direction."
(let ((vector (vector-alloc dimension)))
((pointer->procedure void
(dynamic-func "gsl_vector_set_basis" lib-gsl)
(list '* int))
vector n)
vector))
;; Histogram functions
(define histogram-struct
(list size_t '* '*))
(define (histogram-alloc bins min max)
(let ((histogram
((make-gsl-alloc "gsl_histogram_alloc" (list size_t) "gsl_histogram_free")
bins)))
((pointer->procedure int
(dynamic-func "gsl_histogram_set_ranges_uniform" lib-gsl)
(list '* double double))
histogram min max)
histogram))
(define-public (histogram->points histogram)
(match (parse-c-struct histogram histogram-struct)
((n range bin)
(let ((range (array->list (pointer->bytevector range (1+ n) 0 'f64))))
(map cons
(drop-right range 1)
(array->list (pointer->bytevector bin n 0 'f64)))))))
(define-public (histogram->pdf histogram)
(match (parse-c-struct histogram histogram-struct)
((n range bin)
(let ((range (pointer->bytevector range (1+ n) 0 'f64)))
((pointer->procedure int
(dynamic-func "gsl_histogram_scale" lib-extentsampling)
(list '* double))
histogram
(/ n
(- (array-ref range n)
(array-ref range 0))
(fold + 0 (array->list (pointer->bytevector bin n 0 'f64)))))
histogram))))
(define-public (rhistogram bins min max)
(case-lambda
(() (histogram-alloc bins min max))
((histogram) histogram)
((histogram input)
((pointer->procedure int
(dynamic-func "gsl_histogram_increment" lib-extentsampling)
(list '* double))
histogram input)
histogram)))
;; nd-random
(define-public (random-direction-vector dimension)
(let ((vector (vector-alloc dimension)))
((pointer->procedure void
(dynamic-func "random_direction_vector" lib-extentsampling)
(list '* '*))
%gsl-random-state vector)
vector))
(define-public (subsampling-random-vector mean max-theta)
(let ((vector (vector-alloc (vector-size mean))))
((pointer->procedure void
(dynamic-func "subsampling_random_vector" lib-extentsampling)
(list '* '* double '*))
%gsl-random-state mean max-theta vector)
vector))
(define-public (shifted-gaussian-random-vector mean max-theta standard-deviation)
(let* ((vector (vector-alloc (vector-size mean)))
(cost
((pointer->procedure int
(dynamic-func "shifted_gaussian_random_vector" lib-extentsampling)
(list '* '* double double '*))
%gsl-random-state mean max-theta standard-deviation vector)))
vector))
(define-public (shifted-gaussian-random-vector-cost mean max-theta standard-deviation)
((pointer->procedure int
(dynamic-func "shifted_gaussian_random_vector" lib-extentsampling)
(list '* '* double double '*))
%gsl-random-state mean max-theta standard-deviation (vector-alloc (vector-size mean))))
(define %integration-workspace
((make-gsl-alloc "gsl_integration_workspace_alloc" (list size_t) "gsl_integration_workspace_free")
1000))
(define-public (shifted-gaussian-pdf theta mean max-theta standard-deviation)
((pointer->procedure double
(dynamic-func "shifted_gaussian_pdf" lib-extentsampling)
(list double double double double unsigned-int '*))
theta (vector-norm2 mean) max-theta standard-deviation (vector-size mean) %integration-workspace))
|
