diff options
Diffstat (limited to 'src/extent-sampling.sc')
| -rw-r--r-- | src/extent-sampling.sc | 155 |
1 files changed, 155 insertions, 0 deletions
diff --git a/src/extent-sampling.sc b/src/extent-sampling.sc new file mode 100644 index 0000000..c6220e6 --- /dev/null +++ b/src/extent-sampling.sc @@ -0,0 +1,155 @@ +(sc-include "macros/macros") + +(pre-include "stddef.h") +(pre-include "gsl/gsl_integration.h") +(pre-include "gsl/gsl_rstat.h") +(pre-include "gsl/gsl_vector.h") +(pre-include "nd-random.h") +(pre-include "utils.h") +(pre-include "extent-sampling.h") + +(sc-define-syntax (with-vector var size body ...) + (with-alloc var gsl-vector* + (gsl-vector-alloc size) + gsl-vector-free + body ...)) + +(sc-define-syntax (with-rstats var body ...) + (with-alloc var gsl-rstat-workspace* + (gsl-rstat-alloc) gsl-rstat-free + body ...)) + +(sc-define-syntax (with-integration-workspace var intervals body ...) + (with-alloc var gsl-integration-workspace* + (gsl-integration-workspace-alloc intervals) + gsl-integration-workspace-free + body ...)) + +(pre-define CONFIDENCE-INTERVAL-FACTOR 1.96) + +(pre-let* (VOLUME-MINIMUM-SAMPLES 100) + (define (volume extent-oracle true-volume r dimension rtol stats) + (double extent-oracle-t double (const gsl-rng*) (unsigned int) double gsl-rstat-workspace*) + (declare volume double) + (let* ((vn double (ln-volume-of-ball dimension))) + (with-vector x dimension + (do-while (or (> (/ (* CONFIDENCE-INTERVAL-FACTOR + (gsl-rstat-sd-mean stats)) + volume) + rtol) + (> (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))))) + stats) + (set volume (gsl-rstat-mean stats))))) + (return volume))) + +(define (volume-window extent-oracle true-volume r dimension rtol number-of-samples) + (double extent-oracle-t double (const gsl-rng*) (unsigned int) double (unsigned int*)) + (declare volume double) + (let* ((vn double (ln-volume-of-ball dimension)) + ;; This is the window length used in Volume.m of + ;; Lovasz-Vempala's code. + (window-length int (+ (* 4 dimension dimension) 500)) + (accurate-estimates int 0)) + (with-rstats stats + (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))))) + stats) + (set volume (gsl-rstat-mean stats)) + (cond + ((< (rerror volume true-volume) rtol) + (set+ accurate-estimates 1)) + (else (set accurate-estimates 0))))) + (set (pointer-get number-of-samples) + (gsl-rstat-n stats)))) + (return volume)) + +(pre-let* (INTEGRATION-INTERVALS 1000) + (define (integral-per-direction integrand direction r n radius rtol neval) + (double integrand-t (const gsl-vector*) (const gsl-rng*) (unsigned int) double double int*) + (define (f r params) (double double void*) + (set+ (pointer-get neval) 1) + (return (* (gsl-pow-uint r (- n 1)) + (integrand r direction)))) + (declare gsl-f (struct-variable gsl-function (function (address-of f))) + result double + error double) + (with-integration-workspace w INTEGRATION-INTERVALS + (gsl-integration-qag (address-of gsl-f) + 0 + radius + 0 + rtol + INTEGRATION-INTERVALS + GSL-INTEG-GAUSS15 + w + (address-of result) + (address-of error))) + (return (* (surface-area-of-ball n) + result)))) + +(pre-let* (INTEGRAL-MINIMUM-SAMPLES 100) + (define (integral integrand extent-oracle true-integral r dimension rtol stats) + (double integrand-t extent-oracle-t double (const gsl-rng*) (unsigned int) double gsl-rstat-workspace*) + (declare integral double + error double) + (with-vector x dimension + (do-while (or (> error rtol) + (> (rerror integral true-integral) rtol) + (< (gsl-rstat-n stats) INTEGRAL-MINIMUM-SAMPLES)) + ;; TODO: If this neval is not used anywhere, remove it. + (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)) + stats)) + (set integral (gsl-rstat-mean stats)) + (set error (/ (* CONFIDENCE-INTERVAL-FACTOR (gsl-rstat-sd-mean stats)) + integral)))) + (cond + ((> error rtol) + (GSL-ERROR-VAL "Integral failed to converge" GSL-ETOL integral))) + (return integral))) + +(define (volume-cone extent-oracle r mean omega-min omega-max number-of-samples variance) + (double extent-oracle-t (const gsl-rng*) (const gsl-vector*) double double (unsigned int) double*) + (declare volume double) + (let* ((dimension (unsigned int) (: mean size)) + (vn double (ln-volume-of-ball dimension)) + (theta-min double (solid-angle->planar-angle (* omega-min (surface-area-of-ball dimension)) dimension)) + (theta-max double (solid-angle->planar-angle (* omega-max (surface-area-of-ball dimension)) dimension))) + (with-rstats stats + (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))))) + stats))) + (cond + (variance (set (pointer-get variance) (gsl-rstat-variance stats)))) + (set volume (* (gsl-rstat-mean stats) + (- omega-max omega-min))))) + (return volume)) + +(define (volume-experiment extent-oracle r mean samples-per-cone solid-angle-factor + solid-angle-threshold-exponent-factor number-of-samples) + (double extent-oracle-t (const gsl-rng*) (const gsl-vector*) + (unsigned int) double double (unsigned int*)) + (define + volume double 0 + omega-max double 0.5 + N int 0) + (let* ((dimension int (: mean size))) + (do-while (> omega-max (pow 2 (- (* dimension solid-angle-threshold-exponent-factor)))) + (let* ((omega-min double (/ omega-max solid-angle-factor))) + (set+ volume (volume-cone extent-oracle r mean omega-min omega-max samples-per-cone NULL)) + (set+ N samples-per-cone) + (set omega-max omega-min)))) + (set+ volume (volume-cone extent-oracle r mean 0 omega-max samples-per-cone NULL)) + (set+ N samples-per-cone) + (cond + (number-of-samples (set (pointer-get number-of-samples) N))) + (return volume)) |