challenge-3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## IBM Quantum Challenge Fall 2021\n",
    "\n",
    "# Challenge 3: Classify images with quantum machine learning\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "    \n",
    "We recommend that you switch to **light** workspace theme under the Account menu in the upper right corner for optimal experience."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "Machine learning is a technology that has attracted a great deal of attention due to its high performance and versatility. In fact, it has been put to practical use in many industries with the recent development of algorithms and the increase of computational resources. A typical example is computer vision, where machine learning is now able to classify images with the same or better accuracy than humans. For example, the ability to automatically classify clothing images has made online shopping for clothes more convenient.\n",
    "\n",
    "The application of quantum computation to machine learning has recently been shown to have the potential for even greater capabilities. Various algorithms have been proposed for quantum machine learning, such as the quantum support vector machine (QSVM) and quantum generative adversarial networks (QGANs). In this challenge, you will use QSVM to tackle the clothing image classification task.\n",
    "\n",
    "QSVM is a quantum version of the support vector machine (SVM), a classical machine learning algorithm. There are various approaches to QSVM, some aim to accelerate computation assuming fault-tolerant quantum computers, while others aim to achieve higher expressive power assuming noisy, near-term devices. In this challenge, we will focus on the latter, and the details will be explained later.\n",
    "\n",
    "For this implementation of QSVM, you will be able to make choices on how you want to compose your quantum model, in particular focusing on the quantum feature map. This is motivated by the tradeoff that a more complex feature map would have greater representation power but be more susceptible to noise, which could be especially critical when using noisy, near-term devices.\n",
    "\n",
    "Many of the concepts that appear in this challenge are explained in the 2021 Qiskit Global Summer School (QGSS). The materials and lecture videos are available, and it is recommended that you study them as well. Refer to the links in each part for the corresponding lectures.\n",
    "\n",
    "<center><img src=\"./resources/ecommerce.jpg\" width=\"640\" /></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Challenge\n",
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Goal**\n",
    "\n",
    "Implement a QSVM model for multiclass classification and predict labels accurately. \n",
    "    \n",
    "**Plan**\n",
    "\n",
    "First, you will learn how to construct QSVM for binary classification of a simple dataset. Then apply what you have learned to a more complex problem, 3-class classification of a different dataset.\n",
    "\n",
    "**1. Tutorial - QSVM for binary classification of MNIST:** familiarize yourself with a typical workflow for QSVM and find the best combination of dimentions/feature maps.\n",
    "\n",
    "**2. Challenge - QSVM for 3-class classification of Fashion-MNIST:** implement a 3-class classifier using binary QSVM classifers. Perform similar investigation as in the first part to find the best combination of dimentions/feature maps. Achieve better accuracy with smaller feature map circuits.\n",
    "\n",
    "</div>\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "\n",
    "Before you begin, we recommend watching the [**Qiskit Machine Learning Demo Session with Anton Dekusar**](https://youtu.be/claoY57eVIc?t=1814) and check out the corresponding [**demo notebook**](https://github.com/qiskit-community/qiskit-application-modules-demo-sessions/tree/main/qiskit-machine-learning) to learn how to do classifications using QSVM.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# General imports\n",
    "import os\n",
    "import gzip\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from pylab import cm\n",
    "import warnings\n",
    "\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "# scikit-learn imports\n",
    "from sklearn import datasets\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.preprocessing import StandardScaler, MinMaxScaler\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.svm import SVC\n",
    "from sklearn.metrics import accuracy_score\n",
    "\n",
    "# Qiskit imports\n",
    "from qiskit import Aer, execute\n",
    "from qiskit.circuit import QuantumCircuit, Parameter, ParameterVector\n",
    "from qiskit.circuit.library import PauliFeatureMap, ZFeatureMap, ZZFeatureMap\n",
    "from qiskit.circuit.library import TwoLocal, NLocal, RealAmplitudes, EfficientSU2\n",
    "from qiskit.circuit.library import HGate, RXGate, RYGate, RZGate, CXGate, CRXGate, CRZGate\n",
    "from qiskit_machine_learning.kernels import QuantumKernel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 1: Tutorial - QSVM for binary classification of MNIST\n",
    "\n",
    "In this part, you will apply QSVM to the binary classification of handwritten numbers 4 and 9. Through this tutorial, you will learn the workflow of applying QSVM to binary classification. Find better combinations and achieve higher accuracy.\n",
    "\n",
    "Related QGSS material:\n",
    "- [**Lab 3**](https://www.youtube.com/watch?v=GVhCOTzAkCM&list=PLOFEBzvs-VvqJwybFxkTiDzhf5E11p8BI&index=17)\n",
    "\n",
    "### 1. Data preparation\n",
    "\n",
    "The data we are going to work with at the beginning is a small subset of the well known handwritten digits dataset, which is available publicly. We will be aiming to differentiate between '4' and '9'. \n",
    "\n",
    "There are a total of 100 data in the dataset. Of these, eighty are labeled training data, and the remaining twenty are unlabeled test data. Each data is a 28x28 image of a digit, collapsed into an array, where each element is an integer between 0 (white) and 255 (black). To use the dataset for quantum classification, we need to scale the range to between -1 and 1, and reduce the dimensionality to the number of qubits we want to use (here N_DIM=5)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load MNIST dataset\n",
    "DATA_PATH = './resources/ch3_part1.npz'\n",
    "data = np.load(DATA_PATH)\n",
    "\n",
    "sample_train = data['sample_train']\n",
    "labels_train = data['labels_train']\n",
    "sample_test = data['sample_test']\n",
    "\n",
    "# Split train data\n",
    "sample_train, sample_val, labels_train, labels_val = train_test_split(\n",
    "    sample_train, labels_train, test_size=0.2, random_state=42)\n",
    "\n",
    "# Visualize samples\n",
    "fig = plt.figure()\n",
    "\n",
    "LABELS = [4, 9]\n",
    "num_labels = len(LABELS)\n",
    "for i in range(num_labels):\n",
    "    ax = fig.add_subplot(1, num_labels, i+1)\n",
    "    img = sample_train[labels_train==LABELS[i]][0].reshape((28, 28))\n",
    "    ax.imshow(img, cmap=\"Greys\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Standardize\n",
    "ss = StandardScaler()\n",
    "sample_train = ss.fit_transform(sample_train)\n",
    "sample_val = ss.transform(sample_val)\n",
    "sample_test = ss.transform(sample_test)\n",
    "\n",
    "# Reduce dimensions\n",
    "N_DIM = 5\n",
    "pca = PCA(n_components=N_DIM)\n",
    "sample_train = pca.fit_transform(sample_train)\n",
    "sample_val = pca.transform(sample_val)\n",
    "sample_test = pca.transform(sample_test)\n",
    "\n",
    "# Normalize\n",
    "mms = MinMaxScaler((-1, 1))\n",
    "sample_train = mms.fit_transform(sample_train)\n",
    "sample_val = mms.transform(sample_val)\n",
    "sample_test = mms.transform(sample_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Data Encoding\n",
    "\n",
    "We will take the classical data and encode it to the quantum state space using a quantum feature map. The choice of which feature map to use is important and may depend on the given dataset we want to classify. Here we'll look at the feature maps available in Qiskit, before selecting and customising one to encode our data.\n",
    "\n",
    "### 2.1 Quantum Feature Maps\n",
    "As the name suggests, a quantum feature map $\\phi(\\mathbf{x})$ is a map from the classical feature vector $\\mathbf{x}$ to the quantum state $|\\Phi(\\mathbf{x})\\rangle\\langle\\Phi(\\mathbf{x})|$. This is facilitated by applying the unitary operation $\\mathcal{U}_{\\Phi(\\mathbf{x})}$ on the initial state $|0\\rangle^{n}$ where _n_ is the number of qubits being used for encoding.\n",
    "\n",
    "The following feature maps currently available in Qiskit are those introduced in [**_Havlicek et al_.  Nature **567**, 209-212 (2019)**](https://www.nature.com/articles/s41586-019-0980-2), in particular the `ZZFeatureMap` is conjectured to be hard to simulate classically and can be implemented as short-depth circuits on near-term quantum devices.\n",
    "\n",
    "- [**`PauliFeatureMap`**](https://qiskit.org/documentation/stubs/qiskit.circuit.library.PauliFeatureMap.html)\n",
    "- [**`ZZFeatureMap`**](https://qiskit.org/documentation/stubs/qiskit.circuit.library.ZFeatureMap.html)\n",
    "- [**`ZFeatureMap`**](https://qiskit.org/documentation/stubs/qiskit.circuit.library.ZZFeatureMap.html)\n",
    "\n",
    "The `PauliFeatureMap` is defined as:\n",
    "\n",
    "```python\n",
    "PauliFeatureMap(feature_dimension=None, reps=2, \n",
    "                entanglement='full', paulis=None, \n",
    "                data_map_func=None, parameter_prefix='x',\n",
    "                insert_barriers=False)\n",
    "```\n",
    "\n",
    "and describes the unitary operator of depth $d$:\n",
    "\n",
    "$$ \\mathcal{U}_{\\Phi(\\mathbf{x})}=\\prod_d U_{\\Phi(\\mathbf{x})}H^{\\otimes n},\\ U_{\\Phi(\\mathbf{x})}=\\exp\\left(i\\sum_{S\\subseteq[n]}\\phi_S(\\mathbf{x})\\prod_{k\\in S} P_i\\right), $$\n",
    "\n",
    "which contains layers of Hadamard gates interleaved with entangling blocks, $U_{\\Phi(\\mathbf{x})}$, encoding the classical data as shown in circuit diagram below for $d=2$.\n",
    "\n",
    "<center><img src=\"./resources/featuremap.png\" width=\"1000\" /></center>\n",
    "\n",
    "Within the entangling blocks, $U_{\\Phi(\\mathbf{x})}$: $P_i \\in \\{ I, X, Y, Z \\}$ denotes the Pauli matrices, the index $S$ describes connectivities between different qubits or datapoints: $S \\in \\{\\binom{n}{k}\\ combinations,\\ k = 1,... n \\}$, and by default the data mapping function $\\phi_S(\\mathbf{x})$ is \n",
    "$$\\phi_S:\\mathbf{x}\\mapsto \\Bigg\\{\\begin{array}{ll}\n",
    "    x_i & \\mbox{if}\\ S=\\{i\\} \\\\\n",
    "        (\\pi-x_i)(\\pi-x_j) & \\mbox{if}\\ S=\\{i,j\\}\n",
    "    \\end{array}$$\n",
    "\n",
    "when $k = 1, P_0 = Z$, this is the `ZFeatureMap`: \n",
    "$$\\mathcal{U}_{\\Phi(\\mathbf{x})} = \\left( \\exp\\left(i\\sum_j \\phi_{\\{j\\}}(\\mathbf{x}) \\, Z_j\\right) \\, H^{\\otimes n} \\right)^d.$$\n",
    "\n",
    "which is defined as:\n",
    "```python\n",
    "ZFeatureMap(feature_dimension, reps=2, \n",
    "            data_map_func=None, insert_barriers=False)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3 features, depth 2\n",
    "map_z = ZFeatureMap(feature_dimension=3, reps=2)\n",
    "map_z.decompose().draw('mpl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note the lack of entanglement in this feature map, this means that this feature map is simple to simulate classically and will not provide quantum advantage.\n",
    "\n",
    "and when $k = 2, P_0 = Z, P_1 = ZZ$, this is the `ZZFeatureMap`: \n",
    "$$\\mathcal{U}_{\\Phi(\\mathbf{x})} = \\left( \\exp\\left(i\\sum_{jk} \\phi_{\\{j,k\\}}(\\mathbf{x}) \\, Z_j \\otimes Z_k\\right) \\, \\exp\\left(i\\sum_j \\phi_{\\{j\\}}(\\mathbf{x}) \\, Z_j\\right) \\, H^{\\otimes n} \\right)^d.$$ \n",
    "\n",
    "which is defined as:\n",
    "```python\n",
    "ZZFeatureMap(feature_dimension, reps=2, \n",
    "             entanglement='full', data_map_func=None, \n",
    "             insert_barriers=False)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3 features, depth 1, linear entanglement\n",
    "map_zz = ZZFeatureMap(feature_dimension=3, reps=1, entanglement='linear')\n",
    "map_zz.decompose().draw('mpl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that there is entanglement in the feature map, we can define the entanglement map:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3 features, depth 1, circular entanglement\n",
    "map_zz = ZZFeatureMap(feature_dimension=3, reps=1, entanglement='circular')\n",
    "map_zz.decompose().draw('mpl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can customise the Pauli gates in the feature map, for example, $P_0 = X, P_1 = Y, P_2 = ZZ$:\n",
    "$$\\mathcal{U}_{\\Phi(\\mathbf{x})} = \\left( \\exp\\left(i\\sum_{jk} \\phi_{\\{j,k\\}}(\\mathbf{x}) \\, Z_j \\otimes Z_k\\right) \\, \\exp\\left(i\\sum_{j} \\phi_{\\{j\\}}(\\mathbf{x}) \\, Y_j\\right) \\, \\exp\\left(i\\sum_j \\phi_{\\{j\\}}(\\mathbf{x}) \\, X_j\\right) \\, H^{\\otimes n} \\right)^d.$$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3 features, depth 1\n",
    "map_pauli = PauliFeatureMap(feature_dimension=3, reps=1, paulis = ['X', 'Y', 'ZZ'])\n",
    "map_pauli.decompose().draw('mpl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The [`NLocal`](https://qiskit.org/documentation/stubs/qiskit.circuit.library.NLocal.html) and [`TwoLocal`](https://qiskit.org/documentation/stubs/qiskit.circuit.library.TwoLocal.html) functions in Qiskit's circuit library can also be used to create parameterised quantum circuits as feature maps. \n",
    "\n",
    "```python\n",
    "TwoLocal(num_qubits=None, reps=3, rotation_blocks=None, \n",
    "         entanglement_blocks=None, entanglement='full',  \n",
    "         skip_unentangled_qubits=False, \n",
    "         skip_final_rotation_layer=False, \n",
    "         parameter_prefix='θ', insert_barriers=False, \n",
    "         initial_state=None)\n",
    "```\n",
    "\n",
    "```python\n",
    "NLocal(num_qubits=None, reps=1, rotation_blocks=None, \n",
    "       entanglement_blocks=None, entanglement=None,   \n",
    "       skip_unentangled_qubits=False, \n",
    "       skip_final_rotation_layer=False, \n",
    "       overwrite_block_parameters=True, \n",
    "       parameter_prefix='θ', insert_barriers=False, \n",
    "       initial_state=None, name='nlocal')\n",
    "```\n",
    "\n",
    "Both functions create parameterised circuits of alternating rotation and entanglement layers. In both layers, parameterised circuit-blocks act on the circuit in a defined way. In the rotation layer, the blocks are applied stacked on top of each other, while in the entanglement layer according to the entanglement strategy. Each layer is repeated a number of times, and by default a final rotation layer is appended.\n",
    "\n",
    "In `NLocal`, the circuit blocks can have arbitrary sizes (smaller equal to the number of qubits in the circuit), while in `TwoLocal`, the rotation layers are single qubit gates applied on all qubits and the entanglement layer uses two-qubit gates.\n",
    "\n",
    "For example, here is a `TwoLocal` circuit, with $R_y$ and $R_Z$ gates in the rotation layer and $CX$ gates in the entangling layer with circular entanglement:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "twolocal = TwoLocal(num_qubits=3, reps=2, rotation_blocks=['ry','rz'], \n",
    "               entanglement_blocks='cx', entanglement='circular', insert_barriers=True)\n",
    "twolocal.decompose().draw('mpl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and the equivalent `NLocal` circuit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "twolocaln = NLocal(num_qubits=3, reps=2,\n",
    "               rotation_blocks=[RYGate(Parameter('a')), RZGate(Parameter('a'))], \n",
    "               entanglement_blocks=CXGate(), \n",
    "               entanglement='circular', insert_barriers=True)\n",
    "twolocaln.decompose().draw('mpl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's encode the first training sample using the `PauliFeatureMap`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f'First training data: {sample_train[0]}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "encode_map = PauliFeatureMap(feature_dimension=N_DIM, reps=1, paulis = ['X', 'Y', 'ZZ'])\n",
    "encode_circuit = encode_map.bind_parameters(sample_train[0])\n",
    "encode_circuit.decompose().draw(output='mpl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Challenge 3a**\n",
    "\n",
    "Construct a feature map to encode a 5-dimensionally embedded data, using 'ZZFeatureMap' with 3 repetitions, 'circular' entanglement and the rest as default.\n",
    "    \n",
    "</div>\n",
    "\n",
    "Submission format:\n",
    "```python\n",
    "ex3a_fmap = ZZFeatureMap(...)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##############################\n",
    "# Provide your code here\n",
    "\n",
    "\n",
    "ex3a_fmap = \n",
    "\n",
    "\n",
    "##############################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check your answer and submit using the following code\n",
    "from qc_grader import grade_ex3a\n",
    "grade_ex3a(ex3a_fmap)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Quantum Kernel Estimation\n",
    "\n",
    "A quantum feature map, $\\phi(\\mathbf{x})$, naturally gives rise to a quantum kernel, $k(\\mathbf{x}_i,\\mathbf{x}_j)= \\phi(\\mathbf{x}_j)^\\dagger\\phi(\\mathbf{x}_i)$, which can be seen as a measure of similarity: $k(\\mathbf{x}_i,\\mathbf{x}_j)$ is large when $\\mathbf{x}_i$ and $\\mathbf{x}_j$ are close. \n",
    "\n",
    "When considering finite data, we can represent the quantum kernel as a matrix: \n",
    "$K_{ij} = \\left| \\langle \\phi^\\dagger(\\mathbf{x}_j)| \\phi(\\mathbf{x}_i) \\rangle \\right|^{2}$. We can calculate each element of this kernel matrix on a quantum computer by calculating the transition amplitude:\n",
    "$$\n",
    "\\left| \\langle \\phi^\\dagger(\\mathbf{x}_j)| \\phi(\\mathbf{x}_i) \\rangle \\right|^{2} = \n",
    "\\left| \\langle 0^{\\otimes n} | \\mathbf{U_\\phi^\\dagger}(\\mathbf{x}_j) \\mathbf{U_\\phi}(\\mathbf{x_i}) | 0^{\\otimes n} \\rangle \\right|^{2}\n",
    "$$\n",
    "assuming the feature map is a parameterized quantum circuit, which can be described as a unitary transformation $\\mathbf{U_\\phi}(\\mathbf{x})$ on $n$ qubits. \n",
    "\n",
    "This provides us with an estimate of the quantum kernel matrix, which we can then use in a kernel machine learning algorithm, such as support vector classification.\n",
    "\n",
    "As discussed in [***Havlicek et al*.  Nature 567, 209-212 (2019)**](https://www.nature.com/articles/s41586-019-0980-2), quantum kernel machine algorithms only have the potential of quantum advantage over classical approaches if the corresponding quantum kernel is hard to estimate classically. \n",
    "\n",
    "As we will see later, the hardness of estimating the kernel with classical resources is of course only a necessary and not always sufficient condition to obtain a quantum advantage. \n",
    "\n",
    "However, it was proven recently in [***Liu et al.* arXiv:2010.02174 (2020)**](https://arxiv.org/abs/2010.02174) that learning problems exist for which learners with access to quantum kernel methods have a quantum advantage over all classical learners.\n",
    "\n",
    "With our training and testing datasets ready, we set up the `QuantumKernel` class with the PauliFeatureMap, and use the `BasicAer` `statevector_simulator` to estimate the training and testing kernel matrices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pauli_map = PauliFeatureMap(feature_dimension=N_DIM, reps=1, paulis = ['X', 'Y', 'ZZ'])\n",
    "pauli_kernel = QuantumKernel(feature_map=pauli_map, quantum_instance=Aer.get_backend('statevector_simulator'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's calculate the transition amplitude between the first and second training data samples, one of the entries in the training kernel matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f'First training data : {sample_train[0]}')\n",
    "print(f'Second training data: {sample_train[1]}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we create and draw the circuit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pauli_circuit = pauli_kernel.construct_circuit(sample_train[0], sample_train[1])\n",
    "pauli_circuit.decompose().decompose().draw(output='mpl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The parameters in the gates are a little difficult to read, but notice how the circuit is symmetrical, with one half encoding one of the data samples, the other half encoding the other. \n",
    "\n",
    "We then simulate the circuit. We will use the `qasm_simulator` since the circuit contains measurements, but increase the number of shots to reduce the effect of sampling noise. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "backend = Aer.get_backend('qasm_simulator')\n",
    "job = execute(pauli_circuit, backend, shots=8192, \n",
    "              seed_simulator=1024, seed_transpiler=1024)\n",
    "counts = job.result().get_counts(pauli_circuit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The transition amplitude is the proportion of counts in the zero state:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Transition amplitude: {counts['0'*N_DIM]/sum(counts.values())}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This process is then repeated for each pair of training data samples to fill in the training kernel matrix, and between each training and testing data sample to fill in the testing kernel matrix. Note that each matrix is symmetric, so to reduce computation time, only half the entries are calculated explicitly. \n",
    "\n",
    "Here we compute and plot the training and testing kernel matrices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "matrix_train = pauli_kernel.evaluate(x_vec=sample_train)\n",
    "matrix_val = pauli_kernel.evaluate(x_vec=sample_val, y_vec=sample_train)\n",
    "\n",
    "fig, axs = plt.subplots(1, 2, figsize=(10, 5))\n",
    "axs[0].imshow(np.asmatrix(matrix_train),\n",
    "              interpolation='nearest', origin='upper', cmap='Blues')\n",
    "axs[0].set_title(\"training kernel matrix\")\n",
    "axs[1].imshow(np.asmatrix(matrix_val),\n",
    "              interpolation='nearest', origin='upper', cmap='Reds')\n",
    "axs[1].set_title(\"validation kernel matrix\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "</div>\n",
    "    \n",
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Challenge 3b**\n",
    "\n",
    "Calculate the transition amplitude between $x = (-0.5, -0.4, 0.3, 0, -0.9)$ and $y = (0, -0.7, -0.3, 0, -0.4)$ using the 'ZZFeatureMap' with 3 repetitions, 'circular' entanglement and the rest as default. Use the 'qasm_simulator' with 'shots=8192', 'seed_simulator=1024' and 'seed_transpiler=1024'.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = [-0.5, -0.4, 0.3, 0, -0.9]\n",
    "y = [0, -0.7, -0.3, 0, -0.4]\n",
    "\n",
    "##############################\n",
    "# Provide your code here\n",
    "\n",
    "\n",
    "ex3b_amp = \n",
    "\n",
    "\n",
    "##############################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check your answer and submit using the following code\n",
    "from qc_grader import grade_ex3b\n",
    "grade_ex3b(ex3b_amp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Related QGSS materials:\n",
    "- [**Kernel Trick (Lecture 6.1)**](https://www.youtube.com/watch?v=m6EzmYsEOiI&list=PLOFEBzvs-VvqJwybFxkTiDzhf5E11p8BI&index=14)\n",
    "- [**Kernel Trick (Lecture 6.2)**](https://www.youtube.com/watch?v=zw3JYUrS-v8&list=PLOFEBzvs-VvqJwybFxkTiDzhf5E11p8BI&index=15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 Quantum Support Vector Machine (QSVM)\n",
    "\n",
    "Introduced in [***Havlicek et al*.  Nature 567, 209-212 (2019)**](https://www.nature.com/articles/s41586-019-0980-2), the quantum kernel support vector classification algorithm consists of these steps:\n",
    "\n",
    "\n",
    "<center><img src=\"./resources/qsvc.png\" width=\"1000\"></center> \n",
    "\n",
    "1. Build the train and test quantum kernel matrices.\n",
    "    1. For each pair of datapoints in the training dataset $\\mathbf{x}_{i},\\mathbf{x}_j$, apply the feature map and measure the transition probability: $ K_{ij} = \\left| \\langle 0 | \\mathbf{U}^\\dagger_{\\Phi(\\mathbf{x_j})} \\mathbf{U}_{\\Phi(\\mathbf{x_i})} | 0 \\rangle \\right|^2 $.\n",
    "    2. For each training datapoint $\\mathbf{x_i}$ and testing point $\\mathbf{y_j}$, apply the feature map and measure the transition probability: $ K_{ij} = \\left| \\langle 0 | \\mathbf{U}^\\dagger_{\\Phi(\\mathbf{y_j})} \\mathbf{U}_{\\Phi(\\mathbf{x_i})} | 0 \\rangle \\right|^2 $.\n",
    "2. Use the train and test quantum kernel matrices in a classical support vector machine classification algorithm.\n",
    "\n",
    "The `scikit-learn` `svc` algorithm allows us to [**define a custom kernel**](https://scikit-learn.org/stable/modules/svm.html#custom-kernels) in two ways: by providing the kernel as a callable function or by precomputing the kernel matrix. We can do either of these using the `QuantumKernel` class in Qiskit.\n",
    "\n",
    "The following code takes the training and testing kernel matrices we calculated earlier and  provides them to the `scikit-learn` `svc` algorithm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pauli_svc = SVC(kernel='precomputed')\n",
    "pauli_svc.fit(matrix_train, labels_train)\n",
    "pauli_score = pauli_svc.score(matrix_val, labels_val)\n",
    "\n",
    "print(f'Precomputed kernel classification test score: {pauli_score*100}%')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Related QGSS materials:\n",
    "- [**Classical SVM (Lecture 4.2)**](https://www.youtube.com/watch?v=lpPij21jnZ4&list=PLOFEBzvs-VvqJwybFxkTiDzhf5E11p8BI&index=9)\n",
    "- [**Quantum Classifier (Lecture 5.1)**](https://www.youtube.com/watch?v=-sxlXNz7ZxU&list=PLOFEBzvs-VvqJwybFxkTiDzhf5E11p8BI&index=11)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2: Challenge - QSVM for 3-class classification of Fashion-MNIST\n",
    "\n",
    "In this part, you will use what your have learned so far to implement 3-class classification of clothing images and work on improving its accuracy. \n",
    "    \n",
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Challenge 3c**\n",
    "\n",
    "**Goal**: Implement a 3-class classifier using QSVM and achieve 70% accuracy on clothing image dataset with smaller feature map circuits.\n",
    "\n",
    "**Dataset**: Fashion-MNIST clothing image dataset. There are following three dataset in this challnge.  \n",
    "- Train: Both images and labels are given.\n",
    "- Public test: Images are given and labels are hidden.\n",
    "- Private test: Both images and labels are hidden.\n",
    "    \n",
    "Grading will be performed on both public test and private test data. The purpose of this is to make sure that quantum methods are used, so that cheating is not possible.\n",
    "    \n",
    "</div>\n",
    "\n",
    "### How to implement a multi-class classifier using binary classifiers\n",
    "\n",
    "So far, you have learned how to implement binary classification with QSVM. Now, how can you scale it up to multi-class classification? There are two approaches to do so. One is the One-vs-Rest approach, and the other is the One-vs-One approach.\n",
    "\n",
    "1. One-vs-Rest: In this approach, multi-class classification is achieved by combining classifiers for each class that classifies the class as positive and the others as negative. Since one classifier is required for each class, the total number of classifiers required for N-class classification is N. The advantage is that fewer classifiers are needed, and the disadvantage is that the labels are likely to be imbalanced in each classification.\n",
    "2. One-vs-One: In this approach, multi-class classification is achieved by combining classifiers for each pair of two classes, where one is positive and the other is negative. Since one classifier is required for each label pair, the total number of classifiers required for N-class classification is N(N-1)/2. The advantage is that labels are less likely to be imbalanced in each classification, and the disadvantage is that the number of classifiers required is larger.\n",
    "\n",
    "Both approaches can be used to solve this problem, but here you will be given hints based on the One-vs-Rest approach. Please follow the hints to solve it.\n",
    "\n",
    "<center><img src=\"./resources/onevsrest.png\" width=\"800\"></center>\n",
    "\n",
    "Figure via [cc.gatech.edu](https://www.cc.gatech.edu/classes/AY2016/cs4476_fall/results/proj4/html/jnanda3/index.html)\n",
    "\n",
    "### 1. Data preparation\n",
    "The data we are working with here is a small subset of clothing image dataset called Fashion-MNIST, which is a variant of the MNIST dataset. We aim to classify the following labels.\n",
    "- label 0: T-shirt/top\n",
    "- label 2: pullover\n",
    "- label 3: dress\n",
    "\n",
    "First, let's load the dataset and display one image for each class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load MNIST dataset\n",
    "DATA_PATH = './resources/ch3_part2.npz'\n",
    "data = np.load(DATA_PATH)\n",
    "\n",
    "sample_train = data['sample_train']\n",
    "labels_train = data['labels_train']\n",
    "sample_test = data['sample_test']\n",
    "\n",
    "# Split train data\n",
    "sample_train, sample_val, labels_train, labels_val = train_test_split(\n",
    "    sample_train, labels_train, test_size=0.2, random_state=42)\n",
    "\n",
    "# Visualize samples\n",
    "fig = plt.figure()\n",
    "\n",
    "LABELS = [0, 2, 3]\n",
    "num_labels = len(LABELS)\n",
    "for i in range(num_labels):\n",
    "    ax = fig.add_subplot(1, num_labels, i+1)\n",
    "    img = sample_train[labels_train==LABELS[i]][0].reshape((28, 28))\n",
    "    ax.imshow(img, cmap=\"Greys\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, preprocess the dataset in the same way as before.\n",
    "- Standardization\n",
    "- PCA\n",
    "- Normalization\n",
    "\n",
    "Note that you can change the number of features here by changing N_DIM."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Standardize\n",
    "standard_scaler = StandardScaler()\n",
    "sample_train = standard_scaler.fit_transform(sample_train)\n",
    "sample_val = standard_scaler.transform(sample_val)\n",
    "sample_test = standard_scaler.transform(sample_test)\n",
    "\n",
    "# Reduce dimensions\n",
    "N_DIM = 5\n",
    "pca = PCA(n_components=N_DIM)\n",
    "sample_train = pca.fit_transform(sample_train)\n",
    "sample_val = pca.transform(sample_val)\n",
    "sample_test = pca.transform(sample_test)\n",
    "\n",
    "# Normalize\n",
    "min_max_scaler = MinMaxScaler((-1, 1))\n",
    "sample_train = min_max_scaler.fit_transform(sample_train)\n",
    "sample_val = min_max_scaler.transform(sample_val)\n",
    "sample_test = min_max_scaler.transform(sample_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Modeling\n",
    "Based on the One-vs-Rest approach, you need to create the following three QSVM binary classifiers\n",
    "- the label 0 and the rest\n",
    "- the label 2 and the rest\n",
    "- the label 3 and the rest\n",
    "\n",
    "Here is the first one as a hint.\n",
    "\n",
    "### 2.1: Label 0 vs Rest\n",
    "Create new labels with label 0 as positive(1) and the rest as negative(0) as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "labels_train_0 = np.where(labels_train==0, 1, 0)\n",
    "labels_val_0 = np.where(labels_val==0, 1, 0)\n",
    "\n",
    "print(f'Original validation labels:      {labels_val}')\n",
    "print(f'Validation labels for 0 vs Rest: {labels_val_0}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See only places where the original label was 0 are set to 1.  \n",
    "\n",
    "Next, construct a binary classifier using QSVM as before.  \n",
    "Note that PauliFeatureMap is used in this hint but you can use a different feature map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pauli_map_0 = PauliFeatureMap(feature_dimension=N_DIM, reps=2, paulis = ['X', 'Y', 'ZZ'])\n",
    "pauli_kernel_0 = QuantumKernel(feature_map=pauli_map_0, quantum_instance=Aer.get_backend('statevector_simulator'))\n",
    "\n",
    "pauli_svc_0 = SVC(kernel='precomputed', probability=True)\n",
    "\n",
    "matrix_train_0 = pauli_kernel_0.evaluate(x_vec=sample_train)\n",
    "pauli_svc_0.fit(matrix_train_0, labels_train_0)\n",
    "\n",
    "matrix_val_0 = pauli_kernel_0.evaluate(x_vec=sample_val, y_vec=sample_train)\n",
    "pauli_score_0 = pauli_svc_0.score(matrix_val_0, labels_val_0)\n",
    "print(f'Accuracy of discriminating between label 0 and others: {pauli_score_0*100}%')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see that the QSVM binary classifier is able to distinguish between label 0 and the rest with a reasonable probability.\n",
    "\n",
    "Finally, for each of the test data, calculate the probability that it has label 0. It can be obtained by ```predict_proba``` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "matrix_test_0 = pauli_kernel_0.evaluate(x_vec=sample_test, y_vec=sample_train)\n",
    "pred_0 = pauli_svc_0.predict_proba(matrix_test_0)[:, 1]\n",
    "print(f'Probability of label 0: {np.round(pred_0, 2)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These probabilities are important clues for multiclass classification.  \n",
    "Obtain the probabilities for the remaining two labels in the same way.\n",
    "\n",
    "### 2.2: Label 2 vs Rest\n",
    "Build a binary classifier using QSVM and get the probability of label 2 for test dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##############################\n",
    "# Provide your code here\n",
    "\n",
    "\n",
    "pred_2 = \n",
    "\n",
    "\n",
    "##############################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 Label 3 vs Rest\n",
    "Build a binary classifier using QSVM and get the probability of label 3 for test dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##############################\n",
    "# Provide your code here\n",
    "\n",
    "\n",
    "pred_3 = \n",
    "\n",
    "\n",
    "##############################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Prediction\n",
    "Lastly, make a final prediction based on the probability of each label.  \n",
    "The prediction you submit should be in the following format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_pred = np.load('./resources/ch3_part2_sub.npy')\n",
    "print(f'Sample prediction: {sample_pred}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to understand the method to make predictions for multiclass classification, let's begin with the case of making predictions for just two labels, label 2 and label 3.\n",
    "\n",
    "If probabilities are as follows for a certain data, label 2 should be considered the most plausible.\n",
    "- probability of label 2: 0.7\n",
    "- probability of label 3: 0.2\n",
    "\n",
    "You can implement this with ```np.where``` function. (Of course, you can use different methods.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred_2_ex = np.array([0.7])\n",
    "pred_3_ex = np.array([0.2])\n",
    "\n",
    "pred_test_ex = np.where((pred_2_ex > pred_3_ex), 2, 3)\n",
    "print(f'Prediction: {pred_test_ex}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can apply this method as is to multiple data.\n",
    "\n",
    "If second data has probabilities for each label as follows, it should be classified as label 3.\n",
    "- probability of label 2: 0.1\n",
    "- probability of label 3: 0.6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred_2_ex = np.array([0.7, 0.1])\n",
    "pred_3_ex = np.array([0.2, 0.6])\n",
    "\n",
    "pred_test_ex = np.where((pred_2_ex > pred_3_ex), 2, 3)\n",
    "print(f'Prediction: {pred_test_ex}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This method can be extended to make predictions for 3-class classification.\n",
    "\n",
    "Implement such an extended method and make the final 3-class predictions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##############################\n",
    "# Provide your code here\n",
    "\n",
    "\n",
    "pred_test = \n",
    "\n",
    "\n",
    "##############################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Submission\n",
    "    \n",
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Challenge 3c**\n",
    "\n",
    "**Submission**: Submit the following 11 items.\n",
    "- **pred_test**: prediction for the public test dataset\n",
    "- **sample_train**: train data used to obtain kernels\n",
    "- **standard_scaler**: the one used to standardize data\n",
    "- **pca**: the one used to reduce dimention\n",
    "- **min_max_scaler**: the one used to normalize data\n",
    "- **kernel_0**: the kernel for the \"label 0 vs rest\" classifier\n",
    "- **kernel_2**: the kernel for the \"label 2 vs rest\" classifier\n",
    "- **kernel_3**: the kernel for the \"label 3 vs rest\" classifier\n",
    "- **svc_0**: the SVC trained to classify \"label 0 vs rest\"\n",
    "- **svc_2**: the SVC trained to classify \"label 2 vs rest\"\n",
    "- **svc_3**: the SVC trained to classify \"label 3 vs rest\"\n",
    "\n",
    "**Criteria**: Accuracy of 70% or better on both public and private test data.\n",
    "\n",
    "**Score**: Solutions that pass the criteria will be scored as follows. The smaller this final score is, the better.\n",
    "1. Each feature map gets transpiled with:\n",
    "    - basis_gates=['u1', 'u2', 'u3', 'cx']\n",
    "    - optimization_level=0\n",
    "2. Calculate the cost for each transpiled circuit:   \n",
    "   cost = 10 * #cx + (#u1 + #u2 + #u3)\n",
    "3. The sum of the costs will be the final score.\n",
    "\n",
    "</div>\n",
    "\n",
    "Again, the prediction you submit should be in the following format.\n",
    "- prediction for the public test data (**sample_test**)\n",
    "- type: numpy.ndarray\n",
    "- shape: (20,)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f'Sample prediction: {sample_pred}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check your answer and submit using the following code\n",
    "from qc_grader import grade_ex3c\n",
    "grade_ex3c(pred_test, sample_train, \n",
    "           standard_scaler, pca, min_max_scaler,\n",
    "           kernel_0, kernel_2, kernel_3,\n",
    "           svc_0, svc_2, svc_3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Additional information\n",
    "\n",
    "**Created by:** Shota Nakasuji, Anna Phan\n",
    "\n",
    "**Version:** 1.0.0"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}