def clean_sh(sh):
"""Turns symbolic solid harmonic functions into string representations
to be using in generating basis functions.
Args
sh (OrderedDict): Output from exatomic.algorithms.basis.solid_harmonics
Returns
clean (OrderedDict): cleaned strings
"""
_replace = {'x': '{x}', 'y': '{y}', 'z': '{z}', ' - ': ' -'}
_repatrn = re.compile('|'.join(_replace.keys()))
clean = OrderedDict()
for key, sym in sh.items():
if isinstance(sym, (Mul, Add)):
string = str(sym.expand()).replace(' + ', ' ')#.replace(' - ', ' -')
string = _repatrn.sub(lambda x: _replace[x.group(0)], string)
clean[key] = [pre + '*' for pre in string.split()]
else: clean[key] = ['']
return clean
python类Mul()的实例源码
def _generate_outer_prod(self, arg1, arg2):
c_part1, nc_part1 = arg1.args_cnc()
c_part2, nc_part2 = arg2.args_cnc()
if ( len(nc_part1) == 0 or
len(nc_part2) == 0 ):
raise ValueError('Atleast one-pair of'
' Non-commutative instance required'
' for outer product.')
# Muls of Tensor Products should be expanded
# before this function is called
if (isinstance(nc_part1[0], TensorProduct) and
len(nc_part1) == 1 and len(nc_part2) == 1):
op = tensor_product_simp(nc_part1[0] * Dagger(nc_part2[0]))
else:
op = Mul(*nc_part1) * Dagger(Mul(*nc_part2))
return Mul(*c_part1)*Mul(*c_part2)*op
def _normal_ordered_form_terms(expr, independent=False, recursive_limit=10,
_recursive_depth=0):
"""
Helper function for normal_ordered_form: loop through each term in an
addition expression and call _normal_ordered_form_factor to perform the
factor to an normally ordered expression.
"""
new_terms = []
for term in expr.args:
if isinstance(term, Mul):
new_term = _normal_ordered_form_factor(
term, recursive_limit=recursive_limit,
_recursive_depth=_recursive_depth, independent=independent)
new_terms.append(new_term)
else:
new_terms.append(term)
return Add(*new_terms)
def _normal_order_terms(expr, recursive_limit=10, _recursive_depth=0):
"""
Helper function for normal_order: look through each term in an addition
expression and call _normal_order_factor to perform the normal ordering
on the factors.
"""
new_terms = []
for term in expr.args:
if isinstance(term, Mul):
new_term = _normal_order_factor(term,
recursive_limit=recursive_limit,
_recursive_depth=_recursive_depth)
new_terms.append(new_term)
else:
new_terms.append(term)
return Add(*new_terms)
def _print_Mul(self, expr):
"""Print a Mul object.
:param expr: The expression.
:rtype : str
:return: The printed string.
"""
assert isinstance(expr, _sympy.Mul)
# Get the precedence.
prec = _sympy_precedence.precedence(expr)
# Get commutative factors and non-commutative factors.
c, nc = expr.args_cnc()
# Print.
res = super(_MEXPPrinter, self)._print_Mul(expr.func(*c))
if nc:
res += '*'
res += '^'.join(self.parenthesize(a, prec) for a in nc)
return res
def expr(self):
"""
Process the expression, by returning a sympy expression
"""
if DEBUG:
print("> ArithmeticExpression ")
ret = self.op[0].expr
for operation, operand in zip(self.op[1::2], self.op[2::2]):
if operation == '+':
ret = Add(ret, operand.expr)
else:
a = Mul(-1, operand.expr)
ret = Add(ret, a)
return ret
def convertUnitTo(self, expr, newUnit):
if expr.astUnit.isNull():
raise XXXSyntaxError("Can't convert a null unit!")
elif not expr.astUnit.rawUnit.isCompatibleWith(newUnit):
raise XXXSyntaxError("Incompatible unit conversion requested.")
else:
factor = expr.astUnit.rawUnit.convertTo(newUnit)
return AST(sympy.Mul(factor,expr.sympy), ASTUnit(newUnit, explicit=False))
#def astToVar(self, var, ast):
# self.currentSubsystem().ssa[var] = ast
# self.currentScope().addInstruction(var)
# return (var, ast.astUnit)
#def astToTemp(self, ast):
# return self.astToVar(self.newTempVar(),ast)
#def astToSymbol(self, name, ast):
# (var, astUnit) = self.astToVar(Symbol(name),ast)
# self.currentScope().setSymbol(name, var)
# return (var, astUnit)
def doit(self, **hints):
expr = self.args[0]
condition = self._condition
if not expr.has(RandomSymbol):
return expr
if isinstance(expr, Add):
return Add(*[Expectation(a, condition=condition).doit() for a in expr.args])
elif isinstance(expr, Mul):
rv = []
nonrv = []
for a in expr.args:
if isinstance(a, RandomSymbol) or a.has(RandomSymbol):
rv.append(a)
else:
nonrv.append(a)
return Mul(*nonrv)*Expectation(Mul(*rv), condition=condition)
return self
def _expand_single_argument(cls, expr):
# return (coefficient, random_symbol) pairs:
if isinstance(expr, RandomSymbol):
return [(S.One, expr)]
elif isinstance(expr, Add):
outval = []
for a in expr.args:
if isinstance(a, Mul):
outval.append(cls._get_mul_nonrv_rv_tuple(a))
elif isinstance(a, RandomSymbol):
outval.append((S.One, a))
return outval
elif isinstance(expr, Mul):
return [cls._get_mul_nonrv_rv_tuple(expr)]
elif expr.has(RandomSymbol):
return [(S.One, expr)]
def _collect_factor_and_dimension(expr):
if isinstance(expr, Quantity):
return expr.scale_factor, expr.dimension
elif isinstance(expr, Mul):
factor = 1
dimension = 1
for arg in expr.args:
arg_factor, arg_dim = Quantity._collect_factor_and_dimension(arg)
factor *= arg_factor
dimension *= arg_dim
return factor, dimension
elif isinstance(expr, Pow):
factor, dim = Quantity._collect_factor_and_dimension(expr.base)
return factor ** expr.exp, dim ** expr.exp
elif isinstance(expr, Add):
raise NotImplementedError
else:
return 1, 1
def _generate_outer_prod(self, arg1, arg2):
c_part1, nc_part1 = arg1.args_cnc()
c_part2, nc_part2 = arg2.args_cnc()
if ( len(nc_part1) == 0 or
len(nc_part2) == 0 ):
raise ValueError('Atleast one-pair of'
' Non-commutative instance required'
' for outer product.')
# Muls of Tensor Products should be expanded
# before this function is called
if (isinstance(nc_part1[0], TensorProduct) and
len(nc_part1) == 1 and len(nc_part2) == 1):
op = tensor_product_simp(nc_part1[0] * Dagger(nc_part2[0]))
else:
op = Mul(*nc_part1) * Dagger(Mul(*nc_part2))
return Mul(*c_part1)*Mul(*c_part2)*op
def _normal_order_terms(expr, recursive_limit=10, _recursive_depth=0):
"""
Helper function for normal_order: look through each term in an addition
expression and call _normal_order_factor to perform the normal ordering
on the factors.
"""
new_terms = []
for term in expr.args:
if isinstance(term, Mul):
new_term = _normal_order_factor(term,
recursive_limit=recursive_limit,
_recursive_depth=_recursive_depth)
new_terms.append(new_term)
else:
new_terms.append(term)
return Add(*new_terms)
def eval(cls, arg):
if not arg.is_Atom:
c, arg_ = factor_terms(arg).as_coeff_Mul()
if arg_.is_Mul:
arg_ = Mul(*[a if (sign(a) not in (-1, 1)) else
sign(a) for a in arg_.args])
arg_ = sign(c)*arg_
else:
arg_ = arg
if arg_.atoms(AppliedUndef):
return
x, y = re(arg_), im(arg_)
rv = atan2(y, x)
if rv.is_number:
return rv
if arg_ != arg:
return cls(arg_, evaluate=False)
def _fix_ma(self, species = None):
"""Check the numerator of the reaction rate, and adds species
to reactant and product if they divide the numerator but their
stoichiometry does not match the degree in the rate."""
remainder = self.kinetic_param.as_numer_denom()[0].cancel()
if remainder.func.__name__ == 'Mul':
mulargs = list(remainder.args) + [i.args[0] for i in remainder.args if i.func.__name__ == 'Mul'] \
+ [i.args[0] for i in remainder.args if i.func.__name__ == 'Pow']
while any(sp.Symbol(s) in mulargs for s in species):
for s in species:
if sp.Symbol(s) in mulargs:
if s in self.reactant: self.reactant[s] = self.reactant[s] + 1
else: self.reactant[s] = 1
if s in self.product: self.product[s] = self.product[s] + 1
else: self.product[s] = 1
remainder = (remainder / sp.Symbol(s)).factor()
if remainder.func.__name__ == 'Mul':
mulargs = list(remainder.args) + [i.args[0] for i in remainder.args if i.func.__name__ == 'Mul'] \
+ [i.args[0] for i in remainder.args if i.func.__name__ == 'Pow']
else: mulargs = []
# update the kinetic parameter
self.__kinetic_param = (self.rate / self.reactant.ma()).cancel()
def _fix_denom(self, species):
"""Remove species that are involved in both reactant and product,
if their concentration divides the denominator of the rate."""
remainder = self.kinetic_param.as_numer_denom()[1].cancel()
#if remainder.func.__name__ == 'Mul':
if remainder != 1:
mulargs = [remainder] + list(remainder.args) + [i.args[0] for i in remainder.args if i.func.__name__ == 'Mul'] \
+ [i.args[0] for i in remainder.args if i.func.__name__ == 'Pow']
while any(sp.Symbol(s) in mulargs and s in self.reactant and s in self.product for s in species):
for s in species:
if sp.Symbol(s) in mulargs and s in self.reactant and s in self.product:
if self.reactant[s] == 1: del self.reactant[s]
else: self.reactant[s] = self.reactant[s] - 1
if self.product[s] == 1: del self.product[s]
else: self.product[s] = self.product[s] - 1
remainder = (remainder / sp.Symbol(s)).factor()
if remainder.func.__name__ == 'Mul':
mulargs = list(remainder.args) + [i.args[0] for i in remainder.args if i.func.__name__ == 'Mul'] \
+ [i.args[0] for i in remainder.args if i.func.__name__ == 'Pow']
else:
if str(remainder) in species: mulargs = [remainder]
else: mulargs = []
# update the kinetic parameter
self._kinetic_param = self.rate / self.reactant.ma()
def freeze_expression(expr):
"""
Reconstruct ``expr`` turning all :class:`sympy.Mul` and :class:`sympy.Add`
into, respectively, :class:`devito.Mul` and :class:`devito.Add`.
"""
if expr.is_Atom or expr.is_Indexed:
return expr
elif expr.is_Add:
rebuilt_args = [freeze_expression(e) for e in expr.args]
return Add(*rebuilt_args, evaluate=False)
elif expr.is_Mul:
rebuilt_args = [freeze_expression(e) for e in expr.args]
return Mul(*rebuilt_args, evaluate=False)
elif expr.is_Equality:
rebuilt_args = [freeze_expression(e) for e in expr.args]
return Eq(*rebuilt_args, evaluate=False)
else:
return expr.func(*[freeze_expression(e) for e in expr.args])
def car2sph(sh, cart):
"""
Turns symbolic solid harmonic functions into a dictionary of
arrays containing cartesian to spherical transformation matrices.
Args
sh (OrderedDict): the result of solid_harmonics(l_tot)
cart (dict): dictionary of l, cartesian l, m, n ordering
"""
keys, conv = {}, {}
prevL, mscnt = 0, 0
for L in cart:
for idx, (l, m, n) in enumerate(cart[L]):
key = ''
if l: key += 'x'
if l > 1: key += str(l)
if m: key += 'y'
if m > 1: key += str(m)
if n: key += 'z'
if n > 1: key += str(n)
keys[key] = idx
# TODO: six compatibility
for key, sym in sh.items():
L = key[0]
mscnt = mscnt if prevL == L else 0
conv.setdefault(L, np.zeros((cart_lml_count[L],
spher_lml_count[L]),
dtype=np.float64))
if isinstance(sym, (Mul, Add)):
string = (str(sym.expand())
.replace(' + ', ' ')
.replace(' - ', ' -'))
for chnk in string.split():
pre, exp = chnk.split('*', 1)
if L == 1: conv[L] = np.array(cart[L])
else: conv[L][keys[exp.replace('*', '')], mscnt] = float(pre)
prevL = L
mscnt += 1
conv[0] = np.array([[1]])
return conv
def _print_Mul(self, expr):
prec = precedence(expr)
pows = [i for i in expr.args if i.is_Pow and i.exp < 0]
if len(pows) > 1:
raise NotImplementedError("Need exactly one inverted Pow, not %s" % len(pows))
if not pows:
no_autoeye = self.__class__({**self._settings, 'use_autoeye': False})
num_terms = [no_autoeye._print(no_autoeye.parenthesize(i, prec)) for i in
expr.args if i.is_number]
mat_terms = [self._print(self.parenthesize(i, prec)) for i in
expr.args if not i.is_number]
if len(mat_terms) >= 2 and self._settings['py_solve']:
raise NotImplementedError("matrix multiplication is not yet supported with py_solve")
if num_terms and mat_terms:
return '*'.join(num_terms) + '*' + '@'.join(mat_terms)
else:
if self._settings['use_autoeye']:
if num_terms:
return ('autoeye(%s)' % '*'.join(num_terms)) + '@'.join(mat_terms)
return '@'.join(mat_terms)
return '*'.join(num_terms) + '@'.join(mat_terms)
[pow] = pows
rest = Mul(*[i for i in expr.args if i != pow])
return 'solve(%s, %s)' % (self._print(1/pow), self._print(rest))
def _rearrange_args(l):
""" this just moves the last arg to first position
to enable expansion of args
A,B,A ==> A**2,B
"""
if len(l) == 1:
return l
x = list(l[-1:])
x.extend(l[0:-1])
return Mul(*x).args
def permute(self, pos):
""" Permute the arguments cyclically.
Parameters
==========
pos : integer, if positive, shift-right, else shift-left
Examples
=========
>>> from sympy.core.trace import Tr
>>> from sympy import symbols
>>> A, B, C, D = symbols('A B C D', commutative=False)
>>> t = Tr(A*B*C*D)
>>> t.permute(2)
Tr(C*D*A*B)
>>> t.permute(-2)
Tr(C*D*A*B)
"""
if pos > 0:
pos = pos % len(self.args[0].args)
else:
pos = -(abs(pos) % len(self.args[0].args))
args = list(self.args[0].args[-pos:] + self.args[0].args[0:-pos])
return Tr(Mul(*(args)))
def _hashable_content(self):
if isinstance(self.args[0], Mul):
args = _cycle_permute(_rearrange_args(self.args[0].args))
else:
args = [self.args[0]]
return tuple(args) + (self.args[1], )
def qubit_to_matrix(qubit, format='sympy'):
"""Coverts an Add/Mul of Qubit objects into it's matrix representation
This function is the inverse of ``matrix_to_qubit`` and is a shorthand
for ``represent(qubit)``.
"""
return represent(qubit, format=format)
#-----------------------------------------------------------------------------
# Measurement
#-----------------------------------------------------------------------------
def __new__(cls, *args):
# args should be a tuple - a variable length argument list
obj = Basic.__new__(cls, *args)
obj._circuit = Mul(*args)
obj._rules = generate_gate_rules(args)
obj._eq_ids = generate_equivalent_ids(args)
return obj
def random_circuit(ngates, nqubits, gate_space=(X, Y, Z, S, T, H, CNOT, SWAP)):
"""Return a random circuit of ngates and nqubits.
This uses an equally weighted sample of (X, Y, Z, S, T, H, CNOT, SWAP)
gates.
Parameters
----------
ngates : int
The number of gates in the circuit.
nqubits : int
The number of qubits in the circuit.
gate_space : tuple
A tuple of the gate classes that will be used in the circuit.
Repeating gate classes multiple times in this tuple will increase
the frequency they appear in the random circuit.
"""
qubit_space = range(nqubits)
result = []
for i in xrange(ngates):
g = random.choice(gate_space)
if g == CNotGate or g == SwapGate:
qubits = random.sample(qubit_space, 2)
g = g(*qubits)
else:
qubit = random.choice(qubit_space)
g = g(qubit)
result.append(g)
return Mul(*result)
def test_random_reduce():
x = X(0)
y = Y(0)
z = Z(0)
h = H(0)
cnot = CNOT(1, 0)
cgate_z = CGate((0,), Z(1))
gate_list = [x, y, z]
ids = list(bfs_identity_search(gate_list, 1, max_depth=4))
circuit = (x, y, h, z, cnot)
assert random_reduce(circuit, []) == circuit
assert random_reduce(circuit, ids) == circuit
seq = [2, 11, 9, 3, 5]
circuit = (x, y, z, x, y, h)
assert random_reduce(circuit, ids, seed=seq) == (x, y, h)
circuit = (x, x, y, y, z, z)
assert random_reduce(circuit, ids, seed=seq) == (x, x, y, y)
seq = [14, 13, 0]
assert random_reduce(circuit, ids, seed=seq) == (y, y, z, z)
gate_list = [x, y, z, h, cnot, cgate_z]
ids = list(bfs_identity_search(gate_list, 2, max_depth=4))
seq = [25]
circuit = (x, y, z, y, h, y, h, cgate_z, h, cnot)
expected = (x, y, z, cgate_z, h, cnot)
assert random_reduce(circuit, ids, seed=seq) == expected
circuit = Mul(*circuit)
assert random_reduce(circuit, ids, seed=seq) == expected
def test_generate_equivalent_ids_1():
# Test with tuples
(x, y, z, h) = create_gate_sequence()
assert generate_equivalent_ids((x,)) == set([(x,)])
assert generate_equivalent_ids((x, x)) == set([(x, x)])
assert generate_equivalent_ids((x, y)) == set([(x, y), (y, x)])
gate_seq = (x, y, z)
gate_ids = set([(x, y, z), (y, z, x), (z, x, y), (z, y, x),
(y, x, z), (x, z, y)])
assert generate_equivalent_ids(gate_seq) == gate_ids
gate_ids = set([Mul(x, y, z), Mul(y, z, x), Mul(z, x, y),
Mul(z, y, x), Mul(y, x, z), Mul(x, z, y)])
assert generate_equivalent_ids(gate_seq, return_as_muls=True) == gate_ids
gate_seq = (x, y, z, h)
gate_ids = set([(x, y, z, h), (y, z, h, x),
(h, x, y, z), (h, z, y, x),
(z, y, x, h), (y, x, h, z),
(z, h, x, y), (x, h, z, y)])
assert generate_equivalent_ids(gate_seq) == gate_ids
gate_seq = (x, y, x, y)
gate_ids = set([(x, y, x, y), (y, x, y, x)])
assert generate_equivalent_ids(gate_seq) == gate_ids
cgate_y = CGate((1,), y)
gate_seq = (y, cgate_y, y, cgate_y)
gate_ids = set([(y, cgate_y, y, cgate_y), (cgate_y, y, cgate_y, y)])
assert generate_equivalent_ids(gate_seq) == gate_ids
cnot = CNOT(1, 0)
cgate_z = CGate((0,), Z(1))
gate_seq = (cnot, h, cgate_z, h)
gate_ids = set([(cnot, h, cgate_z, h), (h, cgate_z, h, cnot),
(h, cnot, h, cgate_z), (cgate_z, h, cnot, h)])
assert generate_equivalent_ids(gate_seq) == gate_ids
def test_random_circuit():
c = random_circuit(10, 3)
assert isinstance(c, Mul)
m = represent(c, nqubits=3)
assert m.shape == (8, 8)
assert isinstance(m, Matrix)
def normal_order(expr, recursive_limit=10, _recursive_depth=0):
"""Normal order an expression with bosonic or fermionic operators. Note
that this normal order is not equivalent to the original expression, but
the creation and annihilation operators in each term in expr is reordered
so that the expression becomes normal ordered.
Parameters
==========
expr : expression
The expression to normal order.
recursive_limit : int (default 10)
The number of allowed recursive applications of the function.
Examples
========
>>> from sympy.physics.quantum import Dagger
>>> from sympy.physics.quantum.boson import BosonOp
>>> from sympy.physics.quantum.operatorordering import normal_order
>>> a = BosonOp("a")
>>> normal_order(a * Dagger(a))
Dagger(a)*a
"""
if _recursive_depth > recursive_limit:
warn.warning("Warning: too many recursions, aborting")
return expr
if isinstance(expr, Add):
return _normal_order_terms(expr, recursive_limit=recursive_limit,
_recursive_depth=_recursive_depth)
elif isinstance(expr, Mul):
return _normal_order_factor(expr, recursive_limit=recursive_limit,
_recursive_depth=_recursive_depth)
else:
return expr
def _gates(self):
"""Create a list of all gates in the circuit plot."""
gates = []
if isinstance(self.circuit, Mul):
for g in reversed(self.circuit.args):
if isinstance(g, Gate):
gates.append(g)
elif isinstance(self.circuit, Gate):
gates.append(self.circuit)
return gates
def eval(cls, arg):
from sympy.simplify.simplify import signsimp
if hasattr(arg, '_eval_Abs'):
obj = arg._eval_Abs()
if obj is not None:
return obj
# handle what we can
arg = signsimp(arg, evaluate=False)
if arg.is_Mul:
known = []
unk = []
for t in arg.args:
tnew = cls(t)
if tnew.func is cls:
unk.append(tnew.args[0])
else:
known.append(tnew)
known = Mul(*known)
unk = cls(Mul(*unk), evaluate=False) if unk else S.One
return known*unk
if arg is S.NaN:
return S.NaN
if arg.is_zero: # it may be an Expr that is zero
return S.Zero
if arg.is_nonnegative:
return arg
if arg.is_nonpositive:
return -arg
if arg.is_imaginary:
arg2 = -S.ImaginaryUnit * arg
if arg2.is_nonnegative:
return arg2
if arg.is_real is False and arg.is_imaginary is False:
from sympy import expand_mul
return sqrt( expand_mul(arg * arg.conjugate()) )
if arg.is_Pow:
base, exponent = arg.as_base_exp()
if exponent.is_even and base.is_real:
return arg
if exponent.is_integer and base is S.NegativeOne:
return S.One
