Change extent oracle and true volume prototypes.

This change will make it easier for the scheme wrapper to curry extent
oracles.

* include/extent-sampling.h (extent_oracle_t): Define as a struct that
bundles the parameters required by the actual extent oracle.
* include/oracles.h (bernoulli_params, uniform_params, beta_params,
cube_params, ellipsoid_params, spheroid_params): New types.
(bernoulli_extent_generator, uniform_extent_generator,
beta_extent_generator): Rename to bernoulli_extent_oracle,
uniform_extent_oracle, beta_extent_oracle respectively and change to
the new extent oracle prototype.
(bernoulli_true_volume, uniform_true_volume, beta_true_volume,
cube_extent_oracle, cube_extent_oracle_with_center, cube_true_volume,
ellipsoid_extent_oracle, ellipsoid_true_volume,
spheroid_extent_oracle, spheroid_true_volume): Change to new extent
oracle and true volume prototypes.
* src/oracles.sc: Likewise.
* src/extent-sampling.sc (invoke-extent-oracle): New macro.
(volume, volume-window, integral, volume-cone): Call using new extent
oracle prototype.

4 files changed, 113 insertions, 64 deletions

diff --git a/include/extent-sampling.h b/include/extent-sampling.h
index 8ed3fb2..fd11bc4 100644
--- a/include/extent-sampling.h
+++ b/include/extent-sampling.h
@@ -5,7 +5,11 @@
 #include <gsl/gsl_rstat.h>
 #include <gsl/gsl_vector.h>
 
-typedef double (*extent_oracle_t) (const gsl_vector*);
+typedef struct {
+  double (*oracle) (const gsl_rng*, const gsl_vector*, void*);
+  void *params;
+} extent_oracle_t;
+
 typedef double (*integrand_t) (double, const gsl_vector*);
 
 void init_random (void);
diff --git a/include/oracles.h b/include/oracles.h
index b8421ed..c294b84 100644
--- a/include/oracles.h
+++ b/include/oracles.h
@@ -4,23 +4,48 @@
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_vector.h>
 
-double bernoulli_extent_generator (const gsl_rng* r, double p, double r0, double r1);
-double bernoulli_true_volume (double p, double r0, double r1, unsigned int dimension);
+typedef struct {
+  double p, r0, r1;
+} bernoulli_params;
 
-double uniform_extent_generator (const gsl_rng* r, double a, double b);
-double uniform_true_volume (double a, double b, unsigned int dimension);
+double bernoulli_extent_oracle (const gsl_rng *r, const gsl_vector *x, void *params);
+double bernoulli_true_volume (unsigned int dimension, void *params);
 
-double beta_extent_generator (const gsl_rng* r, double alpha, double beta);
-double beta_true_volume (double alpha, double beta, unsigned int dimension);
+typedef struct {
+  double a, b;
+} uniform_params;
 
-double cube_extent_oracle (const gsl_vector* x, double edge);
-double cube_extent_oracle_with_center (const gsl_vector* x, const gsl_vector* center, double edge);
-double cube_true_volume (double edge, unsigned int dimension);
+double uniform_extent_oracle (const gsl_rng *r, const gsl_vector *x, void *params);
+double uniform_true_volume (unsigned int dimension, void *params);
 
-double ellipsoid_extent_oracle (const gsl_vector* x, const gsl_vector* axes);
-double ellipsoid_true_volume (const gsl_vector* axes);
+typedef struct {
+  double alpha, beta;
+} beta_params;
 
-double spheroid_extent_oracle (const gsl_vector* x, double eccentricity);
-double spheroid_true_volume (double eccentricity, unsigned int dimension);
+double beta_extent_oracle (const gsl_rng *r, const gsl_vector *x, void *params);
+double beta_true_volume (unsigned int dimension, void *params);
+
+typedef struct {
+  double edge;
+  const gsl_vector* center;
+} cube_params;
+
+double cube_extent_oracle (const gsl_rng *r, const gsl_vector *x, void *params);
+double cube_extent_oracle_with_center (const gsl_rng *r, const gsl_vector *x, void *params);
+double cube_true_volume (unsigned int dimension, void *params);
+
+typedef struct {
+  const gsl_vector* axes;
+} ellipsoid_params;
+
+double ellipsoid_extent_oracle (const gsl_rng *r, const gsl_vector *x, void *params);
+double ellipsoid_true_volume (unsigned int dimension, void *params);
+
+typedef struct {
+  double eccentricity;
+} spheroid_params;
+
+double spheroid_extent_oracle (const gsl_rng *r, const gsl_vector *x, void *params);
+double spheroid_true_volume (unsigned int dimension, void *params);
 
 #endif
diff --git a/src/extent-sampling.sc b/src/extent-sampling.sc
index c6220e6..b54c44b 100644
--- a/src/extent-sampling.sc
+++ b/src/extent-sampling.sc
@@ -25,6 +25,9 @@
               gsl-integration-workspace-free
               body ...))
 
+(sc-define-syntax (invoke-extent-oracle extent-oracle r x)
+  ((struct-get extent-oracle oracle) r x (struct-get extent-oracle params)))
+
 (pre-define CONFIDENCE-INTERVAL-FACTOR 1.96)
 
 (pre-let* (VOLUME-MINIMUM-SAMPLES 100)
@@ -40,7 +43,7 @@
                       (> (rerror volume true-volume) rtol)
                       (< (gsl-rstat-n stats) VOLUME-MINIMUM-SAMPLES))
           (random-direction-vector r x)
-          (gsl-rstat-add (exp (+ vn (* dimension (log (extent-oracle x)))))
+          (gsl-rstat-add (exp (+ vn (* dimension (log (invoke-extent-oracle extent-oracle r x)))))
                          stats)
           (set volume (gsl-rstat-mean stats)))))
     (return volume)))
@@ -57,7 +60,7 @@
       (with-vector x dimension
         (do-while (< accurate-estimates window-length)
           (random-direction-vector r x)
-          (gsl-rstat-add (exp (+ vn (* dimension (log (extent-oracle x)))))
+          (gsl-rstat-add (exp (+ vn (* dimension (log (invoke-extent-oracle extent-oracle r x)))))
                          stats)
           (set volume (gsl-rstat-mean stats))
           (cond
@@ -105,7 +108,7 @@
         (let* ((neval int 0))
           (random-direction-vector r x)
           (gsl-rstat-add (integral-per-direction integrand x r dimension
-                                                 (extent-oracle x) rtol (address-of neval))
+                                                 (invoke-extent-oracle extent-oracle r x) rtol (address-of neval))
                          stats))
         (set integral (gsl-rstat-mean stats))
         (set error (/ (* CONFIDENCE-INTERVAL-FACTOR (gsl-rstat-sd-mean stats))
@@ -126,7 +129,7 @@
       (with-vector x dimension
         (for-i i number-of-samples
           (hollow-cone-random-vector r mean theta-min theta-max x)
-          (gsl-rstat-add (exp (+ vn (* dimension (log (extent-oracle x)))))
+          (gsl-rstat-add (exp (+ vn (* dimension (log (invoke-extent-oracle extent-oracle r x)))))
                          stats)))
       (cond
        (variance (set (pointer-get variance) (gsl-rstat-variance stats))))
diff --git a/src/oracles.sc b/src/oracles.sc
index e0e88ba..f0f53a8 100644
--- a/src/oracles.sc
+++ b/src/oracles.sc
@@ -3,36 +3,45 @@
 (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-generator r p r0 r1) (double (const gsl-rng*) double double double)
-  (return (if* (gsl-ran-bernoulli r p) r1 r0)))
+(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 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 (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-generator r a b) (double (const gsl-rng*) double double)
-  (return (gsl-ran-flat r a b)))
+(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 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)))
+(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 (/ (+ alpha r)
-                   (+ alpha beta r))))
+      (set* vol (/ (+ (: params alpha) r)
+                   (+ (: params alpha) (: params beta) r))))
     (return vol)))
 
 (define (infinity-norm x) ((static double) (const gsl-vector*))
@@ -43,45 +52,53 @@
       (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))))
+(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)
-  (/ (- (/ edge 2)
+  (/ (- (/ (: params edge) 2)
         (* (GSL-SIGN (gsl-vector-get x i))
-           (gsl-vector-get center i)))
+           (gsl-vector-get (: params 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)))
+(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 (: center size)
+    (for-i i (: (: params 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-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 x axes) (double (const gsl-vector*) (const gsl-vector*))
-  (let* ((k double 0))
-    (for-i i (: axes size)
+(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 axes i)))))
+                            (gsl-vector-get (: params 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)))
+(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 x eccentricity) (double (const gsl-vector*) double)
-  (let* ((xsub gsl-vector-const-view
+(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) eccentricity))))))))
+                        (gsl-pow-2 (/ (gsl-vector-get x 0)
+                                      (: params eccentricity)))))))))
 
-(define (spheroid-true-volume eccentricity dimension) (double double (unsigned int))
-  (return (* (volume-of-ball dimension) 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)))))