From db6f3dcc8b612f7c86d95c4818fb8d0a996a204b Mon Sep 17 00:00:00 2001
From: Pete Shadbolt <pete.shadbolt@gmail.com>
Date: Tue, 10 May 2016 22:00:00 +0100
Subject: [PATCH] Adding a circuit model simulator for testing

---
 abp/graphstate.py              |  4 +++
 abp/qi.py                      | 48 +++++++++++++++++++++++++++++++++-
 tests/test_measurement.py      |  6 +++++
 tests/test_qi_circuit_model.py | 22 ++++++++++++++++
 4 files changed, 79 insertions(+), 1 deletion(-)
 create mode 100644 tests/test_measurement.py
 create mode 100644 tests/test_qi_circuit_model.py

diff --git a/abp/graphstate.py b/abp/graphstate.py
index f50ef02..44f8b14 100644
--- a/abp/graphstate.py
+++ b/abp/graphstate.py
@@ -106,17 +106,21 @@ class GraphState(object):
     def measure_z(self, node, force = None):
         """ Measure the graph in the Z-basis """
         res = force if force else np.random.choice([0,1])
+
+        # Disconnect
         for neighbour in self.ngbh[node]:
             self.del_edge(node, neighbour)
             if res:
                 self.act_local_rotation_by_name(neighbour, "pz")
 
+        # Set final state as appropriate
         if res:
             self.act_local_rotation_by_name(node, "px")
             self.act_local_rotation_by_name(node, "hadamard")
         else:
             self.act_local_rotation_by_name(node, "hadamard")
 
+        return res
 
 
     def measure_x(self, i):
diff --git a/abp/qi.py b/abp/qi.py
index 3c3f097..14f64b1 100644
--- a/abp/qi.py
+++ b/abp/qi.py
@@ -7,17 +7,20 @@ Exposes a few basic QI operators
 
 import numpy as np
 from scipy.linalg import sqrtm
+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)
 # Operators
 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) / np.sqrt(2)
+ha = 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)
 
@@ -44,3 +47,46 @@ names = "identity", "px", "py", "pz", "hadamard", "phase", "sqz", "msqz", "sqy",
 by_name = dict(zip(names, common_us))
 
 paulis = px, py, pz
+
+class CircuitModel(object):
+    def __init__(self, nqubits):
+        self.nqubits = nqubits
+        self.d = 2**nqubits
+        self.state = np.zeros((self.d, 1), dtype=complex)
+        self.state[0, 0]=1
+
+    def act_cz(self, control, target):
+        """ Act a CU somewhere """
+        control = 1 << control
+        target = 1 << control
+        print control, target
+        for i in xrange(self.d):
+            if (i & control) and (i & target):
+                self.state[i, 0] *= -1
+
+    def act_hadamard(self, qubit):
+        """ Act a hadamard somewhere """
+        where = 1 << qubit
+        output = np.zeros((self.d, 1), dtype=complex)
+        for i, v in enumerate(self.state):
+            if (i & where) == 0:
+                output[i] += v*ir2
+                output[i ^ where] += v*ir2
+            if (i & where) == 1:
+                output[i] += v*ir2
+                output[i ^ where] -= v*ir2
+        self.state = output
+
+
+    def local_rotation(self, qubit, u):
+        """ Act a local unitary somwhere """
+        pass
+
+    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:
+                s += "|{}>: {}\n".format(label, self.state[i, 0])
+        return s
+
diff --git a/tests/test_measurement.py b/tests/test_measurement.py
new file mode 100644
index 0000000..b572a2f
--- /dev/null
+++ b/tests/test_measurement.py
@@ -0,0 +1,6 @@
+from abp.graphstate import GraphState
+
+def test_z_measurement():
+    #TODO
+    pass
+
diff --git a/tests/test_qi_circuit_model.py b/tests/test_qi_circuit_model.py
new file mode 100644
index 0000000..1371e0b
--- /dev/null
+++ b/tests/test_qi_circuit_model.py
@@ -0,0 +1,22 @@
+from abp import qi
+
+def test_init():
+    """ Can you initialize some qubits """
+    psi = qi.CircuitModel(5)
+    assert psi.d == 32
+
+def test_cz():
+    """ What does CZ do ? """
+    psi = qi.CircuitModel(2)
+    #psi.act_hadamard(0)
+    psi.act_hadamard(0)
+    print psi
+    psi.act_hadamard(1)
+    print psi
+    psi.act_cz(0, 1)
+    print psi
+    psi.act_cz(0, 1)
+    print psi
+    #psi.act_cz(0, 1)
+
+