diff options
| -rw-r--r-- | include/extent-sampling.h | 5 | ||||
| -rw-r--r-- | scm/nsmc/wrap.scm | 10 | ||||
| -rw-r--r-- | src/extent-sampling.sc | 80 |
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*) |