{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quickstart Tutorial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`galtab` is a general approach for calculating the expectation value of\n",
    "counts-in-cells statistics for a given halo catalog and HOD model. It pretabulates\n",
    "placeholder galaxies inside each halo to yield rapid, deterministic results,\n",
    "which is ideal for MCMC likelihood evaluations.\n",
    "\n",
    "This [tutorial](https://github.com/AlanPearl/galtab/blob/main/docs/source/notebooks/intro.ipynb)\n",
    "will demonstrate some basic Counts-in-Cylinders (CiC) calculations\n",
    "using the intended `galtab` workflow.\n",
    "\n",
    "To cite `galtab`, learn more implementation details, and explore an example science\n",
    "use case, check out https://arxiv.org/abs/2309.08675."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prerequisites\n",
    "\n",
    "All of the following are `pip` installable\n",
    "\n",
    "- `galtab`\n",
    "    - `numpy`\n",
    "    - `jax`\n",
    "    - `astropy`\n",
    "    - `halotools`\n",
    "- `matplotlib`\n",
    "- `jupyterlab`\n",
    "\n",
    "After installing the above *and downloading the bolplanck z=0 halotools catalog*,\n",
    "you should be able to run the following cell. In this cell:\n",
    "\n",
    "- set our cosmology and CiC parameters\n",
    "- choose an HOD model\n",
    "- load the simulation data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><i>Table length=5</i>\n",
       "<table id=\"table139826825803952\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>halo_vmax_firstacc</th><th>halo_dmvir_dt_tdyn</th><th>halo_macc</th><th>halo_scale_factor</th><th>halo_vmax_mpeak</th><th>halo_m_pe_behroozi</th><th>halo_delta_vmax_behroozi17</th><th>halo_xoff</th><th>halo_spin</th><th>halo_tidal_force</th><th>halo_scale_factor_firstacc</th><th>halo_c_to_a</th><th>halo_mvir_firstacc</th><th>halo_scale_factor_last_mm</th><th>halo_tidal_id</th><th>halo_scale_factor_mpeak</th><th>halo_pid</th><th>halo_m500c</th><th>halo_id</th><th>halo_halfmass_scale_factor</th><th>halo_upid</th><th>halo_t_by_u</th><th>halo_rvir</th><th>halo_vpeak</th><th>halo_dmvir_dt_100myr</th><th>halo_mpeak</th><th>halo_m_pe_diemer</th><th>halo_jx</th><th>halo_jy</th><th>halo_jz</th><th>halo_m2500c</th><th>halo_mvir</th><th>halo_voff</th><th>halo_axisA_z</th><th>halo_axisA_x</th><th>halo_axisA_y</th><th>halo_y</th><th>halo_b_to_a</th><th>halo_x</th><th>halo_z</th><th>halo_m200b</th><th>halo_vacc</th><th>halo_scale_factor_lastacc</th><th>halo_vmax</th><th>halo_m200c</th><th>halo_vx</th><th>halo_vy</th><th>halo_vz</th><th>halo_dmvir_dt_inst</th><th>halo_tidal_force_tdyn</th><th>halo_rs</th><th>halo_nfw_conc</th><th>halo_hostid</th><th>halo_mvir_host_halo</th></tr></thead>\n",
       "<thead><tr><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>int64</th><th>float32</th><th>int64</th><th>float32</th><th>int64</th><th>float32</th><th>int64</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>int64</th><th>float32</th></tr></thead>\n",
       "<tr><td>1001.57</td><td>12810.0</td><td>200800000000000.0</td><td>1.00231</td><td>1001.57</td><td>202700000000000.0</td><td>-0.0026</td><td>0.0257357</td><td>0.02391</td><td>0.11954</td><td>1.00231</td><td>0.47559</td><td>200800000000000.0</td><td>0.28343</td><td>2812606193</td><td>1.002</td><td>-1</td><td>116580000000000.0</td><td>2811042639</td><td>0.41506</td><td>-1</td><td>0.593</td><td>1.190447</td><td>1091.38</td><td>17390.0</td><td>200800000000000.0</td><td>111800000000000.0</td><td>2536000000000000.0</td><td>-474400000000000.0</td><td>-6566000000000000.0</td><td>65777000000000.0</td><td>200800000000000.0</td><td>21.34</td><td>19.2231</td><td>59.7891</td><td>-18.8001</td><td>43.14082</td><td>0.63663</td><td>36.17984</td><td>17.96339</td><td>223780000000000.0</td><td>1001.57</td><td>1.00231</td><td>1001.57</td><td>158240000000000.0</td><td>16.1</td><td>8.51</td><td>-78.88</td><td>17390.0</td><td>0.12244</td><td>0.137953</td><td>8.629367</td><td>2811042639</td><td>200800000000000.0</td></tr>\n",
       "<tr><td>895.2</td><td>13760.0</td><td>179600000000000.0</td><td>1.00231</td><td>895.2</td><td>181000000000000.0</td><td>-0.01065</td><td>0.041987</td><td>0.06297</td><td>0.50587</td><td>1.00231</td><td>0.56181</td><td>179600000000000.0</td><td>0.29862</td><td>2811077105</td><td>1.002</td><td>-1</td><td>100360000000000.0</td><td>2811055606</td><td>0.50618</td><td>-1</td><td>0.627</td><td>1.146849</td><td>969.05</td><td>7324.0</td><td>179600000000000.0</td><td>128700000000000.0</td><td>1.074e+16</td><td>4931000000000000.0</td><td>-1.185e+16</td><td>47026000000000.0</td><td>179600000000000.0</td><td>41.91</td><td>41.2062</td><td>34.6803</td><td>17.8882</td><td>49.54417</td><td>0.8397</td><td>45.36644</td><td>40.01593</td><td>204460000000000.0</td><td>895.2</td><td>1.00231</td><td>895.2</td><td>142290000000000.0</td><td>2.46</td><td>264.77</td><td>-128.08</td><td>7324.0</td><td>0.491</td><td>0.185805</td><td>6.172326</td><td>2811055606</td><td>179600000000000.0</td></tr>\n",
       "<tr><td>853.83</td><td>4666.0</td><td>129800000000000.0</td><td>1.00231</td><td>853.83</td><td>149500000000000.0</td><td>0.00531</td><td>0.026461901</td><td>0.03607</td><td>0.07568</td><td>1.00231</td><td>0.66381</td><td>129800000000000.0</td><td>0.49606</td><td>2810630242</td><td>1.002</td><td>-1</td><td>87766000000000.0</td><td>2809250167</td><td>0.491</td><td>-1</td><td>0.5774</td><td>1.029343</td><td>926.37</td><td>2747.0</td><td>129800000000000.0</td><td>80320000000000.0</td><td>2133000000000000.0</td><td>-3236000000000000.0</td><td>-3111000000000000.0</td><td>39496000000000.0</td><td>129800000000000.0</td><td>23.35</td><td>-17.5268</td><td>38.9596</td><td>24.3626</td><td>13.88261</td><td>0.76149</td><td>22.02318</td><td>9.80153</td><td>141210000000000.0</td><td>853.83</td><td>1.00231</td><td>853.83</td><td>112010000000000.0</td><td>18.49</td><td>124.89</td><td>-35.19</td><td>2747.0</td><td>0.10074</td><td>0.119293995</td><td>8.628624</td><td>2809250167</td><td>129800000000000.0</td></tr>\n",
       "<tr><td>777.64</td><td>4401.0</td><td>103000000000000.0</td><td>1.00231</td><td>777.64</td><td>104800000000000.0</td><td>0.00498</td><td>0.0516998</td><td>0.05031</td><td>0.09677</td><td>1.00231</td><td>0.47302</td><td>103000000000000.0</td><td>0.38469</td><td>2820592816</td><td>1.002</td><td>-1</td><td>57781000000000.0</td><td>2809483946</td><td>0.65806</td><td>-1</td><td>0.6152</td><td>0.952978</td><td>831.17</td><td>2747.0</td><td>103000000000000.0</td><td>64200000000000.0</td><td>1713000000000000.0</td><td>-1488000000000000.0</td><td>4582000000000000.0</td><td>30529000000000.0</td><td>103000000000000.0</td><td>98.45</td><td>24.7744</td><td>-10.3568</td><td>38.9949</td><td>36.67881</td><td>0.7881</td><td>12.29788</td><td>34.18085</td><td>115110000000000.0</td><td>777.64</td><td>1.00231</td><td>777.64</td><td>82069000000000.0</td><td>-281.37</td><td>-115.39</td><td>-391.28</td><td>2747.0</td><td>0.10259</td><td>0.132334</td><td>7.201309</td><td>2809483946</td><td>103000000000000.0</td></tr>\n",
       "<tr><td>748.56</td><td>11480.0</td><td>99470000000000.0</td><td>1.00231</td><td>748.56</td><td>107600000000000.0</td><td>0.05989</td><td>0.0779697</td><td>0.0348</td><td>0.12465</td><td>1.00231</td><td>0.47409</td><td>99470000000000.0</td><td>0.63275</td><td>2809483946</td><td>1.002</td><td>-1</td><td>59100000000000.0</td><td>2809272603</td><td>0.63781</td><td>-1</td><td>0.67</td><td>0.941893</td><td>748.56</td><td>5218.0</td><td>99470000000000.0</td><td>70970000000000.0</td><td>1207000000000000.0</td><td>-2126000000000000.0</td><td>-2677000000000000.0</td><td>26267000000000.0</td><td>99470000000000.0</td><td>118.79</td><td>29.2183</td><td>52.7796</td><td>6.18836</td><td>26.12877</td><td>0.66155</td><td>10.66037</td><td>22.5009</td><td>108110000000000.0</td><td>748.56</td><td>1.00231</td><td>748.56</td><td>84337000000000.0</td><td>-43.87</td><td>292.95</td><td>-171.47</td><td>5218.0</td><td>0.1579</td><td>0.14077</td><td>6.691006</td><td>2809272603</td><td>99470000000000.0</td></tr>\n",
       "</table></div>"
      ],
      "text/plain": [
       "<Table length=5>\n",
       "halo_vmax_firstacc halo_dmvir_dt_tdyn ... halo_hostid halo_mvir_host_halo\n",
       "     float32            float32       ...    int64          float32      \n",
       "------------------ ------------------ ... ----------- -------------------\n",
       "           1001.57            12810.0 ...  2811042639   200800000000000.0\n",
       "             895.2            13760.0 ...  2811055606   179600000000000.0\n",
       "            853.83             4666.0 ...  2809250167   129800000000000.0\n",
       "            777.64             4401.0 ...  2809483946   103000000000000.0\n",
       "            748.56            11480.0 ...  2809272603    99470000000000.0"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import jax\n",
    "from jax import numpy as jnp\n",
    "\n",
    "from astropy import cosmology\n",
    "import halotools.empirical_models as htem\n",
    "import halotools.sim_manager as htsm\n",
    "import halotools.mock_observables as htmo\n",
    "\n",
    "import galtab\n",
    "\n",
    "# Set our CiC parameters (all lengths are in Mpc/h)\n",
    "proj_search_radius = 2.0\n",
    "cylinder_half_length = 10.0\n",
    "cic_edges = np.arange(-0.5, 16)\n",
    "\n",
    "# Set our cosmology and HOD model\n",
    "cosmo = cosmology.Planck13\n",
    "hod = htem.PrebuiltHodModelFactory(\"zheng07\", threshold=-21)\n",
    "\n",
    "# Load Bolshoi-Planck simulation halos at z=0\n",
    "halocat = htsm.CachedHaloCatalog(simname=\"bolplanck\", redshift=0)\n",
    "halocat.halo_table[:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate CiC the standard way with `halotools`\n",
    "\n",
    "- Populate the halocat with galaxies probabilistically from the HOD model\n",
    "- Compute the number of neighbors within a cylinder around each neighbor\n",
    "- Tally up a histogram of the neighbor counts for a given set of CiC bins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.3675878 , 0.25521013, 0.15238918, 0.08963731, 0.05152562,\n",
       "       0.03218192, 0.0193437 , 0.01162925, 0.00420265, 0.00529649,\n",
       "       0.00362694, 0.00270581, 0.00166955, 0.00149683, 0.00080599,\n",
       "       0.00069085])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Choose your HOD parameters (in this case, we will keep them the same)\n",
    "hod.param_dict.update({})\n",
    "\n",
    "# Populated model galaxies and get their Cartesian coordinates\n",
    "hod.populate_mock(halocat)\n",
    "galaxies = hod.mock.galaxy_table\n",
    "xyz = htmo.return_xyz_formatted_array(\n",
    "    galaxies[\"x\"], galaxies[\"y\"], galaxies[\"z\"], velocity=galaxies[\"vz\"],\n",
    "    velocity_distortion_dimension=\"z\", period=halocat.Lbox, cosmology=cosmo\n",
    ")\n",
    "\n",
    "# Compute CiC (self-counting subtracted by the `-1`)\n",
    "cic_counts = htmo.counts_in_cylinders(\n",
    "    xyz, xyz, proj_search_radius, cylinder_half_length) - 1\n",
    "cic_halotools = np.histogram(cic_counts, bins=cic_edges, density=True)[0]\n",
    "cic_halotools"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now let's do it the `galtab` way"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "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"
     ]
    },
    {
     "data": {
      "text/plain": [
       "DeviceArray([0.35819536, 0.2572313 , 0.15485156, 0.09022042, 0.05297909,\n",
       "             0.03188442, 0.01996169, 0.0124573 , 0.00806456, 0.00497749,\n",
       "             0.00318336, 0.00215645, 0.0015641 , 0.00095675, 0.00078415,\n",
       "             0.00053201], dtype=float32)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Give the Tabulator the halo catalog and a fiducial HOD model\n",
    "gtab = galtab.GalaxyTabulator(halocat, hod)\n",
    "\n",
    "# Prepare the CICTabulator to make predictions\n",
    "cictab = galtab.CICTabulator(gtab, proj_search_radius, cylinder_half_length,\n",
    "                            bin_edges=cic_edges)\n",
    "\n",
    "# Choose your HOD parameters (in this case, we will keep them the same)\n",
    "hod.param_dict.update({})\n",
    "\n",
    "# Predict CiC for this model\n",
    "cic_galtab = cictab.predict(hod)\n",
    "cic_galtab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot the `galtab` vs. `halotools` comparison\n",
    "\n",
    "- `galtab` predicts the CiC expectation value (smooth + deterministic)\n",
    "- `halotools` draws a CiC realization (noisy + stochastic)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkUAAAGzCAYAAAAhXWNYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/av/WaAAAACXBIWXMAAA9hAAAPYQGoP6dpAABcsElEQVR4nO3dd3hUVR7G8e/MpIckEHogoYlAKKGFAAEBjWKjKlgQAQVFUEFcLGsBKzaQFkEsVEFUBBSxANIEpCQERHqPlISWCmkzs39kDY4TICGZTMr7eZ48a86599zfoJu8nHvuuQar1WpFREREpIwzOrsAERERkeJAoUhEREQEhSIRERERQKFIREREBFAoEhEREQEUikREREQAhSIRERERQKFIREREBAAXZxdQklgsFk6ePImPjw8Gg8HZ5YiIiEgeWK1WkpOTCQgIwGi88nyQQlE+nDx5ksDAQGeXISIiItchNjaWmjVrXrFfoSgPIiMjiYyMJCsrC8j+Q/X19XVyVSIiIpIXSUlJBAYG4uPjc9XjDHr3Wd4lJSXh5+dHYmKiQpGIiEgJkdff31poLSIiIoJCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCUZ5ERkYSHBxMaGios0sREREpVGvWrMFgMJCQkFAsxnEmhaI8GD58OLt372br1q3OLkVERMShZs2aRfny5Z1dhlMoFBUnF8+DXkUnIiLiFApFxcD4X/bxxnc7MM+8C+b3hQvHnF2SiEiZZbFYOZeS7tQviyXvf0FOTk6mX79+eHt7U716dT788EM6d+7MyJEjAZg7dy6tW7fGx8eHatWq8eCDDxIfH5/rWGvWrGHQoEEkJiZiMBgwGAyMHTs2X+Ns2LCBZs2a4eHhQdu2bdm1a1e+/x04i4uzCyjrjp+7yMdrDzOQpZhcd8OZ3VgjwzB0eRHaDgOTq7NLFBEpUy5czKDVmyudWkPUyxFULOeep2NHjRrFhg0b+O6776hatSqvvvoq0dHRNG/eHIDMzEzeeOMNGjRoQHx8PKNGjWLgwIEsX77cbqz27dszceJEXn31Vfbt2wdAuXLl8jXO6NGjmTRpEtWqVeO///0v3bp1Y//+/bi6Fv/fZwpFTvbOT3uoYonjGbdFOW2GrEuw4lWsO77E0G0SBLZxYoUiIlJcJScnM3v2bObPn88tt9wCwMyZMwkICMg55pFHHsn557p16zJ58mRCQ0NJSUnJCTx/c3Nzw8/PD4PBQLVq1Wz68jrOmDFjuPXWWwGYPXs2NWvWZPHixfTt27fwPriD6PaZE51MuMTKPfEEcJZkvOz6DfG7sX52GywbBZcSir5AEREp1g4fPkxmZiZt2lz+y7Ofnx8NGjTI+T4qKopu3boRFBSEj48PnTp1AuD48eP5ulZex2nXrl3OP/v7+9OgQQP27NmT78/mDApFThRQ3pNVozpRucnN3JL+AXOybsViNdgcY8AK2z7DPDUUdi3SQmwREcmz1NRUunbtiq+vL1988QVbt25l8eLFAGRkZBT5OMWdbp85WaC/F5H9WrL1aG3eWFadb090ZJzrpzQy2iZvU2o8fPMI5uh5mO6eAP51nFSxiEjpVsHLjaiXI5xeQ17UrVsXV1dXtm7dSlBQEACJiYns37+fm266ib1793Lu3DneeecdAgMDAdi2bdtVx3Rzc8NsNtu05Wec33//PaeWCxcusH//fho1apSnz+NsCkXFRGhtf5YMC2dJTG0G/9iQOy8u4RmXRXgZ0m2OMx3+FfPUMIxdXsDQ/iktxBYRKWRGoyHPi5ydzcfHhwEDBjB69Gj8/f2pUqUKY8aMwWg0YjAYCAoKws3NjSlTpjB06FB27drFG2+8cdUxa9euTUpKCqtWrSIkJAQvL698jfP6669TsWJFqlatyksvvUSlSpXo2bOnAz594dPts2LEaDTQu2VNVoy+Be8uz9DNMp5V5hZ2x5ks6RhWvcalqeEQqw0lRUTKsgkTJtCuXTvuvvtuIiIiCA8Pp1GjRnh4eFC5cmVmzZrF119/TXBwMO+88w4ffPDBVcdr3749Q4cO5b777qNy5cq89957+RrnnXfeYcSIEbRq1YrTp0/z/fff4+aWt5kvZzNYrVqkkldJSUn4+fmRmJiIr6+vw693OjGN93/aS+qOxYx1nU01wwW7Y+YEvEJE3+EElPd0eD0iIlL8paamUqNGDcaPH8+jjz7q7HKKhbz+/tZMUTFWzc+D8fc1Z/iwUTxf9VNmZnW1WYi9ztyUVw835Obxa5jwyz5S07OcWK2IiDjD9u3bWbBgAYcOHSI6Opp+/foB0KNHDydXVvIoFJUATWv6MeuJW6h23ySGerzDn5ZapFldeTnrEcBAWqaFyb8epMsHa/h6W2y+dkIVEZGS74MPPiAkJISIiAhSU1NZv349lSpVcnZZJY5un+VDUd8+y016lpnZ6w+yes1KNqXXzvWYJjV8+bDJUep3vA9cSsZ9XBEREUfR7bNCFBkZSXBwMKGhoc4uBXcXE491acCU0YPpFxaE0WB/TJVTa6i/9klOvdua0zt/LfoiRURESiDNFOVDcZgp+rd9p5N584fdrD9wFgBP0ljh/hw1DWdzjtleuTv1HhyPb4UqzipTRETEaTRTVEY0qObDnEfaMHNgKPUqezPC5VubQATQ4sx3ZE1qxW/fRpKVZb7CSCIiImWbQlEpYDAY6NKwCj+NvImaHfqxG/vdrv1JosPO/7JzXBe2bNPeRiIiIv+mUFSKuJqM3H37ndT4zyZ+DhxJitXD7piW5h2EfH8H3374NAdPns1lFBERkbJJoagU8ivnSddHX+P8oN+I9upg1+9uyKR34mz4uCOfzpvH+dTS8zI/ERGR66VQVIoF1a5Py+d+YE/nj4k3VLbrv8FwksEHh7Pm/fvYtDfWCRWKiEhBdO7cmZEjR173+bNmzaJ8+fKFVk9BjB07lubNmzu1BoWiMqBR5/up+Px29tUZgDmXf+W9+RXP+T34fuMOJ1QnIiIlicFgYMmSJc4uwyEUisoIk4cPDQZMJmPQSk57N7Lrr2/4i+nfr+fDFfvRLg0iIlIWKRSVMZ61WlHt2Q0kdH6LSwYvALKsRoZnPs2f1jpMWnWAZ7/aQUaWxcmViogUA6lnr/8r89JVxj2X+znXwWKx8Nxzz+Hv70+1atUYO3ZsTt+ECRNo2rQp3t7eBAYGMmzYMFJSUq463rRp06hXrx5ubm40aNCAuXPn5vTVrl0bgF69emEwGHK+v9Z5AMePH6dHjx6UK1cOX19f+vbtS1xc3BXrWLNmDW3atMHb25vy5csTHh7OsWPH8v4Hcx1cHDq6FE9GE+U7P4mlUQQJn/VgXGo31lha5HR/u/0EpxLTmP5QK/y8XJ1YqIiIk71f7/rPvfMDaDMk977IULh4zr59bGK+LzN79mxGjRrF5s2b2bRpEwMHDiQ8PJxbb70Vo9HI5MmTqVOnDocPH2bYsGE899xzfPTRR7mOtXjxYkaMGMHEiROJiIhg2bJlDBo0iJo1a9KlSxe2bt1KlSpVmDlzJrfffjsmkylP51kslpxAtHbtWrKyshg+fDj33Xcfa9assasjKyuLnj17MmTIEBYsWEBGRgZbtmzBYMjlNQ6FSDta50Nx3NG6wNJTmBN9lrHf/cm/3yN7Q5VyzBwYSqC/l3NqExFxtrF+13/u1ULRe3ULJRR17twZs9nM+vXrc9ratGnDzTffzDvvvGN3/DfffMPQoUM5ezZ7VmrWrFmMHDmShIQEAMLDw2ncuDEzZszIOadv376kpqbyww8/ANlrihYvXkzPnj1zjrnWeStWrOCOO+7gyJEjBAYGArB7924aN27Mli1bCA0NZezYsSxZsoSYmBjOnz9PxYoVWbNmDZ06dcrXn0lutKO15I17OR5uV5tPHm6Np6vJputgfApvT53OztjzTipORESupVmzZjbfV69enfj4eABWrlzJLbfcQo0aNfDx8aF///6cO3eOixcv5jrWnj17CA8Pt2kLDw9nz549V63hWuft2bOHwMDAnEAEEBwcTPny5XMd29/fn4EDB9K1a1e6devGpEmTOHXq1FVrKAwKRQLALY2q8tXj7ajs457Tdq9pLdMsr3H8k4dY+Yce2RcRKY5cXW2XORgMBiwWC0ePHuXuu++mWbNmLFq0iKioKCIjIwHIyCj++9PNnDmTTZs20b59exYuXMiNN97I77//7tBrak2R5Gha04/Fw9rzyKytVD2zkXEunwJwt3EDm77uy/xz03mwc4iTqxQRKUKjD13/uW7eV+4bvhVw7OqVqKgoLBYL48ePx2jMngP56quvrnpOo0aN2LBhAwMGDMhp27BhA8HBwTnfu7q6Yjab83Veo0aNiI2NJTY21ub2WUJCgs3Y/9aiRQtatGjBiy++SLt27Zg/fz5t27bN459A/ikUiY2aFbz4elBTTJMexNV6+T/6dsbd+P/6IBPPTeWpXp0xGR272E1EpFjwruSgcSs6Ztx/uOGGG8jMzGTKlCl069aNDRs2MH369KueM3r0aPr27UuLFi2IiIjg+++/59tvv2XlypU5x9SuXZtVq1YRHh6Ou7s7FSpUuOZ5ERERNG3alH79+jFx4kSysrIYNmwYnTp1onXr1nZ1HDlyhBkzZtC9e3cCAgLYt28fBw4c4OGHHy7cP6R/0e0zseNX3h/3+2eRZvS0aW9g/Iv7dw7irc++4lKG+Qpni4hIcRASEsKECRN49913adKkCV988QXjxo276jk9e/Zk0qRJfPDBBzRu3JiPP/6YmTNn0rlz55xjxo8fz4oVKwgMDKRFixZ5Os9gMLB06VIqVKjATTfdREREBHXr1mXhwoW51uHl5cXevXu55557uPHGG3nssccYPnw4jz/+eKH82VyJnj7Lh1L59NlVWE9u5+LMe/DOtH1CItnqyfvlX+apwUNs1iCJiIgUR3r6TArMENAC72GrSSpX16bdx3CJVxJe5ePJb3IwPtlJ1YmIiBQuhSK5ugq18B22isQqoTbNrgYzL2dO5seP/sOmg9e3C6uIiEhxolAk1+blj9+QZSTXu9uu6ym+5NicISyJOlr0dYmIiBQihaI8iIyMJDg4mNDQ0GsfXFq5euDTby4XWw2167rf+Cu+SwYw7ZcdepmsiIiUWFponQ9lbaH1lWT8FonLypcw/muPjShLfRY2mcFb9zTH1aS8LSIixYMWWovDuHUYjrXPLLIMbjbt87Nu4avoUwyauZWktEwnVSciInJ9FIrkupga98Rl0PekuWa/LPGDzD4sstwEwG8Hz3LvtI2cSLjkzBJFRETyRaFIrl9QWzweX8WhRsP41Njbpmt/XAo9Izew60T+3vgsIiLiLApFUjCV6lPvvnEsfKw9lcrZ3k47k5xO3483sWpPnJOKExERyTuFIikUIYHlWTwsnBuqlLNpb5O1jYNfjGLuxsNOqkxERCRvFIqk0AT6e7FoaHva1vUHoKnhMB+5TuZxl2VU+PEJ3l0Wg8Wihx1FRKR4UiiSQuXn5crsR9rwaGMjn7u9h5chHYC7Tb/TecvjjJ67hrRMvUxWRESKH4UiKXTuLiZebpWBvzHVpj3MuJehh4bz9PSlnEtJd1J1IiIiuVMoEocwBPfA9NA3ZLp427TXN57gzbPP8NzUeRw+k+Kk6kREROwpFInj1LsZ18E/k+5Z1aa5iiGBSZf+y/sffUTUsQtOKk5ERMSWQpE4VrWmuD++ivQKN9o0lzOkMcUyjq8/e4+1+884qTgREZHLFIrE8coH4v7YL2TUbG/T7GKw8I7xIzbPfYXvYk44qTgREZFsCkVSNDwr4DZwCVmNett1PWdaQPw3/2HupiNOKExERCSbQpEUHRd3XPp8hrndk3Zdg12WU275cKau2I3Vqr2MRESk6CkUSdEyGjF1fQtLxBt2Xb1MGzi0eg6vL9utTR5FRKTIuTi7ACmbjB2exlquMpalwzFZszdznJd1C4stHWDDURIvZvLuvc1wNSm3i4hI0dBvHHEaQ/MHMD24kEyjBz+aQ3k1axBgAODb7ScYOjdKu1+LiEiRUSgS56p/K66PrSSr5wyMRpNN16q98Tz82RYSL2U6qTgRESlLFIrE+ao1pVurunzycGs8XG3/k9xy9DwPfbyeM8l6LYiIiDiWQpEUG10aVmHeo2H4elxe6hZoiCPy/GN8OHUisecvOrE6EREp7RSKpFhpXdufhY+3o7KPO5VIZK7rOwQZz/BG2jjmRL7OvtPJzi5RRERKKYUiKXYaVfdl0SNN+cLzfWob4wAwGay8ZJ7Gz9NHE3X0vJMrFBGR0kihSIqloCoVCQoOs2t/mgXs/vwJ1u6Lc0JVIiJSmikUSfFkcsHz3mmktR1p19Xf+BPJXwzgh+ijRV6WiIiUXgpFUnwZDHjc/hoZEW/bdd1t3ITfkn58+dtuJxQmIiKlkUKRFHtuHYaT1etTsv61AXsH4y6Cf3mQz37arPeliYhIgSkUSYngEtIHY7+vSTd62rQ3Mx7h5o39mbJopd6XJiIiBaJQJCWGsf7NuD26nIsu5W3a6xjjuP+PwUyY9w2ZZotzihMRkRKvzIWiXr16UaFCBe69915nlyLXwVCjJV5DV5HsEWDTXsWQwGOHnmLiJ5/pfWkiInJdylwoGjFiBHPmzHF2GVIQlW7AZ9ivJPo2sGn2NVwi+MRXPPz5FpLS9L40ERHJnzIXijp37oyPj4+zy5CC8q2O3xO/cKFyaE7TVsuNPJv5BFuOnOf+j3/X+9JERCRfilUoWrduHd26dSMgIACDwcCSJUvsjomMjKR27dp4eHgQFhbGli1bir5QKR48y1PhsWWcr3U7Bwjk0Yz/kIY7ALtPJdFn+ka9L01ERPKsWIWi1NRUQkJCiIyMzLV/4cKFjBo1ijFjxhAdHU1ISAhdu3YlPj4+55jmzZvTpEkTu6+TJ0/mu5709HSSkpJsvqSYcfXAf8B8LA8vw61cRZuuo+cucu/0jeyP0/vSRETk2gzWYrrBi8FgYPHixfTs2TOnLSwsjNDQUKZOnQqAxWIhMDCQp556ihdeeCHPY69Zs4apU6fyzTffXPW4sWPH8tprr9m1JyYm4uvrm+frSdE4di6Vhz7bTOz5SzbtQZ5pTO7fnuZ1qzupMhERcaakpCT8/Pyu+fu7WM0UXU1GRgZRUVFERETktBmNRiIiIti0aZNDrvniiy+SmJiY8xUbG+uQ60jhqFXRm2+GtqdB1ctrxrxIY7L5LdJn92bjn4edWJ2IiBR3JSYUnT17FrPZTNWqVW3aq1atyunTp/M8TkREBH369GH58uXUrFnzqoHK3d0dX19fmy8p3qr6evDV4+1oGVQeF7KY5jqR5sZDhBl2U+GrnqzcvMPZJYqISDHlcu1DSpeVK1c6uwRxMD8vV+YNDuP3yQ/TKWVnTnsjwzG8l9/Dd6mf0/3mm5xYoYiIFEclZqaoUqVKmEwm4uLibNrj4uKoVq2ak6qS4srLzYWO9z9HkqmCTXuQ4Qzt1/Zj4dLv9b40ERGxUWJCkZubG61atWLVqlU5bRaLhVWrVtGuXTsnVibFlWvN5pR74lfOudWwaa9kSOKO6CF8vuBLvS9NRERyFKtQlJKSQkxMDDExMQAcOXKEmJgYjh8/DsCoUaP45JNPmD17Nnv27OGJJ54gNTWVQYMGObSuyMhIgoODCQ0NvfbBUqwYK9XF/6nVxHnb7379wL4RfPT5p2Rk6X1pIiJSzB7JX7NmDV26dLFrHzBgALNmzQJg6tSpvP/++5w+fZrmzZszefJkwsLCiqS+vD7SJ8VQWhInZ9xDwHnbzT7TrS5Mq/QSgx97mnLuZW6JnYhImZDX39/FKhQVdwpFJVxmGqc/u59qp1fbNGdZjUzyGcXAoc9RsZy7k4oTERFHKXX7FIkUmKsH1YZ8TXytu22aXQwWnkkez+wpY/RaEBGRMkyhSMoWkytVBszh3I332zQbDVZGpU9j8dTn2Htar3MRESmLFIqk7DGaqPjAdBJChth19cz6iQHTf2XLkfNOKExERJxJoSgP9PRZKWQwUL7n+6S0/U9O0ymrPw9mvkRcmiv9P9vMit1xVxlARERKGy20zgcttC6dLq2ZSPq6D7nn0kscsl7e08hogHd6N6NvaKATqxMRkYLSQmuRPPLsPBKPkdHUadjCpt1ihecW7SRy9UHtfi0iUgYoFIkAHr4Vmf5QK/q2rmnX98HPe3jru+3a/VpEpJRTKBL5PxeTkXfvacYTnev9o9XKGy4zuXnbcF5YsFG7X4uIlGIKRSL/YDAYeP72hrxydzBg5QWXBTzksor2pt08sO9pnp61mtT0LGeXKSIiDqBQJJKLRzvUYVnLaIa6LMtpa2E8yIjjI3lixk+cT81wYnUiIuIICkV5oEfyy6YmHXuQ4V7Bpq2R8ThjzvyHJz5ayl8XtPu1iEhpokfy80GP5JdBZ/aRMbM7bhdP2zT/Za3ESLfXeHtwD26s6uOk4kREJC/0SL5IYajcALchP5PpG2TTXNNwlo8yXuKFaV+y7ah2vxYRKQ0UikSupUJtXAf/TKZ/fZvmKoYEPrOOZdyn81m1R7tfi4iUdApFInnhG4Droz+TVbWZTXMFQwqzTG/yybx5fL0t1knFiYhIYVAoEskr74q4DFqGuWaYTbOP4RIzXd5h2bdzmb72kHa/FhEpoRSKRPLDww/Tw4ux1O1i0+xpyOAT1w+I+XkOb/2wR7tfi4iUQApFIvnl5o3xwYVYG95l22ww857rDL7+7Q+e/XoHmWbtfi0iUpIoFOWB9ikSOy7uGPrMhqZ9c5rSra4MzRxJIuVYvP0EQ+Zs42KGdr8WESkptE9RPmifIrFjscDyZ7FEzeHxjJGsMLe06W4eWJ6ZA0Op4O3mpAJFRET7FIkUBaMR7pqA8bFfebD/43i42v5fKiY2gT4fb+JkwiUnFSgiInmlUCRSUAYDVA+hS8MqfDG4LX6erjbdB+NTuOejDRyIS3ZSgSIikhcKRSKFqFWtCnwztB3V/Txy2vqaVjPq4iTum/YbUce0+7WISHGlUCRSyOpX9WHRE+2pV9mbu42beMflU/q4rONN8wQe/mQDP+w85ewSRUQkFwpFIg4QUN6TpRHJTHT7CKMh+1mGO01bmGV8nbHzf+VjbfIoIlLsKBSJOEg5TzdMJhebtlDjfr5zf5kffvqBl5fsIkt7GYmIFBsKRSKOcmNXDA99g9XV26a5uuE8X7u9Tvq2uQyZs43UdO1lJCJSHCgU5YE2b5TrVucmDI/8iNWvpk2zuyGTD1w/5qZDH/DA9PXEJaU5qUAREfmbNm/MB23eKNct9Sx8PRCOrrfr2mQO5nXP0Xz4SAQNq+m/KxGRwqbNG0WKE+9K0H8xhD1h19XOtJtP0kfz8rT5rD9wxgnFiYgIKBSJFB2TK9zxDvSchsXkbtNV03CWubzCt7Mn8tXWWCcVKCJStikUiRS15g9ifORHzOWq2zR7GjL40GUqx5a8zvhf9umRfRGRIqZQJOIMNVphGroOS2Bbm+YMq4nfLY2Y8utBnlkYQ3qW2UkFioiUPQpFIs5SrgrGAd9jbf1oTtOYrIFEWRsAsCTmJA9/toXEi5nOqlBEpExRKBJxJhc3DHdPgG6T2Bv0IAstt9h0bz5ynt7TNhB7/qKTChQRKTsUikSKg1YDafjIND4bEIqXm8mm69CZVHpF/kbMsXNOKk5EpGxQKBIpRro0rMJXj7ejio/t02m9074l8fPerNq+30mViYiUfgpFIsVMkxp+LBkeToOqPgDcZNzB8y5f0skQQ93F3fj2p5VOrlBEpHRSKBIphgLKe/L1E+24p3Y6U1ynYDJkP55fx3ia2zb148s5H2G26JF9EZHCpFCUB3r3mTiDr4cr795RA5OLm017OUMa9x9+kR8mP83F9AwnVSciUvro3Wf5oHefiTNYE2I582kfqqTssevb7BZGvcfmU6lSJSdUJiJSMujdZyKlhKF8IFVGrOZ4ze52fWEZm0mO7MSx/TucUJmISOmiUCRSErh6EvToHI6GvkLWv/5vW8f6FxXm386edd84qTgRkdJBoUikpDAYqH3Xf4jrPp8EfGy6fLlIg1WD+XPhGNAdcRGR66JQJFLC1Gh5B+bBv3LYVNem3Wiw0njPRA5E3oM1PdlJ1YmIlFwKRSIlUMWaN1L9mbVs9e5s11f/7CqOTLqTTL1MVkQkXxSKREooz3K+tBy1mBU1hmO2GnLaLVYDbyR0ZdCsbSSl6WWyIiJ5pVAkUoKZTEZuHfI2K1tGkmj1AmB8Vh9WW1rw28Gz3DttIycSLjm5ShGRkkGhSKQU6NqjHzvvXMInlm5EmnvktO+PS6Fn5AZ2/pXgvOJEREoIhSKRUqJjWBhthkylUjkPm/YzyencO30Ti6L+AnOWk6oTESn+FIpESpGQwPIsHtaeG6qUs2nPyLKwYtEnnP4gjMzzx5xUnYhI8aZQJFLKBPp7sWhoe9rVrZjT1tKwn4mukVS7dJDkqZ25cGibEysUESmeFIpESiE/L1fmPNqGge1rU8twmk/cxuNhyH4Szd9yHve5d3F002InVykiUrwoFOVBZGQkwcHBhIaGOrsUkTxzNRkZ270xz9/ZhAvYvgDRizQCfxpE9KLxTqpORKT4MViteidAXuX1Lbsixc2fh46RNu9BWll32fX9VrUfYUMm4+ri4oTKREQcL6+/vzVTJFIGNK5Xi1ojfmSd5812fR3ivmDrB705l5DohMpERIoPhSKRMqJSeV/aPfsNa6sNtOtrn7aW2Eld2X3waJHXJSJSXCgUiZQhri4mOg2dxNZmr5NpNdn0NbfuwXPu7fy8fqOTqhMRcS6FIpEyKLT3CI7fPotUPG3a6xhO0WrlfXy6YCGZZouTqhMRcQ6FIpEyql677qQ//CNnjZVs2isZknho73AmTx3PuZR0J1UnIlL0FIpEyjD/ui3we2otpzxvsGn3MGSSGB9L96kb2HVCC7BFpGxQKBIp41wr1KT6yDWcqtwhp+2TrDuZY+7KiYRL3DNtI4u3/+XECkVEioZCkYiAuw/Vhy7hXIMH+NXQlrezHszpSs+y8MzCHbz+/W6ytM5IREoxhSIRyWZypeL902g6YhGhtSvZdX++4QgPfbZZ64xEpNRSKBKRywwGKpcvxxdDwhjQrpZdt/nIRp6a/KXWGYlIqaRQJCJ2XE1GXuvRhPfubYabKfvHxA2Gv/jM7QOmpb/AO9M+Zcn2E06uUkSkcCkUicgV9W0dyFdD2xHsc4lZbu/ha7iIn+Ein5veYvXXU3ljmdYZiUjpoVAkIlfVPLA8ixqupqbhbE6bm8HMJLeP8Nj0If0/1TojESkdFIpE5Jo8u72H5cY77dpHu35F99h36TVlrdYZiUiJp1AkItfm5o3x/nkQNtSu6wGX1bxx8Q0GTF/J0hitMxKRkkuhSETyxmiCO96FruOwYrDp6mTayTzDWMZ9uUrrjESkxLruUJSZmUlsbCz79u3j/PnzhVmTiBRn7YZh6DsHq4uHTXMj43EWu49h44Y1PPz5Fs6nZjipQBGR65OvUJScnMy0adPo1KkTvr6+1K5dm0aNGlG5cmVq1arFkCFD2Lp1q6NqFZHiIrg7hoE/YPWy3eSxuuE8X7m9jsuRX+k25TetMxKREiXPoWjChAnUrl2bmTNnEhERwZIlS4iJiWH//v1s2rSJMWPGkJWVxW233cbtt9/OgQMHHFm3iDhbzdYYBq+AirYvk/UxXOJz1/cJT17OvdM3ap2RiJQYBqvVas3LgQ888AAvv/wyjRs3vupx6enpzJw5Ezc3Nx555JFCKdLZIiMjiYyMxGw2s3//fhITE/H19XV2WSLFw8Xz8OWDcHyTXVe/jBfZYGnK453q8lzXhpiMhlwGEBFxrKSkJPz8/K75+zvPoUjy/ocqUuZkpsGSJ+DPb3OavjV3YFTmE/D/Rdkd61di6gMt8fNydVKRIlJW5fX3d4GePjtz5gy7d++2a9+9ezdnzpwpyNAiUpK4esA9n0H4SAAOe7fg+czH4B9Pqa0/cJbukb+xPy7ZOTWKiFxDgULRk08+yYULF+zaL1y4wNNPP12QoUWkpDEa4dbXoPen1H1yMa/1aoGryfZ22bFzF+kZuYGfdp1yUpEiIldWoFB05MgRwsPD7drDw8PZtWtXQYYWkZKqWR/wrMCDYUEsGNKWSuXcbbozMtIZOi+aCb/sw2LR3XsRKT4KFIpymyX626VLlwoytIiUAq1r+7PsqQ6EBJYHoDrnWOE2mluN25j860Eem7uN5LRM5xYpIvJ/BQpFzZo1Y9asWXbtc+bMoWnTpgUZWkRKiWp+Hix8rC0PNS/PTLf3qGOM42PXDxlo+omVe+LpGbmBQ2dSnF2miEjBnj47ceIEPXr0wMfHh5YtWwIQHR1NcnIyS5YsoWbNmoVWaHGgp89ErpM5E+sXfTAcXm3T/FnWHbyV1Q9vdzcmPdCcmxtWdVKBIlKaFekj+atWrcp5Cq1Ro0ZEREQUdMhiSaFI5DqZM+H7kRAzz67rJ3MoIzOHkW5w5z+3NWBY53oYDNrPSEQKj/YpcgCFIpECsFph3fuw+i27ru2WGxic8Szn8OPOptV4/94QvN1dnFCkiJRGRRKK+vfvz9y5cwkNDbX5m53VasVgMLBly5brHbpYUigSKQQ7FsLS4WCxXWB93FKZgZnPc9gaQMNqPszo35qgil5OKlJESpMiCUXHjx8nJSUFb29vm/aUlBRSU1Np06bN9Q5dLCkUiRSSI+vgy4cg3faFsQlWb4ZkPMtWa0P8PF2JfLAlHepXusIgIiJ5UyQ7Wo8ePZoLFy5Qq1Ytm6+EhAQ+/PDDggwtIqVZnZvg0V/AL8imubwhlXlub9PNuJHES5k8/PlmPl1/GN3lF5GioM0bRcQ5qjSEwSuhenObZndDFlPcpvKE6TssVitv/rCHUV/tIC3T7Jw6RaTM0OaNIuI8PlVh0HK48Q67ruddv+Rtl09xIYvF209w7/SNnEjQzxURcRxt3igizuXmDfd/AaFD7LoCDOex/v+lsrtOJNF9ym9sPnyuqCsUkTJCmzfmgxZaiziQ1QqbIuGXlwEryeUbcsuFF4lPd7U5zMVo4NVuwfRvW0v7GYlInmjzRgdQKBIpAruXwsrXYOAyDqb5MGROFEfOptoddl/rQF7v2Rh3F5MTihSRkkSbNzqAQpFIETFngil7hijxUiYjv9zO6n1n7A5rEVSe6Q+1oqqvR1FXKCIliEMeya9Rowbdu3fntdde4/vvv+fEiRMFLlRExI7p8i0zP09XPh0QyvAu9QAoTzJ3GzcBsP14At2m/Eb08Ss/9CEiklf5mimaMmUK0dHRREdHs2fPHsxmM5UrV6Zly5a0atWKli1b0rJlS2rVquXImp1GM0UizvXj9qNUXnIfrQ17mZzVkwlZfQADbiYjb/RszH2hQdccQ0TKHoffPktPTycmJiYnJEVHR7N7924yMzPJysq67sKLM4UiESeyWODbwbBrUU7Tt+YOvJA5hAyyZ5b6t63Fq92CcTUV6MFaESll8vr7+7rfuOju7k5YWBgtW7bk559/JjMzkyNHjuDm5na9Q4qIXNmhX20CEUBv029U5zyPZ44kiXLM/f0Y++KSiXywJZV93J1UqIiUVNf116m0tDQWL15Mv379qFy5MoMGDcJkMjF37lzOnLFfDCkiUmD1I+Cu8WCw/bHVzrSbRW6vUdOQ/bNny5Hz3Dl5PRsPnXVGlSJSguXr9tnChQtZtGgRP/74Iz4+PvTq1YvevXvTuXNnTKbS/1isbp+JFAP7f4GvB0Km7WP6Z6x+PJIxmj+sdQEwGuDpW+rz1M31MRm1n5FIWeaQNUVGo5GAgABefvllBg8ejIvLdd99K5EUikSKiZMxML8vpMTZNF+0uvN05pOstLTKaWtXtyIT72+ux/ZFyjCHhKJOnToRExNDcnIynp6eNGvWzObJsyZNmpTqoKRQJFKMJMTCF33gzB6bZrPVwOjMx/nWclNOW0VvNz68rzk33Vi5qKsUkWLAoU+fHThwgKioKJsnzxISEnB3d6dp06Zs2bKlQMUXVwpFIsVMWiIs7A9H1to0X7K6cXP6eE5R0aZ9WOd6jLr1Rlz0dJpImVLkO1ofOXKEbdu2sX37dt5+++3CGLLYUSgSKYayMuD7EbBjvk3zz8aOPH7xCbvDW9eqwOQHWhBQ3rOoKhQRJyv0UHT8+HGCgvK+MdqJEyeoUaNGno8vCRSKRIopqxWWPgkx82ya364+mRlHKtkdXt7LlfF9QrilUdWiqlBEnKjQX/MRGhrK448/ztatW694TGJiIp988glNmjRh0aJFVzzOWWJjY+ncuTPBwcE0a9aMr7/+2tkliUhhMBjg1tfBw8+m+UXDLF6+swGuJtunzxIuZvLo7G28uWw3GVmWoqxURIqxPM8UnTt3jrfeeovPP/8cDw8PWrVqRUBAAB4eHly4cIHdu3fz559/0rJlS1555RXuvPNOR9eeb6dOnSIuLo7mzZtz+vRpWrVqxf79+/H29s7T+ZopEinmNn0EP7+Y/c9+QXDb6xDck5i/EnlqQTSx5y/ZnRISWJ6pD7Qg0N+riIsVkaLisDVFly5dYvny5axfv55jx45x6dIlKlWqRIsWLejatStNmjQpcPFFJSQkhGXLlhEYGJin4xWKRIo5cyZ8GgGN7oZ2T4Lr5XVDiZcyeWHRTn7cddruNB8PF96/txm3N6lelNWKSBEp9Ntnf4uPj8fNzY0HH3yQxYsX89NPPzFv3jyeffbZAgeidevW0a1bNwICAjAYDCxZssTumMjISGrXro2HhwdhYWHX/aRbVFQUZrM5z4FIREoAkysMWQ03jbYJRAB+nq581K8lb/RojNu/nj5LTsti6LxoxizdRVqmuSgrFpFiJF+haMGCBdx444306NGDdu3a0bp160J9rUdqaiohISFERkbm2r9w4UJGjRrFmDFjiI6OJiQkhK5duxIfH59zTPPmzWnSpInd18mTJ3OOOX/+PA8//DAzZswotNpFpJgwXvnHmsFgoH+72nw7rD11KtnfNp+96Rj3TNvIkbOpuZwtIqVdvm6fNWzYkHbt2vHiiy8SGxvL888/T0hICJ999lnhF2YwsHjxYnr27JnTFhYWRmhoKFOnTgXAYrEQGBjIU089xQsvvJCncdPT07n11lsZMmQI/fv3v+ax6enpOd8nJSURGBio22cipUBKehYvLf6DpTEn7fq83Uy83bspPZqXridoRcoqh9w+O3z4MGPGjOHGG2/klltuYd68eXz55ZcFLjYvMjIyiIqKIiIiIqfNaDQSERHBpk2b8jSG1Wpl4MCB3HzzzdcMRADjxo3Dz88v50u32kRKuLTEnH8s5+7CxPua8+49TXF3sf1RmJphZsSXMbz47U4uZeh2mkhZka9QlJWVhZfX5Sc0GjZsiMVi4fRp+4WLhe3s2bOYzWaqVrXdV6Rq1ap5vv6GDRtYuHAhS5YsoXnz5jRv3pw//vjjise/+OKLJCYm5nzFxsYW6DOIiJOcPZD9SpBZd4HlcsgxGAzcFxrEd0924IYq5exOW7Allp6RGzgYn1yU1YqIk+T7RWWzZ88mPDycZs2aUa5cOVxcXLh48aIjait0HTp0wGLJ+54k7u7uuLu7O7AiEXGojIvw65uw5WOwZGW3bZ8HrQbYHNagmg/fPRnOmKV/8nXUXzZ9++KS6TZlA2/0bMK9rWoWVeUi4gT5minq2LEjb775Jh06dKB8+fLUr1+ftLQ0PvvsM1avXk1ysuP+NlWpUiVMJhNxcbZvxY6Li6NatWoOu66IlGAmNzi8+nIgAlj1us1ttL95ubnwfp8QJvQNwcvNZNN3KdPMf77ewaivYkhNz7I7V0RKh3yForVr15KYmMi+ffuYN28evXr1olOnTkybNo1bbrmFChUq0KhRI4cU6ubmRqtWrVi1alVOm8ViYdWqVbRr184h1xSREs7kAl3/9S7Gi2dh3ftXPKV3y5p892QHGlbzsev7NvoE3af+xp5TSYVdqYgUA/m+fQZQv3596tevz/3335/T9s8Xwl6vlJQUDh48aDNmTEwM/v7+BAUFMWrUKAYMGEDr1q1p06YNEydOJDU1lUGDBl33NfMiMjKSyMhIzGYtuBQpcep1gQZ3wr7ll9t+nw6tBkHFermeckOVciwZHs7ry3Yzf/Nxm75DZ1LpGbmBMd0a80CbQAwGQ65jiEjJk+8drR1pzZo1dOnSxa59wIABzJo1C4CpU6fy/vvvc/r0aZo3b87kyZMJCwsrkvq0o7VICXXuEESGgSXzcluDO+GBBdc89fsdJ3nx2z9IyeW2WbeQAN7u1QQfD9fCrFZECpnDXvNRlikUiZRgv7wMG6fYtvVfkj2TdA3HzqXy5Pzt/HHCfi1SrYpeRD7YkiY1/HI5U0SKA4e95kNEpES6aTR4VbJt++lFMF974XStit5880Q7Bravbdd37NxFen+0kVkbjqC/Y4qUbApFIlI2ePjBLa/Ytp3ZA1Ez83S6u4uJsd0b83H/Vvh62C7HzDBbGPv9bp5csF2bPYqUYApFIlJ2tOgPVZvatq1+Cy6ez/MQXRtXY/mIjrQIKm/X98POU9w7fSMnEy4VsFARcQaFojyIjIwkODiY0NBQZ5ciIgVhNMEd79i2XboAa9/N1zA1K3jx1ePtePymunZ9f55MovvUDUQfv1CQSkXECbTQOh+00FqklPjqYdi99PL3BhMM2wSVG+R7qF/3xjHiyxiS02zXJrmZjIzr3ZR7tAu2iNNpobWIyJXc+gaY/vEKH6sZfv7vdQ11c8OqLBkeTp1K3jbtGWYLz369g7eX78Fs0d89RUoChSIRKXsq1IL2T17+3tUbgtravCw2P+pVLseSYeF0rF/Jrm/GusMMnr2V5LTMXM4UkeJEoUhEyqYOo8CnOjS7H57alv3IvtF07fOuwM/LlZkDQxkUXtuub/W+M/T6aCNHz6YWoGARcTStKcoHrSkSKWXSErMf1S9kX245zitLd5Fptv3x6ufpyrR+LWl/g/2Mkog4jtYUiYhciwMCEcD9bYL4YnBb/L3dbNoTL2XS//MtzN101CHXFZGCUSjKAz2SLyL51aaOP0uHh9Owmo9Nu9li5ZWlf/LS4j/INFucVJ2I5Ea3z/JBt89EyhCLBYwF/3tjanoWzyyM4ZfdcXZ9YXX8mfZQK7sZJREpXLp9JiJyPS5dyH4n2vw+UAh/Z/R2d2H6Q6146uYb7Po2HzlPj8jf2B+XXODriEjBKRSJiED2i2G3fgaTW8LvH8HBlbD3h0IZ2mg08OxtDZj8QAvcXWx/7Maev0SvyA2szGUmSUSKlkKRiAiAOQPWT4BL/3gP2i8vQ1Z6oV2ie0gAXw9tRzVfD5v21AwzQ+Zu46M1B9GKBhHnUSgSEQFw84JbX7Ntu3AEfp9WqJdpVrM83z0ZTkhgeZt2qxXe+2kfzyyMIS3z+jaRFJGCUSgSEflbk3sgsK1t27oPICW+UC9TxdeDhY+1pVeLGnZ9S2JOct+M34lLSivUa4rItSkUiYj8zWCA28fZtmUkw6rXC/1SHq4mJvQN4YU7GmIw2PbtiE2g+9Tf2PlXQqFfV0SuTKEoD7RPkUgZUqMlNO9n27Z9HpzaUeiXMhgMDO1Uj08fbo23m+0rRuKS0ukzfRPf7ThZ6NcVkdxpn6J80D5FImVE8mmY0goyUi63BbWHQcuxm9YpJPvjkhk8exvHz1+06xvepR7P3toAo9Ex1xYp7bRPkYjI9fKpBh1H2bYd3wh/LnbYJW+s6sPS4eG0retv1xe5+hCPz4siJT3LYdcXEYUiEZHctR0O5WvZtq14FTIvOeySFbzdmPtoGA+1DbLrW7E7jnunbSQ2l5kkESkcCkUiIrlx9YDb3rBtS4yFjVMde1mTkTd7NuWNHo0x/et22d7TyfSI3MDmQ2fhRHR2SPvxBTh3yKE1iZQVCkUiIlfSqDvU6mDb9tsESHL84uf+7Woz95E2lPdyzWnz4SJ3pv2Az+yb4ZMusGESbJ8LJterjCQieaVQJCJyJTmP6P9jxibzIqx87YqnFKb2N1Ri6bD2dPeP5X2X6Wx2H86brjMJNh7LOcZ88ytQ3v52m4jkn0KRiMjVVG8GLR+2bTMYweLgXacvnoffp1FrYQSTLz5PH5d1eBlsXzmy1XIjA/8IIfFipmNrESkjXJxdgIhIsXfzK9lPnlVuCLe/AzVbOeY6Visc2wBRs2D3d2C++nvXLBg5cziGblPT+KhfS5rU8LvcmZYEJ7dD3U6OqVWkFNI+RXkQGRlJZGQkZrOZ/fv3a58ikbLo7AGoeINj9ilKOQM75kP0HDh38KqHWgwurDC3ZH5WF9ZbmmL5/4S/m4uRsd0a80CbQAwGA/zwLGz9FFo8BLe9CZ4VCr9ukRIir/sUKRTlgzZvFBGHWDEGNky8+jEV6mTfxmvej+gLbjw2J4qzKfYzSb1b1ODtVsl4zLv7cmO5qnDnBxDcvXDrFikhFIocQKFIRBzi3CGY0tK+3eQGjbpBywFQuyMYLy8DPZV4ieFfRBN9PMHmFHcyWOn1EoGWE/bjNeqeHY58qhbyBxAp3rSjtYhIcWExw/6fs18fkpuK9bJDz98q3Qi3vQWj9sK9n2evCzLa/riu7ufJl4+145HwOjbtGbjwcfptpFo97K+z5zuIbAPbv8hevyQiNjRTlA+aKRIRG+YsiJ4F8XvgrvH2/Qmx2S+T3T4Xkk5Al5eg03O5j7V7KexdDq0GQlDbfK1d+vGPU4z+ZqfNa0ACOMubrp9zsykm95PqdoFuE6FC7TxfR6Sk0u0zB1AoEpEch9fATy9C/O7s7wf9BLXagTkT9v8EUbPh4ErgHz9i/YJgRAwYTYVezpGzqTwxL4q9p5P/0Wqlu3Ejb7rPxdeaZH+Sqxfc8iq0ecwhNYkUFwpFDqBQJCIApCfDh00gLeFyW7VmcEMExHwBKXFXPrffIqgf4ZCy0jLNvLp0F19t+8um3Z8k3vSYy51syP3EmqHQfSpUaeiQukScTWuKREQcxd0HOr9o23Z6Z/YrQK4WiAJa2K0NKkweribeuzeE9+5thrvL5eucx5dhacN5JOM/JLtVtj/xr60wvQOseVdrjaRMUygSEbkeoY9CpQbXPs7dF0IHw+Pr4bE1UO9mh5fWt3Ugi4eFU7uil037r5aWtEsaxy9ed9ufZMmEpL8csw+TSAmh22f5oNtnImLj4EqYd0/ufYFh2Y/SN+4Jbt5FWtbfktIyef6bnfy4y/6pt9u8DzLJ63M8k49mN5SrBsM3g2f5Iq1RpCjo9pmIiKPdEAHN+13+3qM8tB0Gw36HR3+BFv2cFogAfD1c+ahfS169OxgXo+0M0C+pN9Dy7FiiAwdiNZiyn55TIJIyTjNF+aCZIhGxYzHD0fVgyYJaHcA1l/2BioGoYxd4cn40pxLT7Pruq5fFf/vdiZ+Xa+4nH/0NqjZRaJISSzNFhSgyMpLg4GBCQ0OdXYqIFDdGE9TtnD1rVEwDEUCrWhX44emOdKxfya5v4SEX7pqynp1/JdifmHQSFjwAkWGwZ5njCxVxIs0U5YNmikSkpDNbrEz99SATV+23e9DMzWTklW7BPBQWlP1SWasVvnwQ9i2/fFBwT7jzfShXpUjrFikIzRSJiIgdk9HAiIj6zHmkDf7ebjZ9GWYLryzZxciFMaSmZ8Gfi20DEcDuJTA1FGLm6/F9KXUUikREyqCO9Suz/OmOtK5Vwa5vacxJekRu4IhHMNxwq/3JaQmw5An47DY49KvCkZQaCkUiImVUNT8PFjzWliEd69j1HYxP4c7ZR1naeCL0mgGe/vYD/LUF5vaCz7sqHEmpoFAkIlKGuZqMvHRXMNMfaoWPu4tN36VMMyMW7uDlI8GkD90ETa6wJ1Ps5v+Ho9vh0GqFIymxtNA6H7TQWkRKs2PnUnliXjS7T9m/PLZpDT8+6teSwDNr4cfnIeHYlQcKbAtd/gt1OzmwWpG800JrERHJl1oVvfl2WHseaBNo1/fHiUTumryeleaW8FQUdJ8C5YNyHyj2dzi4wsHVihQ+hSIREcnh4WpiXO9mjO8Tgoer7a+IpLQsBs/Zxju/HCIr5CF4Mgq6TQa/f4UjF09oP6IIqxYpHApFIiJi555WNVkyPJy6lexfUzJ97SEe/HQz8Rct0GpA9szRP8NRm8FQrnLuA5uztOZIii2tKcoHrSkSkbImJT2L5xft5Iedp+z6ynu58lr3xnQPCcje7DErA3bMhwZ3XTkUrX4bjqyDzi9AnU5gMOR+nEghyuvvb4WifFAoEpGyyGq1MmfTMd78YTeZZvtfGbcFV+XNXk2o4nON15xcugATm0H6/xdyB7X/fzi6SeFIHEoLrUVEpFAYDAYGtK/NV4+3o0Z5T7v+X3bHcduH61gac4Kr/j3792mXAxHA8Y0wpzvMvBMOr9VtNXE6hSIREcmTFkEVWPZUB+5qVt2uL+FiJiO+jOGxuVHEJ6fZn2yxXPmFsn+Ho1l3Zd9aUzgSJ9Hts3zQ7TMRkWzL/zjFK0t2cS41w67PzzN7rVGP5v9fa/S3rAyI+QLWj4fE2CsPXiv88m01kUKgNUWFKDIyksjISMxmM/v371coEhEBzqWk8+p3f+a6CBsgolFV3u7VhCq+/1prlJUBMfNg3XhI+uvKF6jV4f/hqGMhVi1lkUKRA2imSETE3rVmjcZ2D6Zn8xq2s0YAWenZM0fXCke1O8Jtb0JA88ItXMoMLbQWEZEicWfT6qwY1Ym7c1lrlHgpk2cW7mDInCjik/611sjFHVo/Ak9Hw10TwLdm7hc4uh4yLzqgchFbCkUiIlJg/t5uTH2wJdP6taRSOTe7/pV74rj1w3Us3v6X/RNqLu4Q+ug/wlEN2/46naBWewdWL5JNt8/yQbfPRESu7XxqBmO++5Pvd5zMtT+iURXe7tXUfq3R37LSYftcWD8Bkk7AoB8ViqRAtKbIARSKRETy7qddp3h5yS7OptivNfL1cGFs98b0apHLWqO/ZaXDgV+gUbcrXyQtETz8CqliKa20pkhERJzq9ibV+eWZTnQLCbDrS0rLYtRXOxg8extx/15r9DcX96sHolM7YWJTiJpdSBVLWadQJCIiDuPv7caUB1ow/aHc1xqt2hvPrRPW8m10LmuNriYlHhY8kD1T9P3T8OPz2S+bFSkAhSIREXG425tUZ8Uzneh+vbNG/5SVDgsfsn2Mf/N0+OLe7PeriVwnhSIRESkSFbzdmPxAC6Y/1IpK5dzt+v+eNVoUdY1ZI6MLBLW1bz+8Gj65Bc4eKMSqpSxRKBIRkSJ1e5NqrHjmJno0z33W6NmvrzFrZDTBra9Dr4/B9K9bcucPZQejgysdULmUdgpFIiJS5Cp4uzHp/hZ83P/qs0bfXG3WKOR+GLgcvKvYtqcnwhd9YNNHerms5ItCkYiIOE3XxlefNfrP1zt4dPY2TideYdYoMBQeWw3VQ2zbrRb4+UX47snsNUgieaBQJCIiTnWtWaNf98Zz64dXmTXyqwmDfoLGvez7ts+DOT0g5YwDKpfSRqFIRESKha6Nq7Fy1E30zGXWKPn/s0aPzNrKyYRL9ie7ecG9M6HLS/Z9xzfBJ13g9B8OqFpKE+1onQ/a0VpEpGj88udpXlqyizPJ9re+PFyNPNqhDkM71cPHw9X+5N1LYfFQ+5fIunrDExvAv46DqpbiSjtai4hIiXXb/9ca9WpRw64vLdNC5OpDdHp/DbM3HiXTbLE9ILgHPPIz+Na0bW/WFyrUdlzRUuJppigfNFMkIlL0VuyO47+L/8h11gigTiVvnuvagNubVLN9j1pKfPYmj7GboVYHeHgJmHKZWZJSTy+EdQCFIhER50i4mMH7P+/jy62xmC25/9pqVasC/72zIa1q+V9uzEqH1W9D+6fBu2IRVSvFjUKRAygUiYg418H4FN79aS8rdsdd8ZjbG1fjudsbULdyubwNmnkJXD0LqUIpjrSmqBBFRkYSHBxMaGios0sRESnTbqhSjk8ebs3Cx9oSElg+12N++vM0t324jleX7uJsyjX2KDp3CCY1hx0LC71WKXk0U5QPmikSESk+rFYrP/xxivd+2sfx8xdzPaacuwtDO9Xl0Q518XQz2XamJcKnt8LZfdnfh4+AW8Zkv0ZEShXdPnMAhSIRkeInI8vCvN+PMfnXAyRczMz1mKq+7jx7awPuaVUTk9EAFjPMvw8OrrA9sH5XuOdT8NDP+NJEt89ERKRMcHMx8kiHOqwd3YWhnerh5mL/qy0uKZ3nFu3krsnrWbMvHqvVAv517Qc78DN8diucP1IElUtxo5mifNBMkYhI8Xci4RLjf9nH4u0nrvg+2PAbKvLiHY1ocupbWP4fsGTZHuBZAfrOgTo3Ob5gcTjdPnMAhSIRkZJj14lE3vlxL78dPHvFY3q1qMGLwWepsnwIXDpv22l0gTveg9BHHVypOJpCkQMoFImIlDxr959h3PI97D2dnGu/m4uRka3cePzES5jO7rE/IHQw3P6ONn4swRSKHEChSESkZDJbrHwb/Rfjf9nP6aS0XI8J8MxiQcXPqHV2rX1nnZvgrgngXw+MWo5b0igUOYBCkYhIyXYpw8znG44wbc0hUtKz7PoNWHi93GL6Zy3KfYDKjWD47w6uUgqbnj4TERH5F083E8O73MDa0Z0Z2L42LkaDTb8VI6+k3MPTGcPJIJfbZRXrXXnwP76B3ybCwVWQcqZwC5ci4eLsAkRERIpaxXLujO3emAHta/P+z3tZ/sdpm/7vLOEcS6/KDLcJVDUkXO6o2uTKg8bMh0OrLn9frhpUawLVmv7/q1n2NgDaHLLY0u2zfNDtMxGR0inq2AXeXr6HqGMXbNqrcp5nXb7mZtN2KhmSMPeZh6lxt9wHeb8+pMZf/UKuXlAl2DYoVQ0GN+9C+iSSG60pcgCFIhGR0stqtfLzn3G8+9NejpxN/XcvlUkg5IZAJjzUHl+Pf91aS46D8Tde55UNcNN/4OaXr/N8uRatKRIREckHg8HA7U2q8cszN/FGj8ZU9Hb7Zy9nqMDKgyn0/mgjx879KzRZsqD1I1AzNHs2KF+s4FO9oOVLIdBMUT5opkhEpOxITstkxrrDfLzuMBlZFpu+Cl6ufNy/NW3q+NufaDHD+cNweiec3gWn/8j+Sjltf+zfHl0JgaG59509CKdioOm91/9hyjjdPnMAhSIRkbIn6tgFHp+7jbMpGTbtriYDb/dqSp/WgXkbKOUMxP1xOSSd/gPO7gerFf57Ivd1Rad3wdyecPEc3Ps5NO5V8A9UBikUOYBCkYhI2fTXhYsMnr0t112xH+9Ul+e7NsT4r8f78yTzEpw7lP2Umt1Ft8G83pCWmP290RUeWAD1b83/dco4rSkSEREpJDUrePHNE+25pWEVu76P1x7m8XlRpOayGeQ1uXrmHogADqy4HIgALJmw8CE4+lv+ryN5olAkIiKSB+XcXZjxcGuGdKxj17didxx9pm/iZMKlwrtg5xey37v2T1lpMP9+OBFVeNeRHApFIiIieWQyGnjprmDe6d3Ubjfs3aeS6BG5gZjYhMK5mMEAd7wPze63bc9Ihnn3QNzuwrmO5FAoEhERyaf72wQx59E2+Hna7ld0Jjmd+z7exLKdJwvnQkYj9IiEhnfbtl+6kL0A+9yhwrmOAApFIiIi16V9vUosGR5O3Uq2T42lZ1l4cv52Jq08QKE8y2RyyX7yrN7Ntu0pcTCnJyT+VfBrCKBQJCIict3qVPJm8bBwwm+oaNf34cr9jPgyhrRMc8Ev5OIO982DoHa27YnHYU4PvYC2kCgUiYiIFICflyuzBrWhX1iQXd93O05y/4zfiU9OK/iF3LzhwYVQPcS2/dxBmNsr+5aaFIhCkYiISAG5moy82bMJY7oF8+/timJiE+g5dQN7TiUV/EIefvDQYqjc0LY97g/4og+kpxT8GmWYQpGIiEghMBgMDAqvw2cDQynn7mLTdzIxjXumbWTl7riCX8i7IvRfAhVq27b/tRXWjy/4+GWYQpGIiEgh6tKgCt8Oa0/NCp427RczzAyZu41P1h0u+AJs3+rw8FLbF8k2uBM6PV+wccs4hSIREZFCdmNVH5YOD6d1rQo27VYrvLV8Dy8s+sPuJbP5VqF2djDyqghN+0DfOeDqUbAxyzi9+ywf9O4zERHJj/QsMy9++wffRp+w6wur48/0h1pRwdutYBe5cBT8AsFoKtg4pZjefSYiIuJk7i4mxvcJYXTXBnZ9m4+cp+dHGzgYX8DF0RVqKxAVkjIVihISEmjdujXNmzenSZMmfPLJJ84uSURESjmDwcDwLjcw/aGWeLrahpdj5y7S66MN/HbgrGMunpYEMfMdM3YpVKZun5nNZtLT0/Hy8iI1NZUmTZqwbds2Kla033QrN7p9JiIiBbHrRCKPzt5KXFK6TbvJaGBs98b0b1ur8C528TzM6w0nt0PXcdBuWOGNXcLo9lkuTCYTXl5eAKSnp2O1WgtnC3YREZE8aFLDj++e7EDTGn427WaLlVeW7GLsd3+SZS7gAmyA5NMw887sQATw84sQPafg45ZyxSoUrVu3jm7duhEQEIDBYGDJkiV2x0RGRlK7dm08PDwICwtjy5Yt+bpGQkICISEh1KxZk9GjR1OpUqVCql5EROTaqvp68NXj7bizaTW7vlkbj/Lo7G0kpWUW7CKH18CZPbZt3z0NuxYVbNxSrliFotTUVEJCQoiMjMy1f+HChYwaNYoxY8YQHR1NSEgIXbt2JT4+PueYv9cL/fvr5MnsNxaXL1+eHTt2cOTIEebPn09cXCFspCUiIpIPnm4mpj7QkqduvsGub+3+M9zz0UaOn7t4/RcIuR8ixv6r0QrfPgb7f77+cUu5YrumyGAwsHjxYnr27JnTFhYWRmhoKFOnTgXAYrEQGBjIU089xQsvvJDvawwbNoybb76Ze++9N9f+9PR00tMv3/dNSkoiMDBQa4pERKTQLN7+F89/8wcZ/7pt5u/txsf9WxFa2//6B1/1uv0u1yZ3eOgbqHPT9Y9bwpS6NUUZGRlERUURERGR02Y0GomIiGDTpk15GiMuLo7k5GQAEhMTWbduHQ0a2D8m+bdx48bh5+eX8xUYGFiwDyEiIvIvvVrUZMFjbalUzna/ovOpGfT7ZDORqw9yOvE6Xyh78yvQ5jHbNnM6LHgA/tp2nRWXXiUmFJ09exaz2UzVqlVt2qtWrcrp06fzNMaxY8fo2LEjISEhdOzYkaeeeoqmTZte8fgXX3yRxMTEnK/Y2NgCfQYREZHctKpVgcXDwmlQ1cemPcNs4f2f99HunVX0nb6JOZuOciY5/Qqj5MJggNvfhZAHbdszUmDePXB6VyFUX3q4XPuQ0qNNmzbExMTk+Xh3d3fc3d0dV5CIiMj/Bfp78c0T7RjxZQy/7o236bNaYcvR82w5ep6x3/1J27oVubtZALc3qYb/tXbENhqh+5TsILTnu8vtaQkwtxcM+hEq2a9tKotKzExRpUqVMJlMdguj4+LiqFbNfgW/iIhISePj4conD7fm0Q51rniMxQobD53jv4v/IPStlTz8+Ra+2hZL4qWrPLFmcoF7PoUbImzbU+NhTg9I0J0QKEGhyM3NjVatWrFq1aqcNovFwqpVq2jXrp0TKxMRESk8JqOBV+4O5qvH23Fvq5r4eFz5po7ZYmXd/jM8981OWr+5gkdnbWXx9r9Izu2Rfhd36DsXgtrbtif9lR2MUuLtzyljitXTZykpKRw8eBCAFi1aMGHCBLp06YK/vz9BQUEsXLiQAQMG8PHHH9OmTRsmTpzIV199xd69e+3WGjmCdrQWEZGilp5lZv3+s3y/8yQrd8eRmmG+5jluLka6NKjM3c0CuKVRFbzc/hGs0pJgTvfLGzv+rcGd8MCCQq6+eMjr7+9iFYrWrFlDly5d7NoHDBjArFmzAJg6dSrvv/8+p0+fpnnz5kyePJmwsDCH1hUZGUlkZCRms5n9+/crFImIiFOkZZpZvTeeZTtPsWpvHGmZ19792tPVxM2NqtCtWXU6N6iCh6sp+xUgM++AM3uzD6rUAB5eAr4Bjv0ATlIiQ1Fxp5kiEREpLlLTs1i1N55lO06yZv8ZMrKuHZDKubtwa3BV7m5WnY7VsnCbcye4+0L/xeBdet/woFDkAApFIiJSHCWnZbJidxzLdp5i/YEzZJqv/avd18OF++ob6BRyA2GNauNqKjHLjPNNocgBFIpERKS4S7yYyc9/nub7nSfZeOgcZsu1f81X8HLl9ibV6dasOmF1K2IyGiDuTzj6G9TtDJVuzN7zqIRSKHIAhSIRESlJzqWk89Ofp1m24xS/HzlHXn7jVyrnzp1NqzHC+BUVoyZlN/pUhzqdoG6n7P/1q+HYwguZQpEDKBSJiEhJFZ+UxvI/TrFs5ym2HbtwzeMXuY2llXF/7p0V62fPINXtBLU7gGeFwi22kCkUFSI9fSYiIqXJyYRLLP/jFN/vPMWO2AS7/nJcJMb9MVwM1168jcEI1ZtfDkmBbcHVo7BLLhCFIgfQTJGIiJQ2secvsmznKZbtPMmfJ5MAqMIFRrssJNy0iwDD+fwNaHKHTqPhptEOqPb6KBQ5gEKRiIiUZofPpLBwWyyfrT9ClsUKWKljOE24cVf2l2kPvqRce6DuU6Dlww6vN68UihxAoUhERMqCP08mMvrrnew+lWTTbsRCU+NRRtY7xU0uuzDFboasNPsBRuyECrXs27PSYVp7qNnm8qJt3+oO+hSXKRQ5gEKRiIiUFZlmC9PWHGLKrwdy3feoXmVvPujVgBbshyNr4fCa7FeHlK8FI2JyH/TIeph9t21bpQa2i7Y9/Ar7oygUOYJCkYiIlDV7TiUx+psd7DqRZNdnNMDgjnUZdeuN2a8PuZQAibFQrWnug/36Jqx7/8oX868HT0cXTuH/kNff36V3+8pCFBkZSXBwMKGhoc4uRUREpEg1qu7L4mHh/Oe2G3E12W7gaLHCjHWHuXPSeqKOnQfP8lcORJA9m3Q1dToWuN6C0ExRPmimSEREyrJ9p5P5z9c7+ONEol2fwQCPhtfh2dsa4Olmyn2Asweyg9HhNdm30tL/NU6fWdC4V2GXrdtnjqBQJCIiZV2W2cKM9YeZuOIAGWb7fYzqVPLmvXubEVrb/+oDmbPg1A44siY7JMVuhWf+BO+KhV6zQpEDKBSJiIhkOxCXzH++2Znr5o8GAwxsX5vRXRvg5eaStwGz0sHFvXCL/D+tKRIRERGHqV/Vh0VD2/HCHQ1xc7GNE1YrzNxwlDsmrWfz4XN5G9BBgSg/FIpERETkuriYjAztVI/lT3ekRVB5u/5j5y5y34zfGbN0F6npWUVfYD4pFImIiEiB3FClHN8Mbc9LdzbC3cU+WszedIzbJ61j46GzTqgu7xSKREREpMBMRgNDbqrL8hEdaVWrgl1/7PlLPPjJZl5ZUnxnjRSK8kD7FImIiORNvcrl+OrxdrxydzAervYxY+7vx+g6cR0bDha/WSM9fZYPevpMREQk746cTeW5b3aw9eiFXPsfDAvixTsa4uPh6tA69PSZiIiIOFWdSt4sfKwdY7oF4+lqv6Hj/M3HuX3ietYfOOOE6uwpFImIiIjDGI0GBoXX4aeRHQmrY7+h44mES/T/bAsvfruTpLRMJ1R4mUKRiIiIOFytit4sGNKW13s0xiuX14As2BJL1w/XsWZfvBOqy6ZQJCIiIkXCaDTwcLva/DTiJtrVtX+dx6nENAbO3Mpz3+xwyqyRQpGIiIgUqaCKXnwxOIw3ezbBO5dZo7X7z+CMx8AUikRERKTIGY0GHmpbi59G3kT4DbazRuN6N8XP07FPpOVaU5FfUUREROT/Av29mPdoGG/3ako5dxd6t6zBzQ2rOqWWPL66VkRERMQxDAYDD4YF0alBZcq5OS+aaKYoD7SjtYiIiOPVKO+Jn1fR3zb7m3a0zgftaC0iIlLyaEdrERERkXxQKBIRERFBoUhEREQEUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSERERARSKRERERACFojzRaz5ERERKP73mIx8SExMpX748sbGxes2HiIhICZGUlERgYCAJCQn4+fld8TjnvYq2BEpOTgYgMDDQyZWIiIhIfiUnJ181FGmmKB8sFgsnT57Ex8cHg8FQqGP/nWLLyiyUPm/pps9buunzlm6l8fNarVaSk5MJCAjAaLzyyiHNFOWD0WikZs2aDr2Gr69vqfmPMC/0eUs3fd7STZ+3dCttn/dqM0R/00JrERERERSKRERERACFomLD3d2dMWPG4O7u7uxSioQ+b+mmz1u66fOWbmXt8/6TFlqLiIiIoJkiEREREUChSERERARQKBIREREBFIpEREREAIWiYiEyMpLatWvj4eFBWFgYW7ZscXZJDjFu3DhCQ0Px8fGhSpUq9OzZk3379jm7rCLzzjvvYDAYGDlypLNLcZgTJ07w0EMPUbFiRTw9PWnatCnbtm1zdlkOYzabeeWVV6hTpw6enp7Uq1ePN954g9Ly/Mq6devo1q0bAQEBGAwGlixZYtNvtVp59dVXqV69Op6enkRERHDgwAHnFFsIrvZ5MzMzef7552natCne3t4EBATw8MMPc/LkSecVXEDX+vf7T0OHDsVgMDBx4sQiq88ZFIqcbOHChYwaNYoxY8YQHR1NSEgIXbt2JT4+3tmlFbq1a9cyfPhwfv/9d1asWEFmZia33XYbqampzi7N4bZu3crHH39Ms2bNnF2Kw1y4cIHw8HBcXV358ccf2b17N+PHj6dChQrOLs1h3n33XaZNm8bUqVPZs2cP7777Lu+99x5TpkxxdmmFIjU1lZCQECIjI3Ptf++995g8eTLTp09n8+bNeHt707VrV9LS0oq40sJxtc978eJFoqOjeeWVV4iOjubbb79l3759dO/e3QmVFo5r/fv92+LFi/n9998JCAgoosqcyCpO1aZNG+vw4cNzvjebzdaAgADruHHjnFhV0YiPj7cC1rVr1zq7FIdKTk621q9f37pixQprp06drCNGjHB2SQ7x/PPPWzt06ODsMorUXXfdZX3kkUds2nr37m3t16+fkypyHMC6ePHinO8tFou1WrVq1vfffz+nLSEhweru7m5dsGCBEyosXP/+vLnZsmWLFbAeO3asaIpyoCt93r/++stao0YN665du6y1atWyfvjhh0VeW1HSTJETZWRkEBUVRURERE6b0WgkIiKCTZs2ObGyopGYmAiAv7+/kytxrOHDh3PXXXfZ/Hsujb777jtat25Nnz59qFKlCi1atOCTTz5xdlkO1b59e1atWsX+/fsB2LFjB7/99ht33HGHkytzvCNHjnD69Gmb/679/PwICwsrEz+/IPtnmMFgoHz58s4uxSEsFgv9+/dn9OjRNG7c2NnlFAm9ENaJzp49i9lspmrVqjbtVatWZe/evU6qqmhYLBZGjhxJeHg4TZo0cXY5DvPll18SHR3N1q1bnV2Kwx0+fJhp06YxatQo/vvf/7J161aefvpp3NzcGDBggLPLc4gXXniBpKQkGjZsiMlkwmw289Zbb9GvXz9nl+Zwp0+fBsj159fffaVZWloazz//PA888ECpemnqP7377ru4uLjw9NNPO7uUIqNQJE4xfPhwdu3axW+//ebsUhwmNjaWESNGsGLFCjw8PJxdjsNZLBZat27N22+/DUCLFi3YtWsX06dPL7Wh6KuvvuKLL75g/vz5NG7cmJiYGEaOHElAQECp/cySvei6b9++WK1Wpk2b5uxyHCIqKopJkyYRHR2NwWBwdjlFRrfPnKhSpUqYTCbi4uJs2uPi4qhWrZqTqnK8J598kmXLlrF69Wpq1qzp7HIcJioqivj4eFq2bImLiwsuLi6sXbuWyZMn4+LigtlsdnaJhap69eoEBwfbtDVq1Ijjx487qSLHGz16NC+88AL3338/TZs2pX///jzzzDOMGzfO2aU53N8/o8raz6+/A9GxY8dYsWJFqZ0lWr9+PfHx8QQFBeX8/Dp27BjPPvsstWvXdnZ5DqNQ5ERubm60atWKVatW5bRZLBZWrVpFu3btnFiZY1itVp588kkWL17Mr7/+Sp06dZxdkkPdcsst/PHHH8TExOR8tW7dmn79+hETE4PJZHJ2iYUqPDzcbouF/fv3U6tWLSdV5HgXL17EaLT9MWoymbBYLE6qqOjUqVOHatWq2fz8SkpKYvPmzaXy5xdcDkQHDhxg5cqVVKxY0dklOUz//v3ZuXOnzc+vgIAARo8ezc8//+zs8hxGt8+cbNSoUQwYMIDWrVvTpk0bJk6cSGpqKoMGDXJ2aYVu+PDhzJ8/n6VLl+Lj45Oz7sDPzw9PT08nV1f4fHx87NZLeXt7U7FixVK5juqZZ56hffv2vP322/Tt25ctW7YwY8YMZsyY4ezSHKZbt2689dZbBAUF0bhxY7Zv386ECRN45JFHnF1aoUhJSeHgwYM53x85coSYmBj8/f0JCgpi5MiRvPnmm9SvX586derwyiuvEBAQQM+ePZ1XdAFc7fNWr16de++9l+joaJYtW4bZbM75Gebv74+bm5uzyr5u1/r3++/Q5+rqSrVq1WjQoEFRl1p0nP34m1itU6ZMsQYFBVnd3Nysbdq0sf7+++/OLskhgFy/Zs6c6ezSikxpfiTfarVav//+e2uTJk2s7u7u1oYNG1pnzJjh7JIcKikpyTpixAhrUFCQ1cPDw1q3bl3rSy+9ZE1PT3d2aYVi9erVuf5/dsCAAVarNfux/FdeecVatWpVq7u7u/WWW26x7tu3z7lFF8DVPu+RI0eu+DNs9erVzi79ulzr3++/lYVH8g1WaynZelVERESkALSmSERERASFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhEp5Tp16oTBYGDBggU27VOmTCEgIMBJVYlIcaRQJCKlltVqZfv27VSvXp1FixbZ9EVFRdGyZUsnVSYixZFCkYiUWgcOHCA5OZmXX36ZH3/8kYsXL+b0RUdH06pVKydWJyLFjUKRiJRaUVFReHh4MHjwYHx9ffnxxx8BSEtLY8+ePZopEhEbCkUiUmpFR0fTrFkz3Nzc6NWrF9988w0AO3bsICsrKycU/fXXX/Tu3Zt69erRunVr+vTpQ1xcHCdPnqRfv35XPUZESg+FIhEptaKjo3OCT+/evfnhhx9IT08nOjqaypUrExgYiNVqpUePHtx1110cOnSIbdu28fTTT3PmzBkCAgL44osvrnqMiJQeCkUiUmr9c91Q586dcXV15eeff7ZZZL1q1SrKlSvHo48+mnNex44dadKkCUePHqV169ZXPUZESg+FIhEplQ4fPkxCQkJO+HFxcaF79+4sWrTIJizt3r37mmuL8nKMiJR8CkUiUipFRUXh5uZmM5tzzz338N133/Hnn38q5IiIHYUiESmVoqOjadKkCW5ubjltt956K2azmYyMjJxQ1KhRI7Zv337VsfJyjIiUfApFIlIqjRs3jqioKJs2d3d3kpKSsFqt1KlTB4CIiAiSkpKYNWtWznG//fYbu3btyvk+L8eISMmnUCQiZZrBYGDJkiUsWbKEevXq0bhxY6ZMmULlypXzdYyIlHwGq9VqdXYRIiIiIs6mmSIRERERFIpEREREAIUiEREREUChSERERARQKBIREREBFIpEREREAIUiEREREUChSERERARQKBIREREBFIpEREREAIUiEREREUChSERERASA/wGiDGgh7TltsAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "cic_cens = 0.5 * (cic_edges[:-1] + cic_edges[1:])\n",
    "plt.semilogy(cic_cens, cic_galtab, label=\"galtab\", lw=3)\n",
    "plt.semilogy(cic_cens, cic_halotools, label=\"halotools\", lw=3, ls=\"--\")\n",
    "plt.legend(frameon=False)\n",
    "plt.xlabel(\"$N_{\\\\rm CiC}$\")\n",
    "plt.ylabel(\"$P(N_{\\\\rm CiC})$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **In Development:** Differentiate CiC w.r.t. the HOD parameter $\\log M_{\\rm min}$\n",
    "\n",
    "`galtab` is implemented in JAX, so it is portable to GPU and differentiable\n",
    "(in principal), assuming your HOD model is compatible with JAX. Unfortunately,\n",
    "this requires a few modifications to `halotools` models. For example, let's\n",
    "use the `JaxZheng07Cens` and `JaxZheng07Sats` models, originally implemented\n",
    "for the [JaxTabCorr](https://github.com/AlanPearl/JaxTabCorr) project."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can construct a composite HOD model with our JAX-compatible mean\n",
    "occupation functions, which we call `hod_jax`. This model allows us to\n",
    "differentiate `cictab.predict` with `jax.grad`.\n",
    "\n",
    "*Note:* You shouldn't try using `jax.jit` directly on `cictab.predict`, since it\n",
    "contains some I/O lines that can't be compiled. Rest assured that the primary\n",
    "expensive computations will automatically compile and run on the GPU if available."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkMAAAG0CAYAAAAxRiOnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/av/WaAAAACXBIWXMAAA9hAAAPYQGoP6dpAABacElEQVR4nO3deVxU5f4H8M8wMOACiKLsOuK+IMqaGmlComWaaKmZotd7vXXNJczUezMqfwWuoeHV9LqluaVolksZiWGZJoipKKKiIgguKaAoy8z5/TExNTIDA8wG5/N+veYFnPOcZ76PR5mPZ3mORBAEAUREREQiZWXuAoiIiIjMiWGIiIiIRI1hiIiIiESNYYiIiIhEjWGIiIiIRI1hiIiIiESNYYiIiIhEzdrcBdQHSqUSubm5sLe3h0QiMXc5REREpAdBEFBUVAR3d3dYWek+/sMwpIfc3Fx4eXmZuwwiIiKqhezsbHh6eupczzCkB3t7ewCqP0wHBwczV0NEpqZUKpGSkgI3N7cqf6ESkWUpLCyEl5eX+nNcF4YhPVScGnNwcGAYIhKJ0tJS/PDDD9izZw+++uor9OvXD1u3buWpcqJ6qLp/twxDRER/KCwsxIEDB7Bnzx7s378fhYWFAIBevXph3bp1DEJEDRTDEBGJWl5eHvbu3Ys9e/YgMTERpaWlGutbtmyJPXv2oHHjxgZ9359++gktW7ZEx44dDdovEdUcwxARic6lS5ewe/du7NmzB8eOHYMgCFrb2djYYNeuXWjdurXB3jsrKwuzZ89GZmYmUlJSDNYvEdUewxARNXiCICA1NVUdgM6dO6fXdvHx8QgJCTFIDYWFhYiJicEnn3yCkpISfPXVV1Xe6ktEpsMwREQNmiAIWLRoEf79739DoVDovd2//vUvTJ48uc7vr1AosG7dOrz77ru4desWACAgIAAvvvhinfsmIsOwuP+WrFixAnK5HHZ2dggODsaJEyd0tl2zZg1CQkLg5OQEJycnhIWFaW1//vx5DB06FI6OjmjSpAkCAwNx/fp1Yw6DiCyERCLBO++8gzNnzmDIkCF6bdOvXz/ExcXV+b1/+OEH+Pn5YfLkyeogBAAffvghL8YmsiAWFYa2b9+OqKgoREdHIzU1Fb6+vggPD9f4JfJXSUlJGDNmDA4fPoxjx47By8sLAwcORE5OjrrN5cuX8fTTT6Nz585ISkrCb7/9hnnz5sHOzs5UwyIiC9ClSxds3rwZwcHBVbZr06YNvvzyS9jY2NT6vTIzMzFs2DCEhobit99+01jXu3dvDBo0qNZ9E5HhSQRdVw6aQXBwMAIDAxEfHw9ANdGZl5cXpk6dijlz5lS7vUKhgJOTE+Lj4zF+/HgAwOjRo2FjY4NNmzbVuq7CwkI4OjqioKCA8wwR1VPfffcdJk2ahBs3buhs07hxY/z888/w9fWt1Xvcu3cP8+fPR3x8PMrKyrS2+f777xEaGlqr/omoZvT9/LaYI0OlpaVISUlBWFiYepmVlRXCwsJw7NgxvfooLi5GWVkZmjdvDkAVpvbt24eOHTsiPDwcrVq1QnBwMPbs2VNlPyUlJSgsLNR4EVH99ODBA7zxxhsIDw+vMggBwOeff16rIFRWVob4+Hh06NABn3zyic4g9Mwzz2DAgAE17p+IjMtiwtCdO3egUCjg4uKisdzFxQV5eXl69TF79my4u7urA9WtW7fw4MEDxMbGYtCgQfjuu+8wfPhwRERE4MiRIzr7iYmJgaOjo/rF55IR1Q8KBZCUBGzdqvp6+PCP8PX1xapVqzTaWVlZYeDAgRrL5s2bhxEjRtT4PQ8cOIAePXpg6tSpuHv3bpVtea0QkWVqMHeTxcbGYtu2bUhKSlJfD6RUKgEAw4YNw1tvvQUA6NmzJ37++WesWrUK/fr109rX3LlzERUVpf654tkmRGS5EhKA6dMBzYM/bQH4AriiXtK5c2ds3LgRjRs3ho+PDwDgpZdewvvvv1+j97t9+zYiIyNx4MABvdqHhobq/J1DROZlMUeGnJ2dIZVKkZ+fr7E8Pz8frq6uVW67ePFixMbG4rvvvkOPHj00+rS2tkbXrl012nfp0qXKu8lsbW3VzyHj88iILF9CAjBy5JNBCAA8AOwEMBwSiQRRUVFITU1FUFAQ2rZtCwDo1q0bPv/88xrP+dOyZUv873//w6pVqzBo0CDIZLIq28+fP79G/ROR6VhMGJLJZPD390diYqJ6mVKpRGJiInr37q1zu4ULF2L+/Pk4ePAgAgICKvUZGBiIjIwMjeUXL15EmzZtDDsAIjKYhw8fIiIiAh9//DFOnjypPsqrjUKhOiKk/VYQ1a84qTQeP/xwBEuWLEGjRo0AAE2aNEHXrl3x1VdfVftEa13c3d3xz3/+E5s3b67ysRqDBw+u8vcYEZmZYEG2bdsm2NraChs2bBDS09OFyZMnC82aNRPy8vIEQRCEcePGCXPmzFG3j42NFWQymbBz507h5s2b6ldRUZG6TUJCgmBjYyOsXr1ayMzMFD799FNBKpUKycnJetdVUFAgABAKCgoMN1giESovF4TDhwVhyxbV1/Jy3W1jY2MFAAIAoUWLFsKoUaOEtWvXCtevX9dod/iwIKiiUNWvw4crv8eVK1fqPKb79+8L/v7+6lq1vU6cOFHn96kLpVIpKJVKs9ZAZA76fn5bVBgSBEH49NNPhdatWwsymUwICgoSfvnlF/W6fv36CZGRkeqf27Rpo/UXT3R0tEafa9euFdq3by/Y2dkJvr6+wp49e2pUE8MQUd3t2iUInp6aAcXTU7Vcm4cPHwpubm5a/4136dJFmD59urBv3z5h3bpHeoWhLVsMP6YHDx4Iffv21ahNIpEIrVq1Uv88dOhQw7+xnm7cuCF8/PHHwscff2y2GojMSd/Pb4uaZ8hScZ4horqpuKbnyd82EolqwcqVtzFgQAFKSkpQWlqq/vr5559j3bp1VfYtlYZCofi+2hoOHwb696/tCCp7/PgxhgwZonFqH1DNjF9YWIiZM2cCAE6dOoWePXsa7o2rUVJSgr1792L9+vX49ttv4e3tjdTU1FqfCiSqz/T9/GYY0gPDEFHtKRSAXK7t4uYKSgA3oLrzS/e1QbpZAbgK1cXSlS+DlEgAT08gKwuQSmvR/R8UCiA5Gbh5E3B2LkNc3Ajs3/+1Rptly5Zh2rRpyM3NhZeXF4YPH46dO3fW/k1r4NSpU1i3bh22bNmC33//HQBgY2ODn3/+udL1lERioe/nd4O5tZ6ILFNyclVBCFAFmNYAQgDonv9LNyWA6VDdNaaEZiBS/V8vLk5SpyBU+bZ9GwAroPoVuhsA8PHHH2PatGkAoJ7vrKa369fUnTt38MUXX2D9+vU4ffp0pfUfffQRgxCRHhiGiMiobt7Ut6VbHd5lN4CRsLH5L8rK/pyKw8tLgrg4ICKi9j3rOsX35237I/Gf/3TF3LlzNdauXr3aKHetlpeX49tvv8W6devw9ddf65ztOiwsTH2qjoiqxjBEREblpmfGeemlp+Dt7Q6ZTAZbW1vIZDJcuHCh2ucKBgUFYeTIkRgxYgTatHFVn8pycwNCQup+aqzq2/aVaNp0Ld5/v1mltYYOQhcuXMD69euxadMm3KwmYTo7O9dq7iQisWIYIiKjCglRXbOTk6M9VFRc07Nz5/RKwWX48OFa++zTpw9GjhyJiIiISqHDkBdJ63OK78EDJxw9atj3raBQKLBhwwasXbtW72c0AsD69evhpm8KJSKGISIyLqkUWLZMdapJItEMRBWP6YqLq3wE5/z58+qHKkskEjzzzDMYOXIkhg8fDg8PD5PUru8pPv1PBdaMVCpFcHAwEhMTcfz48Sonn6zw5ptvYsiQIcYpiKiB4jFUIjK6iAhg507gyQyjOiKk/ZqepUuXIjQ0FCtXrkRubi6SkpLw5ptvmiwIAfqf4jPmQZju3btjy5YtOH/+PIKCgqps6+Pjg0WLFhmvGKIGirfW64G31hMZxl9vT6/ump6CggI4OjqatsAnVEwLUN0pvrretl8dpVKJjz76CNHR0dD1K9vOzg4nT55Et27djFcIUT3DW+uJyOJIpfpfW2PuIATU/hSfId29exevvfYaDh48WGW7pUuXMggR1RJPkxERVaE2p/gM5cSJE+jVq1elIOTu7q7x80svvYTXX3/deIUQNXA8TaYHniYjopqc4qtr/66uAs6eXYmZM2dUmkdo5MiRWL16NTw8PPDo0SN4eHjg9OnTaNGiheGKIWogeJqMiMiAanKKr6Yqz3AtAfAigO9RMcO1tbU1Fi9ejGnTpkEikcDX1xfHjx/Hpk2bGISI6ohhiIjIjPSZ4drT81fs2LEDvXv3Vq/18/PDs88+i2effdaE1RI1TAxDRERmos8M17a2q/DrrxK4urbUWBsZGYlevXqZokyiBo8XUBMRmYk+M1yXlLTChQstK60JCgqCjY2N0WojEhOGISIiMzH3DNdEpMIwRERkJpYwwzURMQwREZlNxUNsKyZwfJJEAnh5qdoRkfEwDBERmUnFDNdA5UBkqhmuiYhhiIjIrMw5wzURqfDWeiIiM4uIAIYNM+4M10SkG8MQEZEFMOYM10RUNZ4mIyIiIlFjGCIiIiJRYxgiIiIiUWMYIiIiIlFjGCIiIiJRYxgiIiIiUWMYIiIiIlHjPENEpKZQcOI/IhIfhiEiAgAkJADTpwM3bvy5zNNT9ewsPhKCiBoyniYjIiQkACNHagYhAMjJUS1PSDBPXUREpsAwRCRyCoXqiJAgVF5XsWzGDFW7J92+fRvl5eVGrY+IyNgYhohELjm58hGhvxIEIDsb+PHHymkpNTUVL7/8Mh4/fmzEComIjIthiEjkbt7Ur90///k+9u3bB+Evh5BsbGywZ88eDBw4EPfv3zdOgURERsYwRCRybm76tcvMPIIhQ4agV69e2L59OxQKBWxsbAAAycnJeOaZZ5Cbm2vESomIjINhiEjkQkJUd41JJLpaKAFcB5AMADh9+jRGjx6Nzp07Y9euXepWZ86cQZ8+fZCRkWHskomIDIphiEjkpFLV7fNA5UAkkQgAJHBxiYEqFP3p0qVLWFax4R+uXbuGvn374vjx48YrmIjIwBiGiAgREcDOnYCHh+ZyT08Jdu2S4MaNT/HFF1+ge/fu1fZ19+5dDBgwAAcOHDBStUREhiURBG031NJfFRYWwtHREQUFBXBwcDB3OURGU90M1EqlEt988w0+/vjjao/+WFtbY+3atRg/fryRqyYi0k7fz2+LOzK0YsUKyOVy2NnZITg4GCdOnNDZds2aNQgJCYGTkxOcnJwQFhZWZfvXX38dEokEcXFxRqicqP6TSoH+/YExY1Rfn3wUh5WVFYYOHYqpU6fC2rrqCezLy8sRGRmJxYsXG61eIiJDsKgwtH37dkRFRSE6Ohqpqanw9fVFeHg4bt26pbV9UlISxowZg8OHD+PYsWPw8vLCwIEDkZOTU6nt7t278csvv8Dd3d3YwyBqsEpLSzF16lS89tprek+2OGvWLLz99ttQKpXVNyajUSiApCRg61bVV22TaFb46aef8OWXX3JCTRIPwYIEBQUJU6ZMUf+sUCgEd3d3ISYmRq/ty8vLBXt7e2Hjxo0ay2/cuCF4eHgIZ8+eFdq0aSN88sknVfbz+PFjoaCgQP3Kzs4WAAgFBQU1HhNRQ6FUKoU5c+YIrVq1EgDU+PXaa68JpaWl5h6GKO3aJQienoKgmkJT9fL0VC3X5vbt24JUKhXatm0rLF++XCgqKjJtwUQGUlBQoNfnt8UcGSotLUVKSgrCwsLUy6ysrBAWFoZjx47p1UdxcTHKysrQvHlz9TKlUolx48Zh1qxZ6Natm179xMTEwNHRUf3y8vKq2WCIGiCJRIKYmBjk5+fj999/x7Fjx7B+/XrMmTMHw4cPR5cuXdTzDmmzefNmvPjii3jw4IEJqybdz50TdD53ztnZGWFhYcjKysK0adPQunVr/Oc//8FNfWfoJKpnLOYC6tzcXHh4eODnn39G79691cvfeecdHDlyRK9bdf/1r3/h22+/xblz52BnZwdAFWwOHz6Mb7/9FhKJBHK5HDNmzMCMGTN09lNSUoKSkhL1z4WFhfDy8uIF1ETVKC8vR1ZWFjIyMpCRkYELFy6ov96+fRsAEBgYiH379qFly5ZmrrbhUygAubyqx60oYWd3F5Mm/R86dPCGt7c32rVrh7Zt22L79u2YOHGiRmuZTIbXXnsNM2fORNeuXY1dPlGd6XsBddVXQNYjsbGx2LZtG5KSktRBKCUlBcuWLUNqaiokumeUq8TW1ha2trbGKpWowbK2tkaHDh3QoUMHDBkyRGPdvXv31MFo3759GD9+PKysLObgdINU3XPnACs8ftwSK1acBrBcY02LFi0qtS4tLcW6deuwbt06PP/883j77bfRv3//Gv1+JbJEFvObyNnZGVKpFPn5+RrL8/Pz4erqWuW2ixcvRmxsLL777jv06NFDvTw5ORm3bt1C69atYW1tDWtra1y7dg0zZ86EXC43xjCISAcnJyc89dRTmDBhAiZMmMAgZAL6n9Wq/EyWu3fvVrnF/v37MWDAAAQGBmLbtm282JrqNYv5bSSTyeDv74/ExET1MqVSicTERI3TZk9auHAh5s+fj4MHDyIgIEBj3bhx4/Dbb78hLS1N/XJ3d8esWbPw7bffGm0sRESWQN/nzgG1vxYoJSUFY8aMQbt27RAXF4eioqJa90VkLhZ1miwqKgqRkZEICAhAUFAQ4uLi8PDhQ/V56/Hjx8PDwwMxMTEAgAULFuC9997Dli1bIJfLkZeXBwBo2rQpmjZtihYtWlQ61GtjYwNXV1d06tTJtIMjIjKxiufO5eSo7iF7kkQCuLsrsXfvMly9ehlXrlzB5cuqr8eOHatRsLl+/TreeustvP/++3j99dcxdepUeDw5pTmRhbKoMDRq1Cjcvn0b7733HvLy8tCzZ08cPHgQLi4uAFT/2P56aH3lypUoLS3FyJEjNfqJjo7G+++/b8rSiYgsTsVz50aOVAWfvwaiist8li+3gp+fL/z8fNXrMjMz4evri5qwtbVFx44d0aVLF9jZ2eH48eMYNmwYpE/O3ElkgSzmbjJLxsdxEFF9lpAATJ+ueTG1lxcQF6d6Lt1fKRQKhISE6JzSxMnJCV26dEHnzp01vsrlcgYfsjiiu5uMiIi0i4gAhg2r+rlzFZYsWaKe0b9Lly6Vgk+rVq149xg1ODwypAceGSIisTh79izkcjmaNm1q7lKI6oxHhoiIqMa6d+9u7hKITM5ibq0nIiIiMgeGISIiIhI1hiEiIiISNYYhIiIiEjWGISIiIhI1hiEiIiISNYYhIiIiEjWGISIiIhI1hiEiIiISNYYhIiIiEjWGISIiIhI1hiEiIiISNYYhIiIiEjWGISIiIhI1hiEiIiISNYYhIiIiEjWGISIiIhI1hiEiIiISNYYhIiIiEjWGISIiIhI1hiEiIiISNYYhIiIiEjWGISIiIhI1hiEiIiISNYYhIiIiEjWGISIiIhI1hiEiIiISNYYhIiIiEjWGISIiIhI1hiEiIiISNYYhIiIiEjVrcxdARET1n0IBJCcDN28Cbm5ASAgglZq7KiL9MAwREVGdJCQA06cDN278uczTE1i2DIiIMF9dRPriaTIiIqq1hARg5EjNIAQAOTmq5QkJ5qmLqCYYhojqEYUCSEoCtm5VfVUozF0RiZlCoToiJAiV11UsmzGDf0/J8jEMEdUTCQmAXA48+yzw6quqr3I5/+dN5pOcXPmI0F8JApCdrWpHZMksLgytWLECcrkcdnZ2CA4OxokTJ3S2XbNmDUJCQuDk5AQnJyeEhYVptC8rK8Ps2bPh4+ODJk2awN3dHePHj0dubq4phkJkMDwVQZbo5k3DtiMyF4sKQ9u3b0dUVBSio6ORmpoKX19fhIeH49atW1rbJyUlYcyYMTh8+DCOHTsGLy8vDBw4EDk5OQCA4uJipKamYt68eUhNTUVCQgIyMjIwdOhQUw6LqE54KoIslZtb7dqVlZVBqVQaviCiWpIIgrZfseYRHByMwMBAxMfHAwCUSiW8vLwwdepUzJkzp9rtFQoFnJycEB8fj/Hjx2tt8+uvvyIoKAjXrl1D69at9aqrsLAQjo6OKCgogIODg/4DIjKApCTVKbHqHD4M9O9v7GqI/qRQqE7V5uRoD+uAEsANbN78M8aOHa1eun//fty7dw9jx441UaUkVvp+flvMkaHS0lKkpKQgLCxMvczKygphYWE4duyYXn0UFxejrKwMzZs319mmoKAAEokEzZo109mmpKQEhYWFGi8ic6nLqQhBEFBaWmrYgoj+IJWqbp8HAInkybUVR35mYNy4V7Fu3Tr1mq+++gpz5sxBcXGxKcokqpbFhKE7d+5AoVDAxcVFY7mLiwvy8vL06mP27Nlwd3fXCFR/9fjxY8yePRtjxoypMiHGxMTA0dFR/fLy8tJ/IEQGpu+piB07luHRo0cay27duoXp06cboSoilYgIYOdOwMNDc7m1dR6AkQB2QxAETJo0CfHx8VAqldi7dy9u3LiBpUuXmqNkokosJgzVVWxsLLZt24bdu3fDzs6u0vqysjK88sorEAQBK1eurLKvuXPnoqCgQP3Kzs42VtlE1QoJUU1gV/l/3hWUAK5jz54oBAQEIC0tTb0mIyMDq1atwoYNG4xfKIlWRARw9arqVO2WLaqvWVkSdOt2UaPd1KlTMWHCBPV/cGNjY3lDC1kEiwlDzs7OkEqlyM/P11ien58PV1fXKrddvHgxYmNj8d1336FHjx6V1lcEoWvXruHQoUPVXvdja2sLBwcHjReRueh7KgJQIj09HUFBQVi8eDGUSiUyMjIAAG+88QZOnTplmoJJlKRS1TVrY8aovnp6uiEpKQl+fn4a7TZt2qT+/uHDh3j33XdNWyiRFhYThmQyGfz9/ZGYmKheplQqkZiYiN69e+vcbuHChZg/fz4OHjyIgICASusrglBmZia+//57tGjRwij1ExmTrlMRnp7A4MHrAOxWLysrK8OsWbMQFhaGpKQkAKpTxBEREbh7967piibRc3Z2rvZ3+IYNG5CammrCqoi0ECzItm3bBFtbW2HDhg1Cenq6MHnyZKFZs2ZCXl6eIAiCMG7cOGHOnDnq9rGxsYJMJhN27twp3Lx5U/0qKioSBEEQSktLhaFDhwqenp5CWlqaRpuSkhK96yooKBAACAUFBYYdMFENlZcLwuHDgrBli+preblq+YEDBwQXFxcBQJWv8PBwobxiIyITePjwobB69WrByspK59/Lfv36CUql0tylUgOk7+e3RYUhQRCETz/9VGjdurUgk8mEoKAg4ZdfflGv69evnxAZGan+uU2bNlr/YUVHRwuCIAhZWVk6//EdPnxY75oYhqg+uHXrljBs2LBqA9G8efPMXSqJwLFjx4Rx48YJTZo0qfbvJAAhISHB3CVTA6Tv57dFzTNkqTjPENUXgiDgf//7H6ZNm4bHjx/rbPf1119jyJAhJqyMxKawsBDLly/HkiVLcP/+/Wrbt2vXDunp6ZDJZMYvjkSj3s0zRER1d/HiRRw8eLDKIAQAr732Gi5dumSiqkiMHBwc8O677yIrKwvz5s2Dvb19le0vX76snnCXyNR4ZEgPPDJElu7u3buYO3cu1q1bB4Wez+Xw8fHBL7/8gsaNGxu5OiLV39FFixbh008/1TnZoqOjIy5dugRnZ2cTV0cNFY8MEYlI8+bN8frrr+Ptt9+Gt7e3XtucOXMGkydPBv8/RKbQokULxMbG4sqVK3jrrbe0zgdXUFCA999/3/TFkejxyJAeeGSI6hNBEJCamoodO3Zgx44duHr1apXtly9fjqlTp5qmOKI/5Obm4uOPP8bq1atRVlamXi6VSvHbb7+ha9euZqyOGgp9P78ZhvTAMET1lSAIOHnyJHbs2IEvv/wS165dq9TG2toaSUlJ6Nu3rxkqJLG7fv06/u///g/r169HeXk5AGDw4MHYv3+/mSujhoBhyIAYhqghEAQBv/76q/qI0V8fM+Pm5obU1NRqZ3snMpYrV67gww8/xKZNm6BUKnHw4EGEh4ebuyyq5xiGDIhhiBoaQRBw/PhxfPnll9ixYwdu3LiBkJAQJCYmwsbGxtzlkYhlZGTggw8+wLlz55CSkgJra2tzl0T1GMOQATEMUUOmVCpx/Phx7NixA56enpg5c6a5SyLC2bNnIZPJ0LFjR3OXQvUYw5ABMQyRWCiVSlhZ8SZTImoYeGs9EdUYgxARiRF/8xEREZGoMQwRERGRqDEMERERkagxDBEREZGoMQwRERGRqDEMERERkagxDBEREZGoMQwRERGRqDEMERERkagxDBEREZGoMQwRERGRqDEMERERkagxDBEREZGoMQwRERGRRSssLMTixYuRkpIChUJh8P4ZhoiIiMgsFAogKQnYulX1VVfOcXBwwNmzZxEQEICWLVsiIiIC8fHxOH/+PARBqHMdEsEQvTRwhYWFcHR0REFBARwcHMxdDhERUb2XkABMnw7cuPHnMk9PYNkyICKicvvs7Gx07NgRjx8/1lju5uaGAQMGqF9yuVy9Tt/Pb4YhPTAMERER6VZSUoItW7YgMDAQXbt2hZVV1SeeEhKAkSOBJxOIRKL6unOn9kA0d+5cxMbGVtm3t7c3BgwYgNDQUAQEBKBDhw4MQ4bAMERERFS16dOnY/ny5XB0dERwcDD69OmD3r17Izg4GI6Ojup2CgUgl2seEdKkhEx2G126PI9Hjx7g0aNHKC4uRnFxMR49elSr2hiGDIBhiIiIxEihAJKTgZs3ATc3ICQEkEq1t71//z46duyI27dvayyXSCTo1q2bOhxJpaEYP95Lj3fvD+BIXYcAgGHIIBiGiIhIbGp6TQ8ArFu3DpMmTaqm59EAtupRwRgA2/QrVodmzZrh/v371X5+W9fpXYiIiKjB0XVNz40bAkaMAEaM2IpWrY7i999/V7/u3r2Lu3fv6tH7Tb1qcHR8BEfH1mjcuDEaN26MRo0aoby8HMePH69yO3t7ewwfPhyjRo1CUFAQWrZsWe171ToMlZWVIS8vD8XFxWjZsiWaN29e266IiIh0qsmpGqo7hUJ1REj7eSMJACV27XoawDgAylq8QzKAbAAe0DbDj0QiwNNTgqysPZX2c0xMjNYw1LhxYwwdOhSjRo3CoEGDYGdnB0B1ZkcfNQpDRUVF2Lx5M7Zt24YTJ06gtLQUgiBAIpHA09MTAwcOxOTJkxEYGFiTbomIiLSqzakaqpvk5KoubgZUAaY1gBDU7poeJYDpAHb+8f2fgUh1N5kEcXGVA29JSQmWL1+u/tnW1hbPP/88Ro8ejRdeeAFNmjSpRS0qek+6uHTpUsjlcqxfvx5hYWHYs2cP0tLScPHiRRw7dgzR0dEoLy/HwIEDMWjQIGRmZta6KCIioopTNU9+MOfkqJYnJJinrobupn5nsQC4qb9zcnJCu3btEBgYCFtb2yq36tKlCxYu7I3//a8Anp6aMcTTU/dt9V988QXu3r2LF154AZs2bcKtW7eQkJCAV155pU5BCKjBBdRjxozBu+++i27dulXZrqSkBOvXr4dMJsPf/va3OhVnKXgBNRGRaVV3+7VEovrgzMriKTNDS0oCnn22+nabNt3A4MGN0KxZM0j/2AlXrlxB+/btK80K7ejoiNGjR2PixIkICgqC5I8JhWpyCvTHH3+Ej48PnJyc9B4LJ100IIYhIiLT0vcD+fBhoH9/Y1cjLhVBNCdH+3VDVQXR2bNnY+HChX+0kyA0NBQTJ07E8OHD0ahRI+MX/wR9P795NxkREVkcfU/V6H9Kh/QllaquyRo5UhV8/hqIKmaI1nZNz+PHj7F27Vp4e3tjwoQJiIyMROvWrU1Wd10Y9EGt2dnZDebUGBERmY+bW/VtatKOgMzMTL3vroqIUF274+Ghubyqa3pyc3Oxa9cuZGZmYt68efUmCAEGDkO///47Nm7caMguiYhIhEJCVB+8FUciKlMCuI4LF9ZUWnP//n0jVmZ59H3y+927d9GiRQv07dsX77//Po4ePYqysjKd/UZEAFevqk5Fbtmi+pqVpfsuPm9vb/Tr16/a55JZohpdM7R3794q11+5cgUzZ86EQteeqKd4zRARkelV3E0GPHntSsXcNiMhle7F3r178fzzzwMASktL8fzzz+PQoUPqi3QbsppOPRAREYHdu3erf7a3t0f//v0RFhaG5557Dp07d25Qf256f34LNSCRSAQrKytBIpHofFlZWdWky0ri4+OFNm3aCLa2tkJQUJBw/PhxnW1Xr14tPP3000KzZs2EZs2aCaGhoZXaK5VKYd68eYKrq6tgZ2cnhIaGChcvXqxRTQUFBQIAoaCgoFZjIiKi2tm1SxA8PQVBFYdUL5ksTwCGCwAEAEKTJk2ElJQUQRAEYdu2bQIAYf/+/Wau3Ph27RIEiUTzzwZQLZNIVOufdP78ecHKykr9Z/fky8PDQ4iMjBQ2b94s3Lx50/SDMjB9P79rdGTIw8MD//3vfzFs2DCt69PS0uDv71/rI0Pbt2/H+PHjsWrVKgQHByMuLg5ffvklMjIy0KpVq0rtx44di759+6JPnz6ws7PDggULsHv3bpw7dw4ef5zoXLBgAWJiYrBx40a0bdsW8+bNw5kzZ5Cenq6eobI6PDJERGQ+T95+3bXrXTz9dG+N+excXV1x/PhxjB8/HkeOHMGzzz6LH374wYxVG1f1T34X4OBQgEmTPkJZ2WM8fvwYJSUlePz4MQ4dOqT3qUQfHx8899xzCAsLwzPPPFPn+XxMzShHhl588UVh3rx5OtenpaUJEomkJl1qCAoKEqZMmaL+WaFQCO7u7kJMTIxe25eXlwv29vbCxo0bBUFQHRVydXUVFi1apG5z//59wdbWVti6davOfh4/fiwUFBSoX9nZ2TwyRERkQS5duiQ4OztrHNVo27atxs8nT540d5lGc/hw5SNC2l/9dB4FqulLJpMJ/fv3F7744gtBqVSa+49AL/oeGarRVU6zZs1Cnz59dK5v3749Dh8+XJMu1UpLS5GSkoKwsDD1MisrK4SFheHYsWN69VFcXIyysjL1c9KysrKQl5en0aejoyOCg4Or7DMmJgaOjo7ql5eXV63GRERExtGuXTt8/fXXGkf4s7KyNNosWbLE1GWZTG1mia4LuVyO6dOnY+HChRgzZkyDuq4IqOHdZCEhIRg0aJDO9U2aNEG/fv1qVcidO3egUCjg4uKisdzFxQV5eXl69TF79my4u7urw0/FdjXtc+7cuSgoKFC/srOzazIUIiIygaeeegqff/65zvU7duzAtWvXTFiR6eg/pUDtJ2KSy+V455138Ouvv+LKlStYuHAhAgMDG1wQAhrQpIuxsbHYtm0bkpKS9L4WSBdbW9tqn61CRETmcfbsWSQmJuKHH37AkSO6HxSqUCiwbNkyLF261ITVmUbF1AO6ZokGBLi4lOLAgeVo0sQOtra2sLNTfQ0JCcHZs2e19tu2bVu8/PLLePnll+Hv798gg482FhOGnJ2dIZVKkZ+fr7E8Pz8frq6uVW67ePFixMbG4vvvv0ePHj3Uyyu2y8/Ph9tfYnR+fj569uxpuOKJiMhkbt26hfj4eFy6dKnatmvWrMF7772HZs2aGb8wE6p+lmgJ/vtfW/Tq1UNju++++65SEPL29lYHID8/P9EEoL+ymJmRZDIZ/P39kZiYqF6mVCqRmJiI3r1769xu4cKFmD9/Pg4ePIiAgACNdW3btoWrq6tGn4WFhTh+/HiVfRIRkeUaMGAAfvvtN/znP/+BtXXV/6d/8OABVq9ebaLKTKs2s0RXPDfM29sbs2fPRkpKCi5duoTY2FhRHQmqpLZXaGdnZwsKhaLS93Wxbds2wdbWVtiwYYOQnp4uTJ48WWjWrJmQl5cnCIIgjBs3TpgzZ466fWxsrCCTyYSdO3cKN2/eVL+Kioo02jRr1kz46quvhN9++00YNmyY0LZtW+HRo0d618V5hoiILNPZs2eFPn36VHkXlLu7u1BSUmLuUo2mvFx1d9mWLaqv5eXa22VnZwtz5swRUlJS6s3dYHWl7+d3rcOQvb29cPny5Urf19Wnn34qtG7dWpDJZEJQUJDwyy+/qNf169dPiIyMVP/cpk0brX/xo6Oj1W0qJl10cXERbG1thdDQUCEjI6NGNTEMERFZLoVCIaxcuVJwcHDQGYgqplwhcTHKpIt/ZW9vj9OnT8Pb21vj+4aIky4SEVm+3NxcTJ8+HTt37qy0zsfHB6dPnxbvaSCR0vfz22KuGSIiIqoLd3d3fPnll/j6668rzQ935swZHDp0yCx16fsgVTIfhiEiImpQhgwZgvT0dLz11lsaT1BfvHixyWtJSFA9NuPZZ4FXX1V9lctVy7XJz8/HunXrGtwDzy0dwxARETU4TZs2xdKlS3HixAn06tULAHDo0CGcPn26UltjHblJSFDd+v7k88NyclTLtQWiVq1a4d///jf8/f0b9LPVLA3DEBERNVj+/v44ceIEFi9ejMaNG1d6REdNj9zoS6EApk/XPiGi6qlhAqZMKUVZmVJjnUQiwTPPPIPTp08jNDQUw4YNw8WLF+tWDFWLYYiIiBo0a2trzJw5E+fOnUNhYSFu/HGopjZHbvSVnFzVE+UBQIK8PBmcnIZi8ODB+PDDD/H999+jqKhI47FWe/fuRbdu3TBjxgz8/vvvtS+IqsS7yfTAu8mIiBoGQRBQVFSEJk0cIJfrDiwSiWrywqws1WzPNbV1q+pIU/XGANim/snKygru7u7qwPZXTk5OiI6OxhtvvAGZTFbzokTI6HeT/fvf/1Y/Hf6v3xMREVkqiUQCBweHao/cCAKQna06wlMbtX2QqlKp1BqEAODevXuYMWMGunfvjr1796KWxzJIi1qHoblz56qf9fLX74mIiCzdTT0f5q5vuydVPEhV97RGAho3voM2bbJr3HdmZiaGDRuG0NBQpKWl1a5A0sBrhoiISHT0PXIjkeRVWrZlyxbcvXu3yu0qHqSq6uPJPlVHqDZtcsbVq5eRk5ODnTt3IioqCoGBgfoVBuDw4cPw8/PDpEmTcLO2qY0AMAwREZEIVX/kRgngOiZP7oLt27drrDlw4ADGjRsHpVKpfdM/6PsgVXd3d4wYMQJLlixB3759azQOQRCwbt06dOjQAfPnz0dxcXGNticVhiEiIhKdqo7cqIIQAMxAUdF9jB49GpMmTcLDhw8BAOfOncOBAwcQGxtb7ftERABXrwKHDwNbtqi+ZmVpf6L8+vXrERcXV6vxlJWV4b///S8iIiKQkZFRqz7ErNZ3kwHA7du3cfv2bXTt2lVjeXp6Olq2bImWLVvWuUBLwLvJiIgapoQE1XxAf71muXnzhygunozHj7dotO3UqRO++OILPP3003j8+DGsrKzw/fff49lnn61zHT///DMGDBgAmUyG5s2ba7ycnJyqXdaoUSM+d00LfT+/6xSGRo0ahWnTplU6rPfTTz8hPj4eW7durW3XFoVhiIio4VIoVHeN3bypupYoJAS4evUyxowZg19//VWjrVQq1XhUhouLC06dOgU3/W8f0+rBgwewtbWFjY1NnfohTSZ5UGtWVpbW85t9+/bF2bNn69I1ERGRSUilQP/+wJgxqq9SKdCuXTscPXoU77zzjkbbJ58Zlp+fj9GjR6O8vLxONTRt2pRByIzqFIbu3bunc92jR4/q0jUREZFZyWQyLFiwAN999x1cXFx0tvvxxx8xb948E1ZGhlanMNSjRw9s2LCh0vLPP/8cPj4+demaiIjIIjg5OaFDhw5VtomNjcU333xjoorI0Op0zVBOTg6GDRsGe3t7+Pn5AQBSU1NRVFSEPXv2wNPT02CFmhOvGSIiEp9z585h3rx52L17t17tnZyckJqaCrlcbtzCSG/6fn5b1+VNPDw8cPLkSSQmJiI9PR0AMHjwYISFhdWlWyIiIrMrKCiAm5sbnJ2dcefOnWrb37t3D6+88gqSk5Nha2trggrJUOp0ZEgseGSIiEi8ysrKkJiYiC1btmD37t148OBBle2nTJmC+Ph4E1VHVTHJrfXjxo3Dpk2bEBgYqDG/gSAIkEgkOHHiRG27tigMQ0REBADFxcXYt28ftm7din379qG0tFRru23btmHUqFEmro6eZJIwdP36dTx48ABNmjTRWP7gwQM8fPgQQUFBte3aojAMERHRk+7fv4/du3djy5Yt+OGHHzQez9G0aVOcPHkSnTp1MmOFZJJ5hmbNmoV79+6hTZs2Gq/79+/jk08+qUvXREREFq1Zs2aYOHEiDh06hJycHCxbtgxPPfUUANVBgREjRqgf4UGWrU5HhoKCgnSeCvPx8cGZM2dqXZgl4ZEhqglts9lKpeauiohM5cqVK9i2bRu2bNkCf39/bNiwgY/KMBOTnCbr0KEDMjMzta5r3749Ll26VNuuLQrDEOlL23OOPD1VD4TU9mBGImrYzpw5Aw8PDzRv3tzcpYiSSU6TcdJFoj8lJAAjR2oGIQDIyVEtT0gwT11EZD4+Pj4MQvUAJ13UA48MUXUUCkAurxyEKkgkqiNEWVk8ZUZEZCqcdJHIhJKTdQchABAEIDtb1a5/f5OVRUREeqhTGKoQGhqK0NBQQ3RFVC/dvFn7dgUFBXB0dDRsQUREpLcaXTPk4eGBoUOH4oMPPsDXX3+NnJwcY9VFVK+4uenXztb290rLxo0bh2vXrhm4IiIi0leNrhn69NNPkZqaitTUVJw/fx4KhQItW7aEn58f/P394efnBz8/P7Rp08aYNZscrxmi6lRcM5STozolVpkSwA00a+aHpUsXYcKECepbbd3d3dG8eXP89NNPPEJERGRARr+1vqSkBGlpaepwlJqaivT0dJSVlaG8vLzWhVsihiHSR8XdZMCTgahiVtqRAFRPvw4LC8Nnn30GLy8v2NraQhAEhIeH45tvvoG1tUHOXhMRiZ7Rb623tbVFcHAw/v73v2PYsGHw8fFBo0aNKj2ag0gsIiKAnTsBDw/N5Y6ORbC2HoOKIAQA33//PXx8fDBv3jxU/H/k22+/xbRp08BnJxMRmVatwtDjx4+xe/dujB07Fi1btsTEiRMhlUqxadMm3L5929A1EtUbERHA1avA4cPAli2qr3fvOuLs2Q8REhKi0ba4uBgLFizQWLZy5UosX77chBUTEVGNTpNt374du3btwoEDB2Bvb4/hw4cjIiIC/fv3h7QBT57C02RkCEqlEqtXr8Y777yDoqIine0kEgm++uorvPjiiyasjoio4THKNUNWVlZwd3fHu+++i7///e+iubaBYYgMKTs7G//85z9x4MABnW2aNGmCo0ePomfPnqYrjIiogTHKNUMhISEoKirCv/71Lzg6OqJ3796YMmUK1q1bh7S0tAZ34TSRoQmCgLS0NGRkZFTZ7uHDhxgyZAhyc3NNVBkRkXjV6m6yzMxMpKSkaNxJdv/+fdja2sLHx0fnk+zrKx4ZIkPIzMzE9OnTqzwi9CQ/Pz/8+OOPvDGBiKgWjPo4jg4dOqBDhw4YPXq0ellWVhZOnjyJU6dO1aZLogavcePGiIiIgJ2dHQ4dOoQHDx5Uu01qaipeffVVJCQkNOjr8oiIzEnvI0PXr19H69at9e44JycHHk/eY1xP8cgQGVpJSQmOHj2K/fv3Y//+/bhw4UKV7WfOnInFixebqDoioobB4NcMBQYG4p///Cd+/fVXnW0KCgqwZs0adO/eHbt27apZxX9YsWIF5HI57OzsEBwcXOUpt3PnzmHEiBGQy+WQSCSIi4ur1EahUGDevHlo27YtGjVqhHbt2mH+/Pmcy4XMytbWFqGhoViyZAnOnz+Py5cv49NPP8XgwYNhZ2dXqf2SJUvw2WefmaFSIqKGT+/TZOnp6fjoo4/w3HPPwc7ODv7+/nB3d4ednR3u3buH9PR0nDt3Dn5+fli4cCGef/75Ghezfft2REVFYdWqVQgODkZcXBzCw8ORkZGBVq1aVWpfXFwMb29vvPzyy3jrrbe09rlgwQKsXLkSGzduRLdu3XDy5ElMnDgRjo6OmDZtWo1rJDIGb29vvPnmm3jzzTdRXFyMpKQk7Nu3D/v27VM/t2zKlCnw9vbGc889Z+ZqiYgalhpfQP3o0SPs378fycnJuHbtGh49egRnZ2f06tUL4eHh6N69e62LCQ4ORmBgIOLj4wGo5mXx8vLC1KlTMWfOnCq3lcvlmDFjBmbMmKGxfMiQIXBxccHatWvVy0aMGIFGjRph8+bNetXF02RkLoIg4MKFC9i/fz/27duHM2fO4MiRI+jatau5SyMisnhGu4D61q1bkMlkePXVVxEUFFSnIv+qtLQUKSkpmDt3rnqZlZUVwsLCcOzYsVr326dPH6xevRoXL15Ex44dcfr0aRw9ehRLly7VuU1JSQlKSkrUPxcWFtb6/YnqQiKRoEuXLujSpQtmzpyJwsJCXLp0ydxlERE1KDUKQ1u3bsWECRNQVlYGiUSCXr164cCBA2jZsmWdC7lz5w4UCgVcXFw0lru4uFR7cWlV5syZg8LCQnTu3BlSqRQKhQIfffQRxo4dq3ObmJgYfPDBB7V+TyJjcXBwgJ+fn7nLICJqUGo06eIHH3yAV199FRcuXMB3330HANWevjK3HTt24IsvvsCWLVuQmpqKjRs3YvHixdi4caPObebOnYuCggL1Kzs724QVExERkSnV6MjQlStXcPDgQcjlcnTs2BGbN2+Gv7+/xvU4teXs7AypVIr8/HyN5fn5+XB1da11v7NmzcKcOXPUcyL5+Pjg2rVriImJQWRkpNZtbG1tYWtrW+v3JCIiovqjRkeGysvL0bhxY/XPnTt3hlKpRF5eXp0Lkclk8Pf3R2JionqZUqlEYmIievfuXet+i4uLYWWlOUypVAqlUlnrPomIiKjhqPEF1Bs3bkTfvn3Ro0cPNG3aFNbW1iguLjZIMVFRUYiMjERAQACCgoIQFxeHhw8fYuLEiQCA8ePHw8PDAzExMQBUF12np6erv8/JyUFaWhqaNm2K9u3bAwBefPFFfPTRR2jdujW6deuGU6dOYenSpfjb3/5mkJqJiIiofqvRrfX9+vVDWloaioqKYGVlhbZt2+Lq1at45513EBYWhoCAANjb29epoPj4eCxatAh5eXno2bMnli9fjuDgYABA//79IZfLsWHDBgDA1atX0bZtW611JiUlAQCKioowb9487N69G7du3YK7uzvGjBmD9957DzKZTK+aeGs9ERFR/aPv57dBH9RqZWWFDh064Pz583Uq3tIwDBEREdU/fFArERERkR5qdWRIbHhkiIiIqP4x+INaiYiIiBoihiEiIiISNYYhIiIiEjWGISIiIhI1hiEiIiISNYYhIiIiEjWGISIiIhI1hiEiIiISNYYhIiIiEjWGISIiIhI1hiEiIiISNYYhIiIiEjWGISIiIhI1hiEiIiISNYYhIiIiEjWGISIiIhI1hiEiIiISNYYhIiIiEjWGISIiIhI1hiEiIiISNYYhIiIiEjWGISIiIhI1hiEiIiISNYYhIiIiEjWGISIiIhI1hiEiIiISNYYhIiIiEjWGISIiIhI1hiEiIiISNYYhIiIiEjWGISIiIhI1hiEiIiISNYYhIiIiEjWGISIiIhI1hiEiIiISNYYhIiIiEjWGISIiIhI1hiEiIiISNYsLQytWrIBcLoednR2Cg4Nx4sQJnW3PnTuHESNGQC6XQyKRIC4uTmu7nJwcvPbaa2jRogUaNWoEHx8fnDx50kgjICIiovrEosLQ9u3bERUVhejoaKSmpsLX1xfh4eG4deuW1vbFxcXw9vZGbGwsXF1dtba5d+8e+vbtCxsbGxw4cADp6elYsmQJnJycjDkUIiIiqickgiAI5i6iQnBwMAIDAxEfHw8AUCqV8PLywtSpUzFnzpwqt5XL5ZgxYwZmzJihsXzOnDn46aefkJycrHcdJSUlKCkpUf9cWFgILy8vFBQUwMHBQf8BERERkdkUFhbC0dGx2s9vizkyVFpaipSUFISFhamXWVlZISwsDMeOHat1v3v37kVAQABefvlltGrVCr169cKaNWuq3CYmJgaOjo7ql5eXV63fn4iIiCybxYShO3fuQKFQwMXFRWO5i4sL8vLyat3vlStXsHLlSnTo0AHffvst3njjDUybNg0bN27Uuc3cuXNRUFCgfmVnZ9f6/YmIiMiyWZu7AGNTKpUICAjAxx9/DADo1asXzp49i1WrViEyMlLrNra2trC1tTVlmURERGQmFnNkyNnZGVKpFPn5+RrL8/PzdV4crQ83Nzd07dpVY1mXLl1w/fr1WvdJREREDYfFhCGZTAZ/f38kJiaqlymVSiQmJqJ379617rdv377IyMjQWHbx4kW0adOm1n0SERFRw2FRp8mioqIQGRmJgIAABAUFIS4uDg8fPsTEiRMBAOPHj4eHhwdiYmIAqC66Tk9PV3+fk5ODtLQ0NG3aFO3btwcAvPXWW+jTpw8+/vhjvPLKKzhx4gRWr16N1atXm2eQREREZFEs6tZ6AIiPj8eiRYuQl5eHnj17Yvny5QgODgYA9O/fH3K5HBs2bAAAXL16FW3btq3UR79+/ZCUlKT++ZtvvsHcuXORmZmJtm3bIioqCv/4xz/0rknfW/OIiIjIcuj7+W1xYcgSMQw1LAoFkJwM3LwJuLkBISGAVGruqoiIyND0/fy2qNNkRMaWkABMnw7cuPHnMk9PYNkyICLCfHUREZH5WMwF1ETGlpAAjBypGYQAICdHtTwhwTx1ERGReTEMkSgoFKojQtpOClcsmzFD1e5Jv/32G5RKpVHrIyIi82EYIlFITq58ROivBAHIzla1e9Lu3bsxfPhwFBYWGq9AIiIyG4YhEoWbN/Vrt2jRZly8eFFjmbOzM/bu3Yunnnqq0joiIqr/GIZIFNzc9Gu3f///0KlTJ4SFhWHnzp0oKyuDs7MzAOD8+fMICgrC/v37jVgpERGZGsMQiUJIiOquMYlEVwslgOsAVOfJEhMT8fLLL6N169bYsmWLulVBQQGGDBmCmJgYcFYKIqKGgWGIREEqVd0+D2gLRAIACaTSt6EKRX/Ky8vD3r17NVsLAv7973/jlVdewYMHD4xVMhERmQjDEIlGRASwcyfg4aG53MtLgl27JMjNjUdMTAzkcrle/e3cuRN9+vTBlStXDF8sERGZDGeg1gNnoG5YqpuBWqlU4rvvvsOqVavw9ddfV3tbvZOTE7Zv347nnnvOyJUTEVFN8HEcBsQwJE4//PCD3rfUW1lZYeHChYiKioJE94VJRERkQvp+fvM0GZEWmzdvxqBBg/SeW0ipVOLtt9/Ga6+9huLiYiNXR0REhsQwRPSEr7/+GrNnz4ZC23TU1diyZQuefvppXLt2zQiVERGRMfA0mR54mkycFAoF7ty5g7y8POTl5eHmzZsaX//6fVFRkca2zs7O+PLLL9G/f3/zFE9ERHxqPVFdSaVSuLi4wMXFBb6+vlW2ffjwoTogVYSkpKQktG/fHp6eniaqmIiIaoNhiMgAmjRpgnbt2qFdu3bmLoWIiGqI1wwRERGRqDEMERERkagxDBEREZGoMQwRERGRqDEMERERkagxDBEREZGoMQwRERGRqDEMERERkagxDBEREZGoMQwRERGRqDEMERERkagxDBEREZGoMQwRERGRqDEMERERkagxDBEREZGoMQwRERGRqDEMERERkagxDBEREZGoMQwRERGRqDEMERERkagxDBEREZGoMQwRERGRqDEMERERkahZXBhasWIF5HI57OzsEBwcjBMnTuhse+7cOYwYMQJyuRwSiQRxcXFV9h0bGwuJRIIZM2YYtmgiIiKqtywqDG3fvh1RUVGIjo5GamoqfH19ER4ejlu3bmltX1xcDG9vb8TGxsLV1bXKvn/99Vd89tln6NGjhzFKJyIionrKosLQ0qVL8Y9//AMTJ05E165dsWrVKjRu3Bjr1q3T2j4wMBCLFi3C6NGjYWtrq7PfBw8eYOzYsVizZg2cnJyMVT4RERHVQxYThkpLS5GSkoKwsDD1MisrK4SFheHYsWN16nvKlCl44YUXNPquSklJCQoLCzVeRERE1DBZTBi6c+cOFAoFXFxcNJa7uLggLy+v1v1u27YNqampiImJ0XubmJgYODo6ql9eXl61fn8iIiKybBYThowhOzsb06dPxxdffAE7Ozu9t5s7dy4KCgrUr+zsbCNWSUREROZkbe4CKjg7O0MqlSI/P19jeX5+frUXR+uSkpKCW7duwc/PT71MoVDgxx9/RHx8PEpKSiCVSittZ2trW+U1SERERNRwWMyRIZlMBn9/fyQmJqqXKZVKJCYmonfv3rXqMzQ0FGfOnEFaWpr6FRAQgLFjxyItLU1rECIiIiJxsZgjQwAQFRWFyMhIBAQEICgoCHFxcXj48CEmTpwIABg/fjw8PDzU1/+UlpYiPT1d/X1OTg7S0tLQtGlTtG/fHvb29ujevbvGezRp0gQtWrSotJyIiIjEyaLC0KhRo3D79m289957yMvLQ8+ePXHw4EH1RdXXr1+HldWfB7Nyc3PRq1cv9c+LFy/G4sWL0a9fPyQlJZm6fDIQhQJITgZu3gTc3ICQEIAH8YiIyFgkgiAI5i7C0hUWFsLR0REFBQVwcHAwdzkNWkICMH06cOPGn8s8PYFly4CICO3b7Nq1C8899xz3DRERadD389tirhkiSkgARo7UDEIAkJOjWp6QoH271NRUeHt7Y8mSJXj06JHxCyUiogaFYYgsgkKhOiKk7TilapmA6dMFKBSV14eHh+Pu3bt4++230aFDB6xZswbl5eXGLpmIiBoInibTA0+TGV9SEvDss9W369BhMvr1E+Dn5wc/Pz/06NEDUqkULVq0wIMHD/7SrgPmz5+Pl19+WeM6MyIiEg99P78t6gJqEq+bN/Vrl5lZhMzMbeqfpVIpunTpUmmahMzMTIwePRqxsbH46KOPMHjwYEgkEkOWTEREDQT/y0wWwc1N35aaqUmhUODs2bMoKCjQ2jotLQ0vvPACnnnmGRw9erRuRRIRUYPEMEQWISREddeY7oM3Amxs8mBjc7xW/R89ehQhISF4/vnnkZaWVtsyiYioAWIYIosglapunwcqByKJBJBIJNi2zRUPHhTg1KlTWLt2LaZMmYLevXujUaNGer/PgQMH0KtXL4wePRqZmZkGHAEREdVXvIBaD7yA2nS0zTPk5QXExWmfZ6i4uBi9e/fGb7/9VuP3kkql+Nvf/ob33nsPnp6etS+aiIgsEi+gpnopIgIYNky/GagFQcC//vWvGgUhBwcHuLm5wc3NDe7u7nBwcMDevXsxYcIENG7c2IAjISKi+oJhiCyOVAr07199u//973/YuHEjAMDJyQnu7u7qoFMRdv76s5ubG5o0aWLc4omIqN5hGKJ6y9/fH1lZWXB1dYWdnZ25yyEionqKYYjqLT8/P3OXQEREDQDvJiMiIiJRYxgiIiIiUWMYIiIiIlFjGCIiIiJRYxgiIiIiUWMYIiIiIlFjGCKjKisrq9WjMoiIiEyFYYhqRaEAkpKArVtVXxUK7e1sbGwwffp0+Pn5Yfny5bhz544pyyQiIqoWwxDVWEICIJcDzz4LvPqq6qtcrlquzZQpU3Dq1ClMnz4d7u7uiIiIwNdff42ysjJTlk1ERKQVn1qvBz61/k8JCcDIkcCTf2skEtWCL78ERoyQaKwrKyuDXC5Hbm6uxvJWrVrhtddew4QJE+Dj42PUuomISHz0/fxmGNIDw5CKQqE6AnTjhq4WSgA5aNEiAC1aNIOTkxOaN28OJycnpKSkICMjQ2ff/v7+mDBhAsaMGYMWLVoYoXoiIhIbhiEDYhhSSUpSnRKrXn8AR2r1HjY2Nhg6dCgmTJiAQYMGwdqaj88jIqLa0ffzm9cMiVx5eTkWLFiAtLQ0VJeLb97Ut1e3WtdTVlaGXbt24cUXX4SnpydmzZqFc+fO1bo/IiKi6jAM1QOff/45Tpw4AaVSqfc2+t7tZW1tjfz8fPTq1Qve3t6IiopCcnIyFFo2cNMz47Rr1wReXl5o2rSp3vVqk5+fj+XLl2Py5MnYsWNHtWGNiIioNniaTA/GOE2mUADJyaqjLW5uQEgIIJVqb7tx40ZMmDABrq6uGDJkCIYOHYrQ0FA0btxYa/uEBGD6dM1rezw9gWXLgIiIyu1zc3PRtm1blJaWqpe1atUKw4YNw/DhwzFgwADY2tqqrxnKyal8ATUASCSq98nK+nMspaWluHLlCvz9/VFcXFztn4uLiwv69OmDPn36oG/fvvDz84OtrW212xERET2J1wwZkKHDUE3DSnl5OTp37ozLly+rl9nZ2SEsLAxDhw7FkCFD4PbHYZvq7vbatOkxIiJUR4Ssra0hkaju/Hr99dfx2Wefaa3X3t4eL7zwAoYPH46yshcxblwjAJrv8Uc32Lmz8hjefPNNrFixolK/EokEPj4+6uDTp08ftG3bVl0TERFRXTAMGZAhw5DusKL6qi1MAMCGDRswceJEnf0GBgbihReGYuXKd5CfbwNAW6BQArgBoO0f3wNWVlawtraGIAh6zftja2uL7t3fw5Ur03Dv3p+nwby8gLi4yrWfPn0afn5+UCqVsLe3R+/evdVHfoKDg0V9QToRERkXw5ABGSoM6XNrukSSC0fHnrCyUu2WiqMkgiDg999/r+Yd+gFI0qOS/qjt3V5/soJE0g+dO/fHgAFdMXNmENq2bV2p1erVq6FUKtGnTx9069YNUl3nAomIiAxM389v3rdsQsnJVQUhALCCIHji/v3uqF1Y0fcurtrf7fUnJQThMK5dO45r1wbgwIFbGD9+fKWLpidPnmyA9yIiIjIehiETMv6t6fq9Qd++7eDkNATl5eXq17Vr15CVlaXX9l27dsXgwYMxaNAgPP3007Czs6tlvURERObHMGRC+t6a3rGjAxwcAtS3kguCgMePHyM9Pb2aLZMBZAPwgLZZEyru9jpy5P807lwTBAFPP/20zjBkb2+PsLAwDBo0CIMGDULr1pVPhxEREdVXDEMmFBKiCiPV3Zqenv5ZpdvsP/nkE0RFRWnt187ODqGhoRgyZAisrBrh9ddVQUjb3V5xcZVv4T9w4AB+/vlnjWW+vr4YNGgQBg8ejN69e0Mmk9VkqERERPUGw5AJSaWq2+dHjlSFE33DSnFxMRYsWKCxzMPDA0OGDMGQIUMwYMAAjTmHnJ2137qv7W4vpVKJd999F82aNcNzzz2HwYMHIzw8HO7u7nUfMBERUT3Au8n0YIp5hnTdmg4AS5Yswdtvv42goCB1AOrZs2eV8/HoO6nj3bt3ceHCBQQHB/M5YERE1KDw1noDMvcM1Pv27YO/vz9cXV0N8t5ERERiwDBkQHxqPRERUf3Dp9YTERER6YFhiIiIiETN4sLQihUrIJfLYWdnh+DgYJw4cUJn23PnzmHEiBGQy+WQSCSIi4ur1CYmJgaBgYGwt7dHq1at8NJLLyEjI8OIIyAiIqL6xKLC0Pbt2xEVFYXo6GikpqbC19cX4eHhuHXrltb2xcXF8Pb2RmxsrM6Li48cOYIpU6bgl19+waFDh1BWVoaBAwfi4cOHxhwKERER1RMWdQF1cHAwAgMDER8fD0A1B46XlxemTp2KOXPmVLmtXC7HjBkzMGPGjCrb3b59G61atcKRI0fwzDPPaG1TUlKCkpIS9c+FhYXw8vLiBdRERET1SL27gLq0tBQpKSkICwtTL7OyskJYWBiOHTtmsPcpKCgAADRv3lxnm5iYGDg6OqpfXl5eBnt/IiIisiwWE4bu3LkDhUIBFxcXjeUuLi7Iy8szyHsolUrMmDEDffv2Rffu3XW2mzt3LgoKCtSv7Oxsg7w/ERERWR5RTTk8ZcoUnD17FkePHq2yna2tLWxtbU1UFREREZmTxYQhZ2dnSKVS5OfnayzPz883yMzLb775Jr755hv8+OOP8PT0rHN/RERE1DBYTBiSyWTw9/dHYmIiXnrpJQCq01qJiYl48803a92vIAiYOnUqdu/ejaSkJLRt27ZWfQCqC7GIiIiofqj43K7uXjGLCUMAEBUVhcjISAQEBCAoKAhxcXF4+PAhJk6cCAAYP348PDw8EBMTA0B10XV6err6+5ycHKSlpaFp06Zo3749ANWpsS1btuCrr76Cvb29+vojR0dHNGrUSK+6ioqKAIAXUhMREdVDRUVFcHR01Lneom6tB4D4+HgsWrQIeXl56NmzJ5YvX47g4GAAQP/+/SGXy7FhwwYAwNWrV7Ue6enXrx+SkpIAQOeT3devX48JEyboVZNSqURubi7s7e2rfFJ8TVXcsp+dnd0gb9lv6OMDGv4YG/r4gIY/Ro6v/mvoYzTm+ARBQFFREdzd3WFlpfueMYsLQ2LS0B8A29DHBzT8MTb08QENf4wcX/3X0MdoCeOzmFvriYiIiMyBYYiIiIhEjWHIjGxtbREdHd1g5zRq6OMDGv4YG/r4gIY/Ro6v/mvoY7SE8fGaISIiIhI1HhkiIiIiUWMYIiIiIlFjGCIiIiJRYxgiIiIiUWMYMpAff/wRL774Itzd3SGRSLBnzx6N9QkJCRg4cCBatGgBiUSCtLQ0vfr98ssv0blzZ9jZ2cHHxwf79+83fPF6MsYYN2zYAIlEovGys7MzzgCqUdX4ysrKMHv2bPj4+KBJkyZwd3fH+PHjkZubW22/K1asgFwuh52dHYKDg3HixAkjjkI3Y4zv/fffr7T/OnfubOSR6Fbd39H3338fnTt3RpMmTeDk5ISwsDAcP3682n7rwz4Eajc+S9qH1Y3vr15//XVIJBLExcVV26+l7D/AOGOsT/twwoQJlWodNGhQtf0aex8yDBnIw4cP4evrixUrVuhc//TTT2PBggV69/nzzz9jzJgxmDRpEk6dOoWXXnoJL730Es6ePWuosmvEGGMEAAcHB9y8eVP9unbtmiHKrbGqxldcXIzU1FTMmzcPqampSEhIQEZGBoYOHVpln9u3b0dUVBSio6ORmpoKX19fhIeH49atW8Yahk7GGB8AdOvWTWP/HT161Bjl66W6v6MdO3ZEfHw8zpw5g6NHj0Iul2PgwIG4ffu2zj7ryz4Eajc+wHL2YXXjq7B792788ssvcHd3r7ZPS9p/gHHGCNSvfTho0CCNWrdu3VplnybZhwIZHABh9+7dWtdlZWUJAIRTp05V288rr7wivPDCCxrLgoODhX/+858GqLJuDDXG9evXC46OjgatzRCqGl+FEydOCACEa9eu6WwTFBQkTJkyRf2zQqEQ3N3dhZiYGEOVWiuGGl90dLTg6+tr2OIMRJ8xFhQUCACE77//Xmeb+rwP9Rmfpe5DXeO7ceOG4OHhIZw9e1Zo06aN8Mknn1TZj6XuP0Ew3Bjr0z6MjIwUhg0bVqN+TLEPeWTIgh07dgxhYWEay8LDw3Hs2DEzVWQcDx48QJs2beDl5YVhw4bh3Llz5i5JLwUFBZBIJGjWrJnW9aWlpUhJSdHYh1ZWVggLC6sX+7C68VXIzMyEu7s7vL29MXbsWFy/ft00BdZRaWkpVq9eDUdHR/j6+upsU1/3oT7jq1Bf9qFSqcS4ceMwa9YsdOvWrdr29XH/1XSMFerLPgSApKQktGrVCp06dcIbb7yBu3fv6mxrqn3IMGTB8vLy4OLiorHMxcUFeXl5ZqrI8Dp16oR169bhq6++wubNm6FUKtGnTx/cuHHD3KVV6fHjx5g9ezbGjBmj88GCd+7cgUKhqJf7UJ/xAUBwcDA2bNiAgwcPYuXKlcjKykJISAiKiopMWG3NfPPNN2jatCns7OzwySef4NChQ3B2dtbatj7uw5qMD6hf+3DBggWwtrbGtGnT9GpfH/dfTccI1K99OGjQIHz++edITEzEggULcOTIEQwePBgKhUJre1PtQ2uD9URUC71790bv3r3VP/fp0wddunTBZ599hvnz55uxMt3KysrwyiuvQBAErFy50tzlGFxNxjd48GD19z169EBwcDDatGmDHTt2YNKkScYutVaeffZZpKWl4c6dO1izZg1eeeUVHD9+HK1atTJ3aQZR0/HVl32YkpKCZcuWITU1FRKJxNzlGEVtx1hf9iEAjB49Wv29j48PevTogXbt2iEpKQmhoaFmq4tHhiyYq6sr8vPzNZbl5+fD1dXVTBUZn42NDXr16oVLly6ZuxStKoLCtWvXcOjQoSqPmjg7O0MqldarfViT8WnTrFkzdOzY0WL3HwA0adIE7du3x1NPPYW1a9fC2toaa9eu1dq2Pu7DmoxPG0vdh8nJybh16xZat24Na2trWFtb49q1a5g5cybkcrnWberb/qvNGLWx1H2ojbe3N5ydnXXWaqp9yDBkwXr37o3ExESNZYcOHdI4ktLQKBQKnDlzBm5ubuYupZKKoJCZmYnvv/8eLVq0qLK9TCaDv7+/xj5UKpVITEy0yH1Y0/Fp8+DBA1y+fNki958uSqUSJSUlWtfVt32oTVXj08ZS9+G4cePw22+/IS0tTf1yd3fHrFmz8O2332rdpr7tv9qMURtL3Yfa3LhxA3fv3tVZq6n2IU+TGciDBw80km1WVhbS0tLQvHlztG7dGr///juuX7+unrclIyMDgOroT0W6HT9+PDw8PBATEwMAmD59Ovr164clS5bghRdewLZt23Dy5EmsXr3axKNTMcYYP/zwQzz11FNo37497t+/j0WLFuHatWv4+9//buLRVT0+Nzc3jBw5Eqmpqfjmm2+gUCjU56ubN28OmUwGAAgNDcXw4cPx5ptvAgCioqIQGRmJgIAABAUFIS4uDg8fPsTEiRMbxPjefvttvPjii2jTpg1yc3MRHR0NqVSKMWPGmHx8QNVjbNGiBT766CMMHToUbm5uuHPnDlasWIGcnBy8/PLL6m3q6z6s7fgsaR9W9zvmyYBuY2MDV1dXdOrUSb3MkvcfYJwx1pd92Lx5c3zwwQcYMWIEXF1dcfnyZbzzzjto3749wsPDdY7PJPvQYPelidzhw4cFAJVekZGRgiCobiHXtj46OlrdR79+/dTtK+zYsUPo2LGjIJPJhG7dugn79u0z3aCeYIwxzpgxQ2jdurUgk8kEFxcX4fnnnxdSU1NNO7A/VDW+iukCtL0OHz6s7qNNmzYa4xUEQfj000/VYwwKChJ++eUX0w7sD8YY36hRowQ3NzdBJpMJHh4ewqhRo4RLly6ZfnB/qGqMjx49EoYPHy64u7sLMplMcHNzE4YOHSqcOHFCo4/6ug9rOz5L2ofV/Y55krbbzi15/wmCccZYX/ZhcXGxMHDgQKFly5aCjY2N0KZNG+Ef//iHkJeXp9GHOfahRBAEoY55ioiIiKje4jVDREREJGoMQ0RERCRqDENEREQkagxDREREJGoMQ0RERCRqDENEREQkagxDREREJGoMQ0RERCRqDENEREQkagxDREREJGoMQ0Rksfr3748ZM2aYuwwiauAYhohItPr16weJRIKPP/5YY7kgCAgODoZEIsGHH35opuqIyFQYhohIlARBwKlTp9CmTRucOXNGY93GjRuRm5sLAPDz8zNHeURkQgxDRFQvlJSUYNq0aWjVqhXs7Ozw9NNP49dff9VoU1RUhLFjx6JJkyZwc3PDJ598ovNUW2ZmJoqKihAZGakRhoqKijB37lxMmDABAODv72/MYRGRBWAYIqJ64Z133sGuXbuwceNGpKamon379ggPD8fvv/+ubhMVFYWffvoJe/fuxaFDh5CcnIzU1FSt/aWkpKBx48YYM2YMMjIyUFpaCgCYP38+AgIC0LJlS7i6usLNzc0k4yMi82EYIiKL9/DhQ6xcuRKLFi3C4MGD0bVrV6xZswaNGjXC2rVrAaiO6GzcuBGLFy9GaGgounfvjvXr10OhUGjtMzU1FT169ECnTp1gZ2eHCxcuIDMzEytXrsTSpUuRmppao1Nkubm5GDt2rEHGS0SmxTBERBbv8uXLKCsrQ9++fdXLbGxsEBQUhPPnzwMArly5grKyMgQFBanbODo6olOnTlr7rAg7EokEPXr0wJkzZ/DWW2/hjTfeQIcOHZCSklKjU2Tu7u744osvajlCIjInhiEiEqW/Hvnp2bMn4uLicPLkScybNw+PHz/GhQsX1OuvXr0KX19fjB07Fh06dMAbb7yBPXv2IDg4GN27d0dmZiauXr2KgIAAjfaRkZHo0qULRo0aBUEQzDZWIqoawxARWbx27dpBJpPhp59+Ui8rKyvDr7/+iq5duwIAvL29YWNjo3FRdUFBAS5evFipvytXruD+/fvqsNOrVy+cPHkSMTExsLe3x+nTp1FeXq5xZOj8+fN47733cOHCBSQlJeGnn37C8ePHMXXqVMTHx1d6j/Pnz2P27NlIT09Hfn4+jh49arA/DyIyLGtzF0BEVJ0mTZrgjTfewKxZs9C8eXO0bt0aCxcuRHFxMSZNmgQAsLe3R2RkpLpNq1atEB0dDSsrK0gkEo3+UlJSIJPJ0L17dwBAZGQkXnrpJbRo0QKA6qhRy5Yt4eXlpd6mU6dO6lNuXbp0QVhYGADAx8cH+/fvr1Rzp06d1EGtV69euHr1KkJCQgz8J0NEhsAwRET1QmxsLJRKJcaNG4eioiIEBATg22+/hZOTk7rN0qVL8frrr2PIkCFwcHDAO++8g+zsbNjZ2Wn0lZqaiu7du8PGxgaA6vojZ2dnjfW9evXS2MbW1lb9vZWVlfpnKysrrRdp/7W9VCrVeSE3EZkfwxARWaykpCT193Z2dli+fDmWL1+us729vb3GRcwPHz7EBx98gMmTJ2u0i4mJQUxMjM5+1qxZU/uiiajeYRgiogbj1KlTuHDhAoKCglBQUKB+lMawYcPMXBkRWTKJwFsciKiBOHXqFP7+978jIyMDMpkM/v7+WLp0KXx8fMxdGhFZMIYhIiIiEjXeWk9ERESixjBEREREosYwRERERKLGMERERESixjBEREREosYwRERERKLGMERERESixjBEREREosYwRERERKLGMERERESi9v80SrOHD+c+dgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from galtab.jaxhalotools import JaxZheng07Cens, JaxZheng07Sats\n",
    "\n",
    "# Create JAX-compatible composite HOD model\n",
    "hod_jax = htem.HodModelFactory(\n",
    "    centrals_occupation=JaxZheng07Cens(threshold=-21),\n",
    "    satellites_occupation=JaxZheng07Sats(threshold=-21),\n",
    "    centrals_profile=htem.TrivialPhaseSpace(),\n",
    "    satellites_profile=htem.NFWPhaseSpace()\n",
    ")\n",
    "\n",
    "# Define function that predictions P(N_cic = 1)\n",
    "def calc_cic1(logMmin=12.79):\n",
    "    hod_jax.param_dict.update({\"logMmin\": logMmin})\n",
    "    return cictab.predict(hod_jax, warn_p_over_1=False)[1]\n",
    "\n",
    "# Define the derivative of calc_cic1\n",
    "diff_cic1 = jax.grad(calc_cic1)\n",
    "\n",
    "# Note that we shouldn't make logMmin too much lower than that of our fiducial\n",
    "# model. If desired, make more conservative choices for the fiducial parameters.\n",
    "# i.e., low logMmin / logM1 / logM0 values and large sigma_logM values\n",
    "for logmmin in np.linspace(11.0, 15.0, 20):\n",
    "    value = calc_cic1(logmmin)\n",
    "    derivative = diff_cic1(logmmin)\n",
    "\n",
    "    plt.plot(logmmin, value, \"bo\")\n",
    "    plt.quiver(logmmin, value, 1, derivative, angles=\"xy\")\n",
    "\n",
    "plt.xlabel(\"$\\\\log M_{\\\\rm min}$\")\n",
    "plt.ylabel(\"$P(N_{\\\\rm CiC} = 1)$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## `jax.grad` (the arrows in the above plot) isn't working *yet*...\n",
    "\n",
    "- I actually wasn't expecting the above to work perfectly, because it's using the Monte-Carlo mode, which isn't perfectly continuous\n",
    "- But analytic mode moment derivatives aren't working either...\n",
    "- TODO: Figure out what's going wrong"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "main39",
   "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.13"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}