Notebooks/analytic_error_propagation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analytical uncertainty propagation\n",
    "\n",
    "This notebook allows to calculate an analytical approximated error propagation of the hierarchical glSNe analysis. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "# matplotlib configs\n",
    "from pylab import rc\n",
    "rc('axes', linewidth=2)\n",
    "rc('xtick',labelsize=15)\n",
    "rc('ytick',labelsize=15)\n",
    "#rc('text', usetex=True)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# =============================\n",
    "# MST related error propagation\n",
    "# =============================\n",
    "\n",
    "def sigma_lambda_lens(sigma_lambda_int, sigma_kappa_ext):\n",
    "    \"\"\"\n",
    "    Total relative uncertainty in the MST per individual lens\n",
    "    \n",
    "    :param sigma_lambda_int: relative scatter of the internal MST per lens\n",
    "    :param sigma_kappa_ext: relative external convergence uncertainty per lens\n",
    "    \n",
    "    \"\"\"\n",
    "    return np.sqrt(sigma_lambda_int**2 + sigma_kappa_ext**2)\n",
    "\n",
    "def sigma_mu_lens(sigma_m_app, sigma_mu_macro, sigma_m_std_eff, sigma_m_ml_eff, num_images=4):\n",
    "    \"\"\"\n",
    "    relative magnification uncertainty per lens\n",
    "    \n",
    "    :param sigma_m_app: intrinsic scatter in the SNe apparent magnitude\n",
    "    :param sigma_mu_macro: relative macro model magnification uncertainty (per image)\n",
    "    :param sigma_F_eff: effective relative magnification uncertainties from flux measurement and unresolved mili- and micro-lensing (per image)\n",
    "    :param num_images: number of multiple images per glSNe\n",
    "    \n",
    "    \"\"\"\n",
    "    sigma_F_app = sigma_m_app * 1.086  # change from logarithmic magnitude units to relative linear units\n",
    "    sigma_F_std_eff = sigma_m_std_eff * 1.086\n",
    "    sigma_F_ml_eff = sigma_m_ml_eff * 1.086\n",
    "    return np.sqrt(sigma_F_app**2 + 1 / (1 / (sigma_mu_macro**2 + sigma_F_std_eff**2) + (num_images - 1) / (sigma_mu_macro**2 + sigma_F_ml_eff**2)))\n",
    "\n",
    "\n",
    "def sigma_lambda_int_mean(sigma_m_p_mean, sigma_mu_lens, num_lens, sigma_lum_dist=0):\n",
    "    \"\"\"\n",
    "    ucertainty in the mean of the MSR\n",
    "    \n",
    "    :param sigma_m_p_mean: uncertainty in the mean apparent magnitude of the SNe sample\n",
    "    :param sigma_mu_lens: relative magnification uncertainty per lens\n",
    "    :param num_lens: number of glSNe\n",
    "    :param sigma_lum_dist: relative luminosity distance uncertainty between apparent magnitude measurement and strong lensing sources\n",
    "    \n",
    "    \"\"\"\n",
    "    return 1/2 * np.sqrt(sigma_m_p_mean**2 + 1/num_lens * sigma_mu_lens**2 + sigma_lum_dist**2)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ======================================\n",
    "# fixed population MST error propagation\n",
    "# ======================================\n",
    "\n",
    "def sigma_H0_no_mst(sigma_ddt_model, sigma_lambda_lens, num_lens, sigma_ang_dist=0):\n",
    "    \"\"\"\n",
    "    uncertainty for fixed MST population\n",
    "    \n",
    "    :param sigma_ddt_model: relative time-delay distance uncertainty per len for given model\n",
    "    :param sigma_lambda_lens: uncertainty in the MST per lens (internal and external combined)\n",
    "    :param num_lens: number of glSNe\n",
    "    :param sigma_ang_dist: relative angular diameter distance uncertainty of the time-delay distances relative to z=0\n",
    "    \"\"\"\n",
    "    return np.sqrt(1 / num_lens * (sigma_ddt_model**2 + sigma_lambda_lens**2) + sigma_ang_dist**2)\n",
    "\n",
    "# =================\n",
    "# total uncertainty\n",
    "# =================\n",
    "\n",
    "def sigma_H0_combined(sigma_H0_no_mst, sigma_lambda_int_mean):\n",
    "    \"\"\"\n",
    "    total uncertainty\n",
    "    \n",
    "    :param sigma_H0_no_mst: relative H0 uncertainty at fixed lambda_mean\n",
    "    :param sigma_lambda_int_mean: uncertainty in the population mean of lambda\n",
    "    \n",
    "    \"\"\"\n",
    "    return np.sqrt(sigma_H0_no_mst**2 + sigma_lambda_int_mean**2)  # TODO: subtract uncertainty of population relative to sample\n",
    "\n",
    "\n",
    "def sigma_H0_scenario(num_lens, num_images, sigma_ddt_model, sigma_mu_macro, sigma_m_std_eff, sigma_m_ml_eff, sigma_m_app, sigma_m_p_mean, \n",
    "                      sigma_lambda_int, sigma_kappa_ext, sigma_ang_dist, sigma_lum_dist, no_mst=False):\n",
    "    \"\"\"\n",
    "    combined uncertainties given a forecast scenario\n",
    "    \n",
    "    :param num_lens: number of glSNe\n",
    "    :param num_images: number of multiple images per glSNe\n",
    "    :param sigma_ddt_model: relative time-delay distance uncertainty per len for given model\n",
    "    :param sigma_mu_macro: relative macro model magnification uncertainty (per image)\n",
    "    :param sigma_m_std_eff: effective relative magnification uncertainties from flux measurement and unresolved mili- and micro-lensing (per image)\n",
    "     for standardizable image (1 per lens)\n",
    "    :param sigma_m_ml_eff: effective relative magnification uncertainties from flux measurement and unresolved mili- and micro-lensing (per image)\n",
    "     for non-standardizable image (other images than the standardizable one)\n",
    "    :param sigma_m_app: intrinsic scatter in the SNe apparent magnitude\n",
    "    :param sigma_m_p_mean: uncertainty in the mean apparent magnitude of the SNe sample\n",
    "    :param sigma_lambda_int: relative scatter of the internal MST per lens\n",
    "    :param sigma_kappa_ext: relative external convergence uncertainty per lens\n",
    "    :param sigma_ang_dist: relative angular diameter distance uncertainty of the time-delay distances relative to z=0\n",
    "    :param sigma_lum_dist: relative luminosity distance uncertainty between apparent magnitude measurement and strong lensing sources\n",
    "    :param no_mst: boolean; if True, ignores error in the MST population mean\n",
    "    \"\"\"\n",
    "    sigma_lambda_lens_ = sigma_lambda_lens(sigma_lambda_int, sigma_kappa_ext)\n",
    "    sigma_mu_lens_ = sigma_mu_lens(sigma_m_app, sigma_mu_macro, sigma_m_std_eff, sigma_m_ml_eff, num_images)\n",
    "    if no_mst:\n",
    "        sigma_lambda_int_mean_ = 0\n",
    "        sigma_lambda_int_mean_corr = 0\n",
    "    else:\n",
    "        sigma_lambda_int_mean_ = sigma_lambda_int_mean(sigma_m_p_mean, sigma_mu_lens_, num_lens, sigma_lum_dist)\n",
    "        sigma_lambda_int_mean_corr = sigma_lambda_int_mean_\n",
    "    sigma_H0_no_mst_ = sigma_H0_no_mst(sigma_ddt_model, sigma_lambda_lens_, num_lens, sigma_ang_dist)\n",
    "    \n",
    "    sigma_H0_ = sigma_H0_combined(sigma_H0_no_mst_, sigma_lambda_int_mean_corr)\n",
    "    return sigma_H0_, sigma_lambda_int_mean_\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pantheon no MST (0.006140145673262882, 0)\n",
      "Pantheon 00 (0.011608144016722378, 0.00985127497454409)\n",
      "Pantheon 03 (0.014490765172216535, 0.013125581373311997)\n",
      "Pantheon 10 (0.028537599872474933, 0.027869216307470728)\n",
      "Future no MST (0.004658474953124562, 0)\n",
      "Future 00 (0.0074497656012097145, 0.005813571933336565)\n",
      "Future 03 (0.014508805192086636, 0.013740598211614777)\n",
      "Future 10 (0.02711355023750083, 0.02671035787091952)\n"
     ]
    }
   ],
   "source": [
    "\n",
    "kwargs_joint = {'num_lens': 144, \n",
    "                'num_images': 3,\n",
    "                'sigma_ddt_model': 0.04,\n",
    "                'sigma_mu_macro': 0.05,\n",
    "                'sigma_m_ml_eff': 1.0,\n",
    "                'sigma_m_app': 0.1,\n",
    "                'sigma_lambda_int': 0.03,\n",
    "                'sigma_kappa_ext': 0.025,\n",
    "                \n",
    "               }\n",
    "\n",
    "kwargs_pantheon = {'sigma_m_p_mean': 0.008,\n",
    "                   'sigma_ang_dist': 0.004, \n",
    "                   'sigma_lum_dist': 0.015\n",
    "                  }\n",
    "\n",
    "kwargs_future = {'sigma_m_p_mean': 0.006,\n",
    "                   'sigma_ang_dist': 0, \n",
    "                   'sigma_lum_dist': 0\n",
    "                  }\n",
    "\n",
    "\n",
    "kwargs_patheon_no_mst = {**kwargs_joint, **kwargs_pantheon, 'sigma_m_std_eff': 0.0, 'no_mst': True}\n",
    "kwargs_patheon_00 = {**kwargs_joint, **kwargs_pantheon, 'sigma_m_std_eff': 0.0}\n",
    "kwargs_patheon_03 = {**kwargs_joint, **kwargs_pantheon, 'sigma_m_std_eff': 0.2}\n",
    "kwargs_patheon_10 = {**kwargs_joint, **kwargs_pantheon, 'sigma_m_std_eff': 1.0}\n",
    "\n",
    "kwargs_future_no_mst = {**kwargs_joint, **kwargs_future, 'sigma_m_std_eff': 0.0, 'no_mst': True}\n",
    "kwargs_future_00 = {**kwargs_joint, **kwargs_future, 'sigma_m_std_eff': 0.0}\n",
    "kwargs_future_03 = {**kwargs_joint, **kwargs_future, 'sigma_m_std_eff': 0.3}\n",
    "kwargs_future_10 = {**kwargs_joint, **kwargs_future, 'sigma_m_std_eff': 1.0}\n",
    "\n",
    "print('Pantheon no MST', sigma_H0_scenario(**kwargs_patheon_no_mst))\n",
    "print('Pantheon 00', sigma_H0_scenario(**kwargs_patheon_00))\n",
    "print('Pantheon 03', sigma_H0_scenario(**kwargs_patheon_03))\n",
    "print('Pantheon 10', sigma_H0_scenario(**kwargs_patheon_10))\n",
    "\n",
    "print('Future no MST', sigma_H0_scenario(**kwargs_future_no_mst))\n",
    "print('Future 00', sigma_H0_scenario(**kwargs_future_00))\n",
    "print('Future 03', sigma_H0_scenario(**kwargs_future_03))\n",
    "print('Future 10', sigma_H0_scenario(**kwargs_future_10))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot forecast figure with increasing number of glSNe for Pantheon and Future scenario\n",
    "\n",
    "\n",
    "num_plot = 300\n",
    "\n",
    "num_list = []\n",
    "\n",
    "pantheon_h0_00_list = []\n",
    "future_h0_00_list = []\n",
    "\n",
    "pantheon_h0_03_list = []\n",
    "future_h0_03_list = []\n",
    "\n",
    "pantheon_h0_10_list = []\n",
    "future_h0_10_list = []\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "for n in range(1, num_plot):\n",
    "    num_list.append(n)\n",
    "    kwargs_patheon_plot = {**kwargs_joint, **kwargs_pantheon, 'sigma_m_std_eff': 0.0, 'num_lens': n}\n",
    "    h0, _ = sigma_H0_scenario(**kwargs_patheon_plot)\n",
    "    pantheon_h0_00_list.append(h0)\n",
    "    \n",
    "    kwargs_patheon_plot = {**kwargs_joint, **kwargs_pantheon, 'sigma_m_std_eff': 0.2, 'num_lens': n}\n",
    "    h0, _ = sigma_H0_scenario(**kwargs_patheon_plot)\n",
    "    pantheon_h0_03_list.append(h0)\n",
    "    \n",
    "    kwargs_patheon_plot = {**kwargs_joint, **kwargs_pantheon, 'sigma_m_std_eff': 1.0, 'num_lens': n}\n",
    "    h0, _ = sigma_H0_scenario(**kwargs_patheon_plot)\n",
    "    pantheon_h0_10_list.append(h0)\n",
    "    \n",
    "    kwargs_future_plot = {**kwargs_joint, **kwargs_future, 'sigma_m_std_eff': 0.0, 'num_lens': n}\n",
    "    h0, _ = sigma_H0_scenario(**kwargs_future_plot)\n",
    "    future_h0_00_list.append(h0)\n",
    "    \n",
    "    kwargs_future_plot = {**kwargs_joint, **kwargs_future, 'sigma_m_std_eff': 0.2, 'num_lens': n}\n",
    "    h0, _ = sigma_H0_scenario(**kwargs_future_plot)\n",
    "    future_h0_03_list.append(h0)\n",
    "    \n",
    "    kwargs_future_plot = {**kwargs_joint, **kwargs_future, 'sigma_m_std_eff': 1.0, 'num_lens': n}\n",
    "    h0, _ = sigma_H0_scenario(**kwargs_future_plot)\n",
    "    future_h0_10_list.append(h0)\n",
    "    \n",
    "    \n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOydd3iUVfb4PyeT3kPvVVp6AoRelaIiRRcbrkbErrD8FMvakIVdFVcQ/a6IiihiA3ulSREECSUBEpAaSoCQQHojmdzfHzOZ9GSGhExI7ud55pl531ve884kc+bcc+45opRCo9FoNJr6hoO9BdBoNBqNpiK0gtJoNBpNvUQrKI1Go9HUS7SC0mg0Gk29RCsojUaj0dRLtILSaDQaTb2k3igoEfEXkfUiki0iZ0RkjogYqhnjLCLzReR3EckRkUpj5kVkgojsE5FcEYkTkdtq/y40Go1GU1vUCwUlIn7AOkABE4A5wBPAy9UMdQemAdnAH1XMPxj4CtgAXA/8BHwmIqNrLLxGo9ForghSHzbqisizwFNAR6VUuvncU8BsoFXRuUrGilJKichjwFtKKamgz2rASSk1ssS5nwFvpdTg2r0bjUaj0dQG9cKCwmTVrC6jiD4H3IBhVQ1U1WhYEXEBRgBflmn6HBggIj62i6vRaDSaK42jvQUw0xP4reQJpdRJEck2t/1Qg7m7Ak7AwTLnD2BS0N2BqMoGV+XX0mg0Gk0xFa1g1YT6YkH5AakVnE8xt9V0biqYP6VMeylE5AER2VnDa2s0Go3mMqkvFlS9Qym1BFhSZEEdaWZyVbltXk6bXp3sKJmmsbBr1y569+5tbzE0mmoRqVXDyUJ9saBSgIp8QX4UWzo1mZsK5vcr024VhfnGStsKcvNIPpZAXk6eLVNqNBXSp08fe4ug0diV+qKgDmLyNVkQkfaYwsjL+o5s5SiQX3Z+83EhcMiWyapSUBf/OsX+1z/l3JEztsqo0Wg0mjLUFwX1CzBGRLxKnLsNyAE21WRipVQepv1Pk8s03QZsU0ql2TJfobFyBeXe1BuAnJQM24TUaDQaTTnqiw9qMTAd+FpEXgW6YNoD9UbJ0HMROQJsUkrdV+Lc9YAHEGo+/pu5KUopdcL8+l/ARhFZCHwL3GB+jLVVUGN+QaVtbn5eiEBuSqXbtjQaq3nppZfsLYJGY1fqhYJSSqWIyLXA25hCylOBBZiUVEkcgbLpj94BOpY4Xml+vhdYZp5/i1lxzQUeBo4Ddyql1tgsa0HlFpSDuytOLs7kaQtKUwvMnj3b3iJoNHalXigoAKVUHDCymj6drDlXydhvMVlPNaKwCgUlIhh8vbh0UVtQmprTpk0bzpzR/kxN46XeKKirhaoUFED78UPA2amOpNE0ZM6ePWtvETQau6IVlI1Up6DahHarI0k0Go2mYaMVlI0UFlQeJAGQk5pB0uHTNA/ogpu7Sx1JpWmIhIeH21uEUtz+xlpSsi5V28/Pw5nP/9+oKvsYDAaCgoIoKCigV69efPTRR7i7u9ss07fffkv37t3x9/cHYPjw4bz++ut1tods+/btzJgxg7y8PPLy8rjtttuYPXs2y5YtY+rUqURHRxMcHAxAYGAgP/74I506daoT2Wxh2bJl7Ny5k7ffftveopSivoSZXzVUFSQBkHr4FMc+/JELCcl1JJGmobJr1y57i1AKa5STtf3c3NyIjo5m//79ODs7s3jx4suS6dtvvyUuLu6yxtYG99xzD0uWLLHcy6233mppa9euHfPmzbObbA0BraBspLCgsMp2j2amhBXZSRWlFtRorOeBBx6wtwh1wpAhQzhy5AgAEydOpHfv3gQEBLBkyRJLH09PT5577jlCQkLo378/iYmJ/PHHH3z//ffMmjWL0NBQjh49CsDKlSuJiIige/fu/P777wAYjUZmzZpF3759CQ4O5t133wVAKcWsWbMIDAwkKCiIL774AoCNGzcyfPhw/va3v9GzZ0+mTJlCRYUTzp8/T+vWrQGTVVhkyQGMGzeO2NhY/vrrr3Lj1qxZw4ABAwgPD2fy5MlkZmaW67No0SL8/f0JDg7m9ttvB2DHjh0MGDCAsLAwBg4caJl72bJlTJw4kVGjRtGpUyfefvtt3njjDcLCwujfvz8XL14ETBbmjBkzCA0NJTAwkB07dpS7blJSErfccgt9+/alb9++bN26tcrP70qiFZSNGPPzq2x3b+aLCOToSD5NDXnvvffsLcIVp6CggF9++YWgoCAAli5dyq5du9i5cyeLFi3iwoULAGRlZdG/f39iYmIYOnQo7733HgMHDmT8+PHMnz+f6Ohounbtaplzx44dLFy4kJdfNtU8/eCDD/Dx8SEqKoqoqCjee+89jh8/ztdff010dDQxMTGsW7eOWbNmWYJT9uzZw8KFC4mLi+PYsWMVflHPnDmTHj16MGnSJN59911yc3MtbQ4ODjz11FP8+9//LjUmOTmZuXPnsm7dOnbv3k2fPn144403ys39yiuvsGfPHvbu3WuxMHv27Mnvv//Onj17mDNnDv/85z8t/ffv38/XX39NVFQUzz33HO7u7uzZs4cBAwbw8ccfW/plZ2cTHR3N//73P6ZOnVruujNmzGDmzJlERUXx1VdfMW3aNCs+ySuD9kHZiDGvagVl8HTHycWJ3As2JajQaBoVOTk5hIaGAiYL6r77THvvFy1axDfffAPAqVOnOHz4ME2bNsXZ2Zlx48YB0Lt3b9auXVvp3DfffLOlX3x8PGCyWPbu3cuqVasASEtL4/Dhw2zZsoU77rgDg8FAy5YtGTZsGFFRUXh7exMREUG7du0ACA0NJT4+nsGDS9c3ffHFF5kyZQpr1qzh008/5bPPPmPjxo2W9jvvvJN58+Zx/Phxy7nt27cTFxfHoEGDALh06RIDBgwodx/BwcFMmTKFiRMnMnHiRIvc99xzD4cPH0ZEyC/xg3nEiBF4eXnh5eWFj48PN910EwBBQUHs3bvX0u+OO+4AYOjQoaSnp5OaWnq1Z926daWWTdPT08nMzMTT07PS9/xKoRWUjRhzq15fFxGcfL3IaUAK6u2332bhwoUcPXqUpKQkmjVrBpiWR2bMmMHPP/+Mu7s7y5Ytq1eO/cocv0uXLmXBggWICIWFhcybN48JEyZU6PDu2LEjb775JgBxcXH06NEDg8HA2LFjeeWVV+xxWw2CIh9USTZu3Mi6devYtm0b7u7uDB8+3GKRODk5WTJmGwwGCqoIVnJxcSnXTynFW2+9xZgxY0r1/eWXX6qdp7prdu3alYcffpj777+f5s2bW6w+AEdHR5544gleffVVyzmlFKNGjeKzzz6r9NoAP/30E5s3b+aHH35g3rx57Nu3jxdeeIERI0bwzTffEB8fz/DhwyuU18HBwXLs4OBQSvaymcfLHhcWFrJ9+3ZcXV2rlK8u0Et8NmLMq94B3OXOsQTdcV0dSFM3DBo0iHXr1tGxY8dS53/55RcOHz7M4cOHWbJkCQ8//HCdyWSsIidiVZw+fZp58+axZcsW9u7dy/bt2y1RVhU5vO+9916io6OJjo6mTZs2bNiwgejo6DpRTgkJCVf8GvWJtLQ0/Pz8cHd35+DBg2zfvr3aMV5eXmRkVJ+5ZcyYMbzzzjsWi+PQoUNkZWUxZMgQvvjiC4xGI0lJSWzevJmIiAirZf7pp58svqnDhw9jMBjw9fUt1ScyMpJ169aRlJQEQP/+/dm6davF75aVlcWhQ6VzVhcWFnLq1ClGjBjBq6++SlpaGpmZmaSlpdG2bVvA9APscijys23ZsgUfHx98fEoXehg9ejRvvfWW5bjsD4m6RCsoGynMqV5BNb2mLX6tm9aBNHVDWFhYhaGx3333HXfffTciQv/+/UlNTS23ufTFF19k4cKFluPnnnvOYpHMnz/f4rQumXeuKkf5E088QUhICNu2beOZZ56xOJGffPJJq+7l/PnzeHl5WZYrPD096dy5s6WtMoe3PahvUXx+Hs612q8sY8eOtYSdP/PMM/Tv37/aMbfffjvz588nLCzMEiRREdOmTcPf35/w8HACAwN58MEHKSgoYNKkSQQHBxMSEsLIkSN57bXXaNWqldUyL1++nB49ehAaGsrf//53VqxYgcFQOhubs7Mz06dP5/z58wA0b96cZcuWcccddxAcHMyAAQM4eLB00Qaj0chdd91FUFAQYWFhTJ8+HV9fX5566imeffZZwsLCqrQiq8LV1ZWwsDAeeughPvjgg3LtixYtYufOnQQHB+Pv73/ZEZa1glJKP6p4AApQR5oNVkeaDVYbX3xPVUdWUoo6+Mt2lXYxo9q+9YE33nhDRUREqMGDB6sPPvhAHTp0SM2fP1/98ccfpfp17NhRJSUlWY5vvPFG9fvvv1uOR44cqaKiokqNOX78uAoLC1NKKWU0GlWXLl1UcnKyWr16tbr//vtVYWGhMhqN6sYbb1SbNm1SSil14cIFpZRS2dnZKiAgQCUnJyullALUF198oZRSKjk5WXXv3l0VFhYqpZRKSUkpd18ffvihevTRR0udKygoUKNHj1bt27dXkZGR6vvvv7e0vfzyy8rX11dNnDhRLV68WOXk5FR5/1ca07+nRlN7DBs2rNz/aG1Q9D2pavn7V1tQNqJyqy9GmJt4kbPfbyYp/lwdSFRzEhMT2bp1K++//z4bNmzgpptuIj09nX79+tV47k6dOtG0aVP27NnDmjVrCAsLo2nTpqxZs8ZyHB4ezsGDBzl8+DBg+gVXFE5c5CgHk1Vzyy23AODj44Orqyv33XcfX3/9tdWbPA0GA7/++iurVq2ie/fuzJw505KU9cUXX2Tnzp2MHj2aTz/9lLFjbU52r9FoahEdJGEjhdUESQB4tjQV6806X9NiwHVDkT+lR48eLF++3Opxbdu25dSpU5bj06dPW9bHSzJt2jSWLVvGuXPnLGGtSimeffZZHnzwwVJ9q3KUu7q6WpZPHB0d2bFjB+vXr2fVqlW8/fbb/Pbbb1bJLSJEREQQERHBqFGjuPfeey1KqiKHd9OmDWe5VtO4KRlheDWgLSgbUdWEmQM4+Xnj5GQgJ7lhb9YdP348H3/8MUoptm/fjo+Pj8WHU5JJkybx66+/EhUVZYmiGjNmDEuXLrVsUExISOD8+fNWO8qLHMY33HADCxYsICYmxiqZz5w5w+7duy3H0dHRluAPaxzedUnRZlKNprGiLSgbUZeqt6DEYMDJ14vcBpJNYtGiRbz22mucO3eO4OBgbrjhBt5//31uuOEGfv75Z6655hrc3d358MMPKxzv7OzMiBEj8PX1tVhAo0eP5sCBA5b9H56ennzyySeMHTuWxYsX06tXL3r06FGpozwjI4MJEyaQm5uLUqrCjY5ginT69tviKitbt27lySef5MyZM7i6utK8eXOLE3j58uXMnDkTd3d3HB0dK3R41yWNJZOERlMZUvSLUVMxIqIAjjQzbdA7Nnoko1a8XO24mP9bRdLZFK7917Ry+wwaG4WFhYSHh7Ny5Uq6ddPZ3q1FRCpMr6PR1DeKvuOUUrX6ZactKBsRKywogG53jqGni1OjV05xcXGMGzeOSZMmaeWk0WhsQisoW7HCBwXg7ud1hQW5OvD39+fYsWP2FkNTC8T7T8CYdLHafobmTegU913VfUqU2+jcuTPLly+3q7+viMqyjERGRrJ27VqOHTuGi4sLycnJ9OnTx5JKqb4RGRnJuHHj+Nvf/mZvUWqEDpKwEakmWWwR+Zk57P18Hadj46+sQJoGS1HuufqCNcrJ2n4ly200adKE//u//6upeDWmqiwjYFKqS5cutaOEjQ+toGxELlmnoAyODqRtiSYx7nj1nTWNniXvf8DNk2/nUokl5B9++MGOEtUdAwYMsKR1UlWUvxg2bBgTJkygS5cuPPPMM6xYsYKIiAiCgoIsWSR++OEH+vXrR1hYGNdddx2JiYkAzJ49m6lTpzJ8+HC6dOnCokWLyslRVZYRgH/84x8sWLCgwgwOlWVFKcJoNBIZGWm5rwULFgCmjPV9+/YlJCSEW265hezsbMBkAT388MP079+fLl26sHHjRqZOnUqvXr2IjIy0zOvp6cnMmTMJCAjg2muvtaRTKsmuXbsYNmwYvXv3ZsyYMeWyvdRntIKyEWstKAdXFxx9PMlNtO5Xp6Zxsyc6hu9//JEbx08iL8+0GbwoG3VDxmg0sn79esaPHw9QZfmLmJgYFi9ezIEDB1i+fDmHDh1ix44dTJs2zZI7bvDgwWzfvp09e/Zw++2389prr1mudfDgQVavXs2OHTt4+eWXS2UCBwgJCaFly5Z07tyZe++9t9wPhA4dOjB48OByewXXrFnD4cOH2bFjB9HR0ezatYvNmzeX6hMdHU1CQgL79+9n37593HvvvYAp83pUVBQxMTH06tWrVOqhlJQUtm3bxoIFCxg/fjwzZ84kNjaWffv2WfLjZWVl0adPH2JjYxk2bJilvEgR+fn5PP7446xatYpdu3YxdepUnnvuOds+JDuiFZSNWKugAFxbNiHvfIqOxNJYhVf/W9l9OoOx40zh8z/++KO9RbpiFJXbaNWqFYmJiYwaZSoRX1n5C4C+ffvSunVrXFxc6Nq1K6NHjwZM5SSKfEGnT59mzJgxBAUFMX/+fGJjYy3XvPHGG3FxcaFZs2a0aNHCYl0VUVWWkSKeffZZ5s+fT2FhceHSqrKiFNGlSxeOHTvG448/zq+//oq3tzdgquE0ZMgQgoKCWLFiRSl5b7rpJkSEoKAgWrZsSVBQEA4ODgQEBFju18HBgdtuuw2Au+66iy1btpS67l9//cX+/fsZNWoUoaGhzJ07l9OnT1v9OdkbraBsxJBvfYJG99bNIDWdLCsDKzSNG3FwxPP6J9ibmMeo6+uX/6m2KfJBnThxAqWUVT4oa8pJPP744zz22GPs27evXAFBa8pnFGUZefbZZ/n888/56quvSrV369aN0NBQvvzyS8u5oqwoRVnvjxw5YqlvVYSfnx8xMTEMHz6cxYsXW4oARkZG8vbbb7Nv3z5eeumlCuUtea9l77ci+UuilCIgIMAi2759+1izZk2FY+sjWkHZiC0WlHebZjg4Gsi6UH05AI0GQAyOeI6dSZx5ZTgrK8u+Al1h3N3dWbRoEf/9738pKCiocfmLkuUoPvroI5tkqSrLSEmee+45Xn/9dctxZVlRSpKcnExhYSG33HILc+fOtVwnIyOD1q1bk5+fz4oVK2ySF0x7DIuKMH766aflCir26NGDpKQktm3bBpiW/EpaafUdHWZuI7ZYUC37B9ByQCAODvXzd8CTTz3Df+e/Wn1HTZ3gPcC0VCMOBjzHTEdRyC23T+HXH76tZmTdYGjexOowc1sICwsjODiYzz77jLvuuott27YREhKCiFjKX5QtR1EZs2fPZvLkyfj5+TFy5MhSlWyrIz8/v9IsIyUJCAggPDzcomQqy4rSokULy5iEhATuvfdey9Lgf/7zHwD+9a9/0a9fP5o3b06/fv2sqm1VEg8PD3bs2MHcuXNp0aKFJaikCGdnZ1atWsX06dNJS0ujoKCAf/zjHwQEBNh0HXuhM0lUQ9lMEhf8mhBxqOo9HlcLP/30U6nCZBr7sXr1anyH34tPP1O2dmNmCsmf/IMnZzzOv14uHxWm0YBJGRZZbvbkSmWS0AqqGsoqqFQvb3of+8nq8X99uZ6LlwoZcNeoKyOgpkHw8GPT+eJALt59J1CQcYH0r14gI/GkDrDRVElDV1D1c+2pHuNoYxXLgqQUMg7GU2AsrL6zptFTkJ5M+lfPM+sxnShWUz31QTldSbSCshFHo20KyqNNM5zS0knJyK2+s6ZRU5B+nvSvnufZmY/xwvNXz14VjeZKoRWUjTgWFFBow7KLT8dWiLGQ1NPnq++sadRk7PyOF5+aybNPPwXA999/b2eJNBr7Um8UlIj4i8h6EckWkTMiMkdEqi3GIyI+IvKhiKSISJqIrBCRpmX6OIvIiyJyRERyzM8vi4hLZfOWJc3dAwDHwkIuXbLeivLq0BIHByH9Kin/rrEPU+64naUfLuPJJ/6f5Vzv3r3tKJFGY3/qRZi5iPgB64A4YALQFfgvJgX6fDXDvwS6A9OAQuBV4FtgSIk+rwAPmefaA4QDcwFfYIY1MiZ7+eCTbdqTcik7D1cXJ2uGYWjqi0eHljg714u3WlNPGTxoIIMHDSx1rm3btjpIQtOoqS/fmg8BbsDNSql0YK2IeAOzReQ187lyiMgAYDQwTCm12XwuAfhTRK5TSq0zd70TeEcpVVR2dYOItAWmYKWCcpDiL4q87Fzw87TqxkSE3k/93aq+Go1GoymmvizxXQ+sLqOIPsektIZVMy6xSDkBKKV2AMfNbUU4AWllxqYCVodEGh2KVxvzs20PeCgsLMRYqCP5NBqNxlrqi4LqCZTaKq6UOglkm9usHmfmQJlx7wMPisggEfEUkSHAw8Db1gpYWEpB5Vk7DID0w6fY8tTbHIs9adM4TePm/vvvt7cIGo1dqS9LfH6YLJqypJjbLmdclxLHz2Cyxkqm+v2fUmqOtQIaHYt1eYGNCsq9iReSe4nU42chqJNNYzWNlyVLlthbBI3GrtQXC+pKMwu4C3gc05LhdGCKiFSqoETkARHZWXSsSlpQObYt8Rma+ODm5UbGqcTqO2s0ZnQUn6axU18sqBTAp4Lzfua2qsY1r2qciDTDFLH3qFLqPXP7ZhG5BLwtIm8rpcptUlJKLQGWFKU6KiyR8LUgxzYLSkRwb9eCC2eTuFRgxNmx2uh5jaZUZm2NpjFSXyyog5TxNYlIe8Cdin1MlY4zU9I31QVTkER0mT57MCno8vn0K0AZipWKrQoKwLtDK5xT0zmf0rDLJ2g0Gk1tUV8U1C/AGBHxKnHuNiAH2FTNuFYiYimCIiJ9MCmlX8ynTpifw8uMLVo/ibdGQGUofquMl6GgmgR2xndgMB4OtZpLUdOAad26tb1F0GjsSn1Z4luMyS/0tYi8iknBzAbeKBl6LiJHgE1KqfsAlFLbRGQN8LGIPEnxRt0tRXuglFKJIvIt8KqIuAJ7gVDz/CuVUknWCNgyuzgpY0HuJZtv0L1zO8I6t7N5nKbxcubMGXuLoNHYlXphQSmlUoBrAQPwA/AysAAoWwjH0dynJLdhsrKWAh8Du4BJZfrcgynUfDrwM/Ao8C5wH1ZS4F1s3BlzbbegAAou5ZN0/KzODqCxitmzZ9tbBI3Gruh6UNVQFCSx6m//JHSjaT9wwqyHGPrUFJvnOvrFWk5siyVs7kP4ebrWrqCaBoeI6B8zmqsCXQ/KzmQ4OlteF2RdXumMJt3a43Apn8RjeulGo9FoqkMrKCvxTyrORl6YmX1Zc3h3bYuzowMX/jpVW2JpNBpNg0UrKCu51LWD5fXlKiiDjxcezXzIPHbGpppSmsbJzp07q++k0TRgtIKykjNNm1leq6ycy57Hu2s7HBOTSE6//Dk0Go2mMaAVlJXkZRRbTZJ9+cql1fBwek4bTxNPq2slahopffr0sbcI9Z7IyEhWrVoFwLRp04iLi6u077Jly65I6P7ixYv5+OOPa2UuT0/ryvhYy/fff88rr7xSq3PWJfVlH1S9xyelOCdtTRSUW9sWtG1bGxJpNA0HpRRKKRwcLv838/vvv19l+7JlywgMDKRNmzaXfY2KeOihh6zuW1BQgKNj3X3tjh8/nvHjx9fZ9WobbUFZiUOJVEcONiaLLUvKoVPE/PwnBUZdH0pTPxCRK/aojPj4eHr06MHdd99NYGAgp06dYs2aNQwYMIDw8HAmT55MZqZpg/ycOXPo27cvgYGBPPDAAxWG3w8fPpydO3diNBqJjIwkMDCQoKAgFixYwKpVq9i5cydTpkwhNDSUnJwcoqKiGDhwICEhIURERJCRkVFqvo0bNzJs2DAmTJhAly5deOaZZ1ixYgUREREEBQVx9OhRwLRf7fXXXwfgyJEjXHfddYSEhBAeHs7Ro0fZuHEjQ4YMYfz48fj7+wPwxhtvEBgYSGBgIAsXLqzw/Zk/fz59+/YlODiYl156yfKe9erVi/vvv5+AgABGjx5NTo7pB/OiRYvw9/cnODiY22+/HTAp5cceewwwWZvTp09n4MCBdOnSxWJ5FhYW8sgjj9CzZ09GjRrFDTfcYGmzO0W/XPSj4gegAPXj1HnqSLPB6kizwWpD37+rmnBi5Xq1/vE31OnzaTWaR9Oweemll+rsWkV/51fiURnHjx9XIqK2bdumlFIqKSlJDRkyRGVmZiqllHrllVfUyy+/rJRS6sKFC5Zxd911l/r++++VUkrdc889auXKlUoppYYNG6aioqLUzp071XXXXWfpn5KSUqpdKaXy8vJU586d1Y4dO5RSSqWlpan8/PxS8m3YsEH5+PioM2fOqNzcXNWmTRv14osvKqWUWrhwoZoxY4ZSyvQ5zZ8/XymlVEREhPr666+VUkrl5OSorKwstWHDBuXu7q6OHTumlFJq586dKjAwUGVmZqqMjAzl7++vdu/erZRSysPDQyml1OrVq9X999+vCgsLldFoVDfeeKPatGmTOn78uDIYDGrPnj1KKaUmT56sli9frpRSqnXr1io3N7fUPX/44Yfq0UcftbxXf/vb35TRaFSxsbGqa9euSimlVq5cqa6//nplNBrV2bNnla+vr+U9tZYSn3Wtfv9qC8pKDCUsKMNlZpIoollgZwzGAs4d1AUMNZXTGDJJdOzYkf79+wOwfft24uLiGDRoEKGhoXz00UecOGFKpblhwwb69etHUFAQv/32G7GxsZXO2aVLF44dO8bjjz/Or7/+ire3d7k+f/31F61bt6Zv374AeHt7V7j01rdvX1q3bo2Liwtdu3Zl9OjRAAQFBREfH1+qb0ZGBgkJCUyaZEpk4+rqiru7OwARERF07twZgC1btjBp0iQ8PDzw9PTk5ptv5vfffy8115o1a1izZg1hYWGEh4dz8OBBDh8+DEDnzp0JDQ0FTCVZiuQIDg5mypQpfPLJJ5UuI06cOBEHBwf8/f1JTEy0yDN58mQcHBxo1aoVI0aMqPS9rWu0D8pKHEoki3XKq5mCcuvcFk93Fy7GHUcNDqhyGUTTeGnTpk2d5eNTdjqBrmMAACAASURBVNr24OHhUUqGUaNG8dlnn5Xqk5ubyyOPPMLOnTtp3749s2fPJje38mV2Pz8/YmJiWL16NYsXL+bLL79k6dKllyWfi0txMJODg4Pl2MHBgYKCAqvnKXmf1qCU4tlnn+XBBx8sdT4+Pr6UTAaDwbLE99NPP7F582Z++OEH5s2bx759+6q8H3t95ragLSgrcShRUdfpUs0UlLg44921LerEGdKybU88q2kcnD171t4i1Cn9+/dn69atHDlyBICsrCwOHTpkUUbNmjUjMzOzWv9IcnIyhYWF3HLLLcydO9dSV8vLy8viZ+rRowdnz54lKioKMFk/tiicivDy8qJdu3Z8++23AOTl5ZGdXX7P5JAhQ/j222/Jzs4mKyuLb775hiFDhpTqM2bMGJYuXWrxwSUkJHD+fLmydRYKCws5deoUI0aM4NVXXyUtLc0ytjoGDRrEV199RWFhIYmJiWzcuNHKO77yaAvKSgwlootc8i5hLFQYalA6o3lAFxIOnyEtNQtfDx1yrtE0b96cZcuWcccdd5BnXqWYO3cu3bt35/777ycwMJBWrVpZluUqIyEhgXvvvZfCQlMQ0n/+8x/AFCTw0EMP4ebmxrZt2/jiiy94/PHHycnJwc3NjXXr1tU4zHv58uU8+OCDvPjiizg5ObFy5cpyfcLDw4mMjCQiIgIwhceHhYWV6jN69GgOHDjAgAEDAFP4+SeffFLK1VASo9HIXXfdRVpaGkoppk+fjq+vr1Uy33LLLaxfvx5/f3/at29PeHg4Pj4V1Y+te3Sy2GooShb724w36LDia8BU06P1qfV4ujpXNbRKVEEBRnHA0aCNWE3F9O7dm127dtlbDE0jIDMzE09PTy5cuEBERARbt26lVatWVo+/UslitQVlJUHt/Eh2dMSpoAAHICclC8/Wl6+gxNHR8uYrpbQfSlMOrZw0dcW4ceNITU3l0qVLvPDCCzYppyuJtqCqociCMmZksT/gZjyyTSXbHTevoGOvDlWOrY6M2KPEfL6eltNuplvHZtUP0DQqHnjgAZYsWWJvMTSaatHlNuzM7sQs8tyKazjlplrngKwKV28PHNIyOBt9uMZzaRoe7733nr1F0GjsilZQVnJuzXYulQjRzEvNqKK3dTi2a4lPcx+y4o6Re6lmEUQajUbT0NAKykoGZiQhLsU+p7y0mltQIkLz0G64nE0k/kxq9QM0Go2mEaEVlJX4zbqXQnc3y3F+elbtzBvUDXcHSIjRy3ya0iQkJNhbBI3GrmgFZSVbDpzjkqe75bgg4/KKFpbFsUMrWg4KpkvX1lfFzm5N3aGj+EqTmprK//73P3uLYRW1WeaiU6dOJCcn18pcYCqEOX369Fqb70qio/iqoSiK74f5X+D85x66btwCwIlHIxk5+z67yqZp2IiI/tFSgvj4eMaNG8f+/fvLtdV1GYvaxJIYtZJSI506dWLnzp00a1Z/I311FJ+d8Tx2AocSQRIFtbTEV0TG6fMc3H1EfyFp7MayZcvKPYpSAeXn51fYHh0dDUB2dna5Nmv45JNPiIiIIDQ0lAcffBCj0ciJEyfo1q2bJWXRkCFDWLNmDc888wxHjx4lNDSUWbNmVVjGoqL5wJSJYdasWQQEBHDdddexY8cOhg8fTpcuXfj+++8BUzaGWbNmWUpcvPvuu+XkjY+Pp2fPnkRGRtK9e3emTJnCunXrGDRoEN26dWPHjh2W97KozEViYiKTJk0iJCSEkJAQ/vjjjwpLjXz22WcEBQURGBjI008/bfX7VXR/zz33HCEhIfTv39+SCHblypUEBgYSEhLC0KFDAVMZkXHjxgGmhMRTp061vBeLFi2yXOtf//oXPXr0YPDgwdxxxx2WkiJ1Sa0pKDHRWUSCzc8Nbuepg0dxmLmxFhWUUopzX6zhxHebSUqvWa0pjeZq4cCBA3zxxRds3bqV6OhoDAYDK1asoGPHjjz99NM8/PDD/Pe//8Xf35/Ro0fzyiuv0LVrV6Kjo5k/fz4Au3fv5s033+TQoUOVzgemvH4jR44kNjYWLy8vnn/+edauXcs333zDiy++CMAHH3yAj48PUVFRREVF8d5773H8+PFych85coQnnniCgwcPcvDgQT799FO2bNnC66+/zr///e9y/adPn86wYcOIiYlh9+7dBAQEAHD48GEeeeQRYmNjcXJy4umnn+a3334jOjqaqKgoS06/6t6vovvr378/MTExDB061LJFYc6cOaxevZqYmBiLIi7LwYMHWb16NTt27ODll18mPz+fqKgovvrqK2JiYvjll1/YuXPn5XzENabGNrGIOAOvAZFAybz26SLyIfCMUqpm2VXrCQ7eJTIvp9U8zLwIEaFVvwDOfr6eI4cSaNH3mlqbW3P1UtEv+CtJZGRkpW1OTk5Vtru7u1fZXhHr169n165dltx6OTk5tGjRAjDlp1u5ciWLFy+2WGkVUbKMRVXzOTs7M3bsWMBUKsPFxQUnJ6dSZTPWrFnD3r17Lclo09LSOHz4sGX+Ijp37kxQUBAAAQEBXHvttYhIhSU4AH777TdLSXiDwYCPjw8pKSmlSo1ERUUxfPhwmjdvDsCUKVPYvHkzEydOtOr9cnZ2tlhFvXv3Zu3atYApEWxkZCS33norN998c4Xv4Y033oiLiwsuLi60aNGCxMREtm7dyoQJE3B1dcXV1ZWbbrqp0s/gSlIbi7b/BwQAtwK7gVTAFwgHXgLeAh6ohevYHUcfL8trh/Sah5mXxD24G74/buHEzgPkhXbGxanipJCaxsMDDzSIf5tKUUpxzz33WJK5liQ7O5vTp08DpjxxXl5e5fpA+XIdlc3n5ORk8ZNUVjZDKcVbb73FmDFjqpTbniU4rLk/g8FgkWPx4sX8+eef/PTTT5XmdixbvqOmWd1rk9pY4rsFuEkptUYplayUKjA/rwEmApNr4Rr1Ahef4kzHjhm164MyeHvSvGdH3I6f5Hhieq3Orbk6aYCr5KW49tprWbVqlaWMxMWLFy0FCp9++mmmTJnCnDlzuP/++4HS5TJsnc8axowZwzvvvEN+fj4Ahw4dIiur5v/n1157Le+88w5g8nOlpaWV6xMREcGmTZtITk7GaDTy2WefMWzYsHLz2Hp/R48epV+/fsyZM4fmzZtz6tQpq2QeNGgQP/zwA7m5uWRmZvLjjz9aNa62qQ0FpQCnStocze0NAiff4l88Tlm1a0EBNOkXiBeFZJ6pvZBSjaa+4u/vz9y5cxk9ejTBwcGMGjWKs2fPsmnTJqKioixKytnZmQ8//JCmTZsyaNAgAgMDmTVrltXzWcu0adPw9/cnPDycwMBAHnzwwVqxJt588002bNhAUFAQvXv3Ji4urlyf1q1b88orrzBixAhCQkLo3bs3EyZMqPH9zZo1yxJ4MXDgQEJCQqySuW/fvowfP57g4GCuv/56goKC7FKCo8Zh5iLyJjAMmAfEYFri8wFCgWeBzUqpf9RQTrtRFGa+4eHXaD22P473/ROAZB9f+h35oVavpYxGCvLycXJ3rb6zpsGjw8w19qSoBEd2djZDhw5lyZIlhIeHV9i3Ppfb+H/Ac8DrQHtMFpMAp4D3gfJhLVchDg6Ci5cHRvOxe24OefnGWvUVicGAk7sBpRTZOZfwcNeFDBszRU5vjcYePPDAA8TFxZGbm8s999xTqXK6ktTqRl0R8QG8gAylVPmF1quQIgvqwqsf4DFxJCcH/d2yLur11y+0aFKzCpxlUUYjxxd/RWyhM0PvuwEf98uvOaXRaDR1wVWxUVcplaaUOt1QlFNJ/Gbdi0v3TuS5FefjyziXUuvXEYMB35Z+eB47wYH4pFqfX3P1YK/QXo2mvmCVghKRFBFZLyKvi8idItKztgUREX/zNbJF5IyIzBGRatfPRMRHRD40y5gmIitEpGkF/ZqKyLsick5EckTkoIjcba18Ww6e41xqNnkexfn4Ms9ftPr+bMGrXxBNnYXE7ft1GY5GjL0ipzSa+oK1PqhEYDgwAnNUnohkA9GY9j4VPWKVUoW2CiEifsA6IA6YAHQF/otJgT5fzfAvge7ANKAQeBX4FhhSYn5vYDOQCTwOJAP+gNXrZ5fW/kHu8DAKPDwg+QIAOedr34ICU52olj06cP7AEeJOhhF+TYsrch2NRqOpz1iloJRSPc3+pQigHzAUuA4YZH4UObLyRGQfsEsp9YgNcjwEuAE3K6XSgbVmpTJbRF4znyuHiAwARgPDlFKbzecSgD9F5Dql1Dpz138CLkAfpVSO+dwGG+Sjt2M+bk5wxrvY55SXdGUUlIjgNzScZgdOcCrqACFdmmGoJJGkRqPRNFSs/tYz+5fWKqXmAkXK52FgLPAM8ClwBAgDHrRRjuuB1WUU0eeYlNawiodYxiUWKSeznDuA4+a2Iu4FPiihnGzG95HbcenVBVUim8SlC1fO1ebUtT2dpozl2gn9tXJqpOgQ89rhhhtuIDW15gVBSyaArS2mTZtW4b4ojYnL/eYr+s9JMmeQmK+U+rtSKhjwBHrbOF9P4GCpCyh1Esg2t1k9zsyBonEi0hloAaSKyM8icklEkkTkDXMeQatYE3OK44npiG9xusGCi1dOQYkIvqE9cHd3QSmFsVB/WTU2lixZYm8RGgQ///wzvr6+VvWt6zQ/77//viUTu6Y8tV5ARSl1CZNvyhb8MG3wLUuKue1yxnUxv25lfn4Nk1U2FgjBtD+rAHiqoolF5AFK5BA0fr+BrCHBODYp3k1dmHLlUxLlxB5l9+ooXMePpLf2RTUqHnzwwTrNx2dNiYzu3bszcOBAS//Q0FBCQ0PJzs7myy+/LNW3uuSx8fHxXH/99QwePJg//viDtm3b8t133+Hm5kZ0dDQPPfQQ2dnZdO3alaVLl+LnV/qrIDIyEjc3N/bs2cP58+dZunQpH3/8Mdu2baNfv36W+ylZT+njjz/m9ddfR0QIDg5m+fLlREZG4urqyp49exg0aBDPP/88U6dO5dixY7i7u7NkyRKCg4NLXTspKYmHHnqIkydPArBw4UIGDRrE7NmzOXnyJMeOHePkyZP84x//YPr06WRlZXHrrbdy+vRpjEYjL7zwArfddhvDhw/n9ddfp0+fPnh6ejJjxgx+/PFH3Nzc+O6772jZsiVHjx5lypQpZGVlMWHCBBYuXEhmZu1nsqmPNIa1o6K4/Fil1P1Kqd+UUguA/wDTRcS9okFKqSVKqT5Fx25nE1GpGbi1KP4nUVfQgirCQSk8z53nxJ9xZObmX/HraTR1yeHDh3n00UeJjY3F19eXr776CoC7776bV199lb179xIUFMTLL79c4fiUlBS2bdvGggULGD9+PDNnziQ2NpZ9+/aVy4IeGxvL3Llz+e2334iJieHNN9+0tJ0+fZo//viDN954g5deeomwsDD27t3Lv//9b+6+u3yw74wZM5g5c6alLMW0adMsbRWVr/j1119p06YNMTEx7N+/35JZvSSVlcyYMWMGM2bMYN++fbRr1872N/kqpr6UoEzBlB6pLH7mtqrGNa9mXNFz2aCI34CXMUUM7qtWQjFtovVoW3w5x9Qrr6Cc/bvQumtrkmMPsie0G0MC2l7xa2oaJ7aWyyjZ/3LKbYCpdEVoaChgKhMRHx9PWloaqamplmSp99xzD5MnV5xz+qabbrKUumjZsmWpMhjx8fGWucFU9mLy5MmWyrRNmjSxtE2ePBmDwbSrZcuWLRZFOXLkSC5cuEB6eunVknXr1pXyHaWnp1usmorKVwQFBfHEE0/w9NNPM27cOIYMGUJZKiuZsW3bNkttqDvvvJMnn3yy+je2gWDtPqifReQFERljDgmvbQ5SxtckIu0Bdyr2MVU6zkxJ39RR4BLFlpTlEuZnq8LiHRwEY74Rnw6tLOdc60BBiYMDvsP70FoVcD7qABcydEHDxkJlBeYaEjUt9VCy1EXZMhhXsvRFYWEh27dvJzo6mujoaBISEvD09CwlExTfU/fu3dm9ezdBQUE8//zzzJkzp9yclZXMaMxYu8Q3FpO18TOmPUQbMAVKjBWRkSLSrIZy/AKMEZGSRV9uA3KATdWMayUig4tOiEgfTP6nX8DiE1uLaQ9XSa7FFIRxxBoBxWDAWGDEu2NLyznPjExy841VjKodnLp3olWP9jQ78Bdxx3V2icZC7962xho1DHx8fPDz8+P3338HYPny5eVKT1wOI0eOZOXKlVy4YNrHePFixRvthwwZYqlUu3HjRpo1a4a3t3epPqNHj+att96yHFdVVBHgzJkzuLu7c9dddzFr1ix2795ttdz9+/e3WHSff/651eMaAtYu8flgKkAYjilCLxyTgpoG3AcgIucwZTOPAWKUUra8k4uB6cDXIvIqJgUzG3ijZOi5iBwBNiml7gNQSm0TkTXAxyLyJMUbdbeU2AMFMAfYYq7w+xkQjCk0/l/WVvt1cHaiMC8fx5bFSSq8c7K4kJZN22YVF1OrLUQE7zED6RWeRhO9xNdoaNu2baMNNf/oo48sQRJdunThww8/rPGcAQEBPPfccwwbNgyDwUBYWFiFgSGzZ89m6tSpBAcH4+7uzkcffVSuz6JFi3j00UcJDg6moKCAoUOHsnjx4kqvvW/fPmbNmoWDgwNOTk6W+lDWsHDhQu666y7mzZvH2LFj7VL2wl5cdrJYEfHAVFKjSGH1xrS0ZgCUUsqmNN8i4g+8DQzAFJn3PjBbKWUs0Sce2KiUiixxzhdYAEzCZBH+CExXSpUqqiQiYzAFRgQA54ElwLzqMl8UJYvd+txipE0LBjwyidgOo3HNMW2pyl77EUGhXaqaotbJLzCiAGdHXXW3IaPLbWjAVF3Yzc0NEeHzzz/ns88+47vvvrO3WKWod+U2lFJZwFbzAwARccMUwm1zXnalVBwwspo+nSo4l4ppI+691YxdDay2Va4iHJydMF4yRdHl+vpYFFRG/DmoQwWVtWU3f+46hvv1Q+jfvWX1AzQazVXNrl27eOyxx1BK4evry9KlS+0tUp1hlYISkVOY8tt9i8mCqdDxYs7UsN38aFC4uLtQYP41W9DEF86eAyAr4XydyuEgQouzZzmw+xCdWnjRyrfCKHlNA6Co1LmmcTNkyBBiYmLsLYZdsDZI4jtgIqZgg/MislxEJlW2h6gh0q1bG3p1NoeYNy8OT807U7dBC679gmjTuRUt98ayPS6BfKPNuXk1Vwk6k4SmsWOVglJKPaaUag/0x+S76QN8BSSJyHciEllRiYuGhOeEEXhONK1AGloU36ox8UKdyiGOjnhdP5hOroLDrliijydXP0hzVdJYo/g0miJsyiShlNqhlHpWKdULU7mKuZhSCX0AnBORDSIyXUQ6XAFZ7cqhM6msjTkNgEuJzboq6crUhKoKpw6tadqnF53PnSE1OU3n6Wug2BKKrNE0RC471ZFS6qBS6j9KqX5AB2AmYAReB46LyG4RKZ/P4ypF4o7guXEbSqlSm3UNF65MyY3qcBvZj+6PTebaftdgcKjVwBmNRqOpF9RKLj6lVIJS6m2l1HVAS0wRdfFAYG3MXx9o7+tGQEtPRIQm3YvzYXlduMilgiu/WbcsDq4uuLRujoMIWWeSiI5P1iHJDYzWrVvbWwSNxq5c9j6oxkLRPqiS71PB+YucCJgAQLazC032fkfbpralSqktLh08zskVv7AnIIheAwMI7NCk+kEajUZTi1ypfVCXZUGJSHdzFN+DIvKA+XW32hSsvnEqOZOv/zxORk4+huZ+5Ds5AeB+KY/EU/ZLP+R0TXuadWpF1wMH2HvwDOfTLrsmo6aeMXv2bHuLoNHYFasVlIj0EpE3zSXVDwCrgHcwpSlaBRwUkTMistCcFaJBUXjsFF7friYvNQMRIadZcSTfxb9O2E0ucXTEc9ww2ns60XbXHjbtO02WLsvRIKisxIRG01ioVkGJSFcRWQXsx5R3by+mxLF3AzcAN5pfz8GUh28asE9EVopI3eYAuoI4KoVTWgYFOabUfcbWxZF8WccS7CUWAI4tm+I9ZhDdCnNx33uAqCM6oaxGo7n6sSaTRBymekmRwNfmFEeVYs7R9zdghnmsaw1lrBcYXE1LegW5lwBwaNcKdpvKSOWfPGs3uYpwCe6Ob1IK4Z6e+PXQKZA0Gs3VjzUKarJSyurCNGYF9hHwkYhMuGzJ6hlOLs4A5JsVlFvnElnFExLtIVI53K/tR1Fqj4ICI2dSsunQ/MpmWtdcOXbu3GlvETQau1LtEp8tyqmCsfUr5W4NcHI1KagiC6pJj+K9yM6JSfUqxPvSweMcWfQ5v++J52CCffZpaTQaTU2plX1QjQFHs4Iy5pkUlE/39pa2pikpJKXXn0q3Dl7uNFf5dN8Tw86DZzmWmF79IE29o0+fPvYWQaOxK7WmoMREZxEJNj83qPQGTu4mV1qBWUE5dy22oJqnpXDqXKpd5KoIx7Yt8Ro3jC5yic779vPHgbOcvpBpb7E0Go3GJmqsoETEWUQWAinAUSDa/JwiIgtExKWm16gPOLo4IQLGPFMIt4OnO9lN/AAwFBaSuO+YPcUrh3OPzniMiKB7bjrt/vqL7X8lUqAzn2s0mquI2rCg/g+IAG4FWgDO5udbzeffqoVr2B8nR7zcnHGjREaJDm0sr9Pj6peCAnDpE4BH3yBCu7VieGAbHA16Rfdq4qWXXrK3CBqNXbnsiroluAXoppQqWXciGVgjInuAQ8ADtXAduyIODlwTfg3OHYv3Pzl37wjRsQAYj5y0l2iVIiK4jYzAzfy6MDuXAxey8XZzpn0zT3uLp6kGnUlC09ipjZ/UCnCqpM3R3N4g8L7zRlxDexYf+xfvQ3Y6daZeRfIVISKICMbUDFI/+IrEjbvZHHeW+PMZ9hZNUw1t2rSpvpNG04CpDQX1CfCriEw25+hrISLdRGQy8BPwcS1co16w/VCipSYUQNOQ4vSDTZOS6lUkX1kcvD1w7tCGoIR42p2I5/cDZ4k7rUPQ6zNnz9p/A7hGY09qY4nv/wHPYaoD1R6TxSTAKeB94N+1cI16gd+fu/HML4AQU7kN5+6dLG2tLyRz5PRFWvi0rWS0fREHBzzGDQUHISDuKC6FRnYpRc6lAnp3aV79BBqNRlPH1FhBKaWMmPLwzRERH8ALyFBKpdV07vpGhy6tSh07tmhCno83LmnpuBTkczr6CATUTwUFIAYDHuOGIS5OdNtzEGc/b7zdWlU/UGMXwsPD7S2CRmNXrFJQIpIC7Ab2mJ93K6UOlu1nVkoNTjEV4TYoDGOhQillqX9S2KMz7IgBIDP6L5gyzJ4iVouI4D5qII4tmxHm3xVxMv0JJFzIwsfDGU/XytyJmrpm165d9hZBo7Er1vqgEoHhmJbzlgOxIpIhIr+bS3DcIyJBItKg45iPnkvn082HyCxRzsIjpIfltRw8ag+xbEZEcAnpgTg5UpibR/qPG/kz5gQ/7z7JudRse4unMfPAA1d98KtGUyOsUihKqZ5AE2AM8BKwHvAABgGPA0sxbdDNFJE/ReR/V0Zc++J44DDtvvqR3Kw8y7mmEQGW180Trr6CgcbkVIyHTjIwNgb3rEzW7U1g34kLFNbDiMTGxnvvvWdvETQau2K1xaOUSlNKrVVKzQUeMZ9+GBgLPAN8ChwBwoAHa1vQ+oCzqwtiLCQvvbjiiFtwd8vrdkmJxJ28aA/RLhundi3xuvMGXA1C/73RdM1MITr+Auv3JmAs1EpKo9HYj8tdkiv65kpSSq1RSs1XSv1dKRUMeAK9a0e8+oWLtxsAeRnFCsqxUxsK3E3nvXJzOLbnkF1kqwmOrZrh/febcG7ZlO779tE/L5Umni4YHBpUOkWNRnOVUes+I6XUJaVUdG3PWx9w9TJVW8rLLF7GEwcHVHCxHyrrz/11Lldt4ODlgdcd1+PWP4SuAwLo3dUUep6UnkPUkfM6j58dSEiwb6VmjcbeNOightrG2csDEbiUUTqQwG9QqOW158EjZOXllx16VSAGA+7D+uDg5YFSiqyfN3N+3zEOJqTy464TV51/7WpHR/FpGjtaQdmAwd0VR4MD+Zmlv6i9+gdbXnc+c5p9J64uP1RFqKwcCs5doPXWPxlx8QyF+QWsiT6lrak6ZPz48fYWQaOxK1YpKBH5WUReEJExIuJ3JQQREX8RWS8i2SJyRkTmiIjBinE+IvKhiKSISJqIrBCRplX0nyAiSkRsr6ft7ITB0VBOQbn29keZ90W1vZBETOzpikZfVTh4uuN990249g7A/fAxhsfG0MPJyMGEVI7qAogajaYOsNaCGgu8DPyMKVP5BkyBEmNFZKSINKuJEGalt8485wRMmSmeMF+zOr7EtEdrGhAJ9AW+reQ6rsACTPu6LkdODB5u5GeWXuJz8PLA2K2T6bVSXNxku+6rj4ijI+7X9sPr1rE4GAvotm8/o4Pbck0rHwDOpWaTkXN1LmdqNJr6j7WpjnyAcPOjt/lZYVIK9wGIyDkgpuihlPrcBjkeAtyAm5VS6cBaEfEGZovIa+Zz5RCRAcBoYJhSarP5XALwp4hcp5RaV2bILCABU0HFQBvks+DXzBsP5/J63XdEXzIPHQegaexfnE/LoYWP2+Vcot7h1KkN3pGTKEzPxNfPA1VQQP6ZZP48k0NWXgGBHZrg385P15uqZd599117i6DR2BW53BIRIuIBhFKssHoDPQEDoJRS1S7PlZhrM3BGKXV7iXMdgBPAeKXUD5WMmwM8oJRqVeb8MeAbpdQTZeaLBYYB04FApVQfK2RTmG4IgJytexBnJ1z7ltZvWWu3ce7OpwA43aw56csXcFOfjtVNf1WS8+decjbtRPl3ZX+bjpzMuIS7iyOhnZrSuaU3DqLD0zWaxkRR6jelVK3+8192slilVBaw1fwAQETcgBBMCssWegK/lZn/pIhkm9sqVFDmtnI5AYED5raS/Bf4Uim1W2rwBeo2KIxLBUaMhYUYHIotBrcBghSWCAAAIABJREFUISiDATEaaZecxIpdRxqsgnIN64XKziV3ZywhJ85wTWgAMU4+/PFXIm7OjrRp4mFvERsEIlIva4xpNHVFra7JKKVylFLblVK2pjryA1IrOJ9ibqvROBEZiWkp8J/WCiQiD1QUSJGUnsMXmw9x5mJWqfMOnu44hvWyHBdu3nnVhptXhzg74T4iAu+/34SDtyduW3cyNPk01wa1pbWfaa/YkbNpJOq8fhqNpgZUq6BE5NrLnVxErrvcsbWFiDgCi4B5SimrgyOUUksqWgJ0O5lAxPrf8LpUvjihz9hBltcBxw7z56Hzlyf0VYJjq2Z43TUOjxuH4hrakzZNPFDZueSnprP/1EXWxJxmbcxpnYBWo9FcFtZYUL+KyG8iMs7KsG8nEZkkIpswRf1ZQwqmQIyy+JnbajLufnOfZSLiKyK+gDNgMB/bVF/CrXUzOtwwEG/v8stYHmOKFVSvE8fZHHPSlqmvSkQEl4BrcGxpiuzP2bKbjA++5trUs/Rp7Ulqdh5rY06zJuYUFzPrb8Xh+si4cePsLYJGY1es8UGFAW8A3wNJIrIO2IEpEu4ipuq5TYBuQH/gWsAXWIMpiMIaDlLGZyQi7QF3KvYxlRw3pILzPSkONe8BtKPi0PIU4O+YytZbhaGpL4W9A0k3FuJbps2pRyekfWvUqbO45V8iY9NO0m/ujbe7s7XTX/W4DTR95Hl7/6LV/iO0C+/F6XbtiU3KRjD5/vLyjTg5Ouhgimr44YfKXK8aTeOgWgtKKbVfKTUaU2mNNcBNmPYSfQf8DmzGpAxex+Tn+Q7or5S6XikVZ6UcvwBjRMSrxLnbgBxgUzXjWonI4KITItIH6GJuA3gbGFHmsRo4ZH691koZLWyLPsH26Phy50UE73FDLcehhw6wMfaMrdNf1Th4eeAxZhA+992CU7cO5Eftp8OxY0zq1xk/TxcAdhw5z7c74jlwOoVLBUY7S1x/uemmm+wtgkZjV2wOMzcv8/UG/IHmmPZDJQH7gT1KKZvz4Jg36saZ53gVk4J5A1iolHq+RL8jwCal1H0lzv1/9t47Po7jvvt/z+5ePxw6CIAk2ItIiiqkKImyerOqm5prZDuSSxI/KY4dP3b8o5U8yc92XGLHJUoeR5ZtWbZlW82FVu8iJRaRIkVSbABJ9HLA9bvdneePWRwOIECAIEgcwH2/Xvva3dmZ2Zkr+9lp3+96VOvts4DtpG+XUg7XsupPcx/jnGYO8MY9/01rKMI1f3srQ2cEpjft4Mg7PwlA0ufje5//Aj/4yyuOine6YLZ3owV8aCUhzOZ2sm830TVvDju707T3pvAaGgtqS1lcV3patTTHgjuLz2WqUEzTzJ8Bfi6l/MFEFUJK2eNMxvgP1JTyKKqVtm5IVAO1zqqQ2524P0K1CB9HrXM6afjLS5DdCdI5i4B38EfoO3cZ+uxarEOtBDMZSjZtY3vT2aycM6L1pWmNUVORPzYPt5HeuJ3w6zt4x5mLSJyxkF0xZT5JCFg1vxopJbbEdfXh4uIyLoF6B2os6n+AOUBcSnnCfgGc7sArRokzd5iwKPBRZxvrve48zuINIlhVhtHUSW8ie5RACSEoec+VRL/zMwAu2vEGj77WeNoKVCH+NWfiWdhAeuN2Mtv3YGzbzaqVSzjvsjX0y1Fzd5JX9rSxsC7CghmllASOaw6Li4vLNGK866AuQ00w2Ak0CSGOCCG+KYSom7CSFTHhukqMRJJob2LY65EPDsy+Wt64n7de201XzJ3BBqBXlBJ65zsou/tW/KuWo5WGCXgNfIZGZuc+PEJSEfbxZmM3D288wPqth3i7pfe09O7rdu+5nO6MV6DKgCbgv4GfAL3AXwM7nUWx05pgbSW6Juhr6Rr2umf+LAKXquEtTUoueHMrj73eeCqLWPRoJSGCl68hcL5yVWIebCbx+HN4HniUC3paeNeSKs6eV0U6Z7G9sYv+Hr++ZBb7NHlw33vvvZNdBBeXSWW8AvUssFxK+Qkp5Z1SymWoKeaHgIeFEAsnqoDFiF5Zhs+jk2wdXqAAIh8Z8OWzdsc2Hnt1P32p7Kko3pTEmDeTktveiVFfQ3rDdnL3P8zcVzZy4xk1vPOcBoQQ2FKy/o1D/PrV/Wx8u522aHJai9UnPvGJyS6Ci8ukMh6BsoHfSCnNwkAp5UaUIdYox2FSaCqiV0Tw+Tyk2rpH7IYJXXcxWrWaIFCWiLNg9y4eemX/qSzmlEIIgWduPSXvvYrST9xK4KJzkLaNFgoQ9Blk9zZhtXezZmENNZEAe1t7+dMbh/nNqwc42B6b7OK7uLicBMYjUN2ohblHIaXsQc2mu/ZEClXsCF0nUFWK7O4lkTGHj+MxiHzg+vz5VZs38MiGA0QTmVNVzCmLHgkTuOgcIu+/Pj/VOvnEy8R+/Ahlv3+a85PdvO/MWi4+o46qiB+/R03s7IqleXVPG0e6Eli26/XXxWWqMx6B2gjc5fhrGo4uYNpPWatoqGF5WM8/HIej9KPvBo+a5beg5QhzDuznFy/vO1VFnDYIIYj82bsJXnE+QgiSz2wk8V8PMePAAS5bXk+tY6C2N5nlQHuMp988wi9e2sezbzazt6V3yrqof/TRRye7CC4uk8p4FupeCjyFmsF3l5RyQ8E1P2qdVK2Uct5EFnSyGG6hLkDuYDPSMvEuaDhm+o7P/ht9P34EgLfrZ/Od2z7ED+6+mDnVJcdM5zIyVncv2bf245lbjzFzBmZrJ6lXtuJdPBdt3kza0zaHu+Ic6UqQMW1uWzsfXdNo7UniMTQqwr4psXC6ubmZ+vr6yS6Gi8uonKyFuuNyWCiE+DPgB4AP2IvyoptFjUHVA/8ipfzHCSznpDGSQIHqUoomMiyoHc5erSLX1ELT+e8Hx6TPt977AQIXnc3XPnzBlHhITgWye5tI/ull7HgSNIGnoQ7P4rl4ly0gaUPYr9ZSPb6pkZ54Bp9Hp64sSF252kL+4lxr5VqScJkqnCyBGtcsPinlj4HlwPdRInUL8AFU1953ONoCxLRDSknjG/vYvGHPMWeSeRrqKLnjuvz5Ta88z7aDXTy1/YTXNrs4eBc2UPqp24l86Eb8q1dgReOkntmI0ARhvwfzSBtWTx9XnjmTi5bWMrMiRGuvWhD8+r6OfD7N3QkyOdc2oItLsTBul++DMlEuLEqAVinltPLSN1ILSkpJx3cewDu3nrJ3XX7MPHIHm2m68AP5VtSPrr2Jt885hx/cfTGVJf6TVPLTFyklMp5EK1EuUXr/57dYHT3o5RE882fjmT8LfVYNvRk1NlUe9hFP5/jthgMIoCzsY0ZpQG1lQXzHGGc8mbgtKJepQlG1oIYipYxKKQ9NN3E6FkIIKm6/lshVF4wa1zO3ntK7b82fv/fFZ0hH43z9kTem9TqeyUIIkRcngPB7riJ45floZSVktu4i9qv1JH//IuVhH+VhH3YsQcBrcPVZs1g5txKfofN2Sy/P7WzhiOM5OZHOcbA9dkq9JN91112n7F4uLsXIhLSgpjPHGoMCONDWR2cszXkLa46Zjx1L0HTBB7DauwFYv+oCHrnoMv78qqXceuGCCS61y0jInEmuqQXN78WYOQOrL07vD3+JXlGK0VCHZ04dxuw6pN9HVyxNadCLz6OzuznKxreVh+SQz6AqEqAq4mdhbQSvMTktLBeXYqGoW1CnKzKTJf7yVg68sX9Uv0ZaSYjKL38qf37Vlo3Mbm/lf57ezfbGkS1SuEwswmPgXTAbY+YMdW4YBC9fg1ZWQnbHXuKPPEP0ez/H2n+ImtIAXttCZnMsrC3lunMaWL2gmupIgM5Ymk0F41f7Wnt5bW87B9tjxNO5CemaW7Vq1Qnn4eIylXFbUKNwrBaUNE0Of/3HbPWXseIDVzNvxkhLw5z4tk3zTX9JeuN2AFrKK/n/338ngZIg//6xi6ivONqNvMupQ1oWVmsnucYWfCsXo4WDpDfvJPn0BozaaoyGWoyZMzBm1qD5faSzJn7Hmv2m/R3sPhLNG7X1e3RqSgNculxNE8/kLLyGdlwzN90xKJepQlFNMz+dGK2Lr+/BP7B5x2G44wYuWz76mpXc/sMcuvxjyGQKgKfOPo9fX3IlDVVhvnnnWte9RJFhtnWR3X0As7EFs60TbAmaoPwzH0J4PZgdPQifBz0SxrIl0USGjr403XFlvX7tkloAfr+5iVgqR0XYR0WJj8qwn6qIPz8FfjhcgXKZKhSTw0KXAjxz6qncspcdzd3kltbi0Y/da+qZP4vKe/6Czs/+GwBXbn2NPbMa2M4i/vHnG/mXD55P0Od+LcWCMaMSY4YyjCKzOcyWDqzuXoRXCUvqmQ3kDjajRUIY9TWE6msonTUDwxGmfhbXldIZS9MVS7PrcBRbSuZUl3DJMuWh5vV9HYT9BmUhH+Uhn1qrVXdaeK9xcRkRtwU1CqO1oMzWTprv/TWvz57P2defz/xRuvn682r9wOdIPvkqACmvj6/e/hHayytZOaeCf3r/mmOaUHIpHsz2bszDrZiH2zCPtGHHkhgNdUSctW+pl7eilYQw6qrQKssQQuRbWkJARdhP1rR4eOPBQWuwAl6DlXMqWFxfhmVLepMZSoNedM0dNnYpPtwuvkliNIGSUhL94S/ZErOxrr2Yq8+aNaZ8ra4oh6++C/NQKwCt5ZV87bYPk/b5ObOhgq/cvrpoLRy4jIwdSyAzWfSqcqRlEf3eg8i0MhAsvB702ip8Zy3Bd8b8/G+qvysvlbXoSWSIJjJEE1keffA+vrruH+iJZ3h8UyMCCAc8RAJeSoNeFtRGKAv5kFK6VklcJhVXoCaJ0QQKIPHEKzS98AabLr2Ed61dNOZxpMy2PRy58dPIlHqA7ZnZwH+86zZMw2DBjAj/5wNrKA/7JqAWLpOFlBK7K4rZ0qm6B1s78S5fgH/Vcqy+OH0/+i16dQV6TQXGjEr0mgr06nKErueFK5OzaOlJ0pvM5rdYKsuly+qZWRniSFeCl/e0Uhr0Egl4KQupfXVpYNQuZxeXicAVqEliLAKVO9hM1wO/46U5i5h7wXJWLagec/6xXz9B+yfvyZ9vWbCY/77u3UhNo6Y0wLrbVrOgdvRuQ5eph9UXJ71hO1Z7N1Z7FzKnXLeEbrwU37IFlGkeWl/bil5Vjl5VjhYK5NP2L/DWhKArlmZ3c5Q+R7yyprKQcdPqOZSFfBxsj7G/rY9I0EtJwEOJ30NJwEvYb7gtL5cJwRWoSWIsAiVtm94f/pJmw0/pLVcz8zini/d852d0/9MP8+cbli7n/qtuQGoaPo/OZ28+Kz+Y7jI9kVJi9/RhtXdjzJqBFg6y3Ajzwr98Jx9HBP3olWWE3vkO9PIIdioDAjS/b1A+6ZxFbzJLdcSPrmnsa+1l5+EosVQ2Pw0e4PaLFuA1dPa19tEdT1MS8FLi9xDyG4T9Hgy39eUyRlyBmiTGIlAAqRc3Y/fFCV538XG/lUop6frH79L7n7/Kh21Zuoz/e+UN2LqaLHHdObP55DXL8utuXKY/m15/nXOWLsPqjGJ19qito4fw+65GC/pJvbSF1Etb0EqC+VaWXlmGd/kChH70JJv+ca5YKks8beZb5pv2d7CnebDfLI+ucftFCxBCsOtIlGTGJOw3CPk9hP0ewn7DnbDhkscVqElirALVTzpn8dbhHpbPLj8uEzjStun4u68T++nj+bA9ixfz/StvJOvxAjCrMsQX3nMOC+tGdu/hMn3YtGnTMa1JmK2d5Bqb8wJmd0ZBCMr+5sN5x45mcwd6ZSlaeQS9sgy9ohS94ujfT3/LK5bKkUjnMG3JIud39vzOFg51xgfZjSwLeblp9VwAtjV2YdmSkM8g6Gwhn2fSjOy6nHpcgZokjlegOva3sL6pj0tXzGR2Vfi47iVtm84v/Dt9P/pNPqx7zmy+dtW76AupvHRNcMuF8/nAxYvcqejTnONdqCttGzuWQC9VzjBTG7aR23cIq7sXmVQLh/XyCKV33aKuv7QFmTPRyyNKwMojiHBw2B4AW0pSWZNE2iSWyqEJ8pZT1m89REdfelBZZ1aEuOLMmQC88FYLuhBKuPwegl4jPx7mMj1wBWqSOB6Byh04QuxX6xHXX0r5ivEZgJVS0v3P/0n0Oz/Lh5nVFXzvqpvZXT2w+LO2LMBfXrdiVCO1LlOXibQkYacz2N29yKyJZ66yeBL71XpyTS1Q2LU3byYlt14LqDVcwu9FKy1RIlYaHrbrEAYELJlRm9fQqSsPIqXkT28cJp7OkcqY9NdmUV0pFyyegZSSxzc1EfDqBL0GAZ9BwGtQHfFTWeJXY3NSut2JRY4rUJPE8QiUtCwyW3fhW7EI4fPSFUuP2714730P0/n5b4HtPDw8Bs9ffz0PNpwBBfmtWVjNx65YOqodQJepx6kwdSRtG7svoSZoRPvQAn68S+ep9X3ffSC/hsspEP5Vywhecb7qEty4Hb20BK00jBYJI4L+Y/7WbSlJZUwSGROvoVEW8pGzbF7e1UoiY5LKmqSyFlJKzppbyco5lSQzJr9+dT9eQyfo0wl4lYAtmBGhtjxIzrLpiWeccN2d2DFJuAI1SRxvF18/7b0p1m89xPmLalhcXzaueyeefJX2T92DHY3lw2KXnM+/nXUxHWJgsoQArjprFh+6ZBG1ZcFx3cul+Fi3bh3r1q2btPtLKZGJFFY0hh2NYUf70GdU4l00BzueJPr9BwfFF4ZB4JJV+Fcvx05nyGzZhRYJqRZYaRgRCiBGaQn1r/sSQuDz6KRzFm83R0lmTVIZi1TOJJUxOWtuFQtqI3T2pfnDlqZ8eq+h4fforF5Yw8yKEH2pLAfaYgS8On6Pjt9r4PfohNxJHhOKK1CTxHgEKrf/MOk3dvHq3MV0xDNcd07DuBfc5ppaaPv4l8ls3TVQpqpyXrnjVn6iVVBYKl0TXLFiJrddtICG4xz/cnE5XmQmq8QrlsDujWH3xvEsmI1nTj1mayd99z86OIEmCF1/Cb5lC7B6Y2Tf3IsWDqKVhBAlIbSSIMLnPa4eh6xp0dGXJp01SWYtUhmTdM5i2axyqiJ+DnXGeXZH81Hprl45i9ryIIe74mxr7Fbi1S9gXp35NSX4vQaZnIVp2/g9uitox8AVqEliPAKV2bmPxOPPoV9+Pk+YAQxdcP25DeN2bCczWTq/9B367ntk8IUbLuPBtZfxfFt6cJmBtUtmcPOauZw1p9JdjDlFqa+vp7n56IfrVEFmc9h9cey+BFZfHLs3jveM+Rg1FWT3HSL+mydhyP+q5NZr8cybSe5wG9ltux3hUuKlhUPolaUI4/iWWli2apWlcybprEUqazGzIojfa9DSk+Stwz2kskrY0lkLW0redd5cIkEvOw/35P1+Gbpal+j36Fy+op6Ak76jL4XPo+MzdLX3qO5L7TT637kCNUmMR6CklMR+/nusrijZ227gyd0d1JUHuWxF/Qn9aBN/fJGOv/t63isvgBYJk/r4bfy4bhFbD/celaahKsxNq+dw5cqZhHzurKmpxHR3tyEtCzuexI4llQ3DWALvsgVo4SDZXQdIPrMBO54aJGKlH3sPelU5mR17ybyx2xGvEFoogAgH8S5sQHgMpG2DEONak5izbAxdQxMi7z4lnbPIFGwXL6vDo2tsOdDJm03dR+XzgYsXoWuCTfs7aOyI4zM0R7zUONpqx9pMR1+KbM7G03/d0PAYOro2tcRt2guUEGIZ8F3gQiAK/DfwFSnlMV3VCiFKgW8D70Z5CH4c+IyUssu5rgOfBW4EljnJNgFflFK+NoZyjWsMymzrou/+R/GdtYTmM5fz6p425taUcNHS2hMSKaunj84vfYf4L9cPCvcsnkPiLz7CA6KMjXs7jkoX8Oq8Y2kdV66cyco5lVPuD3A6Mt0FaixI20Ym06obMZbAM28WwmOQ2bmPzNZdStziifxMxPL/9SGEz0vymY1ktuxChAKqG9ERsOAVaxCahtUVRZqWCg8FTqiXwbIlWdNyWmkW2ZxFQ7Wa6r+/rY+WnuTANcfz9rvXzAPg2R3NHOqMD8ov6DN43wXzAdj4djvRZAavruP1aPgMnZKAJz+u3dmnek+8hobXo+M1tElpuU1rgRJClAM7gJ3AV4EFwDeAb0kpvzRK2vXAYpQI2U76Ninlxc71MHAI+B/gSUACfwlcBayVUm4aJf9xCRRA8qkNpDftIPy+q9njCbHlQCcLaiNcuHjGCXe7JZ/aQOcX/53cvkODwv3nrSD76Q/xO72UJ7cdJpU9Wt+rSvxcvqKey5bXs6A24nYBFimrVq1i06Zj/jxdcCZzpLPIRDLv0iS77xBmYzN2IoWdSCHjSaRpUvbJ2wGIP/Ys2bf2qww0gRYMoFeWUXL7OwHI7NiLTKYRQT9aMKD2jthNJIl0jmTWJGvaZHMWGdNGE+QFaPP+TtXKMu28CJaFfFx/bgMAj29qpCeeGZRnfXmQK1cqrwovvtVCzrLx6BpeQ8djaFSEfcxxBLS1J4mmCbyGhsdQcQzt+Fue012gvgB8Dpgjpexzwj4HrANq+8OGSXch8DJwqZTyeSdsDbABuFpK+aTTgopIKXsK0nmBPcAzUsqPjlK2cQuUNE36fvIYdiJF6Z3vZntnim2NXfk1ICeKzObo/a+H6P63+5Dx5KBr/necS+ATt/FiRT2PbWqiachbWj8zSgNcuGQGa5fUsqKh3B0IdjktMDt6sLujSsDiSWQiBbpG6JqLAOUp22xqGZRGry6n9KPvASD+6DPY8eQg8dKryvAuUS0jqzeG8HoR/uOb9DEWbCnzraSuWJpU1iSTs8lZFpmcTdBn5K2APLejmVg6R860yZo2OctmTlWYix3bng++uJdcwTo4gAW1EdYuqUVKyR+2HELvFzBH5OrKg8yuCmNLSVNHHI+hMatSTcqargL1PNAspbyjIKwBaARullI+NkK6e4C7pZS1Q8L3A7+VUv7dMe75O8CQUl47StnGLVAAVmcPffc/hl5fTclt17K7pZeqkgBVEf+48hsOs62L6Ld/Qu/9j0I2N+iad9l8Sj91By0Xruapt9p5dkczvcnssPlEAh7WLKph1fxqzp1fRVnIdfUxmdx9993ce++9k12M0xIpJTKTVd2LyTQymQJNw7tQtVwST76K1dHtXE8h01m1yPmWawCI/vAX2H0J1ToL+BEBP55FDQQvVqarUi9vRRg6IuBHBH1ofh9aacmEt9CGq5ctyXfx97fOco54ZU2L0qCXWZVKgJ7b0eyE2/l4S+pLOXteFZmcxS9f3gfARy5b0p//tBSoduD7Usp1Q8ITwDop5ddHSPdLoEZKedmQ8N8BSClvGCGdDzgIPCSl/KtRynZCAgXK71Pijy/iO3spoWvW5sO3NXZRHQlQVz4xP8rcoVZ6vnEfsQf/CNbgrj29rprIh28ieMd1vJEWPPtmMxv3thNPmyPmN39GhHPnV3HuvCqWzy53DdWeYtwxqKmDtG3ImQifspuZ3XVATQBJppCpDDKZRq+vJnD+SrUI+ls/QZqD/3v9zwdpWUR/8As0v08JWMCHFvDhWTQH78IGpGWR238YEfCj+b0Iv0+11I5zduOJYktJn+PeZYaz/nK6ClQO+Hsp5beHhB8G7pdS/u8R0j0BJKSU7x4S/lNgvpRy7Qjp7gE+D6yUUu4epWwnLFAAqVe34Zk9A2Om6tozLZs/bGmitiw44eaKck0t9P7nr+j76ePqza8QTSN49YVEPnIT3kvPY/uRPl7e3coru9vojKWHzxD1xrWgNsKK2RUsm13O8tnlVIQnrhXocjSuQE1fpJSQzWGnMsiUaqVp4SDGjEpkziT5zEZkKoOdSiNTaWQ6g+/cZQTOX4nVF6f3h788Ks/gFefjX70cqzdG4nfPI/w+JWA+JWCeRXMwairU+rWuqFpz1i9uI5iwGivTfQzqlAmUEOIG4FHg74beb0i8u4G7gVVw4gJViNnSgVFXjWnZCAG6ptHemyJn2cftS+pYWD199N33ML3/9WusjqOnwurVFYRuvoyS916FZ/Vy9rb28freDjYf6OStwz2DfAcNR115kOWzy1lcV8rCulIWzIi4rawJxBUol+GQponVGVUCls4g0xlkKoNn3kyMumqs7l4S619S4ZmsasHlzLwjzFxTC7EH/zAoT+ExCN18Od4FszFbOki9+oZqwRWImHfxXLRwEDuZxk6kHPHzgsdAc8aup6tAtQPfk1J+ZUj4WLr4qqWUlw8JH7aLTwhxHvAM8GMp5V+MsWwT0oLqJ7vvEPFfP0H43VfgXTw3H/709iMc6U5QWxbk3PlVVJZMXOtEZrIk/vAiffc/QuqFzcPGMWbXEn73lYTfcyXeFQtJZS22NXax5UAnWw500tgx/CSLQjQBsyrDLHIEa5EjWkGfK1rjobm5mfr6+skuhss0QDpd/kLXsVMZzOb2wQKWyeI7awl6ZRm5g80kn96AzGSQqWy+KzLyoRsx6mvyQxZ5NEHl339M3WeaCtTzwBEp5fsLwmYDTYw+SeIuKWXdkPB9wMOFkySEEIuBF4FXgPeOtr6qIN2ECpS0bWVQ9qwlg5rVlm2zp7mX7U3dZHIWc6tLWN4w8d1o2X2HiP3scWK/+OOgBb+FGLNrCV6zltC1FxFYezbC56UvleWtwz3saOphx+Eedh+JHjX7ZyRqSgPMqQ4zp7qEudUlzKkO01AVdltbo/DYY49x0003TXYxXE5zpGkiM7l8V6DVG8Nq6URmsqoFl8kSuvQ8FXeaCtQXgL9HTTOPOWGfBe5hbNPML5ZSvuiErQZew5lm7oTVOfFagSullMnh8hvhHhMqUIXY8STJ514nePkatKASoqxpsfNQD28diWJaNrVlQZbNKqe+Yng/PeNFWhapl7YQ/+1TJB5/bpBB2kJEOEjw8jUErzifwKWr8cyuzZfz7ZZedjf3srell7dbejnUGWesn5IAZpQFmF0VZmZFiPqKEDPDQkZ4AAAgAElEQVSdraY04C4kxu3ic5k6TPcxqHLUIt03UQtt5wPfBL5duFBXCLEXeE5K+fGCsPXAIgYv1G0vWKgbQLWa5gIfBLoKbp2RUm4ZpWwnTaCy+5qIP/w0WtBP6IZL8TQMNAQzOSUAu45ESWVNwn4Pl6+oPylTv2U2R/KZjcR/+xTJJ17B7hu5O8+YO5PgJasIXLyKwMXnolcOWGpPZU32t/Xxdksve1vU/lBXfNSxrKPuoQlqy4N5waotC1BTGmRGWYAZpQFC/tPDZJMrUC5ThWktUJA3dfQfDDZ1tK6wK04IcRB4Vkp5Z0FYGfAt4D0MNnXU6VyfCxwY4baNUsq5o5TrpAkUKLfdiceexYrG8J+3gsBF5yA8A11fli1p7IjR2BHjkmV16JrGvtY+LNsetxuPYyFzJqlX3yD5x5dI/OklzIPHNlbqXb4Q//kr8a9Zgf+8FRizawe19HKWzZGuBI0dMQ52xGjsiNPYEaO5OznIhfjxEPYbSrBKA3nRqilVa8sqS/xUhH3TYsGxK1AuU4VpL1DFyskWKHBaME9vILNtD1okTPDqC/AuaBgx/rNvNpO1LK45azYAB9tjRIIeykPjc444YrmkJLfnIIknXiH1/CbSG7blXYePhF5bhf+8FUqwVi3Hu3xhvvuykKxpcbgrwZHuBM3dan+kO0lzd4LuIaZbjhcBlId9VIR9VJX4qSjxU1Xipyriz4dVlvgpCXiK2szTvffey9133z3ZxXBxGRVXoCaJUyFQ/eSaWkg+8QpWV1StOr9k9aAutEKypoXX0MmaFr98eT9SSvwendryIPXlQWpKg4T9xsQKViZLetNOUs+/rgRr81tHLQg+Ck3Du2QuvpWL8Z21BN9ZS5RohQIjJklmzLxoNfckaYsmae9N0dabor1XrXyfCHRNUBr0Uh7yURbyUhbyURryUhb0UR5W+/7wspB33O5SXFymO65ATRKnUqBATV5Iv7aD9CtbwZaUfur2YVsghSQyOVp7krT0JGmJJkk7BmL9Hp3qiOr6aqgKEwl6J7SsdjxJ+vUdpDduJ/3am6Rf33GUTcBh0TQ8ixrwLV+Id+k8vGfMx7t0PkZD7Zg8rkYTWdp6k7RFB0SrrTdFdyxNZyxNNDG8KacTJeg1iAQ9RAJeSgIeSvL7gbCh+5DfM+4JH24Xn8tUwRWoSeJUC1Q/diKFeagV71JlfDK1YRue+bMxqsuPma7/Ad7Rl6KjL01HX4pYKsfaJbUsqI0QTWTY3tjNyrmVlAa92FIiYEJaWtKyyO46kBeszNbd5PY2HeWUbiRE0I93yTxHtObhXTIPz8IGjJk1x7XS3bRsuuMZuhzBUsKlzgfCMiSzI5t5migEEPJ78kIW8nkI+QzCfg9Bv5E/D/kNwr7BYXNn15Ps68HQp/54msv0xhWoSWKyBKoQO56k978ewr/2bGXLy7KUM7YxTgRIZ03HpL5OS0+SV/a0cdWZM4kEvew6EmVbYxelQS+RgJdIULUMIs4D9UQnG9jxJJk3lXO5zLbdZN7YTe7tJrDH3k0nfF6MeTPxLpiNZ/5sPPNnKdfiC2ajV5ePW1wzOYveZJZoIkM0kSWadPaJIXsn/HhnI04UPo+uRMxnEPJ7CPoMAl6DgFd39gPHQZ+B33GKF/AZBDy62nt1gl4Dn0cv6nE3l6mJK1CTRDEIFICdyiA0gfB5yby5l9Rzr+FdtgDfmYvQq47dqjoWrT1JDrTH6Etl6UtmSecGxpT63/4jAQ8XL6vDa+hEExmypk11xD/uB52dSJHduY/srv1k3zqg9rsOYHX0jJ54CFpJCM/8WRgNdRgNtXga6jFm1+KZU4cxq3bU7tGxIqUkls4RS+WIpbLOPkefczywL7yePaYx3slAAAGvgd8RLL9Xz4uab8jmd1yY+7163qX5sPEK9l5DcwXwNMQVqEmiWASqkNzhNtKvbVfOCm2JUVelusMWz0UvKzmhvDM5K//A7UtliSVzxNI53nnObDQheHVPG02dcW5buwCA1/a20x3PEOp/q+9/a3cefgGvMeaHltnRQ273AbJvKcHK7j5Ibv/hYe0IjhW9usIRrjqM2bVKyOprMOqrMepr0MpKTuoD1bIl8fSAaCUyJon0kH0mRyJdeN4fliOZMZmkhtu4KRStQlfnXkMJmNcYcJ7Xfzx471zXC657nL0+NL2Oz6N8FbnCOHm4AjVJFKNA9dPfEsns3IfVptYf6zUVeBfPxbNgNsaMygm/Z9x5aNaUqll42xq7aOlJksiYpLPmsN1gkYCHdzkurrc1diGE4MyGCgCaOuMIBh5qXuehVviwsfri5PYfJrfv0KAtu+/Q2CZlHAMR8GHUVaPXVSvRqqseJGB6bRV6VdkJW3seDzfddBOPPvooqaw1IGIZ9fmnsxbJrPrMk1mLVPbosHTWVHFzlpPGJDNBMyCLkUGC5jjYM3TlKdbQBR69IKzgmjoW+TCvsy+83h+nP60nn+/g60fnq43LQ+1UwxWoSaKYBaoQKxojt+cg2T2NmM3tGA11RO64DoDs3iaM2qpT4gwtZ9mkshapjEkqa5LKWWgCls5U3ZDP72xB1wQXLVUmk3796n6SmcHdYAKUUDliNaMsyDnzqgB4s6mb0qCX2VVhpJQ07W1GHG5Da21HHGmHQy1wqBmzqQXzcPvo0+DHgqahV5Wh11Si11Rg1FSgV1eg11Sgz6jAcML1mgq0SHjCHkYnYxafZUvSWZPUIDEzyeSUN9Z0Th2nnXMVbuXD0gXn+XDHFXk6a43ZPuPphqEJdEesDF1DL9wPF+bE7U+ja45Aahq6PiRNPkzF0YfsDW1oniOkce6na2IgX+dYH3JNG/IbdwVqkpgqAlWIHU9ipzIY1eXY6QzR7z6A/4KVBC9ehczmyO5pxJh58ru3xkI8nSOdtciYBQ+9/LF6QFaU+Fg1vxqAX72yj4aqMOcvmoFlSx544e2j8hSg3m4FaPEks3NJFsd6yDW28Ep7krrWFqoPNJJq7eJguAzDttAta2Bv2ei2hWFZ6LaNxzLRxzqpw+tRAlZTgV5Vjl5RilZZil5Rhl5Zhl5ZOhBWWY4WCY042WUqTjO3bEnWtApEbmDr98qaNa28+/H+4/7wnLPPON5bB65bebflWdMiaw1cV+7OXWE8lQjA0DU0R7Ae/vw7gYkXKNec9DRECwfzrSXh8xL5yM2IgLLhZzZ3kPj98ypeKKC6seqqMWZUos+onLBJBWMl7PcQPg7berdcMD9vkFYIuO6chqMeeP0PMtOyyZUGKa0IEakrxbRsrE2NlMwqp6G+jL5klg0v7sGOJbBjSbXFEwPncbVfsmc3C1sOk/T5eO7Mczlr/9vUd3cSDYXZNn8Rum2h2dLZ2+i2jd6TQ+9qRbObqe/qoDSZIO3x0l5WTk20B38uS9rjpa+kBE8kiFFaglEaxltWglFRilFWwkdK5tHxwO/xl5dglEXQykrQy0rQSksQgYm1GjJR6JrIzyw8ldhSYloFAuaIlmmpff9mWpKcOTh8UBxTxTk6nT2Qzh6cR/99TXsg/8J0kzX782QiUWbMmIAOimPhtqBGYSq2oI6FlBKrowezuR3zcBvmkXbs3gFL5lokjDGjksAlq9Ary9SUdm16DkBLKbHs/oeR2ucfMs7DxbRsKvwGZekk8eZOduxvZ1YiRqS7m66OXt5M2ORiCcxYkmw8hZUzsTQdW9OwnO3sfXuo7+6kM1LKq2ecyYU7t1EZ6+NwZTVbFy4ZtZwX7txOZayXI5XVbJu3kEu2byEkLQ7Pncf+hjnofi+G34se9GEE/BhBtelBP0Y4wLmzyghVldImDZpzsHp5PYbHw5GuBJ2xNJom0ARoQo2V6JpAgBMuaKgOowlBbzJLKmtS67j3jqVy5CwLTah4miaUA04nn/58p4NdxBPBsiWW87sybRtr6N6WmJYTx5ZYltr3i5tpDQ4fNmxoGud+Vn5fkGbItaFlGdgKzq2Bc3MYwf3Tl28E3C6+U850E6jhsFMZrLYuzLZOrLYurLYuwrdei15WQnrTDlLPb6L0k7ejBXyYze3YqbTqniqduPGW6YKdSKkXgLZO7O5erK5erK4oVncvua4oqWgcvaMbuqOkown6TBspNGxNYDt7S2hITcMWAlvTqO3uJJDN0hsMcaSqmoXNh/GaJm1l5TRX1gxKO1weF7z1JoFshgO19bw9czZXbt6IEfCzc9FiDtbPQvg8CK93yN6jvKl6Pdw6K4CntIQtaUFjFm5bMwdREuLFxihN3ceepOI1NG6/aCEAL+xsoTeV5cZVcwB4atthZXexQBw1QX7huBCCSMDDpcuV08YNb7ehaxqrF6ju3pd3t5LJWU46JY6CgbyEgLKQNz/++WZTN2G/h7k1JflzUC3x/L2FGLRwPRLwMMMR5H6blxVhP5YtOdKdcOIOTicgX6egsyjblpIeZ7ar32tg2TbxtJmPK1D37z/uz9PjjA/1v0z1vzT0P48m6//XL2C2I2DhgLJS43bxuUw4WsCHNrcez9yjvbfqNZX4zjkD4Vc/wPTmt8ju3Aco75xaaRitPIJWWoJeXuJ0Q0VGtCE43dFCAbRQYNjPcjhkNofVPSBidldv/vzbX/ln/uIjd2JHY9jRGJ7eGBXRPmxdQ5owI9rDjOjY147Na21mXquyTi+TKc544w2WvvEGthBIocRNCoGtCSTOXgg6UikEEPH5WOTx0hhXLe5wKMz80giEAohAABn0Q8CHDPjB2XS/j66nnkQLBYj4ggR8fmJNb6OFAlSYBj7NAJ8X6fOAzwc+L2gatpRIyaCuQoHz8Hfon7UopVRT8Z19/7lEkjWDLJ2p4u9r7aWmNJgXqK0Hu0Z98VwwI5IXqJd2tbJstnIimrNsnttxbEv/ACsaKjhnXhVZ0+b3m5s4b2E1S2eW05fM8fimxlHTn7ewhqUzy4gmsjy+qZFLltUxp7qE5p4kT28/4gicGEYo1Se2dskMZleFaY0meWlXK5ctr6eyxE9jR4wtBzqPSg8DQi+E4MLFMygP+zjSleDNQ91cfEYdQZ9BY0eMfa19+bgnC7cFNQqnQwvqeLDTGeyuKFan80CNxrCifdjRGDKnZuOJgI/yv/ogAKmXtiBNi+ClqwEwm9vBMNTkAJ/XbYEdg2NNkrBTGezemPP5x4Y/dvZ2NIbVF0fmx9oSYzY/NSl4PXmhF6EAWjCACPrRAj6E3+cc+9U4XMAJD/jRgn7l9TXgHwgfNq4PYRjO2NBgUUNKJGpMCwm6ruH3qCUGvcksPkPD7zXyLSJUEiRKUKWTHpWckM8gEvRi2TbN3UnKQj5KAh6yppV3OSOdBM4t83lICTNKA5SHfaSzJvva+phVGaY0qDxcH2iLFdxXlQFJ3o2NlLCwLkJF2E9PPMNbR3pY0VBBJOClNZrk7ZbefDkLy+2UBinh3PlVlIV8HOlOsPNQD+84o5aA12Bfax97W3vz6a4/d45zT7eL75TiCtTYkFIiEyn1cExn8u5CEutfQpoW4RsuASB670PYUeUgWXgMtEgYLRJCCwfVwygUQK8qxzNHtUBkzhzkH+t04mTN4pO2jUym1WSQvnhetOw+Z4JIvOB4UHhyUPwTXYM2qXiMwcIVVHvh9w0IYcCnXqJ8XoTfi+b35Y+Fr/DYq9L1Hw8JF4Xh0/SFzJ1mPkm4AjWxmC0d2L3xgQddX1w9/BJJ7EQKbIlnUQMl77kKgOj3fo5nUQOhay5CSkni0WfVW3HIeavuf5D0v0EHfJOyqPZk8Nhjj3HTTTdNdjFGRFoWdiKlxC6RRCbS2ImUCkuknLDUkLAUMnl0WOHxhKxdK1IGCZjfd5SgaUMFza/GAQc2L3iNgbFCj5EfK1TnBeOHQ67h8QyMM3qNCf2fnCyBOj1fTV0mDaOuGuqqh70mpUSmMoMMyfrPW4HWP56VzWF1dKuHWWZklxr+NWcSvOw8ZM4k/sjT+M5eindhA3YqQ3b3AfUgKHzL7T82iuvvsGrVqskuwjERuo4eCUMkPGF5SimRmawSq2Ra7eNJZDqjzlMZ7FQamVLHA+eDw/vjylQ6fz2fRzI9aV2cMpNVv92+Sbn9YHR9QOQKJsfgHSJyIwie8AwI38miuP6RLqc1QgjEkHVY/jVnDlz3eSn98/cB6u1dPZQyzkMprY6TaSWCqO5BmUznx8bs3hjJP7088v11HeH3ErjsPHzLF2JFY6Re2KRc2ddWYfXGMA82O39g50/d/6d19hgTZy185syZp13LXQiB8PvA70OfeEtdgNMbks3lfztHiVk6g0xm8mJip7PITAaZzqprTrhMZwfipIYPl2mVzs5kIZs7ORUaL5aFTFoU8y/MFSiXKYnQdUTBguTh0IJ+Ih+5OX+u11RQ9uk78g8Ze9CDZOBYc1oEMp3BbOnMC5zV2kli/UvHLpgmKHnfNXjmzSTX2Ezq+dcJ3XApekUpucZmsrsPFgjaQJcMHh1hGAiPgV5dgfAYGAj1tu0tbtf0Uw0hBPi86D4vnKBx5eNB2jYykxskXHnhyygRk8MJXTar0uVM5zw3sGVykOs/7r9mFqQZiCdzOcgMxJ0KuALlctogNA0RDsIYbRIatVWU3X1L/tyzsIGyT90++EGQG/jz9++1/oeepqnWgOMq3urpI7vrgHpQHMM0T+nH34teWcaZRpief/8pZZ/5IMLvI/XSFjLbdquuSK/HETQdPIY6dsICl6xCeAxyh9uwu6L4zlKLgc3WTtWFaugIQ1ddPIauzvuPPcaY/Yy5HB9C05RFl4APOHXCOBxSSjCtowVvkPipawwVvMzgeORM+OyLJ6WcrkC5uIwRoeuIkhCUhMYU3zO7Fs/s2vy5/+yl+M9eCoA0zQGBy5pgmuoNOWeiOflfdsf7CF52nmphAXplGZ65M1U80xp4M06kkDkrn0fgUjV2ld11gOzOfXmBSr+6jeyeg8euo9dD+V9/GIDEH1/E6ukj8v7r1fn6l7C6oo6oOYPsho7w6AjdAENHKwniP3eZuv/uA6DreBeqGZ25/YeRtq3S6RpC10DXlZAbTpjHkze3JW1brdFxW48TjhBCvYxM1AzZz945MfkMwZ3FNwruLD6XqUr/G25/N6jV04dMKjGTlqneoE1L7S21Rwj8560AIL11FzKRInDROQAkn9qA2dENOTMfX5qWGstw8tEqSym9890A9P3kUYTfR8mt1wIQ/eEvsPsSxyyzZ249Jbe9Mx/f01BH6HpnicL3H1SRdN0RN80RuwHB88yfhX/VcgASf3gBz4LZeBfPRWZzpF7c7Aii441a0/L5oGkITUOfUYlRW4U0TXIHjiijv6UlyGwOs71bxRdClUETR+fnVa3Zybb0cKpxZ/G5uJxmrFq1ik2bNo07ff/U5H708giUR8acvr+110/wyvNHTVP4Ihe+9VoKR+BLbrlmsKBZtmolFZyLUGDg/qtX5LtLpZR4FjaAbSMtW4mjrdJg2apFmcsh0wNjK7mm1rxFE5nNkXljt7qfLUecxRdYe7YSqHSW+G+fInj1hejnnIHV3Uvsgd+NWv/g1RfiP+cMrLYu+u5/lPB7r8K7sIHcwWbijz+bFzIliMI5FqqlqAkCl63BM7sWs7md1IubCV55AXplGbmDzWTe2K3iFqYVg899556BHgljtnaS238Y/6plCJ8Xs6UDs6VjwK5m/32HlMPTUI/wGGqRdyyJMbMGIYRaCpLODlNmrd8ExUnBFSgXlyJl8+bNk12E46awxaD5fYOu6VXlx5WXf/XyQfmGrll7XOnLPnHrQFnCQcr/5iP5cyVUjrhZ6ljaNsKZMi0CPiIfuTnf3aqXRyi59VolirYsENeBtNg2xswZ+fsF1p6tXgoAEfTjXTwHbKnuJ6W6vy2R0hFN21YtNEDaUk2A6C9vOoPV0Y104iFtZQXCHkiLLfEunQeOQKVe3Ixv5WKEz0tu/2FSL20Z/TP79B0Ij0H2zb2kXtpC+d9/FIDUy1uVQJ5i3C6+UXC7+Fwmi6noD8qlOJD9Aui0mGTOVBMapBwiarYSPcdekl5djtB1ZTIrGsvblDRbO7F7YwMC6YisdNIFnHFH15LEKcYVKJfJor6+nubm0Q2SurhMNq6po0nCFSgXFxeXY3OyBMpd8ODiUqSsW7dusovg4jKpuC2oUXBbUC6ThTsG5TJVcFtQLi4uLi6nFa5Aubi4uLgUJUUjUEKIZUKIp4QQSSFEsxDiHiHEqA5LhBClQoj/EUL0CCF6hRA/E0IcZQdZCPEuIcR2IURaCLFTCHH7yamJi8vE8Prrr092EVxcJpWiECghRDnwJGrd+buAe4C/A74yhuS/BC4D/hy4EzgPeHhI/u8Afg08A1wH/A74uRDimgmpgIuLi4vLhFMUkySEEF8APgfMkVL2OWGfA9YBtf1hw6S7EHgZuFRK+bwTtgbYAFwtpXzSCVsPeKSUVxSk/T0QkVK+Y5SyuZMkXCYFd5KEy1Rhuk+SuA5YP0SIHgQCwKWjpGvrFycAKeVG4IBzDSGED7gc1dIq5EHgQiFE6YkX38XFxcVloikWW3xLgacLA6SUTUKIpHPtsWOk2zVM+FvONYAFgGeYeG+hBHox8NpoBTxdrBK7FBfu787ldKZYWlDlQHSY8B7n2omk698Pjdcz5PoghBB3CyHcUWoXFxeXSaJYWlBFh5TyXuBeACHE61LK1ZNcpJOCW7epiVu3qcl0r9tE51ksLageYLixoHIGWjrjTde/HxqvfMh1FxcXF5ciolgEahcDY0YACCFmA0GGH2MaMZ1D4djUPiA3TLylgA3sGUd5XVxcXFxOMsUiUH8ArhVClBSE3Q6kgOdGSVfrrHMCQAixGpjvXENKmUGtf7p1SNrbgVeklL1jKN+9Y4gzVXHrNjVx6zY1cet2HBTLOqhyYCfwJvBVlMB8E/i2lPJLBfH2As9JKT9eELYeWAR8FtUi+irQLqW8uCDOO4Bngf9ALeK93on/Tinln05q5VxcXFxcxkVRtKCklD3AlYCOmlL+FeBbwP83JKrhxCnkdlQr60fA/cAm4D1D8n8RuAW4ClgP3Ax8wBUnFxcXl+KlKFpQLi4uLi4uQymKFlQxMl7jtcWEEOJOIYQcZvtkQRwhhPjfQohDQoiUEOJ5IcTZk1nu4RBCLBRC/KcQYpsQwhJCPDtMnDHVpdi+2zHW7eAw32PrMPGKrW63CiEeFUIcEULEhRCbhBDvHybeXUKItx1jzpuEEFcOE2emEOK3QoiYEKJTCPEfQojgqanJ0YylbkKIZ0f4D/qHxCu2ut0ihHhZCNHlfCe7hRBfEkJ4C+Kc9P+buw5qGMSA8dqdKOO1C4BvoAT9S8dIWqxcgZpw0s/+guN/AP4R+HvUzMe/BZ4UQqyQUh71AJxElqPGDl9FWQYZjlHrUqTf7VjqBvAA8N2C82zhxSKt29+iTI/9DdCJqucDQogqKeV3AZyH+g9RtjdfBD4KPC6EOE9K+aYTx4Pqns8CdwBlqHHqMuBDp7JCBYxaN4dngP89JG2m/6BI61aJsu7zdZSRgzU4tlGBv3TinPz/m5TS3YZswBdQ66MiBWGfA5KFYcW+oay7SyA8wnU/0At8uSAsBHQA/zzZ5R9SVq3g+CHg2fHUpRi/29Hq5oQfBP5tlHyKsW5Vw4Q9ABwoON8N/Kjw8wC2Az8tCHs/YAHzCsJuQ02MWlTEdXsWeGiUfIqubiOU8/+gxEqcqv+b28U3POM1XjvVWAtEKDCkK6VMoCaqXDdZhRoOKaU9SpSx1qXovtsx1G2sFGPdOocJ3gLUAwgh5qPsYRZ+bzbwK47+3l6TUh4oCHsY1ep45wQXe0yMVrfjoOjqNgJdQH8X3yn5v7kCNTxHGaGVUjahVH+4hcHFzj4hhOn0I3+iIHwp6s3t7SHxC43tThXGWpep/N1+XAiRFcox50NCiDlDrk+Vul3IwAL5/nINZ8y5QghRXRBvaN2yqIX4xVq3fq5xxl+SQoj1QoiVQ64Xbd2EELoQIijUUp3PAD+Qqhl0Sv5v7hjU8IzXeG2x0YLqI96Imp5/B/BDIURQSvktVF3iUkprSLoeICiE8Dp/lKnAWOsyVb/bR1BjVIeBM1BLMF4QQpwpBxabF33dnMkP7wY+5gSNxZhzB1OzbqCWwPwY2AvMAb6I+t7OklIedOIUc90SgM85vh813gSn6P/mCtQ0Rkq5HjX42s8fnNlDXxJC/PskFctlHEgp/1fB6QtCiJeBragJBd+enFIdH0KIuagxmkeklPdNamEmmJHqJqUsXMv5ghDiSVSL4q+drdhZizI5twb4MsrYwadP1c3dLr7hGa/x2qnAQ0AFMBdVl/AwUz7LgeQUaj3B2OsyLb5bqWa37QbOLQgu2roJISpQ5scagQ8WXBqrMeepWLejkGp220tMke9NSrlZSvmilPKbqC6+TwkhFnCK/m+uQA3PeI3XTgVkwX4Xqutv4ZA4IzmCLGbGWpfp9N1KBr5PKNK6Oet5HkcNsN8opUwWXO4v13DGnLullB0F8YbWzYsyi1asdRuJsXxvk163Ydjs7Odxiv5vrkANz3iN104FbkGt2WgEXgb6KDCk6/zhbsIxtjuFGGtdpsV3K4RYgfrjbyoILrq6CSEM1Iy8RSjbl+2F16WU+1GTCgq/N805H/q9nTdkYsjNqPGRP56c0h+b0eo2Qppa4B0c/b0VVd1G4CJnf4BT9X+b7Ln1xbihmp8twBMo+313A3GKbG3QGOrxa+DzqKmeNwI/Qb25/VVBnC+gZtT8Bcoe4u9QAjZjsss/pC5BlLjeArwC7Cg4D461LsX43Y5WN+AG4Oeo7qPLgU8BR1ALriNFXrd7nd/cZ4ALhmw+J07/OqAvOfW7z3mArSjIx4MyJr0JtSD2/UArBWuliq1uwErnN3inU68/Q7UauoGGIq/bH1EGta8DrkHZR40DDxbEOen/t0mp/FTYgGWoldQp5wP+Jw8qMsMAAAhRSURBVECf7HIdZx3+BTVOkXTqsQn48JA4AjWz6LAT5wXgnMku+zB1mctA18jQbe7x1KXYvtvR6uY86J5CzWbLOQ+v+4D6KVC3g6N9b068u1Az3TKorqQrh8lrFmp9UBy1Jud7OC8nxVg3YCbwe+d7yDpl/jWwdArU7Z9QohlHzcLbDPwV4CmIc9L/b66xWBcXFxeXosQdg3JxcXFxKUpcgXJxcXFxKUpcgXJxcXFxKUpcgXJxcXFxKUpcgXJxcXFxKUpcgXJxcXFxKUpcgXJxGSNCiLmOu+51k12W8SCEqBJC3O+43ZZiGNfyJ5D3s0KIgxOVn4sLuALl4nI68Q2UmZkfAh9GeUg9qQghbhJCPCGEOCyEyAghWoQQLwshviaEqCqId5kjmlIIcdcIeUkhxOMnu8wuxYPrbsPF5fThapR303tOxc2EEF9FuffeBnwfaEN5mz0T+CTKG+twXmnXCSF+KqVMnYpyuhQvrkC5uBQpQggPyiRMeoKyrEXZgTvpCCFqULbcXgMuklLmhlwPj5D0dWA1ylfSv57UQroUPW4Xn0vRIIS40+nGuUII8VkhxD6nW2iPEOLPhsQdcTxICLHOuTa3IOw+J6zSOe4UQsSEEA87FqYRQtwthHhLCJEWQuwSQrzrGGV9vxBimxO3ybnnUS98Qog6IcQPnDhZZ/znXucBPlyZlwshvimEOAykUYZHj/WZhYQQ/1rwWbU640xzCuKsE0JIlO20PyvoSrtzlLwrhRA/EkJ0CSHiQoinhRDnjHG8aT7q+fL8UHECkFLGpZTxYdL9EmUz8vNCiMpR7tFfztVCiN8632lGCLFbCPHF4b4Pl6mF+wW6FCP/AgSA/0QZD/0UcJ8QYq+U8qUTzPuPKOOWX0b5svkM8FshxG9Qlpb/L0oYPgM8JIRYLKU8MCSPm1EP4O+hDLfejHLBPgfl4RYAIUQDyjq518l3n3PPTwGXCyFWywF37f38DGVU8xsoo6MtI1XEaWGtR7lBeMhJs8jJ/xon/8PAb1CGWH+CMuh5r5PFy8fI2wc8CZyNMky7EWW09knG1grb7+xvFEJ8U0rZPIY0oOr8Dyjr118E/vZYkYUQNzBQv284ZbsQuMcp+60jp3YpeibLWq67udvQDeWWQAJbAG9B+EyUUP28IGyuE3fdMPms42hr2fc5Yd8bEvebTngTg11XrHTC/3WYe1rAuQXhAvitc+2CgvBHgHZg1pB7rgbMwrIXlPlZwBjj53WXk+ZrQ8JvcMJ/MiRcAveNMe9PO/G/OEL4wSHhzw4T9l0nbgZ4HvgayoVI+TD3u8yJ+1nn/E+oF4U5Q8r/eMG5H/WC8PzQzwz4Gyf+ZZP9u3a38W9uF59LMfJ9WeBuXkp5BOXUbtEE5P3tIecvOPv7pZR9BffchnLINtw9n5BSbi6IK1EPX4D3AAghSlE+uB4F0s4U7ypn5tpB1Bv/NcOVT0ppjrEu7wFshozVSCl/B2wF/l975xdiVRWF8d9XVBhaQyQV9JCDZPbQH8Wgh6BSyR5iwgeLMnuQyBD0ybKwp3yZpCAoehEr8ykIBgND+kNEFENCFDU1WRYMFIoyYhiE1OphnZvjnX3n3pnpNufm94NhX/beZ+/FZu5ZZ+313bMHqsP/ZsJ9pCN+qal+N9Ac9bViC7CBjNRuA7aRB/z9KmlQk48Ln8hTZOT53BR9VgNXAa8BfU1rfKDqU1pj0yPYQZk6cqRQdwLoKCcxzbHHq7J5G6/RVprz20LdSFX2V+US8vu1kTzHqflvCXlzbeb7VoYXWAT8EhHjhbZvgAXAlYW26Yx9Tp6oenAordUkInkzIu4CLgNWkNt2p0h135NTXPsF1SGNkm5q0W1pVe5h8vo2jhMvrbHpEZyDMnXkzxb1mvB5qoPMWv5fR0SrsTuZczo0rtsHvNGiT0lG/fsM56s1lWM7BByS9Dbp5DcytVJvB7klOEie7NpMY423kRFjiU5zX6aG2EGZXqWRqL+i0NZfqPs3WVqou7EqGxHaD6QTvTgi3u+SHUeANZL6IuJkwZ5TlH9n1Ak/A6skzZ8YRVXCjEXkKaszIiJGJY2TucWp+v0k6VVgq6Q7C10OV+XpLq6xmUO8xWd6koj4jUyQ3y3pnyhHUj9wf5enXy1p2YQ5xdntqqHKvhNkHmStpElScSULZ2nHEPkd3t409r3ArcD+iPhrhmO/A1wIbG2qfwy4vN3Fkq6WdEuLtjvIB4uRUnsTO0lH+3yh7SApQtkuadKDiqR5khZ0MIepKY6gTC/zMnkDe1fSEPmWgk3A12S+o1t8CXwo6RVSBj4ArCJVc59N6PcE8AnwsaS9pDrxAjLCGwD2kuq9mfI68Cj5m6HrSDXbYlJpdxR4ZhZj7wYeB3ZKWsxZmfk6Mjpsd++4Fvhc0jDwARntXQLcDDwMnOnEvog4LmkXBbFERJyWtIF01KOS9lS29QE3AGtJIclH7eYx9cQOyvQyg+TT/COkTHmEzGssp7sOaj8wCjxNih2OkTfQc26iETEmaTmpSBsA1pPS6TEyQnlrNkZExBlJ95C5mgfIG/JJUim3IyLGZjH2H5JWArsq29cBw8BK0nld2maI74DNpNLuQVKscBHp0IeAFyohRCe8SDrdawp2HpS0gowi1wMLSXHLj9V1X3U4h6khSoWsMca0p5KGHweGI2LNXNtj/t84B2WMKSJpXqF6E7mF9t5/bI45D3EEZYwpImkf+baGT8m3QdwOPERuny2rhCrGdA07KGNMkUqAsBm4HphPCi8OAM9GxNG5tM2cH9hBGWOMqSXOQRljjKkldlDGGGNqiR2UMcaYWmIHZYwxppbYQRljjKklfwPgP6K130XWRAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "\n",
    "save_plot = False\n",
    "\n",
    "color_pantheon = 'steelblue'\n",
    "color_future = 'crimson'\n",
    "    \n",
    "plt.plot(num_list, pantheon_h0_03_list, color_pantheon, label='Pantheon SNe sample', linewidth=3)\n",
    "plt.plot(num_list, future_h0_03_list, color_future, label='Roman SNe sample', linewidth=3)\n",
    "\n",
    "plt.plot(num_list, pantheon_h0_10_list, color_pantheon, linestyle='--', alpha=0.5)\n",
    "plt.plot(num_list, future_h0_10_list, color_future, linestyle='--', alpha=0.5)\n",
    "\n",
    "plt.plot(num_list, pantheon_h0_00_list, color_pantheon, linestyle='-.', alpha=0.5)\n",
    "plt.plot(num_list, future_h0_00_list, color_future, linestyle='-.', alpha=0.5)\n",
    "\n",
    "\n",
    "plt.ylim(0, 0.1)\n",
    "plt.xlim(0, 300)\n",
    "\n",
    "line1, = plt.plot(-10, 10, 'k-', linewidth=3, label='realistc microlensing')\n",
    "line2, = plt.plot(-10, 10, 'k--', alpha=0.5, label='extreme microlensing')\n",
    "line3, = plt.plot(-10, 10, 'k-.', alpha=0.5, label='no microlensing')\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "import matplotlib.patches as mpatches\n",
    "patch_1 = mpatches.Patch(color=color_pantheon, label='Pantheon SNe sample', alpha=1)\n",
    "patch_2 = mpatches.Patch(color=color_future, label='Roman SNe sample', alpha=1)\n",
    "\n",
    "\n",
    "\n",
    "# Create a legend for the first line.\n",
    "first_legend = plt.legend(handles=[line1, line2, line3], loc=(0.55, 0.55), frameon=False)\n",
    "\n",
    "# Add the legend manually to the current Axes.\n",
    "ax = plt.gca().add_artist(first_legend)\n",
    "\n",
    "# Create another legend for the second line.\n",
    "plt.legend(handles=[patch_1, patch_2], fontsize=10, handlelength=1, loc=1, frameon=False)\n",
    "\n",
    "plt.vlines(x=144, ymin=0, ymax=1, color='k', linestyle='--', linewidth=1)\n",
    "plt.arrow(x=50, y=0.085, dx=90, dy=0, width=0.0001, head_width=0.005, head_length=5)\n",
    "plt.text(50, 0.09, '~10 years LSST')\n",
    "\n",
    "#plt.hlines(y=0.01, xmin=1, xmax=num_plot, color='k', linestyle='--', linewidth=1)\n",
    "#plt.text(2, 0.011, '1%')\n",
    "plt.xlabel(r'number of glSNe', fontsize=18)\n",
    "plt.ylabel(r'$\\sigma(H_0) / H_0 $', fontsize=18)\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "    \n",
    "if save_plot:\n",
    "    plt.savefig('/Users/sibirrer/Science/Publications/glSNe/Plots/number_forecast.pdf', dpi=None, facecolor='w', edgecolor='w',\n",
    "            orientation='portrait', papertype=None, format='pdf',\n",
    "            transparent=False, bbox_inches=None, pad_inches=0.1,\n",
    "            frameon=None, metadata=None)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}