telegraphic/speedup_frbpoppy.ipynb Secret

## speedup_frbpoppy.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "available-navigation",
   "metadata": {},
   "outputs": [],
   "source": [
    "from frbpoppy import precalc as pc\n",
    "import healpy as hp\n",
    "import numpy as np\n",
    "from astropy.coordinates import Angle"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "collaborative-pillow",
   "metadata": {},
   "source": [
    "### Basic profiling of CosmicPopulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "id": "increased-grass",
   "metadata": {},
   "outputs": [],
   "source": [
    "from frbpoppy import CosmicPopulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 150,
   "id": "precise-pontiac",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now generate FRB population\n",
    "n_source = 1000\n",
    "n_days   = 1.0\n",
    "cosmic_pop = CosmicPopulation.complex(n_source, n_days=n_days)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 153,
   "id": "realistic-andorra",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.8 µs ± 3.27 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)\n",
      "1.04 s ± 15.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "137 ns ± 1.52 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each)\n",
      "382 µs ± 3 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n",
      "125 µs ± 959 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n",
      "1.44 s ± 17.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "179 µs ± 86.9 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n",
      "75.8 µs ± 341 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n",
      "45.3 µs ± 280 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit cosmic_pop.gen_index()\n",
    "%timeit cosmic_pop.gen_dist()\n",
    "%timeit cosmic_pop.gen_time()\n",
    "%timeit cosmic_pop.gen_direction()\n",
    "%timeit cosmic_pop.gen_gal_coords()\n",
    "%timeit cosmic_pop.gen_dm()\n",
    "%timeit cosmic_pop.gen_w()\n",
    "%timeit cosmic_pop.gen_lum()\n",
    "%timeit cosmic_pop.gen_si()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "overall-johnson",
   "metadata": {},
   "source": [
    "Two slowest functions are `gen_dist()` and `gen_dm()`. Target these for optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "previous-constraint",
   "metadata": {},
   "source": [
    "## Speedup `gen_dm()`\n",
    "\n",
    "### Setup NE2001 lookup table\n",
    "\n",
    "First let's setup the current lookup system, which uses an sqlite table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 189,
   "id": "applicable-mineral",
   "metadata": {},
   "outputs": [],
   "source": [
    "ne2001 = pc.NE2001Table()\n",
    "\n",
    "def dm_func_sql(gl, gb):\n",
    "    z = lambda: ne2001.lookup(gl, gb)\n",
    "    return z()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "running-piece",
   "metadata": {},
   "source": [
    "### Setup healpix-based lookup\n",
    "\n",
    "In lieu of sqlite, let's try using healpix coversions, which can be done by healpy ufuncs.\n",
    "\n",
    "The `hp.ang2pix` will convert a sky coordinate into a healpix pixel ID, and `hp.pix2ang` does the opposite.\n",
    "\n",
    "Using a precomputed NE2001 map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 146,
   "id": "republican-legislature",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NSIDE = 128\n",
      "ORDERING = RING in fits file\n",
      "INDXSCHM = IMPLICIT\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlkAAAFzCAYAAAANJxyKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABlRUlEQVR4nO29fdD1zlnf9710zn0/z+/NNsY4ECBgmsAkJFPakpBOy0sTbJNkhqZQEzyGQEwbJpOkIWk6TFrahrTknSGdZjqEFie8FQoJtJOSYrudYCC8JECBEiChvBhMKMYG8/Pv93ue577P0fYP7UqXVrurlY50jqTz/czco/tIq9VqJe1+dV3XrsQYA0IIIYQQMi3FpQtACCGEELJFKLIIIYQQQmaAIosQQgghZAYosgghhBBCZoAiixBCCCFkBiiyCCGEEEJmgCKLEHISIvL5IvI96rcRkd86Yf5vEpG3JbZ/p4j8R1MdT+X7goh81NT5EkKuB4osQq4YEfl5EbkTkVd563/EiqWPvFDRaowx32CMed0FjvusMeZnz31cQsh2oMgihPwcgDe6HyLyuwA8dbniEELINqDIIoR8HYA/qn5/HoCv1QlE5OUi8rUi8qsi8k4R+RIRSbYfIvIaEXmfSyci/5OIvFtt/3oR+SKV/1eLyC+LyC+JyH8rIju7zXdHvlZEfkpEfkNE/jYA8Y77ZhH5SRH5dRF5q4h8RKR83yEif8pb96Mi8hn2/9rtKSIPRORvisgviMiviMhXishTdts7ROQz7f//rt3vD9rfnyoiP5KqJ0LIdqHIIoR8P4CXichvt8LmjwD4ei/Nfw/g5QA+CsAnoxJlfyyVqTHm5wA8D+DfsKs+EcALIvLb7e9PAvAO+//XADgA+K02/esAdOKsrFvzHwD4EgCvAvAzAP4dtf0PA/jPAXwGgA8C8N0AvjFSxP8ZbQve7wDwEQC+PZD2rwH4aAAfZ8v4oQD+K7vtHQA+RZ3Tz6KqI/8cCSFXBkUWIQRorFmvBfBTAH7JbVDC6y8YY95vjPl5AF8O4HMz8n0HgE8WkQ+2v/++/f0aAC8D8KMi8psA/AEAX2SMedEY824AXwHgswP5/UEAP2GM+fvGmHsAfwvA/6e2fyGAv2KM+UljzAHAXwbwcRFr1rd5294E4FuNMU90IhERAP8xgD9rjPk1Y8z7bb6ufO9AW1T9FfX7k0GRRcjVsr90AQghi+DrAHwXgNfAcxWishjdAninWvdOVNacPt4B4NMBvMvm/52oxNljAN9tjCmtyLkB8MuVngFQvQD+YiC/36zXG2OMiOh0HwHgvxORL1frxJZVlx/GmPeLyLejEkt/zS7/eOCYHwTgaQA/pMonAHb2/+8D8NFWLH6cPd8vtVa332PPmxByhVBkEUJgjHmniPwcKkvRF3ib3wPgHpWA+Qm77rdAWbsSvAPA30Alst4B4HsAfCUqkeUsPL8I4AmAV1nrU4pfBvDh7oe1Mn242v6LAL7MGPMNGWUDKlfify0i34Uq2P8fB9K8B8AjAB9rjOmcszHmJRH5IQB/BsCPG2PuROR7Afw5AD9jjHlPZlkIIRuD7kJCiOMLAPw+Y8yLeqUx5gjgmwF8mYg8Zy1Pfw7duK0OxpifRiVQPgfAdxljngfwKwA+E1ZkGWN+GcDbAHy5iLxMRAoR+ddE5JMDWX47gI8Vkc8QkT2A/wTAB6vtXwngL4jIxwJ1QP0bEkX8R6jE418C8L8YY8rAOZQA/kcAXyEir7b5fqiIvF4leweAP4VGOH6n95sQcoVQZBFCAADGmJ8xxvxgZPOfBvAiqqDu70EVNP6WzKzfAeC9xphfUL8FwP+t0vxRVC7JnwDw66hitz4kUMb3AHgDgL8K4L0AfhuAf6K2fxsq1983icjzAH4cVbxXEBt/9a0APtWeU4wvBvD/Avh+m+//CeBjvHN8Do1r0P9NCLlCxBhz6TIQQgghhGwOWrIIIYQQQmaAIosQQgghZAYosgghhBBCZoAiixBCCCFkBiiyCCGEEEJmIHcyUg5BJIQQQgjpIrENnPGdENLitUVq7k6Sw9vLb7l0EQghCyB3nixasghZORRPy4NijJBNELVkUWQRshEoorYHRRghq4Aii5C1QdFEhkJRRshFoMgiZKlQTJG5ofgiZFYosghZAhRUZClQeBEyGRRZhMwJxRPZKhRjhPRCkUXIlFBUkWuFoouQDhRZhJwCRRUhYSi6CKHIIiQJRRQh80ARRq4AiixCNBRVhFwGii6yQSiyyHVCMUXIOqD4IiuGIotcBxRVhGwDii6yIiiyyHahsCJk21BwkYVDkUXWD8UUIURD8UUWAkUWWScUVoSQHCi4yAWhyCLrgcKKEHIKFFzkzFBkkWVCQUUIOQcUXmRGKLLIcqCwIoRcEgouMjEUWeTyUFwRQpYExRaZCIoscl4oqAgha4TCi4yAIovMD4UVIWRLUHCRTCiyyPRQVBFCrgmKLhKBIotMB8UVIeSaodgiHhRZ5DQorAghpAsFFwFFFhkKRRUhhAyHousqocgieVBcEULI6VBsXRUUWSQNxRUhhEwPxdZVQJFF2lBUEULI+aHo2iQUWaSC4ooQQi4PxdamoMi6ViiqCCGbQmx/5vou/7e/TqS71OT1gbND0bVqKLKuEQosQsgiCYkgf5vG76dCaaYgJMJCx58JCq3VQpF1LVBYEULOQkgc9YkloGuBWjMzii8KrlVBkbVlKKwI2QAht1fOPiGXWCydPo4jZr3p20bmIeDafPvxmy9cKNIDRdYWee3us6Yzb/vm+yEsJKaBkNUTe5bHCB0KpM1BsbVYKLK2xmt3n3XpIoSh4CJkPBRFpAcKrUVCkbV2FiuqYui3aAovQvqhwCIjoOhaBNGHtzhnKcg4ViewgG6HEfrNToUQQk5ilf3DFUFL1oLZ7MPTF5hLyDXBlw0yEbRqXQy6C9fCZoVVLn2TB65BiFEwkiFQZJEZoOA6KxRZS+fqxVUuqWHqS2DJZSPLhCKLzAjF1lmgyFoiFFYz0Td0PWfCxFOPT0gKCityASi4ZoMia2lQYG0UCiySA0UWuRAUWrNAkbUUKK42DkUWiUFhRRYExdakUGRdGoqrK4EiizhCAzcIWRgUW5PAebIuCQUWIVcIBRZZAeyf5oWWrBnhzXuF0JJFAIorskpo1RoN3YXngsLqiqHAuj4opshGoeAaBN2F54AC68phh3td8HqTDcP+bBoosiaCNyQhG8WJKb2kwCJXAPu106G78AR4A5IgdBuuD4omQnqhCzEK3YVTQ4FFyMpxFikKLEKyYL83HFqyBsKbjGRDi9bloYAiZBZo1WrB0YWnQnFFBkORdTkorgg5CxRbACiyToMCi5wMBdc0cBZ1QhYHhRZF1igorsjkXIPY0gLICaLc86Z4ImS1XLHYYuD7UCiwyCxsXUT45xf7rQPOGYBOyCZgv9mFIisAbxQyK1sVFH3n5M83lbMPIWRVsP9sQ3ehgjcHWRSXci3quCdXhtT/hBAS4Irch4zJSkFxRVaH/9w6sZMbA0XRRAg5E1cgtqIN5/6cpSCETIQfXO6vzxFLdNsRQsisXK0li9YrQggh5Lxs1KrF0YUaCixCCCHk/Fxb/3tVlqxru7iEkI0j9j3ZlJctByEj2JBVi4HvFFiEkE0hAUeEKcPr/W2pdLH0+ripbXp/QnrYiNC6bpFFgUUI2RR9AmkOcoRZbj6EKDYgtK5TZL12/9nNDz7YhJC1cAkRtXTYhi8X34I58v59++GbJirQ2bk+kdUSWCH4wBJClgbFVT8xVyWZD99FPCMrFVqcJ6tDLJ6BEELOAQXVOHS9peLPyDD67kfer6PYlCWr13p1Kk7J8wEmhKRgh7R+LtXex46p16fSbIgVWbVoyZoEbTLtu9kJIdsiNKpuY50aUbhrO+QaT9EXhMRdqCy891bBZixZs1uxUlBkkWvlXLExU7v3Yy9KZ4w9IYT0sxJr1nYD3y8qrnKJNeAUZ2TtDLXkDrUC54qd2FxNhJBNsHCxRXfhRaGJl2yVoQJLL/X/p7rh+VwRQhbIai1Zq7BgTQGtYGTJpOKT6HYjhEzMQi1a2/pA9NUILKDqpGIBmHobIUMJWZRi91vu9lj+hBAyAWvr/1dnyVpbBZ+FlBWBFq/tkBrdpq93ajshhGyAhVm0thH4ToE1ExRipzGVmJ3o0xSEEHINLEhorV9kUWCdGQqvMH3CZ6rRcoQQQnpZiNBat8iiwFoYvpDYuluSwogQQgYhRVd3mHIeKbEAobVekUWBtVFOiRnL3S+UjpNNkhUghcCUJthRnQvdIeryzNVRku2Qe99OeS9dWGitU2RRYBEGcpMt4johX7gsQVwNhaKL+Ay5fzcitNYnsiiwCCFrImTlWbuAGkuqHsh1cGVCa10iiwKLEDIlsQZfW5JC1qVrEESXJlb3ZN0MfXZWLrTWI7IosAghKWIBtRRJ28K/jhRe62LsszfVdT6z0FqHyKLAImR9DHULUPiQuaEgWwanPOtTXMMzCq3li6zX3b6xfcDAyBZCyHhSzxGFD9kqoXuefcp5OLVd8eMZx+T/trtvPKkMmUQLsj/H0cfgV55vNqYZmZB83PNEMUWujdg9zz5k+ehrt9a26+KWLN+CNQcc6UK2DOOPCJkXirE4sYEjS2Nmi9b6LFlTknqT0fBBImsgdD8vsVEjZ2LrX1xYACFLcGhqjq32HTnTk+j1pOGilqxzWLFOZesPD1kHbLiuBP1FAn8JpCfhDaXry8sdU/8mJ7H2sJaY8WELzGjNWl7g+xoE1hSs5cEil2crDRmAsBBYWyfOLwqkWdv1JASzCa1liaxrEVh9pEZOrOnNh3QJzfQdSrNJUuIkt2P2rS9Dj0Mui75mdGeShTGD0FqOyKLAOp2YIKMwuxybEEy+aAm5mvz0fd+QLATgPUli9N1jhMzExELrugPft0Yq8Dk1G/ZaYwSWxmZG8/XF+bh1/jKUJoSrHwotEiN0j/n3I4UXWTFntWTRirUsQkKBwqvL6sWUT44FagpcvfGeIlNC0UUmYkJr1uXdhRRY18eUgm0qy9tmrFAplhKrRJFVoe+1nLrItfz15UsLIgUZ6WUioUV3ITk/Y77CnvrIry+0xgqmzQqspYgrx1I7+BzxkUrjtvlLP01s39wyDmHM8VJ1sBWBRrcjuTBnsWS97vaN8TlbCCH56OdmaaJqKmKiJSRsctOS4bh6TdX12mGfdB2k5owz5RTWrMu5C1/34E3jdgxNnhfaRsgWOVfcVArt7osJF4oZAoSv/9YEmYP9zvmJTag7IW978g2n7L5Cd2HOiKbY7MYaPhBkqfTNyD20IRka+5PKR++vRwn6xwmtIx1ErBs876W2tY9myP5nJXT9Q+u2ILxSM+6T0/DduyHL/cos+LNaskZbsc5NarJDPjhkSuZqIMZ2aBRHs+OLpVyhFBJZQ/NYDVsQXyHYf7RZuEA6wZp1fnfhagTWnKS+H0a2ySXcfG5ggG8t2ZgrT0RgjOksU9t0mr58pkg7OYV3H5Xj2o7cc1g81xQbds6+ImeAwBlcdktgpNA6r8iiwBqINodSkK2PSzc6MZG1IBYpYJaOE1ju3GvxPE+bsMT7Jputii/gPH1CztceNiqqQowQWiuMybom9M2bM8v2UC79lrRWhtT9ORqgHmvU1J3k1MImlJ9bF1teLb7Acv8b02ybWGxNVecXEWtDYwPXJMbm6BOWeMyNQpF1DcQemKEPzpZF2RIbkQu4+K5e3CwdJ7QWTOweWpSl7FpGQpKLM7m78HUP3rRdsy0ZzhBhlvsx4iEfLl6L2fsMgirkliMBiqKyEsWWQ9P6cVVAfB+f0DUKtdm55V0Yi44RYx921Qx0GZ4nJut1T31OboHa5ATo8oYna+ZMVikdlzW7iMoRD+Q0ckXWlMSu4QUF22JF2JYnal0TsdHVJ05M/LZHX5+bdOExWTn+9L7K6atQPgDk3JzJOjVk/SiGiCaXhgLrdGLXcG6XYewaxkT1GQjF6l3cCpaaP87Bfmca+kRUap8Lz+W3DJE1BX0VekoF80EhOVxQVE1KqDOlaFoWS4nNSt0XMwuw1GCJxVi8GPuVz9jvdS58eprJ3IWjXYVkOKd8143kc+Zv4M0+z9KULjzXyftLPw1ZLqnrdwkRd2Z35CJdkI61tOOpvmcjZLoMoyfM19M1Ukj7z62LLfvS+HnHtm2ZnHqcARGp/yajKJo/f/0YRJo/9zu29NOS5ZK6fjnbp8a/b/VyBkuqtoDpZ3AR04jktkdzt9M5x12IW26pTGLJohVrw6TeTtZiPUs99Gd86wrFk5zM1J0PxVEYf6TeqaMIQ1bFmCVHz4u1pBGDzgLmr5uLUN2doS4WbfFy9LVjY9txCiYAWdasaEWdLLJe//Tn5uzfzmzkzbr4G51cPZO//U4pok4tW2wKA7INLiHezjFScmp6hC77qPUwpL1+60tfl8wqtuEiLaQ2zw750/u6/0NL/1g5aRdlKiaLJXSfnOzy00JFu0yGChj/+Ke470Ium1PL18MUz3HftlPzOmWG+qH7nr0t8q9t6DrH7oux+Pfo1O7I0Hnkln23a5ZF0Sz1vi6NItZv6SUZzpC+e4immJvVji7cqtDy3Uixb7uF9slJ69KTNqn6CqXVyyz6pkHoa/Bzj9WXbkin2NMZxe6tU54bCq1h53EK+lr1tgmnCK0xE6T2nW/foIsYfWX3RVNR9Lel+pyOx8455l7rVbglJ2BIXxXaVy/71i+Bk9yFY1yFZLmcMhv4khuGMQ/zyUxh5Tm1nBNZmmavr75YklD6WLzI0uIBN8Lkz/cc8VU6RmxMea3AOrkNPB6r5QSuyjULryUKnlNJuAyjJ7taSxaZnlPeBpb+QM3+pjNG0Phv4GMtY6nA66xiTCg+QwMfTgmeje07dP05SE2tcqkyxMozUIzm3BuDhEDIkpQ7ECCGLqNf3ljZXDprsar+F0i8z+xi67K2Cu52VdmdVcyJrhHkWL/OQUzoLb3dXwKjLVmvf/bzhh8t5/tasY6ib5ROzvfDljIqhyyXnO/JDUW7NYaWZfCh2o3hqEZwihFGU38/8pzfnBxbvj766kDnkVN/fh5zfMBdi7EzWA4Hi4axn/6xz7fc7Nv3+Nh7zK97Wy/1+ZRlI7bYF81P32jg3BG/Ad76wteEVi/EkpUKpj1131Seel1fhfYNsR6yb05eMYbeCLnlJun609tPsU75wipH7AzsLGJu0Gxh5TqWkKDKmfMm1CHpjr9vGSuPK1PfseZizLFyxVCqDvx1Q+sxJd5i16qPPoGdMx3AAPx41P7yFenfkedIRIAb2/3FBFbu9/D8fU1Zb5NdUaXZ7WCKoiu02C4PR7uX+wwrof1O0SEj2L67sK+iY+lD24fum5NXjKE3winljjHWknOJD8v6xxta7rGCyhESVKcEoQe2jQnwrvLK7Cj9dGPFRg5DO2/HFHP4jN23Lw4sVwwNPfYQoZIj3vrWa4Zax3JFesY5TeIii4mukLDyXzqagjT/7+zLU+j8/PugjqE/2mykciVqoeWWYwYHbJEh7fCZxdJYRrkLR7kKCVkivsXJPQ9jYw0GmJ1Hu/aGiIxT3XTueDnERnuFrHtjY9GWSm7HH6uLnJFyffV7Ck4gzDG5sL7/xpRrYFlSfZrsnADaeW5CJbr6ypuTRpfBldmUQGmq8jkRdcp1ywm/mYPUsRYmcOYi4DKc2F2oO6SYeyS1fezw262QWyeppSZUf0Ovw9YY0qGnAmZTxKxniYZmsGsvJqj6rE45LqoxsVep+0pv988rtT3WaQ2JS+rbP1Q3Y5c6n1C9tI6VUUdDLKCp+h3TljhyP5PSNzt4SATF7r8h7kq97BFcvdPWuHIU0r0PU+VtHSSQxr8f/PotC6Aoq6B6EZjj8bT2N8eSM+ZLAzlhMgu1Gp1MTr85kNPchWMb1lg6xxCBEdovlVdsn9xj9YmdnH1DdXFKA9x34XPyHHM+sXRLuE6h8zyFXPdu6/ADY6ZyBRUQb9z97akOUm8PuUf6rmPnPvLKNMSKloozGrP0958iz1jeKbEQqqO+OsmtO/+4Y9qSnOdOk/MtvdzpNXKtqv55po6VEl67XbNP6J71j9H3AhITlK68rl5b5XXPZQnBri205njhnSN0ZQ309b8524a8/CS4bExWbueYe5IpATLVsXLETt++S2TM+cTSXfo6xczZEzYSua49F/hauyn6CMwgXeO/vdcN/84VqidvaS9D7g7d6YSyy/7GmZtzaG8F5n7wyDFfpEbTPXwA8/hJb7qTjhEQzM3IscxrWziR7eqifa2bPN16f+lj16fmEutDvKXGlF1xmMMukJkxwK4YJiJCQeWdY6EtRne7+EuHs1ztJPwSopPbdMa07+3Wc7/z7xl1veqy21WuTO5UygLYAWJjtIyL0+p7hv36m0OYTUVO2c7ZJ+a8MM5x2KExWa9/7vNnK8wm0CMfQtvGTkEROkbsOFMcK2Zyjh0jVJ6xZbigrz80ummQSy/H7ZbzBp+dl2eFGOJKix0ztD546PjMzdl113deofNUy/d+6mvwgd/9S8CTu2S6Tt6OmKt9gGWh1Yb2TXmA/jrp+4JDLI5vkMD0y3nKx95j91WuKBsSvuDHjqXK5FsB9/v271r0JuaJ01bdrFPx7gW/nK5OjGnWW4FlTphP6+whIX3HWrJBYSLe+v6/p39GT3i4Jescnd3Q6Q2WyNgRiEN93jkj1E491hBTc99ozgnOW38KZLbZ3N2EhLniqbBvx0WPWyxHNDlCFoE6XcDiVIitB7s8uuOFLSLB2JUi0un4aW1DKnZb9NMh+qVj7OjWxL6v/OFfw/2HfSBufv7dpx1jbFnKsjpX24lF68/h14vLu5UkXJf10uXjdWad46U6WP/esr/r8nvLKrtYR72LPB/qfksJIt+aJggLNGcNc/iXTZfBWRelaPZRz4J/b9d1WuftvejltCvqmhtTWa9blrDSqOe2bK6Bm+7B7jdqNvw+S/8couecx8oh1dZcUFMscwqHLfqIR5JyY1z7bLtZ3wTre+MdM71BTmxRTt71PpliKkTRPXcpKrdMtotSlSFdl6rjSbk1Qx27flZz71v9gd7YoZ5/ETfPv2jdQIF0fW/xufvEyuJ1gk5wRQXq0JiYIVYIvyMe0en55W/dV6l8PNHWabP6noWseC3PIutbJIPB+7vWb4nVSega5XbKgTxC94Gp9ZXpWpatAJTjsZpPyz9+TrD6VqZ9OLVfX1jc2SB34etf/ubJC7A2UjEcOXEdRNHnlsh1X4S2D52LKORW67NIpVxDfW6qWBkGHKPjltNWktAQ8dSbul7fdx56ve60hrreQi7BELHrP9S9NXY6gj7Xi78+Va+p+spx76TKMsZNlAoB6Ou0h5Sls2tz3yZj4vqua8gVrvMpdo3Asi8fnfssR+Tm1mWo3rz9mlg+03YfAjB39+lj9HHp+bZSISAb5K2/8Rb3b7RxH2TJGiIScl04oe19bqC+mIVU3lOQetuf9Hj+5HipUS5jhuOfm1DwbMhSFNpeWxLQjqtw/2uP2Jh5oeoh1oJWQHFIPDlC7jzfjeeXIeOa+m+/rbfvlItId+Ta4hKLkdjt4h2zb4nJur8kvL/GdTAp11ZIwME21OJtd79dFfvb/Q5yJ/EyxMrSt/3ojyYrIucAr968OjamWVdi3DNfBs4tdW3dPaOFgG+1y7meIRLxluK+7wfVnvplBarrCSTclu5ZDT93siva10XT96LgXzvfcubj6ivktirLlou9cSPayUmLKmBfbm8a8aU+zdM7LYVjaLjJHCyhDIohOiRV12NDU2ZzF+bOA5Qzi7Wfpi9mIZX3xThF9OQKkb50Y0iJlDGTW+7UyJux+baETCCdHnKfK1SB7vqUuy3TjReNS7HxUx2h5jXkknIZxPZz/7tGWO8X6jDr85RmqeurnsBR2kuN7sRSz11sW0eMeG/ArgPr6/hCFqO+MsUsjjGh5W/X99UxUg+xey503YB0XYe2u0vYOY5Elmjfj0OuWcr6qOvLF2u7XfsapawcpRKuUG154EW6yss7L/GugX//OHyRmdtnxMS33hZ6dsuyFlp1+XdaKJbVtbxv6jjWj2ULr5Ux1CiTky5Xh/TVda6uae2b6y78tFd8QXamm+aSlqKcmYb9dJci8HbthtoH0+XmOYSQJcPf3vodqddsQRU6RGbnlWtedy6PVF4+qaHzIul6SllEco6dgSl6yuCXZ8p0AxBnWdgVkGNZL2FM85cqQ5/AHFQY77xyjp2TPsbJs76PuFZB92HELYiuwKhd6M6yqI/n47frqTY+d+SlK48uV8CN3yq3G2FoTPX/qfXeKs7peeVOeRJKv1W+431fDWCq0YU5c+Tkuq5S8Syx+Iuhx5iLU+N9hkzSGMpLL6Pl8ywoc9dTqFHXb9vGwNzdNcf339pP6byT9eovPTr3XtFZHxvN2Cmj7khjsQnubT7kqvD/XDncdQu5K2KWGmcxCHS0LZFQGphCqqVaH+ygfYtYU0HKfaef0zLduem8d0XldtPWQ/fbX8ae/9z9NUdPVIe2u3qrDX5SBycbu02OZVWPKn2w3lWbVte7/0zELEc9rta6DPqappaJsnavPeLX3C9LyioUO4a/XV+z0P1VpzvW96X4Ama/b8qRavdqa5fK278PdNtUuDy99Z37MrC/dsHadkHUhKSmtvTbnXPb64zYxKmEzlBrziIFVo6AztEcGSI4X2TFXFJ+YXNdV6mb55LusaHkznKdO4t0H2POdUg9+Y1Y35u4LxT0MpRv6rhStN8HcmfGBsJLJzZL0/1WWevQSkR5jatR5xdtXLSA0ufjBJVeuvQuDdB2y5VeA6/FaquuVF37YgaVtcUmaARUaWD2VRnq5U6saEDVyO+VSyd2zdz2vvupiLhoYnnud+37yZaxtXTpyrJZX6j9Q/n5+WhC6/q2+2LNmFps1fWslgDq+9J4t3EtaFOu1pDADTx39TUtAJSNWHZlMG6ahl1PWZVrUYvvYEyTLpO7L7Uo9M9Hl9l/EdDbXb3oZ6I03f133j3jj2JNxVSGfsfOyy+7n96YeL9UeukcscEQLj50v+++BMQ8GKk+cczcZ+cmN6xjiBFnjGFhqK7JyHuZUzgshaECKEdA9d3kcwjGMW8SOcKptX0G16Rfj5266ekgHbH4p87huucQDCzvJuoO7Y/VX+j6+lY/P4/Yb4vZ75JpzD7QwYfySsagZYjkAZgT3m7NzQ5y703ceM635aKA+K6S/a7qEJ3Fy5ha1HQsMCErpBYXjtpy4t/nXgcawh/44OOEmV9Wj0aY9T3fzf6iRUVyl0jb4gs2P+7M4YsXfc7+i0qsPL7w8SmKxtIINJZIX0D6+bTElBIQoTKoUYC1G64Q4BAQRWMme12SUaKPqQwsCzrXfJHlj+Ya8qHWIcfoi4lJfRn9EqRMi6l9UqRM90M4xS0ZY8h5omk0cpaRDPpdC36nlTOnjO70QsdMHdt/kwq51/T+rhH1LQ/WZeO7kQCvUwt1BIUSnr4bzW13bkpjYIoCUpbV0naqclTLfdEtq/6th4NL4y5rVZs3Z5DfccvRP8dAGXSanVS/1dIUArMrqv42sL1vf3hl6Ajv1HZ7DeRo6k9g1GUXqYRfAcihbCyMXv3X+XiWnPq6uOug6rEue6CMUngCqUQjAO2yPrZ+dvz5nfz791jG3SPunou0U85S1sK3xoWmFvCfbVffvnuzrkMteAqvzEWrflvPXMd9ipaQ6tSJPg19fYwJiy7fmueWOjbSnb7bzwnE47GxrO/3VmCpF4pWu5MxiKguuL7evjn1zNM8bImMPnW4JWvMx1SHxCH1FfoSYirV2Ljtfvpcck3TURGSqGe9Xb/lpUyp7qHNOIfQZJepDyLXSzdZpl42GYTFSuccrHuuLO3nMrzGRjfGvpXJ7+S6J9F1U+ilbjC1ew3oiBu3PmhtUu4544aa6zLokU++e670yiACFI2FSIypLD5OvJSVBQgAjFRlNDe7ZmnTu3qWQ9nsr5atMrrO3ooXc+s1J94Q9jp/AEYAOZYob3f1sqq/ttXA3BStZb3/A1vGwLLFHvZcmnw0OiYqtV0HvRvnBTXqvFwZ7o91fWKn6u3o3ddWqNbCqlQWx86Ep2gEclkC+6aMxtVzfV1MnU/109R5O0FbFb45hpQlzG5X5V2oe8oeu3rWAvf9sWyWoefER7c7vhtWiwUtjoDG/bnf1XFngLQFjvGeTUBZoRq3Z/2s+Ut9bP9lLfRc2jo3hbVu9Vr6PUHquxJrq6V6mQHQie0KWXL8gP1gOVJW6lAb2BOiEZtTcIgRJqUJco0Dqe1jP8A+sSHpPO7CU+OQpmSMn9btF1pqQlaVIWjX2xjLk45D6mStypMaHt7jUus7Lzf/jZ4HJyhk/EYnln/MnO+CWmPD5V1jE2rsU+fguymMQR3gqj0whUrnj/jzBZPreHzcPRQ4Z2chAdB2rbnT8Y/p3n7rt+CqsS4fKnGi8jHiYrNuKhGyb4uY0hNMZicda5tf5i67zv9G91f22DG3by0M9gI5eFYwAHCi8ak95GDqZQeXLpBP+MCZ1gEA2NljOpF4u4ccy0qAGYPjzQ7F40MjVrxrVtaxcVX9lzc7G//U7jSNvfDlbo/i4Dfuu0aQPPDEhz6tGzQddi0UDUonBt3N1Xlh2XWtmUAl9ko0z+Cttbq5Zaeu7FKXw5ZBuwKNSNsqql56jBZf2pJVFJBDMzKv9dLirKoBy1bnvnXbo+5Rl76pq1Y5vHStZybU5tVWP5vvsazjzKQoqo9H7/fNqMMQuqhjAuVD9BlBYn35UCNMbj4xhorHMcc+MZ56GzFZsaC+MUH3QwmJgkH7RwRV7ALmnmNdnIAJO544LXjqYxbtRtcXUX7Ad8zdqd8Y+9KF3tT0uYfeFEMxEzGXnu9O0IHiAfdO14VYtDqDemSUFuWhiR+18HEWCCeudEwPtOgoKmHk8rOCqilzla60HX55u2sFMJc6n2NZH7d8sKvEihUj5U1hRRCC21s4g5XtnAo/YNfloSwqdRmMqdeXe0FxaJbHBwWKY7MsdxJc9qVz+NuDZXTnEBFk2s1ZHEqYvdT7FXdHlA9vUDw+NkLZWfpurLXNc4+WD5vrBIPO9apcZVoANxaduhPXmsS6VKtz8MVKk6aqkLI+JwDQ7kld1toiay1etfv5RrkijanuRb0MCQu33bYdRu/vPKpOwMGKMP3So2OxnOtSnxOkHRIAlR7ovvR04p4iAsyVXQSA62/89Wi9aLs4yM4Lit9+6Jdz184cDlUbeneftpz4ljadT4yUsSBUN0sNmJ+KPkvgCeTPk/WBf3zyg2eRctWNsUiNIdb5D84nYaJMMfBcO99N88VA6LMPYwVVTLx0C5WXLnSuKfHXitGIvK1627OnLvAFlf+7HjkmjcvES6fdd3VRPPEUbJxDabSQKgRiKtFSHMpaEBlrCTHuzd8JqZuiJZzcby2oCivAikPZEj31pdEdvedWcxaq4l6n6V7jwh3rvmzyMKZeX1vV1FKnH7usjx9Z3zpHZz3ah5/P4l4JT10/91V9uu1yXwI7we6lQyVWVKyWFsTGWS6du/Jomuvj9lNxZXJoYpRqV5e6v1qxb0WgszWorGUAxKj7sTRt96gta+vlQMe2FbCWLLSF85ilemkQY5rjunO34qp+jsqy206UprEW3uy7YsEnFktbes92SiwCzQjAkGXeuzYAupZGvd0091/dTh+PMKFRhjkuq1jZpyJnlJ9jRiFzMv59MIDveO9XAZhinqxLRumfMHwyixzxlHQxuUbI80+n0vpp+s7F2x76FFFrvS9EdJzHkM8c6H1jH9ONmt79c3Lpeo4f6tz0cRSdgFWg60oLlM1YF6Cp3+yLxvoVcvn6+/suvH3RTtdyzUkdKN0y7Wt3htNZrYB3eNtVTJS75URQPtw15XHp6mPbxIXAPKzWlc/uURxM/Rti67GO01FTAWgyHpPjU8gbOfhUk7kYg4Pdp7zZWQtS7Re1+VbLw9Pt9bFlebOzgs8/CT99G+lrX0277Pocjk81xyjuDPBUdS6lFVRO/DiR54phCkG5k+oaSXNP12UxRfv3jamP0ymzfhcJiCu33ujz9+KCxNikxqs/g8rVh2rZenG4aS/LB040hrebnRWy7vnbFSqOXR23FoCogvydVUt2LS+zlGV3ZO4u3GbUFGi3E25/K3xbeaml6PUh12PopVDXVW3d8srj8nJVrnzrsivqNr4apZjhsjJluGzR9CPET+4ov751l6TTb02b/TbchUM41RqVzHvg1Rl4s/VO6lY2EyZWfW1iioPYW03fMUJCKmVly7FwAXVgeB1XkRJUoTxDDZu2OIXeqLSFyh/9Z61QbiReO/jULSsx5wKY63J6QiiUf5WnG03myoxaUNUuPy/gG2KF1a218twW1grUduGVN+315QOpOn7XZ1kr0/GBoLi3brZ7YwVXc7jDQ8HuicHxQbXU9Xu8BXZ31U/9v8OtC20DgOMNsLtzlq3meh5vnOACinvgeCvY3ZlqqSxlhweC/ROD40213i2bOmjy292bOp9OOXT+dwZaqOh93P/S2u7OsakfMbBxRVW93rx4xHFfYPekrPN2VrTjg0pcHB9W1+f4sKivkzHVOezuSpQ31Xaza65zLUpQiRVn+XJuTuM+jwS0rIPFfVMOMdU9Uu3jlL4rpH02lTVN32O1IcYbsVdvt89ZZ7tIWxlaq1zp0knRxLU5XNhZWT2TrY3uGTyUbQsU2oK/Izq9dMHtukm3VjtTOhFs6vJX25VVJ2Z91xZ1XR8Jyzt2O5jDoVVfTZkSwmjodA8ht+M1EHLbTki+u/CDvnDyg5+VqcVVxOUnzz4N88JL/fsPEFh94spoE7PNW7QFqu/t5RzCKvMYrZFevssw5BrU23xXni6jv90f2u0LKucK0bPRq+NrYWXcW7M+D0Fn1FG5bxRw1ZlIW1jdVKOhjFi32m3lMqqEFarfVnhV8UqNQKo7YieY7gyOViAdHlYi4/CgwO7O/rbLvdt+VwkqI1X5dncG908L9o+VqHkozWu8QWtbK41dvdeizP57fNBsrwRWI3CKexuLdTQt4VMJrUrM/OqnPsErv+cBfu13GrzqhwXlDbB70hV09fKB3f6gWe8ICr8HTVl3TyLbFa36eVDV1/GBoDhU1/jmxerFZ/+oxPFBYZdWuD4UFE8MDs9U1/PwoKjFoKuL3Z1BuQd27veTpm7KfZXO3YdOiGk3bXHXFmRAYz2CsXFnZZPGuT+rkZOVoHQWuCpI3oq0ohEwvjszuZTIDPT6eSlN/YwA9vi1Bcu6Ft2xTbOtNVDEF0c+RXd963n2Lc6+y89N26GO38KPY9KuQeeG1fn5IyXd/6X9fTxWsWsxclxcp07VMJfLceV8x6/+HQBTuAvXRtRy4rn2svMLubC8da9+ZSV4XnpU/W4NS/ZciRGRkvUJAj9N4KbP/l5UrhCKTVMR2ifk2uu46yL7ue1uEsSnql5RHt/DPLyp3Ab3hyZ9S4hpV5t9Ww680TYNaPvcDWy5XQeiDYFFk65J3y5r6zjS/K9H1NX5qzoo3bnuq+kBjAjKo0F521gSqnT2uAXqv91dZRG6f6ZyGVZu0GZ7bZyo9ylaadz2/UsGd8850Yda7DXWQ3eC6lwgyqWl1hr3WxqXiFHb7O/9Iyu6IHV97e7s+d80B909Bo4Pq5/P/thDPH4lYJ67x0sffAs52nweokW7zu05PgIOT0fS6HN0P8t2utpK7NVB6MDFAZCj4Hhj7LkUkBKQsmjXoUC5rFV+xtadXe6tuMIz1fHvn3ZlUeYsALsbJbIMgDpd+5y0axAAxBmQjQGMTWdU/mWTjx74UN4IisLedzc2P9er7G0ehek8T5XL0t7froLrOK/2C1EtqFxcknMDujxKdw3UCDz9vPsG/doC5a3u3ADxttHsbJ2EQgt03tqy1dpehMWZy0tbyZzI9cIgWgYSvz0Nia5T58dKuRwpvKJsV2SFGDt1ROZ+8t732X+kHahoAv9L29KS/X0n/2EPvd24fI/H8BxSGfnW9Iih3vUpt18sjRZbNzsUj+6rOrXTCog/nDnWGEoTeN4SSRErW1CIRcrm3H2uftoxUeha4VQMVb1dGstBNSJMYOwTWd4UcBYsU1THK+22442gOFaiq9w3aY4PKiuQs1RVeSqBZTub+2cr4XL/MoP9i1UHefeyZiTU/ctK7F+UavlSUVu34N9Grn9x/3sdeS0S1LbDMwY37xf7v9Qd+v1zVVnuXl65CDW1MJJq/90jwf7dNzg8YyAHweFpX/gA989W+R2eMdi/VB3v8QdV+wLA4dnqHKtronZ0l9GPM3Z1kHhMD88Y7F+ojnn76y6GrLLEFYdKGN0+DxyeAW5eUHm7a2Oq/I8PgN0jwOyr6wRTrXP1WN5UFq3yBiju2vWrrVjlvqpLvd2JKC2AxRnkrYV2d28q9+OxEjFSVuUSY1CKvbf2Ur8cHB9WAexiKr1TXwtrkXLipWVNE1H3j7Ug7URN5yBVU2lQW9U6E43Wz7ZBPZVETPS0BEkgZjPURsas/DE3YIzkVAlWLIZch36ZU9M4DDlmfeyRBodgXqqeKLhabE9kJV1jwyYRG0tnJEiqLHOMkKyDEb1YoDmZQWABgDw51qN35PGhHTeREmetN1lPWHnLTqB2wGrV5KvTKcGkylCn8dwOzj1YW5xsoL0RtNKVe8GT5wrcvlC2hBLQWL2c26u0nXN5U3Vah6fs8mkX24Wmkz4qofWyqh7vXm5QHIHjg0qwuIDk++fs8hVHyFFgHhyBe+9l46asRNTtEebxDnJXWNFUlcHcGBRPBOWNAYzUguD+2aqjL28Ndk8q8VgcgLtXlNi/WKC8NShvgXJvUBykWd4YFPeCn/zC/wG/7Wv/BO6fNSjuqnyKO6mXjrtXliieCO4+oERxJzg+XblUndXm/lkDc1sJAbkTlA9KFHfOzGiv2a2+3wC5F5hCWRNuDOSugLmp9j3Yc7t/WVnV597g5vkCZl/Vxd3LreCyFqnSXh+zB+RQXZ/iUAnL3V21LO5g45ZgXclVDJiYyjUpzr12tGUsm/KXt6jPpxJVlWgSg3oiTbNDy2pW7pwgstexqNyKThg5QeBiyeRY7V+6/KzQqkSTu6lN82ICm78V7rXbzQqycGyUNNEZbmJXoG2Vchok1halOv/Q9gyPQXT/jkhSFq1IPn5MW1+e9Sd4lgZdii22GZM1tagYawFrWbACb1jqm3xjLVkwphuT5fIsimYkYV/+C7VkRfP1LUWRtCdbsiL5T2XJMlLl5SxZ9ZyEKnYpZclynaUpqo53kCXrucrK03IfImLJCmHaVqymQ7WnbdC2ZBlryXpB2nkoS5YTGzEO1kJ1/4zBzYuBe8uuOjxjLVk2fVUgtT1iyWrcmaF72ahRoF56tC1ZN89XlVkcqvMpDpXQvXkBODxVuS7rMluhbApYwdu1ZDlRU4lTa8naVwKsbalq8k1bstC2ZBnUFhU3eKCOxVKuQyOC/aO2RaW8aVzQleByCs8JIKlFVmvuMV2m0G/rjnvxQx/i2V98VOVbtgce1C43N71DypIVsxK57fqni2fdF7VodMiT+84nonotY5EJYnX8VWdqh3oKBzsthWfJ6u27c6cimNLocIXi6jpjsnKCuQflN2BOq1h59P9+2fSowL5yhwSCf2NrS5Zfhlj+sTrTDaWfl59f6BiBdZ2gdj9N6o3Q/Xbtgp4HyBdQpomyMHXchukO3XYBtToQvtR1YTrxG84aAjfPjwua1TM9u7xg6rdssVYrOVoXAapTkaOBKSuxhQK4ecHYeCUr6I4GxaESZ8XBWCtEE8xtROrRgbsnxgorg0b0uT+pOm4B9o8aUadjsm5erDrN/Yu77jQO7vxrMeL9tuuCgguohJFpjnl8WBVgb914xT1qC5wTNPtHprLQAZXrT4Ddo8ZKsnsENX1CtZ9zEbp8q/8rgVOfo6/l9Tnp83HptEUzUC8376/O7eb9tmxldR2dC+/mBVPFV7lxMdLOx1l0do+s8FHusv1jG5Nl/4epQo18gVKnQztwvxOTBVTPh1tfx2SVtfANxWQVdcyW2xHYPTni8PSuyVMJkOrzSVasHZqytqxVdRlME5CuyvbsL3gxrloo6bbAbwd8dEB5QgyY3a6eksUNVjFQZb7dVzGvh7L7ofIQMbGTElj+eegJW3M4RyB8nc/1CashbFNkAd0Lf6ro8m/IvvlJknn5PVChNnVv2I7w8tJIUcCooTIiUp1v6LM2sUDL2Da3PTYLcGyfmGiz68R3qYoA8L4j5wsov4yugXXna1RD5I0uFC3A3AzWWlQJWqMLxR9dWEpdRhe7Uje6B6lHF4rOr2znX82BVQkfPbqwkGbaBhd0LkcDKb3RhTfu0yvt0YXFgwKmqDrz/aOyGYUWHF1YRkcX7h8hb3Sh69Afda+9E4ZAZPThAyWO7gyOj9FMEeFGzB0MDg/a0yrsHqM1DcPhodSCospH8OQDBI9ebfDwvYIHv2awewI8/xrBy37e1KMNd4/d/dYdPShlI2QAOzrSO4XDQ2UVK+KjC3UdmMLUowL3j00zCvTWXZfmXPU5+aMLd3fe6EI7mlDKqpzOwrW7a84hNGv9qaML9y/dt6dTQdWRyH0zgWxqlGE9mWpodGEpvaMLq+c3sK1+8Qq0SaEpWPy4u4f7enJYN2dZI7ybT/oY+wkl2ZWQ2z2KF59U7cEx8iFnu3/Qkua3Z/p/7aEoivY0Dp3zm0FU9Vn/SBbbdBf2MZWVK/SBTABBN2GqHMptmDxcTrn70vRNRDrUrTjRPFk565PzZIX29QWXLrPeHopdk4HzZKn9g/NkeW+oqXmyACuoBsyTdfdc1fAXB2NdSRKdJ6ueOgInzJOF7hxS+8j23vmp7DxWvmWpPT9Wd76t421l+XLzab3wYdUIvmf+VVkLUze4IKcc5b4SfNqiFZonq8m7mSvseItO/Yjb31oWd/dd0ebmqHJ1UO5UGQ3658lyFtSimmm+ngcL/nES82R5bjzfeuXPkyX3x2aOLkHzgXE/pshVY2ieLOW2dPNkNYJLWt9/bDJqhFTLNWiX9SegQu5/j9a8du7zRs/e2vP1R+F556SbkkPlypPHh2rUsx9wnnIjxixY7n/3IlqWtauw01/3iauhccgUVIPpcxdep8jymTqGy9FXt1pkTT0xqc4/xdSiq5X3QGHVtz0nviuWNide7JQZ3/102rWUkbYWGDrGS9Dd7gm4KohecHxQVG5I21nWE1KqfN3/zhLW5Ik6Rqx1jid4yAfjRJ83oSiQ7iwjMyjU1qLswxug5cZyt0Ckb6qnXXCW2URHXOdj83cC2O3fOqbW08qK0vpder9duf1jaiGh3fWdjjpynMA5ybHsWKHNTXfS45jgcttqi7bnbuscWwm+4JxUIdEFJNsYXefu00nV9ClF51NKLcGl8tfWPSfuikf3gPpIdffA7TJ26sc/P7fdWrDqkeRzCCsymulisi71vSFtKYktTyVHDKUsJ315BoLS22nKKman76PQ6lxT4rgWYDlWNPetr9i59X270B0nVEdHr9Gr0yVGPBbS/rCy3tcvOwAx6lz9MvixYXb7FN8uFGchAyoXR+DbhfW32OB1eIHvEooro/099NuFuyfNNwv7vl24f/Ew7tuFekb5ffj6tb5vqD9zYq9n3zcFb5z17VhNxrp73Eyy6qw+qW8Udr55aD8GrcunPx7d3EdV/bcm99ylz9EU0vl2oSur+4ZhXTZ1zP2L7e/2FYeymuTTWXPUdam/66fccK1vFxbStvg4IeT66X1zD+q5zoLfLgTsNwTT3y6UsmxcgG4/71uEReJbhcdnbmv3+M2Ld61zdHllfbuwPrjbroSZPxmxSqevqrm9bZ7xQgBVH3V5jJt4VO/Y1Ec9393t3sasmarcddq2daj324VA+3uIAMzBi/3qDUlR+59imUr1tZfSBJcgZ0qMCMv/QHQOvvDS6+cmJBZyTLqh0XE5wfUjzrGOEzI9E5QmxEy7DD3Cy/fj57j3+tK584ydeyh2wLdExYSY+t2JFekTXi4P3zW5K7piUXdWulxAe9Z4FwDsp7HCp/5EijuHfXUdzM2u6ZCPxgol+1uJyVLno6wSIeElqrPwP9wcwljhJAcTvJa1eFATWur1TuwdrUAyIpX42zUiMPWx4dpNameOb30nEOgIM00VE1SVBUD1gefIOUpp6rI0oqxoynpfhs/xvvsB7Fo47T0B5QSzjdGr0QLr4LnV3NLda74YcWVxsYqh2Mg6dqqs67XeTwTF3aG593aB5yi1tGkPz95i9+I96o9B10LJnSPirkF/MI6bFd39+SOqQ22m7S/KZx7Ya2MtWjs7IGInbSujo2zKVgtVdw0OZWXZevFxUNy0LFe+sPKtV2qdORzsOebG+mKYsAoZL0ibRL30fSB6GyJrCHMKr9CNnbpho9MQFBg9onHE+WVPH1El7jm+J8DculA+obeskJsvdI+OcVPm7KsFU0ywAY3LQwuv2DH8mDC/nIFz1p1b0N0YiPNqpXXCa6diulouTLvUHb4qo++i8zvbbBJV7UREdx9nGShQ3HXFii5/5xwOYVFRpevm0zpsq4OLl1tP0xGLgaqP6QmmZoOXUHVyrXTefWN2laWxKafrnF0+CF6nmMvVxRy29gndy/qZLlFZs4DaWgugG+QeY1+0LD2+kKm++xiIa03FL7ky6jKHlm4fEZiHN103oUgtsluf1wphmpGUxf2xuvfsecnjalRF0CXon4P+ffTO23oagkHvQ61T3pQaFFTT0Ceyho0uHNvxax/xKVMhTMGYmyrXVTlEYLnt/lQGzp0G9L+5AOgE3JcD67nwJrQLTdmg0W9mgYbbHA4QkXoJoDtTccpylbJo6Qb9eIxb0tyIQ52ffw7+rNH+du1u0Etr2RL36Y+iqOIwPLdeK//DsTtvFtBMqujSaTfoXfNGXk8tUUgz4WJ5aNyTel+7lLtD5bZ5fGgC95XVQVwZH99XyyeHVodZlI1FTFuLNLUITGx365zlrMXx0O6YPYqXTNUhv1TWorF1rd3vR/f1shO0jKqjM/sC8qjslMG5x6oylu37Rp2DPldX30YAcceOWXbrkXbhc9XWolY9+gamsoQpChRPDu0y2o6+jg2z6cQXIq5MgZF1xV3ZHrzhPxf1+mP7xaO0rjdnsb037aWPe45cgLi/TZ1r6yXCiUlpLHhVuoA1yJ2rSPuZ887NPP3AjiKsPtBtDJr7VwtEv+kyaI86dCJ6b9sF07QL4sdnxQRWzH3ogt3dPFmx+hzap5Uq/1M/tbMEYn1drB/09YifLra9b78I+SJLd/y5Bejb35Haz23380mlmRp/cj29DAmuIbFi/lQIfSZ3n5x6Tt0YsQc0tl7HhfnuMLfembmte7KzLMvgUmUQrovjsV0H7vdRT0KFrpXHF3l6e+izRCJVxxFb6rixw9EKsrIRWP7SNbieEKrcPcegRce4POvyHAE9d1dRBMpgIEVTB3J/tMe2nej90XbAtmO/PzbLvU2vh9q7/fVSH99VlwsC9ucMKoqWhUTu1Nu4SLtMbqmuhymKSiy65V4tnSjR57ATFHdNGdyoMeyk+nKAW2p27QBtsaPFtNutdQ5HAxdTKECwDPXSr/fQuUpPPWv3nasbZy3R4kOVuf62JwA3mKMWX8fAfS9STT+g79vDsd22+ff9UaWLPScaJ37qYO7A/FK+sDPN6Dv3rFTnogRVq32wZStLYKdiFP3pYNzSxV65FwEroOriRJpAPzjfuW1Roj2VQ0hUaaucP8VDa7u1FBp7v8VcejlB8KdsXwOpc4ht8w0ZqX40Z5lg3DxZpx6474T60o9NM5Qc4Ra72euHJ1Eu/bFoR6gxiL0B9WGOzXKohSyG3zbWn7II15Vx3yXLXCbRdVMfINSAdc9PfHEWIscVCjRWALfsc2fqDkeJ0tZHef3dY59Dqq1gWoxIc86dj+Xajsn+rpcH27G4b227Tu8o7d11Pi1yvqN2TNZpp0xenlKX3YqKg7e03zj0lwAql5NzSbn13jcRY+vaI9XCItxZjtwIxlhZ6rLrpTGdc+3eVPBeFAJlL01/21CisYjGLCKpfd3SmPit4BEbVdgRXjllcMKxDKQJxTg5Eekse/a+Fm9+kNKKV3OjLF7GdISVbxntuLdNk06MfdEo0bjq7z03nxZV7pz8dsuJN3eseiqHtoXrZHIsYLmWsqEWtXPESy+E7U5GOgU5AslPP8QtGkqXG8M1lFyB5rvZei1pbvuJIjenPjvxDT31rAJjTYGohbEeGOBb1GCtcaEg4tZxIm7LKgPvbTtQj7psuowxV6e3Tg762M3bdb1dz0HkOg19P7TcT6qsoWPrMvSJUlWe/mSmGSwQGTYPIO2u1+5WZ93LKWtKePgvQNbak4xjrC0q4futlc6dcyo//95qicBw+VtWnNh9BLQHeYSejbIJxO+tQ+1K1OuG7BMoY9Qt5t+PLoRAW7U8928Vi3XbnJ8J3KEGgMDGakm9rrKeuoEd7bagNb3EvgAe3EJeeKlJUHpldf/HBFZ9bmV+2zrHCMDc/ce6LKcg5ZruewZjdTahCBw2hUNObFJO7JJL54idZCxN7jHmJMcMOVZoxci9MVP12VdPoca77205R5D5+/rby2O3wc2xvqXcojo+zcUf+HXoxaTVb49+rBzQcXvWZay/MaYsN1p4ucb/eAR2u058SC28fJcL0BZgdZ2oxtnvRFzVOeHl3Ek2T+2eDE5lYbd3yuhf09AQ+r5r3omda5ex5Xpy+/tlcDFwelnn4wkBJ7SUOzX4bOhy1tNyqE6zrnNp0qTKGLqWfsyWdstpN11MmLipVlx9qzquBZAnOuprHrrWRXu7XvrnUltv+q65Fho6jcavJ/9zMX5b5ddjTFhpsWIMsN9bS5xpiS1ze1O5ZO1oXLOTjspyM9vXMVrukK3pNKCmdgi0k/eH5h718dOHBJYbTRi1Mgb6QLf+2ugLecmJjY6t0yNYY/qnh2GWrFRsUk66vnxT6049Rh9DzaI5oi57MrjMdH2xa45Q8PvUD2Gnvkx76VtnOqibuPT27eSJdH1rt6hewjVSdmm8Bk+k7frQ9arWG/3Zo5B7MyDIXMNZN8DOxH9oxyU1eUSCp30vUmgko5+Xw3WOu/Y+Lranjgmy1iv3Ni5+Z6r/90eO9pUhVJ7QOr9zj+RrjKmtM1KaZv9ATJscy6rDMkYJ4XQRg+cV6kRVWVOTbrZmOo+VtVBxRrn4giNWbldG+7v+KLQ/+i9lUQJGthsmfk1jxzame8/HztN7tvTLkrSsWrYM+x3k7r6ZlsWKGXFCzLZXZldUMXDGdJsku170vg7t0twVMGbXdRf6wkqvg3rJcwLreBwkIlp1IN6gpgvSVxa9fdBo97nRdT9Sc2RP4QAAn/aKLxhQOtJirT7ovm80DhlwMId1r4/Ym0gOqYe9r3wpEzYSDYlvGfCFWLQ8qbIm8swoa8xaFk3fV55cQmXNcT/lph9bBrc+lC6Wvi/fUzpDZ/mIEbvfx75wpcqaOn+9X2qevZhbPnDcYP9VGsjezkLvzZllnnvaWrB2lcXqtj1bvdkX1VxvPe+9ciyb2EBXhDpW8Ag5HiHvf6lb7oC4qs/Bj+M9li0BGQplIBXBMI8ziLXveN9X10WIpRkckxUbMdaXri+/IceIjlo7U8WOYkyDlisQ5nSfTj06Zax1LzTE1v/f4erAuQwdqTeQUL2lLEZ98Xq6wQwIsnpkZuCYrXtdx/8MmZ1fl9cPotfn5btiWq4c07hs/JFQMcuZdnOG8M3wen3oHvbdSDH3XJ1Pj5tTn39KDMQsa54VK5gupwy6HKHjDH3mQ3m4387CklOWVH3EtsXWhVyk2grqLxP1rtv2mLBqFSEwYtk89aASYHYUJ4BuoLtBr8CqEqIzR5oYU09IivtD262sz7M+UChf5SL0LHSTiqrYtcndHkt3IUJ1M7S+Uta0U7TGIJHlX+y+i597c6QqKJbH0LJMhV/RoW2pfQeRMlXG0q7JN9/n227FZhy9dHo6gbKdzxFNx5QdNOrFcQERa5USCbGYQX8+Ieem9N2T/mgtezzjHaP1G6jz9e+3jmsy1OnW5Y24/DoB1qqj9gl1inUHq44ZFBGB8nXcRL4QQ/d8Yh2zCOD7ekL32ZB7sJXOpfWXCJ9vrMypbc6l508UWtd14Nr4eepr4i+B7j0QEkMxUmK0FnjK7QUAx2P7nvXbU7+Mqs5bLvrUlEHGANjBFGimpjDWzVwYoLAftN633ejlTWXdirXhLUFWSDVnmV8POwEOyoWsRwLadPXUNcY0L4MO5S5M9mWxZ7fvmoXy8Je522PpckiJb53mjLRic0dokhiD3IUA8PqXvzkrYzINKeG2WKvdkuiLtRtj8Rsw0Wtefol0Y9ySwUP0pMu5l3Ib8aGN/ZSNachdOPTFK3ef3HSnvPgNKf+QOo8x0Uvq4JfdvpdC/eITenHyj1dUbkDRI1WfelhZs4D6o9b1x62NsW7Eoj3ZaITaXei/aB/Kyl34+A549LhTtuDEz/o7g1Z0mZgF7NrJeTk8E2/9jbfon9O5C8l5STVW1+iTz3Ebtxgz6AJIu2pTE+Jq92DMMuITszb1HStVfu94oXslNAFs1Azu3I5+eUPot3cdzxUTcqkYnL44pz6XxphnJOiOChwjN10OQzrSnHrx02VnHX+m0pYVkxb8Qy3rUSGVuAd1GQ4HG3yuyrTfVW68/a4WVfBFT2EgyW9BmcaV7gss92kdl26/B+7v1a4BgdU6N3vepa37qebDCjHk3hlyr81NnzVtCGPcoyMYLrKoqtOk5lSKmUj7lrFjxI5zyrH6/O6xY4TKc0oZIqbjXLfxWJr4p8Zt4EZkJQVZaZTbMvGMHNFjnTJenrFYGM8lGGsI9UhTzw0aHCHpu2U8F2V1yAGjlnS6Hld6qIOvXaOICEDlhqrKX4Y7iPbB0mWKuSFd3n3H8Ms0kL6YU79e9H6BzMIH8V3d2i1nRwZnTRTccnUBndjJHPquV6vM3vPh46ZqMSVgpP3pGH0PHw2MinkXYwDnPTQmbpew20Xpvc652KW5u2vK3Ypp9GLQ9IjC0uSJqylesPvy8F23uelz8J/BHCE3NWPdowMZ7C4EgNc/9/mTHJwslKHunqURE2sTmpSzXbWDRzQmypiKA2sXLt1wxWJZctygQ+eKGTinTL3PAM416CUodnrKFCtnzLI4mZttSH3HCH36bHReA61quWXvcyPudtVxpYA890w14WpRADf7yj0ozTxX5maH0v++ZoTicTM1Qx0LeSirKRsKgTz/YvNRZ30uWmA56kE69sUuNBDgWugT2wsJkXnr+/+ev2pl7sKYqXKMybJv36HHSm2/lCKfmthonxhLu059gZxj8ARaTmdYmfzL1jJ9DGkC/EP4c3pFRzCaZhkcoefPJ+bw6idUt36wft9Ld266VjG8kaM9U3q0LC45Au0E8ZH7GajQ9s6+3nl11vedf0r05NZ3X3tW/5+RYWpEZGuwwAD6hF19T5eR62onADZH4O4eeOpBdV7GVN/SdDFZIpC7AwqjXIkRtMCqymitinf3VT5P7iuB5VuuvPQdcQV0rZ9DBFbKEzElOV6PsfRZzeYQnGP6uAGM/HZhZofXly6lWGOmyjEmy9wLl3us1PbYtj5RMXQZOmbseGOOkXOuqfO+9HWKnc9Y8Tsi2DLkykxaXUIjQ2OdVs7UGqEZ7mNTJZQGwWkXUvfxWEHc26CV3f9DVpWQ8MgRUH0fT+9bhsoZE4GpvPw8QnWQ873RU55vt79epkjdP267XqYYahnrO2ZqP6ASVvcHmJt9V0g5oWQ/lh6zPsv9MVhPzfc0D1XclxNQrfvFK689/6jASp5TJKziHB6InNGGsbbxwlM9BBnbx2Uy3pKV2+H1pZtDmS6RMULylDobex22cj367keH3+m77TlirC/+LFisruBqinKC8AoGCkfEQaxDzAnUD9VnqAMfep/1iQEJnGvM5Zka3h/KY8wyll9qWyyPnJdOYJhQGvp86/nL+tzBQwVVzlcqYgwQcUFX7NFAdqiOv29GEgJWUClLlqtbOZT201LNsd2nhfxvFqqDN3/unFtLVW5dD6GpHuptPXGqawrp0GWdMch8iYyKyQKA1z/7edOXhiyLsTf/uczWwDwm61xyrDopBpZxdMzRnNNU5Bx3plisVTEkLs/f71KcEoMFJNyGI4sz5AXQawNkV0Ceegq4vanyeuimclB2Bpe8bP43u10111YCeVQFuMv9ATgcYR4/bovZUEyWjr+acxThNXEmofbWF74mtDr6cC8zJmtO+jrglLoeOvttTl655Thl5t2cEYupcozhXG9ZY4b0TiXMcqw6KWGUmnwvUJ6+QOqoCAtNwOlvz3FHZsRHJY87Jn3qSwYx11Us7yFfS8iZkmCogJhTMKXit/wvJExFTh1knvMkn41JjHQ2ZQl5CjYu62HtNtQfeTalvreqhfjBbV71+QLMPHqk0kaszX6A+4RM+YWVvq+5LI7cOfxi+85kVRsvsqaYnXXICLCppiMIlT22PZQuxx+dWp5SjqHHzt12bfTFEpzywIXcM7kBkwOvTSrmy5F0Qfrr+kRGTnxQTizTEPrcRENcV7FtQ9ynOfmNJcftmRMjlnKPDiU2mnVI7JVflETnnV+uyItj4hkyd3eQp5+uhNbtTSW09rt6dKATXO3jxIsg7uPNLtbr8ZP28wM0Sz8w381Ef4IVK/XZub79YuvGTImzis/aOVJt7JgX9gxGuwsB4PVPf+5JBydkM4x5yTglKH/CmIbJGsap3X1D43hOnWZgrrzOQY7YmWJahwA5U1oMZqqXQd/6+9yz1fKph9Vm6z6s47WA/ufwUAkjubejDJ9U7kLzkrVi+fdN6QmW4zF6ftcwwXTWi+BCeetLXxfbRHch6efUN5KlNhBn+Yj4GMtuTlB+3/Fib/ADRFfqs02D6it3vqzsgg3saKcURaeOfBvD3PFXI/PPtWYMK8uI0IX4wdMDBnTS979QCa37A2BHGQJofrv8gPDzZ4WVHJQFar+Def6FcEB7XUTrHgy4CZfabs5FbEJpnzWJrxQnWbIAWrPORUokjPWpX9vDnWJ2c/epFqcpyjVBrMH89TQg5mpMTFZf2tSxVkZufM5s7cApsaAxckdi9lB8wCuqf6ywqgPgbyJ2h/tD6xmsRdb9oYrD0hOPFmp+PC2qrGuQ7e7pXCJWLGHFAq7BkhVrMHo/T9GTLpTXpUj5yYf61Pmgt8mNQRhtFUvFjOgOKCaEcq7XFDFfA6agyGXIC8DgmKsxMVl9aU+IN4qR9akcj9zZ5fuOe8r2oIV0iGjKjQUdIp78uMfcMljk9hbm8WPI7W2dh9zd2yxNcNRwLaqMqd2F5tGjKq3a1008jLJspgQ+eJOXnpu+KTn0ulOtsWdibL+W289Pycki69QCjmmsh+SVKz5yxAmFCem7T0Y9oP4cMpoh1qfY/TmkTKFO8UQrXOoZHCMg/A8W53zAeArBcmqeY1+EFmVpmmIATUxIDRVPMXrKZu7ugLs7oDRV0PuD23qbE1PmZt/EXGm0wAJaAqteRtyCs5Iamau3xwa9+P/3vWAsXITFyO3nJ20nTnUXAsDrnvqcSQpDFkxsuHbO8PdLulr6GoMzvrVN+oY05XDj3NGPYxg6DcrWOXVal0sSsja5e2bOF9Ap3Y0+NzeQ/R7ihNbNDXB/XwkvOwqxFlKPnwD7fTUPli4bUJ+/OR7Hd9B9FiefofPQXZqVCjMAeNujr+9LEj25SUQWQKG1OvpMx2M+7LtmhsyRNAGzxDatRXhpQlO4jJmjLkQs31SnnTPH3JD56KaOSToXE8U+DeZMdaatFcXTT8McDpCnnwLu7iEPbqvfNzdN+sMROBwqK1iTSTMNQ1HADPoszgRt69pEVoyFuygzBBZwDTFZV4c/f03KqqTp++bX0LmCtkLovGcUXiFz9cmCK+bqG0PMdTO1+PI70dy53lIdceoTJLlzzeXse2raS+CuZ8giNZf7LkbfXFcT15nvBnLL8qWX2scrBMZar3B33/zWdeWE1eFg47ASAivWJk/Rtm6hffbdlX6/tVDxlctklizgSqxZ12TdIcM4g/Vrcs7lstvIcOyzseb4zwsIyovHy7I/WA4TtsOZVizgKixZQ4dzn3osQnxCVsMJmWUCv9RoxykFWKoTPIcAS8UT5aYN7ZNKuwVC7tEFWOUuIqp0n8I+YLnErs2ZQ0Ic04usudwsKfeYn+baXF1kmcxo8p51xuSUW24Oy9e5Okzf/TVkmoDUPkPyXSJ9MWRndnn6sYoXF1Tut16S9dGnD2YSYZO6CwHgdQ/eNL40ZNv0zaLtf5cttW1o2qUzsQg762zJ5xBipM2YoP4FEJqb8CKsVSzltoGn0Ne2Xglve/INQ5LPP7pQQ6F1Bazpu245LKnxmOoNKvBlgNi2yRn60fZrZqGC6BQuHiPlszZRtbb2dUnt5wQMFFjAVcRkkX5ib0F+mmskdN6XajhmcDP2fS9sVpdjaBliqQJsgyJoShYnqIBli6o5rVCXItZ+rs2jMAMUWUthqBl47MPp9vGXpMtQ9+YcjcgJMy8P6fxyPkfh0s1mBaOYOQtDPu2zKNYysjunDb+G9je3r9m4CJvFXQjQZdgiN4aItFljvVy6wTjj7PWjvuFIZqfv+62rYG0TIm/JKnVpLhwLNsJVCJw7JsuxaqHFh4bMwSXE1xkn84uJrbPGg62E1QieqVm6YPJhH7AeTmxfRwos4FIiC1iA0Nqi/5tsg0tavWYUXqcMv49ZxlIurljaqbhaMTSUtVmfUrCf2BYZbewJAgu4aOD7HJ1IKn4pJKSuwf9N1oe+L88tuGLzzk3IGHES+txQbFtOWjITfZ+KWYPA4ov39ZCKB0ttn4DZLVkA8LrbN56yOyHXy1KCQVf+/TCSYAuWpxAUUSSTt91946lZXM5d6KDQIuRELi24QpYLiq/ls1URpaGgIiOZQGABnCeLLAkzorGXQkbt5+exai49+3LIFUSxtQxCz0bIjbcF9HOghRUFFlkgZ7NkAbRmrQFTmlrQaFFyqsBZG6sWZJe2eAHpD+le6EOto4h9izW33NdgRZqaKxZLS2h3/bbPL9OWmMiKBSzBXeig0LocvoByD/AUVqJrIiRCF8ml5+waS5+YSa1LiaHYOv+4ZD7o1ou2w0vGL+9q2sAEEwosICGyVtb6khjuQTWlaf2fWufvS/II1eEiCc24vIbOLTRSTf+l1qX2j63zj0umITbj9xruwQlItburaUMsfeVey3lcgrNbsgBas0LQurQtVvmmtzaLFzk/nHcQAEVFDqE+bCnt4cRWLGBJ7kLHtQotPpzr5VThu5QGJgsKLnKlAioG2+7puUSbOIPAApYosoB1Ca2QH3pNPnWyXGj1IpOg4+9ocerFb8vJMpmyXZxJYAFLFVnAsoRWqKPjw0fOyaqEVgiKr/mgiDoZtufrJ2bk6Gs7ZxRYAOfJaugbGcGHkFyS1Y/ciXX6FF/9+HWX+kSY/p8CKwg9DdskNWhrqPA6Bxe3ZDleu//s+v/U8NZYpaXS8yEjW2IJDcekpCwyuVaa1L5T5BsSP7QkXQTdnrN9J328/fBN5zjMuixZqeGtKRGVk56QteN3MqsnZZHJFS+pfafIN2VBosCanVg7zvadLJ3FWLKAtjWLEDKMTQgucvVQOJGpOJMVC1hy4LsPhRYh07Hq+C6yGTiSj5ybMwosYE0iC6DQImQOKLTIXHAgEVkSZxZYwNpEFkChRcjcUHSRHGiFImviAgILWOO3Cy9UUYQQcpX0fWuPAossnSXqhsVashy0aJGa0GSMse0kiz4XD61d24PfRyVb5MICa33uQg2F1kpY41B2f/4jMgiKsMvCCTcJubjAAtYusgAKrcWwRiF1Clp40VLWixZdHNU4DRRQhMRZgMAC1hiT5bOQitw+sUkX3d+14Z/7NdbBAFzsTuhzF1um7/xCsU1brxMygtRnlE7JZ6OsQResxpLloEVrBLEPyy4NXcY1Q2sXGUDqU2FkoaylTfXZUNu0MIG1fkuWY2EVu1x8y9OarDFrbwjWVNeEnMjVuIRDVqa1PefaMn/K56YuzJp0wOosWQ5atBQreTCy8GOgtoL+mDAhHnoeKgayLwDfUnVtLHhA0EIF1voD32NsUmz55ue1mKLJaSysMSNks7CNHc+FXhgXKq4c23EXbpKQ2TkUgE62Da8zIaeTGrzjrwvtQ9KE6nFtbtMzsnpLlmPxFi2+MZFToJWLELada+UE9+PCLViO7boLNYsRWmsabULWBwUX2TJsN6+LRHu2EoEFJETW/pyluAoosMjcLDQglZDRsL28XkLxcRtiU5YszWt3n5VOMPZCsjEgS2ZjDRRZOWwvyUDefvzmSxdhDLRkdeh7+P3hu7ROkTWg71EKLnJu2EYS0mKzliwgw5pFyDVAsUWmguEQZEZWasUCriXwPQbFFiEBKL6ID8UTuQArFleO6xZZAIUWIb1QdG0bWqHIAtmAwAIosiootAjJgJ8A2hYUVWShbERgARRZbSi2CBkBRdf5CIlciiWyETYkrhwUWT4UWoScwIbntZkV/8PD/v8UUmTjbFBgARRZaSi4CJmQKUWXLzp03n2CZcyx+vKjMCJkMBsVVhrOk0UIOROnztWVEi/+uthHfkNliIkhX7iN/Z8QQjxoyVLQokXIjOQILooWQjbDFViwHHQXDoFiixBCCBnHFYkrR1RkMWI1wBXeIIQQQsjJsP9sQ5EVgTcKIYQQkg/7zS50F2ZA9yEhhBAShuKK7sKT4A1ECCGEdGH/mIaWrIHQqkUIIeTaobhqwdGFU0OxRQgh5NqguApCd+HU8EYjhBByTbDfGw4tWRNByxYhhJCtQWGVBS1Zc8MbkRBCyJZgv3Y6FFkTwhuSEELIFmB/Ng10F84IXYiEEELWAoXVaDi68JJQbBFCCFkqFFcnw5isS8IbmBBCyBJh/zQvtGSdGVq1CCGEXBqKq0mhu3BpUGwRQgg5NxRXs0B34dLgjU4IIeScsN85P7RkLQRatgghhEwNhdVZoLtwLVBsEUIIORWKq7NCkbVGKLgIIYTkQmF1MSiy1gzFFiGEZGAMINL9P3eflUJxdXEY+L5m+ABNTN6LRX46Qsg0GNP8xdaH0vj/u99+Wn/p7x/bvmDYPywbWrJWCq1bEfz7Wb+hTtlwioTzc+tX/mZMyKSsSLRE0c/8BZ9viqpFQnfhFqHQQiNo9H2s3QVLhAKMbI2lPmvn4IzPMwXWYqHI2jJXKbZ8UbWGRj5kVaPliywZ/RITi3Vaw7N3DvxneMLnmuJq8VBkXQubFFxrslJNCYUXmZM1vqiskZD4Cq33oLBaFRRZ18hqBVfoTfmaOwGKLTKWmDXlWp+lpRB5pimsVgtF1rWyaqFFwlB0kVz4HC0b9SxTYK0aiixS8driDc2PJXXWIZegv36tnON8Qpa+JV1fMg1beB4IAODt5bdcughkOiiySJuW2EoRczVM2YFv3R0Yq6s5z1cHKzOwfpn0XZetPg9XDsXVJqHIImmyRdcp9AWAXkOncmkrHcXWebmGe5okoai6CiiySB4XFVvXQGg4/KXK4Y5P4TWe0FQc13Q/kygUV1cFRRYZxlnEFsDOaSks0a04ZHbtKb9TF5sy5NLCmKwCiqurhCKLnMbZRBdZDnMJL4oUsiEoqggossiUUHBdCXONWKTIIiuHwop4UGSR6aHYInShkWuC4opEoMgi83MVousS0zEQQi4CRRXJhCKLnJfNCq4hQdOEkNVBYUVGQJFFLs8mhRctW4SsFgoqMhEUWWQ5XIXYosgiZLFQXJGJocgiy2STggvgvF+ELAwKKzIjFFlkPWxWeBFCzgIFFTkzFFlknVBwEUJyoLAiF4Qii6wfCi5CiIbCiiwEiiyyXSi+CNk2FFNk4VBkkeuAgouQbUBhRVYERRa5Tii6CFkHFFVkxVBkEaKh+CLkMlBMkQ1CkUVICoouQuaBoopcARRZhJwCRRhZBW4S3DNOhksRRQhFFiGTQtFFLk7qY+UzCiyKKkI6UGQRMicUXeRshMQVRRUhl4Qii5AlQDFGlgLFEyGTQZFFyFKh8CJzQ0FFyKxQZBGyNii+yFAopgi5CBRZhGwdirLtQdFEyCqgyCLk2qEIWx4UUYRsAoosQkgeFGOnQ/FEyFVBkUUIIYQQMgNRkbU/NQNCCCGEENKluHQBCCGEEEK2CEUWIYQQQsgMUGQRQgghhMwARRYhhBBCyAxQZBFCCCGEzABFFiGEEELIDFBkEUIIIYTMAEUWIYQQQsgMUGQRQgghhMwARRYhhBBCyAxQZBFCCCGEzABFFiGEEELIDFBkEUIIIYTMAEUWIYQQQsgMUGQRQgghhMwARRYhhBBCyAxQZBFCCCGEzABFFiGEEELIDOwvXQByHbxKPtjc4a76IQJxG7r/AKL/D233fvvpvE3Bdck8ABPMp5sulL8JpontC0Ck2ieZxss/M20r/ZAyIXEeGfn0lrHv2EOOG91mBl4Hb1+1PpiNuIXpbvLyav1s5WcC6Zs8/du0tb7ebtR2L129r79f97jtvJsySyd/08mntU3Viz6tWP76GLH8W+cbykMA16r4aav/9bam/fmhH3vyVmPMp4GQmaDIImfhDnf4hOK1kEIAKeolCtciC6Qoqla/XgpErLG1kGqd3i7Nvs32op3W225EKvutt91fX//W28X+79LBigG7ze1Tp1Xbq2V1HvX/ge2t37DrisC21v7dZbVN2us625Heju56vT0n7+D6VL62sw3u729vlcE0aaDS+tvdpWulr/6vxYT7X5q0Iko8SJNe7G+9vah/G3vrmTrvAu11etlsr363/qzYCK2v/i+b7Wrdrv7dbN/Z34UYu11vc+tLFLBLu6/br1mW9T4A2vugOnaVxu1f1sfbuf3tee1cvvU20867/u3Kb9fZy70TYAdBAWAnggJif4vdLihQVOul+g8Adh/y068CITNCdyEhhBBCyAxQZBFCCCGEzABFFiGEEELIDFBkEUIIIYTMAEUWIYQQQsgMUGQRQgghhMwARRYhhBBCyAxQZBFCCCGEzABFFiGEEELIDFBkEUIIIYTMgBjT/e4WIVMjIj8O4PGly0EIIYqHxpjfeelCkO3CbxeSc/HYGPPxly4EIYQ4ROQHL10Gsm3oLiSEEEIImQGKLEIIIYSQGaDIIufiqy5dAEII8WC7RGaFge+EEEIIITNASxYhhBBCyAxQZBFCCCGEzABFFgkiIm8RkXfb+a3cun9dRL5PRP4fEfmHIvKywH4fIyI/ov6eF5EvUtv/tIj8CxH55yLy1+26N3n7lCLycXbbl4nIL4rIC/OfNSFkLYjIQxH5pyLyo7Y9+VK7/g32dykiH6/S34rI37Xt14+KyKd4275KRP6liPyUiHxm4ri/RUReEJE/P+f5kW3AmCwSREQ+CcALAL7WTdYnIv8MwJ83xrxDRN4M4DXGmP8ykccOwC8B+ARjzDtF5N8D8F8A+EPGmCci8mpjzLu9fX4XgP/NGPNR9vfvBfBOAD9tjHl2hlMlhKwQEREAzxhjXhCRGwDfA+DPAPgNACWAv4OqvfpBm/5PAvh4Y8wfE5FXA/g/APxuY0xpBdrOGPMlIlIAeKUx5j2R4/4Dm/8PGGP+5tznSdYNLVkkiDHmuwD8mrf6YwB8l/3/7QCib3uW3w/gZ4wx77S//wSAv2qMeWKP8e7APm8E8I2qHN9vjPnlgcUnhGwcU+Es3Df2zxhjftIY8y8Cu/wOAP+X3ffdAN4HwFm63gzgr9htZUJg/WEAPwvgn090GmTjUGSRIfw4gE+3/78BwIf3pP9sKMEE4KMBfKKI/ICIvENEfndgnz/i7UMIIUFEZCciPwLg3QDeboz5gUTyHwXw74vIXkReA+DfAvDhIvIKu/2/EZEfFpFvEZHfFDjWMwC+GMCXTnoSZNNQZJEhvBnAnxSRHwLwHIC7WEIRuUUlyL5Frd4D+AAAvxfAfwbgm63J3+3zCQBeMsb8OAghpAdjzNEY83EAPgzA7xGR1HcI3wLgXQB+EMDfAvC9AA6o2qUPA/BPjDH/JoDvAxByA34pgK9Q1jNCeuG3C0k2xpifAvA6ABCRjwbwhxLJ/wCAHzbG/Ipa9y4A32qqQMB/KiIlgFcB+FW73bd8EUJIL8aY94nIdwL4NFQW91CaA4A/636LyPcC+GkA7wXwEoBvs5u+BcAXBLL4BAD/oR2w8woApYg8Nsb87YlOg2wQWrJINjZYFDYw9EsAfGUieSu2yvK/Avh9No+PBnAL4D0qzzcA+KZJC00I2SQi8kHO1SciTwH4VAA/lUj/tHX5QUReC+BgjPkJ+9L3DwF8ik36+wH8hL+/MeYTjTEfaYz5SFSWsL9MgUX6oMgiQUTkG1GZzT9GRN4lIl8A4I0i8i9RNWT/CsDftWl/s4j8I7Xv0wBeC+BbvWzfAuCj7LQQ3wTg80wzvPWTALzLGPOzXjn+uoi8C8DTthx/cepzJYSskg8B8I9F5McA/DNUMVn/u4j8B7bN+LcBfLuIvNWmfzWAHxaRn0QVW/W5Kq8vBvAXbV6fC+A/BQAR+XQR+UtnOh+yQTiFAyGEEELIDNCSRQghhBAyAxRZhBBCCCEzQJFFCCGEEDIDFFmEEEIIITNAkUUIIYQQMgMUWYQQQgghM0CRRQghhBAyA/8/3X3MDOF+Z7sAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 612x388.8 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "data = hp.read_map('../data/models/healpix/dm-ne2001-30kpc.fits')\n",
    "hp.mollview(data)\n",
    "\n",
    "nside = hp.npix2nside(len(data))  # This sets map resolution  (NSIDE)\n",
    "npix  = hp.nside2npix(nside)      # Compute number of pixels in healpix map\n",
    "pid   = np.arange(npix)           # Give each pixel an index\n",
    "\n",
    "def dm_func_healpix(gl, gb):\n",
    "    pixloc = hp.ang2pix(nside, gl, gb, lonlat=True)\n",
    "    return data[pixloc]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "beautiful-distributor",
   "metadata": {},
   "source": [
    "## Compare timings\n",
    "\n",
    "Create a grid of galactic (lon,lat) points for testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "id": "binary-duncan",
   "metadata": {},
   "outputs": [],
   "source": [
    "gl = np.linspace(-180, 180, 360)\n",
    "gb = np.linspace(-90, 90, 180)\n",
    "\n",
    "gl_grid, gb_grid = np.meshgrid(gl, gb)\n",
    "gl_grid = gl_grid.ravel()\n",
    "gb_grid = gb_grid.ravel()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "radical-building",
   "metadata": {},
   "source": [
    "Now run timing test:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "id": "significant-translator",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.14 ms ± 16.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit dm_func_healpix(gl_grid, gb_grid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "id": "sunset-kingston",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "746 ms ± 1.85 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit dm_func_sql(gl_grid, gb_grid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "id": "phantom-nancy",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Speedup: 237.58x\n"
     ]
    }
   ],
   "source": [
    "print(\"Speedup: {:2.2f}x\".format(746 / 3.14,))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "compound-domestic",
   "metadata": {},
   "source": [
    "### Compare output\n",
    "\n",
    "These methods are pretty different, let's check they have the same output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "id": "unauthorized-garage",
   "metadata": {},
   "outputs": [],
   "source": [
    "v_hpx = dm_func_healpix(gl_grid, gb_grid)\n",
    "v_sql = dm_func_frbpoppy(gl_grid, gb_grid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 139,
   "id": "discrete-plymouth",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7eff1469c150>"
      ]
     },
     "execution_count": 139,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAscAAAD7CAYAAACG0TnRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABQWklEQVR4nO29a6w0yX3e9/y7Z85533f33V2RS66WuxRvoSLThk3RC0qWbEGBZIkiHNFKYIFGgBAIAyaAhFgIDIgCAZv5ZjmRAgQC5FAwQ8qQRMmxFfGDY10IKwwCSRTv98uSWlmrXe7yvrf3PWem+58PVdVd3VN9m+merp7z/ICDM1PTXV1dXf300/+uqhZVBSGEEEIIIQRI5i4AIYQQQgghsUBzTAghhBBCiIXmmBBCCCGEEAvNMSGEEEIIIRaaY0IIIYQQQiw0x4QQQgghhFgmM8ci8joR+ZyIPCwib51qO4QQQg6Hmk0IIQaZYp5jEUkBfB7A3wXwKIA/BfAPVfXTo2+MEELIQVCzCSGkZKrI8WsBPKyqX1LVSwDvAfCGibZFCCHkMKjZhBBiWU2U7wMA/sL7/iiA72la+EzO9RrumKgohBAyLU/jG19V1RfMXY4DGKTZAHWbELJcbuNZXOqFNP0+lTkObbDSf0NE3gLgLQBwDTfwPckPT1QUQgiZlj/I//Wfz12GA+nUbIC6TQg5Df4k/4PW36fqVvEogBd73x8E8Ji/gKq+Q1UfUtWH1jifqBiEEEJ60KnZAHWbEHI1mMoc/ymAV4rIy0TkDMAbAbx3om0RQgg5DGo2IYRYJulWoapbEflpAL8LIAXwTlX91BTbIoQQchjUbEIIKZmqzzFU9d8B+HdT5U8IIWQ8qNmEEGLgG/IIIYQQQgix0BwTQgghhBBimaxbxeKRxPsYngpPc/W/TF0iQgghQ5BA/IdaTQjpgOa4CSegklRNcNNyhBBC4oL6TAjZA5rjLiiuhBBCCCFXBvY5JoQQQgghxEJzTAghhBBCiIXmmBBCCCGEEAvNMSGEEEIIIRaaY0IIIYQQQiw0x4QQQgghhFhojgkhhBBCCLHQHBNCCCGEEGKhOSaEEEIIIcRCc0wIIYQQQoiF5pgQQgghhBALzTEhhBBCCCEWmmNCCCGEEEIsq7kLQEgvxN7HaT5vOQghhMSPtMT+NOc1hbRCc0yWRZvgddEmgl1Cui9+vhRhQsipcYgm74PT0UO2669bz4c6TUBzTObg2GJ66HYPLS/FlhASI3Np8SFMXeZQ/tTwKwfNMdmPJYrqXPDxHSHkWFCbx4eGeR72acsjHZd4zDENxPHo0+C6joPfZ4v0o2998RzoD9sguSqwrcdFqDsGu9EdxhhtfAx/g5jMMRmHsQSUQnx8KKbD4U3avHTVPdt0P9iGl0/9GNIoV4mpjfcoS3zmeKoBV4TEDiPLhphElBwGzXMJ2/XVxX8yfhWekp9AW4/PHF9FTqAhEbID++mRLoZoX0xtR5JxZk0gV4tTmCXDtf0Tb/enZY6HHqx9GyOjIWROljJ7hr8dV+YTF1QyIbG1ndjKQ5bNSH1lJ9v2IcsvkNMyx0MZqwtHW18jQuZgbBFlmyaEkHmhDh+Nq22OD4GNlMTM2FPgtN0MEkIIIScEzTE5OSQRAIDmOnNJFgZNLyFkYiSRRm3uo91umS6o/+QQaI7JUekrbCH6iJ2f/yHb2nf7hBByCoytn0PyHmPboTyo4aQvNMdkMsYW1ynF+tDtU3QJIbHhorRza2csNNWDqyPqOHHQHJNWKKr98IXVrzOKLSFkCvoaX2p4N66O6nVF/b66HGSOReQRAE8DyABsVfUhEXkegN8E8FIAjwD4SVX9xmHFJGNBoZyOUN22RSpIHFy1x6/U7WXQR6up59NCsxwvQ9v+0KcDY4zA+c9U9dWq+pD9/lYA71PVVwJ4n/1OJkYS6fVH4iDmYxJz2cakbR+vwPl0uG5LuJ5IP064bZ0ssR+jWMs1lCm01386IIkAHVlM0a3iDQB+0H5+N4A/BPCzrWsI+3MO5RROAGJgdGJ8jnV+FNvJjrK5KRmu2w30qfsl9/Gk9kaE/6bCWTYf11OnJZxPS+kHf6g5VgC/JyIK4H9X1XcAuE9VHwcAVX1cRF54aCGHGudTmMor9oYzOmNOIxaak7cuoBG/5XBusxzTedNkoK7c+TEuR9HtNpr6eBYFjKwNXlmWML1jWxln0PHYDPMcnEI/+EPN8fer6mNWSH9fRD7bd0UReQuAtwDANdzYuwBtlRxzI11C4+hNbAIaKs9Yr8eMRGyBeNryGIQGMXYZKLI3s+t253b2OOaHnA9Xoo3FptPHIBIdnzvgMSZX4lzBgeZYVR+z/58Ukd8G8FoAT4jI/Tb6cD+AJxvWfQeAdwDAXcnzjtZS9j2wfmM+6cZxFQV0CKH6mSnKHOvNX9OTm77nzUmfXxGwRN3uw5VoN/VuBO47dXsYM5vmvl2PxthOVz5L6eZwbPY+o0TkDhG56T4D+BEAnwTwXgBvsou9CcDvHFrIGIi5E/4gJGn/I8OJqE7HHuDTlV/bdmIfvHIVuWq6vXhCWhL6PgVX8ZyN6HrYd1DwPhodWoZUOSRyfB+A3xYRl8+vq+q/F5E/BfBbIvJmAP8RwD84vJikN31PaHcyHCPSmMiw7dRP1LHLGBKCKevBPyYR9WUeuPJO2U9SUPuePzMexwOhbsfEUAM2VEv3wT+vl3iOj1k/ET0pJDUmvnnZ2xyr6pcA/I1A+tcA/NCw3GR3R9kAm4k9whsS1ENE9hgC3bYNiu0yytiHsc6d2M/BBibX7eBGT6Tt7MPY7WQMLVqi4R1C0/6NpeP0Ksdhco1tPw/ifUPekIo51cY5ZePoKxSnLqR9GFoHQ0V4qYY5BpouVAs1ryfJqRvoGNoadbqbrjra1zxTv/sTw7nSk3jN8RAinparlbkbyhUTVBGB6pG6kQCHRSoi6YoxC6F973uuzH1Okf2IUcNjaUuR6XRdR20XndkYTdP9ej40ynxVDHMs58gEnIY57uLYwhtLg5lZVOcWzRBDynSw6I7VdzqS6YhGYawp9cjV4pTbwRF1eixNjknb+5RlsJZPMS5lCTp+iufZnv30r4Y57mKJDeLIxjcmMTwWbfu8l3EeW3CPMfl9aOqoq8LUfRdPgakHz54CR9Dqq6jPQ2iqn0E6PtVA7qukqYdw5HFL8ZjjqfoDLZ0JhXV2QU0iEYV8uJEcRWyB6czFmIK7dPEe+xyK7DF3VAypm1PR9Inbw9F0OhY9PpQBen6wjh975qNTInIdjcccd3GKortk43sqQgoM25cO4aXYTkjkYkoG0jVDzDGmLevL0rT6lPR5KF373sM8h44JNTyAO0dPUJuXY46HEFMUekmiekxBnTtqHaKP+DXV0R6m+UqK7QJFtH7sjjKocyFMNsjVtZM+7WWfc2GB/XwBTKvRMWryPnS1x5E0/OAuGcV2I9CT2KdabWDKIGAU5lgwQf/ONhZ0gR7l4E8hqKcipD5d+9TWDvcQ3NENc2W7Iw0mPLGowBjn0+zdkSLB1cJRB7mGOJWBx2PptIjRqmO20yQxWjeVeQ/lXddW93tXZLipDQ7U8Mm61l0x5tLTrq1GYY676Ft5S4zoRCGs7IIRpi6KbfU0RHAHGOa92/RYghuJcPsRSprT5dJ3ZoGjTbvYk8XpdNuyU+nxlDofyruW5o6RdpnoprqZSMNN1vG05bE5VT1ehDnuS4wmenGiOtY2exDTSVWft1NV2/f/EMEdILZddbQE0V3U9FGhC2vXxfaqIrL3I+r2bKXyvwn/ZumQ8yCabg/7lGMkjY5Jiw+lMMmuTfTV8Yk03C9TE/u035393HP9RTGmJ+nY/zjMcZvI1hnhwhRlo1iYqB69DseIYPpdBGrdDur70/nIbGzBHSi2ZfbD6iUkpFGeD/sy9g1dj4gVCTDiINcm/HZ71Da8z/E/YoBivIDMQnUhb3/CVNHAPjp+iIb7+bSw9zFLEkiMN+snoJFxmOMhDK30mBrOvg1m6IkzYDsHC+kxBfSQKcU0N2V1eSQD28UYgjvUNA+JXPdkViO8ZMH0620B0fqjEWpPg/tc9phdwPVpnZpjBimm1OkxdHlpUzi6+dg7Br/1fmlIV8Cib9sfYfaMSdYFlqfJR7x+xWOOh+z0oFGiR4pIRyiqvQX1ECE9poAeJPhp9Ws+sNwdZlo6x4C0dNOYQnC7thNiaULpOLbhP6VI+xSMreWuXc7dPufWaWCYBu6rzUuNGCPtHoicovsFSQ0Gem/DXM1kN23udr0vbuCn//2EiMccD2EKIz11Ax0kgCNGFPoK3T5CeqiIznkypQMjgV1mus08t0QrepnmQ2bR8DmWCC9dJJd6sZqDJrPgp7fdlE0VFNmXPQywpN6+rspLql5c2ix75NlHS/to9L6zzMx1ztb3SfNwWl+crre2lRYT7Z4w1mnT8DHGpgQzPuTtqUc6nrH015+AeMzxIRUzp/ge2JdMRIDUi2omsntHBgCbbUM+HdsfIqhd7LuvsT+eq+9WSIw7zXSL4DYZ531NMzAsWrEEIhNG0gNpGAxXP5Zj9T+eYl73PZ6miAiwXpnyiDHHeted0Otn0HVq9CEHktuXkOduA5cb6K1b7dvp0sg2jQ7VSxpIG6rDoZewTB5VTnumWZo01xWzyVingWus21Yoz6aIc0c3jVbjXOTRUMZjDkCechq+CQjpzuCuRx2Lx2GO5cDRnDGKb0uZin1dr4AkBRKBJIlZRwQ4W0NXKbBKgSwHVCHP3oLevm1OXP8kbRO8JiHrs59DhPQAwYxhMNhu2wqIsS+YfpH9Y9EYtUh383DJDeLd1k2j0zgPYcGiGEPbIYcdh05j3cW+/ZEHdtUQEaNzaQqsVua7CPSO68hu3gBSQX6WQtMEmgpEFfm1FeTu60ieu0TyrTX0mWfDxmuo8S1+20P7m/atiZDJDmdSzq089RzL/nZaymfaVouWj6Hjbd00evRvHlXLD6FhWryYmbqMcZhjSOOjjGKJPSviYPEdgUrZEzGP3tK0aojNgsjuvRu3X3QDmztSbM9NuuSKG1/Z4vyJZ5F8/WkThcgaHv/YfMLp+4npoLo/pMHOIRJ53r5/rv3URLhsV5747iW6exhnBCI6Pek6H05OFOvRr3o7P7GXnByXBt3uQ8cA1zaCg1+nnH3HmmJZr8uZlWz0OL/7BrIbZ8jPjClWvz6uKSRTrNIEmqZIkwT69DPt/TRDGt1Rx5112HvsycA6bBqkegwNaQuW2RulnXoJaHlQxwHTPtt0vK+G19etIW4jMbwlz+eEXwtd0r5vkZhjGFGoN6JDDswB4usITUa/V35uPyQxUYdVWnalkPICk73wHjz18juwuWFNsd2spoJnv32NW8+/G3c+dg1nj3wVevsCyLJyGz1FtbH8nX2XO4Rzn3qJwYh5NyaNghua/iwUqZLqcsHIRZPoDjHO7qchffFcETvfC3SCtOmIP4MJUNZp7F2BYmFoPfWZUaBrk1NMWtEYaKhpttVrXa+Q37yB7Poa2xspNBVAABWBpkCyVUCNWd6sEqRnCZCKOW2fea5z2zs63aSVTbrcV1v7jFmp31zuY5pCg7dCetumw468IzLt10koL/96Duw+dahFpHd0vEvDd7bbouO1RSr5TU3XudvSm+UqEIc5dm1qzAvS0Km6fGrGutFQDhnsZiPGslqZkzeR8iQWwfaFd+Gpl93A9rp3UvrnWgLkqeDinjXW184huQKbTXe/m/r3kJi2PrrrIZ5dHGqCxzDRvYxvj+04gXP1WM+3S3g9A72X6Fa2ETDeZD/cRd/XoJOOmoyAYFgd1eu3L/VBWodoexNNwQVfs1PvKd96hfzGGvl5inydAALkK0F2Zkyy5ECyUSRbNaerJpDra8jN60g2WxPYcOdxvStCl0YPvR717erQtFy9OKm33AhdKDQRSK47/82PoW4oqG4/VP4uvXcGP9ewlrfpeFJbFsnuE7m6FmveouN1OvpYNz1lp14No6O64jDHwLgHdl8RLsrSU3wHDHaT9dp2p7DrpGnlEVx2bYXtNSkNcUN13Hpeihs3zpFkGSQRyDarnbhemdoENVTf+5rkrn5Vc46aDkZ99zSRqtULRYsxbjXRrs9ssayLptXzbBDdehU1meg2+g4+jSG6X8e/INcNRv1zE6ELaz16H+O+x0bfOlIttc997039BnLMIEpL+Z0xXqVFVwq1A6b1bIV8nWB7PYWujMZlZ4Lcu6Jma8HqwpjkLE2M6dueI3nmNrCR8Gw09S54frrf7puMaUO6bzyDDNH3LkPcd3kvXVPzXZ1chnS2nk+fc9xbrmK8XZNKGjR8R3fbgyFtBto8gV51G+gQ/pOskMF29O4bnlRvNpuelNXT69Hs0JP+EyMScyxmYFqdfSt/6DRduxnst11HTQhkvQbObH81ETPYrrbsMw+eIT+rRouB2iNEMSLy3HfcgTs/fwlNEkiamShEU79KX0xbyjhYNPfpnnFMk9wiRNp0eNvaSshwNUUafPENiW3dRA810DtlDZjoNkIGex9Co9mPgQz83CefPumkhvQ3qvU6bbrB6yKk7X1uhirb7oqg2pmEzs6MoU+MsS3O0zSBrhJk5wnytUATY+bcGJEiGwW258BKzWfJc2TXV1itUiBNIJvt7o1Yg1HWkKa37UsgXVPsMU6l38mgjdFsAKFD25QOQEJaWT/GXcc+kO6+aR7Ip563Xb8IpCQN+l3R8vDNn/haXjfQoXwqJN4itfblbix6tP1q99DUW7/p2lG/SA6cRaQPoZlQQl14DqUxn/a2HYk5RsNJO8Jj431FePB2Gi4S3gC8QmB9IUkTY3Iz023C5OUV13VN9vY9Xwk0TU2DTxPIBaonrztpkprgdna56BlNDglvmziG6CO6+xrkNhPrs2dTCIq397mX+GKoiW6JuqmGT/OmEfxjdqfdiXZPwyDzT46HezImgSndBlGLDve5MA69GQrR0IdW1mtrgl0/42pXOLWRYE1NQKPsc1yWxTxossuoIl8nSC+3ZpYiVSgAyfJmnW6LIlcKa/V+Xw12xy1JwlFnMbNv7OTvDG49/3r3hnq0N1R+p1V5aWJ9A+1rritPpez+dkP72airtfxd+Wx5FKi2Q9WgVle0vL5cZbvpbn5NgZda+qH38o3rj/UGyqE3q6Fodz2tYzYSpzt+d9L69658mojDHAt2dyZE57QtfdhjzkRHj8jnzn6s18B6BbURiGIeY8DMgZgASBKcP5Xh1r1J+ThJUKiEqBFfF0VOL9TMp+nG44kA26w8oYq+a95nv/FURHdnB0I7VfnaJJJOYFTEpLVEOOoR8mp+3fXcdCEQ1eBxlNCh7WozdYEK3amHxNX/3/S7n0+9f12T+NbL3HUTMOQtgFOZz33yDYj0oIGwTRccMi4CM+OO+9rHFA25aAb0PjRIundxm9pQWtPIJDGa7bpT2L7GFV1LpDJgWhNrkr0rqmTmvM1TQFSQnwH5ZYL8bGX0IEnKp35Oq32jW1wLJFynaNBit06TEa1koJXrhgrMgENniBMp3apbVe3nyjXFapRI4OlnoDuH2G4erurNHUS1v3GROSrXwyILa9hlR/fUBI1CTw+9m67impDnph7rmu3rkK/F2a5ZVhEzg1QtvdLm/aBIIN/WtD6/ddE4D39gxpeh29m3XAcY8qYxYb2vFR2LxWGOgf2jhMXqPdbvOoBD7i4a79KT6uf1CrpeFVGHQsxSK0apSc/XNsLg+xnvxt4YZGD9jOLsWxtgZaYOgqqJIq+S4qStiLgTdaAqnC37URG3mqhVxNj/mGv52EpsFMU350V+DXkF8qxuuLpMqzkOrVP/zafp/OyIQhTiWoift52aAKsfSfBuJFrNtPdbMIJRL2OtvK19q/uK2RDR67Ns1w2oMw2OpvNsTONbdFkJRE+a0vz1rizSPFvCPvgRrIZIVucg6S76zO7gjLE1x+q02gUzgIqZylPADcJTL5IqmVlOV0CSmUf06UUGXSfG760S0/e4Yk6r2m22hYoO72ifM45+sjWpolaHfaNsgxmaCmS7+zIhXXkGPZB/kW8qlah7kS6l7klWGmBdVeu+srynl5X0bNfQhrZTMdCekfYNtG+8jTEXm7+pD6hWNV1LLa9st6blbjvFegDg5Vu57vhGvM0kh24mGzTP7VPwxsIn6XGD2vRb181trru/1W+OQ6Q9umZMFuRo15BIzLGEKynmyE+f7gar1Iis7TpRmR5mnRbmWEVw/o0Nrt1M8dx9Uh6zmhhBgfNv5pAsNyOk4RmgrUDdCRmIArg7/Yrx7fEYriKQTWlaFYDCHHvzfgYHGobyrC/j9rFmjvOVhNu2YneqJ2/doEEONbOOpueLqCiqZtgJc65WXFE11FruT1GevPq72rybzLJZJ7BPoci0IyScbWLYROC3RlGuLxvsM+2lhfymjdRVX8TSImxDdcPXnpAONQl4H2E/dQ4MalRw9Vn/PyY7N+u7EU1NbVDDabYLYkh1HuNko5DcnLOZ1SPJUZzbri+y5AC21qTmQH5m+3zmagZVW6NZ6DTssqlnUj3jVzfHRZcHqW4bde11+1cJaiSF8TP7h6K7iGQmIFNeh9w+GqNZ13TJ7Y2C/ey0UHIto+veekmmRVdCX0+TzMtHbTkS1EwvwhHpWrrfHSTZqh34510j1Bpwm7fTX8nMsS3TPbOclUaxch3IvUi2NdRQNTpcGOIGvXbaV+uLXei+O+b1dWtPUbSvDvdI25k5pG2dUB/yytt+W3Sis/tUYF1XjgmfEMZhju2JuZt+oPA2NarO8vTYbkBU6991lULP1oVhdZEH08XCRQmsWGSKm49eIF+d49YLrbn0zSOAs6eB1S2169t0e0JJmpQDrayIF/3gEt8o99xHW1/q9YOum2JfsHxD6tJzK4iVCIhn+CufvXxNWilu/gXH1YsT291y28U0lKY7aY6dbhctzcWPVhRiqv53e+HLy887y+Zay8cTX6AayfAi041RaXjLtZnmHcPs71ftt67v9UGOA9fffXyqpkBNItt7KiRvOU5vNx2+bjs9nKO+/X7DXVNatXUzcJ/XK+jKvLjDdXszmuMePxsdTTa5MchrgaTmHKgP9pXcmLJ0o1g/kxVBDcCeb2epia7WtNpFWZ1OOrOXrwXJpnre7Ri9BMjO/KcvKPJyTxvN9lEY7jyVwiTnK7dcmb8rh+kiYpZJtl76CsbU2nTkdr0VIFsgX9vv9oahvh9Oy/JVEkzfXhOsbpU6WOz/Ckg25feiXKnZbpG+NstpYn53+66pucmB2HSru5oIkq1ac+5dg8TWsZYG3dWtmdu61PJiHzLd1XJ4+lc7Z3xdVHjr1LouBgMmxQ67+uupyQEtLmYOqV8cQzqe5+Eumt4TkMY89ppAoRZlaxun0JTeYYHiMMdAr36mw5GGz0Oz6WHcnQH09yNJTPcHEcATO13ZPmwr2TGcdz62QbpZ4fKmYHNHGY2447Ec17+2BXJFfuYGhZTbktSOrHWm1EVtvUiAf9e+vSZINwGD4opuBSo7c9MXeULtCXbRZ8wJtDOvCcpuIpV1bZmkujzqebt8VSCub7W67Vbz3jk0OZCtgXSDUqTU5FXk45btcy76N747Jth9t1GkHKW5S4Dk0kRI3LLi/ud+mpp99IW1yUy3GOlifzwzbdKqArjbHcQKbkckum6+kfYw4EEDbD9KdV0TqZGwyIZm9Oj12G66yALBrm6PreP1Y9x0U9/UT7nFDAdngBAz2FnTFEhdkCEpDWqhX0YT04sc+doYy+ysTAesSVRjGpONjbaeVcvj+t0WGutpdb6y5Sn0Fri8U7B+rjRnACrdDFYXZjvZ2jxdEzVaaDZm/uUrM34lXwkS2+1BcmNYASA/Q0Xz8rX5XhhcKQ1yJYhjo4f5mSLZGHOdOGO8BbJr6hlrLdLzM7Ncdk2xftrst2QCXSmSS4GmimQD3L5XTH7WEOsKEPu/NLqeEbY3KprafKwxNolmH5JLUzGSoTDIlzeBs2+5rjB21zZmP9ONFsfUHGOrqfa708Nk43XbQHmsXDeWyhNG1ep12I9MV6K23jZsN8ZC710EtaK9aDbPdcMdeqLXpOO1ZTTxr62ukB3GuOk8HhIJ9vfbnctN+rOHLsVjjg/FVdKYj/l6mOKdkcb+55V7c1LZfQJAEfnNXSSiYgaBG09ucf0rZT7l3blZsBBF3xym3jyWCeyMFoI8hfdIyzO0YkSk/ojONHQgvTQnd3ZeGlFNynxQ5OeVwcvD/zPLqjW1Wq7jzGQC5GtFfqbQVAujbDIDJJfCMBYmOlXTh6owegKk9ksmuPntT+Ppx2+WeaiUy2aCnfktK2ZdK9/rSCZABsjWXGQkc//tZ/s4Mj8zYu/SksxcDMxntcvb9Kx8RFtEnuvmu6+BRnW5PibaN9C683vVFEtSpit8MZPdiIVoOVAl8Q5s7p2r7vEgcghqfU1dFyH3mDTU7aLtvB+iCb7YdjHh47zl0HDxm2xzHca3Y/kmM1z5n9quFKmY6LF7uucCDVIGCwAAuYkIb2+YO3537qsYPU43xkC7R/q7/YXNcvlaTBcDKY2xe/JWDwRs7tgd9KYrY8S2l4J8DWTn9jqQovivCYoBg2o1udDTBIXRMdcPo0f5dTV6mFoDlSjkdgo9M9NKyLUM6/Mtbly7xLO3znDPzVu4vVkhzxNkWYLVeovLzQrr9RaSJUhTcwKfrzLc3qywTjOcrTKskhyvfv5f4o++/BI8d/sc6/UWm80KaZpjs0mRrDNsL1PkKlAVJGdb5FmCs7MtNpvUPFhKctw43+Bym+LG+SUuNiuIFcQ7zjb48jduYn1tg9u3zpDnCfKtPY5PraErhVyYWaN0pbh4vmD1nDXUl8aUm+4dgtVzxiwnW4VkYoyy1Wr3xFCuOe2u6boL9hyg4+o9VVRgV8ddECTxfm96qig1Dc+l/N0N8rddQ6s675ngPDcDILNm3a40+a4xPl2RZn8bTcv3on29SMzxSCJ7aB5t63dFHPzPhQku7578iIIfeSjvuqtdHowourtdry+T144rIu18RmoGgLjocSGyhcFFKYYojbf5bAQVaqIgEBs5XqEi0u7xYW4Fs8jTlUNKMwxniFMrsE5kcwDbpOxndZ5hfX2D69c2OFtlEFEkViGyPLH9/e2JK4rUGjR/7scsN8KZq0AB3Lz/adx38xncXN/GF772AqgKVIHN5Qrrs60R2UQhYv9gBHaV5khEkSY5UlsOJ7KbLMXFdoXLbYrMXgCyrfmfX6ZmgE0m1jjbKIgaU5xsBMnGGmL7OdmWImtGt6MWYUYpvCv3XUrBdI8v1aXb5pBb0cu1uDGQRKGSFNEq2aonsrWbS1VPSO1vbqJ+Z4KLNuk1Sndc7DKS58bsqgLwRNZ/DF88mk9Kg+wJpSaC/J47kHzzWc8ke2IZiCBrmnT3laszVGyPaQyvMkPruc0U+7/X/9dmjDC/lRpdjt1AZZCcCpDeym10Uuzcx+bteE67q/2JUc3fRai9p3z5SkpD6wUnXITPdW/wuzUAguTcfN9eM+vk56bd52uzTH5mz/fEnM+6sif5WQ5sEnP+JsDq2hYiirvuuI3L7Qpnqy1WaY4b6w2+desa7r3zWXz9ueu45/ptPHDHN/F9d38RX93exF+59hg+8txLAAAX+Qrfdf1xfPq5F+FVNx7Dp597EQAgQ4JX3XgMn3z2Abz02tfwqWdehAevfwNfvrgLP/nSj+Dzz96H+699C48893x8x/Wv49Hb9+AVN76Kzz1zH66nG3xrcw3fdfMJ/OXte/CqOx/D55/9dpwnW5wnG3zfzYfx8MV9eOjGl/CxWy9BIjlSKP7G9T/H+5/5Lty7ehofeeY78Oz2HLezFZ5//iw+9837AABPPnUnLi9X0E0CTRQXt1ZALkhuJUgujAneXldc+2pitHsjSC5hxvoorH5Lqd+Z0d5Sx+13e411/alhTa7TcUns9GTuOp8DukogWzPWaHUrM8fQBTvE6b7re24HIDqT7PyVaqnPzgTXzx3XPW6Tl6Y6y41B9q8Nvn7bwbOFQfaCgPmd15E8d2FmZYFr/wEN989Dt43Qb0cMSERijnW/nR774tR0QPzfpOygXowCdsvWjIJsc+jafs7yss9apmXW7nFabueqdX3OrKEtipJVR8xWy1TeteYo71LV3ZXmgGalwAb7+5pvO+nJVqtR50o02PVlLkW6bvCLbhDWJLs+ysmlNYp2nzRdYXvHGZ65psjP8+Ii1Ejlkbx4V4oqz3zzBpJ1juypNYrHUplga8tRXuzKuuzqgVOJHNuocZoDqwxIttYUusdszuyGIsdbFzVGGXFQX0g9czxS5Fi8c62cHs5fRsvqbYoc19OHGM/adnbWd+2mdhDkclsts09g++JGTw/RlaGR44Ztk5Fpixq1LS/Vi3Bhkv3j7P/PFaZfgBWyvDQeivImUlz/RtclISkjZKKmC4DRZO/8c3oDVHTKnUv5GmZOZFV7w28jt6m5oTc6W878UOyK/Z+Kux6Zba/sYLL8QqqR4yI6DUDM+BXXDa4oYgLkqxUkE3zz+jVAgGdTcwOPHFg9m+Cpa/cAOfD1azn+7I578ZG7HsTFxQp333kbF9vywvWHq1fi9maF96//kyJdVfD/rV+O25sVVskrcLFZYbV6MbIswR+nL8U2S/CKe7+G//ZF78evfvn78Kkv34+PrEz+AKB5gs+c34csS/ChswexsRHiJFH8v+evwGab4v86++u42KxsM1BcW/9NfOWbd+L8fIuL22tkWQLd2sekT5eR42QLpFtTj+tLVCLHyRZYPyVIL213uY2Wuu6NL9nR8eK78wLop+OuGatCNbdPE/NqH2bXnr1l99Jwd76EorchrQvN1QyUHsf9dJ4Cz1bP30YN73NOH5FIzPGeTF1hoQPmm+TMb0g1AQZMBHGTmanW1De3th+batEXuWh/OYBVaaqKkbaZltuzAliZMscKcSJ2cMXWCHbq9zkWs4Jvjht3XcpClcZXSwMsQCX6HfgfTIM1hxszSMU9YslTQXYNyM4TM5ikYrS9ggVNfRtWIP2morX/NcKD+bzv6osfSlF0ArgT/fUH53lCGfq9h2Ae0ue4um8dwun/FliusmyTyLrIQhGFaMg78N0XUXnmdnidPnkNNchjLnfSNNx8uOjr0JuTWjeb9k135Bu4yZHMz1d3gxuAEZY0tVqb2yUTk523GNz0bTCaLaKmK1si5imTpxv+fL6AvUF1fUrtVGVm2jctusQlG/PUqXijuTdOw0WeK32Oc5jgg4cbZJetYaPFNt1eX0yf5HIuY3Xpa5gbhASApsjPYAa3CZDeMv1uJRPTxzpZIbu+QnZ2Dek58Ex+p4lO277Jm0yQrxWbrZj0LYAExfeN7SOc2fRsa5b/7Ofvws/cfAlWTyeQHNhuBJIqkkyQp4p8cwO6Vmw2UsxAsVkpnr4069/eSBFl10Rx60KQpkC+BdYwgXLA9TlG2efYDtpLL62e5rYF2EF/ro+x67vsBvNV+hyr97Q3q2puZ5/jWpeDSleFNgPc0Oc4qM/ue717Q1Irwyo1UeNEdudw7vpsv6dffxayzXbSgwzV7CHBjHAGrb/GY46jHVUeKFdSOyCekBeDimDMo2xsw8hN9wQBipd3GGH1pqMRgcBEDlRK05xeZLsnj7tBc4YmywvTnEjzbBWmXN6+NFyEdoxpce2Q8nPSkC7AjomV6hQ+ziC7/crXgvw5s0C+kmp/aK+IjabY3Z+Efq8dwuAgPP9wVk7S2joa+Ky+cUVplr2Bd33MbmXbuZ+veturCSMQFMVKXm2iEzK6Tcs2RQ5Cy4bSW8xv7zz6pPtEqysnQnB0OKo3R0W6hL8X6Vqmu4ty6/JoHmjTZNq9dKk41tpnTWEccI7i5RbedGuSe4OR7W/WS5bZZN55XWyzqi3FPMM2b03E3Es2zFaxs5tutgpP+9xTPtF8p/ueP1tFocliX1Rio83qbc/PP9kaTUvdbBJqdDtbmwHebhxLMYsFUJnFQjyP5GaxqOMG0xnsNVD9dFenWuSferNb5Kna2SqMES9uIOz3vrNVwDO2odkq/Kd7rk2Ws1iUxx0AhsxWUV/GFL68Jkw2WwVQvlTMlcn91sPcNmp5Vnfg4fV7/da0vK8FQ2ar6CAOcxxqEL3W8+4cjjVARqQa/fLTfRKBXOamcfiRFNcnKNuafkQV0bWmJxckl3kR2ag//vbnmwTsybrJyj7DRT9nYGfKuK5yA6XxhmdQgaoRrnyv7nclTWRnGcm1vKN27fqyGgnJ10kxGX1RlvosBjv70vJbQ9PoPSeyhn5HVSCLtFIgi3mO4S0bNMNaydMtf8g8x7uCGaiEkQxp73mO900H+gkcI7rHo0m3Q/oYSnffD0lv2paPFy2u0NAXWS4uock5ZAMbQTUvAhEVo9dF9LgM4Uou0EyK7mbFuV8rXyWKmOXGHPta7fTTznMcehFH6zzHHm42o3owZCfwoSjHwdTSi3mOFXYAsRb6bca3CNaeMd9eT6rar7mNansvJREU0VLduS5pMf9xdibFNHi78xnb6HQtim7qNZBuzb2Llpfd3qzGun10mj7iPMcVf+JHc+16xX83/qYePc7zwI1dQ7vvq8M9DGrveY67fuvS5H0CGPvqjl+PHZuNwxwDlQ7bs9F0EPv0ba6/ecmLIONyYx7VqX0ZCFzjy816W9tJP/GiBahNqOCdZOKn5YBsMkieW+F2K5uIRjEY0A76aO0E77KtJ3nC1PiyDvGXqa3vr2NP/rJPaJlPOWekvWFYJUHD2/R2vKF0vRCk+mYjfxnPyHrfG/vuumXqZreybsDw7vzWIKr1NFf+oYb1kGX7LtclhGNvL1iGHAe95W7Mt8ItnWPcjPjHy33vwn/DXtPx8jXaRwRysYGuUoimKF7kYd9wKoCdHSgrli/6HouY2ERDOy+7xqkxxraMFa0OBDh8Rn1DnjPWgXzq6cXjfxddTb1AkdXs9FY5tV3VJGv1f7Ez1oD5Mmaj5Rf3pFjfypHermmceDMzeVqsNqpbeaOeNbjl8p4xd10jExT70dSP173pr2J4gWY9b3lDXmmUG7R7z8BE5xvy+gQRd64hHWXyGeN6M8X6e+bVaY5F5J0A/h6AJ1X1r9m05wH4TQAvBfAIgJ9U1W/Y334OwJthgvT/g6r+bo8STyqyO28KG55B688iEjb3Nt1pDvIckiRlxMGNnrPmtYgcu826u1C/DLW+RZIpsNmW312UGihGXwsATRJjwHf6Ge2WeedyURPjHXEO5BGsizo1kam8ZS8tIzZHIXBiB/vn1j93dWHoM4gtb/jNu4PvLEeojQ65I5/i/Ns3z0Ney9x3m4fcjMdwI9/B0XT70Fdo97lgA8Pr3F++bd3gk7MEuNwY3VSFrFLTLS4BpOjO5plZuMHUwXt5u7zdnNM8NU/7Ko+ExUSmARQmuaLndWOJPbS4KYgRMOB+f+Ridc8oArs3AeblKYmXD6pGGaVh3+mLXUtfP7UJpvsGtjCE9YhgcdEty13vkhOM7ALlsapFbJ2Gix+JrFzDGtLdd1eWNj1vOxe6fEguzcsNjc4O1e45tP5g2svcJ3L8LgC/BOBXvbS3Anifqv4zEXmr/f6zIvIqAG8E8FcBvAjAH4jId6pqu7LpCAZ2jL6F2vNASVWBtKNfslxuzAmxXpn99CMhiTWjua8S5nNdaI1g5eXd5zYznd3d99r61b7MaBXKosw7+7qbJoG0xmhunyCbM8a+gNhIRGhbQfou16ed1ZepP+bqa1D9fluB5YIi2WjCW7YTFMMJRKd+g3ZEDtaHq8e7MLVuA5V2oKrmvJ3qwlrvb3wISbOR8G9k5Wxt9idJIJJ7XSmkjCK7SG/TKeffJLtpDbO8GKikiZhbEmnQ66ap6IDKckM1eEdbQ5FdeEbZS0+cXvtGMQFkI0V/6cAGq/XdEEkOfg+tB1S7InjBn0owomKcw9pe37/iOAHNprfJ8NaXq283NNNDSNOH6mzf5f0nKqHP+zLx9cHpS/173/ShdJpjVX2/iLy0lvwGAD9oP78bwB8C+Fmb/h5VvQDwZyLyMIDXAvijzpL06k848gV/3wNYv2aEDoA/+CA3hlhcFMI8d4Ou0t3+MHUS/4T171JzyGZbdngPCKkEDHfT90Kg6/SMAgfNd4iQAa9Pu+VuDrbSPOAmRG37+fU1nn3Jnbj52W+Uib3EI7BMmyHtEsmW5Xs/WvOFK7APgwwJB6mdPEfRbQW0NuAmGCg4hLrmjxW0r+fjBTyKfbAXV9HUTKOWSBnYSLwnbE06Wz/3rc7JNqtEs6VmjOt5Nc7TXPssTcv4NOlp0/LOqIeMRp5XI9+ul8nWjxxbLfdvotquE6ngW3/lbtx44hLrrzy3W5yazlWertpyFvtTN6hFucOaXLmGueh+w7LBz320PGRAe2p3XeOH3IgGl/XLMmYwZcjNrk+PG9+6vrjvrelN+XYUc98+x/ep6uMAoKqPi8gLbfoDAP7YW+5Rm7aDiLwFwFsA4FpyZ7PxHfsOZBJj4N89B058zUwfr1whui4m1JZLI7TBt+wV65qTXLK8fISU52YEaf1xYZKU0YuBUeKKQNfpMr0dv1entwusHhIgkWbD3oUtT7LZ4uYnbu3m3SeLPjMptBnmtvUHGN9dQQuVY6Cw7XtOHXAXPip+dCkUrdgnclHPz087HY6n24ewT/t0UeV6GlBNDy0HVAMe7rhnVc0WOzgPVotdv+DKOqG26Rstp92hdtag0YU2+3mHoq4N3S/8ly4EjXYgWryzjE9IyzJvWW+UszPWjSY+wD0f/kq/c7hu2OvR4lpZG+fXDUWra+sGDbX/vW4uu56o9PUhLefXkLMkuKwkZf61p+E76VOc53Wm6q3WmG97DY49IC/U7oMlUNV3AHgHANy9foE2nghjm9mpD3KG3YbmyLcwk3mvgTQtBKCY4iZNq8vX7nwlh4kUb7fQwJ3ezgCzHbH1B7W0LNeVXuTXwzyI7M4IUafh2Fci38ekLdrQtdyYxjfUVtsuGJNHhGOKONfLoi2/9SBBIJIxPJuFMr5uD2WMtlsYmEC02RkBSXpsyzPQmkEvcki2AlZ2IuBEIFp2rwB2o7s7XaZyM3ORbmw/2oboLwBvsGjgt9DyxXrhdF9HW7W4j1GuE7rm+Pn0/d+Wf1P3iyHl6hVk6BnACPzepeXBHkpT3gT2wS9TUw+qHj2rTpV9zfETInK/jT7cD+BJm/4ogBd7yz0I4LHO3BTtgjWFoZ3SSCQt5d3m5kQSgaxWpUkWAXI3s3gt0gFYcc3MY8x6fVgzrvUIduAEb+yD0yZQbaO9+6x/LELRhkP6PzXesHWb1uDjrn3Nb+Mo4DkHM5wQvmHyIyan1/3kuLo9hDHbcnDe08xqYs/tVIIcCs0uge3W02w7zqOIyjYYNxctdtpd5liWsf64PKQBjYZ4YMCjb75F/rVH0qG5pusGNmTehhq3LvaZlaEp/UDje7CO76zb0Uab/OpQH+tHhZsCeqdMx+HY1xy/F8CbAPwz+/93vPRfF5FfhBnY8UoAH+iV41BxHPHCNcZgn4rp7CybHYjhJrpP0lJk/Ttvv2zbbZn3zpyP4UeC1QIGDDTQOjClWLVv/YxhkIdGLsamx772Mr7AsMjv0O4Se7T/Pu18jIEMURO68QzdzJ7mjcf8ut2HsbS9ns/gfLNqe3GBje3WTM2ZJO3T+mWZWT5r6EZhN1GmNwQ5gOana23ndN9IsBuMFdqX+uHNG9KzDMWbMENN4linU1Pf2YZ62lvLD7kG9FnvKGQNnz3GvqlZEH2mcvsNmEEc94rIowD+KYy4/paIvBnAfwTwDwBAVT8lIr8F4NMAtgB+qteIZ+jBDWTu0exDOsZXohkZgKSsIgWqfYFCtHXd8B8JVpIDBtrltVPIHrNxAAGTPlYUKYKIXed8vAP7yA+N/rZsf9S2Xrsoav2mLGKKaa+GmPlQvZ5ehPhIuo3J626udli0px2DvTW6t92at5imSWlAnW66+Y8b66bhMXjfIIcX6WucKQkotShQh5WbYDd25dApCmee4rBztpShGnyoAe5xbsSssyICzVqeNp84EsPBuTu9V7/3zh+fdiNzzqc34KUBvRtir/6+Hds9YCaI/uv13HdfYPqsM0dkb59+8XtGfzvPy1nnh5yAhb9Y4/ee+j8+pKoPzV2OY3L36l79W3e+ofF3vw3vNcVbF2OdAwPb3iCz0KaxfXSuS6OHlGWsR+cjv6Z3knzadHdMHe+aXSF2Ha+/DGnhOjzkyecfPfM7+Nb2q40LxvOGvL5M3Zia+iX51H/vOhBdZQ5E7po33RDR2MmzR/+6YoRxnxOiIdrRudoe0YRYBgEcOr3gsYQzghvcIEMG+Mx9kfBZ+AXiaGj/yFfUJqHvG/csfbsnmbzbls269bT1KSGAUBeMxkV7TEHah7H6vIZQxe4rVjGOxu37NPAQHd+3XY+l6UO9ydw6fKD2DnryOVGf4/GZ8qUFY6/XNRK2iaaGOod5LvIeMEgF6CHWC2TfKPQYj836tPvRZgTouZ9tfRCH0PJYd1Tqo9gPfQw49wViSRy7ro7ZlnwGaDRw5CCHIx96vjZEf7tmhaiPe9ljvtrB5TqUA8dvjGKA5whiHHObY3S/iEh74zHHTcQYFdvXOHTtyxHMc7mpgSa6sr14GvBgRu4j2ftR8dgCOqWIuLynFqqxIrT73qwO4Yr2u5uUGLW9T5lCbWGq6PMgveoRiW7doJ2hJa0Z35CxbIsSR/Lgr+BYQQyTUb/lhuY7JVPp8JQcQY/jMcdTV+zYjXCf/Loa4VTmObD9of3/hs3GER+T962fIgI8gXAOeiQ8NVNcGKbqEhGjkYuBueql3nbqfSeHMqTdTKnThwQ6gMO02QU9MniD+vJwZHiM7R2JwQGMJCkH+4YGGUag46OyJB12HEF34jHHXcTewPowUBx3mLr7RksZYhi4OSljta8jC+cUx2XKYz258V6i0JPDj9sx1z80yAHMGugIF6c2i5JWo9ZuG6Fo9uKuDV11nOdQ17UMmPxp3uwzskzBlH7tSHochzlWjcb8HjrC+uAGN7BvW5BDxLlPGcgwE3xAfS7uwtPB6HOKHwOeD+1ctfo5NMgBHC/Q0bM8O+dlFi5fr/N3Ke2hbV9CEeMT0/Gxy7ToJ44B4jDHE3Howd9n/ckfW48hzMB+jyWO1fgjFJJe7HHSTiaak74B0uuvPsMk8YfU2VWds3MSZghqHMtkRBHkAA430X3LEwtNA2q7Zo/qu2wXC4v+xsTiAh8d5Y3SHIcq2Z/wf+ns1ZdsCAd0n2jlBOp+ECNeUMaf43XGY+Fv+6A+jssy1iFotg8jVj2PQqPH7Ae9JIbuy5Gmi5ykrcbUZzuCt+HFZLCjMMeK7kqJVUSnJIoodJ1D+/ssJYKxJwe305jEcmoO3dcTEfOl0ke3T51JNXqqIMcx8V8ucax+qLGZ4KVo+gnoMTBg7vWO36Mwx7MydcOduMEMOZFHuaNauLmN5mK+FMGMmSF1GIlwkyPi2kdodoUjtocogxzHZuoyBvKffoYianiFfesjUm0+XXMcS8MdWo4JG0pU03h5RGNYp2CudnjM12sv4aUwe83nHadoLxlfX9pfrKDjvY44lEdkF/LJu3HMxEloeyxewlHXdl9/3W9L0GRHZOeiI25z3PT+9r7LLpEhb7abgJMQs77E3mZ8EZTEfA8J4dyMUY4YxbzP28vIOIxpYKdkIfrcNNNSk7k+ad2PsR05xtDOUB6H5BujFoeYOOARjznuu6Ntj8quEn32+ZQu3qd6jPuKmFsuFkM8Nn32KzbRPtU2OYSZzWK0zKzPTWY3ShPsz3oz5hODOTgFfd5nH0LBmxgZ0K7iMcdDWerJc0z2Eegppuc61WPV9njLX8aln4JwzsnQ+otdqK8SVleC5uxU9aGLJfWZn/oYjTUDztRQw8PsG7yJWKOXa44Poe0O5xiNP6YGsZTHmcdk3zbQtB4FdR72jYCQaajfeB9DZ0JdkZYS5fIZY8pEPyrrp50a1NvlEPETw9M0x30qfM4TiBftYUx5rOr1SmG92kQs1ifBUDM2RZ/MMbsoLaEthKKyMZhi/0bFp+kJXNsy1O3Tpe+xHflcjMQcKxt3F0u+aDdFaWI45jGUoQVtuIhJItBcIQMetzbldSh+GfxtDCnbooi8zURD18j5U6jHJevyIYx17JoGk3XV2Sm0nT1p03GnufVlTlaLfYa0iR7nZCTmmIxCU+MYIs6HTgVzIl0LmkzeVAZzKK4cMZSnqQyh9Cbx9n8jS6FnUGNh5/7ozBT52psYjpd/HRq5PC6oEIN2+kGOMcozRIu72CcAsxg0R9drQGiOrwL7iEsMAjkB+4hEDCJ6KrTVZVuUnJCT51iaO4HhnAxbzpCB7GMom27GY9H0mIIcdYaW7dR0muaYREuMgkGOzz6PEf3fCCEeERjjMYIUffLgNeR49KnrJWkyzTGZBYoWGYN9ItE+SxJrQmIkpm4KJG6WpMk0xxhm1IYeuCU1hrE5abEcM/oSS1/DK06fCDUhVw1GaMkx6dtVZmpOzhxPfZJOkf8pjPA/eXGc8lHk0gbsnBAn325PlGNHK0PbWqpWN3Gy58Lc3Uio26NyrIBjHOZYT/jEHMgU9RDDdF+zMrc4jsUYs5GQvTjJ82ImxjK2xzwmoW1N+cTxGCyyTQ+dOi8G7e8zjzMZlTGiz3GYYzIpixTBPsQgfDHg1wOFlxyLA4IaJ6tJDVy1/d2LMedOjnlGDprl2dFcu2ZyozkmERKrqC2BQ+epJoSQsZhLy5d0DWFwI0pojklcLEnUYmYqk9x1fCjuhFxNqN2H0+ftgOQo0ByT40IBPS5T13f98eXIr/AkhEQE9Xt62gIbdfPc9CZBt0zb8aL+tkJzTMaHAjoM7dEfUeIb0ANgv2O9z4AZCjkhx4U6bnD6XNfgpvTRtttQ/6H0etqYr3X3jfcV0mGaY0JC9DGsx0Q1XoM8lKEX3SGCzMeShOzHKZth38iOoe1+Hl35LV23hz4ZPBH9pTkmh3Nqojq1Md5XoEPrLF14+3Bq7YuQWDilc6uPpu6r7XWdHaLh9eVOXbNPxEDTHJNuhvZlio0Yo8Bj53XqgnsIJyLWhIzKkjTcEZuW74O/D1dVtxfQF5rmeAnE9KripQnqKYhpH65adKKNpbVRQo7FUs6N2HS73q1trPIxuLFLJMa5c0si8k4ReVJEPumlvV1E/lJEPmr/Xu/99nMi8rCIfE5EfnSqgs+G5sf7i6X8MaLa74+QLiKJVIwJdXsiQrp4zGtCjNeUQ1mKbk9Zvpj3OyaO2Lb7RI7fBeCXAPxqLf1/VdX/xU8QkVcBeCOAvwrgRQD+QES+U1Wzzq3EeuLO2Y0g1joh8cNHd/05zfPsXViKbve9Oalvq77esY/jababaaEBbOeUBl4fm5HPx05VUtX3A/h6z/zeAOA9qnqhqn8G4GEArz2gfPMT8x33qcKIMCEHsSjd3jfquYSo6FWHuj0c1lUUHPI88adF5OP28d232bQHAPyFt8yjNo0QQsj8ULfJdNAMjwfrb1b2Nce/DOAVAF4N4HEAv2DTQ88DgkdYRN4iIh8UkQ9ucLFnMQghnfBCRQzUbTI+NMPkBNnLHKvqE6qaqWoO4FdQPoJ7FMCLvUUfBPBYQx7vUNWHVPWhNc73KQZZKuwyQcjRoW6T0aBeHwfW72zsZY5F5H7v608AcCOi3wvgjSJyLiIvA/BKAB84rIiEkFGg0F5pqNvkYGiIyRWhc7YKEfkNAD8I4F4ReRTAPwXwgyLyaphHb48A+O8AQFU/JSK/BeDTALYAfqrXiGdCCCGjQd0mo0EzPC+cwWIWRCNo+HfJ8/R75Ie6F2QDWRYRtC1S4yqeQ0doh3+g/+eHVPWhyTcUEUHdvort65ShhsfBVT2v/JekjNwW/0Tfh6f0640Vu6w35HVVzlVtQIScGrwoL5O+x41aHTc8/8hU7NO2ZmiPyzLHXQytQAo0IceHF15CEx0nPDfJEE64vZyWOR5K04GlIBNiOGHxIwugrf25R63Ua0LGhbp/xc1xE12CTMLwhIofHiNyKri23KdNU7e7oTacDjyWB0NzPBQa52Ym6DRPCCEHQwPdDfU7Xnhcjg7N8ZiM1YCXLNIUWELIEtlHt5as1SGo34QAoDmOkzHEaU7RpsASQq4C+xrq0HqxGG3qNyE0xyfLlOLWJeIUVkIICdOkj4fq5pjmmgaZXHFojslwKJqEEBIXdV3e1yxT3wmhOSaEEEJODppcQvYmmbsAhBBCCCGExALNMSGEEEIIIRaaY0IIIYQQQiw0x4QQQgghhFhojgkhhBBCCLHQHBNCCCGEEGKhOSaEEEIIIcRCc0wIIYQQQoiF5pgQQgghhBALzTEhhBBCCCEW0QheMSkiXwHwLICvzl2WHtwLlnNsllJWlnNcllJOoLusL1HVFxyrMDEgIk8D+Nzc5ejBKbWzWGA5x2Up5QSWU9aDNDsKcwwAIvJBVX1o7nJ0wXKOz1LKynKOy1LKCSyrrMdiKXWylHICyykryzkuSyknsJyyHlpOdqsghBBCCCHEQnNMCCGEEEKIJSZz/I65C9ATlnN8llJWlnNcllJOYFllPRZLqZOllBNYTllZznFZSjmB5ZT1oHJG0+eYEEIIIYSQuYkpckwIIYQQQsiszG6OReR1IvI5EXlYRN46d3l8ROQREfmEiHxURD5o054nIr8vIl+w/79tprK9U0SeFJFPemmNZRORn7N1/DkR+dGZy/l2EflLW68fFZHXR1DOF4vIfxCRz4jIp0TkH9n0qOq0pZwx1uk1EfmAiHzMlvV/sumx1WlTOaOr0xiIWbOBeHWbmj16OReh2R1ljapeqdkeqjrbH4AUwBcBvBzAGYCPAXjVnGWqle8RAPfW0v45gLfaz28F8PMzle0HALwGwCe7ygbgVbZuzwG8zNZ5OmM53w7gHweWnbOc9wN4jf18E8DnbXmiqtOWcsZYpwLgTvt5DeBPAHxvhHXaVM7o6nTuv9g125YxSt2mZo9ezkVodkdZo6pXanb5N3fk+LUAHlbVL6nqJYD3AHjDzGXq4g0A3m0/vxvA35+jEKr6fgBfryU3le0NAN6jqheq+mcAHoap+7nK2cSc5XxcVT9sPz8N4DMAHkBkddpSzibmrFNV1Wfs17X9U8RXp03lbGK2Oo2AJWo2EIFuU7PHZSma3VHWJmLTwqjq9BiaPbc5fgDAX3jfH0V7gzk2CuD3RORDIvIWm3afqj4OmAYP4IWzlW6XprLFWM8/LSIft4/w3COaKMopIi8F8N0wd6PR1mmtnECEdSoiqYh8FMCTAH5fVaOs04ZyAhHW6cwsYd+XpNvRnQstRHsuLEWzgfh1m5ptmNscSyAtpukzvl9VXwPgxwD8lIj8wNwF2pPY6vmXAbwCwKsBPA7gF2z67OUUkTsB/BsAP6OqT7UtGkg7WlkD5YyyTlU1U9VXA3gQwGtF5K+1LD5bWRvKGWWdzswS9v0UdDu2eo72XFiKZgPL0G1qtmFuc/wogBd73x8E8NhMZdlBVR+z/58E8NswYfgnROR+ALD/n5yvhDs0lS2qelbVJ2zDzgH8CsrHG7OWU0TWMML1a6r6b21ydHUaKmesdepQ1W8C+EMAr0OEderwyxl7nc5E9Pu+MN2O9lzwifVcWIpmN5U11nq1ZfsmrrBmz22O/xTAK0XkZSJyBuCNAN47c5kAACJyh4jcdJ8B/AiAT8KU7012sTcB+J15ShikqWzvBfBGETkXkZcBeCWAD8xQPgDFyeX4CZh6BWYsp4gIgH8J4DOq+oveT1HVaVM5I63TF4jIPfbzdQA/DOCziK9Og+WMsU4jIFrNBhap21GdC03EeC4sRbPbyhpbvVKzPfQIIzXb/gC8Hmbk5hcBvG3u8njlejnM6MaPAfiUKxuA5wN4H4Av2P/Pm6l8vwHz2GADc1f05rayAXibrePPAfixmcv5rwB8AsDHbaO9P4Jy/m2YxywfB/BR+/f62Oq0pZwx1ulfB/ARW6ZPAvgnNj22Om0qZ3R1GsNfrJptyxatblOzRy/nIjS7o6xR1Ss1u/zjG/IIIYQQQgixzN2tghBCCCGEkGigOSaEEEIIIcRCc0wIIYQQQoiF5pgQQgghhBALzTEhhBBCCCEWmmNCCCGEEEIsNMeEEEIIIYRYaI4JIYQQQgix0BwTQgghhBBioTkmhBBCCCHEQnNMCCGEEEKIheaYEEIIIYQQC80xIYQQQgghFppjQgghhBBCLDTHhBBCCCGEWGiOCSGEEEIIsdAcE0IIIYQQYqE5JoQQQgghxEJzTAghhBBCiIXmmBBCCCGEEAvNMSGEEEIIIRaaY0IIIYQQQiw0x4QQQgghhFhojgkhhBBCCLHQHBNCCCGEEGKhOSaEEEIIIcRCc0wIIYQQQoiF5pgQQgghhBALzTEhhBBCCCEWmmNCCCGEEEIsNMeEEEIIIYRYaI4JIYQQQgix0BwTQgghhBBiWc1dAAC4V75dL3FpvoiYf+5HqXyrLNP8e1Oa/VD7KZzWlQ+gTesFVgttQ4PLNa0PQMSs07pMYBsDtjN0+cp6A9dx6dqQPmj7Q7fd+JvuvR87e9LUPMT929nzStN2+TU3Xd1Zx8/Tb7476cXvofTdfG1r382rMX9v+Z1tNOTlfqvtT+uyA7ZR2W/4dalFXQlq+lNZXmppZcqHPn7xu6r6Olwh7pX79RIX5ot4tdOl2aFlujR75+cDNDuwSGjZSTQ7tBw1+8C8YtJsk+cQzfbz3T0twtrldLtNkyvl8PPqodlmOW3+DbXf9tDsPtvY1erasp5m+8uNodlRmONLXOJ7kr8LSQQQE8wuPidihDexQW4RIEnsf4G4Zdxv/u/ipdt8EFre+13dOkn191C6ipSxd/dd7OdKulsHlXW0SPOWA0weiVvPExGppdXW0aT6m9mGn4bq8kVaeQa7tK51dtLRkB7Ib3f74fSiPPW8d9bRxnKZP60KcbG+lsugXHbnd1edleW1WEdsugSWF1HbXPzf1DaT8ncASIrvapukFssnKNOcoLjPSfG7FuskLt3+1p2eV9ISyZG636DF76n97PJKRZHA/ean50hg/9vl3brl/7zYfoq8uo7NK7X5uuXLMpnlU5TrpKj/ppVyld/dfqBMB5AKkELsZ7F5iU0TpCJI7MlrvifF9/T+L9yLK8YlLvC9qx8xXyTZ0WyT7LS61GwApW7XNdv8uKvZQHV5L6+KFvu/AzvpdR33NRtwGlrVbD+fumab37Cj2WV6h44nLb/t5FXTbNSW61hnJx0N6Q3LN10TdvS4YZ1qurZup1OzUV22Wianx97yRV5VrS5uqr3lfU32l/M12x6qHc0GPF1GQMdrmu2n17XZ5dWcXtVsk55bTa5qtvutrtkACh329bfIK6DZbvt1zQac9lY1G0CxTV+zTbrNu6bZrlx1zTbb3tVss21Yja5qtvktGazZ7FZBCCGEEEKIheaYEEIIIYQQC80xIYQQQgghFppjQgghhBBCLDTHhBBCCCGEWGiOCSGEEEIIsdAcE0IIIYQQYqE5JoQQQgghxEJzTAghhBBCiIXmmBBCCCGEEAvNMSGEEEIIIRaaY0IIIYQQQiw0x4QQQgghhFhEVecuA0Tk3wO4d+5yBLgXwFfnLkQHLON4LKGcLOM4jF3Gr6rq60bML3qo2wexhDICLOeYLKGMwNUpZ6tmR2GOY0VEPqiqD81djjZYxvFYQjlZxnFYQhnJfizh2C6hjADLOSZLKCPAcjrYrYIQQgghhBALzTEhhBBCCCEWmuN23jF3AXrAMo7HEsrJMo7DEspI9mMJx3YJZQRYzjFZQhkBlhMA+xwTQgghhBBSwMgxIYQQQgghlitrjkXkxSLyH0TkMyLyKRH5Rzb97SLylyLyUfv3em+dnxORh0XkcyLyo0cs6yMi8glbng/atOeJyO+LyBfs/2+bq5wi8p969fVREXlKRH5m7roUkXeKyJMi8kkvbXC9icjftPX/sIj8byIiE5fxfxaRz4rIx0Xkt0XkHpv+UhG55dXnv5ixjIOP7Qxl/E2vfI+IyEdt+iz1SA6Huj1q+aLUbbud6LW7pZzU7/HKOZ+Gq+qV/ANwP4DX2M83AXwewKsAvB3APw4s/yoAHwNwDuBlAL4IID1SWR8BcG8t7Z8DeKv9/FYAPz93Oe32UwBfBvCSuesSwA8AeA2ATx5SbwA+AOBvARAA/zeAH5u4jD8CYGU//7xXxpf6y9XyOXYZBx/bY5ex9vsvAPgnc9Yj/0Y5ztTtacoajW7bbUWv3S3lpH6PVM7a70fV8CsbOVbVx1X1w/bz0wA+A+CBllXeAOA9qnqhqn8G4GEAr52+pK3lebf9/G4Af99Ln7OcPwTgi6r65y3LHKWMqvp+AF8PbLt3vYnI/QDuUtU/UnPm/aq3ziRlVNXfU9Wt/frHAB5sy2OOMrYQTT06bOTgJwH8RlseU5eRHA51ezKi0W1gGdrdVE7q9/jlnEPDr6w59hGRlwL4bgB/YpN+2j4Seaf36OYBAH/hrfYo2kV5TBTA74nIh0TkLTbtPlV9HDAXDAAvjKCcAPBGVBtwbHU5tN4esJ/r6cfiv4G5+3W8TEQ+IiL/j4j8HZs2VxmHHNs56/HvAHhCVb/gpcVUj2QPqNujErtuA8vTboD6PRZH1/Arb45F5E4A/wbAz6jqUwB+GcArALwawOMwoXzAhOjrHGuqj+9X1dcA+DEAPyUiP9Cy7GzlFJEzAD8O4F/bpBjrsommMs1Zn28DsAXwazbpcQDfoarfDeB/BPDrInLXTGUcemznPOb/ENULf0z1SPaAuj0eC9dtIE7NoX6Py9E1/EqbYxFZwwjsr6nqvwUAVX1CVTNVzQH8CsrHRo8CeLG3+oMAHjtGOVX1Mfv/SQC/bcv0hH2E4B4lPDl3OWEuAh9W1SdseaOrSwyvt0dRfSx2lLKKyJsA/D0A/5V9PAT7qOtr9vOHYPqDfeccZdzj2M5VjysA/wWA33RpMdUjGQ51e3SWoNvAQrTblo/6PRJzafiVNce2D8u/BPAZVf1FL/1+b7GfAOBGTr4XwBtF5FxEXgbglTAdv6cu5x0ictN9huns/0lbnjfZxd4E4HfmLKelcncXW1162+5db/bx3dMi8r22zfzX3jqTICKvA/CzAH5cVZ/z0l8gIqn9/HJbxi/NVMZBx3aOMlp+GMBnVbV41BZTPZJhULcnYQm67bYftXYD1O8JmEfDdcTRhkv6A/C3YcLtHwfwUfv3egD/CsAnbPp7AdzvrfM2mDuUz+FIo9gBvBxm9OjHAHwKwNts+vMBvA/AF+z/581czhsAvgbgbi9t1rqEEfzHAWxg7ijfvE+9AXgIRjy+COCXYF+eM2EZH4bp9+Xa5b+wy/6Xtg18DMCHAfznM5Zx8LE9dhlt+rsA/Pe1ZWepR/6Ncpyp2+OWMzrdttuJXrtbykn9HqmcNv1dmEHD+YY8QgghhBBCLFe2WwUhhBBCCCF1aI4JIYQQQgix0BwTQgghhBBioTkmhBBCCCHEQnNMCCGEEEKIheaYEEIIIYQQC80xIYQQQgghFppjQgghhBBCLP8/BqWodEUC79cAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x432 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pylab as plt\n",
    "\n",
    "plt.figure(figsize=(12, 6))\n",
    "plt.subplot(1,2,1)\n",
    "plt.imshow(v_sql.reshape((180, 360)))\n",
    "plt.colorbar(orientation='horizontal')\n",
    "plt.subplot(1,2,2)\n",
    "plt.imshow(v_hpx.reshape((180, 360)))\n",
    "plt.colorbar(orientation='horizontal')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "remarkable-exemption",
   "metadata": {},
   "source": [
    "They look the same, but the two are actually different enough for concern:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 148,
   "id": "loaded-scratch",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7eff14d4e7d0>"
      ]
     },
     "execution_count": 148,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWgAAADsCAYAAABZszyvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnsklEQVR4nO2da8wk11nnf09Vdfc773jGHns8zmB7iYMcaZ1dMGEwrGBR2AAxXoTDB5D5ZO0ieRfBCnaFVs5GWtiVooWwLAitYBWWCHM1WUGUCLGQEHH5suTi4BgnxsnkAh7bseP4Ore3u6qe/VCXPl1vVXe/t+nqd/4/qdVVT5069ZxTdf516jmnus3dEUII0T+iVTsghBCiHQm0EEL0FAm0EEL0FAm0EEL0FAm0EEL0FAm0EEL0lAMTaDO728yeNLOzZvbgQR1HCCEOK3YQ86DNLAY+C3w3cA74OPDD7v6ZfT+YEEIcUg6qB30XcNbdv+DuY+Bh4N4DOpYQQhxKkgPK92bgqWD9HPAtnU5cu+mjU9cekCtCCHGwXDz75Rfc/cb9zvegBNpabDOxFDN7AHgAYHjqOP/4l//VAbkihBAHyyP3/Le/P4h8DyrEcQ64NVi/BXgmTODu73H3M+5+Jjm+eUBuCCHE+nJQAv1x4HYzu83MhsB9wAcP6FhCCHEoOZAQh7unZvbjwJ8CMfBed//0QRxLCCEOKwcVg8bd/xj444PKXwghDjsHJtDrTGQ7mxueuxGZk3vb2KgQQuwOCXSDNnG2OYLtpTiH+0qohRD7gQR6CXxJwZUwCyH2Ewl0A4msEKIv6NfshBCip0ighRCip0ighRCip0ighRCip0ighRCip0ighRCip0ighRCip0ighRCip0ighRCip0ighRCip0ighRCip0ighRCip0ighRCip0ighRCip0ighRCip0ighRCip0ighRCip0ighRCip0ighRCip0ighRCip0ighRCip0ighRCip0ighRCip0ighRCip0ighRCip0ighRCip0ighRCip0ighRCip0ighRCipyR72dnMvgS8BmRA6u5nzOx64PeB1wNfAn7I3V/am5tCiCtFZF4v524r9ETsRw/6O939Tnc/U64/CHzE3W8HPlKuCyHWCAtEWqyOgwhx3As8VC4/BLz9AI4hhDggcjeyPFLvuQfsVaAd+JCZPWJmD5S2m9z9WYDy+9QejyGEEFcle4pBA9/m7s+Y2Sngw2b2d8vuWAr6AwDDU8f36IYQQhw+9tSDdvdnyu/ngfcDdwHPmdlpgPL7+Y593+PuZ9z9THJ8cy9uiF0Qmc8MBgnRha6V1bFrgTazo2Z2rFoGvgd4HPggcH+Z7H7gA3t1UuwvkTlWfq5Ew3PFMg8FEukrz15CHDcB7zezKp/fdfc/MbOPA+8zsx8B/gH4wb27KfaT3K2+M1+JgSDNCDgcaNDwyrNrgXb3LwDf0GL/KvDWvTglDh41tv0lMj+0dXpYy7UO6E1CIYToKRJoIfYB9TLFQSCBFkKIniKBFkKIniKBFkKIniKBFkKIniKBFkKIniKBFqIFvf0o+oAEWggheooEWogWrrbX0/XE0E8k0EIcAOsmeFfbDWldkEALcQBI8MR+sNcf7Bf7hLuR5dNeV5pHuBtmzsYgXaFnQohVIYHeByoh3c1+aRYRRTl5HpG5ceH8BvlWDOPy4cbg/JGUeJAzHE3YHE2uWO+sKlezfF3l7Yt9kd99KY8Qi1CIYx/YaeOrhDnLpw13nMacf3mTfBxD5LCRYZspDHOixEnPD7j40hFevbAx09Pe6XF3miYUtbZ0oT1M1/ZdfRbl28yn+nTl37SHIjwvXdPe5d9O/G6rxy57W30KEaIe9BXE3ch92hjjqPhXk61JwsVXN2ArJjo64Q1f8wInRhc5ObpAhDPxiMdfPM0Lr1zDZCvhgo04trm1dG9tGVFoE7V59nD/RSI9z8ed5NuW/7KivV/2rvLtVLTn5S9ExVUp0FWDuJKNIcuNPI8w81qY3Y3Lk4QL5zdgHBFdM+Eb/tE5btl8mdRjtrKEI/GEUZTyjTc8zYXrhjz2ldO89OI1nDc4tnm5tTxNIVhWNNu2L9qvK5/9Shemnyfmy+Tftv9O89nr/k1bl+i33WQVPrn6uCoF+krjXogzFL3m3AE34qhomPk4ZnjdFt9w89OcGp3nQjripfERco8YxinHki2OxBOOxmNObl7kpa8eYzKJ61BHZNPj7ETE1o2qHLlPyxzSZl/m5rRq5on4sj32tvRty12x965lsVquSoG+khefu9WxZjOvQxyVOKdpDGnRiI4lW1zKBrwy2eByNgBgnMdczgbcMLrAseQyrz/2VT6bv4580j58sF8i1CZ2eVltlT1vVOM8e7hv23JX+i57uL1pr/yoRKnp97z924R/Xvp5fnfltRMWhYnC+Piy9p0cq+1pLLQvCqctk143g256PUh4UAMq+/GYulN7HPm25dwhzw2bRLzp9LMcTbZIooxxnpBEOUmUA8WUu1fGG+QeMYpS3vi1X8YnEXHknWKR++7t4Xdb+tyZGagMhbltADOMu4d5tfkRCupe7Flu27ZV/jVvPs00zfKHZWyWJ6Srnpp12HUO2vKYZ9+vcMsie7PX3eyVN/OoOiPNY4a2rrz3Wp559mU5qHx3Q+8Euu2ENlnmpM7Lf1Esr63n0Xastgsu/HYvZmuEaauGXgmImeOx88r4SJ1vYtl0uRTqyIrBwtRjcgwb5tv8aetdhiLbJVZNe5f4dD1ih/u3PXp3xVnbvnc7iBeGkaonlupmEYpIlb4aE+g6d2F5q3paxv82/5rnIqy3ZW6M1Q2lK32Yfxtd9mVZVrDaQi7hctvNoK2n31xeJpwTstexhoPQnd3SO4He6R11N+mXSdfWcNvsXRcXTBtGKBYwK85R5DDI+fIrxwAYRSmbyYTEMiKcCGcYpZwYXiK2Yv25166p82keq/pu87mrgezmuxrwDPOLrPsRuxLPeeJa3NDiOn1oz/OoDguFYlylb4p4msVFSKkhwmH68LvtuG1+T28Cs+nm1VPzHLTZ284bLBbjZUU+3K/Lvl8sK15dolvR7BiFN9lqfV44p+04exHfZfLZb3oXgw4rt2tAo83edlIWTeHq6kmHx2gT6bb1NioRq3o/1be7EUU5kUGSZESjjGEy7TUfTbaILCfNYyJzjg+K2RoRzjOXjnPhtQ02NsfTAUeKsEkzvBD6HkV5a08yTLusPRS1OMqZZDGxOXm5T5ZHxNG0h5+7EZlvs2d5RGReb6/+eDXNI2JzMjeSMr2Z108jmRtW+mFB+nA/AyZZXKSDenuVvs2eVaLZsBcDunntb2Wv/K6OW61DdVMuy1deB2GdLvukEEXbn5Sq/Nt6klVnoHltzxsXaNIcIwhDWPsRV59HW7tstu9mufdbLOeJ+EEcbx69E+guFlXKojtjm8DPu/D34xGmrXE11zcGKenGhFde3uSxIzfz9dc/XU+ty70QjCTK6sb/xZduIB/HjI5fKmPYxZuIVa+8EpnIplP5cjfIo1oAK5Gs/4m6xR5H+ez+Mz3lQoycqahlbrVIObOvqlei6W7kZU+1KFdep8uDiz8yZ1Kmq19yL/NP84g8rxootd9V+pwplahW9nC9IgcGcVbvX5e5TJ/XPVyr/azLH/jt5nV5qnzCeqntYXlCf8rzW9VBbW/k1+zFhT3xUECa11lz/5ztYjOvl95kvwdDKx/ntbt5Pdhm+z4oIb3SA5q9E+idVkDXyWjau07eogGQvZDnEYMkq3t+OWUcuuw95w6xwdGNMeOtAX9/7iRmzpkb/oEtT4gsr3tkA8t4ebLJOI05euLS9h5SKcAOde+OUljC3mnuhucRaRozCHrtlYhDIcBxVPRY07L3V/U4oRCdLIuK8EywLXNjMkmI47ysT3CPyCMnTSOSZCqfaRqRx8ZkEjMYZHVdj8cJSZKTTmIGw7QWoPE4ru3JIKvzr9ab53UyThgMZ3/DJM0iLl0YsbE5nkk/TuM6feXHZJyQDDLSNJr1b2vQbh8n5IOMLJuWM8+tXnenrq8iNGJ1fVZXWR6cr5lzW16Xccu1mrsRMXsd5w1Rr2zznhib6dtCDoUv0+XmrJ29CvU8Ue2aLdLWrpv5LNKUK90r3gm9E2hgRlCWYdmT2rS3iXrFvB7JMuJdxPmqnqlvi0MDDEohy4GjRy9zgQ2e+soJLoyHnNy8wK1HXwKKuPRfPv1GXnllk43NMUmcTXub5UyPMN5qZe/ZmK3LKOj9VQJR9cDrR/KyN5p5GIc1oshJ03im7iaTQjQnk5goctwhjnPG45g49lo807TYf7yVEMU56SQhTjImkyLd1tagFimLcsbjQuQvXxwSJ0UvNopyti4NiAcZly8MiQZ50QtNstqeXh4QDTOy1wYcv+k8589vECcZaSm2WRYxGKWMtwbFU0caEQ8ysrTwc+vysKiHNCIZply+OCy+Lw0xg3wSkYxSti4PiOO8Tu/AYJDV9TPeKo536ZUNBkfHTCaFfTKOGAzTWrTTtKi3PDeSJCvj7Hkt3HWog+CmyzQEE4p51eOOo/AZoju2uyiU12wbbdMjq5DHImHuemLtCkEuI5Ztabra9W7y6gu9FOidiPNB0dbzbtse0ozdVo2l6ilXNC/oyJzRIGV03Xleu7jBi68c5YWvHuOJ7DS4YUnOYJRy7PglhklWhzPCvDIvxbjyj6k4V4247rGZMwx6l1XvqhDdQhzzPGIyKR77K9E1c7IsmumZb6VFmiw13GGcG1HsZGkRgti6VPQ4ty4kWOR4Xlxy2ThisJFy+eURJDlkBplhWxEcT/GXN8g3cvLS93wrwoeOb43wUU52CYidfGuIj3Ly1wZYbng64NrPR8Rfk5FfjsknA0icycUEYmfyyhAf5GRxme9rAxg4k5eH+DAvKsmN7MViPbs8JN/IMTdw8OeG5DdO8NdGRXoDzMnTEXZ8TH4xgXKGzfCaMZNLA5KNCZOtAVGSM94aYFHOxddGDDZS3L3+sSwLbqbV+ajEufquzm34ZFTdeKvrrk2Am3Rdz9V3W9w5mnMNzyOuZiEF4a0uP/YjtHiY6KVAryuzvZTiOzLHbSqmVf+oboBB7xbg2qOXmGQRkTEzRW80SGdjzGWYpPlIF5WNuopBR7U/s72Vulcd5TPhiyJuXvifJEWvEyDPIuIkq8U5TyNIpvvm5cs2vhXjR1IwyDPDL8dMJhF2IcY3i/1xiC7ETDZihi/ETE4UAh2lMPpqxPhSxOhF49JNgIFHsPF8xNYNORtfibh8E3hUKMfGczGXT8Hg1YjxDRnJ+YhXb88YfvIkdn2GD5wTN73K+cevJz2WM/pqzOQawxPIN3KGX42ZXJcz+krM1qmyMA5Hvhyxdb2x8ZXSD8ByY/PLxiUfMHjFuHwjUNyvGL0UsZUNiSdGdgxSNyzO8csxWZzjmUFsWFzUgV9I8FFxk6wGA8PrqDpT1blsG4Cszm/zWgqvrybVgO4kixkNJvV12UXXtnlPm5Ut9NXNW8cA5uVblTd3YxhnXJokM+8TzGPbOMuashYC3ecYUZOwZ1uFNPKyV1ZvN6/FNow3VrHGyGCUpIzKs5PmUT2tKzxOjuN5RBLneBWOqAeBijmU4UBgVw02L+JpmCS40ZThGCrRT/LpoBYQJWXDjFIsroQE8o2MZJAxiRyLvKgHc/LYiTcyxgCxQ27kiXH5ZE5+XUo+Ssg3pse8fNrxQc6lgeOjyg6XXwc+zBkPHAY5k+sK37bK9filAcPfP8Hga43JTRmXE4fEi18MBMY3Fukvv85hmENe3EAunjZ8kHNxEJGXxzM3LtxsZNenpEeSwA9n60bHjo/JLiYwyLHIiRIn30yJg1i15xEW5UTHJ3ge4VFW1FMZRgqvkZnronEu62ugjEHvtOcZCmrXW5HFsdtFumvwe7chi9ljz5Yd4KULR9gYTpbO4zCIM6yBQFciF6+orudN71lm36o3HE59qn4DOqIUq468qsfaJMrZSpN6il7TpypvZ7b3FQ4MwrR37T6d9hb2tPJS7NNsOqBUxEaN4TAtBhZLsYmCgb8qBg1FDzxNY5IyxjwYFrHpI8e2pjHpSczgmi3yPGJ0aszW1mAmVh67EZ8eM9lKiJIi1hzFGdkkJr62iCnHSdFDio9npOO4iCVPYuLNSVGO41tMxgn+uozn7xrBictEkTM4tkU6SYoYdBaRXDMu/Do2JsvK6XZpRHLtVnGc61OsmiEyiYhumOBpTHxqMp0LThGDznPDRil5FhHFOVkas7E5IcsihhsT8mwagx5uFscrwklWh5WaMejqnIVXSFsMujlDJ7xeq1k5lNeSmW9r+M3QSE4okvNf1w9pm85X3XDSRoija2xnNopeMBqk9Tz4ZQT/MIgzrIFA72VkuGtUd94jWUVb+nB9Gcx85nchwhDH9I3CKIhVT3vcYaO8NElaY4uRgQflSIK5uqE4N2cFVGUJB5RmGnHcsMdFKGMQ/LNLlkV1HHs0CmdnxIE9LUWoEPPhKC38TAqRHyYTxmnMkSPj4li5cfnikCh2PDeOHN2q852ME0ZHJqSTmI3N8TZ7kW89E5vJOGGjzNePTP/kYDJOGG0U9kGZOkmK9KNR2UMbFek2NsekaVSnZ6PavxDd0TDwo5yNkmVhOfNgtkch4sXAbjEwOBikGBCXT/3VtL4wlBFOaYwb5zEMfTQHB2fOc8u1Hs6Lb9unGcNue0W/bRZH1/zqttBGVy+8LV11TXb510y/Lk/ci+i9QO+lotsG+BaJbSWC+zFYEUfFFK7RIC1jt9NtdS/amMYTrRzkc+PSeFD7kuXRzMyNZh5hOaqGWsUmI7dtPa22+FzVyEN7HOX1/jC9YUTmxEk2fWGkKq85lmQzoZQkyYo3IMu8Q/ExczaGk/rfZKIIjhwdFzNJJkkdXgEYjtLCn3K2SLVttDEJYpu+zZ5mUT1lbjhIsdH2vw+Lo3ybPTxe5a+7MdqYFINy0ey1NRqlxTzoaLZ8w2FezwOv6tDMieLtMdlqsK+K24b2NmEycwje6KzHIVrm37fFqUP/m3TZQ1GeN+UutO2Utql/824kzfK3dcrWld4L9KLR6Hn7dY1mh1PRmsdZNP2o2rasgEdWzCeu3vqLrJhxMX2cKxsV07hj8VbgdApcOHOj8jsyIBDnrouyGkhK4uLRNs2mU/GaA4ihvTlY01bmOBCTSnirXl3zTcLqTcOuaWCDOJuZyz0oB0WTslcJxY1gEM/eGGyBPYqc9MIAu8bxAXUIodpeiec8e1Wu6qWTpj2cUVH9wFWzXqreb/M8VfXUZW+KbZsIzQ5Ozxen6WDybBijshXna7refJOwa1ZHM4+9EF5rbWWv7M1yd21bZ7HutUDPe+RZ9GjTvJArusS1KxTStl+zQc0T62YDy73jAtzWayoG5rLcWl9wacYFwzK3rede3gTKfNrStfVEuuxh2eNAlPIgZFPHsaN8xl5tq/wP842C/KvUofi524zYZXm00H7dqdda7RVVaKjL3hTZbfayXM0QRCiyeXBuq/JFgTiHPdy2+m+jS6DCnwCoWNTDDe0zvXFrT7/fr3x3tae2ttKsq9124rrok6Av/LEkM3uvmT1vZo8HtuvN7MNm9rny+0Sw7R1mdtbMnjSzt+3UoebJafYQKntX+kX2tvUqfVe6ZkPoajht+VZT4bLcZi7qOJr+swpMfzEty23mE/oW/jB/1Qir5bBXUx2z6Vcohl1l3M13WCfhizPhU0jUEO3Kn9Bvs9kXb6Jy3cxJ4qzOP0xX2Zt+7dbelr8FflT2uLHezK8KTYX1UuXRVk/hWEVIV2+1Os/hp5m+7SWSrhdLQlvoz25YVtzmtcO26zcU5mp9Xvq29WWefPsizrCEQAO/AdzdsD0IfMTdbwc+Uq5jZncA9wFvKvf5FTPbPjt9DjutnGVFeZl8unrWXXluiwc29gsbX5ZHRcy4pSfTFP3mcijAuU9/GKkS54quhtfWWGtBbGnA8240i2h79OzqGTb9DkWqKX6LevahiJoVMwZeOXuCO04+N2PvSr/sk8Oy39UNplnOsB6XFdcw1DDvzb2d2rvY71e2w7bQ1U7C9eZyVzhj3pN0F8tewzuxHyQLBdrd/wp4sWG+F3ioXH4IeHtgf9jdt9z9i8BZ4K79cXUxe73z7eXmENraLsIkLn7wKOxJN0e/wwaZxPlMg66nGwY953nM6yktavjzhGKeSEQ2+8cEob1L9Jf1u6tX13WTSeKM69/4Ik++eGqbH3u9KS1TnuoGUN1cmzfHnYjpvLK3sRN72/W6H/ausZxm2mU6VftlX5aDync37DYGfZO7Pwvg7s+aWdUKbgb+Okh3rrRtw8weAB4AGJ46vks31odCuLLy0Xe25xv+ZGjEtFdQCzjT9NUFncSzszcqup4EuthLr2sn+y7Kr+l3c4AKZuPubSK7zHq4f9eA16KBsN3Yd8u8J7vmo/6i5fAb9iaOO92nT6K3Tuz3D/a3XZ6tZ8Dd3+PuZ9z9THJ8c5/d6Cdx5CRxXv9SXPMHZ8LBkGZsOfwHlKYILHrk7jtdfs4T2Z2wrMh3HfegWfTo3wydzAuHzVtuO+5efN7vPMV2divQz5nZaYDy+/nSfg64NUh3C/DM7t07nITx1apH3RwYDAcMw8HAcGCxyiv8XmRvbl+0vCj/efkua9/JwM2yfnSxU5Hfj3J3CXD1Xd2Qw/TL5i8ON7sV6A8C95fL9wMfCOz3mdnIzG4Dbgc+tjcXDyeV2FafcMpQmkX1tKwwzW4baVuPatEc2uZAzLLTwfZq77opdM313utjetv2pn9t9dP0t1k/i/Kd9y1ExcIYtJn9HvAW4KSZnQN+GvhZ4H1m9iPAPwA/CODunzaz9wGfofgzjB9z96w1YwFMG2X14/lx1F5d+9F4u3rK4XqbfZkpiHsV55CuGR/L2tvEvMve3H+Zci9jX8bv5v5CNFko0O7+wx2b3tqR/l3Au/bilOgXy4pHM92i0MV+hkna1ndqv1L+CbEsvftXb3E4cLeZgc1VzCEVYt2RQIsDIfxp1KtdnK/28ovd0+vf4hDrS2Tbf4ynyV7e/FwHqqcIWN3vmYv1RgItDgQzn/mFvqsRMy/+bEDiLHaJBFocGIexV7xTlv0PPSHakECLldFXAe/Tz02KqxsNEgoRsO13uoVYIRJoIYToKQpxCBGg0IboE2vdg9YcWyHEYWZtBToUZom0EOIwspYCLUEWQlwNKAYtrnrC/w7MdfMXPWItBVoDOWK/aP6xqxB9Yi1DHELsBxJn0Xck0OLQs2i2TxXWUHhD9I21DHEIsRO6QmKhIEucRR9RD1oIIXqKBFoIIXqKBFoIIXqKBFoIIXqKBFoIIXqKBFoIIXqKBFoIIXqKBFoIIXqKBFoIIXqKBFoIIXqKBFoIIXqKBFoIIXqKBFoIIXqKBFoIIXqKBFoIIXqKBFoIIXrKQoE2s/ea2fNm9nhg+xkze9rMHi0/9wTb3mFmZ83sSTN720E5LoQQh51letC/AdzdYv9Fd7+z/PwxgJndAdwHvKnc51fMLN4vZ4UQ4mpioUC7+18BLy6Z373Aw+6+5e5fBM4Cd+3BPyGEuGrZSwz6x83ssTIEcqK03Qw8FaQ5V9q2YWYPmNknzOwT6asX9+CGEGJZ9E/m68VuBfpXga8D7gSeBX6htLf982brFeHu73H3M+5+Jjm+uUs3hBDi8LIrgXb359w9c/cc+DWmYYxzwK1B0luAZ/bmohBiv9C/l68XuxJoMzsdrP4AUM3w+CBwn5mNzOw24HbgY3tzUQghrk6SRQnM7PeAtwAnzewc8NPAW8zsTorwxZeAfwPg7p82s/cBnwFS4MfcPTsQz4UQ4pCzUKDd/YdbzL8+J/27gHftxSkhhBB6k1AIIXqLBFoIIXqKBFoIIXqKBFoIIXqKBFoIIXqKBFoIIXqKBFoIIXqKBFoIIXqKBFoIIXqKBFoIIXqKBFoIIXqKBFoIIXqKBFoIIXqKBFoIIXqKBFoIIXqKBFoIIXqKBFoIIXqKBFoIIXqKBFoIIXqKufuqfcDMvgJcAF5YtS9LcBL5ud+si6/yc39ZFz9hsa9f6+437vdBeyHQAGb2CXc/s2o/FiE/95918VV+7i/r4ieszleFOIQQoqdIoIUQoqf0SaDfs2oHlkR+7j/r4qv83F/WxU9Yka+9iUELIYSYpU89aCGEEAESaCGE6CkrF2gzu9vMnjSzs2b24Kr9CTGzL5nZ35rZo2b2idJ2vZl92Mw+V36fWJFv7zWz583s8cDW6ZuZvaOs4yfN7G0r9vNnzOzpsl4fNbN7euDnrWb252b2hJl92sx+orT3qk7n+NnHOt0ws4+Z2adKX/9Lae9bnXb5ufo6dfeVfYAY+DzwBmAIfAq4Y5U+Nfz7EnCyYXs38GC5/CDwcyvy7TuANwOPL/INuKOs2xFwW1nn8Qr9/Bngp1rSrtLP08Cby+VjwGdLf3pVp3P87GOdGnBNuTwAPgp8aw/rtMvPldfpqnvQdwFn3f0L7j4GHgbuXbFPi7gXeKhcfgh4+yqccPe/Al5smLt8uxd42N233P2LwFmKul+Vn12s0s9n3f2T5fJrwBPAzfSsTuf42cUq69Td/Xy5Oig/Tv/qtMvPLq6Yn6sW6JuBp4L1c8y/2K40DnzIzB4xswdK203u/iwUjQU4tTLvttPlWx/r+cfN7LEyBFI94vbCTzN7PfCNFD2p3tZpw0/oYZ2aWWxmjwLPAx92917WaYefsOI6XbVAW4utT/P+vs3d3wx8L/BjZvYdq3Zol/Stnn8V+DrgTuBZ4BdK+8r9NLNrgD8AftLdX52XtMV2xXxt8bOXderumbvfCdwC3GVm/2RO8pX52uHnyut01QJ9Drg1WL8FeGZFvmzD3Z8pv58H3k/xGPOcmZ0GKL+fX52H2+jyrVf17O7PlQ0iB36N6ePhSv00swGF6P2Ou/9hae5dnbb52dc6rXD3l4G/AO6mh3VaEfrZhzpdtUB/HLjdzG4zsyFwH/DBFfsEgJkdNbNj1TLwPcDjFP7dXya7H/jAajxspcu3DwL3mdnIzG4Dbgc+tgL/gLpRVvwARb3CCv00MwN+HXjC3f9HsKlXddrlZ0/r9EYzu65cPgJ8F/B39K9OW/3sRZ0e9AjpEiOo91CMRH8eeOeq/Qn8egPFSO2ngE9XvgE3AB8BPld+X78i/36P4rFrQnFH/5F5vgHvLOv4SeB7V+znbwF/CzxGcbGf7oGf307xmPoY8Gj5uadvdTrHzz7W6dcDf1P69Djwn0t73+q0y8+V16le9RZCiJ6y6hCHEEKIDiTQQgjRUyTQQgjRUyTQQgjRUyTQQgjRUyTQQgjRUyTQQgjRUyTQQgjRUyTQQgjRUyTQQgjRUyTQQgjRUyTQQgjRUyTQQgjRUyTQQgjRUyTQQgjRUyTQQgjRU5L9yuikvc7HjKcGs+k/K25fqNPMrrel6bIFK9v+wrHNNs+P2VXv2r8lbW3Ylkdb2q48Krtt/+fJueUIjrUw773vs23/OduXyXOpsi6x3Xfjx9ztvueybStd1yVl4eL2P89oNpHy6pyT7zQP68i72aS2bQs3mzfdrNNbo7xtacJ024/T2McWbKfpi28r48L0c44Vpm3uO5tvwzeDQO069mlun1oeeWzrT939blrYN4EeM+Zbou/GoupsRMWyRVDbDIui6RmrlsvtVqW1afo6zYwtXG/ZJ8gTM9yseFZopPHaD2aO0UzvVrYCM7x65jArhGF69U33C/apxaNOT+1btd6apt5uwXaK4y/ch9l9ZrbbNtv2NNvz2JZnYF92n06/GnlUdG2fzdPnlmUmDY00FNtoHKeyNdNY2z6FZ+X2qbBbyz5mHlyCwT7V4cs0xXaI6n222wCiUqCiervPbjcv03idJiqXoyCPqG074Xreaotn8sjrPOLmujkRxXpMXucRW17nGVte5FkfI6/zKfKc7h+Xtogqj7w+TpFuuk917OoYMXm9X+1rkGdc1kfl53Q9LFu1H2W+EGPBuhFhxOXFEGGlLQrWo3o9Pv25k3SgEIcQQvQUCbQQQvQUCbQQQvQUCbQQQvQUCbQQQvQUCbQQQvQUCbQQQvQUCbQQQvQUCbQQQvQUCbQQQvQUCbQQQvQUCbQQQvQUCbQQQvQUCbQQQvQUCbQQQvQUc9/+I+G7ysjsT4DO3zU9YE4CL6zo2PvBOvu/zr6D/F818h9e6PrB/n0T6FViZp9w9zOr9mO3rLP/6+w7yP9VI//noxCHEEL0FAm0EEL0lMMi0O9ZtQN7ZJ39X2ffQf6vGvk/h0MRgxZCiMPIYelBCyHEoWOtBNrM7jSzvzazR83sE2Z2V7DtHWZ21syeNLO3BfZvMrO/Lbf9sln1x/erwcz+Xenjp83s3YF9LfwvffopM3MzOxnYeu+/mf28mf2dmT1mZu83s+uCbb33v4mZ3V36e9bMHly1P03M7FYz+3Mze6K83n+itF9vZh82s8+V3yeCfVrPwyoxs9jM/sbM/qhcv3L+u/vafIAPAd9bLt8D/EW5fAfwKWAE3AZ8HojLbR8D/hlgwP+t9l+R/98J/BkwKtdPrZP/pT+3An8K/D1wcp38B74HSMrlnwN+bp38b5QlLv18AzAs/b9j1X41fDwNvLlcPgZ8tqzrdwMPlvYHlzkPKy7HfwB+F/ijcv2K+b9WPWjAgePl8rXAM+XyvcDD7r7l7l8EzgJ3mdlp4Li7/z8vavA3gbdfYZ9DfhT4WXffAnD350v7uvgP8IvAf6Q4FxVr4b+7f8jd03L1r4FbyuW18L/BXcBZd/+Cu4+BhynK0Rvc/Vl3/2S5/BrwBHAzhZ8PlckeYlqnrefhijrdwMxuAf4l8L8D8xXzf90E+ieBnzezp4D/DryjtN8MPBWkO1fabi6Xm/ZV8Ubgn5vZR83sL83sm0v7WvhvZt8PPO3un2psWgv/G/xrih4xrKf/XT73EjN7PfCNwEeBm9z9WShEHDhVJutjmX6JokOSB7Yr5n+yl50PAjP7M+B1LZveCbwV+Pfu/gdm9kPArwPfRfH42cTn2A+MBf4nwAngW4FvBt5nZm9gffz/TxRhgm27tdh657+7f6BM804gBX6n2q0l/Ur83wF99m0GM7sG+APgJ9391Tlh/F6Vycy+D3je3R8xs7css0uLbU/+906g3f27uraZ2W8CP1Gu/h+mjx3nKGKjFbdQhD/OMX2MDe0HxgL/fxT4w/Jx+WNmllO8y997/83sn1LE1T5VNrBbgE+WA7W997/CzO4Hvg94a3keoEf+74Aun3uFmQ0oxPl33P0PS/NzZnba3Z8tw0hVqK9vZfo24PvN7B5gAzhuZr/NlfR/1QH4HQbrnwDeUi6/FXikXH4Ts8H5LzAd5Pk4RY+1GuS5Z4X+/1vgv5bLb6R4HLJ18b9Rli8xHSRcC/+Bu4HPADc27Gvhf8PnpPTzNqaDhG9atV8NH40ibv9LDfvPMzvI9u5F52HVH+AtTAcJr5j/Ky/4Divp24FHykr4KPBNwbZ3UoyaPkkw0g6cAR4vt/1PypdzVuT/EPjt0p9PAv9infxvlKUW6HXxn2LQ5ing0fLzv9bJ/5by3EMxM+LzFCGclfvU8O/bKR7xHwvq/B7gBuAjwOfK7+sXnYdVfxoCfcX815uEQgjRU9ZtFocQQlw1SKCFEKKnSKCFEKKnSKCFEKKnSKCFEKKnSKCFEKKnSKCFEKKnSKCFEKKn/H8SIMJfMKMkEgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(v_sql.reshape((180, 360)) - v_hpx.reshape((180, 360)))\n",
    "plt.colorbar(orientation='horizontal')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "settled-inspection",
   "metadata": {},
   "source": [
    "## Speeding up `gen_dist()`\n",
    "\n",
    "Approach here is to use an interpolation function instead of sqlite"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 218,
   "id": "municipal-parcel",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>z</th>\n",
       "      <th>dist</th>\n",
       "      <th>vol</th>\n",
       "      <th>dvol</th>\n",
       "      <th>cdf_sfr</th>\n",
       "      <th>cdf_smd</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.00000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>0.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.00001</td>\n",
       "      <td>0.000044</td>\n",
       "      <td>3.630882e-13</td>\n",
       "      <td>3.630882e-13</td>\n",
       "      <td>2.816478e-17</td>\n",
       "      <td>8.261876e-16</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.00002</td>\n",
       "      <td>0.000089</td>\n",
       "      <td>2.904685e-12</td>\n",
       "      <td>2.541597e-12</td>\n",
       "      <td>2.253219e-16</td>\n",
       "      <td>6.609439e-15</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.00003</td>\n",
       "      <td>0.000133</td>\n",
       "      <td>9.803245e-12</td>\n",
       "      <td>6.898559e-12</td>\n",
       "      <td>7.604724e-16</td>\n",
       "      <td>2.230666e-14</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.00004</td>\n",
       "      <td>0.000177</td>\n",
       "      <td>2.323716e-11</td>\n",
       "      <td>1.343391e-11</td>\n",
       "      <td>1.802626e-15</td>\n",
       "      <td>5.287456e-14</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>649996</th>\n",
       "      <td>6.49996</td>\n",
       "      <td>8.638431</td>\n",
       "      <td>2.699997e+03</td>\n",
       "      <td>3.621527e-03</td>\n",
       "      <td>9.999987e-01</td>\n",
       "      <td>9.999998e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>649997</th>\n",
       "      <td>6.49997</td>\n",
       "      <td>8.638435</td>\n",
       "      <td>2.700000e+03</td>\n",
       "      <td>3.621523e-03</td>\n",
       "      <td>9.999991e-01</td>\n",
       "      <td>9.999999e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>649998</th>\n",
       "      <td>6.49998</td>\n",
       "      <td>8.638439</td>\n",
       "      <td>2.700004e+03</td>\n",
       "      <td>3.621519e-03</td>\n",
       "      <td>9.999994e-01</td>\n",
       "      <td>9.999999e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>649999</th>\n",
       "      <td>6.49999</td>\n",
       "      <td>8.638443</td>\n",
       "      <td>2.700008e+03</td>\n",
       "      <td>3.621515e-03</td>\n",
       "      <td>9.999997e-01</td>\n",
       "      <td>1.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>650000</th>\n",
       "      <td>6.50000</td>\n",
       "      <td>8.638447</td>\n",
       "      <td>2.700011e+03</td>\n",
       "      <td>3.621511e-03</td>\n",
       "      <td>1.000000e+00</td>\n",
       "      <td>1.000000e+00</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>650001 rows × 6 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "              z      dist           vol          dvol       cdf_sfr  \\\n",
       "0       0.00000  0.000000  0.000000e+00  0.000000e+00  0.000000e+00   \n",
       "1       0.00001  0.000044  3.630882e-13  3.630882e-13  2.816478e-17   \n",
       "2       0.00002  0.000089  2.904685e-12  2.541597e-12  2.253219e-16   \n",
       "3       0.00003  0.000133  9.803245e-12  6.898559e-12  7.604724e-16   \n",
       "4       0.00004  0.000177  2.323716e-11  1.343391e-11  1.802626e-15   \n",
       "...         ...       ...           ...           ...           ...   \n",
       "649996  6.49996  8.638431  2.699997e+03  3.621527e-03  9.999987e-01   \n",
       "649997  6.49997  8.638435  2.700000e+03  3.621523e-03  9.999991e-01   \n",
       "649998  6.49998  8.638439  2.700004e+03  3.621519e-03  9.999994e-01   \n",
       "649999  6.49999  8.638443  2.700008e+03  3.621515e-03  9.999997e-01   \n",
       "650000  6.50000  8.638447  2.700011e+03  3.621511e-03  1.000000e+00   \n",
       "\n",
       "             cdf_smd  \n",
       "0       0.000000e+00  \n",
       "1       8.261876e-16  \n",
       "2       6.609439e-15  \n",
       "3       2.230666e-14  \n",
       "4       5.287456e-14  \n",
       "...              ...  \n",
       "649996  9.999998e-01  \n",
       "649997  9.999999e-01  \n",
       "649998  9.999999e-01  \n",
       "649999  1.000000e+00  \n",
       "650000  1.000000e+00  \n",
       "\n",
       "[650001 rows x 6 columns]"
      ]
     },
     "execution_count": 218,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import sqlite3, pandas as pd\n",
    "\n",
    "db = sqlite3.connect('../data/models/universe/h0-67d74-wm-0d3089-wv-0d6911.db')\n",
    "df = pd.read_sql_query(\"SELECT * FROM DISTANCES\", db)\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "accurate-petroleum",
   "metadata": {},
   "source": [
    "### Setup sqlite lookup method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 197,
   "id": "advanced-chuck",
   "metadata": {},
   "outputs": [],
   "source": [
    "dt = pc.DistanceTable(H_0=67.74, W_m=0.3089, W_v=0.6911).lookup\n",
    "m = dt(z=np.array([1.0]))\n",
    "dist_co_max = m[1]\n",
    "vol_co_max = m[2]\n",
    "cdf_sfr_max = m[-2]\n",
    "cdf_smd_max = m[-1]\n",
    "\n",
    "def gen_dist_sql(vol_co):\n",
    "    z, dist_co, = dt(vol_co=vol_co)[:2]\n",
    "    return z, dist_co"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "stretch-beatles",
   "metadata": {},
   "source": [
    "### Setup interpolation method\n",
    "\n",
    "(sidenote: I think when using the interpolation method we probably don't need 650k points to be precalculated)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 201,
   "id": "pressed-baseline",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.interpolate import interp1d\n",
    "\n",
    "vol_co_to_z    = interp1d(df['vol'], df['z'])\n",
    "vol_co_to_dist = interp1d(df['vol'], df['dist'])\n",
    "\n",
    "def gen_dist_interp(vol_co):\n",
    "    z = vol_co_to_z(vol_co)\n",
    "    dist = vol_co_to_dist(vol_co)\n",
    "    return z, dist"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adult-testimony",
   "metadata": {},
   "source": [
    "### Compare timings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "announced-nancy",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_gen = 10000\n",
    "vol_co = vol_co_max * np.random.random(n_gen)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 202,
   "id": "based-purpose",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.41 ms ± 23.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit gen_dist_interp(vol_co)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 203,
   "id": "superb-german",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10.2 s ± 49.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit gen_dist_sql(vol_co)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "regulated-consumer",
   "metadata": {},
   "source": [
    "### Check results are the same"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 207,
   "id": "growing-scott",
   "metadata": {},
   "outputs": [],
   "source": [
    "z_interp, dist_co_interp = gen_dist_interp(vol_co)\n",
    "z_sql, dist_co_sql       = gen_dist_sql(vol_co)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 217,
   "id": "labeled-supplier",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "print(np.allclose(z_interp, z_sql, atol=1e-5))\n",
    "print(np.allclose(dist_co_interp, dist_co_sql, atol=1e-4))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}