From 9266f4ce5213ea3197484ca22cb82ca082382fc6 Mon Sep 17 00:00:00 2001
From: Arun Isaac
Date: Tue, 20 Apr 2021 16:52:35 +0530
Subject: Return the first accurate estimate.

* include/extent-sampling.h (volume_window): Delete function.
* src/extent-sampling.sc (integral, volume): Return the first accurate
estimate.
(volume-window): Delete function.
* scm/nsmc/wrap.scm (volume-window): Delete function.
---
 src/extent-sampling.sc | 80 ++++++++++++++------------------------------------
 1 file changed, 22 insertions(+), 58 deletions(-)

(limited to 'src')

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*)
-- 
cgit v1.2.3