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)
python类symbols()的实例源码
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
