From ec527081eaaec124c070288befc0815b139ee232 Mon Sep 17 00:00:00 2001
From: Arun Isaac
Date: Wed, 30 Jun 2021 14:46:38 +0530
Subject: Implement window stopping criterion for integrals.

* src/extent-sampling.sc (integral): Implement window stopping
criterion.
---
 src/extent-sampling.sc | 24 +++++++++++++++---------
 1 file changed, 15 insertions(+), 9 deletions(-)

(limited to 'src')

diff --git a/src/extent-sampling.sc b/src/extent-sampling.sc
index 0c0708c..d46f31d 100644
--- a/src/extent-sampling.sc
+++ b/src/extent-sampling.sc
@@ -78,15 +78,21 @@
     (return (* (surface-area-of-ball n)
                result))))
 
-(define (integral integrand extent-oracle true-integral r dimension rtol stats)
-  (void (const integrand-t*) (const extent-oracle-t*) double (const gsl-rng*) (unsigned int) double gsl-rstat-workspace*)
-  (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 dimension
-                                             (invoke-extent-oracle extent-oracle r x) rtol)
-                     stats))))
+(pre-let* (WINDOW-LENGTH 1000)
+  (define (integral integrand extent-oracle true-integral r dimension rtol stats)
+    (void (const integrand-t*) (const extent-oracle-t*) double (const gsl-rng*) (unsigned int) double gsl-rstat-workspace*)
+    (define accurate-estimates int 0)
+    (with-vector x dimension
+      (do-while (< accurate-estimates WINDOW-LENGTH)
+        (random-direction-vector r x)
+        (gsl-rstat-add (integral-per-direction integrand x dimension
+                                               (invoke-extent-oracle extent-oracle r x) rtol)
+                       stats)
+        (cond
+         ((rtol? (gsl-rstat-mean stats) true-integral rtol)
+          (set+ accurate-estimates 1))
+         (else
+          (set accurate-estimates 0)))))))
 
 (define (volume-cone extent-oracle r mean omega-min omega-max number-of-samples variance)
   (double (const extent-oracle-t*) (const gsl-rng*) (const gsl-vector*) double double (unsigned int) double*)
-- 
cgit v1.2.3