ClassificationTask.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9a21b8df",
   "metadata": {},
   "source": [
    "# Classifying the iris dataset on Quandela's QPU Ascella"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4bb04ecc",
   "metadata": {},
   "source": [
    "_This notebook was developed by Dario A. Fioretto and Alexia Salavrakos._"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29333a71",
   "metadata": {},
   "source": [
    "In this notebook, we train a model to classify the [iris toy dataset](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.load_iris.html) using Ascella. We use a variational quantum algorithm (or quantum neural network) as explained in the paper.\n",
    "\n",
    "\n",
    "As seen in the image below, we use modes 3 to 7 of Ascella for our ansatz. We input three photons in modes 3, 5 and 7 respectively. In a first block, the variational parameters $\\theta$ are encoded by phase-shifters. In a second block, the features of the data $x$ are encoded. We then have a third block with more variational parameters $\\theta$. Then, we have a block that implements pseudo-PNR on detectors 3 and 7 using modes 1,2 and 8,9 respectively. The lambda parameters are weights that multiply the output probabilities in order to define the observable and build the model."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3ed4687",
   "metadata": {},
   "source": [
    "![alternative text](images/ansatz_classifier.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de58f698",
   "metadata": {},
   "source": [
    "## Imports and configurations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9addd053",
   "metadata": {},
   "outputs": [],
   "source": [
    "import perceval as pcvl\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import skopt as sk\n",
    "\n",
    "import scipy\n",
    "from scipy.stats import poisson\n",
    "\n",
    "from sklearn import datasets\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import accuracy_score\n",
    "from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay\n",
    "\n",
    "from datetime import datetime\n",
    "import os\n",
    "import time\n",
    "import pandas as pd\n",
    "from itertools import combinations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ecf9b13d-b6f8-4bc9-9dc5-20aa4520cf0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "pcvl.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7030f84a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# to display the circuit\n",
    "from perceval.rendering.circuit import DisplayConfig, SymbSkin\n",
    "DisplayConfig.select_skin(SymbSkin)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c625fc8b",
   "metadata": {},
   "source": [
    "Insert token from https://cloud.quandela.com/ here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58c90cae",
   "metadata": {},
   "outputs": [],
   "source": [
    "qpu_token = ''"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70516243",
   "metadata": {},
   "source": [
    "Choose whether to run a noisy local simulation or a real experiment on Ascella:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8d0710e",
   "metadata": {},
   "outputs": [],
   "source": [
    "run_on_qpu = True"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0539834e",
   "metadata": {},
   "source": [
    "## Define quantum classifier class"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28b357a2",
   "metadata": {},
   "source": [
    "We define a class that contains all the relevant functions as methods: building the ansatz, computing the outputs of the experiment, implementing pseudo-PNR and its mapping, as well as building the classifier estimator from the output probabilities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09c37fc7",
   "metadata": {},
   "outputs": [],
   "source": [
    "class QuantumClassifierOnAscella():\n",
    "    def __init__(self, dataset, n_classes, run_on_qpu):\n",
    "        self.dataset = dataset\n",
    "        self.n_classes = n_classes\n",
    "        \n",
    "        # For a test on a very small dataset, use:\n",
    "        #self.small_set_indices = [0,1,2, 51, 52, 53, 101, 102, 103]\n",
    "        #self.X = np.array([self.dataset.data[i] for i in self.small_set_indices])\n",
    "        #self.Y = np.array([self.dataset.target[i] for i in self.small_set_indices])\n",
    "        \n",
    "        # Prepare data \n",
    "        self.X = self.dataset.data\n",
    "        self.Y = self.dataset.target\n",
    "        self.X_angle = self.transform_features_to_phases()\n",
    "        self.X_train, self.X_test, self.y_train, self.y_test = self.split_dataset(self.X_angle, self.Y)\n",
    "        \n",
    "        # Connect to the remote processor and get Ascella circuit\n",
    "        self.processor = pcvl.RemoteProcessor(\"qpu:ascella\", qpu_token)\n",
    "        self.circuit = self.processor.specs[\"specific_circuit\"]\n",
    "        print(\"Received circuit from Ascella.\")\n",
    "        self.circuit_phases = self.circuit.get_parameters()\n",
    "        \n",
    "        # If running a simulation, change processor to a local backend with noisy source\n",
    "        if not run_on_qpu:\n",
    "            self.source = pcvl.Source(losses = 0.92, emission_probability=1, multiphoton_component=0,\n",
    "                                      indistinguishability=0.92)\n",
    "            self.processor = pcvl.Processor(\"SLOS\", self.circuit, self.source)\n",
    "        \n",
    "        # Postselect on observing 3 photon counts (discard photon losses)\n",
    "        self.processor.set_parameter(\"mode_post_select\", 3)\n",
    "        \n",
    "        # Define chip phases and set values of fixed phases\n",
    "        self.theta_phases = ['phi2', 'phi4', 'phi3', 'phi5', 'phi13',\n",
    "                             'phi14', 'phi15', 'phi16', 'phi17', 'phi24',\n",
    "                             'phi25', 'phi26', 'phi27', 'phi35', 'phi36', 'phi38',  # end block 1, 16 params\n",
    "                             'phi57', 'phi58', 'phi59', 'phi60', 'phi61', \n",
    "                             'phi68', 'phi69', 'phi70', 'phi71', 'phi79',\n",
    "                             'phi80', 'phi81', 'phi82', 'phi83', 'phi90', 'phi92'] # end block 2, 16 params\n",
    "        self.data_phases = ['phi37', 'phi39', 'phi47', 'phi49'] \n",
    "        self.set_fixed_phases()\n",
    "        \n",
    "        # Define expected output states and pseudo PNR mapping dictionary\n",
    "        self.possible_output_states = self.create_dict_output_states_with_pseudo_PNR(n_modes = 5, n_photons = 3,\n",
    "                                                                                     starting_mode = 2)\n",
    "        self.correspondance_dict = self.create_correspondance_pseudo_pnr()\n",
    "        \n",
    "        # Choose number of samples per experiment\n",
    "        self.n_samples = 5*1e4\n",
    "\n",
    "        # Define max number of shots if run on QPU, with some margin\n",
    "        if run_on_qpu:\n",
    "            # Define simple circuit to get the estimate\n",
    "            remote_qpu = pcvl.RemoteProcessor(\"qpu:ascella\", qpu_token)\n",
    "            remote_qpu.set_circuit(pcvl.Circuit(9))\n",
    "            remote_qpu.with_input(pcvl.BasicState([0, 0, 1, 0, 1, 0, 1, 0, 0]))\n",
    "            remote_qpu.min_detected_photons_filter(3)\n",
    "            # Get estimated required number of shots\n",
    "            required_shots = remote_qpu.estimate_required_shots(nsamples=self.n_samples)\n",
    "            # Assign with a margin factor of 10\n",
    "            self.required_shots = required_shots * 10\n",
    "        \n",
    "        # Define empty lists for the outcomes and probabilities and for the estimator \n",
    "        self.outcomes = []\n",
    "        self.probabilities = []\n",
    "        self.predicted_class = []\n",
    "        self.model_output = []\n",
    "        self.model_output_error = []\n",
    "        \n",
    "        # Create folder to save data\n",
    "        self.path = datetime.now().strftime(\"%Y%m%d_%H%M%S\") + '/'\n",
    "        os.mkdir(self.path)\n",
    "        \n",
    "        \n",
    "    def split_dataset(self, X, Y):\n",
    "        # Given the data X and the target Y, split in a training and test set\n",
    "        X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size=0.25, random_state=42)\n",
    "        return X_train, X_test, y_train, y_test\n",
    "\n",
    "    \n",
    "    def transform_features_to_phases(self):\n",
    "        # Given the features from data X, reformat them as phases to be encoded by phase shifters\n",
    "        a = []\n",
    "        b = []\n",
    "        for i in range(len(self.X[0])):\n",
    "            a.append((max(self.X[:, i]) - min(self.X[:, i])) / 2)\n",
    "            b.append(((max(self.X[:, i]) + min(self.X[:, i])) / 2))\n",
    "\n",
    "        def sigmoid_Pi(x):\n",
    "            return np.pi / (1 + np.exp(-x))\n",
    "\n",
    "        X_angle = sigmoid_Pi((self.X - b) / a)\n",
    "\n",
    "        return X_angle\n",
    "\n",
    "    \n",
    "    def get_ascella_fixed_phases_dict(self):\n",
    "        # In our approach without transpilation, we set some phases on the chip to pi, or pi/2\n",
    "        # The phases set to pi prevent photons from escaping to extra modes\n",
    "        # The phases set to pi/2 enable pseudo-PNR by redirecting photons to some specific modes\n",
    "        phases_set_to_pi = {'phi6': np.pi, 'phi12': np.pi, 'phi28': np.pi, 'phi34': np.pi, 'phi46': np.pi,\n",
    "                            'phi48': np.pi, 'phi50': np.pi, 'phi56': np.pi, 'phi72': np.pi, 'phi78': np.pi,\n",
    "                            'phi94': np.pi, 'phi102': np.pi, 'phi104': np.pi,\n",
    "                            'phi112': np.pi, 'phi114': np.pi, 'phi122': np.pi, \n",
    "                            'phi124': np.pi, 'phi126': np.pi}\n",
    "        phases_set_to_pi_half = {'phi100': np.pi/2, 'phi110': np.pi/2, 'phi116': np.pi/2, 'phi128': np.pi/2}\n",
    "        return ({**phases_set_to_pi, **phases_set_to_pi_half})\n",
    "\n",
    "    \n",
    "    def get_ascella_theta_phases_dict(self, theta_values):\n",
    "        # Create a dictionary that assigns a value to the phases corresponding to theta parameters\n",
    "        theta_phases_dict = dict((self.theta_phases[i], theta_values[i]) for i in range(len(self.theta_phases)))\n",
    "        return (theta_phases_dict)\n",
    "\n",
    "    \n",
    "    def get_ascella_data_phases_dict(self, data_values):\n",
    "        # Create a dictionary that assigns a value to the phases corresponding to the data encoding\n",
    "        data_phases_dict = dict((self.data_phases[i], data_values[i]) for i in range(len(self.data_phases)))\n",
    "        return (data_phases_dict)\n",
    "\n",
    "    \n",
    "    def set_fixed_phases(self):\n",
    "        # Create a dictionary for the fixed phases (set to 0, pi or pi/2)\n",
    "        # Then set all other phases to 0 for uniformity, with the exception of the data phases\n",
    "        # Some of these 0s will be replaced by theta parameters in set_unitaries_ascella\n",
    "        all_phases_ordered = [('phi' + str(i)) for i in range(132)]\n",
    "        fixed_phases_dict = self.get_ascella_fixed_phases_dict()\n",
    "        puppet_data_phases_dict = self.get_ascella_data_phases_dict([1 for _ in range(len(self.data_phases))])\n",
    "        \n",
    "        for (i, j) in enumerate(all_phases_ordered):\n",
    "            if j in fixed_phases_dict:\n",
    "                self.circuit_phases[i].set_value(fixed_phases_dict[j])\n",
    "            elif j in puppet_data_phases_dict:\n",
    "                # do nothing\n",
    "                # keep the data phases non assigned, as they will be sent in batch to the processor\n",
    "                pass\n",
    "            else:\n",
    "                self.circuit_phases[i].set_value(0)\n",
    "\n",
    "        return 0\n",
    "\n",
    "    \n",
    "    def set_unitaries_ascella(self, theta_values):\n",
    "        # Assign values to the theta parameter phases in the circuit\n",
    "        theta_phases_dict = self.get_ascella_theta_phases_dict(theta_values)\n",
    "        all_phases_ordered = [('phi' + str(i)) for i in range(132)]\n",
    "        for (i, j) in enumerate(all_phases_ordered):\n",
    "            if j in theta_phases_dict:\n",
    "                self.circuit_phases[i].set_value(theta_phases_dict[j])\n",
    "\n",
    "        return 0\n",
    "\n",
    "    \n",
    "    def create_dict_output_states(self, n_modes, n_photons, starting_mode = 0):\n",
    "        # Create a list with the possible output states\n",
    "        possible_states = []\n",
    "\n",
    "        for comb in combinations(np.arange(n_modes), n_photons):\n",
    "            vector_state = [1 if i in comb else 0 for i in range(n_modes)]\n",
    "            possible_states.append(pcvl.BasicState([0]*starting_mode + vector_state + \n",
    "                                                   [0]*(12 - n_modes - starting_mode)))\n",
    "\n",
    "        return possible_states\n",
    "\n",
    "\n",
    "    def create_dict_output_states_with_pseudo_PNR(self, n_modes, n_photons, starting_mode = 0):\n",
    "        # Create a list with the possible output states\n",
    "        possible_states = []\n",
    "\n",
    "        for comb in combinations(np.arange(n_modes), n_photons):\n",
    "            vector_state = [1 if i in comb else 0 for i in range(n_modes)]\n",
    "            possible_states.append(pcvl.BasicState([0]*starting_mode + vector_state + \n",
    "                                                   [0]*(12 - n_modes - starting_mode)))\n",
    "\n",
    "        # Here we hardcode the extra states from pseudo PNR:\n",
    "        # 3 photons detected in modes 3/7\n",
    "        possible_states.append(pcvl.BasicState([0]*2 + [3] + [0]*9))\n",
    "        possible_states.append(pcvl.BasicState([0]*6 + [3] + [0]*5))\n",
    "        \n",
    "        # 2 photons detected in modes 3 + 1 photon elsewhere\n",
    "        for comb in combinations(np.arange(4), 1):\n",
    "            vector_comb = [1 if i in comb else 0 for i in range(4)]\n",
    "            possible_states.append(pcvl.BasicState([0]*2 + [2] + vector_comb + [0]*5))\n",
    "        \n",
    "        # 2 photons detected in modes 7 + 1 photon elsewhere\n",
    "        for comb in combinations(np.arange(4), 1):\n",
    "            vector_comb = [1 if i in comb else 0 for i in range(4)]\n",
    "            possible_states.append(pcvl.BasicState([0]*2 + vector_comb + [2] + [0]*5))\n",
    "\n",
    "        return possible_states\n",
    "    \n",
    "    \n",
    "    def create_correspondance_pseudo_pnr(self):\n",
    "        # Prepare a correspondance dictionary for the pseudo PNR\n",
    "        # The keys are the raw state and the values are the mapped pseudo pnr state \n",
    "        # This allows us to do the coarse graining or reinterpretation of the observed states as PNR states\n",
    "        correspondance_dict = {}\n",
    "        \n",
    "        # \"Normal\" 3 photon coincidences between modes 3 and 7\n",
    "        states_non_PNR = self.create_dict_output_states(n_modes = 5, n_photons = 3, starting_mode = 2)\n",
    "        for state in states_non_PNR:\n",
    "            correspondance_dict[state] = state\n",
    "\n",
    "        # Hardcoded mapping for pseudo PNR\n",
    "        # Three photons in detectors 1,2,3 \n",
    "        # Same in 7,8,9\n",
    "        correspondance_dict[pcvl.BasicState([1, 1, 1] + [0]*9)] = pcvl.BasicState([0]*2 + [3] + [0]*9)\n",
    "        correspondance_dict[pcvl.BasicState([0]*6 + [1, 1, 1] + [0]*3)] = pcvl.BasicState([0]*6 + [3] + [0]*5)\n",
    "        \n",
    "        # Two photons in detectors 1,2,3 and one photon in detectors 4 to 7 \n",
    "        for comb in combinations(np.arange(4), 1):\n",
    "            vector_comb = [1 if i in comb else 0 for i in range(4)]\n",
    "            correspondance_dict[pcvl.BasicState([1, 1, 0] + vector_comb + [0]*5)] = pcvl.BasicState(\n",
    "                [0, 0, 2] + vector_comb + [0]*5)\n",
    "            correspondance_dict[pcvl.BasicState([1, 0, 1] + vector_comb + [0]*5)] = pcvl.BasicState(\n",
    "                [0, 0, 2] + vector_comb + [0]*5)\n",
    "            correspondance_dict[pcvl.BasicState([0, 1, 1] + vector_comb + [0]*5)] = pcvl.BasicState(\n",
    "                [0, 0, 2] + vector_comb + [0]*5)\n",
    "            \n",
    "        # Two photons in detectors 1,2,3 and one photon in detectors 8 or 9\n",
    "        correspondance_dict[pcvl.BasicState([1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0])] = pcvl.BasicState(\n",
    "            [0, 0, 2, 0, 0, 0, 1, 0, 0, 0, 0, 0])\n",
    "        correspondance_dict[pcvl.BasicState([1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0])] = pcvl.BasicState(\n",
    "            [0, 0, 2, 0, 0, 0, 1, 0, 0, 0, 0, 0])\n",
    "        correspondance_dict[pcvl.BasicState([1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0])] = pcvl.BasicState(\n",
    "            [0, 0, 2, 0, 0, 0, 1, 0, 0, 0, 0, 0])\n",
    "        correspondance_dict[pcvl.BasicState([1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0])] = pcvl.BasicState(\n",
    "            [0, 0, 2, 0, 0, 0, 1, 0, 0, 0, 0, 0])\n",
    "        correspondance_dict[pcvl.BasicState([0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0])] = pcvl.BasicState(\n",
    "            [0, 0, 2, 0, 0, 0, 1, 0, 0, 0, 0, 0])\n",
    "        correspondance_dict[pcvl.BasicState([0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0])] = pcvl.BasicState(\n",
    "            [0, 0, 2, 0, 0, 0, 1, 0, 0, 0, 0, 0])\n",
    "                                            \n",
    "        # Two photons in detectors 7,8,9 and one photon in detectors 3 to 6\n",
    "        for comb in combinations(np.arange(4), 1):\n",
    "            vector_comb = [1 if i in comb else 0 for i in range(4)]\n",
    "            correspondance_dict[pcvl.BasicState([0, 0] + vector_comb + [0, 1, 1] + [0]*3)] = pcvl.BasicState(\n",
    "                [0, 0] + vector_comb + [2, 0, 0] + [0]*3)\n",
    "            correspondance_dict[pcvl.BasicState([0, 0] + vector_comb + [1, 0, 1] + [0]*3)] = pcvl.BasicState(\n",
    "                [0, 0] + vector_comb + [2, 0, 0] + [0]*3)\n",
    "            correspondance_dict[pcvl.BasicState([0, 0] + vector_comb + [1, 1, 0] + [0]*3)] = pcvl.BasicState(\n",
    "                [0, 0] + vector_comb + [2, 0, 0] + [0]*3)\n",
    "            \n",
    "        # Two photons in detectors 7,8,9 and one photon in detectors 1 or 2\n",
    "        correspondance_dict[pcvl.BasicState([1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0])] = pcvl.BasicState(\n",
    "            [0, 0, 1, 0, 0, 0, 2, 0, 0, 0, 0, 0])\n",
    "        correspondance_dict[pcvl.BasicState([0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0])] = pcvl.BasicState(\n",
    "            [0, 0, 1, 0, 0, 0, 2, 0, 0, 0, 0, 0])\n",
    "        correspondance_dict[pcvl.BasicState([1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0])] = pcvl.BasicState(\n",
    "            [0, 0, 1, 0, 0, 0, 2, 0, 0, 0, 0, 0])\n",
    "        correspondance_dict[pcvl.BasicState([0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0])] = pcvl.BasicState(\n",
    "            [0, 0, 1, 0, 0, 0, 2, 0, 0, 0, 0, 0])\n",
    "        correspondance_dict[pcvl.BasicState([1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0])] = pcvl.BasicState(\n",
    "            [0, 0, 1, 0, 0, 0, 2, 0, 0, 0, 0, 0])\n",
    "        correspondance_dict[pcvl.BasicState([0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0])] = pcvl.BasicState(\n",
    "            [0, 0, 1, 0, 0, 0, 2, 0, 0, 0, 0, 0])\n",
    "\n",
    "        # One photon in detector 1,2 and 2 photons in detectors 4,5,6,7\n",
    "        for comb in combinations(np.arange(4), 2):\n",
    "            vector_comb = [1 if i in comb else 0 for i in range(4)]\n",
    "            correspondance_dict[pcvl.BasicState([1, 0, 0] + vector_comb + [0]*5)] = pcvl.BasicState(\n",
    "                [0, 0, 1] + vector_comb + [0]*5)\n",
    "            correspondance_dict[pcvl.BasicState([0, 1, 0] + vector_comb + [0]*5)] = pcvl.BasicState(\n",
    "                [0, 0, 1] + vector_comb + [0]*5) \n",
    "            \n",
    "        # One photon in detector 8,9 and 2 photons in detectors 3,4,5,6\n",
    "        for comb in combinations(np.arange(4), 2):\n",
    "            vector_comb = [1 if i in comb else 0 for i in range(4)]\n",
    "            correspondance_dict[pcvl.BasicState([0, 0] + vector_comb + [0, 1, 0] + [0]*3)] = pcvl.BasicState(\n",
    "                [0, 0] + vector_comb + [1, 0, 0] + [0]*3)\n",
    "            correspondance_dict[pcvl.BasicState([0, 0] + vector_comb + [0, 0, 1] + [0]*3)] = pcvl.BasicState(\n",
    "                [0, 0] + vector_comb + [1, 0, 0] + [0]*3)\n",
    "            \n",
    "        \n",
    "        # Finally, 1 photon in 1,2 and 1 photon in 4,5,6 and 1 photon in 8,9\n",
    "        for comb in combinations(np.arange(3), 1):\n",
    "            vector_comb = [1 if i in comb else 0 for i in range(3)]\n",
    "            correspondance_dict[pcvl.BasicState([1, 0, 0] + vector_comb + [0, 1, 0] + [0]*3)] = pcvl.BasicState(\n",
    "                [0, 0, 1] + vector_comb + [1, 0, 0] + [0]*3)\n",
    "            correspondance_dict[pcvl.BasicState([1, 0, 0] + vector_comb + [0, 0, 1] + [0]*3)] = pcvl.BasicState(\n",
    "                [0, 0, 1] + vector_comb + [1, 0, 0] + [0]*3)\n",
    "            correspondance_dict[pcvl.BasicState([0, 1, 0] + vector_comb + [0, 1, 0] + [0]*3)] = pcvl.BasicState(\n",
    "                [0, 0, 1] + vector_comb + [1, 0, 0] + [0]*3)\n",
    "            correspondance_dict[pcvl.BasicState([0, 1, 0] + vector_comb + [0, 0, 1] + [0]*3)] = pcvl.BasicState(\n",
    "                [0, 0, 1] + vector_comb + [1, 0, 0] + [0]*3)\n",
    "        \n",
    "        return correspondance_dict\n",
    "        \n",
    "    \n",
    "    def map_pseudo_PNR_states(self, raw_results):\n",
    "        # Given a list of raw results from the experiment with pseudo PNR\n",
    "        # Return reinterpreted results\n",
    "        # For example: 111000000000 is sent to 003000000000\n",
    "        results_PNR = {}\n",
    "        \n",
    "        for raw_state, pnr_state in self.correspondance_dict.items():\n",
    "            if pnr_state in results_PNR.keys():\n",
    "                previous_count = results_PNR[pnr_state]\n",
    "                results_PNR[pnr_state] = previous_count + raw_results[raw_state]\n",
    "            else:\n",
    "                results_PNR[pnr_state] = raw_results[raw_state]\n",
    "\n",
    "        return results_PNR\n",
    "    \n",
    "    \n",
    "    def prepare_batch_data_phases(self, X_train):\n",
    "        # Given a list of data points to be sent in batch, create a list of dictionaries \n",
    "        # Each dictionary contains the phases (e.g. phi_37) and the corresponding values (e.g. 1)\n",
    "        list_data_phases = []\n",
    "        for x in X_train:\n",
    "            data_phases_dict = self.get_ascella_data_phases_dict(x)\n",
    "            list_data_phases.append(data_phases_dict)\n",
    "        \n",
    "        return list_data_phases\n",
    "\n",
    "    \n",
    "    def get_results_ascella_batch(self, list_data_phases):\n",
    "        # Given a batch of data points, run an experiment on Ascella\n",
    "\n",
    "        # Define initial state: 3 photons\n",
    "        initial_state = pcvl.BasicState([0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0])\n",
    "\n",
    "        # Send computation to remote processor\n",
    "        self.processor.set_circuit(self.circuit)\n",
    "        self.processor.with_input(initial_state)\n",
    "        self.processor.set_parameter(\"parameter_iterator\", list_data_phases)\n",
    "        \n",
    "        if run_on_qpu:\n",
    "            sampler = pcvl.algorithm.Sampler(self.processor, max_shots_per_call=self.required_shots)\n",
    "        else:\n",
    "            sampler = pcvl.algorithm.Sampler(self.processor)\n",
    "            \n",
    "        remote_job = sampler.sample_count\n",
    "        remote_job.name = \"vqc_\" + datetime.now().strftime(\"%Y%m%d_%H%M%S\")\n",
    "        remote_job.execute_async(self.n_samples)\n",
    "\n",
    "        while not remote_job.is_complete:\n",
    "             time.sleep(0.02)\n",
    "\n",
    "        # Retrieve results and assign them to outcomes and probabilities\n",
    "        results = remote_job.get_results()\n",
    "        for j in range(len(results['results_list'])):\n",
    "            result_batch_item = self.map_pseudo_PNR_states(results['results_list'][j]['results'])\n",
    "            counts_list = []\n",
    "            for i in self.possible_output_states:\n",
    "                counts_list.append(result_batch_item[i])\n",
    "            counts_array = np.array(counts_list)\n",
    "            self.outcomes.append(counts_array)\n",
    "\n",
    "            probs = counts_array / np.sum(counts_array)\n",
    "            self.probabilities.append(probs)\n",
    "        return self.outcomes, self.probabilities\n",
    "    \n",
    "    \n",
    "    def get_results_ascella_all_data_batch(self, theta_values):\n",
    "        # For given theta parameters, set phases values, then prepare batch of data phases\n",
    "        # Then call get_results_ascella to run experiment on Ascella\n",
    "        self.set_unitaries_ascella(theta_values)\n",
    "        list_data_phases = self.prepare_batch_data_phases(self.X_train)\n",
    "        self.get_results_ascella(list_data_phases)\n",
    "        \n",
    "\n",
    "    def get_results_ascella(self, data_phases):\n",
    "        # Given a data point, run an experiment (without batch functionality)\n",
    "        \n",
    "        # set value to the right phase in the circuit\n",
    "        data_phases_dict = self.get_ascella_data_phases_dict(data_phases)\n",
    "        all_phases_ordered = [('phi' + str(i)) for i in range(132)]\n",
    "        for (i, j) in enumerate(all_phases_ordered):\n",
    "            if j in data_phases_dict:\n",
    "                self.circuit_phases[i].set_value(data_phases_dict[j])\n",
    "\n",
    "        # define initial state and send circuit to remote simulator\n",
    "        initial_state = pcvl.BasicState([0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0])\n",
    "\n",
    "        self.processor.set_circuit(self.circuit)\n",
    "        self.processor.with_input(initial_state)\n",
    "        \n",
    "        if run_on_qpu:\n",
    "            sampler = pcvl.algorithm.Sampler(self.processor, max_shots_per_call=self.required_shots)\n",
    "        else:\n",
    "            sampler = pcvl.algorithm.Sampler(self.processor)\n",
    "\n",
    "        remote_job = sampler.sample_count\n",
    "        remote_job.name = \"vqc_\" + datetime.now().strftime(\"%Y%m%d_%H%M%S\")\n",
    "        remote_job.execute_async(self.n_samples)\n",
    "\n",
    "        while not remote_job.is_complete:\n",
    "             time.sleep(0.02)\n",
    "\n",
    "        # Retrieve results and assign them to outcomes and probabilities\n",
    "        results = remote_job.get_results()\n",
    "        rr = self.map_pseudo_PNR_states(results['results'])\n",
    "\n",
    "        counts_list = []\n",
    "        for i in self.possible_output_states:\n",
    "            counts_list.append(rr[i])\n",
    "        counts_array = np.array(counts_list)\n",
    "        self.outcomes.append(counts_array)\n",
    "\n",
    "        probs = counts_array / np.sum(counts_array)\n",
    "        self.probabilities.append(probs)\n",
    "        return self.outcomes, self.probabilities\n",
    "    \n",
    "    \n",
    "    def get_results_ascella_all_data(self, theta_values):\n",
    "        # For given theta parameters, set phases values\n",
    "        # Then call get_results_ascella to run experiment on Ascella point by point (no batch)\n",
    "        self.set_unitaries_ascella(theta_values)\n",
    "        data_point_counter = 0\n",
    "        for x in self.X_train:\n",
    "            self.get_results_ascella(x)\n",
    "            data_point_counter = data_point_counter + 1\n",
    "            if (data_point_counter % 100 == 0):\n",
    "                print('100 datapoint processed')\n",
    "        \n",
    "    \n",
    "    def class_prediction(self, lambda_values):\n",
    "        # Given measurement parameters lambda,\n",
    "        # For the results of an experiment, obtain model prediction\n",
    "        thresholds = np.linspace(-1, 1, self.n_classes, False)\n",
    "        \n",
    "        for i in [proba*100000 for proba in self.probabilities]: # factor added for the optimization\n",
    "            model_output_data_point = np.sum(lambda_values*i)\n",
    "            self.model_output.append(model_output_data_point)\n",
    "            #print(i)\n",
    "            if model_output_data_point <= thresholds[1]:\n",
    "                self.predicted_class.append(0.)\n",
    "            elif model_output_data_point > thresholds[-1]:\n",
    "                self.predicted_class.append(float(self.n_classes - 1))\n",
    "            else:\n",
    "                for j in range(1, self.n_classes - 1):\n",
    "                    if thresholds[j] < model_output_data_point <= thresholds[j + 1]:\n",
    "                        self.predicted_class.append(j)\n",
    "                        \n",
    "    \n",
    "    def class_prediction_with_error(self, lambda_values):\n",
    "        # Given measurement parameters lambda,\n",
    "        # For the results of an Ascella experiment, obtain model prediction\n",
    "        # With an error estimated from counts\n",
    "        thresholds = np.linspace(-1, 1, self.n_classes, False)\n",
    "        \n",
    "        for prob_index, proba in enumerate (self.probabilities):\n",
    "            i = proba*100000 # factor added for the optimization\n",
    "            distribution_counts = self.outcomes[prob_index]\n",
    "            model_output_data_point = np.sum(lambda_values*i)\n",
    "            model_output_error_data_point = [(lambda_values[j]**2)*(poisson.std(distribution_counts[j])**2)\n",
    "                                             for j in range(len(lambda_values))]\n",
    "            \n",
    "            self.model_output.append(model_output_data_point)\n",
    "            self.model_output_error.append(np.sqrt(np.sum(model_output_error_data_point)))\n",
    "            \n",
    "            if model_output_data_point <= thresholds[1]:\n",
    "                self.predicted_class.append(0.)\n",
    "            elif model_output_data_point > thresholds[-1]:\n",
    "                self.predicted_class.append(float(self.n_classes - 1))\n",
    "            else:\n",
    "                for j in range(1, self.n_classes - 1):\n",
    "                    if thresholds[j] < model_output_data_point <= thresholds[j + 1]:\n",
    "                        self.predicted_class.append(j)\n",
    "\n",
    "                        \n",
    "    def cost_estimator(self, lambda_values):\n",
    "        # Given measurement parameters lambda and results from Ascella experiment, \n",
    "        # Compute all predictions and compute accuracy and cost\n",
    "        cost = 0\n",
    "        self.class_prediction(lambda_values)\n",
    "        y_predicted = np.array(self.predicted_class)\n",
    "        accuracy = accuracy_score(self.y_train, y_predicted)\n",
    "        cost = 1 - accuracy\n",
    "        self.predicted_class = []\n",
    "        return cost\n",
    "    \n",
    "    \n",
    "    def reset_outcomes(self):\n",
    "        # Reset probabilities, outcomes, predictions before running another experiment\n",
    "        self.probabilities = []\n",
    "        self.outcomes = []\n",
    "        self.predicted_class = []\n",
    "        self.model_output = []\n",
    "        self.model_output_error = []\n",
    "        self.predicted_class = []\n",
    "        \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34ef1f94",
   "metadata": {},
   "source": [
    "## Train the model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9ac95e0",
   "metadata": {},
   "source": [
    "We can now define the classifier using the iris dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51ba5504",
   "metadata": {},
   "outputs": [],
   "source": [
    "Classifier = QuantumClassifierOnAscella(dataset = datasets.load_iris(), n_classes = 3, run_on_qpu = run_on_qpu)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08e4af1d",
   "metadata": {},
   "source": [
    "Print the Ascella circuit below.\n",
    "\n",
    "Notice that several phases are already initialised in the definition of the QuantumClassifierOnAscella class.\n",
    "\n",
    "In order to train the model faster, we do not use transpilation in this notebook. This means we hardcode the phases that we use to encode the data and the chip parameters. This is done in the definition of the QuantumClassifierOnAscella class (see variables theta_phases and data_phases)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc03df80",
   "metadata": {},
   "outputs": [],
   "source": [
    "pcvl.pdisplay(Classifier.circuit)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e03b6ee2",
   "metadata": {},
   "source": [
    "Now we can launch the optimisation.\n",
    "\n",
    "As explained in the article, we optimise separately the chip parameters ($\\theta$) and the measurements parameters ($\\lambda$). \n",
    "\n",
    "We use different optimisers for each set. The theta optimisation is done via Gaussian processes, while the lambda optimisation is done via Nelder-Mead."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85891ef8",
   "metadata": {},
   "outputs": [],
   "source": [
    "theta_optimisation_space = [sk.space.space.Real(0, 2*np.pi) for _ in range(len(Classifier.theta_phases))]\n",
    "\n",
    "bounds=[]\n",
    "for i in range(len(Classifier.possible_output_states)):\n",
    "    bounds.append((-4, 4))\n",
    "\n",
    "theta_optimizer = sk.Optimizer(theta_optimisation_space, \"GP\", acq_func=\"EI\", acq_optimizer=\"sampling\",\n",
    "                               initial_point_generator=\"lhs\")\n",
    "\n",
    "best_accuracy = 0\n",
    "best_params = []\n",
    "best_iteration = 0\n",
    "\n",
    "# Split the data - each batch will be sent as one job to the QPU \n",
    "# The split is needed in order to stay below the runtime limit\n",
    "n_splits = 10\n",
    "split_data = np.array_split(Classifier.X_train, n_splits)\n",
    "\n",
    "for i in range(15):\n",
    "    print(\"Iteration \" + str(i))\n",
    "    next_theta = theta_optimizer.ask()\n",
    "    \n",
    "    start_watch = time.time()\n",
    "    \n",
    "    if run_on_qpu:\n",
    "        Classifier.set_unitaries_ascella(next_theta)\n",
    "        # send data to qpu with split\n",
    "        for data_batch in split_data:\n",
    "            list_data_phases = []\n",
    "            list_data_phases = Classifier.prepare_batch_data_phases(data_batch)\n",
    "            Classifier.get_results_ascella_batch(list_data_phases)\n",
    "    else:\n",
    "        Classifier.get_results_ascella_all_data(next_theta) \n",
    "    \n",
    "    end_watch = time.time()\n",
    "    print('Obtained results from experiment in ' + str(round(end_watch - start_watch)) + ' seconds.')\n",
    "    \n",
    "    start_watch = time.time()\n",
    "    opt_lambda = scipy.optimize.minimize(Classifier.cost_estimator,\n",
    "                                         x0 = [0.0]*len(Classifier.possible_output_states),\n",
    "                                         method = \"Nelder-Mead\", bounds=bounds)\n",
    "    end_watch = time.time()\n",
    "    print('Lambda optimization done in ' + str(round(end_watch - start_watch)) + ' seconds.')\n",
    "\n",
    "    print('Objective function:' + str(opt_lambda.fun))\n",
    "    \n",
    "    # keep best accuracy\n",
    "    if (1 - opt_lambda.fun) > best_accuracy:\n",
    "        best_accuracy = (1 - opt_lambda.fun)\n",
    "        best_params = [*next_theta, *opt_lambda[\"x\"]]\n",
    "        best_iteration = i\n",
    "    \n",
    "    # write to output files\n",
    "    outcome_report = {}\n",
    "    for j in range(len(Classifier.outcomes)):\n",
    "        outcome_report[\"Counts\" + str(j)] = Classifier.outcomes[j]\n",
    "        outcome_report[\"Prob\" + str(j)] = Classifier.probabilities[j]\n",
    "\n",
    "    parameter_report = {\"Parameters\": [*next_theta, *opt_lambda[\"x\"]]}\n",
    "    iteration_report = {\"Time\": [datetime.now().strftime(\"%Y%m%d_%H%M%S\")],\n",
    "                        \"Accuracy\": [1 - opt_lambda.fun], \"Cost\": [opt_lambda.fun]}\n",
    "    df_iteration = pd.DataFrame(data = iteration_report)\n",
    "    df_iteration.to_csv(Classifier.path + \"Iteration\" + str(i) + \".csv\", index = False)\n",
    "    df_outcome = pd.DataFrame(data = outcome_report)\n",
    "    df_outcome.to_csv(Classifier.path + \"Outcomes\" + str(i) + \".csv\", index = False)\n",
    "    df_parameters = pd.DataFrame(data = parameter_report)\n",
    "    df_parameters.to_csv(Classifier.path + \"Parameters\" + str(i) + \".csv\", index = False)\n",
    "    \n",
    "    # feed round to sk.Optimizer and reset\n",
    "    res = theta_optimizer.tell(next_theta, opt_lambda.fun)\n",
    "    Classifier.reset_outcomes()\n",
    "    \n",
    "print('Best iteration is ' + str(best_iteration) + ' with accuracy ' + str(best_accuracy) )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "937b88ca",
   "metadata": {},
   "source": [
    "## Evaluate the model on test set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "719ea52a",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_splits_tests = 2\n",
    "split_data_test = np.array_split(Classifier.X_test, n_splits_tests)\n",
    "\n",
    "# Set theta parameters to the best ones and run\n",
    "if run_on_qpu:\n",
    "    Classifier.set_unitaries_ascella(best_params[:len(Classifier.theta_phases)])\n",
    "    # send data to qpu with split\n",
    "    for data_batch_test in split_data_test:\n",
    "        list_data_phases = []\n",
    "        list_data_phases = Classifier.prepare_batch_data_phases(data_batch_test)\n",
    "        Classifier.get_results_ascella_batch(list_data_phases)\n",
    "else:\n",
    "    Classifier.set_unitaries_ascella(best_params[:len(Classifier.theta_phases)])\n",
    "    for x in Classifier.X_test:\n",
    "        Classifier.get_results_ascella(x)\n",
    "\n",
    "Classifier.class_prediction(best_params[len(Classifier.theta_phases):])\n",
    "y_predicted_test = np.array(Classifier.predicted_class)\n",
    "accuracy_test = accuracy_score(Classifier.y_test, y_predicted_test)\n",
    "print(\"Test set accuracy is: \" + str(accuracy_test))\n",
    "\n",
    "# Write to output files\n",
    "outcome_report = {}\n",
    "for j in range(len(Classifier.outcomes)):\n",
    "    outcome_report[\"Counts\" + str(j)] = Classifier.outcomes[j]\n",
    "    outcome_report[\"Prob\" + str(j)] = Classifier.probabilities[j]\n",
    "\n",
    "parameter_report = {\"Parameters\": best_params}\n",
    "iteration_report = {\"Time\": [datetime.now().strftime(\"%Y%m%d_%H%M%S\")],\n",
    "                    \"Accuracy\": [accuracy_test], \"Cost\": [1 - accuracy_test]}\n",
    "df_iteration = pd.DataFrame(data = iteration_report)\n",
    "df_iteration.to_csv(Classifier.path + \"Iteration\" + \"_test\" + \".csv\", index = False)\n",
    "df_outcome = pd.DataFrame(data = outcome_report)\n",
    "df_outcome.to_csv(Classifier.path + \"Outcomes\" + \"_test\" + \".csv\", index = False)\n",
    "df_parameters = pd.DataFrame(data = parameter_report)\n",
    "df_parameters.to_csv(Classifier.path + \"Parameters\" + \"_test\" + \".csv\", index = False)\n",
    "\n",
    "# Reset\n",
    "Classifier.reset_outcomes()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa19eb42",
   "metadata": {},
   "source": [
    "## Reproduce plots from article and get confusion matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82e779d6",
   "metadata": {},
   "source": [
    "This can be executed independently of the optimisation as it recovers results saved in csv file. \n",
    "\n",
    "Change the path as needed."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74e3219a",
   "metadata": {},
   "source": [
    "**Train set**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7624187",
   "metadata": {},
   "outputs": [],
   "source": [
    "best_params_train_path = Classifier.path + \"Parameters\" + str(best_iteration) + \".csv\"\n",
    "best_outcomes_train_path = Classifier.path + \"Outcomes\" + str(best_iteration) + \".csv\"\n",
    "\n",
    "# If needed: change path and redefine Classifier class\n",
    "# Classifier = QuantumClassifierOnAscella(dataset = datasets.load_iris(), n_classes = 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82399530",
   "metadata": {},
   "outputs": [],
   "source": [
    "best_params_csv = pd.read_csv(best_params_train_path).squeeze(\"columns\")\n",
    "best_params = [i for i in best_params_csv.values.tolist()]\n",
    "\n",
    "outcomes_df = pd.read_csv(best_outcomes_train_path)\n",
    "outcomes_array = np.array(outcomes_df).transpose()\n",
    "\n",
    "outcomes_from_csv = []\n",
    "probs_from_csv = []\n",
    "for i in range(0, outcomes_array.shape[0], 2):\n",
    "    outcomes_from_csv.append(outcomes_array[i])\n",
    "    probs_from_csv.append(outcomes_array[i + 1])\n",
    "    \n",
    "Classifier.outcomes = outcomes_from_csv\n",
    "Classifier.probabilities = probs_from_csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4c558d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "Classifier.class_prediction_with_error(best_params[len(Classifier.theta_phases):])\n",
    "y_predicted_train = np.array(Classifier.predicted_class)\n",
    "accuracy_train = accuracy_score(Classifier.y_train, y_predicted_train)\n",
    "cm = confusion_matrix(Classifier.y_train, y_predicted_train)\n",
    "disp = ConfusionMatrixDisplay(confusion_matrix=cm)\n",
    "disp.plot()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58a9fc6f",
   "metadata": {},
   "outputs": [],
   "source": [
    "f = plt.figure()\n",
    "f.set_figwidth(10)\n",
    "f.set_figheight(8)\n",
    "\n",
    "# Color codes\n",
    "light_purple_code = '#efdeff'\n",
    "dark_purple_code = '#af5cff'\n",
    "light_blue_code = '#cce2ff'\n",
    "dark_blue_code = '#516d91'\n",
    "\n",
    "point_color_dict = {0:dark_purple_code, 1:dark_blue_code, 2:'black'}\n",
    "thresholds = np.linspace(-1, 1, Classifier.n_classes, False)\n",
    "\n",
    "for i in range(len(Classifier.model_output)):\n",
    "    plt.errorbar(x = i, \n",
    "                y = Classifier.model_output[i],\n",
    "                yerr = (Classifier.model_output_error[i]), \n",
    "                c = point_color_dict[Classifier.y_train.tolist()[i]],\n",
    "                fmt = '.')\n",
    "\n",
    "plt.hlines(thresholds, xmin = 0, xmax = len(Classifier.model_output), color = 'grey', linestyle = 'dashed')\n",
    "plt.fill_between([0, 112], [thresholds[0], thresholds[0]], [thresholds[1], thresholds[1]], facecolor = light_purple_code)\n",
    "plt.fill_between([0, 112], [thresholds[1], thresholds[1]], [thresholds[2], thresholds[2]], facecolor = light_blue_code)\n",
    "plt.fill_between([0, 112], [thresholds[2], thresholds[2]], [1, 1], facecolor = 'lightgrey')\n",
    "\n",
    "plt.xlim([0,len(Classifier.model_output)])\n",
    "plt.ylim([-1,1])\n",
    "\n",
    "plt.ylabel(r'Class Estimator', fontsize = 18)\n",
    "plt.xlabel('')\n",
    "plt.xticks([])\n",
    "\n",
    "plt.rcParams.update({\n",
    "    'savefig.dpi': 300,  # to adjust notebook inline plot size\n",
    "})\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4c8c791",
   "metadata": {},
   "outputs": [],
   "source": [
    "Classifier.reset_outcomes()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7f9bba6",
   "metadata": {},
   "source": [
    "**Test set**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69b3d910",
   "metadata": {},
   "outputs": [],
   "source": [
    "best_params_test_path = Classifier.path + \"Parameters\" + \"_test\" + \".csv\"\n",
    "best_outcomes_test_path = Classifier.path + \"Outcomes\" + \"_test\" + \".csv\"\n",
    "\n",
    "# If needed: change path and redefine Classifier class\n",
    "# Classifier = QuantumClassifierOnAscella(dataset = datasets.load_iris(), n_classes = 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e117e3d",
   "metadata": {},
   "outputs": [],
   "source": [
    "best_params_csv = pd.read_csv(best_params_test_path).squeeze(\"columns\")\n",
    "best_params = [i for i in best_params_csv.values.tolist()]\n",
    "\n",
    "outcomes_df = pd.read_csv(best_outcomes_test_path)\n",
    "outcomes_array = np.array(outcomes_df).transpose()\n",
    "\n",
    "outcomes_from_csv = []\n",
    "probs_from_csv = []\n",
    "for i in range(0, outcomes_array.shape[0], 2):\n",
    "    outcomes_from_csv.append(outcomes_array[i])\n",
    "    probs_from_csv.append(outcomes_array[i + 1])\n",
    "    \n",
    "Classifier.outcomes = outcomes_from_csv\n",
    "Classifier.probabilities = probs_from_csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3fe2f30",
   "metadata": {},
   "outputs": [],
   "source": [
    "Classifier.class_prediction_with_error(best_params[len(Classifier.theta_phases):])\n",
    "y_predicted_test = np.array(Classifier.predicted_class)\n",
    "accuracy_test = accuracy_score(Classifier.y_test, y_predicted_test)\n",
    "cm = confusion_matrix(Classifier.y_test, y_predicted_test)\n",
    "disp = ConfusionMatrixDisplay(confusion_matrix=cm)\n",
    "disp.plot()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79cc0d3d",
   "metadata": {},
   "outputs": [],
   "source": [
    "f = plt.figure()\n",
    "f.set_figwidth(10)\n",
    "f.set_figheight(8)\n",
    "\n",
    "# Color codes\n",
    "light_purple_code = '#efdeff'\n",
    "dark_purple_code = '#af5cff'\n",
    "light_blue_code = '#cce2ff'\n",
    "dark_blue_code = '#516d91'\n",
    "\n",
    "point_color_dict = {0:dark_purple_code, 1:dark_blue_code, 2:'black'}\n",
    "thresholds = np.linspace(-1, 1, Classifier.n_classes, False)\n",
    "\n",
    "for i in range(len(Classifier.model_output)):\n",
    "    plt.errorbar(x = i, \n",
    "                y = Classifier.model_output[i],\n",
    "                yerr = (Classifier.model_output_error[i]), \n",
    "                c = point_color_dict[Classifier.y_test.tolist()[i]],\n",
    "                fmt = '.')\n",
    "\n",
    "plt.hlines(thresholds, xmin = 0, xmax = len(Classifier.model_output), color = 'grey', linestyle = 'dashed')\n",
    "plt.fill_between([0, 112], [thresholds[0], thresholds[0]], [thresholds[1], thresholds[1]], facecolor=light_purple_code)\n",
    "plt.fill_between([0, 112], [thresholds[1], thresholds[1]], [thresholds[2], thresholds[2]], facecolor=light_blue_code)\n",
    "plt.fill_between([0, 112], [thresholds[2], thresholds[2]], [1, 1], facecolor='lightgrey')\n",
    "\n",
    "plt.xlim([0,len(Classifier.model_output)])\n",
    "plt.ylim([-1,1])\n",
    "\n",
    "plt.ylabel(r'Class Estimator', fontsize = 18)\n",
    "plt.xlabel('')\n",
    "plt.xticks([])\n",
    "\n",
    "plt.rcParams.update({\n",
    "    'savefig.dpi': 300,  # to adjust notebook inline plot size\n",
    "})\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2dfbe531",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e39dfd84",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}