Python special.factorial方法代码示例

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import factorial [as 别名]
def _K_T_xs(self, temperature, volume, params): # K_T_xs, eq. 3.20 of de Koker thesis
        f = self._finite_strain(temperature, volume, params)
        theta = self._theta(temperature, volume, params)
        K_ToverV=0.
        for i in range(len(params['a'])):
            ifact=factorial(i, exact=False)
            for j in range(len(params['a'][0])):
                if i > 0:
                    jfact=factorial(j, exact=False)
                    prefactor = float(i) * params['a'][i][j] \
                        * np.power(theta, float(j)) / ifact / jfact 
                    K_ToverV += prefactor*self._d2fdV2(temperature, volume, params) \
                        * np.power(f, float(i-1))
                if i > 1:
                    dfdV = self._dfdV(temperature, volume, params)
                    K_ToverV += prefactor * dfdV * dfdV \
                        * float(i-1) * np.power(f, float(i-2))
        return volume*K_ToverV 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import factorial [as 别名]
def _alphaK_T_xs(self, temperature, volume, params): # eq. 3.21 of de Koker thesis
        f = self._finite_strain(temperature, volume, params)
        theta = self._theta(temperature, volume, params)
        sum_factors = 0.
        for i in range(len(params['a'])):
            ifact=factorial(i, exact=False)
            if i > 0:
                for j in range(len(params['a'][0])):
                    if j > 0:
                        jfact=factorial(j, exact=False)
                        sum_factors += float(i)*float(j)*params['a'][i][j] \
                            * np.power(f, float(i-1)) * np.power(theta, float(j-1)) \
                            / ifact / jfact
                            
        return -self._dfdV(temperature, volume, params) \
            * self._dthetadT(temperature, volume, params) \
            * sum_factors 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import factorial [as 别名]
def _get_M(self, delta_t):
        n = len(self.a)
        A = np.zeros(n)
        for i, ai in enumerate(self.a):
            Ae = [self.a[i] * (-1)**j * factorial(i) / factorial(j) / factorial(
                i-j) / ((delta_t)**i) for j in range(i + 1)]  # Elementary A to assemblate in A
            for j, aej in enumerate(Ae):
                A[j] += aej

        n = len(self.b)
        B = np.zeros(n)
        for i, ai in enumerate(self.b):
            Be = [self.b[i] * (-1)**j * factorial(i) / factorial(j) / factorial(
                i-j) / ((delta_t)**i) for j in range(i + 1)]  # Elementary B to assemblate in B
            for j, bej in enumerate(Be):
                B[j] += bej

        Mo = [-x/B[0] for x in B[1:][::-1]]
        Mi = [x/B[0] for x in A[::-1]]
        return (Mi, Mo) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import factorial [as 别名]
def _falling_factorial(x, n):
    r"""
    Return the factorial of `x` to the `n` falling.

    This is defined as:

    .. math::   x^\underline n = (x)_n = x (x-1) \cdots (x-n+1)

    This can more efficiently calculate ratios of factorials, since:

    n!/m! == falling_factorial(n, n-m)

    where n >= m

    skipping the factors that cancel out

    the usual factorial n! == ff(n, n)
    """
    val = 1
    for k in range(x - n + 1, x + 1):
        val *= k
    return val 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import factorial [as 别名]
def __init__(self, xi, yi, axis=0):
        _Interpolator1DWithDerivatives.__init__(self, xi, yi, axis)

        self.xi = np.asarray(xi)
        self.yi = self._reshape_yi(yi)
        self.n, self.r = self.yi.shape

        c = np.zeros((self.n+1, self.r), dtype=self.dtype)
        c[0] = self.yi[0]
        Vk = np.zeros((self.n, self.r), dtype=self.dtype)
        for k in xrange(1,self.n):
            s = 0
            while s <= k and xi[k-s] == xi[k]:
                s += 1
            s -= 1
            Vk[0] = self.yi[k]/float(factorial(s))
            for i in xrange(k-s):
                if xi[i] == xi[k]:
                    raise ValueError("Elements if `xi` can't be equal.")
                if s == 0:
                    Vk[i+1] = (c[i]-Vk[i])/(xi[i]-xi[k])
                else:
                    Vk[i+1] = (Vk[i+1]-Vk[i])/(xi[i]-xi[k])
            c[k] = Vk[k-s]
        self.c = c 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import factorial [as 别名]
def imp_lsc(self, gamma, sigma, w, dz):
        """
        gamma - energy
        sigma - transverse RMS size of the beam
        w - omega = 2*pi*f
        """
        eps = 1e-16
        ass = 40.0

        alpha = w * sigma / (gamma * speed_of_light)
        alpha2 = alpha * alpha

        inda = np.where(alpha2 > ass)[0]
        ind = np.where((alpha2 <= ass) & (alpha2 >= eps))[0]

        T = np.zeros(w.shape)
        T[ind] = np.exp(alpha2[ind]) *exp1(alpha2[ind])

        x= alpha2[inda]
        k = 0
        for i in range(10):
            k += (-1) ** i * factorial(i) / (x ** (i + 1))
        T[inda] = k
        Z = 1j * Z0 / (4 * pi * speed_of_light*gamma**2) * w * T * dz
        return Z # --> Omm/m 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import factorial [as 别名]
def H(n, x):
    """Explicit expression for Hermite polynomials."""
    prefactor = factorial(n)
    terms = tf.reduce_sum(
        tf.stack(
            [
                _numer_safe_power(-1, m)
                / (factorial(m) * factorial(n - 2 * m))
                * _numer_safe_power(2 * x, n - 2 * m)
                for m in range(int(np.floor(n / 2)) + 1)
            ],
            axis=0,
        ),
        axis=0,
    )
    return prefactor * terms 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import factorial [as 别名]
def orbit_cardinality(orbit: list, modes: int) -> int:
    """Gives the number of samples belonging to the input orbit.

    For example, there are three possible samples in the orbit [2, 1, 1] with three modes: [1, 1,
    2], [1, 2, 1], and [2, 1, 1]. With four modes, there are 12 samples in total.

    **Example usage:**

    >>> orbit_cardinality([2, 1, 1], 4)
    12

    Args:
        orbit (list[int]): orbit; we count how many samples are contained in it
        modes (int): number of modes in the samples

    Returns:
        int: number of samples in the orbit
    """
    sample = orbit + [0] * (modes - len(orbit))
    counts = list(Counter(sample).values())
    return int(factorial(modes) / np.prod(factorial(counts))) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import factorial [as 别名]
def test_displaced_squeezed_state_fock(self, r_d, phi_d, r_s, phi_s, hbar, cutoff, tol):
        """test displaced squeezed state returns correct Fock basis state vector"""
        state = utils.displaced_squeezed_state(r_d, phi_d, r_s, phi_s, basis="fock", fock_dim=cutoff, hbar=hbar)
        a = r_d * np.exp(1j * phi_d)

        if r_s == 0:
            pytest.skip("test only non-zero squeezing")

        n = np.arange(cutoff)
        gamma = a * np.cosh(r_s) + np.conj(a) * np.exp(1j * phi_s) * np.sinh(r_s)
        coeff = np.diag(
            (0.5 * np.exp(1j * phi_s) * np.tanh(r_s)) ** (n / 2) / np.sqrt(fac(n) * np.cosh(r_s))
        )

        expected = H(gamma / np.sqrt(np.exp(1j * phi_s) * np.sinh(2 * r_s)), coeff)
        expected *= np.exp(
            -0.5 * np.abs(a) ** 2 - 0.5 * np.conj(a) ** 2 * np.exp(1j * phi_s) * np.tanh(r_s)
        )

        assert np.allclose(state, expected, atol=tol, rtol=0) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import factorial [as 别名]
def test_odd_cat_state(self, tol):
        """test correct even cat state returned"""
        a = 0.212
        cutoff = 10
        p = 1

        state = utils.cat_state(a, p, fock_dim=cutoff)

        n = np.arange(cutoff)
        expected = np.exp(-0.5 * np.abs(a) ** 2) * a ** n / np.sqrt(fac(n)) - np.exp(
            -0.5 * np.abs(-a) ** 2
        ) * (-a) ** n / np.sqrt(fac(n))
        expected /= np.linalg.norm(expected)

        assert np.allclose(state, expected, atol=tol, rtol=0)


# ===================================================================================
# Random matrix tests
# =================================================================================== 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import factorial [as 别名]
def test_loss_channel_on_coherent_states(
        self, setup_backend, T, mag_alpha, phase_alpha, cutoff, tol
    ):
        """Tests various loss channels on coherent states (result should be coherent state with amplitude weighted by sqrt(T)."""

        rootT_alpha = np.sqrt(T) * mag_alpha * np.exp(1j * phase_alpha)

        backend = setup_backend(1)

        backend.prepare_coherent_state(mag_alpha, phase_alpha, 0)
        backend.loss(T, 0)
        s = backend.state()
        if s.is_pure:
            numer_state = s.ket()
        else:
            numer_state = s.dm()
        n = np.arange(cutoff)
        ref_state = (
            np.exp(-0.5 * np.abs(rootT_alpha) ** 2)
            * rootT_alpha ** n
            / np.sqrt(factorial(n))
        )
        ref_state = np.outer(ref_state, np.conj(ref_state))
        assert np.allclose(numer_state, ref_state, atol=tol, rtol=0.0) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import factorial [as 别名]
def test_nongaussian(self, a, phi, setup_backend, cutoff, tol):
        """Tests that probabilities of particular Fock states |n> are
        correct for a nongaussian state."""
        backend = setup_backend(2)

        alpha = a * np.exp(1j * phi)
        n = np.arange(cutoff)
        ref_state = np.exp(-0.5 * np.abs(alpha) ** 2) * alpha ** n / np.sqrt(fac(n))
        ref_probs = np.abs(ref_state) ** 2

        backend.prepare_coherent_state(np.abs(alpha), np.angle(alpha), 0)
        backend.prepare_fock_state(cutoff // 2, 1)
        state = backend.state()

        for n in range(cutoff):
            prob_n = state.fock_prob([n, cutoff // 2])
            assert np.allclose(prob_n, ref_probs[n], atol=tol, rtol=0) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import factorial [as 别名]
def test_pure(self, a, phi, setup_backend, cutoff, batch_size, tol):
        """Tests that the numeric probabilities in the full Fock basis are
        correct for a one-mode pure state."""
        backend = setup_backend(1)

        alpha = a * np.exp(1j * phi)
        n = np.arange(cutoff)
        ref_state = np.exp(-0.5 * np.abs(alpha) ** 2) * alpha ** n / np.sqrt(fac(n))
        ref_probs = np.abs(ref_state) ** 2

        backend.prepare_coherent_state(np.abs(alpha), np.angle(alpha), 0)
        state = backend.state()

        probs = state.all_fock_probs(cutoff=cutoff)
        if isinstance(probs, tf.Tensor):
            probs = probs.numpy()
        probs = probs.flatten()

        if batch_size is not None:
            ref_probs = np.tile(ref_probs, batch_size)

        assert np.allclose(probs, ref_probs, atol=tol, rtol=0) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import factorial [as 别名]
def test_reduced_state_fock_probs(self, cutoff, setup_backend, batch_size, tol):
        """Test backend calculates correct fock prob of reduced coherent state"""
        backend = setup_backend(2)

        backend.prepare_coherent_state(np.abs(a), np.angle(a), 0)
        backend.prepare_squeezed_state(r, phi, 1)

        state = backend.state(modes=[0])
        probs = np.array([state.fock_prob([i]) for i in range(cutoff)]).T

        n = np.arange(cutoff)
        ref_state = np.exp(-0.5 * np.abs(a) ** 2) * a ** n / np.sqrt(fac(n))
        ref_probs = np.abs(ref_state) ** 2

        if batch_size is not None:
            ref_probs = np.tile(ref_probs, batch_size)

        assert np.allclose(probs.flatten(), ref_probs.flatten(), atol=tol, rtol=0) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import factorial [as 别名]
def test_rdm(self, setup_backend, tol, cutoff, batch_size):
        """Test reduced density matrix of a coherent state is as expected
        This is supported by all backends, since it returns
        the reduced density matrix of a single mode."""
        backend = setup_backend(2)
        backend.prepare_coherent_state(np.abs(a), np.angle(a), 0)
        backend.prepare_coherent_state(0.1, 0, 1)

        state = backend.state()
        rdm = state.reduced_dm(0, cutoff=cutoff)

        n = np.arange(cutoff)
        ket = np.exp(-0.5 * np.abs(a) ** 2) * a ** n / np.sqrt(fac(n))
        rdm_exact = np.outer(ket, ket.conj())

        if batch_size is not None:
            np.tile(rdm_exact, [batch_size, 1])

        assert np.allclose(rdm, rdm_exact, atol=tol, rtol=0)