Add state equality checking, global phase fixing.

abp/clifford.py

@@ -26,26 +26,17 @@ def get_name(i):
 
 def find_clifford(needle, haystack):
     """ Find the index of a given u within a list of unitaries, up to a global phase  """
-    needle = normalize_global_phase(needle)
+    needle = qi.normalize_global_phase(needle)
     for i, t in enumerate(haystack):
         if np.allclose(t, needle):
             return i
     raise IndexError
 
-
-def normalize_global_phase(m):
-    """ Normalize the global phase of a matrix """
-    v = (x for x in m.flatten() if np.abs(x) > 0.001).next()
-    phase = np.arctan2(v.imag, v.real) % np.pi
-    rot = np.exp(-1j * phase)
-    return rot * m if rot * v > 0 else -rot * m
-
-
 def find_cz(bond, c1, c2, commuters, state_table, ab_cz_table):
     """ Find the output of a CZ operation """
     # Figure out the target state
     target = qi.cz.dot(state_table[bond, c1, c2])
-    target = normalize_global_phase(target)
+    target = qi.normalize_global_phase(target)
 
     # Choose the sets to search over
     s1 = commuters if c1 in commuters else xrange(24)
@@ -64,7 +55,7 @@ def compose_u(decomposition):
     """ Get the unitary representation of a particular decomposition """
     matrices = ({"x": qi.sqx, "z": qi.msqz}[c] for c in decomposition)
     output = reduce(np.dot, matrices, np.eye(2, dtype=complex))
-    return normalize_global_phase(output)
+    return qi.normalize_global_phase(output)
 
 
 def get_unitaries():
@@ -96,7 +87,7 @@ def get_state_table(unitaries):
     for bond, i, j in tqdm(params, desc="Building state table"):
         state = qi.bond if bond else qi.nobond
         kp = np.kron(unitaries[i], unitaries[j])
-        state_table[bond, i, j, :] = normalize_global_phase(
+        state_table[bond, i, j, :] = qi.normalize_global_phase(
             np.dot(kp, state).T)
     return state_table
 

abp/qi.py

@@ -3,6 +3,7 @@
 
 """
 Exposes a few basic QI operators
+And a circuit-model simulator
 """
 
 import numpy as np
@@ -20,7 +21,7 @@ id = np.array(np.eye(2, dtype=complex))
 px = np.array([[0, 1], [1, 0]], dtype=complex)
 py = np.array([[0, -1j], [1j, 0]], dtype=complex)
 pz = np.array([[1, 0], [0, -1]], dtype=complex)
-ha = np.array([[1, 1], [1, -1]], dtype=complex) * ir2
+ha = hadamard = np.array([[1, 1], [1, -1]], dtype=complex) * ir2
 ph = np.array([[1, 0], [0, 1j]], dtype=complex)
 t = np.array([[1, 0], [0, np.exp(1j*np.pi/4)]], dtype=complex)
 
@@ -48,6 +49,15 @@ by_name = dict(zip(names, common_us))
 
 paulis = px, py, pz
 
+
+def normalize_global_phase(m):
+    """ Normalize the global phase of a matrix """
+    v = (x for x in m.flatten() if np.abs(x) > 0.001).next()
+    phase = np.arctan2(v.imag, v.real) % np.pi
+    rot = np.exp(-1j * phase)
+    return rot * m if rot * v > 0 else -rot * m
+
+
 class CircuitModel(object):
     def __init__(self, nqubits):
         self.nqubits = nqubits
@@ -84,6 +94,13 @@ class CircuitModel(object):
             output[i ^ where] += v*u[not q, q]
         self.state = output
 
+    def __eq__(self, other):
+        """ Check whether two quantum states are the same or not
+            UP TO A GLOBAL PHASE """
+        a = normalize_global_phase(self.state)
+        b = normalize_global_phase(other.state)
+        return np.allclose(a, b)
+
 
     def __str__(self):
         s = ""

tests/test_against_circuit_model.py

@@ -14,7 +14,7 @@ def test_hadamard_only_multiqubit():
         g.act_hadamard(i)
         c.act_hadamard(i)
 
-    assert np.allclose(g.to_state_vector().state, c.state)
+    assert g.to_state_vector() == c
 
     for i in range(100):
         a, b = np.random.randint(0, n-1, 2)
@@ -22,9 +22,7 @@ def test_hadamard_only_multiqubit():
             g.act_cz(a, b)
             c.act_cz(a, b)
 
-    s1 = clifford.normalize_global_phase(g.to_state_vector().state)
-    s2 = clifford.normalize_global_phase(c.state)
-    assert np.allclose(s1, s2)
+    assert g.to_state_vector() == c
 
 
 def test_all_multiqubit():
@@ -34,13 +32,13 @@ def test_all_multiqubit():
     c = CircuitModel(n)
 
     for i in range(10):
-        i = np.random.randint(0, n-1)
-        j = np.random.randint(0, 24)
-        print i, j
-        g.act_local_rotation(i, j)
-        c.act_local_rotation(i, clifford.unitaries[j])
+        qubit = np.random.randint(0, n-1)
+        rotation = np.random.randint(0, 24-1)
+        g.act_local_rotation(qubit, rotation)
+        c.act_local_rotation(qubit, clifford.unitaries[rotation])
 
-    assert np.allclose(g.to_state_vector().state, c.state)
+
+    assert g.to_state_vector() == c
 
     #for i in range(100):
         #a, b = np.random.randint(0, n-1, 2)
@@ -48,7 +46,5 @@ def test_all_multiqubit():
             #g.act_cz(a, b)
             #c.act_cz(a, b)
 
-    #s1 = clifford.normalize_global_phase(g.to_state_vector().state)
-    #s2 = clifford.normalize_global_phase(c.state)
-    #assert np.allclose(s1, s2)
+    assert g.to_state_vector() == c
 

tests/test_cphase_against_anders.py

@@ -30,7 +30,7 @@ def test_cz_table():
 
         # Go and compute the output
         computed_output = np.dot(qi.cz, input_state)
-        computed_output = clifford.normalize_global_phase(computed_output)
+        computed_output = qi.normalize_global_phase(computed_output)
 
         # Now look up the answer in the table
         bondp, c1p, c2p = ab_cz_table[bond, c1, c2]

tests/test_cphase_table.py

@@ -15,7 +15,7 @@ def test_cz_table():
 
         # Go and compute the output
         computed_output = np.dot(qi.cz, input_state)
-        computed_output = clifford.normalize_global_phase(computed_output)
+        computed_output = qi.normalize_global_phase(computed_output)
 
         # Now look up the answer in the table
         bondp, c1p, c2p = clifford.cz_table[bond, c1, c2]

tests/test_normalize_global_phase.py

@@ -7,5 +7,5 @@ def test_normalize_global_phase():
         u = qi.pz
         phase = np.random.uniform(0, 2*np.pi)
         m = np.exp(1j*phase) * u
-        normalized = clifford.normalize_global_phase(m)
+        normalized = qi.normalize_global_phase(m)
         assert np.allclose(normalized, u)

tests/test_qi_circuit_model.py

@@ -1,11 +1,38 @@
 import numpy as np
 from abp import qi
 
-def _test_init():
+def test_init():
     """ Can you initialize some qubits """
     psi = qi.CircuitModel(5)
     assert psi.d == 32
 
+def test_single_qubit_stuff():
+    """ Try some sensible single-qubit things """
+    psi = qi.CircuitModel(2)
+    psi.act_local_rotation(0, qi.px)
+    assert np.allclose(psi.state[1], 1)
+    psi.act_local_rotation(0, qi.px)
+    assert np.allclose(psi.state[0], 1)
+    psi.act_local_rotation(0, qi.pz)
+
+
+def test_more_single_qubit_stuff():
+    """ Try some sensible single-qubit things """
+    psi = qi.CircuitModel(2)
+    psi.act_local_rotation(0, qi.px)
+    psi.act_local_rotation(1, qi.px)
+    psi.act_cz(0, 1)
+
+def test_equality():
+    """ Test that equality succeeds / fails as desired """
+    a = qi.CircuitModel(2)
+    b = qi.CircuitModel(2)
+    assert a == b
+    a.act_local_rotation(0, qi.px)
+    assert a != b
+
+
+
 def test_hadamard():
     """ What does CZ do ? """
     psi = qi.CircuitModel(3)