#	R A N D O M   V A R I A B L E   G E N E R A T O R S
#
#	distributions on the real line:
#	------------------------------
#	       normal (Gaussian)
#	       lognormal
#	       negative exponential
#	       gamma
#
#	distributions on the circle (angles 0 to 2pi)
#	---------------------------------------------
#	       circular uniform
#	       von Mises

# Translated from anonymously contributed C/C++ source.

from whrandom import random, uniform, randint, choice # Also for export!
from math import log, exp, pi, e, sqrt, acos, cos

# Housekeeping function to verify that magic constants have been
# computed correctly

def verify(name, expected):
	computed = eval(name)
	if abs(computed - expected) > 1e-7:
		raise ValueError, \
  'computed value for %s deviates too much (computed %g, expected %g)' % \
  (name, computed, expected)

# -------------------- normal distribution --------------------

NV_MAGICCONST = 4*exp(-0.5)/sqrt(2)
verify('NV_MAGICCONST', 1.71552776992141)
def normalvariate(mu, sigma):
	# 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.

	while 1:
		u1 = random()
		u2 = random()
		z = NV_MAGICCONST*(u1-0.5)/u2
		zz = z*z/4
		if zz <= -log(u2):
			break
	return mu+z*sigma

# -------------------- lognormal distribution --------------------

def lognormvariate(mu, sigma):
	return exp(normalvariate(mu, sigma))

# -------------------- circular uniform --------------------

def cunifvariate(mean, arc):
	# mean: mean angle (in radians between 0 and pi)
	# arc:  range of distribution (in radians between 0 and pi)

	return (mean + arc * (random() - 0.5)) % pi

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

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

	u = random()
	while u <= 1e-7:
		u = random()
	return -log(u)/lambd

# -------------------- von Mises distribution --------------------

TWOPI = 2*pi
verify('TWOPI', 6.28318530718)

def vonmisesvariate(mu, kappa):
	# mu:    mean angle (in radians between 0 and 180 degrees)
	# kappa: concentration parameter kappa (>= 0)
	
	# if kappa = 0 generate uniform random angle
	if kappa <= 1e-6:
		return TWOPI * random()

	a = 1.0 + sqrt(1 + 4 * kappa * kappa)
	b = (a - sqrt(2 * a))/(2 * kappa)
	r = (1 + b * b)/(2 * b)

	while 1:
		u1 = random()

		z = cos(pi * u1)
		f = (1 + r * z)/(r + z)
		c = kappa * (r - f)

		u2 = random()

		if not (u2 >= c * (2.0 - c) and u2 > c * exp(1.0 - c)):
			break

	u3 = random()
	if u3 > 0.5:
		theta = mu + 0.5*acos(f)
	else:
		theta = mu - 0.5*acos(f)

	return theta % pi

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

LOG4 = log(4)
verify('LOG4', 1.38629436111989)

def gammavariate(alpha, beta):
        # beta times standard gamma
	ainv = sqrt(2 * alpha - 1)
	return beta * stdgamma(alpha, ainv, alpha - LOG4, alpha + ainv)

SG_MAGICCONST = 1+log(4.5)
verify('SG_MAGICCONST', 2.50407739677627)

def stdgamma(alpha, ainv, bbb, ccc):
	# ainv = sqrt(2 * alpha - 1)
	# bbb = alpha - log(4)
	# ccc = alpha + ainv

	if alpha <= 0.0:
		raise ValueError, 'stdgamma: alpha must be > 0.0'

	if 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

		while 1:
			u1 = random()
			u2 = random()
			v = log(u1/(1-u1))/ainv
			x = alpha*exp(v)
			z = u1*u1*u2
			r = bbb+ccc*v-x
			if r + SG_MAGICCONST - 4.5*z >= 0 or r >= log(z):
				return x

	elif alpha == 1.0:
		# expovariate(1)
		u = random()
		while u <= 1e-7:
			u = random()
		return -log(u)

	else:	# alpha is between 0 and 1 (exclusive)

		# Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle

		while 1:
			u = random()
			b = (e + alpha)/e
			p = b*u
			if p <= 1.0:
				x = pow(p, 1.0/alpha)
			else:
				# p > 1
				x = -log((b-p)/alpha)
			u1 = random()
			if not (((p <= 1.0) and (u1 > exp(-x))) or
				  ((p > 1)  and  (u1 > pow(x, alpha - 1.0)))):
				break
		return x

# -------------------- test program --------------------

def test():
	print 'TWOPI         =', TWOPI
	print 'LOG4          =', LOG4
	print 'NV_MAGICCONST =', NV_MAGICCONST
	print 'SG_MAGICCONST =', SG_MAGICCONST
	N = 100
	test_generator(N, 'random()')
	test_generator(N, 'normalvariate(0.0, 1.0)')
	test_generator(N, 'lognormvariate(0.0, 1.0)')
	test_generator(N, 'cunifvariate(0.0, 1.0)')
	test_generator(N, 'expovariate(1.0)')
	test_generator(N, 'vonmisesvariate(0.0, 1.0)')
	test_generator(N, 'gammavariate(0.5, 1.0)')
	test_generator(N, 'gammavariate(0.9, 1.0)')
	test_generator(N, 'gammavariate(1.0, 1.0)')
	test_generator(N, 'gammavariate(2.0, 1.0)')
	test_generator(N, 'gammavariate(20.0, 1.0)')
	test_generator(N, 'gammavariate(200.0, 1.0)')

def test_generator(n, funccall):
	import sys
	print '%d calls to %s:' % (n, funccall),
	sys.stdout.flush()
	code = compile(funccall, funccall, 'eval')
	sum = 0.0
	sqsum = 0.0
	for i in range(n):
		x = eval(code)
		sum = sum + x
		sqsum = sqsum + x*x
	avg = sum/n
	stddev = sqrt(sqsum/n - avg*avg)
	print 'avg %g, stddev %g' % (avg, stddev)

if __name__ == '__main__':
	test()