fixed_N_PSs.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import jax\n",
    "import jax.numpy as jnp\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import numpyro\n",
    "import numpyro.distributions as dist\n",
    "from numpyro.infer import MCMC, NUTS, Predictive\n",
    "\n",
    "numpyro.set_host_device_count(4)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x15f67a0a0>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAGdCAYAAABU0qcqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAABFgUlEQVR4nO3df3RU9Zk/8PdMZjIZ8pMgSYgEjIgiIqCAmGLdFlM5rvXomm9r+3V3WeupRzZYkd1Ty56q1dMa156qdQ9CdV2w363L1u2itVatRcVTBESEFUqLilHCjyQqJCExmSQz9/sHa9rhPk+cB276Gabv1zk5Rz+53Pu5M3fmk8l953lCnud5ICIi+hMLu54AERH9eeICRERETnABIiIiJ7gAERGRE1yAiIjICS5ARETkBBcgIiJyggsQERE5EXE9gWOlUikcOHAAxcXFCIVCrqdDRERGnufhyJEjqK6uRjisf87JugXowIEDqKmpcT0NIiI6QS0tLRg/frz6/axbgIqLiwEAFxd9GZFQNO17qSPd4r8JFxdlvP9QxSniuNf+YcbbJ/e8l/HxACBv0mniuHU/Eu3cLY+Vtm2kqlLed0+P6ZiZzmO4fWiPYainVxwfbG074bkkz6nNeB+RD46I432njRbHC947LI53zBorjhft7fONRd9rF7e1Pj+Wa8LyWgOAcGGhOC49P9ZrQrs+xX1PnyJPsHmfOGyZt0Y7H23fqQr5WkmMjYvjkZ5B31j3hIIMZ/e/24+TP50UHUz5xgo+8h8PAKLrtvnGBjGA3+CXQ+/nmqxbgD75tVskFEUklJ/2vdQxC9InwsdsN+z+82LiuKfsQ9o+pMxDk6cc07ofiXbulsdK2zYS1vY9YDpmpvMYbh/qYxhOygcIYC6hSOYv5khevzyu7COinE9eVNleeKUG9fxYrgnLaw0AwsocpefHek1o5y/uW3m8oRzTMm+Ndj7avlPKHJPaNRTxLwja9aPJi8kLUF7UvwBJxwPg+6AAAPjfCqOfdhuFIQQiInKCCxARETkRyrZ2DF1dXSgtLcX84mt9v4KzCBcpv2ft1n4/Lv8O3yIyrsp0TItQlXxvwGv9wLQf6XGxPibhT/m9biaCeLyHoz0XQeg9139T9YMZ8rWaVH4jMuqA/LKL9mb+chy9Rb4HpF0TI/mYW68Jy3WoCeJ8tHlr7x+DB1sz3o82P/V9orJcHA/1+O//AYBX6L+4koXyddhbIf96L9oj/wq79xT/73wjyrVZvNN/vQ0mE1i354fo7OxESUmJ+O8AfgIiIiJHuAAREZETXICIiMgJLkBEROQEFyAiInIi6/4Q9RPhwkLfH2xZEihe0Sh5W+V4WmLFmzfTNxZ9V56HNj8t9WI5n+Tb72a87XCkpJE1TWTdXpqjed4BJJ7UlNUkufRTokK+hgYK8/zbzpL/Wj+2VamyEJf/QC9+yP8HgAAQb0/4xgYr5HSR9qK2XPva86MmI42pMcu+rYnGIM7HmgC1HFMTbjskjmvvZd21/v33lsvPckmL/MfJFtI1CACh7o/9Yyn5D7OPxU9ARETkBBcgIiJyggsQERE5wQWIiIicMIcQ9u/fj1tvvRXPPvssPv74Y5xxxhlYtWoVZs+eDeBoI6I77rgDjzzyCDo6OjBv3jysWLECkydPNh1nsLXNV3nWdCPReBNRChsAQGjDdv+gclM0qFI8QdxEtQYfLPImny6OSzcjAfn8tZvWVpb9JKdPEsfF5xjAwF/NFcdLNr3vG9u3YIK87SFbpSvtJnJcqLqT1yPf6DWXZzLcLFevw5G8xo37lq59S+AHCCZoY2UNMUnXYVQoEwXIwRkAKGpW3j96Mi+DJpUQSiUTQAZvNaZPQIcPH8a8efMQjUbx7LPPYteuXfjBD36A0aP/0Mfi3nvvxYMPPoiVK1di8+bNKCwsxIIFC9DXJ9czIiKiP0+mT0D//M//jJqaGqxatWporLb2D027PM/DAw88gG9/+9u48sorAQA//vGPUVlZiSeffBJf+cpXApo2ERGd7EyfgH7+859j9uzZ+NKXvoSKigqcd955eOSRR4a+39zcjNbWVtTX1w+NlZaWYu7cudi4caO4z0Qiga6urrQvIiLKfaYF6N133x26n/P8889j0aJF+MY3voHHHnsMANDaevSXfpWV6a1yKysrh753rKamJpSWlg591dTIfxRIRES5xbQApVIpnH/++bj77rtx3nnn4YYbbsDXv/51rFy58rgnsGzZMnR2dg59tbS0HPe+iIjo5GG6BzRu3DhMnTo1bezss8/Gz372MwBAVdXRtEZbWxvGjRs3tE1bWxtmzpwp7jMWiyEWU/q1H0NtMiclamZOFbbUy12ElfI6A0I6zmuXf02opsCMiZogUjxa+Q5LWRPtmFrKKmk4T+25HDhdnl/em3vEccvjou0jpKT6pEZbANB14UTfWNUrSmmdD+USKIfPlFNG1T+TSy5JSSOt+ZicdxqGUIpIe51YU5RaYlKilZsayUSnOe2mXLfS6007n6AaA0qNESM9g+K28R37xHHt9RYR3uMOz6kQtx29xb9tKCmX7TmW6RPQvHnzsHv37rSxt956CxMnHn1B1tbWoqqqCuvWrRv6fldXFzZv3oy6ujrLoYiIKMeZPgHdcsst+MxnPoO7774bX/7yl/Haa6/h4YcfxsMPPwwACIVCWLJkCb773e9i8uTJqK2txW233Ybq6mpcddVVIzF/IiI6SZkWoDlz5mDt2rVYtmwZ7rrrLtTW1uKBBx7AtddeO7TNN7/5TfT09OCGG25AR0cHLrroIjz33HMoKPD3Lycioj9f5koIX/ziF/HFL35R/X4oFMJdd92Fu+6664QmRkREuY214IiIyImsbUgnkZJAABCWao0Zmzt5hfKvCA9+xr990X5527JdcjpOmh8AhKrGiuMpIWWmpd20fWhJtZChdloQjcAAOQmlVUjT6ptZEpBWlhp2ABAXGsFp89bk1UTFce0al0hJJQCAck2oTQ33+P/0QW6LZ693qB3Tsm+NpY6bNUmnpfcGA2oMadmH9p4V27jbN6a+TpTnR7uGpPfDksc3idsmpTEvswZ4/AREREROcAEiIiInuAAREZETXICIiMgJLkBERORE1qbg8iadhry89Bpxye27xG2lhItUJwkAek+xnXJM6GgZ/1Cut9RdK6dY4oVyJ06tNpmUZFHTOkriSSPtR0v8WGtzafX3pMSbWmtMSJgBQMfnTxPHR2+RE0LSftTUmLGGn/i8CfXUhlPxSps4riXyJNYOmkF0/wyiriEgz93a+dRCTYFp8zOm3aQ0alA/3WuJVul9IrDHUEhGaqTHJOz1AxlcKvwERERETnABIiIiJ7gAERGRE1yAiIjIiawNIST3vIdQKL1cieWGZkwpaRJrl4+37wtl4vjY//GXWHnv/8qFZAr2yK3A8vrkG+XVG+S7dJYbvZZSJ4D8GEqlfwA9nJBSyhZpIYyiX/vDI4euOEeboijSqxXvyZx2M1ctC2RopoaePnFYLfOjXJ9aSaiUEsCRBNGoTWMtOWNtpBgE8aa4VqLGUN5ruP1ooQVxHwGV1ZIe26AeV+08JdK8Q8kEQwhERJS9uAAREZETXICIiMgJLkBEROQEFyAiInIia1NwFmKyS9k2WZgvj8vhI0R6/GV3Tn0qJmwJfDRN3kfZO1LLJltzL0uzN8CWjlMTQkr6JjldLi3UWy7/PDMoJN60ckZdSqO27lOVfccrxPE+YS5VPXK5nJCSYNNISTVLSg2wv/BCwvOsPceBXBMBpam0uUjpQK2MjJo8M5QiUtN4yjGDKDmkHVPbt/XTgPQcqftWnk9LuS3LvD0vswaN/AREREROcAEiIiInuAAREZETXICIiMgJLkBERORE1qbgwsVFCIfkxNqxui6c6BuLtyfEbXsr5ARb7LC87z3/x594mvisnOAas1PeR/FOOU2mpX4sza20WmNq87EAGlZF35XTPWPflPfdXe9vVHf4TPm57Zok5xfzerXZyI/MuFf9j0vHVLnZXdnTcvMtLX0l1ebS0l6WOl7Hsx+JtZGgJXWp0a43bd5Jw/6158FSG3KkWV5XlmQgYKunZ00vak0kB4VUp6WhYcobyOz4GW1FREQUMC5ARETkBBcgIiJyggsQERE5wQWIiIicyNoUnERLgxTv9Hcc1RIlkXa5O2lveaU4PupA5mu0VgutSOlyGYSg6mdJ1NSU0kUyVCQ/tlIKUKvhpv1M1FN94h1Rta6qWm27UHuXOC7NUKubZ02TWdJuGmsCUpqji4SZ9ZjaeY5kt9UgWJOBlvSiRq0lqdRBtHSVFf+918+OqERElL24ABERkRNcgIiIyAkuQERE5AQXICIiciJrU3DJc2oRiqSnx0IbtovbSqmSkKGOFwBU/lJOzVnEd8jjWmpMTUgJ4y6SPWqtOiXtpiUPpVp9A/GQaS5jdsoJttaL5W6zH/T65zjueTlNpc1bO39LHTON9nxKdQABub6btfOpJTWmnXtQSTWx46axTqGl+6e1U6g6bkiCaazdjaFdn4b0okZ7P5TO01Ne9ynhuWctOCIiympcgIiIyAkuQERE5AQXICIiciJrQwh5v21G3rEN6SwNkYxlLbSbqJabempZHG17ZdwSOLCWepGYy64oxxy4dLY4LpXi6b1YLn00GJcPqSnaI1/CUume3/3DKeK2Z//gQ3Fce96kG7SWMjfaPgCoN5zF58jajNB4I17cRwANDdV9G69lSwjDGuLR5hLE601jbUgnnb9WEkq73rTn84RLdqX6gQx2wU9ARETkBBcgIiJyggsQERE5wQWIiIicMC1A3/nOdxAKhdK+pkyZMvT9vr4+NDY2YsyYMSgqKkJDQwPa2toCnzQREZ38zCm4c845B7/+9a//sIPIH3Zxyy234JlnnsETTzyB0tJSLF68GFdffTU2bNhgnljqSDdSoWjamKVshrUhm7WJl0Qro2JNK0n7MTcqm1QjH3NPi3/fSomNoMr/9NX6SxFVvCL/YCJtCwCtc/PF8VEH5BI9Rfv948Xv94vbDlaUiONQxqPv+q+hIEq0DEe6DkfymNq1bL0O1dSYdH0ak6iDSskhS5kfNR1mTLtZ3oPUxKQxYSil/ayleDQnej4pT36tHcu8AEUiEVRV+U+ys7MTjz76KB5//HHMnz8fALBq1SqcffbZ2LRpEy688ELroYiIKIeZ7wG9/fbbqK6uxumnn45rr70We/fuBQBs3boVAwMDqK+vH9p2ypQpmDBhAjZu3BjcjImIKCeYPgHNnTsXq1evxllnnYWDBw/izjvvxGc/+1ns3LkTra2tyM/PR1lZWdq/qaysRGur/hdJiUQCiURi6P+7urpsZ0BERCcl0wJ02WWXDf339OnTMXfuXEycOBE//elPEY8b/4T9fzU1NeHOO+88rn9LREQnrxOKYZeVleHMM8/EO++8g6qqKvT396OjoyNtm7a2NvGe0SeWLVuGzs7Ooa+WFv9NciIiyj0nVAuuu7sbe/bswd/8zd9g1qxZiEajWLduHRoaGgAAu3fvxt69e1FXV6fuIxaLIRaLZXQ8SzJFS2xoKRFrLS/xmMq4pyXSlO21dI+JkHYD5MdFO561KVesXWnqJ8xlcPokeVtF1WY5VfPBDDkdN/aFTt/Y+1eUmfYd37FPHJdSg1o6zFpTzPKYW5NdFtbzsdafS23flfFcLHXJgGDq1VkbBkqPl3Xe1tRpeOZU/zENj+twx7TU1xT/fYYN6UwL0D/+4z/iiiuuwMSJE3HgwAHccccdyMvLw1e/+lWUlpbi+uuvx9KlS1FeXo6SkhLcdNNNqKurYwKOiIh8TAvQvn378NWvfhUfffQRxo4di4suugibNm3C2LFHfyq4//77EQ6H0dDQgEQigQULFuChhx4akYkTEdHJzbQArVmzZtjvFxQUYPny5Vi+fPkJTYqIiHIfa8EREZETXICIiMiJrO2IKhnJhFAQ3SK1fXs9fRnvQ923MWWkpXWkxJt13xpL59e8N/eI2/bW+5M9gNxVFQCKTqkQxzum+uu4jX9R6TbaLv/xs1YjT+pcKVeks7N0UDXXGjNc49YknbpvIakF2FJw1vpm1vSZRL32DZ2Wg0g6DktIl1q6xFpp+5YSgGGvH8jg7YOfgIiIyAkuQERE5AQXICIicoILEBEROcEFiIiInMjaFFy4uAjhUHqdL1MtImOCayRZEyhSokZL9qgJIaULaxCsj620vTZvLe0mJc8AoOTxTeK4tP/ec8eL2w4qXVhjG3eL4xAShilj8spakzCIZNdIUruWKmk3KQmmXVfWa19Ka5k7ChsFkrwLIL2ovU6C6G6s7ftEEqD8BERERE5wASIiIie4ABERkRNcgIiIyInsDSEUFiIcziyEIJadsZaRMZSj0bYdyRvFWhmMVGGB/A+UhnRiyRBjAy9tLtpNSgstsKHdcNZKvUjlj7QGc9abvxGpLI61IdsINk3TWJrjBRXiscxxJMvIaMzlpgznYymrBOjvH0E0jbM+ttIxLe9vmTak4ycgIiJyggsQERE5wQWIiIic4AJEREROcAEiIiInsjYFN9jaBoSiaWNaGkQqs2FtYhVEKkktR2IsJSLNRSv1oh7TkMBR017avgNIJWkN89TEj/b8aOPK/oMgPj8BNDQcTlD7OVHWtJ9p+4CSd1KyS20CZ00vWpr6BZSWtTbRlGiliFynFPkJiIiInOACRERETnABIiIiJ7gAERGRE1yAiIjIiaxNwUkN6bTklCfUQ0sp+w0JNcIAqGkqMSUyqUbeR9sh5agnTk0AFo0Sx7XEiqV5lLbvIATVIExLCEkN+QZOl1N90XeVGlzavgVqIkm5VlJKozaNJa2kvU60lJmUdtSee2uaypLUCqpJn6XZnfY6ker9AfrjItVBNDeRNLI0rtSuCa1+o7Qf7dOKmAD0+oEMnnp+AiIiIie4ABERkRNcgIiIyAkuQERE5AQXICIiciJrU3ChilMQyouljxk6brZfXCmOxw/JaZ0ipbNosjDfN6alprQEipa0kZJamuT0SeK4NheNlATLK1S6iir7sKaVLEkoa106yzHztG2Nc7Ek0sJKMlI7pmnfhgSgdd+BpRQN6Tjr60f76Vl6rUTau8RtD8+pEMfLX3pPHJcSt9p4WLnerAlDLcFmqXlnSbtZSddhOBVlCo6IiLIXFyAiInKCCxARETnBBYiIiJzI2hDCkamnIBJNv7EX7UmK2+77XNQ3VrlFvs07GA+J41LYALDd5NduAFobNlluFrf9Za04XtIyII4XNJ94uaBUZbn8DUM5I+0Gutp8y9gw0FKOxdrsTdqPNSShlkrSSt1YygKNYBgkKJZrQgvDqI9JT79/H0p4IP7hoDiulW3SwgzS86Zey8aAhxYgCKJJodqoT5i7JbCQ8vzPgXicjPdIREQUIC5ARETkBBcgIiJyggsQERE5wQWIiIicyNoUnMX4l/2JrxdX/6u47dQVfy+OD8Tj4ni8YmLG8yjZ9L44rjYIU0iJmkNT5Pl9XC2n+sp/L6d7+mr9CTYtGSdtCwDxHfvEca28jLitMdWmsaTmrKmhIJqsWRuBJbWkniEFZ9n2eLaXWBu+WUr9aI+VlsaUEq29FTFhS6C3XP4ZXHv9aCW+Kl4R5qck77Sf+oMoixNUetHy+pSOmfLkFO6x+AmIiIic4AJEREROcAEiIiInuAAREZETJ7QA3XPPPQiFQliyZMnQWF9fHxobGzFmzBgUFRWhoaEBbW1tJzpPIiLKMcedgtuyZQt+9KMfYfr06Wnjt9xyC5555hk88cQTKC0txeLFi3H11Vdjw4YNpv13np6HvFhe2ljskLxedn3Bn9g485W/FbfVTvijOXKducn/z1/T6MhEOd2i1khTaPXnpDlGe+WKWKV75HEt9dNxRp5vLHKmXPdq7DYlfaTVfAugjps1xWNpYBfUvMU6c9o5KuMhrRaekryTElLW+nOaIBKJ2ly05njepBrf2OD2XeK2WpKuu1Z+rKSakfH2hLjtR9Pk5nDRXqV23CE563lkmj/tWLxTTvpZrisgmBqGlqZ21n2Ldf28/pFrSNfd3Y1rr70WjzzyCEaPHj003tnZiUcffRT33Xcf5s+fj1mzZmHVqlV49dVXsWnTpuM5FBER5ajjWoAaGxtx+eWXo76+Pm1869atGBgYSBufMmUKJkyYgI0bN4r7SiQS6OrqSvsiIqLcZ/4V3Jo1a/DGG29gy5Ytvu+1trYiPz8fZWVlaeOVlZVobZX/yKqpqQl33nmndRpERHSSM30Camlpwc0334yf/OQnKCiQf0dqtWzZMnR2dg59tbS0BLJfIiLKbqYFaOvWrWhvb8f555+PSCSCSCSC9evX48EHH0QkEkFlZSX6+/vR0dGR9u/a2tpQVSXfpIzFYigpKUn7IiKi3Gf6Fdwll1yCHTt2pI1dd911mDJlCm699VbU1NQgGo1i3bp1aGhoAADs3r0be/fuRV1dnWlipe8mEYmmp1n2XynXZ5JO4uvT5NTdQx2XiONFe+SH4oPz/OPVP5OTI9YuipGePnH88JwKcdyi+1R/2g0AksIH17J35ASgNu/B6ZPE8bCle6yx/pg5HSfVgrPWPBOSWoCcXtQeK421K2Z45lT/PvbIvy1QO4sWyYkvqdaaltTS0lTa+WhzCQnXvnSOANBXIc9bS5kNVvh/iNWSq1Wb5c6dUpdlAKjcIqdOi5r916eWMNMSg9ZacNJjq71OrF2Zpf0E0YH1WKYFqLi4GNOmTUsbKywsxJgxY4bGr7/+eixduhTl5eUoKSnBTTfdhLq6Olx44YXBzZqIiE56gVfDvv/++xEOh9HQ0IBEIoEFCxbgoYceCvowRER0kjvhBejll19O+/+CggIsX74cy5cvP9FdExFRDmMtOCIicoILEBEROZG1HVGjHycRiaSns059Sp5ub7k/lfQI5onbxvfL+5DSYQCQGO1PvWg13/Le3COOS3WvAABKiin+YeY15XpPkc9nUG6gKqZ+Ij1yutBTOjpGlbRbEDXFtJSV9pOSlvoJojOk+tOZcJ5yNmoYhoQdoF9bEkvaDVCeN2V+nnLNquk45RqSUnCtF5WJ2xbtl1OaMWXfBz8jn7+k+1T58R7/spyO07oHS68VtT6ekZYOlGrnaUk1c71DYXtLSo8dUYmIKKtxASIiIie4ABERkRNcgIiIyImsDSFEegYRiaTfHP9ghnzDUAoKnPYv8to6WCjfXOyqkUtvaI3gRMqNW+mGKwBAuXErhQK0G/84d7w4HP9Q3jz6q9d9Y968meK22ry1kkOWG+XWBm5BlO7R9qHeXNXCI0IIRQsPDBbaXmLaTe6+urN8Y1p4JK9HvsZV0mNuLPOjNZ7rq5UDNQOF/pvlVb/pELfVrkNt3+Ne9c9Fe0wSSpmfWLt8PlLjOQAYtXazbyxsbBiohhba5GtCbo0n067xoJoaHi9+AiIiIie4ABERkRNcgIiIyAkuQERE5AQXICIiciJrU3B5Hw8gLy99fRz3qrxtb0XMN6Y1CNv3FTn1Mf5FOfUi0VI5WqkTqzyhuZeWeInv2Gfa94HFn/GNaekjrYlVXuuJl/vQSuVopUSCSOVoJWo0YSV9JO0nqmw7qKQUtdQloKXG/A0GtRScluwqaJavW7GxmVJuqm1WqThe/vtecVybo1RC6v0rysRtJz7dIY5riUGJ+tpUHivtNV70azkdCOG6Va835VrWrnH1tSKV6FHSixotHWdpPicm6VL9QAaVe/gJiIiInOACRERETnABIiIiJ7gAERGRE1yAiIjIiaxNwaXe/D1SofSkUJ6SzCgxpHhOWyNHMzwheQbYapBZU1aWY1pTYNr2Y7cJaSCt7pfyeGvNx7R6YBbWxlkWWqrPXA9LGA8p89ZSYNFe+aWnbS+l4KT0JwAU75SvK63BoKS7Vn7uK15pE8ffU9KlWoKtc5L/tVK6R856anOJt8v196QErPbaVBvMKa9NjVh70LQH/drXxj0hqZc0NmLUrn0pHac1xksJ7x8pL7N6hPwERERETnABIiIiJ7gAERGRE1yAiIjICS5ARETkRNam4MLTpyCcl57y0WpziWklJQVnJXX/1Dp/hpSkjZa+srCmw9T0nqVrqZKoiRg7i0ppOkv6ZjhaUk/av7ZvbXzg0tkZzyO2cbc4rl0rZT1yt1Wts2rHGf4UXNF+OTVmrUkoXZ8lxn1oCbZ9XygTx2OH/WPxD+UEoJYM1Ejnb+0IauiDrNJeP0HtR7r28yafLm6rvQdZ0qWp7bvEcemYoWQCyOD0+QmIiIic4AJEREROcAEiIiInuAAREZETWRtCkErxSE2fAKUMhlJeRiuZAqW8jHQTWb25qN201250KqEFqWTKoHID0NrATZq7WopGOR9rEyupdE9KK0NkeI5HmqnZn1YuRXmOD15UJo7n9cq3v8f+j7+0SfRXr8tzUZ5PLZygPf8S7bqKKPOuedIWKpEcXCDPT2ukKD3m2s35Qa3porK9VqJHegyDKB+l7VvbvzXwZHldqa9v4boKpViKh4iIshgXICIicoILEBEROcEFiIiInOACRERETmRtCi5SVYlIOL0siTV9ZaElhKRsj/V4WhomKZT5AeSGWuZjGhN5kqCOmac1sDPsw0p6zC3nru1Do5U+GqwoEccLDsmla4rf9zcZA5Qma0qDMKlRGaBf49J5ak0HtZ9YtSZ4GinpqTWeE5soQi/NJbEm0rS0m3Z9ikncgF4/I8mSOrXMjw3piIgoq3EBIiIiJ7gAERGRE1yAiIjICS5ARETkRNam4CwsCRRLjbTh9iPRkkNafabou7YGaRKtZpX2k4V0/tZzV5vjKfuRzt+aELLWiLM0pNNY64FJtOe4GLZEnnTMRN1Z4raHzywTx6t/lnmdMGtTO+v2fbX+hpHRnqS4rZQAHI7lGteoDR21Y0rvQQE0bhyO+H5jfC1bHhdLTbpM8RMQERE5wQWIiIic4AJEREROcAEiIiInTAvQihUrMH36dJSUlKCkpAR1dXV49tlnh77f19eHxsZGjBkzBkVFRWhoaEBbW1vgkyYiopNfyPM8uZWh4Omnn0ZeXh4mT54Mz/Pw2GOP4fvf/z62bduGc845B4sWLcIzzzyD1atXo7S0FIsXL0Y4HMaGDRsynlBXVxdKS0vxOVyJyDEdUa2pLIk1sSHt25qk09Ij1lSWxJwmM5yP9rhqXT61tJ80R2vyLqjtLUYyOaQZUOoDvnd53DcWOxwSty3aL9eZ6yuXf96Uaq1pyTNr2k27tqQEl3r9KDXvoHQ9Fvcxwok0y3tQEK8f7ZjW13JQXVuPNej148UjP0FnZydKSuR6iIAxhn3FFVek/f/3vvc9rFixAps2bcL48ePx6KOP4vHHH8f8+fMBAKtWrcLZZ5+NTZs24cILLzyO0yAiolx13PeAkskk1qxZg56eHtTV1WHr1q0YGBhAfX390DZTpkzBhAkTsHHjRnU/iUQCXV1daV9ERJT7zAvQjh07UFRUhFgshhtvvBFr167F1KlT0draivz8fJSVlaVtX1lZidZW/VdNTU1NKC0tHfqqqakxnwQREZ18zAvQWWedhe3bt2Pz5s1YtGgRFi5ciF27dh33BJYtW4bOzs6hr5aWzH+vS0REJy9zKZ78/HycccYZAIBZs2Zhy5Yt+OEPf4hrrrkG/f396OjoSPsU1NbWhqoq/UZsLBZDLBazz5yIiE5qJ1wLLpVKIZFIYNasWYhGo1i3bh0aGhoAALt378bevXtRV1d3whMdjpRkUesWBVATyprIsiZtLLXtgki3BJV2CyLtZ31+LM9FUDXvLOkjLdWmyeuRO0mW7PGn4BL+cmoA9LTb6LfkfX9wnv95HrtN3rf2hmG9JlJCR1Rr2u1PnewCjF1BA6gvOZwgUrRa/Uop7agdT3qOw6kokMFDZVqAli1bhssuuwwTJkzAkSNH8Pjjj+Pll1/G888/j9LSUlx//fVYunQpysvLUVJSgptuugl1dXVMwBERkY9pAWpvb8ff/u3f4uDBgygtLcX06dPx/PPP4wtf+AIA4P7770c4HEZDQwMSiQQWLFiAhx56aEQmTkREJzfTAvToo48O+/2CggIsX74cy5cvP6FJERFR7mMtOCIiciJrG9JFqioRCeenjZluZis3Ir15M8Xx0Ibtpv2I2yo3HbXGZtayJhJzwMFyg1Yr62EMW1j2EVQ5I1OzLmOZFmku6nP/5h5xXDVJ/ju4+CF/eZ1EeZ64rVRaBwB6K+S0qRZOsNCeHy3IIgYLlHPXxgcK88Vx6bUc1HUVRDkwjTWcIF6HxgCG1lwxJJyP5T0l5WV2TfETEBEROcEFiIiInOACRERETnABIiIiJ7gAERGRE1mbghtsbQOOaUinsZRd0dJuWlJNSolYE1yWpAkgJ2qsqRytxIZ4PsYET6pSqQGzXS5KK80xiNQQMEzDNyFhaG4yZkj9WBvPaemwgxeVieNF+5O+sbH/IyeN9s2X953XJ89l3KsJ35jWkE6jpa8S544Xx+PC8+P1yBP0hLI9ADBYKL99xQNIh6klhAJIxQZVikeaS1CN98TXuKERZ8obELf1HT+jrYiIiALGBYiIiJzgAkRERE5wASIiIie4ABERkRNZm4KTWNNnFpaGWpYaYcdD2r81OaOdj8Ta2CsVQIO9IBoDAsOke6TmZtr5aDXvDEk9bR9qGlFJdhUINd8AoFdoMjcYD4nbjjrgieMWgxUl4ri1tt1AoVyvLiakrLpr5esn3u5P6Q17TKEJYPRd+TqxPm/aNWR5zVprRqopWmGO2jNvrnmnNAGUSPP2kgkgg0uFn4CIiMgJLkBEROQEFyAiInKCCxARETnBBYiIiJw4qVJwlq6G1qSaqZunlLDCMLWslBSLJQ2j7UNj6dxoTYFp41p9M+mnHC1lpKX3tOdHTQ4JKR5rbS41BWhI8Gk/4Wn7KIN8bSUq/I9tpGdQ3FbrfKqlyfJ6/DXlpOMBQG+9PL+STe+L49Eefw07APhgVqlv7ONqOdU3GI+L44UH5O3zev01JEsKlZp0O/aJ41ItQQA4opx/0a/9dRC7lW2LdyqpNuWY7V+eJo5X/rLZNyYlAAH9jd7yPqG+RwrXcpK14IiIKJtxASIiIie4ABERkRNcgIiIyImsDSGEi4sQDuWnjWk3bqWbY9aGUtrNOHHf4pb2m9baDXTpRry1rEcQDbUswYzhSHO3lD4Kci4S6/Mj3XS1Pj/qzd+2Q+L4gFCmJvqr18Vti7V5GxQ0y4GamBKG6a07S96+Xb6xXrnDf55tf1krbhtXyhPtWyAXnqn9b//2WtigV2uYp2wvhQ00WgDjva/I13jpHvk8y3/fK45Lc9ceb620kvZeNii8Pi3vKWGvH8ggq8NPQERE5AQXICIicoILEBEROcEFiIiInOACRERETmRtCi51pBupkL+khkRKFFnLy1jKyFiSZMAw6RFDeR21OdqkGnF4cLuc1pEeK60szkj+dGJJHQ5HK1FkaSZn3rexOaBIed6gpOCk9FXImHYLosmclurTUmOHPn+aON53UZlvbFCuuIOuSfKVWPvfcrmX3lP8b2sx5TUrbQvIDfMAoEM5n/iH/rJIsY27xW1H1cildbpPlc+zfY7cvLBkj3/7sUp5pkh7lziupVEHLp3tGwsrCTupQWWKpXiIiCibcQEiIiInuAAREZETXICIiMgJLkBERORE1qbgpFpwWrJJSk5Zm4zBkJrTklpBNTaTUnNqqk8cBZLzZorjISUNEwRLozYtGag9x9bkobi9kjxLKYlBjSVJKNWNA/Sab5bz9ArldJR2PhEotQcNte2s13j5S++J4xZv3yTXiNMSbFIiLaHUqtNIDfMAINor15/7YEa+b6ysUG5Il4zLjfRih+R9xw7J28cP+WvNac0IS96VE2za8ykl+Cx1DVkLjoiIshoXICIicoILEBEROcEFiIiInOACRERETmRvCq6wEOFwerLEVCfMmHiy1PdSt1WOGeqRu0tGDHXpEkrnxoHCPHG8qDnzRJqWprKkDgG9g6hUU20kO5wCSpJwT4tpH5Yup9pzaU2NmZKeyvlYOu1q22vbarRUo1ZjUapLpyW4xr8s1xXrqpFrRfad6U+kdc5OyPPrk7uQFu2R3xoTkBNpRfuFLqzt8jE/miY/JmOEfRzdXn6NJ8r94+Nf6BC33fvX8jUx9n/6xXGp267l9cBacERElNW4ABERkRNcgIiIyAkuQERE5IRpAWpqasKcOXNQXFyMiooKXHXVVdi9O71kQ19fHxobGzFmzBgUFRWhoaEBbW1tgU6aiIhOfqYU3Pr169HY2Ig5c+ZgcHAQ//RP/4RLL70Uu3btQmHh0fTOLbfcgmeeeQZPPPEESktLsXjxYlx99dXYsGGDaWKpnh6kQulJCrWzqJQoUhJCWpJDq+UldS3VkkptSv2oilfkFFxK6boopeZiSjfCgVr5fLqV8aJmcVikJbW05wHd8hylWmNW1pSVNheJWvdMq0snnI+1PqBGTcdJ9ba0xJwxwRYSHivTaw36+Wv7kd54Rm3IvDsnAIzd2imOJyr818RgXE7YafXXkkp31sRoeVyqP6el+jRaaq50j5xSlWiJ23Gvyq+HvB45BSd129U6BJ8I0wL03HPPpf3/6tWrUVFRga1bt+Liiy9GZ2cnHn30UTz++OOYP38+AGDVqlU4++yzsWnTJlx44YXBzZyIiE5qJ3QPqLPz6E8g5eVHf5LfunUrBgYGUF9fP7TNlClTMGHCBGzcuFHcRyKRQFdXV9oXERHlvuNegFKpFJYsWYJ58+Zh2rRpAIDW1lbk5+ejrKwsbdvKykq0tsof0ZuamlBaWjr0VVMj/zEnERHlluNegBobG7Fz506sWbPmhCawbNkydHZ2Dn21tNj+Wp2IiE5Ox1WKZ/HixfjFL36BV155BePH/6FETFVVFfr7+9HR0ZH2KaitrQ1VVfLNyFgshlgss5t1lhud2g1+tRGYUo4mOX2Sf1C5cVf++15xvP3iSnnfSmOqov3+RlMdZ8jlOMre8W873Pa95f6ghDbvaEDN4aTnx9zUT3k+LaWVLA21hiOVrpFu5APHcZ6GJoUjGXzQAg5WltehGjTp8d/gB4Bkob/kDgAcFkrxaA5Pk0MIRbUd4niiWQ4aSQ3ptDI3VZvl16xGCjgAckM+qcTRcLT3Qynco14TQgmycDIBvPnUpx//U7f440l5HhYvXoy1a9fixRdfRG1teqfCWbNmIRqNYt26dUNju3fvxt69e1FXV2c5FBER5TjTJ6DGxkY8/vjjeOqpp1BcXDx0X6e0tBTxeBylpaW4/vrrsXTpUpSXl6OkpAQ33XQT6urqmIAjIqI0pgVoxYoVAIDPfe5zaeOrVq3C3/3d3wEA7r//foTDYTQ0NCCRSGDBggV46KGHApksERHlDtMC5Hny70v/WEFBAZYvX47ly5cf96SIiCj3sRYcERE5kbUN6SzEVJaW1FKSHFoJiyNTpVSJnJhrvVhuKDVmizyVrkny9onR/gTbYFz+9LlvgZLi2SMfs6TF3yhKKxkSaZfL3Gg/tWipLCnZZU1qaaWVgqAlzzQhISFkLUVjnYuYVDMm6SyPeRBJR0BvAihf+bLQhu3iuDdvpjg+dps/kfjBefK1XKG8NjsPyTV3Sg7L249+y5940xrmjXlU/qN87fmR86xA15en+ca0x0rbt6eUIJNSnVqiUXxtenIC0DevjLYiIiIKGBcgIiJyggsQERE5wQWIiIic4AJEREROZG0KLlRxCkJ5x6SzAmhspjUw66tVEh6Cj6bJNdzOmSInft4qk5Mmsd8WieOjDviTbQNXdMj72CCndRKzusXxfaP9KcDTnpFrwWl1paJK3TM1CSWkZ7TUobXplTUJFoSk0PDN2sBNo+1HStlJNekAIE9JNlkeW0tNOkBPAWos6UC1MaBSk1FqxjjueXl+Wp3G8S/K17ilyZzWME9q9gbY6wlK++/4v3LFmfKX3hPHtb/sHDjd//xE2uVWOYNCvczkYB+wSdn5H+EnICIicoILEBEROcEFiIiInOACRERETnABIiIiJ7I2Bee1fwgvlN5l0JJ40hJC1o6BnZP8a3REDo3hd9smiuNFtXIaJiqk3QDgozn+jonRXWXitmP2y1W1ivbLab+BuH8s7025cFxISVOZu80KXUvNteAUlrpn2vlIqbbh9m1hrjOnzNHyeGnno70mxOMp41qKVNu3lryT0oHauWvpMC1JWbxTHpdoddkGLp0tjhc1Z/58at16NdpzHJ45VRyXzr/saTmJq9V8O7hATiNKte2gvHfmCWlEL+mvOSnhJyAiInKCCxARETnBBYiIiJzgAkRERE5wASIiIieyNgVnIaVHtPSNVsfs0F/WiuNSJ9IzP/OeuK2WgjvSISeHipX6blriTVK2S67P1HqRvI+q3/iP+YHQWRHQE0KasFLfS3p+1PpeSu20IBJp2jWhJswm1cjjI9idNQhap9BBZXupxlei7ixx2/iOffIxlXScRnw+jfXntHScJWGn7TuinKd2fUr715KB2rw1WldZrY6bhdQ9FpBr3mkdXgu+5L9+BnsSwNWffnx+AiIiIie4ABERkRNcgIiIyAkuQERE5ETWhhBSR7qRCsk3vY5laW6lNXcq/728j2hvgW/sLZwmbluxU74t2H2q3MQqeihfHC85JJfXkXRMlctjSGEDAEgW+o9Z+ctmcVtPuYmqlXqx0G7mquEEQ8kdAGKAQCuNou5DufmrzdG0b4X22IrlWJRSNFrTtMG43Ejxo8v99Zkm/VQuH6VRb6wrN/+lwIFaKsgQNtBYb/xr1EZ9hmtCew/SwiOhDdvF8SDKWR2aItTmAlD+e3+9sXi7so8n/E39kv2ZlUPiJyAiInKCCxARETnBBYiIiJzgAkRERE5wASIiIieyNgUn0dJuWqrEso+U0FQJACK9/kTRac/Y0kfjXpUTOFIjJ0BOqmnbaqR9WPdjbWyW0pqPGcroqCkj63NvaAYWxDEtzRIBveFZQbPc1E/KV2rNFfctkNOYRXvyxPHYYXFY1HvueHFcm7eWPpMeLzUBqDVqM5RzspTQGW4ulmvf3IzQmHazpDFThf40LwBEe+VrRXqf0BK3cSG1OziQWZKXn4CIiMgJLkBEROQEFyAiInKCCxARETnBBYiIiJzI2hRcuLgI4VB6kiuItJt6vDY5xVMk1ANLTp8kb9ssp15CSs2uvtpycTzSI7QOM9Yl036ykBqHqXXZlPSN2thNm4shraOxztHC2hxPOqY1SYdfvS6PKymr7lr/MaM9SXHbsp3yyzpVL8fdok+X+cb2fcE/BgBF++V0U6xdTllpaTLr61MSSPpVuZY1WqrPE2oPhrXGhdZGh9r2ynuWZd+je+QU4JFp/vHi9+X3sei7/udhMJVZ2pafgIiIyAkuQERE5AQXICIicoILEBEROcEFiIiInMjaFJyFmEpSEkzmtJJArdmkJW2UucSVRI2U7tEqK2kJLi0hZElwWWtQWVJJVloNLi2RJ51TEIk5wFiDS0vSSR1OAXhKWqlIOM+OK87JeB4AMOYhed4Dhf6rq/z3CXHbwULlLUNLaSqPuXStBFHrUTumdR/a9Tao1asTnmc1Xak8Vup7k5J2szyG6jWrvAeNWrtZ3l6ahzTmDWT0b/kJiIiInOACRERETnABIiIiJ8wL0CuvvIIrrrgC1dXVCIVCePLJJ9O+73kebr/9dowbNw7xeBz19fV4++23g5ovERHlCHMIoaenBzNmzMDXvvY1XH311b7v33vvvXjwwQfx2GOPoba2FrfddhsWLFiAXbt2oaBALtchCRcWIhxOL8Vjufmt3USEcqNPu0kp3dQLorQMoN+gluaulTQJokSNtZmadv5quRPLDVrlfLTz9+bNFMe1oIjlmJbzN4dBtLIrCmn/o7e0y8dUGtVpAQKppI/WuFAquwIACCCYoj0PI9mIUp2LsURPUO8JkiBe41IJLkAP8UjvQVoZIlGqH8jg6TEvQJdddhkuu+wy8Xue5+GBBx7At7/9bVx55ZUAgB//+MeorKzEk08+ia985SvWwxERUY4K9B5Qc3MzWltbUV9fPzRWWlqKuXPnYuPGjUEeioiITnKB/h1Qa+vRz1yVlZVp45WVlUPfO1YikUAi8Ye/Oejq6gpySkRElKWcp+CamppQWlo69FVTo5QeJyKinBLoAlRVdfRGX1tbW9p4W1vb0PeOtWzZMnR2dg59tbTYbs4SEdHJKdBfwdXW1qKqqgrr1q3DzJkzARz9ldrmzZuxaNEi8d/EYjHEYjHfuFcYh5eXPq42QhPSIFpiw5r4ErfVUinGhmzafrTzFPdhmDcgp4GsZUqCKLkTxLwBIPXmHnFcKotkKU803Hiq0t9IUGtoaEkGDseS9Iy0K7/GVtJx0vbaNThoLduklBySmjRqr1nrYxVEqSRryS7LHK3vE6YkoXHenrIf6flPGl6zmZbiMS9A3d3deOedd4b+v7m5Gdu3b0d5eTkmTJiAJUuW4Lvf/S4mT548FMOurq7GVVddZT0UERHlMPMC9Prrr+Pzn//80P8vXboUALBw4UKsXr0a3/zmN9HT04MbbrgBHR0duOiii/Dcc8+Z/gaIiIhyn3kB+tznPgfP0z64AaFQCHfddRfuuuuuE5oYERHlNucpOCIi+vPEBYiIiJzI2oZ0XvuH8EL5n76hwlLbDbClsqypNm37UNVYcVyre2ZhqZ9lbdRmbWAn0c7dVG8KACbJfzfmCSmrIBKQAIDtu/xjxn1bG+yJz6eWGhNSegCQpyQGpV+oW69ZTUhJB0qvlVBA9dQsqTHrtR9E/bmgmuNJ7xNBvZaluSenTxK3FesDZlgLjp+AiIjICS5ARETkBBcgIiJyggsQERE5wQWIiIicyNoUXOpIN1KhaEbbSskPawJFS5oMBpBI01hWf2tHR0vdM2saz5omU7vTGqi1tpTElyRRd5Y4Htu4WxzXHhdLrT6Ntg/tMZTG1cST8lgFkTy0JjStqSzTvpU6cykhpajOQ0lRiknHYVhqLKp187R0reF6055j7X1MrREndFDVUpQpacyTO+oei5+AiIjICS5ARETkBBcgIiJyggsQERE5wQWIiIicyNoUnERLj0jJD0tyBLAlvsy1wxRazS7pSbHWj7J0buyrlecRVR4Tay04KWWlnbtXKPeNSpw7Xhz/YIZcLzAp7Gb8i0raS0tCCfXkAGPHTe2xMnbFlFhTVpYEm3b9jGSqzZwa0+rMSWPa87CnRR5XEnZQtpdeV9a0qPV1JY1rz3FQdekynUfG//aEj05ERHQcuAAREZETXICIiMgJLkBERORE1oYQwsVFCB/TkM5SqsJ6Q896Y9DCegN00NIcTzsfQzmW6K9et+3b+JhINzrzhFIfgB5CiPQMiuNF++VLuHOS/2erg5+Rj1n2Tkwcj7fLAYeoFKpQHm9LGGQ4gQQfDCVgrDenreVlLPs3BzkCCEpI5XyGYwkrWV9XQZSEsl5v0nWr7UOad8obyOg4/AREREROcAEiIiInuAAREZETXICIiMgJLkBERORE1qbgkufUIhRJT0SlNmwXt5VSJZayPQAQMaR1giproaVeLM2trEkbKTljTeVYmlgB8k85WhO0wYoScfy9y+PyXHpD4viYnUnfWG+5/PNWtMe/7XCkMkLaT3JBNSWT9mNNe1nSiyNZckfbvzUxp6WyTIlBYzrMwvoYqu8rhhJKQb2WgyjRkwl+AiIiIie4ABERkRNcgIiIyAkuQERE5AQXICIiciJrU3DR99oRCafX4pKrgdlqWVnSbhpzHS9j4smy/6DqtVn2rT1WeZNPF8fF5KGWgiuUL8mSPeIwknI4Tky2xdsTpmNG2rvEcSlJqNXvG8nag9bn2DKXkbyuNNbEoOUxtO5DvZa16zaA1JiW0LW8f1hrElof80yFvX4gg0uFn4CIiMgJLkBEROQEFyAiInKCCxARETnBBYiIiJzI2hTcYGsbEIpmtu0IJtgsRjLtpgkqCWXZt5YQ0jo0WtI6A4V54vjYn+4UxxN1Z4njXTXStSNfT5W/bBbHtdp2lsdcSzYhgNSci6RaUPXntFSWxFq/Ubq2rI9JUqm/pp2/pX6j2iVXef2YavgZU76m623mVPmgSmfnTPATEBEROcEFiIiInOACRERETnABIiIiJ7I2hBAuLkI4lP/pG+JPf9PVeqNPY2kmp978nVQjjys3Bi3NujTaDVqNdJ7ajdiiZvkx0cIGBc2HxPH4Dn/JFC34kBJHAQQQErE+Vhrp8VIbsjkIJ1j3bQkKaD8le6Yj2mhBG+35NJVQCigIJQYFAnh9A/I1lNq+K+NtMz7Ocf9LIiKiE8AFiIiInOACRERETnABIiIiJ0ZsAVq+fDlOO+00FBQUYO7cuXjttddG6lBERHQSGpEU3H/+539i6dKlWLlyJebOnYsHHngACxYswO7du1FRUZHRPsKFhQgf05DOmvoRt9XKd2gJHEMZGWsaxNIkSi3foaTdgkg8WVJ6x7O9uG2bnGqLKecZUp7PgdOF0ihK47nYxt0Zzu4oMSGknbuxfEkQCUv1WjGUxbGWkdHK5VjKM1kFUeLKWhbHUm5LewyDSiOKST2lYZ4miNestI+UN5DRvx2RT0D33Xcfvv71r+O6667D1KlTsXLlSowaNQr/9m//NhKHIyKik1DgC1B/fz+2bt2K+vr6PxwkHEZ9fT02btzo2z6RSKCrqyvti4iIcl/gC9CHH36IZDKJysrKtPHKykq0tvo/kjY1NaG0tHToq6ZG+cNKIiLKKc5TcMuWLUNnZ+fQV0vL8Zf2JiKik0fgIYRTTjkFeXl5aGtrSxtva2tDVZX/hl8sFkMsFhv6f887WmBjMNXv2zbl+ceswim5J8ygdtNshOZhnoswj+Hmot0EDBvmbt2HZXvt3PW5yMdUH8PBPmFMvtzzAng+1XNPJuR/YD2fDG/qHp2MvG/LNW7ZFgC8pNzHyTOcp/V1pT1W0n7Um+LG15VlLqbn7Dh4wrUVsj73Cuk1m2mwAAAGcXTbT97PVd4IuOCCC7zFixcP/X8ymfROPfVUr6mp6VP/bUtLi4ejZZ74xS9+8YtfJ/FXS0vLsO/3IxLDXrp0KRYuXIjZs2fjggsuwAMPPICenh5cd911n/pvq6ur0dLSguLiYhw5cgQ1NTVoaWlBSUnJSEw1K3R1dfE8c8SfwzkCPM9cE/R5ep6HI0eOoLq6etjtRmQBuuaaa/DBBx/g9ttvR2trK2bOnInnnnvOF0yQhMNhjB8/HgAQCoUAACUlJTn95H+C55k7/hzOEeB55pogz7O0tPRTtxmxdgyLFy/G4sWLR2r3RER0knOegiMioj9PWb0AxWIx3HHHHWkpuVzE88wdfw7nCPA8c42r8wx53qfl5IiIiIKX1Z+AiIgod3EBIiIiJ7gAERGRE1yAiIjIiaxegHKtq+orr7yCK664AtXV1QiFQnjyySfTvu95Hm6//XaMGzcO8Xgc9fX1ePvtt91M9jg1NTVhzpw5KC4uRkVFBa666irs3p3e8K2vrw+NjY0YM2YMioqK0NDQ4KsdmO1WrFiB6dOnD/3hXl1dHZ599tmh7+fCOR7rnnvuQSgUwpIlS4bGcuE8v/Od7yAUCqV9TZkyZej7uXCOn9i/fz/++q//GmPGjEE8Hse5556L119/fej7f+r3oKxdgD7pqnrHHXfgjTfewIwZM7BgwQK0t7e7ntpx6+npwYwZM7B8+XLx+/feey8efPBBrFy5Eps3b0ZhYSEWLFiAvj5/Yc1stX79ejQ2NmLTpk144YUXMDAwgEsvvRQ9PX/oUnnLLbfg6aefxhNPPIH169fjwIEDuPrqqx3O2m78+PG45557sHXrVrz++uuYP38+rrzySvz2t78FkBvn+Me2bNmCH/3oR5g+fXraeK6c5znnnIODBw8Off3mN78Z+l6unOPhw4cxb948RKNRPPvss9i1axd+8IMfYPTo0UPb/Mnfg06k6OhIuuCCC7zGxsah/08mk151dXVGBU1PBgC8tWvXDv1/KpXyqqqqvO9///tDYx0dHV4sFvP+4z/+w8EMg9He3u4B8NavX+953tFzikaj3hNPPDG0ze9+9zsPgLdx40ZX0wzE6NGjvX/913/NuXM8cuSIN3nyZO+FF17w/uIv/sK7+eabPc/Lnefyjjvu8GbMmCF+L1fO0fM879Zbb/Uuuugi9fsu3oOy8hOQtatqLmhubkZra2vaOZeWlmLu3Lkn9Tl3dnYCAMrLywEAW7duxcDAQNp5TpkyBRMmTDhpzzOZTGLNmjXo6elBXV1dzp1jY2MjLr/88rTzAXLruXz77bdRXV2N008/Hddeey327t0LILfO8ec//zlmz56NL33pS6ioqMB5552HRx55ZOj7Lt6DsnIBsnZVzQWfnFcunXMqlcKSJUswb948TJs2DcDR88zPz0dZWVnatifjee7YsQNFRUWIxWK48cYbsXbtWkydOjWnznHNmjV444030NTU5Pterpzn3LlzsXr1ajz33HNYsWIFmpub8dnPfhZHjhzJmXMEgHfffRcrVqzA5MmT8fzzz2PRokX4xje+gcceewyAm/egEStGStTY2IidO3em/T49l5x11lnYvn07Ojs78V//9V9YuHAh1q9f73pagWlpacHNN9+MF154AQUFBa6nM2Iuu+yyof+ePn065s6di4kTJ+KnP/0p4vG4w5kFK5VKYfbs2bj77rsBAOeddx527tyJlStXYuHChU7mlJWfgKxdVXPBJ+eVK+e8ePFi/OIXv8BLL7001F4DOHqe/f396OjoSNv+ZDzP/Px8nHHGGZg1axaampowY8YM/PCHP8yZc9y6dSva29tx/vnnIxKJIBKJYP369XjwwQcRiURQWVmZE+d5rLKyMpx55pl45513cua5BIBx48Zh6tSpaWNnn3320K8bXbwHZeUClJ+fj1mzZmHdunVDY6lUCuvWrUNdXZ3DmY2c2tpaVFVVpZ1zV1cXNm/efFKds+d5WLx4MdauXYsXX3wRtbW1ad+fNWsWotFo2nnu3r0be/fuPanOU5JKpZBIJHLmHC+55BLs2LED27dvH/qaPXs2rr322qH/zoXzPFZ3dzf27NmDcePG5cxzCQDz5s3z/UnEW2+9hYkTJwJw9B40ItGGAKxZs8aLxWLe6tWrvV27dnk33HCDV1ZW5rW2trqe2nE7cuSIt23bNm/btm0eAO++++7ztm3b5r3//vue53nePffc45WVlXlPPfWU9+abb3pXXnmlV1tb6/X29jqeeeYWLVrklZaWei+//LJ38ODBoa+PP/54aJsbb7zRmzBhgvfiiy96r7/+uldXV+fV1dU5nLXdt771LW/9+vVec3Oz9+abb3rf+ta3vFAo5P3qV7/yPC83zlHyxyk4z8uN8/yHf/gH7+WXX/aam5u9DRs2ePX19d4pp5zitbe3e56XG+foeZ732muveZFIxPve977nvf32295PfvITb9SoUd6///u/D23zp34PytoFyPM871/+5V+8CRMmePn5+d4FF1zgbdq0yfWUTshLL70k9k1fuHCh53lHY5C33XabV1lZ6cViMe+SSy7xdu/e7XbSRtL5AfBWrVo1tE1vb6/393//997o0aO9UaNGeX/1V3/lHTx40N2kj8PXvvY1b+LEiV5+fr43duxY75JLLhlafDwvN85RcuwClAvnec0113jjxo3z8vPzvVNPPdW75pprvHfeeWfo+7lwjp94+umnvWnTpnmxWMybMmWK9/DDD6d9/0/9HsR2DERE5ERW3gMiIqLcxwWIiIic4AJEREROcAEiIiInuAAREZETXICIiMgJLkBEROQEFyAiInKCCxARETnBBYiIiJzgAkRERE5wASIiIif+P1dyDHd4TZK7AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_xy = 64  # Number of grid points in each direction\n",
    "\n",
    "# Deterministic jax simulator that takes in an array of positions x, y and puts down point sources at those positions\n",
    "# on a grid (-1, 1) x (-1, 1)\n",
    "@jax.jit\n",
    "def simulate(x, y, s):\n",
    "    \"\"\" Simulate a map with point sources at positions x, y with fluxes s.\n",
    "    \"\"\"\n",
    "    psf_std = 0.1  # PSF width\n",
    "    mu_bkg = 1.  # Background level\n",
    "\n",
    "    grid = jnp.linspace(-1, 1, n_xy)\n",
    "    xx, yy = jnp.meshgrid(grid, grid)\n",
    "    mu = jnp.zeros_like(xx)\n",
    "    \n",
    "    for i in range(len(x)):\n",
    "        mu += s[i] * jnp.exp(-((xx - x[i])**2 + (yy - y[i])**2) / psf_std ** 2)\n",
    "    return mu + mu_bkg\n",
    "\n",
    "# Simulate data\n",
    "mu, sigma = 10, 2\n",
    "N = 20\n",
    "\n",
    "\n",
    "x, y = jax.random.uniform(jax.random.PRNGKey(0), (2, N), minval=-1, maxval=1)  # Random positions\n",
    "s = jax.random.normal(jax.random.PRNGKey(42), (N,)) * sigma + mu  # Fluxes drawn from Gaussian with POI mu, sigma\n",
    "\n",
    "# Simulate map and Poisson fluctuate\n",
    "mu_map = simulate(x, y, s)\n",
    "data = jax.random.poisson(key=jax.random.PRNGKey(42), lam=mu_map, shape=(n_xy, n_xy))\n",
    "\n",
    "# Plot\n",
    "plt.imshow(data, origin='lower')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## NumPyro model for fixed $N$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def model(data):\n",
    "\n",
    "    # Parameters of interest\n",
    "    mu = numpyro.sample('mu', dist.Uniform(3, 20))\n",
    "    sigma = numpyro.sample('sigma', dist.Uniform(0.1, 4.))\n",
    "\n",
    "    x_ary, y_ary, s_ary = [], [], []\n",
    "\n",
    "    # Loop over point sources, positions and fluxes are parameters\n",
    "    for i in range(N):\n",
    "\n",
    "        x = numpyro.sample('x_{}'.format(i), dist.Uniform(-1, 1))\n",
    "        y = numpyro.sample('y_{}'.format(i), dist.Uniform(-1, 1))\n",
    "        s = numpyro.sample('s_{}'.format(i), dist.LeftTruncatedDistribution(dist.Normal(mu, sigma), low=0.))\n",
    "\n",
    "        x_ary.append(x)\n",
    "        y_ary.append(y)\n",
    "        s_ary.append(s)\n",
    "\n",
    "    x_ary = jnp.array(x_ary)\n",
    "    y_ary = jnp.array(y_ary)\n",
    "    s_ary = jnp.array(s_ary)\n",
    "\n",
    "    # Sort by brightness (label switching)\n",
    "    idx_sort = jnp.argsort(s_ary)\n",
    "    \n",
    "    x_ary = x_ary[idx_sort]\n",
    "    y_ary = y_ary[idx_sort]\n",
    "    s_ary = s_ary[idx_sort]\n",
    "    \n",
    "    # Deterministic map sim\n",
    "    mu_map = simulate(x_ary, y_ary, s_ary)\n",
    "\n",
    "    # Poisson likelihood\n",
    "    return numpyro.sample('obs', dist.Poisson(mu_map), obs=data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "sample: 100%|██████████| 2500/2500 [01:52<00:00, 22.15it/s, 127 steps of size 2.69e-02. acc. prob=0.90] \n"
     ]
    }
   ],
   "source": [
    "mcmc = MCMC(NUTS(model), num_warmup=500, num_samples=2000, num_chains=1)\n",
    "mcmc.run(jax.random.PRNGKey(0), data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "                mean       std    median      5.0%     95.0%     n_eff     r_hat\n",
      "        mu      9.52      0.75      9.55      8.36     10.76     64.69      1.07\n",
      "       s_0      7.51      2.49      8.24      1.87     10.07     39.77      1.14\n",
      "       s_1      8.02      1.76      8.52      5.16     10.21      9.03      1.17\n",
      "      s_10     17.34      0.81     17.33     16.04     18.71   2013.46      1.00\n",
      "      s_11     11.61      0.89     11.59     10.02     12.89   1319.41      1.00\n",
      "      s_12      9.27      0.60      9.28      8.27     10.26   2175.57      1.00\n",
      "      s_13      7.08      3.62      6.79      2.21     13.46     26.22      1.15\n",
      "      s_14     10.11      1.81     10.53      8.23     12.70     27.99      1.10\n",
      "      s_15      6.84      1.89      7.26      3.42      9.48      9.43      1.06\n",
      "      s_16      9.06      0.62      9.05      8.04     10.03   3643.52      1.00\n",
      "      s_17      8.25      1.19      8.25      6.43     10.09    291.92      1.01\n",
      "      s_18      9.71      1.36      9.83      7.99     11.79     82.62      1.02\n",
      "      s_19      7.67      1.61      7.90      5.87     10.27     10.18      1.32\n",
      "       s_2     10.53      0.71     10.50      9.29     11.64   2131.21      1.00\n",
      "       s_3      8.72      2.42      9.64      4.14     11.45     22.05      1.02\n",
      "       s_4      8.47      2.26      8.81      5.27     13.33     27.34      1.04\n",
      "       s_5      9.56      0.65      9.56      8.53     10.67   2433.03      1.00\n",
      "       s_6      7.57      0.75      7.56      6.32      8.71    796.15      1.00\n",
      "       s_7     10.05      0.63     10.04      8.92     11.00   4038.04      1.00\n",
      "       s_8     14.72      0.96     14.69     13.14     16.23   1200.87      1.01\n",
      "       s_9      8.70      2.63      8.96      4.57     13.44     27.20      1.03\n",
      "     sigma      3.11      0.48      3.12      2.41      3.96     20.73      1.14\n",
      "       x_0      0.46      0.19      0.37      0.35      0.82      8.41      1.26\n",
      "       x_1      0.52      0.07      0.55      0.36      0.57      5.87      1.22\n",
      "      x_10     -0.90      0.00     -0.90     -0.90     -0.89   2704.19      1.00\n",
      "      x_11     -0.58      0.01     -0.58     -0.59     -0.57   1361.44      1.00\n",
      "      x_12      0.18      0.01      0.18      0.17      0.19   2660.53      1.00\n",
      "      x_13      0.55      0.19      0.57      0.40      0.79     18.40      1.06\n",
      "      x_14      0.13      0.01      0.13      0.12      0.14    821.34      1.00\n",
      "      x_15      0.18      0.22      0.05      0.03      0.56      3.84      1.45\n",
      "      x_16     -0.22      0.01     -0.22     -0.23     -0.21   2513.32      1.00\n",
      "      x_17     -0.22      0.01     -0.22     -0.24     -0.21   1819.21      1.00\n",
      "      x_18     -0.30      0.01     -0.30     -0.31     -0.28   1015.93      1.00\n",
      "      x_19      0.54      0.40      0.88      0.03      0.90      2.63      2.74\n",
      "       x_2     -0.01      0.01     -0.01     -0.02     -0.00   2737.65      1.00\n",
      "       x_3      0.78      0.03      0.78      0.75      0.81    159.30      1.00\n",
      "       x_4      0.41      0.06      0.41      0.34      0.46     15.45      1.02\n",
      "       x_5     -0.56      0.01     -0.56     -0.57     -0.55   4061.68      1.00\n",
      "       x_6      0.97      0.01      0.97      0.95      0.99    634.05      1.00\n",
      "       x_7     -0.80      0.00     -0.80     -0.81     -0.80   2136.28      1.00\n",
      "       x_8     -0.74      0.01     -0.74     -0.75     -0.73   1489.21      1.00\n",
      "       x_9      0.51      0.08      0.55      0.41      0.60      3.16      1.98\n",
      "       y_0     -0.07      0.22      0.10     -0.38      0.12      3.37      2.11\n",
      "       y_1      0.23      0.05      0.26      0.11      0.27      6.04      1.21\n",
      "      y_10      0.25      0.00      0.25      0.24      0.25   3522.02      1.00\n",
      "      y_11     -0.45      0.01     -0.45     -0.46     -0.44   2208.75      1.00\n",
      "      y_12     -0.61      0.01     -0.61     -0.62     -0.60   3457.60      1.00\n",
      "      y_13     -0.30      0.22     -0.37     -0.46     -0.20     23.03      1.01\n",
      "      y_14     -0.12      0.01     -0.12     -0.14     -0.11    562.40      1.00\n",
      "      y_15      0.15      0.06      0.12      0.10      0.26      3.38      1.58\n",
      "      y_16     -0.87      0.01     -0.87     -0.88     -0.86   2446.15      1.00\n",
      "      y_17      0.59      0.01      0.59      0.57      0.61   1284.09      1.00\n",
      "      y_18      0.72      0.01      0.72      0.70      0.73    809.63      1.00\n",
      "      y_19      0.58      0.45      0.96      0.09      0.99      2.80      2.79\n",
      "       y_2      0.93      0.01      0.93      0.92      0.94   3052.27      1.00\n",
      "       y_3     -0.31      0.01     -0.31     -0.32     -0.29    319.35      1.00\n",
      "       y_4     -0.14      0.17     -0.23     -0.28      0.13      5.32      1.42\n",
      "       y_5      0.25      0.01      0.25      0.24      0.26   2794.06      1.00\n",
      "       y_6     -0.53      0.01     -0.53     -0.55     -0.52   1148.03      1.00\n",
      "       y_7      0.76      0.01      0.76      0.75      0.77   2265.14      1.00\n",
      "       y_8     -0.49      0.00     -0.49     -0.50     -0.49   2578.49      1.00\n",
      "       y_9     -0.32      0.07     -0.36     -0.40     -0.22      2.74      2.70\n",
      "\n",
      "Number of divergences: 0\n"
     ]
    }
   ],
   "source": [
    "mcmc.print_summary()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Posterior predictive maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Predictive distribution\n",
    "predictive = Predictive(model, mcmc.get_samples())\n",
    "samples_predictive = predictive(jax.random.PRNGKey(321), None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {