Implement spheroid importance sampling experiments.

* src/spheroid.sc: New file.
* CMakeLists.txt: Add spheroid executable.

2 files changed, 62 insertions, 0 deletions

diff --git a/CMakeLists.txt b/CMakeLists.txt
index ad0bc1f..b4fae28 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -73,6 +73,9 @@ add_executable(volume-bodies volume-bodies.c)
 target_link_libraries(volume-bodies nsmc)
 add_executable(integral integral.c)
 target_link_libraries(integral nsmc -lm)
+add_executable(spheroid spheroid.c)
+target_link_libraries(spheroid nsmc)
+
 # Build and install scheme wrapper.
 if(${GUILE_FOUND})
   configure_file(scm/nsmc/load-libs.scm.in scm/nsmc/load-libs.scm)
diff --git a/src/spheroid.sc b/src/spheroid.sc
new file mode 100644
index 0000000..830fa7c
--- /dev/null
+++ b/src/spheroid.sc
@@ -0,0 +1,59 @@
+(sc-include "macros/macros")
+
+(pre-include "stdio.h")
+(pre-include "stdlib.h")
+(pre-include "gsl/gsl_blas.h")
+(pre-include "extent-sampling.h")
+(pre-include "oracles.h")
+(pre-include "utils.h")
+
+(pre-define DIMENSION-START 10)
+(pre-define DIMENSION-STEP 10)
+(pre-define DIMENSION-STOP 100)
+
+(pre-define RTOL 0.1)
+(pre-define SAMPLES-PER-CONE 1000)
+(pre-define SOLID-ANGLE-FACTOR 1.1)
+
+(define (experiment eccentricity solid-angle-threshold-exponent-factor fp)
+  (void double double FILE*)
+  (declare params (struct-variable spheroid-params
+                                   (eccentricity eccentricity)))
+  (declare oracle (struct-variable extent-oracle-t
+                                   (oracle spheroid-extent-oracle)
+                                   (params (address-of params))))
+  (fprintf fp "dimension\tsamples\n")
+  (with-rng r
+    (for-i-step dim DIMENSION-START DIMENSION-STOP DIMENSION-STEP
+      (with-vector mean dim
+        (gsl-vector-set-basis mean 0)
+        (gsl-vector-scale mean (/ (gsl-blas-dnrm2 mean)))
+        (declare samples (unsigned int))
+        (define volume double (* 2 (volume-importance (address-of oracle)
+                                                      r
+                                                      mean
+                                                      SAMPLES-PER-CONE
+                                                      SOLID-ANGLE-FACTOR
+                                                      solid-angle-threshold-exponent-factor
+                                                      (address-of samples))))
+        (let* ((true-volume double (spheroid-true-volume dim (address-of params))))
+          (unless (rtol? volume true-volume RTOL)
+            (fprintf stderr
+                     "Volume estimation failed to converge: dimension=%d, true volume=%g, estimated volume=%g, rtol=%g\n"
+                     dim
+                     true-volume
+                     volume
+                     RTOL)
+            (exit EXIT_FAILURE)))
+        (fprintf fp "%d\t%d\n" dim samples)
+        (fflush NULL)))))
+
+(sc-define-syntax* (run-experiment eccentricity solid-angle-threshold-exponent-factor)
+  `(begin
+     (with-data-file fp ,(format #f "spheroid-eccentricity-~a.dat" eccentricity) "w"
+      (experiment ,eccentricity ,solid-angle-threshold-exponent-factor fp))))
+
+(define (main) (int)
+  (run-experiment 10 2)
+  (run-experiment 100 6)
+  (return 0))