Python sympy 模块,symbols() 实例源码

我们从Python开源项目中,提取了以下50个代码示例,用于说明如何使用sympy.symbols()

项目:pyrsss    作者:butala    | 项目源码 | 文件源码
def chapman_vec(z_vec, Nm_vec, Hm_vec, H_O_vec):
    """
    Vectorized implementation of the Chapman function evaluation
    routine :func:`chapman`. The input arguments must be sequences
    with the same length and the output is an :class:`NP.ndarray` with
    that length.
    """
    try:
        chapman_vec._chapman_sym_f
    except AttributeError:
        sym_vars = SYM.symbols('z Nm Hm H_O')
        chapman_vec._chapman_sym_f = SYM.lambdify(sym_vars,
                                                  chapman_sym(*sym_vars),
                                                  modules='numexpr')
    return chapman_vec._chapman_sym_f(z_vec,
                                      Nm_vec,
                                      Hm_vec,
                                      H_O_vec)
项目:muesr    作者:bonfus    | 项目源码 | 文件源码
def isSymbolic(self):
        """
        True if Fourier components are defined as sympy symbols, 
        False otherwise.

        ALWAYS False for MM objects.
        """        
        return False


#def set_k(mm, kval):
#    '''
#    Sets propagation vector of a MagneticModel object.
#    '''
#    if (type(kval) is np.ndarray):
#        mm.k = kval
#        return True
#        
#    elif (type(kval) is list):
#        mm.k = np.array(kval)
#        return True
#        
#    else:
#        return False
#
项目:muesr    作者:bonfus    | 项目源码 | 文件源码
def set_params(self, values):
            """
            Coverts the symbolic values into numerical values.
            The order of the parameters must be the same as the one given
            in the class instantiaion.

            :param list values: a list containing the values for the sympy symbols.            
            :return: None
            :rtype: None
            :raises: RuntimeError           
            """



            if self._symFClambda != None:
                self.fc_set( np.array(self._symFClambda(*values), \
                                      dtype=np.complex), \
                             self._inputType)
            else:
                raise RuntimeError("Symbolic FC not defined!")
项目:pymoskito    作者:cklb    | 项目源码 | 文件源码
def __init__(self, settings):
        Trajectory.__init__(self, settings)

        # setup symbolic expressions
        tau, k = sp.symbols('tau, k')

        gamma = self._settings["differential_order"] + 1
        alpha = sp.factorial(2 * gamma + 1)

        f = sp.binomial(gamma, k) * (-1) ** k * tau ** (gamma + k + 1) / (gamma + k + 1)
        phi = alpha / sp.factorial(gamma) ** 2 * sp.summation(f, (k, 0, gamma))

        # differentiate phi(tau), index in list corresponds to order
        dphi_sym = [phi]  # init with phi(tau)
        for order in range(self._settings["differential_order"]):
            dphi_sym.append(dphi_sym[-1].diff(tau))

        # lambdify
        self.dphi_num = []
        for der in dphi_sym:
            self.dphi_num.append(sp.lambdify(tau, der, 'numpy'))
项目:texta    作者:texta-tk    | 项目源码 | 文件源码
def _build_logical_expression(self, grammar, terminal_component_names):
        terminal_component_symbols = eval("symbols('%s')"%(' '.join(terminal_component_names)))
        if isinstance(terminal_component_symbols, Symbol):
            terminal_component_symbols = [terminal_component_symbols]
        name_to_symbol = {terminal_component_names[i]:symbol for i, symbol in enumerate(terminal_component_symbols)}
        terminal_component_names = set(terminal_component_names)

        op_to_symbolic_operation = {'not':operator.invert, 'concat':operator.and_, 'gap':operator.and_, 'union':operator.or_, 'intersect':operator.and_}

        def logical_expression_builder(component):
            if component['id'] in terminal_component_names:
                return name_to_symbol[component['id']]
            else:
                children = component['components']
                return reduce(op_to_symbolic_operation[component['operation']],[logical_expression_builder(child) for child in children])

        return simplify(logical_expression_builder(grammar))
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_Tuple():
    t = (1, 2, 3, 4)
    st = Tuple(*t)
    assert set(sympify(t)) == set(st)
    assert len(t) == len(st)
    assert set(sympify(t[:2])) == set(st[:2])
    assert isinstance(st[:], Tuple)
    assert st == Tuple(1, 2, 3, 4)
    assert st.func(*st.args) == st
    p, q, r, s = symbols('p q r s')
    t2 = (p, q, r, s)
    st2 = Tuple(*t2)
    assert st2.atoms() == set(t2)
    assert st == st2.subs({p: 1, q: 2, r: 3, s: 4})
    # issue 5505
    assert all(isinstance(arg, Basic) for arg in st.args)
    assert Tuple(p, 1).subs(p, 0) == Tuple(0, 1)
    assert Tuple(p, Tuple(p, 1)).subs(p, 0) == Tuple(0, Tuple(0, 1))

    assert Tuple(t2) == Tuple(Tuple(*t2))
    assert Tuple.fromiter(t2) == Tuple(*t2)
    assert Tuple.fromiter(x for x in range(4)) == Tuple(0, 1, 2, 3)
    assert st2.fromiter(st2.args) == st2
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_twave():
    A1, phi1, A2, phi2, f = symbols('A1, phi1, A2, phi2, f')
    c = Symbol('c')  # Speed of light in vacuum
    n = Symbol('n')  # Refractive index
    t = Symbol('t')  # Time
    w1 = TWave(A1, f, phi1)
    w2 = TWave(A2, f, phi2)
    assert w1.amplitude == A1
    assert w1.frequency == f
    assert w1.phase == phi1
    assert w1.wavelength == c/(f*n)
    assert w1.time_period == 1/f
    w3 = w1 + w2
    assert w3.amplitude == sqrt(A1**2 + 2*A1*A2*cos(phi1 - phi2) + A2**2)
    assert w3.frequency == f
    assert w3.wavelength == c/(f*n)
    assert w3.time_period == 1/f
    assert w3.angular_velocity == 2*pi*f
    assert w3.equation() == sqrt(A1**2 + 2*A1*A2*cos(phi1 - phi2) + A2**2)*cos(2*pi*f*t + phi1 + phi2)
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_dagger():
    x = symbols('x')
    expr = Dagger(x)
    assert str(expr) == 'Dagger(x)'
    ascii_str = \
"""\
 +\n\
x \
"""
    ucode_str = \
u("""\
\n\
x \
""")
    assert pretty(expr) == ascii_str
    assert upretty(expr) == ucode_str
    assert latex(expr) == r'x^{\dag}'
    sT(expr, "Dagger(Symbol('x'))")
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_scalars():
    x = symbols('x', complex=True)
    assert Dagger(x) == conjugate(x)
    assert Dagger(I*x) == -I*conjugate(x)

    i = symbols('i', real=True)
    assert Dagger(i) == i

    p = symbols('p')
    assert isinstance(Dagger(p), adjoint)

    i = Integer(3)
    assert Dagger(i) == i

    A = symbols('A', commutative=False)
    assert Dagger(A).is_commutative is False
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_UGate():
    a, b, c, d = symbols('a,b,c,d')
    uMat = Matrix([[a, b], [c, d]])

    # Test basic case where gate exists in 1-qubit space
    u1 = UGate((0,), uMat)
    assert represent(u1, nqubits=1) == uMat
    assert qapply(u1*Qubit('0')) == a*Qubit('0') + c*Qubit('1')
    assert qapply(u1*Qubit('1')) == b*Qubit('0') + d*Qubit('1')

    # Test case where gate exists in a larger space
    u2 = UGate((1,), uMat)
    u2Rep = represent(u2, nqubits=2)
    for i in range(4):
        assert u2Rep*qubit_to_matrix(IntQubit(i, 2)) == \
            qubit_to_matrix(qapply(u2*IntQubit(i, 2)))
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_UGate_CGate_combo():
    a, b, c, d = symbols('a,b,c,d')
    uMat = Matrix([[a, b], [c, d]])
    cMat = Matrix([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, a, b], [0, 0, c, d]])

    # Test basic case where gate exists in 1-qubit space.
    u1 = UGate((0,), uMat)
    cu1 = CGate(1, u1)
    assert represent(cu1, nqubits=2) == cMat
    assert qapply(cu1*Qubit('10')) == a*Qubit('10') + c*Qubit('11')
    assert qapply(cu1*Qubit('11')) == b*Qubit('10') + d*Qubit('11')
    assert qapply(cu1*Qubit('01')) == Qubit('01')
    assert qapply(cu1*Qubit('00')) == Qubit('00')

    # Test case where gate exists in a larger space.
    u2 = UGate((1,), uMat)
    u2Rep = represent(u2, nqubits=2)
    for i in range(4):
        assert u2Rep*qubit_to_matrix(IntQubit(i, 2)) == \
            qubit_to_matrix(qapply(u2*IntQubit(i, 2)))
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_cg_simp_sum():
    x, a, b, c, cp, alpha, beta, gamma, gammap = symbols(
        'x a b c cp alpha beta gamma gammap')
    # Varshalovich 8.7.1 Eq 1
    assert cg_simp(x * Sum(CG(a, alpha, b, 0, a, alpha), (alpha, -a, a)
                   )) == x*(2*a + 1)*KroneckerDelta(b, 0)
    assert cg_simp(x * Sum(CG(a, alpha, b, 0, a, alpha), (alpha, -a, a)) + CG(1, 0, 1, 0, 1, 0)) == x*(2*a + 1)*KroneckerDelta(b, 0) + CG(1, 0, 1, 0, 1, 0)
    assert cg_simp(2 * Sum(CG(1, alpha, 0, 0, 1, alpha), (alpha, -1, 1))) == 6
    # Varshalovich 8.7.1 Eq 2
    assert cg_simp(x*Sum((-1)**(a - alpha) * CG(a, alpha, a, -alpha, c,
                   0), (alpha, -a, a))) == x*sqrt(2*a + 1)*KroneckerDelta(c, 0)
    assert cg_simp(3*Sum((-1)**(2 - alpha) * CG(
        2, alpha, 2, -alpha, 0, 0), (alpha, -2, 2))) == 3*sqrt(5)
    # Varshalovich 8.7.2 Eq 4
    assert cg_simp(Sum(CG(a, alpha, b, beta, c, gamma)*CG(a, alpha, b, beta, cp, gammap), (alpha, -a, a), (beta, -b, b))) == KroneckerDelta(c, cp)*KroneckerDelta(gamma, gammap)
    assert cg_simp(Sum(CG(a, alpha, b, beta, c, gamma)*CG(a, alpha, b, beta, c, gammap), (alpha, -a, a), (beta, -b, b))) == KroneckerDelta(gamma, gammap)
    assert cg_simp(Sum(CG(a, alpha, b, beta, c, gamma)*CG(a, alpha, b, beta, cp, gamma), (alpha, -a, a), (beta, -b, b))) == KroneckerDelta(c, cp)
    assert cg_simp(Sum(CG(
        a, alpha, b, beta, c, gamma)**2, (alpha, -a, a), (beta, -b, b))) == 1
    assert cg_simp(Sum(CG(2, alpha, 1, beta, 2, gamma)*CG(2, alpha, 1, beta, 2, gammap), (alpha, -2, 2), (beta, -1, 1))) == KroneckerDelta(gamma, gammap)
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_represent():
    x, y = symbols('x y')
    d = Density([XKet(), 0.5], [PxKet(), 0.5])
    assert (represent(0.5*(PxKet()*Dagger(PxKet()))) +
            represent(0.5*(XKet()*Dagger(XKet())))) == represent(d)

    # check for kets with expr in them
    d_with_sym = Density([XKet(x*y), 0.5], [PxKet(x*y), 0.5])
    assert (represent(0.5*(PxKet(x*y)*Dagger(PxKet(x*y)))) +
            represent(0.5*(XKet(x*y)*Dagger(XKet(x*y))))) == \
        represent(d_with_sym)

    # check when given explicit basis
    assert (represent(0.5*(XKet()*Dagger(XKet())), basis=PxOp()) +
            represent(0.5*(PxKet()*Dagger(PxKet())), basis=PxOp())) == \
        represent(d, basis=PxOp())
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
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
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_kinetic_energy():
    m, M, l1 = symbols('m M l1')
    omega = dynamicsymbols('omega')
    N = ReferenceFrame('N')
    O = Point('O')
    O.set_vel(N, 0 * N.x)
    Ac = O.locatenew('Ac', l1 * N.x)
    P = Ac.locatenew('P', l1 * N.x)
    a = ReferenceFrame('a')
    a.set_ang_vel(N, omega * N.z)
    Ac.v2pt_theory(O, N, a)
    P.v2pt_theory(O, N, a)
    Pa = Particle('Pa', P, m)
    I = outer(N.z, N.z)
    A = RigidBody('A', Ac, a, M, (I, Ac))
    assert 0 == kinetic_energy(N, Pa, A) - (M*l1**2*omega**2/2
            + 2*l1**2*m*omega**2 + omega**2/2)
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_potential_energy():
    m, M, l1, g, h, H = symbols('m M l1 g h H')
    omega = dynamicsymbols('omega')
    N = ReferenceFrame('N')
    O = Point('O')
    O.set_vel(N, 0 * N.x)
    Ac = O.locatenew('Ac', l1 * N.x)
    P = Ac.locatenew('P', l1 * N.x)
    a = ReferenceFrame('a')
    a.set_ang_vel(N, omega * N.z)
    Ac.v2pt_theory(O, N, a)
    P.v2pt_theory(O, N, a)
    Pa = Particle('Pa', P, m)
    I = outer(N.z, N.z)
    A = RigidBody('A', Ac, a, M, (I, Ac))
    Pa.set_potential_energy(m * g * h)
    A.set_potential_energy(M * g * H)
    assert potential_energy(A, Pa) == m * g * h + M * g * H
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_one_dof():
    # This is for a 1 dof spring-mass-damper case.
    # It is described in more detail in the KanesMethod docstring.
    q, u = dynamicsymbols('q u')
    qd, ud = dynamicsymbols('q u', 1)
    m, c, k = symbols('m c k')
    N = ReferenceFrame('N')
    P = Point('P')
    P.set_vel(N, u * N.x)

    kd = [qd - u]
    FL = [(P, (-k * q - c * u) * N.x)]
    pa = Particle('pa', P, m)
    BL = [pa]

    KM = KanesMethod(N, [q], [u], kd)
    KM.kanes_equations(FL, BL)
    MM = KM.mass_matrix
    forcing = KM.forcing
    rhs = MM.inv() * forcing
    assert expand(rhs[0]) == expand(-(q * k + u * c) / m)
    assert KM.linearize() == \
        (Matrix([[0, 1], [-k, -c]]), Matrix([]), Matrix([]))
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_pend():
    q, u = dynamicsymbols('q u')
    qd, ud = dynamicsymbols('q u', 1)
    m, l, g = symbols('m l g')
    N = ReferenceFrame('N')
    P = Point('P')
    P.set_vel(N, -l * u * sin(q) * N.x + l * u * cos(q) * N.y)
    kd = [qd - u]

    FL = [(P, m * g * N.x)]
    pa = Particle('pa', P, m)
    BL = [pa]

    KM = KanesMethod(N, [q], [u], kd)
    KM.kanes_equations(FL, BL)
    MM = KM.mass_matrix
    forcing = KM.forcing
    rhs = MM.inv() * forcing
    rhs.simplify()
    assert expand(rhs[0]) == expand(-g / l * sin(q))
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_rigidbody3():
    q1, q2, q3, q4 = dynamicsymbols('q1:5')
    p1, p2, p3 = symbols('p1:4')
    m = symbols('m')

    A = ReferenceFrame('A')
    B = A.orientnew('B', 'axis', [q1, A.x])
    O = Point('O')
    O.set_vel(A, q2*A.x + q3*A.y + q4*A.z)
    P = O.locatenew('P', p1*B.x + p2*B.y + p3*B.z)
    I = outer(B.x, B.x)

    rb1 = RigidBody('rb1', P, B, m, (I, P))
    # I_S/O = I_S/S* + I_S*/O
    rb2 = RigidBody('rb2', P, B, m,
                    (I + inertia_of_point_mass(m, P.pos_from(O), B), O))

    assert rb1.central_inertia == rb2.central_inertia
    assert rb1.angular_momentum(O, A) == rb2.angular_momentum(O, A)
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_simple_trace_cases_symbolic_dim():
    from sympy import symbols
    D = symbols('D')
    G = GammaMatrixHead(dim=D)

    m0, m1, m2, m3 = tensor_indices('m0:4', G.LorentzIndex)
    g = G.LorentzIndex.metric

    t = G(m0)*G(m1)
    t1 = G._trace_single_line(t)
    assert _is_tensor_eq(t1, 4 * G.LorentzIndex.metric(m0, m1))

    t = G(m0)*G(m1)*G(m2)*G(m3)
    t1 = G._trace_single_line(t)
    t2 = -4*g(m0, m2)*g(m1, m3) + 4*g(m0, m1)*g(m2, m3) + 4*g(m0, m3)*g(m1, m2)
    assert _is_tensor_eq(t1, t2)
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_partial_velocity():
    q1, q2, q3, u1, u2, u3 = dynamicsymbols('q1 q2 q3 u1 u2 u3')
    u4, u5 = dynamicsymbols('u4, u5')
    r = symbols('r')

    N = ReferenceFrame('N')
    Y = N.orientnew('Y', 'Axis', [q1, N.z])
    L = Y.orientnew('L', 'Axis', [q2, Y.x])
    R = L.orientnew('R', 'Axis', [q3, L.y])
    R.set_ang_vel(N, u1 * L.x + u2 * L.y + u3 * L.z)

    C = Point('C')
    C.set_vel(N, u4 * L.x + u5 * (Y.z ^ L.x))
    Dmc = C.locatenew('Dmc', r * L.z)
    Dmc.v2pt_theory(C, N, R)

    vel_list = [Dmc.vel(N), C.vel(N), R.ang_vel_in(N)]
    u_list = [u1, u2, u3, u4, u5]
    assert (partial_velocity(vel_list, u_list, N) ==
            [[- r*L.y, r*L.x, 0, L.x, cos(q2)*L.y - sin(q2)*L.z],
            [0, 0, 0, L.x, cos(q2)*L.y - sin(q2)*L.z],
            [L.x, L.y, L.z, 0, 0]])
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_vector_simplify():
    x, y, z, k, n, m, w, f, s, A = symbols('x, y, z, k, n, m, w, f, s, A')
    N = ReferenceFrame('N')

    test1 = (1 / x + 1 / y) * N.x
    assert (test1 & N.x) != (x + y) / (x * y)
    test1 = test1.simplify()
    assert (test1 & N.x) == (x + y) / (x * y)

    test2 = (A**2 * s**4 / (4 * pi * k * m**3)) * N.x
    test2 = test2.simplify()
    assert (test2 & N.x) == (A**2 * s**4 / (4 * pi * k * m**3))

    test3 = ((4 + 4 * x - 2 * (2 + 2 * x)) / (2 + 2 * x)) * N.x
    test3 = test3.simplify()
    assert (test3 & N.x) == 0

    test4 = ((-4 * x * y**2 - 2 * y**3 - 2 * x**2 * y) / (x + y)**2) * N.x
    test4 = test4.simplify()
    assert (test4 & N.x) == -2 * y
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_simplify():
    f, n = symbols('f, n')

    m = Matrix([[1, x], [x + 1/x, x - 1]])
    m = m.row_join(eye(m.cols))
    raw = m.rref(simplify=lambda x: x)[0]
    assert raw != m.rref(simplify=True)[0]

    M = Matrix([[            1/x + 1/y,                 (x + x*y) / x  ],
                [ (f(x) + y*f(x))/f(x), 2 * (1/n - cos(n * pi)/n) / pi ]])
    M.simplify()
    assert M == Matrix([[ (x + y)/(x * y),                        1 + y ],
                        [           1 + y, 2*((1 - 1*cos(pi*n))/(pi*n)) ]])
    eq = (1 + x)**2
    M = Matrix([[eq]])
    M.simplify()
    assert M == Matrix([[eq]])
    M.simplify(ratio=oo) == M
    assert M == Matrix([[eq.simplify(ratio=oo)]])
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_Matrix_berkowitz_charpoly():
    UA, K_i, K_w = symbols('UA K_i K_w')

    A = Matrix([[-K_i - UA + K_i**2/(K_i + K_w),       K_i*K_w/(K_i + K_w)],
                [           K_i*K_w/(K_i + K_w), -K_w + K_w**2/(K_i + K_w)]])

    charpoly = A.berkowitz_charpoly(x)

    assert charpoly == \
        Poly(x**2 + (K_i*UA + K_w*UA + 2*K_i*K_w)/(K_i + K_w)*x +
        K_i*K_w*UA/(K_i + K_w), x, domain='ZZ(K_i,K_w,UA)')

    assert type(charpoly) is PurePoly

    A = Matrix([[1, 3], [2, 0]])

    assert A.charpoly() == A.charpoly(x) == PurePoly(x**2 - x - 6)
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_jscode_loops_matrix_vector():
    n, m = symbols('n m', integer=True)
    A = IndexedBase('A')
    x = IndexedBase('x')
    y = IndexedBase('y')
    i = Idx('i', m)
    j = Idx('j', n)

    s = (
        'for (var i=0; i<m; i++){\n'
        '   y[i] = 0;\n'
        '}\n'
        'for (var i=0; i<m; i++){\n'
        '   for (var j=0; j<n; j++){\n'
        '      y[i] = x[j]*A[i*n + j] + y[i];\n'
        '   }\n'
        '}'
    )
    c = jscode(A[i, j]*x[j], assign_to=y[i])
    assert c == s
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_dummy_loops():
    # the following line could also be
    # [Dummy(s, integer=True) for s in 'im']
    # or [Dummy(integer=True) for s in 'im']
    i, m = symbols('i m', integer=True, cls=Dummy)
    x = IndexedBase('x')
    y = IndexedBase('y')
    i = Idx(i, m)

    expected = (
        'for (var i_%(icount)i=0; i_%(icount)i<m_%(mcount)i; i_%(icount)i++){\n'
        '   y[i_%(icount)i] = x[i_%(icount)i];\n'
        '}'
    ) % {'icount': i.label.dummy_index, 'mcount': m.dummy_index}
    code = jscode(x[i], assign_to=y[i])
    assert code == expected
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_jscode_loops_add():
    from sympy.tensor import IndexedBase, Idx
    from sympy import symbols
    n, m = symbols('n m', integer=True)
    A = IndexedBase('A')
    x = IndexedBase('x')
    y = IndexedBase('y')
    z = IndexedBase('z')
    i = Idx('i', m)
    j = Idx('j', n)

    s = (
        'for (var i=0; i<m; i++){\n'
        '   y[i] = x[i] + z[i];\n'
        '}\n'
        'for (var i=0; i<m; i++){\n'
        '   for (var j=0; j<n; j++){\n'
        '      y[i] = x[j]*A[i*n + j] + y[i];\n'
        '   }\n'
        '}'
    )
    c = jscode(A[i, j]*x[j] + x[i] + z[i], assign_to=y[i])
    assert c == s
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_ccode_Indexed():
    from sympy.tensor import IndexedBase, Idx
    from sympy import symbols
    n, m, o = symbols('n m o', integer=True)
    i, j, k = Idx('i', n), Idx('j', m), Idx('k', o)
    p = CCodePrinter()
    p._not_c = set()

    x = IndexedBase('x')[j]
    assert p._print_Indexed(x) == 'x[j]'
    A = IndexedBase('A')[i, j]
    assert p._print_Indexed(A) == 'A[%s]' % (m*i+j)
    B = IndexedBase('B')[i, j, k]
    assert p._print_Indexed(B) == 'B[%s]' % (i*o*m+j*o+k)

    assert p._not_c == set()
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_ccode_loops_matrix_vector():
    n, m = symbols('n m', integer=True)
    A = IndexedBase('A')
    x = IndexedBase('x')
    y = IndexedBase('y')
    i = Idx('i', m)
    j = Idx('j', n)

    s = (
        'for (int i=0; i<m; i++){\n'
        '   y[i] = 0;\n'
        '}\n'
        'for (int i=0; i<m; i++){\n'
        '   for (int j=0; j<n; j++){\n'
        '      y[i] = x[j]*A[%s] + y[i];\n' % (i*n + j) +\
        '   }\n'
        '}'
    )
    c = ccode(A[i, j]*x[j], assign_to=y[i])
    assert c == s
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_dummy_loops():
    # the following line could also be
    # [Dummy(s, integer=True) for s in 'im']
    # or [Dummy(integer=True) for s in 'im']
    i, m = symbols('i m', integer=True, cls=Dummy)
    x = IndexedBase('x')
    y = IndexedBase('y')
    i = Idx(i, m)

    expected = (
        'for (int i_%(icount)i=0; i_%(icount)i<m_%(mcount)i; i_%(icount)i++){\n'
        '   y[i_%(icount)i] = x[i_%(icount)i];\n'
        '}'
    ) % {'icount': i.label.dummy_index, 'mcount': m.dummy_index}
    code = ccode(x[i], assign_to=y[i])
    assert code == expected
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_recursive():
    from sympy import symbols, exp_polar, expand
    a, b, c = symbols('a b c', positive=True)
    r = exp(-(x - a)**2)*exp(-(x - b)**2)
    e = integrate(r, (x, 0, oo), meijerg=True)
    assert simplify(e.expand()) == (
        sqrt(2)*sqrt(pi)*(
        (erf(sqrt(2)*(a + b)/2) + 1)*exp(-a**2/2 + a*b - b**2/2))/4)
    e = integrate(exp(-(x - a)**2)*exp(-(x - b)**2)*exp(c*x), (x, 0, oo), meijerg=True)
    assert simplify(e) == (
        sqrt(2)*sqrt(pi)*(erf(sqrt(2)*(2*a + 2*b + c)/4) + 1)*exp(-a**2 - b**2
        + (2*a + 2*b + c)**2/8)/4)
    assert simplify(integrate(exp(-(x - a - b - c)**2), (x, 0, oo), meijerg=True)) == \
        sqrt(pi)/2*(1 + erf(a + b + c))
    assert simplify(integrate(exp(-(x + a + b + c)**2), (x, 0, oo), meijerg=True)) == \
        sqrt(pi)/2*(1 - erf(a + b + c))
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def runtest_ufuncify(language, backend):
    has_module('numpy')
    a, b, c = symbols('a b c')
    fabc = ufuncify([a, b, c], a*b + c, language=language, backend=backend)
    facb = ufuncify([a, c, b], a*b + c, language=language, backend=backend)
    grid = numpy.linspace(-2, 2, 50)
    for b in numpy.linspace(-5, 4, 3):
        for c in numpy.linspace(-1, 1, 3):
            expected = grid*b + c
            assert numpy.sum(numpy.abs(expected - fabc(grid, b, c))) < 1e-13
            assert numpy.sum(numpy.abs(expected - facb(grid, c, b))) < 1e-13

#
# tests of language-backend combinations
#

# f2py
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_functions():
    one_var = (acosh, ln, Heaviside, factorial, bernoulli, coth, tanh,
            sign, arg, asin, DiracDelta, re, Abs, sinh, cos, cot, acos, acot,
            gamma, bell, harmonic, LambertW, zeta, log, factorial, asinh,
            acoth, cosh, dirichlet_eta, loggamma, erf, ceiling, im, fibonacci,
            conjugate, tan, floor, atanh, sin, atan, lucas, exp)
    two_var = (rf, ff, lowergamma, chebyshevu, chebyshevt, binomial,
            atan2, polygamma, hermite, legendre, uppergamma)
    x, y, z = symbols("x,y,z")
    others = (chebyshevt_root, chebyshevu_root, Eijk(x, y, z),
            Piecewise( (0, x < -1), (x**2, x <= 1), (x**3, True)),
            assoc_legendre)
    for cls in one_var:
        check(cls)
        c = cls(x)
        check(c)
    for cls in two_var:
        check(cls)
        c = cls(x, y)
        check(c)
    for cls in others:
        check(cls)

#================== geometry ====================
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_Routine_argument_order():
    a, x, y, z = symbols('a x y z')
    expr = (x + y)*z
    raises(CodeGenArgumentListError, lambda: Routine("test", expr,
           argument_sequence=[z, x]))
    raises(CodeGenArgumentListError, lambda: Routine("test", Eq(a,
           expr), argument_sequence=[z, x, y]))
    r = Routine('test', Eq(a, expr), argument_sequence=[z, x, a, y])
    assert [ arg.name for arg in r.arguments ] == [z, x, a, y]
    assert [ type(arg) for arg in r.arguments ] == [
        InputArgument, InputArgument, OutputArgument, InputArgument  ]
    r = Routine('test', Eq(z, expr), argument_sequence=[z, x, y])
    assert [ type(arg) for arg in r.arguments ] == [
        InOutArgument, InputArgument, InputArgument ]

    from sympy.tensor import IndexedBase, Idx
    A, B = map(IndexedBase, ['A', 'B'])
    m = symbols('m', integer=True)
    i = Idx('i', m)
    r = Routine('test', Eq(A[i], B[i]), argument_sequence=[B, A, m])
    assert [ arg.name for arg in r.arguments ] == [B.label, A.label, m]

    expr = Integral(x*y*z, (x, 1, 2), (y, 1, 3))
    r = Routine('test', Eq(a, expr), argument_sequence=[z, x, a, y])
    assert [ arg.name for arg in r.arguments ] == [z, x, a, y]
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_simple_c_codegen():
    x, y, z = symbols('x,y,z')
    expr = (x + y)*z
    result = codegen(("test", expr), "C", "file", header=False, empty=False)
    expected = [
        ("file.c",
        "#include \"file.h\"\n"
        "#include <math.h>\n"
        "double test(double x, double y, double z) {\n"
        "   return z*(x + y);\n"
        "}\n"),
        ("file.h",
        "#ifndef PROJECT__FILE__H\n"
        "#define PROJECT__FILE__H\n"
        "double test(double x, double y, double z);\n"
        "#endif\n")
    ]
    assert result == expected
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_dummy_loops_c():
    from sympy.tensor import IndexedBase, Idx
    # the following line could also be
    # [Dummy(s, integer=True) for s in 'im']
    # or [Dummy(integer=True) for s in 'im']
    i, m = symbols('i m', integer=True, cls=Dummy)
    x = IndexedBase('x')
    y = IndexedBase('y')
    i = Idx(i, m)
    expected = (
        '#include "file.h"\n'
        '#include <math.h>\n'
        'void test_dummies(int m_%(mno)i, double *x, double *y) {\n'
        '   for (int i_%(ino)i=0; i_%(ino)i<m_%(mno)i; i_%(ino)i++){\n'
        '      y[i_%(ino)i] = x[i_%(ino)i];\n'
        '   }\n'
        '}\n'
    ) % {'ino': i.label.dummy_index, 'mno': m.dummy_index}
    r = Routine('test_dummies', Eq(y[i], x[i]))
    c = CCodeGen()
    code = get_string(c.dump_c, [r])
    assert code == expected
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_output_arg_c():
    from sympy import sin, cos, Equality
    x, y, z = symbols("x,y,z")
    r = Routine("foo", [Equality(y, sin(x)), cos(x)])
    c = CCodeGen()
    result = c.write([r], "test", header=False, empty=False)
    assert result[0][0] == "test.c"
    expected = (
        '#include "test.h"\n'
        '#include <math.h>\n'
        'double foo(double x, double &y) {\n'
        '   y = sin(x);\n'
        '   return cos(x);\n'
        '}\n'
    )
    assert result[0][1] == expected
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_simple_f_code():
    x, y, z = symbols('x,y,z')
    expr = (x + y)*z
    routine = Routine("test", expr)
    code_gen = FCodeGen()
    source = get_string(code_gen.dump_f95, [routine])
    expected = (
        "REAL*8 function test(x, y, z)\n"
        "implicit none\n"
        "REAL*8, intent(in) :: x\n"
        "REAL*8, intent(in) :: y\n"
        "REAL*8, intent(in) :: z\n"
        "test = z*(x + y)\n"
        "end function\n"
    )
    assert source == expected
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_f_code_argument_order():
    x, y, z = symbols('x,y,z')
    expr = x + y
    routine = Routine("test", expr, argument_sequence=[z, x, y])
    code_gen = FCodeGen()
    source = get_string(code_gen.dump_f95, [routine])
    expected = (
        "REAL*8 function test(z, x, y)\n"
        "implicit none\n"
        "REAL*8, intent(in) :: z\n"
        "REAL*8, intent(in) :: x\n"
        "REAL*8, intent(in) :: y\n"
        "test = x + y\n"
        "end function\n"
    )
    assert source == expected
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_simple_f_header():
    x, y, z = symbols('x,y,z')
    expr = (x + y)*z
    routine = Routine("test", expr)
    code_gen = FCodeGen()
    source = get_string(code_gen.dump_h, [routine])
    expected = (
        "interface\n"
        "REAL*8 function test(x, y, z)\n"
        "implicit none\n"
        "REAL*8, intent(in) :: x\n"
        "REAL*8, intent(in) :: y\n"
        "REAL*8, intent(in) :: z\n"
        "end function\n"
        "end interface\n"
    )
    assert source == expected
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_simple_f_codegen():
    x, y, z = symbols('x,y,z')
    expr = (x + y)*z
    result = codegen(
        ("test", expr), "F95", "file", header=False, empty=False)
    expected = [
        ("file.f90",
        "REAL*8 function test(x, y, z)\n"
        "implicit none\n"
        "REAL*8, intent(in) :: x\n"
        "REAL*8, intent(in) :: y\n"
        "REAL*8, intent(in) :: z\n"
        "test = z*(x + y)\n"
        "end function\n"),
        ("file.h",
        "interface\n"
        "REAL*8 function test(x, y, z)\n"
        "implicit none\n"
        "REAL*8, intent(in) :: x\n"
        "REAL*8, intent(in) :: y\n"
        "REAL*8, intent(in) :: z\n"
        "end function\n"
        "end interface\n")
    ]
    assert result == expected
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_dummy_loops_f95():
    from sympy.tensor import IndexedBase, Idx
    # the following line could also be
    # [Dummy(s, integer=True) for s in 'im']
    # or [Dummy(integer=True) for s in 'im']
    i, m = symbols('i m', integer=True, cls=Dummy)
    x = IndexedBase('x')
    y = IndexedBase('y')
    i = Idx(i, m)
    expected = (
        'subroutine test_dummies(m_%(mcount)i, x, y)\n'
        'implicit none\n'
        'INTEGER*4, intent(in) :: m_%(mcount)i\n'
        'REAL*8, intent(in), dimension(1:m_%(mcount)i) :: x\n'
        'REAL*8, intent(out), dimension(1:m_%(mcount)i) :: y\n'
        'INTEGER*4 :: i_%(icount)i\n'
        'do i_%(icount)i = 1, m_%(mcount)i\n'
        '   y(i_%(icount)i) = x(i_%(icount)i)\n'
        'end do\n'
        'end subroutine\n'
    ) % {'icount': i.label.dummy_index, 'mcount': m.dummy_index}
    r = Routine('test_dummies', Eq(y[i], x[i]))
    c = FCodeGen()
    code = get_string(c.dump_f95, [r])
    assert code == expected
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def test_output_arg_f():
    from sympy import sin, cos, Equality
    x, y, z = symbols("x,y,z")
    r = Routine("foo", [Equality(y, sin(x)), cos(x)])
    c = FCodeGen()
    result = c.write([r], "test", header=False, empty=False)
    assert result[0][0] == "test.f90"
    assert result[0][1] == (
        'REAL*8 function foo(x, y)\n'
        'implicit none\n'
        'REAL*8, intent(in) :: x\n'
        'REAL*8, intent(out) :: y\n'
        'y = sin(x)\n'
        'foo = cos(x)\n'
        'end function\n'
    )
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def Simple_manifold_with_scalar_function_derivative():
    Print_Function()
    coords = (x, y, z) = symbols('x y z')
    basis = (e1, e2, e3, grad) = MV.setup('e_1 e_2 e_3', metric='[1,1,1]', coords=coords)
    # Define surface
    mfvar = (u, v) = symbols('u v')
    X = u*e1 + v*e2 + (u**2 + v**2)*e3
    print(X)
    MF = Manifold(X, mfvar)

    # Define field on the surface.
    g = (v + 1)*log(u)

    # Method 1: Using old Manifold routines.
    VectorDerivative = (MF.rbasis[0]/MF.E_sq)*diff(g, u) + (MF.rbasis[1]/MF.E_sq)*diff(g, v)
    print('Vector derivative =', VectorDerivative.subs({u: 1, v: 0}))

    # Method 2: Using new Manifold routines.
    dg = MF.Grad(g)
    print('Vector derivative =', dg.subs({u: 1, v: 0}))
    return
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def derivatives_in_spherical_coordinates():
    Print_Function()
    X = (r, th, phi) = symbols('r theta phi')
    curv = [[r*cos(phi)*sin(th), r*sin(phi)*sin(th), r*cos(th)], [1, r, r*sin(th)]]
    (er, eth, ephi, grad) = MV.setup('e_r e_theta e_phi', metric='[1,1,1]', coords=X, curv=curv)

    f = MV('f', 'scalar', fct=True)
    A = MV('A', 'vector', fct=True)
    B = MV('B', 'grade2', fct=True)

    print('f =', f)
    print('A =', A)
    print('B =', B)

    print('grad*f =', grad*f)
    print('grad|A =', grad | A)
    print('-I*(grad^A) =', -MV.I*(grad ^ A))
    print('grad^B =', grad ^ B)
    return
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def Dirac_Equation_in_Geometric_Calculus():
    Print_Function()
    vars = symbols('t x y z')
    (g0, g1, g2, g3, grad) = MV.setup('gamma*t|x|y|z', metric='[1,-1,-1,-1]', coords=vars)
    I = MV.I

    (m, e) = symbols('m e')

    psi = MV('psi', 'spinor', fct=True)
    A = MV('A', 'vector', fct=True)
    sig_z = g3*g0

    print('\\text{4-Vector Potential\\;\\;}\\bm{A} =', A)
    print('\\text{8-component real spinor\\;\\;}\\bm{\\psi} =', psi)

    dirac_eq = (grad*psi)*I*sig_z - e*A*psi - m*psi*g0
    dirac_eq.simplify()

    dirac_eq.Fmt(3, r'%\text{Dirac Equation\;\;}\nabla \bm{\psi} I \sigma_{z}-e\bm{A}\bm{\psi}-m\bm{\psi}\gamma_{t} = 0')

    return
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def derivatives_in_spherical_coordinates():
    Print_Function()
    X = (r, th, phi) = symbols('r theta phi')
    curv = [[r*cos(phi)*sin(th), r*sin(phi)*sin(th), r*cos(th)], [1, r, r*sin(th)]]
    (er, eth, ephi, grad) = MV.setup('e_r e_theta e_phi', metric='[1,1,1]', coords=X, curv=curv)

    f = MV('f', 'scalar', fct=True)
    A = MV('A', 'vector', fct=True)
    B = MV('B', 'grade2', fct=True)

    print('f =', f)
    print('A =', A)
    print('B =', B)

    print('grad*f =', grad*f)
    print('grad|A =', grad | A)
    print('-I*(grad^A) =', (-MV.I*(grad ^ A)).simplify())
    print('grad^B =', grad ^ B)
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def MV_setup_options():
    Print_Function()

    (e1, e2, e3) = MV.setup('e_1 e_2 e_3', '[1,1,1]')
    v = MV('v', 'vector')
    print(v)

    (e1, e2, e3) = MV.setup('e*1|2|3', '[1,1,1]')
    v = MV('v', 'vector')
    print(v)

    (e1, e2, e3) = MV.setup('e*x|y|z', '[1,1,1]')
    v = MV('v', 'vector')
    print(v)

    coords = symbols('x y z')
    (e1, e2, e3, grad) = MV.setup('e', '[1,1,1]', coords=coords)
    v = MV('v', 'vector')
    print(v)

    return
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def main():
    Format()
    a = Matrix(2, 2, (1, 2, 3, 4))
    b = Matrix(2, 1, (5, 6))
    c = a*b
    print(a, b, '=', c)

    x, y = symbols('x, y')

    d = Matrix(1, 2, (x**3, y**3))
    e = Matrix(2, 2, (x**2, 2*x*y, 2*x*y, y**2))
    f = d*e

    print('%', d, e, '=', f)

    xdvi()
    return
项目:zippy    作者:securesystemslab    | 项目源码 | 文件源码
def Simple_manifold_with_scalar_function_derivative():
    Print_Function()
    coords = (x, y, z) = symbols('x y z')
    basis = (e1, e2, e3, grad) = MV.setup('e_1 e_2 e_3', metric='[1,1,1]', coords=coords)
    # Define surface
    mfvar = (u, v) = symbols('u v')
    X = u*e1 + v*e2 + (u**2 + v**2)*e3
    print('\\f{X}{u,v} =', X)
    MF = Manifold(X, mfvar)
    (eu, ev) = MF.Basis()
    # Define field on the surface.
    g = (v + 1)*log(u)

    print('\\f{g}{u,v} =', g)

    # Method 1: Using old Manifold routines.
    VectorDerivative = (MF.rbasis[0]/MF.E_sq)*diff(g, u) + (MF.rbasis[1]/MF.E_sq)*diff(g, v)
    print('\\eval{\\nabla g}{u=1,v=0} =', VectorDerivative.subs({u: 1, v: 0}))

    # Method 2: Using new Manifold routines.
    dg = MF.Grad(g)
    print('\\eval{\\f{Grad}{g}}{u=1,v=0} =', dg.subs({u: 1, v: 0}))
    dg = MF.grad*g
    print('\\eval{\\nabla g}{u=1,v=0} =', dg.subs({u: 1, v: 0}))
    return