notebooks/tutorials/intro.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Author: Ismael Mendoza*\n",
    "\n",
    "*Last updated: 05/22/23*\n",
    "\n",
    "# MultiCAM intro notebook \n",
    "\n",
    "In this tutorial, we will show you the basics of using MultiCAM for connecting the mass accretion history (MAH) of dark matter halos to their present-day properties.\n",
    "\n",
    "We will use pre-processed data of $10,000$ halos in a narrow mass bin of $[10^{12}, 10^{12.2}]h^{-1}, M_{\\odot}$ from the ROCKSTAR catalog of the Bolshoi n-body dark matter simulation suite. We also pre-processed their merger tree information to extract the main line progenitor branch of each halo and compute the MAH. The data is available in the `data/` directory of the MultiCAM repository at: https://github.com/ismael-mendoza/multicam/tree/main/data/processed/bolshoi_m12. \n",
    "\n",
    "The MultiCAM paper has more details about the data and the algorithms used in this tutorial and is available at: https://arxiv.org/abs/2302.01346. If you use the MultiCAM algorithm in your research, please follow the instructions in the README to cite: https://github.com/ismael-mendoza/multicam#how-to-cite."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Installation\n",
    "\n",
    "Follow these steps to ensure you can run this notebook. First, pip install required packages. Then, install the MultiCAM package from the github repository.\n",
    "\n",
    "```bash\n",
    "pip install -r https://raw.githubusercontent.com/ismael-mendoza/multicam/main/requirements.txt\n",
    "pip install git+https://github.com/ismael-mendoza/multicam.git\n",
    "```"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data can be downloaded by running: \n",
    "\n",
    "```bash\n",
    "wget https://raw.githubusercontent.com/ismael-mendoza/multicam/main/data/processed/bolshoi_m12/final_table.csv\n",
    "wget https://raw.githubusercontent.com/ismael-mendoza/multicam/main/data/processed/bolshoi_m12/z_map.json\n",
    "```\n",
    "\n",
    "or you can directly navigate to the `data/processed/bolshoi_m12/` directory at this [link](https://github.com/ismael-mendoza/multicam/tree/main/data/processed/bolshoi_m12) and download both files manually.\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note:** If you are working in google colab, you can directly run these commands in a code cell by appending `!` at the beginning of each line. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Importing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load and understand data\n",
    "\n",
    "First we load the data that will be using to produce our results. The MultiCAM repository has a set of functions to load the data from the pre-processed files.\n",
    "\n",
    "Feel free to ask me if you have questions about the details, but essentially the data is loaded into an astopy table from csv files, and cuts are applied corresponding to a resolution scale of 50 particles per halo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from multicam.mah import get_mah\n",
    "mah_data = get_mah('.', cutoff_missing=0.05, cutoff_particle=0.05)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start by understanding the data. First let's extract the MAH of a single halo and plot it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10000, 165)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# extract mah of all haloes \n",
    "ma = mah_data['ma_peak']\n",
    "ma.shape # the first dimension is the number of haloes, the second is the number of snapshots."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first dim corresponds to the number of haloes and the second one corresponds to the number of snapshots. The snapshots are ordered from the earliest to the latest. Just like in the MultiCAM paper, we use the definition of $m(a)$ corresponding to $M_{\\rm peak}$ which means that our MAH definition always monotonically increases. \n",
    "\n",
    "Next, we extract the scales corresponding to the snapshots used. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(165,)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scales = mah_data['scales']\n",
    "scales.shape"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again the dimension corrsponds to the number of snapshots we decided to use. The snapshots don't encompass all of the simulation snapshots, but rather we make a cut on the earliest snapshot consider based on particle resolution. \n",
    "\n",
    "Now we can plot the MAH of multiple haloes. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj0AAAG2CAYAAACUDjeHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABpTElEQVR4nO3dd3hUddrG8e9MeiAJgZBCCITeewepIoiKoqsiNsS2SrGwr4ACYlnFilhQXATZdXVxEdcGAhpARJr0XkOHBEIgvc6c949jgkiR1DOTuT/XlSszZ86ZeQ4tN79qMwzDQERERKSCs1tdgIiIiEh5UOgRERERj6DQIyIiIh5BoUdEREQ8gkKPiIiIeASFHhEREfEICj0iIiLiEbytLsBVOJ1Ojh8/TlBQEDabzepyRERE5AoYhkFaWho1atTAbr98W45Cz2+OHz9OTEyM1WWIiIhIMRw5coSaNWte9hyFnt8EBQUB5i9acHCwxdWIiIjIlUhNTSUmJqbw5/jlKPT8pqBLKzg4WKFHRETEzVzJ0BQNZBYRERGPoNAjIiIiHkGhR0RERDyCQo+IiIh4BIUeERER8QgKPSIiIuIRFHpERETEIyj0iIiIiEdQ6BERERGPoNAjIiIiHsElQ8/y5csZOHAgNWrUwGaz8dVXX/3pNcuWLaNt27b4+flRv359Zs+eXeZ1ioiIiPtwydCTkZFBq1atmDZt2hWdf+DAAa6//np69+7Npk2beOKJJ3jwwQdZtGhRGVcqIiIi7sIlNxwdMGAAAwYMuOLzp0+fTp06dXjzzTcBaNKkCStWrOCtt96if//+ZVWmiIh4AIfTwcnMkxgYF7yWnefgTGauBVW5EUc2XllnAAj0r0zD2q0tK8UlQ09RrVq1ir59+553rH///jzxxBOXvCYnJ4ecnJzC56mpqWVVnoiIuKHdybuZHz+f+QfmczLzpNXlVAgNc+zMe3izZZ9fIUJPQkICERER5x2LiIggNTWVrKwsAgICLrhm8uTJPP/88+VVooiIuImVx1fyzoZ32H56e+Exb5s3XnYv8h0G+U6nhdW5BlvhdwM7zsJjNozfvoDzWsbMI16FV1qjQoSe4nj66acZPXp04fPU1FRiYmIsrEhERKyS48hh7u65fLP/G3Ym7wTAx+5Dz5o9uaHuDXSv2Z1PVh3jxe92AFCraiChgT7cf1UdBrasgd1u7Q/zcpGXBdv/B7/OhGPr/vz8KrWh3X3Q+AYIawA263+NKkToiYyMJDEx8bxjiYmJBAcHX7SVB8DPzw8/P7/yKE9ERFzY/rP7GbN8DHvO7AHMsDO40WAeavkQoX6hHE/J5pkvd/LF+qMAjL6mIY9d3cDKkstHXhYcXgWHV8Pp/bA/Dn4bm1MoIBTq9ISQmmD3gpAYqFoHQn/7srvWfKkKEXq6dOnCggULzjv2ww8/0KVLF4sqEhERV3c66zQfbvmQubvnkm/kU9W/Kn9t+VcG1BlAqH8oX244yssLNpCUbo7/tNvgib4NGdWnvsWVlxHDgJM7YP8S8+vQSsjPPv+ckFrQ/j5ocw9Uqm4ec4EWnCvlkqEnPT2dffv2FT4/cOAAmzZtomrVqtSqVYunn36aY8eO8a9//QuARx55hPfee48xY8Zw//33s2TJEv773/8yf/58q25BRERcVI4jh1nbZjF722wy8zMB6FmzJ891fY6th5w8/tleUjJz2Xw0BQBvu40mUcFMuL4JnepWs7L0suF0wp6F8PObF3ZbBdWAuj0hvAlEtjBbdexe1tRZClwy9Kxbt47evXsXPi8YezN06FBmz57NiRMnOHz4cOHrderUYf78+Tz55JO8/fbb1KxZk48++kjT1UVE5Dx/7MpqVq0ZT7Z7kk5Rndh05Cx//WQVuQ5zYK6X3cYTVzfg4Z518fN23x/0hZIPQOJ2MJyQnggpR8FwwL44s4UHwMsP6vSAen3Mr+qN3Kol58/YDMO4cOEBD5SamkpISAgpKSkEBwdbXY6IiJSyE+kn+Mu3fyEtN42q/lV5uuPT9I/tj81m43R6Dje8u4ITKdn0alSdv7StSZOoYOqHV7a67JJx5JutOOtmml1Wl+IbBB0fhM7DoXJ4+dVXCory89slW3pERERKk9NwMv6X8aTlptG0WlOmXT2NsICwwtffWLybEynZ1K1eiXeHtCHI38fCaovJ6TDH4GSnwsGfzZCzLw4yCtYYskFUK/D2g8AwqBIDXj4QFAWt7zQHJVdwCj0iIlKhpeSk8PaGt/k14VcCvAN4vcfr5wWe3QlpfP7rEQBe+0tL9wo8Tgcc2wBbPofNcyA37cJzAqtBm7uh3TBzZpUHU+gREZEKa0PiBkYtGUVqrrnqfk3nHbzy7SngVOE5uxLScBowoHkk7WOrWlTpJWSnQvJ+cybV0V/BkWceNwxzXE7SXshJufC6yJbnxuXU6my27ohCj4iIVCwOp4MjaUc4nnGcvy37G+l56ThzIslOHMCGjEZAwgXX+HjZGHtt4/Iv9mIMwww4Pz4Ph1b8+fn+IVDvamh7L8R0BLu3Qs4lKPSIiEiF8uLqF5m3d17h8/zMWLIOP8D1zWvRqe7FW3KaR4cQG1apvEq8uAPLYelkOL4R8rPOHQ8Mg6iW5nRx/5DfHa8KVetC9SbgpR/nV0K/SiIiUmGsPbG2MPA48wNxZtXCkTiEkb2a8Ld+DbG5yvRrpwNObIaTO+HULohfBglbzr1u94GWg6HXOHPAsZQKhR4REakQdiUkM3zReLBBbnJnchIHMfqahtx7f22qBPpaV5jTCWknzC0dDv0Cp/fBiS2Qffb88+ze0P5+6PgwhMaaM6ukVCn0iIiI21t3KIn75z+JUSkBZ35lfFKv5/UhbbixVQ1rCjp7xJxRtX+JObvq991VBfxCoEZrs4uqVheo19vt1shxNwo9IiLitpLSc3j1++18e/w9fKpsAsOLCR1f4Nb7rsbHq4w3u8xIgi3/hZQj4Mw3VzhOTzQfJ2w1Vz4uYPMyt3Ko2wsimkNYQ3PNHI3FKVf61RYREbe0/1Q6d89eQErQLHyqHANsvNjtJQY1uKZ0Pygn3QwxZw6Yj/MyzSnk8UvBkXvp62K7Q7Obze9V66i7ygUo9IiIiNvZePgMw/79HXnhH+Dlk0Yl72Be7v4ifWr1Kfqb5eeaQcbpgNSjkH7SbKVJ3Ab7l8Lh1eDMu/i1Ndqc24QzKAqCa5itOmENoFq9kt2klDqFHhERcStxOxMZMfd7vKI/xO6dTt3g+nzY7wMiK0Ve/sK0BNj0qbnxZkHAST5ofv99V9TFBEdDtfrmVg12L7Nrqn5fiGhWavclZU+hR0RE3EJuvpM3F+9mxupVBNSegd07nfpVGjCr/0xC/S+yb9TZI3BsPZzea66Bc2ilOd7mcgLDzNYauxcE1TAHF9frYw42dpXp7lJsCj0iIuLy4k+l8+h/F3Iw73sC62zAZs+nUWgjPur3EVX8q0BeFmSehvwcOLIGtn0J+34EjPPfKKaz2UJj94KgSAitY463CagKNrsGFldw+t0VERGXdeh0Bu8v38h3x9/FXnkHBavttI9oz9RuLxGy/WvYNs9cA+dig4prtDG7pWp2hPpXa5yNh1PoERERlxN/Kp3XF+1m8Z6t+Md8jL1yMhg2ukVexcM1r6LNoY3Y3ut8/mabXr5ma014UzPgtBqikCPnUegRERGXErczkSfmbCLbto/KtWfj9M6mBv5My/Wl/tr/wupPz50cWgfaDYXGN5gtOhp3I5eh0CMiIi7BcDqY891CDq79lpGBG/koKp1su52mOTm8n3CUas7fZlgFVIXYq6DdfVC3N9jLeBFCqTAUekREpPwYBiTtMXcSP7bB3HQzPREj+yxGVgpDcPJzkD9PhFcn126na3YeU/wbUqnLEHN8To02UKWWWnSkWBR6RESkbGWchn0/mPtQ7V8KGScvOMX229eigCCejgglzwZXR3Tk9avfx8fHr9xLlopJoUdEREqfYZgrGa+bCTu+Pn9mlXeAudFmjTYYUa34106DTzancqbyCRxRX+Mw8rmm9jW82uNVfOzaukFKj0KPiIiUHkc+bP8SVkyFk9vPHY9oDg2uMRf6i+kE3mbrzftL9/H6pt3YfHOoEvMdDiOffrX78UqPVxR4pNQp9IiISPGknoBTu849T94PK9+FMwfN594B0OJWaH8/RLe94PIv1h/l9UW7AQf1mn5DYk4OnaI68WqPV/G268eTlD79qRIRkSvndJq7i6+bBbu/B8Nx4TmB1aDzcOjwgLlX1UWs3J/E019uAaBb+81sydhLkG8Qf+/2dwUeKTP6kyUiIheXlw07v4Fd30FGEmSdhfREyEw6d061BoVdVXj7QYvboO1Q8A285NvuO5nOI5+sJ89h0KN5Jlsy5wEwsfPEP980VKQEFHpEROR8SXthw79g478hK/nC1/2CzdWO298P4Y2L9Nan03O4f/avpGbn07pWIKcDP8CR5mBAnQEMqDOglG5A5OIUekREBM4ehq1fmIOQE7aeOx5cE9rcDdUbQUAV8K8C1RtftiXncl5duIvDyZnEVA2ga7vNfLLrEOGB4YzvNL5UbkPkchR6REQ8mWHAmumweCI488xjdm9zllW7YdCgX6ntPJ6Slcc3m48D8PTACCau+8R83PFpQvxCSuUzRC5HoUdExFNlnYGvR5pjdgBqdYFWd0CTGyGwaql/3NebjpGd56RRRBCLE2aQ68ylU1Qnrq51dal/lsjFKPSIiHiio+vhi/vMbi0vX+j/MnR4sMy2dzAMg8/WHAagS/Mk5h2Ow8vmxbgO47BpSwkpJwo9IiKe5I/dWaGxcNtsc0+rMrTpyFl2JaTh522wPm02AEMaD6F+aP0y/VyR31PoERHxFFln4esR57qzmtwIN70H/mU/nmbO2iMANG+yjT1pBwj1C+XR1o+W+eeK/J5Cj4hIRefIgw3/hGWvQMYpsPuY3VkdHyqX3crTss0BzDavdI4YXwHweNvHCfYNLvPPFvk9hR4RkYrK6TQXF1zyIpzeZx6rVh9umXHRbSHKytebjpOV56B6nTiyHBk0rdaUQfUHldvnixRQ6BERqWgc+bBtHqyYcm5vrMAw6DUO2t0HXuW3kWfBAGa7/1Gy/VcD5hR1L7tXudUgUkChR0SkokjaC9u+hE2fwtlD5jG/YOj0CHQdBf7l3520Yl8SO06cpVLst4DBwLoDaR3eutzrEAGFHhER9+Z0miFn7Yfnr6RcsOlnx4fKZaDypby/dD/ewZuwBxwiwDuAJ9o9YVktIgo9IiLuKmEbfPckHF1rPrd7Q93e0PwWaHoT+FaytLwNh8+wKj6JSvUXA/DXln8lPDDc0prEsyn0iIi4i+QDcGo3nDlgBp7N/wHDAb6VoecYaHNPmaykXFzvL92P3e8Edp+zBHgHcHfTu60uSTycQo+IiKtz5MHCcfDrRxe+1vQm6D8ZQqLLv67L2J2Qxo87E/GttgeATpGd8PPys7gq8XQKPSIirizjNMwdCgd/BmwQ2RxC60DVOuamoHV7WV3hRU3/aT8A1cMPkAp0i+5mbUEiKPSIiLiuhG0wZ4i5P5ZvZfjLR9BogNVV/aljZ7PM3dTt2WTY9oEBV0VfZXVZIgo9IiIux+mAle/C0pfAkWu27Az5D4Q3sbqyK7JgywkcToPG9RI4ZjiIDY6lZlBNq8sSUegREXEpp/fDV4/CkTXm84bXwqAPXGqA8p9ZvCMBgKph8RxLVSuPuA6FHhERV2AY5kDlH56FvEzwDYJrJ0Obu8tlf6zSkpSew7pDZwCDE7mbAIUecR0KPSIiVks5Cl+PhPil5vPY7jDofahSy9q6iiFuZyKGAY1iMjiefRI/Lz/aRbSzuiwRQKFHRMQ6+Tlm686yVyEnBbz9oe/z0PFhsNutrq5YFm9PBCAm+jDHU6BDZAf8vf0trkrEpNAjIlLenE7Y/iXEPW/OzAKIbg83T4ewBtbWVgIZOfn8vC/JfOy1HVDXlrgWhR4RkfJ0YgvMHw1HfzWfB0VB72eg1Z3g5d7/JC/fc4rcfCe1wuzsOrsZUOgR1+Lef8NERNxFThosfRnWTAfDaa67c9WT5qagvoFWV1cqFu8wu7Za1UtlWWo+0ZWjqRXkfuOSpOJS6BERKWv74uCbUZB6zHze7GZz64jgKGvrKkV5DidxO83QUycqj2WpEBsci82NZp5JxafQIyJSVnIzYPFEWDfTfB4aC9e/CfX7WlpWWVh7IJnU7HyqVfIlwD8dgIhKERZXJXI+hR4RkbJwdD3Me8DcER3MGVl9nwPfSpaWVVZ+/K2V5+om4ZzMWglARKBCj7gWhR4RkdJkGOa4ncUTwZkHwTXhpvegXm+rKytTa+KTAejeoDrfnzIDkEKPuBqFHhGR0pKZbI7d2fWd+bzJjXDjuxBQxdKyylpqdh47E1IB6FinKh8fMrehiKwUaWVZIhdQ6BERKan0k7D6fVj7EeSmgZcv9HsJOj7kVltIFNf6Q2cwDKhVNZCIYH8SM9TSI65JoUdEpLhyM8xp6L9+BPnZ5rHwZmZ3VnRba2srR78eMLu2OsRWJTMvk9Rcs9VHA5nF1Sj0iIgUx+E18NUjkBxvPo9uB93/z9wV3U23kCiuXw+aoadjnVBOZp4EINA7kMo+la0sS+QCCj0iIkWRlgCrpsGq98xFBoOj4Yap0OAaj+jK+qPsPAebj6QAZktPYqa5/UREpQit0SMux2X/OzJt2jRiY2Px9/enU6dOrF279rLnT506lUaNGhEQEEBMTAxPPvkk2dnZ5VStiFRo2SmwbhbMvgHebAwr3zEDT6sh8OhKaNjPIwMPwJajKeQ6nIRV9qNOWCUSMn4bxByoQczielyypefzzz9n9OjRTJ8+nU6dOjF16lT69+/P7t27CQ8Pv+D8zz77jHHjxjFr1iy6du3Knj17uO+++7DZbEyZMsWCOxCRCiH9FKz5ANbOgJzUc8ej20P30dD4eutqcxG/79qy2WwkZv42iFnjecQFuWTomTJlCg899BDDhg0DYPr06cyfP59Zs2Yxbty4C85fuXIl3bp148477wQgNjaWIUOGsGbNmkt+Rk5ODjk5OYXPU1NTL3muiHiYk7vMwckb/w35WeaxsEbQ+k5zC4nQ2tbW50LW/m4QM6CZW+LSXK57Kzc3l/Xr19O377ll2u12O3379mXVqlUXvaZr166sX7++sAssPj6eBQsWcN11113ycyZPnkxISEjhV0xMTOneiIi4l6yzsGUufHwdvN8Jfp1hBp4abeGOz2D4arjqCQWe33E4DTYcOgP8LvSopUdcmMu19CQlJeFwOIiIOP8vTEREBLt27broNXfeeSdJSUlcddVVGIZBfn4+jzzyCM8888wlP+fpp59m9OjRhc9TU1MVfEQ8zandsG0e7F8Cx9ab43QAbF7QaIC5zk6dnh47XufP7DyRSlpOPkF+3jSJCgbOhR6N6RFX5HKhpziWLVvGyy+/zPvvv0+nTp3Yt28fjz/+OC+++CITJ0686DV+fn74+fmVc6UiYrn8XNj1Lfw6Cw6tOP+1sIbQ7BZoey+ERFtTnxspGM/TtnYoXnYzGBYMZFZLj7gilws9YWFheHl5kZiYeN7xxMREIiMv/j+HiRMncs899/Dggw8C0KJFCzIyMnj44YcZP348dg9bM0NE/sDpgIM/w7YvYee3kGX+sMZmN9fVaXSduTdWSE1r63Qz5wYxm11b2fnZnM05C2hMj7gmlws9vr6+tGvXjri4OAYNGgSA0+kkLi6OkSNHXvSazMzMC4KNl5cXAIZhlGm9IuLCDAO2fA4/PgdpJ84drxwJ7YZC26Fq0SmBTYfPAtC+dihA4cKEAd4BBPsGW1WWyCW5XOgBGD16NEOHDqV9+/Z07NiRqVOnkpGRUTib69577yU6OprJkycDMHDgQKZMmUKbNm0Ku7cmTpzIwIEDC8OPiHiYU7th/t/MFh6AgFBzA9Dmt0Dtq8DLJf/5cxtOp0FimjkDNjasEvC7QcyBWphQXJNL/q0fPHgwp06d4tlnnyUhIYHWrVuzcOHCwsHNhw8fPq9lZ8KECdhsNiZMmMCxY8eoXr06AwcO5KWXXrLqFkTEKpnJ8PObsOZDcOaBdwD0HANdRoC3xvGVljOZuTicBjYbVK3kC2g8j7g+lww9ACNHjrxkd9ayZcvOe+7t7c2kSZOYNGlSOVQmIi4pLxvWfmgGnmxzWwQaDoABr2qaeRk4lW628oQG+uLjZf4n9PctPSKuyGVDj4jIFTl7GNbPhg2fQIY5poTwZnDNC9Cg72UvleI79VvXVvXK51rPClt6FHrERSn0iIh7Sj4Ai8bDnu/Pra8TXBP6jIeWg8Gu8XxlKem3lp7qQedCT+EaPZW0Ro+4JoUeEXEvhgHrP4ZFEyAvwzxWpwe0f8DcC8vLx9r6PERBS09YZd/CYwVbUCj0iKtS6BER95F8ABb8H+z70Xxeuxvc8BZUb2RtXR6osHvrIi096t4SV6XQIyKuL3EHrHgLtn1hdmV5+UHfSdDpUdDio5ZISs8FzoWeHEcOydnmYoUKPeKqFHpExHU58uG7x83dzgvU7wv9X1brjsXOdW+ZoadgYUI/Lz9C/EIsq0vkchR6RMQ1OfLhfw+bG4JigyYDofvfoEZrqysTLuze+v14Hi1MKK5KoUdEXI/TAV+PMAOP3QcGf2Luei4u44+ztzSeR9yBOsNFxLXk58I3j8GWOWDzgts+VuBxMXkOJ8mZ5piegu4trdEj7kAtPSLiOg7+AvNHw6ld5g7of/nI7NYSl5KckYthgJfdRmigOWX9SNoRQFtQiGtT6BER62Umw+IJsOlT83lgGAycqsDjogrG81Sr5IuX3YbTcPLzUXNj17bhba0sTeSyFHpExFqH18AXwyD1GGCDdveZ09EDQq2uTC6hYN+tgq6tHad3cDLrJIHegXSK6mRlaSKXpdAjItZwOmHVu/Dj82A4oFp9GDQdYjpYXZn8iT/O3FpyeAkA3aK74evle8nrRKym0CMi5S83A758GHZ9Zz5vfqvZneUXZGlZcmX+OHNr2dFlAPSO6W1VSSJXRKFHRMpXRhJ8NhiOrTNXVh7wqtmlpbVd3MbvFyY8mnaUvWf24mXzokfNHhZXJnJ5Cj0iUn6S4+HffzG/B4TCkM+hlsaAuJvfd28tPbIUgLYRbbUSs7g8hR4RKXuGYS40+P0YyDwNVWrB3V9CWAOrK5Ni+H331le/hR51bYk7UOgRkbKVtM9ce+fAT+bzqFZw51wI0nou7qqgpSfAL5sNiRsAhR5xDwo9IlJ2tv/PHLDsyDXH7/R4Cro9Bt5+VlcmJVAQeg5lrcdhOGgQ2oCaQTUtrkrkzyn0iEjZOLUHvhphBp56feD6N6FqXaurkhLKyXeQmp0PwObkXwC18oj7UOgRkdKXmwlzh0JeBsR2h7u+ALuX1VVJKUhKN/fc8vHOY03CSgD6xPSxsiSRK6YNR0WkdGUkmYHn5A6oFA5/manAU4EUdG1VqRZPVn4W0ZWjaVqtqcVViVwZtfSISOnZ+R18+zhkJoHdx9wwVAOWK5Sk30KPV9BmAPrF9sOmNZbETSj0iEjJ5aTDgqdg82fm8/CmcPN0c6aWVCin0nPAlkuW9zYA+sf2t7gikSun0CMiJZO4Hf47FE7vBZsduj4GvZ/RDK0K6lRaDt6Vd+K05RITFEPTquraEveh0CMixWMYsPHfZgtPfhYE1YBbZ0LtrlZXJmUoKT0H7+AtgNnKo64tcScKPSJSPEtfguWvm4/r94WbP4RKYdbWJGXuROpZvCvvBuDa2GstrkakaBR6RKTofnrtXODpPQG6/w3smgzqCQ5l/YrNL58wv5o0DG1odTkiRaLQIyJF8/MUs5UHoN9L0HWktfVIuUpiLQBdIvqoa0vcjkKPiPw5pwMOroDNc87N0Lp6kgKPh0nLTSPXdyc2oJ9mbYkbUugRkUvLzTC7sjZ9Bhknzx3v9TR0H21dXWKJeXu+wmZz4MiOoH0NzdoS96PQIyIXd2Qt/O+vkBxvPg8IhSYDoeVgiL3K2tqk3DkNJ//Z9R/zSWoXKvlqlW1xPwo9InJOfq65fcS2ebDqPTCc5lT0Aa9Ao+vAy8fqCsUiK46t4HjGEQyHP9XopvE84pYUekQ8ldMBp3bD8Y1wfIP5PWEbOHLOndPyDjPwBIRaV6e4hH/v+DcAeWc70DZGW4uIe1LoEfFEuxbAwrFw9vCFr/mHQI020OEhaHJD+dcmLmf/2f2sOrEKDBu5Z7rQo2d1q0sSKRaFHhFPcvYwfD8Wdi8wn/tUghqtzZBT8FW1LqjrQn7ns53mjL289CYYeVXp3lCLUIp7UugR8RRbv4BvRkFeJti9oeso6PEU+FayujJxYSk5KXwb/y0AecndaBIVTHiQv8VViRSPQo9IRed0wJIXYcVb5vNaXeGGKRDexNq6xOXlO/N5Y90bZOVnUdkWQ1pmXXp0UCuPuC+FHpGKLDsFvnwY9iw0n3d7Aq5+FuyabiyXl5mXyZjlY/jp6E/YsJFzqi9go2cDjecR96XQI1IROfJg/WxY9gpkJoG3P9z4LrS83erKxA0kZSUxMm4k209vx8/LjxHNJvL8594E+HjRLlYz+cR9KfSIuLvcDNg6FzKTzedOB2z+DyTvN59XrQd/mQHR7ayrUdxGZl4mDy56kP0p+wn1C+WdPu+wemcQsIvOdavi561WQnFfCj0i7sowzEUEf3gWUo9d+HpgGPQaB+3u06KCcsVe+/U19qfsp3pAdWZfO5tawbV4/ZvVAPRoqK4tcW8KPSLu6MRm+H4cHF5pPg+pBXV6nHu9Wj3o+BD4BVlTn7ilRQcXMW/vPGzYmNx9MrWCa5GZm8+vB84ACj3i/hR6RNxJxmlzJtb62YAB3gHmxp9dR4FPgNXViRs7nn6c51c+D8CDLR6kU1QnANbEJ5PrcBJdJYC6YVreQNybQo+IqzMMSNoDO76BVe+aM7IAmv8FrnkBQmpaW5+4vTxHHmOXjyUtL42W1VvyaOtHC1/7/NcjgNnKo/22xN0p9Ii4guxUSNwOidvMXc0Np3k8PwcOLD83KBkgogUMeBViu1lTq1Qo+c58xv48lk2nNlHZpzKvdn8VH7s5BuynPadYuD0BL7uNoV1rW1ypSMkp9IiUp+xUWDcLfp0Jmad/O2iYqyRfjt3HHLPT/C/Q6g6tsyOlwmk4mfjLRH449AM+dh/e7PkmNYPMlsOcfAfPfbMdgPu6xtI4MtjKUkVKhUKPSHnIOA1rpsPaD891T/1RcDRENIfqDcHLzzxms0FEM6h3Nfjrh46UHsMweHH1i3wX/x1eNi/e6PkGXaO7Fr7+0c8HOJCUQfUgP57o28DCSkVKj0KPSFlKPQ4r34P1H59rzQlrCFeNhtpdzp3nFwyBVa2pUTzSOxvf4Ys9X2C32Xml+yv0qdWn8LWjZzJ5d8leACZc34Qgfy15IBWDQo9IWchJhyV/h3UzwZFrHotsCT3+DxoPBLvd2vrEo32972s+2voRAM91eY5r61xb+Fq+w8mEr7aRneekU52q3NiqhlVlipQ6hR6R0nZwBXw1HM4eMp/X6gLd/w/qX212V4lYaH3iep5b9RwAD7V4iJsb3Fz4mmEYjP/fNpbtPoWvl50XBzXXjC2pUBR6REqLIx8WjzfH7gAE14SBb0ODvtbWJfKbI6lHeGLpE+Q787mm9jWMbDPyvNdf+X4Xn687gt0G7wxpTcMILW4pFYva2EVKy5IXzwWetkNh+CoFHnEZy48uZ+jCoZzNOUvTak156aqXsNvO/Qj4YNl+PlweD8Art7Tk2uZRVpUqUmbU0iNSGnYtgF+mmo9v/ge0GmxpOSIF0nLTeO3X1/hq31cA1Ampw7t93iXA21zB2+k0eG3Rbqb/ZK4FNf66JtzeIcaqckXKlEKPSEklx8P/HjEfdx6uwCMu40z2GYbMH8Kx9GPYsHFP03sY1WYU/t7+AGTm5vPk55tYtD0RgNHXNOShHnWtLFmkTCn0iJREXhZ8fi/kpEBMJ3NbCBEX8c7GdziWfoyoSlG80v0V2ka0LXztQFIGo/6zgW3HUvH1svParS0Z1CbawmpFyp5Cj0hxOZ3mLK3ErRAYBrfNBi+tZyKuYXvSdubtmQdwXuBJSs/hnbi9fLbmMPlOg2qVfPnHve1oV1vrREnFp9AjUhyGAQvHwfYvzS0ibvsYgrWeibgGp+Hk5TUvY2BwQ90baBvRlnyHk1m/HOCduH2k5+QDcHXjcJ67sRkxVQMtrlikfLjs7K1p06YRGxuLv78/nTp1Yu3atZc9/+zZs4wYMYKoqCj8/Pxo2LAhCxYsKKdqxeOsmGJuKQFw83RzXywRF/H1vq/ZkrSFSj6VGN1uNHsT0/jL9FW8vGAX6Tn5tIgO4bOHOjHzvg4KPOJRXLKl5/PPP2f06NFMnz6dTp06MXXqVPr378/u3bsJDw+/4Pzc3FyuueYawsPD+eKLL4iOjubQoUNUqVKl/IuXim/jpxD329ida1+BFrdaW4/I76TmpjJ1w1QAHm31KF/+msqbi9eR63AS5O/NxOubcmu7mtjtWnRQPI9Lhp4pU6bw0EMPMWzYMACmT5/O/PnzmTVrFuPGjbvg/FmzZpGcnMzKlSvx8THHVMTGxl72M3JycsjJySl8npqaWno3IBXX4dXw7ePm425PQOdHLS1H5PdOZ53msSWPkZydTN2QumQmdebV73cB0LtRdSbf0pLIEH+LqxSxjst1b+Xm5rJ+/Xr69j23qJvdbqdv376sWrXqotd88803dOnShREjRhAREUHz5s15+eWXcTgcl/ycyZMnExISUvgVE6N1KeRPpB6Hz+8BZx40HQR9n7O6IpFC8Snx3LXgLrYkbSHEL4RB0aN5beE+AMZe25hZ93VQ4BGP53KhJykpCYfDQURExHnHIyIiSEhIuOg18fHxfPHFFzgcDhYsWMDEiRN58803+fvf/37Jz3n66adJSUkp/Dpy5Eip3odUMPk5ZuDJOAnhzWDQ+9pHS1zG6hOruXvB3RxLP0ZMUAwvdviQN77NwjBgSMdaPNKzrvbQEsFFu7eKyul0Eh4ezj/+8Q+8vLxo164dx44d4/XXX2fSpEkXvcbPzw8/P79yrlTckmHA/L/BsXXgXwXu+BR8K1ldlQhZ+VlMXT+Vz3Z9BkDr6q15qvVkRvx7Dxm5DjrVqcrzNzZT4BH5jcuFnrCwMLy8vEhMTDzveGJiIpGRkRe9JioqCh8fH7y8vAqPNWnShISEBHJzc/H19S3TmqWC2/iJ+WWzw62zoGodqysSYdPJTYxfMZ7DaYcBuKbmTeSdHMjN723B4TSoVTWQ6Xe3w9fb5Rr0RSzjcn8bfH19adeuHXFxcYXHnE4ncXFxdOnS5aLXdOvWjX379uF0OguP7dmzh6ioKAUeKZmErbDgKfNxnwlQ/2pr6xGPZxgG/9z+T4YuHMrhtMME+4RR1/EEX/7QhW83J+FwGvRoWJ1/3d+R0Er690/k91yupQdg9OjRDB06lPbt29OxY0emTp1KRkZG4Wyue++9l+joaCZPngzAo48+ynvvvcfjjz/OqFGj2Lt3Ly+//DKPPfaYlbch7i47Ff47FPKzoUE/6Pak1RWJh8vKz+K5lc+x4MBva5Clt+XYsYEccwZgs8GA5pEM71Wf5tEh1hYq4qJcMvQMHjyYU6dO8eyzz5KQkEDr1q1ZuHBh4eDmw4cPY7efa6SKiYlh0aJFPPnkk7Rs2ZLo6Ggef/xxxo4da9UtiLszDPhmFCTvh5AYuPlDsLtcw6h4kN3Ju3lmxXj2nNmNYdjJSbyevDNdia4SyF/aRvOXdjWpXU1jzUQux2YYhmF1Ea4gNTWVkJAQUlJSCA4OtrocsdrPb5oLENp94P6FULO91RWJh0rISOC9je/xzf5vMDBw5lci+9hdNA1tw//1b0T3+mFaaFA8WlF+frtkS4+IpVZPP7ficv+XFHikXB1MOcjWpK3Ep8QTfzaeFcd+IddpLqSal9KSoMybeeH6ztzcJlphR6SIShx6HA4Ha9asYcOGDSQmJnLmzBlCQ0OJiIigXbt2dOzY8bxZVSIubf1sWPhbt2iPMdDpr5aWI54hx5HDD4d+YO7uuWw4ueGC1/MzY3Ek3cCDHXoysk99An31/1WR4ij235wVK1Ywbdo05s+fT0ZGRuFxwzDOWxOicuXKXH/99YwYMYJu3bqVrFqRsrRlLnz7hPm46yjo/Yyl5Yhn2HF6B4/++CjJ2ckA2LBjZNUmNzscZ044zuxo2ke14aW/tqB+eJDF1Yq4tyKP6fn555958skn2bhxI4ZhYLfbad68Oc2aNaNatWoEBweTkpLC6dOn2bZtGzt27MDpdGKz2Wjbti1Tpkyhe/fuZXU/xaYxPR7uzEGY1hnys6DDQ3Dd61pxWcqc03AyZP4QdpzeQZB3NdKT2pN2qi1Gfghhlf3o2ySc/s0i6dWouhYYFLmEMhvTc8cddzB37ly8vb258cYbue++++jTpw9BQZf+30dqaipxcXHMnj2bhQsX0qtXL26//Xb+85//FOWjRcqOYcB3T5qBJ7Y7DHhNgUfKxTf7vmXH6R3YDX9O7HgEwxFE06hgnh3YlI6xVTVmR6SUFamlx8/Pj4cffpgJEyZcsDfWlUhMTOSFF15g5syZZGdnF/n6sqSWHg+25b/w5UPg5QeProSw+lZXJBXcnsQ0Zq3czfzkJ8A7hZyT10JKH57o24CHutfFx0vLI4hcqaL8/C5S6Dl48CCxsbElrY8DBw5Qp45rLeWv0OOhMpPhvfaQedpccbnHU1ZXJBVYWnYeU3/cy+yVB/Gq+gN+1X+E/FCuC32Lv3ZvRN3qla0uUcTtlFn3VmkEHsDlAo94sEXjzcAT3hS6Pm51NVKBzd9ygue+3c6ptBxs3qlUrr4cJ/BKr3FcX6+d1eWJeATNexTPdWw9bP4MsMHAd8Bb+xRJ2ZixPJ6XFuwEoHa1AGo1XsCm5FxaVW/FdXUHWFydiOdQ6BHP9fMU83urOyCmg7W1SIX13pK9vLF4DwAPXlWHwMjFfLx9Od42b8Z2GKtZWSLlqFRCz7x585g7dy67d+8mNTWViw0Tstls7N+/vzQ+TqTkTu2GXd+Zj7s9YWkpUjEZhsFbP+zhnSX7ABh9TUOiYjbxwqqZAEzqOokW1VtYWaKIxylR6DEMg9tvv50vv/zyokEHzLDzxwULRSz3y9vm98Y3QHhja2uRCic5I5cJX21lwdYEAMYNaEzTescYteTvAAxvNZxB9QdZWKGIZyrRvMgZM2Ywb948WrZsyaJFi7jllluw2Wzs3r2b7777jsGDBwMwYcIE4uPjS6VgkRI7ewS2fG4+vupJa2uRCiduZyL93lrOgq0JeNttPDkglCNes3hsyWM4DSc31buJR1o9YnWZIh6pRC09n3zyCX5+fnz//fdERkby2WefAdCgQQMaNGjAddddR+/evRk+fDi9evWidu3apVK0SImsmgbOfHMhQm0mKqUkPSefv3+3gzm/HsHmlU5MzRM0a3CYWYficBpOAK6NvZZJXSep5VvEIiUKPdu2baNLly5ERkYCFP5F/n131sMPP8zUqVN5/fXX6dOnTwnLFSmhjNOw4Z/m4+6jra1FKoy1B5L529xNnMhbTWCdJXj5J3IW+MXs3aJHzR480vIRjeERsViJQk9WVhZRUVGFz/38/IBzCwUVaN26NT/88ENJPkqkdKx6D/IyIaoV1O1tdTXi5nLyHUxZvId//ByPd5VfCKj5beFrDUIb0CGiAzfWu5FmYc0srFJECpQo9ERERHDq1KnC5+Hh4QDs27ePdu3OLbaVnJzscttOiAdKPQGrPzAf9xij/bWkRHYcT+XJzzexOzEV32pL8QtfDMBdTe7i4ZYPU9W/qsUVisgflWggc/369c8boNyhQwcMw2D69OmFx3bu3MmyZcuoV69eST5KpOR+etXcVLRmR2h8vdXViJvKzXcybek+bpq2jH2pmwmJ+aow8AxvNZyxHcYq8Ii4qBK19PTr149nnnmGHTt20LRpU/r160dMTAyzZs1i06ZN1KpVi7i4OPLy8rj33ntLq2aRokvaBxv+ZT6+5nm18kiR5TuczNtwlCkrviXF5yf86u/HZs/F+dvrT7V/inub6d85EVdWotAzZMgQ8vPzycrKAsDX15fPP/+cQYMGsX79etavXw/ATTfdxOOPa18jsdCSF8BwQIP+ULur1dWIGziTfYaTmScB2HTkDFOX/0SK7494VU3A57dzQv1D6RzVmRvq3kCPmj2sK1ZErkiRdlm/UllZWSxfvpzk5GSaNGlC69atS/sjSp12WXdDTidknIQ/+yN8ahd8MgiwwaO/QIQGlcqF8hx5HEg9wIpjK1h6eCmbT23G4MI/Wz42f25teCu3NLyJhqENsdtKNEpAREqozHZZv1IBAQH079+/LN5aBJLjYdNn5lfqsSu/rtUdCjwCmMtq7Dmzh+VHl7P+5HoOpRzieMbxwvV0ClTyqkJajsN87B3IfS1v466mdxDiF3KxtxURF6cNR8V9pCXA1yNh3++XP7DBlfxPOzgaeo8vs9LEenmOPPan7Cc9N52MvAwy8zMLX3MaTpKzkzmZeZLEjEQ2ntpIQkbCBe8R4B1A2/C29IrpxZFjsUz74TQAd3aqxYs3NcfLrrFgIu6sSKFnw4YNtG3btsQfWlrvIx7k8Br4772QngDYoF4faHM3NLoOfPytrk4sdjb7LPd8fw8HUw9e8TV+Xn50jupMt+huNKjSgFrBtageUJ3MXAfjvtzKt5uPA/BIz3qMvbaRVlEWqQCKFHo6dOjAbbfdxnPPPUfjxkXfpHH79u08//zzzJs3D4fDUeTrxUOt+xgWPAXOPKjeBAZ/AmENrK5KXITTcDJuxTgOph4kwDuAiMAIKvtUxt/b/7zxNqH+oYQHhhMRGEGdkDp0iOxAgHfAee+172Q6j/x7PftOpuNttzH++iYM61anvG9JRMpIkULPY489xrRp05g7dy5dunThvvvu4+qrr6ZOnUv/oxAfH88PP/zA7NmzWbt2LV5eXprJJVfm5E5YNB72x5nPm94EN70PfpWtrUtcyodbPuSXY7/g7+XPJwM+oVHVRkV+j+w8Bx//cpB3l+wlM9dBRLAf0+5sS/tYrbcjUpEUefbWzp07eeqpp/j+++8Lj1WvXp0mTZpQrVo1goODSU1N5fTp0+zcufO8FZuvu+46XnvtNZo0aVJ6d1BKNHvLRRgGnN5nrpy8frY5zdzuA72fMXdEVxeD/M4vx37h0R8fxcDg2U4v0Dm86BMo1h86w+uLdnPsrLn0Rpe61XhnSBuqB/mVdrkiUgaK8vO72FPW9+7dy3vvvcfXX3/N4cOHL3lerVq1GDRoECNGjKBBA9ftklDoKWeGYbbknD0EaSfMLSIStsKRNZCVfO68xjfANS9ANa3oLefbfXov93w/lCxHGv5Z3Ug6NPBPVy+4nKgQf57q34hBraOxa8CyiNsol9Dze/Hx8WzcuJHExERSUlKoUqUK4eHhtG3b9rJdX65EoaccHV4DcS/AoRUXf93LD2I6Qs+xUKd7+dYmbmHdifU8tHg4+WTiyIom89AjYPgQ4ONV5MbASn7e3Nc1lgeuqoO/j1fZFCwiZabc1+mpW7cudevWLY23korszCH4fgzsWWg+9/KF8KYQXAOCIqFafYjpBJEtwdvX2lrFZS0++CNP/TQGJ3k4s2pxV+zzdO1di+Y1QtQlJSKXpXV6pPx8MQyOrQeblzndvOcYCKlpdVXiJvIceczYMpMPtrwPGDjSm/Bmz9cZ0Ly21aWJiJtQ6JHycXyjGXi8fOGRFVC96DNsxHP9mvArE1e8wLGMgwA4Ujrybr+/c3XjKGsLExG3UuLQk5+fz9y5c4mLi+P48eNkZ2df9DybzUZcXFxJP07c1bqPze9Nb1LgkSuSmZfJz8d+Zn78ApYeWQKAM78yXmdu5MMbH6R7w+oWVygi7qZEoefUqVP069ePLVu28GfjobWaqQfLToWtX5iP299vbS3isrLys9hxegdbT21l48mNrDy+kmzHuf9E5Z7pxFVV72HyXzsRHqRVuEWk6EoUesaMGcPmzZupX78+jz76KA0aNCAoKKi0apOKYut/IS8DwhpBrS5WVyMuIt+Zz7akbaw6voqVx1eyNWkrDuMPK7XnVSMnpTmBue158/p+XN8iSv+BEpFiK1Ho+e6774iIiGD16tVUraqVS+UiDONc11b7+7W4oAAQfzaev/30N/ad3Xfe8eoB1WkR1gJHdgyLfg0mL6sGTaNCmPFAe6KrBFzi3URErkyJQk9WVhbXXnutAo9c2tF1kLgNvP2h1WCrqxEXsCB+Ac+teo6s/CyCfILoUqMLXWp0oXNUZxJOB/LK97tYd+gMANe3iOL121oS6Ks5FyJSciX6l6RBgwZkZWWVVi1SEa2bZX5v/hcICL2iS/YkprFwWwLOkq+bKS4kx5nOutRP2ZWxGIAo3+b0rvokAflVOHoYXlyTwKLtiQD4+9h5/OqGPNKzrrqzRKTUlCj0PPDAA4wdO5ajR49Ss6bWW5E/yDoD2780H7cbdkWXpGbnce/MtSSkXnwWoLghWw6+VVfiW+0nbF7m72tOUm/2nLqGPZwCzu3PZ7fBbe1iePKahkSGaLCyiJSuEoWekSNH8vPPP9OnTx/effddrrnmGux2e2nVJu5u8xzIz4aIFlCz/RVd8ur3u0hIzSa6SgC9G2tKsjvLMzI5mBvH/pz55BipAATZY2jufxfh9VrCH7ZT8/f24rb2MTSK1GQIESkbJe4o//DDD+nZsyfXXXcd3t7eREVFXTT42Gw29u/fX9KPE3dx3gDm+65oAPPq+NN8usbcvPbN21vRuW61MixQytLnuz7n7Y1vk5abBkCtoFoMbz2ca2Ovxcuu/a1ExBolCj1Hjhyhe/fuHDlyBMMwyMvLu+SO6+qX9zCHVkLSbvCpBC1u/9PTs/McjJu3BYA7O9VS4HFThmHw9oa3mbltJgB1Q+pyf/P7ua7udfjYfSyuTkQ8XYlCz9ixYzl8+DBXXXUVo0ePpkGDBlSuXLm0ahN3tv63Vp4Wt4L/n+9aP/XHvRw8nUlksD/jBjQu4+KkLOQ783l+1fN8te8rAB5r8xgPtHgAu01d3iLiGkoUen788Udq167NDz/8gJ+fdjeW32Schh1fm4/b//kA5qT0HGb8HA/A3wc1J9hfLQLuJt+Zz+hlo1l6ZCl2m53nujzHzQ1utrosEZHzlOi/YFlZWXTs2FGBR8636VNw5EKNNubXn1gdfxqH06BJVDB9m0aUQ4FS2t7Z+A5LjyzFz8uPqb2mKvCIiEsqUehp2rQpycnJpVWLVARO57murSucpr46/jQAXTSOxy3FHYrj423m7/lLV71E71q9La5IROTiShR6Ro0axU8//cS2bdtKqx5xV04nxC+DL4ZBcjz4BZsLEl6BVfvN0NO5rlb2djcHUw4y/pfxANzb9F76x/a3uCIRkUsr0Zieu+++mx07dtCnTx9efPFFBgwYQK1atUqrNnEXuxfCgqcg5Xcz97o9Bn5/Pqj9ZFo2+09lYLNBxzoKPe4kMy+TJ5c9SUZeBm3D2/JEuyesLklE5LJKFHq8vM6ttzF8+PDLnmuz2cjPzy/Jx4krOrXHbN3JywS/EHO2Vpu7r2gsD8CaeLN7tElkMFUCfcuyUilFhmEw8ZeJ7Du7j7CAMN7o+YampIuIyytR6DGKsDdSUc4VN5GXDV/cbwaeOj3hzs/Bp2g7YReO56mn8Tzu5KOtH7H40GK87d682fNNqgdq9WwRcX0lCj1Op7O06hB39OMkSNwKgWFwyz+KHHjgXOjRYoTu46cjP/HuxncBeKbTM7SNaGtxRSIiV0arhknx7P4e1kw3Hw/6AIIii/wWJ1N/N54nVuN53EH82XjG/jwWA4PBjQZzW8PbrC5JROSKKfRI0eWkw9cjzMedR0DDfsV6m9UHzPE8TaOCCQnUeBBXt/L4SoYuHEpGXgbtItoxtsNYq0sSESmSEm84Kh5o06eQeRqq1oO+k4r9Nuracg9Ow8msbbN4d+O7OA0nTao24c2eb+LjpaAqIu5FoUeKxumA1e+bj7uMAO/ir8at0OP6MvIyeObnZ1hyZAkAN9e/mfGdx+PnpVXYRcT9KPRI0eyaD2cOQkAotBpS7Lc5mZpNvMbzuLTj6ccZtWQUe87swcfuwzOdnuHWhrdaXZaISLEp9EjRrJpmfm//APgGFvtt1mg8j0vbdHITjy99nOTsZKr5V+PtPm/Tqnorq8sSESkRhR65ckfXwZHV4OULHR8q0VsdSMoAoHmNkNKoTErR6hOrGf7jcPKceTQKbcS7fd4lqnKU1WWJiJSYS8/emjZtGrGxsfj7+9OpUyfWrl17RdfNmTMHm83GoEGDyrZAT7PqPfN7i9uKNUX99xJSswGICPEvaVVSig6nHuZvy/5GnjOPXjV78a8B/1LgEZEKw2VDz+eff87o0aOZNGkSGzZsoFWrVvTv35+TJ09e9rqDBw/yf//3f3Tv3r2cKvUQZw/Djq/Nx50vv+XIlUhMMUNPZLBCj6vIyMvgsSWPkZqbSsuwlrzR6w0CfYrfhSki4mpcNvRMmTKFhx56iGHDhtG0aVOmT59OYGAgs2bNuuQ1DoeDu+66i+eff566deuWY7Ue4Jd3wHBC3V4Q2bzEb1fQ0hMZollArsBpOHn656fZn7Kf6gHVeav3W5qhJSIVjkuGntzcXNavX0/fvn0Lj9ntdvr27cuqVasued0LL7xAeHg4DzzwwJ9+Rk5ODqmpqed9ySWkJcCGf5mPu/+tVN4ysaB7Sy09LuH9Te+z9MhSfO2+vN37bcIDw60uSUSk1Llk6ElKSsLhcBAREXHe8YiICBISEi56zYoVK5g5cyYzZsy4os+YPHkyISEhhV8xMTElrrvCWvkuOHIgpjPElrzbMDffSVJ6LqDuLVew6OAiPtzyIQCTuk6iRfUWFlckIlI2XDL0FFVaWhr33HMPM2bMICws7Iquefrpp0lJSSn8OnLkSBlX6aYykmDdb12KPZ4Cm63Eb3kqPQcAHy8bVSv5lvj9pPh2J+9m4i8TAbi36b3cWO9GiysSESk7LjllPSwsDC8vLxITE887npiYSGTkhbOG9u/fz8GDBxk4cGDhsYId4L29vdm9ezf16tU77xo/Pz/8/DRm4U+teg/yMqFGG6h/dam8ZcJvg5jDg/yxlUKIkuJJzk7msSWPkZWfRZeoLjzZ7kmrSxIRKVMu2dLj6+tLu3btiIuLKzzmdDqJi4ujS5cuF5zfuHFjtm7dyqZNmwq/brzxRnr37s2mTZvUdVVcmcmw9rfuwh5jSqWVB86N54nUdHXL5Dny+Nuyv3E84zi1gmrxes/X8ba75P+BRERKjcv+Kzd69GiGDh1K+/bt6dixI1OnTiUjI4Nhw4YBcO+99xIdHc3kyZPx9/enefPzZxRVqVIF4ILjUgSrP4DcdIhoDo0GlNrbJmi6uqWchpOJKyeyLnEdlXwq8U6fdwjx0yKRIlLxuWzoGTx4MKdOneLZZ58lISGB1q1bs3DhwsLBzYcPH8Zud8mGqorhzCFY+Y75uJTG8hTQzC1rvbPhHebHz8fb5s2UnlOoV6Xen18kIlIBuGzoARg5ciQjR4686GvLli277LWzZ88u/YI8ycKnIT/bnK3V9KZSfWut0WOdObvmMHPbTACe6/ocXaO7WlyRiEj5UVOJXGjPItg9H+zecN0bpdrKA+e6t9TSU76WH13O5LWTARjZeiQ31S/dMCsi4uoUeuR8eVmw4CnzcefhEN641D+icCCzQk+5OZx6mHHLx+E0nNzS4BYebvmw1SWJiJQ7hR453y9vw9lDEFQDeo4p9bc3DON33VsKPeUhMy+Tx5c+TlpeGq2qt2JCpwlaKkBEPJJCj5yTehxWvGU+7v8S+AWV/kdk5ZOdZ66hpO6tsmcYBs+ufJZ9Z/cRFhDGlF5T8PHysbosERFLuPRAZilnP71mDl6u1QWa3VwmH1HQylMl0Ad/H68y+QwxOQ0n7258l0UHF5kztXpN0Z5aIuLRFHrElBwPGz8xH1/9bKkPXi5QEHoigtTKU5Yy8zKZ8MsEfjj0AwDjOo6jTXgbi6sSEbGWQo+Ylr0CznyodzXULrtpzIVr9Gg8T5k5ln6Mx5Y8xp4ze/C2e/Ns52e5uUHZtNyJiLgThR6Bkzthy3/Nx1dPLNOPSixcjVlr9JSFhIwE7llwD6eyTlHVvypTe09VC4+IyG8UegSW/B0woMmN5saiZShB09XLTGZeJqOWjOJU1inqhdRj+jXTiax04Qa9IiKeSqHHk53eD6vfh13fATboPb7MP1LdW2XDaTh5ZsUz7EreRVX/qrzf930FHhGRP1Do8UQpR2HhONj5HWCYxzr9tUwWIvwjtfSUjfc2vkfc4Th87D683fttalSuYXVJIiIuR6HHE817CA6vNB836AddRkKdHuXy0QkpOYDW6ClNy48uZ8bWGQA83/V5Woe3trYgEREXpdDjaQ6vNgOPly88GAdRLcvto/McTk5nmKFHqzGXjqz8LF5e8zIAdze5m4H1BlpckYiI69KKzJ5mxVTze6s7yjXwAJxMy8EwwMfLRtVA33L97IpqxpYZHEs/RlSlKEa1GWV1OSIiLk2hx5Mk7oA93wM26Pp4uX98we7q4UH+2O3a+6mk4lPi+Xj7x4C5+GCgT6DFFYmIuDaFHk+y8h3ze9MbIax+uX984cwtrdFTYoZh8NLql8h35tOzZk96x/S2uiQREZen0OMpzh6GrXPNx92esKSEgpYejecpue8PfM/ahLX4efkxruM47ZouInIFFHo8xapp5jYTdXtBdFtLSkhMK2jpUegpqX/v/DcAD7R4gJpBNS2uRkTEPSj0eAKnEzaaPyStauWB329BodBTEsfSj7E1aSs2bNzW8DaryxERcRsKPZ4g7QTkpoPdG2K7W1bGCXVvlYrFBxcD0D6yPWEBYRZXIyLiPhR6PMGZg+b3kBjwsm5ppqNnsgCoGRpgWQ0VQUHo6V+7v8WViIi4F4UeT1AQekJjLSshz+HkREpB6NHU6uI6mnaUbae3YbfZubr21VaXIyLiVhR6PIELhJ4TZ7NxGuDrbad6ZU1ZL67Fh8xWng4RHdS1JSJSRAo9nsAFQs/RM5mA2bWlhQmLr6Brq19sP4srERFxPwo9nsAFQs+R30JPjLq2iu1I2hG2n95udm3VUteWiEhRKfR4AlcIPckaxFxSBa08HSI7UC2gmsXViIi4H4Weii43AzJOmo9doaWnqlp6isMwDBYeXAhA/1jN2hIRKQ6Fnoru7GHzu38VCKhiWRkF09XVvVU8iw4tYlfyLnztvuraEhEpJoWeis4FurYAjiQXtPSoe6uo0nPTeW3ta4C57URV/6oWVyQi4p4Ueio6Fwg92XkOTqblAFqjpzje3/w+p7JOERMUwwMtHrC6HBERt6XQU9G5QOgp6Nqq5OtFaKCPZXW4o93Ju/ls52cAjO80Hj8vrXEkIlJcCj0VnUuEnnODmG02rdFzpZyGkxdXv4jDcNCvdj+6RXezuiQREbem0FPRFYae2paVcER7bhXLR1s/YvOpzQR6BzKmwxiryxERcXsKPRWZYbhGS09ywWrMGs9zpZYcXsK7G98FYGzHsURUirC4IhER96fQU5GlJ0J+Ntjs5g7rFimcrq41eq7I3jN7efrnpwEY0ngItzS4xeKKREQqBoWeiqyglSekJnhZN4D43BYU6t76M2ezz/LYksfIzM+kY2RHnurwlNUliYhUGAo9FZkLdG3BuTV61L11ebmOXJ5Y9gRH048SXTmaN3u+iY9ds91EREqLQk9F5gKhJz0nnzOZeYAWJrwcp+FkwooJrE9cT2Wfyrzb512q+FexuiwRkQpFoacic4HQUzBdvUqgD0H+arW4lHc2vMP3B7/H2+bNW73fokFoA6tLEhGpcBR6KjIXCD0Fu6trz61L++/u/zJz20wAnu/2PJ2jOltckYhIxaTQU5G5ROgpGM+jrq2LWXNiDS+veRmAEa1HcGO9Gy2uSESk4lLoqajysiDthPk4tI5lZRz53WrMcr6jaUf5209/w2E4uKHuDfy15V+tLklEpEJT6Kmozh4xv/sFQ0CoZWUUrtGjlp7zZOZl8tjSx0jJSaFZtWZM6jJJW3SIiJQxhZ6K6vfbT1j4w1TT1S+U78xnwi8T2HtmL9X8qzG191T8vf2tLktEpMLztroAKSMuMJ7HMIzfrcaslh6AHad38NzK59iZvBNvuzdTe08lslKk1WWJiHgEhZ6KygVCT0pWHuk5+YBaerLys3h/0/v8a8e/cBpOgn2Deb7r87QOb211aSIiHkOhp6JygdBTMF29epAf/j5eltVhtZXHV/LCqhc4ln4MgAGxAxjTcQxhAWEWVyYi4lkUeiqqgtBTJdayEgpmbnnqdPWz2Wd5fd3rfLP/GwAiK0UysfNEetTsYXFlIiKeSaGnIjIMF2npKdho1PO6trae2soTS5/gZNZJbNi4s8mdjGoziko+lawuTUTEYyn0VEQZSZCXAdigSoxlZXjqIOav933NC6teINeZS52QOrzY7UVaVW9ldVkiIh5PoaciKmjlCY4Gbz/LyihcmNBDWnrynfm8ue5N/r3z3wD0junN5O6T1bojIuIiFHoqIhfo2gLPWqMnMy+TMcvH8NPRnwB4pNUjPNrqUew2LYUlIuIqFHoqIhcIPZ60Rk9SVhIj4kaw4/QO/Lz8ePmql+kX28/qskRE5A8UeioiFwg9p9JzyMl3YrdBVEjFDT37z+5n+I/DOZ5xnFC/UN7p847W3hERcVEKPRWRC4SegjV6IoP98fWumF08y44sY9zP48jIy6BWUC0+6PsBtYJrWV2WiIhcgkJPReQCoedowRo9FXB3dcMwmLltJu9seAcDg/YR7ZnSawqh/tZt7CoiIn9Ooaeiyc+FVHPlX2tDT8Hu6hUn9BiGwZqENczeNptfjv8CwOBGgxnbcSw+dh+LqxMRkT+j0FPRpBwBDPCpBJWs2+agcGHCCjCI2TAMvt7/Nf/c/k/2nd0HgLfNm6c7Pc3tjW63uDoREblSLj3YYtq0acTGxuLv70+nTp1Yu3btJc+dMWMG3bt3JzQ0lNDQUPr27XvZ8yusMwfM76GxYLNZVsa5LSjcu6XHMAymrJ/CxF8msu/sPgK8A7ij0R3876b/KfCIiLgZlw09n3/+OaNHj2bSpEls2LCBVq1a0b9/f06ePHnR85ctW8aQIUNYunQpq1atIiYmhn79+nHs2LFyrtxiLjCeB84NZI5x4323DMPg7Q1vM3v7bACGtxpO3G1xjO88ntiQWEtrExGRonPZ0DNlyhQeeughhg0bRtOmTZk+fTqBgYHMmjXroud/+umnDB8+nNatW9O4cWM++ugjnE4ncXFx5Vy5xVwg9DicBsfPFqzR474tPdM2TWPmtpkAPNPpGR5t/ShBvkEWVyUiIsXlkqEnNzeX9evX07dv38Jjdrudvn37smrVqit6j8zMTPLy8qhatepFX8/JySE1NfW8rwrBBUJPQmo2+U4DHy8bEcH+ltVRXHmOPCavmcyHWz4EYGyHsQxpPMTiqkREpKRcMvQkJSXhcDiIiIg473hERAQJCQlX9B5jx46lRo0a5wWn35s8eTIhISGFXzEx1m3MWaoKQ09ty0ooGMRco0oAXnbrxhUVR0JGAsMWDeOzXZ8B8H/t/4+7m95tcVUiIlIaKuTsrVdeeYU5c+awbNky/P0v3tLw9NNPM3r06MLnqamp7h98DAPOHDIfW7owoWttNJqWm8bZ7LMXfc3A4GzOWRIzEzmefpxZ22aRnJ1MkE8QL3d/mV4xvcq1VhERKTsuGXrCwsLw8vIiMTHxvOOJiYlERkZe9to33niDV155hR9//JGWLVte8jw/Pz/8/KzbgbxMZJ2BnN+66apYtzKwlXtu5TnyWHl8JYsOLmLv2b0cSz9GWm5akd6jcdXGTOk1hZggNw/BIiJyHpcMPb6+vrRr1464uDgGDRoEUDgoeeTIkZe87rXXXuOll15i0aJFtG/fvpyqdSEF09WDosDHullTVkxXP5l5kg83f8iiQ4tIyUm54PUA74BL7nge5BtERGAE4YHhNK7amHub3ou/t/uNRRIRkctzydADMHr0aIYOHUr79u3p2LEjU6dOJSMjg2HDhgFw7733Eh0dzeTJkwF49dVXefbZZ/nss8+IjY0tHPtTuXJlKleubNl9lCsXGMQMcPS36eo1y2m6+toTa3lq+VMkZycDUM2/GgPqDKBzVGeiK0dTo3INAn1co6tNRESs47KhZ/DgwZw6dYpnn32WhIQEWrduzcKFCwsHNx8+fBi7/dz/3D/44ANyc3O59dZbz3ufSZMm8dxzz5Vn6dZxldBzpmA15rINGoZh8PH2j3l7w9s4DScNQhvwVPun6BjZES+7V5l+toiIuB+XDT0AI0eOvGR31rJly857fvDgwbIvyNWdjje/V7Fu5lZuvpMTqdlA2Q5k3pC4gfc2vcevCb8CcGO9G5nQeQIB3u67GKKIiJQtlw49UkSHVpjfo9taVsLxs1kYBvj72Amr7Fuq753nzGND4gY+2voRq0+sBsDX7su4TuO4tcGt2CzcdkNERFyfQk9FkXzA7N6y+0DtbpaV8ftBzCUNIam5qcSfjWfPmT2sPL6SNSfWkJ6XDpgbfg5qMIiHWjxEjco1Sly3iIhUfAo9FUX8MvN7TEfws27g9qLt5gDy2sUczxN/Np5/7fgXy48u51TWqQter+JXhb61+/JgiweJrhxdolpFRMSzKPRUFPFLze91e1lWws97T/Hv1YcBuK9b7BVfl+vIZW3CWv6z6z8sP7r8vNfCA8OpF1KPthFt6VajG02rNdUgZRERKRaFnorA6YADv4UFi0JPSmYeT83dAsC9XWrTvUH1y56fnZ/Nd/Hf8dPRn1hzYg1Z+eY0dxs2esf0ZkiTITSr1kwbfIqISKlR6KkITmw2V2P2C4Ya1gxifu7b7SSkZlMnrBLjBjS+5HmGYfDj4R9549c3OJ5xvPB4WEAYfWv15e6md1M72LrZZyIiUnEp9FQEBeN5YruDV+n9luY6ctlxegf5zvzLnrcq/jTf7NqDdyA81K8525M3kufMI9eRS44jp/B7jiOHpUeWsubEGgAiAiO4vdHtdI/uTqOqjS65YrKIiEhpUOipCArG89TrXWpvuSFxA8+ufJZDqYeu6PzAWPP75I1/fq6v3ZdhzYdxf/P7tVKyiIiUG4Ued5ebCYfNNWtKYzxPZl4m72x8h892foaBQbBvMFX9q1703PScfE6m5mAAlf28iQg+t4Grr5cvfl5++Hr54u/lX/i8WkA17m5yNzWDapa4VhERkaJQ6HF3R1aDIxeCa0K1+iV6q18TfuXZX57laPpRAG6ufzP/1+H/CPYNPu+8PIeTf648yEsLdmIYcGu7mrz6l5Z42bU4oIiIuC6FHne3/3dT1Yu5GGBGXgZvrX+Lz3d/DkBkpUie6/Ic3aLPX+QwIyefOb8eYdaKAxw7a862GtqlNpMGNsOuwCMiIi5OocedOZ2w70fzcTHH8+w9s5dRS0ZxLP0YALc1vI3R7UZT2ffcAofbjqUwd90R/rfxGKnZ5qDmsMq+PNqrPvd3i9X2DyIi4hYUetzZ8tfg5A7w9i/WeJ7dybt5aPFDnMk5Q3TlaJ7r+hydozoXvr5ibxKvLNzJtmOphcfqhFXioe51uaVtNP4+WiRQRETch0KPu9q1AJZNNh9fPwUqhRXp8h2nd/DwDw+TkpNC02pN+cc1/yDEL6Tw9S/WH2XsvC04nAa+XnauaRbB7e1juKp+mMbuiIiIW1LocUen9sCXD5uPOz4Mbe664ksNw+Cnoz/xzIpnSMtNo2VYSz645oPCwcqGYfDBT/t5beFuAG5qXYNJA5tRtVLp7pguIiJS3hR63E36Sfj8LshNM3dT7//yFV2W78xn0cFFzNo2iz1n9gDQunprPuj7QeH4nf2n0nl/6X7mbTBnb/21R13GXttYg5RFRKRCUOhxJ8c3wpy7IPUYBEfDbf8EL58/vSwzL5P7F93P9tPbAQj0DuS2hrcxvPVwArwDmL/lBP9efYhV8acLr5l4Q1MeuKpOmd2KiIhIeVPocRdb5sI3IyE/G6o1gCFzoPLlN/UEs7vq2ZXPsv30doJ8gxjadCh3NL6DEL8QkjNyGfnpOpbsOgmA3QZ9Godz/1V16FqvaGOEREREXJ1Cjzv4eQrEPW8+btAf/jID/EMuf81v/rn9nyw6uAhvmzfTrp5Gm/A2AKyOP83jczaSmJqDr7edh7vX5c5OtahRJaCs7kJERMRSCj2ubt2sc4Gn2xNw9bNgv7Kp4qtPrOatDW8BMLbj2MLA89HP8by8YCdOA+pWr8R7Q9rStEbw5d5KRETE7Sn0uLLt/4PvRpuPezwFfSZc8aXH04/z1E9P4TSc3FTvJgY3GoxhGLwdt5epP+4F4C9ta/LCTc2o5Kc/BiIiUvHpp52r2vcjzHsIMKD9/dB7/BVfmp2fzRNLn+BszlmaVG3ChM5mWHp5wU5m/HwAgKf6N2JE75Lt1SUiIuJOFHpcTW4mLH0JVr8PhhOa3QzXvXHF+2oZhsHfV/+dnck7qeJXham9p+Lv7c9z32xn9sqDgGZmiYiIZ1LocSUHlsM3o+DMQfN5yzvgxneueAwPwH93/5ev93+N3WbntR6vUaNyDTYePsPslQex2eDlm1swpGOtsqlfRETEhSn0uIpt8+CLBwDDXIPnhregYf8ivcXmU5t55ddXAHii7RN0qdEFgCk/mIsR/qVtTQUeERHxWAo9Vju5C36dARv+BRjQcrDZneVf9NlUU9dPJd+ZzzW1r+G+ZvcBsCb+ND/vTcLbbuPxqxuUbu0iIiJuRKHHKo58WDwe1kw/d6zZLTDogyJ1ZxWIT4lnXeI67DY7YzqMwWazYRgGby42W3lu7xBDTNXA0qpeRETE7Sj0WCHrLHwxDPYvMZ83vsHcOLROjysesPxHX+z5AoAe0T2IrBQJwIp9Saw9mIyvt51RfTRTS0REPJtCT3k7vR8+Gwyn94JPINzyD2gysERvmePI4Zv93wBwW6PbAHMW1xu/tfLc1akWUSFaaVlERDybQk95OnsYPh4A6YnmYOUh/4GoViV+28UHF5OSk0JUpSi61egGQNzOk2w+cpYAHy8e7VWvxJ8hIiLi7hR6ykt2Cnx6uxl4wpvCPV9BUESpvHVB19YtDW7By+6F02nw5m8ztoZ2jSU8yL9UPkdERMSd2a0uwCPk58Ln98CpnVA5Eu6aW2qBZ9+ZfWw4uQEvmxe3NLgFgIXbE9h5IpXKft78tUfdUvkcERERd6fQU9YMA757Eg78BD6V4K7/QkjNUnv7L/aarTw9a/YkPDAch9MoXJfngavqEFrJt9Q+S0RExJ0p9JS1jZ/Apn+DzQ63zS6VMTwFVh1fVdi1VTCA+ZvNx9h3Mp2QAB8e6K6tJkRERApoTE9Za3Eb7F0MdXtBw36l9rYrjq3g8SWPk+vMpVfNXnSt0ZU8h7NwB/W/9qxLsL9PqX2eiIiIu1PoKWs+AXD7J8Vef+dilh1Zxuhlo8lz5tE7pjdv9nwTu83O7FUHOHQ6k7DKvtzXNbbUPk9ERKQiUOgpD6UUeBxOBx9v/5hpG6eRb5jbTbza41V87D58/uthXpy/A4ARvesT6KvfWhERkd/TT0Y3cTj1MONXjGfTqU0AXF/3el7s9iI+dh/+teogz369HYB7OtdmaJdY6woVERFxUQo9bmB+/HyeX/U8WflZVPapzNiOY7mp3k2k5eTzzvLdvLtkHwAPXlWH8dc3wVaKXWkiIiIVhUKPCzIMg+SMXFKz8vhi/yf8a/f7ALSs1p7HW02gklcYz3+7g7nrjpCR6wBgZO/6/K1fQwUeERGRS1DosdDp9ByW7j7F4eRMjp7J5NiZLE6kZJOQkk2uIx+/yG/wDV0NQO7pHvyy81p+WbHrvPdoEF6ZR3vV45a2pbf2j4iISEWk0GOBbcdS+OfKg3y9+Ti5+c6LnGHgHz0Hn+AtYNjg9I34pnbH1++3l23QrnYo93erQ/cGYWrdERERuQIKPWXMMAyOnsli/aEz/HowmV8PJrMnMb3w9WY1gmkVU4XoKgHUDA2gRpUAVid9xT+2b8HH7sOrPV7lmtrXWHgHIiIiFYNCTxlbtD2RR/69/rxj3nYb17WI4r5usbSJqXJeS83O0zv5+Of3APi/9v+nwCMiIlJKFHrKWOuYKnjbbTSPDqF97VDax1alY52qVL3InliZeZmMWT6GPGcevWJ6MaTxEAsqFhERqZgUespYZIg/257vj7+P12XPy3Pk8dKalziYepDwwHBe7PqixuqIiIiUIoWecnCpwGMYBptPbea7+O9YdHARZ3POYrfZeaX7K1Txr1K+RYqIiFRwCj0WWZ+4nrfWv8XmU5sLj1Xzr8aoNqPoENnBwspEREQqJoWechZ/Np631r/FsqPLAPD38qdfbD+ur3M9HaM64m3Xb4mIiEhZ0E/YcpLvzGf29tm8v+l98px5eNm8uKXBLTza6lGqB1a3ujwREZEKT6GnHMSnxDNxxUS2JG0BoHt0d57q8BR1QupYXJmIiIjnUOgpY3GH4xi7fCw5jhyCfIIY23EsN9a7UTOzREREyplCTxlrVq0ZvnZf2ke057muzxFZKdLqkkRERDySQk8Zi6wUyX9u+A+1gmqpdUdERMRCCj3loHZwbatLEBER8Xh2qwsQERERKQ8KPSIiIuIRFHpERETEIyj0iIiIiEdw6dAzbdo0YmNj8ff3p1OnTqxdu/ay58+dO5fGjRvj7+9PixYtWLBgQTlVKiIiIq7OZUPP559/zujRo5k0aRIbNmygVatW9O/fn5MnT170/JUrVzJkyBAeeOABNm7cyKBBgxg0aBDbtm0r58pFRETEFdkMwzCsLuJiOnXqRIcOHXjvvfcAcDqdxMTEMGrUKMaNG3fB+YMHDyYjI4Pvvvuu8Fjnzp1p3bo106dPv+D8nJwccnJyCp+npqYSExNDSkoKwcHBZXBHIiIiUtpSU1MJCQm5op/fLtnSk5uby/r16+nbt2/hMbvdTt++fVm1atVFr1m1atV55wP079//kudPnjyZkJCQwq+YmJjSuwERERFxOS4ZepKSknA4HERERJx3PCIigoSEhItek5CQUKTzn376aVJSUgq/jhw5UjrFi4iIiEvy2BWZ/fz88PPzs7oMERERKScu2dITFhaGl5cXiYmJ5x1PTEwkMvLiG3ZGRkYW6XwRERHxLC4Zenx9fWnXrh1xcXGFx5xOJ3FxcXTp0uWi13Tp0uW88wF++OGHS54vIiIinsVlu7dGjx7N0KFDad++PR07dmTq1KlkZGQwbNgwAO69916io6OZPHkyAI8//jg9e/bkzTff5Prrr2fOnDmsW7eOf/zjH1behoiIiLgIlw09gwcP5tSpUzz77LMkJCTQunVrFi5cWDhY+fDhw9jt5xqqunbtymeffcaECRN45plnaNCgAV999RXNmze/os8rmLmfmppa+jcjIiIiZaLg5/aVrMDjsuv0lLejR49q2rqIiIibOnLkCDVr1rzsOQo9v3E6nRw/fpygoCBsNlu5fW7BoohHjhypkIsiVuT7q8j3BhX7/nRv7qsi319Fvjcou/szDIO0tDRq1KhxXg/Qxbhs91Z5s9vtf5oQy1JwcHCF/ENeoCLfX0W+N6jY96d7c18V+f4q8r1B2dxfSEjIFZ3nkrO3REREREqbQo+IiIh4BIUei/n5+TFp0qQKuzp0Rb6/inxvULHvT/fmviry/VXkewPXuD8NZBYRERGPoJYeERER8QgKPSIiIuIRFHpERETEIyj0iIiIiEdQ6CkH06ZNIzY2Fn9/fzp16sTatWsvee6MGTPo3r07oaGhhIaG0rdv38ue7wqKcn9ffvkl7du3p0qVKlSqVInWrVvzySeflGO1RVOUe/u9OXPmYLPZGDRoUNkWWEJFub/Zs2djs9nO+/L39y/HaoumqL93Z8+eZcSIEURFReHn50fDhg1ZsGBBOVVbNEW5t169el3w+2az2bj++uvLseKiKerv3dSpU2nUqBEBAQHExMTw5JNPkp2dXU7VFk1R7i0vL48XXniBevXq4e/vT6tWrVi4cGE5Vnvlli9fzsCBA6lRowY2m42vvvrqT69ZtmwZbdu2xc/Pj/r16zN79uwyrxNDytScOXMMX19fY9asWcb27duNhx56yKhSpYqRmJh40fPvvPNOY9q0acbGjRuNnTt3Gvfdd58REhJiHD16tJwrvzJFvb+lS5caX375pbFjxw5j3759xtSpUw0vLy9j4cKF5Vz5nyvqvRU4cOCAER0dbXTv3t246aabyqfYYijq/X388cdGcHCwceLEicKvhISEcq76yhT13nJycoz27dsb1113nbFixQrjwIEDxrJly4xNmzaVc+V/rqj3dvr06fN+z7Zt22Z4eXkZH3/8cfkWfoWKen+ffvqp4efnZ3z66afGgQMHjEWLFhlRUVHGk08+Wc6V/7mi3tuYMWOMGjVqGPPnzzf2799vvP/++4a/v7+xYcOGcq78zy1YsMAYP3688eWXXxqA8b///e+y58fHxxuBgYHG6NGjjR07dhjvvvtuufwsUOgpYx07djRGjBhR+NzhcBg1atQwJk+efEXX5+fnG0FBQcY///nPsiqxREp6f4ZhGG3atDEmTJhQFuWVSHHuLT8/3+jatavx0UcfGUOHDnXp0FPU+/v444+NkJCQcqquZIp6bx988IFRt25dIzc3t7xKLLaS/p176623jKCgICM9Pb2sSiyRot7fiBEjjD59+px3bPTo0Ua3bt3KtM7iKOq9RUVFGe+99955x2655RbjrrvuKtM6S+pKQs+YMWOMZs2anXds8ODBRv/+/cuwMsNQ91YZys3NZf369fTt27fwmN1up2/fvqxateqK3iMzM5O8vDyqVq1aVmUWW0nvzzAM4uLi2L17Nz169CjLUousuPf2wgsvEB4ezgMPPFAeZRZbce8vPT2d2rVrExMTw0033cT27dvLo9wiKc69ffPNN3Tp0oURI0YQERFB8+bNefnll3E4HOVV9hUpjX9TZs6cyR133EGlSpXKqsxiK879de3alfXr1xd2E8XHx7NgwQKuu+66cqn5ShXn3nJyci7oQg4ICGDFihVlWmt5WLVq1Xm/FgD9+/e/4j/HxaUNR8tQUlISDoeDiIiI845HRESwa9euK3qPsWPHUqNGjQv+cLiC4t5fSkoK0dHR5OTk4OXlxfvvv88111xT1uUWSXHubcWKFcycOZNNmzaVQ4UlU5z7a9SoEbNmzaJly5akpKTwxhtv0LVrV7Zv327pZr1/VJx7i4+PZ8mSJdx1110sWLCAffv2MXz4cPLy8pg0aVJ5lH1FSvpvytq1a9m2bRszZ84sqxJLpDj3d+edd5KUlMRVV12FYRjk5+fzyCOP8Mwzz5RHyVesOPfWv39/pkyZQo8ePahXrx5xcXF8+eWXLhfGiyMhIeGivxapqalkZWUREBBQJp+rlh4X9sorrzBnzhz+97//ufSA0aIKCgpi06ZN/Prrr7z00kuMHj2aZcuWWV1WiaSlpXHPPfcwY8YMwsLCrC6nTHTp0oV7772X1q1b07NnT7788kuqV6/Ohx9+aHVpJeZ0OgkPD+cf//gH7dq1Y/DgwYwfP57p06dbXVqpmjlzJi1atKBjx45Wl1Jqli1bxssvv8z777/Phg0b+PLLL5k/fz4vvvii1aWV2Ntvv02DBg1o3Lgxvr6+jBw5kmHDhmG360d3camlpwyFhYXh5eVFYmLieccTExOJjIy87LVvvPEGr7zyCj/++CMtW7YsyzKLrbj3Z7fbqV+/PgCtW7dm586dTJ48mV69epVluUVS1Hvbv38/Bw8eZODAgYXHnE4nAN7e3uzevZt69eqVbdFFUJI/mwV8fHxo06YN+/btK4sSi6049xYVFYWPjw9eXl6Fx5o0aUJCQgK5ubn4+vqWac1XqiS/bxkZGcyZM4cXXnihLEsskeLc38SJE7nnnnt48MEHAWjRogUZGRk8/PDDjB8/3mUCQnHurXr16nz11VdkZ2dz+vRpatSowbhx46hbt255lFymIiMjL/prERwcXGatPKCWnjLl6+tLu3btiIuLKzzmdDqJi4ujS5cul7zutdde48UXX2ThwoW0b9++PEotluLe3x85nU5ycnLKosRiK+q9NW7cmK1bt7Jp06bCrxtvvJHevXuzadMmYmJiyrP8P1Uav3cOh4OtW7cSFRVVVmUWS3HurVu3buzbt68wqALs2bOHqKgolwk8ULLft7lz55KTk8Pdd99d1mUWW3HuLzMz84JgUxBeDRfaWrIkv3f+/v5ER0eTn5/PvHnzuOmmm8q63DLXpUuX834tAH744Yci/ewoljIdJi3GnDlzDD8/P2P27NnGjh07jIcfftioUqVK4VTfe+65xxg3blzh+a+88orh6+trfPHFF+dNM01LS7PqFi6rqPf38ssvG4sXLzb2799v7Nixw3jjjTcMb29vY8aMGVbdwiUV9d7+yNVnbxX1/p5//nlj0aJFxv79+43169cbd9xxh+Hv729s377dqlu4pKLe2+HDh42goCBj5MiRxu7du43vvvvOCA8PN/7+979bdQuXVNw/l1dddZUxePDg8i63yIp6f5MmTTKCgoKM//znP0Z8fLyxePFio169esbtt99u1S1cUlHvbfXq1ca8efOM/fv3G8uXLzf69Olj1KlTxzhz5oxFd3BpaWlpxsaNG42NGzcagDFlyhRj48aNxqFDhwzDMIxx48YZ99xzT+H5BVPWn3rqKWPnzp3GtGnTNGW9onj33XeNWrVqGb6+vkbHjh2N1atXF77Ws2dPY+jQoYXPa9eubQAXfE2aNKn8C79CRbm/8ePHG/Xr1zf8/f2N0NBQo0uXLsacOXMsqPrKFOXe/sjVQ49hFO3+nnjiicJzIyIijOuuu84l1wspUNTfu5UrVxqdOnUy/Pz8jLp16xovvfSSkZ+fX85VX5mi3tuuXbsMwFi8eHE5V1o8Rbm/vLw847nnnjPq1atn+Pv7GzExMcbw4cNdMhgYRtHubdmyZUaTJk0MPz8/o1q1asY999xjHDt2zIKq/9zSpUsv+rOr4H6GDh1q9OzZ84JrWrdubfj6+hp169Ytl7WjbIbhQu1/IiIiImVEY3pERETEIyj0iIiIiEdQ6BERERGPoNAjIiIiHkGhR0RERDyCQo+IiIh4BIUeERER8QgKPSIiIuIRFHpERETEIyj0iIiIiEdQ6BERERGPoNAjIiIiHkGhR0QqpLVr1zJmzBg6duxIZGQkvr6+REREMHDgQH788UeryxMRC2iXdRGpkPr27cvSpUtp1qwZMTExVKpUif3797NhwwYApk6dyuOPP25xlSJSnhR6RKRC+v7772ndujVRUVHnHV+1ahXXXnstWVlZHDhwgOjoaIsqFJHypu4tEamQBgwYcEHgAejSpQsjRowgLy+Pr7/+2oLKRMQq3lYXICJSVk6fPs38+fPZtm0bZ86cIS8vD4C9e/cCsHv3bivLE5FyptAjIhXSjBkzePLJJ8nIyLjkOampqeVYkYhYTd1bIlLhrF+/nr/+9a/k5OTw6quvsmPHDtLT03E6nRiGwYcffgiAhjSKeBa19IhIhTN37lwMw2DUqFGMGTPmgtcLurdExLOopUdEKpzk5GQAateufcFr2dnZzJs3r7xLEhEXoNAjIhVOkyZNAPjnP/9JWlpa4fHs7GyGDx/OgQMHrCpNRCyk0CMiFc6wYcOoXbs2GzdupE6dOtx8883ceuut1K5dmy+++EKLEop4KIUeEalwqlSpwrp16xg+fDhVqlTh+++/Z9WqVfTr148NGzbQunVrq0sUEQtoRWYRERHxCGrpEREREY+g0CMiIiIeQaFHREREPIJCj4iIiHgEhR4RERHxCAo9IiIi4hEUekRERMQjKPSIiIiIR1DoEREREY+g0CMiIiIeQaFHREREPIJCj4iIiHiE/wfrb8wRdnV65gAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot the mass assembly history of a random halo\n",
    "plt.plot(scales, ma[0])\n",
    "plt.plot(scales, ma[10])\n",
    "plt.plot(scales, ma[20])\n",
    "plt.xlabel('a', fontsize=16)\n",
    "plt.ylabel('m(a)', fontsize=16)\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see all our monotonically increasing and by the definition end up at $m(a) = 1$ when $a = 1$ or $z = 0$. \n",
    "\n",
    "The dataset also contains ROCKSTAR present-day properties of the haloes. This includes concentration, centre-of-mass offset, virial ratio, spin, and triaxiality. Let's access these properties:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((10000,), (10000,))"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cat = mah_data['cat'] # astropy table with present-day halo properties\n",
    "cvir = cat['cvir']\n",
    "xoff = cat['x0']\n",
    "t_over_u = cat['t/|u|']\n",
    "spin = cat['spin_bullock']\n",
    "c_to_a = cat['c_to_a']\n",
    "\n",
    "cvir.shape, spin.shape # 1 dimensional quantities, one per halo"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup the model and train\n",
    "\n",
    "For the purposes of this tutorial we will use the MAH information and halo properties as features and targets for prediction respectively. In this section, we will walk you through how to setup the MultiCAM model and train it.\n",
    "\n",
    "First we import relevant classes from `mah.models` module:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "from multicam.models import LinearRegression"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `LinearRegression` model implements the option of using the MultiCAM framework with a linear regression algorithm for connecting features and targets within the gaussianized space (see Fig 1 of the MultiCAM paper).\n",
    "\n",
    "First we need to create an instance of the model. For now, let's use the model to predict concentration from MAH only."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_features = ma.shape[1] # number of features is the number of scales\n",
    "n_targets = 1 # number of targets is 1, the present-day halo property\n",
    "model = LinearRegression(n_features, n_targets, use_multicam=True)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lastly before training, let's split the data into training and testing: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "from multicam.models import get_tt_indices\n",
    "\n",
    "n_haloes = len(cat)\n",
    "train_idx, test_idx = get_tt_indices(n_haloes, test_ratio=0.25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "ma_train = ma[train_idx]\n",
    "ma_test = ma[test_idx]\n",
    "cvir_train = cvir[train_idx]\n",
    "cvir_test = cvir[test_idx]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last flag means that we are using the MultiCAM framework which includes a gaussianization step and matching quantiles. If this flag is turned off, then the model will simply be a linear regression model.\n",
    "\n",
    "Now we can proceed with training, the only thing to be careful with is that all features and targets need to be reshaped to be 2D arrays if they are not already. To train we use the `.fit()` method of the model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.fit(ma_train, cvir_train.reshape(-1, 1)) # reshape to 2D array"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and that's it! We are ready to predict with this model."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prediction and results"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To predict, simply use the `.predict()` method of the model on the test features: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "cvir_pred = model.predict(ma_test)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compare them to the test concentration via a histogram: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkUAAAG2CAYAAABxkPLrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABg2ElEQVR4nO3deVwVZf//8ddhFURAXECUrTTFcsVS0tJbDTQzS/LObvdKuw0zM5f85Z5KWbnepq1qmbelpablvtCt4oZaLuQWhneyZAa4gnDm94e35xshioeDh+X9fDzOQ89c11zzOeMJ3s1cM2MyDMNAREREpJxzsHcBIiIiIiWBQpGIiIgICkUiIiIigEKRiIiICKBQJCIiIgIoFImIiIgACkUiIiIiADjZu4CSzGw2c+bMGSpVqoTJZLJ3OSIiIlIIhmFw/vx5/P39cXAo/PEfhaKbOHPmDAEBAfYuQ0RERKxw+vRpatWqVej+CkU3UalSJeDaTvX09LRzNaXHv7r/h/OZUMkTBi15yN7liIhIOZOZmUlAQIDl93hhKRTdxPVTZp6engpFtyHyH3eRfSkHF3cn7TcREbGb2536olAkNhfWM9TeJYiIiNw2XX0mIiIigkKRiIiICKDTZ1IMzp85j2E2MDmYqOR/e5PcRERuxGw2k52dbe8ypIRwdnbG0dHR5uMqFInNfTgwnswM8PSCoSvb2LscESnlsrOzSUxMxGw227sUKUG8vb3x8/Oz6X0EFYpERKTEMgyD5ORkHB0dCQgIuK0b8UnZZBgGly5dIi0tDYAaNWrYbGyFIhERKbFycnK4dOkS/v7+uLu727scKSHc3NwASEtLo3r16jY7labILSIiJVZubi4ALi4udq5ESprrIfnq1as2G1OhSERESjw9f1L+qji+EwpFIiIiImhOkYiIlEZJSXD27J3bXtWqEBh457YndqFQJCIipUtSEoSGwqVLd26b7u6QkKBgVMaVyFCUm5vL+PHjWbRoESkpKfj7+9O3b19Gjx5tOYdoGAbjxo3jww8/JD09nZYtWzJ37lzq1KljGefcuXO89NJLrFq1CgcHB6Kiopg5cyYeHh72+mgiIlJUZ89eC0SLFl0LR8UtIQF69ry23dsIRW3atKFx48bMmDHDJmX07duX9PR0VqxYYZPxJL8SGYreeust5s6dy8KFC7n33nvZu3cv/fr1w8vLi8GDBwMwdepUZs2axcKFCwkJCWHMmDFERkZy5MgRKlSoAECPHj1ITk5mw4YNXL16lX79+jFgwAAWL15sz48nIiK2EBoKTZvauwopQ0pkKNqxYwddunShU6dOAAQHB/Pvf/+b3bt3A9eOEs2YMYPRo0fTpUsXAD799FN8fX1ZsWIF3bt3JyEhgbVr17Jnzx6aNWsGwOzZs3n00Ud555138Pf3t8+HK+uSkujdxwFzroGDown27bNuHJ2/F5FSrG/fvsTGxhIbG8vMmTMBSExM5MKFCwwfPpz//Oc/VKxYkYiICKZPn07VqlUBWLZsGRMmTODEiRO4u7vTpEkTVq5cydtvv83ChQuB/7vqasuWLbRp08Yun6+sKpGh6MEHH+SDDz7g2LFj3HPPPfzwww9s27aNadOmAde+WCkpKbRv396yjpeXF82bNycuLo7u3bsTFxeHt7e3JRABtG/fHgcHB3bt2sWTTz6Zb7tZWVlkZWVZ3mdmZhbjpyyD/neev6otzvPr/L2IlGIzZ87k2LFj3HfffUycOBG49ryuBx54gOeff57p06dz+fJlRo4cyd///nc2b95McnIyzzzzDFOnTuXJJ5/k/Pnz/Oc//8EwDIYNG0ZCQgKZmZnMnz8fAB8fH3t+xDKpRIai1157jczMTOrVq4ejoyO5ublMnjyZHj16AJCSkgKAr69vnvV8fX0tbSkpKVSvXj1Pu5OTEz4+PpY+fxUTE8OECRNs/XHKD1ud57fy/L2ISEnh5eWFi4sL7u7u+Pn5ATBp0iSaNGnClClTLP0++eQTAgICOHbsGBcuXCAnJ4euXbsSFBQEQIMGDSx93dzcyMrKsowntlciQ9GXX37J559/zuLFi7n33ns5cOAAQ4YMwd/fnz59+hTbdkeNGsXQoUMt7zMzMwkICCi27ZVZOs8vIpLPDz/8wJYtW254sc/JkyeJiIigXbt2NGjQgMjISCIiInjqqaeoXLmyHaotn0pkKBo+fDivvfYa3bt3B64l5V9++YWYmBj69OljScmpqal5HgSXmppK48aNAfDz87M8LO66nJwczp07V2DKdnV1xdXVtRg+UflycFMaVxMScK7gSIOoe+xdjohIiXDhwgU6d+7MW2+9la+tRo0aODo6smHDBnbs2MH69euZPXs2r7/+Ort27SIkJMQOFZc/JTIUXbp0Kd+TkB0dHTGbzQCEhITg5+fHpk2bLCEoMzOTXbt2MXDgQADCw8NJT08nPj6esLAwADZv3ozZbKZ58+Z37sOUFra4EVpCAkkEsPTfOVw0pVLRA66GWBGKEtyoSgA6cSYipZmLi4vl2W0ATZs25auvviI4OBgnpxv/+jWZTLRs2ZKWLVsyduxYgoKCWL58OUOHDs03ntheiQxFnTt3ZvLkyQQGBnLvvfeyf/9+pk2bxrPPPgtc+9IMGTKESZMmUadOHcsl+f7+/jzxxBMAhIaG0qFDB/r378+8efO4evUqgwYNonv37rry7K9sdCO0JAIIJYFm+/dgAgzgjTBrRgrFnQQSkn9WMBKRgiUklOjtBAcHs2vXLk6dOoWHhwfR0dF8+OGHPPPMM4wYMQIfHx9OnDjBkiVL+Oijj9i7dy+bNm0iIiKC6tWrs2vXLn777TdC/zdHMzg4mHXr1nH06FGqVKmCl5cXzs7Otvyk5V6JDEWzZ89mzJgxvPjii6SlpeHv788LL7zA2LFjLX1GjBjBxYsXGTBgAOnp6bRq1Yq1a9da7lEE8PnnnzNo0CDatWtnuXnjrFmz7PGRSjYbTZA+m+DGpZ4VadoEHExQ0QOmT7/9cRK+S6TnmBDOpjspFIlIflWrXrtCtWfPO7dNd/dr270Nw4YNo0+fPtSvX5/Lly+TmJjI9u3bGTlyJBEREWRlZREUFESHDh1wcHDA09OT77//nhkzZpCZmUlQUBDvvvsuHTt2BKB///5s3bqVZs2aceHCBV2SXwxMhmEY9i6ipMrMzMTLy4uMjAw8PT3tXU7x2bcPwsIgPr5IE6SvDzPm4a04mMDTC4aubHP743yeQFjPUOIXJdC0xx24W62IlFhXrlwhMTGRkJCQPP/Tq2efSYHfDaz//V0ijxSJiIjcVGCgQorYnMOtu4iIiIiUfQpFIiIiIigUiYiIiAAKRSIiIiKAJlpLMXCr6ICrsxmPSsrcIiJSeigUic1FTnpYjz4TEZFSR/8rLyIiIoKOFMn/JBHA2QS3Io1xp+64LyKiezf+n+DgYIYMGcKQIUPsXUqpp1AkJCU7E0oCl3pWLPJYVtwJX0TkttjocY23xd392v/4ldRgJLahUCScTXfiEhVZ9EYioY+GFGmsqlXhhw93s/NCDm4eTnR+4wEbVSkico2NHtdYaAkJ1x6zdvZs8YWi7OxsXFxcimdwKTSFIrEIDblikwnSy368RGYGeHplF30wEZEChIYW6XGNxapNmzbcd999AHz22Wc4OzszcOBAJk6ciMlkIjg4mOeee47jx4+zYsUKunbtyoIFC9i2bRujRo1i7969VK1alSeffJKYmBgqVrx2JD8tLY3nnnuOjRs34ufnx6RJk+z5McscTbQWEREpBgsXLsTJyYndu3czc+ZMpk2bxkcffWRpf+edd2jUqBH79+9nzJgxnDx5kg4dOhAVFcWPP/7IF198wbZt2xg0aJBlnb59+3L69Gm2bNnCsmXLeO+990hLS7PHxyuTdKRIRESkGAQEBDB9+nRMJhN169bl4MGDTJ8+nf79+wPQtm1bXn31VUv/559/nh49elgmTNepU4dZs2bRunVr5s6dS1JSEmvWrGH37t3cf//9AHz88ceE3olziOWEjhSJiIgUgxYtWmAymSzvw8PDOX78OLm5uQA0a9YsT/8ffviBBQsW4OHhYXlFRkZiNptJTEwkISEBJycnwsLCLOvUq1cPb2/vO/J5ygMdKRIREbGD6/OErrtw4QIvvPACgwcPztc3MDCQY8eO3anSyi2FIhERkWKwa9euPO937txJnTp1cHR0vGH/pk2bcuTIEWrXrn3D9nr16pGTk0N8fLzl9NnRo0dJT0+3ad3lmU6fiYiIFIOkpCSGDh3K0aNH+fe//83s2bN5+eWXC+w/cuRIduzYwaBBgzhw4ADHjx9n5cqVlonWdevWpUOHDrzwwgvs2rWL+Ph4nn/+edzcinbjXfk/OlIkIiKl0p26i7612+nduzeXL1/mgQcewNHRkZdffpkBAwYU2L9hw4bExsby+uuv89BDD2EYBnfffTdPP/20pc/8+fN5/vnnad26Nb6+vkyaNIkxY8ZYV6Dko1AkIiKlStWq1+4w3bPnndumNXfrd3Z2ZsaMGcydOzdf26lTp264zv3338/69esLHNPPz4/Vq1fnWdarV6/bK0wKpFAkNndfi0pcuZBDBQ99vUTE9gIDrx290bPPxNb0W0tsLmJU2K07iYgUQWCgQorYnkKRiIiIjW3dutXeJYgVdPWZiIiICApFIiIiIoBOn0kx+Fe3WM5nGlTyNDFoaWt7lyMiIlIoOlIkNpedbZCVde1PERGR0kKhSERERASFIhERERFAoUhERESs1KZNG4YMGWJ5HxwczIwZM+xWT1FporWIiIiNtWnThsaNG9s0IPTt25f09HRWrFhhszFtbc+ePVSsWNHeZVhNoUhERERsolq1avYuoUh0+kxERMSG+vbtS2xsLDNnzsRkMmEymSwPgD106BAdO3bEw8MDX19fevXqxdk/PcRt2bJlNGjQADc3N6pUqUL79u25ePEi48ePZ+HChaxcudIyZkF3zW7Tpg0vvfQSQ4YMoXLlyvj6+vLhhx9y8eJF+vXrR6VKlahduzZr1qzJs96tart48SK9e/fGw8ODGjVq8O677+bb9l9Pn02bNo0GDRpQsWJFAgICePHFF7lw4YKlfcGCBXh7e7Nu3TpCQ0Px8PCgQ4cOJCcnW7Hni65EhqLg4GDLP/qfX9HR0QBcuXKF6OhoqlSpgoeHB1FRUaSmpuYZIykpiU6dOuHu7k716tUZPnw4OTk59vg4IiJSjsycOZPw8HD69+9PcnIyycnJBAQEkJ6eTtu2bWnSpAl79+5l7dq1pKam8ve//x2A5ORknnnmGZ599lkSEhLYunUrXbt2xTAMhg0bxt///ndLYEhOTubBBx8ssIaFCxdStWpVdu/ezUsvvcTAgQPp1q0bDz74IPv27SMiIoJevXpx6dIlgFvWBjB8+HBiY2NZuXIl69evZ+vWrezbt++m+8LBwYFZs2Zx+PBhFi5cyObNmxkxYkSePpcuXeKdd97hs88+4/vvvycpKYlhw4ZZu/uLpESePtuzZw+5ubmW94cOHeKRRx6hW7duALzyyit8++23LF26FC8vLwYNGkTXrl3Zvn07ALm5uXTq1Ak/Pz927NhBcnIyvXv3xtnZmSlTptjlM4mIiA3FxV173UqNGvDMM3mX/fvfUJgjEeHh1163ycvLCxcXF9zd3fHz87Ms/9e//kWTJk3y/B765JNPCAgI4NixY1y4cIGcnBy6du1KUFAQAA0aNLD0dXNzIysrK8+YBWnUqBGjR48GYNSoUbz55ptUrVqV/v37AzB27Fjmzp3Ljz/+SIsWLW5Zm7+/Px9//DGLFi2iXbt2wLXgVatWrZvW8ddJ2JMmTeKf//wn7733nmX51atXmTdvHnfffTcAgwYNYuLEibf8jMWhRIaiv56TfPPNN7n77rtp3bo1GRkZfPzxxyxevJi2bdsCMH/+fEJDQ9m5cyctWrRg/fr1HDlyhI0bN+Lr60vjxo154403GDlyJOPHj8fFxcUeH6vceOzFIK5eycW5gqO9SxGRsiorCzIzb93P0zP/sosXC7duVtbt13UTP/zwA1u2bMHDwyNf28mTJ4mIiKBdu3Y0aNCAyMhIIiIieOqpp6hcufJtb6thw4aWvzs6OlKlSpU8AcvX1xeAtLS0QtV2+fJlsrOzad68uWW5j48PdevWvWkdGzduJCYmhp9++onMzExycnK4cuUKly5dwt3dHQB3d3dLIAKoUaOGpa47rUSGoj/Lzs5m0aJFDB06FJPJRHx8PFevXqV9+/aWPvXq1SMwMJC4uDhatGhBXFwcDRo0sPyjA0RGRjJw4EAOHz5MkyZNbritrKwssv70H0FmYf6jkXzuiQyxdwkiUta5ut448PzVja6EqlixcOu6ut5+XTdx4cIFOnfuzFtvvZWvrUaNGjg6OrJhwwZ27NjB+vXrmT17Nq+//jq7du0iJOT2fq46OzvneW8ymfIsM5lMAJjN5kLVduLEidvaPsCpU6d47LHHGDhwIJMnT8bHx4dt27bx3HPPkZ2dbQlFN6rVMOzzRIQSH4pWrFhBeno6ffv2BSAlJQUXFxe8vb3z9PP19SUlJcXS58+B6Hr79baCxMTEMGHCBNsVLyIixcPKU1tA/tNpxcDFxSXPNBCApk2b8tVXXxEcHIyT041//ZpMJlq2bEnLli0ZO3YsQUFBLF++nKFDh95wTFu5VW133303zs7O7Nq1i8DAQAD++OMPjh07RuvWN37GZXx8PGazmXfffRcHh2tTmL/88stiqd9WSuRE6z/7+OOP6dixI/7+/sW+rVGjRpGRkWF5nT59uti3KSIiZU9wcDC7du3i1KlTnD17FrPZTHR0NOfOneOZZ55hz549nDx5knXr1tGvXz9yc3PZtWsXU6ZMYe/evSQlJfH111/z22+/ERoaahnzxx9/5OjRo5w9e5arV6/arN5b1ebh4cFzzz3H8OHD2bx5M4cOHaJv376WsHMjtWvX5urVq8yePZuff/6Zzz77jHnz5tms5uJQokPRL7/8wsaNG3n++ecty/z8/MjOziY9PT1P39TUVMvkMz8/v3xXo11/f7MJaq6urnh6euZ5ye07sy+F07vOcGZfwUflRETKsmHDhuHo6Ej9+vWpVq0aSUlJ+Pv7s337dnJzc4mIiKBBgwYMGTIEb29vHBwc8PT05Pvvv+fRRx/lnnvuYfTo0bz77rt07NgRgP79+1O3bl2aNWtGtWrVLBcX2cKtagN4++23eeihh+jcuTPt27enVatWhIWFFThmo0aNmDZtGm+99Rb33Xcfn3/+OTExMTaruVgYJdi4ceMMPz8/4+rVq5Zl6enphrOzs7Fs2TLLsp9++skAjLi4OMMwDOO7774zHBwcjNTUVEuf999/3/D09DSuXLlS6O1nZGQYgJGRkWGDT1NyxS86YsC1P23h3ce3GONabzHefXxLiahHREqvy5cvG0eOHDEuX75s71KkhLnZd8Pa398ldk6R2Wxm/vz59OnTJ8/5TS8vL5577jmGDh2Kj48Pnp6evPTSS4SHh9OiRQsAIiIiqF+/Pr169WLq1KmkpKQwevRooqOjcbXxxDkREREpG0psKNq4cSNJSUk8++yz+dqmT5+Og4MDUVFRZGVlERkZmeeeB46OjqxevZqBAwcSHh5OxYoV6dOnj93ueyAiIiIlX4kNRREREQVeklehQgXmzJnDnDlzClw/KCiI7777rrjKExERkTKmRE+0FhEREblTFIpEREREUCgSERERARSKRERERACFIhERERFAoUhEREQEUCiSYhD9WTijvgkn+jMrH9YoIiKlSnBwMDNmzLC8N5lMrFixwm71WKvE3qdISi9XT901XETKtzZt2tC4ceM8QaE8SU5OpnLlyvYu47YpFImIiNiBYRjk5ubmeZRVWXGzh6+XZDp9JiIiYkN9+/YlNjaWmTNnYjKZMJlMnDp1iq1bt2IymVizZg1hYWG4urqybds2+vbtyxNPPJFnjCFDhtCmTRvLe7PZTExMDCEhIbi5udGoUSOWLVt20zqCg4OZNGkSvXv3xsPDg6CgIL755ht+++03unTpgoeHBw0bNmTv3r151tu2bRsPPfQQbm5uBAQEMHjwYC5evGhpT0tLo3Pnzri5uRESEsLnn3+eb9t/PX02cuRI7rnnHtzd3bnrrrsYM2YMV69etbSPHz+exo0b89lnnxEcHIyXlxfdu3fn/PnzhdjjtqNQJDYX98FBts44QNwHB+1diojIHTdz5kzCw8Pp378/ycnJJCcnExAQYGl/7bXXePPNN0lISKBhw4aFGjMmJoZPP/2UefPmcfjwYV555RV69uxJbGzsTdebPn06LVu2ZP/+/XTq1IlevXrRu3dvevbsyb59+7j77rvp3bu35bFaJ0+epEOHDkRFRfHjjz/yxRdfsG3bNgYNGmQZs2/fvpw+fZotW7awbNky3nvvPdLS0m5aR6VKlViwYAFHjhxh5syZfPjhh0yfPj1Pn5MnT7JixQpWr17N6tWriY2N5c033yzU/rGVsnfMTuwu7tvfycwATy8IH2DvakSkLIqLu/a6lRo14Jln8i77978hOfnW64aHX3vdLi8vL1xcXHB3d7/haaSJEyfyyCOPFHq8rKwspkyZwsaNGwn/X0F33XUX27Zt4/3336d169YFrvvoo4/ywgsvADB27Fjmzp3L/fffT7du3YBrR3DCw8NJTU3Fz8+PmJgYevTowZAhQwCoU6cOs2bNonXr1sydO5ekpCTWrFnD7t27uf/++wH4+OOPCQ0NvelnGD16tOXvwcHBDBs2jCVLljBixAjLcrPZzIIFC6hUqRIAvXr1YtOmTUyePLnQ+6qoFIpERKTUycqCzMxb9/P0zL/s4sXCrZuVdft1FUazZs1uq/+JEye4dOlSviCVnZ1NkyZNbrrun49E+fr6AtCgQYN8y9LS0vDz8+OHH37gxx9/zHNKzDAMzGYziYmJHDt2DCcnJ8LCwizt9erVw9vb+6Z1fPHFF8yaNYuTJ09y4cIFcnJy8PzLP05wcLAlEAHUqFHjlkegbE2hSERESh1X1xsHnr+qWPHGywqzrmsxXUhb8S9FOTg4WE5fXffn+TYXLlwA4Ntvv6VmzZp/qfHmRTo7O1v+bjKZClxmNpst23rhhRcYPHhwvrECAwM5duzYTbd3I3FxcfTo0YMJEyYQGRmJl5cXS5Ys4d133y2w1uu1Xa/rTlEoEhGRUsfaU1uQ/3RacXBxcSE3N7dQfatVq8ahQ4fyLDtw4IAlJNSvXx9XV1eSkpJueqrMFpo2bcqRI0eoXbv2Ddvr1atHTk4O8fHxltNnR48eJT09vcAxd+zYQVBQEK+//rpl2S+//GLTum1FE61FRERsLDg4mF27dnHq1CnOnj170yMebdu2Ze/evXz66accP36ccePG5QlJlSpVYtiwYbzyyissXLiQkydPsm/fPmbPns3ChQttWvfIkSPZsWMHgwYN4sCBAxw/fpyVK1daJlrXrVuXDh068MILL7Br1y7i4+N5/vnncXNzK3DMOnXqkJSUxJIlSzh58iSzZs1i+fLlNq3bVhSKREREbGzYsGE4OjpSv359qlWrRlJSUoF9IyMjGTNmDCNGjOD+++/n/Pnz9O7dO0+fN954gzFjxhATE0NoaCgdOnTg22+/JSQkxKZ1N2zYkNjYWI4dO8ZDDz1EkyZNGDt2LP7+/pY+8+fPx9/fn9atW9O1a1cGDBhA9erVCxzz8ccf55VXXmHQoEE0btyYHTt2MGbMGJvWbSsm468nMsUiMzMTLy8vMjIy8k0IK0v2fZ5AWM9Q4hcl0LTHza8gKIxpXbZarj4burKN3esRkdLrypUrJCYmEhISQoUKFexdjpQgN/tuWPv7W0eKRERERLAyFOXm5pKZmUlOTk6e5ZcvX2bChAk8+eSTvPLKK5w5c8YmRYqIiIgUN6uuPps4cSKTJk1i69atPPTQQ8C1+xi0adOGvXv3YhgGJpOJr7/+mgMHDpTKh8KJ9WoEuuCZmUtFT0d7lyIiIlJoVoWiTZs24efnZwlEAKtWrWLPnj3cc889vPjii6xZs4b169fz4Ycf5rljpZR9z8x+0N4liIiI3DarTp8lJiZSr169PMtWrlyJyWTi888/Z/DgwaxatYpq1ard8oF1IiIit6JrguSviuM7YVUo+v333/M9z2X79u3UrFnTcutvJycnWrRocdPLEEVERG7G0fHaafjs7Gw7VyIlzaVLl4D8d8IuCqtOnzk5OXHx4kXL+z/++IPjx4/z97//PU+/SpUqkZGRUbQKRUSk3HJycsLd3Z3ffvsNZ2dnHBx00XR5ZxgGly5dIi0tDW9vb0twtgWrQtFdd93Fzp07MZvNODg4sHr1agzDoFWrVnn6paWlUa1aNZsUKqXHv1/awcX/TbTW/CIRKQqTyUSNGjVITEwssY+GEPvw9vbOd9aqqKwKRY8//jhTpkyhS5cutG/fnrfeegtHR0c6d+5s6WMYBvv37yc0VDffK2+Sk7L/d/PGwj33R0TkZlxcXKhTp45OoYmFs7OzTY8QXWdVKBoxYgQrV67k22+/5dtvvwXgtddeIzAw0NJn27ZtnD17Nt/RIxERkdvl4OCgO1pLsbMqFHl6erJ7926WLVtGamoq999/f74n9/7++++8/PLLPP300zYpVERERKQ4WRWKANzc3OjVq1eB7U888QRPPPGEtcOLiIiI3FE2mcafnZ1NcnIy586ds8VwIiIiIndckULRokWLeOCBB6hYsSK1atVi2LBhlrbly5fzj3/8g8TExCIXKSIiIlLcrA5Fzz//PH369GHv3r24ubnlu7PkPffcw5IlS/jqq6+KXKSIiIhIcbMqFH3++ed88skn3HfffezZs+eGN2i89957qVWrFmvWrLGqsF9//ZWePXtSpUoV3NzcaNCgAXv37rW0G4bB2LFjqVGjBm5ubrRv357jx4/nGePcuXP06NEDT09PvL29ee6557hw4YJV9YiIiEjZZlUo+uCDD/Dw8GD16tWEhYVhMplu2K9BgwZWnT77448/aNmyJc7OzqxZs4YjR47w7rvvUrlyZUufqVOnMmvWLObNm8euXbuoWLEikZGRXLlyxdKnR48eHD58mA0bNrB69Wq+//57BgwYcPsfWERERMo8q64+++GHH2jevDkBAQE37efj40Nqauptj//WW28REBDA/PnzLctCQkIsfzcMgxkzZjB69Gi6dOkCwKeffoqvry8rVqyge/fuJCQksHbtWvbs2UOzZs0AmD17No8++ijvvPMO/v7+t12XFE54pypkXcrF1d32N9YSEREpLlYdKcrKysLLy+uW/X777Ter7jj5zTff0KxZM7p160b16tVp0qQJH374oaU9MTGRlJQU2rdvb1nm5eVF8+bNiYuLAyAuLg5vb29LIAJo3749Dg4O7Nq1q8DPlZmZmeclty98QAPaDGlM+IAG9i5FRESk0KwKRTVr1iQhIeGmfQzD4MiRI3mO8BTWzz//zNy5c6lTpw7r1q1j4MCBDB48mIULFwKQkpICgK+vb571fH19LW0pKSlUr149T7uTkxM+Pj6WPn8VExODl5eX5XWrI2EiIiJSdlgVitq1a8dPP/3EypUrC+zz2Wef8d///pdHHnnktsc3m800bdqUKVOm0KRJEwYMGED//v2ZN2+eNeUW2qhRo8jIyLC8Tp8+XazbExERkZLDqlA0bNgwXF1d+cc//sGMGTM4c+aMpe3cuXPMmzePF198kYoVKzJ48ODbHr9GjRrUr18/z7LQ0FCSkpIALE/F/et8pdTUVEubn58faWlpedpzcnI4d+5cgU/VdXV1xdPTM89Lbl9WZpblJSIiUlpYFYrq1KnDwoULMZvNvPrqqwQEBGAymVi4cCHVqlUjOjqanJwcFixYkOchsYXVsmVLjh49mmfZsWPHCAoKAq5Nuvbz82PTpk2W9szMTHbt2kV4eDgA4eHhpKenEx8fb+mzefNmzGYzzZs3t+ZjSyHN6RVHzONxzOkVZ+9SRERECs3qmzd269aNPXv20K1bNypVqoRhGBiGQYUKFejcuTNxcXFERUVZNfYrr7zCzp07mTJlCidOnGDx4sV88MEHREdHA2AymRgyZAiTJk3im2++4eDBg/Tu3Rt/f3/L89ZCQ0Pp0KED/fv3Z/fu3Wzfvp1BgwbRvXt3XXkmIiIi+Vj9QFiA++67jyVLlmAYBr///jtms5mqVavi4FC0R6rdf//9LF++nFGjRjFx4kRCQkKYMWMGPXr0sPQZMWIEFy9eZMCAAaSnp9OqVSvWrl1LhQoVLH0+//xzBg0aRLt27XBwcCAqKopZs2YVqTYREREpm4oUiq4zmUxUrVrVFkNZPPbYYzz22GM33ebEiROZOHFigX18fHxYvHixTesSERGRsqnIoSg7O5v4+Hh+/fVX4Nrl+mFhYbi4uBS5OBEREZE7xepQlJOTw4QJE5g9ezbnz5/P01apUiUGDx7M2LFjcXKyycEoERERkWJlVWIxm808/vjjrFu3DsMwqFy5suUmjYmJifzxxx9MnjyZ+Ph4Vq1aVeQ5RiIiIiLFzaq08tFHH7F27VqCgoJYtmwZv//+O3v37mXv3r38/vvvfPXVVwQFBbF27Vo+/vhjW9csIiIiYnNWhaJPP/0UNzc3Nm/eTNeuXfO1P/nkk2zatAlXV1fLozlERERESjKrQtGhQ4do06YNwcHBBfYJCQmhbdu2HDp0yNraRERERO4Yq+YUZWVl4eXldct+lSpVIitLj3oob7qPq0fuVTOOzppLJiIipYdVoSggIIC4uDhyc3NxdHS8YZ/c3Fx27txJrVq1ilSglD7+TW/8bDkREZGSzKr/lY+MjCQpKYmXX36Zq1ev5mvPzs5m8ODBJCUl0bFjxyIXKSIiIlLcrDpS9Nprr7F48WLmzp3LypUr6d69u+WS/J9//pkvvviCM2fO4OPjw8iRI21asIiIiEhxsCoU1axZk7Vr19KtWzeSkpKYNm1annbDMAgMDGTZsmXUrFnTJoVK6XFsXSJXr+TiXMGReyJD7F2OiIhIoVh9u+n777+fY8eOsXTpUrZu3ZrnMR9t2rShW7duetRHObX6vV/IzABPLxiqUCQiIqVEkZ7B4eLiQo8ePfI8vV5ERESkNNI10yIiIiIoFImIiIgAhTx9dtddd1m9AZPJxMmTJ61eX0REROROKFQoOnXqlNUbMJlMVq8rIiIicqcUKhQlJiYWdx0iIiIidlWoUBQUFFTcdYiIiIjYlSZai4iIiKBQJMXAxcWEq+u1P0VEREqLIt28ce/evSxbtoyjR4+SmZmJYRj5+phMJjZt2lSUzUgpM2hpa3uXICIictusDkXDhg1j+vTpliBkMpnyhKLr73X1mYiIiJQGVp0+W7p0KdOmTaNmzZq8//77REREALBu3Tr+9a9/ER4ejmEYvPbaa2zevNmmBYuIiIgUB6tC0QcffICjoyObNm2if//+1KhRA4BHHnmEF198ke3bt/P6668zbdo0vLy8bFqwiIiISHGwKhTt37+f5s2bU6dOnQL7TJgwgRo1ajBp0iSri5PSaX1MPN+8vov1MfH2LkVERKTQrJpTdP78eQIDAy3vXVxcALhw4QIeHh4AODg40Lx5c2JjY21QppQmh3aeJzMDPL0gwt7FiIiIFJJVR4qqVatGenq65X3VqlWB/I8DuXjxIpmZmVYXJyIiInKnWBWKgoOD+eWXXyzvmzRpgmEYLF682LIsJSWF2NhY3Q1bRERESgWrQlG7du346aefLEeGOnbsiI+PD2+99RbdunXj1VdfpXnz5ly8eJGoqChb1isiIiJSLKyaU9S9e3fOnDnD6dOnCQ4OpmLFisyfP5/u3bvz1VdfWfqFhYUxatQomxUrIiIiUlysCkWhoaF8+OGHeZZ17tyZ48ePs2rVKs6dO0doaCidO3fG0dHRJoWKiIiIFKciPebjr/z9/XnhhRdsOaSIiIjIHVEiHwg7fvx4TCZTnle9evUs7VeuXCE6OpoqVarg4eFBVFQUqampecZISkqiU6dOuLu7U716dYYPH05OTs6d/igiIiJSShT5SNGvv/7Kr7/+ypUrVwrs8/DDD9/2uPfeey8bN260vHdy+r9SX3nlFb799luWLl2Kl5cXgwYNomvXrmzfvh2A3NxcOnXqhJ+fHzt27CA5OZnevXvj7OzMlClTbrsWERERKfusDkUrV67ktdde49ixYzftZzKZrDpC4+TkhJ+fX77lGRkZfPzxxyxevJi2bdsCMH/+fEJDQ9m5cyctWrRg/fr1HDlyhI0bN+Lr60vjxo154403GDlyJOPHj7fcbFKKR52G7ly+kIObh03PzoqIiBQrq35rrVmzhqioKMxmM15eXtx11114enratLDjx4/j7+9PhQoVCA8PJyYmhsDAQOLj47l69Srt27e39K1Xrx6BgYHExcXRokUL4uLiaNCgAb6+vpY+kZGRDBw4kMOHD9OkSZMbbjMrK4usrCzLe9140jqd33jA3iWIiIjcNqtC0eTJkzGbzYwfP57XXnvN5kdemjdvzoIFC6hbty7JyclMmDCBhx56iEOHDpGSkoKLiwve3t551vH19SUlJQW4duPIPwei6+3X2woSExPDhAkTbPpZREREpHSwKhQdOHCAxo0bM3bsWFvXA1y7GeR1DRs2pHnz5gQFBfHll1/i5uZWLNsEGDVqFEOHDrW8z8zMJCAgoNi2JyIiIiWHVVefOTo65rkarLh5e3tzzz33cOLECfz8/MjOzs7z7DWA1NRUyxwkPz+/fFejXX9/o3lK17m6uuLp6ZnnJSIiIuWDVaGoYcOG/Pe//7V1LQW6cOECJ0+epEaNGoSFheHs7MymTZss7UePHiUpKYnw8HAAwsPDOXjwIGlpaZY+GzZswNPTk/r169+xusurD3p+z7QuW/mg5/f2LkVERKTQrApFQ4YMYfv27ezdu9fW9QAwbNgwYmNjOXXqFDt27ODJJ5/E0dGRZ555Bi8vL5577jmGDh3Kli1biI+Pp1+/foSHh9OiRQsAIiIiqF+/Pr169eKHH35g3bp1jB49mujoaFxdXYulZvk/F86bycy49qeIiEhpYVUoioqKYsyYMURGRvLee++RlJRk06L++9//8swzz1C3bl3+/ve/U6VKFXbu3Em1atUAmD59Oo899hhRUVE8/PDD+Pn58fXXX1vWd3R0ZPXq1Tg6OhIeHk7Pnj3p3bs3EydOtGmdIiIiUnYUaqL1zZ5f9tJLL/HSSy8V2G7NfYqWLFly0/YKFSowZ84c5syZU2CfoKAgvvvuu9varoiIiJRfhQpFhmFYvYGirCsiIiJypxQqFJnNmhsiIiIiZVuJfCCsiIiIyJ2mUCQiIiKCQpGIiIgIoFAkIiIiAlj57DORm3mktz9Xr+TiXKHgWzmIiIiUNApFYnMNou6xdwkiIiK3TafPRERERChkKPr+++85duxYcdciIiIiYjeFCkVt2rThzTfftLxv27YtU6dOLbaipHQ7e/R30g7/xtmjv9u7FBERkUIr9JyiPz+uY+vWrQQHBxdHPVIGfDriIJkZ4OkFQ1e2sXc5IiIihVKoI0WVKlUiOTm5uGsRERERsZtCHSlq2LAhmzdvZuzYsdSuXRuAEydO8OmnnxZqI71797a+QhEREZE7oFChaMSIETz11FNMnjzZsmz79u1s3769UBtRKBIREZGSrlChqHPnzuzevZsVK1bwyy+/sGDBAu6++25atmxZ3PWJiIiI3BGFnmjdqFEjGjVqBMCCBQto1aoVn3zySbEVJiIiInInWXVH63HjxtGkSRNb1yIiIiJiN1aHIhEREZGypEjPPsvJyWHZsmVs2bKFX3/9FYCaNWvyt7/9jaeeegonJz1aTUREREoHq1PLgQMHeOqpp0hMTMxzY0eAjz76iDFjxrB06VIaN25c1BpFREREip1VoejMmTNERERw9uxZfH196d69O3fffTcAP//8M0uWLOHkyZNERkZy4MABatSoYdOipWTrPzcMw2xgcjDZuxQREZFCsyoUvfXWW5w9e5bnn3+emTNn4ubmlqd9ypQpDB48mI8++oipU6cyffp0mxQrpUMl/0r2LkFEROS2FeoxH3+1Zs0aAgMDmTt3br5ABFChQgXee+89AgMD+fbbb4tcpIiIiEhxsyoUnT59mgcffBBHR8cC+zg5OREeHs7p06etLk5ERETkTrHq9JmrqyuZmZm37Hf+/HlcXV2t2YSUYvGLEsi+lIOLuxNhPUPtXY6IiEihWBWK6tevz5YtWzh9+jQBAQE37JOUlMSWLVt09Vk5FLs0lcwM8PRCoUhEREoNq06f9e7dm8uXL9O+fXu+++67fO2rV6/mkUce4cqVK3oYrIiIiJQKVh0p6t+/P1999RWbNm2ic+fO+Pj4EBISAkBiYiLnzp3DMAzat29P//79bVqwiIiISHGw6kiRo6Mj3377LSNGjKBixYr8/vvv7N27l7179/L7779TsWJFRo4cyerVq3FwsGoTIiIiIneU1Xe0dnFx4c0332TChAns3bs3z2M+mjVrpgnWIiIiUqoU+eFkrq6utGzZ0ha1iIiIiNiNzm2JiIiIUApC0ZtvvonJZGLIkCGWZVeuXCE6OpoqVarg4eFBVFQUqampedZLSkqiU6dOuLu7U716dYYPH05OTs4drl5ERERKixIdivbs2cP7779Pw4YN8yx/5ZVXWLVqFUuXLiU2NpYzZ87QtWtXS3tubi6dOnUiOzubHTt2sHDhQhYsWMDYsWPv9EcQERGRUqLEhqILFy7Qo0cPPvzwQypXrmxZnpGRwccff8y0adNo27YtYWFhzJ8/nx07drBz504A1q9fz5EjR1i0aBGNGzemY8eOvPHGG8yZM4fs7OwCt5mVlUVmZmael9y+KtWdqFbdRJXqRZ6yJiIicseU2FAUHR1Np06daN++fZ7l8fHxXL16Nc/yevXqERgYSFxcHABxcXE0aNAAX19fS5/IyEgyMzM5fPhwgduMiYnBy8vL8irobt1yc30+bEX0l63p82Ere5ciIiJSaCUyFC1ZsoR9+/YRExOTry0lJQUXFxe8vb3zLPf19SUlJcXS58+B6Hr79baCjBo1ioyMDMtLD7MVEREpP6w6v5GUlITJZCqWIymnT5/m5ZdfZsOGDVSoUMHm49+Mq6ur7q8kIiJSTll1pCg4OJju3bvbuhbg2umxtLQ0mjZtipOTE05OTsTGxjJr1iycnJzw9fUlOzub9PT0POulpqbi5+cHgJ+fX76r0a6/v95HRERE5M+sCkWenp6WZ53ZWrt27Th48CAHDhywvJo1a0aPHj0sf3d2dmbTpk2WdY4ePUpSUhLh4eEAhIeHc/DgQdLS0ix9NmzYgKenJ/Xr1y+WuuX/fDU8js/+uY2vhsfZuxQREZFCs+r0Wf369Yttvk2lSpW477778iyrWLEiVapUsSx/7rnnGDp0KD4+Pnh6evLSSy8RHh5OixYtAIiIiKB+/fr06tWLqVOnkpKSwujRo4mOjtbpsTvgl2NZZGaAp5fuCyUiIqWHVUeK+vfvz/bt29mzZ4+t6ymU6dOn89hjjxEVFcXDDz+Mn58fX3/9taXd0dGR1atX4+joSHh4OD179qR3795MnDjRLvWKiIhIyWfVkaJ+/fqxf/9+IiIiGD58OFFRUQQHBxfbUZitW7fmeV+hQgXmzJnDnDlzClwnKCiI7777rljqERERkbLHqlDk6Oho+fuYMWMYM2ZMgX1NJpMeryEiIiIlnlWhyDCMYukrIiIiYi9WhSKz2WzrOkRERETsqkTe0VpERETkTlMoEhEREaGIoejkyZOMGDGCVq1aUbduXUaMGGFp27VrFx988AEZGRlFLlJERESkuFk1pwhg4cKF/POf/yQrKwu4dpXZ2bNnLe2XLl1i4MCBuLi40Ldv3yIXKqVH0795k3UxB9eKVn+9RERE7jirfmvt3LmT559/Hnd3d9544w1at25N8+bN8/Rp3bo1Xl5erFq1SqGonGkzpLFNxklIrAD7ijZG1aoQGGiTckREpIyzKhRNnToVwzD49ttvadWq1Q37ODg40LhxY44cOVKkAqX8qeqdgzsX6TkmBAq+BVahuLtDQoKCkYiI3JpVoWj79u088MADBQai6/z8/Ni3r4j/qy/lTmCNqyQQytlF6yA01OpxEhKgZ084e1ahSEREbs2qUJSenk5gIX7LXL58mezsbGs2IeVcIKcJDL0MTe1diYiIlBdWhaIqVarwyy+/3LLfiRMn8PPzs2YTUopN67KVzAzw9IKhK9vYuxwREZFCsSoUtWjRglWrVnH48GHuvffeG/bZvn07hw8fpmfPnkUqUMqxhIQiru8GWH/6TUREyherQlF0dDQrVqwgKiqKJUuW0Lhx4zztCQkJPPvss5hMJl588UVb1CkFSEq6NmemKBISK9imGFupWvXaDOkiB+omwD5ITgZq2KAwEREpy6wKRe3atWPo0KFMmzaNsLAw7r77bkwmE+vWraNhw4YcOXIEs9nMiBEjaNGiha1rlv9JSro2D/nSpaKOFII7F6nqnWOLsoouMPDaUaKipr3vUq5dvZaejkKRiIjcitV313vnnXeoW7cu48eP58SJEwAkJyeTnJxM1apVGTduHNHR0TYrVPI7e/ZaIFq0qEgXaUFCAlV7RhJYY4WtSiu6wMCiXzJW1NNvIiJSrhTplsP9+/fn+eefZ//+/fz888+YzWYCAgK4//77cXLS3YzvlNBQaFqkq7QuA6dtVI2IiEjpVOTkYjKZaNq0KU2L9ltZRERExK5scjjHMAx+//13DMOgSpUqODgU6TmzIiIiIndckdLLhg0b6NChA5UqVcLX1xc/Pz8qVapEhw4dWLduna1qFBERESl2Voei4cOH06FDB9avX8+lS5cwDAPDMLh8+TLr16/n0Ucf5dVXX7VlrSIiIiLFxqpQtGjRIt59910qVKjAq6++yo8//sj58+c5f/48Bw8eZNiwYbi5uTFjxgwWLVpk65qlhOs6vDY9x4TQdXhte5ciIiJSaFbNKZo9ezaOjo6sXbuWhx56KE/bvffey9SpU3n88cdp06YN//rXv3RX63ImuFUte5cgIiJy26w6UnTo0CFatWqVLxD92fX2Q4cOWV2ciIiIyJ1iVSiqUKEC/v7+t+zn7++Pi4uLNZsQERERuaOsOn0WFhbGjz/+eMt+P/74I82aNbNmE3I7EhK4dgPGoqxvO6e2/ZecrFycXB11Kk1EREoNq0LR66+/Trt27Zg6dSojRoy4YZ+3336bhIQEZs+eXaQC5SauP+i0Zw9gf9HGcne/9iBWG/j67RNkZoCnFwxVKBIRkVKiUKHo+++/z/PeZDIxaNAgRo0axdKlS+nVqxchISEAJCYmsmjRIuLj4xk8eLBu5Ficrj/o9I1J8Khf0caqWrXozxoTEREpxQoVitq0aYPJZMq33DAM4uPj2bdvX77lALNmzWL27Nnk5JSQp6+XVSEh0LQoT4QVERGRQoWihx9++IahSERERKSsKFQo2rp1azGXISIiImJfmvAjIiIigkKRiIiICGDlJfnXXblyhb1793LmzBmuXLlSYL/evXvf1rhz585l7ty5nDp1Crj26JCxY8fSsWNHy3ZfffVVlixZQlZWFpGRkbz33nv4+vpaxkhKSmLgwIFs2bIFDw8P+vTpQ0xMDE5ORfrIIiIiUkZZnRDefvttpkyZQmZm5i373m4oqlWrFm+++SZ16tTBMAwWLlxIly5d2L9/P/feey+vvPIK3377LUuXLsXLy4tBgwbRtWtXtm/fDkBubi6dOnXCz8+PHTt2kJycTO/evXF2dmbKlClWfV4REREp26wKRf/6178YOXIkAA0aNKBOnTpUqlTJZkV17tw5z/vJkyczd+5cdu7cSa1atfj4449ZvHgxbdu2BWD+/PmEhoayc+dOWrRowfr16zly5AgbN27E19eXxo0b88YbbzBy5EjGjx9f4KNHsrKyyMrKsrwvTOATERGRssHqUOTk5MRXX32VL8DYWm5uLkuXLuXixYuEh4cTHx/P1atXad++vaVPvXr1CAwMJC4ujhYtWhAXF0eDBg3ynE6LjIxk4MCBHD58mCZNmtxwWzExMUyYMKFYP095MHRlG3uXICIictusmmh96tQpHn744WINRAcPHsTDwwNXV1f++c9/snz5curXr09KSgouLi54e3vn6e/r60tKSgoAKSkpeQLR9fbrbQUZNWoUGRkZltfp06dt+6FERESkxLLqSFH16tWpVq2arWvJo27duhw4cICMjAyWLVtGnz59iI2NLdZturq64urqWqzbEBERkZLJqiNFHTt2JC4uDrPZbOt6LFxcXKhduzZhYWHExMTQqFEjZs6ciZ+fH9nZ2aSnp+fpn5qaip/fted/+fn5kZqamq/9epuIiIjIX1kVisaNG0d2djaDBw8mOzvb1jXdkNlsJisri7CwMJydndm0aZOl7ejRoyQlJREeHg5AeHg4Bw8eJC0tzdJnw4YNeHp6Ur9+/TtSb3m2dcYB1k3ey9YZB+xdioiISKFZdfrM39+fbdu28fjjj1O3bl3+9re/ERgYiIND/oxlMpkYM2bMbY0/atQoOnbsSGBgIOfPn2fx4sVs3bqVdevW4eXlxXPPPcfQoUPx8fHB09OTl156ifDwcFq0aAFAREQE9evXp1evXkydOpWUlBRGjx5NdHS0To/dAfu2pJOZAZ5e0GaIvasREREpHKtCkWEYzJw5k59++gmz2cyCBQvy9TGZTBiGYVUoSktLo3fv3iQnJ+Pl5UXDhg1Zt24djzzyCADTp0/HwcGBqKioPDdvvM7R0ZHVq1czcOBAwsPDqVixIn369GHixInWfFwREREpB6wKRW+//TazZ8/GycmJxx57jDp16uDh4WGzoj7++OObtleoUIE5c+YwZ86cAvsEBQXx3Xff2awmERERKdusCkUfffQR7u7u/Oc//ynwnj8iIiIipYlVE61Pnz7NQw89pEAkIiIiZYZVocjPz8+mj/UQERERsTerQtGTTz7Jf/7zH65cuWLrekRERETswqpQNH78eHx8fHjmmWc4e/asrWsSERERueOsmmg9ZMgQ6taty4oVK9i8eTNhYWE3vU/Rra4mExEREbE3q0LRggULMJlMAJw/f56tW7cW2FehqPwJuseVS+dzca/kaO9SRERECs2qUDR//nxb1yFlSNTb4fYuQURE5LZZFYr69Olj6zpERERE7MqqidYiIiIiZY1CkYiIiAhWnj579tlnC91XE63Ln4X9t3EhIxcPL0f6fNjK3uWIiIgUitVXn93M9SvTDMNQKCqHfk/LITMDsrJy7F2KiIhIodn06jOz2cwvv/zCd999x969exkyZAiNGjUqUoEiIiIid0KxXH02fvx4RowYwYcffsi+ffusKkxERETkTiq2idZTpkyhUqVKjB07trg2ISIiImIzxRaKnJycaNq0KRs3biyuTYiIiIjYTLFekn/58mX++OOP4tyEiIiIiE0UWyhKSEhg27ZtBAQEFNcmRERERGzGqonWn376aYFt58+fJyEhgc8++4wrV67wj3/8w+riRERERO4Uq0JR3759LfciuhHDMADo0qULo0ePtq4yERERkTvIqlDUu3fvAkORi4sLNWvWpH379jz44INFKk5Kp9bdfMm+lIOLu1VfLxEREbsoljtaS/kW1jPU3iWIiIjcNj0QVkRERASFIhERERGgkKfPbna1WWH07t27SOtL6XL+zHkMs4HJwUQl/0r2LkdERKRQChWKbnW12a0oFJUvHw6MJzMDPL1g6Mo29i5HRESkUAoVitq2bXvboSguLo5Lly4VKUyJiIiI3CmFCkW38/yy//znP4wYMYLLly8D0KBBA+sqExEREbmDbDbR+tChQ3Tu3Jk2bdqwa9cuAgICWLBgAfv377fVJkRERESKTZHvrnf69GnGjBnD559/Tm5uLlWqVOH//b//R3R0NC4uLraoUURERKTYWR2K/vjjDyZPnsx7773HlStXcHd35+WXX2bkyJF4enraskYRERGRYnfbp8+uXLlCTEwMd911F9OnTycnJ4cBAwZw4sQJJk+ebJNAFBMTw/3330+lSpWoXr06TzzxBEePHs1XR3R0NFWqVMHDw4OoqChSU1Pz9ElKSqJTp064u7tTvXp1hg8fTk5OTpHrExERkbKn0KHIbDbzwQcfULt2bUaPHk1GRgZPPvkkhw4dYt68efj5+dmsqNjYWKKjo9m5cycbNmzg6tWrREREcPHiRUufV155hVWrVrF06VJiY2M5c+YMXbt2tbTn5ubSqVMnsrOz2bFjBwsXLmTBggWMHTvWZnWKiIhI2VGo02dff/01r7/+OseOHcMwDFq3bs1bb73FAw88UCxFrV27Ns/7BQsWUL16deLj43n44YfJyMjg448/ZvHixbRt2xaA+fPnExoays6dO2nRogXr16/nyJEjbNy4EV9fXxo3bswbb7zByJEjGT9+vOY7iYiISB6FCkVPPfUUJpPJMm/o0UcfJScnhx07dhRqIw8++GCRiszIyADAx8cHgPj4eK5evUr79u0tferVq0dgYCBxcXG0aNGCuLg4GjRogK+vr6VPZGQkAwcO5PDhwzRp0iTfdrKyssjKyrK8z8zMLFLdIiIiUnrc1kTrS5cuERMTQ0xMTKHXMZlMRZrHYzabGTJkCC1btuS+++4DICUlBRcXF7y9vfP09fX1JSUlxdLnz4Hoevv1thuJiYlhwoQJVtcq1/Se2gBzjhkHJz1aT0RESo9ChaLAwEC73Zk6OjqaQ4cOsW3btmLf1qhRoxg6dKjlfWZmJgEBAcW+3bKmat0q9i5BRETkthUqFJ06daqYy7ixQYMGsXr1ar7//ntq1aplWe7n50d2djbp6el5jhalpqZaJnz7+fmxe/fuPONdvzqtoEnhrq6uuLq62vhTiIiISGlQIs9vGIbBoEGDWL58OZs3byYkJCRPe1hYGM7OzmzatMmy7OjRoyQlJREeHg5AeHg4Bw8eJC0tzdJnw4YNeHp6Ur9+/TvzQURERKTUKPIdrYtDdHQ0ixcvZuXKlVSqVMkyB8jLyws3Nze8vLx47rnnGDp0KD4+Pnh6evLSSy8RHh5OixYtAIiIiKB+/fr06tWLqVOnkpKSwujRo4mOjtbRoGJ28KtjXL2Si3MFRxpE3WPvckRERAqlRIaiuXPnAtCmTZs8y+fPn0/fvn0BmD59Og4ODkRFRZGVlUVkZCTvvfeepa+joyOrV69m4MCBhIeHU7FiRfr06cPEiRPv1McotzZ8eobMDPD0QqFIRERKjRIZigzDuGWfChUqMGfOHObMmVNgn6CgIL777jtbliYiIiJlVImcUyQiIiJypykUiYiIiKBQJCIiIgIoFImIiIgACkUiIiIigEKRiIiICKBQJCIiIgKU0PsUSenmUckBMP/vTxERkdJBoUhsbsCih+1dgoiIyG3T/8qLiIiIoFAkIiIiAigUiYiIiACaUyTFYNWY3Vy+kIObhxOd33jA3uWIiIgUikKR2NzxHy+RmQGeXtn2LkVERKTQdPpMREREBIUiEREREUChSERERARQKBIREREBFIpEREREAIUiEREREUChSERERARQKBIREREBdPNGKQb3tajElQs5VPDQ10tEREoP/dYSm4sYFWbvEkRERG6bTp+JiIiIoFAkIiIiAigUiYiIiACaUyTF4F/dYjmfaVDJ08Sgpa3tXY6IiEih6EiR2Fx2tkFW1rU/RURESguFIhEREREUikREREQAhSIRERERQKFIREREBCihoej777+nc+fO+Pv7YzKZWLFiRZ52wzAYO3YsNWrUwM3Njfbt23P8+PE8fc6dO0ePHj3w9PTE29ub5557jgsXLtzBTyEiIiKlSYkMRRcvXqRRo0bMmTPnhu1Tp05l1qxZzJs3j127dlGxYkUiIyO5cuWKpU+PHj04fPgwGzZsYPXq1Xz//fcMGDDgTn0EERERKWVK5H2KOnbsSMeOHW/YZhgGM2bMYPTo0XTp0gWATz/9FF9fX1asWEH37t1JSEhg7dq17Nmzh2bNmgEwe/ZsHn30Ud555x38/f3v2GcRERGR0qFEhqKbSUxMJCUlhfbt21uWeXl50bx5c+Li4ujevTtxcXF4e3tbAhFA+/btcXBwYNeuXTz55JM3HDsrK4usrCzL+8zMzGL7HElxv3L256KNn7D9nI2qERERkVIXilJSUgDw9fXNs9zX19fSlpKSQvXq1fO0Ozk54ePjY+lzIzExMUyYMMHGFeeXFPcroQ96c4maRR7LnYtUvcvTBlXZzmMvBnH1Si7OFRztXYqIiEihlbpQVJxGjRrF0KFDLe8zMzMJCAiw+XbO/pzJJWqyaOB2Qlv6FGmsqnd5Ehhe9HBlS/dEhti7hLwSE2Hf5aKNUbUqBAbaph4RESmRSl0o8vPzAyA1NZUaNWpYlqemptK4cWNLn7S0tDzr5eTkcO7cOcv6N+Lq6oqrq6vtiy5AaEsfmvYIvWPbK3e8va/9OWY0jNlftLHc3SEhQcFIRKQMK3WhKCQkBD8/PzZt2mQJQZmZmezatYuBAwcCEB4eTnp6OvHx8YSFhQGwefNmzGYzzZs3t1fpcqddD82LPofQIhwpSkiAnj3h7FmFIhGRMqxEhqILFy5w4sQJy/vExEQOHDiAj48PgYGBDBkyhEmTJlGnTh1CQkIYM2YM/v7+PPHEEwCEhobSoUMH+vfvz7x587h69SqDBg2ie/fuuvLsDjizL4Xcq2YcnR3wb1rwkbk7JYGiHo1zoyoBKA6JiJRtJTIU7d27l7/97W+W99fn+fTp04cFCxYwYsQILl68yIABA0hPT6dVq1asXbuWChUqWNb5/PPPGTRoEO3atcPBwYGoqChmzZp1xz9LebRkwk9kZoCnFwxdab9QVLXqtbNePXsWdaRQ3EkgIflnBSMRkTKsRIaiNm3aYBhGge0mk4mJEycyceLEAvv4+PiwePHi4ihPSonAwGtnvs6eLdo4Cd8l0nNMCGfTnRSKRETKsBIZikRsJTDQBtOAEq7cuo+IiJR6JfIxHyIiIiJ3mkKRiIiICApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoDuUyTFIPqzcHuXICIictsUisTmXD1d7V2CiIjIbdPpMxEREREUikREREQAnT6TYhD3wUGyLuXi6u5I+IAG9i5HRESkUBSKxObivv2dzAzw9ILwAfauRkREpHB0+kxEREQEhSIRERERQKFIREREBFAoEhEREQEUikREREQAhSIRERERQKFIREREBFAoEhEREQF080YpBjUCXfDMzKWip6O9SxERESk0hSKxuWdmP2jvEkRERG6bTp+JiIiIoFAkIiIiAigUiYiIiACaUyTF4N8v7eDi/yZal6X5RQmJFWBf0caoWhUCA21Tj4iI2JZCkdhcclI2mRng6ZVr71Jsoqp3Du5cpOeYEBhTtLHc3cwk/OSgYCQiUgIpFIncQmADLxIqNOXslYpFGieBUHpe/pyzB5MJDKxho+pERMRWFIpEbiUwkMCjGwg8e7Zo43yXcu1IU3o6oFAkIlLSKBSJFEZgYNEnAyUk2KYWEREpFgpFIndYwvZzQNECUtW7PAkMr2mbgkREBCgHoWjOnDm8/fbbpKSk0KhRI2bPns0DDzxg77KkHKp6l+e1CdtzW8Lcoo3lzkUSdvyqYCQiYkNlOhR98cUXDB06lHnz5tG8eXNmzJhBZGQkR48epXr16vYuT8qZwPCaJOz4lbM/JxVpnITt5+g5tyX/2ZRGqGvR69JtAkRErinToWjatGn079+ffv36ATBv3jy+/fZbPvnkE1577TU7VyflUWB4zSIf3anqfRD3uba5RQCAe4Vcvn77Z6pVzin6YLbg7Q01ij4RvWrWrwS6phZ5nKQsX866Fv2InK3qsRmlYSmKpCQo6sUnUOK+h2U2FGVnZxMfH8+oUaMsyxwcHGjfvj1xcXE3XCcrK4usrCzL+4yMDAAyMzNtWtuFSxeATC5cumDzsUuCK1cvkpUDV67aft8JeAc7sNu1Eb9nFe0WAQBnqUrPK4vo8JKvDSqzpaJ/b9xwYBHDqYr1P7jPUpWeLOJyCanHplwrwBuTroVQkduRng5jRkPWlSIP5VchE7+9qyEgoOh1/cn13z2GYdzeikYZ9euvvxqAsWPHjjzLhw8fbjzwwAM3XGfcuHEGoJdeeumll156lYHX6dOnbys7lNkjRdYYNWoUQ4cOtbw3m82cO3eOKlWqYDKZbmuszMxMAgICOH36NJ6enrYutdTR/shL+yMv7Y+8tD/y0v74P9oXeRW0PwzD4Pz58/j7+9/WeGU2FFWtWhVHR0dSU/Oew09NTcXPz++G67i6uuLqmnfmqncRDy17enrqi/sn2h95aX/kpf2Rl/ZHXtof/0f7Iq8b7Q8vL6/bHsfBVgWVNC4uLoSFhbFp0ybLMrPZzKZNmwgPD7djZSIiIlISldkjRQBDhw6lT58+NGvWjAceeIAZM2Zw8eJFy9VoIiIiIteV6VD09NNP89tvvzF27FhSUlJo3Lgxa9euxde3+K+0cXV1Zdy4cflOx5VX2h95aX/kpf2Rl/ZHXtof/0f7Ii9b7w+TYdzu9WoiIiIiZU+ZnVMkIiIicjsUikRERERQKBIREREBFIpEREREAIWiYjNnzhyCg4OpUKECzZs3Z/fu3fYu6Y74/vvv6dy5M/7+/phMJlasWJGn3TAMxo4dS40aNXBzc6N9+/YcP37cPsUWs5iYGO6//34qVapE9erVeeKJJzh69GiePleuXCE6OpoqVarg4eFBVFRUvhuOlhVz586lYcOGlpushYeHs2bNGkt7edoXN/Lmm29iMpkYMmSIZVl52ifjx4/HZDLledWrV8/SXp72xXW//vorPXv2pEqVKri5udGgQQP27t1raS9PP0+Dg4PzfT9MJhPR0dGA7b4fCkXF4IsvvmDo0KGMGzeOffv20ahRIyIjI0lLS7N3acXu4sWLNGrUiDlz5tywferUqcyaNYt58+axa9cuKlasSGRkJFeuFP3BgiVNbGws0dHR7Ny5kw0bNnD16lUiIiK4ePGipc8rr7zCqlWrWLp0KbGxsZw5c4auXbvaseriU6tWLd58803i4+PZu3cvbdu2pUuXLhw+fBgoX/vir/bs2cP7779Pw4YN8ywvb/vk3nvvJTk52fLatm2bpa287Ys//viDli1b4uzszJo1azhy5AjvvvsulStXtvQpTz9P9+zZk+e7sWHDBgC6desG2PD7Yc3DVuXmHnjgASM6OtryPjc31/D39zdiYmLsWNWdBxjLly+3vDebzYafn5/x9ttvW5alp6cbrq6uxr///W87VHhnpaWlGYARGxtrGMa1z+7s7GwsXbrU0ichIcEAjLi4OHuVeUdVrlzZ+Oijj8r1vjh//rxRp04dY8OGDUbr1q2Nl19+2TCM8vf9GDdunNGoUaMbtpW3fWEYhjFy5EijVatWBbaX95+nL7/8snH33XcbZrPZpt8PHSmysezsbOLj42nfvr1lmYODA+3btycuLs6OldlfYmIiKSkpefaNl5cXzZs3Lxf7JiMjAwAfHx8A4uPjuXr1ap79Ua9ePQIDA8v8/sjNzWXJkiVcvHiR8PDwcr0voqOj6dSpU57PDuXz+3H8+HH8/f2566676NGjB0lJSUD53BfffPMNzZo1o1u3blSvXp0mTZrw4YcfWtrL88/T7OxsFi1axLPPPovJZLLp90OhyMbOnj1Lbm5uvrtm+/r6kpKSYqeqSobrn7887huz2cyQIUNo2bIl9913H3Btf7i4uOR76HBZ3h8HDx7Ew8MDV1dX/vnPf7J8+XLq169fLvcFwJIlS9i3bx8xMTH52srbPmnevDkLFixg7dq1zJ07l8TERB566CHOnz9f7vYFwM8//8zcuXOpU6cO69atY+DAgQwePJiFCxcC5fvn6YoVK0hPT6dv376Abf9bKdOP+RApKaKjozl06FCeORLlUd26dTlw4AAZGRksW7aMPn36EBsba++y7OL06dO8/PLLbNiwgQoVKti7HLvr2LGj5e8NGzakefPmBAUF8eWXX+Lm5mbHyuzDbDbTrFkzpkyZAkCTJk04dOgQ8+bNo0+fPnauzr4+/vhjOnbsiL+/v83H1pEiG6tatSqOjo75Zr2npqbi5+dnp6pKhuufv7ztm0GDBrF69Wq2bNlCrVq1LMv9/PzIzs4mPT09T/+yvD9cXFyoXbs2YWFhxMTE0KhRI2bOnFku90V8fDxpaWk0bdoUJycnnJyciI2NZdasWTg5OeHr61vu9smfeXt7c88993DixIly+f2oUaMG9evXz7MsNDTUckqxvP48/eWXX9i4cSPPP/+8ZZktvx8KRTbm4uJCWFgYmzZtsiwzm81s2rSJ8PBwO1ZmfyEhIfj5+eXZN5mZmezatatM7hvDMBg0aBDLly9n8+bNhISE5GkPCwvD2dk5z/44evQoSUlJZXJ/3IjZbCYrK6tc7ot27dpx8OBBDhw4YHk1a9aMHj16WP5e3vbJn124cIGTJ09So0aNcvn9aNmyZb5beBw7doygoCCg/P08vW7+/PlUr16dTp06WZbZ9Pth4wnhYhjGkiVLDFdXV2PBggXGkSNHjAEDBhje3t5GSkqKvUsrdufPnzf2799v7N+/3wCMadOmGfv37zd++eUXwzAM48033zS8vb2NlStXGj/++KPRpUsXIyQkxLh8+bKdK7e9gQMHGl5eXsbWrVuN5ORky+vSpUuWPv/85z+NwMBAY/PmzcbevXuN8PBwIzw83I5VF5/XXnvNiI2NNRITE40ff/zReO211wyTyWSsX7/eMIzytS8K8uerzwyjfO2TV1991di6dauRmJhobN++3Wjfvr1RtWpVIy0tzTCM8rUvDMMwdu/ebTg5ORmTJ082jh8/bnz++eeGu7u7sWjRIkuf8vTz1DCuXckdGBhojBw5Ml+brb4fCkXFZPbs2UZgYKDh4uJiPPDAA8bOnTvtXdIdsWXLFgPI9+rTp49hGNcuIx0zZozh6+truLq6Gu3atTOOHj1q36KLyY32A2DMnz/f0ufy5cvGiy++aFSuXNlwd3c3nnzySSM5Odl+RRejZ5991ggKCjJcXFyMatWqGe3atbMEIsMoX/uiIH8NReVpnzz99NNGjRo1DBcXF6NmzZrG008/bZw4ccLSXp72xXWrVq0y7rvvPsPV1dWoV6+e8cEHH+RpL08/Tw3DMNatW2cAN/yMtvp+mAzDMIpwJEtERESkTNCcIhEREREUikREREQAhSIRERERQKFIREREBFAoEhEREQEUikREREQAhSIRERERQKFIREREBFAoEimTNmzYQL9+/bjnnnvw9PTE1dWVGjVq8MgjjzB9+nR+++03e5coVtq6dSsmk4k2bdrYuxSRMkehSKQMOXv2LI888ggREREsWLCAq1ev8re//Y2oqChCQ0PZsWMHQ4cO5a677mLXrl32LrdE6tu3LyaTiQULFthl+yaTCZPJZJdti5R3TvYuQERsIyMjg1atWnH06FHq1avHBx98wEMPPZSnT1ZWFgsXLmTcuHEkJyfbqVIpigceeICEhATc3d3tXYpImaNQJFJGvPTSSxw9epTg4GC2b9+Oj49Pvj6urq4MGDCALl26kJ6efueLlCJzd3enXr169i5DpEzS6TORMuDnn39m8eLFAEybNu2GgejPfH19qVu3bp5lS5YsoV27dvj4+ODq6kpQUBDPPvssx44du+EYwcHBmEwmTp06xZYtW4iIiKBy5cq4ubnRtGlTPv300wK3bxgGX3/9NY899hh+fn64uLjg5+dHq1ateOutt7h8+XK+deLj4+nRoweBgYG4urri4+NDZGQk3333nU3qO3XqFCaTiYULFwLQr18/y6ksk8nE+PHjLX3/fIpr/vz5hIeH4+XlZdkewC+//MJbb71F27ZtLTV7e3vTqlUr3n//fcxmc57tjx8/Ps9psz9v+8/j3mpO0U8//US/fv0ICgqy7Kd27drx5Zdf3rD/9e2OHz+e3377jejoaAICAnBxcSEgIICXXnpJAVrKD0NESr2ZM2cagOHt7W3k5OTc1rpms9no3bu3ARhOTk5G27Ztje7duxv33HOPARju7u7GmjVr8q0XFBRkAMaYMWMMk8lkhIWFGd27dzdatGhhAAZgTJ8+Pd962dnZRteuXQ3AcHBwMFq0aGE888wzxiOPPGLUrFnTAIzExMQ868yYMcNwcHAwAKNx48bGU089ZbRq1cpwcXExAGPChAlFru+3334z+vTpY9x9990GYLRs2dLo06eP5bV8+XJL3+vrDxo0yHBwcDBatWplPPPMM0bz5s2NU6dOGYZhGG+88YYBGCEhIUa7du2M7t27G61bt7bU3LVrV8NsNlvGXL58udGnTx/L2H/edp8+fYzffvvNMAzD2LJliwEYrVu3zveZV69ebVSoUMEAjLp16xrdu3c32rZtazg6OhqA8eyzz+ZbZ9y4cZa2WrVqGb6+vkbXrl2NRx991PDy8jIA4/777zeys7PzrStS1igUiZQBvXr1MgCjbdu2t73u3LlzDcCoWrWqsX//fstys9ls+YXp7e1tpKWl5VnveuhwdnY2Vq1aladt/vz5BmB4eXkZly5dytM2dOhQAzCCg4ONAwcO5Gkzm83Gxo0bjfT0dMuytWvXGiaTyahataoRGxubp/+PP/5o1KpVywCMrVu32qS+68Fk/vz5Be6z68HF09PTiIuLu2Gf3bt3GwcPHsy3/NdffzUaNWpkAMaXX35Z4NgFKSgUpaSkWELMpEmT8gSuPXv2GJUrVzYA44MPPsiz3vV/Y8Do27evceXKFUtbUlKSJaguXry4wJpEygqFIpEyoEOHDgZgdO/e/bbXvX5kZNasWfnazGaz0bBhQwMwJk+enKfteugYOnToDcetV6+eARjff/+9ZVlqaqrlSMnevXsLVV/z5s0NwFi2bNkN27/88ksDMKKioopcn2HcXiiaOHFioT7DX61bt84AjG7duhU4dkEKCkXXj0yFhYXdcL133nnHAIw6derkWX49FNWqVcu4ePFivvXefPPNAo8yiZQ1mlMkUo7997//5eTJkwD06dMnX7vJZKJfv34AbNmy5YZjdO7c+YbLQ0NDAfj1118ty7Zs2UJ2djZhYWGEhYXdsr6zZ8+ye/du3NzcCtzO9bk1O3bsKHJ9t+upp566aXtWVharVq1i7Nix/POf/6Rfv3707duX999/H4CjR49ave2/2rp1K3Djf0eA5557DoDjx49z5syZfO3t2rW74RVttthPIqWFrj4TKQOqVasGQFpa2m2td/0XXZUqVfD09Lxhn7vvvjtP378KDAy84fLr4125csWy7JdffgEo9NVTiYmJGIbB5cuXcXV1vWnfgm5IeTv13a7g4OAC23bu3MnTTz9NUlJSgX0yMzOt3vZfXf/3CQkJuWG7t7c3Pj4+nDt3jv/+97/4+/vnaS/O/SRSWigUiZQBYWFhfPbZZ+zbt4/c3FwcHR3v2LYdHIrvgPP1K7Q8PDyIioqyaozirM/Nze2Gyy9dusQTTzxBamoq/fr1Y+DAgdSuXRtPT08cHR05duwYdevWxTCMYqvtdhXnfhIpLRSKRMqAxx57jKFDh5Kens4333zDk08+Waj1atasCcDvv/9OZmbmDY8W/fzzz3n6FsX1oxE//fRTofoHBAQA107jffLJJ6XmF/f3339PamoqTZs25ZNPPsnXfvz4cZtvs2bNmvz000+Wf6+/ysjI4Ny5c5a+IpJf6fgJIyI3dffdd/PMM88A8Oqrr1p++RUkLS2No0ePUqtWLcvpsRs91sIwDMvyv/3tb0Wus23btri4uBAfH8++fftu2d/f35+GDRty/vx51q5dW+TtF4aLiwsAOTk5Vo9xff8XdEpq0aJFBa7r7Oxs1favz626fp+lv7oezurUqaNQJFIAhSKRMmL27NnUrl2bxMREWrVqxbZt2/L1yc7O5pNPPqFJkyYkJCQAMGzYMADeeOMNfvjhB0tfwzCYNGkSBw4cwNvbm/79+xe5xurVqzNw4EAAunXrxqFDh/K0G4bB5s2bycjIsCybNGkScO1miqtWrco3pmEY7Nq1i/Xr1xe5PoBatWoBcPjwYavHuD45edOmTRw5ciRP2wcffMAXX3xh8+33798fT09P9u3bx5QpU/Kcmtu/f79lPw4fPvy2xhUpT3T6TKSMqFy5Mtu3b+fpp59m69atPPTQQ4SEhNCwYUPc3d1JTU1l9+7dXLhwAU9PT8tE2xdeeIEdO3bw2Wef0axZM1q3bk316tXZt28fR48exc3NjcWLF1smcxfV1KlTSUxM5JtvvqFRo0Y0b96ckJAQzp49y+HDh/n1119JTEzEy8sLuHb12MyZM3n11Vd5/PHHqV27NnXr1sXLy4vffvuNH374gbS0NEaOHElERESR63viiSeYMGECs2bN4tChQwQEBODg4MDjjz/O448/XqgxmjRpQpcuXVi5ciVNmjShTZs2+Pj4cODAAY4ePcr/+3//j8mTJ99w3aioKN555x3at29P27ZtqVSpEgBvvfUWVapUKXCbvr6+fP7553Tr1o3XX3+dzz77jCZNmpCWlkZsbCw5OTn069fPJuFWpMyy4+0ARKSYrFmzxujdu7dRu3Ztw8PDw3B2djb8/PyMRx55xJgxY4bx+++/51tn8eLFRps2bQxvb2/D2dnZCAgIMPr27Wv89NNPN9zG9fsA/fXu09fd7H4/ZrPZWLx4sREREWFUqVLFUt9DDz1kvP3228bly5fzrXPw4EFjwIABRp06dYwKFSoY7u7uxl133WVERkYas2bNMn799Veb1bd8+XKjZcuWRqVKlQyTyWQAxrhx4yzt3OJeQoZx7c7db7/9ttGgQQPD3d3d8PHxMSIiIoz169cbiYmJBmAEBQXlW+/y5cvGiBEjjNq1a1vu6fTnz3GzO1obhmEcOXLE6NOnj1GrVi3D2dnZ8Pb2Nv72t78ZS5YsuWH/6/cp+vPn+7NbbU+kLDEZRgm6/EFERETETjSnSERERASFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSERERAeD/A5GMa7zHkfB4AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "_, bins, _ = plt.hist(cvir_test, label='test', histtype='step', bins=25, color='r')\n",
    "plt.hist(cvir_pred, label='pred', histtype='step', bins=bins, color='b')\n",
    "plt.axvline(np.mean(cvir_test),  ls='--', label='test median', alpha=0.5, lw=2, color='r')\n",
    "plt.axvline(np.mean(cvir_pred), ls='--', label='true median', alpha=0.5, lw=2, color='b')\n",
    "plt.xlabel('Concentration', fontsize=16)\n",
    "plt.ylabel('Number of haloes', fontsize=16)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also examine their spearman correlation coefficient as in the paper."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.7841946739737133"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from multicam.correlations import spearmanr\n",
    "\n",
    "spearmanr(cvir_test, cvir_pred)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multivariable prediction \n",
    "\n",
    "MultiCAM enables the prediction of multiple targets at once. Let's try to predict concentration, virial ratio, centre-of-mass offset, spin, and triaxility all at once. \n",
    "\n",
    "For this we create new model, using the same class as before. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_features = ma.shape[1]\n",
    "n_targets = 5 # <---- note the change since we are predicting multiple targets\n",
    "mv_model = LinearRegression(n_features, n_targets, use_multicam=True)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now to carefully stack the corresponding variables to use. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((10000, 165), (10000, 5))"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = ma \n",
    "y = np.vstack([cvir, xoff, t_over_u, spin, c_to_a]).T # just need to remember the order of the targets\n",
    "x.shape, y.shape"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and doing the train/test split again... "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_idx, test_idx = get_tt_indices(n_haloes, test_ratio=0.25)\n",
    "x_train = x[train_idx]\n",
    "x_test = x[test_idx]\n",
    "y_train = y[train_idx]\n",
    "y_test = y[test_idx]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit and predict as before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "mv_model.fit(x_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = mv_model.predict(x_test)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute spearman correlation coefficients:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cvir: 0.8022848241735719\n",
      "xoff: 0.6351004548960729\n",
      "t_over_u: 0.7500029134232726\n",
      "spin: 0.4891788779504358\n",
      "c_to_a: 0.4948907902178588\n"
     ]
    }
   ],
   "source": [
    "# compute spearman correlation for each target\n",
    "target_names = ['cvir', 'xoff', 't_over_u', 'spin', 'c_to_a']\n",
    "for i in range(n_targets):\n",
    "    print(f'{target_names[i]}:', spearmanr(y_test[:, i], y_pred[:, i]))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, more complicated analysis and or plots can be done from `y_pred` and `y_test` as desired."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "nbody-relaxed-P7T5KfjX-py3.9",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.16"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}