Implement offcenter volume experiments.

* src/volume-bodies.sc: Implement offcenter volume experiments.

1 files changed, 85 insertions, 35 deletions

diff --git a/src/volume-bodies.sc b/src/volume-bodies.sc
index f4ee2a3..c53d0b3 100644
--- a/src/volume-bodies.sc
+++ b/src/volume-bodies.sc
@@ -2,26 +2,36 @@
 
 (pre-include "stdio.h")
 (pre-include "stdlib.h")
+(pre-include "unistd.h")
 (pre-include "extent-sampling.h")
+(pre-include "nd-random.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)
+(pre-define TRIALS 10)
+(pre-define DIMENSION-START 10)
+(pre-define DIMENSION-STOP 100)
+(pre-define DIMENSION-STEP 10)
+(pre-define CUBE-OFFSET-MAXIMUM 0.03125)
+(pre-define ELLIPSOID-AXES-START 0.25)
+(pre-define ELLIPSOID-AXES-END 0.5)
+(pre-define ELLIPSOID-OFFSET-MAXIMUM 0.05)
 
-(sc-define-syntax* (with-cube-oracle edge body ...)
+(sc-define-syntax* (with-cube-oracle rng edge center dimension body ...)
   (let ((params (sc-gensym)))
     `(begin
        (declare ,params (struct-variable cube-params
                                          (edge ,edge)
-                                         (center NULL)))
+                                         (center ,center)))
        (declare oracle (struct-variable extent-oracle-t
-                                        (oracle cube-extent-oracle)
+                                        (oracle ,(if (eq? center 'NULL)
+                                                     'cube-extent-oracle
+                                                     'cube-extent-oracle-with-center))
                                         (params (address-of ,params))))
        ,@body)))
 
-(sc-define-syntax* (with-ellipsoid-oracle axes-start axes-end dimension body ...)
+(sc-define-syntax* (with-ellipsoid-oracle axes-start axes-end center dimension body ...)
   (let ((params (sc-gensym))
         (axes (sc-gensym)))
     `(with-vector ,axes ,dimension
@@ -29,45 +39,85 @@
          (gsl-vector-set ,axes i (+ ,axes-start (* i (/ (- ,axes-end ,axes-start)
                                                         (- ,dimension 1))))))
        (declare ,params (struct-variable ellipsoid-params
-                                         (axes ,axes)))
+                                         (axes ,axes)
+                                         (center ,center)))
        (declare oracle (struct-variable extent-oracle-t
-                                        (oracle ellipsoid-extent-oracle)
+                                        (oracle ,(if (eq? center 'NULL)
+                                                     'ellipsoid-extent-oracle
+                                                     'ellipsoid-extent-oracle-with-center))
                                         (params (address-of ,params))))
        ,@body)))
 
-(sc-define-syntax* (experiment-dimension dimension true-volume filename)
+(sc-define-syntax* (experiment-dimension dimension trials true-volume filename)
   `(begin
+     (fprintf stderr "Dimension %d: " dimension)
      (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))
+       (for-i trial ,trials
+         (fprintf stderr "%d " trial)
+         (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)))
+       (fprintf stderr "\n")
        (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 (log-to-data-file filename body ...)
+  (begin
+    (with-data-file fp filename "w"
+      (fprintf fp "dimension\tsamples\n"))
+    body ...))
 
-(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))))))
+(sc-define-syntax (with-center rng center dimension maximum-offset body ...)
+  (with-vector center dimension
+    (random-vector-in-sphere rng CUBE-OFFSET-MAXIMUM center)
+    body ...))
+
+(sc-define-syntax (cube-experiment filename)
+  (log-to-data-file filename
+    (for-i-step dimension DIMENSION-START DIMENSION-STOP DIMENSION-STEP
+      (with-cube-oracle rng EDGE NULL dimension
+        (experiment-dimension dimension 1 cube-true-volume filename)))))
+
+(sc-define-syntax (offcenter-cube-experiment filename)
+  (log-to-data-file filename
+    (for-i-step dimension DIMENSION-START DIMENSION-STOP DIMENSION-STEP
+      (with-center rng center dimension CUBE-OFFSET-MAXIMUM
+        (with-cube-oracle rng EDGE center dimension
+          (experiment-dimension dimension TRIALS cube-true-volume filename))))))
+
+(sc-define-syntax (ellipsoid-experiment filename)
+  (log-to-data-file filename
+    (for-i-step dimension DIMENSION-START DIMENSION-STOP DIMENSION-STEP
+      (with-ellipsoid-oracle ELLIPSOID-AXES-START ELLIPSOID-AXES-END NULL dimension
+        (experiment-dimension dimension 1 ellipsoid-true-volume filename)))))
+
+(sc-define-syntax (offcenter-ellipsoid-experiment filename)
+  (log-to-data-file filename
+    (for-i-step dimension DIMENSION-START DIMENSION-STOP DIMENSION-STEP
+      (with-center rng center dimension ELLIPSOID-OFFSET-MAXIMUM
+        (with-ellipsoid-oracle ELLIPSOID-AXES-START ELLIPSOID-AXES-END center dimension
+          (experiment-dimension dimension TRIALS ellipsoid-true-volume filename))))))
 
 (define (main) (int)
   (with-rng rng
-    (cube-experiment)
-    (ellipsoid-experiment)))
+    (ellipsoid-experiment "volume-ellipsoid.dat")
+    (offcenter-ellipsoid-experiment "volume-offcenter-ellipsoid.dat")
+    (cube-experiment "volume-cube.dat")
+    (offcenter-cube-experiment "volume-offcenter-cube.dat")
+    (with-data-file fp "volume-body-parameters.tex" "w"
+      (fprintf fp "\\\\newcommand{\\\\bodyrtol}{%.1f}\n" RTOL)
+      (fprintf fp "\\\\newcommand{\\\\cubeedge}{%.1f}\n" EDGE)
+      ;; The code uses semi-axes (radius) whereas we wish to report
+      ;; axes (diameter). So, we double.
+      (fprintf fp "\\\\newcommand{\\\\ellipsoidaxesstart}{%.1f}\n" (* 2 ELLIPSOID-AXES-START))
+      (fprintf fp "\\\\newcommand{\\\\ellipsoidaxesstop}{%.1f}\n" (* 2 ELLIPSOID-AXES-END))
+      (fprintf fp "\\\\newcommand{\\\\cubeoffsetmaximumpercent}{%.2f}\n" (* 100 (/ (* 2 CUBE-OFFSET-MAXIMUM) EDGE)))
+      (fprintf fp "\\\\newcommand{\\\\ellipsoidoffsetmaximumpercent}{%.0f}\n" (* 100 (/ ELLIPSOID-OFFSET-MAXIMUM
+                                                                                        ELLIPSOID-AXES-END)))
+      (fprintf fp "\\\\newcommand{\\\\offsettrials}{%d}\n" TRIALS))))