{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Bayesian signal reconstruction with Gaussian Random Fields (Wiener Filtering) - CMB example"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Florent Leclercq,
\n",
"Imperial Centre for Inference and Cosmology, Imperial College London,
\n",
"florent.leclercq@polytechnique.org"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This notebook requires the CLASS code (https://class-code.net) and its python wrapper classy, as well as the Planck colormap (https://raw.githubusercontent.com/zonca/paperplots/master/data/Planck_Parchment_RGB.txt).\n",
"Optionally, pre-computed data are available:
\n",
"* http://www.florent-leclercq.eu/data/Bayes_InfoTheory/sqrtsignalcovarPix.npy (2.4 GB)\n",
"* http://www.florent-leclercq.eu/data/Bayes_InfoTheory/sqrtCovWF.npy (2.4 GB)\n",
"* http://www.florent-leclercq.eu/data/Bayes_InfoTheory/CovWFinvnoisecovarmat.npy (2.4 GB)\n",
"* http://www.florent-leclercq.eu/data/Bayes_InfoTheory/invnoisecovar.npy (98.4 kB)"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import scipy.linalg\n",
"import os.path\n",
"import matplotlib.pyplot as plt\n",
"from classy import Class\n",
"import healpy as hp\n",
"np.random.seed(123457)\n",
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# Download and define the Planck color map\n",
"from matplotlib.colors import ListedColormap\n",
"if not os.path.isfile(\"data/Planck_Parchment_RGB.txt\"):\n",
" !wget https://raw.githubusercontent.com/zonca/paperplots/master/data/Planck_Parchment_RGB.txt --directory-prefix=data/\n",
"planck = ListedColormap(np.loadtxt(\"data/Planck_Parchment_RGB.txt\")/255.)\n",
"planck.set_bad(\"gray\") # color of missing pixels\n",
"planck.set_under(\"white\") # color of background"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup Cls using CLASS"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# Define cosmology (what is not specified will be set to CLASS default parameters)\n",
"paramsLCDM = {\n",
" 'h':0.702,\n",
" 'n_s':0.9619,\n",
" 'Omega_b':0.045,\n",
" 'Omega_cdm':0.272-0.045,\n",
" 'output':'tCl lCl',\n",
" 'l_max_scalars': 1000,\n",
" 'ic':'ad',\n",
"}\n",
"# Create an instance of the CLASS wrapper\n",
"cosmoLCDM = Class()\n",
"# Set the parameters to the cosmological code\n",
"cosmoLCDM.set(paramsLCDM)\n",
"# Run the whole code.\n",
"cosmoLCDM.compute()"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"# Compute the power spectrum for the current set of cosmological parameters\n",
"res=cosmoLCDM.raw_cl(1000)\n",
"ell=res['ell']\n",
"Cl=res['tt']"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"# Clean CLASS\n",
"cosmoLCDM.struct_cleanup()"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"# Define a function of l from the arrays\n",
"from scipy.interpolate import InterpolatedUnivariateSpline\n",
"Clfunc = InterpolatedUnivariateSpline(ell, Cl, k=2)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAa4AAAEWCAYAAAA+bHOCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xd8VfX9+PHX+97sDQkrgwTCDATCCsvBkqpF3OKsipXab9Xaoa1tv2r9dVhr/TqrOHEianFA3QhSkA1hbwgQSCAJZO/cz++PG2gIGfdmnXuT9/PxyEPvOZ/z+bzDgbzzOeczxBiDUkop5S1sVgeglFJKuUMTl1JKKa+iiUsppZRX0cSllFLKq2jiUkop5VU0cSmllPIqmriUUkp5FU1cSimlvIomLqWUUl5FE5dSSimv4mN1AB1RVFSUSUhIsDoMpZTyKhs2bMgxxnRrqpwmrjaQkJDA+vXrrQ5DKaW8iogccqWcPipUSinlVTRxKaWU8iqauJRSSnkVfcfVTiorK8nIyKCsrMzqUJQbAgICiI2NxdfX1+pQlFI1NHG1k4yMDEJDQ0lISEBErA5HucAYQ25uLhkZGfTp08fqcJRSNfRRYSsSkctE5KX8/PxzzpWVlREZGalJy4uICJGRkdpLVsrDaOJqRcaYRcaYOeHh4fWe16TlffSeqcZUVZSzd+Mytiz9kKzDe60Op9PQxNWJZGVlcf3115OYmEhSUhKXXnope/bsIT09naFDh9Z7TVVVFVFRUTz44INnHV+8eDEjRoxg+PDhJCUlMXfuXAB2797NpEmTSElJYfDgwcyZM6dVYn/qqacoKSlp9vVpaWl89tlnrRKLUsbhYM0H/yD3L4Po/+nlDPvuDnq+Npq0xy8h59hBq8Pr8DRxdRLGGK688komTZrE/v372bFjB3/5y184fvx4o9d99dVXDBw4kPfffx9jDOAcaDJnzhwWLVrE5s2b2bRpE5MmTQLg3nvv5Re/+AVpaWns3LmTe+65p1Xi18SlPEVZSREbnryKsdsf5ZRvT9aPeZKdl7zP93FzGFC8gaqXpnH0wE6rw+zQNHF1EkuXLsXX15e77rrrzLGUlBTOP//8Rq+bP38+P//5z+nduzerV68GoLCwkKqqKiIjIwHw9/dn4MCBAGRmZhIbG3vm+uTk5HPqNMZw//33M3ToUJKTk1mwYAEAy5YtY8aMGWfK3X333cybN49nnnmGY8eOMXnyZCZPngxASEgIv/rVrxg5ciRTp04lOzsbgEmTJp1ZtSQnJ4eEhAQqKip46KGHWLBgASkpKWfaU8pdFeXl7H72akYWLmNN4r0MfHAlo394B4PH/oAJd/ydzKs/IZBSHG9dRVFhntXhdlg6qtACf1y0nR3HClq1zqToMB6+bEiD57dt28aoUaPcqrO0tJQlS5Ywd+5c8vLymD9/PuPHj6dr167MnDmT+Ph4pk6dyowZM7jhhhuw2Wz84he/YMqUKUyYMIHp06dz++23ExERcVa9CxcuJC0tjc2bN5OTk8OYMWO44IILGozj3nvv5cknn2Tp0qVERUUBUFxczMiRI/nHP/7Bo48+yh//+Eeee+65eq/38/Pj0UcfZf369Q2WUcoVaS/NIbV0NWuG/J6x1z1wzvnEYePZWfwKA7+4kfWv303qvW9bEGXHpz0u1aDFixczefJkgoKCuPrqq/noo4+orq4G4JVXXmHJkiWkpqbyxBNPMHv2bABuv/12du7cybXXXsuyZcsYN24c5eXlZ9W7YsUKbrjhBux2Oz169ODCCy9k3bp1bsVms9mYNWsWADfffDMrVqxohe9YqYatW/Qyqbkfs6rXzfUmrdMGj7+Utb1uJPXkInZvWNZ+AXYi2uOyQGM9o7YyZMgQPvzwQ7eumT9/PitXruT0Sve5ubksXbqUadOmAc7HgMnJydxyyy306dOHefPmARAdHc3s2bOZPXs2Q4cOPae3d/pdWV0+Pj44HI4zn90Zhn569F/tOnQYu2ot2UcPMnjD/7LLN4nRtz/ZZPmhN/6Jk08uourL/4VR/2mHCDsX7XF1ElOmTKG8vJyXX375zLF169bx3Xff1Vu+oKCAFStWcPjwYdLT00lPT+f5559n/vz5FBUVsWzZsjNl09LSiI+PB+CLL76gsrIScI5izM3NJSYm5qy6L7jgAhYsWEB1dTXZ2dksX76c1NRU4uPj2bFjB+Xl5eTn57NkyZIz14SGhlJYWHjms8PhOJOI3333Xc477zzAuTL/hg0bAM5K1HWvV8odR9/9GXZTTfD1r+Dr599k+ZCwruwZ8BOGVGxhx9pv2iHCzkUTVychInz00Ud8/fXXJCYmMmTIEB555BGio6MB5zD22NjYM19z585lypQp+Pv/9x/p5Zdfzqeffkp1dTWPP/44AwcOJCUlhYcffvhMb+urr75i6NChDB8+nB/84Af8/e9/p2fPnmfFcuWVVzJs2DCGDx/OlClTePzxx+nZsydxcXFcd911DBs2jJtuuokRI0acuWbOnDlccsklZwZnBAcHs337dkaNGsW3337LQw89BMCvf/1rXnjhBSZMmEBOTs6Z6ydPnsyOHTt0cIZyW9o380kpXkla4l3EJbr+tGTYzLspIJiy5c+0YXSdkzT02EY13+jRo03d/bh27tzJ4MGDLYqo4wkJCaGoqKhd2tJ713lVVVZw9K8jEAy9frvJpd5Wbavm/ozUY++QO2cz3WPi2yjKjkNENhhjRjdVTntcSinVgA0fP0u8I4Pccb9zO2kB9J5yJ3Yx7Fs6r/WD68Q0cSmv1F69LdV5lRTlk7j9GXb6JpEy7cZm1RHTP4W9Pv3pdvCTVo6uc9PE1YoaW2RXKeVdtnz0D6LIQy56FLE1/0flycQr6F+9n/SdG1oxus5NE1cramqRXaWUdygrLabf/jfY6j+SQakXtaiuxMm34jBC5ur3Wyk6pYlLKaXq2LzoeaLIw3bBr1tcV1TPOPb59qfr0aWtEJkCTVxKKXWWyopy4na+zC6fwSSNv6RV6syNmUL/yj3kHs9olfo6O01cncif//xnhgwZwrBhw0hJSWHNmjWNln/kkUd44oknAHjooYf45pvWmUh5//33M2TIEO6///5WqQ/OXf39008/5bHHHmu1+lXnsfmL14k2Jygbf1+L3m3V1m3kZdjEcGDVx61SX2enSz51EqtWrWLx4sVs3LgRf39/cnJyqKiocPn6Rx99tNVimTt3LtnZ2WdNbm6ptLQ01q9fz6WXXgrAzJkzmTlzZqvVrzqPsC2vctgWw7BJ17ZanYnJE8j+qAu2fV8Dd7davZ2V9rg6iczMTKKios4ki6ioqDOrZiQkJPCb3/yG1NRUUlNT2bdv3znX33bbbWeWUEpISODhhx9m5MiRJCcns2vXLsC5Yvvs2bMZM2YMI0aM4JNPzh0CPHPmTIqLixk7diwLFiw4q15wTiwG5xYnkyZN4pprrmHQoEHcdNNNZ9Y4XLduHRMmTGD48OGkpqaSn59/zrYl8+bN4+67nT8gDh06xNSpUxk2bBhTp07l8OHDZ76ne++9lwkTJtC3b1+313JUHc/u9d8yoGoPmQNuwWa3t1q9YrNxKGw0CUWbMLXW41TNoz0uK3z+W8ja2rp19kyGSxp+NDZ9+nQeffRRBgwYwLRp05g1axYXXnjhmfNhYWGsXbuWN998k/vuu4/Fixc32lxUVBQbN27kn//8J0888QSvvPIKf/7zn5kyZQqvvfYaeXl5pKamMm3aNIKDg89c9+mnnxISEkJaWhoAn3/+eYNtbNq0ie3btxMdHc3EiRNZuXIlqampzJo1iwULFjBmzBgKCgoICgo6Z9uS00tQgXNfrx/96EfceuutvPbaa9x77718/LHzkU1mZiYrVqxg165dzJw5k2uuuabR71t1bIXLn6fIBDLk0ruaLuwmR/xEIrd+zaE9acQPGtnq9Xcm2uPqJEJCQtiwYQMvvfQS3bp1Y9asWWf9cL/hhhvO/HfVqlVN1nfVVVcBMGrUKNLT0wHnOoWPPfYYKSkpTJo0ibKysjO9m+ZITU0lNjYWm81GSkoK6enp7N69m169ejFmzBjAmXB9fBr//WvVqlXceKNzAuktt9xy1hYoV1xxBTabjaSkpCZ3g1YdW07WYYblL2Vb9xmEhHVp9fqjhzuH1WdtWdJESdUU7XFZoZGeUVuy2+1MmjSJSZMmkZyczBtvvMFtt90G/HdbkLr/35DTjxztdjtVVVWAc7uSf/3rX2d2Q3ZF7W1IjDFnvXer/Q7sdDvGGJfia0zt62u3oet2dm57P3uO8VJNzPR72qT+mL5JnKArPkdWAq03MKkz0h5XJ7F792727t175nPtrUiAMyumL1iwgPHjxzerjR/84Ac8++yzZxLApk2bmrym9jYkn3zyyZktURoyaNAgjh07dmbjycLCQqqqqhrdtmTChAm89957ALzzzjtntkBR6rTKinL6HX6fLQFjiOs/vE3aEJuNw2EjiS/U91wtpT2uTqKoqIh77rmHvLw8fHx86NevHy+99NKZ8+Xl5YwdOxaHw8H8+fOb1cb//u//ct999zFs2DCMMSQkJDT5ruzOO+/k8ssvJzU1lalTp571Pqw+fn5+LFiwgHvuuYfS0lICAwP55ptvmDx58pnHlA8++OBZ1zzzzDPMnj2bv//973Tr1o3XX3+9Wd+f6ri2Ln2fkZzi6Ogft2k71TFjiCr4hqyM/fTs3b9N2+rIdFuTBohIX+D3QLgx5hoRCQb+CVQAy4wx7zR0rbdta5KQkMD69euJioqyOhSP5Mn3TrWOzX+7iF6l++j6+934+Pq1WTt7Nn7HgE9nsnHsU4y85PY2a8dbed22JiISISIfisguEdkpIs16XiUir4nICRHZVs+5i0Vkt4jsE5HfNlaPMeaAMeaOWoeuAj40xtwJ6AQhpTqIExkHGFqyjv0xl7dp0gKITxpDhbFTfnh904VVgzzpUeHTwBc1vRs/IKj2SRHpDpQaYwprHetnjKk76Wge8BzwZp3r7cDzwEVABrBORD4F7MBf69Qx2xhzos6xWOD0GPZqN783j3Z6VKBSndGBr+fSXQy9p85p87b8A4LY45tIaO6WNm+rI/OIHpeIhAEXAK8CGGMqjDF5dYpdCHwiIgE119wJnLMntjFmOXCynmZSgX01PakK4D3gcmPMVmPMjDpfdZMWOJNdbM3/e8Sfm1KqZRzV1fQ+vJCtfinE9E1qlzbzuiTTp3wPVU0MRFIN85QfwH2BbOB1EdkkIq/UvFM6wxjzAfAF8J6I3ATMBq5zo40Y4Eitzxk1x+olIpEi8iIwQkQeBBYCV4vIC8CiBq5pdD8ufZ/offSedWw7v19MtDlB+bCb261NW+wogqWMw3s3t1ubHY2nJC4fYCTwgjFmBFAMnPMOyhjzOFAGvADMNMa4sw1ufZN/GvypZIzJNcbcZYxJNMb81RhTbIy53Rjz04YGZjS2H1dAQAC5ubn6g9CLGGPIzc0lICDA6lBUGylbO488Qhg6tXk7HDdH98ETAcjZ/X27tdnReMo7rgwgwxhzernyD6kncYnI+cBQ4CPgYdxbrTIDiKv1ORY41qxomyE2NpaMjAyys7Pbq0nVCgICAoiNjW26oPI6+SezSS5YzsbuVzIusPFpGK0pJjGZYuOPyWzlZd86EY9IXMaYLBE5IiIDjTG7ganAjtplRGQE8DLwQ+Ag8LaI/MkY8wcXm1kH9BeRPsBR4Hqg3X7N8vX1pU+fPu3VnFKqCTuXvMU4qSJq4o/atV273c4R3z6E5O1q13Y7Ek95VAhwD/COiGwBUoC/1DkfBFxrjNlvjHEAtwKH6lYiIvOBVcBAEckQkTsAjDFVOHtoXwI7gfeNMdvb7LtRSnm0kL0fccQWQ+Kw9l9JJT9sADEVB3QFjWbyiB4XgDEmDWhw4pkxZmWdz5U4e2B1y93QSB2fAZ81dF4p1TlkHt5HUvlW1iXMIa6VNot0S/ckIk5+SnbmIbrF6JMYd3lSj0sppdrFwWVvYhND7wtvtaT90PgUAI7t0YnIzaGJSynVqRhj6JH+Kbt9B9Gr7xBLYogZ5Hy4VHJEJyI3hyYupVSnsn/7OhIdB8nvd6VlMYR36UYWUfhk72i6sDqHJi6lVKdyYuVbVBkb/SffYmkcxwMTiSze23RBdQ5NXEqpTqO6upo+mZ+xI2g0Xbo3uHBOuyiJGEBcdQZVFeWWxuGNNHEppTqNHWu+ohc5VA65xupQsHcfhK9Uk3Vot9WheB1NXEqpTqN43XxKjD9Jk663OhTC4px7vOWkn7MDk2qCJi6lVKdQVlbKoFNL2Bl+PoEh564n2t569R0GQFmWrqDhLk1cSqlOYfvyhURQhN/IBtcoaFfhXbuRQwS2XB2g4S5NXEqpzmHzAk4SxuCJnrOB+Qm/OEKL060Ow+to4lJKdXj5eScZUvQ9+7pNx8fXz+pwzigK6UPPyiO63ZGbNHEppTq8Xd++TYBU0mXcTVaHchYT2Z8uFHIyO9PqULyKJi6lVIcXtGshR6Un/UZMsjqUswRFO0cWZh7QvbncoYmrFYnIZSLyUn5+vtWhKKVqnDh6iKTyNI7E/hCxYiX4RkT1GQpAUYYu/eQOz7qLXs4Ys8gYMyc83Pqhtkopp/1L52EXQ+wF1qwE35gesf0pN76QoyML3aGJSynVoUUd/IS9Pv2J7T/c6lDOYfPx4Zi9FwEFB60Oxato4lJKdVjpuzbRv3o/J/tebnUoDcrzjyWiLMPqMLyKJi6lVId17D9vUG2ExCme95jwtLLQOHpUZ2EcDqtD8RqauJRSHZJxOIg/9hk7AkcQ1bO31eE0SLr2JVAqyD1+xOpQvIYmLqVUh7Rr/bfEmOOUD7ra6lAaFdijHwA5h3XNQldp4lJKdUj5a9+hzPgyaMqNVofSqK6x/QEoyNxncSTeQxOXUqrDqSgvZ2DO1+wIO4+QsK5Wh9OoHr0H4jBCdc4Bq0PxGpq4lFIdzvb/fEQXCrEPn2V1KE3y8w/ghC0Kn/x0q0PxGpq4lFIdTvXmBeQRQtIFV1odikty/aIJLT1qdRheQxOXUqpDKSrMI6lgJbsjp+HrF2B1OC4pCe5Nt6pjVofhNTRxKaU6lB3fzidIyglP9exBGbWZiHgiySf/1EmrQ/EKmriUUh2K/84PyZRuDBwzzepQXObXLRGA44d3WxyJd2hR4hKR/iJyo4hMFBHfes4HtaR+pZRyR3bWYYaUbuRQ9A8Rm93qcFwWGu0cEl+UqYvtuqKlPa4ngMnAL4GDIvKqiKTUOj9dRP7YwjaUUsol+799Cx9xEH3+LVaH4paouEEAVOiQeJe0NHE9D0wBcoDfAGnAMyLykYj0Bz4DLmphG0op5ZKuBz5mv70vvQeNtjoUt4R1iSLfBGPPS7c6FK/QosRljPkKGAJ8DkwA7gCSgEuA7UAmsLCFMSqlVJMO7d3KgKo95PTx3JXgGyIiZPv0IKBYh8S7wqelFRhjyoCPa74AEJHImrovA6a2tA2llGrKse/mEWeExMmeuxJ8Ywr8exFZdtjqMLyCSz0ucbpLRMa5Ut4Yk2uMOW6MeQW4vUURWkRE+ta8s/uw5nOwiLwhIi+LyE1Wx6eU+i/jcBB3dDE7AlKIiuljdTjNUhEcQ7fqE7q9iQtcfVR4I873WRPdbaCmR+YSEbGLyCYRWexuO7XqeE1ETojItnrOXSwiu0Vkn4j8trF6jDEHjDF31Dp0FfChMeZOYGZz41NKtb7dG74l1mRRNtizV4JvVEQcQVJO4alsqyPxeK4mrpuA3cCTjRUSkZ+IyMMi0qWZ8fwc2NlA3d1FJLTOsX71FJ0HXFzP9XacyfcSnO/hbhCRJBFJFpHFdb6611NvLHB6w5xql78jpVSby1vzNqXGj8FTvPdhiF9UAgAnMnRIfFNcTVyjgQ+MMaaJcl8D/wtc424gIhIL/BB4pYEiFwKfiEhATfk7gWfqFjLGLAfqm36eCuyr6UlVAO8BlxtjthpjZtT5OlHP9Rk4kxfoxG2lPEZFeRmDalaCD/bwleAbE9azLwAFWQctjsTzufoDOBxIb6qQMeYA8A0woxmxPAU8ANT7gNcY8wHwBfBezTum2cB1btQfw397TOBMRDENFRaRSBF5ERghIg/iHB15tYi8ACxq4JrLROSl/Px8N8JSSrXEjuULiaAIn5E3WB1Ki3SLdT5AqsjRxNUUV0cVngQiXSy7HLjNnSBEZAZwwhizQUQmNVTOGPO4iLwHvAAkGmOK3GmmviobaSsXuKvO4UYHmhhjFgGLRo8efacbcSmlWsCx+T1OEsaQid43DL62sC7dKTH+kH+k6cKdnKs9rt04H9W5IhPo5WYcE4GZIpKO8xHeFBF5u24hETkfGAp8BDzsZhsZQFytz7GALseslBcryMthSOH37On2A3z8/K0Op0XEZiPb3h3/Ip3L1RRXE9fHwCUiMtaFst3cDcIY86AxJtYYkwBcD3xrjLm5dhkRGQG8DFyOs+fTVUT+5EYz64D+ItJHRPxq2vnU3ViVUp5j95K38JdKuo73riWeGpLv34vQ8iyrw/B4riauuTh7JwtrEkhjJgP7WhRV/YKAa40x+40xDuBW4FDdQiIyH1gFDBSRDBG5A8AYUwXcDXyJc+Ti+8aY7W0Qp1KqnQTvXshhiaF/yvlWh9IqyoJjiKo+bnUYHs+ld1zGmFIRuQbnwIvvReQ54EVjzP7a5UTkXmA68GhzAzLGLAOW1XN8ZZ3PlTh7YHXLNfiG1hjzGc71E5VSXu744T0kVWxhVfxd9LZ1jIG+JjyOiJwi8vNOEh7hvSMk25rLd9sYsxbne6504FfA7prJwu+JyL9EZA/wf8ABmpjvpZRSLZW+7A0Ael/olYvz1MsvMh6AnIy2eGjVcbj1a4oxZhMwHPgfYCswDOeQ9CuBvsBi4EJjTGErx6mUUmcYh4Oe6Z+ww3cIMX0HWR1Oqwnt4VyuKi9LtzdpjNuL7NZM3n0ReLFmMd14wI5zcu+pVo5PKaXOsXfzSgY4jrBm8ENWh9KqusY453KVZ6dbG4iHa9Hq8DVznXJbKRallHLJyZWvUW58GXzRj6wOpVV16R5LhfHRuVxNaFbiEpEBOPfh6o5zEm82sM0Yo4tsKaXaVFlpMYNzvmRb2PmMinB79o1HE5udE/Zu+Opcrka5nLhEZDDOlSSuBXqcPlzzX1NT5jjwPjDXGFPvYrlKKdUS27+dzyiK8R/jnftuNSXftyehZbo2QmOaTFwikgj8DecAjFLgPzjnSe3H+ZhQgK5AP2Ac8GPgHhFZCPymZv1CpZRqFb5b55NFFIMnNGdJVM9XGhxN/MmVTRd0wbYViyhf8xokzWTUJR1n9KUrPa4dOEcQ3gYsNMYUN1ZYRIJxrg5/b821AS2MUSmlADiRcYChpRtYE3c7PX1avIG7R6oOi6PbyTyKi4sIDg5pdj0Htq1mwNe3Yaca+5pv2RocQfIFV7ZipNZxZTj8dcaY0caYt5pKWgDGmGJjzBvGmFHArJaHqJRSTge+eQWbGHpP+bHVobQZ367OuVwnMvY3UbJxhYv/QIkEkDNnI0elB37L/9JhdlduMnEZYz5pbuUtuVYppWozDgexhxay3S+ZmL5DrA6nzQT3SAAgP7P5ievo/m0ML1vHzvhb6BHTl4zBdzCwag/7N69opSit1THWSVFKdXi7131FrMmkaPD1VofSprrG9AegtAVzuQ4ve51qIyRO/wkAA6fcRqWxk7Pug2bVV1ZazKbHL2XzY9MozK9vn9721eqJS0RuFpFvW7tepVTnVrjqDYpNAEOndYyV4BsS2TOBKmPDcepws+vodmwpu/2G0D3GuRJHRFQPdgUMo1fW0mbVt+nDxxlRspLhZevY/q/Hmh1Xa2mLHlc8ru/dpZRSTSopymPIqSVs6zKF4NBwq8NpUzYfX7JtkfgWZjTr+uysw/Sr3k9BzNk/hotjzifecYScLPcnN8ce/ICdvkPY5p9CTMbiZsXVmvRRoVLK4+34+i2CpJzgcbdZHUq7yPPtSXAz53Klr/k3AFEpl5x1vOvQqQAc2viVW/UdObCTOMdR8vteRlHvqcQ5jpKVcbBZsbUWl8aTiog7c7E69q9DSql2F7TjPQ5LNENSL7I6lHZRHBRNbN6GZl1rDq2kgGD6Jk8463jf5AkUfxJA1YEVwB0u15e16UvigOiRF1NWnA97/0HmtmX0jO3TrPhag6s9rgScCanYha/KVo9SKdVppe/aSFLlNjL7XoN0kH23mlIdGkc3k0tZWZnb10blbyU9YDA2u/2s4z6+fhzy70943g636rNlrCKXcOL6D6d30jgqjJ2Kw5vcjqs1ufq34CCw3hiT3NQX8FwbxquU6mSylrxAhbHTv2aEXGdg7xqPXQzZx9x7JFdcmEd81SGKu6XUe76gyxDiKw9QVVnhcp3BBQfI8u+L2GwEBASQaY/G95S1y9K6mrg2ACNdLGuaGYtSSp2lrKSIpOx/szXsArr2iLU6nHYT1D0BgFPH3NtQ8uCWFdjFENR3bL3nfWNHECgVHN7tWo+ptLyK6KoMKrv0O3PsZGACUWXpbsXV2lxNXJuASBFJcKHsIWB5cwNSSqnTtn71BmEUEzCu466UUZ8uvRIBKDmR7tZ1BftWAxA/9Px6z3cbOA6A7L3rXKrv0OGDhEkJvj0GnDlWGtGPaEcWVeWlbsXWmlxKXMaYvxpjbMaYdBfKvm2MmdziyJRSnV7o9rc5ItEkjb/U6lDaVbeYvjiM4Dh1yK3rAo5vIkN6EdGtV73nYxOTKTYBcMy1HldO+jYAwuP+u1KJPaofPuLgRIZ166d3jjedzSAifUXkVRH5sOZzsIi8ISIvi8hNVsenVEd3cPtaBlXu4Gji9Z1mUMZpPv6B5EoX7AXuzeXqVbqH48EDGzxvs9s54ptASL5r76hKMncB0L1P8pljYTVLUmUfbdlaii3hEX8bRCRARNaKyGYR2S4if2xBXa+JyAkR2VbPuYtFZLeI7BOR3zZWjzHmgDGm9pjRq4APjTF3AjObG59SyjXZy16k3Pgy8OLOMyijtlzfngSXuj6XqzAvl17mBBVRSY2WKwjpS48K13py9ty9lOJmBhEyAAAgAElEQVSPX5f/vl+M6NUXgNJs6+ZyeUTiAsqBKcaY4UAKcLGIjKtdQES6i0honWP9ONc84OK6B0XEDjwPXAIkATeISJKIJIvI4jpf3eupNxY4PeW82s3vTynlhtKiAgZnf8aW8El0ieppdTiWKAnsRZfKLJfLH93jnPcVEDu80XKOqIFEkUdeTtN1hxYd5IRfHNTq8Ub26oPDCCaveSt7tAaPSFzGqajmo2/NV93RiRcCn4hIAICI3Ak8U09dy4H6VoFMBfbV9KQqgPeAy40xW40xM+p8najn+gycyQs85M9NqY5q21evEUopQeM716CM2ipD4+juyKGy0rWpsfnpaQD0Gji60XKBMc73Vcf2bW68/WoHPasyKA49e6KxX0AguRKBT9FRl+JqCx7zA1hE7CKSBpwAvjbGrKl93hjzAfAF8F7NO6bZwHVuNBHDf3tM4ExEMY3EEykiLwIjRORBYCFwtYi8ACxq4JrLROSl/Px8N8JSStVmjCF8+zuk2+JIGjvd6nAsY+sSh69Uk33MxQEax7dTQDA9Yvo2WqxHorNHVnDknLcpZzl8PIdYsiFqwDnncn26E1TSvCWpWoPHJC5jTLUxJgVnryZVRIbWU+ZxoAx4AZhZq5fmCqmv2UbiyTXG3GWMSawZVVlsjLndGPNTY8w7DVyzyBgzJzxcV71Sqrl2bVjGgOo9nBhwY6cblFFbYLcEAE66uC9XeMFuMnz7Nvln1iO2n3Nk4YmdjZY7fnAHNjEERQ8+51yxX3dCKnNciqsttOhvhTj1FhG/1grIGJMHLKP+91TnA0OBj4CH3aw6A4ir9TkWsO5XBqVUvYq+e44iE8iQS39qdSiWioh2vsIvPt70IAhHdTVxFQcpCD+3d1SX2Gwc9e1NcEHjk5uLa0YURvY+d9POyqAowh2nmmyrrbT015muOJeDOq8llYhINxGJqPn/QGAasKtOmRHAy8DlwO1AVxH5kxvNrAP6i0ifmkR7PfBpS+JWSrWunGOHGF6wlG09LiM4rIvV4VgqKsY5Cbn6ZNOPCrMO7yFYyrD1POdBVb3yQ/rSszy98ULZziHzoTGDzjklwd2JoIiyMmsmIbdGP7y+R3Du6gUsFZEtOBPM18aYupu+BAHXGmP2G2McwK04V+k4OxiR+cAqYKCIZIjIHQDGmCrgbuBLYCfwvjFmeyvErpRqJfs+fwYfHMRM/7nVoVguICiUXMKxFzS9f9bxvc4RheF9RrhUd3XXAXTjFPmnchssE1iwnxO27uAXdM45e2gPAE4et2aAhkvbmrQ1Y8wWoNE/cWPMyjqfK3H2wOqWu6GROj4DPmtmmEqpNlRRVkr/Ix+wJWgsKf1c6zl0dCd9ehDowiCIsowtOIwQN9C1JWUDY5LgAGTuSyN8zNR6y0SWHeZkUDz1zQ3yjXBOUSg6mQnx9c1Kalud982nUsqjbPvqdSLJh3F3WR2KxygM6EVERWaT5fxzd3LU1ougENcGhnVLGAZA/uH6RxaWVVQR5zhKWXhivecDahJX6Skv6HGJyAV1Dp3+UxomIlW1T9TMp1JKqaYZQ9iWV0mXWIadd7nV0XiM8tA4ehaupLqqCrtPwz+uu5XsIzso8azRZ43pET+QcuOLo4GRhUcPHyBRypBu9Q/2COnqnElUkXfcxRZbl7uPCpdR/xDyf9T6f6kpY6+nnFJKnWPfhm/oV7WPlYN+T4JdHwSdZovsh19mFZlH99Mrvv41CEuK8olxZJIROcPleu0+PhzyiSUwv/6h9rmHtpMIhNQzMAMgvHs0AI6i+tZqaHvuJq66q75H4Byafj/OPbuUUsptBd89T74JJvnSOVaH4lFCYgbBNshJ395g4srYk8YAMfjHDHOr7lNBfehVVP+jwvIsZ0+s9uK6tQUHh1JkArEVe0HiMsZ8V/uziETW/G9a3XNKKeWKzCP7GVbwHet6Xc/4sAirw/EoPeKdC+YWZ+5psEzeQecWJd0TXRtReFpl1wH0LFhKcWE+waFnvxuzndxHEYGERtX/8FFEyJcwbGXWzOXSPrlSylIH/v1/CIY+l9xndSgeJ7JnnHOVi9yGJws7jm+nxPgT3efcFS4a499rMDYxHNu/5ZxzoUUHyfKJBWl4tlOxPQzf8jy32mwtmriUUpbJzztJcua/2Bw2iZ4NPArrzMRmI8snmsDC9AbLhObtJsM3HpvdvWEFXWtGFp46dO7jwu7lh8kL7nPO8drKfMMJqLJmXVZNXEopy+xY9DRhUkLY1F9aHYrHyg+Kp2t5/VuIGIeD6IoD5IU2vdRTXdF9k6gyNqqOn7VIESVFefQkh6qIxhfrrfQLJ7jaOxPXSaAPsLKpgkopVVtFeRmJ+99km18K/VLqzrRRp1VG9KGX43i9yyvlZh2hC4U4up+7nmBTfP0COGaPJiDv7N2Qj+11Pjr069V4nY6ALoQ43FnnvPW0KHHV7KN1yBhT3loBKaU6h61fvEx3TlI9/l6rQ/Fovt0H4iMOMg/sOOfcsT3rAQjt3fjmkQ3JDexDVGn6WcdOHXImrqi+jY9SNAFdCJdiyiva/8e/PipUSrU746gmavNL7LclkHzBlVaH49FC41MAyDu48ZxzxYecs5BiB6c2q+7yLv2JdmRSVlpy5lh11g4qjA/RfZIavzjYOai84FR2s9puCU1cSql2t3XpB8Q7DpM97C5sOuG4Ub0HDqfC2Ck/eu7ov4DjaRyRaMK7dmtW3X49B+MjDo7u2/rfOvP3kWGPxce38d2q7MFdASjWxKWU6uiMMfisfoYsujHq0tlWh+Px/P0DybDHEXhq1znnYkp2cjzU/fdbp0X2HwPAyb2rzxzrUXqAU8GND8wA8A1x9rjK8jVxKaU6uG2rvySpcjtHBt2Or5+/1eF4hZOhA+hRevbyTCeOHqQ7J6nq5dqK8PWJ6zeMAoIxGesAyD2eQS+yqeje9DuzgDBnL6+8sOGtUdqKS4lLRAaJyI0iMklEQpsoq2sUKqUaZJY9zknCSL7sHqtD8RrV3ZPpSS45x9LPHDu6zbmOeUT/sc2u12a3cyhgEFF5zkeFh7f+x1nngPFNXhsU7kxclUU5zW6/uVztca0Anga+BnJFZKmI3Cwi9f26tEBE3J9UoJTq8HatW8Kw8g3s63c7AcFhVofjNUIHOqcLZGz+5syxir3LKDH+9E2e2KK6i7uNJL76EPk5mZTtX0GlsdMneUKT14VEOBOXo/hki9pvDlcTVxbwAtAb50K73wG/A9JF5B4R8a1V9mXg9VaNUinVIVR8+xinCGXoFTrh2B39hk2g2ARQvm/FmWMxuavYEzQCP/+AFtXddeTl2MWwZ/kCoo8vY3dAMgFBjT5YAyA0IpIqY4NSz01cdwK3AvuBPwCFwBzgDmAasENELq4puxVIaN0wlVLebv+m7xhWupZdCT8iKEQX03WHn58fBwKH0vOU811Uxv7txJpMSuNaPnG7//CJHJMeDNv8J+IdGRQnurY9it1uo0CCsZW1/3qFLiUuY8wqIBG4CTiBcxuT74BFQAoQBvxbRDYAe4GFbRKtUsprlXzzF/IIYciVv7Y6FK9UnDCNeEcGB7etImPZaziM0PeCWS2uV2w2Mkf9Gn+p5LDEMHyG6ztQF0sItoqCFsfgLpe3NTHGVOHce+sjgJr3WOOAJCCqplgPYCiwpHXDVEp5s4NbVpBcvJoVcT/hvPCuVofjlQZOm03FzicoWvQ7BpTvY2vgaIbH9muVukfNmMOxoRcSGdnDpceEp5XaQvCpbP9ln9zdSPIMY8we4JxNYkTkBpw9s49bEJdSqgPJ++LPFBBM8tW/sToUr9Ulqgdrev+IsUdeo8z4EnzJI61af3SC+6vzl/uE4FdV2KpxuKLZiashxpj5IvKv1q5XKeWddm9awYiS71kT/xPGRkQ2fYFqUOrt/2DLf84nrEcf+g1yb+PItlDpG0pYiQdOQBaRqe5WaoypqLl2WnOCUkp1HGVfPkI+IQy5SntbLSU2G8MuvIoED0haANW+YQQ5itu9XVcGZ3whIt+KyAxXJheLiK+IXCki3wGftTxEpZS32rHqc4aXrWN34mxCwrW31dE4/MMIMcUYY9q1XVceFY4AngQ+BXJE5GtgLc6h8ScBAboC/XEO1pgCdAG+wjniUCnVCRmHA9u3j5JNF5KvesDqcFRbCIggUCooKSslKDCo3ZptMnEZY7YB00VkPPA/wOXADUDdFCtAAc6h8C8YY9a1cqxKKS+yffmHDK3cweqk3zMu2PWRasp72IPCASjMO+lZieu0mrlcq2oeF47COQy+G84Elg1sAzYZYxxtEahSynsYRzVB//krGdKTEZfrRpEdlU9QFwCK83OhV2z7tevuBcaYapyPCte2fjhKqY5gyxevM7z6AKuGP0ZsC5ckUp7LN9jZ4yorat9ln1xKXCIyCBgJHAM2GGMaHLgvIvaa5KaU6oQqK8qJWvd39tsSGHPZnVaHo9qQf6hzMnlZwal2bbclq8PfoqvDK6XqWv/xM8SYLAonPIiPT6tPFVUeJDDEmbgqij0zcdW3OvyD6OrwSqlaCvNP0n/Hc+z2TWL4lOusDke1seCaKQ5VJe270K6uDq+UajXb33+EKPKwXfIXxKYbrHd0ITXrTjpK2zdxudSPN8asEpFE4DLgCpyrw3erOZ0BBOBcHT4NGAS81gaxtisR6Qv8Hgg3xlwjIsHAP4EKYJkx5h1LA1TKw5w4vIsRGe+yLvwixoycbHU4qh34BIZRbQQpa98V4l3+lcgYU2WM+cgYc6sxpicwGLgdmI9ze5PXcQ7e8MHN1eFFJK7mvdlOEdkuIj935/o6db0mIidEZFs95y4Wkd0isk9EfttYPcaYA8aYO2odugr40BhzJzCzufEp1VEd/eA3OBBirn7M6lBUexGhWIKR8vx2bdZTVoevAn5ljNkoIqHABhH52hizo1a93YHS2iMaRaSfMWZfnbrmAc8Bb9aJyw48D1yEs5e4TkQ+BezAX+vUMdsYc6LOsVicj0EBdNSkUrVsX/U5IwqX8X3vOUyIb52tNpR3KLYFY69s3xXiW/0htDFmPs7E5c41mcaYjTX/XwjsBGLqFLsQ+EREAgBE5E7gmXrqWo5zKaq6UoF9NT2pCuA94HJjzFZjzIw6X3WTFjiT3ekZdvrwXqkaVZWV+H39e44TyYjrH7I6HNXOymwh+Hl74oL/rg7fHCKSgHN9xDV16vwA+AJ4T0RuAmYD7gxbigGO1PqcwbnJsXYckSLyIjBCRB7EuZTV1SLyAs5Ho/Vdc5mIvJSf377dZqWstPZfT9LfsZ+s1AcJ1KWdOp1yn1D823lPLo+aZCEiIcC/gPuMMee87TPGPC4i7+Ecmp9ojHFn602p51iDSxobY3KBuntY395YA8aYRcCi0aNH66xL1SnkZB1h6M6n2B6QwrCL72j6AtXhVPqGElR+qF3b9JhHXjVzwf4FvGOMWdhAmfOBocBHwMNuNpEBxNX6HItzMIlSqpnS5/+KAMoJvfppHf7eSVX7tf+eXB7xN01EBHgV2GmMebKBMiNwTm6+HGfPp6uI/MmNZtYB/UWkj4j4Adfj3KpFKdUMe9Z8zuj8L1kXczO9B+gORp2Vwz+cEEqoqm6/9dU9InEBE4FbgCkiklbzdWmdMkHAtcaY/TUr0N8KnNM/FZH5wCpgoIhkiMgd4BzOD9wNfIlz8Mf7xpjtbfctKdVxVVWU4f/VAxyjO8NvdOf3R9XRSEA4oVJKQUlZu7XpEe+4jDErqP8dVO0yK+t8rsTZA6tb7oZG6vgM3ZVZqRbbuODPpFYfZv2EF4gOCbM6HGUhW1AEAEV5uXQNbZ89uTylx6WU8hKZh3YzdN9cNgZNZNRFDf6eqDoJn5rEVVzQflubaOJSSrnMGEPWe/cB0GvWUzhfT6vOzK9ma5PSgtx2a1MTl1LKZRu/epcRpd+ztd9d9IrX3YsUBJzZ2kR7XEopD1OQl03cqj9w0J7AqFm/szoc5SECa7Y2qSxqvz25NHEppVyy94176GryqJjxHD5+AVaHozxEcJgzcVW3455cmriUUk3atvR9Rp36nNUxtzJwxPlWh6M8SEDNOy7TjntyaeJSSjWqIC+bHt/9hgO2eEbf8herw1EeRvyCqcLerlubaOJSSjVq97x76GLyqLzseQIC22eejvIiIhRKCHZNXEopT7D52/cYk/c5a2P1EaFqWKktGJ/K9tsFWROXUqpeOVlHiF3+AAdtCYzSR4SqEaX2MPw1cSmlrGQc1RybdzvBpgS59lX8A/QRoWpYhU8oAdXu7DLVMpq4lFLnWLvgrwwrW8fmpPtJGDza6nCUh6v0CyPIoYlLKWWRA9vWMGLX/5EWOJ7Ua++3OhzlBRz+YYSYIoxpcG/eVqWJSyl1RllJEbaFd1AgIfS+/VXdHFK5xAREEEYxxeVV7dKe/q1USp2x+bV7SHAc4djkp+jaPcbqcJSXkIBw/KSawsL2GaChiUsBcOzADjZ+8QZl5eVWh6IssmHxy4zNWciqHjcw7MIrrQ5HeRGfoC4AFOXntEt7mrgUABX/+gkjV99L9l+T+eyNx8nIab/lW5T1Du9cx+B1v2eH71BG3/G01eEoL+MT4kxcpe20J5cmLgXH0kgo3sIn5gKq/SO49OCf4dlRvPXcw6zcdbTdXrgqaxTnn8T2wY8oliCibp+Pr5+/1SEpL+Nfs7VJWWH77MmliUvBmrmU2wJ51u9OEn67hpwr3sUW2pNbcp6iz/zzePaxB3jrP7spLKu0OlLVyoyjmv0v30yP6uNkTX+R7tG9rQ5JeaGgmhXiK4vb50mNJq7Origbtn3IhohLKJJgECEq5YdE/2oFFTcuxC+yD/eWv8T0b6bzz7/8kkcXrmfv8UKro1atZMO7jzCsaCWr+91H8oSLrQ5HeamgmhXiq4rbZ08un3ZpRXmuDa9DdQUrI6+G2vMHRfAbMJWoAVMhfQVBX/6Z32S+Sc7mT3hpw6Xs7T2LWRMHM21wd3zs+vuPN9r83UeM2Pssa0MmM/HGP1gdjvJiwRFRQPvtyaU/cTqzqgpY9yr0m0aOfyOPiBLOI/Qnn8PsLwnrM5Lf+c7n/zJvYcv8h7j4b//m+aX7yC3S0YjeJH33ZhK+/RlH7HEMuWseNv3lQ7WALTAcACnTxKXa2o5PoCgLxt7lWvne4/C77WP48RLC+o/nAd8FfFx5F+Xf/JXpf13ELxekkXZERyN6ulM5Wdjfm4VD7ATe+gHBoRFWh6S8nd2XEgKwa+JSbW7NixDZDxKnundd7GhsN30Ac5YRMuBCfun7ISv972XAjqe57fkvuG7uKrYfa7+9eZTrystLyZx7Dd0dOWTPeJ0e8YOsDkl1EPm2CPwrdFShaksZ6+Hoekj9CTR3WZ/oEXDDu3DXCgIGTeMuWci64F9wcdZcbn32M/7w8VbySipaN27VbI7qarY+dzNJlVvZNuavDBg9zeqQVAdSYO9KUKXO41Jtac2L4B8GKTe0vK6eyXDdm/A/q/EdfCm3m4/5PvA+4tf/lav//hFvrz5EtUPnglnJOBysfvGnjC78htV9fsaoGXdaHZLqYEr8uhJa1T6jCjVxdUYFmbD9IxhxM/iHtl693QfDNa8iP1uL39Ar+LHv53xu7qZi8QPc9swnrE9vn9/G1LlWvvUwE7IXsLb7dYy95U9Wh6M6oDL/SMId+o5LtZX1r4GjGsb8uG3q7zYArpqL3L0e35TruM33G17J+zE7XpnDH9/+khMFZW3TrqrX9x8+xXkHn2FT2BRG/+RFXfFdtQlHYDfCTSGOqrZfqED/Bnc2VeXOxDXgBxCZ2LZtRSYiVzyP7d4N2FJu5CbfpTy49waW/eNG3v1iORVVjrZtX7Hx3y8xbusjbA0YzdCfvYvNbrc6JNVB+YT1wCaGk9nH2rwtTVydzbaFUJLj+hD41tAlAd8rnsH+8zTKht3ClbKc61ZdzpLHrmHN+nXtF0cns/GLNxi29jfs8E+m3z0f4+sfaHVIqgPz79ILgJMnMtq8LU1cnYkxsOYF6DYI+k5q//Yj4gi7+ml8f7mFzIG3MKXqP4xedBHfP3ENxw7uav94OrANX77N0FW/YJ/fQOJ/9gmBwa34LlOpeoRGOfdvK8zWxKVa05E1kLkZxv4ERKyLIyyauBufgfs2s7X3zaQULqfLvPNZ+8bvqKrQ918tlfbZywz//h4O+SYS87PFhIZ3tTok1QmER/cDoCrnQJu3pYmrM1n9AgSEw7BZVkcCgH9ENCl3PEfBj1exPXgsqQef59hjo9i39nOrQ/Na6xY+zbA197PbL4me93xJaM0ackq1ta7dYikx/phT6W3eliauziI/A3YugpE/Ar9gq6M5S8+4REbdv4j1E+did1TQ77Pr2fzMLIpPZlodmlf5/u0/MmbLQ2wLHEXCzz/XnpZqVza7jV2+SWSXtkNbbd+E8gjrXgUMjPHMiaciwuiLrif0V+v5rsetDM79mupnRrHj06fAoaMPG2Mc1ax54SdM2PckaSHnM+i+xQSHhFkdluqESq//kLhr/9bm7Wji6gwqS2HDPBh4KXSJtzqaRoWFhnPhT59h79VfcsAnkaSND3PgbxPI2aujD+tTVlrMpievZOzx91jd7RqG3fcxfgE6elBZY2K/KFLi2n7RZk1cncHWD6D0JIz7qdWRuGzIsDEM+e0yvhn8J8LKjtHl7YvY8dr/UF2qi/eeln08g/1PXsTIou/4PvEXjP3py9h8dIs91fFp4urojIHVL0KPoRA/0epo3OLrY2farHsonbOKpSE/ZNChd8l7PIUDy95yfl+d2MHNy6l+4UISK/aweeyTTLjlEV0RQ3Ua+je9o0tfASe2Wz8EvgXiYmKY+uu3WTnpPbKJoO+yu9n9j+mcPLLT6tAssf6jp4lZeCUGIeOqjxl+yR1Wh6RUu9LE1dGteRECu0LytVZH0iIiwvmTLybugdV8Ff9LYgq3Evzq+Wx687dUlbfDMCYPUFRUwKqnbmb05ofY5Z+M/aff0W/4eVaHpVS708TVkZ06BLs/g1G3gW/HeGEfHOjP9NsfJvu2FWwKnMCIAy+Q9bdR7Pp+kdWhtakDad9x6h9jGZ+3iHWxP2LIA9/QvUeM1WEpZQlNXB3ZupcBabtV4C3Up08/xj7wCevOewVjHAz66mY2PHkVWYc61tJRFeXlrH7tfnp/dAV+ppwdF73FmB8/i10HYahOTBNXR1VRDBvfhKSZEN4xfzMXEcZMu5ao+zfwfdydDM1fTrfXxrHl6WvI3b/R6vBabPemFRx4/DzGHX6JTeFT8btnDUkTZ1odllKW08TVUW1+D8ry23cVeIsEBgUz4Y4nOHXnOlZ2u56+J/9D5FuT2fvUpeTv/o/V4bmt+NRxNjx/K/0/nkH36iy2jn+KMb/8kC6R3awOTSmPoImrIzIG1syFXikQN9bqaNpNz9g+XHD3i5yas4HPu80m8tRmwufP4MiTkyje/oXHD6Gvrqpkw/uPUf30CIaf+JRV3a/F975NJP/gdqtDU8qj6IPyjujAUsjZDVe86LVD4FsiLiaWuJ/9H/uP/pYlHz/NxBPzCf5gFpmBA7Bf8Eu6j70ObJ6zoWJVeQlbP3uJyK0vMcpxlC1+Kfj88O9MHJ5qdWhKeSRNXB3RmrkQ3A2GXmV1JJZKjOlB4s/+wo4jv2Tpv19ifOab9P3yLrKW/D/yRvyMAdN/jM3X37L4ik5mseffT5Ow/x1GkM9eWyIbxz3LiOk362RipRqhiaujyd0Pe76ECx8AH+t+KHuSpLgoku76HSfyfs7iL94kcddLDF73O06s/wcH+9xA7zEz6DUwtV16Yaaqgr3rviJv3QKSc79gpFSw3m8M+8ffzagLZmKza8JSqimauDqatS+DzQdGz7Y6Eo/TPSKYGdf/lIrKOaxaupDQdc8w9sBzcOA5CiSE4xGjCBw4mejhF2HrkQSt1OupKi/hwNrFlG/5mN7ZyxlAIaXGjy1dphEx9ReMTtZHgkq5QxNXR1JeCJvehiFXQmhPq6PxWH6+dsZPvxamX8vRQ/vZu/ZzzMHvSDyZRuzq72D1IxTawsnsMprq3hPpmjyN7gnJTT++MwZH4QnyM3Zw/OA2So7twufUXvqVpDGAcvJNMFuCx2FPupzkC69gbGh4+3zDSnUwmrg6krR3oaKwUwyBby0x8YnExN8N3E1eSQVfbErj5LYlROWsYWhOGtG5S2DTo1QaOxXiR6X4U2Xzp8ruT4X4UyF+VOCHb3UJPSuPEEoJXYAuQJnx5ag9mi1dLsY+5DIGjLuU80M8axNPpbyRJq6OwuFwDsqIHQOxo6yOxitFBPlx8cRUmOh8dFdaXsW2vdso3PktjpMHqSovpbqiBFNZiq2qjAAqCZAK/CmjzBZIWsR0ysL7QmQ/ohKG0q//IBID/Um0+PtSqqPRxNVR7F8CJ/fD5N9ZHUmHEejvw9ChKTA0xepQlFK16BCmjmL1CxDSEwbrkkBKqY5NE1dHkL3H2eMa82Pw8bM6GqWUalOauDqCtXPB7ufcvkQppTo4TVzerjQP0uY7N4oM0UVYlVIdnyYub5f2DlQWQ+ocqyNRSql2oYnLmzmqnUPge4+HaB35ppTqHDRxebM9X0LeIZ1wrJTqVDRxebM1L0JYDAyaYXUkSinVbsR4+OZ63khE8oG99ZwKB/JdOBYF5LRBaE2pL5b2qsfVa5oq19B5d4570j2B1rkvnnhPGjvn6ffF0/+teOs9iTfGND3KzBijX638Bbzk6vEGjq33pLjbox5Xr2mqnDt/9t5wT1rrvnjiPfHm++Lp/1Y6+j3RR4VtY5Ebxxsqa4XWiqU59bh6TVPl3Pmzb+i4J90TaJ14PPGeNHbO0++Lp/9b6dD3RB8VeiARWW+MGW11HOq/9J54Jr0vnqc97on2uDzTS1YHoM6h98Qz6X3xPG1+T7THpZRSyqtoj/184tAAAAOISURBVEsppZRX0cSllFLKq2jiUkop5VU0cXk4EQkWkTdE5GURucnqeJSTiPQVkVdF5EOrY1FOInJFzb+TT0RkutXxKCcRGSwiL4rIhyLy09aoUxOXBUTkNRE5ISLb6hy/WER2i8g+EfltzeGrgA+NMXcCur1xG3LnvhhjDhhj7rAm0s7DzXvycc2/k9uAWRaE22m4eV92GmPuAq4DWmWYvCYua8wDLq59QETswPPAJUAScIOIJAGxwJGaYtXtGGNnNA/X74tqH/Nw/578oea8ajvzcOO+iMhMYAWwpDUa18RlAWPMcuBkncOpwL6a3+QrgPeAy4EMnMkL9H61KTfvi2oH7twTcfob8LkxZmN7x9qZuPtvxRjzqTFmAv+/vbtnjSKKozD+HFBTieBrES0sJKKVVtoJvqSxEaPEMljY6AfwC/gB0kQbWzFFBEFBsIiFiATBJlgLwc7CQkQUrkVCYpEF2Sx75zrPr1l2ip0Df2YOcxfuwEj+7vBG2B2TbD1ZwXphTQJLwPUkC3Rry5u+2HYuSQ4keQicSXK/TrTeGnSt3AMuATNJfNfP+A26Vi4kmU/yCHg5ihPtGsWPaCSyzbFSSvkOzI07jDYNmstXwJtjHYNmMg/MjzuMNg2ayzKwPMoT+cTVHWvAsb++HwW+VMqiLc6le5xJN41tLhZXd6wAJ5IcT7IHmAWeV84k59JFzqSbxjYXi6uCJE+Ad8BUkrUkt0spv4G7wCvgE7BYSlmtmbNvnEv3OJNuqj0XN9mVJDXFJy5JUlMsLklSUywuSVJTLC5JUlMsLklSUywuSVJTLC5JUlMsLqlHkpxPsprkq69nUassLqknkuwGngLPNg6drRhHGprFJfXHVeAwsADsw41p1SiLS+qPa8Ab4PjG9w8Vs0hDs7ik/pgGXgOXgbellG+V80hDsbikHkgyxfoy4XvgJvC4biJpeBaX1A/nNj73AvuBxYpZpB2xuKR+OMX6G2rvAA9KKT8q55GGtqt2AEljcRCYAE4DNypnkXbE4pL64xAwV0r5WTuItBMuFUr9cAT4WEp5UTuItFMWl/SfS3ILuAKcTDKZZDrJ0sZOGlJzXCqU/mNJJoAZYBa4CKwAn1lfMvxVM5s0rJRSameQJOmfuVQoSWqKxSVJaorFJUlqisUlSWqKxSVJaorFJUlqisUlSWqKxSVJaorFJUlqisUlSWrKH+ROejQIVOK0AAAAAElFTkSuQmCC\n",
"text/plain": [
"