我们从Python开源项目中,提取了以下50个代码示例,用于说明如何使用sympy.symbols()。
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)
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 #
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!")
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'))
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))
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
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)
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'))")
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
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)))
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)))
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)
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())
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_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)
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
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([]))
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))
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)
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)
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]])
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
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)]])
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)
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
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
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
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()
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
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
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))
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
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 ====================
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]
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
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
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
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
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
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
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
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
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' )
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
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
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
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)
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
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
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