RX_Tomo_Armonk.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook we use Quantum State Tomography (QST) to obtain fidelities of Rx(θ) interactions, using:\n",
    "\n",
    "1. a standard decomposition using 2 X90 pulses\n",
    "2. an optimized implementation using a scaled X180 pulse\n",
    "\n",
    "The state on which we do tomography is produced by the following circuit and produces a sinusoid probability of reading out 1.\n",
    "```\n",
    "        ┌────────┐\n",
    "q_0: |0>┤  Rx(θ) ├\n",
    "        └────────┘\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from copy import deepcopy\n",
    "\n",
    "from qiskit import IBMQ, pulse, schedule, transpile, assemble, QuantumCircuit, QuantumRegister, ClassicalRegister\n",
    "from qiskit.pulse import pulse_lib as pulse_library, reschedule\n",
    "from qiskit.ignis.measurement.discriminator.iq_discriminators import LinearIQDiscriminator\n",
    "from qiskit.ignis.measurement.discriminator.filters import DiscriminationFilter\n",
    "import qiskit.ignis.mitigation.measurement as mit\n",
    "from qiskit.ignis.verification.tomography import process_tomography_circuits, ProcessTomographyFitter\n",
    "from qiskit.ignis.verification.tomography import state_tomography_circuits, StateTomographyFitter\n",
    "from qiskit.ignis.verification.tomography.basis import PauliBasis\n",
    "from qiskit.visualization.state_visualization import plot_state_city\n",
    "from qiskit.quantum_info import state_fidelity\n",
    "from qiskit.circuit import Gate\n",
    "from qiskit.scheduler import schedule_circuit, ScheduleConfig\n",
    "from qiskit.quantum_info import state_fidelity, process_fidelity\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import qiskit as q\n",
    "import qiskit.tools.jupyter\n",
    "from qiskit.pulse import pulse_lib as _  # only want to call pulse_lib as q.pulse.pulse_lib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pulse_compiler_helper_fns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<AccountProvider for IBMQ(hub='ibm-q', group='open', project='main')>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q.IBMQ.load_account() # Load account from disk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "provider = q.IBMQ.get_provider(group='open')\n",
    "backend = provider.get_backend('ibmq_armonk')\n",
    "defaults = backend.defaults()\n",
    "circ_inst_map = defaults.circuit_instruction_map\n",
    "measure = circ_inst_map.get('measure', qubits=backend.configuration().meas_map[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### experiment parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qubit_t1(qubit):\n",
    "    return backend.properties().qubits[qubit][0].value\n",
    "\n",
    "def qubit_t2(qubit):\n",
    "    return backend.properties().qubits[qubit][1].value\n",
    "\n",
    "def sort_qubits_by_t1():\n",
    "    \"\"\"in ascending order\"\"\"\n",
    "    qubits = [i for i in range(backend.configuration().n_qubits)]\n",
    "    qubits = sorted(qubits, key=qubit_t1)\n",
    "    return qubits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sort_qubits_by_t1()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.         0.14279967 0.28559933 0.428399   0.57119866 0.71399833\n",
      " 0.856798   0.99959766 1.14239733 1.28519699 1.42799666 1.57079633\n",
      " 1.71359599 1.85639566 1.99919533 2.14199499 2.28479466 2.42759432\n",
      " 2.57039399 2.71319366 2.85599332 2.99879299 3.14159265]\n"
     ]
    }
   ],
   "source": [
    "shots = 2000\n",
    "thetas = np.linspace(0, np.pi, 23, endpoint=True)\n",
    "print(thetas)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Direct RX(θ) Gate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "rx_schedules_std = []\n",
    "rx_schedules_opt = []\n",
    "\n",
    "for theta in thetas:\n",
    "    circuit = q.QuantumCircuit(1)\n",
    "    circuit.rx(theta, 0)\n",
    "    \n",
    "    transpiled_circuit = q.transpile(circuit, backend, optimization_level=0)\n",
    "    schedule = q.schedule(transpiled_circuit, backend=backend, inst_map=circ_inst_map)\n",
    "    rx_schedules_std.append(schedule)\n",
    "\n",
    "    with qiskit.PulseBackedOptimizationContext():\n",
    "        decomposed_circuit = circuit.decompose()\n",
    "        pulse_compiler_helper_fns.update_basis_gates_and_circ_inst_map(decomposed_circuit, backend, circ_inst_map)\n",
    "        transpiled_circuit = q.transpile(decomposed_circuit, backend, optimization_level=0)\n",
    "        schedule = q.schedule(transpiled_circuit, backend=backend, inst_map=circ_inst_map)\n",
    "        rx_schedules_opt.append(schedule)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set up State Tomography"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qst_circuits_and_schedules(schedule, theta):\n",
    "    \"\"\"\n",
    "    create the state tomography circuits\n",
    "    the state preparation circuit is (H $otimes$ H).ZZ.(H $otimes$ H)\n",
    "    \"\"\"\n",
    "    num_qubits = 1\n",
    "    qubit = [0]\n",
    "    qr = QuantumRegister(num_qubits)\n",
    "    my_gate = Gate(name='rx_{}'.format(theta), num_qubits=num_qubits, params=[])\n",
    "\n",
    "    # create the tomography circuits\n",
    "    qc = QuantumCircuit(qr)\n",
    "    qc.append(my_gate, qr[:])\n",
    "    qst_circuits = state_tomography_circuits(qc, measured_qubits=qr)\n",
    "\n",
    "    # transpile so that the singles are in terms of the u gates\n",
    "    qst_circuits = transpile(qst_circuits, backend, initial_layout=qubit,\n",
    "                             basis_gates=['u1','u2','u3',my_gate.name])\n",
    "\n",
    "    # add the custom gate to the circ_inst_map\n",
    "    circ_inst_map.add(my_gate.name, qubits=qubit, schedule=schedule)\n",
    "\n",
    "    # schedule!\n",
    "    schedule_config = ScheduleConfig(circ_inst_map, meas_map=backend.configuration().meas_map)\n",
    "    qst_schedules = [schedule_circuit(circ, schedule_config) for circ in qst_circuits]\n",
    "    \n",
    "    return qst_circuits, qst_schedules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "tomo_circuits_std = []\n",
    "tomo_schedules_std = []\n",
    "for rx_schedule, theta in zip(rx_schedules_std, thetas):\n",
    "    tomo_circs, tomo_scheds = qst_circuits_and_schedules(rx_schedule, theta)\n",
    "    tomo_circuits_std.append(tomo_circs)\n",
    "    tomo_schedules_std.append(tomo_scheds)\n",
    "    \n",
    "tomo_circuits_opt = []\n",
    "tomo_schedules_opt = []\n",
    "for rx_schedule, theta in zip(rx_schedules_opt, thetas):\n",
    "    tomo_circs, tomo_scheds = qst_circuits_and_schedules(rx_schedule, theta)\n",
    "    tomo_circuits_opt.append(tomo_circs)\n",
    "    tomo_schedules_opt.append(tomo_scheds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUgAAAEvCAYAAAAn5VjRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAed0lEQVR4nO3deXRc533e8e/vDlaSIAECIrgA3EGJOwWB+ypKlV3bkeseu7ay2HFjs0lPe+y0ieMlsZulTay2dpKe2rJiJ6d23DS2qtiuY0d1XNnWsVZKFGlZEkVwExeRIrhTFAjivm//eO8AQ0qvBJDAzAB4PufMmZk7g7k/3HfuM+/dzXuPiIi8VlLqAkREypUCUkQkQgEpIhKhgBQRiVBAiohEKCBFRCIqijkyM9sGbAMYP378LQsWLCjm6EVkDNixY0eX9/6GofgsK9V+kO3t7f7HP3moJOMWkdFrYt2EJ733HUPxWVrEFhGJUECKiEQoIEVEIhSQIiIRCkgRkQgFpIhIhAJSRCRCASkiEqGAFBGJUECKiEQoIEVEIhSQIiIRCkgRkQgFpIhIhAJSRCRCASkiEqGAFBGJUECKiEQoIEVEIhSQIiIRCkgRkQgFpIhIhAJSRCRCASkiEqGAFBGJUECKiEQoIEVEIhSQIiIRRQ1IM9tmZtvNbHtXV1cxRy0iMmhFDUjv/b3e+w7vfUdTU1MxRy0iMmhaxBYRiVBAiohEKCBFRCIUkCIiEQpIEZEIBaSISIQCUkQkQgEpIhKhgBQRiVBAiohEKCBFRCIUkCIiEQpIEZEIBaSISIQCUkQkQgEpIhKhgBQRiVBAiohEKCBFRCIUkCIiEQpIEZEIBaSISIQCUkQkQgEpIhKhgBQRiVBAiohEKCBFRCIUkCIiEQpIEZGIogakmW0zs+1mtr2rq6uYoxYRGbSiBqT3/l7vfYf3vqOpqamYoxYRGTQtYouIRCggRUQiFJAiIhEKSBGRCAWkiEiEAlJEJEIBKSISoYAUEYlQQIqIRCggRUQiFJAiIhEKSBGRCAWkiEiEAlJEJEIBKSISoYAUEYlQQIqIRCggRUQiFJAiIhEKSBGRCAWkiEiEAlJEJEIBKSISoYAUEYlQQIqIRCggRUQiFJAiIhEKSBGRiKIGpJltM7PtZra9q6urmKMWERm0ogak9/5e732H976jqampmKMWERk0LWKLiEQoIEVEIhSQIiIRCkgRkQgFpIhIhAJSRCRCASkiEqGAFBGJUECKiEQoIEVEIhSQIiIRCkgRkQgFpIhIhAJSRCRCASkiEqGAFBGJUECKiEQoIEVEIhSQIiIRCkgRkQgFpIhIhAJSRCRCASkiEqGAFBGJUECKiEQoIEVEIhSQIiIRCkgRkYiiBqSZbTOz7Wa2vaurq5ijFhEZtKIGpPf+Xu99h/e+o6mpqZijFhEZNC1ii4hEKCBFRCIUkCIiEQpIEZEIBaSISIQCUkQkQgEpIhKhgBQRiVBAiohEKCBFRCIUkCIiEQpIEZEIBaSISIQCUkQkQgEpIhKhgBQRiVBAiohEKCBFRCIUkCIiEQpIEZGI6wpIMztgZj8aolpERMrKsPQgzWy1mf2jmZ03s3Nm9g9mtmI4xiUiMlyGPCDNbA3wY2AO8GngM0Ab8JCZLR3q8YmIDJfh6EH+OdADbPLef957/3lgE+CB/1r4xt/86EeYWDfhitt73v3uYShJypHaf+waKW1fMZA3mVkrIdzeAhihh/jR13nffGAl8Jfe+yP54d77I2Z2P/B+M5vqvT8GsHv3bgByuRy//uu/Qd3EOubPb7vOf0nKUXd3NzU1NVcMU/uPDSO57d80IM2sHvgJ0ArcAzwLbAYeBGqvevvK7P6Rqz5jLvALhHC9Bfj7/GuNjY2cPXuWzs5Ovnnffdf4b0g5279/P+94+9v46tf+mltuueWK19T+o9tIb/uBLGJ/DJgNfNh7/2+891/w3r8XuB9ouuq907P7I1cN/2OgMns8o/CFhQsX8sd/8ic88MA/8MADDwyqeBkZfv8/fIYLF16hpaXlNa+p/Ue3kd72AwnIfwYcB7561fDPvs57x2X3l/IDzKwSuBP4/lXv6fOhD32YxsZG/s93vjOAcqQseA/OQZpCb2+4Xb4cbj094XbpEpcvXOB73/sev3zXXTTX18OlS+FWQO0/Ag2g/aNt39vb9zHl3vYDWQc5F3jCe58WDvTev2RmZ65678XsvrpgWBNQAxwrfI+dOMG2mTOhtpaae+5h/oQJHH30ESq//BeQ5PC5XFggl4FxHktTSHvDlzbNvryD4X3BzWX3hBmBbDjhYXjs+/+u7zOufNB14QLd3d3ctLeTyrvvxnp6oPcyjBvHtqnToKKCmi99KbT/Y49S+ZWv4JMEkkTtP1Ce0PYuzQIre+z8m/9t4Yfk2975/vZ1BcPzjesKHkO0/a9ue7zHLlyAqkr+VUMDALX/6T/SVl3N0Z/8mKo//INrngSFZsDUIfkgAO/9G96AbuChyGungR8VPL+LMIk+VDBsHJAC38pe2wNszx7rpptuug3prRWOvlmuDfQ2kB7kPqDNzHKFvUgzmwbUX/XeJ7L7tcCXAbz3F83sp8Ct2T+w0Xt/rH35cv+TL95D8sJufvTgg7z9q1/lnne/h/dv2Rx6P85hLs268S77dXThl/Ly5axbn92n2fucC79uzmG+4LnPf57DY1nPxEI5nnBvCZiBGd4IzwGSfDcmvPa6+n5Bs8/L98AAu+IXmP5ekffZZ5L1lhLIhXuf9D/uu6+shMoKqKqGygp8rgIqcpDkIJf1uHO58P4k3PskPywXxpHLxf+H4dLdzR2f+TTHz5/n8c/eTdXyZbiWVjh7lmTvXpLdu3nwkUd4+ze/wT13/SLvX7cO0hRzLkzDwt6wc9jly1m7F3wPfH+753tO1jcs5Ype0dX/fr7pEgMMn93nvwv93xXepP2zDyrofVm+p+Wuavu+v7Ew3nyPOUnwZn3t1zc8l4OKSqiqhMoqfGVF1q4Vfe3d3/6FbZ/0t33+eTG75d5zx+9+iuOnT/PEJz9J1U0LcfPn46dN65uWDz74IO+88xf4whe+yC//yq8MyWgn1k04OiQfxMAWsb8NfBx4P/BXBcN/5+o3eu87zWw78B4z+z3vfb7QPwPuA7oIvckwUy9bxvbLl/ngb/0WN86YwV3tN2NnzuJaW6Gmpu+7+6a8z2ak7OayGawgWHEpluZnOlcQqP3hGRZPXf+M5um/L1zkvKKygpkpsRBmSfYlTyz74lbg+wIv/4UN7/H5GSDX/4UOwyquGDbiljd7e7GjR7FL3fzu+97H2+++m1/85jf4wtatNE2YABMm4GbM4MnKSj74R3/IjVOnctf8eeAcvrk5TIMB8Ve2c5qGH9HCds6Ha354YdunBT/EzmWLpwXfk77VDVwZhEDfD5xl7W5kP0JJ1vYVWcAVBFZh6FkSvhe5/h86kiR8Z/I/fhUVxf9RGwrnzpG8dJTfe8tbeNsXv8hd3/17vvCud9F0w5S+t+zYsYMPf+jXWLBgAf/ive8tYbFx5v0bx5CZNQBPE7Y+3wP8HNhC6CXWAj/33m8peP86wi5Ah4H/lg3+t9nfVwBngPubm5s/vGLFzfzgB/+XmTNn8nf338/8JEfyzM9IDhzE19Tgp0/PwkFGDO+xkyexrhP4qdNwS5bgFizgL772VT7227/NpEmTuPPOO2lsbOKZZ57pb/+//jptly6RPPcsdvIkbtp0mDSp1P+NDFZPD3bkMIbh2ubjFi3m3u9/743b/lvfZt68eUNWwsS6CU967zuG4rPeNCABzGwm8DngjmxQfkfxHwIHCgMye/9a4I+A1YSf3IeBTxC2mv8OsMXMmubOnce7/vm7+MhHPsqk/MzQ3U2yZw/Jrp3Y8eO4qdOg/uoleSlL3d0kL76InzAev3AR6eLFV7TdU089xZ9+/nM89NBDnD9/ntbWma9pf3v5OMnOXSSde/Cpw7e0QFVVqf4jGSjvsa4u7GQXrnUmftky3Pz5fR2cgbT9UCl6QA6H9vZ2/+OfPBR/w6lT5HbsIHlhN955fGtrWNyQ8uM9dvw4du4cbs5sXPstIdiuddHQOZI9e7CdO0kOH8JPacZPnjykJcsQunSJ5OBB/IQJuCWLcUuWwoQJJStnKAOyfBNn8mTSW2/Ft7SQ7HiKpHMPbvoMmDix1JVJofzMMbGOdN063LJlcNVhZYOWJLgbb4Tp0+GpJ0mefx7OncPPnBnW30nZsK4urKsLN3s2rqMjdGRG4jrTiPINSOibUVxzM7knHid5fjf+/PmwFUwzSumdPk1y/Bhu9hzcypVh5hhKdXWkmzbjm6eSPPkkyQu7cTNnwbjXHGsgxZam2MGDWFUlbs0a0hUrRmW7lHdA5tXXk269DT+lmWTHDpLOTtzs2Vo3VSrOYYcPY97hbrmFtGMljB8/POMyw910E76pieTxx0leeAHf0IC/4YbhGZ+8uYsXSV48iJ82jd6Olfj580dVr7HQyAhIgFwOt3x5mFEefZRk396wP10J13WMSZcvk+zbi29sIr2lHbdwUVF6876pifS22/CNk8nt3An79+NnzdKSRJGFPRS6cDctJF29GhobS13SsBo5AZnxM2aQ3nEHPPooybPP4uvr1ZsolgsXSA4fws2ajVu3LqzqKKbqatzKVVDfQPLEEyR79uDmzNGSRDF4jx06FJYaVq8mbW+//nXNI8CIC0ggrJvasgU/cSK5p3fAgQOhNzFKu/nlwLq6sFOncEuWkK5ZW7qNZdl6aT9pEskjj2hJohh6e0kOHMDX1ZGuWlm0pYZyMHL/y8pK3MqVpJs2wbhxJHv3XnGWEBki3mNHjmDnz5GuXk26eUtZ7Engp04l3boVt3QZyeHD2KlTpS5pdMr2S/bNzWF6L14yZsIRRmoPMi9JcAsX4ceNJ/fIIyR7O3Hz5mt/yaHiHHbwIFRXk27YiFu6tLxmjkmTwg9kTQ3Jzp1w+TK+ubnUVY0er7xCcuhFXFsbbu0afMFhgmNFGX3br52fNYv01ltxrTNJOjuhu7vUJY18vb2hVz5xIm7jRtzy5eUVjnnV1aTr1pGuXo1d7sEOHSp1RaPDuXMkhw7hFi0Oq7PGYDjCKAlIICwC3HYbrq2N5MB+uHCh1CWNXL29JJ2d+OYppJs2hZ22y1kuh+voIF23HiorsH37snNYyrWwkydJjh/HrVhBumVLWaxSKZXRtSza0BAatLqK5JlndOTNtbh0iWT/flxrK27zJvyUEbLIaoZbsgRfU0Pu4YdJ9u7FzZtXnr3eMmbHj2OvvBLW769aNeZXV42+/76ujnTjJqiqInl6J845nexioF59leTgAdy8eWEajsDjn/38+aSVleQe/mk4oGDu3DE/kw+UHT2K9fSQrl0bVqlouo3CgASorQ2LW2YkT+/EO6eTHbyZ7OgId9NC0nXrRvSPip81i7SigtxPs5Ccrw13byi/pwL0H0+vnjcwitZBvkZVFen6DbiVK7HTp7GXXy51ReUrv7Vy0WLSjRtHdDjm+RkzSLdsxrW2hA13PT2lLqls2cGDkFh/z1Hh2Gd0T4mKCtJVq0hXrcLOn1dIvp4LF0I4LlwUdpmpqyt1RUPGT2nG3boVN3sWyb69r7ma4pjnPXbgANRU4zZsCLtx6WCLK4zugASoqAgrnNeswV69iB0/XuqKyseFCyRHDuOWLgs9x1F4Nhbf1ES6eQtu3jyS/fu0C1iec9j+/TBuHG79BtxNC0tdUVka/QEJYYfy7Kwz1t2NHR2ya/qMXOfOkRw+jFu8mHTDhuE7G085aGwk3bgJNz/bBWys9ySdww4egIkTSdevL//duEpobAQkhJDs6CBduxbrvYwdOVLqikrn7FmSo0fDfm6bt0BtbakrGn6TJ5Nu3hxCciz3JJ0j2bcvhOO6dfi2tlJXVNbGTkBC2Fdu2TLSteswH85pOOacPUty7FjoUa9bB9XVpa6oeCZNCkdc3Xhj6ElevFjqioorC0ffUE+6cVM4j6O8obEVkBB6ksuWka5Zi5mNrZA8c4bk5eO49ptJ16wZE6ereo1sP1m3aDHJiwfHTkimKcnevfjGxnCW9tmzS13RiDD2AhL6e5Jr1kBiYTeHUc5OnSI5cQK34uawj+hY6jlebeJE0g1hw8SYCEnnwkmOp0wh3bgxnBpQBmRsBmTGLV2KW78eKivD7g4lusLjcAuX4zzZv1hdWVnqkkpvwgTSLVtCT/LQi6P32P3e3nC6snw4zpxZ6opGlDEdkABu4SLchvVQWxt2exhlJzmwEyews2dIV68Kp8hXOPYbNy5suFm6jOTIYTh3rtQVDa3e3rDOcdo00lu3hkvxyqDo+CvA3bQQn+TIPfpoOFHDnDmj4mgCO34cu3iRdOXKcKmC7CLuUqC2lnT9eqjIkex4evQcu58/I9P06aSbNuGnTi11RSNSUQPSzLYB2wBah/oSodfJL1gQjt99+OFw4t2580Z0oNjhw1iakq5ejWtvHxWBP2zyx+5jJLt2jfxj93t6SPbtC2dk2rJ5zJ7LcSgUda7x3t/rve/w3nc0NTUVc9QD4ufODScHbW4Ox+9evlzqkgbP+3C9YggbIhSOA1NVFaZXR0c4dv/YsVJXdG0uXgzX6Zk3l3TrVoXjddKccxXf0kJ661Zca+vIOzt5moadgGtrSTduwC0ZW9cPuW4VFeG6O2vXYpcuhV3ARtKGuwsXwhmZFi0OBwCM8kuyFoPmntfhm5txW7b071A8ErZwXr4cVsg3NpJu3BiuPKcTDwxeLodrbw+HX+ZyI2fD3enT4eiopUtHzRmZyoE20kSEkxxshooKkuefxzc1le96qe5ukgP7ca0zcevX46dPL3VFI1uS4JYuDWcnf/yxcE7J2bPL9vrbdvQo1t0dDqVdtaps6xyJFJBvpK6OdOtWGD+eZNdOePXVED7l1DM7c4bk2Eu4tgWh19PQUOqKRg3f1kY6blz/9bdnziqvMx45F65XXVMd1p8uXTqiNyyWIy1iv5mqKtL160k3bIRcLqyXLIfrb2dngU5OnsTd0kF6++0Kx2HgZ8wgvf32cNTNoUPYyZOlLinIX6+6oZ50y63hRLcKxyGnHuRAJAlu+XJ8fX3YV7JzD66lFSZMKE09PT0kBw/g6yaGHZ0XLdIlBYZTfT3pbbdBfX1Ykjh3LhyRUqJAspMnsRMv49oW4FavGjkXVhuBNFcNgp81i96JE8ltf4LkuefxNTVhkbuIW4rt1Cns5eO42bNxK1fhy2x/0lGrupp0zRp8YyPJ9u0ke14o/o9kby926BCWJLjVa0hXrBjd5/EsAwrIwWpoIN1yK755KsnTT5O8sBs3o2X4L1XQ0xOOGa6uxq1cFWaOUXR5hBEhSXA33ohraiL3xBNhEffUqfAjOcw9+PwPo58+g9729nCqMu3CNewUkNeishK3bBl+ajPJUztI9nbiT5wIM8pQn0IsTcMhg+fO4WbOxN18M36UHAo5YjU2kt5+O761hWTnLpLOPfimG8JeDkPdLq+8QnLkCIyrxXWsJF22TLvwFJEC8jr4Kc1hRpkzB/vZz0hefBFfVRWOe73eoOzt7TvRhL9hSljXuHBheW1FHcsqKnCLl+BaWsnt2oW98ALJC7uHLigvXCA5dgxyCW7evPCDPHNmee1BMQYoIK9XRUW4psesWfjOTuz558MvvkvxDZPx9fUDX/xyDs6fD4e6dXeHU1St3xA+X72G8jRpEumGDVhbG8lzz2EHDoSgnFAX2n4w6yh7ekLbnzkTrjQ4ezbuppvCEoPOwlQSCsihUlMTDu2bPx936BDJ/n3Y0ZdI9u8PYVk7Dmpr8VVV/Vs/ncN6euDSJeyVV0JA1tXhZkzHzZuPa22FSZNK+3/JmzPDT51KOnUqduJl7OCL2L59JKdOwuFDUFmFHzcOqqvxFRWhd+l9WH3S0wMXL2Ld3VBZia+vx7XfjJs9Bz9jhvZOKDFN/aFWUxN2MG5rg9OnsRMnSE52wUvHsIsXw6Vn03Domk+SMFNMmoRf0BaO1mm6AX/DDVrHOEL5G6aEE0QsX447fjzskvNyuOdST7iqpnMhVCsqoKoK39KCa24O7d/cXLrdx+Q1FJDDqaEB39BAyoLQY+juhu5uLNvR3Ody4dIHtbUKxNGmshLf0tJ/ktreXnj11bAY7VLA8JWVoe3H8uUvypwCsljMwsxQW8sIOj+MDJWKir7dstT+I4e6LSIiEQpIEZEIBaSISIQCUkQkQgEpIhKhgBQRiVBAiohEKCBFRCIUkCIiEQpIEZEIBaSISIQCUkQkQgEpIhJR1IA0s21mtt3Mtnd1dRVz1CIig1bUgPTe3+u97/DedzQ1NRVz1CIig6ZFbBGRCAWkiEiEAlJEJEIBKSISoYAUEYlQQIqIRCggRUQiFJAiIhEKSBGRCAWkiEiEAlJEJEIBKSISoYAUEYlQQIqIRCggRUQiFJAiIhEKSBGRCAWkiEiEAlJEJEIBKSISoYAUEYlQQIqIRCggRUQiFJAiIhEKSBGRCAWkiEiEAlJEJEIBKSISUdSANLNtZrbdzLZ3dXUVc9QiIoNW1ID03t/rve/w3nc0NTUVc9QiIoOmRWwRkQgFpIhIhAJSRCRCASkiEqGAFBGJUECKiEQoIEVEIhSQIiIRCkgRkQgFpIhIhAJSRCRCASkiEqGAFBGJUECKiEQoIEVEIhSQIiIRCkgRkQgFpIhIhAJSRCRCASkiEqGAFBGJUECKiEQoIEVEIhSQIiIRCkgRkQgFpIhIhAJSRCRCASkiEmHe+9KM2Ow8sLskIx+YJqCr1EW8gXKur5xrA9V3vcq9vhu993VD8UEVQ/Eh12i3976jhON/Q2a2XfVdm3KuDVTf9RoJ9Q3VZ2kRW0QkQgEpIhJRyoC8t4TjHgjVd+3KuTZQfddrzNRXso00IiLlTovYIiIRRQ9IM3urme02s04z+3ixx5/V0GpmD5rZs2b2czP7SDZ8spn9wMz2ZPcN2XAzsz/Pat5lZu1FqjNnZjvM7LvZ8zlm9lhWx9+aWVU2vDp73pm9PrsItdWb2X1m9ryZPWdma8tp+pnZb2Zt+4yZ/Y2Z1ZRy+pnZX5rZy2b2TMGwQU8vM/tA9v49ZvaBYaztP2dtu8vM/s7M6gte+0RW224ze0vB8GGZt1+vvoLX/r2ZeTNryp4P7bTz3hftBuSAvcBcoArYCSwqZg1ZHdOA9uxxHfACsAi4G/h4NvzjwGezx28Dvg8YsAZ4rEh1/jvgfwLfzZ5/A3hf9vge4Deyx/8auCd7/D7gb4tQ2/8APpQ9rgLqy2X6ATOA/UBtwXT71VJOP2AT0A48UzBsUNMLmAzsy+4bsscNw1TbHUBF9vizBbUtyubbamBONj/nhnPefr36suGtwAPAQaBpOKbdsM5Er/OPrgUeKHj+CeATxawhUte3gX9C2HF9WjZsGmFfTYAvAXcVvL/vfcNYUwvwQ2Ar8N2swbsKvrR90zL7kqzNHldk77NhrG1SFkB21fCymH6EgDyUzQwV2fR7S6mnHzD7qhAa1PQC7gK+VDD8ivcNZW1XvfYu4OvZ4yvm2fy0G+55+/XqA+4DlgMH6A/IIZ12xV7Ezn9x8w5nw0omW5y6GXgMaPbev5S9dAxozh6Xou4/BT4GuOx5I3DGe9/7OjX01Ze9fjZ7/3CZA5wA/ipbBfBlMxtPmUw/7/0R4L8ALwIvEabHk5TP9Msb7PQq1fzzLwm9srKpzczeCRzx3u+86qUhrW9Mb6QxswnA/wY+6r0/V/iaDz8zJdnEb2bvAF723j9ZivEPQAVhkeeL3vubgVcIi4h9Sjz9GoB3EoJ8OjAeeGspahmoUk6vN2JmnwJ6ga+XupY8MxsHfBL49HCPq9gBeYSw3iCvJRtWdGZWSQjHr3vv788GHzezadnr04CXs+HFrns9cKeZHQD+F2Ex+8+AejPLHx5aWENffdnrk4CTw1jfYeCw9/6x7Pl9hMAsl+l3O7Dfe3/Ce38ZuJ8wTctl+uUNdnoVdTqa2a8C7wB+KQvwcqltHuHHb2c2j7QAT5nZ1KGur9gB+QTQlm1NrCKsEP9OkWvAzAz4CvCc9/5zBS99B8hv3foAYd1kfvj7sy1ka4CzBYtGQ857/wnvfYv3fjZhGv0/7/0vAQ8C747Ul6/73dn7h6034r0/BhwysxuzQbcBz1Im04+waL3GzMZlbZ2vryymX4HBTq8HgDvMrCHrJd+RDRtyZvZWwiqeO733F6+q+X3Zlv85QBvwOEWct733P/PeT/Hez87mkcOEja7HGOppN1QrUQexsvVthK3Ge4FPFXv8WQ0bCIszu4Cns9vbCOudfgjsAf4RmJy934D/ntX8M6CjiLVuoX8r9lzCl7ET+CZQnQ2vyZ53Zq/PLUJdK4Dt2TT8FmHLYNlMP+D3geeBZ4CvEba6lmz6AX9DWB96OZuhf+1aphdhfWBndvvgMNbWSVhnl58/7il4/6ey2nYD/7Rg+LDM269X31WvH6B/I82QTjsdSSMiEjGmN9KIiLwRBaSISIQCUkQkQgEpIhKhgBQRiVBAiohEKCBFRCIUkCIiEf8ffEFKBIfqmoYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from qiskit.visualization.pulse.matplotlib import SchedStyle\n",
    "style = SchedStyle(d_ch_color=['r','r'], use_table=False, figsize=(5, 5), remove_spacing=True, dpi=800)\n",
    "tomo_schedules_std[22][2].filter(channels=[backend.configuration().drive(0)]).draw(style=style, plot_range=(0,1400), scaling=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUgAAAEvCAYAAAAn5VjRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAbx0lEQVR4nO3da2xj+X3e8e+PpEhdqNuMNJLmsp5Z726aBYrG6zFqN43hNqnjGkHcFjbqNECcpsUCLVC4bdDUrtEUBfrGaZukKZo4mzRGG7hpYjetXQOtUTtOvIbtXc/seHZukkb3+4UiKeouXv59wSOtVqv/jjRD8ZDi89kVSB2e4XnmSHzm3I855xARkbeKhB1ARKRWqSBFRDxUkCIiHipIEREPFaSIiIcKUkTEI1bNiZnZi8CLAG1tbe9+7rnnqjl5EWkAt27dSjnneivxXhbWcZAvvPCC+7NvvRzKtEXk7OpoT950zl2vxHtpFVtExEMFKSLioYIUEfFQQYqIeKggRUQ8VJAiIh4qSBERDxWkiIiHClJExEMFKSLioYIUEfFQQYqIeKggRUQ8VJAiIh4qSBERDxWkiIiHClJExEMFKSLioYIUEfFQQYqIeKggRUQ8VJAiIh4qSBERDxWkiIiHClJExEMFKSLioYIUEfFQQYqIeKggRUQ8VJAiIh4qSBERDxWkiIiHClJExEMFKSLiUdWCNLMXzeyGmd1IpVLVnLSIyIlVtSCdcy85564756739PRUc9IiIiemVWwREQ8VpIiIhwpSRMRDBSki4qGCFBHxUEGKiHioIEVEPFSQIiIeKkgREQ8VpIiIhwpSRMRDBSki4qGCFBHxUEGKiHioIEVEPFSQIiIeKkgREQ8VpIiIhwpSRMRDBSki4qGCFBHxUEGKiHioIEVEPFSQIiIeKkgREQ8VpIiIhwpSRMRDBSki4qGCFBHxUEGKiHioIEVEPFSQIiIesbADNIqSK7G8uUxmO81WYRvnHM2xBJ2JLi60XqAp2hR2RBE5RAV5yoqlIuOr4wyuDLKwPs/67nr5BQPnSrQ2tXG+pYfne36Yp7veSXOsOdzAIrKvqgVpZi8CLwJcuXKlmpMORXorzWuLrzG8Mky+lKe/rY+B5AARK2/ZcM6xkd9geWOJr6/Ncq1rlBf6382l9kshJxcRqHJBOudeAl4CeOGFF1w1p11tU7kpvjvzHWbXZrnUfomORMdbxjEzkvEkyXiS7cI2o5lR0ltp3nPxPfy58z+8X6QiEg59Ak/BSGaEP538Josbizx77tkjy/Gw5lgzz517ju3CNi9Pv8ythVs4d6b/DRGpeSrIChvJjPDt6W+ztrvGM93PEIscfyHdzHiq8ymaIk28OvcKry28ppIUCZEKsoLm1uZ4ZfZ7bOTXudp5FTN7rPcZSA6QiCV4bfEGgysPKpxSRI5LBVkhme0M3575NstbKa52Xn3i7YcDyQGKpRKvzL3KdG66QilF5CRUkBWwXdjm5vwNZnLTvLPr6YrtXHmq4ymy21m+P/cq2e1sRd5TRI5PBVkBd5fvci91jysdV060zfFRzIxrXVeZWJ3g5sJNiqVixd5bRB5NBfmEpnPT3Fm+Q3u8nWQ8WfH3j0fjXGq/xODKA22PFKkyFeQT2MpvcXPhJqvbq/S19Z3adDoSHcQiMW4t3mJlM3Vq0xGRN1NBPoH7qftMrk7wjs6nTv2g7kvJSyxtLnFn+Q6FUuFUpyUiZSrIx7SwvsDry6/THm+vyvnT0UiUi8mL3EvdYywzeurTExEV5GMpuRJ3U3fJbKdPddX6sM5EJwCvL99hu7BdtemKNCoV5GMYST9keGWIgbaBqp8vfTF5kencFK8v3a7qdEUakQryhHYKO9xL3aPkSnQ1d1V9+s2xZjoTnQyuDLGytVL16Ys0EhXkCQ2nh5nKTYd6SbK+tj5Sm8vcX75PyZVCyyFy1qkgT2B9d507y3doiTWHemHbiEXoa+tjKD3IwvpCaDlEzjoV5AmMZEZY3Figv60/7Ch0N3eT211jcGVQV/wROSUqyGNKb6W5u3yXZFOyJu4fY2b0tV5gNPuQydxk2HFEziQV5DGNZUZZ2lzkYvvFsKPsO9dyjo38Fg9XhrUtUuQUqCCPYW13jaH0EJ2Jzpq7DUJf6wUmchPMr8+HHUXkzKmtT3uNur98j6XNRXpbesOO8hadiU42dje4vfgDLUWKVJgK8hFyOzlGMiN0JDprYtvjYWZGX1sfs+szzK3NhR1H5ExRQT7CxOoES5tLXGi9EHYUr67mLtZ2NxjLjmqPtkgFqSDfxtruGneWXifZlKzohXBPQ1/rBYbSQ8ytaylSpFJUkG9jIjtBaitV1QtSPK5zLedY211jNDMSdhSRM0MF6bFT2GE0O0I8Eq/JbY9HOd98nonVSdJb6bCjiJwJKkiPsewoU6vTdbH0uOd8y3lSm8sMpYfCjiJyJqggj1AoFRjJjGIGLU0tYcc5tmgkSldzF6OZUXI7ubDjiNQ9FeQR5tfnmV+fq6ulxz29rb2sbKWYzk2FHUWk7qkgD3HOcW/5Lpv5TZJNlb9L4WmLRWIkognuBn8HEXl8KshDljaXmF2fo7e1FzMLO85j6W/rL/891mbCjiJS11SQh0ytTrK6nQ3lauGVkoglcM4xmhnT6YciT0AFeUB6K83gyiCdzbV3UYqT6m3tZWJ1nCltixR5bPXdAhU2szZDaitVExfEfVJdzV1sFbaYWp3U6Ycij0kFGSiUCgyvDNESa6n7pcc951vOM5YdJ7uTDTuKSF2qahOY2YtmdsPMbqRSqWpO+pHGs2Msbi7S09ITdpSK6W7uJrO9wmhapx+KPI6qFqRz7iXn3HXn3PWentopovIOjVFKrkRbvC3sOBUTi8ToSHQynHnIVn4r7DgidedsrEs+ocWNRWbWZs7U0uOe883nSW2lmNJ9a0ROTAUJTKyOs7qTozPRGXaUits75Gc8O06hVAg7jkhdafiCzG5nGU4Pc77lXN0eGP4oF5MDjGfHda1IkRNq+IKczk2T3kpzvuV82FFOTTKeZKu4zUR2POwoInWloQuyUCowsTpOPBqv+SuGP6lzzd1M5aZ0lR+RE2jogpxYnWAqN0Vva+3drbDSzrecZ3lzieH0cNhRROpGwxakc46p1SnypQLJeP1dteekYpEYrU1tTGTH2SnshB1HpC40bEFmtjNM5iboOcPbHg/rbellcXNRO2tEjqlhC/JhepjsdoauRP1eteekWppayBfzDK0M6So/IsfQkAW5ld9iNDtGR6KTaCQadpyq6k/2M7M2zfLmcthRRGpeQxbkdG6K1FaK882Ns3q9pz3eTm4nx2R2IuwoIjWv4QqyWCoylh3HuRKJWCLsOFUXsQjt8XbGV8fYyG+EHUekpjVcQc6tzzGxOs5AciDsKKHpa+tjfn1BB46LPELDFeTk6iSbhU3a4+1hRwlNU7SJpmgTo9lRnZ8t8jYaqiDXd9eZXJ2gu7k77Cih62npYXFjgaWNpbCjiNSshirI4fQwS5tLDblz5rD2eDtruxvcX7kfdhSRmtUwBblb3GUiO05rrJWmaFPYcUJnZvS29DCbmyG9lQ47jkhNapiCnF+fZ2FjgZ7Ws3dR3MfV3dxNejvNjO6fLXKkhijIkisxuDLITnGX1qbWsOPUjGgkSnOsmaGVQTbzm2HHEak5DVGQqc1lZnJT9Lf1hR2l5lxMXmR+fYHZtdmwo4jUnIYoyMnVSVZ3cnQkOsKOUnPi0TgOx8P0sO6fLXLImS/I9d11xrKjtMfbz8z9riutp+U88xsLLG3qkB+Rg858Y0ysTjC/vkCfVq+9upu7yW1nmc5Nhx1FpKac6YIslAqMZEb2zxyRo5kZHYlOBlcekN3Ohh1HpGac6YKcX59ncWOBC60Xwo5S8y60XWBla4VZHfIjsu/MFqRzjnvL99jKb9HW1BZ2nJoXi8SIR+I8zIxQLBXDjiNSE85sQS5vLjG3Pktva++Zvd91pfW09jC3NstYdizsKCI14cwW5PTaDNntLF3NjXNLhSeVjCfJl/KMZ8d0SwYRzmhBZrYzPEjdpzPRqUN7Tqi3tZfp3JRuySDCGS3ImdwMy5vLOrTnMXQmOlndyfEw/TDsKCKhO3MFWSgVGEoP0RJrabgbclVCxCL0tPYwkhkhs50JO45IqM5cQY5mRlhYn6e3tTfsKHWru7mbzPYKE6sTYUcRCdWZKshiqcjDzAiArtrzBGKRGK2xNkbSD3VjL2loZ6og59fnmV2b0dJjBfQn+5ldm9WNvaShVbUgzexFM7thZjdSqVRF37vkSjxI3Wcjv6mr9lRAPBqnKdrEUHqIfDEfdhyRUFS1IJ1zLznnrjvnrvf0VPbK3ksbS0zmJhlo66/o+zayC60XWFhfZDI3GXYUkVCcmVXssewYOV3zsaKS8SRFV+B+6p4OHJeGdCYKMrOdYSQzQndLtw4Mr7ALrReYW5vTpdCkIZ2JNhnNjLKytUxvi3bOVFpHooOtwjYPMw91EQtpOHVfkLmdXHBaYZcODD8lF5MDjGZGmFufCzuKSFXVfUGOZcdIbaV0aM8p6kh0sJHfZDg9HHYUkaqq64LM7eQYTg/R2tRKLBILO86Z1tvSw8TqhO5+KA2lrgtyLDvG7NosF5MXw45y5p1vOc/qzirDK7r7oTSOui3IrfwWD1L3aWtq09JjFZgZfa0XGMuO6u6H0jDqtiDvpe6xuLGo+81UUVdzF7ndHHeW7+i4SGkIdVmQa7trDKeHScaTJGKJsOM0jIhF6G/rZyI7zvz6fNhxRE5dXRbkaGaUxQ3d6zoMnYlOcrtrPEwPaylSzry6K8j0Vpo7S6/THm/XtscQmBn9bX0MrgwyuapztOVsq7uCHE4Psby1zEByIOwoDau7uZvt4g6DqQdaipQzra4KcnlzieH0Q7oTOuc6bANt/UzmJhnNjIYdReTU1E3LlFyJe8v3SG0tc6FNe67D1pHooFAqcHf5DjuFnbDjiJyKuinIubU5RjIj9LX2aemxRlxqv8T02gyj2ZGwo4icirpomnwxz63F11jbXaO7uTvsOBJojjXTFInxg8XbZLezYccRqbi6KMjR7AiTq5M81fEUZhZ2HDngSvsVFjfmGVx5EHYUkYqr+YJc3Vnl7tJdzEx3KqxB0UiU7uZzPFgZ1IUs5Myp+YK8n7rP1NoUT3U8FXYU8bjQeoHsdoY7S3d0gy85U2q6IOfW5hhcecC55nM6KLyGmRmX2y8zmh1hJPMw7DgiFVOzBblb3OXmwk1Wt1d1Mdw6kIwniViEWwu3tMNGzoyaLcjBlUHGV8e50nFFh/XUicvtl5nfWODu8l2dYSNnQk02z+LGIj9YvEVLtFk7ZupILBKjv62PO0uvM5LWqrbUv5oryHwxz62FW6S2VrjYriuF15tzLecoUeIHS7fJ7eTCjiPyRJ6oIM1swsz+tEJZALi9dJvh9CBX2i9r1bpOXW6/zExumlfnX9WtYqWunUoDmdlfNLOvm9mameXM7P+a2Y886s/Nrs1yd/kuyXg7yXjyNKJJFcSjcS62X2R4ZYih9GDYcUQeW8UL0szeC/wZcA34ZeBfAc8CL5vZn/f9ufXddb47+12y2xn62/orHUuqrDPRSdSi3Ji/qauPS906jSXI3wB2gfc7537NOfdrwPsBB/z7o/5AvpjntYWbTOcmudZ1TacTnhGXOy6T2kpxc+EmG/mNsOOInNixCtLMrpjZH5nZarDK/L/N7J1HjPcM8B7gi865/fPOgudfBH7CzN6yeHg3dZfbS7fpbxsgHo0//t9GakrEIlzteAcP00PcnL+h7ZFSdx5ZkGbWBXwL+FvA7wOfAjaBbwJth0Z/T/D43SPe6nuAAe8+OHA0M8qthddojbXS1dx1svRS81qaWuhvG+D15TvcXroddhyREznO+Xu/BFwFfsE59/lg2G+a2a8Dnzw07t5xOUddtWBv2KW9AUsbi7wy+z0285s83fX08VMLJVeiUCpQLBUpuiIlV6LoihRL5eclSjjncM5RogQOXPDfQRb8V/7fMDMiRMqPFsEoPx7+ilq0/BiJPvJog67mLrYKW7y2cJOORAfPdD9zmrNGpGKOU5B/A1gE/uuh4Z/lrQW5d1T3UZeY3j44TqFU4FvTL7O4ucSz3c9ouyNQLBXJl/LsFncplArslnYpFAvkS/n94TuFHXaLO+RLhf1SLLkSJVfCUaJUcuVHF9ShK2FmOOcA4/Bsdg7A7Y9jFqE8ihExw4KyLBdnUJoWIYKVH638elOkiVgkRjQS3X8ei8TeKFUizK3P8ZXhL/OXLv8o7+p7l27ZKzXvOAX5NPB959ybNiA55+bN7PBJt5vB41G/+c0Hx8lsZ/j+3Kt0Jbp4sHL/JJnrjuON8iuU8uRLBfLFPI4SzhEUXYFCqRgsERYoutIbS4fB48Glv6OW6vaXAA88BzBnbw7jCwk4VwjKlKB0HTjeWCLFlYcf8fywNxUpBmas764zuDLIU53vIBFVQcopaKa9Um9V6UvkzAWPl454bW/YPzOzFwF+6yO/WeHJi0jD66huQY4Bz5pZ9OBSpJkNAIf3qnw/eHwf8LuHXnsv5eWUH3POLVx7/pr71O99+sAqYpFisKq4vx1tb7uaK5VXFSFYUokA7q1LUUT2VwHNLNh+ZnDg0XjzOmb5/Rwl98bj3pLTfq7S3qpsEczAlf9ULBKjyaLEIk3EojFikSaaIjGaok0km9pJNrXR3NRCIpqgJdZCIhY8RhMkoomGPVOoRIndwm7YMeSMutx7ce7RYx3PcQryy5T3XP8c8PkDw//54RGdcyNmdgP4mJn9S+fcHICZXQQ+BvyJc24ByvdW/tvPf7y8o8EVKZQKFEoFSq5IoRR874KdEPvPSxRcMF6pyG4xz25pl3xxb9U1H7zHXsGVt8GVt8XtrQQeXhUsb3PbK9TydrXyY8xixKNNxINCS0TjNEXjNEWaiEebaIrGiUfjxCPlx0SsPJ6uXflozbHmR48kErLjfJJ/Bfg7wO+Y2buBe8AHKC8lpo4Y/5OUDwF62cz+YzDsH1Fe7PvFgyNGI1GikejjJX8bzrn9JdCDS4N7RXnQfjEesbdWO45EGtsjC9I5lzGzHwN+lfJSJJRPJfwrwDeOGP87ZvYB4N8EXw74DvAx51xVDoQzM6IWJUrly1dEGsex1gWdc1PAR4946apn/O8CP/74sUREwteYewlERI5BBSki4qGCFBHxUEGKiHioIEVEPFSQIiIeKkgREQ8VpIiIhwpSRMRDBSki4qGCFBHxUEGKiHioIEVEPFSQIiIeKkgREQ8VpIiIhwpSRMRDBSki4qGCFBHxUEGKiHioIEVEPFSQIiIeKkgREQ8VpIiIhwpSRMRDBSki4qGCFBHxUEGKiHioIEVEPFSQIiIeVS1IM3vRzG6Y2Y1UKlXNSYuInFhVC9I595Jz7rpz7npPT081Jy0icmJaxRYR8VBBioh4qCBFRDxUkCIiHipIEREPFaSIiIcKUkTEQwUpIuKhghQR8VBBioh4qCBFRDxUkCIiHipIEREPFaSIiIcKUkTEQwUpIuKhghQR8VBBioh4qCBFRDxUkCIiHipIEREPFaSIiIcKUkTEQwUpIuKhghQR8VBBioh4qCBFRDxUkCIiHipIEREPFaSIiIcKUkTEQwUpIuKhghQR8VBBioh4VLUgzexFM7thZjdSqVQ1Jy0icmJVLUjn3EvOuevOues9PT3VnLSIyIlpFVtExEMFKSLioYIUEfFQQYqIeKggRUQ8VJAiIh4qSBERDxWkiIiHClJExEMFKSLioYIUEfFQQYqIeKggRUQ8VJAiIh4qSBERDxWkiIiHClJExEMFKSLiYc65cCZstgYMhTLx4+kBavnGObWcr5azgfI9qVrP90POufZKvFGsEm/ymIacc9dDnP7bMrMbyvd4ajkbKN+Tqod8lXovrWKLiHioIEVEPMIsyJdCnPZxKN/jq+VsoHxPqmHyhbaTRkSk1mkVW0TEo+oFaWYfMrMhMxsxs09Ve/pBhitm9k0zu29m98zsk8Hwc2b2/8zsYfDYHQw3M/uNIPPrZvZClXJGzeyWmX01+P6amb0S5PhDM4sHwxPB9yPB61erkK3LzL5kZoNm9sDM3ldL88/M/knws71rZn9gZs1hzj8z+z0zWzKzuweGnXh+mdkngvEfmtknTjHbvw1+tq+b2f80s64Dr306yDZkZj95YPipfLaPynfgtV80M2dmPcH3lZ13zrmqfQFRYBR4GogDt4Hnq5khyDEAvBA8bweGgeeBXwE+FQz/FPDZ4PmHgf8DGPBe4JUq5fynwH8Dvhp8/0fAx4PnnwP+QfD8HwKfC55/HPjDKmT7L8DfD57Hga5amX/AJWAcaDkw334+zPkHvB94Abh7YNiJ5hdwDhgLHruD592nlO2DQCx4/tkD2Z4PPrcJ4FrweY6e5mf7qHzB8CvA14BJoOc05t2pfoiO+Iu+D/jage8/DXy6mhk8ub4M/DXKB64PBMMGKB+rCfDbwM8cGH9/vFPMdBn4BvBXga8GP/DUgV/a/XkZ/JK8L3geC8azU8zWGRSQHRpeE/OPckFOBx+GWDD/fjLs+QdcPVRCJ5pfwM8Av31g+JvGq2S2Q6/9TeALwfM3fWb35t1pf7aPygd8CfgLwARvFGRF5121V7H3fnH3zATDQhOsTr0LeAXoc87NBy8tAH3B8zBy/zrwS0Ap+P48kHXOFY7IsJ8veH01GP+0XAOWgc8HmwB+18zaqJH555ybBf4dMAXMU54fN6md+bfnpPMrrM/PL1BeKquZbGb2EWDWOXf70EsVzdfQO2nMLAn8D+AfO+dyB19z5X9mQtnFb2Y/BSw5526GMf1jiFFe5fkt59y7gA3Kq4j7Qp5/3cBHKBf5RaAN+FAYWY4rzPn1dszsM0AB+ELYWfaYWSvwL4BfPu1pVbsgZylvN9hzORhWdWbWRLkcv+Cc++Ng8KKZDQSvDwBLwfBq5/5R4KfNbAL475RXs/8D0GVme6eHHsywny94vRNYOcV8M8CMc+6V4PsvUS7MWpl/PwGMO+eWnXN54I8pz9NamX97Tjq/qjofzezngZ8CfjYo8FrJ9k7K//jdDj4jl4HXzKy/0vmqXZDfB54N9ibGKW8Q/0qVM2BmBvxn4IFz7lcPvPQVYG/v1icob5vcG/5zwR6y9wKrB1aNKs4592nn3GXn3FXK8+hPnHM/C3wT+Kgn317ujwbjn9rSiHNuAZg2sx8KBv04cJ8amX+UV63fa2atwc96L19NzL8DTjq/vgZ80My6g6XkDwbDKs7MPkR5E89PO+c2D2X+eLDn/xrwLPAqVfxsO+fuOOcuOOeuBp+RGco7XReo9Lyr1EbUE2xs/TDlvcajwGeqPf0gw1+mvDrzOvCD4OvDlLc7fQN4CHwdOBeMb8B/CjLfAa5XMesHeGMv9tOUfxlHgC8CiWB4c/D9SPD601XI9SPAjWAe/i/KewZrZv4B/xoYBO4Cv095r2to8w/4A8rbQ/PBB/rvPc78orw9cCT4+runmG2E8ja7vc/H5w6M/5kg2xDw1w8MP5XP9lH5Dr0+wRs7aSo673QmjYiIR0PvpBEReTsqSBERDxWkiIiHClJExEMFKSLioYIUEfFQQYqIeKggRUQ8/j8+dem88U4R+gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "style = SchedStyle(d_ch_color=['g','g'], use_table=False, figsize=(5, 5), remove_spacing=True, dpi=800)\n",
    "tomo_schedules_opt[22][2].filter(channels=[backend.configuration().drive(0)]).draw(style=style, plot_range=(0,1400), scaling=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run tomo experiments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "import itertools\n",
    "flat_schedules_std = list(itertools.chain(*tomo_schedules_std))\n",
    "flat_schedules_opt = list(itertools.chain(*tomo_schedules_opt))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/pranavgokhale/Developer/qiskit/qiskit-terra/qiskit/providers/models/pulsedefaults.py:155: UserWarning: `qubit_freq_est` and `meas_freq_est` now have units of Hertz(Hz) rather than gigahertz(GHz).\n",
      "  warnings.warn('`qubit_freq_est` and `meas_freq_est` now have units of '\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "job_std is 5e375652804aa70018a115a8\n",
      "job_opt is 5e375658fcf365001889c107\n"
     ]
    }
   ],
   "source": [
    "job_std = q.execute(flat_schedules_std, backend=backend, shots=shots)\n",
    "job_opt = q.execute(flat_schedules_opt, backend=backend, shots=shots)\n",
    "print('job_std is %s' % job_std.job_id())\n",
    "print('job_opt is %s' % job_opt.job_id())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "job_std.status()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job_opt.status()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run a mitigation cal job"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "import qiskit.ignis.mitigation.measurement as mit\n",
    "meas_qcs, meas_labels = mit.complete_meas_cal(qubit_list=[0])\n",
    "meas_qcs_transpiled = q.transpile(meas_qcs, backend, basis_gates=['x'])\n",
    "meas_schedules = q.schedule(meas_qcs_transpiled, backend=backend, inst_map=circ_inst_map)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "job_mit is 5e37565fee0f3b0019d3d1f6\n"
     ]
    }
   ],
   "source": [
    "job_mit = q.execute(meas_schedules, backend)\n",
    "print('job_mit is %s' % job_mit.job_id())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "job_mit.status()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job_std = backend.retrieve_job('5e375652804aa70018a115a8')\n",
    "job_opt = backend.retrieve_job('5e375658fcf365001889c107')\n",
    "job_mit = backend.retrieve_job('5e37565fee0f3b0019d3d1f6')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-----\n",
    "### Fit tomography data to get statevector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit.result import marginal_counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "result_std = job_std.result()\n",
    "result_opt = job_opt.result()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.96679688, 0.09960938],\n",
       "       [0.03320312, 0.90039062]])"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# mitigate results\n",
    "result_mit = marginal_counts(job_mit.result(), [0])\n",
    "readout_fitter = mit.CompleteMeasFitter(result_mit, meas_labels)\n",
    "readout_fitter.cal_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "result_std_marg = marginal_counts(result_std, [0])\n",
    "result_opt_marg = marginal_counts(result_opt, [0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "from copy import deepcopy\n",
    "result_cpy_std = deepcopy(result_std_marg)\n",
    "result_cpy_opt = deepcopy(result_opt_marg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "standard_statevectors = []\n",
    "optimized_statevectors = []\n",
    "\n",
    "for i, theta in enumerate(thetas):\n",
    "    result_std_mit = readout_fitter.filter.apply(result_cpy_std)\n",
    "    result_opt_mit = readout_fitter.filter.apply(result_cpy_opt)\n",
    "    \n",
    "    # remove results for other thetas\n",
    "    # TODO: allow selecting subset of result\n",
    "    single_qst_size = len(tomo_schedules_std[0])\n",
    "    del result_std_mit.results[:i*single_qst_size]\n",
    "    del result_std_mit.results[single_qst_size:]\n",
    "    del result_opt_mit.results[:i*single_qst_size]\n",
    "    del result_opt_mit.results[single_qst_size:]\n",
    "    assert len(result_std_mit.results) == single_qst_size\n",
    "    assert len(result_opt_mit.results) == single_qst_size\n",
    "\n",
    "    qst_fitter = StateTomographyFitter(result_std_mit, tomo_circuits_std[i])\n",
    "    sv = qst_fitter.fit(method='cvx')\n",
    "    standard_statevectors.append(sv)\n",
    "    \n",
    "    qst_fitter = StateTomographyFitter(result_opt_mit, tomo_circuits_opt[i])\n",
    "    sv = qst_fitter.fit('cvx')\n",
    "    optimized_statevectors.append(sv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### plot state fidelities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "# what are the ideal statevectors\n",
    "from qiskit.quantum_info.states import Statevector\n",
    "ideal_statevectors = []\n",
    "for theta in thetas:\n",
    "    circuit = QuantumCircuit(1)\n",
    "    circuit.rx(theta, 0)\n",
    "    ideal = Statevector.from_instruction(circuit).data\n",
    "    ideal_statevectors.append(ideal)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "standard_fidelities = []\n",
    "for sv1, sv2 in zip(standard_statevectors, ideal_statevectors):\n",
    "    fidelity = state_fidelity(sv1, sv2)\n",
    "    standard_fidelities.append(fidelity)\n",
    "    \n",
    "optimized_fidelities = []\n",
    "for sv1, sv2 in zip(optimized_statevectors, ideal_statevectors):\n",
    "    fidelity = state_fidelity(np.abs(sv1), np.abs(sv2))\n",
    "    optimized_fidelities.append(fidelity)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x11b8b75c0>"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEGCAYAAABy53LJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOydeXhU1f2H35M9rElIgJCQBZJAQNkhiYosgrJErFtdarWt1tra2tbaWkVttSpotdSlbq2t1f7UWpVFBDGyqCgEgrKIZEhIAiRhCSEJIWTP+f1x5iaTZGZyJzN3JpD7Ps88IXfuPfcAyXzu+S6fI6SUmJiYmJiY6MXP1xMwMTExMTm7MIXDxMTExMQlTOEwMTExMXEJUzhMTExMTFzCFA4TExMTE5cI8PUEvEFkZKRMSEjw9TRMTExMzip27NhxQkoZ1fF4rxCOhIQEcnJyfD0NExMTk7MKIcRBe8fNUJWJiYmJiUuYwmFiYmJi4hKmcJiYmJiYuIQpHCYmJiYmLmGocAgh/imEOC6E+MbB+0II8awQIl8IsVsIMcnmvVuEEHnW1y02xycLIfZYr3lWCCGM/DuYmJiYmLTH6BXHa8A8J+/PB5Ktr9uBFwGEEBHAH4A0YBrwByFEuPWaF4Ef21znbHwTExMTEw9jqHBIKT8DTjo55QrgdanYCoQJIaKBy4AsKeVJKWUFkAXMs743QEq5VSpb39eB7xj5dzAxMTExaY+vcxwxwGGb74utx5wdL7Zz3Bg++AD+/W/DhjcxMfEMJ2tPUlBR4Otp9Bp8LRyGIYS4XQiRI4TIKSsrc30AKeHll+GHP4TXXvP4/ExMmluaqayr9PU0zlqklGws3MiN793IsKeHMe7FcZxpPOPrafUKfC0cJcBwm+9jrcecHY+1c7wTUspXpJRTpJRToqI6dcx3jRDwv//BnDnwox/Bv/7l+hgmJk5YsnkJw5cNJ6fUdDVwhWOnj/HkF08y6vlRzH59Nmvy1nDJiEuoaaxhy+Etvp5er8DXwrEKuNlaXZUOVEkpjwDrgEuFEOHWpPilwDrre6eEEOnWaqqbgZWGzS40FFauhLlz4dZbTfEw8SgbCjdwuuE0C/5vAXnleb6eTo+mRbbw8YGPueada4hdFsu9n9zLkH5D+Pd3/k3pb0p5++q38Rf+bCza6Oup9goM9aoSQrwFzAQihRDFqEqpQAAp5UvAGmABkA+cAX5ofe+kEOJPwHbrUI9IKbUk+89Q1VqhwFrryzhCQ2HFCrjySiUeUqoViImJGzS3NLO9dDvzk+azrWQbl/3nMr689UuG9hvq66n1KEpOlfCvnf/i1a9fpaiyiEGhg7hr2l3cNuk2UqNS2507ZdgUNhVt8sq8lm1ZRtzAOK4ec7VX7tfjkFKe86/JkydLt6mtlfKyy6QUQspXX3V/PF9z7JiUy5ZJ2dLi65n0SvYc2yP5I/LfO/8ttx7eKvs81kdOfGmirKqr8vXUfE5jc6NclbtKXv7m5dLvYT/JH5Gz/z1bvrXnLVnXWOfwunuz7pWBjwTK0/WnDZ3fmYYzMvhPwTLyyUjD7+VrgBxp5zPV16Gqs4eQELXyuPRStfJ49VVfz8g9fvEL+PWvYe9eX8+kV5JdnA1AWkwaabFpvHvtu+w+tpur/nsVDc0NPp6d7zhw8gAjnhnBorcXsa1kG7+74Hfk/SKP9Tev5/rzric4INjhtbMSZtHY0siXh780dI6bD22mvrmeE2dO8MqOVwy9V0/FFA5X0MRj3jy47Tb4xz98PaPu8cUX8M476s9FRT6dSm8luySbsJAwkgclAzA/eT6vLnqV9YXr+cGKH9AiW3w8Q9+wIncFh08d5u2r3+bwrw+zZM4SkiKSdF17YdyFBPgFGJ7nyCrIItAvkIzYDJ7a8hR1TXWG3q8nYgqHq4SEwPLlMH8+/PjHZ594tLSolUa4tRH/oF27fRODyS7JZlrMNPxE26/gLRNuYcklS3jrm7f4zbrfoCIFvQtLuYVBoYO47rzrCPQPdOnafkH9mDpsquF5jqyCLC4YfgGPzn6U0upSXtv5mqH364mYwtEdQkLg/fdhwQIlHn//u69npJ8334Tt2+Evf4HgYFM4fMDphtN8c/wb0mLSOr1374X3cte0u/hr9l956sunfDA732IptzAqclS3r5+ZMJPtpds53XDag7Nq43jNcXYe3cncEXOZlTCLjNgMlm5eSmNzoyH366mYwtFdbMXj9tvhlbMg1nnmDNx3H0yeDDffDPHxZqjKB+wo3UGLbLErHEIIls1bxnfHfpffffI73tj1hg9m6DssJyyMGtR94ZiVMIumlia+OPSFB2fVxvqC9QDMHTkXIQSLpy/mYNVB/m/P/xlyP4D6pnr+vuPvnKo/Zdg9XMUUDncIDlbisXAh/OQnqtO8J/PUU1BcDMuWgZ+fEg5zxeF1sktUYnxazDS77/sJP17/zuvMTpzNj1b9iI/yP/Lm9HxGVV0Vx2qOuSUcFwy/wNA8R1ZBFuEh4UyOngzAguQFTBg6gSWbl9Dc0mzIPZ/84kluX307D2962JDxu4MpHO4SHAzvvafE4447eq54lJTAE0/ANdfA9OnqmCkcPiG7JJsR4SOI6uvY0SA4IJjl1y1nbNRYrnnnGraXbHd47rmCpdwC4Faoqm9QX6bFTDMkzyGlJKsgi9mJs/H38wfUCvGB6Q+wv3w/7377rsfvmX8yn8c+f4xg6c+LOS9y9PRRj9+jO5jC4Qk08cjMVOLx0ku+nlFnFi+GpiYlHhrx8XDsGNTW+m5evZDs4my7YaqODAgewNrvrSWqbxQL3jz3u8stJ6zC4caKA1S4Kqc0h+r6ak9MqxVLuYXiU8XMHTG33fErU68kNTKVRz9/1KPVcFJK7lxzJ0HSj49fa6a+qY4/f/Fnj43vDqZweIrgYHj3Xbj8cvjpT1USuqewY4dy+f3Vr2DEiLbjCQnq66FDuod6Oedlxr80vldW/HiCklMllFSX6BIOgOj+0ay7aR0Al/3nsh7zxGkElnIL/sKfkREj3RpnZsJMmmUzmw9t9tDMFFkHsgCV37DFT/ixePpivjn+DR9YPvDY/d7Z+w4fH/iYx1pmcvFBuMl/Ii/mvMix08c8do/uYgqHJwkOVsaIY8f2nGS5lKr8NipKrTpsiY9XX10IV3126DN2H9tN2ZluOA6btOY30mL1CQdAyqAU1ty4hmM1x5j/f/N7VJLUk1jKLSSGJxLkH+TWOBcMv4BAv0CP5zmyCrIYET6CEeEjOr133XnXMSJ8BI9+/qhHHqqq6qr41bpfMTl6Mj/LVuM9UHE+9c31/PlL3686TOHwNMHBcMklquS1qcnXs1HJ+88/hz/9CQYMaP9eN4SjqLIIULFXE9fJLs4m0C+QCUMnuHTd1JipvPfd9/jm+Ddc+d8rqW+qN2iGvsPdiiqNPoF9SItN82ieo7G5kU1FmzqFqTQC/AK476L7yCnNIasgy+37PbDhAY6dPsZLC1/Ef5vKbyWX1PG987/HC9tf4HjNcbfv4Q6mcBhBRoYqfd2927fzqK+H3/4WzjtP2aR0ZNgwCAhwqSS3sKIQMIWju2SXZDNh6ARCAkJcvnZe0jz+ueifbCjcwF1r7zJgdr6jRbaQdzLPI8IBKs+x48gOquqqPDLetpJtVDdUOxQOgJvH30zsgFge/exRt+6VU5rDCzkvcOfUO5lSGw7l5eqNI0d44OIH1KrDx7kOUziMICNDfd3i470BnnkGCgtVs1+AHSPkgACIjdW94qhtrOXI6SOAKRzdobmlmZzSHN35DXt8f/z3uWX8Lby99+1zypbkUNUh6prq3KqosmVmwkxaZIvH8hxZBVkIBLMTZzs8J8g/iHsvvJfPD33OZwc/69Z9mluauWP1HQzuO5hHZz8K2Sq0yahRcPQoKYNSuPH8G/nb9r/5dNVhCocRxMVBdLRvheP4cXj0UVUmPNfxU5IrJbkHq9rOM4XDdfaW7aWmscal/IY9ZsTP4FT9qdYqpHMBT1VUaWTEZhDkH+SxcFVWQRZThk0hPDTc6Xm3TryVIX2HdHvV8WLOi+w4soNlly1jYMhA2LYN+vRRG8odUQ9tD0xXqw5fOguYwmEEQqhVhy+F46GHVJntU138cLnQPa7lN/oH9TeFoxvYOuK6Q3psOgBbi7e6Paeegid6OGwJDQwlPTbdIwnyqroqsouznYapbO/7m4zfkFWQxbaSbS7dp7S6lPvX38/cEXO5bux16mB2NkyZoiID1dVQU8OoyFHccN4NPl11mMJhFBkZUFCgnvy9zZ49yj/rpz+F0aOdn5uQAKWl0Ni1146W35iVOMsUjm6QXZJNRGiEbrdXR4yKHMXA4IGtFVrnApYTFgYED2BI3yEeG3NWwiy+Pvq12/u6byraRLNs7lSG64g7ptxBRGgEj33+mEv3uXvd3TQ0N/C3BX9DCKFylF9/DdOmqQgGwFFVjv3AxQ9Q11TH018+7dI9PIUpHEbhqzyHlPCb38DAgfCHP3R9fny8cswtLu7y1MLKQoL8g5geN52KugpO1p7s8hqTNjRHXLXrcffxE35Mi5l2zq04Rg0a5fa/jS1anuPzg5+7NU5WQRZ9AvuQEZuh6/z+wf35VdqvWGVZxa6ju3Rd8/GBj/nv3v9y//T7W6322bULGhogLa1NOKzhqtGRo7n+vOt5fvvzlNV4vzTeFA6jmDwZAgO9Lxxr1kBWlhKNQYO6Pt+FktzCykLiB8aTMigFMPMcrlBdX83e43vdDlNppMWksef4Hmoaajwynq9x1xXXHumx6QT7B7ud58gqyGJG/Aynm0h15Bdpv2BA8AAe3/x4l+fWNtbysw9/RsqgFO698N62N7TEeFoaDLVuKWwVDoAHL36Q2sZan+Q6TOEwipAQmDjRu8LR2KhWGykp8LOf6btG6x7XkecoqiwiMTyxNdRyrltgeJKc0hwk0mPCkR6bTotsIac0xyPj+ZKahhqKTxV7LDGuERIQQsbwDLfyHIeqDrG/fL+u/IYtYSFh3Dn1Tv6393/knsh1eu7SzUs5UHGAFxa80F6ctm1TK43Y2E4rDlCrjhvOV7kOb686TOEwkowM1QioI3/gEV56CSwWlRAP1LkJzvDhKpmvZ8VRUUhiWCIjwkcgEOaKwwW6csR1FW2ccyHPsb98P+C5iipbZiXMYufRnVTUVnTrekc2I3r4dfqvCQkIYenmpQ7PsZywsPSLpXzv/O9xyYhL2r+Zna3yG0Ko6EFAQGuOQ+OB6Q9wpvEMT2/xbq7DFA4jychQlU3eaASsqIA//lF1rWdm6r8uKEg9zXQhHNX11ZTXlpMYlkhIQAjDBw4nv8IUDr1kl2STFJHEoD46woc6iOobxcjwkedEnsPTFVW2zEyYiUR2u68iqyCL6H7RjI0a6/K1UX2juGPKHfxn939aC0tskVLyszU/IzQglKcv7fDBf/Ik5OWpMBWobRCGDm234gBIjUpVuY5tz3PizAmX59hdTOEwEm8myB95BCorVbOfqwnGhIQuQ1VaKW5CWAIASRFJ5opDJ1JK3Y64rpAWm3ZOrDgsJywIBMkRyR4fOy0mjZCAkG7lOVpkC+sL1zNnxJxuJ+3vueAe/P38eeKLJzq99+aeN9lQuIEllyxhSL8O1WTbrKW8aTY/M3aEA1Su40zjGa9WWBkqHEKIeUIIixAiXwjxezvvxwsh1gshdgshNgkhYm3ee0II8Y31dZ3N8deEEIVCiJ3Wl2umP95k+HBl62G0cOzfD88/r2xFxo1z/XodTYCFleqJKTE8EYCkcFM49FJ8qpgjp494XDjSY9IprS6l+FTXFXE9GUu5hbiBcYQGhnp87OCAYC4YfkG38hw7j+7kxJkTLuc3bBnWfxg/mvAj/rXzX+3+nypqK7j747uZFjON2yff3vnC7Gz1ADhlStux6Gi7wpEalcp1513Hc9ue89qqwzDhEEL4A38D5gNjgBuEEGM6nPYU8LqUchzwCLDEeu1CYBIwAUgD7hFC2Dr0/VZKOcH62mnU38FtvNUI+NvfQmioMjLsDvHxcPgwNDvewUxbaieGWYUjIokTZ064XSPfG+iOI64etPHO9nCVERVVtsyMn8nuY7tdLh/X8htzRsxx6/73XnQvzS3N7aqfFm9YzIkzJ3hp4Uutm0K1Y9s2GDOmvTFpdHSnHIeGtur4y5a/uDVXvRi54pgG5EspC6SUDcDbwBUdzhkDbLD+eaPN+2OAz6SUTVLKGmA3MM/AuRpHRobyizpmkIf+V1/BqlVw//0wpJvNU/HxKoFv52lGo7CykD6BfYjsEwnQWll14OQBfff4z3/g2We7N7+znOzibIL8gxg/ZLxHxx0/ZDxB/kGtHelnI1JK9pfvNyQxrjErcRYSyadFn7p0XVZBFucNPo/o/tFu3T8hLIHvj/8+r+x4heM1x8kuzualnJe4a9pdTIye2PkCKdsS47ZER0NZmV3X7TFRY/ju2O/y3LbnKD9T7tZ89WCkcMQAh22+L7Yes2UXcJX1z1cC/YUQg6zH5wkh+gghIoFZwHCb6x6zhreWCSHsFlcLIW4XQuQIIXLKyny4d4TReY6VK1Xi7Lbbuj+GVpLrJFxVVFlEYlhia6xXEw7d4apnn4Unn+z+HM9iskuymTh0okt9AHoIDghmUvQktpacvSuO0upSTjecNlQ4pg6bSmhAqEt5jtrGWjYf2uxWmMqW+y66jzrrDn53fHgHw/oP45FZj9g/uaBAOeKmdVihDh2qRMXBQ+iDFz9ITUONV1Ydvk6O3wPMEEJ8DcwASoBmKeXHwBrgS+AtYAugxVHuA0YDU4EI4N6OgwJIKV+RUk6RUk6JinK8t7PhTJpkbCPg6tVKnCIjuz+GjibAwsrC1vwG0LqZjS7hkBJyc9W+56fOzU2IHNHU0sSOIzs8nt/QSItJY0fpDhqbvVTy7WGMrKjSCA4I5sK4C13Kc3x+6HPqm+s9Jhwpg1K47rzreGrLU+w8upNn5j1D/+D+9k+2bfyzxU4vhy1jB4/l2rHX8uy2Zw1fdRgpHCW0XyXEWo+1IqUslVJeJaWcCCy2Hqu0fn3MmsOYCwhgv/X4EamoB/6FCon1XEJClHgYIRwlJSpU5Ur5rT3i4tRXB8IhpWzt4dDoG9SXYf2H6SvJLS1VBm2gBKQX8c3xbzjTeMbj+Q2N9Nh0aptq2XN8jyHj1zbWMnzZcF7b+Zoh43vaFdcRM+Nnsuf4Ht3J408KPiHIP4iL4y/22Bzuv+h+ABYkL+Cq1Kscn5idrRxxzzuv/fEOflX28Naqw0jh2A4kCyEShRBBwPXAKtsThBCRQghtDvcB/7Qe97eGrBBCjAPGAR9bv4+2fhXAd4BvDPw7eIaMDMjJ8Xwj4Jo16uvll7s3Tt++amtZByW5FXUVVDdUt5biauguybUVi337uj/PsxBPOeI6QhvXqDzHpqJNFJ8q5sO8Dw0Z31JuoU9gH2IGdIxie5ZZibMAdOc5sgqyuGD4BfQN6uuxOZw/5Hy23LqFt65+y3l577ZtyrKo4x46Xaw4AM4bfB7XjLnG8FyHYcIhpWwCfg6sA/YB70gp9wohHhFCLLKeNhOwCCH2A0MAzU4yEPhcCPEt8Apwk3U8gP8TQuwB9gCRgHvbbXkDrRFwlz7DM92sXq3yE2M6Fqt1AycluR0rqjSSI5L12Y70ZuEoySayT6Tdfao9QUJYAoP7Djasn2NNnno4Mapyy1JuIWVQCn7C2Kj5lGFT6BPYR1ee43jNcXYe3emxMJUt6bHpDAge4PiEhoY2R9yOaMUvToQD4KEZD1HdUM2yrcvcmKlz7GwL5zmklGtQuQrbYw/Z/Pld4F0719WhKqvsjel4C66eim2C3LYu2x1qa+GTT+BHP3K94c8e8fGwd6/dtzr2cGgkRSRxrOYY1fXVjuO1oISjf3/V19ILhcMTjriOEEKQFpNmyAe7lJI1+WvwE34UnyqmtLqUYf2HefQelhMWj9mwOCPIP4iL4i7SledYX7AewBDh6JJdu5Sdesf8BiiXh0GDuhQObdXxbPaz3J1xNxGhER6fpq+T472D4cMhJsazeY5Nm9S+5u7mNzQSEtSKQ8pOb2krDnuhKoADFV2U5Obmqn1BUlN7lXCcqj/FvrJ9hoWpNNJj07GUW7rtx+SIvJN5FFQU8L3zvwd4PhxW11RHUWWR4fkNjZnxM9lbtrfLzY+yCrIIDwlnUvQkr8yrHY4S4xpOejlseehi66pjizGrDlM4vIWnGwE/+EDlJmbM8Mx48fFqFWOndLmosoiwkDDCQsLaHdddkmsrHAcOqCeqXsD2ku0edcR1hDa+qzvOdYUWprp/+v0E+gV6fFWTfzIfiTS0osoWPXkOKSVZBVnMTpxtvzHPaLKzVdnt8OH233dgO9KR84ecz9WpV/NM9jOG7JtjCoe3yMhQyWcdTwtdIqXKb8ydq6q2PIGTktzCysJO+Q2AkeEjgS6Eo7pabRI1erR6tbRAfu+wKvG0I64jpsZMRSA8nudYm7+W0ZGjGR05mglDJ3h8fG9VVGlMjp5M38C+TvMclnILxaeKfROmApUYT0tzHH52YDtij4dmPMT4oeMNsSExhcNbeLIRcM8eZRHiqTAVtAmHncqqjj0cGv2D+zOk7xDnwrFfWWa3rjig14SrskuySRmUQnhouKH3GRA8gDFRYzy6IqhpqGFT0SYWJC0AVDgspzSH5hbHtjSuovVwaBuDGU2gfyDT46c7zXO4Y6PuNhUV6vfFXmJcQwtV2Qkpd2TckHF8/sPPDfn3NYXDW0yapJJbnhCO1avV1wUL3B9Lw8GKQ0pJUWURCQMT7F7WZUmuVlE1ejSMsj5Z9gLhMMoR1xFpMcopV+r4QNHDxqKNNDQ3sCB5Qev4NY017C2zX0DRHSzlFob1H+a8sMLDzIyfyb4T+zh22n73dVZBFiPCRxhWBecUe464HYmOVpVXFZ7NZ7mKKRzeIjjYc42Aq1er6qxo9zx02hEWpvYp7yAcx2qOUddUZ3fFATqFw98fRo5UOZn4+F4hHIeqDnGs5pjXhCM9Np2TtSc95li8Jm8NfQP7clHcRYAxhoqWExavhak0ZibMBODTg53zHI3NjWwq2uS7MJXmiDt1quNz7Gwh6wtM4fAmWiNgQ0P3xygrg61bPRum0rDTy+Goh0MjKSKJkuoSzjSesT9mbi6MGKGEE3pNZZVRjriO0O7jiTyElJI1eWuYM2JOq7/WyPCRDAod5LHKKimlcsX1snBMHjaZfkH92FjYOVyVXZJNdUO1b4UjNbW9I25HdDQBegNTOLxJRgbU1bnXCLh2rYpvGiEcdjZ0ctTDoaFVVhVUFNgfU6uo0khNVdvbtrS4O9seTXZxNsH+wYwb0o39UbrB2Kix9A3s65EP9twTuRysOsj8pPmtx4QQHt04SrPk91ZFlUaAXwDT46az6eCmTu9lHcjCT/gxO9EHrWJStiXGnaHDdsQbmMLhTTyRIF+9Wm0ONcmAGnM7Kw5t57/4gfF2L3FaktvcrJJ9HYWjtlbXHudnM9kl2UyKnkSQf5BX7ufv58/UmKkeccrVynDnJ89vdzwtJo1vy77lVL37RpWt5oZeXnGA2oc890QuR6rbP7VnFWQxZdgUw4sZ7FJYCCdOOE+Mg7ni6JXExqpXd4WjoQHWrYOFCz3TLd6R+HjlXlvZtjlTYUUhg/sOdujZ41Q4CgvVnG2FQ/vzOWx22NjcaKgjriPSY9LZeXQntY21bo2zNn8tY6PGEjcwrt3xtJg0JJLtJdvdGh9sSnG9vOIA+3mOqroqtpVs822YCrpecfTrpwwQTeHoZbjTCLh5s/pgNyJMBXZLch31cGiEhYQR2SfSvmeVJg5aGa7tn8/hPMee43uoa6rzWn5DIy02jaaWJr4++nW3x6iur+azg5+1VlPZovWjeCJcZSm3EOwf7HAlayQToycyIHhAuzzHxqKNNMtm3wpHaCicf77z84RwqZfDKEzh8DYZGSpM053/+NWrVZL5kks8Py+wu6FTYWVhJ6uRjiRFJNm3V9eEY5TNU2VkpHqdw8JhtCOuIzzhlLuhcAONLY3t8hsa4aHhjBo0ymPCkRSR5JPubHt5jqwDWfQN7EvG8AyvzwdQwmHPEdceOm1HjMQUDm/T3TyHlMpmZPZsVdZqBB16OZpbmjlUdcjpigOclOTm5sLgwRDRwWTtHK+syi7JJqpPVJeC62mi+0cTNzDOrTzHmrw19A/qz4VxF9p9Py1WGSq62y9iOWHsPuNdMSthFvvL91NaXQqo/MaMhBley0m1Q3PE7SpMpWGuOHohEyd2rxFw/35l1WFUmArUSiA0tDVUVVJdQlNLk8OKKo2k8CQOVx2mrqmu/RsdK6o0NOHwULNaTyO7JJu02DTDHHGdkRaT1u0Vh+aGO3fkXIcfoGkxaRyvOc7Bqu4XNzQ2N3Kg4oBPEuMaWp5jU9EmDlYeJO9knu/CVLt3K/+2rhLjGjr9qozEFA5vExyslqRffunadVq3+MKFnp+ThhBtLrl03cOhkRSRhES2nt+KM+E4edKuoeLZTmVdJbkncr0eptJIj03nYNVBjp52PZSxt2wvxaeK7YapbMcH98JhhZWFNLU0+VQ4JgydwMDggWws3MgnBZ8APrJRB/2JcY3oaJXrPOOgd8oLmMLhCzIyYMcO1xoBV69WibN4g5OJNiW5WimunhwHdKisOnECyssdCweck+EqreLIV8LhTp6jtQzXiXCcP/h8QgJC3Mpz+LKiSsPfz5+L4y9m08FNZBVkMaz/MMZEeWBDtO6Qna02aYqL6/pc6BG9HKZw+IKMDLU0/Vpn9UtlJXz+ubFhKo34+NZQVWFlIQLRqSyzI3aFw9ajqiPnsHBoH6hTY5zYRhjIpOhJBPgFdOuDfW3+WsYNGed0G9dA/0AmR092y3rElz0ctsxKmEX+yXxW71/NnBFzfBJaBJRwOHPE7UgP6OUwhcMXuJogX7dONdN5QzgSEtRKoaaGwspCYgbEtNpOOCIiNIKwkDD9whEbq2rRz8FejuySbEZHju60d4m3CA0MZfyQ8S5/sFfVVbH50OZWN1xnpMWk8dWRr2ho7p51juWEhag+Ub5ptLNBy3PUNAY9BVUAACAASURBVNb4LkylxxG3Iz3Ar8oUDl8QE6M2atErHKtXq8S13hioO9hUVhVWOO/h0BBCdC7Jzc1Ve4XYW377+SlBOcdWHN52xHVEemw620u3u2SB/knBJzS1NHXqFnc0fn1zPbuP7e7W/Czlvq2o0hg/dDzhIUq85oyY45tJbLc2U7ryu22uOHoxehsBm5thzRploe7vhZp3G+EoqizSXVKaHJHcecWRkuJ4zudgSW5RZRFlZ8p8LhxpMWmcbjjNt2Xf6r5mbf5aBgYPJCO26z6GVkPFbibIfWFuaA8/4ceC5AWkxaQxtN9Q30xCjyNuRyIjVb+HmePohWRkqM2YSkqcn7d1q6pA8kaYClqFo6Ewn+JTxbpWHKDyHEWVRW3hC0cVVRqpqervf/q0uzPuMXjbEdcRrjrlSilZm7+WuSPnEugf2OX5wwcMZ2i/od3qF6msq+R4zfEeIRwAry56lQ23bPDdBLKz1e/JwIH6r/HzU8l0c8XRC9Gb51i9Wj1dXHqp8XMCtQwODOTQ4W+QyC57ODSSIpJokS2qEquuTvlU2VqNdER77xzKc2QXZxMSEML5g7uwjTCY5IhkwkPCdec5dh/bTWl1qa78BqjQZHpserdWHD2hosqW4IBg+gT28c3NpWxLjLuKj3s5DBUOIcQ8IYRFCJEvhPi9nffjhRDrhRC7hRCbhBCxNu89IYT4xvq6zuZ4ohAi2zrmf4UQPmj1dM6x08e6NpqbOFH1dOgRjosvdu2JxB38/CAujqLjastXvaGqdpVVeXnKNr2rFQecU+EqzRFXz1O7kbhqga6V4c5Lmqf7HmkxaeSdzKP8TLlLc+spFVU9gqIifY649vCx7YhhwiGE8Af+BswHxgA3CCE6Fko/BbwupRwHPAIssV67EJgETADSgHuEENruJk8Ay6SUSUAFcKtRfwdXyS7O5rp3r2PYX4aR+VYmLdLJnhNBQaoR0JlwFBXBN994L0ylER9P4SnVy+FKqAqswuGsoqr1giS1kjpHhKOxuZGvjnzl8/yGRnpMOnuP79Vlgb42fy0Th04kur/+HSW1v+e2km0uzctywkKAX4Bvtmbtabja+GeLj21HjFxxTAPypZQFUsoG4G3gig7njAG0AONGm/fHAJ9JKZuklDXAbmCeUIXWs4F3ref9G/iOgX+HLmluaWb5vuVc9M+LSH81nXX568hMyWRD4QZeznnZ+cVaI2B9vf33tW5xbwtHQgKFjccJ8AsgdkBs1+cDUX2i6B/Uv71wpKQ4viAwUG0ne46EqnYf2019c32PEY60WGWBnlOa4/S8itoKvjz8pV03XGdMGTYFgXC5X8RSbmFE+Aifr8p6BNnZqvKwK0dcewwdCsePQ1OT5+elA13CIYQIF0KMFUKMEELoFZsY4LDN98XWY7bsAq6y/vlKoL8QYpD1+DwhRB8hRCQwCxgODAIqpZRNTsbU5ny7ECJHCJFTZoC1RU1DDc9ve55Rz4/iqneuoqS6hL9e9lcO//owK65bwdwRc/lt1m8723DYkpHRZnBmj9Wr1YdvcrLH5++U+HgKA2qIGzBct3tpa0muJhzx8apXwxnnUGVVT0mMa7RaoHeRh8gqyKJZNjvtFrdH/+D+nDf4vG4JhxmmsqI54gZ2Q0Sjo1WO5Phxz89LBw5FQAgxUAhxvxBiD7AVeBl4BzgohPifEGKWB+5/DzBDCPE1MAMoAZqllB8Da4AvgbeALYD+onRASvmKlHKKlHJKVFSUB6aqKK0u5b5P7mP4suH8Yu0viOobxf+u/R95v8jjl+m/pH9wf4QQ/GPRP/ATfty66lbHIStnCfLTp2HjRu+vNgDi4ykKg4TgIS5d1k44nIWpNFJTlXFjY2M3J9pzyC7JZnDfwT7ZX8IeEaERpAxK6bLyaW3+WsJDwrsleJqhol6n3OaWZvLK80zhAPXA+NVX3e/N8rHtiLPVw7uoFcN0KeUoKeVF1g/i4cBS4AohhLP8QglqlaARaz3WipSyVEp5lZRyIrDYeqzS+vUxKeUEKeVcQAD7gXIgTAgR4GhMo9h1dBe3rLiFhL8m8OSXTzI7cTZf/OgLtty6hWvGXEOAX3sf/biBcTx96dNsLNrIKztesT/osGGqQc6ecHzyifrh8oVwJCRQGA6J0rXu56SIJGVgZ9mnXziampR4nOVojX8+s62wg1b55OiDvUW2sDZvLZeOvLTTz68e0mLTqKirIO+knU287HCo6hD1zfU9pqLKp+zZ45ojbkd83AToUDiklHOllG9oH+Qd3tshpfyVlPJVJ2NvB5KtVVBBwPXAKtsThBCRNqGv+4B/Wo/7W0NWCCHGAeOAj6X6DdgIXGO95hZgpZ6/aHfQfrHmvD6HCS9P4L1v3+OnU35K3i/yePe773LB8AucXn/bpNuYM2IOv836bathYCccNQKuXq0qqS66yP2/iIuciRnMsX6QWBfi0nVJEUk0tTRxKKhWv3DAWR+uqqitwFJu6TH5DY20mDSO1RxzaIG+8+hOjtUcczm/oeGqU65ZUWWDO4lx8LntiNN8hbVc9jdCiGesr7uFELosHK15iJ8D64B9wDtSyr1CiEeEEIusp80ELEKI/cAQ4DHr8UDgcyHEt8ArwE02eY17gbuFEPmonIcz8XKLzDczWfDmAnJP5PLEnCc4/OvDPDP/Gd0VIUII/nH5PxAIbl11q/0nv4wMKC5WL42WFvjwQ5g3r3vxTzc52FeFjhI6PTI4JzlC5WLyI9AnHNo5Z7lwaJVFPSW/odGVU253ynBtSY1MpV9QP915jp7Ww+FTsrPVJmfddbvuqcIhhPgF8BpwBsgEcoAQ4AMhxE/1DC6lXCOlTJFSjpRSPmY99pCUcpX1z+9KKZOt59wmpay3Hq+TUo6xvtKllDttxiyQUk6TUiZJKa/VrjGCm8ffzBtXvkHBLwv43YW/65YpW3xYPE9d+hQbCjfYD1nZy3N89ZWKXfoiTAUUnlYilljaRS9KB1pLcvUKR79+yrPLiXC8tvM1Pj7wsUvz8DYf5X9EsH9w6xN4T2HckHGEBIQ4bARcm7+WKcOmMLjv4G6N7+/nz9RhU3U3GlrKLYSFhBHVx3M5x7MWVx1xOxIcrHbW7IE5jp8Bc6SULwKnrGGrx1F9FXd7ZXY+5vrzruemcTe5vZ3kjyf9mDkj5nBP1j2dQ1YTJqiSPFvhWL1aNeLN696ToLu0buBU5NqSY2i/ofSRAeQPDVKWCHpwYnZ4qv4Ut39wO4veWuRyv4C3kFKyPHc5c0bMoV9QP19Ppx2aBbq9FUH5mXK2Fm/V3S3uiLSYNHYd29V1wyttFVU9KQ/kEyoqwGJx37TUh70czoSjCRUKArCNsUQDZ38ZjBfRQlYAt626rX3Iyl4j4OrVaiUSGenlmSoKKwsJafFnaJ5rP5RCCJJqQsgb3kf/k1RqqqrCaulcefZR/kc0tjQSEhDClf+9kiPVvt0u0x67ju3iYNVBrhx9pa+nYpf02HS7FugfH/iYFtmiyw23q/GbWpr4+mjXe8v4ep/xHkOOtbemu4lxDR/ajjgTjp8DHwohXgdihRCvCCHWAMuBH3tlducQ8WHxPDX3KdYXrufvX/29/ZsZGSo8VV8PpaWqKdBHYSpQLq/xIgxRXOJyg1FSWRP5YS7sJZ6aqrbAtM3xWFlpWUlkn0g23LKByrpKrn7nauqbDItMdovl+5bjJ/y4fNTlvp6KXdJi0qhvrmfX0V3tjq/NX8ug0EFMHebehlN6nXJPN5ympLrETIxDW2LcFUdce/TEFYeU8lMp5VTgWeAO4CPgIWuJ7BfemuC5xO2Tb2d24mx+8/FvOFhpU+miNQJ+9ZVKioNPhaOwspDEkGhl6d6Ve68tp06RVFJHQdBp/XtBOKisamxu5MP9H5KZksmk6Em8dsVrbCnewp1r7tTdN+ANVlhWcOHwC7udJzCa1sonm3BVi2zho/yPuCzpMt0Nno4Y2m8ocQPjuuwX2V+uvM9M4aDNETfMzc2+NL8qH/w+OEuOJwBIKXOklO9bXzk27wtbU0KTrhFC8OoiVQT24w9+3PYBaJsgX71aVVqMHeujWWLdwClBfXPQfimnXXJzSToJDTRTfKrzCsIuDoTj04OfUlVfxXdGKUeZa8dey+Lpi3n161d5YfsL+udkIAUVBew+tpvvjPap641TYgfEEt0vul0Ce0fpDsrOlLmd39DQGgGdYVZUWXHHEbcj0dEqSlHpYvmjB3AWqvqzEOI9IcTNVruRwUKIOCHEbCHEn4AvACe+2Sb2SAhL4M9z/0xWQRb/+ErlPYiOVmKxcaNq/MvM7H61hZtU1VVRUVdBwlBrVZR1/3FdWIUDOuw/7oyoKFUd0kE4VuauJDQglLkj27b0fGTWI1yecjm//OiXbCrapH9eBrEidwVAjxaOVgt0mxXHmrw1CASXJV3mkXukx6ZzsOogx04fc3iOpdyCQLRW3vVaioqgrMwzwuHDklxnoaprgQeBUSiX289RzXa3ARZgtpQyyxuTPNewDVkdqjqkDmZkqNXGmTM+z28AJCZMUAdcXXGcUqEP3cIhRCfPKiklKy0rmTtybru9EvyEH/+56j+kDErh2v9d67ip0kssz13OuCHjerzTa1pMGvkn8zlx5gSg8hvTYqYR2cczxRet/SJO+jks5RYSwhIICXCtqfScY5u1OtDdxDj4tHvcaQOglPJbKeViKeVMq+3IRCnljVLK/0gp67w1yXMNP+HHPy7/By2ypa3KSgtX9ekDM2f6bG6FldZS3MGj1BONi8IRMziJYP9g/cIBnYRj59GdHD51mCtGdTRThgHBA1hx/Qoamxv5ztvfoaahRv99PMjxmuN8ceiLHltNZYuWwN5Wso2ymjK2lWzrdre4PSZFTyLAL8BpP4dZUWVFc8QdN879sXzoV9Wl060Q4k/WvTW07wcIIf5l7LTOfRLDE1tDVq9+/WqbcMydq36wfERrD0dYogqfuSgcfqNTGRkxkvwKF4Rj9Gi1oc0J9US8IncFAkFmiv2VV8qgFN6+5m32HN/DD1f+0CfJ8lWWVUhkjw5TaUwZNgU/4cfW4q2sO7AOiXTZDdcZoYGhjB8y3uGKQ0rJ/vL9ZmIclHBMmuQZR4ieuuKwEgBsE0KME0LMRXlQ7TB2Wr2Dn0z5CbMSZnH3urs5lBgBs2fDHXf4dE5FlUX0C+pHRGiEEg69OY7GRmVWOHo0yRHJrq84oHVvjpWWlVwY57xSaV7SPJZespT/ffs/lm5eqv9eHmJF7goSwhIYP2Q8fPaZ6oA/edLr89BDv6B+rRboa/PXEtUnisnDJnv0HmkxaWwv2W63mq6kuoSaxhpTOOrq3HPE7Uj//hAa2jOFQ0p5H/A7IBu1cdJCKeXzRk+sN+An/Hh10au0yBZu/+hO5Cef+KxbXKOwspDEsETV3ZuQAIcO2W3O63xhoRKP0aNJikjiwMkDzndAtMWmsqqosohdx3bZDVN15J4L7uGG825g8YbFrN6/Wt+9PEB1fTVZBVl8Z9R31L/TqlWqD+Xbb702B1dJj0lnW8k21uWvY17SPPx0b6ujj7TYNKobqtl3orMLgFlRZSUrS4nHpZd6ZjwhfNbLoSdUdTGql+MRYBPwnBBimMHz6jUkhify5NwnWXdgHf/8+p++no4SjnDrdrHx8aq/5JjjaplWbLaLTYpIorapVn+nd3y8enLat49VFmWgrEc4tH1PJkZP5Mb3biT3hHd2E/wo/yMamhu4MtWa39C6/g8fdnyRj0mLTaOyrpLy2nKP5jdax3diqGi64lp57z3leD17tufG9NHe43oeO54CrpVSLpFS3gj8nbbtXk08wB1T7mBmwkzu/vhuDlf57sNHSmnt4bARDtAXrtKEY9So1pJLvfs04OcHo0bBvn2stKwkNTKV5EH6dj3sE9iH5dctJzQwlCvevoLKOuNr2pfnLieyTyQXDr9QCesOa+S2BwuH1gjoJ/y4dKSHnnhtSB6UTHhIuN08h+WEhX5B/RjWvxc/bzY2qpXp5ZcrmyFP0VNXHECGlLJ1DS6lfB+YbtyUeh9ayKq2sZbntj3ns3mU15ZT01hDgtb8l2D9qidBnpurqrDCwtpccl3Mc5w88A2fFn2qa7VhS9zAON699l0KKgq48b0b9Xetd4OG5gY+zPuQRSmLVNf1rl1te8b3YOEYHTmaAcEDSI9NV/krD+Mn/JgWM82+cJRbSBmU0rvNDTdtUuaGV1/t2XF95FfV1X4c/kC4zfdBQojbgc+MnlhvY0T4CGYmzPRqrL4j7SqqoG3FoVc4rLmK4QOGE+gX6LJwrAktplk2c8Vo14QDYHr8dJ6f/zxr89fywIYHXL5eLxsLN3Kq/lRbNdVWawnqoEEqH9RD8RN+/HPRP3n60qcNu0daTBrfHP+G0w2n2x039xkH3n9fldp7Kr+hER0NVVVQ69oWCO7izHLkeuAksFsI8akQ4lKgAJgPfM9L8+tVZKZksu/EPg6cPOCT+7f2cGg5jn79VFd3V6EqKdvtM+7v58+I8BGuCcfo0awcBUODI5kW073mqJ9M+Ql3TL6DpV8s5a09b3VrjK5YkbuCvoF92zrat26FmBjV0NWDVxwAV4+52tA9Q9Jj02mRLeSUtjoTUdtYy8HKg71bOJqbYflyWLBAiYcn8VEvh7MVxwPAZCnlMODXwAfAT6WUV0opv/LK7HoZC5MXAvBh3oc+ub/Wid0aqgJ9vRxlZWoZbrN5U1JEkkvCUZ8yko+SYFHwOLcqfp6Z/wwXxV3Eratu7eQI6y4tsoWVlpXMT57f1gG9dSukp6ty3B4uHEajCb5tgjz/ZD4S2bsrqrZsUQUmV13l+bF9ZDvi7De0QUqZD2AVijwp5QfemVbvZGTESFIjU30WriqsKCQiNIIBwQPaDiYkdC0cWte3HeHQ25y3IbCY08FwRYV7LrNB/kG899336B/cn7s+usujzYHbSrZx5PSRVuNFjh+HggIlHHFxqoHRyyEDl3jnnTbLCwMY1GcQSRFJ7ZxyzYoqVDVVUBAsXOj5sX3UBOhMOAZb9xi/WwhxNxDW4XsTA8hMyWRT0Saq66u9fm+th6Md2orD2QewTSmuRlJEEjWNNRyr0VHKC6w48CF9GwWz97nvZDO472D+MOMPfHbwM9bmr3V7PI3l+5YT4BfAwhTrB4C2r4K24oCeu+rYvx9uvBEeftjQ22hOuZpgaz0cKYNSDL1vj0VKld+YOxcGDOj6fFfpgaGqvwP9bV4dvzcxgMyUTBpbGskq8L5/ZLseDo34eKipgfJyxxfm5qrYbWyby35yhCqn1ROuapEtrNq/inmnhxLy7f5uzb0jP570Y0aGj+T3n/zeI1VW2haxsxJmERZi3Udh61YICFAWEj1dOB58UMXanezv7gnSY9M5cvpIq62+pdxC7IBY+gb1NfS+PZYdO1TRhKerqTQiI8Hfv+esOKSUDzt7eXOSvYkLhl9AWEiY18NVLbKFg5UHSRiY0P4NPSW5ubmqD8Ov7cfJlZLc7SXbOXr6KFf0nQR5eS7vOtiJxkYCb7udR2NvZs/xPby55033xgP2ndhH3sm89qaGW7fC+PFKNHuycOzYocJUUVGq0KHGOGPIjk65vb6i6v331Qf7okXGjO/vD4MH9xzhMPENAX4BzE+az4d5H+q37PAAR08fpb653v6KA7oWDpswFaitcgP8AnQJx0rLSvyFPwtHzleNUgfcrCp76y147TW+u3w/E4dO5MGND7q95ay298aiUdYPgOZmlS9It1Ypaautnigc99+vyoWXLlWhE4vFsFuNHzqeYP9gthZvRUqpXHF7q3BIqfIbM2eqf3+j8EEToKHCIYSYJ4SwCCHyhRC/t/N+vBBivRBitxBik+2OgkKIJ4UQe4UQ+4QQzwpr95D1PIsQYqf11TP37HSDzJRMjtccb1fWaDSdejg0uuoeP3NGiUoH4QjwCyAhLEG3cFwcfzER51n3YHYnnNLcDEuWAOD3yXqWXrKEg1UHeSnnpe6PiRKOtJg0YgbEqAN798Lp023CERKinvx6mnBs2AAff6zEQ5urgZ5aQf5BTIyeSHZJNsdrjlNVX9V7K6q+/VblloyoprLFB7Yjzvo4fmn9emF3BrY2D/4N1fcxBrhBCDGmw2lPAa9LKcehvLCWWK+9ALgQGAecB0wFZthc9z3r3ucTpJTHuzO/noxmQufNcJXdUlyA8HDlwuloxZGXp56sOggHqHBVV7YjeeV5fFv2reoWH2X9gMl1w3Pq/ffV9fPmwdGjzK2LYXbibB79/FFO1Z/q1pDFp4rZXrq9vYW61vin2eGDClf1pCZAKeG++9Rq6Gc/g6QklZMx2IwxPSadHaU72Fu2F+jFFVXvv6+MCK80eM+WHrbi+KH1a3c9MKYB+VLKAillA/A20LEleAxtvlcbbd6XQAgQBAQDgYC+8pxzgIjQCC4cfqFXhUNr/uskHEI47+WwU1GlkRTedUnuSstKANUtPnAgDBvW/RWHlPDYY0qAXlD7kov161l6yVJOnDnB0192r2t6Za6aY6f8RmQkjLDZ/a+n9XIsX67CaQ8/rFZEQUGQnGx4gjwtNo3aplre/fZdoBe74r73nnqw0CqfjGLoUFUa3myc1U5HnAnHPiFEHjDKGkrSXnuEELt1jB0D2P4WFVuP2bIL0NZxVwL9hRCDpJRbUEJyxPpaJ6W0/Wn/lzVM9aAWwuqIEOJ2IUSOECKnrKxMx3R7FpkpmXx99GtKTpV45X6FFYUM7TeU0MDQzm8625cjN1eJS3JnU8KkiCRO1Z9q3bLUHistKxk3ZFybYHXYDdAlPvxQeUfddx8kJqo5ZWUxNWYq1465lqe3PO10X2xHLM9dzujI0e0/ALXGP9sfv54kHE1NsHixEvSbb247PmaM4SsOLUH+9jdvExIQQtzAOEPv1yM5cED9LBpVTWVLdLTa+uC494IvzqqqbkCZGeYDl9u8Mq1fPcE9wAwhxNeoUFQJ0CyESAJSgViU2MwWQmjGit+TUp5vndt04PsO5v+KlHKKlHJKVFSUh6brPbTd77zVRW63h0OjqxVHYqKyRe9AV5VVZTVlfHn4y7aGOlDCkZvrvG/EHlLCo4+qKrAbb1TH5sxR5nKNjTw6+1Hqmup49LNHXRq2oraCTUWb2q82KiuVuKV3sO+Ii4PqauUd5Gtef139Oz7+uApPaYwZozbcqnevWMAZCWEJDO47mIq6CpIjkj2+98dZwfvvq69Gh6nAJ70cXe05flRKOR711K/1b5RKKfXsJ1oCDLf5PtZ6zHb8UinlVVLKicBi67FK1Opjq5TytJTyNLAWyLC+X2L9Wg28iQqJnXOkRqaSGJbotXBVUWVR5zCVRkKC+rC094Fop6JKoyvhWL1/NS2ypb2pYWqq+vAtcXGltWGDasi79962bTnnzFGlp9nZpAxK4bZJt/HyjpcpqCjQPezq/atpls3t8xta93VH4dBKcn2d56irgz/8QflnfafD1rZjxqin0/2e6ZexhxCiddXRa8NU77+v+nsSHTyMeRIfdI/r2chpBpCHSnS/AOy3bu7UFduBZCFEohAiCLgeWNVh7EghWh9H7gO0nYwOoVYiAUKIQNRqZJ/1+0jrtYGo1c83OuZy1iGE2nP7k4JPqG001saiqaWJQ1WHnK84oPOqo6VFlXY6EI6EsAT8hJ9D4VhpWcnwAcOZOHRi20Gb3QBd4rHHVH7kBz9oOzZrluot+eQTAB6a8RABfgE8uPFB3cOusKwgpn8MU4ZNaTu4dasKUU2d2v7kntLL8cILakfCpUvbh9JACQd4LVzVKxPjxcXqZ8ToaioNH/hV6VlD/gW4VEo5Q0p5MXAZsKyri6SUTcDPgXXAPuAdKeVeIcQjQgitG2YmYBFC7AeGAI9Zj78LHAD2oPIgu6w+WcHAOmuOZSdqBfN3XX/Ts5DMlExqm2rZWLTR0PsUn1J25p16ODQcCcfhw8qbyYFwBAcEEzcwjvyKzsJxpvEMHx/4mEWjFrXfp6E7wvHll7BxI9xzj0oCa4SHw5QpastOYFj/Yfwq/Ve8uedNvj7ydZfD1jbW8lH+R1wx6or24ZatW2Hs2M4WEj1BOKqqVHjq0kuVcHYkJUWJqdGVVVYX3tGR9n82zmlWqJ4fr+Q3oMcKR6CUsrVjSEq5H1Xl1CVSyjVSyhQp5Ugp5WPWYw9JKVdZ//yulDLZes5tUsp66/FmKeVPpJSpUsoxUsq7rcdrpJSTpZTjpJRjpZS/lFJ6r5TAy8yIn0HfwL6Gh6sc9nBoOOoed1JRpZEckWx3xZF1IIvaptrOmzYNGaKqq1wpyX3sMVXhdPvtnd+bM0eFsE6pUtzfXfg7wkPCuW/9fV0Om1WQxZnGM21bxILKpWiJ8Y5ER6tOXl8Kx9NPK3uYxx+3/35ICIwcabhwzEyYyUsLX+KqVC89dfck3ntPPQA5+b3wKCEh6iGpp+Q4rOQIIf4hhJhpff0d8F5nWi8mOCCYS0deyur9qz3q8toRhz0cGoMHqx/OjsJhxxW3I47s1VdaVjIgeAAzEma0f0MI1yqrvvoK1qyBX/8a+trxQ5ozR5UpfvopAGEhYdw//X7WHVjHxkLnK7nlucsJCwljRrzNHPfvVxby9oQjIECFy3wlHMeOwV/+At/9Lkye7Pi81FTDhcPfz5+fTPkJfQI9vP9ET6esDD77zHurDQ0v93LoEY6fAt8Cd1lf31qPmXiBzJRMDp86zJ7jewy7R2FlIX7Cz3HZpBCqYqhjSW5urtroKTLS4dhJEUmcrD3JydqTrceaW5pZvX81C5MXEuRvZ/9lV4Tj8cfVCuXOO+2/f8EFquIrq8008ufTfk7sgFh+v/73DgW5qaWJDywfkJmSSaC/zQJba/yzJxzg2ybAxx5TifE//cn5eWPGKAFsbPTOvHoTK1eq3J+38hsaXt5CtkvhkFLWSyn/Yq1+z+2gQgAAIABJREFUukpKuUwLKZkYz4LkBQCGhqsKKwuJHRDb/gOyI/b25dAqqpzsJa1VVtnuarileAtlZ8oc7y2emqqenisqnE/8229VWOAXv1DiYY/gYLj44tYEOUBIQAgPz3yYbSXbeH/f+3Yv23xoM+W15e1LhUEJx4ABbbmYjviql6OwEF56CW69VeUxnDFmjOrzyHdhh0YTfbz/vqqkmjDBu/f1su1ILyywPrsY2m8o02KmGSocTktxNez1cjgpxdXQhMPWemRF7goC/QKZnzzf/kV6E+RLlqjw1C9/6fy8OXPUWDYlvjePv5kxUWNYvGExTS2d3XhX5K4gJCCEeUnz2r+xdSukpbVzAm5HXJyqqjEwtGiXP/xB5Vceeqjrc7XKKoM7yHsdlZXqAeWqq5w+TBmCFqry0s+dKRxnAZnJmWwt3srxGmM6QwsrnDT/acTHq87UM2fU9xUValXQhXCMCB+BQLTmOaSUrLSsZFbirPY7DdqiRzgOHIA334Q77nAaKgOUcACsX996KMAvgMdnP46l3MK/vv5Xu9O1vTfmjpjbfh+JmhrYvdtxmArUiqO+XsW6vcXu3fCf/ygBjelozmAH7f/M4DxHr+PDD1X4z9v5DVDCUVfnteZTPX0c53tjIiaOyUzJRCJZm+e53ew06pvqKa0u1Scc0Ba/16y5uxCOkIAQYgfEtgrHvhP7yD+Z7zhMBSosFhzsXDieeEI1+v3mN87nDTBunNqLIqv95liLRi3iguEX8MdP/8iZxjOtx3ce3cmhqkPtu8UBcnJU/Lor4QDv5jkWL1ahunvv1Xd+377q39gUDs/y3nuqOCItzfv39nJJrp4VxwtCiG1CiJ8JIRwEkk2MZMLQCQzrP4zVeZ4PVx2sOohEOu7h0OhYkquVyzqK9dtgW1mlGQa27mthD39/ZVToSDgOH4bXXlPxfD0Gcn5+cMklKoxgs5QXQrD0kqWUVpfyXHabl+fy3OX4Cb9W25dWtMS4sw8Gb/dybN4Mq1cr0QgP13+dFzyrehU1NfDRR8pixFEY00i8bDuiJzk+Hfgeyj5khxDiTSHEXMNnZtKKEILM5EzW5a+jobnBo2N3WYqr0XFfjtxc5baqCYoT2gmHZSVThk0hdkCs84tGj3bcy/HUU0oAfve7Lu/dypw56peqw4fl9PjpLExeyJLNS1orv1bkrmB63HSi+nbwONu6VRknOtuUx5vCISX8/vfqQ+Ouu1y7dswY9e/rRUfVc5qPPlLNsN6uptLwsu2ILmmUUuYBDwD3ouw/nhVC5AohemF3j2/ITMmkuqGazw9+7tFxu2z+0xg2TPUp2K44kpPbG+g5ICkiibIzZeSeyCW7JNt5mEojNVVVCtV2sFs5dgxeeQW+//02MdODluewqa7SWHLJEk7Vn2Lp5qUcOHmAPcf3tPemAueNf7ZERakwmzeEY80a+OILlRDv42K/xJgxKhdTWGjM3Hob77+vHigu1uPGZAA9TTiEEOOEEMtQtiGzgcullKnWP3dpPWLiGS4ZcQkhASEer64qrCwk0C+QYf2HOT/R3189TdsKh87O2OQIZbm+bIv6cdEtHFJ2NuNbtgwaGtSTtivEx7farHfk/CHnc9O4m3hu23M8v+15gM7CceiQWrF0JRxCeKckt6VF2cePHKlCdq7iJc+qXkF9vQoXXnGFrgcpQxgwQDXp9hThQG3k9BUwXkp5p5TyK1DOtqhViIkX6BPYh9mJs/lg/wce7SIvqiwibmAc/n7+XZ+sleRq+4LrFA6tJPf13a+TGJbIeYPP6/oie5VVJ0/C3/6mOqO76lWwh43NekcemfUILbKFv2b/lQlDJ3QO3W3Zor52JRzgnSbAt96CPXuUlXygLgeg9piVVZ5j/XplaeOLaioNIbzay6FHOJZLKd+QUrbGDLRtZaWUbxg2M5NOZCZncqDiAPvLPWeJXVhZ2HViXEPb0OnAAdVAplM4RoSrXfLqmuq4YtQV7U0NHaGZ8dkKx3PPqX2+779f33w7YmOz3pGEsAR+OkUZInSqpgIVpgoNhfN1FBkaveJoaIAHH1RNZt/9bvfGGDhQle6awuE+772nnvgvucS38/Ci7Yge4bjZzrEfeHgeJjpYmLIQ8GwXua4eDo34eCgtVX0DoFs4+gb1bQ2FdQoBOSIkRHXgasJRXQ3PPKPCAXo+vO3RwWa9Iw9e/CA3jbuJH038Uec3t25VTrt6nu7j4tS/U1PnxkKP8Pe/q9zEkiXuVfCYlVXu09SkbEYyM1Vuy5d40XbE4U+dEOIGIcQHQKIQYpXNayNw0tF1JsYRNzCOcUPGeaws93TDacrOlOkXjoQElXfQ8gSj9O+1kBSRpPZSj7tQ/wRtPatefFE1HS5erP/6jnSwWe/IoD6DeOPKNzpXfNXXw9dfq/2j9TB8uMpBGPVL/L//qdXGZZe5N45WWdXS4pl59UY++0y5EfuqmsoWL644nGVyvkTt/BcJPG1zvBrQs+e4iQFkJmfyxBdPUFlXSVhImFtjHaxUie4uS3E1tCqmjz5SYY7+/XXf67HZj1FVV0WAnwvJw9Gj4eOPVXjq6afVHhMdN09ylTlzVPPgqVOd99NwxNdfq/CQnvwGtG8CHD7c+bndYd8+9YTrrq3FmDEqdHf4sGsVaiZtvP++CmHOm9f1uUYTHa1sT+rq2u9LYwDO9hw/KKXcJKXMkFJ+avP6yrpJk4kPyEzJpFk2sy5/ndtjFVZaS3FdyXGA8mJyca+Bi+Iuag216SY1VX1gL16s7E7cWW1odLBZ14Wexj9bjOzlOHlS/VvoaLzsErOyyj1aWpRwzJtn39Lf23ixCdBZqGqz9Wu1EOKUzataCHHK8JmZ2GVazDQi+0R6JFylu4dDY/jwtqdcb2xSo304PvccTJ/umRp5zWbdQZ7DLlu3qrzFsC5KljWMFA4tdOcJ4dDGMIWje2Rnq9CQL6upbPGi7YizFcdF1q/9pZQDbF79pZQ61/gmnsbfz58FyQtYk7eG5hb3un6LKosIDQhlcN/B+i4ICmr78PSmcEjpmdUGtNmsO8hz2EVP458tAweqMJgRwuGC1UuXDBqkdlzsrcIhpapKi4yEG26A119Xqzm9vPeeKpZY6OJK2ii82ATobMUR4exl+MxMHJKZnMnJ2pNsLd7q1jiFlYUkhCXoK4/V0MJVnvjg6oqwMJVLmTJF5Tc8hR2bdYccOaJ6V1wRDjCuJHffPhW/9lROojdXVr3yiio0OP982LABbrlFCemUKarc+csvHVuySKnCVHPmqJ/TnkBPCFUBO1BbxO6w8zK3jvUhl468lAC/ALfLcl3q4dDQvKm8tZ/yypXw7rue3d/Ajs26Q7Sej+4IhxFNgPv2qWo2fx0Nm3rQhMPb+4f4GotFbTc8Z476OThyRLkf/+lPalX6+ONw4YXKQub66+Hf/1Z2Nxo7d6qS6J5QTaURFaXKs30cqkqUUo6wfu34GmH4zEwcMjBkIBfHX+x2nsOlHg6NqVOVeOiN97vL5Mmer/jRbNb15Dm2blXhiIkTXbuHkSsOT4p2aqqqMCst9dyYPZ3GRrjpJpXreu019WHr56d+1h54QPl/lZXB22/DokXKbeAHP1A5BO2cZcvUNVfosM/xFv7+MHiwb4VDQyhuEkI8aP0+TggxzfCZmTglMzmTb45/0+pu6yqVdZVU1VfpL8XV+OUvIS/P+zuceRIHNut22bJFiYar5Y1xcerDp66u+/PsSG2t6tz3ZJiwN1ZWPfywWl28/LLjja8iIuC665SwlJbCjh3K3iU0VDVevvGGypVFRdm/3ld4yXZE134cQAZwo/X7auBvegYXQswTQliEEPlCiE6udEKIeCHEeiHEbiHEJiFErM17Twoh9goh9gkhnhXWQLwQYrIQYo91zNbjvQ1tr4gP93/YretdrqjSEMJ3Rm6eZM4c9WTm7AOzqQm2b3c9TAVtlVXFxd2bnz0sFiV0RghHb9lGdvNm9cH/gx/ANdfou8bPDyZNUgUamzfDiRMqv/GPfxg61W7hpSZAPcKRJqW8E6gDkFJWAEFdXSSE8EcJzHxgDHCDEGJMh9OeAl6XUo4DHgGWWK+9ALgQGAecB0xF2bkDvAj8GEi2vnpA5433SR6UTMqglG6Fq74t+5Zntz0LuNDDca4x17qljLNw1Z496ilfb8e4LUbsBOjJiiqNwYPV03VvWHFUVSk7/oQEePbZ7o8THq42bBo50mNT8xg9SDgarSIgAYQQUYAej4JpQL6UskBK2QC8DXQMCI4BNlj/vNHmfQmEoAQqGAgEjgkhooEBUsqtUlnEvg7oND8698hMzmRD4QZON5zu8tzS6lL+suUvTHp5EmNfGMvru17n8pTLGRs11gsz7YHExSmbdWfCoTX+ubPi8GSeY98+9fTbHWdgRwjReyqr7rpLCfkbb7jkenBWMXSoSuIbvEGXHuF4FlgODBZCPAZsBh7XcV0MYPtbU2w9ZssuQCtLuBLoL4QYJKXcghKSI9bXOinlPuv1tmt/e2MCIIS4XQiRI4TIKSsr0zHds4/MlEwamhtYX2C/Oqi6vpp//397Zx4fVXn9//chbLKqAS1IYiJVSCgYIEVQcCtqEIpspWxVqKhfXKnIT3yxI7S2UOSLYltUiAsuFASjZSv+ImgVJUDCliCoqAFkq4rsAZ/vH8+dMAwzIZPMzJ1Jzvv1mlfuPHPvc8+9M5nPPOc55zy5L3HLK7eQ8HQCI1aMoGqVqvxvxv+y+9HdZPXPokZVlwuzuUkJZdYBKxyXXlq2yfkmjtc11MJxxRWhL6aXmgpbtlTsyKr5822expgxNgm0otKokc1oP3AgrKc5r7PaGDNPRNYBvwIE6OF8iYeCx4BnRWQwsBrYBZwWkZ8DKYBnzuPfItIJOOa3F/92zwZmA6Snp1fI/4iOiR2pV6Me7372Lnc0t4O1otNFLP98OfM2zePtgrc5duoYyRcmM7rTaAa2HEizBqUvTFjh6dzZFk/85BPo2PHc1z2Jf2WZRrvgAjtxGmrhCEcYdGqqLWWyf791XVU0CgvhvvugXTsrHBUZ7yTASy8N22kCCodPkt8+4HXv14wx56uQuwu7TrmHJk5bMc5iUL2cPusAvY0x34vIPcAaY8xh57Wl2An6VzgjJn77rExUi6tGxs8zeHf7u3z8zcfM2zSPN7e8yYGjB4i/IJ4haUMY1GoQ7Zu0Dy7Jr7LgXWbdVzgOHrSrDw4ZUvb+QxmSe+qUtadLl9D05413ZFVFE46ffrKJfUVF8OqrZVv0KpbwFo60tLCdprQJgPuBz4Dtzva6UvS9FrhSRJJFpDrQD8jy3kFEGoiIx4YngDnO9tfADSJSVUSqYSfG840xe4BDItLeiaa6E3i7FLZUWLpd2Y1vD3/LtXOu5cUNL3Jz8s1k9cti94jdzOo6iw4JHVQ0AlFSmfVPP7V/yzK/4SGUSYBffmkLPoYjY78ih+Q+/bTNCp8xw85pVXQiVK8q4IjDGJMMICLPY1cBXOI870IpJqSNMadE5EFgORAHzDHGbBGRSUCOMSYLuBH4k4gYrKvqAefwBdg1zTdhJ8qXGWPecV67H8gELgCWOo9KS8+Unjyw6wHSG6fTK6UX9WpoGbGgCFRmfc0aOxpJTy973wkJkJ1dfhshtMUNfWnc2F57RROOvDy7WmSPHmVblz0WiVDZkdIE5Lc3xtzjeWKMWSoifylN547YLPFpG+e1vQArEr7HnQbuC9BnDjZEVwHqVK/Ds7c/67YZsUvnzra8xKpV8Otfn2lfs8bWMKpTp+x9JyZaQQpm7Y9AhCMU14OI7bciCcexYzBwoA01fv752E5YDYaaNW3trDCPOEoTVbVbRMaISJLzGA1UovoESoXGX5n1n36yE+blcVNBaENy8/Ptr8n69cvflz8qWkjuqFE2Uiwz01a/rUxEIJejNMLRH2iIDcldBFzitClK7OOvzHpBgU0WC5VwhGKeIz8/vBWJU1Nt/P9/K8Cq0MuX2wS/hx8u//K6sUgE1h4/r3AYY/5rjHnEGNPaeTxSiogqRYkdfMusexL/ypIx7k2oRhzGhC8U10NFKT1y4IAtJ9KiBTz1lNvWuEME6lWVtB7HDOfvOyKS5fsIq1WKEkk85Uc8ZdbXrLERV+WNwmnc2E6wl1c49uyx8yThHnFAbLurjIF77rGjpnnzrAuyMuJxVYUxobOkyfFXnL/TwnZ2RYkGWrY8U2b9zjutcFxzjf3SLw9Vq1rxKK9whDOiykNiItSqFdvCMWcOLF4MU6fC1Ve7bY17/OxnNjjg0KGwzYmVJBz7AYwxq8JyZkWJFrzLrB86BJs3h24d6VAkAUZCOKpUie3Iqh07bMn/m26CRx912xp38U4CDJNwlPSTarFnQ0QWhuXsihIteMqsv/SSHeKXd2LcQyiSAAsKbDiv5wshXMRqZJUxcP/9doT30kvlHynGOhHI5SjpDnsHPuuKf0rFxjPP8RcnRaldiNYqS0y0tZLK42/2RFSFOxchNdXaeuhQeM8Tav71LxsVN2HCmYCEyoz3iCNMlCQcJsC2olQ8PGXWCwtt9NJFF4Wm34QEuwpgeaqVhjuiyoPHFRZLkVUnT8KIEXYd9gceOP/+lYEIlB0pSTiuFpFDIvIj0MrZPiQiP4pIjP0kUZRS0Lmz/RsqNxWUPyT3hx/sF0A45zc8xGJk1axZtvjjX/9a8QsYlpYLL7T5SW4IhzEmzhhTzxhT1xhT1dn2PNeCSErFI5zCUdZ5jkhMjHtITrZfOLEiHAcO2PXDb7sNbr/dbWuiB5Gw53JU8lkkRfGia1frJ//tb0PXZ3lHHJEUjqpVrcsnVlxV48bB4cMwfXrlqUVVWsJcdkSFQ1E81KgB48fboX6oaNjQ9ltW4SgogOrV7WggEsRKZNXmzfCPf8CwYWdcbMoZwlx2RIVDUcJJlSp2GdnyjDiuusqOBiJBairs3AlHjkTmfGXBGPjDH2yOwoQJblsTnairSlFinPIkAYa7uKEvqan2i3nbtsidM1jefdcma06YAPHxblsTnTRqZEuvnDgRlu5VOBQl3JQ1CfD4cfjii8iE4nqI9sgqT/ht8+bWTaX4J8xJgBEa/ypKJSYxEXbvhtOnIS6u9Mdt327XBonkiOPnP7dusWgVjmeftfdlyRINvy0J71yOyy8Pefc64lCUcJOQYEUj2MnKSEZUeahWzc6pRKNw7N8PkyZBRgZ06eK2NdFNmEccKhyKEm7KGpKbn2/DTJs1C71NJRGtxQ69w2+Vkglz2REVDkUJN2VNAiwogKSkyK8rkZoKn39u51iihY0bYfZsW8wwkiOwWOWSS2xEnwqHosQoiYn2b1lGHG58Saam2rmVzz6L/Ln9YYwtla7ht6UnLs7mEKlwKEqMUr8+1K0bnHCcPm1DYiMZUeUh2paRzcqyqzNOnAgXX+y2NbFDGHM5wiocIpIhIttEZIeIjPLz+uUi8p6IbBSR90WkidN+k4jkej2Oi0gP57VMEfnS67W0cF6DooSEYHM5vvrKuorcGHFcdZV1c0TDPMeJEzb8NiUF/ud/3LYmtghj2ZGwheOKSBwwC7gFKATWikiWMcb70zgNeNkY85KI3Az8CfidMSYbSHP6uRjYAazwOm6kMWZBuGxXlJATrHC4EVHloWZNaNo0OoTjmWfsfMvSpRp+GyyNGkFeXli6DueIox2wwxjzhTHmJPAGcIfPPqnA/3e2s/28DtAHWGqMORo2SxUl3ASbBOimcEB01Kzatw+efNJWvs3IcNeWWORnP4O9e63bM8SEUzguA7x/YhU6bd7kAb2c7Z5AXRHxrSHQD3jdp22K4956WkRq+Du5iNwrIjkikrN///6yXYGihIrERPtFWNoSEAUFNjLGLZ9+aqqdHC8qcuf8AGPHwtGjdq0NJXgaNbKicfBgyLt2e3L8MeAGEdkA3ADsAorlUUQaAS2B5V7HPAE0B34JXAw87q9jY8xsY0y6MSa9YcOGYTJfUUqJJyS3sLB0+7sVUeUhNRVOnYIdO9w5f14evPCCXdXPjQCBikCjRjYPaN++kHcdTuHYBXgvANzEaSvGGLPbGNPLGNMaGO20fe+1S19gkTGmyOuYPcZyApiLdYkpSnQTTBKgMdEhHOCOu8pT/fbCC23Sn1I2une3I9xf/CLkXYdTONYCV4pIsohUx7qcsrx3EJEGIuKx4Qlgjk8f/fFxUzmjEEREgB7A5jDYriihJZgkwH374Lvv3P2l7clWd0M43n4bsrNteRENvy071aqFLaAgbMJhjDkFPIh1M+UD840xW0Rkkoh0d3a7EdgmIp8BlwJTPMeLSBJ2xLLKp+t5IrIJ2AQ0ACaH6xoUJWQEM+Jwe2IcoHZtm7UeaeE4cQIee8yOeO67L7LnVkpNWKvjGmOWAEt82sZ5bS8A/IbVGmN2cu5kOsaYm0NrpaJEgAsugAYNYkc4wJ3IqpdftuG3y5ZFbvEqJWjcnhxXlMpDaXM58vOhTh27cqCbpKba7PUwhHMGZM4caNECbr01cudUgkaFQ1EiRWmFo6DAzm+IhN+mkkhNta6jL7+MzPkKCmDNGhg82P1rV0pEhUNRIkVpkwDdjqjyEOnIqpdessX5Bg2KzPmUMqPCoSiRIjERfvgBfvwx8D4//mhzPaIhd8EjXpEQjtOn7fxGly5nVq9TohYVDkWJFKWJrCoosH+jYcRRr56dZ4mEcPz733Z53cGDw38updyocChKpCiNcERLRJWHSEVWzZ0L8fHw61+H/1xKuVHhUJRIUZokwPx8G4batGlkbDofqanWpp9+Ct85vvsOFi+GAQOgevXwnUcJGSocihIpGje261ycb8Rx5ZXRU0I8JcUWGgx22dtgeP11OHkShgwJ3zmUkKLCoSiRolo1W3jufHMc0eKmgshEVmVmQqtWkKZrssUKKhyKEklKyuU4edJWo41G4fjkk/D0v2ULrF1rRxuauxEzVNqc/qKiIgoLCzl+/LjbpihhpGbNmjRp0oRq0eL6SUgIvCrbjh02LDUaQnE9XHyxXUTp2Wfh0Uft+umhJDPTzukMHBjafpWwUmmFo7CwkLp165KUlIToL50KiTGGgwcPUlhYSHJystvmWBIT4Z13bOlw389dtEVUeZg8GdLTYfp0mDgxdP0WFcErr0C3bqBr5sQUldZVdfz4ceLj41U0KjAiQnx8fHSNKhMS4Phx/6uyeYQjmkYcAG3bQu/eVjhCuZrm8uV2aVPN3Yg5Kq1wACoalYCoe49LyuXIz7cjktq1I2tTaZg0yUZX/fnPoetz7lw70rj99tD1qUSESi0cihJxShKOaIuo8iY1FX73OzvXsWvX+fc/HwcOWJfdoEHRE3qslBoVjihjxowZHD16NGT9JSUlceDAgTIf//7779OtWze/r/Xv359WrVrx9NNPM27cOFauXBnU8f72ycrK4qmnngJg8eLFbHVjBbpwEigJ8Kefols4AMaPt3ZODsHaaa+/buc41E0Vk1TayfFoZcaMGQwaNIhatWq5cv7Tp08TFxd33v2+/fZb1q5dy44dO0J6/u7du9O9u10gcvHixXTr1o1UT0hoReCSS2x2tO+I45tvrCso2uY3vElOhnvugdmzYeRIuOKKsvc1dy60aWPzN5SYQ4UDYPhwyM0NbZ9paTBjRsCXjxw5Qt++fSksLOT06dOMHTuWvXv3snv3bm666SYaNGhAdnY2w4YNY+3atRw7dow+ffow0YlqSUpK4q677uKdd96hqKiIf/7znzRv3pyDBw/Sv39/du3aRYcOHTDGFJ+zR48efPPNNxw/fpxHHnmEe++9F4A6depw3333sXLlSmbNmsXhw4cZPnw4tWrVomPHjn7tv/XWW9m1axdpaWk888wzvPjii3Tr1o0+ffqwbNkyv8cfOXKEhx56iM2bN1NUVMSECRO44447zuo3MzOTnJwcBgwYQFZWFqtWrWLy5MksXLiQ3/zmN6xfvx6A7du389vf/rb4ecxQpYotHOgrHNEaUeXL6NF2saUJE2w127KQlwcbNsDMmSE1TYkc6qpyiWXLltG4cWPy8vLYvHkzGRkZPPzwwzRu3Jjs7Gyys7MBmDJlCjk5OWzcuJFVq1axcePG4j4aNGjA+vXrGTZsGNOmTQNg4sSJdOzYkS1bttCzZ0++9nKJzJkzh3Xr1pGTk8PMmTM56ET2HDlyhGuuuYa8vDzS09O55557eOedd1i3bh3ffvutX/uzsrJo2rQpubm5dOrUqbj9+PHjAY+fMmUKN998M59++inZ2dmMHDmSI0eO+O3/2muvpXv37kydOpXc3FyaNm1K/fr1yXUEfu7cuQyJ1RIV/pIAY0U4GjeGhx6CV1+1yXtlITPTzmsMGBBS05TIoSMOKHFkEC5atmzJiBEjePzxx+nWrdtZX77ezJ8/n9mzZ3Pq1Cn27NnD1q1baeUM73v16gVA27ZteeuttwBYvXp18XbXrl256KKLivuaOXMmixYtAuCbb75h+/btxMfHExcXR+/evQEoKCggOTmZK6+8EoBBgwYxe/bsUl9XScevWLGCrKysYpE7fvz4WcJ2PoYOHcrcuXOZPn06b775Jp9++mmpj40qEhLggw/ObsvPt9VhYyGf4fHH4e9/h3HjYOHC4I49edKKTvfu9nqVmERHHC5x1VVXsX79elq2bMmYMWOYNGnSOft8+eWXTJs2jffee4+NGzfStWvXs3ISatSoAUBcXBynTp0q8Xzvv/8+K1eu5OOPPyYvL4/WrVsX91WzZs1SzWuUF2MMCxcuJDc3l9zcXL7++mtSgviF3bt3b5YuXcq7775L27ZtiY/VL57ERLtYk/da3tGy6l9piI+HESPgrbdsuZBgWLLERlTF6mhRAcIsHCKSISLbRGSHiIzy8/rlIvKeiGwUkfdFpInTfpOI5Ho9jotID+e1ZBH5xOnzTRGJyTrMu3fvplatWgwaNIjeoOS0AAAPjElEQVSRI0cW++rr1q3Lj84KcYcOHaJ27drUr1+fvXv3snTp0vP2e/311/Paa68BsHTpUr777jsAfvjhBy666CJq1apFQUEBa9as8Xt88+bN2blzJ59//jkAr7/+elDXVdLxt912G88880zxvMuGDRtK7Mv7XoAVuNtuu41hw4bFrpsK7Ijj9GnwdgNGe0SVL3/4gxWQMWOCOy4z067wd9ttYTFLiQxhEw4RiQNmAV2AVKC/iPiGx0wDXjbGtAImAX8CMMZkG2PSjDFpwM3AUWCFc8yfgaeNMT8HvgPuDtc1hJNNmzbRrl070tLSmDhxImOcf8B7772XjIwMbrrpJq6++mpat25N8+bNGTBgANddd915+x0/fjyrV6+mRYsWvPXWWyQmJgKQkZHBqVOnSElJYdSoUbRv397v8TVr1mT27Nl07dqVNm3acMkllwR1XSUdP3bsWIqKimjVqhUtWrRg7NixJfbVr18/pk6dSuvWrYuFaODAgVSpUoVbb701KLuiCt9cjgMH7COWhKNePRg1ClasgFWrSnfMvn3wr3/ZfJCq6iWPaYwxYXkAHYDlXs+fAJ7w2WcLkOBsC3DITz/3AvO89jkAVPV3jkCPtm3bGl+2bt16TpsS/UydOtWMGTMmqGOi7r3OyzMGjJk/3z5fvdo+X7LEXbuC5ehRYxo3Nua664z56afz7z99ur3OLVvCb5sSEoAc4+c7NZyuqssA79CRQqfNmzygl7PdE6grIr6O636Ax98RD3xvjPE49P31qVRQevbsycsvv8wjjzzitinlwxkFFicBxkpElS8XXGBdVf/5DyxbVvK+xtjcjXbtzpRqV2IWtyfHHwNuEJENwA3ALqB4xlBEGgEtgeXBdiwi94pIjojk7A9lYTbFNRYtWsTGjRtp0KCB26aUj/r1oU6dM66q/HyoVeuMoMQSd99tEwNHjy55edkNG2DTJs0UryCEUzh2AQlez5s4bcUYY3YbY3oZY1oDo52277126QssMsYUOc8PAheKiMdBek6fXn3PNsakG2PSG8ZCiKNSeRA5O5cjPx+aNbPJgbFG9eq21PqGDTbKKhCZmVCjBvTrFzHTlPARzk/qWuBKJwqqOtbllOW9g4g0EBGPDU8Ac3z66M8ZNxWOzy0b6OM03QW8HQbbFSW8eAtHrEVU+TJggHU/jR17doixhxMnYN486NEDvPKKlNglbMLhzEM8iHUz5QPzjTFbRGSSiHR3drsR2CYinwGXAlM8x4tIEnbE4huy8TjwqIjswM55vBiua1CUsOERjiNH4KuvYls44uLgySetAL766rmvv/su/Pe/mrtRgQhrTJwxZgmwxKdtnNf2AmBBgGN34mfi2xjzBdAupIYqSqRJTLR5HJs22efRXNywNPTsaRd8Gj/euqOc5FTAToo3bgydO7tnnxJSYtCpWnnxLbl+++238/3335dwxNl4lywvDzfeeCM5OTnntH/wwQe0aNGCtLQ0du3aRZ8+ffwcHfj4QPt4rvP777/nueeeK7f9UYEnl8NTij6WRxxg522mTLGjpxdeONO+Z4+NuLrzTjsyUSoEKhwxhK9wLFmyhAsvvLDUx3fv3p1Ro85J4A8Z8+bN44knniA3N5fLLruMBQv8DiaDxnOdFVI4VqywX6hOba+Y5tZboVMnu16H53P66qt23kOjqSoUmr4JDF82nNxvQ1tWPe1naczIKLl44vTp05kzx8YDDB06lOHDh7Nz504yMjJo27Yt69evp0WLFrz88su88MIL55RcT0pKIicnh8OHD5ORkUH79u356KOP+OUvf8mQIUMYP348+/btY968ebRr1664ZPmzzz5LWlpasR3btm1j2bJlpKen+y17fuzYMYYMGUJeXh7Nmzfn2LFj51zLCy+8wPz581m+fDlLly5lypQpdOvWjc2bN5d4/IoVKxg/fjwnTpygadOmzJ07lzp16pzVt+c6R40axeeff05aWhq33HILe/fupVevXvTo0QOwWeV9+/Y9p1R7VOIRjo8/hqZNbXRSrOMZdVx/vV0pcORIG03VoYONGlMqDDricIl169Yxd+5cPvnkE9asWcPzzz9fXLtp27Zt3H///eTn51OvXj2ee+45vyXXvdmxYwcjRoygoKCAgoICXnvtNT788EOmTZvGH//4x3P29xQafPLJJ0lPT+faa68NWPb8b3/7G7Vq1SI/P5+JEyeybt26c/obOnRocRn0efPmnfVaoOMPHDjA5MmTWblyJevXryc9PZ3p06cHvGdPPfVUcSn3qVOncvfdd5OZmQnYWlwfffQRXbt2LfV74Coe4Th1KvbdVN506gRduti1yd97D7Zu1UnxCoiOOOC8I4Nw8OGHH9KzZ09q164N2BLpH3zwAd27dychIaG4LtWgQYOYOXMmjz32WIn9JScn07JlSwBatGjBr371K0SEli1bsnPnTr/HbN++nZEjR5KdnU21atUClj1fvXo1Dz/8MACtWrUqLuteWgIdv2bNGrZu3Vp8rSdPnqRDhw6l7veGG27g/vvvZ//+/SxcuJDevXtTNVZqINWqZYsEHjxYsYQDrKuqbVvo29dml/ft67ZFSoiJkf+yyoWIlPjcHzW8oliqVKlS/LxKlSp+S64fPnyYvn378vzzz9OoUSPgTNnzZhFyKxhjuOWWW4KuwOvNnXfeyauvvsobb7zB3LlzQ2hdBEhIqJjC0aYN9OkDCxbAwIE2U16pUKiryiU6derE4sWLOXr0KEeOHGHRokXFizl9/fXXfPzxxwC89tprxcuv+pYZLw+///3vGTJkyFkLSAUqe+5dqn3z5s1nrUJYGgId3759e/7zn/8Ur1t+5MgRPvvss4D9+Lv+wYMHM8NZiCvm1ib3uKtiPRTXH5Mn2+t78EG3LVHCgAqHS7Rp04bBgwfTrl07rrnmGoYOHUrr1q0BaNasGbNmzSIlJYXvvvuOYcOGAWeXXC8PX331FQsWLGDOnDmkpaWRlpZGTk5OwLLnw4YN4/Dhw6SkpDBu3Djatm0b1PkCHd+wYUMyMzPp378/rVq1okOHDhQUFATsJz4+nuuuu45f/OIXjBw5EoBLL72UlJSU2Fyfw1ObqiIKR7NmtohjgPL9Smwjnl+XFZn09HTjmzeQn58f1OpzkWLnzp3F0UjK+Tl69CgtW7Zk/fr11A/gEonW95oNG+wE8nnmrxTFLURknTEm3bddRxxKzLJy5UpSUlJ46KGHAopGVNO6tYqGEpPo5HiUkZSUpKONUtK5c2e++uort81QlEpHpR5xVAY3XWVH32NFCT2VVjhq1qzJwYMH9YulAmOM4eDBg9SsWdNtUxSlQlFpXVVNmjShsLAQXR2wYlOzZk2aNGnithmKUqGotMJRrVo1kpOT3TZDURQl5qi0ripFURSlbKhwKIqiKEGhwqEoiqIERaXIHBeR/UBZA/4bAAdCaE64iBU7IXZsVTtDS6zYCbFja7jtvNwY09C3sVIIR3kQkRx/KffRRqzYCbFjq9oZWmLFTogdW92yU11ViqIoSlCocCiKoihBocJxfma7bUApiRU7IXZsVTtDS6zYCbFjqyt26hyHoiiKEhQ64lAURVGCQoVDURRFCQoVjhIQkQwR2SYiO0RklNv2eBCRBBHJFpGtIrJFRB5x2ieIyC4RyXUet0eBrTtFZJNjT47TdrGI/FtEtjt/L3LZxmZe9yxXRA6JyPBouZ8iMkdE9onIZq82v/dQLDOdz+xGEWnjsp1TRaTAsWWRiFzotCeJyDGve/t3l+0M+F6LyBPO/dwmIrdFys4SbH3Ty86dIpLrtEfunhpj9OHnAcQBnwNXANWBPCDVbbsc2xoBbZztusBnQCowAXjMbft8bN0JNPBp+wswytkeBfzZbTt93vdvgcuj5X4C1wNtgM3nu4fA7cBSQID2wCcu23krUNXZ/rOXnUne+0XB/fT7Xjv/V3lADSDZ+U6Ic9NWn9f/CoyL9D3VEUdg2gE7jDFfGGNOAm8Ad7hsEwDGmD3GmPXO9o9APnCZu1YFxR3AS872S0APF23x5VfA58aYqFla0BizGvivT3Oge3gH8LKxrAEuFJFGbtlpjFlhjDnlPF0DuF7jPsD9DMQdwBvGmBPGmC+BHdjvhohQkq0iIkBf4PVI2eNBhSMwlwHfeD0vJAq/nEUkCWgNfOI0Pei4Bea47QJyMMAKEVknIvc6bZcaY/Y4298Cl7pjml/6cfY/YrTdTw+B7mE0f25/jx0NeUgWkQ0iskpEOrlllBf+3utovp+dgL3GmO1ebRG5pyocMYyI1AEWAsONMYeAvwFNgTRgD3YY6zYdjTFtgC7AAyJyvfeLxo6xoyImXESqA92BfzpN0Xg/zyGa7mEgRGQ0cAqY5zTtARKNMa2BR4HXRKSeW/YRI++1D/05+0dOxO6pCkdgdgEJXs+bOG1RgYhUw4rGPGPMWwDGmL3GmNPGmJ+A54ngkDoQxphdzt99wCKsTXs97hPn7z73LDyLLsB6Y8xeiM776UWgexh1n1sRGQx0AwY6Iofj+jnobK/Dzh1c5ZaNJbzXUXc/AUSkKtALeNPTFsl7qsIRmLXAlSKS7PwS7QdkuWwTUOzbfBHIN8ZM92r39mX3BDb7HhtJRKS2iNT1bGMnSjdj7+Ndzm53AW+7Y+E5nPULLtrupw+B7mEWcKcTXdUe+MHLpRVxRCQD+H9Ad2PMUa/2hiIS52xfAVwJfOGOlSW+11lAPxGpISLJWDs/jbR9fugMFBhjCj0NEb2nkYoOiMUHNkLlM6xyj3bbHi+7OmJdExuBXOdxO/AKsMlpzwIauWznFdiIlDxgi+ceAvHAe8B2YCVwcRTc09rAQaC+V1tU3E+smO0BirA+9rsD3UNsNNUs5zO7CUh32c4d2DkCz+f0786+vZ3PRC6wHvi1y3YGfK+B0c793AZ0cfu9d9ozgf/x2Tdi91RLjiiKoihBoa4qRVEUJShUOBRFUZSgUOFQFEVRgkKFQ1EURQkKFQ5FURQlKFQ4FMUlRKSviOQ7j9+4bY+ilBYNx1UUF3ASI98HbsDmXqwCrjfGHHbTLkUpDTriUBR3yACWGmMOG1vheCW25ImiRD1V3TZAUSopicAQEenmPG8I7HXRHkUpNTriUBT3+JsxJs0Ykwb8w21jFKW0qHAoijt8zdnrOiQAUbN4lKKUhE6OK4oLOGup5ADXYSfHPwZa6+S4EgvoHIeiuIAx5rCIjMEupwrwhIqGEivoiENRFEUJCp3jUBRFUYJChUNRFEUJChUORVEUJShUOBRFUZSgUOFQFEVRgkKFQ1EURQkKFQ5FURQlKP4Pjtgy4zDWsJAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(thetas * 180 / np.pi, standard_fidelities, 'r')\n",
    "plt.plot(thetas * 180 / np.pi, optimized_fidelities, 'g')\n",
    "plt.xlabel('θ')\n",
    "plt.ylabel('Fidelity of Rx(θ)')\n",
    "plt.legend(['standard fidelity', 'optimized fidelity'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.9842197286344213, 0.9908958547858515)"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(standard_fidelities), np.mean(optimized_fidelities)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### plot standard vs. optimized rotation paths"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit.quantum_info.synthesis.one_qubit_decompose import OneQubitEulerDecomposer\n",
    "from qiskit.quantum_info import Operator\n",
    "from qiskit.extensions.standard import U3Gate, RXGate\n",
    "\n",
    "standard_paths = []\n",
    "optimized_paths = []\n",
    "for theta in thetas:\n",
    "    u3 = U3Gate(theta, -np.pi/2, np.pi/2)              # RX(θ) in terms of U3\n",
    "\n",
    "    decomposer = OneQubitEulerDecomposer(basis='U1X')  # standard way\n",
    "    standard_circ = decomposer(Operator(u3))\n",
    "    standard_paths.append(standard_circ)\n",
    "\n",
    "    decomposer = OneQubitEulerDecomposer(basis='ZXZ')  # optimized way\n",
    "    optimized_circ = decomposer(Operator(u3))\n",
    "    optimized_paths.append(optimized_circ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit.quantum_info.operators.pauli import pauli_group, Pauli\n",
    "from qiskit.quantum_info.states import DensityMatrix\n",
    "\n",
    "paulis = pauli_group(1)\n",
    "X = paulis[1]\n",
    "Y = paulis[2]\n",
    "Z = paulis[3]\n",
    "\n",
    "ideal_rhos = [DensityMatrix(sv) for sv in ideal_statevectors]\n",
    "ideal_Xexps = [np.real(np.trace(np.dot(X.to_matrix(), rho.data))) for rho in ideal_rhos]\n",
    "ideal_Yexps = [np.real(np.trace(np.dot(Y.to_matrix(), rho.data))) for rho in ideal_rhos]\n",
    "ideal_Zexps = [np.real(np.trace(np.dot(Z.to_matrix(), rho.data))) for rho in ideal_rhos]\n",
    "\n",
    "std_rhos = [DensityMatrix(sv) for sv in standard_statevectors]\n",
    "std_Xexps = [np.real(np.trace(np.dot(X.to_matrix(), rho.data))) for rho in std_rhos]\n",
    "std_Yexps = [np.real(np.trace(np.dot(Y.to_matrix(), rho.data))) for rho in std_rhos]\n",
    "std_Zexps = [np.real(np.trace(np.dot(Z.to_matrix(), rho.data))) for rho in std_rhos]\n",
    "\n",
    "opt_rhos = [DensityMatrix(sv) for sv in optimized_statevectors]\n",
    "opt_Xexps = [np.real(np.trace(np.dot(X.to_matrix(), rho.data))) for rho in opt_rhos]\n",
    "opt_Yexps = [np.real(np.trace(np.dot(Y.to_matrix(), rho.data))) for rho in opt_rhos]\n",
    "opt_Zexps = [np.real(np.trace(np.dot(Z.to_matrix(), rho.data))) for rho in opt_rhos]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x11b973f28>]"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeXhTZfrw8e/TlJR9axHZQXBBBxVpdaKCjVUEdJS6RHQQlxmrzPjqDCJanaLTKpWK+/xEqqOCC5gByjCK49ImskWmgCiu7LILssnadLnfP04KKZSuadPl/lxXruRsOU/T5Nzn2Y2IoJRSqvGKCHcClFJKhZcGAqWUauQ0ECilVCOngUAppRo5DQRKKdXIRYY7AVURExMjPXv2DHcylFKqXlm2bNkvItLh+PX1MhD07NmTpUuXhjsZSilVrxhjfiptvRYNKaVUI6eBQCmlGjkNBEop1chpIFBKqUZOA4FSSjVyIQkExpg3jDE7jDHfnGS7Mca8ZIxZY4z52hhzQdC2240xqwOP20ORHqWUUhUXqhzBW8CQMrYPBU4PPJKAyQDGmPbA48BFwIXA48aYdiFK0wlWZvrwXpXOykxfTZ1CKaXqnZD0IxCR+caYnmXsch0wTawxr78wxrQ1xnQC4oFPRWQ3gDHmU6yAMj0U6Qq2MtNH73sS6Isf/yd2VpJNvyRHhY/dNctL9A3xFT5GKaXqi9rqUNYF2BS0vDmw7mTrT2CMScLKTdC9e/dKJ2DXLC97u+axsGcRl27Io3CWFypwUa9OAFFKqfqg3vQsFpFMIBMgNja20rPpbP5dNElxRfhtYC8sIvPU6Aodt2uWl774iaQQwc+uCgYQpZSqL2qr1dAWoFvQctfAupOtD7lN5+0ir0kEhRHgbxLBpvN2Vei46BviWdDVxpOXGhZ2tRF9Q3yFz+nzQXq69ayUUnVVbeUI5gL3GWNmYFUM7xORbcaYj4EJQRXEg4HkmkhAfM94oiKj8Bf6sdvsxPeMr9BxB4bCsG0Gv4DdGHKGVux8Ph8kJIDfD3Y7ZGeDQzMSSqk6KCSBwBgzHaviN8YYsxmrJVATABF5FZgHDAPWAIeAOwPbdhtj0oDcwFulFlcch5qjm4PsUdl4N3iJ7xmPo1vFrsreDV7yTQFFCPmmAO8Gb4WO9XrhgjwfA4u8LMiLx+t1aCBQStVJoWo1dEs52wX480m2vQG8EYp0lMfRzVHhAFAsvmc8dpu90jmJa6J9/LZzPIt75pO2oQkx0V5AI4FSqu6pN5XF4VLVnMSBndO4+nZ/oHLaT/bOaWggUErVRRoIKqAqOQlvT/CvwqqcFmu5ou/g81lFS/HxWq+glKp5OtZQDYkfNAp7ZBQ2DPYmUcQPGlWh43w+SI73ceCxdJLjfdriSClV4zRHUEMc3Rxk3+GpdJHS6mk+5vkTsOPH77czc1o2Ds0WKKVqkAaCGlSVIqXL8GIP6sB2GV60bkEpVZO0aKiO6TEqHhNlp9DYiIiy02NUfMUO1N5rSqkq0hxBXeNwYPNkV662WHuvKaWqQXMEdZCvK6Rfaj1XiNdrBYHCQuvZ663B1CmlGhrNEdQxvk0+EqYlHO3Alj0qu/x6hvh4CiPtUOSHSDu2+PhaSatSqmHQHEEd493gxV/op1AK8Rf68W7wlnuMDwcJks140kiQbHxauayUqgTNEdQxVRnSwuuFhYUOPhcHtkJrWasIlFIVpYGgjqnKkBbx8VYdcXFdsZYMKaUqQwNBHVTZ/gcOB7wwy8esZV5uGBCvHdCUUpWigaAB8G3y8ZflCfiL/CxYbqffbypQwayUUgFaWdwAVKWCWSmlimkgaACKK5htxlapORMA7ZGslArZDGVDgBcBG/C6iDx93PbnAWdgsTlwioi0DWwrBFYGtm0UkWtDkabGpKpzJuDzUehMwPj9iN1u9WjW+gWlGp1qBwJjjA34P+BKYDOQa4yZKyLfFe8jIn8N2v//Af2D3uKwiJxf3XQ0dlUZ4O6naV665PmxUUh+np/N07z00ECgVKMTiqKhC4E1IrJORPzADOC6Mva/BZgegvOqavqcePzYycdGPnY+Jz7cSVJKhUEoAkEXYFPQ8ubAuhMYY3oAvYCcoNVNjTFLjTFfGGOGn+wkxpikwH5Ld+7cGYJkq9NHORhmz+YJk8Ywezanj9LcgFKNUW03Hx0BzBSRwqB1PURkizHmNCDHGLNSRNYef6CIZAKZALGxsVI7yW3YHA64dQ7MWga3DtDqAaUaq1AEgi1At6DlroF1pRkB/Dl4hYhsCTyvM8Z4seoPTggEKvS0/4FSCkJTNJQLnG6M6WWMsWNd7Ocev5Mx5iygHeALWtfOGBMVeB0DXAJ8d/yxqmZo/wOlFIQgRyAiBcaY+4CPsZqPviEi3xpjUoGlIlIcFEYAM0QkuFinLzDFGFOEFZSeDm5tpGpWVQa4U0o1PKbkdbl+iI2NlaVLl4Y7GQ2Cb5Ov8v0PwOqAVplZ1JRSYWeMWSYiscev17GGGrmq9D/QqTGValh0iAlVeV4vkmdNjSl5OjWmUvWdBgJVaSuj4zlcZHVEO1xkZ2V0fLiTpJSqBi0aUpX2wS4HH0ZkM7DIy4KIeK7e5aBfuBOllKoyDQSq0uLjIS3KwRd+B3Y7PBMf7hQppapDA4GqNIfDqh/WRkNKNQwaCFSVOBwaAJRqKLSyWFWJb5OP9AXp+DbphDZK1XeaI1CV5tvkI2FawtEeydmjdIwipeozzRGoStMxipRqWDQQqEqr1hzJSqk6R4uGVKVVeY5kpVSdpIFAVUmVxihSStVJWjSkatXKTB/eq9JZmamtjZSqKzRHoGrNykwfve9JoC9+/J/YWUk2/ZI0V6FUuGmOQNWaXbO82PETSSFN8LNrljfcSVJKEaJAYIwZYoz50RizxhjzSCnb7zDG7DTGrAg8/hi07XZjzOrA4/ZQpEfVTdE3xOPHGrU0HzvRN8SHO0lKKUJQNGSMsQH/B1wJbAZyjTFzS5ly8n0Rue+4Y9sDjwOxgADLAsfuqW66VN3TL8nBSrLZNctL9A3xWiykVB0RijqCC4E1IrIOwBgzA7iOik1CfxXwqYjsDhz7KTAEmB6CdKk6qF+SAzQAKFWnhKJoqAuwKWh5c2Dd8W4wxnxtjJlpjOlWyWMxxiQZY5YaY5bu3LkzBMlWSikFtVdZ/B+gp4icC3wKTK3sG4hIpojEikhshw4dQp5ApZRqrEIRCLYA3YKWuwbWHSUiu0QkL7D4OjCgoseqhkVHLVWq7glFHUEucLoxphfWRXwEcGvwDsaYTiKyLbB4LfB94PXHwARjTLvA8mAgOQRpUnWQjlqqVN1U7RyBiBQA92Fd1L8H3CLyrTEm1RhzbWC3+40x3xpjvgLuB+4IHLsbSMMKJrlAanHFsWp4dNRSpeqmkPQsFpF5wLzj1o0Pep3MSe70ReQN4I1QpEPVbcWjlhbnCCo1aqnPp3NjKlVDdIgJVWuqPGqpz0ehMwHj9yN2OzZPtgYDpUJIA4GqVVUZtfSnaV665PmxUUh+np/N07z00ECgVMjoWEOqzvuckkNTfE58uJOkVIOigUDVeaePcjDMns0TJo1h9mxOH6W5AaVCSYuGVJ3ncEC614HX6yA9XqsHlAo1zRGo+qGrDy5Nt56VUiGlOQJV52lHNKVqluYIVJ2nHdGUqlkaCFSdV9wRzWZsle+IppQqlxYNqTqvyh3RlFIVooFA1QtV6YimlKoYLRpSDZvPB+np1rNSqlSaI1ANl88HCQng94PdDtk6RpFSpdEcgWq4vF4kzw+Fhdaz1xvuFClVJ2kgUA3Wyuh4DhdZYxQdLrKzMjo+3ElSqk4KSSAwxgwxxvxojFljjHmklO1jjDHfBSavzzbG9AjaVmiMWRF4zA1FepQC+GCXg8ER2YwnjcER2XywS4uFlCpNtesIjDE24P+AK4HNQK4xZq6IfBe025dArIgcMsaMBjKAmwPbDovI+dVNh1LHi4+HtCgHX/gd2O3wTHy4U6RU3RSKyuILgTUisg7AGDMDuA44GghExBO0/xfAyBCcV6kyORzwwiwfs5Z5uWFAPA6tKFaqVKEIBF2ATUHLm4GLytj/D8BHQctNjTFLgQLgaRGZU9pBxpgkIAmge/fu1Uqwahx8m3z8ZXkC/iI/C5bb6fcbHaNIqdLUamWxMWYkEAs8E7S6h4jEArcCLxhjepd2rIhkikisiMR26NChFlKr6jsdo0ipiglFINgCdAta7hpYV4Ix5grgMeBaEckrXi8iWwLP6wAv0D8EaVJKxyhSqoJCUTSUC5xujOmFFQBGYN3dH2WM6Q9MAYaIyI6g9e2AQyKSZ4yJAS7BqkhWqtp0jCKlKqbagUBECowx9wEfAzbgDRH51hiTCiwVkblYRUEtgX8ZYwA2isi1QF9gijGmCCt38vRxrY2UqhYdo0ip8hkRCXcaKi02NlaWLl0a7mQopVS9YoxZFqiTLUF7FiulVCOngUCp4+mIpaqR0dFHlQqmI5aqRkhzBEoF0xFLVSOkgUCpIDpiqWqMtGhIqSAf7HLwYUQ2A4u8LIiI5+pdDvqFO1FK1TANBEoF0RFLVWOkgUCpIA6HVT/s9VpBQeuJVWOggUCp4zgcGgBU46KVxUop1chpIFDqOL5NPtIXpOPbpB3KVA1YuhRGjYLc3HCn5CgNBEoF8W3ykTAtgRRPCgnTEupUMBARVmxfwUtLXjq6Luk/Sfwt52/sO7IvjClTFbJnD/zpT3DhhfD221anxTpCA4FSQeraZDYiwvJty0n+LJkz/nEG/af0Z8zHY9i6fyuFRYUc8B/gqQVP0evFXmQsyuBQ/qGwprdeOXLEuiufMgVGj4aPPir/mKp67z0480zrXPffD7t2wSWXWNs++QSKimru3BWglcV1RVERLF4MM2fCokVwww1w333QsmW4UxZSeQV52CJsREbUza9e8WQ2/kJ/2CazEREKigpoYmvC9G+m8/vZv8dmbFze63LGXTyO4WcNp0MLa5a+9254j4cufojHch7j4c8e5oUvXuD9G99nYI+BtZ7uEr76Cr7/HrZtsy56djtER8Of/2xt//xz2LkTmjWDpk2t5zZt4JxzrO27d1vHtGgB1tD11XP4MOTlQdu2sGoVuFzw7bdQUGBtb97culsH+OkniIyELl2qf95i338PvXtbF/3zzz+2/n//g6uugt/9Dt55B1q3Dt05K0NE6t1jwIABUikTJ4rk5JRcl5NjrQ+nggKR7dut1zt3ithsIlFRIv37i4BIhw4iS5aE5FRFRUWy9detsvCnhfL2V2/L371/l/vn3X90+2vLXpMHP35Q/rHkHzJv1Tz5YecPcjj/cJXOVVBYIGt3r5X1e9aLiMiew3vkymlXSo/ne4h5wkhkaqSc/tLp8vqy10VE5HD+Yfngxw/kh50/SF5BXrX/1upavHGxTJg/QRZvXFxr5ywsKhTfJp88+PGD0uP5HvLykpdFRGTXoV3y+rLX5ZeDv5T7HvM3zJcrp10p2/db36mfD/wshUWFoUlgQYHIqlUiRUXW8rx5ImPGiIwYIXLZZSJnnCHSpcux/W+9VYpAlnRBnncYWdYJka5dj20fMsT6jgc/zjrr2PZLL7XWRUVZx/XvLzJ69LHtb74pkpkpkpUlsnChlbZ9+46ldeFCkZdeErn9dpF+/azf1vjx1va9e0UGDxZJThaZOVNk/fpjf1cg7dKkicgdd4h8+23VPq/9+0Ueesj6nERE8vJECkv5XxQVibz8spW+c84RWbu2auerIKw5Yk64ptbN27JQi4uz7gDcbvj4Y/j1V3j3XSubJhKaO46KKiiA+fOtO//Zs+HssyEnB2Ji2PPhTGa32sLHWz4nYnd7Wn63lrEdbZwFrFr2Cdn5P9KyWVtaRbWipb0lreyt6NexH82bNCe/MJ/DBYfZuG8j6/asO/qYNHgSdpudv378V15c8mKJpHRt1YVJWQdpsuxLci/YwrROOzliO5ZF7da6Gxt33w7NmvGiLZd9TYo4rXlnep1xEafFXUXrJi1psW4T0qoVj658ge9/XcuqXatYu2ct/kI/d19wN5m/y6R1VGsO+A9wSfdLuLP9nfgL/azevZr2zdoDsGrXKq6Zfg0AESaC7m2606d9Hx699FGcvZz8mvcrG/dtpHOrzrRt2pYIU3MlmiJC3w596dK6C11bdw35e+88tJPNv26mSIqI7WwNC3/vB/fy4eoP2fzrZppENGFw78H0ad8HgPbN2vOHC/5Qofcf2GMgn9z2ydFzJb6faBUdXf4UV59+NaYy3/ONG8HjgWXLrMeKFXDoEOzbZ921ejzw6qvQqRN07mzd5XbqBIWF7MrbyzsjevHPuNNZuW81ILj63sT717wFwJwf5jDw5QyiD2dYd+pHjljPdvux8//lL3DttfDLL7Bjh5V7sNmObX/ySVi7tmSahw6FefOs3/OQIXDgAHTsCAMGwHXXWevAynl8/PHJ//Ynn7RyL6+/Dm+9BddcA8nJHLnwAv675r+s2rXq6GP93vX89bd/ZYxjDIEPHubMgQcegE2brJzO0KEl/7Zgxlg5/7594aabrPqDmTOtTizFMjKsa5jTeWydx2MVa40bV8Y/seJCMjGNMWYI8CLWDGWvi8jTx22PAqYBA4BdwM0isiGwLRn4A1AI3C8iZfyHLJWdmCYjA66M9NA/3WV9mQ8FlaO2aAHJyfDYY1BYCB98YF2cTzut5BcvFJ59FiZOtL7UzZvD1VdbAerGG/Ft8nHZW5eRX5RPt9bdaNakGfvz9pN1cxYXnTqAt67pyp2On094yxX3rOC8U8/jldxX+PO8P5fY1sreim/+9A3d23Tn8w2f882ObzjtSDNO23KIHjf9kaa2KLj4YuuHXVCAHDrI9sJ9rB92Meuud5LvP8KdF90DRUVcMQpyeoEEXUuG9riCeXd+BsDZfwYjcMavkZx59iDOGPJ7Ypv15tyxz1jv36bNscfgwdaPM+Bw/mFWbF/Bmt1rWLN7DWv3rGXN7jVMSJjA5b0u58NVHx4NFDZjI7p5NB2ad+DN694krkscX23/ijk/zKFDiw50aN7h6HOf9n2IiozigP8Am3/dzM6DO9lxcAc7D+1k58GdPHzpw9htdv7vf/9H5vJMdh7cyS+HfiG/KB+AvL/lYbfZeTT7UWZ9P4voZtG0b9ae9s3ac0qLU5g0eBIA/9vyP/Ye2Uv7Zu2xGRubf91MflE+1/e9HoA/zv0j3g1eNv+6mbxCa7rui7tdzKK7FgFw9XtX0ySiCTf0vYHfnfk72jZtW+2vmogw45sZpHhSWLtnLZd0u4QJCRMY1GNQyR3z8+G776yL/fLl8NBD0KMHTJ5sFZW0aGFd5AcMgP794eabrYtbQYH1+wgEFxE5GmgS309kzg9ziOscxx/6/4Gr+lyFiNCrXS/W7VlH75d6E2EiuKjLRQw7fRjDTh/G+aeeX7kAf/iw9TsKfnTocOxiv3Ah9OplBakKBkAR4ad9Px270G/+ilUrvVy+aCvjzvsThyb8nRYTWgBwSotTOCP6DLq06sJ9F97Hpd0v5dvl/+W110ZzZ9YGzjuln/UZFtcDVMSaNVbwGzsW7rrr6OoZ93i4/n0X9iy3FQw8HvyJLmbf7GbEFGcZb3iik01MU+1AYIyxAauAK4HNWHMY3yJBU04aY/4EnCsi9xpjRgCJInKzMeZsYDpwIdAZ+Aw4Q0QKyzpnZQOBx2Ndb3OHjqfn22lsvGkM931yHc/c9T1nFn0Pl10GiYnWHUYf604Mu92q3Onb1yrXHDTIKu+cOxciIo49bDYrenfvbm1fsODYtvx8+OwzeOYZq2xyyhTwesm7/jr+e0YEM9ZkMaDTAMZePJYjBUd43PM4N51zEwM6DSh59yZC3n8/YM9TKez//isO9OnG/qQ72H9pHIN6xtMqqhVLty7Fs95D9zbdOa3dafRq14voZtHW++zZY+WGpk2z6iFat7busqKiys8RiVitGw4d4si+Xfz0yxrWF+1mLXtoY2vOyDXNYd8+ZN8+zK+/WneMV19tlXtu2gTDh1vrih/5+dad1h/+YJXRjhwJ550H55577LlDhxJJ2H5gO59v+JztB7YfvYjvOLSDpxOe5syYM3nzyze5a+5dJyR95eiV/OaU3zBx4UQeyX7khO1bxmyhc6vOvLXiLbJ+yLKCSCCQtGvajjv73wnAG1++wcdrP2b34d1HHwbDugfWgQjXv3stWWs/KPHe3dt056eOT8OePTy8byabCnbTtagl3U45nW7O6zit3WmcmzEV9u61Pl+bzfqOXHKJdWcI1l1ry5bQrt2xR9Om5X/hwbpQb91K/r7dvPHDDFLXvM5W/y7ev/RFXAn3W2X4d98NX39tlZ0DtGpl5VKvuMK6sP7yC5xxRpk3RBv3beTNL99k6ldTybk9h55te/LV9q+IMBH063jiKE2FRYUs3bqUj9Z8xLzV88jdajWhnH7DdEb8ZgT7juxDkJAEw1I/lqICNu7byNrda4/eeHRq1YmxF48FIDojmt2HdwPQ0t6SM6PP5JYzb+DB80ZD27Z8Ofv/6DX+edr+5RHruxv0/5j60l0k7XwTfyT073g+d/a/i1v73Up08+iKJzAvz/pdgnXHf/75eBY24cXhHtzGhf3+0fhfmoxL3Dwwx1kik1ARNRkIHMATInJVYDkZQETSg/b5OLCPzxgTCWwHOgCPBO8bvF9Z56x0jmBRBm0XRXL9I+m822o0t/46mdsdyfzQp4DuG8fhclk3P4f25PFgwgp6HPqeHge/o8eh74ku+IZ1t01g6D9uYce/c2h9YwJNC0q+/4Kxcxj4zHXsnDqPDndcXWLbQVsrVjz5IZc8MpCp83N45L132Bkzm8Im+4j0R9N18xheu+1RrrjCyn3/5S8npn/CBOvGffEi4cN75nLXhhR6H1zJ386ZzcKYRF54wbpp++wzK1cbbMj2txi37h4i8v3s73Y270beTnbHW9kZdazY4+23oVs3eP996ybmeDNnQkyMlUt+660Tt8+bZ2VwXnnFijfHKx7JedIzwqdzDyMY/LZm9Nn/JfduTCY28ivYvv3o/g/2+4Rl7a+k+6EfGCBLeeCN8+DMM3l0fCRfL/yVA7Y2YAzdDv1IbNNveGDUHvJ3/8LsuVvZtncHb/a5m3z7L9y3/GsSd81mz7nw1amGfRvbUrD/FN7t/ApN8mO4aJeH2A4/cXOSlVP5+/OtWbunPRubnwVAp8PruLz3T9x97c/w88/MeuVnfslvw/TuDwPw3FeXc+7+RWxp6WdLa9jdDDbYOpK5/0mKulzNp19dTKcjG0p8FvNjEhl/zmwApub2Jdq+n5btoyjy53Nk2x6yT7mFZ8/IxEgRn81vgo2SrUlmdnmAf/R5gSZFeTz3VQKdThViovZTsHc/B7ft573uD/N+t4foeOQn3l/S8+hxhyPhtQFQsGciczuN45fOD9OseSantezGeX36EtXyt7w300GrA78tcb6j373F8Oij1roi4+eXmDls7/xP9rb/FID+ba/AfPIsLQ+WvPhPmWLdT/3nP1aG+HjPZ+7gm8Mfk/ftMN7JjGZz1xdZ2/tBWv/qoM2+gdx+c2tOaduC5t+N5u2pkRxs/h1++3ZsRS2wFbbgrdeaE9O6BbOmdTz63Ssyfo403cDhZmu4++E1FBQVULRoDB98AMsHxLG/9bFrR7PIZvzuzN/xm+/fJzsbfomZQ2RBO5odOoNOrU5l9izrJik52Zqq4vy9XkavfZAzDyxnb9NTaTvmLujXj798MYKvviyiOT+wqmsO2zu9yYFWy2nXtB3bx27nvtF2Vq0q+beffz688IL1euRI2Lz52LbovK1MX9oH+8Dfwr/+xaDEaK5aNJ7HitKY1CyFAR+mVjoIwMkDQSjqCLoAm4KWNwMXnWwfsSa73wdEB9Z/cdyxpVbVG2OSgCSA7t27VyqBcTsice0ey5o/TuKZ18Yw5/y2rLx0LJctm8SuoP1+JY/5p23hYMv1HGi5moMtv+NIsw38vw7NGAp4Tj3IiL9BE38bmuWdYj2OdOCBM7oyEPipXx9uGzyNFkeiaX6kPU0Ko/j+FD+PnmN9HNNWvcjODl5idiZyyo4RtN2TQIQ0qfgfYgyLYq5jcfTvuHjXf/BFW8Ul7ea9C1vbQZOhnLF/OYN/noa3g4tv2lzCjy0HsG/EvbR7YBTeLRcw47larA8pJf15tuZHF9e06k9K7H+tVns7dvDOuK/Zt+Br1rS0WlU4dv2H0evGwbmAzcaTRUKEFHHtxb/wa5Norto+lZGb0mEBNAFuNDb2R7ZjXrs38duacShiH1va9SO21SHO3rqPzav2kZd3gJkxnQG4etvrxK+cCTlWeh4Hdto7c5NjCwAPrrqH2P99ZuVZgetMJN+2dhwNBEvbDebwby7kkus78tMHv3JFTjqR7OQu7ue2wnP4c38fVwyO4E9/sUNUFJcPsVNkjt1d3x73/dGbkCOHYNgwrBxYwJ/7+7jeuYfBF+7lwKY9vPfKHla1tIrU7EVH8EdEUWA30PsU8iJa4fmiJetbWBfivU06kHHG61w+vBWx8a3YubslXz/fip9jrN/OYVtHNrXuy9qYdbh3rYRdbpqcdwoXL7KKH9f2fojDzVbz0o+nsTzyNI7s7MWRqHNpmteNgsi9/HD277H7O5F05niSh9zJ6qU9eDKr8l+JmGancNsZt/H+D9Zy2z2X033jI+yOnsem7k/z5BfW5/HP7laLni1dX2Zbl1ePHv+b16BpZFOebXcYgB/P/CPbO70JxgqgD/wXerfrzb1YZfhdNv+FoogjNDvch3ZFfcj+dyciTARpaYH0/DL86HuX9ktZ0Taeey5YygV7c7hjZwZtJ0ywyu8dNyMmgoOcTZctZ9Nly31En/01N4z+BrvNqh9Yee4wWhzox6nb76T5obPK/Fx2RXVmzlWv4vr0bjLuPIPzm/6Be3mTVFJ40Ezmf3vbkrGogHGX1JE6AmPMjcAQEfljYPk24CIRuS9on28C+2wOLK/FChZPAF+IyDuB9f8EPhKRmWWds9KT12dkkEkk9+5O59Ko0SzMe4WxbW6hh2zi634duarPVVzf93pW7VrFmf84kwgTwYX3MEoAACAASURBVBnRZ3Bux3M5r+N53ND3Bs6MOZPVu1bj/tbNtgPb2HZgG1v3b2Xb/m3Mvnk2F3S6gNeXv87d/7n76GmjbFHkFeax/oH19Gzbky2/biG6eTRNIyuYva8IEbjoIisb2aGDlaWPioLnnjvWHK6+ys+3mvp9/TV8841V3NaunVWk0aqVVfS0Z4+1rm1bqxilMhWiBw9axxcXae3bZzXjHTbM2u7zWWXRHTvCqada54kovRz7p9HpdHk1hUgKycfG1nvT6DE5OQQfQs076D/Ihr0b2H1499Fmp2M+HsMnaz9h3Z51HC6wLrKDegzi8zs+B2DF9hX0O6UftogQ16MFEREOFxzmoP/g0eay6/asY9O+TRzKP8TB/IMc9B+koKjgaIX621+9zerdq+nTvs/RR4fmHSpXUV4ZW7daFcvFxTkncdB/kFtn38qHqz6kUAqJ7RxLr7a9+OMFf2Rw78H8fOBnXl36Ki3sLWjRpAUt7C1o3qQ5F/3chDUPjOSmIQcY5L2BPkP+yaLZ/2TV1WNx/3YSzsQxlUpuoy4aKq4jGDzxMd7bNKHEtnZN2/HYwMd48OIHKSwqZPm25Zxzyjk0b9L8JO92cj8f+Jmvf/76aJDYcXAHAzoNYPhZw2lhb1Hp96uw/Hx4803473+tyrKbbrIuWqr2+HwsHBnP513yuWxLEy59x9sgRq4TEX4++DPr96xn+4HtDD9reM1dVBuBnw/8zDtfv8PsH2az+/BuUgalcGu/W1m2dRmxr51wfebd69/l4qe/5bX1mUy46heG9hnK4p9yOee9ZKZdUkDvKZXLEdRkIIjEqixOALZgVRbfKiLfBu3zZ6BfUGXx9SLiMsacA7zHscribOD0UFcWZ2RA5Oke0te46NW2F9/u/JZbuvyNU7bfxlPjuugXW1Wbb5OPhLec+Iv82CPsZN/hwdGt/gcCFVoZizKI6xyHs9exAn7Peg+5W3MZ4xjD4fzDHMw/aOV4/Afp0roLr7/cntPPWs87h8Yw+8c5pAxKwUlqlVqPniwQVLtBtogUAPcBHwPfA24R+dYYk2qMuTaw2z+BaGPMGmAMx3IC3wJu4Dvgv8CfywsCVRF3kxUE3De6+d/d/+ODWz7g3zue40rXag0CKiS8G7z4pYBCBL8UhH1oClU3xXWOwzXThWe9B7CCgGumi7jOcURGRNIqqhWntjyV09qdRr+O/WjfrD3jxkHbfhuYv2khKYNSmLx0MvT0hKoLgaW0XmZ1/VHZnsUTF06UnHUlexbnrMuRiQvD3LNYNRiLNy6WZk82E9vfbdLsyWa12itZ1S8563IkJiNGUnJSJCYj5oRr08n2L97v+OXK4CQ9i0PSoay2VbqyWKla4Nvkw7vBS3zPeC0WUmUa7xlP2vw0UgalkOpMLXPfsoqTKttqqMbqCMKh1gJBLXTtVko1LsXFQaNjRzN56WTcN7pLXORrUo3VETRoxWMUeazyvKPNj+Liwpsu1bD4fJCebj2rBq04CLhvdJPqTMV9o7tEnUG4NI5B56rK6bS6yrpc1njlkydby1Xp0qdUaXw+SEiwhpmw2yE7u0E0O1Wly92aWyIH4OzlxH2jm9ytubWWKyiN5gjK43RaQSAtzXrWIKBCyeu1gkBhofVcPB6HapDGLQLnhpLrnBus9eGkgaAMGRnw5XMeKyeQkgKTJ/Plcx4yMsKdMtVgxMdbOQGbzXoOHn5YNTgz1sbhTyxZ3OxPdDFjbXiLmzUQlOHKSA/dxrr4MtkNqal8meym21gXV0aGtzxPNSAOBytfyMabkMbKF7RYqKHrOMKJS9xWMBg/Hn+iC5e46TgivCUNWkdQhv4FuXw5yc3gdCej98LkyU4+meSmf0EuoEVEqvp8Pkj4iwO/34F9AWT301jQkDmdwBwnL109mrFpabzULIUHPqz8cNKhpjmCMmRcAnsTS1YR7E201isVClpF0Pg48TDaTCaVFEabyTgJfwmDBoIyxHWOI3G6ixfnekhJgRfnekicbnUHVyoUtIqgkQnUCfw+0k1BSiq/j3SXrDMIEw0EZdngRNxujMsFzvEYlwtxu2GDFgup0HA4rBajaWnacrQxWDsj9+jsYqmp8MAcq85g7YzcsKZL6wjKkJsLc5534mH00e7gToeT3FxtRapCx+HQANBYzOo9jgfmHLt+FNcZzMp1Es6xCnSIiXKEszu4ahx0jCJVW2pyqsoGK7g7uLOXE2dPZ4llparLt8lHwrQE/IV+7DY72aOyNRjUBw1sHDKtIyhDWd3BlQoF7wYv/kI/hVKIv9Cv8xjUFw1sHDLNEZShtCFenb2cmhtQIRPfMx67zX40RxDfM77iB/t8VnvT+HitZKhtDWwcsmoFAmNMe+B9oCewAXCJyJ7j9jkfmAy0BgqBp0Tk/cC2t4DLgH2B3e8QkRXVSZNS9Ymjm4PsUdmVryPQwerCL3gcspSUehsEoPpFQ48A2SJyOtZ8w4+Uss8hYJSInAMMAV4wxrQN2v6QiJwfeDSMIJCRcWK7YI8HHaRIlcbRzUHywOTK1Q1oT7SwyliUgSfruRLjkHmyniNjUf38jVc3EFwHTA28ngoMP34HEVklIqsDr7cCO4AO1Txv3dbAyg9VHaQ90cIqbkckriVj8WQmQ2oqnsxkXEvGErejfpa2VzcQdBSRbYHX24GOZe1sjLkQsANrg1Y/ZYz52hjzvDEmqoxjk4wxS40xS3fu3FnNZNesjFynNVCdyxpYCpc1cF1Gbv3NOqo6RnuihZV5rYCn2k/CtSad8Z7xuNak81T7SZjXCsKdtCopNxAYYz4zxnxTyuO64P0CEyOftFOCMaYT8DZwp4gUBVYnA2cBcUB74OGTHS8imSISKyKxHTrU7QzF2lMziJ8GG4Za5Ycbho4mfpq1XqlQ8eEgnWR8aBCobfLQOB57ZgxDY6zOpkNjRvPYM2OQh+pf01GoQGWxiFxxsm3GmJ+NMZ1EZFvgQr/jJPu1Bj4EHhORL4Leuzg3kWeMeRMYW6nU11EjBsbh3pDIireEDQNT2LvwRSLueIkRA7PCnTTVQGhdcXg5nZCc6WHskskMjErhnR8nMynTibOeVhhXt2hoLnB74PXtwL+P38EYYweygGkiMvO4bZ0CzwarfuGbaqanTnBugNmzhdtGGpyRcNtIw+zZcsLMREpVldYVh5dnvYf0NS5G2t0sSE1lpN1N+prwzz1cVdUNBE8DVxpjVgNXBJYxxsQaY14P7OMCBgF3GGNWBB7nB7a9a4xZCawEYoAnq5meuiE3l7aj5lC47H64LI3CZffTdtQcq9ehUiGgdcXhlbs1l+Q+bj6a7CQlBT6a7CS5T/3tbKpjDdUAjweG/9WDcbm43zGal3yTEbebOc+HfwIK1XBof7LwKW4IWNyH7PjluupkYw3pEBM1YMYXVhDIusVNqjOVrFusoaxnfFE/s42qjurqg0vTrWdVq3JzS170izsa19dMv+YIakDGogziOseVGIrCs95D7tbcUoetUKqydLA6VRU6+mgt0jGKVE0rbbA6DQSqqrRoSKl6qHiwOpuxVX6wOqWOozkCpeqhKg9Wp1QpNBAoVU85ujk0AKiQ0KKhukRHLVVKhYEGgjpkxto4/IklRy31J7qYsVZHLVWqJgx7KoPnskrefD2X5WHYU43r5ksDQR3ScYQTl7itYDB+PP5EFy5x03GEtjZSoePzQXq69dzYXXF2HGOXuI4Gg+eyPIxd4uKKsxvXzZcGgjok157BoLfgJb81aulL/tEMestar1Qo+HyQHO/jwGPpJMf7Gn0wGLM6l1fbJzN2iYtBj49n7BIXr7ZPZszqetozrIo0ENQhcZ3jSP8+kbN7vEgqKZzd40XSv08krnPjujtRNWf1NB/z/An8XVKY509g9bRGHgni4kh6Jp1rDg9lQUQa1xweStIz6Y1uEikNBHWIcwO8+56QeK3hs/GQeK3h3fd01FIVOpfhxY6fSAppgp/L8IY7SeHldJL5UDKLo97htvUDWRz1DpkPJdftAYNqgAaCOmTtjFxe2TSHm8+6nwURadx81v28smkOa2c0rmyqqjk9RsVjouwUGhsRUXZ6jIoPd5LC6rksD/fuTuflvSOZNnUBL+8dyb2700+oQG7oNBDUIbN6j2PQW/DRL5NJGZTCR79MZtBb1nqlQsLhwObJxvZUGjaPzmbz2XdWHcEtWR9BSgq3ZH3Eq+2T+ey7xnXzpR3K6pC4mzy4Zrpw3+i2xibq6bSWb3IDjSurqmqOryt4L4X4rjT6SS7nXRxXcvxop5Mkl4sktzvcSatV1coRGGPaG2M+NcasDjy3O8l+hUGT0swNWt/LGLPEGLPGGPN+YDazRit3a+7RIADWQHXuG8uZ7EI7oalKKB61NMWTQsK0BHybGnllcUMbT7qKqls09AiQLSKnA9mB5dIcFpHzA49rg9ZPBJ4XkT7AHuAP1UxPvTbuknEnjFDq7OUse+jquMAdTVAnNFyuRtfqQVVMaaOWVlhD7IAwbtyJFcNOp7W+EaluILgOmBp4PRVr3uEKCcxTfDlQPI9xpY5XAcV3MC6rE1q9mCZJhU2VRy31+SAhAVJSrOe6Fgw0Z1wt1Q0EHUVkW+D1dqDjSfZraoxZaoz5whhTfLGPBvaKSEFgeTPQ5WQnMsYkBd5j6c6dO6uZ7IYjY1EGnp7AaKsTGqNH4+lprVfqeMWjlqY50yo3mY3XC34/FBZaz15vTSaz0jJOXYvnr8NL5Iw9fx1Oxqlrw5uweqLcymJjzGfAqaVseix4QUTEGHOy6c56iMgWY8xpQE5gwvp9lUmoiGQCmWDNUFaZYxuyuM5xuKYn4p4rOFNS8Mx9EVfTl3DfkhXupKk6qkqjlsbHg91uBQG73VquQ+IGjsC1+X3cfx2O89oHrN+By+AeOCLcSasXyg0EInLFybYZY342xnQSkW3GmE7AjpO8x5bA8zpjjBfoD8wC2hpjIgO5gq7Alir8DY3az0/Du4sE162G0Q6Y3NTqhPbzT8CUcKdONRgOB2RnWzmB+Pg61+w0919Okvtm4Rp+DaPnpzF5eHOS+35A7r+cOBtXcX+VVLdoaC5we+D17cC/j9/BGNPOGBMVeB0DXAJ8J9ZkyR7gxrKOV2WLw+qENrTT/aTNT2NoJ6sTWhyNq9WDqnk+HKSTjK8ONjqNi4P5d0DSEiHtMut5/h3aZqKiqhsIngauNMasBq4ILGOMiTXGvB7Ypy+w1BjzFdaF/2kR+S6w7WFgjDFmDVadwT+rmZ5Gp/cUqxPaOz9OZmBRCu/8aHVC6z1Fb4NU6NT1umInHv7UbTiTzrUzsCiFSefa+VO34ThpXD2Eq6pagUBEdolIgoicLiJXiMjuwPqlIvLHwOvFItJPRM4LPP8z6Ph1InKhiPQRkZtEJK96f07j41nvIX2Ni5F2NwtSUxlpd5O+xoVnvf4AVOh4vZDXwUehI528Dr66VleM54sZ/P5Ww81RWSxITeXmqCx+f6vB88WMcCetXtAhJuq53K25JPdx89FkJykp8NFkJ8l9yumEBtrcTlVK9Pk+im5LgMtTKLotgejz61aWIHdQb5L7ZpX8HfTNIndQ73AnrV7QISbquTj/OFxJJXrI43I5cbvL6UcQd1zX+uKOaI2sa72qmF0tvUQ08VNEIRERfna19FKXBqio8u9AAZojqPeq3ENeO6KpSojvGU9UpNURLSqyEh3RaomOFFE9xmq8U7/ExsbK0qVLw52MhmH8eKsjWkoKpKaGOzWqDvNt8uHd4CW+Z3zl+yFUVEaGlVsNviHxeKwreiMb9qEmGGOWiUjs8es1R9CYeTwwebIVBCZPPrHOQKkgjm4OkgcmVy4IVHZ8Ih07Kyy0jqCxCq4TOFaoqsVDKnSK25wW90bOLn/+gwx7LnGZyThdLmvYlMmT8WQmk2vPZZwOxV5jNEfQWGmhqqppVRifKK5zHK416XhGD4W0NDyjh+Jak67zdtcwDQSNVAbj8Bx3h+XBSQZaDqtCpHh8IputwuMT5f7LyVOSjMv/DuPHD8Tlf4enJJncf2luoCZpIGiktChW1TiHg5UvZONNSGPlCxWbFvPKSA/XP5JOYpuRpEUsILHNSK5/JJ0rI7X+qiZpq6FGrPjiHyiK1eoBFVJVqCKAjAwyieTe3elcGjWahXmTebV9MkkUaKuhENBWQ6qEjEUZ0NMTPI0B9PToPAYqZKoyhYHnpjgeM+klhkx5zKTjuUmzqjVJA0EjFdc5jsTpLl6c6yElBV6c6yFxuksr5VTIVKGKoOpDpqhq0eajjdUGJ+J2Y1wucIzGNJ2MuN3gcEKvcCdO1XcZizKI6xxHdrbz6BQGR071kLEot8w5uHWoiPDQHEEjlZsLc553cr9jNGnz07jfMZo5zzvLbz2qg9WpCoibvxbX9ESW2f4Bl6azzPYPXNMTiZtf9tSR2qo5TESk3j0GDBggqvpy1uVITEaMpOSkSExGjOSsy6nAQTkiMTHWc2nLSomI5OTIy87mwhOIecIITyAvO5vr9yTMgKVSyjVVcwSNlGe9B9dMF+4b3aQ6U3Hf6MY1swLzGOhgdaoinE72334rRkAQjMD+22+t+PekskNTqGqpViAwxrQ3xnxqjFkdeG5Xyj5OY8yKoMcRY8zwwLa3jDHrg7adX530qIrL3ZqL+0Y3zl7WD9PZy4n7xgpWyjmdlGhupEFAlaJV3/MQwBSBBJYrpK5Ph9YAVbey+BEgW0SeNsY8Elh+OHgHEfEA54MVOIA1wCdBuzwkIjOrmQ5VSaVV2Dl7OY8GhpMZ9lQGw5tEkhQ0WF1my7bMyS9g3mPazltZ7kn34M7/Gy/Pb87+gRfRav4SUiL+xsrsc5iSXM6NQ2ntTivQGU1VXXWLhq4DpgZeTwWGl7P/jcBHInKomudVYTK8SSSP/jKWzIeSITWVzIeSefSXsQxvog3Q1DFNmcGbbwmXXPMByX/P4ZJrPuDNt4SmVGDqyPh4CiPtFBobhZEVbHeqqqW6v96OIrIt8Ho70LGc/UcAzx237iljzHggG3hETjJvsTEmCUgC6N69e9VTrKoliQKImcS9u9N55/G9Vs/PmEnWeqUCXrT15sv75jA43cnovTB5spNPkufwYkH5RY8+HPyp8wu07TqLvZtv4BUcdWgutIap3CEmjDGfAaeWsukxYKqItA3ad4+InFBPENjWCfga6Cwi+UHrtgN2IBNYKyLlzo6iQ0yE36DHx7MgIo2BRSnM/7tOaNOQFfcJCC429Kz3kLu17D4BULV5j0Y/5ePVwwlg80OhnXubZTP5MQ0FoVDlISZE5AoR+U0pj38DPwcu5sUX9R1lvJULyCoOAoH33hZo1ZQHvAlcWNk/TNW+57I8LMybzMCiFBbmTea5rAoMCKb9D+qt4j4BxS3KPOs9FeoTUOV5j3p6rSAQUQgRfmtZ1ajq1hHMBW4PvL4d+HcZ+94CTA9eERREDFb9wjfVTI+qYc9leRi7xMWki9zM/3sqky5yM3aJq/xgoMOd1lvO347A7RZc0xMZ7xmPa3oibrfg/O2Ikx4TPO9RauqxFscVCQajBlnzIxtsRDWxM2pQfOj+GFW60joXVPQBRGOV7a8GPgPaB9bHAq8H7dcT2AJEHHd8DrASKwC8A7SsyHm1Q1n4DH1yojw7u2SnoGdn58jQJyeWf3Bx57OUFO2EVt/k5EjKsGbCE0jKsGbl/u8mTjxxl5wca31FLN64WCbMnyCLNy6uYoJVaThJhzIdhlrVrqoUGquw86z34Jp6DaPnH2LyoOa4b/+g3KbGYeHzcXRwI21yeoKT1RFomz9VKzIyrElH+gcVGn/Z1smnBU4dZr6Os/oEJDJ7TiTOa1NwznmR6yMScTXJKr9PQG3y+fDdFo+3Sz7x/2yC422vBoMK0iEmVK24MtJDt7Euvky2Co2/THbTbaxLZ56qTVWssC/uE9B21BxITaXtqDkV7xNQi3w500gY4SclXkgY4ceXMy3cSao3NBCoWtG/IJdNk9wMTncyfjwMTneyaZKb/hVoV65CpIoV9i/aetMj0Ceg+H/X4745vGjrXQuJrrgZzcBvg8II8EdYy6piNBCoWpFxCexNLDlE0d5Ea33ZB2qz01DJsOfiyUwuMWCgJzOZDHs5wXjcOPqPcZb43/Uf46xzU0f6D49CCqMwhQYpisJ/eFTFD27sg9yVVoNc1x/aaqj+yVmXI22eipHW5+VISopI6/Os5XKHvtZhr0OmeNjxnJTbREByUm6r8PDj9aHB1+LFIvbei8UMmiD23otlcUUbHC1eLNKsmYjNZj1X+MD6h5O0Ggr7Rb0qDw0E9U9OzrGLf0pOytGgUKELSn24CtUTObOflZiHjaSMHygxDxvJmf1s+cfUo1i8eLHIhAmVvJZPmGAFAbCeJ0yosfSF28kCgRYNqVpR1RnRMhZl4OlJiTIlT09rfdkHapHSCTwenEnpjD5zJGkRCxh95kicSenl9vKqT7OGORyQnFzJxkJVmVy5oSktOtT1h+YI6qeqzIiWNCFH2v69jeSc11okJUVyzmstbf/eRpImaJFSZXmGTpQpE58t8T+YMvFZ8QytYC+vhqxKWYn6B80RqHCq6oxo90bBm28J119vGO+E6683vPmWcG9UOSdswDOpZSzKOOFz86z3lJtLWn53HPfuTie5j/U/SO7j5t7d6Sy/W4f5yNwLVxVaz42RBgJVK6o6I1r/glx63DeHvCX3kzY/jbwl99PjvjnlNjutcpFSPRDXOa5EEC0OsnGdy76gF5ySy6SL3KQnWc1A05OcTLrITcEpdbCMpxZlfuTjnkUJfJKfwj2LEsj8qBG2HCotm1DXH1o01LjkrMuR5o/HCM4Uaf54BVu5rMuRmKdKFinFPNWmQsdWWnUH1qnC+XJmlyziyZn9bIXPl5Ji1YumpNRM8uqbwWkThPE24QmEFJsMTtPKYqXqFM96D4nTXUTOcZNyaSqRc9wkTi+/SMm5AWvETJdVpORyGWvEzA01kMjaHlk1Ls6q9I0ZStr8NEbHDLUqfStwvioPDd2A3TAgHgrtUGiDIru13NiUFh3q+kNzBI1H0tSJJZqZFjdDTZpa9t2vZ+hEWf5sjqTkpFgjZuakyPJnc8quGK3OnX1VmrhW8XwTJ4pV6RvUDHTKxGfLTabWn5/clHmLZXDaBJkyr3KVxYvnTZEJaYNl8bwpNZSy4JNVv0Ib7Ueg6qOqXpur1G+hulfKSpa5TJyaZBVdBZ0v57zWMnFqUpnHPTs7R8zDMfLe3VbHsPfuvk3MwzEnDA9+wvlquQSroVs8b4o0ewyxjUeaPUbNBoMQdXo7WSDQoiFVp40bd2JDH2dFRjfo6cG4XIjbDZ5UxO3GuFzQ8+RlIVUeggGqVOYSN3AELpfB89fhMH48nr8Ox+UyxA08+YQvYFX6vto+mStf/wjvwBSufP0jXm2fXG6lb5U/S1WqrOxZJcY2ysqeVXMn83rB74fCQuvZ6w3p22sgUA1S7tZcsm5x88C1TtLS4IFrnWTdUnYrpbXz40j8Ph3P6KGQloZn9FASv09n7fxyyt6rOB1X7r+cJPfNwjW8gPEL03ANLyC5bxa5/yq7ies4fxxJz6Tzn5FunAtS+c9IN0nPpDPOr81Aa1NT/w3YC8FWCPYia7nG1HSnt9KyCRV9ADcB3wJFQGwZ+w0BfgTWAI8Ere8FLAmsfx+wV+S8WjSkKqKyxfY5OSIXn/OsxIwLlL2PM3LxOc+We1zSXUNlysSSQzVMmfisJN01tNzzXdc6Rx4dYs389eiQZnJd6woMuzHRqv8I/tuWP6tlPLVt8WKR83tNkfiBg+X8XlMqVVrz9ZTF4hk8Qb6eUvGDQlEfQU3UEQB9gTMB78kCAWAD1gKnAXbgK+DswDY3MCLw+lVgdEXOq4FAladKxf05OZLXJkbOGHKb8ARyxpDbJK9N+RGkuMy+uIz++OWyzvfxOa3FPq6NDByfIvZxbeTjc1qXez6t9K07qlJ/+/WUxZLT1S5PXmokp6u9QsFg8cbFEpXaTMwTNolKbVblKTxPFgiqVTQkIt+LyI/l7HYhsEZE1omIH5gBXBeYsP5yYGZgv6lYE9grVW1VGh8nN5dFbyaz+aKP4PMUNl/0EYveTC53UJ0xiVbHrLFLXAx6fDxjl7iYdJGbMYllF/F4vpjB72813ByVxYLUVG6OyuL3txo8X5Q94Ut9GvunoavK2EbzP57G1bf7efxy4erb/cz/uPwJdKbN95JX4EcoJC/fz7T53qonuhS1MVVlF2BT0PJm4CKsie/3ikhB0PouJ3sTY0wSkATQvXv3mkmpajBKqwB1OsseYcJzU9yxPgvXOnlxjpNEu4usW9yUNzDFmEQnc1aMZkFEGgOjUsoNAgC5g3qTvCOL9CRnoI7ZSXJmFrmn5JZ5vqr8baruWHUJ+PcFKpnFWi7Xhnirr4P4ochuLYdQuTkCY8xnxphvSnlcF9KUlENEMkUkVkRiO3ToUJunVo3EjAW5iNvNnOedpKZao6WK282MBeXfaj+X5WFh3mQGFqWwMG8yz2VVoNWQfxzpSc4SdczpSU7i/NqMpyEbcdMomkREEVFkaGKLYsRN5U+gM+pyB/YZ2ZjP07DPyGbU5aGdi7ncHIGIXFHNc2wBugUtdw2s2wW0NcZEBnIFxeuVCove28cx5/mSRS7WUNll32o/l+UpURz0XJaTsUtcQNnFQ2UV8ejdfcPl6OYg5y4P3g1e4nvG4+hW/kXd4QDv2w68Xgfx8ZUcZrsCjFV/UM03McYLjBWRpaVsiwRWAQlYF/pc4FYR+dYY8y9glojMMMa8CnwtIq+Ud77Y2FhZuvSEUykVFsOeyuCKs+NKXPSfy/Lw2Xe5zHtM7+5V3WGMWSYisSesr04gMMYkAi8DHYC9wAoRucoY0xl4XUSGBfYbBryA1YLoDRF5KrD+NKzK4/bAl8BIEckrTO5UlAAABSlJREFU77waCJRSqvJqJBCEiwYCpZSqvJMFAu1ZrJRSjZwGAqWUauQ0ECilVCOngUAppRq5ellZbIzZCfxUxcNjgF9CmJyaUl/SCfUnrZrO0Kov6YT6k9aaTmcPETmhR269DATVYYxZWlqteV1TX9IJ9Setms7Qqi/phPqT1nClU4uGlFKqkdNAoJRSjVxjDASZ4U5ABdWXdEL9SaumM7TqSzqh/qQ1LOlsdHUESimlSmqMOQKllFJBNBAopVQj16gCgTFmiDHmR2PMGmPMI+FOTzFjTDdjjMcY850x5ltjzAOB9U8YY7YYY1YEHsPqQFo3GGNWBtKzNLCuvTHmU2PM6sBzuzCn8cygz2yFMeZXY8xf6srnaYx5wxizwxjzTdC6Uj9DY3kp8J392hhzQZjT+Ywx5odAWrKMMW0D63saYw4HfbavhjmdJ/1fG2OSA5/nj8aYq2ornWWk9f2gdG4wxqwIrK+9z7S0iYwb4gNrCOy1wGmAHfgKODvc6QqkrRNwQeB1K6z5G84GnsCa5yHsaQxK6wYg5rh1GcAjgdePABPDnc7j/u/bgR515fMEBgEXAN+U9xkCw4CPAAP8FlgS5nQOBiIDrycGpbNn8H514PMs9X8d+F19BUQBvQLXBFs403rc9meB8bX9mTamHMGFwBoRWScifqx5EGp1us2TEZFtIrI88Ho/8D1lzN9cB10HTA28ngoMD2NajpcArBWRqvZEDzkRmQ/sPm71yT7D64BpYvkCa1a/TuFKp4h8IsfmGf8Ca2bBsDrJ53ky1wEzRCRPRNYDa7CuDbWirLQaYwzgAqbXVnqKNaZA0AXYFLS8mTp4sTXG9AT6A0sCq+4LZMPfCHeRS4AAnxhjlhljkgLrOorItsDr7UDH8CStVCMo+cOqa59nsZN9hnX5e3sX/7+9s3eNIoqi+O8gahE/QLEQVNhIrFUsLBIrC1c0oDYRwQg2go1Y2Oz/YCcKIggSQUTFrbWwVYzRRPyIWBmWDaTQwsaPa/HehNlNNuW8J3N/MOzs3Vk4nPfgzLvz2A2rlYKGpDeSXkgaSyWqxGpjnbOfY0DXzOZLtUo8rVMQZI+kTcAj4IqZ/QBuAnuB/UCHsGxMzaiZHQSawGVJR8ofWljTZrEnWdIGYBx4GEs5+rmCnDwchKQW8BuYiqUOsMfMDgBXgfuStqTSx38y1n2cpfempTJP6xQEC8Du0vtdsZYFktYTQmDKzB4DmFnXzP6Y2V/gNhUuYQdhZgvxdRF4QtDULdoV8XUxncIemsC0mXUhTz9LDPIwu3kr6QJwAjgXQ4vYalmK568Jvfd9qTSuMdbZ+QnL/+1+GnhQ1Kr0tE5B8AoYkdSId4oTQDuxJmC5N3gH+GBm10v1ci/4FDDX/90qkTQkaXNxTnhwOEfwcTJeNgk8TaNwBT13WLn52ccgD9vA+bh76DDwvdRCqhxJx4BrwLiZ/SzVd0haF8+HgRHgaxqVa451G5iQtFFSg6DzZdX6VuEo8NHMvhWFSj2t6ml5DgdhB8ZnQrK2Uusp6RoltALeATPxOA7cA2ZjvQ3sTKxzmLDj4i3wvvAQ2A48B+aBZ8C2DDwdApaAraVaFn4SwqkD/CL0qC8O8pCwW+hGnLOzwKHEOr8QeuzFPL0Vrz0T58QMMA2cTKxz4FgDrejnJ6CZeuxj/S5wqe/ayjz1n5hwHMepOXVqDTmO4zir4EHgOI5TczwIHMdxao4HgeM4Ts3xIHAcx6k5HgSO4zg1x4PAcRyn5vwDVeOnRoZOKboAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(thetas * 180 / np.pi, ideal_Xexps, 'b--')\n",
    "plt.plot(thetas * 180 / np.pi, ideal_Yexps, 'bx')\n",
    "plt.plot(thetas * 180 / np.pi, ideal_Zexps, 'b.')\n",
    "\n",
    "plt.plot(thetas * 180 / np.pi, std_Xexps, 'r--')\n",
    "plt.plot(thetas * 180 / np.pi, std_Yexps, 'rx')\n",
    "plt.plot(thetas * 180 / np.pi, std_Zexps, 'r.')\n",
    "\n",
    "plt.plot(thetas * 180 / np.pi, opt_Xexps, 'g--')\n",
    "plt.plot(thetas * 180 / np.pi, opt_Yexps, 'gx')\n",
    "plt.plot(thetas * 180 / np.pi, opt_Zexps, 'g.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### plot trajectories on Bloch sphere"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from qutip import Bloch\n",
    "import numpy as np\n",
    "from math import pi\n",
    "b = Bloch()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "b.clear()\n",
    "b.zlabel = ['$z$', '']\n",
    "b.point_size = [1, 32, 20, 20]\n",
    "b.add_points([0,0,0])\n",
    "b.add_points([std_Xexps, std_Yexps, std_Zexps], meth='l')\n",
    "b.add_points([opt_Xexps, opt_Yexps, opt_Zexps], meth='l')\n",
    "#b.add_points([ideal_Xexps, ideal_Yexps, ideal_Zexps], meth='l')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### for a sample point, plot gate vs. pulse trajectory of Rx(θ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "theta:  122.72727272727273\n",
      "        ┌──────────┐┌──────────┐┌────────────┐┌──────────┐┌──────────┐\n",
      "q_0: |0>┤ U1(pi/2) ├┤ Rx(pi/2) ├┤ U1(5.2836) ├┤ Rx(pi/2) ├┤ U1(pi/2) ├\n",
      "        └──────────┘└──────────┘└────────────┘└──────────┘└──────────┘\n",
      "        ┌───────────┐\n",
      "q_0: |0>┤ Rx(2.142) ├\n",
      "        └───────────┘\n"
     ]
    }
   ],
   "source": [
    "sample_point = 15\n",
    "print('theta: ', 180 * thetas[sample_point] / np.pi)\n",
    "print(standard_paths[sample_point])\n",
    "print(optimized_paths[sample_point])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# standard trajectory\n",
    "b = Bloch()\n",
    "b.zlabel = ['$z$', '']\n",
    "b.point_size = [1, 32, 20, 20]\n",
    "b.point_marker = ['o', 's', '.', '^']\n",
    "b.add_points([0,0,0])\n",
    "\n",
    "z1_path = [[0],[0],[0]]\n",
    "\n",
    "x1_path = [ideal_Xexps[:21],\n",
    "           ideal_Yexps[:21],\n",
    "           ideal_Zexps[:21]]\n",
    "\n",
    "z2_angle = standard_paths[sample_point][2][0].params[0]\n",
    "z2_path = [list(np.sin(np.linspace(0,z2_angle,20))),\n",
    "           list(-np.cos(np.linspace(0,z2_angle,20))),\n",
    "           list(np.zeros(20))]\n",
    "\n",
    "z2_end_x = z2_path[0][-1]\n",
    "x2_path = [list(np.full(20, z2_end_x)),\n",
    "           list(np.sqrt(1-z2_end_x**2)*np.cos(np.linspace(0,np.pi/2,20))),\n",
    "           list(np.sqrt(1-z2_end_x**2)*np.sin(np.linspace(0,np.pi/2,20)))]\n",
    "\n",
    "x2_end_x = x2_path[0][-1]\n",
    "x2_end_y = x2_path[1][-1]\n",
    "x2_end_z = x2_path[2][-1]\n",
    "\n",
    "z3_path = [list(np.sqrt(1-x2_end_z**2)*np.cos(np.linspace(np.pi, 3*np.pi/2,20))),\n",
    "           list(np.sqrt(1-x2_end_z**2)*np.sin(np.linspace(np.pi, 3*np.pi/2,20))),\n",
    "           list(np.full(20, x2_end_z))]\n",
    "trajectory = [x1_path[i] + z2_path[i] + x2_path[i] + z3_path[i] for i in (0,1,2)]\n",
    "b.add_points(trajectory, 'l')\n",
    "\n",
    "# optimized trajectory\n",
    "z1_path = [[0],[0],[0]]\n",
    "x1_path = [ideal_Xexps[:sample_point], ideal_Yexps[:sample_point], ideal_Zexps[:sample_point]]\n",
    "z2_path = [[0],[0],[0]]\n",
    "trajectory = [x1_path[i] for i in (0,1,2)]\n",
    "b.add_points(trajectory, 's')\n",
    "\n",
    "b.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}