;;;; 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)
  (only type-checks check-natural-integer check-vector check-procedure)
  (only type-errors error-vector)
  random-source
  (srfi 27)
  (srfi 27 distributions)
  srfi-27-vector-support
  srfi-27-distributions-support)

;;;

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

;;;

#; ;incomplete
(define:-type (vector-iota (n fixnum)) -> (vector-of fixnum)
  (import (only vector-lib vector-unfold))
  (vector-unfold values n) )

(: vector-iota! (vector fixnum #!optional fixnum -> void))
;
(define (vector-iota! vec n #!optional (m #f))
  #;(unless (vector? vec) (error 'vector-iota! "not a vector" vec))
  (let ((start (if (not m) 0 n)) (end (if (not m) (vector-length vec) m)))
    #;(unless (<= start end) (error 'vector-iota! "limits swapped" start end))
    #;(unless (<= end (vector-length vec)) (error 'vector-iota! "limits exceeded" end (vector-length vec)))
    (do ((i start (fx+ i 1)))
        ((fx= i end) vec)
      (vector-set! vec i i) ) ) )

;;;

;;

;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