diff options
| -rw-r--r-- | include/extent-sampling.h | 6 | ||||
| -rw-r--r-- | include/oracles.h | 51 | ||||
| -rw-r--r-- | src/extent-sampling.sc | 11 | ||||
| -rw-r--r-- | src/oracles.sc | 109 |
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))))) |