Repository: rigetticomputing/pyquil
Pull request #357 (2018-03-14 20:28:05)

PauliTerms remember their qubit ordering
See changelog and #355 for motivation.

In particular, this will be important for using QubitPlaceholder's cc #354 as well as for tailoring pyquil programs to physical devices. 

This PR attempts to keep as much behavior the same as possible. Pay close attention to the changes in the tests for what behaviour changes. It should only be in cases where a string is being generated... either the str() of a paulisum or a program resulting from exponentiating a pualisum

willzeng: @mpharrigan it looks like it might be a semaphore problem with the build? is that the case
mpharrigan: It's failing on flake8 style check and not running the actual tests
mpharrigan: @ncrubin @stevenheidel review?

Patch:
diff --git a/docs/source/changes.rst b/docs/source/changes.rst
new file mode 100644
index 0000000..444e7da
--- /dev/null
+++ b/docs/source/changes.rst
@@ -0,0 +1,15 @@
+Changelog
+=========
+
+v1.9
+----
+
+ - :py:class:`PauliTerm` now remembers the order of its operations. ``sX(1)*sZ(2)`` will compile
+   to different Quil code than ``sZ(2)*sX(1)``, although the terms will still be equal according
+   to the ``__eq__`` method. During :py:class:`PauliSum` combination
+   of like terms, a warning will be emitted if two terms are combined that have different orders
+   of operation.
+ - :py:func:`PauliTerm.id()` takes an optional argument ``sort_ops`` which defaults to True for
+   backwards compatibility. However, this function should not be used for comparing term-type like
+   it has been used previously. Use :py:func:`PauliTerm.operations_as_set()` instead. In the future,
+   ``sort_ops`` will default to False and will eventually be removed.
diff --git a/docs/source/index.rst b/docs/source/index.rst
index 85f8218..f51dc05 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -56,6 +56,7 @@ Contents
    compiler
    noise
    modules
+   changes
 
 
 Indices and Tables
diff --git a/pyquil/paulis.py b/pyquil/paulis.py
index bb44bb4..52762d7 100644
--- a/pyquil/paulis.py
+++ b/pyquil/paulis.py
@@ -24,7 +24,7 @@ import copy
 from .quil import Program
 from .gates import H, RZ, RX, CNOT, X, PHASE, STANDARD_GATES
 from numbers import Number
-from collections import Sequence
+from collections import Sequence, OrderedDict, defaultdict
 import warnings
 from six import integer_types as six_integer_types
 from six.moves import range
@@ -50,6 +50,12 @@ class UnequalLengthWarning(Warning):
 integer_types = six_integer_types + (np.int64, np.int32, np.int16, np.int8)
 """Explicitly include numpy integer dtypes (for python 3)."""
 
+HASH_PRECISION = 1e6
+"""The precision used when hashing terms to check equality. The simplify() method
+uses np.isclose() for coefficient comparisons to 0 which has its own default precision. We
+can't use np.isclose() for hashing terms though.
+"""
+
 
 class PauliTerm(object):
     """A term is a product of Pauli operators operating on different qubits.
@@ -66,30 +72,44 @@ class PauliTerm(object):
         assert op in PAULI_OPS
         assert isinstance(index, integer_types) and index >= 0
 
-        self._ops = {}
+        self._ops = OrderedDict()
         if op != "I":
             self._ops[index] = op
         if not isinstance(coefficient, Number):
             raise ValueError("coefficient of PauliTerm must be a Number.")
         self.coefficient = complex(coefficient)
-        self._id = None
 
-    def id(self):
+    def id(self, sort_ops=True):
         """
-        Returns the unique identifier string for the PauliTerm (ignoring the coefficient).
-         Used in the simplify method of PauliSum.
+        Returns an identifier string for the PauliTerm (ignoring the coefficient).
+
+        Don't use this to compare terms. This function will not work with qubits that
+        aren't sortable.
 
-        :return: The unique identifier for this term.
+        :param sort_ops: Whether to sort operations by qubit. This is True by default for
+            backwards compatibility but will change in pyQuil 2.0. Callers should never rely
+            on comparing id's for testing equality. See ``operations_as_set`` instead.
+        :return: A string representation of this term's operations.
         :rtype: string
         """
-        if self._id is not None:
-            return self._id
+        if sort_ops and len(self._ops) > 1:
+            warnings.warn("`PauliTerm.id()` will not work on PauliTerms where the qubits are not "
+                          "sortable and should be avoided in favor of `operations_as_set`.",
+                          FutureWarning)
+            return ''.join("{}{}".format(self._ops[q], q) for q in sorted(self._ops.keys()))
         else:
-            s = ""
-            for index in sorted(self._ops.keys()):
-                s += "%s%s" % (self[index], index)
-            self._id = s
-            return s
+            return ''.join("{}{}".format(p, q) for q, p in self._ops.items())
+
+    def operations_as_set(self):
+        """
+        Return a frozenset of operations in this term.
+
+        Use this in place of :py:func:`id` if the order of operations in the term does not
+        matter.
+
+        :return: frozenset of strings representing Pauli operations
+        """
+        return frozenset(self._ops.items())
 
     def __eq__(self, other):
         if not isinstance(other, (PauliTerm, PauliSum)):
@@ -97,7 +117,15 @@ class PauliTerm(object):
         elif isinstance(other, PauliSum):
             return other == self
         else:
-            return self.id() == other.id() and np.isclose(self.coefficient, other.coefficient)
+            return (self.operations_as_set() == other.operations_as_set() and
+                    np.isclose(self.coefficient, other.coefficient))
+
+    def __hash__(self):
+        return hash((
+            round(self.coefficient.real * HASH_PRECISION),
+            round(self.coefficient.imag * HASH_PRECISION),
+            self.operations_as_set()
+        ))
 
     def __ne__(self, other):
         # x!=y and x<>y call __ne__() instead of negating __eq__
@@ -135,8 +163,7 @@ class PauliTerm(object):
     def get_qubits(self):
         """Gets all the qubits that this PauliTerm operates on.
         """
-        # sort the keys to get a deterministic iteration order over qubits
-        return sorted(self._ops.keys())
+        return list(self._ops.keys())
 
     def __getitem__(self, i):
         return self._ops.get(i, "I")
@@ -238,7 +265,7 @@ class PauliTerm(object):
         :rtype: PauliTerm
         """
         assert isinstance(other, Number)
-        return self + other
+        return PauliTerm("I", 0, other) + self
 
     def __sub__(self, other):
         """Subtracts a PauliTerm from this one.
@@ -260,7 +287,7 @@ class PauliTerm(object):
 
     def __str__(self):
         term_strs = []
-        for index in sorted(self._ops.keys()):
+        for index in self._ops.keys():
             term_strs.append("%s%s" % (self[index], index))
 
         if len(term_strs) == 0:
@@ -404,7 +431,8 @@ class PauliSum(object):
         elif len(self.terms) != len(other.terms):
             warnings.warn(UnequalLengthWarning("These PauliSums have a different number of terms."))
             return False
-        return all([term == other.terms[i] for i, term in enumerate(self.terms)])
+
+        return set(self.terms) == set(other.terms)
 
     def __ne__(self, other):
         """Inequality testing to see if two PauliSum's are not equivalent.
@@ -560,28 +588,7 @@ class PauliSum(object):
         """
         Simplifies the sum of Pauli operators according to Pauli algebra rules.
         """
-        def coalesce(d):
-            terms = []
-            for k in sorted(d):
-                term_list = d[k]
-                if (len(term_list) == 1 and not
-                        np.isclose(term_list[0].coefficient, 0.0)):
-                    terms.append(term_list[0])
-                else:
-                    coeff = sum(t.coefficient for t in term_list)
-                    if not np.isclose(coeff, 0.0):
-                        terms.append(term_with_coeff(term_list[0], coeff))
-            return PauliSum(terms)
-
-        like_terms = {}
-        for term in self.terms:
-            id = term.id()
-            if id not in like_terms:
-                like_terms[id] = [term]
-            else:
-                like_terms[id] = like_terms[id] + [term]
-
-        return coalesce(like_terms)
+        return simplify_pauli_sum(self)
 
     def get_programs(self):
         """
@@ -595,6 +602,37 @@ class PauliSum(object):
         return programs, coefficients
 
 
+def simplify_pauli_sum(pauli_sum):
+    # You might want to use a defaultdict(list) here, but don't because
+    # we want to do our best to preserve the order of terms.
+    like_terms = OrderedDict()
+    for term in pauli_sum.terms:
+        key = term.operations_as_set()
+        if key in like_terms:
+            like_terms[key].append(term)
+        else:
+            like_terms[key] = [term]
+
+    terms = []
+    for term_list in like_terms.values():
+        first_term = term_list[0]
+        if len(term_list) == 1 and not np.isclose(first_term.coefficient, 0.0):
+            terms.append(first_term)
+        else:
+            coeff = sum(t.coefficient for t in term_list)
+            for t in term_list:
+                if list(t._ops.items()) != list(first_term._ops.items()):
+                    warnings.warn("The term {} will be combined with {}, but they have different "
+                                  "orders of operations. This doesn't matter for QVM or "
+                                  "wavefunction simulation but may be important when "
+                                  "running on an actual device."
+                                  .format(t.id(sort_ops=False), first_term.id(sort_ops=False)))
+
+            if not np.isclose(coeff, 0.0):
+                terms.append(term_with_coeff(term_list[0], coeff))
+    return PauliSum(terms)
+
+
 def check_commutation(pauli_list, pauli_two):
     """
     Check if commuting a PauliTerm commutes with a list of other terms by natural calculation.
@@ -606,6 +644,7 @@ def check_commutation(pauli_list, pauli_two):
     :returns: True if pauli_two object commutes with pauli_list, False otherwise
     :rtype: bool
     """
+
     def coincident_parity(p1, p2):
         non_similar = 0
         p1_indices = set(p1._ops.keys())
@@ -730,6 +769,7 @@ def _exponentiate_general_case(pauli_term, param):
     :returns: A Quil program object
     :rtype: Program
     """
+
     def reverse_hack(p):
         # A hack to produce a *temporary* program which reverses p.
         revp = Program()
@@ -789,7 +829,7 @@ def suzuki_trotter(trotter_order, trotter_steps):
               type: o=0 is A and o=1 is B.
     :rtype: list
     """
-    p1 = p2 = p4 = p5 = 1.0 / (4 - (4**(1. / 3)))
+    p1 = p2 = p4 = p5 = 1.0 / (4 - (4 ** (1. / 3)))
     p3 = 1 - 4 * p1
     trotter_dict = {1: [(1, 0), (1, 1)],
                     2: [(0.5, 0), (1, 1), (0.5, 0)],
@@ -848,7 +888,7 @@ def trotterize(first_pauli_term, second_pauli_term, trotter_order=1,
     if not (1 <= trotter_order < 5):
         raise ValueError("trotterize only accepts trotter_order in {1, 2, 3, 4}.")
 
-    commutator = (first_pauli_term * second_pauli_term) +\
+    commutator = (first_pauli_term * second_pauli_term) + \
                  (-1 * second_pauli_term * first_pauli_term)
 
     prog = Program()


Test patch:
diff --git a/pyquil/tests/test_parametric.py b/pyquil/tests/test_parametric.py
index 0b7fc12..7ec1e81 100644
--- a/pyquil/tests/test_parametric.py
+++ b/pyquil/tests/test_parametric.py
@@ -163,4 +163,10 @@ def test_exponentiate_paraprog():
     xterm = PauliTerm("X", 2) * PauliTerm("X", 1)
     paraprog = exponential_map(xterm)
     prog = paraprog(1)
-    assert prog.out() == "H 1\nH 2\nCNOT 1 2\nRZ(2.0) 2\nCNOT 1 2\nH 1\nH 2\n"
+    assert prog.out() == ("H 2\n"
+                          "H 1\n"
+                          "CNOT 2 1\n"
+                          "RZ(2.0) 1\n"
+                          "CNOT 2 1\n"
+                          "H 2\n"
+                          "H 1\n")
diff --git a/pyquil/tests/test_paulis.py b/pyquil/tests/test_paulis.py
index f0c37e3..a0b2ef0 100644
--- a/pyquil/tests/test_paulis.py
+++ b/pyquil/tests/test_paulis.py
@@ -15,20 +15,21 @@
 #    limitations under the License.
 ##############################################################################
 
-import pytest
-from pyquil.paulis import (PauliTerm, PauliSum, exponential_map, exponentiate_commuting_pauli_sum,
-                           ID, UnequalLengthWarning, exponentiate, trotterize, is_zero,
-                           check_commutation, commuting_sets, term_with_coeff, sI, sX, sY,
-                           sZ, ZERO)
-from pyquil.quil import Program
-from pyquil.gates import RX, RZ, CNOT, H, X, PHASE
-
 import math
-from itertools import product
+import warnings
 from functools import reduce
+from itertools import product
 from operator import mul
-from six.moves import range
+
 import numpy as np
+import pytest
+from six.moves import range
+
+from pyquil.gates import RX, RZ, CNOT, H, X, PHASE
+from pyquil.paulis import PauliTerm, PauliSum, exponential_map, exponentiate_commuting_pauli_sum, \
+    ID, UnequalLengthWarning, exponentiate, trotterize, is_zero, check_commutation, commuting_sets, \
+    term_with_coeff, sI, sX, sY, sZ, ZERO, is_identity
+from pyquil.quil import Program
 
 
 def isclose(a, b, rel_tol=1e-10, abs_tol=0.0):
@@ -93,9 +94,9 @@ def test_simplify_term_xz():
 
 
 def test_simplify_term_multindex():
-    term = PauliTerm('X', 0, coefficient=-0.5) * PauliTerm('Z', 0, coefficient=-1.0) \
-        * PauliTerm('X', 2, 0.5)
-    assert term.id() == 'Y0X2'
+    term = (PauliTerm('X', 0, coefficient=-0.5) *
+            PauliTerm('Z', 0, coefficient=-1.0) * PauliTerm('X', 2, 0.5))
+    assert term.id(sort_ops=False) == 'Y0X2'
     assert term.coefficient == -0.25j
 
 
@@ -160,7 +161,22 @@ def test_getitem():
 def test_ids():
     term_1 = PauliTerm("Z", 0, 1.0) * PauliTerm("Z", 1, 1.0) * PauliTerm("X", 5, 5)
     term_2 = PauliTerm("X", 5, 5) * PauliTerm("Z", 0, 1.0) * PauliTerm("Z", 1, 1.0)
-    assert term_1.id() == term_2.id()
+    with pytest.warns(FutureWarning) as w:
+        assert term_1.id() == term_2.id()
+    assert 'should be avoided' in str(w[0])
+
+
+def test_ids_no_sort():
+    term_1 = PauliTerm("Z", 0, 1.0) * PauliTerm("Z", 1, 1.0) * PauliTerm("X", 5, 5)
+    term_2 = PauliTerm("X", 5, 5) * PauliTerm("Z", 0, 1.0) * PauliTerm("Z", 1, 1.0)
+    assert term_1.id(sort_ops=False) == 'Z0Z1X5'
+    assert term_2.id(sort_ops=False) == 'X5Z0Z1'
+
+
+def test_operations_as_set():
+    term_1 = PauliTerm("Z", 0, 1.0) * PauliTerm("Z", 1, 1.0) * PauliTerm("X", 5, 5)
+    term_2 = PauliTerm("X", 5, 5) * PauliTerm("Z", 0, 1.0) * PauliTerm("Z", 1, 1.0)
+    assert term_1.operations_as_set() == term_2.operations_as_set()
 
 
 def test_pauliop_inputs():
@@ -214,7 +230,7 @@ def test_ps_adds_pt_2():
     assert str(b + 1.0) == "(3+0j)*I"
     assert str(1.0 + b) == "(3+0j)*I"
     b = sX(0) + 1.0
-    assert str(b) == "(1+0j)*I + (1+0j)*X0"
+    assert str(b) == "(1+0j)*X0 + (1+0j)*I"
     b = 1.0 + sX(0)
     assert str(b) == "(1+0j)*I + (1+0j)*X0"
 
@@ -233,12 +249,11 @@ def test_ps_sub():
     b = 1.0 - sX(0)
     assert str(b) == "(1+0j)*I + (-1+0j)*X0"
     b = sX(0) - 1.0
-    assert str(b) == "(-1+0j)*I + (1+0j)*X0"
+    assert str(b) == "(1+0j)*X0 + (-1+0j)*I"
 
 
 def test_zero_terms():
-    term = PauliTerm("X", 0, 1.0) + PauliTerm("X", 0, -1.0) + \
-        PauliTerm("Y", 0, 0.5)
+    term = PauliTerm("X", 0, 1.0) + PauliTerm("X", 0, -1.0) + PauliTerm("Y", 0, 0.5)
     assert str(term) == "(0.5+0j)*Y0"
 
     term = PauliTerm("X", 0, 1.0) + PauliTerm("X", 0, -1.0)
@@ -254,14 +269,13 @@ def test_zero_terms():
     term4 = PauliSum([])
     assert str(term4) == "0j*I"
 
-    term = PauliSum([PauliTerm("X", 0, 0.0), PauliTerm("Y", 1, 1.0) *
-                     PauliTerm("Z", 2)])
+    term = PauliSum([PauliTerm("X", 0, 0.0), PauliTerm("Y", 1, 1.0) * PauliTerm("Z", 2)])
     assert str(term) == "0j*X0 + (1+0j)*Y1*Z2"
     term = term.simplify()
     assert str(term) == "(1+0j)*Y1*Z2"
 
 
-def test_exponentiate():
+def test_exponentiate_1():
     # test rotation of single qubit
     generator = PauliTerm("Z", 0, 1.0)
     para_prog = exponential_map(generator)
@@ -269,47 +283,57 @@ def test_exponentiate():
     result_prog = Program().inst(RZ(2.0)(0))
     assert prog == result_prog
 
+
+def test_exponentiate_2():
     # testing general 2-circuit
-    generator = PauliTerm("Z", 1, 1.0) * PauliTerm("Z", 0, 1.0)
+    generator = PauliTerm("Z", 0, 1.0) * PauliTerm("Z", 1, 1.0)
     para_prog = exponential_map(generator)
     prog = para_prog(1)
     result_prog = Program().inst(CNOT(0, 1)).inst(RZ(2.0)(1)).inst(CNOT(0, 1))
     assert prog == result_prog
 
+
+def test_exponentiate_bp0_ZX():
     # testing change of basis position 0
-    generator = PauliTerm("Z", 1, 1.0) * PauliTerm("X", 0, 1.0)
+    generator = PauliTerm("X", 0, 1.0) * PauliTerm("Z", 1, 1.0)
     param_prog = exponential_map(generator)
     prog = param_prog(1)
-    result_prog = Program().inst([H(0), CNOT(0, 1), RZ(2.0)(1), CNOT(0, 1),
-                                  H(0)])
+    result_prog = Program().inst([H(0), CNOT(0, 1), RZ(2.0)(1), CNOT(0, 1), H(0)])
     assert prog == result_prog
 
+
+def test_exponentiate_bp1_XZ():
     # testing change of basis position 1
-    generator = PauliTerm("X", 1, 1.0) * PauliTerm("Z", 0, 1.0)
+    generator = PauliTerm("Z", 0, 1.0) * PauliTerm("X", 1, 1.0)
     para_prog = exponential_map(generator)
     prog = para_prog(1)
-    result_prog = Program().inst([H(1), CNOT(0, 1), RZ(2.0)(1), CNOT(0, 1),
-                                  H(1)])
+    result_prog = Program().inst([H(1), CNOT(0, 1), RZ(2.0)(1), CNOT(0, 1), H(1)])
     assert prog == result_prog
 
+
+def test_exponentiate_bp0_ZY():
     # testing change of basis position 0
-    generator = PauliTerm("Z", 1, 1.0) * PauliTerm("Y", 0, 1.0)
+    generator = PauliTerm("Y", 0, 1.0) * PauliTerm("Z", 1, 1.0)
     para_prog = exponential_map(generator)
     prog = para_prog(1)
     result_prog = Program().inst([RX(math.pi / 2.0)(0), CNOT(0, 1), RZ(2.0)(1),
                                   CNOT(0, 1), RX(-math.pi / 2)(0)])
     assert prog == result_prog
 
+
+def test_exponentiate_bp1_YZ():
     # testing change of basis position 1
-    generator = PauliTerm("Y", 1, 1.0) * PauliTerm("Z", 0, 1.0)
+    generator = PauliTerm("Z", 0, 1.0) * PauliTerm("Y", 1, 1.0)
     para_prog = exponential_map(generator)
     prog = para_prog(1)
-    result_prog = Program().inst([RX(math.pi / 2.0)(1), CNOT(0, 1), RZ(2.0)(1),
-                                  CNOT(0, 1), RX(-math.pi / 2.0)(1)])
+    result_prog = Program().inst([RX(math.pi / 2.0)(1), CNOT(0, 1),
+                                  RZ(2.0)(1), CNOT(0, 1), RX(-math.pi / 2.0)(1)])
     assert prog == result_prog
 
+
+def test_exponentiate_3cob():
     # testing circuit for 3-terms with change of basis
-    generator = PauliTerm("X", 2, 1.0) * PauliTerm("Y", 1, 1.0) * PauliTerm("Z", 0, 1.0)
+    generator = PauliTerm("Z", 0, 1.0) * PauliTerm("Y", 1, 1.0) * PauliTerm("X", 2, 1.0)
     para_prog = exponential_map(generator)
     prog = para_prog(1)
     result_prog = Program().inst([RX(math.pi / 2.0)(1), H(2), CNOT(0, 1),
@@ -317,10 +341,13 @@ def test_exponentiate():
                                   CNOT(0, 1), RX(-math.pi / 2.0)(1), H(2)])
     assert prog == result_prog
 
+
+def test_exponentiate_3ns():
     # testing circuit for 3-terms non-sequential
-    generator = PauliTerm("Y", 3, 1.0) * PauliTerm("Y", 2, 1.0) * PauliTerm("I", 1,
-                                                                            1.0) * PauliTerm("Y", 0,
-                                                                                             1.0)
+    generator = (PauliTerm("Y", 0, 1.0) *
+                 PauliTerm("I", 1, 1.0) *
+                 PauliTerm("Y", 2, 1.0) *
+                 PauliTerm("Y", 3, 1.0))
     para_prog = exponential_map(generator)
     prog = para_prog(1)
     result_prog = Program().inst([RX(math.pi / 2.0)(0), RX(math.pi / 2.0)(2),
@@ -433,23 +460,31 @@ def test_check_commutation():
     assert check_commutation([term2], term3)
     assert not check_commutation([term1], term3)
 
+
+def _commutator(t1, t2):
+    with warnings.catch_warnings():
+        warnings.filterwarnings('ignore',
+                                message=r"The term .+ will be combined with .+, "
+                                        r"but they have different orders of operations.*",
+                                category=UserWarning)
+        return t1 * t2 + -1 * t2 * t1
+
+
+def test_check_commutation_rigorous():
     # more rigorous test.  Get all operators in Pauli group
     p_n_group = ("I", "X", "Y", "Z")
     pauli_list = list(product(p_n_group, repeat=3))
     pauli_ops = [list(zip(x, range(3))) for x in pauli_list]
     pauli_ops_pq = [reduce(mul, (PauliTerm(*x) for x in op)) for op in pauli_ops]
 
-    def commutator(t1, t2):
-        return t1 * t2 + -1 * t2 * t1
-
     non_commuting_pairs = []
     commuting_pairs = []
     for x in range(len(pauli_ops_pq)):
         for y in range(x, len(pauli_ops_pq)):
 
-            tmp_op = commutator(pauli_ops_pq[x], pauli_ops_pq[y])
+            tmp_op = _commutator(pauli_ops_pq[x], pauli_ops_pq[y])
             assert len(tmp_op.terms) == 1
-            if tmp_op.terms[0].id() == '':
+            if is_identity(tmp_op.terms[0]):
                 commuting_pairs.append((pauli_ops_pq[x], pauli_ops_pq[y]))
             else:
                 non_commuting_pairs.append((pauli_ops_pq[x], pauli_ops_pq[y]))
@@ -563,8 +598,40 @@ def test_from_list():
         pterm = PauliTerm.from_list([("X", 0), ("Y", 0)])
 
 
+def test_ordered():
+    term = sZ(3) * sZ(2) * sZ(1)
+    prog = exponential_map(term)(0.5)
+    assert prog.out() == "CNOT 3 2\n" \
+                         "CNOT 2 1\n" \
+                         "RZ(1.0) 1\n" \
+                         "CNOT 2 1\n" \
+                         "CNOT 3 2\n"
+
+
 def test_numpy_integer_types():
     idx_np, = np.arange(1, dtype=np.int64)
     assert isinstance(idx_np, np.int64)
     # on python 3 this fails unless explicitly allowing for numpy integer types
     PauliTerm("X", idx_np)
+
+
+def test_simplify():
+    t1 = sZ(0) * sZ(1)
+    t2 = sZ(0) * sZ(1)
+    assert (t1 + t2) == 2 * sZ(0) * sZ(1)
+
+
+def test_dont_simplify():
+    t1 = sZ(0) * sZ(1)
+    t2 = sZ(2) * sZ(3)
+    assert (t1 + t2) != 2 * sZ(0) * sZ(1)
+
+
+def test_simplify_warning():
+    t1 = sZ(0) * sZ(1)
+    t2 = sZ(1) * sZ(0)
+    with pytest.warns(UserWarning) as e:
+        tsum = t1 + t2
+
+    assert tsum == 2 * sZ(0) * sZ(1)
+    assert str(e[0].message).startswith('The term Z1Z0 will be combined with Z0Z1')