(sc-include "macros/macros")

(pre-include "stdio.h")
(pre-include "stdlib.h")
(pre-include "unistd.h")
(pre-include "extent-sampling.h")
(pre-include "nd-random.h")
(pre-include "oracles.h")

(pre-define EDGE 1.0)
(pre-define RTOL 0.1)
(pre-define TRIALS 10)
(pre-define DIMENSION-START 10)
(pre-define DIMENSION-STOP 100)
(pre-define DIMENSION-STEP 10)
(pre-define CUBE-OFFSET-MAXIMUM 0.03125)
(pre-define ELLIPSOID-AXES-START 0.25)
(pre-define ELLIPSOID-AXES-END 0.5)
(pre-define ELLIPSOID-OFFSET-MAXIMUM 0.05)

(sc-define-syntax* (with-cube-oracle rng edge center dimension body ...)
  (let ((params (sc-gensym)))
    `(begin
       (declare ,params (struct-variable cube-params
                                         (edge ,edge)
                                         (center ,center)))
       (declare oracle (struct-variable extent-oracle-t
                                        (oracle ,(if (eq? center 'NULL)
                                                     'cube-extent-oracle
                                                     'cube-extent-oracle-with-center))
                                        (params (address-of ,params))))
       ,@body)))

(sc-define-syntax* (with-ellipsoid-oracle axes-start axes-end center dimension body ...)
  (let ((params (sc-gensym))
        (axes (sc-gensym)))
    `(with-vector ,axes ,dimension
       (for-i i ,dimension
         (gsl-vector-set ,axes i (+ ,axes-start (* i (/ (- ,axes-end ,axes-start)
                                                        (- ,dimension 1))))))
       (declare ,params (struct-variable ellipsoid-params
                                         (axes ,axes)
                                         (center ,center)))
       (declare oracle (struct-variable extent-oracle-t
                                        (oracle ,(if (eq? center 'NULL)
                                                     'ellipsoid-extent-oracle
                                                     'ellipsoid-extent-oracle-with-center))
                                        (params (address-of ,params))))
       ,@body)))

(sc-define-syntax* (experiment-dimension dimension trials true-volume filename)
  `(begin
     (fprintf stderr "Dimension %d: " dimension)
     (with-rstats samples
       (for-i trial ,trials
         (fprintf stderr "%d " trial)
         (with-rstats stats
           (volume (address-of oracle)
                   (,true-volume dimension (struct-get oracle params))
                   rng
                   dimension
                   RTOL
                   stats)
           (gsl-rstat-add (gsl-rstat-n stats)
                          samples)))
       (fprintf stderr "\n")
       (with-data-file fp ,filename "a"
         (fprintf fp "%d\t%g\n" dimension (gsl-rstat-mean samples))))))

(sc-define-syntax (log-to-data-file filename body ...)
  (begin
    (with-data-file fp filename "w"
      (fprintf fp "dimension\tsamples\n"))
    body ...))

(sc-define-syntax (with-center rng center dimension maximum-offset body ...)
  (with-vector center dimension
    (random-vector-in-sphere rng CUBE-OFFSET-MAXIMUM center)
    body ...))

(sc-define-syntax (cube-experiment filename)
  (log-to-data-file filename
    (for-i-step dimension DIMENSION-START DIMENSION-STOP DIMENSION-STEP
      (with-cube-oracle rng EDGE NULL dimension
        (experiment-dimension dimension 1 cube-true-volume filename)))))

(sc-define-syntax (offcenter-cube-experiment filename)
  (log-to-data-file filename
    (for-i-step dimension DIMENSION-START DIMENSION-STOP DIMENSION-STEP
      (with-center rng center dimension CUBE-OFFSET-MAXIMUM
        (with-cube-oracle rng EDGE center dimension
          (experiment-dimension dimension TRIALS cube-true-volume filename))))))

(sc-define-syntax (ellipsoid-experiment filename)
  (log-to-data-file filename
    (for-i-step dimension DIMENSION-START DIMENSION-STOP DIMENSION-STEP
      (with-ellipsoid-oracle ELLIPSOID-AXES-START ELLIPSOID-AXES-END NULL dimension
        (experiment-dimension dimension 1 ellipsoid-true-volume filename)))))

(sc-define-syntax (offcenter-ellipsoid-experiment filename)
  (log-to-data-file filename
    (for-i-step dimension DIMENSION-START DIMENSION-STOP DIMENSION-STEP
      (with-center rng center dimension ELLIPSOID-OFFSET-MAXIMUM
        (with-ellipsoid-oracle ELLIPSOID-AXES-START ELLIPSOID-AXES-END center dimension
          (experiment-dimension dimension TRIALS ellipsoid-true-volume filename))))))

(define (main) (int)
  (with-rng rng
    (ellipsoid-experiment "volume-ellipsoid.dat")
    (offcenter-ellipsoid-experiment "volume-offcenter-ellipsoid.dat")
    (cube-experiment "volume-cube.dat")
    (offcenter-cube-experiment "volume-offcenter-cube.dat")
    (with-data-file fp "volume-body-parameters.tex" "w"
      (fprintf fp "\\\\newcommand{\\\\bodyrtol}{%.1f}\n" RTOL)
      (fprintf fp "\\\\newcommand{\\\\cubeedge}{%.1f}\n" EDGE)
      ;; The code uses semi-axes (radius) whereas we wish to report
      ;; axes (diameter). So, we double.
      (fprintf fp "\\\\newcommand{\\\\ellipsoidaxesstart}{%.1f}\n" (* 2 ELLIPSOID-AXES-START))
      (fprintf fp "\\\\newcommand{\\\\ellipsoidaxesstop}{%.1f}\n" (* 2 ELLIPSOID-AXES-END))
      (fprintf fp "\\\\newcommand{\\\\cubeoffsetmaximumpercent}{%.2f}\n" (* 100 (/ (* 2 CUBE-OFFSET-MAXIMUM) EDGE)))
      (fprintf fp "\\\\newcommand{\\\\ellipsoidoffsetmaximumpercent}{%.0f}\n" (* 100 (/ ELLIPSOID-OFFSET-MAXIMUM
                                                                                        ELLIPSOID-AXES-END)))
      (fprintf fp "\\\\newcommand{\\\\offsettrials}{%d}\n" TRIALS))))