blob: 64f2df31ac5c08bcd275fe02b6f0d60f4deda013 (
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
60
61
62
63
64
65
66
|
(sc-include "macros/macros")
(pre-include "gsl/gsl_blas.h")
(pre-include "gsl/gsl_vector.h")
(pre-include "extent-sampling.h")
(pre-include "nd-random.h")
(pre-include "oracles.h")
(pre-define DIMENSION 100)
(pre-define CUBE-EDGE 1.0)
(pre-define ELLIPSOID-AXES-START 0.25)
(pre-define ELLIPSOID-AXES-END 0.5)
(pre-define SAMPLES-PER-ITERATION 5000)
(pre-define ITERATIONS 10)
(sc-define-syntax (invoke-extent-oracle extent-oracle r x)
((struct-get extent-oracle oracle) r x (struct-get extent-oracle params)))
(define (find-centroid oracle center filename) (void extent-oracle-t gsl-vector* (const char*))
(with-rng rng
(with-vector new-center DIMENSION
(with-vector x DIMENSION
(with-file fp filename "w"
(fprintf fp "%g\n" (gsl-blas-dnrm2 center))
(for-i iteration ITERATIONS
(gsl-vector-set-zero new-center)
(for-i trial SAMPLES-PER-ITERATION
(random-direction-vector rng x)
(gsl-vector-scale x (invoke-extent-oracle oracle rng x))
(gsl-vector-add new-center x)
(gsl-vector-scale x (/ -1.0 (gsl-blas-dnrm2 x)))
(gsl-vector-scale x (invoke-extent-oracle oracle rng x))
(gsl-vector-add new-center x))
(gsl-vector-scale new-center (/ 1.0 2 SAMPLES-PER-ITERATION))
(gsl-vector-memcpy center new-center)
(fprintf fp "%g\n" (gsl-blas-dnrm2 center))))))))
(define (main) (int)
(with-vector center DIMENSION
(gsl-vector-set-basis center 0)
(gsl-vector-scale center (* 0.25 CUBE-EDGE))
(declare cube-params (struct-variable cube-params
(edge CUBE-EDGE)
(center center)))
(declare cube-oracle (struct-variable extent-oracle-t
(oracle cube-extent-oracle-with-center)
(params (address-of cube-params))))
(find-centroid cube-oracle center "cube.dat")
(gsl-vector-set-basis center 0)
(gsl-vector-scale center (* 0.5 ELLIPSOID-AXES-START))
(with-vector axes DIMENSION
(for-i i DIMENSION
(gsl-vector-set axes i (+ ELLIPSOID-AXES-START (* i (/ (- ELLIPSOID-AXES-END
ELLIPSOID-AXES-START)
(- DIMENSION 1))))))
(declare ellipsoid-params (struct-variable ellipsoid-params
(axes axes)
(center center)))
(declare ellipsoid-oracle (struct-variable extent-oracle-t
(oracle ellipsoid-extent-oracle-with-center)
(params (address-of ellipsoid-params))))
(find-centroid ellipsoid-oracle center "ellipsoid.dat")))
(return 0))
|
