diff options
Diffstat (limited to 'src/centroid.sc')
-rw-r--r-- | src/centroid.sc | 66 |
1 files changed, 66 insertions, 0 deletions
diff --git a/src/centroid.sc b/src/centroid.sc new file mode 100644 index 0000000..64f2df3 --- /dev/null +++ b/src/centroid.sc @@ -0,0 +1,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)) |