From e1538eee83ba5d7ce9a84453434473e7a4ed8f8a Mon Sep 17 00:00:00 2001
From: Arun Isaac
Date: Wed, 30 Jun 2021 14:54:32 +0530
Subject: Implement offcenter volume experiments.
* src/volume-bodies.sc: Implement offcenter volume experiments.
---
src/volume-bodies.sc | 120 ++++++++++++++++++++++++++++++++++++---------------
1 file changed, 85 insertions(+), 35 deletions(-)
(limited to 'src')
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))))
--
cgit v1.2.3