{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Copyright (C) 2024 Richard Stiskalek\n",
    "# This program is free software; you can redistribute it and/or modify it\n",
    "# under the terms of the GNU General Public License as published by the\n",
    "# Free Software Foundation; either version 3 of the License, or (at your\n",
    "# option) any later version.\n",
    "#\n",
    "# This program is distributed in the hope that it will be useful, but\n",
    "# WITHOUT ANY WARRANTY; without even the implied warranty of\n",
    "# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General\n",
    "# Public License for more details.\n",
    "#\n",
    "# You should have received a copy of the GNU General Public License along\n",
    "# with this program; if not, write to the Free Software Foundation, Inc.,\n",
    "# 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from tqdm import trange\n",
    "from joblib import dump\n",
    "from h5py import File\n",
    "\n",
    "import csiborgtools\n",
    "\n",
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "fpath = \"/mnt/extraspace/rstiskalek/catalogs/PV_compilation.hdf5\"\n",
    "\n",
    "with File(fpath, 'r') as f:\n",
    "    RA_2MTF = f[\"2MTF/RA\"][...]\n",
    "    DEC_2MTF = f[\"2MTF/DEC\"][...]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/shared/python/3.11.7/lib/python3.11/pty.py:89: RuntimeWarning: os.fork() was called. os.fork() is incompatible with multithreaded code, and JAX is multithreaded, so this will likely lead to a deadlock.\n",
      "  pid, fd = os.forkpty()\n",
      "  0%|          | 0/1 [00:00<?, ?it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Assigning calibration to 62/1248 galaxies.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1/1 [00:18<00:00, 18.57s/it]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saving to  /mnt/extraspace/rstiskalek/csiborg_postprocessing/flow_mock/Carrick2MTFmock_seed0.hdf5\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "!rm /mnt/extraspace/rstiskalek/csiborg_postprocessing/flow_mock/*\n",
    "\n",
    "\n",
    "reader = csiborgtools.read.Carrick2015Field()\n",
    "velocity = reader.velocity_field()\n",
    "boxsize = csiborgtools.simname2boxsize(\"Carrick2015\")\n",
    "\n",
    "mean_mag = 10.31\n",
    "std_mag = 0.83\n",
    "h = 0.7\n",
    "nrepeat_calibration = 1\n",
    "\n",
    "kmax = 1\n",
    "seed = 53\n",
    "\n",
    "for k in trange(kmax):\n",
    "    mock, truths = csiborgtools.flow.mock_Carrick2MTF(\n",
    "        velocity, boxsize, RA_2MTF, DEC_2MTF,\n",
    "        mean_mag=mean_mag, std_mag=std_mag, h=h,\n",
    "        nrepeat_calibration=nrepeat_calibration,\n",
    "        calibration_max_percentile=10, calibration_rand_fraction=1.,\n",
    "        seed=seed, verbose=True)\n",
    "    seed += 1\n",
    "\n",
    "    fname = f\"/mnt/extraspace/rstiskalek/csiborg_postprocessing/flow_mock/Carrick2MTFmock_seed{k}.hdf5\"  # noqa\n",
    "    print(\"Saving to \", fname)\n",
    "    with File(fname, 'w') as f:\n",
    "        for key in mock.keys():\n",
    "            f.create_dataset(key, data=mock[key])\n",
    "\n",
    "    fname_truths = fname.replace(\".hdf5\", \"_truths.pkl\")\n",
    "    dump(truths, fname_truths)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = np.isfinite(mock[\"mu_calibration\"][0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGdCAYAAACyzRGfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAip0lEQVR4nO3de3CU1eH/8c8GyAKS3RggCSkJd0HlokYIqYoIKRAd5eaIyIygDBYarJhaIY6K2O8YKlPxMoi2FdAqYukIVBmgECRoDSipKaIlJZkgWJKg0GQhyILJ+f3hjy1LwmWT3bPZ+H7N7AzZ5+yz5+wTdt+z2YvDGGMEAABgSVS4JwAAAH5ciA8AAGAV8QEAAKwiPgAAgFXEBwAAsIr4AAAAVhEfAADAKuIDAABY1TrcEzhXXV2dDh06pJiYGDkcjnBPBwAAXAJjjI4dO6akpCRFRV34uY1mFx+HDh1ScnJyuKcBAAAa4eDBg+ratesFxzS7+IiJiZH0w+RdLleYZwMAAC6Fx+NRcnKy73H8QppdfJz5U4vL5SI+AACIMJfykglecAoAAKwiPgAAgFXEBwAAsIr4AAAAVhEfAADAKuIDAABYRXwAAACriA8AAGAV8QEAAKwiPgAAgFXEBwAAsIr4AAAAVhEfAADAKuIDAABY1TrcEwAAIJJ1n7c+3FMI2P6Ft4X1+nnmAwAAWEV8AAAAq4gPAABgFfEBAACsIj4AAIBVxAcAALCK+AAAAFYRHwAAwCriAwAAWEV8AAAAq4gPAABgFfEBAACsIj4AAIBVxAcAALCK+AAAAFYRHwAAwCriAwAAWBVQfCxdulQDBw6Uy+WSy+VSenq6NmzY4Nt+8uRJZWVlqWPHjurQoYMmTpyoysrKoE8aAABEroDio2vXrlq4cKEKCwu1a9cujRgxQmPHjtUXX3whSXr44Yf13nvvafXq1crPz9ehQ4c0YcKEkEwcAABEJocxxjRlB3FxcVq0aJHuvPNOde7cWStXrtSdd94pSdq7d6+uvPJKFRQUaOjQoZe0P4/HI7fbrerqarlcrqZMDQCAkOs+b324pxCw/QtvC/o+A3n8bvRrPmpra7Vq1SrV1NQoPT1dhYWFOn36tDIyMnxj+vXrp5SUFBUUFJx3P16vVx6Px+8EAABaroDj4/PPP1eHDh3kdDo1c+ZMrVmzRldddZUqKioUHR2t2NhYv/EJCQmqqKg47/5yc3Pldrt9p+Tk5IAXAQAAIkfA8dG3b18VFRVp586dmjVrlqZOnaovv/yy0RPIyclRdXW173Tw4MFG7wsAADR/rQO9QHR0tHr37i1JSk1N1aeffqoXXnhBkyZN0qlTp1RVVeX37EdlZaUSExPPuz+n0ymn0xn4zAEAQERq8ud81NXVyev1KjU1VW3atFFeXp5vW3FxsQ4cOKD09PSmXg0AAGghAnrmIycnR5mZmUpJSdGxY8e0cuVKbdu2TZs2bZLb7db06dOVnZ2tuLg4uVwuPfjgg0pPT7/kd7oAAICWL6D4OHz4sO69916Vl5fL7XZr4MCB2rRpk372s59JkhYvXqyoqChNnDhRXq9Xo0eP1ssvvxySiQMAgMjU5M/5CDY+5wMAEEn4nI8fWPmcDwAAgMYgPgAAgFXEBwAAsIr4AAAAVgX8IWMAAIRKJL54E4HjmQ8AAGAV8QEAAKwiPgAAgFXEBwAAsIr4AAAAVhEfAADAKuIDAABYRXwAAACriA8AAGAV8QEAAKwiPgAAgFXEBwAAsIr4AAAAVhEfAADAKuIDAABYRXwAAACriA8AAGAV8QEAAKwiPgAAgFXEBwAAsIr4AAAAVhEfAADAKuIDAABYRXwAAACriA8AAGAV8QEAAKwiPgAAgFXEBwAAsIr4AAAAVhEfAADAKuIDAABYRXwAAACriA8AAGAV8QEAAKwiPgAAgFXEBwAAsIr4AAAAVhEfAADAKuIDAABYFVB85ObmavDgwYqJiVF8fLzGjRun4uJivzHDhw+Xw+HwO82cOTOokwYAAJEroPjIz89XVlaWduzYoc2bN+v06dMaNWqUampq/MbNmDFD5eXlvtOzzz4b1EkDAIDI1TqQwRs3bvT7ecWKFYqPj1dhYaGGDRvmO799+/ZKTEwMzgwBAECL0qTXfFRXV0uS4uLi/M5/66231KlTJ/Xv3185OTk6ceLEeffh9Xrl8Xj8TgAAoOUK6JmPs9XV1WnOnDm64YYb1L9/f9/599xzj7p166akpCTt3r1bc+fOVXFxsd59990G95Obm6sFCxY0dhoAgPPoPm99uKcANMhhjDGNueCsWbO0YcMGffTRR+ratet5x23dulUjR45USUmJevXqVW+71+uV1+v1/ezxeJScnKzq6mq5XK7GTA0AIOID57d/4W1B36fH45Hb7b6kx+9GPfMxe/Zsvf/++9q+ffsFw0OS0tLSJOm88eF0OuV0OhszDQAAEIECig9jjB588EGtWbNG27ZtU48ePS56maKiIklSly5dGjVBAADQsgQUH1lZWVq5cqXWrVunmJgYVVRUSJLcbrfatWun0tJSrVy5Urfeeqs6duyo3bt36+GHH9awYcM0cODAkCwAAABEloDiY+nSpZJ++CCxsy1fvlzTpk1TdHS0tmzZoueff141NTVKTk7WxIkT9fjjjwdtwgAAILIF/GeXC0lOTlZ+fn6TJgQAAFo2vtsFAABYRXwAAACriA8AAGAV8QEAAKwiPgAAgFXEBwAAsIr4AAAAVhEfAADAKuIDAABYRXwAAACriA8AAGAV8QEAAKwiPgAAgFXEBwAAsIr4AAAAVhEfAADAKuIDAABYRXwAAACriA8AAGAV8QEAAKwiPgAAgFXEBwAAsIr4AAAAVhEfAADAKuIDAABYRXwAAACriA8AAGAV8QEAAKwiPgAAgFXEBwAAsIr4AAAAVhEfAADAKuIDAABYRXwAAACriA8AAGAV8QEAAKwiPgAAgFXEBwAAsIr4AAAAVhEfAADAKuIDAABYRXwAAACriA8AAGAV8QEAAKwKKD5yc3M1ePBgxcTEKD4+XuPGjVNxcbHfmJMnTyorK0sdO3ZUhw4dNHHiRFVWVgZ10gAAIHIFFB/5+fnKysrSjh07tHnzZp0+fVqjRo1STU2Nb8zDDz+s9957T6tXr1Z+fr4OHTqkCRMmBH3iAAAgMrUOZPDGjRv9fl6xYoXi4+NVWFioYcOGqbq6Wq+99ppWrlypESNGSJKWL1+uK6+8Ujt27NDQoUODN3MAABCRmvSaj+rqaklSXFycJKmwsFCnT59WRkaGb0y/fv2UkpKigoKCBvfh9Xrl8Xj8TgAAoOVqdHzU1dVpzpw5uuGGG9S/f39JUkVFhaKjoxUbG+s3NiEhQRUVFQ3uJzc3V26323dKTk5u7JQAAEAEaHR8ZGVlac+ePVq1alWTJpCTk6Pq6mrf6eDBg03aHwAAaN4Ces3HGbNnz9b777+v7du3q2vXrr7zExMTderUKVVVVfk9+1FZWanExMQG9+V0OuV0OhszDQAAEIECeubDGKPZs2drzZo12rp1q3r06OG3PTU1VW3atFFeXp7vvOLiYh04cEDp6enBmTEAAIhoAT3zkZWVpZUrV2rdunWKiYnxvY7D7XarXbt2crvdmj59urKzsxUXFyeXy6UHH3xQ6enpvNMFAABICjA+li5dKkkaPny43/nLly/XtGnTJEmLFy9WVFSUJk6cKK/Xq9GjR+vll18OymQBAEDkCyg+jDEXHdO2bVstWbJES5YsafSkAABAy8V3uwAAAKuIDwAAYBXxAQAArCI+AACAVcQHAACwivgAAABWER8AAMAq4gMAAFhFfAAAAKuIDwAAYBXxAQAArCI+AACAVcQHAACwivgAAABWER8AAMAq4gMAAFhFfAAAAKuIDwAAYBXxAQAArCI+AACAVcQHAACwivgAAABWER8AAMAq4gMAAFhFfAAAAKuIDwAAYFXrcE8AwI9L93nrwz2FgO1feFu4pwC0KDzzAQAArCI+AACAVcQHAACwivgAAABWER8AAMAq4gMAAFhFfAAAAKuIDwAAYBXxAQAArCI+AACAVcQHAACwivgAAABWER8AAMAq4gMAAFhFfAAAAKuIDwAAYBXxAQAArAo4PrZv367bb79dSUlJcjgcWrt2rd/2adOmyeFw+J3GjBkTrPkCAIAIF3B81NTUaNCgQVqyZMl5x4wZM0bl5eW+09tvv92kSQIAgJajdaAXyMzMVGZm5gXHOJ1OJSYmNnpSAACg5QrJaz62bdum+Ph49e3bV7NmzdKRI0fOO9br9crj8fidAABAyxX0+BgzZozeeOMN5eXl6be//a3y8/OVmZmp2traBsfn5ubK7Xb7TsnJycGeEgAAaEYC/rPLxdx9992+fw8YMEADBw5Ur169tG3bNo0cObLe+JycHGVnZ/t+9ng8BAgAAC1YyN9q27NnT3Xq1EklJSUNbnc6nXK5XH4nAADQcoU8Pr7++msdOXJEXbp0CfVVAQCACBDwn12OHz/u9yxGWVmZioqKFBcXp7i4OC1YsEATJ05UYmKiSktL9eijj6p3794aPXp0UCcOAAAiU8DxsWvXLt1yyy2+n8+8XmPq1KlaunSpdu/erddff11VVVVKSkrSqFGj9Jvf/EZOpzN4swYAABEr4PgYPny4jDHn3b5p06YmTQgAALRsfLcLAACwivgAAABWER8AAMAq4gMAAFhFfAAAAKuIDwAAYBXxAQAArCI+AACAVcQHAACwivgAAABWER8AAMAq4gMAAFhFfAAAAKuIDwAAYBXxAQAArCI+AACAVcQHAACwivgAAABWER8AAMAq4gMAAFhFfAAAAKuIDwAAYBXxAQAArCI+AACAVcQHAACwivgAAABWtQ73BACgues+b324pwC0KDzzAQAArCI+AACAVcQHAACwivgAAABWER8AAMAq4gMAAFhFfAAAAKuIDwAAYBXxAQAArCI+AACAVcQHAACwivgAAABWER8AAMAq4gMAAFhFfAAAAKuIDwAAYBXxAQAArAo4PrZv367bb79dSUlJcjgcWrt2rd92Y4yefPJJdenSRe3atVNGRob27dsXrPkCAIAIF3B81NTUaNCgQVqyZEmD25999lm9+OKLeuWVV7Rz505ddtllGj16tE6ePNnkyQIAgMjXOtALZGZmKjMzs8Ftxhg9//zzevzxxzV27FhJ0htvvKGEhAStXbtWd999d9NmCwAAIl5QX/NRVlamiooKZWRk+M5zu91KS0tTQUFBg5fxer3yeDx+JwAA0HIFNT4qKiokSQkJCX7nJyQk+LadKzc3V26323dKTk4O5pQAAEAzE/Z3u+Tk5Ki6utp3OnjwYLinBAAAQiio8ZGYmChJqqys9Du/srLSt+1cTqdTLpfL7wQAAFquoMZHjx49lJiYqLy8PN95Ho9HO3fuVHp6ejCvCgAARKiA3+1y/PhxlZSU+H4uKytTUVGR4uLilJKSojlz5uj//u//1KdPH/Xo0UNPPPGEkpKSNG7cuGDOGwAARKiA42PXrl265ZZbfD9nZ2dLkqZOnaoVK1bo0UcfVU1NjR544AFVVVXpxhtv1MaNG9W2bdvgzRoAAEQshzHGhHsSZ/N4PHK73aqurub1H0AL1H3e+nBPAfjR27/wtqDvM5DH77C/2wUAAPy4EB8AAMAq4gMAAFhFfAAAAKuIDwAAYBXxAQAArCI+AACAVcQHAACwivgAAABWER8AAMAq4gMAAFhFfAAAAKuIDwAAYBXxAQAArCI+AACAVcQHAACwivgAAABWER8AAMAq4gMAAFhFfAAAAKuIDwAAYBXxAQAArCI+AACAVcQHAACwivgAAABWER8AAMAq4gMAAFhFfAAAAKuIDwAAYBXxAQAArCI+AACAVcQHAACwivgAAABWER8AAMAq4gMAAFhFfAAAAKuIDwAAYBXxAQAArCI+AACAVcQHAACwivgAAABWER8AAMAq4gMAAFhFfAAAAKuCHh9PPfWUHA6H36lfv37BvhoAABChWodip1dffbW2bNnyvytpHZKrAQAAESgkVdC6dWslJiaGYtcAACDCheQ1H/v27VNSUpJ69uypKVOm6MCBA+cd6/V65fF4/E4AAKDlCvozH2lpaVqxYoX69u2r8vJyLViwQDfddJP27NmjmJiYeuNzc3O1YMGCYE8D+FHoPm99uKcAAAFzGGNMKK+gqqpK3bp103PPPafp06fX2+71euX1en0/ezweJScnq7q6Wi6XK5RTAyIe8QGgMfYvvC3o+/R4PHK73Zf0+B3yV4LGxsbqiiuuUElJSYPbnU6nnE5nqKcBAACaiZB/zsfx48dVWlqqLl26hPqqAABABAh6fDzyyCPKz8/X/v379fHHH2v8+PFq1aqVJk+eHOyrAgAAESjof3b5+uuvNXnyZB05ckSdO3fWjTfeqB07dqhz587BvioAABCBgh4fq1atCvYuAQBAC8J3uwAAAKuIDwAAYBXxAQAArCI+AACAVcQHAACwivgAAABWER8AAMAq4gMAAFhFfAAAAKuIDwAAYBXxAQAArCI+AACAVcQHAACwivgAAABWER8AAMCq1uGegG3d560P9xQCtn/hbeGeAgAAQcMzHwAAwCriAwAAWEV8AAAAq4gPAABgFfEBAACsIj4AAIBVxAcAALCK+AAAAFYRHwAAwCriAwAAWEV8AAAAq4gPAABgFfEBAACsIj4AAIBVxAcAALCK+AAAAFYRHwAAwCriAwAAWEV8AAAAq4gPAABgFfEBAACsah3uCaBl6j5vfbinAABopnjmAwAAWEV8AAAAq4gPAABgFfEBAACsIj4AAIBVxAcAALAqZPGxZMkSde/eXW3btlVaWpo++eSTUF0VAACIICGJj3feeUfZ2dmaP3++/vGPf2jQoEEaPXq0Dh8+HIqrAwAAESQk8fHcc89pxowZuu+++3TVVVfplVdeUfv27bVs2bJQXB0AAIggQf+E01OnTqmwsFA5OTm+86KiopSRkaGCgoJ6471er7xer+/n6upqSZLH4wn21CRJdd4TIdlvKIXqtgilSLydAeDHIhSPK2f2aYy56Nigx8e3336r2tpaJSQk+J2fkJCgvXv31hufm5urBQsW1Ds/OTk52FOLWO7nwz0DAEBLEsrHlWPHjsntdl9wTNi/2yUnJ0fZ2dm+n+vq6nT06FF17NhRDofjkvbh8XiUnJysgwcPyuVyhWqqYdXS18j6Ihvri3wtfY2sL/SMMTp27JiSkpIuOjbo8dGpUye1atVKlZWVfudXVlYqMTGx3nin0ymn0+l3XmxsbKOu2+VytchfqrO19DWyvsjG+iJfS18j6wutiz3jcUbQX3AaHR2t1NRU5eXl+c6rq6tTXl6e0tPTg311AAAgwoTkzy7Z2dmaOnWqrr/+eg0ZMkTPP/+8ampqdN9994Xi6gAAQAQJSXxMmjRJ33zzjZ588klVVFTommuu0caNG+u9CDVYnE6n5s+fX+/PNy1JS18j64tsrC/ytfQ1sr7mxWEu5T0xAAAAQcJ3uwAAAKuIDwAAYBXxAQAArCI+AACAVWGPj9zcXA0ePFgxMTGKj4/XuHHjVFxc7DemtLRU48ePV+fOneVyuXTXXXfV+xCzc3Xv3l0Oh6PeKSsryzdm+PDh9bbPnDkzqOtbunSpBg4c6Pvgl/T0dG3YsMG3/eTJk8rKylLHjh3VoUMHTZw48aJrM8boySefVJcuXdSuXTtlZGRo3759fmOOHj2qKVOmyOVyKTY2VtOnT9fx48eDurZQrO/06dOaO3euBgwYoMsuu0xJSUm69957dejQIb9xDR3fhQsXNvv1SdK0adPqzX3MmDF+Y2wdv1CtsaH/ew6HQ4sWLfKNaS7H8Pe//72GDx8ul8slh8OhqqqqS9rvkiVL1L17d7Vt21ZpaWn65JNP/LY35nZrjFCs71Lul23cf0qhWd9TTz1Vb+79+vXzGxPJx6+5PP5dkAmz0aNHm+XLl5s9e/aYoqIic+utt5qUlBRz/PhxY4wxx48fNz179jTjx483u3fvNrt37zZjx441gwcPNrW1tefd7+HDh015ebnvtHnzZiPJfPDBB74xN998s5kxY4bfuOrq6qCu769//atZv369+fe//22Ki4vNY489Ztq0aWP27NljjDFm5syZJjk52eTl5Zldu3aZoUOHmp/+9KcX3OfChQuN2+02a9euNf/85z/NHXfcYXr06GG+++4735gxY8aYQYMGmR07dpgPP/zQ9O7d20yePDmoawvF+qqqqkxGRoZ55513zN69e01BQYEZMmSISU1N9RvXrVs38/TTT/sduzO/M815fcYYM3XqVDNmzBi/uR89etRvjK3jF6o1nr228vJys2zZMuNwOExpaalvTHM5hosXLza5ubkmNzfXSDL//e9/L7rPVatWmejoaLNs2TLzxRdfmBkzZpjY2FhTWVnpG9OY2625rO9i98vG2Ln/DNX65s+fb66++mq/uX/zzTd+YyL5+DWXx78LCXt8nOvw4cNGksnPzzfGGLNp0yYTFRXld6NUVVUZh8NhNm/efMn7feihh0yvXr1MXV2d77ybb77ZPPTQQ0Gb+6W6/PLLzR//+EdTVVVl2rRpY1avXu3b9q9//ctIMgUFBQ1etq6uziQmJppFixb5zquqqjJOp9O8/fbbxhhjvvzySyPJfPrpp74xGzZsMA6Hw/znP/8J0ar+pynra8gnn3xiJJmvvvrKd163bt3M4sWLgzntS9bU9U2dOtWMHTv2vNvDffyMCf4xHDt2rBkxYoTfec3hGJ7tgw8+uOQ79yFDhpisrCzfz7W1tSYpKcnk5uYaY0zQbrfGaur6znXu/bIx4bv/NKbp65s/f74ZNGjQebe3tOPXnB7/zgj7n13OVV1dLUmKi4uTJHm9XjkcDr8PTmnbtq2ioqL00UcfXdI+T506pTfffFP3339/vS+re+utt9SpUyf1799fOTk5OnEidF8FX1tbq1WrVqmmpkbp6ekqLCzU6dOnlZGR4RvTr18/paSkqKCgoMF9lJWVqaKiwu8ybrdbaWlpvssUFBQoNjZW119/vW9MRkaGoqKitHPnzhCtLjjra0h1dbUcDke97/xZuHChOnbsqGuvvVaLFi3S999/H6ylNCiY69u2bZvi4+PVt29fzZo1S0eOHPFtC9fxk0JzDCsrK7V+/XpNnz693rZwH8PGOHXqlAoLC/1uk6ioKGVkZPhuk2D97gcqGOtryLn3y2fYvP+Ugru+ffv2KSkpST179tSUKVN04MAB37aWdPyay+PfucL+rbZnq6ur05w5c3TDDTeof//+kqShQ4fqsssu09y5c/XMM8/IGKN58+aptrZW5eXll7TftWvXqqqqStOmTfM7/5577lG3bt2UlJSk3bt3a+7cuSouLta7774b1HV9/vnnSk9P18mTJ9WhQwetWbNGV111lYqKihQdHV3vQTUhIUEVFRUN7uvM+ed+WuzZl6moqFB8fLzf9tatWysuLu68+22KYK7vXCdPntTcuXM1efJkvy9L+uUvf6nrrrtOcXFx+vjjj5WTk6Py8nI999xzwVyapOCvb8yYMZowYYJ69Oih0tJSPfbYY8rMzFRBQYFatWpl/fhJoT2Gr7/+umJiYjRhwgS/85vDMWyMb7/9VrW1tQ3+H9y7d6+kH/4PNvV2C0Qw13euhu6XJXv3n1Lw15eWlqYVK1aob9++Ki8v14IFC3TTTTdpz549iomJaVHHL9yPf+fTrOIjKytLe/bs8XtGo3Pnzlq9erVmzZqlF198UVFRUZo8ebKuu+46RUVd2hM3r732mjIzM+t9ze8DDzzg+/eAAQPUpUsXjRw5UqWlperVq1dwFiWpb9++KioqUnV1tf7yl79o6tSpys/PD9r+wy1U6zt9+rTuuusuGWO0dOlSv23Z2dm+fw8cOFDR0dH6+c9/rtzc3KB/vHCw13f33Xf7/j1gwAANHDhQvXr10rZt2zRy5MhgTDlgofwdXbZsmaZMmaK2bdv6nd8cjmGw7uDDLZTra+h+WbJ3/ykFf32ZmZm+fw8cOFBpaWnq1q2b/vznPzf4DF2ohfL4hfvx73yaTXzMnj1b77//vrZv366uXbv6bRs1apRKS0v17bffqnXr1oqNjVViYqJ69ux50f1+9dVX2rJlyyXVXFpamiSppKQkqDd+dHS0evfuLUlKTU3Vp59+qhdeeEGTJk3SqVOnVFVV5VfYlZWVSkxMbHBfZ86vrKxUly5d/C5zzTXX+MYcPnzY73Lff/+9jh49et79NkUw13fGmfD46quvtHXr1ot+RXRaWpq+//577d+/X3379m3yms4WivWdrWfPnurUqZNKSko0cuRI68dPCt0aP/zwQxUXF+udd9656NhwHMNXX3014H116tRJrVq1qvfOh7Nvk8TExKD8blyqYK7vbBe6Xz5XqO4/pdCt74zY2FhdccUVKikpkdRyjl9zePw7n7C/5sMYo9mzZ2vNmjXaunWrevTocd6xnTp1UmxsrLZu3arDhw/rjjvuuOj+ly9frvj4eN12220XHVtUVCRJfg/qoVBXVyev16vU1FS1adNGeXl5vm3FxcU6cODAef/e16NHDyUmJvpdxuPxaOfOnb7LpKenq6qqSoWFhb4xW7duVV1dne8XLJSasj7pf+Gxb98+bdmyRR07drzodRYVFSkqKqrenytCoanrO9fXX3+tI0eO+H7vwn38pOCt8bXXXlNqaqoGDRp00bHhOIaNER0drdTUVL/bpK6uTnl5eb7bJFi/G43VlPVJgd0vn2Hr/lNq+vrOdfz4cZWWlvrmHunH74zm+PjnE7aXuv5/s2bNMm6322zbts3vLT8nTpzwjVm2bJkpKCgwJSUl5k9/+pOJi4sz2dnZfvsZMWKEeemll/zOq62tNSkpKWbu3Ln1rrekpMQ8/fTTZteuXaasrMysW7fO9OzZ0wwbNiyo65s3b57Jz883ZWVlZvfu3WbevHnG4XCYv/3tb8aYH97OlZKSYrZu3Wp27dpl0tPTTXp6ut8++vbta959913fzwsXLjSxsbFm3bp1vrceN/RW22uvvdbs3LnTfPTRR6ZPnz4heatmsNd36tQpc8cdd5iuXbuaoqIiv98Jr9drjDHm448/NosXLzZFRUWmtLTUvPnmm6Zz587m3nvvbfbrO3bsmHnkkUdMQUGBKSsrM1u2bDHXXXed6dOnjzl58qTvMraOXyjWeEZ1dbVp3769Wbp0ab3rbE7HsLy83Hz22WfmD3/4g5Fktm/fbj777DNz5MgR3z7OvX9ZtWqVcTqdZsWKFebLL780DzzwgImNjTUVFRW+MZdyuzXX9V3sftnW/Weo1verX/3KbNu2zZSVlZm///3vJiMjw3Tq1MkcPnzYNyaSj58xzePx70LCHh+SGjwtX77cN2bu3LkmISHBtGnTxvTp08f87ne/83vLkDE/vG1v/vz5fudt2rTJSDLFxcX1rvfAgQNm2LBhJi4uzjidTtO7d2/z61//Oujvc77//vtNt27dTHR0tOncubMZOXKk75fKGGO+++4784tf/MJcfvnlpn379mb8+PGmvLzcbx/n3h51dXXmiSeeMAkJCcbpdJqRI0fWW+ORI0fM5MmTTYcOHYzL5TL33XefOXbsWFDXFor1lZWVnfd34sx71AsLC01aWppxu92mbdu25sorrzTPPPOM34N3c13fiRMnzKhRo0znzp1NmzZtTLdu3cyMGTP8HrSMsXf8QrHGM1599VXTrl07U1VVVe86m9MxnD9//kXvgxq6f3nppZdMSkqKiY6ONkOGDDE7duzw234pt1tzXd/F7pdt3X+Gan2TJk0yXbp0MdHR0eYnP/mJmTRpkikpKfG73kg+fsY0j8e/C3EYY0zwn08BAABoWNhf8wEAAH5ciA8AAGAV8QEAAKwiPgAAgFXEBwAAsIr4AAAAVhEfAADAKuIDAABYRXwAAACriA8AAGAV8QEAAKwiPgAAgFX/D8kkk+pn4/lwAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.hist(mock[\"mu_calibration\"][0, m], bins=\"auto\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAGdCAYAAADqsoKGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAABNfElEQVR4nO3df1iUZb4/8PeAzKDCDIHiwAqKP0IRsaTEOZWZQqIdM8XdMv2q5Wq64CZuu0ZnW7M9LbRe3zV32+jHMWu/RuzakVzriMdf0GaoiLGCFKssauUMrBIDQgzI3N8/2JllZAbmmRlgeHi/rmuunOe5n3vum9GLT/ePz60QQggQERERDXA+/d0AIiIiIk9gUENERESywKCGiIiIZIFBDREREckCgxoiIiKSBQY1REREJAsMaoiIiEgWGNQQERGRLAzp7wb0FbPZjKtXryIwMBAKhaK/m0NEREROEEKgsbER4eHh8PHpfixm0AQ1V69eRURERH83g4iIiFzw1VdfYfTo0d2WGTRBTWBgIICOH4pare7n1hAREZEzGhoaEBERYf093p1BE9RYppzUajWDGiIiogHGmaUjXChMREREssCghoiIiGSBQQ0RERHJAoMaIiIikgUGNURERCQLDGqIiIhIFhjUEBERkSwwqCEiIiJZGDTJ94iIiMh17WaB09V1qG1sQWigP2ZEBcPXx7vOUmRQQ0RERN3KL9dj24EK6I0t1mvBw5V45I5wJMVovSbAUQghRH83oi80NDRAo9HAaDTymAQiIiIn5ZfrsWHPWXQXLIRp/LF1YQySY8M8/vlSfn9zTQ0RERHZ1W4W2HagotuABgD0xhas33MWvzxwHkVV19Fu7p/xEk4/ERERkV2nq+tsppx6suvEJew6calXR266w5EaIiIisqu20fmApjODsQUb9pxFfrnewy3qHoMaIiIisis00N+l5yyTT9sOVPTpVBSDGiIiIrJrRlQwwjT+cGVfk0DHWpvT1XWebpZDDGqIiIgGmXazQFHVdewv/abbhb2+PgpsXRgDAC4FNoDrU1iukBTUZGdnIy4uDmq1Gmq1GjqdDgcPHrTef/PNNzF79myo1WooFArU19f3WOfYsWOhUCi6vFJTU61lZs+e3eX++vXrpTSdiIiI0LFF+96Xj2HZWyfxdG4plr11Eve+fMzh+pfk2DBkr5gOzTA/lz7P1SksV0gKakaPHo2srCyUlJTgzJkzmDNnDhYtWoTz588DAJqbm5GcnIznnnvO6TqLi4uh1+utr8OHDwMAvv/979uUW7t2rU25X//611KaTkRENOhZcs7cuqNJ78TC3vrmNkmfpUBH/poZUcGuNNUlkrZ0L1y40Ob9Sy+9hOzsbJw8eRJTpkzBpk2bAAAFBQVO1zly5Eib91lZWRg/fjzuv/9+m+vDhg2DVquV0lwiIiL6p55yzggAGfvKkBSjtckObHlOCsvTWxfG9GmmYZfX1LS3tyM3NxdNTU3Q6XQeaUxrayv27NmDJ598EgqF7Q/hvffew4gRIxAbG4uMjAw0Nzd3W5fJZEJDQ4PNi4iIaLByJufMt81tePXYBcnP3Uqr8Uf2iul9nqdGcvK9srIy6HQ6tLS0ICAgAHl5eYiJifFIYz788EPU19dj9erVNtcff/xxjBkzBuHh4Th37hy2bNmCyspK7Nu3z2FdmZmZ2LZtm0faRURENNA5u2D37U+rkTZnonWExdnn0h4Yj4mjAvv1sEvJQU10dDRKS0thNBrxwQcfYNWqVSgsLPRIYLNr1y7Mnz8f4eHhNtfXrVtn/fPUqVMRFhaGuXPnoqqqCuPHj7dbV0ZGBjZv3mx939DQgIiICLfbSERENBA5u2DX2HITrx67iKcTJ0p67p4JI6EbH+Jy+zxBclCjVCoxYcIEAEB8fDyKi4uxc+dOvPHGG2415PLlyzhy5Ei3oy8WCQkJAICLFy86DGpUKhVUKpVbbSIiIhpI2s0Cp6vrUNvY0mXEZEZUMIKG+qH+u54X/O448jcAQNqcCdZcNd1NQfX1gmBH3D77yWw2w2Qyud2Q3bt3IzQ0FA899FCPZUtLSwEAYWF9O1dHRETkrfLL9dh2oMIm+Oh8BpOvjwJP3BNlDVh6suPI3/D+6ct44eEpeHhaGN74pNph2YenhfXLdNOtJAU1GRkZmD9/PiIjI9HY2IicnBwUFBTg0KFDAACDwQCDwYCLFy8C6Fh/ExgYiMjISAQHd0Rwc+fOxeLFi5GWlmat12w2Y/fu3Vi1ahWGDLFtUlVVFXJycrBgwQKEhITg3LlzSE9Px6xZsxAXF+dW54mIiOTAslX71p1NljOYsldMR1KMFmYh7cgCQ4MJ6/ecxXClb7fl3jv1FSZp1dBqhvbbehpAYlBTW1uLlStXQq/XQ6PRIC4uDocOHUJSUhIA4PXXX7dZnDtr1iwAHaMwlsW/VVVVuHbtmk29R44cwZUrV/Dkk092+UylUokjR47glVdeQVNTEyIiIpCSkoKf//znkjpKREQkR91t1Rbo2F6dsa8Mz/73OdR/d9Olz2hqbe/2/g3TTaT/6a8A0G8ndAOAQgiJYdsA1dDQAI1GA6PRCLVa3d/NISIi8oiiqutY9tbJ/m6GlWWMxlNbuqX8/ubZT0RERAOYwfhdfzfBRn+d0A0wqCEiIhqw8sv1+OXHX/R3M7rojxO6AQ/sfiIiIqK+52hxsKdZppNc+Zy+PKEb4EgNERHRgNPTOU6eYglo1s2Kcun5vjyhG2BQQ0RENOA4ex6Tws2d1ZYznDIWxOD1FdOhGercBE9/nNANMKghIiIacJyd1nF3f/N3rf/aAp4Uo8VQv56Dmv46oRtgUENERDTg9NW0Tv13N7F+z1nkl+txuroOhoaeg6ng4cp+OaEbYFBDRETkldrNAkVV17G/9BsUVV232R5tOY+pr8ZBtu4vx6cX/+FU2Z8/NLlfAhqAu5+IiIi8jjPnOG1dGIMNe85CAdd2JklR09iK3x+vcqqsVjO0l1vjGEdqiIiI+kh3oy8Wlq3aty4EtpzjlF+uBwAkx4Yhe8V0aDV9u8PIkf5aHNwZR2qIiIj6QE+jL4Bz5zhtO1CBpBgtfH0USI4NQ1KMFqer6/B6YRUK/+bcFJGn9efi4M44UkNERNTLnB196Wmrtr1Mvb4+CsyICkbZN0an2/PUrCiEeXCEx7L1u7/W0lhwpIaIiKgXSRl9cXar9tX677DrL3/H5bpmjAkehttHBaKuqbXH5xQAfv/4dCyIC8PPkifjdHUdjlQYsOvEJQk9spX2wASkJ93eryM0FgxqiIiIepGU0ZdL15qdqvOZD/5qk4PG2XBi9b+NwYK4jtEUXx8FdONDoBsfgrujgvHsvjLUN7c5WdO/3DNhhFcENACDGiIiol7l7OiLwfgd3j99xamytybVc3b3U9Awpd3rlrU5J/9+HUVV1wEIJIwNwU//+xxqGlrs1q9Ax7RTfy4MvhXX1BAREfUiZxPlXbvR6lRyO3e8f/qK3R1XQMfIzT0TRuCZedF4Zt4k3Bc9Ei88HAOg60iQtywMvhWDGiIiol7SbhYwC4GgoX4OyygABA3zw6vHL/Z6ewwNJptFxj1xtG3cWxYG34rTT0RERL3A3hbuW1kS57mylsVVzk6HWXTeNl7b2ILQwI4pJ28aobFgUENERORhli3cPa11GaVWoeWmuU+DGlfOjbIsKvZ2nH4iIiLyoO62cFsMU/riPxZMxval0yQFNO6MjXhDxt/exqCGiIjIg3rawg0Aza3teOl/vsDG9z+XVHdiTGi395NiQqHAwFnY62kMaoiIiDxIypqV+u+kjdKUf9OAtfdF4da4xEfRkSX4rZV3D6iFvZ7GNTVEREQe5MqaFWdYkvTNmTQKP503Cf+v6JI1o/D/0Y2FckjHOMVAWtjraQxqiIiIPGhGVDDCNP4wGO0nrXNXbWMLlEN8sOa+cQ7LDJSFvZ7G6SciIiIP8vVRYOtC+0nrPKG3RoLkgEENERGRhzlKWueu4SpfWe9echeDGiIiol6QHBuGT7fMQers8R6rs8nUjsMVBo/VJzcMaoiIiHqJr48C904c6dE6tx2ocHh+02DHoIaIiKgXzYgKRvBw+6dju0JvbJF0ftNgwqCGiIioG+1mgaKq69hf+g2Kqq5LHiXx9VHgPxfFerRNUs9vGiy4pZuIiMgBe4dShmn8sXVhDJJjw9BuFk7lg1kQF4anvo7CG59Ue6Rd3AFln6SRmuzsbMTFxUGtVkOtVkOn0+HgwYPW+2+++SZmz54NtVoNhUKB+vr6Hut84YUXoFAobF6TJk2yKdPS0oLU1FSEhIQgICAAKSkpqKmpkdJ0IiIiSSyHUt565IHB2IINe84iLecs4n95GMveOomnc0ux7K2TuPflY8gv19ut787I29xu02A4v8kdkoKa0aNHIysrCyUlJThz5gzmzJmDRYsW4fz58wCA5uZmJCcn47nnnpPUiClTpkCv11tfn376qc399PR0HDhwAHv37kVhYSGuXr2KJUuWSPoMIiIiZ3V3KKX45+ujc/ouxxxYAp5bAxtLfZ4g9/Ob3CFp+mnhwoU271966SVkZ2fj5MmTmDJlCjZt2gQAKCgokNaIIUOg1Wrt3jMajdi1axdycnIwZ84cAMDu3bsxefJknDx5EjNnzpT0WURERD1x5lBKewQ6RlO2HahAUozWGny4Wl9nwcP98KvFU2V/fpM7XF4o3N7ejtzcXDQ1NUGn07nViAsXLiA8PBzjxo3D8uXLceXKFeu9kpIStLW1ITEx0Xpt0qRJiIyMRFFRkcM6TSYTGhoabF5ERETOcGchruWMps47lI64mVtG7T8EJzMSGdD0QHJQU1ZWhoCAAKhUKqxfvx55eXmIiYlxuQEJCQl45513kJ+fj+zsbFRXV+O+++5DY2MjAMBgMECpVCIoKMjmuVGjRsFgcPyXJDMzExqNxvqKiIhwuY1ERDS4eGIhriUwyi/XY9eJS27V9f340dYDK8kxyT+h6OholJaW4tSpU9iwYQNWrVqFigrX5wnnz5+P73//+4iLi8O8efPwP//zP6ivr8ef/vQnl+sEgIyMDBiNRuvrq6++cqs+IiIaPCyHUrqzcuXStWaPraVJjLG/RINsSd7SrVQqMWHCBABAfHw8iouLsXPnTrzxxhseaVBQUBBuv/12XLx4EQCg1WrR2tqK+vp6m9Gampoah+twAEClUkGlUnmkTURENLhYDqVcv+esy3W8cuRvAODWWhoFAC13OznN7bEss9kMk8nkibYAAG7cuIGqqiqEhXXMG8bHx8PPzw9Hjx61lqmsrMSVK1fcXstDRETkSHJsGJ68Z6xbdez+zPm8NLeOClnec7eT8yQFNRkZGfjkk09w6dIllJWVISMjAwUFBVi+fDmAjvUvpaWl1lGWsrIylJaWoq7uX4ul5s6di1dffdX6/plnnkFhYSEuXbqEzz77DIsXL4avry+WLVsGANBoNFizZg02b96M48ePo6SkBE888QR0Oh13PhERUa9KcmPaRwCob27rsZyFZpifzXutxh/ZK6ZzcbAEkqafamtrsXLlSuj1emg0GsTFxeHQoUNISkoCALz++uvYtm2btfysWbMAdGzBXr16NQCgqqoK165ds5b5+uuvsWzZMly/fh0jR47Evffei5MnT2LkyH8dALZjxw74+PggJSUFJpMJ8+bNw2uvveZyp4mIiAD0mBHYsrbGYGyxm7PGGcOUvmhube+2jAKA/xAfvPfDBFy7Yeo2OzE5phBCDIqjPhsaGqDRaGA0GqFWq/u7OURE1IfsBS+HKwzdHoFgkV+ud2ttzXCVL5pM3Qc1Fu+vnQnd+BCXP0uOpPz+5tlPREQka/bObwoa5md3asiSEdiT0z5NpnbMj9XiYHnPuWp4UKV7uOmdiIhky9H5TY7WulimLrYdqEC7WXhsS/b4kcOdKseDKt3DoIaIiGSn3Sxw4sI1PPvfZZLXwnTOCOyJ4w0AQDduRLd5b3hQpWdw+omIiGTF3nSTKzw1FRSm8cfM8SHYujAGG/achQKwCbS4ddtzOFJDRESy4Wi6yRWhgf4emQ6yBCvJsWHIXjEdWo1tndy67TkcqSEiIlmwrH/x1Jbeb5taMS9W69aW7oVxWptgJTk2DEkx2m63kZPrOFJDRESy4Kn1Lxa//LhjgfDWhR2HNrsSdhw4Z0B+ud7mmq+PArrxIVh0x/egGx/CgMaDGNQQEZEseHo7tN7Ygh2H/wbNUCV+//idXaaNgm7JAGyPAv/aSUW9j9NPREQ0YLWbBU7+/TqKqq7j62+bPV7/q8cv4tXjFxGm8cfzD8XgtuFKm2mjV49dxI5/HlxpT+edVEyq1/sY1BAR0YCUX67Hs/vKJJ2v5CqDsQWpOR1J+Rbd8T3rdeN3rU49z6R6fYPTT0RENOBYji7oi4AG6JqUD+gYJfqw9KpTzzOpXt9gUENERF6h3SxQVHUd+0u/QVHVdYfrUNrNAi/82f0sv1J1nkoCOhYm1zX1PFITPNyPSfX6CKefiIio39lLmGfvcEmgI5gwNPTfdI5lKsnZKaXFd3yPO5z6CEdqiIioXzlKmGc5XPLWLdHurk+5bdgQl7ZnW1imkpydUkqM0brxaSQFgxoiIuo33SXMs7eOBXB/fcq2f48FID3vzK3nM82ICuZ5Tl6GQQ0REfWbnhLm3bqOBegIJrRq1wObfzSZ7B5X0B175zP5+igcJubjeU79g0ENERH1G2enkjqX8/VR4IWHY1z+zMt1zUiODcOnW+YgPfF2p55xdD4Tz3PyLlwoTERE/cbZqaRbyyXHhuH1FdOx+U9/RXNru6TPHBM8zPrn3OIr3ZYNGuqH3y+fjpnjHB9nwPOcvAeDGiIi6jeWdSmODoxUoGPUo/O6lHazwOnqOnzX2g5/P19JQY2PAvg/urEAnDsrqv67NvgoFD0GKJbznKh/MaghIqJ+4+ujwMPTwvDGJ9UOy3Rel2Jv67cUa++LgnJIx8oLV6a+yLsxqCEion6TX67Hm90ENOtmRVnXpVi2frtyNKSPoiOgyVjwr7U4rk59kffiQmEiIuoX3W3ntvjjma/RbhZOlXVkmJ8Pzm9LtgloAG7JliMGNURE1C+cWtPS3IYtH/wVJ/9+3eUpp+Y2M85e/rbLdW7Jlh8GNURE1C+cXavywdlv8KM9JW59VmpO18zEALdkyw3X1BARUb+QslbF2HLTrc+q/64NG/acdZhrhluy5YFBDRER9YsZUcEIGuqH+u/a+uwztx2oQFKMtkvAwi3Z8sDpJyIi6he+Pgqs/rcxffZ59o5cIHlhUENERP0iv1yP3OKvJD8XNNTP5n2Yxh9PzYrqct0R5p2RL04/ERERgH9l6nVmXYmUsvbKf9vUitQc13LO/P7x6fDxUXT57Fm3h2L5f53q8XnmnZEvBjVERGQ3U2+Yxh9bF8Z0WVgrpayj8j4KSA5oLEcmzBxv/xymmeNCJB+5QPLC6SciokHOkqn31jwwBmMLNuyx3QotpWx35c0SIxpn8sYw7wxJCmqys7MRFxcHtVoNtVoNnU6HgwcPWu+/+eabmD17NtRqNRQKBerr63usMzMzE3fffTcCAwMRGhqKRx55BJWVlTZlZs+eDYVCYfNav369lKYTEZEd3WXqtVzbdqCix6y+t5btqW6pnM0bw7wzg5uk6afRo0cjKysLEydOhBAC7777LhYtWoTPP/8cU6ZMQXNzM5KTk5GcnIyMjAyn6iwsLERqairuvvtu3Lx5E8899xwefPBBVFRUYPjw4dZya9euxYsvvmh9P2zYMHvVERGRBD1l9b11x5CzZXXjQ5zKGOyM5x+ajNX3RDk9wsK8M4OXpKBm4cKFNu9feuklZGdn4+TJk5gyZQo2bdoEACgoKHC6zvz8fJv377zzDkJDQ1FSUoJZs2ZZrw8bNgxarVZKc4mIqAe9cVK1payndhmNCFRJDkiYd2ZwcnlNTXt7O3Jzc9HU1ASdTuexBhmNRgBAcLDtQq733nsPI0aMQGxsLDIyMtDc3NxtPSaTCQ0NDTYvIiKyJeWkaqmnWntqlxF3K5GzJO9+Kisrg06nQ0tLCwICApCXl4eYmJieH3SC2WzGpk2bcM899yA2NtZ6/fHHH8eYMWMQHh6Oc+fOYcuWLaisrMS+ffsc1pWZmYlt27Z5pF1ERHJlOana2R1DUsr2VHdPuFuJpJI8UhMdHY3S0lKcOnUKGzZswKpVq1BRUeGRxqSmpqK8vBy5ubk219etW4d58+Zh6tSpWL58Of7whz8gLy8PVVVVDuvKyMiA0Wi0vr76SnqCJyIiuZOyY8iV3UWP3R3hMABSAFh731i77eJuJXKF5KBGqVRiwoQJiI+PR2ZmJqZNm4adO3e63ZC0tDR89NFHOH78OEaPHt1t2YSEBADAxYsXHZZRqVTWXVqWFxERdSVlx5CzZfPL9bj35WPYceSC3c+0lP+Ph6bg9RXTEcbdSuQBbiffM5vNMJlMLj8vhMDGjRuRl5eHgoICREVF9fhMaWkpACAsjH/ZiYg8QcqOoZ7KWnLTOJpySk+8HWlzJljLc7cSeYqkoCYjIwPz589HZGQkGhsbkZOTg4KCAhw6dAgAYDAYYDAYrCMoZWVlCAwMRGRkpHXh79y5c7F48WKkpaUB6JhyysnJwf79+xEYGAiDwQAA0Gg0GDp0KKqqqpCTk4MFCxYgJCQE586dQ3p6OmbNmoW4uDiP/SCIiAY7KTuGbi3bbhYoqroOg/E7/PLjLxwGNAoAucVXkDZngsufTeSIpKCmtrYWK1euhF6vh0ajQVxcHA4dOoSkpCQAwOuvv26zONeyJXv37t1YvXo1AKCqqgrXrl2zlsnOzgbQkWCvM8szSqUSR44cwSuvvIKmpiZEREQgJSUFP//5zyV3loiIPM/eMQiO3JrLhsiTFEIITyR79HoNDQ3QaDQwGo1cX0NE5CE9TTU5svOxO7Doju/1SptIXqT8/ubZT0RE5BJ3jkFg7hnqDTylm4iIXOLKMQjMPUO9iUENERH1qN0suuxOknoMAnPPUG9jUENERN2ytxA4TOOPx+6OkFSPVuOPrQtjmHuGeg2DGiIicsjRQmCDsQU7jlxA0DA/GJvbHK6rCR7uh+f/fQq0auaeod7HoIaIaJCwN4XUXZDR3UJgAdujEhT/vNb5PQD8avFUjsxQn2FQQ0Q0CDiaQupuOqinhcACQH1zG9ITb0du8RWbspxqov7AoIaISOa6m0LasOeswzOWnF0IPHbEMHy6ZQ6POaB+x6CGiEjGeppCAoBn95UhUOWHmeNDbAKREcNVTn1GaKA/jzkgr8Dke0REMuZMLpn65jYs33UK9758DPnlegAdozs/2fvXbp9ToGMKizlnyFtwpIaISMak5JKxTEetmxWFNz+p7jZTMHPOkDfiSA0RkYxJOY5A/PP11l+6D2gAYJRa5XAtDlF/YVBDRCRT7WYBs1kgaKifpOfMThzm9H9/cAcDGvI6nH4iIpIhe1u4PenaDVOv1EvkDgY1REQy42gLtyfxlG3yRpx+IiKSke62cFv0tKzXR+G4DHc8kTdjUENEJCPObOG2BDy3Bi6Kf77W3hfl8D7AHU/kvRjUEBHJiLNbuNfcMxZaje0Uklbjj+wV05GxIAbZK6Y7vM8FwuStuKaGiEhGnF3rkhijxXMPxTg82iA5NgxJMVoefUADCoMaIiIZmREVjDCNPwzGFrvrahToGHFxZk0Mjz6ggYZBDRHRANJuFj2Onjx2dyR2HPlbl2c7r4k5XGGQfGo3kbdjUENENEDYyz1jCUSSYrR49dgF7D5xCfXftdl9XvvPsgBcOrWbyNsxqCEiGgAc5Z4xGFuwfs9ZDFP6orm13eHz6YkTkTZnIgDg3pePOTy1WwFg24EKJMVouX6GBhzufiIi8nLd5Z6xXOsuoFEAyC3+CkDPW74FAL2xBaer61xuL1F/YVBDROTlnMk9053OgYqzW76lnO5N5C0Y1BAReTlPBRi1jS24dK3ZqbI8BoEGIgY1RERezlMBxojhKrx/+kqP5bRqFY9BoAGJQQ0RkZez5J5xZ9lumMYfUACGhp5HfZbNiOQiYRqQGNQQEXk5Xx+FdSu2q6HG1oUxuHbD5FTZsSOGu/gpRP2LQQ0RUR9oNwsUVV3H/tJvUFR1He3m7s7R7io5NszueUw9CRrmh9f/mXfG2WksrqehgYp5aoiIell3SfOkJLnrfB6ToaEFv/zoPOqa7CfaA4AAlS9OP5cI5ZCO/3/15BEKRN5I0khNdnY24uLioFaroVarodPpcPDgQev9N998E7Nnz4ZarYZCoUB9fb1T9f7+97/H2LFj4e/vj4SEBJw+fdrmfktLC1JTUxESEoKAgACkpKSgpqZGStOJiPqFJWnerVuyLdl788v1kuqznMekVft3G9AAwA1TO0ouf2vzrKNprM5HKHA9DQ1UkoKa0aNHIysrCyUlJThz5gzmzJmDRYsW4fz58wCA5uZmJCcn47nnnnO6zj/+8Y/YvHkztm7dirNnz2LatGmYN28eamtrrWXS09Nx4MAB7N27F4WFhbh69SqWLFkipelERH3OmaR52w5USJ6KApzf5n1rOUfTWFqNP49HoAFPIYSQ/q+pk+DgYGzfvh1r1qyxXisoKMADDzyAb7/9FkFBQd0+n5CQgLvvvhuvvvoqAMBsNiMiIgIbN27Es88+C6PRiJEjRyInJwdLly4FAHz55ZeYPHkyioqKMHPmTKfa2dDQAI1GA6PRCLVa7VpniYgkKKq6jmVvneyx3PtrZ0o+DdvZutMemIB7JozocvClMwdjEnkDKb+/XV4o3N7ejtzcXDQ1NUGn07lUR2trK0pKSpCYmPivBvn4IDExEUVFRQCAkpIStLW12ZSZNGkSIiMjrWXsMZlMaGhosHkREfUlZ0dTjlQYJNft7DbvV49fxLK3TuLel4/ZTHVZprEW3fE96MaHMKAhWZAc1JSVlSEgIAAqlQrr169HXl4eYmJiXPrwa9euob29HaNGjbK5PmrUKBgMHf/IDQYDlEpllxGfzmXsyczMhEajsb4iIiJcaiMRkas7l5zN3vvHM19JnoKSus3b1TU8RAOJ5KAmOjoapaWlOHXqFDZs2IBVq1ahoqKiN9rmloyMDBiNRuvrq6++6u8mEdEAlF+ux70vH8Oyt07i6dxSu6Me9oKe/HI9XjnyN6c+44apHa8euyC5bVK2ebu7hodoIJC8pVupVGLChAkAgPj4eBQXF2Pnzp144403JH/4iBEj4Ovr22UnU01NDbRaLQBAq9WitbUV9fX1NqM1ncvYo1KpoFKpJLeJiMjCsnPp1hDAMuqRvWI6AHTZrq1Vq9By02x3gbAju09cQtqciZKngTpv8z5x8R949XiVw7KdD7aUuoaHaCBwO/me2WyGyeRclspbKZVKxMfH4+jRozb1HT161LpOJz4+Hn5+fjZlKisrceXKFZfX8hAR9cSZnUsZ+8qw3t527QYT6pu73259q/rv2vDOiWqXkvNZ1sdMHBXoVHmewE1yJWmkJiMjA/Pnz0dkZCQaGxuRk5ODgoICHDp0CEDH+heDwYCLFy8C6Fh/ExgYiMjISAQHdyRzmjt3LhYvXoy0tDQAwObNm7Fq1SrcddddmDFjBl555RU0NTXhiSeeAABoNBqsWbMGmzdvRnBwMNRqNTZu3AidTuf0ziciIqlOV9d1CVY6EwC+lRi49OSXH39h/bMryfmYMZgGO0lBTW1tLVauXAm9Xg+NRoO4uDgcOnQISUlJAIDXX38d27Zts5afNWsWAGD37t1YvXo1AKCqqgrXrl2zlnn00Ufxj3/8A7/4xS9gMBhwxx13ID8/32bx8I4dO+Dj44OUlBSYTCbMmzcPr732msudJiLqSX+PZnSe4nI2sGHGYBrs3M5TM1AwTw0RSeFsHpjeZAlCPt0yx+m1NpZ1QABsAhvL00ywRwNNn+SpISKSs/gxt6G/U7d0XtjrLGYMpsGMB1oSEdlRcvlbeMvOZ6lTYZ13RDFjMA0mDGqIiOyQEkgoAEnbty3U/kPQ0HKzx3KuLOy17IgiGkw4/UREZIezgUR64kSnkt91lvrAeLy/dibO/Dyp26MOFOjYBcWFvUTOYVBDRGRHT2crWQKOtDkT8emWOXh/7UysmBnpVN0TRgZANz4EyiE+Do86sLzfujCG00ZETmJQQ0Sy4eoZTfZ0d7bSrQGHZaonKmS4U3XXNbVa/8yFvUSewzU1RCQL+eX6LscVuJLArjNLwNHlGAQH9QYHOHc0y63luLCXyDMY1BDRgOfMGU3uBDbOBhxatXNra+yV48JeIvcxqCGiAa2nM5oU6DhwMilG6/LIh7MBh2UdTnfHK3DhL1Hv4ZoaIhrQnDmjSWoCO1dZ1uEoYH8djgJc+EvUmxjUENGA5mw+mb46y4kLf4n6D6efiGhA88aTqbnwl6h/MKghogHNW0+m5sJfor7H6SciGtCk5JMhInljUENEXkdqEj131rF4MmEfEfUvTj8RkVdxNYmelHUs7WaB09V1OFJhQF7pN6hrapP0WUTknRRCiEHxvyUNDQ3QaDQwGo1Qq9X93RwissNREj2L1x6/Ewviwt3+jFuDps4sIRB3KhF5Bym/vzn9REReobskehZp73+O/zmnd/kzLEFTT3ltgI6EfZyKIhpYGNQQkVfoKYkeAJgF8KOcs8gvlx7YOBM0WfRlwj4i8hwGNUTkFaQkx3NlFMWZoMmdNhFR/2NQQ0ReQUpyPFdGUVwJUPoyYR8RuY9BDRF5BUsSPWdJDVKkBCgK8OBJooGIQQ0ReYXOSfScIXUUxRI0OZuCjwn7iAYeBjVE5DWSY8Pw2uN3ortYwtVRlO4yD3cWxoMniQYsJt8jIq+yIC4cr0KBH+Wc7XLP3WMPLJmHb81TEzJciUV3hCMpRsuDJ4kGMCbfIyKv5GpmYWdYMgrzBG0i7yfl9zeDGiLyWgw+iIgZhYlowGNAQ0RScU0NEXmd3px6IiL54kgNEXkVR+czGYwt2LDHtSMSiGhwYFBDRF6ju/OZeNAkEfVEUlCTnZ2NuLg4qNVqqNVq6HQ6HDx40Hq/paUFqampCAkJQUBAAFJSUlBTU9NtnQqFwu5r+/bt1jJjx47tcj8rK0tiV4nI2/V0PhMPmiSi7kgKakaPHo2srCyUlJTgzJkzmDNnDhYtWoTz588DANLT03HgwAHs3bsXhYWFuHr1KpYsWdJtnXq93ub19ttvQ6FQICUlxabciy++aFNu48aNErtKRP2l3SxQVHUd+0u/QVHVdYcjLc4efcCDJonIHkkLhRcuXGjz/qWXXkJ2djZOnjyJ0aNHY9euXcjJycGcOXMAALt378bkyZNx8uRJzJw5026dWq3W5v3+/fvxwAMPYNy4cTbXAwMDu5QlIu8nZdGvs0cfXGs0YX/pN9wVRUQ2XF5T097ejtzcXDQ1NUGn06GkpARtbW1ITEy0lpk0aRIiIyNRVFTkVJ01NTX4+OOPsWbNmi73srKyEBISgjvvvBPbt2/HzZs3u63LZDKhoaHB5kVEfUvqot/4Mbd1e0SCxS8//gJP55Zi2Vsnce/Lx7h4mIgAuBDUlJWVISAgACqVCuvXr0deXh5iYmJgMBigVCoRFBRkU37UqFEwGAxO1f3uu+8iMDCwy5TVj3/8Y+Tm5uL48eN46qmn8Ktf/Qo/+9nPuq0rMzMTGo3G+oqIiJDUTyJyjyuLfksufwupa4C5K4qILCTnqYmOjkZpaSmMRiM++OADrFq1CoWFhR5pzNtvv43ly5fD3992CHrz5s3WP8fFxUGpVOKpp55CZmYmVCqV3boyMjJsnmtoaGBgQ9SHpCz61Y0PAeDaWhmBjjOhth2oQFKM1vrZTNpHNPhIDmqUSiUmTJgAAIiPj0dxcTF27tyJRx99FK2traivr7cZrampqXFqLcxf/vIXVFZW4o9//GOPZRMSEnDz5k1cunQJ0dHRdsuoVCqHAQ8R9T5XFv06u6bmVpYA6dVjF5FbfIVJ+4gGKbfz1JjNZphMJsTHx8PPzw9Hjx613qusrMSVK1eg0+l6rGfXrl2Ij4/HtGnTeixbWloKHx8fhIaGutV2Iuo9zgYoncvNiApGmMYfro6r7DjyNybtIxrEJAU1GRkZ+OSTT3Dp0iWUlZUhIyMDBQUFWL58OTQaDdasWYPNmzfj+PHjKCkpwRNPPAGdTmez82nSpEnIy8uzqbehoQF79+7FD3/4wy6fWVRUhFdeeQV//etf8fe//x3vvfce0tPTsWLFCtx2220udpuIeltPAYoCHaMoM6KCrdd8fRTYujDGet8TmLSPaPCQFNTU1tZi5cqViI6Oxty5c1FcXIxDhw4hKSkJALBjxw78+7//O1JSUjBr1ixotVrs27fPpo7KykoYjUaba7m5uRBCYNmyZV0+U6VSITc3F/fffz+mTJmCl156Cenp6XjzzTel9pWI+lB3AYrl/daFMV3WuyTHhiF7xXRoNbYjPe4si2HSPqLBQSGEGBT/6yLl6HIi8hxXD6e89ZTub5tMSM35HABsdlQpbnnfnZ2P3YFFd3xPeieIqN9I+f3NU7qJqFclx4YhKUYreUeSr4/CuivKIttH0SVA0mr88djdEdhx5EKPbXF1ITIRDQwMaoio19kLUFzhKEACgNzir2AwttgdtVGgI/jpvH6HiOSHQQ0RDSiOAqStC2OwYc/ZLtNR3a3fISJ5cXtLNxGRN3C0wFir8Uf2iunMU0M0CHCkhohkw9X1O0QkDwxqiEhWPLV+h4gGHk4/ERERkSwwqCEiIiJZYFBDREREssCghoiIiGSBQQ0RERHJAoMaIiIikgUGNURERCQLDGqIiIhIFhjUEBERkSwwqCEiIiJZYFBDREREssCghoiIiGSBQQ0RERHJAoMaIiIikgUGNURERCQLQ/q7AUTkvHazwOnqOtQ2tiA00B8zooLh66Po72YREXkFBjVEA0R+uR7bDlRAb2yxXgvT+GPrwhgkx4b1Y8uIiLwDp5+IBoD8cj027DlrE9AAgMHYgg17ziK/XN9PLSMi8h4Maoi8XLtZYNuBCgg79yzXth2oQLvZXgkiosGDQQ2RlztdXddlhKYzAUBvbMHp6rq+axQRkRdiUEPk5WobHQc0rpQjIpIrBjVEXi400N+j5YiI5IpBDZGXmxEVjDCNPxxt3FagYxfUjKjgvmwWEZHXYVBD5OV8fRTYujAGALoENpb3WxfGMF8NEQ16DGqIBoDk2DBkr5gOrcZ2ikmr8Uf2iunMU0NEBIlBTXZ2NuLi4qBWq6FWq6HT6XDw4EHr/ZaWFqSmpiIkJAQBAQFISUlBTU1Nt3WuXr0aCoXC5pWcnGxTpq6uDsuXL4darUZQUBDWrFmDGzduSGk60YCXHBuGT7fMwftrZ2LnY3fg/bUz8emWOQxoiIj+SVJG4dGjRyMrKwsTJ06EEALvvvsuFi1ahM8//xxTpkxBeno6Pv74Y+zduxcajQZpaWlYsmQJTpw40W29ycnJ2L17t/W9SqWyub98+XLo9XocPnwYbW1teOKJJ7Bu3Trk5ORIaT7RgOfro4BufEh/N4OIyCsphBBuZewKDg7G9u3bsXTpUowcORI5OTlYunQpAODLL7/E5MmTUVRUhJkzZ9p9fvXq1aivr8eHH35o9/4XX3yBmJgYFBcX46677gIA5OfnY8GCBfj6668RHh7uVDsbGhqg0WhgNBqhVquld5SIiIj6nJTf3y6vqWlvb0dubi6ampqg0+lQUlKCtrY2JCYmWstMmjQJkZGRKCoq6raugoIChIaGIjo6Ghs2bMD169et94qKihAUFGQNaAAgMTERPj4+OHXqlMM6TSYTGhoabF5EREQkX5KDmrKyMgQEBEClUmH9+vXIy8tDTEwMDAYDlEolgoKCbMqPGjUKBoPBYX3Jycn4wx/+gKNHj+Lll19GYWEh5s+fj/b2dgCAwWBAaGiozTNDhgxBcHBwt/VmZmZCo9FYXxEREVK7SkRERAOI5FO6o6OjUVpaCqPRiA8++ACrVq1CYWGhyw147LHHrH+eOnUq4uLiMH78eBQUFGDu3Lku15uRkYHNmzdb3zc0NDCwISIikjHJQY1SqcSECRMAAPHx8SguLsbOnTvx6KOPorW1FfX19TajNTU1NdBqtU7XP27cOIwYMQIXL17E3LlzodVqUVtba1Pm5s2bqKur67ZelUrVZcExERERyZfbeWrMZjNMJhPi4+Ph5+eHo0ePWu9VVlbiypUr0Ol0Ttf39ddf4/r16wgL69imqtPpUF9fj5KSEmuZY8eOwWw2IyEhwd3mExERkUxICmoyMjLwySef4NKlSygrK0NGRgYKCgqwfPlyaDQarFmzBps3b8bx48dRUlKCJ554Ajqdzmbn06RJk5CXlwcAuHHjBn7605/i5MmTuHTpEo4ePYpFixZhwoQJmDdvHgBg8uTJSE5Oxtq1a3H69GmcOHECaWlpeOyxx5ze+URERETyJ2n6qba2FitXroRer4dGo0FcXBwOHTqEpKQkAMCOHTvg4+ODlJQUmEwmzJs3D6+99ppNHZWVlTAajQAAX19fnDt3Du+++y7q6+sRHh6OBx98EL/85S9tpo7ee+89pKWlYe7cudb6f/vb37rbdyIiIpIRt/PUDBTMU0NERDTw9EmeGiIiIiJvwqCGiIiIZIFBDREREckCgxoiIiKSBQY1REREJAsMaoiIiEgWGNQQERGRLDCoISIiIllgUENERESywKCGiIiIZIFBDREREckCgxoiIiKSBQY1REREJAsMaoiIiEgWGNQQERGRLDCoISIiIllgUENERESywKCGiIiIZIFBDREREcnCkP5uANFA0W4WOF1dh9rGFoQG+mNGVDB8fRT93SwiIvonBjVETsgv12PbgQrojS3Wa2Eaf2xdGIPk2LB+bBkREVlw+omoB/nlemzYc9YmoAEAg7EFG/acRX65vp9aRkREnTGoIepGu1lg24EKCDv3LNe2HahAu9leCSIi6ksMaoi6cbq6rssITWcCgN7YgtPVdX3XKCIisotBDVE3ahsdBzSulCMiot7DoIaoG6GB/h4tR0REvYdBDVE3ZkQFI0zjD0cbtxXo2AU1Iyq4L5tFRER2MKgh6oavjwJbF8YAQJfAxvJ+68IY5qshIvICDGqIepAcG4bsFdOh1dhOMWk1/sheMZ15aoiIvAST7xE5ITk2DEkxWmYUJiLyYgxqiJzk66OAbnxIfzeDiIgckDT9lJ2djbi4OKjVaqjVauh0Ohw8eNB6v6WlBampqQgJCUFAQABSUlJQU1PjsL62tjZs2bIFU6dOxfDhwxEeHo6VK1fi6tWrNuXGjh0LhUJh88rKypLYVSIiIpIzSUHN6NGjkZWVhZKSEpw5cwZz5szBokWLcP78eQBAeno6Dhw4gL1796KwsBBXr17FkiVLHNbX3NyMs2fP4vnnn8fZs2exb98+VFZW4uGHH+5S9sUXX4Rer7e+Nm7cKLGrREREJGcKIYRb+d2Dg4Oxfft2LF26FCNHjkROTg6WLl0KAPjyyy8xefJkFBUVYebMmU7VV1xcjBkzZuDy5cuIjIwE0DFSs2nTJmzatMnldjY0NECj0cBoNEKtVrtcDxEREfUdKb+/Xd791N7ejtzcXDQ1NUGn06GkpARtbW1ITEy0lpk0aRIiIyNRVFTkdL1GoxEKhQJBQUE217OyshASEoI777wT27dvx82bN7utx2QyoaGhweZFRERE8iV5oXBZWRl0Oh1aWloQEBCAvLw8xMTEoLS0FEqlskswMmrUKBgMBqfqbmlpwZYtW7Bs2TKbaOzHP/4xpk+fjuDgYHz22WfIyMiAXq/Hb37zG4d1ZWZmYtu2bVK7R0RERAOU5KAmOjoapaWlMBqN+OCDD7Bq1SoUFha63ZC2tjb84Ac/gBAC2dnZNvc2b95s/XNcXByUSiWeeuopZGZmQqVS2a0vIyPD5rmGhgZERES43U4iIiLyTpKDGqVSiQkTJgAA4uPjUVxcjJ07d+LRRx9Fa2sr6uvrbUZrampqoNVqu63TEtBcvnwZx44d63HOLCEhATdv3sSlS5cQHR1tt4xKpXIY8BAREZH8uJ1R2Gw2w2QyIT4+Hn5+fjh69Kj1XmVlJa5cuQKdTufweUtAc+HCBRw5cgQhIT3nASktLYWPjw9CQ0PdbT4RERHJhKSRmoyMDMyfPx+RkZFobGxETk4OCgoKcOjQIWg0GqxZswabN29GcHAw1Go1Nm7cCJ1OZ7PzadKkScjMzMTixYvR1taGpUuX4uzZs/joo4/Q3t5uXX8THBwMpVKJoqIinDp1Cg888AACAwNRVFSE9PR0rFixArfddptnfxpEREQ0YEkKampra7Fy5Uro9XpoNBrExcXh0KFDSEpKAgDs2LEDPj4+SElJgclkwrx58/Daa6/Z1FFZWQmj0QgA+Oabb/DnP/8ZAHDHHXfYlDt+/Dhmz54NlUqF3NxcvPDCCzCZTIiKikJ6errNehkiIiIit/PUDBTMU0NERDTw9EmeGiIiIiJvwqCGiIiIZIFBDREREckCgxoiIiKSBQY1REREJAsMaoiIiEgWGNQQERGRLDCoISIiIllgUENERESywKCGiIiIZIFBDREREckCgxoiIiKSBQY1REREJAsMaoiIiEgWGNQQERGRLDCoISIiIllgUENERESywKCGiIiIZIFBDREREckCgxoiIiKSBQY1REREJAsMaoiIiEgWGNQQERGRLDCoISIiIllgUENERESywKCGiIiIZIFBDREREckCgxoiIiKSBQY1REREJAsMaoiIiEgWGNQQERGRLEgKarKzsxEXFwe1Wg21Wg2dToeDBw9a77e0tCA1NRUhISEICAhASkoKampquq1TCIFf/OIXCAsLw9ChQ5GYmIgLFy7YlKmrq8Py5cuhVqsRFBSENWvW4MaNG1KaTkRERDInKagZPXo0srKyUFJSgjNnzmDOnDlYtGgRzp8/DwBIT0/HgQMHsHfvXhQWFuLq1atYsmRJt3X++te/xm9/+1u8/vrrOHXqFIYPH4558+ahpaXFWmb58uU4f/48Dh8+jI8++giffPIJ1q1b50J3B692s0BR1XXsL/0GRVXX0W4W/d0kIiIij1IIIdz67RYcHIzt27dj6dKlGDlyJHJycrB06VIAwJdffonJkyejqKgIM2fO7PKsEALh4eH4yU9+gmeeeQYAYDQaMWrUKLzzzjt47LHH8MUXXyAmJgbFxcW46667AAD5+flYsGABvv76a4SHhzvVzoaGBmg0GhiNRqjVane6PODkl+ux7UAF9MZ/BYphGn9sXRiD5NiwfmwZERFR96T8/nZ5TU17eztyc3PR1NQEnU6HkpIStLW1ITEx0Vpm0qRJiIyMRFFRkd06qqurYTAYbJ7RaDRISEiwPlNUVISgoCBrQAMAiYmJ8PHxwalTpxy2z2QyoaGhweY1GOWX67Fhz1mbgAYADMYWbNhzFvnl+n5qGRERkWdJDmrKysoQEBAAlUqF9evXIy8vDzExMTAYDFAqlQgKCrIpP2rUKBgMBrt1Wa6PGjXK4TMGgwGhoaE294cMGYLg4GCH9QJAZmYmNBqN9RURESG1qwNeu1lg24EK2BuKs1zbdqCCU1FERCQLkoOa6OholJaW4tSpU9iwYQNWrVqFioqK3mibWzIyMmA0Gq2vr776qr+b1OdOV9d1GaHpTADQG1twurqu7xpFRETUS4ZIfUCpVGLChAkAgPj4eBQXF2Pnzp149NFH0draivr6epvRmpqaGmi1Wrt1Wa7X1NQgLCzM5pk77rjDWqa2ttbmuZs3b6Kurs5hvQCgUqmgUqmkdk9WahsdBzSulCMiIvJmbuepMZvNMJlMiI+Ph5+fH44ePWq9V1lZiStXrkCn09l9NioqClqt1uaZhoYGnDp1yvqMTqdDfX09SkpKrGWOHTsGs9mMhIQEd5sva6GB/h4tR0RE5M0kjdRkZGRg/vz5iIyMRGNjI3JyclBQUIBDhw5Bo9FgzZo12Lx5M4KDg6FWq7Fx40bodDqbnU+TJk1CZmYmFi9eDIVCgU2bNuE///M/MXHiRERFReH5559HeHg4HnnkEQDA5MmTkZycjLVr1+L1119HW1sb0tLS8Nhjjzm982mwmhEVjDCNPwzGFrvrahQAtBp/zIgK7uumEREReZykoKa2thYrV66EXq+HRqNBXFwcDh06hKSkJADAjh074OPjg5SUFJhMJsybNw+vvfaaTR2VlZUwGo3W9z/72c/Q1NSEdevWob6+Hvfeey/y8/Ph7/+v0YP33nsPaWlpmDt3rrX+3/72t+70e1Dw9VFg68IYbNhzFgrAJrBR/PO/WxfGwNdHYedpIiKigcXtPDUDBfPUME8NERENPFJ+f0teKEy22s0Cp6vrUNvYgtDAjqkcbxv5SI4NQ1KM1uvbSURE5A4GNW4YSCMgvj4K6MaH9HcziIiIeg1P6XYRM/USERF5FwY1LmCmXiIiIu/DoMYFzNRLRETkfRjUuICZeomIiLwPgxoXMFMvERGR92FQ4wJLpl5HG6IV6NgFxUy9REREfYdBjQssmXoBdAlsmKmXiIiofzCocVFybBiyV0yHVmM7xaTV+CN7xXSvy1NDREQkd0y+5wZm6iUiIvIeDGrcxEy9RERE3oHTT0RERCQLDGqIiIhIFhjUEBERkSwwqCEiIiJZYFBDREREssCghoiIiGSBQQ0RERHJAoMaIiIikgUGNURERCQLDGqIiIhIFnhMgoe0mwXPgCIiIupHDGo8IL9cj20HKqA3tlivhWn8sXVhDE/rJiIi6iOcfnJTfrkeG/actQloAMBgbMGGPWeRX67vp5YRERENLgxq3NBuFth2oALCzj3LtW0HKtButleCiIiIPIlBjRtOV9d1GaHpTADQG1twurqu7xpFREQ0SDGocUNto+OAxpVyRERE5DoGNW4IDfT3aDkiIiJyHYMaN8yICkaYxh+ONm4r0LELakZUcF82i4iIaFCSFNRkZmbi7rvvRmBgIEJDQ/HII4+gsrLSpkxVVRUWL16MkSNHQq1W4wc/+AFqamq6rXfs2LFQKBRdXqmpqdYys2fP7nJ//fr1Uprvcb4+CmxdGAMAXQIby/utC2OYr4aIiKgPSApqCgsLkZqaipMnT+Lw4cNoa2vDgw8+iKamJgBAU1MTHnzwQSgUChw7dgwnTpxAa2srFi5cCLPZ7LDe4uJi6PV66+vw4cMAgO9///s25dauXWtT7te//rXU/npccmwYsldMh1ZjO8Wk1fgje8V05qkhIiLqI5KS7+Xn59u8f+eddxAaGoqSkhLMmjULJ06cwKVLl/D5559DrVYDAN59913cdtttOHbsGBITE+3WO3LkSJv3WVlZGD9+PO6//36b68OGDYNWq5XS5D6RHBuGpBgtMwoTERH1I7fW1BiNRgBAcHDHmhGTyQSFQgGVSmUt4+/vDx8fH3z66adO1dna2oo9e/bgySefhEJhGxS89957GDFiBGJjY5GRkYHm5mZ3mu9Rvj4K6MaHYNEd34NufAgDGiIioj7m8jEJZrMZmzZtwj333IPY2FgAwMyZMzF8+HBs2bIFv/rVryCEwLPPPov29nbo9c5l1v3www9RX1+P1atX21x//PHHMWbMGISHh+PcuXPYsmULKisrsW/fPrv1mEwmmEwm6/uGhgbXOkpEREQDgssjNampqSgvL0dubq712siRI7F3714cOHAAAQEB0Gg0qK+vx/Tp0+Hj49xH7dq1C/Pnz0d4eLjN9XXr1mHevHmYOnUqli9fjj/84Q/Iy8tDVVWV3XoyMzOh0Wisr4iICFe7SkRERAOASyM1aWlp+Oijj/DJJ59g9OjRNvcefPBBVFVV4dq1axgyZAiCgoKg1Woxbty4Huu9fPkyjhw54nD0pbOEhAQAwMWLFzF+/Pgu9zMyMrB582br+4aGBgY2REREMiYpqBFCYOPGjcjLy0NBQQGioqIclh0xYgQA4NixY6itrcXDDz/cY/27d+9GaGgoHnrooR7LlpaWAgDCwuzvLlKpVDZre4iIiEjeJAU1qampyMnJwf79+xEYGAiDwQAA0Gg0GDp0KICOwGTy5MkYOXIkioqK8PTTTyM9PR3R0dHWeubOnYvFixcjLS3Nes1sNmP37t1YtWoVhgyxbVZVVRVycnKwYMEChISE4Ny5c0hPT8esWbMQFxfncueJiIhIPiQFNdnZ2QA6EuF1tnv3buvC3srKSmRkZKCurg5jx47Ff/zHfyA9Pd2mvGV6qrMjR47gypUrePLJJ7t8rlKpxJEjR/DKK6+gqakJERERSElJwc9//nMpzSciIiIZUwghRH83oi80NDRAo9HAaDRac+gQERGRd5Py+5tnPxEREZEsMKghIiIiWXA5+d5AY5llYxI+IiKigcPye9uZ1TKDJqhpbGwEAOaqISIiGoAaGxuh0Wi6LTNoFgqbzWZcvXoVgYGBXc6UGugsiQW/+uqrQbEIejD1dzD1FWB/5Www9RUYXP3t7b4KIdDY2Ijw8PAeTycYNCM1Pj4+XbIfy41arZb9P57OBlN/B1NfAfZXzgZTX4HB1d/e7GtPIzQWXChMREREssCghoiIiGSBQY0MqFQqbN26ddCcdTWY+juY+gqwv3I2mPoKDK7+elNfB81CYSIiIpI3jtQQERGRLDCoISIiIllgUENERESywKCGiIiIZIFBTT/Lzs5GXFycNWmRTqfDwYMHrfdbWlqQmpqKkJAQBAQEICUlBTU1Nd3WKYTAL37xC4SFhWHo0KFITEzEhQsXbMrU1dVh+fLlUKvVCAoKwpo1a3Djxo1e6aOFp/va1taGLVu2YOrUqRg+fDjCw8OxcuVKXL161abc2LFjoVAobF5ZWVm91k+L3vhuV69e3aUvycnJNmXk8N0C6NJPy2v79u3WMt763b755puYPXs21Go1FAoF6uvrnar397//PcaOHQt/f38kJCTg9OnTNvdd+Tm6qzf6mpmZibvvvhuBgYEIDQ3FI488gsrKSpsys2fP7vLdrl+/3tPd66I3+vvCCy906cukSZNsyvTHdwv0Tn/t/btUKBRITU21lum171dQv/rzn/8sPv74Y/G3v/1NVFZWiueee074+fmJ8vJyIYQQ69evFxEREeLo0aPizJkzYubMmeLf/u3fuq0zKytLaDQa8eGHH4q//vWv4uGHHxZRUVHiu+++s5ZJTk4W06ZNEydPnhR/+ctfxIQJE8SyZcsGVF/r6+tFYmKi+OMf/yi+/PJLUVRUJGbMmCHi4+Ntyo0ZM0a8+OKLQq/XW183btzo1b4K0Tvf7apVq0RycrJNX+rq6mzKyOG7FULY9FGv14u3335bKBQKUVVVZS3jrd/tjh07RGZmpsjMzBQAxLfffttjnbm5uUKpVIq3335bnD9/Xqxdu1YEBQWJmpoaaxlXfo7u6o2+zps3T+zevVuUl5eL0tJSsWDBAhEZGWnz3d1///1i7dq1Nt+t0WjsrW5a9UZ/t27dKqZMmWLTl3/84x82ZfrjuxWid/pbW1tr09fDhw8LAOL48ePWMr31/TKo8UK33Xab+K//+i9RX18v/Pz8xN69e633vvjiCwFAFBUV2X3WbDYLrVYrtm/fbr1WX18vVCqVeP/994UQQlRUVAgAori42Frm4MGDQqFQiG+++aaXemWfO3215/Tp0wKAuHz5svXamDFjxI4dOzzZbJe5299Vq1aJRYsWObwv5+920aJFYs6cOTbXvPG77ez48eNO/yKYMWOGSE1Ntb5vb28X4eHhIjMzUwghPPZz9AR3+3qr2tpaAUAUFhZar91///3i6aefdrOlnuFuf7du3SqmTZvm8L43fbdCeP77ffrpp8X48eOF2Wy2Xuut75fTT16kvb0dubm5aGpqgk6nQ0lJCdra2pCYmGgtM2nSJERGRqKoqMhuHdXV1TAYDDbPaDQaJCQkWJ8pKipCUFAQ7rrrLmuZxMRE+Pj44NSpU73UO1ue6Ks9RqMRCoUCQUFBNtezsrIQEhKCO++8E9u3b8fNmzc91RWneLK/BQUFCA0NRXR0NDZs2IDr169b78n1u62pqcHHH3+MNWvWdLnnbd+tK1pbW1FSUmLzM/Lx8UFiYqL1Z+SpfyPu8ERf7TEajQCA4OBgm+vvvfceRowYgdjYWGRkZKC5udljn+kMT/b3woULCA8Px7hx47B8+XJcuXLFes8bvlugd77f1tZW7NmzB08++WSXw6R74/sdNAdaerOysjLodDq0tLQgICAAeXl5iImJQWlpKZRKZZdf0KNGjYLBYLBbl+X6qFGjHD5jMBgQGhpqc3/IkCEIDg52WK+neLKvt2ppacGWLVuwbNkym0PVfvzjH2P69OkIDg7GZ599hoyMDOj1evzmN7/xZNfs8nR/k5OTsWTJEkRFRaGqqgrPPfcc5s+fj6KiIvj6+sr2u3333XcRGBiIJUuW2Fz3xu/WFdeuXUN7e7vdf7dffvklgI5/t+7+HF3lyb7eymw2Y9OmTbjnnnsQGxtrvf74449jzJgxCA8Px7lz57BlyxZUVlZi3759Hvnc7ni6vwkJCXjnnXcQHR0NvV6Pbdu24b777kN5eTkCAwP79bsFevf7/fDDD1FfX4/Vq1fbXO+t75dBjReIjo5GaWkpjEYjPvjgA6xatQqFhYX93axe0Vt9bWtrww9+8AMIIZCdnW1zb/PmzdY/x8XFQalU4qmnnkJmZmavp/X2dH8fe+wx65+nTp2KuLg4jB8/HgUFBZg7d64nmuyy3vx7/Pbbb2P58uXw9/e3ue6N362nfhl4k97sa2pqKsrLy/Hpp5/aXF+3bp31z1OnTkVYWBjmzp2LqqoqjB8/3u3P7Y6n+zt//nzrn+Pi4pCQkIAxY8bgT3/6k93Rx77Wm9/vrl27MH/+fISHh9tc763vl9NPXkCpVGLChAmIj49HZmYmpk2bhp07d0Kr1aK1tbXLavOamhpotVq7dVmu37pqvvMzWq0WtbW1Nvdv3ryJuro6h/V6iif7amEJaC5fvozDhw/bjNLYk5CQgJs3b+LSpUtu9qZnvdHfzsaNG4cRI0bg4sWLAOT33QLAX/7yF1RWVuKHP/xhj2W94bt1xYgRI+Dr69vjv1tP/J1xhSf72llaWho++ugjHD9+HKNHj+62bEJCAgBY/673pt7qr0VQUBBuv/12m3+3/fXdAr3X38uXL+PIkSNO/9sF3P9+GdR4IbPZDJPJhPj4ePj5+eHo0aPWe5WVlbhy5YrD+c6oqChotVqbZxoaGnDq1CnrMzqdDvX19SgpKbGWOXbsGMxms/UvVl9xp6/AvwKaCxcu4MiRIwgJCenxM0tLS+Hj49NlmqYvuNvfW3399de4fv06wsLCAMjru7XYtWsX4uPjMW3atB7LesN36wqlUon4+Hibn5HZbMbRo0etPyNP/Z3xBHf6CnSknUhLS0NeXh6OHTuGqKioHp8pLS0FAOvf9b7kbn9vdePGDVRVVVn74k3fLeC5/u7evRuhoaF46KGHeizrse/X40uPSZJnn31WFBYWiurqanHu3Dnx7LPPCoVCIf73f/9XCNGxzS8yMlIcO3ZMnDlzRuh0OqHT6WzqiI6OFvv27bO+z8rKEkFBQWL//v3i3LlzYtGiRXa3dN95553i1KlT4tNPPxUTJ07s9W2/nu5ra2urePjhh8Xo0aNFaWmpzdZAk8kkhBDis88+Ezt27BClpaWiqqpK7NmzR4wcOVKsXLmyV/vaG/1tbGwUzzzzjCgqKhLV1dXiyJEjYvr06WLixImipaXF+owcvlsLo9Eohg0bJrKzs7t8pjd/t3q9Xnz++efirbfeEgDEJ598Ij7//HNx/fp1ax1z5swRv/vd76zvc3NzhUqlEu+8846oqKgQ69atE0FBQcJgMFjLOPNzHAh93bBhg9BoNKKgoMDm321zc7MQQoiLFy+KF198UZw5c0ZUV1eL/fv3i3HjxolZs2b1al97q78/+clPREFBgaiurhYnTpwQiYmJYsSIEaK2ttZapj++297qrxAdu/ciIyPFli1bunxmb36/DGr62ZNPPinGjBkjlEqlGDlypJg7d671L5MQQnz33XfiRz/6kbjtttvEsGHDxOLFi4Ver7epA4DYvXu39b3ZbBbPP/+8GDVqlFCpVGLu3LmisrLS5pnr16+LZcuWiYCAAKFWq8UTTzwhGhsbB1Rfq6urBQC7L0s+hJKSEpGQkCA0Go3w9/cXkydPFr/61a9sgoCB0t/m5mbx4IMPipEjRwo/Pz8xZswYsXbtWptfekLI47u1eOONN8TQoUNFfX19l8/05u9269atdv9edu7fmDFjxNatW23q/d3vficiIyOFUqkUM2bMECdPnrS578zP0dN6o6+O/t1anrly5YqYNWuWCA4OFiqVSkyYMEH89Kc/7ZM8Nb3R30cffVSEhYUJpVIpvve974lHH31UXLx40eZz++O7FaL3/i4fOnRIAOjyu0eI3v1+FUII4d5YDxEREVH/45oaIiIikgUGNURERCQLDGqIiIhIFhjUEBERkSwwqCEiIiJZYFBDREREssCghoiIiGSBQQ0RERHJAoMaIiIikgUGNURERCQLDGqIiIhIFhjUEBERkSz8f5TsShC908s9AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "# plt.hist(mock[\"mu_TFR\"] - mock[\"mu_calibration\"][1])\n",
    "plt.scatter(mock[\"mu_true\"], mock[\"mu_calibration\"][0])\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([125])"
      ]
     },
     "execution_count": 127,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum(np.isfinite(mock[\"mu_calibration\"]), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAhOUlEQVR4nO3dfXBU5f2/8feGkBADuyEREiKJZJQRrAgKilEHETI8DoKmtVTaWmWIVUAhM2Li8FD7VYMMRQaLoFYjVKhKp6BCjbVRodYYIIgoRR4qSDRssIPZhSAhkvv3h8P+XIgJkLOce5PrNXNmmnNOls/tQnL17JPHGGMEAABgkRi3BwAAADgVgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOrFuD3AuGhoaVFVVpU6dOsnj8bg9DgAAOAPGGB0+fFjp6emKiWn6GklUBkpVVZUyMjLcHgMAAJyDyspKde/evclzojJQOnXqJOn7BXq9XpenAQAAZyIYDCojIyP0e7wpURkoJx/W8Xq9BAoAAFHmTJ6ewZNkAQCAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgnVi3BwDQ+vQoWNei7983d7RDkwCIVlxBAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1uHTjAFYh09DBsAVFAAAYJ2zDpQNGzZozJgxSk9Pl8fj0Zo1a0LH6uvr9dBDD6lPnz5KTExUenq6fv3rX6uqqirsNg4dOqQJEybI6/UqKSlJEydO1JEjR1q8GAAA0DqcdaDU1taqb9++Wrx48WnHjh49qi1btmjWrFnasmWL/va3v2nnzp265ZZbws6bMGGCtm/frrfffltr167Vhg0blJeXd+6rAAAArYrHGGPO+Zs9Hq1evVrjxo370XM2bdqka6+9Vl988YUyMzO1Y8cOXX755dq0aZMGDBggSSopKdGoUaP05ZdfKj09vdk/NxgMyufzKRAIyOv1nuv4ACKkpc8haSmegwLY6Wx+f0f8OSiBQEAej0dJSUmSpLKyMiUlJYXiRJJycnIUExOj8vLyRm+jrq5OwWAwbAMAAK1XRAPl2LFjeuihh/SLX/wiVEp+v19du3YNOy82NlbJycny+/2N3k5RUZF8Pl9oy8jIiOTYAADAZRELlPr6et1+++0yxmjJkiUtuq3CwkIFAoHQVllZ6dCUAADARhF5H5STcfLFF1/onXfeCXucKS0tTQcPHgw7/7vvvtOhQ4eUlpbW6O3Fx8crPj4+EqMCAAALOX4F5WSc7N69W//85z+VkpISdjw7O1s1NTWqqKgI7XvnnXfU0NCggQMHOj0OAACIQmd9BeXIkSPas2dP6Ou9e/dq69atSk5OVrdu3fTTn/5UW7Zs0dq1a3XixInQ80qSk5MVFxen3r17a8SIEZo0aZKWLl2q+vp6TZkyRePHjz+jV/AAAIDW76wDZfPmzbr55ptDX+fn50uS7rzzTv3ud7/T66+/Lknq169f2Pe9++67Gjx4sCRpxYoVmjJlioYOHaqYmBjl5uZq0aJF57gEAADQ2px1oAwePFhNvXXKmbytSnJyslauXHm2fzQAAGgj+CweAABgHQIFAABYh0ABAADWIVAAAIB1IvJGbQCim9sf9gcAXEEBAADWIVAAAIB1CBQAAGAdAgUAAFiHJ8kCrRBPcgUQ7biCAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsE6s2wMAOF2PgnVujwAAruIKCgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALDOWQfKhg0bNGbMGKWnp8vj8WjNmjVhx40xmj17trp166aEhATl5ORo9+7dYeccOnRIEyZMkNfrVVJSkiZOnKgjR460aCEAAKD1OOtAqa2tVd++fbV48eJGj8+bN0+LFi3S0qVLVV5ersTERA0fPlzHjh0LnTNhwgRt375db7/9ttauXasNGzYoLy/v3FcBAABalbP+LJ6RI0dq5MiRjR4zxmjhwoWaOXOmxo4dK0lavny5UlNTtWbNGo0fP147duxQSUmJNm3apAEDBkiSnnrqKY0aNUrz589Xenp6C5YDAABaA0efg7J37175/X7l5OSE9vl8Pg0cOFBlZWWSpLKyMiUlJYXiRJJycnIUExOj8vLyRm+3rq5OwWAwbAMAAK2Xo4Hi9/slSampqWH7U1NTQ8f8fr+6du0adjw2NlbJycmhc05VVFQkn88X2jIyMpwcGwAAWCYqXsVTWFioQCAQ2iorK90eCQAARJCjgZKWliZJqq6uDttfXV0dOpaWlqaDBw+GHf/uu+906NCh0Dmnio+Pl9frDdsAAEDr5WigZGVlKS0tTaWlpaF9wWBQ5eXlys7OliRlZ2erpqZGFRUVoXPeeecdNTQ0aODAgU6OAwAAotRZv4rnyJEj2rNnT+jrvXv3auvWrUpOTlZmZqamTZumRx99VD179lRWVpZmzZql9PR0jRs3TpLUu3dvjRgxQpMmTdLSpUtVX1+vKVOmaPz48byCBwAASDqHQNm8ebNuvvnm0Nf5+fmSpDvvvFMvvviiZsyYodraWuXl5ammpkY33nijSkpK1KFDh9D3rFixQlOmTNHQoUMVExOj3NxcLVq0yIHlAACA1sBjjDFuD3G2gsGgfD6fAoEAz0dBq9SjYJ3bI0S1fXNHuz0CgEacze/vqHgVDwAAaFsIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWiXV7AKA16lGwzu0R0AItvf/2zR3t0CRA28UVFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADW4X1QALQ6vA8NEP24ggIAAKxDoAAAAOsQKAAAwDqOB8qJEyc0a9YsZWVlKSEhQZdccon+7//+T8aY0DnGGM2ePVvdunVTQkKCcnJytHv3bqdHAQAAUcrxQHniiSe0ZMkS/fGPf9SOHTv0xBNPaN68eXrqqadC58ybN0+LFi3S0qVLVV5ersTERA0fPlzHjh1zehwAABCFHH8VzwcffKCxY8dq9OjvP82zR48e+stf/qKNGzdK+v7qycKFCzVz5kyNHTtWkrR8+XKlpqZqzZo1Gj9+vNMjAQCAKOP4FZTrr79epaWl2rVrlyTp448/1vvvv6+RI0dKkvbu3Su/36+cnJzQ9/h8Pg0cOFBlZWVOjwMAAKKQ41dQCgoKFAwG1atXL7Vr104nTpzQY489pgkTJkiS/H6/JCk1NTXs+1JTU0PHTlVXV6e6urrQ18Fg0OmxAQCARRy/gvLqq69qxYoVWrlypbZs2aJly5Zp/vz5WrZs2TnfZlFRkXw+X2jLyMhwcGIAAGAbxwPlwQcfVEFBgcaPH68+ffroV7/6laZPn66ioiJJUlpamiSpuro67Puqq6tDx05VWFioQCAQ2iorK50eGwAAWMTxQDl69KhiYsJvtl27dmpoaJAkZWVlKS0tTaWlpaHjwWBQ5eXlys7ObvQ24+Pj5fV6wzYAANB6Of4clDFjxuixxx5TZmamfvKTn+ijjz7SggULdPfdd0uSPB6Ppk2bpkcffVQ9e/ZUVlaWZs2apfT0dI0bN87pcQAAQBRyPFCeeuopzZo1S/fdd58OHjyo9PR03XPPPZo9e3bonBkzZqi2tlZ5eXmqqanRjTfeqJKSEnXo0MHpcQAAQBTymB++xWuUCAaD8vl8CgQCPNwDK/Fpum3bvrmj3R4BsNLZ/P7ms3gAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1Yt0eALANn0QMAO7jCgoAALAOgQIAAKzDQzwA4DAnHibcN3e0A5MA0YsrKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsE6s2wMATutRsM7tEQAALcQVFAAAYJ2IBMpXX32lX/7yl0pJSVFCQoL69OmjzZs3h44bYzR79mx169ZNCQkJysnJ0e7duyMxCgAAiEKOB8o333yjG264Qe3bt9ebb76p//znP/rDH/6gzp07h86ZN2+eFi1apKVLl6q8vFyJiYkaPny4jh075vQ4AAAgCjn+HJQnnnhCGRkZKi4uDu3LysoK/W9jjBYuXKiZM2dq7NixkqTly5crNTVVa9as0fjx450eCQAARBnHr6C8/vrrGjBggH72s5+pa9euuuqqq/Tcc8+Fju/du1d+v185OTmhfT6fTwMHDlRZWVmjt1lXV6dgMBi2AQCA1svxQPn888+1ZMkS9ezZU2+99Zbuvfde3X///Vq2bJkkye/3S5JSU1PDvi81NTV07FRFRUXy+XyhLSMjw+mxAQCARRwPlIaGBl199dV6/PHHddVVVykvL0+TJk3S0qVLz/k2CwsLFQgEQltlZaWDEwMAANs4HijdunXT5ZdfHravd+/e2r9/vyQpLS1NklRdXR12TnV1dejYqeLj4+X1esM2AADQejkeKDfccIN27twZtm/Xrl26+OKLJX3/hNm0tDSVlpaGjgeDQZWXlys7O9vpcQAAQBRy/FU806dP1/XXX6/HH39ct99+uzZu3Khnn31Wzz77rCTJ4/Fo2rRpevTRR9WzZ09lZWVp1qxZSk9P17hx45weBwAARCHHA+Waa67R6tWrVVhYqN///vfKysrSwoULNWHChNA5M2bMUG1trfLy8lRTU6Mbb7xRJSUl6tChg9PjAACAKOQxxhi3hzhbwWBQPp9PgUCA56PgNHwWD1qDfXNHuz0C4Liz+f3NZ/EAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArBPxQJk7d648Ho+mTZsW2nfs2DFNnjxZKSkp6tixo3Jzc1VdXR3pUQAAQJSIaKBs2rRJzzzzjK688sqw/dOnT9cbb7yhVatWaf369aqqqtJtt90WyVEAAEAUiVigHDlyRBMmTNBzzz2nzp07h/YHAgE9//zzWrBggYYMGaL+/furuLhYH3zwgT788MNIjQMAAKJIxAJl8uTJGj16tHJycsL2V1RUqL6+Pmx/r169lJmZqbKyskZvq66uTsFgMGwDAACtV2wkbvTll1/Wli1btGnTptOO+f1+xcXFKSkpKWx/amqq/H5/o7dXVFSkRx55JBKjAoCVehSsa9H375s72qFJAHc4fgWlsrJSDzzwgFasWKEOHTo4cpuFhYUKBAKhrbKy0pHbBQAAdnI8UCoqKnTw4EFdffXVio2NVWxsrNavX69FixYpNjZWqampOn78uGpqasK+r7q6WmlpaY3eZnx8vLxeb9gGAABaL8cf4hk6dKg++eSTsH133XWXevXqpYceekgZGRlq3769SktLlZubK0nauXOn9u/fr+zsbKfHQRRq6aVtAED0czxQOnXqpCuuuCJsX2JiolJSUkL7J06cqPz8fCUnJ8vr9Wrq1KnKzs7Wdddd5/Q4AAAgCkXkSbLNefLJJxUTE6Pc3FzV1dVp+PDhevrpp90YBQAAWMhjjDFuD3G2gsGgfD6fAoEAz0dphXiIB2g5XsUDG53N728+iwcAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHVi3R4ArU+PgnVujwAAiHJcQQEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANbh04wBoBVq6aeK75s72qFJgHPDFRQAAGAdAgUAAFiHQAEAANZxPFCKiop0zTXXqFOnTuratavGjRunnTt3hp1z7NgxTZ48WSkpKerYsaNyc3NVXV3t9CgAACBKOR4o69ev1+TJk/Xhhx/q7bffVn19vYYNG6ba2trQOdOnT9cbb7yhVatWaf369aqqqtJtt93m9CgAACBKOf4qnpKSkrCvX3zxRXXt2lUVFRUaNGiQAoGAnn/+ea1cuVJDhgyRJBUXF6t379768MMPdd111zk9EgAAiDIRfw5KIBCQJCUnJ0uSKioqVF9fr5ycnNA5vXr1UmZmpsrKyhq9jbq6OgWDwbANAAC0XhENlIaGBk2bNk033HCDrrjiCkmS3+9XXFyckpKSws5NTU2V3+9v9HaKiork8/lCW0ZGRiTHBgAALotooEyePFmffvqpXn755RbdTmFhoQKBQGirrKx0aEIAAGCjiL2T7JQpU7R27Vpt2LBB3bt3D+1PS0vT8ePHVVNTE3YVpbq6WmlpaY3eVnx8vOLj4yM1KgAAsIzjgWKM0dSpU7V69Wq99957ysrKCjvev39/tW/fXqWlpcrNzZUk7dy5U/v371d2drbT4+ActPQtsgEAaCnHA2Xy5MlauXKlXnvtNXXq1Cn0vBKfz6eEhAT5fD5NnDhR+fn5Sk5Oltfr1dSpU5Wdnc0reAAAgKQIBMqSJUskSYMHDw7bX1xcrN/85jeSpCeffFIxMTHKzc1VXV2dhg8frqefftrpUQAAQJSKyEM8zenQoYMWL16sxYsXO/3HAwCAVoDP4gEAANYhUAAAgHUIFAAAYB0CBQAAWCdib9QGAIheLX0/pH1zRzs0CdoqrqAAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA7vgwIAcBzvo4KW4goKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDp9m3Aq19FNEAQBwG1dQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB3e6h4A0Oq09CM/9s0d7dAkOFdcQQEAANYhUAAAgHV4iKcRXBoEAMBdXEEBAADWIVAAAIB1eIjHQi19iAkAoh0/B1vGif9+bj9dwdUrKIsXL1aPHj3UoUMHDRw4UBs3bnRzHAAAYAnXAuWVV15Rfn6+5syZoy1btqhv374aPny4Dh486NZIAADAEq49xLNgwQJNmjRJd911lyRp6dKlWrdunV544QUVFBS4NRYAAC3GQ1Qt50qgHD9+XBUVFSosLAzti4mJUU5OjsrKyk47v66uTnV1daGvA4GAJCkYDEZkvoa6oy36/pbO1dI/HwDQMvwcj8zv2JO3aYxp9lxXAuV///ufTpw4odTU1LD9qamp+uyzz047v6ioSI888shp+zMyMiI2Y0v4Fro9AQCgJfg5Htn/BocPH5bP52vynKh4FU9hYaHy8/NDXzc0NOjQoUNKSUmRx+NxZaZgMKiMjAxVVlbK6/W6MoMb2uq6JdbeFtfeVtcttd21t9V1S+dn7cYYHT58WOnp6c2e60qgXHjhhWrXrp2qq6vD9ldXVystLe208+Pj4xUfHx+2LykpKZIjnjGv19vm/hJLbXfdEmtvi2tvq+uW2u7a2+q6pcivvbkrJye58iqeuLg49e/fX6WlpaF9DQ0NKi0tVXZ2thsjAQAAi7j2EE9+fr7uvPNODRgwQNdee60WLlyo2tra0Kt6AABA2+VaoPz85z/X119/rdmzZ8vv96tfv34qKSk57YmztoqPj9ecOXNOe+iptWur65ZYe1tce1tdt9R2195W1y3Zt3aPOZPX+gAAAJxHfFggAACwDoECAACsQ6AAAADrECgAAMA6BEoTlixZoiuvvDL0pjXZ2dl68803Q8ePHTumyZMnKyUlRR07dlRubu5pbz4XrZpb+7PPPqvBgwfL6/XK4/GopqbGvWEd1NS6Dx06pKlTp+qyyy5TQkKCMjMzdf/994c+GyraNXef33PPPbrkkkuUkJCgLl26aOzYsY1+NEW0aW7dJxljNHLkSHk8Hq1Zs+b8DxoBza198ODB8ng8Ydtvf/tbFyd2zpnc72VlZRoyZIgSExPl9Xo1aNAgffvtty5N7Iym1r1v377T7u+T26pVq877rARKE7p37665c+eqoqJCmzdv1pAhQzR27Fht375dkjR9+nS98cYbWrVqldavX6+qqirddtttLk/tjObWfvToUY0YMUIPP/ywy5M6q6l1V1VVqaqqSvPnz9enn36qF198USUlJZo4caLbYzuiufu8f//+Ki4u1o4dO/TWW2/JGKNhw4bpxIkTLk/eMs2t+6SFCxe69tEakXIma580aZIOHDgQ2ubNm+fixM5pbu1lZWUaMWKEhg0bpo0bN2rTpk2aMmWKYmKi+9dmU+vOyMgIu68PHDigRx55RB07dtTIkSPP/7AGZ6Vz587mT3/6k6mpqTHt27c3q1atCh3bsWOHkWTKyspcnDByTq79h959910jyXzzzTfuDHUeNLbuk1599VUTFxdn6uvrz/NU50dTa//444+NJLNnz57zPFXknbrujz76yFx00UXmwIEDRpJZvXq1e8NF2A/XftNNN5kHHnjA3YHOox+ufeDAgWbmzJkuT3R+NPXvvF+/fubuu+8+zxN9L7pT8Dw6ceKEXn75ZdXW1io7O1sVFRWqr69XTk5O6JxevXopMzNTZWVlLk7qvFPX3lacyboDgYC8Xq9iY6PiczfPWHNrr62tVXFxsbKysqz9VPFz0di6jx49qjvuuEOLFy9u9LPCWosfu89XrFihCy+8UFdccYUKCwt19OhRF6eMjFPXfvDgQZWXl6tr1666/vrrlZqaqptuuknvv/++26M6qrl/5xUVFdq6dat7V4ldyaIosm3bNpOYmGjatWtnfD6fWbdunTHGmBUrVpi4uLjTzr/mmmvMjBkzzveYEfFja/+h1ngF5UzWbYwxX3/9tcnMzDQPP/zweZ4wcppb++LFi01iYqKRZC677LJWc/WkqXXn5eWZiRMnhr5WK7uC0tTan3nmGVNSUmK2bdtmXnrpJXPRRReZW2+91cVpnfVjay8rKzOSTHJysnnhhRfMli1bzLRp00xcXJzZtWuXy1O33Jn+jLv33ntN7969z/N0/x+B0oy6ujqze/dus3nzZlNQUGAuvPBCs3379jYRKD+29h9qjYFyJusOBALm2muvNSNGjDDHjx93aVLnNbf2mpoas2vXLrN+/XozZswYc/XVV5tvv/3WxYmd8WPrfu2118yll15qDh8+HDq3tQXKmfx9P6m0tLRVPaz3Y2v/97//bSSZwsLCsPP79OljCgoKXJrWOWdynx89etT4fD4zf/58l6YkUM7a0KFDTV5eXugf6qm/mDMzM82CBQvcGS7CTq79h1pjoJzq1HUHg0GTnZ1thg4d2ip+OTelsfv8pLq6OnPBBReYlStXnuepIu/kuh944AHj8XhMu3btQpskExMTY2666Sa3x4yIpu7zI0eOGEmmpKTkPE91fpxc++eff24kmT//+c9hx2+//XZzxx13uDRd5DR2ny9fvty0b9/eHDx40KWpeA7KWWtoaFBdXZ369++v9u3bq7S0NHRs586d2r9/f6t9nsbJtbc1P1x3MBjUsGHDFBcXp9dff10dOnRwebrIauo+N9//H5xW+Xfi5LoLCgq0bds2bd26NbRJ0pNPPqni4mJ3h4yQpu7zk+vv1q3beZzo/Dm59h49eig9PV07d+4MO75r1y5dfPHFLk0XOY3d588//7xuueUWdenSxaWpXPw042hQWFiokSNHKjMzU4cPH9bKlSv13nvv6a233pLP59PEiROVn5+v5ORkeb1eTZ06VdnZ2bruuuvcHr3Fmlq7JPn9fvn9fu3Zs0eS9Mknn6hTp07KzMxUcnKym6O3SFPrPhknR48e1UsvvaRgMKhgMChJ6tKli9q1a+fy9C3T1No///xzvfLKKxo2bJi6dOmiL7/8UnPnzlVCQoJGjRrl9ugt0tS609LSGn1ibGZmprKyslyY1llNrf2///2vVq5cqVGjRiklJUXbtm3T9OnTNWjQIF155ZVuj95iTa3d4/HowQcf1Jw5c9S3b1/169dPy5Yt02effaa//vWvbo/eIs39bJekPXv2aMOGDfr73//u4qTiSbJNufvuu83FF19s4uLiTJcuXczQoUPNP/7xj9Dxb7/91tx3332mc+fO5oILLjC33nqrOXDggIsTO6e5tc+ZM8dIOm0rLi52b2gHNLXukw9nNbbt3bvX3cEd0NTav/rqKzNy5EjTtWtX0759e9O9e3dzxx13mM8++8zlqVuuub/rp1Ireg5KU2vfv3+/GTRokElOTjbx8fHm0ksvNQ8++KAJBAIuT+2MM7nfi4qKTPfu3c0FF1xgsrOzzb/+9S+XpnXOmay7sLDQZGRkmBMnTrg05fc8xhjjThoBAAA0juegAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArPP/AH77U2J8Ne9jAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "m = np.isfinite(mock[\"mu_calibration\"])\n",
    "\n",
    "plt.figure()\n",
    "plt.hist(mock[\"mu_calibration\"][m], bins=\"auto\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "31.757872992956177"
      ]
     },
     "execution_count": 101,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.percentile(mock[\"mu_true\"], 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAhHElEQVR4nO3df1BVdf7H8ReIIKL3EiQgCcmUk7ZZlhZRjbnKqOiYFrutmzvrlqNtqaXMZNCorfu1MMctRzOptkg32cqdzUo3Whdb3XYJFTPLNdRNkw0vtmPcq5hI8vn+0Xinq4Q/OHA+F5+PmTuznHO4vj9di+d+7g8ijDFGAAAAFol0ewAAAIDTESgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArBPl9gAXoqmpSTU1NerevbsiIiLcHgcAAJwDY4yOHDmi1NRURUa2vEcSloFSU1OjtLQ0t8cAAAAXoLq6Wr169WrxmrAMlO7du0v6boEej8flaQAAwLkIBAJKS0sL/hxvSVgGyqmndTweD4ECAECYOZeXZ/AiWQAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWCfK7QEA2Kd3/rpWff/+BaMdmgTAxYodFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1+KA2ANbhg+IAsIMCAACsQ6AAAADrECgAAMA65x0omzZt0pgxY5SamqqIiAitWbMmeK6xsVGPPvqo+vfvr7i4OKWmpuqXv/ylampqQu7j8OHDmjBhgjwej+Lj4zVp0iQdPXq01YsBAAAdw3kHSn19va677jotW7bsjHPHjh3Ttm3bNGfOHG3btk1//vOfVVVVpTvuuCPkugkTJmjnzp1av3691q5dq02bNmnKlCkXvgoAANChnPe7eHJycpSTk9PsOa/Xq/Xr14cce/bZZ3XTTTfpwIEDSk9P165du1RaWqotW7Zo0KBBkqSlS5dq1KhRWrRokVJTUy9gGQAAoCNp89eg+P1+RUREKD4+XpJUXl6u+Pj4YJxIUnZ2tiIjI1VRUdHsfTQ0NCgQCITcAABAx9Wmn4Ny/PhxPfroo/r5z38uj8cjSfL5fEpKSgodIipKCQkJ8vl8zd5PYWGh5s2b15ajAnBQaz/HBADabAelsbFRd999t4wxWr58eavuq6CgQH6/P3irrq52aEoAAGCjNtlBORUnX3zxhTZs2BDcPZGklJQUHTp0KOT6b7/9VocPH1ZKSkqz9xcTE6OYmJi2GBUAAFjI8R2UU3GyZ88e/e1vf1NiYmLI+aysLNXV1amysjJ4bMOGDWpqalJmZqbT4wAAgDB03jsoR48e1d69e4Nf79u3T9u3b1dCQoJ69uypn/zkJ9q2bZvWrl2rkydPBl9XkpCQoOjoaPXr108jR47U5MmTVVRUpMbGRk2bNk3jx4/nHTwAAEDSBQTK1q1b9eMf/zj4dV5eniRp4sSJ+s1vfqO3335bkjRgwICQ73v//fc1ZMgQSdKqVas0bdo0DRs2TJGRkcrNzdWSJUsucAkAAKCjOe9AGTJkiIwxP3i+pXOnJCQkqKSk5Hz/aAAAcJHgd/EAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA60S5PQAA5/XOX+f2CADQKuygAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDp81D1gIT6qHsDFjh0UAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdc47UDZt2qQxY8YoNTVVERERWrNmTch5Y4zmzp2rnj17KjY2VtnZ2dqzZ0/INYcPH9aECRPk8XgUHx+vSZMm6ejRo61aCAAA6DjOO1Dq6+t13XXXadmyZc2eX7hwoZYsWaKioiJVVFQoLi5OI0aM0PHjx4PXTJgwQTt37tT69eu1du1abdq0SVOmTLnwVQAAgA4l6ny/IScnRzk5Oc2eM8Zo8eLFmj17tsaOHStJWrlypZKTk7VmzRqNHz9eu3btUmlpqbZs2aJBgwZJkpYuXapRo0Zp0aJFSk1NbcVyAABAR+Doa1D27dsnn8+n7Ozs4DGv16vMzEyVl5dLksrLyxUfHx+ME0nKzs5WZGSkKioqmr3fhoYGBQKBkBsAAOi4HA0Un88nSUpOTg45npycHDzn8/mUlJQUcj4qKkoJCQnBa05XWFgor9cbvKWlpTk5NgAAsExYvIunoKBAfr8/eKuurnZ7JAAA0IYcDZSUlBRJUm1tbcjx2tra4LmUlBQdOnQo5Py3336rw4cPB685XUxMjDweT8gNAAB0XI4GSkZGhlJSUlRWVhY8FggEVFFRoaysLElSVlaW6urqVFlZGbxmw4YNampqUmZmppPjAACAMHXe7+I5evSo9u7dG/x637592r59uxISEpSenq4ZM2Zo/vz56tOnjzIyMjRnzhylpqZq3LhxkqR+/fpp5MiRmjx5soqKitTY2Khp06Zp/PjxvIMHAABIuoBA2bp1q3784x8Hv87Ly5MkTZw4Ua+88opmzZql+vp6TZkyRXV1dbrttttUWlqqLl26BL9n1apVmjZtmoYNG6bIyEjl5uZqyZIlDiwHAAB0BBHGGOP2EOcrEAjI6/XK7/fzehR0SL3z17k9Qljbv2C02yMAaMb5/PwOi3fxAACAiwuBAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOtEuT0A0BH1zl/n9ggXtdb+89+/YLRDkwC4UOygAAAA6xAoAADAOgQKAACwDoECAACsw4tkAeA0vMgWcB87KAAAwDoECgAAsA6BAgAArEOgAAAA6zgeKCdPntScOXOUkZGh2NhYXXHFFfq///s/GWOC1xhjNHfuXPXs2VOxsbHKzs7Wnj17nB4FAACEKccD5amnntLy5cv17LPPateuXXrqqae0cOFCLV26NHjNwoULtWTJEhUVFamiokJxcXEaMWKEjh8/7vQ4AAAgDDn+NuN//etfGjt2rEaP/u5tdr1799Yf//hHbd68WdJ3uyeLFy/W7NmzNXbsWEnSypUrlZycrDVr1mj8+PFOjwQAAMKM4zsot9xyi8rKyrR7925J0scff6wPPvhAOTk5kqR9+/bJ5/MpOzs7+D1er1eZmZkqLy9v9j4bGhoUCARCbgAAoONyfAclPz9fgUBAffv2VadOnXTy5Ek98cQTmjBhgiTJ5/NJkpKTk0O+Lzk5OXjudIWFhZo3b57TowIAAEs5voPyxhtvaNWqVSopKdG2bdu0YsUKLVq0SCtWrLjg+ywoKJDf7w/eqqurHZwYAADYxvEdlEceeUT5+fnB15L0799fX3zxhQoLCzVx4kSlpKRIkmpra9WzZ8/g99XW1mrAgAHN3mdMTIxiYmKcHhUAAFjK8R2UY8eOKTIy9G47deqkpqYmSVJGRoZSUlJUVlYWPB8IBFRRUaGsrCynxwEAAGHI8R2UMWPG6IknnlB6erp+9KMf6aOPPtLTTz+t++67T5IUERGhGTNmaP78+erTp48yMjI0Z84cpaamaty4cU6PAwAAwpDjgbJ06VLNmTNHDz74oA4dOqTU1FTdf//9mjt3bvCaWbNmqb6+XlOmTFFdXZ1uu+02lZaWqkuXLk6PAwAAwlCE+f5HvIaJQCAgr9crv98vj8fj9jjAGXrnr3N7BLho/4LRbo8AWOl8fn7zu3gAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGCdKLcHAGzTO3+d2yMAwEWPHRQAAGAdAgUAAFiHp3gAwGFOPE24f8FoByYBwhc7KAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsw+egoMPho+oBIPyxgwIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA67RJoHz55Zf6xS9+ocTERMXGxqp///7aunVr8LwxRnPnzlXPnj0VGxur7Oxs7dmzpy1GAQAAYcjxQPn666916623qnPnznr33Xf173//W7/73e90ySWXBK9ZuHChlixZoqKiIlVUVCguLk4jRozQ8ePHnR4HAACEIcc/qO2pp55SWlqaiouLg8cyMjKC/9sYo8WLF2v27NkaO3asJGnlypVKTk7WmjVrNH78eKdHAgAAYcbxHZS3335bgwYN0k9/+lMlJSXp+uuv14svvhg8v2/fPvl8PmVnZwePeb1eZWZmqry83OlxAABAGHI8UD7//HMtX75cffr00XvvvacHHnhADz30kFasWCFJ8vl8kqTk5OSQ70tOTg6eO11DQ4MCgUDIDQAAdFyOP8XT1NSkQYMG6cknn5QkXX/99fr0009VVFSkiRMnXtB9FhYWat68eU6OCQAALOb4DkrPnj119dVXhxzr16+fDhw4IElKSUmRJNXW1oZcU1tbGzx3uoKCAvn9/uCturra6bEBAIBFHA+UW2+9VVVVVSHHdu/ercsvv1zSdy+YTUlJUVlZWfB8IBBQRUWFsrKymr3PmJgYeTyekBsAAOi4HH+KZ+bMmbrlllv05JNP6u6779bmzZv1wgsv6IUXXpAkRUREaMaMGZo/f7769OmjjIwMzZkzR6mpqRo3bpzT4wAAgDDkeKDceOONevPNN1VQUKDf/va3ysjI0OLFizVhwoTgNbNmzVJ9fb2mTJmiuro63XbbbSotLVWXLl2cHgcAAIShCGOMcXuI8xUIBOT1euX3+3m6B2fonb/O7RGAVtu/YLTbIwCOO5+f3/wuHgAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHWi3B4AAHCm3vnrWvX9+xeMdmgSwB3soAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACs0+aBsmDBAkVERGjGjBnBY8ePH9fUqVOVmJiobt26KTc3V7W1tW09CgAACBNtGihbtmzR888/r2uvvTbk+MyZM/XOO+9o9erV2rhxo2pqanTXXXe15SgAACCMtNlH3R89elQTJkzQiy++qPnz5weP+/1+vfTSSyopKdHQoUMlScXFxerXr58+/PBD3XzzzW01EsJEaz/iGwAQ/tpsB2Xq1KkaPXq0srOzQ45XVlaqsbEx5Hjfvn2Vnp6u8vLythoHAACEkTbZQXnttde0bds2bdmy5YxzPp9P0dHRio+PDzmenJwsn8/X7P01NDSooaEh+HUgEHB0XgAAYBfHd1Cqq6v18MMPa9WqVerSpYsj91lYWCiv1xu8paWlOXK/AADATo4HSmVlpQ4dOqQbbrhBUVFRioqK0saNG7VkyRJFRUUpOTlZJ06cUF1dXcj31dbWKiUlpdn7LCgokN/vD96qq6udHhsAAFjE8ad4hg0bpk8++STk2L333qu+ffvq0UcfVVpamjp37qyysjLl5uZKkqqqqnTgwAFlZWU1e58xMTGKiYlxelQAAGApxwOle/fuuuaaa0KOxcXFKTExMXh80qRJysvLU0JCgjwej6ZPn66srCzewQMAACS14duMW/LMM88oMjJSubm5amho0IgRI/Tcc8+5MQoAALBQhDHGuD3E+QoEAvJ6vfL7/fJ4PG6PA4fxOShA6+1fMNrtEYAznM/Pb34XDwAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALBOlNsDoOPpnb/O7REAAGGOHRQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh4+6xxn4qHog/LX23+P9C0Y7NAlwYRzfQSksLNSNN96o7t27KykpSePGjVNVVVXINcePH9fUqVOVmJiobt26KTc3V7W1tU6PAgAAwpTjgbJx40ZNnTpVH374odavX6/GxkYNHz5c9fX1wWtmzpypd955R6tXr9bGjRtVU1Oju+66y+lRAABAmHL8KZ7S0tKQr1955RUlJSWpsrJSgwcPlt/v10svvaSSkhINHTpUklRcXKx+/frpww8/1M033+z0SAAAIMy0+Ytk/X6/JCkhIUGSVFlZqcbGRmVnZwev6du3r9LT01VeXt7sfTQ0NCgQCITcAABAx9WmgdLU1KQZM2bo1ltv1TXXXCNJ8vl8io6OVnx8fMi1ycnJ8vl8zd5PYWGhvF5v8JaWltaWYwMAAJe1aaBMnTpVn376qV577bVW3U9BQYH8fn/wVl1d7dCEAADARm32NuNp06Zp7dq12rRpk3r16hU8npKSohMnTqiuri5kF6W2tlYpKSnN3ldMTIxiYmLaalQAAGAZx3dQjDGaNm2a3nzzTW3YsEEZGRkh5wcOHKjOnTurrKwseKyqqkoHDhxQVlaW0+MAAIAw5PgOytSpU1VSUqK33npL3bt3D76uxOv1KjY2Vl6vV5MmTVJeXp4SEhLk8Xg0ffp0ZWVl8Q4eAAAgqQ0CZfny5ZKkIUOGhBwvLi7Wr371K0nSM888o8jISOXm5qqhoUEjRozQc8895/QoAAAgTEUYY4zbQ5yvQCAgr9crv98vj8fj9jgdDh91D6C1+Kh8NOd8fn7zywIBAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdaLcHgDO652/zu0RAABoFXZQAACAdQgUAABgHZ7iAQA4rrVPNe9fMNqhSRCu2EEBAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdfioewCAdfiofLCDAgAArEOgAAAA6/AUDwCgw+EpovDHDgoAALAOgQIAAKzDUzzNYGsQAAB3sYMCAACsQ6AAAADr8BSPhVr7FBMAILzxUgOXd1CWLVum3r17q0uXLsrMzNTmzZvdHAcAAFjCtR2U119/XXl5eSoqKlJmZqYWL16sESNGqKqqSklJSW6NBQBAq7m9E+7En+/2LoxrOyhPP/20Jk+erHvvvVdXX321ioqK1LVrV7388stujQQAACzhyg7KiRMnVFlZqYKCguCxyMhIZWdnq7y8/IzrGxoa1NDQEPza7/dLkgKBQJvM19RwrFXf39q5WvvnAwBah/+Ot83P2FP3aYw567WuBMr//vc/nTx5UsnJySHHk5OT9dlnn51xfWFhoebNm3fG8bS0tDabsTW8i92eAADQGvx3vG3/GRw5ckRer7fFa8LiXTwFBQXKy8sLft3U1KTDhw8rMTFRERER7TpLIBBQWlqaqqur5fF42vXPdhtrZ+2s/eLB2ll7W6zdGKMjR44oNTX1rNe6EiiXXnqpOnXqpNra2pDjtbW1SklJOeP6mJgYxcTEhByLj49vyxHPyuPxXHR/cU9h7az9YsPaWfvFpi3Xfradk1NceZFsdHS0Bg4cqLKysuCxpqYmlZWVKSsry42RAACARVx7iicvL08TJ07UoEGDdNNNN2nx4sWqr6/Xvffe69ZIAADAEq4Fys9+9jN99dVXmjt3rnw+nwYMGKDS0tIzXjhrm5iYGD3++ONnPOV0MWDtrP1iw9pZ+8XGprVHmHN5rw8AAEA74pcFAgAA6xAoAADAOgQKAACwDoECAACsQ6D8gOXLl+vaa68NflhNVlaW3n333eD548ePa+rUqUpMTFS3bt2Um5t7xgfPhaOzrfuFF17QkCFD5PF4FBERobq6OveGdVhLaz98+LCmT5+uq666SrGxsUpPT9dDDz0U/L1Q4e5sj/v999+vK664QrGxserRo4fGjh3b7K+lCEdnW/spxhjl5OQoIiJCa9asaf9B28DZ1j5kyBBFRESE3H7961+7OLFzzuVxLy8v19ChQxUXFyePx6PBgwfrm2++cWli57S09v3795/xmJ+6rV69ul3nJFB+QK9evbRgwQJVVlZq69atGjp0qMaOHaudO3dKkmbOnKl33nlHq1ev1saNG1VTU6O77rrL5alb72zrPnbsmEaOHKnHHnvM5Umd19Laa2pqVFNTo0WLFunTTz/VK6+8otLSUk2aNMntsR1xtsd94MCBKi4u1q5du/Tee+/JGKPhw4fr5MmTLk/eemdb+ymLFy9u91+t0dbOZe2TJ0/WwYMHg7eFCxe6OLFzzrb28vJyjRw5UsOHD9fmzZu1ZcsWTZs2TZGR4f9js6W1p6WlhTzeBw8e1Lx589StWzfl5OS076AG5+ySSy4xv//9701dXZ3p3LmzWb16dfDcrl27jCRTXl7u4oRt49S6v+/99983kszXX3/tzlDtpLm1n/LGG2+Y6Oho09jY2M5TtY+W1v7xxx8bSWbv3r3tPFX7OH3tH330kbnsssvMwYMHjSTz5ptvujdcG/v+2m+//Xbz8MMPuztQO/r+2jMzM83s2bNdnqj9tPTv+4ABA8x9993XzhMZE/4p2A5Onjyp1157TfX19crKylJlZaUaGxuVnZ0dvKZv375KT09XeXm5i5M66/R1X0zOZe1+v18ej0dRUWHxOzfP2dnWXl9fr+LiYmVkZFj7G8UvVHNrP3bsmO655x4tW7as2d8V1lH80OO+atUqXXrppbrmmmtUUFCgY8eOuThl2zh97YcOHVJFRYWSkpJ0yy23KDk5Wbfffrs++OADt0d13Nn+fa+srNT27dvd2S1u9yQKIzt27DBxcXGmU6dOxuv1mnXr1hljjFm1apWJjo4+4/obb7zRzJo1q73HdNwPrfv7OuoOyrms3RhjvvrqK5Oenm4ee+yxdp6w7Zxt7cuWLTNxcXFGkrnqqqs61O5JS2ufMmWKmTRpUvBrdbAdlJbW/vzzz5vS0lKzY8cO8+qrr5rLLrvM3HnnnS5O66wfWnt5ebmRZBISEszLL79stm3bZmbMmGGio6PN7t27XZ7aGef637oHHnjA9OvXr52n+w6B0oKGhgazZ88es3XrVpOfn28uvfRSs3Pnzg4fKD+07u/rqIFyLmv3+/3mpptuMiNHjjQnTpxwaVLnnW3tdXV1Zvfu3Wbjxo1mzJgx5oYbbjDffPONixM754fW/tZbb5krr7zSHDlyJHhtRwuUc/k7f0pZWVmHemrvh9b+z3/+00gyBQUFIdf379/f5OfnuzSts87lcT927Jjxer1m0aJFrsxIoJyHYcOGmSlTpgT/JT39h3N6erp5+umn3RmuDZ1a9/d11EA53elrDwQCJisrywwbNqzD/HD+Ic097qc0NDSYrl27mpKSknaeqn2cWvvDDz9sIiIiTKdOnYI3SSYyMtLcfvvtbo/ZJlp63I8ePWokmdLS0naeqn2cWvvnn39uJJk//OEPIefvvvtuc88997g0Xdtq7nFfuXKl6dy5szl06JArM/EalPPQ1NSkhoYGDRw4UJ07d1ZZWVnwXFVVlQ4cONAhX6txat0Xo++vPRAIaPjw4YqOjtbbb7+tLl26uDxd22rpcTff/Z+bDvv34tTa8/PztWPHDm3fvj14k6RnnnlGxcXF7g7ZRlp63E+tv2fPnu04Ufs5tfbevXsrNTVVVVVVIed3796tyy+/3KXp2lZzj/tLL72kO+64Qz169HBlpo716j4HFRQUKCcnR+np6Tpy5IhKSkr097//Xe+99568Xq8mTZqkvLw8JSQkyOPxaPr06crKytLNN9/s9uit0tK6Jcnn88nn82nv3r2SpE8++UTdu3dXenq6EhIS3By91Vpa+6k4OXbsmF599VUFAgEFAgFJUo8ePdSpUyeXp2+dltb++eef6/XXX9fw4cPVo0cP/fe//9WCBQsUGxurUaNGuT16q7W09pSUlGZfGJuenq6MjAwXpnVWS2v/z3/+o5KSEo0aNUqJiYnasWOHZs6cqcGDB+vaa691e/RWa2ntEREReuSRR/T444/ruuuu04ABA7RixQp99tln+tOf/uT26K12tv/OS9LevXu1adMm/eUvf3FvUFf2bcLAfffdZy6//HITHR1tevToYYYNG2b++te/Bs9/88035sEHHzSXXHKJ6dq1q7nzzjvNwYMHXZzYGWdb9+OPP24knXErLi52b2iHtLT2U09pNXfbt2+fu4M7oKW1f/nllyYnJ8ckJSWZzp07m169epl77rnHfPbZZy5P7Yyz/Z0/nTrQa1BaWvuBAwfM4MGDTUJCgomJiTFXXnmleeSRR4zf73d5amecy+NeWFhoevXqZbp27WqysrLMP/7xD5emdda5rL2goMCkpaWZkydPujSlMRHGGONOGgEAADSP16AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACs8/8fu0/PYQGRKAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.hist(mock[\"mu_true\"], bins=\"auto\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "from jax import numpy as jnp\n",
    "from jax.scipy.special import logsumexp\n",
    "\n",
    "def normal_logpdf(x, loc, scale):\n",
    "    \"\"\"Log of the normal probability density function.\"\"\"\n",
    "    return (-0.5 * ((x - loc) / scale)**2\n",
    "            - jnp.log(scale) - 0.5 * jnp.log(2 * jnp.pi))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu_true = np.copy(mock[\"mu_true\"])\n",
    "mu_calibration = np.copy(mock[\"mu_calibration\"])\n",
    "e_mu = np.copy(mock[\"e_mu_calibration\"])\n",
    "\n",
    "mu_calibration = np.stack([mu_calibration, mu_calibration])\n",
    "e_mu_calibration = np.stack([e_mu, e_mu])\n",
    "\n",
    "mu_calibration[0, 100:] = np.nan\n",
    "e_mu_calibration[0, 100:] = np.nan\n",
    "\n",
    "mu_calibration[1, 50:] = np.nan\n",
    "e_mu_calibration[1, 50:] = np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "h = 0.7\n",
    "\n",
    "# Now, the rest of the code except the calibration likelihood\n",
    "# uses the distance modulus in units of h\n",
    "mu_true_h = mu_true + 5 * jnp.log10(h)\n",
    "\n",
    "# Calculate the log-likelihood of the calibration, but the\n",
    "# shape is `(n_calibrators, n_data)`.\n",
    "ll_calibration = normal_logpdf(\n",
    "    mu_calibration, mu_true[None, :],\n",
    "    e_mu_calibration)\n",
    "\n",
    "# Create a mask for valid (non-NaN) log-likelihoods\n",
    "calibration_mask = ~jnp.isnan(ll_calibration)\n",
    "\n",
    "# Replace NaN values with zero (or another neutral value) for safety\n",
    "ll_calibration_clean = jnp.where(calibration_mask, ll_calibration, 0.0)\n",
    "\n",
    "# Count the number of valid calibrators for each galaxy (non-NaN entries)\n",
    "counts = jnp.sum(calibration_mask, axis=0)\n",
    "\n",
    "# Now apply logsumexp only to the valid log-likelihoods\n",
    "ll_calibration_sum = jnp.where(\n",
    "    counts > 0,\n",
    "    logsumexp(ll_calibration_clean, axis=0) - jnp.log(counts),\n",
    "    0.0  # Return zero likelihood if no valid calibrators\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {},
   "outputs": [],
   "source": [
    "from jax.lax import cond\n",
    "def ll_calibration(mu_calibration, mu_true, e_mu):\n",
    "    # Use jnp.where to apply element-wise conditional logic\n",
    "    return jnp.where(\n",
    "        jnp.isfinite(mu_calibration),  # Check for finite values\n",
    "        normal_logpdf(mu_calibration, mu_true, e_mu),  # Use valid values\n",
    "        0.0  # Return 0 for invalid (non-finite) values\n",
    "    )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "33.15014277245038"
      ]
     },
     "execution_count": 121,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mu_calibration[0, 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Array([[1.5757915, 1.3846258, 2.0763617, ..., 0.       , 0.       ,\n",
       "        0.       ],\n",
       "       [1.5757915, 1.3846258, 2.0763617, ..., 0.       , 0.       ,\n",
       "        0.       ]], dtype=float32)"
      ]
     },
     "execution_count": 127,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ll_calibration(mu_calibration, mu_true[None, :], e_mu)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([33.15014277, 33.36291487, 31.68284461, ...,         nan,\n",
       "                nan,         nan]),\n",
       " array([33.10009268, 33.30408597, 31.68137438, ..., 31.3564011 ,\n",
       "        32.63194483, 33.72852162]))"
      ]
     },
     "execution_count": 125,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mu_calibration[0], mu_true"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1.5757915 1.3846258 2.0763617 ...       nan       nan       nan]\n",
      " [1.5757915 1.3846258 2.0763617 ...       nan       nan       nan]]\n",
      "\n",
      "[[1.5757915 1.3846258 2.0763617 ...      -inf      -inf      -inf]\n",
      " [1.5757915 1.3846258 2.0763617 ...      -inf      -inf      -inf]]\n",
      "\n",
      "[1.5757914 1.3846259 2.0763617 ... 0.        0.        0.       ]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "ll = normal_logpdf(mu_calibration, mu_true[None, :], e_mu)\n",
    "\n",
    "print(ll)\n",
    "print()\n",
    "\n",
    "\n",
    "mask = ~jnp.isnan(ll)\n",
    "ll = jnp.where(mask, ll, -jnp.inf)\n",
    "\n",
    "print(ll)\n",
    "print()\n",
    "\n",
    "counts = jnp.sum(mask, axis=0)\n",
    "ll = jnp.where(counts > 0, logsumexp(ll, axis=0) - jnp.log(counts), 0.)\n",
    "\n",
    "print(ll)\n",
    "print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Array(100, dtype=int32)"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "jnp.sum(~jnp.isnan(ll))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "venv_csiborg",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}