{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "unlimited-saying",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "demographic-final",
   "metadata": {},
   "source": [
    "## Load the block of simulation data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "mediterranean-throat",
   "metadata": {},
   "outputs": [],
   "source": [
    "from load_tng_parametric_scatter_data import load_tng_mah_data\n",
    "res = load_tng_mah_data()\n",
    "halo_ids, log_mah_500c, log_mah_200m, log_mah_mvir, tng_t, log_mah_fit_min = res\n",
    "n_halos = halo_ids.size"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cardiovascular-medicare",
   "metadata": {},
   "source": [
    "## Load the best-fit parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "outstanding-mayor",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import h5py\n",
    "from collections import OrderedDict\n",
    "\n",
    "drn = \"/Users/aphearin/work/random/0112/FIT_DATA\"\n",
    "\n",
    "fit_data_200m = OrderedDict()\n",
    "with h5py.File(os.path.join(drn, \"M200mean_Msun.h5\"), 'r') as hdf:\n",
    "    for key in hdf.keys():\n",
    "        fit_data_200m[key] = hdf[key][...]\n",
    "\n",
    "fit_data_500c = OrderedDict()\n",
    "with h5py.File(os.path.join(drn, \"M500crit_Msun.h5\"), 'r') as hdf:\n",
    "    for key in hdf.keys():\n",
    "        fit_data_500c[key] = hdf[key][...]\n",
    "\n",
    "fit_data_vir = OrderedDict()\n",
    "with h5py.File(os.path.join(drn, \"Mvir_Msun.h5\"), 'r') as hdf:\n",
    "    for key in hdf.keys():\n",
    "        fit_data_vir[key] = hdf[key][...]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "minus-conversation",
   "metadata": {},
   "source": [
    "## Calculate the best-fit smooth MAHs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "residential-legislature",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:absl:No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)\n"
     ]
    }
   ],
   "source": [
    "from diffmah.individual_halo_assembly import calc_halo_history\n",
    "\n",
    "dmhdt, log_mah = calc_halo_history(tng_t, tng_t[-1], fit_data_vir['logmp_fit'], \n",
    "                  10**fit_data_vir['mah_logtc'], fit_data_vir['early_index'], fit_data_vir['late_index'])\n",
    "fit_data_vir['log_mah'] = log_mah\n",
    "\n",
    "dmhdt, log_mah = calc_halo_history(tng_t, tng_t[-1], fit_data_500c['logmp_fit'], \n",
    "                  10**fit_data_500c['mah_logtc'], fit_data_500c['early_index'], fit_data_500c['late_index'])\n",
    "fit_data_500c['log_mah'] = log_mah\n",
    "\n",
    "dmhdt, log_mah = calc_halo_history(tng_t, tng_t[-1], fit_data_200m['logmp_fit'], \n",
    "                  10**fit_data_200m['mah_logtc'], fit_data_200m['early_index'], fit_data_200m['late_index'])\n",
    "fit_data_200m['log_mah'] = log_mah"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "blond-custody",
   "metadata": {},
   "source": [
    "## Plot a randomly selected halo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ideal-survey",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuEAAAEjCAYAAACciN+wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABcXElEQVR4nO3dd1gc57XH8e/QQYUFVFFHvVgFkHsXcou7JdfYiZNYSnHKTZGiNDvlxpGTm3vTI9mJHTsusnCvsZB7F0KyegNUkYQELL3ve/+YBa8QiLbsLOzv8zz7oJ2ZnT1gn9mz77zFMsYgIiIiIiKBE+Z0ACIiIiIioUZFuIiIiIhIgKkIFxEREREJMBXhIiIiIiIBpiJcRERERCTAVISLiIiIiASYinARERERkQCLcDoACR2WZS3y/nM8kALcZYxxt9hf7H2aYoy5v5XXt7lfRNpnWdYC4CbgPsANLADcxpiVPscoF0WCiGVZywGXMWax07GI/6gIl4CwLGtRiw/5BcB67IK8+UPdGJPpfZ5iWdaKpgtOe/tFpFNSsfPPDaz0LaKViyJBaZXTAYj/WVoxU3qaZVkpwGJjzNIW20uwW8MzLctab4xJa7G/eVt7+0WkYyzLWtBUQLexX7koIhIA6hMugbKolW3FQKJlWS7s7ikn7bcsK6O9/f4LUSS0KRdFRAJH3VGkxxlj8oCEVnalANnen8Wt7Hf77Gtzv2VZqcBy7/P7gETABcw1xiz1KQ5SgbxTtQKKhAJvTriwcybVpztKt3LR5/y+fcYTfbuiichnvLmyFDuHFno/L7Esaw12Ti0G8oAVAMaY+d7PvAewPz9Xe4+bb4xZGPBfQLpFRbg4wnvhyTLG5LTTgubCLqrb3O89x3K8hbgxJsf7HjdZlrXcpxtMlrcLjIpwCWV52HnS9GFfbFnWGmPMfNrJtQ7sbxpAltuiz/giFeIiJzPGrLQsCyCtKSe9VmB/RroBLMtaiv0Zh/czr+n5cuxi3BXAsMVPVIRLwPn0EfdnH9Ji7ILc9yKW19pxlmW5fGdlEQklTV9SfZ9blpXubV3rFm93lUXGGN87XwuApO6eW6Sv8hbiJdit3r7b3T5Pff8NJ3/maYaiXkh9wsUJy4F5HTjO1cn97laOKerA+4iEujwg/RT7Xe28vml/Bi2+/Bpj7m85KFtETvJU0zS+3i+z7g68prWGJulFVIRLQHlvVS9t8Q0/m9ZvcycCOR3Y36StvqoiQnPXkJJTHOKvXBSRzlnBZy3hGcaYrA68xt1z4UggqAiXgPF+y1/h22XEsqwMb0Fe7P3278tljMlqb38PhizSF93XyrYUPut/2p1czKGV2VNaOV5EfHi7ibm83TUlRKgIl4DwDr7MblGA+/ZBXY7PNIbefVmd2A+tt9C5uh61SN/izT+37zbvwllP+eRml3PRe47m2+o+NH2hSPuWY8920tHGJVfPhSKBoMV6pMd5v9nntrE7wWf0d5eWyvYWAcuwP+jvM8bc7y0smqYtXGqMybIsa4l320pgeYtBnCIhw6dIdoHdb7uV/V1ett6ba27sPqsuTQsq0j7vHaPlLVef9X6GLsf+jFuK3S2s+TMPe9Vbd0CDFb9QES4iIiIiEmDqjiIiIiIiEmAqwkVEREREAkxFuIiIiIhIgKkIFxEREREJsJBctn7QoEFm7NixTochEnL27t2Lck8k8JR7Is5Yv379cWPM4Nb2heTsKMnJyebw4cPcc8893HvvvU6HIxIykpOTUe6JBJ5yT8QZlmWtN8akt7YvJFvCk5OTKSgocDoMkZCj3BNxhnJPJPioT7iIiIiISICpCBcRERERCTAV4SIiIiIiAaYiXEREREQkwFSEi4iIiIgEmIpwEREREZEAC8kivKCgAMuyNFeqSGfVVkBRLg1573L4vUfZmvlr1q34Buue/t8OvVy5J+IM5Z5I8NE84SJiqymD0oP2o+wQlBV4H/a/PWUFhNWVA/aFY7j3UWsi2VSd0aG3UO6JOEO5JxJ8QrIIFwlJ1W5w74OSfeDe/9mj9AC4D0Bt6YnHW2HUxw3hsElkb10iuTVjKTQJVEQNZsSoFFJSUhg1ZhzjRoxgbpQuJSIiIp2hT06RvsIYqDgKxXn2oygXSvKhOB9K9kKN+8TjowdC/ChwjYbRZ0H8SHCNgoEjIX4E7xwO5+tPbAJg7tgE5o5LZH5KErNGuggPswL+64mIiPQlKsJFepv6aijaA8d22j+P77J/FuVCXcVnx4VF2AV2wjgYkQYJYyFhDLjG2D9jXGC1Xkw/te4AP3p2AxOG9OehO+cyPD42IL+aiIhIqFARLhKsGmrtQrtwOxzbDoU77J8l+wDjPciyW6+TJsKoMyFpAiSlQOJ4u5U7vGMpXlhWwyd7i/kk337sOFLOeRMH8dfbUhkQE9ljv6KIiEhvVFPfSExkeLfOoSJcJBiUH4Ejm+HIJjiyBY5utVu3TaO9PyzSLrCT58DMm2HwJBg0GZLGQ2THW6nLa+rZeaSc/cVVHCqpJr+okpx9JewtqgIgLiqc1NEJLLt8BF86dxyR4SE5gZKIiISixnqoPA6VhVB5zPvvY95HEabyGOXFR6gtO0ZMvZuG726j/8CELr+dinCRQDLGnnGkYAMc3ggFG+Hwp3bCN3GNhqEzYOpVMHQaDJluF9vhHWuRbmj0cLi0hgPFVewrrmJvUSV5xyrZcaSMA8XVJxw7eEA0s0a6uO2MMcwdl8j05IEqvEVEpG+pr4Hyw3aDV9PPiiPen0ehotD+WVXU+uvDo6iNTmRfdSxHGvpTEZ5C0pDhTGpo7FZYKsJFelJNKRzKgUPZcHA9FOTYiQ5ghcPgKTAhA4bPhGEzYeh0iHW1fbr6Rg65qzlQXMXh0hqKK+soqazjeEUtBe4aCkqrOVxaQ6PHNL8mKjyMUYmxzBzh4ua5o5kybABjB/VjhCu227fSREREHNVYbzdulR7wTrN7AEpbTLNbXXzy68KjoP8wGDAUElNg9JnQfyj0Gwz9h0C/IdBvEPQbzMs7K/ivpz5lVGIs38mYxCXThxId0f3PTxXhIv5ijD0LyYGPYf9HmAMfQ+F2LG//7er48ZQkncnRMdMo6DeNkgGTaAyPsV/rAXPIULzrKPuK8zlaVkNFbQOVtY1U1DZQU9dIdX0jDT7FdZPYyHCS+keR7Ipl7thEkl0xjE6MY1RCHKMS40h2xWo2ExER6Z08HrvVummmr5ZT7ZYXgPGc+JrYRIgfYT9GzYWByTBguM9jGMQmtDk5QZO6Bg+Pf7yPn7+0jbTRCTz4hXRccVF++9VUhEtIq65rJGd/CXsKKygsr6GwrJaqug7eXjKGoXV7mVj9KROrNzG+6lMSPPa37XJi2eCZSHbjDeSYiWzypFBW0w+O+p4g96RThlmQ7IpleHwMQwbE0H9QBP2iI4iNDCc2Koy4qAiSXTGMTLCL66R+UWrNFhGR3q2pq2bRHijOtWf78p1mt7H2s2OtMBg4wp58YNx53ql2R9k/40fZBXdUXCfe2lBe20BJZR35xyvZfbSCXUfL2Xa4jN1HK6hr9JAxdSh/vnWO3z9vVYRLr1Tf6GFrQRnr95Ww+2g5tQ0e6ho8J3TDaM/xilo+PeimvtF+TXiYxeD+0fSLDsdq49vxMM8RUhs+ZXbDJuY0biLB2AvcHLOS2BI9k4MDZnEkfjblAycQGx1FXFQEl8VGclNsJPGxkSTEReGKiyQuqvX36B8dQVSE+mSLiEgf1FhvF9fHdtjT6x7f5Z1uNxfqK5sPMxGxeFxjqXelUD/qIiriRlEUNYKj4UMpDBtMeb1FRW0jJZV1FB2upWh3HZV1DVTWHqS6bh8eY2irGjDGLrwbjaGx0VBZ10DL0mFQ/2imDh/Al84dx8yR8VwybSgRPTBeKiSL8IKCAizL4p577uHee+91OhzpoF1Hy1m7vZAPco+TvbeE6nq7xXpQ/yhio8KJCg8jIqzjSdI/JoIvnTuOM1OSmJ48kEH9oglr2W2jthzy34E9WZD7JpTle188DFIug7HnwdhzGJwwjsHt3NYS5Z6IU5R7vV9VXQO7jlawp7CC/OMV1NR72n9RABgDHmPwGEODx+DxGOoaGulXc4QhVbtJrs1nRF0eI+r3Mrz+IBE0NL/2WPhQDoSNZLfnQrY0DGWPZzh7PcM4QgKmIgwOtny3Yu/DlhAXSVL/aBLjohgyIIZ+gyKIiQgjPMzy9jRp/XM5zLIb3sIsi/7REbji7Iay0YlxTBo6gIR+/utyciohWYQnJydTUFDgdBjSQfuKKvntf3by0qbDAEwa2p+b5o5i7thE0sYkMCw+xn9vZgwc3w27/wO7/gP7PwJPPUT1twvuM78GKRfCoEnt9iWTkyn3RJyh3Oud9hRW8PzGQ3yYW3TCnduIMCtouiKG4WGcVcAMK59p7GUKeUwy+xjIZ4vHHQ0bwv7wMXwSO5f9EaPZHz6GI5EjaQiPJToinGHxMQwdGM05keGcHxZGZLhFeJhFZLj974ExkST2iyKpfxQDYiLpFx1BXGT4yQ1nvUxIFuES/IwxfJJfTOb6gzy74RCR4WF86+IJfP7MMQwZ6MeiG8DTaA+m3PEy7HzFXvIdYMg0OOvrMGE+jDoDIgLzzVhERELb/qIq/m/tLp7bcAjLspgxIp4vn5vC7FEuJg7tz+jEOGemkzXG7qd9KAcOrbcfRzZDfZW9hlxEjD3F7rCFMGwGDD0NhkxlaMxAhgJzAx9xUFMRLkGlodHDE+sO8MA7eewvrqJfVDi3nD6ab148wb/Fd0Od3c1k+wt28V113F4QZ9z5cObXYdJl9kAPERGRACiurOP1rUd4dcsR3t9znPAwiy+fO46vXjCepP7RzgRVV2UX2gc+hoPr7EfTXNoRsTB8FqR90f45fJa9enMHV2oWFeESRNZuP8qvX9lO7rFK0sYk8J2MiVw2YxhxUX7637SxAfa+A1uege0vQo3b7mYy6VKYcqU9X3fMQP+8l4iISAdU1Dbw1zf38OB7+dQ1eBiVGMuXzxvHl84Zx1B/3/ltT3WJ3Q1z3/uw7wN7MTmPtw/3oMkw6XIYmW4/Bk9Vwd1N+utJUPjrW3u4/7WdpAzuxwN3pJMxdUibM5R0ijH2bbPNT9nFd2UhRA2AKVfAtGth/MUQGeCLnIiICPDK5sPc88JWjpXXct2cEXz53HFMTx7on8+/jqitgP0fQt5bsPddOLwJMPZCNiPS4exv2YvYjJwLcYmBiSmEqAgXxz3y4V7uf20n18xO5ncLZ/mnn1vpIfj0Cfj0SSjaDeHRdov3aQth4nyIjO3+e4iIiHTRa1uO8I3HczhtRDwrbk8jdXRCz7+pxwNHPv1sxq8Dn9iTD4RHw6jT4cIfwthz7QJcDVQ9TkW4BFxNfSPbDpdRWFbLtoJS/vjGHuZPG9r9AryhDna9CjmPQu5aewWtMefA2d+Eadeccjl4ERGRQFm3t5hvPbmBWSNdPHHXmcRG9eBMJzWlsGct7H7dLr4rj9nbh820Jx9Iuchu7VbjVMCpCJeAOlJaw60PfETe8c8m5b9w8mD+dMucrhfgJfsg51928V1ZCAOS4dzvwpzbIDHFT5GLiIh0TWFZDf94P5+aukYM8PzGAka4YvnnF+f2TAHuPmDP9rXjZbt/t6fBXqZ9QoY949f4i6D/EP+/r3SKinAJmEPuam594COKKur4v5tmM2FIf4YMjGZw/+jO938zBvLehI9Xwq7X7Dm7J11mj9KekAFhwTF/qoiIhLbS6nru+Ocn7CmsoH+MXXYNj49h5e3pJPpzUZjifNj6rD3rV8EGe9ugyXDW3TD5crtftz4bg4qKcAmIw6XV3LTiQ0qr63nky6d3ve9bfbXdz/ujv8HxnRA3CM7/vl18x4/0a8wiIiLdUVPfyFf+tY7cYxU8fOfpnDtxkH/foPQgbHnaLr6bCu/kVMi4F6ZcBYMm+Pf9xK9UhEtA/PXNXI6V17L6q2cxc6Sr8yeoLIJPVsK6B+w5SofPgutWwPTrIMKh+VNFRETaYIzhW09sIHtfCX+6ZY7/CvCaUtj6HGx6yu5qgoHkOTD/lzD9WnCN9s/7SI9TES49rqHRwyubD5MxbWjnC3D3Afjwz7D+X9BQbc9Revbd9oBLLRsvIiJB6tODpby+7ShLLpvMlTOTu3cyj8de52LDY/Y6Fw3VkDQBLlwGpy2ApPH+CVoCSkW49LgPcosoqqzjqs5chEr2wru/h42PAwZm3gTnfBsGT+6pMEVERPzmxU8LiAoP47YzxnT9JOVHYcOj9qNkL8TEw+xb7ceINDVG9XIqwqXHvfhpAQOiI7hw8uD2D3bvh7fvt4vvsHC7r/c539YS8iIi0mt4PIaXNhVwweTBxMdGdu7FxtjdTNY9aLd6expg7Hlw8U/t1Z01f3efoSJcelRtQyOvbT3CJdOHERN5ilHZFYXwzm8h+yH7m/3cr8C534GB3byFJyIiEmDr9hZztKyWq2d14jOsvho2r4aPV8DRLRDjgjO+CulfUneTPsoPSxP2PgUFBViWxb333ut0KH3e2zuPUV7TwNWz27gQ1ZbDm7+GP8yG7H/ac3t/awNccb8K8D5IuSfiDOVeYL24qYDYyHDmTe3AXNyVRfDWcvjf6fDCN+1tV/8ZvrcDLv1vFeB9WEi2hCcnJ1NQUOB0GCHhhU8LSOwXxdnjk07c4Wm0+7i98St79a5p18K8n+li08cp93qHRo9hzbYjfJRX7HQoIS/SU0NsQylxDaXENpYR21DW/DOmsZzwqDjO+vLv2j2Pci9w7MkIjpAxbShxUacos9wH4IM/Qc4j9kDLiZfaEw+MPU99vUNESBbhEhiVtQ2s3V7IDWkjTlwNM/8deG2Zfbtt1JlwyyoYmeZcoCIC2P1YH/lwL/94P58DxdXERYV3fSVbOZkx9KeKJEpJpJREyrw/S0mgjATKcFGOizJcVOCijBjq2zxdHREcCNd0dMHmg9wiiivruGrm8NYPKM6zJx749An7+cyb4OxvwZApgQtSgoKKcOkxf387l+r6Rm5I9S6iU3oIXv+xvahA/GhY8JA9z7e+8YsEhWc3HOLeF7eRNiaBH10+lUumDyM8TPnZLo8Hqo5DWQGUH4Hyw/bPiiP2eJeKo96fhdBY2/o5ouOhXxLEJUHcSO/PBIhNhLhEn58J9iPGRVRkLON1/Qw6z28sYEBMBBe0nIygZJ899mnj4xAWAWl3wjnf0rzeIUxFuPSI/OOVrHg7j2tnJzNnxAD7ltub94FptOc1PefbEBnrdJgi4uPJdftJGdSPzK+ehaXizmaMvUCYe7+9OmHZoc9+lhXYjQsVR+wZLFqKGwT9h8KAoTBoEvQbDP2HQL8h0G+Q/e+4QXbBHeHH5cvFMVsOlfLcxkPcdsZooiO8kxFUFNqzfq1/2Gfigf+CgW20lEvIUBEufmeM4ecvbiUqIoyfpdbAAxfCkc12f7fLl0PiOKdDFJEW9hSWs25vCcsunxJ6BXhNmT0Hc9PDvc9utXTvh9IDUF914vERMTBwBMSPgLHn2oPIBwy3i6oByTBgmF1gh3dyajrp1RoaPfzwmU0kxEXxvfmTobbCboD64E/QUAOpt8P5S+z/b0RQES494PVtR/lo50FWT3qTxCf+bbcE3fgoTL1KXU9EgtSqdQeICLO4vqn7WF9TVwVFe7yPXPtncS4U59tdSXzFuOwuAoMmwoQMe52C+FEQP9L+GZeoa5mc5J/v57PlUBl/uWU28TtXwdpf2ndJpl0DF/8MBk1wOkQJMirCxa8K3NU88+xq1sb9hRH7D9t93ub/3F7lS0SCUl2Dh6dzDpExdSiDB0Q7HU73VLvh2A4o3AbHdsHxnXB8t92i7WvgCEhMgSmfs38mjLXv0rnGQKzLgcClN9tfVMXv1+xi0bgirvjoVji8EUakw02PwqjTnQ5PgpSKcPGbIncZ7/7lm/yt/nkaBo6E616AlAucDktE2pG1/SjFlXXcdHovWpm2sR6O7bRnWTq6BY5uhcLt9qDIJpH97Nbs0WfBoC/YLZFJE+2iOyrOudilTyksq+E7/3yd+8Ie4rrDb9ldk65/EE5boDsmckoqwsUvKvdvpOxfX+Cmxr0UTr6VITf8FqL7Ox2WiHTAk+sOkBwfw/kTB7d/sBPqa+xCu2ADHP7UfhzbAY119v7waBg8GcZdAEOm2o/BU+yuI2GaYlF6TmFZNY/85Vc8VPMQA8Lr4Kxv2/2+9fknHaAiXLrHGGo++DuRa35KfxPHpgtWMvPim5yOSkQ6aH9RFe/uPsa3Lp4YHNMRejxwfBccyoaD2XBovd21pGn2kdhEGD7LXs57+CwYOgOSJkC4Ps4ksIr2baPgkbv4fuMWyobOJWzhX+wvgyIdpKuWdF1VMfVPf5WY3P/wpmc2nqv/yrz06U5HJSKd8O+P9xFmWdxyukNzFddV2YX2/g9h/0d28V1Tau+LjocRc+wpTZPnwPDZ9uBI3eIXJ3kaKX7jj/R/79eMM5HknX0fKfO/qrsu0mkqwqVrDq6ncdUdUH6EXzXewRk3/5j504c5HZWIdEJ1XSOr1h3gsunDGBYfE5g3rau0C+6978He96Egx9vKbdndSKZfByNPh5Fz7RZuFTYSTIrzqXrqLhKPrONt0ki65W/MmKLWb+kaFeHSOcbAugfxvLaMox4X32z4OXfffiMXTR7idGQi0kkvflpAaXU9t581pufepLHBbunOfQPy34aD6+yiOywCklPh7G/C6LNh1Fx7JUiRYGQMZsOjNLy8lIYGuDfim9x21xImDhvodGTSi6kIl46rr6H62W8Ru20VbzbO4W8JP+C+285n0tABTkcmIp1kjOHhD/YyeegAzhiX6N+TlxXA7jWw+3XIfwdqy8AKs7uTnHU3jDsfRp8JUf38+74iPaG6hKrMrxOX+wrrGqfxWPIylt2cwcgEzbAj3aMiXDokd89OIjLvYEzNDv6v4QYqzvguj10+9bNleUWkV8nZX8K2w2X893Uzur9CpjH2vMg7X7UfRzbZ2weOhOnXwvh5duEd5+diX6QHGWP49KM1jFx7N/H1x/mt+TzDrvg+fzpjLGHBMIhZej0V4dKutWtf47R3FtOPGp6csJwbPncHoxLVAiDSm63OPkj/6Aiund3FJbQ9jXbf7u0v2o+yQ3Zr98jTIeNemHip3cdbgyilFyqvruO1B3/GtcdXctRK4rXTHuT2+VcEbuyEhAQV4dKm+kYPTz/2d67JvYeKiATqb3+Bm8fOcjosEfGDjQfcpI9NoF90Jz4GPB678N76DGx7ASoLISLGbum+6Mcw6VLoN6jnghYJgN0Hj3Dg4S+zsOE98gdfRPIX/8nn++sujvifinBpVUNDI6v/uISbS/9BQf/pDFv0NBHxmv1EpC+oqW9kd2EF86cN7dgLjm6FTatg89NQdhAiYu2Ce/q1MPES9e2WPuOj7E9IevFOLrAK2D/nB4y7+se6myM9plcW4ZZluYBlwCpjTI7P9gygGMgAsnz3SSd4Gtn/2N3cWvY4e4dfwtgvPwKRsU5HJSJ+sv1wGY0ew/Tk+LYPqiqGTU/Bp4/bK1SGRdgt3hn3wuTLtSKg9DkNe95k2ku3QVgYpdc9yeiZlzodkvRxvbIIB9KBFN8NlmWlAEuNMfMty0rELtIXOhFcr1Zfg3l2ESn5z7M66lpuuOufEK7BlyJ9yZaCMgBOG9miCPd4YO87kPOI3c+7sc5elfLy+2HGDepqIn3Xun8Q9soPOOwZTvHVj3DWzDSnI5IQ0CuLcGNMlmVZC1tsywPme5+mAOsCHlhvV1sBT96Klf82v6y/jWlX/YgwFeAifc6Wg6UkxEWS3DTIrNoNGx+H7H9A0R6IcUH6l2DO7TBshpOhivQsjwey7oEP/sjH4en8Ln4JmWmpTkclIcKxItzbpWQRkGSMWdrK/iVAHpAIYIxZ2cHzLgDGt3ZOOYXqEnjsRsyhbP404Lu8Vn8eP5yd7HRUItIDthSUMmNEPFbRHvj473YBXl9lz2xy3QqYdo26oEnfV18Dz30Vtj7LvpRbuW3bFfzh6tO6P2WnSAc5UoR7+267gPFt7F8OrDPGZDY9tyxrQdPzUzHGZFqWlWJZVoYxJsufcfdZlUXw6DVQuIN9F/+V3788kJ9fnUJkuJaLFulrausbGFi4jnsGrYU/vwfhUXDajXDGIrvriUgoqCmDJ26Bfe/B/F/y3U/TGJ1UxxWnDXc6MgkhjhThTcWxZVlzsYvxlha1aMleBSwH2izCvS3rGGPcQBawmjaKfPFReRz+dTUU58KtT5KZO4rwsFyuUSu4SN9iDOx8lcY3fssTETnUViXAhcvsbif9hzgdnUjgVByDx26wZ/25/gEKx17N+hfX8oNLJxOuRXgkgIKuT7hlWa11xnJjz3hyKouAJGCp93hN6tmeimPwyNVQnAe3PAnjLyLrpXdIG5OAKy7K6ehExB88HtjxIrzzWziyGU/sCH5a/0UWffWnjBqmgZYSYkoPwiPX2D9vfgImXcLaT/YDkDG1g1N2ivhJ0BXh2MVzcYttJzz3dmdJB9yWZeGdinAlkOHdtxDNjHJqVcX2hag4H25dBSkXcrCkih1HyvnRFVOcjk5EussY2PUfeONXcHQzJI6Ha//Gb/Km8sKmQn4xNMnpCEUCy70fHr7SHgN1+3Mw5iwA1mw7yqjEWCYN1bSbEljBWIS72tphWZbLGOP2dmc5Yf4gbzeUpu4q6gt+KjWl8O/r7VkQvAU4wBs7CgGYp9YAkd5t3wew5mdwcB0kjIPrVsJpCyAsnM3vvWcPytTgMwklxfl218vaUrsAH2mXEFV1Dby35zi3nTFaOSEBF4wj79yc3JXEr11LCgoKsCyr+XHvvff68/TBra4SHrsRjmyGGx+B8Rc178raXsi4Qf0YP1itAdIzQjr3AuHYLnuw2UOXQ+khuOoPcPc6mHUThIVT3+hh+5FyZow4xSI90ieFdO6V7LNbwGvL4I4XmgtwgHd3H6euwcN8NT6JA4KxJbyYk1vDXdDc2t1tycnJFBQU+ONUvUtjPTz1BTj4CSz4J0y+rHlXRW0DH+UWccdZYxwMUPq6kM29nlZdAm8th09WQmQczPsZnPE1iIo74bDdRyuoa/CoCA9BIZt7ZQX22Ke6cvjCiyfNAJS17SgDYyKYO07DyCTwgq4IN8bkWJblbrE5EXUx6R6PB56/G/assVvHpl93wu73dh+jrtGjrigivYinsZH8rBUkZy8npr6UrcOv5+OxX6OmIQHeP7ng2nGkHIAZyQMDHapI4FUU2l1QKovgjudPKsAbPYY3dhRy0ZQhmpJXHBF0RbjXUy3mBZ8PrHAyoF4v62ew6Um46CeQ9sWTd28vZGBMBOljEwIfm4h0Sl2Dh8eef4m0Lb9iptnJx54p/Lx+Cdvyx0J+IVDY5mvHJsUxNqlfwGIVcURNKTx6PZQdgs8/c0IXlCYbD5RQVFmnWVHEMU4t1pOKPeXgAu/zJUCWd5YTjDGLLcta4p3pJAXI7chCPR3V1DfunnvuCY1+cR+vhA/+BHO/Aud//6TdxytqeX3rEbUGSI8LudzrCfXVbHh4KbcffJTq8AHkzPk1U+ffxXORHbucR4RZhGku5JATUrlXXwNP3ArHtsOtTzXPgtLSox/uIyo8jAsmDw5wgCK2Nq/almVdD5wOmG6+xzpjzDO+G7zFdg5wf1svMsa0ua+7Qqpv3K7/wGtLYdLlcPn90Mro71++tI3q+kbuvmiCAwFKKAmp3OsJ+z6g5plvcEZpHh8nXMEZi/5Capz6skr7Qib3PI3wzFfslTCvfxAmzGv1sDd3FPLcxgK+NW8iA2MiAxykiO1UTSfzgR/64T1WAs+0e5T43+FNsPpOGHYa3PAghIWfdMibOwt5fmMB3543kYlDBzgQpIi0q74G3vgl5sO/UGIN5tcR9/Crxd+EWBUPIif4z49g+4tw6X0ws/XlQspr6vnxs5uZOKQ/37hIC2uLc05VhOcZY0q7+waWZWV39xzSBRWF8MTNEOuCW1ZB9MnTDlbWNvCTZ7cwYUh/vq4LkUhwOrwJnlkEx7azYch1fH7/Vfz1zvOIVwEucqKPV8LHf4czvw5nfb3Nw+5/bSeHy2p45mtnEx1xcuOUSKC02QHYGPNbf7yBv87jT0194/psv7iGWlh1u70q5i1PwsDhJx1S29DIksxNHHJX85vrT9OFSAKiz+eePxkDH/0NHpwH1SXsnv8wCw4s5Oq5E7lw8hCno5Neps/nXlPXy8lXwCW/avOwRz7cy6Mf7eNL54xjzmhNRCDO8svATMuyxgGpQAmQAJiW/cCDSZ/uG2cMvPJ9OPARLHgIhs886ZDSqnoWPZrNx/nF/OiKKaSPVZ9SCYw+nXv+VFkEz38Ddr0Kky6j6oo/ctcD20h2GX5y5TSno5NeqE/nXuF2yPzSKbteNjR6+OVL2/jXh/vImDqUH1w62YFARU7U7SLcsqzZAMaYp1tsvwFYY4wp6+57SCdk/xNyHoHzvgczrj9pd2F5Dbc+8DH7i6r4w82zuWb2CAeCFJE2HcyGp+6AymNw2XI4YzH3v7iNvUVVPHHXmfSPDtaZZUUcUF0CT94KUf3sO79RJ0+/6fEYvv5YDq9vO8pd543jh5dPJVwzBEkQ8MfVfHzLAhzsotw7w0rQtoj3OQfXw2s/hAnz7fnAWzDG8P3VmzhYUsUjXz6dM1OSHAhSRFplDKx7EF5bBgOT4ctrIHk2H+Qe5+EP9nLnOWM5a7xyVqSZp9FuAXcfgC++bOdNK/72di6vbzvKT6+cxpfPHRfgIEXa1q0i3LKssdhTDTY9/z6QY4x5w7up2wM7pYMqi2D1F6D/MLh+JYSd3N3/Xx/s5Z1dx/jltTNUgIsEk4Y6eOV79l2siZfYORxr91f909o9jHDFsuTSKQ4HKRJk1v4Cct+Aq/4Io89o9ZBP8ov5n9d3cvWsZL50ztjAxifSju6uzJIAFPk8vxl7asMm3Z1jXDrC44Fn7oKKo3Djv6CVeYN3HS3nvld3cPGUIXz+jNEOBCkirao4Bo9c7e1G9n17NiNvAX6wpIoP84q4ae4oYqM0eFqk2Y5X4P3/g7Q7Ie0LrR5SXFnHt57YwOjEOP77uhlYrayTIeKkbhXhxpgNQLrP83RjzDKfQ4Jy6HGfGyX+wR8gdy1c9hsYkXrS7roGD995ciP9oyNYfsNMXYjEMX0u97rr2E548GIo2Ag3/APm/fSEu1jPbTgEwHVzNHZDuqdP5V7JXnjuqzB8lv251wpjDD96ZjPFlXX8+dZUBmhBHglC/ugT3mqh7Z0xJc8P5/e7PjVK/GA2vPErmHYtpH+p1UP+9lYu2w6XsfL2NAYPiA5sfCI++lTuddfe9+wBZeHRcOfLMCLthN3GGJ7JOcTp4xIZlRjnUJDSV/SZ3GuohdVftO+zL/wXRMa0etiLmw7z2tYj/PDyKcwYER/QEEU66pQt4ZZlDWzvBN4BmDdYlnWx9zXx3plR4r0t5dJTqt2QeScMSIar/tDqkvQ7j5Tz5zd3c/WsZC6ZPizwMYrIybY8DY9ca4/h+ErWSQU4wMYDbvKOV7IgdWTg4xMJVmt/AQUb4Nq/QmLrgyyPlddyz/NbmD3KxV3npQQ4QJGOa68lfJn3cUpNs6NYljUPKGptthTpAS9/D0oPwZdes1fGbKGh0cOSzE8ZEBPJPVdpbmGRoLDuQXj5+zD6LLjl8eb+3y09nXOQ6IgwLj9NX55FAMh7Cz78M6R/GaZe2eohxhh++twWKusa+d3CmZqKUIJae33CF1iWdVFHT2aMWWuM2di9kKRDtjwDWzLhgqUw6vRWD3nwvXw+PVjKz6+eTlJ/dUMRcZQx8PZv7S/Pky6F259pswCvbWjkxU8Pc+n0YerLKgL2CtDPfg2SJp5yRczV2Qd5besR/itjEhOGDAhggCKd115L+CAgyzuQLwfIAtYA2a0twmNZ1kAtzhMA5UftD/LkVHtRnlas31fM7/6zk8tnDOPKmScvWy8iAWQMvPFLePd/YOZNcM1fILzt4vqN7YWUVtdzQ5q6oohgDLz0X1BZCLdkQVTrYyQ2HXTzk+e3cM6EJO46T/OBS/BrryU8FfgqcCPwFJCGXYiXWJa127Ksv1mWdb13vnDoQNeVYNCrR4kbAy9+C+qr4LoVEH7y96jiyjrufnwDya5YfqPZUCSI9Orc6ypjIOseuwBP/QJc+/dTFuAAT2UfYMiAaM6dMChAQUpf16tz79MnYdtzcNGPIXlOq4cUVdTy1UfXM7h/NH+6JZWI8O7OwCzS807ZEm6MyQce8Pb1zjPGXAJgWVYqdkE+H3gQiLcsy+19WdAX4r16lPjGx2DXa3Dpr2HwpJN2ezyG76zaSFFlHc987WziY3UrW4JHr869rjAGXv/JZ/1Yr/hdqwtp+Xpl82He3HmM782fpP6s4je9NvdK9sErP4Ax58A53271EGPsz73jlXU8/dWzSewXFeAgRbqmQ1MUGmPWAnhnPck1xuRgd095wLt9HHZB/oMeilMA3Pvh1R/CmHPhjK+1eshjn+znnV3H+PV1p2laJhGnvfnfdgE+9y644retzmDk60hpDcue2cyskfF89cLxAQpSJEgZAy/cbf/7ur9DWOsLVr265Qjv7j7OL66Zzmkj9bknvUen7td4Zz0p9XZBGeizPd8YsxLQrCg9xRh44ZuAgWv/0mprWmVtA3/I2s3p4xK55fRRgY9RRD7z7u/hnd/CnNvh8vvbLcCLKmr59pMbqGvw8L83zSZSt9Ml1K1/CPLfgUt+Ca7WV3qua/Dwm1d3MHnoAG47Y0yAAxTpnk4v1uPtopJvWdY8y7LijTHP+Oxe4b/Q5AQbH7OnZ/rc/0DC2FYP+ed7+RyvqGXlHWnqBy7ipHUPwtqfw4wF9hz+p+iCUt/o4V8f7OUPa3dTXdfI8htmkjK4fwCDFQlCZQWw5h4YdwGkfbHNwx75cC/7i6t4+M656r4lvU6XV8xs0UWlxBjzhrdAF3+rKIT//BhGnw1pra+KWVRRy4p38rh0+lBSR7c+7ZmIBMC2F+x5wCdddspb6ADVdY184aFP+CS/mAsmDeanV05jwhAV4CK8ugQa6+Cq/2vzLpK7qo4/vbGH8yYO4sLJQwIbn4gfdOt+p7dLSi6QYlnWKp9ZUsSfXv8J1FWeskXtL2/mUlXXwA8unRzg4ESk2d734emvwMi5sOChU86CUt/o4e7Hc1i3t5j/WTiLh++cqwJcBGD7S7D9RXsdjMTWV7ysbWjkx89uoaymnh9dMTXAAYr4R7vL1luWNdvbB/z73ikJ/+OdnrARKMEeoHk/MB5YHICYu61XTdWU9zZsWgXnfqfV2VCMMTzwTh7/+nAvN6aP0uIEEtR6Ve511rGd8OQtdt/VW1e1OZcx2LMYLX16E2t3FPLLa2ZwQ9pIdSGTHtVrcq+mzJ4NZch0OPubrR5SUlnH7f/4hJc3H2bJpVOYOnxgq8eJBLv2uqO4AeP9mQ/kARuATO+/83pjF5ReM1VTQy28/F1IGNfqojzlNfX88OnNvLz5MJdNH8ZPrtTS9BLcek3udVZlETx+I4RHweefhrjENg+tqmtgSeYmXtp0mO/On8Tnz9RgMul5vSb33voNlB+Gmx494U5So8ew/XAZH+cX8+iHeykoreEPN8/mmtkjHAxWpHvaK8LzsAdbuoGs3lhw92rv/wGK9tgf6pGxuKvqeHvXMXL2lbDhgJttBWV4jGHZ5VNYdH6KWtJEnNBQC6s+D2WH4YsvQ0LbRfX+oioWPZrNzqPlLLlsMl+7QNMQijQr3A4f/x1S74CR6QBsPljKn97YzUd5RZTVNACQMrgfT9x1Bmlj2v6yK9IbtFeEZxpjfgtgWdYcy7K+4rMvzxjzhu/BlmVd32K2FOkq9wF7irNp12LGz+O5DQf5xYvbKKmqJy4qnJkj41l0fgqXTB/G7FEup6MVCU1Ny2nv/wBu+AeMmtvmoYVlNVzzl/fwGHj4ztO5YNLgAAYqEuSMsbuhRA+AeffQ0Ojht6/vZOU7eSTGRXHFacM5MyWJM1ISGR4f63S0In7R3oqZP/T59wbsriiAvUCPZVl3AS7sLiulwBJARbg/rPkpADUX/4KfPb2Jp7IPkjraxT+unMbMEfFaklckGKx70J4+9PwlcNqCUx76f2t3U17TwKvfPo+JQzV2Q+QEW5+Bve/C5/6HY57+3P3gx3ycX8wtp49m2RVTGBij1Z+l7+nOFIX5eFfMhOal7JP8EVTIy38Xtj6L+4zv8/kn9rPlUBnfungC387QMtYiQWP/R/DaD+2pCC9cdspDc49VsGrdAT5/xmgV4CIt1VXBf36CZ9hMnjYZ/O5P71JaXc/vb5zF9akjnY5OpMd0uQhvyRiTY1nWSn+dL2Q1NlD70g+ojBzGee+eRlhUFQ/ekU7GtKFORyYiTcoOw1N32DOhXLfilIvxAPzuPzuJiQjjm/MmBihAkV7kwz9DeQFfq/4a/3lmK9OTB/LQF09nWrJmPZG+rc0i3LKsscaYvZ05mW/3le6cJ5Stf+4PpBVtZ0njf3HjWRNZdH6K+r+JBBNPoz0XeG053P4cxLpOeXjO/hJe3XKE/8qYxKD+0QEJUaTXKD9C/Tu/J6txLiXD5vKvGydw/sRBmmhAQsKpWsKXAl/zw3ssB27yw3n6vE25Bxmz6Q9sj57BPd9aRqI+sEWCz7u/h33vwTV/haHtTwv629d2Mqh/NF85b1wAghPpXY4+/1MSGupYO/IbPP6VMzTeSULKqYpwy7Ks+7p5fsv7CCpNixbcc889QbNwwfGKWrIfv5cvWaVE3vgU8SrApQ8KxtzrlP0fwVv3wYwFMPvWdg9fv6+ED/OK+MnnptIv2m+9/0Q6LRhzr2DnOobtWc2z0Vfxsy9cqQJcQk57LeF9UrAtWmCM4WePZvG7hhdwj78K14QznQ5JpEcEW+51SnWJ3Q3FNQqu/F/owO3yv721B1dcJLecPjoAAYq0LRhzr+i5H9GPOE7/wn2a/URCUptFuDGmNJCBhLIth8o4/9BKoiMbibvyl06HIyItGQMvfNNeye9Lr0NM+wPGdh4pJ2t7Id/JmKhWcJEWire9xWnVn/DmqLu5aIRmQJHQpHs/QeC9D99nYfjbNKR+GRLVb1Qk6OT8C7a/CPN+BiPTOvSSv7+dS1xUOF84a2zPxibS2xhDzWv3UGhcpFzxHaejEXGMinCHNXoMk7b9kfqwGKIv7rM9gER6r+J8eO1HMO58OOubHXrJgeIqXvi0gFtOH01Cv6geDlCkd/HsWkNy2UZecn2eMcO1cqyELhXhDtuy/h3mmQ85MPlO6Ke1jkSCiscDz38DwsLt2VDamQ+8yX2vbifcsjQjikhLHg+Vr93Dfs9gBl94l9PRiDhKRbjDot65j1LTj1Gf+4HToYhIS9n/gH3vw6W/tgdkdsArmw/zyuYjfDtjoub4F2lpx0sMKNnGA+E3cclpHcspkb5KRbiD6vZ+yNTyD3lv6K3EDEh0OhwR8eU+AFn3QsqFMOfzHXpJcWUdP31uC6eNiGfx+Sk9Gp5Ir2MM9W/+hjwznJjUm4mOCHc6IhFHaci+g8pf+QUeM5CBF3asn6mIBEZVbT01T32D+MZG1s+4h5rdxzv0usc+3kdZTT2PLdSiIyItmZ2vEHlsK39v/DpfO0tfUkVUhDukfu9HJBV+wP+G38E3J+uWnEiw+DiviP88+Ud+Vvc299bfwcNPHQYOd/j1379kElOGtT+FoUhIMYaSV39FuWcIkzLuZNygfk5HJOI4vxXhlmXNAxZ4n64xxjzjr3P3RXlP38MgM4DJn/u2WsxEgoAxht+8uoMn3t3MW9EP4U6cxVVX/Yyrwjp+yzwmMpxpw1WAi7R0NOdFhpZuY3Xid7nr/IlOhyMSFPxShFuWdQOQAmR6N6ValvUVY8yD/ji/vzm9fO+aNa8wv/wj3hz9da5ImxDw9xdxitO5dyqfHixlxTt5PJb8Cgkl5VgL/0Sapk+TPsLJ3CupqKXk1f+mgcFcdft/ERbW/mqzIqHAXy3hbmPMb32er/W2jAclJ5fvfXNHIda7v6MiYgDn37rMkRhEnBKMS2c3+SS/iBlWHmcXP4d1+iIYPsvpkET8xqnc219Uxe8ffIj/a9jBjtR7mJKkO0UiTfxVhJtWtsX76dx9xsPv5/PUy6/yStR6qs/6IeGxuhiJBIt1eUX8NvZfWLGD4eIfOx2OSK+3taCUL/zzE37fkEl9TBJTLv+a0yGJBBV/FeEJlmXdB6zzPp8L5Prp3H3CL17cxj/fz2fVoCxMXT9iz9XFSCRYeDyG5H3PMtXshktWQozaEES6o6a+kbsf38AkDnA+OXDWjyFS8+aL+PJLEW6MedqyrDzgJu+mVcaYDf44d1+wZttR/vl+Pt9KjeL07W9hnfk1iHU5HZaIeOUdOsI3PI9zPGEWg2be6HQ4Ir3e/67ZRf7xSp6c9i4ciIO5X3E6JJGg0+Ui3LKssS02lQB/99l/nzEm5Ds9V9U1cO8LW5k0tD/f7vcSlhUGZ37d6bBExEfNW//DBMvN4YxHwdKgMZHu2HjAzQPv5rF4djRDd71gF+BxWpBOpKXutIRnAesB30+spr7hFjAHCPki/A9Zuznkrua5L04mPPNRmHkjxI9wOiwRaeLez+S8h3nFOp/Lp5/ndDQivVpdg4clmZ8ydGAM341/A4xRw5NIG7pThC82xqxta6dlWXO6ce4+YfvhMh58L5+b545i9pGnoaEaztbqmCJBJevnNBqL98d+nSvUCi7SLS9+WsCuoxX845apRL/6KEy9ChLGOB2WSFDq8iox7RTgs4FxXT13X/GbV3cwMCaCpfPGwCcrYOKlMGSq02GJSJOCjbAlkwcarmDCxClORyPSqxljeOiDfCYO6c/FdW9BTSmcqUkIRNril6UaLcsaZ1nW65ZlrbMsKxt4EHuGlJC15VApb+86xlfOSyEh73moKoKz73Y6LBHxtfYX1EXGs7LhSuaOVZ9Vke5Yv6+ELYfK+OLZY7A+WWHPtT/qDKfDEgla/pqicLH3AVDs/Znip3P3Sn97O5f+0RF8/ozR8PAdMHQGjFV/U5Ggkf8O5K4la/g3MHUDmarl5kW65aEP9jIwJoIbEvbAsR1w7d800FnkFPzSEo49JWG+MSYfmGeMKQUS/HTuXifvWAWvbD7M7WeNIf7oR1C4Fc74qi5GIsHCGLsveP9kfnX0bOaOTSBcS2mLdFmBu5rXthzh5tNHE7P+QYgbBNOvdzoskaDmryI80bKspoV6Tvcu3LPQT+fudVa8nUdUeBhfOmccfPR3iEuC00L2zyESfHa8DIeyWRl+IyX1EfzoCo3VEOmOf3+0D2MMX5wK7HoN0r8EkTFOhyUS1Py1WM9avH3AjTE/tCzrBiDHH+fubQ6WVPHMhoPcPHc0g+sLYOcrcN73dDESCRbGwFu/oThmFL87msbvbprBxKEDnI5KpNc6WFLFvz7Yy6XTh5G850mwwiDti06HJRL0/NUSfgJjzNN8Nmd4yDDGsOyZzUSGh/HVC8fDugchLFwrhYkEkZINz8HRzfyq7HPcePo4rpsz0umQRHotYww/fnYLBvjRpeNh42Mw6TKthyHSAX5pCfeunrkQaJpewALmEWIzpDy57gDv7j7OL6+dwYh+wIZ/w5QrYeBwp0MTCXkNjR5++PQmvrjlXkqtodRPW8Cvr5rmdFgivdozOYd4e9cx7r1qGqOOvgGVxyD9TqfDEukV/DU7yg+BNZzYBcXlp3P3Cofc1fz3y9s5KyWJ204fDZuegBo3zP2y06GJCPDWzmOUbHyBGVF7KZr3v/zpvJBqIxDxu2PltfzipW2kjUngjrPGwiPfhPjRMP5ip0MT6RX81R1lhTHmaWPM2qYHsNRP5/a7goICLMvi3nvv9ds5f/bcFjzGcP+CmYSFWZD9T0iaqGkJRXz0RO511Ors/Xw36jmMawxJZ98e8PcXcVJP5N7Kd3KprG1g+Q0zCSvOhb3vQtoddjdMEWlXl1vCLcvynVQ3wbKsi4E8n22LgWVdPX9PSk5OpqCgwG/nK3BXs3ZHId+eN5FRiXFweBMcXAeX3qdpCUV8+Dv3OqqoopaanWuZHpkL5/0BwiMDHoOIk/yde3UNHp7JOUTG1KFMGNIf/vMQhEXAHH3BFemoUxbhlmUNNMaUtbE7B1iP3f+7NXMI0iLc317aZF/YrpvjHYiS/U+IiIHZtzgYlYg0eX5jAV8Je5GGuKFEzFJeinTXGzsKKaqs48a5I6GhDj59AiZfDgOGOR2aSK/RXkv4MtoupBd7u520yrKsOV2Oqpd54dMCZo2MZ+ygflBbDptXw4wbIDZk1ysSCSobP36LL4VvgbN/DhHRTocj0uutzj7A0IHRnD9xMOx6BaqK1Aou0knt9QlfYFnWRa3tOFUB7t2/octR9SK5xyrYcqiMq2Yl2xu2PA11FZCm0eEiwWBrQSnz3U9SF9FfszaI+MHRshre3FnIDakjiQgPg42PQ/+hMH6e06GJ9CrttYQPArIsu19zDpCFPQtKdmvdVNrpvtInvbCxAMvisyJ8w79h8FQYme5sYCIhxhjDP97LZ29R5Qnbiw7s5M9hH1OfejfExDsUnUjf8XTOQTwGFqaPgopjsPs/cObXIdxfE66JhIb2MiYVyACKgRRgPvasJ8ayrDw+K8pzjDF7OXX3lT7HGMOLnxZw5rgkhg6MgcId9oDMS/5bAzJFAuzfH+3jVy9vxxUXSbhP/i3xrIawcGLO/YaD0Yn0DcYYVmcf5PRxiYwb1A8+fBg8DTD7NqdDE+l1TlmEG2PygQcsy5oH5BljLgGwLCsVSMMuyh8E4i3LcntfFjJF+NaCMvKOV3LX+Sn2ho3/tkeHz7zJ2cBEQkz+8Ur++5XtnD9pMP+6cy5WUxFe7Ybf3wHTb9KiWSJ+sONIOfnHK1l8fgoYAxsegxFpMGSK06GJ9DodunfU1P/bsqwbgFxjTA5295QHvNvHYRfkP+ihOIPSK5sPExFmcfmMYdBYD58+aS/X23+w06GJhIxGj+F7T20kKjyM+2+Y+VkBDnb3sPoqOGOxcwGK9CEf5BYBcN6kwXD4UyjcCp/7H4ejEumdOrVYjzHmaaDUsqzrfecJN8bkG2NWAk/7O8Bg9t6e48wZ7cIVFwW7X7eX69XocJGAeuj9fHL2u/nltTMYFh/z2Q5PI3yyEkafDcNnORegSB/yYe5xxibFMcIVC5tWQXiUPRuYiHRap1fM9BbczwBzLcu6vsXuFf4JK/iVVNax+VAp507wtnpv+Lc9OnxChrOBiYSYV7ccYdYoF1c3DY5usus1cO9TK7iInzQ0evg4r5izxg+yv+RueRomXqLpeEW6qMvL1nuXp3/GsqwbvKtlNvUhDwkf5hVhDJw7MQkqj9st4TNv0uhwkQDyeAzbD5cxZ5TrxG4oAB//HQaOgClXOhOcSB+zpaCM8toGzh6fZC9RX3EUTlvodFgivVaXi3BoXro+F0ixLGuVZVlj/RJVL/Du7uMMiI5g1kgXbH3WHh0+62anwxIJKfuKq6iqa2Ta8IEn7ji6DfLfgblf0RdjET/5IPc4AGemJNmL0kUNgEmXOhyVSO/V7rL12FMTNj3Gt3gO9rL1biAPWEyIzI7y3p5jnDk+yV6oYNMqGDoDhk53OiyRkLK1oBSAacktivDsf0J4NKR+wYGoRPqmD3OLmDx0AINjDGx7EaZeBZGxTocl0mu110TkBoz3Zz52ob0ByPT+Oy+UuqA02V9UxYHiar5ybgoU5dpzg8//hdNhiYScbQVlRIRZTBza/7ONdZX2F+Np10C/JOeCE+lDahsaWbe3mJvnjoY9a6C2FE5b4HRYIr1ae0V4HvZgSzeQFYoFd2ve3XMMgHMnDoLNfwYsmKGLkUigbS0oY8KQ/kRHhPtsfBZqy7REvYgfbdzvpqbeY/cH37wc+g2GcRc4HZZIr9ZeEZ5pjPktgGVZcyzL+orPvjxjzBu+B1uWdb135pQ+7b3dx0mOjyElKc5ucRt3HsSPcDoskZCz7XAZ509sMS9/9kMwaDKMPsuZoET6oA9yiwiz4IwRUfDMa5D2BY23EOmm9lbM/KHPvzdgd0UB7AV6LMu6C3Bhd1kpBZYAfboIb/QYPsgt4tLpQ7EKcqA4D877ntNhiYScwvIajpXXntgf/MhmOJQNl94HLWdLEZEu+zCviOnJ8cTvXwuNtTC95QzFItJZXf4a27SkfdNz71L2fb4D5u7Cckqr6zlrfBJs+qM9+GvqVU6HJRJythWUATDdtwhf/7Cdk5qpSMRvahsa2XjAze1njoHtf4f+w2DUGU6HJdLrdWuKQl/epexX+ut8p2JZlsuyrOXewt93+yLvY4VlWSltvb47thyyP/hPGz4Atj0PE+dDTHxPvJWInMK2w3YuTm2anrCuCjY9BdOvhbhE5wIT6WM2HyylrsHDmSOjYfcau+EpzG/lg0jI8msW+XZf6WHpfDZFItDcEp9tjFkJrKaHVu/ccqiUuKhwxlVvhoojMP26nngbEWnH1oIyRiXGEh8baW/Y8bI9IHPO7c4GJtLHfLK3GIAzG3OgocaeeUhEuq1XfpU1xmQBxS02p2DPUw6QjV2o+93WglKmDR9I+LbnICIGJl3WE28jIu3YXlB24iI9Gx+D+NEw5hznghLpg9blFzN+cD8G5L0McYNgzNlOhyTSJzhWhHu7lCyxLGt5G/uXWJa1oKmLSXvnM8ZkGmOaivB07ELcrzwew7aCMmYM7+ftinIJRPdv/4Ui4leVtQ3kF1Uybbi3K1jpIch7y+4LrtvkIn7T6DFk7yvh7DH9YNfrMPVKCAtv/4Ui0i5HPq0sy8oAMrBX4HS1sn859hSImd7uJeMty+rMRNyLgYX+iNXX3qJKKusauSBmD1QWwgyNDhdxwrbDZRjjMyhz0yrAwOxbHI1LpK/ZeaSc8poGrojZBvWV6ooi4keOFOHGmCxjTCb2IkCtWeTd32QVn3U1OSVvq/lSY0xb5+6yLd7ZGGaVvQmRcXZLuIgE3EufFhAVHkb62AQwBjY+bs8Lntgj47FFQtY6b3/wmeVvQWwCjD3P0XhE+pKgu2/bcsYTLzd2y3l7r83AXtkzz/tvv9p6qJTYcEPCvtdg0qUQ1c/fbyEi7aipb+TZDYe4bMYwXHFRcGg9FO2GWWoFF/G3T/YWM2JABHH71sKkyyE80umQRPqMYFzuKpGTB12e8NxbYKcDbsuyMMbkeIv31UCxZS/SkQNk+TOwrQVlXJu0H6vsmG7JiTjktS1HKKtp4Oa5o+wNGx+3B0lPv9bRuET6GmMM2XuLuW3YIawDpTD5cqdDEulTgq4lnFb6iDexLMsFzd1Z0owxS73zk2OMyTHGJBhjxnsfbfYJLygowLKs5se9997bblDGGLYUlHJl5Hp7MZAJ8zv7e4mEvK7kXktPrtvP6MQ4zkxJgsZ62PacXRxovn6RNnUl9w4UV3O0rJaMsPUQHgXjL+75QEVCSDC2hLuxW8N9+XXljeTkZAoKCjr1mkPuatxVdcyKeh/GX6RZUUS6oCu55yv/eCUf5RXzg0snExZmwe63oaoIZnRm3LZI6OlK7uXsLwEME0rehXEX6HNPxM+CsSW8mJNbw10APTHYsqO2FpQx3dpL/5rDMOVKp8IQCWlPZR8gzIIFaSPtDVuehuh4e+VaEfGrHUfKmRpeQFTZPnVFEekBQVeEe7uXuFtsTsTP/bs7a+uhUi6LWI+xwnQxEnHICxsLuHDyEIYOjIH6Gtjxkj1vcUS006GJ9Dm7j5azoP9m+4kWphPxu6Arwr2eajEv+Hz8uAx9U9+4zvRH3VJQxpWR67FGnwX9BvkrFJGQ0pXca1JV18AhdzVpYxLsDXvW2MvUz7jBv0GK9EFdyb1dheVcZK2H4bMhfkSPxSYSqhzpE+6dySQDWOB9vgR7asGmQZaLvStmZmAvR5/bYt7wbulK37iyQzsZ59kHU77qrzBEQk53+oQfLKkGYGRCrL1hy9P2EtrjLvBXeCJ9Vmdzr6qugeriI4yL2QaTl/VgZCKhy5Ei3Fts5wD3n+KYNvcFWkllHXOq3odIYMrnnA5HJCQdKK4CYFRiHNRWwM7XYM5tEB6M48tFerc9hRVcFL4BCwOT1RVFpCcEa3eUoLL9cBnzw9dTkTAVEsY4HY5ISGouwhPiYNdr0FCtrigiPWTX0QouCNtEQ9xQGDbT6XBE+iQV4R2Qu38/adYuwjQgU8QxB0qqiY0MZ1D/KNj+IvQfCqPOdDoskT5pz5ESzgvbTNjEDLAXwBMRPwvJIryzA1Ss3DcItwxx06/o2cBE+rjuDMzcX1zFyIRYrIYa2L3G7hoWFpKXMJFO62zuNe7PJt6qtItwEekRIdmZsrMDVEYce4eysHgGjkjtwahE+r7uDMw8UFxl9wfPfRPqK2HqVX6OTqTv6vTn3vH38RBG2PiLejAqkdCmZqR21NfXM7t2PfsSzoawcKfDEQlJxhgOllQzKiHW7ooSEw9jz3M6LJE+qbK2gTn16zkycAbEJjgdjkifpSK8HYe3vEuCVUHNON2SE+fk5OSQlpbG0qVLnQ7FEe6qeipqGxjtioSdr8DkKyA80umwJASEYu7l79vLrLA8KkapFVyc1dfzT0V4O6q3vUKDCcN1mqZoklNbuHAhK1eu7JFzp6amsnjx4k69prVYejLGnnSgxJ4ZZWbjVqhxqyuKnEC5518V214HIG7apQ5HIr2B8q/rQrII78wAlYSDb7DeTGbsyOSeD0x6tcWLF5OR0XN3TBITEzt8bF5eHm63+6TtPR1je7o6MPNAsb1Qz/jjb0BkHIy/uAeik95Kude+zuRev/1vUWQGMnyKZh+S9in/uk4DM0/FfYAh1bm82O9OzggPye8r0gnBlODLly9n/PjxJ213OsauDsw8UFKFhQfX/tdhQgZExvZAdNJbOf3/ta9en3seD2PcH5EdncbF4RoHJe1z+v9tX8Gaf21RZXkqe9YAUJx8obNxSFBZuXIlWVlZZGZmNt8ma9lvLSsri7S0NO6//34yMzPJyspi8eLF5OXlkZmZ2fzavLy8E473Pd/8+fNPeRsuJyeHnJyc5nM1ffvPysoiLy+PNWvWNMfaWowAbre7OUbfY1vG3zLeQDtQXMVZsQcJqziqVWtDmHIvAI5uZqCnlIKkswLzftJrKP/8LyRbwjuqdscajpskEsdqtbCe9PMXt7KtoMyR956WPJB7rpre4ePvv/9+MjIySE21p6tsSsymfmu5ubmA/a178eLFrF69mjVr1jQfu3DhQtavX998vuXLl7NixYrm45v2NZ1v1apVbcaycOFCVq9ezYIFCwBYunRp87lycnIAWLRoUfPxLWMEmDdv3gnxLFy4kMTExFbjLy4uZsWKFSxfvrzDfy9/OVBSzZUxm6DaggnzA/7+fZVyT7nXUuOeNwkHqkZq9qGepvxT/qklvC2N9YTve4d3GmcyNXmg09FIkEhJSeGuu+5i5cqVuN3uExK9tX5rLpfrhP2+x7hcrm59u16/fn3zBTElJaVD5/J9/8zMTFJSUk7Yf9NNN3HfffedEKPva51qCT9YXMU5nmwYdTr0S3IkBnGWci8wudew5012ekbSb9CIgLyf9A7KP7WEB9bBbCLqK3jHM5OfD+7vdDR9Wme+jTut6Zv3ihUrWLx4MYsWLWLFihVtHt/y4uSb2P6wdOlS5s6dS3FxMcXFxZ167bp161qNz/di05kBMT3F4zHUlRxiTOQumPQzp8PpU5R7Xdcnc6++hsiDH/K+5yJG9I8O/PuHGOVf1/WV/AvJlvAOjRLPXYuHcN43M0jsFxWw2CS4ZWVlsWDBAtasWYMxhuzs7B77hnyqC4vb7SYtLY1ly5axYMEC0tPTT9jnKzMzs9VzjB8//qT3cLvdJ7UQ+FNXZkc5Wl7DOWywn0zSVKGhSrnXPR3KvQMfE9ZYy3ueGQxSES4+lH89IySL8OTkZIwx7RThb3Cg3zQi4hKI0Mwo4rVmzZoTLjy+I66Li4tbnRqpo1JSUk64MDT1R2vt/NnZ2bhcrubWhaaY8vLyyMvLIyUlhaKiopPew/ccixYtau4/12TVqlUsW7asy79DezqUey0cKK5mXlgONXHJMGRaj8UmwU251z0dyr2IaA4NvZCPPVMZrCJcfCj/eoa6o7TltkwefvItktxqBZfPjB8/nqysLFJSUnC73cydO5eUlBRycnJYvXo1eXl5ZGVlkZiYeNLzFStWNI8QT01NZcWKFWRnZ7Ny5UoWLVpERkYGq1evbv72Pn/+fJYuXcrKlStJT08/4XwZGRmkp6ezcuVKUlJSSE1NJT09nczMTJYsWUJqaiqrVq1q3g+cFGNGRgZr1qxpvq2Xl5fH4sWLSU1NbfX38Y2/6dZkIBw6VsKlYVuoTbmZGMsK2PtKcFHuBSD3Rp/J81P/h8p9Oxk0QJ998hnlX8/kn2WM8esJe4P09HSTnZ3d7nEL//4B4WEWTy7SVE0i/pCenk5Hcs/Xs6sf5rqt36b+5lVETlF3FJGu6Gju/eLFbaxat5+tv1CuifiDZVnrjTHpre1TP4tTKKqoU784EYcNOvw2NUQROf4Cp0MR6fOOV9QyaIA+90QCQUX4KRyrqFURLuKw8aUfsy16llbJFAmAY+X63BMJFBXhbahtaKS8poEkzYwi4pySfSQ3HuJAwplORyISEo5X1DKovz73RAJBRXgbiivrAHRbTsRBdbvWAlA+4nyHIxEJDccrahmszz2RgAjJIrwj86UeL7eLcLWEi/hPZ+cJr92ZRYFJZMBITU0o0h0dyb36Rg8lVfXqjiISICE5RWFycjIFBQWnPOZ4ZS0ASboYifhNR3KvmaeRmAPv8kpjKpMHadVake7oSO4VVXjvAOtzTyQgQrIlvCOOl9tFuBYsEHHIoRwi68t413MaoxPjnI5GpM87XmF/7qkIFwkMFeFtKPL2CU/SABURZ+S+gQeLjZGzSYiLdDoakT7vmLcIH6yFekQCIiS7o3REUUUtMZFhxEWFOx2KSGjKfYO9UZMY2G8ollbKFOlx9Q0ehgyIZnD/GKdDEQkJKsLbcNy7UI8+/EUcUFMKB9fxQeT1jElSVxSRQLhk+jAumT7M6TBEQoa6o7TheEWtBmWKOGXv+2Aaea1qCqNVhIuISB+kIrwNC9JGcvuZY5wOQ4JcVlYWaWlpLF26tN3jVq5cSVZWFvfffz9Lly4lMzOThQsXOhJP0Mt/BxMewycNEzQoU1ql3BNxjvLPP0KyCO/IfKnXzB7BgrSRgQtKeqWMjAwWL158wraFCxeycuXK5uc5OTmsWLGCRYsWkZGRQVFRETk5OSxYsIC5c+f2eDzBpMPzhOe/Q+ngNOqIZExiv4DEJr2Lcq9zOjtHv8ipKP/8IyT7hHdqrmKRTlq8eDEpKSnNz7Oysk644IwfP77530uWLAlobE7rUO5VHIPCreyf/G0A9QmXDlPutU2fe9LTlH+dF5JFuEhPysjIaPeYpKSkAETSS+19F4BNkTOJCLMYHq+ZGqRjlHsizlH+dZ6KcHHeqz+EI5udee9hp8Hlv+nUS9xuN/fdd1/zN/zc3NzmfTk5Odx1111kZGSwfPlysrKyWLNmDQAulwuANWvW4Ha7cblc3HjjjWRnZ7N06VJuuukmUlJScLlcrF69mqVLl5KTk9P8mqVLlza3MjRtz8vLY82aNSxfvrz5/E0yMzNbfW3Q2/cBRA3gk9oxjEioICI8JHvNBYZyT7knzlH+hXz+qQgX6aR58+axdu3a5sRft25d877U1FQWL17cfHHKyMhovmgsWrSo+bjc3Nzm50192VavXt180crLy2PhwoWsX7+++TXLly9nxYoVgN33bvXq1SxYsACApUuXNu8D+zbg8uXLASguLmbFihXNz4PeZffB6YvYu6pQgzLlBMo9Eeco//xPRbg4r5Pfxp2UlZUFcMI37/Hjx5/QIpCYmHjC847yPWdiYiKJiYkn7MvLy2t+vn79+ubjU1JSTtjXtM33XE0XuF4hPBIGT2J/8V4+d9pwp6Pp25R7J50zpHNPAkv5d9I5Qy3/VISLdEJOTk6P3dryvfAAJ91ia2np0qXMnTuX4uJiiouLT3mu3qa0uh53Vb0GZUoz5Z6Ic5R/PUOdLUU6ITU19aRv3oHmdrtJS0tj2bJlLFiwgPT09BP29QX7iioB1B1Fmin3RJyj/OsZagkX6YSm0d9Ng0uAE/qugd0PrScvCNnZ2bhcrub3b7ow5uXl9doLUUt7CisAmDCkv8ORSLBQ7ok4R/nXM1SEi3TS2rVrue+++5g/fz5g3zrLzMwkMzOTlJQUVq9eTV5eXnMfulWrVgF2S0JiYmLz/vHjx7No0SJycnJOeE1iYiIrVqwgLy+PzMxMUlNTWbFiBdnZ2axcuZJFixaRnp7OypUrSUlJITU1lfT0dDIzM8nIyDjluZoGswS73YUVRIRZjEnSQj3yGeWeiHOUf/5nGWOcjiHg0tPTTXZ2ttNhiISc9PR0OpJ7dz2STf7xSrK+e0EAohLp+zqaeyLiX5ZlrTfGpLe2T33CRSTo5BZWMFFdUUREpA8LySK8oKAAy7K49957nQ5FJKR0JPdqGxrZW1Sp/uAifqTPPZHgE5J9wpOTkykoKHA6DJGQ05Hc23u8Co/RoEwRf9LnnkjwCcmWcBEJXpoZRUREQoGKcBEJKrsLy7EsGD9YRbiIiPRdKsJFJKhMGTaAO84cQ0xkuNOhiIiI9JiQ7BMuIsHrshnDuWzGcKfDEBER6VFqCRcRERERCTAV4SIiIiIiAaYiXEREREQkwFSEi4iIiIgEmIrwEKAV0nqO/rZyKvr/o+fobyunov8/eo7+tv5jGWOcjiHg0tPTTXZ2ttNhBIxlWYTif+dA0N+2c9LT01HuiT/ob9s5yj3xF/1tO8eyrPXGmPTW9qklXEREREQkwFSEi4iIiIgEWEh2R7Es6xiwz+k4AigZKHA6iD5Kf9vOSQVynA4igPT/R8/R37ZzlHviL/rbds4YY8zg1naEZBEuIiIiIuIkdUcREREREQkwFeEiIiIiIgGmIlxEREREJMAinA5AeoZlWYuANGC1d9NCYLkxJs+5qHony7JcwCIgyRiztJX9S4A8IBHAGLMyoAFKUFHu+Y9yTzpDuec/yr3AUBHet92InUQ5wF26EHWeZVkZgAsY38b+5cA6Y0xm03PLshY0PZeQpdzrJuWedJFyr5uUe4GjIrwPM8YkOB1Db2eMyQKwLGsu9kWppUUtWglWAcsBXYxCmHKv+5R70hXKve5T7gWO+oSLdJFlWamtbHYDGQEORSSkKPdEnKHc8y+1hPdh3v5xxajPVk9JxP77+mr5XEKQcq/HKfekVcq9Hqfc8yMV4X1XNuBu6g9nWdZqy7KK1WfLr1xt7bAsy2WMcQcuFAkiyr2e52prh3IvpCn3ep6rrR3Kvc5Td5Q+yhiT02JAyjpgmVPx9FFuvK0tPlo+lxCj3AsIN8o9aUG5FxBulHt+oyK8j/KObvaVB7TWl0u6rpiTWwVcAGoNCF3KvYBQ7slJlHsBodzzIxXhfZBlWSnAGu88n740VZMfGWNysFsFfCUCWYGPRoKBci8wlHvSknIvMJR7/qUivA/y3o5b2uJb6U3YUwiJfz1lWdYCn+fzgRVOBSPOUu4FlHJPmin3Akq55yeWMcbpGKQHeFsFmpIkCcjVKPHO807HlAEs9m5aAWR5WwOajlmCvTBECmg0fqhT7vmHck86S7nnH8q9wFERLiIiIiISYOqOIiIiIiISYCrCRUREREQCTEW4iIiIiEiAqQgXEREREQkwFeEiIiIiIgGmIlykHZZluSzLWuF9pAT4vZd433dB+0eL9C3KPRFnKPcCI8LpAER6CbcxZmnLjd65UpOAIuxVxIqNMZnei0dWd5fxNcbc710BbhmQ2Z1zifRSyj0RZyj3epiKcOk2y7KWAy5jzOJ2D27/XK6WCezP8/uL9wKxGlhhjLnfd7tlWYuApUCaQ+FJiFDuKffEGco95Z4/qAgXf1jlx3PdCLRcecuf5/eXtdhLJGf5bjTGuC3Legot4SuBodzzUu5JgCn3vJR7XaciXLrNdylbP5hPi4uRn8/fbd5v/MUtL0RNvBekPn0LTYKDcu9Eyj0JFOXeiZR7XaOBmRI0mm6/OR1HByzGviV3KsHYiiHSKuWeiDOUe6FNLeG9WNM3U+/TRGPMyk7sy/M+dQGJQDawHHuQxX3ebS5grjFmqWVZGd7jU4E8Y0ym91wpeG9BGWPmd+T92/hdFnjfL8U76APsloFE3/NblpXalTh93mcJkNP0e7cXVxtS+ezv1yrf9/XpK+cGFhpj8rzb13jjv8t76APY/x1WAynAfGPMwi7EJz1MuafcE2co95R7fYoxRo9e+MBOyEU+z1OanuP9n9lnnwtY4/33olb2Lff+OwNY32L/6qb9PttKWjxPbTp/e7G18zudcJ5TnL+rca4GUlvEuaADcfn+jVIA43ueDv73WoQ9mMV32wLsgTcn/V7e91zS8v31cP6h3FPu6eHMQ7mn3OtrD3VH6YW8I5QXmRO/zS4Axnu/MacY77dOsPtqAXneb6Zg31by3dd0C6kYOzl8v+229s232BtDE3dHYuvAr9YWd4vnnY7T23KRak7sZ7cKn79FR7R4z868biX24JuW290+T5t/L2OM2/iMPpfgoNxT7okzlHvKvb5I3VF6pwxaJF/T/7gtbrn5ygXSjDGLLctabVmWAbKA1S0uHO5WXlvkj9j8zN3KtlPFmQG4fW7bgfc2YBfeOw9Ix76916zpQox94V2EfVtxuc8FLMuyrAXGnk/VReu/Q5cudhIwyj3lnjhDuafc63PUEt73uE6xL9E7H+lCIAH7ttRCy7J8pxUqbuV1bv+F13HtrNLV2Thd2H3lsnwemcaYrrRUZAIn9Vnz+RbftLjB8hYtCPfxWQtEhml9lLm7C/FIcHCdYp9yT7knPcd1in3KPeVe0FIR3jvl0Mo3We+3zKzW9mF/S12DvQJVU+JkGXtQiT+XpD1VbJ2V2u1oPtNqXF1h7BXEUry3QE/lhAum95ZgSjsXWQluyr3OU+6JPyj3Ok+5F+RUhPdC3m+ZT/n0dWuS4f0fPs83UbwXgvSm22+tvM73W2tiK2/p8kdsHXh5HideMNynOLZTcXq/fRe3vIC0EmdHzQeWe0e3t3Qjbd9eW+F9tDUHrKuL8UgAKPcA5Z44QLkHKPf6nvZGbuoRvA9gCXYfrAxajHb27lvgfSzBOxrZ++8Mn8cC7AtAKvYo6hI+G528ALtP3XrsC13T6w12QqV4H02vW9SR2Dr4Oy3wPj/h/F2Ns7XzdzQuTjFK23u+5T5/70U+f09XG+dqayR88+/l+9pTvb8ezjyUe8o9PZx5KPeUe33pYXl/WRFpg7dFZZmxb8f543wLTIt5XAP5/iK9hXJPxBnKvcBQdxSRALIsK6MzFyIR8Q/lnogzlHttUxEu0sMsy1rh0yfP5WQsIqFEuSfiDOVex6gIF+l5q7GnyerU7TgR6TblnogzlHsdoMV6RDomxbKs1cBS08nVw0zr86J2iGVZS4C5fLa6m0ioUe6JOEO518M0MFNEREREJMDUHUVEREREJMBUhIuIiIiIBJiKcBERERGRAFMRLiIiIiISYCrCRUREREQCTEW4iIiIiEiA/T/VtM5A1nDDTQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "rng = np.random.RandomState(999)\n",
    "ih = rng.randint(0, n_halos)\n",
    "\n",
    "fig, (ax0, ax1, ax2) = plt.subplots(1, 3, sharex=True, figsize=(12, 4), sharey=True)\n",
    "fig.subplots_adjust(wspace=0)\n",
    "\n",
    "for ax in ax0, ax1, ax2:\n",
    "    xlabel = ax.set_xlabel(r'${\\rm cosmic\\ time\\ [Gyr]}$')\n",
    "    yscale = ax.set_yscale('log')\n",
    "ylabel = ax0.set_ylabel(r'$M_{\\rm halo}\\ [M_{\\odot}]$')\n",
    "title = ax0.set_title(r'${\\rm 200m}$')\n",
    "title = ax1.set_title(r'${\\rm 500c}$')\n",
    "title = ax2.set_title(r'${\\rm vir}$')\n",
    "\n",
    "__=ax0.plot(tng_t, 10**log_mah_200m[ih, :], label=r'${\\rm simulation}$')\n",
    "__=ax0.plot(tng_t, 10**fit_data_200m['log_mah'][ih, :], label=r'${\\rm diffmah}$')\n",
    "__=ax1.plot(tng_t, 10**log_mah_500c[ih, :], label=r'${\\rm simulation}$')\n",
    "__=ax1.plot(tng_t, 10**fit_data_500c['log_mah'][ih, :], label=r'${\\rm diffmah}$')\n",
    "__=ax2.plot(tng_t, 10**log_mah_mvir[ih, :], label=r'${\\rm simulation}$')\n",
    "__=ax2.plot(tng_t, 10**fit_data_vir['log_mah'][ih, :], label=r'${\\rm diffmah}$')\n",
    "\n",
    "for ax in ax0, ax1, ax2:\n",
    "    leg = ax.legend()\n",
    "fig.savefig('example_tng_fit.png', bbox_extra_artists=[xlabel, ylabel], bbox_inches='tight', dpi=200)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "pressing-island",
   "metadata": {},
   "source": [
    "## Calculate fitting residuals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "likely-pantyhose",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats.mstats import trimmed_std, trimmed_mean\n",
    "\n",
    "\n",
    "def measure_fitting_residuals(tarr, log_mah_sim, log_mah_fit, lgm_min, dlgm_min, pcut=0.01):\n",
    "    assert tarr.size == log_mah_sim.shape[1], \"shape mismatch between log_mah and tarr\"\n",
    "\n",
    "    nh, nt = log_mah_sim.shape\n",
    "    avg_error = np.zeros(nt)\n",
    "    std_error = np.zeros(nt)\n",
    "\n",
    "    gen = zip(range(nt), tarr)\n",
    "    for i, t in gen:\n",
    "        log_mahs_at_t = log_mah_sim[:, i]\n",
    "        log_mah_errs_at_t = log_mah_fit[:, i] - log_mahs_at_t\n",
    "\n",
    "        low_mass_msk = log_mahs_at_t > lgm_min\n",
    "        dlgm_msk = log_mahs_at_t > log_mah_sim[:, -1] - dlgm_min\n",
    "        msk = low_mass_msk & dlgm_msk\n",
    "\n",
    "        if msk.sum() >= 10:\n",
    "            avg_error[i] = trimmed_mean(log_mah_errs_at_t[msk], limits=(pcut, pcut))\n",
    "            std_error[i] = trimmed_std(log_mah_errs_at_t[msk], limits=(pcut, pcut))\n",
    "    return avg_error, std_error\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "australian-james",
   "metadata": {},
   "outputs": [],
   "source": [
    "mu_200m, std_200m = measure_fitting_residuals(tng_t, log_mah_200m, fit_data_200m['log_mah'], 10, 2.5)\n",
    "mu_500c, std_500c = measure_fitting_residuals(tng_t, log_mah_500c, fit_data_500c['log_mah'], 10, 2.5)\n",
    "mu_mvir, std_mvir = measure_fitting_residuals(tng_t, log_mah_mvir, fit_data_vir['log_mah'], 10, 2.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "excessive-encoding",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEQCAYAAAB80zltAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA6gUlEQVR4nO3deZxcZZ3o/8+39qreqruzsyUdQBCQpOmo4IJKt8I4jgpJuMM4KChpHddxSYPMjMx1NCRy7w/nitjB13gHRwU6orhDOojr/JROo4K4JKmAhJCll+p9qeW5f5xzitOV6kov1emu6u/79apX9znPqVNPVXWf73l2McaglFJKFZJnvjOglFKq9GhwUUopVXAaXJRSShWcBhellFIFp8FFKaVUwfnmOwMLwZIlS8zq1avnOxtKKVVU9u7d22WMWZorTYMLsHr1ajo6OuY7G0opVVRE5NnJ0rRaTCmlVMFpcFFKKVVwGlyUUkoVnAYXpZRSBafBRSmlVMFpcFFKKVVwGlyUUkoVnAYXpZRSBafBRSmlVMFpcFFKKVVwGlyUUkoVnAYXpZRSBafBRSmlVMFpcFFKKVVwGlyUUkoVnK7nopRS86izs5P29nYAHn/8cZqbm2lsbAQgHo+zc+dO6urqiMViNDY2Ul9ff9K0hUCDi1JKzaP29na2bt0KWAFjzZo17Nmzh/r6ejZt2kRrayt1dXUANDU10dbWRjQazZu2EGi1mFJKzZPOzk62bduW2Y5GozQ0NNDe3k48HicWi2WCB0BdXd1J0xYKLbkopUrKv3739zx9uH9eXvulqyr51FsumPLx9fX1tLW1TdgXi8WIRqN0dHScUAqJRqPs3r2baDSaN62lpYVrr72Wuro6otEobW1ttLS00NnZCcDu3btpaWmZEJwKTYOLUkrNI6d9BazA0tPTw+bNm2lvb6empmbCsbW1tcRiMeLx+KRpjY2NNDc309bWxu7duzPn3bRpE3v37s0cv337dlpbW+fsfWlwUUqVlOmUHBaa5uZm9uzZkymV9PT0THpsvjRgQsmmpqZmQjCKRqPEYrFZ5fVktM1FKaUWgB07dtDS0pLp8RWNRonH4xOO6e7upqamJm+aI7tkc6ob+rXkopRS82zXrl0TuhLHYjEaGhpOKJ3E43Gamprypi0UWnJRSql51N7eTjQanTB+pbOzM9NzzF191dHRQWNjY960hUJLLkopNU9isVjO0obT8N7W1pYZKNnT08M999yTqd6aLK2zs5O2tjZisVimU0BrayuxWIxdu3ZRX19Pa2srHR0d7Ny5ky1btszJexNjzJycuJg0NDSYjo6O+c6GUkoVFRHZa4xpyJWm1WJKKaUKToOLUkqpgpuT4CIiV8/FeZVSShWHWTfoi8hNQAvQ6+wC1gPe2Z5bKaVUcSpIbzFjzNnubRG5phDnVUopVZwKUS2Wq5vV3hz7lFJKLRIFKbmIyP2AM5pHgCuADYU4t1JKqeJTiODSDDwAxF37ogU4r1JKqSJViODSZozZ494hIjoiUSmlFrFCtLkYEXmDiKx2HsDNBTivUkqpIlWIkstOrAZ8ce1bD9xSgHMrpZQqQgVpc8lRLba+AOdVSqmSt2PHDrq7u7n22mvp6emhra0ts0JkPB7PTE7prDLpnj15srSFoBDB5YBdFea2GXiiAOdWSqmSt3PnTnbu3EljYyP33HNPZv+mTZtobW3NrHXf1NREW1sb0Wg0b9pCMKPgIiL7gXpjTD/QjlaLKaXUjESjUXp7e0/YH4/HicVimeABUFdXR3t7O42NjZOmbdy48ZTk+2RmFFyyRuRvMsZMKKXMRbWYiGzFGktTY+dhZ55jo4CzSMEGYHe+45VSJeYrbz5x3wVvg5ffBOPD8LVNJ6avuw7W/x0MdcMD15+YvuFGuPAa6DsEDzafmH7ZB+AlV804y84CYU7A6OjoOKEUEo1G2b17N9FodNK0jRs3EovFaG1tZcOGDUSjUWpqak55ldmse4u5A4vdW+wN2cFmtkRkOxAzxuyyg8RaEckXnm8xxuywH5uAFhGZmxVxlFJqlnbt2kVdXR2dnZ20tLQAVsmlpqZmwnG1tbX09PTkTQOriuyWW25h48aN1NTUcNNNN52aN+JSiIkrtwG7sarFbse6kL/HGPPl2Z7bZYsxpsW1fT+wHdiVIz9RoC5rdyvW5JpaelFqMbjh+5OnBSL508tq86dXnZ4/fZrcK0Fu3LiRlpaWzOqUTrDIZbK0Xbt2TSjZ1NfXs2fPnpzHzqVCjHNpN8Y8ilUNdbP9+8ECnBcAEclVlosD+RaLbhQRd4CJc2LAUUqpedfZ2Tlhu76+PlP1FY/HJ6R1d3dTU1OTNy0Wi51QqpmPRv5Jg4uIVE7xHFUisga4xNUluWrWOXtRDZAdoicN58aYuDGm2hgTc+1uwup4oJRSC0ZnZydXXHHFhH3xeJy1a9fS0NBwQukkHo/T1NSUN62+vj5viedUyVdymWpvryew5hfbJCJVIvIJCltKiE6WYFeB5WUf04hVLZbT4cOHEZHM47bbbpt+LpVSaprq6+vZvn37hH2xWIzNmzcTjUZpaGggFnvxPrmjo4PGxsa8aY2NjZnzOHbuPPUtAmKMyZ0gsg+rrePHpzZLJ+SjEWv+smrXvjrgAFBtjImf5PltQKsxZtKSS0NDg+no0OnQlFKnXmdnJ+3t7USjUQ4cOMC1116bc6BkT08PDQ0NU07btm0bGzZYk9PX19dP6LZcKCKy1xjTkDMtT3DpBZyqsU6saqXdQIc9viX7+Mpc+2fLbnPZa4yRfPsmee5WoDNfYAENLkopNRP5gku+arF64L1Yo+0fAC7BCjC9IrJPRO4Wkatdo/PnZNCkMaaTidP5g9UOkzdg2F2VM4HFLgEppZQ6BSYNLsaYg8aYe7Au7DFjzBuNMR6sQYk7gFrgy1jTv3Tz4qDFufBA1riWJqzuxYBVTeZOtwNJDdAhIlG7Gm3hTLqjlFIl7qTjXJweYCJyDXDALkl0AvfY+9dgXew/MVeZNMY0i8hWO2jU2flwj3HZaOdhl92Av9ve3+o65oQxMUoppebGpG0uOQ+2Asl6rLEt/VlptxtjinIdF21zUUqp6cvX5jKtEfrGmIPAQRG5QkSqjDEPupJbJ3ueUkqpxWWmE1e6q8p6jTGP2oFHKaWUmvn0L/YI/gNAnYjcn2NNF6WUUovUpCUXO3jUuR5rs7bBmqwyjjUVfjO6hotSSinyV4vFAWP/PIgVQJ7A6nUVw+qerFVhSimlTpAvuMSwGunjWL3DNJAopZSaknzBZZcx5nNgrSwpIu9xpcXsqfUzROTqrN5jSimlFqlJg4t7zIq9sqR7xck1InIT1ozFBugDtgIaXJRSSs24K/JB7BH6kJlIsrZQmVJKKWXZtGkTTU1NE1asLAazXuYYrMklRUSXEFZKqQJrbm6ek+ny51pBggtMrEZTSilVGM7iX8WmYMFFKaUWiht+dMNJj7n89Mt514Xvyhz/1rPfytvOfhu9o7189LGPnvT52ce/84J38rozXjetfO7atYuWlhbq6upoa2sjGo2yadMmOjs7aWtrA+Cmm26isbGR7du3097eTnNzM83NzUSjUVpbW9mzZw/RaHRar3sqaHBRSql5snHjRnp6eti9e3cmQDjVYE5VWHNzMwcOHACsUkxzczP3338/e/fuBViQgQU0uCilStBXrvzKjI+vDlVP6/nTPT7bli1baGlpyWzH4/EJbSw1NTWZ4OJw0hdyI/+M5xZTSilVGI2NjezaZS05NZWSyIYNG+Y4R7M365KLPQdZI9YKlVGsEf2Pk2PNF6WUUidqbm5m+/btgFVVVgpmVXKxp9zfgTWpZQxrhcqYvb1dRK6edQ6VUqrENTY20tHRQU9PzwlpPT09xOPxU5+pWZp1ycUY897J0kTkitmeXymlFoNbbrmFhoaJizo6vcZisRjt7e0A3H///YBVfbaQ21ymtczxCU8WucYY88086UUx35guc6yUUtNXsGWOc4iJyANYi4Z12/ucaWDq0KWPlVJqUZpVm4sx5gljzGagA2vhMAF6gA5jzLXZMycXi+7uboaHh+c7G0opVbQKNbfYpFVjxSiRSDA6OkokEpnvrCilVFGa03EuxdxbbHh4mEQiMd/ZUEqpojSrkouI3A5UTZYMXEKRrvGSTqcZGBigpqZmvrOilFJFZ7bVYs6cBLHZZmSh8Xq99PX1UV1djYjMd3aUUqqozDa4PAA0GGP25Eos5ouyx+MhlUoxPDxMWVnZKXtdYwypVIp0Oo0xBmNM5vdUKkUqlSKRSJBKpTLpDhHB4/Hg8Xjwer14PJ7Md+D8dNKcdK/Xi4gU9XellFp4ZhVcjDF9QM7AYqdPmlYMPB4PR48eZdWqVYRCoYKdN5VKMT4+TiKRYHx8nLGxMZLJZCaoTHahdwJJvmDgHDPZ+KXs57nPmZ3u3ud+OEHLCUxOoHIe7mPcP53jlVKlT2dFzsPr9ZJKpTh8+DArV64kHA7P6DzGGEZHRxkeHmZoaCjTUcAYc8JF213aOFXyDaR1Ske5gtZkgSxf4PN6vfh8Pnw+34TSU3aAylX6UkoVj2kFFxG52xjzvjzptwPrgd3GmDtmm7mFwB1gli9fTnl5+UmfY4whkUgwMjLC8PAwIyMjmQu0cyFdSBfMfHkpZD6dzyCZTJJIJE4ITrley/2ZOYEpEAhkfs8OSPlKZgvpMz8Zd1B3qkWzf2aXKLOPz1d6zVVSzXVMvhJtLs53MN3Pu5i+GzU10y257AIQkXVAzD3rsYhsA7qMMW8SkSuKZeqXqfB6vaTTaY4cOcKSJUsIh8P4/X48HqsntxNMnNLJ8PBw5p/bXX202M3kAu9cIFOpFMlkktHR0UmD0cmmMsp10cs+h3s7V8ks+7nZJc/strJ8eXCe4xzrDh7u5xbibyfXZzNXf5MzmVIquxrVXap1/8z30P+xhWXa1WIist/+tVpEtrlKKBuNMeeA1dZiz5hcMpxA0tXVlblL9Hq9+P1+xsfHJwQT559Bzd7J7rCzTSfw5Lq7z3e3n338dKsFs4+dLGgttoul+3NwgmoymZyQlq9066Rntwl6PJ5MSdfv92eqY50Sr5pb0w0u9cAldkM+InKTiFTaJZjsbz1egPwtKM4/Pbx4YRofH19UF4KFLl81z1SPne651exMtcotn+yg7ZQAx8fHGRoamnBud9ufO/BkV7/q9z070w0uMSew2B4AGoBHgd6sY2c+3XIRKLY6fKVK2WQBKlcJxQlETm9Np9STfb5gMEg4HCYYDOL3+/H7/fo/Pw3TDS41IvIGrIkqa4BmYJuIVAHVWcfWYQUdpZRaMJwAka/q2hjD2NgYIyMjEwKKz+cjGAwSDAYJBAIEAgEt5UxiWsHFGHOPiHwJq2H/ALATaMJa5rhZRD5upzVSgqP2lVKLg9OpwB2AnKq2oaEhBgcHJwSUQCBAKBQiEokQCoW0zZUZNOjbK09mrz75TQAR6QFuBh4p1un2lVIqF3eHATent+j4+Dj9/f0YYwgEApSVlWV6li7G0s2MBlGKyHuwJqWswRrT8mWw1nfhxMCjlFIlyynlOJypmnp7e4nH45n9TummrKyMUChU8j3WZtIV+RGsKq8D9s8GEWkGrnCPe1FKqcUoVwknu3QDUFZWRkVFBeFwuCQDzXRH6L8H2JTVYwwRiQJbgJIYla+UUoWUq3QzNDSU6SZdXl6eCTSlUn023ZJLb3ZgATDGxEXkYIHypJRSJU1E8Pmsy68xhsHBQQYHB/F4PFRWVlJWVkYwGCzqQDPd4JJv7EpJj2tRSqm54A406XSaeDxOPB5HRAiHwxM6BhST6QaXWhFZbYx5xr3TnmtsbaEypZRSi1H2LCDOTOrOoM5oNEokEimKNpqZjHN5QETW8OI4ljqskfvXFjx3Sim1SGVXnY2Pj3P06FEAKisrqaioWNBVZzMZ57JZRNZjTfsSBW63uyArpZSaA+4OAcYY+vr66O/vx+PxUFZWlqk6W0glmhmNc7GDyYSAIiIfL5U1XJRSaqESkUz7SzqdZmBggIGBAQAikQjRaJRQKDTvJZpJg4uIPDyN8wjWoEoNLkopdYrkaqMZHh7G5/MRjUapqKiYt9JMvpKLAC1Mbep8AW4vRIaUUkpNn9NG48wQcPz4cbq6uigrK6O8vJxwOHxK5zzLF1xaptOWIiItBciPUkqpWXBPuukM1hwcHAQgFAplBmzOdaCZNLhMt5HeGKODKNWiNZxIMTKeJhLwEvLpWj9qYcjV46yrq4uenh6qqqqoqqrKpBfa3JxVqQVgJJHmQPcoY8k0ybQhmTak0gavCH6vEPB5CHiFkM9DddhHZciL12MFhVTa8HzfOPu7R9nfNULPSJKA10PIJwR9Hrwe4dhAgkN9YxzqG6d7OJl5Xa9AJOClLOAhbWA8mWY8ZRhLpSkPeDl3aZhzl4Q4Z2mYNdVBEmnD4FiKwfE0g2MpBlwPZ3s0mWY0aRhLpBlJpqkJ+1h/WhmXnF7OhSsiBH3569XTxtA9lKRvNMkZ0eBJj1elx12icQ/WdIJMoQdpFk1wEZGtWGNragCMMTsLebxaOIwxxEdSmQv3833jHOobI5mGi1ZEeNmqCOcsCePzTCwdpI0h1j3Kr58b5PHnBnnyhSESaQMIyDjiHcUkywAv1oQSWasPAlUhL5UhH8cGxxlNWpNOeD1QHUkznoJEwsdoMo0BomFhZXSUl5wxSHlZPz5vEtIhUskIqWSIZCJCSKoJeD2kvF14PUkSY8vY1zXKfb/tIpWe/DMQoCzgoSLopSLoJez3UBn0Ei73E/R5ONw/zn2/7eJrT3QR8ArnLwtTGfIR9AlBr4egXxhJpDk6kLAegwmS6Rffz9m1YS5YEeaMqiA+j+D1gNcjeEXwCHg8glesC5JHmLDfI1bwHU8axlKG8WSaRNpQHvQSDXmJhn1EQz4iASt4u0txg2MpnouP8Zf4GM/FxxlLpYmGfETDXqpCPqJhH+VBDxUBL+VBrwbBOeJ0BDDGZIJMKBTKTD1TiCozca89vVCJyHbgcWPMrlzbsz2+oaHBdHR0ZLaPHDnC8PCwLvgzx5Jpw+G+cZ7pHeNgzzAH+o7z3OBRjo0cIyndlPmPUTtawUsHqzgr3MV9qzoYO/pG4n2vIeRLEwyMk0pDIm1IpiGVFjyBLrzhZzit6mnGgoe4LbmGc7xR9nh7+Xf+zD+f+wWq/cv55cHP89DwfxPET1CChKUCr6kmaZaQTpYT8vUjvsPcWv1mzvL28eX+/58HR/bx4OXfxhhD65++yA+PPEyKySNEja+C+869GU9imH869HW6zDh3vrIVjKH16f+PvmQAf/o0vB7wySB4kqS9adLpARIjz7PGX8M7KuuRdIIPHv02F1ZfzDvO/yAy0sPNv72VVRXnEOYC4r1nEjsWYDiRti/4aUYTaYI+Dysq/CyvCLC83MvyCj9VoQD7u0b4/dER/nhsOBNA54qAFfDsINE3msqkeQV8XmEsTx4CXqEm4uPcJWHOWxZmabk/U3oM+Tz2T8HrEfpGU/SNJImPJjOvs6TMz5KIz/pZ5iMSyP8/nUobPMKiq9Z0FkJLp9OZaWeqqqqIRCJ5PwsR2WuMaciVViwlly3GGHeHgfuB7VirXhbi+JLnGY3jHXqBVGQZ6VANnOSfRxLD9D/9FZ4ZOMCrQmdhQlH2SZoj4SWcf+ab8ZgU/cf3kvT5KI8sJxJZhXhe/Mfd393FY/sfo7frGVYfMwwT4omawwzIaiLpJtKjcfrD94Kvl0RggCHfOMYLVIFUgR8IpNP8Vf8AHxnrwyTh+HAt5196NsnUGTwf+0/uDzyCAAH74RYwPi4d6mP1cA9nJZJcIUnKq5Zx0bJlRHxhIh2/pWKsl2GP0O/x0OX1ciwQoTcQYphhytKGs4bGWRP7MUtTad4YCrJq+YXWZyPCRX/6FsvNGKuSKVYmk6xKJknVvZnnL7qBwUQf4d3vZVS6WLnv7wH4YDBAb91V1ncx3k/Ps99nbyjIaFY3UZ/4CHsCREf7WDIyytLubwBwYW01p4frABjrO0Blz5/5+fCzDHmsNflWLolQG15OyF9JBKF8fJg3hdfwilSAZ4ee5V2jv6Pl9Pdx6RlvZrXvJ+xLPMDVZ6yhJlKHeCOItwyPJ4gYq2rEGKGu7DzCnkriiV4Oj/yFM8Ln4iNIfyLOQPo4hnHSMo6RBCmTgHQQk46QTIQZHw9jkuWMp2DMrtJLG8OqygCByAuMep9l2BxhPD1G0BMGE4R0iHQqCKkQ6VSEivTZDCUMLwyM8+djo/z04OxX9Aj7PNSW+aiN+Kgt8xP2e+gZTtI9lKB7OEnvSJKqkJeXrSzj4lVlXLQiQnXYZ1Wjeq3qVKfEPJJMMzKeZmg8zXAiTXnQw4ryAD5v8QWm7E4AIyMjmS7N1dXVM+rSPOPgkmuOsbkgIvU5dsexllKe9fGlzDt8jO59u+g78kte88JTiEnz79VVVK77MK9fvZHIM4/gGR9kqO6vIJ0gdPhXHE/28QMfPPz8YxxNWH00On+3B78YHq6t5nuRao49vIYoA7zmzFt5rCwCgM8YqlIGY6qIEyHtt6apOCuQ4Au8AMDfB5fjkSFSo2/kNcED/DL8NGVpw2ljSVaMBVjhKyfykr+jetklrBgdorb7jyTLVvB82QpSZct5XzpJOlgFwHjPMGd0eTHiwYhgEIx48K/7ABdUXcCK8RHS/gjpcC2HsCrCLjMmE1TPuvzfWZ0YRJKjeBJDeMb6MN4go6e/mvHUONGDjyAiJC9ezqGylazwhVhpUjj33Vee/U7S/nJSkaWkIsusoB0op87jg3SK4OvuxjM+wFGPD+OLsMQfoTq8hDRgPD4+s/Z9MHiErtEX8PrKCYRqkFWXkV56EZIcxd/7Z4wnwGGvHyM+3iVCKlSNAQJLXsa/rvssnq4nebb7CZ4cOsjTHKcrvJKh9DgDI0f5y/Ax1j/3K6L9AxCK8vfRKKd5ywHw9PyJ8sEjPJroYah/8n47dw942ZAO8nDQw5dDg9zdcBdnlNfwvSe/xJe6f3zSv7+d6/4Xp5WfyYOHf8hXn/kvvvaqrxH2hWl94gt8t28vAYQy42FY0oxlzXvrBX7mWYd4POwIHWZgTZrWS+6hbzTJo7FWjo13I8aDMQLGw7gnxJDPMGJ6GR+OsRwvn5OzGTEB/nf6GZK+VSwv/wA9QwkGBnbSn/RzfKCM0ZSPymAf/kqhegks8fSyZLSX5EAd3/nFBu4yy0i5qlG9JIlKL2HPMFFPnEpvDxWePn7DmRxPL6c2KbwmeJxAWS2h8iq8qRF8Y308xVqOJyMsM0d5hf9ZvJXLCFatpLxmOSuqy62g610YVYDZE2keP36c7u5uqqqqKC8vJxAITKlkN+NqMRF52Bjzphk9eXqv0wi0GmPWuvZFsab/P+EdTvd4gMrKSvP5z3+eG264gVQqxeWXX87b3/52Nm7cyOjoKDfccAPXXXcdb3nLWxgcHOSmm27i+uuv56qrriIej/O+972PG2+8kaamJo4fP86HPvQhtmzZwutf/3qOHDnCP/7jP/IP//APvOY1r+G5555j69atfOhDH+LSSy/l4MGDfPKTn+RjH/sYDQ0N7N+/n3/553/iE1tbWL9+PU8//TSf/vSnufXWW7nwwgt58skn+exnP8unPvUpXrp6Bcd/9hV+8cvvce4rzmEskOTIUDffPZRg1YbLeHb0aZ5PxlmSSPO3hxp4Yvwsnj5rN929Z/LK6BbeSRvfTj3C6/tH6At6+VFZhN+EggCkRk6jeuylRMdO5y+/Oc6mqy7H732BPzz3DE90DnLt1W8lNPw9nu2O0T98jOplEUYZZF8qzJFxPxtWrudvfMN0/3GY7/3kae668w68I8f5r+//gkd+3sF/fOFzBLr/yNd/+At+vPfP/J+7rSaxL33pSzz11FN84QtfAOCuu+5i37593HnnnQDceeedHDp0iDvusMbr3nHHHXR3d7Nt2zYAbr/9dgYHB/m3f/s3AD796U+TTqf51Kc+BcBtt92G3+/n1ltvBeDWW2+lsrKSlharoNvS0sLy5cv56Ec/CsDHPvYxzjzzTD784Q8D8OEPf5hzzz2X97///QC8//3v5+KLL2bLli0AbNmyhVe+8pXceOONALz73e/m8ssv5/rrrwfg+uuv58orr+S6664D4LrrruNtb3sbmzdvJpVK8Y53vINNmzZx9dVXT+tv742vexXHe/r44Ec+yj+853pe33A+x3r6+Mitn+Wm933ghL+9yy65iONP7eFbD3yJxte9kuVLq/hNuIE77/p3Wq5exxrvEarigxz+8x9ZsfYMusPC86d/gju23cGdf1fJ+OAThIwhnDaETBpvoIrvnPkR/mvXf7H1ijDpof1cNThEhTE8Hgryk6rlLC37OK1fauXe/5GmrO/PVJkQx/oTrFi+nLFlF/CfY6/i67u+zhff7ic51sX5/UMMDQ3zyxUVPFt9Ft7+q7nvvvtY/84efuv3kBIhCSRFMCIEUxHGukd5bfkIZ6fghq5BSIyyY1U14ehLSB54JXt/+Rg9m48Qz6ry9iL4E0FkKElZ+RiXj4xwW1cP40ZoWHM6TVzEn556ObWpJ+lc93iuy0iGzxg+3BPnXf0D/MXn44aVyzhz6Cp+d/B8Lq16lL1n/Y4lyRRLU9ajIunlB2NN9Iwu4+J0N39dM0K5CXD80CEk4MW7NMJT0ev4xc+f4OW1g7xl3Ur84Uq+853vsmrlcs4771x+U34J3/rBd7hseZAr159HuLKWO+9q5eJ1DVx51ZvoNeX8z9v/F69dv4p1Lz2TsnANX/7yvWx4+WVc9ddvZXA8zdZPfZrXv+alnHVWlIoRPz/5yWOc/tIzWXH2mQTDG7jr7rtpuhyojDOQ7mF0KMG9n/hVZzqdviTn55D3U8rvVJX9opNmQCRqjInP8njGxsa48cYbMxeE0047DXfQTaVSpNPpaW07z3cGNGWfzxgD6SRLjv6MD5z/PBc8eye+w2PUDR3l3fXCUGoIk07zsv1fYMvaI7z02f/LssODnN2zn4tfbhVdPWO9fHP0Ib5zeTlwyDp5BXAeHOr/HWZsFeG+9TzTV88jy1Zy0Rk1vPDnizn8l2P8aniMR3k5VWc+wSOrhgEIjlaTOl5Pen+Yf7nm7Vx6zjJ+8P3vc++BR9l4yd9TVnY+3/72t/nTMz9i8/pmgsH309bWRtv327jvvv/E4/HwjW98g4ceeoh/vc+62H71j1+ld+wPpMpXkCpfwZDpIJlMkg5VM3rapXSbpxhJvPjZpdNpUqnUhM/Kve3UC+dLz7edSqUmFO+ne75kMnnCd5l9fvf5kslk3r8V9/km+1uZyt+e8UcwniFSqRQJCZKsWs3Y8AuMJk3Ovz0TKKc/ehG7YxWsu+avWHPJBpbu2wfHQaqvpKa+nqeffprbvn4b/3zVP3HRRRdx7He/I5VK0VP3Ic4//3z27t3Ljh07+MxnPsvZa+tY8XgH8hdh7aq/5czwCPv+/Ht++fPHeMtVb+TGJav47nHr9YcuuoXgqtXc94tOWltb+eIXv8iSJUsI/uhHmD6D71V3UVtVxb3f+x73fv9evvIf/0FDyMe3vvcjUqkUN7/0Xwh4Pfz0pz/lpz/9KZ+89ZOYyBLufbiThx56iI99/UHw+Pj8V7/KI48+wle/+lUwhq/s+78MjqX52rr/zcjwMXY/+m2e+8tB3v2uDxGqPo+77v02Tz75JK13bifdu4/7HvwKxJ/h+vPquGDVZSSf2E/8+BGay9fj9QZ56sn9DAyk+Ju3/i3+0FK+uus+Rujl0lefx8pVK9n5ncfoBuovPJ2m9a8l1fEgg8k0VzS8lt7R4/z5hQPsCxmGytIk5Gf4gaftx/85cpy3l43waDjMh6uXUvn8Bo7WXoypvZtPjn6XshGDeTX8TP7EWPrn0P8f8Gr4CfCLgYf4+p+O8N11Cb5T/jiffvLr/OnAZ+C8qxmq+jd2DA/CMPjeZtiT3sfPf30vR6UcNo3wiDwF4/Dr557jhjWG7eEo9/ZXcOzxNXDBJrpDt/J4MEVVKs2qlI/ly5evZBLFUnJpM8ZUu/bVYS2zXJ0dLKZ7PJziBv3xIY51dVK5fAMhT4CfffdKtldmtxhYfOJjSTLBkmSCz3cNUFVxFl+uiPAlOca7ltzNc13wTN8e9o0NMzq6ApOsxCTLMakINeEAL10R5vVrq7jsrIqcDZmH+sZ48Mlj7D7UwfLQMt567vm84ewqyk7S6KmUKhxjDAPJAbrGuuge66Z3rIt10YtZFl5J93gPv+/7PfXV9UR8ZfzuhV/xy6M/ZXB8AC9CxBehzB8mVHE6EX+YwYHD9A0e5g1mLWXDI/w++Ry/8h3j8tW3sTRSyVDfT3l24PcMJQYZGh9kODXMSDrBcOX5VAeqOD01wDKvl9dFL0aMcGD4CC+kRli27GrS6TRlqeeojS5j6dIznJH/k5ZciiG41AN73VVaufbN9Hg4NcHF13eQij/cz1PP/YDmpZV85OybkZELefKZX9KROEg6HcKkfaTTPgxCwD+CLzCI+AYwnn5Wjr+Nw91L+ctAF0bGMYlqygM+zl4S5tylIdbWhlhZEWBZudUrxr9A6m+VUsXLKbk7K2RWVFRMaHMp6t5ixphOEYln7a4B2gtx/FwaTg4Te+5HPBr7BnX9R3hf3xChildzfu8qPvNDQzp5iLLAas5Zcr7VndJv9UQRscYD9MVT9I+m6BtNMhTycXZtgNfW1XHOEiuYrKqcWsOaUkpNldMt2RiD1+tl6dKllJeXT/tme8EHF9sDIrLRNU6lCWh1Eu1qr3pXet7j50I6McJz+77Bb488xv5kL38oq+KFRByAcjEM+S7i5SN/R+9YFfWnlXNTQxnrTyvLORhQKaVOJfc4FyAzzqWsrGzGN7BFEVyMMc0istVuT6kDDmQNiNyIFUB2TfH4gpHEEN/6dQsPjcY47vUgYljh9RMYqaVi7FXE40s5NnAWT1ZUcfUrqrnyvCi1keJaC1spVZqMMSST1tRFwWAwM0K/EPONFUVwATDG7DhJ2o4c+6btDy/08832X3Dty1+C11uR85jHj/83Dz3bxs0v28b+owlSXfs53V9OMH4Jfxx4A/3pMpaV+zkzGuDSs0K8enUF604rw6NVWEqpBcAJKiJCVVUV0Wh08c4tdqrsfeoR2sf/hcivXsE1l9+W2f/7+O8Z7kqz5E8PExj6ISNVAd5+738znlhC2LONc5dX8aqXlPPB08pZUxsk4tceV0qphSN7ipdoNEo0Gl2QsyKX5G345te9jbav/E++y6/5q7EBwsEKvn/oh7Tu+yJXDQ5ze9dx/lvWcVngLVxw0bm8ZFk5604r0+67SqkFxxk75fweDAapqKiY3/VcpmB7wXKxgAT8Pq6MXsMXR3bx1V99hv6aM3is64e8amSEzSOns/9v7mFVbR3XzHdGlVJqEk4XYhHJrEQZCoXmrJSSy4xfyRizp5AZWUj+5lU38evvfJPvRJ6Crqcw8ddwzen11F7VeNIJH5VSaj64Z3jw+Xwz7kJcKNrmkoPX66Exei3xwXvxjP41W694Lysqco+iV0qp+eKu9gKrC3E0GiUcDs/7GDgNLpN45cv/B2c+u55VZ74En18/JqXUwuDu6SUilJeXZ6q9FtIaVHrVnISIcObaC+Y7G0opdUJPr8rKSiorKwkGg/NeQpmMBhellFqAskfNh0KhzKj56S7cNR80uCil1ALh7uUFVhtKWVkZ4XCYQKC42n1nHVxE5GpjzIP272uw1lNZ4+xTSik1OXejvDO4sby8fEFXeU3FbJY5rrR/rXH93g30sAiXFFZKqanK7uVVVlZGZWXlgujlVSizKbk0Ac3AGmATL47Y72WOZyBWSqli4p7G3hGJRKioqCASiRRFG8p0zWYQ5TeBb4rIFaU8oFIppabDGJMpmbhLIYFAgGAwmGlDKcWA4jbrNpdcgUVEVhtjnpntuZVSaqFzgolTMhER/H4/lZWVhEIhAoEAPp+vZKq7pmpGwUVEqoArjDEPisjHs5Ox2lzmfAlkpZQ6VXIFEYff76esrIxIJHLK5/BaqGb6CTRitbUAvBy4Pyu9b8Y5UkqpeZY9aNHh8/kIhUKZEonf78fv9y+6UslUzDS49NptLgD3u34HQERO+Xr1Sik1FU6jenYjuztAiAiBQIBwOEwoFMLv9+Pz+Uq+naSQZhpcLhGRHmPMb4BqV1dkx03AHbPK2QLnFJGNMXg8Hr1zUWoeuXthubedAOL+/3T+Z51ShxM4vF4vXq8Xn8+3KNtICm3S4CIilcaY/lxpxpjPicjtIrIDa416pyuysX+uoQSDi3sqBmMMXq8Xv99PIpHIBBoRwev16h/mPMmuC3fvz/ecfLLPl+t4977sO+DpPCfXc0ud+3NwB4VcP2Hy78OZyNF5OAHEqb5ygoaWQE6NfCWXW+xHTsaYmwFEZL0x5gl3mohcUZjsLRzOtAy1tbWEw2H8fn9mBlJjDOPj44yOjjI8PMzIyMiEUo2WbGbOfXFxP6YbQBzui0/2/sleP/s13T+d353v2F3V4jyyz++++GU/x7mByfW87Hzl+5vK9Vm4A91c/z1O5btwc38ezv+Mx+PJlCa8Xu+E/bke+j+2sOQLLhtF5BFjzI/znSA7sNj7Smrci/MPv2LFCsrLy09IFxGCwSDBYJCqqqoJwWZoaIiRkRHgxdJOrovbYpArULj3T3ZH6r7IOHegTnDPvhBlV384iu0zz+6ZlP0z192687zsY3KZLMi6093Hufe5f8/1GTsBIvs52e8v+3mqtOQLLkuAdvtL7wTagd1AR67qsnzVaMXMWdlt5cqVRCKRKT0nO9ik02nGxsYYHh5mcHCQRCKRuWPNLso7zz+VTnaXma8+e7JjJitZuOu1c92V5vp9pp9HMV+wsi/SpaaYvxs1NfmCSz1Wl+MerHaVJqAFMCIS48Vg02kPmMxbjVaMnLl/TjvtNEKh0IzP4/F4CIfDhMNhamtrSSaTJBIJEokEY2NjjI2NkUwmc3Z9dMu+W51Mrou/Y7KL/mR1/dnBz33Ryy41ZFdrZB+rbVFKLR6TBhdjzEHgHrv9JGaMeSOAiNQDl2AFmy8DVSISt59WUsHFGMOKFStmFVhycRoVw+HwCWlOgHHXvbvr4pPJZOaRXXJwLuTOhd5pE3KXiLLTtb5aKTUXTtoV2Wk/EZFrgAPGmE6sarJ77P1rsALNJ+Ywn6dcOp3G6/VOuSqsUJyLvVJKFbMpX8XsgZJ9InK1e1yLMeagMWYn8M3Jn118UqkUVVVVekevlFIzMK1BlHZV2UERuUJEqszEBcFKapp9j8dDRUXFfGdDKaWK0oxG6GdVlfUaYx61A0/JiEQi+P3++c6GUkoVpRlX7ttVYweAOhG5X0RWFyxX88zv91NdXT3f2VBKqaKVd/oXrC7IzmNt1jZYU73EgRjWqpQl0VustrZ2vrOglFJFLV+1WBxrrrA4cBArgDwB7LJ/j5VaVZhSSqnCyBdcYliN9HGgXQOJUkqpqcoXXHYZYz4H1uSUIvIeV1rMGPOo+2ARuTqr95hSSqlFKt8I/Ztdvz+BVSUGWAMnReQmIIpVddYHbAU0uCillJpxV+SD2CP0ITMljLaCK6WUAmbRFdnNnhJmZyHOpZRSqvgVbBIrdzWaUkqpxU1nSFRKKVVwGlyUUkoVnAYXpZRSBafBRSmlVMHNqCuymz0HWSOwAWvcSxx4HGtUf/9sz6+UUqr4zCq42FPuN2HNjhxzJa0FmkRkt47aV0qpxWfWJRdjzHsnSxORK2Z7fqWUUsVnrttcqub4/EoppRag2ZZcYiLyAFa1WLe9z5kGpo4SW/pYKaXU1Myq5GKMecIYsxnowFo4TIAeoMMYc232zMlKKaUWh1m3uQAYY75ZiPMopZQqDXPa5iIiV8/l+ZVSSi1Ms+2KfDuTN9oLcAm6xotSSi06s60WO2D/jOU9Siml1KIy2+DyANBgjNmTK1FEZnn6zHm2YgWwGgBjzKRrx4hIFNhib24Aduc7XimlVOHNKrgYY/qAnIHFTp80bapEZDvwuDFml7MtIhud7RxuMca0uJ5/QETyBiSllFKFVQwTV27JCiT3A825DrRLLXVZu1uBlhOPVkopNVemFVxE5O6TpN8uIg+LyMdnl63M+epz7I5jTZQ5mUYRcQeYOCcGHKWUUnNouiUXp2pqnT0bcoaIbAO6jDFvAp4oUDfkGqxBmW7Z2xnGmLgxptoY4+5g0AS053uRw4cPIyKZx2233TbjDCullJpBm4uI7Ld/rRaRbcaYO+ztjcaYc8Bqa7FnTJ6taJ58RI0x8ZPkNYpVysk7geaqVas4fPjwDLKnlFIql+mWXOqBS4wxZxtjaoE+Vwkmu2tYfLaZs89Rk7Uvezufe4BNxpjOAuRFKaXUFE235BKze4g5HgAagEeB3qxjTa4TiMhG4NqTvE6PMaYZqwosmpUWBasKLN8J7O7LrcaYvFViSimlCm+6waVGRN6ANVFlDVavrW0iUgVUZx1bhxV0JrB7fk3WjTj72E4RiWfngZO0odgBrNMJLCLSqEFGKaVOnWlVixlj7gE2A88AbVgDG5uA24FmEfm4iKwWkfdQuFH7D9jBwtGEayp/Ealzp4tII1YA6hCRqN1zLFevM6WUUnNEjMlZezWzk4msxyrNPFLI5Y3tKq5O7C7F7gGRdlqTMabJbsDPrp4D2GWM2TTZ+RsaGkxHR0ehsquUUouCiOw1xjTkTJtJcLFLJpdglRB2G2O+PLsszi8NLkopNX35gsu0R+iLyCNYjfgxrLaXBhF5PHvci1JKqcVrWg36dollU1aPMfdkkXfkep5SSqnFZboll97swAKZbsEHC5IjpZRSRW+6wSVfA03hegYopZQqatMNLrUisjp7p4isA9YWIkNKKaWK37TaXIwx94jIAyKyhhfHsdRhjdw/2ah7pZRSi8S0J640xmy2x7M0YE3Fcrsx5olCZ0wppVTxmtFKlHYwmRBQROTjrhmSlVJKLWKTBhcReXga5xGsQZUaXJRSSuUtuQjW8sDxKZxHsOYXU0oppfIGl5bptKWIiK5Tr5RSCsjTFXm6jfTGGB1EqZRSCpjB3GJKKaXUyWhwUUopVXAaXJRSShWcBhellFIFp8FFKaVUwWlwUUopVXAaXJRSShWcBhellFIFp8FFKaVUwWlwUUopVXAaXJRSShWcBhellFIFp8FFKaVUwWlwUUopVXAaXJRSShWcBhellFIFp8FFKaVUwWlwUUopVXAaXJRSShWcBhellFIFp8FFKaVUwWlwUUopVXAaXJRSShWcBhfg8OHD852FOXfbbbfNdxbmXKm/x1J/f6DvsZSIMWa+8zDvRMSU+ucgIuh7LG6l/v5A32OxEZG9xpiGXGlaclFKKVVwGlyUUkoVnG++M7BADIvIH+Y7E3NslYiUeuNSqb/HUn9/oO+x2Jw1WYK2uSillCo4rRZTSilVcBpclFJKFZwGF6WUUgWnwUUppVTBLdreYiISBbbYmxuA3caYnfOXo7knIq3GmOb5zkeh2d/lLcABe1eHMaZz/nJUWCKyFYjbm1FjzI55zM6suf73ao0xLTnStwIxoAagGP8v873HxXLtWbTBBbjF/aWLyAF75GzJfckAIrIdqJvvfBSa/Y/aZoxpsre3YAWaTfOZr0IRka3uYCIi9dn7iomINAJRYO0k6duBx40xu5xtEdnobBeDk71HFsm1Z1FWi9kXpOwLbStwwl1UKRCR+vnOwxy6B+u7czxAaX2P17o37BLZhnnKy6wZY9rtQBGf5JAtWYHkfqCoStv53uNiuvYsyuBiaxQR95ccpwTv7G0NwO75zsQc2Qi0i0idiNQbY+LGmNh8Z6qAekSkzdmwS2b3z2N+5swkN0FxoPEUZ2WuLYprz6IMLvYFqDrrItQEtM9XnuaKiGzEupsvOa6LUYNrX5t9d1gqmrEuRr12W0RPMVURTVMN0JO1L3u7qC2ma8+iDC7Z7ItRIyVWNLXfV9wYE5/nrMyVzN2eMSZmVxndj1VVVhLsi9A2oAPYThFXiU1BdLKEErthyCjVaw9ocHHcA2wqpR5Gts3GmJK7I3KJ2z87XPtiWFVlJUFEWoF2u8NCE7DFXU1WYuLYPcRcsrdLTaleexZ1bzEg0+2xtdQuwnaVUUm9pxxiYFU1uPbFwbojLPYSm/0dxp0LjzGmXUTWAAfnN2dzpocTSy9ROOE7Lgmleu1xLOrgYrdHdDpfrog0ltAXXYNVV+9sbwDq7D/oXaXQ6G2MiYlIPCuQRCmdqsAaoNu9wxgTF5FS+RudwBjTKSLxrN01lOBNUolfe4BFXC1m90WvATpEJGr33iiZLrt2d8gdzgOrt1jc3i76wOKyDdjs2r7W3lf07ItNk3ufXUdfSt9ftgfsC6+jiYldzYteqV97HItyyn37H7Q3R9IuY0xJDL5zs7uvbsLqVbUN2Fkid/ZApnoho1gHGOZiX3iaeXH2gaIebGdX9TXy4tgVp02p03XMVqATu8NGsb3ffO9xMV17FmVwUUopNbcWbbWYUkqpuaPBRSmlVMFpcFFKKVVwGlyUUkoVnAYXpZRSBafBRSk1gT32otV+nNLZekVkq/26JTOFz2K1qEfoK6UmFc+zSmQt1swBcexZmu1g0D7b8VPGmB2ulUVLdfbnRUGDi1LTZK+WGC3EktG55kAr5PkLxVnxE2suLPfKmFF7kG4LcMk8ZU8tQBpclJq+Qi7WtRnIHoG+EBcD2wO0ZM9/Zc919gAlNkWLmj0NLkpNU4GnR28iK7gstOnX7ZJJz2QTK9oBRquw1ATaoK/UPHGqv+Y7H1PQjFUlls9CLG2peaQlF1U0nDtoe7PGPaHhFNKcmYSj2DPSYq3sGMeazLPGTttgjGmxZ64Fa7bamLO0sN17qhXAXsDrpK8/yXvZaL+eswwCWCWYGvf57UkQp51P1+s4k0BGp5KvSdRzkpmY3a/raoOJYy2EFbP377bzf5N96D1Y30Mb1iSVTaU2eeOiZozRhz4W/APrArvFtV3nbGNfnFxpUWC3/fuWHGnb7d8bgb1Z6W1Oumtfb9Z2vXP+k+XtJO9pwnnynH+m+WwD6rPyuXEK+XJ/RnWAcZ9nit/XFqzGf/e+jVgdFU54X/Zrbs1+fX0U70OrxdSCZ/dU2mIm3nVvBNbad/Z1xrVGjbF6X8XsO2h4cepzJ82pwunButi578pz3aH3ZK3hHp9K3qbw1iYTz9qedj6dNULMxPab+3F9FlNhZrj2j/15bM6xP+7azLwvY0zclNBSCUqrxVRxaCTrYupciLKqvNwOAJcYY5pFpE1EDNaKhm1ZgSCe47ndOfZNO28FFs+xL18+G4G4q9oM7Gq4Gbx2DGstoAkdDZzAihVIt2BV6213BaR2EdlorHEwUXK/h1Je+GxR05KLKnbRPGk19jiSTUA1VrXQJhFxd5vtyfG8eOGyN3UnGQ0/3XxGsdpg2l2PXcaYmZSodmEtNjeBq7ThDLbcnlXS2caLJaXJlvGNzyA/qghocFHFILMqoZt9N9yeKw3rbno31khv50LYbqxG+EJOaZIvb9NVyKVuc+ZrJow1Ur/OroLMZ0IAtKvk6k71FDJqYdDgohY8+274AVcbiqPRvoDF3Bc++8Le4FR/5Xie++66JsdLRguRtyk8PcbEABDPc+y08mmXEnqyA0KOfE5VE7B9kjm/NjN59Var/Zhs7E50hvlRC5y2uaiiYLedbHW1sUSN3f3VGLPJTnMu1HXAFfbv3VjBx7nYR7EukvVYpZo6EdlqrDmtNmI1xsdFpNMY02535a2zn7PdPsd2oEFEthhjdubL20neU9yepNEZpLjLfg+Z82N11Z12Pu1G8iY7Xw3YpQozw/Xo7SDqnG+787liBb12O5+57MTqYjwh+Lg+/wY77zvNLOclUwuLGKvrn1JKAZmS3y0mx8SVMzzfxqkE27l6fTU/tFpMKTVnRKRxOoFFlQ4NLkqpgrKr+py2nuh85kXNHw0uSqlCa8PqBj6t6jBVWrRBXymVS52ItGFNsz+tgY6TjGeZErtxfwM6EWbR0wZ9pZRSBafVYkoppQpOg4tSSqmC0+CilFKq4DS4KKWUKjgNLkoppQru/wEHfW5P/8a0xwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1)\n",
    "ylim = ax.set_ylim(-0.35, 0.35)\n",
    "xlim = ax.set_xlim(2, 13.8)\n",
    "__=ax.plot(np.linspace(0, 20, 5000), np.zeros(5000), ':', color='k')\n",
    "\n",
    "__=ax.plot(tng_t, mu_200m, label=r'${\\rm 200m}$')\n",
    "__=ax.plot(tng_t, mu_500c, '--', label=r'${\\rm 500c}$')\n",
    "__=ax.plot(tng_t, mu_mvir, '-.', label=r'${\\rm vir}$')\n",
    "leg = ax.legend()\n",
    "\n",
    "__=ax.fill_between(tng_t, -std_mvir, std_mvir, color='lightgray', alpha=0.7)\n",
    "\n",
    "xlabel = ax.set_xlabel(r'${\\rm cosmic\\ time\\ [Gyr]}$')\n",
    "ylabel = ax.set_ylabel(r'$\\log_{10}M_{\\rm fit} - \\log_{10}M_{\\rm sim}$')\n",
    "fig.savefig('tng_fit_residuals.png', bbox_extra_artists=[xlabel, ylabel], bbox_inches='tight', dpi=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "skilled-style",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
