blob: 96804793e9d61247bfb811986528d331a7ce92b0 (
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
|
(sc-include "macros/macros")
(pre-include "math.h")
(pre-include "gsl/gsl_blas.h")
(pre-include "gsl/gsl_randist.h")
(pre-include "utils.h")
(define (bernoulli-extent-generator r p r0 r1) (double (const gsl-rng*) double double double)
(return (if* (gsl-ran-bernoulli r p) r1 r0)))
(define (bernoulli-true-volume p r0 r1 dimension) (double double double double (unsigned int))
(return (* (volume-of-ball dimension)
(+ (* p (gsl-pow-uint r1 dimension))
(* (- 1 p) (gsl-pow-uint r0 dimension))))))
(define (uniform-extent-generator r a b) (double (const gsl-rng*) double double)
(return (gsl-ran-flat r a b)))
;; TODO: Verify the accuracy of this function for non-trivial a, b.
(define (uniform-true-volume a b dimension) (double double double (unsigned int))
(return (- (exp (+ (ln-volume-of-ball dimension)
(* dimension (log b))
(- (log (+ dimension 1)))))
(exp (+ (ln-volume-of-ball dimension)
(* dimension (log a))
(- (log (+ dimension 1))))))))
(define (beta-extent-generator r alpha beta) (double (const gsl-rng*) double double)
(return (gsl-ran-beta r alpha beta)))
(define (beta-true-volume alpha beta dimension) (double double double (unsigned int))
(let* ((vol double (volume-of-ball dimension)))
(for-i r dimension
(set* vol (/ (+ alpha r)
(+ alpha beta r))))
(return vol)))
(define (infinity-norm x) ((static double) (const gsl-vector*))
(let* ((max double (fabs (gsl-vector-get x 0))))
;; TODO: Start this loop from i = 1, not i = 0. That would be
;; slightly faster.
(for-i i (: x size)
(set max (GSL-MAX max (fabs (gsl-vector-get x i)))))
(return max)))
(define (cube-extent-oracle x edge) (double (const gsl-vector*) double)
(return (/ edge 2 (infinity-norm x))))
(sc-define-syntax (compute-cube-extent-oracle-minimizand i)
(/ (- (/ edge 2)
(* (GSL-SIGN (gsl-vector-get x i))
(gsl-vector-get center i)))
(fabs (gsl-vector-get x i))))
(define (cube-extent-oracle-with-center x center edge) (double (const gsl-vector*) (const gsl-vector*) double)
(let* ((min double (compute-cube-extent-oracle-minimizand 0)))
;; TODO: Start this loop from i = 1, not i = 0. That would be
;; slightly faster.
(for-i i (: center size)
(set min (GSL-MIN min (compute-cube-extent-oracle-minimizand i))))
(return min)))
(define (cube-true-volume edge dimension) (double double (unsigned int))
(return (gsl-pow-uint edge dimension)))
(define (cube-maximum-extent edge dimension) (double double (unsigned int))
(return (/ (* edge (sqrt dimension))
2)))
(define (ellipsoid-extent-oracle x axes) (double (const gsl-vector*) (const gsl-vector*))
(let* ((k double 0))
(for-i i (: axes size)
(set+ k (gsl-pow-2 (/ (gsl-vector-get x i)
(gsl-vector-get axes i)))))
(return (/ (sqrt k)))))
(define (ellipsoid-true-volume axes) (double (const gsl-vector*))
(let* ((vol double (volume-of-ball (: axes size))))
(for-i i (: axes size)
(set* vol (gsl-vector-get axes i)))
(return vol)))
(define (spheroid-extent-oracle x eccentricity) (double (const gsl-vector*) double)
(let* ((xsub gsl-vector-const-view
(gsl-vector-const-subvector x 1 (- (: x size) 1))))
(return (/ (sqrt (+ (gsl-pow-2 (gsl-blas-dnrm2 (address-of (struct-get xsub vector))))
(gsl-pow-2 (/ (gsl-vector-get x 0) eccentricity))))))))
(define (spheroid-true-volume eccentricity dimension) (double double (unsigned int))
(return (* (volume-of-ball dimension) eccentricity)))
|
