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))