def b_mat(ind_mat):
    r""" Calculates the B coefficients from [1]_ Eq. 27.

    Parameters
    ----------
    index_matrix : array, shape (N,3)
        ordering of the basis in x, y, z

    Returns
    -------
    B : array, shape (N,)
        B coefficients for the basis

    References
    ----------
    .. [1] Ozarslan E. et. al, "Mean apparent propagator (MAP) MRI: A novel
    diffusion imaging method for mapping tissue microstructure",
    NeuroImage, 2013.
    """

    B = np.zeros(ind_mat.shape[0])
    for i in range(ind_mat.shape[0]):
        n1, n2, n3 = ind_mat[i]
        K = int(not(n1 % 2) and not(n2 % 2) and not(n3 % 2))
        B[i] = K * np.sqrt(factorial(n1) * factorial(n2) * factorial(n3)
                           ) / (factorial2(n1) * factorial2(n2) * factorial2(n3))

    return B

def radon_n_bspline_general( y , theta , m , n ):
    exponent = 2*m - n + 1
    y_plus = positive_power( exponent )
    

    ##  Calculate  m+1 fold derivative (analytically)
    deriv_positive_power = scalar_product( ( misc.factorial( exponent ) / \
                   misc.factorial( exponent - (m +1) ) ) , positive_power( m - n ) )

    if theta == 0.0 or theta == np.pi/2.0 or theta == np.pi:
        y_plus = deriv_positive_power
        
    
    ##  Consider special case of theta = 0 , pi
    if np.abs( theta - np.pi/2.0 ) > eps:
        for i in range( m + 1 ):
            y_plus = finite_difference( y_plus, np.cos(theta) )
    

    ##  Consider special case of theta = pi/2
    if np.abs( theta ) > eps  and np.abs( theta - np.pi ) > eps:
        for i in range( m + 1 ):
            y_plus = finite_difference( y_plus , np.sin(theta) )
    
    return y_plus(y)/float( misc.factorial( exponent ) )  

    def B(l,ll,r,rr,i,l1,l2,Ra,Rb,Rp,g1,l3,l4,Rc,Rd,Rq,g2):
        """
        Expansion coefficient B.

        Source:
            Handbook of Computational Chemistry
            David Cook
            Oxford University Press
            1998
        """

        b = 1
        b *= (-1)**(l) * theta(l,l1,l2,Rp-Ra,Rp-Rb,r,g1)
        b *= theta(ll,l3,l4,Rq-Rc,Rq-Rd,rr,g2)
        b *= (-1)**i * (2*delta)**(2*(r+rr))
        b *= misc.factorial(l + ll - 2*r - 2*rr,exact=True)
        b *= delta**i * (Rp-Rq)**(l+ll-2*(r+rr+i))

        tmp = 1
        tmp *= (4*delta)**(l+ll) * misc.factorial(i,exact=True)
        tmp *= misc.factorial(l+ll-2*(r+rr+i),exact=True)

        b /= tmp

        return b

    def normal_reference_constant(self):
        """
        Constant used for silverman normal reference asymtotic bandwidth
        calculation.

        C  = 2((pi^(1/2)*(nu!)^3 R(k))/(2nu(2nu)!kap_nu(k)^2))^(1/(2nu+1))
        nu = kernel order
        kap_nu = nu'th moment of kernel
        R = kernel roughness (square of L^2 norm)

        Note: L2Norm property returns square of norm.
        """
        nu = self._order

        if not nu == 2:
            msg = "Only implemented for second order kernels"
            raise NotImplementedError(msg)

        if self._normal_reference_constant is None:
            C = np.pi ** (0.5) * factorial(nu) ** 3 * self.L2Norm
            C /= 2 * nu * factorial(2 * nu) * self.moments(nu) ** 2
            C = 2 * C ** (1.0 / (2 * nu + 1))
            self._normal_reference_constant = C

        return self._normal_reference_constant

    def wdm_linearity(self, domain):
        # Calculate dispersive terms:
        if self.beta is None:
            self.factor = (0.0, 0.0)
        else:
            if self.centre_omega is None:
                self.Domega = (domain.omega - domain.centre_omega,
                               domain.omega - domain.centre_omega)
            else:
                self.Domega = (domain.omega - self.centre_omega[0],
                               domain.omega - self.centre_omega[1])

            terms = [0.0, 0.0]
            for n, beta in enumerate(self.beta[0]):
                terms[0] += beta * np.power(self.Domega[0], n) / factorial(n)
            for n, beta in enumerate(self.beta[1]):
                terms[1] += beta * np.power(self.Domega[1], n) / factorial(n)
            self.factor = (1j * fftshift(terms[0]), 1j * fftshift(terms[1]))

        # Include attenuation terms if available:
        if self.alpha is None:
            return self.factor
        else:
            self.factor[0] -= 0.5 * self.alpha[0]
            self.factor[1] -= 0.5 * self.alpha[1]
            return factor

def arcsin(n,x):
    #construct a backwards order of your n's
    integers = np.arange(n,-1,-1, dtype=np.float64)
    coeff = factorial(integers*2)/((2*integers+1)*(factorial(integers)**2)*(4**integers))
    allcoeff = np.zeros(2*n+2)
    allcoeff[::2]=coeff[:]
    return np.polyval(allcoeff,x)

def get_n_combos(n,k):
    """
    Returns the number of combinations of n choose k.
    Uses Stirling's approximation when numbers are very large.
    """

    result = None
    try:
        fac_n = factorial(n)
        fac_k = factorial(k)
        fac_nk = factorial(n-k)
        summ = fac_n + fac_k + fac_nk
        if summ == np.inf or summ != summ:
            raise ValueError("Values too large. Using Stirling's approximation")
        result = fac_n/fac_k
        result /= fac_nk
    except: # Catch all large number exceptions.  
        pass
    if not result or result == np.inf or result != result: # No result yet.  
        if n == k or k == 0:
            result = 1
        else:
            x = stirling(n) - stirling(k) - stirling(n-k) # Use Stirling's approx.
            result = np.exp(x) 
        
    return result

    def update(self, rho):
        """
        Calculate the probability function for the given state of an harmonic
        oscillator (as density matrix)
        """

        if isket(rho):
            rho = ket2dm(rho)

        self.data = np.zeros(len(self.xvecs[0]), dtype=complex)
        M, N = rho.shape

        for m in range(M):
            k_m = pow(self.omega / pi, 0.25) / \
                sqrt(2 ** m * factorial(m)) * \
                exp(-self.xvecs[0] ** 2 / 2.0) * \
                np.polyval(hermite(m), self.xvecs[0])

            for n in range(N):
                k_n = pow(self.omega / pi, 0.25) / \
                    sqrt(2 ** n * factorial(n)) * \
                    exp(-self.xvecs[0] ** 2 / 2.0) * \
                    np.polyval(hermite(n), self.xvecs[0])

                self.data += np.conjugate(k_n) * k_m * rho.data[m, n]

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
    n = int(n)
    sch = np.array([1.0])
    sch2 = np.array([(-1.0) ** m * np.sqrt(
        (2.0 * factorial(n - m)) / factorial(n + m)) for m in range(1, n + 1)])
    sch = np.append(sch, sch2)
    if isinstance(x, float) or len(x) == 1:
        leg = lpmv(np.arange(0, n + 1), n, x)
        return np.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 = np.array([sch * leg]).T
            else:
                out = np.append(out, np.array([sch * leg]).T, axis=1)
    return out

def series_problem_a():
	c = np.arange(70, -1, -1) # original values for n
	c = factorial(2*c) / ((2*c+1) * factorial(c)**2 * 4**c) #series coeff's
	p = np.zeros(2*c.size) # make space for skipped zero-terms
	p[::2] = c # set nonzero polynomial terms to the series coeff's
	P = np.poly1d(p) # make a polynomial out of it
	return 6 * P(.5) #return pi (since pi/6 = arcsin(1/2))

def cumulant_from_moments(momt, n):
    """Compute n-th cumulant given moments.

    Parameters
    ----------
    momt: array_like
        `momt[j]` contains `(j+1)`-th moment.
        These can be raw moments around zero, or central moments
        (in which case, `momt[0]` == 0).
    n: integer
        which cumulant to calculate (must be >1)

    Returns
    -------
    kappa: float
        n-th cumulant.

    """
    if n < 1:
        raise ValueError("Expected a positive integer. Got %s instead." % n)
    if len(momt) < n:
        raise ValueError("%s-th cumulant requires %s moments, "
                         "only got %s." % (n, n, len(momt)))
    kappa = 0.
    for p in _faa_di_bruno_partitions(n):
        r = sum(k for (m, k) in p)
        term = (-1)**(r - 1) * factorial(r - 1)
        for (m, k) in p:
            term *= np.power(momt[m - 1] / factorial(m), k) / factorial(k)
        kappa += term
    kappa *= factorial(n)
    return kappa

def GL_mode(u,p):
	r = len(u[:,0])
	z = len(u[0,:])
	L = zeros((r,z))
	for i in range(p+1):
		L += float(factorial(p))*(-u)**i/( float(factorial(p - i))*(float(factorial(i)))**2. )
	return L

    def update_rho(self, rho):
        """
        calculate probability distribution for quadrature measurement
        outcomes given a two-mode density matrix
        """

        X1, X2 = np.meshgrid(self.xvecs[0], self.xvecs[1])

        p = np.zeros((len(self.xvecs[0]), len(self.xvecs[1])), dtype=complex)
        N = rho.dims[0][0]

        M1 = np.zeros((N, N, len(self.xvecs[0]), len(self.xvecs[1])), dtype=complex)
        M2 = np.zeros((N, N, len(self.xvecs[0]), len(self.xvecs[1])), dtype=complex)

        for m in range(N):
            for n in range(N):
                M1[m,n] = exp(-1j * self.theta1 * (m - n)) / \
                    sqrt(pi * 2 ** (m + n) * factorial(n) * factorial(m)) * \
                    exp(-X1 ** 2) * np.polyval(hermite(m), X1) * np.polyval(hermite(n), X1)
                M2[m,n] = exp(-1j * self.theta2 * (m - n)) / \
                    sqrt(pi * 2 ** (m + n) * factorial(n) * factorial(m)) * \
                    exp(-X2 ** 2) * np.polyval(hermite(m), X2) * np.polyval(hermite(n), X2)

        for n1 in range(N):
            for n2 in range(N):
                i = state_number_index([N, N], [n1, n2])
                for p1 in range(N):
                    for p2 in range(N):
                        j = state_number_index([N, N], [p1, p2])
                        p += M1[n1, p1] * M2[n2, p2] * rho.data[i, j]

        self.data = p

def beta(p,a):
    """Returns the pth coefficient beta.

    inputs:
    p -- (int) pth coefficient
    a -- (float) uncorrected fractional distance of MC

    output:

    beta_p -- pth beta coefficient for given vel. or accel.
    """

    beta_p = a ** (p+1) / factorial(p+1)

    if p > 0:

        for q in range(p):

            if a >= 0:

                beta_p -= beta(q,a) / factorial(p + 1 - q)

            else: # a < 0

                beta_p -= (-1) ** (p+q) * beta(q,a) / factorial(p + 1 - q)

    return beta_p

 def eval_expx(lamb,i,j):
     '''get the i,j element of exp(lamb(a+a^dag))'''
     res=0.
     r=min(i,j)+1
     for m in xrange(r):
         res+=lamb**(i+j-2*m)*sqrt(factorial(j)*factorial(i))/factorial(m)/factorial(i-m)/factorial(j-m)
     return exp(lamb**2/2)*res #it is a pending issue whether '-' sign make sense here.

def shore_odf_matrix(radial_order, mu, smoment, vertices):
    """
    Eq. 33 (choose ux=uy=uz)
    """

    ind_mat = shore_index_matrix(radial_order)

    n_vert = vertices.shape[0]

    n_elem = ind_mat.shape[0]

    odf_mat = np.zeros((n_vert, n_elem))

    rho = mu

    for i in range(n_vert):

        vx, vy, vz = vertices[i]

        for j in range(n_elem):

            n1, n2, n3 = ind_mat[j]
            f = np.sqrt(factorial(n1) * factorial(n2) * factorial(n3))

            k = mu ** (smoment) / np.sqrt(2 ** (2 - smoment) * np.pi ** 3)

            odf_mat[i, j] = k * f * _odf_cfunc(n1, n2, n3, vx, vy, vz, smoment)

    return odf_mat

    def time(self, t, s=1.0):
        """
        Complex Paul wavelet, centred at zero.

        Parameters
        ----------
        t : float
            Time. If `s` is not specified, i.e. set to 1, this can be
            used as the non-dimensional time t/s.
        s : float
            Scaling factor. Default is 1.

        Returns
        -------
        out : complex
            Value of the Paul wavelet at the given time

        The Paul wavelet is defined (in time) as::

            (2 ** m * i ** m * m!) / (pi * (2 * m)!) \
                    * (1 - i * t / s) ** -(m + 1)

        """
        m = self.m
        x = t / s

        const = (2 ** m * 1j ** m * factorial(m)) \
            / (np.pi * factorial(2 * m)) ** .5
        functional_form = (1 - 1j * x) ** -(m + 1)

        output = const * functional_form

        return output

def beta_m(a, B, N):
    A = np.zeros([N,N])
    alpha = np.zeros([N,1])
    for i in range(N):
        alpha[i,0] = a ** (i+1) / factorial(i+1)
        for j in range(i+1):
            A[i,j] = B[i-j]/factorial(i-j)
    return A.dot(alpha)

 def integrand(zbar, T=T, data=data):
     fromb_times = data[1][data[2]=="B"]
     #total = len(fromb_times)
     numbefore = np.sum(fromb_times<=zbar)
     numafter = np.sum(fromb_times>zbar)
     pbefore = zbar**numbefore*np.exp(-zbar)/float(factorial([numbefore])[0])
     pafter = (1.4*(T-zbar))**numafter*np.exp(-1.4*(T-zbar))/float(factorial([numafter])[0])
     return pbefore*pafter*.2*np.exp(-.2*zbar)

def likelihood(theta, n, x):
    """
    :param theta: probability of infection
    :param n: n testing
    :param x: x positive
    :return likelihood probability
    """
    return (factorial(n) / (factorial(n-x) * factorial(x))) * (1 - theta) ** (n - x) * (theta ** x)