diff options
| -rw-r--r-- | src/volume-bodies.sc | 120 |
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)))) |