aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--include/extent-sampling.h5
-rw-r--r--scm/nsmc/wrap.scm10
-rw-r--r--src/extent-sampling.sc80
3 files changed, 22 insertions, 73 deletions
diff --git a/include/extent-sampling.h b/include/extent-sampling.h
index 9ceddf5..edf8b19 100644
--- a/include/extent-sampling.h
+++ b/include/extent-sampling.h
@@ -43,11 +43,6 @@ double volume
const gsl_rng* r, unsigned int dimension, double rtol,
gsl_rstat_workspace* stats);
-double volume_window
-(extent_oracle_t *extent_oracle, double true_volume,
- const gsl_rng* r, unsigned int dimension, double rtol,
- unsigned int* number_of_samples);
-
double integral
(integrand_t *integrand, extent_oracle_t *extent_oracle, double true_integral,
const gsl_rng* r, unsigned int dimension, double rtol,
diff --git a/scm/nsmc/wrap.scm b/scm/nsmc/wrap.scm
index 33c1bee..005dd86 100644
--- a/scm/nsmc/wrap.scm
+++ b/scm/nsmc/wrap.scm
@@ -369,16 +369,6 @@
true-volume %gsl-random-state dimension rtol stats)
(rstat-n stats)))
-(define-public (volume-window extent-oracle true-volume dimension rtol)
- (let ((samples (make-c-struct (list unsigned-int) (list 0))))
- ((pointer->procedure double
- (dynamic-func "volume_window" lib-nsmc)
- (list '* double '* unsigned-int double '*))
- (maybe-procedure->extent-oracle extent-oracle)
- true-volume %gsl-random-state dimension rtol samples)
- (match (parse-c-struct samples (list unsigned-int))
- ((samples) samples))))
-
(define-public (integral integrand extent-oracle true-integral dimension rtol)
(let ((stats (rstat-alloc)))
((pointer->procedure double
diff --git a/src/extent-sampling.sc b/src/extent-sampling.sc
index 4ed6d33..971f196 100644
--- a/src/extent-sampling.sc
+++ b/src/extent-sampling.sc
@@ -49,45 +49,17 @@
(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 (invoke-extent-oracle extent-oracle r 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*))
+(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))
- ;; 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 (invoke-extent-oracle extent-oracle r 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))))
+ (let* ((vn double (ln-volume-of-ball dimension)))
+ (with-vector x dimension
+ (do-while (> (rerror volume true-volume)
+ rtol)
+ (random-direction-vector r x)
+ (gsl-rstat-add (exp (+ vn (* dimension (log (invoke-extent-oracle extent-oracle r x)))))
+ stats)
+ (set volume (gsl-rstat-mean stats)))))
(return volume))
(sc-define-syntax (invoke-integrand integrand r x)
@@ -116,26 +88,18 @@
(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))
- (random-direction-vector r x)
- (gsl-rstat-add (integral-per-direction integrand x r dimension
- (invoke-extent-oracle extent-oracle r x) rtol)
- 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 (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)
+ (with-vector x dimension
+ (do-while (> (rerror (gsl-rstat-mean stats) true-integral)
+ rtol)
+ (random-direction-vector r x)
+ (gsl-rstat-add (integral-per-direction integrand x r dimension
+ (invoke-extent-oracle extent-oracle r x) rtol)
+ stats)
+ (set integral (gsl-rstat-mean stats))))
+ (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*)