1 files changed, 90 insertions, 0 deletions
diff --git a/src/oracles.sc b/src/oracles.sc
new file mode 100644
index 0000000..9680479
--- /dev/null
+++ b/src/oracles.sc
@@ -0,0 +1,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)))