Not there yet

abp/graphstate.py

@@ -118,20 +118,27 @@ class GraphState(object):
     def measure(self, node, basis, force=None):
         """ Measure in an arbitrary basis """
         basis = clifford.by_name[basis]
-        old_basis = basis
         ha = clifford.conjugation_table[self.node[node]["vop"]]
         basis, phase = clifford.conjugate(basis, ha)
-        assert phase in (-1, 1)  # TODO: remove
 
-        # TODO: wtf
-        force = force ^ 0x01 if force != -1 and phase == 0 else force
+        #TODO: wtf
+        #force = force ^ 0x01 if force != -1 and phase == 0 else force
+        if force != None and phase == 0+0j: 
+            force = not force
+
+        result = random.choice([0, 1])
+        if which == clifford.by_name["px"]:
+            result = self.measure_x(node, result)
+        elif which == clifford.by_name["py"]:
+            result = self.measure_y(node, result)
+        elif which == clifford.by_name["pz"]:
+            result = self.measure_z(node, result)
+        else:
+            raise ValueError("You can only measure in PX, PY or PZ")
 
-        which = {1: self.measure_x, 2:
-                 self.measure_y, 3: self.measure_z}[basis]
-        res = which(node, force)
-        res = res if phase == 1 else not res
+        if phase == 1:
+            result = not result
 
-        # TODO: put the asserts from graphsim.cpp into tests
         return res
 
     def toggle_edges(a, b):
@@ -142,14 +149,11 @@ class GraphState(object):
                 done.add((i, j), (j, i))
                 self.toggle_edge(i, j)
 
-    def measure_x(self, node, force=None):
+    def measure_x(self, node, result):
         """ Measure the graph in the X-basis """
         if len(self.adj[node]) == 0:
             return 0
 
-        # Flip a coin
-        result = force if force != None else random.choice([0, 1])
-
         # Pick a vertex
         friend = next(self.adj[node].iterkeys())
 
@@ -179,11 +183,8 @@ class GraphState(object):
 
         return result
 
-    def measure_y(self, node, force=None):
+    def measure_y(self, node, result):
         """ Measure the graph in the Y-basis """
-        # Flip a coin
-        result = force if force != None else random.choice([0, 1])
-
         # Do some rotations
         for neighbour in self.adj[node]:
             # NB: should these be hermitian_conjugated?
@@ -197,11 +198,8 @@ class GraphState(object):
         self.act_local_rotation(node, "msqz" if result else "msqz_h")
         return result
 
-    def measure_z(self, node, force=None):
+    def measure_z(self, node, result):
         """ Measure the graph in the Z-basis """
-        # Flip a coin
-        result = force if force != None else random.choice([0, 1])
-
         # Disconnect
         for neighbour in self.adj[node]:
             self.del_edge(node, neighbour)
@@ -268,8 +266,10 @@ class GraphState(object):
                 output[a][b] = "Z"
             else:
                 output[a][b] = "I"
+        # TODO: signs
         return output
 
     def __eq__(self, other):
         """ Check equality between graphs """
         return self.adj == other.adj and self.node == other.node
+

abp/qi.py

@@ -9,12 +9,13 @@ And a circuit-model simulator
 import numpy as np
 import itertools as it
 
+
 def hermitian_conjugate(u):
     """ Shortcut to the Hermitian conjugate """
     return np.conjugate(np.transpose(u))
 
-# Constants 
-ir2 = 1/np.sqrt(2)
+# Constants
+ir2 = 1 / np.sqrt(2)
 # Operators
 id = np.array(np.eye(2, dtype=complex))
 px = np.array([[0, 1], [1, 0]], dtype=complex)
@@ -22,23 +23,29 @@ py = np.array([[0, -1j], [1j, 0]], dtype=complex)
 pz = np.array([[1, 0], [0, -1]], dtype=complex)
 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)
-
-sqx = np.array([[ 1.+0.j, -0.+1.j], [-0.+1.j, 1.-0.j]], dtype=complex)*ir2
-msqx = np.array([[ 1.+0.j, 0.-1.j], [ 0.-1.j, 1.-0.j]], dtype=complex)*ir2
-sqy = np.array([[ 1.+0.j, 1.+0.j], [-1.-0.j, 1.-0.j]], dtype=complex)*ir2
-msqy = np.array([[ 1.+0.j, -1.-0.j], [ 1.+0.j, 1.-0.j]], dtype=complex)*ir2
-sqz = np.array([[ 1.+1.j, 0.+0.j], [ 0.+0.j, 1.-1.j]], dtype=complex)*ir2
-msqz = np.array([[ 1.-1.j, 0.+0.j], [ 0.+0.j, 1.+1.j]], dtype=complex)*ir2
+t = np.array([[1, 0], [0, np.exp(1j * np.pi / 4)]], dtype=complex)
+
+sqx = np.array(
+    [[1. + 0.j, -0. + 1.j], [-0. + 1.j, 1. - 0.j]], dtype=complex) * ir2
+msqx = np.array(
+    [[1. + 0.j, 0. - 1.j], [0. - 1.j, 1. - 0.j]], dtype=complex) * ir2
+sqy = np.array(
+    [[1. + 0.j, 1. + 0.j], [-1. - 0.j, 1. - 0.j]], dtype=complex) * ir2
+msqy = np.array(
+    [[1. + 0.j, -1. - 0.j], [1. + 0.j, 1. - 0.j]], dtype=complex) * ir2
+sqz = np.array(
+    [[1. + 1.j, 0. + 0.j], [0. + 0.j, 1. - 1.j]], dtype=complex) * ir2
+msqz = np.array(
+    [[1. - 1.j, 0. + 0.j], [0. + 0.j, 1. + 1.j]], dtype=complex) * ir2
 
 # CZ gate
 cz = np.array(np.eye(4), dtype=complex)
-cz[3,3]=-1
+cz[3, 3] = -1
 
 # States
-zero = np.array([[1],[0]], dtype=complex)
-one = np.array([[0],[1]], dtype=complex)
-plus = np.array([[1],[1]], dtype=complex) / np.sqrt(2)
+zero = np.array([[1], [0]], dtype=complex)
+one = np.array([[0], [1]], dtype=complex)
+plus = np.array([[1], [1]], dtype=complex) / np.sqrt(2)
 bond = cz.dot(np.kron(plus, plus))
 nobond = np.kron(plus, plus)
 
@@ -60,11 +67,12 @@ def normalize_global_phase(m):
 
 
 class CircuitModel(object):
+
     def __init__(self, nqubits):
         self.nqubits = nqubits
-        self.d = 2**nqubits
+        self.d = 2 ** nqubits
         self.state = np.zeros((self.d, 1), dtype=complex)
-        self.state[0, 0]=1
+        self.state[0, 0] = 1
 
     def act_cz(self, control, target):
         """ Act a CU somewhere """
@@ -86,19 +94,18 @@ class CircuitModel(object):
         output = np.zeros((self.d, 1), dtype=complex)
         for i, v in enumerate(self.state):
             q = i & where > 0
-            output[i] += v*ha[q, q]
-            output[i ^ where] += v*ha[not q, q]
+            output[i] += v * ha[q, q]
+            output[i ^ where] += v * ha[not q, q]
         self.state = output
 
-
     def act_local_rotation(self, qubit, u):
         """ Act a local unitary somwhere """
         where = 1 << qubit
         output = np.zeros((self.d, 1), dtype=complex)
         for i, v in enumerate(self.state):
             q = i & where > 0
-            output[i] += v*u[q, q] # TODO this is probably wrong
-            output[i ^ where] += v*u[not q, q]
+            output[i] += v * u[q, q]  # TODO this is probably wrong
+            output[i ^ where] += v * u[not q, q]
         self.state = output
 
     def __eq__(self, other):
@@ -108,12 +115,10 @@ class CircuitModel(object):
         b = normalize_global_phase(other.state)
         return np.allclose(a, b)
 
-
     def __str__(self):
         s = ""
         for i in range(self.d):
             label = bin(i)[2:].rjust(self.nqubits, "0")
-            if abs(self.state[i, 0])>0.00001:
+            if abs(self.state[i, 0]) > 0.00001:
                 s += "|{}>: {}\n".format(label, self.state[i, 0].round(3))
         return s
-

tests/test_measurement.py

@@ -1,7 +1,10 @@
+import numpy as np
 from abp import GraphState
+from abp import qi
 from anders_briegel import graphsim
 
-def test_measurements():
+def test_single_qubit_measurements():
+    """ Various simple tests of measurements """
 
     # Test that measuring |0> in Z gives 0
     g = GraphState([0])
@@ -20,7 +23,9 @@ def test_measurements():
     g.act_local_rotation(0, "pz")
     assert g.measure(0, "px") == 1, "Measuring |-> in X gives 1"
 
-    # Test random outcomes
+
+def test_random_outcomes():
+    """ Testing random behaviour """
     ones = 0
     for i in range(1000):
         g = GraphState([0])
@@ -28,6 +33,22 @@ def test_measurements():
         ones += g.measure(0, "pz")
     assert 400 < ones < 600, "This is a probabilistic test!"
 
+def test_projection():
+    """ Test that projection works correctly """
+    g = GraphState([0])
+    g.act_local_rotation(0, "hadamard")
+    g.measure(0, "pz", 0)
+    print g.to_state_vector()
+    assert np.allclose(g.to_state_vector().state, qi.zero)
+
+    g = GraphState([0])
+    g.act_local_rotation(0, "hadamard")
+    print g.to_state_vector()
+    g.measure(0, "pz", 1)
+    print g.to_state_vector()
+    assert np.allclose(g.to_state_vector().state, qi.one)
+
+
 
 
 def test_z_measurement_against_ab():