(sc-include "macros/macros")

(pre-include "stddef.h")
(pre-include "gsl/gsl_integration.h")
(pre-include "gsl/gsl_rstat.h")
(pre-include "gsl/gsl_vector.h")
(pre-include "nd-random.h")
(pre-include "utils.h")
(pre-include "extent-sampling.h")

(sc-define-syntax (with-vector var size body ...)
  (with-alloc var gsl-vector*
              (gsl-vector-alloc size)
              gsl-vector-free
              body ...))

(sc-define-syntax (with-rstats var body ...)
  (with-alloc var gsl-rstat-workspace*
              (gsl-rstat-alloc) gsl-rstat-free
              body ...))

(sc-define-syntax (with-integration-workspace var intervals body ...)
  (with-alloc var gsl-integration-workspace*
              (gsl-integration-workspace-alloc intervals)
              gsl-integration-workspace-free
              body ...))

(sc-define-syntax (invoke-extent-oracle extent-oracle r x)
  ((struct-get extent-oracle oracle) r x (struct-get extent-oracle params)))

(pre-define CONFIDENCE-INTERVAL-FACTOR 1.96)

(pre-let* (VOLUME-MINIMUM-SAMPLES 100)
  (define (volume extent-oracle true-volume r dimension rtol stats)
    (double extent-oracle-t double (const gsl-rng*) (unsigned int) double gsl-rstat-workspace*)
    (declare volume double)
    (let* ((vn double (ln-volume-of-ball dimension)))
      (with-vector x dimension
        (do-while (or (> (/ (* CONFIDENCE-INTERVAL-FACTOR
                               (gsl-rstat-sd-mean stats))
                            volume)
                         rtol)
                      (> (rerror volume true-volume) rtol)
                      (< (gsl-rstat-n stats) VOLUME-MINIMUM-SAMPLES))
          (random-direction-vector r x)
          (gsl-rstat-add (exp (+ vn (* dimension (log (invoke-extent-oracle extent-oracle r x)))))
                         stats)
          (set volume (gsl-rstat-mean stats)))))
    (return volume)))

(define (volume-window extent-oracle true-volume r dimension rtol number-of-samples)
  (double extent-oracle-t double (const gsl-rng*) (unsigned int) double (unsigned int*))
  (declare volume double)
  (let* ((vn double (ln-volume-of-ball dimension))
         ;; This is the window length used in Volume.m of
         ;; Lovasz-Vempala's code.
         (window-length int (+ (* 4 dimension dimension) 500))
         (accurate-estimates int 0))
    (with-rstats stats
      (with-vector x dimension
        (do-while (< accurate-estimates window-length)
          (random-direction-vector r x)
          (gsl-rstat-add (exp (+ vn (* dimension (log (invoke-extent-oracle extent-oracle r x)))))
                         stats)
          (set volume (gsl-rstat-mean stats))
          (cond
           ((< (rerror volume true-volume) rtol)
            (set+ accurate-estimates 1))
           (else (set accurate-estimates 0)))))
      (set (pointer-get number-of-samples)
           (gsl-rstat-n stats))))
  (return volume))

(pre-let* (INTEGRATION-INTERVALS 1000)
  (define (integral-per-direction integrand direction r n radius rtol neval)
    (double integrand-t (const gsl-vector*) (const gsl-rng*) (unsigned int) double double int*)
    (define (f r params) (double double void*)
      (set+ (pointer-get neval) 1)
      (return (* (gsl-pow-uint r (- n 1))
                 (integrand r direction))))
    (declare gsl-f (struct-variable gsl-function (function (address-of f)))
             result double
             error double)
    (with-integration-workspace w INTEGRATION-INTERVALS
      (gsl-integration-qag (address-of gsl-f)
                           0
                           radius
                           0
                           rtol
                           INTEGRATION-INTERVALS
                           GSL-INTEG-GAUSS15
                           w
                           (address-of result)
                           (address-of error)))
    (return (* (surface-area-of-ball n)
               result))))

(pre-let* (INTEGRAL-MINIMUM-SAMPLES 100)
  (define (integral integrand extent-oracle true-integral r dimension rtol stats)
    (double integrand-t extent-oracle-t double (const gsl-rng*) (unsigned int) double gsl-rstat-workspace*)
    (declare integral double
             error double)
    (with-vector x dimension
      (do-while (or (> error rtol)
                    (> (rerror integral true-integral) rtol)
                    (< (gsl-rstat-n stats) INTEGRAL-MINIMUM-SAMPLES))
        ;; TODO: If this neval is not used anywhere, remove it.
        (let* ((neval int 0))
          (random-direction-vector r x)
          (gsl-rstat-add (integral-per-direction integrand x r dimension
                                                 (invoke-extent-oracle extent-oracle r x) rtol (address-of neval))
                         stats))
        (set integral (gsl-rstat-mean stats))
        (set error (/ (* CONFIDENCE-INTERVAL-FACTOR (gsl-rstat-sd-mean stats))
                      integral))))
    (cond
     ((> error rtol)
      (GSL-ERROR-VAL "Integral failed to converge" GSL-ETOL integral)))
    (return integral)))

(define (volume-cone extent-oracle r mean omega-min omega-max number-of-samples variance)
  (double extent-oracle-t (const gsl-rng*) (const gsl-vector*) double double (unsigned int) double*)
  (declare volume double)
  (let* ((dimension (unsigned int) (: mean size))
         (vn double (ln-volume-of-ball dimension))
         (theta-min double (solid-angle->planar-angle (* omega-min (surface-area-of-ball dimension)) dimension))
         (theta-max double (solid-angle->planar-angle (* omega-max (surface-area-of-ball dimension)) dimension)))
    (with-rstats stats
      (with-vector x dimension
        (for-i i number-of-samples
          (hollow-cone-random-vector r mean theta-min theta-max x)
          (gsl-rstat-add (exp (+ vn (* dimension (log (invoke-extent-oracle extent-oracle r x)))))
                         stats)))
      (cond
       (variance (set (pointer-get variance) (gsl-rstat-variance stats))))
      (set volume (* (gsl-rstat-mean stats)
                     (- omega-max omega-min)))))
  (return volume))

(define (volume-experiment extent-oracle r mean samples-per-cone solid-angle-factor
                           solid-angle-threshold-exponent-factor number-of-samples)
  (double extent-oracle-t (const gsl-rng*) (const gsl-vector*)
          (unsigned int) double double (unsigned int*))
  (define
    volume double 0
    omega-max double 0.5
    N int 0)
  (let* ((dimension int (: mean size)))
    (do-while (> omega-max (pow 2 (- (* dimension solid-angle-threshold-exponent-factor))))
      (let* ((omega-min double (/ omega-max solid-angle-factor)))
        (set+ volume (volume-cone extent-oracle r mean omega-min omega-max samples-per-cone NULL))
        (set+ N samples-per-cone)
        (set omega-max omega-min))))
  (set+ volume (volume-cone extent-oracle r mean 0 omega-max samples-per-cone NULL))
  (set+ N samples-per-cone)
  (cond
   (number-of-samples (set (pointer-get number-of-samples) N)))
  (return volume))