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