{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quickstart Tutorial\n",
    "\n",
    "MultiGrad allows you to easily implement data-parallelized, differentiable models in the Jax framework. All you need to provide is two Jax traceable functions: (1) a function that accepts model parameters and predicts a set of summary statistics that multigrad will automatically **sum over all MPI ranks**, and (2) a function that accepts the summary statistic predictions, compares them to a set of targets, and returns a loss function. All of the specific details of these functions are up to you!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mpi4py import MPI\n",
    "import jax.scipy\n",
    "from jax import numpy as jnp\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from diffopt import multigrad"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example Model\n",
    "\n",
    "Let's imagine we have a model that probabilistically predicts the stellar mass of a galaxy residing in a halo based on two parameters: a stellar-to-halo mass ratio $f$ and a scatter of $\\sigma$ (both of which must be positive, and could reasonably span a range of orders of magnitude, so we will actually parameterize their logarithmic counterparts $\\log f$ and $\\log \\sigma$). We will define the probability of a stellar mass given the halo mass using a log-normal distribution $P(\\log M_\\star | M_h; f, \\sigma) = \\mathcal{N}(\\log f + \\log M_h, \\sigma)$. Imagine we have a simulated dataset of 10,000 halos with masses ranging from $10^{10}$ to $10^{11}$ $M_\\odot$. We want to fit our model parameters to reproduce an observed stellar mass function, which is the number density of galaxies as a function of stellar mass $\\Phi(\\log M_\\star)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate fake halo masses between 10^10 < M_h < 10^11 as a power law\n",
    "def load_halo_masses(num_halos=10_000, comm=MPI.COMM_WORLD):\n",
    "    quantile = jnp.linspace(0, 0.9, num_halos)\n",
    "    mhalo = 1e10 / (1 - quantile)\n",
    "\n",
    "    # Assign halos evenly across given MPI ranks (only one rank for now)\n",
    "    return np.array_split(mhalo, comm.size)[comm.rank]\n",
    "\n",
    "\n",
    "# Compute one bin of the stellar mass function (SMF)\n",
    "@jax.jit\n",
    "def calc_smf_bin(params, logsm_low, logsm_high, volume, log_halo_masses):\n",
    "    # Unpack model parameters f and sigma\n",
    "    log_f, log_sigma = params\n",
    "    mean_logsm = log_f + log_halo_masses\n",
    "    sigma = 10 ** log_sigma\n",
    "\n",
    "    # Integrating the log-normal PDF over the bin, we get erf:\n",
    "    erf_high = 0.5 * (1 + jax.scipy.special.erf(\n",
    "        (logsm_high - mean_logsm)/(jnp.sqrt(2) * sigma)))\n",
    "    erf_low = 0.5 * (1 + jax.scipy.special.erf(\n",
    "        (logsm_low - mean_logsm)/(jnp.sqrt(2) * sigma)))\n",
    "    prob_in_bin = erf_high - erf_low\n",
    "\n",
    "    # Sum probabilities, convert to number density, divide by bin width\n",
    "    return jnp.sum(prob_in_bin) / volume / (logsm_high - logsm_low)\n",
    "\n",
    "\n",
    "# Compute the stellar mass function over all bins (loop over calc_smf_bin)\n",
    "@jax.jit\n",
    "def calc_smf(params, smf_bin_edges, volume, log_halo_masses):\n",
    "    smf = []\n",
    "    logsm_low = smf_bin_edges[0]\n",
    "    for logsm_high in smf_bin_edges[1:]:\n",
    "        smf_bin = calc_smf_bin(\n",
    "            params, logsm_low, logsm_high, volume, log_halo_masses)\n",
    "        smf.append(smf_bin)\n",
    "        logsm_low = logsm_high\n",
    "    return jnp.array(smf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inspecting the \"truth\"\n",
    "\n",
    "We've now defined a model that will predict the stellar mass function, given two parameters and a set of halo data. Now we just have some auxiliary data to define, such as the volume spanned by our 10,000 halos and the stellar mass bins over which we are measuring the stellar mass function. Let's set a volume of 1.0 for simplicity (in whatever units happens to make this a reasonable assumption) and ten bins evenly dividing the range over $9 < \\log M_\\star < 10$.\n",
    "\n",
    "Let's also set fiducial \"truth\" values for our stellar-to-halo mass fraction $\\log f = -2$ and $\\log \\sigma = -0.5$ (roughly corresponding to $f = 0.01$ and $\\sigma = 0.316$). Then, let's go ahead and \"load\" our simulated halos and compute the \"true\" stellar mass function, which will be our target summary statistic."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "An NVIDIA GPU may be present on this machine, but a CUDA-enabled jaxlib is not installed. Falling back to cpu.\n"
     ]
    }
   ],
   "source": [
    "volume = 1.0\n",
    "smf_bin_edges = jnp.linspace(9, 10, 11)\n",
    "\n",
    "true_params = jnp.array([-2.0, -0.5])\n",
    "log_halo_masses = jnp.log10(load_halo_masses(10_000))\n",
    "true_smf = calc_smf(true_params, smf_bin_edges, volume, log_halo_masses)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1oAAAF8CAYAAADWyMqcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABdiElEQVR4nO3deVxV1f7/8fc5IOAAKCqDM5qpiIIizmODaWVzWqlpal2VLO12u/azMhuuZWVmouZQWpqZ2XArG2xyTnEAB8wBcUZJUMABFNi/P/zCjQRlOOdszuH1fDzOo9hnnb0/a4t7+zlr7c+yGIZhCAAAAABgM1azAwAAAAAAV0OiBQAAAAA2RqIFAAAAADZGogUAAAAANkaiBQAAAAA2RqIFAAAAADZGogUAAAAANkaiBQAAAAA25m52AOVdbm6ujh8/Lm9vb1ksFrPDAYAKxTAMZWRkqE6dOrJa+W4wD/cmADBHSe5LJFrXcPz4cdWvX9/sMACgQjty5Ijq1atndhjlBvcmADBXce5LJFrX4O3tLenyyfTx8TE5GgCoWNLT01W/fv38azEu494EAOYoyX2JROsa8qZk+Pj4cDMDAJMwPa4g7k0AYK7i3JeY8A4AAAAANkaiBQAAAAA2RqIFAAAAADZGogUAAAAANkaiBQAAAAA2RqIFAAAAADZGeXcAACqQnFxDmxJTlZyRKX9vL7UP9pOblfL5AGBrJFoAAFQQ3+9M0qSv45WUlpm/LcjXSxP7hahPaJCJkQGA62HqIAAAFcD3O5M0atHWAkmWJJ1Iy9SoRVv1/c4kkyIDANdEogUAgIvLyTU06et4GYW8l7dt0tfxysktrAUAoDRItAAAcHGbElOvGMn6K0NSUlqmNiWmOi4oAOVSTq6hDQkp+ir2mDYkpDj9FzC//fabLBaLzpw54/Bj84xWEaKjoxUdHa2cnJwy76vR+G+LfO/ga7eVef8AAFxNckbRSVZp2gFwTWY8x9mzZ0+Fh4dr2rRp5WpftsCIVhGioqIUHx+vmJgYs0MBAKBM/L29itVuyabD2p+cYedoAJRH5fU5TsMwlJ2dbcqxy4pECwAAJxEdHa2QkBBFRkaW6HPtg/0U5OulaxVx//1Aqnq/vVrjlsbq4KlzpQ8UgOkMw9D5i9nFemVkXtLE/+666nOcL/43XhmZl665L8Mo/lTDoUOHatWqVXrnnXdksVhksVi0YMECWSwW/fDDD2rXrp08PT21Zs0aDR06VHfddVeBz48dO1Y9e/Yscl8HDx7Mb7tlyxa1a9dOVapUUefOnbVnz54Snc/SYOogAABOIioqSlFRUUpPT5evr2+xP+dmtWhivxCNWrRVFqnAP6bykq8Jt7XQpsRU/Rh/Ul9sO6b/xh3XvW3raswNTVXfr4otuwHAAS5cylHICz/YZF+GpBPpmWr14o/XbBv/0i2q4lG8FOOdd97R3r17FRoaqpdeekmStGvXLknSM888ozfffFONGzdW9erVS7Wv2rVr5ydbEyZM0FtvvaXatWtr5MiRGjZsmNatW1esOEuLRAsAgAqgT2iQZg1qe8XzF4F/ef5iRLfG2nE0TVNX7tGve/7Up5uP6ottx9S/XX09fsN1CvKtbGIPALgaX19feXh4qEqVKgoMDJQk/fHHH5Kkl156STfffHOZ9vVXr776qnr06CFJGj9+vG677TZlZmbKy6t4U6tLg0QLAIAKok9okG4OCdSmxFQlZ2TK39tL7YP95Gb936TCVvV89cEj7bXl0Gm9vXKv1u4/pcUbD2vZlqN6qH0Dje7VpNjPfAEwT+VKbop/6ZZitd2UmKqhH1y7LsGCRyLVPtjvmse1hXbt2tlkP3lat26d//9BQZcLeyQnJ6tBgwY2Pc5fkWiZrKiKhFQjBADYg5vVok5Nal6zXUTDGlo0ooM2HkjRWyv3alNiqhasP6hPYg7r4U6N9I/ujVWzmqcDIgZQGhaLpdhT+Lo1ra0gXy+dSMss9Dktiy6PfndrWrvAFzP2VLVq1QI/W63WK57/unTpUrH3V6lSpfz/t1gu9yE3N7cMEV4bxTAAAECROjSuqaWPddSi4R3UpkF1ZV7K1ZzVB9Rtyq9644c/dOb8RbNDBFBGec9xSrqiaE7ezxP7hdglyfLw8CjWckq1a9dWUlLByoexsbGl2pejkGgBAICrslgs6tq0lj4f1VkfDI1UaF0fnb+Yo+hfE9Tt9V817ae9Ss8s/jfLAMqfvOc4A30LTg0O9PXSrEFt7baOVqNGjbRx40YdPHhQp06dKnKU6YYbbtDmzZv14Ycfat++fZo4caJ27txZqn05CokWAAAoFovFol7N/fX141313uAINQ/0VkZWtqb9tE/dXv9V0b/u17ks51zvBsDlZGvtv2/Qkkc76p0HwrXk0Y5a++8b7JZkSdLTTz8tNzc3hYSEqHbt2jp8+HCh7W655RY9//zzeuaZZxQZGamMjAw9/PDDpdqXo1iMkhS7r4DySuimpaXJx8enVPso6jmsq+EZLQCwzTXYFZWX85Kba2jFziS9vXKvEv68vO5WzaoeGtmjiQZ1bKjKHrZ5KB4AyouSXH8phlFOXS05IwkDAJQHVqtFt7euo76hQfpv3DG989M+HUw5r1dX7NacNQcU1bOJHuzQQJ7uJFwAKh6mDgIAgDJxs1p0d5t6+umpHppyb2vVrV5Zf2Zk6cWv49Xzjd+0eOMhXcw291kJAHA0Ei0AAGAT7m5W9Y+sr1+f7qlX7gpVoI+XktIyNeGLnbrhrd/06eYjys4pmHDl5BrakJCir2KPaUNCinJyeaIBgGtg6iAAALApD3erBnVsqPsi6unjjYc187cEHT19Qc98tl2zfkvQkzc2Vb+wOloZf0KTvo5XUlpm/meDfL00sV+IXR++BwBHYEQLAADYhVclNw3rGqw1z/TSs32bq0aVSko8dU5jl8aq62u/aOSirQWSLEk6kZapUYu26vudSUXsFQCcA4kWAACwq8oebvpHjyZa8+8b9HTv6+Xt6aak9MxC2+ZNHJz0dTzTCAE4NaYOOqGiKhJSjRAAUJ5V83TX4zc0VfMgH41YuLnIdoakpLRMbUpMVacmNR0XIADYECNaAADAoYq7qHFyRuGjXgDgDEi0AACAQ/l7exWrXc2qHnaOBADsh0QLAAA4VPtgPwX5eslyjXYvfxOv3w+kOCQmALA1Ei0AAOBQblaLJvYLkaQrkq28n6t4uGnPybN6YM7vevzjrTp+5oJDYwSAsiLRAgAADtcnNEizBrVVoG/BaYSBvl6aPait1v37Bg3s0EAWi/TN9iTd+NYqzfhlnzIv5ZgUMQCUjMUwDGqnXkV6erp8fX2VlpYmHx+fUu2jqCqBjkI1QgDOyhbXYFfkSuclJ9fQpsRUJWdkyt/bS+2D/eRm/d84185jaZr09S7FHDwtSWrgV0XP3x6im1r4y2K51uRDALCtklx/Ke8OAABM42a1XLWEe2hdX336j076b9xx/WfFbh1OPa9HP9ys7tfX1gu3h+g6/2oOjBYAio+pgwAAoFyzWCy6M7yufvlnT43q2UQeblat3vun+kxbrf+s2K2MzEtmhwgAVyDRAgAATqGqp7v+3ae5fhjXXTc291d2rqE5qw+o15ur9NmWo8rN5WkIAOUHiRYAAHAqwbWqav7QSH0wNFLBtarq1NksPb0sTvfOXq/tR8+YHR4ASCLRKlJ0dLRCQkIUGRlpdigAABeTkZGhyMhIhYeHq1WrVpo7d67ZITmlXs399f3Ybhrft7mqerhp2+EzujN6nf792XadOptldngAKjiqDl6DK1QdvBoqEgIoz1yput5f5eTkKCsrS1WqVNH58+cVGhqqmJgY1axZdFGIv3LV81IWJ9Mz9dp3f+iLbcckSd5e7hp30/Ua3KmhKrnxvTIA2yjJ9ZcrDwAADubm5qYqVapIkjIzM5WTkyO+9yybAB8vvT0gXJ+N7KSWdXyUkZmtl76J123T12j9/lNmhwegAiLRAgCghFavXq1+/fqpTp06slgs+vLLL69oM3PmTAUHB8vLy0sRERFas2ZNgffPnDmjsLAw1atXT88884xq1arloOhdW7tGfvrv4131n7tbqUaVStp78qwemrdRoxZt0dHT580OD0AFQqIFAEAJnTt3TmFhYZoxY0ah7y9dulRjx47VhAkTtG3bNnXr1k19+/bV4cOH89tUr15dcXFxSkxM1Mcff6yTJ086KnyX52a16KEODfTb0700tHMjWS3SdztP6Ma3VmnaT3uVeSnH7BABVAAkWgAAlFDfvn31yiuv6J577in0/alTp2r48OEaMWKEWrRooWnTpql+/fqaNWvWFW0DAgLUunVrrV69usjjZWVlKT09vcAL1+ZbpZJevKOlVjzZTR2C/ZSVnatpP+3TjW+t0vc7k5iuCcCuSLQAALChixcvasuWLerdu3eB7b1799b69eslSSdPnsxPltLT07V69Wo1a9asyH1OnjxZvr6++a/69evbrwMuqHmgjz55rKNmPNRGQb5eOnbmgkYu2qpB8zdq78mM/HY5uYY2JKToq9hj2pCQohzW5QJQBu5mBwBzFVURkWqEAFA6p06dUk5OjgICAgpsDwgI0IkTJyRJR48e1fDhw2UYhgzD0OOPP67WrVsXuc9nn31WTz31VP7P6enpJFslZLFYdHvrOrqhub9m/Zag91Yf0Lr9Ker7zhoN6dRILev66M0f9igpLTP/M0G+XprYL0R9QoNMjByAsyLRAgDADiwWS4GfDcPI3xYREaHY2Nhi78vT01Oenp62DK/CquLhrn/2bqb7I+rrlW/j9WP8Sb2/LrHQtifSMjVq0VbNGtSWZAtAiTF1EAAAG6pVq5bc3NzyR6/yJCcnXzHKBfM0qFlFcx5upwVDI+VmtRTaJm/i4KSv45lGCKDESLQAALAhDw8PRUREaOXKlQW2r1y5Up07dzYpKhTFs5LbVZMoQ1JSWqY2JaY6LigALoGpgwAAlNDZs2e1f//+/J8TExMVGxsrPz8/NWjQQE899ZQGDx6sdu3aqVOnTpozZ44OHz6skSNHlum40dHRio6OVk4O5cltJTkj89qNStAOAPKQaKFQRRXJkCiUAQCbN29Wr1698n/OK1QxZMgQLViwQAMGDFBKSopeeuklJSUlKTQ0VCtWrFDDhg3LdNyoqChFRUUpPT1dvr6+ZdoXLvP39ipWu+pVKtk5EgCuhkQLAIAS6tmz5zXXYBo9erRGjx7toIhQWu2D/RTk66UTaZm62p/oi1/t0pv9KymiYQ2HxQbAufGMFgAAqLDcrBZN7BciSfp7SYy8n3283JWYcl73z16vySt2K/MSUzcBXBuJFgAAqND6hAZp1qC2CvQtOI0w0NdLswe11ZpnbtA9besq15DeW31At01fo22HT5sULQBnwdRBAACcBMUw7KdPaJBuDgnUpsRUJWdkyt/bS+2D/fJLv0/tH66+oUH6f1/sUMKf53TvrPX6R48mGntTU3m6u5kcPYDyyGJca5J5BZf3wHFaWpp8fHxKtY+rFZZwRhTDAOAotrgGuyLOi3lOn7uoF7/epa9ij0uSmvpX01v9w9S6XnVzAwPgECW5/jKihRIrKnEkAQMAuLoaVT30zgNt1Dc0SM99uUP7ks/q7pnrNapHE4258TpGtwDk4xktAACAEuoTGqgfx/VQv7A6ysk1NOPX/brj3XXaeSzN7NAAlBMkWgAAAKXgV9VD7z7YRrMGtlXNqh7aczJDd0av09SVe3UxO9fs8ACYjEQLAACgDPq2CtKP47rr1laBysk1NP3nfbozep12HWd0C6jISLQAAHAS0dHRCgkJUWRkpNmh4G9qVvPUzIERmvFQG9WoUkm7k9J154x1euenfbqUw+gWUBFRdfAaqDpYdhTJAFBaVNcrHOelfPszI0vPfblDP+w6KUlqWcdHb/UPU/NA/qwAZ1eS6y8jWgAAADZU29tTswdF6J0HwlW9SiXtOp6ufu+u1Yxf9imb0S2gwiDRAgAAsDGLxaI7w+vqx3HddXNIgC7lGHrzx726e+Z67TmRYXZ4AByARAsAAMBO/L29NGdwhN4eECYfL3ftOJamfu+uVfSv+xndAlwciRYAAIAdWSwW3d2mnlY+1UM3NvfXxZxcvfHDHt07a732JzO6BbgqEi0AAAAHCPDx0rwh7fTm/WHy9nJX3NE03Tp9rd5blaCcXGqTAa6GRAsAACdBeXfnZ7FYdF9EPa0c10M9m9XWxexcTf7uD903e70S/jxrdngAbIjy7tdAeXf7ovQ7gKuhjHnhOC+uwTAMLdt8VC9/E6+MrGx5ulv1dO9mGtY1WG5Wi3JyDW1KTFVyRqb8vb3UPthPblaL2WEDFVpJrr/uDooJAAAAf2GxWNQ/sr66Nq2lfy/frjX7TunVFbv1/a4T6hcWpPdWHVBSWmZ++yBfL03sF6I+oUEmRg2guJg6CAAAYKI61Svrw2HtNfmeVqrm6a4th07rxf/GF0iyJOlEWqZGLdqq73cmmRQpgJIg0QIAADCZxWLRg+0b6NsnusrDrfB/nuU96zHp63iKZwBOgEQLAACgnDh+JlMXr7K+liEpKS1TmxJTHRcUgFLhGS2YqqhCIRTJAABURMkZmdduVIJ2AMzDiFYRKKELAAAczd/by6btAJiHRKsIUVFRio+PV0xMjNmhAAAgiS8BK4L2wX4K8vXS1Yq4WyzS+axsh8UEoHRItAAAcBJ8Cej63KwWTewXIklFJluGIQ3/cLOmfP+Hsq/yPBcAc5FoAQAAlCN9QoM0a1BbBfoWnB4Y5Ouldx8M18OdGkqSZv6WoEHzNyo5nee1gPKIYhgAAADlTJ/QIN0cEqhNialKzsiUv7eX2gf7yc1qUb+wuops5Kfxy7fr9wOpunX6Wr37YBt1alLT7LAB/AWJFsqloqoRSlQkBABUDG5WS5HJU7+wOgqp46PRi7Zqz8kMDZz3u/7Zu5lG9Wgiq/VqT3gBcBSmDgIAADihJrWr6cuoLrovop5yDemNH/Zo2MIYnT530ezQAIhECwAAwGlV9nDTm/eHacq9reXpbtVve/7UbdPXaOvh02aHBlR4JFoAAABOrn9kfX0xuouCa1XV8bRM9Z+9Qe+vTZRhGGaHBlRYJFoAAAAuIKSOj/77eBfd2ipQ2bmGXvomXqMXb1V65iWzQwMqJIphwOkUVSiDIhkAgIrO26uSoh9qq4XrD+rVFbv13c4T2p2UruiBbdWyjq/Z4QEVCiNaAAAALsRisWhol2B9+o9Oqlu9sg6mnNfdM9frk02HmUoIOBCJFgAATiI6OlohISGKjIw0OxQ4gTYNauibMV3Vq1ltXczO1fjPd+ify+J0/mK22aEBFQKJFgAATiIqKkrx8fGKiYkxOxQ4iRpVPTR/SKSe6dNMVov0+dZjuit6nfYnnzU7NMDlkWgBAAC4MKvVotE9r9PHj3ZUbW9P7T15VnfMWKuvYo+ZHRrg0ki0AAAAKoCOjWvq2ye6qlPjmjp/MUdPfhKr577coazsHLNDA1wSiRYAAEAF4e/tpUUjOmjMDddJkhb9flj3zdqgwynnTY4McD2Ud4fLoOw7AADX5ma16J+9mymiYQ2NWxqrHcfSdNu7a/TW/WHq3TLQ7PAAl8GIFgAAQAXUs5m/vn2im9o0qK6MzGw99tEW/WfFbl3KyTU7NMAlkGgBAABUUHWqV9bSxzppeNdgSdKc1Qf04JzfdSIt0+TIAOdHogUAAFCBebhb9fztIZo9qK28Pd21+dBp3Tp9jdbs+9Ps0ACnRqIFAAAA9QkN0tdjuiokyEep5y7q4fc36e2Ve5WTaygn19CGhBR9FXtMGxJSlJNrmB0uUO5RDAMur6giGRKFMgAA+KtGtarq89GdNenreC3ZdFjv/LxPP+w6odRzF5WckZXfLsjXSxP7hahPaJCJ0QLlGyNaAAAAyOdVyU2T72mltweEycPNqj9OZBRIsiTpRFqmRi3aqu93JpkUJVD+kWgBAADgCneE1ZVP5cInP+VNHJz0dTzTCIEikGgBAOAkoqOjFRISosjISLNDQQWwKTFVp85eLPJ9Q1JSWqY2JaY6LijAiZBoAQDgJKKiohQfH6+YmBizQ0EFkJxRvBLvxW0HVDQkWgAAALiCv7eXTdsBFQ1VB1GhFVWRkGqEAICKrn2wn4J8vXQiLVNFPYVlkWQYPKMFFIYRLQAAAFzBzWrRxH4hki4nVIUxJA35YJOWbznqsLgAZ0GiBQAAgEL1CQ3SrEFtFehbcHpgkK+Xpj8QrttaBelSjqF/LovT1B/3MLoF/AVTBwEALuXnn3/WL7/8ovXr1+vo0aM6deqUqlSpotq1a6tVq1bq0aOHbr/9dgUGBpodKuAU+oQG6eaQQG1KTFVyRqb8vb3UPthPblaLbm9dRw1rVtHM3xI0/Zf9OpR6Xq/f21peldzMDhswHYkWAMDpnT17VtOnT9fcuXN1+PDh/G/Vvby85OfnpwsXLmjnzp3avn27Fi9eLHd3d91xxx0aN26cunTpYnL0QPnnZrWoU5OaV2y3Wi16pk9zNapZVf/vix36Kva4jp2+oDkPt5NfVQ8TIgXKDxItoBBFFcmQKJQBlDezZ8/Wiy++qOTkZIWFhemxxx5Tp06d1K5dO1WrVi2/nWEY2rdvnzZu3Kgff/xRX331lb744gvdeeedeuuttxQcHGxiLwDn1j+yvurWqKyRi7Zo86HTunvmOr0/NFJNale79ocBF2UxyjCZtiJMz0hPT5evr6/S0tLk4+NTqn1c7R/tcD4kWoDjFOcaXKlSJQ0cOFD/+te/1LJly2Lv+8KFC1qyZIkmT56swYMH64UXXrBV2HZni3sTYA/7TmbokQUxOnr6gnwrV9J7gyPUsfGVI2GAsyrJ9bfEiVZxpmekpaUpNzdXkpx+egaJFv6ORAtwnOJcgxMSEtSkSZNSHyMnJ0dHjx5Vw4YNS70PRyPRQnl26myWHv1ws7YdPqNKbha9dk9r3RtRz+ywAJsoyfW3RFUHZ8+ereuuu07PPfecqlevrldeeUW//PKL0tPTdf78eR09elQpKSm6dOmS/vjjDy1cuFADBgzQjz/+qO7du+uee+5RYmJimToHAMBflSXJkiQ3NzenSrKA8q5WNU8tebQjFQlR4ZUo0RozZoz69OmjHTt2aNu2bXr22WfVs2fPAnPgJclisej666/X4MGD9dFHH+nkyZOaO3euduzYoY8++simHQAAAED54lXJTe8+2Eaje17+ImT6L/s1dmmsMi/lmBwZ4DglKobxxx9/lOqbw8qVK2vYsGEaMmSIjh5lQTsAAABXR0VCVHQlSrSYngEU/cwdz24BAHAlKhKioirR1EEAAACgpLpcV0ufj+qsejUq61DKed0zc71+P5BidliAXdkl0Tp48KBGjx6t/v3769lnn9Unn3yi+Pj4/EqEAACYJS0tTVu3bjU7DKDCaRrgrS+juqhNg+pKu3BJg+dv1PItPFIC12WXROu+++7T+vXr1aBBA+3bt0/PPfecWrVqpWrVqikiIsIehwQAoFiys7P1wAMP6P3339d7771ndjhAhUJFQlQkJXpGq7h2796tzZs3q0WLFvnbMjIyFBsbq+3bt9vjkAAAXNM333yjN998U/v379c///lPHThwwOyQSiQ6OlrR0dHKyaFyG5xXXkXChjWraOZvCZr+y34dSj2v1+9tLa9KbmaHB9iMXUa02rZtq7S0tALbvL291a1bN0VFRdnjkAAAXNPtt9+u5557ThaLRTVq1NCPP/5odkglEhUVpfj4eMXExJgdClAmeRUJp9zbWu5Wi76KPa5B8zYq9dxFs0MDbMZi2GisdtiwYWrdurXCwsJ0/vx5RUdHa9myZapataotdm+akqz+XJSiqtShYqAaIVB6trgG/11KSoqmTJmi119/Xfv379d1111nk/06kj3OC2CWdftPaeSiLcrIzFbDmlWoSIhyrSTXX5uNaHl6emrZsmW68847dccdd2jlypVq3ry5nn76aS1btkx79+5l/i0AwHS+vr569dVXJckpkyzA1VCREK7KZs9ozZo1K///ExISFBcXl/9avny5Dh06pCpVqqhly5bauHGjrQ4LAECJuLvb5fFkAGXQNMBbX4zuokc/3KzYI2c0eP5GvXZPa90bUc/s0IBSs8vdpkmTJmrSpInuueee/G3p6ekUwwAAAEChant76pPHOuqfn8bp2x1J+ueyOB1KOadxN18vi8VidnhAiZV46uD+/ftLdSAfHx91795djz/+eKk+DwAAANeWV5FwdM8mkqTpv+zX2KWxyrxEpU04nxKPaF1//fXy8fFReHi4IiIi1LZtW7Vt21bNmzd3qW8bKKELAADgeHkVCRvWrKIJX+zUV7HHdez0Bc15uJ38qnqYHR5QbCWuOmi1WmWxWPILW+QlV1WqVFFYWFh+4hUREaGQkBC5uTn3eghUHYQ9UZEQuLriXoOdtXpgaVF1EBVFURUJc3INbUpMVXJGpvy9vdQ+2E9uVtf5wh/lV0muvyUe0fLz81NaWppuvfVW3X333Tp48KC2bNmirVu3av369Vq/fn1+8uXp6alWrVopIiJCM2fOLF1vAAC4hooy2wKoaPIqEj6yICa/IuHwrsFasumwktIy89sF+XppYr8Q9QkNMjFaoKASj2idOXNGzz//vGbPni0fHx+99NJLGjVqlKxWq06ePJmfdOX998iRI7JYLE47BY8RLdgTI1rA1RX3GsxsC8C1/ZmRlV+RsDB5X6fMGtSWZAt2VZLrb6kXLN65c6fGjh2rX375RS1bttS0adN04403XtEuJSVFW7ZsUe/evUtzGNORaMGeSLSAqyvuNbhWrVpFzrZISkqSJJeabUGihYroXFa2Il5ZqcxLuYW+b5EU6Oultf++gWmEsBuHLFgcGhqqn376SZ999pnOnTun3r1765577lFiYmKBdjVr1nTaJAsA4Bz279+vkSNHasWKFfrnP/+p2rVr66uvvtKxY8eUlJSkb775RpMmTdIdd9yh2rVrKyYmRu+9957ZYQMoge1H04pMsiTJkJSUlqlNiamOCwq4ilInWnnuuece/fHHH3rppZe0cuVKhYSEaMKECTp37pwt4gMA4JqqV6+ud999V9u2bVObNm00ZswYhYWF6eeff1ZAQIBuvfVWPffcc/riiy906NAh/fnnn/ruu+/MDhtACSRnZF67UQnaAfZW5kRLkjw8PDRhwgTt2bNH9957ryZPnqxmzZpp165dttg94LIajf+20BeA0mG2BeC6/L29bNoOsDebJFpHjx7Vd999p0WLFslisahGjRpKSkpSQkKCLXYPAECJMNsCcD3tg/0U5Oulqz19FehzudQ7UB6UuLz7unXrtGPHjgKv9PT0/EpPtWvXVkREhMLDwxUWFmbzgAEAKI682RaPPPKInnnmGU2ePFkLFy7UDz/8oJYtW5odHoAScrNaNLFfiEYt2iqLLj+T9Xdelaw6m5kt3yqVHB0ecIUSJ1rdunWTxWKR1WpV48aN1bt3b4WHh+cnVnXq1LFHnAAAlMjRo0fzvxD8+2wLEi3AOfUJDdKsQW016ev4Auto1armofMXc3Qw5bwenPu7PhreXjWreZoYKVCKREuS3N3d1bdvX3Xv3j1/bRLKywIAzMJsC6Di6BMapJtDArUpMVXJGZny9748XXBfcoYGzduo+KR0PTDndy0e0UH+PjyvBfOUeB2t1q1b648//lB2dnb+miSS1Lhx4/zFIPOSLz8/558jyzpaMAPrawGXlXTB4rzZFnkzLVx1tgXraAGFS/jzrAbO3agT6ZlqVLOKFj/aUXWrVzY7LLgQuy9YnJWVpbi4OG3dujV/QcidO3fq0qVLl3f6fwlYgwYN8pOv//f//l8pumI+Ei2UNyRhqEhKkmhVqlSpwsy2INECinYk9bwemve7jqReUN3qlbV4RAc1qlXV7LDgIuyeaBXm0qVL2r59e4Hka8eOHcrKypLFYlFOTo4tDuNwJFoob0i0UJEU9xrMbAsAf5WUdkED527UgVPn5O/tqcUjOqhpgLfZYcEFmJJoFSY7O1u7du3Sli1bNGzYMHsdxq5ItFDekGihIinJNZjZFgD+6s+MLA2ev1F/nMiQX1UPfTS8vVrW8TU7LDi5cpNouQISLZQ3JFqoSMp6DWa2BVCxnT53UUM+2KTtR9Pk4+WuhcPaq02DGmaHBSdWkutvqaoOAgDgDCpVqqSIiAhFRETo0UcflVRwtgUA11ajqocWjeigYR/EaPOh0xo0b6PmD41Ux8Y1zQ4NFQAjWtfAiBacBSNdcEWM3BSO8wKUzPmL2RqxcLPWJ6TIq5JV7w1upx7X1zY7LDihklx/rQ6KCQAA/J8jR46oZ8+eCgkJUevWrbVs2TKzQwJcWhUPd70/NFK9mtVW5qVcPbpws37cdcLssODiyjR1sDgFLqxWq3x8fNSsWTPdfvvtqlu3blkOCQCA03N3d9e0adMUHh6u5ORktW3bVrfeequqVqUENWAvXpXc9N7gdnryk236bucJjVq8VdMGhKtfmGuts4fyo0yJ1oIFC/KrOBU2A9FisRTYPmbMGL3wwgt67rnnynJYAACcWlBQkIKCgiRJ/v7+8vPzU2pqKokWYGce7la9+2Ab/euz7fpi2zE9+ck2ZV7K0f3t6psdGlxQmRKthIQEjR07VjExMXryySfVuXNnBQQE6OTJk1q3bp2mT5+u9u3ba8KECYqLi9Mrr7yiiRMnqmnTphowYICt+gAAQD5HzLZYvXq13njjDW3ZskVJSUn64osvdNdddxVoM3PmTL3xxhtKSkpSy5YtNW3aNHXr1u2KfW3evFm5ubmqX59/6AGO4O5m1Vv3h8mrklVLNh3Rvz7brsxLORrcqZHZocHFlCnRWrp0qTZt2qS4uDj5+/vnb7/++uvVrVs3DR06VOHh4fr111/1zDPPqG/fvgoJCdHMmTNJtAAAduGI2Rbnzp1TWFiYHnnkEd17771XvL906VKNHTtWM2fOVJcuXfTee++pb9++io+PV4MGDfLbpaSk6OGHH9a8efOuerysrCxlZWXl/5yenl7sWAFcyWq16D93t5JXJTd9sO6gnv9qly5cytFj3ZuYHRpcSJmqDjZt2lR9+/bV9OnTi2wzZswYff/999q3b58kaeDAgfr222915syZ0h7Woag6CGdHNUI4s9JcgxMTE0s82+LIkSP6+OOPS/UloMViuWJEq0OHDmrbtq1mzZqVv61Fixa66667NHnyZEmXk6ebb75Zjz76qAYPHnzVY7z44ouaNGnSFdupOgiUjWEYevPHPYr+NUGSNO6m6/XEjdflf1kD/J3Dqg4ePXpUnp6eV23j5eWlo0eP5v/coEEDZWZmluWwAAAUKW+2RWxsrP7973+rW7du+TMtxo8fr61bt+r333/Xr7/+qhEjRmjdunWqVq2aZs6caZPjX7x4UVu2bFHv3r0LbO/du7fWr18v6fI/7oYOHaobbrjhmkmWJD377LNKS0vLfx05csQmsQIVncVi0b9uaa6ne18vSXr7p7167fs/Ch0NB0qqTIlW3bp19dVXXxWYzvBXWVlZ+uqrrwrMfU9OTlaNGqzIDQCwj/nz5+v+++8vMKX9rwIDA3X//fdr7ty5ki7fy26//XbFxcXZ5PinTp1STk6OAgICCmwPCAjQiROXy0mvW7dOS5cu1Zdffqnw8HCFh4drx44dRe7T09NTPj4+BV4AbOfxG5rq+dtDJEnvrTqgF/+7S7m5JFsomzIlWsOHD9f+/fvVo0cPffvtt0pNTZUkpaam6ptvvlH37t2VkJBQ4MHkNWvWKCwsrGxRAwBQhPIy2+LvU48Mw8jf1rVrV+Xm5io2Njb/1apVK5seH0DJDO8arFfvDpXFIi3ccEjjP9+uHJItlEGZimE888wz2r17txYtWqQ77rhD0uVKTrm5uZIu31QGDhyo8ePHS5JOnjyp2267TX369Clj2AAAFC5vtsUrr7xSaMJl79kWtWrVkpubW/7o1V+P8fdRLgDly8AODVW5kpueXhanTzcfVealXL3VP0yV3Mo0NoEKqky/NW5ubvrwww+1cuVKPfzwwwoPD1ejRo0UHh6uIUOGaOXKlfroo49ktV4+TEBAgN5++23dcsstNgkeAIC/M3u2hYeHhyIiIrRy5coC21euXKnOnTuXad/R0dEKCQlRZGRkmfYDoGj3tK2nGQ+1lbvVov/GHVfU4q3Kys4xOyw4oTJVHawIqDoIV0ZFQpR3pbkG5+Tk6JFHHtGiRYvyp+oVNtti4cKFslqtOnnypF577TX16dOn2F8Enj17Vvv375cktWnTRlOnTlWvXr3k5+enBg0aaOnSpRo8eLBmz56tTp06ac6cOZo7d6527dqlhg0bluJMFGSLexOAq/vlj5MauWirLmbnqsf1tTV7UIQqe7iZHRZMVpLrL4nWNZBowZWRaKG8K8s1+Oeff9aiRYu0fft2paeny8fHR2FhYRo4cKBuvPHGMsX122+/qVevXldsHzJkiBYsWCDp8oLFU6ZMUVJSkkJDQ/X222+re/fuZTpuHhItwDHW7T+lEQs368KlHHVs7Kd5QyJVzbNMT97AyTk80Vq/fr0WLFig2NjY/IPmTR/s2rVrWXdvKhItuDISLZR3JBSF47wAjhNzMFWPfBCjs1nZatOguhY80l6+lSuZHRZM4rB1tCTp6aefVrdu3TRv3jxt3rxZCQkJ2rJli+bPn68ePXroqaeeKushAAAAAFNENvLT4hEd5Fu5krYdPqOH5v6u1HMXzQ4LTqBMidaHH36oqVOnqlmzZlqyZImSkpKUnZ2tEydO6JNPPlHz5s31zjvv6MMPP7RVvAAAFMv69ev12GOPqX379mrWrJkiIyP16KOPau3atWaHBsDJhNWvrk8e66ha1Ty063i6HpizQcnptl0SAq6nTFMHO3XqpOPHj2vnzp3y9va+4v309HS1atVKQUFB+v3338sUqFmYOoiKiCmFKC9Kew1++umn9fbbbyvvFvfXYhgWi0VPPvmkpk6dapeY7Sk6OlrR0dHKycnR3r17mToIONj+5LMaOO93nUzPUnCtqlo8ooPqVK+snFxDmxJTlZyRKX9vL7UP9pOb1XLtHcLpOGzq4M6dO3XvvfcWmmRJko+Pj+655x7t2rWrLIcBAKDYXHm2RVRUlOLj4xUTE2N2KECFdJ1/NS37R2fVq1FZiafO6f7ZG/TRhoPq+vovenDu73ryk1g9OPd3dX39F32/M8nscGGyMj+jda0BsbzSugAAOMKsWbNUv359bdy4UQMGDMhfJNjf31/9+/fXhg0bVK9ePc2cOdPkSAE4owY1q+jTf3RS41pVdezMBT3/1S4lpRWcRngiLVOjFm0l2argypRohYaGavny5Tp79myh72dkZGj58uVq2bJlWQ4DAECxMdsCgL3VqV5Zix/tIPcipgfmDUNM+jpeObmspFRRlSnRGjlypI4ePapOnTpp+fLlOnXqlCTp1KlT+uyzz9S5c2cdPXpUo0aNskmwAAAUB7MtANjbwVPnlX2VJMqQlJSWqU2JqY4LCuVKmVZcGzJkiGJjY/XOO++of//+kgo+cGwYhsaMGaMhQ4aUPVIADnO1Ai4UykB5lzfb4uWXX1a1atWueJ/ZFgBsITmjeFUHi9sOrqfMz2i9/fbbWr16tYYOHarw8HA1atRI4eHheuSRR7Rq1Sq98847togTAIBiceXZFtHR0QoJCVFkZKTZoQAVnr+3l03bwfWUaUQrT9euXdW1a1db7AoAgDJx5dkWUVFRioqKyi8vDMA87YP9FOTrpRNpmSpqAmGQ7+VS76iYbJJoAQBQnrz99tu699579cEHHyg2Nlbp6eny8fFRmzZtNGTIEHXr1s3sEAE4OTerRRP7hWjUoq2ySIUmW0/dfD3raVVgJFoAAJfEbAsA9tYnNEizBrXVpK/jC5R4t1qkXEP6cMMh3RIaKB+vSiZGCbOUKNEaNmxYqQ5isVg0f/78Un0WAAAAKK/6hAbp5pBAbUpMVXJGpvy9veRX1UMPzf1dO46lacSCzVo4rL0qe7iZHSocrESJ1oIFC0p1EBItAAAAuCo3q0WdmtQssG3hsPZ6cM7v2nQwVaMWb9Gcwe3k4V7mOnRwIiVKtBITE+0VBwAnUVTpd8q+wyzMtgBQHoXW9dX7j0Rq8PyN+m3Pnxr3aaymP9CGZ7YqkBIlWg0bNrRXHAAAlEpFmm0RHR2t6Oho5eTkmB0KgGKIbOSn9wa304iFMfp2e5K8Pd01+Z5WLJpeQVAMAwDg1CrSbAvKuwPOp8f1tTX9gTaK+nirPok5omqe7ppwWwuSrQqgRInW7bffrkmTJikiIqLEB7pw4YKio6NVtWpVp1wkEgBQPjHbAkB517dVkF67t7We+Wy75q1NlE/lSnrixqZmhwU7K9ETeUeOHFH79u114403asGCBUpPT7/mZzZv3qyxY8eqYcOGeuGFF1SrVq1SBwsAAAA4o/7t6uuF20MkSVNX7tUH6yrOaHxFVaIRrdjYWH3wwQd66aWXNGzYMI0YMULNmzdX27ZtFRAQoBo1aujChQtKTU3Vvn37tHnzZqWlpclqtap///569dVX1ahRIzt1BYCZKJIBszDbAoCzGNY1WBmZ2Xr7p72a9HW8qnm66/529c0OC3ZSokTLYrFo2LBhGjp0qL799lstWLBAq1at0qJFi65oa7Va1bp1a911110aMWKE6tSpY7OgAQDIkzfbomfPnho8eLDuuece+fj4XPUzmzdv1qJFi/Txxx/r7NmzWrhwoYOiBVDRPXHjdUrPvKT5axP17+Xb5e3lrj6hQWaHBTuwGIZhlHUnu3fv1tGjR5WSkqLKlSurdu3aatmypUs8qJv3wHFaWto1b9xFKeqbfqAiYEQLZVGca7BhGPmzLQ4fPiyr1erysy1scW8CYB7DMDR++Q4t3XxEldwsmj8kUt2vr212WCiGklx/bZJouTISLaBsSLRQFiW5Bufm5haYbZGamnpFG1eZbUGiBTi/nFxDTyzZpm93JMmrklWLhndQu0Z+ZoeFayjJ9dcm5d3Pnj2rffv2yTAMNW3aVN7e3rbYLQAAxWa1WtWvXz/169dPkmvOtmAdLcB1uFktentAuM5mZWvV3j/1yIIYLXm0o0LrOu81CgWVaUQrOTlZ48aN0/Lly3Xp0iVJkru7u+69915NnTpVgYGBNgvULIxoAWXDiBbKgpGbwnFeANdx4WKOhry/SZsOpqpmVQ99OrKTmtSuZnZYKEJJrr8lKu/+V8nJyerSpYuWLFkii8WiiIgItW3bVm5ubvrkk0/UrVs3JScnl3b3AACUydmzZ7Vt2zZt3bpVGRkZZocDAIWq7OGmeUPbKbSuj1LOXdSgeRt19PR5s8OCDZQ60XriiSeUkJCg4cOH6+TJk9q0aZNiYmJ04sQJDRs2TAkJCRozZowtYwXghBqN/7bIF2APycnJGjhwoGrVqqV27dopMjJStWrV0kMPPaQTJ06YHR4AXMHHq5IWPtJeTWpXVVJapgbN26g/M7LMDgtlVKpE69ixY1q2bJl69OihuXPnFhg28/Hx0bx589SjRw8tX75chw8ftlmwAABcDbMtADirmtU8tXhER9WrUVkHU85r8PyNSjt/yeywUAalSrR++eUXSdK//vWvItuMHz9eubm5+vnnn0sXGQAAJcRsCwDOLNDXS4uGd1Btb0/9cSJDQxds0rmsbLPDQimVuOrgsGHDtH37dknShx9+qM8++6zQdhcuXJAkvfvuu1q7dq3mz59fhjABALi6v8+2+Ku82RYJCQn5sy0aNGhgUqQAULRGtapq0fAO6v/eBm07fEaPfbRZ84dEyquSm9mhoYRKnGgtWLBAkmSxWPTpp59eta3FYlFsbKzi4uJItAAAdlXc2RarVq3Szz//rEceecRRoQFAiTQL9NbCYe01cO7vWrc/RWOWbNOsgW3l7lbq8gowQYkTrQ0bNmjx4sWKjo7WkiVL1KhRo0LbHTp0SA888IDGjBmjhx56qKxxAnBBRRXEoCQ8SorZFgBcTXj96po7pJ2GfhCjlfEn9cxn2/Xm/WGyWi1mh4ZiKnGi1aFDB6WmpmrGjBk6efKkBgwYUGi7mJgYWSwW9e3bVx06dChzoAAAFIXZFgBcUecmtTTzobYauWiLPt92TNW83DXpjpayWEi2nEGpxh979eolPz8/vfLKK0pKSrri/ePHj2vSpEny8/NTr169yhwkAABXs2HDBj3++OOSpCVLlmjDhg2Fvj755BMZhqExY8Zo/fr1JkddctHR0QoJCVFkZKTZoQBwkJtCAvRW/zBZLNKHGw7pzR/3mB0SiqlUiZaXl5eef/55nTp1Sh07dtS8efO0e/du7d69W3PnzlXHjh2VmpqqF154QZ6enraOGQCAAjp06KC+ffvKMAydPHlSHTp0KPSVnJzs1LMtoqKiFB8fr5iYGLNDAeBAd4bX1St3hUqSon9N0OxVCSZHhOKwGIZhlPbDjz/+uGbOnHnF8KVhGBo9erRmzJhR5gDNlp6eLl9fX6WlpRVYL6wkWJgVKBme0UKeklyDMzMzVa9ePVmtVsXFxSkoKKjA+8ePH1dYWJgk6ejRo079RaAt7k0AnM/sVQl67bs/JEmv3h2qgR0amhxRxVOS62+Jn9H6qxkzZqhfv356//33FR8fL8Mw1LJlSw0bNky33HJLWXYNoAK72pcTJGEoSt5si3Hjxqljx456/vnn1aVLF0nS2rVr9fLLLys1NVXTpk1z6iQLQMU1skcTZWReUvSvCXruy52q5umuO8Prmh0WilCmREuSbrnlFpIqAEC58OSTT2rfvn2aOXOm/vGPfxR4L2+2BQsWA3BmT/dupozMbH244ZD++Wmcqnm668YWAWaHhUJQjB8A4FJmzJih7777Tvfdd59CQkLUokUL3Xffffruu+9cYko7gIrNYrHoxX4tdXebusrONTRq8VatTzhldlgoRJlHtAAAKG+YbQHAlVmtFr1xX2udzcrWyviTenThZi1+tKPC61c3OzT8BSNaAAAAgJNxd7Pq3QfbqMt1NXXuYo6GvL9Je05kKCfX0IaEFH0Ve0wbElKUk1vquncoI0a0AAAAACfkVclNcwa306D5G7Xt8BndP3u9PN3d9OfZrPw2Qb5emtgvRH1Cg66yJ9gDI1oAAACAk6rq6a4FQ9urbnUvpWdmF0iyJOlEWqZGLdqq73cmmRRhxcWIFgCnUlTpd8q+AwAqqmpe7rqUU/gUQUOSRdKkr+N1c0ig3KyWQtvB9hjRAgAAAJzYpsRUJWdkFfm+ISkpLVObElMdFxRItAAAAABnlpyRadN2sA0SLQAAnER0dLRCQkIUGRlpdigAyhF/by+btoNtkGgBAOAkoqKiFB8fr5iYGLNDAVCOtA/2U5Cvl6729FWQr5faB/s5LCZUkGIYR44c0eDBg5WcnCx3d3c9//zzuv/++80OC4ANUSQDAFBRuVktmtgvRKMWbZVFl5/J+rt72talEIaDVYgRLXd3d02bNk3x8fH66aefNG7cOJ07d87ssAAAAACb6BMapFmD2irQt+D0wMqV3CRJC9Yd1O6kdDNCq7AqxIhWUFCQgoIuL9Lm7+8vPz8/paamqmrVqiZHBgAAANhGn9Ag3RwS+H9VCDPl7+2l8PrV9ciCTfr9QKqGL4jRl4934VktB3GKEa3Vq1erX79+qlOnjiwWi7788ssr2sycOVPBwcHy8vJSRESE1qxZU+i+Nm/erNzcXNWvX9/OUQMAAACO5Wa1qFOTmrozvK46Nampyh5umj0oQo1rVdXxtEw9unCzLlzMMTvMCsEpEq1z584pLCxMM2bMKPT9pUuXauzYsZowYYK2bdumbt26qW/fvjp8+HCBdikpKXr44Yc1Z84cR4QNAAAAmK56FQ+9PzRSNapUUtzRND31aaxycwtf4Bi24xSJVt++ffXKK6/onnvuKfT9qVOnavjw4RoxYoRatGihadOmqX79+po1a1Z+m6ysLN1999169tln1blz5yKPlZWVpfT09AIvAAAAwJk1qlVV7w1uJw83q77beUJTfthjdkguz+mf0bp48aK2bNmi8ePHF9jeu3dvrV+/XpJkGIaGDh2qG264QYMHD77q/iZPnqxJkybZLV4AjlVUNUKJioQAgIqlfbCfXr+vlcYtjdPsVQkKrlVFAyIbmB2Wy3KKEa2rOXXqlHJychQQEFBge0BAgE6cOCFJWrdunZYuXaovv/xS4eHhCg8P144dOwrd37PPPqu0tLT815EjR+zeBwAAAMAR7m5TT0/c2FSSNOGLnVq//5TJEbkupx/RymOxFFwXwDCM/G1du3ZVbm5usfbj6ekpT09Pm8cHAAAAlAfjbmqqg6fO6b9xxzVy0RZ9PrqLrvOvZnZYLsfpR7Rq1aolNze3/NGrPMnJyVeMcgEAAAAVncVi0ZT7WiuiYQ2lZ2Zr2IIYpZ67aHZYLsfpEy0PDw9FRERo5cqVBbavXLnyqkUvAAAAgIrKq5Kb5gyOUH2/yjqcel6PfbhZmZco+25LTpFonT17VrGxsYqNjZUkJSYmKjY2Nr98+1NPPaV58+bp/fff1+7duzVu3DgdPnxYI0eONDFqAAAAoPyqWc1THwyNlLeXuzYfOq1/L98uw6Dsu604xTNamzdvVq9evfJ/fuqppyRJQ4YM0YIFCzRgwAClpKTopZdeUlJSkkJDQ7VixQo1bNjQrJABOIGiKhJSjRAAUFFc5++t2YMiNOT9Tfoq9riCa1XV2JuuNzssl+AUiVbPnj2vmV2PHj1ao0ePdlBEAAAAgGvocl0tvXJXqMZ/vkPTftqnRjWr6q42dc0Oy+k5xdRBAAAAAPbzQPsG+kePxpKkZz7brs0HU02OyPmRaAEAAADQv29prj4tA3UxJ1ePfbRFh1LOmR2SUyPRAgAAACCr1aK3B4SrdT1fpZ67qEcWxCjt/CWzw3JaTvGMFgA4UlFFMiQKZcBc0dHRio6OVk4OJZgB2EdlDzfNe7id7oxepwN/ntPIRVu0cFh7ebgzPlNSnDEAAJxEVFSU4uPjFRMTY3YoAFyYv4+X5g+JVFUPN204kKLnv9xJ2fdSINEqQnR0tEJCQhQZGWl2KAAAAIBDhdTx0YyH2spqkZZuPqLZqw6YHZLTIdEqAt8aAgAAoCLr1dxfE/u1lCS9/v0fWrEjyeSInAuJFgAAAIBCDencSEM7N5IkjVsaq9gjZ0yNx5mQaAEAAAAo0nO3tVCvZrWVlZ2rEQs369iZC2aH5BSoOggAJVBURUKqEQIAXJW7m1XvPtRW981arz9OZGj4ghgtG9lJ3l6VzA6tXGNECwAAAMBVVfN01/yhkart7ak/TmTo8Y+3KTsn1+ywyjUSLQAAAADXVLd6Zc0f0k5elaxatfdPTfo6nrLvV0GiBQAAAKBYWterrmkD2shikT76/ZAWrD9odkjlFokWAAAAgGLrExqo8X2aS5Je/iZeP+8+aXJE5ROJFgAAAIASeax7Yz0QWV+5hjRmyTbFH083O6Ryh6qDAGADRVUjlKhICABwPRaLRS/fFaojp89r3f4UDV8Yoy+juijAx8vs0MoNRrQAAAAAlFglN6tmDoxQk9pVlZSWqeELY3T+YrbZYZUbJFoAAAAASsW3ciV9MLS9/Kp6aOexdD35SaxycqlEKJFoAQAAACiDBjWraO7DEfJwt2pl/Em9/v0fZodULpBoFSE6OlohISGKjIw0OxQAAACgXIto6Kc37mstSZqz+oA+3njY5IjMR6JVhKioKMXHxysmJsbsUAAAAIBy787wuhp30/WSpOe/2qk1+/5UTq6hDQkp+ir2mDYkpFSoaYVUHQQAOyuqIiHVCAEAruaJG6/TwZRz+mLbMT26cLOqebnr1NmL+e8H+XppYr8Q9QkNMjFKx2BECwAAAIBNWCwWvXZvKzWpXVWZ2bkFkixJOpGWqVGLtur7nUkmReg4JFoAAAAAbMbdalVGZuFl3vMmDk76Ot7lpxGSaAEAAACwmU2JqUrOyCryfUNSUlqmNiWmOi4oE5BoAQAAALCZ5IxMm7ZzVhTDAACTUCQDAOCK/L29bNrOWTGiBQAAAMBm2gf7KcjXS5artAny9VL7YD+HxWQGEi0AAExw9913q0aNGrrvvvvMDgUAbMrNatHEfiGSVGSy9dxtLeRmvVoq5vxItAAAMMETTzyhDz/80OwwAMAu+oQGadagtgr0LTg9MC+12n40zfFBORjPaAEAYIJevXrpt99+MzsMALCbPqFBujkk8P+qEGbK39tLKWez9PiSbXpv9QG1aVDdpRcuZkQLAIASWr16tfr166c6derIYrHoyy+/vKLNzJkzFRwcLC8vL0VERGjNmjWODxQATOZmtahTk5q6M7yuOjWpqdvD6mhE12BJ0tPLtivhz7MmR2g/JFoAAJTQuXPnFBYWphkzZhT6/tKlSzV27FhNmDBB27ZtU7du3dS3b18dPnzYwZECQPnz777N1b6Rn85mZWvUoi06l1X44sbOjqmDAFDOFFX2XaL0e3nRt29f9e3bt8j3p06dquHDh2vEiBGSpGnTpumHH37QrFmzNHny5BIfLysrS1lZ/1v8Mz09veRBA0A5UcnNqhkD2+j26Wu19+RZjf98h6Y/EC6LxbWKYzCiBQCADV28eFFbtmxR7969C2zv3bu31q9fX6p9Tp48Wb6+vvmv+vXr2yJUADCNv7eXoge2lbvVoq/jjmvB+oNmh2RzJFoAANjQqVOnlJOTo4CAgALbAwICdOLEifyfb7nlFt1///1asWKF6tWrp5iYmCL3+eyzzyotLS3/deTIEbvFDwCOEtnIT8/e2kKS9Oq3u7X5YKrJEdkWUweLEB0drejoaOXk5JgdCgDACf19CoxhGAW2/fDDD8Xel6enpzw9PW0WGwCUF8O6NNK2w6f1zfYkRX28Vd+M6aba3q5xvWNEqwhRUVGKj4+/6jeMAAD8Xa1ateTm5lZg9EqSkpOTrxjlAoCKzmKx6PV7W+s6/2o6mZ6lxz/equycXLPDsglGtADAiRRVKIMiGeWHh4eHIiIitHLlSt19993521euXKk777zTxMgAoHyq6umu2YMidOeMtdqYmKopP+zR//u/KYXOjBEtAABK6OzZs4qNjVVsbKwkKTExUbGxsfnl25966inNmzdP77//vnbv3q1x48bp8OHDGjlyZJmOGx0drZCQEEVGRpa1CwBQrlznX01v3B8mSZqz+oC+25FkckRlx4gWAAAltHnzZvXq1Sv/56eeekqSNGTIEC1YsEADBgxQSkqKXnrpJSUlJSk0NFQrVqxQw4YNy3TcqKgoRUVFKT09Xb6+vmXaFwCUN7e2CtKj3YI1d02i/vXZdl0f6K0mtauZHVapkWgBAFBCPXv2lGEYV20zevRojR492kERAYBr+Hef5oo7mqZNiaka+dEWfRnVRVU9nTNlYeogAAAAgHLB3c2qGQ+1kb+3p/Yln9W/l2+/5hdb5ZVzpocAgAKKKpIhUSgDAOBc/L29NHNgWz0w53d9sz1JbRvU0LCuwWaHVWKMaAEA4CQohgGgomjXyC+/8uB/VjjnYsYkWgAAOAnWeARQkTzSpZFubx2k7FxDoxdvVXJGptkhlQiJFgAAAIByJ28x46b+1ZSckaUxH29zqsWMSbQAAAAAlEtVPd01e3CEqnm65y9m7CxItAAAAACUW01qV9Mb97WWdHkx4xVOspgxVQcBwMUVVZGQaoQAAGfRt1WQHuveWHNWH9C/lsXp+gBvXedfvhczZkQLAAAnQdVBABXZM7c0U4dgP527mKORi7boXFa22SFdFYkWAABOgqqDACqyy4sZt1WAj6f2O8FixiRaAAAAAJxCbW9PzRzYVu5Wi77ZnqQP1h00O6QikWgBAAAAcBoRDf004bb/LWYcU04XMybRAgAAAOBUhnZupDvC6ig711BUOV3MmKqDRYiOjlZ0dLRycnLMDgUA7IJqhAAAZ2WxWDT5nlbanZSufcln9fjH27R4RAdVcis/40jlJ5JyhgeOAQAAgPLrr4sZb0pM1ZTv/zA7pAJItAAAcBKUdweAgprUrqY377+8mPHcNYnlajFjEi0AAJwEsy0A4Ep9QoP0j+6NJUn/Whan/ckZJkd0GYkWAAAAAKf2r1uaqWPjvMWMt5aLxYxJtAAAAAA4NXc3q9598H+LGT9TDhYzJtECAAAA4PT+upjxt9uT9L7JixlT3h0AUEBRZd8lSr8DAMq3iIZ+eu62Fnrx63hNXrFbrer6qn2wnymxMKIFAAAAwGUM+etixh9vVXK6OYsZk2gBAAAAcBkWi0Wv3dtK1wdU058ZWXr84226lJPr8DhItAAAAAC4lCoe7po96P8WMz6Yqte/c/xixjyjBQAotqKe3+LZLceIjo5WdHS0cnJyzA4FAMq9xrWr6c37wzRy0RbNW5uosHrVVcvbU8kZmfL39lL7YD+5WS12Oz6JFgAATiIqKkpRUVFKT0+Xr6+v2eEAQLnXJzRQ/+jRWO+tOqAnPtmmvxZ8D/L10sR+IeoTGmSXYzN1EAAAAIDLal338hdTf19V60RapkYt2qrvdybZ5bgkWgAAAABcUk6uoVe+3V3oe3mJ16Sv45WTa/vFjUm0AAAAALikTYmpSkorury7ISkpLVObElNtfmye0QIAlBmLHAMAyqPkjOKtoVXcdiXBiBYAAAAAl+Tv7WXTdiVBogUAAADAJbUP9lOQr5eKKuJu0eXqg+2D/Wx+bBItAAAAAC7JzWrRxH4hknRFspX388R+IXZZT4tECwAAAIDL6hMapFmD2irQt+D0wEBfL80a1NZu62hRDAMAACcRHR2t6Oho5eTkmB0KADiVPqFBujkkUJsSU5WckSl/78vTBe0xkpWHRAsAACcRFRWlqKgopaeny9fX1+xwAMCpuFkt6tSkpsOOx9TBIkRHRyskJESRkZFmhwIAAADAyZBoFSEqKkrx8fGKiYkxOxQAAAAAToZECwAAAABsjEQLAAAAAGyMRAsAAAAAbIxECwAAAABsjPLu12AYhiQpPT291PvIzTpvq3AAwOmU5fqZ99m8azEus8W9CQBQciW5L5FoXUNGRoYkqX79+iZHAgDOyXda2feRkZHBulF/wb0JAMxVnPuSxeBrwqvKzc3V8ePH5e3tLYul5CtHp6enq379+jpy5Ih8fHzsEGH5Rv/pf0Xuv8Q5KGv/DcNQRkaG6tSpI6uV2e55ynpvckUV/e/atXB+isa5KRrn5koluS8xonUNVqtV9erVK/N+fHx8KvQvKP2n/xW5/xLnoCz9ZyTrSra6N7miiv537Vo4P0Xj3BSNc1NQce9LfD0IAAAAADZGogUAAAAANkaiZWeenp6aOHGiPD09zQ7FFPSf/lfk/kucg4refzgOv2tXx/kpGuemaJybsqEYBgAAAADYGCNaAAAAAGBjJFoAAAAAYGMkWgAAAABgYyRaAAAAAGBjJFolsHr1avXr10916tSRxWLRl19+WeB9wzD04osvqk6dOqpcubJ69uypXbt2XXO/y5cvV0hIiDw9PRUSEqIvvvjCTj0oG3v0f+7cuerWrZtq1KihGjVq6KabbtKmTZvs2IvSs9eff55PPvlEFotFd911l20DtyF7nYMzZ84oKipKQUFB8vLyUosWLbRixQo79aL07NX/adOmqVmzZqpcubLq16+vcePGKTMz0069KL1r9f/zzz/XLbfcolq1aslisSg2NrZY+3WWayDMlZGRobFjx6phw4aqXLmyOnfurJiYmCLbf/7557r55ptVu3Zt+fj4qFOnTvrhhx8cGLHjlPTc/NW6devk7u6u8PBw+wZpotKcn6ysLE2YMEENGzaUp6enmjRpovfff99BETtOac7N4sWLFRYWpipVqigoKEiPPPKIUlJSHBSxcyHRKoFz584pLCxMM2bMKPT9KVOmaOrUqZoxY4ZiYmIUGBiom2++WRkZGUXuc8OGDRowYIAGDx6suLg4DR48WP3799fGjRvt1Y1Ss0f/f/vtNz344IP69ddftWHDBjVo0EC9e/fWsWPH7NWNUrNH//McOnRITz/9tLp162brsG3KHufg4sWLuvnmm3Xw4EF99tln2rNnj+bOnau6devaqxulZo/+L168WOPHj9fEiRO1e/duzZ8/X0uXLtWzzz5rr26U2rX6f+7cOXXp0kWvvfZasffpTNdAmGvEiBFauXKlPvroI+3YsUO9e/fWTTfdVOT9YvXq1br55pu1YsUKbdmyRb169VK/fv20bds2B0dufyU9N3nS0tL08MMP68Ybb3RQpOYozfnp37+/fv75Z82fP1979uzRkiVL1Lx5cwdG7RglPTdr167Vww8/rOHDh2vXrl1atmyZYmJiNGLECAdH7iQMlIok44svvsj/OTc31wgMDDRee+21/G2ZmZmGr6+vMXv27CL3079/f6NPnz4Ftt1yyy3GAw88YPOYbclW/f+77Oxsw9vb21i4cKEtw7U5W/Y/Ozvb6NKlizFv3jxjyJAhxp133mmnqG3LVudg1qxZRuPGjY2LFy/aM1ybs1X/o6KijBtuuKHAtqeeesro2rWrzWO2pb/3/68SExMNSca2bduuuR9nvQbCsc6fP2+4ubkZ33zzTYHtYWFhxoQJE4q9n5CQEGPSpEm2Ds9UZTk3AwYMMJ577jlj4sSJRlhYmB2jNE9pzs93331n+Pr6GikpKY4I0TSlOTdvvPGG0bhx4wLbpk+fbtSrV89ucTozRrRsJDExUSdOnFDv3r3zt3l6eqpHjx5av359kZ/bsGFDgc9I0i233HLVz5RHpe3/350/f16XLl2Sn5+fPcK0m7L0/6WXXlLt2rU1fPhwe4dpV6U9B//973/VqVMnRUVFKSAgQKGhofrPf/6jnJwcR4RtM6Xtf9euXbVly5b8KbMHDhzQihUrdNttt9k95vLAVa6BsK/s7Gzl5OTIy8urwPbKlStr7dq1xdpHbm6uMjIynO7+ci2lPTcffPCBEhISNHHiRHuHaKrSnJ///ve/ateunaZMmaK6devq+uuv19NPP60LFy44ImSHKc256dy5s44ePaoVK1bIMAydPHlSn332WYW5Z5WUu9kBuIoTJ05IkgICAgpsDwgI0KFDh676ucI+k7c/Z1Ha/v/d+PHjVbduXd100002jc/eStv/devWaf78+cV+lqU8K+05OHDggH755RcNHDhQK1as0L59+xQVFaXs7Gy98MILdo3Zlkrb/wceeEB//vmnunbtKsMwlJ2drVGjRmn8+PF2jbe8cJVrIOzL29tbnTp10ssvv6wWLVooICBAS5Ys0caNG9W0adNi7eOtt97SuXPn1L9/fztH61ilOTf79u3T+PHjtWbNGrm7u/Y/BUtzfg4cOKC1a9fKy8tLX3zxhU6dOqXRo0crNTXVpZ7TKs256dy5sxYvXqwBAwYoMzNT2dnZuuOOO/Tuu+86OHrnwIiWjVkslgI/G4ZxxTZbfKa8KktfpkyZoiVLlujzzz+/4tsVZ1GS/mdkZGjQoEGaO3euatWq5YjwHKKkvwO5ubny9/fXnDlzFBERoQceeEATJkzQrFmz7B2qXZS0/7/99pteffVVzZw5U1u3btXnn3+ub775Ri+//LK9Qy03XOkaCPv56KOPZBiG6tatK09PT02fPl0PPfSQ3NzcrvnZJUuW6MUXX9TSpUvl7+/vgGgdqyTnJicnRw899JAmTZqk66+/3oRoHa+kvzu5ubmyWCxavHix2rdvr1tvvVVTp07VggULXG5Uq6TnJj4+Xk888YReeOEFbdmyRd9//70SExM1cuRIB0fuHFz7awwHCgwMlHT529mgoKD87cnJyVd8W/v3z/39m9trfaY8Km3/87z55pv6z3/+o59++kmtW7e2W5z2Upr+JyQk6ODBg+rXr1/+ttzcXEmSu7u79uzZoyZNmtgxatsq7e9AUFCQKlWqVOCi3qJFC504cUIXL16Uh4eH/YK2odL2//nnn9fgwYPzHyRu1aqVzp07p8cee0wTJkyQ1era34e5yjUQ9tekSROtWrVK586dU3p6uoKCgjRgwAAFBwdf9XNLly7V8OHDtWzZMqebLVFcJTk3GRkZ2rx5s7Zt26bHH39c0uV7j2EYcnd3148//qgbbrjB0V2wq5L+7gQFBalu3bry9fXN39aiRQsZhqGjR48WexTVGZT03EyePFldunTRv/71L0lS69atVbVqVXXr1k2vvPJKgfsfGNGymeDgYAUGBmrlypX52y5evKhVq1apc+fORX6uU6dOBT4jST/++ONVP1Melbb/kvTGG2/o5Zdf1vfff6927drZO1S7KE3/mzdvrh07dig2Njb/dccdd6hXr16KjY1V/fr1HRW+TZT2d6BLly7av39/fpIpSXv37lVQUJDTJFlS6ft//vz5K5IpNzc3GYYhwzDsFm954SrXQDhO1apVFRQUpNOnT+uHH37QnXfeWWTbJUuWaOjQofr4448rxDMkxTk3Pj4+V9x7Ro4cqWbNmik2NlYdOnQwIXLHKO7vTpcuXXT8+HGdPXs2f9vevXtltVpVr149R4XrUMU9N0XdsyRViHtWiTm6+oYzy8jIMLZt22Zs27bNkGRMnTrV2LZtm3Ho0CHDMAzjtddeM3x9fY3PP//c2LFjh/Hggw8aQUFBRnp6ev4+Bg8ebIwfPz7/53Xr1hlubm7Ga6+9Zuzevdt47bXXDHd3d+P33393eP+uxR79f/311w0PDw/js88+M5KSkvJfGRkZDu/ftdij/39X3qsO2uMcHD582KhWrZrx+OOPG3v27DG++eYbw9/f33jllVcc3r9rsUf/J06caHh7extLliwxDhw4YPz4449GkyZNjP79+zu8f9dyrf6npKQY27ZtM7799ltDkvHJJ58Y27ZtM5KSkvL34czXQJjr+++/N7777rv8vydhYWFG+/bt8yuWjh8/3hg8eHB++48//thwd3c3oqOjC9xfzpw5Y1YX7Kak5+bvXLnqoGGU/PxkZGQY9erVM+677z5j165dxqpVq4ymTZsaI0aMMKsLdlPSc/PBBx8Y7u7uxsyZM42EhARj7dq1Rrt27Yz27dub1YVyjUSrBH799VdD0hWvIUOGGIZxubzzxIkTjcDAQMPT09Po3r27sWPHjgL76NGjR377PMuWLTOaNWtmVKpUyWjevLmxfPlyB/WoZOzR/4YNGxa6z4kTJzquY8Vkrz//vyrviZa9zsH69euNDh06GJ6enkbjxo2NV1991cjOznZQr4rPHv2/dOmS8eKLLxpNmjQxvLy8jPr16xujR482Tp8+7biOFdO1+v/BBx9c8++zM18DYa6lS5cajRs3Njw8PIzAwEAjKiqqQNI0ZMgQo0ePHvk/9+jR46q/r66kpOfm71w90SrN+dm9e7dx0003GZUrVzbq1atnPPXUU8b58+cdHLn9lebcTJ8+3QgJCTEqV65sBAUFGQMHDjSOHj3q4Midg8UwGOcDAAAAAFviGS0AAAAAsDESLQAAAACwMRItAAAAALAxEi0AAAAAsDESLQAAAACwMRItAAAAALAxEi0AAAAAsDESLQAAAACwMRItAAAAALAxEi3ABAcPHpTFYtHQoUPNDgUAAO5LgB2QaAEVWOvWrWWxWGSxWLRmzZpC25w+fVq1atXKb7d3714HRwkAqCi4L8GVkGgBFVRmZqZ2794td3d3SdKOHTsKbffCCy/o9OnTkiRvb281bdrUYTECACoO7ktwNSRaQAUVFxen7Oxs9e3bV1artdAb2s6dOzVr1izdeuutkqTw8HBZLBZHhwoAqAC4L8HVkGgB5czChQvVsWNHVatWTdWqVVPHjh21cOHCQttmZ2dr8uTJatKkiby8vHTddddp8uTJOnDgwDXn2m/dulWS1L17dzVp0qTQG9rYsWPl7e2tu+++W5IUERFR9g4CAJwK9yWgdNzNDgDA/4wbN07Tpk1T3bp1NXz4cFksFi1fvlxDhw5VXFycpk6dWqD9sGHD9NFHH6lJkyaKiopSVlaWpk2bpg0bNlzzWFu2bJEktW3bVuHh4frhhx8KvP/555/r559/1jvvvKNDhw7ltwUAVBzcl4AyMAA4XGJioiHJGDJkSP621atXG5KMFi1aGGfOnMnffubMGaN58+aGJGPNmjX523/66SdDktGuXTvj/Pnz+duTkpKMwMDAK/b/d23atDEkGadPnzZeffVVQ5Jx6NAhwzAMIzMz0wgODjZatGhhXLp0ybj99tsNScbOnTttdxIAAOUG9yXA9pg6CJQTCxYskCS9+OKL8vX1zd/u6+uriRMnFmgjSYsWLZIkPf/886pcuXL+9sDAQD355JNXPdbFixe1c+dONW7cWNWrV1d4eLik/z14/OabbyoxMVHTpk2Tu7u7tmzZoipVqqh58+YF9jNy5EhKAQOAi3LG+xJQnpBoAeXEtm3bJEk9e/a84r28bbGxsfnb4uLiJEmdO3e+on1h2/5q+/btunTpUv6Ui7CwMEmXb2jHjh3T5MmTdccdd6h37946ceKEkpKSFBYWJjc3tyv207p162L1DwDgXJzxvgSUJyRaQDmRnp4uq9Wq2rVrX/FeQECArFar0tLSrmhfs2bNQttfTd4Dx3kPEdetW1e1a9fWjh079Mwzz+jSpUt66623JBWcM/9XhmFox44d+TdDAIBrcbb7ElDekGgB5YSPj49yc3P1559/XvFecnKycnNz5ePjc0X7lJSUK9qfPHnyqscq7CYVFham7777Th9//LHGjRun6667TtL/bn5/v6ElJCTo7NmzunDhgrp27aoqVaooIiJCBw4cKGaPAQDlmbPdl4DyhkQLKCfatGkjSfrtt9+ueG/VqlWSlD9nXfrftIr169df0b6wbX9V2E0qPDxcp0+fVmBgoCZMmJC/Pe/m9/cSunFxcXJzc9OMGTM0ZcoUxcTEKDs7W1OmTLnqsQEAzsHZ7kuFMQxDX3/99TXbAfZAogWUE0OGDJEkTZo0Senp6fnb09PTNWnSpAJtJGngwIGSpJdfflmZmZn520+cOKF33nmnyONcunRJO3bsUIMGDVSrVq387aNHj9YXX3yhH3/8Ud7e3vnbt27dKk9PT4WEhBTYT1xcnAICAvTpp5+qc+fOatmypW666SadOnWqNN0HAJQzznZf+vs+T506pZ9//llvvvmmTp06pTNnzhSz54BtsI4WUE50795dY8aM0bvvvqvQ0FDde++9MgxDn3/+uY4cOaInnnhC3bt3z29/0003aeDAgVq8eLFatWqlO++8U1lZWfr000/VoUMHff3117Jar/wuZdeuXcrKyrpiykVwcLCCg4MLbDt16pSOHDmidu3aqVKlSgXei4uL04ABAwpMGzl48GD+1A4AgHNztvvS302ZMkVvvvmmDMNQx44d9cknn6hdu3ZlOCNAyTCiBZQj06dP1/vvv6/AwEDNmTNHc+fOVWBgoN5///1Cvw1csGCBXn75ZeXk5Ojdd9/VihUrNHbsWD333HOSVCAJylOSh4iv1jYuLk4dOnQosC02NpbiGADgQpzpvvRXlSpV0pQpU9SwYUP5+Pho/PjxJFlwOIthGIbZQQCwrXnz5unRRx/VzJkzNWrUKJvvPy0tTdWrV9fu3bvz1zBJT09X9erVtX37doWGhtr8mAAA52Xv+1Jh9u3bp82bN6t379769NNPHXZcIA+JFuDETpw4oYCAAFkslvxtx44dU5cuXXT06FElJiaqfv36Nj/u6tWr1adPH2VkZOSvYbJmzRrdfPPNOnv2rNzdmZUMABWRWfcloDziX0OAE3vttdf07bffqlu3bvL399fhw4f1zTffKCMjQy+++KLdbmZxcXFq2bJlgYUi4+Li1KJFC5IsAKjAzLovAeURI1qAE/v+++81depUxcXF6fTp0/Ly8lLr1q01evRoPfTQQ2aHBwCoYLgvAf9DogUAAAAANkbVQQAAAACwMRItAAAAALAxEi0AAAAAsDESLQAAAACwMRItAAAAALAxEi0AAAAAsDESLQAAAACwMRItAAAAALAxEi0AAAAAsDESLQAAAACwMRItAAAAALCx/w9gzHNdCgK+HAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def plot_hmf_and_smf(smf, plotarg=\"C0o-\", label=None, axes=None):\n",
    "    if axes is None:\n",
    "        _, axes = plt.subplots(ncols=2, figsize=(10, 4))\n",
    "    axes[0].hist(log_halo_masses, bins=50, color=\"C0\")\n",
    "    axes[0].semilogy()\n",
    "    axes[0].set_xlabel(\"$\\\\log M_h$\", fontsize=14)\n",
    "    axes[0].set_ylabel(\"$\\\\Phi(\\\\log M_h)$\", fontsize=14)\n",
    "\n",
    "    smf_bin_cens = 0.5 * (smf_bin_edges[:-1] + smf_bin_edges[1:])\n",
    "    axes[1].semilogy(smf_bin_cens, smf, plotarg, label=label)\n",
    "    axes[1].legend(frameon=False)\n",
    "    axes[1].set_xlabel(\"$\\\\log M_\\\\star$\", fontsize=14)\n",
    "    axes[1].set_ylabel(\"$\\\\Phi(\\\\log M_\\\\star)$\", fontsize=14)\n",
    "    return axes\n",
    "\n",
    "plot_hmf_and_smf(true_smf, label=\"truth\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introducing MultiGrad\n",
    "\n",
    "To fit such a model with multigrad, we have to define a class that inherits from an abstract base class called `OnePointModel`. To make this class useable, you must define two methods: `calc_particle_sumstats_from_params` (which simply calls the `calc_smf` function we already defined above) and `calc_loss_from_sumstats` (which in this example simply computes a mean-squared error of the log of the stellar mass function in each bin). **Important:** The sumstats that we return **must** be summable over all ranks, so be careful not to return log quantities here. You can always take the logarithm of your sumstats in the loss function!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "class MySMFModel(multigrad.OnePointModel):\n",
    "    def calc_partial_sumstats_from_params(self, params):\n",
    "        # Accessing global variables is fine, but I prefer to store them in\n",
    "        # the `aux_data` attribute, which we will define during construction\n",
    "        bin_edges = self.aux_data[\"smf_bin_edges\"]\n",
    "        volume = self.aux_data[\"volume\"]\n",
    "        log_halo_masses = self.aux_data[\"log_halo_masses\"]\n",
    "\n",
    "        return calc_smf(params, bin_edges, volume, log_halo_masses)\n",
    "\n",
    "    def calc_loss_from_sumstats(self, sumstats):\n",
    "        # Add 1e-10 so that log values always remain finite\n",
    "        target_sumstats = jnp.log10(self.aux_data[\"target_sumstats\"] + 1e-10)\n",
    "        sumstats = jnp.log10(sumstats + 1e-10)\n",
    "        # Reduced chi2 loss function assuming unit errors (mean squared error)\n",
    "        return jnp.mean((sumstats - target_sumstats)**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "aux_data = dict(\n",
    "    log_halo_masses=log_halo_masses,\n",
    "    smf_bin_edges=smf_bin_edges,\n",
    "    volume=volume,\n",
    "    target_sumstats=true_smf  # SMF at truth: params=(-2.0, 0.2)\n",
    ")\n",
    "model = MySMFModel(aux_data=aux_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "At true_params:\n",
      "Loss = 0.0\n",
      "Grad(loss) = [0. 0.]\n",
      "At true_params + [0.1, 0.1]:\n",
      "Loss = 0.44032094\n",
      "Grad(loss) = [2.6187496 4.2603974]\n"
     ]
    }
   ],
   "source": [
    "# Since this notebook kernel is only a single MPI rank, the partial sumstats\n",
    "# are identical to the total sumstats:\n",
    "assert np.all(model.calc_partial_sumstats_from_params(true_params) ==\n",
    "              model.calc_sumstats_from_params(true_params))\n",
    "\n",
    "# Using the `calc_loss_and_grad_from_params` method:\n",
    "\n",
    "# At the true set of parameters, the loss is at a global minimum of 0 :)\n",
    "loss, gradloss = model.calc_loss_and_grad_from_params(true_params)\n",
    "print(\"At true_params:\")\n",
    "print(\"Loss =\", loss)\n",
    "print(\"Grad(loss) =\", gradloss)\n",
    "\n",
    "# Perturbing these parameters yields a descendable gradient\n",
    "loss, gradloss = model.calc_loss_and_grad_from_params(\n",
    "    true_params + jnp.array([0.1, 0.1]))\n",
    "print(\"At true_params + [0.1, 0.1]:\")\n",
    "print(\"Loss =\", loss)\n",
    "print(\"Grad(loss) =\", gradloss)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Recovering the truth with gradient descent\n",
    "\n",
    "MultiGrad provides several built-in gradient descent methods, notably Adam and BFGS. Adam is particularly useful for stochastic problems, but BFGS will likely be more powerful in this case. Running it is as simple as..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "BFGS Gradient Descent Progress:  15%|█▌        | 15/100 [00:00<00:05, 16.69it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BGFS has converged: True\n",
      "Initial guess = [-3.5         0.19999999]\n",
      "True params = [-2.  -0.5]\n",
      "Converged params = [-1.99999865 -0.50000094]\n",
      "\n",
      "Full results info:\n",
      "  message: CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH\n",
      "  success: True\n",
      "   status: 0\n",
      "      fun: 5.930545117494024e-12\n",
      "        x: [-2.000e+00 -5.000e-01]\n",
      "      nit: 15\n",
      "      jac: [-1.169e-05 -2.370e-05]\n",
      "     nfev: 31\n",
      "     njev: 31\n",
      " hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "# Initial guess for our parameters. If it's too far off, there is always a\n",
    "# risk of getting stuck in local minima or other zero-valued gradients\n",
    "init_params = true_params + jnp.array([-1.5, 0.7])\n",
    "\n",
    "# Run gradient descent using the BFGS method powered by scipy\n",
    "params, loss, results = model.run_bfgs(init_params)\n",
    "print(\"BGFS has converged:\", results.success, flush=True)\n",
    "print(\"Initial guess =\", init_params, flush=True)\n",
    "print(\"True params =\", true_params, flush=True)\n",
    "print(\"Converged params =\", results.x, flush=True)\n",
    "print(\"\\nFull results info:\", flush=True)\n",
    "print(results, flush=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1oAAAF8CAYAAADWyMqcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAB3yklEQVR4nO3deVxN+f8H8Ne97ahIWiiViBKlTdn3bewMhkEYM0ZjHWP4MpZhhsEYjDB2M9axz9i3sRYllSVLpUWUVLRRWs7vD7/uuFNpu3W6t9fz8bgP7jmfe877c+Sc3vfzOe8jEQRBABERERERESmMVOwAiIiIiIiIVA0TLSIiIiIiIgVjokVERERERKRgTLSIiIiIiIgUjIkWERERERGRgjHRIiIiIiIiUjAmWkRERERERArGRIuIiIiIiEjB1MUOoLLLzc3Fs2fPoKurC4lEInY4RERViiAISE1NRd26dSGV8rvBPLw2ERGJoyTXJSZaRXj27BnMzc3FDoOIqEp78uQJzMzMxA6j0uC1iYhIXMW5LjHRKoKuri6AdwdTT09P5GiIiKqWlJQUmJuby87F9A6vTURE4ijJdYmJVhHypmTo6enxYkZEJBJOj5PHaxMRkbiKc13ihHciIiIiIiIFY6JFRERERESkYEy0iIiIiIiIFIyJFhERERERkYIx0SIiIiIiIlIwJlpEREREREQKxvLuREREVUhOrgC/iCTEp2bASFcbblYGUJOyfD4RkaIx0SIiIqoiTt2NxcK/QxCbnCFbZqqvjfl97NDD3lTEyIiIVA+nDhIREVUBp+7G4sudt+SSLACIS87Alztv4dTdWJEiIyJSTUy0iIiIVFxOroCFf4dAeG9Z9czXACBbtvDvEOTkCvk+S0REpcNEi4iISMX5RSTJjWRJc3NwaOcMbN2/APZxYRAAxCZnwC8iSbwgiVRATq4A3/BEHA16Ct/wRJX58iIyMhISiQRBQUFih6JUeI9WIby9veHt7Y2cnJwyb8ty1vFC10Uu/ajM2yciIvqQ+FT56YIOsaGwToxB44RodHp8E6cbueOXNiMQn+ooToBEKqCi74GUSD5cxGb06NHYvn17ibfr6emJV69e4ciRI6ULjGQ4olUILy8vhISEwN/fX+xQiIiIysRIV1vufWC9Jug8fgMONu2IHIkU3UOv49S2STD6bDSirt4UKUoi5SXGPZCxsbGy16pVq6Cnpye3bPXq1XLts7KyFB4DfRgTLSIiIiXh7e0NOzs7uLq6luhzblYGMNXXxvvff0fVqouve3+NbmO98XeTtgAAj4DzMGvXEsuX7UNkQroCIydSLoIg4PXb7GK9UjOyMP+veyhokmDesgV/hSA1I6vIbQlC8acampiYyF76+vqQSCSy9xkZGahZsyb+/PNPdOjQAdra2ti5cycWLFgAR0dHue2sWrUKlpaW7+JcsAA7duzA0aNHIZFIIJFIcPHiRVnbx48fo2PHjqhWrRocHBzg6+tbouNa1XDqIBERkZLw8vKCl5cXUlJSoK+vX+zPqUklmN/HDl/uvAUJIPcL4WNDc0zu9y2yvp0Fq1+XIzcuDt6J1bFh5SUMcqqHye71YGZmqPC+EFVmb7JyYDfvtEK2JQCIS8lAswVnimwb8n13VNNU3K/n3377LX7++Wds27YNWlpa2Lhx4wfbz5gxA/fv30dKSgq2bdsGADAwMMCzZ88AAHPmzMGKFSvQqFEjzJkzB5988gnCwsKgrs6UoiAc0SIiIqoCetibYv2nTjDRl59GaKKvjfWfOmGgZy+0CPgHWhf/QccmRsjJFXDq6gPoNG6IGz2G4HlIqEiRE1FpTZ06FQMHDoSVlRXq1q1bZPsaNWpAR0cHWlpastExTU1N2foZM2bgo48+go2NDRYuXIioqCiEhYWVZxeUGtNPIiKiKqKHvSm62pnALyIJ8akZMNLVhpuVAdSk/04qtG9kgm2NTBAQ9RIB839G7dfJqH16PzKbH8H1nkPQ8JfFMGxoKV4niCqAjoYaQr7vXqy2fhFJ8NxW9D3928e4ws3KoMj9KpKLi4tCt9e8eXPZ301N3xX4iI+PR5MmTRS6H1XBREtkhVUkZDVCIiIqD2pSCTysaxfZztmiFpy3L0ZIj5bAvHmwCw2C+7FdeHPqAK73HQ6bFd/DwMqsAiImqngSiaTYU/jaNqoDU31txCVnFHiflgTvRo7bNqoj96VGRahevbrce6lUmu8+sJIUydDQ0JD9Pa/qYW5ubhkiVG2cOkhERESFshvWB7YPAnBn2wE8tGoKnexMuB/aBnXbJlh9yB+vXr8VO0QiUeXdAwkA/02j8t7P72NX4UlWQerUqYO4uDi5ZOu/z8bS1NRUyOONiIkWERERFUEilaKZ5yDYhN1G8IZdCDVvjNONPPCLXzza/vQPVp17hJTkNLHDJBJNUfdAlsdztEqjQ4cOePHiBZYtW4bw8HB4e3vj5MmTcm0sLS1x+/ZtPHz4EAkJCSwLXwZMtIiIiKhYJFIpHL4YjoaRIai19Tc0MdFFamY2/t57Adl1zeA7fgbSE1+JHSaRKHrYm+Lqt52wZ7w7Vg9zxJ7x7rj6badKk2QBgK2tLdatWwdvb284ODjAz88PM2bMkGszfvx4NG7cGC4uLqhTpw6uXbsmUrTKTyKUpGB/FZRXQjc5ORl6enql2kZh92F9CO/RIiJSzDlYFVWW45KbK+DE3VhkfDUFg68cAAAkVdfHo9ET4fDDbOjU1BUtNiKi8lCS8y8TrSKIlWh9CJMwIqoqKktCUdlUtuOSk5WNwJ/WwWTVMpglPgUAJNQwQNhnk9Bi8UxoVa8mcoRERIpRkvMvpw4SERFRmahpqMNl7mSYPH0M/7nLEVvTGIZpSXBftRBRjZpj1/VIvM1mZTIiqlqYaBEREZFCqGtpwnXRDNR+Gokb3yxGvJ4hDtm0xZwj99Dp54v40z8a2ZnyVQpzcgX4hifiaNBT+IYnIieXE22ISDXwOVpERESkUJrVtNFy2RxkfDcVpn5PYOgbg5iXb3Bm6WZ4XN6G51NnosXML3H24Qss/DsEsckZss+a6mtjfh+7SlVAgIioNDiiRUREROVCW7c6RndugiszO2J2zyYYH/Q3zBNi4DJ3MqLMGuKvuWsQ9+q13GfikjPw5c5bOHU3VqSoiYgUg4kWERERlSsdTTV80d4aTW9cgO/YaUjRroEG8VFYd3Qpjm+fgq6h14H/r82VN3Fw4d8hnEZIREqNVQeLUBmrDhaG1QiJSNVUtup6lYWyH5d//EIRPH0+xvkfge7bNwCAg0074uveX8u12zPeHR7WtcUIkYioQKw6SERERJVWimY1rGozAm0nbME698F4raGF403a5msXn5pRwKeJiJQDEy0iIiKqUEa62gCAVzp6WNbeE22/2IIL1q6y9R/fPoORt46hthZ/TSFSNIlEgiNHjnywjaenJ/r371/sbUZGRkIikSAoKKhMsakaVh0kIiKiCuVmZQBTfW3EJWdAAJBYvaZsXa3XyZh7YQv0M9PxuPNZ3Pt5JZp+0le0WIkqM09PT7x69arIxOl9sbGxqFWrFoB3CZKVlRUCAwPh6Ogoa7N69Wrw7qKy41dFREREVKHUpBLM72MHAJD8Z12qdg0sbzcKydo10CA2HE2H98PNVj0Qdze04gMlUkEmJibQ0tL6YBt9fX3UrFmzYgJSYUy0iIiIqML1sDfF+k+dYKKvLbe8Tq3qaLPyOwgPH+FG94+RCwlcfE9Dz6kZfMdNR0ZqukgRE1V+HTp0wOTJkzFz5kwYGBjAxMQECxYskGvz/tRBKysrAECLFi0gkUjQoUMHAPmnDp46dQpt2rRBzZo1Ubt2bfTu3Rvh4eElii02NhYfffQRdHR0YGVlhd27d8PS0hKrVq0CUPD0w1evXkEikeDixYuyZSEhIejVqxdq1KgBY2NjjBw5EgkJCbL1Bw4cQLNmzaCjo4PatWujS5cuSE9/d964ePEi3NzcUL16ddSsWROtW7dGVFRUifpREpw6qEIKq27IaoRERFQZ9bA3RVc7E/hFJCE+NQNGutpwszKAmvTdOFfLU38i7PRlZHlNgm34bbhtW41xek0xYnxvdLE1gkTy3/EwIsXK+wW9IGpqatDW1i5WW6lUCh0dnQ+2rV69eimjlLdjxw5Mnz4dN27cgK+vLzw9PdG6dWt07do1X1s/Pz+4ubnh3LlzaNq0KTQ1NQvcZnp6OqZPn45mzZohPT0d8+bNw4ABAxAUFASptHjjNqNGjUJCQgIuXrwIDQ0NTJ8+HfHx8SXqW2xsLNq3b4/x48dj5cqVePPmDb799lsMGTIEFy5cQGxsLD755BMsW7YMAwYMQGpqKq5cuQJBEJCdnY3+/ftj/Pjx2LNnD96+fQs/P79yPY8w0SIiIiLRqEklHyzh3rB7OwiPAnHzp/W4dd4PF7VMcPH3m2hnUwcLWhmjQROLCoyWqpoaNWoUuq5Xr144fvzfL7mNjIzw+vXrAtu2b99eblTG0tJSbhQGgMLuiWrevDnmz58PAGjUqBHWrl2L8+fPF5ho1alTBwBQu3ZtmJiYFLrNQYMGyb3fsmULjIyMEBISAnt7+yJjevDgAc6dOwd/f3+4uLgAADZv3oxGjRoVu18AsH79ejg5OeHHH3+ULdu6dSvMzc3x6NEjpKWlITs7GwMHDoSFxbtzQ7NmzQAASUlJSE5ORu/evWFtbQ0AsLW1LdH+S4pTB4mIiKhSk0ilcJnthRHHt+DLDtbQVJPimU8ATJs1xvVBY5Eanyh2iESVRvPmzeXem5qalnjk6L/Cw8MxfPhwNGjQAHp6erIph9HR0cX6/MOHD6Gurg4nJyfZsoYNG8qKchRXQEAA/vnnH9SoUUP2atKkiSxGBwcHdO7cGc2aNcPHH3+MTZs24eXLlwAAAwMDeHp6onv37ujTpw9Wr16N2NjYEu2/pDiiRUREREqhupY6vu3RBENczHH7s8PQyc6E+6FtSDhzBA9mfAfnOZMhVVcTO0xSIWlpaYWuU1OT/1n7UDLz3+l1kZGRZYrrQzQ0NOTeSyQS5Obmlmmbffr0gbm5OTZt2oS6desiNzcX9vb2ePv2bbE+X9ho3fvL847R+8uysrLk2ufm5qJPnz746aef8m3L1NQUampqOHv2LHx8fHDmzBn8+uuvmDNnDm7cuAErKyts27YNkydPxqlTp7Bv3z7MnTsXZ8+ehbu7e7H6UVIc0SIiIiKlYmVYHf2ObETwhl14YmgGw7SXcF0wHaE2Dgg9/o/Y4ZEKqV69eqGv9+/PKqrt+/dnFdZWDHn3ZOXk5BTaJjExEffv38fcuXPRuXNn2NraykaJiqtJkybIzs5GYGCgbFlYWBhevXole583jfH9Uab/PpfLyckJ9+7dg6WlJRo2bCj3yjuGEokErVu3xsKFCxEYGAhNTU0cPnxYto0WLVpg9uzZ8PHxgb29PXbv3l2ivpQEE61CeHt7w87ODq6urkU3JiIiKoHU1FS4urrC0dERzZo1w6ZNm8QOSSk5fDEcRlGPcP2LmUjX1EHjiHuw7t0ZV3qPREJaptjhEVV6RkZG0NHRwalTp/D8+XMkJyfna1OrVi3Url0bGzduRFhYGC5cuIDp06eXaD9NmjRBly5d8Pnnn8PPzw+BgYH4/PPPoaOjIytGoaOjA3d3dyxduhQhISG4fPky5s6dK7cdLy8vJCUl4ZNPPoGfnx8eP36MM2fOYOzYscjJycGNGzfw448/4ubNm4iOjsahQ4fw4sUL2NraIiIiArNnz4avry+ioqJw5swZPHr0qFzv02KiVQgvLy+EhITA399f7FDKzHLW8UJfRERU8apVq4ZLly4hKCgIN27cwJIlS5CYyPuMSkOrmg7cN/yE17fv4mabjyCFgEuvNdFxxUVsvRqBrJyyTZkiUmXq6upYs2YNfvvtN9StWxf9+vXL10YqlWLv3r0ICAiAvb09pk2bhuXLl5d4X7///juMjY3Rrl07DBgwAOPHj4eurq7cyODWrVuRlZUFFxcXTJkyBYsXL5bbRt26dXHt2jXk5OSge/fusLe3x5QpU6Cvrw+pVAo9PT1cvnwZvXr1go2NDebOnYuff/4ZPXv2RLVq1fDgwQMMGjQINjY2+Pzzz/HVV1/hiy++KPmBKyaJwMc+f1BKSgr09fWRnJwMPT29Um2jMic0LP1ORJWZIs7BlV1SUhJatGiBgIAAGBoaFuszVeG4lNb9w2cwK0IdwfFvAAAD0x9jrJsZ7EcOEDkyInpfTEwMzM3Nce7cOXTu3FnscIqtJOdfjmgRERGV0OXLl9GnTx/UrVtX7uGf71u3bh2srKygra0NZ2dnXLlyRW79q1ev4ODgADMzM8ycObPYSRZ9mO2Abjg0tSN+HNAMdTQBr93LYD9qIG65dUbs7Qdih0dUZV24cAF//fUXIiIi4OPjg2HDhsHS0hLt2rUTO7Ryw0SLiIiohNLT0+Hg4IC1a9cWuH7fvn2YOnUq5syZg8DAQLRt2xY9e/aUK4Vcs2ZNBAcHIyIiArt378bz588rKnyVpyaVYHjL+jj3VSskeLRDjkQKJ/8LqOXsAN/Rk5GRUnglOSIqH1lZWfjf//6Hpk2bYsCAAahTp47s4cWqilMHi8Cpg0RE4lGGKXISiQSHDx9G//79ZctatmwJJycnrF+/XrbM1tYW/fv3x5IlS/Jt48svv0SnTp3w8ccfF7iPzMxMZGb+W9whJSUF5ubmlfq4VCYRF3zwesJXaBr6ruJZbE1jxM3/AY6Tx0Ai5XfORFR8nDpIREQkkrdv3yIgIADdunWTW96tWzf4+PgAAJ4/f46UlBQA7y7aly9fRuPGjQvd5pIlS6Cvry97mZubl18HVJBVp1awe3ATAcvWI06/DkxfPUeLaZ9h8bcb8Oh5qqxdTq4A3/BEHA16Ct/wROTk8rtoIio9PrC4iitstI0jXUREpZOQkICcnBwYGxvLLTc2NkZcXByAdzeBjxs3DoIgQBAEfPXVV2jevHmh25w9e7ZcOeW8Ea2yGDFiBDQ0NGBmZgZzc3O5l76+vqzksqqQSKVw/mYCXn/2CXy/+h9e+wdgi7Q+tq++gtEelmhqWgMrzoYiNjlD9hlTfW3M72OHHvamIkZORMqKiRYREVE5+G+iIgiCbJmzs3O+B3F+iJaWFrS0tBQWW25uLvbv34+srKwC13fs2BEXLlyQvV+5ciX09PTkkjFdXV2FxVORqtXSh8cub0QnpKPbifs4E/Ich84F45Nd36Kt6wDsb94FguTdhJ+45Ax8ufMW1n/qxGSLiEqMiRYREZECGRoaQk1NTTZ6lSc+Pj7fKJdYcnNzsXXrVjx58iTfKykpCbVr15ZrO2vWrHxJWd4Uxm7duuHnn3+WLb927Rrq1KkDc3Nz6OjoVFifSqq+YXVsHOWCiw/iETx+GholPsGyU2swIugkFnT5AoH1mkAAIAGw8O8QdLUzgZpUtUb5iKh8MdEiIiJSIE1NTTg7O+Ps2bMYMODfZzedPXu2wIeBikFdXR2ffvppgetev36N169fy95nZmZizJgxcslYcnKy7GVraytrm5ubi44dO8qSstq1a8tGwMzMzODh4YGRI0fK2r99+xaamprl1Mvi0dJQw6/uQ5CioYOp13bDIS4UB3d+g+3OfbCs/ShkaGgjNjkDfhFJ8LCuXfQGiYj+HxMtIiKiEkpLS0NYWJjsfUREBIKCgmBgYID69etj+vTpGDlyJFxcXODh4YGNGzciOjoaEyZMKNN+vb294e3tjZycnLJ2oVDVqlVDtWrVZO91dHTw22+/ybVJTU1FTEwMnjx5An19fdnylJQUWFtb48mTJ0hPT0diYiISExNl0yQTEhJkiVZubi50dXVRs2ZNmJubo379+nIvW1tbNG3atNz6mSc+NQPZaurY4jYAf9l1wLeXtmPw3fMYG/AXOjy+iW96TUWAmR3iUzOK3hgR0XtY3r0Iql7evTRYKIOIKkplLe9+8eJFdOzYMd/y0aNHY/v27QDePbB42bJliI2Nhb29PX755ReFPZizsh6XPIIg4NWrV7JkLO/VtGlTfPLJJwCA2NhY1K1bt9BtfPzxx/jzzz8BvEvKevTogXr16uVLyMzNzeUSw5LyDU/EJ5uuyy3rEH4TS079CtO0ROy374JvPpqKHWNd0d7GqNT7IaoKFixYgPXr1yM+Ph6HDx/GkSNH8OrVqwIf6q6sSnL+ZaJVBCZa+THRIqKKUtkTCrGownERBAGJiYlyiVh0dLTs1atXL8yZMwdA0UnZmDFjsHXrVgDvkrJVq1bJJWNGRkaQFvK8rJxcAW1+uoC45Ay8/wuRXkYaplzbg9WtP0GKdg1Y1a6GFYObw9mK0wep8vD09MSOHTtk7w0MDODq6oply5Z9sJJpSSxYsABHjhwpsoDP/fv3YWdnh8OHD8Pd3R21atVCRkYGBEFAzZo1AQAdOnSAo6MjVq1apZDYxFCS8y+nDhIREVGFk0gkMDQ0hKGhIVq0aPHBtjVq1MAff/whl4hFR0cjKioKaWlpMDAwkLV9/vw5vv76a7nPa2pqyqYnDh48GBMnTgTwLikLCw3Ft10sMe3gA0gAWbKVol0DizuPhwBAT1sdEQnpeNmtF6472MNx2xpo61ZX4NEgKr0ePXpg27ZtAIC4uDjMnTsXvXv3RnR0dIXGER4eDgDo16+frMKqIqulKiMmWkRERFSp6erqFli8QxAEJCcny92zlp2djWHDhsmSsWfPnuHt27cIDw9HeHg4nJycZG2fP3+OJk2aAAD0ahkgR6c2cqvXhrpeHajrGaFuYwes8BoMjwaG+GPZ7+gS5geE+SHqyjlkbNqCxn27lH/niYqgpaUFExMTAICJiQm+/fZbtGvXDi9evECdOnUAAE+fPsX06dNx5swZSKVStGnTBqtXr4alpSWAd9OhZ86ciXv37kFDQwNNmzbF7t278c8//2DhwoUA/n1kxbZt2+Dp6SkXw4IFC2Tt8kaPBUGAp6enbOqgp6cnLl26hEuXLmH16tUA3t3fmheDKmKiRUREpCQqohiGMpFIJLIpSXnMzc2xZ88e2fusrCw8e/ZMlnjZ2NjI1sXHx6NGjRpIS0tDyssk4GUSgFDZ+tGO02TPz/r4y94Y91td/PDiBSzio5HTrxuOuneEwQ9z0cjO9oPTE0mJpacXvk5NDdDWLl5bqRR4/3EHBbWtXvZR0rS0NOzatQsNGzaUPabh9evX6NixI9q2bYvLly9DXV0dixcvRo8ePXD79m1IpVL0798f48ePx549e/D27Vv4+flBIpFg6NChuHv3Lk6dOoVz584BgFwBnDwzZsyApaUlxowZg9jY2AJjW716NR49egR7e3t8//33ACBLBFUVEy0iIiIl4eXlBS8vL9k9AlQ0DQ0NWFhYwMLCIt86BwcHpKSkIDk5GVFRUXJTEqOiotDKw0PWNiIiAlvjnuEwgDUAPoWAftcv4G7nC+gNoOOMGVi+fDkA4OXLl1i3bp1sv5aWlqhbty7U1NQqptOkODVqFL6uVy/g+Hv34RsZAe89GkFO+/bAxYv/vre0BBIS5NuUsmzCsWPHUOP/40xPT4epqSmOHTsmS/z37t0LqVSKzZs3y41K1axZExcvXoSLiwuSk5PRu3dvWFtbA4DcYxtq1KgBdXV12ahZQWrUqCH70qOwdvr6+tDU1ES1atU+uC1VwkSLSqyw4h4skkFERMomb1SsZs2acHBwKLSdra0tTp8+jaioKDyIisK8U+fwVeBN2Ofm4A8Ac1I1kJmdAy11NTx8+BBz586V+7y6ujrMzMxgaWmJL774AsOGDQMAZGRk4OnTpzA3Nxf9mWKknDp27Ij169cDAJKSkrBu3Tr07NkTfn5+sLCwQEBAAMLCwqCrqyv3uYyMDISHh6Nbt27w9PRE9+7d0bVrV3Tp0gVDhgyBqampGN1RKUy0iIiIiIqgr6+Pbt26/btg8WK8jHqGm8PHYUn9drhZvTH6/noNPw9xgK6uLsaMGYPIyEhERUXhyZMnyMrKQmRkJCIjI+UeZB0cHAx3d3dIJBLUrVtXbhTMwsIC7dq1kxtdoAqWllb4uv+OUMbHF972v9NKIyNLHdJ/Va9eHQ0bNpS9d3Z2hr6+PjZt2oTFixcjNzcXzs7O2LVrV77P5k3d27ZtGyZPnoxTp05h3759mDt3Ls6ePQt3d3eFxVkVMdEiIiIiKoVaFnXhcu0kPrsTi8gjd/HweSoOjJmF7sZq2PDbOmhWe3f/Tk5ODmJjY2VTEl1cXGTbSEhIgLa2tmxk6+nTp/Dx8ZGtX7VqlSzRCgoKgpeXFywtLWFlZQUrKyvZ383NzaGhoVGxB6AqKMl9U+XVtoQkEgmkUinevHkDAHBycsK+fftgZGT0wXLkLVq0QIsWLTB79mx4eHhg9+7dcHd3h6ampsLuC1XktpQBEy0iIiKiMujZzBRuVgZYsf0ffPvzNuhkZyL8wmlg21ZYd2sLNTU1mJmZwczMDK1bt5b77EcffYTXr18jPj4eUVFRslGwvL+//yykhw8fwsfHRy4RyyOVSrFhwwaMHz8eABATE4Pz58/LEjLeI6a6MjMzERcXB+Dd/YFr165FWloa+vTpAwAYMWIEli9fjn79+uH777+HmZkZoqOjcejQIXzzzTfIysrCxo0b0bdvX9StWxcPHz7Eo0ePMGrUKACApaUlIiIiEBQUBDMzM+jq6pa6bLulpSVu3LiByMhI1KhRAwYGBipdRIaJFhERkZJg1cHKq3YNLSz5qgcC3qxCgwXfwvpZGLJ6dsT1TyfC+bfl0NAu/BdTiUQCY2NjGBsbw83NrdB2bdu2xZ9//omIiAhERkbK/oyMjERGRgaMjIxkbX18fORKcGtoaKB+/fqyxOuzzz6T7Ss7OxtqamqyQgmkXE6dOiW7n0pXVxdNmjTB/v370aFDBwBAtWrVcPnyZXz77bcYOHAgUlNTUa9ePXTu3Bl6enp48+YNHjx4gB07diAxMRGmpqb46quv8MUXXwAABg0ahEOHDqFjx4549epVgeXdi2vGjBkYPXo07Ozs8ObNG5Uv7y4RhFKWOKkiSvL058IUVjyiqmCRDCIqLUWcg1URj0vllhAejehhnnC6+Q8AIMysEdR2bIdVp1blsr/c3FzEx8dDV1cX1f9/StqpU6ewYsUK2QhZdna23GcOHz6M/v37AwD279+PUaNGwdLSssBpiU2aNJFVtSOq6kpy/uWIFhEREZECGVrXR+0b53Bz2QY0/H4WGsaE4nWPTti87zI8+7pAXU2xU6WkUmm+ctk9evRAjx49ALy7R+zp06eyUbCIiAg4OjrK2uaNiD148AAPHjzIt/1Dhw7JCnj4+vpi3759aNCggexlZWUFnfefEUVEAJhoERERESmcRCqFy6yJSBjcC4FDPXFRpx5W33iBozE+WPGxAxqb6Ba9EQVRU1ND/fr1Ub9+fbRr1y7f+ilTpmDgwIH5piTm/WllZSVre+3aNaxevTrfNkxNTdGgQQOsWrVKVuwjMTERGRkZMDU1Ven7cIgKw0SLiIiIqJwYNrREbf8LiLz1BHrHHuDO02TMnLsDsySRcF37I9S1xH92lqamJqytrWUPq/0QNzc3zJgxA48fP5a9UlJSEBsbi9jYWLmCG9u3b8eMGTOgpaUFKysruVGwBg0aoH379rKH3BKpIiZaREREROVIIpVigIsFWtkY47s/AzF12yTYxUfg0ZkT0Nq5AxZtXcUOsdjatWsnNyomCAJevnwpS7psbGxk6169egU1NTVkZmYWOC0xICAATk5OAIB9+/bh6NGj+ZKxevXqsVoiKS0mWkREREQVwFhPG7+NbQn/p5ORsmQubKLvI7Nja1z/bBpcf/0BahrK92uZRCKBgYEBDAwM5J4PBgCLFi3CvHnz8OTJEzx+/BgRERFyI2ENGjSQtb169Sr27NmTb/saGhqwtLTE8ePH0ahRIwBAWFgYUlNT0aBBA+jr65dvB4nKQPn+RxMREVVRLO+u/CRSKdzmT0X8kD6IGDYKDrd94P7bMjw8cxw6O3egfitnsUNUKA0NDdno1IcMGzYMFhYWcolYZGQksrKyEBoaitq1a8varlmzBr/++isAwMDAQDbtMe81ePBg6OpW3D1wRIVhefcisLx7+WLpdyL6EJYxLxiPi2oQcnPhv+AX2C6bB93M18hU08DxDQfQb2wfqEklyMkV4BeRhPjUDBjpasPNygBq0qrzrKucnBw8e/YMERERaNu2rew5XzNmzMAff/yB+Pj4Aj/34sULGBoaAgCWL1+Oq1evypKwhg0bwtraGhYWFtDQ0KiwvpDqYHl3IiIiokpOIpXC7fuvETekLyKGjcKb5FR8HSrBrt980cfBFL9deozY5AxZe1N9bczvY4ce9qYiRl1x1NTUYG5uDnNzc7nlK1aswIoVK5CWlobHjx8jPDxc9nr27Jnc6Nfly5dx7NixArddv359BAYGyqYf3r17F7m5uWjQoAGfG0YKwRGtInBEq3xxRIuIPoQjNwXjcVE9Qm4uDv5zDwsuP0NaZjZ03magX8hF7HPoBkHyrjR63ljW+k+dqkyyVVbXrl1DcHAwwsLCZMnY48eP8ebNG9SoUQMpKSmykbJ+/frhr7/+AgAYGxvnm5I4bNgwqKtzjKKq44gWERERkRKRSKUY3LkZXB0boOvKy5h5eQfGBPyNAff+wTe9piK6likEvEu2Fv4dgq52JlVqGmFptW7dGq1bt5ZbJgiCrBx9XpIFANWrV4eBgQGSkpLw/PlzPH/+HD4+PrJ1I0aMkLWdMWMGIiMj5RKxhg0bwszMjFUSSYaJFhEREVEl8exVBt7m5OKRoQXSNbTRMuYejm+fjG96TcWpxq0hAIhNzoBfRBI8rGsXuT3KTyKRoG7duqhbt67c8t27dwN4V5b+/emI4eHhss/lOXPmDO7cuZNv25qammjatCkCAgJk7W/fvg09PT2Ym5szCatimGiRqAqbVskphUREVBXFp767J2uPYw9csXTEyuMr4RYTgg1HlmCrc18s6TgGWWoasnakeDVr1oSzszOcnQuvALlixQrcv39fLhmLiIjA27dvkZWVJZeUjRo1CsHBwbIKjI0aNULDhg3RsGFDNGnSBJ07d66IbpEImGgVgiV0iYiIqKIZ6WrL/h5T0wSffLIE31z+HRNuHMTYgL/Q4tlDTOw/S64dVbxu3bqhW7ducstycnLw5MkTJCcnyy3X1NSEpqYm3r59i4cPH+Lhw4eydba2tggJCZG9nzJlCgDIkrFGjRrBwsKC94YpKf6rFcLLywteXl6yG96IiIjExi8BVZ+blQFM9bURl5wBAUCOVA1LO4zBzXp2+Pn4Sli9fAo1CHidmS12qPQfampqsLS0zLfcz88POTk5iImJQVhYGEJDQxEWFoawsDDUr19f1k4QBGzbtg2pqalyn1dXV4eVlRU6deqEDRs2yJY/efIEJiYmLFNfibHqYBFYdVAcnDpIRACr6xWGx0W1nbobiy933gIAvP9Lmlnyc5glP8f1+s0BABM7WGN6l0ZQV+d9P6ogOzsbW7duzZeMZWS8myb60UcfyUrVC4IAAwMDpKamwtLSUm4ErGHDhrC1tYWVlZWY3VFZrDpIREREpKR62Jti/adOWPh3iNxztHLqW2BEr+6wiXyJ332jcH/LXjyYegwmfx+AYSNL8QImhVBXV8fnn38utyw3NxfPnj1DaGgotLX/nS6akpKCzMxM5OTkyBXsyNOjRw+cPHlS9n727NkwNzeHjY0NGjVqBHNzc0il0vLtEDHRIiIiIqpsetiboqudCfwikhCfmgEjXW24WRlATSpBH4d6cKunixZdxqJecjwSnJxwb8M2NB3RT+ywScGkUinMzMxgZmYmt1xfXx9paWmIjY2VGwHL+7N58+ayti9fvsTSpUvlPq+trQ1ra2vY2Njgo48+wrhx42TrBEGQK+ZBpcepg0Xg1MHKh9MKiaoOTpErGI8LAUC07y3kDBwMq7gI5Eik8B87FW4blkHKqYT0noSEBCxduhSPHj1CaGgowsPDkZWVJVs/ceJEeHt7A3h3bjE3N0ejRo1gY2MjGwHL+7NmzZoi9aLy4NRBIiIiIhVX38MJb+4Hw3/Ap3C9+Bfct6xEsJ8vLI4dRM36pmKHR5WEoaEhVqxYIXufnZ2NqKgohIaG4tGjR2jWrJlsXWhoKFJSUhAQEICAgIB825o2bRpWrlwJAHjz5g1OnDgBGxsbNGzYEDo6OuXfGSXDRIuIiIhISenU1IXrP0fh/90KNFs6Bw53fBHX3AHBl3zh4GAtdnhUCamrq8Pa2hrW1tbo0aOH3LrmzZvj3r17stGv9/+MjY1FvXr1ZG0fPXqEwYMHy96/fw+YjY0NOnXqBAcHhwrrV2XERIuIiIhIybkumoHw9q2g+clQXKnXFPP2PcT/UqUY09qS99tQsWloaMDOzg52dnb51qWmpuL9O44yMzPh5uaGR48e4dWrV3jy5AmePHmC8+fPAwB++uknWaIVGhqKadOmoXHjxnIvY2Njlf75ZKJFREREpAKsu7RC6r0g3Dj+ANkPX+H7YyEIufMY8wY6QM/YUOzwSMnp6urKvXdzc8ONGzcgCAISExPlRr8ePXoEV1dXWds7d+7g+PHjOH5cvm6Bnp4ebGxsMHfuXPTr966YS141xWrVqpV/p8oZEy1SOoUVF2GRDCIiqup0jWpjlWcrtPCJxJJjd9HvpxlI+e4FXuzeA+tubcUOj1SQRCKBoaEhDA0N0apVqwLbODk5Yf369Xj48KHsFRkZiZSUFNy8eVOuOMfp06fRr18/1K9fHzY2NvlGwZSpND0TLSIiIiIVIpFI4NnaCi7SNBhueAaTV/HI7NUZft8uguuibyBRkl9SSXVYWlpiwoQJcssyMzMRHh6Ohw8fyiVojx8/BgBER0cjOjoa586dk/vcnj17MGzYMADAw4cP4e/vL0vCKlsVViZaRERESsLb2xve3t7IyckROxRSAvYezfDqdjCCeg+G4+1rcPtxFvyvXEHTo7tQrZa+2OFRFaelpVXg/WBTpkzBp59+Kjf6lfcKDw+HjY2NrO2JEycwffp02XsTExO50a+hQ4fKFfCoaHyOVhH4HC3lwamDRKqHz4sqGI8LlURudg5uTJwFt80roSbkItLEEpIDB2DR2lns0IhKJDs7GxKJBGpq754Vt2vXLmzcuBEPHz7E8+fP87X38/OT3Su2Z88e7Nu3D40bN8bw4cNLXRGRz9EiIiIiIgCAVF0NHhuX417HtjD+fAws4yIRMngojp64iH4tzMQOj6jY1NXlU5cRI0ZgxIgRAIDk5GQ8evRIbgTs/dEvHx8fHD16FADQsmXLCik9z0SLiIiIqApo+klfvHC7hVsDhmOOyzDc3xcM/6iX+K63HbTU1cQOj6hM9PX14erqKlft8H2enp6wsbHBw4cP0aJFiwqJiYkWERERURVRx9oCBkFX0eXcI9y/EIad16Ohe/QwPp00GPWc8j87iUhVODs7w9m5YqfLMtEilcGy70REREVTk0rwdbfGcLaohc0/78W03xchY+9yBC5fixZfjRY7PCKVwfqeRERERFVQh8ZGWDG5Jx5bNIFeRhpaTPLE9UFjkZWRKXZoRCqBiRYRERFRFWVi3wgN7t3E9X6jAADuh7YhzN4N8ffDRY6MSPkx0SIiIiKqwjR1tOB+ZAcCf9mMVK1qsA2/DQ0XZ9zZvl/s0IiUGhMtIiIiIkKLqeOQfOU6wswaodbrZBz7/SR+OfsIObkCcnIF+IYn4mjQU/iGJyInl49hJSoKi2GQyvvQA6NZKIOIiOhfZq7NkBESiKNf/4jfankA50Nx+l4cktLfIj7133u3TPW1Mb+PHXrYm4oYLVHlxhEtIiIiIpLR1q2Ofht/wC/DHKGpJkV0dDx+3vQNXGLuydrEJWfgy523cOpurIiRElVuTLSIiIiIKJ++DvWgp6OOST770DYqCHt3z8bnNw4CgoC8iYML/w7hNEKiQjDRIiIiUhLe3t6ws7ODq6ur2KFQFeAXkYSEtLf4tdVQHLFrD3UhF/+7uA2rjq2AVvZbCABikzPgF5EkdqhElRITLSIiIiXh5eWFkJAQ+Pv7ix0KVQHxqRkAgNeaOpjaewbmdJuILKka+odcws69c1HrdbJcOyKSx0SLiIiIiPIx0tX+941Egl0temH0xwuRolUdrk9DcPiPGbBKeirfjohkWHWQqrTCKhKyGiEREVV1blYGMNXXRlxyhuyeLB9LRwz4dAW2HViAalkZeKuuAUHgPVpEBeGIFhERERHloyaVYH4fOwCA5L3l4YbmGDDyZ4wa8j2e6hlh9DY/HAyIESdIokqMiRYRERERFaiHvSnWf+oEE3356YGadU0wcfIAfNTMFFk5As78sB7XR3hByM0VKVKiyodTB4mISKWcP38eFy5cgI+PD2JiYpCQkIBq1aqhTp06aNasGdq3b4/evXvDxMRE7FCJlEIPe1N0tTOBX0QS4lMzYKSrDTcrA6hJJejdvC6aIg1jfv4ZOtmZCIh8jKanDkBbt7rYYROJjokWEREpvbS0NKxZswabNm1CdHS07J4RbW1tGBgY4M2bN7h79y5u376NXbt2QV1dHX379sW0adPQunVrkaMnqvzUpBJ4WNfOt1wqlWDiiHbwf7AYjj/OgrPPKdxv0Qom50+ilkVdESIlqjyYaBEVoLAiGQALZRBVNhs2bMCCBQsQHx8PBwcHfP755/Dw8ICLiwtq1KghaycIAkJDQ3Hjxg2cOXMGR48exeHDh9GvXz/8/PPPsLKyErEXRMrNddEM3G3UAPU/HwXb8NuIcXZF6t9/o76Hk9ihEYmmTIkWp2cQEZHYJk2ahBEjRuCbb75B06ZNC20nkUhgY2MDGxsbjBw5Em/evMGePXuwZMkS/PHHH5g3b14FRk2keuxHDUSkVX2k9esDs8RnSO7UHiFbdsJueD+xQyMShUQoYU3O4kzPSE5ORu7/3wyp7NMzUlJSoK+vj+TkZOjp6ZVqGx8aHSHlwxEtoopTnHNweHg4rK2tS72PnJwcxMTEwMLCotTbqGiKuDYRlZfEiCdI6NwTjSPuYXPLAajlvQaDnM3EDotIIUpy/i1R1cENGzagYcOGmDt3LmrWrInFixfjwoULSElJwevXrxETE4PExERkZWXhwYMH2LFjB4YOHYozZ86gXbt2GDhwICIiIsrUOSIioveVJckCADU1NaVKsogqu9pW5rAIuo6Dw6fhx3ae+Hp/MFaeecjnbVGVU6JEa9KkSejRowfu3LmDwMBAzJ49Gx06dJCbAw/8Oz1j5MiR+OOPP/D8+XNs2rQJd+7cwR9//KHQDhARERFR5aKtVwMD/vgZEzrZAADWn72PY8OnICM1XeTIiCpOie7RevDgQam+OdTR0cHYsWMxevRoxMTwgXZEREREqk4qlWBmjyawrF0dwvjx6BN8Gvf9L7EiIVUZJUq0OD2DqPB77njvFhERUX5DXM1xd/rnSPn86r8VCY8dQ333FmKHRlSuSjR1kIiIiIiopOxHDUTS2Yt4VssEZonPoN+xHe7t+UvssIjKVbkkWpGRkZg4cSKGDBmC2bNnY+/evQgJCZFVIiQiIhJLcnIybt26JXYYRFWOZVsXaAX44aFVU+hnpKHRp4Pgv3CV2GERlZtySbQGDx4MHx8f1K9fH6GhoZg7dy6aNWuGGjVqwNnZuTx2SUREVCzZ2dkYNmwYtm7dit9++03scIiqlLyKhLdadoFmbjYaL5mLdfuvsyIhqaQyPbC4MPfv38fNmzdha2srW5aamoqgoCDcvn27PHZJRERUpGPHjmHFihUICwvD119/jcePH4sdUol4e3vD29sbOTk5YodCVGraejXgePUUfD/1wlqY41pAIh7mBOGnQc2hraEmdnhEClMuI1pOTk5ITk6WW6arq4u2bdvCy8urPHZJRERUpN69e2Pu3LmQSCSoVasWzpw5I3ZIJeLl5YWQkBD4+/uLHQpRmUjV1eCxdwP6fT0K6lIJjgY9w4K52/Ay6pnYoREpjERQ0Fjt2LFj0bx5czg4OOD169fw9vbG/v37Ub16dUVsXjQlefpzYQqrUkdVA6sREpWeIs7B/5WYmIhly5bhp59+QlhYGBo2bKiQ7Vak8jguRGK5FpaApauO4I8t05BWXQ+5rEhIlVhJzr8KG9HS0tLC/v370a9fP/Tt2xdnz55FkyZNMGPGDOzfvx+PHj3i/FsiIhKdvr4+fvjhBwBQyiSLSNW0bmiItcMc8VqnBisSkkpR2D1a69evl/09PDwcwcHBstfBgwcRFRWFatWqoWnTprhx44aidktERFQi6urlcnsyEZWBRRsXJAb44WHnnmgccQ86nw6Cf+gKuM6bInZoRKVWLlcba2trWFtbY+DAgbJlKSkpLIZBRERERAWqbWWO6kHXcavbADjdOAfX+VPh++gR3H//FRIpH/1KyqfEP7VhYWGl2pGenh7atWuHr776qlSfJyIiIiLVJqtIOPRzAIDHrnXYP/47ZGSx0iYpnxKPaNnY2EBPTw+Ojo5wdnaGk5MTnJyc0KRJE0gkkvKIURQsoUtERERU8d5VJPwN/o0aIeuPnZhbyxV/br6BjaNcYFBdU+zwiIqtxFUHpVIpJBKJrLBFXnJVrVo1ODg4yBIvZ2dn2NnZQU1NuZ+HwKqDVJ5YkZDow4p7DlbW6oGlxaqDVFVcexSPCbsDkZqRDcta2tjevR4sHZsgJ1eAX0QS4lMzYKSrDTcrA6hJVecLf6q8SnL+LfGIloGBAZKTk9GrVy8MGDAAkZGRCAgIwK1bt+Dj4wMfHx9Z8qWlpYVmzZrB2dkZ69atK11viIiIilBVZlsQVTWtbYxweGIreG7zx7BD3qi14AwO/rABKzJNEJucIWtnqq+N+X3s0MPeVMRoieSVeETr1atX+O6777Bhwwbo6enh+++/x5dffgmpVIrnz5/Lkq68P588eQKJRKK0U/A4okXliSNaRB9W3HMwZ1sQqbaEhGQkurVG44h7eCtVx6yek3DIvrNsfd7XKes/dWKyReWqJOffUj+w+O7du5g6dSouXLiApk2bYtWqVejcuXO+domJiQgICEC3bt1KsxvRMdGi8sREi+jDinsONjQ0LHS2RWxsLACo1GwLJlpUFaW/SsXlVr3Q8/5VAMDqVp/glzbDgf//vy0BYKKvjavfduI0Qio3FfLAYnt7e5w7dw4HDhxAeno6unXrhoEDByIiIkKuXe3atZU2ySIiIuUQFhaGCRMm4MSJE/j6669Rp04dHD16FE+fPkVsbCyOHTuGhQsXom/fvqhTpw78/f3x22+/iR02EZXA7cS3mNhnJrzdPwYATPHZg1+O/Qz1nGwAgAAgNjkDfhFJIkZJ9K8yP5Rg4MCBePDgAb7//nucPXsWdnZ2mDNnDtLT0xURHxERUZFq1qyJX3/9FYGBgWjRogUmTZoEBwcHnD9/HsbGxujVqxfmzp2Lw4cPIyoqCi9evMDJkyfFDpuISiA+NQOCRIrl7Ufjm56TkSVVw4CQi1jz1zLgvQla8akZH9gKUcVRyNPfNDU1MWfOHDx8+BCDBg3CkiVL0LhxY9y7d08RmydSWZazjhf4IqLS4WwLItVlpKst+/v+5t0wfuBcvNbQwsnGrWXTB//bjkhMCkm0YmJicPLkSezcuRMSiQS1atVCbGwswsPDFbF5IiKiEuFsCyLV42ZlAFN9bVnhi4vWrmj7xRb8bdde1sZE712pd6LKoMTl3a9du4Y7d+7IvVJSUmSVnurUqQNnZ2c4OjrCwcFB4QETEREVR95sizFjxmDmzJlYsmQJduzYgdOnT6Np06Zih0dEJaQmlWB+Hzt8ufMWJHh3T1Zi9Zqy9capCfjx1A6kjW4G/XpGYoVJJFPiRKtt27aQSCSQSqVo0KABunXrBkdHR1liVbdu3fKIk4iIqERiYmJkXwj+d7YFEy0i5dTD3hTrP3XCwr9D5J6jZVhdA+t2L4PzkxCEt2yDnCsXYGBlJmKkRKVItABAXV0dPXv2RLt27WTPJmF5WSIiEgtnWxBVHT3sTdHVzgR+EUmIT82Ake676YLRbTYjoe9HsH4aikiPtsi9fAGGNlZih0tVWIkTLXt7ezx48AB//fUX/v77b9nyBg0ayB4GmZd8GRhwjiwREZU/zrYgqlrUpBJ4WNeWW2bV0QPRZ84hvkd3WD6PREyrtoj75wJMmtmIFCVVdaV6YHFmZiaCg4Nx69Yt2QMh7969i6ysrHcb/f/KL/Xr15clX//73/8UG3kF4QOLqbLhQ46pKinuOVgqlUJDQ6PKzLbgA4uJCvcsMATo3Bl1X8YhrqYRsk+fhZlbc7HDIhVRkvNvqRKtgmRlZeH27dtyydedO3eQmZkJiUSCnJwcReymwjHRosqGiRZVJcU9Bzdv3hwPHjxAdna27Ms+QHVnWzDRIvqw5yGheNuhM8xfPIFfgxao5XMJjYx1xQ6LVEBJzr+lukerIBoaGnB2doazszPGjx8PAMjOzsa9e/cQEBCgqN0QERHlc/v27UJnW4SHh2P//v0qNduCiD7M2K4REnyv4tqAEZjWdjyyN17HH+Pc0LSuvtihURWisBEtVcURLapsOKJFVUlZz8GcbUFUtb1Mf4vR2/xwOyYZetrq2DmgIZo7WIsdFikxUUa0iIiIKhvOtiCq2mpV18TOz1pi7DZ/mJ04BMulg3Bv+x40HdZH7NCoCmCiRaRkChsh5UgXUfGoq6vDwcGBZd6Jqgg9bQ38PtYVod7ToZeZjgYjP8bt17+j+dghYodGKk4qdgBERERVzZMnT9ChQwfY2dmhefPm2L9/v9ghEam0aloaaOxzDsHNW0MnOxNNPh+BwF+3ix0WqbgyjWiNHTu2yDZSqRR6enpo3LgxevfujXr16pVll0REREpPXV0dq1atgqOjI+Lj4+Hk5IRevXqhevXqYodGpLK09WrA1vccbnX4CE7+F9BsyjgEvHkD55lfih0aqagyJVrbt2+XVXEqqKaGRCKRWz5p0iTMmzcPc+fOLctuiYiIlJqpqSlMTU0BAEZGRjAwMEBSUhITLaJypllNG82vnMTNLgPhcvU4HGd9Bb/0N3BbOF3s0EgFlSnRCg8Px9SpU+Hv748pU6agVatWMDY2xvPnz3Ht2jWsWbMGbm5umDNnDoKDg7F48WLMnz8fjRo1wtChQxXVByIiIpmKmG1x+fJlLF++HAEBAYiNjcXhw4fRv39/uTbr1q3D8uXLERsbi6ZNm2LVqlVo27Ztvm3dvHkTubm5MDc3L1EMRFQ66lqacPrnKPx6DoXbuYO4deoaHvYYiJEelmKHRiqmTInWvn374Ofnh+DgYBgZGcmW29jYoG3btvD09ISjoyP++ecfzJw5Ez179oSdnR3WrVvHRIuIiMpFRcy2SE9Ph4ODA8aMGYNBgwblW79v3z5MnToV69atQ+vWrfHbb7+hZ8+eCAkJQf369WXtEhMTMWrUKGzevPmD+8vMzERmZqbsfUpKSrFjJaL8pOpqcD39J/6cswZLcxsCR+/hTVYOPm/H0u+kOGV6jlajRo3Qs2dPrFmzptA2kyZNwqlTpxAaGgoAGDFiBI4fP45Xr16VdrcVis/RImXHaoSkzEpzDo6IiCjxbIsnT55g9+7dpfoSUCKR5BvRatmyJZycnLB+/XrZMltbW/Tv3x9LliwB8C556tq1K8aPH4+RI0d+cB8LFizAwoUL8y3nc7SIykYQBPx85hHW/hMGrey3WKMWhm5LvoZEynpxVLAKe45WTEwMtLS0PthGW1sbMTExsvf169dHRkZGWXZLRERUKLFnW7x9+xYBAQGYNWuW3PJu3brBx8cHwLtf7jw9PdGpU6cikywAmD17NqZP//cekpSUFE41JFIAiUSCGd0bQ0ddgsZffIou4f7wffwQ7vs2MtmiMivTT1C9evVw9OhRuekM78vMzMTRo0fl5r7Hx8ejVq1aZdktERFRobZs2YKPP/5YLsl6n4mJCT7++GNs2rQJwLtrWe/evREcHKyQ/SckJCAnJwfGxsZyy42NjREXFwcAuHbtGvbt24cjR47A0dERjo6OuHPnTqHb1NLSgp6entyLiBTHq7MNavTqDgDwOLAFfn0+RW52jshRkbIrU6I1btw4hIWFoX379jh+/DiSkpIAAElJSTh27BjatWuH8PBwuRuTr1y5wodEEhFRuakssy3y7hPLIwiCbFmbNm2Qm5uLoKAg2atZs2YK3T8RlYz7mkW4MWsJciFByxN7cLPbYORkZYsdFimxMk0dnDlzJu7fv4+dO3eib9++AN5VcsrNzQXw7qIyYsQI2fSJ58+f46OPPkKPHj3KGDYREVHB8mZbLF68uMCEq7xnWxgaGkJNTU02evX+Pv47ykVElUvLJbNws5oOWsyfDrd/juBmhz5wOH8EGtof/vKGqCBlGtFSU1PD77//jrNnz2LUqFFwdHSEpaUlHB0dMXr0aJw9exZ//PEHpP8/x9XY2Bi//PILunfvrpDgiYiI/kvs2RaamppwdnbG2bNn5ZafPXsWrVq1KtO2vb29YWdnB1dX1zJth4gK5/LdFAQvX48sqRpcfE7Bv+sgZHIaIZVCmUa08nTu3BmdO3dWxKaISME+VPWSFQlJFVXEbIu0tDSEhYXJ3kdERCAoKAgGBgaoX78+pk+fjpEjR8LFxQUeHh7YuHEjoqOjMWHChDL1zcvLC15eXrKqV0RUPpy+/hzB1XRgPuMrLLfsCL3fA7DhU2foaKqJHRopkTKVd68KWN6dVBkTLarsynIOPn/+PHbu3Inbt28jJSUFenp6cHBwwIgRI8r85eDFixfRsWPHfMtHjx6N7du3A3j3wOJly5YhNjYW9vb2+OWXX9CuXbsy7TePIq5NRFQ03+BIjD3wAG+ycuDewACbR7uihpZCxilISZXk/KuQRMvHxwfbt29HUFCQbKd50wfbtGlT1s2LiokWqTImWlTZMaEoGI8LUcXxj0zCmG3+sIh6gO/996Lhub+gb1pH7LBIJBX2HC0AmDFjBn755Rfk5Wt50zMCAgKwdetWTJkyBStXrizrboiIiIiIKpyrpQF2eTqjpss4WCQ9RZhba+RevYhaFnXFDo0quTIVw/j999+xcuVKNG7cGHv27EFsbCyys7MRFxeHvXv3okmTJli9ejV+//13RcVLRERULD4+Pvj888/h5uaGxo0bw9XVFePHj8fVq1fFDo2IlIyDlSFy/vwTSdX10TAmFMnurZHwKELssKiSK9PUQQ8PDzx79gx3796Frq5uvvUpKSlo1qwZTE1Ncf369TIFKhZOHaSqiFMKqbIo7Tm4sNkWwLvnWynrbAtvb294e3sjJycHjx494tRBogoWdS0AOj27wyg1EU8MzaDxzwWY2DdCTq4Av4gkxKdmwEhXG25WBlCTSoreICmdklyXyjSidffuXQwaNKjAJAsA9PT0MHDgQNy7d68suyEiIio2VZ5t4eXlhZCQEPj7+4sdClGVZNHaGVn/XMKzWiYwT4iB0LYtDu+/hDY/XcAnm65jyt4gfLLpOtr8dAGn7saKHS6JrEyJFgAUNSAmkTCbJyKiirN+/XqYm5vjxo0bGDp0qOwhwUZGRhgyZAh8fX1hZmaGdevWiRwpESmjes5NIb18GU/qmMP01XPkLlyI2OQMuTZxyRn4cuctJltVXJkSLXt7exw8eBBpaWkFrk9NTcXBgwfRtGnTsuyGiIio2DjbgojKm4l9I2hevYS9jt0xp7tXvvV5wxAL/w5BTi6fpFRVlSnRmjBhAmJiYuDh4YGDBw8iISEBAJCQkIADBw6gVatWiImJwZdffqmQYImIiIqDsy2IqLw9VtPDrO6TkKGhLVtWJ+2l7O8CgNjkDPhFJIkQHVUGZSrvPnr0aAQFBWH16tUYMmQIAPkbjgVBwKRJkzB69OiyR0pEFeZDBVxYKIMqu7zZFosWLUKNGjXyredsCyJShPjU96YLCgJmXPkDIwJPYsjwpQitY1FwO6pSynyP1i+//ILLly/D09MTjo6OsLS0hKOjI8aMGYNLly5h9erVioiTiIioWFR5toW3tzfs7Ozg6uoqdihEVZ6R7r8jWVrZb9E6Mhi1MlKx88/vYP4qrsB2VLWUqbx7VcDy7kTyOKJFFam05+Bp06Zh9erVsimCBc22UOYvAhVxbSKissnJFdDmpwuIS86AAED/TSr27pkN2xeRiNY3xscjfoLUzAxXv+3EUu8qpMLKuxMREVVGnG1BROVNTSrB/D52AAAJgGQdXYwasgiRNU1RP/k5du77DrNcDZlkVWFlukeLiIiosmrTpg3atGkjdhhEpMJ62Jti/adOWPh3CGKTM/CiRi18Omwx/tz1LRolPgHGDkXKzWvQM6otdqgkghIlWmPHji3VTiQSCbZs2VKqzxIRERERVVY97E3R1c4EfhFJiE/NgJGuNjKHOyCpZ2c0evIQP8/dgIlrZ0FHU03sUKmClSjR2r59e6l2wkSLiIiIiFSVmlQCD+v3Rq2sayPs4F9YvuFv7DFwxJ1dAdg40gWa6rxrpyopUaIVERFRXnEQkZIorLgLi2SQWDjbgogqo4bd22FgY3sc3nIDFx++wOztV7FsdCuoafDOnaqiRP/SFhYWRTciIiKqQFVptoW3tze8vb2Rk5MjdihEVAyulgb4baQLZnmfxmf/+wq39raAy5kDkEg5slUVMKUmIiKlVpVmW3h5ecHLy0tWXpiIKr/2NnXwa6Nc2CREQ+18JK4PGouWB7cy2aoCSpRo9e7dGwsXLoSzs3OJd/TmzRt4e3ujevXqSvmQSCIiqpw424KIKjuXr0bBL/4F3BbNgPuRHfAdVxMe21aJHRaVsxKl0k+ePIGbmxs6d+6M7du3IyUlpcjP3Lx5E1OnToWFhQXmzZsHQ0PDUgdLRERERKSM3L7/GtcnfwcA8Ni+GjemzBM5IipvEkEQhOI2FgQB27Ztw/fff4/o6GhIpVI0adIETk5OMDY2Rq1atfDmzRskJSUhNDQUN2/eRHJyMqRSKYYMGYIffvgBlpaW5dgdxSvJ058LU1jxAKKqgEUyqCyKcw6uirMtFHFtIiJx+I6eAo/f1wAA/Ob9DLeF00WOiEqiJOffEk0dlEgkGDt2LDw9PXH8+HFs374dly5dws6dO/O1lUqlaN68Ofr374/PPvsMdevWLVkviIiIiiFvtkWHDh0wcuRIDBw4sMiL382bN7Fz507s3r0baWlp2LFjRwVFS0RVnfu2X3A9+RXcj/4O09XLcKZPf3RzaSB2WFQOSjSiVZj79+8jJiYGiYmJ0NHRQZ06ddC0aVOVuFGXI1pEZcMRLSqL4pyDOduCI1pEykbIzcX5oRMxr05LvKhljC2jXdHOpo7YYVExlOT8q5BES5Ux0SIqGyZaVBYlOQfn5ubKzbZISkrK10ZVZlsw0SJSfjm5AibvCcTxO7HQ1pBi9yf2cLIzFzssKkK5TR0sTFpaGkJDQyEIAho1agRdXV1FbJaIiKjYpFIp+vTpgz59+gBQzdkWfI4WkepQk0rwy1BHpGVmQ+fvI7BwG46wg0fRsHs7sUMjBSnTiFZ8fDymTZuGgwcPIisrCwCgrq6OQYMGYeXKlTAxMVFYoGLhiBZR2XBEi8qCIzcF43EhUh1vMrMR3rwl7B/dQlJ1faSdvYD6Hk5ih0WFKMn5t9RPSouPj0fr1q2xZ88eSCQSODs7w8nJCWpqati7dy/atm2L+Pj40m6eiIioTNLS0hAYGIhbt24hNTVV7HCIiAqko6WO+pdPI8zMBgbpydDq2QOxtx+IHRYpQKmnDk6ePBnh4eEYN24cfv75Z1lGl5KSgunTp2Pr1q2YNGkS9u3bp7BgiUj5fGhEl6NdVB6qwmwLIlItesaGyL58HlHurWERH42YTl2QcMMHhtb1xQ6NyqBUI1pPnz7F/v370b59e2zatElu2ExPTw+bN29G+/btcfDgQURHRyssWCIiog/hbAsiUlYGVmbQuXgBz2qZwCzxKVLadULyU56vlFmpEq0LFy4AAL755ptC28yaNQu5ubk4f/586SIjIiIqofdnWzx//hx+fn7w9/dHXFwcxo4di/DwcEyaNEnsMImICmRka43c02eQUMMADZ6F4+8vv0N6ZrbYYVEplXjq4NixY3H79m0AwO+//44DBw4U2O7NmzcAgF9//RVXr17Fli1byhAmERHRh/13tsX78mZbhIeHy2Zb1K/PKTlEVPmYuTZDxN8ncPS7X7DYtjdO/nETW0a7QltDTezQqIRKnGht374dACCRSPDnn39+sK1EIkFQUBCCg4OZaBERUbkq7myLS5cu4fz58xgzZkxFhUZEVCJWHVrCefcGVNt0HdfCEjF5102s+6QF1LU0xQ6NSqDEiZavry927doFb29v7NmzB5aWlgW2i4qKwrBhwzBp0iQMHz68rHESkQoqrFAGi2RQSXG2BRGpGkfzmtg02gXjN/ug79IZCNqiD6d/jkKqzpEtZVHiRKtly5ZISkrC2rVr8fz5cwwdOrTAdv7+/pBIJOjZsydatmxZ5kCJiIgKw9kWRKSKWlkbYru9FI7LfKHxMAc3+o6A27HdkEhL/YQmqkCl+lfq2LEjDAwMsHjxYsTGxuZb/+zZMyxcuBAGBgbo2LFjmYMkIiL6EF9fX3z11VcAgD179sDX17fA1969eyEIAiZNmgQfHx+Roy45b29v2NnZwdXVVexQiKiCuI7ojeDFq5ELCVqe3Icbw78UOyQqplIlWtra2vjuu++QkJAAd3d3bN68Gffv38f9+/exadMmuLu7IykpCfPmzYOWlpaiYyYiIpLTsmVL9OzZE4Ig4Pnz52jZsmWBr/j4eKWebeHl5YWQkBD4+/uLHQoRVSCX2V7w/3YxAMB930Zc/2KmyBFRcUgEQRBK++GvvvoK69atg0QikVsuCAImTpyItWvXljlAsaWkpEBfXx/JyclyzwsriQ89sJWI8uM9WpSnJOfgjIwMmJmZQSqVIjg4GKampnLrnz17BgcHBwBATEyMUn8RqIhrExEpn+sTvoX7b8sAAH4zF8PtpzkiR1T1lOT8W+J7tN63du1a9OnTB1u3bkVISAgEQUDTpk0xduxYdO/evSybJqIq7ENfTjAJo8LkzbaYNm0a3N3d8d1336F169YAgKtXr2LRokVISkrCqlWrlDrJIqKqy33DT/BNTobH3t9g/8sinOraCz26tBA7LCpEmRItAOjevTuTKiIiqhSmTJmC0NBQrFu3Dl988YXcurzZFnxgMREpM/dd63Al8y3Watsg4EIsfqtXF51tjcUOiwrAkiVERKRS1q5di5MnT2Lw4MGws7ODra0tBg8ejJMnT6rElHYiqtokUilaH9iCuv16IDtXwJe7bsHnYZzYYVEByjyiRUREVNlwtgURqTKpVILlg5sjLTMb0ZduoG6r8Xi0fQds+nQWOzR6D0e0iIiIiIiUjLqaFL9+0gILgw/BMukpjIf0R8TFG8jJFeAbnoijQU/hG56InNxS172jMuKIFhERERGREtLWUEOz0wfx0LUNGkfcw9s+PTFo/CoEadaWtTHV18b8PnboYW/6gS1ReeCIFhERERGRkqpeuyZMrlxAmIkV6qS9xK9bZ8IkJUG2Pi45A1/uvIVTd2NFjLJq4ogWESmVwkq/s+w7ERFVVTVM62DY6KVYt3EqrF7GYue+uRgy4ickVdOHAEACYOHfIehqZwI1qaSozZGCcESLiIiIiEiJ+UUk4T6q49OhP+Cpbh00TIrBzEs7ZOsFALHJGfCLSBIvyCqIiRYRERERkRKLT80AADzVN8LIoYtwrHEbLOr0WaHtqGIw0SIiIlIS3t7esLOzg6urq9ihEFElYqSrLfv749pm+Kr/LKRrVftgOyp/TLSIiIiUhJeXF0JCQuDv7y92KERUibhZGcBUXxv57r4SBHzmdwgjbx2Dqb423KwMxAivyqoSxTCePHmCkSNHIj4+Hurq6vjuu+/w8ccfix0WESkQi2QQEVFVpSaVYH4fO3y58xYkeHdPFgC0j7iFuf9sRY5ECoe2jlCT8oHGFalKjGipq6tj1apVCAkJwblz5zBt2jSkp6eLHRYRERERkUL0sDfF+k+dYKL/7/TAS1ZOOOjQFWpCLnounobH531EjLDqqRIjWqampjA1ffeQNiMjIxgYGCApKQnVq1cXOTIiIiIiIsXoYW+KrnYm8ItIQnxqBox0teH4XWfcc26Dpo9uodqg/ki46QfDhpZih1olKMWI1uXLl9GnTx/UrVsXEokER44cyddm3bp1sLKygra2NpydnXHlypUCt3Xz5k3k5ubC3Ny8nKMmIiIiIqpYalIJPKxro59jPXhY14ZODR2YnT+OJ3XMYZL8Akmde+HNq1Sxw6wSlCLRSk9Ph4ODA9auXVvg+n379mHq1KmYM2cOAgMD0bZtW/Ts2RPR0dFy7RITEzFq1Chs3LixIsImIiIiIhKdvpkJJMeO4WU1PdhE38f9bv2Rm5MrdlgqTykSrZ49e2Lx4sUYOHBggetXrlyJcePG4bPPPoOtrS1WrVoFc3NzrF+/XtYmMzMTAwYMwOzZs9GqVatC95WZmYmUlBS5FxERERGRMjNza47YrbvxRl0LR/QbYdmZR2KHpPKU/h6tt2/fIiAgALNmzZJb3q1bN/j4vLvhTxAEeHp6olOnThg5cuQHt7dkyRIsXLiw3OIloopVWDVCgBUJiYioarEb+hFO1PbF7+eeAZfCYWVYDUNd64sdlspSihGtD0lISEBOTg6MjY3llhsbGyMuLg4AcO3aNezbtw9HjhyBo6MjHB0dcefOnQK3N3v2bCQnJ8teT548Kfc+EBERERFVhF5dWmBy50YAgJ92++L2wTMiR6S6lH5EK49EIv+INkEQZMvatGmD3NzizUPV0tKClpaWwuMjIiIiIqoMpnVphFcPH2P0/Bmos/4loowvwqKNi9hhqRylH9EyNDSEmpqabPQqT3x8fL5RLiIiIiKiqk4ikeB/o9siq2Yt6GWmQ71fX7yMfCp2WCpH6RMtTU1NODs74+zZs3LLz549+8GiF0REREREVZW2bnUYnTuBZwamqJcUi7guvZCRmi52WCpFKRKttLQ0BAUFISgoCAAQERGBoKAgWfn26dOnY/Pmzdi6dSvu37+PadOmITo6GhMmTBAxaiIiIiKiysvAygxZR/9CilZ12Ibfxt2eH0Mo5u02VDSluEfr5s2b6Nixo+z99OnTAQCjR4/G9u3bMXToUCQmJuL7779HbGws7O3tceLECVhYWIgVMhEpgcIqErIaIRERVRUWbVxwd+PvaDxmCFyuncR1zylw//1XscNSCUoxotWhQwcIgpDvtX37dlmbiRMnIjIyEpmZmQgICEC7du3EC5iIiIiISEnYjxqIwFk/AADq/bUfx64+FDki1aAUI1pERERERFR+3H74FqeSX2OupBFSTj6GiVkduFgaiB2WUlOKES0iIiIiIipf3dYsgIubLd7m5OLzPwIQ9TxZ7JCUGhMtIiIiIiKCVCrBL0Md0dxMHx2vn0R2CyckP40XOyylxamDRET/UViRDICFMkhc3t7e8Pb2Rk5OjtihEJGK0tFUw+ZBtsiduxMmyS9wr3Mv6Ny6Cs1q2mKHpnQ4okVERKQkvLy8EBISAn9/f7FDISIVZmRaG+kHjiBdUwdNHwYgqM9wln0vBSZahfD29oadnR1cXV3FDoWIiIiIqEJZd2mFsDWbkSORwu3CYVz/crbYISkdJlqF4LeGRERERFSVOXwxHDenzgMAeGxchlsrN4sckXJhokVERERERAVquXI+bnw0HABg960X7l+6KXJEyoPFMIiIiIiIqFDOB7chyC0KF3XrY+fllzjq8Ab1auqIHValx0SLiKgECqtIyGqERESkqtS1NNHQ9zxmbbyBhLhUjNvuj/0TPKCrrSF2aJUapw4SEREREdEH1aimhS2erqijq4XIJy9w4tNpyM58K3ZYlRpHtIiIiIiIqEj1aupgyyhnvGnbAS2jbuNGegLcju+GRMqxm4LwqBARERERUbE0N68FjWlTkAsJWp7ahxvTF4odUqXFRIuIiIiIiIrNacpY+E2YCQBwXbMYQev/EDmiyomJFhERERERlUhL7x/h13kg1IRc2Ez5HOFnr4kdUqXDe7SIiBSgsGqEACsSEhGR6pFIpWhxbDfuOrWB/f2bqPHxALy4cR11GjcQO7RKgyNaRERERERUYhraWjA/fwJRRvWh/jYTP2w+j9dvs8UOq9JgokVERERERKWib1oH6idOwPOLX3FErS6m7A1CTq4gdliVAhMtIiIiIiIqtXrOTbFgam9oqktxNuQ51u6+InZIlQITrUJ4e3vDzs4Orq6uYodCRERERFSpOVsYYPng5ugU5ofPxnbHjdlLxQ5JdEy0CuHl5YWQkBD4+/uLHQoRERERUaXXz7Eexld/hepZGXD+aQ7ubD+AnFwBvuGJOBr0FL7hiVVqWiGrDhIRlbPCKhKyGiEREaka922/4ObjcLhcPQ6LCZ4YdusV/KuZytab6mtjfh879LA3/cBWVANHtIiIiIiISCEkUimanfoTd6yaQS8zHSt/n4va6a9k6+OSM/Dlzls4dTdWvCArCBMtIiIiIiJSGHUdHUwfNh+RNU1hnvwcmw4tglZWJgAgb+Lgwr9DVH4aIRMtIiIiIiJSGL+IJITmamPs4PlI1qoOp2cP8Zn/Edl6AUBscgb8IpJEi7Ei8B4tIiIiIiJSmPjUDADA49pmmDBgDjqF++G3loMKbaeqmGgREYmERTKIiEgVGelqy/7ua9EcvhbNi2ynijh1kIiIiIiIFMbNygCm+tqQ/Ge5ek42pl3ZBcP0lzDV14ablYEo8VUUJlpEREQiGDBgAGrVqoXBgweLHQoRkUKpSSWY38cOAOSSrR9Oe2OKzx78+tcyfNe9EdSk/03FVAsTLSIiIhFMnjwZv//+u9hhEBGVix72plj/qRNM9P+dHrix5UCkaerAI/oOav34vYjRVQzeo0VERCSCjh074uLFi2KHQURUbnrYm6KrnQn8IpIQn5oBI113PGwshfM3X8Dj4BYErvJAi6njxA6z3HBEi4iIqIQuX76MPn36oG7dupBIJDhy5Ei+NuvWrYOVlRW0tbXh7OyMK1euVHygREQiU5NK4GFdG/0c68HDujacZ3yO6/1HAwAafTsZ0b63RI6w/DDRIiIiKqH09HQ4ODhg7dq1Ba7ft28fpk6dijlz5iAwMBBt27ZFz549ER0dXcGREhFVPs57fkNIQwfUePsaOQMHIT0pWeyQygWnDhIRVTKFlX0HWPq9sujZsyd69uxZ6PqVK1di3Lhx+OyzzwAAq1atwunTp7F+/XosWbKkxPvLzMxEZmam7H1KSkrJgyYiqiQ0tLVgdOIwEpxcUCfpOTb8ehjT542GRKJaxTE4okVERKRAb9++RUBAALp16ya3vFu3bvDx8SnVNpcsWQJ9fX3Zy9zcXBGhEhGJxrCRFeJ/34OBnr/g1zd1sN0nUuyQFI6JFhERkQIlJCQgJycHxsbGcsuNjY0RFxcne9+9e3d8/PHHOHHiBMzMzODv71/oNmfPno3k5GTZ68mTJ+UWPxFRRbEb0A1DR777UuqH4/dxMyJR5IgUi1MHC+Ht7Q1vb2/k5OSIHQoRESmh/06BEQRBbtnp06eLvS0tLS1oaWkpLDYiospibGtLBEa/ROLfp1Gj9RS8uHIOdawtxA5LITiiVQgvLy+EhIR88BtGIiKi/zI0NISamprc6BUAxMfH5xvlIiKq6iQSCX7q3xSLL21Bk9gwxPcaiOzMt2KHpRAc0SIiUiKFFcpgkYzKQ1NTE87Ozjh79iwGDBggW3727Fn069dPxMiIiCqn6tW0oL5/H9I7tUXTR7fgO+JLeBzYInZYZcYRLSIiohJKS0tDUFAQgoKCAAAREREICgqSlW+fPn06Nm/ejK1bt+L+/fuYNm0aoqOjMWHChDLt19vbG3Z2dnB1dS1rF4iIKhWL1s54+ONqAIDHwa24tXKzyBGVnUQQBEHsICqzlJQU6OvrIzk5GXp6eqXaxodKNRMRKYKqjmgp4hxcHi5evIiOHTvmWz569Ghs374dwLsHFi9btgyxsbGwt7fHL7/8gnbt2ilk/5X1uBARldX1gWPgfng70jSrIeniFdT3cBI7JDklOf9y6iAREVEJdejQAUV9Tzlx4kRMnDixgiIiIlINLnt+w71mwWgaGoiEAYOQfi8Q1WvXFDusUuHUQSIiIiIiqhTUtTRhfOIQXugaIKi2Bb47eq/IL7YqK45oERGpgA9NUVbVaYVERKSaDBtaIvjcFcw4HIHsR69gfy0SY9tYiR1WiXFEi4iISEmwGAYRVRUObnb430d2AIAlx+/h9rXbIkdUcky0iIiIlASf8UhEVcmY1pYY3FAX3gd/gGnPjkgIixQ7pBJhokVERERERJWORCLBwgHNYZMShzqpSXjeayCylOhhxky0iIiIiIioUqpeuyakhw4gTbMamoYG4ubwsj2PsCIx0SIiIiIiokqrfitnPFry/w8zPrQNt37eJHJExcOqg0REKq6wioSsRkhERMrCafpnuH7NB+6HtqHx/6Ygyt0JFq2dxQ7rgziiRUREpCRYdZCIqjKX3Rtwz8YJ1d++QdaQYUh/U7nv12KiRUREpCRYdZCIqrJ3DzM+jDvmtpjR+UvMPHy3Uj/MmIkWEREREREpBUPr+nh7+QrumjXB8dux2HotUuyQCsVEi4iIiIiIlIazZW3M+cgWAHBk2zGEHDotckQFY6JFRERERERKxbOVJaZrxmL/7zNQx3N4pXyYMasOFsLb2xve3t7IyckROxQionLBaoRERKSsJBIJxk0ZjLiNP8DyeSRCeg2A/u3r0NDWEjs0GY5oFYI3HBMRERERVV7Va9eE9PBBpGlWg11oEAI++ULskOQw0SIiIlISLO9ORCSvvocTQn9aAwBwP7IDASs2ihzRv5hoERERKQnOtiAiyq/F1HG4PmgsAKDJnKmIunpT5IjeYaJFRERERERKzWXXetnDjIO+/QFpmdlih8REi4iIiIiIlFvew4zXdh+Haa3H4NsDt0V/mDETLSIiIiIiUnqG1vXh8dsySNXVcfxOLLZcjRA1HpZ3JyIiOYWVfQdY+p2IiCo3ZwsDzP3IFksOB6HalK8Q8q0X7IaKc+3iiBYREREREamM0a0ssSriJIYHnoTRZ6OQ8EickS0mWkREREREpDIkEgnab16OCBMrGKYlIf6jAcjKyKzwOJhoERERERGRSqlWSx/qhw8iVasa7MKCcVOEhxnzHi0iIiq2wu7f4r1bFcPb2xve3t7IyckROxQiokrP3L0Fbv30K5ymjoPHkR24ucwdWYM+RnxqBox0teFmZQA1qaTc9s9Ei4iISEl4eXnBy8sLKSkp0NfXFzscIqJKz2nKWPhe9YHHgS2wnTsNfR8B4YbmAABTfW3M72OHHvam5bJvTh0kIiIiIiKV9XL2PPjUb44sNXUYpb+ULY9LzsCXO2/h1N3YctkvR7SIiIiIiEgl5eQKWHQ6FG/7zoROVgZiaprI1gkAJAAW/h2CrnYmCp9GyBEtIiIiIiJSSX4RSYhNzkBi9ZpySVYeAUBscgb8IpIUvm+OaBERUZnxIcdERFQZxadmKLRdSXBEi4iIiIiIVJKRrrZC25UEEy0iIiIiIlJJblYGMNXXRmF3X0nwrvqgm5WBwvfNRIuIiIiIiFSSmlSC+X3sACBfspX3fn4fu3J5nhYTLSIiIiIiUlk97E2x/lMnmOjLTw800dfG+k+dyu05WiyGQUREpCS8vb3h7e2NnJwcsUMhIlIqPexN0dXOBH4RSYhPzYCR7rvpguUxkpWHiRYREZGS8PLygpeXF1JSUqCvry92OERESkVNKoGHde0K2x+nDhbC29sbdnZ2cHV1FTsUIiIiIiJSMky0CuHl5YWQkBD4+/uLHQoRERERESkZJlpEREREREQKxkSLiIiIiIhIwZhoERERERERKRgTLSIiIiIiIgVjefciCIIAAEhJSSn1NnIzXysqHCIipVOW82feZ/POxfSOIq5NRERUciW5LjHRKkJqaioAwNzcXORIiIiUk/6qsm8jNTWVz416D69NRETiKs51SSLwa8IPys3NxbNnz6CrqwuJpORPjk5JSYG5uTmePHkCPT29coiwcmP/2f+q3H+Ax6Cs/RcEAampqahbty6kUs52z1PWa5Mqqur/14rC41M4HpvC8djkV5LrEke0iiCVSmFmZlbm7ejp6VXpH1D2n/2vyv0HeAzK0n+OZOWnqGuTKqrq/9eKwuNTOB6bwvHYyCvudYlfDxIRERERESkYEy0iIiIiIiIFY6JVzrS0tDB//nxoaWmJHYoo2H/2vyr3H+AxqOr9p4rDn7UP4/EpHI9N4XhsyobFMIiIiIiIiBSMI1pEREREREQKxkSLiIiIiIhIwZhoERERERERKRgTLSIiIiIiIgVjolUCly9fRp8+fVC3bl1IJBIcOXJEbr0gCFiwYAHq1q0LHR0ddOjQAffu3StyuwcPHoSdnR20tLRgZ2eHw4cPl1MPyqY8+r9p0ya0bdsWtWrVQq1atdClSxf4+fmVYy9Kr7z+/fPs3bsXEokE/fv3V2zgClRex+DVq1fw8vKCqakptLW1YWtrixMnTpRTL0qvvPq/atUqNG7cGDo6OjA3N8e0adOQkZFRTr0ovaL6f+jQIXTv3h2GhoaQSCQICgoq1naV5RxI4kpNTcXUqVNhYWEBHR0dtGrVCv7+/oW2P3ToELp27Yo6depAT08PHh4eOH36dAVGXHFKemzed+3aNairq8PR0bF8gxRRaY5PZmYm5syZAwsLC2hpacHa2hpbt26toIgrTmmOza5du+Dg4IBq1arB1NQUY8aMQWJiYgVFrFyYaJVAeno6HBwcsHbt2gLXL1u2DCtXrsTatWvh7+8PExMTdO3aFampqYVu09fXF0OHDsXIkSMRHByMkSNHYsiQIbhx40Z5daPUyqP/Fy9exCeffIJ//vkHvr6+qF+/Prp164anT5+WVzdKrTz6nycqKgozZsxA27ZtFR22QpXHMXj79i26du2KyMhIHDhwAA8fPsSmTZtQr1698upGqZVH/3ft2oVZs2Zh/vz5uH//PrZs2YJ9+/Zh9uzZ5dWNUiuq/+np6WjdujWWLl1a7G0q0zmQxPXZZ5/h7Nmz+OOPP3Dnzh1069YNXbp0KfR6cfnyZXTt2hUnTpxAQEAAOnbsiD59+iAwMLCCIy9/JT02eZKTkzFq1Ch07ty5giIVR2mOz5AhQ3D+/Hls2bIFDx8+xJ49e9CkSZMKjLpilPTYXL16FaNGjcK4ceNw79497N+/H/7+/vjss88qOHIlIVCpABAOHz4se5+bmyuYmJgIS5culS3LyMgQ9PX1hQ0bNhS6nSFDhgg9evSQW9a9e3dh2LBhCo9ZkRTV///Kzs4WdHV1hR07digyXIVTZP+zs7OF1q1bC5s3bxZGjx4t9OvXr5yiVixFHYP169cLDRo0EN6+fVue4Sqcovrv5eUldOrUSW7Z9OnThTZt2ig8ZkX6b//fFxERIQAQAgMDi9yOsp4DqWK9fv1aUFNTE44dOya33MHBQZgzZ06xt2NnZycsXLhQ0eGJqizHZujQocLcuXOF+fPnCw4ODuUYpXhKc3xOnjwp6OvrC4mJiRURomhKc2yWL18uNGjQQG7ZmjVrBDMzs3KLU5lxREtBIiIiEBcXh27dusmWaWlpoX379vDx8Sn0c76+vnKfAYDu3bt/8DOVUWn7/1+vX79GVlYWDAwMyiPMclOW/n///feoU6cOxo0bV95hlqvSHoO//voLHh4e8PLygrGxMezt7fHjjz8iJyenIsJWmNL2v02bNggICJBNmX38+DFOnDiBjz76qNxjrgxU5RxI5Ss7Oxs5OTnQ1taWW66jo4OrV68Waxu5ublITU1VuutLUUp7bLZt24bw8HDMnz+/vEMUVWmOz19//QUXFxcsW7YM9erVg42NDWbMmIE3b95URMgVpjTHplWrVoiJicGJEycgCAKeP3+OAwcOVJlrVkmpix2AqoiLiwMAGBsbyy03NjZGVFTUBz9X0GfytqcsStv//5o1axbq1auHLl26KDS+8lba/l+7dg1btmwp9r0slVlpj8Hjx49x4cIFjBgxAidOnEBoaCi8vLyQnZ2NefPmlWvMilTa/g8bNgwvXrxAmzZtIAgCsrOz8eWXX2LWrFnlGm9loSrnQCpfurq68PDwwKJFi2BrawtjY2Ps2bMHN27cQKNGjYq1jZ9//hnp6ekYMmRIOUdbsUpzbEJDQzFr1ixcuXIF6uqq/atgaY7P48ePcfXqVWhra+Pw4cNISEjAxIkTkZSUpFL3aZXm2LRq1Qq7du3C0KFDkZGRgezsbPTt2xe//vprBUevHDiipWASiUTuvSAI+ZYp4jOVVVn6smzZMuzZsweHDh3K9+2KsihJ/1NTU/Hpp59i06ZNMDQ0rIjwKkRJfwZyc3NhZGSEjRs3wtnZGcOGDcOcOXOwfv368g61XJS0/xcvXsQPP/yAdevW4datWzh06BCOHTuGRYsWlXeolYYqnQOp/Pzxxx8QBAH16tWDlpYW1qxZg+HDh0NNTa3Iz+7ZswcLFizAvn37YGRkVAHRVqySHJucnBwMHz4cCxcuhI2NjQjRVryS/uzk5uZCIpFg165dcHNzQ69evbBy5Ups375d5Ua1SnpsQkJCMHnyZMybNw8BAQE4deoUIiIiMGHChAqOXDmo9tcYFcjExATAu29nTU1NZcvj4+PzfVv738/995vboj5TGZW2/3lWrFiBH3/8EefOnUPz5s3LLc7yUpr+h4eHIzIyEn369JEty83NBQCoq6vj4cOHsLa2LseoFau0PwOmpqbQ0NCQO6nb2toiLi4Ob9++haamZvkFrUCl7f93332HkSNHym4kbtasGdLT0/H5559jzpw5kEpV+/swVTkHUvmztrbGpUuXkJ6ejpSUFJiammLo0KGwsrL64Of27duHcePGYf/+/Uo3W6K4SnJsUlNTcfPmTQQGBuKrr74C8O7aIwgC1NXVcebMGXTq1Kmiu1CuSvqzY2pqinr16kFfX1+2zNbWFoIgICYmptijqMqgpMdmyZIlaN26Nb755hsAQPPmzVG9enW0bdsWixcvlrv+EUe0FMbKygomJiY4e/asbNnbt29x6dIltGrVqtDPeXh4yH0GAM6cOfPBz1RGpe0/ACxfvhyLFi3CqVOn4OLiUt6hlovS9L9Jkya4c+cOgoKCZK++ffuiY8eOCAoKgrm5eUWFrxCl/Rlo3bo1wsLCZEkmADx69AimpqZKk2QBpe//69ev8yVTampqEAQBgiCUW7yVhaqcA6niVK9eHaampnj58iVOnz6Nfv36Fdp2z5498PT0xO7du6vEPSTFOTZ6enr5rj0TJkxA48aNERQUhJYtW4oQecUo7s9O69at8ezZM6SlpcmWPXr0CFKpFGZmZhUVboUq7rEp7JoFoEpcs0qsoqtvKLPU1FQhMDBQCAwMFAAIK1euFAIDA4WoqChBEARh6dKlgr6+vnDo0CHhzp07wieffCKYmpoKKSkpsm2MHDlSmDVrluz9tWvXBDU1NWHp0qXC/fv3haVLlwrq6urC9evXK7x/RSmP/v/000+CpqamcODAASE2Nlb2Sk1NrfD+FaU8+v9flb3qYHkcg+joaKFGjRrCV199JTx8+FA4duyYYGRkJCxevLjC+1eU8uj//PnzBV1dXWHPnj3C48ePhTNnzgjW1tbCkCFDKrx/RSmq/4mJiUJgYKBw/PhxAYCwd+9eITAwUIiNjZVtQ5nPgSSuU6dOCSdPnpT9P3FwcBDc3NxkFUtnzZoljBw5UtZ+9+7dgrq6uuDt7S13fXn16pVYXSg3JT02/6XKVQcFoeTHJzU1VTAzMxMGDx4s3Lt3T7h06ZLQqFEj4bPPPhOrC+WmpMdm27Ztgrq6urBu3TohPDxcuHr1quDi4iK4ubmJ1YVKjYlWCfzzzz8CgHyv0aNHC4Lwrrzz/PnzBRMTE0FLS0to166dcOfOHblttG/fXtY+z/79+4XGjRsLGhoaQpMmTYSDBw9WUI9Kpjz6b2FhUeA258+fX3EdK6by+vd/X2VPtMrrGPj4+AgtW7YUtLS0hAYNGgg//PCDkJ2dXUG9Kr7y6H9WVpawYMECwdraWtDW1hbMzc2FiRMnCi9fvqy4jhVTUf3ftm1bkf+flfkcSOLat2+f0KBBA0FTU1MwMTERvLy85JKm0aNHC+3bt5e9b9++/Qd/XlVJSY/Nf6l6olWa43P//n2hS5cugo6OjmBmZiZMnz5deP36dQVHXv5Kc2zWrFkj2NnZCTo6OoKpqakwYsQIISYmpoIjVw4SQeA4HxERERERkSLxHi0iIiIiIiIFY6JFRERERESkYEy0iIiIiIiIFIyJFhERERERkYIx0SIiIiIiIlIwJlpEREREREQKxkSLiIiIiIhIwZhoERERERERKRgTLSIiIiIiIgVjokUkgsjISEgkEnh6eoodChEREa9LROWAiRZRFda8eXNIJBJIJBJcuXKlwDYvX76EoaGhrN2jR48qOEoiIqoqeF0iVcJEi6iKysjIwP3796Gurg4AuHPnToHt5s2bh5cvXwIAdHV10ahRowqLkYiIqg5el0jVMNEiqqKCg4ORnZ2Nnj17QiqVFnhBu3v3LtavX49evXoBABwdHSGRSCo6VCIiqgJ4XSJVw0SLqJLZsWMH3N3dUaNGDdSoUQPu7u7YsWNHgW2zs7OxZMkSWFtbQ1tbGw0bNsSSJUvw+PHjIufa37p1CwDQrl07WFtbF3hBmzp1KnR1dTFgwAAAgLOzc9k7SERESoXXJaLSURc7ACL617Rp07Bq1SrUq1cP48aNg0QiwcGDB+Hp6Yng4GCsXLlSrv3YsWPxxx9/wNraGl5eXsjMzMSqVavg6+tb5L4CAgIAAE5OTnB0dMTp06fl1h86dAjnz5/H6tWrERUVJWtLRERVB69LRGUgEFGFi4iIEAAIo0ePli27fPmyAECwtbUVXr16JVv+6tUroUmTJgIA4cqVK7Ll586dEwAILi4uwuvXr2XLY2NjBRMTk3zb/68WLVoIAISXL18KP/zwgwBAiIqKEgRBEDIyMgQrKyvB1tZWyMrKEnr37i0AEO7evau4g0BERJUGr0tEisepg0SVxPbt2wEACxYsgL6+vmy5vr4+5s+fL9cGAHbu3AkA+O6776CjoyNbbmJigilTpnxwX2/fvsXdu3fRoEED1KxZE46OjgD+vfF4xYoViIiIwKpVq6Curo6AgABUq1YNTZo0kdvOhAkTWAqYiEhFKeN1iagyYaJFVEkEBgYCADp06JBvXd6yoKAg2bLg4GAAQKtWrfK1L2jZ+27fvo2srCzZlAsHBwcA7y5oT58+xZIlS9C3b19069YNcXFxiI2NhYODA9TU1PJtp3nz5sXqHxERKRdlvC4RVSZMtIgqiZSUFEilUtSpUyffOmNjY0ilUiQnJ+drX7t27QLbf0jeDcd5NxHXq1cPderUwZ07dzBz5kxkZWXh559/BiA/Z/59giDgzp07soshERGpFmW7LhFVNky0iCoJPT095Obm4sWLF/nWxcfHIzc3F3p6evnaJyYm5mv//PnzD+6roIuUg4MDTp48id27d2PatGlo2LAhgH8vfv+9oIWHhyMtLQ1v3rxBmzZtUK1aNTg7O+Px48fF7DEREVVmynZdIqpsmGgRVRItWrQAAFy8eDHfukuXLgGAbM468O+0Ch8fn3ztC1r2voIuUo6Ojnj58iVMTEwwZ84c2fK8i99/S+gGBwdDTU0Na9euxbJly+Dv74/s7GwsW7bsg/smIiLloGzXpYIIgoC///67yHZE5YGJFlElMXr0aADAwoULkZKSIluekpKChQsXyrUBgBEjRgAAFi1ahIyMDNnyuLg4rF69utD9ZGVl4c6dO6hfvz4MDQ1lyydOnIjDhw/jzJkz0NXVlS2/desWtLS0YGdnJ7ed4OBgGBsb488//0SrVq3QtGlTdOnSBQkJCaXpPhERVTLKdl367zYTEhJw/vx5rFixAgkJCXj16lUxe06kGHyOFlEl0a5dO0yaNAm//vor7O3tMWjQIAiCgEOHDuHJkyeYPHky2rVrJ2vfpUsXjBgxArt27UKzZs3Qr18/ZGZm4s8//0TLli3x999/QyrN/13KvXv3kJmZmW/KhZWVFaysrOSWJSQk4MmTJ3BxcYGGhobcuuDgYAwdOlRu2khkZKRsagcRESk3Zbsu/deyZcuwYsUKCIIAd3d37N27Fy4uLmU4IkQlwxEtokpkzZo12Lp1K0xMTLBx40Zs2rQJJiYm2Lp1a4HfBm7fvh2LFi1CTk4Ofv31V5w4cQJTp07F3LlzAUAuCcpTkpuIP9Q2ODgYLVu2lFsWFBTE4hhERCpEma5L79PQ0MCyZctgYWEBPT09zJo1i0kWVTiJIAiC2EEQkWJt3rwZ48ePx7p16/Dll18qfPvJycmoWbMm7t+/L3uGSUpKCmrWrInbt2/D3t5e4fskIiLlVd7XpYKEhobi5s2b6NatG/78888K2y9RHiZaREosLi4OxsbGkEgksmVPnz5F69atERMTg4iICJibmyt8v5cvX0aPHj2Qmpoqe4bJlStX0LVrV6SlpUFdnbOSiYiqIrGuS0SVEX8bIlJiS5cuxfHjx9G2bVsYGRkhOjoax44dQ2pqKhYsWFBuF7Pg4GA0bdpU7kGRwcHBsLW1ZZJFRFSFiXVdIqqMOKJFpMROnTqFlStXIjg4GC9fvoS2tjaaN2+OiRMnYvjw4WKHR0REVQyvS0T/YqJFRERERESkYKw6SEREREREpGBMtIiIiIiIiBSMiRYREREREZGCMdEiIiIiIiJSMCZaRERERERECsZEi4iIiIiISMGYaBERERERESkYEy0iIiIiIiIFY6JFRERERESkYEy0iIiIiIiIFIyJFhERERERkYL9HyZ0/blQc9NvAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "axes = plot_hmf_and_smf(true_smf, label=\"Truth\")\n",
    "axes = plot_hmf_and_smf(model.calc_sumstats_from_params(init_params),\n",
    "                        \"k--\", label=\"Initial guess\", axes=axes)\n",
    "axes = plot_hmf_and_smf(model.calc_sumstats_from_params(results.x),\n",
    "                        \"r--\", label=\"Best fit\", axes=axes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The same thing, but parallel!\n",
    "\n",
    "The power of multigrad is that it sums summary statistics over all given MPI ranks. Note that our halo loading function will automatically distribute the halos evenly across all available ranks, which is what makes this sum mathematically valid. So let's redo this whole fit, but this time in the script `smf_grad_descent.py` which is copied below (this is essentially just copy-and-pasting everything we've done in this notebook):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```python\n",
    "\"\"\"\n",
    "smf_grad_descent.py\n",
    "\"\"\"\n",
    "\n",
    "from mpi4py import MPI\n",
    "import jax.scipy\n",
    "from jax import numpy as jnp\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from diffopt import multigrad\n",
    "\n",
    "\n",
    "def load_halo_masses(num_halos=10_000, comm=MPI.COMM_WORLD):\n",
    "    # Generate fake halo masses between 10^10 < M_h < 10^11 as a power law\n",
    "    quantile = jnp.linspace(0, 0.9, num_halos)\n",
    "    mhalo = 1e10 / (1 - quantile)\n",
    "\n",
    "    # Assign halos evenly across given MPI ranks\n",
    "    return np.array_split(mhalo, comm.size)[comm.rank]\n",
    "\n",
    "\n",
    "# Compute one bin of the stellar mass function (SMF)\n",
    "@jax.jit\n",
    "def calc_smf_bin(params, logsm_low, logsm_high, volume, log_halo_masses):\n",
    "    # Unpack model parameters f and sigma\n",
    "    log_f, log_sigma = params\n",
    "    mean_logsm = log_f + log_halo_masses\n",
    "    sigma = 10 ** log_sigma\n",
    "\n",
    "    # Integrating the log-normal PDF over the bin, we get erf:\n",
    "    erf_high = 0.5 * (1 + jax.scipy.special.erf(\n",
    "        (logsm_high - mean_logsm)/(jnp.sqrt(2) * sigma)))\n",
    "    erf_low = 0.5 * (1 + jax.scipy.special.erf(\n",
    "        (logsm_low - mean_logsm)/(jnp.sqrt(2) * sigma)))\n",
    "    prob_in_bin = erf_high - erf_low\n",
    "\n",
    "    # Sum probabilities, convert to number density, divide by bin width\n",
    "    return jnp.sum(prob_in_bin) / volume / (logsm_high - logsm_low)\n",
    "\n",
    "\n",
    "# Compute the stellar mass function over all bins (loop over calc_smf_bin)\n",
    "@jax.jit\n",
    "def calc_smf(params, smf_bin_edges, volume, log_halo_masses):\n",
    "    smf = []\n",
    "    logsm_low = smf_bin_edges[0]\n",
    "    for logsm_high in smf_bin_edges[1:]:\n",
    "        smf_bin = calc_smf_bin(\n",
    "            params, logsm_low, logsm_high, volume, log_halo_masses)\n",
    "        smf.append(smf_bin)\n",
    "        logsm_low = logsm_high\n",
    "    return jnp.array(smf)\n",
    "\n",
    "\n",
    "def plot_hmf_and_smf(smf, logmh_per_rank=None, plotarg=\"C0o-\",\n",
    "                     label=None, axes=None):\n",
    "    if axes is None:\n",
    "        _, axes = plt.subplots(ncols=2, figsize=(10.5, 4))\n",
    "    if logmh_per_rank is not None:\n",
    "        colors = [f\"C{i}\" for i in range(len(logmh_per_rank))]\n",
    "        axes[0].hist(logmh_per_rank, bins=np.linspace(10, 11, 50),\n",
    "                     color=colors, stacked=True)\n",
    "        for i in range(len(colors)):\n",
    "            axes[0].hist([], bins=np.linspace(10, 11, 50),\n",
    "                         color=f\"C{i}\", label=f\"MPI Rank = {i}\")\n",
    "        axes[0].legend(frameon=False)\n",
    "        axes[0].semilogy()\n",
    "        axes[0].set_xlabel(\"$\\\\log M_h$\", fontsize=14)\n",
    "        axes[0].set_ylabel(\"$\\\\Phi(\\\\log M_h)$\", fontsize=14)\n",
    "\n",
    "    smf_bin_cens = 0.5 * (smf_bin_edges[:-1] + smf_bin_edges[1:])\n",
    "    axes[1].semilogy(smf_bin_cens, smf, plotarg, label=label)\n",
    "    if label is not None:\n",
    "        axes[1].legend(frameon=False)\n",
    "    axes[1].set_xlabel(\"$\\\\log M_\\\\star$\", fontsize=14)\n",
    "    axes[1].set_ylabel(\"$\\\\Phi(\\\\log M_\\\\star)$\", fontsize=14)\n",
    "    return axes\n",
    "\n",
    "\n",
    "class MySMFModel(multigrad.OnePointModel):\n",
    "    def calc_partial_sumstats_from_params(self, params):\n",
    "        # Accessing global variables is fine, but I prefer to store them in\n",
    "        # the `aux_data` attribute, which we will define during construction\n",
    "        bin_edges = self.aux_data[\"smf_bin_edges\"]\n",
    "        volume = self.aux_data[\"volume\"]\n",
    "        log_halo_masses = self.aux_data[\"log_halo_masses\"]\n",
    "\n",
    "        return calc_smf(params, bin_edges, volume, log_halo_masses)\n",
    "\n",
    "    def calc_loss_from_sumstats(self, sumstats):\n",
    "        # Add 1e-10 so that log values always remain finite\n",
    "        target_sumstats = jnp.log10(self.aux_data[\"target_sumstats\"] + 1e-10)\n",
    "        sumstats = jnp.log10(sumstats + 1e-10)\n",
    "        # Reduced chi2 loss function assuming unit errors (mean squared error)\n",
    "        return jnp.mean((sumstats - target_sumstats)**2)\n",
    "\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    volume = 1.0\n",
    "    smf_bin_edges = jnp.linspace(9, 10, 11)\n",
    "    true_params = jnp.array([-2.0, -0.5])\n",
    "\n",
    "    log_halo_masses = jnp.log10(load_halo_masses(10_000))\n",
    "    logmh_per_rank = MPI.COMM_WORLD.allgather(log_halo_masses)  # for plotting\n",
    "\n",
    "    # We must sum calc_smf over all MPI ranks this time\n",
    "    # Could equivalently use model.calc_sumstats_from_params(true_params)\n",
    "    true_smf = multigrad.reduce_sum(\n",
    "        calc_smf(true_params, smf_bin_edges, volume, log_halo_masses))\n",
    "\n",
    "    aux_data = dict(\n",
    "        log_halo_masses=log_halo_masses,\n",
    "        smf_bin_edges=smf_bin_edges,\n",
    "        volume=volume,\n",
    "        target_sumstats=true_smf  # SMF at truth: params=(-2.0, -0.5)\n",
    "    )\n",
    "    model = MySMFModel(aux_data=aux_data)\n",
    "\n",
    "    # Initial guess for our parameters. If it's too far off, there is always a\n",
    "    # risk of getting stuck in local minima or other zero-valued gradients\n",
    "    init_params = true_params + jnp.array([-1.5, 0.7])\n",
    "\n",
    "    # Run gradient descent using the BFGS method powered by scipy\n",
    "    _, _, results = model.run_bfgs(init_params)\n",
    "\n",
    "    init_smf = model.calc_sumstats_from_params(init_params)\n",
    "    final_smf = model.calc_sumstats_from_params(results.x)\n",
    "    # Print and plot results from the root rank only\n",
    "    if not MPI.COMM_WORLD.rank:\n",
    "        print(\"BGFS has converged:\", results.success, flush=True)\n",
    "        print(\"Initial guess =\", init_params, flush=True)\n",
    "        print(\"True params =\", true_params, flush=True)\n",
    "        print(\"Converged params =\", results.x, flush=True)\n",
    "        print(\"\\nFull results info:\", flush=True)\n",
    "        print(results, flush=True)\n",
    "\n",
    "        axes = plot_hmf_and_smf(\n",
    "            true_smf, logmh_per_rank, label=\"Truth\")\n",
    "        axes = plot_hmf_and_smf(\n",
    "            init_smf, None, \"k--\", label=\"Initial guess\", axes=axes)\n",
    "        axes = plot_hmf_and_smf(\n",
    "            final_smf, None, \"r--\", label=\"Best fit\", axes=axes)\n",
    "\n",
    "        plt.savefig(\"smf_gradient_descent.png\", bbox_inches=\"tight\")\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's run this with three MPI ranks. Executing `mpiexec -n 3 python smf_gradient_descent.py` yields approximately the same results down to machine precision, but with each rank only loading 1/3 of the halo data. There is no speedup observed here due to the computation being very cheap with a small dataset of 10,000 halos, but if we increased this significantly, we would see a 3x speedup.\n",
    "\n",
    "```txt\n",
    "BFGS Gradient Descent Progress:  16%|█▌        | 16/100 [00:03<00:15,  5.26it/s]\n",
    "BGFS has converged: True\n",
    "Initial guess = [-3.5         0.19999999]\n",
    "True params = [-2.  -0.5]\n",
    "Converged params = [-2.00000006 -0.49999987]\n",
    "\n",
    "Full results info:\n",
    "  message: CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL\n",
    "  success: True\n",
    "   status: 0\n",
    "      fun: 4.862954657014473e-12\n",
    "        x: [-2.000e+00 -5.000e-01]\n",
    "      nit: 16\n",
    "      jac: [ 5.082e-06  9.755e-06]\n",
    "     nfev: 29\n",
    "     njev: 29\n",
    " hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>\n",
    "```\n",
    "![Results plot](./smf_gradient_descent.png)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "main311",
   "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.11.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}