aboutsummaryrefslogtreecommitdiff
path: root/src/oracles.sc
blob: f0f53a8594b38f01a336111635bf4fc7c9430326 (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
(sc-include "macros/macros")

(pre-include "math.h")
(pre-include "gsl/gsl_blas.h")
(pre-include "gsl/gsl_randist.h")
(pre-include "oracles.h")
(pre-include "utils.h")

(define (bernoulli-extent-oracle r x params) (double (const gsl-rng*) (const gsl-vector*) void*)
  (let* ((params (const bernoulli-params*) (convert-type params bernoulli-params*)))
    (return (if* (gsl-ran-bernoulli r (: params p))
                 (: params r1)
                 (: params r0)))))

(define (bernoulli-true-volume dimension params) (double (unsigned int) void*)
  (let* ((params (const bernoulli-params*) (convert-type params bernoulli-params*)))
    (return (* (volume-of-ball dimension)
               (+ (* (: params p) (gsl-pow-uint (: params r1) dimension))
                  (* (- 1 (: params p)) (gsl-pow-uint (: params r0) dimension)))))))

(define (uniform-extent-oracle r x params) (double (const gsl-rng*) (const gsl-vector*) void*)
  (let* ((params (const uniform-params*) (convert-type params uniform-params*)))
    (return (gsl-ran-flat r (: params a) (: params b)))))

;; TODO: Verify the accuracy of this function for non-trivial a, b.
(define (uniform-true-volume dimension params) (double (unsigned int) void*)
  (let* ((params (const uniform-params*) (convert-type params uniform-params*)))
    (return (- (exp (+ (ln-volume-of-ball dimension)
                       (* dimension (log (: params b)))
                       (- (log (+ dimension 1)))))
               (exp (+ (ln-volume-of-ball dimension)
                       (* dimension (log (: params a)))
                       (- (log (+ dimension 1)))))))))

(define (beta-extent-oracle r x params) (double (const gsl-rng*) (const gsl-vector*) void*)
  (let* ((params (const beta-params*) (convert-type params beta-params*)))
    (return (gsl-ran-beta r (: params alpha) (: params beta)))))

(define (beta-true-volume dimension params) (double (unsigned int) void*)
  (let* ((params (const beta-params*) (convert-type params beta-params*))
         (vol double (volume-of-ball dimension)))
    (for-i r dimension
      (set* vol (/ (+ (: params alpha) r)
                   (+ (: params alpha) (: params 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 r x params) (double (const gsl-rng*) (const gsl-vector*) void*)
  (let* ((params (const cube-params*) (convert-type params cube-params*)))
    (return (/ (: params edge) 2 (infinity-norm x)))))

(sc-define-syntax (compute-cube-extent-oracle-minimizand i)
  (/ (- (/ (: params edge) 2)
        (* (GSL-SIGN (gsl-vector-get x i))
           (gsl-vector-get (: params center) i)))
     (fabs (gsl-vector-get x i))))

(define (cube-extent-oracle-with-center r x params) (double (const gsl-rng*) (const gsl-vector*) void*)
  (let* ((params (const cube-params*) (convert-type params cube-params*))
         (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 (: (: params center) size)
      (set min (GSL-MIN min (compute-cube-extent-oracle-minimizand i))))
    (return min)))

(define (cube-true-volume dimension params) (double (unsigned int) void*)
  (let* ((params (const cube-params*) (convert-type params cube-params*)))
    (return (gsl-pow-uint (: params edge) dimension))))

(define (ellipsoid-extent-oracle r x params) (double (const gsl-rng*) (const gsl-vector*) void*)
  (let* ((params (const ellipsoid-params*) (convert-type params ellipsoid-params*))
         (k double 0))
    (for-i i (: (: params axes) size)
      (set+ k (gsl-pow-2 (/ (gsl-vector-get x i)
                            (gsl-vector-get (: params axes) i)))))
    (return (/ (sqrt k)))))

(define (ellipsoid-true-volume dimension params) (double (unsigned int) void*)
  (let* ((params (const ellipsoid-params*) (convert-type params ellipsoid-params*))
         (vol double (volume-of-ball (: (: params axes) size))))
    (for-i i (: (: params axes) size)
      (set* vol (gsl-vector-get (: params axes) i)))
    (return vol)))

(define (spheroid-extent-oracle r x params) (double (const gsl-rng*) (const gsl-vector*) void*)
  (let* ((params (const spheroid-params*) (convert-type params spheroid-params*))
         (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)
                                      (: params eccentricity)))))))))

(define (spheroid-true-volume dimension params) (double (unsigned int) void*)
  (let* ((params (const spheroid-params*) (convert-type params spheroid-params*)))
    (return (* (volume-of-ball dimension)
               (: params eccentricity)))))