aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/extent-sampling.sc11
-rw-r--r--src/oracles.sc109
2 files changed, 70 insertions, 50 deletions
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)))))