presentation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bayesian Framework for Quantum Algorithms in Qiskit\n",
    "\n",
    "## _Team: Baye's Baes_\n",
    "\n",
    "**Team members:** _Conrad Haupt, Annie Ho, Shawal Kassim, Eric-Muthemba Kiarie, Roman Peters_\n",
    "\n",
    "**Coach:** _Ismail Akhalwaya_\n",
    "\n",
    "\n",
    "This project was developed for the Qiskit Camp Africa in 2019. It shows that Bayesian Inference can be integrated into Qiskit using the QInfer toolkit. This notebook showcases a proof-of-concept for the Iterative Quantum Phase Estimation circuit where the phase of the unitary is determined through a Sequential Monte-Carlo sampling scheme over an arbitrary prior distribution. The experimental results used for updating the posterior come from a QASM simulation of the IQPE circuit at iteration 1.\n",
    "\n",
    "[1] N. Wiebe and C. E. Granade, “Efficient Bayesian Phase Estimation,” Phys. Rev. Lett., vol. 117, no. 1, p. 010503, Jun. 2016.\n",
    "\n",
    "[2] C. Granade et al., “QInfer: Statistical inference software for quantum applications,” Quantum, vol. 1, p. 5, Apr. 2017.\n",
    "\n",
    "[3] H. Abraham et al., Qiskit: An Open-source Framework for Quantum Computing. 2019."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/conrad/.anaconda3/lib/python3.7/site-packages/qiskit/providers/ibmq/ibmqfactory.py:181: UserWarning: Credentials are already in use. The existing account in the session will be replaced.\n",
      "  warnings.warn('Credentials are already in use. The existing '\n"
     ]
    }
   ],
   "source": [
    "# Import major packages\n",
    "import qiskit as qk\n",
    "import qinfer as qf\n",
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "\n",
    "# Qiskit imports for IQPE\n",
    "from qiskit import IBMQ\n",
    "from qiskit.aqua.algorithms.single_sample import IQPE,QPE\n",
    "from qiskit.aqua.components import iqfts\n",
    "from qiskit.aqua.operators import WeightedPauliOperator\n",
    "from qiskit.aqua.algorithms import ExactEigensolver\n",
    "from qiskit.aqua.components.initial_states import Custom\n",
    "from qiskit import Aer,execute\n",
    "from qiskit.aqua import QuantumInstance\n",
    "from qiskit.providers.aer import noise\n",
    "\n",
    "# QInfer imports for modelling the IQPE algorithm\n",
    "from qinfer import distributions\n",
    "from qinfer.smc import SMCUpdater         # The QInfer Sequantial Monte-Carlo Sampler\n",
    "from qinfer import FiniteOutcomeModel\n",
    "\n",
    "# Tertiary imports\n",
    "import matplotlib.pyplot as plt\n",
    "import time\n",
    "\n",
    "# Intial setup\n",
    "IBMQ.load_account()\n",
    "ibmqx_provider = IBMQ.get_provider(hub='ibm-q-wits', group='internal', project='default')\n",
    "mpl.rc('text', usetex=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {
    "hide_input": false
   },
   "outputs": [],
   "source": [
    "# Function definitions for the entire notebook\n",
    "EPS = 1e-15\n",
    "def chop(value):\n",
    "    if np.abs(value) <= EPS:\n",
    "        return 0\n",
    "    else:\n",
    "        return value\n",
    "def radiandiff(a,b):\n",
    "    FULL_ROTATION = 2 * np.pi\n",
    "    return np.min(np.abs([a - b, b - a, a - b + FULL_ROTATION, b - a + FULL_ROTATION]))\n",
    "def complexfromangle(angle):\n",
    "    real = chop(np.real(np.exp(1.j * angle)))\n",
    "    imag = chop(np.imag(np.exp(1.j * angle)))\n",
    "    return real,imag"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# IQPE on ZZ Unitary with specified phase"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we setup the IQPE circuit using Aqua with an arbitrary phase"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set-up unitary with specified phase\n",
    "angle = 4 * np.pi / 5\n",
    "real,imag = complexfromangle(angle)\n",
    "\n",
    "PAULI_DICT_ZZ = {\n",
    "    'paulis': [\n",
    "        {\"coeff\": {\"imag\": imag, \"real\": real}, \"label\": \"ZZ\"}\n",
    "    ]\n",
    "}\n",
    "QUBIT_OP_ZZ = WeightedPauliOperator.from_dict(PAULI_DICT_ZZ)\n",
    "eigen_results = ExactEigensolver(QUBIT_OP_ZZ,4).run()\n",
    "eigen_value = eigen_results['eigvals'][0]\n",
    "eigen_vector = eigen_results['eigvecs'][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigenvalues [-0.80901699+0.58778525j -0.80901699+0.58778525j  0.80901699-0.58778525j\n",
      "  0.80901699-0.58778525j]\n",
      "Phase is 2.5133 radians\n"
     ]
    }
   ],
   "source": [
    "print('Eigenvalues',eigen_results['eigvals'])\n",
    "print('Phase is %.4f radians' %(angle))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setup initial state as eigenvector\n",
    "state_in = Custom(QUBIT_OP_ZZ.num_qubits, state_vector=eigen_vector)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use QASM Simulator for testing\n",
    "simulator = Aer.get_backend('qasm_simulator')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Construct IQPE circuit and simulate\n",
    "iqpe_builder = IQPE(\n",
    "    QUBIT_OP_ZZ,\n",
    "    state_in,\n",
    "    num_iterations=1,\n",
    "    num_time_slices=1,\n",
    "    expansion_mode='suzuki',\n",
    "    expansion_order=1,\n",
    "    shallow_circuit_concat=True)\n",
    "\n",
    "# This is required to setup iqpe_builder so we can call construct_circuit() later on\n",
    "iqpe_results = iqpe_circ = iqpe_builder.run(simulator,shots=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# QInfer Models and Definitions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we define our IQPE QInfer model as well as the required analytical representations for the inference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extend the QInfer FiniteOutcomeModel as IQPE has only one bit per measure\n",
    "class IQPEModel(FiniteOutcomeModel):\n",
    "    \n",
    "    def __init__(self):\n",
    "        super(IQPEModel, self).__init__()\n",
    "\n",
    "    # Number of parameters overwhich to conduct the Bayesian Inference\n",
    "    @property\n",
    "    def n_modelparams(self):\n",
    "        # phi\n",
    "        return 1\n",
    "\n",
    "    # Number of parameters that are needed for the simulation but are known.\n",
    "    # i.e. do not need to be inferred\n",
    "    @property\n",
    "    def expparams_dtype(self):\n",
    "        return [('theta','float64'), ('m','int64')]\n",
    "\n",
    "    @property\n",
    "    def is_n_outcomes_constant(self):\n",
    "        return True\n",
    "\n",
    "    # This function will verify whether the new inferred model parameters are valid for our model.\n",
    "    @staticmethod\n",
    "    def are_models_valid(modelparams):\n",
    "        # NOTE: If T2 is included as a model parameter, restrict it to positive\n",
    "        return np.logical_and(modelparams[:] >= -2*np.pi, modelparams[:]<=2*np.pi, modelparams[:] > 0).all(axis=1)\n",
    "\n",
    "    def n_outcomes(self, expparams):\n",
    "        return 2\n",
    "\n",
    "    # This function implements our probability for reading a 0 or 1 on our IQPE circuit.\n",
    "    # These equations are taken from [1] and are used by the updater/sampler to modify the posterior\n",
    "    def likelihood(self, outcomes, modelparams, expparams):\n",
    "        super(IQPEModel, self).likelihood(outcomes, modelparams, expparams)\n",
    "        # Probability of getting 0 for IQPE is\n",
    "        # P(0) = (1 + cos(M(theta - phi)))/2\n",
    "        # P(1) = (1 - cos(M(theta - phi)))/2\n",
    "        \n",
    "        # This is the probability of getting a 0 with T2 depolarisation noise from [1].\n",
    "        # Not currently used.\n",
    "        # pr0 = np.exp(-expparams['m'] / expparams['T2']) * (1 + np.cos(expparams['m'] * (expparams['theta'] - modelparams[:]))) / 2 + (1 - np.exp(-expparams['m'] / expparams['T2']))/2\n",
    "        \n",
    "        pr0 = (1 + np.cos(expparams['m'] * (expparams['theta'] - modelparams[:]))) / 2\n",
    "\n",
    "        return FiniteOutcomeModel.pr0_to_likelihood_array(outcomes, pr0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 136,
   "metadata": {},
   "outputs": [],
   "source": [
    "# As our model parameter is an angle, ensure we have the circular mean for the distribution so we allow\n",
    "# for wraparound. This is done by converting angles to x,y points on a unit circle, averaging over\n",
    "# x and y independently, and then combining the resulting average x and y values to create a new\n",
    "# vector. If the resulting magnitude is 0, there is no circular average as all points are equi-distant.\n",
    "# The resulting angle is approximately the circular average of the collection of angles/phases.\n",
    "class PhaseSMCUpdater(SMCUpdater):\n",
    "    @staticmethod\n",
    "    def particle_mean(weights,locations):\n",
    "        locs = [complexfromangle(a) for a in locations.flatten()]\n",
    "        loc_x = []\n",
    "        loc_y = []\n",
    "        for x,y in locs:\n",
    "            loc_x.append(x)\n",
    "            loc_y.append(y)\n",
    "        avg_x = np.average(loc_x,weights=weights)\n",
    "        avg_y = np.average(loc_y,weights=weights)\n",
    "        avg_angle = np.angle(avg_x + 1.j*avg_y)\n",
    "        avg_mag = np.sqrt(avg_x**2 + avg_y**2)\n",
    "        return avg_angle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bayesien Inference Parameters\n",
    "N_PARTICLES = 5000               # number of samples for PhaseSMCUpdater\n",
    "N_EXP = 25                       # Number of times to choose a new theta\n",
    "USE_QISKIT_QASM_SIMULATOR = True # Toggle parameter if we don't want to run any circuit simulations\n",
    "SIM_N_SHOTS=1024                 # number of shots for the circuit execution\n",
    "\n",
    "# DType list for experiment metric logging\n",
    "performance_dtype = [\n",
    "    ('outcome', 'i1'),\n",
    "    ('est_mean', 'f8'), ('est_cov_mat', 'f8'),\n",
    "    ('true_err', 'f8'), ('resample_count', 'i8'),\n",
    "    ('elapsed_time', 'f8'),\n",
    "    ('like_count', 'i8'), ('sim_count', 'i8'),\n",
    "    ('bme', 'f8'),\n",
    "    ('var', 'f8'),\n",
    "    ('bme_err', 'f8')\n",
    "]\n",
    "performance = np.empty((N_EXP), dtype=performance_dtype)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 138,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set-up initial prior and QInfer model\n",
    "# prior = distributions.UniformDistribution([0, 2 * np.pi])\n",
    "prior = distributions.NormalDistribution(angle,np.pi/4)\n",
    "model = IQPEModel()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 139,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use IBMQX Rochester as the base device for the partially ideal noise-model. We \n",
    "backend_device = ibmqx_provider.get_backend('ibmq_rochester')\n",
    "partially_ideal_model = noise.device.basic_device_noise_model(\n",
    "    backend_device.properties(),gate_error=False,readout_error=False,\n",
    "    thermal_relaxation=False,standard_gates=True)\n",
    "data_backend = Aer.get_backend('qasm_simulator')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 163,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 / 25 [ 1 2.2247863572512836 2.54531063001566 ]\n",
      "1 / 25 [ 1 [2.94084272] 2.5576314486878684 ]\n",
      "2 / 25 [ 1 [2.81167142] 2.5628875896183163 ]\n",
      "3 / 25 [ 1 [2.27997105] 2.5656032435855853 ]\n",
      "4 / 25 [ 1 [3.25967239] 2.369671933681773 ]\n",
      "5 / 25 [ 1 [2.64024075] 2.392361257076481 ]\n",
      "6 / 25 [ 1 [2.51572721] 2.4106386198086036 ]\n",
      "7 / 25 [ 1 [2.01676131] 2.42560295657491 ]\n",
      "8 / 25 [ 1 [2.4287447] 2.4380449342646493 ]\n",
      "9 / 25 [ 1 [1.99341921] 2.4812820612574646 ]\n",
      "10 / 25 [ 1 [2.37077158] 2.4880674161607748 ]\n",
      "11 / 25 [ 1 [2.79005046] 2.4939565631141196 ]\n",
      "12 / 25 [ 1 [2.81977423] 2.499101227847673 ]\n",
      "13 / 25 [ 1 [2.20683687] 2.5036271556645726 ]\n",
      "14 / 25 [ 1 [2.44050937] 2.5076365743992572 ]\n",
      "15 / 25 [ 1 [2.64024075] 2.5112121219360324 ]\n",
      "16 / 25 [ 1 [2.53116409] 2.5144206473790987 ]\n",
      "17 / 25 [ 1 [2.59477881] 2.517316440536718 ]\n",
      "18 / 25 [ 1 [2.43948248] 2.5199438390105096 ]\n",
      "19 / 25 [ 1 [2.54266406] 2.5223392865542396 ]\n",
      "20 / 25 [ 1 [2.60721876] 2.524532943365652 ]\n",
      "21 / 25 [ 1 [2.62169114] 2.526549942717673 ]\n",
      "22 / 25 [ 1 [2.39060614] 2.528411372774796 ]\n",
      "23 / 25 [ 1 [2.59980335] 2.5301350462472034 ]\n",
      "24 / 25 [ 1 [2.51685617] 2.531736106498142 ]\n",
      "Inference completed\n"
     ]
    }
   ],
   "source": [
    "# Create a Bayesian Inference Updater\n",
    "updater = PhaseSMCUpdater(model, N_PARTICLES, prior,zero_weight_policy='skip')\n",
    "# Set-up initial experimental parameters\n",
    "THETA = prior.sample()[0][0]\n",
    "M = 1\n",
    "T2 = 100\n",
    "\n",
    "posterior_marginal = []\n",
    "memory_list = []\n",
    "theta_list = []\n",
    "\n",
    "##############################################\n",
    "# Run for each experiment we defined earlier #\n",
    "##############################################\n",
    "for idx_exp in range(N_EXP):\n",
    "    # Define experimental parameters\n",
    "    expparams = np.array([(M,THETA)], dtype=model.expparams_dtype)\n",
    "    \n",
    "    # Simulate IQPE circuit and get results for inference\n",
    "    circuit = iqpe_builder.construct_circuit(k=M,omega=THETA,measurement=True)\n",
    "    results = execute(circuit,simulator,shots=SIM_N_SHOTS,memory=True,noise_model=ideal_model)\n",
    "    counts = results.result().get_counts()\n",
    "    memory = results.result().get_memory()\n",
    "    memory_list.append(memory)\n",
    "    theta_list.append(THETA.flatten())\n",
    "\n",
    "    # Start by simulating and recording the data. #\n",
    "\n",
    "    # Retrieve the outcome of the simulation either from the circuit simulation or a model simulation\n",
    "    if USE_QISKIT_QASM_SIMULATOR:\n",
    "        outcomes = np.array([[int(m) for m in memory]])\n",
    "    else:\n",
    "        outcomes = model.simulate_experiment(np.array([[angle]]),expparams,repeat=SIM_N_SHOTS)\n",
    "        outcomes = outcomes.reshape((1,outcomes.shape[0]))\n",
    "        \n",
    "    model._sim_count = 0\n",
    "    model._call_count = 0\n",
    "\n",
    "    # Time the update process\n",
    "    #     tic = toc = None\n",
    "    #     tic = time.time()\n",
    "\n",
    "    # Update the posterior particles using the result of the circuit simulation\n",
    "    updater.batch_update(outcomes, expparams)\n",
    "    \n",
    "    #     performance[idx_exp]['elapsed_time'] = time.time() - tic\n",
    "\n",
    "    # Record the performance of this updater.\n",
    "    est_mean = updater.est_mean()\n",
    "    performance[idx_exp]['est_mean'] = est_mean\n",
    "    performance[idx_exp]['true_err'] = radiandiff(est_mean,angle) ** 2\n",
    "    performance[idx_exp]['est_cov_mat'] = updater.est_covariance_mtx()\n",
    "    performance[idx_exp]['resample_count'] = updater.resample_count\n",
    "    performance[idx_exp]['like_count'] = model.call_count\n",
    "    performance[idx_exp]['sim_count'] = model.sim_count\n",
    "    \n",
    "#     Log\n",
    "    print(idx_exp,'/',N_EXP,'[',M,THETA,est_mean,']')\n",
    "    \n",
    "    \n",
    "    # Re-evaluate experiment parameters\n",
    "    uniform_draw_01 = np.random.uniform(low=0.05,high=0.95)\n",
    "#     uniform_draw_01 = np.random.uniform()\n",
    "    cumsum_particles = np.cumsum(updater.particle_weights)\n",
    "    draw_index = (cumsum_particles<= uniform_draw_01).argmin()\n",
    "    THETA = updater.particle_locations[draw_index]\n",
    "#     current_variance = updater.est_covariance_mtx()[0][0]\n",
    "#     M=idx_exp % 6 + 1\n",
    "    \n",
    "    posterior_marginal.append((updater.posterior_marginal()))\n",
    "\n",
    "print('Inference completed')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Posterior Marginal for each Iteration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 212,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.5132741228718345\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAESCAYAAAAR2wXeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeXxcV334/c+5d/YZSaORZC3e5XiLCYm3JCU7kaHQQiG14/Jja3mKVWj70P54sAP9QdtXWlKnlAeelh/YYSs/aElilrAmWEkUQgKJl2xOnMSRvMjWvoxmn7nLef64I0WWJXkkjZbI5/166aVZ7j33nHs1X50599zvFVJKFEVRlEuDNtcVUBRFUWaPCvqKoiiXEBX0FUVRLiEq6CuKolxCVNBXFEW5hKigv4AIIQaEEAeFEPfnf++dhW3WCyFailBOixDi/jFe3yuEmPYUs+nUUwgRFkIcGee9Tfl9fXB6NZw9xTpm82U7yuQINWVz4RBCDAArpZTR/PO9QJ+U8u5JlFEvpWydqTpOsN0WICql3Dzq9SNAvZSyfLbrNKIOYeDh0XXLv3cQ2DG0z9+IinnM5+rvRymc6ukvbPcC2ya5zgW97Vl0rxBi+9ATIcQmoGkO61OI+jdywM8r5jGfy78fpQAq6C9sO4GDAEKIXUPDEEKIhvxr9SNeGxoOGhqu2J1fZrcQ4kj+tXD+52D+9ftHD32Ms53z1pmgvgfydR5Z/3tHLjCi7CNj1Se/zN4Ry+3L12m4nkOP8+8dGTkMNrr8iXZufh8N7cNNY7X/Ym0ffQwusg+P5Ns29Hv3UF1H1mnk8Rqn3iP3xaSP+Xj7aXRZk/jbuOBYjLVflCKRUqqfBfIDDOAE+aGfvfnXNwH3j1juCBAG9gINo8o4MuLxJuBg/nE9sC+/ngR2518PD60zwXbOW2ecureMsf2hcgfGWH7XOPWpH6rDyPaNqufQOuGLlT963fHqPUH7V0zU9tHHoIB9OFRnCWzPPz6YX++C4zXONs9rz2SP+Xj7aYyyJvO3cd6xGOtvU/0U58eFstCMNb7cyPk95iagASdY7BNCHMD54I4ei90JREb0UOvzv6Ny7PME422naYJ1Rhsa4mlljKEdIUR9vsyRw1ajyx7q4faNeDxa64j91H+R8gs1VvtvHqN+I513DMYpY2gfjqzzyP3TCkTydR7reE1GQcd8Cvup0HYNHYuL/W0qU6SGdy5dESllE69/aMecnQLskVLuyP8MncjsH2fZMbczyXWGhnjGGtrZhDNmfBgnQA4ZLjsfHCL5IYGtUsoD42zngnH4CcqfjjATtL3AYzC0D8+r8xj/3GHs4zVZEx7zIu6nMdsFBe8XZQpU0L803I/T0xqyHbhP5GdaSCn3AIdHjM0O9Y7vHblevncHY3xIJ9rORdY5Tz5o1wObpJRHR729Bbg3//rIXuxw2fk6Nkkpt0kpdxSyzQLKL9RY7f8xE7R99DEAfjlGGfeNufKFxjteFzXJYz7hfhrnXMJEfxtj1WfMv01l+tTwziVAStkknPnkQz2mPVLKqBDidiHEUGBszb/WBDwshDgspWwcdSLtILB/CtuZ7Ad2aBx5tCbgoBBiG04QuqAHLaVszZ9E3J5fpikfOApx0fInMlb7ufg/u4ZRx+DHQog1U9mHUsqjYxyvQobUJnvMJ9pPw2XhtH+obpP92xi9X97oM6TmDTVPX1lQ8rM/WqSU+/MB5R6cXul4wzyKcklRwzvKQjM8hp7vHQ6d5FQUBdXTVxaYEb37oWGDVill4wSrKMolRQV9RVGUS4ga3lEURbmEqKCvKIpyCZn3UzYrKyvlihUr5roaE0omkwSDwbmuRlEslLaodswvC6Ud8MZpy5EjR3qllFWjX5/3QX/FihUcPnx4rqsxoebmZm6++ea5rkZRLJS2qHbMLwulHfDGaYsQ4vRYr6vhHUVRlEuICvqKoiiXEBX0FUVRLiHzfkxfURRlKgzD4OzZs2QymaKWW1ZWxvHjx4ta5nT4fD6WLFmC2+0uaHkV9BVFWZDOnj1LSUkJK1asQAhRtHLj8TglJSVFK286pJT09fVx9uxZVq5cWdA6anhHUZQFKZPJUFFRUdSAP98IIaioqJjUtxkV9BWlCEzDor8jiWXYc10VZYSFHPCHTLaNanhHUaZI2pJHv/sy505EifemkRK2/sEKrn7XVO6/oiizQ/X0FWWKuk/HOf5kB6UVPja/cwX+EjeDPem5rpaiTEgFfUWZopPP9SA0wds/+iaueVc9ZVV+UrHcXFdLmSei0SibN79+m+LW1qnd233Pnj3s2LGDPXsKvQHcxNTwjqJMUetzvSwqyxH9t3/B7OuD1JUkIivmulrKGP7xpy/yUnusKGVZloWu61xeV8rfv2tDwevt2bOH+++/f1LbampqoqKigr1793L33Xdz9OhRNm3aNNkqn0f19BVlCqJdKQY6kpQ+/SMGH3iA7IkT6N1tqqevjOnuu++mqamJxsZGolHndr87duxgx44d7N+/n9bWVhobG9mxY8d53wgOHjw4HOTr6+tpamqadl1UT19RpuDkc70AVEZfYtXBX+GKRGj5wOfJSTeWaaO7VH9qPplMj/xipjJPf/fu3Rw8eJB9+/YBzj+BnTt3sn37dhobG2loaKCpqYkjR44QDr9+r/hoNEok4tztMxwO09fXN+36q6CvKFPQ+mw3JekOSres5hW7g1z3afBkAUjFcpREfHNcQ2U+a2lpoaWlhUOHDg0H+YaGhvMCPjiBvr+/H3D+AVRUVEx72yroK8okpWI5OltjrOw8whe2HuGJn/8JAI259QjPTSroKxe1efNmIpEI27dvB5yTvKMDPsC2bds4evQoDQ0NHDp0iG3btk172+o7qKJM0qkXnKGdcPJlDtXb7L1hL/sa9mGGTABSUTVtUxnb0Jj9rl27OHjw4PA4/ngaGhro6+ujsbFx+Pl0qZ6+okxS6+FOfJk+nlnRzk31Dbyz/p0AHCp3wQAkOqJwVfUc11KZa+FwmCNHjgw/P3jw4HnvD43vD9m7d++Y5Yz3+lSpoK8ok2DkLNpeGWBxz3N86zqTf+xsg2//IbgD6No5ABKd0TmupaKMTw3vKMokJAey2LbA5hyl9dW8+eVfgZmBRCc+OnAbCZJ9ybmupqKMq+g9fSFEGBgaeNoqpbzgMjIhxADQCjSN9b6izFfJTmcmxUt1ST7U3gqr3w7vvw+A4Bc3kz0TIxbV57KKijKhmejp3w5EpJQHAIQQu8ZYZoeUcrMK+MobTTof9HNhg5sTcdqv+Swvtg/Sn8wRrqjEk4uRShhzXEtFGV/Re/pSyv0jntYD+8ZYLCyEqJdSTi0ZhaLMkd7uLkDncqOLjrUf5sZvnMGWZwD4eKnGUitGJrtkbiupKBOYsRO5Qoh6oH+cwB4B+oUQ+6SUjWOsuwvYBVBdXU1zc/NMVbMoEonEvK9joRZKW2aqHe2vvgaspcyT4bZj11IbFLznMg/RjCR3pgKLGDnby6OPPlqUXO7qeExdWVkZ8Xi86OValjUj5U5HJpMpeP/O5Oyd7WMFdHj924AQIiqE2D40FDTq/f0AW7ZskTfffPMMVnP6mpubme91LNRCactMteMnza8xADznv5yMKOX+j13P8oogAF/64uPktH6kcHPdtTfg8U//46WOx9QdP358Rm5rOJ9ulzjE5/OxcePGgpadkdk7+UB+d/5xw6j3dgkhtuefTj+RhKLMokwii2blOG5G+NedV/JfgzE+39LODzr7yQbqMNxOJkeVeE0pVmrl1tbWolyJO2QmZu80AHuFEJ/Ov7Qn//pBKeU24D5gy9A/g9G9fEWZz7IpG5eZZsOGN3FnbIDWVBZNgClhnRnidvcgAImBDOHqwBzXVhn2yzug84WiFOW3TNBdUHMFvONfCl5vKqmVo9EoBw4cmPI/jLEUvacvpWySUq7Kz87ZLKVsyr++Lf87ml+maejbgKK8UViGhman+JkVJmpY3H/VKlpvfDPNV68lUVoLXufCrETnwBzXVJlPpppaORwOs3v37jHz8kyVuiJXUSbBsjwIO0VpxeV8f+taarxuANYF/YQjEfANQjafikGZPybRI7+Y9CymVp4JKugryiTY0otLDvCHS3TiZ/6VU/2/xuOuoKxsE1WBa7C9ObS0SaKnOHdpUhamQlMrzwQV9BVlEqTwI2kn3fl52tyvUh6+GsMc5PSZfQSyUQyvn5JcjFifd66rqsxjhaZWngkq946iTIKl+UGkuHLZ27jh+qfZuPE7XL31AW668VmWB3OkvGV4cjEScTV7R7nQZFMrg3MCuLW1lT179gyfD5gO1dNXlALZtsR2BbCNNP5MLb29ScJhN36/H10PsNTvpkWrRDdjZDNlc11dZY4VM7VyMdMrq6CvKAXKJp3bIdpalkcffIJMxpkCWFFRwerVqykPSxKiGkmMnOWZy6oqyrhU0FeUAiW7nGRrpp7lxhvfS1nZIgYGBjh58iSHDh2iqvY1BtxXYYg+TPzYtkTTpp+KQVGKSQV9RSnQQHsnAKbHos6oZvn69QhNcP3115PL5fjef/0tfb4Iln4ShEY6niNYpk7oKvOLOpGrKAUaaD8LQM5jYf6yk+6vPEumxTmx5vF4CAaWkQj4sNzOa6nB7JzVVVHGo3r6ilKgaFc3ECLjhYrrlpN5sY/ee17AXRvEf2UVAfcSRCCD9ORTMXREqVpWOreVVpRRVE9fUQoU73MuuMr6dfrKeqj55BbC716FcGvEHjyFflhS6etHep3bJSba++eyuooyJhX0FaVAqbgzXGMEPfz0q3fxytOPE3pLHYs+fhU1e7YS9peySO/D9GYAGOxS+XeU+UcFfUUpkJG2EbaF5Q9QumgRP//3L3DsUWfutavcR3llBRErRtrvQTfTDPYOznGNlblUrNTKQ4nZGhvHvD3JpKkxfUUpkDR0dFIECdDw0b/i0APf4KGvfZm2F5/nilvfTvmiCKFMmoSvggojQXJQznWVlby9T+/l5f6Xi1KWZVnous66yDr2XF34bb6nklp5//793HPPPYTDYRobG2lqaqKhoeHiK05A9fQVpUDS8qBbKcozgsf++xzv/OvPsPEd7+K1w7/j3n+4gxee+CWlSZtBTw3CzmGkM3NdZWWemGpq5V27dg3n5Fm1apVKw6Aos0r6EDJFOAHxviy//n4Lb//oLm74kw/z6lNP8OJ/PUQgFSShRwADK6d6+vPFZHrkFzOV2yUWI7XywYMHJ/1NYSwq6CtKoYQfZBJ3zM2mbUs4evAsLz5ezptuXMyGm27ltUeexE4HiAdKQeawVCoGZRyTTa3c2NjIvn37ipKJUwV9RSmQFH4EvciElxXHvk/v5e/mN/edIFDqYeWbKwnURaAD0pUCQQ5L+ue6yso8NZnUynv27GHPnj3U19cXZdtqTF9RCmRpAaRIo2VyRO+9l+tv8BIs9/LLr73Af/3jU2TSAdxpPzJoAjkk7rmusjLPTDa18shzAdu2bePAgenfUlz19BWlALZtY7n8aHYGEZdowSCxL+3lfV//Bq3P9PDcI2c5+bzF2nKLclcciVBB/xJXjNTKu3fvZvfu3UWtl+rpK0oBEj19IHQsPYduxJG73kfqqadIPfowa66uYccdWyhfvBiZTVEl+gEDW6igr8w/KugrSgHaThwDwHSZJPV+vrT0ON7Vq+neezd21rlSd9HSOnLZGGEjiRQ5pKZO5Crzjwr6ilKA0685c6ctt8VAsJcnun5L6q/eh3HuHGc+9GESjz1GeV0JKcugPCmRwsDWPEippm0q84sK+opSgL6OcwCYHklPZZwSTwn3uH9H7T//M0ZPN22Nf4Hx3a+SsCWRuB+pGSA0bFMFfWV+UUFfUQqQ6k8AIN0arkiW9y++hUfaHqHnrVdw2UMPUXvXXfjaXiRle/ClK5HCBCCXNeey2opyARX0FaUAMuUEb126KPf08f4zx/G7/Hz9ha8j3G7C730PkTW1oJXiSpZha87yRiI9l9VWlAuooK8ohcg5H5VAykW9bhB+7RF2Lv99Hjz1IK1RZ7zfv3oVPkvDlXZh6/mefiw1Z1VWlLGooK8oBdBMD0gb30COdXopCI0PJzL4XX52/mwnXzzyRewVSwil05DNYGsWAInB6SfIUt6YipVaubGxkcbGRvbsKU7+IHVxlqIUQLO96FaahOghF9hCpsqm8rn7uf+jB/nfL36Dbx/7Noe7PLwv+VbMXCmm28ZtQnowNtdVV4DOz3+e7PHipFY2LYt+Xce7fh01n/lMwetNJbXy0aNHh1Mw7Nixg6NHj7Jp06bJVvk8Re/pCyHCQojt+Z8LLzFzltkuhGgQQhT3UjNFmSFC+tCsFNlAJ989U8O9vB3S/Sw9c4i7briLH7z7B/RV+wmmukibFpbuzNpJxVRPX5l6auVNmzZRX19PNBqltbV12gEfZqanfzuAlHK/EGKrEGKXlHL/0JtCiO3595uEEPVCiAYpZdMM1ENRikcG0GQab1kHR+zbOdhay4cqViMOfR2ueh+ry1ezZulVWForSbse0yUhC5mY6unPB5PpkV/MbKdWHurtT/fmKUOKHvRHBnigHtg3apGtwL35x63AJuC8oC+E2AXsAqiurqa5ubnY1SyqRCIx7+tYqIXSlqK3Q/hBpigryXA2U0cmZXDXa24+U3+Yow98lVjZenwxH53BdtK2BzP/yTrTegpjGvVQx2PqysrKiMfjRS/XsqyCyo3H49i2PbzsyPVefvllXn75ZZ544gkCgQCJRIKbbroJXdcvKHv16tX88Ic/5HOf+xw/+clPuOWWWy7YViaTKXj/ztiYvhCiHuiXUo4+ezE6f2jF6HXz/zj2A2zZskXefPPNM1LHYmlubma+17FQC6UtxWyHZVm89u12kDEMbRF/dsNq7j/cxo997+bjuVfZ9Nxn4eqPYl5+HYcqf06Z5cJyOyOnQX9wWvVQx2Pqjh8/PukeeSEK7elbloWmacPL6ro+/Pjaa6+9ILVyVVXVBeWOHMOvra2ls7NzzG37fD42btxYUP1ncvbOdinlWHfyjQKRGdyuohTVS62vYGsBbC1DT2YtD8g0W9ZV0R3ewJ8+cyVc9T546musfeBvaKsSBLMGpkcAkE3l5rj2ynwy2dTK9fX1w8scOnSIXbt2TbsOM9LTF0Jsl1LenX88esz+EK/39uuBg6PXV5T55Oizv8F0rUA3s7zMelo9sDTkBU3nTOmb4d1fgms+xuLvvJueKo3S4wkstw6AkVFX5F6qipFaORwOX7DcdM3E7J0GYK8Q4ogQ4siI1w8CSCkPAPX55VAncZX5rqftBFJzYes5TlRcAZrgaTuHO9VDonKDs1D15WjLriVUJgmlejDdTlplM2fNYc0V5UIzcSK3CVg1xuvbRjy+u9jbVZSZYsacE2tSMzgZLuMt4RBPRhOUawOkA2t4rTvBZYtCUHsVl3X/hqzehfTUo1k5dXN0Zd5RV+QqykWItDMub2lZ7EovfxV6kmqPjlEfAdvkw998ml+80IGsvYq1uRzdoU404UWzc9iWmOPaK8r5VNBXlItw5c/FJl02y4LdyLZ/5EZ/J8l161n06o8o8bn4+PeO0thksTZrcLYigdf2otkGWOojpswv6i9SUSYgpcRlOB+TjlAVH6jfgsdTxXXaM0i3B3NxOT/76+u584828Ng5m0oinK0SuC0Xws4hpT7HLVCU86mgrygTiMWew2U4gfu1ilW8p7qcUGgdy7O/wRUfJLbhKly6xgd/bwXX1lfQYq/Eigg8BgiZA1ult1LmFxX0FWUCXd0/QzedmTiuFZex1OehJLSedOoEpcefJbn6chKmM0Nn47Iwv80spdKXIpiznaCPujm6Mr+ooK8o45DSoqvr52iWF4Cy+FkGBgYIhdYjZY6l5w4jXW6+cqYbS0o2LivnBXslq00Df9oEVNC/lBUrtTI4uXomuohrMtR3T0UZRzR6mEy6G016EbZJ1WAX3/zmN9m+/S0ArJav0PLKMf5f4OH+GJ9bWsML9ko+mstxKiMxpIEkNLeNUAB4/L5X6W1LFKUsy7LQdZ3KpSFuuH1NwetNJbUyOP88+vr6Jr3eeFRPX1HG0dX1U9LtlWi2F81O87d/4uRJ+e//fhQhPFRX51j8/a/ztcuX05k12PHSSdzl1SyyytCyEsghhWduG6HMC1NNrQzOP4vGxrEy2kyN6ukryhhs26Cr+5cMvLYepA/NzrCkpoaPfOQjfPnLX0bKWmpqOhHAe6rLeWtFKW8//Apm2E1bx3JkToJQQX++mEyP/GJmM7Xy0aNHCYfD1NfXF63+qqevKGMYHHwG04ySy5YDAYTMAFBeXk4kEiGTrqS62gCcK25LXTo3lJfQF9A5aqzANg2kMJCaGtNXLtTS0sLBgwfZs2fPcJBvaGi4IJf+nj17hr8FHD16lP37949V3KSonr6ijKGv7wUA0v4SQsIHMk3mtQF8l5VTXV3NQDRIba1NKGQPr3NtOMT/KXVxTK6k2jqHFB5sTfX0lQtt3rz5gtTKowM+vJ6kLRqN8tGPfrQoWTZVT19RxtDd/QKWpZP06UjND2QZOHACO2tSU1NDT7czo6e6+vXUydeUBZEhN8e1VVhaComJ1DxIW+XfURyTTa08E1RPX1HGkEqdJJMuJeO1sDU/aAbWYJbBX56iel01yWQ5QH6Ix7HY52FJwIMRWYR5LoXUnPeMnIXHpz5ql5pipFYeWdZUZv6MRfX0FWUMlnkOMxUhRxJL96O7LUJvqSP5uw7Kc35M00si4Tkv6ANcWxYiGXJhutJI4eTSz8WSc9EERRmTCvqKMopt59DdUfTUIoSRxHL5cPk0St++Aj3iw36wC6/Xy+BgyXnDOwDXhIMkS13YuomtO1fqpgfVzdGV+UMFfUUZpa/vZYSQaLlaXBmnJ+8LetA8OuW3XYbVn6UyUE4qGaGi0sSyssPrXlMWwi7zYLkNLM0J+vGB/jlph6KMRQV9RRmls+NZANriEk/WmZ3jyU/C8a4KI3wuKrRSstlFaBokk68Or7s64CVS4sF2SyyXcwI32d87uw1QlAmooK8oo/R1HAMg2t2BJ+ME/fZXj/HCI79CCIG7NkB5LkAmUwlAInF8eF0hBNeUhch4NSyXs25iYGCWW6Ao41NBX1FGSaZOYuR8uLMWHsO581UgHOZ3P7wXyzTx1IYID3pIp0tIpVz09f/mvPWvKQuSdnuwNaenH4uqMX1l/lBBX1FGsUUHZqqcYMaNy3CmWq677kZiPV28/MRjuGuDhHMBpISzbZX09TVjWZnh9a8Jh0h5SzFdzj+MVFzN3lHmDxX0FWUE0zRx+Qaw4mX4zJLhXPqrr7mGqhX1PPWj+9CrfbjRkTmL7u4lWFaS/v7Hh8u4IuQn5S3DdjtBP5NKz0lblLlVjNTKR48eZdWqVcPJ2YaStU2HumJEUUboOPEKHk+GwXgZLi2EaSWxXeCvLOXa23by0y/exanTzxMSAq/hImYvxeU6QXf3g1RVbQPApQlsbyWWuwcAI52daJPKLHj02/vpPj31fPYjWaaF7tJZtLyeW/608LQIU02tvH379gkv3Jos1dNXlBHaTjwFwGnPBlw+1/ANVLxBD6u3/h6RxUt56oF7cVX6qbBL8PmCRCJvpbfvYWz79eBueyJI3enpmxlr9huizCvTSa3c2trKnj17OHr0aFHqonr6ijJCb88LBEPQlaygriSFNuhD2Dl0l9M/uva9t/OL//g3siuzLNOq6BRxXK6rMc0f0d//JJWVtwAgfFXgcj68tqFy78y1yfTIL2Y2UysPPf70pz/Nrbfeyv333z/tNMuqp68oI2S1dqQtyGbiJEM6SD+aleH5R9sAWPN716PpLgaNHla6awCIx5ag6yG6ex4aLkcEwwh8IG2kKeakLcr8VWhq5fr6eu6//37C4TA7d+6kqalp2tsuKOgLIVZMe0uKMs8ZsTQEejEyJbi1swwEXQjpR5NpHr/3BKde6EV3ualYspSewTOEpA8za3DixEmqKm+lp+cgtu1cwesJRXDbbjTbAEv1rZTzbd68mW3btrF3797h8fqxUiuPHNJpaWkpys1UCv1r3CuEmNz3GUV5g2n99VE8/kFy8QpKcx0kPAKED0GWYNjLC81nAVi0op62jpcQCNwJ58NYVnYLphklGn0aAG+oHK+ho9k5sPW5bJYyj0w2tXI4HB4e+w+HwzQ0NEy7DoWO6Q8AR4UQTcAgIKWUn5721hVlHjl38gz+dTES/RsIxrvprTSRwo8QOdZfV8vhX5wi1pumavlKXnzsYRKhFHW5MOfsJF1dlWiaj97eR4hEriNSEsCXtdDsHEKqu2ddioqRWnloeKeYCu7pA28D7gb2ARPes0sI0SCEODjB+wNCiCNCiOLNQ1KUaYomzqHrFnUrNuJLCQyRRWp+dM3g8uvqEMBLT7RTtdz5it2R7uIyVx3l5eW8+OIJSkLricdfBKA66MGXziLsHEg1X0KZPwoK+lLKk0ADsBv44/zziZa/2NmGHVLKzVLKPYVVU1Fmnu3uAqBi1SbcaQ+WmcLW/Ohum5KIj+VvquD4Ex1ULF0OQG+mi6WBGjZs2MDJkyfx+dYQTxxHSpuakA8tl0OQA1RPX5k/Cj2R+1WgD6enPyiEuHea2w0LIYp3e3dFKQLpGwRgwLsEzfDhySYwXX48PudjsuGGxaRiOTpbs5RUVBEzevDpHtYtWY2UksHBUiwrQTrdxiKfGztngDQAdZ9cZf4o9HtnREr5w/zje4QQ072pYwToF0Lsk1I2jn5TCLEL2AVQXV1Nc3PzNDc3sxKJxLyvY6EWSlsm2w4pJdLr5Mj5+aFTLLPL8KRj2LoHQ2Zobm5G2hKXHx5/4BhaqIREn/PNoP2ZVgKBAMeODXDZanjqqfvp4Gos0mgyB1pgyvv0Uj0exVBWVkY8Hi96uZZlzUi505HJZArev4UGfSGEeCtwGNgGTCsBhJRyf77QqBBiu5TywBjv7wfYsmWLvPnmm6ezuRnX3NzMfK9joRZKWybbjmQ0ys/PJrAML75wJR5Rii/j3PxkUV31cFnB1EkO/ewk6665imdP/5icbbA6sIzkVvj1rx/hshARU1MAACAASURBVNU6y5bBFUvfwr3f/w98IockPOV9eqkej2I4fvz4pC+iKsRULs6aaT6fj40bNxa0bKFj+rcDm4F7gJX551MihNglhNief9o31XIUpZi6j7chPAlsI0i0+wwuVxBvxrmSNhApG15u/VtqATCNCGBzLPoq6ed72bD+cqTU0UQd8cRxwi6djMcADKRQY/rK/DFh0BdC3JX//S9ABXASqBx6fYL1tgNbRgR3RszmuQ+ICiEaAEb38hVlLvS+cg7dkwRZRjZ6CuFy4805V9KGIuXDy5VEfATLPBg557Vj0RewkwYlg27Ky8tJpSuJx1/CpQkMD0gMpKbG9JX542LDO/flfz8BnC200HwgHz1ksy3/OwoMze6Z/jXFilIEA+29eKrSeFwrkfFTiEAYd9a5qKq0InLespXLShjsSWPako5sC1rAReqZbhYvXkxfXwCfr5tsrpes3wPCQAoV9C9F0WiUW2+9dXiufmtr65SuqG1qahqeq793794xr9ydjAmDvpTymfzDDwB/LqWcX2cvFKVI0vE4fk+aYKgOb+ocBILoSTe2C4IVZectW7W0hDPH+ohbkhJd4L+iktTRbureVstThwMsXgyJ+ItkfD6kMLBVT3/ORX/aQq69ODezsSyTtO7CUxck/K5VBa83ldTKra2t7Nu3r6gXaBV6cVYU54rcrwoh7rrY8I6ivNEYpNB1C1doCaFUD9mQhZ5Pq+yvKD1v2aqlJUgJGdtHqUvDf1UV0rCpyIZIJJxvBfH4S+SC5diaidR0TMOc9TYp88dUUysfOHCASCRCY2Mj27ZtK0pdCp298y/5nyEqV6yysPgSACT9ywmmU6TLJbrt9NC9Jb7zFq1cGgJAauW4tS4ygTR6uZeSUxLb9gKVxBMvYYZqkJozVyEdjVFSdf4wkTJ7JtMjv5jZTK3c19fHqlWr2L179/A/junm3ym0p38rsAfYlb8aV6VPUBaMXDaH9KUAiIoavBmdQb+OZvsB8PjOT5hWUuHDG3DhctcB0HPmFIGrFmG3xKmqrCKdriIefwk9WI2tOTdQifV2z2KLlPmu0NTKFRUVbNq0CXDy8BTjRiqFBv1tUsq/GPF8YNpbVpR5ouNEO7bfGe9tTYbQs34GXC6E7UOzMmj6+R8TIQSVS0OEvPXYUnL2pWMENi4CCdXucvr7gqTTp3D7K7DyQb+/q2PW26XMX4WmVt60adNworapnggerdCgL4QQt+GkT7ht2ltVlHmk45U28DrDO6diGlo2SEx3gQygycyY61QuLaHMV0dX1uL4b5rRIh5cVX4qsiGiUeervyvgxnY5I6F9PT2z0xhlXptsauWGhgai0Sg7duygpaWF7du3j7tsoQoa05dS3i6E+BQgcFIy/MXF1lGUN4rY6Q7sSAxp6wz0R5FmAMu2EATR5Ng3Na9aWoKuuek1AtTGBml95hAVyyuIvNhHwu2M3cvBNmzdBhMGB9R1iJeaYqRWHmu56So04VqplPJfpZR/IaX8urqTlrKQZHtiCE8CKUuwBs7hskuQtgXCjybGD/oAtquOYHmEY48exLu8lHDah2UEsdMuvIkWbM3p6ScHB2etPYoykcnk05/ouaK8YVmJFC5vCl2U446dw6uVIM00UvjQ9bGnWoZrApi2QbmvjstvfCsnnzmMFQENjUjaJpf1U0oXVv4ccDaensUWKcr4LpaGYaUQ4j5gpxDiISHEr4QQvwLKJ1pPUd4opJRIaePxpNE8VfgzZ3FrfjQzjdT8uFz2mOtpmmAw20XYV8Obbm5A2javHHsC4ZZUplxk7ABhrQvyJ4GNdG42m6Uo47rYFbkngduFEJ+SUv7rLNVJUWZNajAHLoHHnSblXUHIPoPmrsSbiWG6/ATz8/fHEs10srT0TZTXLqZuzXqONR9kUe4aFgUW02X6CASTyHxP386qi7OU+aHQ4Z2DAEKIFfkrcq+awTopyqzpORsFN7g9WWKe5QSNXqRbJ5B2bqDiDYyfIXMg04lH9xHrzbDh5gb6288S6zjGItdiDMOHP5BA4CRtk4aYrSYpyoQKDfpDNzrZh5MkTY3pKwtC52tt2AFn6GXQVUcwG0X6Jf5sBqm58JUGx103mu0EoLctzpprr0eX0OlLUCr9SDuA35cA6QR7Yamgr8wPhQb9inzvXkgpHwbUVARlQeg+eRaZvzCrMxvGk02j+bP4ss5YfrA0NO66g9kuTNvgzEv9cOYM5YkUbcE0wqWh26VoAjTbKUdY+rjlKMpsKjTo78Pp7e/KPz80M9VRlNk12N6J6YsB0J70oWc1bH8OTz6Xfkl52bjr2tKiLXaMVw910ffLJioSGQZi/eg1frSs8w3BZdoI20TIQj9qykIRjUbZvHnz8PORidQKtX//frZt28a2bdvYvHkzd99997TrVejFWQ/nb2S+RwjRok7qKguBlBIrGkOrd7649sQlWjZEMiRxtztBOlQ58US1luhhVoY38spTXSxetpJXjEFS3iR6IuAsYBpodg4h1d2z5tIvf/lLOjs7i1KWZVnouk5NTQ3veMc7Cl5vKqmVd+3axa5dTl97x44dw4+no9CLs76Kc2vDu4FBIcS9096yosyxxEAWj51Dep3bRFixKDJbQtIHLtPJsOmrmDib4kCmg6oaD2c861my7e14g0G64q24ss56GWE6Qd8uNKGtshBNNbXykAMHDrB169Zp30AFCk+tHJFS/jD/+B4hxPjJIhTlDaK/PYlfkyQ9caQM4o13oeVKGfDq1Jl+cIG/avzhnSEr3G0cCtaSWrucpT1neLHlN9RWXwlA2msjpAGonv5cmkyP/GJmM7XykLvuuuu8lA7TMZmEa28VQpQKIf4Y56YqivKG1nNmALfbg9uTQtPKCSW7EKafAZcXzXLG5P2hi9/1KnL0AVx2lpePpVm64c3EejvJmU5/KuvVETKHkOruWcrrCk2tDMXLrjlksgnX/gJ4Wkp5e9FqoChzpKPlDG53CI8nA67llOQ6cdl+ErjQLGdM3huc+CMSsSyMF59j1c40J57pwV7rzNbRpAvL0jGCHoSdAeGd8fYobxybN28mEokMZ81sbW0dd+jm6NGjbN26tWjbvlgahlIhxH35Mfx9UsrbpZRfKNrWFWUO9badQ3i8eDxpUp7lhGQ3Hi3kZNiUQTQ7gz4ql77Z18fpD3yQ7n/7IgHb5qqMk3r5zX90BbYlaTr0IlmfJIgbw/BiB31ACol/DlqozDeTTa0ME/9DmIqL9fTvwZmuOfR4Z9G2rChzyLJs4r0d2NXVeDxpot5leDiOW1+JZphACN0+P0malJKOz/09qWefJXXkCP8EZITAu3Ytizatxld7hmB7LZVrV2P3JjANH3pQQ5JBarVz0k5l7hQrtfLu3buLWq+LjekLKeXD+Quy1N2ylAVjsCuNZsbIeFNomk3cXYvbjiLcHnQzi0YAXTs/rfLgj35M4uGHWfTJT7LygR/zssdDuW1T+gd/gJSSY8HfsSixnE1XNmCkYhimF7c3hxQZLF319JX54WJBv2XEYxX0lQWjrz1BQDPJeJ2bm/TkyvEZCfC40a0kUgvgGpFW2Th3jq5//mcCW7cS+fCH8K1Zw/7yMHdUVVLxkT/juZ7neCbwOAKBpi/DMFMYOS9+/yC2nsXS/RjZsXPzK8psuljQbxRCHBJCHM4/PiGEeE0IcWI2KqcoM6W/PUlAE2S9/QB0pILo2Sy4XbhzMWw9iNvrnJSVtk37pz8DQO1ddyG01z82g7qOcLn4actPSZb0EShz031GoPl1MAIs8vVi6zkQGmdaX5n9hirKKBdLrRyZrYooymzqPt1LSPeT8HYD0JVycVlWgtcikIlhuIOUBpygH/vpT0k9/TS1/3QnniWLLygrZ+V48NSDvHX5W1khqjhxuIu6Ui8YQSq9HXTqBgJoOX6MVevfPJvNVJQLqIQgyiWp53QbAXcpeJ0UDLF4FjJBNF8WfzqB4QrgL/FjZzJ0f+nL+N70Jspuu23Msh4/+zixXIx3rXoXK66owMhYWC4XwggQJImVvxFL5+mTs9Y+RRmPujZcueTkMiaJ/i700hBeXxIhgmiDPZAtA38GX9wCoREqDzHw3e9idnRQN2pYZ6Sftv6UCl8F19Zei6wU6C6NlCFwGSV4yGG7nJujJ3vUaTFl7qmevnLJ6e9IIq0BpNuH15vC7arGm+xBy5SQCVj4svlkawEvvfv2E7rpJoLXXgOALW2aTjfRl3ZOANsem8fOPsY7Vr4Dl+bC7dVZvDZMrDeNns+0aeezKlsJdSJXmXsz0tMXQjQAe6SU28Z5fztOKodNUsrp5wpVlEnoP5fEtgcw3fV4vUncvpWUpntwmSUkAgJ3zo0E5PNHsZNJqj75PwEwbZO/f/Lv+UnLTwh7w5j1JtItMW2Td61613D5K2p6aT1ajsg6Q0e2ls+pn1M3UrmURKNRbr311uG5+lNNpzB08VYkErlgbv9UzEjQl1I2CSH2jPVePuAPLVMvhGiQUjbNRD0UZSwdLVECWo6klsPrTWL5LyNsduOyAkR9UGb6MF2Qe+JRam97L741a8haWXY/tptH2h7hg5d/kOd6nuP5m54HE1aVrWJ9ZL1TuG2x/Ny/IbQdaPmrde38jVQ0y4207XGHiZSZ8+qrdxJPHC9KWZZlousuSkLrWbPmswWvN5XUykP59Hft2sX+/fs5cODAcOqGqZqLv76twFDu0FZg0xzUQblESSk582IfAQ2S+iBud46kZzkB2YMmvQy4POimcyGVOxuj6mMfI22m+cuH/5JH2h7hjqvvYPfW3Xzn979D4CknP8/OdTsRIt+Lf+b/UBp9khK/gIzTp0rqTvDXLA+Z48UJPMoby1RTK2/ZsoUjR44QjUY5cuQImzZNP1zOxYnc0UkkKkYvIITYRf4uXdXV1TQ3N89CtaYukUjM+zoWaqG0Zbx2ZKKSZDRBnRYk6usF4Ni5HD7Rjy4WM6B5EPkMm3LJIh5/9RW++eQ3eT71PB+o+ACLuxbT3OWUm30yi/Y7jZr6Gpo7m9HNFNc89fekS9dDthxSA5imi2jERV1/Gs32cuw/v0PynYWn+V3ox2MmlZWVEY8790qorf0bipUIY+gmKsBw+WOJx+PYtk08HudjH/sYDz74IF/4gpO67M477+Td734373nPe/jEJz7Btddey69+9St+/etfEw6Hh8tdvXo1hmFw2223ceWVV1JVVTXmNjOZTMH7dy6CfhSYcP6/lHI/sB9gy5Yt8uabb56Fak1dc3Mz872OhVoobRmvHc82neGE1UzQFabXfwaA4OIrydnfQHP7SdkaIp9hs37ne/nPkmd57sxzfGrLp/jQhg+dV9ZQEqxbbrnFeeGRfwIjiufDP2DJr/wkDx/HNL0YkSCanUbYXqraXmXrJPbvQj8eM+n48eOTzntfiELz6VuWhaZpw8vquj78+OzZs5w9e5Zjx45RVVVFKBTibW97G0uXLj2vjP379/POd76T7du3s3//fh566KExh3d8Ph8bN24sqP5zEfQP8Xpvvx44OMGyilJUbS/14y9JEjRLsX3O1+z2VICwmUJ4PAjTRpNBNCtN8+oM33rxW+xcu5MPXv7BiQtO9cNvvwIbboMlWwhFWoiZJhhegiUuNDuFwE/6ueewolH0ImZNVN54Ck2tPDQUNPrxdMzImH7+ZO2WoZO2+dcOAkgpDwD1+Rk+qJO4ymwxcxbnTkTxBgbRvaV4vUlA0JHQcGWyCLcL3bRAhtDsFHe+/GVuWHwDd1x9x+tj9uN55rtgpOAGZ6ZPSbkXw5KYho+qQBxkGvCDbZN44okZb6syP002tfLQHbd27NjBoUOHinKP3JmavXMAODDqtW0jHqtpmsqs63htEMuwyaXOgn81Xl8SXa+kM5HjzVkL4dJx5dJoIoyQKYLuIHtv3ItLu8jHxLbg0D2w7C1QcwUAwXIfOSnQDD91wTNk5QakiKCHwwz+8EeEbrwRfQaGHpT5o1iplSc74+di1Nwx5ZJx5ng/QjdJdHeQERZebxKvt5ZUNAqZALpLoNtJpAhgiATvX/9+SjwFBOZXH4LoGbjm9V5YKOwlJzUwAkT0dhAZbC1AbvvbSD7xBC3b3kbfN7+FrTJvKrNMBX3lktH2Uj/lVQmCWgkJkcHrTRIMLsUV74ZsGM1t4zFi2K4gpivF+9e//8JCbAtOPQG51OuvPb0PSupg3R8OvxSKeMlKN+SC+FwppJbB1v3cv0Wy4sABfBs20H333bR9dPpf1xVlMlTQVy4JycEsfecSeAM9BN1hEiKN15vC71tCMNmDlisDf45gehDDHcQX8FLmLXu9AMtwxu2/cjV8+51wzy2sDKZZFshAazNs/Qjo7uHFfUE3tu5HM5zpn5aWw3L5eeXY8+jrV7PsG1+n8uMfJ/X00xidnbO8N5RLmQr6yiWh7biTN99InWNRxQrS3n40zcbwLCaS7UXPlYAvhz+dxHAFqKmqcVa0bXjuXvj/NsIDfwkuP2y7E1L9fG3TK9yx7jToHtj0p+dtTwiBv6wckXUu9LJ1A4C6gSCPtD1CV7KL9mudS/ITzY/Nzk5QFFSWTeUS0f5qFE9Ap7ftBJev2MmJwG8AiOt1hK3f4bEC5Pw5vDEnw2Z5RTmc/i089BloPwq1V8IffBFWbwMh4M07ef7vNrM1Eoc3vQ9CVRdsM1QeQWR9ANi6iZAQyfj51GOfchaQkv9drhNsfpTyP1G3n1Zmhwr6yiWh+3ScSHWOU+filMgwts/p+cfEIgKyF01WkQhIfFnnSstSsx2+9X87Y/Xv3QdX3A4jc+aUVLP7+VXcWBXlH3d/fsxtllVVkO7yAmDpJi4TfDkXf37Fn1MTqKE/28+h+n+n4rdPYmcyaD7fzO4ERUEFfeUSYOQs+juSLFndgy7c5NImbk8CgH5ZhkcMAHXEAxY+w8mwGWq9Dy6/Fj74Q/AExyxXInispxwCF15g/pPn2rH9AWTK+YhZLhuXCboh+cSmTzj1sg3++vL/5B1HBkk99RShm26akfYrykhqTF9Z8PrOJpC2xMicY1HZMhI4M3eE8NOW1vAwiO4OMugBl+H0tv1BA3Z8a9yAP5HDp/r5xPefobktip12Tu7aXucm65r5+kfOrblZ+9bbyLihp+nBIrRUmU+i0SibN28efj4ykdpkNDY20tjYyJ49YyYunjTV01cWvO7TMQAGu1tZW7uZRC6D15vC46mhZSDFWiOF5vUyoNkstfzgBv8f3AGldZPeVsaw2POD55ESXowmeYuhYZputGoBg+AyRnzkuo/zRxXreXyl4MpHH0ZKefErf5Up+eyJsxxLpItSlmVa6C6dN4X83Ll6ScHrTSW18oEDB9i8eTO7du1iz549HD16dNqZNlVPX1nwuk/H8YVsBtrbqAwuJeHL4fUl8fvrOB1No+ey6C4XUXzDGTb9l18/pW39xyOv0dKT5JPb1jCIjW2BYXipWuZ81HTTA9mEk5ztq9dRf38jPVcuxtsbJ/PKK0VrszK/TDW18qFDh9iyZQsAW7du5fDhw9Oui+rpKwte9+k4obIBom0QMENE/eeI+FIEAkvoiKXRMibC7SJni+EMm97A5D8aL7XH+NpjLfzxpiX85S2Xcd8TpzEGBIbhoyJ4hpi00WwPfPlKSPXChvfCyz9nzQpnSOm1X3yfK9b9QzGbruRNpkd+MYVm2RxpKIfOUOqFu+++m507d7J9+3YaGxtpaGigqamJI0eOnJd4bSjQb9q0iUOHDlFRcUEm+klTQV9Z0HIZk4HOJIsWd+HWvIhBm75IH9XuNB5PHbnBAex0EM0lEDmcDJt2Gk2f3Jdg25bc8cPnCQfcfPYP16Npgus3VJE55cKXKSHo70a3nJz6aC543/dh7TvgV/+LW3/7FZpr6yh5uAn+5z/MzI5Q5pWWlhZaWlo4dOjQcJBvaGi4INPm9u3b2bZtGy0tLbS2trJq1appb1sFfWVB621LgIRcuoOlizdgYpGSzhWwcddSKpI9kA1DiUTL2ggZQpeZSW/nJ8+18/zZQb608yrCAQ8ADZfX8Owv3JRmytD0M2hWGiF98LfHXr9698ZPEXz2v0itdLPit31kY1G8pSrt8kJXaGpleD1R244dO7j99tunvW01pq8saEMncWO9bSyuWEO/SOD1OtM1e0QNi9J9wykYXEYKIQK4xOSSoGVNiy/86hU21JXy7itfP/l73WWVxHQfIl2O0Cw0mQbpPy9dA74yuOXvqC3pRZPwwqMHxtiCslBMNrVyNBodHvvfuXPnuP8YJkP19JUFrft0HH+JycCZfso91ZwORPF4nWRp3TJCVaYPV64UAmk8RhwpluFyWReUY1g23fEsHdE0/ckcb7mscvi97/3uDGcH0tx12xVo2uuzb3xuHXdpOXoq/0GVKRB+MokEvlDo9cI3fZirfrOPUyJOxxNN8Ed/PjM7Q5lVxUitHA6Hi55aWQV9ZUHrORMnFI4zcAZ82QADgQ5KQmlAo93wEbG7cFsBTH+GkmwvlitIwGfS9FIX+3/dSlc8w0AyRyxjnlfu8ooAuWA17swA//7ICa6/rJIbVl+YiqFmySK0WGn+WQYpqjh+6Ek23vK21xfSXYRu/X/o/vHn0J97eeZ2hqKggr6ygGXTJtGuFLUr+3AJDyJm01cVZ0l5gkCgnpZBg1V6J7pcRjxkUd7Zi+kO0pnp447vHGZlZZArl4QpD7gpD3qoKfVRU+bDtCSfe+AYHRvejzfRTjZlsOf3141Zh8vXLKHrEWd2jhQZpObnpRdfOD/og5OWufZ/seSFLOf6TrG4YsUM7x3lUqWCvrJg9ZyJA2Bku1i8aC02kt7UAJf5eykpuZ7WjgRXunpBrGYwIAlmYhiuAL1mD5955zr+7LqVuMeZxbNlRTlv+ZuvkIqs4V1X1nHFkrIxl6upK6EzLTEMD1LLYukBumK92Lbksw8c43Rfiv/4HxsJBwIsvWI95rOv8swj32Pxjr+bsf1yKbkULniTUk5qeXUiV1mwhk7iJnrbWFy1lkGRApFEiAFCofW0D6TxigFcbj8DHg2XaYLQuOXK5ey6cdW4AR8gHPBQ9eoDVB+/j39+75vGXS5U7sXMGWQyIWwth617SVsJ/u7HL/C9p87w29Y+/mT/7+iJZ1n5ro8D0Pv4Q8XdEZcon89HX1/fpIPiG4mUkr6+PnyTSNanevrKgtV9MkaoXKP3ZAeRJTV0kSYYHADAF1xHdDCLMFIIj49eoVObc6ZaLqotL6h8AfgHT1Pqc4+7TKjcR8oSaJkSbFcOAGkY/PfTbfzVLZdxTX2Ej37nMDv3/5bv/V/XM1gJgRO95KwcHt0zvR1wiVuyZAlnz56lp6enqOVmMplJBdmZ5vP5WLKk8IvPVNBXFiQpJedORKlcnKK3VeI3QwwEBygpGQQg5l5DZep3JAxBwOWmX/rQsz5sAYFFhQX9QvhDbtJ4CGbC4E1DFnTT5CPXreSTb1uDEILvfOQaPvLtQ7z/64e4e81iLjt8jiMnfsHvrXtP0epxKXK73axcubLo5TY3N7Nx48ailztb1PCOsiBlByGTMPD4oujChRaX9IkYFRVp3O4KOq0Qy1PtJDIuhFsnZ7jQsvm8O9XFC/pCE7hKKhDpMO6gM9yk5wSf/cP1w2PNV6+M8B//YyOtvUk6LmvAn4PnD36jaHVQlJFU0FcWpGS389vKdbIovBwpJT3pAYKhfkpC63lpsI9lmQ5y2RDCpeHLmM5FWkCg1FvUupRV1yGSJfhLnW8Z7hwXnFy8ee0ibllbxb8OOrdQHHyp1blVo6IUmQr6yoKU7JaURHwMdJxmcfV6UmTJGCl0vYtQyXqeOPEkS8xOPIkweEyCmRRu07mIKlhe3KBfs3I5JIL4y5zzCbo59sfuM+9czzk9RKzcR81ZyYtPfbmo9VAUUEFfWYCkLUl1Q93qEvraTlMZXEyfN4XfPwiYlITW09qeJOxqpzxZggim8RsDuO1ydJnG7dGLWp/KZWXYMT++iHMlsMtw05eIX7Dc6uoS3nf1Un5Xdjnr2iQ/ePbraNbkUkIoysWooK8sOP0dSawclFaksUyTkCwjWpIhGHJ62jJpETcjuD29lKRD2P40PrMHoYXx6pNPtnYxFXUhzJQbLWgAoFseHv7Nz8Zc9m8a1vBi7RWUZOBEjyRy9sdFr49yaVNBX1lwzr3qBHdN9CIQ6AmNPj1BZUUaTXh4vOlOUpT+/+3deXxU1d348c+5d/Ytk2SyAwkQCLLJKiIUqEJdQatotdXq0/aBWmvt9tMuv7Z2edqnti9a29qqLa22WquFWrdaBQQBRfZVCAQSyL5nkpkks917nj8ShCIoS0KY5LxfL14JuffMnHPnznfOnHvu92BaW7EbHjo9Bq54I3GrH4er5+uTmuMiFrMQM6wI00Az7NRt2HjSfQMeO9M/tYCIbmHifni36VVob+r5SikDlgr6Sr9TtT+I1Q1tTRWkuLIgLmmItZDiC+LoSPDXTknctNApIyQsXtrcJt6OFqL2VHz+nh3PB7A5LJguL2bMh2Z0oEkHzvr3J3U76q4rLqJp9GQuLRY877Yj33x/Ii5FOVsq6Cv9ijQlVSUtuDOhqngvBYPGESFOW2cIq1bO4YhJWmMmee21tHdYiDi8tDg07J0x4jYv/szeyWXvycnB7ExBMzsR0onF6TzlvkIIpt55M6ntJokGK3t2PIlsONAr9VIGnl4J+kKIhUKIuUKI+0+xvUUIsVUIobowSo9qqm4n2p7AkRql/nAp2WmFNFnCWG2dYE2wXlopqspjSGcNsbCOtLto0JxYO7vGdVJye26O/vGyh+dDhw8hOwEnjZnOD0wP4JkzB+x2LtsHz3q9HPjLlzDN/ptOQDl/ejzoCyEWAkgpVwJBIcTck+x2s5RyspTygZ5+fmVgOzqeL80qkBKPSKVCNOLpTr+QWuXCsPnIj9Wgr5mTGQAAHjhJREFUhWzoVo0m6UKPdvXwfUMCp3zsc5E1NAvZ7gatA/DQ4JYcbG495f66x43nIzOZWaLzmsdNXugdnnzycRX4lXPWGz39qcDR5dxLgUkn2ccvhBjWC8+tDHDVB4L4Ag46myvQrTba6xN0UInX23W31vWbszhYOIossxpfKIBVixGLW9HjXT18b6bngx7+rKXneTBDXjRrE3FbgESikfXbV39gGd+VV+FtizGk2uCptMHMKvsF3/r71n6dQEzpfb2Re+fEQdGTLd+eBjQLIR6TUi4+caMQYhGwCCArK4s1a9b0eCV7UjgcvuDreLqSuS1SSo7sk3hzoPrQEcL2IfhMQYs1QqqvgXDUQu6Ww+xYMIpL9i/D3RlA90ns0fh7QX/b7k3oxaeXijcYDAKc1vEyDUki5MPmrSYSduNtaaLy7ddYYzn1cJKw6AQsFq4q8fLM5ZLPNe3Fu2spSwgyOSu50mYl83l1omRvS2+cOUG6gvopSSkfBxBCBIUQC6WUy06y/XGAKVOmyDlz5vRCNXvOmjVruNDreLqSuS0tte3sfXYje4HUtkbSB00kRoIWYWOIpwV/rZeODDe1wolubcYVH4GwWXFG29FNP7rZyRXzrj3t5zu6XunpHq9nXi7GkbOdtjD4ggLd7PzQstv/9ATTD5bz6zkhVg2fzpcPPc9na67jq7dckVR54pP5vDpRsrelN4Z3NnOstz8M+I+FIYUQi46O+wNqArLSY6pLunree6v2I4Ab0qM0izAWSwSfLYZra4jKO+5CdBok9BCm8GFxadjjTWj4sWodvVo/R3oGloyupGuudhe+6IfPFIpOmYzeGGRWQ4CnfG7cdHJRw6u8UVzfq3VV+q8eD/rdvfZhRy/gdl/QRQhxNPg/x3EXeE/s5SvK2dqytZZ2Ibl6cARNEzhDggZ7E15vV9/CWmnhwGWzcIdCRCMG7Q4Hwh3HGW/EtKTicJx67nxPSB8+CJwmSBNrLBXd5yJqfPBzRiZMQPN6ufVgBjuCB9iRN57bbW/yq5UH1Ni+clZ6ZcqmlPIhKeXKo8M43X+b1/0z2L1tpZTyod54fmXgSRgmtYdaaXULHLV7yXO0ERMTKXM04/TVAZA24nK2G1AYrKaz3ULEqtHubsfX0UjU7sft691FS/JGFWCGM7AkgmhmBg1ZNrbWtnxwIZsN33XXEninhAzDzdMZWQyX5YjqbawtaezV+ir9k7o5S+kX/r7uCJ44jBtlp6GqmoKUDKRppzreSJq3GUstpC64la3NYbI7qpBtVmx2Gw3uKO6OdmI2LykZ3l6tY8bgFGRDLprZCCJAvRu27X3nQ8v5b7wJGY1yd8NYVrQdpMbh5TPOtTysevvKWVBBX0l6kbjBCyu7ZgkXNj8LCGy5C0hgoBkWUlzN2Brc1E6YSKipk8HU4Qin49EFLe4Y9k4bCI203JNNNOs53jQH8ZoAmrWRuC2DeKKRpuIPnrYJ4Bg7BntREePfacCQBv8eOomrxdsUl9fyTmlzr9ZZ6X9U0FeS3tMby3G1GWgWQVvZW1h0gc3MY4+jFIctguaI4g9MZVsogmiJMdRagr8jA7vFQtBmxRrp6uGnDO6dG7OOEprAcGRjcTUQt3nxNDfhOlz54eWEwH/TjZj7DjA3VsjLehyr0cEnnJtZur6sV+us9D8q6CtJLRSJ88jqgxTpVnJTm6lo9zGo8CJcrTp7beV4ui/iZky7ma1t7dha4+Q5ynDH0tHsVlqEG717xSxv3gfONO4RQyZORPd0LZCeEgRvq/20hmh88+cjrFZu2O/jQHslBzJH8N+e9awqrqOssb23q630IyroK0nt8bWlhMMxPBFJeudKGqNuAnljsUgdi9eF31UPJqQOnc2W1jCDWqto18NIfBgOK6G4Fy3RfTdumqPX6zvqsuFY9a7ncbc7cZHFkY4PXyjFkpqK54oryFpXjN3QeGXQKHJDuxmlVfOnt1RvXzl9KugrSau+LcIf1pXx8cEBkCAjxQA0hNqIiThGKEGqvRYnuXRKG/trw0w29lGesNJpc5BICWJ2WNEMP5oZwebs/btcM4Z4MSKZAFhjqWipNlYdrDitsv6bbkS2tnFbSxGvRKoxNSvfyNrE37dU0toR781qK/2ICvpK0np4VQlxw+TyNB8Cg1bNhcXtJCOURpm7gfZYFGdaK/6c6ewIdUAwxiXaPo5EXJgOO1FvDdb2dnSRipXweamzEILUQROxxFvQjAB1mRZK9755WmXd06ejp6cze6+grrOBrSNmMaN9BUY8wl83lfdyzZX+QgV9JSkdagjzt80VfHLaENoPVhKwHKIi7KIloDE8OphmLYjDEUY4DHz+CWxt7UBriTLWdoBoWMdus9Hoj5LSWYXUUrHZYuet7mPnTu2etplBvVugHX7jtMoJiwXftdfg3rSPQMLJy/50LNEW7snZz5NvHyZumL1cc6U/UEFfSTqGKfnRy3txWDQWTcuntloQcOyjo70TzdV1UbautZwsrWvYxOcdx9qWEDnBRhoczeQ3puC0atR5HAzuOETU7sfpOX8JzPJGpiK0FuK2DGKJJtKqWzBPc759yvwFEI/z6bpCXm/ZQzRlCHfY3qS2LcK1v1rHSzurMVT6ZeUDqKCvJBXDlHztuR2s3t/A168somPrNiQaojslcqH9IqLWOE12nZzBBkLYEM4RbKprY2L0XTY5HBTWZ6DbLJSLbDIjDcRsPnzpvbA47ilouoZujRK3+fA0N+JuTWN/qPO0yjrGjsE2dChTdnYQjod5Y+RM0ureYun1GZgS7n1mO1f9ci1vH1J36yonp4K+kjQShslXn9vBP3dU8/+uLOK/Zgzl8Nt7ceptVLZ3EvTEuSR6MRXRYqSm4cprIiVlIhtaoxgtUaZpxewznLii6UinoKIjBxnTQGhk5Gae17ZkDR0BQErQxC2GsubA4dMqJ4QgZcF8LDv3MyaRxXJCIDSu6Hyd1748i1/fNpG4YfLJ32/kBy/tJRLv3XxCSvJRQV9JCgnD5CvP7eSFHdXcf1UR93y0EOPwRsqb8xg8JEpdSQm2jBwcEStHQofwOVuIJyrIyryGN5pD2FpjjLLuw1ul0ZiSg5HShqslQWek6y7clCEZ57U946+cAYAr7EJLt1O+b8WHlDjGd911AHy6qoCNjTspHz4Ltj+NLg3mX5zLq/fN4s7p+fzxrTKu+/V6NY9f+Q8q6CsXvFjC5N5ntvPSzmq+cfUovjCnEIDaV54hKj1Y8p0IQ3JJ+iwMGac61U5efhMgCATmsaqxlfT6ZlpdjYwvg1ZfgA5/M8PaDuNq7Qr6voLz29PPLMoFwB7z0zykE0fZutMuaxs8GOekSRRurEZHY3nmIAhVw5s/BSlx2nS+f/1Y/vyZS2gMR/nyszvUOL/yHhX0lQtaJG5w91NbeXVPLd+5bjSfnz28a0PNTg4fMtE0k92Vb2AIyUXtw6hv2UnUZsOfVoY/ZQrVpp/Khg7GRfayxW5jeGMBdj1Otd8kW+4htSMHpElKhvO8tsvmtGKJt6IZ6bzr85BaEyVxBoE5ZcF8jENl3Cgn8s+WPcTH3wprH4IV34Hui8KzRmbwg+vHsrMiyB9Vugalmwr6ygXFNCUv76rm4ZUl/P9/7ubmRzewqrieH94wls/OHNq1k5Sw+iccjk0ja7iH5n378ebko0V0ykU1LlcbFkstmZlXsbq5Da2+k+naXirCTjpdo3HRSrkjBVeigrhzBCnuKDbH+V9+0CqCCDIxO1rwtOeyu/X0h2F811yDcLmYv12jOdLMmok3wdTPwdu/hle+BmbX9M3543OYNzqLn7++Xw3zKIAK+soFpLk9xn89sZkv/nU7v1h5gFd21RBLmCy55WLuuDT/2I7rfk5w7w6CiRyq7VvxhnUuzZyNaSaoGORj+PBWADIyr2J1cwhnQ4SZzs1kVGjUpw3B7nRSEcvC1hSlzVfAoBGnXqe2N6Vk2Yg488hsKMbpzGXN3v2nXVb3+fDfcD3ON7ZQKDNYdnA5XPNzmHEfbFkKzy8GI44Qgh/dMBabReOB5bsw1TDPgKeCvnJB2Hy4mWseXseGQ0388PoxHPjR1Wz/7sd47SuzuHHSoGM77lkOb/yIw+mfA6Ck7F+YFkGgxUd1+7sEDYP0wBGgECyZrK8OMj70LmW2EGPK/YQtIdy+ALZWE0vLcEzNypBpw/qkzaNmjMGwOMmu1YgMilNb/PoZlU+9/XZkPM7nSvPZUL2B0rYymPt9uPw7sPs5ePZ2iHeS5XPwnWtHs6msmd+9eaiXWqMkCxX0lT63oyLIbY+/g92q8Y8vXMYd0wuwWU5yapZvhOfvhiGXUZqYg5kSJqvWZNSIqQjNw8HMCD5fBMMoQ4jJbGptJ17bwXx9A2ttHtIjo/HQRIe/gdxQOd6OQpCS3KLezaN/KsNmjQZp4gwNpzjbjqds4xmVtw8bhnvmTIa/UYJf93DvqntpijTDrK/DtUvgwGvwlxsh0srNUwZx3fgcfvbafh5TgX9AU0Ff6VOxhMn9y3YS8Nh58Z6ZjM1LOfmODQfgb7dBSh7lE35HzaE2DoqXsCd0RrRlEJFRyrRORo+JACCYzIqmVux17cy2b6KuyU5D2nAGO7MoS60hLb4XoRfisoZwuK3nscXHOD027GYtmhhFKBYktV6nLnpmidPSPn0HsrGJ34hPUd9Rz90r7yYcC8PUz8LCpVC5CV74IgL4xScmcO34HH7yajGPqsA/YKmgr/Sp3645yIG6MD++cSwprlME32AF/OUGEDrGrX9n/Yu1mL4IergMl8eHK57NIfsRDCOO27WNlJRJhAjwdHkjl7VuZ607weQSnWaXxJKSS2O6ILW5hjbfMPIK3Oe3wSfIyJW0u4eSXV1CSjyfx3YfOaPy7pkzsRUU4H1xLUvmLKGkpYT7Vt9H1IjC2Jvgiu/Cvhdh+1NYdY2HPzGB+Rfn8r+vFvPt53dT3xbppZYpFyoV9JU+s782xCOrD3L9hFwuH5V18p3CDV0BPxqGO/7Bnt12Wmo7WJ3+BLmNDsYbQxDudA6khCgqChKLV5E/ZDGvYCda28F87S1etXoZXzWSFBGm0ltBGbk4mrMxdTtDpxWe30afYMzcKUhNJ6POjQh4qNnw3BmVF5pG6u23E9m5iwkH4vxgxg/YVLuJL6z8AqFYCKbfC0NnwasPQNMhLLrGL265mLsuK+BvmyuY9bPV/Phf+2gKf3hOf6V/UEFf6RMJw+T+5bvwOqx897rR79/BNKH4FXjyOmitgk8+S6e7iLdfPECFv5gRmo4mBXmeKTS4gzSFmsnO2YbHcxH4ZvNvaSOvJkyuewdFuzVqAxeT7Q4QGbyV5qZM7O0jAci7OO88t/w/FVw6EmHGcHQOpzRfZ1jNC+w6g6mbAP4bP45j9GiqvvJVLm/M4sczf8y2um3c9e+7qI80wg2Pgm6F5Z8DI45F13hwwRhWfXU214zN4Q/rSpn10GqWrDhAKKLy8vd3Kugr551pSu5fvoudFUEeXDCGdI/9uI0GbH0SHpkKf/skxDvgtr9C/nReeGYDiahJy5hdDDtkYZS1EIsni4O5YXJyqjDNGoYW3MujlQ0kajqZ0Pw2r7mszNmXQ73PTdTvoThrMEMbtmFlBA6acflsfXcgAItVxy4rkPooGmnFd3goT647swu6msvF4KV/wDZkMBV3383c0GAeueIRKkOV3P6v2zlMHOY/DNXb4MV7wegK7AUBN0s+MYHXvzKL2UUZ/GpVCR95aDVPbzxyWks4KslJBX3lvDJNybee380/tlXxtXkjWXBx7rGNTYfgT9fAS18CmwcW/hHu3Q7DL+f511bRtM2gsmAXlx9Jo6OlhVH+6dT7QuytKGZ44T487iJEykdZWlaPZ38zi9yvEi53UJ/5MdJFDQ25m/l3/VwuDa8j7B1GIEv03YE4TlaBjU5XLplVR0hxjMex43ennWr5KEtqKoOXLsWSEaB80WIm1Nj401V/ImpEWbRiEbUFl8Kcb8HOZ+DphRBpfa9sYaaX335qMi99cSajc3x8+/k9fPMfu4klVH7+/kgFfeW8SRgmD770Ln/bXMG9lxdy7xVdmSaREjb8Fn53GTTsg48/BovWdF2I1C08te7vHHkhRltqHXMznVTu3sW0ziw0TybrbAfIz29EiDoKhn6RJUfqMUqCLDb+yX5HJXN2D6LWEyE7bRT7hzhYUPY6tcFsDIuT4VNH9OXheM+E62YBkFGXQmR0O4NLI7xaeeapka2ZmeT/8Y/oqX6O3HkXmcvf4tHLf0tbrI3FKxYTvHQx3PA7OLwell4JzaX/UX7coBSe+uw0vvjRwq4Fan7/Dg0hNdbf36igr/SquGHy5oEGHli2i6n/s5I/bzjC4lnD+Oq8kd07RLrGml/7JgybA1/YCBffCkLQGm3lZ+uXULFMgs3kY5c42PXvVxjaLsifdAc7vRW0x6oZOmwbbvcIVsSn8sT+GqZU7GCe+yVWhNIIp3yUQrugPLuE9fXTydZfYUTVVADyZ4zss+NyvLxxBeiJMLZoIa+lefBrl7Dq5T+c1WNZ8/IYunw53o/No2HJEtzfephfT/oxlaFK7ll1Dx1jboDbl0NbNfzmkq6UDW0175XXNMHXryzi17dNZE91K1f+ci1PbzyiErb1IyroKz0uljBZubeOrz23k8k/XMGdf9zEK7trmD0yg6V3TuEbV49CCAHtjfDn62HPMrjie3Db38CXQ0lLCd/f8H2uf/omgi+68cbTmJh1gHV//QM5+Jk29X6aZITdcj9Tpm4A2mnJe4gHdpeTt7ua7zh/y/cdGdy2Np+ou51IdiplRY3Mr36GwJbJ1ObOZeR4L95UR18fKgCEJvA5G4k5xnLp7lcpuSRC/q51Z3xB9yjd4yFvyRKyH/weHRs34v/8D/lFxufZ07SHhS8t5C27Fb6wASbdAVufgF9NgJe/CjW73nuM+Rfn8s97ZlCY6eHbz+/h2l+t443iOhX8+4Hzn2VK6bf214Z4dnMF/9xRRXN7DJ/DwtzRWVw9NoePjAjgsOpdO0oJ+1/tmkYYroObnyBx0XxWl6/kmeJn2F61k8k1H2Nh9QNoUiMQ+jtbdlRS4BzCtIJPUSNaWe/bx5jh69H1OuzDl/K5nXHyNpfwY37DYw7BohcKOJQ/l4IsO9snv8P6stnccOQJwoHF+L0xPvrfk/r2YJ1g3n0LWPbTTfgaPkOT9xUC1in8/ncPMHPqTWf1eEIIUm+9Fce4cVTd92Uyvv4rnlj8Cb7reYfPr/w88/Ln8aWP3EfBZV+CtT+H7U915ezJnQiT74JxtzAq28eziy7l33tq+Z9/7eMzT2wh2+fghol5XD8hl6IsL5p2YVwXUU6fCvrKOSltCPOv3TW8sruWfTVtWHXBvNFZLJw8iJmFGf+ZTiERhdI1sOYnUL0dUofScNvT/CN8iGX/uJqOpjhTWuby3zW3YEYs+CObSDS+SbM/wBX5t+DW81hrK6bKvovCYXtx+OpYH3iEJzbamLvjZW7wPMFLHXauenM2h4cOoiBV582p+2goHsz8pl/iiX2DuFMy/xuzsRz9ALpAZBQGmH1LHquX2Rh2cDpbph3hks0udkSWcvllM8hwnt0sI+eYMQxdvozqb34LHnmaX44qYu9lH+FnsTeZf2QFo9NHc82Ia5h92WLyS99GbPszvHQfrPgeTPo0YuIdXD12BJdflMnKvfUs31bJ79eV8uibh/C7rEzJT+PSYWl8ZEQGI7M8Xd/glAua6I2pWUKIhUAQmCSlfOhMtx9vypQpcsuWLWdchz/c+TMQWQgziCnaMPUOTC2KqccxLAkSFhPTIonbBUK3oGkmvqwcMm2ppDpScNp8GHFJPBKnM9RJZzhIItqBaSaQho6ZsAE60gIdRgzNY0fYdXSHDd3pwJ2WijcjgCcjE4vdjSksmAgaWkPU1tbR2FBLZzRCAh0TnbjmwLA4Ma1uhGbBkPK9aXMWTcOqa1h1gaYJNAF2i06G106m106G106K04rPYcXrsOCw6tgsGhZNnPGbcM2aNcyZM+ek2xKGSU1rhN1Vrbx9qJG3DzZR2p2ud3J+KteNz+H68dmk0QptVdBcBi2Hoekg1O6GhmIi0qA4bTDvjrqCTUJSvqOMEbXDGdIxDhlzYibqsHfsx2GJ4nVmk+bKw0xJpcZbQnvKQfxZh2nz2dmXuJgt1dOZUL2WwsRG4nUagdqxWN0FaEMSvDssTqVmYdzud0irHkfMNZO41cX8RUUMmTzkjM+nkzl6nNasWdMjjwfw+i+eo2R/AFdoHdWDSjACg3FXlNE58TruuflGUj3vX8tXmpKOUIxQYyctdUGCdS3EInESUYNE3MDmsOPwuhGlJSQ2b4DDB7CKGJ35Xvb7g2z2V1EZMDHSfQzPHccYm5+imr0Ulb7DoHgM3ZMNBTNh0FRIL6TJnsfqWjubykNsKmvmcFMHAJleOzMKA4zJ9VGU7WVklpeAx46uiQ88r86FlJLOuEFHzCCaMDEMScI0iSZM2qMJQtEE0biJVRdYdA1dCBKmiWFK4t37JgyJYUrcdh2fw4rHYcFm6XrP2XTtvffV0W81q95YzaRpM4ibJtI0kYaB02HH67BeUN98hBBbpZRT3vf3ng763QEdKeUyIcQioFRKufJ0t5/obIP+o5+9H2s8Dyn8GBY/Cav31DtLA1usBc1oBmLd/xJIYQImEgMEwPHHSp7wUyCFGu88X4QUdB37o28y0f03DWHaEXiRwkPUmY1EEPC1MuPOGQwam91jdeiNoC+l5Kkv/4a26BgArNEWdOMwEEWKBFIYHN9udc5duLrORzgWIzQEGkgd0BFSp+uy6qm/dUaKarnn2/97ds9/iqDfG8M7U4Fnu38vBSYBK89gO90fBosAsrKyzu5NlVqPEaoArMR1K1EzgpAd6FYwNSeadCGkFS1hxxr3IEw/UqSA8CGxIoUNhIYUOlLoHAsucOwN9/6/KedT15tJSAlIhDQAE2HG0GQIoQXx2BvJmjMae1o6BxuLObimuMee/cEHHwR6NugDDP74ODpr2qjbuodExI6p5yKFFSmsmNrx+YnUOXfhO/ahLGTXOSpkAqSJINF1zspTL14faS/t8fOrN4K+/4T/n5i39sO2I6V8HHgcunr6Z/O1sDe+Sp5Kb3117Qv9pS39oh23Legf7aCfvB7dzm9bPt3jj9gbUzaDQNo5bFcURVF6SW8E/c0c680PA1ac4XZFURSll/R40JdSLgOGCSHmdv9/JYAQYsUHbVcURVF6X6/M0z/ZNEwp5bwP2q4oiqL0PpWGQVEUZQBRQV9RFGUAUUFfURRlAFFBX1EUZQDpldw7PUkI0QAc6et6fIgAcOarXlyY+ktbVDsuLP2lHZA8bcmXUmac+McLPugnAyHElpPluEhG/aUtqh0Xlv7SDkj+tqjhHUVRlAFEBX1FUZQBRAX9nvF4X1egB/WXtqh2XFj6SzsgyduixvQVRVEGENXTVxRFGUBU0FcURRlAVNA/R0KIhUKIuUKI+/u6Lueiuw1JneZaCOHvfj0WCiF+2tf1OVfdr8nc/tAWgGRvhxCiRQixNdnboYL+OThuvd+VQPBouuhk1E9SXN8CpHWn7z667GZSEkJMAuZ1vy6ThBDD+rpO56L7vZHUbQBullJOllI+0NcVORcq6J+bqXSt8wvH1vtV+oiU8vHupTahK8Ak7QeZlHKblPIBIYQfKJVSln5ooQtU9wdW0tb/OP5k//AFFfTP1Yeu96ucf91vzOZkDpTHmULXEqPJbFg/eS3SgGYhxGN9XZFzoYL+uVHr/V6YFkopF/d1JXpC9/CO/+hQYrIRQsztJ0OHR79JBukayk3K1wNU0D9Xar3fC4wQYuHRldmS+RqLEOKnx12TSObORXP3xeiFdC2TmpRDoEKIRccF+qY+rcw5UkH/HPSn9X67T+gpydyD6X4dfto9w2JrX9fnHD0GlHa3yX/ctYqk0n1tYiVdH1onDocmk+c4brLG0ckCyUjdkasoijKAqJ6+oijKAKKCvqIoygCigr6iKMoAooK+oijKAKKCvqIoygCigr6inINknXeuDFwq6CvKWRJCTJJSbuvreijKmVDz9BXlDHXn9pnEsbtkn+u+PV9RLniqp68oZ6A76+Xi7jsyS+m6U/ObfVsrRTl9Kugrypm5BXi2+/fm7p/JnF5AGWBU0FeUMxc8bjx/EZDUKykpA4sa01eUM9Sd/dIPbKMrb7+6mKskDRX0FeUs9Kc88crAooZ3FOXsNH/4Lopy4VE9fUVRlAFE9fQVRVEGEBX0FUVRBhAV9BVFUQYQFfQVRVEGEBX0FUVRBpD/A9mxa5q2jc/AAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot for each iteration\n",
    "print(angle)\n",
    "fig_iter = plt.figure()\n",
    "ax_iter = plt.subplot(111)\n",
    "for x in posterior_marginal:\n",
    "#     print(x[0],x[1])\n",
    "    _ = ax_iter.plot(x[0],x[1])\n",
    "ax_iter.vlines(angle,ymin=0,ymax=2)\n",
    "_ = ax_iter.legend(['Iter ' + str(x) for x in range(len(posterior_marginal[0:10]))])\n",
    "_ = plt.grid()\n",
    "_ = plt.ylabel('Posterior')\n",
    "_ = plt.xlabel(r'$\\phi$')\n",
    "_ = plt.title('Posterior Marginal for some iterations')\n",
    "plt.savefig('iterative_posterior.png',dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 213,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAESCAYAAAAR2wXeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAeu0lEQVR4nO3dXXBb6Xkf8P8Dgh8gJRKiPqi1Xa+W2m13PW7SkaikjT9Cr6m4k2knmSy1zkzaXnRsKtNe1DdZTWZ60U4z3Wid6fiqjrS9yU3H2tX4Ism0tcmNkXidNJWo9cfaq/VK3A99WFpJxKEEAiRA4OnFeQ91hAWBA+Icnhc4/98MhwAOgPO+A+qPV895z3tEVUFERMmQirsBRES0cxj6REQJwtAnIkoQhj4RUYIw9ImIEiQddwOIghKRPIAnVNUx92cBnAQwDmAJwFe9bb7nXwTgAMgCuKSqp3zblgAs+3ZxSlUvbbHfhu8TFRGZBDCvqoej3A8lj3DKJnULf+h7ga+qx822IwBe9Ydkgy+J0wDuqepL9duC7rf+fdpo+6SqLgXvLVE0WN6hbnUawAnvjhmhnzeBvJVzAI6HsO/tvM+rIeyXqGMMfeo6ZlTvNBilnwMw0+SlXwYwH0ITNt9HROZEZN78zJjHJn2PzZsvoiPm9gvmOS+IyKJ5LGt+5s3jr5r7i74+N9rPI68JoV+UAKzpUzeahFuPr7dktvm9KiLe7Ut1JRn/Nniloi185H3Ml89xX4lpUUS+CPc4w2lVXfBeICIzdaWo46p61NTuTwM4BfcLa15VT4hI1vfarfYD/2uatJ1oE0OfulGjcAcafxmcaFK3b7YtyHNPwv3fhWcBJoQBnBGR8wDONKjlfxnAuG907vXF2eI4wVb7WWjyGqKGWN6hrmPq95PeaNhXx/eCME7jZoTv/a9hcYvnnVLVE+bnqHlseYvnNtzPNl5DxNCnrvVVPDw4ek9E5uGOiF/cwTa8avbpmQXwijdTx0zrvOj7cvJKNuf8rzMlHsCdEhp4Py1eQ9QQyzvUlVT1vIjAhD0ATAF4sY1yDQAsioj/+S+q6vk22rAgIkd8B1xPmemkz4uIV2NfMo8tAHhNRC6q6knvoKx5zjyAs9vYT3ar1xBthfP0qSeY0fIigKOcD0+0NYY+9QyvTMLQJ9oaQ5+IKEF4IJeIKEEY+kRECcLQJyJKEOunbO7bt08PHToUdzOaWl1dxcjISNzNCEWv9IX9sEuv9APonr4sLi7eVdX99Y9bH/qHDh3CxYsX425GU7lcDtPT03E3IxS90hf2wy690g+ge/oiIu83epzlHSKiBGHoExElCEOfiChBGPpERAnC0CciSpBIZu94l3ODe7WfUw22z8JdEvYILwBBRLRzQh/p+y7ttgD3uqCTddtnAXe5WACO7wuCiIgiFnroq+olVT1l1vpearDi4TE8vKTdEoAjYbeBiIgai/LkrCk0vqpP/YUf9tY/QUTmAMwBwMTEBHK5XOiNC1OhULC+jUH1Sl92oh9f+9rXAADf+MY3ItsHPw/7dHtfIgt9c7WfEyIyW3c1IgcPr++51WvPwlxJaGpqSm0/+61bztALolf6shP9yGbd8UuU++HnYZ9u70sUNf3TZqQONA74C3g42p+Ee6k4IiLaAVFM2TwDYMkcoM2aUfvmtUzNqH/SO4BrDugSEdEOCL28Yw7cegdqF3yPH/fd5jRNIqIY8OQsIqIEYegTESUIQ5+IKEEY+kRECcLQJyJKEIY+EVGCMPSJiBKEoU9ElCAMfSKiBGHoExElCEOfiChBGPpERAnC0CciShCGPlFIVDXuJhC1FOXlEol6mgK4/9iv4Ct/dgHv3l3F9XwJ//FffAr/+p8+HnfTiLbEkT7RNlWG9yP/+K/jyocFPHVgNwbSKbzxfj7uZhE1xZE+0TYV9j4D1Kr49r/7DMZHBvA7//0HuHV/Le5mETXFkT7RNqgqVvc9jczK+xgfGQAAHBwbYuiT9Rj6RNtw6QMH1cExjNx7a/OxA7uHcHuFoU92Y+gTbcNf/OgmpFbBcP7K5mMHx4awWq6isL4RY8uImmPoE7Vpo1rDX/74F8jkl5CqljcfPzg6BAC4xdE+WYyhT9Sm/7u0jLuFdYzcu/zI4xMm9G+zrk8WY+gTtekvfnQTuwbTyOSXHnn84BhH+mQ/hj5RG8obNfzvN3+B3/jUBFL6aO1+s7zDkT5ZjKFP1IZfrJRwf20D/+zw3o9sywz0YXQozfIOWS30k7NEJAtgxtw9pqqnGjwnD2AJwEKj7US2yhcrAIC9uwYabj84NsTyDlktipH+8wDGVfU8AIjIXIPnnFDVowx86jb5ojtbZyzTOPQnRoc40ierhT7SV9WzvruTAM40eFpWRCZVdanBNu+LYg4AJiYmkMvlwm5mqAqFgvVtDKpX+hJVP/7uplvHf+fNN+A4DgA8up/iOt6/Ww1t3/w87NPtfYls7R0RmQSwvEWwjwNYFpEzqnqyfqP54jgLAFNTUzo9PR1VM0ORy+VgexuD6pW+RNWPpdffBX78M3xp+rP402wWAB7Zz2L5bfzg5hV89nOfR7qv8/9I8/OwT7f3JcoDubONAh1wQ11VHQCOiMxG2AaiUDmlCkSA0Ux/w+0To0OoKXC3UG64nShukYS+iMyq6kvm9kzdtjlf0N+LYv9EUXGKZYwO9aMvJQ23c9om2S700Dchf1pEFkVk0ff4vLn5CtwR/gwAeAd8ibpBvljBnuHGo3zg4QlaPJhLtoriQO4CgMMNHj9ufjsAFszDC/XPI7KZUywjO9x45g4AHBgdBMDQJ3vx5CyiNjjFCrJNRvr7RgaRTgnn6pO1GPpEbcgXy9jTZKSfSgkO7B5kTZ+sxdAnakOrkT4ATIzxBC2yF0OfKKDyRg2F9Y2mI33AncHD8g7ZiqFPFNBKyV13p+VIf3QIt++v70STiNrG0CcKyDHr7jSbvQO40zYL6xu8bCJZiaFPFJC3wmazefoAL5tIdmPoEwXkrbDZqqbPyyaSzRj6RAE5m8sqtxjp87KJZDGGPlFAjlfeGWk9ewfg+jtkJ4Y+UUD5YgX9fYKRgb6mz+NlE8lmDH2igLx1d0Qar7Dpt2dkYHOKJ5FNGPpEAeWLZWRb1PM9mf4+FMvViFtE1D6GPlFATrHScuaOZ3igD8Uy5+mTfRj6RAEFWXfHMzKY5kifrMTQJwqo1Qqbfpn+PpQY+mQhhj5RAKra1kjfLe8w9Mk+DH2iAIrlKsrVWst1dzyZAZZ3yE4MfaIAnFKwdXc8wwN9KPFALlmIoU8UQH412AqbnuGBPhQrVahqlM0iahtDnygAbwmGoDX9zEAfVIG1Si3KZhG1jaFPFEDQFTY9IwNpAOBcfbIOQ58ogHZr+hmzPg8P5pJtGPpEATimpj/WxoFcAChVGPpkF4Y+UQD5YgXDA30YTDdfYdMzzJE+WSod9huKSBbAjLl7TFVPNXjOLAAHwBFVfSnsNhCFzWnjbFwAyPSzpk92imKk/zyAcVU9DwAiMuffaAIfqroAwBGRmY++BZFd8sVy4Jk7gK+8w5E+WSb0kb6qnvXdnQRwpu4pxwCcM7eXABwBsOB/gvmimAOAiYkJ5HK5sJsZqkKhYH0bg+qVvoTdjw9ulTCUxiPv6TgOADTcz42CO1Xzwg9/gr7bb217v/w87NPtfQk99D0iMglgWVWX6jZl6+7vrX+t+eI4CwBTU1M6PT0dSRvDksvlYHsbg+qVvoTdj/98MYdDHxvF9PSRzceyWfdPudF+rueLwOvfwxNP/kNMH/vktvfLz8M+3d6XKA/kzqrqyQaPOwDGI9wvUejcFTaDl3ceztNneYfsEknoi8isd4C2Qc3+Ah6O9icBzEfRBqKwVGuKlVLwC6gAnKdP9go99E3InxaRRRFZ9D0+DwDmAO+k92VgDugSWevBWgWqwdfdAYDBdAop4YFcsk8UB3IXABxu8Phx321O06SukffW3Ql4fVwAEBEMc3llshBPziJqwSl6K2wGD33ALfGUKpynT3Zh6BO1sGLW3RlrY6QPuHP1V9c50ie7MPSJWthu6Gf6eclEsg9Dn6iF+x2M9FneIdsw9Ila8Eb6o22G/sggD+SSfRj6RC2slCoYTKcw1B9shU1Ppr+PUzbJOgx9ohZWSpW2SzuAuU4uQ58sw9AnamGlVGl7uiYAZDhPnyzE0CdqobORPg/kkl0Y+kQtrJQ2th36pUoVqhpBq4i2J1Doi8ihaJtBZK/7pUrbM3cA94xcVWCtUougVUTbE3Skf1pEdkfaEiJLbbe883B5ZZZ4yB5BF1zLA7gkIgsAVgCoqv5hdM0issNGtYbC+vbKO/7llT9ypSCimAQN/dPmhyhR7q+5o/Tt1vQBoFThDB6yR6Dyjqq+C2AGwAsAnjP3iXredtfdAR6GPqdtkk2CHsj9JoB7AF4CsCIi51q8hKgndBL6mX7W9Mk+QQ/kjqvqt1X1XVV9GcCeKBtFZItQRvpcXpksEjT0RUSeFZFREXkO7sXNiXpeKKHPmj5ZJGhN/3kARwG8DOAJc5+o53VU3vEO5LK8QxZpGvoi8qL5/ccA9gJ4F8A+73GiXnd/m8sqA/55+hzpkz1aTdl8xfz+AYDrEbeFyDrbXVYZeHSePpEtmoa+qr5hbv4rAF9R1QfRN4nIHivF7Z2NCwCD6RRSAq6pT1YJenKWg4dn5DoAwDNyKQm2uwQDAIgIhrm8MlkmaOj/sfnxcNlASoROQh9wSzycp082CTpl84sATgGYM2fjckkGSoROQ59XzyLbBA3946r6+777+WZPFpEZEZlvsj0vIosiwi8PslrHI/1+hj7ZpZ2Ts34HQNb8bkpVF1o85YSqHlXVUwH3TxSL7a6l73EvpMLyDtkjUE1fVZ8XkT8AIHCXZPj9Vq9pISsik6q61GijiMwBmAOAiYkJ5HK5DncXrUKhYH0bg+qVvoTRj5oqHqxvIH/7BnK5Ox/Z7jjuienN9rO+WsLKSvPnNMPPwz7d3pdAoS8io6r6dd/9Q6r6Xgf7HQewLCJnVPVk/UZVPQvgLABMTU3p9PR0B7uKXi6Xg+1tDKpX+hJGP/KrZeA78/jlZ57C9Gef+Mj2bDYLAE338z8/uIgPlouYnv78ttrAz8M+3d6XwFfOanG/Lap6VlUdAI6IzHbyXkRR6WQJBg8P5JJtmo70ReQJuAE/IyKTcMs7HTGlm2VVPQ93uWYiK4UR+pmBNKdsklWajvTNUsrPA3hRVb+kqr/h/TR7nRm9T/lH8b7ZPK/AHeHPmH2c76wLRNHYDP1hjvSpdwQ9OWsecGv5AE4COKeqP9zqySbIz9c9dtz8dgB4s3tazfIhik1Y5Z1SpQpVhUjH/1Em6ljQmr53sPUM3KDm/HrqeeGUd/qgCqxVamE1i6gjQUN/r4j8EwCiqq8BWImwTURWCGWk3++ttMm6PtkhaOifgTvanzP3L0TTHCJ73C9VMLDNZZU9w4NcU5/sEvTkrNfM7J1TInLVP2efqFd1ugQD8PCSiSVeMpEsEWikLyLfhDu98iUAKyJyLtJWEVkgzNBfXWd5h+wQtLwzrqrfNlM4XwawJ8pGEdkgjNDP9Lv/meaFVMgW7Sy49qyIjIrIczAXUiHqZWGO9FnTJ1sECn1zgtZRAP8DwBPmPlFPCzX0WdMnS7RahmEUbtArgK/yAC4lSSjlHe9ALqdskiVazd55Ge50Te/2l6NtDpEdqjXFg7WNjtbSB4CRAU7ZJLu0Ku+Iqr5mTshqerUsol7yYM09MSsb0kifoU+2aBX6V323GfqUGGGcjQsAg+kUUsIzcskerco7J81qmAJg0qyaKQBUVZ+KvHVEMQkr9EUEwwNpjvTJGk1DX1XHd6ohRDZxip0vq+zJDPRxnj5ZI+g8faJEyRfLAIDxkYGO34tr6pNNGPpEDdwrmNAfDiP0efUssgdDn6iBfLGMlHRe0weAsUwa90sMfbIDQ5+ogXurZewZHkAq1fnVrrKZATilcgitIuocQ5+ogfxqGXtCqOcDQHa4f/PAMFHcGPpEDSyvlkM5iAu4M4CcUgWqGsr7EXWCoU/UwPJqOZSDuIBb3ilv1HidXLICQ5+ogXwx3PIOANb1yQoMfaI6tZoiX6xgb1ihb2YAsa5PNmDoE9W5v1ZBtaahjfS9s3oZ+mQDhj5RneVVtwwT3kjffZ8VlnfIApGEvojMiMh8k+2z5jkvRLF/ok54oR96TZ8jfbJAJKGvqgtbbTMrdXrPccwqnkTW8EI/tNk7mwdyGfoUv1ZLK0fhGIBz5vYSgCMAHvmSEJE5AHMAMDExgVwut5Pta1uhULC+jUH1Sl866cffXXfD+ec/WcS9K1uPixzHAYCW+1FVpAX48eWryOm1ttrCz8M+3d6XOEI/W3d/b/0TVPUsgLMAMDU1pdPT0zvQrO3L5XKwvY1B9UpfOunHz3JXgDffxm9+8dc3r3zVSDbr/ikH2c+ev13A2P4DmJ7+pbbaws/DPt3elzgO5DoAuE4/WSu/WsZQf6pp4Lcrm+FSDGSHOEL/Ah6O9icBbHnAlygOy6sV7B0ZDPU9uf4O2SKq2TuzAKa8g7bmsXkAUNXzcC+9OGPub3nQlygOy6vr2DPS+ZLKfmOZAR7IJStEUtM3wX6+7rHjvtsvRbFfojAsFysYj2Ck/9ObK6G+J9F28OQsojrLq+sYD+HauH6s6ZMtGPpEdfKrldBOzPJkh/tRqlSxVuG1cileDH0in/WNKgrrG6EtweAZMyd63Wddn2LG0Cfyya+6oRz6SD/Ds3LJDgx9Ip+wF1vzcP0dsgVDn8gnXzSLrYW07o7HW2nTKXKlTYoXQ5/I55632FpUI32WdyhmDH0in3xEoe9dSGWF5R2KGUOfyGd5tQwRYCwT7jz93YNp9KWE18ml2DH0iXyWV8sYy/Qj3RfuPw0RwRhP0CILMPSJfJaL5dBLO57scD9r+hQ7hj6Rz3KhHNoVs+plM/2s6VPsGPpEPvlIR/oDrOlT7Bj6RD7LqxGGPmv6ZAGGPpGhqsgXy6EvweAZG2Z5h+LH0CcyHqxvoFLV0Jdg8GQzA2YftUjenygIhj6R4Z2YFfYSDB7vrFyutElxYugTGVEtweDhUgxkA4Y+kbFcMCP9qGr6Ga60SfFj6BMZt+6vAQAOjg5F8v5ZUzZa4bRNihFDn8i4ni+hv09wYHe4F0X3ZDnSJwsw9ImMG04JH8tmkEpJJO/PC6mQDRj6RMaNfBEfz2Yie//dQ/0Q4YFcihdDn8i4ni9FGvp9KcHoUD9WePUsihFDnwjA+kYVHz5Yxyf2DEe6H660SXGLJPRFZFZEZkTkhS2250VkUUROR7F/onb9wnFn7nx8T3QjfYDr71D80mG/oYjMAoCqLojIpIjMqOpC3dNONHiMKDbX8yUAiLS8AwB7dw3iplOKdB9EzYQe+gCOAThnbi8BOAKgPuCzIjKpqkuN3kBE5gDMAcDExARyuVwEzQxPoVCwvo1B9Upf2u3HX193R9/X3/4RcteC/QfYcRwAaGs//WvruPrhBv7qe99DSlrPEkrq52Gzbu9LFKGfrbu/t8FzxgEsi8gZVT1Zv1FVzwI4CwBTU1M6PT0deiPDlMvlYHsbg+qVvrTbj0vffRspuYLf/tI0+gNeKjGbdf/U29nPreEP8J33foKnfvlX8Q/GWx8/SOrnYbNu70sUNX0HbqhvSVXPqqoDwPHKQURxup4v4eDoUODA367DB3YBAK7cKUS6H6KtRPEXfgEPR/uTAOb9G0Vkzhf09yLYP1HbrjulyGfuAMDh/W7oX/2QoU/xCD30VfU8gEkRmTH3FwBARLzwfwXuCH/G93yiWN3IlyKfuQO4K3iOjwzgKkf6FJMoavpQ1ZcaPHbc/Hbw8MAuZ/BQ7DaqNdy6v4ZP7EDoA8Dh/SO4+uHqjuyLqB5PzqLEu3V/DdWaRj5d03N4/y6O9Ck2DH1KvBveHP0dGuk/eWAX7q2WN6/URbSTGPqUeN6JWTtxIBfwHczlaJ9iwNCnxLthzpB9bCyai6fUY+hTnBj6lHg38iXs3z2Iof6+Hdnfx/dkMJhO4QqnbVIMGPqUeNed4o7N3AHcJZaf2DeCq3c4g4d2HkOfEu9GxOvoN/LkAc7goXgw9CnRajXFTWdtx2bueA7v34Vry0WsVao7ul8ihj4l2p3COsrV2o7N3PEcPrALNQXeu8cSD+0shj4l2uZ0zZ0u72yuwcPQp53F0KdEu54vAti5E7M8T+wbgQg4g4d2HEOfEu2H1xwMplP4ZIC17cOUGejDx7MZHsylHcfQp0T7/jt38auTe3dsjr7fkwd24ee3H+z4finZGPqUWDedEq58WMDnn9oXy/6nHt+Dy7ce4MMHa7Hsn5KJoU+J9fo7dwEAn40p9J99egIAkLt8J5b9UzIx9Cmxvn/lLvbvHsQ/mtgdy/6feWw3HhsbwmuXb8eyf0omhj4lUq2meP2dO/jcU/sgIrG0QUTw7NMH8P137mJ9gydp0c5g6FMi/fTmfeSLFXz+qf2xtuOLzxxAsVzF3y8tx9oOSg6GPiXS37zj1tE/82Q89XzPrx3eh6H+FP7q8oextoOSg6FPifT9d+7gmcdGsX/3YKztGOrvw2cO78Nrl29DVWNtCyUDQ58SZ3V9A4vv52Obqlnv2WcO4NpyiSdq0Y5g6FPi/P2791CpKj4Xcz3f8+zTBwAAr73FEg9Fj6FPibK8WsYf/eVb2DPcj6lDe+JuDgDgsbEMPvXYKL7z01ss8VDkGPqUGKVyFV/5swu47pTw8r+ZimXpha08d/QTuPSBg9P/520GP0UqHXcDiHZCtab4D996A29cc/DN3zuCqUPjcTfpEf/2M4ewdKeAP/3rq9g9lMa//8KTcTeJelQkoS8iswAcAEdU9aV2t4fhv/6vt/De3VVMjA7hwO5BjO8awK7BNIYH0tg1mMbuoTRGBtPIZvqRHe6P7QQditZKsYI///FNnLvwAd68cR//6V9+Cv/804/F3ayPEBH8l9/6NIrlKr7+HXe0/7u/8sm4m0U9KPTQN4EOVV0QkUkRmVHVhaDbw7JWqeL9e0X8v/eW4RQrTZ87PNCHj2UzODg6hKH+Pgz1pzCY7kN/nyDdJ0inUqbtpg+QzfsC4Nr1dbxe+Bl64Xvj2rV1/GD1Z3E3oyV/BUQB1FRRrbk/xXIVV66t4b+9+Tou33qA8kYNTx/cjdPP/WN8+Zi9QZpKCb4++0solav4k+/+HH/y3Z9jYlgwdfMSdg2kMZBOYSCdgv/PrFv+5rrl7yqIoH2pr9JVVbFRVWzUaqhUFRvVGjbM3+xW/ui3P429u8KdVixh1w9F5DSAc6p6SURmUDeab7XdPGcOwBwATExMHP3Wt77VUZvKVcVqRbFeBdY2FKUNYK3q/n5QVtwt1XCvpHDWFZUaUKm6v6vqflDVmhssfo+GjgLokn99Lahq1/yvpz78UgL0CTDQJxjpq2FsKI2DI4Jf+1gaj4+muqZfNVVcdWp4x6ni8p0ybq+lUK4ClZpio/bwed1U+e+mv6tW2ulL/d9onwB9KXF/m59Uauv0+IOpIWSHtnfo9Qtf+MKiqk7VPx5FeSdbd39vm9uhqmcBnAWAqakpnZ6eDq1xUcjlcrC9jUH1Sl+6vR/Pmt/d3g9Pr/QD6P6+RDF7xwHQ7ChZq+1ERBSRKEL/Ah6O5icBzLe5nYiIIhJ66KvqeQCTpl4P7yCtiMw3205ERNGLZMpmo2mYqnq82XYiIooez8glIkoQhj4RUYIw9ImIEoShT0SUIKGfkRs2EbkD4P2429HCPgB3425ESHqlL+yHXXqlH0D39OVxVf3IRSOsD/1uICIXG53u3I16pS/sh116pR9A9/eF5R0iogRh6BMRJQhDPxxn425AiHqlL+yHXXqlH0CX94U1fSKiBOFIn4goQRj6REQJwtDvkIjMisiMiLwQd1s6YfrQ1ctci0jWfB6z5gptXc18JjO90Bdg86p5XUtE8iKy2O39YOh3wH+9XwCOt1x0N+qRJa6fBzBulu/2LrvZlUTkCIDj5nM5IiKTcbepE+bfRlf3AcAJVT2qqqfibkgnGPqdOQZgydxeAnAkxrYknqqeNZfaBNyA6dovMlW9pKqnRCQLYElVl1q+yFLmC6tr2++T7fYvX4Ch36mW1/ulnWf+YS53c1D6TMG9xGg3m+yRz2IcwLKInIm7IZ1g6HeG1/u106yqnoy7EWEw5Z2sV0rsNiIy0yOlQ+9/kg7cUm5Xfh4AQ79TvN6vZURk1rsyWzcfYxGR075jEt08uFg2B6Nn4V4mtStLoCIy5wv6e7E2pkMM/Q700vV+zR/0VDePYMzncNrMsFiMuz0dOgNgyfQp6ztW0VXMsYkFuF9a9eXQbvIKfJM1vMkC3Yhn5BIRJQhH+kRECcLQJyJKEIY+EVGCMPSJiBKEoU9ElCAMfaIOdOu8c0ouhj7RNonIEVW9FHc7iNrBefpEbTJr+xzBw7NkXzGn5xNZjyN9ojaYVS9PmjMyl+CeqfmH8baKKDiGPlF7ngdwztxeNr+7eXkBShiGPlH7HF89fw5AV19JiZKFNX2iNpnVL7MALsFdt58Hc6lrMPSJtqGX1omnZGF5h2h7lls/hcg+HOkTESUIR/pERAnC0CciShCGPhFRgjD0iYgShKFPRJQg/x/L71WOeZ6KvQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax_qinfer = updater.plot_posterior_marginal(smoothing=0)\n",
    "plt.vlines(angle,ymin=0,ymax=3)\n",
    "plt.grid()\n",
    "_ = plt.ylabel('Posterior')\n",
    "_ = plt.xlabel(r'$\\phi$')\n",
    "_ = plt.title('IQPE Posterior')\n",
    "plt.savefig('final_posterior.png',dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 207,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAESCAYAAAAWtRmOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de1xUdf4/8NeHu9wFYRS8ICAoXkG8pokFpVtWtnTbtsu2hW3Xrb5pu9WvvfT9ttpt22pLar/ftv1WWpaW9W0NVOxiogJe8QKCoqAgNxWQ63x+f3CoicswDDNz5sx5PR8PHs4czjnz/jA4L87nc87nCCkliIiITLmpXQARETkfhgMREfXAcCAioh4YDkRE1APDgYiIemA4EBFRDx5qF0DaJIRYDiBNeRqt/Fui/PuRlDLThq9Vp+y71mTxCillvq1eQ0uEEMEANkspp1uwbh2AsVLKejvWEw0gS0oZY6/XIMcTvM6BBksIsRJAjZRylZ32b/cPOC2xNhyEENFSypL+trGwBpvti5wTu5WI9OMjJ90XOSGGA9mcECJYCJElhFguhPhIeZ5n8v1jJo+XCyHylPWDbfA6PzxX1slQlmUJIVJ7W6eX/f5kG2XZR0KIpO5t6F6/Jfvvq90mr5ln+rNQ6slTvpYri4OFEKuVZSst+FmtBJDUVZultfdWU/d99fL+9vbzC1b2YXHNpDIpJb/4NagvACsBLDd5HgxAdi1TnueZfP+Y8m8SOvuqgc5xi9V97L8OQJbpl5nXMX2ehM7xj6795AGIMl2nl9fqbZtgAOkAVpqu01v93Wsw8xp9thtARtcyZV3Tn12wyWsEd/18zLxWncl6eeZqMFe7aU297OuH99fMz8/imvnlHF8ckCZ7qZf9j0HcBCDE5C/saDPr3iB7H3Po/jqmz5cBWGvyvWwAKf3U1ts2qVLKdcpfuyuUuteaqb+/tve6nTKwm4ofB/q71v2hHtk5dhAMoMTk52E6UG8pi2rvoyZzev35Kf8OtmZyIIYD2Yul//lXSCmzbfg6/b1usAXrdBei/FvS9WEppVwhhJiBbvUrH9yW7L/7dkkA3gJwDzrPzFrWz/a2GJw3W7sVNfWl6+fHEwo0hGMOZC/dPwhMxxO6/kpdC5MPHOWDd7CvY/r8I/z0Ay0dwIZetjHV2zYfKo9XK9/rOkunr/r7+xDsbbtkAGtl5+m50d3WvamX17CKyViGJbX3VVP3fZky9/MjDeGRA9md0hVSIoTIApCvfEFKmd81cKmsmgWgr+6YPCGE6QfXc+jsqjD3utlCiCSTwdIV6OeDu7dtTLpCstH54XeDmfr7vb6jj+3WAcgSQqQpNdaarLvapJ5spe3WyAawWQixW0q5zILas3urqfu+0Plz7Wpbrz+/gZ5sQOrjdQ5ERNQDu5WIiKgHhgMREfXAcCAioh4YDkRE1IOmz1YaNmyYjIqKsmrbxsZG+Pn52bYgDdFz+/XcdkDf7WfbO9uel5dXLaUMM7e+psMhKioKu3fvtmrbnJwcpKSk2LYgDdFz+/XcdkDf7WfbUwAAQogT/a2vyW4lIcQSIUTmuXPn1C6FiMglaTIcpJQbpZQZQUFBapdCROSSNBkORERkXwwHIiLqgeFAREQ9MByIiKgHhoPGNbd14OyFFrXLICIXw3DQuGe/KMTs5zbjqQ37Ud3AkCAi29BkOPA6h04dRokv95+BIcAbH+w8iZTnc/D61mI0t3WoXRoRaZwmw4HXOXQqKKtDTWMrfvezCfjqkUsxJyYUz286goUv5OCT/FMwGnmvDiKyjibDgTplHaqEp7vAgvgwxIT5463bk7EmYzaG+Xvj0Q/3Yslr32J7cbXaZRKRBjEcNCyrsBKzo0MR6OP5w7LZ0aH49P5L8MrN01Df1IZfvJ2LX7+zC8VVF1SslIi0huGgUcfONqDkbCPSEgw9vufmJnDttEhsfmwBnlg8HjtLa3HlX7/Bk+v388wmIrKIpmdl1bPswkoAwOUTeoZDFx9Pd9y7IAY3TB+Jv20uwnu5ZdhQUI7fpMRgHO8dTkRmMBw0KquwEhMjAhEZPKTfdUP9vfHHayfh9rlRWPnlYbzw1VEEewtU+5/Ajcmj4OnOA0gi+il+KmhQTUML8srqeu1SMicmzB+Ztyfjw2VzEDZE4Mn1B5D20jZ8uqecZzYR0U8wHDRo8+EqSAmkmulSMmfm2BD8fpYP/nFHMnw83fHwmj246tVvsfVwFSS7m4gIDAdNyi6sRESQDyZGBFq9DyEELp9gwP89NB9/vWkaGlva8at3duHG1d9j1/FaG1ZLRFrkdOEghIgWQiQJIZYLIaLVrsfZNLd14JuiaqQmGCCEGPT+3NwErkuMRPajC/Dn6ybheE0Tbnjze9z1zi4UVpy3QcVEpEVOFw4AkgCUAMgGkK5yLU7n26JqXGzrGPB4Q3+8PNxw2+wx2PZ4CpYvisfu47W46tVv8PCaApyoabTpaxGR87NLOAghUoUQWd2WpSvLl5vbVkq5TkpZDyAVwDp71Kdl2YcqEeDtgVljQ+2yf18vD9yXEotvll+GexfEYNPBM7j8xW14cv1+nDnXbJfXJCLnY5dwkFJmmz4XQqSbLK8XQqR2Le/2Faws7woGdn6bMBolsg9VYUF8GLw87HvQF+TriRWLxuPrxxfilpmjsXbXSVz6/Fb8ceNBVF1gSBC5OmGvs1OEEFlSyjTl8UoAa6WU+coHf5KUclUf26UCWIHOrqUsKeW6bt/PAJABAAaDYfqaNWusqq+hoQH+/v5WbauW4voOPLujGfdO8cbsiMFdojLQ9p9tMuKzY234rqIdHgK4bLQHfjbWC4Hegx/3cDQtvve2pOf2s+2dbV+4cGGelDLZ3PqOugguuNvzPvtElKOLbDPfzwSQCQDJyckyJSXFqoJycnJg7bZqyf33YXi4leA3SxcgaIhn/xuYYU37bwBwvLoRf9tShA0F5dhWLnHH3CgsuzQaQ/28BlWPI2nxvbclPbefbU+xeH1HDUjXAwix1c70ej+HrMJKzIoOGXQwDEbUMD+8dOM0fPXIAqQlGLD662OYt3ILXvzqCM41talWFxHZlqPCYRd+PHqIBpBlZt1+6fF+DqXVjSiuarD6wjdbiw33x99uScSm316KlPhwvLqlGPNWbcEr2UU438yQINI6e52tlA4g2WQgeh2A6K6B6O4D1tS/ron2bH0K62DFGQLw+q1J+PLh+ZgTHYqXs49i/sqteH1rMRpa2tUuj4isZJcxByUM1nVb1usAtDWEEEsALImNjbXVLp1e1qFKTBgRiJFDfdUupVcTRgQi8/Zk7D91Dn/NPornNx3B29+U4O750bh9zhgE+KjXFUZEA+eMF8H1S2/dSrWNrdh9vBZpE8LVLqVfk0cG4R93zsD6++Zi6qhgPL/pCOat3Iq/bWZ3E5GWaDIc9Gbr4SoYJZCWMFztUiyWOHoo3vnVTHx6/yWYETUUL2UdxSV/2YKXs45y4JpIAzQZDno7WymrsBLDA30wKdL6ifbUMnVUMN6+YwY+f3Ae5kSH4pXNRT+c3VTX2Kp2eUTUB02Gg566lZrbOvB10VmkJoTbZKI9tUyKDELm7cn44qF5mDduWOfZTSu3YNW/D6OWIUHkdHgnOCf3/bEaNLV2aKpLyZyJEUF445fTceTMBby6pQhvbDuGd7Yfx21zxuCe+dEY5u+tdolEBI0eOeipW+mrwkr4e3tgdrTNriF0CvHDA/DaL5Lw1W8vRVqCAW99XYL5K7diQ0G52qURETQaDnrpVjIaJTYfqsSCuDB4e7irXY5djDME4JWbE5H16AKMHeaHF7OO8G50RE5Ak+GgF/vKz6HqQgtSE5z/FNbBignzx6/njcXJ2ovYfaJO7XKIdI/h4MSyCs/A3U1gYbzrhwMALJo0HEM83fFJPruWiNSmyXDQy5hDdmEVZkaFINhXOzOeDoaftweunGjAF/sq0NLeoXY5RLqmyXDQw5hDWU0TjlReQKqTzaVkb0uTRuJ8czu2Hq5SuxQiXdNkOOhB1iFloj0nmYXVUS6JCUVYgDe7lohUxnBwUlmFZxBvCMDoUOecaM9ePNzdcM3UCGw9UsUrqIlUpMlwcPUxh/qmVuw6Xud003M7ytLESLR1SHy+/7TapRDplibDwdXHHLYeqUKHUepuvKHLxIhAxBn8eUEckYo0GQ6uLquwEuEB3pgS6Zrh1x8hBJYmjkTeiTqcqGlUuxwiXWI4OJmW9g5sO3IWqQkGuLlpd6K9wbp2WgSEANbz6IFIFQwHJ/P9sRo0tnbo7iyl7iKCh2D22FCsLyjndBpEKmA4OJnsQ5Xw9XLHnJhQtUtR3dKkSJyoaUJ+Wb3apRDpjibDwVXPVpJSIruwCpeOC4OPp2tOtDcQiycNh7eHGwemiVSgyfs5SCk3AtiYnJx8j9q1DMa5pjYcrbqAI2cu4GjlBRw6fR5nzjfr9hTW7gJ8PHHFxOHYuK8CT1+dAC8PTf4tQ6RJmgwHrWlqbUdxVcMPIXCksgFHz1zAmfPNP6wT4O2BuOEBuH3OGCye7Bo39rGFpYkR2Li3AjlHqnDFRP5ciByF4WAHre1GbDt6Fp/trcC+U/Uoq21C15iql4cbxoX7Y25MKOKHByBueADiDQEYEeSj6duA2sv8cWEI9fPC+oJyhgORAzEcbERKifyyemwoKMfn+ypQ19SGUD8vzI4OxfWJIxE/3B9xhgCMCfWDu45PUR0oT3c3LJkagfdzy3CuqQ1Bvp5ql0SkCwyHQTpe3Yj1BeXYsKccJ2qa4O3hhismDsfSxAjMHxcGT3f2kw/W9UmReGf7cfzfgdO4ZeZotcsh0gWGgxVqG1vxxb4KfFJQjoKyeggBzIkOxQMLY7Fo0nAE+PCvW1uaHBmEmDA/rM8vZzgQOQjDwUIt7R3ILqzC+oJy5BypQrtRIt4QgCcWj8e10yIwImiI2iW6rM7pNCLxwldHcbK2CaNC9DVTLZEaNBkOQoglAJbExsY65PUaW9px2z9ykV9WD0OgN+6aNxbXTYtEQkSgQ16fgGundYbDhoJyPHj5OLXLIXJ5mgwHR17n0NzWgbv/uRt7T53DizdMxXWJkRxQVsGoEF/MHBuC9QXleOCyWJ7ZRWRnHC01o7XdiPvey8eO0hq8cMMU/Hz6SAaDiq5PjERJdSP2nXKtK+OJnBHDoQ/tHUY8snYPthyuwrPXTcLSxJFql6R7iyePgJeHG2dqJXIAhkMvjEaJFR/vxxf7T+Opqybg1llj1C6JAAQN8UTqhHBs3FuBtg6j2uUQuTSGQzdSSjzz2UF8nH8Kj6TG4e750WqXRCaWJo5ETWMrvj56Vu1SiFwaw8GElBJ/+fdh/GvHCSy7NBoPXe6Ys6HIcgviwjDU15NdS0R2xnAw8dqWYqzeVoJfzh6NJxaP5xkxTsjLo3M6jazCSpxvblO7HCKXxXBQvP1NCV7MOorrEyPxp2smMRic2NLESLS0G/Hv/WfULoXIZTEcALyfW4ZnvziExZOGY1X6FF3fu1kLpo0Kxthhfvik4JTapRC5LKcLByFEtBAiSQixXAgRbO/X21BQjic37EdKfBheuTkRHpwoz+kJIXDdtEjsKKlFef1FtcshcknO+ElYC6BEeRxizxfadPAMHvtoL2aNDcGbv5zOO41pyNLESADAp3s4ME1kD3b5NBRCpAohsrotS1eWLze3rZTSIXeT33+2HQ++X4ApI4Pw9h0zeM9mjRkd6ovkMUOxPr8csutOSkRkM3YJBylltulzIUS6yfJ6IURq1/JuX8FCiAwlILIBpNujvp2ltXi1oAWx4f54586Z8PfW5BRTurc0KRJFVQ3IPlSldilELkfY668uIUSWlDJNebwSwFopZb4SDElSylV9bJekPEwFkNn9SEIIkQEgAwAMBsP0NWvWDLi2igYj/nmgCfcn+iHQW5+Dzw0NDfD391e7jEG52C7x7I6LqGiQWDrOE1dHe8LNgrPMXKHtg6Hn9rPtnW1fuHBhnpQy2dz6jvqTufvAcmhfK0op85WH+X18PxNAJgAkJyfLlJQUqwqK8M+Btdu6gpwc12j/pfPb8fv1+/HJngrUuQXj5ZumIdjXy+w2rtJ2a+m5/Wx7isXrO2oEth52HlwmffLz9sBfb5qGP183Cd8V1+Cqv32LfaccMmxF5NIcFQ678OPRQzSALDPr9ksIsUQIkXnuHKdups5TW2+bPQYf3TsHAJD+xvf4144THKgmGgR7na2UDiDZZCB6HYDoroHo7gPWAyWl3CilzAgKChp8seQypo4KxucPzsOcmFA8veEAHlm7B02t7WqXRaRJdhlzUMJgXbdlvQ5AW8PRtwkl7Rjq54X/uXMGXttajJezj6Lw9Hm88cvpiAnT5yAkkbU0edUXjxzIHDc3gYcuH4d375qJ6oZWXPPqt/hi32m1yyLSFE2GA5El5o8Lw+cPzkPc8ADc/34+/rjxIFrbeZMgIktoMhw4IE2WiggegrUZc3Dn3Cj8z3fHcXPm96htZkAQ9UeT4cBuJRoILw83/OGaiXj1lkQcOXMBz2y/iPyyOrXLInJqmgwHImssmRqBTx+4BD7uArdk7uA4BJEZmgwHdiuRtWLDA/D0nCGYFBmE+9/Pxxs5x3g9BFEvNBkO7FaiwQj0Enjv7llYMjUCK/99GE98vB9tHRyHIDLF6UhJl3w83fHKTdMQFeqLV7cU41R9E/5+63QEDfFUuzQip6DJIwciW3BzE3jsing8nz4FuSW1SH9jO07WNqldFpFT0GQ4cMyBbOmG5FF4966ZqDzfjKV//w57TnLiPiKLwkEIEWXfMgaGYw5ka3Njh+GT++ZiiJc7blr9Pb7czzOZSN8sPXJYKYQIsGslRCqLDQ/A+vsuwcSIQNz3fj5Wb+OZTKRfloZDHYB8IcQbQojnhBDP2bMoIrUM8/fG+/fMxs8mj8BzXx7G79cf4JlMpEuWnq20UvlyCpyVlezJx9Mdr96ciDEhvvh7zjGcqmvC67cmIdCHZzKRflh05CClLO3+Ze/C+qmHYw5kV25uAssXjceqn0/B98dqcEvmDtQ3tapdFpHDWHy2khDiHqVb6T/sWRCRM7lxxii8dUcyiqoacOvbuQwIUl2HUaKgrA7FVQ12fR1Lz1Z6A0ANgFUAzgkh1tq1KiInsjA+HJm3TWdAkGrOXWzDxr0VeHTtHsz4z2ws/ft2vLPdvh04lo45hEgpP1EevyWEuMFeBRE5oxQlIDL+lYdb387Fe3fPQrCvl9plkYuSUqKoqgFbDldhy+Eq5J2oQ4dRItjXEylxYVg4PhwL4sLsWoOl4SCEEJcB2A0gDQCvEiLdYUCQPTW3dWD7sWpsOVyFrYfPorz+IgBgwohA3LsgGpeND8e0UUPh7iYcUo+l4bACQDqAewHslFLeaL+S+sezlUgtDAiypXNNbfjywGlsOngG24/VoKXdCF8vd1wSOwwPXBaLlPgwjAgaokptlobDXwDcI6V83p7FWEpKuRHAxuTk5HvUroX0hwFBg9HU2o6swkps3FuBbUfPoq1DYnSIL26ZORqXjQ/HrOgQeHu4q12mxeFQByBPCJENpUtJSvk7u1VF5OQYEDQQre1GbDt6Fp/trUB2YSUutnVgeKAP7pwbhWumRmJSZCCEcEx3kaUsDYcP4EQXwRE5AwYEmdNhlMgtqcGneyrw5YHTON/cjqG+nrg+KRLXTI3AjKgQuDlo/MAalobDfQDullJesGcxRFrDgCBTUkrsPXUOn+4pxxf7TqPqQgv8vNxx5cThWDItAvNih8HTXRuTYQ+kWymf3UpEPTEg6GJrBzbsKcc/tx/H4TMX4OXhhoXxYbhmaiQunxAOH0/1xxAGSpNzKxE5GwaEPpXVNOFfO45j7a6TON/cjgkjAvHc9ZNx1ZQRmp+Ly2w4CCGmSSn3dJ9LSQhxPQBV51cicjYMCH0wGiW+Ka7Gu9uPY8uRKrgJgUWThuPOuVFIHjPU6QaWrdVf59cPXUfKFBpdltmnHMvwTnDkrFLiw/HW7ckoqmzAgx8UoJ3TfbuM881t+J/vSpH60jbc8d87sfdUPR5cGIvvVlyG13+RhBlRIS4TDED/3UrCgscOx+scyJktiAvDs0snYfm6fVi16Qh+/7MJapdEg1BUeQHvfn8Cn+SfQmNrBxJHB+OVm6dh0aThTnE9gr30Fw6yj8d1dqiFyGXcmDwKB8vPIfPrEiSMCMR1iZFql0QDVFBWh5eyjuKbomp4ebhhyZQI3DF3DKaMDFa7NIfoLxzShBC70HmkEC2ESFYeJ9q9MiKNe+rqBBw6cwErPt6H2HB/TIrk/Ue0oLjqAp7fdASbDlYi1M8Lj18Zj5tnjEKov7fapTlUf+Ew1iFVELkgT3c3/P3WJFzz6rfIeHc3PntwHobp7ANGS8rrL+KvWUfxcf4p+Hp54NG0OPx63lj4eVt6UqdrMdtqKSVHfIkGYZi/N1bfloz0N7fjvvfy8d7dszRzEZRe1Da24u9bi/HujhOABO66ZCzuWxiLED99n2mmz0gkcqDJI4Pwl59PxiNr9+LZzwvxx2snqV0SAWhsacc/vi1F5tclaGptR/r0kXg4NQ6RwerMgupsGA5EDrA0cSQOlp/H29+WYmJEEG6cMUrtknSrtd2I93NP4LWtxahuaMWVEw14/Mp4xIYHqF2aU2E4EDnIE4vH4/CZC3hqwwGMM/gjcfRQtUvSlQ6jxPaKdjz1Yg5O1V3E7OgQvHX7eL4PfbD0HtJ3CyHeFEI8pzznPaSJBsjD3Q2v3pIIQ5A37v3fPFSdb1a7JN0orW7Ez9/Yjsx9LQga4ol375qJD+6ZzWAww9KRsTQp5b0mz+1+nYMQgnM5kcsZ6ueFzNuScf5iO+793zy0tHeoXZJLk1LivdwT+Nkr36C0uhEZU7yx8YF5uDQuzKWuZrYHS8NBKPMpBSv/2pUQIgmAPq40Id2ZMCIQL9wwFfll9Xjm04OQUva/EQ3Y2QstuPufu/Hk+gNIjhqKTb+9FHMjPJz6HgrOxKJwUO4ZHYPOC+BCuh1F9CCESBVCZHVblq4sX25hbfUWrkekOVdNGYH7UmKwZtdJvJdbpnY5LiersBKL/vo1vimuxjNLEvDPX83E8CAftcvSFIsGpIUQlwHIA5APQAohLpNSbulrfSllthBihcn26SbLo4UQqcrj9G6bZgNIBVCCziuyo6WUJQNsE5EmPHZFPApPn8cfPjuIOEMAZo4NUbskzWtsacefPy/Eml0nkTAiEGtunoZxBp6FZA1Lz1aabvI4BsBQAH2GQy9mAOgaxC4BkAQgW0q5rpd11wkhgsFuJXJx7m4Cr9yciOte/w73vZeHzx6YhwieY2+1vBN1ePTDPSirbcJvUmLwSGocvDx4waG1hDX9nUKI5/q7E5wQIktKmaY8Xg1gtZQyXwiRis4B7hXmtjez3wwAGQBgMBimr1mzxprdoKGhAf7+/lZt6wr03H5na3tFgxF/+v4iRvi54cnZPvCwc5+4s7V/sNqNEp8da8PGY20I8RHImOKN+JDeZ0t1tbYPhGnbFy5cmCelTDa3vqXdSm/ip7OyJsPkXg8WqAdgk2NmKWUmgEwASE5OlikpKVbtJycnB9Zu6wr03H5nbHvg6Ao88H4BTnpH4VeX2HdKM2dsv7VKzjbgkbV7sPdUE36eNBJ/uCYBAWbuwOZKbR+ogbbd0m6lLwCc6noipfzNwMrCLvzYTRQNIMvMuv0SQiwBsCQ2NnYwuyFyGldNHoE1sSfxctZRXDstUvfz+vSn8xTVMjz7RSF8PN3xxq1JWDx5hNpluRRLO+R+CaBYSlkgpSzob2VloDnZZCB6HToHmFOV59nWFqxsv1FKmREUxCmQyTUIIfD01QloaGnHy1lH1S7HqRmNEk9tOICnNhzAzLGh2PTbSxkMdmDpkUMdgHwhRDaUU0zNjTkoYbCu27JV1hbZHY8cyBXFDw/AL2ePwf/uOIFbZ4/G+OGBapfkdIxGiSc3HMAHO8uwbEE0nlg0nhez2YnZIwflFFYAWAngCgCr0Nnfn2nnuszikQO5qkdS4xDg44k/bSzkxXHdGI0Sv1+/Hx/sLMP9C2MYDHbWX7fSMgCQUpZ2/3JAbUS6M9TPC4+mxWH7sRp8VVipdjlOw2iUWPHxPqzZdRIPXRaL/7ginsFgZ5o8CVgIsUQIkXnuHO9FRK7n1lmjMS7cH//5xSHOvYTO2VQfX7cPH+WdwsOXj8OjDAaH6C8c0oQQRd2+ioUQRQ6prg/sViJX5uHuhv+3JAFltU3472+Pq12OqjqMEo9/tBcf55/CI6lxeCQtTu2SdKO/cMiWUo7r9hUrpRznkOqIdGr+uDCkTjDgtS1Fup3au73DiMc+3INPCsrxWFocHk7lx44jsVuJyEk9edUEtHYY8fymI2qX4nDtHUY8+uFebNhTgcevjMeDlzMYHK2/cBjUxWr2wm4l0oOxw/xw1yVj8VHeKew9qZ9Jits7jPjt2j34bG8FViwaj/sX8pR1NZgNBynlW44qhIh6euCyWAzz98KfPtfHqa1tHUY8vGYPPt93Gr9bPB6/SYlRuyTd0mS3EpFeBPh4YvmV45F3og6f7a1Quxy7ausw4qEPCvDF/tN46qoJWLaAwaAmTYYDxxxIT9Knj8SkyED85cvDaGptV7scu2htN+KB9/Px5YEzePrqBNw9P1rtknRPk+HAMQfSEzc3gWeWTMTpc814c5vr3fvKaJR48IN8bDpYiWeWJODX8+w7Ky1ZRpPhQKQ3M6JCsGRqBFZvO4by+otql2NTa3adxKaDlXjqqgl2n66cLMdwINKIJxaPhxDAc/93SO1SbKa6oQV/+fIQZkeH8IjByWgyHDjmQHoUGTwEyy6Nwef7TmNnaa3a5djEf31xCBfbOvDsdZM5JYaT0WQ4cMyB9OreBTEYEeSDP248iA6jtk9t3V5cjU8KynHvghjEhuvz1p3OTJPhQKRXQ7zc8bufTcDBivNYl3dS7XKs1tLegac2HMCYUF9e5OakGA5EGrNkyggkjxmK5zcdwYXmNrXLscqbOSUoqW7En66dBB9Pd7XLoV4wHIg0RojOU1trGlvx+tZjapczYKXVjXg9pz8TCWsAAA3USURBVBhXTxmBBXFhapdDfWA4EGnQ5JFBuDJhOD7OPwWjhsYepJR4esMBeLu74f9dnaB2OWSGJsOBZysRAYsmDcfZCy3Yc0o7k/J9trcC3xZX4/FF8QgP9FG7HDJDk+HAs5WIgIXx4XB3E8jSyO1Ez11sw58/P4QpI4Nw66wxapdD/dBkOBAREOTriVljQ5CtkXB4ftNh1Da24L+WToa7G69pcHYMByINS0swoKiqAcerG9UuxayCsjq8l1uGO+ZGYVIkj/i1gOFApGFpCQYAcOqupfYOI36//gDCA7zxKO8BrRkMByINGznUFxNGBDp1OLyz/TgOnT6PPyyZiAAfT7XLIQsxHIg0Li3BgN0nalHb2Kp2KT1U1F/ES1lHsTA+DIsmDVe7HBoAhgORxl2RYIBRApsPOd/Rwx83HoRRSvzp2kmcWE9jGA5EGjcxIhAjgnycrmspu7ASmw5W4qHLx2FUiK/a5dAAaTIceBEc0Y+EEEidYMA3RdVobutQuxwAQFNrO5757CDGhfvj7nm85acWaTIceBEc0U+lJRhwsa0D3xZVq10KAOCV7CKU11/Efy6dDC8PTX7M6B7fNSIXMDs6FAHeHk7RtXT4zHm8/W0pbkweiZljQ9Quh6zEcCByAV4eblgQH4bNhytVn4jvhU1HEejjgd8tnqBqHTQ4DAciF5GWYEB1QysKTqo3EV97hxHfH6vGVVNGYKifl2p10OAxHIhcREp8ODxUnojvYMV5NLZ2YNbYUNVqINtgOBC5iKAhnpgdHYqswjOq1ZBbWgMAmBXNsQatYzgQuZDUCeE4drYRJWcbVHn93JJaRA/zQ3gA79WgdQwHIheSquJEfB1GiZ3Ha3nU4CKcLhyEENFCiI+EEMvVroVIa0YO9UWCShPxHTp9Hhea2zne4CKcLhwU90gpV6ldBJEWpSUYkFdWh+qGFoe+bm5pLQCON7gKu4SDECJVCJHVbVm6sry/I4JaANHK+rzunmiA0hIMkBLYcqjKoa+bW1KD0SG+GBE0xKGvS/Zhl3CQUmabPhdCpJssrxdCpHYt7/YVDCBZSpkPoARAuj3qI3JlEyMCERk8BFkOnKXVqIw38Ipo1+HhoNeZAWCt8rgEQBKAbCnluu4rCiFKhBBJAJIBZDqoPiKX0TkRXzjW7j6Ji60dGOLlbvfXPFp1AfVNbZjFcHAZjgqH4G7P+xyxklKWKA/ze/u+ECIDQAYAGAwG5OTkWFVQQ0OD1du6Aj23Xw9tD2/rQHObEW9u2IrE8J/+N7dH+7NPtHU+qCpCTs4xm+7blvTw3vdloG13VDjUA7DJnxRSykwoRxTJyckyJSXFqv3k5OTA2m1dgZ7br4e2z2034s39WTjtFoZHUqb+5Hv2aP+H7+UhIqge6YsXOvVNffTw3vdloG131NlKu/Dj0UM0gCwz6/aL93MgMs/Lww0p48Ox+VAVOuw8EZ+UEjtLazErOtSpg4EGxl5nK6UDSDYZiF6HzjOQUpXn2ea27w/v50DUv7QEA2oaW1FQVmfX1zl2tgHVDa0cb3AxdulWUsJgXbdlNrtuQQixBMCS2NhYW+2SyOWkxIfB071zIr7kKPt9cO8o6bq+gRe/uRJnvQjOLB45EPUv0KdrIj77ntKaW1qL8ABvRIXyPtGuRJPhQESWSUswoKS6EcVV9pmIT0qJ3JIajje4IE2GAwekiSyTOqFzIr5sO10Qd7ymCVUXWjje4II0GQ7sViKyTETwEEyKtN9EfLklnfdvmM35lFyOJsOBiCyXOsGA/LI6nL1g+4n4cktrMczfCzFh/jbfN6lLk+HAbiUiy/0wEd9h2x49dI03zBwbwvEGF6TJcGC3EpHlEkYoE/HZuGvpVN1FVJxr5v0bXJQmw4GILCeEQFqCAd8UVaOptd1m+91RwvtFuzKGA5EOpCUY0NJuxDdF1TbbZ25pLYJ9PREXHmCzfZLz0GQ4cMyBaGBmjg1BoI+HTbuWcktrMDMqBG5uHG9wRZoMB445EA2Mp7sbFo4Px5bDVTDKwU/EV1F/ESdrL3LKDBemyXAgooFLSzCgtrEVRXXGQe9rZ9f9onnxm8tiOBDpxIK4MPh7e+CrrhvzDEJuaQ0CfDwwYUSgDSojZ6TJcOCYA9HABfh44u75Y5FX2YG9J+sHta/cklrMjAqBO8cbXJYmw4FjDkTWuXt+NAI8gVWbDlu9j6rzzSipbuQprC5Ok+FARNbx9/bA1TFe+K64Bt8VW3daa+4P4w0cjHZlDAcinVk4ygMRQT5YtekIpBVnLuWW1sDf2wMTIzje4MoYDkQ64+Uu8NvUOOw9WY9NBwd+3UNuSS2mjxkKD3d+fLgyvrtEOnR9UiRiwvzw4ldH0GG0/OihpqEFRVUNHG/QAU2GA89WIhocD3c3PHZFPIqqGrC+oNzi7XZyvEE3NBkOPFuJaPAWTxqOyZFBeDnrKFraOyzaJre0FkM83TFlJP/vuTpNhgMRDZ4QAssXxaO8/iLezy2zaJsdJTWYPmYoPDne4PL4DhPp2LzYYZgTHYrXthSjscX8dN71Ta04UnmBU2boBMOBSMeEEHh8UTxqGlvx39+Wml13Z2ktpAQn29MJhgORziWNHoq0BAMyvy5BXWNrn+vlltbC28MNU0dxvEEPGA5EhP+4Ih4Nre14c9uxPtfJLa1B4uhgeHu4O7AyUgvDgYgQPzwAS6dF4p3tx3HmXHOP759vbkNhxXmewqojDAciAgA8khYHo5R4ZXNRj+/tPl4Lo+T9ovVEk+HAi+CIbG9UiC9+MXM0Ptx9EqXVjT/5Xm5JLbzc3ZA0eqhK1ZGjaTIceBEckX08cNk4eLm74aWsoz9ZvqO0FlNHBcHHk+MNeqHJcCAi+wgL8MZd86KwcW8FDlZ0Hpk3tLTjQPk5jjfoDMOBiH4i49IYBA3xxAubjgAA8k7UocMoOd6gMwwHIvqJoCGeuHdBDLYeOYudpbXILamBu5vgeIPOMByIqIc750YhPMAbq/59GLmltZgcGQQ/bw+1yyIHYjgQUQ9DvNzx4OXjsPtEHfJO1LFLSYcYDkTUq5tnjMKYUF8AwGwORusOw4GIeuXp7oanr0pAnMEfMzgTq+44ZSeiECIDQAmAWillvtr1EOlVaoIBqQkGtcsgFThdOHQFg5QyW+1aiIj0yi7dSkKIVCFEVrdl6cry5f1sPh1AsLJ+kj3qIyIi8+wSDt3/6hdCpJssrxdCpHYt7/YVDKAeQDaAfAA32aM+IiIyz1HdSjMArFUelwBIApAtpVzXfUUhxHMAbgRQC2C1g+ojIiITjgqH4G7P+zwvTkpZDyCzr+8rYxIZAGAwGJCTk2NVQQ0NDVZv6wr03H49tx3Qd/vZ9hyL13dUONQDsMm5cFLKTCjhkZycLFNSUqzaT05ODqzd1hXouf16bjug7/az7SkWr++o6xx24cejh2gAWWbW7Rfv50BEZF/2OlspHUCyyUD0OgDRXQPRgz1NlfdzICKyLyGlVLsGqwkhzgI4YeXmwwBU27AcrdFz+/XcdkDf7WfbO42RUoaZW1nT4TAYQojdUspktetQi57br+e2A/puP9tueds5txIREfXAcCAioh70HA59XkuhE3puv57bDui7/Wy7hXQ75kBERH3T85EDERH1geFAREQ96DIcBjB9uMsRQtQJIfKEECvVrsVRBjmFvOb10X6X/z0QQnRN/Z9u2k49vPdm2m7x+667cOhr+nAduUFKOV1KuULtQhzF0inkXVUfMxLo4ffgRgAhXbM/CyEydPTe92i7stzi91134YDO6cNLlMdd04frSbAQIlrtIlSm998BQAe/B1LKTGWiTqBzTrds6OS976PtwADedz2Gg8XTh7uoEAC1Qgg93ytD778DgI5+D5QPw1opZQl09t53azswgPddj+Fgs+nDtUj5i6IenYfU6WrXoxJd/w4Auvs9SJdSLlMe6+29N237gN53R93PwZnYdPpwLVH6HWuVfsgatetRkW5/BwB9/R4IIdKllKuUx6nQ0XvfS9u7jiIset91d+Rg6+nDNeZDmAzC9XabVldk7ynknV339kMnvwdK+1YqZ+fkAfp573trOwb4vvMKaSIi6kF3Rw5ERNQ/hgMREfXAcCAioh4YDkRE1APDgYiIemA4kEtTJiDLM3k+qCkjTLcXQiQ5y8R1rj4VBjkew4H0ZrAf5j9sL6XMd6KJ65wipMh1MBxIN5QpmlOFEKuFEMHKso+UrwwhRLTyvY+EECkmjzN62970yMFkP6uV50lCiCwhxErTI5du9azutk1v++h6vFyZZrrHfntrF9Fg8SI4cmnKh+VmKeV05XmWlDJNebwcQImUcp3yIbwSndMpTFfmn+naxzEpZUwv2ycBuAmdUxHkSymzlauQQwDsBrBSSpmmBEiW6dW4XYHTNXOmUktv+1gmpVzW9X0Atb3t17QuIlvgkQPpWQyArg/ZrjDI7goG5a/0jD63/tEMdH6QA50f4NOVx10zYfY2j810k23M7aM35vZLZBMMB9KzPHT+5b3CZOygKxgyAKRCmY+mn/3sUtaF8u8xC177GIDkAewjxoJ9EtkMw4F0R+nXj1a6dNK6+v67rVaCzi6j3/W1fddzZebLm0TnrTjTumbCNEdZZ3rXOENv+5BS5qNzwryVsGCa6e51EQ0GxxyIiKgHHjkQEVEPDAciIuqB4UBERD0wHIiIqAeGAxER9cBwICKiHhgORETUw/8HNS06icQGVhwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig_true_err = plt.figure()\n",
    "ax_true_err = plt.subplot(111)\n",
    "ax_true_err.semilogy(performance['true_err'])\n",
    "ax_true_err.grid()\n",
    "_ = plt.xlabel('Iteration count')\n",
    "_ = plt.ylabel('True Error')\n",
    "_ = plt.title('True Error over each Iteration')\n",
    "plt.savefig('true_error.png',dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 206,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAESCAYAAADnvkIDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de3wU53kv8N+ju0BaCXGRsbguF2NsbJCEL7Fjgy2cq53YFjhNYzdNY5EmaZu2McTpOWnS5iQRaXOSNqmL0uOTJk1jA03TnDiJiwBhJ3ZsI3ExNncBAsxFIK2QBLrtPuePmZWG1Ui7Wu3s7OX3/Xz2s7uzM7PPuwPz6H3fmfcVVQUREVGoDLcDICKixMQEQUREtpggiIjIFhMEERHZYoIgIiJbTBBERGSLCYJGJCLtItIoIlstj/IItvOKyLF4xGjz3aN+r4gcE5HNNstrRcSxa75FpDz4Gzr1HWG+v1hEGiNct11Eih2Ox7V/IxS5LLcDoIR3v6r6wq0kIl5VbQYA83leLL7cut8Y8tosqwIQtpzjUAtgdSS/ZSKJ5e/v1L8Rcg5rEBQrw/4qT+D9Pi8i1cE3Zq2o3oHvsfImW3IwxfL3d+rfCDmECYLGzGwesDY71QIINqGsszZnBF+bTTjB53Xmuo2WfQb31Rhs3gjdr7lsnaXZK7jeRnPZRgAlERRhC4DHLO8fA/B8SBntvscuxmD5gjHU2vxe6wAEf7Nyc1mNZX9Vln0Ff0O7ZrCIYrLsv9F8rDMXF48Wp833hf397WKO5FiGNnmN8nuM+tuSw1SVDz5sHwDaAWy1PszltQCqQtZttLwuDr43XyuAYvO9Aqg2X28FUB6ynxoAG0fYb7klBi+AjTCahrZaYw5TpmM2+220bmv3PSPFaFM+2+8Pfq9l/5ut32/uJ7ivdTbbjyWmcpvjEVGcluMeXC/c7z9azKMdS+u/kXC/R9iY+XDmwT4ICseu3XwrgI0isgXGCSBcG3WzZR/NGGrOaYb5F7+IeGGc7FeNsp/HAJRY/rr2wug3sP613RYmlqBgM5M1ntG+Z7QYreWL5PvX4toaS72533oAPlXdYLPNWGK6pkakqj7zL/mxxhlRDKExR3gsrUb7PcYbM40Dm5hozFS1HkP/+SO5MuaaBBOacMxml80AdsH4q3Q061V1tfmoiDBkO8FmpmHNS3bfEybGWPQtBJvGRjsJjiUmO7GI0+73H4w5iphGEvw9krHfJmUwQdCYBa9GUdX1AHZZ2+Oj3GUlgOdVtQk2VxhZ9vs8jL82B+OAUZtZbVkvkj4ImLUeL4wmrqaQj+2+Z9QYx2izdf8AqgFsMl+PdEIcS0zPw9LHYq4btTC/f2jMkR5Lq9F+D3IRm5gonEYRsZ4Avg6jmWG1+b7ZbMKoB7BNRHYBWD/G76gHsFVEVsE42bSFfLZNRHap6lq59l6Craq6QURWmx2eu2A0GUUq2IZ+DVVtCv0eGDWOkWIcE1WtF+O+iGDta72lGWikbSKOyVx3o2X/9TCOWzRG/f0B1NmsH/ZYwvJvJJrfg+JDVDkfBBERDccmJiIisuVIgjCvX642H7bXLsvQMA615vtyMYZBaAxZfs16REQUH071QawBAFWtE5HlIlKjqqFtlavNq2GCSlR1HjB4JYRvhPWIiCgOHKlBqGqdJSF4YT+MQbH16oqQJGAd/+Wa9YiIKD4cvYrJPLG3jXAjVQmANhHZqKrWS+dCaxu261nXh3HHJvLz8ytmzpwZVayBQAAZGenVJcMypweWOfWNt7yHDx++qKpTh33g5G3asLn93madWphDL5jvN0eynt2joqJCo7Vjx46ot01WLHN6YJlT33jLC2CX2pxTHUuxIlKt5u33wcG3LJ/VyNBompcsy4sjWY+IiJzn1FVMVQCCo3deM2Kn+XITAF8wcajqFnN5Ca69sWak9YiIyGGO9EGo0eE8bDIQVV1lPvsw1HFdb/m8GZZb7kdaj4iInJc+vThERDQmTBBERGSLCYKIiGwxQSSg3gE/mlu7EAhwIEUicg+H+04AFzp70HTSh6aWdjSdbMe+Mx3oGwjg2Y9X4r5FpW6HR0Rpigkizgb8ARw814mmlnY0nmxHU0s7TrVdBQDkZGZgyYwiVFfMwL+/1oIzvh6XoyWidMYEEQfHWrvws91nsOtEO/ae9uFKnx8AMK0wFxWzJ+EP7pyDZbMm4eYyD3KzMtHT78e/v9aCzp5+lyMnonTGBOEQVcXLRy7i2d8eR8OhVmRmCBZP92BN5Uwsm1WMitmTUFacDxEZtm1uVgZyMjNw+eqAC5ETERmYIGLsSt8Aftp0Bj945QSOXujC1MJc/HnVQnz09lmYWpgb0T5EBIV5WbjMGgQRuYgJIkbO+K7ih6+ewHOvn0LH1X4sKSvCt9bcig/cMh25WZlj3p8nPxuXrzJBEJF7mCDGQVXR1NKOZ39zAr9+6xxUFe+9+Tp84q65qJg9ybb5KFKevCx09rCJiYjcwwQRhX5/AC/sO4tnf3sc+053wJOXhU/ePReP3zkbMyZNiMl3ePKz2cRERK5ighij3gE/PvWjRuw41Arv1In42w/fjEfLyzAhJ7Y/ZWFeFt7xXY3pPomIxoIJYgz6BgL4zI+bsONQK77y0E14/I7ZyMiIvhlpNJ68bFxmExMRuYgJIkL9/gD+9Ce7UX/gAv7mQzfhiTvnOPp9nvxs3gdBRK7iWEwRGPAH8Lnn9+DXb53Dlz642PHkABid1D39AfQO+B3/LiIiO0wQYfgDir/cvBcv7DuLL75/ET5x99y4fK8nPxsAeCUTEbmGCWIUgYBi3ZZ9+K897+Cp99yAmnuGTZLnmMI8o/WP90IQkVucmpO6WESqzUftCOu0m3NW14ZZVi0iVSKyzolYRxIIKL74n2/iP5pO43NVC/CZlfPj+fXw5Bk1CHZUE5FbnKpBrAFQoqpbAEBEamzWWa2qFaq6fqRlIlINDM5x7RORKofivYaq4ks/34/n3jiFz66cjz+7f0E8vvYaQ01MrEEQkTscuYpJVessb70ANtqsViwiXlVtHmXZcgDPm6+bAZQDqLfuxEw+NQBQWlqKhoaGqGLu6upCQ0MDVBU/PtCH+pYBvH9uNipy3sHOnWej2ud4nO4MAABebdwL/xlnLjYLljmdsMzpId3K7FR5Hb3MVUS8ANpCkkBQCYA2EdmoqmtHWFYcss3k0J2YyagOACorK3XFihVRxdrQ0IB7770X/+uFA6hvOY5P3DUX//ODN45ruIzxeMd3Ffjtdsz0LsSK22Y58h0NDQ2I9vdKVixzeki3MjtVXqc7qastJ/9rqGqdqvpgNB1Vj7DMByNpOE5VseHFQ/iX3xzHH9w529XkAAw1MbGTmojc4liCEJFqVd1gvq4K+awmmBQAXBppGYA3MFSL8ALY6lS8/3m0H880HMNHb5+FLz90k6vJAQAm5mQiQ3iZKxG5x6mrmKoA1JpXJDValgdP8Jtg6XQ2O7OHLTOXey3Lrul/iJV/2HYEPz/WjzWVM/DVD93senIAjDkhOGAfEbnJqU7qegDDbhpQ1VXmsw9Dnc31Iy0zl29wIkargCredX0Wvv7ILY6NrRSNwrwsNjERkWs4FhOAz1UtxPbMM8hMoOQAcMA+InIX76Q2ZSRAs1IoTx4H7CMi9zBBJDBPfhYuX2UNgojcwQSRwArz2ElNRO5hgkhgnrxsdlITkWuYIBKYJz8L3X1+DPgDbodCRGmICSKBBUd07eplPwQRxR8TRAIbGm6DCYKI4o8JIoENThrEjmoicgETRAIbmjSICYKI4o8JIoF58oPTjrKJiYjijwkigbEGQZR+VBVt3X1uhwGAYzEltMEEwXshiFJKT78fp9uvoKXtClouXUFL21W0tF3BqTZjWb8/gIN/+15kZbr7NzwTRAIrGOykZhMTUbJp7+7D8UvdOHmpGycvXbkmAZy/3HvNuhNyMjGrZAJmTZ6AuxdMwaySCRgIKLIyXQrexASRwDIzBIW5WRywjygBqSp8V/oHk8Dxi1dw8lI3TlzsxolLV9BhqfmLANM9eZhZMgH3LJg6mAxmlkzArJIJmDwxJyHmoQnFBJHgPPnZ7KQmclFPvx/HL3bjWGsXmluN5xMXu3H8Yvc1tXsRoKw4H3OnTMSDt07HnMkTjceUiZhZko9ct6sDUWCCSHCFeVnspCZymKriQmcvjl3owrGL3Whu7cKxVuP5jO8qVIfWDSaBh5ZenxJJYDSOJAgRKQYQnId6uaqut1mnHUAzgHpVXT/SNqHrORFvIuOAfUSxEwgo3um4iiPnu3D4fCcOn+/C0QudONbafc2QNhNyMuGdOhHlsyZhdcVMeKdOxLypBZg7ZSLyc1IrCYzGqRrEGgBQ1ToRWS4iNapaF7LO6pA5pkfaJnS9tOLJz8I7vh63wyBKKqqKM76rOHy+E0fMRHDkfCeOXuhCd59/cL1phblYUFqA6ooZg0nAO3UirvPkJWSfQLyJWutOTnyByGYA61W1OWR5NYCm0OWh24y2nrluDYAaACgtLa147rnnooqzq6sLBQUFUW3rpLp9vTjc7sff3Tsh5vtO1DI7iWVOPd39itOdAZyyPM50+dHrHzrBF+UKygoEZQUZKCvIwPXm88Ts1EgC4z3GK1eubFTVytDljvZBiIgXQNsIJ/cSAG0islFV146yje16QWYtow4AKisrdcWKFVHF2tDQgGi3dVLD5bewv+20I7ElapmdxDInL39AcfJSNw6c7cTBc5dx4OxlHDjbiTO+q4PrFE/Ixo3XTYK3qAMryhdhYWkhFpYWoHhCjouRO8+pY+x0J3W13UkdGDyxQ0R8IlKtqlvsthllvbRQmJeFzt4BBAKKjIzU+GuHKJyefj8OnuvEm2c68PY7HXj7bCcOn+vE1X6jeSgzQ+CdMhEVsyfh9++YhRune3DjdR6UenIhIsYJ847ZLpci+TmWIMyT+QbzdZW1H8FsFmozT/aXRtoGgNduvXTiycuGKtDdN4BC885qolTS0+/H22cvY/+ZDrx5ugP737mMI+c7MRAwmr+NWoEHv3fbLCyaXojF0z2YP60Aednp01nsFqeuYqoCUCsiT5uLglckbVXVVQA2Aag014Oqbhlhm2HrORFvIhscsK+HCYKS39U+P94+aySCN88YSeFoaxf8ZjIomZiDm8uKcN+iqVhSVoSby4pQVpzPDmOXOJIgzNrCPJvlq8xnH4BgjaJ+tG1C10s31vGYyorzXY6GKHKqiuMXu7HnlA+7W3zYfaodB88O1QymFBjJ4IGbSnFzWRGWlBVhehGvHkokvFEuwRVywD5KEh1X+rHntA97zGSw55QPvivGv9uC3CzcMqMIa+/14tYZxbhlRvFgfwElLiaIBBdsYurkgH2UQFQVx1q78frxNjS1tGN3SzuOtXYDMIacWDitEO9ZfB2WzSrGslmTMH9aATJ5kUXSYYJIcJwTghLBgD+AA2c78fqJNrx+/BJ2nWjHJXPOgskTc7BsVjEeXlaGZbMm4ZYZRewvSxFMEAnOk88mJoq/nn4/9p3uwBsn2vDa8TY0nWwfHIpiZkk+7r1hKm6bU4Lb5pZg7pSJbCpKUUwQCa6Qc0JQHPQO+LG7xYffHr2I15rbsOe0D30DAQDAwtICfGjp9bhtrpEQphfxYol04fhQG/FUWVmpu3btimrbpUuXori4OMYRxcbJ5X+GwvN7UdLSENP9+ny+hC2zU1hmgwLonzAVV4tmo6doNnoKZ0AzcwANIKf7PPIun0Ze52nkdp5B5sBV+x0nsHQ7zl/+8pfHdSe1iMR/qA2KjQx/LwJZuW6HQUluIMczmBCuFs1CIHsiACD76iUUtO5HfsdJ5F4+hUx/b5g9UbpggjB9+9vfTtjxalZ9ayfmly/AM88+HX7lMUiVMXrGIp3K3NnTj98evYhNO/eh+UoOTl+6AgCYWpiLB+ZPwV3zp+Cu+ZNTsskonY4zYJTXCUwQScCTn82rmCgiJy91Y9uBC9h+8AJeO34J/X5FXiZw14ISPHHnHNy9YAoWTCtgpzJFhAkiCRTmZeFSV5/bYVAC6vcHsOtEO7YfPI/tBy8M3oswf1oBPnHXXKxcNA1dJ/ah6r7lLkdKyYgJIgl48rJx4mK322FQgmjr7sPOwxew7cAF7Dzcis6eAWRnCu7wTsbH7piN+xZNw+zJEwfXb2hhbYGiwwSRBDz5WbzMNc2d8V3Fr948i1/vP4emlnYEFJhSkIv33Xwd7ltUirsXTEFBLv87U2yF/RclIk8CWA1gEoA2AJNhDL29WVX/xdnwCBial1pV2XacRoJJ4YU3z2J3iw8AcON0Dz573wLcv2galpQVcY4QctSoCUJE/hnAJlV9wOaz+0XkGVX9Y8eiIwDGgH0DAcXVfj8m5PCvxFR2uv0KfvXmObzw5lnsOWUkhcXTPXjqPTfg/UumY+6UiWH2QBQ7I55tRKQIxrzQHXafq+o2ANtExKOql50KkK4dsI8JIvWcaruCX+0/ixfePIe9ZlK46XojKXxgyXTMYVIgl4x4tglNDCJyX8jn281nJgeHWeeEKPXkuRwNxUJbdx/+a88Z/GzPO4NJ4eYyD9a910gK1k5mIrdE9OeoiOwCsBVGHwRg3Km/3amg6FqDA/bxXoik1u8PYOehVmxuPIXtBy+g36+46XoP1r93Ed6/5DomBUo4kbZX1KtqxLfxikgxgCrz7XJVXW+zTjuAZnPfwSlJqwH4AJRb5qYetizdDA7Yd5VXMiWjA2cv4z8aT+Nne87gYlcfphTk4OPvmoNHK2Zg0XUet8MjGlGkCaJcRF6EcUIHAITpnF5jrlMnIstFpEZV60LWWW1OMwpgMBFAVetFxGvOQ10cusy6TbrgnBDJp627Dz/fcwZbmk5j/5nLyM4UVN1YiuqKGbhn4VRkZ2a4HSJRWJEmiGE1gNGEJAMvgI02qxWLiFdVg0lnOYDnzdfNAMphXFIbuiz9EkQ+h/xOBsEmpC2Np7Ht4Hn0+xVLyorwlYduwkO3Xo9JE3PcDpFoTCJKEKq6O5qdi4gXQJslCViVAGgTkY2quhZmbcFi8gjLQr+jBkANAJSWlkY9aFVXV5djA16NV5/fGJJ971uHMLPneMz2m8hldooTZb7cq2g43Y8dLQNo71V4coD7Z2bhrrJszCwcAPpOYO8bJ2L6nWPB45z6nCpvpJ3UT8KoRSiADgCqqpEM7lJtnvyHCdYyRMRn6WcoCVnNbpndfuoAYz6IaEdwTPTRH3N2/AqTr5+JFStujNk+E73MTohlmfed9uEHr5zAL/aeRZ8/gHcvmIIn7pyDFTckVhMSj3Pqc6q8kTYxVanqfDNRbAIQtsNaRKotHc1VIf0NNTBqFltg3JUNAG9gqMbghXHVVLHNsrTkyctCJ5uYXNc3EMCv9p/FD145gd0tPkzMycRHbpuJJ+6cg/nTCtwOjyimIk0Q7ebzLhjNOfePtrLZwVwrIsFEErxKaauqroKRZCrN9WAmCojIOsuy+pGWpaPgcBvkjguXe/Dj11rw76+3oLWzF3OnTMRfP7gY1RUzUGheRECUaiJNELWA0RchIpUI02ltnsjn2SxfZT77MNTZXG/5fNhlrOl6aWuowvxsdlLHmapi9ykf/vWVE/jlm2fR71esuGEqPv6uObhnwVSOg0QpL9IEcb+ZGNpV9WkReR68US6uPHlZrEHEiapix6EL+MftR7G7xYfC3Cx87I7ZeOLOORwLidJKpAlilao+JiLfMN+3j7o2xZwnLxtnfMk3eXwyCQQU//32eXx3xxHsP3MZZcX5+MpDN+HRihkcSpvSUqT/6kVEHgFQZD5TnHny2UntFH9A8cKbZ/G97Udx6Hwn5kyegA3Vt+DhZWUJdTUSUbxFeh/EGhF5CoAAKFHVTzkbFoViJ3XsDfgD+K897+B7DUfR3NqN+dMK8O3HluKDt0xHFhMDUeQzyqnqN0XEAwAc4jv+PPnZ6B0IoKffj7zsTLfDSWp9AwH8tOk0/qnhGFrarmDRdYX43kfL8b6br2PHM5FFpDfKbQJQBOPGNYFxw9xjDsZFIYID9nX2DDBBRKmn349tLf34q79rwBnfVSwpK0Ld4xWourGUiYHIRqQ1iDZVXeNoJDSq4IB9nT39mFqY63I0ySUQUPxszxl888VDONvRh/JZxfjqwzdjxcKpnMKVaBSRJohGEfk8rh3N9afOhER2OGBfdF5rvoSvvnAAb57pwC0zivD4QuCPH3kXEwNRBCJNEGth3NDG/1Uusc4qR+GduNiNr//qAF586zymF+Xh248txUO3Xo+XXtrJ5EAUobFMGPQFRyOhURVyToiIdFzpxz9sP4IfvnoC2ZkZ+PwDC/FHd3uRn8N+G6KxcmrCIIqxYBMT74Ww1+8P4N9+dxLf2XYEHVf78VjlTPzFqoWYxjm8iaLmyIRBFHtsYrKnqqg/cAFf/+UBNF/sxl3zJ+Ov3r8Yi6/nVJ5E4+XohEEUOxNyMpGZIWxisnjrnQ589RcH8GrzJXinTsSzH6/EyhumsY+BKEZGvF1URIpE5BkRmTPC53PNz/mnWhyIiDlgH5uYAOD85R488k+v4OC5y/jKQzfhxc/dg/sWlTI5EMXQiDUIVe0QkS8AeFpElmGo/yE4FegxAF/gHdXxU5iXjU7WIAAADYcuoHcggP/89F1sTiJyyKhNTKraAeALgFFjgJEYms3lFGee/CzeB2HaebgVpZ5c3Di90O1QiFLWWMZiOu5kIBQeB+wzDPgD+M2Ri3jPTdexSYnIQY4MWSkixSJSbT5qw6xbaz6Xi8gxEWk0H8Hl7db36cyTl81OagB7T/twuWcA994w1e1QiFKaU7OgrAEAVa0TkeUiUqOqdaErmXNNe823Jao6z1xeDmNgQABYnc5zUVsV5nFOCADYefgiMgS4e/4Ut0MhSmkR1SDMK5q+MdpVTVaqWmdJCF5Y5p227NOLa2+8s67jVdXgZ8XmumnPk88mJsDof7h1ZjGKJ+S4HQpRShu1BiEiS1V1D4DVAL4GYyymGhHZai4flXlib7Oc7K28qlof2oZsU9soAdAmIhtVda3Nd9QAqAGA0tJSNDQ0hAvLVldXV9TbxkvbuT509/mxbfsOZMZgeOpkKHOorj7FvlNX8NC87KhiT8YyjxfLnPocK6+qjvoAsAzAkwDuA7DUXHZfuO3M9daNsLzK8npzyGebR9imFkD1aN9XUVGh0dqxY0fU28bLs79p1tnrf6FtXb0x2V8ylDnUz/ec0dnrf6GNJ9ui2j4ZyzxeLHPqG295AexSm3Nq2CYmNe6iblPV7aq6x7wnolxE7hORpSNtJyLVqrrBfF0V8nGbiFSJSDUAr9nnABEpDtlHjbkOAFwKF2uq44B9RvNSUX42bp1RHH5lIhqXSK9iahKRr5sJwQegI5gw7FY2E0Jt8Ioky/KtAKCqTWr0OZTAuLciqARAm+X9JgC+YIJR1S2RFiwVefLSe8A+VcVLh1tx94IpMWliI6LRRToW03ER+QaMq5NUVb8fZv16APNslq8KeV8HoM7yvhnG3BPB9z4MdXCn/ZVMnvz0HrDv4LlOXOjsxb0LeHkrUTyM5Ua5DgCjJgZylifNm5heOtwKALhnIRMEUTw4cqMcOWNw2tE0HbBv5+FW3FBaiOuKOMcDUTwwQSSRdO6k7u4dwK4T7bx7miiOwt0Hcd9In6nq9tiHQ6MpzM2CCNJywL7fNV9Cnz+Ae9j/QBQ34foggh3N1RjqJJ4MYC4AJog4y8gQFORmpWUn9UuHW5GXnYHKOZPcDoUobYQb7vv7gHHZqqp+M7hcRJ5xOjCyl64D9u083Io7vZORl53pdihEaSPSPohJIvKIiMwRkUcBLHcyKBpZOg7Y13LpCk5cusKrl4jiLNIEsRpGc9N6GM1LI/ZNkLPSccC+nUeMy1vvZYIgiquxJgifqv4deD+Ea4wmpvSqQew81IoZk/Ixd8pEt0MhSiuRJohVqvopy/t2J4Kh8Dz56dVJ3TcQwKvHLuLehVM5exxRnEWaIEREHoExN8MjTgZEo/PkZaMzjTqpG0+2o7vPz/4HIhdElCBUdQ2MJiaBMfPbp8JsQg7x5GWhs3cAgYC6HUpcvHSkFVkZgnfNm+x2KERpJ6KxmETEE3KZ6xxVPeFYVDQiT342VIGuvoHBsZlS2c5DrSifPWnwLnIiip9Im5hqw7ynOBkcsC8N+iEudPbg7bOXefUSkUvCDbUxF0YyqLLMC82eQhcFB+xLh3shXj58EQAvbyVyS7g7qY8DWCMiT1mbmMg9hWlUg3jpSCumFORg8XSP26EQpaVIO6mZHBLE0JwQqV2DCAQULx+5iHcvmIoMzh5H5IqIEoSIPCkiR0XkiIjsEpE3wqxfLCLV5mPU/grr5yLSbk5Tal1Wbc5fvS6SWFPd0JwQqV2D2P9OB9q6+9i8ROSiSDupq1R1PoANAO4HsC3M+mtgXA67BQBEpMZuJXOuaa9l0WpVrVDV9ebn1cDgFKaDc1Ons3SZVW7nIWN4jbsXTHE5EqL0FWmCCN45vQtADYwkMSJVrTPnmwaMBDBsPmmz07s5ZHGxpTMcMAYFDK7TDKA8wnhTVkFeenRSv3SkFTeXeTClINftUIjSVqRzUtcCgKruFpFKAF+IZCPzZN+mqqGJAAC8qlofMnxCCYA2EdmoqmsBFIdsM+xuKbN2UgMApaWlaGhoiCS0Ybq6uqLeNt5yM4G3jjSjIfPMuPaTqGW+0q9oPHkF75+bHfP4ErXMTmKZU59T5Y0oQZhXM0FEPACeH8P+q80T/TXM+SWG1SqCtQ4R8ZnNSz4YSWO02OoA1AFAZWWlrlixYgzhDWloaEC028bbpFe2oWjKFKxYceu49pOoZf71/rMIaBMer6rA7d7Y3kGdqGV2Esuc+pwqb6R3Um8CUASgw1ykAB4Ls021qm4wX4cmhDazP6EYgFdEygFUwqhtbAFwyVzvDQzVIrwAtkYSb6ozBuxL3SamnYcvoiA3C+WzOXsckZsi7YNoU9X3qOoa8xEuOVQBqDWvSGq0LN8KAKraZCaMEgwlgE2wdESr6hYzWXgty4bVOtJRYV42OntTs5NaVfHS4Va8a95kZGdG+s+TiK+yiFAAABFaSURBVJwQaR9Eo4h8HpZOZVX96UgrmyfyeTbLV4W8H2weMtWHPCNYC6EhnrwsXOzqczsMRxxr7cYZ31V8euWwfz5EFGeRJoi1ME7avGMpAXjys9F8sdvtMByx87Bxees9C3j/A5HbIk0Q9aoa0ZVL5DxPXupOO/rS4VZ4p07EzJIJbodClPYiTRDlIvIirm1i+mNnQqJwPPlZ6OwZgKqm1CxrPf1+/K75En7vtlluh0JEiDxBrHc0ChqTwrxsDAQUV/v9mJAT6SFMfK8fb0PvQAD33sDmJaJEEG6476+r6tMwLmkNncJst2NR0aiG5oQYSKkEsfNwK3KyMnDHXM4eR5QIwp1dNpnPY7k5jhw2OGBfTz+uK8pzOZrYeelwK26fW4L8nEy3QyEihJ8PIlhLOKaql4PLRWSOgzFRGMEaRGcKDdh3uv0KjlzowprKmW6HQkQmTjmahArzgkN+p8bd1Jd7+vHpHzchJzMDqxaXuh0OEZnGOuVo6lwyk8Q8+akz5HdX7wA+/uzrePudy9j4eAXmTJnodkhEZOKUo0nIkyLTjl7t8+OPfvAG9p7uwHd/bxnuv5G1B6JEEmkT01bA6HsQka+LyFIHY6IwBpuYknhOiJ5+P2p+tAuvn2jDt9bcivctme52SEQUItIEERyyeyOMITfYB+GivOxM5GRlJG0TU99AAJ/+cRNePnIRGx69BR9aWuZ2SERkI9IEMdmsNYiqbsPQsN/kEmO4jeSrQfT7A/iTnzRh+8EL+NrDS7CaVy0RJaxIE8RGAJ+COXMbjHkayEWe/Kykq0H4A4o/f34PXnzrPP76wcX46O0cUoMokYW7immpqu4xaw3bLB8dczYsCifZBuwLBBRPbdmLX+w7i6fftwh/eNdct0MiojDC1SCeDr4QkWcsy4dNI0rxVZhnDNiXDAIBxV/97E38tOkM/mLVQqy9l3M9ECWDcAlCInhNLvDkZydFE5Oq4iv/7y385PVT+OzK+fjT+xe4HRIRRShcgtARXrc7EAuNQTJ0UqsqvvbLA/jXV0/iyXfPxV8+sNDtkIhoDMIN1rdKRN6AUWPwikil+XrZaBuJSDGAKvPtclUdcbhwEalV1fUjbSMi7TDmoagfbT/pJhk6qb+19TC+//JxPHHnbHzx/Tem1NwVROkgXIKItidxDWDMOS0iy0Wkxpx/+hoiUgXAG2ab1eYc12ThyctG30AAPf1+5GUn3uin//JyM/5x+1F8ZPlMfPnBm5gciJKQqIZO8xDjLxDZDGC9qjaHLA8mhlpVXT3SNiJSDaApdHvLujUwL78tLS2teO6556KKs6urCwUFBVFt64ZtLf340dt9+M7KCSjKje7k61SZ918cwN/v6kV5aSY+szQXGQmUHJLtOMcCy5z6xlvelStXNqpq5bAPVNWxB4zawcYRPqsynzePtg2Mk3/xSPuxPioqKjRaO3bsiHpbN/xs92mdvf4XevRCZ9T7cKLMJy526S1fflEf+NZO7erpj/n+xyvZjnMssMypb7zlBbBLbc6pkd4oF61qVR12SayIVOnIzUbXbKOqdarqA+AzaxOExBywr7t3ADU/bAQA1D1RgYm5qTPbHVE6cixBiEi1qm4wX1eFfNwmIlXmCd8rIuV224hIjSUpXHIq1mQUnFUuUe6FUFV8fvNeHLnQie9+dBlmT+aw3UTJzpEEYSaEWhFpFJFGy/KtAKCqTWYNogRG89FI22yCUXOoMrfb4kS8yagwL7HmhPjejqP41f5zePp9N+LdC6a6HQ4RxYAjbQDmyX/Y7bKquirkfR0A69VNdrfY1oc8E6xNTO7XILYdOI+/33oYH156PT75bg6hQZQqnO6DIIcEm5jcrkEcvdCFzz23Bzdd78E3Hr2Fl7MSpRAmiCSVn52JrAxBp4sJ4nJPP2p+tAs5WRnY+HhlQt6PQUTR42UmSUpEUJiX5VoTUyCg+Nxze9By6Qp+/MnbUVac70ocROQc1iCSmJsD9n1r62FsP3gBX3pwMW73TnYlBiJyFhNEEnNrTohfvnkW391xFI9VzsTjd8yO+/cTUXwwQSQxT37854Q4eO4yPr95L5bNKsbffJhjLBGlMiaIJObJi28Tk+9KH5784S4U5Gbhnz9WgdwsdkoTpTJ2UiexeHZSD/gD+JOf7Mb5jl48t/YOlHry4vK9ROQeJogkFs8axP+uP4yXj1zEhkdvQfmsSXH5TiJyF5uYkpgnPxtX+vzo9wcc/Z49p3x4puEY1lTOwJrlMx39LiJKHEwQScyTZ1QAuxzsqO4d8OOpzXtR6snD//jgYse+h4gSDxNEEovHgH3/uO0ojlzowtceWTI4/hMRpQcmiCTmyXd2wL79ZzrwzM5jeLR8BlbeMM2R7yCixMUEkcSCTUxO1CD6BgL4/Oa9mDwxB19i0xJRWuJVTEksWINwYsC+f2o4ioPnOvH9JypRNIFNS0TpiDWIJFYYrEHEuInpwNnL+O72o/jQ0uuxanFpTPdNRMmDCSKJDfZBxLAG0e8P4Kkte1E8IRtffvCmmO2XiJKPU1OOFotItfmoDbNureV1tTkX9brRlpGhICcLIojpgH11LzVj/5nL+NsP3YxJE3Nitl8iSj5O1SDWACgJziEtIjV2K5lzTXvN19XA4HSlPjMpDFvmULxJKSNDUJibhcsxug/i8PlOfKf+CD6wZDret2R6TPZJRMnLkQShqnXmfNOAkQCGzSctIl4AzZZFyy3vmwGUj7CMLGI1J8SAP4CntuxDQV4WvvIhNi0RkcNXMZlJoE1Vm20+9qpqvWW46OKQzyePsCz0O2oA1ABAaWkpGhoaooq1q6sr6m3dJAO9OH76XFSxW8v8y+N92HuqH5+6NRf7d70a2yATSLIe5/FgmVOfU+V1+jLXalVdG7pQRKrMZiMrH4CSCJZdw6yp1AFAZWWlrlixIqpAGxoaEO22brr+0KtQACtW3DnmbYNlPtbahZ/Vv4wHFpdi/UcqUnqOh2Q9zuPBMqc+p8rrWIIQkWpV3WC+Dk0IbWZ/QjEAr4iUA3gDQzUGL4Ct5vvQZWThyc/G6farUW/vDyjWbdmH/OxMfPXhm1M6ORDR2Dh1FVMVgFoRaRSRRsvyrQCgqk1mwiiBmQDMDm1vsCNaVevtljkRbzIb77SjP3jlBBpPtuOvH1yMaYWc44GIhjhSgzBP5PNslq8KeT/YPGS+32CzzbBlNKQwLyvqTurz3QF883cHcd+iaXh4WVmMIyOiZMcb5ZKcJz8bXb0DCAR0TNsFAopn9/ciOzMDX3t4CZuWiGgYJogk58nLgirQ1Te2eyH+7bWTONQewP/8wGJcV8SmJSIajgkiyQ0N+R15M5PvSh+++eIh3DQ5A6srZzgVGhElOSaIJOeJYsC+Z3YeQ1fvAD6yKJdNS0Q0IiaIJOcZ46xy5zp68IPfnsDDy8ows5CHn4hGxjNEkhtrE9N3th1GQBV/XrXQybCIKAUwQSS5YA2iM4IB+461dmHTrtP4/dtnY2bJBKdDI6IkxwSR5Dz5kU87+q3/PozcrAx89r75TodFRCmACSLJFeRG1km977QPL7x5Fp98txdTCnLjERoRJTkmiCSXlZmBiTmZYWsQ33zxECZNyMaT754bp8iIKNkxQaQAT342OkdJEL89ehEvH7mIz6ycj0Kzz4KIKBwmiBRgDNhn38Skqtjw64MoK87Hx+6YHefIiCiZMUGkgNEG7Pv1/nPYe7oDn6tagLzszDhHRkTJjAkiBYw07eiAP4Bv/vchLJhWgEfKOaQGEY0NE0QK8ORl2d4H8R9Np9Hc2o3Pv+cGZGZwSA0iGhsmiBTgyR8+aVBPvx/frj+CpTOL8cDiUpciI6JkxgSRAow+iAGoDs0J8aNXT+JsRw/Wv3cRB+Qjoqg4NeVosYhUm4/aEdapMh+15vtyETkWnKbUsrzd+p6G8+Rlwx9QXOnzAzDuqv5ew1Hcs3Aq7pw32eXoiChZOVWDWAOgxJxTGiJSY/1QRMoBrDKnJi0XEa+5/jxVrQDwJICN5uqrVbVCVdc7FGvSCw7YF+yHqNvZDN+Vfqx7zw1uhkVESc6pOanrLG+9GDrZBz9vAtAkIsUAmlW1GUCzdZtgcgFQLCJecx2yYR3yOyMD+D+/OY4Hb70eN5cVuRwZESUzRxJEkFkzaBvl5F4JwBeyTU1IgikB0CYiG1V1rc131ACoAYDS0lI0NDREFWtXV1fU27rt+EWj5rDzldfxu7MD6Bvw425Pe9jyJHOZo8Uyp4d0K7Nj5VVVxx4A1kWwzkYA1Zb3m0dYr9a6nt2joqJCo7Vjx46ot3Xb7pZ2nb3+F/p/f9Os855+Qb/4030RbZfMZY4Wy5we0q3M4y0vgF1qc0517ComEalW1Q3m66qQz2ot/RI+GLUEmE1O1vVqRKTafHvJqViTXXDa0W9vO4KsTMGf3r/A5YiIKBU4dRVTFYDa4BVJluVbzZcbATSb6xXrUJNSCYA2y642AfAFE4wO9UuQRbCT2nelH39411yUevJcjoiIUoFTndT1AObZLF9lPls7pestnzcDWGt577N8PrgeXavQrEF48rLwqXuG/exERFFxtJOa4iM3KxO3zS3Bh5eWoWgCh/MmothggkgRm9be6XYIRJRiONQGERHZYoIgIiJbTBBERGSLCYKIiGwxQRARkS0mCCIissUEQUREtpggiIjIlqhlmspkJyKtAE5GufkUABdjGE4yYJnTA8uc+sZb3tmqOjV0YUoliPEQkV2qWul2HPHEMqcHljn1OVVeNjEREZEtJggiIrLFBDGkLvwqKYdlTg8sc+pzpLzsgyAiIlusQRARkS0mCCIissUEAUBEqkWkSkTWuR1LPIhIuzlfeK3bsTjNPK5bQ5al9PEeocwpe8xFpNg8ptXW8qXqcR6lvDE/xmmfIESkGhicR9snIlUuhxQPq1W1QlXXux2I08zjOigdjndomU2pfMzXAChR1S0AICI1KX6ch5XXXB7zY5z2CQLAcgDN5utmAOUuxhIvxSLidTsIl6Tj8QZS+Jirap2qBq/i8QKoRwof5xHKCzhwjJkggOKQ95NdiSK+SgC0ichGtwNxQToebyANjrl5cmxT1WakwXEOKS/gwDFmggB8MH7YtGH+BeKDUfWudjueOEu74w2kzTGvVtW15ut0OM7W8jpyjLNisZMk9waG/trwAtg6yrpJz2yvbDPbLy+5HY8L0up4A+lxzEWkWlU3mK+rkOLH2aa8wdpETI9x2tcgzB/UG+zEGqGDL5VsgqXTLtjRlarMv6QqLZ2WKX+8Q8uMFD/mZrlqzSt4GoHUPs525YVDx5h3UhMRka20r0EQEZE9JggiIrLFBEFERLaYIIiIyBYTBBER2WKCoLRhDnLWaHk/rmEJrNuLSHmiDISXqkNqUPwxQVA6G+8JfXB7VW1KoIHwEiJRUfJjgqC0ZA4BXSUiG0Wk2Fy22XzUiIjX/GyziKywvK6x295ag7DsZ6P5vlxEtopIrbUGExLPxpBt7PYRfL3OHMZ62H7tykUULd4oR2nDPGFuU9UK8/1WVV1lvl4HoFlVt5gn4loYwzNUmOPbBPdxTFXn2WxfDuAxGMMcNKlqvXkncwmAXQBqVXWVmUS2Wu/sDSad4AidZix2+1irqmuDnwNos9uvNS6i8WANgsgwD0DwRBtMCPXB5GD+tV4z4tZDlsM4mQPGSbzCfB0ccdNunJwKyzaj7cPOaPslGhcmCCJDI4y/wNdb+hKCyaEGQBXM8W7C7OcNc12Yz8ci+O5jACrHsI95EeyTaNyYICitme38XrN5Z1WwLyBktWYYzUdPj7R98L05wuZjYkz5uSo44uZozHUqgv0OdvtQ1SYYA/DVIoJhrEPjIooG+yCIiMgWaxBERGSLCYKIiGwxQRARkS0mCCIissUEQUREtpggiIjIFhMEERHZ+v8SG+kVNwV+1gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig_est_mean = plt.figure()\n",
    "ax_est_mean = plt.subplot(111)\n",
    "ax_est_mean.plot(performance['est_mean'])\n",
    "ax_est_mean.grid()\n",
    "ax_est_mean.hlines(angle,xmin=0,xmax=N_EXP)\n",
    "_ = plt.xlabel('Iteration count')\n",
    "_ = plt.ylabel('Estimated $\\phi$ (mean)')\n",
    "_ = plt.title('Estimated Mean for each Iteration')\n",
    "plt.savefig('mean.png',dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Presentation plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 198,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the naive prediction of the IQPE phase, phi\n",
    "E_i = memory_list\n",
    "# Convert to integers\n",
    "theta_i = theta_list\n",
    "phi_i = []\n",
    "for i,E in enumerate(E_i):\n",
    "#     print(i,E)\n",
    "    e = [int(v) for v in E]\n",
    "\n",
    "    # Naive sum\n",
    "    sum = 0\n",
    "    for v in e:\n",
    "        sum += v\n",
    "    p_i = sum / len(E)\n",
    "\n",
    "    # Inverse naive phi approximation using expected value from experimental results\n",
    "    phi = theta_i[i][0] - (np.arccos(1 - (2 * p_i))) / M\n",
    "    phi_i.append(phi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 205,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAESCAYAAAASQMmzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de1gc530v8O/LXdzvCwgJhASSdUESINmyZQvZ4NiRZbc2cnKOG7uNa9SkSeo+PZXac5rT5DR1KqWnJz2nbSq5cV2ndmxLbn2Jr2AL2bEdC5B1vwtxkQQCsSyw3GHf88fOjFbL7rL3ZZnv53l4gJ2Zd9533tn97bzvO+8IKSWIiIgAICLUGSAiotmDQYGIiDQMCkREpGFQICIiDYMCERFp5nxQEEL0CSGa7V67KIRIdbFNkRDiop/zUebn9Pyex0CmO9v2GUqhKq8QItX+vRCEfZYJIeqEEHXB3K/N/kN+boXiuPtizgcFhdGTD2UpZYuUcrGf85DuawJCiCL1b3/mMVDpuisU+/QX22Pn7nrhXF4v7AKwTUpZHYqdB/pYu1v/4UQvQWEXgD8PdSb8YF+YpasH7h47vR7jIimlKdSZCKC5V69Syjn9A6BP+X0R1hNU/TtV+btO+Wm2eS0VQLPy9z4AZTbpXVR+71C2qVO3c7L/GgBVyvq1Lta7KT0ARTZ5q4M1sEnl7x12eUxVtt1l83uHWi6bfTgqq9N0leW1NttV2e1vj7o/J2WyT0s9dvZlc1SWaWkreVW32ePseDo4ltq6AMrU8tvs5yKAPa7SsDkWzcrPDvtjF6Jj7Oz8nOm8dlg37pyLTvLhqAw7bMpd5mQ7n+rKwfapNsd5n5Nza6b3iaM0p9WFB/V/07F28B5xmF+7+nSWhxnrxuPPTH8kMpt/cCMo1KonE2yCgt1JrS63PZFqbCqgTKm4MrUClErZ42Tfu5TlRTZvrh0O1puWnrJtld16tieL/ckubU5CCaDG5iQtc1bWGdItA7DPdj2bN4nt/vpcnPCOTvCbyjZDWfpsjss+R9vPdCxt8m77IZNqm2/l9Spnaajb2ubZ/tiF6BhPOz/dPK+d1Y1b56KDYz6tDLbpBqKunORVPW47nJTZ5ftkhjSn1YU79W9/rB28Rxzm165enL0vXNaNNz9R0Akp5V5HHcxKm2AVAIdtnlLK/UKIXQB2AvgagFeU3+lCCPXScVq7opJuqpSyRQhRpfxOB+CofdNRejsB7BFC7If15GqZoYgt8sZleguAepu/090pqwPbYS2vql7Zvt5uf0Y301PVwaZsDrZ3lrZad702f9tzVjfbYP0weUpKaVLOA9u01Q/+eidpmGBzLKSTJpFgH2Mn56c3+VB5cy46K8N+L/YFuF9XGU62N0kpdzvZ50zvk2onabp1vnt53F3l11ne1Tzc9F5y43NiRroJCor9sEZwANqIoGcBPAXrSbHdyXYtamVLKXcKIdYB2CmlrHeyPmA9MexHXKjffhyZlp4QolrJUzOANBf7AqwfWhr7Dy0PyjoTtcPc63ZiKWW9XdkW2a0yLW01qCqjWGZ6EzmqG/WDxd2OwZvSUD54XQrhMbY/P33Nh6/nIuD+wApf68q+nuwDiD2X7xMXac5YFz4cd3e/VDl6X9i/l9ypG5f00tGs+jFurqgKAK9IKQ/D9Qm4R9lOjcKv2KbjZASC9g3dRrWUcq+DdaelJ4QoktaREzsBNKlXOK6G0s7AZVmdpLsPNx+vGgCverhf23SLlH3dVDY4/9Zvm78iAPVSymop5TYXqzqrm2cBlANYZzMSzbZ+tuPGt0ZHaahXiPbp2h67UB1j+/PT3fN6Wt3Ag3PRD2Xwta6cbe9Lx7Y77+2buFv/Ltjn11G9ONv3THXjMV1dKSiXovW4cbVQD6BOibQmOI/Y9bCe+NuUdA7bjb2uA3DTN1clgtcKIWpgvRytlVI6/ObgJD2TEEL98GuxyfuHQogmWC/pPeGqrA7TVcpQZjPGeqfN5fyMlHVblHIdVn4AoMq2bHDjTaxcKexQjqcJ1gAx7Rg4OpY2V3YmIcRTAD4E8AgAKM0EasA57CwNKeVuIcQem2NRD+uxsj12uxDkY2yTtnZ+wo3z2lnduHsu2qXlsAwzZdofdeUgr46+dLnNSfldpelu/bu7f2fvGWeqXNWNN4TSWUEBovQnuGpmIjcozTcXlb6hVFi/Tb4ipZyp3dpZemqHnl7uFwhbrKvg0lvzEYUvra1Y+TbkqHmOiHykq+YjCms7ATwrhFCb4Fqc9M8QkQ/YfERERBo2HxERkYZBgYiINGHfp5CZmSkLCwu92nZoaAgJCQn+zdAsp7cy6628AMusF76Uubm5+bqUMsvRsrAPCoWFhWhqavJq24aGBlRWVvo3Q7Oc3sqst/ICLLNe+FJmIUSbs2VsPiIiIg2DAhERaRgUiIhIw6BAREQaBgUiItIwKBARkSboQ1KVGQ+rlH/XOZr+WAjRB+WpSI6WExFRYITiPoVHAe3xmOuU5wzYT2y2LRjTTT/99NNITfX5mRRhxWQy6arMeisvwDLrxQ9+8IOApBv0oGAXANQHY9tLVZ8o5CgNIUQtlAflGAwGNDQ0eJWXqakpmEw+P5MirOitzHorL8Ay64XZbPb6s88lKWVIfqAEBCfLamGdP9/hctuf8vJy6a0DBw54vW240luZ9VZeKVlmvfClzACapJPP1FB2NNdI54+n3CutD1IxKY9fJCKiIAhJUBBC1Egpdyt/V9ktq7UJBL1BzxwRkY4FPSgoQWCXEKLZ5kHfsHlQ9quwXiFUAYD08hm8RETkuVB0NNcDmPYAbilltfLbBEAdecQH3hMRBRFvXiMiIg2DAhERaRgUiIhIw6BAREQaBgUiItIwKBARkYZBgYiINAwKRESkYVAgIiINgwIREWkYFIiISMOgQEREGgYFIiLSMCgQEZGGQYGIiDQMCkREpGFQICIiDYMCERFpGBSIiEjDoEBERBoGBSIi0jAoEBGRhkGBiIg0UcHeoRAiFUCV8u86KeVOB+vUADABKJNS7g5m/oiI9CwUVwqPAkiXUu4HACFEre1CJSBASlkPwCSEqJqeBBERBYKQUoZu50LsA7BTStli89ouAK9IKQ8rAWHa1YISSGoBwGAwlL/88ste7d9sNiMxMdHr/IcjvZVZb+UFWGa98KXMmzdvbpZSVjhaFvTmI5UQogiA0TYgKFLt/s+w31ZKuRfAXgCoqKiQlZWVXuWhoaEB3m4brvRWZr2VF2CZ9SJQZQ5lR3ONlHK7g9dNANKDnRkiIgpRUBBC1KhNQg76DBpx42qhCEBdMPNGRKRnQQ8KShDYJYRoFkI027xeBwBKB3SRGiyUDmciIgqCoPcpKB/yix28Xm3zN4ehEhGFAG9eIyIiDYMCERFpGBSIiEjDoEBERBoGBSIi0jAoEBGRhkGBiIg0DApERKRhUCAiIg2DAhERaRgUiIhIw6BAREQaBgUiItIwKBARkYZBgYiINAwKRESkYVAgIiINgwIREWkYFIiISMOgQEREGgYFIiLSMCgQEZGGQYGIiDQhCQpCiCohRJ2L5X1CiGYhxK5g5ouISO9CEhSklPUzrLJNSlkupdwZlAwRERGA2dt8lCqEKAp1JoiI9EZIKUOzYyHqpJTVTpbVAngVwC4p5XYny2sBwGAwlL/88ste5cFsNiMxMdGrbcOV3sqst/ICLLNe+FLmzZs3N0spKxwulFKG5AdAnRvr7AJQ42qd8vJy6a0DBw54vW240luZ9VZeKVlmvfClzACapJPP1FnXfCSEqBVC1Cj/9oY0M0REOhOq0Uc1ACpsPvxhMxrpVQAmIUQVAEgp94cgi0REuhQVip0qH/T77V6rVn6bAKijk2YapURERH4065qPiIgodBgUiIhIw6BAREQaBgUiItIwKBARkYZBgYiINAwKRESkYVAgIiINgwIREWkYFIiISMOgQEREGgYFIiLSMCgQEZGGQYGIiDQMCkREpGFQICIiDYMCERFpGBSIiEjDoEBERBoGBSIi0jAoEBGRhkGBiIg0DApERKQJSVAQQlQJIepcLK9R1tkRzHwREeldSIKClLLe2TIhRI3NOiYhRFXQMkZEpHOzsfloHYAW5e8WAGUhzAsRka5EhToDDqTa/Z9hv4IQohZALQAYDAY0NDR4tSOz2ez1tuFKb2XWW3kBllkvAlXm2RgUTADSXa0gpdwLYC8AVFRUyMrKSq921NDQAG+3DVd6K7PeyguwzHoRqDK71XwkhEgRQvyNEOJnQohCv+fiZo24cbVQBMBphzQREfmXy6AghFij/LkNwDMA/gzANpvXvaJ0JleoncrKa3UAIKXcD6BI7WB21SlNRET+5bL5SEp5RAixFoAAUAHAKKX8iRDibl92qnzw77d7rdrm792+pE9ERN6ZsU9BSvmlEKJISvkRAChBokwIAViDxJEA55GIiILE3Y7mw0KIHwN4BdaO4H41SBAR0dzhVlCQUl4SQvwNgEet/8pnA5stIiIKBbeHpEop+wEwGBARzWGz8Y5mIiIKEQYFIiLSMCiEsdGJqVBngYjmGAaFMNXc1odVP3gfl64PhTorRDSHMCiEqeOXTZiYkjjS0RfqrBDRHMKgEKbajMMAgHPXzCHOCRHNJQwKYapDCQrnrw2GOCdENJcwKISptl5eKRCR/zEohCEpJdqNw4iMEOjoG8bIOEchEZF/MCiEoZ7BMYxNWrC+MB1SAhd7eLVARP7BoBCG1E7me27JBgCcY78CEfkJg0IYalf6EzaVZCE6UrBfgYj8hkEhDLUZhyEEUJCRgEWZCRyBRER+w6AQhjqMw8hLmYeYqAgUG5JwrptBgYj8g0EhDLX1DmFhejwAoCQ7CR3GEQyPT4Y4V0Q0FzAohKF248iNoGBIBABc6Ga/AhH5jkEhzAyPT+K6eQwLM6xBodiQBAA4z85mIvIDBoUw064MR1WvFAoz4hETGcF+BSLyCwaFMKMOR1WDQlRkBIqyEnilQER+EZKgIISoEUJUCSF2OFneJ4RoFkLsCnbeZjv1SqFAaT4CgCXZibyBjYj8IuhBQQhRAwBSynoAJiFElYPVtkkpy6WUO4Obu9mv3TiMpLgopMyL1l4rMSThct8IhsY4AomIfBOKK4V1AFqUv1sAlDlYJ1UIURS8LIWPduMwCjLiIYTQXuMIJCLyl6gQ7DPV7v8MB+ukAzAKIfZIKbfbLxRC1AKoBQCDwYCGhgavMmI2m73eNlTOdAwjPynipnz3mS0AgDcPNqIvP9rJllbhWGZf6K28AMusF4EqcyiCggnWD32npJR7AUAIYRJC1Egp9ztYvhcAKioqZGVlpVcZaWhogLfbhsKURcJY9x4eWleAyspbtNcnpyz4y8/fR1R6/k2vOxJuZfaV3soLsMx6Eagyh6L5qBE3rhaKANTZLhRC1Kr9DgB6g5mx2a5rYBTjUxYUpCfc9Lo6AomdzUTkq6AHBeVbf5Hawax0OEMIoQaHV2HTAW1/laBn9sNRbZUYkjhbKhH5LBTNR5BS7nbwWrXy2wSgXnm53n49PetwMBxVVWJIxJtHr2JobBIJsSGpViKaA3jzWhhpMw4hMkIgNyVu2rIl2cp0FxyBREQ+YFAII+3GEcxPnYeoyOnVpg5LZb8CEfmCQSGMtPcOOWw6AqwP3ImJiuADd4jIJwwKYaTdOIwFDjqZASAyQmBxViI7m4nIJwwKYWJgdAJ9wxMocBIUAGsTEu9qJiJfMCiECVfDUVUlhiRcMY3AzDmQiMhLDAphQh2O6qz5CACKs62dzexXICJvMSiEiTb14TpOOpoB65UCwKewEZH3GBTCRLtxGGnx0UiOcz7h3YL0eMRGRXBYKhF5jUEhTHQYh7EwI8HlOtoIJHY2E5GXGBTCRFvvsMtOZlWJIZF9CkTkNQaFMDA5ZcEV0wgWps+bcd1iQxI6+0cxODoRhJwR0VzDoOCmiSkLJqcsIdn3VdMopixy2pTZjmidzWxCIiIvMCi46bFnv8BTLzRBShn0fbe7MRxVpc6BxCYkIvIGg4IbBkYn0NhmxIGzPdjffDno+28zDgFwPGW2vQVp8YiLjuB0F0TkFQYFNzS39UFKICspFn/1q1PoHhwN6v7bjcOIiYyAIXn6lNn2IiIElmQnclgqEXmFQcENTa1GREUI/OvvrsPopAX/8/WTQd1/h3EY+WnzEBkh3Fq/ODuJN7ARkVcYFNzQeKkPK+anYOX8FPxxVQneO9mFd493Bm3/bb3DLu9ktldsSETXwCj6RzgCiYg8w6Awg7HJKRy5bML6wjQAwFN3LsLK+cn4/hsnYRoeD/j+pZRod/MeBVWJ8hS2C91sQiIizzAozOD45X6MT1pQUZgOAIiKjMCuR0phGh7Hj94+HfD9m4YnMDg26VlQ4BxIROQlBoUZHGo1AgAqCtK011bkpeAPNi3G/ubLOHiuJ6D7V4ejehIU8tPmYV50JEcgEZHHGBRm0NTah8VZCchIjL3p9e/cvQSLsxLw3//jeECfX9Duxuyo9tQRSOfZfEREHmJQcMFikWhqNWL9ovRpy+KiI7HrkVJc7R/B375/NmB58OZKAbB2NnNYKhF5KiRBQQhRI4SoEkLs8GZ5sJzrHsTA6CTWFU4PCgBQUZiOJzYU4t8+b0WT0szkb+29w8hMjEV8TJRH25UYknBtYIwjkIjII0EPCkKIGgCQUtYDMAkhqjxZHkyNl6wf9M6CAgD86VeWIi9lHna8dgyjE1N+z0ObccitO5nt8SlsROSNUFwprAPQovzdAqDMw+VBc6i1DznJcchPcz47aUJsFH788Cq09Azh/3103u956DCOeNx0BNwYgcTOZiLyhGdtEv6Ravd/hofLIYSoBVALAAaDAQ0NDV5lxGw2O91WSolfnxlBcVoEDh48OGNaG+dH4WcNF5E9dgUFyZFe5cfepEXiqmkEloFuj8tokRIxkcCBw6eRN9Kive6qzHNRIMo7NiVxuncKx65P4dT1KWxeGI2vFDp/Il6w6a2OAZbZn0IRFEwAnLfHzLwcUsq9APYCQEVFhaysrPQqIw0NDXC2bYdxGH3vH8ADty5D5YbCGdNas34cVX/3MV5tjcHrf3gHoiN9vwhr6TFDfnAQd5UtR2V5vsfbLzv5a4xER6Oy8lbtNVdlnov8Vd623iEcONONA2d78HlLL8YnLYiPiUR8TAw+7orAM09sghDuTUMSaHqrY4Bl9qdQBIVG3LgaKAJQ5+HyoGhqm7k/wVZqfAz+6qEV+NaLh/HsJy34duUSn/PgzXBUW8XZSfjkfGDvo5irxian0HipDwfOduPA2W609Fhnqi3KTMDv3FqAu5dlY92iNPzqaCf+ZN9RfNlhQtnCtBlSJZr9gh4UpJT7hRA71A5kpUMZQog6KWW1s+XBduhSH5LiorS2eXfcvyoX96/MwU/rz+MrK3KwOCvRpzyoQaHAiz4FwPpshdcOX0b/8ARS4mdP88Zs9s7xTrz+5RV8euE6hsanEBMVgQ1FGXj8tgJULs1GYebNDzqqXmFAzH9G4K2jVxkUaE4IxZUCpJS7HbxW7Wp5sDW2GlFRkOb2zKSqHz60Ap9euI4/e+0YXqndgAgPt7fV3juM2KgIZCXFzryyA1pnc/eg21c8enbiSj++/eJh5KXE4bfL5mPz0mxsWJzhcjhwclw0Ni/NwtvHOvEXW5Z7fL4QzTa8ec0B49A4LnSbsc7BTWszyU6Kw/cfWI7G1j78+xdtPuWj3WidCM/btupi5SlsvInNPfubLyMmMgLv/tFd+NFvrcI9txjcuj9k6+o8dA+O4dClwNyrQhRMDAoOqDeiefvtuqY8H3csycBP689jwofnOrcbh726R0GVlzIP8TGRnBjPDeOTFrxx5AqqVxg8bmq7e1k25kVH4q1jVwOUO6LgYVBwoKmtDzFRESjNT/FqeyEEnthQCOPQOD69cN2rNKSUaDcOu/VcZmciIgSK+RQ2t3x0pht9wxOo8WKUV3xMFKqWG/DeiS6fvgQQzQYMCg4cumTE6vwUxEZ5f7/BpqVZSI6LwptHvPv2eN08juHxKa87mVXFhiSc7+aVwkz2N19GVlIs7lyS6dX2W0tzYRwax2cXe/2cM6LgYlCwMzw+iRNX+n3umI2NisT9K3Px/skur6a/8HU4qqrEkIiewbGgPBAoXF03j6HhbDceXjsfUV7eX7JpaRaS4qLw1lE2IVF4Y1Cwc6TDhEmL9KqT2d6Da/IwND6FD093e7xthzY7asIMa7pWzOkuZvTGkauYtEg84kXTkSo2KhL3Ls/B+ye7MDbp/zmwiIKFQcFO46U+CAG/jDm/rSgDWUmxeOPIFY+3beu1BgVX8y6548YcSOxXcGZ/82Wszk/x6J4UR7auzsXg6CQ+PuddPxLRbMCgYKex1YhlOclImef7zV6REQIPlOai4WyPx1NYtxuHkZMch7ho3+ZRykuJQ2JsFGdLdeLk1X6c7hzwqoPZ3h1LMpEWH80mJAprDAo2JqcsONzeh3WF/rsz9aE18zE+ZcH7J7o82q7dOORzfwJgHQm1JDuRzUdOqPcmbF2d53Na0ZERuH9VLupOXcPweOCexkcUSAwKNk51DmB4fMqvd/+uzk9BQUY83vTw26N645o/FPPRnA5Z7024iqrl2UiNj/FLmltL8zAyMYWPznjej0Q0GzAo2Ghs7QPg/U1rjggh8ODqPHx28Tq6B0fd2mZ0YgrXBsZ8Ho6qKjEk4bp5HMYhjkCy1XC2G8ahcb80HanWL0pHVlIsm5AobDEo2Gi8ZMSC9HnISYnza7oPrcmDRQJvH+t0a/0OPw1HVanTXbBf4Wb7my8jMzEWdxVn+S3NyAiBLatyceBsDwZH+ShUCj8MCgopJRpbjQGZOG5JdhJuyU3GG27eyKbeo+DL3cy2bkyMx34FVa95DB+d6cbDZd7fm+DM1tV5GJ+0oO7UNb+mSxQMDAqKS9eH0Ds0HrDZRB9ak4cjHSa0K0NNXVGHo/qr+Sg3JQ5JHIF0E+3ehDL/NR2pyhamYn7qPDYhUVhiUFA0+jgJ3kzU0S1vHp35noV24zASYiKRnuCfzk8hBJYYOAeSrf3Nl1Gan4KlOb7dm+CIENahyJ+cv44+9uNQmGFQUDS29iE9IQaLs3y7g9iZ+anzsK4wDW8cuQoppct1243DWJiR4NfHO5ZkJ3G2VMWpqwM41TkQkKsE1dbVeZi0SLx30rOhyEShxqCgUB+qE8jn7D64Og/nu8040+X6G7t1OKpvdzLbKzYkondoHAPjrgOSHrx2+DKiI62jwgJlRV4yFmUm4FecTpvCDIMCgO6BUbT1DmO9H+Y7cuWrq3IRGSFc3rNgsUh0+PEeBZU6B9KVQX1P7TwxZcHrX15B1S0GpPmpec4RIQS2lubi84u9bg9FJpoNGBQQmPsTHMlIjMXGJZl400UTUvfgGMYmLViY4d9mrJV5yYiJisBLZ8Zx3Tzm17TDScPZHvT6+d4EZ7autg5Ffvc4m5AofDAowNp0NC86EsvzkgO+r4fW5OGKaQSH2/scLtemzPbzlUJGYiyefbwC14Ys+Nqez9HVr89vr/ubO6z3JpT4794EZ4oNSVhqSOIoJAorDAqwBoWyglRE+3m8uiP3rshBbFSE03sW2nqHAPhvOKqtTSVZ+JOKOHT1j+LRPZ9rN8nphXFoHB+d6cZvr80LSl0D1plTm9r6cNU0EpT9EflK90FhcHQCpzsHUFEQ2KYjVWJsFKpuMeDtY52YdPDoxg7jMCIEkJfq345m1dL0SLz41G0wDY/j0T2fo6VHPyOS3jxyBRNTvj03wVMPlFo7s929m50C64ppBDv3H0PFj+q9flTuXKf7oNDc1geLRMA7mW1tXZ2H3qFxfOrg0Y3txmHkpsxDTFTgqmbNglS8XLsB45MWPLrnNzjTNRCwfc0m+w9fxsr5yViWE/hmQlVhZgJK81PwFkchhVSveQz/661T2PyTBvznl1cQEynwzecb8evzDAz2dB8Umlr7EBkhsHZhatD2Wak8utHR85vbjMMo8NOcR64sz0vGK9tvQ2QE8PW9v8Gxy6aA7zOUTncO4MSVAdQE8N4EZx4ozcWxy/1ovT4U9H3r3cDoBP7ug7O4a/cBPP/ZJfzW2jx89N824a3vbsSizAQ8+W+N+OR8T6izOauEJCgIIWqEEFVCiB1OlvcJIZqFELsCnZdDrUaszEtGfExUoHeliYuOxH0rchw+vzkQw1GdWZKdhH3bb0dibBQee/YLNCl3dc9FrzUr9yasmR/0fW9RmpB4z0LwjE5MYe/HF3HX7gP4vx9dwKalWfjgjzdhd81q5KfFIyMxFi89dRuKshLx5L814eA5BgZV0IOCEKIGAKSU9QBMQogqB6ttk1KWSyl3BjIvExaJox2mgA9FdeShNfNhHpvEAZt5981jk7huHvfb7KjuWJgRj1e3b0BWUiy+8fNDc7KddWLKgtePXME9ywx+mzrEE/NT56GiIA2/CnC/gsUi8d6JLrx8ZlxXfUW2JqYsePGLNmz6yQE8884ZlOan4q3vbMQ/PVaOJdmJN62bnhCDl37/VizJSsRTLzSh4SyfgQEAYqYpF/y+Q+u3/1eklIeVgFAmpdxtt04NgMNSyhYnadQCqAUAg8FQ/vLLL3uVl2NXzfi7YwLfXRuLckPwrhQAwCIlnj4wguK0CHx3rXWq7o5BC77/6Qi+vToW63MDkx+z2YzExMRpr5vGLPjbxlF0DUt8Z00s1mQH93gEitlsxoXhOPz08Bj+qCwWa0NUrvq2Cfz76XH89R3zMD/Jv9/FJiwSn1+dxDuXJtA1ZH0/Rwrg7oVReHBxDJJiAneX/mwxMGjGycE4/OeFcXQPSyxJjUBNSQyWpc/8OFvzuMRPmkZxxWzB99bGojQrPM59Z+9ld2zevLlZSlnhaFkoSm/feJ/hYJ10AEYhxB4p5Xb7hVLKvQD2AkBFRYWsrKz0KiNvP/cBgAn87pY7kZEY61Uavnh48CReOtSOstvuQHJcNN4/2QV82oz77qxAaX5g+jgaGhrg7Hht2jiOx587hH84MoC///oKbCnNDUgegqmhoQFnTUnITDTiO4/cHbShqPaWD47ipSS5+U8AAA3+SURBVGc+xLXY+Xiscqlf0jSPTeKXX7TjXz5rwbWBcSzPTcb/eHAxJjvP4NBwFl5pbMdvrgHfvXsJnri9ELFRvj3vezbqGxrHB6e68A+fnkTH4BiW5SThbx5dis1Lsz2asmbjxnH8zs+/wD8cMWPPN0qxeVm2T/nqH5mAEEBSbFTAps5x9V72RUCCgvJN3l6L2mQE64e+U8qHPoQQJiFEjZRyfwCyiXN9FizOSghJQACAB9fk4fnPWvH+iS5sq1igTatdkB6YSflmkpYQgxefuhXf/NdGfPeXhzE6sTqowzcDwTwuUX/6Gh7fUBiygAAA2UlxuK0oA28d68QfV5f49EFx3TyG5z9txQuft2JgdBIbijLwk5rVuLM4E0IINPSdw4/vW4Xfu6MQz7xzGs+8cwYvfN6GnfctwwOluV7vu613CO+e6MIn53tQmp+Kb96xCFlJwX/v9JrH8P7Ja3j3RCc+u9iLKYtETrzA3399DbaW5iEiwvPypcbH4N+fvBXf+PkhbP9FM/75G2W4e5nBozSklDjcbsLzn7Xi3eOdmLRIxMdEwpAcB0NyLHKS42BIiYMhKQ45KXEwJFt/ZyfFhvTctBeQoKB+qDvRiBtXC0UA6mwXKgHFqASC6WM2/cRikbhgmsLWNcHvT1CtXZCKBenz8ObRq9agYBxGclwUUuKjQ5an5LhovPDketS+0Iw/2XcUhy4Z8cDqXNxWlOH3E3dwdAIjE1PITvLvk+5s/aZzEhNTMijTWsxk6+o8/Pl/HMfJqwNYOT/F4+07jMN49pMWvNLYgfEpC76yPAd/ULkYaxY4vqosMSTh+d9bj0/O9+Cv3z6N7/7ySzz36SX8xZZbUO7mfTkXugfx7vEuvHOiC6c7rUOXi7IS8NnFi3ju15fwtXUL8NSdRX57IJQzPYNjeP9kF9490YnftBgxZZEoyIhH7V1F2LIqFz3nDmOzj4MItMDw3BfY/otm/OyxclQtnzkwjE1O4Z3jnXj+01YcvdyPpLgoPHF7IXKS49A1MIqugVFc6x9FU1sfugfGMG53f5IQQEZCDO5dkYNnfnuVT2Xwh6A3H0kp9wshdqgdzMrVA4QQdVLKagCvAqiwWR6Qq4Tz3WYMTQR+viNX1Oc3/6zhInoGx5ThqKG5SrAVHxOFf3miAj986xTeOHIFrzR1IC0+Gvcuz8H9q3Jw++JMr+6jsFgkTl4dwMFz3fj43HU0t/dhyiKxdmEqtqzKxZbSXOSm+PemvV9fmcSKvGTckhu8exOcuW9FDr7/+gn8/YfnUX2LAXExkZgXHYn4mEjERVv/nhdz8+/oSIEzXYP454MX8atjnYgQwMNr81G7qQiLs9xrT76zOAtvfy8TrzVfxt9+cBaP/OxzbFmVix33LZ12vkkpcbpzEO+e6MS7J7pwQXlaX3lBGv5iyy34yoocLEiPR0uPGXsOtuCXh9rx4hfteGh1Hr5VuVibeNEfugdG8d7JLrxzvBOHLhlhkUBRZgK+tWkx7l+Vg+W5ydpVT8N5/zTRpMRH4xdP3orHf/4FvvViM/7psXJUOwkM3QOjePELa/mvm8ewOCsBf/VbK/Hw2vlIiHX80SqlRN/wBLr6R3FNDRjKT1Fm6N/7QAg6mv2toqJCNjU1ebzdL37Thu+/fgIf/+nmoI72sXfu2iDu/T8f44cPrsDzn7VieW4y/vGxsoDtz9N2yNGJKXx8rgfvHO9E/elumMcmkRwXherlOfjqqhxsLM502VbdMziGT8734ONzPfjk/HX0Kg+dWTk/GXcVZyE+JhLvHO/CKeVbaEVBGraU5uKrq3JhSPb8CmJscgonrvSjua0Pja19qDt1DX+5dTl+745FHqcVCH/40mGP7m6OjBCYskgkxETisdsK8M07Fs34DHFXdTw0NolnP2nBnoMtmLRY8MSGQnzn7iVo7R3Guyc68d6JLrT1Wu+qX78oHV9dlYuvrMhxWhed/SN49uNL+OWhdoxMTOHe5QZ8e/MSp1cvrvSPTOD45X4cvWzCwbM9aGwzQkpgSXYivroqF19dlYOlhiSHzV/+bl/vH5nA488dwqmr/fjH/1qGe1fkaMuOdJjw/KeX8LbSRHT30mz87h2F2LgkM6BT79vzpcxCCKcdzboNCt/75Zf4+EwnvvzB/UGtSEfu++nHiI2OxKmr/XhyYxH+7P5lAduXLyfS2OQUfn3+Ot453oW6U10YGJ1EUmwU7rklG/evysWmkixERgg0t/Xh43M9OHiuByevWj/sMxJicFdJFu4qycTGJVnT2qJbesx453gnfnWsE2e6BiEEsK4gHVtKc3H/yhxkO/lQ6jWP4XC7CU1tRjS39uHYlX6MT1ovzwsy4rEofhz/+Pv3OP3mFmwWi4RpxNpsNjI+hdGJKQyPT930v/q3+js1Phrbyhe43azoTh1fGxjF//7gLPY1X4YAYJFAVITA7Usycf/KHFQvNyDTg74249A4nv/0Ep7/zNrPcceSDHy7cgluX5zh8P01OjGFk1f7cbSjH8cum3D0cj8u2dzct9SQhPtX5WDLqly3rj4C0ek6MDqBx39+CCeu9OOnX1+DKYvEv37aiiMdJiTGRmFbRT6e2FCIwhB9w2dQcMLboHD7jz9E/rwJvPr0fQHIlWf+qeECdr93FgDw44dX4b+sXxiwffnrzTM+acFnF6/j3eNdeP9UF0zDE0iIiYQQAuaxSURFCJQVpGFTSRY2lWRheW6y2x2AF7qtAeLtY504e80aINYXpuOB0lysXZiGk1f70dTah+a2PrQoHyTRkQIr56egoiAN5QXpKC9IQ1ZSbMBGaMxmnpT5dOcAXm3qwIq8FFTfYvC5P8s8NomXvmjDs59cQs/gGFYvSMW3KxcjP20ejl22BoAjHf04d20QUxbrZ09OchxK81OwekEqVuenYtX8FI/zEah6HhidwBPPHcKX7dY7/osyE/DE7YV4pDwfiSH+ohGooDA7vj4FmXlsElnJcViWNDseOLO1NE8LCsG6m9lXMVERqFyajcql2fjR1Ep80WLEeyc7ISVwV0kWbl+cgaQ47z5glmQn4nv3FON79xTj/LVBvK1cQXz/jZPaOqnx0ShfmIaainxUFKSjND8FcdFzb8hloN2Sm4y/3LrCb+klxkah9q7FeHxDIV47fBn/fPAitv+iWVueMi8apfkpuGfZYiUIpDi9CpwNkuOi8cI31+NnDRexflE67irO8mp0UzjRZVBIjI3CG394BxoaGkKdFQDAgvR4lBekobmtL2yCgq3oyAhsLM7ExuJMv6ddbEjC04YkPF1VgnPXBnHq6gBWzk9GUWbinH9zhrO46Eg8dmsBvlaxAPWnr2Fs0oI1C1KxMD0+5M21nkqKi8aO+wLXpDvb6DIozEZPblyESYtE7gydiHpWYkhCiR9Ht1DgRUVG4L6V4X8TpJ4wKMwS1tEVfPMQUWjNntvoiIgo5BgUiIhIw6BAREQaBgUiItIwKBARkYZBgYiINAwKRESkYVAgIiJN2E+IJ4ToAdDm5eaZAObek+pd01uZ9VZegGXWC1/KXCClzHK0IOyDgi+EEE3OZgqcq/RWZr2VF2CZ9SJQZWbzERERaRgUiIhIo/egsDfUGQgBvZVZb+UFWGa9CEiZdd2nQEREN9P7lQIREdlgUCAiIo0ug4IQokYIUSWE2BHqvASLEKJPCNEshNgV6rwEmlK3dXavzdk6d1LeOVvfQohUpT5rbMs3x+vYWZn9Xs+6CwpCiBoAkFLWAzAJIapCnKVg2SalLJdS7gx1RgJNqVvNXK9z+/Iq5nJ9PwogXUq5HwCEELVzvY7hoMzK636vZ90FBQDrALQof7cAKAthXoIpVQhRFOpMhIge63zO1reUcq+UUh15UwSgHnO8jp2UGQhAPesxKKTa/Z8RklwEXzoAoxBiT6gzEgJ6rPM5X9/Kh6FRStkCndSxXZmBANSzHoOCCdYDqSvKNw0TrJfWNaHOT5Dprs51Ut81Usrtyt96qWPbMgeknqP8kUiYacSNbxVFAOpcrDsnKO2PRqU9sjfU+QkBXdW5HupbCFEjpdyt/F0FHdSxgzKrVw1+rWfdXSkoB7BI7Yhy0kk317wKm843tbNqrlK+MVXYdD7O6Tq3Ly/meH0r5dqljLppBnRRx9PKjADVM+9oJiIije6uFIiIyDkGBSIi0jAoEBGRhkGBiIg0DApERKRhUKA5S5lErNnmf5+mA7DdXghR5uskZEKIPUKIfcpvnydxm6vTWlBwcUgqzVlCiFQAH0opy5X/90kpt/mQnk/b26VVB2CXP8fT+zN/pF+8UiBdUL6JVynfylOV1/YpP7VCiCKbb+6VNn/XOtre9krBJp09yv9lQog6IcQu2ysVm7wUAY5vsHKSlvr3DmVq6GnpOyofkTcYFEgXlOkBmqSU26WUJuVD9BXlm3W5sloVgKeklA3KetsA7HS0vZquks4eZd06mymNoUxnXO9gGuci3JjRU+MqLSdl0tJ3lj8iTzEokF4tBlCtfNtXP0Tr1Q9U5du4yw9lxToATcrfh3EjwKgf+o7mpGmBNTC4m5YjrtIn8hqDAulVM4A6KeVOmweUqAGhFtarhldxI2A406isC+X3xZl2rE577OAKYqa0Fs+UNpGvGBRIV5T2+iLlgSXVat+B3WotAL4G4M+dba/+rzTbfE3pOK5WZ7GciZSyGsA2mz6EXY7SklIehnWyu11wY2po+/wReYqjj4iISMMrBSIi0jAoEBGRhkGBiIg0DApERKRhUCAiIg2DAhERaRgUiIhI8/8BpaSR+kpK5v8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the naive phi estimate\n",
    "plt.plot(phi_i)\n",
    "plt.grid()\n",
    "plt.hlines(angle,xmin=0,xmax=25)\n",
    "_ = plt.xlabel('Iteration Count')\n",
    "_ = plt.ylabel(r'$\\phi$')\n",
    "_ = plt.title('Naive $\\phi$ estimation using expectation values of experimental runs')\n",
    "plt.savefig('naive.png',dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Draw the Circuit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 221,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Disable TeX text rendering\n",
    "mpl.rc('text', usetex=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 224,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Render the IQPE circuit for M=1 and the last Theta used\n",
    "circuit = iqpe_builder.construct_circuit(k=1,omega=THETA,measurement=True)\n",
    "circ_fig = circuit.draw(output='mpl')\n",
    "circ_fig.savefig('iqpe_1.png',dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 225,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Render the IQPE circuit for M=2 and the last Theta used\n",
    "circuit = iqpe_builder.construct_circuit(k=2,omega=THETA,measurement=True)\n",
    "circ_fig = circuit.draw(output='mpl')\n",
    "circ_fig.savefig('iqpe_2.png',dpi=300)"
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}