neutra_nuts_neutrino_experiments.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3268b88a-3e13-46b0-9619-15a91d2925fc",
   "metadata": {},
   "source": [
    "# NeuTra NUTS experiments for neutrino NSI problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "15807d83-ec52-4f1b-b7e8-4f244d391a07",
   "metadata": {},
   "outputs": [],
   "source": [
    "import arviz as az\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import stats\n",
    "import pickle\n",
    "import pandas as pd\n",
    "import optax\n",
    "\n",
    "import json as js\n",
    "from tqdm import tqdm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "76fb6d44-80dd-4b49-97a1-e1587f633199",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "os.environ[\"CUDA_DEVICE_ORDER\"]=\"PCI_BUS_ID\"\n",
    "os.environ[\"CUDA_VISIBLE_DEVICES\"]=\"1\" #select GPU, -1 means use CPU\n",
    "\n",
    "from jax import config\n",
    "config.update(\"jax_enable_x64\", True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ec94a8be-9a5f-4b8f-b677-b19d182566e1",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2024-05-19 18:36:49.602800: W external/xla/xla/service/gpu/nvptx_compiler.cc:679] The NVIDIA driver's CUDA version is 11.2 which is older than the ptxas CUDA version (11.8.89). Because the driver is older than the ptxas version, XLA is disabling parallel compilation, which may slow down compilation. You should update your NVIDIA driver or use the NVIDIA-provided CUDA forward compatibility packages.\n"
     ]
    }
   ],
   "source": [
    "from snuddd import config\n",
    "from snuddd.targets import Nucleus, electron_xe, Electron\n",
    "from snuddd.models import GeneralNSI\n",
    "from snuddd.jinterp import interp1d\n",
    "\n",
    "import numpyro\n",
    "from numpyro import optim\n",
    "from numpyro.diagnostics import hpdi\n",
    "import numpyro.distributions as dist\n",
    "from numpyro import handlers\n",
    "from numpyro import infer\n",
    "from numpyro.infer import MCMC, NUTS, autoguide, reparam\n",
    "\n",
    "import jax.numpy as jnp"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f5c7728-1d85-42c0-8bd3-eb75ea1874c6",
   "metadata": {},
   "source": [
    "## Preparation for SNuDD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "6157ab0f-3cfa-4dfe-bda4-07cca4473d4d",
   "metadata": {},
   "outputs": [],
   "source": [
    "Xe_nucleus = Nucleus(54, 132, mass=131.9041535 * config.u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "11862eba-36bd-4852-9b1c-42a7aa8e9384",
   "metadata": {},
   "outputs": [],
   "source": [
    "SM_matrix = np.array([[0, 0, 0],\n",
    "                      [0, 0, 0],\n",
    "                      [0, 0, 0]])\n",
    "\n",
    "SM_eta = 0\n",
    "SM_phi = 0\n",
    "\n",
    "sm_model = GeneralNSI(SM_matrix, SM_eta, SM_phi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "cfde6386-4b0a-4b28-87c8-c27636989426",
   "metadata": {},
   "outputs": [],
   "source": [
    "Xe_nucleus.update_model(sm_model)\n",
    "Xe_nucleus.prepare_density()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "14828b57-a23c-4242-ac9c-2f540b647b46",
   "metadata": {},
   "outputs": [],
   "source": [
    "E_Rs = np.logspace(-1, 1, 100) / 1e6 # Recoil Energy in GeV \n",
    "\n",
    "sm_spec=Xe_nucleus.spectrum(E_Rs, nu = '8B')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "31809d25-9ca1-483a-a424-e97e77bd6669",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAHTCAYAAAD7+WZaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABxfklEQVR4nO3deVhUZf8G8PvMwLAPuyAIiisoaoq4K2CKaZmaVrYQroW+Zcqrlm2mlZSWy1uSWQaaubRpWZaZu6L+3DDNFVMBAdn3dYbz+wOZGlkEzsDMwP25rrl0znnmPF+WcW6f85znCKIoiiAiIiIinZHpuwAiIiKi5oYBi4iIiEjHGLCIiIiIdIwBi4iIiEjHGLCIiIiIdIwBi4iIiEjHGLCIiIiIdIwBi4iIiEjHGLCIiIiIdIwBi4iIiEjHGLCIdEAQhDo9Dhw4gLfffhuCICA9PV3fZTdITEwM3n77bWRnZ2ttj46OhiAIuHnzZqP0W9/j//nnn5gyZQq8vLxgbm4Oa2tr9O7dG8uWLUNmZmaV49b2M2tInebm5rh161aV/YGBgfD19dXaduLECYwfPx6enp4wMzODi4sLBgwYgP/+97/16ruhavqZStUYvxNTp06FmZkZzp8/X2Xf+++/D0EQsHPnTp31R9RQJvougKg5OHbsmNbzd955B/v378e+ffu0tnft2rXeH9iGJiYmBosXL8bkyZNhZ2en2f7www/j2LFjaN26tf6Ku+vzzz/HrFmz0KVLF8yfPx9du3ZFWVkZTp06hbVr1+LYsWPYvn271muioqLg7e1d5Vhdu3ZtUA0lJSV444038NVXX9Xa7pdffsGjjz6KwMBALFu2DK1bt0ZycjJOnTqFrVu34qOPPmpQ//VR089Uqsb4nVi1ahX27t2L0NBQnDhxAqampgCA8+fPY9GiRZg8eTLGjBmjs/6IGooBi0gH+vfvr/Xc2dkZMpmsynZDVFhYCEtLS8nHcXZ2hrOzsw4qkubYsWOYOXMmRowYgR07dsDMzEyzb8SIEfjvf/+L3377rcrrfH190adPH53V8dBDD2Hz5s2YN28eevbsWWO7ZcuWwcvLC7t374aJyT//JE+aNAnLli1rUN+6+plK1Ri/E0qlEuvXr0dwcDDeffddLF68GGVlZQgJCYGLiwtWrVql0/6IGoqnCIn05M6dO3jqqadga2sLFxcXTJ06FTk5OVptrl27hqeffhqtWrWCmZkZfHx8sGbNmirHOnLkCB588EHY2NjA0tISAwcOxC+//FKlXeXpyTNnzmDixImwt7dHhw4d6tzf22+/jfnz5wMAvLy8tE6j1XQ66PLly3jqqafg4uICMzMzeHp64rnnnkNJSQkAIC4uDlOmTEGnTp1gaWkJd3d3jBkzptpTQHWxdOlSCIKAdevWaYWrSgqFAo8++mi9j7ts2TIIgoCrV69qtl27dg0mJibYtGlTlfYLFiyAo6MjXnnllVqPm5GRAScnJ61wVUkmu/8/0bX9TOvyva3tZ/rvr7Muv4f3qu53orLev/76676//zUZPnw4wsLCsHTpUpw+fRpvv/02zp07h/Xr18PW1larbUNrJ5KKAYtITyZMmIDOnTvj+++/x6uvvorNmzdj7ty5mv0XL16Ev78/Lly4gI8++gg///wzHn74YcyePRuLFy/WtDt48CCGDRuGnJwcrF+/Hlu2bIGNjQ3GjBmDbdu2Vdv3Y489ho4dO+Lbb7/F2rVr69zf9OnT8dJLLwEAfvjhBxw7dgzHjh1D7969q+3n3Llz8Pf3x/Hjx7FkyRL8+uuviIiIQElJCUpLSwEASUlJcHR0xPvvv4/ffvsNa9asgYmJCfr164crV67U63uqVquxb98++Pn5wcPDo96vValUWg+1Wq3ZX3m8+Ph4zbZOnTqhT58+2LFjR5Xj2djY4I033sDu3burnCr+twEDBuDEiROYPXs2Tpw4gbKysnrVXam6n2ldvrf3+5nW9fewvu73+38/y5cvh6enJyZOnIgPPvgAYWFhGDFihFabxqqdqE5EItK50NBQ0crKqtp9ixYtEgGIy5Yt09o+a9Ys0dzcXCwvLxdFURRHjhwptmnTRszJydFq9+KLL4rm5uZiZmamKIqi2L9/f7FVq1ZiXl6epo1KpRJ9fX3FNm3aaI73777feuutKnXVtb/ly5eLAMQbN25otYuKiqqyfdiwYaKdnZ2Ymppa7feiOiqVSiwtLRU7deokzp07t9bj3yslJUUEIE6aNKnO/VUet7qHXC7XtDt8+LAIQFy/fr3W68eNGyf26dOnyvFOnjwplpSUiO3btxf79Omj+TkEBASI3bp107RPT08XBw8erOnT1NRUHDhwoBgREaH1M61JbT/Te9X0va3pZyqKdf+9qE51P7O6/v7XxebNm0UAoqura7XfKym1E0nFESwiPbn3NFWPHj1QXFyM1NRUFBcXY+/evRg/fjwsLS21RlVGjx6N4uJiHD9+HAUFBThx4gQmTpwIa2trzbHkcjlCQkKQmJhY7SjQhAkTtJ7Xtb/6KCwsxMGDB/HEE0/UOg9HpVJh6dKl6Nq1KxQKBUxMTKBQKHDt2jVcunSpXn1KsXHjRpw8eVLrceLECc3+yhGshIQErdddunQJ7u7u1R5ToVDg3XffxalTp/DNN99U28bR0RGHDx/GyZMn8f7772Ps2LG4evUqFi5ciO7du9f5atN7f6aA9O9tY/xeVKrt978uysvL8fHHH0MmkyE1NRXnzp1rstqJ6oIBi0hPHB0dtZ5XzhcqKipCRkYGVCoVPv74Y5iammo9Ro8eDQBIT09HVlYWRFGs9iotNzc3ABVzfO51b/u69lcfWVlZUKvVaNOmTa3twsPD8eabb2LcuHHYuXMnTpw4gZMnT6Jnz54oKiqqV59OTk6wtLTEjRs36vU6APDx8UGfPn20Hn5+fpr97u7ukMlkWqcIv/32W1y5cgUhISE1HnfSpEno3bs3Xn/99VpP//Xp0wevvPIKvv32WyQlJWHu3Lm4efNmnSe6V/c7IPV72xi/F5Vq+/2viw8//BDHjh3D5s2b0alTJ0ydOlXrtY1ZO1Fd8CpCIgNkb2+vGYX6z3/+U22byvWdZDIZkpOTq+xPSkoCUBE67iUIQoP6qw8HBwfI5XIkJibW2m7Tpk147rnnsHTpUq3t6enp9V4yQC6X48EHH8Svv/6KxMTE+4a7+jAxMYGrq6smYMXFxWHmzJkYM2ZMtaNHlQRBwAcffIARI0Zg3bp1derL1NQUixYtwsqVK3HhwoU6vebenykg/XvbGL8XunDx4kW89dZbeO655/Dkk0+ibdu2GDRoEF5//XWsWLECgOHWTi0HAxaRAbK0tERQUBDOnj2LHj16QKFQ1Ni2X79++OGHH/Dhhx/CwsICQMXpk02bNqFNmzbo3LmzTvur60iDhYUFAgIC8O233+K9996rNugBFcHg3qv9fvnlF9y+fRsdO3a8b+33WrhwIXbt2oUZM2bgxx9/rPK1lJWV4bfffmvQWkkeHh5ISEhAcnIyRo8eja5du2LLli33fd3w4cMxYsQILFmypMrk++Tk5GpHnypP4VWORDZEXb+3Nf1M6/N70VRUKhVCQ0Ph5OSE1atXA6hYJiU8PBwrVqzAhAkTMGjQIIOsnVoWBiwiA7V69WoMHjwYQ4YMwcyZM9GuXTvk5eUhLi4OO3fu1FyZFhERgREjRiAoKAjz5s2DQqFAZGQkLly4gC1btlQ7siGlv+7du2vah4aGwtTUFF26dKn2mCtWrMDgwYPRr18/vPrqq+jYsSPu3LmDn376CZ999hlsbGzwyCOPIDo6Gt7e3ujRowdOnz6N5cuXN3j0acCAAfj0008xa9Ys+Pn5YebMmejWrRvKyspw9uxZrFu3Dr6+vlUC1oULF6BSqaocr0OHDpo5ZB4eHjh37hyGDBkCT09P7Ny5E1ZWVnWq64MPPoCfnx9SU1PRrVs3zfaRI0eiTZs2GDNmDLy9vVFeXo7Y2Fh89NFHsLa2xssvv9yg7wOAOn9va/qZ2tjY1Pn3oqlERETg1KlT+PXXX7VG4d555x3s3LkTU6dORWxsLCwsLAyudmph9D3Lnqg5qstVhGlpaVrbq7vi6saNG+LUqVNFd3d30dTUVHR2dhYHDhwovvvuu1qvPXz4sDhs2DDRyspKtLCwEPv37y/u3Lmzzn3Xt7+FCxeKbm5uokwmEwGI+/fvr/Eqv4sXL4qPP/646OjoKCoUCtHT01OcPHmyWFxcLIqiKGZlZYnTpk0TW7VqJVpaWoqDBw8WDx8+LAYEBIgBAQG1fn9qExsbK4aGhoqenp6iQqEQraysxF69eolvvfWW1lWNtV1FCED8/PPPNW3nzp0rAhCffvppTf33+vdVhPd6+umnRQBaVxFu27ZNfPrpp8VOnTqJ1tbWoqmpqejp6SmGhISIFy9evO/XWdvPtK7fW1Gs/mdaqa6/FzV9L6q7irAuv//3io2NFU1NTcUZM2ZUu//YsWOiTCbTukKyobUTSSWIoig2aaIjIiIiauZ4FSERERGRjjFgEREREekYAxYRERGRjjFgEREREekYAxYRERGRjjFgEREREekYFxrVgfLyciQlJcHGxqbOizoSERGRfomiiLy8PLi5uUEm0+2YEwOWDiQlJVW5/QUREREZh4SEBJ3euxRgwNIJGxsbABU/IKVSqedqiIiIqC5yc3Ph4eGh+RzXJQYsHag8LahUKhmwiIiIjExjTO/hJHciIiIiHWPAIiIiItIxBiwiIiIiHeMcLCIiA6RWq1FWVqbvMoiMmqmpKeRyuV76ZsAiIjIgoigiJSUF2dnZ+i6FqFmws7ODq6trk69TyYBFRGRAKsNVq1atYGlpycWLiRpIFEUUFhYiNTUVANC6desm7Z8Bi4jIQKjVak24cnR01Hc5REbPwsICAJCamopWrVo16elCTnInIjIQlXOuLC0t9VwJUfNR+X5q6jmNDFhERAaGpwWJdEdf7ycGLCIiIiIdY8AiIiIyEtHR0bCzszO6Y9fH5MmTMW7cuDq3P3DgAARBMLgrbxmwiIhIstTUVLzwwgvw9PSEmZkZXF1dMXLkSBw7dkzTpl27dhAEAcePH9d67Zw5cxAYGFiv/gIDAyEIAgRBgJmZGdzd3TFmzBj88MMPVdpWthMEAVZWVujUqRMmT56M06dP37efs2fP4pFHHkGrVq1gbm6Odu3a4cknn0R6ejoAw/1wbyyV38d7f4YlJSVwdHSEIAg4cOCAfoozMAxYREQk2YQJE3Du3Dls2LABV69exU8//YTAwEBkZmZqtTM3N8crr7yikz5nzJiB5ORkxMXF4fvvv0fXrl0xadIkPP/881XaRkVFITk5GX/99RfWrFmD/Px89OvXDxs3bqzx+KmpqRg+fDicnJywe/duXLp0CV9++SVat26NwsJCnXwN+lBaWirp9R4eHoiKitLatn37dlhbW0s6bnPDgEVERJJkZ2fjyJEj+OCDDxAUFIS2bduib9++WLhwIR5++GGtti+88AKOHz+OXbt21Xi8wMBAzJkzR2vbuHHjMHnyZK1tlpaWcHV1hYeHB/r3748PPvgAn332GT7//HP88ccfWm0rF5ts164dgoOD8d133+GZZ57Biy++iKysrGrriImJQW5uLr744gv06tULXl5eGDZsGFatWgVPT0/cvHkTQUFBAAB7e3sIgqCp8bfffsPgwYNhZ2cHR0dHPPLII7h+/brm2Ddv3oQgCPjhhx8QFBQES0tL9OzZU2vED6g4befp6QlLS0uMHz8eGRkZWvuvX7+OsWPHwsXFBdbW1vD396/ytbdr1w7vvvsuJk+eDFtbW8yYMaNOx65JaGgotm7diqKiIs22L7/8EqGhoVXanj9/HsOGDYOFhQUcHR3x/PPPIz8/X7NfrVYjPDxc831asGABRFHUOoYoili2bBnat28PCwsL9OzZE999912datUnBiwiIgMmiiIKS1V6edz7QVcTa2trWFtbY8eOHSgpKam1bbt27RAWFoaFCxeivLxcF98iLaGhobC3t6/2VOG95s6di7y8POzZs6fa/a6urlCpVNi+fXu13wsPDw98//33AIArV64gOTkZq1evBgAUFBQgPDwcJ0+exN69eyGTyTB+/PgqX/Prr7+OefPmITY2Fp07d8ZTTz0FlUoFADhx4gSmTp2KWbNmITY2FkFBQXj33Xe1Xp+fn4/Ro0fjjz/+wNmzZzFy5EiMGTMG8fHxWu2WL18OX19fnD59Gm+++Wadjl0TPz8/eHl5ab72hIQEHDp0CCEhIVrtCgsL8dBDD8He3h4nT57Et99+iz/++AMvvviips1HH32EL7/8EuvXr8eRI0eQmZmJ7du3ax3njTfeQFRUFD799FP89ddfmDt3Lp599lkcPHiwTvXqCxcaJSIyYEVlanR9a7de+r64ZCQsFff/mDAxMUF0dDRmzJiBtWvXonfv3ggICMCkSZPQo0ePKu0rPzC//vrrKh/KUslkMnTu3Bk3b968b1tvb28AqLFt//798dprr+Hpp59GWFgY+vbti2HDhuG5556Di4sL5HI5HBwcAACtWrXSmiA+YcIErWOtX78erVq1wsWLF+Hr66vZPm/ePM0o3+LFi9GtWzfExcXB29sbq1evxsiRI/Hqq68CADp37oyYmBj89ttvmtf37NkTPXv21Dx/9913sX37dvz0009aQWbYsGGYN2+e5vlbb71132PXZsqUKfjyyy/x7LPPIioqCqNHj4azs7NWm6+//hpFRUXYuHEjrKysAACffPIJxowZgw8++AAuLi5YtWoVFi5cqPl+rV27Frt3//P7XlBQgBUrVmDfvn0YMGAAAKB9+/Y4cuQIPvvsMwQEBNSpXn3gCBYREUk2YcIEJCUl4aeffsLIkSNx4MAB9O7dG9HR0VXaOjs7Y968eXjrrbckzweqjiiKdVr7qHJUqra27733HlJSUrB27Vp07doVa9euhbe3N86fP1/rsa9fv46nn34a7du3h1KphJeXFwBUGVn6dwCtvJVL5a1dLl26pAkVle59XlBQgAULFqBr166ws7ODtbU1Ll++XKWfPn36aD2vy7Fr8+yzz+LYsWP4+++/ER0djalTp1Zpc+nSJfTs2VMTrgBg0KBBKC8vx5UrV5CTk4Pk5GStfk1MTLRqvXjxIoqLizFixAjNSKm1tTU2btyodcrVENVrBKt9+/aSO5wzZw5mz54t+ThERC2BhakcF5eM1Fvf9WFubo4RI0ZgxIgReOuttzB9+nQsWrSoytwpAAgPD0dkZCQiIyOr7JPJZFVOydV1FW61Wo1r167B39//vm0vXboEAJrwUxNHR0c8/vjjePzxxxEREYFevXrhww8/xIYNG2p8zZgxY+Dh4YHPP/8cbm5uKC8vh6+vb5VAaWpqqvl7ZdCrPI1Yl1O08+fPx+7du/Hhhx+iY8eOsLCwwMSJE6v08++QU9dj16ZyXtm0adNQXFyMUaNGIS8vr0ofNYXXui7+Wfm9+OWXX+Du7q61z8zMrAGVN516jWDdvHkTWVlZEEWxQY9bt261mEtZiYh0QRAEWCpM9PKQugJ2165dUVBQUO0+a2trvPnmm3jvvfeQm5urtc/Z2RnJycma52q1GhcuXKhTnxs2bEBWVlaVU3TVWbVqFZRKJYYPH16nYwOAQqFAhw4dNF+XQqHQ1FgpIyMDly5dwhtvvIEHH3wQPj4+NU6kr03Xrl2rLIdw7/PDhw9j8uTJGD9+PLp37w5XV9c6nR6ty7HvZ+rUqThw4ACee+65au/x17VrV8TGxmr9Dhw9elRzGtfW1hatW7fW6lelUmktn9G1a1eYmZkhPj4eHTt21Hp4eHjUq96mVu85WHPnzsVbb73VoM5kMp6RJCJqbjIyMvD4449j6tSp6NGjB2xsbHDq1CksW7YMY8eOrfF1zz//PFauXIktW7agX79+mu3Dhg1DeHg4fvnlF3To0AErV66s9j/nhYWFSElJgUqlwu3bt/HDDz9g5cqVmDlzpubqvkrZ2dlISUlBSUkJrl69is8++ww7duzAxo0ba1xc8+eff8bWrVsxadIkdO7cGaIoYufOndi1a5dmmYK2bdtCEAT8/PPPGD16NCwsLGBvbw9HR0esW7cOrVu3Rnx8vGauU33Mnj0bAwcOxLJlyzBu3Dj8/vvvVeZIdezYET/88APGjBkDQRDw5ptv1unigboc+34eeughpKWlQalUVrv/mWeewaJFixAaGoq3334baWlpeOmllxASEgIXFxcAwMsvv4z3338fnTp1go+PD1asWKH1s7axscG8efMwd+5clJeXY/DgwcjNzUVMTAysra2rvXLRYIj1IAiCuHjx4vq8RKevN1Q5OTkiADEnJ0ffpRCRESsqKhIvXrwoFhUV6buUeikuLhZfffVVsXfv3qKtra1oaWkpdunSRXzjjTfEwsJCTbu2bduKK1eu1Hrt5s2bRQBiQECAZltpaak4c+ZM0cHBQWzVqpUYEREhjh07VgwNDdW0CQgIEAGIAESFQiG2bt1afOSRR8QffvihSn2V7QCI5ubmYocOHcTQ0FDx9OnTtX5d169fF2fMmCF27txZtLCwEO3s7ER/f38xKipKq92SJUtEV1dXURAETY179uwRfXx8RDMzM7FHjx7igQMHRADi9u3bRVEUxRs3bogAxLNnz2qOk5WVJQIQ9+/fr9m2fv16sU2bNqKFhYU4ZswY8cMPPxRtbW01+2/cuCEGBQWJFhYWooeHh/jJJ5+IAQEB4ssvv1zr970ux67Ov7+Ge1VX/59//ikGBQWJ5ubmooODgzhjxgwxLy9Ps7+srEx8+eWXRaVSKdrZ2Ynh4eHic889J44dO1bTpry8XFy9erXYpUsX0dTUVHR2dhZHjhwpHjx4UBRFUdy/f78IQMzKyqq2rtreV435+S2IYt1PxF6/fh0ODg6wt7dvUJiT+npDlZubC1tbW+Tk5NSY5ImI7qe4uBg3btyAl5cXzM3N9V0OUbNQ2/uqMT+/63WKsEOHDpI6k/p6IiIiImPASVFEREREOsaARURERKRjklZy37t3L/bt24eYmBgkJiYiPT0dlpaWcHZ2Rvfu3REQEIBHHnkErq6uuqq30ZiYmGhW1+3Tpw+++OILPVdERERExqreASs/Px//+9//8PnnnyM+Pl6zWJm5uTkcHBxQVFSECxcu4M8//8TXX38NExMTPProo5g7dy4GDRqk8y9AV+zs7BAbG6vvMoiIiKgZqNcpwrVr16Jjx4544403YGdnh3fffRf79u1Dbm4uCgsLkZiYiIyMDJSVleHy5cvYsGEDnnzySfz+++8YOnQoHnvsMdy4caOxvhYiomahMW6CTNRS6ev9VK9lGkxNTfHMM89g/vz56NatW507KSoqwpYtWxAREYGQkJAGL1Rak0OHDmH58uU4ffo0kpOTsX37dowbN06rTWRkJJYvX47k5GR069YNq1atwpAhQzT7FQoFunfvDgsLC7z33nv1uoEkl2kgIl0oLy/HtWvXIJfL4ezsDIVCIXk1daKWShRFlJaWIi0tDWq1Gp06daqy4LnBLNNw+fLlBi21YGFhgalTpyI0NBSJiYn1fv39FBQUoGfPnpgyZUq1t0fYtm0b5syZg8jISAwaNAifffYZRo0ahYsXL8LT0xNAxW2A3NzccOHCBTz88MM4f/48wxIRNSmZTAYvLy8kJycjKSlJ3+UQNQuWlpbw9PRs8rvJ1GsEC6i4T5CJiaS58Y1KEIQqI1j9+vVD79698emnn2q2+fj4YNy4cYiIiKhyjFGjRuGdd96pcvfxSiUlJSgpKdE8z83NhYeHB0ewiEgnRFGESqXSur8dEdWfXC6HiUnN99U0mBEsAHBzc8Nzzz2HqVOnomvXrjotpjGUlpbi9OnTVe4DFRwcjJiYGABAVlYWLC0tYWZmhsTERFy8eBHt27ev8ZgRERFYvHhxo9ZNRC2XIAgwNTWFqampvkshogaq93hZTk4OVqxYge7du2PgwIFYv3498vPzG6M2nUhPT4dardbcWLKSi4sLUlJSAACXLl1Cnz590LNnTzzyyCNYvXo1HBwcajzmwoULkZOTo3kkJCQ06tdARERExqXeI1jJycn46quvEBUVhePHj+PEiROYM2cOnnjiCUyZMgWDBw9ujDolu3d4UBRFzbaBAwfi/PnzdT6WmZkZzMzMdFofERERNR/1HsFycHDAyy+/jNjYWJw6dQozZ86EQqFAVFQUAgIC4OPjg+XLl+POnTuNUW+9OTk5QS6Xa0arKqWmplYZ1SIiIiLSBUlT6nv37o1PPvkEycnJ2Lx5Mx588EFcu3YNr776Kjw8PDB+/Hj8/PPPel3TRaFQwM/PD3v27NHavmfPHgwcOFBPVREREVFzppPLARUKBSZNmoRJkyYhMTERUVFRiI6Oxo8//oiffvoJrq6uuH37ti66qlZ+fj7i4uI0z2/cuIHY2Fg4ODjA09MT4eHhCAkJQZ8+fTBgwACsW7cO8fHxCAsLa7SaiIiIqOWq9zINdXXlyhVMnz4dR48ehSAIjXq58YEDBxAUFFRle2hoKKKjowFULDS6bNkyJCcnw9fXFytXrsTQoUN10j8XGiUiIjI+jfn5rdOAVVBQgG+++QZffvklYmJiIIoiLC0tMXHiRE3QaY4YsIiIiIyPQa2DVZ3Dhw/jyy+/xHfffYfCwkKIogh/f39MmzYNTz31FGxsbHTRDREREZFRaHDAun37NjZs2IDo6Ghcv34doijC0dER06dPx7Rp0+Dr66vLOomIiIiMRr0D1jfffIOoqCj88ccfUKvVkMlkCA4OxtSpUzFu3DiuPExEREQtXr0D1qRJkwAA7dq1w5QpUzBlyhS0adNG54URERERGasGBaxp06bhwQcfbIx6iIiIiIxevQPW5s2bG6MOIiIiomZD0kruAKBSqbBy5Ur07dsXSqUSJib/ZLbY2FjMmjULV69eldoNERERkdGQtExDUVERgoODERMTAycnJyiVShQUFGj2e3l5ISoqCg4ODnj33XclF0tERERkDCSNYC1duhRHjx5FREQEUlJSMH36dK39tra2CAgIwO7duyUVSURERGRMJAWsbdu2ITAwEAsWLIAgCBAEoUqb9u3bIz4+Xko3REREREZFUsCKj4+Hv79/rW2USiVycnKkdENERERkVCQFLBsbG6SlpdXa5vr163B2dpbSDREREZFRkRSw+vfvj507d9Y4QpWYmIhdu3Zh6NChUrohIiIiMiqSAtb8+fORmZmJ4cOHIyYmBiqVCgBQWFiIvXv3Ijg4GGVlZQgPD9dJsURERETGQNIyDUOHDsWaNWswe/ZsDBkyRLPdxsYGACCXyxEZGQk/Pz9pVRIREREZEUEURVHqQS5duoS1a9fixIkTyMzMhFKpRL9+/TBr1ix069ZNF3UatNzcXNja2iInJwdKpVLf5RAREVEdNObnt04C1v2Ul5dDJpO8aLzBYsAiIiIyPo35+S0p9XzxxRf3baNWq/H0009L6YaIiIjIqEgKWDNnzsSPP/5Y435RFPHMM8/g22+/ldINERERkVGRvEzDU089hcOHD1fZVxmuvvnmG4SFhUnphoiIiMioSApYP//8Mzp06IBHH30U58+f12wXRREhISHYunUrXnjhBaxZs0ZyoURERETGQlLAsrW1xe7du2Fra4uHHnoIN2/ehCiKePbZZ7F582Y8//zz+PTTT3VVKxEREZFRkHxpn5ubG37//XeUlpYiODgYTz31FLZs2YLp06dj7dq1uqiRiIiIyKjoZO2Ezp0749dff0VKSgq+/fZbTJs2DevWrdPFoYmIiIiMTr1Wcl+yZEmt+/v27YvY2Fi4u7trtRUEAW+++WbDKiQiIiIyMvVaaLShi4UKggC1Wt2g1xoDLjRKRERkfBrz87teI1j79+/XaedEREREzVG9AlZAQEBj1UFERETUbDTfGwQSERER6Um9RrBqEx8fDwBwdHSElZWVrg5LREREZHR0NoLVrl07tG/fHl9//bWuDklERERklHQ2guXp6QlBEGBjY6OrQxIREREZJZ0FrJs3b+rqUERERERGjZPciYiIiHRMUsCSy+V45plndFULERERUbMgKWAplUp4eHjoqhYiIiKiZkFSwOrbty/OnTunq1qIiIiImgVJAWvx4sXYt28fNmzYoKt6iIiIiIyepKsIf//9dwQGBmLq1Kn4+OOP0bdvX7i4uEAQBK12giDgzTfflFQoERERkbEQRFEUG/pimaxuA2CCIECtVje0G4PXmHfjJiIiosbRmJ/fkkaw9u/fr6s6iIiIiJoNSQErICBAV3UQERERNRtcaJSIiIhIx3R2qxy1Wo309HSUlJRUu9/T01NXXREREREZNMkB6/Tp03jttddw6NAhlJaWVttGEASoVCqpXREREREZBUkBKzY2FkOGDIGJiQmCg4Oxc+dO9OzZE66urjhz5gzS0tIQGBiItm3b6qpeIiIiIoMnaQ7WO++8AwA4ceIEfvzxRwDA+PHj8euvv+LmzZsICwvDhQsXsGjRIumVEhERERkJSQHryJEjePTRR+Hj46PZVrmsloWFBT755BO4ubnhtddek1YlERERkRGRFLBycnLQvn17zXNTU1Pk5+f/c3CZDIGBgdi7d6+UboiIiIiMiqSA1apVK2RlZWmeu7q64tq1a1ptiouLUVhYKKUbIiIiIqMiaZJ7165dceXKFc3zQYMGYceOHTh+/Dj69++PS5cu4ZtvvoG3t7fkQo1B+LazUCptYWYqh7mpDOamcpib/PN3C1M5zExlsDCVVzxXyKv83cJUDjMTGWQy4f4dEhERkUGSFLAefvhhzJ07F8nJyWjdujVeeeUVbN++HYMGDYKDgwOysrJQXl7eYuZg/X4xFTKz/Ps3rAPzu0HMUmECc1MZLBUmmhBmqagIZJaKiv0WpnJYmclhoTCB1b+2W5nd/VNhAkszOazuHuvem3ETERGRbkm62XNZWRkyMzNhb28PhUIBAIiJicF7772Hv//+G23btsVLL72Ehx9+WGcFG6LKm0Wu3fMnBIUlisvKUVymrvhTpUZxqbriz7vbi8rUKCpVo0RVjqK7+yqfNzZBQEXgUshhbfZP8LI2M4HV3Ye1mRzWZqawMpPDxrxymwlszE3+td0U1mYmkHOkjYiIjFRj3uxZUsBqLvLy8jBs2DCUlZVBrVZj9uzZmDFjRp1fr6sfkLpcRIlKjcLSisBVXHb373cDWWGpGoWlKhSVVf5djaJSlaZ9YakaBaUqFJWqUXB3X0GpGgUlFW0ag5VCDmtzE03gsjE3gdLcFDbm//xdaWEKpYUJbMz++buthSmU5qawVMg5okZERHrBgNXI1Go1SkpKYGlpicLCQvj6+uLkyZNwdHSs0+sb8wekK+XlIorKKgJYYYka+SUqFJSoUFCqQkFJRQir2FbRJr9EhfxilWa75lGsQl6xCqVq3Yy2mciEitBlfjd0WZjC9u7DztIUdhaKiueWprCzMIWdpQL2lhXPzUzkOqmBiIhapsb8/JZ8qxyVSoWPP/4YW7ZsweXLl1FYWKi5LU5sbCzWrVuHOXPmoHPnzpKLbSxyuRyWlpYAKq56VKvVaG65UyYTNKcAYSP9eCUqNfKLK0JXXrEKucVlmvCVW1yGvGIV8u7+Wfk8t6gMuXf/zCkqg6pchKpcRGZBKTILqr/NUm0sFXJN6HKwUsDO0hT2lgrYWyngYGla8aeVAvZ39ztYKWBuylBGRESNT1LAKioqQnBwMGJiYuDk5ASlUomCggLNfi8vL0RFRcHBwQHvvvuu5GJrcujQISxfvhynT59GcnIytm/fjnHjxmm1iYyMxPLly5GcnIxu3bph1apVGDJkiGZ/dnY2AgICcO3aNSxfvhxOTk6NVm9zYGYih5m1HI7WZg16vShWjKjlFqmQU1SG3OIy5BRWBK+cojJkF5Vpglh2YSmyi8qQXfjP83IRmtOkSTnFde7XUiGHg5UCjncDl6O1GRytFXCyqvjT0doMjlYKONuYwcFKAVO5pJVMiIiohZIUsJYuXYqjR4/i/fffx/z587F48WLN7XMAwNbWFgEBAdi9e3ejBqyCggL07NkTU6ZMwYQJE6rs37ZtG+bMmYPIyEgMGjQIn332GUaNGoWLFy/C09MTAGBnZ4dz587hzp07eOyxxzBx4kS4uLg0Ws0tnSAIsFSYwFJhAldb83q9trxcRF6JCtmFpcgqLENWYSmyC0uRWVB2989SZBeWIbOgFFmFpcgoKEVWQSlU5eLdUFaExKyiOvVlb2kKJ2szOFmbwdnmn0fl81Z3H/aWCi6tQUREGpIC1rZt2xAYGIgFCxYAQLWTldu3b4+zZ89K6ea+Ro0ahVGjRtW4f8WKFZg2bRqmT58OAFi1ahV2796NTz/9FBEREVptXVxc0KNHDxw6dAiPP/54tccrKSlBSUmJ5nlubq4OvgqqK5lM0MzTalu3aXIQRRG5xaq7pyNLkJFfEcQyCkqRkV+KjIISZBaUIi2vBBl3T1mqy8W7Aa4M11JrX37DVC7AydoMrZTmcLExg4vSHC7Ku8+V5nC9+1BamHBSPxFRCyApYMXHx2P8+PG1tlEqlcjJyZHSjSSlpaU4ffo0Xn31Va3tlac2AeDOnTuwsLCAUqlEbm4uDh06hJkzZ9Z4zIiICCxevLhR6ybdEoR/QpmXk9V925eXi8gqLEV6finS80uQnl+CtLx/Pe4+T82rCGZlahHJOcVIvs/pSnNTWUXYsjWHm60FXG3N0drOAq2V5mhtZw53OwvYWpgyhBERGTlJAcvGxgZpaWm1trl+/TqcnZ2ldCNJeno61Gp1ldN9Li4uSElJAQAkJiZi2rRpEEURoijixRdfRI8ePWo85sKFCxEeHq55npubCw8Pj8b5AkgvZDLh7vwsM3S5z1UBpapypOdXhK3U3GLcqfwztxgpuSW4k1OMlNxi5BSVobisHDczCnEzo+bbR1mYyuFmZw43Owu42VrA3d4C7nb//Nna1hwmnBtGRGTQJAWs/v37Y+fOncjJyYGtrW2V/YmJidi1a1eVCef6cO+IgCiKmm1+fn6IjY2t87HMzMxgZtawyd3U/ChMZBVhyM6i1nZFpWrcya0Y5ar8Mzmn6J8/s4uRUVCKojI1rqcV4HpaQbXHkcsEuCrN4eFggTb2lmhjbwEPe0t4OlrC08ESztZmnA9GRKRnkgLW/PnzERQUhOHDh2P16tWa5RkKCwtx7NgxvPTSSygrK9Ma7WlqTk5OkMvlmtGqSqmpqZzETk3KQiFHOycrtKvlFGVxmboicGUXITG7CEl3H7ezi3A7qwhJ2cUoVZdXPM8uApBZ5RhmJjJ4OFSELU8HS7R1tEQ7Ryt4OlrCw94SChOOfhERNTZJAWvo0KFYs2YNZs+erbXkgY1NxSkVuVyOyMhI+Pn5SatSAoVCAT8/P+zZs0drvtiePXswduxYvdVFVB1zUzm8nKxqnCdWXi4iPb8ECVmFSMwqQkJmIRIyi5CYXYj4zEIkZRejRFWOuNR8xFUzMV8mAO72FmjnWNFHO0creDlboYOTNdztLXjrIyIiHZG80GhYWBgCAgKwdu1anDhxApmZmVAqlejXrx9mzZqFbt266aLOWuXn5yMuLk7z/MaNG4iNjYWDgwM8PT0RHh6OkJAQ9OnTBwMGDMC6desQHx+PsLCwRq+NSJdkMgGtlOZopTSHX9uq+8vU5UjOLkZ8ZkXgupVRgFsZhbiZUYD4zEIUlqqRkFmEhMwiHL6WrvVahVyGto6WaO9shfbO1ujobI0OrazRwdkKNuamTfQVEhE1D83iVjkHDhxAUFBQle2hoaGIjo4GULHQ6LJly5CcnAxfX1+sXLkSQ4cO1Un/xnCrHCJRFJGWV4Ib6QW4mVGAmxmFuJFWgBvpBbiRUYDSWm427qI0Q8dW1ujUygadXKzR2cUGnVvZwNaSwYuIjBfvRWjgGLDI2JWXi7idXYS/0wvwd1o+rqfl43pqAeLS8pGWV1Lj61rZmKGLqw26uNigs6sNvF1t0KmVDSwUvCURERk+gw1Y48ePx/PPP4+HHnqoRa/bw4BFzVlOURmup+Uj7k4+rqXm4VpqPq7dyb87yb4qQQC8HK3g01oJn9Y28HZVoqubEq1tzVv0vxNEZHgMNmDJZDIIggB3d3dMnToVU6dO1dx6piVhwKKWKL9EhWt38nAlJQ9XKv9MyUNGDTfutrc0RVc3Jbq52aLb3T/bO1lxSQki0huDDViXLl3C559/jk2bNiE9PR0ymQwjRozAjBkz8Oijj8LERPIceqPAgEX0j9S8YlxOzsOl5Ny7jzzEpeVDXV71nxprMxN0dVOih7sturexRXd3W3g5WXGki4iahMEGrEplZWXYsWMH1q9fjz/++AOiKMLJyQmTJ0/GtGnT0LlzZ13UarAYsIhqV1ymxtU7efgrKRd/JeXgr6RcXEzKRUk1E+uV5ibo6WGHBzzs0LONHXp62MHZhgv7EpHuGXzA+reEhASsX78e0dHRSEhIAAAMGTIEM2bMwMSJE5vlCugMWET1p1KXIy4tH+cTc3D+dg7+TMzBxeTcaq9m9HCwQG9Pe/T2tEcvTzv4tFbClLcLIiKJjCpgAYBarcaOHTvw8ssvIykpqaIjQYCDgwNeeeUVhIeHQyZrPv84MmAR6UapqhxX7+QhNiEb5xKyEZuQjbi0fNz7r5SFqRwPeNjBv509/No5oLenHdfqIqJ6M5qAFRcXhy+++AIbN27EnTt3oFAoMGHCBISGhuLs2bNYs2YNEhMTER4ejuXLl+uqW71jwCJqPLnFZTiXkI0zt7JxJj4LZ+OzkFus0mojEwBvVyX6tXdAPy9H9PNygL2VQk8VE5GxMOiAVVJSgm+//RZffPEFDh8+DFEU0aVLF8yYMQOTJ0+Gg4ODVtsRI0bgypUruHPnjuTiDQUDFlHTKS8XcT0tHydvZuHUrUycupmF+MzCKu26uNigf3sHDOjghP7tHWBnycBFRNoa8/Nb0mV+L730Er7++mvk5ORAoVBg0qRJeP755xEQEFBtezMzM4wcORJHjx6V0i0RtWAymYBOLjbo5GKDp/tVLAuTmluMEzcyceJGBk78nYlrqfkVS0fcycOGY7cgCEA3NyUGdXDCwI5O6NvOgYuhElGjkrwOVufOnTWjVY6Ojvd9zdGjR/HHH39g0aJFDe3W4HAEi8iwZOSX4P9uZOLY3xmIuZ5R5cbXChMZ/NvZY0gnZwzp5AQfVyXX4yJqgQz2FOGBAwcQGBiow3KMEwMWkWFLzS3Gsb8zcDQuHUeupSMpp1hrv5O1GQI6OyPI2xlDOjnD1oIT5olaAoMNWFSBAYvIeIiiiL/TC3D4ahoOX0vHsb8zUFiq1uyXywT4edoj0NsZI3xc0LGVNRc+JWqmjCJgxcfHAwAcHR1hZWWli0MaDQYsIuNVqirHqVuZOHAlDfsvp+LaPacT2zpaYriPC4b7uMC/nT1MuP4WUbNhFAFLJpNBJpMhMjISzz//vC4OaTQYsIiaj4TMQhy4koq9l1MRE5eBUvU/C5/aWZpiuI8LHurmisGdnGBuyonyRMbMYK8i/DdPT08IggAbGxtdHZKIqMl5OFgiZEA7hAxoh4ISFQ5fS8Oei6nYd/kOsgrL8N3pRHx3OhFWCjmG+bhglK8rgrq04lWJRKSFc7B0gCNYRM2fulzEyZuZ+O1CCn67kIKU3H8mylsq5Bju44JHerRGQBdnmJkwbBEZA6M4RdiSMWARtSzl5SLOJWbjtwsp+OV8MhKzijT7bMxMENzNFY/1dkf/9o6Qc/kHIoPFgGXgGLCIWi5RFBGbkI2f/0zGL38ma41suSrNMfYBN4zv7Q5vV/7bQGRoDDpglZaWYseOHTh58iSys7OhVqurtBEEAevXr5fSjUFjwCIioGJk69StLOyIvY1f/kxGTlGZZl/X1ko80acNxj7gzvskEhkIgw1Yt27dwogRI3D9+nXUdhhBEKoNXs0FAxYR3atEpcb+y6nYfvY29l1ORZm64t9IhVyGEd1c8GQfDwzq6MRTiER6ZLAB67HHHsOOHTsQEhKCqVOnok2bNjAxqf7CxLZt2za4SEPHgEVEtckuLMWPsUnYdjIBF5NzNdvd7SzwVF8PPOHvgVY25nqskKhlMtiAZWdnB39/f+zZs0eXNRkdBiwiqqsLt3PwzakE7Dh7G7nFKgCAiUzASF9XPNuvLfq3d+DK8URNxGDXwSovL0evXr10VQsRUbPn624LX3dbvDbaB7vOJ2PT8Vs4E5+NX+5Oku/YyhqTB7bDY73dYanQ2VKFRNTEJI1gjRw5EgqFAjt37tRlTUaHI1hEJMXFpFxsOnELO87e1twX0dbCFE/19cRzA9rCzc5CzxUSNU8Ge4rw7NmzGDJkCKKjozFx4kRd1mVUGLCISBfyiitWio86ehPxmYUAKm4+PcrXFWEBHeDrbqvnComaF4MNWEuWLMHJkyexa9cuBAQEoFevXrC1rfoPgCAIePPNNyUVasgYsIhIl9TlIvZdTsWXR27g2N8Zmu1DOjkhLKADBnZw5DwtIh0w2IAlk9XtrvJcpoGIqGEuJuXi88N/46dzSVCXV/xz3d3dFjMDO+Chbq6QcZkHogYz2IB18ODBOrcNCAhoaDcGjwGLiBpbQmYhvjj8N7adSkBxWTkAoLOLNV5+sDNG+TJoETWEwQYsqsCARURNJSO/BNExNxEdcxN5d5d5YNAiahgGLAPHgEVETS2nqAxRR29g/ZEbWkFrXnAXjOjqwjlaRHVg8AErJiYG0dHRiI2N1RTZq1cvPPfccxg8eLAu6jRoDFhEpC/VBS2/tvZ4dZQ3/Ns56Lk6IsNm0AFr3rx5WLlypeZehDKZDOXlFfMDBEHAyy+/jBUrVkiv1IAxYBGRvuUUleGzg9fx5dEbmjlaD3q3woKHvNHF1UbP1REZpsb8/K7bZYA12LhxI1asWIEuXbpgy5YtSE5OhkqlQkpKCrZu3Qpvb2+sXr0aGzdu1FW9RERUDVsLUyx4yBsH5wfhqb6ekMsE7L2cilGrD+H17eeRkV+i7xKJWhRJI1gDBgxAUlISLly4ABubqv9Dys3NRffu3dG6dWscP35cUqGGjCNYRGRorqfl48PdV/DrhRQAgI25CeYM74znBrSFqVzS/62Jmg2DHcG6cOECJkyYUG24AgClUonHHnsMf/31l5RuiIionjo4W+PTZ/2w7fn+6NpaibxiFd75+SJGrjqE/VdS9V0eUbMn+b8x9xsA45UsRET606+9I3a+NBjvP9YdjlYK/J1WgClRJxH21Wmk5BTruzyiZktSwPL19cX333+P/Pz8avfn5eXh+++/R7du3aR0Q0REEshlAib19cT++YGYMcQLcpmA3/5KwfAVBxF19IZmhXgi0h1JASssLAyJiYkYMGAAvv/+e6SnpwMA0tPT8d1332HgwIFITEzEzJkzdVIsERE1nNLcFK8/3BU/vzQYvTztkF+iwuKdFzFuzVGcT8zRd3lEzYrkZRrmzp2L1atXa04F/nuZBlEU8dJLL2H16tXSKzVgnORORMamvFzElpPxeP/Xy8grVkEmAC8EdMCc4Z1gZiLXd3lETcKg18ECgCNHjiAqKgqxsbHIzc3VLDQaGhqKIUOG6KJOg8aARUTGKjWvGO/8fAk7zyUBALq42OCjJ3rC191Wz5URNT6DD1gtHQMWERm73y6kVKyXVVAKuUzAf4I64sWgjlCYcEkHar4MdpmG+Ph45Obm1tomLy8P8fHxUrohIqJG9pCvK36fOxQPd28NdbmI/+29hnFrjiIuNU/fpREZJUkBy8vL677zqyIjI+Hl5SWlGyIiagKO1mZY80xvfPJ0L9hbmuJici7GfHwU35xMuO+SPESkTVLAEkXxvm86vimJiIzLIz3csHvuUAzp5ISiMjUWfP8nZm+NRW5xmb5LIzIajX5yPTExscaV3omIyDC1sjHHhil98cpD3jCRCdh5LgkP/+8wYhOy9V0akVEwqe8LlixZovX8wIED1bZTq9VITEzE1q1b0a9fvwYVR0RE+iOTCZgZ2AH92jtg9pazSMgswsRPY7BoTFc8278t79RBVIt6X0Uok/0z6CUIwn1PAbq5uWH79u3w9/dvWIVGgFcRElFzl1NUhoU//Ild5ytuHv24Xxu8M84X5qZcM4uMl0Et03Dw4EEAFXOrhg0bhsmTJyM0NLRKO7lcDgcHB3h7e2uFsuaIAYuIWgJRFLHu0N/44LfLKBeBnh52WPtsb7S2tdB3aUQNYlAB698WL16MoKAgDB06VJc1GR0GLCJqSQ5fS8NLW84iu7AMTtYKRD7jh75eDvoui6jeDDZgUQUGLCJqaRIyCzFj4ylcTsmDiUxAxGPd8XgfD32XRVQvBrPQ6COPPILTp083qKOioiJ8+OGH+PTTTxv0eiIiMhweDpb4YdZAPNKjNVTlIuZ/9ydW/XGVS/MQ3VWvgJWQkIC+ffviwQcfRHR09H1XcQeAU6dOYc6cOWjbti3eeustODk5NbhYIiIyHJYKE/xvUi/8J6gDAGDVH9cw/7s/UaYu13NlRPpXr1OEoigiKioKS5YsQXx8PGQyGby9vdG7d2+4uLjA3t4eRUVFyMzMxLVr13Dq1Cnk5ORAJpPhiSeewHvvvYd27do14pejHzxFSEQt3eYT8XjzxwtQl4sY0skJkc/0ho25qb7LIqqVwc3BKi8vxy+//ILo6GgcPHgQmZmZVdrIZDL06NED48aNw/Tp0+Hm5qaTghvL+PHjceDAATz44IP47rvv6vVaBiwiImD/5VT8Z/MZFJaq4e1qg41T+6KV0lzfZRHVyOAC1r0uXbqExMREZGRkwMLCAs7OzujWrRtsbW11UWOT2L9/P/Lz87FhwwYGLCKiBjqfmIMp0SeRnl8CLycrfD29H9zsuIwDGabG/Pyu90ru1fHx8YGPj48uDqU3QUFBNa5KT0REddO9jS1+mDkQT31+HDfSC/DEZ8ewZUZ/eDhY6rs0oiZV7xVAo6KikJaW1hi1NNihQ4cwZswYuLm5QRAE7Nixo0qbyMhIeHl5wdzcHH5+fjh8+HDTF0pE1AJ4Olrim7ABaOdoicSsIjzx2TH8nZav77KImlS9A9a0adPg5uaGIUOG4KOPPsK1a9cao656KSgoQM+ePfHJJ59Uu3/btm2YM2cOXn/9dZw9exZDhgzBqFGjEB8f38SVEhG1DO52Ftj2wgB0cLZCck4xnlx3HNfu5Om7LKImU++AFRMTg3nz5iEzMxPz58+Ht7c3fHx88Nprr+H48eONUeN9jRo1Cu+++y4ee+yxavevWLEC06ZNw/Tp0+Hj44NVq1bBw8OjwWtylZSUIDc3V+tBRETaXJTm2PbCAHi72iAtrwRPrjuOqwxZ1ELUO2D1798fERER+Ouvv3Dt2jV88MEHcHZ2xrJlyzBo0CC0bt0aL7zwAnbt2oWSkpLGqLleSktLcfr0aQQHB2ttDw4ORkxMTIOOGRERAVtbW83Dw4OrFxMRVcfJ2gxbn++P7u62yCwoxbNfnEB8RqG+yyJqdJLuwtyhQwfMmzcPhw4dwp07d/DFF1+gX79++PrrrzFmzBg4OTlhwoQJ+Oqrr6pdyqEppKenQ61Ww8XFRWu7i4sLUlJSNM9HjhyJxx9/HLt27UKbNm1w8uTJGo+5cOFC5OTkaB4JCQmNVj8RkbGzs1Tgq2l90cXFBql5JXh2/Qmk5hbruyyiRiUpYP2bo6MjpkyZgh07diA9PR3bt2/HE088gaNHjyI0NBQuLi4IDAzUVXf1JgiC1nNRFLW27d69G2lpaSgsLERiYiL8/f1rPJaZmRmUSqXWg4iIalYZsjwdLBGfWYiQ9f+H7MJSfZdF1Gh0FrD+zdzcHI8++ijWr1+P5ORkHDlyBHPnztUaMWoqTk5OkMvlVfpOTU2tMqpFRESNp5XSHF9P74dWNma4cicPk6NOoqBEpe+yiBpFowQsURRx7do1JCYmQhAEDBw4EMuWLcPly5cbo7taKRQK+Pn5Yc+ePVrb9+zZg4EDBzZ5PURELZmHgyU2Te8HO0tTxCZkY8bGUyhRqfVdFpHOSQpYP/74I6ZOnYqsrCzNtps3b6J79+7w9vZG27Zt8cwzz6C8vHFv/Jmfn4/Y2FjExsYCAG7cuIHY2FjNMgzh4eH44osv8OWXX+LSpUuYO3cu4uPjERYW1qh1ERFRVZ1dbBA9pS+sFHLEXM/Awh/OQwc3FSEyKJJWcl+7di0SExNhb2+v2TZnzhxcvHgRw4YNQ0ZGBrZu3Yphw4Zh2rRpkoutyalTpxAUFKR5Hh4eDgAIDQ1FdHQ0nnzySWRkZGDJkiVITk6Gr68vdu3ahbZt2zZaTUREVLMHPOzw6bN+mBJ9Ej+cuY1OrWwwM7CDvssi0hlJ9yL09PTEiBEjsH79egBATk4OnJ2dMWHCBGzZsgVlZWXo1asXlEplg5dEMAa8FyERUcNsPHYTb/34FwQBWPusH0Z2c9V3SdSCNObnt6RThGlpaWjdurXm+ZEjR6BSqfDUU08BAExNTTFixAjExcVJq5KIiJql5wa0w3MD2kIUgbnbYvFXUo6+SyLSCUkBS6lUIiMjQ/P8wIEDkMlkGDJkiGabqakpCgoKpHRDRETN2FuPdMXgjk4oLFVjxoZTSM3jGllk/CQFLG9vb+zcuROZmZnIycnB1q1b0bt3b605Wbdu3eJyCEREVCMTuQxrnu6N9s5WSMopxvMbT6O4jFcWknGTFLBmz56NpKQkuLu7w8PDA0lJSVpX5qnVahw5cgQ9e/aUXCgRETVftpamWB/qD1uLiuUb3v3lor5LIpJEUsCaMGEC1qxZg27duqFz586IiIjA1KlTNfv37t2LwsJCPPTQQ5ILJSKi5s3LyQofP9ULALDpeDx+/jNJzxURNZykqwipAq8iJCLSnWW/XUbkgeuwMTPBz7MHo62jlb5LombKYK8ijI+PR25ubq1t8vLyNAt+EhER3U/4iM7o09YeeSUqvLj5LFd6J6MkKWB5eXlh9erVtbaJjIyEl5eXlG6IiKgFMZHL8L+nesHO0hTnb+fg/V+b/jZrRFJJCliiKN739gY8A0lERPXlZmeBDydWXCAVdfQmdv+VoueKiOqnUW72/G+JiYmwsbFp7G6IiKiZGd7VBdMHV5wBmf/tOdzOLtJzRUR1V+97ES5ZskTr+YEDB6ptp1arkZiYiK1bt6Jfv34NKo6IiFq2BQ954+TNTJxLzMHCH85jwxR/CIKg77KI7qveVxHKZP8MegmCcN9TgG5ubti+fTv8/f0bVqER4FWERESNJy41H6P/dxilqnIsn9gDj/fx0HdJ1Ew05ud3vUew9u/fD6BibtWwYcMwefJkhIaGVmknl8vh4OAAb29vrVBGRERUHx1bWWPu8M744LfLeOfniwjo7IxWSnN9l0VUK0nrYC1evBhBQUEYOnSoLmsyOhzBIiJqXCp1OcZHxuD87RyM6OqCdSF+PFVIkjXm5zcXGtUBBiwiosZ3KTkXYz4+AlW5iI+f6oUxPd30XRIZOYM6RXiv0tJS7NixAydPnkR2djbU6qoLwgmCgPXr10vtioiIWjCf1kr8J6gjVu+9hkU//YWBHRzhaG2m77KIqiVpBOvWrVsYMWIErl+/Xutkd0EQqg1ezQVHsIiImkapqhxjPj6CK3fy8GhPN/zv7r0LiRrCYEew5s6di7i4OISEhGDq1Klo06YNTEwkD4oRERFVS2Eiw/LHe2DcmqP46VwSxvdyR5B3K32XRVSFpBEsOzs7+Pv7Y8+ePbqsyehwBIuIqGm998tFfH74Bto7WeG3OUOhMOHV6lR/Bnuz5/LycvTqxeFZIiJqWrMf7AQnawX+Ti/AxmM39V0OURWSAtaAAQNw6dIlXdVCRERUJzbmppgX3AUAsHrvNWTkl+i5IiJtkgLW+++/j/379+O7777TVT1ERER18ngfD3RtrUResQor/7iq73KItEiakb5z504EBQXhySefREBAAHr16gVbW9sq7QRBwJtvvimlKyIiIi1ymYC3xnTFpHXHsflEPJ7t3xberpwHS4ZB0iT3ut4Ch8s0EBFRY5m56TR+vZCCQR0dsWlaP67wTnVmsMs0VN6XkIiISF9eG+2DvZdTcTQuA3su3kFwN1d9l0QkLWAFBAToqg4iIqIG8XCwxPTBXog8cB3v7bqEgC7OMDOR67ssauG4cAgRERm9WUEd4WxjhlsZhdh8Il7f5RBJD1gqlQorV65E3759oVQqtVZyj42NxaxZs3D1Kq/uICKixmNtZoI5wzsBANYevI7isuY775eMg6SAVVRUhKCgIMybNw+3bt2CUqnUuiehl5cXoqKisHHjRsmFEhER1eZxPw+42ZrjTm4Jtp1M0Hc51MJJClhLly7F0aNHERERgZSUFEyfPl1rv62tLQICArB7925JRRIREd2PwkSGWUEdAQCRB+I4ikV6JSlgbdu2DYGBgViwYAEEQaj20tj27dsjPp7nw4mIqPE93qcNWt8dxfrmFEexSH8kBaz4+Hj4+/vX2kapVCInJ0dKN0RERHViZiL/ZxRrP+dikf5IClg2NjZIS0urtc3169fh7OwspRsiIqI6e+LuKFZKbjFHsUhvJAWs/v37Y+fOnTWOUCUmJmLXrl0YOnSolG6IiIjqzMxEjlmBHQBUjGKVqDiKRU1PUsCaP38+MjMzMXz4cMTExEClUgEACgsLsXfvXgQHB6OsrAzh4eE6KZaIiKgunvD3gKvy7igWrygkPZB0L0IAWLt2LWbPnl3tvQblcjkiIyOrXF3Y3PBehEREhmfjsZt468e/0NrWHAfmB3J1d6qiMT+/JS80GhYWhnPnzuHFF1+Ev78/OnTogF69eiEsLAxnz55t9uGKiIgM0xN9KkaxknOK8cOZ2/ouh1oYySNYxBEsIiJD9cXhv/HuL5fQxcUGv80ZUu1yQtRyGewI1ogRI/DVV1+hoKBAV/UQERHpzON9PGCpkOPKnTwcu56h73KoBZEUsPbt24fJkyfDxcUFzz77LHbv3o3y8nJd1UZERCSJrYUpJvRuAwCIirmp32KoRZG80OjSpUvRvn17bN68GaNHj4a7uzv++9//4syZM7qqkYiIqMFCB7YDAPxx6Q7iMwr1Wwy1GJIClru7O1555RX8+eefiI2Nxdy5cyGXy7Fy5Ur4+/ujW7du+OCDD5CQwEtkiYhIPzq2ssbQzs4QRWDDsZv6LodaCJ1PchdFEXv37sWmTZuwfft25OfnQyaToaysTJfdGBROciciMmz7L6diSvRJ2JiZ4PhrD8LKzETfJZEBMNhJ7tURBAHDhw/H4sWL8dJLL8HExITzsoiISK8COjvDy8kKeSUqfH8mUd/lUAug04CVnZ2NdevWYejQoejQoQMiIiJgYWGByZMn67IbIiKiepHJBIQOaAsAiI65ifJyrlBEjUvyGGlpaSl27tyJTZs24ddff0VZWRnkcjkeeughhISEYOzYsTA3N9dFrURERA02sY8HPvz9Kv5OK8Cha2kI7NJK3yVRMyYpYM2YMQPff/89cnJyIIoi+vTpg5CQEEyaNAnOzs66qpGIiEgyazMTPN6nDaKO3kR0zE0GLGpUkgLW+vXr0a5dO/znP/9BSEgIOnfurKu6iIiIdC50QDtEx9zEgStpuJ6Wjw7O1vouiZopSQHr4MGDGDJkiK5qISIialTtnKwQ1KUV9l1OxTcnE7BwtI++S6JmStIkd4YrIiIyNk/6ewAAfjh7Gyo1r3KnxiH5KkKVSoWVK1eib9++UCqVMDH5Z1AsNjYWs2bNwtWrV6V2Q0REpBNBXVrBwUqBtLwSHLqWpu9yqJmSFLCKiooQFBSEefPm4datW1Aqlfj3uqVeXl6IiorCxo0bJRdKRESkCwoTGcY94A4A+O4018SixiEpYC1duhRHjx5FREQEUlJSMH36dK39tra2CAgIwO7duyUVSUREpEsT/SpuAP3HxVRkFZTquRpqjiQFrG3btiEwMBALFiyAIAgQBKFKm/bt2yM+Pl5KN01i/PjxsLe3x8SJE/VdChERNbKubkp0ba1EqbocO/9M0nc51AxJCljx8fHw9/evtY1SqUROTo6UbprE7NmzeSqTiKgFqRzF4mlCagySApaNjQ3S0mqfIHj9+nWjWHQ0KCgINjY2+i6DiIiayNgH3GAiE/BnYg6upOTpuxxqZiQFrP79+2Pnzp01jlAlJiZi165dGDp0qJRucOjQIYwZMwZubm4QBAE7duyo0iYyMhJeXl4wNzeHn58fDh8+LKlPIiJq3hytzTDMu2I1d94AmnRNUsCaP38+MjMzMXz4cMTExEClUgEACgsLsXfvXgQHB6OsrAzh4eGSiiwoKEDPnj3xySefVLt/27ZtmDNnDl5//XWcPXsWQ4YMwahRo7Tmfvn5+cHX17fKIymJ596JiFqqytOEP5zhmlikW5JWch86dCjWrFmD2bNnay06WnmqTS6XIzIyEn5+fpKKHDVqFEaNGlXj/hUrVmDatGmaqxhXrVqF3bt349NPP0VERAQA4PTp05Jq+LeSkhKUlJRonufm5urs2ERE1HSCvFvB0UqB9PwSHLyahgd9XPRdEjUTkhcaDQsLw7lz5/Diiy/C398fHTp0QK9evRAWFoazZ89WWbpB10pLS3H69GkEBwdrbQ8ODkZMTEyj9BkREQFbW1vNw8PDo1H6ISKixmUql2FcL66JRbonaQSrko+PD1avXq2LQ9Vbeno61Go1XFy0/9fh4uKClJSUOh9n5MiROHPmDAoKCtCmTRts3769xiskFy5cqHXaMzc3lyGLiMhITfRrg/VHbuCPS3eQVVAKeyuFvkuiZkAnAcsQ3LsGlyiK1a7LVZP6LIZqZmYGMzOzOrcnIiLD5dNaiW5uSvyVlIuf/0xCyIB2+i6JmgHJpwj1zcnJCXK5vMpoVWpqapVRLSIiouqMfcANALDrfN3PfBDVpl4B6+mnn8YPP/zQ4M6kvr46CoUCfn5+2LNnj9b2PXv2YODAgTrti4iImqdRvq0BACduZCA9v+Q+rYnur14Ba+vWrbhw4UKDO2vo6/Pz8xEbG4vY2FgAwI0bNxAbG6tZhiE8PBxffPEFvvzyS1y6dAlz585FfHw8wsLCGlwrERG1HB4OlujRxhblIrD7L45ikXT1noMVGxvb5LeUOXXqFIKCgjTPKyeYh4aGIjo6Gk8++SQyMjKwZMkSJCcnw9fXF7t27ULbtm2btE4iIjJeo7u3xp+JOdh1PhnP9OPnB0kjiKIo1rWxTFYx4FWfyeOVKrt5++238dZbb9X79YYsNzcXtra2yMnJgVKp1Hc5RETUAPEZhRi6fD9kAnDy9eFwtObFTM1dY35+12sEKyoqSnKHDzzwgORjEBER6ZqnoyV83ZW4cDsXv1+8g6f6euq7JDJi9QpYoaGhjVUHERGR3o3u3hoXbudi1/lkBiySxOiXaSAiItKV0XevJoy5noHMglI9V0PGjAGLiIjornZOVujaWgl1uYg9F3k1ITUcAxYREdG/PNyjYhTrFy46ShIwYBEREf3LKF9XAEBMXDqyC3makBqGAYuIiOhf2jtbw9vVBqpyEb9fvKPvcshIMWARERHd4+HuFacJd51P1nMlZKwYsIiIiO4x6m7AOhqXjpzCMj1XQ8ZIcsBSqVRYuXIl+vbtC6VSCROTf5bWio2NxaxZs3D16lWp3RARETWZjq2s0cXFBmVqEX9c4mlCqj9JAauoqAhBQUGYN28ebt26BaVSiX/fecfLywtRUVFNfu9CIiIiqYK7uQAA9l1O1XMlZIwkBaylS5fi6NGjiIiIQEpKCqZPn66139bWFgEBAdi9e7ekIomIiJpakHcrAMChq2koU5fruRoyNpIC1rZt2xAYGIgFCxZAEIRqbwLdvn17xMfHS+mGiIioyfVsYwcHKwXySlQ4fStL3+WQkZEUsOLj4+Hv719rG6VSiZycHCndEBERNTm5TEBgZ2cAwH6eJqR6khSwbGxskJaWVmub69evw9nZWUo3REREelF5mpDzsKi+JAWs/v37Y+fOnTWOUCUmJmLXrl0YOnSolG6IiIj0YmhnZ8hlAq6l5iMhs1Df5ZARkRSw5s+fj8zMTAwfPhwxMTFQqVQAgMLCQuzduxfBwcEoKytDeHi4ToolIiJqSrYWpvBraw8A2H+Fo1hUdyb3b1KzoUOHYs2aNZg9ezaGDBmi2W5jYwMAkMvliIyMhJ+fn7QqiYiI9GSYdyv8341M7LuciucGtNN3OWQkJAUsAAgLC0NAQADWrl2LEydOIDMzE0qlEv369cOsWbPQrVs3XdRJRESkF8O8W+H9Xy/j2PUMFJWqYaGQ67skMgKSAxYA+Pj4YPXq1bo4FBERkUHp1Moa7nYWuJ1dhJjr6XjQx0XfJZERkLxMQ25ubq1t8vLyuA4WEREZLUEQMIxXE1I9SQpYXl5e9x25ioyMhJeXl5RuiIiI9KoyYO2/nKp1SziimkgKWKIo3vcXjb+IRERk7Pq3d4SZiQxJOcW4cidP3+WQEZAUsOoiMTFRc1UhERGRMbJQyDGwgyMAniakuqn3JPclS5ZoPT9w4EC17dRqNRITE7F161b069evQcUREREZimHerbD/ShoOXE7DrMCO+i6HDJwg1vMcnkz2z6CXIAj3PQXo5uaG7du33/eehcYsNzcXtra2yMnJgVKp1Hc5RETUCBKzCjH4g/2QywSceWMEbC1N9V0SSdSYn9/1HsHav38/gIq5VcOGDcPkyZMRGhpapZ1cLoeDgwO8vb21QhkREZExamNvic4u1rh6Jx8Hr6Xh0Z5u+i6JDFi9A1ZAQIDm74sWLUJQUBDvNUhERC1CQGdnXL2Tj5i4dAYsqpWkoaVFixYxXBERUYsxsIMTAODY3xl6roQMnU5WcgcqJrWnp6ejpKSk2v2enp666oqIiEgv/L0cIJcJuJVRiNvZRXC3s9B3SWSgJAes06dP47XXXsOhQ4dQWlpabRtBEKBSqaR2RUREpFfWZibo0cYWZ+Ozcex6Bib6tdF3SWSgJJ0ijI2NxZAhQ3Ds2DEEBwdDFEX06NEDwcHBcHJygiiKCAgIQEhIiK7qJSIi0qsB7SvWw4q5nq7nSsiQSQpY77zzDgDgxIkT+PHHHwEA48ePx6+//oqbN28iLCwMFy5cwKJFi6RXSkREZAAq52Edv57Bu5VQjSQFrCNHjuDRRx+Fj4+PZlvlL5uFhQU++eQTuLm54bXXXpNWJRERkYHwa2sPU7mApJxi3Moo1Hc5ZKAkBaycnBy0b99e89zU1BT5+fn/HFwmQ2BgIPbu3SulGyIiIoNhoZCjl6c9AF5NSDWTFLBatWqFrKwszXNXV1dcu3ZNq01xcTEKC5nwiYio+fhnHhYDFlVPUsDq2rUrrly5onk+aNAg/P777zh+/DgA4NKlS/jmm2/g7e0trUoiIiIDUnnj52Och0U1kBSwHn74YRw6dAjJyckAgFdeeQWiKGLQoEFwdnZG9+7dkZ2dzTlYRETUrDzgaQczExnS80sQl5p//xdQiyMpYIWFheH27dtwdKxI8j179sTevXvx0EMPwcnJCcOHD8fOnTsxfvx4nRRLRERkCMxM5PBv5wCA87CoepIWGjU1NYWLi4vWtoEDB+KXX36RVBQREZGhG9DBEUfi0hETl4HnBrTTdzlkYCSNYMnlcjzzzDO6qoWIiMhoDLg7D+v4jQyUl3MeFmmTFLCUSiU8PDx0VQsREZHR6O5uCyuFHNmFZbiUkqvvcsjASApYffv2xblz53RVCxERkdEwlcvg73V3HhaXa6B7SApYixcvxr59+7BhwwZd1UNERGQ0/r1cA9G/SZrk/vvvvyMwMBBTp07Fxx9/jL59+8LFxQWCIGi1EwQBb775pqRCiYiIDM2A9hX3Jfy/G5lQqcthIpc0bkHNiCBKWCFNJqvbL5IgCFCr1Q3txuDl5ubC1tYWOTk5UCqV+i6HiIiaiLpcRK8lvyO3WIUd/xmEBzzs9F0S1UNjfn5LGsHav3+/ruogIiIyOnKZgH7tHbHn4h2c+DuDAYs0JAWsgIAAXdVBRERklPza2mPPxTuITcjWdylkQHiymIiISILKUauz8dl6rYMMCwMWERGRBD3a2EImACm5xUjJKdZ3OWQgGLCIiIgksFSYoItrxQTp2IQsPVdDhoIBi4iISCKeJqR7MWARERFJ1MvTDgBwlhPd6S5JASs+Ph65ubXffykvLw/x8fFSuiEiIjJove6OYJ1PzIFKXa7fYsggSApYXl5eWL16da1tIiMj4eXlJaWbRpeXlwd/f3888MAD6N69Oz7//HN9l0REREakg7M1bMxMUFSmxpU7efouhwyApHWwRFHE/RaCl7BQfJOxtLTEwYMHYWlpicLCQvj6+uKxxx6Do6OjvksjIiIjIJMJ6OlhhyNx6YhNyEY3N1t9l0R61uhzsBITE2FjY9PY3Ugil8thaWkJACguLoZarTaKYEhERIaDE93p3+o9grVkyRKt5wcOHKi2nVqtRmJiIrZu3Yp+/fo1qLhKhw4dwvLly3H69GkkJydj+/btGDdunFabyMhILF++HMnJyejWrRtWrVqFIUOG1LmP7OxsBAQE4Nq1a1i+fDmcnJwk1UxERC1L5UR3ruhOQAMC1ttvv635uyAIOHDgQI0hCwDc3NzwwQcfNKQ2jYKCAvTs2RNTpkzBhAkTquzftm0b5syZg8jISAwaNAifffYZRo0ahYsXL8LT0xMA4Ofnh5KSkiqv/f333+Hm5gY7OzucO3cOd+7cwWOPPYaJEyfCxcVFUt1ERNRyVI5gxaXmI6eoDLYWpvotiPRKEOt5LuzgwYMAKuZWDRs2DJMnT0ZoaGiVdnK5HA4ODvD29oZMprszkYIgVBnB6tevH3r37o1PP/1Us83Hxwfjxo1DREREvfuYOXMmhg0bhscff7za/SUlJVphLTc3Fx4eHo1yN24iIjIeQ5btQ0JmEb6a1hdDOjnruxy6j9zcXNja2jbK53e9R7D+fYPnRYsWISgoCEOHDtVpUfVRWlqK06dP49VXX9XaHhwcjJiYmDod486dO7CwsIBSqURubi4OHTqEmTNn1tg+IiICixcvllQ3ERE1P7087JGQWYTY+GwGrBZO0tDSokWL9BquACA9PR1qtbrK6TwXFxekpKTU6RiJiYkYOnQoevbsicGDB+PFF19Ejx49amy/cOFC5OTkaB4JCQmSvgYiImoeKk8Tch4WSVqmwZAIgqD1XBTFKttq4ufnh9jY2Dr3ZWZmBjMzs/qUR0RELcAD/1rRvT6fQ9T8SA5YpaWl2LFjB06ePIns7Gyo1eoqbQRBwPr166V2VS0nJyfI5fIqo1WpqamcpE5ERE2qm5sSCrkMmQWlSMgsgqejpb5LIj2RFLBu3bqFESNG4Pr167WuG9WYAUuhUMDPzw979uzB+PHjNdv37NmDsWPHNkqfRERE1TEzkcPHTYlzCdk4m5DFgNWCSQpYc+fORVxcHEJCQjB16lS0adMGJia6P+uYn5+PuLg4zfMbN24gNjYWDg4O8PT0RHh4OEJCQtCnTx8MGDAA69atQ3x8PMLCwnReCxERUW16edhVBKz4bIx9wF3f5ZCeSEpD+/btw4MPPogNGzboqp5qnTp1CkFBQZrn4eHhAIDQ0FBER0fjySefREZGBpYsWYLk5GT4+vpi165daNu2baPWRUREdK9ennaIjuFE95ZOUsAqLy9Hr169dFVLjQIDA+9765pZs2Zh1qxZjV4LERFRbSqvJLyYlIsSlRpmJnL9FkR6IWmZhgEDBuDSpUu6qoWIiMjoeTpYwsFKgVJ1OS4m5eq7HNITSQHr/fffx/79+/Hdd9/pqh4iIiKjJgiCZhTrHE8TtliSThHu3LkTQUFBePLJJxEQEIBevXrB1ta2SjtBEPDmm29K6YqIiMhodG2txL7LqbickqfvUkhP6n0vwn+r6z0GBUGodn2s5qIx72VERETGZ+e5JLy05Sx6edph+6xB+i6HamBQ9yL8t/379+uqDiIiomaji6sNAOBqSh7Ky0XIZFzRvaWRFLD+feNnIiIiquDlZAVTuYCCUjVuZxfBw4ELjrY0kia5ExERUVWmchk6OFsDAK5wHlaLJDlgqVQqrFy5En379oVSqdRayT02NhazZs3C1atXpXZDRERkVCpPE165w4DVEkk6RVhUVITg4GDExMTAyckJSqUSBQUFmv1eXl6IioqCg4MD3n33XcnFEhERGQtNwOIIVoskaQRr6dKlOHr0KCIiIpCSkoLp06dr7be1tUVAQAB2794tqUgiIiJj482A1aJJCljbtm1DYGAgFixYAEEQIAhVr5Jo37494uPjpXRDRERkdDq7VASs62n5KFWV67kaamqSAlZ8fDz8/f1rbaNUKpGTkyOlGyIiIqPjbmcBazMTqMpF3EgvuP8LqFmRFLBsbGyQlpZWa5vr16/D2dlZSjdERERGRxAEdHapuJLwcgrvSdjSSApY/fv3x86dO2scoUpMTMSuXbswdOhQKd0QEREZpS6uFauDcx5WyyMpYM2fPx+ZmZkYPnw4YmJioFKpAACFhYXYu3cvgoODUVZWhvDwcJ0US0REZEwqJ7pf5VINLY6kZRqGDh2KNWvWYPbs2RgyZIhmu41NxS+UXC5HZGQk/Pz8pFVJRERkhConuvOmzy2PpIAFAGFhYQgICMDatWtx4sQJZGZmQqlUol+/fpg1axa6deumizqJiIiMTuUIVmJWEfJLVLA2k/yxS0ZCJz9pHx8frF69WheHIiIiajbsrRRoZWOG1LwSXL2Th96e9vouiZqIpDlYI0aMwFdffaW1ejsRERH9gyu6t0ySAta+ffswefJkuLi44Nlnn8Xu3btRXs7F1IiIiCp1cWHAaokkLzS6dOlStG/fHps3b8bo0aPh7u6O//73vzhz5oyuaiQiIjJaHMFqmSQFLHd3d7zyyiv4888/ERsbi7lz50Iul2PlypXw9/dHt27d8MEHHyAhIUFX9RIRERkVTcC6kwdRFPVcDTUVQdTxT1sURezduxebNm3C9u3bkZ+fD5lMhrKyMl12Y1Byc3Nha2uLnJwcKJVKfZdDREQGpKhUja6LfoMoAidfHw5nGzN9l0R3Nebnt86vFxUEAcOHD0enTp3Qpk0bLF++XLMAKRERUUtjoZCjrYMlbmYU4kpKHgNWC6HTgJWdnY1vvvkGmzZtQkxMDERRhI2NDSZMmKDLboiIiIxKF1cb3MwoxOWUXAzu5KTvcqgJSA5YpaWl2LlzJzZt2oRff/0VZWVlkMvleOihhxASEoKxY8fC3NxcF7USEREZpS6uSuz+6w5vmdOCSApYM2bMwPfff4+cnByIoog+ffogJCQEkyZNgrOzs65qJCIiMmpcqqHlkRSw1q9fj3bt2uE///kPQkJC0LlzZ13VRURE1Gx00dz0OR/l5SJkMkHPFVFjkxSwDh48qHWTZyIiIqqqnaMlFCYyFJWpkZBViLaOVvouiRqZpHWwGK6IiIjuz0QuQ0dnawDAZZ4mbBHqHbB69+6NdevWaW3bvXs3wsPDq22/ePFimJjw7uFERNSyVZ4mjEvN13Ml1BTqHbBiY2ORkpKite348eNYvXp1ja/hyrVERNTSeTpYAgASMgv1XAk1BUmnCImIiKhuPCoDVhYDVkvAgEVERNQEPOwtAAAJmUV6roSaAgMWERFRE6gcwUrKLoK6nFNnmjsGLCIioibgojSHqVyAqlxEcg5HsZo7BiwiIqImIJcJcLfjacKWokHrJ2zatAnHjx/XPI+LiwMAjB49ukrbyn1EREQtnYeDJW5mFCIhsxADOjjquxxqRA0KWHFxcdUGp99++63a9oLAWwIQERHxSsKWo94B68aNG41RBxERUbPnYc+1sFqKegestm3bNkYdREREzZ6Hw905WFmcg9XccZI7ERFRE+EIVsvBgEVERNREKudgpeaVoLhMredqqDExYBERETURe0tTWCnkAIBEniZs1hiwiIiImoggCLySsIVgwCIiImpCmoDFeVjNGgMWERFRE+JE95aBAYuIiKgJaZZq4O1ymjUGLCIioiakGcHiHKxmjQGLiIioCXEOVsvAgEVERNSE2thXnCLMLVYhp6hMz9VQY2HAIiIiakJWZiZwtFIA4ChWc8aARURE1MR4mrD5Y8C6y8TEBA888AAeeOABTJ8+Xd/lEBFRM8bFRps/E30XYCjs7OwQGxur7zKIiKgF8LDnUg3NHUewiIiImhhHsJo/owhYhw4dwpgxY+Dm5gZBELBjx44qbSIjI+Hl5QVzc3P4+fnh8OHD9eojNzcXfn5+GDx4MA4ePKijyomIiKriau7Nn1GcIiwoKEDPnj0xZcoUTJgwocr+bdu2Yc6cOYiMjMSgQYPw2WefYdSoUbh48SI8PT0BAH5+figpKany2t9//x1ubm64efMm3NzccOHCBTz88MM4f/48lEplo39tRETU8lSu5p6YVQRRFCEIgp4rIl0TRFEU9V1EfQiCgO3bt2PcuHGabf369UPv3r3x6aefarb5+Phg3LhxiIiIqHcfo0aNwjvvvIM+ffpUu7+kpEQrrOXk5MDT0xMJCQkMZUREdF9l6nL4vbMH5SKw/78BcFaa67ukFik3NxceHh7Izs6Gra2tTo9tFCNYtSktLcXp06fx6quvam0PDg5GTExMnY6RlZUFS0tLmJmZITExERcvXkT79u1rbB8REYHFixdX2e7h4VG/4omIqMXruErfFVBGRgYD1r3S09OhVqvh4uKitd3FxQUpKSl1OsalS5fwwgsvQCaTQRAErF69Gg4ODjW2X7hwIcLDwzXPs7Oz0bZtW8THx+v8B9Tc+Pv74+TJk/ouo870WW9j963r40s9npTX1/e19Wlf+T9cjlDfH9/fhtM33991a195Bqq2z/yGMvqAVene89f1Oac9cOBAnD9/vs59mZmZwczMrMp2W1tb/gN8H3K53Ki+R/qst7H71vXxpR5Pyuvr+9qG9KVUKo3qd1cf+P42nL75/q5fXzKZ7q/5M4qrCGvj5OQEuVxeZbQqNTW1yqgW6d9//vMffZdQL/qst7H71vXxpR5Pyuvr+1pj+z00Fsb2feX7u+mO1xLf381mkrufnx8iIyM127p27YqxY8c2aJJ7feXm5sLW1hY5OTlG9b83Iro/vr+Jmq/GfH8bxSnC/Px8xMXFaZ7fuHEDsbGxcHBwgKenJ8LDwxESEoI+ffpgwIABWLduHeLj4xEWFtYk9ZmZmWHRokXVnjYkIuPG9zdR89WY72+jGME6cOAAgoKCqmwPDQ1FdHQ0gIqFRpctW4bk5GT4+vpi5cqVGDp0aBNXSkRERGQkAYuIiIjImBj9JHciIiIiQ8OARURERKRjDFhEREREOsaA1cTGjx8Pe3t7TJw4Ud+lEJFEP//8M7p06YJOnTrhiy++0Hc5RKRDUj+vOcm9ie3fvx/5+fnYsGEDvvvuO32XQ0QNpFKp0LVrV+zfvx9KpRK9e/fGiRMnGuWWG0TU9KR+XnMEq4kFBQXBxsZG32UQkUT/93//h27dusHd3R02NjYYPXo0du/ere+yiEhHpH5eM2D9y6FDhzBmzBi4ublBEATs2LGjSpvIyEh4eXnB3Nwcfn5+OHz4cNMXSkSSSX2/JyUlwd3dXfO8TZs2uH37dlOUTkT3YQif5wxY/1JQUICePXvik08+qXb/tm3bMGfOHLz++us4e/YshgwZglGjRiE+Pl7Txs/PD76+vlUeSUlJTfVlEFEdSH2/Vze7oq43mCeixqWLz3PJRKoWAHH79u1a2/r27SuGhYVpbfP29hZfffXVeh17//794oQJE6SWSEQ60pD3+9GjR8Vx48Zp9s2ePVv8+uuvG71WIqofKZ/nUj6vOYJVR6WlpTh9+jSCg4O1tgcHByMmJkZPVRFRY6jL+71v3764cOECbt++jby8POzatQsjR47UR7lEVA9N9XluFDd7NgTp6elQq9VwcXHR2u7i4oKUlJQ6H2fkyJE4c+YMCgoK0KZNG2zfvh3+/v66LpeIJKjL+93ExAQfffQRgoKCUF5ejgULFsDR0VEf5RJRPdT181zq5zUDVj3dO8dCFMV6zbvgVUZExuN+7/dHH30Ujz76aFOXRUQ6cL/3t9TPa54irCMnJyfI5fIqo1WpqalVUjARGTe+34mar6Z6fzNg1ZFCoYCfnx/27NmjtX3Pnj0YOHCgnqoiosbA9ztR89VU72+eIvyX/Px8xMXFaZ7fuHEDsbGxcHBwgKenJ8LDwxESEoI+ffpgwIABWLduHeLj4xEWFqbHqomoIfh+J2q+DOL93aBrD5up/fv3iwCqPEJDQzVt1qxZI7Zt21ZUKBRi7969xYMHD+qvYCJqML7fiZovQ3h/816ERERERDrGOVhEREREOsaARURERKRjDFhEREREOsaARURERKRjDFhEREREOsaARURERKRjDFhEREREOsaARURERKRjDFhEREREOsaARURERKRjDFhE1CJFR0dDEATNY9KkSVr7Dxw4AEEQ8Pbbb+unwHvExcVp1duuXTt9l0REtWDAIiKD0qNHD60gUd0jLi5OZ/2NHTsWixYtwsSJE3V2zNo89dRTEAQBW7durbVdRkYGzMzM4OTkhNLSUjg4OGDRokVYtGgRbG1tm6RWImo4E30XQERUqbi4GJcuXUKrVq0wc+bMatvI5XJ06NBBZ32OGzcOkydP1tnx7mfatGnYunUroqKiqoya/dumTZtQWlqKkJAQKBQKODg4aEbToqOjm6ZYImowBiwiMhjnzp2DSqXC0KFDDebUnK49+OCDaNeuHf744w8kJCTAw8Oj2nZRUVEAKgIZERkfniIkIoNx5swZAECfPn30XEnNsrOzMWTIEMjlcqxbt05r36FDhzBmzBg4OTnBzMwMnTp1whtvvIHCwkJNG0EQMGXKFJSXl2PDhg3V9nH69GmcO3cOffv2ha+vb6N+PUTUOBiwiMhgnD59GoDhBqykpCQMGTIEJ0+exLfffovnn39es2/t2rUIDAxETEwMHnnkEcyePRvu7u547733MGLECJSWlmraTpkyBTKZDNHR0RBFsUo/HL0iMn48RUhEBqNyBOuXX37B4cOHq+z38PDQW+i4evUqgoODkZWVhd9++w2BgYGafRcvXsRLL72EBx54AH/88QccHBw0+95//30sXLgQH3/8Mf773/8CqPg6RowYgd27d+PQoUMICAjQtC8pKcHmzZthaWlZ6xwtIjJsDFhEZBBKS0tx4cIFAMDKlSurbTNlyhS9BKyTJ09i9OjRkMvlOHjwIB544AGt/Z999hlUKhX+97//aYUrAFiwYAFWrFiBLVu2aAIWAEydOhW7d+/Gl19+qRWwtm/fjqysLISGhkKpVDbq10VEjYcBi4gMwp9//omysjI8+eST913CoJKfn59m1MvMzAxdu3bF8uXL8eCDD+qsrsOHD+Ojjz6Ci4sLdu/eXe0VjMePHwcA/Pbbb/jjjz+q7Dc1NcXly5e1to0bNw6Ojo747rvv8Mknn8DGxgYA8OWXXwKoCGBEZLwYsIjIIFQGpe7du9epfVlZGS5cuIAvv/wSo0aNQl5eHt544w1MmDABqampUCgUOqnr7NmzyM/Px6hRo2pc3DMzMxMA8N5779X5uAqFAs8++yxWr16Nb775BtOmTUNCQgL27t2LTp06YejQoboon4j0hJPcicggVE5w79GjR53a//XXXygtLcXIkSPh6uqKTp064emnn0ZOTg6SkpJ0VteLL76IqVOn4ttvv0VISAjUanWVNpWn8nJzcyGKYo2Pe1We7qwctYqOjkZ5eTlHr4iaAQYsIjIIlSNYdQ1YZ86cQbt27eDm5gYAuHbtGt5//324uLigTZs2OqtLJpPhiy++wPTp07Fly5ZqQ1a/fv0A/HOqsK66d+8Of39/xMTE4PLly4iOjoZcLkdoaKjO6ici/WDAIiK9Kysrw/nz52Fra4u2bdvW6TVnzpxBfHw8rK2tYW5ujs6dO+PKlSv46quvYGKi29kPgiBg3bp1mDFjBrZs2YJnnnlGK2TNmjULJiYmeOmll5CQkFDl9dnZ2Th79my1x64cxZo+fTr+/vtvjB49Gq1bt9Zp/UTU9DgHi4j07q+//kJJSQns7e1rXMHdwsICr7zyiub52bNnMX/+fEyfPh1paWmYNWsWZs6ciREjRjRKjYIg4LPPPtOELVEU8fXXX8PExAS+vr6IjIzEzJkz0aVLF4wePRodOnRAbm4u/v77bxw8eBCTJ0/G2rVrqxz3qaeeQnh4OI4ePQqAa18RNRcMWESkd5Xzr1JSUrB48eJq2wwYMEATsMrLy3Hu3Dm8+uqr6NixIzp27IiIiAhMmjQJzz33HMzNzRulTkEQsHbtWshkMqxduxaiKGLz5s0wMTHBjBkz8MADD2DFihU4dOgQfvrpJ9ja2sLT0xNz586t8bSfUqnExIkTsXHjRri4uODhhx9ulNqJqGkJYnUzL4mIDNjly5fh4+OD69evo3379gAqTjPa29tj48aNeOyxx+57jOjoaEyZMgVRUVFNerNnXai8mvHmzZt6rYOIasY5WERkdM6cOQNLS0utZRNMTU0xdOhQ/Pjjj/U61pQpUyAIgsGvmh4XFwdBECAIAm7duqXvcojoPniKkIiMztmzZ+Ht7Q2ZTPv/iMOHD8fSpUuhVqshl8trPcYDDzyARYsWaZ4b+k2VHRwctOq1s7PTXzFEdF88RUhERESkYzxFSERERKRjDFhEREREOsaARURERKRjDFhEREREOsaARURERKRjDFhEREREOsaARURERKRjDFhEREREOsaARURERKRjDFhEREREOsaARURERKRjDFhEREREOvb/VuvznvMXQYoAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.loglog(E_Rs*1e6, abs(sm_spec), label=\"SNuDD Standard Model\") # abs because of the residual imaginary part from the density matrix. All elements have zero imaginary part anyway.\n",
    "plt.legend()\n",
    "plt.title(r\"Theoretical CE$\\nu$NS rate in Xe\")\n",
    "plt.ylabel(r'${\\rm d}R/{\\rm d}E_R\\,\\,\\left[{\\rm ton}^{-1}\\,{\\rm yr}^{-1}\\,{\\rm keV}^{-1}\\right]$', size=14)\n",
    "plt.xlabel(r'$E_{R}\\,\\,\\left[{\\rm keV}\\right]$', size=14)\n",
    "\n",
    "plt.ylabel(r\"Event rate [events/(ton$\\cdot$year$\\cdot$keV)]\")\n",
    "\n",
    "plt.xlim(xmin=1e-1, xmax=1e1)\n",
    "plt.ylim(ymax=1e6, ymin=1e-5)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "41c813f9-a5d8-4da2-ad1e-9277757df7d7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "\"\"\" IN SI UNIT \"\"\"\n",
    "\"\"\" Reduced Planck Constant in J·s \"\"\"\n",
    "HBAR = 1.054571800E-34\n",
    "\"\"\" Speed of Light in m/s \"\"\"\n",
    "C = 2.99792458E8\n",
    "\"\"\" Avogaro's Number in 1/mol \"\"\"\n",
    "A_V = 6.022140857E23\n",
    "\"\"\" Charge of Electron in C \"\"\"\n",
    "CH_E = 1.6021766208E-19\n",
    "\n",
    "\"\"\" keV in J \"\"\"\n",
    "KEV = 1e3 * CH_E\n",
    "\"\"\" GeV in J \"\"\"\n",
    "GEV = 1e9 * CH_E\n",
    "\n",
    "GeV = 1.0\n",
    "MeV = GeV / 1e3\n",
    "keV = GeV / 1e6\n",
    "eV = GeV / 1e9\n",
    "\n",
    "# GeV ^ -1\n",
    "m = 1.0 / (HBAR * C / GEV) / GeV\n",
    "cm = m / 1e2\n",
    "km = m * 1e3\n",
    "mm = m / 1e3\n",
    "um = m / 1e6\n",
    "nm = m / 1e9\n",
    "pm = m / 1e12\n",
    "fm = m / 1e15\n",
    "\n",
    "# GeV\n",
    "kg = 1.0 * C**2 / GEV * GeV\n",
    "ton = kg * 1e3\n",
    "g = kg * 1e-3\n",
    "\n",
    "# GeV ^ -1\n",
    "sec = 1.0 * C * m\n",
    "dy = 24.0 * 3600.0 * sec\n",
    "yr = 365.0 * dy\n",
    "\n",
    "# Fermi's const\n",
    "Gf = 1.1663787e-5 / GeV**2\n",
    "\n",
    "# Neutron mass\n",
    "mn = 0.9314940954 * GeV\n",
    "\n",
    "# Mixing angle\n",
    "sin2the = 0.2223"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "e3f0fd69-1d5a-4e29-aeb3-126b330ea26a",
   "metadata": {},
   "outputs": [],
   "source": [
    "NSI_matrix = np.array([[0.0, 0.0, 0.0],\n",
    "                       [0.0, 0.0, 0.0],\n",
    "                       [0.0, 0.0, 0.0]])\n",
    "\n",
    "NSI_eta = 0\n",
    "NSI_phi = 0\n",
    "\n",
    "NSI_model = GeneralNSI(NSI_matrix, NSI_eta, NSI_phi)\n",
    "\n",
    "Xe_nucleus.update_model(NSI_model)\n",
    "Xe_nucleus.prepare_density()\n",
    "\n",
    "NSI_spec=Xe_nucleus.spectrum(E_Rs, nu = '8B')\n",
    "\n",
    "acc_data = pd.read_csv('acceptance.csv').drop(columns='Unnamed: 0')\n",
    "\n",
    "accept = np.array(acc_data['acceptance'])\n",
    "energy = np.array(acc_data['energy'] * keV)\n",
    "accept = interp1d(energy, accept)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "b329ac42-973a-42ec-9baa-b538fac95553",
   "metadata": {},
   "outputs": [],
   "source": [
    "Xe_electron = Electron(Xe_nucleus, None, None)\n",
    "Xe_electron.update_model(NSI_model)\n",
    "Xe_electron.prepare_density()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "91bece62-534c-4037-8f1d-d72926aef461",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "90% UL B8 solar neutrino rate: 13.12 per ty\n"
     ]
    }
   ],
   "source": [
    "er = energy\n",
    "tot_rate = (np.abs(Xe_nucleus.spectrum(er, nu = '8B')) * accept(er).astype(\"float64\") * (yr**(-1) * keV**(-1) * ton**(-1)))\n",
    "\n",
    "delta_energy = er[1] - er[0]\n",
    "total_rate = np.sum(tot_rate) * delta_energy\n",
    "upper_lim_B8_flux = 5.25e6 * 2.3 * cm**-2 * sec**-1\n",
    "lower_lim_B8_flux = 5.25e6 * 0.9 * cm**-2 * sec**-1\n",
    "tot_flux = 5.46e6 * cm**-2 * sec**-1\n",
    "upper_lim_total_rate = total_rate / tot_flux * upper_lim_B8_flux\n",
    "lower_lim_total_rate = total_rate / tot_flux * lower_lim_B8_flux\n",
    "print(\"{:} {:.2f} per ty\".format(\"90% UL B8 solar neutrino rate:\", upper_lim_total_rate / (yr**(-1) * ton**(-1))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "b44e8ba8-af12-4746-954e-f551def04de5",
   "metadata": {},
   "outputs": [],
   "source": [
    "er_energies = np.linspace(24, 50, 27)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "82dddd3a-3dfb-4d4b-ba8e-907d7e05bae6",
   "metadata": {},
   "outputs": [],
   "source": [
    "pp_sum = jnp.real(jnp.sum(Xe_electron.spectrum(er_energies, nu = 'pp')))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e153576a-e0fd-4092-affd-0280aefecdda",
   "metadata": {},
   "source": [
    "## Define the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "deff05a9-1996-46d9-88c3-1d4155ede6e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import jax\n",
    "import jax.numpy as jnp\n",
    "\n",
    "import numpyro\n",
    "from numpyro import distributions as dist, infer\n",
    "from numpyro import handlers\n",
    "\n",
    "\n",
    "# @handlers.reparam(config={\"phi\": infer.reparam.CircularReparam()})\n",
    "def nsi_rate(rate = None):\n",
    "    eps_ee = numpyro.sample(\"eps_ee\", dist.Uniform(-5.0, 5.0))\n",
    "    eps_em = numpyro.sample(\"eps_em\", dist.Uniform(-5.0, 5.0))\n",
    "    eps_et = numpyro.sample(\"eps_et\", dist.Uniform(-5.0, 5.0))\n",
    "    eps_mm = numpyro.sample(\"eps_mm\", dist.Uniform(-5.0, 5.0))\n",
    "    eps_mt = numpyro.sample(\"eps_mt\", dist.Uniform(-5.0, 5.0))\n",
    "    eps_tt = numpyro.sample(\"eps_tt\", dist.Uniform(-5.0, 5.0))\n",
    "    sineta = numpyro.sample(\"sineta\", dist.Uniform(-1., 1.))\n",
    "    phi = numpyro.sample(\"phi\", dist.Uniform(-0.5 * jnp.pi, 0.5 * jnp.pi))\n",
    "    \n",
    "    NSI_matrix = jnp.asarray([[eps_ee, eps_em, eps_et],\n",
    "                       [eps_em, eps_mm, eps_mt],\n",
    "                       [eps_et, eps_mt, eps_tt]])\n",
    "    NSI_eta = jnp.arcsin(sineta)\n",
    "    NSI_phi = phi\n",
    "\n",
    "    NSI_model = GeneralNSI(NSI_matrix, NSI_eta, NSI_phi)\n",
    "    Xe_nucleus.update_model(NSI_model)\n",
    "    Xe_nucleus.prepare_density()\n",
    "\n",
    "    tot_rate = jnp.sum(jnp.abs(Xe_nucleus.spectrum(er, nu = '8B')) * accept(er) * (yr**(-1) * keV**(-1) * ton**(-1))) * delta_energy\n",
    "    tot_lambda = tot_rate * 0.6 * ton * yr + 5.38\n",
    "    \n",
    "    # Xe_electron = Electron(Xe_nucleus, None, None)\n",
    "    Xe_electron.update_model(NSI_model)\n",
    "    Xe_electron.prepare_density()\n",
    "    pp_sum_NSI = jnp.real(jnp.sum(Xe_electron.spectrum(er_energies, nu = 'pp')))\n",
    "    \n",
    "    measured_pp_ratio = numpyro.sample(\"measured_pp_ratio\", dist.TruncatedNormal(1.72, 1.72, low=0.0), obs=pp_sum_NSI/pp_sum)\n",
    "    \n",
    "    \n",
    "    numpyro.sample(\"rate\", dist.Poisson(tot_lambda), obs=rate)\n",
    "\n",
    "    \n",
    "    # numpyro.factor('rate_lh',jax.scipy.stats.poisson.logpmf(rate, tot_lambda))\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8db712ec-64ad-4117-8e14-1cd76710473f",
   "metadata": {},
   "source": [
    "## Load NeuTra Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "08bb3d71",
   "metadata": {},
   "outputs": [],
   "source": [
    "guide = autoguide.AutoBNAFNormal(\n",
    "    nsi_rate,\n",
    "    hidden_factors=[10, 10],\n",
    "    num_flows=3\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "128f2f22-3bf8-486f-aeaa-5e536b35d346",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsIAAAH5CAYAAACLR1gqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAoOklEQVR4nO3df5BV9X3w8c/dBdS1sESp4IogtsagCygLMeAPCvShMorjmFgTLdJoxmIhFailJrYxj08TzMxTx7ZBqJbRmIwDzvhjnJYkwpCAEVIVwahU0UhFIyzFR1iUKmb3PH/UvbgC6uI53O/lvl4zzIS7Z3e/y4nJm6+f+z2lLMuyAACAGlNX6QUAAEAlCGEAAGqSEAYAoCYJYQAAapIQBgCgJglhAABqkhAGAKAm9aj0AqpNR0dHvP7669G7d+8olUqVXg4AAB+SZVns2rUrmpqaoq7uwPu+QribXn/99TjxxBMrvQwAAD7Gq6++GgMHDjzgx4VwN/Xu3Tsi/ucPtk+fPhVeDQAAH9bW1hYnnnhiudsORAh3U+c4RJ8+fYQwAEDCPm6M1ZvlAACoSUIYAICaJIQBAKhJQhgAgJokhAEAqElCGACAmiSEAQCoSUIYAICaJIQBAKhJQhgAgJokhAEAqElCGACAmiSEAQCoSTUbwrt27YrRo0fHGWecEcOGDYs777yz0ksCAOAQ6lHpBVRKQ0NDrFy5MhoaGmL37t3R3Nwcl1xySRx77LGVXhoAAIdAze4I19fXR0NDQ0REvPPOO9He3h5ZllV4VQAAHCrdDuHf/va38Td/8zcxZMiQOOqoo+Lkk0+Om2++OTo6OnJb1KpVq2LKlCnR1NQUpVIpHnroof1ed/vtt8eQIUPiyCOPjJaWlnj00Ue79X127NgRI0aMiIEDB8bcuXOjX79+OaweAIBq0O3RiO9973uxcOHC+MEPfhCnn356PPnkk/HVr341Ghsb47rrrtvn+sceeyw+//nPR8+ePbu8/vzzz0ffvn1jwIAB+3zO22+/HSNGjIivfvWr8cUvfnG/61iyZEnMmjUrbr/99jj77LPjn//5n2Py5MmxYcOGGDRoUEREtLS0xLvvvrvP5z7yyCPR1NQUffv2jaeffjpaW1vjkksuiS996UvRv3//7v6RFGbLzv+Or9+7rtLLAICKOqJnXVw/6dQ4c9BnKr0UDjOlrJvzABdeeGH0798/Fi1aVH7ti1/8YjQ0NMQPf/jDLtd2dHTEyJEj45RTTonFixdHfX19RERs3Lgxxo0bF7Nnz465c+d+9AJLpXjwwQfj4osv7vL6WWedFSNHjowFCxaUXxs6dGhcfPHFMW/evO78SBERce2118aECRPi0ksv3e/H58+fH/Pnz4/29vbYuHFj7Ny5M/r06dPt79Md/7n97fiD//vzQr8HAFSDy0adGN/70vBKL4Mq0dbWFo2NjR/ba93eET7nnHNi4cKFsXHjxvjsZz8bTz/9dPziF7+I2267bZ9r6+rqYunSpXHeeefFlVdeGT/84Q9j06ZNMWHChLjooos+NoIPZM+ePbF27dq44YYburw+adKkWL169Sf6Gq2trXHUUUdFnz59oq2tLVatWhXXXnvtAa+fMWNGzJgxo/wHeyj8bu8jYuGfjDwk3wsAUvTIhtZ44KnfxHs5jmBCp26H8F//9V/Hzp0743Of+1zU19dHe3t7fOc734mvfOUr+72+qakpVqxYEeedd15cfvnlsWbNmpg4cWIsXLjwoBe9ffv2aG9v32eMoX///rF169ZP9DVee+21uPrqqyPLssiyLGbOnBnDh6f1N82jj+gR5zcfX+llAEDFbP5/uyu9BA5j3Q7hJUuWxI9+9KO499574/TTT4/169fHrFmzoqmpKaZNm7bfzxk0aFDcc889MW7cuDj55JNj0aJFUSqVPvXiP/w1siz7xF+3paUl1q9f/6nXAAAUpxTv//+6g50oQLdPjfirv/qruOGGG+LLX/5yDBs2LKZOnRqzZ8/+yLnc1tbWuOaaa2LKlCmxe/fumD179qdadL9+/aK+vn6f3d9t27Yl9WY3ACAfOpgidDuEd+/eHXV1XT+tvr7+gMenbd++PSZOnBhDhw6NBx54IFasWBH33XdfXH/99Qe34ojo1atXtLS0xLJly7q8vmzZshg7duxBf10AIC05/AtkOKBuj0ZMmTIlvvOd78SgQYPi9NNPj3Xr1sWtt94aV1111T7XdnR0xPnnnx+DBw+OJUuWRI8ePWLo0KGxfPnyGD9+fJxwwgn73R1+66234qWXXir/ftOmTbF+/fo45phjykejzZkzJ6ZOnRqjRo2KMWPGxB133BGbN2+O6dOnd/dHAgAS56FXFKHbIfxP//RP8bd/+7fx53/+57Ft27ZoamqKP/uzP4tvfetb+1xbV1cX8+bNi3PPPTd69epVfn3YsGGxfPnyAz7O+Mknn4zx48eXfz9nzpyIiJg2bVrcfffdERFx2WWXxRtvvBE333xzbNmyJZqbm2Pp0qUxePDg7v5IAEDiZDBF6PY5wrXuk55LBwB8eot+sSn+z79uiItGNMU/fuXMSi+HKvFJe63bM8IAAIeKEWGKJIQBgOT519cUQQgDAMnqPDXCJCdFEMIAQLKMRlAkIQwAJKvzibH2gymCEAYA0qeEKYAQBgCS5clyFEkIAwDJ6uzgzJYwBRDCAEDyHBpBEYQwAJAusxEUSAgDAMkqj0bYEaYAQhgASJ4ZYYoghAGAZJmMoEhCGABIVun94QijERRBCAMAydPBFEEIAwDJ6hyNsCNMEYQwAJAsI8IUSQgDAFXAljD5E8IAQLKMRlAkIQwAJKtkOIICCWEAIF2dO8KVXQWHKSEMACQvMxtBAYQwAJAsgxEUSQgDAMkqvf9uOfvBFEEIAwDJMxlBEYQwAJAsoxEUSQgDAMkqOTWCAglhACB5To2gCEIYAEhWyWwEBRLCAECyPFmOIglhACB5JiMoghAGAJK1981ySpj8CWEAAGqSEAYAkmc0giIIYQAgWeVHLAthCiCEAYBkOTOCIglhACBZ3ixHkYQwAJA8oxEUQQgDAMnyQA2KJIQBgGTtHY2A/AlhACB9SpgCCGEAIFkGIyiSEAYAkuXUCIokhAGA5Dk1giIIYQAgYe8/Wa7Cq+DwJIQBgGSVDAlTICEMACQvMxtBAYQwAJCszg1hGUwRhDAAkKyS2QgKJIQBgOSZjKAIQhgASJbRCIokhAGAZJmMoEhCGABIVjmEzUZQACEMACRPBlMEIQwAJKsUZiMojhAGANL1fgebjKAIQhgASF5mOIICCGEAIFneK0eRhDAAkCxPlqNIQhgASJ4dYYoghAGAZHmyHEUSwgBAskxGUCQhDAAkLzMbQQGEMACQLA/UoEhCGABIltEIiiSEAYDkmYygCEIYAEjW3lMjlDD5E8IAQLqMRlAgIQwAJKvzzXJGIyiCEAYAkqeDKYIQBgCS5dQIiiSEAYBkld8sZzaCAghhACB5MpgiCGEAIFmlztkIJUwBhDAAkCwzwhRJCAMAybMhTBGEMACQLG+Wo0hCGABIltEIiiSEAYDk2Q+mCEIYAEiYRyxTHCEMACTLaARFEsIAQPIywxEUQAgDAMnae2pERZfBYUoIAwDJKpmNoEBCGABIlh1hiiSEAQCoSUIYAEhW52SEJ8tRBCEMACSrFGaEKY4QBgCSZz+YIghhACBZe0cjKrsODk9CGACAmiSEAYDkebIcRRDCAECyjEZQJCEMACTLqREUSQgDAMmzIUwRhDAAkCyjERRJCAMAySqZjKBAQhgAqAK2hMmfEAYAktX5ZjmjERRBCAMAyTIaQZGEMACQrM4OtiFMEYQwAJC8zGwEBRDCAECyysenVXYZHKaEMACQMEPCFEcIAwDJMxlBEYQwAJCsvU+WU8LkTwgDAMkyGEGRhDAAkDz7wRRBCAMAySo5NoICCWEAIFlGIyiSEAYAkmdDmCIIYQAgWU6NoEhCGABIVslwBAUSwgBA8uwHUwQhDAAka+9oRGXXweFJCAMAycvsCVMAIQwAJKtkRJgCCWEAIHlGIyiCEAYAktX5ZDkdTBGEMACQLJMRFEkIAwDpsyVMAYQwAJCs8vFpSpgCCGEAIFmeLEeRhDAAkDynRlAEIQwAJGvvaATkTwgDAMkyGEGRhDAAkLzMbAQFEMIAQLqMRlAgIQwAJKvz1AgbwhRBCAMAUJOEMACQrJJ3y1EgIQwAJOuDHewNc+RNCAMAySrZEqZAQhgAqAo2hMmbEAYAktVlNKJiq+BwJYQBgGSZjKBIQhgAqAreLEfehDAAkKzSB4YjZDB5E8IAQLqMRlAgIQwAVAWTEeRNCAMAyfrgm+UywxHkTAgDAMkyGUGRhDAAUBWMRpC3mg3hXbt2xejRo+OMM86IYcOGxZ133lnpJQEAH+IRyxSpR6UXUCkNDQ2xcuXKaGhoiN27d0dzc3Nccsklceyxx1Z6aQDA+7o8Wc6OMDmr2R3h+vr6aGhoiIiId955J9rb2x3UDQBQQ7odwieddFKUSqV9fs2YMSO3Ra1atSqmTJkSTU1NUSqV4qGHHtrvdbfffnsMGTIkjjzyyGhpaYlHH320W99nx44dMWLEiBg4cGDMnTs3+vXrl8PqAYC8ODWCInU7hJ944onYsmVL+deyZcsiIuLSSy/d7/WPPfZYvPfee/u8/vzzz8fWrVv3+zlvv/12jBgxIr7//e8fcB1LliyJWbNmxY033hjr1q2Lc889NyZPnhybN28uX9PS0hLNzc37/Hr99dcjIqJv377x9NNPx6ZNm+Lee++N1tbWT/znAAAUr8uT5XQwOStln3IeYNasWfGv//qv8eKLL+4z0N7R0REjR46MU045JRYvXhz19fUREbFx48YYN25czJ49O+bOnfvRCyyV4sEHH4yLL764y+tnnXVWjBw5MhYsWFB+bejQoXHxxRfHvHnzuv1zXHvttTFhwoQDBn2ntra2aGxsjJ07d0afPn26/X0AgE/unffa43N/+5OIiHjuf/9RHH1Ezb69iW74pL32qWaE9+zZEz/60Y/iqquu2u+7Ouvq6mLp0qWxbt26uPLKK6OjoyN+/etfx4QJE+Kiiy762Aj+qO+7du3amDRpUpfXJ02aFKtXr/5EX6O1tTXa2toi4n/+sFatWhWnnnrqAa+fP39+nHbaaTF69OiDWjMA8OnYECZvn+qvVQ899FDs2LEj/vRP//SA1zQ1NcWKFSvivPPOi8svvzzWrFkTEydOjIULFx70992+fXu0t7dH//79u7zev3//A45bfNhrr70WV199dWRZFlmWxcyZM2P48OEHvH7GjBkxY8aM8t8wAIBDy5vaydunCuFFixbF5MmTo6mp6SOvGzRoUNxzzz0xbty4OPnkk2PRokW5nAv44a+RZdkn/rotLS2xfv36T70GAKA4jhGmSAc9GvHKK6/E8uXL42tf+9rHXtva2hrXXHNNTJkyJXbv3h2zZ88+2G8bERH9+vWL+vr6fXZ/t23bts8uMQBweLAfTN4OOoTvuuuuOO644+KCCy74yOu2b98eEydOjKFDh8YDDzwQK1asiPvuuy+uv/76g/3W0atXr2hpaSmfWNFp2bJlMXbs2IP+ugBAWpwaQZEOajSio6Mj7rrrrpg2bVr06HHgL9HR0RHnn39+DB48OJYsWRI9evSIoUOHxvLly2P8+PFxwgkn7Hd3+K233oqXXnqp/PtNmzbF+vXr45hjjolBgwZFRMScOXNi6tSpMWrUqBgzZkzccccdsXnz5pg+ffrB/EgAQIKMRlCkgwrh5cuXx+bNm+Oqq676yOvq6upi3rx5ce6550avXr3Krw8bNiyWL19+wMcZP/nkkzF+/Pjy7+fMmRMREdOmTYu77747IiIuu+yyeOONN+Lmm2+OLVu2RHNzcyxdujQGDx58MD8SAJA6O8Lk7FOfI1xrnCMMAIfOb9s74vdv/HFERKz/1v+Kvg29PuYz4BCdIwwAUKQPngZl6468CWEAAGqSEAYAkvXB98rZECZvQhgASNYHT43wtibyJoQBAKhJQhgASFaXN8tVcB0cnoQwAFAVTEaQNyEMACTN0+UoihAGAKpCZjiCnAlhACBp5Q1hHUzOhDAAkLSS2QgKIoQBgKpgQ5i8CWEAIGmd+8FOjSBvQhgASJrJCIoihAGAquDUCPImhAGApJXeH44wGkHehDAAkLb3RyN0MHkTwgAA1CQhDAAkbe+pEfaEyZcQBgCS1nlqhA4mb0IYAICaJIQBgKSVwkHCFEMIAwBJMxpBUYQwAJA0+8EURQgDAFXBk+XImxAGAJJWKnmyHMUQwgBA0oxGUBQhDABUBRvC5E0IAwBpK58aIYXJlxAGAJJWfsRyRVfB4UgIAwBQk4QwAJA0p0ZQFCEMACStVD42QgmTLyEMAEBNEsIAQNLKb5azIUzOhDAAkLTyjHCF18HhRwgDAFCThDAAkDSjERRFCAMASes8NSIzHEHOhDAAkLjSx18CB0EIAwBVwWgEeRPCAEDSyqMRQpicCWEAIGkGIyiKEAYAqoI3y5E3IQwAJM1oBEURwgBA0kqGIyiIEAYAoCYJYQAgaUYjKIoQBgCSVn7EsjfLkTMhDABATRLCAEDSSu/PRhiNIG9CGACoCjqYvAlhAABqkhAGAJK299QIe8LkSwgDAEkrh3Bll8FhSAgDAEnzZDmKIoQBgKpgMoK8CWEAIGml8oawEiZfQhgASFr5yXI6mJwJYQAAapIQBgCSVn6yXIXXweFHCAMASTMaQVGEMAAANUkIAwBp82Q5CiKEAYCklUcjKroKDkdCGACAmiSEAYCklU+NsCVMzoQwAJC0vaMRSph8CWEAAGqSEAYAklbybjkKIoQBgKSVwpPlKIYQBgCSVt4RhpwJYQCgKjg1grwJYQCgKjg1grwJYQAgac4RpihCGACAmiSEAYCkOT2NoghhACBpnadGZGYjyJkQBgCgJglhACBp5R3hyi6Dw5AQBgCS1vlkOSVM3oQwAAA1SQgDAEnbOxphS5h8CWEAIGnl49N0MDkTwgAA1CQhDACkzSOWKYgQBgCS5slyFEUIAwBJ82Q5iiKEAQCoSUIYAEia0QiKIoQBgKSVvFmOgghhAABqkhAGAJJWKv8nW8LkSwgDAEnbe2pEZdfB4UcIAwBQk4QwAJC00vvDETaEyZsQBgDSZjSCgghhAABqkhAGAJK294EatoTJlxAGAJLm1AiKIoQBAKhJQhgASJpTIyiKEAYAkrZ3NEIKky8hDAAkrVT6+GvgYAhhAABqkhAGAJJWnhE2GUHOhDAAkLTyjLC3y5EzIQwAQE0SwgBAVTAaQd6EMACQtFLJjDDFEMIAANQkIQwAJK3zGGEbwuRNCAMASfNkOYoihAEAqElCGABImtEIiiKEAYCklfY+UQNyJYQBgKrgyXLkTQgDAEkrffwlcFCEMACQtL2nRlR2HRx+hDAAkLj3nyxX4VVw+BHCAADUJCEMACTNaARFEcIAQNL2niOshMmXEAYAoCYJYQAgaUYjKIoQBgCSVnJqBAURwgAA1CQhDAAkrVR+t5w9YfIlhAGApJVnhCu7DA5DQhgAgJokhAGApJXfLGdLmJwJYQAgbeXj05Qw+RLCAEBVkMHkTQgDAEkrffwlcFCEMACQtFLJjDDFEMIAQNLKxwhXdBUcjoQwAAA1SQgDAEkrOTWCgghhACBp3ixHUYQwAAA1SQgDAElzagRFEcIAQNL2nhqhhMmXEAYAoCYJYQAgbeVTIyq7DA4/QhgASFrp/RLWweRNCAMAVcGOMHkTwgBA0koOEqYgQhgASJpTIyiKEAYAqoLRCPImhAGApBmNoChCGABIWimUMMUQwgBA0krlc4TNRpAvIQwAQE0SwgBA0kqeLEdBhDAAkDhPlqMYQhgAgJokhAGApBmNoChCGABImifLURQhDABATRLCAEDSjEZQFCEMACSt5NQICiKEAYDqYEuYnAlhACBpnaMRkDchDAAkbe+pEZAvIQwAVAWTEeRNCAMASSuZjaAgQhgAqAoeqEHehDAAkDTnCFOUmg3hXbt2xejRo+OMM86IYcOGxZ133lnpJQEAcAj1qPQCKqWhoSFWrlwZDQ0NsXv37mhubo5LLrkkjj322EovDQD4AA/UoCg1uyNcX18fDQ0NERHxzjvvRHt7e2T+nQsAJMdoBEU5qBD+zW9+E3/yJ38Sxx57bDQ0NMQZZ5wRa9euzW1Rq1atiilTpkRTU1OUSqV46KGH9nvd7bffHkOGDIkjjzwyWlpa4tFHH+3W99mxY0eMGDEiBg4cGHPnzo1+/frlsHoAAKpBt0P4zTffjLPPPjt69uwZP/7xj2PDhg3x93//99G3b9/9Xv/YY4/Fe++9t8/rzz//fGzdunW/n/P222/HiBEj4vvf//4B17FkyZKYNWtW3HjjjbFu3bo499xzY/LkybF58+byNS0tLdHc3LzPr9dffz0iIvr27RtPP/10bNq0Ke69995obW3txp8EAHAo7H2ghi1h8tXtGeHvfe97ceKJJ8Zdd91Vfu2kk07a77UdHR0xY8aMOOWUU2Lx4sVRX18fEREbN26M8ePHx+zZs2Pu3Ln7fN7kyZNj8uTJH7mOW2+9Na6++ur42te+FhERt912W/z0pz+NBQsWxLx58yIiPvEudf/+/WP48OGxatWquPTSS/d7zfz582P+/PnR3t7+ib4mAJCPkkfLUZBu7wg//PDDMWrUqLj00kvjuOOOizPPPPOAJy7U1dXF0qVLY926dXHllVdGR0dH/PrXv44JEybERRddtN8I/iT27NkTa9eujUmTJnV5fdKkSbF69epP9DVaW1ujra0tIiLa2tpi1apVceqppx7w+hkzZsSGDRviiSeeOKg1AwCfjg4mb90O4ZdffjkWLFgQp5xySvz0pz+N6dOnx1/8xV/EPffcs9/rm5qaYsWKFfHYY4/F5ZdfHhMmTIiJEyfGwoULD3rR27dvj/b29ujfv3+X1/v373/AcYsPe+211+K8886LESNGxDnnnBMzZ86M4cOHH/SaAIBieLIcRen2aERHR0eMGjUqvvvd70ZExJlnnhnPPfdcLFiwIK688sr9fs6gQYPinnvuiXHjxsXJJ58cixYtyuW/1B/+GlmWfeKv29LSEuvXr//UawAAilWejHBsBDnr9o7w8ccfH6eddlqX14YOHdrlTWof1traGtdcc01MmTIldu/eHbNnz+7+Sj+gX79+UV9fv8/u77Zt2/bZJQYADg86mLx1O4TPPvvseOGFF7q8tnHjxhg8ePB+r9++fXtMnDgxhg4dGg888ECsWLEi7rvvvrj++usPbsUR0atXr2hpaYlly5Z1eX3ZsmUxduzYg/66AECCTEZQkG6PRsyePTvGjh0b3/3ud+OP//iP4/HHH4877rgj7rjjjn2u7ejoiPPPPz8GDx4cS5YsiR49esTQoUNj+fLlMX78+DjhhBP2uzv81ltvxUsvvVT+/aZNm2L9+vVxzDHHxKBBgyIiYs6cOTF16tQYNWpUjBkzJu64447YvHlzTJ8+vbs/EgCQME+WoyjdDuHRo0fHgw8+GN/4xjfi5ptvjiFDhsRtt90WV1xxxT7X1tXVxbx58+Lcc8+NXr16lV8fNmxYLF++/ICPM37yySdj/Pjx5d/PmTMnIiKmTZsWd999d0REXHbZZfHGG2/EzTffHFu2bInm5uZYunTpAXemAYDqZjSCvJUyk+fd0tbWFo2NjbFz587o06dPpZcDAIe97/3k+Vjw81/HVWcPiW9NOe3jP4Ga90l77aAesQwAcKh4shxFEcIAQNI6T0b177DJmxAGAKAmCWEAIGkl56dRECEMACRt72iE2QjyJYQBgKogg8mbEAYAkmYwgqIIYQAgbe/PRpiMIG9CGACoCs4RJm9CGABImtEIiiKEAYCkeaAGRRHCAEBV0MHkTQgDAEnzQA2KIoQBgKQZjaAoQhgAqBJKmHwJYQAgaQYjKIoQBgCSZjSCoghhAKAqCGHyJoQBgKSVSoYjKIYQBgCqgkcskzchDAAkzYwwRRHCAEBV0MHkTQgDAEnzZDmKIoQBgKQZjaAoQhgAqAreLEfehDAAkDSDERRFCAMASSsfI2xDmJwJYQCgKuhg8iaEAYCkOTWCoghhACBpe0+NsCdMvoQwAFAVZDB5E8IAANQkIQwAJK30/myEyQjyJoQBgKqgg8mbEAYAklY+RtiWMDkTwgBA0kpOT6MgQhgASJoHy1EUIQwAVAclTM6EMACQtJLZCAoihAGApJWfLGdLmJwJYQCgKjg0grwJYQAgaQYjKIoQBgDS5slyFEQIAwBVwYwweRPCAEDSjEZQFCEMACStfGqEDWFyJoQBgKqgg8mbEAYAklYyHEFBhDAAkDSjERRFCAMAVUIJky8hDAAkrXMwwo4weRPCAEDSSkaEKYgQBgCS1vlmORvC5E0IAwBVITMbQc6EMACQNqMRFEQIAwBJK79ZrqKr4HAkhAGAqmAygrwJYQAgaSXHRlAQIQwAJM1oBEURwgBAVXBqBHkTwgBA0kxGUBQhDAAkTQhTFCEMAFQFkxHkTQgDAEnb+4hlJUy+hDAAkDSjERRFCAMAVcFoBHkTwgBAVRDC5E0IAwBJ82Q5iiKEAYCk7X2ynC1h8iWEAYCqYDSCvAlhACBpJiMoSo9KLwAA4KN0niP83Ott8aUFqyu8Gg7W+c0D4mvnnlzpZXQhhAGApA1oPDIiIt5697fx5CtvVng1HKzTmvpUegn7EMIAQNJGDuobD/752Ghte6fSS+FTOPGYhkovYR9CGABIWqlUijMHfabSy+Aw5M1yAADUJCEMAEBNEsIAANQkIQwAQE0SwgAA1CQhDABATRLCAADUJCEMAEBNEsIAANQkIQwAQE0SwgAA1CQhDABATRLCAADUJCEMAEBNEsIAANQkIQwAQE0SwgAA1KQelV5AtcmyLCIi2traKrwSAAD2p7PTOrvtQIRwN+3atSsiIk488cQKrwQAgI+ya9euaGxsPODHS9nHpTJddHR0xOuvvx69e/eOUqlU+Pdra2uLE088MV599dXo06dP4d+P/LmH1c89rH7uYfVzD6vfobyHWZbFrl27oqmpKerqDjwJbEe4m+rq6mLgwIGH/Pv26dPHP/hVzj2sfu5h9XMPq597WP0O1T38qJ3gTt4sBwBATRLCAADUJCGcuCOOOCJuuummOOKIIyq9FA6Se1j93MPq5x5WP/ew+qV4D71ZDgCAmmRHGACAmiSEAQCoSUIYAICaJIQBAKhJQhgAgJokhBN2++23x5AhQ+LII4+MlpaWePTRRyu9pJq1atWqmDJlSjQ1NUWpVIqHHnqoy8ezLItvf/vb0dTUFEcddVT8wR/8QTz33HNdrnn33Xfj61//evTr1y+OPvrouOiii+K1117rcs2bb74ZU6dOjcbGxmhsbIypU6fGjh07Cv7pDn/z5s2L0aNHR+/eveO4446Liy++OF544YUu17iHaVuwYEEMHz68/ESqMWPGxI9//OPyx92/6jNv3rwolUoxa9as8mvuY9q+/e1vR6lU6vJrwIAB5Y9X5f3LSNLixYuznj17ZnfeeWe2YcOG7LrrrsuOPvro7JVXXqn00mrS0qVLsxtvvDG7//77s4jIHnzwwS4fv+WWW7LevXtn999/f/bMM89kl112WXb88cdnbW1t5WumT5+enXDCCdmyZcuyp556Khs/fnw2YsSI7Le//W35mvPPPz9rbm7OVq9ena1evTprbm7OLrzwwkP1Yx62/uiP/ii76667smeffTZbv359dsEFF2SDBg3K3nrrrfI17mHaHn744ezf/u3fshdeeCF74YUXsm9+85tZz549s2effTbLMvev2jz++OPZSSedlA0fPjy77rrryq+7j2m76aabstNPPz3bsmVL+de2bdvKH6/G+yeEE/X5z38+mz59epfXPve5z2U33HBDhVZEpw+HcEdHRzZgwIDslltuKb/2zjvvZI2NjdnChQuzLMuyHTt2ZD179swWL15cvuY3v/lNVldXl/3kJz/JsizLNmzYkEVE9stf/rJ8zZo1a7KIyJ5//vmCf6rasm3btiwispUrV2ZZ5h5Wq8985jPZv/zLv7h/VWbXrl3ZKaecki1btiwbN25cOYTdx/TddNNN2YgRI/b7sWq9f0YjErRnz55Yu3ZtTJo0qcvrkyZNitWrV1doVRzIpk2bYuvWrV3u1xFHHBHjxo0r36+1a9fGe++91+WapqamaG5uLl+zZs2aaGxsjLPOOqt8zRe+8IVobGx033O2c+fOiIg45phjIsI9rDbt7e2xePHiePvtt2PMmDHuX5WZMWNGXHDBBfGHf/iHXV53H6vDiy++GE1NTTFkyJD48pe/HC+//HJEVO/965H7V+RT2759e7S3t0f//v27vN6/f//YunVrhVbFgXTek/3dr1deeaV8Ta9eveIzn/nMPtd0fv7WrVvjuOOO2+frH3fcce57jrIsizlz5sQ555wTzc3NEeEeVotnnnkmxowZE++88078zu/8Tjz44INx2mmnlf/P0f1L3+LFi+Opp56KJ554Yp+P+ecwfWeddVbcc8898dnPfjZaW1vj7/7u72Ls2LHx3HPPVe39E8IJK5VKXX6fZdk+r5GOg7lfH75mf9e77/maOXNm/OpXv4pf/OIX+3zMPUzbqaeeGuvXr48dO3bE/fffH9OmTYuVK1eWP+7+pe3VV1+N6667Lh555JE48sgjD3id+5iuyZMnl//zsGHDYsyYMfF7v/d78YMf/CC+8IUvRET13T+jEQnq169f1NfX7/M3n23btu3zNy0qr/Mdsx91vwYMGBB79uyJN9988yOvaW1t3efr/9d//Zf7npOvf/3r8fDDD8fPfvazGDhwYPl197A69OrVK37/938/Ro0aFfPmzYsRI0bEP/zDP7h/VWLt2rWxbdu2aGlpiR49ekSPHj1i5cqV8Y//+I/Ro0eP8p+x+1g9jj766Bg2bFi8+OKLVfvPoRBOUK9evaKlpSWWLVvW5fVly5bF2LFjK7QqDmTIkCExYMCALvdrz549sXLlyvL9amlpiZ49e3a5ZsuWLfHss8+WrxkzZkzs3LkzHn/88fI1//7v/x47d+503z+lLMti5syZ8cADD8SKFStiyJAhXT7uHlanLMvi3Xffdf+qxMSJE+OZZ56J9evXl3+NGjUqrrjiili/fn2cfPLJ7mOVeffdd+M//uM/4vjjj6/efw5zf/sdueg8Pm3RokXZhg0bslmzZmVHH3109p//+Z+VXlpN2rVrV7Zu3bps3bp1WURkt956a7Zu3brycXa33HJL1tjYmD3wwAPZM888k33lK1/Z75ExAwcOzJYvX5499dRT2YQJE/Z7ZMzw4cOzNWvWZGvWrMmGDRvmyJ8cXHvttVljY2P285//vMuxP7t37y5f4x6m7Rvf+Ea2atWqbNOmTdmvfvWr7Jvf/GZWV1eXPfLII1mWuX/V6oOnRmSZ+5i6v/zLv8x+/vOfZy+//HL2y1/+Mrvwwguz3r17l9ukGu+fEE7Y/Pnzs8GDB2e9evXKRo4cWT7qiUPvZz/7WRYR+/yaNm1almX/c2zMTTfdlA0YMCA74ogjsvPOOy975plnunyN//7v/85mzpyZHXPMMdlRRx2VXXjhhdnmzZu7XPPGG29kV1xxRda7d++sd+/e2RVXXJG9+eabh+inPHzt795FRHbXXXeVr3EP03bVVVeV//fwd3/3d7OJEyeWIzjL3L9q9eEQdh/T1nkucM+ePbOmpqbskksuyZ577rnyx6vx/pWyLMvy32cGAIC0mREGAKAmCWEAAGqSEAYAoCYJYQAAapIQBgCgJglhAABqkhAGAKAmCWEAAGqSEAYAoCYJYQAAapIQBgCgJv1/XE47pDisqcYAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 800x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "N_1 = 4000\n",
    "N_2 = 1000\n",
    "epochs = N_1 + N_2\n",
    "high_lr = 0.8e-2\n",
    "lr_multiplier = 0.7\n",
    "schedule = optax.join_schedules([optax.constant_schedule(high_lr), optax.constant_schedule(high_lr*lr_multiplier)], [N_1])\n",
    "x = np.arange(0, epochs, 10)\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.plot(\n",
    "    x,\n",
    "    [schedule(i) for i in x],\n",
    ")\n",
    "ax.set(yscale='log')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "5052f44d-f679-4e0a-8037-fe41cab5f019",
   "metadata": {},
   "outputs": [],
   "source": [
    "svi = infer.SVI(nsi_rate, guide, optax.adamw(schedule), infer.Trace_ELBO(num_particles=30))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "ed537395-9b6b-46f7-941d-690209b08972",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█| 1/1 [00:21<00:00, 21.81s/it, init loss: 227.8295, avg. loss [1-1]: 227.8\n"
     ]
    }
   ],
   "source": [
    "svi_result = svi.run(jax.random.PRNGKey(0), 1, rate=6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "158b018c-a61a-4b02-a91e-ca533e5225a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# with open(\"params_with_sineta_pp.pickle\", \"wb\") as f:\n",
    "#     pickle.dump(svi_result.params, f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "495b4dfc-56c6-4a53-9791-0a9eba629176",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"params_with_sineta_pp_2.pickle\", \"rb\") as f:\n",
    "    loaded_dict = pickle.load(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "783ab07d-9255-4dba-b963-1e12fd2618ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "neutra = infer.reparam.NeuTraReparam(guide, loaded_dict)#svi_result.params)\n",
    "neutra_model = neutra.reparam(nsi_rate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "5d229227-faf6-41d2-9d6f-3fe76f305334",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "sample: 100%|█| 7000/7000 [1:07:03<00:00,  1.74it/s, 31 steps of size 1.20e-01. \n",
      "sample: 100%|█| 7000/7000 [1:06:27<00:00,  1.76it/s, 31 steps of size 1.20e-01. \n",
      "sample: 100%|█| 7000/7000 [1:08:17<00:00,  1.71it/s, 31 steps of size 1.20e-01. \n",
      "sample: 100%|█| 7000/7000 [1:04:58<00:00,  1.80it/s, 31 steps of size 1.20e-01. \n",
      "sample: 100%|█| 7000/7000 [1:08:43<00:00,  1.70it/s, 31 steps of size 1.20e-01. \n",
      "sample: 100%|█| 7000/7000 [1:07:00<00:00,  1.74it/s, 31 steps of size 1.20e-01. \n",
      "sample: 100%|█| 7000/7000 [1:07:41<00:00,  1.72it/s, 31 steps of size 1.20e-01. \n",
      "sample: 100%|█| 7000/7000 [1:04:32<00:00,  1.81it/s, 31 steps of size 1.20e-01. \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 6h 38min 28s, sys: 2h 19min 33s, total: 8h 58min 2s\n",
      "Wall time: 8h 55min 59s\n"
     ]
    }
   ],
   "source": [
    "kernel = infer.NUTS(neutra_model, target_accept_prob=0.98, step_size=0.12, adapt_step_size=False)\n",
    "\n",
    "sampler = infer.MCMC(kernel,\n",
    "    num_warmup=2000,\n",
    "    num_samples=5000,\n",
    "    num_chains=8,\n",
    "    progress_bar=True,\n",
    "    chain_method='sequential'\n",
    ")\n",
    "\n",
    "%time sampler.run(jax.random.PRNGKey(22), rate = 6, extra_fields = ['num_steps', 'accept_prob'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "2f98a375-d6dd-4955-a6c9-a9c9c72e6bd0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "                           mean       std    median      5.0%     95.0%     n_eff     r_hat\n",
      "auto_shared_latent[0]      0.07      1.04      0.07     -1.59      1.68  32720.95      1.00\n",
      "auto_shared_latent[1]     -0.07      1.01     -0.07     -1.77      1.51  33248.67      1.00\n",
      "auto_shared_latent[2]     -0.08      1.12     -0.09     -1.78      1.65  23459.66      1.00\n",
      "auto_shared_latent[3]     -0.07      1.02     -0.08     -1.72      1.56  29834.50      1.00\n",
      "auto_shared_latent[4]      0.06      1.08      0.06     -1.54      1.91  32626.17      1.00\n",
      "auto_shared_latent[5]      0.03      1.05      0.03     -1.63      1.73  32164.95      1.00\n",
      "auto_shared_latent[6]     -0.05      1.91     -0.03     -2.11      2.05   2814.58      1.00\n",
      "auto_shared_latent[7]     -0.02      1.60     -0.01     -2.03      1.99   6078.68      1.00\n",
      "\n",
      "Number of divergences: 3\n"
     ]
    }
   ],
   "source": [
    "sampler.print_summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "838208f7-30f4-4937-9f92-520b00717f93",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'nrer_az/ntr_ddl_fixed_size_0.98_5k_8.json'"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "arviz_obj = az.from_numpyro(sampler)\n",
    "# az.to_json(arviz_obj, \"nrer_az/ntr_ddl_fixed_size_0.98_5k_8.json\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "fa0f0630-d1ad-4148-a3a6-460ef7dd6b4c",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "sample: 100%|█| 7000/7000 [1:04:19<00:00,  1.81it/s, 31 steps of size 1.20e-01. \n",
      "sample: 100%|█| 7000/7000 [1:06:07<00:00,  1.76it/s, 31 steps of size 1.20e-01. \n",
      "sample: 100%|█| 7000/7000 [1:06:36<00:00,  1.75it/s, 15 steps of size 1.20e-01. \n",
      "sample: 100%|█| 7000/7000 [1:05:14<00:00,  1.79it/s, 15 steps of size 1.20e-01. \n",
      "sample: 100%|█| 7000/7000 [1:06:55<00:00,  1.74it/s, 31 steps of size 1.20e-01. \n",
      "sample: 100%|█| 7000/7000 [1:07:39<00:00,  1.72it/s, 31 steps of size 1.20e-01. \n",
      "sample: 100%|█| 7000/7000 [1:06:45<00:00,  1.75it/s, 31 steps of size 1.20e-01. \n",
      "sample: 100%|█| 7000/7000 [1:07:15<00:00,  1.73it/s, 31 steps of size 1.20e-01. \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 6h 36min 37s, sys: 2h 17min 33s, total: 8h 54min 11s\n",
      "Wall time: 8h 52min 4s\n"
     ]
    }
   ],
   "source": [
    "sampler = infer.MCMC(kernel,\n",
    "    num_warmup=2000,\n",
    "    num_samples=5000,\n",
    "    num_chains=8,\n",
    "    progress_bar=True,\n",
    "    chain_method='sequential'\n",
    ")\n",
    "\n",
    "%time sampler.run(jax.random.PRNGKey(20), rate = 6, extra_fields = ['num_steps', 'accept_prob'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "ef23ab9f-bdd5-4790-a07c-3210a9ae6fc8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "                           mean       std    median      5.0%     95.0%     n_eff     r_hat\n",
      "auto_shared_latent[0]      0.07      1.05      0.08     -1.56      1.75  31233.41      1.00\n",
      "auto_shared_latent[1]     -0.08      1.00     -0.08     -1.75      1.49  33422.67      1.00\n",
      "auto_shared_latent[2]     -0.07      1.12     -0.08     -1.79      1.65  25083.05      1.00\n",
      "auto_shared_latent[3]     -0.07      1.01     -0.09     -1.65      1.60  35140.81      1.00\n",
      "auto_shared_latent[4]      0.04      1.08      0.05     -1.67      1.79  29003.16      1.00\n",
      "auto_shared_latent[5]      0.03      1.03      0.03     -1.63      1.70  35469.87      1.00\n",
      "auto_shared_latent[6]     -0.04      1.62     -0.02     -2.07      2.17   7851.00      1.00\n",
      "auto_shared_latent[7]      0.00      1.54     -0.01     -2.07      1.98   8837.32      1.00\n",
      "\n",
      "Number of divergences: 2\n"
     ]
    }
   ],
   "source": [
    "sampler.print_summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "5306d82c-6415-4f80-b4cb-7a0e735e7254",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'nrer_az/ntr_ddl_fixed_size_0.98_5k_8_1.json'"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "arviz_obj = az.from_numpyro(sampler)\n",
    "# az.to_json(arviz_obj, \"nrer_az/ntr_ddl_fixed_size_0.98_5k_8_1.json\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "9f166511-caef-49db-8ed1-450efa5942ec",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "sample: 100%|█| 7000/7000 [1:06:18<00:00,  1.76it/s, 31 steps of size 1.20e-01. \n",
      "sample: 100%|█| 7000/7000 [1:05:09<00:00,  1.79it/s, 31 steps of size 1.20e-01. \n",
      "sample: 100%|█| 7000/7000 [1:05:21<00:00,  1.79it/s, 31 steps of size 1.20e-01. \n",
      "sample: 100%|█| 7000/7000 [1:06:32<00:00,  1.75it/s, 31 steps of size 1.20e-01. \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 3h 18min 4s, sys: 1h 7min 5s, total: 4h 25min 9s\n",
      "Wall time: 4h 24min 3s\n",
      "\n",
      "                           mean       std    median      5.0%     95.0%     n_eff     r_hat\n",
      "auto_shared_latent[0]      0.06      1.05      0.07     -1.64      1.68  14010.25      1.00\n",
      "auto_shared_latent[1]     -0.08      1.01     -0.07     -1.75      1.53  15763.03      1.00\n",
      "auto_shared_latent[2]     -0.07      1.17     -0.06     -1.76      1.73   9603.74      1.00\n",
      "auto_shared_latent[3]     -0.08      1.01     -0.10     -1.69      1.61  17913.14      1.00\n",
      "auto_shared_latent[4]      0.05      1.08      0.06     -1.70      1.79  14873.02      1.00\n",
      "auto_shared_latent[5]      0.03      1.02      0.04     -1.63      1.69  16602.41      1.00\n",
      "auto_shared_latent[6]     -0.02      1.72      0.00     -2.07      2.11   2849.19      1.00\n",
      "auto_shared_latent[7]     -0.01      1.62     -0.02     -1.98      1.96   5103.52      1.00\n",
      "\n",
      "Number of divergences: 2\n"
     ]
    }
   ],
   "source": [
    "sampler = infer.MCMC(kernel,\n",
    "    num_warmup=2000,\n",
    "    num_samples=5000,\n",
    "    num_chains=4,\n",
    "    progress_bar=True,\n",
    "    chain_method='sequential'\n",
    ")\n",
    "\n",
    "%time sampler.run(jax.random.PRNGKey(18), rate = 6, extra_fields = ['num_steps', 'accept_prob'])\n",
    "sampler.print_summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "fdcffe4a-ae34-46e4-81e8-5c722e18f167",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'nrer_az/ntr_ddl_fixed_size_0.98_5k_4.json'"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "arviz_obj = az.from_numpyro(sampler)\n",
    "# az.to_json(arviz_obj, \"nrer_az/ntr_ddl_fixed_size_0.98_5k_4.json\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:numpyro_gpu]",
   "language": "python",
   "name": "conda-env-numpyro_gpu-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}