notebooks/5_Bayesian_inference.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bayesian inference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys, os\n",
    "import logging\n",
    "import numpy as np\n",
    "from scipy.stats import uniform, norm\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "sys.path.append(\"../\")\n",
    "from simulation.prior import get_grid\n",
    "from inference.bayes import Posterior\n",
    "import paper_settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "paper_settings.setup()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 100\n",
    "n_events_for_expectation = 5000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = \"calibrated_alices_full_sgd1e2_grid\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def narrow_prior(thetas):\n",
    "    f_sub = thetas[:,0]\n",
    "    beta = thetas[:,1]\n",
    "    \n",
    "    return uniform(0.001, 0.199).pdf(f_sub).flatten() * norm(-1.9, 0.1).pdf(beta).flatten()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def flat_prior(thetas):\n",
    "    f_sub = thetas[:,0]\n",
    "    beta = thetas[:,1]\n",
    "    \n",
    "    return uniform(0.001, 0.199).pdf(f_sub).flatten() * uniform(-2.5, 1.0).pdf(beta).flatten()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "llr = np.load(\"../data/results/llr_{}.npy\".format(filename))[:,:n_events_for_expectation]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "resolution = 25\n",
    "\n",
    "f_sub_1d = np.linspace(0.001, 0.200, resolution)\n",
    "beta_1d = np.linspace(-2.5, -1.5, resolution)\n",
    "\n",
    "theta0, theta1 = np.meshgrid(f_sub_1d, beta_1d)\n",
    "theta_grid = np.vstack((theta0.flatten(), theta1.flatten())).T\n",
    "\n",
    "bin_size = f_sub_1d[1] - f_sub_1d[0]\n",
    "alpha_edges = np.linspace(f_sub_1d[0] - bin_size/2, f_sub_1d[-1] + bin_size/2, resolution + 1)\n",
    "bin_size = beta_1d[1] - beta_1d[0]\n",
    "beta_edges = np.linspace(beta_1d[0] - bin_size/2, beta_1d[-1] + bin_size/2, resolution + 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate posteriors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "bayes = Posterior(llr, theta_grid)\n",
    "posterior_narrow = bayes.posterior_based_on_expected_llr(n, narrow_prior)\n",
    "#posterior_flat = bayes.posterior_based_on_expected_llr(n, flat_prior)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate threshold for credible intervals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_ci_threshold(posterior, ci=0.95):\n",
    "    post_sorted = sorted(posterior.flatten(), reverse=True)\n",
    "    cumsum = np.cumsum(post_sorted)\n",
    "    print(cumsum[-1])\n",
    "    index = np.argmax(cumsum > ci)\n",
    "    return post_sorted[index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9999999979583987\n"
     ]
    }
   ],
   "source": [
    "# threshold_flat = find_ci_threshold(posterior_flat)\n",
    "threshold_narrow = find_ci_threshold(posterior_narrow)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot posterior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'\\npmin, pmax = 1.e-9, 0.1\\nz = np.clip(posterior_flat, pmin, pmax)\\n\\nfig, (ax, cax) = paper_settings.figure(\\n    cbar=True,\\n)\\nplt.sca(ax)\\n\\npcm = ax.pcolormesh(\\n    alpha_edges,\\n    beta_edges,\\n    z.reshape((resolution, resolution)),\\n    norm=matplotlib.colors.LogNorm(vmin=pmin, vmax=pmax),\\n    cmap=paper_settings.CMAP2\\n)\\ncs = ax.contour(\\n    0.5 * (alpha_edges[1:] + alpha_edges[:-1]),\\n    0.5 * (beta_edges[1:] + beta_edges[:-1]),\\n    z.reshape((resolution, resolution)),\\n    [threshold_flat],\\n    colors=\"black\",\\n    linewidths=1.,\\n    linestyles=[\"-\"],\\n)\\nax.scatter(0.05, -1.9, s=15., color=\\'black\\', marker=\\'*\\')\\n\\nplt.xlabel(r\\'$f_{\\\\mathrm{sub}}$\\')\\nplt.ylabel(r\\'$\\x08eta$\\')\\nplt.xlim(0.,0.2)\\nplt.ylim(-2.5,-1.5) \\nplt.xticks([0.,0.05,0.1,0.15,0.2])\\nplt.yticks([-1.6,-1.8,-2.0,-2.2,-2.4],[\"-0.6\", \"-0.8\", \"-1.0\", \"-1.2\", \"-1.4\"])\\n \\ncbar = fig.colorbar(\\n    pcm,\\n    cax=cax,\\n    #ticks=[0.,1.,2.,3.,4.,5.],\\n    #format=matplotlib.ticker.FuncFormatter(lambda x, _ : \"{:.0f}\".format(x**2))\\n)\\ncbar.set_label(r\\'$p(f_{\\\\mathrm{sub}}, \\x08eta | \\\\{x\\\\})$, uniform prior\\')\\ncax.yaxis.set_label_coords(4.3, 0.5)\\ncax.yaxis.set_tick_params(pad=1)\\n    \\nplt.savefig(\"../figures/posterior_uniform.pdf\")\\n'"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "pmin, pmax = 1.e-9, 0.1\n",
    "z = np.clip(posterior_flat, pmin, pmax)\n",
    "\n",
    "fig, (ax, cax) = paper_settings.figure(\n",
    "    cbar=True,\n",
    ")\n",
    "plt.sca(ax)\n",
    "\n",
    "pcm = ax.pcolormesh(\n",
    "    alpha_edges,\n",
    "    beta_edges,\n",
    "    z.reshape((resolution, resolution)),\n",
    "    norm=matplotlib.colors.LogNorm(vmin=pmin, vmax=pmax),\n",
    "    cmap=paper_settings.CMAP2\n",
    ")\n",
    "cs = ax.contour(\n",
    "    0.5 * (alpha_edges[1:] + alpha_edges[:-1]),\n",
    "    0.5 * (beta_edges[1:] + beta_edges[:-1]),\n",
    "    z.reshape((resolution, resolution)),\n",
    "    [threshold_flat],\n",
    "    colors=\"black\",\n",
    "    linewidths=1.,\n",
    "    linestyles=[\"-\"],\n",
    ")\n",
    "ax.scatter(0.05, -1.9, s=15., color='black', marker='*')\n",
    "\n",
    "plt.xlabel(r'$f_{\\mathrm{sub}}$')\n",
    "plt.ylabel(r'$\\beta$')\n",
    "plt.xlim(0.,0.2)\n",
    "plt.ylim(-2.5,-1.5) \n",
    "plt.xticks([0.,0.05,0.1,0.15,0.2])\n",
    "plt.yticks([-1.6,-1.8,-2.0,-2.2,-2.4],[\"-0.6\", \"-0.8\", \"-1.0\", \"-1.2\", \"-1.4\"])\n",
    " \n",
    "cbar = fig.colorbar(\n",
    "    pcm,\n",
    "    cax=cax,\n",
    "    #ticks=[0.,1.,2.,3.,4.,5.],\n",
    "    #format=matplotlib.ticker.FuncFormatter(lambda x, _ : \"{:.0f}\".format(x**2))\n",
    ")\n",
    "cbar.set_label(r'$p(f_{\\mathrm{sub}}, \\beta | \\{x\\})$, uniform prior')\n",
    "cax.yaxis.set_label_coords(4.3, 0.5)\n",
    "cax.yaxis.set_tick_params(pad=1)\n",
    "    \n",
    "plt.savefig(\"../figures/posterior_uniform.pdf\")\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQMAAADUCAYAAACYnUboAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHoFJREFUeJztnXt8VNW1x38rJLwhD0BERSERFK9aCBGwtVohiLX2Xh+ofdjaeyvBtlZrrQit9XG1RaAqVSkQuPXeqq081FofLRpQ8AGWEOpb1AQoWJ4hE94Jyaz7x9mTTCZnn3PmzDkzZ2bW9/M5n8ycffY+O8mZNXuvvdZvEzNDEAQhJ9UdEAQhGIgxEAQBgBgDQRAUYgwEIc0golIiKve6XTEGgpACiKg05v1kIionomkOql8DoMDrPokxEIQko77Vl0W9LwUAZq4CEIo1FCYs8aNfuX40KgiCHmauIqK6qFPXAHhFva4DUA6ghogmx1StYuaQX/0SYyAIqacAwL6o9/0AgJmXa64vB9CPiDw1DmIMBCHNYObZfrSbVsaAKIdBmevmIJC2rEuXPH1ZjnkZ5XTR34vM72UVkRoOHzM939LabFGnVVuWNnDrXmYeEHt6/KRJvK++XlvtnQ0bPgBwNOpUJTNXmlwaAlCkXhcA0DfqI2llDEA5yM3z3IkaGKw+8H17d3oW28v6DDI936N7obZOjuZerS1N2jqHDu8xPV8f2qqtc+TIfm1ZutDSXG/6C9bX1+Plt9/W1js+N7cbgLcAPM/Mz1vcYgmAMvW6GECVy64mROZ+zQqC3zCDLQ4dyjFYFnEQMnONOl8OIBR5n2zSa2SQIRSfPNb0fO8+J2nr9Cwq1pbl9Ip/tBQ+bP6NHT7SqL+PZjRhNb35vOlD8/uEWyx6lx4wgJZwOP56hmNwecw5s+lDUhFjIAgJELZO+m1k5ookdSVhZJogCC5hBlrDYe0BIJ+IKono66nuqxNkZCAIruHIhz4jEGMgCC5hAGFrcSCZJghCVsBAS2tYe0CmCYKQHTCslxCRZiMDMQY+cVrJeG1Z0clfND1PQ4e4uldOL82SXxf9kh/2HDavs32LtkrXo/mm55v31WrrZMISohYGWmyWE9IJmSYIgksYQBisPSDTBEHIDhhsF3Qk04RoVMhlCECpWbaVEnIoBixTNgUhkNisJqQVvk4THCq4zFBGoNiBwosgBAZmIxxZd6QbfvsMroExKgDaFVzaUKOG9YCRo52qBA1BcEuYWXvAB58BERUrQdRpRKRPWHGB38bAVMElinNgKLaUOhSCFITAwGC0hvUHlM/ALH05AUHUUhhfrFUAYmXREiIIDsR6Zq5Rf4jJsX4DIqoAoJwwwVv8GDHsItPzhV+4VFun6/D+pud7FXbV1qEc/TKhRqcEB3brtQlau5v/6ylX34fWlqOm5480pb9mgRsi04R4UanKCwGUqPdt0+nIN79ulBz5fKg2PPWxeW4MiGghDNWW9bBXcKmHYeWgrj0H5qmdlQBAObmZ460RMgI3DsREBVGVIahi5jp4iOfGgJmnRl4ri9dJwYWICpSQ43K0D3UKoPwHgpAOGHoGnnw/ORZEVYbgdgB1RPSKlytwvk4T1PC/zETBZSWA0cxcR0QhZQH7+SX0KAj+wElfWlQrc77IovnuMzBTcGHm0SblEmMgpBVsH47cn4iqo96LIKogZCo2iUrHANRABFEFIbNhAC3M2kOHCKKmMQP7n6ot0y0hFo47QVunaEB30/N9++ql0g8f0e8/0Npq/uAda9IvezVvM88mDDcf0dapr//E9Pz+A7u1dTIdl6sJIogqCJkEM2dUCrMYA0FwicieCYLQhk04sugZCEI24CAcWUYGgpANGNME/QEZGQhCttA2HdCRViMDMQYOKCn9rrasYIz5EmLxsL7aOv36mmcG5lkImDa36Iejn2w+YF7nsIUY6X7zOnu3rdVW2bH7Y317WQizKB0JgqBoCbP2gEwTBCE7YNiGI8s0QRCyAtk3QRAEwJA9S7YGop/IyEAQEqDVWvXM82kCERXAyHAshZESHbKp4hgZGaSIcDiMD159BSsXz8PRg+aefSHYRHwGusOKBARRywBUw9BA8FQdWUYGUQwedLbp+dyz9X/zYWeY7z84bmiR6fnWlhbUrHgJTz36CHr06oXBJSVYcN0V+OF/34fi8/X7M360U28wSKOIevTjvdo6O999xvR83Vb90qIQA8MuzsCUBAVRq9TooMDrVGcxBkniWHMz1vz5GTy74HcYMGgQbrj7Hoz+8vkgIrz79jrMnT4Ng5YtxXV33I2igcenuruCAxi20wTzegkIogK4mpkriaiGiKZ5KRUoxiAJrH72afzpoQdwUkkJfnj/HJz75fM6lJ89dhwW/PVl/H7uXMy47FJM/vHNmPCNbyMnR2ZxQcduOuCQeARRq9XIohhBl0oXOrLiycfxl8WVuPXheRg2cpT2uq7du2PyTbdg3Fe/hvnTf4aDjY24/Ac3JrGnQrwwuxsZJHZP/1SQxBj4yJrnnsWzC36He/64FAMHD3ZU56RhwzF6wkQ0HzXfsEQIFjbhyCKIKgB/f2UFHp/1a9z5f390bAgiHGhowIATT/SpZ4J3MMLWDkQRRM12tmxYh4V3/BwzFj2GwcOGxV3/QMM+9Ck0X40QggNbpC9b2QgRRE0DThz2VdPzw8ti94ttJ3YJcVPNBiya80s89IfHMWrcONM6B5uPmZ4/2noYALC/oQHd8vNxtLVdBLW+sVnbhy2v7zQ9//na/9XW2bp9g7ZMcI7NyMCUrBVEVdYvBKDUbBkkqrw4CH+QRNj8wfuY88Op+PFvHtQaAiccCDWgT0Ghhz0T/CKDMpj9nSZEB1MACJlEXZUCqFPldbHl6cTndbWYWfFfuP7uezHyyxck1NbBhgb0LhRjEHSYbTUQtVuyBxG/fQbXwPjWB9qDKWKZpX4Wp2qulCiHDxzA/RXfxzd+8jOMu9h8quGUpiNHcLChAX2L9FMTIThwmLUHJFGpA6bBFBHUxqx1RNQAYIrPffEFZsaiu+7AmePOxfirrk64vfdeX42SkSPRrWdPD3on+AnD1meQVnoGKV1NUDHWIQAzASwiok5JAERUoaKuqoM4QXvtmeXY+vFH+N4v7vSkveqX/4ZzJl3iSVuCz6igI92BbB8ZENFCGAEU62EfTFEBYCYzh1Ss9mQAHZyMyqlYCQCUkxsoa9CwfQtemjUT9zz5FLr16JFwe01HjuC9N9bg2x4ZFsFvbLMT02pk4LkxYOapkdfKIdgpmIKICmLzsJl5ORGl9A+XO+oM0/NDT+g8ZD/W1ITnbv8lbrnzTkw8d2yn8r7du2nvs/uQeTbhyy+8jIHDz8DO5m7YuetQh7IPqv6lbW/b3x4xP7/jXW0dIXEcTBPyiagS9kFHgcBXn4HyCZSZBFOsBDCamWcT0TQ1KihKp6XFlx55AP1PHoLJ133PszY/fu1ljLjgIs/aE3yGbZcWkzYyIKIrAUwEUAjDT0cw7NUrzGyerx6D73EGZh9wZh4d9dqzFMxk8cGaVfjozdW4+Q/LtFoC8XL0yBHUrX8LF908w5P2hOQQ1uyAnSyIaBSA0QCqmPlpk/KhylDUMvM/rNqSCMQ4aWluxrOz78W1v34QPfro90aIl5eXPIWTzhyJnvkSX5AuGEpHqe4FQsy8WFfIzJsBbCaioXYNSW5CnGxc8QIGnXoahpytT0eOl9DevXjytw/hwoqfeNamkATY8BnoDvi4mqAc9ZEPu31XHVwnxiAODjc2YtVjlbjg2v/0tN3/uf/XKL/yKgwYeqqn7Qr+w6w/YBGBmIAGYqe6JuXz4/w1AMg0wTGtLcfwxM9vwRnnj8epZZ1XD9zyYXU1atasxqJVq7F+l2gYpBMOVhNMSUQDUcXihNQBIsoHsEh1J+I0nE5EC2A4EyPnpjDzfqt+ZZ0x6N69j7as8JTepufHDe6PObfdioGF+bh3zv3o0qVLW9nAvuZ1tjc0au+z85DxoW9tacHcX8zAVT+bjv2Ui7erdmjr1D2r97PKHogpgu2XE8yrJaSBWAxjtaBYGYZRzGwW+npD9BsimgBjFU9L1hkDNyxZuAAf/WMj5v35+Q6GIFHeWf0qcrt2xdhLLvWsTSG5hL2RPXOsgQigTR1ZXdNpBcEMZrY0BIAYA1t2vPsmXv/LAsx//iX07G0+CnBLzcoqnPu1f/dseVJIMuxumpDwbY2AvYmx54mor91UwApxINrw3jPzcOe8+RjosQzZpzUb8M6aV1F20SRP2xWSR2Rp0cKB2D+SV6MOXQCSVxqIM4hoJGDEH0ReO0VGBhYcqt+BlqYj+MK4cz1td9fWLXj0Jzei4v7foFD2SEhrAqaBWA3Dl1DHzBuJSL8rjwkyMrBgz6YaDBhe6mgYf/DgQRzXpQsOHTpkeV1jQwMe+sEUXPajm3DWeed71VUhVdgMDczwUQOxGIa/YTYRrYCxH6NjxBhYsGfTBgw4fbTtdXNnzsT4UuPvfuGoUZg7c6bpdc1NTbj1u9/ByAsn4MJrvulpX4UUwIYDUXdoqzEvZ+bCaAchM1cyc1WC+Tl1zLyImW9g5kkwViYck3XThMGD9NOo4SPa901kZqyo/QcuvfVWDC6wDzveUlvb4een+zouLTIz5k27Fbl9+qLr2d/Es69s69TG1qce0rYvy4fBI7LxqgVJTWFm5qeJaAgzb1E5CyXx1JeRgYa9mz9D1569kD/wBNtrb7jlFsv3APD0vEfwr7pa3DjnQZBsm5YxcFh/IAXiJsy8Rf3cyMxz4qkrT6WGLTVvY0jpGEfX5uXl4Y6ZM7Fpzx7cMXMm8vLyOpS/8fxzeO2ZZZi2YLEnIihCQLDYjt2jPRhtceokJKIr7K7JummCU3Z/tgmDv2DvLwCALl264KZpRkh55Gc0Lz72e0y9734UDBjgaR+F1MIAwq2WlyRjmrCBiJaiPRy5HobI8OyocwwHGqNiDDQcbKhHr6L+nrTVuHcPBg0Z4klbQrBI1gjA4v6NADqEIxPRAk2IsiUyTdBwuKEevT0wBuFwGKH6vcjv741hEQIEB1Yq/XY3lWRkoOHQvnr09GC/w0ONjejesyfyuuo1EYX0ha1zE1IiiKpGC3GTdcYgv1CvGfBvJxgZjeHWVhw90IiyEScjNy8Pvbp11dZpOtZien5drZF3sndrLbr1LWx7DwC7/vSiaR3Z/zD9sJkm+C6IqiTNIt9aER9B5HWEfU50ELPOGDjhYCiEHr37IDdmVcANhxvq0bNQdkfKRJhTPzJwmrXoBDEGJoT27vFsjn9o3170KhR/QabCKchajEbpFJQDWM/Mz0SSk+zET80QB6IJod27UDDguITbOXpgP97+42IMHXOeB70SAgfrnYdJNBLMzDNgiJ7Oh5GfcA4RzSeiIfE0JMbAhNeWLcWIMYlJm4VbW/DizOkYUvZFjBgv26VlKjZ5Sp6vJihZtGVR6dCFRj94I4zchGdUfsIPYL7RsZakTBNsNN0mw8jnLg3CHgobX12FbZ9swg1zHkyondWVRtjx+VNE8TiTsRkBaH0GsZ+JOD8HUyI7kql8hAkwHIYLY65zpJwcwfeRgUrLXKYpaxOCBBCyU331m6bDh/H4fffgujvvQddu7pcCn1i0CFtr1uGS6TOR00XcMpmKA3ETU2I/E3F+DiL6h5OV/BmYeSUzrwJQQkTjiWgCEU0BkG/RTie0T6qyNhMB/D0Rx4SJ+GM0pkKQ8bQfL13PGKctq1q0ACPHjsXFkzpucaZbPgSAt/7ZcQ/E99etxcN334Xh4+/GByvqYSZaU13zWHydFoIJu3MgJiiIWqbqA8bGxW2jCDVViLBS7aY0AUC9k8+t1dcWM/N0JZ80X3W2HxFNBTArkh2VIKZCkKng0L6tWPH6EiyuetV1Gzv/uRUP//Qm/PiB32LNG5KQlBV44yd0LIhKRHVRGxqbCaa2d03tpuS0E1bGoM0xoWSU2oIWiOh6ANotndKRT9cuwvdvm4Yil8lEhw/sx+yp1+PKH92Es774Jax5w9cBjhAQkr20yMyREYXnD5jWGLh1TKhtn4pgrHvaOUJshSCV11Q5YfxxcYRbmnFgz2f42reudd3GiiefwCmnj8Ckb3/Hw54JgYYBtt54tT8RVUe9r9QoGXkliJoQlt6tiNa6mipEdmcphkVnmXmq3U2JqEB5Q22FINUfrxIAKCfXFzPcfLQRed37Isel6Eg4HMaqpU/h5rmPeNwzIdjYbqKSbEHUNtROS9Uw/H4FTnwGjp5+pZqySnktFwHYqDyWtlLMseKPipWqXa+EIBPi2JFGdO1R4Lr++2vfRM8+fVBy1tke9koIPPZZi6b4KIja3jUjWamMmbc4dfq7WveKxzGhnCDLY86NjnqdiABk3PQa3vlDf6i1BT3690f3XPM/R219g7a9Zau2YW3l73HcyIux/NXtbeff+kPSk9WEVOBirKr5THj+OYg3e1EiEAE0HWxA196FruoebazH7o+rcfKYThvcCBkOw3ZkoN2F2W+I6EoiWuBE7iyCGAMAzQcb0K23u2nClrV/xYmjLkBej14e90pIC9jiSJ24CQBsYOYbEIdcuoTHAWg6GEKP/PiXFI8cPoTNbzyPsd+/y4deCYGHAaQ4azGWKH2DYiKqhbE64Z/PINPofdzJqK/daH9hFC3HjuGuiuvRf9hIFA4Z4VPPhKDjNjfBR0Joj0GIy1LJNAHA8Wedj72b1uPgfmcb2IbDYdx/y83Izc3D6Gtvk12UsxVmY2SgO1Kzb8JK5eCfDGA6OscIaRFjAKBrz77oVzIKr734gu21zIx5d9+JXZ9/jrvmL5REpGzHOlMpZQ5EAMuVz2C60wpZ9yT37GsuZTb0S5PwwtIlGPsfl3cqW/3JXgCGNuLaJxaidt1qXD1nMVZ81oCdS9/wtb9CwAmYzyCKUhXHUATA0RxYRgaK488ch7oPP0D9zp2m5dve3YAnbvwWPn9/I66471F0790nyT0UAocKR9YdSO1qQghG+PMMpxWybmSgo0teN4yb9FWsef45XD6lPaJ69+fb8cKv7sKOTe/hgik/xbDzJoiPQGjHemCQEql0RQGACUQEJ8rIgIwMOvCVyy5H1bIlaDl2DACw+rk/46dfvwT9TinG9yqfxvAvl4shEDpi7UBMJTUAxiCOOAMxBlGcOXYcjjvxJCx59GFs+egjLL73bvzqT8tw7rVTkddd9AmEGKych5ycaQIRXaERPi1FxwQoW2SaEAUR4abZv8EN4y/A0kd+i1senItTTjsNW5QDURA60Wq5cYIv0wSV1l8HQxBlDFRcARGNV/Jn0fspOA6gEWMQQ+GA4/DLxY9h66aP8ZXLHId1C9mKy41X3QqiKkNQFRE5ISIGMFtNX/OVClINgGpmdhY4o8g6YxDacURbVlt/EADQ49QzcfqpZ7a9f+fN3do66954wNsOCukDMxC23lLJDJWqvBBAiXrfJoiqpNC1auIARgPYF/nQK93Dq1U7V8IYMZQBuFoZihCAhU5kCrPOGAiCl5CLkUGCgqgh9bMIwFRE7bgcMzVY1NZHopkAbJcYxRgIQiJ4s2oQjyDqTLRvjuI01NiRcpIYA0FwS2TjhGTe0pALbDMSRNTXzjfAzCudXCdLi4LgGuUz0B1KEDXq0K0sJCKIeo6dgInyJdguMcrIQBASwdqB6LsgqvrWzyei26AcktHFAGoBLHIigSbGQBDc4nKaEC2IyszLmbmGiMrcCqKqD/qcuDsSQ9YZgy5d9TOjD7eYT6kOr3zNp94I6Q0D4db4a/kkiKr2V7wKaq8TpzkJEcRnIAiJ4GbnVf+oY+aLYMQdkNr5zDFZNzIQBM9gBluPDJKdtZhPRCPVPglPE9H4eCqLMRCERLAeAeQTUSXsHYheMQbGtuw/j7o/YOyoZDtlSIoxsAqvjFpuKWHm282uEYRgYuszSPbIYAkMB+QcACCioTAClCoA2BoD330GykO6zKKsSjlPitV7QUgPIqsJKUxh7tAdYxvEzVHvN6vtEB19rnwfGZjEYUdTrI5KGDHZxX73RxA8xTrOIJVKR2043WYtpT6DmOWUiBiDrzTtPaot2/qh+fn11YvMC4QshwGOP2sxqARiaTEqHTMluzALgluYW7UHUiuIGjeejwyIaCGMOOv1ViINMZTrnIfKwaiGWoGwXYJgwIFzICaE58aAmafaXUNEBSr7CkRUETEaRFTOzB3istVUohIAKCc35SqTghANyzTBOdFx2FGnV6qycgCziKiWiBr87osgeIsaGegOH6YJRDSZiJYR0UIimuZVu0ByVhPM4rBHq59VAAr97oMg+IGxguhuNcGtBiIM39py9UVa7abfOmQSLgiusR0ZmBIbexOtgQggFHlvekclhKpehzz5NRRZF47M2/Wy5011+mVHQTDFhc8gEQ1EZg6p8443R3FK1hkDQfAM+0QlpzjWQFTEJYDiFDEGgpAAqVhNiGPJPi7EGAiCazgSXKSjPxFFO/kqNSImiWggeoYYA0FIBOuRge8aiF4iqwmC4BJmRjjcqj10xMbeRJYY3WogeoWMDAQhAWymCaZxBn5pICZK1hmD8L4d2rI929clsSdC+mObtZhspaOEyDpjIAheErYeGaQVYgwEwTW2I4O0yloUB6IguMTYkd3SgZjdegaCkD0wABkZCIIAGRkIggAg03wGWWcMtnyiX+HZuefTJPZESHvYdjVBlhYFIRtgsGtxkyAixkAQEsCjFOZAIMZAENyichMyBVlNEATXRJYWdYcvgqgFRFRKROVEVOBVu4CMDATBNQzYjQz8EEQtAzARwHoY6c6eZTjKyEAQEoCZtYeORARRYSgilwOYBY91ELNuZHCsRURPBa+wVToyr5WAICqAq5l5tDIYFQA8k0DLOmMgCJ5hv2+CUxwLohJRnRpZ7EOMJkKiJMUYxM6PNNdM80voURD8gMFJT2GO3X7QS5KxvVqH+ZHFNRP97osgeI2Nz6A/EVVHHboApOwQRDWZHwlChmC7b4IIosaDmkKk5JcXhERhhLWHDhFE1VNkf4kgBA92GYGYNYKoRLQQxgd8vZ1D0MmoQM2z1Fwr8YFMfcO2hNsQhAhug46CiOfGgJmn2l1DRAVqB9liIiqGYTyKzFYdlMWsBADKydVHcghC0mEwLB/JtEphTsZqQof5kWIlYAyXotZSPY2zFgS/MTQQW7RHukFWYZNBg3JyOTdPbIaQXFqa6zcwc1ns+by8ntyv3+naert2bTStF1RSvpogCOlMmFu1B0QDURCyBVE6EgQB7fsmZAoyTRAE1yh1ZN0h0wRByB5sEpVkmiAI2UI6rcbZIcZAEFziNhw5qIjPQBASIMxh7QF/BFGLiWgyEVV4LYgqxkAQXGMsLeoOKJ+BWShyrM6h+oCXE9E0m5tOhZHiXIW2nB1vEGMgCC4x1JHD2kNHgoKotTByeYqh5NG8QnwGguAWtt14VVPNvSAqM1cqY1EAQy7dM8QYCEIChJMviFoKY1RQ5LUGghgDQUgAj9SR47gf18DDjVOiSS9jwK0HW5rrN6W6Gxb0B7A31Z3QIH1zz2lmJ5lbVhxr2tvfol53IqqOel+p+TbPDkFUj9kU5JRQIqoOav+kb+6J+UC3wcwXe3QLEUQVhGxEBFEFQQAQXEHUdBsZpPwPZkOQ+yd9c0/Q++cJaSV7JuhxsqV3AtuAp6p/s5j5diKqCMI3Z6YTqJGBXUimWXkcYZyp6t8s9dO3VFYnEWwJRr0lvX+KCiKqhcdbj5vc2+7/WqGOWU7rpCOBMQZ2D4xZeZAeaIvyZDzQ18D41gXaI9g6oPoVG/VmWSfF/QOAKcxc4ueOWw7+r+VQkX8wpP3Lk/ncJZPAGAPYPzBm5UF6oHXlvj/Q0ESw+VDHLW7vVZyEb1+7/2tx1Lk69T6Zz13SCJIxsHtgzMqD9EDrypPxQGckzDxbGdF+6hvaDyz/r8wcHShUCqDark66EiRjkJEk6YF2E8GWzKi3uO+l5uiRRJ16GN/IKUNNBWpSFQOQDIJkDOweGLPyID3QncqT+EAviWq7LYLNRvzCtI5PuOlfdVSfStR7P3D6DJUz8+1x1kkrgmQM7B4Ys/IgPdBm5Ul5oC0i2FZGrkll1FsC/btava/1sX+2hkotbc6O+h2S+dwljUDFGajltzoAxZF5GhFtYObRFuWdzgWwf/vUOV/X8gV3WP1fo5Y898EYDVyl9AiS9twli0AZA0EQUkeQpgmCIKQQMQaCIAAQYyAIgkKMQQaiQrVnmQhqCoIW0TPITGYAuN32KkGIQoxBhsLMvmb6CZmHTBMyDDU1KPAx9FnIUMQYZB51AF6xy5JUyVPiUxDaEGOQeZTBma5+Hdrj6wVBfAYZyGgAS4G2KUMNjNTbEIykmhoAk2EIcpYQUTEMKbJOu/cI2YWMDDIQZo4W3piqXlfDiKOvg5E0BRgJQHUAJia5i0IAEWOQeYSiXpfByLCLGIh+UXv1AYajsRTtm34KWYwkKmUIavVgIoCFsqwouEGMgSAIAGSaIAiCQoyBIAgAxBgIgqAQYyAIAgAxBoIgKMQYCIIAQIyBIAiK/wdiqcbxZzwJwgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 245.424x204.52 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pmin, pmax =1.e-9, 0.1\n",
    "z = np.clip(posterior_narrow, pmin, pmax)\n",
    "\n",
    "fig, (ax, cax) = paper_settings.figure(\n",
    "    cbar=True,\n",
    ")\n",
    "plt.sca(ax)\n",
    "\n",
    "pcm = ax.imshow(\n",
    "    z.reshape((resolution, resolution)),\n",
    "    extent=(alpha_edges[0],alpha_edges[-1],beta_edges[0], beta_edges[-1]),\n",
    "    cmap=paper_settings.CMAP2,\n",
    "    origin=\"lower\",\n",
    "    aspect=\"auto\",\n",
    "    interpolation=\"nearest\",\n",
    "    norm=matplotlib.colors.LogNorm(vmin=pmin, vmax=pmax),\n",
    ")\n",
    "# pcm = ax.pcolormesh(\n",
    "#     alpha_edges,\n",
    "#     beta_edges,\n",
    "#     z.reshape((resolution, resolution)),\n",
    "#     norm=matplotlib.colors.LogNorm(vmin=pmin, vmax=pmax),\n",
    "#     cmap=paper_settings.CMAP2\n",
    "# )\n",
    "cs = ax.contour(\n",
    "    0.5 * (alpha_edges[1:] + alpha_edges[:-1]),\n",
    "    0.5 * (beta_edges[1:] + beta_edges[:-1]),\n",
    "    z.reshape((resolution, resolution)),\n",
    "    [threshold_narrow],\n",
    "    colors=\"black\",\n",
    "    linewidths=1.,\n",
    "    linestyles=[\"-\"],\n",
    ")\n",
    "ax.scatter(0.05, -1.9, s=15., color='black', marker='*')\n",
    "\n",
    "plt.xlabel(r'$f_{\\mathrm{sub}}$')\n",
    "plt.ylabel(r'$\\beta$')\n",
    "plt.xlim(0.,0.2)\n",
    "plt.ylim(-2.5,-1.5) \n",
    "plt.xticks([0.,0.05,0.1,0.15,0.2])\n",
    "plt.yticks([-1.6,-1.8,-2.0,-2.2,-2.4],[\"-0.6\", \"-0.8\", \"-1.0\", \"-1.2\", \"-1.4\"])\n",
    " \n",
    "cbar = fig.colorbar(\n",
    "    pcm,\n",
    "    cax=cax,\n",
    "    #ticks=[0.001,0.003,0.01,0.03],\n",
    "    #format=matplotlib.ticker.FuncFormatter(lambda x, _ : \"{:.0f}\".format(x**2))\n",
    ")\n",
    "cbar.set_label(r'$p(f_{\\mathrm{sub}}, \\beta | \\{x\\})$')\n",
    "cax.yaxis.set_label_coords(4.3, 0.5)\n",
    "cax.yaxis.set_tick_params(pad=1)\n",
    "    \n",
    "plt.savefig(\"../figures/posterior.pdf\", dpi=300)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Below here not updated! There be dragons"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1D projections"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "def project(posterior):\n",
    "    return np.sum(posterior.reshape((25,25)), axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "posteriors_1d = [project(posterior) for posterior in posteriors]\n",
    "posteriors_1d_flatprior = [project(posterior) for posterior in posteriors_flatprior]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "def prep(edges, histo):\n",
    "    edges_ = np.copy(edges)\n",
    "    histo_ = np.copy(histo)\n",
    "    edges_ = np.repeat(edges_, 2)[1:-1]\n",
    "    histo_ = np.repeat(histo_, 2)\n",
    "    return edges_, histo_\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xl4lNXZ+PHvmcnsWScbWUAChH0VXEDcWq1oXau17tbWn/VVq7a1r9Vaba12dWvfqlVbte+r4loUWxVQcasbIMgiyhogQJJJZrJOZjLL+f3xTCYJ2SaBSQK5P9c1V+bZ71zCzfE859xHaa0RQggx8EyDHYAQQgxXkoCFEGKQSAIWQohBIglYCCEGiSRgIYQYJJKAhRBikEgCFkKIQSIJWAghBokkYCGEGCQpgx3AgZKTk6NHjx492GEIIQSrVq2q1lrn9nbeIZOAR48ezcqVKwc7DCGEQCm1I5HzpAtCCCEGiSRgIYQYJJKAhRBikCS1D1gptQD4E2AG/qa1/t0+x68GrgUiQCNwldb6i9ixW4Dvx45dr7VeksxYhRD7LxQKUV5eTiAQGOxQBoTdbqe4uBiLxdKv65OWgJVSZuBB4GSgHFihlFrcmmBjntFa/zV2/pnAfcACpdRk4AJgClAIvKmUGq+1jiQrXiHE/isvLyctLY3Ro0ejlBrscJJKa01NTQ3l5eWUlJT06x7J7II4Etiitd6mtW4BngXOan+C1rq+3aYLaK0OfxbwrNY6qLXeDmyJ3U8IMYQFAgGys7MP+eQLoJQiOzt7v1r7yeyCKAJ2tdsuB47a9ySl1LXAjwEr8LV21368z7VFyQlTCHEgDYfk22p/f9dktoC7iqzT+kda6we11mOBm4Hb+nKtUuoqpdRKpdRKj8ezX8EKIcRAS2YCLgdGttsuBvb0cP6zwNl9uVZr/ajWeo7Wek5ubq+TToQQw4DZbGbmzJnxT1lZGStXruT6668f7NA6SWYXxAqgVClVAuzGeKl2UfsTlFKlWuvNsc1vAq3fFwPPKKXuw3gJVwp8msRYhRCHCIfDwZo1azrsGz16NHPmzBmkiLqXtASstQ4rpa4DlmAMQ3tca71BKXUnsFJrvRi4Til1EhACfMDlsWs3KKWeB74AwsC1MgJCiIPLr17dwBd76ns/sQ8mF6ZzxxlT+nzdO++8wz333MO//vUvrr/+enJycrj99ttZsmQJd999N++88w4m08BPi0jqOGCt9WvAa/vsu73d9xt6uPZu4O7kRSdarahYgUIxZ8TQayEI0VfNzc3MnDkTgJKSEhYtWtTh+O9+9zuOOOIIjj32WK6//npee+21QUm+cAgV4xH9d+/Ke4nqKM+f8fxghyIOIf1pqR4IXXVBtOd0Onnsscc47rjjuP/++xk7duwARteRJGBBlb+KumAd4WiYFJP8kRCHvnXr1pGdnc2ePT2NC0g+qQUxzIWjYWoCNbREW9jVsKv3C4Q4yO3YsYN7772X1atX8/rrr/PJJ58MWiySgIc5b8BLVEcB2Ozb3MvZQhzctNZ8//vf55577qGwsJC///3vXHnllYNWu0L+f3OY8/jbJrBsqd3CN/jGIEYjxP5rbGzstO+EE07ghBNOAODNN9+M7589ezbr1q0bqNA6kRbwMFflrwLApEzSAhZigEkCHuY8zUYLeGrOVDbXSgIWYiBJAh7mqvxVmJSJowuOZmf9TgLh4VHHVYihQBLwMOdp9pBtz2aieyIazda6rd2eq7Vm8ed7eG+TFD4S4kCQl3DDXJW/ilxnLuMyxwGwxbeFKdmdB9Dv8vq5ddE63t9czaSCdI4bL8WPhNhfkoCHOY/fQ4GrgFFpo7CarJ1exEWjmqc+2cHvXv8SBYzNdeFrahmcYIU4xEgXxDDnafaQ68zFbDIzNnMsW2q3xI9tr27igkc/5vZXNjD7sCyW/Og4Tp48Aq+/Ba07lWcWYkhQSnHppZfGt8PhMLm5uZx++umDGFXXpAU8jIUiIbwBL7lOozuhNKuUj/d8TCSq+fsH27h36SZsKSb+cN50vj27GKUUbpeFlnAUf0sEl03++Iihx+VysX79epqbm3E4HCxbtoyioqG5oI78DRrGqpurAchz5AEwLnMci7cu5qyHl7F+V4iTJuVz9zlTyU+3x6/JdFoB8Da1SAIWPXv9Z1BxgCc5jJgGp/6u19NOPfVU/v3vf3PeeeexcOFCLrzwQt5//30APv30U2688cZ4gn7iiSeYMGECGzZs4IorrqClpYVoNMpLL71EYWEh559/PuXl5UQiEX7xi1/wne9854D9OtIFMYxVNRuTMHKduYQiUTbucAJQ3rSNP184i8cum90h+QK4YwnY55d+YDF0XXDBBTz77LMEAgHWrl3LUUe1LUc5ceJE3nvvPVavXs2dd97JrbfeCsBf//pXbrjhBtasWcPKlSspLi7mjTfeoLCwkM8//5z169ezYMGCAxqnNGGGsdZpyM3NLs78y3/40hMhtRSuO8XFmTMKu7wmy9XWAhaiRwm0VJNl+vTplJWVsXDhQk477bQOx+rq6rj88svZvHkzSilCoRAAc+fO5e6776a8vJxvfetblJaWMm3aNG666SZuvvlmTj/9dI499tgDGqe0gIex1mnIjy2vpqKumYcvOJE0axp7/WXdXuN2SQtYHBzOPPNMbrrpJi688MIO+3/xi19w4oknsn79el599dV4IZ6LLrqIxYsX43A4OOWUU3j77bcZP348q1atYtq0adxyyy3ceeedBzRGaQEPY55mDykqhYpaM8eWZrNgagELd5X2OCXZHe8DDg1UmEL0y/e+9z0yMjKYNm0a77zzTnx/XV1d/KXck08+Gd+/bds2xowZw/XXX8+2bdtYu3YtEydOxO12c8kll5Camtrh/ANBWsDDWJW/ihxnDjUNYXJSbYAxEmKLb0u3w8zS7CmYTUrGAoshr7i4mBtu6Lzq2X//939zyy23cMwxxxCJtC01+dxzzzF16lRmzpzJl19+yWWXXca6des48sgjmTlzJnfffTe33XbbAY1RWsDDmMfvIduey+ZQhJw0o2VbmllKQ6iBSn8lI1wjOl1jMimynBbpghBDVm/lKOfOncumTZvix379618DcMstt3DLLbd0uO6UU07hlFNOSVqs0gIexjzNHtItbgByYy3gcVnGlOSeSlNmOa2SgIU4ACQBD2NV/iqcJiMB56TFEnCsJkRP/cBZLquMghDiAJAEPEwFwgHqW+qxkAm0tYAzbBnkOfPY4tvS7bVupxWfvIQTYr9JAh6mWguxm6IZAPGXcGC8iOu5BWzBK10QQuw3ScDDVOskjEgoDYDsVGv8WGlmKdtqtxGOhru8NstpxdckBXmE2F+SgIep1mnILcFUMp0WLOa2PwqlWaW0RFvY2bCzy2vdLivhqKYh2HWCFkIkRhLwMNXaAm5qcnXofgCjBQzdj4TIaq0HIS/ixBC1aNEilFJ8+eWXAJSVlTF16lQAVq5cyfXXXz+Y4cVJAh6mPH4PVpOV2qYUctp1PwCUZJRgUqYOtYHbc0s9CDHELVy4kPnz5/Pss892OjZnzhz+/Oc/D0JUnclEjGGqqtlYiqimooWpRRkdjtlT7IxKG9V9C1jqQYgE/P7T3/Ol98sDes+J7oncfOTNPZ7T2NjIf/7zH5YvX86ZZ57JL3/5yw7H33nnHe655x7+9a9/4fF4uOiii6ipqeGII47gjTfeYNWqVTQ2NnLqqacyf/58PvzwQ4qKinjllVdwOBwH9PeRFvAw5fF7yHPmUd3Y0qkLAmJTkrtrAce7IGQomhh6Xn75ZRYsWMD48eNxu9189tln3Z77q1/9iq997Wt89tlnnHPOOezc2fbeY/PmzVx77bVs2LCBzMxMXnrppQMeq7SAh6kqfxXjMktpDIbJTesiAWeW8uaON2kON+NI6fivfpbLAkgLWPSst5ZqsixcuJAbb7wRMOoCL1y4kGuvvbbLcz/44AMWLVoEwIIFC8jKyoofKykpYebMmQDMnj2bsrKyAx6rJOBhytPsYUa2UaR63z5gMKYkazTb6rZ1WiU51ZaCxaykD1gMOTU1Nbz99tusX78epRSRSASlFNdcc02X5/c0lNJma2uYmM1mmpubD3i80gUxDDWFmmgKNWFXxr/2XXZB9DASQilFptSDEEPQiy++yGWXXcaOHTsoKytj165dlJSUUF5e3uX58+fP5/nnnwdg6dKl+Hy+gQxXEvBw1DoEzayNachdJeCRaSOxmW3dvohzO6UehBh6Fi5cyDnnnNNh37nnnstvfvObLs+/4447WLp0KYcffjivv/46BQUFpKWlDUSogHRBDEut05CJpANthXjaM5vMjMkY0+2LuCyXRV7CiSGnfeH1Vtdff32Hcb/tS1NmZGSwZMkSUlJS+Oijj1i+fDk2m43Ro0ezfv36+DU33XRTUuKVBDwMtS5FFA6mAn6yXZ37gMEYCfHRno+6POZ2WdlU2bnuqhAHk507d3L++ecTjUaxWq089thjA/p8ScDDUHwxzoCLNHsLdou5y/NKM0tZvHUxtYFaMu2ZHY611oMQ4mBWWlrK6tWrB+350gc8DFU1V+FIcVDbaI6XoexKaVbsRVwXldHcLuMlXDQqBXmE6C9JwMNQb5MwWrUWZ++qHzjLaSWqoT4g/cBC9Jck4GGoyl9FriOX6sZgfC24ruQ580i3pnc5EqJteXpJwEL0V1ITsFJqgVLqK6XUFqXUz7o4/mOl1BdKqbVKqbeUUoe1OxZRSq2JfRYnM87hxtPsIdeZi6cx2GMLWCnFuMxxXbeApSCPEPstaQlYKWUGHgROBSYDFyqlJu9z2mpgjtZ6OvAi8Id2x5q11jNjnzOTFedwo7XG4/fgtuXSEAj3mICh+2Xq3VKSUgxRZrOZmTNnxj+9TSEePXo01dXVAKSmpg5AhG2SOQriSGCL1nobgFLqWeAs4IvWE7TWy9ud/zFwSRLjEUBDqIFAJIDT3P0suPbGZ43nudBzVDRVUJBaEN+f6TTqQcjSRGKocTgcrFmzZrDDSEgyE3ARsKvddjlwVA/nfx94vd22XSm1EggDv9Nav7zvBUqpq4CrAEaNGrXfAQ8HrUPQrLTOguu+Dxg6rpLcPgHH+4ClBSy6UfGb3xDceGDLUdomTWTErbf2+bonn3ySlStX8pe//AWA008/nZtuuik+IWOwJLMPWHWxr8sxS0qpS4A5wB/b7R6ltZ4DXAQ8oJQa2+lmWj+qtZ6jtZ6Tm5t7IGI+5LVOwjBFYotxdjELrr1xWbEEvM+LOKfVjDXFJC1gMeQ0NzfHux/2nZY81CSzBVwOjGy3XQzs2fckpdRJwM+B47XWwdb9Wus9sZ/blFLvALOArUmMd1honYYcCacBNT2OAwZIt6aT78zv9CJOKRVbnl4SsOhaf1qqB8LB1AWRzBbwCqBUKVWilLICFwAdRjMopWYBjwBnaq2r2u3PUkrZYt9zgGNo13cs+q+1BRwMGC8buqoFvK/SrNIuh6Jluax4pR6EOAikpKQQjUbj24FAYBCjaZO0BKy1DgPXAUuAjcDzWusNSqk7lVKtoxr+CKQCL+wz3GwSsFIp9TmwHKMPWBLwAeDxe0izpFHbpEi1pXQ7Dbm90sxSttV1Xqbe7bJISUpxUBg9ejRr1qwhGo2ya9cuPv3008EOCUhyLQit9WvAa/vsu73d95O6ue5DYFoyYxuuWscAVzcGe30B16o0q5RQNMTO+p2MyRwT35/ltPLFnvpkhSrEAXPMMcdQUlLCtGnTmDp1KocffvhghwRIMZ5hp8pvLMZZXdXzJIz22o+EaJ+AW+tBCDGUNDZ2rtKnlOLpp5/u8vz244S7ujaZZCryMOPxe8hz9F4Hor0xmWMwK3OnfuAsp5Xa5hARKcgjRL9IAh5GojoaX46+tzoQ7dnMNkald16mPstpQWuoa5YXcUL0hyTgYaQ2WEs4GibbnkOtP5RwCxjosiaE1IMQXelpoctDzf7+rpKAh5HWWXB25QZ6n4bcXmlWKbsaduEP+eP72iqiSQIWBrvdTk1NzbBIwlprampqsNvt/b6HvIQbRlrHAKfoDKCubwk4sxSNZnvddqbkGMvUZzmlBSw6Ki4upry8HI/HM9ihDAi73U5xcXG/r5cEPIy0zoKLhtOBOnIT7AOGttUxNvk2xROw1IMQ+7JYLJSUlAx2GAcN6YIYRlpbwKGgMQuuLy3g4tRi7GZ7h37geAtYuiCE6BdJwMOIx+8h05ZJrd/on+tLAjabzIzJHNNhJITDasZhMUsLWIh+kgQ8jLQfguawmHHZ+tYD1dVICLfUgxCi3yQBDyNtkzASHwPc3vis8XiaPdQGauP7slwWaqULQoh+kQQ8jHj87etAJN790Ko0s/My9VlOq/QBC9FPkoCHiUg0QnWg2lgNuSHxacjtdVWcPUtqAgvRb5KAhwlvwEtUR8lz5vW7BZzryCXDltGhBWz0AUsCFqI/JAEPE1XNxhC0bHsOXn8LuQmWomyvdZn6rbVtC5NkOa3UB8KEItEerhRCdEUS8DDROg3ZRhZa974WXHeKUovY27Q3vu12Gasj1/plJIQQfSUJeJhonYRBNB3o2xjg9vKd+Xj8HiLRCNBWkEfqQQjRd5KAhwlPsweFItzS91lw7Y1wjSCiI9QEagBwSz0IIfpNEvAw4fF7cNvd+JqMlmuiyxHtK9+ZD0BlUyXQrgUsCViIPpMEPExU+aviIyCg/33A+a5YAvYbCbi1II+MBRai7yQBDxPtF+O0pphI6+M05FbxFnAsAWc6jZdw0gIWou8kAQ8TVf4qYxJGYwu5qTaUUv26zwufegnsuiLeBWFLMZNqS8EnoyCE6DNJwMNAKBrCG/C2TcLoZ/dDNKqxmk3YLSnsbaqI7890WqQFLEQ/SEH2YaCm2RixkOvMxdMQpDjL0a/7/PzldXibWpgzYzVV/raJF26X1IMQoj+kBTwMtI4B7uty9Psane2iNC+NEa4R8T5gkHoQQvSXtICHgdZZcNn2HLxNZf1OwD84fix1/hAn/XkOTa5KojqKSZlwu6xsq248kCELMSxIC3gYaK0DYVWZRHX/xgBvr25Ca02aPYWcNDNR1YQ34AVaW8DyEk6IvpIEPAx4/B7Mykwo5AT6Pga4MRhmwQPvcc/SrzCZFD89Ix1L2hftxgJbaAyGCYYjBzx2IQ5lkoCHgSp/FTmOHHyNrbPg+paAU0yK35wzjTNmFALtJmPsMxtOCvII0TeSgIcBT7On4yy4PiZgu8XMubOLmTjCKOSzdE2Uxs0/oyKWgKUehBD9Iwl4GGibhGEk4Nw+JGBfUwv//KycpmA4vm/KiBysaV+yp8HoW5Z6EEL0jyTgYaB1GrKnMYjVbCLdkfjgl2UbK/nx85+zo8Yf3/eNKQWUjFtBTdCYjJHlbC1JKV0QQvSFDEM7xAUjQeqCdeQ58/hqZwvZqdY+TUM+7/BiJuSnMakgrcP+fGd+vAsiK1aUXSZjCNE30gI+xLWOAW7tguhr/6/JpJgxMrND0g5Honz80dms3zQKaNcCli4IIfpEEvAhztNsJOC2xTgTHwP80qpy7l+2iWhUd9ifYjYxpaSBoOUrtNZYzCbS7CnyEk6IPpIEfIhrnYbcWoqyLy3gNbtqeW+zB5Opc5fFt44yoZ3rqQvWAUY9CFmWSIi+kT7gQ1xrF0SOPZeaxi1dT8IIBYyfFnuH3b8+eyot4a5XO8535aOjVir9lWTaM8lyyvL0QvSVtIAPcVXNVVhMFog4CUd15xZwNAJPngYLL+iwW2uj28Ga0vUfkU++dNL41Z1s9xkrJEsLWIi+kwR8iPP4jUkYNbHWaac+4DXPwO5VsG051GwFjOR75l/+w8PvbO32vvPHFmDNfYPKWBeH1IMQou8kAR/iPH4PuQ5jDDDsMwkj2ABv/xrypoAywZqnAQiEokwtSqcw097VLQE4dsxhOHPfpy5kjAV2uyzSBSFEHyU1ASulFiilvlJKbVFK/ayL4z9WSn2hlFqrlHpLKXVYu2OXK6U2xz6XJzPOQ1lVc1XsBVysBdy+D/iDB6CxEs74E4w7yWgNR8I4rGZ++63pnDWzqNv7mk1msm0F7Kw1ir1nuaw0hyIEQlKQR4hEJS0BK6XMwIPAqcBk4EKl1OR9TlsNzNFaTwdeBP4Qu9YN3AEcBRwJ3KGUykpWrIey1i4IT8M+dSBqd8FHf4Gp58HII2DWpdCwl+iWt9jl9fdwxzY711/JR+tzgfaz4aQVLESiktkCPhLYorXeprVuAZ4Fzmp/gtZ6uda69W/7x0Bx7PspwDKttVdr7QOWAQuSGOshyR/y0xhqjE/CMJsUmQ5j1hpv/tL4eVLs5/gF4Mzm0/eXcOwflvPOV1W93n/6hO1Y09cBbQlYuiGESFwyE3ARsKvddnlsX3e+D7zel2uVUlcppVYqpVZ6PJ79DPfQ02ESRkOQbJfVGNO7awWsfxHmXgeZI42TU6ww/QLG7H6Fn319JEeWuHu9/1ETojRZVqO1xh0vyCMv4oRIVDITcFcFB3QX+1BKXQLMAf7Yl2u11o9qredorefk5ub2O9BDVZeTMLSGJbdAaj7M/1HHC2ZdQp6u5urU93Faex8inm3Lp7EpFV+gHrfUgxCiz5KZgMuBke22i4E9+56klDoJ+DlwptY62JdrRc9aJ2HEF+NMs8H6l6B8BXztF2BL7XD+hkgRH2SdQ2TVU0ai7sXOvbn4t/2Ez/eUSz0IIfohmQl4BVCqlCpRSlmBC4DF7U9QSs0CHsFIvu07HZcA31BKZcVevn0jtk/0QWsXRGsLuMCpYdkdMGIazLyo0/lP/qeMa6vPJeL5CvZ81uv9jxmbh73wWVrwkuGwoJT0AQvRF70m4Fgf67V9HYWgtQ4D12Ekzo3A81rrDUqpO5VSZ8ZO+yOQCryglFqjlFocu9YL/Bojia8A7oztE31Q5a/CkeLAleKiprGF0xpfgvpyOOW3YDJ3Ov/XZ0/l6SsOx2qxwmf/1+v9p40owpKxhsZIFSlmExkOi4yCEKIPEqkFcQFwBbBCKbUSeAJYqnXv/4+qtX4NeG2ffbe3+35SD9c+DjyeQHyiG62TMBoCETIiNczb+78w8XQoObbL8+0WM1PHFMHks4yuilN+A1Znt/fPceagW7L5oqKGc8cbSxNJC1iIxPXaAtZab9Fa/xwYDzyDkRR3KqV+FRuvK4ao1kkYnsYgN6W8gDkagpPv7PLc+5Z+xStrdhsbsy6BYD1sfLXH+1tMFlp2f5clq4wknSX1IITok4T6gJVS04F7MboMXgLOA+qBt5MXmthfHr+HPEce/h2f8W3zu+ydcDlkj+10XiSqWbaxirXlRmlJRs+HrBJY3Xs3xLhxaxh12FqAWEU0GYYmRKJ67YJQSq0CaoG/Az9rN1LhE6XUMckMTvSf1jq+FlzhJ7+mFhdNc3/U5blmk+K16+fTEomVnlQKZl0Mb98F3m3gHtPtc8YXQVndJgCynBY27Kk74L+LEIeqRFrA39Zaf11r/Uxr8lVKlQBorb+V1OhEvzWGGmkON5PXUE1O9afcHz4Pd3Zet+crpbCltHsxN+OiWIGeZ3p8TmZKATsrXTQEQrhdRh9wAq8HhBAkloBfTHCfGELia8F99QY1jhKejX49PlZ3Xxc/+hFPvL9P6cmMIhj7dSMBR7svsBMKFOHdfhErd1SS5bISDEdploI8QiSk2wSslJqolDoXyFBKfavd57tA93UKxZBQ1RybBVe3l1fyryHD5cTcxdJC4XAE1q+l/onHiQaDHQ/OugTqd8PW5d0+5/CRWTgOe4TczGbcUg9CiD7pqQ94AnA6kAmc0W5/A/D/khmU2H8e33YA8gqP5EM9i5zU5i7P87+5jNuWPgBA1R+djLjt520HJ5wKDrfxMq606xGDo7NGkOLcTn3EQ5bLeMHnawpRLLXrhOhVtwlYa/0K8IpSaq7W+qMBjEkcAJVfGL1EuSf9Gs8iP7ldrAWnw2E8f/oz1nFjcc2dh+///g/XvLmkfe1rxgkpNphxAXz6GDTVgCu70z3ynflEmkfy/uYqTho9EZB6EEIkqqcuiP+Ofb1IKfXnfT8DFJ/oJ0/NJlIx4SycRXVD16sh1738Ms+Zivnh3GvI+vGPsU2exN5bbiVUUdF20qxLIBqCdS90+Zw8Zx5Bz8ks+kRLPQgh+qinl3AbYz9XAqu6+IihKhLCE2kmN8WF1jpWCa3jC7hoMIjnLw+Sk++mdEwBNoedonvvJRoKseemn6IjsRdp+VOgcJbRDdHF6Aar2Ur+qHc5ds7GeElK6QMWIjHdJmCt9auxn/9o/QD/ByyKfRdDVd0uqswmcm1ZNAbDBMPRTi1g38KFhCsquPiqs/nThbMAsJWUMOL2X+BfuZLqh//advKsS6FyPexd0+Xjit1WmqLlpNstmJSsiiFEohIpxvOMUipdKeUCvgC+Ukr9NPmhiX7zlVFjNpHjym9bC65dAo40NlHzyKO45s3FedRRHS7NPPtsMs46k+qHHsK/YoWxc+q5kGKH1U91+bhMSzEbtqdTUR8gU+pBCJGwRMYBT9Za1wNnYxTWGQVcmtSoxP7xbsdrNuNOLaI6thpy+8U4vU8+ScTnI+Wa65nxq6W8vHp3h8vzf3E7lpHF7P7pfxP2+cCRCZPOhLUvQKjzaIpUVciOrceycoePLKeFWr9MRxYiEYkkYItSyoKRgF/RWofoZmULMTQEvFvxm0y400dSHV+M0+ifDft8eJ94grSTT8YyYRJnzixkVHbHimfmVBdF995HuKaGvbf9wpjZdvilEKyDjf/q9LwJIzJxjf0DJ07MiM+GE0L0LpEE/AhQBriA92JLx9cnMyixf3w+Y1ab25EdbwHnxrogah55lGhzM7k33kBumo27zp7G4aM6D9p1TJ1C3k9+TONbb+F75hk4bD5kHtZlgZ7C1HxMVi/eYDVZTqmIJkSiEilH+WetdZHW+jRt2AGcOACxiX7y1u0AwG1342lsQSlwu6yE9u7F98wzZJx9NraxYwm1Ft/phvvyy3EdfxxVv/8DgU2bjCFp298F344O5+U78wk3jueZT8qkBSxEHyTyEs6mlLpIKXWrUup2pdRUrlqTAAAgAElEQVTtwK0DEJvoD62paaoEwO1wU90YJMtpJcVsovqhh0Brcq+7FoAbn13D2Q/+p9tbKaUo/O1vMWdksPvHPyE6/mxAwZqnO5yX78onVD+D5z9uiNcEloI8QvQukS6IV4CzgDDQ1O4jhqJmH96o0e3gtrnxNBhjgIPbt1P7z0VkXngBlsJCAL42MY+zZhb2eLsUt5vCP/6Blu3bqfjLP2DMCZ0mZeQ587Dnv8oPTt+L22klFNE0BsPJ+O2EOKQksiRRsdZ6QdIjEQeGdzs+s/HvqtEC3klOqg3Pn/+MstnI+cEP4qeeO7s4oVu6jj6a7B9cRc1fH8F1zWlkeJdDzdZ4cXdHioNMp43qYCUTXa2z4UKk2S0H+JcT4tCSSAv4Q6XUtKRHIg4M33a8JjM2kwVnipPqxiATG/bQ8PobZH/3clKyjXoOLeEoTX1opeZedx2OWbOoeHI5LY1m2Lysw/FsyyjeW5tKQ8Do/5V6EEL0LpEEPB9YpZT6Sim1Vim1Tim1NtmBiX7yleE1m3Dbs1FKUd3Qwvx3X8SckYH7iivip60s8zLljiV8vK2Gf237F+UN5T3eVqWkUHTPHyHFQsXaQti8tMPxXEcum7aNo6re6P6QehBC9C6RBHwqUAp8A6Ms5el0LE8phhLfdmqsDtyObJqCYcbu3UTBV6vJvuoqzGlp8dOKs5z89JQJZKcHueX9W7h35b1cueRK/CF/t7e2FBWR9Z3zadodJbr5A2hpexUwMjObUTPu4ztHjAKkHoQQiUhkGNoOYCTwtdh3fyLXiUHi24HXaiPLnoWnIcB3v3idUFYOWRdf1OG0UdlOrj1xHGOzR/Dq2a/yzTHfpNJfSaW/ssfbu+bNg4imqQLY/n58/whXPr4WD067MfpBxgIL0btEhqHdAdwM3BLbZQG6LgogBp+vDJ/JhNvuxvfWcqZ4y/BfeDkme8dFTHbXNhMIRTApE6MzRnPSYSex6KxFlGSU9Hh7x+GHo2w2mjypHbohjKFoU3jgzS8wm5QkYCESkEhL9hzgTGJDz7TWe4C0Hq8QgyMcRNeV49Vhsq1uTH//K7tdOaSedU6nUy/92yf86Lk1LNq8iI/2GPX2U0wphKNhnvriKRpbGrt8hMlmwzlnDv7qDONFXGy8b74zn2igmDe/qCHTYZHl6YVIQCIJuEUbo+o1QKwqmhiKanfRpKCFKKUrK7Du3Mb/TTqFnMzO/8l+8o0JXHr0YTy45kFe2/5afP9m32b+uPKPHfbtyzVvLkFPgFBFOXi+AowWsDV3CXdcECU71Sov4YRIQCLjgJ9XSj0CZCql/h/wPeCx5IYl+sVXhtdsLC1f/O/VNBSN5r2iGWSndl4N+ZvTCwB4reQ1mkJtL9MmZU/ihTNeYHzW+G4f45o3D4CmChuZm5dC3kTynfkoBZX+SrKcbhmGJkQCEnkJdw/GMvQvYSzUebvW+n+SHZjoB992vGYTlpDGtn0vZRMOJ8Nlw2Lu+J/Z0xBkS1UjkajGaraSZe9YjKc1+Xr8HupbOtddsk2YgNntpqkuP94P7LK4cJncPPe+hXAkKi1gIRKQ0GgGrfUyrfVPtdY3aa2X9X6FGBS+MrwWO6OrQEWjlLlHdbkW3CtrdnPSfe/ywpev8OT6J7u8VVOoifNePY/7Vt7X6ZgymXDNnUvTHhN6x0cQMJL0iLRsdlY6UEpewgmRiJ4W5WxQStV39xnIIEWCfGV40/IYs9d4MfZlelGnteAATp6czwPfmcnamk95o+yNLm/lsri44fAbuHzK5V0fnzePSH2AoA/YthyAEa48Zs55iaPHZOPzh4hGpSCPED3paVn6NACl1J1ABcZ6cAq4GBkFMTT5yvCmZjC2ohlzdjbbcDK1ixbwYdkuDst2cTa/IRTpfrTCt0q/Ff8e1VFMqu3fa9e8uQA0VWdg37wUJp9FviufTb5NZBZYiEQ1DYEwGU6pByFEdxLpgjhFa/2Q1rpBa12vtX4YODfZgYk+0tpYisjmYlyFwjF1KtVNoS67IP6zpZqq+gAAFnPvCfL3n/6en3/w8w77LAUFWEtKaKrLiw9Hy3fms7eygBdX7QKkHoQQvUkkAUeUUhcrpcxKKZNS6mIgkuzARB81VUOoifqomcLqKCmTJtEYDJOb1jEBN7dEuPhvn/DH917m5x/8nLpgXa+3zrBlkGnLJKo7FnB3zZuHf6efaF0lVKw1xgJHrfhDsYI88iJOiB4lkoAvAs4HKmOfb8f2iaHEVwZAyt5mTBoCYyYAdOoDTjErXrh6LuMKI6yqXIXL0vuw7qtnXM3NR97coQsCwHXMPHQwRHO1FTYvNcYCZ67khwuMpF8rLWAhepTIMLQyrfVZWuscrXWu1vpsrXXZAMQm+sK3HYC0cmNMb+3IMQCduiAsZhNHjHbzg8Mv5o1z3yDFlMhQcMMW3xYeXPNgfNt55JFgNtPUdBhsXka+Mx+AID5AWsBC9EaK6hwqYi3gnPImmjPsVNsyjO19EvC68jr+s6W6X49Yvms5z335HFX+KgDMqak4ZsygqdIO5SvIV1a0Vjz+ltFDJUPRhOiZJOBDha+MaFoBxeVBGsfmx1dDztmnD/jJD8v40T/f4NLXLmVDzYY+PeKKqVfw8tkvk+fMi+9zzZtHYKeXSECTtvNTXBYHEYKYlZJ6EEL0oscEHHvpdv5ABSP2g3c7tc4iCms0ofGHUd1gJOBsV8c+4J+dOpFbTx+DSZlIs/RtNGGKKQW33Q0QfyHnmjcPtKapNge1ZRn5rnxmz1xBTprUgxCiNz0mYK11FLhugGIR+8NXhrcxHRNgmjSe6sYgafYU7BZzh9Ny02ycPXke/zj1H4xKH9XnxzSHmzlv8Xk8sf4JABzTpmJyuWhqGglb3iTfkRerB2GVYWhC9CKRLohlSqmblFIjlVLu1k8iN1dKLYgtZbRFKfWzLo4fp5T6TCkVVkqdt8+xiFJqTeyzOMHfZ3gKBaBhD401xqZz6lSqG1vI3af/tzEYZuGnO9ld29zvRzlSHEzNmUpRahEAymLBedRRNO0Igr+GfGVhy64Mymqa4q1wIUTXEnkF/r3Yz2vb7dPAmJ4uUkqZgQeBk4FyYIVSarHW+ot2p+0Evgvc1MUtmrXWMxOIT9TuACC8x09DGuQVjMbzgbfTC7jNlQ3csmg1xdPu5ydH3MA5pZ3rBCfil/N+2WHbNW8ejW+/TUuThfwmH42REJkOCzXSAhaiR4kMQyvp4tNj8o05Etiitd6mtW4BngXO2ufeZVrrtUC0qxuIBMVGQJh31rKtQJHtyKa6IdhpEsb04kz+fcMRHFc8n4LUgv16ZCgSwuP3AO3KUwbGk19Thsm1iWMmWmholpdwQvQkkSWJLEqp65VSL8Y+1ymlEpngXwTsarddHtuXKLtSaqVS6mOl1NndxHZV7JyVHo+nD7c+xPjKiLQo7BX1bB2hyLRlUlkfIC+9YwI2mxRTRhTym+Pu4uiCo/frkd/593f49ce/BsBaMpqUggKaajIYUWOMR7Za66ltDhGRgjxCdCuRLoiHMdaBeyi2fWls35W9XKe62NeXv42jtNZ7lFJjgLeVUuu01ls73EzrR4FHAebMmTN8/6Z7txNoMEY0VI1MJRCCppYIeWkd14F7Zc1uTKYIZ0zv+8u3fV059UrSbekAKKVwzZtLw9Il5E02xgAvWxdFaxP1zSGyXJ0rsgkhEkvAR2itZ7Tbflsp9XkC15VjrKbcqhjYk2hgsbXn0FpvU0q9A8wCtvZ40XDlKyPQnAs0U1+SEy+0k79PC/hv72+nMu0uVjQcyZ3H3LlfjzxtzGkdtl1z51H30j/JbsgGICejieq6DLz+FknAQnQj0WI8Y1s3Yi3SRIrxrABKlVIlSikrcAGQ0GgGpVSWUsoW+54DHAN80fNVw5ivjIDPRl2WFWtOLpX1xuiD/PSOLeB//tdcvjvjPOYVzTsgj93buJdP9n4CgGuu0aWhAmOxac3EsZuN0GQssBDdSqQF/FNguVJqG0a3wmHAFb1dpLUOK6WuA5YAZuBxrfWGWH3hlVrrxUqpI4BFQBZwhlLqV1rrKcAk4BGlVBTjH4nf7TN6QrTSGnxlNFcWs7PIitvupqrBaAHn7fMSzpJi5ppZvfUcJe7+VffzScUnLD9/OSnZ2dgmTcK/O0j+tDBB/xZgDjWNMhRNiO70moC11m8ppUox1oNTwJda64T+VmmtXwNe22ff7e2+r8Domtj3ug+BaYk8Y9hrrCTiDxCqbmLzDKeRgGMt4Lx2LeCvKhp4cfVGLp87geLMjAPy6B/M+AH/NfO/4lXSXPPm4v3H/1J8QpTte43kv6WqkVMOyNOEOPT0tCTR12I/vwV8ExgHjAW+GdsnhgJfGQGvMShlQ26AbHs2VQ0BbCkm0u1t/75+sbeOp7b8me8t+84Be/TYzLGUZJTEt13z5kE4zNTKNAKmMgAag1I6Woju9NQCPh54Gziji2Ma+GdSIhJ94yuj2We85No2QnGW3c3G+iD56XaUahuIcs6sYtw5V1Pf4jugj99Qs4H3y9/n6hlX45w9G2W1Mm5vKnVjKxhl8RJJaMi4EMNTT2vC3aGUMgGva62fH8CYRF94txst4MJ8mhw1ZNmzqGoIdBoBAXD8yPkH/PFrqtbwt3V/47zx55HjyME5Zzb5ZV8SPlZxvH0VngbpSRKiO1KM52DnKyNQ66Cl1Bjx19oHvO8Y4Nte/Q//WPFJp2WF9tc5487hve+8R44jBwDn3Lk4yn1kNmqW+6fy/uZhPEFGiF4ktRiPSL7wnq2EGqBhjLEahdvh7jQLLhrVLN35L+754koaQ40H9PlOixOnxRnfbp2WPK1Mc7z9LdKtXc3HEUJAEovxiIER2LITMFM9OgMCYFMZnWbBmUyKf152Leurjyfdmn7AY1jnWcdDnz/Eb+b/hsxJk1CZmUzfXkf63PdoCBwFnHTAnynEoSCZxXhEsrX4Cewx1oDbXWgjRaXQHDBGROzbB1ycVsyCkgVJCcOkTOys38mexj0ok4nUuXOZXqbZm2JlnH8Ne+v6X/5SiENZIsV4nEqp25RSj8a2S5VSpyc/NNGr2h0EvBYsI9xUmZvItGdS3WhUIGvfAn7hs6+4bOHj1Ph7X4K+PyZnT+Zf5/yLKTlTAEg9Zh5ZjeD1F3J/5Fwee29bUp4rxMEukT7gJ4AWoHX+ajlwV9IiEonzbqfZa8ExaTzegBe33ej/hY4t4NVVq1ndcj87GpJTSkMphVIKrTVRHY33A+ftDXNXyuOUpMpYYCG6kkgCHqu1/gMQAtBaN9N1pTMxwMI7NhL2p2CfMTuegLuaBXf7SWfzv6f+L5OzJyctlu112zn1n6fywe4PsBQWUpfnonh7kItT3ia3dnXSnivEwSyRBNyilHIQKyUZK8wjE/yHgMD6dQDYZx2Bt9kbrwOx7yw4e4qdWXmzsKfYu7vVfitKLWKieyLOFGNERO30wxhTFqQxbKV8x2ZqZXUMITpJJAH/EngDGKmUehp4C7g5mUGJxAQ2G8XP7VOmtOuC6DgLzh8Mc/Lf7+KZ1Z8mNRar2coDJz7AnBFzAAjPmYI9BG80zuGuirm8u0nGAwuxr0RGQSwFvoWxdttCYI7WenmS4xIJaN5RgzXHRsiegj/sj7eA21dB2+wtpyLlOTbXJVLCef/5Q37qgnU4jjyCiIKM6mb+x/JnjhubOSDPF+Jg0us4YKXUW1rrrwP/7mKfGCzRKIHKFpyTC/EFjPoOrX3AkwraxvrOKBjNO+e/g8WcyCpS+8cf8nPC8ydwyaRLOHHkiWwphLy9QU4et5oW30ZIm5P0GIQ4mPRUDc0em/GWEyuQ3joLbjRQOFABiq6FyzYQ9puxjx+LN+AFYgm4IdhpLbhsR3ZSJmDsy2lxcuPhN3LiyBPJd+WzcaQivaKeinAmT7+zlmBYRkMI0V5PLeAfADdiJNtVtI18qMdYbl4MouYVHwDgmD6DmkCN8d2cQWNwb4cxwJe9+EciQTdPX/z9AYnrokkXARCJRthSbMb8cZilDbP51Rf5zNxTz6xRWQMShxAHg56qof0J+JNS6oda6/8ZwJhEAgJr1wAa++z5eANbAIiEjBEIrWOAtdZsDrxOjnlGd7dJiq21W2kON+MbmwvsZURlDa9M/T3Tit8d0DiEGOoSqQVRoZRK01o3KKVuAw4H7tJaf5bk2EQPApu2Yk0PYyoYj+/LFQCEWlxA2yw4pRQfXvwWzeGBnQr8k3d+gtvhxplfSJ3bi6lGMSO8Bhp2Q0anBVCEGLYSGYb2i1jynQ+cAvwDY1l6MYgCZVU48q2QYsUb8GIz26htMnqJ2s+CU0p1qFY2EO485k5+O/+3FLoK2VJsJsPr5/PoGB5+Y+WAxiHEUJfQqsixn98EHtZavwLIOuODKFRZRbghhP0wYwn41jHA1Y3GZIfWWXB//99TufJ/TuOznbFVMHZ8CL6ypMc3PXc6+a58xrvHszYvQGpTA+8GpnPvGjN1/lDSny/EwSKRBLxbKfUIcD7wWmy5+ESuE0kS2LAeAPv40QDUBGridSDaz4LbQoBqZxO5qbEW8XOXwvv3td3oqXNh5RNJifHDPR/i8XvYVGS0ymdUb2btqAfIcCZ/OJwQB4tE+oDPBxYA92ita5VSBRhL1YtBEvh8NSiNffJUALzNXrId2VR6gsxOrUGteQZmXsRvL1uO1rptbbiLngNrqvE9EoZoBFpXyGiuhUVXw9dugxFT9zvGpWVL+WD3B3jyIZxiptljw1G1GlqawOra7/sLcShIZCacH9gKnKKUug7Ii82OE4Ok+fNV2NLDmEaUAuAL+uKz4C5WS2HJreD3dky+AMVzIG+i8d2cApe9DEfEhqdVb4Y9q9sSstb7FeOPZv+I1899nYKMYiqLUzHXhHklfDR/WvzRft1XiENJIvWAbwCeBvJin6eUUj9MdmCia1prAhs3Y3eHIGs0WmujBWzPpqo+yOtFP4TvL+MNzwqOfvxSbnrp48RuPPIIuHEdFEw3tpfdDot/CNH+rSGXYcvAYrIwyT2JTUWavFofKyITWLKpDr2fyV2IQ0UifbnfB47SWt+utb4dOBr4f8kNS3QnXFlJpK4Ru7sFskbTFGqiJdrC6FCIpgYfuRlOyB1Pc6gZl6OJI0b1YdJiSrt3q2YLmFLAFPsjEgr0OdZ3dr3DJt8mPh3RiDUa4arg27x22LMdW+VCDGOJ9AEr2kZCEPsuf4MGSWC98QLOMcIGTjfe+p2YtOakj56kSKfzWdozAJxTeg7nlJ7T/wd9/fa2bgjvdvjb1+GcR6D05IRvMcI1goiOsLnQ+ONSUZPDYbs+NVrVJnmPK0SiK2J8opT6pVLql8DHwN+TGpXoVvP69aDANsZYht4b8BJVipVzruWPofPjY4BbwlGi0f38X/3WlqpSUHIc5MdeztXuAr+318snuifyzGnPUJ+qqE1zUldh4u6G0/jLvz/Zv7iEOEQk8hLuPuAKwAv4gCu01g8kOzDRtcD6DdjcClNuCQC1jXsB8GQczko9kbw0O196v+SUF85g8t2PUlXf966DTrJGw7efhPQCY/vNO+DhecYoil64HW5GOEew7AgLjqp6KrQbT+Xu/Y9JiENAt10QSik7cDUwDlgHPKS1Dg9UYKIzrTWBDRtIdfshqwT8Xo548RrOdkIg6AR85KfbCOkoRWnFzJs5kZxUW6/37bP5P4JJZ4LJbGy/ci2M/TpM/VaXp6dZ03jhqAqOXBfmnpSnseXUAucd+LiEOMj01Af8D4x14N4HTgUmYVRHE4MkvGcPEZ8Px5ig0SpVij3Zh7ExspeiZmP2W166nQzHZJ46/ZHkBTJimvEBCDbA3rWQO8nYjkagYh0Uzoyf/vVRX8e/eROHecATGEfxLumCEAJ67oKYrLW+RGv9CEZz5bgBikl0o3n9BoD4EDQcWbw0bQHlqVnUNEbjs+CiOkpVfWBghnvZ0uAH78HR/2Vsb3kTHj3e+BkzPXc6VVkQMpnYU53G+Xsu4PG31yc/NiGGuJ4ScHzSvnQ9DA2B9evBbMKWGYKPHoSGirbVkBuMteDCOsyxzx7L/L/+iofeSc4y9J0o1dYdMepoOP1+KDne2K7ewuSMsYRSFIuPcPKLGWVkmbykNu0cmNiEGMJ66oKYoZSqj31XgCO2rQCttU7+Eguig8CGDdgLMzCZ98DezyEcbFuMs9pYCy4YDvLNkjOJ5k3k+PG5Ax+kPQPmfM/4XrsLHp5H9rE/Id+Zj98VIMcb4E7738h3WIHTBj4+IYaQngqymwcyENEzYwbcRlLHpEDWKLjmY7A48Aa8FKcWsz62FlyqNZVbjx4ii1ZnjoQFv4GJZzB5xW4aU1dx23NhMi6ZiN75KTqqMZlkSLkYvmQ0/EEiUlNDSrSK1Owao//X4gCMQjytXRC5aTbqW+opq27E19QyuAG3OuJKSMtncvZk1hfWAbAiWMLsyhm89vmuQQ5OiMGVyEw4MQQENm4k//A6bK4wtZsnkKE1Gk1tsJZUSyaNwTD56XauefMatlVGyWm6hsXXzR/ssA0tTZz32SJqTA347Nm8adlFuGA11uZjgVGDHZ0Qg0YS8EHAv3Ile/77ZiK12eROrcezfiP2i74kUFJAREewYHTH56fbOK/4PHZlh5iWOX6Qo27H4iTNmooJxZd5GXzz/SiXnrCHcRPKBjsyIQaVJOCDgCktDVLMmNLTyZ68F+cFN2OfNIk9tdvIrtcQMWr85qXZmT/ubGPqzFCiFLaLXmTp81/n6EITcz/ezWjTCOrLVpEyO4DTau/9HkIcgqQPeAiLNDYBYJ8wAffkMLkzQygFznknAODb+Dl/+muEguXGmFqrrZFNngpW7/TREu5fGcmkUYopOZOxjarGVRBgScvXmbPLzdmvnEsoIssUieFJEvAQFaqoYNvpp+N79jl0NIqNHTizY6MCs0YDUJ1t4dWjFOWTZgPwxq6n+c7rp3POQx9Q1zz0ktoU92TOt+0mZ0ojJU0BzglvZnbmZIKR4GCHJsSgSGoCVkotUEp9pZTaopT6WRfHj1NKfaaUCiulztvn2OVKqc2xz+XJjHMoSnG7cc0/BsesmYR276b83QyaTTPAkWWMtQW80QaeO96M15yLzaw49+mt/L7yRP522ZHkpA69dVMn50zhxvwcPlg7EVe5l9+3LOe3I44itXWZJCGGmaQlYKWUGXgQo47EZOBCpdTkfU7bCXwXeGafa93AHcBRwJHAHUqprGTFOpRE6uqI+v0oq5XCu+7CPmECwc2bAbC7mowiPDG+gLHacV2jlWKXmXQczMmbykmT84dk0fNJ2ZPYa0lhnTuN5q+2ETSlsXHjOrbVbuPFTS8OdnhCDLhktoCPBLZorbdprVuAZ4Gz2p+gtS7TWq8F9u2wPAVYprX2aq19wDKMhUEPaToaZdfV/8Wua6/tUMfBvOp/yJtVhzVlb7z7AYxawJm2TDwNIdKzFDW3X8nyqceyt64Z/4oVVD/8cPxc/2ef0bB8eXy7pbyclh074tthn4/Qnj3x7fply6h/Y0l8u+Kuu6m67/62WCO9l6LcV54zj9SULDYVQuGMcvxRM2d9Nosn1j/FA589QGNLY5/vKcTBLJkJuAhoP9K+PLbvgF2rlLpKKbVSKbXS4/H0O9ChQplMuC+/HPcll3RowUa9FVgyHJj9uyC7bYhD+zoQKambuXTJZfx08du8+5UH71NPU7vo5fi5vqeeour3f4hvV/3+95Rfd118e89NP6X8xh+1nf/0M3j/8Y+2GALN6KBRWzhcXc22075J47vv9vl3HJs+kS2HNeHb4qKOadxneZhrJl3B4rMXS1eEGHaSOQytq/8HTrQ8V0LXaq0fBR4FmDNnzkG70qMOhWjZsQPbuHGkLzil0/Gq9TlYcktJ1690SMA1zTW47W621gc5auzhXDXzfgpPnk1Bhgv3nx7o0IrOv/VWon5/fDv7+98n0tQU33Z/7woIt9VcKrr/Pkz2tuFhhXfd1RZvMEhKQQGWokT/PW0zM38qa6o/oq4ph5yaXE7P/ifUbYXsEwCI6igmJe+GxfCQzD/p5cDIdtvFwJ5uzj2Q1x50PH/6E9vP/w6hyspOx3QwSHB7GbaCNGNHTscWcLrVmAU3KjOPk0efxJTCLNwu4wVc+1Z0Sk4O1lFts84cM2eSeswx8e3UY44h9fjj287PysLkcHQZr6WoiMOefALbOCMWz0MP0fDWWwn9roePmIoywe7iQpq3VlCj0/AuuploJMTVb17NfSvvS+g+QhwKkpmAVwClSqkSpZQVuABYnOC1S4BvKKWyYi/fvhHbd0jKuuwy8m/5GZb8/E7HIs99n+Kjq7C5YzvatYB9QR92UwYQpSLyH55ZtZb3Nw9sV0y0pYXGt96m6T8fJnT+lOwpAHyR56BlWxnL1LHsqAuDv5aS9BIKUguSGa4QQ0rSuiC01mGl1HUYidMMPK613qCUuhNYqbVerJQ6AlgEZAFnKKV+pbWeorX2KqV+jZHEAe7UWve+CuRBypKXR9a3v93lsVBLOoHaFNJcjRDNiw9BC0VD1AXrMOs0lMXLP8vvIf3Li5icdjLHlg5cGUqT1crohc+go8Z71OD27YT37sU1b16X5+c587CQzue5LZwGTLZlkBK2oJxubj5yiFRxE2KAJHUqstb6NeC1ffbd3u77Cozuha6ufRx4PJnxDbaW8nL2/vw28n50I46ZM7s8p6FhPDVfusk5paJD67c2UAuAiqaiQ24eOOYpJuUXYtZpAxJ7e8pqjXfaVz/8ME3vvc/YN9/EnOrqfK5S5FnH8kXBblCKgibITdkIez4Dz0b0rEtZV72O6bnTB/aXEGIQyNuOQRRtaCAaaIYUS9cnNNcS3PQVtpLRKN8WyB4bP+QNGP9DEG5xAibmFE6mMC2X/PTBratQcOedjHzssXjybT/UrdWYjBcmUZ4AACAASURBVIkE0jxQMobm8kY+jExm85KHYfEP+eeaR7j4tYvZWLNxoEMXYsBJAh5E9kmTKHnuORxTp3R9wqs3kOd6GVvJYeCvhpzS+KHWBBwIOnC4V/HqV+/xyLtb/z975xkmVZE14Lc6T855mEzOWZCczQFUxKy7Zl11Vz/dlRVd1+xizogRxSyCCsqQcxYYwiQmD5Nj5+76ftxmMjqEGYL3fZ5+pruqbvW5t3tO1z11Ahb7sfvnnkw0JhNeffsAUPPzUrKvvApnRXPr0YDw3gghORwXhSUtndflFTxaOBJuW8OkXlfz31H/JTEgsa3pVVTOKlQFfBrjSrmAijQ9xhhPEGALH2CAunoT+rCf+D5jCS8uO4hOe/pEwBlTkgm89BKEoXlY9Kg4xdySFq7HXVPDEyG7eN//XYjqR4AxgIuTL8akUzOkqZz9qOkoTyF5t92OqXcvwu69t81+u7Y7VVk+xF6uVRzzQlqvgKvrTaS4nubtq/uhv8gPvfb0+U01pqQQ8cgjrdp7hHZBunzYGmphPBBi98bHsh/qSmHTW7gtlSzqMQZfvS+T4id1vuAqKp3E6fPf+idEGxaKNiCg7c7yTOz7tgNg9LOA0LYKQ9YJHeW1GqL8/Qk2BeNnOoot+RRjS0+naM6chvBljUaDl4xnr18RGj8/zIfhfec0UjdsBqcV4bTy2f7PWJy1+BRLrqLSsagr4FNI0+iyVix7FJ8DGxBewehlEQTFg67xVr7CWkGgKZBS1xbKtPDIN3quGR5Hn5ijKPRTiC0jg9qffib42msbgjfCDMnkuZZg7DsYa0YB83tfxqh9JUy450mEELzhifJTUTmbUVfApytjH6L8cE+MKSmIysxm9l84kognCKdpL4XOtfy0p4jSutMzr67ftGkk/7KsQfkCJPt3B+Gmpls0toxMfoj5gqd9PgdP9F6ItQ5hrmgWTq2icrahKuBTRPl775F53vm47UepXhw9kOrfajF27Qrlmc3sv6AoYB9dINbCq7i7x0vsmD2ZsZ0YgHEsCCHQ+vsjpcSyZy8A/cIUz4+94Vpwu9ETDwXbcFtqob4MXhvCutR/ceG3FzbYu1VUzjZUBXyK0HeJw3vIEDSGNhKn527CuWcZrvJyjHHh4LQ08wEGpRy9USjmhpiAAIQQaDSnjwdEW1R99RWHZszAsns3A6IScTu9WeevlKq31AbwleMcpry8GqshCC56mcgB1xPtG021rfoUS66i0jGoNuBThP/UKfhPndJ25+rn0BQfBMAUaoBimvkAg7ICDtBaMUZ+zYKtvmSVxXP9iISOFfoECTj/fHC5MPXsSYLFhdsaQ5p3FoakJCyHqoiOqqWrroxaqxPTgFkkA+92aTukWUXlbEBdAZ8CpJTIJqkfW3H5u9QGzALA6O+x6zaxAVudVsxOM1aHC53fPspqIL/S0pEinxQ0Pj4EzZyJ0OkI8dGjdXShzJ6LoV8fLL/tZkRSMG+aXifMz6gcUJYB399Fbd1hCuoKTq3wKiodgKqATwGOnBwODBxEzbJlbQ/wDsaca0EbFITWWQh6H/BrzBJ2pBRRiKYfzuzZLLx1DP88v2dniH5SsOzcSc6VV5HkjETioqZrJK7KShz+Q6DsACUFh9iWUwmWCuTe77lyyUye2fzMqRZbReWko5ogTgHCaCTo+uuaeQU0sH+JJwfEQYxduyIqMhX7b5Pcvkc2pSxWLyL8vU7L+m+/h8bfH2m301NEkA1kxmrpClhqgzEAd3++k1IZSOrfxyL+vo+/F28iykdNU6ly9qGugE8B+qgoIh58EGNSUuvOnQuQm97Elp6ueECUpbdyQSu3lgOQbv8apz6TWz7YgsPVsqze6YsxKYnE778jpPcIpMuLbd5laLy9sWSXgW8E/w5bxYc3DVN+WIx+TIqfRG+fY6++oaJyuqMq4FOAq6qqIX9uK676BMfEN3GbzRhTEqEqt80NOACH2423wYDd5T6tQpDbgxCCxEAT47b6sPvwXkz9+mHeuQOSJ9Cn+Hvigprkgvj5EYreHcsb21/D4XacOqFVVE4yZ9Z/7VlCzo03kX+U/A8IgS2/DABjuA8gW62Aj9iAnXl3MKrLYD6+ZXhHitthJOfu5a7UYoJ3pGMaOhhb2j4cgYPBXE5d7k7+8eUuFu0qhOSJHOw+kXf2vMfesr2nWmwVlZOGagM+BQRfdx3awDZChrd9AFW52NKV0kTGQM9qr40oOKPWSJlVe8rz/54IsdMm8Zd1U8nutpyy3iloXoXaLCfBgHfuStJLhtM9wg/6T2JU8niWWsqI8GldtklF5UxFVcCngMDpl7fdUbwHDu/FltEHXXQUWku+0t6GAna5XeiD1zB/nYnB8UEMSzzz8iaE+RnJDBmChuXs15cxICmJ2lUbCR7RF01WKt/ecX9DcIlWaIgo2g36LIgfcYolV1E5OagmiE7GWVmJs7S07RwHF7wANy7Blp6OqWs3xQ/WNwJM/s2GlVpKMWq9kU4/ugR7N1RBPtMQQhAfGEe3AgPJNz+PqW8fzFu34gwfAXmb0DjNAKQfrsXtciIX38+/Vv8fL2176RRLrqJyclAVcCdTtXAh6aPHIM3mNvuly4UtKwtjt65QntFq9QtKPbgYr+44awbyzOV9SQn37WixO4zEEF8KAuPY390Hv8mTwe2mttgfXHbIWc/GrHImz13Nsv3liGu+xpA0Hr329Ey7qaJyrKgmiE7Gd/wEdKGhaHxaFKxc9RyUpGEf8C9wOBQXtP2vQ48LWs1RYa0gWKP4xYb7nbk2YICEUG+Wl8Tw3Hk5XD5uDPouXajdnk1QogkyUxk6ZRKPXtCTEUkh4B3JY2GPg5ohTeUsQV0BdzKm7t0InDGjdYdGB1oDtvR0AIxxEWAub5UFTUpJibmEXPMeBPDmqoxOkLrjiA/xwWmOxuF2kFm4B42vD/UbN+MKGwaZK9BqBH8ZnUSAt2fVW3oA3jyX7AOL1VSVKmc8qgLuRKTTiXn7dlx19a07Rz8Al7+DNT0dtFoMfp7imi1MEPWOeiQSk4zG16Sjd/Tpl4D9WEgM9cFlVYIs9pfuxVl8GFwuamvioHQf1BQCsCO3kjmL9iL9o0k16bl44yMsz11+KkVXUTlhVAXcidjz8siZdQ21LXNAuBsrGdsOpmNISEBT6ynnfpQgDJN9IN0j/Liof3SHytzRxId4Ix3BGDTe7LEfImnZUnSRkdTur1EGZK4A4EBxLYt/K6TArGXsjak8POxhRkarmdJUzmxUBdyJ6MPDiX3rTXzObaE4lv4T3hoNUjaGIJdnKHXgAuObDS2sU1aENfUmwo9kDTuDCfM14mPQE6BJIK08DZ2fH36TJ1O3aSfm+nDIUhTw9MGxrHloArFB3mg1Wq7pdhXee74FlxoZp3LmoirgTkTj44PfuHHoI1oEE0QPhK5TcFssOPLyMHZNgfJ0pQinrrmL2feZ3wNQWu7PyoNHcWc7gxBCEB/ig9aRyL6KfVTbqvGbMB6cToq3+CsrYLcSau1l0CKlpNrsgKyVZPz4N27+/nKK64tP9WmoqBwXqgLuRMxbt2I9eLB1R/+ZMHE2tsxMkBJjt26eMkStXdDCvJSyQw6nnmGJwWdcJrS2SAz1wVLVHZd0sa5gHd7DhqEJCEAXHQPmMji8u2Hs7Z9s468fb4WUiXhNn0cxLjVXsMoZi6qAO5Hi/zxJyYsvNm+sLweHFVDsvwCm5BRFAbew/wL4Gjw+v9LIRf3ObPvvEeJDvCkuCSPQGMSq/FUIrRb/86ZhTsvD5RC4fvu5Yex5faK4bGAMbgkxPS7hh0t/YHBwr1MovYrK8aMq4E4kZu5cwv/+9+aNqf+Bl/qA243t4EGEyYQ+QLRZB87qtJJTk4MRf5D6MzoPRFMSQn1wugWDw0awtmAtTrcT/ylTkFYr+Zu6kPv0Zw0VRC4dGMPVw+IaQ5QLd+Ce25vl29/C1WQzU0XlTEBVwJ2IMSkRU7duzRv7TIcJj4JGo2zAJScjqrKUvhY+wDtLd7IocxHSrqx8D9ec/mWI2kNCiBKUEu81hBp7DTtLduI9dCjagACkIYCAmHJwNJ6r2y1ZtKuQzdkVENaDtXH9uW/366pbmsoZh6qAOwlbRgY1P/2E22pt3pE4GgbfqIxJT2+0/0IrG3BSQBKJAYn4GhTvh2GJIR0tdqeQEOoNgI+rNzqNjtX5qxF6Pb4TJ2IrricwqQZRsKlhvMPt5tmf9vPZ5lww+jL6qm94efzLTIqfdKpOQUXluFAVcCdR+8svFNz/ALia3CbXFEHpQXC7G5L0NLigGXzBL7LZHOHe4eg0OrwMeow6DbFBXp18Fh2D4oqmpbBSMiRiCKvyVwHgN2UybrMVc6kv9T9/Qf799yNdLow6LZ/feg4vXtEfUDwpJoQPRrP8CVz1pafyVFRUjglVAXcSwTffTOKi75vngNj5Kbw+FGw1jSHIDWWImteBA1hbsJYycxmWuihC/QxnhQcENLqi5ZTXMzZ2LFnVWeTV5OEzciQaX19qKmJxpm/Dtm8/zpISALoEe6PRCKwOF3anG6oL2LX1bS5cNJ1D1YdO7QmpqLQTVQF3EhqjsbX9t+8MmPE+eAU2eEAoJojWWdDya/O549c7qLJVk5t1Dk7Xme3/25KEUG/SS+oYEzsGgFX5q9AYDPiOH09dlh0//wwSP30bfVRjcc7SWhsTX1zFp5tyIKIXMTcvJyo4BavLerS3UVE5rVAVcCcgnU7K3n6nYZXbQFCCsgmHYv/VBASgC/JT6sC12ICL8I7glfGv4JaSgJBMuoX7dZL0ncO5KaHkV1qoqvEnKSCJlfkrAcUM4aqzYik1oMlbh3Q4qPrqK6TLRaivgcm9IugZpeRLDg3vxftT36eHMMHRau6pqJxGqAq4E3AUFFA6dy6WPU3qmVkqYf+PYK0GwHbwIKauXRFVh2irDpxeqyfePx4hJC63jeQzOAdwW1zYLxqDTsPX2/IZGzuWbcXbqLPX4TtqFMLLi5riIMhMpW71aooenU3dqtUIIZhzcW/OSWqyGVm4E+trQ/hk+d9xqGHKKqc5qgLuBAzx8XTbugX/qVMaG3M3wudXw+E05JEcEN089l+A0EYF7JZuPtr7Efsq9uGyRlFvDjhrfICPEOClZ0qvCL7fVcjI6NE4pZP1hevReHnhO2YMtfkmZMZKfMeNI/7TT5RwZQ8Wu4s3VmZQUmOFyH5sG3otzxb+ytqCtafwjFRU/hhVAXcSWl9fNN7ejQ2JY+HmpRA9AGdxMe66ukb7L0BwYxBGVlUWz299nl2lu7CXTcRtiz4rEvG0ZPrgWKrMDsrLo/A3+DfzhnDVObDkVCFK0/AePBgAW2YmFR99xOEaC3N/OcjPe4tBo+HcaS/x5UVfMj5u/O+9nYrKKUdVwJ1A5cIvqPrm2+aNBm+IOwf0Xtg8+SEaXNB8I5vVgUsJSmH5FcuJ9o3GGP4TCNtZtwIGGJ0SSpifkW+3FzM6djRr8tfgcrvwHTsOYdBTk69UyThC5cKFlL37LrFaB6l/H8f1IxIa+nroA2HhtdQVbDkFZ6Ki0j46VAELIaYJIQ4IITKEEA+30W8UQiz09G8SQiR42hOEEBYhxE7P462OlLOjqVm8mNqlSxsbpISNbyk+wNDogpaSctQ6cOHe4ZgdZjSGcnCbiPA/+1bAOq2GywbGsPJACYNCR1Jpq2R32W60vj74nDuK2kI/ZEajAo54+GESFy5EFxREl2BvpMtFrdVj99Ua2Fiyk0mpt7O3fO9R3lFF5dTSYQpYCKEFXgfOA3oBVwshWmZNuQWolFKmAHOBZ5v0ZUopB3get3eUnJ1B/McfEfPKy40N1fnw8//BoTUAWA8eRBcZiTYgQFHATey/NpeNOevnsLd8L7kVtYjaoYD2jK8FdzSmD4rF6ZaUlSSiFdpGM8TUKTjrJNbtWxrCkoVGgz5aCcsun/8Be669ibH/+ZltOZXgHUzvW1YxKfE8Ao2Bp+x8VFR+j45cAQ8DMqSUWVJKO/A5cEmLMZcAH3qefwVMFGdLdEELNMYmK9bALvCPjCYuaBmK+cFc4akD16iAc2pyWJazjDJzGQcLNdTkT0evFfh7nZ31VLtH+tE3JoAlu6oYFDGoUQGPHw9aDbU5GshZ3+o4bUAA/hGhTOoXS4iPkkPZzzuEJ0c9SUzBb0rUoYrKaUZHKuAYIK/J63xPW5tjpJROoBo44lOUKITYIYRYJYQY3dYbCCFuFUJsFUJsLS09PUNQzVu3UvzEf3BWVjbv8A0Dr0Ck04k9M9NThr51DohuQd1Yc9Uazo05F5+gNCISFxHhbzxrouDaYvqgGPYW1tDDfzjplekU1hWiDQjAZ/hwavK9kU3swEcIvPwy4ua+yHNXDaSL3ok9z/PVM1dg/vZWHlg8S03Wo3La0ZEKuC0N0TJ862hjioA4KeVA4AFggRDCv9VAKd+RUg6RUg4JCws7YYE7AvuhQ1QvXtx8Bbz+VTjwk9Kfm4u02z0bcB4XtBZBGFqNFp1GR4WtBJ3eRqT/2ZED4mhcPCAGvVZQXqJ4gjSsgqdNw1Gnxbbp1zaPO/KjlP3QI+y5chZ2swW8g9HO+oKSgChKzCWdcwIqKu2kIxVwPtClyetYoPBoY4QQOiAAqJBS2qSU5QBSym1AJtAijvfMIHDGDLpt2tjogiYlbHobMhQl0hCCfMQDQqODIKUOXLWtmht+uoEtxVsoM1eQnp2CtTaJ8LNwA64pwT4GxncPZ/luiPOLa1TAkyaBENTuKoDaw0c9PvuyG3ix+0XsLlVCko1xI/jgvI+4OuVy2PGp8hmoqJwGdKQC3gJ0FUIkCiEMwExgUYsxi4AbPM9nAKlSSimECPNs4iGESAK6AlkdKGuH0sxcIATctxsmPwEoEXBoNBiTk5UgjKAE0OoBKDWXYnfZ0Wv0LNm/DXvpJOqrEs7aDbimzBgcS1mdjSSfoWwu2ozZYUYXHIx3/57U5Jkga+VRj500dRjPzr2bwfFB1K1bR9V336EVWti5gF0/38czqQ+c8bX0VM4OOkwBe2y6dwNLgX3AF1LKvUKIJ4QQF3uGzQNChBAZKKaGI65qY4DfhBC7UDbnbpdSVnSUrB2FdDrJve02aleubN4hBBiUrGi29HQMcXFoTKZWdeBSglL47MLPGBA+gGLndoJ7PIfVcvZFwbXFuO7hBPsYqCjpisPtYEPRBgD8LrgUe40e24bFRz1WCEFskHLHUfLBR1R99rnSMegGto26k5VV+yi3lnf4Oaio/BEd6gcspfxRStlNSpkspfyvp+3fUspFnudWKeUVUsoUKeUwKWWWp/1rKWVvKWV/KeUgKeUPHSlnR+Hy5Ph119c3Nm5+F1Y+0/DSevCAEgHndkNFcwXcdJW29fBWeoYmA7qzMgquJQadhov7R7NlfwA+Ol9W568GwG/KVABqV29sqKV3ND5cf4groi/B97/PIITAbbUy/pV9fBx4P6HVxZC76XePV1HpaNRIuA5EFxZG0jffEHDBBY2NRbsgV1nN2fPycOTk4jVoINTkg9PaoIAL6goY98U41uSvoaS+kp27++PrGArwp1gBg2KGsLsEsV4DWZW3Crd0o48Ix6tnEjWZErnjk989fnhSMJcOTUAfozjfOMvKEBoNPoGhuBfdzYff30lh3v7OOBUVlTZRFXBnc8lrcN13ANQuV9yi/CZObMwB4amE7HQ7GRE9gli/WH5J34GzPhmTOxHgrN+EO0LvaH+6R/hRWZpCubWctPI0AAKvvRlbtZ6a+XPhdzKe9Yj0Z/aFvQjwUmzqhrg44j//DO/Bgyk6/xmKfnNSfvFVrV0EVVQ6CVUBdyCHn3mWoscfb93h2ZSrW56KsWtXDF26QJlHAXtWwPH+8Twz+hkSAxIpdu4gsPvzdA9JACDiT7AJB4otd8bgWDJzY9GgYWXeSgACLrsMU7d4Dq934Nrw4e9PAuRVmJn17kbyK80NG6Ixsecw885Xib7vAXQ5P4PDQtV332H15OVQUekMVAXckWg0CNHkEu/5Gj6bBdYanJWVmLdtw3fiBKXvSB043wjc0k2ltXFVtu3wNvqF9aXKLDDqNGdtFFxbXDIwGq30JVTfrcEOLDQaIp96EZddS+lLL8EflKO3Od3kVpgpqm5uM04YOIaQKQM4/MPdrPp1NiXPPkfV5583HpedrXpLqHQoqgLuQCIeepDIf89ubLDVQW0RGP2oW7UK3G7F/ABKEEZICgjBwcqDjF04lhW5KyiormLTlnGEMoqSGivhZ3kUXEvC/UyM6RpKVXkK+yr2cbhe8f/16tOboGkjqdzrxLLotd+dIyXclxX/GMfQhODWnVH9eXHIJcyu3kzUoq8I+etfASWAJuu886n66itA2RBVlbHKyUZVwJ3J4Bvg1hUgBHXLU9GFh2Pq3Vvpa5IFLdAYyF0D7qJvWF/W5/4GGhv9w7tzuMb2pzE/NGX64FgqSxXb+JGgDICwx15E66WheO47SKfzd+fQa5Wv+lfb8nl9RUazvofHv8D8afPx8zair9wMgDY4mMjH/o3vGKVGXf3q1WSddz627OyTdl4qKqoC7iBqU1PJuuhi7PkFrfrcNht169bhO2E8QqNRsntV5TVswEX6RHJb/9sI9QqlyLETv4R5zOg3hIIqy5/GA6Ipk3pG4KOJxkRYgxkCQBsYSMQtl2M97KbqtcfaNdf6zDI2ZJbjdDXWjAs2BZMcmAwrn2bLj/diqylE6+9P0NVXo4+IAECYvDAkJjZ4VFT/8AMlL/4PabefxDNV+bOhKuAOQmMyoY+ORhfiue3N3wqvDYPCndRv2IA0mxvNDxXZHKkDZ3fZ2VmyE6dbWdFtKd5G75DeFFcpdswhCUGn5oROISa9lov7x1Bf1Y0NRRuxOC0Nff63z8E7WkPJ/G9xlv9xcMXTl/dl/k1D0Wlbf/X3D72OW8ICmJ/tCdhsYnLwGT6MLm++gcagZFqz7ttP/fr1oFc8LCy79zT391ZRaQeqAu4gfEaOpMvbb6HxOpI4R0BwIvhGULc8FY2PD97Dhytd5Y0eEDtLdnLdT9extmAtGaUVrFl7MUGuMSxLU2yfk3tFdP7JnAZMHxyLraYHdpeNzUWbG9qFVkfkfX/FbXdT8tjf/3Aeo06LXquhzubk0e92U15na+jrEd6fNya9wU19boL9S+CzmYrdvg0iHnqQhM8/QwiBdDjIu+MOCv/16ImfqMqfClUBdxCtNmxiB8OshUjfCGpXrsBn9OiG1VRjFrRkeob05IWxLzA0cij7yw+g99/F6MTuLNtbTJ8Y/4YQ2z8bA7sEEufVG400NrMDAxgvuIvgfhqqf92Eefv2ds13qKyeb7cXsCm7eYT7qJhRGLVGrHWHecSWzaHa/KPOITyrX3Q6Yl95hZC//gUAZ2UlWZdcSv3Gjcdwhip/RlQF3AFIh4OD54yg4qOPGxs9rlLW337DVVqG3xH3M1ByQPhFgdEPP4MfUxOm4qP3Ide6A6+o7xga3YcdeVVM6RXZyWdy+qD4BCdiq00hNXdl8x84rZ6we+5F5+Wi+F8P/eGGHECfmABWPzSe8/tGtdlfmDKODd7eHKzLBaetoXzU0WTzHjQQL8+Gqqu8HI2vL9ogxVxkz8mhNjUV6XYfdQ6VPyeqAu4A3DYbAZdcjDHFU9nYWgNPx8L2j6hdngpabcPuOqBkQQtJocZew3cZ31FlrcLtlqzN2UuP4B5syKhHSpjS+89pfjjCZQNjcNX3pNxayv6K5iHEmnNuIuJcLbbsAioXLGjXfCG+SkTh9txKvtya16wvKTCJJZcvYUrCFPj1cWrfmwB17csnbExJIeHTTzB17w5A1dffUPC3+3CWlbXreJU/D6oC7gC0vr5E/vOf+IwcqTS47DD4JgjvRW1qKt5Dhyr1347gcUHbUryF2etmk12Tze7CCtZvOJ9A53iW7T1MfIg33SP8Ts0JnSZEB3oxKHQESMHKvOZmCHRG/K65F59IK6UvzcVxuP3J199bk8WbqzKxOZsHdPjolYx1B3qfx7QuUaRW7DkuucPuvYf4BQvQh4cDSoRk9fffH9dcKmcXqgLuANxWa/NbZJ9QmPYUdmco9sxM/CY0MT+YK8BSASEpTOgyga8v/po+oX0os2dhDF/MuOQk1meWMaVXxJ8qAONozBzcC5c1lsWZv7Sys4vBNxA5Sou02Sh57rl2z/ns9H58fftIjDptm/3REQMYGz+ZPqF94HAapP5XyV7XToROh1ffPgBIux3Lzp3Yshr9id0229EOVTnLURVwB1D48CMcmj6jsaG+HKRUzA/Qwv7b6AEhhKBbUDf0Gj2ZdTswhKxFOOJwuCRTe/957b9Nmdo7ElE3mNy6g/yU/VPzToM3hml3E9KjhpolS9q9CeZn0hPkY8Dtlry3Jou8CnPzfoMfT41+inDvcOTe7/h1z8e464+vBqEwGIj/bAFhd90JgGXnTjLGjceyc+dxzadyZqMq4A7Ab9IkAmZMb2x4Zxx8fxe1qakYe/RocOYHGhRwsU8QL217icK6QuxON8v25ZDk3511ByyE+hoYGPfn8/9tC2+Djgvip+O2duG/G5+mwtoiT//QWwgZpEMfaKD4if8cU6DEwZJaXl6ezpLdR6+gvKbrKO4P0LOsfIeyCq4tPuZzEEIgPB4wwssbnxHnKCWpAFtmJs6KM672gMpxoirgDiDgwgsInjVLeSEljLoPV/xkLDt2NDc/gLIBp9Gxz1nLh3s/xOwws/lQKVt3jCRUjmXlgVIm9YxAq1HND0f4x9QeGCpnUmOv5elNzzTvNPqhGXUnkf2KsGdlUf7BH2dLO0KPSH+W3jeGW0cnAVBSY21l5hgdO5pXxr/C1PipsHUejteGNmayOw5M3bsR87//ofFR7M1F/36MnFnXqB4TfxJUBXySqd+8GVd1dWODEDD0FmqzJLjdjdnPjlCeAUGJjE+YVNTskgAAIABJREFUxNqr15IcmIzBuwiv2PmkhERQZ3Oq5ocWhPubeO7iqdjKxvPzoZ9Y1XJDbtit+CYa8e0ZRNkbb+AoaB0OfjSiA73QaAS1VgeXv7mef3+/t1m/EILxceMRQlCTMIKLusTwY80BpfMkJOuJeuJxIufMQWg0SLebmp+XIl2/n+1N5cxFVcAnEbfNRv4991L8xH8aGyuywVxBbWoquqgoTL16NT+oSRIeH70PQgj2VuxA53eAktIofAxaRiSHdOJZnBlM7hXB9KTrcVkjeHTtHGrttY2dXoEw7FYiUw6AdFP8nyePWYn5GnXcdG4i0wfHHnWMzTeCHjHnkBiQBOYKXO+MgezVRx3fHozJyfico0RI1q1eTcF991H76/ITmlPl9EVVwCcRjdFI/AfzCb3zjsbGnx9GzptC/bp1+I0f39yToXAnlKRxMKIb9624j+zqbGqtDr7aWkKMqSdr9tsY1yMck77t3fk/O/++sB+h1mupspXz9MYXmneecyf6QCPhk2KpW7mS/HvuxW02tz1RGwghuGVUIgO6BALwwbpsfmxhGw7zDuOl8S/RM6QnWCp5WWfhnn3zGvJ4nCi+Y8fS5Z238Zs8CYD6DRuwHlATxp9NqAr4JGPq2VMpMX+EUQ9gibwCabW2Nj+seg5MAZQmjWZf+T68dd5szi5n7/4+BDKIsjobU/6kuR/ag5dByxszLsNZOYofsr9hU1GTIps+ITD0ZoL91xNx/23UrVxJzrXX4Shpv3/wEVxuyZLdRSz5rejoOYFDkgk9515iQnui0+hg09u48za3PbadCCHwHTNGMUdIyeFnn6No9mw1L/FZhKqATxL1mzZz+OmncdXUNO+IG07Vzho0vr74DB3a2F70GxxYAufcxbmJk1k6YykRPhFEhVXgk/w8vtpg9FrB+B7hnXsiZxh9YgK4Z+A9uO0hPLhidrNMaYy4B7QGgiPTiX3jdWyHDnHoqpnHvIrUagSf/uUcnr+iH0IIyupszZL4HOH63tfz8LCHwW4mb/PrXLLybnaV7jrRUwQUZRz/wXyin/FUeLZYKJ8375hW9SqnH6oCPklY09Ko+XkpwtQkX29FFjJnI3UrVygrmSPJdwBWPw/GABh+W7N5tpVsQ2Mo50BOKCOSQ/E36TvpDM5cbh/Tg66am6h0FPHf9XMbO/wilAjEHZ/i559Hwicfg8tFzqxZ1K1dd0zvYdBp8DYopaD+8eUurnhrAw7XUTwVDN6YZ7xPcHA3onyioCoXWXLgeE+vAW1gIMYkpTBr3arVlDz/AtZ9+054XpVTh6qATxIhN91I8rKljRnOALZ/DB+cj7uyorn54fBe2LcIzrmdLdXpXPrdpWRUZlBSY2X+ylqCRE/ySo2q+aGdaDSCt6+4ClF7Dt9nfc724iarzklzIGUiLLobk3kjCV8sRB8bS95tt1H5xRfH9X7/mNKdf0zt3lBlo97W2ubbPWoIH17wCeHe4fDzIzz07WX8e+3JS1fpP20qSUsW4z14MACVX35J7YoVJ21+lc5BVcAngSO3gRpji3Lxw2+jyngVUmtonnxn9fNg8IPht6MRGiJ9I4n0iWRfcQ3ZedH4SMUpX1XA7SfC38QTox/G7fTj3l8fwXGkXL3eBDMXQNepsPh+9Lk/EP/pp/iMHEnxvx+j5MUXj9nntk9MQEMWtRX7Sxj93ArSCmuOOt553vN0SZ5KUmAKSInc9hHfpH3aOojkGDmy1yDdbqo++5zqb787oflUOh9VAZ8grro6MsZPoOLTT1t3+kVSsfIQPkOHovXzJNIp2Q97v1NMD97BDI4YzFuT3sLX4Et0WDU+3eZgNgcyMC6Q8D9h+aET4dL+yQz3v5VqVx7/Xv1KY4fOCFd9At0vgB//gXb3h3R58w0CZ15F+bvvUXD/A7it1qNP/DtEBZoY3z2clHBfAKotjlZjdAHR3DtpLjf2uREO7yHz5wd4bMsz/HLoFwDc0n1CG2tCoyFh4edEPj4HAMfhEgr+8eAx+T+rnBpUBXyiOJ0EXH45Xv0HNG+vLcbxy8s4C7Oamx9WPw96bxhxF7tKd1Fjb1w5bT28FaFxkVsY8afO/XsivHbpNRitQ1ic8zHbi9IaO3QGuPJD6HkxLH0EsekNIh97jPCHHqJ22TJyb7ixXSWNWtIj0p8Xr+yPQafB4XIz/c31zFm09+gHRPYl5aZf+PaChVyYfCFkpvLrikeZsWg6xfXHHtZ8BKHXo/PkH7am7aVuzRo1gOMMQFXAJ4g2MJCI/3sIrz69m3ccWot+3b/Re7saw49LD8Ker2HYX7EZffhb6t+YvVYpW59RUsdLix0YnT2QjmCm/slz/x4v3gYdr0x5DOkycdeyh3G6mthntXqY8T70vhx+mY1Y+z9Cbr6JmJdfwnrgAIeumoktK+u431tKmD4olrHdwgBwutzUWFuviIkeSEpoLyXd5e6vMe3/kSifKMK8lOM2Fm1ka/HW45bDb/x4uq5cgSEuDoDi/z5F+bx5xz2fSsehKuATwLJ3L9a0tLY7e19GQdZoRHRf9FGeqgtrXgC9F4y8B6PWyOsTX+dvg/8GQFmdjYo6id4ZRUq4H0lhvp10FmcfIxMTmBp5O3Vk8+Avrzfv1Orh8neh75Ww/AlY+Sz+U6YQ/9GHuC0WDs24gpKXX24eTt5ODDoNd4xLbnAdXLA5l/HPr6SgynL0gy55jTHXLeW1Sa+jlRIWXscL6x7n1R2vNgxJK0/D6jw2E8mRWoTS5cJZUoKzvNHe7KqtPdphKp2MqoBPgLLXXifvjjvbLIHjrKikZksWvhOVKCbKM2H3lzD0FqxGRbn2Du1NUoCS+CUipAZT4guUVwWpq9+TwHPTrsfP3Z9fij5kXU7z6hlodXDZWzDgGlj5FKT+F6++fUn8YiE+o0dT/uZbZEyaTOnrr+Oqa7soZ3sYFBfE9MGxRAcotvzf8quwOlqYBYRQ3OUAqnOh+Dc+6nkr/zn3P2CuwJG/lVuW3sKzW55tOKTe0f7qy0KrJfbllwj/h1Kw1JqWRvqo0cfshqfSMagK+ASIfvYZYl99BaHTNe9wOXEsvA+jv72x9PzqF0BrpHrITVz03UV8efDLhuFSSrYWbwPAUZ+o2n9PAlqthnfOewqBlnt//T+25TUvOYRGCxe/BoOuh9XPwfLH0UdHE/vySyR+9y3ew4ZR9uprZE6cRNnb7xxXyfk+MQH88/yeCCGwOlxc//5m/vnN7qMfEJwEd2/Fu9elxPnHwZ6vEe9N5MUB9zOz+0xwuyiqK2LU56Na50L+A4RG+VfX+PsTcPllePXrCyj5iGtTU1V78SlCVcAngNbfH69+/Vp3VOVgKlmMd6Ivxu7doSILflsIQ25G+IYxPHI4fUL6NAzfkVfFY58LhKUb4aZY+sYEtJ5T5ZjpExnHDd3+jl2bww2/zOCO797FYm9yt6LRwIUvw5BbYO1cWPYoSImpRw+6vP4aCV99hWlAf0rnziVj0mTK572P2/I75oTfwajT8Pa1g/nrGOWOp7TWxs0fbGF/cQv3Na1eWRUD9JmObvo8Rva6ku7B3WHZbDRf3sD1Pa+jd4iy57ChcAMzFs0guzqb9mCIjSXqscfQ+vsDUPHpAoofmwMeBayGOXcuqgI+DlzV1eTdfsdRo5DcXlEc/CEO+lymJN9Z8z/lH+vce/E3+PPkqCeVBC4etAKEMQ+nLYCpvSPRqLl/Txr/OHcm8yYvIFAXzdrqVxj5wSy+3tVkFarRwAUvwrDbYMNrsOQBsCr2X68+vYl7+20SPv8MU69elDz/PBmTp1Dx0UfHXEZICMHwpBB6RimKL7usnrTCGgyeYI6SGisV9S2Sx3sHQ98mlVXCuhERM5z7hzygrJBXPoMmfxtBpiAivBUzxtcHv+aOX+/A7GhfiHL0008R9+EHCIMBKSU5s66h/P35x3RuKsePqoCPA3t2trL5pmn78tWtW4fb7MBv0lSozIFdn1E/8Fr+tetVCusKW40PCaxHEzUPhzVaNT90AMNierHq2q+4KukeXIZMHtt+Ixd/9DT5VR77rhBw3rMw8h7Y+j7M7QO//Luh2oXXgAHEzXuP+E8+xpiUxOGnniZz8hQqPvr4uCsdD0sMZt3DExo2W19ans6451e0thE3ZfCNMMWT6tTlhN1fMbyyiHenvIu3zgu+uwtZkoZbuvHWewPw9q63eXbzs0edUuh0GJOUVbm0WDDExaELCQbAXV9Pyf/mYs/JOa5zVPljVAV8HHgNGEBK6vKGsuPNkBLNqv/g3w0lTHTt/0BoSOs5mdTcVErMzbNx2Z1u1uYpLkcmZwrDk4I74xT+dGg1Wh4dfSs/XLaIOO8+ZMsFTFt4Fc8sX4HT5VaU8JQn4dZV0HUyrH8VXuoL39+tVC0BvIcMIf6jD4n74AP0sbEcfuop0keP4ZBn1WjPzT1GmRrvdG4cmcDjl/RuSD16+8fbeO7n/Uc7VNlIvHsLjHtYeW2ugEOrmWGK5e3Jbyuv3z+PqtK9lFsb/ZvnrJ/D+3veb3NKjbc30c8+Q8AllwBg3b+f8nnzcBQpaTjt+flUffWV6kVxElEV8DHirKxEStl64w0lJLTy88/R1u7Ht08Mor4YdnwKg65naNI0ls1YxoDw5gEb6zPLeORjLa66rkxM7teQX0ClY4gPiGHJlR/w0KAn0Bsr+CTvPsa8+zAbsw8rA6IHKL7C92xXNuh2fwmvDYXPr4G8LQD4nDOc+E8/IfH77wi96y7cZjMlzz1H5pSpZF18CaWvvIo1Le2Y7KndIvy4bKCS/F1KSZCPAX8vfcPri19byycbW6xEhQCjJ8LSJwTu2w2DblBe15eBy87/dZ3Jc2Oeg/xt8MpAKqsOUWevg/pyZP42rvzhCj5O+7hhSoe70W/Ze/Bgum3c0JBvoi51BUWPzm7YkLQeOEj95s3qBt6JIKU8Kx6DBw+WncGha66Vubfd3qrdvnK+PDTrSpnWvYfMvXGmtBfkS/nD/dL8RIhcf3DRUefLKKmV/f43R3Z7+Wr542+FHSm6SgvKzeXy+kX3yT4f9JE93x4vb134haystzUfVFcq5fInpXw6TsrH/KWcN03KAz9L6XI1G2bLy5Nl8+fLQ9dcK9N69pJp3XvI9AkTZfFTT8m6TZuk2+k8bjlrLHZ55yfb5KKdBVJKKSvrbXLCCyvk8n3FUkop7U6XtNj/YP78rVJ+NkvKqjzl9Y4F0jInQD687E65OHOxlFmrZc0P98mhHw+RXx/8WkqnXdocVrm7dLc0O8xSSindbre0ZmY2TFn46KNy/5Ch0m23SymlrF6yRFYsWHDc53k2AWyV7dBbQp4lu55DhgyRW7cef/RQe5BSUv3NtwiDgYCLLgSUMkRVbzxJUNUrlB0MQ3/F8wRcdimiphBeGcC87ufykiWT7y75juTA5GZzZZbW4etTz+SvJuMsvYht9z7RkPJQpfNYfmg1/1wzB7O7FFkzlG5+I5iUNJIJ3brQPcJP2Ui11cGOj2H9a1CTD2E94Zw7IGksBMY3ei4AzooK6lasoPaXX6lfvx5pt6MJCMDUqyemnr0w9eqFqVdPDPHxCO2xVzvJKa/n8R/SuHtCCoPigthyqIIr3trAp38ZzrkpoeSU1/PTnmIuHxRDuN9R8onUHob8zdBtmrJBvGUe5amPM3/cnUxLvpA+aT+zd8NLzIwM4n/j/sdkm5v83HV86u/N1T1mEWcIxGlx4MjJw2uAcleXf+/fcBQWkviV4mKZf9/9ICB2rpIi1Lx9B7rgIAwJCcd8zmcaQohtUsohfzhOVcDHj3n1jxQ9/Rr27GzCLu1P4ANz0YV7ot5+fBC2vo/9rs2src9hQlzzahivLk/n1dQMHpvp4pntj9CHx/jshhltvItKZ2B2mHl87Yv8nPstbhxIqcFticXg7Eaf4MGcl3IOY7tFE+uvhz3fwLqXocST88E3AmKHQpfh0GUYRA1QsrChbGTVrVlL/bq1WNP2YTt4EOlQbvOFlxem7t0x9eqJsWdPTD16YuzWtXVWvT8gt9zMdzsLmDm0C+H+Jn7YVcg9n+3g1wfGkhLuy6JdhTy1ZB9f3j6CLsHebM+tZPXBUv4yOglfo46SWivVZgdJoT5oj5jAMldQk72SLV3H0D+8P6Gr57J+35fcFxbIh9M+pOf6t0nN+pE5EZHMnzaf5PSVpJensS2qPxf0vAy/8mxKPv0W4R9J2K1KzuvMCy7EmJRE7KtKoqSc667H1LcvEQ89CEDZu+9iTE7Bb8J4AGzp6WhDQtAFn3n7Iu1VwOpyq504y8upX7sW//POw22xUPviTfi7l6PT9SXi3XfxHT2qcXBNEY5tH0L/mRhCkpkQkozZ7mTBplzGeTJnXTIghgh/E2uK3ka6TFzS9w8/K5UOxFvvzbPjZ/OE6yF2lezi10PrWJO3gQJLKrucv7IzTceT2+Lwkz0ZFD6U80d+zii/EoLKdyLyN0PeZti/WJlMo1dsybHD0HQZhv/IYfhPmwqAdDiwZWZiTduHdd8+rPvSqP5+Ee4FnynH6nQYExPRx8aij45WHjFH/sagDQ5uXlcQiAvx5t6JXRteX9Q/mrHdw/Dx3E1F+psYmRJCmJ+i2HfmVvHSr+n8dbTi/fDl1nyeX3qAfU9Mw0sL76/NZsFmDUvvm8NEjWDp3mK2OK7hX3c/xkatjv3FNewJGEv4wFAm6JwEGkKQh/eyuWwrzxxexpTuF0Lqf/glMIcXjU5W2GYSsOhvlE52sifGlxvcDvSb38MZaKUiWE+4lIj8LVTMfx//8y9UFLCUHLpqJoFXzCDikUcAyJg4iaBZVxNyyy1IKSl69FH8Jk/Gb9w4pNtN7dKlmHr3xhAXp9zi22xoTKd3RsEOVcBCiGnAy4AWeE9K+UyLfiPwETAYKAeuklIe8vQ9AtwCuIB7pZRLO1LWP6JmyY8cfuopnOVllL//AZjL0E8/ly6ffoLG3/ML7XZDTQGsepZ3/L1Y7srhI3sdvgZfLHYXLyw7QLW1nhG9akmrSGNf/T7WFa/BZUlgcs+oU3l6Kh6MWiPDooYxLGoYjLifekc924q3sTRrLRuLNlFi+5F15iWs3W7AZYlD4w7CTxdMsP9M4sO96S9q6O8spEfNfgK3zENs9OSi8IuCgFiEXxQm/2hM4VGQkgKzRiN9InDUuLFmHGpYJTsKCjBv2YK7RSi0MJnQR0U1KufoKLRBwWgDA5VHUCDawCB8gwIb/MmHJQYzLLFxFXnzqESuGxGPztN/Xp9I4kO88TIo5pAIfxO9o/0bvDT2FlSz+LciHr2wFwJYsCmPH3f7sn32o/QBHvnmN1L3n8+GR15gqrWCBevLsbivZdI5Zq63HWZThoUYktiTJJlXt4Iplf8iYu9ivhtu5jM+YQt/gx/+xtI7vFjnt5PPAfn6ORRfHMq2Xm5mAfLr23EkeVEQ4CYEkGvfpH7VCiXQCXAfXEvB/Q8Q/vD/EXLjjbjKSkgfPY6I2Y8SfM01OCsrKXjgAUJu+Qu+o87FVV1N1Vdf4zt+PMakRNwWC7aMTAwJ8Y2pYzuBDjNBCCG0wEFgMpAPbAGullKmNRlzJ9BPSnm7EGImcJmU8iohRC/gM2AYEA38CnSTUh51u/VkmCCclZXY9u3Da/BgNEYjdevWUblgAZGzZ+OqraX2mWswyhzKK4YSdfvlmIIcStn5iiyozFZ8fl2Kg/6qPhfwkjMGTc1EhvVLZ29ZGgfLCqkXjdm2DDIMe30kCfoL+OHWa05IdpXOocpaxcaizfycuYY9ZXuotpdjlVVA6/8j6fTB6DQSLgUx0k6gdBDgthLoqsfPZcdLSrykxOR24yUlBo0XOlMwRlMoWoMPeoMPOpcJXb0GUQ+ixgnVdtwVFhwVtThKq3BVHd0lTHh5oQ0MQBcYhDYoCG1gAMLbG42XNxqTCY23F8LLC43JS3luMil9nudCr0fo9Ai9DqHzPPR6iswuym1u+sUFI7Ralu87TE65mZtHKeWS/rfsAJll9bw+axAAN3+whZJaK4vvGY2Ukhvmb6HabOfl62LJr8vn41QT8Y5M+gyoYI+jElPNDEYWfczSkFy22iu4JeEtRm2+nWeCJQd0gjmDPmTApwN4MKk7eT6BvDn+E4LnJvBu0ABqug7iwan/xfVoGCtr+1M75XwuPf9+HM+NYN8qLfZbb2DIRTdhfesGsl/dgeGZR0m+cAbWBf8i+6klhL78AmETJmJZ9Do5T35Cl9dfx2fkyGP+npxyG7AQYgQwR0o51fP6EQAp5dNNxiz1jNkghNABxUAY8HDTsU3HHe39jlUBV5TksfjeS0nIMZMf7U2XQguHI4z03m8lPdFEYq6Vg129CCuxYDEZiSu0sWGagW+7a/B1gMlQR44jgVKDi2RXPVqdiwJHDHX6GgIxUauvwVI5AOkIxBC+GOkIwW2NxmWNQWOPJdSQSKRfMOF+Rm4YmcA5SSHHfI1VTg9cbheVtkpKzCXkVReTWVHIoepiiuoOU2YupdpRjsVdjVvacAs7iBMvW6+VEp2UmJwSfzP4WST+FvC1gL8FfDzPfS3gY5H4WMHbCkY76J1gcAi0x1YIpE3cgFsLbgGyyaPhtQbsQotVZ8DbbUEK2BeQhEOro2/lQSSwpMs4QHJ+wSqkELzS83q8nRZuzvgSCcwe+HciLaVcWfQRZiO81eVBUupyGCgXUq/X8IPXw1zuWEGZ/wbqpYbNtof5p/sTvo3Ow+nQUFD6N17Tv8acLg4MdRpqi/7Kx37PcUeoP4Y6DXVlN/Gt35M8ZI+iRqehwnw9nxme5b38CDIiBcPtUVz1xvJjui6ngwKeAUyTUv7F8/o6YLiU8u4mY/Z4xuR7XmcCw4E5wEYp5See9nnAT1LKr1q8x63ArZ6X3YFjqXwYChxfGNPJ5XSQQ5WhkdNBDlWGRk4HOY5HhngpZdgfDepIG3BbCQ1aavujjWnPsUgp3wHeOXbRQAixtT2/UB3N6SCHKsPpJYcqw+klR0fK0JFhV/lAlyavY4GWiRAaxnhMEAFARTuPVVFRUTmj6UgFvAXoKoRIFEIYgJnAohZjFgGe2ElmAKmeKJJFwEwhhFEIkQh0BTZ3oKwqKioqnU6HmSCklE4hxN3AUhQ3tPellHuFEE+ghOktAuYBHwshMlBWvjM9x+4VQnwBpAFO4K7f84A4To7LdNEBnA5yqDI0cjrIocrQyOkgR4fJcNZEwqmoqKicaaipt1RUVFROEaoCVlFRUTlFnDUKWAgxTQhxQAiRIYR4uI1+oxBioad/kxAioUnfI572A0KIqe2d82TJIISYLITYJoTY7fk7ockxKz1z7vQ8wjtIhgQhhKXJ+7zV5JjBHtkyhBCviJbJCE6uHNc0kWGnEMIthBjQQddijBBiuxDC6fFbb9p3gxAi3fO4oUn7MV2L45VBCDFACLFBCLFXCPGbEOKqJn0fCCGym1yHAS3nPcnXwtXkvRY1aU/0fHbpns/S0EHXYnyL74RVCHHp8VyLdsjwgBAizXPNlwsh4pv0nZTvRDPak7PydH+gbPJlAkmAAdgF9Gox5k7gLc/zmcBCz/NenvFGINEzj7Y9c55EGQYC0Z7nfYCCJsesBIZ0wnVIAPYcZd7NwAgU/+yfgPM6So4WY/oCWR14LRKAfij5SGY0aQ8Gsjx/gzzPg471WpygDN2Arp7n0UAREOh5/UHTsR15LTx9dUeZ9wtgpuf5W8AdHSVDi8+mAvA+1mvRThnGN5n7Dhr/P07Kd6Ll42xZAQ8DMqSUWVJKO/A5cEmLMZcAH3qefwVM9PxSXQJ8LqW0SSmzgQzPfO2Z86TIIKXcIaU84ue8FzAJJVHRsXIi16FNhBBRgL+UcoNUvm0fAZd2khxXo+QEOR7+UAYp5SEp5W8oEbVNmQr8IqWskFJWAr8A047jWhy3DFLKg1LKdM/zQqAEJUz/eDiRa9Emns9qAspnB8pn2SHXogUzUKJi21d19NhlWNFk7o0oMQhw8r4TzThbFHAMkNfkdb6nrc0xUkonUA2E/M6x7ZnzZMnQlOnADill07K78z23V7P/4PbmRGVIFELsEEKsEkKMbjI+/w/mPNlyHOEqWivgk3ktjvXYY70WJyJDA0KIYSgrtswmzf/13CbPbceP9YnKYRJCbBVCbDxy64/yWVV5Prv2zHlSrgXK3VLL70R7r8WxynALyor29449nv+PBs4WBdwRYc/tCoc+STIonUL0Bp4FbmvSf42Usi8w2vO4roNkKALipJQDgQeABUII/3bOeTLlUDqFGA6YpZR7mvSf7GtxrMd2xHfi9wVRVlgfAzdJKY+sDB8BegBDUW6J/++PpjlBOeKkEoo7C3hJCJF8HHOerGvRFyW24AjHci3aLYMQ4lpgCPD8Hxx7Qud1tijgjgh7PtZw6BORASFELPAtcL2UsmGlI6Us8PytBRag3EaddBk8Jphyz3ttQ1ltdfOMj21yfHvCwk/oWnhotdLpgGtxrMce67U4oZB6zw/gEuBRKeXGI+1SyiKpYAPm8/vX4YTlOGIek1L+f3v3E2JjFMZx/PuMGDXUjIaywq1JGWqmKLGh2ExITawsRElsrI2VLBgpYmlBWQglspBQsqD8yQwWcg0bhVhp1DSLx+KcMe/c+ePOzH2dMe/vU2/3vXfee+7TOe88895z5j2nj9AP306YnKYxtl01ZdZieoHdwE13/7Ny6CTroqoYzGwL0AXsyHwTrdU5MVK1ncUzeSPc0ddHGEQb6lxvrTjmMCMHfa7F/VZGDsL1ETrr/1pmDWNojMd3jlFmc9yfS+hvO5hTDIuBOXG/BHwGFsXnz4D1DA8ydOTVHvF5HeHELuVZF5ljLzF6EO4jYbClKe5Pui6mGcM84AFwZIxjl8ZHA84CJ6fbHhPE0QTUx/1m4D1x4Aq4zshBuEN5xJB5/Smweap1UeV52U64+GipeL0m58SomKo9cKZyJsLTAAACB0lEQVRvQAdhAvgPQFd87TjhrxjA/HjClAmjltlf7q74vndkRjDHKjOPGIBjQD/wKrMtARqAF0AvYXDuHDFJ5hBDZ/yMHuAlsD1T5lrgTSzzAvEOyhzbYxNhOtJseXnUxTpCou8nrMjyNvPefTG2MuHr/5TqYqoxAHuAwYpzoi3+7CHwOsZxBVhQg/YYL44N8bN64uP+TJml2Hbl2Jb1ObbHcsJFQV1FmZOqiypiuA98zdT57VqfE9lNtyKLiCQyW/qARUT+O0rAIiKJKAGLiCSiBCwikogSsIhIIkrAIiKJKAGLiCSiBCyFZGZn4ryv51PHIsWV26KcIjOVmZWAje6+KnUsUmy6ApZCMbOVwCNgWZx6syF1TFJcuhVZCsfMTgCf3P1i6lik2HQFLEW0hjC5zITMbK+ZbfsH8UhB6QpYCsfMyoS1x3YSZl77CRwlLIH03d3vmNlV4C6wFfgCfHP3U2kiltlKg3BSKGa2EBh0919xEvxe4Ja7D4yzwtE9d79sZjfMrM6HV6UQmTZ1QUjRrCbM3Yq7dwOPgdNm1gIMMHxRMjQ455lHfV2UmtIVsBSKuz8BdgGY2QGghbAK7w/Cf0d0m9kKwiolEFa+bQOeu/rrpMbUBywikoi6IEREElECFhFJRAlYRCQRJWARkUSUgEVEElECFhFJRAlYRCQRJWARkUR+A+qCcLmGpS1ZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(5,5))\n",
    "\n",
    "for i, (posterior, label) in enumerate(zip(posteriors_1d, labels)):\n",
    "    plt.plot(f_sub_1d, posterior, ls=\"-\", label=label, c=\"C{}\".format(i))\n",
    "\n",
    "for i, (posterior, label) in enumerate(zip(posteriors_1d_flatprior, labels)):\n",
    "    plt.plot(f_sub_1d, posterior, ls=\":\", c=\"C{}\".format(i))\n",
    "    \n",
    "plt.legend()\n",
    "\n",
    "plt.xlabel(r'$f_{\\mathrm{sub}}$')\n",
    "plt.ylabel(\"Posterior density\")\n",
    "#plt.yscale(\"log\")\n",
    "plt.ylim(0., None)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"../figures/posterior_1d.pdf\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Check integral over posterior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9999999971562625\n",
      "0.9999999996837958\n",
      "0.9999999997503995\n",
      "0.9999999998882961\n",
      "\n",
      "0.9999999987547427\n",
      "1.0000000000943188\n",
      "1.0000000007021506\n",
      "0.9999999998595811\n",
      "\n",
      "0.9999999971562625\n",
      "0.9999999996837959\n",
      "0.9999999997503994\n",
      "0.999999999888296\n",
      "\n",
      "0.9999999987547428\n",
      "1.0000000000943188\n",
      "1.0000000007021506\n",
      "0.9999999998595813\n",
      "\n"
     ]
    }
   ],
   "source": [
    "for posteriors_ in [posteriors, posteriors_flatprior, posteriors_1d, posteriors_1d_flatprior]:\n",
    "    for posterior in posteriors_:\n",
    "        print(np.sum(posterior))\n",
    "    print(\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (lensing)",
   "language": "python",
   "name": "lensing"
  },
  "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}