Added an option for non unitary gates

the reasons for this are just for experimentation with circuit partitioning

benchmark_csv.py

@@ -1,56 +0,0 @@
-# This script is to run the benchmarks.
-import numpy as np
-import sys
-import csv
-import time
-import qcgpu
-import os.path
-
-results_file = 'benchmark_results.csv'
-
-def bench_qcgpu(n, depth):
-    state = qcgpu.State(n)
-
-    h = qcgpu.gate.h()
-    x = qcgpu.gate.x()
-    sqrt_x = qcgpu.gate.sqrt_x()
-
-    start = time.time()
-
-    for level in range(depth):
-        for q in range(n):
-    
-            state.apply_gate(h, q)
-            state.apply_gate(sqrt_x, q)
-
-            if q != 0:
-                state.apply_controlled_gate(x, q, 0)
-        
-    runtime = time.time() - start
-    return {'name': 'qcgpu', 'num_qubits': n, 'depth': depth, 'time': runtime}
-
-def create_csv(filename):
-    file_exists = os.path.isfile(filename)
-    csvfile = open(filename, 'a')
-   
-    headers = ['name', 'num_qubits', 'depth', 'time']
-    writer = csv.DictWriter(csvfile, delimiter=',', lineterminator='\n',fieldnames=headers)
-
-    if not file_exists:
-        writer.writeheader()  # file doesn't exist yet, write a header
-
-    # writer.writerow({'TimeStamp': dic['ts'], 'light': dic['light'], 'Proximity': dic['prox']}
-
-    return writer
-
-def write_csv(writer, data):
-    print("Qubits: " + str(data['num_qubits']) + ", Depth: " + str(data['depth']) + ", Time: " + str(data['depth']))
-    writer.writerow(data)
-
-writer = create_csv(results_file)
-for i in range(2, 25):
-    for d in [5,10,15,20]:
-        data = bench_qcgpu(i, d)
-        write_csv(writer, data)
-
-

examples/deutsch-jozsa.py

@@ -0,0 +1,41 @@
+import qcgpu
+
+# 3 qubits, f(x) = x_0 NOT x_1 x_2
+# Balanced
+balanced_state = qcgpu.State(3)
+
+balanced_state.apply_all(qcgpu.gate.h())
+
+# Oracle U_f
+balanced_state.h(2)
+balanced_state.z(0)
+balanced_state.cx(1, 2)
+balanced_state.h(2)
+
+balanced_state.apply_all(qcgpu.gate.h())
+
+outcomes = balanced_state.measure(samples = 1000)
+
+if int(max(outcomes, key=outcomes.get)) == 0:
+    print('constant')
+else:
+    print('balanced')
+
+
+# 3 qubits, f(x) = 0
+# Constant
+constant_state = qcgpu.State(3)
+
+constant_state.apply_all(qcgpu.gate.h())
+
+# Oracle is equivalent to the identity gate, 
+# thus has no effect on the state
+
+constant_state.apply_all(qcgpu.gate.h())
+
+outcomes = constant_state.measure(samples = 1000)
+
+if int(max(outcomes, key=outcomes.get)) == 0:
+    print('constant')
+else:
+    print('balanced')

notebooks/Benchmarking.ipynb

@@ -1,188 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "code",
-   "execution_count": 1,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "import sys\n",
-    "import time\n",
-    "import pandas as pd\n",
-    "import numpy as np\n",
-    "\n",
-    "# QCGPU Imports\n",
-    "import qcgpu\n",
-    "\n",
-    "# ProjectQ Imports\n",
-    "from projectq import MainEngine\n",
-    "import projectq.ops as ops\n",
-    "from projectq.backends import Simulator\n",
-    "from projectq.types import Qureg\n",
-    "\n",
-    "# Qiskit Imports\n",
-    "from qiskit import ClassicalRegister, QuantumRegister\n",
-    "from qiskit import QuantumCircuit, execute"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 2,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "def bench_qcgpu(n, depth):\n",
-    "    state = qcgpu.State(n)\n",
-    "\n",
-    "    h = qcgpu.gate.h()\n",
-    "    x = qcgpu.gate.x()\n",
-    "    t = qcgpu.gate.t()\n",
-    "\n",
-    "    start = time.time()\n",
-    "\n",
-    "    for level in range(depth):\n",
-    "        for q in range(n):\n",
-    "    \n",
-    "            state.apply_gate(h, q)\n",
-    "            state.apply_gate(t, q)\n",
-    "\n",
-    "            if q != 0:\n",
-    "                state.apply_controlled_gate(x, q, 0)\n",
-    "        \n",
-    "    runtime = time.time() - start\n",
-    "    return runtime\n",
-    "    \n",
-    "def bench_projectq(n, depth):\n",
-    "    eng = MainEngine(backend=Simulator(gate_fusion=True), engine_list=[])\n",
-    "    qbits = eng.allocate_qureg(n)\n",
-    "\n",
-    "    start = time.time()\n",
-    "\n",
-    "    for level in range(depth):\n",
-    "        for q in qbits:\n",
-    "            ops.H | q\n",
-    "            ops.T | q\n",
-    "            if q != qbits[0]:\n",
-    "                ops.CNOT | (q, qbits[0])\n",
-    "\n",
-    "    runtime = time.time() - start\n",
-    "\n",
-    "    for q in qbits:\n",
-    "        ops.Measure | q\n",
-    "    eng.flush()\n",
-    "    return runtime\n",
-    "\n",
-    "def bench_qiskit(n, depth):\n",
-    "    q = QuantumRegister(n)\n",
-    "    c = ClassicalRegister(n)\n",
-    "    qc = QuantumCircuit(q, c)\n",
-    "\n",
-    "    \n",
-    "    \n",
-    "    for level in range(depth):\n",
-    "        for i in range(n):\n",
-    "            qc.h(q[i])\n",
-    "            qc.t(q[i])\n",
-    "\n",
-    "            if i != 0:\n",
-    "                qc.cx(q[i], q[0])\n",
-    "\n",
-    "    start = time.time()\n",
-    "    \n",
-    "    job_sim = execute(qc, \"local_statevector_simulator\")\n",
-    "    job_sim.result()\n",
-    "    \n",
-    "    runtime = time.time() - start\n",
-    "    return runtime\n"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 3,
-   "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "0.016415119171142578"
-      ]
-     },
-     "execution_count": 3,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "bench_qcgpu(20,3)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 4,
-   "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "0.20844078063964844"
-      ]
-     },
-     "execution_count": 4,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "bench_projectq(20,3)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 5,
-   "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "2.865260124206543"
-      ]
-     },
-     "execution_count": 5,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "bench_qiskit(20,3)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3",
-   "language": "python",
-   "name": "python3"
-  },
-  "language_info": {
-   "codemirror_mode": {
-    "name": "ipython",
-    "version": 3
-   },
-   "file_extension": ".py",
-   "mimetype": "text/x-python",
-   "name": "python",
-   "nbconvert_exporter": "python",
-   "pygments_lexer": "ipython3",
-   "version": "3.6.6"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}

notebooks/Verification.ipynb

@@ -1,134 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# Verification\n"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 49,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "import numpy as np\n",
-    "import qcgpu\n",
-    "from projectq import MainEngine\n",
-    "import projectq.ops as ops\n",
-    "from projectq.backends import Simulator, CircuitDrawer"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 50,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "def qcgpu_amplitudes(n, depth):\n",
-    "    state = qcgpu.State(n)\n",
-    "\n",
-    "    h = qcgpu.gate.h()\n",
-    "    x = qcgpu.gate.x()\n",
-    "    t = qcgpu.gate.t()\n",
-    "\n",
-    "    for level in range(depth):\n",
-    "        for q in range(n):\n",
-    "    \n",
-    "            state.apply_gate(h, q)\n",
-    "            state.apply_gate(t, q)\n",
-    "\n",
-    "            if q != 0:\n",
-    "                state.apply_controlled_gate(x, q, 0)\n",
-    "        \n",
-    "    return state.amplitudes()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 58,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "def projectq_amplitudes(n, depth):\n",
-    "#     drawing_engine = CircuitDrawer()\n",
-    "#     eng = MainEngine(drawing_engine)\n",
-    "    sim = Simulator(gate_fusion=True)\n",
-    "    eng = MainEngine(backend=sim, engine_list=[])\n",
-    "    qbits = eng.allocate_qureg(n)\n",
-    "\n",
-    "    for level in range(depth):\n",
-    "        for q in qbits:\n",
-    "            ops.H | q\n",
-    "            ops.T | q\n",
-    "            if q != qbits[0]:\n",
-    "                ops.CNOT | (q, qbits[0])\n",
-    "\n",
-    "\n",
-    "    eng.flush()\n",
-    "    amplitudes = sim.cheat()[1]\n",
-    "    for q in qbits:\n",
-    "        ops.Measure | q\n",
-    "   \n",
-    "    return amplitudes\n",
-    "#     return drawing_engine.get_latex()\n"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 70,
-   "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "True"
-      ]
-     },
-     "execution_count": 70,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "def verify(n, depth):\n",
-    "    amps_qcgpu = qcgpu_amplitudes(n,depth)\n",
-    "    amps_projectq = np.transpose(np.array([projectq_amplitudes(n,depth)]))\n",
-    "#     print(amps_qcgpu)\n",
-    "#     print(amps_projectq)\n",
-    "    return np.allclose(amps_qcgpu, amps_projectq, rtol=1e-04)\n",
-    "                        \n",
-    "verify(25,15)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 2",
-   "language": "python",
-   "name": "python2"
-  },
-  "language_info": {
-   "codemirror_mode": {
-    "name": "ipython",
-    "version": 2
-   },
-   "file_extension": ".py",
-   "mimetype": "text/x-python",
-   "name": "python",
-   "nbconvert_exporter": "python",
-   "pygments_lexer": "ipython2",
-   "version": "2.7.15rc1"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}

qcgpu/backend.py

@@ -86,52 +86,41 @@ __kernel void apply_controlled_gate(
     {
         amplitudes[one_state] = cfloat_add(cfloat_mul(D, one_amp), cfloat_mul(C, zero_amp));
     }
-
-
 }
 
 /*
  * Applies a controlled-controlled single qubit gate to the register.
- * NOT MIGRATED
  */
 __kernel void apply_controlled_controlled_gate(
-    __global cfloat_t *const amplitudes,
-    __global cfloat_t *amps,
-    int control1,
-    int control2,
+    __global cfloat_t *amplitudes,
+    int control,
+    int control_2,
     int target,
     cfloat_t A,
     cfloat_t B,
     cfloat_t C,
     cfloat_t D)
 {
-    int const state = get_global_id(0);
-    cfloat_t const amp = amplitudes[state];
+    int const global_id = get_global_id(0);
+    int const zero_state = nth_cleared(global_id, target);
+    int const one_state = zero_state | (1 << target); // Set the target bit
 
-    int const zero_state = state & (~(1 << target));
-    int const one_state = state | (1 << target);
+    int const control_val_zero = (((1 << control) & zero_state) > 0) ? 1 : 0;
+    int const control_val_one = (((1 << control) & one_state) > 0) ? 1 : 0;
+    int const control_val_two_zero = (((1 << control_2) & zero_state) > 0) ? 1 : 0;
+    int const control_val_two_one = (((1 << control_2) & one_state) > 0) ? 1 : 0;
 
-    int const bit_val = (((1 << target) & state) > 0) ? 1 : 0;
-    int const control1_val = (((1 << control1) & state) > 0) ? 1 : 0;
-    int const control2_val = (((1 << control2) & state) > 0) ? 1 : 0;
+    cfloat_t const zero_amp = amplitudes[zero_state];
+    cfloat_t const one_amp = amplitudes[one_state];
 
-    if (control1_val == 0 || control2_val == 0)
+    if (control_val_zero == 1 && control_val_two_zero == 1)
     {
-        // Control is 0, don't apply gate
-        amps[state] = amp;
+        amplitudes[zero_state] = cfloat_add(cfloat_mul(A, zero_amp), cfloat_mul(B, one_amp));
     }
-    else
+
+    if (control_val_one == 1 && control_val_two_one == 1)
     {
-        // control is 1, apply gate.
-        if (bit_val == 0)
-        {
-            // Bitval = 0
-            amps[state] = cfloat_add(cfloat_mul(A, amp), cfloat_mul(B, amplitudes[one_state]));
-        }
-        else
-        {
-            amps[state] = cfloat_add(cfloat_mul(D, amp), cfloat_mul(C, amplitudes[zero_state]));
-        }
+        amplitudes[one_state] = cfloat_add(cfloat_mul(D, one_amp), cfloat_mul(C, zero_amp));
     }
 }
 
@@ -289,6 +278,23 @@ class Backend:
             self.dtype(gate.d)
         )
     
+    def apply_controlled_controlled_gate(self, gate, control1, control2, target):
+        """Applies a controlled controlled gate (such as a toffoli gate) to the quantum register"""
+
+        program.apply_controlled_controlled_gate(
+            self.queue,
+            [int(2**self.num_qubits / 2)],
+            None,
+            self.buffer.data,
+            np.int32(control1),
+            np.int32(control2),
+            np.int32(target),
+            self.dtype(gate.a),
+            self.dtype(gate.b),
+            self.dtype(gate.c),
+            self.dtype(gate.d)
+        )
+
     def seed(self, val):
         random.seed(val)
         
@@ -311,6 +317,23 @@ class Backend:
             results[np.binary_repr(i, width=self.num_qubits)] += 1
         
         return dict(results)
+
+    def measure_first(self, num, samples):
+        probabilities = self.probabilities()
+        # print(probabilities)
+        # print(np.sum(self.amplitudes()))
+        choices = np.random.choice(
+            np.arange(0, 2**self.num_qubits), 
+            samples, 
+            p=probabilities
+        )
+        
+        results = defaultdict(int)
+        for i in choices:
+            key = np.binary_repr(i, width=self.num_qubits)[-num:]
+            results[key] += 1
+        
+        return dict(results)
        
 
     def qubit_probability(self, target):

qcgpu/gate.py

@@ -14,7 +14,7 @@ def memoize(func):
     return memoized_func
 
 class Gate:
-    def __init__(self, gate):
+    def __init__(self, gate, unitary=True):
         gate = np.array(gate)
         
         if gate.shape != (2, 2):
@@ -25,8 +25,9 @@ class Gate:
             )
 
         # Check the gate is unitary
-        if (not np.allclose(np.eye(gate.shape[0]), np.dot(gate.conjugate().transpose(), gate))):
-            raise ValueError("gate is not unitary.")
+        if unitary:
+            if (not np.allclose(np.eye(gate.shape[0]), np.dot(gate.conjugate().transpose(), gate))):
+                raise ValueError("gate is not unitary.")
 
         self.a = complex(gate[0, 0])
         self.b = complex(gate[0, 1])

qcgpu/state.py

@@ -89,6 +89,20 @@ class State:
         # TODO: Check that gate is correct
 
         self.backend.apply_controlled_gate(gate, control, target)
+    
+    def apply_controlled_controlled_gate(self, gate, control, control2, target):
+        if not isinstance(target, int) or target < 0:
+            raise ValueError("target must be an int >= 0")
+        
+        if not isinstance(control, int) or control < 0:
+            raise ValueError("control must be an int >= 0")
+
+        if not isinstance(control2, int) or control2 < 0:
+            raise ValueError("control must be an int >= 0")
+
+        # TODO: Check that gate is correct
+
+        self.backend.apply_controlled_controlled_gate(gate, control, control2, target)
 
     def measure_qubit(self, target, samples=1):
         return self.backend.measure_qubit(target, samples)
@@ -99,6 +113,9 @@ class State:
     def measure(self, samples=1):
         return self.backend.measure(samples)
 
+    def measure_first(self, num=1, samples=1):
+        return self.backend.measure_first(num, samples)
+
     def amplitudes(self):
         return self.backend.amplitudes()[0]
     
@@ -140,6 +157,9 @@ class State:
     def cnot(self, control, target):
         self.apply_controlled_gate(qcgpu.gate.x(), control, target)
 
+    def toffoli(self, control, control2, target):
+        self.apply_controlled_controlled_gate(qcgpu.gate.x(), control, control2, target)
+
     def u(self, target, theta, phi, lda):
         a = np.exp(-1j * (phi + lda) / 2) * np.cos(theta / 2)
         b = - np.exp(-1j * (phi - lda) / 2) * np.sin(theta / 2)