about summary refs log tree commit diff
diff options
context:
space:
mode:
authorArun Isaac2021-04-20 16:52:35 +0530
committerArun Isaac2021-04-20 17:01:55 +0530
commit9266f4ce5213ea3197484ca22cb82ca082382fc6 (patch)
treec96f47466097bac01a63bded1e6ff78a83e49fff
parent2eacc2841e34c926b9bb739c9800783b61e1d01a (diff)
downloadnsmc-9266f4ce5213ea3197484ca22cb82ca082382fc6.tar.gz
nsmc-9266f4ce5213ea3197484ca22cb82ca082382fc6.tar.lz
nsmc-9266f4ce5213ea3197484ca22cb82ca082382fc6.zip
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.
-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*)