2 - Speedup Jit.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f5649ca4-c075-48f6-9e20-bbbdd6877d9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "\n",
    "import h5py\n",
    "import numba\n",
    "import numpy as np\n",
    "\n",
    "from utils.measure_and_compare import measure_and_compare, reset_measure_and_compare\n",
    "\n",
    "reset_measure_and_compare()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "10c7081e-0e8e-4a2e-a71a-d124772343cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_data_hdf5(path):\n",
    "    segments = []\n",
    "    counts = []\n",
    "    centers = []\n",
    "    covs = []\n",
    "    for hdf5_path in sorted(path.glob(\"*.hdf5\"), key=lambda p: int(p.stem)):\n",
    "        with h5py.File(hdf5_path, \"r\") as f:\n",
    "            segments.append(f[\"segments\"][:])\n",
    "            counts.append(f[\"counts\"][:])\n",
    "            centers.append(f[\"centers\"][:])\n",
    "            covs.append(f[\"covs\"][:])\n",
    "    return segments, counts, centers, covs\n",
    "\n",
    "\n",
    "def combine_stats_per_segment(covs, means, sizes):\n",
    "    reshaped_sizes = sizes.reshape((-1, 1, 1))\n",
    "    sizes_sum = sizes.sum()\n",
    "\n",
    "    # bias correction\n",
    "    corrected_covs = covs * ((reshaped_sizes - 1) / reshaped_sizes).astype(covs.dtype)\n",
    "\n",
    "    e = corrected_covs + means.reshape(\n",
    "        (means.shape[0], means.shape[1], 1)\n",
    "    ) * means.reshape((means.shape[0], 1, means.shape[1]))\n",
    "    new_e = (e * reshaped_sizes).sum(axis=(0)) / sizes_sum\n",
    "    new_mean = (means * sizes.reshape((-1, 1))).sum(0) / sizes_sum\n",
    "    new_cov = new_e - np.outer(new_mean, new_mean)\n",
    "\n",
    "    # bias correction\n",
    "    new_cov = new_cov * sizes_sum / (sizes_sum - 1)\n",
    "\n",
    "    return sizes_sum, new_mean, new_cov\n",
    "\n",
    "\n",
    "def combine_stats(chunked_stats):\n",
    "    concated_stats = tuple(\n",
    "        map(lambda x: np.concatenate([i for i in x if len(i) > 0]), chunked_stats)\n",
    "    )\n",
    "    sort_idx = np.argsort(concated_stats[0])\n",
    "    all_segments = concated_stats[0][sort_idx]\n",
    "    all_counts = concated_stats[1][sort_idx]\n",
    "    all_centers = concated_stats[2][sort_idx]\n",
    "    all_covs = concated_stats[3][sort_idx]\n",
    "    g_segments, g_idx = np.unique(all_segments, return_index=True)\n",
    "    g_counts = np.split(all_counts, g_idx[1:])\n",
    "    g_centers = np.split(all_centers, g_idx[1:])\n",
    "    g_covs = np.split(all_covs, g_idx[1:])\n",
    "\n",
    "    def empty_array(other):\n",
    "        return np.ndarray((len(other),) + other[0].shape[1:], dtype=other[0].dtype)\n",
    "\n",
    "    segments = empty_array(g_segments)\n",
    "    counts = empty_array(g_counts)\n",
    "    centers = empty_array(g_centers)\n",
    "    covs = empty_array(g_covs)\n",
    "    for i, (segment, count, center, cov) in enumerate(\n",
    "        zip(g_segments, g_counts, g_centers, g_covs)\n",
    "    ):\n",
    "        count_i, center_i, cov_i = combine_stats_per_segment(cov, center, count)\n",
    "        segments[i] = segment\n",
    "        counts[i] = count_i\n",
    "        centers[i] = center_i\n",
    "        covs[i] = cov_i\n",
    "    return np.array(segments), np.array(counts), np.array(centers), np.array(covs)\n",
    "\n",
    "\n",
    "@measure_and_compare()\n",
    "def load_and_combine_stats(path):\n",
    "    chunked_stats = load_data_hdf5(path)\n",
    "    return combine_stats(chunked_stats)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "2fb150e3-95d8-489b-880b-70d857f3b185",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARQAAAD8CAYAAAC2EFsiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA++klEQVR4nO2dd3wcxdnHv7N3OvUuWZYs94ILLoDBmBZjm2bTe+cloSeU0AIkoeZNXiDUEAimhd4JLVSDwRiwjXsvki0X9V7uJF3Zef/YlSwXWSfd3u2etN/PZz+7N7u389zq7qeZZ555RkgpsbGxsTECxWwDbGxseg+2oNjY2BiGLSg2NjaGYQuKjY2NYdiCYmNjYxi2oNjY2BhGxAVFCHGiEGKjEKJACHFHpOu3sbEJHyKScShCCAewCTgO2An8AlwgpVwXMSNsbGzCRqRbKIcBBVLKLVJKL/AWcFqEbbCxsQkTzgjXNwDY0eH1TmBKxwuEEFcBV+kvD4mQXTY2fZkqKWW2ETeKtKB0iZRyDjAHQAhhzwuwsQk/24y6UaS7PMXAwA6v8/UyGxubXkCkWyi/ACOFEEPRhOR84MII22ATBgSCZFca8XFZyHGjkTGC1PE1TExp2OvaFfWpNKxJx1HrRd26kcbmMloCLUjsBmm0E1FBkVL6hRC/A74EHMCLUsq1kbTBxhgcwkFO1ngCM3O5emgV56b9L9XHTyV9kI9BSbHEKA5cwolTcez1Xr8awCv9uP1eSpr8NK52kL5oAW823MOzy/NwLFhLVVMRARkw4ZPZhEJEh427i+1DsRbxjnjSRh3NOWdUcvnwt4k7K4thyUm4lBjD6nD7W9hQ3YLr3S38peRyfngvi4qC721xCS9LpZSTjbiRLSg2+0UgyEkfzZgr4/jb4V8y8CQXeXGpEau/sKkOz7+3cM6m66l5sZhKt2H+Q5td2IJiE14EgtzBR/KvCwaTe+NTHJyTgiLMm6nhVX0sLmzh55dO4vHXoGTHj6bZ0guxBcUmPHQUkiG3/4Px6elmm7QXq2pq2fTwldz4722UlC0x25zegC0oNsbTL2kot16fx4m3fWJJIdmTn7c38Oyzx/DxIyXUtlaabU40Y5igIKW07AZIewv/5hAO2e/cWfL7xY0yoAZkNOEL+OWi98plvyNOkQJh+rOM0m2JNOo3a9SNwrFZ4EH3+i0nebj87x2XytLmesN+5GawqaFG3nfXFJkem236M43CzRYUewt9y5k0S85f0BB1rZLO8AX8cvlzG2TekGNMf7ZRthkmKHaCpT6IUzgZcvU0vpz7GkcfmWzq6I2ROBUHk644gDf/+1+yzzjZbHP6JL3jm2QTNE7h5Io7D2TJk58wMdP6jteecMzYJL554RWGX3c0DrF3pK5N+LAFpQ8R74jn6jsP5K/3LiDTlWS2OfumsNCQ24xPT2fRY59x8Z0H2qISSYzqO4Vjw/y+Za/ZnMIpv77zKtnkazbMZ2E4y5drbr2AcT6dqtZGefFdE6RDOEz/G1h4s52y9hb85hROec1dk2Sd123YDzUsTJigfyWNxRaVLjdbUOwtuE0g5JV3TZQ13ibDfqBhA6S85JKw3HqHp07mXHSi6X8Pi262oNhbcFv/i0+Q29x1hv0ww8aGDdrXsaoqbFWsqqmR2dNONf1vYsHNHja26Zrc0cfxwf++w6CEVLNN6Zpf/1rbZ2aGrYrx6em8/5fXyEkZGbY6+jq2oPRS0mOz+fCBV5g6KMVsU4Ljp5/gtNPCXs2RRyTyxU3XE6fEhb2uvogtKL0Qh3Dw+M2HMunM8P23N5Rteo6Tp58Oe1WKUBh191UcedvosNfVF7EFpReSMXM2M+9+w9BMamHlmmu0fV5eRKpLcMTy4O/nkTf82IjU15ewBaWXkRyTyn+v+FdEs6qFzBdfwPTpEa3ykJw0vvr1DbgUV0Tr7e3YgtLLOPKWIYw7K8NsM4KnvFzbP/dcxKvOu+04si87JuL19mZsQelFZCcO5q7LfiDBEWu2KcFzww3aftiwiFedHpPIB7NfJ8mZHPG6eyu2oPQijr0hlamjEsw2o3u88w4ceqhp1Y87LYXUSw8zrf7ehi0ovYTsxMH89uIf97kOjmWprdX2L71kmgmJzjg+mPWG3UoxCFtQegkzb0jliNHxZpvRPW69VduPG2eqGXYrxThsQekFpLjSuebsKGudALz4IhxwgNlWkOiM4/ODHsIpIr0yb+/DFpReQPw5U5g8MUpiTtpwu7X9q6+aa4dO6q+Hkz3yV2abEfXYgtIL+HTo0dE1sgPwpz9pexMdsh0ZlJDKr8+uN9uMqMcWlCgnO3Ew3it+Z7YZ3efxxyMWGRssx588j+SYKAoItCC2oEQ5yuljODh/j4luDQ3mGBMsra3a/vXXzbVjD8YcLEkYd6TZZkQ1tqBEOe8PHk+co0P4eEMDpKaCEJCfD99/b55xnfHXv2r7adNMNWNPsmOTOXVWsdlmRDW2oEQxCc5E1Fl/3r0wJQUaG+HGG6G4WPvRCqGVf/qpKXbuxf33Q7I14z6uz33RTmodAragRDFpAw5i4KTA3ieSkjQfhZTQ3Kw5QBsb4ZRTNHFRFHjrLe18pPH7tf0bb0S+7iDwnDGc9ARr+XaiCVtQohgxPY78+C7+08fFwQMPaOLh88FDD2nHF1ygCYsQMGcOqGpkjH7sMW1/sjUX4hqWpeAYOd5sM6IWW1CimN/lursXzOZ0wm23aYISCOxKaHT11eBwaOLyyCPauXBx++3hu7cBZLoS6TetxmwzohZbUKKYU5Lv7PmbFQWuvVYTF1WF117Tym+9VRMeIeDee8HrNcRWYFcryKLdHdAyuv0uPtFsM6KWHguKEGKgEGKeEGKdEGKtEOJGvTxDCPG1EGKzvk/Xy4UQ4kkhRIEQYpUQ4mCjPkRfxCmcuEcYFBQmBFx0kSYuUsJHH0FiItx3H8TGaudvvhk8ntDqact5cv75odscRo7IvtdsE6KXnqbLB3KBg/XjZGATMBZ4CLhDL78DeFA/ngV8DgjgcGBREHWYvbyAZbf02Gz58/YILI8xd66U/fu3SY22XXmllHU9qLvt/RZn2YsFfW1RMPOX0ZBSlkopl+nHjcB6YABwGvCyftnLwOn68WnAK/rfbCGQJoTI7Wn9fR2HMw5njAx/RTNmQGmpJgULF8KoUVpLIy1Na7mcfz5UVHR9n7YRpWefDau5RuD5VTaJMVGyWoDFMMSHIoQYAhwELAJypJSl+qkyIEc/HgDs6PC2nXrZnve6SgixRAixxAjbeivOIWMZmhHh2bFTpsDGjZo4rFoFhxwCb78NOTmauMyeDTt27Pu9b76p7a+4InL29hRhtgHRS8iCIoRIAt4HbpJS7hbzLaVsa1IFjZRyjpRyspRycqi29WoUgVOY6FMfPx6WLNHEZfNmLYDus89g0CBNXI45Ritv46KLtL1ijwP0ZkL66wohYtDE5HUp5Qd6cXlbV0bft7WHi4GBHd6er5fZRDsjRsC8eZq47NwJp54KP/ygdY+EgEmTtOsefthUM4MlOTVAQmJ/s82ISkIZ5RHAC8B6KeWjHU59DFymH18GfNSh/FJ9tOdwoL5D18imtzBggDZKJKXmW7nwQli5Ujt3220wfDj8/LO5NnZBVoKCkpptthnRSU+9ucBRaN2ZVcAKfZsFZALfAJuBuUCGfr0A/gkUAquByUHUYbb327Jb//Enyjqv26BxjTCTlKSN7lxzze6jRdnZUn79tdnW7cWPRfUy1ZVp+t84gpthozyG3CRcmwUetGW3qBIUkPLuu3e9drulvP323cUlLk7KDz4wz8YOLNhSL1Nc6ab/jSO4mT9sbGMuSnkFpZ5Ws83omnnztP1dd+0qS0iABx/UpKS1VZt93NICZ56p+VyEgJdfpn2o2SZqsAUlSmmu30F9UxT8+S6+WNvHdpKi0uWCP/9ZEw+/X5slDfA//7Nr8uI//xm5yYuA3Aa+QEvE6utNRME30mZf+FUvreVRkLejpAR+//vgrnU4tDwuUmoC8vzzWvnvfrdr8uLf/qbNmg4jqas30RxoDmsdvRVbUKKUJl8Dqb9sM9uM/bN4sbZ/4IHuv1cI+M1vdonLO+9okxbvuktr1QihHbcY35KoK33TXkS9h9iCEqVIJFuK7jfbjP1zySXaPjHE2btCwDnnaC0TKbUAuowMrbUSH6+dv/56LYmUAciL72Ppfa/y4p0nMP6mw+k/7gRSXOmG3Lu3I6SFHV9CCOsaZwHG3HAEax7/AcXMiNn9IYQWat82yzgc/PCD5qfZvn1X2aWXaomcMjJ6dMsftzQQqN71TBW3StqKNfxQ8SfuWZuE/GYF1e6dSHrN13OpNCgy3RaUKMUhHGRMn83Gz98gPcaC+TvWrNHC82trtYmEkWDZMrjsMq3uNs48E556CnKDm4fqCbSyZJ0POulJCa8keWMZ/ym8iGc/y6Z65Vf4pd8A403FFpS+iEDQL2UE6rTRnHXgDi7jZAb9+Vby4lLNNm1vDj1011wfM1i/Hn79a22GdBszZ2rpLocO7fRt5a2NbFwtCKbx4aj14/j+EU5duwL3f+fh8bsNMNwUbEHpSyQ4E0kZczSHzS7lFjGNlJjhAAjhoPk3F3H4QAsKihBw7rnabGSzKSrS0lx+9dWussMOg5degrFjd7t0dW0ttQXdW9ZVtEiS5/6b41d+Rf0nX+FVDcxyFxlsQekLJDgTSTxiBjcfVclxymk4lIS9rmmZfCGHn5JpgnX7YcsWbc5OWZmW2sBKlJVpDtz33tNGjbxeTfwAVar8XNREoKpnPinFrVL75XVc/VElZeu/jiYfiy0ovZmOQnK84ywU0fkQZkzyOAbdNJVkZ3wELeyCGTPg22+jLtK13udh5VoVQgxzidvu4c+/nMCSJ0qpaSk3xrjwYpigRDhDj83+EAj6j5rOlWc3cJrjmP0KSRuBpgIK6g7moCwLCcq338JJJ5ltRbcp9bSCr3vdnX3RMiiBP+b/QLLreWa/9zml67/q+k29BLuFYhEy4nIYcM1Q5iTMJM7Zr1vv9Rx+MUecZJE4ieJibQnU7dth4MCur7cIoXZ3OiN+ayPnzbuA7c8uICDDuDxJaBjWQrFoAEPfIm/Yr/jL78fySsqF3RYTgNSVC6hsNSaoK2Suu07bR5GYANT6PARqjf85NA9N5pXzPub3d04kOcaCznODsQXFRBzCQfZZs3nzvElMjTmjx/fxNRexfoNFEqF+/DEcdZTZVnSbbXV+CFMDQk1UOOXk73nq1iPJSR4enkosgi0oJpHkTObgmw/iszFHtA8Dh0L6/A24/SbPkK2q0vYvvGCuHd3E7W+hqTLM7kQnDDnjbb68+Wb6p48Jb10mYguKCWQnDua424bybMJFOJUu1iYOErV6McuLwjsLt0vaZhWPGmWuHd2ksKEZQlzDLFjqZ1/KFzf8rteKii0oESYneTg335jH3TFXIIRx6QekDJD6+UZzWymvvQYTJphX/75oatpvLhW3v4W60tBHdrpDm6hkJw6OaL2RwBaUCNIvaSg3Xd+PE50XhOX+avVilhWY1Epp0FdQeeUVc+rvjEsu2a9PZ1NdS8RaJx2pP+lS/nLDiF7nqLUFJUIkOZO56Ib0sIkJaK2UjI9/oqyloeuLjeZOfeH2iRMjX/f+2E+W/Vqfm8Zik0KxFBh9yoecc9sQHAa2VM3GFpQIkOBM5Khbh3GJ45Kw1+Vzb2LLdyqqjFzKRACefhqGDYtsncHwpz9p+8LC3Yr9aoANpQEwc9qNEy4/bj79/+dYE40wFltQwoxDOMj+9WH8NeZ/DPWZ7I/kxR+yvCKCrZRmPV3iq69Grs5gaUudcPvtuxUXNDbgqzD/668mKnw48yXyhv3KbFMMwfwn2svJOHo272YfGzExAQgEGol9Y82urs/69ZrDNFzce6+2P+KI8NURCmPGwAcftL+s9jZRsT0mqBQFkcAzMo3/XDCNVJfFJnn2AFtQwkj/9DG8cEQyMY60iNftr1tO0ZdevKpPm1dzSRi7Ww89BFlZ4bt/qDz5pLYPBGgJeNmwU+00gZJZtMy+lSG/HWG2GSFjC0qYSHAmcsq1Lga6pphmQ/zyd1i4ohW2bdNmAIeDtgz0r78envsbwcyZAKgvvcjaKs9u6R0tgxOen/weeYOjL8q4IxZ8sr2DxNOmcU0EnLD7Q8oA6Z++p7246abwVPLQQ9r++OPDc38jufEm3DusO8HeMzKNey6MxSmsa2NX2IISBnKSh/PCyLSg0g+EG9fOdQCsnjo1PBX86U+dL+JlIWp/ey2Kx2MZv0lnjJv+PulHRl/qhzZsQTEYgWDY1WmmdnU6kr3oOwCU5xexurbW2JsH9Nl0b75p7H0NpqylgfXnaXEyccVF5hrTBYFUB+8dPSxq1wWyBcVgcnMP5WHXdLPNaCdj1S+oigOfpxBlzs+sqjFQVNqyy5/R85nS4WZncz2bNiv447R8MUP/8WeTLeoaZeb9pB5xgtlm9AhbUAxEIDjgEh+JTmvlAqmcMg0An6cQx3M/s7LaIFGZOBE8JsStB8mmhlq2bHS0j+h4Bo0g+9sPTbUpGNRkhaeOciKwSEqKbmALioHk5h7KA86ZZpvRjvBrIzBVhx3dXubzFOJ45nMWbG40Jpo23kKpJ3VUqbK2rpaygpjd8sMW3vKwdhCwbOa0dgYdPofcwUeabUa3sQXFIASCERZrnaSvXgKAL3X3FfT8vmpS33yfHxe4zc+hYjCeQCtLyxqoLojZK2FS7VRN7HM+e8MEy7qHNy+OKReaMCcrRGxBMYh+ycP4s+NQs83YjexF33d6LhBoJOWbV9n0ajFb3XWRMyqMlLc2srTAS/NO535Hc4Y/cnvnJy3EAxmPRd1sZFtQDEKdPoZMp7WS5sRXlODOG9TpeSkDiKLP8DzxHfPXNdESiLoFqgDwqj5W1dSycb1ANuzf71B83rU4PU0Rsiw0aidPJmFidAW6hSwoQgiHEGK5EOJT/fVQIcQiIUSBEOJtIbRgDCFErP66QD8/JNS6rYJLcXHtuPKIztcJljaH7P7wte4g9Z1XWPdaMevqayM/U7mHqFKlpKWeX4qaqdsSE9R6OtuuuAOA2NLtXVxpPmqiwvTjdphtRrcwooVyI7C+w+sHgceklCOAWuA3evlvgFq9/DH9ul5BVr8JnIS1hvlc1RUA1I+dFNT1UgZQtnyK+o/PWfhFvTk5VbpBtbeJX0obKVjn0ELpgwxY8+v+pKH/vCeM1hnHHa7fEqfEmW1G0IQkKEKIfGA28Lz+WgDTAT3em5eB0/Xj0/TX6Odn6NdHPb4jc3q0/EU4yf5lPgDS0b0wbr+vmoSFr1H9+Pf88GMjW911lmmxqFLVhaSOteugtdjRo0z1zfnD6PfVe11faAE8h5xFRtbYri+0CKG2UB4HbgfavnGZQJ2U0q+/3gkM0I8HADsA9PP1+vW7IYS4SgixRAixJETbIoJAcPoBO802Yy+yF88P6f2+5iKSv/o37sc+Y/GntSyrqMMTaDXIuu7hVX3sbK5n0Y4m1q5Dc7qGkOmy4BZ9/tF+cs1aBW+/WPzTLbY+9H7o8SwkIcTJQIWUcqkQYppRBkkp5wBz9DosPvMC0mKzOF21WNpDnaqDQ89P4vdVE7fkDZSlsWxJn0jt5LFkHuhnSFICCY7wzeHxqj6qvB7KGiRNVU7wOAybh1N7hDaRsd8Xb1MxK3wpOQ1BgZOGlbQ37a1OKNMajwROFULMAuKAFOAJIE0I4dRbIflAsX59MTAQ2CmEcAKpQHUI9VuC+Nwx5Dgt1iRVtX5A1WHGZQFTZStqzWKSv1qMf24CW1IOoGnYaPxjY8nNDZAV6yLREYtT6b5jWpUqXtVPo7+VOq+fmkaF1jqHoSKyG3pPe/gjt1tfUICL5RReYRXS6jMbCUFQpJR3AncC6C2UW6WUFwkh3gXOBt4CLgM+0t/ysf76Z/38t9LKCysHiW9iCg4lwWwzdiN142oAWvrlhuX+qupBrVtO3LLlsAw8jmR2ujLxDJsMR8fjjIEU3Y+oCEiJcaIIgSolDT4/qv5Xr/doQautbgdqswCfgEBklrQoPucqBrw7JyJ1hUr/fpcR73wdj99ttildEo7EC38A3hJC/AVYDrQtI/cC8KoQogCoAc4PQ90RZ8TIUrNN2Ivshd9FtL5AoJFAcyOu+hHU79QWLtttpWUHtDc1AtbI9bHtyjsZ8O4cYst20NrfOtHN+8I9eCSpKUPw1Kw125QuMeSvK6X8DvhOP94CHLaPa1qAc4yozyo4hIMjXdZrZCXt2EJLZuRHnTodUbLg1Bl/mpaycug/72XDA9ZeOtWfGIOITTTbjKCwI2VDwKW4mKTmmW3GPgkmoM1ovKnRFSbenDeEfl++Y7YZXaLGO/APslZYQmfYghICiXFZDLBYwK+zsR6A2gmRn1cklegKKyqMkuFjGSdgoPVawvvCFpQQcLqSiBNpZpuxG1m//ACA6rJ+WkazqTnqRAD6ffWuyZb0HmxB6WW0pXy0CYK24eO/32bYLR2NdQx7/E6OOSyZYw5LxtFk7SkMRmMNl3uUoqRk4lSslWDI4fNSO+5gs82IGkrOuoK8958P6R7pC79h+KO3k1C0qb0sEJfA5j88RiApJVQTowpbUELAn5eOU0k224xd6GE9lVN6x7KWkWDbVXeR9/7zuMqL8eYM6PoNgLO+hsHP/x8D3n5mt/LK6aez9fr7aRkwNBymRgW2oPQikrZq/yE9+X33C91dfOnZAAx9+l423vdcp9el//QVIx65nfgduxZd9yckUXjLQ5TPuhAc1ktdYQa2oISAc2spvkCdKUuN7gvbf9IzWnIHkfP5W7sJirOumsHP/Y0B7z6727UVM89k6+/upzVvcOQM9IPqiQ53py0oIaC2uglIL5EJFu+a1M1r8SWa1wVz1TXQmm8tn1IwFNzyEAfeej4ZP3zO8Ef/QHzx1vZz/sQUvRVyASjm/KgdjX5cKypNqbu72IJiUfot+Jq8bz8BoPqgw6k8bBotOV0H0ZkR0NaG4o+uFJLOuiqGzPkree9pLZMDbzkXgIrjz2brb++jNbfz9JmRxNEaQHqbzTYjKGxBCQFfSx0NailxGB/FWDNpCjHuBrIXfU/m8oVkLl/Yfq41PYvKKdOomTSlPd7E0aytj1NtQMqCnuJoqkabUG5RpCRjweeM+PvtxJVuay/2paQT01CLNy2ThV8VmWdfJyQWbaWxbovZZgSFLSgh4PHWUUl5GOQE/EkpFJ9wFsUnnKUVqCqpm9aQveg7krYVkP/Fe+R/sSvrmKef1noJxJs381k0V4EfS32rYmoqGfLsA+T+56XdystPOp+i6+6hNScfgJF/u2Gva6xCbelruP3RkVjbQn/66KMl0MICZSfjIlGZolA/egL1oye0Fzkb68lasoDsRd+RUFECwKQHbgTAl5BE5ZRpVB9yBIGEpEhYiOqvQ/gk0mliCL6UZH73CSMeuZ3YiuL2Ym9aJoW3PEzl8We3B7R1pOjqP5P7n5dwVZbizQ5P2oee8r6YHxW5UMAWlJCQSH5ujuVqk/yg/uRUyo6dTdmxs5l0/w00jBhDxdQZZC/6jtRNa8ib9yl58z5tv75h2AFUTplG44gxIIx3MKqqh5jaVrzxkU2qHFNdwZB/3U/uR7vnNSufdQFbr7sXb7+ufU++DG34eMgz97Hp7n+Fxc4eocJ7xVlmWxE0tqCESPG6HORhAUssoVE5ZRpNQ0fRNHRUe5mjxUPG8oVkL/qelC0bSdmysf1c2dEnUHbsbMPqlzKAq64Rb16YBUVKsuZ9zPBHbiO2clc+Gm9GttYKmXnmPlshXdHabwD9P33dUoLirPHBN9HR3QFbUELGtbYS/+RGU2NREnYWAdA4bPRe5wJxCVROnU7l1OntZfEl28le/D3+MAwxx1Ztp4lsw+/bkQnXziJt2QIAyk6+iKJr7zGkm1Jwy0OM+8NFWsSxRRZkSF29mLqyVWabETS2oIRIY/k6tgXyGOEwL9w9e7G+5GiQP4LmvEFsP/2SsNiiureiuA9CTQxfzMbqJ/6DdMUa/qOvnnYKAFnffEjVzDMMvXdPebP+T1GR+rGN6Ai/szBNvgbeEptNtSF9zVJUpzXC6wKBRlzV4V2AXcbGhacF0Z682rjZx6EgWiRzFlgzgVdn2IISIhLJ9yuGmG2GqQFte5Kwc2vXF1mU0tMuI7a63GwzAEhbsRHvwh/NNqNb2IJiAM6Fm2n0FZlSt/Bp0alVh1pnUW1Zvx7Fbe0saJ1RdM3dALiqyky2BJ4tu5J6b3StNGMLigFU1m3i3+IbU+rOWPULoEV7WoVAoJGEHbVmm9EjfHpy7yHP3G+qHc5qPx98EF7ndjiwBcUAAjLAJwsHI2Xk07tbdYZxbMlaLWo2CmnNzqX/J6+aakPMgn9Qvel7U23oCbagGETgp0UU+X6OeL1xVeW484dEvN6u8DcXEVccPfETHSm85WHtwKR16ESzZPbqX/Cq0TXZEmxBMYyalnLuaCru+sIwYCWHbBtSBkgsXBaVrZSqY08FIGvex6bUnzT3BdyffWtK3aFiC4qBNLxaRZVvTcTqi63SRiPqxlhzsfaobaWYOHwsWiQnrPw2qmJPOmILioGU1W/kd02RE5T2gLYeLFAeCaK5lVJ2yiW7hfVHCuWLB2j679yI12sUtqAYTN3LlRFrpWQtWRCRekLB31xEUkF0DX0CFF2rDR/HVEUuJsVV0sLZn6ymORAdyZT2hS0oBlPesJkrqjdHbMSn8tCjI1JPT5EygGvbjzjqLbjA8X7wZvUHYMizD0SmQhVeXnwc5au/jEx9YcIWlDDQ8PIalvo/C28lAe0HWnWY9ZfM8PuqSV29GqIs1q01M2evlAjhIvP7H3nvieqoyXvSGbaghIHa1kru+iEOv9oYtjrS1q8AoDUzOhbRDtSvIqGwzmwzukXhzfrax2EePnaVtDD9g6eocu8Iaz2RwBaUMFH74+f80Re+NXOtGtDWGVIGiC2Yi7PKZ7YpQdM24zjzu0/CVofwSq6bP5vKRf8NWx2RxBaUMOGXfpY9Vc0SX3i+jInF22i2WKrCrggEGkld9TOiOUqa9frw8YhHbg/P/VUo/fQyNjyxOuq7Om3YghJGqt07+eO3sfgCdWG5fzQuOerzFJL23ftRM3mwbPaFu+WmNZKYj/+X3z65OKpHdfbEFpQwkORMJnfUDHIuPp7RvyrmJ3WuoaM+MfXaxLva8ZMNu2ekqPNt4qRtr/HR18dGhagUXXcvoGXPN5LM+Us5+6VvqfKYE10dLuyMbQYhEPRLHkbgiHEMPbCEfFdlu1r/RB3p/q8Z7zzOkNyzWb/8AICMcYV8r0jiC9RxTdMyvJsyebkgETiW046bF9bsbqHSllpy8Jz/peCOxw25Z+b8pRz/xB8pKV1syP2sREh/SSFEmhDiPSHEBiHEeiHEVCFEhhDiayHEZn2frl8rhBBPCiEKhBCrhBAHG/MRzEUgyEkeTtZZszjwumR+dXARg1ze3R5sAPiCnez0LTGkzmhzyIImJlc3f0rNSj2Duyp5eUE+H0RBS8WblkneBy8Ycq/M7xdrYrIjuhInBUuo/xqeAL6QUo4GJgLrgTuAb6SUI4Fv9NcAJwEj9e0q4JkQ6zadrMSBZJ0xmwOvS+aYUTtIUzr/YXiBd5TVlHhDFxUl4Kcmiro7vkAd1zZ/ws4l/UHt4HxUJa8tyOejr4+1dOCbIbOPVUj97ysc/+TdvVZMIARBEUKkAscALwBIKb1SyjrgNKAtGuhl4HT9+DTgFamxEEgTQkTXMIVOnBJHv6mnMPLafhwzevt+haQjzRLeVFaEJipSq8uKM4z3RZuYbF+8h5i0obdU7vhqOnHbPZE3MAgqjz8bgIwFn/fo/cIrafnwJmY+/NdeLSYQWgtlKFAJvCSEWC6EeF4IkQjkSCnbZlWVATn68QCgY+TOTr1sN4QQVwkhlgghjOkfGEz/1APIuWIqR07bRl5M92Mq2kRlvbdnjtrkwg3affKssZD3/mjxV2hisqSL/xuqZOPPAzh/7rlkzl9qvYjatuHjv3d/+NhR6+eFz47mrL+/0+scsPsiFEFxAgcDz0gpDwLc7OreACCllNC9AXYp5Rwp5WQppaXa9E7hJGfSLMZcncDkzOqQHlyzhA9FESv8XxBQu/dfOVr8JxXelZxX+VPnLZN90LQmhWPnP4Lnw+st51cpP+n83RZYD4b0Xwo4/dHLeO1va2jyhy9q2kqE8rvYCeyUUi7SX7+HJjDlbV0ZfV+hny8GBnZ4f75eZnmSnMlknHUCU04sJsNhTF8/APyXUj5W36fFX9Hl9QCumkpSCjfgTUkzxIZwIGWAJb5PuHhzGc3rE7t/gxq4dGENf/x8GqnLuvcDDidb9eFjZ11Vl9cqbpWi/5zHCXdfStWXHxMwITWoWfRYUKSUZcAOIcQBetEMYB3wMXCZXnYZ8JF+/DFwqT7aczhQ36FrZFmyEgeScfmhHD1yBzHC+GjGtbKVF/mEMu+yLrtAY5/SZr5a1X/iC9TxaOAl7lwWD+Uh/K/yS9b/lM9xX/+BHR9cgLPa/IQq3hytdz5kzl87v0iFzO9+5synLuTKB+dSWr06QtZZByFD8FwLISYBzwMuYAtwOZpIvQMMArYB50opa4QQAngKOBHwAJdLKffrJxEiDL/gbpCdOJh+V+QxOiH82bNcwBFkMEWZvs9lTR0tHsY/pPUoV93+IGpcfNhtChYpA2z2zee26nqaNyYF3cUJCkWQMLqRVzIn4Dz6HvxZ5i1oNnXmIGIaapm/eI/uiwoZP6/i3K23UfHcVqqbLf9/ck+WGuViCElQwo2ZghJJMdmtXqFwojqafOckHEpCe3nmLz8w8HNtsuGKu5+MqE37o8FXyP3q96xaNkhzDoWLjsJyzL34MyMfk9nvi7cZffcV7YIivJL0X1ZHs5C0YZig2JGy+8AsMQGolCpviHUMUzczI3AImc4xCOEgc/lPAARi4yJu075w+3fwMl/x0eZ8ZNlAuul77z6qxLMuibOVrcRXncW9WTEcPPQxGsbmI12RWdi84oRz8QwahaukBeWXpzhl4xJ8nyymttXYsPxoxm6h7IGZYrInLmCkiGWyOopT/u91WjKy2XTV7aiuWNNsqvGt53XxE58WDECWmJvLViYLxo7fwT2Ow8kcej0No/ORcWEQFz/ElXhg+dNc2vA1je96qShfgV+a79sxCLvLEw6SY1LJvfpAxidba4jPAbz9t1UsnT4W7+GnEufIMWROULD41UZWB+bz70AF69flQ13Eqg4amSzIHlnF+UkBThQn4hpyAd70dFqz45Hx3RAZVRulSdxWiWP7d7zg/Sdvb8tB+XwHlbUbe5OIdMQWFKNxCieZ5x/HUUOsN5Id1xrgtUfXcvmNY/EnOhmMi3FqHvliBInOPBRh/CRBv9pIqX8dL4nVLCzuj397HPit+13ZE5kocFT6ET4VcZBKSmoT5yV3HohYipfPSrPxNsQSs7SRlrKNNHlro3KxrR5g+1CMJvOQE5gyeKfZZuyTI9fXAdCY4AQJG/CyQRThoojsgINhMpkBMpsskUuCkolTSUbgCKoVo0ovUvrxqU00yQoqZRlblAq20UyVkPy0PQ91Wwyh+EhiVBWfEtkZxcItkTEKvnVLqVldTilg7mrFfQNbUIDczPGMnVEcljgTI5i5omaf5V6gWAYopg5EHQ4245KQFFDoRwwuHLikgzy5+0LqAaFSImoJIKnGSzMqTah4gYCgXTsU4MD8MlbVD+xxN2dYcxMLl86j35Ent4ewRwoZI3CNPIi4td/RorZEtO6+Sp8XlARnIgnnpZGm1JttSqeMLG2mNK3rbk0AbeS2GZVKWnedEPvwCQWpnWmKSt7ICkqW5/Soy7MlTouWPbuymPf65Xf7/aESSI0hY+BkSrf92GvSLFoZ62a2iRCJR8xgQqp1xaSNbyZlmFb3sLhmHPmtXV+4L/RWydOblhtoUffw56WTmZBnWv19iT4tKDnJwxlz+A5LPwSXT5sk99349C6uDB8KMCavArozWtKBO4cdqB2YNAAgHeAcPh6XEl0Z7qIRK/+WwopDOODUkWQbNNkvXEzdoLWe6pLMCzkHyHAESBres5bci7lDADitqsRAi7qHmhxDZtpI0+rvK/RZQemXMZaDBlo/VHrmCuusC3xAWh0yufutFKl3e57buMxok4K3QQF12BASnD2YAW0TNH1SUBzCgXriABItOqrTkTE7PVSmmNs6aSNRSFKH1vbovXcPHWuwNd1HjVNIyxpjthm9mj4pKNHSOmnjm4nmOWT3ZGRKQ49aKXPyhgEwu8rc5x7Iz7FbKWGkTwpK4NhBUdE6ifFrDtlvJ5jnkN2TRCFJyG/q9vtUvdvz7MalRpvUPTviFFJTh5lqQ2+mzwlKVuJARg8zzznYHQ7b1ABATYq1RieGp9dCbPdbKQ8MHo3LAlM9ZN4AzSlvYzh9TlAcYyZYfmSnjc4iZM0mzRFAZHZ/ktzT+cMBOK6m3GiTuoWaHENaXD9Tbeit9ClBcQgH8YeYn5YgWMZva6I20XrBzAqQnlPX7fcFhPZ1m7PB3G6PVCA2c4ipNvRW+pSgZCTkMyqtzmwzuoWVHLIdyU/w9Kjb8+CgA0hUzW8hBvpnEKdYI1lVb6JPCYoy5kCSg1yUy2ycAc1OqwpKiqL2qNvzZP4IAKbVBpfpP1xIl0JirHWc3b2FPiMoAoFjXPTktjikQJvQVxnEpEAzUIDkrO4nompLY/CsiUFuoHd70gZ2faFNt+gzgpLoTGJAVs+CsszAqg7ZjvRL8IDS/W7PowNHku7v/qqLRiPT0hBENqVCb6fPCEpyyiD69WDpULM4aEsjjXHWHtpMdfqRPVjN49GB2pyao4JYNCucqEku4p0JXV9oEzR9RlDkgIHERUEwW0fMTFkQDPFCIpK672D1KppQ/svsbk+MID4mxVQbeht9RlD8edb0RewLRV8oa65FHbJtKIArpWd5Up4aMJx+vh7mWDEIqYAr3trPONroM4ISkxM9KQAPKtScnWUZ5i2XESxx8T0ThYcGaSvYTqk3dza1SEwztf7eRp8QFJfiIiPdWktj7I+ZK63vkG0jLqZny0q0OLRuzzObzO32EGO9wMFopk8IioKCK0rC7QEO3dxAc0x0/GlSYrw9GukBeDZvKPmt5rYc1cR4e6THQKLjWxsiCa5U0pzRM8ID1nfIthGjqD3+Fv1t8GgADmkwcTg/wpn4ezt9QlAcjjhiDYiQza9s4e3/W8WT/9rAqQsrSHEbv4qc0Gfjzo0SQQkFj0PrbjxtYrdHxjqIVazvq4oW7A5kN6hMdbF+YCIHbndz6bwyLp1X1n6uIjWGuRMzmDchg9rknmdYG1+k5RrZmdU35pm81H8wl5dtM61+qYCiOCA6ZmRYHltQukGrS+Hei4a3v3b5VKaur+O4FTWMLvZw4fxyLpy/a2p+TZKTuRMz+HZiBlWpwQ1bHxcFEbIdSVBUZCyIHjbW7h8ylsvLtjGhqY5VSWldv0FKTq8qYY4ew3LuuCl8l26nIrAKtqCEgDdG4fsJGXw/YVf3JMavctjGemaurGH8Njfn/ljBuT/umgjXEO9g7qQMvpmYQXn63k3tqRvq8fXQyWkGraqCCCGcxO3Uuz0bl3PUIcfu+yIpOadyJ//ctGK34ruHjmV+WnbPK7cxHFtQDMbnVPhxXDo/jts1k9UZUDlkcyPHrahm0tYmzvy5kjN/rmw/73EpzJ2U0e43iRaHLGirFYbKqzmDuKR8++6FUnJh+Q4eL1i5W/Edww7kpdwh7Zn0bayFkBZIydcZQhgTK5+VMICDr8+2VOoCRZUcVKiJzOSCvWNkvA6htWQmZbCtXw8mzESIyoCD9YsH9miZ0jaS/T4KF37BsZOO4dDGWh4qXL3b+VtGTOC1nEFhERHFE6B2xVy8avTMRA8DS6WUk424UUgtFCHE74Er0FbKXQ1cDuQCbwGZwFLgEimlVwgRC7wCHAJUA+dJKYtCqT9YmlqqqfblkRxrbqh3R1RFsHRkCktH7ppLMq6okfve3MrCA1I4fGMDs5ZWM2vprkhSFdpFprB/vCWGPFUZmg1CSs6r2AHAvBXz28tvGDmRt/oNDPtnFH6JT42ukAIr02NBEUIMAG4Axkopm4UQ7wDnA7OAx6SUbwkh/gX8BnhG39dKKUcIIc4HHgTOC/kTBIGKil9af4R85kotHuPvZw5pLxNSMna7m5krajh6XR3Hr6jh+D0ct9+OT+ebSRlsHJAQcZGp97m63TpRpOSa4i3cW7Rur3P9jjw5sp/Bwi30aCRUH4oTiBdC+IAEoBSYDlyon38ZuBdNUE7TjwHeA54SQggZgT6XT/VRVZfMiP7N4a4qJI5eV7dXmRSCtYOTWDs4iSdOG6QXSg4o9jBzRQ3Hrq5lur515IexacydlMHaQYlh/YG2+oNLseBUVX5bXMgft21oL1OBK0cfwidZeaT5vGxa9CWjmpvYlJAcJmv3Rmn0ILFFxSh6LChSymIhxN+B7UAz8BVaF6dOStk2iLgTGKAfDwB26O/1CyHq0bpFuyXFEEJcBVzVU7v2aSsSb00c9DfyruEhqBnGQrAxP5GN+Yn882Q965iUDC9tZubKGo7TWzN7CtTPo1OZOzGDVUOTDPNHuN2d+3diVJWbdmzmth2b2stahcLVBxzMZ1m5u11bF6MNq/9j03JOmHSMIbYFhcmh/72NULo86WitjqFAHfAucGKoBkkp5wBz9DoM+9ehFEswfzXMzgk1QlYICvMSKMxL4NmT8tuLh5Q3M2NFDTNX1jB1Qz1TN9Tz8BmDWTQ6NWSTVcDXuPfQ98AWD0uXfNP+usnh4KoDDmFuRs5+7/dBVh5nRnhBddXTswXgbfZNKF2emcBWKWUlgBDiA+BIIE0I4dRbKflAsX59MTAQ2CmEcAKpaM7ZiODYXopPCmIsmmRpVLEHgIJcY0d0inLieeGEAbxwwoCuL+4mXilQ3Q7EHl2GFL+PhSkZPDpwZLeCzu4afiBnVpUw3NNEYUKS0ebuhVChudncrHG9jVA8lduBw4UQCUIIAcwA1gHzgLP1ay4DPtKPP9Zfo5//NhL+kzbcdduoCVg3pWJ7DlkLjNwES43fidiHW2ptUiqnTjiy2xGsNTFaa+fJzSsMsK5rREuAZn/0pLWIBnosKFLKRWjO1WVoQ8YKWlflD8DNQogCNB/JC/pbXgAy9fKbgTtCsLvbNPnq2FEfejM/XOzpVI0GqlviQoo/2RefZOZyaGNknoWjvrmvx58YTkijPFLKe4B79ijeAhy2j2tbgHNCqS8UAjJA88pEmG5uhrD98d2B0bVOTH218flY7xx+IKdUlzKk2U1RfKLh9++INDlJdm/E+sEZBuLYvJWWEAOxwkk0pSxwS0GgpuezqjujwqXNsn5s88ourgwNEQB3486w1tEX6VOCUl23mS2e8P7X6wnDSjWH7Ib86FnSobglHuEOjwvsi4wcjmwIb0tScfvx+Gz/idH0KUHxql6q1mSZbcZeRKNDtrIqfP6o24ePByC/xRO2OkRJCX5pfIKsvk6fEhQAsXg1jaq1PnY0rBLYkTpVIVAevixnZbHa0PmjBeHp9gifpKm+KCz37utY65cVAard21lbYa1WigIsGGPdEag9KWpMhubwjvh/m5bNtDA5TR01Hpp8DWG5d1+nzwlKQAaQ87z4LOacnTsp02wTgsInBQ3F4Re/W0ZMACC31dj5V0KFQOlWe/5OmOhzggJQtX0hq+vSzDYDgMEV2g9m7WDrOYv3RWFzAtSG/2tTHKc5qB8uWGXofZVGHzXuHYbe02YXfVJQvKqX1i+kJVopM3T/STRkIPNJQUVRJqiR+e/+Q2omx9dWdH1hkAgV5LZC2xkbRvqkoIB1WikzomiVwM2exIi0Ttr4/chJAPTzGjMjWKn3Ut241ZB72eybPisoXtVL66de3Ca3UmL9koUHGB9xajRuKagsiFzrBGC73u15cI+UkD1BBMC/dY3dOgkzfVZQAMqLF7Nsi/GzcLvL3InWd8iur8lAmDAwsiglndnVZV1f2AXO0nqqPZFNjdAX6dOCEpAB/B8vpajVnJXj8qq1pvyqoeGfqh8K270uPAWRy6LWkRv0bk+mr+f5gJVWSWPxSntkJwL0aUEBqGkpp+STZFMctG0BbaqF1+FpkYKiLf3Ba86PcWu8JrZ/LVzTo/cLFZTCrTT67ERKkaDPCwpARcH3/LxuUMTrtbpDVgVWV2dCpbmCtywpjTN6mMnNWVxLWc3eybBtwoMtKOipDT77jiXVkfVlJLaqLBlhTlciGDa4k2jeaH537PpRkwBI93Uvd4mj3kfNjiV2VyeC2IKi0+RvxPP6Fra0RnaRcqumLCjxxVC1NiuiozqdsVnPgv/AlrVBv0dpUfFuXk6LaiehjiS2oHSg0r2N8rcUKiOQKjKnVnMyLhtuvSHjcr+TzevyoNV8MWljTWIK51YGl79EBEBuWEdNS3nXF9sYii0oe1BSspCCN+LCLipt/hOrOWTdUrCxIM+UIeL98btRBwFaAuz9IQKgbCigoskOYDMDW1D2QcnOnyl4MyGsojLDgikL3FJocTkmO2H3xbpErSV339bOuz1tYlJWuz5SZtnsgS0onVCy48ewikpqc4AVFoo/qQw4WLo5H1li3ZUBNiQkcVH5vif22WJiDWxB2Q8lO36k4NUYClqMXSunDas4ZMv9TtatyYdya38drht1MABJe3R7lFaJsm6TLSYWwNrfIAtQUrqY0ucqWFpj3I8/q14b/lwy0lyHrAoUtMSzYVW+5Xwm+2JNUiqv5gzC2WE5J6UpgH/NUsrqN5pomU0btqAEQUXTVqpfXM78DYMMyZo/fZXmP/E7zHv8LVKwvCaDkpU5YUs2HQ5uGTmRuhgXQgVnuRv3mh+obi412ywbnZDW5elLNPrq8Xz4OctGTmPQqXXkx+x/tGF/mO2QLfc72bQ1F1mmQBQGfQmfxFG4g8rqtfbsYYthC0o3CMgApZu+wff0ALaePJFDhxUT14O1kjOb/KwZFPkMbXWqwqplg3As3ogjz0cgPRYZRW1UoYKjpoXWolVU2DEmlsQWlB5Q5SnG8W4Zy/sfguvUBA7KqOl23zGSDtkWKVhRnk3gw2Kqaj9BInE2FJKdfgDqkMGoCdYd2WlD8QRQirZRWbvRbpVYGFtQekhABigpXUzcc3EsHHU0Kce5GZvU1KWwpDdqXaVFB4Q/0XOLFKyoyEJ+Vk95+TwCMtB+zi/9lNasJb5+C+k5EwkMyEKNtV78idIqcewoo65qPR6/22xzbLpASGndPrQQPehPmES8I570/EMQM2IZl11NsqLu8zqXT2XGyho+nxy+pTzqVIW1xTkoc6spL1++m5B0RpIzmeS0Yaj5A5AJDlO7QkIF4QngKC6hrnazLSThZ6mUcrIRN7IFxWAcwkFW8jDUyWPIGFfFiER3j/ws3aVFCgrcidSszIKlK6lyb+/RLNs4JY7UxAEo/QcTyEhAxkSu1aK0Shw1bvzlRdR5ivGq3ZtdbNNjbEGJBlyKi8y0kfgnDSdpVD0DUxpIU1RiDPhYPimo9DvZXpOGf1UcYnMBNY1bDfsRCgTxzgRSEvNRsvIIpCcgXYqhLRehgvCqOGo9BCp30ugpsVsj5mALSrThEA5SXBnEpw4kkN8fdaBCUm4jcS4f2XHa2jzJioqrw0d2q0p73EtlSzzu5jg8xYk4t3lRKspprNtKk68hIvk+4pQ4EmPTcSVkIZJSUZMSkfFOZIwCeiNmX2Ij2np+AYkISBS3D+HxIBtr8TbX4G6ttVMMmI8tKL0BgSBGiSHOqQ0hx8dlobgS2s/7m+to1ZfMbPG78ak+yyQLarM91hGPU3EhhIIrPgMhdqmKDPjxttYhpYov0IJXbbXUZ7BpxxYUGxsbwzBMUKIorMnGxsbq2IJiY2NjGF0KihDiRSFEhRBiTYeyDCHE10KIzfo+XS8XQognhRAFQohVQoiDO7znMv36zUKIy8LzcWxsbMwkmBbKv4ET9yi7A/hGSjkS+EZ/DXASMFLfrgKeAU2AgHuAKcBhwD1tImRjY9N76FJQpJTzgT2nx54GvKwfvwyc3qH8FamxEEgTQuQCJwBfSylrpJS1wNfsLVI2NjZRTk/n8uRIKduSUJQBOfrxAKBjjr6delln5XshhLgKrXUD0Ar0bMk4c8gCqsw2IkiiyVaILnujyVaAA4y6UciTA6WU0sjhXSnlHGAOgBBiiVHDWZEgmuyNJlshuuyNJltBs9eoe/V0lKdc78qg7yv08mJgYIfr8vWyzsptbGx6ET0VlI+BtpGay4CPOpRfqo/2HA7U612jL4HjhRDpujP2eL3MxsamF9Fll0cI8SYwDcgSQuxEG635P+AdIcRvgG3AufrlnwGzgALAA1wOIKWsEUI8APyiX3e/lDKYPIhzgv8oliCa7I0mWyG67I0mW8FAey0dem9jYxNd2JGyNjY2hmELio2NjWFYVlCEECcKITbqYfx3dP2OsNszUAgxTwixTgixVghxo17e7WkIEbTZIYRYLoT4VH89VAixSLfpbSGESy+P1V8X6OeHmGBrmhDiPSHEBiHEeiHEVIs/29/r34M1Qog3hRBxVnm+pk6XkVJabgMcQCEwDHABK4GxJtuUCxysHycDm4CxwEPAHXr5HcCD+vEs4HO09EOHA4tMsPlm4A3gU/31O8D5+vG/gGv14+uAf+nH5wNvm2Dry8AV+rELSLPqs0ULytwKxHd4rv9jlecLHAMcDKzpUNatZwlkAFv0fbp+nN5l3ZH+4gT5QKYCX3Z4fSdwp9l27WHjR8BxwEYgVy/LBTbqx88CF3S4vv26CNmXjzbPajrwqf6FqQKcez5jtCH8qfqxU79ORNDWVP0HKvYot+qzbYv8ztCf16do00ss83yBIXsISreeJXAB8GyH8t2u62yzapcn6FB9M9CbrAcBi+j+NIRI8ThwO9oSxgCZQJ2U7YvadLSn3Vb9fL1+faQYClQCL+ldtOeFEIlY9NlKKYuBvwPbgVK057UU6z5fCON0mY5YVVAsixAiCXgfuElKudsS41KTctPH4YUQJwMVUsqlZtsSJE60JvozUsqDADe7ZrAD1nm2ALr/4TQ0IcwDEomiya7hfJZWFRRLhuoLIWLQxOR1KeUHenF3pyFEgiOBU4UQRcBbaN2eJ9Bmf7cFM3a0p91W/XwqUB0hW0H777dTSrlIf/0emsBY8dkCzAS2SikrpZQ+4AO0Z27V5wsRmi5jVUH5BRipe81daI6sj800SAghgBeA9VLKRzuc6u40hLAjpbxTSpkvpRyC9uy+lVJeBMwDzu7E1rbPcLZ+fcRaA1LKMmCHEKJt1usMYB0WfLY624HDhRAJ+veizV5LPt992BC+6TKRcmT1wKk0C20kpRD4owXsOQqtmbgKWKFvs9D6wt8Am4G5QIZ+vQD+qdu/Gphskt3T2DXKMwxYjDY14l0gVi+P018X6OeHmWDnJGCJ/nw/RBtZsOyzBe4DNqCl13gViLXK8wXeRPPt+NBaf7/pybMEfq3bXABcHkzddui9jY2NYVi1y2NjYxOF2IJiY2NjGLag2NjYGIYtKDY2NoZhC4qNjY1h2IJiY2NjGLag2NjYGMb/A1bDPPHNMkSyAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 9.32s (17.98s wall)  load_and_combine_stats(data/jit)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<a href='profiles/load_and_combine_stats__data_jit.svg' target='_blank'>profiles/load_and_combine_stats__data_jit.svg</a><br>"
      ],
      "text/plain": [
       "/home/jonathan/code/data-analysis-speedup/profiles/load_and_combine_stats__data_jit.svg"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from pathlib import Path\n",
    "\n",
    "from utils.preparation import get_segmentation, plot\n",
    "\n",
    "combined_stats = load_and_combine_stats(Path(\"data/small\"))\n",
    "\n",
    "segmentation = get_segmentation()\n",
    "plot(segmentation, combined_stats)\n",
    "\n",
    "combined_stats = load_and_combine_stats(Path(\"data/jit\"), max_time=25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "3fa620ae-afed-4337-a3c0-c7295ec31520",
   "metadata": {},
   "outputs": [],
   "source": [
    "@numba.njit\n",
    "def combine_stats_per_segment_jit(\n",
    "    covs: np.ndarray, means: np.ndarray, sizes: np.ndarray\n",
    ") -> np.ndarray:\n",
    "    reshaped_sizes = sizes.reshape((-1, 1, 1))\n",
    "    sizes_sum = sizes.sum()\n",
    "\n",
    "    # bias correction\n",
    "    corrected_covs = covs * ((reshaped_sizes - 1) / reshaped_sizes).astype(covs.dtype)\n",
    "\n",
    "    e = corrected_covs + means.reshape(\n",
    "        (means.shape[0], means.shape[1], 1)\n",
    "    ) * means.reshape((means.shape[0], 1, means.shape[1]))\n",
    "    new_e = (e * reshaped_sizes).sum(axis=(0)) / sizes_sum\n",
    "    new_mean = (means * sizes.reshape((-1, 1))).sum(0) / sizes_sum\n",
    "    new_cov = new_e - np.outer(new_mean, new_mean)\n",
    "\n",
    "    # bias correction\n",
    "    new_cov = new_cov * sizes_sum / (sizes_sum - 1)\n",
    "\n",
    "    return sizes_sum, new_mean, new_cov\n",
    "\n",
    "\n",
    "def combine_stats_jit(chunked_stats):\n",
    "    segments = []\n",
    "    counts = []\n",
    "    centers = []\n",
    "    covs = []\n",
    "    concated_stats = tuple(\n",
    "        map(lambda x: np.concatenate([i for i in x if len(i) > 0]), chunked_stats)\n",
    "    )\n",
    "    sort_idx = np.argsort(concated_stats[0])\n",
    "    all_segments = concated_stats[0][sort_idx]\n",
    "    all_counts = concated_stats[1][sort_idx]\n",
    "    all_centers = concated_stats[2][sort_idx]\n",
    "    all_covs = concated_stats[3][sort_idx]\n",
    "    g_segments, g_idx = np.unique(all_segments, return_index=True)\n",
    "    g_counts = np.split(all_counts, g_idx[1:])\n",
    "    g_centers = np.split(all_centers, g_idx[1:])\n",
    "    g_covs = np.split(all_covs, g_idx[1:])\n",
    "    for segment, count, center, cov in zip(g_segments, g_counts, g_centers, g_covs):\n",
    "        count_i, center_i, cov_i = combine_stats_per_segment_jit(cov, center, count)\n",
    "        segments.append(segment)\n",
    "        counts.append(count_i)\n",
    "        centers.append(center_i)\n",
    "        covs.append(cov_i)\n",
    "    return np.array(segments), np.array(counts), np.array(centers), np.array(covs)\n",
    "\n",
    "\n",
    "@measure_and_compare()\n",
    "def load_and_combine_stats_jit(path):\n",
    "    chunked_stats = load_data_hdf5(path)\n",
    "    return combine_stats_jit(chunked_stats)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "aa2fb5eb-77c2-4d0e-890b-356b03a87bc3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "13.45s (25.87s wall) load_and_combine_stats_jit(data/jit), -30.73% faster (load_and_combine_stats was  9.32s)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<a href='profiles/load_and_combine_stats_jit__data_jit.svg' target='_blank'>profiles/load_and_combine_stats_jit__data_jit.svg</a><br>"
      ],
      "text/plain": [
       "/home/jonathan/code/data-analysis-speedup/profiles/load_and_combine_stats_jit__data_jit.svg"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 8.78s (16.23s wall) load_and_combine_stats_jit(data/jit), 6.17% faster (load_and_combine_stats was  9.32s)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<a href='profiles/load_and_combine_stats_jit__data_jit.svg' target='_blank'>profiles/load_and_combine_stats_jit__data_jit.svg</a><br>"
      ],
      "text/plain": [
       "/home/jonathan/code/data-analysis-speedup/profiles/load_and_combine_stats_jit__data_jit.svg"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "_ = load_and_combine_stats_jit(Path(\"data/jit\"), max_time=25)\n",
    "_ = load_and_combine_stats_jit(Path(\"data/jit\"), max_time=20)"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}