def alpha(k,m,n):
     tau = t_intervals[n-1]
     i = np.arange(m+1)[:,np.newaxis]
     the_sum = np.sum((1j*k*omega)**(-m+i-1) * \
                          tau[:,np.newaxis]**i / factorial(i), axis=-2)
     integral = -factorial(m) * np.exp(-1j * k * omega * tau) * the_sum
     return integral[:,1] - integral[:,0]

def Combinations(values, k):
    """This function outputs all possible combinations of k elements from the column vector values"""
    n = len(values)
    try:
        values = sp.row_stack(values)
    except:
        raise ValueError, "I need a 2d column array"

    if k > n:
        raise ValueError, "k must be <= %d" % n
    elif k<=0 or k%1 != 0:
        raise ValueError, "k must be > 0"

    #out = sp.array([],ndmin=2)
    if k == 1:
        return values
    else:
        #This loop iterates through all the elements of the values that have at least
        #k elements.  For each element it then calls Combinations(values[i+1:], k-1) which
        #returns combinations of size k-1 for the elements succeeding the current element
        #We do not want to get repeats of combinations
        #print "for i in range(%d)" % (n-(k-1))
        for i in range(n-(k-1)):
            #Calculate the number of possible combinations (to allow proper concatenation
            #in the recursive call
            numCombs = sp.factorial(n-i)/(sp.factorial(k-1)*sp.factorial(n-i-(k-1)))
            combs = Combinations(values[i:], k-1)
            ones = values[i]*sp.ones((numCombs,1))
            #print "ones: %s \t\t combs: %s" % (str(ones.shape), str(combs.shape))
            print combs

def prehamiltonian( genlags, EC, EJ, EL ):
    ### NB: omits a factor of EJ (on top of the flux dependency and the 
    ### diagonal terms) for use as a derivative later
    size = len(genlags)
    hbar_w0 = sqrt( 8. * EL * EC )
    phi0 = ( 8. * EC / EL ) ** .25 
    arg = phi0**2/2
    
    genlags = [[ f(arg) for f in row ] for row in genlags]
    ret = [ range(size) for i in range(size) ] #values set below

    for row in range(size):
        for col in range(size):
            #the nonzero cosine elements
            if (col-row)%2==0:
                n = min(row,col)
                m = abs(col-row)/2 # because of Hermitianness
                ret[row][col] = -(-2)**-m \
                    * sqrt(factorial(n)/factorial(n+2*m)) \
                    * phi0**(2*m) * exp(phi0**2/-4) \
                    * genlags[n][2*m]
            #the nonzero sine elements
            else:
                ### IS THIS PART RIGHT?
                n = min(row,col)
                m = (abs(col-row)-1)/2
                ret[row][col] = -(-2)**(-m) * 2**-.5 \
                    * sqrt(factorial(n)/factorial(n+2*m+1)) \
                    * phi0**(2*m+1) * exp(phi0**2/-4) \
                    * genlags[n][2*m+1] ## Check overall signs
    return array(ret)

def Combinations(values, k):
    """This function outputs all the possible combinations of k elements from the vector values"""
    
    if int(k) < 0:
        raise ValueError("k must a positive integer")
    
    #Make input vectors column vectors
    if values.shape == (1,values.size):
        values = sp.atleast2d(values).T.copy()
    
    out = sp.array([]).reshape(0,1)
    n = max(values.shape)
    if k == 1:
        out = values
    else:
        #the following loop interates through all the elements of the vector values that have at least k elements after them.  For each element it then calls Combinations(values[i+1:], k-1) which returns combinations of size k-1 for the elements succeeding the current element.  This is so that we do not get repeats of combinations
        #nck = sp.misc.comb(n,k, exact=True)
        #out = sp.zeros((nck, k))
        for i in range(n-(k-1)):
            #Calculate the number of possible combinations (to allow proper concatenation in the recursive call
            nCombs = int(sp.factorial(n-i)/(sp.factorial(k-1)*sp.factorial(n-i-(k-1))))
            #This is the recursive call
            
            print Combinations(values[i+1:], k-1).reshape((-1,1))
            out = sp.r_[out, sp.c_[values[i]*sp.ones((nCombs,1)), Combinations(values[i+1:], k-1).reshape(-1,1)]]
            
    return out

def _sch_lpmv(n,x):
	'''
	Outputs array of Schmidt Seminormalized Associated Legendre Functions S_{n}^{m} for m<=n.
	
	Parameters
	----------
	n : int 
	    Degree of polynomial.
	
	x : float 
	    Point at which to evaluate
	
	Returns
	-------
	array of values for Legendre functions. 
	
	'''
	from scipy.special import lpmv
	sch=array([1.0])
	sch2=array([(-1.0)**m*sqrt((2.0*factorial(n-m))/factorial(n+m)) for m in range(1,n+1)])
	sch=append(sch,sch2)
	if isinstance(x,float) or len(x)==1:
		leg=lpmv(arange(0,n+1),n,x)
		return array([sch*leg]).T
	else:
		for j in range(0,len(x)):
			leg=lpmv(range(0,n+1),n,x[j])
			if j==0:
				out=array([sch*leg]).T
			else:
				out=append(out,array([sch*leg]).T,axis=1)
	return out

 def delta(a,b,c):
     """ Calculate delta """
     fac = zeros(4,long)
     fac[0] = factorial(a+b-c)
     fac[1] = factorial(a-b+c)
     fac[2] = factorial(-a+b+c)
     fac[3] = factorial(a+b+c+1)
     return sqrt(prod(fac[0:3])/fac[3]);

def basis2d(n0,n1,beta=[1.,1.]):
    """2d dimensionless Cartesian basis function"""
    b=hermite2d(n0,n1)
    b[0]*=((2**n0)*(n.pi**(.5))*scipy.factorial(n0))**(-.5)
    exp0=lambda x: beta[0] * b[0](x) * n.exp(-.5*(x**2))
    b[1]*=((2**n1)*(n.pi**(.5))*scipy.factorial(n1))**(-.5)
    exp1=lambda x: beta[1] * b[1](x) * n.exp(-.5*(x**2))
    return [exp0,exp1]

def bernstein(n):
    """ Bernstein-Polynome (n!/(i!(n-i)!)* t^i (1-t)^(n-i)"""
    t = linspace(0,1,20)

    bmat = zeros((20,n+1))
    for i in range(0,n+1):
        bmat[:,i] = sp.factorial(n)/(sp.factorial(i)*sp.factorial(n-i))*t**i*(1-t)**(n-i)
    return bmat

def plot_messung(n, t, b_len, y_lim, lens, rate):
    mu, sigma = np.mean(n), np.std(n)
    print "Mittelwert: %.3f" % mu
    print "Varianz: %.3f Poisson: %.3f" % (sigma ** 2, mu)
    print "Standardabweichung: %.3f Poisson: %.3f" % (sigma, np.sqrt(mu))
    print "Standardabweichung des Mittelwerts: %.3f Poisson: %.3f" % (
        sigma / np.sqrt(len(n)),
        np.sqrt(mu) / np.sqrt(len(n)),
    )
    print "Gesamtanzahl der Ereignisse:", np.sum(n)
    print "mittlere Zaehlrate: %.3f 1/s" % (np.mean(n) / (rate))
    print "Standardabweichung der Zaehlrate: %.3f 1/s" % (np.std(n) / rate)
    print "Standardabweichung der mittleren Zaehlrate: %.3f 1/s" % (np.std(n) / rate / np.sqrt(len(n)))
    print "Schiefe: %.3f Poisson: %.3f" % (np.mean((n - n.mean()) ** 3) / sigma ** 3, 1 / np.sqrt(mu))
    print "Kurtosis: %.3f Poisson: %.3f" % (np.mean((n - n.mean()) ** 4) / sigma ** 4 - 3, 1 / (mu))
    print ""
    fig = plt.figure(figsize=(16, 12))
    ax = fig.add_subplot(111)

    # the histogram of the data
    n, bins, patches = ax.hist(n, lens, normed=1, facecolor="yellow", alpha=0.75)

    bincenters = 0.5 * (bins[1:] + bins[:-1])
    b = np.linspace(0, b_len, 1000)
    b2 = np.arange(0, b_len, 1) + 0.5

    # add a 'best fit' line for the normal PDF
    # y = mlab.normpdf( b, mu, sigma)
    # y2 =mlab.normpdf( b2, mu, sigma)
    # l = ax.plot(b, y, 'b--', linewidth=1)
    poisson = lambda k: 1.0 / (sc.factorial(k) * np.exp(mu)) * mu ** k
    poisson2 = lambda k: 1.0 / (sc.factorial(k) * np.exp(mu + sigma)) * (mu + sigma) ** k
    poisson3 = lambda k: 1.0 / (sc.factorial(k) * np.exp(mu - sigma)) * (mu - sigma) ** k
    normal_k = lambda k: 1.0 / np.sqrt(2 * np.pi * mu) * np.exp(-(k - mu) ** 2 / (2 * mu))
    normal_k2 = (
        lambda k: 1.0 / np.sqrt(2 * np.pi * (mu + sigma)) * np.exp(-(k - (mu + sigma)) ** 2 / (2 * (mu + sigma)))
    )
    normal_k3 = (
        lambda k: 1.0 / np.sqrt(2 * np.pi * (mu - sigma)) * np.exp(-(k - (mu - sigma)) ** 2 / (2 * (mu - sigma)))
    )

    nk = ax.plot(b, normal_k(b), "g--", linewidth=1)
    nk = ax.plot(b, normal_k2(b), "g--", linewidth=1)
    nk = ax.plot(b, normal_k3(b), "g--", linewidth=1)
    p = ax.plot(b, poisson(b), "r--", linewidth=1)
    p = ax.plot(b, poisson2(b), "r--", linewidth=1)
    p = ax.plot(b, poisson3(b), "r--", linewidth=1)
    l = ax.scatter(b2, normal_k(b2), marker="x", c="b")
    p = ax.scatter(b2, poisson(b2), marker="x", c="b")

    ax.set_xlabel("Anzahl der Impulse")
    ax.set_ylabel("Wahrscheinlichkeit")
    plt.xlim(0, b_len)
    plt.ylim(0, y_lim)
    ax.grid(True)

    plt.show()

def dimBasis2d(n0,n1,beta=[1.,1.],phs=[1.,1.]):
    """2d dimensional Cartesian basis function of characteristic size beta
    phs: additional phase factor, used in the Fourier Transform"""
    b=hermite2d(n0,n1)
    b[0]*=(beta[0]**(-.5))*(((2**n0)*(n.pi**(.5))*scipy.factorial(n0))**(-.5))
    exp0=lambda x: b[0](x/beta[0]) * n.exp(-.5*((x/beta[0])**2)) * phs[0]
    b[1]*=(beta[1]**(-.5))*(((2**n1)*(n.pi**(.5))*scipy.factorial(n1))**(-.5))
    exp1=lambda x: b[1](x/beta[1]) * n.exp(-.5*((x/beta[1])**2)) * phs[1]
    return [exp0,exp1]

def C(n, r):
    if n - r > r:
        num = np.prod(np.arange(n - r + 1, n+1))
        den = factorial(r)
        return num / den
    else:
        num = np.prod(np.arange(r + 1, n+1))
        den = factorial(n - r)
        return num / den

def plot_pdf(order, N, iterations):
	order_stats = []
	for it in range(iterations):
		numbers = [np.random.uniform(0,1) for i in range(N)]
		numbers.sort()
		order_stats.append(numbers[order-1])
	plt.figure()
	n, bins, patches = plt.hist(order_stats, iterations/20, normed=1, facecolor='green')
	y = lambda x: int(sp.factorial(N))/(int(sp.factorial(N-order))*int(sp.factorial(order-1))) * x**(order-1) * (1-x)**(N-order)
	plt.plot(bins, y(bins), 'r--', linewidth=3)

def factorialQuotient(numerator, denominator):
    """
    result=numerator!/(denominator!)
    """
    diff = numerator - denominator
    if diff > 0:
        result = factorial(diff)
    else:
        result = 1.0 / factorial(-diff)
    return result

def Wigner6j(j1,j2,j3,J1,J2,J3):
#======================================================================
# Calculating the Wigner6j-Symbols using the Racah-Formula                
# Author: Ulrich Krohn                                            
# Date: 13th November 2009
#                                                                         
# Based upon Wigner3j.m from David Terr, Raytheon                         
# Reference: http://mathworld.wolfram.com/Wigner6j-Symbol.html            
#
# Usage: 
# from wigner import Wigner6j
# WignerReturn = Wigner6j(j1,j2,j3,J1,J2,J3)
#
#  / j1 j2 j3 \
# <            >  
#  \ J1 J2 J3 /
#
#======================================================================

    # Check that the js and Js are only integer or half integer
    if ( ( 2*j1 != round(2*j1) ) | ( 2*j2 != round(2*j2) ) | ( 2*j2 != round(2*j2) ) | ( 2*J1 != round(2*J1) ) | ( 2*J2 != round(2*J2) ) | ( 2*J3 != round(2*J3) ) ):
        print 'All arguments must be integers or half-integers.'
        return -1
    
# Check if the 4 triads ( (j1 j2 j3), (j1 J2 J3), (J1 j2 J3), (J1 J2 j3) ) satisfy the triangular inequalities
    if ( ( abs(j1-j2) > j3 ) | ( j1+j2 < j3 ) | ( abs(j1-J2) > J3 ) | ( j1+J2 < J3 ) | ( abs(J1-j2) > J3 ) | ( J1+j2 < J3 ) | ( abs(J1-J2) > j3 ) | ( J1+J2 < j3 ) ):
        print '6j-Symbol is not triangular!'
        return 0
    
    # Check if the sum of the elements of each traid is an integer
    if ( ( 2*(j1+j2+j3) != round(2*(j1+j2+j3)) ) | ( 2*(j1+J2+J3) != round(2*(j1+J2+J3)) ) | ( 2*(J1+j2+J3) != round(2*(J1+j2+J3)) ) | ( 2*(J1+J2+j3) != round(2*(J1+J2+j3)) ) ):
        print '6j-Symbol is not triangular!'
        return 0
    
    # Arguments for the factorials
    t1 = j1+j2+j3
    t2 = j1+J2+J3
    t3 = J1+j2+J3
    t4 = J1+J2+j3
    t5 = j1+j2+J1+J2
    t6 = j2+j3+J2+J3
    t7 = j1+j3+J1+J3

    # Finding summation borders
    tmin = max(0, max(t1, max(t2, max(t3,t4))))
    tmax = min(t5, min(t6,t7))
    tvec = arange(tmin,tmax+1,1)
        
    # Calculation the sum part of the 6j-Symbol
    WignerReturn = 0
    for t in tvec:
        WignerReturn += (-1)**t*factorial(t+1)/( factorial(t-t1)*factorial(t-t2)*factorial(t-t3)*factorial(t-t4)*factorial(t5-t)*factorial(t6-t)*factorial(t7-t) )

    # Calculation of the 6j-Symbol
    return WignerReturn*sqrt( TriaCoeff(j1,j2,j3)*TriaCoeff(j1,J2,J3)*TriaCoeff(J1,j2,J3)*TriaCoeff(J1,J2,j3) )

def multivariate_polya_vectorized(x,alpha):
    """Multivariate Pólya PDF. Vectorized implementation.
    """
    x = np.atleast_1d(x)
    alpha = np.atleast_1d(alpha)
    assert(x.size==alpha.size)
    N = x.sum()
    A = alpha.sum()
    likelihood = factorial(N) / factorial(x).prod() * gamma(A) / gamma(N + A)
    likelihood *= (gamma(x + alpha) / gamma(alpha)).prod()
    return likelihood

def plot_messung3():
    H = np.array([6.45 * 10 ** 2, 3.45 * 10 ** 2, 9.5 * 10, 14, 2])
    mu = np.sum(H * np.arange(0, 5, 1)) / np.sum(H)
    n = np.arange(0, 5, 1)
    sigma = np.sqrt((np.sum(H * (np.arange(0, 5, 1) - mu) ** 2) / (np.sum(H) - 1)))
    print "Mittelwert: %.3f" % mu
    print "Varianz: %.3f Poisson: %.3f" % (sigma ** 2, mu)
    print "Standardabweichung: %.3f Poisson: %.3f" % (sigma, np.sqrt(mu))
    print "Standardabweichung des Mittelwerts: %.3f Poisson: %.3f" % (
        sigma / np.sqrt(np.sum(H)),
        np.sqrt(mu) / np.sqrt(np.sum(H)),
    )
    print "Gesamtanzahl der Ereignisse:", np.sum(H * n)
    print "mittlere Zaehlrate: %.3f 1/s" % (mu / (0.5))
    print "Standardabweichung der Zaehlrate: %.3f 1/s" % (sigma / 0.5)
    print "Standardabweichung der mittleren Zaehlrate: %.3f 1/s" % (sigma / 0.5 / np.sqrt(np.sum(H)))
    print "Schiefe: %.3f Poisson: %.3f" % ((np.sum(H * (n - mu) ** 3) / np.sum(H)) / sigma ** 3, 1 / np.sqrt(mu))
    print "Kurtosis: %.3f Poisson: %.3f" % ((np.sum(H * (n - mu) ** 4) / np.sum(H)) / sigma ** 4 - 3, 1 / (mu))

    print ""

    fig = plt.figure(figsize=(16, 12))
    plt.scatter(np.arange(0, 4 + 1, 1), H / 1101, marker="^", s=80)
    poisson = lambda k: 1.0 / (sc.factorial(k) * np.exp(mu)) * mu ** k
    poisson_min = lambda k: 1.0 / (sc.factorial(k) * np.exp(mu - sigma)) * (mu - sigma) ** k
    poisson_max = lambda k: 1.0 / (sc.factorial(k) * np.exp(mu + sigma)) * (mu + sigma) ** k
    normal_k = lambda k: 1.0 / np.sqrt(2 * np.pi * mu) * np.exp(-(k - mu) ** 2 / (2 * mu))
    normal_k_min = (
        lambda k: 1.0 / np.sqrt(2 * np.pi * (mu - sigma)) * np.exp(-(k - mu - sigma) ** 2 / (2 * mu - 2 * sigma))
    )
    normal_k_max = (
        lambda k: 1.0 / np.sqrt(2 * np.pi * (mu + sigma)) * np.exp(-(k - mu + sigma) ** 2 / (2 * mu + 2 * sigma))
    )
    b = np.linspace(0, 5, 100)
    b2 = np.arange(0, 5, 1)
    # plt.plot(b, normal_k(b), 'g--', linewidth=1)
    plt.plot(b, normal_k_min(b), "g--", linewidth=1)
    # plt.plot(b, normal_k_max(b), 'g--', linewidth=1)
    # plt.plot(b, poisson(b), 'r--', linewidth=1)
    plt.plot(b, poisson_min(b), "r--", linewidth=1)
    # plt.plot(b, poisson_max(b), 'r--', linewidth=1)
    plt.scatter(b2, normal_k(b2), marker="x", c="b")
    plt.scatter(b2, poisson(b2), marker="x", c="b")
    plt.title("Messung 3")
    plt.xlabel("Anzahl der Pulse")
    plt.ylabel("Wahrscheinlichkeit")
    plt.grid(True)
    plt.xlim(0, 5)
    plt.ylim(0, 0.7)
    plt.show()

def multivariate_polya(x, alpha):
    """Multivariate Pólya PDF. Basic implementation.
    """
    x = np.atleast_1d(x).flatten()
    alpha = np.atleast_1d(alpha).flatten()
    assert(x.size==alpha.size)
    N = x.sum()
    A = alpha.sum()
    likelihood = factorial(N) * gamma(A) / gamma(N + A)
    # likelihood = gamma(A) / gamma(N + A)
    for i in range(len(x)):
        likelihood /= factorial(x[i])
        likelihood *= gamma(x[i] + alpha[i]) / gamma(alpha[i])
    return likelihood

def _wigner_laguerre(rho, xvec, yvec, g, parallel):
    """
    Using Laguerre polynomials from scipy to evaluate the Wigner function for
    the density matrices :math:`|m><n|`, :math:`W_{mn}`. The total Wigner
    function is calculated as :math:`W = \sum_{mn} \\rho_{mn} W_{mn}`.
    """

    M = prod(rho.shape[0])
    X, Y = meshgrid(xvec, yvec)
    A = 0.5 * g * (X + 1.0j * Y)
    W = zeros(np.shape(A))

    # compute wigner functions for density matrices |m><n| and
    # weight by all the elements in the density matrix
    B = 4 * abs(A) ** 2
    if sp.isspmatrix_csr(rho.data):
        # for compress sparse row matrices
        if parallel:
            iterator = (
                (m, rho, A, B) for m in range(len(rho.data.indptr) - 1))
            W1_out = parfor(_par_wig_eval, iterator)
            W += sum(W1_out)
        else:
            for m in range(len(rho.data.indptr) - 1):
                for jj in range(rho.data.indptr[m], rho.data.indptr[m + 1]):
                    n = rho.data.indices[jj]

                    if m == n:
                        W += real(rho[m, m] * (-1) ** m * genlaguerre(m, 0)(B))

                    elif n > m:
                        W += 2.0 * real(rho[m, n] * (-1) ** m *
                                        (2 * A) ** (n - m) *
                                        sqrt(factorial(m) / factorial(n)) *
                                        genlaguerre(m, n - m)(B))
    else:
        # for dense density matrices
        B = 4 * abs(A) ** 2
        for m in range(M):
            if abs(rho[m, m]) > 0.0:
                W += real(rho[m, m] * (-1) ** m * genlaguerre(m, 0)(B))
            for n in range(m + 1, M):
                if abs(rho[m, n]) > 0.0:
                    W += 2.0 * real(rho[m, n] * (-1) ** m *
                                    (2 * A) ** (n - m) *
                                    sqrt(factorial(m) / factorial(n)) *
                                    genlaguerre(m, n - m)(B))

    return 0.5 * W * g ** 2 * np.exp(-B / 2) / pi

def Binomial(x, n, p):
    """
    Binomial Log-Likelihood 

    :param x: data
    :param n:
    :param p:

    >>> Binomial([2,3],6,0.3)
    -2.81280615454
    """
    x = array(x)
    like = sum(x * log(p) + (n - x) * log(1. - p) + log(scipy.factorial(n)) - log(scipy.factorial(x)) - log(
        scipy.factorial(n - x)))
    return like

def Negbin(x, r, p):
    """
    Negative Binomial Log-Likelihood 

    :param x: data
    :param r:
    :param p:

    >>> Negbin([2,3],6,0.3)
    -9.16117424315
    """
    x = array(x)
    like = sum(r * log(p) + x * log(1 - p) + log(scipy.factorial(x + r - 1)) - log(scipy.factorial(x)) - log(
        scipy.factorial(r - 1)))
    return like