aboutsummaryrefslogtreecommitdiff
path: root/src/spheroid.sc
blob: 830fa7c10818cdf1049b43e534a3866e093a9c6b (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
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))