Implement volume of bodies experiments.

* src/volume-bodies.sc: New file.
* CMakeLists.txt: Add volume-bodies executable.

2 files changed, 75 insertions, 0 deletions

diff --git a/CMakeLists.txt b/CMakeLists.txt
index e56757d..ad0bc1f 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -69,6 +69,8 @@ install(TARGETS nsmc LIBRARY
   PUBLIC_HEADER DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/nsmc)
 
 # Build executables
+add_executable(volume-bodies volume-bodies.c)
+target_link_libraries(volume-bodies nsmc)
 add_executable(integral integral.c)
 target_link_libraries(integral nsmc -lm)
 # Build and install scheme wrapper.
diff --git a/src/volume-bodies.sc b/src/volume-bodies.sc
new file mode 100644
index 0000000..f4ee2a3
--- /dev/null
+++ b/src/volume-bodies.sc
@@ -0,0 +1,73 @@
+(sc-include "macros/macros")
+
+(pre-include "stdio.h")
+(pre-include "stdlib.h")
+(pre-include "extent-sampling.h")
+(pre-include "oracles.h")
+
+(pre-define EDGE 1.0)
+(pre-define RTOL 0.1)
+(pre-define ELLIPSOID-AXES-START 0.5)
+(pre-define ELLIPSOID-AXES-END 1.0)
+
+(sc-define-syntax* (with-cube-oracle edge body ...)
+  (let ((params (sc-gensym)))
+    `(begin
+       (declare ,params (struct-variable cube-params
+                                         (edge ,edge)
+                                         (center NULL)))
+       (declare oracle (struct-variable extent-oracle-t
+                                        (oracle cube-extent-oracle)
+                                        (params (address-of ,params))))
+       ,@body)))
+
+(sc-define-syntax* (with-ellipsoid-oracle axes-start axes-end dimension body ...)
+  (let ((params (sc-gensym))
+        (axes (sc-gensym)))
+    `(with-vector ,axes ,dimension
+       (for-i i ,dimension
+         (gsl-vector-set ,axes i (+ ,axes-start (* i (/ (- ,axes-end ,axes-start)
+                                                        (- ,dimension 1))))))
+       (declare ,params (struct-variable ellipsoid-params
+                                         (axes ,axes)))
+       (declare oracle (struct-variable extent-oracle-t
+                                        (oracle ellipsoid-extent-oracle)
+                                        (params (address-of ,params))))
+       ,@body)))
+
+(sc-define-syntax* (experiment-dimension dimension true-volume filename)
+  `(begin
+     (with-rstats samples
+       (with-rstats stats
+         (volume (address-of oracle)
+                 (,true-volume dimension (struct-get oracle params))
+                 rng
+                 dimension
+                 RTOL
+                 stats)
+         (gsl-rstat-add (gsl-rstat-n stats)
+                        samples))
+       (with-data-file fp ,filename "a"
+         (fprintf fp "%d\t%g\n" dimension (gsl-rstat-mean samples))))))
+
+(sc-define-syntax* (cube-experiment)
+  (let* ((filename "volume-cube.dat"))
+    `(begin
+       (with-data-file fp ,filename "w"
+         (fprintf fp "dimension\tsamples\n"))
+       (for-i-step dimension 10 100 10
+         (with-cube-oracle EDGE (experiment-dimension dimension cube-true-volume ,filename))))))
+
+(sc-define-syntax* (ellipsoid-experiment)
+  (let* ((filename "volume-ellipsoid.dat"))
+    `(begin
+       (with-data-file fp ,filename "w"
+         (fprintf fp "dimension\tsamples\n"))
+       (for-i-step dimension 10 100 10
+         (with-ellipsoid-oracle ELLIPSOID-AXES-START ELLIPSOID-AXES-END dimension
+                                (experiment-dimension dimension ellipsoid-true-volume ,filename))))))
+
+(define (main) (int)
+  (with-rng rng
+    (cube-experiment)
+    (ellipsoid-experiment)))