{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e89bd8d4",
   "metadata": {},
   "source": [
    "# 4. Predict observables and generate mocks\n",
    "\n",
    "This notebook defines its own JAX model and physical parameters below.\n",
    "It builds on the [model tutorial](models.ipynb#compose-the-jax-model). Always distinguish a\n",
    "projected observable from an intrinsic quantity and a dispersion from a variance.\n",
    "\n",
    "\n",
    "This notebook is self-contained. Install JeansPy with the `numpyro_cpu` and\n",
    "`plotting` extras as described in the installation guide, then select that\n",
    "environment as your Jupyter kernel and run cells from top to bottom.\n",
    "Saved outputs are an example run; timings and short-chain results can vary."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "923a4db1",
   "metadata": {},
   "source": [
    "## Setup\n",
    "\n",
    "Run this setup in a fresh kernel. It defines all objects used below.\n",
    "\n",
    "The generalized NFW halo uses `ZhaoModel` with fixed `alpha=1`, `beta=3`\n",
    "and an untruncated cutoff. Its inner slope `gamma` is sampled in the\n",
    "inference examples; `truth[\"gamma\"]=1` generates standard-NFW mocks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "0a75e27d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "os.environ.setdefault(\"JEANSPY_JAX_PLATFORM\", \"cpu\")\n",
    "os.environ.setdefault(\"JEANSPY_JAX_ENABLE_X64\", \"true\")\n",
    "\n",
    "# The example uses no progress widgets; ignore only their optional-import warning.\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\", message=\"IProgress not found.*\")\n",
    "\n",
    "# Import JeansPy before JAX so its runtime settings take effect.\n",
    "from jeanspy.model_jax import (\n",
    "    DSphModel, PlummerModel, ZhaoModel, ConstantAnisotropyModel,\n",
    ")\n",
    "from jeanspy.sampler_numpyro import JeansLikelihoodModel, ParameterSpec\n",
    "import jax\n",
    "import numpyro.distributions as dist\n",
    "from numpyro.infer import MCMC, NUTS, init_to_value\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "model = DSphModel(submodels={\n",
    "    \"StellarModel\": PlummerModel(),\n",
    "    \"DMModel\": ZhaoModel(),\n",
    "    \"AnisotropyModel\": ConstantAnisotropyModel(),\n",
    "})\n",
    "fixed = dict(re_pc=200., alpha=1., beta=3., r_t_pc=np.inf)\n",
    "truth = dict(**fixed, rs_pc=500., rhos_Msunpc3=.1, gamma=1., beta_ani=0.)\n",
    "options = dict(solver=\"kernel\", n_u=128, n_kernel=32,\n",
    "               dm_mass_method=\"numeric\", dm_mass_n_steps=128)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a83a124c",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(123)\n",
    "u = rng.uniform(size=32)\n",
    "R_pc = 200. * np.sqrt(u / (1. - u))  # Projected Plummer radii.\n",
    "e_vlos_kms = np.full(32, 2.)\n",
    "variance = np.asarray(model.sigmalos2(R_pc, params=truth, **options))\n",
    "vlos_kms = rng.normal(0., np.sqrt(variance + e_vlos_kms**2))\n",
    "data = dict(R_pc=R_pc, vlos_kms=vlos_kms, e_vlos_kms=e_vlos_kms)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c247b63c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import jax.numpy as jnp\n",
    "params = dict(truth)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b8ed0cf",
   "metadata": {},
   "source": [
    "## Calculate the LOS dispersion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "3ad41af3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import jax.numpy as jnp\n",
    "\n",
    "variance = model.sigmalos2(R_pc, params=params, solver=\"kernel\",\n",
    "                          n_u=256, n_kernel=64)\n",
    "sigma_kms = jnp.sqrt(variance)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2a602f3",
   "metadata": {},
   "source": [
    "`sigmalos2` returns the intrinsic LOS second central moment in (km/s)².\n",
    "The mean velocity does not belong inside its square root. NumPy users can\n",
    "call `model.sigmalos(R_pc, n=256, n_kernel=64)` directly on a NumPy/SciPy model;\n",
    "JAX users take the square root as above. The different keyword names are\n",
    "documented in the individual\n",
    "[NumPy](https://gomeshun.github.io/jeanspy/dev/api/all.html) and\n",
    "[JAX](https://gomeshun.github.io/jeanspy/dev/api/all.html) API entries.\n",
    "\n",
    "## Plot a radial profile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "7dec8a8f",
   "metadata": {
    "mystnb": {
     "image": {
      "alt": "Predicted line-of-sight dispersion versus projected radius in pc."
     }
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAG1CAYAAADwRl5QAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAT1NJREFUeJzt3XlcFPXjBvBndtnlvm8EEUVFkVtRNG/N1ExLTc0DUyuPNKys+HZoh1l9O/x1mWZhnqV5puaZ5n2h4IUoguDBoSIs5wK78/vD2m/kBbgwezzv12tftbMzs8/qBk+fmfmMIIqiCCIiIiITJJM6ABEREVF9YdEhIiIik8WiQ0RERCaLRYeIiIhMFosOERERmSwWHSIiIjJZLDpERERkslh0iIiIyGRZSB1ASlqtFteuXYO9vT0EQZA6DhEREdWAKIooKiqCj48PZLL7j9mYddG5du0a/Pz8pI5BREREdXD58mX4+vredx2zLjr29vYAbv9BOTg4SJyGiIiIakKlUsHPz0/3e/x+zLro/H24ysHBgUWHiIjIyNTktBOejExEREQmi0WHiIiITBaLDhEREZksFh0iIiIyWSw6REREZLJYdIiIiMhksegQERGRyWLRISIiIpPFokNEREQmi0WHiIiITBaLDhEREZksgyk6RUVFiIuLg7+/P6ytrdGxY0ccPXr0nuvv3r0bgiDc8cjJyWnA1ERERGTIDOamnhMmTMDp06exZMkS+Pj4YOnSpejVqxfOnj2LRo0a3XO71NTUajfk9PDwaIi497X3wnV89UcamrnbIsDNFk3d7NDU3RZ+LjZQyA2mWxIREZk8gyg6ZWVlWL16NdavX48uXboAAGbNmoXffvsN8+bNwwcffHDPbT08PODk5NRASWvmzDUVjmTk40hGfrXlFjIBAW62aO5ph0APezT3sENLL3sEuNmyABEREdUDgyg6VVVV0Gg0sLKyqrbc2toa+/btu++24eHhUKvVaNOmDWbNmoVOnTrdc121Wg21Wq17rlKpHi74PfRr4w0vByukXy/GxRslyLhegowbJSir1OBCXjEu5BUD+N8hNqVchmYedmjlZY8gb3sE+zgi2McBTjbKeslHRERkLgRRFEWpQwBAx44doVQqsXz5cnh6emLFihWIjY1FYGAgUlNT71g/NTUVu3fvRtu2baFWq7Fw4UIsWbIEhw8fRmRk5F3fY9asWXj33XfvWF5YWFjt8Fd90GpFXCssQ1pese5xPrcIqTlFKKnQ3HUbX2drBPs4INTXCWG+TgjxdYSjtaJecxIRERk6lUoFR0fHGv3+Npiic/HiRYwbNw579uyBXC5HZGQkWrRogcTERKSkpNRoH127dkXjxo2xZMmSu75+txEdPz+/Bik696LVirhaUIZzOUU4l63C2WwVzlxTISu/9K7rN3W3RbivEyL9nRHl74wWnvaQy4QGTk1ERCQdoyw6fyspKYFKpYK3tzeGDRuG4uJibNq0qUbbzpgxA/v27cPBgwdrtH5t/qAaWmFZJc5eU+H01UKcvFqI5MsFdy0/9pYWiPB3Rjt/Z3Ro5opQX0dYWsglSExERNQwavP72yDO0fknW1tb2Nra4tatW9i6dSs++eSTGm+blJQEb2/vekzXcBytFYhp5oqYZq66ZfklFUi+UoATmbeQmHULJ7IKUKSuwp7z17Hn/HVgO2BpIUOUvzM6NHVFp0BXhPk6wYInOhMRkZkymBGdrVu3QhRFtGzZEmlpaZgxYwasrKywd+9eKBQKxMfH4+rVq1i8eDEAYO7cuQgICEBwcDDKy8uxcOFCfPXVV9i2bRt69uxZo/c05BGdmqjSaHEupwiJmbdwJCMfhzNu4kZxRbV17C0t0L6pKzo3d0OXFu4IcLOVKC0REZF+GOWITmFhIeLj43HlyhW4uLhg8ODBmD17NhSK2yffZmdnIysrS7d+RUUFXnnlFVy9ehU2NjYIDQ3Fjh070L17d6k+QoOzkMvQppEj2jRyRGzHJhBFERevF+NQej4OXryJAxdv4FZpJXak5GJHSi4AoImrDbq19ECPIA9EB7jASsHDXEREZLoMZkRHCsY+ovMgWq2IM9dU2Jt2Hfsu3MDRS/mo1Pzvr9taIUeXFm7oE+yFHkEevJydiIiMglGfjNyQTL3o/Fuxugr7025g17k87ErNQ67qf1egyWUC2ge44LE2XnisjRc87K3usyciIiLpsOjUkLkVnX8SxdujPdvO5GDb2VycyynSvSYTgPYBrng8zBt923jDxZYjPUREZDhYdGrInIvOv126UYJtZ3Ow+VQOki4X6JbLZQK6NHfDU5G+6N3ak+f0EBGR5Fh0aohF5+4u55di06ls/JZ8DWeu/e82GfaWFugX4o2nIhshOsAFgsCJComIqOGx6NQQi86DXbxejHUnrmLN8au4WlCmWx7gZoth7fwwJMoXbnaWEiYkIiJzw6JTQyw6NafVijhyKR9rjl/BppPZuvtzWcgE9G7tiZHt/dEp0JWjPEREVO9YdGqIRaduStRV+C35GlYcvYzkf5zP08zdFrEdm+CpSF/YWRrMFE1ERGRiWHRqiEXn4aVkq7D8cBbWHL+iG+Wxt7TA4ChfjOsUgMauNhInJCIiU8OiU0MsOvpTVF6J1YlXsPhgJtJvlAC4fZl63zbeeL5LU4T5OUkbkIiITAaLTg2x6OifVitib9oN/LgvA3+ev65b3j7ABRO7NUO3Fu48j4eIiB4Ki04NsejUr3M5KizYk44NSddQpb39NQv1dcSL3QPRu7UnCw8REdUJi04Nseg0jOzCMvywNwPLDmehrPL2eTxBXvaY2qM5+rbxgkzGwkNERDXHolNDLDoN62axGj/sy8BPBy7pTlwO8rLHjD4t0SPIgyM8RERUIyw6NcSiI42C0gr8uP8SEvZloEhdBQCI8nfGjD4t0aGpq8TpiIjI0LHo1BCLjrQKSisw78+L+OnAJZRXagEAXVq44z/9ghDkxb8PIiK6OxadGmLRMQy5qnJ8/UcaVhzJQpVWhEwAhrVrjJd7t4C7PW8vQURE1bHo1BCLjmHJvFmCj7ecw+ZTOQAAW6Uck7sHYvwjAbxrOhER6bDo1BCLjmE6eikfH2w8i+QrhQAAPxdrzHw8GL1ae0qcjIiIDAGLTg2x6BgurVbEbyev4aPfzyG7sBwA0CPIAzMHtIa/q63E6YiISEosOjXEomP4StRV+OqPNPywLx2VGhFKCxkmdm2Gyd2a8XAWEZGZqs3vb1kDZSKqE1tLC7zRNwi/v9QFjwS6oaJKiy93XkC/L/fiSEa+1PGIiMjAseiQUQj0sMOS8dH4+pkIuNtbIv16CZ6efxD/WXsKqvJKqeMREZGBYtEhoyEIAh4P9cGO6V0xvJ0fAGD54Sz0/vxPbD+bK3E6IiIyRCw6ZHQcbRT4aHAoVjzXAQFutshVqfHc4mN4dVUyR3eIiKgaFh0yWjHNXPH7S53xQpemEATg18Qr6PPFHuy7cEPqaEREZCBYdMioWSnkiO/XCqteiIG/qw2yC8sx6ofDeHvdaZT9deNQIiIyXyw6ZBLaNnHB7y91xpgYfwDAkkOZGPD1Ppy9ppI4GRERSYlFh0yGjdIC7w1sg6Xj28PD3hJpecUY9M1+JOzPgBlPF0VEZNZYdMjkPNLcDVviuqBXKw9UaLR497ezGLfoKG4Uq6WORkREDYxFh0ySi60S349pi/cGBkNpIcOu1Ovo/+VeHL3ESQaJiMwJiw6ZLEEQMCamCTa82AmBHnbIVakxfMEhLNhzkYeyiIjMBIsOmbwgLwesn9IJA8N9oNGK+HDzObywJBGFZZxzh4jI1LHokFmwtbTA3GHh+GBQGyjlMmw7m4sBX+1DSjavyiIiMmUsOmQ2BEHAqA7++HVSDBo5WSMrvxRPfXsAm05mSx2NiIjqCYsOmZ1QXydsmvYIOjd3Q1mlBlOWH8enW1Oh1fK8HSIiU8OiQ2bJyUaJhLHt8FznAADA17vS8NziY7xXFhGRiWHRIbNlIZfhzf6t8cWwMCgtZNh5Lg9PfXsAWTdLpY5GRER6wqJDZu/JCF/8OjEGXg5WSMsrxpPf7kdiJufbISIyBSw6RLh93s76FzuhTSMH3CypwIjvD2ND8jWpYxER0UNi0SH6i6eDFVa+EINerTxRUaXFtBUn8NXOC5xckIjIiLHoEP2DjdIC80dHYcIjt09S/mz7eby++iSqNFqJkxERUV2w6BD9i1wm4K3HW+ODQW0gE4CVx67ghSWJKKvQSB2NiIhqiUWH6B5GdfDH/NFtYfnXFVnPLDyE/JIKqWMREVEtsOgQ3Ufv1p5YNqE9HK0VOJFVgCHfHcDlfF5+TkRkLFh0iB6gbRMX/DoxBt6OVki/XoLB8w7gQm6R1LGIiKgGWHSIaqC5pz3WTO6IFp52yCtS4+n5B3HqSqHUsYiI6AEMpugUFRUhLi4O/v7+sLa2RseOHXH06NH7brN7925ERkbC0tISgYGBWLRoUcOEJbPk7WiNX56PQaivI26VVuKZ7w/hSAYnFiQiMmQGU3QmTJiA7du3Y8mSJTh16hQeffRR9OrVC1evXr3r+hkZGejfvz+6d++OpKQkxMXFYcKECdi6dWsDJydz4myrxLIJ7REd4IIidRXG/HgYu1PzpI5FRET3IIgGMBtaWVkZ7O3tsX79evTv31+3PCoqCn379sUHH3xwxzavv/46Nm3ahNOnT+uWDR8+HAUFBdiyZUuN3lelUsHR0RGFhYVwcHB4+A9CZqO8UoNJSxOxK/U6FHIBX42IxGNtvKSORURkFmrz+9sgRnSqqqqg0WhgZWVVbbm1tTX27dt3120OHjyIXr16VVvWp08fHDx48J7vo1aroVKpqj2I6sJKIcf80W3RP8QblRoRU5Yfx+ZT2VLHIiKifzGIomNvb4+YmBi8//77uHbtGjQaDZYuXYqDBw8iO/vuvzxycnLg6elZbZmnpydUKhXKysruus2cOXPg6Oioe/j5+en9s5D5UFrI8OWICDwV0QgarYipK07gN94fi4jIoBhE0QGAJUuWQBRFNGrUCJaWlvjyyy8xYsQIyGT6ixgfH4/CwkLd4/Lly3rbN5knuUzAf4eGYXCkLzRaES/9fALrk+5+XhkRETU8C6kD/K1Zs2b4888/UVJSApVKBW9vbwwbNgxNmza96/peXl7Izc2ttiw3NxcODg6wtra+6zaWlpawtLTUe3Yyb3KZgE+GhEImAKsSr2D6L0kQRWBQRCOpoxERmT2DGdH5m62tLby9vXHr1i1s3boVAwcOvOt6MTEx2LlzZ7Vl27dvR0xMTEPEJKpGLhPw8eBQDG/nB60ITF+ZxJEdIiIDYDBFZ+vWrdiyZQsyMjKwfft2dO/eHUFBQXj22WcB3D7sNGbMGN36EydORHp6Ol577TWcO3cO3377LVauXInp06dL9RHIzMlkAj58MgQjohtDFIGXVyZjy2meoExEJCWDKTqFhYWYMmUKgoKCMGbMGDzyyCPYunUrFAoFACA7OxtZWVm69QMCArBp0yZs374dYWFh+Oyzz7Bw4UL06dNHqo9ABJlMwOxBbXTn7ExdcQI7U3IfvCEREdULg5hHRyqcR4fqi0YrIu6XJPyWfA1KuQwLY9uiSwt3qWMREZkEo5tHh8jUyGUCPn86DH2CPVGh0eL5JcdwKP2m1LGIiMwOiw5RPVHIZfhqRCR6BHmgvFKL8YuO4uSVAqljERGZFRYdonqktJDh25GR6NjMFSUVGsT+eARpeUVSxyIiMhssOkT1zEohx4IxbRH2113PRy08giu3SqWORURkFlh0iBqAnaUFEp6NRnMPO+SoyjH6hyO4XqSWOhYRkclj0SFqIC62SiwZ3x6NnKyRcaMEY348gsKySqljERGZNBYdogbk5WiFZRPaw83OEinZKryw5BjUVRqpYxERmSwWHaIG1sTNFj+Nawc7SwscSs/Hy78kQ6s12+msiIjqFYsOkQSCfRwxf3QUFHIBm05l472NZ2HGc3cSEdUbFh0iiXQKdMNnT4cDABYduIT5e9KlDUREZIJYdIgk9ESYD97q3woA8NHv57A68YrEiYiITAuLDpHEJnRuiuc6BwAAXl99EvvTbkiciIjIdLDoEBmA+L6t8HioN6q0IiYuTcT5XM6eTESkDyw6RAZAJhPw6dAwtPV3RlF5FZ5NOIq8onKpYxERGT0WHSID8fetIgLcbHG1oAzjFx1DaUWV1LGIiIwaiw6RAXGxVSJhbDs42yhw6mohpq04AQ3n2CEiqjMWHSID08TNFgtj20JpIcOOlDzM3pQidSQiIqPFokNkgKL8XfDFX3Ps/Lg/A8sOZ0obiIjISLHoEBmo/qHeeKV3CwDAO+vP8LJzIqI6YNEhMmAv9gjEoHAfaLQiJi1NxMXrxVJHIiIyKiw6RAZMEAR8NDgUUf7OUJVXYfyio7hVUiF1LCIio8GiQ2TgrBRyzB8dBV9na1y6WYpJyxJRqdFKHYuIyCiw6BAZATc7S/wQ2w52lhY4lJ6P9347K3UkIiKjwKJDZCRaetlj7rBwCAKw5FAmlh7ilVhERA/CokNkRHq19sSMPi0BALM2nMGh9JsSJyIiMmwsOkRGZlLXZngizAdVf12JdTm/VOpIREQGi0WHyMgIgoBPhoQipJEjbpVW4rnFx1Ci5j2xiIjuhkWHyAjdvgFoFNztLXEupwivrEyGKPKeWERE/8aiQ2SkvB2t8d2oKCjkAracycE3u9KkjkREZHBYdIiMWJS/M94f2AYA8Nn28/jjXK7EiYiIDAuLDpGRGx7dGCPbN4YoAi+tSEI6bxNBRKTDokNkAmYOCEZbf2cUqavw/JJEFJVXSh2JiMggsOgQmQClhQzfjoqEl4MV0vKK8fLKZGi1PDmZiIhFh8hEeNhb4bvRUVDKZdh+Nhff7ubJyURELDpEJiTczwnvDwoGcPvk5N2peRInIiKSFosOkYkZ1q4xRkT/dXLyz0nIusmZk4nIfLHoEJmgWU+0RpifEwrLKvHC0kSUVWikjkREJAkWHSITZGkhx3ejIuFqq0RKtgrxa05y5mQiMkssOkQmytvRGl8/Ewm5TMC6pGv46cAlqSMRETU4Fh0iExbTzBXxfYMAAB9sSkFiZr7EiYiIGhaLDpGJG/9IAPqHeqNKK2LysuO4XqSWOhIRUYNh0SEycYIg4OPBoQj0sEOuSo2pK46jSqOVOhYRUYNg0SEyA3aWFvhuVBRslXIcSs/Hf7elSh2JiKhBsOgQmYlADzv8d2gYAGD+n+nYcjpb4kRERPWPRYfIjPQL8caERwIAAK+uOsk7nRORyWPRITIzr/cNQnSAC4rVVZi09DhKK6qkjkREVG9YdIjMjEIuw9cjIuBub4nU3CK8ufY0JxMkIpNlEEVHo9Hg7bffRkBAAKytrdGsWTO8//779/3hu3v3bgiCcMcjJyenAZMTGScPByt8PSICcpmAtSeuYtnhLKkjERHVCwupAwDAxx9/jHnz5uGnn35CcHAwjh07hmeffRaOjo6YNm3afbdNTU2Fg4OD7rmHh0d9xyUyCe2buuL1x1riw83n8N5vZxHSyBFhfk5SxyIi0iuDKDoHDhzAwIED0b9/fwBAkyZNsGLFChw5cuSB23p4eMDJyameExKZpuc6N8XxzAJsOZODycuOY+PUR+Bsq5Q6FhGR3hjEoauOHTti586dOH/+PAAgOTkZ+/btQ9++fR+4bXh4OLy9vdG7d2/s37//vuuq1WqoVKpqDyJzJggCPhkaiiauNrhaUIa4X5Kg1fJ8HSIyHQZRdN544w0MHz4cQUFBUCgUiIiIQFxcHEaOHHnPbby9vfHdd99h9erVWL16Nfz8/NCtWzccP378ntvMmTMHjo6Ouoefn199fBwio+JgpcC8UVGwUsjw5/nr+OqPNKkjERHpjSAawOUWP//8M2bMmIH//ve/CA4ORlJSEuLi4vD5558jNja2xvvp2rUrGjdujCVLltz1dbVaDbX6f/f5UalU8PPzQ2FhYbXzfIjM0erEK3hlVTIEAfjp2Wh0aeEudSQiortSqVRwdHSs0e9vgxjRmTFjhm5UJyQkBKNHj8b06dMxZ86cWu0nOjoaaWn3/r9RS0tLODg4VHsQ0W2Do3wxIroxRBF46ecTuFZQJnUkIqKHZhBFp7S0FDJZ9ShyuRxabe1uPJiUlARvb299RiMyKzMHtEabRg64VVqJycuOo6KKN/8kIuNmEEVnwIABmD17NjZt2oRLly5h7dq1+Pzzz/Hkk0/q1omPj8eYMWN0z+fOnYv169cjLS0Np0+fRlxcHP744w9MmTJFio9AZBKsFHLMGxkFR2sFki4XYPams1JHIiJ6KAZxeflXX32Ft99+G5MnT0ZeXh58fHzwwgsv4J133tGtk52djays/01qVlFRgVdeeQVXr16FjY0NQkNDsWPHDnTv3l2Kj0BkMvxcbPDFsDCMW3QMPx3MRFQTFzwR5iN1LCKiOjGIk5GlUpuTmYjMzX+3nsM3uy7CRinH+imd0NzTXupIREQAjPBkZCIyPC/3bolOga4ordBg4tJEFKt5808iMj4sOkR0V3KZgP8bHgEvBytcvF6CN1af5M0/icjosOgQ0T252Vnim5ERsJAJ2HgyG4sOXJI6EhFRrbDoENF9Rfm74D/9WgEAZm9KQWLmLYkTERHVHIsOET3Qs52aoH+oN6q0IqYsO44bxeoHb0REZABYdIjogQRBwMeDQ9HM3RY5qnJMXX4CVRpOJkhEho9Fh4hqxM7SAvNHR8FWKcfB9Jv477ZUqSMRET0Qiw4R1Vighz0+GRIGAJj/Zzq2nM6ROBER0f2x6BBRrfQP9caERwIAAK+uSkb69WKJExER3RuLDhHV2ut9gxDdxAXF6ipMXJqIEk4mSEQGqlb3utqwYUOt36B3796wtrau9XZEZLgUchm+fiYC/b/ah/O5xXhjzSl8OTwcgiBIHY2IqJpa3etKJqvdAJAgCLhw4QKaNm1a62ANgfe6Ino4Ry/lY8SCQ6jSinirfytM6GyY/60TkWmp13td5eTkQKvV1uhhY2NT5w9BRIavXRMXvP14awDAnN/P4cDFGxInIiKqrlZFJzY2tlaHoUaNGsWREiITNybGH09FNIJGK2Lq8hO4VlAmdSQiIp1aHboyNTx0RaQfZRUaDJ53AGezVQjzdcQvL8TASiGXOhYRmah6PXRFRPRv1ko55o+OgpONAslXCjFrwxne6ZyIDEKdi05ZWRlKS0t1zzMzMzF37lxs27ZNL8GIyLj4udjgy+ERkAnAz0cvY+nhLKkjERHVvegMHDgQixcvBgAUFBSgffv2+OyzzzBw4EDMmzdPbwGJyHh0aeGO1x4LAgC8u+EMjmTkS5yIiMxdnYvO8ePH0blzZwDAr7/+Ck9PT2RmZmLx4sX48ssv9RaQiIzLC12aYkCYD6q0IiYvS+TJyUQkqToXndLSUtjb2wMAtm3bhqeeegoymQwdOnRAZmam3gISkXERBAGfDA5Fa28H3CiuwPNLjqG8UiN1LCIyU3UuOoGBgVi3bh0uX76MrVu34tFHHwUA5OXl8QomIjNnrZRjwZgouNgqcfqqCm+sPsmTk4lIEnUuOu+88w5effVVNGnSBNHR0YiJiQFwe3QnIiJCbwGJyDj5Otvgm2ciIZcJWJd0DQv2pEsdiYjMUK2LzjvvvIPExEQMGTIEWVlZOHbsGLZu3ap7vWfPnvjiiy/0GpKIjFNMM1e889fMyR9tOYedKbkSJyIic1PronPlyhX07dsXvr6+ePfdd5GTk4Oqqv/duTg6OhpBQUF6DUlExmtMjD+ead8YoghMW3ECqTlFUkciIjNS66Lz448/IicnBytWrIC9vT3i4uLg5uaGwYMHY/HixcjP5+WkRPQ/giDg3SeC0aGpC0oqNBj/01HcLFZLHYuIzIRebgGRkpKC3377DevXr0diYiKio6PxxBNPYMSIEWjUqJE+ctYL3gKCqOHcKqnAwG/2Iyu/FNEBLlg6vj2UFpycnYhqr8FvAdGqVSu89tpr2L9/Py5fvozY2Fjs3bsXK1as0MfuicgEONsq8UNsW9hbWuBIRj7eXneaV2IRUb3jTT05okPUoHal5mH8oqPQisAbfYMwsWszqSMRkZGpze9vi4d5o/Lycpw8eRJ5eXnQarXVXnviiSceZtdEZKK6t/TAW/1b472NZ/HR7+fg52yD/qHeUsciIhNV56KzZcsWjBkzBjdu3LjjNUEQoNFwJlQiurtxjwQgK78Uiw5cwvSVSfBytEKUv7PUsYjIBNX5HJ2pU6di6NChyM7OhlarrfZgySGiB3n78dbo1coDFVVaPL/4GLJulkodiYhMUJ2LTm5uLl5++WV4enrqMw8RmQm5TMD/DY9Am0YOuFlSgbGLjqCwtFLqWERkYupcdIYMGYLdu3frMQoRmRtbSwv8GNsOPo5WSL9egucW8wagRKRfdb7qqrS0FEOHDoW7uztCQkKgUCiqvT5t2jS9BKxPvOqKyDCcy1Fh6HcHUVRehT7Bnvh2ZBTkMkHqWERkoGrz+7vOReeHH37AxIkTYWVlBVdXVwjC/34oCYKA9HTDv4Efiw6R4TiUfhNjfjyCiiotnmnfGLMHtan2c4WI6G8NUnS8vLwwbdo0vPHGG5DJjHN2UxYdIsPy+6lsTF5+HKIITO/VAi/1ai51JCIyQA0yM3JFRQWGDRtmtCWHiAxP3xBvvPdEMADgix3nsfxwlsSJiMjY1bmlxMbG4pdfftFnFiIijI5pgqk9AgEAb607hU0nsyVORETGrM4TBmo0GnzyySfYunUrQkND7zgZ+fPPP3/ocERknl7u3QI3itVYceQyXvr5BKyVMvQI4lQWRFR7dS46p06dQkREBADg9OnTegtERCQIAj4YFIIStQYbkq9h4tLjWDS2HToGukkdjYiMDG/qyZORiQxWpUaLycuOY/vZXNgo5Vgyvj1vFUFEDXMy8ooVK+752owZM+q6WyIiHYVchq9GRKBzczeUVmgwNuEITl8tlDoWERmROhedSZMm4ffff79j+fTp07F06dKHCkVE9DcrhRzzR0ehXRNnFJVXYdQPh1l2iKjG6lx0li1bhhEjRmDfvn26ZVOnTsXKlSuxa9cuvYQjIgIAG6UFfhjbDuF+TigorcTIhYdx6grLDhE9WJ2LTv/+/fHtt9/iiSeeQGJiIiZPnow1a9Zg165dCAoK0mdGIiI4WCmwZHw0Ihs7obCsEs8sPITkywVSxyIiA1fnq64A4JlnnkFBQQE6deoEd3d3/PnnnwgMDNRXNiKiauytFFg8vj3G/ngExzJvYdTCw1g8PhoRjXmCMhHdXa2uunr55ZfvunzVqlWIjIxEs2bNdMtqM4+ORqPBrFmzsHTpUuTk5MDHxwdjx47FW2+9dd973ezevRsvv/wyzpw5Az8/P7z11lsYO3Zsjd+XV10RGadidRXGJRzFkUv5sLO0wA+xbdG+qavUsYiogdTm93etRnROnDhx1+WBgYFQqVS612t7I76PP/4Y8+bNw08//YTg4GAcO3YMzz77LBwdHe95F/SMjAz0798fEydOxLJly7Bz505MmDAB3t7e6NOnT63en4iMi52lBRaNa4dxi47iUHo+xvx4BF8/E4nerTmpIBFVZxDz6Dz++OPw9PTEDz/8oFs2ePBgWFtb3/MKrtdffx2bNm2qNlnh8OHDUVBQgC1bttTofTmiQ2Tcyis1eHH5CexIyYVcJuCjp0IwtK2f1LGIqJ7V2zw6J0+ehFarrfH6Z86cQVVV1QPX69ixI3bu3Inz588DAJKTk7Fv3z707dv3ntscPHgQvXr1qrasT58+OHjw4D23UavVUKlU1R5EZLysFHJ8NyoSQ6J8odGKmPHrScz/86LUsYjIgNSq6ERERODmzZs1Xj8mJgZZWQ+++/Abb7yB4cOHIygoCAqFAhEREYiLi8PIkSPvuU1OTg48PasPU3t6ekKlUqGsrOyu28yZMweOjo66h58f/8+PyNhZyGX475BQvNClKQBgzu/nMHvTWWi0kg9WE5EBqNU5OqIo4u2334aNjU2N1q+oqKjReitXrsSyZcuwfPlyBAcHIykpCXFxcfDx8UFsbGxtIt5XfHx8tROqVSoVyw6RCRAEAfH9WsHFVok5v5/D93szcOlmKeYOC4et5UNdXEpERq5WPwG6dOmC1NTUGq8fExMDa2vrB643Y8YM3agOAISEhCAzMxNz5sy5Z9Hx8vJCbm5utWW5ublwcHC453taWlrC0tKyxvmJyLi80LUZvBytMOPXk9h+NhdDvzuIH8a2hbfjg38OEZFpqlXR2b17d72EKC0thUxW/SiaXC6/7/lAMTEx2Lx5c7Vl27dvR0xMTL1kJCLjMDC8EXydbfDCkmM4m63CwK/34/sxbRHm5yR1NCKSQJ1nRtanAQMGYPbs2di0aRMuXbqEtWvX4vPPP8eTTz6pWyc+Ph5jxozRPZ84cSLS09Px2muv4dy5c/j222+xcuVKTJ8+XYqPQEQGJMrfGeumdEJLT3vkFanx9PyDWHP8itSxiEgCBnF5eVFREd5++22sXbsWeXl58PHxwYgRI/DOO+9AqVQCAMaOHYtLly5VG1XavXs3pk+fjrNnz8LX1xdvv/02JwwkIp2i8kpMW3ECu1KvAwBGRDfGzAGtYaWQS5yMiB5GbX5/G0TRkQqLDpHp02hFfPXHBfzfzgsQRSDYxwHzRkahsWvNLqogIsNTb/PoEBEZG7lMQFyvFvjp2Wg42yhw5poK/b/ai61ncqSORkQNgEWHiMxClxbu2DStMyIbO6GovAovLEnEjFXJKCqvlDoaEdUjvRWdESNG6GYa3rBhA1atWqWvXRMR6YWPkzV+fj4GL3RpCkEAViVewWNz9+LAxRtSRyOieqK3onPmzBk4ODjg7Nmz+M9//oPdu3dj6tSp+to9EZFeKC1kiO/XCr88HwM/F2tcLSjDM98fxru/nUFZhUbqeESkZ3orOgqFAqIoIiEhAfHx8fjmm29w4MABfe2eiEivogNc8PtLXTAiujEAIGH/JfT6/E9sPZMDM75Gg8jk6K3oTJo0CZGRkfj1118xaNAgAEBJSYm+dk9EpHd2lhaY81QIEp5th0ZOt0d3XliSiHGLjiLzJn9+EZkCvV5eXlBQAAsLC9jZ2SEtLQ0ffPABFi1apK/d6x0vLyeiv5VWVOGbXWlYsCcdlRoRSgsZXujSFM93aQp7K4XU8YjoHySZR6eyshLLli3D9evX0bp1a/Tt2/eO2zoYGhYdIvq3i9eLMXP9GexLu32CsrONApO6NcOYmCacaJDIQEgyj87w4cNx7NgxWFtbY+PGjYiMjMT58+f1tXsiogbRzN0OS8ZH47tRkWjqbotbpZX4cPM5dP3vLiw9lImKqnvfg4+IDI/eRnQiIiJw4sQJ3fOkpCRMmzYNe/bs0cfu6wVHdIjofqo0Wqw9cRVzd1zA1YIyAICHvSViOzbBM9GN4WyrlDghkXmSZETH3t4eaWlpuufh4eG4deuWvnZPRNTgLOQyDG3rhz9e7Yp3nwiGp4Ml8orU+O/WVMR8tBNvrj2FtLwiqWMS0X3obUTn5MmTGDFiBPr164fWrVsjJSUFZ8+excaNG/Wx+3rBER0iqo2KKi02nbqGhXszcOaaSrc8orETnor0xYBQbzjZcJSHqL5JdlNPtVqNdevWISUlBd7e3hg9ejRsbAz3xnksOkRUF6Io4nBGPhbuzcCu1DxotLd/jCrlMvQI8kD/UG90aeEOR2terUVUHxq06MTFxSEsLAyhoaFo06YNLC0tH2Z3DYpFh4geVl5ROTYkXcPq41eRkv2/UR65TEBbf2f0bOWB7i09EOhhB0EQJExKZDoapOgUFxfDzs4OGzZswMmTJ3Hy5EmcOXMGgiAgODgYoaGhCA0NxYABA+r0IRoCiw4R6dPZayqsT76KP1LycCGvuNprzjYKRPk7I9LfGW39XdCmkQNslBYSJSUybg1SdORyOVauXInBgwdXW15eXo7Tp0/j5MmTOHXqFL744ou67L5BsOgQUX3JulmKP87lYue5PBzJyIf6Lpel+7lYo4WHPVp42aO5hx18nW3g42QFTwcrKOSGPQ8ZkZQapOjIZDL06tULJSUlEAQBbdu2xciRI9GuXbs6hZYCiw4RNYSKKi3OXCtEYuYt3SOvSH3P9WUC4OlgBXd7SzhaK+Bso4STjQJO1gpYKeWwVshhpZDDSiGDQi6DXBAgCALkMgEyAdCKgFYUodWK0IpAlVYLjVZElUZElVZElVaLSo0IzV//rNKI0OjWv/3voggIAATh9j5v71uApUIGK4vb729pIYOdlQUcrRW6h7ONEtZKTqxI9avBio6rqyuGDx8OW1tbJCYmYu/evXjxxRfx6aef1il4Q2PRISKp5JdU4Hxuke5xMa8E1wrLkF1QjgqNcU9KaG9lAS8HK3g5WsHLwQqNnK3RzN0OTd1t0dTNjkWIHlptfn8/1AHi5cuXo3fv3rrnJ0+exMCBA9GoUSNMnz79YXZNRGTSXGyV6NDUFR2aulZbrtWKuFGixrWCctwsVqOgtBK3SitQWFaJwrJKlFVoUFapQXmlFuWVGlRqtBBF3B6REW+P4MgEQPbXSMztfwqwkAtQyGWwkN3+dwuZ7K9/CrD4a1RILhMgCPhrhOivPH+NDol/jQypK7Uor7r93uWVGhSrq1BYVgnVX/kqNSKKyqtQVF58x3lKf2vkZI3WPg4I83VEqK8TQn0deVk+1Zs6j+i4ublh3759CAoKqrZ806ZNmD59ulHc/oEjOkRE+iOKIorVVchVqZFTWI4cVTlyCsuQebMU6TdKcPF6MQpKK++6bVM3W3QMdMUjgW6IaeoGRxtemk/31iCHrnr16oWoqCh8/PHH1ZanpqYiLCwM5eXlddltg2LRISJqWPklFbiQW4RTVwuRfKUQJ68UIPNmabV1BAEIaeSIHkEe6Bfijea8NJ/+pUGKzqFDh9C9e3cMGTIEkydPRmhoKEpKSjBjxgwcOXIEKSkpdQrfkFh0iIikV1BagaOXbmF/2g3sS7uBtH8d8mrmbou+bbzRP9Qbrbz5s5oacMLA5ORkvPTSS9i3bx/+3o2VlRVWrVqFfv361XW3DYZFh4jI8OQUlmPPhevYejoHey/cqHZydptGDhjerjEGhvvA3oqHt8xVg98CIi8vD4mJidBqtWjfvj3c3NwedpcNgkWHiMiwqcorsetcHjadzMbu1Ou60mOtkOPxUG+MiWmCEF9HiVNSQ5PsXlfGhkWHiMh45JdUYM3xK/j56OVqh7c6NnPFxK7N0Lm5G8/lMRMsOjXEokNEZHxEUURi5i0sO5yF35Kvoeqvm6q28nbAxK5N8XioD+QyFh5TxqJTQyw6RETG7WpBGX7cl4EVR7JQWqEBALT0tMeMPi3Rs5UHR3hMFItODbHoEBGZhoLSCiw5mImF+zJQWHZ7rp62/s54o28Q2jZxkTgd6RuLTg2x6BARmZbCskp89+dFJOzPQHnl7ROXH23tiXcGtIavs43E6UhfWHRqiEWHiMg05arKMXfHBaw8dhkarQhrhRzTejbH+EcCoLTgneGNHYtODbHoEBGZtvO5RXh73WkczsgHAAR62OH9gW0Q08z1AVuSIavN72/WWiIiMlktPO3x8/Md8PnTYXCzUyItrxgjvj+E+DUnUayukjoeNQAWHSIiMmmCIOCpSF/sfLkbRnVoDEEAVhy5jL7/twdH/hrpIdPFokNERGbB0UaBDwaFYPmEDmjkZI3L+WUYtuAgPtycgvJKjdTxqJ6w6BARkVmJaeaKLXGd8XRbX4gisGBPOgZ+vf+Om4mSaWDRISIis2NvpcAnQ8KwcExbuNkpkZpbhIFf78OG5GtSRyM9Y9EhIiKz1au1Jza/1BkdmrqgpEKDaStO4J31p6Gu4qEsU8GiQ0REZs3D3gpLx7fHlO7NAACLD2bi6e8O4sqtUomTkT6w6BARkdmzkMswo08QEsa2g5ONAslXCjHom/1IzLwldTR6SCw6REREf+ke5IFN0zqjtbcDbhRXYMT3h7DuxFWpY9FDYNEhIiL6h0ZO1lg1MQaPtvZERZUWcb8k4bNtqdBqzfZGAkaNRYeIiOhfbC0t8N2oKEzsevu8na/+SMPUFSc4344RYtEhIiK6C5lMwBt9g/DfIaFQyAVsOpWNsQlHUFReKXU0qgUWHSIiovsY2tYPi8e1h52lBQ6l52P4gkO4XqSWOhbVEIsOERHRA8Q0c8XPz3eAm50SZ66pMPS7A7icz8vPjQGLDhERUQ20aeSIVRM7wtfZGpdulmLwvAM4l6OSOhY9gMEUnSZNmkAQhDseU6ZMuev6ixYtumNdKyurBk5NRETmJMDNFqsndURLT3vkFakxYsEhnL3GsmPIDKboHD16FNnZ2brH9u3bAQBDhw695zYODg7VtsnMzGyouEREZKY8Hayw8oUYhPk54VZpJUYuPIQz1wqljkX3YDBFx93dHV5eXrrHxo0b0axZM3Tt2vWe2wiCUG0bT0/PBkxMRETmytFGgSXjo/9Rdg7j9FWWHUNkMEXnnyoqKrB06VKMGzcOgiDcc73i4mL4+/vDz88PAwcOxJkzZ+67X7VaDZVKVe1BRERUFw5Wt8tOuJ8TClh2DJZBFp1169ahoKAAY8eOvec6LVu2xI8//oj169dj6dKl0Gq16NixI65cuXLPbebMmQNHR0fdw8/Prx7SExGRufi77EQ0dkJh2e2yw3N2DIsgiqLBzWndp08fKJVK/PbbbzXeprKyEq1atcKIESPw/vvv33UdtVoNtfp/cx+oVCr4+fmhsLAQDg4OD52biIjMU1F5JWJ/PILjWQVws7PEqokxCHCzlTqWyVKpVHB0dKzR72+DG9HJzMzEjh07MGHChFptp1AoEBERgbS0tHuuY2lpCQcHh2oPIiKih2VvpUDCs9F/3QxUjVELDyO7sEzqWAQDLDoJCQnw8PBA//79a7WdRqPBqVOn4O3tXU/JiIiI7s3RWoGfxkUjwM0WVwvKMGrhYdws5gzKUjOooqPVapGQkIDY2FhYWFhUe23MmDGIj4/XPX/vvfewbds2pKen4/jx4xg1ahQyMzNrPRJERESkL+72llg6oT18HK1w8XoJxiYc5b2xJGZQRWfHjh3IysrCuHHj7ngtKysL2dnZuue3bt3Cc889h1atWqFfv35QqVQ4cOAAWrdu3ZCRiYiIqmnkZI0lE9rD1VaJU1cL8dziY1BX8a7nUjHIk5EbSm1OZiIiIqqN01cLMXzBIRSrqzAw3Adzh4Xfd8oUqjmjPhmZiIjIFLRp5Ih5oyJhIROwPukaPt2WKnUks8SiQ0REVE86N3fHh0+FAAC+2XURK45kSZzI/LDoEBER1aOn2/phWs/mAIC31p3GrtQ8iROZFxYdIiKieja9V3M8FdkIGq2IKcuOc/bkBsSiQ0REVM8EQcBHT4WiU6ArSis0eG7xMdzgHDsNgkWHiIioASgtZPj2mSg0cbXB1YIyTFqaiIoqrdSxTB6LDhERUQNxtFFgYWw72Fta4OilW3hn/WmY8SwvDYJFh4iIqAEFetjhy2ciIAjAz0cvY/HBTKkjmTQWHSIiogbWvaUH4vsGAQDe23gW+9NuSJzIdLHoEBERSeC5zk3xVMRfV2ItP44rt0qljmSSWHSIiIgkIAgCPnwqBKG+jigorcTkZcdRXsl7Yukbiw4REZFErBRyfDsyEk42Cpy8Uoj3Np6VOpLJYdEhIiKSkK+zDf5v+O2Tk5cfzsKviVekjmRSWHSIiIgk1rWFO+J6tgAAvLn2FGdO1iMWHSIiIgMwtUcgurV0h7pKi0nLElFYVil1JJPAokNERGQAZDIBc4eFw9fZGpk3S/H6ryc5maAesOgQEREZCCcbJb4dGQmFXMCWMzlYejhL6khGj0WHiIjIgIT6OuH1x25PJvj+xrM8X+chsegQEREZmPGPBKBHkAcqqrR4ccVxlFZUSR3JaLHoEBERGRhBEPDp0DB4Olgi/XoJ3ll/RupIRotFh4iIyAC52Crxf8MjIBOAXxOvYO0Jzq9TFyw6REREBqpDU1dM69kcAPDW2tPIusn7YdUWiw4REZEBm9qjOaKbuKCkQoPpK5NQpdFKHcmosOgQEREZMLlMwGdPh8He0gKJmbcwb/dFqSMZFRYdIiIiA+fnYoN3BwYDAP5v5wUkXy6QNpARYdEhIiIyAk9GNEL/UG9UaUVM/yWJl5zXEIsOERGRERAEAbMHtYGXgxXSb5Tgw80pUkcyCiw6RERERsLJRonPng4DACw9lIU/zuVKnMjwsegQEREZkU6Bbhj/SAAA4I3Vp1BYyruc3w+LDhERkZGZ0aclmrrbIq9IjXd/46zJ98OiQ0REZGSsFHJ8OjQMMgFYc+Iqtp3JkTqSwWLRISIiMkKRjZ3xfJdmAID/rD2NWyUVEicyTCw6RERERiquV3M097DDjWI1Zm7gIay7YdEhIiIyUn8fwpLLBGxIvoYtp7OljmRwWHSIiIiMWJifEyZ2bQoAeHPtaeTzEFY1LDpERERGblrP5mjpaY+bJRX4YONZqeMYFBYdIiIiI2dpIcfHQ0J1V2HtTs2TOpLBYNEhIiIyAeF+Tni20+2JBN9cexrFat4LC2DRISIiMhmvPNoCvs7WuFpQhk+3pkodxyCw6BAREZkIG6UFPnwyBADw08FLSMy8JXEi6bHoEBERmZAuLdwxONIXogi8vvok1FUaqSNJikWHiIjIxLz9eCu42SmRlleMb3ZdlDqOpFh0iIiITIyTjRKznggGAHy3+yIuXi+WOJF0WHSIiIhMUP8Qb3Rr6Y4KjRZvrj0FURSljiQJFh0iIiITJAgC3h/YBlYKGQ6l52PN8atSR5IEiw4REZGJ8nOxwbSezQEAszenmOUdzg2m6DRp0gSCINzxmDJlyj23WbVqFYKCgmBlZYWQkBBs3ry5ARMTEREZvuc6N0ULTzvkl1Tgo9/PSR2nwRlM0Tl69Ciys7N1j+3btwMAhg4detf1Dxw4gBEjRmD8+PE4ceIEBg0ahEGDBuH06dMNGZuIiMigKeQy3dw6vxy7jKOX8iVO1LAE0UDPToqLi8PGjRtx4cIFCIJwx+vDhg1DSUkJNm7cqFvWoUMHhIeH47vvvqvRe6hUKjg6OqKwsBAODg56y05ERGRo3lh9Ej8fvYwWnnbYOLUzlBYGM9ZRa7X5/W2Qn7KiogJLly7FuHHj7lpyAODgwYPo1atXtWV9+vTBwYMH77lftVoNlUpV7UFERGQO3ugbBFdbJc7nFiNhf4bUcRqMQRaddevWoaCgAGPHjr3nOjk5OfD09Ky2zNPTEzk5OffcZs6cOXB0dNQ9/Pz89BWZiIjIoDnZKPFG3yAAwP/tvIDswjKJEzUMgyw6P/zwA/r27QsfHx+97jc+Ph6FhYW6x+XLl/W6fyIiIkM2ONIXUf7OKK3QYPamFKnjNAiDKzqZmZnYsWMHJkyYcN/1vLy8kJubW21Zbm4uvLy87rmNpaUlHBwcqj2IiIjMhUwm4L2BwZAJwMaT2TiQdkPqSPXO4IpOQkICPDw80L9///uuFxMTg507d1Zbtn37dsTExNRnPCIiIqMW7OOI0R38AQDvbDiDiiqtxInql0EVHa1Wi4SEBMTGxsLCwqLaa2PGjEF8fLzu+UsvvYQtW7bgs88+w7lz5zBr1iwcO3YML774YkPHJiIiMiovP9oSrra3b/pp6icmG1TR2bFjB7KysjBu3Lg7XsvKykJ2drbueceOHbF8+XIsWLAAYWFh+PXXX7Fu3Tq0adOmISMTEREZHUdrhdmcmGyw8+g0BM6jQ0RE5kqrFTF0/kEkZt7C46He+PqZSKkj1ZjRz6NDRERE9evfJyYfTr8pdaR6waJDRERkpoJ9HDE8ujEAYNZvZ6HRmt5BHhYdIiIiM/bqoy3hYGWBlGwVfj6aJXUcvWPRISIiMmMutkpM790CAPDp1lQUllZKnEi/WHSIiIjM3KgO/mjuYYdbpZWYu/O81HH0ikWHiIjIzCnkMswcEAwAWHwwExdyiyROpD8sOkRERIRHmruhd2tPaLQi3tt4FqYy+wyLDhEREQEA3urfCkq5DHsv3MCOlDyp4+gFiw4REREBAPxdbTG+cwAA4MPNKSZxHywWHSIiItKZ3K0Z3OyUyLhRgiWHMqWO89BYdIiIiEjH3kqBVx5tCQD4cucFFJRWSJzo4bDoEBERUTVPt/VDkJc9CssqMXfHBanjPBQWHSIiIqpGLhPwVv/WAIClhzJx8XqxxInqjkWHiIiI7vBIczf0DPJAlVbEnM0pUsepMxYdIiIiuqv4fq1gIROwIyUP+9NuSB2nTlh0iIiI6K4CPewwqoM/AOD9jcZ5d3MWHSIiIrqnl3o2h72VBc7lFGHN8StSx6k1Fh0iIiK6J2dbJV7sHggA+GzbeZRVaCROVDssOkRERHRfsR2boJGTNXJU5fhxf4bUcWqFRYeIiIjuy0ohx4w+tycRnLf7Im4UqyVOVHMsOkRERPRAT4T5oE0jBxSrq/DlTuOZRJBFh4iIiB5IJhPwn36tAADLD2ch3UgmEWTRISIiohrp2MwNPf6aRPDjLeekjlMjLDpERERUY/F9gyATgK1ncnH0Ur7UcR6IRYeIiIhqrLmnPYa18wMAfPT7OYiiYU8iyKJDREREtRLXqwWsFDIkZt7C9rO5Use5LxYdIiIiqhVPByuM6xQAAPhkayqqNFqJE90biw4RERHV2gtdm8HJRoG0vGKsOX5V6jj3xKJDREREteZorcCUbrdvDfH59vMorzTMW0Ow6BAREVGdjI7xh4+jFXJU5fjpwCWp49wViw4RERHViZVCjpcfvX1riG92paGwtFLiRHdi0SEiIqI6ezKiEVp62kNVXoVv/0yTOs4dWHSIiIiozuQyAa89dntUZ9H+S8gpLJc4UXUsOkRERPRQegR5oF0TZ6irtPjyD8O64SeLDhERET0UQRAwo08QAGDl0cu4dKNE4kT/w6JDREREDy06wAXdWrqjSiviix3npY6jw6JDREREevHqX1dgbUi+hpRslcRpbmPRISIiIr1o08gR/UO9IYrAp1tTpY4DgEWHiIiI9OiV3i0glwnYeS4PiZn5Usdh0SEiIiL9aepuhyGRvgCAT7akQhRFSfOw6BAREZFevdSrOZQWMhzOyMfeCzckzcKiQ0RERHrl42SN0R38AQAL9qRLmsVC0ncnIiIikzS5WzPYKuUY37mppDlYdIiIiEjvXO0sdTf8lBIPXREREZHJYtEhIiIik8WiQ0RERCbLYIrO1atXMWrUKLi6usLa2hohISE4duzYPdffvXs3BEG445GTk9OAqYmIiMiQGcTJyLdu3UKnTp3QvXt3/P7773B3d8eFCxfg7Oz8wG1TU1Ph4OCge+7h4VGfUYmIiMiIGETR+fjjj+Hn54eEhATdsoCAgBpt6+HhAScnp3pKRkRERMbMIA5dbdiwAW3btsXQoUPh4eGBiIgIfP/99zXaNjw8HN7e3ujduzf2799/33XVajVUKlW1BxEREZkugyg66enpmDdvHpo3b46tW7di0qRJmDZtGn766ad7buPt7Y3vvvsOq1evxurVq+Hn54du3brh+PHj99xmzpw5cHR01D38/Pzq4+MQERGRgRBEqe+2BUCpVKJt27Y4cOCAbtm0adNw9OhRHDx4sMb76dq1Kxo3bowlS5bc9XW1Wg21Wq17rlKp4Ofnh8LCwmrn+RAREZHhUqlUcHR0rNHvb4MY0fH29kbr1q2rLWvVqhWysrJqtZ/o6GikpaXd83VLS0s4ODhUexAREZHpMoii06lTJ6SmplZbdv78efj7+9dqP0lJSfD29tZnNCIiIjJiBnHV1fTp09GxY0d8+OGHePrpp3HkyBEsWLAACxYs0K0THx+Pq1evYvHixQCAuXPnIiAgAMHBwSgvL8fChQvxxx9/YNu2bVJ9DCIiIjIwBlF02rVrh7Vr1yI+Ph7vvfceAgICMHfuXIwcOVK3TnZ2drVDWRUVFXjllVdw9epV2NjYIDQ0FDt27ED37t2l+AhERERkgAziZGSpFBYWwsnJCZcvX+b5OkREREbi74uJCgoK4OjoeN91DWJERypFRUUAwMvMiYiIjFBRUdEDi45Zj+hotVpcu3YN9vb2EARBt7xdu3Y4evTofbd90Dp/t01THC2qyZ+Psb6/vvZd1/3UdrvarM/v9f1J/b2uzwz8XvN7bWoZRFFEUVERfHx8IJPd/7oqsx7Rkclk8PX1vWO5XC5/4Je9JusAMMnL2Gv62Y3x/fW177rup7bb1WZ9fq/vT+rvdX1m4Pea32tTzPCgkZy/GcTl5YZmypQpelnHVEn92evz/fW177rup7bb1WZ9fq/vzxA+e31l4Pda+r9bqRjCZ5c6g1kfuqpPtZm1kchY8HtNpojfa9PGEZ16YmlpiZkzZ8LS0lLqKER6w+81mSJ+r00bR3SIiIjIZHFEh4iIiEwWiw4RERGZLBYdIiIiMlksOkRERGSyWHSIiIjIZLHoSODJJ5+Es7MzhgwZInUUIr24fPkyunXrhtatWyM0NBSrVq2SOhLRQysoKEDbtm0RHh6ONm3a4Pvvv5c6EtUBLy+XwO7du1FUVISffvoJv/76q9RxiB5adnY2cnNzER4ejpycHERFReH8+fOwtbWVOhpRnWk0GqjVatjY2KCkpARt2rTBsWPH4OrqKnU0qgWO6EigW7dusLe3lzoGkd54e3sjPDwcAODl5QU3Nzfk5+dLG4roIcnlctjY2AAA1Go1RFEExwaMD4tOLe3ZswcDBgyAj48PBEHAunXr7ljnm2++QZMmTWBlZYX27dvjyJEjDR+UqBb0+b1OTEyERqOBn59fPacmuj99fK8LCgoQFhYGX19fzJgxA25ubg2UnvSFRaeWSkpKEBYWhm+++eaur//yyy94+eWXMXPmTBw/fhxhYWHo06cP8vLyGjgpUc3p63udn5+PMWPGYMGCBQ0Rm+i+9PG9dnJyQnJyMjIyMrB8+XLk5uY2VHzSF5HqDIC4du3aasuio6PFKVOm6J5rNBrRx8dHnDNnTrX1du3aJQ4ePLghYhLVSl2/1+Xl5WLnzp3FxYsXN1RUohp7mJ/Xf5s0aZK4atWq+oxJ9YAjOnpUUVGBxMRE9OrVS7dMJpOhV69eOHjwoITJiOquJt9rURQxduxY9OjRA6NHj5YqKlGN1eR7nZubi6KiIgBAYWEh9uzZg5YtW0qSl+qORUePbty4AY1GA09Pz2rLPT09kZOTo3veq1cvDB06FJs3b4avry9LEBm0mnyv9+/fj19++QXr1q1DeHg4wsPDcerUKSniEtVITb7XmZmZ6Ny5M8LCwtC5c2dMnToVISEhUsSlh2AhdQBztGPHDqkjEOnVI488Aq1WK3UMIr2Kjo5GUlKS1DHoIXFER4/c3Nwgl8vvOFktNzcXXl5eEqUiejj8XpMp4vfafLDo6JFSqURUVBR27typW6bVarFz507ExMRImIyo7vi9JlPE77X54KGrWiouLkZaWprueUZGBpKSkuDi4oLGjRvj5ZdfRmxsLNq2bYvo6GjMnTsXJSUlePbZZyVMTXR//F6TKeL3mgDw8vLa2rVrlwjgjkdsbKxuna+++kps3LixqFQqxejoaPHQoUPSBSaqAX6vyRTxe02iKIq81xURERGZLJ6jQ0RERCaLRYeIiIhMFosOERERmSwWHSIiIjJZLDpERERkslh0iIiIyGSx6BAREZHJYtEhIiIik8WiQ0R3NWvWLISHh0sd464aOluTJk0wd+5c3XNBELBu3bp6fT9BECAIAgoKCvS2327duun2y7tyk7lg0SEycmPHjtX98lIqlQgMDMR7772Hqqqqh9rvq6++Wu2Ghw/LkItTbWVnZ6Nv3771+h7vvfcesrOz4ejoqLd9rlmzBkeOHNHb/oiMAW/qSWQCHnvsMSQkJECtVmPz5s2YMmUKFAoF4uPj71i3oqICSqXygfu0s7ODnZ1dfcRtcKIoQqPRwMJCPz/yvLy89LKf+7G3t9f7+7i4uEClUul1n0SGjiM6RCbA0tISXl5e8Pf3x6RJk9CrVy9s2LABwO0Rn0GDBmH27Nnw8fFBy5YtAQCnTp1Cjx49YG1tDVdXVzz//PMoLi7W7fNuIzALFy5Eq1atYGVlhaCgIHz77bfVXr9y5QpGjBgBFxcX2Nraom3btjh8+DAWLVqEd999F8nJybrRp0WLFgEACgoKMGHCBLi7u8PBwQE9evRAcnJytf1+9NFH8PT0hL29PcaPH4/y8vL7/nns3r0bgiDg999/R1RUFCwtLbFv3z5cvHgRAwcOhKenJ+zs7NCuXTvs2LGj2rZ5eXkYMGAArK2tERAQgGXLlt2x/38euvr7vf55iCkpKQmCIODSpUsAgMzMTAwYMADOzs6wtbVFcHAwNm/efN/P8G+LFi2Ck5MT1q1bh+bNm8PKygp9+vTB5cuXq63322+/oV27drCysoKbmxuefPLJWr0PkanhiA6RCbK2tsbNmzd1z3fu3AkHBwds374dAFBSUoI+ffogJiYGR48eRV5eHiZMmIAXX3xRV0D+bdmyZXjnnXfw9ddfIyIiAidOnMBzzz0HW1tbxMbGori4GF27dkWjRo2wYcMGeHl54fjx49BqtRg2bBhOnz6NLVu26IrF34dkhg4dCmtra/z+++9wdHTE/Pnz0bNnT5w/fx4uLi5YuXIlZs2ahW+++QaPPPIIlixZgi+//BJNmzZ94J/DG2+8gU8//RRNmzaFs7MzLl++jH79+mH27NmwtLTE4sWLMWDAAKSmpqJx48YAbhfDa9euYdeuXVAoFJg2bRry8vIe5q8DU6ZMQUVFBfbs2QNbW1ucPXu2TqNlpaWlmD17NhYvXgylUonJkydj+PDh2L9/PwBg06ZNePLJJ/Hmm29i8eLFqKioqHWhIjI5Et89nYgeUmxsrDhw4EBRFEVRq9WK27dvFy0tLcVXX31V97qnp6eoVqt12yxYsEB0dnYWi4uLdcs2bdokymQyMScnRxRFUZw5c6YYFhame71Zs2bi8uXLq733+++/L8bExIiiKIrz588X7e3txZs3b94157/3J4qiuHfvXtHBwUEsLy+vtrxZs2bi/PnzRVEUxZiYGHHy5MnVXm/fvv0d+/qnXbt2iQDEdevW3XOdvwUHB4tfffWVKIqimJqaKgIQjxw5ons9JSVFBCB+8cUXumUAxLVr11Z7r1u3buleP3HihAhAzMjIEEVRFENCQsRZs2Y9MMvf/P39q72fKIpiQkKCCEA8dOjQHdkOHz4siuLtP6uRI0fed98ZGRkiAPHEiRM1zkNkzHjoisgEbNy4EXZ2drCyskLfvn0xbNgwzJo1S/d6SEhItfNyUlJSEBYWBltbW92yTp06QavVIjU19Y79l5SU4OLFixg/frzu3B07Ozt88MEHuHjxIoDbh2siIiLg4uJS49zJyckoLi6Gq6trtf1mZGTo9puSkoL27dtX2y4mJqZG+2/btm2158XFxXj11VfRqlUrODk5wc7ODikpKcjKytK9l4WFBaKionTbBAUFwcnJqcaf6W6mTZuGDz74AJ06dcLMmTNx8uTJOu3HwsIC7dq1uyNbSkoKgNt/Bz179nyorESmhoeuiExA9+7dMW/ePCiVSvj4+Nxx0u0/C01d/H3uzvfff39H6ZDL5QBuHy6ry369vb2xe/fuO1572HIB3Pm5X331VWzfvh2ffvopAgMDYW1tjSFDhqCioqLO7yGT3f7/RVEUdcsqKyurrTNhwgT06dMHmzZtwrZt2zBnzhx89tlnmDp1ap3f927q8ndAZOo4okNkAmxtbREYGIjGjRvX6MqiVq1aITk5GSUlJbpl+/fvh0wm052s/E+enp7w8fFBeno6AgMDqz0CAgIAAKGhoUhKSkJ+fv5d31OpVEKj0VRbFhkZiZycHFhYWNyxXzc3N13Ww4cPV9vu0KFDD/yMd7N//36MHTsWTz75JEJCQuDl5aU7YRi4PUJSVVWFxMRE3bLU1NT7zmXj7u4O4PYl53+72xw1fn5+mDhxItasWYNXXnkF33//fa3zV1VV4dixY3dka9WqFYDbfwf6nBKAyBSw6BCZoZEjR8LKygqxsbE4ffo0du3ahalTp2L06NHw9PS86zbvvvsu5syZgy+//BLnz5/HqVOnkJCQgM8//xwAMGLECHh5eWHQoEHYv38/0tPTsXr1ahw8eBDA7UnwMjIykJSUhBs3bkCtVqNXr16IiYnBoEGDsG3bNly6dAkHDhzAm2++qfuF/tJLL+HHH39EQkICzp8/j5kzZ+LMmTN1+tzNmzfHmjVrkJSUhOTkZDzzzDPQarW611u2bInHHnsML7zwAg4fPozExERMmDDhviMlgYGB8PPzw6xZs3DhwgVs2rQJn332WbV14uLisHXrVmRkZOD48ePYtWuXrpzUhkKhwNSpU3XZxo4diw4dOiA6OhoAMHPmTKxYsQIzZ85ESkoKTp06hY8//rjW70NkSlh0iMyQjY0Ntm7divz8fLRr1w5DhgxBz5498fXXX99zmwkTJmDhwoVISEhASEgIunbtikWLFulGdJRKJbZt2wYPDw/069cPISEh+Oijj3SHtgYPHozHHnsM3bt3h7u7O1asWAFBELB582Z06dIFzz77LFq0aIHhw4cjMzNTV7iGDRuGt99+G6+99hqioqKQmZmJSZMm1elzf/7553B2dkbHjh0xYMAA9OnTB5GRkdXWSUhIgI+PD7p27YqnnnoKzz//PDw8PO65T4VCgRUrVuDcuXMIDQ3Fxx9/jA8++KDaOhqNBlOmTEGrVq3w2GOPoUWLFndcml8TNjY2eP311/HMM8+gU6dOsLOzwy+//KJ7vVu3bli1ahU2bNiA8PBw9OjRgxMEktkTxH8eWCYi+kt8fDz27t2Lffv2SR3F7DRp0gRxcXGIi4vTLVu0aBHi4uIe+pYQly5dQkBAAE6cOGEyM1UT3Q9HdIioGlEUcfHiRezcuRPBwcFSxzFbr7/+Ouzs7FBYWKi3ffbt25d/p2R2eNUVEVVTWFiI1q1bo127dvjPf/4jdRyz9Oeff+qu3LK3t9fbfhcuXIiysjIA0E2QSGTqeOiKiIiITBYPXREREZHJYtEhIiIik8WiQ0RERCaLRYeIiIhMFosOERERmSwWHSIiIjJZLDpERERkslh0iIiIyGSx6BAREZHJ+n8RX7uPvZVWJwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "R_grid = np.geomspace(10., 3000., 100)\n",
    "sigma_grid = np.sqrt(model.sigmalos2(\n",
    "    R_grid, params=params, solver=\"kernel\", n_u=256, n_kernel=64))\n",
    "fig, ax = plt.subplots()\n",
    "ax.semilogx(R_grid, sigma_grid)\n",
    "ax.set(xlabel=\"Projected radius [pc]\", ylabel=r\"$\\sigma_{los}$ [km/s]\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1eb3f47",
   "metadata": {},
   "source": [
    "These host-side NumPy conversions are for plotting, after evaluation; keep\n",
    "JAX operations inside a differentiated likelihood. To check numerical\n",
    "stability, evaluate the same physical parameters and radii with larger\n",
    "`n_u` and `n_kernel` and compare the predictions. A finer radius grid makes\n",
    "the plot smoother but does not refine the underlying integration.\n",
    "\n",
    "## Inspect the components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d89464e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "tracer = model[\"StellarModel\"]\n",
    "halo = model[\"DMModel\"]\n",
    "surface_density = tracer.density_2d(R_grid, re_pc=params[\"re_pc\"])\n",
    "r_grid = np.geomspace(10., 3000., 100)\n",
    "spatial_density = tracer.density_3d(r_grid, re_pc=params[\"re_pc\"])\n",
    "halo_density = halo.mass_density_3d(r_grid, params=params)\n",
    "mass = halo.enclosed_mass(r_grid, params=params)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d6f8fa7",
   "metadata": {},
   "source": [
    "The tracer densities are number-weighting functions; the halo density and\n",
    "mass set the potential. Halo density and mass methods can apply different\n",
    "cutoff conventions. Read the selected halo's API entry before using its\n",
    "density in a custom integral. J/D factors additionally require a distance\n",
    "and an angular aperture; follow the [factor guide](../guides/factors.md).\n",
    "\n",
    "## Generate mock observations\n",
    "\n",
    "At a fixed observed radius, the Gaussian likelihood assumes\n",
    "\n",
    "$$\n",
    "v_i \\sim \\mathcal{N}\\!\\left(v_{\\rm mem},\n",
    "\\sigma_{\\rm los}^2(R_i)+\\epsilon_i^2\\right).\n",
    "$$\n",
    "\n",
    "In this expression the second argument is the **variance**. NumPy's\n",
    "`rng.normal` expects its square root as `scale`. The Quickstart implements\n",
    "this explicitly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "2412f1d3",
   "metadata": {
    "mystnb": {
     "image": {
      "alt": "Predicted LOS dispersion and 32 mock velocities with measurement errors."
     }
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA48AAAFpCAYAAAA86x25AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAcpdJREFUeJzt3XdUFNffBvBnll5XkKqgIKgooghYsGuw9xY1GnvsMWhMIkmsUdEUNWpiiYpRozFGoyb2XrGB2LsIiIAi0gWEnfcPX/cnAlLcZRb2+ZwzR2Z2ZvbZdeX63blzryCKoggiIiIiIiKid5BJHYCIiIiIiIg0H4tHIiIiIiIiKhSLRyIiIiIiIioUi0ciIiIiIiIqFItHIiIiIiIiKhSLRyIiIiIiIioUi0ciIiIiIiIqlK7UAaSkUCjw+PFjmJmZQRAEqeMQEVEpEkURKSkpqFSpEmQyfpfKNpGISHsVtU3U6uLx8ePHcHR0lDoGERFJKCoqCg4ODlLHkBzbRCIiKqxN1Ori0czMDMCrN8nc3FziNEREVJqSk5Ph6OiobAu0HdtEIiLtVdQ2UauLx9fdcszNzdlQEhFpKXbRfIVtIhERFdYm8iYPIiIiIiIiKhSLRyIiIiIiIioUi0ciIiIiIiIqFItHIiIiIiIiKhSLRyIiIiIiIioUi0ciIiIiIiIqFItHIiIiIiIiKhSLRyIiIiIiIiqUxhSPKSkp8Pf3R9WqVWFkZIQmTZrgwoULBe5/7NgxCIKQZ4mNjS3F1ERERERERNpBY4rHkSNH4uDBg9iwYQOuXr2Kdu3awc/PD9HR0e887vbt24iJiVEuNjY2pZT4lXtPUvDXhSiERCTgeVpWqT43ERGVf4GBgWjQoAHMzMxgY2ODHj164Pbt27n2ycjIwPjx41GxYkWYmpqid+/eiIuLkygxERGVV7pSBwCAFy9eYNu2bdi5cydatGgBAJg5cyb+/fdfLF++HHPmzCnwWBsbG1SoUKGUkuZ19NZTzN1zU7luYawHVxtTVLc1Qw0bU9SwNUNNOzNUNDWQLCMREZVdx48fx/jx49GgQQNkZ2fj66+/Rrt27XDjxg2YmJgAACZNmoTdu3dj69atkMvlmDBhAnr16oXTp09LnJ6IiMoTjSges7OzkZOTA0NDw1zbjYyMcOrUqXce6+npiczMTNSpUwczZ85E06ZNC9w3MzMTmZmZyvXk5OT3Cw7AxtwAzVyt8OBpKh4nZeB5+ktcePgcFx4+z72fmQFq2ZujdiVz1Kkkh0dlORwtjSAIwntnICKi8mvfvn251tetWwcbGxuEhISgRYsWSEpKwpo1a7Bp0ya0adMGABAUFIRatWrh7NmzaNy4sRSxiYhIzdKzslF7+n4AwI3Z7WGsr/7STiOKRzMzM/j6+uK7775DrVq1YGtri82bNyM4OBiurq75HmNvb48VK1bAx8cHmZmZWL16NVq1aoVz587By8sr32MCAwMxa9YslWbv7lkZ3T0rA3j1F/jgaRruPUnF3ScpuBOXijtxKYh4lo4nKZl4kvIUx+88VR4rN9JDXQc56jtWgGeVCvB0tIClib5K8xERUfmSlJQEALC0tAQAhISE4OXLl/Dz81Pu4+bmhipVqiA4OLjA4lEdX6gSEVH5phHFIwBs2LABw4cPR+XKlaGjowMvLy8MGDAAISEh+e5fs2ZN1KxZU7nepEkT3L9/H4sWLcKGDRvyPSYgIACTJ09WricnJ8PR0VFlr8FYXxd1KstRp7I81/a0zGzcik3BzZhkXH+cjGvRSbgdm4KkFy9x8m48Tt6NV+7rYm2CBk6W8HGyRCNnSzhaGqssHxERlW0KhQL+/v5o2rQp6tSpAwCIjY2Fvr5+nls4bG1t3zmInDq+UCUiovJNY4pHFxcXHD9+HGlpaUhOToa9vT369euHatWqFfkcDRs2fGc3VwMDAxgYlP69hyYGuvCuagHvqhbKbVnZCtyJS8GlqESERSbiUtRzPHiahvv/v/x5IQoAULmCEXxdKqKJS0U0c7WCjblhQU9DRETl3Pjx43Ht2rVCb+koCnV/oUpEROWPxhSPr5mYmMDExATPnz/H/v378f333xf52LCwMNjb26sxnero68qUVyk/blwVAPA8LQshEc9xISIBF8ITcOVREqITX+DvkEf4O+QRAMDNzgwtalijVQ1rNHC2hJ6OxgyYS0REajRhwgT8999/OHHiBBwcHJTb7ezskJWVhcTExFxXH+Pi4mBnZ1fg+aT6QpWIiMoujSke9+/fD1EUUbNmTdy7dw9ffPEF3NzcMGzYMACvviGNjo7G+vXrAQCLFy+Gs7Mz3N3dkZGRgdWrV+PIkSM4cOCAlC/jvViY6MOvti38atsCeNXd9WLEcwTff4bT9+Jx7XESbsWm4FZsCladeAAzA120qGmND9xs8IGbLeTGehK/AiIiUjVRFPHpp5/in3/+wbFjx+Ds7JzrcW9vb+jp6eHw4cPo3bs3gFfTWEVGRsLX11eKyEREVE5pTPGYlJSEgIAAPHr0CJaWlujduzfmzp0LPb1XBVFMTAwiIyOV+2dlZeHzzz9HdHQ0jI2NUbduXRw6dAitW7eW6iWonImBLlrWsEbLGtYAgGepmTh9/xlO3HmKY7efID41C7uvxGD3lRjoygT4ulREhzp2aO9uBytODUJEVC6MHz8emzZtws6dO2FmZqa8j1Eul8PIyAhyuRwjRozA5MmTYWlpCXNzc3z66afw9fXlSKtERKRSgiiKotQhpJKcnAy5XI6kpCSYm5tLHadYFAoRV6KTcPhmHA7eiMOt2BTlYzIBaOpqha71KqFDHTuYG/KKJBHR28pKG1DQlE5BQUEYOnQoACAjIwOff/45Nm/ejMzMTLRv3x6//vrrO7utvq2svB9ERPSKKqfqKGobwOKxnDSU4fFp2HctFvuuxeDyoyTldn1dGdrWskUfbwc0r24FXd4jSUQEoHy1AarA94OIqGyRonjUmG6r9H6crUwwtpULxrZyQcSzNPx7+TF2XX6MO3Gp2H01BruvxsDazAC9vCrjo4ZVULWiidSRiYiIiIhIBWKTMlDN2lTtz8Mrj+X8W9brj5Pwd8gj7Ax7jIS0LOX25tWt8FHDKmhb25ZXI4lIK2lDG1AcfD+IiMqWDcEPMW3ndQCvblsL7OWBfg2qlOhc7LZaBNrUUGZlK3Dk1hNsPh+JE3ef4vXfur3cEIN9nTCgoSMqGOtLG5KIqBRpUxtQFHw/iIjKjpikF2g6/wgUb1RyOoKAU1Nbw15uVOzzsdsq5aKvK0OHOnboUMcOUQnp2Hw+ElsuRCEmKQML9t3Cz4fvoLeXAz5pXg1OVuzSSkRERESkqcLj03IVjgCQI4p4GJ9eouKxqNhfUQs5Whrjyw5uOD21DX7oUxe17c2R8VKBP85Fos1PxzB+UyiuRScVfiIiIiIiIip1zlYmkL01GLeOIMDJylitz8viUYsZ6umgr48jdk9shj9HNUbrmtZQiMDuKzHosvQUhgadR1hUotQxiYiIiIjoDfZyI8zq5q5clwnAvF511HrVEWDxSHg1h1jjahURNKwh9n7WHD08K0FHJuDY7afo8ctpDF93AVceJUodk4iIiIiI/l9vbwflz4cmtyzxYDnFweKRcqllb47F/evj8OSW6OPtAB2ZgCO3nqDbstMYuzEE95+mSh2RiIiIiIjeYCc3LJXnYfFI+XKyMsGPfevh8OSW6OVVGYIA7L0Wi3aLTiBg+1U8Sc6QOiIREREREZUiFo/0Tk5WJlj4oSf2fdYCfrVskKMQsfl8JFr/eAy/HL2HjJc5UkckIiIiIqJSwOKRiqSmnRlWD2mArWN84elYAWlZOfhh/220XXQce6/GQIunCyUiIiIi0gosHqlYGjhZYvvYJljczxN25oaISniBsX+EYkjQBUQ8S5M6HhERERERqQmLRyo2mUxAj/qVcWRKS0xs4wp9HRlO3HmKtotO4OdDd5GZza6sRERERETlDYtHKjFjfV1MblcT+ye1QPPqVsjKVmDRoTvo9PNJhEQ8lzoeERERERGpEItHem/OViZYP7whlg6oD2szA9x/moY+K85g9r83kJ6VLXU8IiIiIiJSARaPpBKCIKBrvUo4NOnV/JCiCKw9HY4Oi0/ifHiC1PGIiIiIiOg9CaIWD5OZnJwMuVyOpKQkmJubSx2nXDl+5ykCtl3B46QMCAIwpqULJvnVgL4uv68gIs3ANiA3vh9ERNqrqG0A/ydPatGyhjX2T2qBvv9/FXL5sfvo8ctp3IlLkToaEREREWmB9KxsOE3dDaepu3krlYqweCS1MTPUww9962HFIC9YGOvhRkwyui49hc3nIzkvJBERERFRGcPikdSuQx177J/UAi1qWCMzW4GA7Vfx6eZLSMl4KXU0IiIiIiIqIhaPVCpszAyxbmgDTO3oBl2ZgP+uxKDL0lO4Fp0kdTQiIiIiIioCFo9UamQyAWNauuCvMb6oXMEIEc/S0Xv5GWwLeSR1NCIiIiIiKgSLRyp1XlUssGdic7Su+aob6+dbL2PajmvIylZIHY2IiIiIiArA4pEkITfWw5ohDeDvVx2CAGw4G4H+q4LxJCVD6mhERERERJQPXakDkPaSyQT4+9VAPYcK+OzPSwiNTESPZafx2xAfuFeSSx2PiAi7du0q9jFt27aFkZGRGtIQERFJi8UjSa61mw12TmiGEb9fwIOnaeizPBiL+tVDhzr2UkcjIi3Xo0ePYu0vCALu3r2LatWqqScQERGRhNhtlTSCs5UJ/hnXFM2rW+HFyxyM2RiKX4/d43yQRCS52NhYKBSKIi3GxsZSxyUiIlIbFo+kMeRGegga2gBDmzgBAL7fdxvf7riG7BwOpENE0hgyZEixuqAOGjQI5ubmakxEREQlEZvEcTVUgcUjaRRdHRlmdnPHjK61IQjAH+ciMWpDCNIys6WORkRaKCgoCGZmZkXef/ny5bCyslJjIiIiKqo3p4PzW3gcWy5ESpimfGDxSBppWFNnLB/oDQNdGY7ceoIBv53Fs9RMqWMRESklJydjx44duHnzptRRiIjoLTFJLzBj13XlukIEvt5+DTFJLyRMVfaxeCSN1aGOHTaPagxLE31ceZSEviuCEZ3If/BEJI0PP/wQy5YtAwC8ePECPj4++PDDD1G3bl1s27ZN4nRERPSm8Pg0KN4aOiNHFPEwPl2aQOUEi0fSaF5VLLB1jC8qVzDCg/g09P71DO7GpUgdi4i00IkTJ9C8eXMAwD///ANRFJGYmIglS5Zgzpw5EqcjIqI3OVuZQCbk3qYjCHCy4sBm74PFI2k8F2tT/D3WF642pohNzkDflcG4HJUodSwi0jJJSUmwtLQEAOzbtw+9e/eGsbExOnfujLt370qcjoiI3mQvN8Ksbu7KdZkAzOtVB/ZyzsP7Plg8UplgLzfC1tG+qOdYAYnpLzFw9TlcfJggdSwi0iKOjo4IDg5GWloa9u3bh3bt2gEAnj9/DkNDQ4nTEZG2Sc/KhtPU3XCauhvpWRxYMD+9vR2UPx+a3BL9GlSRME35oDHFY0pKCvz9/VG1alUYGRmhSZMmuHDhwjuPOXbsGLy8vGBgYABXV1esW7eudMKSJCxM9LFpZCM0rmaJ1MxsDF57Hmfux0sdi4i0hL+/PwYOHAgHBwdUqlQJrVq1AvCqO6uHh4e04YiI6J3s5PySTxU0pngcOXIkDh48iA0bNuDq1ato164d/Pz8EB0dne/+4eHh6Ny5M1q3bo2wsDD4+/tj5MiR2L9/fyknp9JkYqCLoKEN0by6FdKzcjAs6AKO3X4idSwi0gLjxo1DcHAw1q5di1OnTkEme9WEVqtWjfc8EhGRVtCI4vHFixfYtm0bvv/+e7Ro0QKurq6YOXMmXF1dsXz58nyPWbFiBZydnfHTTz+hVq1amDBhAvr06YNFixaVcnoqbUb6OvhtsA/8atkgM1uBUetDcJQFJBGpSfPmzfHjjz/i7t278PHxQc+ePWFqaqp8vHPnzmjatKmECYmIiEqHRhSP2dnZyMnJyXPPiJGREU6dOpXvMcHBwfDz88u1rX379ggODi7weTIzM5GcnJxrobLJUE8Hvw70Rgd3O2TlKDB6QwiO33kqdSwiKoc++eQTBAcHw8vLC7Vq1cJXX32F06dPQxTFwg8mIiIqRzSieDQzM4Ovry++++47PH78GDk5Odi4cSOCg4MRExOT7zGxsbGwtbXNtc3W1hbJycl48SL/uQADAwMhl8uVi6Ojo8pfC5UefV0Zln5UH+1q2yIrW4FR6y/i1F3eA0lEqjV48GBs27YN8fHx+Omnn5CYmIi+ffvCzs4Ow4cPx44dOwpsd4iIiMoTjSgeAWDDhg0QRRGVK1eGgYEBlixZggEDBijvKVGFgIAAJCUlKZeoqCiVnZukoacjw7KPvOBXyxaZ2QqM+P0CB9EhIrUwMDBAp06dsHLlSjx+/Bi7du2Cvb09pk2bhooVK6JLly44ffq01DGJiIjURmOKRxcXFxw/fhypqamIiorC+fPn8fLlS1SrVi3f/e3s7BAXF5drW1xcHMzNzWFklP/8LQYGBjA3N8+1UNmnryvDLwPro43bq3sgR/5+EaGRz6WORUTlXKNGjTB37lxcvXoVV69exQcffFBgbxkiIqLyQGOKx9dMTExgb2+P58+fY//+/ejevXu++/n6+uLw4cO5th08eBC+vr6lEZM0jIGuDpYP8kIz11ejsA5dex43Y3hPKxGpR2pqaq77562trTFp0iT06dNH6mhERERqozHF4/79+7Fv3z6Eh4fj4MGDaN26Ndzc3DBs2DAAr7qcDh48WLn/mDFj8ODBA3z55Ze4desWfv31V/z111+YNGmSVC+BJGagq4NVg73hXdUCyRnZ+HjNOdx/mip1LCIqJ15PEWViYgK5XA4LCwtYWFigQoUKsLCwkDoeERG9xVhfFw/nd8bD+Z1hrK8rdZxyQWPexaSkJAQEBODRo0ewtLRE7969MXfuXOjp6QEAYmJiEBkZqdzf2dkZu3fvxqRJk/Dzzz/DwcEBq1evRvv27aV6CaQBjPV1sXZoAwxYdRY3YpIxaPU5bBvbBJUq5N+VmYioqAYNGgRRFLF27VrY2tpCEASpIxERAQBikzJQzdq08B2J3pMgavFY48nJyZDL5UhKSuL9j+XMs9RMfLgyGPefpsHVxhRbR/vCwkRf6lhEpEGK2waYmpoiJCQENWvWLIV0pY9tIlHZsiH4IabtvA4AkAlAYC8P9GtQRdJM6VnZqD19PwDgxuz2vNpXhhS1DdCYbqtEqlTR1ADrRzSCvdwQ956kYvjvF5CelS11LCIqwxo0aMBRuolII8QkvcCMXdeV6woR+Hr7NcQkcdogUi9+HUDlVuUKRlg/vCH6rAjGpchEjPsjFL8N9oGeDr8zIaLiW716NcaMGYPo6GjUqVNHeVvFa3Xr1pUoGRFpm/D4NCje6juYI4p4GJ8Oezlv1SH1YfFI5Vp1WzOsHdoAA1efxbHbTzF121X82Lcu71UiomJ7+vQp7t+/rxzIDQAEQYAoihAEATk5ORKmIyJt4mxlApmAXAWkjiDAycpYulCkFXgJhso976oW+HWgF3RkAraFPsLiQ3eljkREZdDw4cNRv359BAcH48GDBwgPD8/1JxFRabGXG2FWN3flukwA5vWqw6uOpHa88khaoY2bLeb0qIOA7Vfx8+G7qFzBCB82cJQ6FhGVIREREdi1axdcXV2ljkJEhN7eDsoBcw5NbsnRVqlU8MojaY0BDatgQutX/+kL+Ocqjt95KnEiIipL2rRpg8uXL0sdg4goDzu5odQRSEvwyiNplc/b1UB04gv8cyka4zaG4O+xTVDLnkPSE1HhunbtikmTJuHq1avw8PDIM2BOt27dJEpGRERUOnjlkbSKIAhY0LsufKtVRFpWDkasu4AnKRlSxyKiMmDMmDF49OgRZs+ejb59+6JHjx7KpWfPnmp97hMnTqBr166oVKkSBEHAjh07cj0uiiKmT58Oe3t7GBkZwc/PD3fv8v5uIiJSLRaPpHX0dWVYMcgb1axM8DgpA5+sD0HGS46SSETvplAoClzUPdJqWloa6tWrh19++SXfx7///nssWbIEK1aswLlz52BiYoL27dsjI4NfjhERkeqweCStJDfWw9qhDVDBWA+XoxLx+V+XoXh7wiQiojc8evSowMfOnj2r1ufu2LEj5syZk+8VTlEUsXjxYnz77bfo3r076tati/Xr1+Px48d5rlASEZWW2CR+eVUesXgkreVkZYIVg7yhpyNg99UYLDp0R+pIRKTB2rVrh4SEhDzbT58+jQ4dOkiQ6JXw8HDExsbCz89PuU0ul6NRo0YIDg4u8LjMzEwkJyfnWoiI3se2kP99yea38Di2XIiUMA2pA4tH0mqNq1XEvJ4eAIClR+7h38uPJU5ERJqqcePGaNeuHVJSUpTbTpw4gU6dOmHGjBmS5YqNjQUA2Nra5tpua2urfCw/gYGBkMvlysXRkdMXEVHJxSS9wIxd15XrChH4evs1xCS9yHf/9KxsOE3dDaepu5GelV1aMcs8qd83Fo+k9fr6OGJUi2oAgC/+voyrj5IkTkREmmj16tWoUqUKunbtiszMTBw9ehSdO3fG7NmzMWnSJKnjFVtAQACSkpKUS1RUlNSRiKgMC49Pw9t3AOWIIh7Gp0sTiNSCxSMRgK86uKFVTWtkvFRg1IaLHIGViPKQyWT4888/oaenhzZt2qBbt24IDAzEZ599JmkuOzs7AEBcXFyu7XFxccrH8mNgYABzc/NcCxGVHcb6ung4vzMezu8MY33pZ99ztjKBTMi9TUcQ4GRlLE0gUgsWj0QAdGQClgyoj2rWJohJysDoDSHIzOYIrETa7sqVK7mWW7duYebMmYiKisKgQYPQokUL5WNScXZ2hp2dHQ4fPqzclpycjHPnzsHX11eyXET5kbrLHamPvdwIs7q5K9dlAjCvVx3Yy40kTEWqJv3XFEQawtxQD2uGNED3ZadwKTIR03dcx/zeHhAEofCDiahc8vT0hCAIEMX/9cV6vb5y5UqsWrUKoihCEAS1TteRmpqKe/fuKdfDw8MRFhYGS0tLVKlSBf7+/pgzZw6qV68OZ2dnTJs2DZUqVUKPHj3UlomI6G29vR0wbeer+x4PTW6JatamEiciVWPxSPQGZysTLP3IC8OCzmPLxSjUcZDj48ZVpY5FRBIJDw+XOgIA4OLFi2jdurVyffLkyQCAIUOGYN26dfjyyy+RlpaGUaNGITExEc2aNcO+fftgaGgoVWQi0nJ2cv7+KY9YPBK9pWUNa3zR3g0L9t3CrF3X4WZnhgZOllLHIiIJVK2qGV8etWrVKtfVz7cJgoDZs2dj9uzZpZiKiIi0De95JMrHmJbV0LmuPbIVIsZuDC1wmGkiKt927dqFly9fFnn/PXv24MUL/r4gIqLyicUjUT4EQcAPferCzc4M8amZGLMxlAPoEGmhnj17IjExscj79+/fHzExMeoLREREJCF2WyUqgLG+LlZ97IOuy07hclQiZv97A3N7ekgdi4hKkSiKGDp0KAwMDIq0f0YGp/khIqLyi8Uj0TtUqWiMn/t7Yti6C/jjXCQ8HSugr4+j1LGIqJQMGTKkWPsPHDiQ8yUSERVTbFIGR2YtASneNxaPRIVoVdMGk/xqYOHBO/hmxzXUsjdHncpyqWMRUSkICgqSOgIRUbm0LeSR8me/hccR2MsD/RpUkTBR2SD1+8Z7HomKYEJrV3zgZoOsbAVGbwjB87QsqSMRERERlUkxSS8wY9d15bpCBL7efo0DFBZCE943Fo9ERSCTCVjYzxNVKxojOvEFJv0VBoWi4GHziYiI6N1ik3iPcHljrK+Lh/M74+H8zjDWL7iDY3h8Gt7+b1SOKOJhfLrasqVnZcNp6m44Td2N9KxstT2POknxvr2NxSNREcmN9LB8oDcMdGU4dvsplh29J3UkIiKiMuXtLndbLkRKmIak4mxlApmQe5uOIMDJyliaQGWEJrxvLB6JiqF2JXPM6VEHALDo0B2cvPtU4kRERERlgyZ0uSPNYC83wqxu7sp1mQDM61UH9nIjCVNpPk1431g8EhVTXx9H9G/gCFEEPvszDI8T2egRaYMHDx5IHYGoTNOELnekOXp7Oyh/PjS5JQfLKSKp3zcWj0QlMLObO9wrmSMhLQvjN4UiK1shdSQiUjNXV1e0bt0aGzdu5HyORCWgCV3uSDPZyQ2ljlAmSfG+sXgkKgFDPR0sH+gNc0NdXIpMxIJ9t6SORERqFhoairp162Ly5Mmws7PD6NGjcf78ealjEZUZmtDljojeD4tHohKqUtEYP/atBwBYcyoc+67FSpyIiNTJ09MTP//8Mx4/foy1a9ciJiYGzZo1Q506dbBw4UI8fcp7oIkKI3WXOyJ6Pyweid5DO3c7jGpRDQDwxdbLiHiWJnEiIlI3XV1d9OrVC1u3bsWCBQtw7949TJkyBY6Ojhg8eDBiYmKkjkhUJrCrIlHZw+KR6D190b4mfKpaICUzG+P+CEXGyxypIxGRGl28eBHjxo2Dvb09Fi5ciClTpuD+/fs4ePAgHj9+jO7du0sdkYiISC1YPBK9Jz0dGZZ+VB+WJvq4/jgZs/+7IXUkIlKDhQsXwsPDA02aNMHjx4+xfv16REREYM6cOXB2dkbz5s2xbt06hIaGSh2ViIhILTSieMzJycG0adPg7OwMIyMjuLi44LvvvoMoigUec+zYMQiCkGeJjeV9Z1T67OVGWNzPE4IAbDoXiZ1h0VJHIiIVW758OT766CNERERgx44d6NKlC2Sy3M2ojY0N1qxZI1FCIiIqitgkjphdUrpSBwCABQsWYPny5fj999/h7u6OixcvYtiwYZDL5Zg4ceI7j719+zbMzc2V6zY2NuqOS5SvFjWsMaG1K5YeuYevt19FncpyuFibSh2LiFTk4MGDqFKlSp6CURRFREVFoUqVKtDX18eQIUMkSkhERAXZFvJI+bPfwuMI7OXBAZtKQCOKxzNnzqB79+7o3LkzAMDJyQmbN28u0hDoNjY2qFChgpoTEhWNv18NXHiYgLMPEjD+j1D8M64pjPR1pI5FRCrg4uKCmJiYPF9SJiQkwNnZGTk5vN+ZiEgTxSS9wIxd15XrChH4evs1tKhhXeamijHW18XD+Z0le36N6LbapEkTHD58GHfu3AEAXL58GadOnULHjh0LPdbT0xP29vZo27YtTp8+/c59MzMzkZycnGshUiUdmYAl/evDytQAt2JTMPONX1REVLYVdCtFamoqDA05aiQRUXG8LoIezu8MY331Xs8Kj0+D4q1f4TmiiIfx6Wp93vJII648Tp06FcnJyXBzc4OOjg5ycnIwd+5cDBw4sMBj7O3tsWLFCvj4+CAzMxOrV69Gq1atcO7cOXh5eeV7TGBgIGbNmqWul0EEALAxN8SS/p4YuOYctlyMQqNqlujl5VD4gUSkkSZPngwAEAQB06dPh7GxsfKxnJwcnDt3Dp6enhKlIyKiwjhbmUAmIFcBqSMIcLIyLvggypcgvmtUmlLy559/4osvvsAPP/wAd3d3hIWFwd/fHwsXLizWvSMtW7ZElSpVsGHDhnwfz8zMRGZmpnI9OTkZjo6OSEpKynXfJJEqLD50B4sP3YWRng7+/bQpXG3MpI5ERG9ITk6GXC4vtA1o3bo1AOD48ePw9fWFvr6+8jF9fX04OTlhypQpqF69utozq1NR3w8iorJoQ/BDTNv5qkeYTADveXxLUdsAjbjy+MUXX2Dq1Kno378/AMDDwwMREREIDAwsVvHYsGFDnDp1qsDHDQwMYGBg8N55iYri0zbVceFhAk7fe4Zxf4Ri5/hmvP+RqAw6evQoAGDYsGH4+eefWVgREZVBvb0dlMXjocktUY2DGpaIRtzzmJ6enmf0Oh0dHSgUimKdJywsDPb29qqMRlRiOjIBi/vVh7WZAe7EpWL6zmtSRyKi9xAUFMTCkYioHLCT8z71ktKIK49du3bF3LlzUaVKFbi7u+PSpUtYuHAhhg8frtwnICAA0dHRWL9+PQBg8eLFcHZ2hru7OzIyMrB69WocOXIEBw4ckOplEOVhbWaAn/t7YtDqc9ga8giNqlVEH2/e/0hUVvTq1Qvr1q2Dubk5evXq9c59t2/fXkqpiIiIpKERxePSpUsxbdo0jBs3Dk+ePEGlSpUwevRoTJ8+XblPTEwMIiMjletZWVn4/PPPER0dDWNjY9StWxeHDh1S3ptCpCmauFjB368GFh68g2k7rqGegxzVbXn/I1FZIJfLIQiC8mciIiJtphED5kiFgwNQaclRiBiy9jxO3YtHdRtT7JzQVO3DUhPRu7ENyI3vBxGVZ+lZ2ag9fT8A4Mbs9vx/2FvUMmDOrl27ih2kbdu2MDIqW5NvEqmajkzA4v6e6PTzSdx9kopvd1zDT33rKa9oEJHmCw8PR3Z2dp5RVe/evQs9PT04OTlJE4yIiKiUFKt47NGjR7FOLggC7t69i2rVqhXrOKLyyMrUAEsG1MdHv53F9tBoNHauiA8bOEodi4iKaOjQoRg+fHie4vHcuXNYvXo1jh07Jk0wIiKiUlLs0VZjY2OhUCiKtLw5kTIRAY2rVcTn7WoCAKbtvIabMckSJyKiorp06RKaNm2aZ3vjxo0RFhZW+oGIiIhKWbGKxyFDhhSrC+qgQYN43wTRW8a2dEHLGtbIzFZg/B+hSMl4KXUkIioCQRCQkpKSZ3tSUhJycnIkSERERFS6ilU8BgUFwcys6KNELl++HFZWVsUORVSeyWQCFvXzhL3cEA/i0zB1+1Vo8bhVRGVGixYtEBgYmKtQzMnJQWBgIJo1ayZhMiIiotLBYYaIJGBpoo9lH3mh38pg7L4Sg0bOlhjs6yR1LCJ6hwULFqBFixaoWbMmmjdvDgA4efIkkpOTceTIEYnTEVF5wpFBVc9YXxcP53eWOkaZV+x7Hl978eIF0tPTlesRERFYvHgxDhw4oJJgROWdd1ULBHSqBQD47r8buByVKG0gInqn2rVr48qVK/jwww/x5MkTpKSkYPDgwbh16xbq1KkjdTwiIiK1K/HXGN27d0evXr0wZswYJCYmolGjRtDT00N8fDwWLlyIsWPHqjInUbk0vKkTLoQnYN/1WIz7IxT/fdoMFib6UsciogJUqlQJ8+bNkzoGUbnAq2tEZU+JrzyGhoYqu+38/fffsLW1RUREBNavX48lS5aoLCBReSYIAr7vWxdOFY0RnfgC/lvCoFDw/kciTZWYmIiffvoJI0eOxMiRI7Fo0SIkJSVJHYuIiKhUlLh4TE9PVw6ec+DAAfTq1QsymQyNGzdGRESEygISlXfmhnpYPsgbhnoyHL/zFEuP3JM6EhHl4+LFi3BxccGiRYuQkJCAhIQELFy4EC4uLggNDZU6HhERkdqVuHh0dXXFjh07EBUVhf3796Ndu3YAgCdPnnB6DqJiqmVvjjk9PAAAiw/fwfE7TyVORERvmzRpErp164aHDx9i+/bt2L59O8LDw9GlSxf4+/tLHY+IiEjtSlw8Tp8+HVOmTIGTkxMaNmwIX19fAK+uQtavX19lAYm0RR9vBwxoWAWiCPj/eQmPnqcXfhARlZqLFy/iq6++gq7u/+7L0tXVxZdffomLFy9KmIyIiKh0FLt4nD59OkJCQtCnTx9ERkbi4sWL2L9/v/LxDz74AIsWLVJpSCJtMaNrbXhUluN5+kuM2RiCjJeceJxIU5ibmyMyMjLP9qioqGLNgUxERFRWFbt4fPToETp27AgHBwfMmjULsbGxyM7OVj7esGFDuLm5qTQkkbYw1NPB8kFesDDWw7XoZHy74xpEkQPoEGmCfv36YcSIEdiyZQuioqIQFRWFP//8EyNHjsSAAQOkjkekkdKzsuE0dTecpu5GelZ24QcQkUYr9pjIa9euhUKhwOnTp/Hvv//C398fMTExaNu2Lbp3744uXbrA0tJSHVmJtIKDhTGWDvDC4LXn8HfII3g6VsCgxlWljkWk9X788UcIgoDBgwcrvzTV09PD2LFjMX/+fInTERERqV+J7nmUyWRo3rw5vv/+e9y+fRvnzp1Do0aNsHLlSlSqVAktWrTAjz/+iOjoaFXnJdIKzapb4csOr67gz/r3OkIinkuciIj09fXx888/4/nz5wgLC0NYWBgSEhKwaNEiGBgYSB2PSKto0xXN2KSMAh/TpvehKPh+qF+JB8x5U61atfDll1/i9OnTiIqKwpAhQ3Dy5Els3rxZFacn0kqjW1RDxzp2eJkjYuzGEMQlF9x4EFHpMTY2hoeHBzw8PGBsbCx1HCIqh7aFPFL+7LfwOLZcyHu/NZEUit1ttTDW1tYYMWIERowYoepTE2kVQRDwQ996uP80FXfiUjF6Qwj+HNUYhno6Ukcj0hq9evUq8r7bt29XYxKi8i02KQPVrE2ljqERYpJeYMau68p1hQh8vf0aWtSwhr3cSMJkRO9ZPGZkZODKlSt48uQJFApFrse6dev2XsGICDA10MVvg33QbdlphEUl4tsd1/BDn7oQBEHqaERaQS6XSx2BqNx6++paYC8P9GtQRcJEmiE8Pg2Kt8bKyxFFPIxPZ/FIkitx8bhv3z4MHjwY8fHxeR4TBAE5OZxigEgVqlY0wbKP6mPI2vP4O+QR3CuZY1hTZ6ljEWmFoKAgqSMQlUu8ulYwZysTyATkKiB1BAFOVuW7m3x6VjZqT381/d+N2e1hrK/yDpKkAiW+5/HTTz9F3759ERMTA4VCkWth4UikWs2rW+PrTrUAAHN238SJO08lTkSknbKzs3Ho0CGsXLkSKSkpAIDHjx8jNTVV4mREZcu7rq5pO3u5EWZ1c1euywRgXq86Wl9Uk2YocfEYFxeHyZMnw9bWVpV5iKgAI5o5o7eXA3IUIsZvCsW9JylSRyLSKhEREfDw8ED37t0xfvx4PH366kucBQsWYMqUKRKnIypbXl9de5M2XF0rqt7eDsqfD01uye68pDFKXDz26dMHx44dU2EUInoXQRAwr1cdNHCyQEpGNoavu4iEtCypYxFpjc8++ww+Pj54/vw5jIz+dwWgZ8+eOHz4sITJiMoeXl0rOju5odQRiJRK3Jl42bJl6Nu3L06ePAkPDw/o6enlenzixInvHY6IcjPQ1cHKj33Q/ZdTiExIx5gNIdgwsiEMdDkCK5G6nTx5EmfOnIG+vn6u7U5OTpzXmKgI3h5Rtbe3A6btfHXf46HJLTnaKlEZUOLicfPmzThw4AAMDQ1x7NixXKM/CoLA4pFITSxN9LF2SAP0+vUMzj9MQMD2q/ipbz2OwEqkZgXd0//o0SOYmZlJkIhI8xV1RNX3ubrGaT5e4fuQG98P9Shxt9VvvvkGs2bNQlJSEh4+fIjw8HDl8uDBA1VmJKK3VLc1wy8DvaAjE7A9NBqLDt6ROhJRudeuXTssXrxYuS4IAlJTUzFjxgx06tRJumBv+OWXX+Dk5ARDQ0M0atQI58+flzoSabGCRlSNSXrx3ud+uyh1mrobTlN3Iz0r+73PXZa8/T5suRApYRrp8f1QvxIXj1lZWejXrx9kshKfgojeQ4sa1pjXsw4AYMmRe/jzPH9BEqnTTz/9hNOnT6N27drIyMjARx99pOyyumDBAqnjYcuWLZg8eTJmzJiB0NBQ1KtXD+3bt8eTJ0+kjkZaSl0jquZXlGojdRbnZRHfj9JR4spvyJAh2LJliyqzEFEx9WtQBRPbuAIAvtlxDUdv8z+JROri4OCAy5cv4+uvv8akSZNQv359zJ8/H5cuXYKNjY3U8bBw4UJ88sknGDZsGGrXro0VK1bA2NgYa9eulToaaSl1jaiaX1GqjTjdSW58P0pHie95zMnJwffff4/9+/ejbt26eQbMWbhw4XuHI6LCTWpbA9GJGdgW+gjj/wjFpk8aw9OxgtSxiMqdjIwMGBoaYtCgQVJHySMrKwshISEICAhQbpPJZPDz80NwcHC+x2RmZiIzM1O5npycrPacpF1ej6j6elAcVY2o+roo1fYCMr/3QZunO+H7UTpKfOXx6tWrqF+/PmQyGa5du4ZLly7lWoiodAiCgMBeHmhe3QrpWTkYGnQed+M4BySRqtnY2GDIkCE4ePAgFAqF1HFyiY+PR05OTp65l21tbREbG5vvMYGBgZDL5crF0dGxNKKSllHHfIX5TfOhjcrzdCexSRnFPqY8vx+apMTF49GjR9+5EFHp0deVYfkgb9RzrIDE9JcYtOYcohLYTYNIlX7//Xekp6eje/fuqFy5Mvz9/XHx4kWpY5VYQEAAkpKSlEtUVJTUkaicU+V8hW8Wpf9+2lRl59UUxvq6eDi/Mx7O7wxj/YI7CqqjOJeKKga7KU/vh6YqcfG4efPmAh/74osvSnpaIiohUwNdrBvaANVtTBGXnImP15zD05TMwg8koiLp2bMntm7diri4OMybNw83btxA48aNUaNGDcyePVvSbFZWVtDR0UFcXFyu7XFxcbCzs8v3GAMDA5ibm+daiMoiW3PVFaVlmSqL89KmjsFuyvL7oclKXDyOHTsWe/fuzbN90qRJ2Lhx43uFIqKSsTDRx4YRjeBgYYSHz9IxaPU5JKRlSR2LqFwxMzPDsGHDcODAAVy5cgUmJiaYNWuWpJn09fXh7e2Nw4cPK7cpFAocPnwYvr6+EiYjKlhRr66pQnpWttZO51EWcLCbgmnaZ7fExeMff/yBAQMG4NSpU8ptn376Kf766y92WyWSkJ3cEBtHNIK1mQFux6Xgo9/O4jkLSCKVycjIwF9//YUePXrAy8sLCQkJGtHjZvLkyfjtt9/w+++/4+bNmxg7dizS0tIwbNgwqaMREb2TukbmJdUrcfHYuXNn/Prrr+jWrRtCQkIwbtw4bN++HUePHoWbm5sqMxJRMTlZmWDzJ41hZWqAW7EpGLj6HBLTWUASvY/9+/djyJAhsLW1xdixY2Fra4sDBw4gIiIC8+fPlzoe+vXrhx9//BHTp0+Hp6cnwsLCsG/fvjyD6BCVZyUZaIWkx8Fuyo4SF48A8NFHH2HOnDlo2rQp/v33Xxw/fhw1atQo9nlycnIwbdo0ODs7w8jICC4uLvjuu+8giu8eg/nYsWPw8vKCgYEBXF1dsW7duhK+EqLyx9XGFJs/aQQrU33ciEnGoDUsIIneR8+ePfHixQusX78esbGxWLlyJVq0aCF1rFwmTJiAiIgIZGZm4ty5c2jUqJHUkYjUbselx8qfSzrQCkmPg92UDcXqYD558uR8t1tbW8PLywu//vqrcltx5nlcsGABli9fjt9//x3u7u64ePEihg0bBrlcjokTJ+Z7THh4ODp37owxY8bgjz/+wOHDhzFy5EjY29ujffv2xXlZROVWdVszbPqkMQasOotr0cnov+os1o9oCBsz3kROVFxxcXEwMzOTOgaRVkrPykbt6fsBADdm5/5/3rw9N5U/vx5opUUNa161ek9vv+fqvi/1TRzsRnMV61NQ0PyNrq6uSE5OVj4uCMWbcOfMmTPo3r07OnfuDABwcnLC5s2bcf78+QKPWbFiBZydnfHTTz8BAGrVqoVTp05h0aJFLB6J3lDD1gybRzXGoNXncCs2BX1XBGPjiEZwtOR9BETFwcKRSDMVNNAKi0ci1StW8aiugXCaNGmCVatW4c6dO6hRowYuX76MU6dOvfPqZXBwMPz8/HJta9++Pfz9/Qs8JjMzE5mZ/5u6IDk5+b2zE5UFNWzN8PeYJhi45iwinqW/KiBHNoSrDf8zTERE6vN6RFV1kgm5C0gOtEKkPsW65/HKlStQKBRF3v/69evIzi58SNmpU6eif//+cHNzg56eHurXrw9/f38MHDiwwGNiY2PzDAJga2uL5ORkvHiR/5wwgYGBkMvlysXR0bHIr4WorKtS0RhbRzdBdRtTxCZnoM+KYJwPT5A6FhER0Xv5ulMt5c8caIVIvYpVPNavXx/Pnj0r8v6+vr6IjCz8puW//voLf/zxBzZt2oTQ0FD8/vvv+PHHH/H7778XJ16hAgICkJSUpFyioqJUen4iTWcnN8SW0b7wdKyAxPSXGLT6HHaGRUsdi4iIqFjenCPyo0b/G1hF2wZaKc25MssCdb8fmjbnohSK9a6Koohp06bB2LhoXQGysoo2suMXX3yhvPoIAB4eHoiIiEBgYCCGDBmS7zF2dnaIi4vLtS0uLg7m5uYwMsr/2yYDAwMYGBgUKRNReWVpoo/NnzTGpC1h2Hc9Fp/9GYbIZ+mY0Ma12PcrE2mziIgIpKWlwc3NDTLZew1eTkQqUpSBVmKTMlDN2rQU0qiWlAPYEL1WrE9dixYtcPv27SLv7+vrW2Ah96b09PQ8Da+Ojs47u8j6+vpiz549ubYdPHgQvr6+Rc5HpK2M9HXw60AvzN93C6tOPMBPB+/gdlwKFvSuCxMDNkZEb1q7di0SExNzjTg+atQorFmzBgBQs2ZN7N+/n7dCEGmwbSGPlD/7LTyOwF4eWnWFksoHTfjio1j/Szx27JhaQnTt2hVz585FlSpV4O7ujkuXLmHhwoUYPny4cp+AgABER0dj/fr1AIAxY8Zg2bJl+PLLLzF8+HAcOXIEf/31F3bv3q2WjETljUwm4OtOtVDF0hgzd13Hf1dicCcuBSsGeUv+i4lIk6xatQqjR49Wru/btw9BQUFYv349atWqhQkTJmDWrFlYvXq1hCmJqCAxSS8wY9d15Tqn86CyRNO++NCIfjZLly5Fnz59MG7cONSqVQtTpkzB6NGj8d133yn3iYmJyXX/pLOzM3bv3o2DBw+iXr16+Omnn7B69WpO00FUTIMaV8WfoxrDxswAd+JS0W3Zaey7FiN1LCKNcffuXfj4+CjXd+7cie7du2PgwIHw8vLCvHnzcPjwYQkTEtG7hMenFTidB5Emy++Lj4DtVxGTlP/goKVBI/qnmZmZYfHixVi8eHGB+6xbty7PtlatWhU49yQRFZ2PkyX+m9gME/64hPMPEzBmYyj6N3DEtC612Y2VtN6LFy9gbm6uXD9z5gxGjBihXK9WrRpiY2OliEZEReBsZcLpPIrpzfsrS0tpTOtS1uT3xYdChKTzmGrElUcikp6NmSH++KQRRresBkEA/rwQhU5LTiI08rnU0YgkVbVqVYSEhAAA4uPjcf36dTRt2lT5eGxsLORyuVTxiLRObFJGsfa3lxthVjd35bq2T+dRkhFDi/uek2q8/uLjTTIBkn7xweKRiJT0dGQI6FgLm0Y2RiW5ISKepaPvimDM23MTaZnaOSQ10ZAhQzB+/Hh899136Nu3L9zc3ODt7a18/MyZM6hTp46ECYnKv7fv+9pyofCp4N7U29tB+bO2TeehCiV5z+n9vf3FBwDM6uYu6RcfLB6JKA9fl4rY698CPTwrIUchYtWJB2i78DgOXGfXPNI+X375JT755BNs374dhoaG2Lp1a67HT58+jQEDBkiUjqj8K2jAm5Le91WU6Ty03dtXGt/3PaeSe/OLj/zWS5vKiscBAwYgOTkZALBr1648jSsRlS1yIz0s7l8fa4f6wMHCCI+TMjBqQwhGrLuAO3EpUscjKjUymQyzZ8/GpUuXsHfvXtSqVSvX41u3bs11DyQRqRYHvCl9Ec/S8mzje56btnblVdlIGNevX4e5uTlu3LiBr7/+Gi1btsSJEyewdOlSVT0FEUmgjZstfKtZYcmRu/jtxAMcvvUER28/QS8vB0xqWwOVK2jnPSOkfV68eIGDBw/izp07AIAaNWqgbdu2RZrPmIhKrrABbzjQiupVrWiSZ1t5GGTozYGAbsxuD2P94pVCmjZthhRUVjzq6elBFEUEBQUhICAAAwcOzHVPCBGVXUb6Oviqgxt6ezngx/23se96LP4OeYRdlx+jr7cDRjavBmervA0NUXmxa9cujBw5EvHx8bm2W1lZYc2aNejatatEyYjKt/SsbPgGHsm1TdsHvAHUP1n82117+Z5zvtDXVNZtdezYsfDy8sLff/+NHj16AADS0vJe8iaissvVxhQrPvbGP+OaoHE1S2RlK/DHuUi0+ekYRm+4iPPhCRBFsfATEZUhZ86cQZ8+fdCiRQucPn0aCQkJSEhIwKlTp9C8eXP06dMHZ8+elTomkdbQ1gFv3nfQoPehre/5m9h9+hWVXXkcOXIk+vTpA11dXZiYmODevXto3Lixqk5PRBqkfhULbP6kMc6HJ2DV/3dl3X89Dvuvx8GpojH6eDugl5cDKrFLK5UDc+bMwbBhw7By5cpc25s0aYImTZpg9OjRmD17Nvbs2SNRQiLtoo0D3kh91Usb3/O3acp8oeq+6lwYlV15fPnyJXbs2IHly5dj9+7dqFatGtatW6eq0xORhhEEAY2qVcSaoQ1waHIL9G/gCGN9HTx8lo4fD9xB0wVH0Hv5GSw/dh9341J4RZLKrLNnz2LChAkFPj5+/HgEBweXYiIi0ja86iU9TZkvVOppU1RWPPbv3x8XL16EkZER/vvvP3h5eSkHFSCi8s3Vxgzze9fFhW/88GPfemjkbAlRBEIinmPBvltou+gEmn9/FJO3hGHTuUjcjUtBztutIJGGevHiBczNzQt8XC6XIyNDO0fdI5JKSSa6L8vymyxenVe93uwim9+6tpJivlBNmzZFZd1WHzx4gG3btinXw8LCMHLkSJw4cUJVT0FEGs7EQBd9vB3Qx9sBjxNf4PCtJzh0Iw7B95/h0fMXePQ8GtsvRQMADPVkqGFrhhq2ZqhuY4rKFkaoXOHVYmmiD10dTkNLmqF69eo4cuQIhg0blu/jhw8fRvXq1Us5FRFpk9dXvabtfNV1VVVXvfLrAvl2F1kAmLHrOvxq22rVwDCFKa2uvO+aNkWKvw+VFY9mZma4d+8eXF1dAQCenp54/vy5qk5PRGVMpQpG+LhxVXzcuCrSMrMREvEcFx8m4MLD57gU9RwZLxW48igJVx4l5Xu8uaEuLE30ITfSg4GeDoz0dGCoJ4Oujgw6ggCZAMgEAQpRRI6IV3/miMhWKPAyR0SOQsTLHAWyFSKycxTIEUUoFK/2E0VAEF51vZUJgK6ODIa6Mhj+//OYG+migvGr565oog87uSHs5UawkxtCbqRXyu8kSW3YsGGYMmUKbG1t0alTp1yP7d69G19++SW+/vpridIRUVGUh+k8ens7KIvHQ5Nblvi+t8Kmm8ivi6xChGTFirZzszfXiHstX1NZ8bhs2TJ0794dnTp1Qu3atXHz5k1UrVpVVacnojLMxEAXLWpYo0UNawBAjkJExLM03I5Nwe24FDx4mobHiS8QnfgCcckZUIhAckY2kjM0ryuSpYk+XKxN4Gpjiuo2ZqjnWAHulcxhqKcjdTRSk88++wxnzpxBly5dULNmTdSqVQuiKOLmzZu4e/cuevToAX9/f6ljEpEWKelVr6IMvJPfwDAyAWV+jseySl1XnUtKZcVj3bp1ERoaih07duDmzZtwcXHBzJkzVXV6IipHdGQCqlmbopq1KTp62Od6LDtHgcQXL5GYnoXn6S+RlP4SGdk5eJGVg4xsBXJyFMgRAVEUoRBFyATh/xdAR0cGPZkAXR0ZdGUCdHUE6Mpe/ayj82o/HUGAIADi66uVooiX2QpkZiuQ8TIHGS9zkJyRjcT0LCSmv8SztCzEJGUgNukFnqe/REJaFhLSsnDh4f96VujKBLjZm6GRc0U0q26FRs6WxZ54mDSXTCbD1q1bsWXLFmzevBm3bt0CALi5uWHmzJno37+/xAmJiIrmXQPvvC5G3i5WAGBWN3dedZSQqq46q8J7/+/G398f9erVQ926dVGnTh3069dPFbmISEvp6shgZWoAK1MDqaPkkZ6VjQdP03D/aSruPUnFzZhkhEUlIj41C9eik3EtOhlrToVDX0eGBs4W6FjHHh3q2Gnka6Hi69evX75t3IMHDzBmzBgcOHBAglREREVX1Okm3ixWXq+XN1JPeVFSUk+bUuLiMTU1FaampmjTpg2uXLmCvXv34vr16xAEAe7u7qhbty7q1q2Lrl27qjIvEZFkjPV1UaeyHHUqy5XbRFHE46QMhEY8x+l78Th5Nx7RiS9w+t4znL73DNN3XkPjahXRo35ldK1bCUb67N5a3qSkpODw4cNSxyAql968P08mAIG9PGCsr6sVI6yqQ0m6QN6Y3b7c9KYp7H5PKlyJPwlyuRx//fUXevfujW7duim3Z2Rk4Nq1a7hy5QqOHDnC4pGIyjVBEJSjxHatVwmiKCI8Pg2HbsZh95UYXH6UhDP3n+HM/Wf47r8b6FW/Mj5qVBU17cykjk5EpNHedX8eBy8rOU3qAlmainK/JxWuxMWjKIpYuXIlFi5cCEEQ4OPjg4EDB6JBgwbw8fGBj4+PKnMSEZUJgvDqfs5R1qYY1cIFUQnp2HX5Mf68EImohBf4PTgCvwdHoFVNa4xr5YqGzpZSRyYikkR6VjZqT98PIP+rW++6P6+eoxz0/qTuAlmainK/JxXuvSZSu3TpEry8vNCsWTNcv34dzZs3x5QpU1SVjYiozHO0NMb41q44PqU11g9viPbutpAJwLHbT/HhymD0XXEGJ+8+lTomEZHGeX1/3puknKKAyjZ+nlTjvTowb9q0CW3btlWuX7lyBd27d0flypUxadKk9w5HRFReyGSCcrqSh/FpWHniAbaFPMKFh8/x8ZrzaF7dCgEda6F2JXOpo9Jb6tevD0EQCnw8PT29FNMQaY933Z/35j2PZXXgEypdqpjyojzMF/q+Slw8WlpawtHRMde2unXrYtmyZZg0aRKLRyKiAjhZmSCwlwf8/apjxfH72Hg2AifvxuPUvZPoWb8ypnZwg4259nQl0nQ9evSQOgKR1iro/jxtHPiEhcv709b7PVWpxMWjp6cngoKCsGDBglzbXV1dERkZ+d7BiIjKO1tzQ8zo6o5hTZzxw4Hb+PfyY2wPjcbBG3H4soMbBjasAtnbfWyo1M2YMUPqCESE/92fx4FPSBW06X5PVSrxPY9z5szBkiVL8PHHHyM4OBhpaWl48uQJ5s2bB2dnZ1VmJCIq16pUNMbSAfWxc3xT1HOQIyUjG9N2XEPvFWdwKzZZ6nhERBrlXQOfEJF6lfjKY+PGjXH27Fl89tlnaN68OUTx1b9iQ0NDbN26VWUBiYi0RT3HCtg+rik2no3AD/tv41JkIrouPYUp7WpiZPNq0OFVSCKiIk90T1ReaFKX5fcaMKdevXo4duwYnjx5gpCQECgUCjRq1AhWVlaqykdEpFV0ZAKGNHFCe3c7fLvjGg7djEPg3ls4evsJfvrQE5UrsEsWEWk3VQx8Qu+mScUKaZb3mqrjNRsbG3Ts2BGdO3dm4UhEpAJ2ckP8NtgbC3p7wFhfB2cfJKDD4hPYdy1W6mhERJLr7e2g/PnQ5JblfrAcVXpdGD6c3znP3JpEhVFJ8UhERKonCAL6NaiCPRObw9OxAlIysjFmYwjm772F7ByF1PGIiDQCBz4hKj38uoGISMM5WZlg6xhfLNh7C6tPhWPF8fu48igRSwbUh5WpgdTxyrUlS5YUed+JEyeqMQlR+ca5GonKBhaPRERlgJ6ODN92qQ3PKhXw5d9XcOb+M3RbegprhzWAm5251PHKrUWLFuVaf/r0KdLT01GhQgUAQGJiIoyNjWFjY8PikaiYtHGuRqKyjt1WiYjKkC51K2Hn+KaoZmWCx0kZ6LM8GMfvPJU6VrkVHh6uXObOnQtPT0/cvHkTCQkJSEhIwM2bN+Hl5YXvvvtO6qhEZUpBczXGJL2QMBURFYbFIxFRGVPd1gz/jGuKxtUskZqZjeHrLuCPcxFSxyr3pk2bhqVLl6JmzZrKbTVr1sSiRYvw7bffSpiMqOwp6lyNHNyFVImfp/fH4pGIqAySG+th/fBG6OVVGTkKEd/8cw0/7L+lnHOXVC8mJgbZ2dl5tufk5CAuLk6CRERl1+u5Gt/EuRqJNB+LRyKiMkpfV4af+tbDJL8aAIBfjt7H9J3XoXj763xSiQ8++ACjR49GaGiocltISAjGjh0LPz8/CZMRlT2v52p8jXM1EpUNLB6JiMowQRDwmV91fNejDgQB2HA2Ap9vvcypPNRg7dq1sLOzg4+PDwwMDGBgYICGDRvC1tYWq1evljoeUZnDuRqJyh6N6ezr5OSEiIi89+yMGzcOv/zyS57t69atw7Bhw3JtMzAwQEZGhtoyEhFpqo8bV4WZgS4+33oZ/1yKRmpmNn75yAv6uvyOUFWsra2xZ88e3LlzB7du3QIAuLm5oUaNGhInIyr7ijtX4+t714g0WXpWNmpP3w8AuDG7fbm4z1JjXsGFCxeQk5OjXL927Rratm2Lvn37FniMubk5bt++rVwXBKHAfYmIyrse9SvD1EAX4zaF4uCNOEzYFIpfBnpBT4cFpCo5OTlBFEW4uLhAV1djmlEiIiK105j/UVhbW8POzk65/Pfff3BxcUHLli0LPEYQhFzH2NralmJiIiLN41fbFqsH+0BfV4YDN+Lw2Z+X2IVVRdLT0zFixAgYGxvD3d0dkZGRAIBPP/0U8+fPlzgdERGR+mlM8fimrKwsbNy4EcOHD3/n1cTU1FRUrVoVjo6O6N69O65fv17gvgCQmZmJ5OTkXAsRUXnTooY1Vg7yhr6ODHuuxmLSX7wHUhUCAgJw+fJlHDt2DIaG/+ti5+fnhy1btkiYjIiIqHRoZPG4Y8cOJCYmYujQoQXuU7NmTaxduxY7d+7Exo0boVAo0KRJEzx69KjAYwIDAyGXy5WLo6OjGtITEUmvtZsNfh3oBT0dAf9efoyp269yGo/3tGPHDixbtgzNmjXL9cWmu7s77t+/L2EyIiKi0qGRxeOaNWvQsWNHVKpUqcB9fH19MXjwYHh6eqJly5bYvn07rK2tsXLlygKPCQgIQFJSknKJiopSR3wiIo3gV9sWSwd4QUcm4O+QR5i/75bUkcq0p0+fwsbGJs/2tLQ03nNPRERaQeOKx4iICBw6dAgjR44s1nF6enqoX78+7t27V+A+BgYGMDc3z7UQEZVnHerYYX4vDwDAyuMPsOoEr5CVlI+PD3bv3q1cf10wrl69Gr6+vlLFIiIiKjUaN0xcUFAQbGxs0Llz8YZfzsnJwdWrV9GpUyc1JSMiKpv6+jgiIS0LgXtvYd6eW6hoYpBrfjUqmnnz5qFjx464ceMGsrOz8fPPP+PGjRs4c+YMjh8/LnU8IiIitdOoK48KhQJBQUEYMmRInuHPBw8ejICAAOX67NmzceDAATx48AChoaEYNGgQIiIiin3FkohIG4xu6YJRLaoBAL7cdgUn7z6VOFHZ06xZM4SFhSE7OxseHh44cOAAbGxsEBwcDG9vb7U979y5c9GkSRMYGxujQoUK+e4TGRmJzp07w9jYGDY2Nvjiiy+QnZ2ttkxERKSdNOrK46FDhxAZGYnhw4fneSwyMhIy2f9q3efPn+OTTz5BbGwsLCws4O3tjTNnzqB27dqlGZmIqMyY2sENT1My8c+laIzbGIrt45qguq2Z1LHKFBcXF/z222+l+pxZWVno27cvfH19sWbNmjyP5+TkoHPnzrCzs8OZM2cQExODwYMHQ09PD/PmzSvVrETFYayvi4fzi9fTjIikJYhaPPxecnIy5HI5kpKSeP8jEWmFzOwcfLz6PM4/TICDhRF2jG8KK1MDqWNJorhtgJ+fHwYNGoRevXpJ0masW7cO/v7+SExMzLV979696NKlCx4/fqyc73jFihX46quv8PTpU+jr6xfp/GwTiYhUKz0rG7Wn7wcAHPm8JapZm0qcqGBFbQM0qtsqERGpl4GuDlZ+7A2nisZ49PwFPll/ERkvc6SOVSa4u7sjICAAdnZ26Nu3L3bu3ImXL19KHQvBwcHw8PBQFo4A0L59eyQnJ79z/mPOfUxEpF7bQv43haDfwuPYciFSwjSqweKRiEjLWJjoY83QBpAb6eFSZCK+/PsK54Asgp9//hnR0dHYsWMHTExMMHjwYNja2mLUqFGSDpgTGxubq3AEoFyPjY0t8DjOfUxEpD4xSS8wY9f/vsBTiMDX268hJumFhKneH4tHIiIt5GJtihWDvKErE7Dr8mOsPhkudaQyQSaToV27dli3bh3i4uKwcuVKnD9/Hm3atCnWeaZOnQpBEN653Lql3nk5OfcxEZH6hMenQfHW97I5ooiH8enSBFIRjRowh4iISo+vS0VM61IbM3ZdR+Dem6hlb45m1a2kjlUmxMbG4s8//8TGjRtx5coVNGzYsFjHf/755xg6dOg796lWrVqRzmVnZ4fz58/n2hYXF6d8rCAGBgYwMNDO+12JiNTN2coEMgG5CkgdQYCTlbF0oVSAVx6JiLTYYN+q6OPtAIUITNgciqiEsv2NqDolJycjKCgIbdu2haOjI5YvX45u3brh7t27OHv2bLHOZW1tDTc3t3cuRR3oxtfXF1evXsWTJ0+U2w4ePAhzc3OOQE5EJBF7uRFmdXNXrssEYF6vOrCXG0mY6v3xyiMRkRYTBAFzetTB3bgUXH6UhE/WX8Q/45rCSF9H6mgax9bWFhYWFujXrx8CAwPh4+NTKs8bGRmJhIQEREZGIicnB2FhYQAAV1dXmJqaol27dqhduzY+/vhjfP/994iNjcW3336L8ePH88oiEZGEens7YNrOV/c9Hpqs2aOtFhWvPBIRaTlDPR2s+NgbVqb6uBWbgmk7r0kdSeOIooglS5bg3r17WLRoUakVjgAwffp01K9fHzNmzEBqairq16+P+vXr4+LFiwAAHR0d/Pfff9DR0YGvry8GDRqEwYMHY/bs2aWWkYiI3s1Obih1BJXgPI+c04qICABw5n48Bq0+B4UIfN+nLj70Kd+jbxanDVAoFDA0NMT169dRvXr1UkpYutgmEhGp1pvzPN6Y3R7G+prb6ZPzPBIRUbE0cbHC5LY1AADTdlzDzRjO+/eaTCZD9erV8ezZM6mjEBERSYbFIxERKY1r5YoWNayRma3A+D9CkZqZLXUkjTF//nx88cUXuHaN3XqJiEg7sXgkIiIlmUzA4n6esDM3xIP4NHy9/Sq0+O6GXAYPHozz58+jXr16MDIygqWlZa6FiIiovNPcjrdERCQJSxN9/DKwPj5ceRa7Lj9GyxrW6O3tIHUsyS1evFjqCERERJJi8UhERHl4V7WE/wfV8dPBO5i+8xq8q1rAycpE6liSGjJkiNQRiIiIJMVuq0RElK9xrV3R0MkSaVk5+GxLGF7mKKSOJLn79+/j22+/xYABA/DkyRMAwN69e3H9+nWJkxEREakfi0ciIsqXjkzAov6eMDfUxeWoRCw+dEfqSJI6fvw4PDw8cO7cOWzfvh2pqakAgMuXL2PGjBkSpyMiIk1jrK+Lh/M74+H8zho9TUdxsHgkIqICVa5ghPm96wIAfj12H+ceaO9UFVOnTsWcOXNw8OBB6OvrK7e3adMGZ8+elTCZZkjPyobT1N1wmrob6VkcpZeIqDxi8UhERO/UycMeH/o4QBSBz7de1trpO65evYqePXvm2W5jY4P4+HgJEhEREZUuFo9ERFSoaV1qw8HCCI+ev8Dc3TekjiOJChUqICYmJs/2S5cuoXLlyhIkIiIiKl0sHomIqFBmhnr4sW89CAKw+XwUjtyKkzpSqevfvz+++uorxMbGQhAEKBQKnD59GlOmTMHgwYOljkdERKR2LB6JiKhIGleriOFNnQEAX227iudpWRInKl3z5s2Dm5sbHB0dkZqaitq1a6NFixZo0qQJvv32W6njERERqR2LRyIiKrIv2teEq40pnqZkYtrOa1LHKVX6+vr47bff8ODBA/z333/YuHEjbt26hQ0bNkBHR0fqeERERGrH4pGIiIrMUE8HCz+sBx2ZgP+uxGDftVipI5U6R0dHdOrUCb1790ZaWhqeP38udSQiIqJSweKRiIiKpa5DBYxuUQ0AMG3nNSSma0f3VX9/f6xZswYAkJOTg5YtW8LLywuOjo44duyYtOGIiIhKAYtHIiIqtokfVIeLtQmepmTiu/9uSh2nVPz999+oV68eAODff//FgwcPcOvWLUyaNAnffPONxOmIiIjUj8UjEREVm6GeDr7v82r01W2hj3D09hOpI6ldfHw87OzsAAB79uzBhx9+iBo1amD48OG4evWqxOk0S2xShtQRiIhIDVg8EhFRiXhXtVCOvvr19qtIyXgpcSL1srW1xY0bN5CTk4N9+/ahbdu2AID09HQOmANgW8gj5c9+C49jy4VICdMQEZE6sHgkIqISm9KuJqpWNEZMUgZ+2H9b6jhqNWzYMHz44YeoU6cOBEGAn58fAODcuXNwc3OTOJ20YpJeYMau68p1hQh8vf0aYpJeSJiKiIhUjcUjERGVmJG+Dub19AAAbDgbgdDI8jvy6MyZM7F69WqMGjUKp0+fhoGBAQBAR0cHU6dOlTidtMLj06AQc2/LEUU8jE+XJhAREamFrtQBiIiobGvqaoVeXpWxPTQaX2+/in8/bQY9nfL53WSfPn3ybBsyZIgESTSLs5UJZAJyFZA6ggAnK2PpQhERkcqVz9adiIhK1beda8PCWA+3YlOw+mS41HHU5vDhw+jSpQtcXFzg4uKCLl264NChQ1LHkpy93Aizurkr12UCMK9XHdjLjSRMRUREqsbikYiI3puliT6+7VwbAPDz4TuIfFb+uiv++uuv6NChA8zMzPDZZ5/hs88+g7m5OTp16oRffvlF6niS6+3toPz50OSW6NegioRpiIhIHVg8EhGRSvTyqoymrhWR8VKBb3degyiKhR9UhsybNw+LFi3C5s2bMXHiREycOBGbNm3CokWLMG/ePKnjaRQ7uaHUEYiISA1YPBIRkUoIgoA5PTygryvDiTtPsfdarNSRVCoxMREdOnTIs71du3ZISkqSIBEREVHpYvFIREQq42xlgrEtXQAAs/+9gdTMbIkTqU63bt3wzz//5Nm+c+dOdOnSRYJEREREpYujrRIRkUqNbeWCfy5FIzIhHUsO38XXnWpJHanElixZovy5du3amDt3Lo4dOwZfX18AwNmzZ3H69Gl8/vnnUkUkIiIqNRpz5dHJyQmCIORZxo8fX+AxW7duhZubGwwNDeHh4YE9e/aUYmIiIsqPoZ6OcuTNtafCcTs2ReJEJbdo0SLlsmbNGlhYWODGjRtYs2YN1qxZg+vXr6NChQpYu3at1FGJiIjUTmOuPF64cAE5OTnK9WvXrqFt27bo27dvvvufOXMGAwYMQGBgILp06YJNmzahR48eCA0NRZ06dUorNhER5aO1mw3a1bbFgRtxmLbzGraMagxBEKSOVWzh4eV32hEiIqLi0pgrj9bW1rCzs1Mu//33H1xcXNCyZct89//555/RoUMHfPHFF6hVqxa+++47eHl5YdmyZaWcnIiI8jO9a20Y6slwPjwBO8KipY6jMvHx8YiPj5c6BhERUanTmOLxTVlZWdi4cSOGDx9e4DfVwcHB8PPzy7Wtffv2CA4OLvC8mZmZSE5OzrUQEZF6OFgY49M21QEAgXtulenBcxITEzF+/HhYWVnB1tYWtra2sLKywoQJE5CYmCh1PCIiolKhMd1W37Rjxw4kJiZi6NChBe4TGxsLW1vbXNtsbW0RG1vw0PCBgYGYNWuWqmISEVEhRjZ3xl8XoxDxLB3LjtzD1I5uUkcqtoSEBPj6+iI6OhoDBw5ErVqvBgC6ceMG1q1bh8OHD+PMmTOwsLCQOCkREZF6aeSVxzVr1qBjx46oVKmSSs8bEBCApKQk5RIVFaXS8xMRUW4GujqY3qU2AGDNqQcIj0+TOFHxzZ49G/r6+rh//z5WrlwJf39/+Pv7Y9WqVbh37x709PQwe/ZsqWNKzlhfFw/nd8bD+Z1hrK+R300TEdF70rjiMSIiAocOHcLIkSPfuZ+dnR3i4uJybYuLi4OdnV2BxxgYGMDc3DzXQkRE6tXGzQatalrjZY6I7/67IXWcYtuxYwd+/PHHPL1dgFdt0ffff5/v/I9ERETljcYVj0FBQbCxsUHnzp3fuZ+vry8OHz6ca9vBgweVc28REZFmEAQB07rUhp6OgCO3nuDIrbjCD9IgMTExcHd3L/DxOnXqvPOWCSIiovJCo4pHhUKBoKAgDBkyBLq6ubu8DB48GAEBAcr1zz77DPv27cNPP/2EW7duYebMmbh48SImTJhQ2rGJiKgQLtamGN7UGQDw3X83kZWtkDhR0VlZWeHhw4cFPh4eHg5LS8vSC0RERCQRjSoeDx06hMjISAwfPjzPY5GRkYiJiVGuN2nSBJs2bcKqVatQr149/P3339ixYwfneCQi0lAT2rjCytQA4fFpWB/8UOo4Rda+fXt88803yMrKyvNYZmYmpk2bhg4dOkiQjIiIqHQJoiiKUoeQSnJyMuRyOZKSknj/IxFRKdhyIRJfbbsKM0NdHP+iNSxN9CXLUtQ24NGjR/Dx8YGBgQHGjx8PNzc3iKKImzdv4tdff0VmZiYuXrwIR0fHUkyvemwTiYi0V1HbAA6HRkREpaaPtyN+PxOBGzHJWHzoDmZ31/zeIg4ODggODsa4ceMQEBCA19+5CoKAtm3bYtmyZWW+cCQiIioKFo9ERFRqdGQCvu1SCx/9dg5/nIvEx42rorqtmdSxCuXs7Iy9e/fi+fPnuHv3LgDA1dWV9zoSEZFW0ah7HomIqPxr4mKFtrVtkaMQMXfPTanjFIuFhQUaNmyIhg0bsnAkIiKtw+KRiIhK3dedakFPR8Cx209x7PYTqeMQERFREbB4JCKiUudsZYIhvk4AgMA9t5Cj0Nqx24iIiMoMFo9ERCSJT9tUh9xID7fjUvB3SJTUcYiIiKgQLB6JiEgScmM9fNrGFQDw04E7SM/KljgRkWqlZ2XDaepuOE3dzc83EZULLB6JiEgyH/tWhaOlEZ6kZGL1yXCp4xAREdE7sHgkIiLJGOjq4Mv2bgCAFcfv40lKhsSJiIiIqCAsHomISFJd6trD07EC0rNysPjQXanjEBERUQFYPBIRkaQEQcA3nWsBAP48H4l7T1IlTkRERET5YfFIRESSa+Bkiba1baEQgR/335Y6DhEREeWDxSMREWmEL9vXhEwA9l2PRWjkc6njEBER0VtYPBIRkUaobmuGPt4OAID5e25BFEWJExEREdGbWDwSEZHG8PerAQNdGc4/TMDR20+kjiO5hw8fYsSIEXB2doaRkRFcXFwwY8YMZGVl5drvypUraN68OQwNDeHo6Ijvv/9eosRERFSesXgkIiKNUamCEYY2dQIAfL/vNnIU2n318datW1AoFFi5ciWuX7+ORYsWYcWKFfj666+V+yQnJ6Ndu3aoWrUqQkJC8MMPP2DmzJlYtWqVhMnpbbFJnIaGiMo+Fo9ERKRRxrV0hbmhLm7FpuCfS9FSx5FUhw4dEBQUhHbt2qFatWro1q0bpkyZgu3btyv3+eOPP5CVlYW1a9fC3d0d/fv3x8SJE7Fw4UIJkxMAbAt5pPzZb+FxbLkQKWEaIqL3x+KRiIg0itxYD+NauwIAFh28g8zsHIkTaZakpCRYWloq14ODg9GiRQvo6+srt7Vv3x63b9/G8+cFDzyUmZmJ5OTkXAupTkzSC8zYdV25rhCBr7dfQ0zSCwlTERG9HxaPRESkcYb4OsHGzADRiS+w5UKU1HE0xr1797B06VKMHj1auS02Nha2tra59nu9HhsbW+C5AgMDIZfLlYujo6N6Qmup8Pg0vN3rOkcU8TA+XZpAREQqwOKRiIg0jpG+Dj79oDoAYMnhe0jPypY4kWpNnToVgiC8c7l161auY6Kjo9GhQwf07dsXn3zyyXtnCAgIQFJSknKJimKRrkrOViaQCbm36QgCnKyMpQlERKQCulIHICIiyk8/H0f8duIBIhPSse7MQ4xr5Sp1JJX5/PPPMXTo0HfuU61aNeXPjx8/RuvWrdGkSZM8A+HY2dkhLi4u17bX63Z2dgWe38DAAAYGBsVMTkVlLzfCrG7umLbzVddVmQDM61UH9nIjiZMREZUci0ciItJI+royTGpbHZO2XMaKY/cxsGFVyI31pI6lEtbW1rC2ti7SvtHR0WjdujW8vb0RFBQEmSx3pyFfX1988803ePnyJfT0Xr0/Bw8eRM2aNWFhYaHy7FR0vb0dlMXjocktUc3aVOJERETvh91WiYhIY3WrVxk1bE2RnJGNVSfvSx2n1EVHR6NVq1aoUqUKfvzxRzx9+hSxsbG57mX86KOPoK+vjxEjRuD69evYsmULfv75Z0yePFnC5PQ2O7mh1BGIiN4bi0ciItJYOjIBU9rVBACsPfUQT1MyJU5Uug4ePIh79+7h8OHDcHBwgL29vXJ5TS6X48CBAwgPD4e3tzc+//xzTJ8+HaNGjZIwORERlUcsHomISKO1rW2Leo4V8OJlDpYf066rj0OHDoUoivkub6pbty5OnjyJjIwMPHr0CF999ZVEiYmIqDxj8UhERBpNEAR83rYGAGDjuQjEJmVInIiIiEg7sXgkIiKN17y6FRo6WSIrW4Ffjt6TOg4REZFWYvFIREQaTxAETG736urjnxci8eg5J1onIiIqbSweiYioTGhcrSKauVrhZY6IpYd59ZGIiKi0sXgkIqIy4/XVx79DH+FhfJrEaYiIiLQLi0ciIiozvKpYoHVNa+QoRPx8+K7UcYiIiLSKrtQBiIiIimNy25oIjUxEdVtTiKIIQRCkjkSUL2N9XTyc31nqGEREKsPikYiIyhQPBznOBnwAI30dqaMQERFpFXZbJSKiMoeFIxERUenTmOIxOjoagwYNQsWKFWFkZAQPDw9cvHixwP2PHTsGQRDyLLGxsaWYmoiIiIiISDtoRLfV58+fo2nTpmjdujX27t0La2tr3L17FxYWFoUee/v2bZibmyvXbWxs1BmViIiIiIhIK2lE8bhgwQI4OjoiKChIuc3Z2blIx9rY2KBChQpqSkZERERERESAhnRb3bVrF3x8fNC3b1/Y2Nigfv36+O2334p0rKenJ+zt7dG2bVucPn36nftmZmYiOTk510JERERERESF04ji8cGDB1i+fDmqV6+O/fv3Y+zYsZg4cSJ+//33Ao+xt7fHihUrsG3bNmzbtg2Ojo5o1aoVQkNDCzwmMDAQcrlcuTg6Oqrj5RAREREREZU7giiKotQh9PX14ePjgzNnzii3TZw4ERcuXEBwcHCRz9OyZUtUqVIFGzZsyPfxzMxMZGZmKteTk5Ph6OiIpKSkXPdNEhFR+ZecnAy5XM424P/x/SAi0l5FbQM04sqjvb09ateunWtbrVq1EBkZWazzNGzYEPfu3SvwcQMDA5ibm+daiIiIiIiIqHAaUTw2bdoUt2/fzrXtzp07qFq1arHOExYWBnt7e1VGIyIiIiIiImjIaKuTJk1CkyZNMG/ePHz44Yc4f/48Vq1ahVWrVin3CQgIQHR0NNavXw8AWLx4MZydneHu7o6MjAysXr0aR44cwYEDB6R6GUREREREROWWRhSPDRo0wD///IOAgADMnj0bzs7OWLx4MQYOHKjcJyYmJlc31qysLHz++eeIjo6GsbEx6tati0OHDqF169ZFft7Xt3ty1FUiIu3z+ne/Btz6rxHYJhIRaa+itokaMWCOVB49esQRV4mItFxUVBQcHBykjiE5tolERFRYm6jVxaNCocDjx49hZmYGQRAAvLoKeuHChXz3L+ix16O2RkVFafwgPO96fZr2HCU5T3GPKcr+he3Dz4zmPEdJz1Oc49T9mSno8bL0mQHKxudGFEWkpKSgUqVKkMk0YggASeXXJmqK0vg8lfbzqvLc73uu0vjdWZz9y9rvOylI9W+iuPhvV73nUuW/3aK2iRrRbVUqMpksT2Wto6NT4C+qdz0GoEyM4FrYa9Ck5yjJeYp7TFH2L2wffmY05zlKep7iHKfuz0xhj5eFzwxQdj43crlcRWnKvvzaRE1RGp+n0n5eVZ77fc9VGr87S7J/Wfl9JwWp/k0UF//tqvdcqv63W5Q2kV+1vmX8+PEleqysKI3XoKrnKMl5intMUfYvbB9+ZjTnOUp6nuIcp+7PTHHzaKqy9LkhzSfV37U6n1eV537fc5XG7873eR7Kq6y8l/y3q95zlda/3TdpdbdVVeHEylRc/MxQcfEzQ0Tagr/viDQXrzyqgIGBAWbMmAEDAwOpo1AZwc8MFRc/M0SkLfj7jkhz8cojERERERERFYpXHomIiIiIiKhQLB6JiIiIiIioUCweiYiIiIiIqFAsHomIiIiIiKhQLB6JiIiIiIioUCwe1axnz56wsLBAnz59pI5CZUBUVBRatWqF2rVro27duti6davUkagMSExMhI+PDzw9PVGnTh389ttvUkciIlI5tpFE0uNUHWp27NgxpKSk4Pfff8fff/8tdRzScDExMYiLi4OnpydiY2Ph7e2NO3fuwMTEROpopMFycnKQmZkJY2NjpKWloU6dOrh48SIqVqwodTQiIpVhG0kkPV55VLNWrVrBzMxM6hhURtjb28PT0xMAYGdnBysrKyQkJEgbijSejo4OjI2NAQCZmZkQRRH8XpCIyhu2kUTSY/H4DidOnEDXrl1RqVIlCIKAHTt25Nnnl19+gZOTEwwNDdGoUSOcP3++9IOSxlDlZyYkJAQ5OTlwdHRUc2qSmio+N4mJiahXrx4cHBzwxRdfwMrKqpTSExEVDdtIorKPxeM7pKWloV69evjll1/yfXzLli2YPHkyZsyYgdDQUNSrVw/t27fHkydPSjkpaQpVfWYSEhIwePBgrFq1qjRik8RU8bmpUKECLl++jPDwcGzatAlxcXGlFZ+IqEjYRhKVAyIVCQDxn3/+ybWtYcOG4vjx45XrOTk5YqVKlcTAwMBc+x09elTs3bt3acQkDVLSz0xGRobYvHlzcf369aUVlTTI+/yueW3s2LHi1q1b1RmTiOi9sI0kKpt45bGEsrKyEBISAj8/P+U2mUwGPz8/BAcHS5iMNFVRPjOiKGLo0KFo06YNPv74Y6mikgYpyucmLi4OKSkpAICkpCScOHECNWvWlCQvEVFJsI0kKhtYPJZQfHw8cnJyYGtrm2u7ra0tYmNjlet+fn7o27cv9uzZAwcHBxaWWqwon5nTp09jy5Yt2LFjBzw9PeHp6YmrV69KEZc0RFE+NxEREWjevDnq1auH5s2b49NPP4WHh4cUcYmISoRtJFHZoCt1gPLu0KFDUkegMqRZs2ZQKBRSx6AypmHDhggLC5M6BhGRWrGNJJIerzyWkJWVFXR0dPIMShEXFwc7OzuJUpEm42eGSoKfGyLSBvxdR1Q2sHgsIX19fXh7e+Pw4cPKbQqFAocPH4avr6+EyUhT8TNDJcHPDRFpA/6uIyob2G31HVJTU3Hv3j3lenh4OMLCwmBpaYkqVapg8uTJGDJkCHx8fNCwYUMsXrwYaWlpGDZsmISpSUr8zFBJ8HNDRNqAv+uIygGph3vVZEePHhUB5FmGDBmi3Gfp0qVilSpVRH19fbFhw4bi2bNnpQtMkuNnhkqCnxsi0gb8XUdU9gmiKIqlW64SERERERFRWcN7HomIiIiIiKhQLB6JiIiIiIioUCweiYiIiIiIqFAsHomIiIiIiKhQLB6JiIiIiIioUCweiYiIiIiIqFAsHomIiIiIiKhQLB6JiIiIiIioUCweiYiIiIiIqFAsHolKaObMmfD09JQ6Rr5KO5uTkxMWL16sXBcEATt27FDr8wmCAEEQkJiYqLLztmrVSnnesLAwlZ2XiKi8Y5v4P2wTqTxj8Ujl3tChQ5W//PT19eHq6orZs2cjOzv7vc47ZcoUHD58WEUpNbvhLa6YmBh07NhRrc8xe/ZsxMTEQC6Xq+yc27dvx/nz51V2PiIiTcM2sfSxTaTyRFfqAESloUOHDggKCkJmZib27NmD8ePHQ09PDwEBAXn2zcrKgr6+fqHnNDU1hampqTriljpRFJGTkwNdXdX8SrCzs1PJed7FzMxM5c9jaWmJ5ORklZ6TiEjTsE18N7aJr7BNpPzwyiNpBQMDA9jZ2aFq1aoYO3Ys/Pz8sGvXLgCvvoXt0aMH5s6di0qVKqFmzZoAgKtXr6JNmzYwMjJCxYoVMWrUKKSmpirPmd+3oqtXr0atWrVgaGgINzc3/Prrr7kef/ToEQYMGABLS0uYmJjAx8cH586dw7p16zBr1ixcvnxZ+Y3wunXrAACJiYkYOXIkrK2tYW5ujjZt2uDy5cu5zjt//nzY2trCzMwMI0aMQEZGxjvfj2PHjkEQBOzduxfe3t4wMDDAqVOncP/+fXTv3h22trYwNTVFgwYNcOjQoVzHPnnyBF27doWRkRGcnZ3xxx9/5Dn/m110Xj/Xm11pwsLCIAgCHj58CACIiIhA165dYWFhARMTE7i7u2PPnj3vfA1vW7duHSpUqIAdO3agevXqMDQ0RPv27REVFZVrv3///RcNGjSAoaEhrKys0LNnz2I9DxFRWcc2MTe2iWwTqeh45ZG0kpGREZ49e6ZcP3z4MMzNzXHw4EEAQFpaGtq3bw9fX19cuHABT548wciRIzFhwgRlA/a2P/74A9OnT8eyZctQv359XLp0CZ988glMTEwwZMgQpKamomXLlqhcuTJ27doFOzs7hIaGQqFQoF+/frh27Rr27dunbJhedz3p27cvjIyMsHfvXsjlcqxcuRIffPAB7ty5A0tLS/z111+YOXMmfvnlFzRr1gwbNmzAkiVLUK1atULfh6lTp+LHH39EtWrVYGFhgaioKHTq1Alz586FgYEB1q9fj65du+L27duoUqUKgFf/sXj8+DGOHj0KPT09TJw4EU+ePHmfvw6MHz8eWVlZOHHiBExMTHDjxo0SfYOdnp6OuXPnYv369dDX18e4cePQv39/nD59GgCwe/du9OzZE9988w3Wr1+PrKysYjfIRETlDdvEV9gmsk2kIhCJyrkhQ4aI3bt3F0VRFBUKhXjw4EHRwMBAnDJlivJxW1tbMTMzU3nMqlWrRAsLCzE1NVW5bffu3aJMJhNjY2NFURTFGTNmiPXq1VM+7uLiIm7atCnXc3/33Xeir6+vKIqiuHLlStHMzEx89uxZvjnfPp8oiuLJkydFc3NzMSMjI9d2FxcXceXKlaIoiqKvr684bty4XI83atQoz7nedPToURGAuGPHjgL3ec3d3V1cunSpKIqiePv2bRGAeP78eeXjN2/eFAGIixYtUm4DIP7zzz+5nuv58+fKxy9duiQCEMPDw0VRFEUPDw9x5syZhWZ5rWrVqrmeTxRFMSgoSAQgnj17Nk+2c+fOiaL46r0aOHDgO88dHh4uAhAvXbpU5DxERGUF28S82CYWjG0ivY3dVkkr/PfffzA1NYWhoSE6duyIfv36YebMmcrHPTw8ct3TcfPmTdSrVw8mJibKbU2bNoVCocDt27fznD8tLQ3379/HiBEjlPd9mJqaYs6cObh//z6AV91S6tevD0tLyyLnvnz5MlJTU1GxYsVc5w0PD1ee9+bNm2jUqFGu43x9fYt0fh8fn1zrqampmDJlCmrVqoUKFSrA1NQUN2/eRGRkpPK5dHV14e3trTzGzc0NFSpUKPJrys/EiRMxZ84cNG3aFDNmzMCVK1dKdB5dXV00aNAgT7abN28CePV38MEHH7xXViKiso5tYv7YJhIVjt1WSSu0bt0ay5cvh76+PipVqpTnJvg3G8SSeH3fx2+//Zan0dLR0QHwqltQSc5rb2+PY8eO5XnsfRsnIO/rnjJlCg4ePIgff/wRrq6uMDIyQp8+fZCVlVXi55DJXn1HJYqictvLly9z7TNy5Ei0b98eu3fvxoEDBxAYGIiffvoJn376aYmfNz8l+TsgIipv2Cbmj20iUeF45ZG0gomJCVxdXVGlSpUijZ5Wq1YtXL58GWlpacptp0+fhkwmUw4e8CZbW1tUqlQJDx48gKura67F2dkZAFC3bl2EhYUhISEh3+fU19dHTk5Orm1eXl6IjY2Frq5unvNaWVkps547dy7XcWfPni30Nebn9OnTGDp0KHr27AkPDw/Y2dkpb+AHXn1rmZ2djZCQEOW227dvv3NeKWtrawCvhip/Lb/5ohwdHTFmzBhs374dn3/+OX777bdi58/OzsbFixfzZKtVqxaAV38HqhxKnoioLGKbWDRsE4nyYvFIlI+BAwfC0NAQQ4YMwbVr13D06FF8+umn+Pjjj2Fra5vvMbNmzUJgYCCWLFmCO3fu4OrVqwgKCsLChQsBAAMGDICdnR169OiB06dP48GDB9i2bRuCg4MBvJrkNzw8HGFhYYiPj0dmZib8/Pzg6+uLHj164MCBA3j48CHOnDmDb775RtkgfPbZZ1i7di2CgoJw584dzJgxA9evXy/R665evTq2b9+OsLAwXL58GR999BEUCoXy8Zo1a6JDhw4YPXo0zp07h5CQEIwcOfKd3166urrC0dERM2fOxN27d7F792789NNPufbx9/fH/v37ER4ejtDQUBw9elTZuBWHnp4ePv30U2W2oUOHonHjxmjYsCEAYMaMGdi8eTNmzJiBmzdv4urVq1iwYEGxn4eISJuwTWSbSPQai0eifBgbG2P//v1ISEhAgwYN0KdPH3zwwQdYtmxZgceMHDkSq1evRlBQEDw8PNCyZUusW7dO+S2rvr4+Dhw4ABsbG3Tq1AkeHh6YP3++sgtP79690aFDB7Ru3RrW1tbYvHkzBEHAnj170KJFCwwbNgw1atRA//79ERERoWyw+/Xrh2nTpuHLL7+Et7c3IiIiMHbs2BK97oULF8LCwgJNmjRB165d0b59e3h5eeXaJygoCJUqVULLli3Rq1cvjBo1CjY2NgWeU09PD5s3b8atW7dQt25dLFiwAHPmzMm1T05ODsaPH49atWqhQ4cOqFGjRp4h3YvC2NgYX331FT766CM0bdoUpqam2LJli/LxVq1aYevWrdi1axc8PT3Rpk0bToBMRFQItolsE4leE8Q3O10TUZEFBATg5MmTOHXqlNRRtI6TkxP8/f3h7++v3LZu3Tr4+/u/s7tQUTx8+BDOzs64dOlSnjnLiIgof2wTpcM2kUoTrzwSFZMoirh//z4OHz4Md3d3qeNora+++gqmpqZISkpS2Tk7duzIv1MiomJgm6gZ2CZSaeFoq0TFlJSUhNq1a6NBgwb4+uuvpY6jlY4fP64cnc7MzExl5129ejVevHgBAMoJoImIqGBsE6XHNpFKE7utEhERERERUaHYbZWIiIiIiIgKxeKRiIiIiIiICsXikYiIiIiIiArF4pGIiIiIiIgKxeKRiIiIiIiICsXikYiIiIiIiArF4pGIiIiIiIgKxeKRiIiIiIiICvV/Pv7vB8jnkWcAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 900x350 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rng = np.random.default_rng(123)\n",
    "u = rng.uniform(size=32)\n",
    "R_pc = 200. * np.sqrt(u / (1. - u))  # Projected Plummer radii.\n",
    "e_vlos_kms = np.full(32, 2.)\n",
    "variance = np.asarray(model.sigmalos2(R_pc, params=truth, **options))\n",
    "vlos_kms = rng.normal(0., np.sqrt(variance + e_vlos_kms**2))\n",
    "data = dict(R_pc=R_pc, vlos_kms=vlos_kms, e_vlos_kms=e_vlos_kms)\n",
    "\n",
    "R_grid = np.geomspace(10., 3000., 100)\n",
    "sigma_kms = np.sqrt(model.sigmalos2(R_grid, params=truth, **options))\n",
    "prediction, axes = plt.subplots(1, 2, figsize=(9, 3.5), layout=\"constrained\")\n",
    "axes[0].semilogx(R_grid, sigma_kms)\n",
    "axes[0].set(xlabel=\"Projected radius [pc]\", ylabel=r\"$\\sigma_{los}$ [km/s]\")\n",
    "axes[1].errorbar(R_pc, vlos_kms, yerr=e_vlos_kms, fmt=\".\")\n",
    "axes[1].set(xscale=\"log\", xlabel=\"Projected radius [pc]\",\n",
    "            ylabel=\"Observed LOS velocity [km/s]\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0ba74b9",
   "metadata": {},
   "source": [
    "The setup above supplies `fixed`, `truth` and `options` for this mock.\n",
    "The inverse Plummer CDF gives projected radii; real surveys may need spatial\n",
    "selection, membership and measurement-error models. Generating Gaussian\n",
    "velocities tests the assumed likelihood workflow, not a self-consistent\n",
    "phase-space distribution. For an observed dataset the velocity likelihood\n",
    "conditions on the supplied radii, so it does not also fit their distribution.\n",
    "\n",
    "Next: [turn predictions into an MCMC likelihood](inference.ipynb)."
   ]
  }
 ],
 "metadata": {
  "jeanspy": {
   "execution": {
    "code_sha256": "f66a1541e8f03375ff485a36d2214da26d44d575eb9e947bceba61a7bc11938b",
    "jax_enable_x64": true,
    "lock_sha256": "a30f366620e23a49d6d7466055223a761824454169d8fc0364c24ac25b626dc7",
    "packages": {
     "arviz": "1.0.0",
     "corner": "2.2.3",
     "jax": "0.9.1",
     "jeanspy": "0.1.0",
     "matplotlib": "3.10.8",
     "numpy": "2.4.3",
     "numpyro": "0.20.0",
     "scipy": "1.17.1"
    },
    "platform": "cpu",
    "source_sha256": {
     "src/jeanspy/__init__.py": "61671d182af04c0056cfb8d273db499d8774694711506698dbe7dc4a63d3b403",
     "src/jeanspy/_axisymmetric_components.py": "4b203c616c87e55182d105b37f9fa49e44fe1df6bc430ccccb09edc5cca42530",
     "src/jeanspy/_axisymmetric_params.py": "85da8de45bf5ba8cd9a40a137f7ccf971f7ce9534f312ef49fcd33539a230802",
     "src/jeanspy/_jax_env.py": "2ef46eaf5cbd8e72a0c17351490207c29586979082cabc87128c9bd16fc689a4",
     "src/jeanspy/_numpy/__init__.py": "de2d19d61e00702f13996582369f6acbe92b74628e6f569f91c3215a53580117",
     "src/jeanspy/_numpy/core.py": "d073c6409d5f8cf527947c1ddf3916cb39cef3f6d4e43173cad59d148e463296",
     "src/jeanspy/_numpy/inference.py": "ab58728400b6e56189ffba9ee805d9edbf59ca13dfebe05498801592314552db",
     "src/jeanspy/_numpy/jfactor.py": "823f58a66a615ba59e057022ce51e2adfab577c25a4ced2f2033304c2cbd6cfe",
     "src/jeanspy/_numpy/profiles.py": "e77fda6fb068392a63f0479e4c33bf80a22de3d2ab2cd101ea0c1f33e579491d",
     "src/jeanspy/_numpy/solver.py": "b419796e2a912a46f65404cb5026ccb73799e1dbb8f80052f404148dd1f5ef3c",
     "src/jeanspy/_sampling_identity.py": "b0e70b77e38f25e2148ab1f108a714512a8b6054405c846c3a41730b07081cec",
     "src/jeanspy/_sersic_deprojection.py": "004c5566448698ecfe536096c3fc8d5c348a8d4e1610042763c77568c0f3d821",
     "src/jeanspy/_zhao.py": "df83c7427dc583123c7464a03d2d45272e90f5a6724c9729eee64fad92e0c41a",
     "src/jeanspy/axisymmetric.py": "683d67b91f8eb67e8d898eb23e6c8ca0e5c78020a65a7d05e183c19b076bd506",
     "src/jeanspy/axisymmetric_factors.py": "754b256f7d3d43cc72571cde2b22bc9419917cf13b17342cd5f52406c4111336",
     "src/jeanspy/axisymmetric_inference.py": "57078226b3643f8babda4d942b96ae546622da28802119d3d99ff6a0c5b45695",
     "src/jeanspy/axisymmetric_jax.py": "de5f78937104ce7008242d55c321c8a3c6d9491f6aa25d4be3bdab8de076ed54",
     "src/jeanspy/baes_eta2.py": "c274a16ca8c4b4d8d76ec30a2e705b5a5f6ce7659b31e31b5c685f92ab7a7ece",
     "src/jeanspy/dequad.py": "6260e7ee2d8da12265573e86ff81ba0a3e8912df5fec8588a84bfda9d3b1133d",
     "src/jeanspy/hyp2f1_jax.py": "03d8df3b48fcd9d418e289ea2ffaf174ffbdf14dffe483a150403c092bbaeea8",
     "src/jeanspy/model.py": "72af25df8f34816f27d373fdecf41408c79eb505cb4cb52c830822d5d8f12f2b",
     "src/jeanspy/model_jax.py": "355ba42ad16b2cadd00906679dd06ee12fd85fe4f1fc48630c779e68e4fabd64",
     "src/jeanspy/parameters.py": "7760f4c85063105eb1328910867f3c57516f885dc577533bf9b91ac08a648cb2",
     "src/jeanspy/sampler.py": "42ca0bbcef84b7c1620fd977d71810284a5d6abd70ff6047ae3175822d63f10e",
     "src/jeanspy/sampler_numpyro.py": "ccd367d576bc55173d85939ad4bf7eeb8a44130d693f40ee56f4eade278c2195",
     "src/jeanspy/sersic.py": "487f473351f9a309ad7d8dad1b1405e496196413e7606d0b0a18c5f71bfcf4a3"
    },
    "utc": "2026-09-18T00:32:47.065310+00:00"
   },
   "mcmc": false
  },
  "kernelspec": {
   "display_name": "Python 3 (JeansPy)",
   "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.12.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}