{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Deblending with *Scarlet*\n",
"
Owner(s): **Fred Moolekamp** ([@fred3m](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@fred3m))\n",
"
Level: **Intermediate**\n",
"
Last Verified to Run: **2021-03-26**\n",
"
Verified Stack Release: **v21.0.0**\n",
"\n",
"The purpose of this tutorial is to familiarize you with the basics of using `scarlet` to model blended scenes, and how tweaking various objects and parameters affects the resulting model. This tutorial use `scarlet` as a stand-alone package, and does not depend on the LSST DM Science Pipelines. This notebook was developed from an early version of the [scarlet quickstart tutorial](https://pmelchior.github.io/scarlet/0-quickstart.html). A tutorial that is more focused on using `scarlet` in the context of the LSST DM Science Pipelines is available at [Deblending/LsstStackDeblender.ipynb](../Deblending/LsstStackDeblender.ipynb).\n",
"\n",
"### Learning Objectives:\n",
"\n",
"After working through this tutorial you should be able to: \n",
"1. Configure and run `scarlet` on a test list of objects;\n",
"2. Understand its various model assumptions and applied constraints.\n",
"\n",
"More documentation is available on in the `scarlet` [docs](https://pmelchior.github.io/scarlet/) and the core concepts of the code are discussed [here](https://pmelchior.github.io/scarlet/1-concepts.html).\n",
"\n",
"### Logistics\n",
"This notebook is intended to be runnable on `lsst-lsp-stable.ncsa.illinois.edu` from a local git clone of https://github.com/LSSTScienceCollaborations/StackClub.\n",
"\n",
"## Setup"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-23T20:41:31.149585Z",
"iopub.status.busy": "2021-04-23T20:41:31.145847Z",
"iopub.status.idle": "2021-04-23T20:41:33.141633Z",
"shell.execute_reply": "2021-04-23T20:41:33.142859Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"nb-kadrlica-r21-0-0\r\n",
"lsst_distrib 21.0.0+973e4c9e85 \tcurrent v21_0_0 setup\r\n"
]
}
],
"source": [
"# What version of the Stack are we using?\n",
"! echo $HOSTNAME\n",
"! eups list -s | grep lsst_distrib"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-23T20:41:33.157121Z",
"iopub.status.busy": "2021-04-23T20:41:33.155775Z",
"iopub.status.idle": "2021-04-23T20:41:34.055532Z",
"shell.execute_reply": "2021-04-23T20:41:34.056774Z"
}
},
"outputs": [],
"source": [
"# Import the necessary libraries\n",
"import os\n",
"\n",
"%matplotlib inline\n",
"import matplotlib\n",
"import matplotlib.pyplot as plt\n",
"# don't interpolate the pixels\n",
"matplotlib.rc('image', interpolation='none')\n",
"\n",
"import numpy as np\n",
"from astropy.visualization.lupton_rgb import AsinhMapping\n",
"\n",
"import scarlet\n",
"import scarlet.display\n",
"from astropy.visualization.lupton_rgb import AsinhMapping, LinearMapping"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-23T20:41:34.066039Z",
"iopub.status.busy": "2021-04-23T20:41:34.064580Z",
"iopub.status.idle": "2021-04-23T20:41:34.069246Z",
"shell.execute_reply": "2021-04-23T20:41:34.070380Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"scarlet: 1.0.1+unknown\n"
]
}
],
"source": [
"print(\"scarlet: %s\"%scarlet.__version__)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Display functions\n",
"\n",
"Below are several useful functions used throughout this tutorial to visualize the data and models."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-23T20:41:34.092582Z",
"iopub.status.busy": "2021-04-23T20:41:34.091161Z",
"iopub.status.idle": "2021-04-23T20:41:34.094061Z",
"shell.execute_reply": "2021-04-23T20:41:34.095182Z"
}
},
"outputs": [],
"source": [
"def show_psfs(psfs, filters, norm=None):\n",
" rows = int(np.ceil(len(psfs)/3))\n",
" columns = min(len(psfs), 3)\n",
" figsize = (45/columns, rows*5)\n",
" fig = plt.figure(figsize=figsize)\n",
" ax = [fig.add_subplot(rows, columns, n+1) for n in range(len(psfs))]\n",
" for n, psf in enumerate(psfs):\n",
" im = ax[n].imshow(psf, norm=norm)\n",
" ax[n].set_title(\"{0}-band PSF\".format(filters[n]))\n",
" plt.colorbar(im, ax=ax[n])\n",
" plt.show()\n",
"\n",
"def display_diff_kernels(psf_blend, diff_kernels):\n",
" model = psf_blend.get_model()\n",
" for b, component in enumerate(psf_blend.components):\n",
" fig = plt.figure(figsize=(15,2.5))\n",
" ax = [fig.add_subplot(1,4,n+1) for n in range(4)]\n",
" # Display the psf\n",
" ax[0].set_title(\"psf\")\n",
" _img = ax[0].imshow(psfs[b])\n",
" fig.colorbar(_img, ax=ax[0])\n",
" # Display the model\n",
" ax[1].set_title(\"modeled psf\")\n",
" _model = np.ma.array(model[b], mask=model[b]==0)\n",
" _img = ax[1].imshow(_model)\n",
" fig.colorbar(_img, ax=ax[1])\n",
" # Display the difference kernel\n",
" ax[2].set_title(\"difference kernel\")\n",
" _img = ax[2].imshow(np.ma.array(diff_kernels[b], mask=diff_kernels[b]==0))\n",
" fig.colorbar(_img, ax=ax[2])\n",
" # Display the residual\n",
" ax[3].set_title(\"residual\")\n",
" residual = psfs[b]-model[b]\n",
" vabs = np.max(np.abs(residual))\n",
" _img = ax[3].imshow(residual, vmin=-vabs, vmax=vabs, cmap='seismic')\n",
" fig.colorbar(_img, ax=ax[3])\n",
" plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Load and Display the data\n",
"\n",
"The `file_path` points to a directory with 147 HSC blends from the COSMOS field detected by the LSST pipeline. Changing `idx` below will select a different blend."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"execution": {
"iopub.execute_input": "2021-04-23T20:41:34.111813Z",
"iopub.status.busy": "2021-04-23T20:41:34.110495Z",
"iopub.status.idle": "2021-04-23T20:41:35.695150Z",
"shell.execute_reply": "2021-04-23T20:41:35.696421Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"--2021-04-23 20:41:34-- https://github.com/pmelchior/scarlet/raw/master/data/hsc_cosmos_35.npz\r\n",
"Resolving github.com (github.com)... 140.82.112.4\r\n",
"Connecting to github.com (github.com)|140.82.112.4|:443... connected.\r\n",
"HTTP request sent, awaiting response... 302 Found\r\n",
"Location: https://raw.githubusercontent.com/pmelchior/scarlet/master/data/hsc_cosmos_35.npz [following]\r\n",
"--2021-04-23 20:41:35-- https://raw.githubusercontent.com/pmelchior/scarlet/master/data/hsc_cosmos_35.npz\r\n",
"Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 185.199.111.133, 185.199.109.133, 185.199.108.133, ...\r\n",
"Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.199.111.133|:443... connected.\r\n",
"HTTP request sent, awaiting response... 200 OK\r\n",
"Length: 242498 (237K) [application/octet-stream]\r\n",
"Saving to: ‘hsc_cosmos_35.npz’\r\n",
"\r\n",
"100%[======================================>] 242,498 --.-K/s in 0.02s \r\n",
"\r\n",
"2021-04-23 20:41:35 (15.4 MB/s) - ‘hsc_cosmos_35.npz’ saved [242498/242498]\r\n",
"\r\n",
"Background RMS: [0.35903206 0.85931915 1.004924 1.2481077 1.1298395 ]\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAeMAAAJOCAYAAACEHW+VAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABLVElEQVR4nO3de7RlaV3e++c3L2utfatbV9E03QiIKJBEQfsQlIgcjB68whjxGBAVM4zkdk40mhjiwCgmRj1x4CWGE/ESOh4VCInCIJhIWtDEGBCEINAQuXTb0E3f6rKv6zbne/7Yqzubctd+n6qu2m931fczRo3al7nn+653vnP+1txrvc+OlJIAAEA5VekOAABwraMYAwBQGMUYAIDCKMYAABRGMQYAoDCKMQAAhVGMgSsoIp4XEZ8q3Q8Aj2wUYzxqRMTtEfGXS/cDV15EPD0i3hMRZxb//nNEPH3P9384ImYRsbnn3+eW7DPwcFCMATwS3SXpmySdkHRS0lskvf68bd6QUlrd8+8Th91J4HKhGONRKSK+IyJ+PyJ+KiLORsQnIuLLFl+/MyLujYiX7dn+6yLifRGxvvj+D5+3v2+PiDsi4oGI+MG9d+ERUUXEKyLi44vvvzEiTlxkf79v0ae7I+Kv7fn610bEhyNiIyI+HRF/f8/3XhgR71/0+eMR8YJMGyci4l9HxF2Lu8nf3PO974qIj0XE6Yh4S0Q8bvH1WIzhvYt2/jgi/vzie6+LiNdExG8t7jx/PyIeGxE/vdj/RyLimXvaeFpEvHNxPD4UEd/oPM79pJTOppRuT7sRgSGpk/R57nifNy7fERH/NSJ+ctHvT0bE1+z5/jsj4sci4t2LMXjzg8c3In7uvLvv+flzB7gsUkr849+j4p+k2yX95cXH3yFpLumvSaol/VNJfyrpX0oaSvpqSRuSVhfbP0/SX9DuE9AvlHSPpBctvvd0SZuS/pKkgaSflDTb09Z3S/rvkm5a7PvnJf36nn59QNK3XKDPz1v080cktZK+VtK2pOOL798t6csXHx+X9MWLj58l6Zykr1r0+UZJT82Mz3+Q9IbFflpJX7H4+vMl3S/pixf9/xeSfm/xvf9D0nslHdNu0XuapBsW33vd4ue+RNJI0u9I+qSkb98z5u9YbNtK+pikH1iM4fMX4/8FBz1O45ifXYxfL+mVe77+w4vxOS3pQ5L+1gH7+I7F8fyuRb//lnbvvGPx/XdK+rSkPy9pRdK/k/T/7bOfZ0i6T9IzS58L/Lv6/hXvAP/45/7Tny3Gf7Lne39BUpJ0/Z6vPSDpGRfY109L+qnFx/9Yn11clyVN97R1m6Sv3PP9GxYX98bo8/Mk7ezdVtK9kp69+PhPJf0NSUfO+7mff7B/5tjcsChYx/f53i9J+n/2fL666P8TF0Xzf0p6tqTqvJ97naRf2PP5/y3ptvPG/Ozi4y+X9Jm9+5D065J++KDHaT62FUl/W9LX7fna0yU9blFcv0y7xf4lF/j575D0sfOOb5L02MXn75T04+fteyqp3vO1U4v59+LS5wH/rs5//Joaj2b37Pl4R5JSSud/bVWSIuIvRsQ7IuK+iDgn6W9q97VIafeifueDP5RS2tZuIX/QEyT9xuLXr2e1W5w7Sdeb/XwgpTTf8/n2g/2S9Fe0e7d8R0T8bkR86eLrj5f0cXP/D25/OqV0Zp/vPU7SHQ9+klLa1O7juzGl9DuSfk67v1G4NyJeGxFH9vzs+eO57/gu2rgzpdTv+f4d2r2jly78OLNSSluS/pWkfxMRj1l87cMppbtSSl1K6b9J+hntvsZ8IZ/Zs7/txYere75/556P79Dunf5JSYqIVtKbJP1aSun8162By4JijGvFr2n3TUCPTykd1e7FPRbfu1u7v4KWJEXEkqTr9vzsnZK+JqV0bM+/UUrp0w+3UymlP0wpvVDSYyT9pqQ37mnzyRexqzslnYiIY/t87y7tPqGQJEXEinYf36cXffjZlNKXaPeO8PMl/YOLexQPtfH4iNh7TfmcPW1c6HG6Ku3e0d54ge8/+NrypXr8no8/R7u/Obh/8fm/kLQu6ZUPY//AgSjGuFasaffOcRwRz5L0LXu+9yZJ37B4A9hAu69H7r2w/ytJPxoRT5CkiDgVES98uB2KiEFEvDQijqaUZtq94D94Z/lLkv5aRHzl4g1kN0bEUy+0r5TS3ZJ+S9JrIuJ4RLQR8dzFt399sa9nRMRQ0j+T9K6U0u0R8b8tfmvQStqSNN7Th4vxLu3e8X//ou3nSfoGSa/PPM4Ljc1XRcQzI6Je3Km/WtIZ7f5W4sE3tx1fvAHtWZL+rqQ3X0K/H/Stsbucalm7r++/KaXURcTfkPQVkl563l0/cFlRjHGt+NuSfiQiNrT7GvFDd2YppQ9p9/XQ12v3LnlTu6/rThab/Ix276p/e/Hz/13SX3zw5xfvHH7pJfbr2yTdHhHr2v3V+UsXfXq3dt+c9lPafaPS72rP3e0B+5pJ+sii/9+z2Nd/lvSD2n1j0t3aveN+8eJnjkj6Be0Wuju0++vrf36xDyKlNNVu8f0a7d5RvkbSt6eUPnLQ4zzAMe0+iTin3V/XP1nSC1JK48X3X6zdN4xtSPo3kn4ipXTLgz+8eOfzl1/EQ/gV7b5G/hntvlnt7y6+/hJJnyvprj3vqP6Bi9gvYHnw3YQAFiJiVbvv4n1KSumThbuDKywi3qndd0//Yum+4NrFnTEgKSK+ISKWF6+n/qSkP9buu2cB4IqjGAO7XqjdNyHdJekp2l3C8oj7tdF5ARR7/13Mr2QfESLiBy7wWH6rdN+Aw8avqQEAKIw7YwAACmsOtbG6Tm3bHrhNPKylgp8tRf6uv0+Xr70w2pOkZDxGZ1chr70u8u3Vzr7MX6IYzcleEmr85qb3GvSaM7YJY4FLVN5gOcfZ/t2VMQ7J2FvtHmdjG2stkHn4eqtf3s6cu5DO2FVlHp3KuM50zr4u31SXOUXlXCKdY+Oept61L78z99jML+P1yrmObo8n96eUTp3/9UMtxm3b6slPeOKB29S9e8Ty203beXabycwbgt44e0d1vj1JmqQ6u01rXMWG3tVJZwf5sTrWddlterO9tnYudeYvZbr8QGxV+X05T4B2t8trp/mtYugtSR0ZYzo2+/7ZeRv7m9b5fh3zrk5q81NmNxYto2q8sdro8o+vrry+D/r8ObhujNVaPbXaG40PvgmRpPUq315q3Ccb+YOz5DzbkDSp89uNjX0NzGd51Sw/DqM6f90e9t71+HSTn1e9cc2WpCPG088/+shH79jv6w/r19QR8YKI+Gjs/jWYVzycfQEAcK265GIcEbV282y/Rrsxei+JPX/8GwAAeB7OnfGztPuXUD6xSN95vXaXhwAAgIvwcIrxjfrsv3TyKe0T4h4RL4+I90TEezrjdUkAAK41V3xpU0rptSmlm1NKN9e19yI4AADXkodTjD+tz/6zYzctvgYAAC7CwynGfyjpKRHxpMWfnXuxdv+yDQAAuAiXvM44pTSPiP9L0n+SVEv65cWforugPkLb7cH1fzm815XbmbG2NOUf3sBc+b5urGs+NvV+DT8f5NucG+215m/9l4whPWPs7IS5bm/LSAmYmuN+nbEUcmTMmZkxFyRJXf4xNpk5LElTN8jCWDtbDb19OUcnjPNmbK4nnxvHcDO/vFaTbmC1d6rPH+eJMZ6SNKwm2W1GkZ8zk7l3cFJrnPPG2tmpGWRxfJLfbmKuyQ5jPsya/DE8Ps2PuSTtGP0aG8k7tXmN2Yn8tW9ins8TK2Nhfw8r9COl9DZJb3s4+wAA4FpHNjUAAIVRjAEAKIxiDABAYRRjAAAKoxgDAFAYxRgAgMIoxgAAFEYxBgCgsIcV+nEpuu7gKJONxuvSUPk0ng0jgWVWeTFWp7pZdpvt2ogbktTO831vayNPyUyVGhtBNPXciJjpvfYmTb7BtvfG/S4jHee6Pp/s04aXxjMzTolRyh+/2ozscZLWjFCwXclJBjOef8fYbC4/3xtjXg3zp5YkadzmU5fmZkLVxBj3Jhnz2ExmWjKuRavKP74HknfebDvz2DwnnLS8k8Y5OK+9wWqNBL9tY1f3D7xkt1mXf3yP66fWvqbTSy+p3BkDAFAYxRgAgMIoxgAAFEYxBgCgMIoxAACFUYwBACiMYgwAQGEUYwAACjvU0I8qJa1kAhNmTviEpDAW9y/3+ecayQgPkaQwAjaayC/al6Sx8gv3t6r8OPTmWJ0wHuOmEX6y5T08nerz/ZqZ4QxryocJOMM+NtsbVfmx2jaOX28cP0kazfL9OlN7z5lPJic9I7+vceu1V/f5gR8b4S7dyDs2y7N8v0ZGMIgkzYwAimk/ym5zzJifkrTR5PsexuNbDi8BpnPSSLz8EM1m+Q3PGAEi88a8Xk3yx9A5b2bepV2NEQCzqaG1r9ZtdB/cGQMAUBjFGACAwijGAAAURjEGAKAwijEAAIVRjAEAKIxiDABAYRRjAAAKO9TQD4VU1Qcv1g4jqEOSxvUgu83cCCVYS1OrvUhGSIAZODCc5hes5+MGpLr22tvo8n2fd/l9DY3F8ZK0ZKRwTIzxlKRZlx+rNuXnzLYR8iBJM2P+tco/PjecYTPyR/o6YwwkqTPGfanKP75JeGkQrZFvUDfG8TPGU5J2jDCSM2Y4w8kuf3y2javjmc5rb2jMh+1Bfo4eG3vzeGIE/Yw77zg3Vf741Cm/r9R5gRibAyNAZJI/OLV53jSDy3dPOjT6fiHcGQMAUBjFGACAwijGAAAURjEGAKAwijEAAIVRjAEAKIxiDABAYYe7zhgA8JD56pru/KFXaefJn6eUpMe+6h9r6QMfKN0tFEAxBoBCPv39r9Daf/t9PekffK+2moH60VLpLqGQQy3GKYXm84OTWtYqLxGrNpKS+pR/eEvyEm2mRnpMO/f2dY/y6WHX95PsNitG4pIkPdC02W1G8/zjS/ndSJLuMY5NbSZiDY3jU/f51Jtjc2+seiOUaGDsahDeK0D3GO11ZvJZE/n5vmmcE8fNpKSuznd+bZ4/NlPzxTInCepo7SWfbRhz9HF9/lo0k3dSpH3a61ZXtfnFN+sJr/xBhSqd2J5I2wef9/cbKWSStGocwtWZlxZ1n3GcR5FvcGrMBUmqjJS4JSPZzUlNlKTGSGoczLx97Qwv/ZVf7owBoIDJjTeqOXNad/zIP9HOF3y+Vj/4IT3xx39c9c5O6a6hAN7ABQAFpLrW9lOfplP/9o162l/9q6p3dnTXX//rpbuFQijGAFDA4J57NLjnHq388R9Lkk789m9r62lPK9wrlEIxBoAC2gceUHvPPRo/4YmSpHPPfraWPv7xsp1CMbxmDACFPP7Hf0y3/9iPqW9bLd15p578yleW7hIKoRgDQCHLH/2onvotL5EkDYwVDbh68WtqAAAKoxgDAFDYof6aOpQUmaCKs7X3/GBorB/fNhIcqvAWoqvJLwxfm3rBEmry6wg7Y435xAwZWTXCEhrjedncTGfoamNBvtn3qTHuTrfWOu84p1l+Z+tL+W2Omo/vSJX/1WSYz5mdrQb9LLuNmSshGaEzXZ8fhyp5581A+fN5xzwF68YIuTHm3nzuXULzET5So3yf1mcjq72llA8/6Wp33PPnzrzKb3O6Hlrt3TAdZ7d5YJCfC23vHZvWODrj1kjnkbzUoAvgzhgAgMIoxgAAFEYxBgCgMIoxAACFUYwBACiMYgwAQGEUYwAACqMYAwBQGMUYAIDCDjWBq4rQcmTqfz44RpLUG6kwJ4xEm5m8pKQwIrHub1prX2tGv6LPP74HqoHV3nLK7+voPJ/MdNpIJJKk48ZYjSvvQDfz/JgmI0Vt2njH2UldGhhzZj28udDENL+RmRI3NtKuhrnzT1L03h8smFb59qbDfHteLpM0m+XTjZbNVClnrE5P8+fXwDi3JKk1EqruNtIHb0r5dCpJmhqX9qFxHZKkc21+XzdM8/vqWu+c3075c2fVmKMD5a9pkrTT54/z2JjrktSYbe6HO2MAAAqjGAMAUBjFGACAwijGAAAURjEGAKAwijEAAIVRjAEAKIxiDABAYYca+pGSNEsHL55ujFACSVKXX9y/nc8I0NDYjyRNjNCIiRkA0BmBEEuz/KL21cZbtH/EeIwby/m+r3lDpbExDBN5oRgyAkQmmTklSaPkdX5U5zu/Oc+3t2KGTwwifwqOzSCcHSOYIKp8v44aOSSSlIxxP21sI/Ocd7Yam+fgyBjTvslvtGoem23jSnvECCJxIyXOKn+cl8MLslhJ+QfpzL2BcU2TpC2jX6vGPF7esZpT3+TnTDKvH7VxvboQ7owBACiMYgwAQGEUYwAACqMYAwBQGMUYAIDCKMYAABRGMQYAoDCKMQAAhR1q6Ick5XI4xuai6XubfNdr5RernwgvJGBmLO4/GV4CQO+MepNfZB5GEIkkzY316qNJ/nnZtPbam7ZGgMjMiy/Yyc4YaV7nt0lmuMvI6JaR+WGNuSTJCByw5oukofHUetTl5+i2vOM8NkJ1Vo3Qj6kRtCJJinxoxNwIzpCkHSM0YqvND+jcHCv1+YlVpfyA1sZ8kaSlWX67eWUcQElhBKk0Vf7YnDMncmucO50RapKWrObUdPlxWDKn6NaQ0A8AAB61KMYAABRGMQYAoDCKMQAAhVGMAQAojGIMAEBhFGMAAAqjGAMAUBjFGACAwg41gSuU1GSSaLar1trXsSqfaDNL+ecag9p7PrJtJEFthrevI0Ya1FLkk1zuDy9BZ8NJeTL6NLKDkvL9mhjHRpJGRt+ddK0jMy8Z54HW2S7fpxhazanp88lFJztv4O+P/Okcxjxeb71zcGwkra0ax2bHvAwty0h5MtO8ThjJYKvGnJmayW5HjTmzGc51wXt88yq/3UZ+OCVJa0ba1dSYV0M3odC4Js/TILvNaaPfkrTS5AeinXvj3hvz6kK4MwYAoLBsMY6IX46IeyPig3u+diIi3h4Rf7L4//iV7SYAAFcv5874dZJecN7XXiHp1pTSUyTduvgcAABcgmwxTin9nqTT5335hZJuWXx8i6QXXd5uAQBw7bjUN3Bdn1K6e/HxZyRdf6ENI+Llkl4uSQPjzx4CAHCtedhv4EopJR3w9tKU0mtTSjenlG5ujL87CwDAteZSi/E9EXGDJC3+v/fydQkAgGvLpRbjt0h62eLjl0l68+XpDgAA157si7gR8euSnifpZER8StIPSfpxSW+MiO+UdIekb3Ya6yO03R78q+rWXKjdGM8jossvwJ7V3kL0U3W+X7PkLQx3FshvRX4hejX3xqo1xqrt8tvsmIEKoy7f984MSHHCC6ZG32eNtxjf6dbY2NcJYwwkaXmWb/CB1nuvxWrKz4dxys+9duCN1aoRZLHTOGEQ3jyu+vxYnTSDLJrZNN+ecS2aVd483jACWfo+v809M+96dbzPH5sTZj7FeJbv1zSfwaFaXoPXd/nHOK/y+5qb15ixsa/OvH7MH8YLv9lRTim95ALf+spLbxYAADyIBC4AAAqjGAMAUBjFGACAwijGAAAURjEGAKAwijEAAIVRjAEAKIxiDABAYYf6Z5R6SZuZFJ1MQNdDui6fjuMk+6TOiI6RtKl8KkxXe8O5bCTMGMFTCiN5SpKWq/w4bNT5hJkmme0Z477hPEBJk1l+G+cIVkaimSQ1xpgeMXY1NBN71kf5CT+X1/dZlZ9/UeXbG5hP0edGItY88sd56AVwaWWeH4dxPlhLkjRr88dn2UjU2zbPCWMa67o+3/nPVN71asNIFhw1ZoLfxLjWdsZY1V7f14f5OXrE2M/UmOuS1BvXorr3xmrbSH28EO6MAQAojGIMAEBhFGMAAAqjGAMAUBjFGACAwijGAAAURjEGAKAwijEAAIUdauhHpaTl6uDl7314XQojyMLIBbEW9ktSF/nnLWeMEARJGhj7Wkn5gINNMzhjamzX9Ma+Om8R/XhoBGfkMwkkSVNjDX2b8jt7YNBa7bVGvsYRI+BgbD7PnUzy+6pGXt8ny/njMzLOryXz2NSD/IZTIwCm3fbaS/N830fGeSNJM+P4bPX5yTey4jykM0YgUF3n+5TmXviEE/Sz3nnzamAEAtVGuMuJ3js2O8ax2Wzz42lchiRJa12+X/c0Xl1qzYCe/XBnDABAYRRjAAAKoxgDAFAYxRgAgMIoxgAAFEYxBgCgMIoxAACFUYwBACjsUEM/kkLzdHCTZypvIfpKyi9Ebyb5bWa1t4heg/wmN/QTa1fO2vfTxqL23lhoL0mNkZzRp3xgxKDx0iC2Jvl9bZiBJa2x2REZfU/e885R5B+js6x/7HRcUtXnj/OaMdclKRkhFauD/DZLZt8nxrjHND/uVTe12uvq/MjPzXk1MEJunFCdaWNer4z2OuWvRQPj+EnSejXM76v3zueZcW1YNc6brdYMRTLCk+pZfi4MGm+swiiDJ70pqqERIPLxC3ydO2MAAAqjGAMAUBjFGACAwijGAAAURjEGAKAwijEAAIUd6tImAI9eXTvQu37m59W3A6W61snfvVWfe8trS3cLuCpQjAFYqtlUz/rev61mvKO+rvUHP/sLuu7d/01Hb/tg6a4Bj3r8mhqAJSQ14x1JUmoa9U0jJS9kA8DBDvXOuErSciaJpluc7Dmt8sk3x9pZvk/h5ClJm5Efqo3aS+MJJ2knjKQkI9VHksZOKpGRoLN7Oc5r2/zjG3tDpROT/L5OG9N4OPPGamY8xGQ8vmPG8ZOkTePpcJ28lLh6ZpzORgpSXV24vVRV+i+vuUVbN96kz/mNf6tTHz74rriZG+lvlXdPcKzKJ9zdOfEm1rLR5KjOz5m5mcxU50OsVBmn4KqRQiZJm0bfw0x2uy7y11Fn6s3NJ24T47p23Lhu98bck6RzdX4yzM1r7TEz0HE/3BkDsEXf67l/89v0l1/8DTr3tD+njSc9uXSXgKsCxRjARWu3NnXife/R/c/60tJdAa4KFGMAlsnRY5qtrEqSusFQD9z8F7Xyp7eX7RRwleDd1AAskxMn9f5/+I+VqkqKSte/4+16zB/819LdAq4KFGMAliOf/Jie+ze//aHP53PvzUQA8vg1NQAAhVGMAQAojGIMAEBhh/qacYqk/oBAAUnq08Da18xY/N7l16rr/sZYjS9pUDmvj3nPbcIICphZi/a9FeZN5Pe13ucXyM/ljdWa0a958vY1c8ISlH98yck0kTQ3AgBqYxwmnbn63wheSI23r1ZGOMM83/dpZZw4kmKeH/dmmn98Z4y5J0nR5QM96sYLZ9g2tqkjPxeGVmtSMo7zuaX8ce6n3nkzMQJ6KrPzg538tW/DOIYD4zokSZWx3dl5fi6sWtdsaTnl5/u2GZBSm2E/++HOGACAwijGAAAURjEGAKAwijEAAIVRjAEAKIxiDABAYRRjAAAKoxgDAFDYoYZ+RMoHXiw3Y29ffb7rm5WxaH/mLdJOfX5fg8pc1K78YvQdo7154z2Xqvp8mMCa0fVtc9H+lrHNcuWlcGwZqR9dmx/Pygx3aY2HODeO85Yx9yTppBEOsmHtyQuWqPt8wMFk7B3n3pgP02QcG++UV9Pk58xKMs/BLj9WM2NXEzNkZDgzgjNS/pp2nfn4nJiYauzN0WmV71fM8/vaNoKMJGnJCLlZ7vJz4T6zvK0YoR/JCBaSpHVr5PfHnTEAAIVRjAEAKIxiDABAYRRjAAAKoxgDAFAYxRgAgMIoxgAAFEYxBgCgMIoxAACFHWoCV6fQVrQHbjOtvaSk6/p8ok1jJBJN5LW3M8y3F2b4ytH8rjS3InS89LDeSMdZro1OGalFkjQb5bfrx95gnTASnNbn+fai9dobGZv1RmpbMlLWJEmR39eo9sa9TfnOr07zaUPb4Z0TUyNtqDa2uS68sdpK+Xm8YSafnTKuH8mYe23njVUYx3DZOOl3Gm8uPG6WiTqUtGGez5tLRpKhEajXy2tvMsu31xmBWMudN69aI9VsapbKS8/f4s4YAIDiKMYAABRGMQYAoDCKMQAAhVGMAQAojGIMAEBhFGMAAAqjGAMAUNihhn5UlTRYOnghdjXxFtHPmvxC7RVjYXhXeQvDl4wAkWriPbe5fznfZpfy43BwfMr/MjDGaqsfZLc5mSZWextGyMjEnHoPGH1fNY7zGfNp59QIZzg+zSccDM3gjNTnO1Y3XpRA9PmBWG/z7c16r715n3+Mc6NPO603F8IIljhu9n3DCKCojECggdnexLjvaY3jvGZmyWwrfz6HOa9a5Y/h0EjhmBvnsiStd/mxWnGu217GiJIRvFMrH5YjSVsP4/6WO2MAAAqjGAMAUBjFGACAwijGAAAURjEGAKAwijEAAIVRjAEAKIxiDABAYRRjAAAKO9QErpSk+SyTdmImYhmhKZoaSVBd8lJhmjq/XQy9lJajXb7zEyPxq62nVnuOdSN5amykN0nSRp+fVkdqL/2navJtjvr8nLlu6iVibRqb9Ub02frMG6u6yTe4Zu1JGlb5MZ3kzj95qWCS1MpIIuvy22xbrUnJmDKRvL5XA+N8NpKgTo+8S+jA6HsyrkWfMRO4Bk3+OC+ZEVXLxgHaaY3rqHnrd8KoAdPeSM2qvGu7jGv7LLy8wyq869q+P5vbICIeHxHviIgPR8SHIuK7F18/ERFvj4g/Wfx//JJ7AQDANcx5rjKX9H0ppadLerakvxMRT5f0Ckm3ppSeIunWxecAAOAiZYtxSunulNIfLT7ekHSbpBslvVDSLYvNbpH0oivURwAArmoX9ZpxRDxR0jMlvUvS9Smluxff+oyk6y/wMy+X9HJJappDfYkaAIBHBfvd1BGxKunfSfqelNL63u+l3Xce7PsqeErptSmlm1NKN9cUYwAA/gyrGEdEq91C/KsppX+/+PI9EXHD4vs3SLr3ynQRAICrm/Nu6pD0S5JuSym9es+33iLpZYuPXybpzZe/ewAAXP2c3xs/R9K3SfrjiHj/4ms/IOnHJb0xIr5T0h2SvvmK9BAAgKtcthinlP6rdMHV4V95MY0lSTMdHHIwmHvBGb0RGpGMAIe28xaGz4ynLUaOgCRp1Vjc31f5sIQuM5YP7Svyr0YMmnyfpmZ7oz5/bLrKCxzYMUIcBkYQQmM8Pkl6zCzf92mb73s19NpbSfnjPDODLGbGq05Nlz+/3CCcfm6EjBihLWZzqvd/W8pnmZnHecXYLBmhEZURLCRJZ435vmI8vo1qYLV3ztjXYyM/9ySpGxjnoDFWs947NtvG9eqkka0xMUNNJsZmXecFLB0zr5H7IQ4TAIDCKMYAABRGMQYAoDCKMQAAhVGMAQAojGIMAEBhFGMAAAqjGAMAUNih/uWGKiUtp4NDB9olYzW3pNZYP76VjNXc5hrtLSd4IbxF7eeqfErFZpPv2MrEa28p8mN6X7TZbbyIAKke5Ntbm3p7WzKCAmbGWG2EFwCgUf44n1D++G2bE2u9NsIgnHks6Uz+ECoG+VM+JkaKiiRnSDe7/DjUZlDHwEkHMc9BJ+ih7/IbNa13vVoz5vFOn597R6odq71JnR/3c50xYSSN6vy5OjDCkyoz6GfZCCyZ1vmxmhnXPUlaNUKKGrPvg85rcz/cGQMAUBjFGACAwijGAAAURjEGAKAwijEAAIVRjAEAKOxQlzYBuLze/da3qd7aUvS9+n6uz/+Ol5buEoBLQDEGHuW+8G98l9qzZ7XuLRsF8AjEr6kBACjskO+MQ1UmkWdqpN5I0sBIQTJCaLQT3vORJSMVZtVM/2mn+cc4rvKHZn3gtXdKB6eeSdJa5FN2jNAbSZIRwKVkJto4SWtN5OdCY6ZYrRtJawMjKWlunlmt0a2oD3p8vT74mtdISTr1G2/UY37jTQfuazgx0oacpCtJZ40IrqPGsRkmL41t3bhcOfNFkgbGXJ7V+bGKuXdS1F1+u6Gxqz68ZDcj5E/DyhusRvnjnIxt5sY2khRzI83LSNfadAqAZPVqnrx93WcG/e2HX1MDj2LP+K6XaXjfvZoeP6EP/IvXanT7J3Xkfe8t3S0AF4lfUwOPYsP77pUkDc6c1vF33qqtP/cXCvcIwKWgGAOPUt1oSfPl5Yc+PvfsL9PSx/+kcK8AXAp+TQ08Sk1PXKcP//OfliSlutaJ//g2HfuD3y/bKQCXhGIMPEot3fUpfclLv+mhzyfmG68APPLwa2oAAAqjGAMAUBjFGACAwg7/NePMy1or+XwKSdJ6m8/+c55prCZjdbykqbPY3gyymBqjPqynRnPe4Rsb2zl7GiUjzUNSZYQJRPKyG1snSMUICZibUZEjo73kvDZrhBJIXojDzAickaTBJD//dozwk8YI85CkoRG8s2rMmS3znmDFmKQz8zjXfb7vbZ8/NltmCMd1yp/PnREmc9oYc0kKJ03Gy1rR3Ag/mTnXmM47J8ZGENNane/8ytwbq+06P1Zrvdf3WX/pJZU7YwAACqMYAwBQGMUYAIDCKMYAABRGMQYAoDCKMQAAhVGMAQAojGIMAEBhhxr6kSRNqoMXTw/MpwdDI7ygz7QlSZp5DdatEQbReeEMO8YC+WQ8T1oxwxl6I4RjkvJ9Ghn7kaSuyY/76szrexhBAZuDyzeNT1T5MIHNzgg4MIJIJGnFSF7YkhcsMW3ybS6nfHszJzBC0pFJfpvG2NfmwJ0L+e1G8vZVe0OatWSET0jSwAgEWjfmzKAxr1fGMGyZYzUwAmyWZvltts3r1VqdD+vYNg5gb4YwOVkk095LoxqZj3E/3BkDAFAYxRgAgMIoxgAAFEYxBgCgMIoxAACFUYwBACiMYgwAQGEUYwAACqMYAwBQ2OEmcIXUDTLJKZ2XdCIjPWZqJNokM9Hmhsj3a92JvZG0ZCTMhHFoHpPy+5GkjbrNbnPSSA8LMwlKyRjTkRF7I6nt823GNL+fsdMnSX3kt6uMYRgZaVGSNDP21ffevtrIp0GNjYi71kySS8pvNzcSiY5arUkTZ9x775wYGPchDxhjdWpqTD5J28Y5uGyMZ29eY4Zdfi60tXcO3lcNsttUg/y+Gi+sTDKuM1WXP85HzHOwq/PjvtN4176Ree7shztjAAAKoxgDAFAYxRgAgMIoxgAAFEYxBgCgMIoxAACFUYwBACiMYgwAQGGHGvoRkgaZxdpnKu/5wWqdX0G+PM8vVh8bARyS9EDKD9VaeMPZKt/m0FivvtUOrfZk9CvymQRqjRAESVpKTniBtzi+7/NhAq3RraGxH0na6vNjVRuBM31jhtdMjRAOc45OjX4Z2QxqjaAOyZsz00F+XyfkpUEsO2E5ZtBD3eS3OyVjsGovDGJoTPd2mn987dy7xoy7fL+SmeGzFsY5OMtvs1F5fe+N8JptI5xnI3lz4egs39760Lv2xcxrcz/cGQMAUBjFGACAwijGAAAURjEGAKAwijEAAIVRjAEAKIxiDABAYRRjAAAKoxgDAFDYoSZwJUmT6uCEkmUzNWVqbGcFCc2NGCFJQycFyUyVCuWjb+atEY9Te8+llkf5gVir8+MwarzInuXIj0MynwZuy0jzMpKLurGXYjVP+e3qLp82VFXuXMjva5o5Zx407PPb1caxGZhXhUmbb2/VSGOrG28yGM1pxWhPksJIlZoY56nZnGadMQ7OBWvgzeN+lh/T1HnjPuiNZLA2P55b8lLpnBS140Zo245xjCXvWrQ69/a1PbI22xd3xgAAFEYxBgCgMIoxAACFUYwBACiMYgwAQGEUYwAACqMYAwBQGMUYAIDCDjX0o1LSUibkYGgldUhL0/x2m0ZGRVMbq8clVZEPxRgPvASAtst3bBD550nLAy+EY21o7Gs4zG4zWvICUlqjvdoIu5CkY9P88RmP82ECW403r2InfwzHxn763nue64SDjMzwgqbJ76uL/JxpWm8eHxvkj83AOJ+XzHM+jKtV5Z0SGhgBKU5ERdd7YzVRvr3taf4BLnVe6Edb5+fMzAjLkaSzXb5fY2O+r8699sbz/L7aKv/4to1jLEnjOv/4rpt7gSUbRlDMhXBnDABAYRRjAAAKoxgDAFAYxRgAgMIoxgAAFEYxBgCgMIoxAACFUYwBACjsUEM/QlIue2FqBhxURv5EGGECy7W3MHxY5xesV8ZidUnqjb63bb5fwyWvveHycnab42uD7DbLx0dWe6PV/ML3KnlhCd12PlhiYz2/IL8+t2O1l5TvV5pMs9usz7xTa8kYhqjN58xNfo6uyQnq8JpbM86dGOYfYG2egyvGMFTmvmSE6nQpv6/eDJaYdvntaiO0ZXvmzQXn9Bon71q7bGR1TJzzxgzEGNX5fe0YXR8ZfZKkiPwD3DZDg5bMGrAf7owBACgsW4wjYhQR746I/xERH4qIVy2+/qSIeFdEfCwi3hAR+VsrAADwZzh3xhNJz08pfZGkZ0h6QUQ8W9JPSPqplNLnSToj6TuvWC8BALiKZYtx2rW5+LRd/EuSni/pTYuv3yLpRVeigwAAXO2s14wjoo6I90u6V9LbJX1c0tmU0oPvBvmUpBsv8LMvj4j3RMR75uZfHAEA4FpiFeOUUpdSeoakmyQ9S9JT3QZSSq9NKd2cUrq5qS/9z0sBAHC1uqh3U6eUzkp6h6QvlXQs4qG/MHqTpE9f3q4BAHBtcN5NfSoiji0+XpL0VZJu025R/qbFZi+T9OYr1EcAAK5qTjLBDZJuiYhau8X7jSmlt0bEhyW9PiL+qaT3SfqlK9hPAACuWtlinFL6gKRn7vP1T2j39WNbn0LbmYSSFfNlZWdRc9vnE1g6MxVm20jHWQkv8WXN2G7VSF1aGnh9P7Kcj/xaum4lu83Jxy5Z7a2eyD/H6ysv0Wb7XD7tSndvZzeZme8drLv8sRnM8ylWa8qngklSJ2McjLknSbUxpkeN5KJ+6L16NTD6tWRMmcZNNzIen5se1hsPsTdirOadN1aTWX7gqya/TWOkCkpSb8zjJnl9P2ME7w3H+W3OmcemrvLjUBvpWmvmNeasMaYjI41NkkbzibXdfkjgAgCgMIoxAACFUYwBACiMYgwAQGEUYwAACqMYAwBQGMUYAIDCKMYAABTmJHBdNlUkrbYHL+g+O/e6tKp8ioMTErCTz8OQJM27fL+Oy0hUkFS3+bCOqbHNY8IL/Rgs5fu+eiQ/EEdPGav/JR197DC7Taq9RfTt0JgPm/lxn5zLB3VI0nadb6+L/Fj15lwYGOEgThCJJNVVfruZERQzzB8+SVLb5o/hcJB/vn9k4I1V0+T35YZ+OFkrcyM0aOJ1XY0x32OWPzbGIZYkjY05MzFCTSRpqc/3fXuQvx7Pwru2T3by26wZfZ/seO0Njfk3Nm9b57UTR7U/7owBACiMYgwAQGEUYwAACqMYAwBQGMUYAIDCKMYAABR2qEubgJz3Xvc0/eJT/4q6qPTVn/oDffV9by3dJQC44rgzxiNGp9DPP+3/1A/90f+rf/n7P6rfu+FL9OkjN5TuFgBccYd/Z5xZrH1y5oUznBnmF8jHKL8wfDA1F743+YXhUXnPbcIIChgp36+5EbogScM6v6/a2Fc98kJG6qX8wvf98ko+tvYEPW7ygB6vdWlU6Xn3v1/vf9IX60kf/U8H72uQn8Z15U31QeTHYWIM+8xJlZDUGe3VZjDNijGxoskfw5XGOyfWjH4NjZyYJSMYRJKWnZAR84pWGefXrDPCLibeWE2M4Awng2Nuhn7sJGNepXxQhyQN5vl5tdPn59VaeJ3vjPl3xriPbM32mj7/+E7PzTAPN3RmH9wZ4xHj/uFRnZqcfejzk9OzOr10tFyHAOCQUIwBACiMYoxHjJOTc7pveOyhz+8fHNOJnXPlOgQAh4RijEeMz9+4U3eNTukzwxOaRa3fPfVM3XzXB0p3CwCuOJY24RGjVq+/9fF/r1f++ZfvLm265916/MbdpbsFAFccxRiPKM86c5ue9d7bHvp8o2BfAOCw8GtqAAAKoxgDAFAYxRgAgMIO/zXjODj1pWqMyB5JjZFctLaTT5ip5SV+9U0+bqjKPLYHzY2nQGEk2sydyB5JMyNhZmYkn23veGNVb83y25hJNeOtfJvdTv7xdZ0ReyZp3ufH1Bn2MNN/VowUpGHnjftslJ9Yq0a60cALWtPAOIirVb69I2aS3NpqfrslMyXOuHxoMs3PmXbLO+fXjQSuiZEKVhvJWpI07PL7mpiJga2R7LaU8tuMw2svGWVp2TiftxpvrKpZfruT5rX9rDOxLtSPS/5JAABwWVCMAQAojGIMAEBhFGMAAAqjGAMAUBjFGACAwijGAAAURjEGAKCwQw/9iO7g+r9RewEHVWY/krRtrP+f197zESekwlmILklz4znQ2AifGBlBHZK0MclvN9wYZ7cZ3OuNVT+ZZrepkrevnfX8Yvv19Ul+G2MMJGna5dubVPltZkaYx678xJq2XpDFdUY4Q2sM+9DMLWiMk2JYGUEdA28uHFnJX65W18zQD6NfW0aYTOq8wZqMjbAVI5zHDWSZXMbjvNHkd1Y7ITe1F4Sznc9X0tLMGE/zerxk9GvdCDWRpNrbbF/cGQMAUBjFGACAwijGAAAURjEGAKAwijEAAIVRjAEAKOzw/54xADwCvfIrXqnfe8JzdGLnjH7z335L6e7gGsOdMQBIetH/fKv+1du+p3Q3cI2iGAOApJvvfr+OjtdLdwPXqEP9NXWn0Nk4OF5lZCRPSVJlJNG0RrpRMlNo2nH+ecvWyNtZM893fmo8TdoemwkzRrrWZmWM+9Rrb3LamFbJixKa7+STrLY284lfk418Spck9dv57eZGso8bwGVMBS17p4TqlN/QOeHnzsklqZLRnjGPq4HVnOpB/vxqh2YC1wXSw5phpYhQO6rVzPOPr3Ki+SRVRspTivxgReNNhqGRxrZp7UmqjWvyljEOjXltHzqXGePxTTqvvE0jP2e2zMRAZ6wuhDtjAAAKoxgDAFAYxRgAgMIoxgAg6Xuf+yq9+Gtfq08e/Rw995t+U7/5tG8s3SVcQ1hnDACSXv17P/RZn29uen92E7gcuDMGAKAwijEAAIVRjAEAKOxQXzOulLQUB6chtL23aL/t86/nVHX+ucbqzGpOUV++14/CSHqoU75jM/PondnOL5CfzfOL9uc73hhstAcHu0hSY6atpGl+Ef32JD9Wk22v7zuz/FhN+3yiR1V5ASlOcEZvBgnMjeAWJ1iicRNLjGPYGV3vjHANSZrO8tuNzWCaqPJ9n03y7c3MoTK6rtTnj83cbK9TfsOZMQaSVBkBG4Muv691MzjDmX+TTHiUJLVmQMrUOL9W5F0/murS72+5MwYAoDCKMQAAhVGMAQAojGIMAEBhFGMAAAqjGAMAUBjFGACAwijGAAAUdrihH0lazaxaPz3IL+aWpBPGIvPOWMy9MbSa05F5fuF71XnDmRojhMPYz2zsJZaEkQ4yHY6z2/QT87mbFZBihn6k/DGcGkEP1XxqtbdtzKt5nQ8lqI0wD0kaGOEMdXj7MoZKvdGvZIaMOEEWEyNMZjz2gjq2tvLzvTfOU0lS5Ps13TECZ7a8sRpPjfaMcXdCcCRpq8sHC82SF2TRdfnzPjojCMc4tyRp0wgjcY5ybbY3Ms6vuXNySVLvtbkf7owBACiMYgwAQGEUYwAACqMYAwBQGMUYAIDCKMYAABRGMQYAoDCKMQAAhVGMAQAo7FATuJJCk3Rwk40m1r7WG6Pr+VAYjd1glTr/vGVkpPpI0rDPd2xqJLnUZtjQjrHNfJJP7BnMjAGVNDY6NjBTnmpjs85Ixxn0XtrQrMr3veny7Q2TN1ZOuFbnhdJpaqSaDfv845uY6Wgj4yHudMZ4jt00tvw22403rypjs/ks3y8jFEyStGnsa+qkdM29x9dP8uM+n3j7ckK/xsaAxo5572ecg87lozXvNbea/HabvVcqu0QCFwAAj1oUYwAACqMYAwBQGMUYAIDCKMYAABRGMQYAoDCKMQAAhVGMAQAo7FBDP/pK2h4evM3IXNQ+MxIAJvUgu83Rmbdqvxrk26vMIItG+e02ayPAIczUDyNAZNnIxHBCFyQvOCOZTwNr5fe1FfnAkqG80I/tNr+vNSMFwRkDSQojKKaZe/tK+a5bYTKjqdWcpkb4idGcbW6MQ1V5DTqhHxOj8xMv20VjI9Bjx7gU7ZjXxy1ju25mBu9M8w9yY54/oWf5y7EkKXX5sZo755cRcCNJ3Sx/4iyHd6CnZujMfrgzBgCgMLsYR0QdEe+LiLcuPn9SRLwrIj4WEW+ICPN5DwAA2Oti7oy/W9Jtez7/CUk/lVL6PElnJH3n5ewYAADXCqsYR8RNkr5O0i8uPg9Jz5f0psUmt0h60RXoHwAAVz33zvinJX2/9NC7aa6TdDal9OC7Yj4l6cb9fjAiXh4R74mI93Rz7000AABcS7LFOCK+XtK9KaX3XkoDKaXXppRuTindXDt/9hAAgGuMUx2fI+kbI+JrJY0kHZH0M5KORUSzuDu+SdKnr1w3AQC4emXvjFNK/yildFNK6YmSXizpd1JKL5X0DknftNjsZZLefMV6CQDAVezhrDP+h5K+NyI+pt3XkH/p8nQJAIBry0W9iJtSeqekdy4+/oSkZ13cz0vzeSZdpW+tfVVG0smS8qkpYT4fmUa+vdpMxFoy3sc2qJ19eWlDrZGgU6V8p87W3nQxwsM0NRJ7JEmR365WProotd5Yrc3yYzU3Up5mtff4Nvv8do8xE7icRzg1oqc2zLGqjSS5NMnvqzaTkuZGWpkzXyQpjDi5rsvva2xGjE26/GOcGYlY28Z4StJsbszjmbev8TSfUDUyrmljM8IvjCTDkdH1yp1XxgWr6bw3H0d/6e+LIoELAIDCKMYAABRGMQYAoDCKMQAAhVGMAQAojGIMAEBhFGMAAAqjGAMAUNih/uWGUFKdDl6IvdN4zw/mxgLyI3V+4XtvhC5I0mCaX0DetN6i9sYIhKiMbo3No9cZcRBbfX5hf197j6+dGwEp+UMjSYp8tyQjAGaezICUQX67VeXnQjgBFZKcPI9p5+3LyHZRY4TXGFNdkvTAMN+vI8aU2TaDMypjvtdGYIQkyZijRgaHkhGWI0mbRoDIdJpvsJ95B2fH6NaOGYpRGZsl4zrqXmvHxjwedPmNNrxppWiN89kMKaqN8+tCuDMGAKAwijEAAIVRjAEAKIxiDABAYRRjAAAKoxgDAFAYxRgAgMIoxgAAFHa4oR8ptDQ7eCV2X3lpEEeMRfvn2nxixDEzfGI58gvDKzNwYFLnhz2MlfYjc9H+etVmt9k2QgmWZzOrvcoIHKjDC0to6/zK/ejyx7lqvANdGXPG6JKSGWQhI5BlMneST6TeSFJxwlaOjr151XX5+X7WCPExTgdJUjPJt+dGLsyUH9N2np+jM/M4z4xAj6kRJnO2y5/LkhRGOMjUDBDZMgJzZsa8mss75weRnxDbxrxy8zeaLt+vsRlY0jbuef9ncWcMAEBhFGMAAAqjGAMAUBjFGACAwijGAAAURjEGAKAwijEAAIVRjAEAKIxiDABAYYeawJUkzTLlvzeSoCTpuJGasjPNp6F0rZdCMzaSktzBPFPlo2HWjGSfuvPSXtZSfqzCeF5WGf2WpGnk+1UZaUOSNDLGYWxEYlVObJaktS4fJbRkBGJNzLFq+nx722bft4yxWjWSks6a7S31+XmVjJS4wdRqzrpzCGPuSdLEGIehsY3ZdY2TMR+MgKrpzM0Yy497Z3beCepy0q6WzXOi6vMd640EtbF5fexSfl8Ts05URlLjBX/2kn8SAABcFhRjAAAKoxgDAFAYxRgAgMIoxgAAFEYxBgCgMIoxAACFUYwBACjsUEM/KkmjTPm/p/aeH0yG+a6vTPOLzJMZPrFjPG85MrN2pdZ4iGNj8fhc3iL6uXGYjxgBDgOzvfU+/wCnziBIOtc4YSv5bdZ6r+8DI/AiDfL7qs3gjG6SDxxITuqCpOVxPqWiNoYhGeMpSeutkX5iBGe4+so4B43QFklKKb/d3BgHI7NFkrTT58dqboTzqPcuMhtG+EkyxlOSpkZgSW9MrDp57e0YARuDaX48p4133jTG41s2w6g689zZD3fGAAAURjEGAKAwijEAAIVRjAEAKIxiDABAYRRjAAAKoxgDAFAYxRgAgMIONfQjlNRkQjaOmWumIxmhGEb4hBtkMTVCHO6O1trXiVl+cf+4yffdDSwZGUEBm12+704wiCSNjLHaNBfHN0auxHEjwGE49J53njKCAmJodMqbVuqMcJeZHZaQ325zmh/3iXeYFUaQRdtO89uYgSzOVhvWnqRt49IXlRE+0Xvn4Njo/LjLj+fQmOuSN1bTzhv3ythsasyFoTGekjRPw3x7TX4cTsy8x7cV+e0aZxAkzWfGteECuDMGAKAwijEAAIVRjAEAKIxiDABAYRRjAAAKoxgDAFDYoS5tAq5Wm9ddr1u/+59p59h1Ukp6+m+/SV/41l8t3S0AjxIUY+AyiL7Tl73uJ3XqE7dpOlrWm37yDbrp/X+gE5/6ROmuAXgU4NfUwGWwcuZ+nfrEbZKkwXhbxz/1SW1dd33hXgF4tDjUO+NeoXF/cJPDLp/YI0ndLN/1zkh52grv+Uid8vtaSvmkK0maZMZAkhojYWxkJnDNjHFYVr7vW5UZj9blx/Tk3It5qoz0qeEwPw61+bRzPspvuDw6eBzWTz5O9z/5qbrp9g9oMDh4X+tzI7loxxur6Tx/fMZGKt2OmTa0XOX71RmJX7Pw5pUzj1sz2a0y0uTmxnyfJC9xqUr5Me2NcT9tprEtG9ermHljtZ2Zw5JUG9er2hgDSWrn+XStSZPv+525k2/h6CxfcyK85DO518j9fvSSfxLAnzEdLuk//L1X6yte9xMa7myV7g6ARwmKMXCZdHWjt/29V+sLfv8/6PPefWvp7gB4FKEYA5dBknTry1+lE3d9Ul/8tl8p3R0AjzK8mxq4DO7+gmfqI8/9Bl33p/9Tv/Zjb5Qkfdmv/6ye9L7/UrhnAB4NKMbAZfC4j75Pf/clX/jQ52H+eToAkPg1NQAAxVGMAQAojGIMAEBhh/qacQppXB+8KLqbegu1e2OzEym/sH+n8Rbt9/P8a4A78vZ13OjXKOUXok9Sa7UnY4H87cY4HPcyRrRkBCpsGKEEkrRkPF0MI3ghenMxvhHCkbr8QHgxHdLYeGl5boYlzI3js2O8lj0K70A3RoMTI9DDnMWqjW6Na2+sZpVxrqZ80ENvvjVgOM9faneMwIhR482s6IywlZF3vXKOs4zwpLNmwJJzqrbG9bg2QmkkKYxjODDDVip3Quz3s5f8kwAA4LKgGAMAUBjFGACAwijGAAAURjEGAKAwijEAAIVRjAEAKIxiDABAYRRjAAAKO9wELklJB6farMlLMJmmfNe7Kr+vZn75UlpWjCQoSXpgkI+YWenz+1rr8glBkrTd5TOOTjppQ97D09R4jjeUl4g1VT79pzZSnqp8oJkkaVLn58PUOM4peSlW02l+YiUjbUiSGmOsBkbeVWUkjEnSunEOztv8vDphjIEkjY2Eu2nm+vKgfpafM12dj/lbjZnV3nadPycaI3qqnnp5ZXNjs865qEkaGuf9jnHOJ/PaHkaKVW2k/B2deRes2qgTEzOpsU3efNgPd8YAABRGMQYAoDCKMQAAhVGMAQAojGIMAEBhFGMAAAqjGAMAUBjFGACAwg419KNOSSv9wYECO633/GAlJvmNUn4R/bjyhqA1FplHeIEDp+b5BeT3GyEB/dwLzlhq8n3vjOdljRkSkIxDuGOEjEjSwAjr2DLCLpZ7b151RuBAmub7Xhv7kaSY5/u+PvGO88Q4PjNjjtad1945I2xlNMmP+9gM6hgb833WeOfzUeXDGZKRcrPVeu31db7vS8ZcOKN8EIkkLdVGcIa1J+mcMZVHRt8r8/oxnhvnqhGcVHnZNWqMnI4zRkiMJDW1F8qyH+6MAQAozHpaFxG3S9qQ1Emap5RujogTkt4g6YmSbpf0zSmlM1emmwAAXL0u5s74f08pPSOldPPi81dIujWl9BRJty4+BwAAF+nh/Jr6hZJuWXx8i6QXPezeAABwDXKLcZL02xHx3oh4+eJr16eU7l58/BlJ1+/3gxHx8oh4T0S8Z27+lSEAAK4l7rup/1JK6dMR8RhJb4+Ij+z9ZkopRez/VrmU0mslvVaSlkcj7+10AABcQ6w745TSpxf/3yvpNyQ9S9I9EXGDJC3+v/dKdRIAcHmkqtIHf+X1+pNX/2zprmCPbDGOiJWIWHvwY0lfLemDkt4i6WWLzV4m6c1XqpMAgMvjMy/+Fo1u/2TpbuA8zq+pr5f0G7G7wL+R9Gsppf8YEX8o6Y0R8Z2S7pD0zdk9hVRlFr9vGEECktRM84urB8b66+W59zq2tX68MYMljPCCuRFYsjPw2quN1e87kY8A8OIGpNYJ4Zh7r5DM6nxgSW0EiGx1XgJAjPPbTI25d8wMspDRr/uTN/KtcZzrWX6bjeQdm9U6/xhnxnjeVXvxE8eafHutGeKz3ee3q4xrUTfx2kvD/L6mKf8q3sDYRpImFzifp495jM4857l6zL/+Zd330pdqZoz9aJZv07mOnjWCOiRpo8pvt9zlz8HWDBYKo040RqiJJA3c836/NnIbpJQ+IemL9vn6A5K+8pJbBgAcqrv+3vfphp/7GfXLK6W7gvOQwAUA14D153y5mtNntPyRj+Q3xqE71GxqAEAZW1/0RVp/7nP14S97jtJwoG5lVXf80D/RE171g6W7BlGMAeCacMNrfk43vObnJEmbX/wluvdbv41C/AjCr6kBACiMO2MAuMas/tF7tfr+95buBvbgzhgAgMIoxgAAFEYxBgCgsEN9zThJ6jJhLtf3XoLJ3Ei7Cic3a/+/b7EPIyUoeSktW8auIh88paWRsZGkYZdv0Emx2vAenqo+f2wa5wFK2jL2tRb5vu8YYyBJ4RxDY5vtTNLcg+6NfLpW8oZKjfEQ+8iP57z2zolxm3+MI2M/q2a60STyl6uBmbpUGeOwacyZUcys9objfHtnRvltjvfesVma5o/NpnkrtmSM6aaRlDcYe+fEspHAdbTPtzc2jrEkzY3z66hZl8bVpd/fcmcMAEBhFGMAAAqjGAMAUBjFGACAwijGAAAURjEGAKAwijEAAIVRjAEAKOxQQz96hdYzYQjHcqkgCyNju1Vjm88MvDCIeeS3O2oGiCzJ2K7JL3xf3fHa2zIW0YexTdt5i/ZndX6sNuTtayAjcGCaD86IyhurqsqfElXK72s29dobtPnHN6m9UIxunu97b+zrATMsYXnHCEIwxiqMbSQrdkdd8ubVsJtmtxkM2uw2M/OcnxpDOujyx2bLPG9aJxTDzDty+j4ysk8q99gYKRznjGvMihmQMnEmltn3/JXowrgzBgCgMIoxAACFUYwBACiMYgwAQGEUYwAACqMYAwBQGMUYAIDCKMYAABR2qKEflaSV/uDF09tm4MBA+QXyO6P8fmozDGKQ8gEH29aepLERFHBikm8vMgEqDxqGEc5gLJBf8vJRtNPk2ztXecf5JiM8Y914eIO5GcjS5ufVtnH8ZskbrDTLH8OV5IV+VEZAyt1GeM3ADJaYGUEIc2Mclirv8a3O8ts5ARWSdLoZZrcZpnySxcAMltg0tjmRz7rQRuM9wHPG9fGoGfoxNM6v4Ty/0X1WbIs0NQKPhsb53NduKFK+X31jhlHNzEHdB3fGAAAURjEGAKAwijEAAIVRjAEAKIxiDABAYRRjAAAKoxgDAFAYxRgAgMIoxgAAFHaoCVx9kna6g9NVTjppUZLuHbbZbY5qkt/R2GpOmyMjJagzk5L6QX4bYxhOt95zqROR71dK+X3NzGPT7hgJY7U39daNlLGBERiVKi9VatsI0Knr/Hiezgc8SZK6Lt9gYyRdSdLQSM6KPt/3tjKioCTNq/w5OJjn2/MenbRhHMNh4/V9xThXJ30+maky0uYkacVI6toxzsHaTPw6Zpzzdyt//CRp1UiJa42EqjBSwSTJCMFTGOla5zov8as1rmtHOu84n20uvaRyZwwAQGEUYwAACqMYAwBQGMUYAIDCKMYAABRGMQYAoDCKMQAAhVGMAQAo7FBDP6pIWs0syq9n3r6GRpjAJPLPNbbNETje5xe+T+UtMh+l/IMct/l9pUyAyoM2jc06Y0H+srdmX+vDfN/XYmrtayPl0zOO9vkF+a25aD+cABEjBMENCdg2QgKSGYvhjOgxY5suefO4Ms7VkTGv5o13TzA2xqE3z8GBEfoxb/PtzZ0JI2nVOHfmg3yQRTXzQj9OG4Eeq2YIx1EjQGRiXB/H5q1faxznNhnhLkvesZnM8+dgZwY6OcE0F8KdMQAAhVGMAQAojGIMAEBhFGMAAAqjGAMAUBjFGACAwijGAAAURjEGAKAwijEAAIUdagJXUmgSByfk7CTv+cHQSLGaGElJlZmYkqy0mnzqjSTNjPSYe/p8e4/1wobUJSNdy3h8UXmPr0r5Md12k4uM6TCr8qlEM2MuSJKR6yP1+dNmtfeSkmrj2ER4+9o2EtkaY1dd5hx9UNUaiVFGEtlo5qUbTav8uBshVpKkDSP5bGAkTx2deUlrZ41rw9Dpe+9FFI5qJxHLnFd9/iSsq/w4rE6s5nR6kJ9/AyPxa2Se89vG9WrceqXy+ip/Bbn7Al/nzhgAgMIoxgAAFEYxBgCgMIoxAACFUYwBACiMYgwAQGEUYwAACqMYAwBQWKRkrpK/HI1F3CfpjvO+fFLS/YfWCUiMeSmMexmM++FjzC/sCSmlU+d/8VCL8X4i4j0ppZuLduIaw5iXwbiXwbgfPsb84vFragAACqMYAwBQ2COhGL+2dAeuQYx5GYx7GYz74WPML1Lx14wBALjWPRLujAEAuKZRjAEAKKxYMY6IF0TERyPiYxHxilL9uNpFxC9HxL0R8cE9XzsREW+PiD9Z/H+8ZB+vNhHx+Ih4R0R8OCI+FBHfvfg6434FRcQoIt4dEf9jMe6vWnz9SRHxrsW15g0RMSjd16tRRNQR8b6IeOvic8b9IhQpxhFRS/qXkr5G0tMlvSQinl6iL9eA10l6wXlfe4WkW1NKT5F06+JzXD5zSd+XUnq6pGdL+juL+c24X1kTSc9PKX2RpGdIekFEPFvST0j6qZTS50k6I+k7y3Xxqvbdkm7b8znjfhFK3Rk/S9LHUkqfSClNJb1e0gsL9eWqllL6PUmnz/vyCyXdsvj4FkkvOsw+Xe1SSnenlP5o8fGGdi9QN4pxv6LSrs3Fp+3iX5L0fElvWnydcb8CIuImSV8n6RcXn4cY94tSqhjfKOnOPZ9/avE1HI7rU0p3Lz7+jKTrS3bmahYRT5T0TEnvEuN+xS1+Vfp+SfdKerukj0s6m1KaLzbhWnNl/LSk75fULz6/Toz7ReENXNe4tLu2jfVtV0BErEr6d5K+J6W0vvd7jPuVkVLqUkrPkHSTdn8D99SyPbr6RcTXS7o3pfTe0n15NGsKtftpSY/f8/lNi6/hcNwTETeklO6OiBu0exeByygiWu0W4l9NKf37xZcZ90OSUjobEe+Q9KWSjkVEs7hL41pz+T1H0jdGxNdKGkk6IulnxLhflFJ3xn8o6SmLd9sNJL1Y0lsK9eVa9BZJL1t8/DJJby7Yl6vO4vWyX5J0W0rp1Xu+xbhfQRFxKiKOLT5ekvRV2n29/h2SvmmxGeN+maWU/lFK6aaU0hO1ey3/nZTSS8W4X5RiCVyLZ1E/LamW9MsppR8t0pGrXET8uqTnafdPmt0j6Yck/aakN0r6HO3+SctvTimd/yYvXKKI+EuS/oukP9b/eg3tB7T7ujHjfoVExBdq941CtXZvNN6YUvqRiPhc7b5J9ISk90n61pTSpFxPr14R8TxJfz+l9PWM+8UhDhMAgMJ4AxcAAIVRjAEAKIxiDABAYRRjAAAKoxgDAFAYxRgAgMIoxgAAFPb/A28KXHg8kxiGAAAAAElFTkSuQmCC\n",
"text/plain": [
"