Python sympy.S Examples

def test_multiply(self):
        coords = x, y, z = symbols('x y z', real=True)
        ga, ex, ey, ez = Ga.build('e*x|y|z', g=[1, 1, 1], coords=coords)

        p = Pdop(x)
        assert x * p == Sdop([(x, p)])
        assert ex * p == Sdop([(ex, p)])

        assert p * x == p(x) == S(1)
        assert p * ex == p(ex) == S(0)
        assert type(p(ex)) is Mv

        # These are not defined for consistency with Sdop
        for op in [operator.xor, operator.or_, operator.lt, operator.gt]:
            with pytest.raises(TypeError):
                op(ex, p)
            with pytest.raises(TypeError):
                op(p, ex) 

def _merge_terms(terms1, terms2):
    """ Concatenate and consolidate two sets of already-consolidated terms """
    pdiffs1 = [pdiff for _, pdiff in terms1]
    pdiffs2 = [pdiff for _, pdiff in terms2]

    pdiffs = pdiffs1 + [x for x in pdiffs2 if x not in pdiffs1]
    coefs = len(pdiffs) * [S(0)]

    for coef, pdiff in terms1:
        index = pdiffs.index(pdiff)
        coefs[index] += coef

    for coef, pdiff in terms2:
        index = pdiffs.index(pdiff)
        coefs[index] += coef

    # remove zeros
    return [(coef, pdiff) for coef, pdiff in zip(coefs, pdiffs) if coef != S(0)] 

def test_integrate():
    moments = quadpy.tools.integrate(lambda x: [x ** k for k in range(5)], -1, +1)
    assert (moments == [2, 0, sympy.S(2) / 3, 0, sympy.S(2) / 5]).all()

    moments = quadpy.tools.integrate(
        lambda x: orthopy.line_segment.tree_legendre(x, 4, "monic", symbolic=True),
        -1,
        +1,
    )
    assert (moments == [2, 0, 0, 0, 0]).all()

    # Example from Gautschi's "How to and how not to" article
    moments = quadpy.tools.integrate(
        lambda x: [x ** k * sympy.exp(-(x ** 3) / 3) for k in range(5)], 0, sympy.oo
    )
    S = numpy.vectorize(sympy.S)
    gamma = numpy.vectorize(sympy.gamma)
    n = numpy.arange(5)
    reference = 3 ** (S(n - 2) / 3) * gamma(S(n + 1) / 3)
    assert numpy.all([sympy.simplify(m - r) == 0 for m, r in zip(moments, reference)]) 

def test_boson_states():
    a = BosonOp("a")

    # Fock states
    n = 3
    assert (BosonFockBra(0) * BosonFockKet(1)).doit() == 0
    assert (BosonFockBra(1) * BosonFockKet(1)).doit() == 1
    assert qapply(BosonFockBra(n) * Dagger(a)**n * BosonFockKet(0)) \
        == sqrt(prod(range(1, n+1)))

    # Coherent states
    alpha1, alpha2 = 1.2, 4.3
    assert (BosonCoherentBra(alpha1) * BosonCoherentKet(alpha1)).doit() == 1
    assert (BosonCoherentBra(alpha2) * BosonCoherentKet(alpha2)).doit() == 1
    assert abs((BosonCoherentBra(alpha1) * BosonCoherentKet(alpha2)).doit() -
               exp(-S(1) / 2 * (alpha1 - alpha2) ** 2)) < 1e-12
    assert qapply(a * BosonCoherentKet(alpha1)) == \
        alpha1 * BosonCoherentKet(alpha1) 

def test_eval_trace():
    up = JzKet(S(1)/2, S(1)/2)
    down = JzKet(S(1)/2, -S(1)/2)
    d = Density((up, 0.5), (down, 0.5))

    t = Tr(d)
    assert t.doit() == 1

    #test dummy time dependent states
    class TestTimeDepKet(TimeDepKet):
        def _eval_trace(self, bra, **options):
            return 1

    x, t = symbols('x t')
    k1 = TestTimeDepKet(0, 0.5)
    k2 = TestTimeDepKet(0, 1)
    d = Density([k1, 0.5], [k2, 0.5])
    assert d.doit() == (0.5 * OuterProduct(k1, k1.dual) +
                        0.5 * OuterProduct(k2, k2.dual))

    t = Tr(d)
    assert t.doit() == 1 

def test_booleans():
    assert sympify(sympy.S.true) == true
    assert sympy.S.true == true._sympy_()

    assert sympify(sympy.S.false) == false
    assert sympy.S.false == false._sympy_() 

def test_Reaction_string():
    from sympy import S
    r = Reaction({'A': 1, 'B': 2}, {'C': S(3)/2}, checks=[
        chk for chk in Reaction.default_checks if chk != 'all_integral'])
    assert r.string() == 'A + 2 B -> 3/2 C'
    assert r.string(Reaction_arrow='-->', Reaction_coeff_space='') == 'A + 2B --> 3/2C' 

def test_xk(k):
    n = 10

    moments = quadpy.tools.integrate(
        lambda x: [x ** (i + k) for i in range(2 * n)], -1, +1
    )

    alpha, beta = quadpy.tools.chebyshev(moments)

    assert (alpha == 0).all()
    assert beta[0] == moments[0]
    assert beta[1] == sympy.S(k + 1) / (k + 3)
    assert beta[2] == sympy.S(4) / ((k + 5) * (k + 3))
    quadpy.tools.scheme_from_rc(
        numpy.array([sympy.N(a) for a in alpha], dtype=float),
        numpy.array([sympy.N(b) for b in beta], dtype=float),
        mode="numpy",
    )

    def leg_polys(x):
        return orthopy.line_segment.tree_legendre(x, 19, "monic", symbolic=True)

    moments = quadpy.tools.integrate(
        lambda x: [x ** k * leg_poly for leg_poly in leg_polys(x)], -1, +1
    )

    _, _, a, b = orthopy.line_segment.recurrence_coefficients.legendre(
        2 * n, "monic", symbolic=True
    )

    alpha, beta = quadpy.tools.chebyshev_modified(moments, a, b)

    assert (alpha == 0).all()
    assert beta[0] == moments[0]
    assert beta[1] == sympy.S(k + 1) / (k + 3)
    assert beta[2] == sympy.S(4) / ((k + 5) * (k + 3))
    points, weights = quadpy.tools.scheme_from_rc(
        numpy.array([sympy.N(a) for a in alpha], dtype=float),
        numpy.array([sympy.N(b) for b in beta], dtype=float),
        mode="numpy",
    ) 

def test_gauss_sympy():
    n = 3
    a = 0
    b = 0
    _, _, alpha, beta = orthopy.line_segment.recurrence_coefficients.jacobi(
        n, a, b, "monic", symbolic=True
    )
    points, weights = quadpy.tools.scheme_from_rc(alpha, beta, "sympy")

    assert points == [-sympy.sqrt(sympy.S(3) / 5), 0, +sympy.sqrt(sympy.S(3) / 5)]
    assert weights == [sympy.S(5) / 9, sympy.S(8) / 9, sympy.S(5) / 9] 

def test_chebyshev(dtype):
    alpha = 2

    # Get the moment corresponding to the weight function omega(x) =
    # x^alpha:
    #
    #                                     / 0 if k is odd,
    #    int_{-1}^{+1} |x^alpha| x^k dx ={
    #                                     \ 2/(alpha+k+1) if k is even.
    #
    n = 5

    if dtype == sympy.S:
        moments = [sympy.S(1 + (-1) ** kk) / (kk + alpha + 1) for kk in range(2 * n)]

        alpha, beta = quadpy.tools.chebyshev(moments)

        assert all([a == 0 for a in alpha])
        assert (
            beta
            == [
                sympy.S(2) / 3,
                sympy.S(3) / 5,
                sympy.S(4) / 35,
                sympy.S(25) / 63,
                sympy.S(16) / 99,
            ]
        ).all()
    else:
        assert dtype == numpy.float
        tol = 1.0e-14
        k = numpy.arange(2 * n)
        moments = (1.0 + (-1.0) ** k) / (k + alpha + 1)

        alpha, beta = quadpy.tools.chebyshev(moments)

        assert numpy.all(abs(alpha) < tol)
        assert numpy.all(
            abs(beta - [2.0 / 3.0, 3.0 / 5.0, 4.0 / 35.0, 25.0 / 63.0, 16.0 / 99.0])
            < tol
        ) 

def stroud_1966_b(n):
    s = 1 / sqrt(3)
    data = [(frac(4, 5 * n + 4), z(n))]
    for k in range(1, n + 1):
        r = sqrt(frac(5 * k + 4, 15))
        arr = numpy.full((2 ** (n - k + 1), n), S(0))
        arr[:, k - 1 :] = pm((n - k + 1) * [1])
        arr[:, k - 1] *= r
        arr[:, k:] *= s
        b = frac(5 * 2 ** (k - n + 1), (5 * k - 1) * (5 * k + 4))
        data.append((b, arr))

    points, weights = untangle(data)
    return CnScheme("Stroud 1966b", n, weights, points, 5, _source, 3.393e-14) 

def test_entropy():
    up = JzKet(S(1)/2, S(1)/2)
    down = JzKet(S(1)/2, -S(1)/2)
    d = Density((up, 0.5), (down, 0.5))

    # test for density object
    ent = entropy(d)
    assert entropy(d) == 0.5*log(2)
    assert d.entropy() == 0.5*log(2)

    np = import_module('numpy', min_module_version='1.4.0')
    if np:
        #do this test only if 'numpy' is available on test machine
        np_mat = represent(d, format='numpy')
        ent = entropy(np_mat)
        assert isinstance(np_mat, np.matrixlib.defmatrix.matrix)
        assert ent.real == 0.69314718055994529
        assert ent.imag == 0

    scipy = import_module('scipy', __import__kwargs={'fromlist': ['sparse']})
    if scipy and np:
        #do this test only if numpy and scipy are available
        mat = represent(d, format="scipy.sparse")
        assert isinstance(mat, scipy_sparse_matrix)
        assert ent.real == 0.69314718055994529
        assert ent.imag == 0 

def test_sets():
    x = Integer(2)
    y = Integer(3)
    x1 = sympy.Integer(2)
    y1 = sympy.Integer(3)

    assert Interval(x, y) == Interval(x1, y1)
    assert Interval(x1, y) == Interval(x1, y1)
    assert Interval(x, y)._sympy_() == sympy.Interval(x1, y1)
    assert sympify(sympy.Interval(x1, y1)) == Interval(x, y)

    assert sympify(sympy.S.EmptySet) == EmptySet()
    assert sympy.S.EmptySet == EmptySet()._sympy_()

    assert sympify(sympy.S.UniversalSet) == UniversalSet()
    assert sympy.S.UniversalSet == UniversalSet()._sympy_()

    assert FiniteSet(x, y) == FiniteSet(x1, y1)
    assert FiniteSet(x1, y) == FiniteSet(x1, y1)
    assert FiniteSet(x, y)._sympy_() == sympy.FiniteSet(x1, y1)
    assert sympify(sympy.FiniteSet(x1, y1)) == FiniteSet(x, y)

    x = Interval(1, 2)
    y = Interval(2, 3)
    x1 = sympy.Interval(1, 2)
    y1 = sympy.Interval(2, 3)

    assert Union(x, y) == Union(x1, y1)
    assert Union(x1, y) == Union(x1, y1)
    assert Union(x, y)._sympy_() == sympy.Union(x1, y1)
    assert sympify(sympy.Union(x1, y1)) == Union(x, y)

    assert Complement(x, y) == Complement(x1, y1)
    assert Complement(x1, y) == Complement(x1, y1)
    assert Complement(x, y)._sympy_() == sympy.Complement(x1, y1)
    assert sympify(sympy.Complement(x1, y1)) == Complement(x, y) 

def test_conv9b():
    x = Symbol("x")
    y = Symbol("y")
    assert sympify(sympy.I) == I
    assert sympify(2*sympy.I+3) == 2*I+3
    assert sympify(2*sympy.I/5+sympy.S(3)/5) == 2*I/5+Integer(3)/5
    assert sympify(sympy.Symbol("x")*sympy.I + 3) == x*I+3
    assert sympify(sympy.Symbol("x") + sympy.I*sympy.Symbol("y")) == x+I*y 

def test_conv9():
    x = Symbol("x")
    y = Symbol("y")
    assert (I)._sympy_() == sympy.I
    assert (2*I+3)._sympy_() == 2*sympy.I+3
    assert (2*I/5+Integer(3)/5)._sympy_() == 2*sympy.I/5+sympy.S(3)/5
    assert (x*I+3)._sympy_() == sympy.Symbol("x")*sympy.I + 3
    assert (x+I*y)._sympy_() == sympy.Symbol("x") + sympy.I*sympy.Symbol("y") 

def __eq__(self, A):
        if isinstance(A, Pdop) and self.pdiffs == A.pdiffs:
            return True
        else:
            if len(self.pdiffs) == 0 and A == S(1):
                return True
            return False 

def __call__(self, arg):
        # Ensure that we return the right type even when there are no terms - we
        # do this by adding `0 * d(arg)/d(nonexistant)`, which must be zero, but
        # will be a zero of the right type.
        dummy_var = Dummy('nonexistant')
        terms = self.terms or ((S(0), Pdop(dummy_var)),)
        return sum([coef * pdiff(arg) for coef, pdiff in terms]) 

def __init_from_symbol(self, symbol: Symbol) -> None:
        self.terms = ((S(1), Pdop(symbol)),) 

def _latex(self, print_obj):
        if len(self.terms) == 0:
            return ZERO_STR

        self = self._with_sorted_terms()

        s = ''
        for coef, pdop in self.terms:
            coef_str = print_obj.doprint(coef)
            pd_str = print_obj.doprint(pdop)
            if coef == S(1):
                if pd_str == '':
                    s += '1'
                else:
                    s += pd_str
            elif coef == S(-1):
                if pd_str == '':
                    s += '-1'
                else:
                    s += '-' + pd_str
            else:
                if isinstance(coef, Add):
                    s += r'\left ( ' + coef_str + r'\right ) ' + pd_str
                else:
                    s += coef_str + ' ' + pd_str
            s += ' + '

        s = s.replace('+ -', '- ')
        return s[:-3] 

def _consolidate_terms(terms):
    """
    Remove zero coefs and consolidate coefs with repeated pdiffs.
    """
    new_coefs = []
    new_pdiffs = []
    for coef, pd in terms:
        if coef != S(0):
            if pd in new_pdiffs:
                index = new_pdiffs.index(pd)
                new_coefs[index] += coef
            else:
                new_coefs.append(coef)
                new_pdiffs.append(pd)
    return tuple(zip(new_coefs, new_pdiffs)) 

def test_shorthand(self):
        coords = x, y, z = symbols('x y z', real=True)
        ga, ex, ey, ez = Ga.build('e*x|y|z', g=[1, 1, 1], coords=coords)

        assert Sdop(x) == Sdop([(S(1), Pdop({x: 1}))]) 

def Product_of_Rotors():
    Print_Function()
    (na,nb,nm,alpha,th,th_a,th_b) = symbols('n_a n_b n_m alpha theta theta_a theta_b',\
                                    real = True)
    g = [[na, 0, alpha],[0, nm, 0],[alpha, 0, nb]] #metric tensor
    """
    Values of metric tensor components
    [na,nm,nb] = [+1/-1,+1/-1,+1/-1]  alpha = ea|eb
    """
    (g3d, ea, em, eb) = Ga.build('e_a e_m e_b', g=g)
    print('g =',g3d.g)
    print(r'%n_{a} = \bm{e}_{a}^{2}\;\;n_{b} = \bm{e}_{b}^{2}\;\;n_{m} = \bm{e}_{m}^{2}'+\
        r'\;\;\alpha = \bm{e}_{a}\cdot\bm{e}_{b}')
    (ca,cb,sa,sb) = symbols('c_a c_b s_a s_b',real=True)
    Ra = ca + sa*ea*em  # Rotor for ea^em plane
    Rb = cb + sb*em*eb  # Rotor for em^eb plane
    print(r'%\mbox{Rotor in }\bm{e}_{a}\bm{e}_{m}\mbox{ plane } R_{a} =',Ra)
    print(r'%\mbox{Rotor in }\bm{e}_{m}\bm{e}_{b}\mbox{ plane } R_{b} =',Rb)
    Rab = Ra*Rb  # Compound Rotor
    """
    Show that compound rotor is scalar plus bivector
    """
    print(r'%R_{a}R_{b} = S+\bm{B} =', Rab)
    Rab2 = Rab.get_grade(2)
    print(r'%\bm{B} =',Rab2)
    Rab2sq = Rab2*Rab2  # Square of compound rotor bivector part
    Ssq = (Rab.scalar())**2  # Square of compound rotor scalar part
    Bsq =  Rab2sq.scalar()
    print(r'%S^{2} =',Ssq)
    print(r'%\bm{B}^{2} =',Bsq)
    Dsq = (Ssq-Bsq).expand().simplify()
    print('%S^{2}-B^{2} =', Dsq)
    Dsq = Dsq.subs(nm**2,S(1))  # (e_m)**4 = 1
    print('%S^{2}-B^{2} =', Dsq)
    Cases = [S(-1),S(1)]  # -1/+1 squares for each basis vector
    print(r'#Consider all combinations of $\bm{e}_{a}^{2}$, $\bm{e}_{b}^{2}$'+\
          r' and $\bm{e}_{m}^2$:')
    for Na in Cases:
        for Nb in Cases:
            for Nm in Cases:
                Ba_sq = -Na*Nm
                Bb_sq = -Nb*Nm
                if Ba_sq < 0:
                    Ca_th = cos(th_a)
                    Sa_th = sin(th_a)
                else:
                    Ca_th = cosh(th_a)
                    Sa_th = sinh(th_a)
                if Bb_sq < 0:
                    Cb_th = cos(th_b)
                    Sb_th = sin(th_b)
                else:
                    Cb_th = cosh(th_b)
                    Sb_th = sinh(th_b)
                print(r'%\left [ \bm{e}_{a}^{2},\bm{e}_{b}^{2},\bm{e}_{m}^2\right ] =',\
                      [Na,Nb,Nm])
                Dsq_tmp = Dsq.subs({ca:Ca_th,sa:Sa_th,cb:Cb_th,sb:Sb_th,na:Na,nb:Nb,nm:Nm})
                print(r'%S^{2}-\bm{B}^{2} =',Dsq_tmp,' =',trigsimp(Dsq_tmp))
    print(r'#Thus we have shown that $R_{a}R_{b} = S+\bm{D} = e^{\bm{C}}$ where $\bm{C}$'+\
          r' is a bivector blade.')
    return