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

  (require-library
    type-checks type-errors
    random-source
    srfi-27-uniform-random srfi-27-distributions
    srfi-27-vector-support)

;;;

; (in case special processing needed near limits TBD)
(define (*reciprocal n) (/ 1.0 n))
(define (*-reciprocal n) (/ -1.0 n))

(define (vector-iota-set! vec n)
  (do ((i 0 (fx+ i 1)))
      ((fx= i n))
    (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 randoms)
  (let ((n (vector-length vec)))
    (vector-iota-set! vec n)
    (do ((k n (fx- k 1)))
        ((fx= k 1) vec)
      (let* ((i (fx- k 1))
             (j (randoms k))
             (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 (make-uniform-random-integers)))
  (lambda (n)
    (check-cardinal-integer 'make-random-permutations n 'length)
    (*random-permutation! (make-vector n 0) randoms)) )

(define (random-permutation! vec #!key (randoms (make-uniform-random-integers)))
  (check-vector 'random-permutation! vec)
  (*random-permutation! vec randoms) )

;;

(define (make-random-vector #!key (randoms (make-uniform-random-reals)))
  (lambda (n)
    (check-cardinal-integer 'random-vector n 'length)
    (vector-filled! (make-vector n) randoms)) )

(define (random-vector! vec #!key (randoms (make-uniform-random-reals)))
  (check-vector% 'random-vector! vec)
  (vector%-filled! vec randoms) )

;;; Normal vectors

;;

; Fills vect with inexact real random numbers the sum of whose
; squares is equal to 1.0.  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)))) )

(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.0) (sigma 1.0) (randoms (make-uniform-random-reals)))
  (let-values (((norms pl) (make-random-normals #:mu mu #:sigma sigma #:randoms randoms)))
    (lambda (n)
      (check-cardinal-integer 'random-hollow-sphere n 'length)
      (**random-hollow-sphere! (make-vector n) norms) ) ) )

(define (random-hollow-sphere! vec #!key (mu 0.0) (sigma 1.0) (randoms (make-uniform-random-reals)))
  (check-vector% 'random-hollow-sphere! vec)
  (*random-hollow-sphere! vec mu sigma randoms) )

;;

; Fills vect with inexact real random numbers the sum of whose
; squares is less than 1.0.  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)))) )

(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.0) (sigma 1.0) (randoms (make-uniform-random-reals)))
  (let-values (((norms pl) (make-random-normals #:mu mu #:sigma sigma #:randoms randoms)))
    (lambda (n)
      (check-cardinal-integer 'random-solid-sphere n 'length)
      (**random-solid-sphere! (make-vector n) randoms norms) ) ) )

(define (random-solid-sphere! vec #!key (mu 0.0) (sigma 1.0) (randoms (make-uniform-random-reals)))
  (check-vector% 'random-solid-sphere! vec)
  (*random-solid-sphere! vec mu sigma randoms) )

) ;module srfi-27-vector