Dynamic_Traffic_Management_with_Qutrits.ipynb

{
 "cells": [
  {
   "source": [
    "# Dynamic Traffic Management with Qutrits\n",
    "\n",
    "The implementation that we conceptualized for the use of three distinct energy states in a Quantum Computer was the use of variational hybrid algorithms to optimize traffic grids. While it was at first quite challenging to think of a real-world use for three orthogonal basis states (the |0>, |1>, and |2> states), we eventually came to the conclusion that four-way intersections behave with three distinct settings that can be generalized as nodes in quantum systems, affecting the energy of the system like the settings of traffic lights affect the congestion in the grid."
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "source": [
    "### Importing the necessary packages"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "centered-bachelor",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "configrc.store_credentials:WARNING:2021-04-11 08:31:01,089: Credentials already present. Set overwrite=True to overwrite.\n",
      "ibmqfactory.load_account:WARNING:2021-04-11 08:31:01,286: Credentials are already in use. The existing account in the session will be replaced.\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import networkx as nx\n",
    "import matplotlib.pyplot as plt \n",
    "import random\n",
    "from matplotlib import cm\n",
    "from matplotlib.ticker import LinearLocator, FormatStrFormatter\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "from qiskit import Aer, IBMQ\n",
    "from qiskit import QuantumRegister, ClassicalRegister, QuantumCircuit, transpile, assemble\n",
    "from qiskit.providers.ibmq import least_busy\n",
    "from qiskit.tools.monitor import job_monitor\n",
    "from qiskit.visualization import plot_histogram\n",
    "from scipy.optimize import curve_fit\n",
    "import qiskit.pulse as pulse\n",
    "import qiskit.pulse.library as pulse_lib\n",
    "from qiskit.compiler import assemble\n",
    "from qiskit.pulse.library import SamplePulse\n",
    "\n",
    "from qiskit import IBMQ \n",
    "IBMQ.load_account() \n",
    "provider = IBMQ.get_provider(hub=\"ibm-q\", group=\"open\", project=\"main\") \n",
    "backend = provider.get_backend(\"ibmq_armonk\") \n",
    "backend_config = backend.configuration()\n"
   ]
  },
  {
   "source": [
    "## Data Structure to model traffic intersections\n",
    "\n",
    "We implemented a custom data structure in Python: `TrafficGrid` that would allow the user to create an n-by-m grid of nodes (intersections), and to declare which segments of road between the nodes are being used as roads. To better map this situation to a traditional graph-based combinatorial optimization problem, the user can also add “weights” to each edge in the graph, which we conceptualized as some function of the cars held up on that street times the amount of time they have been waiting."
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "radio-chocolate",
   "metadata": {},
   "outputs": [],
   "source": [
    "class TrafficGrid:\n",
    "    def __init__(self, x_length, y_length, p_left, p_right):\n",
    "        \"\"\"Data structure for a 2-dimensional traffic grid that can be\n",
    "        represented as a (x_length * y_length) grid of nodes (intersections),\n",
    "        (x_length - 1) * y_length North/South roads, and\n",
    "        (y_length - 1) * x_length East/West roads.\n",
    "\n",
    "        p_left is the probability of a car turning left at any 3-pronged intersection,\n",
    "        and p_right is the probability of a car turning right at any 3-pronged\n",
    "        intersection. These are used probabilistically to optimize how much\n",
    "        pressure the next time period will result in for the system. In large-scale\n",
    "        implementation, these could be implemented as an array with specific\n",
    "        approximations for each intersection in the grid for more accuracy.\n",
    "        Probability of proceeding straight is assumed 1 - p_left - p_right.\n",
    "        \"\"\"\n",
    "\n",
    "        self.x_length = x_length\n",
    "        self.y_length = y_length\n",
    "\n",
    "        self.p_left = p_left\n",
    "        self.p_right = p_right\n",
    "        self.p_straight = 1 - p_left - p_right\n",
    "\n",
    "        # Lists of street values from grid, vehicle weights & pedestrian weights\n",
    "        self.ns_streets = [[] for i in range(x_length)]\n",
    "        self.n_weights = [[] for i in range(x_length)]\n",
    "        self.s_weights = [[] for i in range(x_length)]\n",
    "\n",
    "        self.ew_streets = [[] for i in range(y_length)]\n",
    "        self.e_weights = [[] for i in range(y_length)]\n",
    "        self.w_weights = [[] for i in range(y_length)]\n",
    "\n",
    "    def add_ns_streets(self, streets):\n",
    "        \"\"\"Method to populate the traffic grid with North/South streets.\n",
    "        Streets is a list of lists, where each list is a column, and the\n",
    "        contents include which segments of the column are valid streets.\n",
    "\n",
    "        E.g. streets = [[0, 1], [0, 1, 2]] indicates two blocks in the first\n",
    "        column, and three blocks from the second column.\n",
    "        \"\"\"\n",
    "        self.ns_streets = streets\n",
    "        self.n_weights = [[0 for i in j] for j in streets]\n",
    "        self.s_weights = [[0 for i in j] for j in streets]\n",
    "\n",
    "    def add_ew_streets(self, streets):\n",
    "        \"\"\"Method to populate the traffic grid with East/West streets.\n",
    "        Streets is a list of lists, where each list is a column, and the\n",
    "        contents include which segments of the column are valid streets.\n",
    "        \"\"\"\n",
    "        self.ew_streets = streets\n",
    "        self.e_weights = [[0 for i in j] for j in streets]\n",
    "        self.w_weights = [[0 for i in j] for j in streets]\n",
    "\n",
    "    def add_n_weight(self, column, street, weight):\n",
    "        self.n_weights[column][street] += weight\n",
    "\n",
    "    def add_s_weight(self, column, street, weight):\n",
    "        self.s_weights[column][street + 1] += weight\n",
    "\n",
    "    def add_e_weight(self, row, street, weight):\n",
    "        self.e_weights[row][street] += weight\n",
    "\n",
    "    def add_w_weight(self, row, street, weight):\n",
    "        self.w_weights[row][street] += weight\n",
    "\n",
    "    def add_n_pweight(self, column, street, weight):\n",
    "        self.n_pweights[column][street] += weight\n",
    "\n",
    "    def add_s_pweight(self, column, street, weight):\n",
    "        self.s_pweights[column][street + 1] += weight\n",
    "\n",
    "    def add_e_pweight(self, row, street, weight):\n",
    "        self.e_pweights[row][street] += weight\n",
    "\n",
    "    def add_w_pweight(self, row, street, weight):\n",
    "        self.w_pweights[row][street] += weight\n",
    "\n",
    "    def intersections(self):\n",
    "        \"\"\"Produces list of intersections\"\"\"\n",
    "        intersections = set()\n",
    "        for i, column in enumerate(self.ns_streets):\n",
    "            for street in column:\n",
    "                n_node = (self.x_length * street) + i\n",
    "                s_node = (self.x_length * (street + 1)) + i\n",
    "                intersections.add(n_node)\n",
    "                intersections.add(s_node)\n",
    "\n",
    "        for i, row in enumerate(self.ew_streets):\n",
    "            for j, street in enumerate(row):\n",
    "                e_node = (self.x_length * i) + street\n",
    "                w_node = (self.x_length * i) + street + 1\n",
    "                intersections.add(e_node)\n",
    "                intersections.add(w_node)\n",
    "\n",
    "        return intersections\n",
    "\n",
    "    def edges(self):\n",
    "        \"\"\"Produces list of edges/weights\"\"\"\n",
    "        edges = []\n",
    "        for i, column in enumerate(self.ns_streets):\n",
    "            for j, street in enumerate(column):\n",
    "                n_weight = self.n_weights[i][j]\n",
    "                s_weight = self.s_weights[i][j]\n",
    "\n",
    "                n_node = (self.x_length * street) + i\n",
    "                s_node = (self.x_length * (street + 1)) + i\n",
    "                edges.append((n_node, s_node, s_weight))\n",
    "                edges.append((s_node, n_node, n_weight))\n",
    "\n",
    "        for i, row in enumerate(self.ew_streets):\n",
    "            for j, street in enumerate(row):\n",
    "                e_weight = self.e_weights[i][j]\n",
    "                w_weight = self.w_weights[i][j]\n",
    "\n",
    "                e_node = (self.x_length * i) + street\n",
    "                w_node = (self.x_length * i) + street + 1\n",
    "                edges.append((e_node, w_node, w_weight))\n",
    "                edges.append((w_node, e_node, e_weight))\n",
    "\n",
    "        return edges\n",
    "\n",
    "    def graph(self):\n",
    "        \"\"\"Creates & plots a bidirectional graph object of intersections\n",
    "        (as nodes) connected by two-directional streets (as edges).\n",
    "        \"\"\"\n",
    "        G = nx.DiGraph()\n",
    "        G.add_nodes_from(self.intersections())\n",
    "        G.add_weighted_edges_from(self.edges())\n",
    "\n",
    "        # Generate plot of the Graph\n",
    "        colors  = [\"r\" for node in G.nodes()]\n",
    "        default_axes = plt.axes(frameon=True)\n",
    "        pos = nx.spring_layout(G)\n",
    "\n",
    "        nx.draw_networkx(G, node_color=colors, node_size=600, alpha=1, \n",
    "                        ax=default_axes, pos=pos)\n",
    "        return G\n",
    "\n",
    "    def meshgrid(self, step_size=0.1):\n",
    "        \"\"\"Plots a 3D map to indicate ground state graphically.\n",
    "        Source: IBM Qiskit QAOA Documentation.\n",
    "        \"\"\"\n",
    "        # Evaluate the function\n",
    "        a_gamma = np.arange(0, np.pi, step_size)\n",
    "        a_beta = np.arange(0, np.pi, step_size)\n",
    "        a_gamma, a_beta = np.meshgrid(a_gamma, a_beta)\n",
    "\n",
    "        F1_1 = np.sin(2 * a_beta)**2 * np.sin(2 * a_gamma)**2 \n",
    "        F1_1 -= 0.5 * np.sin(4 * a_beta) * np.sin(4 * a_gamma)\n",
    "        F1_2 = 1 + np.cos(4 * a_gamma)**2\n",
    "        F1 = 3 - (F1_1 * F1_2)\n",
    "\n",
    "        # Grid search for the minimizing variables\n",
    "        result = np.where(F1 == np.amax(F1))\n",
    "        a = list(zip(result[0], result[1]))[0]\n",
    "\n",
    "        gamma = a[0] * step_size\n",
    "        beta = a[1] * step_size\n",
    "\n",
    "        # Plot the expectation value F1\n",
    "        fig = plt.figure()\n",
    "        ax = Axes3D(fig)\n",
    "\n",
    "        surf = ax.plot_surface(a_gamma, a_beta, F1, cmap=cm.coolwarm, \n",
    "                                linewidth=0, antialiased=True)\n",
    "\n",
    "        ax.set_zlim(1,4)\n",
    "        ax.zaxis.set_major_locator(LinearLocator(3))\n",
    "        ax.zaxis.set_major_formatter(FormatStrFormatter(\"%.02f\"))\n",
    "        plt.show()\n",
    "\n",
    "        # The smallest parameters and the expectation can be extracted\n",
    "        print(\"\\n --- OPTIMAL PARAMETERS --- \\n\")\n",
    "        print(f\"The maximal expectation value is: M1 = {np.amax(F1):.03f}\")\n",
    "        print(f\"This is attained for gamma = {gamma:.03f} and {beta:.03f}\")\n"
   ]
  },
  {
   "source": [
    "### Initializing `TrafficGrid`"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "ultimate-posting",
   "metadata": {},
   "outputs": [],
   "source": [
    "grid = TrafficGrid(3, 4, 0.25, 0.25)\n",
    "grid.add_ns_streets([[0, 1], [0, 1, 2], [0, 1, 2]])\n",
    "grid.add_ew_streets([[0, 1], [0], [0, 1], [0, 1]])\n",
    "grid.add_n_weight(1, 2, 1)\n",
    "grid.add_e_weight(3, 0, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "victorian-disney",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<__main__.TrafficGrid at 0x7f51cc80a950>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "signal-queue",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<networkx.classes.digraph.DiGraph at 0x7f51cc410890>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADnCAYAAAC9roUQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8+yak3AAAACXBIWXMAAAsTAAALEwEAmpwYAABe6klEQVR4nO2deXhMZ/vHv2f2LYkgkYhEhEgIEmKp2KldLcXPXvtaFKUvVUVfSynVl1arVLVVrddWXVS1qkrVUkUpqmrfUiUiIstk8v39cWReIcskc2bOJDmf6zqXmDnP89znZPKd5zz3/dy3QBIKCgoKCu5BJbcBCgoKCiUJRXQVFBQU3IgiugoKCgpuRBFdBQUFBTeiiK6CgoKCG9Hk9WbZsmUZGhrqJlMUFBQUigeHDx/+h6RfTu/lKbqhoaH45ZdfXGOVgoKCQjFFEISLub2nLC8oKCgouBFFdBUUFBTciCK6CgoKCm5EEV0FBQUFN5KnI03BBZDA+fPA4cPAmTPA/fuA0QiEhwOxsUDlyoAgyG2lgoKCi1BE113cvg289x7w+uvA3buARgPcuwdkZgIqFWA2Azab+O+ECcDw4YBfjhEnCgoKRRhlecHVkMCHHwIVKwIzZwI3boiz27t3RcEFxH+TksTXb94E5swBKlUCVqwQ2ysoKBQbFNF1JSkpQMeOwJgx4qw2JcXxdsnJwPPPAy1bim0VFBSKBYrouorUVODJJ4EffhAFtDAkJwP79wNNm4qzYAUFhSKPIrquYsQI4MgRx2e3uZGaCpw6BQwYII1dCgoKsqI40lzBjh3Apk3OC24WqanA9u3AZ58BXbtK02dRhQSOHQMOHgR++gm4fl18LSAAiIsD6tUTo0CUCBAFD0XIq1xP3bp1qeReKCCZmUBwMHDtmvR9lykjOuI0JfC7MjUVWL0aWLgQ+OcfUWgfXXIxGgG1GvDxASZPFiNAzGZ57FUo0QiCcJhk3ZzeU5YXpObbb8VIBFeQng58/rlr+vZkDhwAqlYFpkwBLl4U17pzWuNOSRGdjlevAtOni7HPe/a4314FhTxQRFdq/vOfAonupwCqATADqAwgT4lISgKWLHHKvCLH0qVAixbA5csFcybevy8uPbRrB8yb5zr7FBQKSAl8TnUxP//s8KnfAvgXgPUA6gO47kijw4fFR+uSsGa5dCkwbZpza+P37wNz54rLPi+9JJ1tCgqFRJnpSsmNGwUSiJkAXgbwBMRfRNCDI08EATh3rrAWFh327QOmTpUmVO7+fWD+fGDnTuf7UlBwEkV0peTKFcBgcOhUG4BfANwEUAVABQBjAeQr2VqtOE5xJiUF+L//ky76AxCFt08f1623Kyg4iCK6UpKR4fCp8QCsADZCXMc9CuAIgDkSj1MkWbECSEiQvt+SuCau4HEoa7pS4u39v3wK+WB88O84AIEPfp4EUXTn5tWQFMcprpDAokUOLSu8CWANgOMA+jz4OU9SU8V14unTxdCykozNBvzyi+gjOHAAuHVLvCdVqgD16wMNGwIhIXJbWSxRRFdKwsPFP2wH8IW4pPCwO8wh19j9+0BUVMFtKyocPAgkJjp0ankALwH4Bg4sy2SRni5uzW7VqlDmFXlu3waWLxejbNLSxKemh5dxBAGwWACrFahdW1xXf+qpkuG4dRPK8oKUaLVAWJjDpw8GsAzA3wASACwB0Cm/RuXLAyZTYS30fA4edHj55GkAXQGUKUj/qanizK4ksmWL+PmcO1fcYJKU9Pi6OSm+npoqRuL07Svm/rh8WR6biyGK6ErNM8+IO6McYAaAegCqQozVrQ1gel4N9Hrxj6A4s2ePw08LhcJqBX780XX9eyKZmcCwYUD//uJTREHub1bSperVge+/d52NJQhFdKVm+HCHc+BqASwHcAfADQBLAeQZ+yAIwLPPOmmgh+OK7dOPct2hiOjiAQkMHAh88knhw+8yMsSdfk89BezaJa19JRBFdKXGz0+c7ToYOuYwej3w9NNAhQrS9utpuCNpe0lKDP/mm8DmzdLFO3fpIsajKxQaRXRdweuvA6VKSdunxQK89Za0fcrMqVOncOHCBWRLuuSOEkVly7p+DE/g/HnpNphkkZoqTipK0heXxCii6wrMZjENo1QOL6MR2LhReiGXmUGDBqFKlSooXbo02rVrh4EDB2JbUhJsWq1D7TMApELcaGJ78HO+Lji1GmjSxBmziw6TJ4sRClJitYq7BUvauriEKKLrKho0EDOCmUzOhduYTMCGDUDz5pKZ5gncv38fpUqVgs1mw507d/DNN9/gww8/xOZLl8QoEAeYAzHe+VUAax/8nO/mErNZjEMt7vz9N7BtmxiPmw+nALQE4ANxd+SW/Brcvw+89prTJpZUFNF1Ja1aiZ7f8PCCz3pNJrE45Y8/inXWiji3b9/Gf/7zH7Rq1QplypSB2WzGrl27oFKJH0GDwYBNmzZh1alTUDsY/TELAB85ZuXXiBTLKBV31q936Ms+A0AXiKGKtwG8C6A/gDN5NSLFFKZ370pgaMlDEV1XU7MmcOKEmOHK1xfw8srz9DSdDvTxAV54QSzTExvrJkOl5dq1a5g7dy4aNWoEHx8flClTBtOnT8fdu3cxZswYnD17Fv/88w8EQUDp0qXx008/oVu3bkjPzMS5Tp2Q4eBst0DodMDIkaJTsrizc6dDuStOA7gGYCIANcQZbyMAH+XX0GgUy1EpFBySuR6xsbFUkBCrldy8mRw8mIyIoE2vp1UQSIOBDA8nn3mGY8uV44DeveW2tMCcOXOG06ZNY2xsLE0mEwHQx8eHjRs35vz583n9+vUc27Vt25Ymk4k6nY6CIBAAS6tUTLNYSHFOJd3h7U3Gx7v5zshEhQoO3ZPjAM0AMx967UmAXfNrq9eTixfLfZUeC4BfmIuuKqIrI8uXL6dGo6HNZrO/VqFCBQLgunXrZLQsf3799VeOHz+eNWrUoF6vJwCWKVOGTz75JJctW8aEhASH+vnmm2/s7QFQpVJx9+7d5JdfkiaTdIJrMpHr17v2pngSDt67dICVAC548PM3ALUA2zhyTydNkvsqPRZFdD2U8PBwAuCWLVvsr4WFhREADQYD33//fdlsexibzcbdu3dz2LBhDA8Pp1arpSAIDAgIYKdOnbhmzRqmpKQUuN/Lly+zUaNGdsE1Go185ZVX/nfC889LI7wmEzlihIR3pAhQgPt2DGBTgKUfiG0/gEMcaTtxotxX6bEoouuB/Pjjj9TpdATAiIgIZmZmkiQtFotdhNRqtSwzXpvNxi+++IJ9+/ZlxYoVqVarqVKpGBwczJ49e3Ljxo20Wq2F7j8pKYndu3enSqVi5cqVuWbNGqpUKtauXZsZGRn/OzEzU/zDdkZ4TSZy2DDyoaeJEoG/f6HvWUOA7+R3nkZDzp0r91V6LIroehiZmZmsXbu2XVzNZjO/+eYb3rt3jwCo0WgIgHFxcdyzZ4/L7UlLS+O6devYrVs3li9fnoIgUKPRMCwsjAMHDuSOHTuyLYEUFpvNxueee44ajYZly5blhg0b7O8tX76cFy5cyLnhpk1kqVKkTue4eOh0pJcX+fHHoniXNJ58skAz3RSAyQBfAxgKMDW/dj4+5Ndfy32VHosiuh7GuXPn7DPZLJF9+umnmZSUxDZt2nDhwoUEwF9++cUl4yclJXHFihVs164d/fz8CIA6nY6RkZEcPXo09+3bJ/mYS5Ysodlspslk4muvvVbwDm7eFJcbvLzEQ6V6XAgEgal6PTOMRnLsWDIX512JYO5ch7+kJgMsBdGh1g7gn46ItV5P/v233FfpsSii64HcvXuXs2fPpre3N69cucKkpKRs71epUoVdunSRZKybN29y8eLFbN68OX19fZm1fhodHc3JkyfzxIkTkoyTE5s3b6afnx81Gg3Hjh3r1LIESTI1ldy6lZw2jYyLI6tWFSM/nniCnDKFS1u0YOXy5SWZmRdpzp4Vo2IKuyyT39GggdxX6NEoouuhzJs3j6VLl871PZPJVKh+L126xNmzZ/OJJ56gl5cXAdBisbB+/fqcOXNm7o/xEnLw4EGGh4dTEAR27dqViYmJLh+TJEeOHEkAHDBggH2dvMTSpIlrBNdiEUMfFXIlL9FVNkfIiNVqte/IepSJEyciNTUV3333Xb79nD59Gv/6179Qp04dmEwmhISEYMmSJdDpdJgxYwbi4+ORlJSEAwcOYNasWahYsaLUl2Ln0qVLaNiwIRo0aICyZcviwoUL2LJlC7zdVGIo6UHhyU8//RTDhw8XZxYllUWLHM7t7DCCIJbxeeopafstQSiiKyN5ia7BYEBUVBTmzZv32Hu//PILxo0bh6ioKBgMBlSrVg2rV6+Gv78/Xn/9dSQmJiIhIQG7d+/GlClT4O/v7+pLwd27d9GtWzeEhobi1q1bOHDgAPbt24cQN9fZOnHiBADx3n7wwQf417/+5dbxPYr69YFRo6StNGIwiLlANEqlr8Ki3DkZSU9PhzqPAoljx47FuHHj8N1332HdunXYs2cPLly4AJvNhoCAANSvXx8vvvgiunfvDoPU+XsdJCMjAxMnTsQ777yD0qVLY9OmTejWrZsstgDAhQsXoNVqYbVaUatWLTRs2FA2WzyC+fPFrGDHjjlfkcNkEmurVa8ujW0lldzWHais6bqciRMnMigoKNtrVquVn332GXv37m3fnQaAISEh7NWrF7du3eq8M0oiFi1aRJPJRLPZzCVLlshtDkly9OjRXLlyJfV6PZcvXy63OZ5BUpLoaHQm3tloJN96S+4rKTJAcaR5JmPHjmVwcDA/+ugjdunShYGBgfYY2SpVqnDIkCGsXr0669WrJ7ep2diwYQPLli1LjUbDCRMmeGSkQLNmzVi7dm25zfAc0tOZNGkSU1UqUhAKJraBgeSPP8p9BUUKRXQ9iMTERC5fvpxt2rSx5xzQ6/WsXr06x44dy4MHD2Y7f+PGjVSr1R4xu/35559ZuXJlqlQqdu/e/bEwN09i8+bNVKvVHvmFIAcnT56k0Whk35o1ya5dxXAyszlnoVWpxFjo0qXJ2bPJ5GS5zS9yFD3RzcgQj2JAfHw8Fy5cyKZNm7JUqVIEQJPJxNq1a7NmzZoMCQnJtw+9Xs8333zTDdbmzPnz51m/fn0KgsBGjRrx8uXLstniKDabjVqtlqtXr5bbFFlJTU3liy++SK1WSwBcunSp+EZ8PLlqlZjxLiqKDAkhQ0PJRo3EGOht28SseAqFwvNF9/ffycmTybp1s687GY1k7drkhAnk0aPuscVJLly4wJkzZ7J+/fr2PApeXl5s2LAhX3nlFV66dMl+7oABA1i1atV8+2zdujWjoqJcaXaOJCQksFOnThQEgRERES7bIecq4uLiWL9+fbnNkI1Lly4xJCSERqPRvuvwa2XrrlvwXNHdv5+MjRXFVaPJfV1JrRbFuGZNcvdu19pUQE6cOMHJkyczOjra/uH29fVl8+bN+frrr/PmzZu5tu3Tpw+rVauW7xi7d++mIAhMdtNjntVq5ahRo6hWq1muXDlu3brVLeNKzccff/xY6sySRHx8POvXr2/fbu7t7f3Y8pWCa/A80U1NJcePF8W2MF7UESNkW2fav38/R48ezcjISHuWMD8/P7Zv357vvvtugdY5e/TowZo1azp0rsVi4ezZswtrtsPMnz+fRqORZrP5f4+iRRSbzUaNRsNPP/1UblNk4/jx4/YdiYIg8Ny5c3KbVCLwLNFNThb3zBdGcB9ddnDx1lKbzcYdO3Zw4MCBDAsLo0ajoSAILF++PLt168ZPPvmEaWlphe6/a9eujImJcejc7t27s1KlSoUeKz8+/fRTlilThlqtlpMmTSo2s8N69eqxcePGcpshG1WqVGFsbCyTkpK4cuXK7KkzFVyG54huRgbZrJk0iTj0erJePTI9XTLzrFYrN27cyJ49ezI4OJgqlYpqtZqhoaHs168fv/zyS0nFqGPHjnT0HmfNWPJarigMe/fuZaVKlahSqdizZ0+3LWG4i1WrVlGn08lthiy88847VKlUvHLlitymlDg8R3QXLZK+BMusWYU2JyUlhWvWrGGnTp0YEBBAQRCo1WpZtWpVDh8+nD/++KNLZ3xt27ZlgwJkaypdujQnTJggydhnz55lbGwsBUFg06ZNefXqVUn69TSsVivVanWRXZcuLGlpaTQajXz22WflNqVE4hmi+9dfzi0p5LXU8Pvv2Ybav38/e/To8ViWqYSEBC5btoxPPvkky5QpY4+RrVGjBp977jn++uuv0l2vA7Rq1YpxcXEOnz948GAGBgY6NWZCQgI7dOhAQRBYrVo1HjlyxKn+igIxMTFs0aKF3Ga4lZ49e9LHx6fYLBMVNTxDdEeMyDtCobCHSkU+VD133bp1NBqN1Ol03LdvH+fPn8/GjRvTx8eHWTGysbGxfPHFF3nmzBnprq8QNGvWjE2bNnX4/IsXLxJAoVIzpqWlcfjw4VSr1QwMDORXX31V4D6KKsuWLaPBYJDbDLdx4sQJCoKQrfaegnuRX3Tv3XNoWSEVYkG8EIAWgNEAtzkivAYDM//5h88++6y91E3W4e3tzbi4OM6dO9fjHqEbN27Mli1bFqhNYGAgBw0a5PD5NpuN//73v2kwGGixWEpkPoK0tDSqVCru2LFDblPcQnh4uMO+AgXXkJfouifL2K5dDqWCywAQDGA3gBAA2wD8H4DjAELzaGcVBIypVAmrHuRSzaJz587YunVrIY12PRkZGTAVMO1enz59sGbNGofO/fjjjzFu3Djcu3cPEydOxPz583NNJVmc0el0iIyMxOLFi9G6dWu5zXEp7777Lv766y9cvHhRblMUcsE9f4GHDgHJyfmeZgYwC6LAqgB0AlAJwOF82mlSUzGrfXvs3bsXb7zxBtq3bw9fX1+cP3/eObtdTEZGBjQFzEs6ffp0JCQk4NixY7me8+OPPyI0NBTPPPMM2rVrhzt37mDBggUlUnCzGDRoEPbs2SO3GS4lPT0dEyZMwKhRo1ChQgW5zVHIjdymwJRyeaFNm0Kt194AqAd4ypHzGzbMNmRmZibTJQwncwW1a9cuVB20sLAwduvW7bHXz5w5wzp16lAQBDZv3pzXS3JhxkdISUmhIAhuqa4sF4rzzHOA7OV6HnnsdwQrgH4ABgKIdKTB3bvZ/isIArRabYHHdSeFmekCwODBg7Fjxw77/2/fvo127dohIiICaWlpOHbsGHbt2oWAgAApzS3SGAwGhIeH47XXXpPbFJdw8uRJbNy4Ee+//36JfqIpCrjnt1NAYckEMACADsCbjjbycIHNCZvNVijRnTRpEu7fv48dO3ZgyJAh8Pf3x4kTJ/D111/jxIkTqFmzpgusLfr0798fu3btktsMl9C1a1fUrl1b1qodCo7hHtEND3f4VAIYCiAewCYADktp5coFNktuCiu6BoMBfn5+aNeuHTZt2oTly5fjypUraNu2rQusLD5MnDgR9+7dwy+//CK3KZKS5TzzZKexwv9wj+g2agSYzQ6dOhrAKQBfAHC4jqnBADRrZv+vzWbDuXPncPv27QIa6l4yMzMLvATywQcfoEyZMrh16xZUKhVu3bqFESNGuMjC4oXFYkGlSpWwYMECuU2RjCzn2ciRIxXnWRHBPaLbuDGQmZnvaRcBrABwFEAAAMuD4+N82qVmZKDjvHkICgpCQECAff1u7ty5ThruWmw2m8Oi+/333yMkJARDhgxBx44dcefOHQDAZ5995joDiyG9e/fGt99+K7cZkvHMM89Ap9PhzTcdXohTkBn3iG7VqkBk/u6wihCXF1IB3Hvo6JdPO21ICO5HRODatWuIj49HRkYG1Go1GjRo4KzlLsUR0T116hRiYmLw5JNPIjw8HNevX8fatWthsVhQt25dLF682E3WFg8mT56MxMRE/P7773Kb4jSnTp3Cf//7X8V5VsRw329q2jSHlxgKhNkM9UsvYdcPP2Dx4sX2UuQZGRno1asXypcvj2effRbXrl2Tfmwnsdls0Ol0Ob73zz//oE2bNoiKikJmZiaOHz+OnTt3wt/f337OpEmTcOjQIdhsNneZXOTx9fVFSEgI5s+fL7cpTtOlSxfFeVYEcZ/o9ugB1KoFqNXS9alSAVWqAM88A0AUoXnz5kEQBPTv3x/nzp1Dhw4dsGHDBgQFBSEwMBBjxozBlStXpLPBCXJypKWmpmLQoEEoV64cTp06hR07duC3335DVFTUY+179OgBjUaDlStXusvkYkGPHj3w9ddfy22GU6xatUpxnhVVcgvgpdQJb0jywgXSYpEu2Y3JRJ4+/dgw69at4+lHXr9w4QKHDx9Of39/AmC5cuU4cuRIXrx4UdprLADlypXjiy++SFLMkTBjxgzq9Xr6+Phw1apVDvXRqlUrh6tPKIjcuHGDAPjnn3/KbUqhyErbOHr0aLlNUcgFyJ7w5mF++EGanLpGI1nIInsXL17kyJEjWa5cOQKgv78/hw8fXqjsXc7g5+fHmTNnctWqVfTx8aFer+eMGTMKtKNo586dFASBKSkpLrS0+FG+fHkOGTJEbjMKRe/evZWdZx6OZ4kuSe7ZQ5YqJVZ/KKjYarWklxcpUcaoS5cucfTo0QwICCAg1jsbOnSoW2pJeXt709vbmyqVigMHDiy0cJrNZs6ZM0di64o3Y8aMob+/v9xmFJiTJ09SEARu3LhRblMU8sDzRJck//mH7NJFnPWq1fmLrUpFms1k27bkjRsuMeny5cscM2YMAwMD7QI8ePBgnj17VtJxTpw4wZo1axIAw8PDnS7B061bN1auXFki60oGly9fJgBZl5cKQ9WqVVm7dm25zVDIB88U3Sx+/ZXs31+c9VosorBqteJhNouzWp2O/L//E0u2u4mrV69y3LhxLF++PAGwTJkyHDhwoFOJz+Pj49myZUsKgsDo6GhaLBa+8cYbTtt65MgRCoLAW7duOd1XScLf379IrYuuXLmSKpWqyH1RlEQ8W3SzsFrJ334j16wRa6ktXEiuXk0eOSJp8cnCcP36dY4fP55BQUF2AX7mmWcec9blRkpKCvv160eVSsWQkBDu3LmTJOnl5SVZUnFfX18+//zzkvRVUhg6dKjT5Y/cRZbzbNSoUXKbouAARUN0iwjx8fGcOHEiK1SoQAAsXbo0+/Xrx5MnTz52rs1m47Rp06jT6ViqVCmuWbMm2/tms9nhKIX8GDRoEMuXLy9JXyWFs2fPEkCRSIHZu3dvent7F2/nWUYG+d135Jw5YjrYGjXIqCixgviLL5Kff04WEYexIrouIj4+npMmTWJwcDAB0NfXl3379uWJEye4YsUKenl5Ua/Xc/bs2Tn+sRiNRn744YeS2HL+/PkiuUYpN2XKlOHEiRPlNiNPTp06VbydZ0lJ5CuvkGXKiMuJOdVSVKtJb29xCfK550gP/6JURNcN3Lx5k1OmTKGfn5+9PlulSpV4+PDhXNsYDAZ++umnktkQEBDAoUOHStZfSaB///4MDg6W24w8qVq1KmNiYuQ2wzXs3En6+5MGg+MRTDqdKM4ffUQ+UvHbU8hLdJUN2xJx48YNbNu2Df/88w9atmyJcePGITMzE7GxsfD19UWvXr3w22+/ZWtDMtdtwIWhd+/eSgKcAjJt2jRcvnzZYzPSvffeezh79mzx3Hn26qvAU08Bf/8NpKY63i49XSyMMGoUMHgwUNS2weemxlRmug5x/fp1Nm/enIIgsE6dOo9FN9y+fZvTpk1jaGgoBUGgj48Pe/TowSNHjlCn00laCv3mzZsEwOPHj0vWZ0mgVKlSnDZtmtxmPEZaWhpNJlPxdJ4tWCDNJimTiXzmGY+b8UJZXpCe5ORk9unThyqVihUrVuTu3bvzbZOQkMAXX3yRlSpVoiAIBMDGjRvzl19+kcyuSpUqsUePHpL1VxLo2bMnK1WqJLcZj9GnT5/i6Tzbu1fcUSpVOgCzmXz/fbmvKhuK6EqIzWbjCy+8QK1WS19fX3700UeF6icxMZGCILB8+fIUBIFeXl7s1q0bDx486JR9s2fPpsVicaqPksbhw4cpCALv3r0rtyl2Tp8+TUEQuGHDBrlNkZb798ny5aUT3KzDYiGvXJH76uwooisRy5cvp5eXFw0GA+fMmeP0DESlUvHgwYNMTEzkzJkzWblyZbsAd+3alfsLsRkkOTmZgiA4NPNW+B9eXl6cPXu23GbYiYiIKJ7OszfflGZZIaf0ACNHyn11dhTRdZKvvvqKAQEBVKvVHDFihGSl3VUqFY8ePZrttaSkJM6ePZvh4eEUBIEWi4WdO3fmvn37HO63evXqbNOmjSQ2lhS6dOnC8PBwuc0gSa5atap47jzLzCQrVnRIRM8DbA+wFMByAJ8FaHVkfffePbmvkqQiuoXmyJEjjIyMpCAI7NixIxMSEiTtXxCEHDdVZJGUlMR///vfdgE2m83s1KkT9+zZk2e/y5Yto16vl9TW4s5PP/1EQRB4//59We3Icp6N9KBZm2QcPSquvzoguu0BDgSYAvA6wBoA/5NfOy8vcv16ua+SZN6iq4SM5cC1a9fQtGlT1KlTBxaLBX/++Se+/PJLlCpVStJxSNorXeSExWLBSy+9hDNnzuDevXuYPn06/vzzTzRt2hQWiwUdO3bEjz/++Fi7UaNGISMjA5s3b5bU3uJMXFwcjEYjli1bJqsdgwcPhkajwfLly2W1wyUcOuTwqecB/B8AA8R6ie0A5Ftg6d49YN++wlrnPnJTY5bAmW5ycjJ79uxJlUrFSpUq5TujdBYAvHr1aoHbJScnc/78+fZZuMlkYvv27blr1y77OfXq1WNcXJyE1hZ/2rVrx+rVq8s2frF1nmUxeLDDa7TvABwAMBngFYBRADc70rZuXbmvkqSyvJAvNpuNkyZNolarZZkyZfjJJ5+4ZVwATmcGS0lJ4cKFC1mtWjUKgkCj0ci2bdty+vTp1Gg0xS/cyIV89913VKlUTEtLk2X8yMhIRkdHyzK2W2jXzmHRPQmwDkA1xN2dAwFmOtLWQ0L/FNHNg6VLl9JsNtNoNHL+/PluHRsAk5OTJesvJSWFixYtYvXq1YkHH9bq1avz22+/lWyM4o7BYODSpUvdPu7q1aupUqncXr3ErXTo4JDg2gCGAJwDMBXgPwA7A5ziiOiGhcl9lSQV0c2RrVu3sly5clSr1Rw1ahStVqtbx7fZbATgsnHT0tJYpUoV6vV6qlQqGgwGPvnkk9y+fbtLxisutGzZ0u2zzSzn2fDhw906ris5cuQIDx8+nP3zPXy4Q6J788GE4c5Dr215sMSQb/v69eW76IdQRPchDh8+zIiICAqCwKeeekryiARHSUtLo7ik7jq+/fZbqlQqJiYm8o033mDNmjXtAtyqVStu27bNpeMXRT7//HOq1Wq3fgn369ev2O08a926NbVaLXU6HaOjo9m6dWu+17gxrQ6W6KoEcD7EMLEEgF0B9smvnUpFekhO6aIjuhkZZFqaS/ZRX7p0iXFxcRQEgfXr1+f58+clH6Mg3L171+WiS5ImkynbsklaWhqXLVvGWrVqUaVSUa/Xs0WLFvziiy9cbktRQafTceXKlW4Z648//qAgCFzvIaFOUnDz5k22b9/evtU96/i/qChmOLj99wjAZhDjdMsA7AnwRn7tvLxID0l/6bmie/06OW8e2bSpWKgy69tKrSarVCEHDBAr/joxA0hKSuLTTz9NlUrFypUr8+eff5bwAgrP33//7RbR7dKlC6tUqZLje1arlW+99Rajo6PtAty8eXNu3brV5XZ5Mo0bN2ZdN3nBi4Pz7PLly5w9ezafeOIJWiwWAqDJZKJarSYAGo1GfvbZZ+JkqkoVh0S3UIfFIm4z9gA8T3QvXiQ7dxbrouWXR9NiIQMCyHffLdAM2Gazcfz48dRoNCxbtqzHheFkFUZ0NVl5BfJbRrFarXz77bdZu3ZtqtVq6nQ6Nm3alJs3by5Wj72OsH79erdEfrz//vtF0nl25swZTps2jXXq1KHJZCIAlipVik2bNuVrr73GmzdvMj09nXq9nmXKlOGRI0f+1/jddx3eIFGgQ6cjx4+X7Z48iueIbmYmuWKFeNNzyg6f12E2k3Fx5KVLj3V769YttmnThjceVAlesmQJzWYzTSYTX3vtNWmvQSLOnj1LQRDcMlapUqX4wgsvOHy+1WrlihUrWKdOHbsAN2nShBs2bCgRAmyz2ajVarl27VqXjVGUnGfHjh3jhAkTWLNmTer1egJincAnn3ySb731FhMTE3Nst379+sfj0NPSxAgDqUXX25uMj3fD3XAMzxDdzExy3Djnkl1oNGTp0uSpU9m67ty5M9VqNVu1akU/Pz9qNBqOHTvW7REJBeHEiRNuE90BAwYwKCioUG0zMjK4atUqxsbGUq1WU6vVslGjRly/fn2xFuAGDRq4dHNJv3796OXl5ZH38Oeff+aoUaMYERFBrVZLQRBYrlw5duzYke+//77TW6Xv7d7NtIJOuvI6TCbyv/+V6OqlwTNEd/p0abILCYIovJcvkyTXrl1Lg8HArMX61q1b5/rN60kcPnyYKpXKLWP99ddfBMDLD+5ZYbHZbFy9ejXr1q1rF+C4uDiuW7fOI8XDGdasWUOtVuuS6zpz5ozHOM9sNhu/++47Dho0iGFhYdRoNBQEgUFBQezevTvXr1/v9GaR9PR0Hjx4kIsXL2ZsbCwBcElEBDOlyKlrMpFjx0p0N6RDftH96SdpkxZrNGSTJjywf/9jHtIOHTpIY7OL2bdvH9VqtdvGK1eunKSPsjabjWvWrGH9+vWp0Wio1Wr5xBNPcO3atcVCgG02G9VqNTdt2iR535GRkaxVq5bk/TqCzWbj1q1b2bt3b4aEhFClUtkT8fft25dffPGFpL+/zMxMlitXLptjTaPRiLPl5cud0wWTSSxS6WFVI0i5RTctjaxQQTrBfXBkGAwcodNRo9EwNDSUderUYePGjTl69GjnbXYDu3btokajcdt448aNY9myZV3St81m40cffcQGDRpQo9FQo9GwQYMG/OCDD4q0ANepU4fNmjWTtM81a9a41XmWlpbGTz75hN26dbMnzNdoNKxcuTIHDx7M77//3uW/ozfffNMuuGq1mjNmzPjfmwcOkCEhBXsKNhhIX1+xJLuHIq/ofvKJGIEgsegSYHq5ch75LecI33zzDbVardvGi4+PJwD+/vvvLh3HZrPx448/ZsOGDanVaqnRaFivXj2uWbOmyAnw22+/LWmKTKvVSrPZzGHDhknW56MkJydz1apV7NChA/39/SkIAnU6HSMjIzlmzJhCJcZ3hvT0dNavX5+CIFCr1dJgMDD+UYdXSgq5bJkovhZLjrNfm07HNINBLNP+yivk7dtuvY6CIq/oxsQUSEjPANQD7OfI+RYLWUTzCnzxxRfU6XRuHbNixYr8v//7P7eNZ7PZ+OmnnzIuLs4uwHXr1uWqVauKhABbrVaqVCrJdu4NGDCAXl5ekjp4ExISuGzZMrZq1YqlS5cmABoMBtaqVYuTJk2StUjp1atXWa5cOfr4+PDEiRNs1KgR+/Tpk3uDzEzy0CHyrbfI3r3JJk3IRo3Irl15acIEtjEYmCZhrhJXIp/o3rolxs8VQHRbA2zsqOgC5IgRztkoExs3bnR7ovGXX36ZXl5ebh0zC5vNxv/+979s3LgxdTod1Wo1Y2NjuXLlSo8W4Jo1a7J169ZO95PlPPv000+d6ic+Pp4LFixgkyZN6OPjw6yNCLGxsXzxxRd59uxZp22Vgj179tBgMDAiIoJJSUkkxUiYwlZd2b9/PwGwVatWzCwCT7fyie6335I+Pg4L7icQt/vNLIjoRkU5Z6NMrFu3jgaDwa1jJiUlURAE7t27163jPorNZuPGjRvZpEkTuwDXqVOH77zzjseF+b3++us0Go2Fapuens5du3YxMzOT1apVK5Tz7OLFi5w1axYbNGhg3+3l5eXFhg0bcs6cObziQcUYs3j77bepUqnYtWtXyb5Q33rrLfua8IsvvihJn65EPtF9/XVx15kD4pkIMBzg5YKKbhEtS7NmzRqaTCa3j1utWjW2a9fO7ePmRmZmJrds2cJmzZpRr9dTrVYzJiaGy5cv9wgBTk1NpSAI2RLEO8oPP/xAAPaCo444z/744w9OnTqVtWvXptFoZNZur2bNmnHRokW8efNmIa7CfQwfPpyCIEhe5LNHjx7MilAymUxctWqVpP1LjXyi+8orYlytA+I5HuCrD34ukOi6aYOB1KxYsUKWUulvvPGG22fYBWHr1q1s3ry5PSVldHQ0ly1bJqsAR0ZGsmPHjgVut3HjRvvsVK1W87nnnmNGRka2c44cOcLx48ezRo0a2XZ7tWnThsuXLy8SMeekOKvPil5xRe6OChUq2EVXo9GwX79+ko8hJfKJ7rx5YvKafITzCMDqANMKI7pujHWVkmXLltHb29vt42Y5h4pCUpsvvviCLVu2pMFgoEqlYs2aNbl06VK3V3aYO3duob4gV6xYYQ+V0ul0DAgI4LZt2zhy5EhWrVo1226vTp06cc2aNUxJSXHBFbiWa9euMSAggD4+Pjx9+rRLxujWrRufe+45qlQqbt682SVjSIl8ort2rZhuLR/hXALQBLHUcjmAZoAGgLUdEd2AAOdslInXX3+dPj4+sowdGxvLRo0ayTJ2Ydm2bRtbtWqVTYDfeOMNtwhwcnIyBUHggQMHRA/7r7+Sb79N9u9PNm9ONmtG9uxJvvGGuBHowWx2woQJBEBBEOjt7U21Wm3f7dWjRw/+97//9YglFGfYu3fvYw4zV1K5cmX26tXL5eM4i3yie/q0QzG6yRDLLGcdzwPsDvBvR0S3iOxAe5QFCxbQ19dXlrHXrl1bpOunbd++na1bt6bRaKRKpWJUVBQXL17sUgGuFhbGt6OjydBQMflSDsH8mQYD041G3tHr+aq3N70fPA4HBgayX79+/Oqrr4rsPc8JVzjM8mPkyJGFziPiTuQTXZvNoZnuo4fDywsmEylDPSsp+Pe//80yZcrIMnZWFq333ntPlvGlZMeOHWzTpg2NRiMFQWD16tW5aNEiaR/TDx3ibR8f3ivAZzhVreYtQeAdJ0PEPBVXOczy49ChQxQEweOfEOTdHDF5coFjdR0+DAaP35mSGzNmzKC/v79s4zdr1owxMTGyje8KvvvuO7Zr144mk4mCIDAyMpILFixwToDfeYc0Gh2rRJvbxGDy5CK7c/JRXO0wcwSdTsePP/5YlrEdRV7RvXBB2mQ3WYdWK66pFVGmTp3KABnXo7/++mtZy427ml27drF9+/bZBHj+/PkFqr6cvGQJ07Va5z+rJhM5aZILr9Y9uMNh5gg1atQoVDSJO8lLdFVwNRUrAhMnAiaTpN3SYAAWL5a0T3eSnp4OtVot2/jt2rWDwWDAG2+8IZsNrqR58+bYtm0bkpOTsXv3blSuXBlz5syBxWJBREQE5s6di/v372dr8/fff2PBggVo0qQJmloswMSJ0Fqtzhtz/z7wzjvAli3O9yUTP/30E8LCwuDj44MrV64gIiJCNlvatWuHAwcOyDa+0+SmxpRqpkuKmcaqVnUofMyRI12r5YhSpThjxgyePn26SDonxo0bx+DgYFlteOqppxgeHi6rDe5m79697NSpE81mMwGwdOnSDAoKsv/f29ubTZ94grd9fQu/pJDb4eND/vOP3LegwLzzzjtud5jlxZ9//kkAbomWKCyQdaYLADod8P33QNmygEbjXF8mEzTTp+O70qXx73//G7Vq1YLFYkG9evXw5ZdfSmOvG7BarVCp3HP7c2PGjBk4e/Ys7t69K6sd7uKPP/7A559/jitXriAzMxMAkJSUhBs3biA5ORmhoaGYNGkSvunfH77p6RCkNiAlBZg3T+peXcqIESMwevRozJw5E1u2bJH9MwsAVapUgdFoxIcffii3KYUjNzWmlDPdLK5cISMiCldBQhCYqlLxrzFjSJLnzp2jRqMhIIblqNVqybJBuYNhw4YxLCxMbjPo4+PDqVOnym2GS/j11185btw4RkVF2Xd7lS1blm3btuU777yTbaa0f/9+du3alV4WC8/k8hlcBjAWoA7gwEfe+w5gBEAjwOYAL+T2OfbyIlNTZbwrjvGww+yzzz6T25zHqFu3Llu0aCG3GbkC2StHPEx6OjljhuhcczAvA728yKpV+c6YMQTAoUOH8p9//mG3bt2oUqkIgBaL5fEieB7MwIEDPeLRvl+/frIvc0jF3r17OXz48Gy7vQICAti5c2d+8MEHjkUxHD/OjFwqVG8CuAXgqEdE9yZAb4D/BZgCcDLABnl9lj1QxB7GUxxmefHKK6/ItrnIETxLdLO4ckWsm+brKwqwt7cYWqbRiDNhb2/x58aNyS++IDMyeOfOHfvM1mg0cuDAgQTAJk2asGrVqtTpdEVmttu3b19GRkbKbQbPnDlDAEXqC4sUY423b9/OZ555hpUqVbLv9qpQoQJ79OjBjRs3Fi6W87338i0RPv0R0V0BsOFD/78HcUflqZzaq1RkASozu5usHWZVq1b16DXT69evEwD//vtvuU3JEc8U3SwyM8nz58mNG8klS8iFC8Utlj/9ROYQ3lOmTBm78KpUKsbFxTEhIYGkmCRaEAROmzbN9XY7Sc+ePVmjRg25zSBJ+vv7c9SoUXKbkSdWq5WbN29mz549GRwcTJVKRbVazdDQUPbv35/btm2TxskzdGi+T16Piu74B7Pfh8+JArgxtz6eeMJ5O13AihUrqFKp2KVLF49wmOWHt7c3582bJ7cZOeLZoltAmjdvbhfd+vXr8969e9neX7VqFdVqNRs3buzRMajdunVjdHS03GaQJJ999ln6+fnJbUY20tLSuHbtWnbp0oWBgYH22l7h4eEcNmwYd+/e7RphaNeuwKI7BOC/HjknDuD7ufVRubL0djvJiBEjZNlh5gxNmzZl/fr15TYjR/ISXfldkQWkYcOG0Gq16NChAw4fPoz4+Phs7w8dOhRHjx7F8ePHERQUhL/++ksmS/PGarXKGqf7MDNmzMDNmzdx6tQp2Wy4f/8+3n33XbRv3x7+/v7Q6/UYMmQI/vzzT3Tv3h0HDx6E1WrFmTNnsHLlSjRt2tRpT/r58+dx79697C+SBe7HAuDR+I+7ALxya1CIMVxFRkYGnnjiCaxevRpbtmzByy+/LLdJDtOjRw+cOHFCbjMKTm5qTA+d6d66dcueDDomJoYBAQE5zniSk5NZp04dajQarl+/3t1m5kv79u1Zr149uc2wExISwt69e7ttvISEBL7xxhts0aIFfX197ev00dHRnDx5Mk+cOOFyGyIjI6lWq1mvXj0uWrSIv//+O+9261bgme6KBzPbh9d0jbmt6QKkh/zer1+/7vEOs7y4e/cuAfDcuXNym/IYKE7LCw+TlJREk8nELl265HrO2LFjKQiCx5Vmb926NRs2bCi3GXamT5/u0vy+169f5/z589moUSN6e3sTAM1mM+vVq8cZM2a4/Q/nzp07rF27tn2pKusYr9fTmktUjRVidMJUgP0f/GyFmA3P+8EabgrAF5B79IIN4D0PqOtXVBxm+VGmTBn+61//ktuMxyi2okuKJVEEQeDq1atzPWfDhg3UaDSMiYkp0N57V9KiRQs2adJEbjPsZNVP27dvnyT9XbhwgTNmzGC9evWy7fZq1KgR582bx2vXrkkyjiOcPHmSCxYsYKdOnRgaGmqP2c1KMA6AWq2WU6ZMYeYvv+SajnTmIwKNB68R4LcQ43QNAJsBPJ+L6N6FWAcwKCiI48eP5/Xr1912H7Ioag6zvGjfvj1r1qwptxmPUaxFlySnTJlCjUaT52zp3Llz9PPzo7e3N48dO+ZG63KmSZMmHhfcHRERwQ6FzE988uRJTpkyhdHR0fbaXr6+vmzRogWXLFnCW7duSWzt46Snp3P79u2cMGEC4+Li6OfnZ4/j9vHxYXR0NAcPHsyPP/6Yd+7c4fbt26nRaFi2bFn+/PPPYieZmWSFCvkuMRT6MJn459GjHDRoEMuWLUtAzLc7duxYt4TtZTnMZs2a5fKx3MEHH3xAnU4ntxmPUexFlySjo6MZGBiY5ze31Wpls2bNqFaruWLFCjda9zgNGzaUpLS3lBSk8u3hw4c5btw4Vq9e3T5z9PPzY7t27bhixQqXP7LGx8dz5cqV7NOnD6tXr26vRabRaBgYGMjmzZtz6tSp3LNnT66ficTERD733HP859F8CP/5T+F2TeZ3aLXkyJHZhjp79iwHDx5MPz8/uwCPGTOGly9flvR+Wa1Wj95hVljS0tIoCAKPHDkitynZKBGim5iYSJPJxG7duuV77ksvvURBENinTx/ZHq/q1avH9u3byzJ2bqSlpVGlUvGLL77I9rrNZuOePXs4bNgwhoeHU6PR2Hd7denShR999JHLanvZbDYePnyYs2bNYps2bRgcHEytVktArAobHh7O7t27c+nSpQ5V23WIe/fIwEDpRddiETcF5cK5c+c4dOhQ+vv7EwADAgI4atQopwW4qDvM8iMwMNDjfDYlQnRJMYeqIAhcs2ZNvufu2LGDer2e4eHh9s0V7qROnTp86qmn3D5uftSpU4eNGjXitm3bOGDAAIaGhlKtVlOlUjE4OJg9e/bkpk2bXJK5Pzk5mZs3b+aoUaNYt25dli5dmoIgUBAEli5dmnXr1uWoUaO4efNm16/N794tbR5os5nMw+/wKBcuXODw4cPtAlyuXDmOHDmywF8sxcVhlhc9evTwiC31D1NiRJckn3/+eWo0Goc+nNevX2dwcDBNJpNkDiRHiY6OdmhW7g6sVis3btzInj172sO31Go1K1WqxAEDBnD79u2SPxFcvHiRS5cuZffu3RkeHm5fB9ZqtQwODmabNm04a9YsHj58WD5nzyuvSLPMYDKJCfcLWT3i4sWLHDlyJAMCAgiA/v7+HDZsWL6f8eLkMMuLzz//nGq12qOusUSJLknWrFmT5cuXd+iXYLPZ2KlTJ6pUKi5cuNAN1olERUWxZ8+ebhvvYVJSUvjRRx+xc+fODAgIoCAI1Gq1rFq1KocNG0a1Wu3Q04IjZC1NTJ06lc2bN2dgYKA9O5zFYmH16tXZt29frly5kvHx8ZKMKSmzZzsnvGYzOWCAvUKws1y+fJmjR4+2C7Cfnx+HDh36mBN55MiRxcphlhc2m40qlYrff/+93KbYKXGim5iYSKPRyO7duzvcZtGiRVSpVOzQoYNbvjEjIyPZt29fl49DiuFg77zzDtu2bWv3mOv1ekZFRXHcuHE8fPhwtvObNm3KOnXqFHicxMRErl27loMHD2Z0dDR9fHwoCAJVKhX9/PwYFxfHCRMmcPv27R69RftRbF9/zTtaLW2OZsXLcppZLOQHH7isPtrly5f57LPPMjAwkICYtvKZZ55hTExMsXOY5UdoaCgHDBggtxl2SpzokuTOnTspCAI/+OADh9vs27ePJpOJFSpUcHkcaXh4OAcOHOiSvm/dusUlS5awefPm2XZ7xcTE8IUXXuDJkyfzbP/VV1/lWz/t9OnT9tjXSpUq2SMY9Ho9K1WqxE6dOnHBggX5juXpxMfHs2LFivQ3GGh76SUxK56Xl5gB7xGhzVSpxPfMZnLsWNKNmduuXr3KIUOGZAuRGzhwIM+cOeM2G+Rk0KBBDAkJkdsMOyVSdEly4sSJDq/vZpGQkMDw8HDq9Xpu377dZbaFhYVx2LBhkvR17do1zps3j3Fxcdl2e9WvX58zZ84slFffaDRy8eLFtFqt+ca+Dho0iGvXrpXFIekqMjMz+d5779k3dtgTq1it5Pbt5KxZZIsWZPXqZLVq3KdS8UDLluSWLTlmx3M1DzvM/vzzTz733HMMCgoiIJYkGjBgQLGMXMhi7969FATBY9Z1S6zokmLl0KCgoAL9Mmw2G/v160dBEPjiiy+6xK6KFSsWOszl3LlzfOmll1i3bt1su70aN27M+fPnF3qXU3x8PFetWsW+ffvSYrFQEAQWNPa1OJCRkcHGjRvTZDIxa+fZc889l2ebJ598UradUVkOs86dOz/2e4mPj+fEiRNZoUIFuwD369evyD+B5IRGo+HGjRvlNoNkCRfdhIQEGo3GQjmtVq5cSbVazSZNmki+BhkcHMzx48c7dO7vv//OyZMnMzo6mgaDwf7H07JlSy5durTAM0ybzcZff/2Vs2fPtse+6nQ6+zJEeHi4PYXmb7/9VoirK9pkZmby+eeft8cDq1Qqzp8/P882WSXtXRFKlxcFcZj9/fffnDRpEoODg5m1Y7BPnz5uSS7kDiIjI9m1a1e5zSBZwkWXJL/77jsKgsC1a9cWuO3x48fp4+NDPz8/SZOylC9fnpMmTcrxvUOHDnHs2LGsVq2aXQz9/PzYvn17rly5skDxlvfv3+eWLVs4evToAse+ent7u2ym7+lYrVaaTCaWLVuWgiBw5cqV+bbR6/V8++233WCdaN8TTzxRaIfZzZs3OWXKFIaEhBAAS5Uqxd69exfpL9nx48ezXLlycptBUhFdkuIvRKPR8OLFiwVum5yczNq1a1Oj0fC///2vJPYEBARw6tSptNls3L17N4cOHZptt1dgYCC7du3KtWvXOjzLfjj2tWrVqvbH46zY19atW3PmzJkOx7726dPHo5wT7mTgwIG0WCxMT0/nl19+yRs3buTbplmzZoWK+igoWTvMvL29JVmnvXXrFv/1r3+xYsWK9rX6nj17ekSOEkexWq385JNPCICxsbGsWLGirPYoovuAqKgoVqhQodDrkWPGjKEgCBzzoCJxYbDZbNy2bRv1ej19fHyy7fbq1asXP/vss3wfUbNiX6dNm5Zr7GufPn2cjn09ffo0AciSCUtOzp49S5VKVeAno02bNrk8SP+nn36iwWBgeHi4S3aYJSQkcNq0aQwNDaUgCPTx8WGPHj08LrfBo7Rs2dIeQQOAlSpVktUeRXQfkLW+26tXr0L3UdA0kVarlRs2bGCPHj0YFBREQRDsQhsTE8MdO3bk+UeamJjIjz/+mIMHD2ZMTMxjsa8NGzZ0aeyrn5+fU18yRZGoqKhC1a+z2WzUarWSbSx5lLwcZq4gISGB06dPZ1hYGAVBoLe3N59++mkeOnTI5WMXlP3792dzfPbr109WezxfdO/fJ9etE4sCRkaKsZDe3mRAAPnkk+TcuaRE3tYdO3ZQEASuW7eu0H1kpYn08fHh8ePHs72XkpLCDz74gJ06dXpst9fIkSO5d+9ekqSvry8XLFiQre3p06e5cOFCe+xrltNMr9czNDSUHTt2dHvs66hRo+jv7++28eTmo48+okqlKvT6fVxcHBs0aCCxVf9zmM2cOVPyvh0hMTGRM2bMsAuwl5cXu3XrxoMHD8piT05s3rzZ/sTnyBq8K/Fc0U1IIJ97Tty54+XFPHf3GI1k3brkd985Pey4ceOo1Wqdyt6Unp7OJk2aUBAE9uvXj23atMm226tGjRocP358jo9lWU6aZs2aMS4ujv7+/tliX2vVquUxsa/Xrl0jgBIRZG+1Wmk2mzlo0KBC97F27VpqtVrJZqLOOsxcQWJiImfOnMnKlStTEARaLBZ26dKF+/fvl9s0zpgxgwDcnkvlUTxTdL/+mixdmizI1sqs5CEDBpCJiU4NX61aNQYHBxf4j+Off/7h4sWL2bx5c5YqVcr+OOPr68sXXnjhMcfG33//bY99fTjva1ZsbbNmzTh16lTXVbeVgODgYLdtWZaTLOeZM2FfNptNsrp8DzvMTp065XR/riApKYmzZ89meHi4XYA7d+4sn+hducJZvr78Oy5OTEbv7S0eoaFkz57k22+Tj+ZPdgGeJ7qvveZcEhGDgQwLI51w8Ny6dYsGg4F9+vTJ87yrV69yzpw5bNiwIb28vOzOqgYNGnDWrFm8ePEit2/fTq1Wy5CQEE6bNo1t27bNMfb16aef5htvvMELFy7QYrHw3XffLbT97mTatGn08fGR2wyXUljnWU7UrVuXjRs3dqqPffv2udRh5gru3bvHOXPmsGrVqhQEgWazmR07duSePXtcP/jx42TbtqTBQJtOl/ekzWAQBfj8eZeZ41mi+9Zb0qTL02jISpXEJYpCsn37dgqCwE8//dT+2rlz5/jiiy8yNjbWvjDv4+PDJk2acMGCBYyPj2dKSgo/++wzjh49mvXq1bPHvmbNYCMjIzly5Ehu3LgxV2ebyWRymcNFahITEykIAg8cOCC3KS4jKiqKUVFRkvS1cuVKp0rIuNth5gqSk5M5b948RkRE2AW4Q4cO3L17t7QDWa3ilmyjkRQEx/VDrRZ16K23XJKQyHNE9/hxaRND6/WkE5EIpBiLqlKpWK1atWy7vVq1asVly5bxt99+47Jly9ijR4/HYl8rVKhgj309dOgQrVYrO3ToQJVKxUWLFuU5rtFodMqZ526qVq3KTp06yW2GS8hynv3111+S9Ge1WqlWq/n5558XuK3cDjNXkJyczFdffZWRkZEUBIEmk4nt2rXjrl27nOs4JYVs2dL51Jt9+0qWejMLzxDdjAwxMqEg30aOrvF+9VW2oW7fvs0RI0bw9u3bj5lx8OBBjhkzhpGRkfbHf7VaTb1ezxdeeIGTJ09+LPbVbDazWrVq7N27N1esWJFv3OrChQupUqnYsWPHXGcqer3eY/aJO8Jrr73mcP20ooQUzrOciI6OZsuWLQtkR8OGDT3KYeYKUlJSuHDhQlarVo2CINBoNLJt27bcuXNnjudnZmYyM6eZaEaGGNkkxSTOZCIl/v17huh+9lmupa2dPqpVsw9z5swZli9fnoIgcMqUKdy1axcHDx7MKlWqZKvtVbduXTZr1ozR0dH2zFxZ220bNmzI8ePHc9u2bYWOff3pp59oMpkYHByco0jrdLrHapF5Mln107Zt2ya3KZIyaNAgms1myXMmLF26lAaDwaFzH3aYFcdENLmRkpLCRYsWsXr16lSpVDQajWzdujV37NhhP+ett95i7dq1eefOneyNnfUL5SS8Eu02JT1FdBs2zPFilwGMBagDOPCh19MAdgdY8YEY7srvhh06xO+//z7brhQAFASBpUqVYvny5VmuXLlseV9DQ0PZoUMHzp8/n2+//TYFQZDE65zFw2kiH/4gkaRWq33sNU8nJiaGzZo1k9sMyZDSefYoWV9S3+UT4lgUHWauIC0tjUuWLGGNGjWoUqloMBj45JNPMjg4mBqNhjVq1Pif8P75p7TLlFmHtzf599+SXI/8opuQIMba5nChmwBuATgqB9FdAnAPwIB8RDdTpeKnD3KH5nR4e3uzVq1aHDhwID/66KNcY19Hjx5NrVbLqxImn7bZbOzbty8FQeBLL71kf12j0fCHH36QbBx3sHr1akljUOWmRo0akjnPcqJ69eps165dru9nOcyeeuqpYnNPpSAtLY3/+c9/GBERkW3yVLFiRVF4+/cXHWFSi65eT86YIck1yC+6O3eK3yJ5XPD0R0T34SMov5kuwF8eih54+Cjo7CEiIoIVK1aU/I8gK01k06ZN7Y6Wn3/+WdIxXE1WDOpHH30ktylOI7XzLCcWLFhAk8mU43ujRo0qdg4zqVm8eDHVarVddAHQX6ejNZcJ3KNHM4B6gOYHR1VHhNfXl0xPd9r2vERXBXfw229AaqpLh6hjMsFqtWLv3r2YMmUKKleuDAC4c+dOgfrZu3cvbty4gUGDBklq37Bhw/Drr7/i6NGjKF++PDIzM6HT6SQdw9WoVCo88cQTeOONN+Q2xSkyMjIwatQoPPPMMwgLC3PZOOPGjUNKSgp++umnbGPHxcVh1apV2LJlC2bNmuWy8Ys6Wq0WLVq0wPPPP4+VK1di586d+GrsWFCjcbiPNwHce3D84UgDmw34+efCGewouakxpZzpzpqV7zeMszNdqtWPDVvYNbKvvvqKgiBww4YNzl75YyQnJzMmJoYAuHjxYsn7dzWff/45VSoV0yWYDciFq5xnOREeHs4uXbqQLLkOM0kZO9bhCKhmAFcWZonh9dedNhOyz3RVKkAQXDtGDv1bLJZCddWhQweMHDkS/fr1w7Vr15y1LBsmkwlHjhwBAEyePBnjxo2TtH9X89RTT0Gn0+Gtt96S25RC8ddff+HDDz/EihUroCnAjKmw9O/fH99//z1+/vlnVKpUCV5eXrh69SqqVavm8rGLJfv3i/LoINMAlAXQCMAPjjRISwP27SuUaQ6TmxpTypnuqlX5hnc4PdP185PG1ocIDw9naGioS5wcALh8+XJqNBrWqVPHoTSRnkL79u0ZGRkptxmFwtXOs0e5e/cuAbHkT6dOnRSHmbOEhzs8a90P8C7AVIBrAFoAnnWkbYsWTpsJ2We6deoAanWOb2UASAVge3CkPngNANIe/B8A0h/8nOt3XJ06Ehn7P/bt24fr169jyJAhkvcNAD179sTp06dx6dIllC9fHidOnHDJOFLz0ksv4Y8//sC9e/fkNqVAfPzxxzh58iS2bt3qtjFfeOEFAEBERAS++OILqFTu+ZMrthTg/jUA4AVAD2AgxNnuNkcauvgJyD2fgKgowGrN8a05AIwAXgWw9sHPcx68F/Hg/1cBtH3w88WcOjEYgNatpbUZQNmyZbFx40Z8+OGH2LRpk+T96/V6VK5cGVevXkXNmjURExOD9957T/JxpCYuLg4WiwULFy6U2xSHedh5luVkdfV4WQ6zp59+GtevX3f5mCWCChUK3VRAHpO2hwkNLfQYDpHbFJhSLi+Q5ODBromtA8SsQQ7UsCosw4cPp06nk7RsDYDHnFHTpk2jIAjs37+/ZOO4il69esleh6oguNN5dv36dQYGBtodZrdu3SIA/v777y4fu9gze7aY7CofTUgAuB1gCkArwLUATQD/yE9LLBby/fedNhOyx+mS0ie7yTrUarJzZ+nszIXw8HBWqlRJkjU5q9VKADnuKd+2bRt1Oh0jIiJkT2CeFydPniQAp2qwuYtz5865bOfZo+S2wyw4OLhIfJl6PDt35l3w4MHxN8C6D9ZxfQA2ALjDET0xmUgJin16huiS5MCB4qxUStE1m8kLF6S1Mwfi4+Op1+slSYySnJxM8SEjZ65evcqgoCCazWaPyMafG2XLluXYsWPlNiNfatSowerVq7t8nKwdZjk5zCZOnMgyZcq43IZiT0YG6e/vmidmgKxVSxIzPUd0796V9oaZzeTy5dLamAdbt26lIAjcvHmzU/1kPW7mhc1ms6eJ9NR43uHDh7NcuXJym5Ena9eupUql4tmzZ106TtYOsxm5bCO9fv06Abh0B1yJYf581zw1m83kxx9LYqLniC4pLjPksyXYocNkIocNc0kC4rwYOnQodTodbzixhnz16tV8RTeLBQsWeOz+/MuXLxOAywWtsFitVlosFg4cONClY2SlZMzvyzgwMJBDhw51mS0lhvv3yaAgaQVXpSJjYiTLq+tZokuSx46RZcsWvD7aw4I7bpzbBTeLypUrs3LlyoVuf+7cOQqC4PD5e/fupdFoZHBwsFNi7wqCgoI4YMAAuc3IEVc7zx51mOXHmDFjSlRlZZfy00/SznZNJlLC4queJ7qkmHmsd2/xYlUqx6f/ZcuS33zjOrscID4+njqdrtCzlpMnTxZIdElxSaJKlSrU6/X89ttvCzWuK3jhhRdYqlQpuc14jCznmauS8zzsMLt7965DbS5dukQATlWhVniIt9+WJqeu0UgWospHXnim6GZx4ADZo4c46/X2zj77FQTRU2kykRUrivWMPCTn6GeffUZBEAqV5f/IkSNUqVQFbmez2di7d28KgsCXX365wO1dQUJCAgVB4KFDh+Q2JRuudJ7l5TDLD39/f44ZM8YldpVIVq0qeH20rEOtFidyLigm4Nmim8Xdu+QPP5CLF5NjxojrtS+8QH7yiTjtl2kpIS8GDx5MnU5X4LCpAwcOFEp0s8j6o2/WrJlb4k7zIzw8nJ3dELbnKK50no0ePTpPh1l+DBkyhOXLl5fYqhLO8eNiKbCCVKaxWMTCCi6KfCoaolsEsdlsDAsLY5UqVRw6Pzk5mZUrV2bp0qUJgCEhIezbt2+hxj527Bi9vb3p7+/PC24ImcuLV199Nde8se4my3n2zDPPSN5vXFycQw6zvDh79iyLSnxzkcJqJTdsIGNjxbBUb+/sm7G0WtLHR3ySbt6c3LbNpRM5RXRdyPXr1x1e383MzGStWrUIiAnWdTodp0yZUuixk5OTGR0dTa1WK2uRy9TUVKpUKo8oPzR48GDJnWcFdZjlR5kyZThp0iQJLFPIkatXxTXamTPJkSPJ0aPJOXPIr7+WrBxPfiii62I2b95MQRC4detWnj59mmvWrMn13G3bttlLvZvNZkl2nWXFiI4bN87pvgpLdHQ0W0iQnckZXOE8e9hhlpiYKEmf/fr1Y0hIiCR9KXgmiui6gUGDBlGtVtNgMNBiseRcNpribDcwMJAAOGvWLMnG//TTT6nRaBgbG8uUlBTJ+nWUlStXUqfTyRpLLLXzzBmHWV4cP36cAHj79m3J+lTwLPISXSXPnASkpKQgPT0dmZmZSE1NRUZGBm7cuJHjuYIgYMCAAQCA559/XjIbevXqhdOnT+PChQsIDAzE77//LlnfjjBkyBBkZmZi/fr1bh03i6y0jZ9//rkk/Y0ZMwajRo3C9OnTJU/JWKNGDfj4+GDRokWS9alQhMhNjanMdB1m37591Gg01Gg09rXa7du3/++EzEzyl1/IN98ke/ViUnQ0j5pMZKdO5Kuvikk8JCp/k5aWxkaNGlGtVnP16tWS9OkocXFxrFu3rlvHJKV1nknlMMuPHj16MCwszGX9K8gLlOUF13PlyhUOGzaMWq2WAMStp6mpYmxxxYpiPGBOO2h0OtHTWqoU+eKL5M2bktgzdepUCoIguRc/L7Zs2UK1Wu32MDapnGdSO8zy4pdffqEgCIWu46fg2Sii60YuX77MqlWrsqnFQltoaMF2zOj1YliLRJEIWWkiIyMjJXMC5YfBYODSpUvdMhb5P+fZhx9+6FQ/+/fvp9FolNRhlh9eXl585ZVX3DKWgntRRNfNZK5axXSNhpnO7AMfPpyUwHnzcJrIgwcPSnB1edO2bVu3pFHMombNmk6Pt3LlSqpUKnbs2NGtjsDOnTuzatWqbhtPwX0ooutOVq+WZj+4ySRW25AggNtms7Fdu3ZUqVRcsmSJ89eYBz/++CMFQXBLoc1169Y5vfPM2R1mzrBnzx4KgiBLtImCa1FE11389pu0mY/MZklKh2Qxf/58qlQqdunSxaUzOi8vL0nD4XLCWefZww4zOTeWmEwmLly4ULbxFVyDIrruID2djIgoXOKN/PaIX7kimZl79uyh0WhkSEiIy7ai9ujRg6GhoS7pOwtnnGfudJjlR7t27dxaEl7BPeQlukqcrlSsXw9cuSJKpZSkpgIvvyxZd40bN8aVK1eg1WoREhKCnTt3StZ3FjNnzsSFCxfwzz//SN43AJw/fx4ffPAB3n77bWgKWC77wIEDCAsLg8ViweXLl1GtWjWX2Ogozz//PE6dOoX09HRZ7VBwI7mpMZWZbsGoUcOhmWs/gAEAvQCGA1zpyGzXaBSzsEmIzWZjr169KAgCZ86cKWnfpJhfYPz48ZL3SxbeeSaXwyw/DAYDly1bJrcZChICZXnBxZw96/Ba7gmAqQ9+PgWwHMBfHFlicFEy7nfeeYcqlYotWrSQNL526NChDAgIkKy/LNatW0dBEArsPBszZoxsDrP8aNGiBWNiYuQ2Q0FCFNF1NevWOVQW+tHj9INZ73pHzh850mXmHzlyhF5eXpKmicyqknDu3DlJ+iML5zzzFIdZXmzdulWWTSUKriMv0VXWdKXg4EHg3j2HTx8DwAQgEkAggA6ONNq3r1CmOUJMTAyuXbuGgIAAhIeHY8uWLU73GRwcjPLly2PWrFnOG/iAESNGgCTee+89h86/ceMGQkJCcPz4cRw9ehTdu3eXzBYp6dy5M9RqNT744AO5TVFwA4roSsHVqwVyoC0HkARgD4CnAegdaXTrVrb/3rlzB19++aX4uCIBFosFx44dw5AhQ9C9e3dMnDjR6T779u2LL774QgLrCu48y3KYmc1mXLlyBVFRUZLY4Spq1KiBl19+GVWrVsW0adPkNkfBleQ2BaayvOA43bsXOiRsJMD/OHBeZmAgU1NTuXnzZrZp04Z6vZ4AXLJ3f926ddRoNKxbt65Tgfu3bt2iIAj89ddfnbapIM4zT3WY5URiYiIjIiKo0+kIgIIguDzGWcH1QFlecDFBQYVumgHgLwfOOxUfD4PBgKeffho7duxAWloa1Go1li1bhq+//lrS8Kw+ffrg5MmTOH/+PAIDA3Hq1KlC9VO6dGmEhYXhlVdeccqeTz75BL///rtDaRufffZZjBgxAtOnT8eXX34paUpGV2CxWFC9enX77F2v1yMyMlJmqxRcSm5qTGWm6zhr1zpUFC8e4CcAkwBmANwO0ARwqwMz3YzBg/nqq6/SaDRSpVIRADUaDX19fe2ZzQRBoMlkYvny5VmnTh1269aNU6dO5ccff8wzZ84UeNaXlpbGuLg4qtXqPKth5MXcuXNpMpn4+eefc/r06bkmd88NR51nRcFhlhtWq5VPPfUUBUGgSqXikSNH5DZJwUmgRC+4mDNnHMq38DfApgB9HsTp1gD4riPLEBYL+UD0EhIS2LdvX6rVajZt2tRugtVq5dGjR7lq1SpOmDCBHTp0YM2aNVmuXDl7eSAA1Ov1LFu2LKtVq8Y2bdpwzJgxfOutt7h//36mpaXleHkvvPACBUEQ01UWgCNHjrB37972cY1Go8Nt169fz/j4eA4dOjTfnWfx8fEMDAykl5cXT5w4USAbPQWr1Wqvn2dfMkpKEmt63bnjkdWwFXInL9EVmIcjpm7duvzll1/cMN8uBkRGAn/84Zq+DQbgxg3Ax8f+0rfffgubzYZ27do53M2lS5dw6NAhHDt2DKdPn8bFixdx/fp1JCQk4P79+8jMzIRGo4HZbEaZMmVQvnx5hIWFITIyElarFXPmzEGVKlWwf/9+eHt75zves88+ixUrVsBmswEAwsLC8Ndf+S+mkIROp4NOp8P9+/fx/vvvY9CgQTmee+DAAbRo0QJBQUE4fPiwQ3Z5Kmm3buGFoCDMjYmB5c8/gbt3AY0GyMwEtFogKgro1AkYPhwICJDbXIU8EAThMMm6Ob6ZmxpTmekWjNWrxQQ1hXSo5Xqo1eSAAW65hMTERO7cuZOLFi3i4MGD2axZM4aHh9PX19deFSPr8PPzY2xsLJ9++mlOnTqV69at49mzZ7MtYVitVvbo0cPuJGrZsqVDdty9e9c+niAIDA8P5+nTpx87b9WqVVSpVOzQoYPHO8zy5P598vnnSZOJVoMh78+DwSDmXe7Wjbx+XW7LFXIByvKCG0hLIytVkl50zWZSog0LzmK1Wnn48GFGRUVREARWq1aNNWvWpL+//2NLGH5+fqxevTrbtGnDKlWqEAAbNmyY6xLGw5w/f94enaHT6ejt7c2ff/6Z1x8SmawdZi+99JIrL9n17N9PBgUVPDudVituyFm3Tu4rUMiBvETXs127RQmdDtiwATAapevTbAZeew2oWFG6Pp1Ao9GgTp06OHHiBObOnYs//vgDlStXxvXr1zFkyBAsXboUFy9exNq1azFy5EhUr14dCQkJSE1NBQDYfv4ZL+n1+Fqlwlm1Gpe1Wpwzm/FTWBi2d+uG3R9+iDt37uDGjRv26Izu3bvj/PnzyMjIQFBQEL7++ms0atQI7777LjZs2IB///vfMt8VJ/j6a6BlSzHOOyWlYG2tViApCRg2DJg71zX2KbiG3NSYyky3cLz5pnRJzHv18mgHSlaayNKlS1On09HHx4f379/PflJGBvnuu0yrUIGpWi0zdbocrzcZ4H2A3wFs9GDGXLp0acbGxrJ79+4MCwuzLzeYzeYi6zCzs2+fNJ+TrM/K8uVyX5HCQ0BxpLmZpUuBqVMLPnvJwmwGnnoK+Ogj0ZHiwXz77bdo27YtSEKv1+P111/HmDFjxDf/+APo1Qs4exZITna4z1SVCidq1sSnTzyBk5cu4a+//sKZM2eynaPT6eDj4wN/f38EBwejSpUqqFGjBmJjYxETE1PglI9uJTkZCA8Hrl+Xrk+TCTh6VOxXQXbycqQpousq9u4VBSchwXHx1WoBvR5YtgwYOBAQBNfaKAEhISG4evUqMjMzAQBmsxl37tyB5ocfgK5dxWt/8F6BMBiAwEDgp5/Qqn9/fP/991CpVDCbzUhLS8PEiROhVqtx5swZXLhwATdu3EBCQgJSUlLsURgWiyVbFEb16tVRu3Zt1KtXD6VKlZLyNhSMceOA994r/JdyTqhUQHQ0cPhwkfjcFHcU0ZWL+/eBN98EXn9d/DktDXg0WbVaLc5sbTZgyBBxhly+vDz2FoIDBw7ghx9+wK5du7B3714kJyfjXw0bYv6xYxDu33euc40G9319UfGff2CpWBHjxo1DkyZNEBMTA61Wm2uzO3fu4NChQzhy5AhOnTqFc+fO4erVq7h16xbu3buHjIwMqFQqGI1G+Pr6IiAgABUrVkRERASio6NRv359hISEuGY3W2KiGO71YJ1bUsxm4NtvgYYNpe9boUAoois3mZnAnj3A/v3A7t3AtWvia2XKAI0aAfXrA61bS+uEk4lju3ahctu2sFitkvRnU6uR3LAhvH/8UbIZXHp6On777Tf8+uuvOHHiBP78809cvnwZf//9N+7evYu0tDQA4pZcb29v+xJGeHi4fQkjOjo63yWMli1bonv37hg9evT/BHzpUmDaNPFLOB9uAxgKYAeAsgDmA+ibVwNBEJ8uNm/Ot28F16KIroL76NED+PJLcVYvFWaz+Djeq5d0feZBZmYmLl26hAMHDuC3337DH3/8gYsXL+a5hBEUFISwsDBUq1YNderUQa1atRAYGAij0YioqCh88sknCAsLA554AjhwwCE7+gDIBPAegKMAOgLYByDPfGkmk5hmVFlikBVFdBXcw6lTQGystGuVWQQGijXoPCSBzZ07d3Dw4EEcPXoUJ0+exLlz53Dt2rVsSxgPIwgCmjZpgp0HD0LtwNJCMgBfACcAVH3w2gAAQQBezauhyQQcPw6EhRXmshQkIi/R9WAXr0KR4403xPjRfLA88v8UiIndl+XVKCkJ2LlTXIbxAEqVKoU2bdqgTZs2Ob7/zTffoGvXrkhNTbWvP6dduoRMqxVqB/o/A/GPs+pDr0UD2J1fQ40GOHFCEV0PxjOmDQpFHxJYtw54ZIaXE/ceOm4AMALomW+je8CqVc5a6TZSU1OhUqnQt29f7NixA6mpqfj5u++gdXDd/h6AR7NI+EBMfp8nZIHC8xTcjzLTVZCGixfFCIwCsgmAP4Amjpy8f3+B+5eLLl264O7du1CrH5rXarUOh89ZANx95LW7ALzyaygI4u5IBY9FmekqSMPRo6KoFJAPADwDwCG3z7Vrrgm1chHZBBcQ16UdjOqoCjHB/Z8PvXYM+TjRsggNdWgMBXlQRFdBGm7fLvBM9yLENcqBjjbQaMR0h0UVrRaoVMmhU80Q6+e9DNGp9hOArRCdaXmSkgLUqOGEkQquRhFdBWkoRIjSRwAaA3BMhgo/jkfRrp3DTwTLIToZ/SGGj70NB2a61auLuxoVPBZFdBWkwd9f3F1XAD5EAWa5gDiTfiiRe5FkzBiH82mUBvAZxJnuJeSzMQIAvLyAKVOcsU7BDSiiqyANdeoUaEPEPgBX4UDUwsNUrFj0nUQREeK9ckW8sVotbk5R8GgU0VWQhsBAwPJoBG7ufABxzTJfb3wWKhXQxKEYB89n9WrplwBMJtf0qyA5iugqSMfQoQ7/0a+AuKbrMEYjMHJkYazyPKpWBebMEYVSCoxGoE0boFs3afpTcCmK6CpIx7PPus7RFRQkJgYqLkycCAwY4LzwGgxAzZrixhSFIoEiugrSERIilo+ROlua0QisWFH0IxceRhCAt98Wxbew98tsFrdF//BDschQV1JQRFdBWhYuFFNWSiWQBgPQty/QvLk0/XkSgiAuM+zeLW5o8HJwhdtsBry9gfffBz7/XBHcIoYiugrSYjQCO3aIouAsBoMY6L8sz1Q4RZ969cSSRuvWic5CrVa8fxaLuPyQJbJarRj9sHSpWOqnZ4FiPxQ8BCX3goL0VKsmlitq3lzMDvZotQxHMJuBmBhg+/aSMZNTq4FOncQjLU3MFHbqlJjsXKcDKlcW74ejs2EFj0URXQXXUKMGcOYMMHw48M03jme+0mrFzQOzZwMTJhR4w0WxQK8X8xLHxsptiYILUJYXFFxH6dLApk1iJYk2bcTlAm/vxzcG6HTi62azGBb2++/A88+XTMFVKPYoM10F19O8uXjcuAHs2yeWqzl5UnyM9vISQ8Hq1hULKkoVu6qg4KEooqvgPgICgKefFg8FhRKKsrygoKCg4EYU0VVQUFBwI4roKigoKLgRRXQVFBQU3IgiugoKCgpuRCCZ+5uCcBNiKSsFBQUFBcepSNIvpzfyFF0FBQUFBWlRlhcUFBQU3IgiugoKCgpuRBFdBQUFBTeiiK6CgoKCG1FEV0FBQcGN/D/w+pDsk/6JSgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "grid.graph()"
   ]
  },
  {
   "source": [
    "\n",
    "## Setting Up Qiskit Pulse to use higher energy states"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "pharmaceutical-interval",
   "metadata": {},
   "outputs": [],
   "source": [
    "dt = backend_config.dt\n",
    "backend_defaults = backend.defaults()"
   ]
  },
  {
   "source": [
    "### Figuring out the Transition Frequencies for the transition from |0> to |1> state and from |1> to |2> state"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "scenic-brake",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Qubit 0 has an estimated frequency of 4.971857898749946 GHz.\n"
     ]
    }
   ],
   "source": [
    "GHz = 1.0e9 # Gigahertz\n",
    "MHz = 1.0e6 # Megahertz\n",
    "us = 1.0e-6 # Microseconds\n",
    "ns = 1.0e-9 # Nanoseconds\n",
    "qubit = 0 # qubit we will analyze\n",
    "default_qubit_freq = backend_defaults.qubit_freq_est[qubit] # Default qubit frequency in Hz. \n",
    "print(f\"Qubit {qubit} has an estimated frequency of {default_qubit_freq/ GHz} GHz.\")\n",
    "\n",
    "# scale data (specific to each device)\n",
    "scale_factor = 1e-14\n",
    "\n",
    "# number of shots for our experiments\n",
    "NUM_SHOTS = 1024\n",
    "\n",
    "### Collect the necessary channels\n",
    "drive_chan = pulse.DriveChannel(qubit)\n",
    "meas_chan = pulse.MeasureChannel(qubit)\n",
    "acq_chan = pulse.AcquireChannel(qubit)\n",
    "# Drive pulse parameters (us = microseconds)\n",
    "# This determines the actual width of the gaussian\n",
    "drive_sigma_us = 0.075\n",
    "\n",
    "# This is a truncating parameter, because gaussians don't have a natural finite length\n",
    "drive_samples_us = drive_sigma_us * 8\n",
    "\n",
    "def fit_function(x_values, y_values, function, init_params):\n",
    "    \"\"\"Fit a function using scipy curve_fit.\"\"\"\n",
    "    fitparams, conv = curve_fit(function, x_values, y_values, init_params)\n",
    "    y_fit = function(x_values, *fitparams)\n",
    "    \n",
    "    return fitparams, y_fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "incoming-antenna",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_job_data(job, average):\n",
    "    \"\"\"Retrieve data from a job that has already run.\n",
    "    Args:\n",
    "        job (Job): The job whose data you want.\n",
    "        average (bool): If True, gets the data assuming data is an average.\n",
    "                        If False, gets the data assuming it is for single shots.\n",
    "    Return:\n",
    "        list: List containing job result data. \n",
    "    \"\"\"\n",
    "    job_results = job.result(timeout=120) # timeout parameter set to 120 s\n",
    "    result_data = []\n",
    "    for i in range(len(job_results.results)):\n",
    "        if average: # get avg data\n",
    "            result_data.append(job_results.get_memory(i)[qubit]*scale_factor) \n",
    "        else: # get single data\n",
    "            result_data.append(job_results.get_memory(i)[:, qubit]*scale_factor)  \n",
    "    return result_data\n",
    "def get_closest_multiple_of_16(num):\n",
    "    \"\"\"Compute the nearest multiple of 16. Needed because pulse enabled devices require \n",
    "    durations which are multiples of 16 samples.\n",
    "    \"\"\"\n",
    "    return (int(num) - (int(num)%16))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "hairy-ending",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Drive pulse parameters (us = microseconds)\n",
    "drive_sigma_us = 0.075                     # This determines the actual width of the gaussian\n",
    "drive_samples_us = drive_sigma_us*8        # This is a truncating parameter, because gaussians don't have \n",
    "                                           # a natural finite length\n",
    "\n",
    "drive_sigma = get_closest_multiple_of_16(drive_sigma_us * us /dt)       # The width of the gaussian in units of dt\n",
    "drive_samples = get_closest_multiple_of_16(drive_samples_us * us /dt)   # The truncating parameter in units of dt\n",
    "# Find out which measurement map index is needed for this qubit\n",
    "meas_map_idx = None\n",
    "for i, measure_group in enumerate(backend_config.meas_map):\n",
    "    if qubit in measure_group:\n",
    "        meas_map_idx = i\n",
    "        break\n",
    "assert meas_map_idx is not None, f\"Couldn't find qubit {qubit} in the meas_map!\"\n",
    "# Get default measurement pulse from instruction schedule map\n",
    "inst_sched_map = backend_defaults.instruction_schedule_map\n",
    "measure = inst_sched_map.get('measure', qubits=backend_config.meas_map[meas_map_idx])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "eight-bradley",
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_ground_freq_sweep_program(freqs, drive_power):\n",
    "    \"\"\"Builds a program that does a freq sweep by exciting the ground state. \n",
    "    Depending on drive power this can reveal the 0->1 frequency or the 0->2 frequency. \n",
    "    Args:\n",
    "        freqs (np.ndarray(dtype=float)): Numpy array of frequencies to sweep.\n",
    "        drive_power (float) : Value of drive amplitude.\n",
    "    Raises:\n",
    "        ValueError: Raised if use more than 75 frequencies; currently, an error will be thrown on the backend \n",
    "                    if you try to do this.\n",
    "    Returns:\n",
    "        Qobj: Program for ground freq sweep experiment.\n",
    "    \"\"\"\n",
    "    if len(freqs) > 75:\n",
    "        raise ValueError(\"You can only run 75 schedules at a time.\")\n",
    "    \n",
    "    # print information on the sweep\n",
    "    print(f\"The frequency sweep will go from {freqs[0] / GHz} GHz to {freqs[-1]/ GHz} GHz \\\n",
    "using {len(freqs)} frequencies. The drive power is {drive_power}.\")\n",
    "    \n",
    "    # Define the drive pulse\n",
    "    ground_sweep_drive_pulse = pulse_lib.gaussian(duration=drive_samples,\n",
    "                                                  sigma=drive_sigma,\n",
    "                                                  amp=drive_power,\n",
    "                                                  name='ground_sweep_drive_pulse')\n",
    "    # Create the base schedule\n",
    "    schedule = pulse.Schedule(name='Frequency sweep starting from ground state.')\n",
    "    \n",
    "    schedule |= pulse.Play(ground_sweep_drive_pulse, drive_chan)\n",
    "    schedule |= measure << schedule.duration\n",
    "    \n",
    "    # define frequencies for the sweep\n",
    "    schedule_freqs = [{drive_chan: freq} for freq in freqs]\n",
    "\n",
    "    # assemble the program\n",
    "    # Note: we only require a single schedule since each does the same thing;\n",
    "    # for each schedule, the LO frequency that mixes down the drive changes\n",
    "    # this enables our frequency sweep\n",
    "    ground_freq_sweep_program = assemble(schedule,\n",
    "                                         backend=backend, \n",
    "                                         meas_level=1,\n",
    "                                         meas_return='avg',\n",
    "                                         shots=NUM_SHOTS,\n",
    "                                         schedule_los=schedule_freqs)\n",
    "    \n",
    "    return ground_freq_sweep_program"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "saved-sweden",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The frequency sweep will go from 4.951857898749947 GHz to 4.991857898749947 GHz using 75 frequencies. The drive power is 0.3.\n"
     ]
    }
   ],
   "source": [
    "num_freqs = 75\n",
    "ground_sweep_freqs = default_qubit_freq + np.linspace(-20*MHz, 20*MHz, num_freqs)\n",
    "ground_freq_sweep_program = create_ground_freq_sweep_program(ground_sweep_freqs, drive_power=0.3)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "homeless-airport",
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_ground_freq_sweep_program(freqs, drive_power):\n",
    "    \"\"\"Builds a program that does a freq sweep by exciting the ground state. \n",
    "    Depending on drive power this can reveal the 0->1 frequency or the 0->2 frequency. \n",
    "    Args:\n",
    "        freqs (np.ndarray(dtype=float)): Numpy array of frequencies to sweep.\n",
    "        drive_power (float) : Value of drive amplitude.\n",
    "    Raises:\n",
    "        ValueError: Raised if use more than 75 frequencies; currently, an error will be thrown on the backend \n",
    "                    if you try to do this.\n",
    "    Returns:\n",
    "        Qobj: Program for ground freq sweep experiment.\n",
    "    \"\"\"\n",
    "    if len(freqs) > 75:\n",
    "        raise ValueError(\"You can only run 75 schedules at a time.\")\n",
    "    \n",
    "    # print information on the sweep\n",
    "    print(f\"The frequency sweep will go from {freqs[0] / GHz} GHz to {freqs[-1]/ GHz} GHz \\\n",
    "using {len(freqs)} frequencies. The drive power is {drive_power}.\")\n",
    "    \n",
    "    # Define the drive pulse\n",
    "    ground_sweep_drive_pulse = pulse_lib.gaussian(duration=drive_samples,\n",
    "                                                  sigma=drive_sigma,\n",
    "                                                  amp=drive_power,\n",
    "                                                  name='ground_sweep_drive_pulse')\n",
    "    # Create the base schedule\n",
    "    schedule = pulse.Schedule(name='Frequency sweep starting from ground state.')\n",
    "    \n",
    "    schedule |= pulse.Play(ground_sweep_drive_pulse, drive_chan)\n",
    "    schedule |= measure << schedule.duration\n",
    "    \n",
    "    # define frequencies for the sweep\n",
    "    schedule_freqs = [{drive_chan: freq} for freq in freqs]\n",
    "\n",
    "    # assemble the program\n",
    "    # Note: we only require a single schedule since each does the same thing;\n",
    "    # for each schedule, the LO frequency that mixes down the drive changes\n",
    "    # this enables our frequency sweep\n",
    "    ground_freq_sweep_program = assemble(schedule,\n",
    "                                         backend=backend, \n",
    "                                         meas_level=1,\n",
    "                                         meas_return='avg',\n",
    "                                         shots=NUM_SHOTS,\n",
    "                                         schedule_los=schedule_freqs)\n",
    "    \n",
    "    return ground_freq_sweep_program"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dense-thesis",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We will sweep 40 MHz around the estimated frequency, with 75 frequencies\n",
    "num_freqs = 75\n",
    "ground_sweep_freqs = default_qubit_freq + np.linspace(-20*MHz, 20*MHz, num_freqs)\n",
    "ground_freq_sweep_program = create_ground_freq_sweep_program(ground_sweep_freqs, drive_power=0.3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "controlled-relationship",
   "metadata": {},
   "outputs": [],
   "source": [
    "ground_freq_sweep_job = backend.run(ground_freq_sweep_program)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "revolutionary-alloy",
   "metadata": {},
   "outputs": [],
   "source": [
    "def apply_sideband(pulse, freq):\n",
    "    \"\"\"Apply a sinusoidal sideband to this pulse at frequency freq.\n",
    "    Args:\n",
    "        pulse (SamplePulse): The pulse of interest.\n",
    "        freq (float): LO frequency for which we want to apply the sweep.\n",
    "    Return:\n",
    "        SamplePulse: Pulse with a sideband applied (oscillates at difference between freq and cal_qubit_freq).\n",
    "    \"\"\"\n",
    "    # Time goes from 0 to dt*drive_samples, sine arg of form 2*pi*f*t\n",
    "    t_samples = np.linspace(0, dt * drive_samples, drive_samples)\n",
    "    # No amplitude for the sine\n",
    "    sine_pulse = np.sin(2 * np.pi * (freq - cal_qubit_freq) * t_samples) \n",
    "\n",
    "    # Create sample pulse w/ sideband applied\n",
    "    # Note: need to make sq_pulse.samples real, multiply elementwise\n",
    "    sideband_pulse = SamplePulse(np.multiply(np.real(pulse.samples), sine_pulse), \n",
    "                                    name='sideband_pulse')\n",
    "    \n",
    "    return sideband_pulse  \n",
    "\n",
    "# More would be implemented above to create a distinguisher for the three\n",
    "# orthogonal quantum states, but we unfortunately don't have the time nor\n",
    "# resources to do a thorough implementation at the moment."
   ]
  },
  {
   "source": [
    "## Implementation of QAOA in Qiskit"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "juvenile-reality",
   "metadata": {},
   "outputs": [],
   "source": [
    "class QAOA():\n",
    "    def __init__(self, grid, p, epsilon, eta):\n",
    "        \"\"\"Quantum Appoximate Optimization Algorithm implementation for\n",
    "        Traffic Grid data structure with 3-level quantum state options.\n",
    "        p is the number of QAOA layers in the circuit.\n",
    "        Epsilon is the stepsize parameter for calculating the gradient,\n",
    "        and Eta is the stepsize parameter for jumps in gradient descent.\n",
    "        \"\"\"\n",
    "        self.grid = grid\n",
    "        self.G = grid.graph()\n",
    "        self.nodes = grid.intersections()\n",
    "        self.edges = grid.edges()\n",
    "        self.n = len(nodes)\n",
    "        self.p = p\n",
    "        self.epsilon = epsilon\n",
    "        self.eta = eta\n",
    "\n",
    "    def initialize(self, circuit):\n",
    "        \"\"\"Apply the layer of Hadamard gates to all qubits.\"\"\"\n",
    "        circuit.h(range(self.n))\n",
    "        circuit.barrier()\n",
    "\n",
    "    def hamiltonians(self, circuit, gammas, betas):\n",
    "        \"\"\"Apply the Ising type gates with learned gamma/beta, \n",
    "        for p number of layers of QAOA.# Get default measurement pulse from instruction schedule map\n",
    "inst_sched_map = backend_defaults.instruction_schedule_map\n",
    "measure = inst_sched_map.get('measure', qubits=backend_config.meas_map[meas_map_idx])\n",
    "        \"\"\"\n",
    "\n",
    "        for i in range(self.p):\n",
    "            gamma = gammas[i]\n",
    "            beta = betas[i]\n",
    "            for edge in E:\n",
    "                j = edge[0]\n",
    "                k = edge[1]\n",
    "                circuit.cp(-2 * gamma, j, k)\n",
    "                circuit.p(gamma, j)\n",
    "                circuit.p(gamma, k)\n",
    "    \n",
    "            # Apply the single qubit X rotations\n",
    "            circuit.barrier()\n",
    "            circuit.rx(2 * beta, range(self.n))\n",
    "\n",
    "        # Finally measure the result in the computational basis\n",
    "        circuit.barrier()\n",
    "        circuit.measure(range(self.n), range(self.n))\n",
    "\n",
    "    def draw(self, circuit):\n",
    "        \"\"\"Draw the circuit for comparison\"\"\"\n",
    "        circuit.draw()\n",
    "\n",
    "    def cost_function(self, x):\n",
    "        \"\"\"Cost of a certain allocation of states, specifically:\n",
    "        |0> : intersection open for North/South traffic.\n",
    "        |1> : intersection open for East/West traffic.\n",
    "        |2> : intersection red lighted, open for pedestrians.\n",
    "\n",
    "        Cost is for a single street connecting two intersections,\n",
    "        considering the increase of traffic on that street at the\n",
    "        next time interval, as well as the time people have been\n",
    "        waiting already on that street.\n",
    "        \"\"\"\n",
    "\n",
    "        E = self.grid.edges()\n",
    "        if (len(x) != len(self.intersections)):\n",
    "            return np.nan\n",
    "\n",
    "        C = 0\n",
    "        for index in E:\n",
    "            Nj = index[0]\n",
    "            j = self.intersections.index(index[0])\n",
    "            Nk = index[1]\n",
    "            k = self.intersections.index(index[1])\n",
    "\n",
    "            if Nj == Nk - self.grid.x_length:\n",
    "                direction = \"n\"\n",
    "            elif Nj == Nk + self.grid.x_length:\n",
    "                direction = \"s\"\n",
    "            elif Nj == Nk + 1:\n",
    "                direction = \"e\"\n",
    "            else:\n",
    "                direction = \"w\"\n",
    "\n",
    "            # If destination intersection is closed\n",
    "            if (direction in [\"n\", \"s\"] and x[j] != 0):\n",
    "                C += self.G[Nj][Nk][\"weight\"] * 2\n",
    "\n",
    "            elif (direction in [\"e\", \"w\"] and x[j] != 1):\n",
    "                C += self.G[Nj][Nk][\"weight\"] * 2\n",
    "\n",
    "            # If source intersection is open\n",
    "            if (direction == \"n\" and x[k] == 0):\n",
    "                # Check if there is northbound traffic\n",
    "                source = index[1] + self.grid.x_length\n",
    "                for incoming in E:\n",
    "                    if incoming[0] == index[1] and incoming[1] == source:\n",
    "                        C += self.G[Nk][source][\"weight\"] * self.grid.p_straight\n",
    "            \n",
    "            elif (direction == \"s\" and x[k] == 0):\n",
    "                # Check if there is southbound traffic\n",
    "                source = index[1] - self.grid.x_length\n",
    "                for incoming in E:\n",
    "                    if incoming[0] == index[1] and incoming[1] == source:\n",
    "                        C += self.G[Nk][source][\"weight\"] * self.grid.p_straight\n",
    "\n",
    "            elif (direction == \"e\" and x[k] == 1):\n",
    "                # Check if there is eastbound traffic\n",
    "                source = index[0] - 1\n",
    "                for incoming in E:\n",
    "                    if incoming[0] == index[1] and incoming[1] == source:\n",
    "                        C += self.G[Nk][source][\"weight\"] * self.grid.p_straight\n",
    "\n",
    "            elif (direction == \"w\" and x[k] == 1):\n",
    "                # Check if there is westbound traffic\n",
    "                source = index[0] + 1\n",
    "                for incoming in E:\n",
    "                    if incoming[0] == index[1] and incoming[1] == source:\n",
    "                        C += self.G[Nk][source][\"weight\"] * self.grid.p_straight\n",
    "\n",
    "            # If vehicles are turning onto this road\n",
    "            if (direction == \"n\" and x[k] == 1):\n",
    "                sourcel = index[1] - 1\n",
    "                sourcer = index[1] + 1\n",
    "\n",
    "                # Check if there are vehicles turning left/right into this road\n",
    "                for incoming in E:\n",
    "                    if incoming[0] == index[1]:\n",
    "                        if incoming[1] == sourcel:\n",
    "                            C += self.G[Nk][sourcel][\"weight\"] * self.grid.p_left\n",
    "                        elif incoming[1] == sourcer:\n",
    "                            C += self.G[Nk][sourcer][\"weight\"] * self.grid.p_right\n",
    "\n",
    "            elif (direction == \"s\" and x[k] == 1):\n",
    "                sourcel = index[1] + 1\n",
    "                sourcer = index[1] - 1\n",
    "\n",
    "                # Check if there are vehicles turning left/right into this road\n",
    "                for incoming in E:\n",
    "                    if incoming[0] == index[1]:\n",
    "                        if incoming[1] == sourcel:\n",
    "                            C += self.G[Nk][sourcel][\"weight\"] * self.grid.p_left\n",
    "                        elif incoming[1] == sourcer:\n",
    "                            C += self.G[Nk][sourcer][\"weight\"] * self.grid.p_right\n",
    "\n",
    "            elif (direction == \"e\" and x[k] == 0):\n",
    "                sourcel = index[1] - self.grid.x_length\n",
    "                sourcer = index[1] + self.grid.x_length\n",
    "\n",
    "                # Check if there are vehicles turning left/right into this road\n",
    "                for incoming in E:\n",
    "                    if incoming[0] == index[1]:\n",
    "                        if incoming[1] == sourcel:\n",
    "                            C += self.G[Nk][sourcel][\"weight\"] * self.grid.p_left\n",
    "                        elif incoming[1] == sourcer:\n",
    "                            C += self.G[Nk][sourcer][\"weight\"] * self.grid.p_right\n",
    "\n",
    "            elif (direction == \"w\" and x[k] == 1):\n",
    "                sourcel = index[1] + self.grid.x_length\n",
    "                sourcer = index[1] - self.grid.x_length\n",
    "\n",
    "                # Check if there are vehicles turning left/right into this road\n",
    "                for incoming in E:\n",
    "                    if incoming[0] == index[1]:\n",
    "                        if incoming[1] == sourcel:\n",
    "                            C += self.G[Nk][sourcel][\"weight\"] * self.grid.p_left\n",
    "                        elif incoming[1] == sourcer:\n",
    "                            C += self.G[Nk][sourcer][\"weight\"] * self.grid.p_right\n",
    "\n",
    "        return C\n",
    "\n",
    "    def energy_from_params(self, gammas, betas, shots=10000):\n",
    "        circuit = QuantumCircuit(self.n, self.n)\n",
    "        self.initialize(circuit)\n",
    "        self.hamiltonians(circuit, gammas, betas)\n",
    "\n",
    "        backend = Aer.get_backend(\"qasm_simulator\")\n",
    "        TQAOA = transpile(circuit, backend)\n",
    "        qobj = assemble(TQAOA, shots=shots)\n",
    "        QAOA_results = backend.run(qobj).result()\n",
    "\n",
    "        # Evaluate the data from the simulator\n",
    "        counts = QAOA_results.get_counts()\n",
    "        avr_C = 0\n",
    "        max_C = [0, 0]\n",
    "\n",
    "        for sample in list(counts.keys()):\n",
    "            # Use sampled bit string x to compute C(x)\n",
    "            x = [int(num) for num in list(sample)]\n",
    "            tmp_eng = self.cost_function(x)\n",
    "\n",
    "            # Compute the expectation value and energy distribution\n",
    "            avr_C = avr_C + counts[sample] * tmp_eng\n",
    "\n",
    "        return avr_C/shots\n",
    "\n",
    "    def gradient_energy(self, gammas, betas):\n",
    "        \"\"\"Uses a symmetric difference to calulate the gradient.\"\"\"\n",
    "\n",
    "        # Gamma-component of the gradient\n",
    "        gs = [gamma + self.epsilon for gamma in gammas]\n",
    "        bs = [beta + self.epsilon for beta in betas]\n",
    "        gs_neg = [gamma - self.epsilon for gamma in gammas]\n",
    "        bs_neg = [beta + self.epsilon for beta in betas]\n",
    "\n",
    "        grad_g = self.energy_from_params(gs, betas)\n",
    "        grad_g -= self.energy_from_params(gs_neg, betas)\n",
    "        grad_g /= 2 * self.epsilon\n",
    "\n",
    "        # Beta-compoonent of the gradient\n",
    "        grad_b = self.energy_from_params(gammas, bs)\n",
    "        grad_b -= self.energy_from_params(gammas, bs_neg)\n",
    "        grad_b /= 2 * self.epsilon\n",
    "        \n",
    "        return grad_g, grad_b\n",
    "\n",
    "    def gradient_descent(self, iterations=100):\n",
    "        \"\"\"Randomly initializes gamma and beta, learns through gradient descent.\"\"\"\n",
    "\n",
    "        angles = []\n",
    "        for i in range(self.p):\n",
    "            gamma = random.random()\n",
    "            beta = random.random()\n",
    "            print(f\"Initialized Gamma {i + 1}: {gamma}\\nInitialized Beta {i + 1}: {beta}\\n\")\n",
    "            angles.append(gamma)\n",
    "            angles.append(beta)\n",
    "\n",
    "        for i in range(iterations):\n",
    "            grad_g, grad_b = self.gradient_energy(gammas, betas)\n",
    "            gammas = [g - self.eta * grad_g for g in gammas]\n",
    "            betas = [b - self.eta * grad_b for b in betas]\n",
    "            if not i % 25:\n",
    "                energy = self.energy_from_params(gammas, betas)\n",
    "                print(f\"Step {i} Energy: {energy}\")\n",
    "\n",
    "        for i in range(self.p):\n",
    "            print(f\"\\nLearned Gamma {i + 1}: {gammas[i]}\\nLearned Beta {i + 1}: {betas[i]}\\n\")\n",
    "\n",
    "        return gammas, betas\n",
    "\n",
    "    def simulate(self, shots=10000):\n",
    "        \"\"\"Run QAOA circuit on local simulator\"\"\"\n",
    "        circuit = QuantumCircuit(self.n, self.n)\n",
    "        self.initialize(circuit)\n",
    "        gammas, betas = self.gradient_descent()\n",
    "        self.hamiltonians(circuit, gammas, betas)\n",
    "\n",
    "        backend = Aer.get_backend(\"qasm_simulator\")\n",
    "        TQAOA = transpile(circuit, backend)\n",
    "        qobj = assemble(TQAOA, shots=shots)\n",
    "        QAOA_results = backend.run(qobj).result()\n",
    "\n",
    "        plot_histogram(QAOA_results.get_counts(), figsize = (8, 6), bar_labels = False)\n",
    "\n",
    "        # Evaluate the data from the simulator\n",
    "        counts = QAOA_results.get_counts()\n",
    "        avr_C = 0\n",
    "        max_C = [0, 0]\n",
    "        hist = {}\n",
    "\n",
    "        for k in range(len(grid.edges()) + 1):\n",
    "            hist[str(k)] = hist.get(str(k), 0)\n",
    "\n",
    "        for sample in list(counts.keys()):\n",
    "            # Use sampled bit string x to compute C(x)\n",
    "            x = [int(num) for num in list(sample)]\n",
    "            tmp_eng = cost_function(x, G)\n",
    "\n",
    "            # Compute the expectation value and energy distribution\n",
    "            avr_C = avr_C + counts[sample] * tmp_eng\n",
    "            hist[str(round(tmp_eng))] = hist.get(str(round(tmp_eng)), 0) + counts[sample]\n",
    "\n",
    "            # Save best bit string\n",
    "            if (max_C[1] < tmp_eng):\n",
    "                max_C[0] = sample\n",
    "                max_C[1] = tmp_eng\n",
    "\n",
    "        M1_sampled = avr_C/shots\n",
    "\n",
    "        print(\"\\n --- SIMULATION RESULTS --- \\n\")\n",
    "        print(f\"Sampled mean value = {M1_sampled:.02f}\\n\")\n",
    "        print(f\"The approximate solution is x* = {max_C[0]} with C(x*) = {max_C[1]}\\n\")\n",
    "        print(\"The cost function is distributed as: \\n\")\n",
    "        plot_histogram(hist, figsize = (8,6), bar_labels = False)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "name": "python385jvsc74a57bd0bb76021927b5b2fe67fa472802372d33b40b9e48188a8eba32b6a9e8190a9a16",
   "display_name": "Python 3.8.5 64-bit ('base': conda)"
  },
  "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}