init

1 files changed, 551 insertions, 0 deletions

diff --git a/194-impl.scm b/194-impl.scm
new file mode 100644
index 0000000..b175f3a
--- /dev/null
+++ b/194-impl.scm
@@ -0,0 +1,551 @@
+; SPDX-FileCopyrightText: 2020 Arvydas Silanskas
+; SPDX-FileCopyrightText: 2020 Bradley Lucier
+; SPDX-License-Identifier: MIT
+
+;;
+;; Parameters
+;;
+(define current-random-source (make-parameter default-random-source))
+
+(define (with-random-source random-source thunk)
+  (unless (random-source? random-source)
+    (error "expected random source"))
+  (parameterize ((current-random-source random-source))
+                (thunk)))
+
+;;
+;; Carefully return consecutive substreams of the s'th
+;; SRFI 27 stream of random numbers.  See Sections 1.2 and
+;; 1.3 of "An object-oriented random-number package with many
+;; long streams and substreams", by Pierre L'Ecuyer, Richard
+;; Simard, E. Jack Chen, and W. David Kelton, Operations Research,
+;; vol. 50 (2002), pages 1073-1075.
+;; https://doi.org/10.1287/opre.50.6.1073.358
+;;
+
+(define (make-random-source-generator s)
+  (if (not (and (exact? s)
+                (integer? s)
+                (not (negative? s))))
+      (error "make-random-source-generator: Expect nonnegative exact integer argument: " s)
+      (let ((substream 0))
+        (lambda ()
+          (let ((new-source (make-random-source))) ;; deterministic
+            (random-source-pseudo-randomize! new-source s substream)
+            (set! substream (+ substream 1))
+            new-source)))))
+
+;;
+;; Primitive randoms
+;;
+
+(define (make-random-integer-generator low-bound up-bound)
+  (unless (and (integer? low-bound)
+               (exact? low-bound))
+    (error "expected exact integer for lower bound"))
+  (unless (and (integer? up-bound)
+               (exact? up-bound))
+    (error "expected exact integer for upper bound"))
+  (unless (< low-bound up-bound)
+    (error "upper bound should be greater than lower bound"))
+  (let ((rand-int-proc (random-source-make-integers (current-random-source)))
+        (range (- up-bound low-bound)))
+    (lambda ()
+      (+ low-bound (rand-int-proc range)))))
+
+(define (make-random-u1-generator)
+  (make-random-integer-generator 0 2))
+(define (make-random-u8-generator)
+  (make-random-integer-generator 0 256))
+(define (make-random-s8-generator)
+  (make-random-integer-generator -128 128))
+(define (make-random-u16-generator)
+  (make-random-integer-generator 0 65536))
+(define (make-random-s16-generator)
+  (make-random-integer-generator -32768 32768))
+(define (make-random-u32-generator)
+  (make-random-integer-generator 0 (expt 2 32)))
+(define (make-random-s32-generator)
+  (make-random-integer-generator (- (expt 2 31)) (expt 2 31)))
+(define (make-random-u64-generator)
+  (make-random-integer-generator 0 (expt 2 64)))
+(define (make-random-s64-generator)
+  (make-random-integer-generator (- (expt 2 63)) (expt 2 63)))
+
+(define (clamp-real-number lower-bound upper-bound value)
+  (cond ((not (real? lower-bound))
+         (error "expected real number for lower bound"))
+        ((not (real? upper-bound))
+         (error "expected real number for upper bound"))
+        ((not (<= lower-bound upper-bound))
+         (error "lower bound must be <= upper bound"))
+        ((< value lower-bound) lower-bound)
+        ((> value upper-bound) upper-bound)
+        (else value)))
+
+(define (make-random-real-generator low-bound up-bound)
+  (unless (and (real? low-bound)
+               (finite? low-bound))
+    (error "expected finite real number for lower bound"))
+  (unless (and (real? up-bound)
+               (finite? up-bound))
+     (error "expected finite real number for upper bound"))
+  (unless (< low-bound up-bound)
+     (error "lower bound must be < upper bound"))
+  (let ((rand-real-proc (random-source-make-reals (current-random-source))))
+   (lambda ()
+     (define t (rand-real-proc))
+     ;; alternative way of doing lowbound + t * (up-bound - low-bound)
+     ;; is susceptible to rounding errors and would require clamping to be safe
+     ;; (which in turn requires 144 for adjacent float function)
+     (+ (* t low-bound)
+        (* (- 1.0 t) up-bound)))))
+
+(define (make-random-rectangular-generator
+          real-lower-bound real-upper-bound
+          imag-lower-bound imag-upper-bound)
+  (let ((real-gen (make-random-real-generator real-lower-bound real-upper-bound))
+        (imag-gen (make-random-real-generator imag-lower-bound imag-upper-bound)))
+    (lambda ()
+      (make-rectangular (real-gen) (imag-gen)))))
+
+(define make-random-polar-generator
+  (case-lambda
+    ((magnitude-lower-bound magnitude-upper-bound)
+     (make-random-polar-generator 0+0i magnitude-lower-bound magnitude-upper-bound 0 (* 2 PI)))
+    ((origin magnitude-lower-bound magnitude-upper-bound)
+     (make-random-polar-generator origin magnitude-lower-bound magnitude-upper-bound 0 (* 2 PI)))
+    ((magnitude-lower-bound magnitude-upper-bound angle-lower-bound angle-upper-bound)
+     (make-random-polar-generator 0+0i magnitude-lower-bound magnitude-upper-bound angle-lower-bound angle-upper-bound))
+    ((origin magnitude-lower-bound magnitude-upper-bound angle-lower-bound angle-upper-bound)
+     (unless (complex? origin)
+       (error "origin should be complex number"))
+     (unless (and (real? magnitude-lower-bound)
+                  (real? magnitude-upper-bound)
+                  (real? angle-lower-bound)
+                  (real? angle-upper-bound))
+       (error "magnitude and angle bounds should be real numbers"))
+     (unless (and (<= 0 magnitude-lower-bound)
+                  (<= 0 magnitude-upper-bound))
+       (error "magnitude bounds should be positive"))
+     (unless (< magnitude-lower-bound magnitude-upper-bound)
+       (error "magnitude lower bound should be less than upper bound"))
+     (when (= angle-lower-bound angle-upper-bound)
+       (error "angle bounds shouldn't be equal"))
+     (let* ((b (square magnitude-lower-bound))
+            (m (- (square magnitude-upper-bound) b))
+            (t-gen (make-random-real-generator 0. 1.))
+            (phi-gen (make-random-real-generator angle-lower-bound angle-upper-bound)))
+       (lambda ()
+         (let* ((t (t-gen))
+                (phi (phi-gen))
+                (r (sqrt (+ (* m t) b))))
+          (+ origin (make-polar r phi))))))))
+
+(define (make-random-boolean-generator)
+  (define u1 (make-random-u1-generator))
+  (lambda ()
+    (zero? (u1))))
+
+(define (make-random-char-generator str)
+  (when (not (string? str))
+    (error "expected string"))
+  (unless (> (string-length str) 0)
+    (error "given string is of length 0"))
+  (let* ((int-gen (make-random-integer-generator 0 (string-length str))))
+   (lambda ()
+     (string-ref str (int-gen)))))
+
+(define (make-random-string-generator k str)
+  (let ((char-gen (make-random-char-generator str))
+        (int-gen (make-random-integer-generator 0 k)))
+    (lambda ()
+      (generator->string char-gen (int-gen)))))
+
+;;
+;; Non-uniform distributions
+;;
+
+(define PI (* 4 (atan 1.0)))
+
+(define (make-bernoulli-generator p)
+  (unless (real? p)
+    (error "expected p to be real"))
+  (unless (<= 0 p 1)
+    (error "expected 0 <= p <= 1"))
+  (let ((rand-real-proc (random-source-make-reals (current-random-source))))
+   (lambda ()
+     (if (<= (rand-real-proc) p)
+         1
+         0))))
+
+(define (make-categorical-generator weights-vec)
+  (define weight-sum
+    (vector-fold
+      (lambda (sum p)
+        (unless (and (number? p)
+                     (> p 0))
+          (error "parameter must be a vector of positive numbers"))
+        (+ sum p))
+      0
+      weights-vec))
+  (define length (vector-length weights-vec))
+  (let ((real-gen (make-random-real-generator 0 weight-sum)))
+   (lambda ()
+     (define roll (real-gen))
+     (let it ((sum 0)
+              (i 0))
+       (define newsum (+ sum (vector-ref weights-vec i)))
+       (if (or (< roll newsum)
+               ;; in case of rounding errors and no matches, return last element
+               (= i (- length 1)))
+           i
+           (it newsum
+               (+ i 1)))))))
+
+;; Normal distribution (continuous - generates real numbers)
+;; Box-Muller algorithm
+;; NB: We tested Ziggurat method, too,
+;; only to find out Box-Muller is faster about 12% - presumably
+;; the overhead of each ops is larger in Gauche than C/C++, and
+;; so the difference of cost of log or sin from the primitive
+;; addition/multiplication are negligible.
+
+;; NOTE: this implementation is not thread safe
+(define make-normal-generator
+  (case-lambda
+    (()
+     (make-normal-generator 0.0 1.0))
+    ((mean)
+     (make-normal-generator mean 1.0))
+    ((mean deviation)
+     (let ((rand-real-proc (random-source-make-reals (current-random-source)))
+           (state #f))
+       (unless (and (real? mean)
+                    (finite? mean))
+         (error "expected mean to be finite real number"))
+       (unless (and (real? deviation)
+                    (finite? deviation)
+                    (> deviation 0))
+         (error "expected deviation to be positive finite real number"))
+       (lambda ()
+         (if state
+             (let ((result state))
+              (set! state #f)
+              result)
+             (let* ((r (sqrt (* -2 (log (rand-real-proc)))))
+                    (theta (* 2 PI (rand-real-proc))))
+               (set! state (+ mean (* deviation r (cos theta))))
+               (+ mean (* deviation r (sin theta))))))))))
+
+(define (make-exponential-generator mean)
+  (unless (and (real? mean)
+               (finite? mean)
+               (positive? mean))
+    (error "expected mean to be finite positive real number"))
+  (let ((rand-real-proc (random-source-make-reals (current-random-source))))
+   (lambda ()
+     (- (* mean (log (rand-real-proc)))))))
+
+(define (make-geometric-generator p)
+
+  (define (log1p x)
+    ;; Adapted from SRFI 144
+    (let ((u (+ 1.0 x)))
+      (cond ((= u 1.0)
+             x) ;; gets sign of zero result correct
+            ((= u x)
+             (log u)) ;; large arguments and infinities
+            (else
+             (* (log u) (/ x (- u 1.0)))))))
+
+  (unless (and (real? p)
+               (> p 0)
+               (<= p 1))
+          (error "expected p to be real number, 0 < p <= 1"))
+  (if (zero? (- p 1.))
+      ;; p is indistinguishable from 1.
+      (lambda () 1)
+      (let ((c (/ (log1p (- p))))
+            (rand-real-proc (random-source-make-reals (current-random-source))))
+        (lambda ()
+          (exact (ceiling (* c (log (rand-real-proc)))))))))
+
+;; Draw from poisson distribution with mean L, variance L.
+;; For small L, we use Knuth's method.  For larger L, we use rejection
+;; method by Atkinson, The Computer Generation of Poisson Random Variables,
+;; J. of the Royal Statistical Society Series C (Applied Statistics), 28(1),
+;; pp29-35, 1979.  The code here is a port by John D Cook's C++ implementation
+;; (http://www.johndcook.com/stand_alone_code.html )
+
+;; NOTE: this implementation calculates and stores a table of log(n!) on first invocation of L >= 36
+;; and therefore is not entirely thread safe (should still produce correct result, but with performance hit if table
+;; is recalculated multiple times)
+(define (make-poisson-generator L)
+  (unless (and (real? L)
+               (finite? L)
+               (> L 0))
+    (error "expected L to be finite positive real number"))
+  (let ((rand-real-proc (random-source-make-reals (current-random-source))))
+   (if (< L 30)
+       (make-poisson/small rand-real-proc L)
+       (make-poisson/large rand-real-proc L))))
+
+;; private
+(define (make-poisson/small rand-real-proc L)
+  (lambda ()
+    (do ((exp-L (exp (- L)))
+         (k 0 (+ k 1))
+         (p 1.0 (* p (rand-real-proc))))
+        ((<= p exp-L) (- k 1)))))
+
+;; private
+(define (make-poisson/large rand-real-proc L)
+  (let* ((c (- 0.767 (/ 3.36 L)))
+         (beta (/ PI (sqrt (* 3 L))))
+         (alpha (* beta L))
+         (k (- (log c) L (log beta))))
+    (define (loop)
+      (let* ((u (rand-real-proc))
+             (x (/ (- alpha (log (/ (- 1.0 u) u))) beta))
+             (n (exact (floor (+ x 0.5)))))
+        (if (< n 0)
+            (loop)
+            (let* ((v (rand-real-proc))
+                   (y (- alpha (* beta x)))
+                   (t (+ 1.0 (exp y)))
+                   (lhs (+ y (log (/ v (* t t)))))
+                   (rhs (+ k (* n (log L)) (- (log-of-fact n)))))
+              (if (<= lhs rhs)
+                  n
+                  (loop))))))
+    loop))
+
+;; private
+;; log(n!) table for n 1 to 256. Vector, where nth index corresponds to log((n+1)!)
+;; Computed on first invocation of `log-of-fact`
+(define log-fact-table #f)
+
+;; private
+;; computes log-fact-table
+;; log(n!) = log((n-1)!) + log(n)
+(define (make-log-fact-table!)
+  (define table (make-vector 256))
+  (vector-set! table 0 0)
+  (do ((i 1 (+ i 1)))
+      ((> i 255) #t)
+      (vector-set! table i (+ (vector-ref table (- i 1))
+                              (log (+ i 1)))))
+  (set! log-fact-table table))
+
+;; private
+;; returns log(n!)
+;; adapted from https://www.johndcook.com/blog/2010/08/16/how-to-compute-log-factorial/
+(define (log-of-fact n)
+  (when (not log-fact-table)
+    (make-log-fact-table!))
+  (cond
+    ((<= n 1) 0)
+    ((<= n 256) (vector-ref log-fact-table (- n 1)))
+    (else (let ((x (+ n 1)))
+           (+ (* (- x 0.5)
+                 (log x))
+              (- x)
+              (* 0.5
+                 (log (* 2 PI)))
+              (/ 1.0 (* x 12.0)))))))
+
+
+
+(define (gsampling . generators-lst)
+  (let ((gen-vec (list->vector generators-lst))
+        (rand-int-proc (random-source-make-integers (current-random-source))))
+
+    ;remove exhausted generator at index
+    (define (remove-gen index)
+      (define new-vec (make-vector (- (vector-length gen-vec) 1)))
+      ;when removing anything but first, copy all elements before index
+      (when (> index 0)
+        (vector-copy! new-vec 0 gen-vec 0 index))
+      ;when removing anything but last, copy all elements after index
+      (when (< index (- (vector-length gen-vec) 1))
+        (vector-copy! new-vec index gen-vec (+ 1 index)))
+      (set! gen-vec new-vec))
+
+    ;randomly pick generator. If it's exhausted remove it, and pick again
+    ;returns value (or eof, if all generators are exhausted)
+    (define (pick)
+      (let* ((index (rand-int-proc (vector-length gen-vec)))
+             (gen (vector-ref gen-vec index))
+             (value (gen)))
+        (if (eof-object? value)
+            (begin
+              (remove-gen index)
+              (if (= (vector-length gen-vec) 0)
+                  (eof-object)
+                  (pick)))
+            value)))
+
+    (lambda ()
+      (if (= 0 (vector-length gen-vec))
+          (eof-object)
+          (pick)))))
+
+
+;;; Code for binomial random variable generation.
+;;; Written by Brad Lucier, lucier@math.purdue.edu
+
+;;; binomial-geometric is somewhat classical, the
+;;; "First waiting time algorithm" from page 525 of
+;;; Devroye, L. (1986), Non-Uniform Random Variate
+;;; Generation, Springer-Verlag, New York.
+
+;;; binomial-rejection is algorithm BTRS from
+;;; Hormann, W. (1993), The generation of binomial
+;;; random variates, Journal of Statistical Computation
+;;; and Simulation, 46:1-2, 101-110,
+;;; DOI: https://doi.org/10.1080/00949659308811496
+;;; stirling-tail and BTRD (mentioned in the comments)
+;;; are also from that paper.
+
+;;; Another implementation of the same algorithm is at
+;;; https://github.com/tensorflow/tensorflow/blob/master/tensorflow/core/kernels/random_binomial_op.cc
+;;; That implementation pointed out at least two bugs in the
+;;; BTRS paper.
+
+(define (stirling-tail k)
+  ;; Computes
+  ;;
+  ;; \log(k!)-[\log(\sqrt{2\pi})+(k+\frac12)\log(k+1)-(k+1)]
+  ;;
+  (let ((small-k-table
+         ;; Computed using computable reals package
+         ;; Matches values in paper, which are given
+         ;; for 0\leq k < 10
+         '#(.08106146679532726
+            .0413406959554093
+            .02767792568499834
+            .020790672103765093
+            .016644691189821193
+            .013876128823070748
+            .01189670994589177
+            .010411265261972096
+            .009255462182712733
+            .00833056343336287
+            .007573675487951841
+            .00694284010720953
+            .006408994188004207
+            .0059513701127588475
+            .005554733551962801
+            .0052076559196096404
+            .004901395948434738
+            .004629153749334028
+            .004385560249232324
+            .004166319691996922)))
+    (if (< k 20)
+        (vector-ref small-k-table k)
+        ;; the correction term (+ (/ (* 12 (+ k 1))) ...)
+        ;; in Stirling's approximation to log(k!)
+        (let* ((inexact-k+1 (inexact (+ k 1)))
+               (inexact-k+1^2 (square inexact-k+1)))
+          (/ (- #i1/12
+                (/ (- #i1/360
+                      (/ #i1/1260 inexact-k+1^2))
+                   inexact-k+1^2))
+             inexact-k+1)))))
+
+(define (make-binomial-generator n p)
+  (if (not (and (real? p)
+                (<= 0 p 1)
+                (exact-integer? n)
+                (positive? n)))
+       (error "make-binomial-generator: Bad parameters: " n p)
+       (cond ((< 1/2 p)
+              (let ((complement (make-binomial-generator n (- 1 p))))
+                (lambda ()
+                  (- n (complement)))))
+             ((zero? p)
+              (lambda () 0))
+             ((< (* n p) 10)
+              (binomial-geometric n p))
+             (else
+              (binomial-rejection n p)))))
+
+(define (binomial-geometric n p)
+  (let ((geom (make-geometric-generator p)))
+    (lambda ()
+      (let loop ((X -1)
+                 (sum 0))
+        (if (< n sum)
+            X
+            (loop (+ X 1)
+                  (+ sum (geom))))))))
+
+(define (binomial-rejection n p)
+  ;; call when p <= 1/2 and np >= 10
+  ;; Use notation from the paper
+  (let* ((spq
+          (inexact (sqrt (* n p (- 1 p)))))
+         (b
+          (+ 1.15 (* 2.53 spq)))
+         (a
+          (+ -0.0873
+             (* 0.0248 b)
+             (* 0.01 p)))
+         (c
+          (+ (* n p) 0.5))
+         (v_r
+          (- 0.92
+             (/ 4.2 b)))
+         (alpha
+          ;; The formula in BTRS has 1.5 instead of 5.1;
+          ;; The formula for alpha in algorithm BTRD and Table 1
+          ;; and the tensorflow code uses 5.1, so we use 5.1
+          (* (+ 2.83 (/ 5.1 b)) spq))
+         (lpq
+          (log (/ p (- 1 p))))
+         (m
+          (exact (floor (* (+ n 1) p))))
+         (rand-real-proc
+          (random-source-make-reals (current-random-source))))
+    (lambda ()
+      (let loop ()
+        (let* ((u (rand-real-proc))
+               (v (rand-real-proc))
+               (u
+                (- u 0.5))
+               (us
+                (- 0.5 (abs u)))
+               (k
+                (exact
+                 (floor
+                  (+ (* (+ (* 2. (/ a us)) b) u) c)))))
+          (cond ((or (< k 0)
+                     (< n k))
+                 (loop))
+                ((and (<= 0.07 us)
+                      (<= v v_r))
+                 k)
+                (else
+                 (let ((v
+                        ;; The tensorflow code notes that BTRS doesn't have
+                        ;; this logarithm; BTRS is incorrect (see BTRD, step 3.2)
+                        (log (* v (/ alpha
+                                     (+ (/ a (square us)) b))))))
+                   (if (<=  v
+                            (+ (* (+ m 0.5)
+                                  (log (* (/ (+ m 1.)
+                                             (- n m -1.)))))
+                               (* (+ n 1.)
+                                  (log (/ (- n m -1.)
+                                          (- n k -1.))))
+                               (* (+ k 0.5)
+                                  (log (* (/ (- n k -1.)
+                                             (+ k 1.)))))
+                               (* (- k m) lpq)
+                               (- (+ (stirling-tail m)
+                                     (stirling-tail (- n m)))
+                                  (+ (stirling-tail k)
+                                     (stirling-tail (- n k))))))
+                       k
+                       (loop))))))))))