notebooks/meps-cqr.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/aangelopoulos/conformal-prediction/blob/main/notebooks/meps-cqr.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import json\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.image import imread\n",
    "import seaborn as sns\n",
    "!pip install -U --no-cache-dir gdown --pre"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load cached data\n",
    "if not os.path.exists('../data'):\n",
    "    os.system('gdown 1h7S6N_Rx7gdfO3ZunzErZy6H7620EbZK -O ../data.tar.gz')\n",
    "    os.system('tar -xf ../data.tar.gz -C ../')\n",
    "    os.system('rm ../data.tar.gz')\n",
    "    \n",
    "data = np.load('../data/meps/meps-gbr.npz')\n",
    "X, labels, upper, lower = data['X'], data['y'], data['upper'], data['lower']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Problem setup\n",
    "n=1000 # number of calibration points\n",
    "alpha = 0.1 # 1-alpha is the desired coverage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Split the softmax scores into calibration and validation sets (save the shuffling)\n",
    "idx = np.array([1] * n + [0] * (labels.shape[0]-n)) > 0\n",
    "np.random.shuffle(idx)\n",
    "cal_labels, val_labels = labels[idx], labels[~idx]\n",
    "cal_upper, val_upper = upper[idx], upper[~idx]\n",
    "cal_lower, val_lower = lower[idx], lower[~idx]\n",
    "val_X = X[~idx]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conformal prediction happens here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get scores. cal_upper.shape[0] == cal_lower.shape[0] == n\n",
    "cal_scores = np.maximum(cal_labels-cal_upper, cal_lower-cal_labels)\n",
    "# Get the score quantile\n",
    "qhat = np.quantile(cal_scores, np.ceil((n+1)*(1-alpha))/n, interpolation='higher')\n",
    "# Deploy (output=lower and upper adjusted quantiles)\n",
    "prediction_sets = [val_lower - qhat, val_upper + qhat]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The empirical coverage before calibration is: 0.7378499927462643\n",
      "The empirical coverage after calibration is: 0.9305092122443058\n"
     ]
    }
   ],
   "source": [
    "# Calculate empirical coverage (before and after calibration)\n",
    "prediction_sets_uncalibrated = [val_lower, val_upper]\n",
    "empirical_coverage_uncalibrated = ((val_labels >= prediction_sets_uncalibrated[0]) & (val_labels <= prediction_sets_uncalibrated[1])).mean()\n",
    "print(f\"The empirical coverage before calibration is: {empirical_coverage_uncalibrated}\")\n",
    "empirical_coverage = ((val_labels >= prediction_sets[0]) & (val_labels <= prediction_sets[1])).mean()\n",
    "print(f\"The empirical coverage after calibration is: {empirical_coverage}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEQCAYAAABfiGi4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOydeXxU5fm3rzNbQvYEskFYA4QoIksA/ckaKlhKCW4F0WoVi1oX1KJSUAioCMgrVSsVq9UiKMq+qMVW9lBRhAghLAIhkBWykQUy+/vHZCaZrJNllsx5rn4smTMz59xz5szzPffy3I9kNpvNCAQCgUAAKNxtgEAgEAg8ByEKAoFAILAhREEgEAgENoQoCAQCgcCGEAWBQCAQ2BCiIBAIBAIbLhOF3bt3M2XKFJKSkpg8eTLffvstABkZGUydOpUJEyYwdepULly44CqTBAKBQFALyRXzFMxmM8OGDWPt2rX07duXU6dOcd999/HTTz/xhz/8gbvvvpukpCS2bt3Kxo0bWb16tbNNEggEAkE9uMxTUCgUlJWVAVBWVkZERATFxcWkp6czadIkACZNmkR6ejpFRUWuMksgEAgENVC54iCSJPHXv/6VP/3pT/j5+VFRUcEHH3xAbm4ukZGRKJVKAJRKJREREeTm5hIWFtbkfg0GA3l5eURFRaFSueSjCAQCgVfjkpHUYDCwatUqVq5cyZAhQ/jpp5949tlnWbZsWav2m5eXx7hx4/juu++IiYlpI2vbnrfffpuysjKeeeYZgoKC3G2OR/Daa68B8PLLL7vZEs/Aej6effZZAgIC3GyN+zl9+jTr16/H19eX2bNnu9sct3PhwgXWrFkDwIABA5g8ebLTjuWS8NHJkye5fPkyQ4YMAWDIkCF06NABHx8f8vPzMRqNABiNRi5fvkx0dLQrzBIIBIJ2Qc3Ur7PTwC4RhaioKPLy8jh//jwA586do7CwkO7duxMfH8+OHTsA2LFjB/Hx8Q6FjtojovegQOAY4rdijytFwSXho/DwcJKTk5k1axaSJAGwePFiQkJCSE5OZs6cOaxcuZKgoCCWLl3qCpPcgvWzCwQCQXPwOlEAmDx5cr1xsNjYWNavX+8qMwQCQTtAeAr2mEwm299eIwoCgRXxg28YTz83JpOJrKwsKioqnHocvV7Pr371KyRJ4uTJk049VnvAej4A1Gq1w+fE39+fmJgYFArHMwVCFAQCgcMUFBQgSRJxcXHNGmiay/Xr1ykuLkahUBAVFeW047QXKisrbfO3fH19Hcq7mkwmsrOzKSgoICIiwuFjid5HAoHAYUpKSoiMjHSqIAjqYvMgm5GXVCgUREZGcvXq1WYdS3yzAoHAYYxGI2q12t1myJbmlqqo1WoMBkOz3iNEQSAQNAtRRedempN3asl3JURBIPAgPD3RLFcOHTrEXXfdBUB+fj6///3vbc/FxcW1SeL9v//9L8eOHav3uabCR5s2beKZZ55ptQ0gREEgEAiaRWRkJJ9++mmz39dUGKcxUbDiCh9NiILA5Yi7YUFTOHqNZGVlMXz48DqPrf+uWLGCKVOmMGHCBA4fPmx73e7du7nrrruYPHkyU6ZM4dSpUwD8+c9/5q677uK3v/0tTz75ZL1J2trHBPjoo49ISkpiwoQJ7Ny507Y9Li6Od999l7vvvpu//e1vnD59munTp3PnnXcyceJEPvnkEwD279/Prl27+OCDD0hKSmLLli0AbN68mXvvvZfp06fz/PPPc/HiJQB0Oh3z589n/PjxTJ06tUkxaQ6iJFUgEHglJSUlDBw4kOeee45t27axfPly1q1bR0ZGBi+//DJr166lR48e6HQ6dDodAPPmzbOVe65YsYJ//OMfDjXkUygUbN26lfPnz3PfffeRkJBAx44dAfDx8WHjxo0AlJeX88knn6DRaKioqODee+9l5MiRjBw5ksTERPr3788DDzwAwOHDh/nmm29Yu3YtBoOBnTt3snTpElatWsUXX3xBVlYWX331FQaDgfvvv7/NmoIKUXAh4g5Z4FUcS4a0hU7ZtV/Vf+U9/gxRy1u2Dz8/xo4dC8DAgQNtLXQOHjzIqFGj6NGjBwAajQaNRgPA1q1b2b59O3q9nmvXrtle0xT33nsvAL169eKGG24gNTWVcePGAXDnnXfaXldZWUlycjKnT59GkiQuX77MqVOniI2NrbPPXbt2cerUKe69915MJhMGg4Hy8nLAkuOYMmUKarUatVrN5MmTOXLkSPNPUj0IUXAhompD0BTt6sZhQLLlPydw7do1SkpKUCgUNNZIXKVS2Z0zrVZr+9s60IPlTr6pmP7hw4f5/PPPWbduHWFhYWzfvp0vv/yyxZ/Bip+fn+3vt956i/DwcJYsWYJKpeKRRx6xs7kmZrOZu+++m1mzZlFRUcHVq1dRKp0/ZIucgkAgaLd06tQJvV5PZmYmgK3jcmPcdttt7Nu3z7YevE6no7y8nNLSUgICAggJCUGn09lCPo5gfe2FCxdIT09n4MCB9b6urKzMtijYmTNn7PIcAQEBttUpARITE9m6dSt5eXmAZY7I6TOnAbjlllvYunUrBoOByspKhz63owhPQSDwINqVp+ABqFQq5s2bx8MPP0xYWBhjxoxp8j09evTg1Vdf5bnnnsNoNKJUKlmyZAkjR45k27ZtTJgwgdDQUBISEjh+/LhDdhiNRqZMmcL169dZtGiRLZ9QmyeeeIIXX3yRDRs20LNnT4YOHWp7bvLkyfzlL3/h3//+Nw8//DBTpkzh2Wef5YknnkCv16PVahk7NpH4fv343e9+x+nTp5k4cSKhoaHcdNNNFBYWOmRrU0jmdnwVZmVliZXX2iEmk4nFixcDYuU1K9aV15566ilCQkLcbE3DnDx5kvj4eKcfp2b4SPQ+wubJKFUaFJKZ8PBwh9/b3O9MhI8EAoGgneCKtKQQBYHLacfOqdMR58aCOA8NIXnHcpwCgUAgaAOEpyAQCAQCq3cguUAVXFJ9lJWVxZNPPml7XFZWRnl5OT/88AMZGRnMmTOHkpISQkJCWLp0qcMTRgQCb0OETQRN4ewrxCWiEBMTw9atW22PX3/9dYxGIwALFixg+vTpJCUlsXXrVubPn8/q1atdYZbATYiBTyBoIZIE3pZT0Ol0bN++nbvvvpvCwkLS09OZNGkSAJMmTSI9Pd227Jy3IQZDgUDQEsxms8s6IrhcFHbt2kVkZCQ33ngjubm5REZGolQqAVAqlURERJCbm+tqs1yCaHMhELiX2h1Ok5KSqKysBCwziM+cOdPqYxw6dIgDBw60+L3WdRvchctFYePGjdx9992uPqzAgxAeU8OIc+Natm7diq+vb7Pe01QPpR9++IGUlJTWmNUgkgtKUl3a5iI/P58ff/yRZcuWARAdHU1+fr5tqrnRaOTy5ctER0e70iyBQOBhNGfgO3r0KMuWLbOtfvbiiy+SkpLCDz/8gF6vJzQ0lMWLF9OlS5c6742Li+PIkSP4+/sDsG3bNg4ePEhZWRkPPfSQrY11YmIiEydO5Pvvv6dv374899xzPP/881RUVKDVahk9ejQvvvgip0+fZt26dZhMJg4ePMhvfvMbZs6cyd69e/n73/+OTqdDrVbzl7/8xdYfacWKFXz99dcEBQUxbNiwRs6H5JKSVJeKwubNmxk9ejShoaEAdOzYkfj4eHbs2EFSUhI7duwgPj7e1s9cIJAb7clTSN6TzMK9zmmdbeXPQ/7M8kkNt84uKSnhqaee4t1332Xw4MEYjUbKy8u54YYbeOmllwBYv349y5cvZ8WKFU0er7CwkE2bNlFQUMCUKVNISEigX79+gKXVxIYNGwBLN9b3338ff39/9Ho9M2bMYN++fYwaNYpp06Zx7do12/EvXrzIypUr+eijjwgICOCXX37hj3/8I3v27GHXrl3s2rWLLVu24Ovra1elWQcXRZ9dLgrz5s2z25acnMycOXNYuXIlQUFBtp7nAoHAs0kek0zymGSn7Nva60ehaDzCnZqaSmxsLIMHDwYsecng4GC2bNnCZ599xrVr15oM99TknnvuASzdV8eMGcMPP/xgE4UpU6bYXmc0Glm2bBlHjx7FbDZTUFDAqVOnGDVqVJ197t+/n4sXL3L//ffbthkMBgoKCjh06BATJ060eSr33HMPK1eubNRGrwof1VymzkpsbCzr1693pRkCN9Oe7oZdjTg3rSc7O5s33niDDRs20LVrV44cOeLQ6mlNUXNdhI8//pjS0lLWr1+Pj48Pr7zySoPrIgCMHDnSFjZvKVKN/3cmYkazQCBotwwcOJBz585x9OhRwHIHn5ubi1qtJjw8HJPJxLp16xze3+bNmwEoKipi7969ddZitlJWVkZ4eDg+Pj7k5+fz3Xff2Z6rvS7Cbbfdxv79+/nll19s26xrKt9yyy188803XLt2DaPR2OAaDracggsQ6ykIBB6E8BSaR0hICO+++y5Llizh2rVrKBQKXnrpJe644w7bWgOjR4+2W8ymMUJDQ7nrrrsoKyvjscceIy4urt7X/f73v2fWrFlMmjSJyMhIbr31Vttzv/rVr9iyZQtJSUm2RPObb77JvHnzqKysRK/XM3jwYAYMGMDYsWNJTU0lKSnJlmjOz89v0D7L3DXnXiNiPQUXINZTsEen09lcabGeggXregozZ84kIiLCzdY0jKvWU6iZUxDrKVgS6pVaLWpNB7TXy4mOjnZ43pNYT0Hg8bTj+xCBixDXiD2uPB9CFAQCD0IMhoKG8dI2FwKBGPjaN+L7cw+S7f8c/w5a8l0JURAIPAhPH3CVSiV6vd7dZsgOs9lsrwoOotfrUamaV08kREEgEDhMSEgI+fn5mEwmd5siQ5q3xI7JZCI/P5/g4OBmHUWUpApcjqffDbsTTz83nTp1Iisri9OnTzv1OJWVlbbJYMXFxU49VnugoqICo8mEUqXGoNNSXFzc5GxvAH9/fzp16tSsYwlREAgEDqNQKOjWrZvTj7N7925SUlLQaDS8+OKLTj+ep7N27VoKSyuJ6hHPmZ92MWvWLAIDA51yLBE+ErgcT78bdifi3FiwngdxPiyYTCYkSbLNTXBm+E6IgkAg8DiEKNhjMpmQFAqQFLbHzkKIgsDliB+6oCmEKNhj8xSq8gjOPC9CFAQCD0IMghaEKNhjMpuRJIUIHwm8E/FDFzSFEAV7LJ6CAskF4aNGq48MBgO7du1iz549nDp1irKyMgIDA+nXrx+jRo3iV7/6VbMnRsgZcYELmkJcIxaEKNhjMpmQ1NWegjPPS4Mj+ueff86qVauIjY1l6NChjB07Fn9/fyoqKjh37hzr169nyZIlPPbYY9x3331OM9CbcLSrobcjfuiCpqh5jZjNZtn/dmrnFNziKVy8eJH169cTHh5e57nbb7+dxx9/nMuXL/Pxxx87zTiBQG4IwbQgzoM9ZpMJhYtyCg2KgnXR6caIiIhw6HUCQU3ED17QFMJTsMdkMqOSJJfkFJqdaDaZTJw4cYLc3NxmvU+r1bJgwQLGjx/Pb3/7W1555RUAMjIymDp1KhMmTGDq1KlcuHChuSYJBAIvo6YoiD5L1fMU3J5oBnjvvfcYNmwYQ4cOxWQy8fvf/57Tp0+j1+uZO3cuU6dOdehAb775Jj4+PuzcuRNJkigoKABgwYIFTJ8+naSkJLZu3cr8+fNZvXp16z6VwKMRnkLDiHNjobanIHc8akbz+vXr6d27NwD79++ntLSUgwcPsmHDBlatWuXQQSoqKtiyZQuzZs2yfahOnTpRWFhIeno6kyZNAmDSpEmkp6dTVFTU0s8jEAi8ACEK9pjNVSWpLpi81qCn8Je//AWAgoIC23q6J0+eRJIkFixYAFi6F1pf98YbbzR4kEuXLhESEsLf/vY3Dh06hL+/P7NmzcLX15fIyEiUSiVg6dUeERFBbm4uYWFhbfMJBYJ2hBgALQhRsKfaU3Bj+Mg6yB86dIiZM2fSs2dP7rzzTubNm8fw4cMB+N///teoGFgxGo1cunSJG264gZdeeomff/6Zxx9/nLfffruNPoZAIPAmag56IqdgSTRLCqV7q4+sTJ06lWnTphEaGoq/vz/Dhg0D4NixYw630I2OjkalUtnCRDfffDOhoaH4+vqSn5+P0WhEqVRiNBq5fPky0dHRrfhIAoGgvSM8BXs8Yp6Clccee4zhw4dTUFDAbbfdZlMqhUJhCx01RVhYGMOHDyclJYURI0aQkZFBYWEhPXr0ID4+nh07dpCUlMSOHTuIj48XoSOBbBEDoAUhCvaYTEbPaHNhZeDAgXW29e/fv1kHWrhwIXPnzmXp0qWoVCqWLVtGUFAQycnJzJkzh5UrVxIUFMTSpUubtd/2hLi4BQLHECWp9nhESerq1auZNm0aGo2mwTfrdDrWrVvHgw8+2OSBunbtyqefflpne2xsLOvXr3fQ3PaN3CfgCASOUnPQk/vNlNlsrprAp3Bv+KigoIDbb7+d0aNHM3ToUHr27GnrfXThwgV++OEH9u3bR1JSktOMEwjkhtwHQCsifFSN9fNbPAU3Jpqff/55/vCHP7B582Y2bNjAmTNnKCsrIygoiLi4OEaPHs1zzz1HaGio04wTCATyRISPqrF+fo9INIeFhTFjxgxmzJjhNAMEAoGgNiJ8VE21KLgmpyAW2RG4HLn/yAVNI8JH1dhEwUWJZiEKLkTuF7egacQ1YkGIQjVGoxGo8hQUHtD7SCAQCFyNmNFcjV1OQXgK3oXc73gEAkcRnkI1VgFQKFxTkuqQKBQVFVFRUQFYXJmNGzeyefNm2Su4QCBwDsJTqMYjcwqPPfYYmZmZAKxYsYJ//vOffPLJJyxZssRphgkEckTud8VWhKdQjW2egnU5TklyvyhcuHCB+Ph4ALZt28Y//vEP/vWvf/H11187zTBvRO4Xt0DgKEaTCVwwUas9ULMkFUAhKdzf+0ihUKDX68nIyCAwMJDOnTtjMplsISWBoDkIcRQ0hclkQqFUYTLohSjYwkcWkZSc7Ck4JAqjRo1i1qxZlJSUMHHiRADOnj1LZGSk0wwTCATyxWQyoVAoMSFEwVqSSpWnICk8wFN4/fXX2bx5MyqVytbrqLi4mKefftpphnkj4g5ZIHAMq6dg/VvO2KqPrKLgCeEjjUbD1KlTMZlMFBQUEBERYVt9TSBoLkIcBU1hEQWl7W85U7P6yPqv2xPNpaWl/PnPf2bAgAGMHz8egO+++44VK1Y4zTCBQI4IwbRgCR8JTwHqJpolSVEdUnICDonCggULCAgIYNeuXajVagAGDRrEN9984zTDvBHxgxcIHEN4CtXU9hQUnpBT+N///sf+/ftRq9W2ft5hYWEUFhY6zTCBQCBfzCYTSiEKQP2egtvDR4GBgRQXF9tty8nJITw83ClGCQQCeSPCR9XUKUlVeEBJ6r333sszzzzDs88+i8lk4ujRo7z11ltMmzbN4QMlJiai0Wjw8fEBYPbs2YwcOZLU1FTmz5+PVqulS5cuvPnmm3Ts2LFln8bDEeEjC+I8NIw4NxZMJrMIH1Vh/fxFuZlcKy32jOqjP/7xj/j4+LBo0SIMBgNz585l6tSpPPTQQ8062DvvvEPfvn1tj00mEy+88AJvvPEGCQkJrFy5kuXLl/PGG28071MIBAKvQpSkVmNNKhfknEepUnuGKEiSxEMPPdRsEWiKtLQ0fHx8SEhIAGDatGmMGzfOa0VB3AUKBI5hMlsmr4EQBVv4CEv4CE9INH///fd06dKFrl27cuXKFZYvX45CoeD5559vVl5h9uzZmM1mhgwZwvPPP09ubi6dO3e2PR8WFobJZKKkpISQkJDmfxpBu0CIo6ApzDU8BblfLzYBqNIEj0g0L1y40FYJsGTJEgwGA5Ik8corrzh8oLVr17Jt2zY2btyI2Wxm0aJFLbNYIPBi5D4AguUc1CxJdWZNfnugWgCsiWbnzlNwyFPIz8+nc+fOGAwGDhw4YJuvMHLkSIcPFB0dDVhmR0+fPp0nnniCBx98kJycHNtrioqKUCgUwksQCGSMdRCsrCgFhChUl6RS9a8HeAoBAQEUFBTw448/Ehsbi7+/PwAGg8Ghg1y7do2ysjLAchfw9ddfEx8fT//+/amsrOTw4cMArFu3jjvuuKMln0PQjhB3w4LGsIpAxdUiQOQU6vMU3J5TeOCBB7jnnnvQ6/XMnTsXgCNHjtCrVy+HDlJYWMjTTz+N0WjEZDIRGxvLggULUCgULFu2jAULFtiVpHorYjAUNIW4RmouVC/ZPZYrtQVAkiTLehNOwiFRmDlzJrfffjtKpZJu3boBEBkZyWuvvebQQbp27cqWLVvqfW7w4MFs377dQXMF3oAY+ASNIUTBnupEs+V8KBRKTHo3iwJAz549gWoDu3fv7hyLBF6PEAVBY9QWAbmLgk0kqx5LkoTJ3YnmEydOsGjRIk6fPo1WqwUsP2xJkjh58qTTjBMI5IYQzJqLyljWI5a7KNQNHykwmZx3nTgkCnPmzGHs2LEsXrwYX19fpxnj7YgfvAVxHgSNUXOyloR3i0LynmQW7l3o2Ivz7R8+tfBJh962YPQCksckO2yTQ6KQnZ3Nc889Z4vxCQStQYiCoDGqPQXLf14tCmOSmxywd+/ezcGDBwnqZJno28E/mPLCSzw7a5ZTbHKoJPX222/nwIEDTjFAIBBUIwSzZgzd8j9RkmqyraUAlpJUs7urj7RaLU899RRDhgyhU6dOds8tW7bMKYYJvBcx8Akaw95T8O7wkSOYTCbbWgpQtfKau0Whd+/e9O7d22lGCOSFEIWGEeemZrhIlKSCdW0Je0/B7ZPXnnrqKacZICfED17QFOIasZ+nIAlPAaPRWMdTMBk9YJ5CSkoKX331FUVFRbz//vscP36c8vJybr31VqcZJ/BOxMDXMN5+bppVbWNd7fc03LfwPoeP0dxqG0+ndk5BoVBgMrtZFD799FNWr17Nvffey86dOwHw9fXl9ddfF6IgaDbePvC1Bm8/N45U25w8eZKNGzcSEBrB9fISukRH8eCDD7rGQA+kvpyCMyevOVR99K9//YuPP/6YmTNn2mJbvXr1IiMjw2mGCQRyxNtFwRFqVx/JPXxUX/UROO9acUgUKioqbK2vrXMVDAYDarXaKUYJvBsx8DWMODeg1+stf1TNaLY9lil1PQXnJuAdEoWhQ4fywQcf2G1bvXo1w4cPd4pRAu9GDHwNI85N3USzoy36vZW6noJzlyl1KKfw8ssv8/jjj7N+/XoqKiqYMGEC/v7+rFq1yilGeSviBy9oCnGN1PIUkNDLXBTqVh9ZPAW3ikJERAQbN27k+PHjZGdnEx0dzYABA+xqZwUCRxEDn6AxrJ6BZe6ahFHmouBqT8HhUd1kMtkU3Gg0ih+2oMWIa6dhxLnBtga8Jafg+AqP3oqrcwoOeQqnTp3iySefRKfTERkZSV5eHj4+Prz33nv069fPKYYJvBcx8DWMODcWEVAoLUOTJMJH9c5TsG53Bg6Jwty5c7n//vt5+OGHkSQJs9nMJ598wty5c9m0aZNTDPNGxA/egjgPDSPOjb0oIFka4tVu9SAnLDkFyXZtSE4WBYfO8oULF3jooYdsboskSTz44INcuHCh2Qf829/+RlxcHGfOnAEgNTWVyZMnM2HCBB555BEKCwub2INA4L3IvSMoWERBqbTEzWuWwMsVi6egtD22hpLcKgqjR49m165ddtt2797NmDFjmnWwEydOkJqaSpcuXQDLh3rhhReYP38+O3fuJCEhgeXLlzdrn4L2h7gbbhghCpbqo5qeAshbFIwmk91aNs72FBwKHxmNRp577jn69+9PVFQUeXl5pKWlMW7cOF588UXb6xpro63T6Vi0aBH/7//9P9uU9bS0NHx8fEhISABg2rRpjBs3jjfeeKM1n8ljEYOhBXEeGkacG8t4o6i6M5YQomAymVColRixFPo421NwSBT69u1L3759bY979+7NiBEjmnWgt99+m8mTJxMTE2PblpubS+fOnW2Pw8LCMJlMlJSUEBIS0qz9C9oPYuBrGOEpWDwFyZpoFp6CRRQ8zVNobevso0ePkpaWxuzZs1u1H4F3IEShYYQo1E00W7fJFZPJhLJWQzzrdmfgUE5h5cqVdX7I169fZ/78+Q4d5Mcff+TcuXOMGzeOxMRE8vLymDFjBpmZmeTk5NheV1RUhEKhEF6ClyNEwZ6a50Ocm1olqVWiIOf+RyZjrclrntD7aP/+/dx3331cunQJgCNHjjB58mTKy8sdOsjMmTM5cOAAu3btYteuXURFRfHRRx/x6KOPUllZyeHDhwFYt24dd9xxRws/iucjfvAWxHmwp+Ydn/AUQG8w2HIKiJxC3clrntD7aO3ataxatYp77rmHUaNGceDAAebNm8ekSZNadXCFQsGyZctYsGABWq2WLl268Oabb7ZqnwLPR4iCPUIU7NHr9Gh81Rj0WuEpUE+bC0/ofaRQKBg/fjxbt25l586djB07lnHjxrX4oDXLWwcPHsz27dtbvC+BoL0jwkf26PQ6Oqiq2vJXDYA6nc6NFrkXk7n+klS3rqewZs0apk+fzrRp09i7dy+SJJGUlERqaqpTjPJWxA/egjgP9tS84xPnxuIpKKtEwRo2kbWnYKy78hq4uffRhg0bWLNmDX369AHgr3/9K1u2bOHxxx/n+++/d4ph3oj4wVsQ58EeET6qxmw2o9fraoiC8BQ8svfR+vXr66yyNmXKFLHITjMRg6EFcR7sqfnjlvvSk9aEslJpLwqy9hTM9XsKbi1JVavVnDt3jvfee4+FCxcCcO7cOa5eveoUowTejdzvhmtTUwjkfm6sHoGyZk5BkmTtKZjrJJo9QBS++eYbHnjgAfLz89m6dSsA165dY8mSJU4xylsRd8gW5D7w1UaEj6qxDv4KVXVkQqVSy1YUzGZzVUlqzUSzB1QfvfPOO3z88cf069ePb775BoB+/fpx6tQppxjlrQhRsCD3ga82InxUjTVMpKwhCkqVWrbho+p22R7WJbWoqIi4uLgqg6rbZ9dUL0HTCFGwIETBHuEpVFMnfFT1t1w9Bev14MouqQ6Jwo033mgLG1n56quvGDBggFOM8laEKFiQ+8BXGyEK1dQnCn87uN4AACAASURBVAqlfD0Fq+foykSzQ+GjefPmMWPGDDZs2MC1a9eYMWMGGRkZ/POf/3SKUQLvpnZdvtw9TpForqah8JHsPQUXJpodEoXY2Fi++eYb28I60dHRjBkzBn9/f6cY5a0IT8GCEAV7RE6hmno9BSEKtXofecDkNYAOHTowceJEpxghF4QoWKgdLpHr2rtWhKdQTYM5hcoKd5nkVqzXg6IeT8GtbS4EbYMQBQsihm6POB/V2MJHyhqioFSj08kzp2C7Hmommp3cEE+IggsRomBBDIL22O4GlSoRPqpnnoJSpUanl2f4yFaSWk/4SIiCFyAGQAs1xVGck+rwkVKtwSBzUdBqtSiVKrtwiUqtQafVyvKmqv5Es/AUvAY5XtT1ITwFe6z9flQqNUaDEAWVWmO3TanWYDKZZOlF1TdPwfJY4frqo9GjRztUFbJnz562tMerEQOgBSEK9tg8BZUGo1G+K4xBlaeg9rHbZhWJyspKAgIC3GGW26gvfASWVhfOuslsUBTECmhth/XLE56CBSEK9tiFj3SOLXHrrVhEobanYBEJOZalepSnMGzYMKccUI44OwbY3hCiYE+1p6BGd11+IZKa1Bc+qukpyI365ilYH7t18hrAyZMnOXz4MMXFxXZ3vLNmzXLo/X/605/IyspCoVDg5+fHK6+8Qnx8PBkZGcyZM4eSkhJCQkJYunQpPXr0aPYHaQ8IT8GCEAV7rKKgUmu4JsO4eU0qtVqUPsF221QqiyhotVp3mORWqhPNtTwFd4SPavLFF1/wxhtvcNttt7Fv3z5GjRpFSkpKs9ZpXrp0KYGBgQD897//Ze7cuWzevJkFCxYwffp0kpKS2Lp1K/Pnz2f16tUt+zQeiggf2WNNrIJ3z+BN3pPMwr0LHX9DhuWfFxbOdvgtC0YvIHlMcvMM82C0lVoCAuxzCtbwkaxFwdM8hQ8//JAPP/yQhIQEhg4dynvvvcfevXv5+uuvHT6QVRAAysvLkSSJwsJC0tPT+fjjjwGYNGkSr776KkVFRYSFhTXzo3g+4q7YglzaOiSPSXZowN6/fz979+4lpu9gcs79zNy//MX5xnkoOl3dnII1fCRvUaidU5DcKwqFhYUkJCQAlunWJpOJ0aNH88ILLzTrYPPmzSMlJQWz2cyHH35Ibm4ukZGRKJWWXuFKpZKIiAhyc3O9UhSEp2ChpqdQ82+5Up1TUGEyGmXbD8pkMqHT6RrMKXilKBxLhrSGvcnuwMu9gEuvgPW++vAjjIgGSoHPHDhG/wUwINlhkxwShaioKLKysoiJiaFHjx589913hIaG1lm3uSlef/11ALZs2cKyZcsczke0d0Si2Z6a3oE3ewqOYjQaUSiUKKoWUjGZTLYbJTlhrS5SqWuHj7xYFAYkNzpgnz9/ns8++4xBifdy/vhBAAaNvYdDX39CbI+uTJkypc1Ncmjy2qOPPsq5c+cAS8L4hRde4KGHHuLJJ59s0UGnTJnCoUOHiIqKIj8/3zYwGI1GLl++THR0dIv266lY/QMhChaMRqOtl4sQhSpRUCptq2vJ1XuyDvq1w0cKhRKFUuWdotAE9fU+sjx0c/jorrvusv09evRofvjhB/R6vcOtsysqKigtLbUN9rt27SI4OJiOHTsSHx/Pjh07SEpKYseOHcTHx3tf6KhKFYQoWDAajahUGgx6rRAFLE3glEoVSpXl52gwGPDx8WniXd6HddCvHT6ybpOzKNSXaHZr9dGBAwfo0qULPXv2BECj0ZCdnU1OTg633XZbk++/fv06s2bN4vr16ygUCoKDg3n//feRJInk5GTmzJnDypUrCQoKYunSpa37RB6MGAAtGAwGlGq1EIUq9Ho9CqUKhVJle+y1NBJDj6Aqfp7/CuTDiJox9M7AVZwSQ/dkGpzR7G5PYdGiRaxZs8Zum5+fH4sWLWLnzp1Nvr9Tp058+eWX9T4XGxvL+vXrHTGjHWP5YoWnYMFoNKKsqj0XomARSYVKbWsX7dWi0EgM/ezZs6xbt45Bib8juFM0R3dvACwx9J/+8zkRYYHcd999rrPVA2h4noLzSlIdyikUFhYSERFhty0iIoIrV644xShvw1j15YkB0IKhhijINX5eE5unoJKBp9AIjYWPlGofET6yw3megkOi0LVrV/73v//ZbTt06BAxMTFOMcrbsIqB8BQsGI1GVFWVa+KcVIuCUg7ho0awtrFQaermU1QaH67LsM1Fdfiobu8jt+YUnnrqKZ5++mnuueceunbtyqVLl9i0aROLFy92ilHehsloGfi8fQBs9gxe4JUdr8AOx1/vbTN4AXR6PUqlBkVV+Eiu3lO1p1CPKKh9KC+RnyhUjxm15q1Ibm5z8atf/Yp//vOfbNiwgb179xIVFcWHH37IgAEDnGKUt2EyVZfcejOOzuD9f2+9RVB4d/IzT5KYmMj//d//Od84D0av16Pw9bdVH8nZU5AUClvCvSYqtQ+VlfILHzXsKbg50QwwYMAAIQItwGw22748b/cUHMVoMKKsCh95u1A6gkFvwNdfJtVHjVBZWYla7VPvbG6VxgeDQW8pUpDRxL6G2ly4xVP4+9//zhNPPAHA22+/3eAO5DIruaXU/OLEAGjBaDLaKm3EOQG9QY+fXEpSG6GysrLefAJUh5QqKysdnh/lDXhUSWpeXl69fwuaR81BT3gKlovcaDBUzVJVClEADNbJazLPKVRWVtabTwBQa6o7pcpRFOrOaFZgNjvnt9OgKCxcWJ0wfOONN5xycDkg+vzYYxVGhdLSukCud8VWzGZznZJUOa4wBhZRqL0Up5WanoKcaCingCRhcnX46NKlSw7toGvXrm1mjDciOoLaY+tvo9KgUmlkOwBa0ev1mM1mVGofFAolkkIh23NSWVmJyr9jvc9Zw0pyE4XGWmebXR0+uv322y0HrtXGt/bjkydPOsUwb0F4CvZYBzylSo1SpZa9p1CzNl+SJNQaX65fv+5mq9xDpVZLcIjwFGpSHT6qL6fgYk/h1KlTtr83btzIwYMHefrpp+ncuTM5OTm899573HrrrU4xypuoKQTCU6ghCmoNCpValrNUa1ItCr6Wf9U+shv4rGgbySnIVRQa9hQUmF0tCjV5++23+fbbb/H1tVy4PXr0YNGiRUyYMMGug6qcaMlELTLgoYUPOfxyr5yoVctTkGuoxIrVK7CGR1QaeYqCwWCwzHRvqPpIpuGjxnMKbpynYDKZyM7OJjY21rYtJydH1tU0jk7Uys3N5aOPPgIgMjKSP/7xj062zLOxFwUNWl25my1yL9ZBTl11J6xS+8gyfGTzmBrwFBRKFZKkEKJQheSORHNN/vCHP/DQQw9x1113ERUVRV5eHps2beKhhxy/65Ur1pCRUqXGIHIKNVbX0qBSa6iskLenUCd8pPGlsrTUnSa5hcb6HgFV+Rb5NcVrcJEdJPeGjx599FH69u3Lv//9b9LT0wkPD2fx4sWMGjXKKUZ5E7b1d9UajAYhCtXVR9bwkbx+5LWpPRiqND6UVgpPoT7kmG9paPKa28NHAKNGjRIi0AKsoqBSaUSimerZukqVBqVKg05UHwHVg6Fa7YtWq61T5eftNNYMz4pSrZGdKNQM0Wuvl2PU68g+ewyc6Ck41Dpbp9OxYsUKxo0bx5AhQwDLamy1F94R1MUWPlJrMBqFKNScp6BUqTEaDLLOTVln8VoFQKXxwWw2yy5MUi2OdddSsCLHNRWsNwc5545TWX4VvfY6vxzZTWV5idM8BYdEYfHixZw5c4bly5fbLt4+ffrw+eefO8Uob8LmKag1Yp4CVTkFSUKhVKJUqau3yZTr16/bxdHlWnpZ22Mym81or5dzrbSI7LPHqib4aWS3poLZbAZJojDnvN127fVy9+YU/vvf//Ltt9/i5+eHQmHRkcjISPLz8x06SHFxMS+++CIXL15Eo9HQvXt3Fi1aRFhYGKmpqcyfPx+tVkuXLl1488036dix/lmN7RFruETt42ebvSqnsEBtdDodKpUGSZJQVt0VarVaW7mz3LB2BrUi19JLW/io6vNb74wBfjmy2/Kc2ofyq/L0FGpHGcxmE3XWWGgjHPIU1Gp1nbvcoqIiQkJCHDqIJEk8+uij7Ny5k+3bt9O1a1eWL1+OyWTihRdeYP78+ezcuZOEhASWL1/e/E/hwVjvgjW+/pZmcDL3FnQ6nc1DEJ4ClJeXo/b1sz3WVP1dXi6vUl2tVoskSbZOsbXvjAtzzqOSYfWRJdFc/+DvrNbZDonCHXfcwUsvvWTrh3T58mUWLVrEb37zG4cOEhISwvDhw22PBw4cSE5ODmlpafj4+JCQkADAtGnT+Pe//93cz+DRVIuCn91juaLT6WxLcaqq1mmW8zkpLinB1z/I9tjXPxiAkpISd5nkFrRarV1upfadsdFoQKXWoNfrZJeDqj+yIOEcSXBQFJ577jliYmKYPHkypaWlTJgwgYiICJ588slmH9BkMvH555+TmJhIbm4unTt3tj0XFhaGyWTyqh+ENYau9vGtfixjdDodiioxsIaP5HpOKisr0VZW4usfTN6Fk+RdOInG1w+FUulVvwFHsCTcG04yA7YOqnLyFhoLN7t1OU6NRsPcuXOZO3cuRUVFhIaGtjgu/uqrr+Ln58cDDzzAf/7znxbtoz2h1+tRVc3eBfkOgFZE+Kiaq1ctMXNf/yCyz/4MQFSPeHz9gmQnClqttsG22VZUNXJQHTp0cIVZbqexgd8topCTk1Pv9tzcXNvfNe/0m2Lp0qVkZmby/vvvo1AoiI6OtjtGUVERCoXC4VxFe8A6CMp9ALRSVlaOOrATAGofyw/ba+Pnx5IhreH+WJHAy72A/FcgsGrj4UcYEQ7ogM8cOEb/BTAguZWGup/GFtixIsfKLGv1UR0kwElRtEZFITEx0eYR1KdKkiQ53Dr7rbfeIi0tjQ8++ACNxqL4/fv3p7KyksOHD5OQkMC6deu44447mvsZPBqLKGiEKGApz716tYSu0b0B8OkQgEKpoqioyM2WOYkByY0O2IcOHeI///kPtyXNJO3gVwAMGnsPZ37aRWH2L7wwe7Zr7PQAKrValOqARl9T01OQC2azGcnFieZGRaFfv35UVlZy5513MnnyZCIiIlp0kF9++YVVq1bRo0cPpk2bBkBMTAzvvfcey5YtY8GCBXYlqd5EbU9BzusHlJSUYDab8Qu0eIKSJOEXEExhYaGbLXMPJSUllh5QGvtyXF//YLSVlVRWVsqmVFdbqaVDaOOl6HLNKdSvCZJ7RGHLli2cOXOGzZs3c9999xEbG0tSUhLjx49v1sXap08fTp8+Xe9zgwcPZvv27c2zuh1hSayK8BFg8wg6BITatvkGhHivp9AExcXF+PoH1cnPWauRSkpKiIqKcodpLker1RLYRKLZ6inIKXwENOgp4KT6oyarj/r27ctLL73Erl27+MMf/sCePXsYMWIEJ06ccIpB3ka1pyASzTZRCKzOGfkFhlJSUiK7MkOAkpKr+PgF1dleUxTkgNlsRqfTOpxTkJ+nUFcUJNw8TwHgwoUL/Pjjj6SmphIfH09QUN2LWVAXnU6HUik8BbCIgkrjg7pGuKRDQIjXlSE7gtlspqTWHAUrHWQmCjqdDrPZbCtRbgj55hQafs4ZNBo+Kikp4auvvmLz5s1UVFSQlJTEmjVrmlVxJHfKyysID+uKUqVGoVR5b6UNNFlt82vg1zHATy/aP9EL+Pcrjh3DS6ptrl27hsGgtwlATVQaX1Rqja1k1dtpai0FKwqlCoVSKavwUYPVRzgreNSEKIwcOZKYmBiSkpK4+eabAcjMzCQzM9P2GrFOc8NUVlai02nx8QtEkiR8/QIp9eYFVJqotnn7nXfoEBLNDbdUV5jpKis4uO1DJkyYwNChQ51vo4dw5coVwD6UZkWSJDoEhHD58mVXm+UWbAsvqRr3FKyvkZO33eDkNUkCd3gK4eHhaLVavvzyS7788st67JL47rvvnGKYN2AVAF8/SxG6j1+AbO7+amMwGCgrLaVj13i77WofP1RqjewqkDIzM0GSCOpYv9cd3KkzWRlpGAwGVCqHlz1pl1gr8hQOfE65revtjt5HjX4Lu3btcspB5YJVFHz8LPXXPh0CKS246E6T3EZxcTEAHQJD7bZb74rlVoGUmZlJYEg46gZCJiERMWT9kkpOTg7dunVzsXWupXrhJXWTr1Wo1PJbrKqBpILbE82C5mP1CnxsnkIg5eXlsuyUavUEOgTUDZf4BoRQKCNRMBgMZGdnExwe0+BrgsO7ANiFar0VmygoHRAFpUp2nkLDJanOEQYhCk6ktLQUSZLw8fUHwLfKYygrK3OnWW7B6gn41SMKfoEhlF69KhuxzMrKwmg0EhrRsCioNb4EhIRz4cIF1xnmJqyDvEPhI6VaVhNAG0s0OwvvDlY6kyYqbQDGAmN7Akdert7YC/haXpU2YBEFja+frcIk74KlPUpUj3g6BIRgNpspLi6mU6dO7jTTJVjzCcGdGq/iC4mIIev8ca/PKzTLU1Cp0OllVn3UxPNtvWiX915pzqaJShuA1atXc/WagUGJ9wJQUVrEj//+lKSkJG666Sbn2+hBFBUV0SEg2PY4N8My+TGqR7ytAqeoqEgWonChKp/QVAlmSHgMWWeOkp2dTffu3V1knetpTk5BqVKjr/DiCr56aCx85AxE+MiJXL1aaksyQ3UVkleXpTZAYWGRXXuLmvhVbZdDslmv15OdnU1II6EjK8HhFk/C2/MKzRIFpRqdjMJH7kCIgpMwm82UlZXaksxguejVGl/ZiYJWq6WiorzemnwAtY8vao0vBQUFLrbM9WRnZ2MyGglpJMlsxZpX8HZRsOUUlE0HLhQqlaxyCu5AiIKTKC8vx2Qy2bwDKz5+gbJpX2Dl558tC8g0NhAGR8Rw8uQpr29hYMsnhDvWFSAkIoas7GyvLsM0GAwolCqHYuNKpRqD3EShgS6pzkKIgpM4e/YsAIGhkXbbA8MiuXjxomzK6oxGIwcP/o+Q8C4Ed4pu8HXd+g1Bq63kp59+cqF1rufChUwCQyOabP5mJSQiBmNVCau3otPpUDkQOgKLt200GmXZQNFVCFFwEsfT0ugQEExgmL0oRHaLQ6/Xc+bMGTdZ5lqOHTtGeXkZ3eIbb2ERFBZFaGQ3vj90yGvDA5WVlWTnZDsUOrIS3KkLSBLnzp1zomXuRa/XOxQ6guoQk7deI7VxpPqorRGi4ARKS0u5mJlJRLd+dVzi4PAu+PoFcvz4cTdZ5zpMJhMpKQcJDI0gNLLpWbnd44dyraLCFm7yNtLS0jAZjUR06+vwe9QaH8KiunPs2HGvvTvW6/UOJZkBmS5WVX/rbGchRMEJWNeaiOweV+c5SZII79qX8+fPU1FR4WrTXEp6ejolJcV0u2GYQ/Hi4KoQ08GDB71yItvRo0cJCAknMLR5KxhG97yR8vIyr/UWmuMpiBb0zkeIghM4fvw4QWGR+AXWX4IZ2b0fZrPZ4fWt2yNms5mUlBT8gzvSqXMvh94jSRLd4odSWlrqdZ5Ubm4u+fn5RPe6sdnv7di5JxpfP44ePeoEy9yPdXVCR1Cq5BU+cgcuEYWlS5eSmJhIXFycXSw9IyODqVOnMmHCBKZOneoVU/qvXLnC5cuXiejWr8HXBIR0IiC4o9cNfDU5c+YMV65coVu/hGbNuAyL6kFASDgHDx70qnBJamoqCqWy0euiIRQKJZHd4zl79qxXrsdhNBpRKJQOvVaqep03epKegktEYdy4caxdu5YuXbrYbV+wYAHTp09n586dTJ8+nfnz57vCHKdy/PhxJEkiolufRl8X0a0f2dnZtu6h3oTZbOZASgodAoIJ7+p4/ByqvYWioiJOnTrlJAtdi16vJy0tjfCYPg12RW2K6J43YjKZOHbsWBtb535MJpPDNw6SpLC9R+AcXCIKCQkJREfblyMWFhaSnp7OpEmTAJg0aRLp6entelar2WwmLe0EIZHd0FQ1wWsIa7IxLS3NFaa5lIyMDHJzcugaNwSFovmXWHiXWPwCQzlw4IDT2gO7kpMnT6LVaonu2fzQkRW/oFCCwzuTmprqFeekJiaz2TbYN4VVPLztHHgSbssp5ObmEhkZiVJpcQeVSiURERHk5ua6y6RWc+nSJUpLrxLZzT7BnHfhpK0BnBVf/yCCw7twPC3N6y7wlJQUfDr4E9UjvukX14OkUNAtPoHLly/b5nu0Z1JTU+kQEGxrh91SonveSFFREZcuXWojyzwDs8nkcCdQqygIT8GCKEn1cNLS0lAoVXTqEmu3PTfjhK0BXE0iu8VRVFhIXl6eq0x0OllZWWRmZhLTd7DDFSX1EdEtDl//wHbvLRQWFnLx4kWiet7Y6m6W4TF9UKk1XpdwNpma0emzyqNoz9dEc6n31DixJtVtohAdHU1+fr4tYWQ0Grl8+XKdMFN74fr166Snn6RTl16o1E2vNQsQ3rUPCoXSq37kBw4cQO3jS+fY1nWBVSiUdI0bQnZ2drvu/fPzzz8jSRJRPW5o9b6UKjUR3eI4efKkVy1ebzKbkBwMMwpPwfm4TRQ6duxIfHw8O3bsAGDHjh3Ex8cTFhbmLpNazLVr11izZg1anZYuvW92+H1qjS9RPW/gyJEj7N+/34kWuoa8vDzOnj1Llz6DHJ6M1BhRPW5E4+tHSkpKG1jneoxGIz///DNh0T3x6dB4jslRonreiMFgsM2F8QbMzfAUJBl6Cq7GJaLw2muvMWrUKPLy8nj44Yf5zW9+A0BycjJr1qxhwoQJrFmzhoULG1+0xhMpLy/n00/XcOVKAf1v+22TC6fUps+gMUR2j2fv3r3s3r27XV/sKSkpqNSaZgljYyhVKmL6DiYjI4OcnJw22acrOXfuHBUVFS2am9AQgaERBIR08irv0mQ2NTvRLDwF5+GSRXZefvllXn755TrbY2NjWb9+vStMcAplZWV8umYNJSVXuWnkZIdaOdRGUijoN+x2FAoFKSkpGI1Gxo0b1+arKTmbwsJCTp48Sbd+CS0uu6yPzrE3cenUjxw4cIDf/e53bbZfV3D8+HE0Ph0Ii+rR6OvMZjPa6+UY9Tqyzx6jc+xNDX7/kiQR2T2ecz/vp7CwkI4dOzrBctdiNjm+5KQcq4/q/ahO/Pgi0dxCrl69yr/+tZqrV0sZMGpKiwTBiiRJ9E0YR5feN/P999+zc+fOdnfRp6SkoFCqiOk7qE33q1Jr6Nx7IGfOnOHy5cttum9notPp+OWXX+gU06fJstycc8epLL+KXnudX47sJudc45MareXM3hJCas48BWQ6T8Gg16JWqxk2bBhqtRqT0dJK3Rk3j0IUWkBxcTGrV6+mvOIaA0bfSUgrSw3B8uX2HjSamL6DOXz4MF999VW7EYaioiKOHz9O51790fj6tfn+Y/oMRKXWsHfv3jbft7M4ffo0BoPBoeZ3hTnnG31cG58OAYSEx5CWdqLdXCONYW5B+MgbPrfjmDHodQwaNIjx48czcOBAjCbnzegWotBMCgsLWb16NdeuV3Lz6DsJ7th21VKSJBF78wi6xQ8lNTWVbdu2tYs7on379qFQKOkWn+CU/at9fOnSZxCnT59uN+W7J06cwNcv0KEck9FoaPRxfUR060tRUSH5+fktttFTMNWz+HztO2OD3rL4krVKqT38LtoCSZLAjK0UeefOnZaWKVXtPoSn4GauXLnC6tWr0eoM3Dzm7jprJbQFkiTR66b/o0f/Wzl+/Dhbtmzx6D4vV65cIS0tjc69b25yFndr6Np3EGqND3v27HHaMdqK69evc/78ecK79nFabig8pjeSpPCKEJLZVNdTqH1nbNBbuqLKLdEsSRJmzKjUPuj1en788ceqrrJCFNxOfn4+q1evRm80c/OYuwkICXfq8XrcMIxeA0aQnp7Opk2bPFYY9u3bh1Klpmu/IQ6/x5pYvVZaRPbZYw6FAlQaH2L6Dubs2bMevwrZyZMnMZlMRHSr2zq9rVD7dCA0qhsnTrT/EJKpnkRz7Ttj69wfWYaPXPxZhSg4QG5uLqs//RSTpGTg2HvwD3ZNxUe3fkPoPWg0p0+fZv369R63Tm9+fj4nT54kps8gND4dHH5fcxOrVrr0GYjapwN7PDy3cOLECfwCQ5x+4xDRtS+lpaVkZWU59TjOxpJTqC0K9nfGtuVLZZZobsoTEJ6CG8jKymLNmjVISg0Dx9zT4BoJziKmz0D6Dknk7NmzfPHllx7VR37v3r2o1D7ExDWv4qi5iVUrKrWGrnFDyDh/nosXLzbrmK6irKyMzMxMwrvGOb2suFOXWBRKZbsPIUmSwvE7f7NFDFrSaLG9YnZm/Wk9yOfMtoBLly6xdu1nKNQdGDjmbjoEBLvFjs6xN9Fv2O1kZGTw+efrPEIYcnJyOHPmDDF9B6HW+DbrvS1JrFrp0nsAGl8/9uzZ45EhBOvCSc1ZcrOlqNQaOkb3JD09vV3fOSuVCktTPAewfk5rI01vx5podiVCFBqguLiYL778ErWvHwPH3o2vf5Bb7YnqcQPxwydw8WIm27dvd/uAuGfPHtQaX2L6DnTpcZUqNd36DeXixYseuShTWloaASHh+Ae5pl1LRLc4rl275pHnwlGUSiVmB0ssra+Tiyg0hrM8USEK9aDVavniiy8wGs30H5GET4eAFu+rJUnVhojsFmdLPruzV9KlS5c4f/48XfsNqY71upDo2P74+AV4nLdQXFxMTk6OS7wEK2FRPVCpNe06hKRQKh32dGTpKYjwkXsxmUxs2rSJwsIibvi/ifgFhrRqfy1NqjZE17jBRPaIZ9++faSnp7dqXy1lz969aHz92qzHUXNRKlV0jx9Gdna2Ry1mX1hYCEBwp9ZPZnQUpUpFQGgkBQUFLjtmW9MST0FWOYV621w0o914M5HPmXWQ7777jnPnztFn8BhCI7q2en8tTao2hCRJxA1JJLhTqhP4fwAAFuhJREFUNNu2bXP5okQXLlwg88IFuvZLaJNOqC0lqscN+PoHsWfvXo/xFqp/pK62x9yuB0mlUonJQVEQnoLzab9XkhM4evQohw4dokvvm1u9HoCV1iRVG0KhVHHj/01C5dOBL774grKyslbv0xHMZjN79uzFp4N/m52flqJQKul+wzDycnM5c+aMW22xYh2YXS1SZnP7FwVHE80ip1AD4Sk4l8zMTL755htCI7sRO3CUu81pEo2vH/1v+y3XK7UuK1U9f/48WVmX6BY/FGUrVlVrKyK7x+MXGOIx3oJNFJpZCdRQSwdHMZtM7VsUFAqHPQWzHD2FBq5tET5yIsXFxazfsAFf/2BuuHViu/mBBYSEE3/LHeTl5jq9IslsNrNn7158/QJbtQB9W6JQKOh+w3CuXL5sKwV1tz1gmYzVHBpq6eAw5nYuCs3wFEwy8xQsbS4aeM5Jx2y/V1IbUV1pZKL/iN+26VoArqBT5170GnCb0yuSzp49S25ODt1uGNaqtZfbmoiuffEPCmPv3n1ur9VvqSg01NLBUcxOTDq6AqVSicksPIUGqedmz7JFeAptjl2l0a2/cfls5baia9wQIrs7tyLp0KEf8OkQQFSP+Fbvq7XhkppICgXd4odSWFhARkZGq21rDdXho+Z5bA22dHAQsww9hfb8eZuDtSFeQ885A3mc2QawVhr1HjSa0MjWVxq5C0mSiEtIJKij8yqSwsJC0VVeo7KitNX7anW4pBalhXlIkkRISOvKh1tLSz2F1tLucwr1lKTWzllZH8vNU7DmFGqfD0mShCi0NampqbZKoy69B7jbnFajUKrof9skVBpLRVJ5eXmb7n/UqFGoVErOHz/Y6n21NlxSk2tlxeScO87AgQPdvjSl20TBC6uPOnbuVe9jueUUFAoFZnPd8+HrH4Sk8GJRyMjIYOrUqUyYMIGpU6e6ZMq+NTEZEhHj1OO0ZaikKVRqH/xDwikvL2/zzpkBAQHceuutFGSd5WpBTqv21dpwSU3OH0tBpVIyevToVtnUFrS0+qi1tPfwkaKe6qPOsTfhGxCC2qcDfQaPtZVAy9FTMJtNdc6HX2CY0zwFj8gYLliwgOnTp5OUlMTWrVuZP38+q1evduoxJ0+ezLovvuDE/76mz8DRdOnjnNm5NUMlZrOZY+mnnXMcnZYTB7+i+PIlxowZQ1xc2/fyv+WWW/jpyBHO/XyAQYn3uj25ebUgh4Lsc4waNYqAgJa3Imkr3NXrv12Fj44lQ9pCu013AkQAh5+x2z4youqPEuCnGk/0Aja/0vAx+i+AAcmtNNRF1HM+anI7cHs34Kc59ucDIBr4bHbTx2jm+ZDMbi7wLiwsZMKECRw6dAilUonRaGT48OF8++23hIU13lQsMzOT8ePHs3btWqKiopp9bL1ez86dO7lw4QLRvfrTrV+CwwNd57z36Hzl702+7rI+ht3aGQwYOIRjqT8x1ucjItSO3cXnhD9BTtSTTb5Oe72C0z9+S2VFKYmJifTr18+h/beEEydOsHv3bvoMGUvHqB52zzl6TlqD9ZyYzWZOHPwao66cB3//e9Rq982utlJZWcmHH35Il9430zVusEvOR1bHx9j2cxA33XQTI0eOdOqxnMW+ffs4cfIUQ8ffb7f9xP++BuDGWyfatl06c4TsX37mySefdPtNiSs4ePAgqampDPv1Q3bn4/zxFMquXGLGjBmt2n9UVBQqVa18hbtFIS0tjZdeeomvvvrKtm3ixIm8+eab3Hhj4/Xwhw8f5v7772/0NZ5Ap06diIqKIi8vr133qBEIBN7Fd999R0yMfQjdI8JHLaV///6sXbuW8PBw2cQYBQKBoK2oL8LidlGIjo4mPz8fo9FoCx9dvnyZ6OjoJt/r6+tLQkKCC6wUCAQCeeD27FTHjh2Jj49nx44dAOzYsYP4+Pgm8wkCgUAgaHvcnlMAOHfuHHPmzKG0tJSgoCCWLl1Kr169mn6jQCAQCNoUjxAFgUAgEHgGbg8fCQQCgcBzEKIgEAgEAhtCFAQCgUBgQ4iCQCAQCGy4fZ6Ct7N06VJ27txJdnY227dvp2/fvu42ya0UFxfz4osvcvHiRTQaDd27d2fRokWyL0H+05/+RFZWFgqFAj8/P1555RXi41u/dkV7529/+xvvvvuu+O0AiYmJaDQafHwsDSRnz57tlNYmQhSczLhx43jwwQfbRTsOVyBJEo8++ijDhw8HLKK5fPlyFi9e7GbL3MvSpUsJDAwE4L///S9z585l8+bNbrbKvZw4cYLU1FS6dOniblM8hnfeecfp4ijCR04mISHBodnZciEkJMQmCAADBw4kJ6d1rbi9AasgAJSXl8ui2Vtj6HQ6Fi1aRHJysrtNkR3CUxC4DZPJxOeff05iYqK7TfEI5s2bR0pKCmazmQ8//NDd5riVt99+m8mTJ9dp1iZ3Zs+ejdlsZsiQITz//PMEBQW1+TGEpyBwG6+++ip+fn488MAD7jbFI3j99dfZs2cPzz33HMuWLXO3OW7j6NGjpKWlMX36dHeb4lGsXbuWbdu2sXHjRsxmM4sWLXLKcYQoCNzC0qVLyczM5K9//Wv7WSDGRUyZMoVDhw5RXFzsblPcwo8//si5c+cYN24ciYmJ5OXlMWPGDA4cOOBu09yKNQyt0WiYPn06R44cccpxRPhI4HLeeust0tLS+OCDD9BoWr4+s7dQUVFBaWmp7Ue/a9cugoODCQkJcbNl7mHmzJnMnDnT9jgxMZH3339f1tVH165dw2g0EhgYiNls5uuvv3ZadZoQBSfz2muv8e2331JQUMDDDz9MSEiI3YJCcuOXX35h1apV9OjRg2nTpgEQExPDe++952bL3Mf169eZNWsW169fR6FQEBwczPvvvy/7ZLOgmsLCQp5++mmMRiMmk4nY2FgWLFjglGOJhngCgUAgsCGCuQKBQCCwIURBIBAIBDaEKAgEAoHAhhAFgUAgENgQoiAQCAQCG0IUBO2Wbdu28cgjj7TZ/nQ6HRMnTuTy5csAzJkzhxUrVjT4+kGDBnHp0iWH9h0XF0dmZmab2FmbnJwcBg0ahNFodMr+BRaWLFnCZ5995m4znI4QBQ9g+/bt3HXXXQwaNIgRI0bw6KOPcvjwYXeb5fFMnjyZf/7zn7bHrR14v/jiCxISEoiIiHDo9UePHqVr164tPl5b0blzZ44ePYpSqXS3KV7NI488wqpVq9DpdO42xakIUXAzH3/8MYsXL+bxxx8nJSWF3bt3M336dL777jt3m2bDYDC424Q6OMOmdevWkZSU1Ob7bUs88buoj/ZipxVH7I2IiKBXr17s2rXLBRa5DyEKbqSsrIx33nmH+fPnM378ePz8/FCr1SQmJvLSSy8BcOzYMaZOnUpCQgIjRoxg0aJFdncqcXFxfP7554wfP56EhAQWLlxIzfmIX375Jb/+9a8ZNGgQEydO5MSJEwDk5+fz9NNPc8stt5CYmMjq1att73n33Xd55plnmD17NoMHD663r//evXuZOHEigwYNYuTIkXz00Ue253bv3k1SUhIJCQlMmzaNU6dOAfDBBx/wzDPP2O3ntdde47XXXrOdj7lz5zJixAhGjhzJihUrbCGRTZs2MW3aNBYvXszw4cN599132bRpE/fddx+Abb2KpKQkBg0axNdff82kSZPsfsB6vZ7hw4eTnp5e5/Pk5ORw6dIlbr75ZrvtpaWlzJw5k0GDBnHvvfdy8eJFu3Nv9UyKi4t5/PHHGTx4MHfffTcrVqyw2Wbl4MGDDX5PGzZs4Ne//jVDhw5lxowZZGdn2x1n7dq1jB8/nvHjx9exPSsri7i4ONvAtmnTJsaNG8egQYNITExk27Ztdd4D9t/zoEGD+O1vf0tGRgarVq3i1ltvZfTo0Xb9hvLz83n8/7d3/jFNXW0c/wKFIcWhwhAUXRY3GUw72kHLBjWCLWIHBOIQNHFTpoxNVJRtsLluRjOCis4wpgvJlJgZwaGRhbBVSLT8YcyazPFDNEC3ohQEGb/bQm/b5/3DcGOhSPV9ja/kfv7qvX3Oc5/znNP7nPuc23OysiAWiyGXy3H+/HmHuib6jM1mQ2lpKWQyGSQSCXbv3o3BwUGHtvT39+Ojjz5CeHg4xGIxNm3aBJvNxl7XUV/t6emBQCCw09nS0gKJRAKGYZ7Ir1qtFlu3boVYLMbatWtRU1NjZ6dYLIZarXZYh1kDcTwz1Go1hYSEEMMw08o0NTXRjRs3iGEYunv3LsXHx9Pp06fZ75cvX06ZmZk0NDREer2eJBIJqdVqIiKqqamh6OhoamhoIJvNRjqdjjo7O8lqtVJKSgp9//33ND4+Tnfu3KHY2Fiqr68nIqLi4mIKDQ2l2tpaslqtZDKZptgVFRVFGo2GiIgGBwepubmZiIhu3rxJkZGR9Ndff5HFYqGLFy9STEwMjY+PU2dnJwkEAhoZGSEiIovFQlFRUXTjxg0iIvrkk09IqVSSwWCgvr4+Wr9+PZ07d46IiC5cuEAhISF05swZYhiGTCYTXbhwgdLT0+18odPp2OPS0lLavXs3e1xbW0sJCQkO/XzlyhVSKBR25/Ly8kgsFlNDQwMxDEN79+6lnJwch9fLycmhnJwcMhqN1NbWRqtWrZpi23TtVFtbSzKZjNrb24lhGPrhhx8oLS3NruyWLVtoYGDAYVvcvXuXli9fTgzDkMFgIKFQSFqtloiIenp6qLW11WGdi4uLacWKFVRfX08Mw9Bnn31GMTExdOLECTKbzVRRUUExMTGs/KZNm+ibb76hsbExamlpIYlEQteuXWN1Te4zZWVllJqaSt3d3TQ+Pk5KpZL27Nnj0JaioiJSKpVkNpvJbDaTRqMhm802Y1/dvHkzVVRUsHoKCwtJqVQ+kV8NBgOtWrWKKisriWEYunnzJonFYmpra2PLqFQqSk5OdliH2QL3pPAMGRwcxPz588HjTb8E1YoVKxAWFgYej4egoCCkpaVBo9HYyWzfvh0vvvgiFi1aBIlEwo7MKysrsW3bNggEAri4uODll1/G4sWL0dTUhP7+fmRnZ8PDwwNLlizBhg0b7EZFYWFhkMlkcHV1haen5xS7eDwe2tvbMTo6Ch8fH7zxxhsAHuTl09LS8Oabb8LNzQ0pKSlwd3dnd9AKDQ1FXV0dAOD69evw9PREWFgY+vr6oFar8eWXX8LLywu+vr7YsmWL3TpR/v7+2Lx5M3g8nkObJpOUlAS1Wo3R0VEADyamk5KSHMoODw+Dz+dPOS+TySAQCMDj8ZCUlIRbt25NkbFarbh8+TJ27tyJOXPm4NVXX0VycvIUuenaqby8HJmZmVi2bBl4PB6ysrJw69Ytu1FtZmYm5s2b51S9XV1d0dbWhrGxMfj7++O1116bVjY8PBxSqRQ8Hg/x8fEYGBhAZmYm3N3doVAooNfrMTw8jO7ubvz555/49NNP8cILLyAkJASpqamoqqpidU3uM+Xl5dizZw8CAgLg4eGB7OxsqFQqh6kaHo+H+/fvo6urC+7u7ggPD4eLi8uMfTUxMRHV1dUAwC4Ul5iY+ER+vXr1KhYvXoz169eDx+MhNDQUa9euxe+//87K8/l8DA8Pz9gGzzPcgnjPkHnz5mFgYAAWi2XawPDPP/+gsLAQzc3NMJlMsFqt7A14gpdeeon9PGfOHBgMBgBAd3c3li5dOkWnXq9Hb28vwsPD2XNWq9XuOCAg4JG2FxcX4+TJkzh69CiCg4ORm5sLoVCIrq4uXLp0CT///DMryzAM+0ZPQkICqqurkZycjOrqaiQkJAB4kL6xWCyIjo5my9lsNrtd62ayaTILFy6ESCSCSqWCXC5HfX099u3b51DWx8eH9dvD+Pn5sZ89PT1hNBqnyPT398NisdjZ6mi3venaqaurCwUFBTh06BD7PRGhp6eH3YrS2d37vLy88N133+HUqVPYt28fRCIR8vLysGzZMofyvr6+dvWbP38+O2E9EYCMRiN6e3vh4+MDb29vVn7RokVobm5mjye3T1dXF3bs2GG3NLqrqyv+/fdfLFy40E72ww8/RElJCfs2WVpaGjIzM2fsq3FxcTh48CB6e3uh0+ng6urKfve4ftXr9WhsbJxyrYcHEgaD4alsbPP/BBcUniFCoRAeHh6oq6tDfHy8Q5n9+/cjNDQUR48ehbe3N8rKyqBSqZzSHxgYaJcDf/h8UFAQLl++PG3ZmVboFAgEOHnyJBiGwdmzZ5GTkwO1Wo3AwEBkZWXh448/dlhu3bp1OHToEO7du4fa2lpUVFQAADuavH79+rQB8klWDU1JScEvv/wCq9WKsLCwKTejCYKDg9HZ2fnIAD0dCxYsAI/Hw7179/DKK68AeBCQnWXCZ9M9xQCPV3epVAqpVIqxsTEcP34cSqXyv36V0t/fH0NDQxgdHWUDQ3d3t50/J9sYEBCAgoICvPXWWzPq9/b2Rn5+PvLz89Ha2ooPPvgAK1eunLGv+vj4ICoqCjU1Nfj777+hUChYOx7Xr4GBgYiIiMDp06enlddqtXj99ddnrM/zDJc+eobMnTsXu3btwoEDB1BXVweTyQSGYaBWq9mdtwwGA/h8Pvh8PrRaLc6dO+e0/vfeew+nTp1Cc3MziAgdHR3Q6/UQCATg8/koLS3F2NgYrFYrWltb0djY6JRes9mMX3/9FSMjI3B3dwefz2dHg6mpqSgvL0dDQwOICEajEVevXmVTOAsWLIBYLMYXX3yBoKAgdgTr7++PqKgoFBYWYnR0FDabDXfu3MEff/zhdH39/Pym/G9AJpOhpaUFZ86ccZjSmSAgIABLly512gcP4+bmBrlcjpKSEphMJmi1Wru0ykykp6ejtLQUbW1tAB5MuP/222+PbQcA9PX1oa6uDkajER4eHvDy8vqfbGIUGBgIoVCIY8eOYXx8HLdv30ZlZeUjb7gbN27E8ePH2XRNf38/mzqczJUrV9DR0QEiwty5c+Hm5gYXFxen+mpiYiKqqqqgUqnY1BHw+H5dvXo1dDodLl26BIZhwDAMGhsbodVqWRmNRgOpVOqc055TuKDwjMnIyEB+fj5OnDiBt99+G6tXr8bZs2chk8kAAHl5eaiuroZIJIJSqYRCoXBa97p165CVlYXc3FyIRCLs2LEDQ0NDcHNzw48//ojbt29jzZo1iIyMxFdffcXeuJ2hqqoKsbGxEIlEKC8vx5EjRwAAK1euxMGDB3HgwAFEREQgLi4OFy9etCubkJCAa9eusamjCQ4fPgyGYaBQKBAREYFdu3bh/v37TtuUnZ2N/Px8hIeHszlnT09PxMXFobOzE3K5/JHl09PTH+tm/jBff/01RkZGEBUVhc8//xzvvvuu0xsIyeVybNu2DXv37oVIJEJCQgLq6+ufyA6bzYaysjJIpVKIxWJoNBrs37//iXRN5tixY9Dr9ZBKpcjOzsbOnTvxzjvvTCv//vvvIzY2FhkZGRAKhdiwYcO0QbejowNbt26FUChEWloaNm7ciMjISKf6amxsLHQ6Hfz8/OxG8Y/rV29vb/z000+oqamBVCpFdHQ0ioqK2Lf9ent70d7ezv42Zyvcfgocs56SkhLodDoUFRU9Us5sNiM5ORllZWVO/4FtOo4cOYK+vj67fDbH801hYSGWLFnCvv48W+GCAsesZnBwECkpKTh8+DAiIiKe2nW0Wi0YhkFwcDCampqwfft2fPvtt7N+VMkx++AmmjlmLefPn0dBQQGSkpKeakAAHsz95Obmore3F76+vsjIyMCaNWue6jU5OJ4G3JMCBwcHBwcLN9HMwcHBwcHCBQUODg4ODhYuKHBwcHBwsHBBgYODg4ODhQsKHBwcHBwsXFDg4ODg4GD5DyN7ItsUuumxAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot results stratifying across cancer severity (one of many covariates)\n",
    "cancer_vecs = val_X[:,40:45]\n",
    "for col in range(cancer_vecs.shape[1]):\n",
    "    one_val = cancer_vecs[:,col].max()\n",
    "    cancer_vecs[:,col] = (cancer_vecs[:,col] == one_val).astype(int)\n",
    "cancer_dx = cancer_vecs.dot(np.arange(5)+1).astype(int)\n",
    "\n",
    "# Trim missing labels\n",
    "val_labels_trimmed = val_labels[cancer_dx > 0]\n",
    "val_lower_trimmed = val_lower[cancer_dx > 0]\n",
    "val_upper_trimmed = val_upper[cancer_dx > 0]\n",
    "cancer_dx_trimmed = cancer_dx[cancer_dx > 0]\n",
    "valid_cancer_dx = np.arange(5)+1\n",
    "\n",
    "lower_uncal = [(np.e**(val_lower_trimmed[cancer_dx_trimmed == dx]) - 1).mean() for dx in valid_cancer_dx]\n",
    "lower_cal = [(np.e**(val_lower_trimmed[cancer_dx_trimmed == dx]-qhat) - 1).mean() for dx in valid_cancer_dx]\n",
    "upper_uncal = [(np.e**(val_upper_trimmed[cancer_dx_trimmed == dx]) - 1).mean() for dx in valid_cancer_dx]\n",
    "upper_cal = [(np.e**(val_upper_trimmed[cancer_dx_trimmed == dx]+qhat) - 1).mean() for dx in valid_cancer_dx]\n",
    "\n",
    "sns.set(style='white')\n",
    "sns.set_palette('pastel')\n",
    "plt.figure()\n",
    "plt_idx = np.random.choice(np.arange(val_labels_trimmed.shape[0]), size=(val_labels_trimmed.shape[0],), replace=False)\n",
    "sns.violinplot(x=cancer_dx_trimmed[plt_idx], y=np.e**val_labels_trimmed[plt_idx]-1, color='#BBD9FA')\n",
    "plt.plot(valid_cancer_dx-1, lower_uncal, marker='_', markersize='25', linewidth=0, color='orange', label='uncalibrated')\n",
    "plt.plot(valid_cancer_dx-1, upper_uncal, marker='_', markersize='25', linewidth=0, color='orange')\n",
    "plt.plot(valid_cancer_dx-1, lower_cal, marker='_', markersize='25', linewidth=0, color='green', label='calibrated')\n",
    "plt.plot(valid_cancer_dx-1, upper_cal, marker='_', markersize='25', linewidth=0, color='green')\n",
    "plt.legend(loc='upper right')\n",
    "plt.ylim(0,80)\n",
    "plt.xlabel('Cancer severity (higher is more severe)')\n",
    "plt.ylabel('Medical expenses ($)')\n",
    "sns.despine(top=True,right=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This method was first introduced in https://arxiv.org/abs/1905.03222."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}