(define-module (language python module _random)
  #:use-module (oop pf-objects)
  #:use-module (language python string)
  #:use-module ((language python module python) #:select (int))
  #:use-module (language python def)
  #:use-module (language python try)
  #:use-module (language python exceptions)
  #:export (Random))

(define-syntax-rule (aif it p . l) (let ((it p)) (if p . l)))

(define PI    (* 4 (atan 1)))
(define TWOPI (* 2 PI))
(define LOG4  (log 4.0))
(define _e    (exp 1))

(define NV_MAGICCONST (/ (* 4 (exp -0.5)) (sqrt 2.0)))
(define SG_MAGICCONST (+ 1.0 (log 4.5)))

(define-python-class Random ()
  (define seed
    (lambda (self s)
      (fastset self '_state (seed->random-state (format #f "~a" s)))))

  (define setstate
    (lambda (self s)
      (fastset self '_state s)))

  (define getstate
    (lambda (self)
      (aif it (fastref self '_state)
           it
           (let ((ret (random-state-from-platform)))
             (fastset self '_state ret)
             ret))))

  
  (define getrandbits '(no))
  
  (define random
    (lambda (self)
      (let lp ()
        (set! *random-state* (getstate self))
        (let ((x (random:uniform)))
          (if (= x 1.0)
              (lp)
              (begin
                (fastset self '_state *random-state*)
                x))))))


  (define randrange
    (lambda* (self start #:optional (stop None) (step 1) (_int int))
      (define (fallback)
        ((rawref self '_randrange) self start stop step _int))
      
      (if (number? start)
          (if (eq? stop None)
              (_randbelow self start)
              (if (number? stop)
                  (begin
                    (if (<= stop start)
                        (raise (ValueError "zero range in randrange")))
                    (if (equal? step 1)
                        (+ start (_randbelow self (- stop start)))
                        (if (number? step)
                            (let* ((width (- stop start))
                                   (n (cond
                                       ((= step 0)
                                        (raise
                                         (ValueError
                                          "step of 0 is invalíd in randrange")))
                                       ((> step 0)
                                        (floor-quotient
                                         (+ width  step - 1) step))
                                       (else
                                        (floor-quotient
                                         (+ width  step + 1) step)))))
                              (+ start (* step (_randbelow self n))))
                            (fallback))))
                  (fallback)))
          (fallback))))
  
              
  (define randint
    (lambda (self a b)
      "Return random integer in range [a, b], including both end points.
      "
      (randrange self a b)))
  
  (define _randbelow
    (lambda (self n)
      "Return random integer in range [a, b], including both end points.
      "
      (set! *random-state* (getstate self))
      (let ((x ((@ (guile) random) n)))
        (fastset self '_state *random-state*)
        x)))
      

  ;; -------------------- triangular --------------------

  (define triangular
    (lambda* (self #:optional (low 0.0) (high 1.0) (mode None))
        "Triangular distribution.

        Continuous distribution bounded by given lower and upper limits,
        and having a given mode value in-between.

        http://en.wikipedia.org/wiki/Triangular_distribution

        "
        (let ((u   (random self))
              (c   (if (eq? mode None)
                       0.5
                       (let ((den (- high low)))
                         (if (= den 0)
                             low
                             (/ (- mode low) 1.0 den))))))
              
          (if (> u c)
              (let* ((u    (- 1.0 u))
                     (c    (- 1.0 c))
                     (t    high)
                     (high low)
                     (low  high))
                (+ low (* (- high low) (sqrt (* u c)))))
              (+ low (* (- high low) (sqrt (* u c))))))))

  ;; -------------------- normal distribution --------------------

  (define normalvariate
    (lambda (self mu sigma)
      "Normal distribution.

       mu is the mean, and sigma is the standard deviation.

      "
      ;; mu = mean, sigma = standard deviation

      ;; Uses Kinderman and Monahan method. Reference: Kinderman,
      ;; A.J. and Monahan, J.F., "Computer generation of random
      ;; variables using the ratio of uniform deviates", ACM Trans
      ;; Math Software, 3, (1977), pp257-260.

      (let lp ()
        (let* ((u1 (random self))
               (u2 (- 1.0 (random self)))
               (z  (/ (* NV_MAGICCONST (- u1 0.5)) u2))
               (zz (/ (* z z) 4.0)))
            (if (<= zz (- (log u2)))
                (+ mu (* z sigma))
                (lp))))))
  
  ;; -------------------- lognormal distribution --------------------

  (define lognormvariate
    (lambda (self mu sigma)
      "Log normal distribution.

       If you take the natural logarithm of this distribution, you'll get a
       normal distribution with mean mu and standard deviation sigma.
       mu can have any value, and sigma must be greater than zero.

      "
      (exp (normalvariate self mu sigma))))

  ;;## -------------------- exponential distribution --------------------

  (define expovariate
    (lambda (self lambd)
      "Exponential distribution.

       lambd is 1.0 divided by the desired mean.  It should be
       nonzero.  (The parameter would be called \"lambda\", but that is
       a reserved word in Python.)  Returned values range from 0 to
       positive infinity if lambd is positive, and from negative
       infinity to 0 if lambd is negative.

       "

       ;; lambd: rate lambd = 1/mean
       ;; ('lambda' is a Python reserved word)

       ;; we use 1-random() instead of random() to preclude the
       ;; possibility of taking the log of zero.
       
       (- (/ (log (- 1.0 (random self))) lambd))))

  ;;## -------------------- gamma distribution --------------------

  (define gammavariate
    (lambda (self alpha beta)
        "Gamma distribution.  Not the gamma function!

        Conditions on the parameters are alpha > 0 and beta > 0.

        The probability distribution function is:

                    x ** (alpha - 1) * math.exp(-x / beta)
          pdf(x) =  --------------------------------------
                      math.gamma(alpha) * beta ** alpha

        "

        ;; alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2
        ;; Warning: a few older sources define the gamma distribution in terms
        ;; of alpha > -1.0
        
        (if (or (<= alpha 0.0) (<= beta 0.0))
            (raise (ValueError "gammavariate: alpha and beta must be > 0.0")))

        (cond
         ((> alpha 1.0)
            ;; Uses R.C.H. Cheng, "The generation of Gamma
            ;; variables with non-integral shape parameters",
            ;; Applied Statistics, (1977), 26, No. 1, p71-74

            (let* ((ainv (sqrt (- (* 2.0 alpha) 1.0)))
                   (bbb  (- alpha LOG4))
                   (ccc  (+ alpha ainv)))

              (let lp ()
                (let ((u1 (random self)))
                  (if (or (< u1 1e-7) (> u1 .9999999))
                      (lp)
                      (let* ((u2 (- 1.0 (random self)))
                             (v  (/ (log (/ u1 (- 1.0 u1))) ainv))
                             (x  (* alpha (exp v)))
                             (z  (* u1 u1 u2))
                             (r  (+ bbb (* ccc v) (- x))))
                        (if (or (>= (+ r SG_MAGICCONST (- (* 4.5 z))) 0.0)
                                (>= r (log z)))
                            (* x beta)
                            (lp))))))))

         ((= alpha 1.0)
          ;; expovariate(1)
          (let lp ((u (random self)))
            (if (<= u 1e-7)
                (lp (random self))
                (- (* (log u) beta)))))

         (else
          ;;alpha is between 0 and 1 (exclusive)
          ;; Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle
          
          (let lp ()
            (let* ((u  (random self))
                   (b  (/ (+ _e alpha) _e))
                   (p  (* b u))
                   (x  (if (<= p 1.0)
                           (expt p (/ 1.0 alpha))
                           (- (log (/ (- b p) alpha)))))
                   (u1 (random self)))
              (if (> p 1.0)
                  (if (<= u1 (expt x (- alpha 1.0)))
                      (* x beta)
                      (lp))
                  (if (<= u1 (exp (- x)))
                      (* x beta)
                      (lp)))))))))

  ;; -------------------- beta --------------------
  ;; See
  ;; http://mail.python.org/pipermail/python-bugs-list/2001-January/003752.html
  ;; for Ivan Frohne's insightful analysis of why the original implementation:
  ;;
  ;;    def betavariate(self, alpha, beta):
  ;;        # Discrete Event Simulation in C, pp 87-88.
  ;;
  ;;        y = self.expovariate(alpha)
  ;;        z = self.expovariate(1.0/beta)
  ;;        return z/(y+z)
  ;;
  ;; was dead wrong, and how it probably got that way.

  (define betavariate
    (lambda (self alpha beta)
      "Beta distribution.

       Conditions on the parameters are alpha > 0 and beta > 0.
       Returned values range between 0 and 1.

      "
      
      ;; This version due to Janne Sinkkonen, and matches all the std
      ;; texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution").
      (let ((y (gammavariate self alpha 1.0)))
        (if (= y 0)
            0.0
            (/ y (+ y (gammavariate self beta 1.0)))))))

  ;; -------------------- von Mises distribution --------------------

  (define vonmisesvariate
    (lambda (self mu kappa)
      "Circular data distribution.

      mu is the mean angle, expressed in radians between 0 and 2*pi, and
      kappa is the concentration parameter, which must be greater than or
      equal to zero.  If kappa is equal to zero, this distribution reduces
      to a uniform random angle over the range 0 to 2*pi.

      "
      ;; mu:    mean angle (in radians between 0 and 2*pi)
      ;; kappa: concentration parameter kappa (>= 0)
      ;; if kappa = 0 generate uniform random angle

      ;; Based upon an algorithm published in: Fisher, N.I.,
      ;; "Statistical Analysis of Circular Data", Cambridge
      ;; University Press, 1993.

      ;; Thanks to Magnus Kessler for a correction to the
      ;; implementation of step 4.
      
      (if (<= kappa 1e-6)
          (* TWOPI (random self))
          (let* ((s (/ 0.5 kappa))
                 (r (+ s (sqrt (+ 1.0 (* s s))))))
            (let lp ()
              (let* ((u1 (random self))
                     (z  (cos (* PI u1)))
                     (d  (/ z (+ r z)))
                     (u2 (random self)))
                (if (or (<  u2 (- 1.0 (* d d)))
                        (<= u2 (* (- 1.0 d) (exp d))))
                    (let* ((q  (/ 1.0  r))
                           (f  (/ (+ q  z) (+ 1.0 (* q z))))
                           (u3 (random self)))
                      (if (> u3 0.5)
                          (floor-remainder (+ mu (acos f)) TWOPI)
                          (floor-remainder (- mu (acos f)) TWOPI)))
                    (lp))))))))

  (define uniform
    (lambda (self a b)
      "Get a random number in the range [a, b) or [a, b] depending on rounding."
      (+ a (* (- b a) (random self)))))

  
  ;; -------------------- Pareto --------------------
  (define paretovariate
    (lambda (self alpha)
      "Pareto distribution.  alpha is the shape parameter."
      ;; Jain, pg. 495

      (let ((u (- 1.0 (random self))))
        (/ 1.0 (expt u (/ 1.0 alpha))))))

  ;; -------------------- Weibull --------------------
  (define weibullvariate
    (lambda (self alpha beta)
      "Weibull distribution.

      alpha is the scale parameter and beta is the shape parameter.

      "
      ;; Jain, pg. 499; bug fix courtesy Bill Arms

      (let ((u (- 1.0 (random self))))
        (* alpha (expt (- (log u)) (/ 1.0 beta))))))


  (define gauss
    (lambda (self mu sigma)
      "Gaussian distribution.

        mu is the mean, and sigma is the standard deviation.  This is
        slightly faster than the normalvariate() function.

        Not thread-safe without a lock around calls.

        "

        ;; When x and y are two variables from [0, 1), uniformly
        ;; distributed, then
        ;;
        ;;    cos(2*pi*x)*sqrt(-2*log(1-y))
        ;;   sin(2*pi*x)*sqrt(-2*log(1-y))
        ;;
        ;; are two *independent* variables with normal distribution
        ;; (mu = 0, sigma = 1).
        ;; (Lambert Meertens)
        ;; (corrected version; bug discovered by Mike Miller, fixed by LM)

        ;; Multithreading note: When two threads call this function
        ;; simultaneously, it is possible that they will receive the
        ;; same return value.  The window is very small though.  To
        ;; avoid this, you have to use a lock around all calls.  (I
        ;; didn't want to slow this down in the serial case by using a
        ;; lock here.)

        (let ((z (fastref self 'gauss_next)))
          (fastset self 'gauss_next #f)
          (if (not z)
              (let ((x2pi  (* (random self) TWOPI))
                    (g2rad (sqrt (* -2.0 (log (- 1.0 (random self)))))))
                (set! z (* (cos x2pi) g2rad))
                (fastset self 'gauss_next (* (sin x2pi) g2rad))))

          (+ mu (* z sigma))))))