Egan_single_qubit_errors.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Single qubit errors validated against Egan paper\n",
    "========================================\n",
    "\n",
    "This is a simulation of single qubit errors in an ion trap quantum computer to compare the Qiskit noise model with experiments by Egan [1].  The noise model will be used on a simulation of the Steane code.\n",
    "\n",
    "This simulation uses the randomised benchmark function from Qiskit to simulate Egan's characterization of single qubit errors.  A circuit is constructed with N-1 sequential random Clifford gates and a final Clifford gate that inverts the others.  In the circuit diagram below the Clifford gates are separated by a barrier.\n",
    "\n",
    "The circuit is transpiled into the native ion computer gate set before computing the noise.\n",
    "\n",
    "With no noise the qubit will always finish in the starting |0> ground state.  A simulation was run and the ground state population is shown plotted against the number of Clifford gates.  For comparison a calculation of the expected ground state population using a error per Clifford gate of 3.4(8)×10−4 as quoted by Egan is plotted.\n",
    "\n",
    "We use a depolarising error with a parameter which represents the chance that a X, Y, Z and I error occurs.  Since a Z and I error won't be detectable we have doubled Egan's stated error of 1.8(3)×10−4 per native Pauli gate to use 3.66 1.8(3)×10−4 as the depolarising error parameter.\n",
    "\n",
    "Egan shows two results for the fidelity per Clifford gate, an upper bound for the 1 ion chain, and a lower bound for the less accurate but more realistic 15 ion chain.  Egan's upper bound for the 1 ion chain is 0.990(3) and the lower bound from the 15 ion chain is 0.988(0).  The simulation shows 0.9877, slightly below Egan's lower band.  It should be noted that the simulations are random, and some noise is inevitable.  The single gate error is not a key parameter.  It is an order of magnitude less that the two gate error and is often neglected in simulations by other workers. \n",
    "\n",
    "The intercept on the y axis in our simulation is 0.995.  This is in agreement with Egan's results and shows that 0.046 is a reasonable parameter for the measure noise.\n",
    "\n",
    "[1] Egan L, Debroy DM, Noel C, Risinger A, Zhu D, Biswas D, et al. Fault-Tolerant Operation of a Quantum Error-Correction Code. arXiv Available from: http://arxiv.org/abs/2009.11482\n",
    "\n",
    "This is part of an MSc project on the Steane code.  Please see https://goldsmdn.github.io/ for full details."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import the necessary function modules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit import(\n",
    "    QuantumCircuit,\n",
    "    QuantumRegister,\n",
    "    ClassicalRegister,\n",
    "    execute,\n",
    "    Aer\n",
    "    )\n",
    "\n",
    "# Import the randomized benchmarking Function and other qiskits modules\n",
    "import qiskit.ignis.verification.randomized_benchmarking as rb\n",
    "from qiskit.providers.aer.noise import NoiseModel\n",
    "from qiskit.providers.aer.noise.errors import pauli_error, depolarizing_error\n",
    "from qiskit.compiler import transpile\n",
    "\n",
    "import random\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "from helper_functions import (\n",
    "    get_noise, \n",
    "    mean_of_list,\n",
    "    calculate_standard_error\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define constants so the process flow can be controlled from one place:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "SINGLE_GATE_SET = ['id', 'ry', 'rx']\n",
    "TWO_GATE_SET = ['rxx']\n",
    "BASIS_GATE_SET = SINGLE_GATE_SET + TWO_GATE_SET\n",
    "\n",
    "SHOTS = 10000                            #Number of shots to run \n",
    "ITERATIONS = 25                          #Number of iterations for each simulation\n",
    "MEASURE_NOISE = 0.0046                   #Measurement noise\n",
    "SINGLE_GATE_DEPOLARISING = 0.000366      #Parameter for single gate depolarising channel\n",
    "TWO_GATE_DEPOLARISING = 0.0              #Two gate noise not relevant\n",
    "SINGLE_CLIFFORD_GATE_NOISE = 0.000348    #Used for calculation\n",
    "NOISE = True                             #Test with noise?\n",
    "\n",
    "SIMULATION_GATE_DEPTH_LIST = [1, 3, 7, 11, 15, 19]  #Number of Clifford gate for simulation\n",
    "CALCULATION_GATE_DEPTH_LIST = [0, 4, 6, 8, 12, 16, 20]  #Number of Clifford gate for calculations\n",
    "\n",
    "SIMULATOR = Aer.get_backend('qasm_simulator')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set up circuits in a list."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate randomised benchmark circuit \n",
    "rb_opts = {}\n",
    "#Number of Cliffords in the sequence\n",
    "rb_opts['length_vector'] = SIMULATION_GATE_DEPTH_LIST\n",
    "# Number of seeds (random sequences)\n",
    "rb_opts['nseeds'] = ITERATIONS\n",
    "rb_circs, xdata = rb.randomized_benchmarking_seq(**rb_opts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Draw an example of random circuit for eight Clifford gates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAyMAAAB7CAYAAABn5l/CAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAARFklEQVR4nO3dfVRV9Z7H8c95SFQ0Q620RKz0qhDqMnq4tAK8eq2xogTGh7nLriup1PSqi2U5Y4qRWWo+XFOrW63RWa1ZhVKrSZu6VoKg1Q1GrDRKM0rN8onJRCWFM3+ITITC4XDO+e29eb/+Ou7HL7/1/bnX5+yzz3H5fD6fAAAAACDM3KYLAAAAANA6EUYAAAAAGEEYAQAAAGAEYQQAAACAEYQRAAAAAEYQRgAAAAAYQRgBAAAAYARhBAAAAIARhBEAAAAARhBGAAAAABhBGAEAAABgBGEEAAAAgBGEEQAAAABGEEYAAAAAGEEYAQAAAGAEYQQAAACAEYQRAAAAAEYQRgAAAAAYQRgBAAAAYITXdAHAxZSVlTW6fuXKlZoyZUqj2/Tr1y+YJQGtUlNzUWI+wn6C0df0NNBy3BmBba1atcp0CQBqMR/hRPQ1EHqEEQAAAABGEEYAAAAAGEEYgW2tX7/edAkAajEf4UT0NRB6hBEAAAAARhBGYFsZGRmmSwBQi/kIJ6KvgdDjq30v4PVi6UBF+M97dZSUlhDYvqZqllpWtx1lfblLO37+2ci5B3bsqCV9YwPal75unkDrtmPNdmXXsWYuNg99HT527JGW1GzX6/n06dNVWloa3IL8MGjQIC1fvjzoxyWMXMCBCunrQ6araB471mxXO37+WVsqjpkuo9ns2CPUjMbYdaztWLcda7Yru461Heu26/W8tLRUBQUFpssIGj6mBdt6+OGHTZcAoBbzEU5EXwOhRxiBbTX1a88Awof5CCeir4HQI4zAtpKSkkyXAKAW8xFORF8DoUcYgW0dPnzYdAkAajEf4UT0NRB6hBEAAAAARhBGYFuxsYF9JR6A4GM+wonoayD0CCOwrby8PNMlAKjFfIQT0ddorTp16hS2cxFGYFtz5841XQKAWsxHOJET+vr44XL9/YXxpsuAIX379tXMmTP16quvqqSkRJ9++qmKioq0atUq3XffferQoUODfeLj4/XVV1/poYceCkuNhBHY1rp160yXAKAW8xFORF/Drm6++Wa99957Kisr06JFizR69GgNHjxY8fHxuvXWWzV58mStXbtW33//vZYvX153JyQ+Pl7vv/++rrjiCt11111yuVwhr9Uxv8B+9uxZuVwueTyesJ97/fwU9bx+mG669zG/llvJT4f2qujVR/X9l4U6c/qEIiKjdOU1Cfqnqa/J421jujxHqPnsc1XPvsC7a9XV0pkz8ixZJHf89eEvzA926w/mIvxht7GmrxGo/LVT9ePeT3T8SLnWz0/RbX9aoiuvucF0WZLs2dd2uJ57vV49+eSTysrKksfjUWVlpXJzc7Vlyxbt3LlTVVVV6tKliwYPHqzU1FQlJSVp2rRpysjIUE5OjubPn6/LL79cGzduVEZGhnw+X+hrDvkZAvDGG29ozpw52r17t2JiYvTII4+oqKhI+fn5Ki8vV3l5ua655hotW7ZMR44cqUt2ZWVl6tOnj+nybeXNxSPUM3647lv8pdq0u1Qnjh3QN6UbwtJ8rYU7/nq5/+v1est8v/yi6qxHpMsukyvOug9I0h/hw1iHD2MdPox16PgT9FL+/KyOHy7XR6/P0/CH1pgt2AGsfj33er3Kzc3VyJEjVV1drUWLFmnBggX66aefGmy7efNmLVmyRAMHDtTq1auVmJio559/Xi6XSxs3blR6erqqqqrCU3dYztIM77zzjtLT0zVs2DAtWLBAlZWVys7O1qlTpxrc9Vi8eLEGDRqklStXyuv1qkuXLoaqtqdTPx9VxcEvdef01xXR/tztuY5demjA0ImGK/NPQUGB6RICVr1kmXy//CLvvz0ql9uan5a0e3/YiRPG2i7z0QljbRdOGGsr9zVBzxqsdD1fsWKFRo4cqWPHjik1NVVbt25tcp8dO3Zo8uTJ2rZtm9q3b6+amhotXbo0bEFEsuAzI3PnzlWvXr309ttvKzU1VWPHjtWmTZv0448/Ntg2KipKb731lu655x7deeed6ty5s4GK7atdxy7q0iNO77+UqS8K/0NHD+yy1X9iO3fuNF1CQKpf+U/5tu+QNydbrnbtTJdzUXbvDztxwljbZT46YaztwgljbdW+Ph/04odOVET7TnK5XHVBz3tJhOnyWg0rXc+HDx+uSZMmqaqqSrfffrtfQUQ694zIpk2b1L59e+3du1dut1svvPCC2oXx77HUnZHKykoVFxcrKytLXu//lxYTE6PExESVl5fX2z41NVXuZqRQfx/CSZ+9WT36p/h9XEn6x5tPquTtZ+otO3P6hHpeP8zvYxQU5Osvw4c067znBVLzuf3y9T9vL9X2d5br6P7PFRF5mQb+capuuvcxv8erJXU3ZsaMGY2uX7ZsmV/bBJtn8dNyDxwQ0L41WwpV89o6eRYukOvKK5u9f35+vlw33hLQuQPpEdP9YWouSoHX3RrnomRmPtp1rLnGOKevQ3GN8Wesfx304v8wUVdcm6DOV/Vv8cPGra2v7Xo9/y2Xy6UVK1ZIOvemfnFxsV/7nX9Y/fwzImPGjNHWrVs1YMAATZ8+XU899VS97QsKCprVY/6++WCpMFJRUSGfz6du3bo1WNetW7cGYaR79+5hqqxpN90z+4IPYVldu45ddevoBbp19AKdqTqp3R/n6v2XH1CHzlcrLvl+0+U5Ss1XX6l68VJ5pv9F7tj+psvxix37g7mIpthxrOlr/Ja/Qc9zSVtFde9nsNKLs2tfW+16PnToUPXt21f79u3TkiVL/Nrnt0Hk/DMiM2fO1LvvvquJEydq4cKFqqmpCXH1FgsjUVFRcrlc+uGHHxqsu9Cy5r4D4G9Ce3aT9PWhZh06KJKTU7R+fmC3sINR8yUR7RWbNF47/v6sDn9b6vd+Lam7MWVlZY2uX7ZsmR588MFGt1m6dGkwS5IkDSv+WFsqjjVrH9+RI6rOfkLu9JFyDw38Hb6UlBS9F+DHHFraIyb6w9RclAKvuzXORcnMfLTrWHONcU5fh+Ia4+9Y+xv0Ii/rphtTZ/l1btM9EoiW1GzX63lKSkq955nGjh0rSXrxxRdVXV3d5P4XCyKStGnTJu3Zs0e9e/dWYmKiioqK6vZLTk5Wfn5+QDU3xlLPjERGRiohIUF5eXk6e/Zs3fJvv/1W27ZtM1iZM52urNDW1/5VR/Z9ruqzZ1RTfVa7/5Gno/s/19V9bzNdXpMef/xx0yX4xXf6tKqzc+SK7S/3n8eZLsdvdu8PO3HCWNtlPjphrO3CCWNtl74+H/S6Rg9oVtBD81j1ep6QkCDpXJBoSmNBRDr3xv0HH3xQ77ihZqk7I5KUk5OjESNGaMSIEZo6dapOnDihefPm6coAPo+Hxnk8bXTy+CFt/GuaKv/3oNxury69vJeSx61Qn5v/2XR5TRo1apTpEvziK9oq3+490nf7dPae9AbrPdOmtujdlVCxe3/YiRPG2i7z0QljbRdOGGur9vXpygqVbFikvol/UlT3vnK5XPq65E0d3f+5Eu727y4Ims+q1/PY2HNfKbxjx45Gt2sqiJxXWloqSYqLiwt6rRdiuTByxx13aP369ZozZ47S09MVExOjWbNmqbCwMCS3hoIh47H8Zi23ikvaRuqPD7xsuoyA9e/fX1988YXpMprkHjZU7mFDTZfRbHbsD+aiOXaZj3Yca/raHKv2tROCnh372orXc7fbrZycHLVp00anTp1qdLvc3Nwmg4gkffTRR5o/f75KSkpCVXY9lgsjkpSWlqa0tLR6ywoLC+te9+rVy3ZfDwgAABAMTgh6CI6amho98cQTfm03duxYZWVlKTMzs9HfEdm+fbu2b98ezDIbZckwAgAAACB4SktLNW6cdZ51Oc9SD7ADzZGSkmK6BAC1mI9wIvoaCD3b3BlZs2aN6RJgMc8995zpEgDUYj7CiehrIPS4MwLbmjRpkukSANRiPsKJ6Gsg9AgjsC2rfrsa0BoxH+FE9DUQeoQRAAAAAEYQRgAAAAAYQRiBbVnxh6iA1or5CCeir4HQI4zAtnJzc02XAKAW8xFORF8DoWebr/YNp6uj7HdeUzWbPHd2drZGjRoV9vMO7Ngx7OcMxrnp6/Cc2441B4OJ+WjXsWYu2ufc9HV49m2JlpzXrtfzQYMGNXufvd8dlCRd27N7vdehPq8/XD6fzxeSIwMtVFZW1uj6/v37N3kLvV+/fsEsCWiVmpqLEvMR9hOMvqanYRezFv5NkvT0ow/We20FfEwLAAAAgBGEEdjW6tWrTZcAoBbzEU5EXwOhRxiBbcXFxZkuAUAt5iOciL4GQo8wAttKTk42XQKAWsxHOBF9DYQeYQQAAACAEYQRAAAAAEbwOyOwrKa+MjE7O5uvVQTCwJ95xnyE3dDXgDVwZwS2NW/ePNMlAKjFfIQT0ddA6BFGAAAAABhBGAEAAABgBGEEAAAAgBGEEQAAAABGEEYAAAAAGEEYAQAAAGAEYQQAAADABeXn5ysuLk69e/dWZmamqqurg3p8wggAAACABmpqapSZmal169Zpz549On78uF555ZWgnoMwAgAAAKCBTz75RFdddZViY2MlSRMmTFBeXl5Qz0EYAQAAANDA/v37FR0dXffvnj17at++fUE9hzeoRwMAAABg1Nfffq8NH3zYYPlf/z2vwevIdm11X/rtanNJw1jg8/lCV2Qt7owAAAAADnJtz+6K6tRBBw8d1cFDR+uW//b1wUNH9fvBcRcMIpIUHR1d707Id999px49egS1VsIIAAAA4CAul0tptyepQ/t2jW53Q/zvFPe7Xhddn5CQoP3792vXrl2SpJdffllpaWnBLJUwAgAAADhNh8h2Srvjtouuj+rUUXcPTWz0GB6PRy+99JIyMjJ03XXXqUOHDho3blxQ63T5wvFhMAAAAABht/6/C1T86Zf1lrkkPfAvd+va6O5mivoV7owAAAAADnX3H36vzp061lt2200DLBFEJAeFkaqqKtMlAAAAAJYSEdFGo+4aIlftv7td3lnDb7vRaE2/Ztkw8tlnnyk9PV1du3ZV27Zt1adPH82ePVuSNH78eHXt2lXFxcVKTk5WZGSk7r//fsMVAwAAANbTq0c3Jd08UB6PW6PvGiKv12O6pDqW/J2RkpISJSUlKTo6WosXL1ZMTIy++eYbbdu2rW6bkydPKiMjQ9OmTdO8efPk9Tb9p8xa+LdQlg0AAABY2q9/aySUnn70Qb+2s2QYycrKUmRkpD7++GN16tSpbvmECRPqXp86dUoLFy7U6NGjTZQIAAAAoIUsF0ZOnjypwsJCTZkypV4Q+S2Xy6V77723Wcf2N6EBAAAACD3LhZGKigrV1NQ0+euOUVFRioiIaNax+ZgWAAAAEHr+3gSw3APsUVFRcrvdOnDgQKPbuVyuRtcDAAAAsDZL/ujhkCFDtGvXLu3evVuXXnppg/Xjx4/Xhg0bdOTIEQPVAQAAAAgGy90ZkaRnnnlGJ06c0C233KI1a9Zo8+bNWrt2rTIzM02XBgAAACBILPfMiCTdcMMN+vDDDzVnzhzNmDFDp0+fVnR0tMaMGWO6NAAAAABBYsmPaQEAAABwPkt+TAsAAACA8xFGAAAAABhBGAEAAABgBGEEAAAAgBGEEQAAAABGEEYAAAAAGEEYAQAAAGAEYQQAAACAEYQRAAAAAEYQRgAAAAAYQRgBAAAAYARhBAAAAIARhBEAAAAARhBGAAAAABhBGAEAAABgBGEEAAAAgBGEEQAAAABGEEYAAAAAGEEYAQAAAGAEYQQAAACAEYQRAAAAAEYQRgAAAAAYQRgBAAAAYARhBAAAAIARhBEAAAAARhBGAAAAABhBGAEAAABgBGEEAAAAgBH/B+0bhQJFbtNeAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1045.18x144.48 with 1 Axes>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rb_circs[0][1].draw(output = 'mpl', filename = 'Single_Clifford_gates_circuit.jpg', fold = 60)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Transpile the circuit and draw the equivalent transpiled circuit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "qt = transpile(rb_circs[0][1], basis_gates = BASIS_GATE_SET)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm0AAACFCAYAAAAJg894AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAZp0lEQVR4nO3de1RVdf7/8SegoqAohqOigKhpCCqKonnFKDLrZ4iYeZmRRmpSHK3RqWY0zSYrs28yJtqkefmO98tPzb5WP8N0Wajf1DAvWN7wguYlLyRgxeX3hwtGlMtBD2efj74ea7nWYZ99eZ1P77P3u33O2duloKCgABERERFxaq5WBxARERGR8qlpExERETGAmjYRERERA6hpExERETGAmjYRERERA6hpExERETGAmjYRERERA6hpExERETGATU3b2bNnGTx4ME2bNiUsLIwHH3yQNWvWAFCzZs0yl01PTyckJKTCwUpbr5ubG6GhoYSEhDBgwACys7Nvexv2MH36dIKDgwkJCWHQoEFcu3bNkhwiIiJydyu3aSsoKCA6OpoePXpw9OhRdu3axbJlyzh16pQj8t2iRo0apKamsm/fPqpVq8YHH3xgSQ6AjIwMZsyYwc6dO9m3bx95eXksW7bMsjwiIiJy9yq3adu0aRPVqlXj+eefL5oWEBDAn//851vmfe+99wgJCSEkJITExMSi6bm5uQwZMoSgoCBiY2PJzs4GIDo6mrCwMIKDg/nwww8rHL579+4cPnwYgLy8PJ599lmCg4OJiooiJyen1G1kZWXx+OOP07ZtW0JCQli+fDkAixYtIjw8nNDQUP70pz+Rl5dXbobc3FxycnLIzc0lOzsbX19fjhw5Qr169WjSpAmhoaHUrVuXZs2akZmZWeHXKCIiIgI2NG379++nffv25a5o165dzJ8/nx07drB9+3bmzJnDt99+C8D333/PyJEjSUtLw8vLi1mzZgEwb948du3axc6dO5kxYwY//fSTzcFzc3P59NNPad26NQCHDh0iISGB/fv3U6dOHVavXl3qNj777DN8fX3Zs2cP+/bto3fv3qSlpbF8+XK+/vprUlNTcXNzY/HixQD06dOH06dP35KhUaNGjBs3Dn9/fxo2bEjt2rWJioqiWbNmdOvWjX//+9+kpqbSpk0b1q5di5eXl82vT0RERORGFf4hQkJCAm3btqVjx47Fpn/11Vf069cPT09PatasSUxMDFu3bgXAz8+Prl27AjB06FC++uorAGbMmEHbtm3p3LkzJ0+e5NChQ+VuPycnh9DQUDp06IC/vz/Dhw8HIDAwkNDQUADCwsJIT08vdRutW7dm48aNvPzyy2zdupXatWuTnJzMrl276NixI6GhoSQnJ3P06FEANmzYgK+v7y1ZLl26xLp16zh27BinT58mKyuLRYsWAdeb3cLv2aWlpdGyZcuKDLOIiIhIMVXKmyE4OLjorBVAUlISFy5coEOHDjZvxMXF5Za/N2/ezBdffMG2bdvw8PAgIiLCpi/xF36n7Wbu7u5Fj93c3MjJySl1Gy1atGD37t1s2LCBCRMmEBkZibe3N8OGDeOtt96y+XV98cUXBAYGUq9ePQBiYmJISUmhf//+XLt2DW9vb06ePImPjw/VqlWzeb0iIiIiNyv3TNtDDz3EtWvXmD17dtG0wu+k3ah79+6sXbuW7OxssrKyWLNmDd27dwfgxIkTbNu2DYAlS5bQrVs3rly5gre3Nx4eHhw8eJDt27fb6zUVKW0bp0+fxsPDg6FDh/LXv/6V3bt3ExkZyapVqzh37hwAFy9e5Pjx42Wu39/fn+3bt5OdnU1BQQHJyckEBQVx4MABgoKCgOtn2Qofi4iIiNyucps2FxcX1q5dy5YtWwgMDCQ8PJxhw4YxderUYvO1b9+euLg4wsPD6dSpE/Hx8bRr1w6Ali1bkpSURFBQEJcuXWLEiBH07t2b3NxcgoKCeOWVV+jcubPdX1xp29i7d2/RDw4mT57MhAkTaNWqFW+88QZRUVG0adOGRx55hDNnzgClf6etU6dOxMbG0r59e1q3bk1+fj7PPfdcsY9Ga9Sowe7duzl48KDdX5+IiIjcO1wKCgoKrA4hIiIiImXTHRFEREREDKCmTURERMQAatpEREREDKCmTURERMQAatpEREREDKCmTURERMQAatpEREREDKCmTURERMQAatpEREREDKCmTURERMQAatpEREREDKCmTURERMQAVawOIFKagwcPlvn8zJkzGTVqVJnzPPDAA/aMJHJHyqtpUF3bi8Za7kY60ybGSkpKsjqCiN2prh1HYy2mUdMmIiIiYgA1bSIiIiIGUNMmxlq1apXVEUTsTnXtOBprMY2aNhEREREDqGkTY8XGxlodQcTuVNeOo7EW0+iSHyX4fhP8fM7x2631O2j50O0ta1VmuLPcJjJ1rFXXFaO6dox7bZzBzLE29b34wgsvkJqaatc8tggNDSUxMdHu61XTVoKfz8HlU1anqBgTM5vK1LE2MbeJmU2lsXYcE8faxMwAqampbNmyxeoYdqOPR8VYCQkJVkcQsTvVteNorMU0atrEWOVdyVzERKprx9FYi2nUtImxevToYXUEEbtTXTuOxlpMo6ZNjHX+/HmrI4jYneracTTWYhr9EOE2jZ0dQdrxbbi5VcXV1Y0G3oEMjhxPz7YDrI5WJlNzm8bUcTYxt4mZTaWxdhwTx9rEzKZR03YHhjz8KkMenkBeXi7rUmby1pLBNG/UjkY+za2OViZTc9+sVatWVkcok6njbGJuEzOXRnXtOBpr+zMxs0n08agduLlV4bFOz5KXn8uR06lWx7GZqbkLrV692uoINjF1nE3MbWLmm6muHUdjXXlMzHy7vLy8cHFxcci21LTZwW+5v/JJymwAGvu0sDiN7UzNXWjixIlWR7CJqeNsYm4TM99Mde04GuvKY2Lm+vXrM3LkSObNm8eOHTv47rvv2LFjB/PnzychIYEGDRrcsoyPjw9bt24lKSnJIY2bmrY7sCR5CtGv1uGJv9dg/ucT+MuAuTT1bQNAxoXDjEwM47fcXwFYsXkaCz53jh1EWbnfXDyY7Qc+KZp30oJodn7//6yKWqaVK1daHaFMqg/HMXWsS6K6dhyTx1rvRfsJCAhgyZIlnDx5kqSkJJ555hnCw8Np3bo14eHhxMXFMXPmTE6cOMGyZcsIDAwErjdsycnJtGnThoiICGrXrl3pWe+api03N5e8vDyHbnNw5HjW/uMyq167QPgDfdhz+Mui5xr5NKdb6/4s2/QWZy4eY3PqMgZHjndovtKUlXvEk4ks+PxVcn65yta9/xfP6rXp0DLKwrTmUn04jqljbSKNtePovVj54uPj2bdvH4MGDcLV1ZW1a9eSkJBAt27daNOmDV27dmXkyJGsWbMGFxcXBg4cyN69exk3blxRw5aWlkavXr24fPlyped1yqZtzZo1hISE4O7uTosWLZg7dy5xcXE0adIEgPT0dFxcXEhMTGTChAn4+fnh7u7O0aNHLclby8ObvwyYy46D/0PKvnVF05+K+Cvb0z7hzcWDGNE3kWpV3C3JV5qScnvX/B39uo0had1oliS/wfN9p1uc0nyqD8cxdaxNpLF2HL0XK8fkyZOZM2cONWvWZMWKFQQGBtKvXz9mzZrF119/zd69e0lJSWH27NnExMTQpEkTli5diqenJ9OmTSvWsJ09e9YhmZ2uafvss8/o378/vr6+rFy5ksmTJ/POO++QnJx8y7zTpk3j22+/ZebMmXz88cfcd999FiS+zsujLv27/4V5n/2d/Px8AKq4VaV10x5czblESGA3y7KVpaTcj3aM49T5H4juOhovj7oWJyydSfeTU304jqljXUh17Timj7Xei7cvPj6eiRMnkpubS1xcHAMHDuTkyZNlLpORkcHo0aM5ceJE0bSFCxc6rGEDJ2zaJk6cSJMmTdiwYQN9+/Zl0KBBbNy4scRB8fb2Zv369Tz55JM8/vjj1K1rbdH26z6Gi5ln2LjrvwFI/3E/+9O/pl3zh9mwY46l2cpyc24A3/uaO/1PtPfv3291hApRfTiOqWMNqmtHMn2sQe/F2xEQEMD06dfPTMbHx7Nw4UKbliv8Dpu/vz8ZGRkAvPrqqzRt2rTSst7MpaCgoMBhWytHVlYWtWrVYuzYsUybNq3YcxEREaSnpxf9CwwM5G9/+xtvvvmmzeu39Zcd7z7/JW2bRVQk+i3y8/MZ+0FPRvRNpLFPC8YkdeGd577Au1b9UpfZc2Qz4z7odVvbs0fmG72zLI4+neJt+j+hO8ldlhdffLHM56dPn27TPPZmVX2A89SIo+rD1LEuTXn1CvdeXVu1/wCzx7qQ3osVt2TJEgYNGsSKFSsYOHCgTcvc+KODwo9E33vvPQYPHsyqVasYMODOLiBsayvmVGfaLl26REFBQYk/qy1pWsOGDR0R67as3zab+xuF0aJxGB7VaxH36D+Y9fELVscSJ6H6cByNteNorKUszlAf9evXJzY2lry8PMaOHWvTMiU1bGfPnmXcuHH89ttvREdH4+vrW8nJrzP6TNv777/PqFGj7J5j5zK4fMruqy1XncbQ4enbW9aqzHBnucty8ODBMp8PCgoiLS2tzHkeeOABe0YCzB1r1XXFVEZdl1fTcO/VtVX7D9BYV4Sp78WIiIhi310cOXIkSUlJrF27ln79+pW7fGkNW6GVK1cSGxvLmDFjmDFjRtH0nj17snnz5tsLXQanOtPm6elJhw4dWL16Nbm5uUXTjx8/TkpKioXJxBlNnjzZ6ggidqe6dhyN9b2nQ4cOAGzcuLHcectr2G5cT+F6K5tTNW0Ar7/+Ounp6fTp04f169ezdOlSoqKiqF+/7M+85d7z1FNPWR1BxO5U146jsb73BAcHA7Bnz54y57OlYQNITU0ttt7K5nRNW+/evVm1ahUZGRn079+fiRMn8tJLLxEZGWl1NHEyQUFBVkcQsTvVteNorO898+bNY8qUKRw5cqTM+WbOnGnTddiOHTvGlClTmDt3bmXEvUUVh2ylgmJiYoiJiSk2bevWrUWPmzRpYvMvLUQcLe3EDj74+EVcXFxp6deRETdc9PLCldNMXTqUX3OvMSzqddq3eBi4/suh56e348muo+jTKd6q6MY59uM+Elc9h6urG773NWfcU/OKfiWemX2Rf65+nsysC4TeH8mQyPEsTp7CxylJ9O74R57p/YbF6UXE0f71r3/ZNN/o0aOpUqUKCQkJZV6H7fz580yYMMFe8crllE2bMyvrIHGjC1cyWLF5Gr3aDSr1AO6sLlzJ4L2V8WRdu2JUbmdRv04A0/60iWpVq/PWkiEcO7OXwIatAVj+5dsMe/QfNPNty4R5TxQ1bdsOrKdOzXpWxjaSX72W/HPU9e+7Tlv+DD+c2klLv44A/HvjZIY9+jr+v/vPF8n7hMcTHNCFbw/ferFuKd29sN8zkTPuq22tFXDuejl37hyxsbFWx7iF03086uwKDxLTR14/8/fDqZ0lzrfrh42EtXik6ACemPAVl6+e49iZvY6Me1t2/bCRyPZDnT53RESE1RFKVNerAdWqVgfAzbUqrq5uRc8d+3EvwU26UMO9Jh7utci6lgnAl98uISK0En5CZ0d7jmym30Rvxs6OYMiUACbOf9LqSFRxq1r0uGoVd+rV9iv6O/3HfSxNfpNxH/TiQPo2ALxr1bf5eo1Wcca6vlv3e8441hXhjPtqW2sFzKsXZ2DMmbYFCxZYHQEo+SCx58hmXlvYj6YN2/LjxWM08w2llkddRkW/Tw33mkXz33wAt1pJuV9/Zh3fHd3CqOj3S208nMXs2bOtjlCmo6e/40rWeQLqtyqalp+fV9Q0eFavTVbOZdKOb6dN0564uriRV5Bb2uos1zqwBy39wnn72c+ZuvQPDO/zltWRAEjZ/zHzP/07jXzux8vzP7eyO5CewqwXduPlUZfJ/92fxISvLExpO2es67tpv3cjZxzrkpi0r7a1Vm7Mb0q9OAOdabsNKfs/5tl3Q7j881m8PO8rOpj914jNtGnak9Exs7j2a1axQizpAG61knIXFBQUy+6MuQuNGDHC6gilysy+yMy1oxg74KNi011c/vOWy/olE88adfj0f+fyaMdnHB2xVBczf2Ts7Ihi/6YsepozF4/SsO7127Wcv3ISn9qNLE56XZfgvswZtw+fOo3ZfuCToumN67UgoH4Q3rXq4+pizq7OWev6btnv3chZx/pmpu2rbamVm/ODc70GZ2XMmTZn0iW4L12C+zJz7Z/ZfuATAhu2LnYwu3z1PE192xbNX3gAnzB0hSV5L2b+yJTFxT96q1urAXG937jlIHw4I7Uou9W5y1MZFy60h7y8XN5eOpTnnniXul7F7+TRtGEbDqRvI9C3DdnXMvGs7sWp8z8waUE0FzIzoKCAkMBuxb6H5Wh1vRrwXyM23zI9Zd86AhoEk5efV6z5tNKvub9QrYo7AB7uXrhXrVH0XKN6Lfgp8wwe1b3Iy3feM5g3c9a6Nm2/ZwtnG2uT9tWlZR0/dFm5tXJzfqteg4nUtFVQSQeJ4z/uL3Yw231oI2H3PwKUfQB3lIochAuzO0NuU235biU/nPyGOf/zEgDDH3uLTalLGBX9Pk9FvMQ7y/7AL7/l8Ieo6xf2/NdfUgH4/JsF5OXnWtqwlSX97H5aBTzIb7m/cPnqOX7KPMN9XtbeSm7nwc9YtfU9ABr53E9gwzYsTp7CkMjxDIuazJuLB/HrbzkMfWQSAJ/+70esT5nFz9kX+Tn7EqNjkqyMbwwT93smMmlfXVpWW2rlxvygeqkINW0VdPNBIqxFFMu+fLvYwWz3oS+I7XH9nmYlHcBbNXnQsvw3Kukg/MOpncT2GMvmPcudNreze6jdIB5qN6jYtMKxq1enMdOe31Tico92jKvsaHdkcOTfix7PGescXxTuEvIkXUKK/yBiSOR4AALqt7rloPJY+HAeCx/uqHh3jbtpv2cik/bVttTKjflB9VIRTnXvUWdxp/dY27JnBT3bVvxK285wX7jbyW7lvQPLczfdNxCco0YqysTMYN29R23hjHVtxX6vLHfzWIPj99X2fi9WJL897z3qKPfEvUfvFrez43IWJmVfsULffZC7j6l1bdK+o5CpYw1mjveNTM9vFX08WoJavzNvu1ZltnLbkyZNsuTegaaOterajG3fa3Vt5X9jjbVjlr1Td7Lt0NDQCi9z9MQZAJr6Nyz2uLK3aws1bSVo+ZDVCSrOxMymMnWsTcxtYmZTaawdx8SxNjEzQGJiYoWXeWXqhwC8/fJzxR47A308KiIiImIANW1irFmzZlkdQcTuVNeOo7EW06hpE2MFBwdbHUHE7lTXjqOxFtOoaRNj9ezZ0+oIInanunYcjbWYRk2biIiIiAHUtImIiIgYQJf8EKdV3pXIJ02aVClXKxepLLbUq+raPjTWcjfSmTYx1muvvWZ1BBG7U107jsZaTKOmTURERMQAatpEREREDKCmTURERMQAatpEREREDKCmTURERMQAatpEREREDKCmTURERMQONm/eTHBwMM2bNyc+Pp68vDy7rl9Nm4iIiMgdys/PJz4+npUrV3L48GEyMzNZtGiRXbehpk1ERETkDn3zzTf4+vrSqlUrAIYPH87q1avtug01bSIiIiJ36NSpU/j5+RX97e/vz8mTJ+26Dd17VERERO5ZlzOvsmjNRvLy84tN/+f81bc8dnV1YUj0I9StXeuW9RQUFFRuUHSmTURERO5hdbxq8kBzf86c+4kz534qmn7z4zPnfqJlU/8SGzYAPz+/YmfWTpw4QePGje2aVU2biIiI3NN6PdgOv4b1ypynUQMfIru0L/X5Dh06cOrUKQ4cOADARx99RExMjF1zqmkTERGRe5qbqytPPdGLqlXcSny+ShU3Bj7xEG5upbdNbm5uzJ07l9jYWJo1a0bNmjX5/e9/b9ecLgWO+BBWRERExMlt272fdRu/vmX6/3m4C13DQixIVJzOtImIiIgAndu1okVg8e+hNQ9oxIPtgy1KVNxd07T98ssvVkcQERERg7m4uBD7WE88qrsDUN29GgP69MTVxcXiZNc5bdO2d+9e+vfvj4+PD9WrV+f+++9n/PjxAMTFxeHj48POnTvp2bMnnp6e/PGPf7Q4sYiIiJjOq5Yn0VHdAIiO6kZtr5oWJ/oPp7xO265du+jRowd+fn5MmzaNgIAAjh07RkpKStE82dnZxMbGMmbMGF577TWqVCn/pbwy9cPKjC0iIiJ3kWXrN7Fs/aZK387bLz9n03xO2bSNHTsWT09PduzYQe3atYumDx8+vOhxTk4OU6dOZeDAgVZEFBEREXEop2vasrOz2bp1K6NGjSrWsN3MxcWF6OjoCq3b1k5WRERExNk4XdN26dIl8vPzy72KsLe3N+7u7hVatz4eFREREWdj60klp/shgre3N66urmRkZJQ5n4uT/JJDRERExBGc8uK6vXr14sCBAxw6dAgvL69bno+Li+OTTz7hwoULFqQTERERcTynO9MG8O6773L16lU6d+7MggUL+PLLL1m4cCHx8fFWRxMRERGxhNN9pw0gLCyMbdu28eqrr/Liiy9y7do1/Pz8ePrpp62OJiIiImIJp/x4VERERESKc8qPR0VERESkODVtIiIiIgZQ0yYiIiJiADVtIiIiIgZQ0yYiIiJiADVtIiIiIgZQ0yYiIiJiADVtIiIiIgZQ0yYiIiJiADVtIiIiIgZQ0yYiIiJiADVtIiIiIgZQ0yYiIiJiADVtIiIiIgZQ0yYiIiJiADVtIiIiIgZQ0yYiIiJiADVtIiIiIgZQ0yYiIiJiADVtIiIiIgZQ0yYiIiJiADVtIiIiIgZQ0yYiIiJiADVtIiIiIgZQ0yYiIiJiADVtIiIiIgZQ0yYiIiJiADVtIiIiIgb4/57N4MtJLtynAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 804.38x144.48 with 1 Axes>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qt.draw(output = 'mpl', filename = 'Single_Clifford_gates_circuit_transpiled.jpg', fold = 30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "noise_model = get_noise(MEASURE_NOISE, SINGLE_GATE_DEPOLARISING, TWO_GATE_DEPOLARISING, SINGLE_GATE_SET, TWO_GATE_SET)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "For 2 gates mean simulated ground state population probability is 0.9944\n",
      "The standard error is 0.0002\n",
      "For 4 gates mean simulated ground state population probability is 0.9932\n",
      "The standard error is 0.0002\n",
      "For 8 gates mean simulated ground state population probability is 0.9918\n",
      "The standard error is 0.0002\n",
      "For 12 gates mean simulated ground state population probability is 0.9899\n",
      "The standard error is 0.0002\n",
      "For 16 gates mean simulated ground state population probability is 0.9887\n",
      "The standard error is 0.0002\n",
      "For 20 gates mean simulated ground state population probability is 0.9870\n",
      "The standard error is 0.0002\n"
     ]
    }
   ],
   "source": [
    "# initialise lists for graphs\n",
    "mean_probability = []\n",
    "standard_error_list = []\n",
    "number_of_gates_simulation = []\n",
    "\n",
    "loop_number = 0\n",
    "for gate_depth in SIMULATION_GATE_DEPTH_LIST:\n",
    "    simulation_probability = []\n",
    "    for iteration in range(ITERATIONS):\n",
    "        qt = transpile(rb_circs[iteration][loop_number], basis_gates = BASIS_GATE_SET)\n",
    "        if NOISE:\n",
    "            result = execute(qt, SIMULATOR, noise_model=noise_model, shots=SHOTS).result()\n",
    "        else:\n",
    "            result = execute(qt, SIMULATOR, shots=SHOTS).result()     \n",
    "        counts = result.get_counts(qt)\n",
    "        probability = counts.get('0') / SHOTS\n",
    "        simulation_probability.append(probability)      \n",
    "    number_of_gates_simulation.append(gate_depth + 1)\n",
    "    # need to add one gate to allow for the reversal\n",
    "    mean = mean_of_list(simulation_probability)\n",
    "    mean_probability.append(mean)\n",
    "    standard_deviation, standard_error = calculate_standard_error(simulation_probability)\n",
    "    standard_error_list.append(standard_error)\n",
    "    true_gate_depth = (gate_depth + 1)    \n",
    "    print(f'For {true_gate_depth} gates mean simulated ground state population probability is {mean:.4f}')\n",
    "    print(f'The standard error is {standard_error:.4f}')\n",
    "    loop_number = loop_number + 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "For gate_depth 0\n",
      "The best fit line for Egan 1 ion chain is 0.9939 and for Egan 15 ion chain is 0.9951\n",
      "A calculation based on the number of Clifford gates gives 0.9954\n",
      "For gate_depth 4\n",
      "The best fit line for Egan 1 ion chain is 0.9932 and for Egan 15 ion chain is 0.9937\n",
      "A calculation based on the number of Clifford gates gives 0.9940\n",
      "For gate_depth 6\n",
      "The best fit line for Egan 1 ion chain is 0.9928 and for Egan 15 ion chain is 0.9930\n",
      "A calculation based on the number of Clifford gates gives 0.9933\n",
      "For gate_depth 8\n",
      "The best fit line for Egan 1 ion chain is 0.9925 and for Egan 15 ion chain is 0.9923\n",
      "A calculation based on the number of Clifford gates gives 0.9926\n",
      "For gate_depth 12\n",
      "The best fit line for Egan 1 ion chain is 0.9918 and for Egan 15 ion chain is 0.9909\n",
      "A calculation based on the number of Clifford gates gives 0.9913\n",
      "For gate_depth 16\n",
      "The best fit line for Egan 1 ion chain is 0.9910 and for Egan 15 ion chain is 0.9895\n",
      "A calculation based on the number of Clifford gates gives 0.9899\n",
      "For gate_depth 20\n",
      "The best fit line for Egan 1 ion chain is 0.9903 and for Egan 15 ion chain is 0.9881\n",
      "A calculation based on the number of Clifford gates gives 0.9885\n"
     ]
    }
   ],
   "source": [
    "theoretical_probability = []\n",
    "egan_one_ion_noise = []\n",
    "egan_fifteen_ion_noise = []\n",
    "number_of_gates_calculation = []\n",
    "\n",
    "for gate_depth in CALCULATION_GATE_DEPTH_LIST:\n",
    "    gate_noise = ((1 - SINGLE_CLIFFORD_GATE_NOISE) ** gate_depth) \n",
    "    theory = gate_noise * (1 - MEASURE_NOISE)  \n",
    "    egan_one_ion = 0.99388 - gate_depth * 0.000177\n",
    "    egan_fifteen_ion = 0.9951 - gate_depth * 0.000348\n",
    "    theoretical_probability.append(theory)\n",
    "    egan_one_ion_noise.append(egan_one_ion)\n",
    "    egan_fifteen_ion_noise.append(egan_fifteen_ion)\n",
    "    number_of_gates_calculation.append(gate_depth) \n",
    "    \n",
    "    print(f'For gate_depth {gate_depth}')  \n",
    "    print(f'The best fit line for Egan 1 ion chain is {egan_one_ion:.4f} and for Egan 15 ion chain is {egan_fifteen_ion:.4f}')\n",
    "    print(f'A calculation based on the number of Clifford gates gives {theory:.4f}')   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot graph:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABiv0lEQVR4nO2dd3iUVfbHP4ck9CodQq8hhYRepRcF6QgsKNiVte6KruuK6Kpr/dnL2sCCgokKKLoKCgIi0gSk1xB6CRAIJKSd3x93MpkkM+mTxv08Tx4yd+77vuedDPOdc8v3iKpisVgsFktOKVPUAVgsFoulZGGFw2KxWCy5wgqHxWKxWHKFFQ6LxWKx5AorHBaLxWLJFVY4LBaLxZIrrHBYSjwiMllEfiygcy0XkVsL4ly5vG6kiAz08FxvEdlV2DFZLJ6wwmEpEYhILxFZLSIxInJGRH4Vkc4AqjpXVQcXdYzeQlVXqmqb1MdZiYzFUhj4FnUAFkt2iEhV4FvgLuALoCzQG7hclHGVBkTEV1WTsmvL7TkspRubcVhKAq0BVPVzVU1W1ThV/VFVtwCIyDQRWZXaWURURO4UkT0ick5E3hQRcTznIyIvichpETkgInc7+rv9EiUiN4vIDhE5KyI/iEgTT0GKyA0iclBEokXkUdfMQETmiMhTLn37isjhDKfoLCLbHdeaLSLlM/YVkU+AxsA3IhIrIg95iGW4iGxy3P9qEQlxeS5SRB4WkS3ARRFp6XgNbhGRKOBnESkjIv9y3M9JEflYRKo5jm/qpn95EfnUce/nRGSdiNT19FpZSjZWOCwlgd1Asoh8JCLXiEiNHBwzHOgMhADXA0Mc7bcB1wChQAdglKcTiMhI4J/AGKA2sBL43EPfdsDbwA1AA6Am4J+DOF2Z7IizBUYs/5Wxg6reAEQB16lqZVV93k0sYcCHwB2OOP4LLBKRci7dJgHDgOpAarbQBwhwxDDN8dMPaA5UBt7IcCnX/lOBakAjxzXvBOJyfuuWkoQVDkuxR1XPA70ABd4DTonIomy+0T6rqudUNQpYhhEKMCLyqqoeVtWzwLNZnONO4D+qusMxFPMMEOoh6xgHfKuqK1T1MvAYkJKL2wR4Q1UPqeoZ4GnMh3teuB34r6r+7sjQPsIM63Vz6fOa41quH+6zVPWio20y8H+qul9VY4FHgIkZMjPX/okYwWjpuOYGx9/NUgqxwmEpETg+vKepqj8QhPlW/0oWhxx3+f0S5hszjuMOuTzn+ntGmgCvOoZezgFnAAEauumb7ryqehGIzuLc7nCN5aDjnHmhCfD31LgdsTfKcD539+3a1sARg2s8vkBdD/0/AX4A5onIURF5XkT88hi/pZhjhcNS4lDVncAcjIDklmOkH0JqlEXfQ8Adqlrd5aeCqq72cF7nuUSkIuYbeCoXgYouj+u5OYdrLI2Box7iys7S+hDwdIa4K6qq6zCbu3O4th3FCJBrPEnACXf9VTVRVZ9Q1XZAD8xQ4Y3ZxGkpoVjhsBR7RKStiPxdRPwdjxthhnHW5OF0XwD3iUhDEakOPJxF33eAR0Qk0HHdaiIy3kPfCGC4Y9lwWeBJ0v//2gRcKyJXiUg94H435/iriPiLyFXAo8B8D9c6gZl38MR7wJ0i0lUMlURkmIhUyeKYjHwOPCAizUSkMmaYbr6n1VMi0k9EgkXEBziPGbrK7VCdpYRghcNSErgAdAV+F5GLGMHYCvw9D+d6D/gR2AL8AXyH+SadnLGjqn4NPIcZfjnvuOY17k6qqtuAvwKfYbKPs4DrqqlPgM1ApOP67kThM8dz+4F9wFNu+gD8B/iXYxjqQTexrMcsAnjDEcdezER3bvjQEfMK4AAQD9yTRf96GPE8D+wAfnEcbymFiC3kZLmSEZFrgHdU1eMy23ycOxK4VVWXFvS5LZaixGYclisKEakgIteKiK+INAQeB74u6rgslpKEV4VDRD50bB7a6uF5EZHXRGSviGwRkQ4uz00Vs4Frj4hMdWnvKCJ/Oo55TcRs7LJYcogAT2CGcP7ADKvMLNKILJYShleHqkTkaiAW+FhVM62AEZFrMeOm12LGsF9V1a6OycH1QCfMyo0NQEdVPSsia4F7gd8x49Ovqer3XrsJi8VisaTDqxmHqq7ArH33xEiMqKiqrgGqi0h9zE7UJap6xrFJawkw1PFcVVVdo0bxPiaLnb8Wi8ViKXiK2uSwIek3ER12tGXVfthNeyZE5HbMDloqVarUsW3btgUXtcVisVwBbNiw4bSq1s7YXtTC4TVU9V3gXYBOnTrp+vXrizgii8ViKVmIyEF37UW9quoI6XfL+jvasmr3d9NusVgslkKiqIVjEXCjY3VVNyBGVY9hPG8Gi0gNhxPqYOAHx3PnRaSbYzXVjcDCIoveYrFYrkC8OlQlIp8DfYFaYuoJPA74AajqO5hVUddidrZeAm5yPHdGRP4NrHOc6kmHYyjAdIxPUQXge8ePxWKxWAqJK2LnuJ3jyD2JiYkcPnyY+Pj4og7FYrF4mfLly+Pv74+fX3pDYxHZoKqdMvYvtZPjlvxx+PBhqlSpQtOmTbF7LC2W0ouqEh0dzeHDh2nWrFmOjinqOQ5LMSU+Pp6aNWta0bBYSjkiQs2aNXM1umCFw+IRKxoWy5VBbv+vW+GwWCwWS66wwmEpthw/fpyJEyfSokULOnbsyLXXXsvu3buzPKZy5cpZPu+J5cuXM3z48Cz7bNq0ie+++875eNGiRTz7bFYly3NGZGQkQUF5KWaYN/r27YtdLJKe77//nk6dOtGuXTvCwsL4+99NqZdZs2bx4osvAjBz5kyWLjUO+StXriQwMJDQ0FDi4uKYMWMGgYGBzJgxI88x5PW9m5EFCxawffv2AjmXJ+zkuKVYoqqMHj2aqVOnMm/ePAA2b97MiRMnaN26dZHEtGnTJtavX8+1114LwIgRIxgxYkSRxFISSUpKwtfX1+PjnB5X0HFs3bqVu+++m8WLF9O2bVuSk5N59913Mx335JNPOn+fO3cujzzyCFOmTAHg3Xff5cyZM/j4+OQphoJkwYIFDB8+nHbt2nnl/GAzDksxZdmyZfj5+XHnnXc629q3b0/v3r2JjY1lwIABdOjQgeDgYBYudL8H9LnnniM4OJj27dvzj3/8A0j/bfv06dM0bdo003Fr166le/fuhIWF0aNHD3bt2kVCQgIzZ85k/vz5hIaGMn/+fObMmcPdd98NmKyhf//+hISEMGDAAKKiogCYNm0a9957Lz169KB58+ZERES4jTUpKYnJkycTEBDAuHHjuHTpEmA+rDp37kxQUBC33347qcvnX3vtNdq1a0dISAgTJ04E4OLFi9x888106dKFsLAw5+sSFxfHxIkTCQgIYPTo0cTFxbmN4aeffiIsLIzg4GBuvvlmLl++DEDTpk15/PHHna/3zp07Mx2bnJzMjBkz6Ny5MyEhIfz3v/8FTCbXu3dvRowYQbt27TI9jo+P56abbiI4OJiwsDCWLVsGwJw5cxgxYgT9+/dnwIAB6a4VGRlJ27Zt3b5eGzZsoE+fPnTs2JEhQ4Zw7Ngx59/9/vvvp1OnTrz66qvpzvf888/z6KOPkupn5+Pjw1133ZXpHqdNm0ZERATvv/8+X3zxBY899hiTJ09mxIgRxMbG0rFjR+bPn5/le+HOO++ka9euPPTQQxw4cIDu3bsTHBzMv/71L7d/E4B///vftGnThl69ejFp0iRnBvTee+/RuXNn2rdvz9ixY7l06RKrV69m0aJFzJgxg9DQUPbt28e+ffsYOnQoHTt2pHfv3m7/frlGVUv9T8eOHdWSO7Zv3572YP19qkv6FOzP+vuyvP6rr76q999/v9vnEhMTNSYmRlVVT506pS1atNCUlBRVVa1UqZKqqn733XfavXt3vXjxoqqqRkdHq6pqnz59dN26dc5jmzRpoqqqy5Yt02HDhqmqakxMjCYmJqqq6pIlS3TMmDGqqjp79mz961//6ozD9fHw4cN1zpw5qqr6wQcf6MiRI1VVderUqTpu3DhNTk7Wbdu2aYsWLTLdz4EDBxTQVatWqarqTTfdpC+88EK6uFVVp0yZoosWLVJV1fr162t8fLyqqp49e1ZVVR955BH95JNPnG2tWrXS2NhYfemll/Smm25SVdXNmzerj4+P8zVIJS4uTv39/XXXrl2qqnrDDTfoyy+/rKqqTZo00ddee01VVd9880295ZZbMt3Df//7X/33v/+tqqrx8fHasWNH3b9/vy5btkwrVqyo+/fvd77Oro9ffPFFZ2w7duzQRo0aaVxcnM6ePVsbNmyY7v6ze70SEhK0e/fuevLkSVVVnTdvnvPcffr00bvuuivTuVRVw8LCdNOmTW6fe/zxx51/i6lTp2p4eHim31XT3neqWb8Xhg0bpklJSaqqet111+lHH32kqqpvvPFGunOksnbtWm3fvr3GxcXp+fPntWXLls54Tp8+7ez36KOPOv9GGWPr37+/7t69W1VV16xZo/369XN7r+n+zzsA1qubz1SbcVhKHKrKP//5T0JCQhg4cCBHjhzhxIkT6fosXbqUm266iYoVKwJw1VVX5fj8MTExjB8/nqCgIB544AG2bduW7TG//fYbf/nLXwC44YYbWLVqlfO5UaNGUaZMGdq1a5cpzlQaNWpEz549AZgyZYrz+GXLltG1a1eCg4P5+eefnbGEhIQwefJkPv30U+eQx48//sizzz5LaGgoffv2JT4+nqioKFasWOEcUgkJCSEkJCTT9Xft2kWzZs2cw4BTp05lxYoVzufHjBkDQMeOHYmMjMx0/I8//sjHH39MaGgoXbt2JTo6mj179gDQpUuXdPsDXB+vWrXKGVvbtm1p0qSJcx5r0KBBHv9u7l6vXbt2sXXrVgYNGkRoaChPPfUUhw+nmWlPmDDB7bkKmqzeC+PHj3cOZ/36669MmjTJ2c8dv/76KyNHjqR8+fJUqVKF6667zvnc1q1b6d27N8HBwcydO9ft+zQ2NpbVq1czfvx4QkNDueOOO5xZWH6wcxyW7On4SqFfMjAw0OOwzty5czl16hQbNmzAz8+Ppk2b5ngNuq+vLykpKQAej3nsscfo168fX3/9NZGRkfTt2zdP95BKuXLlnL+rB6eGjMshRYT4+HimT5/O+vXradSoEbNmzXLGvHjxYlasWME333zD008/zZ9//omq8uWXX9KmTZt8xZvVPfj4+JCUlJTpeVXl9ddfZ8iQIenaly9fTqVKldK1ZXzsiaz6uXu9VJXAwEB+++23XJ0vMDCQDRs20L59+xzFlR8yxpCfJe/Tpk1jwYIFtG/fnjlz5rB8+fJMfVJSUqhevTqbNm3K83XcYTMOS7Gkf//+XL58Od0k5ZYtW1i5ciUxMTHUqVMHPz8/li1bxsGDmZ2fBw0axOzZs51j32fOGKuzpk2bsmHDBgCPwhQTE0PDhqbMy5w5c5ztVapU4cKFC26P6dGjh3MSf+7cufTu3TtX9xsVFeX8wPvss8/o1auXUyRq1apFbGysM96UlBQOHTpEv379eO6554iJiSE2NpYhQ4bw+uuvO8Xpjz/+AODqq6/ms88+A8y31C1btmS6fps2bYiMjGTv3r0AfPLJJ/Tp0yfH8Q8ZMoS3336bxMREAHbv3s3FixezPa53797MnTvXeUxUVFSOhM/d69WmTRtOnTrlbE9MTMxRtjhjxgyeeeYZZ6aTkpLCO++8k+1xnsjpe6Fnz57p+nnq88033xAfH09sbCzffvut87kLFy5Qv359EhMT0x3v+j6tWrUqzZo1Izw8HDACv3nz5jzfWypWOCzFEhHh66+/ZunSpbRo0YLAwEAeeeQR6tWrx+TJk1m/fj3BwcF8/PHHuCvSNXToUEaMGEGnTp0IDQ11Tig++OCDvP3224SFhXH69Gm3137ooYd45JFHCAsLS/ftul+/fmzfvt05Oe7K66+/zuzZswkJCeGTTz7JNAGbHW3atOHNN98kICCAs2fPctddd1G9enVuu+02goKCGDJkCJ07dwbMRPSUKVOcE8r33nsv1atX57HHHiMxMZGQkBACAwN57LHHALjrrruIjY0lICCAmTNn0rFjx0zXL1++PLNnz2b8+PEEBwdTpkyZdAsTsuPWW2+lXbt2dOjQgaCgIO644w63mUlGpk+fTkpKCsHBwUyYMIE5c+aky9By83qVLVuWiIgIHn74Ydq3b09oaCirV6/O9lwhISG88sorTJo0iYCAAIKCgti/f3+O7tsdOX0vvPrqq7z55psEBwdz5Ij76hCdO3dmxIgRhISEcM011xAcHEy1atUAM2netWtXevbsme7/wMSJE3nhhRcICwtj3759zJ07lw8++ID27dsTGBjocTFJbrAmhxa37Nixg4CAgKIOw2LJRGRkJMOHD2fr1q1FHUqhEBsbS+XKlbl06RJXX3017777Lh06dCjw67j7P29NDi0Wi6UEcvvtt7N9+3bi4+OZOnWqV0Qjt1jhsFgsJYqmTZteMdkG4JyfKk7YOQ6LxWKx5AorHBaLxWLJFVY4LBaLxZIrrHBYLBaLJVdY4bAUW3x8fAgNDXX+FISFeUbcmRwWV1wtvguLpk2but3vEh4eTkBAAP369cvzuW+++Wbq1KmTpaX8O++8w8cff5zna2RHbi3mjx49yrhx47wWDxScXb83sauqLMWWChUqFLhVQknCm9bb+eWDDz7gvffeo1evXjnq7+5epk2bxt13382NN97o8bjcbEIsDBo0aODRcaCgKAl2/TbjsJQ4vvvuO9q2bUvHjh259957nQWY3Nmhg7ENGTNmDEOHDqVVq1Y89NBDznPVrl0bMJbkw4YNo3379gQFBWXaGQ6eLdnnzJnDyJEj6du3L61ateKJJ54ACt7+G0xNku7du9OqVSvee+89wNhIzJgxg6CgIIKDg52xZyxOdffddzstVDxZpUdHRzN48GACAwO59dZb3XprPfnkk6xatYpbbrmFGTNm5MkaHYwVSnbmk65Z1qZNm+jWrRshISGMHj2as2fPOl+zhx9+mC5dutC6dWtWrlzp9lzubPbBZE8Zj42MjKR379506NCBDh06OHeguxbdyup95Yqn1/rMmTOMGjWKkJAQunXr5rSCcbXrDw8PJygoiPbt23P11VcDni3sC5Pi+XXGUqy4//77C/ybf2hoKK+88kqWfeLi4ggNDXU+fuSRRxg5ciR33HEHK1asoFmzZk53UTDuqitXrsTX15elS5fyz3/+ky+//BIwHzp//PEH5cqVo02bNtxzzz00atSIdevWAfC///2PBg0asHjxYsD4VeWGtWvXsnXrVipWrEjnzp0ZNmwYtWrVYteuXXzwwQf07NmTm2++mbfeeov77ruPe+65h4ULF1K7dm3mz5/Po48+yocffghAQkKCx+GTLVu2sGbNGi5evEhYWBjDhg3jt99+Y9OmTWzevJnTp0/TuXNn54dMVtSqVYuNGzfy1ltv8eKLL/L+++/zxBNP0KtXL2bOnMnixYv54IMPMh03c+ZMfv75Z1588UU6derESy+9hIjw559/snPnTgYPHuz0fdq4cSNbtmzJlTuxJ2688UZef/11+vTpw8yZM3niiSec76GkpCTWrl3Ld999xxNPPOGs1JfK999/z8KFC/n999+pWLGi07vM07F16tRhyZIllC9fnj179jBp0iS3fxNP76uMuHutH3/8ccLCwliwYAE///wzN954Y6b/Z08++SQ//PADDRs25Ny5c4DJ9qpVq8a6deu4fPkyPXv2ZPDgwekciL2NFQ5LscXdUNWmTZto3ry58z/JpEmTnEaIMTExTJ06lT179iAiTsM9gAEDBjg9ftq1a8fBgwfT/QcPDg7m73//Ow8//DDDhw/PtUnhoEGDqFmzJmAsyFetWsWoUaMy2X+/9tprDB061Gn/DeYbZP369Z3nysr+e+TIkVSoUIEKFSrQr18/1q5dy6pVq5g0aRI+Pj7UrVuXPn36sG7dOqpWrZplzK5W6V999RUAK1ascP4+bNgwatSoke29r1q1invuuQfInTV6boiJieHcuXNO48WpU6cyfvx4t/fizvY9K5t9d8cmJiZy9913s2nTJnx8fDyWLM7ufeXuGqmv76pVq5xfbPr37090dDTnz59Pd1zPnj2ZNm0a119/vfMcP/74I1u2bHEOmcXExLBnzx4rHJbiRXaZQXEhKzt0V+M8d9bgrVu3ZuPGjXz33Xf861//YsCAAcycOTNdn6ws2d3ZfHtqz6v9d1bXcYdrvO5izs4qvSDIqYV6fsnPvbg79uWXX6Zu3bps3ryZlJQUypcvn+Wx2V07r/G98847/P777yxevJiOHTuyYcMGjxb2hYmd47CUKNq0acP+/fud3wxd5yI82aHnhKNHj1KxYkWmTJnCjBkz2LhxY6Y+WVmyL1myhDNnzhAXF8eCBQucWUZB2n8DLFy4kPj4eKKjo1m+fDmdO3emd+/ezJ8/n+TkZE6dOsWKFSvo0qULTZo0Yfv27Vy+fJlz587x008/ZXt+Vwv277//3jmPkBV5tUbPDdWqVaNGjRrOOYjc2r57stn3RExMDPXr16dMmTJ88sknJCcn5z14D7i+bsuXL6dWrVqZssR9+/bRtWtXnnzySWrXrs2hQ4fybGFfkNiMw1JsyTjHMXToUJ599lneeusthg4dSqVKlZxW42Ds0KdOncpTTz3FsGHDcnWtP//8kxkzZlCmTBn8/Px4++23M/V58MEHuf7663n33Xcznb9Lly6MHTuWw4cPM2XKFDp16kRkZKTT/vvmm2+mXbt26ey/7733XmJiYkhKSuL+++8nMDAw2zhDQkLo168fp0+f5rHHHqNBgwaMHj2a3377jfbt2yMiPP/889SrVw+A66+/nqCgIJo1a0ZYWFi253/88ceZNGkSgYGB9OjRg8aNG2d7zPTp07nrrrsIDg7G19c3x9bokyZNYvny5Zw+fRp/f3+eeOIJbrnlFo/9P/roI+68804uXbpE8+bNmT17drbXSGXo0KFs2rSJTp06UbZsWa699lqeeeaZLO9p7NixfPzxx873WkEza9Ysbr75ZkJCQqhYsSIfffRRpj4zZsxgz549qCoDBgygffv2hISEEBkZSYcOHVBVateuzYIFCwo8vqywtuoWtxRnW/VUm2lV5a9//SutWrXigQceKLJ45syZw/r163njjTfStV9p9t+Wkk1ubNXtUJWlxPHee+8RGhpKYGAgMTEx3HHHHUUdksVyRWEzDotbinPGYbFYCh6bcVgsFovFa1jhsFgsFkuusMJhsVgsllxhhcNisVgsucIKh6XYUtxs1adNm+Z1Z1R3uJreFRae7MZXrlxJYGAgoaGhxMXF5encjz76KI0aNaJy5crp2ufMmUPt2rWdf+/333/f7fE9evTI03VzQkZjyJwwc+bMTN5YBc21117r9KoqDtgNgJZiy5Vqq16c7dTnzp3LI488wpQpU3LU3929XHfdddx99920atUqU/8JEyZk2g+TkVSn2uLCk08+6fVrfPfdd16/Rm6wGYelxFFUtuqu/PTTT4SFhREcHMzNN9/M5cuXWbdundOIbuHChVSoUIGEhATi4+Np3rw5YCwkhg4dSseOHendu7fTYnvatGnceeeddO3a1a0996FDhzLZtgP83//9H0FBQQQFBTk9xVytvwFefPFFZs2aBXi2II+Li2PixIkEBAQwevRot9nE+++/zxdffMFjjz3G5MmTs7Rz7927NyNGjKBdu3aZztOtW7d0po65JTVTyer6ffv2Zdy4cU5be3fbDvbu3cvAgQNp3749HTp0YN++fYDZYOru2CeffJLOnTsTFBTE7bff7mx3zUQ9Wai7ktX78fPPPyc4OJigoCAefvhhZ3tqQS1P71NPNv1eQ1VL/U/Hjh3Vkju2b9+e9uC++1T79CnYn/vuyzaGMmXKaPv27Z0/8+bN07i4OPX399f9+/erqurEiRN12LBhqqoaExOjiYmJqqq6ZMkSHTNmjKqqzp49W5s1a6bnzp3TuLg4bdy4sUZFRaW7VkREhN56663Ox+fOncsUz9SpUzU8PNwZw65du1RV9YYbbtCXX35ZExMTtVmzZqqq+ve//107deqkq1at0uXLl+vEiRNVVbV///66e/duVVVds2aN9uvXz3nuYcOGaVJSUqbrzp49W+vVq6enT5/WS5cuaWBgoK5bt07Xr1+vQUFBGhsbqxcuXNB27drpxo0b9cCBAxoYGOg8/oUXXtDHH39cVVX79Omjf/vb31RVdfHixTpgwABVVX3ppZf0pptuUlXVzZs3q4+Pj65bt87ja5D6mg0cOFCTkpL0+PHj2qhRIz169KguW7ZMK1as6PwbeaJSpUpu7zM4OFjHjh2b6W+U8bisrl+1alU9dOiQJicna7du3XTlypWZztOlSxf96quvVFU1Li5OL168mOWx0dHRzmOnTJmiixYtyvSaNGnSRF977TVVVX3zzTf1lltuyXRdT+/HI0eOaKNGjfTkyZOamJio/fr106+//tp53lOnTrl9nyYkJGj37t315MmTqqo6b948598yN6T7P+8AWK9uPlNtxmEptqQOVaX+TJgwgZ07d2ayVU8lJiaG8ePHExQUxAMPPJDOODDV/rp8+fJO+2tXgoODWbJkCQ8//DArV650WmW7Y9euXTRr1ozWrVsDxuJ7xYoV+Pr60qJFC3bs2MHatWv529/+xooVK1i5ciW9e/cmNjaW1atXM378eEJDQ7njjjvSfTMcP348Pj4+bq+ZatteoUIFp237qlWrGD16NJUqVaJy5cqMGTPGYxEjV9zZiK9YscI5/BQSEkJISEi25/Fk5w7Guyu3Nt/XXXcdkZGRbNmyhUGDBjF16tR8Xd/f358yZcoQGhqayWr9woULHDlyhNGjRwNQvnx5p+W6p2OXLVtG165dCQ4O5ueff/ZoTJmdxTu4fz+uW7eOvn37Urt2bXx9fZk8eTIrVqxId5y79+muXbucNv2hoaE89dRTHD58OMvXLr8Uz4FUS/HC2qrnmKuvvprvv/8ePz8/Bg4cyLRp00hOTuaFF14gJSWF6tWre5y3udLt1FPrmQDceuutHivq5YSc2p3n9Nj4+HimT5/O+vXradSoEbNmzcr0umY8Pic267mNz937dPTo0Vna9HsDm3FYShRFaavuGkNkZCR79+4F0lt89+7dm1deeYXu3btTu3ZtoqOj2bVrF0FBQVStWpVmzZoRHh4OmGHizZs35yg+d7btvXv3ZsGCBVy6dImLFy/y9ddf07t3b+rWrcvJkyeJjo7m8uXLfPvtt9me39VOfevWrc4yplnhyc49r7hmX4sWLcrW8iY/169SpQr+/v5OV9nLly87LdfdkSoStWrVIjY21iur67p06cIvv/zC6dOnSU5O5vPPP89kHe/ufZofm/68YjMOS7GluNmqp1K+fHlmz57N+PHjSUpKonPnztx5550AdO3alRMnTjhLt4aEhHD8+HFnhjB37lzuuusunnrqKRITE5k4cSLt27fPNj53tu1gJmZTPyxvvfVWp3X6zJkz6dKlCw0bNqRt27bZnv+uu+7ipptuIiAggICAADp27JjtMZ7s3N1NCLvy0EMP8dlnn3Hp0iX8/f259dZbmTVrFq+99hqLFi3C19eXq666Klvxz+v1U/nkk0+44447mDlzJn5+fk5Bd0f16tW57bbbCAoKol69eunedwVF/fr1efbZZ+nXrx+qyrBhwxg5cmS6Pu7ep/mx6c8r1uTQ4pbibHJY3GzVLZbSgDU5tJRqrK26xVK0ZDtUJSI+qlrwdRMtljzywAMP2AzDYilCcpJx7BGRF0Qk804ei8VisVxx5EQ42gO7gfdFZI2I3C4iVbM7yGKxWCylk2yFQ1UvqOp7qtoDeBh4HDgmIh+JSEuvR2ixWCyWYkW2wiEiPiIyQkS+Bl4BXgKaA98Axct5y1Kk9O1rfiwWS+kmR3McwEjgBVUNU9X/U9UTqhoB/M+74VmuZJ5++mkCAwMJCQkhNDSU33//HTD7FbZv314g10g1j8uKZ555Jtfn9ZYVel7Om9G+PBVXO3BXG/VUC+9z587x1ltv5S9gS6kkJxsAb1TVVa4NItJTVX9V1Xu9FJflCue3337j22+/ZePGjZQrV47Tp0+TkJAA4LFOg7d45pln+Oc//1lo1yssW3VPduCpFt6RkZG89dZbTJ8+3euxWEoWOck4XnPT9npBB2KxuHLs2DFq1arl9PSpVasWDRo0ANJ/O65cuTIzZswgMDCQgQMHsnbtWvr27Uvz5s1ZtGgRkPlb+vDhw1m+fHmma44aNYqOHTsSGBjIu+++C8A//vEP5w72yZMnA/Dpp5/SpUsXp1FhcrJZrT579mxat25Nly5d+PXXX93e16xZs7jhhhvo3r07rVq14r333gMyW5HHx8dz0003ERwcTFhYGMuWLXOew5PFurv4U3nggQcIDAxkwIABnDp1CvBcmCo1C/vHP/7Bvn37CA0NZcaMGdx4441Oiw6AyZMns3DhQrf3aSndePxaIyLdgR5AbRH5m8tTVQH3Fp6WUsn990NO6iml9snJPEdoaNbeiYMHD+bJJ5+kdevWDBw4kAkTJmTy7QFTR6N///688MILjB49mn/9618sWbKE7du3M3XqVEaMGJF9MA4+/PBDrrrqKuLi4ujcuTNjx47l2Wef5Y033nAaE+7YsYP58+fz66+/4ufnx/Tp05k7dy6DBg3i8ccfZ8OGDVSrVo1+/fo57T8ysmXLFtasWcPFixcJCwtz2qNs3LiRrVu30qxZM1566SVEhD///JOdO3cyePBgdu/eDZi6I1u3bqVixYp07tyZYcOG0alTJ7fx16xZk4sXL9KpUydefvllnnzySZ544olsiyUBPPvss2zdutV577/88gsvv/wyo0aNIiYmhtWrV/PRRx/l+PW1lB6yyjjKApUx4lLF5ec8MM77oVmuZCpXrsyGDRt49913qV27NhMmTHDrXVS2bFmGDh0KGMvpPn364OfnR3BwsEdLa0+89tprtG/fnm7dunHo0CH27NmTqc9PP/3Ehg0b6Ny5M6Ghofz000/s37+f33//3WmJXbZsWSZMmODxOiNHjqRChQrUqlWLfv36sXbtWiC9FfmqVaucNudt27alSZMmTuFwZ7GeVfxlypRxxjNlyhRn/9zSp08f9uzZw6lTp/j8888ZO3Zssa1UaPEuHv/qqvoL8IuIzFHVg576WUo/OXVVT8003IwC5QkfHx/69u1L3759CQ4O5qOPPmLatGnp+vj5+TkNBMuUKeMc2ipTpozTqjo7m3ET83KWLl3Kb7/9RsWKFenbt6/bfqrK1KlT+c9//pOu3XUIJzs8WaTn1Irc3fE5jd/d8bnhxhtv5NNPP2XevHnMnj07z+exlGw8Zhwi8orj1zdEZFHGn8IJz3KlsmvXrnTf+Ddt2kSTJk3ydK6mTZuyadMmUlJSOHTokPMbvisxMTHUqFGDihUrsnPnTtasWeN8zs/Pj8TERMAU4ImIiODkyZMAnDlzhoMHD9K1a1d++eUXoqOjSUxMzNJpdeHChcTHxxMdHc3y5cvdOq327t2buXPnArB7926ioqJo06YN4N5iPav4U1JSnHMZn332Gb169crR61alShUuXLiQrm3atGnOErXuysJargyyyjM/cfz7YmEEYrG4Ehsbyz333MO5c+fw9fWlZcuWmSZ8c0rPnj1p1qwZ7dq1IyAggA4dOmTqM3ToUN555x0CAgJo06YN3bp1cz53++23ExISQocOHZg7dy5PPfUUgwcPJiUlBT8/P9588026devGrFmz6N69O9WrV09nB5+RkJAQ+vXrx+nTp3nsscdo0KCBcxgqlenTp3PXXXcRHByMr68vc+bMcWZT7izWg4ODPcZfqVIl1q5dy1NPPUWdOnWyraeeSs2aNenZsydBQUFcc801vPDCC9StW5eAgABGjRqVo3NYSifWVt3ilrzYqhf0UFVpZNasWVSuXJkHH3ywqEPJE5cuXSI4OJiNGzdmWV7XUvLIja16Vquq/gQ8qoqqZl+U2HJFYQWjdLN06VJuueUWHnjgASsaVzhZDVUNL7QoLJYrhFmzZhV1CHlm4MCBHDxo18lYsl5VZd8hW2ZB3FFoPB7q9oUyfkUcUOGiqvlagWOxWEoGuZ2yyGpV1SrHvxdE5HzGf/MZZ8kg8Rwc/AyWDYav6sHvt8LR/0FyQlFH5nXKly9PdHR0rt9QFoulZKGqREdHU758+RwfYyfHsyMpDo79AIci4PAiSLoAZWuA/0hoNA7qDQKfsgUbcDEgMTGRw4cPe9wLYLFYSg/ly5fH398fP7/0oyqeJsdzJBwi0gHohZksX6WqfxRQvIVCga2qSo6HY0scIrIQEmPArxo0HGGGs+oPAp+cq7bFYrEUZzwJR07qccwEPgJqArWAOSLyrxxedKiI7BKRvSLyDzfPNxGRn0Rki4gsFxF/l+eeE5Gtjp8JLu1zROSAiGxy/ITmJJY8sXQpLFsGDhM7fMqD/3XQ/SMYcwL6LIZGo+HIN7BiBHxZB1ZPgUMLTKZisVgspZBsMw4R2QW0V9V4x+MKwCZVbZPNcT6YkrODgMPAOmCSqm536RMOfKuqH4lIf+AmVb1BRIYB9wPXAOWA5cAAVT0vInMcx2S29fRAnjOOq6+GlSuhTh0YMwbGjzdtGf15khPgxDI4FA6HvoaEM+BbGRoON8NZDa4B34q5v77FYrEUIXnOOICjgOv4SzngSA6O6wLsVdX9qpoAzMMUhHKlHfCz4/dlLs+3A1aoapKqXgS2AENzcM2C5fvvYf586NMHPv4YBgyABg3gzjvhp5/A4YWET1loMAS6vg9jjkO/H6HpX+D4Ulg1Dr6sDauuh4NfQNLFQr8Ni8ViKUiyWlX1uoi8BsQA2xxDRLOBrcC5HJy7IXDI5fFhR5srm4Exjt9HA1VEpKajfaiIVBSRWkA/oJHLcU87hrdeFpFyHuK/XUTWi8j61PoDuaZSJbj+evjiCzh5EsLDoV8/+OQTGDgQ6teH22+HJUvSRKSMn5nr6PJfGH0M+v8EzafCyV/g1wlGRFaOhch5kHgh6+tbLBZLMcTjUJWITM3qQFXN0ohfRMYBQ1X1VsfjG4Cuqnq3S58GwBtAM2AFMBYIUtVzIvIoMB44BZwE1qnqKyJSHziOsX1/F9inqu5LmTkocMuRS5dMNhIRAd98AxcvQs2aMHq0Gc7q1w8yrE4gJRlOrYKocDj0JcQfN3Mm9Yea4Sz/68CvasHFaLFYLPkkX6uq8njB7sAsVR3iePwIgKr+x0P/ysBOVfV389xnwKeq+l2G9r7Ag6qa5S53r3pVxcXB//5nspFvvoHYWLjqKiMi48aZ4S13InJ6NURFGBGJOwJlykL9IWZ1VsProGx178RrsVgsOSTPwiEirYD/YOYdnHMdqto8m+N8MZPjAzBzIuuAv6jqNpc+tYAzqpoiIk8Dyao60zGxXl1Vo0UkBPgMCFXVJBGpr6rHxGxpfhmIV9VMK7ZcKTSTw7g4+OEHk4ksWgQXLkCNGjBqlMlEBgyAshn2fGgKnF7jEJEIuHTIDHfVG+QQkRFQ7irvx26xWCwZyI9wrAIex3xIXwfcBJRR1Zk5uOi1wCuYUrMfqurTIvIksF5VFzmGs/6D2R+yAvirql4WkfLARsdpzgN3quomxzl/BmoDAmxyPBebVRxF4o4bHw8//mgykUWL4Px5qF7diMi4cTBokHsRiV7nGM6KgIsHQXyh3kBoPA78R0G5moV7HxaL5YolP8KxQVU7isifqhrs2ualWAucIrdVv3w5TUQWLjQiUq0ajBxpMpFBg6Bchjl+VTiz3mQiUeFw8QCID9Tt7xCR0VC+dtHcj8ViuSLIj3Csxuwaj8AsnT0CPJvdPo7iRJELhyuXL5uNheHhsGABxMQYERkxIk1EMnrGqMLZP9JEJHYvSBmo09cMZ/mPhgp1i+JuLBZLKSY/wtEZ2AFUB/4NVAOeV9U1WR1XnChWwuFKQkJ6ETl3DqpUYUmFESyrNZ5nNgxxLyLntjjmRMLh/C4jIrWvNplIozFQoX5R3I3FYill5HtVlYhUBVRVS9zmg2IrHK4kJJhNhRERbP5oI9OSP+SJCs9y3Shf5PrxMGQIVKiQ/hhViNmWNicSsx0QqN3LISJjoWLGrTMWi8WSM/KTcXQCZgNVHE0xwM2quqHAo/QSJUI4XOgUmsT+HZc5m1CJXr6/8VzS3+lR+U+47jozsX7NNZlFBIxwpK7OOvenaavVwwxnNRoLlRplPsZisVg8kB/LkQ+B6araVFWbAn/FCInFS1Su7ktQ10q88w7srdWNnqxmdM0V7PjuAIwdC7Vrw8SJ8OWXZjNiKtXaQfBMuHYLDN8JIU8Zi5OND8DCxvBDN9jxEsRGFtm9WSyWkk9OMo4/VDUsQ9tGVe3g1cgKkJKWcbhy8SK88go89xxcvKjcPOQos656jYY/fAinT0PFijB8uMlErr3W2KRk5Pwes9EwKhzOOlY5X9XZDGc1HgeVs9ySY7FYrlByPVTlqMEBcCNQAfgcs99iAmbT3d+8FGuBU5KFI5XTp+Hpp+HNN8HHB+6/N4WHe6yk+vefw1dfwalTRkSuvdaszho2zL2IxO53rM6KgDPrTFuNDmY4q/E4qNKycG/MYrEUW/IiHMuyOJ+qav+CCs7blAbhSOXAAZg5E+bONfsJH30U/npHEuXXrTSrs7780hgyVqiQXkQqV858sthIRyYSAdGORXLV26eJSNUSs+LaYrF4gUL3qipOlCbhSGXTJnjkEWOT1agR/PvfMGUK+JBsaoikisiJE2ZJ7zXXGBEZPhyqVMl8wotRcOgrM5x1erVpqxaUJiLV2hXq/VkslqInP6uqqmEsR652NP0CPKmqMQUepZfIq3A8/fTT7N+/n8DAQIKCgggMDKRBgwYYm6ziwc8/w8MPw/r1EBQEzz5rEg0RTOXCVavSROT4cSMiQ4emiUhVN468l46kicipVYBC1QAXEQlyXMBisZRm8iMcX2JqcKTaqN+AqQg4xvNRxYu8Csf06dOJiIjAtZ5HtWrVCAwMdP6kCkrdunWLTFBUja/iP/8Je/eaIoXPPQfdurl0Sk6G1avTROToUWNzMmSIEZERI9yLSNwxU9UwKhxOrTB+WlVap4lI9fZWRCyWUkp+hGOTqoZm11acye9Q1alTp9i2bRvbtm1j69atzt/PnDnj7HPVVVe5FZTatQvPTyoxEd5/H554woxQjRkDzzwDbTJOVaSkpBeRI0eM4aKriFSrlvkCcSfg8AIjIieXGRGp3DJtdVaNDlZELJZSRH6E4zdghqqucjzuCbyoqt29EqkX8MYch6py4sQJp4i4ikpMTNooXu3atd0KylVXec8qPTYWXn4Znn/eOL3fcgs8/ripepuJlBRYs8aISEQEHD5sRGTw4DQRqV4983Hxp+DwQiMiJ34CTYZKzRw71sdBzc5WRCyWEk5+hKM98DHGowrgLDBVVbcUeJReojAnx1WVo0ePphOU1J8LF9LcWurVq+dWUKq5+6afR06eNEt4334bfH3hgQfgoYfcJxOAEZHff08TkUOHTBGqQYOMiIwcaeqLZORytENEIuDEUkhJhIqN00SkVlfjp2WxWEoUeRIOR0Gl51T1QYdXFap63ntheofisKpKVTl06JBbQbnksvu7YcOG6QQlMDCQdu3aUdXd/EMO2bcPHnsMPv/cVLh99FGYPj2zk3uGgGHtWiMi4eEQFWVEZODANBFxlzUlnIXD35hM5PiPkJIAFRoay5PG46F2DysiFksJIT8ZxxpV7ZZlp2JOcRAOT6SkpHDw4MFMYrJ9+3bi4+Od/Ro3buxWUCq52+TngY0bzRLeH3+EJk3MEt6//MVsKMwSVVi3Lk1EDh40KcyAAUZERo0yipSRhBg48q1x8T36P0i5bJx7/cc4RKQXlMnu4haLpajIj3C8DTQEwoGLqe2q+lVBB+ktirNweCI5OZkDBw5kEpSdO3dy+fJlZ79mzZplEpSAgAAquDNBdLB0qVnCu3EjhISYJbxDh+ZwSkLVrP1NFZHISCMi/funiUitWpmPS7wARxY7ROQ7SI6H8nXSRKTO1VDGN9evk8Vi8R75EQ53hoaqqjcXVHDepiQKhyeSkpLYt29fJkHZtWsXiYmJAIgILVq0yCQobdu2pZxjfColxXzu//OfsH8/9O1rJtNnzDDXWb48B8GoGvVJFZH9+0360r+/8c4aPdoYMmYkMRaOfW+Gs44shuRLUK6WKUjVeBzU7WfqrlssliLF7hwvJcLhicTERPbu3ZtuufC2bdvYvXs3ycnJAPj4+NCyZct0YtK6dRC//NKGZ57x5dQpkyw0aWJGpnK1KEoV/vgjTUT27TMi0revyURGj4Y6dTIfl3QJjv3PISLfQlIslL3K1FdvPN6UyvUpm/k4i8XidfKTcTQHXgW6YUwOfwPuV9UD3gjUG1wJwuGJhIQEdu/enUlQ9u7dS0pKCgC+vr60aBGKj89D7NgxEtWyNG9+mb/8xZe//MWHgIBcXlQVNm9OE5E9e6BMmfQiUtdNqdukODOhHhUBRxZB4nnwqw7+I00mUm8Q+GQ1o2+xWAqSfE2OA29i3HEBJgL3qGrXAo/SS1zJwuGJ+Ph4du3alWlT475954CxmD9zH6AM1asfpFOnfYwZk8iAAc1o0aIFPtnOqDtQhS1b0kRk924jIldfbURkzBioVy/zccmX4fgSIyKHF0BiDPhVhYYjjIjUHwI+5TMfZ7FYCoz8CMcWVQ3J0LZZVdsXcIxewwpHzund+xKXLu3k/vu3sWZNFMuX12Lfvk5cvtzR0WMtvr5f0br1JsLCaqXbh9K0aVPKlMliqa0q/Pmn2SMSHg47d5rxMFcRqe+mXnpygtlkmCoiCWfAtzI0vM4hIteAr+fFABaLJW/kRziew2z6m0daPY4awAsAqnrG89HFAysc+Wf79ou8884ZvvmmApGRZtVU2bLrSEj4CIgATlCxYkUCAgIybWxs3LhxZh8vVdi2LS0T2bHDiEivXkZExo51v9U9JRFOLHOIyNdw+TT4VoIGw4yINLjWPLZYLPkmP8KR1VyGqmqxLx9nhaNg2bsX5s+HefNg61YoU0Zp3foYDRuuIilpPnv2rOHo0aPO/pUrV6Zdu3aZBKVhw4ZpguIqItu3GxHp2TNNRBo2zBxIShKc/MUhIl9B/EnwqWDEo9E4aDgM/NxYyFsslhxhV1VZ4fAK27aliciePWZLx6BBMHz4RVq0+JODB7ekm0M5ceKE89iqVau6tV2pd+YM8uWXRkS2bjWdXUXE3z9zICnJxgI+KtwUp4o/buZA6g81IuJ/nZkjsVgsOcYKhxUOr6JqikulisjBg8Yr8ZprYOJEuO46U8k2OjrardPw6dOnneeqUaOGU0x61apFj6NHabxmDb47dpgO3bsbERk3zlSxyhRMCpxaDYccJXLjjkCZsmZCvdE48B8BZasXyutisZRkrHBY4Sg0Um2u5s2DL74wpT8qVDDiMWGCEZOMG9tPnjyZacnwtm3bOHv2rLNP1+rVubV6dYZeuIB/dDQAiZ064TdpkhGRxo3dBJMCp39PE5FLUWZzYd2BZp+I/0go5z2nYoulJGOFwwpHkZCSYooQzptnFlOdOmUq144caTKRQYNMZuIOVeX48eNuM5S6588zDhgPhDn6R9aty7FevfCbNImWAwZQPaMdvCpEr3OISDhcjATxNZsMG483mw7Lu7FLsViuUPIzOS7AZKC5qj4pIo2Beqq61juhFjxWOIoHSUmwbJkZzvrySzh3zri0jxljRKRvXzNHkh2qypEjR5wicnL1ahr//js9jh4l1LGpcS3wQ9Wq7A4JoVanTs45FKfTsCqc3WiykKhwiN0H4gN1+prVWf6joYKbTYoWyxVEfk0OU4D+qhogIjWAH1W1s3dCLXiscBQ/EhJgyRKTiSxYYIpP1a5tpi4mTDCrcrPaEuIOVeXoypVcmD2b6kuWUO/IEQA2lCnD/JQUwoFIoFGjRul9vNq1o12DFCqfWWxE5MJuY/1e+2pHTZExxtXXYrnCyI9wbFTVDiLyh6qGOdrsBkBLgREXB99/b0Tk22/N44YNjYhMnAhduuSxmOD+/RARgYaHI46//7GGDfmlTh3mXr7M0v3701nXN23a1AhJizoE1jlLUOUttK26n4rlBGr3NBPrjcdCRTeruiyWUkh+hON3oAewziEgtTEZR1iWBxYjrHCUHGJj4ZtvzHDW99+bzKRpU5OFTJgAoaF5FJHIyLQd62vNKKuGhRHdrx9/tGzJ2jNn0lnXJyQkAMZpuLl/DQIbJBFY9zyB/hAYFELb7n+hfKuJUKlJQd26xVLsyI9wTMbsFu8AfASMAx5T1S+8Eag3sMJRMjl3DhYuNJnIkiWQnAytWxsBmTgR2rXL44kPHkwTkd9/N22hoSbFGT+epGbN2Lt3r1vr+qSkJADKCLSsB4HNqhMY3JHAbtcR2LE/bdq0oayn2X6LpYSRr1VVItIWGAAI8JOq7ij4EL2HFY6Sz+nT8NVXJhNZtszMbQcFGQGZMAFatszjiaOi0kRkzRrT1r592j6RNm2cXRMSEtizZ48Rko0r2bZxOdt27mXPkXiSzZw8Pj5laN2yGYHBYenmUVq1aoWfn60xYilZ5Cfj+ERVb8iurThjhaN0ceyY+ayfPx9+/dW0dexoBOT66009kTxx6JBZ7hUeDqtXm7bgYGcmQtu2bg+7HL2LXSvfY9uahWzbuZdtR2Db0bLsPZZI6v8vPz8/2rRpk6m4VosWLfDNyVIyi6UIyPfkuMtjH+BPVc3rQEGhY4Wj9BIVZT7n580zFW0BevQwIjJ+vHuz3Rxx+HCaiKSqU1BQWibiaZzs4iE49BUciiDuyCp2HoVtZxqy7WxTth31Y+uugxw4kGb/Vq5cOdq0aeO0W0n9adasWc6t6y0WL5Fr4RCRR4B/AhWAS5hhKoAE4F1VfcRLsRY4VjiuDPbtM1nI/PmmBIgI9OljhrPGjnVfCj1HHDlixsnCw81uRlUjHKmZSGCg++MuHXWKCCdXAApV23Kx1gh2XGrPtqgEtm3f7tzcGBUV5Ty0fPnyTqdhV1Fp0qRJ1tb1FksBkp+M4z8lSSTcYYXjymPHjjTfrF27TBXbgQNNJjJ6NGTcVJ5jjh5NE5GVK42IBASkFxF3y77ijptaIlERcHKZsUKp3NLsE2k8HmqEcSE2lu0OIXHdLX/EsR8FoFKlSm4FpVGjRpmt6y2WfJLfyfEaQCvAWXJNVVcUaIRexArHlUtqFdtUEYmMNBYnQ4aYTGTECBg+3PRdvjyXJz92DL7+2ojIihXGX6VtWzOUNX68mR9x92EefypNRE78BJoMlZo5NhuOg5qd0x137tw5t4Jy/PhxZ58qVaqks65PFZUGDRpYQbHkmfxkHLcC9wH+wCZM7fHfVLW/F+L0ClY4LGBEZN26NPPFI0egfHmoXNnMhWzenMc9IgAnTqRlIr/8YkSkdeu0TCQkxP3JL0fD4UVmx/qJpaZQVcXGaSJSq6vZxe6GMy57T1wF5dSpU84+1apVc2tdX7duXSsolmzJj3D8CXQG1qhqqGNp7jOqOsY7oRY8VjgsGUlJMXPe8+fDu+9CYqKZ+/7732HSJChXLh8nP3kyLRNZtsxcrFWrtEzE0y7GhLNw+BszJ3LsB0hJgAoNodFYIyS1e3oUEVdOnTqVSUy2bdvGmTNpxTqvuuoqt4JSu3btfNy4pbSRH+FYp6qdRWQT0FVVL4vINlX1MCNY/LDCYcmKq682n/Vly5qS6PXqwT33wJ13wlX5dVw/dSq9iCQnQ4sWaZlIWJh7EUk8nyYiR7+HlMvGL8t/jENEekOZnK+6UlVOnDjh1mk4JibG2a927dqZlgwHBgZSs2bNfL4QlpJIfoTja+Am4H6gP6b+uJ+qXuuFOL2CFQ5LVvTta/5dtszsUH/pJfjxR6hYEW6+Ge6/33zW55vTp42IRETATz8ZEWnePC0T6djRg4hcgKPfmTmRo4shOQ7K1zEOvo3HQ50+UCZve0FUlaNHj2bKULZv386FCxec/erVq+dWUDJZ11tKFQVSj0NE+gDVgO9VNbEA4/MqVjgsuWXLFvi//4PPPjN28KNHw4MPmuKDBUJ0tLEFDg83IpKUBM2apYlIp07uRSTposlAoiLg6LfmcbmaRkQajYN6/U2hqnyiqhw6dMitoFy6dMnZr2HDhpnExGldbynx2J3jVjgseeDoUXjjDXj7beOd1b27EZCRI80S3wLhzJk0EVm61IhIkyZpIuLJHjgpzsyFHIowE+xJF6BsDVOQqtE4qDcQfArWNyslJYWDBw9mmpTfvn17Oqfhxo0buxWUSpUqFWg8Fu9id45b4bDkg9hYmD0bXn4ZDhwwQ1f33w833WRqqRcYZ88aZ8fwcDNulphoSuKmikjXru5FJDkeji1xiMhCSIwBv2rQcISZE6k/GHzKZz6ugEhOTubAgQOZBGXHjh1Op2Ew1vUZd8kHBARQIWMtYUuxoCB2joPZPW53jluuWJKTzTTFSy8ZT8QaNeCuu+Duu/Nhb+KJs2dh0SIjIj/+aESkUSOzDX78eOjWzX21q+QEOL7UISILzGot3yrQ8DqHiAwF38L5oE5KSmLfvn1unYYTE81ot4jQvHnzTILSpk0bypf3nthZssfuHLfCYSlgVq+GF180o0x+fjB5Mvztb2ZZb4Fz7pwpVBIeDj/8YAqV+PuniUj37u5FJCURTiwzcyKHv4bLp8G3EjQYZkSkwbXmcSGTmJiY5jTs8rN7926Sk5MBKFOmDC1btsy0S75169bWur6QyI9w9AQ2qepFEZmCqcvxqqoe9E6oBY8VDos32bvXDGHNnm2qFw4ZYuZBBgzIx4bCrIiJSS8ily+bkoljx5ohrZ49PYhIEpz8xSEiX0H8SfCpYMSj0ThoOAz8qngh4JyTkJDA7t270y0X3rZtG3v37iXFUU/e19eXVq1aZRKUli1bWuv6AiY/wrEFaA+EAHOA94HrVbWPF+L0ClY4LIVBdDS88w68/rrZSB4SYjYUTpxo9oh4hfPnTb3d8HBTMvHyZTNmlpqJ9OzpfhY/JRlOrTLDWYe+hLhjUKYcNBjqEJHroGw1LwWde+Lj49m1a1cmQdm/f79H6/pUUWnRooV1Gs4jBVFzfCZwRFU/yDhhXtyxwmEpTC5fhrlzzTzI9u3QoAHcey/ccUc+zBVzwoUL6UUkPt7sZkwVkV693IuIpsDp30wmcigCLh02S3rrDTbDWQ1HQLn87oT0DpcuXWLnzp2ZNjVGRkY6+5QrV462bdtmEpRmzZpZp+FsyI9w/AL8D7MJ8GrgJLBZVYO9Eag3sMJhKQpUzUjSiy+arRqVKsGtt5rVWE2bevnisbGweLERke++M2NodevCmDFGRK6+2rOIRK8z3lmHIuDiQRBfqDfAZCL+o6B8Xv3pC4/Y2Fh27NiRaaf8oUOHnH0qVKjg1mm4cePGVlAc5Ec46gF/Adap6koRaQz0VdWPvRNqwWOFw1LUbNpkMpB584x11bhxZhirS5dCuHhsrBGP8HAjJnFxUKdOehFxV4VQFc5sMAISFQGx+0B8oE5fk4n4j4YKdQvhBgqO8+fPO52GXTOUo0ePOvtUqlTJ6TTsKij+/v5XnDFkgewcL6lY4bAUFw4fNnMg//2vmePu1ctMpF93nfv57ALn4sX0InLpEtSubbbGjx9v/Fc8ici5zUZAosLhwm5juFi7t8lEGo2Big0K4Qa8w9mzZ9NZ16eKyokTJ5x9qlat6lZQ6tevX2oFxQqHFQ5LMeLCBfjgA3jlFTh40Jjn/u1vcOONxiOrULh0ycyFhIebuZGLF02ZxFQR6dfPs4jEbEubE4nZBohx700VkUqNCukmvEt0dLRbp+HTp087+1SvXj3T/ElgYCB16tQp8YJihcMKh6UYkpRkyni8+KKpFVKzJkyfDn/9q5mSKDQuXYL//c+IyDffGBGpWdOIyLhx0L+/2azijpgdZmVWVDic22LaanZz1BQZC5WbFtptFBYnT550Kyhnz5519qlZs2YmMQkMDKRWnmsYFz5WOKxwWIoxqqac+Usvmc3iZcvCDTeYLCQgoJCDiYszIhIRYYKJjTX+8qNGmUxkwADPInJ+t0NEIuDsRtN2Vac0EanSstBuo7BRVY4fP+62uNb58+ed/erUqeNWUGrUqFGE0bsnL5YjfwIeVUVVQwouPO9ihcNSkti922wonDPHrKi99lozD9K3r5c2FGZFfLxZGhYebkTkwgXjszJqlMlEBg70vEkldj9EfWmGs6LXmrYaoWY4q/E4qNqmsO6iSFFVjhw5kklQtm3bRmxsrLNf/fr1MwlKu3btqFat6PbT5EU4mjh+/avj308c/04GUNV/FHSQ3sIKh6UkcuqUceV94w3ze1iYEZDx4z1/4fcq8fHGMysiwhgxnj9vNqaMHGmCGjTIs4hcPAiHvjKZyOnVpq1aUFqJ3Oolpi5cgaGqREVFuXUadrWu9/f3dysolStX9nqM+VmO+4eqhmVosxsALZZCIj4ePv3UDGPt3Gksqu67D267DYrsy+jly8a9NzzciEhMjAnGVUQ81d+9dMSIyKEIOLkSUKga4CIiwUWQWhUfUlJSiIyMdOs07Gpd36RJk0yCEhAQQMUCXF2RH+HYBPxVVX91PO4BvKWqoQUWnZexwmEpDaSkmEVQL71kqhVWqWLE4777jPN6kXH5sqkjkioi585B1aowYoQRkcGDwZPLbdwxOPS1Q0R+MRsQq7RKG86q4aG07hVIcnIy+/fvzyQoO3fudFrXiwjNmjVLJybDhg3Lc6XG/AhHR+BDTOU/wZSOvVlVN+YpkiLACoeltLFxoxGQ+fPN4+uvNxsKO3Ys2rhISDAiEhFhbIPPnjUKlyoiQ4Z4FpH4k8YGPioCTvwMmgyVmqVlIjU7WxFxQ1JSEnv37nVrXZ+UlMSuXbto3bp1ns6d71VVIlINQFVjsutb3LDCYSmtREXBa6/Bu++aees+fYyADBtWSBsKsyIhAX7+2WQiCxaYSodVqpjdjuPGwdCh4KmA0+VoU5AqKgJOLDX28BUbm5VZjcdBrW5mA6LFI6nW9W3atMmzyWN+Mo5ywFigKeDcDaSqT+bgokOBVwEf4H1VfTbD800w2Uxt4AwwRVUPO557Dhjm6PpvVZ3vaG8GzANqAhuAG1Q1gSywwmEp7Zw/D++/bzYUHjoEbdqYpbw33OD5s7lQSUw0IhIRYSphRUdD5cowfLjJRK65xnOgCWfh8DdmOOvYD5CSABUauIhITyhj3W+9QX6E439ADOZDOjm1XVVfyuY4H2A3MAg4DKwDJqnqdpc+4cC3qvqRiPQHblLVG0RkGHA/cA1QDlgODFDV8yLyBfCVqs4TkXcwhotvZxWLFQ7LlUJiovlsfvFFM5xVu7bZTDh9uvm9WJCYCMuXm0zkq6+MiFSqZERk3Diz/tjTBG/ieTjyrclEjn1vSuaWr2t2qzcebyxQyrjZ7W7JE/kRjq2qmuuaZiLSHZilqkMcjx8BUNX/uPTZBgxV1UNi9ubHqGpVEZkBlFfVfzv6fQD8AIQDp4B6qpqU8RqesMJhudJQhV9+MfMg335rphWmToUHHjDZSN++pt/y5UUZJWbrvKuInD5tRGPYMJOJXHut56LuiRfg6HdGRI4uhuQ4KFcbGo02cyJ1+xp7eEue8SQcORkkXC0iebFQbwgccnl82NHmymZgjOP30UAVEanpaB8qIhVFpBbQD2iEGZ46p6pJWZwTABG5XUTWi8j6U6dO5SF8i6XkImLE4ZtvTE2QG24wGwrbtjXz1OfOGXEpcnx9zSbC//4Xjh0z/vM33mhU7/rrjYvv+PHwxRfGBsUVvyrQZAL0Doexp6BXhLF/j/wMlg2Gr+rBmlvg6PemDrulwMhJxrEdaAkcAC5jVlZpdjvHRWQcJpu41fH4BqCrqt7t0qcB8AbQDFiBmUsJUtVzIvIoMB6TYZzEDHV9CqxR1ZaO4xsB32eXEdmMw2KBkyfhzTfhrbfMF/ty5eCmm8yX+/79C9FcMSckJ8OKFWmZyIkTZg7k2mvNcNbw4WaOxB1JcXD8R5OJHFlkhrf8qoH/SJOJ1B8EPh5WdlnSkZ+hqibu2rOrOZ6ToaoM/SsDO1XV381zn2FE43vsUJXFki/i4swu9NOnzebCixfNUFa/fkZEhg0rhEJTuSE5GVauNCLy5ZdGRMqXNxPq48cbEanioVZ68mU4vtRRIncBJJ4D3yqmNG7jcVB/KPgWh9UDxZP8CIfbrUWqGpXNcb6YyfEBwBFMxvAXVd3m0qcWcEZVU0TkaSBZVWc6Jtarq2q0iIQAnwGhDrEIB750mRzfoqpvZRWLFQ6LJT2pcxw//GC+2C9ebH727jXt7dqliUiPHkVkceKO5GT49dc0ETl2zKROriJStaqHYxPgxDIjIoe/Nkt+fStBg2FGRBpcax5bnORHOFLNDgUojxlW2qWq2ZrLiMi1wCuY5bgfqurTIvIksF5VFzmGs/7jOP8KzA71yyJSHkjdYHgeuFNVNznO2RyzHPcq4A/MEt7LWcVhhcNiSY+nyfHdu81k+uLFRlCSkowd1ZAhRkSuucaU7CgWpKSkF5GjR42IDBliROS66zx7sqQkmZ3qURFw+Cuz+dCnAjS4xgxnNRwGfh4E6AqiwGzVRaQDMD117qIkYIXDYsk9588bO6rFi03RwBMnzKR7165p2UhoaDHZzJ2SAr/9ZkQkIgKOHDGGi64i4sl2IyUZTq1yDGd9aWxQypSD+kNMJtLwOijr4dhSToHW4xCRP1U1LyutigQrHBZL/khJMftCUoe01q0z7Q0amPnqYcPM4qhCMGzNnpQUWLMmTUQOHzZjbYMHGxEZOdKziGgKnP4trbrhpcNmSW+9QWbDof9IKFezUG+nKMnPUNXfXB6WAToANbObkC5OWOGwWAqW48eN4eLixcZp/cIF8wW/b9+0bKRFi6KOEiMiv/9uBCQ83Gyr9/Mz7r2pIuKpgJKmQPQ6IyBREXAxEsQH6vY3mYj/KChfpzDvptDJj3A87vIwCYjETE7Huz+i+GGFw2LxHgkJpnphajaya5dpb9MmTUR69fJcqqPQUIW1a9MykYMHjYgMHGiW+I4aZSodejr27EYjIFERELvXeGXV6eOosz4aKtQv1NspDArC5LAygKrGZte3uGGFw2IpPPbtSxOR5cuNsFSpkn6CvVDrqbtD1Yy3pYpIZKTZjDhggMlERo0yNdc9HXtuS9pw1vmdgEDtXg4n3zFQMdOughJJfjKOIEz1v1QpPg1MVdWtBR6ll7DCYbEUDbGxZjP4t9+aCfajR017585p2UiHDkXs5KsKGzYYEQkPhwMHwMcnvYhktZTs3La04awYx8dire6OmiJjoZLbrXAlgvwIx2rgUVVd5njcF3hGVXt4IU6vYIXDYil6VGHTprRs5PffTVu9eiYLGTbMTD142oZRaEFu3JgmIvv3GxHp18+IyOjRWbtFnt9lVmZFRcDZP0zbVZ0dmchYqFIcJn5yTn6EY7Oqts+urThjhcNiKX6cOgX/+58RkR9+MP5Zfn7Qu7fZxzdsGOSx/lDBkKp0qSKyd68Rkb5900SkThaT4xf2pYnIGccytBphaYWpqhblzeWM/AjH15jNeJ84mqYAHVV1dIFH6SWscFgsxZukJFi9Om3z4XZH8YWWLdOGtK6+2nMZc6+jCps3p4nInj1mfK1PHyMiY8ZkPXETG5lWZ/30b6atenBaidxq7QrlNnJLfoSjBvAE0Auzw3sl8ISqnvVGoN7ACofFUrKIjEwb0vr5Z1PWvHJlswBq2DCzd6RBgyIKThW2bEkTkd27jYhcfXWaiNSr5/n4S4eNiERFmI2HKFQNSMtEqgcXk12VeRQOh2fUUlXt583gvI0VDoul5HLpkhGPVCE55CjWEBaWNqTVuXMRTbCrwtataSKyc6f50O/d24jI2LFQP4tlunHH4NDXJhM5+YvZO1KlVVomUiOsSEUkPxnHT8CYklhrPBUrHBZL6SD1czp1SOu338wev9q10ybYBw92vzHc68WrVGHbtjQR2bHDfOj36pUmIlmlSfEn4fACR531n0GToVKztEykZudCF5H8CMdCIAxYAjgrqajqvQUdpLewwmGxlE6io83E+uLFZqL9zBkzf92rV9rcSEBAWmErKMSqh9u3p4nItm0miJ49zWbDceOgodsadIbL0XB4oUNElkJKIlRslLbEt1Z3swHRy+RHOKa6a1fVjwooNq9jhcNiKf0kJZklvqlDWlu2mPamTY2ArFhh3EV++aUIgtuxI01Etjr2evToYTKRcePAP4sNgwln4fA3Zjjr2A+QkgAVGpiNho3GmY2HZXy8EnaBmhyWNKxwWCxXHlFRZtPh4sVmE2JcnLE9mTEDbr21CItV7dyZ5p2Vqm7du6eJSKNGno9NPA9HvjWZyLHvITne+GX5jzFDWnX6QBnfAgs118IhIiMBf1V90/H4dyB158tDqhpRYNF5GSscFsuVTVycmUA/diyt3vrgwXDbbaYGe5EVqtq9Oy0T2bzZtHXtmiYiTbLYdZ4YC0e/M5nIkcWQfMk49/qPNplIvf7G2Tcf5EU4fgUmquohx+NNmGp+lYDZqjogXxEVIlY4LBZL6hzHJ5/Ahx/C++8bx/W6dU3t9VtvLWJH3z170jKRPxy7zrt0SRORrFKkpEtmGCsqAo58A0kXoGwNYwMfPCvPtieehCOr2ZWyqaLhYJWqRjtKxtr6ihaLpUTSqBE8/rjZK/Ltt9CtG7zwgtlsOHAgfPGF2TdS6LRqBY88YixP9uyB//zHTNzMmAHNmpmU6fnnjZdWRnwrGofennNh7Em4epEpQHV4EZQpX+ChZpVx7FXVlh6e26eqJcZ0xWYcFoslK44cgdmzTRZy8KDxNJw2zQxlFantCRi74dRMZMMG09axo8lExo+H5s09H5uSlK85j7wMVc0Flqvqexna7wD6quqkPEdTyFjhsFgsOSE5GZYuhXffhUWLzBf+vn2NgIwZA+UL/st77jhwIE1EUsswduiQNpzV0u13/TyTF+GoAywALmO8qgA6AuWAUap6okAj9CJWOCwWS245fhzmzIH33jMmuVddBVOnGhEJCCjq6DBjbakisnataQsNTctEWrXK9yXys4+jPxDoeLhNVX/OdzSFjBUOi8WSV1JSjOXJe+/B119DYqLZYHj77eZLfoUKRR0hZnztyy+NiKxZY9ratzcCctttWbv4ZoHdx2GFw2Kx5JOTJ+Hjj81Q1p49xtrkhhvMZ3NwcFFH5yAqKr2I7NtnJtfzgBUOKxwWi6WAUDU70N9913xGJySY1Vm33w7XXw+Visu60xMn8lWnNy/LcS0Wi8XihlTvq88+Myuy/u//ICYGbr7Z+BhOn25qQBU5XiruboXDYrFY8kGtWvDAA8bHcOVKGDnSLO0NCzP7995/39ReL01Y4bBYLJYCINVB/eOP4ehReO01Y3Vy222mJMcdd6RtwyjpWOGwWCyWAqZGDbjnHuNh+NtvZnHTJ59Ap05m794778D580UdZd6xwmGxWCxeQsRMmn/4oTFYfPNNs8nwrrtMFnLrrcYKvqStUbLCYbFYLIVAtWpm0vyPP8x+vb/8BebNM8ISGmpE5dy5oo4yZ1jhsFgslkJExPgVvveeyUL++19j63733WZF1rRpsHp18c5CrHBYLBZLEVGlitn7sX69mTifOhW++spUmA0KgldfNeVwixtWOCwWi6UY0KEDvP22WZH1wQdGVO6/32QhU6aY0rfFJQuxwmGxWCzFiMqVzUbCNWtMUcDbbjN1Q/r0MeaKL70Ep08XbYxWOCwWi6WYEhICr79uspA5c8xmwwcfhIYNYdIkWLasaLIQKxwWi8VSzKlY0cx/rFoFW7ea5bw//AD9+5tCU88/bwwYXenbN61cbkFjhcNisVhKEIGB8MorJgv59FMzB/Lww+DvbzYaLllirOC9iRUOi8ViKYGULw+TJxuX3h074N57zdDV4MGmEGBUlPdqp1vhsFgslhJO27bw4ovGqXfePFN+48ABM8HuDX8sKxwWi8VSSihXDiZMgJ9+MpsMmzQxu9ILGiscFovFUgqpWBGaNgUfn4I/txUOi8ViseQK36IOwGKxWCwFz/Ll3ju3zTgsFovFkiuscFgsFoslV1jhsFgsFkuusMJhsVgsllxhhcNisVgsucIKh8VisVhyhRUOi8ViseQKKxwWi8ViyRVWOCwWi8WSK6xwWCwWiyVXWOGwWCwWS66wwmGxWCyWXGGFw2KxWCy5wgqHxWKxWHKFFQ6LxWKx5AorHBaLxWLJFVY4LBaLxZIrrHBYLBaLJVdY4bBYLBZLrrDCYbFYLJZcYYXDYrFYLLnCCofFYrFYcoVXhUNEhorILhHZKyL/cPN8ExH5SUS2iMhyEfF3ee55EdkmIjtE5DUREUf7csc5Nzl+6njzHiwWi8WSHq8Jh4j4AG8C1wDtgEki0i5DtxeBj1U1BHgS+I/j2B5ATyAECAI6A31cjpusqqGOn5PeugeLxWKxZMabGUcXYK+q7lfVBGAeMDJDn3bAz47fl7k8r0B5oCxQDvADTngxVovFYrHkEF8vnrshcMjl8WGga4Y+m4ExwKvAaKCKiNRU1d9EZBlwDBDgDVXd4XLcbBFJBr4EnlJVzXhxEbkduN3xMFZEduXxPmoBp/N4bEnF3vOVgb3n0k9+77eJu0ZvCkdOeBB4Q0SmASuAI0CyiLQEAoDUOY8lItJbVVdihqmOiEgVjHDcAHyc8cSq+i7wbn4DFJH1qtopv+cpSdh7vjKw91z68db9enOo6gjQyOWxv6PNiaoeVdUxqhoGPOpoO4fJPtaoaqyqxgLfA90dzx9x/HsB+AwzJGaxWCyWQsKbwrEOaCUizUSkLDARWOTaQURqiUhqDI8AHzp+jwL6iIiviPhhJsZ3OB7XchzrBwwHtnrxHiwWi8WSAa8Jh6omAXcDPwA7gC9UdZuIPCkiIxzd+gK7RGQ3UBd42tEeAewD/sTMg2xW1W8wE+U/iMgWYBMmg3nPW/fgIN/DXSUQe89XBvaeSz9euV9xM69ssVgsFotH7M5xi8ViseQKKxwWi8ViyRVWOLIgO8uU0oiI+IjIHyLybVHHUhiIyAMOa5utIvK5iJQv6pgKGhH5UEROishWl7YXRGSnw+7naxGpXoQhFjju7tnRfo/jvreJyPNFFZ83EJFGIrJMRLY77u8+R/tVIrJERPY4/q2R32tZ4fBADi1TSiP3YRYzlHpEpCFwL9BJVYMAH8zqv9LGHGBohrYlQJDD7mc3ZlVjaWIOGe5ZRPph3Cnaq2ogxvKoNJEE/F1V2wHdgL86PrP+Afykqq2AnxyP84UVDs/kxDKlVOEwmRwGvF/UsRQivkAFEfEFKgJHizieAkdVVwBnMrT96Fj5CLCGtM22pQJ39wzcBTyrqpcdfUqVz52qHlPVjY7fL2C+ADbEfG595Oj2ETAqv9eywuEZd5YpDYsolsLiFeAhIKWI4ygUHJtJX8TsGzoGxKjqj0UbVZFwM2aTbWmnNdBbRH4XkV9EpHNRB+QtRKQpEAb8DtRV1WOOp45jtj7kCyscFgBEZDhwUlU3FHUshYVjrHck0AxoAFQSkSlFG1XhIiKPYoY45hZ1LIWAL3AVZhhnBvBFarmG0oSIVMbYMd2vquddn3P4+uV7D4YVDs9ka5lSyugJjBCRSMywXH8R+bRoQ/I6A4EDqnpKVROBr4AeRRxToeHwiBuO8X+7EjZ0HQa+UsNaTGZdq4hjKlAcjhpfAnNV9StH8wkRqe94vj6Q7yE6KxyeydYypTShqo+oqr+qNsXc68+qWtq/fUcB3USkouOb5wCunIUBQzHDkiNU9VJRx1NILAD6AYhIa0zZhlLjlOt4D38A7FDV/3N5ahEw1fH7VGBhfq9V1O64xRZVTRKRVMsUH+BDVd1WxGFZChBV/V1EIoCNmOGaPyiFlhQi8jnG3qeWiBwGHsesoiqHcZ4GYyp6Z5EFWcB4uOcPgQ8dS3QTgKmlLNPqiXEL/1NENjna/gk8ixmWuwU4CFyf3wtZyxGLxWKx5Ao7VGWxWCyWXGGFw2KxWCy5wgqHxWKxWHKFFQ6LxWKx5AorHBaLxWLJFVY4LMUGEVERecnl8YMiMquAzj1HRMYVxLmyuc54EdkhIsvcPNdaRL5zuJRuFJEvRKSuiPRNdSMWkRGpTswiUtthj/GHiPTO6ty5iG+5iHRy0+4rIs84Ytvk+Hk0B+f7Z15jsZRcrHBYihOXgTGpdeWLCw4DxJxyC3CbqvbLcI7ywGLgbVVtpaodgLeA2q79VHWRqj7reDgA+FNVw1R1padzF1DcT2FsV4JVNRToDfjl4DgrHFcgVjgsxYkkzAa8BzI+kTFjEJFYx799HYZ1C0Vkv4g8KyKTRWStiPwpIi1cTjNQRNaLyG6HN1dq/ZEXRGSdozbFHS7nXSkii4DtbuKZ5Dj/VhF5ztE2E+gFfCAiL2Q45C/Ab6r6TWqDqi5X1Yz1IqaJyBsiEgo8D4x0fPt/3PXcIlJeRGY7YvjDYRmeevwiEfkZ+ElEKojIPEem8jVQwc29VARuA+5R1XhHbBdUdZZLnwUiskFMnYfbHW3PYpyFN4nIXEfbFMdrv0lE/ut4fX0cf7+tjngz/X0tJQu7c9xS3HgT2CK5K7LTHgjA2GjvB95X1S5iCtncA9zv6NcUY5ffAlgmIi2BGzGuuJ1FpBzwq4ikOuR2wNSsOOB6MRFpADwHdATOAj+KyChVfVJE+gMPqur6DDEGATk2kFTVTQ4h6qSqdzuu2y/13CLyd9NNg0WkrSOG1i5xh6jqGRH5G3BJVQNEJASzSz4jLYEohxW3J252nK8CsE5EvlTVf4jI3Y4MBREJACYAPVU1UUTeAiYD24CGjponSCkrGnUlYjMOS7HC4eb5MabAUk5Z56hFcBnYB6R+8P+JEYtUvlDVFFXdgxGYtsBg4EYxFg2/AzWBVo7+azOKhoPOwHKHOWKqs+zVuYi3IOgFfAqgqjsxVhKpwrFEVVNrUVzt0m8LsCW7E4vITY6M4ZCIpBp93isimzG1OxqR9hq5MgAjpuscr+cAoDnmtW4uIq+L8cg67+ZYSwnCZhyW4sgrmG/Gs13aknB80RGRMhiDulQuu/ye4vI4hfTv8Yz+OgoIZojmB9cnRKQvcDEvwXtgG9CnAM+XFbmNey/QWESqOIaoZgOzxXg6+Thei4FAd1W9JCLLAXcldgX4SFUzVRMUkfbAEOBOjFfSzbmM0VKMsBmHpdjh+Lb8BWYyOJVIzLdZgBHkbOI2I+NFpIxj3qM5sAtjYnmXGDvq1JVPlbI5z1qgj4jUElNieBLwSzbHfAb0EJFhqQ0icrWIBOXhPgBWYoaBUp1eG2PuJyMrMPMrOK4VkrGDwx33A+ANxyR+aunkVHGuBpx1iEZbTD2LVBJTXztMWdJxIlLHcY6rRKSJY7FDGVX9EvgXZijNUoKxGYeluPIScLfL4/eAhY7hkv+Rt2wgCvOhXxW4U1XjReR9zHDWRhER4BTZlNZU1WNilswuw3zLXqyqWVpVq2qcY0L+FRF5BUjEDBvdR95qQrwFvC0if2KysWmqelky1yV6G5M97MBYxnuaZ3kU+DewVUQuAHGYMqNHMdUR73ScYxdmuCqVdzFzUhtVdbKI/Asz31LGcY9/dZxrtqMNSl998ysO645rsVgsllxhh6osFovFkiuscFgsFoslV1jhsFgsFkuusMJhsVgsllxhhcNisVgsucIKh8VisVhyhRUOi8ViseSK/wfF10Tgk96j9gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#graph constants\n",
    "YMIN = 0.985\n",
    "YMAX = 1.000\n",
    "YTICK = 0.005\n",
    "\n",
    "#graph\n",
    "title_text = \"Ground state probabilty against number of Clifford Gates\"\n",
    "#plt.scatter(number_of_gates_simulation, simulation_probability, label = \"Simulation results\" )\n",
    "#plt.plot(number_of_gates_simulation, mean_probability, label = \"Simulated probabilty\" , color=\"blue\")\n",
    "plt.errorbar(number_of_gates_simulation, mean_probability, yerr = standard_error_list,\n",
    "              label = \"Simulated probability\" , color=\"blue\")\n",
    "plt.plot(number_of_gates_calculation, theoretical_probability, label = \"Calculation based on error per Clifford gate\" , color=\"orange\")\n",
    "plt.plot(number_of_gates_calculation, egan_one_ion_noise, label = \"Egan's upper bound for 1 ion chain noise\" , color=\"black\")\n",
    "plt.plot(number_of_gates_calculation, egan_fifteen_ion_noise, label = \"Egan's lower bound for 15 ion chain noise\" , color=\"red\")\n",
    "plt.title('Single qubit errors')\n",
    "plt.xlabel('Number of Clifford Gates')\n",
    "plt.ylabel('Ground state probability')  \n",
    "plt.legend()\n",
    "plt.ylim(YMIN, YMAX)\n",
    "plt.xticks(np.arange(0, 21, 4))\n",
    "plt.yticks(np.arange(YMIN, YMAX, YTICK))\n",
    "fname = \"Single_Clifford_gate_errors_graph.png\"\n",
    "plt.savefig(fname)"
   ]
  },
  {
   "attachments": {
    "Egan%20Figure%20S1.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMMAAACDCAYAAAAwJ6/GAAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAAAJcEhZcwAADsMAAA7DAcdvqGQAABojSURBVHhe7Z1ZbFRXmsd5yMO0NJEijdSdoZWZkWYi8TJSq6NMIpFMEimdhxll6YlEMh1aZJRJp5sJISik2c1udgM2tjE2GIPB7FvYN0MAEwhg9sVgDMYG433DeP+mfsf3kMIp21W37r12lc9PunLVrXJVuXz+93zb+c4AMRgMCiMGg8HCiMFgsDBiMBgsHBdDa2ur3Lt3T/3U9wsKCtTt4uJiqaurU+fu3r2rzhkMfQXHxdDU1CRxcXFy69Ytdf/w4cOSlJQku3btkszMTNm+fbt899136tz169fVc86fPy8bNmyQixcvys2bNyUvL88c5rB9cNG1gytm0t69e+XKlSvS3t4u69atk3PnzsnEiRPl9OnTsmXLFlm6dKmcOHFC9u/fr54/Y8YMGThwoAwePFjmzJkjsbGx5jCHrWPmzJkyZMgQKS0tVWMrFFwxk+Lj4+XgwYPKFDp79qwkJyfLtWvXlDCOHDmiRME5/w9cWFioRGEwhAvjSJvmoeDKzIAgONra2tR9/ZOZQuN/G+7fvy9Tp0617hkM9pk2bZotn9QVMdgBp9uIweAERgwGg4URg8FgYcRgMFgYMRgMFkYMBoOFEYPBYGHEYDBYGDEYDBZGDAaDhRGDwWBhxGAwWBgxGAwWRgwGg4URg8FgYcRgMFgYMRgMFkYMBoOFEYPBYGHEYDBYeC4Gul1s3LhRampqrDPhYcRgcArPxUCzsE2bNsmHH34oc+fOlbKyMusRexgxGJzCczHcvn1bUlNT5fjx45KdnS3Hjh2zHrGHEYPBKTwVA32Qdu/erVpB0vzrxx9/tB6xjxGDwSk8FcPjx49Vu8iFCxfKsmXLJD8/33rEPkYMBqfwVAwMXJoF40DTO/XBgwfWI/YxYjA4hadioI8lZhKds9euXWvrjTtjxGBwCs/FsH79etm8ebP6id8QLkYMBqfwVAx37tyR3NxcZSJhKhkxGPoSnoqBDtpVVVXKeebgdrjc8wlq2syZ1j2DwT6eigHYfYfNRsgxYC6Fy4PSUon5ywhpe1hhnTEY7OGpGNh+asKECTJ//nx1sAuPpra2VjIyMiQnJ0fdP3PmjGRlZanyDTLWRKFwvjGxNMw0d0ruy4SPh0ndlGSpnJ0mjw79IO0Nj61nGAzB46kYGhsbpaKiQu3Tlp6e/lRolYz0gQMH1C49wAzCTj7Tp09XMwj7ve3YsUPS0tLk4cOH6jlTpkyRX/3yV/KXkV+p+43nr0vVkiwpG7dIquIzpfHMZXXeYAgGz80kBv0XX3yhBvKePXuss6Ky0uzZpvMPlG3ExMQ82c8tISFB7dLDzwsXLqjfYdZYtWqVjBs3Tt3XtDc2ScPRM1I5d4WUjV8sNSu2SlPeHevR8Hh0IEea79jbCM/Qt/FcDGxaSD3S4sWLf1aXxI6LVLNyUMCHIPR5CvyIHCEAf6qrq9Xs0RWtFdVStz1byqcmS3nMEqnbtF9aS8qtR0OnatFqeWxmnKjEczE0NzerKz1iaGhosM7aJ5TQKlf0mtXfSfnEeKmYuUxd5dvqQ/sM1cnrpTH3aUEaogPPxcD+ztj++Af4DuFiN8/QeOGGVCet6/AvuNqfvoRHbj3aNUYM0YunYsAvwL7H7mdmQBDhEm7Srb25RRq+PyuV89KVf1Gdtlmarne9/Wn1so0+IXVsym6ILjwVA74AlarLly9XSTevzaSeaKuqlfqdR6Ri2lLlX9Su3yst9zv2nG71HU3SLtWLM6Xx3FV1zhBdeG4msbk5pRiES/ft22edtQ9i4I9wmpbCB1K7ZqfyL+rnZ8jpJemSvWOnlP7bH5S/YYg+PBcDibYVK1ao49ChQ9ZZ+yCGSZMmWffcoflKvlyaFi/Zr3wod//2Fcn/5esqSmXoe5CIJfqoN9QPBc/FQN6ATHLnEKldCMEOGTJElixZogRGsu7s2bMqJ9HS0mI9K3xy8/NkDzNZ6lYp+TxGOd7VqZuk6Wr4C5QMzsH/fMGCBSpqGSqei4EyC/ICDNxbt25ZZ+2jfQbKORAYsw1rJRITE5WTTiKP9RPff/+9er/6+nrrN4OHgkJeb/bSRLk9N1Xaikulre6R1O/+Xiqm+/yLSQlSu26PtBR1ZMYNvQdi4MLY58XAsk9UixDIMnM/XLpzoCn/oEz8hx9+kK1bt6pGBJR48GVRB0UG/OLFiz126CgvL5dXX31V/uU3/yo5IyaLXH/6C2spKpHarN1KFBXTU6R+zzFpqw1ddIbwiRgxHDx4UEWSCK2G2yJGYyeaxOBGBHv37lXlHHyeRYsWSUpKikoIIh5EhJg01ESNme57n6x90tRNnqHpyi2pTtkgZWMXSuWClfL45AWR1tDtV4M9IkYMlGJQk/Ttt9+qxgBXrlyxHrGPHTEEAvMJM+ro0aPKrMK8wsxKSkpS0S9+ZmzeKJVL1krLxTzrt7qmvbVVHp/I9QkiQ8rJXyzbaPwLD4gYMUCrb5AAH9qOx98Zp8QQCD4rjjgOObPH6L/+VRJnzpH4ufN/5qzrvysQbTV1Pv/imFTMsPwLn0mFaWVwnogSA2B+OCEEcFMM/jDod1vlI9WdnHVmDcws7axTgEhyMaCzXlIhdT5nW/kXM3z+hc8JRywGZ4goMdy4cUPGjx+vapQiqYkYnxUfoysQOF/kyZMnnzjrmFn8Y/BL+F38lPJqv6WuRWXSnLlLamISpXrhKmVWtbf8NMMQsXIT/DYnEp99iYgSA040g5fFOt0NrmDxSgyYQXzRocKA83fW+UcR0UIsW3ftlNPXLsvdggKpPnJaGpLWS2Pscmlbu0fac29ITcJa61XcgX88nyeaiCgxPHr0SEWUiNw4EVHySgxOgvmEGUXuA7MqZdkyiU9NkeTMDMlMSZW9MXPl5Cej5Oq//1EaNviu3BW11m86C99dWuZqab+cL49PXbTORjYRIwZMCa6K2No0BWCWCJdIFEMgmHnIvZw9nyujp8TIoqREmffeJzL3dx/K+tf+S66NmCa1O7Kl7v5D+SngGx58d8s3rJO2vTnKqY8GIkYMVKkShVm9erVa5E/WOFyiRQz+ZPi+n8qCQmmIWy210iY3C+/K5Z0HpGhRhtweMV32ffR/kj5yrKzLypLjl85LftE9qbORWVdiWJ8lbYdOS93m8Mvp+wIRZSZVVlbKjBkzVJcM6pTCJRrFkOozI+9duS61c9OtM0/Teu6a3J4cL7t/N1Qy/uMTSfpmvMT7ZpLE1StlVdZalVlnnXhPZugTMRz+Ueq2hD9L9wUiSgxEL9ikhFDkzp07rbP2iVoxXL0hdXGrrDNdwMK8k5ekafFaKZ+YILcXrZRzO/fJ/mNHJXPbFklcmqzMUtaOUDJPpIt/ti6DoeZqte95cuy81G07rM5FOhElBlrFlJSUqPg8fZTCJSrNJJ9Tez/vllRPTrLO9Ezrwwpp2HpI6mPTpMlnXkn2GZG6BmHYFzx8IMdOnOhw1lNSVMiXn0S0psbHScmoOVI6IrbjhVyCWcqprcu6I6IcaAr1iMMTanRqPUM0iaG4uFhiExbJ2bhlUvh3r3W7BLUrmvPuSvXyLR1lIPNXSmtOR2sdTVt7mxQVFcmyjHQZ/vshcvUXv5W8v/mtpE6aLqmrMtQadaqL+SwMLifAR6RNkNtEjBho94LzTKOwlStXqvvhEm1iWJ2ZKc//4z/Iml+/IsW/eEnKpwQ/OwSi8cwVqVqc2dH4IDFLNVrTDPvzn+Q3//yiHH7/T9IwI02q827LjfxbKiFKN8Pk5GQ1i/CTrunUbbGOva4u9Iw5a1jI4rtNxIiBLC4Vo/gNFOtRuBcu0SYGchBvv/227F+RKY1zfVf1SmdMC1pu0nqzYlaalE2Il8fr9srWpFR585Mh8vBQjjQdChzM0P2qTp06Jdu2bVPVu/gh1GrpMnic9dLSUrXKLBA8hqi4EDoRQeyOiBEDIuBLxIGmLIMvM1yi0WdI8dny93w+Q48OtE2Uf+FzmG+PmS/bh42U+kmJUrsxtLIMLmqXLl1S/1MGuV5MhS+CSUTPXAYWpjG+4aBBg+TNN9903W+IGDEQxUAMlG6TYDIOdGCWp6dL0dUbXYZWneJBdaVkxiX4xJAkxe99qcyphiM/qvacdiCPRGadQkXMKgoXmUWoOHj99dfl888/V2aWE1sRdEXEiMENolEMPeUZnILvLm3DOpGC+9JwPFc1VquKX9PhXyxZ60izNDLrRA8RBMJgFuGCqEWinXWn1qzzfsxSdl7LUzHwpRDrJqqAM0af1HAxYrAP313qsp+3vWlvaFT+BS3+y33+Rc3K7dJ8u8h61B7klNhWQIO5hGWAs06YXTvrzCbMKnrNeqjOOn4LeZWu/Jfu8FQMfPm0dWEa4+rAVSFcolYMl69JzdRk64w7KDH4/JPuaC2tVNnp8smJKrJVt/WQtJZVWo8GjxrgNJpu9l2xuxin2lnXa9btOOtsq/zyyy+r/T9CxXMzCceKbChC4KoQLtEohgTfxaLwZr405fx0JXWDYMTgT/OtQqlJ36bacFbMWS6PDp9Ss0gw3C4okJKyMqlLXq8EFgqdnXUupswiOrPu76wzpgYOHGirCLRXfAaK9SjHcGRPtygUA+aEncK7UAlVDP40nrsmVQlrpWys7//o+9l4Nrj17HQjbH0Qfuk+SwHYsgBnnTXqRLGYRfAXPvjgAzWLXL16NaQx5rkY8BVIuKFgMpz+oHD+QUD6nu57OELEuPE3KPLD3/CfIvlSutufwdA1dADBRg+H9sdN0pB9Wip9MwUZb2YOZpCu4Hnh7I/RHSwlxhFnULM+nXHm76xv375d5bq6yqwzjjwVA44UzhIDnOlPw2ow7EIyn4CNqadDZhIOzjFN6uWixLOHDh0qX3/9tbpvCA0GBd+nU7SWV0n99sPKt8DHqNt8UOU0/KHbeedzTtNZ4Djr5LUCZdaphuACTUiYji36YhwKtsVATQyOEUsg/d9YZzgZ8FztmeJwtvmwqHzy5Mkq6sBtmhcDAnnnnXdkxIgR6r4hNLiSYme7AdGnmlU7pGyCz7+ITZVHB3/wzSKNPjMpVdocyqoHItg8Q+fMOjPHRx999GS/wFCwLQauRrw5BXv+hVv8U4ggkJBDMEzhFPMhDJSLfVjgc8KovPRvZU/ozZhJfRvqodjkpezbBZL/6zdUK063CFYMgWCdjadmko4v07mOI1xQd7Q50NEK/sL1AS9KwYv/qdZfsPFks0MbT2qeiCGAT9ATnjvQTEN8WKYlU7Xav2Djl7KxcWqmUP7Fd0ekQm88uXGfI461EkOiTww+kyzUshLPxcCsQDgPv4HYcbgYMUQWVXEZaockf1ru3u/YeHKCtfHk/hPSVm+vZxRiSEz1vcblm1KTttk6Gxyei4HwFjFgEig4y+FixBBZ9BRabbyY92TjyUo2nqSFTQilFVoM9eevSfWi0Kp+PRcDDjKzAokSJyIZRgyRRcWUJGm933PSTW08eexc0BtPap6I4YLPaU9YY50NDs/FQBkGvUnJFVBjEi5GDJGF2rsixN6yHRtPHu3IX0zq8C9aushiR5QYyCwTQiXBZiem2xkjhv5Fs8+/qM3c6edf5DzVk7aFEu5IEANZ5tzcXJU4o28St8PFiKH/ovwL8hf4F3GrpPHUJWlpbpbFazLk0fX8vi0GZgTWM5A0Iz1OFjpcjBgMyr84fk6q5q+U0rELZfv/jpbqrQelJnWT9Yzg8FQMZJCp/1izZo0qrzAzg8FpmitrZMPXE6Vk1Gy5/96XUrfpwJON7XvCc5+Btc+EVimvcAIjBoM/Le1tsnhVujQUlUjVrLSOhUkTrY1hlPPedWm852JwGiMGgz9Ek5JWZ0j1rqNSMnScddbaeHLphg7/YsFKacg5r/bd88eIwRBVtLS1qWhSwVv/I9cH/JPq9uGP2njSJwQEgTCqUzYqocD0uXPk7p3Qa6WMGAx9klafGBLmLZAHY+Kk+L3hyrHuCkym+j3H1cb29dNSZOJ/fyqFD0PfeLLPioGObayIY7lfJB18ZicKF/s7Kum2PE0aCu5J7Yqt1tmeaSt6KDGfD5fCTqsvg6FPioEBxVoIOn13PlgXwaAL9Jj/oQdnV/d7Omh/w+cI9rw+SEayhqO7LXQNPRNOBnr6rNjoMZNYJ633HugMPXiCofNgpH6KLzhYulpUQm6lp14+rN81s0N4RFQ5htP4i4HyDmaAzrBn86hRo9T66+4GNkKYOXPmU9W0lJkjMhYl6fIRrvIMXH8IGXN1J5FIbx+WVCLMO9aVho4X1GIxC3RVhmLEED78f+OWJBgxdCUGFoE/++yzMmzYsB7NEFodUlWryc7OVgvJ6dFDpwU6vSEYmlTppCHraOfNm6dars+ZM0dmzZqllqcikClTpqgVfSxx3bVrl1qiyu8GmqmMGMKH2feG77t9fPmmVMVnWmeDo1+IAdPl3XffVT12eoJqWnrwaFh/ERsbqwYyV3e+MM4xY9CsAOjVo4VBRzc+B50X6N6BsGJiYtRtFjWxjoM13YH2wDZicI6mq/mqkXIo9AsxAEtNMXd6AiGwd4Q2ZejlxEDH/GFmwASinxOtRXgMuM3MwKzBlr60wGFXIhof0LOH8hNaJnKO16JgkdfpjBGDcxgxdCMGrupsBNIT2skN5Ox25wDjIwQayKF0aTBicI6mSzelcs4K615w9BsxRAJGDM7BgqCmy8FFEDX9RgwtxQ+lvSn0Xjr+MDN0NzsES1evY8TQu/QbMVQuXCUtPay9JdKEb6F3ISWCRCRKm1eESQNV2wYa2JhHdIL2/3L1QGctB8nBzhgx9C79RgzFv/+6x51oGNSUl+to0vDhw5UYtN2P30EYlUVKhEdpe4NguK0dYj4DfTwRAj0/x44dq55PboIIEw42GzsSaiUKFRcXp9Z5gBFD79IvxEAz3KsD/l4e/HEcI946GxgGeGZmRxSCBBkDWHf+I99AHoHoEAOYfAPCIcyqW7sTbkVQgBiIHCEY8hV0hWa5K0k9+sSyDpyIlE7yGTH0Lv1DDDuPSuFrf1CdoXvyGxjcX331lRqY5AIY/KdPd2wLy8Ali80g5urPoKfbBzVFtL4BGiazmyk5DXILfD7WfNNBkFmAWYTXRWDMHrTOYY9kMGLoXfqNmfRoX05QDjRZY+x5Sim4+nNfo0s5qDNCAMDgJazqH0Llc/AaPF8/xk++aH6X85zjJ36K/l0jht4lqsTQXaFeJGDE0LtElRi6K+Hu64cp4e59+pQYqAz1HwwMEiASwxUTc6LzTo8IgL76GswQFviwb0MkHXxmTClD7+H5NlZdgYlA9EU7k4QgidDgiNKOEkdVb2elYVaggG78+PFPFuEgKLcOBit2fqDHnDj0AiB+InjOcV//bU4fvDYXDy5A3HbrvfRr8/3h07n1Phy8NhcXvj99v/NzAh0875tvvnlSch8KjouBIjcOSp/1fXb3IfqCGNLT09W6BKIyN2/eVM9hmytKs1944QV1282DPk9smfX++++rsuxAz3HioMJ1zJgx8tJLL8nIkSPV/UDPc+Lgb/r444/ljTfeULcDPcepg79j8ODBqozezffitckPsRc0F8lgvz9C3uwNyAwdKo6LgagNg54SZ2xnMr1UiTIzsAkfYcq8vDw1W1AlCjjLJMWee+65p0wnt3jrrbfk008/te65B1WvAwYMUDkKt2HwDBo0yLrnLs8//7zK27gNe/4988wzngVTXPEZKHtgUOMbMJ0y3QFqxTzhZ2flMnCI/XshBkq3yR67Df4RayYwk9yGJCPJRC8gv6Kz7W7CmMCs9soHc0UMBkMk0mfEQCkETjUOkNtgplFn5DaYiPhGZKfdBBOT92J2oFTEP8HoJARHeH1md2ahYFYc2gHrgLIX/kc4wqwy5H/mNn1GDPxDKWsItIzSSTDhWNKJU++2ScZgYbWd238TAwfTD6iPckt8mCv4fQhv9OjRyqx1y57HxCYIw8pEltfqpblu0mfEgG1NSJbF+G5CzoO6JSJJbtuilJDTVAC71024kLA0FYdz9uzZKszqBvh/xPBZ8kpUkECJE7s2BYLXpk6Mxg4USOpQvZv0GTFwxXZi19BgsRN6CxV6NTFwdADBLZgJMCNwagkMBLMs1g68LsWOxPMpmQm0JsQpiMBdvnxZzRCMCy8iSsaBNhgsjBj6MMws/mUtOK4abmszj6snzwulONCrcGUkYcRgE1rKMJVjNuhMuj96O2AGKo/jN3AbcNz97Xp/s0aXrvPa2P/cZ+DiSOL44yzr3k+YRfg/+Ft8HooE9e+DNi14LxxQfR/fgucTGdJlD7wH+ZADBw6o23xG/Xx+um3q9QWMGGxC/QvVkQwyDpaREtrEaabHEtEqFg7h/LE0dMSIEWrgYt+zqIjyFHow4cgTLQGcRQSAM0zWHucbiBSx0IirP4OS7D3vg8B4Xd6DkgXeh0FO9AqHE/HwHD4PJSEMcux9lrDyOnwWhEgIk8/L5/jss8/UeRxkfo8IH5+DzxrtGDHYhAFJFGfo0KFKDFytEQRlCgx2ZgwGGYWKREYY4Di5hCSpWSJSwrJTBr2G1Xg4psTVeX0dLuU1iLRpGKiEODdv3qzem1ovwpy8N4Od5zKAaVjA63P1p5sgphXRH+p8gOWufA5W8PGZeF1EhBi+/PJLFdpEWMwg/K3+Zlo0YsRgE9ZXY1ZwJSUczOAl7IjJwtWUWYKrKQeDjccQBSLhSs4gpvkAItHoGYBcATkKxAWYPiw1ZWbh6s0yVQapNpkww7jNa/P7iJDHSZJxnhmCga/9BAY2nwXxICI6iSAeZhaExPvq0CYRIz4Pj3sR0elNjBhs0jlhh12ua5D0Y/zUtxnQlHQDV9ji4uKnnGMN9Tj+v+8PC560L6Jfu/Nz/eP++hy+hBaChiy1Di/zubRPwOvjiPMZERN+Dj8DfdZow4jBYLAwYjAYLIwYDAYLIwaDQSHy/4I4MjoR3HKTAAAAAElFTkSuQmCC"
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the graph above the offset corresponds to the state preparation and measurement (SPAM) error of 0.0046 reported by Egan.\n",
    "\n",
    "The slope gives the error per Clifford gate.  Since the error is small the graph is approximated linear.\n",
    "\n",
    "For comparise actual Egan's results are shown below.\n",
    "\n",
    "![Egan%20Figure%20S1.png](attachment:Egan%20Figure%20S1.png)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}