Parallel Qubit Clustering/Parallel Qubit Clustering Batch.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Parallel Qubit Clustering Experiment\n",
    "\n",
    "In this notebook, all the code that is used for the Parallel Qubit Clustering Experiment can be found. The code is organised in such a way that the experiment can be run multiple times with various different configurations. Towards the bottom of the file, you are presented with an array of experiment objects that contain different configurations. These are the configurations that is also referred to in the report. Running this code will perform each experiment, the specified number of times, and calculate the average runtime of each experiment. The experiments are performed sequentially.\n",
    "\n",
    "The comments in this code should provide a good explanation of the workflow. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# =========================================================================\n",
    "# Import all necessary libraries\n",
    "# =========================================================================\n",
    "\n",
    "# Special imports (necessary for backrefrencing to own class-type)\n",
    "from __future__ import annotations\n",
    "\n",
    "# Standard Python libraries\n",
    "import time\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "\n",
    "# Numpy and SciPy\n",
    "import numpy as np\n",
    "from scipy.optimize import linear_sum_assignment as linear_assignment\n",
    "\n",
    "# Qiskit\n",
    "from qiskit import Aer\n",
    "from qiskit import QuantumCircuit\n",
    "from qiskit import QuantumRegister\n",
    "from qiskit import ClassicalRegister\n",
    "from qiskit import transpile\n",
    "from qiskit import execute\n",
    "from qiskit.algorithms.optimizers import SPSA\n",
    "from qiskit.algorithms.optimizers import ADAM\n",
    "from qiskit.algorithms.optimizers import GradientDescent\n",
    "from qiskit.algorithms.optimizers import COBYLA\n",
    "from qiskit.quantum_info import state_fidelity\n",
    "\n",
    "# Sklearn\n",
    "from sklearn import datasets\n",
    "\n",
    "# Other libraries\n",
    "from tqdm import tqdm\n",
    "\n",
    "\n",
    "\n",
    "# =========================================================================\n",
    "# Utility functions\n",
    "# =========================================================================\n",
    "\n",
    "# Square a number\n",
    "def sqr(x: float):\n",
    "    return x * x\n",
    "\n",
    "# Map a number 'a', which is the range between 'b' and 'c',\n",
    "# to become in the range between 'd' and 'e'.\n",
    "def map(a: float, b: float, c: float, d: float, e: float) -> float:\n",
    "    return d + ((a - b)/(c - b)) * (e - d)\n",
    "\n",
    "\n",
    "\n",
    "# =========================================================================\n",
    "# DataPoints\n",
    "# =========================================================================\n",
    "\n",
    "# DataPoint class\n",
    "# -> Provide framework around data-points in the dataset that helps\n",
    "#    track changes throughout the algorithm.\n",
    "class DataPoint:\n",
    "    x: float\n",
    "    y: float\n",
    "    cluster: int\n",
    "    fidelity: float\n",
    "\n",
    "    def __init__(self, x: float, y: float):\n",
    "        self.x = x\n",
    "        self.y = y\n",
    "        self.cluster = 0\n",
    "        self.fidelity = 0\n",
    "\n",
    "    def dist(self, o: DataPoint) -> float:\n",
    "        return sqr(o.x - self.x) + sqr(o.y - self.y)\n",
    "    \n",
    "    def reset_fidelity(self):\n",
    "        self.cluster = 0\n",
    "        self.fidelity = 0\n",
    "        \n",
    "        \n",
    "        \n",
    "# =========================================================================\n",
    "# Experiments\n",
    "# =========================================================================\n",
    "\n",
    "# ParallelQubitExperiment class\n",
    "# -> Provide an easy way to execute multiple experiments in a row\n",
    "class ParallelQubitExperiment:\n",
    "    name: str\n",
    "        \n",
    "    # Experiment configuration\n",
    "    dataset: [DataPoint]\n",
    "    num_qubits: int\n",
    "    num_shots: int\n",
    "    num_iter: int\n",
    "    optimizer = None\n",
    "        \n",
    "    # Information about the dataset\n",
    "    xrange: [float]\n",
    "    yrange: [float]\n",
    "        \n",
    "    # The number of optimization parameters\n",
    "    num_params: int\n",
    "    \n",
    "    # Keep track of the number of circuits that is run under this experiment\n",
    "    parallel_circuit_counter: int\n",
    "    \n",
    "    # The selected qiskit backend\n",
    "    qiskit_backend = None\n",
    "    \n",
    "    # Initialization function: provide all necessary information\n",
    "    # before starting the experiment.\n",
    "    def __init__(self, name: str, data: [DataPoint], num_qubits: int, num_shots: int, num_iter: int, optimizer, num_params: int):\n",
    "        self.name = name\n",
    "        \n",
    "        self.dataset = data\n",
    "        \n",
    "        # Determine the xmin/max and ymin/max\n",
    "        self.xrange = [self.dataset[0].x, self.dataset[0].x]\n",
    "        self.yrange = [self.dataset[0].y, self.dataset[0].y]\n",
    "        \n",
    "        for dp in self.dataset:\n",
    "            if dp.x < self.xrange[0]:\n",
    "                self.xrange[0] = dp.x\n",
    "            if dp.x > self.xrange[1]:\n",
    "                self.xrange[1] = dp.x\n",
    "            if dp.y < self.yrange[0]:\n",
    "                self.yrange[0] = dp.y\n",
    "            if dp.y > self.yrange[1]:\n",
    "                self.yrange[1] = dp.y\n",
    "        \n",
    "        self.num_qubits = num_qubits\n",
    "        self.num_shots = num_shots\n",
    "        self.num_iter = num_iter\n",
    "        self.parallel_circuit_counter = 0\n",
    "        \n",
    "        # Install the optimizer\n",
    "        self.optimizer = COBYLA(maxiter=self.num_iter, tol=0.00001)\n",
    "        \n",
    "        self.num_params = num_params\n",
    "        \n",
    "        # Select the Aer_Simulator backend\n",
    "        self.qiskit_backend = Aer.get_backend(\"aer_simulator\")\n",
    "    \n",
    "    # \n",
    "    # Methods\n",
    "    #\n",
    "    \n",
    "    # [Private] Construct and run the variational circuit for the provided set of DataPoints.\n",
    "    # Extract the 'counts': the number of times a qubit ended up in zero or one state.\n",
    "    # Compute the fidelities from these counts.\n",
    "    # Return the computed fidelities\n",
    "    def __run_parallel_variational_circuit(self, dps: [DataPoint], params):\n",
    "        # Declare a QuantumCircuit with the required number of qubits\n",
    "        registers = [QuantumRegister(1) for i in range(min(len(dps), self.num_qubits))]\n",
    "        \n",
    "        q_circuit = QuantumCircuit(*registers)\n",
    "        \n",
    "        # Name the circuit\n",
    "        q_circuit.name = \"parallel_circuit_\" + str(self.parallel_circuit_counter)\n",
    "        self.parallel_circuit_counter += 1\n",
    "        \n",
    "        # Install the circuit\n",
    "        # Initialize all qubits to state zero\n",
    "        for i in registers:\n",
    "            q_circuit.initialize([1, 0], i)\n",
    "        \n",
    "        # Install the unitary transformations\n",
    "        for i in range(len(registers)):\n",
    "            # 1: unitary transformation for data point properties (x, y)\n",
    "            q_circuit.u(\n",
    "                map(dps[i].x, self.xrange[0], self.xrange[1], 0, math.pi),\n",
    "                map(dps[i].y, self.yrange[0], self.yrange[1], 0, math.pi),\n",
    "                0,\n",
    "                registers[i]\n",
    "            )\n",
    "            \n",
    "            # 2: unitary transformation for the optimization parameters\n",
    "            if len(params) == 1:\n",
    "                q_circuit.u(params[0], 0, 0, registers[i])\n",
    "            elif len(params) == 2:\n",
    "                q_circuit.u(params[0], params[1], 0, registers[i])\n",
    "            elif len(params) >= 3:\n",
    "                q_circuit.u(params[0], params[1], params[2], registers[i])\n",
    "                \n",
    "        # Measure all qubits\n",
    "        q_circuit.measure_all()\n",
    "        \n",
    "        # Run the circuit and obtain the result\n",
    "        job_run = execute(q_circuit, backend=self.qiskit_backend, shots=self.num_shots)\n",
    "        job_counts = job_run.result().get_counts(q_circuit)\n",
    "        \n",
    "        # Count the number of times the individual qubits read zero or one as\n",
    "        # measurement outcome\n",
    "        counts_zero = []\n",
    "        counts_one  = []\n",
    "\n",
    "        for i in range(len(dps)):\n",
    "            counts_zero.append(0)\n",
    "            counts_one.append(0)\n",
    "\n",
    "        # Go over all the outcomes that qiskit encountered\n",
    "        for key in job_counts:\n",
    "            cnt = job_counts[key]\n",
    "            \n",
    "            # For each data point (that is connected to a qubit):\n",
    "            # -> If, during this encounter, this qubit is '0', add the count to\n",
    "            #    its zero-counter in the array. Likewise with '1'.\n",
    "            for i in range(len(dps)):\n",
    "                char = key[i]\n",
    "                if char == '0':\n",
    "                    counts_zero[i] += cnt\n",
    "                elif char == '1':\n",
    "                    counts_one[i] += cnt\n",
    "                else:\n",
    "                    print(\"UNEXPECTED MEASUREMENT OUTCOME!\")\n",
    "\n",
    "        # Calculate fidelities with reference states zero and one\n",
    "        fidelities_with_zero = []\n",
    "\n",
    "        # This is done by dividing the counts by the number of shots\n",
    "        for i in range(len(dps)):\n",
    "            fidelities_with_zero.append(float(counts_zero[i])/float(self.num_shots))\n",
    "\n",
    "        # Before returning the fidelities for each data-point with zero state,\n",
    "        # reverse the order of this array. Qubit follows a different order.\n",
    "        fidelities_with_zero.reverse()\n",
    "        return fidelities_with_zero\n",
    "    \n",
    "    \n",
    "    # [Private] Calculate the cost-function for the provided set of parameters\n",
    "    def __parallel_objective(self, params) -> float:\n",
    "        # Start by pulling every data point through the circuit.\n",
    "        # Create queues of data points, with max length the number of qubits\n",
    "        circuit_queues = []\n",
    "        \n",
    "        queue = []\n",
    "        for dp in self.dataset:\n",
    "            queue.append(dp)\n",
    "            if len(queue) >= self.num_qubits:\n",
    "                circuit_queues.append(queue)\n",
    "                queue = []\n",
    "        if not len(queue) == 0:\n",
    "            circuit_queues.append(queue)\n",
    "            \n",
    "        # For every data point, compute the fidelity with the zero-state\n",
    "        # by using the queues we just made\n",
    "        \n",
    "        # Start by resetting the fidelities and cluster-assignments of all datapoints\n",
    "        for dp in self.dataset:\n",
    "            dp.reset_fidelity()\n",
    "        \n",
    "        datapoint_fid_with_zero = []\n",
    "        \n",
    "        for queue in circuit_queues:\n",
    "            # Compute fidelities for this queue by using class function\n",
    "            queue_fid_with_zero = self.__run_parallel_variational_circuit(queue, params)\n",
    "            \n",
    "            # In the order they're presented in, append these fidelities to the array\n",
    "            # of data point-fidelities\n",
    "            for fid in queue_fid_with_zero:\n",
    "                datapoint_fid_with_zero.append(fid)\n",
    "                \n",
    "        # Now that we have the fidelities for each data point with zero state, compute\n",
    "        # the value of the objective function\n",
    "        total_cost = 0\n",
    "    \n",
    "        for i in range(len(self.dataset)):\n",
    "            for j in range(len(self.dataset)):\n",
    "                if not i == j:\n",
    "                    point_dist = self.dataset[i].dist(self.dataset[j])\n",
    "\n",
    "                    # Compute the objective function for reference state zero\n",
    "                    cost_zero = point_dist * datapoint_fid_with_zero[i] * datapoint_fid_with_zero[j]\n",
    "\n",
    "                    # Compute the objective function for reference state one\n",
    "                    cost_one = point_dist * (1-datapoint_fid_with_zero[i]) * (1-datapoint_fid_with_zero[j])\n",
    "\n",
    "                    # Add both costs to the total cost function\n",
    "                    total_cost += cost_zero\n",
    "                    total_cost += cost_one\n",
    "                    \n",
    "        # Iterate over all datapoints and update the cluster they belong to\n",
    "        for i, dp in enumerate(self.dataset):\n",
    "            if datapoint_fid_with_zero[i] > dp.fidelity:\n",
    "                dp.fidelity = datapoint_fid_with_zero[i]\n",
    "                dp.cluster = 0\n",
    "            if 1-datapoint_fid_with_zero[i] > dp.fidelity:\n",
    "                dp.fidelity = 1-datapoint_fid_with_zero[i]\n",
    "                dp.cluster = 1\n",
    "        \n",
    "        # We're done coputing the final cost! Return it!\n",
    "        return total_cost\n",
    "    \n",
    "    # [Public] Render experiment name\n",
    "    def render_name(self):\n",
    "        return self.name + \" [\" + str(self.num_qubits) + \" qubits, \" + str(self.num_shots) + \" shots, \" + str(self.num_iter) + \" iterations]\"\n",
    "    \n",
    "    # [Public] Run the experiment a specified amount of times.\n",
    "    #\n",
    "    # - Render a scatter plot of the resulting cluster and save it on disk.\n",
    "    # - Record the time each experiment takes. Average the runtimes and return it.\n",
    "    def run(self, times: int) -> float:\n",
    "        average_runtime = 0\n",
    "        \n",
    "        print(\"\\n\\n\")\n",
    "        print(\"Starting experiment \" + self.render_name())\n",
    "        \n",
    "        for i in tqdm(range(times)):            \n",
    "            # Record start time\n",
    "            start_time = time.time()\n",
    "        \n",
    "            # Prepare random initial optimization parameters\n",
    "            params_init = np.random.rand(self.num_params)\n",
    "\n",
    "            # Perform the optimization and obtain the optimal parameters\n",
    "            params_optimal = self.optimizer.minimize(fun=self.__parallel_objective, x0=params_init).x\n",
    "            \n",
    "            # Calculate the run-time and add it to the average\n",
    "            run_time = time.time() - start_time\n",
    "            average_runtime += run_time\n",
    "            \n",
    "            # Separate the points based on the cluster they're in\n",
    "            clusters = [[], []]\n",
    "            \n",
    "            for dp in self.dataset:\n",
    "                clusters[dp.cluster].append(dp)\n",
    "            \n",
    "            # Make a scatter plot of this result\n",
    "            for c in clusters:\n",
    "                plt.scatter(\n",
    "                    [dp.x for dp in c],\n",
    "                    [dp.y for dp in c]\n",
    "                )\n",
    "            \n",
    "            # Give the plot the name of this experiment, together with the run index\n",
    "            plt.title(self.render_name() + \" run \" + str(i))\n",
    "            \n",
    "            # Save the figure to disk\n",
    "            plt.savefig(\"./Images/\" + self.render_name() + \".png\")\n",
    "            \n",
    "            \n",
    "        # Finally, divide the sum of runtimes by the number of runs to obtain the average\n",
    "        # and return it.\n",
    "        average_runtime /= times\n",
    "        return average_runtime\n",
    "\n",
    "    \n",
    "# =========================================================================\n",
    "# Setup and run the experiments\n",
    "# =========================================================================\n",
    "\n",
    "# Select a dataset\n",
    "selected_dataset = datasets.make_blobs(n_samples=100, centers=2, random_state=8)\n",
    "\n",
    "data_points = []\n",
    "for point in selected_dataset[0]:\n",
    "    x = point[0]\n",
    "    y = point[1]\n",
    "    data_points.append(DataPoint(x, y))\n",
    "\n",
    "# Setup the experiments\n",
    "experiments = [\n",
    "    ParallelQubitExperiment(\n",
    "        name=\"Parallel Qubit\",\n",
    "        data=data_points,\n",
    "        num_qubits=1,\n",
    "        num_shots=256,\n",
    "        num_iter=100,\n",
    "        optimizer=None,\n",
    "        num_params=3\n",
    "    ),\n",
    "    ParallelQubitExperiment(\n",
    "        name=\"Parallel Qubit\",\n",
    "        data=data_points,\n",
    "        num_qubits=2,\n",
    "        num_shots=256,\n",
    "        num_iter=100,\n",
    "        optimizer=None,\n",
    "        num_params=3\n",
    "    ),\n",
    "    ParallelQubitExperiment(\n",
    "        name=\"Parallel Qubit\",\n",
    "        data=data_points,\n",
    "        num_qubits=3,\n",
    "        num_shots=256,\n",
    "        num_iter=100,\n",
    "        optimizer=None,\n",
    "        num_params=3\n",
    "    ),\n",
    "    ParallelQubitExperiment(\n",
    "        name=\"Parallel Qubit\",\n",
    "        data=data_points,\n",
    "        num_qubits=4,\n",
    "        num_shots=256,\n",
    "        num_iter=100,\n",
    "        optimizer=None,\n",
    "        num_params=3\n",
    "    ),\n",
    "    ParallelQubitExperiment(\n",
    "        name=\"Parallel Qubit\",\n",
    "        data=data_points,\n",
    "        num_qubits=5,\n",
    "        num_shots=256,\n",
    "        num_iter=100,\n",
    "        optimizer=None,\n",
    "        num_params=3\n",
    "    ),\n",
    "    ParallelQubitExperiment(\n",
    "        name=\"Parallel Qubit\",\n",
    "        data=data_points,\n",
    "        num_qubits=6,\n",
    "        num_shots=256,\n",
    "        num_iter=100,\n",
    "        optimizer=None,\n",
    "        num_params=3\n",
    "    ),\n",
    "    ParallelQubitExperiment(\n",
    "        name=\"Parallel Qubit\",\n",
    "        data=data_points,\n",
    "        num_qubits=7,\n",
    "        num_shots=256,\n",
    "        num_iter=100,\n",
    "        optimizer=None,\n",
    "        num_params=3\n",
    "    ),\n",
    "    ParallelQubitExperiment(\n",
    "        name=\"Parallel Qubit\",\n",
    "        data=data_points,\n",
    "        num_qubits=8,\n",
    "        num_shots=256,\n",
    "        num_iter=100,\n",
    "        optimizer=None,\n",
    "        num_params=3\n",
    "    )\n",
    "]\n",
    "    \n",
    "# Run all experiments 20 times\n",
    "for e in experiments:\n",
    "    runtime = e.run(20)\n",
    "    print(\"Runtime for experiment \" + e.render_name() + \": \" + str(runtime))\n",
    "\n",
    "\n"
   ]
  },
  {
   "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}