###########################################################################
#                                                                         #
# C O P Y R I G H T   N O T I C E                                         #
#  Copyright (c) 2003-10 by:                                              #
#    * California Institute of Technology                                 #
#                                                                         #
#    All Rights Reserved.                                                 #
#                                                                         #
# Permission is hereby granted, free of charge, to any person             #
# obtaining a copy of this software and associated documentation files    #
# (the "Software"), to deal in the Software without restriction,          #
# including without limitation the rights to use, copy, modify, merge,    #
# publish, distribute, sublicense, and/or sell copies of the Software,    #
# and to permit persons to whom the Software is furnished to do so,       #
# subject to the following conditions:                                    #
#                                                                         #
# The above copyright notice and this permission notice shall be          #
# included in all copies or substantial portions of the Software.         #
#                                                                         #
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,         #
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF      #
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND                   #
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS     #
# BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN      #
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN       #
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE        #
# SOFTWARE.                                                               #
###########################################################################
#
# special functions that are helpful for our statistics
# gammaln, gamma, factorial, factln, bico, and beta are inspired by:
#           Numerical Recipes in C
#           by Press, Flannery, Teukolsky, Vetterling
#
from math import log, exp, floor

def gammaln(Z):
    cof = [ 76.18009173, -86.50532033, 24.01409822, -1.231739516, 0.120858003e-2, -0.536382e-5]

    x = float(Z) - 1.0
    temp = x + 5.5
    ser = 1.0
    j = 0

    temp -= (x + 0.5) * log(temp)
    for j in range(6):
        x += 1.0
        ser += cof[j] / x

    return -temp + log(2.50662827465 * ser)


def gammaln2(Z):
    cof = [676.5203681218835, -1259.139216722289, 771.3234287757674, -176.6150291498386, 
           12.50734324009056, -0.1385710331296526, 0.9934937113930748e-5, 0.1659470187408462e-6]

    ser = 0.9999999999995183
    x = float(Z) - 1.0

    for j in range(7):
        x += 1.0
        ser += cof[j] / x

    result = log(ser) - 5.58106146679532777 - float(Z) + (float(Z) - 0.5) * log(float(Z) + 6.5)

    return result


def gamma(Z):
    return round(exp(gammaln(Z)), 8)


def gamma2(Z):
    return round(exp(gammaln2(Z)), 8)


def factorial(N):
    ntop = 4
    arr = [1.0, 1.0, 2.0, 6.0, 24.0]

    if N < 0:
        print "error: factorial only works for non-negative integers"
        return -1
    elif N > 32:
        return gamma(N + 1)
    else:
        while ntop < N:
            j = ntop
            ntop += 1
            arr.append(arr[j] * ntop)
        return arr[N]


def factln(N):
    if N < 0:
        print "error: factorial only works for non-negative integers"
        return -1
    if N <= 1:
        return 0.0
    else:
        return gammaln2(N + 1.0)


def bicoln(n, k):
    return factln(n) - factln(k) - factln(n - k)


def bico(n, k):
    return floor(0.5 + exp(bicoln(n, k)))


def beta(z, w):
    return exp(gammaln(z) + gammaln(w) - gammaln(z + w))


def hyperGeometric(N, r, n, x):
    """ hyperGeometric(N, r, n, x):
                N   is population size
                r   is the number in the population identified as a success
                n   is the sample size
                x   is the number in the sample identified as a success
    """
    N = float(N)
    r = float(r)
    n = float(n)
    x = float(x)
    return exp(bicoln(r, x) + bicoln(N - r, n - x) - bicoln(N, n))