;; [float] -> float -> int
(define w-index
  (lambda (w n)
    (find-index 
     (lambda (e) 
       (< n e)) 
     (tail (differentiate w)))))

;; [a] -> [float] -> a
(define w-choose
  (lambda (l w)
    (list-ref l (w-index w (random 0 1)))))

;; [float] -> [float]
(define normalize-sum
  (lambda (l)
    (let ((n (foldl1 + l)))
      (map 
       (lambda (e) 
	 (/ e n)) 
       l))))

;; [a] -> [float] -> a
(define p-choose
  (lambda (l p)
    (w-choose l (normalize-sum p))))

;; #[a] -> a
(define vector-choose
  (lambda (v)
    (vector-ref v (i-random 0 (vector-length v)))))

;; () -> float
(define random-linear
  (lambda () 
    (min (random 0 1) 
	 (random 0 1))))

;; () -> float
(define random-inverse-linear
  (lambda ()
    (max (random 0 1) 
         (random 0 1))))

;; () -> float
(define random-triangular
  (lambda () 
    (* 0.5 
       (+ (random 0 1) 
          (random 0 1)))))

;; float -> float
(define random-exponential
  (lambda (l)
    (let ((u (random 0 1)))
      (if (zero? u)
          (random-exponential l)
          (/ (- (log u)) l)))))

;; float -> float
(define random-bilinear-exponential
  (lambda (l)
    (let ((u (* 2 (random 0 1))))
      (if (zero? u)
          (random-bilinear-exponential l)
          (* (if (> u 1) -1 1) 
             (log (if (> u 1) 
                      (- 2 u) 
                      u)))))))

;; float -> float -> float
(define random-guassian
  (lambda (sigma mu)
    (let* ((n 12)
           (s (sum (replicate-m n (random 0 1)))))
      (+ (* sigma 
            (/ 1 (sqrt (/ n 12)))
            (- s (/ n 2))) 
         mu))))

;; float -> float
(define random-cauchy
  (lambda (alpha)
    (let ((u (random 0 1)))
      (if (= u 0.5)
          (random-cauchy alpha)
          (* alpha (tan (* u pi)))))))

;; float -> float -> float
(define random-beta
  (lambda (a b)
    (let ((u1 (random 0 1)))
      (if (zero? u1)
          (random-beta a b)
          (let ((u2 (random 0 1)))
            (if (zero? u2)
                (random-beta a b)
                (let* ((y1 (expt u1 (/ 1 a)))
                       (y2 (expt u2 (/ 1 b)))
                       (sum (+ y1 y2)))
                  (if (> sum 1)
                      (random-beta a b)
                      (/ y1 sum)))))))))

;; float -> float -> float
(define random-weibull
  (lambda (s t)
    (let ((u (random 0 1)))
      (if (or (zero? u)
              (= u 1))
          (random-weibull s t)
          (let ((a (/ 1 (- 1 u))))
            (* s (expt (log a) (/ 1 t))))))))

;; float -> float
(define random-poisson
  (lambda (l)
    (let loop ((v (exp (- l)))
               (u (random 0 1))
               (n 0))
      (if (>= u v)
          (loop v (* u (random 0 1)) (+ n 1))
          n))))

;; [a] -> [a]
(define shuffle
  (lambda (l)
    (let ((q (map (lambda (e) (cons (random 0 1) e)) l))
          (c (lambda (a b) (compare (car a) (car b)))))
      (map cdr (sort-by c q)))))