Python scipy.special 模块，factorial() 实例源码

def FF_Yang_Dou_residual(vbyu, *args):
    """
    The Yang_Dou residual function; to be used by numerical root finder
    """
    (Re, rough) = args

    Rstar = Re / (2 * vbyu * rough)
    theta = np.pi * np.log( Rstar / 1.25) / np.log(100 / 1.25)
    alpha = (1 - np.cos(theta)) / 2
    beta = 1 - (1 - 0.107) * (alpha + theta/np.pi) / 2
    R = Re / (2 * vbyu)

    rt = 1.
    for i in range(1,5):
        rt = rt - 1. / np.e * ( i / factorial(i) * (67.8 / R) ** (2 * i))

    return vbyu - (1 - rt) * R / 4. - rt * (2.5 * np.log(R) - 66.69 * R**-0.72 + 1.8 - (2.5 * np.log(
        (1 + alpha * Rstar / 5) / (1 + alpha * beta * Rstar / 5)) + (5.8 + 1.25) * (alpha * Rstar / (
        5 + alpha * Rstar)) ** 2 + 2.5 * (alpha * Rstar / (5 + alpha * Rstar)) - (5.8 + 1.25)
        * (alpha * beta * Rstar / (5 + alpha * beta * Rstar)) ** 2 - 2.5 * (alpha * beta * Rstar / (
        5 + alpha * beta * Rstar))))

def approx_expm(self,M,exp_t, scaling_terms): 
        #approximate the exp at the beginning to estimate the number of taylor terms and scaling and squaring needed
        U=np.identity(len(M),dtype=M.dtype)
        Mt=np.identity(len(M),dtype=M.dtype)
        factorial=1.0 #for factorials

        for ii in xrange(1,exp_t):
            factorial*=ii
            Mt=np.dot(Mt,M)
            U+=Mt/((2.**float(ii*scaling_terms))*factorial) #scaling by 2**scaling_terms


        for ii in xrange(scaling_terms):
            U=np.dot(U,U) #squaring scaling times

        return U

def approx_exp(self,M,exp_t, scaling_terms): 
        # the scaling and squaring of matrix exponential with taylor expansions
        U=1.0
        Mt=1.0
        factorial=1.0 #for factorials

        for ii in xrange(1,exp_t):
            factorial*=ii
            Mt=M*Mt
            U+=Mt/((2.**float(ii*scaling_terms))*factorial) #scaling by 2**scaling_terms


        for ii in xrange(scaling_terms):
            U=np.dot(U,U) #squaring scaling times

        return U

def sph_harm(m, n, az, el, type='complex'):
    '''Compute sphercial harmonics

    Parameters
    ----------
    m : (int)
        Order of the spherical harmonic. abs(m) <= n

    n : (int)
        Degree of the harmonic, sometimes called l. n >= 0

    az: (float)
        Azimuthal (longitudinal) coordinate [0, 2pi], also called Theta.

    el : (float)
        Elevation (colatitudinal) coordinate [0, pi], also called Phi.

    Returns
    -------
    y_mn : (complex float)
        Complex spherical harmonic of order m and degree n,
        sampled at theta = az, phi = el
    '''
    if type == 'legacy':
        return scy.sph_harm(m, n, az, el)
    elif type == 'real':
        Lnm = scy.lpmv(_np.abs(m), n, _np.cos(el))

        factor_1 = (2 * n + 1) / (4 * _np.pi)
        factor_2 = scy.factorial(n - _np.abs(m)) / scy.factorial(n + abs(m))

        if m != 0:
            factor_1 = 2 * factor_1

        if m < 0:
            return (-1) ** m * _np.sqrt(factor_1 * factor_2) * Lnm * _np.sin(m * az)
        else:
            return (-1) ** m * _np.sqrt(factor_1 * factor_2) * Lnm * _np.cos(m * az)
    else:
        # For the correct Condon–Shortley phase, all m>0 need to be increased by 1
        return (-1) ** _np.float_(m - (m < 0) * (m % 2)) * scy.sph_harm(m, n, az, el)

def test_factorial():
    if scipy_factorial is None:
        return
    for i in range(20):
        assert_equal(factorial(i), scipy_factorial(i))

def factorial(N):
    """Compute the factorial of N.
    If N <= 16, use a fast lookup table; otherwise use scipy.special.factorial
    """
    if N < len(FACTORIALS):
        return FACTORIALS[N]
    elif scipy_special is None:
        raise ValueError("need scipy for computing larger factorials")
    else:
        return int(scipy_special.factorial(N))

def fun(self, x, *args):
        self.nfev += 1

        return (prod(x) - factorial(self.N)) ** 2.0

def weights(dim, degree):
        # 1D sigma-points (x) and weights (w)
        x, w = hermegauss(degree)
        # hermegauss() provides weights that cause posdef errors
        w = factorial(degree) / (degree ** 2 * hermeval(x, [0] * (degree - 1) + [1]) ** 2)
        return np.prod(cartesian([w] * dim), axis=1)

def run(self, t,seed=None,plot=False):
            """
            Gillespie Implementation
            """
            l = self.reactions[:,0]
            r = self.reactions[:,1]
            change=r-l
            time =0
            times =[]
            output =[]
            times.append(time)
            output.append(self.population)
            while time<t:           
                lam = np.prod(comb(self.population,l)*factorial(l),axis=-1)*self.rates
                lam_sum = np.sum(lam)
                dt = np.log(1.0/np.random.uniform(0,1))*(1.0/(lam_sum))
                react = np.where(np.random.multinomial(1,lam/lam_sum)==1)[0][0]
                self.population = self.population + change[react]
                time += dt
                output.append(self.population)
                times.append(time)
            output = np.array(output)
            times = np.array(times)
            if plot:
                for i in xrange(output.shape[1]):
                    label = self.species[i]
                    plt.plot(times, output[:, i], label=label)
                plt.legend(loc='best')
                plt.xlabel('t')
                plt.title('Population Profile')
                plt.grid()
                plt.show()
            return times,output


# def main():
#   reactions = [[[1,0],[0,1]], [[0,1],[1,0]]]
#   rates = [1,1]
#   system = MassActionSystem(reactions,rates)
#   system.set_concentrations([2.,4.])
#   system.run()
#   print system.current_concentrations()
#   print system.dydt()
#   system = StochasticSystem(reactions,rates)
#   population = [10,0]
#   system.set_population(population)
#   t,o = system.run(10,plot=True)
#   return