(define-module (extent-sampling wrap)
  #:use-module (srfi srfi-1)
  #:use-module (srfi srfi-26)
  #:use-module (ice-9 match)
  #:use-module (system foreign))

(include "load-libs.scm")

(define (make-gsl-alloc allocator-name allocator-args freer-name)
  (lambda args
    (let ((obj
           (apply (pointer->procedure
                   '* (dynamic-func allocator-name lib-gsl) allocator-args)
                  args)))
      (set-pointer-finalizer! obj (dynamic-func freer-name lib-gsl))
      obj)))

(define-public pi
  ((pointer->procedure double
                       (dynamic-func "pi" lib-extentsampling)
                       (list))))

(define-public volume-of-ball
  (pointer->procedure double
                      (dynamic-func "volume_of_ball" lib-extentsampling)
                      (list int)))

(define-public surface-area-of-ball
  (pointer->procedure double
                      (dynamic-func "surface_area_of_ball" lib-extentsampling)
                      (list int)))

(define-public angle-between-vectors
  (pointer->procedure double
                      (dynamic-func "angle_between_vectors" lib-extentsampling)
                      (list '* '*)))

(define-public planar-angle->solid-angle
  (pointer->procedure double
                      (dynamic-func "planar_angle_to_solid_angle" lib-extentsampling)
                      (list double int)))

;; Random state

(define %gsl-random-state
  ((make-gsl-alloc "gsl_rng_alloc" (list '*) "gsl_rng_free")
   (dereference-pointer
    (dynamic-pointer "gsl_rng_default" lib-gsl))))

(define-public (set-gsl-random-state! seed)
  ((pointer->procedure void
                       (dynamic-func "gsl_rng_set" lib-gsl)
                       (list '* unsigned-long))
   %gsl-random-state seed))

(define-public random-flat
  (cut (pointer->procedure double
                           (dynamic-func "gsl_ran_flat" lib-gsl)
                           (list '* double double))
    %gsl-random-state <> <>))

;; Vector functions

(define vector-alloc
  (make-gsl-alloc "gsl_vector_alloc" (list int) "gsl_vector_free"))

(define-public (vector-size vector)
  (match (parse-c-struct vector
                         (list size_t size_t '* '* int))
    ((size _ _ _ _) size)))

(define-public (vector-scale! vector k)
  "Scale VECTOR by a factor K."
  ((pointer->procedure int
                       (dynamic-func "gsl_vector_scale" lib-gsl)
                       (list '* double))
   vector k)
  vector)

(define-public (vector-norm2 vector)
  ((pointer->procedure double
                       (dynamic-func "gsl_blas_dnrm2" lib-gsl)
                       (list '*))
   vector))

(define-public (vector-add! u v)
  "Return the vector sum of vectors U and V. The sum is stored in U."
  ((pointer->procedure int
                       (dynamic-func "gsl_vector_add" lib-gsl)
                       (list '* '*))
   u v)
  u)

(define-public (basis n dimension)
  "Return the basis vector in the Nth canonical direction."
  (let ((vector (vector-alloc dimension)))
    ((pointer->procedure void
                         (dynamic-func "gsl_vector_set_basis" lib-gsl)
                         (list '* int))
     vector n)
    vector))

;; Histogram functions

(define histogram-struct
  (list size_t '* '*))

(define (histogram-alloc bins min max)
  (let ((histogram
         ((make-gsl-alloc "gsl_histogram_alloc" (list size_t) "gsl_histogram_free")
          bins)))
    ((pointer->procedure int
                         (dynamic-func "gsl_histogram_set_ranges_uniform" lib-gsl)
                         (list '* double double))
     histogram min max)
    histogram))

(define-public (histogram->points histogram)
  (match (parse-c-struct histogram histogram-struct)
    ((n range bin)
     (let ((range (array->list (pointer->bytevector range (1+ n) 0 'f64))))
       (map cons
            (drop-right range 1)
            (array->list (pointer->bytevector bin n 0 'f64)))))))

(define-public (histogram->pdf histogram)
  (match (parse-c-struct histogram histogram-struct)
    ((n range bin)
     (let ((range (pointer->bytevector range (1+ n) 0 'f64)))
       ((pointer->procedure int
                            (dynamic-func "gsl_histogram_scale" lib-extentsampling)
                            (list '* double))
        histogram
        (/ n
           (- (array-ref range n)
              (array-ref range 0))
           (fold + 0 (array->list (pointer->bytevector bin n 0 'f64)))))
       histogram))))

(define-public (rhistogram bins min max)
  (case-lambda
    (() (histogram-alloc bins min max))
    ((histogram) histogram)
    ((histogram input)
     ((pointer->procedure int
                          (dynamic-func "gsl_histogram_increment" lib-extentsampling)
                          (list '* double))
      histogram input)
     histogram)))

;; nd-random

(define-public (random-direction-vector dimension)
  (let ((vector (vector-alloc dimension)))
    ((pointer->procedure void
                         (dynamic-func "random_direction_vector" lib-extentsampling)
                         (list '* '*))
     %gsl-random-state vector)
    vector))

(define-public (subsampling-random-vector mean max-theta)
  (let ((vector (vector-alloc (vector-size mean))))
    ((pointer->procedure void
                         (dynamic-func "subsampling_random_vector" lib-extentsampling)
                         (list '* '* double '*))
     %gsl-random-state mean max-theta vector)
    vector))

(define-public (shifted-gaussian-random-vector mean max-theta standard-deviation)
  (let* ((vector (vector-alloc (vector-size mean)))
         (cost
          ((pointer->procedure int
                               (dynamic-func "shifted_gaussian_random_vector" lib-extentsampling)
                               (list '* '* double double '*))
           %gsl-random-state mean max-theta standard-deviation vector)))
    vector))

(define-public (shifted-gaussian-random-vector-cost mean max-theta standard-deviation)
  ((pointer->procedure int
                         (dynamic-func "shifted_gaussian_random_vector" lib-extentsampling)
                         (list '* '* double double '*))
   %gsl-random-state mean max-theta standard-deviation (vector-alloc (vector-size mean))))

(define %integration-workspace
  ((make-gsl-alloc "gsl_integration_workspace_alloc" (list size_t) "gsl_integration_workspace_free")
   1000))

(define-public (shifted-gaussian-pdf theta mean max-theta standard-deviation)
  ((pointer->procedure double
                       (dynamic-func "shifted_gaussian_pdf" lib-extentsampling)
                       (list double double double double unsigned-int '*))
   theta (vector-norm2 mean) max-theta standard-deviation (vector-size mean) %integration-workspace))