;;;; srfi-27-vector.scm
;;;; Kon Lovett, Feb '10

; Chicken Generic Arithmetic!

(module srfi.27.vector

(;export
  ;
  make-random-permutations
  make-random-vector
  make-random-hollow-sphere
  make-random-solid-sphere
  ;
  random-permutation!
  random-vector!
  random-hollow-sphere!
  random-solid-sphere!)

(import scheme
  (chicken base)
  (chicken type)
  (chicken fixnum)
  (type-checks-numbers integer)
  type-checks-structured
  type-errors-structured
  random-source
  (srfi 27)
  (srfi 27 distributions)
  srfi-27-vector-support
  srfi-27-distributions-support)

;;;

(include-relative "srfi-27-common-types")

(: make-random-permutations (#!rest -> random-vector-function))
(: random-permutation! (srfi-27-vector #!rest -> void))
(: make-random-vector (#!rest -> random-vector-function))
(: random-vector! (srfi-27-vector #!rest -> void))
(: make-random-hollow-sphere (#!rest -> vector))
(: random-hollow-sphere! (srfi-27-vector #!rest -> void))
(: make-random-solid-sphere (#!rest -> vector))
(: random-solid-sphere! (srfi-27-vector #!rest -> void))

;;;

(: *vector-iota! (vector fixnum #!optional (or false fixnum) fixnum -> void))

(define (*vector-iota! vec n #!optional (m #f) (s 1))
  (let ((start (if (not m) 0 n)) (end (if (not m) (vector-length vec) m)))
    (do ((i start (fx+ i 1))
         (j start (fx+ j s)))
        ((fx= i end) vec)
      (vector-set! vec i j) ) ) )

;;;

;;

;Knuth's "The Art of Computer Programming", Vol. II, 2nd ed., Algorithm P of
;Section 3.4.2
;
(define (*random-permutation! vec rndint)
  (let ((n (vector-length vec)))
    (*vector-iota! vec n)
    (do ((k n (fx- k 1)))
        ((fx= k 1) vec)
      ;random-swap
      (let* ((i (fx- k 1))
             (j (rndint n))
             (xi (vector-ref vec i))
             (xj (vector-ref vec j)) )
        (vector-set! vec i xj)
        (vector-set! vec j xi) ) ) ) )

(define (make-random-permutations #!key (randoms (current-random-integer)))
  (lambda (n)
    (*random-permutation!
      (make-vector
        (check-natural-integer 'make-random-permutations n 'length)
        0)
      (check-procedure 'make-random-permutations randoms 'randoms))) )

(define (random-permutation! vec #!key (randoms (current-random-integer)))
  (*random-permutation!
    (check-vector 'random-permutation! vec)
    (check-procedure 'random-permutation! randoms 'randoms)) )

;;

(define (make-random-vector #!key (randoms (current-random-real)))
  (lambda (n)
    (vector-filled!
      (make-vector
        (check-natural-integer 'random-vector n 'length))
      (check-procedure 'make-random-vector randoms 'randoms))) )

(define (random-vector! vec #!key (randoms (current-random-real)))
  (vector%-filled!
    (check-vector% 'random-vector! vec)
    (check-procedure 'random-vector! randoms 'randoms)) )

;;; Normal vectors

;;

; Fills vect with inexact real random numbers the sum of whose
; squares is equal to 1.  Thinking of vect as coordinates in space
; of dimension n = (vector-length vect), the coordinates are
; uniformly distributed over the surface of the unit n-sphere.

(define (**random-hollow-sphere! vec norms)
  (vector%-filled! vec norms)
  (vector%-scale! vec (*reciprocal (sqrt (vector%-sum-squares vec))))
  vec )

(define (*random-hollow-sphere! vec mu sigma randoms)
  (let-values (((norms pl)
                 (make-random-normals #:mu mu #:sigma sigma #:randoms randoms)))
    (**random-hollow-sphere! vec norms) ) )

(define (make-random-hollow-sphere #!key
                                   (mu 0)
                                   (sigma 1)
                                   (randoms (current-random-real)))
  (let-values (((norms pl)
                  (make-random-normals #:mu mu #:sigma sigma #:randoms randoms)))
    (lambda (n)
      (**random-hollow-sphere!
        (make-vector (check-natural-integer 'random-hollow-sphere n 'length))
        norms) ) ) )

(define (random-hollow-sphere! vec #!key
                                   (mu 0)
                                   (sigma 1)
                                   (randoms (current-random-real)))
  (*random-hollow-sphere!
    (check-vector% 'random-hollow-sphere! vec)
    mu sigma randoms) )

;;

; Fills vect with inexact real random numbers the sum of whose
; squares is less than 1.  Thinking of vect as coordinates in
; space of dimension n = (vector-length vect), the coordinates are
; uniformly distributed within the unit n-sphere.

(define (**random-solid-sphere! vec randoms norms)
  (**random-hollow-sphere! vec norms)
  (vector%-scale! vec (expt (randoms) (*reciprocal (vector%-length vec))))
  vec )

(define (*random-solid-sphere! vec mu sigma randoms)
  (let-values (((norms pl)
                  (make-random-normals #:mu mu #:sigma sigma #:randoms randoms)))
    (**random-solid-sphere! vec randoms norms) ) )

(define (make-random-solid-sphere #!key
                                  (mu 0)
                                  (sigma 1)
                                  (randoms (current-random-real)))
  (let-values (((norms pl)
                  (make-random-normals #:mu mu #:sigma sigma #:randoms randoms)))
    (lambda (n)
      (**random-solid-sphere!
        (make-vector (check-natural-integer 'random-solid-sphere n 'length))
        randoms norms) ) ) )

(define (random-solid-sphere! vec #!key
                                  (mu 0)
                                  (sigma 1)
                                  (randoms (current-random-real)))
  (*random-solid-sphere!
    (check-vector% 'random-solid-sphere! vec)
    mu sigma randoms) )

) ;module srfi.27.vector