aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--CMakeLists.txt3
-rw-r--r--src/spheroid.sc59
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))