zero_inflated_negative_binomial_regression_multi_gp.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "d6ef6310-9746-4bc9-b262-40e1a96493c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.special import expit  # Sigmoid function for logistic link\n",
    "\n",
    "from typing import Any, Optional\n",
    "\n",
    "import torch\n",
    "from torch import Tensor\n",
    "from torch.distributions import Laplace\n",
    "\n",
    "from gpytorch.constraints import Interval, Positive\n",
    "from gpytorch.distributions import base_distributions\n",
    "from gpytorch.priors import Prior\n",
    "from gpytorch.likelihoods.likelihood import _OneDimensionalLikelihood\n",
    "\n",
    "from typing import Any, Optional\n",
    "\n",
    "import torch\n",
    "from torch import Tensor\n",
    "from torch.distributions import Laplace\n",
    "\n",
    "from gpytorch.constraints import Interval, Positive\n",
    "from gpytorch.distributions import base_distributions\n",
    "from gpytorch.priors import Prior\n",
    "from gpytorch.likelihoods.likelihood import _OneDimensionalLikelihood\n",
    "\n",
    "\n",
    "import urllib.request\n",
    "import os\n",
    "from scipy.io import loadmat\n",
    "from math import floor\n",
    "import math\n",
    "import torch\n",
    "import gpytorch\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "from gpytorch.models import ApproximateGP\n",
    "from gpytorch.variational import CholeskyVariationalDistribution, DeltaVariationalDistribution, MeanFieldVariationalDistribution\n",
    "from gpytorch.variational import VariationalStrategy\n",
    "from sklearn.preprocessing import MinMaxScaler, StandardScaler\n",
    "\n",
    "\n",
    "N = 1000\n",
    "r = 1\n",
    "pi=0.35\n",
    "X = torch.normal(mean=0, std =1, size=(1,N)).double().cuda().flatten()\n",
    "mu = torch.exp(X)\n",
    "probs = r / (mu + r)\n",
    "y = base_distributions.ZeroInflatedNegativeBinomial(\n",
    "    total_count = r, probs=probs, gate=torch.tensor([pi]).double().cuda()\n",
    ").sample().double().cuda().flatten()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "dbb5cd07-c229-49c6-b941-745f86a095af",
   "metadata": {},
   "outputs": [],
   "source": [
    "from typing import Any, Optional\n",
    "\n",
    "import torch\n",
    "from torch import Tensor\n",
    "from torch.distributions import Laplace\n",
    "\n",
    "from gpytorch.constraints import Interval, Positive\n",
    "from gpytorch.distributions import base_distributions\n",
    "from gpytorch.priors import Prior\n",
    "from gpytorch.likelihoods.likelihood import _OneDimensionalLikelihood\n",
    "\n",
    "class ZeroInflatedNegativeBinomialLikelihood(_OneDimensionalLikelihood):\n",
    "    def __init__(\n",
    "        self,\n",
    "        batch_shape: torch.Size = torch.Size([])\n",
    "    ) -> None:\n",
    "        super().__init__()        \n",
    "        \n",
    "    def forward(self, function_samples: Tensor, *args: Any, **kwargs: Any) -> base_distributions.ZeroInflatedNegativeBinomial:\n",
    "        gate_logits = function_samples[..., 0]\n",
    "        total_count_log = function_samples[..., 1]\n",
    "        probs_logits = function_samples[..., 2]\n",
    "\n",
    "        r = torch.exp(total_count_log)\n",
    "        n = torch.sigmoid(probs_logits)\n",
    "        pi = torch.sigmoid(gate_logits)\n",
    "\n",
    "        return base_distributions.ZeroInflatedNegativeBinomial(total_count = r, probs=n, gate=pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c0a3edc3-5ddc-40b5-a4d8-ad128df3dac9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "50\n",
      "1.1518954381578825\n",
      "100\n",
      "1.1496353102893404\n",
      "150\n",
      "1.1488953317649913\n",
      "200\n",
      "1.1487279331499645\n",
      "250\n",
      "1.1486820693474034\n",
      "300\n",
      "1.1486571388788025\n",
      "350\n",
      "1.1486358772188006\n",
      "400\n",
      "1.1486158004266862\n",
      "450\n",
      "1.1485966260710303\n",
      "500\n",
      "1.1485783361754751\n"
     ]
    }
   ],
   "source": [
    "class GPModel(ApproximateGP):\n",
    "    def __init__(self, inducing_points, num_tasks):\n",
    "        variational_distribution = gpytorch.variational.NaturalVariationalDistribution(inducing_points.size(0), batch_shape=torch.Size([num_tasks]))\n",
    "        variational_strategy = gpytorch.variational.IndependentMultitaskVariationalStrategy(\n",
    "            gpytorch.variational.VariationalStrategy(\n",
    "                self, inducing_points, variational_distribution, learn_inducing_locations=True\n",
    "            ), \n",
    "            num_tasks=num_tasks\n",
    "        )\n",
    "        \n",
    "\n",
    "        super(GPModel, self).__init__(variational_strategy)\n",
    "        self.mean_module = gpytorch.means.ConstantMean(batch_shape=torch.Size([num_tasks]))\n",
    "        self.covar_module = gpytorch.kernels.ScaleKernel(\n",
    "            gpytorch.kernels.LinearKernel(batch_shape=torch.Size([num_tasks])),\n",
    "            batch_shape=torch.Size([num_tasks])\n",
    "        )\n",
    "\n",
    "    def forward(self, x):\n",
    "        mean_x = self.mean_module(x)\n",
    "        covar_x = self.covar_module(x)\n",
    "        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)\n",
    "\n",
    "inducing_points = X.clone()\n",
    "M = 100\n",
    "inducing_points = torch.normal(mean=0, std =1, size=(1,N)).double().cuda().flatten()\n",
    "model = GPModel(inducing_points=inducing_points, num_tasks=3).double()\n",
    "likelihood = ZeroInflatedNegativeBinomialLikelihood().double()\n",
    "\n",
    "if torch.cuda.is_available():\n",
    "    model = model.cuda()\n",
    "    likelihood = likelihood.cuda()\n",
    "\n",
    "model.train()\n",
    "likelihood.train()\n",
    "likelihood.initialize()\n",
    "\n",
    "variational_ngd_optimizer = gpytorch.optim.NGD(model.variational_parameters(), num_data=y.size(0), lr=0.1)\n",
    "\n",
    "hyperparameter_optimizer = torch.optim.Adam([\n",
    "    {'params': model.hyperparameters()},\n",
    "    {'params': likelihood.parameters()},\n",
    "], lr=0.01)\n",
    "\n",
    "mll = gpytorch.mlls.VariationalELBO(likelihood, model, num_data=y.size(0))\n",
    "\n",
    "n_epochs = 500\n",
    "for i in range(1, n_epochs + 1):\n",
    "    if i % (n_epochs // 10) == 0:\n",
    "        print(i)\n",
    "        print(loss.item())\n",
    "    variational_ngd_optimizer.zero_grad()\n",
    "    hyperparameter_optimizer.zero_grad()\n",
    "    output = model(X)\n",
    "    loss = -mll(output, y)\n",
    "    loss.backward()\n",
    "    variational_ngd_optimizer.step()\n",
    "    hyperparameter_optimizer.step()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "18751fd8-1f60-4a16-9c27-e1c17f9a0966",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Axes: xlabel='X1', ylabel='y'>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjIAAAGwCAYAAACzXI8XAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAOmNJREFUeJzt3X90VIWd///XnfwOyUx+EAKRBMIvgVoUkB8Ri4KsrO26/kA/LtuvUuX0VDdiEd21bD+t+/G0xdN2l1bqr+66qFsBpf2gtf22VqlCiyAIpqJWBAkkEBMSIBNIIAnJ/fwBGTPJ/M5M7tyZ5+OcnAP3zr3zvvPzdWbue96GaZqmAAAAbMhhdQEAAACRIsgAAADbIsgAAADbIsgAAADbIsgAAADbIsgAAADbIsgAAADbSrW6gFjr7u5WXV2dcnNzZRiG1eUAAIAQmKapU6dOqaSkRA6H/89dEj7I1NXVqbS01OoyAABABGprazVy5Ei/6xM+yOTm5ko6f0M4nU6LqwEAAKFoaWlRaWmp533cn4QPMj1fJzmdToIMAAA2E+y0EE72BQAAtkWQAQAAtkWQAQAAtkWQAQAAtkWQAQAAtkWQAQAAtkWQAQAAtkWQAQAAtkWQAQAAtkWQAQAAtpXwIwoSxcHG0zp8ok2jC4eofOgQq8sBACAuEGTiXHNbh+5bX6Wt+xs9y+aOL9KaxVPlyk6zsDIAAKzHV0tx7r71Vdp2oMlr2bYDTVq2/j2LKgIAIH4QZOLYwcbT2rq/UV2m6bW8yzS1dX+jqptaLaoMAID4QJCJY4dPtAVcf+g4QQYAkNwIMnFsVEF2wPWjCznpFwCQ3AgycWxMUY7mji9SimF4LU8xDM0dX0T3EgAg6VkaZP7t3/5NhmF4/U2cONGz/uzZs6qsrFRhYaFycnK0aNEiNTQ0WFjx4FuzeKrmjBvqtWzOuKFas3iqRRUBABA/LG+//sIXvqA33njD8//U1M9Luv/++/Xb3/5WGzdulMvl0r333qubb75Z27Zts6JUS7iy0/T80pmqbmrVoeOt/I4MAAC9WB5kUlNTNXz48H7L3W63nnnmGa1bt07z58+XJK1du1aTJk3Sjh07NHv27MEu1VLlQwkwAAD0Zfk5Mvv371dJSYnGjBmjr371q6qpqZEk7d69W52dnVqwYIHnshMnTlRZWZm2b9/ud3/t7e1qaWnx+gMAAInJ0iAza9YsPfvss/r973+vJ598UtXV1frSl76kU6dOqb6+Xunp6crLy/Papri4WPX19X73uWrVKrlcLs9faWlpjI8CAABYxdKvlq677jrPv6dMmaJZs2Zp1KhReumll5SVlRXRPleuXKkVK1Z4/t/S0kKYAQAgQVn+1VJveXl5mjBhgg4cOKDhw4ero6NDzc3NXpdpaGjweU5Nj4yMDDmdTq8/AACQmOIqyJw+fVqffvqpRowYoenTpystLU2bN2/2rN+3b59qampUUVFhYZUAACBeWPrV0oMPPqjrr79eo0aNUl1dnR5++GGlpKRo8eLFcrlcWrp0qVasWKGCggI5nU4tW7ZMFRUVSdexBAAAfLM0yBw5ckSLFy/W8ePHVVRUpCuvvFI7duxQUVGRJGn16tVyOBxatGiR2tvbtXDhQj3xxBNWlgwAAOKIYZp9RisnmJaWFrlcLrndbs6XAQDAJkJ9/46rc2QAAADCQZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2RZABAAC2FTdB5tFHH5VhGFq+fLln2dmzZ1VZWanCwkLl5ORo0aJFamhosK5IAAAQV+IiyOzatUtPP/20pkyZ4rX8/vvv16uvvqqNGzdqy5Ytqqur080332xRlQAAIN5YHmROnz6tr371q/rP//xP5efne5a73W4988wz+o//+A/Nnz9f06dP19q1a/X2229rx44dfvfX3t6ulpYWrz8AAJCYLA8ylZWV+spXvqIFCxZ4Ld+9e7c6Ozu9lk+cOFFlZWXavn273/2tWrVKLpfL81daWhqz2gEAgLUsDTIbNmzQnj17tGrVqn7r6uvrlZ6erry8PK/lxcXFqq+v97vPlStXyu12e/5qa2ujXTYAAIgTqVZdcW1trb75zW/q9ddfV2ZmZtT2m5GRoYyMjKjtDwAAxC/LPpHZvXu3jh07pmnTpik1NVWpqanasmWLHnvsMaWmpqq4uFgdHR1qbm722q6hoUHDhw+3pmgAABBXLPtE5pprrtHevXu9lt15552aOHGiHnroIZWWliotLU2bN2/WokWLJEn79u1TTU2NKioqrCgZAADEGcuCTG5uri655BKvZUOGDFFhYaFn+dKlS7VixQoVFBTI6XRq2bJlqqio0OzZs60oGQAAxBnLgkwoVq9eLYfDoUWLFqm9vV0LFy7UE088YXVZAAAgThimaZpWFxFLLS0tcrlccrvdcjqdVpcDAABCEOr7t+W/IwMAABApggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALCtVKsLAEJxsPG0Dp9o0+jCISofOsTqcgAAccLST2SefPJJTZkyRU6nU06nUxUVFfrd737nWX/27FlVVlaqsLBQOTk5WrRokRoaGiysGIOtua1DdzyzU/P/fYvuXLtL8378lu54ZqfcbZ1WlwYAiAOWBpmRI0fq0Ucf1e7du/Xuu+9q/vz5uuGGG/Thhx9Kku6//369+uqr2rhxo7Zs2aK6ujrdfPPNVpaMQXbf+iptO9DktWzbgSYtW/+eRRUBAOKJYZqmaXURvRUUFOhHP/qRbrnlFhUVFWndunW65ZZbJEkff/yxJk2apO3bt2v27Nk+t29vb1d7e7vn/y0tLSotLZXb7ZbT6RyUY0B0HGw8rfn/vsXv+jcfvJqvmQAgQbW0tMjlcgV9/46bk327urq0YcMGtba2qqKiQrt371ZnZ6cWLFjguczEiRNVVlam7du3+93PqlWr5HK5PH+lpaWDUT5i4PCJtoDrDx1vHaRKAADxyvIgs3fvXuXk5CgjI0N33323Nm3apMmTJ6u+vl7p6enKy8vzunxxcbHq6+v97m/lypVyu92ev9ra2hgfAWJlVEF2wPWjC/k0BgCSneVdSxdffLGqqqrkdrv1y1/+UkuWLNGWLf6/TggmIyNDGRkZUawQVhlTlKO544u07UCTunp9A5piGJozbihfKwEArA8y6enpGjdunCRp+vTp2rVrl37605/qtttuU0dHh5qbm70+lWloaNDw4cMtqhaD6WDjad12+Uid6TinXYdPepbPGTdUaxZPtbAyAEC8sDzI9NXd3a329nZNnz5daWlp2rx5sxYtWiRJ2rdvn2pqalRRUWFxlYil5rYO3be+Slv3N3qWzRidryVXjNYXSlx8EgMA8LA0yKxcuVLXXXedysrKdOrUKa1bt05vvfWWXnvtNblcLi1dulQrVqxQQUGBnE6nli1bpoqKCr8dS0gMvlqu9xxuVlbaET2/tMSiqgAA8cjSIHPs2DHdcccd+uyzz+RyuTRlyhS99tpr+pu/+RtJ0urVq+VwOLRo0SK1t7dr4cKFeuKJJ6wsGTF2sPG01ycxPbpMU1v3N6q6qZVPZAAAHnH3OzLRFmofOuLDm/uO6c61u/yuX3vnDM27eNggVgQAsILtfkcGkGi5BgCEhyCDuNLTcp1iGF7LUwxDc8cX8bUSAMALQQZxZ83iqZozbqjXMlquAQC+xF37NeDKTtPzS2equqlVh463anThED6JAQD4RJBB3CofSoABAATGV0sAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2wg4yS5Ys0datW2NRCwAAQFjCDjJut1sLFizQ+PHj9YMf/EBHjx6NRV1IMgcbT+vNfcdU3dRqdSkAABsxTNM0w92osbFR//M//6PnnntOH330kRYsWKClS5fqhhtuUFpaWizqjFhLS4tcLpfcbrecTqfV5aCP5rYO3be+Slv3N3qWzR1fpDWLp8qVHV+PJQDA4An1/Tuic2SKioq0YsUK/eUvf9E777yjcePG6fbbb1dJSYnuv/9+7d+/P+LCkVzuW1+lbQeavJZtO9CkZevfs6giAICdDOhk388++0yvv/66Xn/9daWkpOjLX/6y9u7dq8mTJ2v16tXRqhEJ6mDjaW3d36iuPh8Kdpmmtu5v5GsmAEBQYQeZzs5O/epXv9Lf/d3fadSoUdq4caOWL1+uuro6Pffcc3rjjTf00ksv6ZFHHolFvUggh0+0BVx/6DhBBgAQWGq4G4wYMULd3d1avHixdu7cqcsuu6zfZebNm6e8vLwolIdENqogO+D60YVDBqkSAIBdhR1kVq9erVtvvVWZmZl+L5OXl6fq6uoBFYbEN6YoR3PHF2nbgSavr5dSDENzxg1V+VCCDAAgsLC/Wrr99tsDhhigr0Ct1WsWT9WccUO9ls0ZN1RrFk8drPIAADYW9icyQKhCaa12Zafp+aUzVd3UqkPHWzW6cAifxAAAQsaIAsRMOK3V5UOHaN7FwwgxAICwEGQQE7RWAwAGA0EGMUFrNQBgMBBkEBO0VgMABgNBBjHR01qdYhhey1MMQ3PHF3EuDAAgKggyiBlaqwEAsUb7NWKG1moAQKwRZBBz5UMJMACA2OCrJQAAYFsEGQAAYFsEGQAAYFsEGQAAYFsEGQAAYFsEGQAAYFsEGQAAYFsEGQAAYFsEGQAAYFsEGQAAYFsEGQAAYFsEGQAAYFsEGQAAYFsEGQAAYFsEGQAAYFsEGQAAYFsEGQAAYFuWBplVq1ZpxowZys3N1bBhw3TjjTdq3759Xpc5e/asKisrVVhYqJycHC1atEgNDQ0WVQwAAOKJpUFmy5Ytqqys1I4dO/T666+rs7NT1157rVpbWz2Xuf/++/Xqq69q48aN2rJli+rq6nTzzTdbWDUAAIgXhmmaptVF9GhsbNSwYcO0ZcsWzZ07V263W0VFRVq3bp1uueUWSdLHH3+sSZMmafv27Zo9e3a/fbS3t6u9vd3z/5aWFpWWlsrtdsvpdA7asQAAgMi1tLTI5XIFff+Oq3Nk3G63JKmgoECStHv3bnV2dmrBggWey0ycOFFlZWXavn27z32sWrVKLpfL81daWhr7wgEAgCXiJsh0d3dr+fLlmjNnji655BJJUn19vdLT05WXl+d12eLiYtXX1/vcz8qVK+V2uz1/tbW1sS4dAABYJNXqAnpUVlbqgw8+0J///OcB7ScjI0MZGRlRqgoAAMSzuPhE5t5779VvfvMbvfnmmxo5cqRn+fDhw9XR0aHm5mavyzc0NGj48OGDXCUAAIg3lgYZ0zR17733atOmTfrjH/+o8vJyr/XTp09XWlqaNm/e7Fm2b98+1dTUqKKiYrDLBQAAccbSr5YqKyu1bt06vfLKK8rNzfWc9+JyuZSVlSWXy6WlS5dqxYoVKigokNPp1LJly1RRUeGzYwkAACQXS9uvDcPwuXzt2rX62te+Jun8D+I98MADWr9+vdrb27Vw4UI98cQTIX+1FGr7FgAAiB+hvn/H1e/IxAJBBgAA+7Hl78gAAACEI27ar4FwHGw8rcMn2jS6cIjKhw6xuhwAgEUIMrCV5rYO3be+Slv3N3qWzR1fpDWLp8qVnWZhZQAAK/DVEmzlvvVV2nagyWvZtgNNWrb+PYsqAgBYiSAD2zjYeFpb9zeqq8/56V2mqa37G1Xd1OpnSwBAoiLIwDYOn2gLuP7QcYIMACQbggxsY1RBdsD1ows56RcAkg1BBrYxpihHc8cXKaXPDymmGIbmji+iewkAkhBBBrayZvFUzRk31GvZnHFDtWbxVIsqAgBYifZr2IorO03PL52p6qZWHTreyu/IAECSI8jAlsqHEmAAAHy1BAAAbIwgAwAAbIsgAwAAbIsgAwAAbIuTfZHUmKINAPZGkEFSYoo2ACQGvlpCUmKKNgAkBoIMkg5TtAEgcRBkkHSYog0AiYMgg6TDFG0ASBwEGSQdpmif/3rtzX3H+BoNgO3RtYSktGbxVC1b/55X11IyTNGmWwtAojFMs88ZjwmmpaVFLpdLbrdbTqfT6nIQZ5JtivYdz+zUtgNNXic6pxiG5owbqueXzrSwMgDwFur7N5/IIKkl0xTtnm6tvnp3ayXLbQEgcXCODJAk6NYCkIgIMkCSoFsLQCIiyABJgm4tAImIIANESSQtzYPdBr1m8VTNGTfUa1kydGsBSFyc7AsMUCQtzVa1Qbuy0/T80plJ160FIHHxiQwwQJEMoLR6aGX50CGad/EwQgwA2yPIAAMQyQBKhlYCQPQQZIABiKSlmTZoAIgeggwwAJG0NNMGDQDRQ5BB0gilQyicLqKDjad1+ESbZozO99nSPGN0vg4db+23L9qgASB6mLWEhBdKh1A4XUS+LpufnaaTbZ1+/993X+62zn5DKxneCACfC/X9myCDhBfKoMRwhin6u+y0sjz90/xxeuKPB7SnpjmkfdEGDQC+MTQSUGiDEs0L/w50mZ6QEWh/uw6fVIphaNfhkyHtS0quoZUAEAucI4OEFkqHUDhdRMEu+15t/xDjb18AgIHjExkktFA6hIJ9u9q7iyjY/qaW5oe8LwDAwPGJDBJaKB1C4XQRBbvs3AlFdCQBwCAiyCDhhTIoMZxhisEua9VgxsEeQAkA8YCuJSSNUDqEwukiCnbZwepIsmoAJQDEEu3XFxBkkOjCaR0HALsI9f2br5YAG2MAJYBkR5ABbIwBlACSHUEGsDEGUAJIdgQZQP47fuK9E4gBlACSHT+Ih6Tmr+Pnezd+Qf/75Q9t0Qm0ZvHUfgMoB6PdGwDiAV1LSGr+On6cWalqOXPOVp1ADKAEkEgYGgkEEWgA5Mm2Tp/LfQ1+jBcMoASQjDhHBkkrWMePP3QCAUD8IMggaQXr+PGHTiAAiB8EGSStQB0/+dlpdAIBgA1YGmS2bt2q66+/XiUlJTIMQy+//LLXetM09d3vflcjRoxQVlaWFixYoP3791tTLBKSvwGPv6680pLBj+GK9/ZwAIg1S0/2bW1t1aWXXqq77rpLN998c7/1P/zhD/XYY4/pueeeU3l5ub7zne9o4cKF+uijj5SZmWlBxUg0ruw0Pb90ps+OH3/L4wGDIgHgvLhpvzYMQ5s2bdKNN94o6fynMSUlJXrggQf04IMPSpLcbreKi4v17LPP6h/+4R9C2i/t10hEDIoEkOhsPzSyurpa9fX1WrBggWeZy+XSrFmztH37dr/btbe3q6WlxesPSCQMigSAz8VtkKmvr5ckFRcXey0vLi72rPNl1apVcrlcnr/S0tKY1gkMNgZFAsDn4jbIRGrlypVyu92ev9raWqtLAqKKQZEA8Lm4DTLDhw+XJDU0NHgtb2ho8KzzJSMjQ06n0+sPSCQMigSAz8VtkCkvL9fw4cO1efNmz7KWlha98847qqiosLAyJIpotC7Hqv052H79tY0Haw+nXRtAorG0/fr06dM6cOCA5//V1dWqqqpSQUGBysrKtHz5cn3ve9/T+PHjPe3XJSUlns4mIBLRaF2OVftzqPsN1DY+mPUCgNUsbb9+6623NG/evH7LlyxZomeffVamaerhhx/Wz3/+czU3N+vKK6/UE088oQkTJoR8HbRfo69otC7Hqv3ZbvsFgFixxfTrq6++WoFylGEYeuSRR/TII48MYlVIZIEmXoc62Toa+0iE/QJAPIjbc2SAWIhG63Ks2p/ttl8AiAcEGSSVaLQux6r92W77BYB4QJBBVPR0w2z9pDGuu2Ki0bocq/bnaO+35z4xLmxPuzaARBQ3s5ZihZN9Y8tXN0yPeO2Kcbd1atn69wbUwRONfcRqv77uk4oxhTIM6e1Pj0e1XgCIlVDfvwkyGBBf3TA94r0rJhqTrWM1HXsg+w3UofR/bvhCXE7zBoC+bNG1BHvz1w3TI967YsqHDvzNPBr7iOZ+g3UoSdK8i4cNuD4AiBecI4OIBeuG6UFXzOChQwlAsiHIIGLBumF60BUzeOhQApBsCDKImL8umx50xQw+BkoCSDYEGfgVyoBBX8MLe4QyxDAe9D3OLfuO6aebP9GfApz/E8n+/bWmv7izRstffE8b360Nqb5g633dJ5NG5OrBa0Mf7RGqgQyh7LstAy0BRIKuJfQTyYDBni6bVIehc92mLbpifB1nT/098rPT9OvKK1VaGNrXaMH232Pu+CLdc9UY3f7fO72uL9Vh6NeVczT5IlfQ+yHY+r/UNuvbL+/VB0dbfK4fiIEMofS1bX52mk62dUa9TgD2Rfv1BQSZ8CXLgMFAreO95Wen6b3vXhvV/acYht/rTXUYOvCDLwe9Hwa6fiAGsu9QbvdEfLwBCE+o7998tQQvPe27fd9kerdSJwJ/x+nLybbOsL9mCrb/QNd7rtvUE28eCHg/bP2kMcj6YzG7HwfyGAn1dk+0xxuA2CHIwEuytO+G2jreY0/NyZjuv69gwem92sD1vFfbHHD9QO7HgTxGwr1dEuXxBiB2CDLwkiztu6G2jveYVpYf0/339aXxRQHXTy0NXM/U0ryA6wdyPw7kMRLu7ZIojzcAsUOQgRcr23cHs2slWOt4b/nZaT6DxefdSMf61R1Ka7o/qQ5D/zRvnM/tHZIuKXGqtCA74P00d8KwqNyP/u6TS0qc/V48Qtl3qLf7QB5vdD8ByYWTfdFPrAYi+jOQDpiB8HWcoXQthToo09f+e1+u8uqx+uoz7/jtWgq0vRR8EORA7sdQB0/6O/ZAgh1XOPsKVjPdT4B90bV0AUEmcrEaiNiX1V1SfY/zT/sbtafmpKaV5fv8JCbcQZnBWtM3vlurbZ82ac7Yobr18lKf9S1bv0cf1bWoV+YJeRBkJPdjqB1dDkOaXOLUmsXTwn6M3PrU29p9+KTXMTkMafqofG28+4qw9uWvZrqfAPtiaCQGLFYDEXsLNuRwMAZO9j3OL40v8nuOSiSDMoPdjrdeXuozwPQwTdPrt2D6XpcUeBBkuPdjsGPsrduUz9pCuY5dh/qfsNxtSrsOnQz7fo+HxxEAa3CODCxlty4pKwZlDvZtFEnHVbg1RPuY7PY4AhA9BBlYym5dUlYMyhzs2yiSjqtwa4j2MdntcQQgeggysJTdhhxaMShzsG+jcDq6Iq0h2sdkt8cRgOghyMByvoYcxnLg5EAHHf6vGSM1uSTX5/pAdYc7/LE3X7fRtFF5MbuN1iyeqqlleV7LLhuZp0nDvY/b1/GGevuGc78fbDyt9TtrtGFnjd/9hnsb0aYNJAa6lhA3Yt0lFe1Bhz2+WOLU3VeN1eSLXD7rHujwx977+fpz72rX4ZMBLzdQvupxZqaq5ew5z/8vucipH9z0RU0ZmRfycfoT6H5vbuvQPb/Yo+0HvVu+rxhbqCe/Or3ffkO5jWjTBuyBWUuwnfKhQzTv4mEx+xrgvvVV2nagyWvZtgNNWrb+vYi27fHRZ6f04rtH/NYd7HpDreu+9VXaU9McUf3h8FVP7xAjSX+tO6Ufv/ZJ0O1CqS/Q/X7f+qp+IUY6/1s2vvYbym00kMcBgPhDkEFSiOWgw0D7CHa9wYY/9uxzsIZ5RjrUMRb1BWsD77vfUGpIlqGoQDIhyCApDMagQ1/7CLZtsOGPPfscrPbiSIc6xqK+UGrpvd9QaqBNG0g8/CAeksJgDDr0tY9g2wYb/tizz8FqL450qGMs6gullt77DaWGYKcE0qYN2A+fyCApDKQ9dyAt18Gud+6EopDqGqz24kiHOsaivp59+tN3v6HUQJs2kHgIMkgaA2nz9rVtqPsIdr2B1m/Zd0w/3fyJ/rS/0eflppbl6bbLR3rO//DXThxs3fqdNVqzeb827KzRgwsn9LseZ6b3h7fTRuXp6ouLvNqhw22nDqX1+YFrJ/Rr+ZbOdy352q+vtvHeNRxsPK3bLh+paX0uM60sdq3sgdACDgwc7ddIOpG0eftq2b2k5EILcmleVK6393qHId34+DadbOv0rO+ZxH3ONPVBnVvPv33I57wi6fN2YlOm31ZjU6bP1mZJmjk6Xx1dpqpqm0M6tt7t0MHaqUNtNe97ufHFQ/SPM0fpaj8dTr62mTE6X/91xwyft0PflvLBbMGmBRwIjunXFxBkEA2DPVl56iN/8AoxPfKz0/Ted68NOp26pzZJfuuWFPJwyFDMHV8U9LYI9XaM5PYOtI3U/3boazAnZTOpGwiO6ddAlAz2ZOUt+475DDGSdLKtUy+9WxM0gPSejB3OuoEIdluEejtGcnsH2yYUgzUpm0ndQHRxjgwQxGC37FYdaQ64/u1P+38VFC8G0sY+kFbuSCZ2h7P/aKIFHIguggwQxGBPVr6s18/++3LF2MKoXl80DaSNfSCt3JFM7A5n/9HEpG4guggyGLBIOy9i1bERbL+9O4FCqcFfy67D6N8CHKierZ80hjQ0sqxwiPL9nPDpzExVUW6mZozOD9gi3VPb3PFF/Z7kDp1fd0lJdM8Zmzu+SKZpBjzGSy5yytGnbEPSqIIs7Tp0XG/uOybjQiu0r9v7ixc5deh4a78uLb/3kc6flD1jVODbS4p+C7a/x5aVLeB0SSERcbIvIhZp50WsOjaC7ffw8dZ+nUC9BarB3dape17Y3e9rnYoxhXrq/+s/vNBfPb6uy9flZo4u0L6GFrnPfN5Vk+owdK7786drfnaa32ORpBmj8mUYhnYeOtFvXYpD6ur2u2nYZo4uUIrD8OqACtQ5FcwVYwtlmvLZUeXL3PFF+v6Nl+jbL3/g93qC3V7OzFT9dtmXVFo4sE93Qnl8u9s6tWz9e4PWtUSXFOyIrqULCDKxE2nnRaw6NoLt118nUKg13PHMTv15f6O6Q9wmUGdR7+0C1f31ueXaU3NSb3zYoI8+O9XvMtPK8vRP88fpiT8e0O7DJzXQbJKdnqKznV3qDuFVwZA0bliOfn7H5Xr4lQ/D6hhySDIMqcvP9fRse6bjXEjH1fv2rG5q1bL1e/RRXYvXcaQYhmSYAQNcKJ1XwYTz+I71xPdIagLiBdOvEVORDt+L1dC+YPt9cWdNwBATrIae/fd9D/S3TaiDJoMNjRyZn63rp5Rob12Lz8vsOnxSKYahXVEIMZLU1hFaiJEkU9L+Y6dVe6It8GBMH+u65T/E9N421OPqfT+YpqkPjrb0O44uM3CIkfoPogxXuI/vWE98j6QmwG4IMohIpJ0XserYCLbf7dWhd/pE0hXTd5tQu2hCGRo50MGTsWb19fcWyu0Vyj4iFY8dSfFYExBNBBlEJNLOi1h1bATbb0V56J0+kXTF9N0m1C6aUIZGDnTwZKxZff29hXJ7hbKPSMVjR1I81gREE0EGEYm08yJWHRvB9nvbzDK/nUCh1BBu3aEOmgxlaGSkgycjkZ+dFvJ+Qj0Gf+tS+7YvhbhtoMsHu70CXacUWhdaIPE4lDIeawKiiSCDsFsyey7va7jgnHFD9cC14wPubyDDGwPxtd9po/L0v2acH6r468orA4aZnhr83R7h1v3AtRM0aUT/gYc92/XcTr5ux0kluXpw4YSA1z12WLbn2Hxd19RSl8YW+X6TSvHxfj6heIj+9bpJfmvuq/dt66u+kfmZ+seZpT7XTRqRq+98ZZJyMlJ87ntySa4euHa8zwGP/mp54NrxWr/zsNb88RNdfXGRJpd4H0dhTpruuWqshqT7vs4Jw3L04LUTfK7r3bLvT6DnxeSSXM9t1fuyPf/vu/9IB4D6rela38/VgTznaOVGvKBrKYmF25Lp7/IPXjtBx9s6VJCdpn//w/6Q9xerjo3qplZ9eNSt594+pF2HPz9/o6eW9482a0/NSU0ry9fI/GxPDfnZaSHdHsHqDtR2LUmXleYpOz3Fq5V77vgifeOqcq36/z/WB3UtPq//8PFWXb/mz16DDsM1YViOMtIc2nu0JfiF/Zg0PFeZqQ69d8Tdr869R5t19y9263R7l2ddfnaafrF0pv7t1Y/8DrmUzndBhfNilJOR4nU94Zg8wqkxQ4eotf2cqo+36tDxz88jCday3zO8s6dN29f9XTGmUGc7u/Sej6GbfdvA+7bV9/1/KANA+z6//D5XF07Q8daOAT3naOXGYKH9+gKCjH/htmQGu3w8tXhGe+hgOPUHG+joS4phyJmVqpYz5yJuH7dST517jzb7rDHVYcg0FdZtEmtzxxdJ8j9UM9Bt3jO8U4rs/g5HKANA+z4+Y/lcjKfnORIb7dcIKNyWzGCX3/rJsbhp8Yyk3TRaLarB2q796TJNnWzrHFD7uJV66vRX47luM65CjCS/reGe23yX/9v8ZFun55ehI7m/wxGojd3X4zOW7da0ciMeEWSSVLTbiX19hB5of7EUi6GDodYfzeGFvYXTPo7oCParwntqTsbs/g5X78dnLNutaeVGPCLIJKlotxNPLc0La3+xFIuhg6HWH83hhb2F0z6O6KgYE/g2n1aWH7P7O1y9H5+xbLemlRvxiCCTpKLVTvx5C+6wkPY3kE6HUDs5whny2LOdv0GF/urf+skxr1r67ifcRmhDUmaqo99AxRTD0OWj8jXMlem32yYeGJK+WOJUZprvlxSHobBvk1grdmZo/LCcfre5JBXnpqskL8tvl1tGqqE//vWYjjaf0SUlzpgeW8/j9pKL+l9PimFowrAc/XTzJ9r4bq2kwM/VGaPzPUM3+/L12O478DTU5wkwmDjZN4mFO7gu2OUDrQ+n46KvQF0S/vb7/Rsv0UP/932/Qx59bedrUGEogw+DDSPsq2JMoQxD/WrzJS8rTc1n4vfcmESXk+6Qw+EYUKdYNKQ6pHMhzGpIdRj6deUcXZSf3e+52PdxGulQz0DPE7qWEE10LV1AkAku3DboYJf3tX4gnQ6BtpX8d3JI8jvkMdB2/+eGL4RUfyRSHNKnP/iKJGnKv71m+RskgsvPTtPI/KwBtawPplSHoQM/+LKkz5+LT/zxgPbUNIf8HAok0PMEiKZQ379TB7EmxKnyoeG9EAW7fN/1PZ0OffXudPC3v2Db+hKNdfMuHha0hkh0dUsb363VsNwMQoxNnGzrjOuOsb7OdZva+G6tbr28VOVDh8i8MFy0r0DPhUD8PU8Aq3CODGJuIJ0OVnSFRDoAMlTbPm1S1ZHmqO4T6G3bp02ef8fqOUSHEuIFn8gg5gbS6WBFV0ikAyBDNWfsUA3LzYjqPoHe5oz9fBxBrJ5DdCghXvCJDGJuIEPrgnZLxWBduAMgw5HqMHTr5aW66uJhQYdYIj7kZ6d5fgXYDnoeYz0ieQ4FQocS4o0tPpF5/PHH9aMf/Uj19fW69NJLtWbNGs2cae1PYR9sPK3DJ9oG9US3QNcZbj3R3Fco19MzALH3d/Jjiobo6ouHqrqpVaZp6vCJNh1rOav6lrMa4cxUl3m+ZdfftvMnFulYy1mdbG3X3l7ziSaX5HqG4fXt3Jg8Ilcj8zN1tPmMxg4bok8aTnvWTRuVp5H5GfrHn2/XpBFOXXXxMFXVntShpjaV5mdqTFG29h/r/3F6ODN/rho/VH/a36iL8rI0c3S+Xv/omPw1o+RlparjXLfaOkNoV0HMzC7Pl8PhkDMjRS0RznYKl8OQuoOce5tqSOd8XMaQqdk/eEMlziw1tp5VTkaa/mFmqdo6zundXufKZKUbuqzUqfwh6fpL7Qm5z4Z2bFPL8lSan6GbfrZNRc50XToyTwcaT6v1bJdys1LVcqZTuZlpumnaRao7eUYvV9XpbOc5lQ/N0U3TLtJFeVk6fKJNKYZ0tPmsDEmzxhTKNE29U33C6/+/ef8zHWpq1ZiiIZoy0qWjzWckGbooL1NHm8+q6VS7TrS1yzQNXTNpmC7Ky9I71Sf0Sf0pHWluU25Gqm6aNlJfuhBEt+w7pqojzRrhylRRbqZSDKnLPP/pUt/rD/W172Djab1TfVySodlhbhfsddaK95lwxUONcd+19OKLL+qOO+7QU089pVmzZuknP/mJNm7cqH379mnYsOAnmkW7a8mKgWmRtB+HO/hxoC3SoV7PN+aW63u//av+Wn8qrP3NHJ2vMx1dXoElkJ7BfrlZqfr68+/6HVbokPTPf3uxXv1LnT76LLyaANhDTrohhyMlrBPsrxhbqCe/Ot3va19zW4f+6YU9fn/iIdB2wV5n7TCYczBqTJj261mzZmnGjBn62c9+Jknq7u5WaWmpli1bpm9961tBt492kLFiYFqk7cfhDn4Md1+R1OzMSh20DpD87DR98aK8mA70A5C45o4vCjhk1l/XV7Dtgr3O2mEw52DUmBBDIzs6OrR7924tWLDAs8zhcGjBggXavn27z23a29vV0tLi9RctVgxMCzqsMZqDH6N0bIGuZzDbWE+2dcZ8oB+AxBVsyGyk2wV6nbXDYM54qzGug0xTU5O6urpUXFzstby4uFj19fU+t1m1apVcLpfnr7S01OflImHFwLSBtE5Gs404nGOLl0F6ADBQkQyZjXS7Q8dbbTGYM95qjOsgE4mVK1fK7XZ7/mpra6O2bysGpg2kdTKabcThHFu8DNIDgIGKZMhspNuNLhxii8Gc8VZjXAeZoUOHKiUlRQ0NDV7LGxoaNHz4cJ/bZGRkyOl0ev1Fy0DaiGN1nVEd/BilYwt0PYPZctzTNhuNtmkAycffa1/Pa1yk2wV6nbXifSZc8VZjXAeZ9PR0TZ8+XZs3b/Ys6+7u1ubNm1VRUWFJTWsWT/WcGNtjzrihnnbfwb7OcOuJ5r4iqfnXlVfqirGFYe9v5ugCzRiVH/Lle7qWfNXRW4ohvfj12ZpWmhd2TQDsISfdkDMzvF8buWJsYcDXvjWLp/p8LasYE3y7YK+zVrzPhCueaoz7rqUXX3xRS5Ys0dNPP62ZM2fqJz/5iV566SV9/PHH/c6d8SVWQyPDHbQY6+uM5uDHaB6bv31VN7XqnYPHZUqaPeb8i8Gh461qOtWuOvcZlbiydK7b7PebDtVNrdpx8LgMSSkOQ3XuM5pWlq8vjS/Sn/Y3ak/NSc//fdWR6jD01sfHdPhEmxZ+YbjXD4dVN7XqF9sP6a+ftWjSCKeunjhMf6lt1sHGVuVnp2lcca6n1p7aR+afr7PpVLs+qmuRYUgFQ9Jl6vxv4BTmZGj2mEIdOdmmTXuOSpJumnaRRuZn62d/3K+P6lrkkNQtUw7DUFqqQ9dOHq6i3Ay99mG9Pms+o+Ot7SpxZWvaqHxtO9CkevcZSVLBkIwL07E7ZBiG0lIcKs3Plis7TXsOn9CJ1g51nOtWelqKUozzBQ3PzVRuVppcWWn6sM6tYy3tchhSWopD2RkpGleUq2HOTBmGVO8+qxOtHTIkHTvdLveZTpndptJTDTkMQ6Z5/vjSUxyqPdmqzm4pI+V8HWc7u+RwSHnZ6SpxZetsZ5dOtHWoraNL7Z1dchhSZnqqznacU3uXqTSHNNyZpdbOczrT3iXjQk1nO7vUZUrOzFRlp6fqM/cZdXZLaYY0YXiuMlJTVN10Wm0dXUp1GCrIydCZjnM60dopU1JuRopG5mcrf0i6znZ06eDxVnWe61JXtykZhlxZaTrXbcrd1qGubnl+0ydFUla6Q/lD0jV7zFB1dZtqaDmjk62dOtLcprYLvyfTZUqmzv/uy8i8LM0oL5BpSnnZaWo506m39h3TqbOd6r5wuYzUFOVmpqj9wm8DdV34sZiMVIdMQ+ro7FZ7V7ey01I0ucSpnIw0HTh2WvXuM2o/1y0ZUm7G+Tfkto4uZaenqGBIhk60tqvjXLe6TVOmzt/uF7mydOx0u3IyUnXXleW6fHSBfvN+naobW3Ws5axOd5zzPNY27TkqwzA1wpWlnMw0faHEqYaWszrY2HrhN21MjSnK0VemlEiSfrH9kPbUNKsoN12XlebrwLHTOn32nJxZqWo5c045mam6adpFqnef1aY9R3W285xGX/gdmZH52Z7n4pGTZzzPcUme53bP/3vqHVM0RJeW5unoyTOe592Rk2d0/HS7TrZ2qMuUrpk0TCPzs7Xj4HEdaDil2pNtyunzOzI9rxElriwNzc1QqsPQuW7T81VI7+sP9bWv72tZONsFe5214n0mXLGsMWHaryXpZz/7mecH8S677DI99thjmjVrVkjbMv0aAAD7SaggMxAEGQAA7CchfkcGAAAgEIIMAACwLYIMAACwLYIMAACwLYIMAACwLYIMAACwLYIMAACwLYIMAACwLYIMAACwrfCmaNlQzw8Xt7S0WFwJAAAIVc/7drABBAkfZE6dOiVJKi0tDXJJAAAQb06dOiWXy+V3fcLPWuru7lZdXZ1yc3NlGIbV5cRUS0uLSktLVVtbmzRzpZLxmKXkPG6OOTmOWUrO4+aY+x+zaZo6deqUSkpK5HD4PxMm4T+RcTgcGjlypNVlDCqn05k0T4QeyXjMUnIeN8ecPJLxuDlmb4E+ienByb4AAMC2CDIAAMC2CDIJJCMjQw8//LAyMjKsLmXQJOMxS8l53Bxz8kjG4+aYI5fwJ/sCAIDExScyAADAtggyAADAtggyAADAtggyAADAtggyCerv//7vVVZWpszMTI0YMUK333676urqrC4rpg4dOqSlS5eqvLxcWVlZGjt2rB5++GF1dHRYXVpMff/739cVV1yh7Oxs5eXlWV1OzDz++OMaPXq0MjMzNWvWLO3cudPqkmJq69atuv7661VSUiLDMPTyyy9bXVJMrVq1SjNmzFBubq6GDRumG2+8Ufv27bO6rJh78sknNWXKFM+PwlVUVOh3v/ud1WUNqkcffVSGYWj58uURbU+QSVDz5s3TSy+9pH379ulXv/qVPv30U91yyy1WlxVTH3/8sbq7u/X000/rww8/1OrVq/XUU0/pX//1X60uLaY6Ojp066236p577rG6lJh58cUXtWLFCj388MPas2ePLr30Ui1cuFDHjh2zurSYaW1t1aWXXqrHH3/c6lIGxZYtW1RZWakdO3bo9ddfV2dnp6699lq1trZaXVpMjRw5Uo8++qh2796td999V/Pnz9cNN9ygDz/80OrSBsWuXbv09NNPa8qUKZHvxERSeOWVV0zDMMyOjg6rSxlUP/zhD83y8nKryxgUa9euNV0ul9VlxMTMmTPNyspKz/+7urrMkpISc9WqVRZWNXgkmZs2bbK6jEF17NgxU5K5ZcsWq0sZdPn5+eZ//dd/WV1GzJ06dcocP368+frrr5tXXXWV+c1vfjOi/fCJTBI4ceKEXnjhBV1xxRVKS0uzupxB5Xa7VVBQYHUZGICOjg7t3r1bCxYs8CxzOBxasGCBtm/fbmFliCW32y1JSfX87erq0oYNG9Ta2qqKigqry4m5yspKfeUrX/F6bkeCIJPAHnroIQ0ZMkSFhYWqqanRK6+8YnVJg+rAgQNas2aNvvGNb1hdCgagqalJXV1dKi4u9lpeXFys+vp6i6pCLHV3d2v58uWaM2eOLrnkEqvLibm9e/cqJydHGRkZuvvuu7Vp0yZNnjzZ6rJiasOGDdqzZ49WrVo14H0RZGzkW9/6lgzDCPj38ccfey7/z//8z3rvvff0hz/8QSkpKbrjjjtk2vCHnMM9bkk6evSo/vZv/1a33nqrvv71r1tUeeQiOWYgUVRWVuqDDz7Qhg0brC5lUFx88cWqqqrSO++8o3vuuUdLlizRRx99ZHVZMVNbW6tvfvObeuGFF5SZmTng/TGiwEYaGxt1/PjxgJcZM2aM0tPT+y0/cuSISktL9fbbb9vuI8twj7uurk5XX321Zs+erWeffVYOh/3yeiT39bPPPqvly5erubk5xtUNro6ODmVnZ+uXv/ylbrzxRs/yJUuWqLm5OSk+aTQMQ5s2bfI6/kR177336pVXXtHWrVtVXl5udTmWWLBggcaOHaunn37a6lJi4uWXX9ZNN92klJQUz7Kuri4ZhiGHw6H29navdcGkxqJIxEZRUZGKiooi2ra7u1uS1N7eHs2SBkU4x3306FHNmzdP06dP19q1a20ZYqSB3deJJj09XdOnT9fmzZs9b+Td3d3avHmz7r33XmuLQ9SYpqlly5Zp06ZNeuutt5I2xEjnH992fK0O1TXXXKO9e/d6Lbvzzjs1ceJEPfTQQ2GFGIkgk5Deeecd7dq1S1deeaXy8/P16aef6jvf+Y7Gjh1ru09jwnH06FFdffXVGjVqlH784x+rsbHRs2748OEWVhZbNTU1OnHihGpqatTV1aWqqipJ0rhx45STk2NtcVGyYsUKLVmyRJdffrlmzpypn/zkJ2ptbdWdd95pdWkxc/r0aR04cMDz/+rqalVVVamgoEBlZWUWVhYblZWVWrdunV555RXl5uZ6zn9yuVzKysqyuLrYWblypa677jqVlZXp1KlTWrdund566y299tprVpcWM7m5uf3Ofeo5nzOic6Ki10iFePH++++b8+bNMwsKCsyMjAxz9OjR5t13320eOXLE6tJiau3ataYkn3+JbMmSJT6P+c0337S6tKhas2aNWVZWZqanp5szZ840d+zYYXVJMfXmm2/6vF+XLFlidWkx4e+5u3btWqtLi6m77rrLHDVqlJmenm4WFRWZ11xzjfmHP/zB6rIG3UDarzlHBgAA2JY9TyAAAAAQQQYAANgYQQYAANgWQQYAANgWQQYAANgWQQYAANgWQQYAANgWQQYAANgWQQYAANgWQQZA3Orq6tIVV1yhm2++2Wu52+1WaWmpvv3tb0uS7rvvPk2fPl0ZGRm67LLLLKgUgFUIMgDiVkpKip599ln9/ve/1wsvvOBZvmzZMhUUFOjhhx/2LLvrrrt02223WVEmAAsx/RpAXJswYYIeffRRLVu2TPPnz9fOnTu1YcMG7dq1S+np6ZKkxx57TJLU2Nio999/38pyAQwyggyAuLds2TJt2rRJt99+u/bu3avvfve7uvTSS60uC0AcIMgAiHuGYejJJ5/UpEmT9MUvflHf+ta3rC4JQJzgHBkAtvDf//3fys7OVnV1tY4cOWJ1OQDiBEEGQNx7++23tXr1av3mN7/RzJkztXTpUpmmaXVZAOIAQQZAXGtra9PXvvY13XPPPZo3b56eeeYZ7dy5U0899ZTVpQGIAwQZAHFt5cqVMk1Tjz76qCRp9OjR+vGPf6x/+Zd/0aFDhyRJBw4cUFVVlerr63XmzBlVVVWpqqpKHR0dFlYOYDAYJp/PAohTW7Zs0TXXXKO33npLV155pde6hQsX6ty5c3rjjTc0b948bdmypd/21dXVGj169CBVC8AKBBkAAGBbfLUEAABsiyADAABsiyADAABsiyADAABsiyADAABsiyADAABsiyADAABsiyADAABsiyADAABsiyADAABsiyADAABs6/8BwL+XYzp5neMAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "model.eval()\n",
    "likelihood.eval()\n",
    "with torch.no_grad():\n",
    "    preds = likelihood(model(X))\n",
    "\n",
    "results_df = pd.DataFrame(dict(\n",
    "    y=y.cpu().numpy(),\n",
    "    X1=X.cpu().numpy(),\n",
    "    preds=preds.sample().median(axis=0).values.cpu().numpy()\n",
    "))\n",
    "\n",
    "results_df.plot.scatter(\"X1\", \"y\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "8e213a12-ada9-410c-ae2e-28beae4d8a29",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Axes: xlabel='X1', ylabel='preds'>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjIAAAGwCAYAAACzXI8XAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAM5pJREFUeJzt3X90VPWd//HXnclvyExCMiQiCeFHEMSCEOWHUCpK63p6uuuPursct1LL2VNZiD+7W2m3tae7bnDbo1ZUtJVC+10BtS3aqm1tqcD3CwFdNPUHQvkRJICB/CAz+QH5MbnfPzBThkySyWSSm0/yfJwz52TuvZ/P533v3Lnz4s6di2Xbti0AAAADuZwuAAAAIFYEGQAAYCyCDAAAMBZBBgAAGIsgAwAAjEWQAQAAxiLIAAAAYyU4XUB/a29v18mTJ5Weni7LspwuBwAARMG2bdXX12vMmDFyubo+7zLkg8zJkyeVl5fndBkAACAGFRUVGjt2bJfzh3yQSU9Pl3R+Q3g8HoerAQAA0QgEAsrLywt9jndlyAeZjq+TPB4PQQYAAMP0dFkIF/sCAABjEWQAAICxCDIAAMBYBBkAAGAsggwAADAWQQYAABiLIAMAAIxFkAEAAMYiyAAAAGMRZAAAgLEcDTLf+973ZFlW2GPKlCmh+efOndOKFSuUlZWlkSNH6tZbb9WpU6ccrBiDxZGqBr154LTKqxudLgUA4CDH/6+ladOm6Y9//GPoeULCX0u677779Nprr+mll16S1+vVypUrdcstt2jnzp1OlIpBoK6pRXdvKtOOg1WhaQsLfVqzZKa8aYkOVgYAcILjQSYhIUG5ubmdpvv9fq1bt04bN27UddddJ0lav369pk6dqt27d2vu3LkDXSoGgbs3lWnnoeqwaTsPVat407v6+bLZDlUFAHCK49fIHDx4UGPGjNGECRN0++2369ixY5KkvXv3qrW1VYsXLw4tO2XKFOXn56u0tLTL/pqbmxUIBMIeGBqOVDVox8EqBW07bHrQtrXjYBVfMwHAMORokJkzZ442bNig3/3ud1q7dq3Ky8v12c9+VvX19aqsrFRSUpIyMjLC2uTk5KiysrLLPktKSuT1ekOPvLy8fl4LDJSPa5u6nX+0hiADAMONo18t3XjjjaG/p0+frjlz5mjcuHF68cUXlZqaGlOfq1at0v333x96HggECDNDxLhRad3OL8gaMUCVAAAGC8e/WrpQRkaGJk+erEOHDik3N1ctLS2qq6sLW+bUqVMRr6npkJycLI/HE/bA0DDBN1ILC31yW1bYdLdlaWGhT+OzCTIAMNwMqiDT0NCgw4cP65JLLlFRUZESExO1devW0PwDBw7o2LFjmjdvnoNVwklrlszU/EnZYdPmT8rWmiUzHaoIAOAkR79a+sY3vqEvfelLGjdunE6ePKmHHnpIbrdbS5Yskdfr1bJly3T//fdr1KhR8ng8Ki4u1rx58/jF0jDmTUvUz5fNVnl1o47WNKogawRnYgBgGHM0yBw/flxLlixRTU2NfD6fFixYoN27d8vn80mSHnvsMblcLt16661qbm7WDTfcoKefftrJkjFIjM8mwAAAJMu2L/ot6xATCATk9Xrl9/u5XgYAAENE+/k9qK6RAQAA6A2CDAAAMBZBBgAAGIsgAwAAjEWQAQAAxiLIAAAAYxFkAACAsQgyAADAWAQZAABgLIIMAAAwFkEGAAAYiyADAACMRZABAADGIsgAAABjEWQAAICxCDIAAMBYBBkAAGAsggwAADAWQQYAABiLIAMAAIxFkAEAAMYiyAAAAGMRZAAAgLEIMgAAwFgEGQAAYCyCDAAAMBZBBgAAGIsgAwAAjEWQAQAAxiLIAAAAYxFkAACAsQgyAADAWAQZAABgLIIMAAAwFkEGAAAYiyADAACMRZABAADGIsgAAABjEWQAAICxCDIAAMBYBBkAAGAsggwAADAWQQYAABiLIAMAAIxFkAEAAMYiyAAAAGMRZAAAgLEIMgAAwFgEGQAAYCyCDAAAMBZBBgAAGIsgAwAAjEWQAQAAxiLIAAAAYxFkAACAsQgyAADAWAlOFwDzHKlq0Me1TSrIGqHx2SPitmw8xx2Iekww3NYXwPAzaILM6tWrtWrVKt1zzz16/PHHJUnnzp3TAw88oM2bN6u5uVk33HCDnn76aeXk5Dhb7DBV19SiuzeVacfBqtC0hYU+rVkyU960xJiXjee4/dmHSYbb+gIYvgbFV0tvv/22nn32WU2fPj1s+n333aff/OY3eumll7R9+3adPHlSt9xyi0NV4u5NZdp5qDps2s5D1Sre9G6flo3nuP3Zh0mG2/oCGL4cDzINDQ26/fbb9ZOf/ESZmZmh6X6/X+vWrdOjjz6q6667TkVFRVq/fr127dql3bt3d9lfc3OzAoFA2AN9d6SqQTsOVilo22HTg7atHQerVF7dGNOy8Ry3P/swyXBbXwDDm+NBZsWKFfriF7+oxYsXh03fu3evWltbw6ZPmTJF+fn5Ki0t7bK/kpISeb3e0CMvL6/fah9OPq5t6nb+0Zq/fjj2Ztl4jtuffZhkuK0vgOHN0SCzefNmvfPOOyopKek0r7KyUklJScrIyAibnpOTo8rKyi77XLVqlfx+f+hRUVER77KHpXGj0rqdX5D11wtJe7NsPMftzz5MMtzWF8Dw5liQqaio0D333KPnn39eKSkpces3OTlZHo8n7IG+m+AbqYWFPrktK2y627K0sNAX9ouY3iwbz3H7sw+TDLf1BTC8ORZk9u7dq9OnT2vWrFlKSEhQQkKCtm/frieeeEIJCQnKyclRS0uL6urqwtqdOnVKubm5zhQ9zK1ZMlPzJ2WHTZs/KVtrlszs07LxHLc/+zDJcFtfAMOXZdsXXRE4QOrr6/Xxxx+HTbvzzjs1ZcoUffOb31ReXp58Pp82bdqkW2+9VZJ04MABTZkyRaWlpZo7d25U4wQCAXm9Xvn9fs7OxEl5daOO1jRGdW+S3iwbz3EHoh4TDLf1BTB0RPv57ViQieTaa6/VlVdeGbqPzPLly/X6669rw4YN8ng8Ki4uliTt2rUr6j4JMgAAmCfaz+9Bc0O8SB577DG5XC7deuutYTfEAwAAkAbZGZn+wBkZAADME+3nt+P3kQEAAIgVQQYAABiLIAMAAIxFkAEAAMYiyAAAAGMRZAAAgLEIMgAAwFgEGQAAYCyCDAAAMBZBBgAAGIsgAwAAjEWQAQAAxiLIAAAAYxFkAACAsQgyAADAWAQZAABgLIIMAAAwFkEGAAAYiyADAACMRZABAADGIsgAAABjEWQAAICxCDIAAMBYBBkAAGAsggwAADAWQQYAABiLIAMAAIxFkAEAAMYiyAAAAGMRZAAAgLEIMgAAwFgEGQAAYCyCDAAAMBZBBgAAGIsgAwAAjEWQAQAAxiLIAAAAYxFkAACAsQgyAADAWAQZAABgLIIMAAAwFkEGAAAYiyADAACMRZABAADGIsgAAABjEWQAAICxCDIAAMBYBBkAAGAsggwAADAWQQYAABiLIAMAAIxFkAEAAMYiyAAAAGMRZAAAgLEIMgAAwFgJThcA9NaRqgZ9XNukgqwRGp89Im5tOpZxW5aCth33/mNt091yF86zbbvXNQCA6RwNMmvXrtXatWt19OhRSdK0adP03e9+VzfeeKMk6dy5c3rggQe0efNmNTc364YbbtDTTz+tnJwcB6uGU+qaWnT3pjLtOFgVmraw0Kc1S2bKm5YYc5tIy8Sz/1jbdLecLbvLmqOpAQCGCsu2bdupwX/zm9/I7XarsLBQtm3rZz/7mX7wgx/o3Xff1bRp07R8+XK99tpr2rBhg7xer1auXCmXy6WdO3dGPUYgEJDX65Xf75fH4+nHtUF/u2PdW9p5qFrBC3ZZt2Vp/qRs/XzZ7JjbRFomnv3H2qa75SR1WXM0NQDAYBft57ejZ2S+9KUvhT1/+OGHtXbtWu3evVtjx47VunXrtHHjRl133XWSpPXr12vq1KnavXu35s6dG7HP5uZmNTc3h54HAoH+WwEMmCNVDRHPPgRtWzsOVqm8ujHi1y49tbE//bsrfe0/1jY9LdeT7moAgKFk0FzsGwwGtXnzZjU2NmrevHnau3evWltbtXjx4tAyU6ZMUX5+vkpLS7vsp6SkRF6vN/TIy8sbiPLRzz6ubep2/tGaxpja9LRMX/uPtU20dfUkUg0AMJQ4HmTef/99jRw5UsnJybrrrru0ZcsWXX755aqsrFRSUpIyMjLCls/JyVFlZWWX/a1atUp+vz/0qKio6Oc1wEAYNyqt2/kFWZ3POkTTpqdl+tp/rG2irasnkWoAgKHE8SBz2WWXqaysTHv27NHy5cu1dOlS7du3L+b+kpOT5fF4wh4w3wTfSC0s9MltWWHT3ZalhYW+iF+fRNOmq2Xi1X+sbXparruae6oBAIYSx4NMUlKSJk2apKKiIpWUlGjGjBn60Y9+pNzcXLW0tKiuri5s+VOnTik3N9eZYuGoNUtmhi507TB/UrbWLJnZpzaRloln/7G26W657mqOpgYAGCoc/dVSJNddd53y8/P1ox/9SD6fT5s2bdKtt94qSTpw4ICmTJmi0tLSLi/2vRi/Whp6yqsbdbSmsVf3S4mmTccyCS5Lbe29u49Mf9XU03IXzpPU6xoAYLCK9vPb0SCzatUq3XjjjcrPz1d9fb02btyoRx55RL///e/1+c9/XsuXL9frr7+uDRs2yOPxqLi4WJK0a9euqMcgyAAAYB4jfn59+vRp3XHHHfrkk0/k9Xo1ffr0UIiRpMcee0wul0u33npr2A3xAAAApEH41VK8cUYGAADzRPv57fjFvgAAALEiyAAAAGMRZAAAgLEIMgAAwFgEGQAAYCyCDAAAMBZBBgAAGIsgAwAAjEWQAQAAxiLIAAAAY8UlyASDQZWVlenMmTPx6A4AACAqMQWZe++9V+vWrZN0PsR87nOf06xZs5SXl6dt27bFsz4AAIAuxRRkfvGLX2jGjBmSpN/85jcqLy/X/v37dd999+nb3/52XAsEAADoSkxBprq6Wrm5uZKk119/XbfddpsmT56sr33ta3r//ffjWiAAAEBXYgoyOTk52rdvn4LBoH73u9/p85//vCSpqalJbrc7rgUCAAB0JSGWRnfeeaf+/u//Xpdccoksy9LixYslSXv27NGUKVPiWiAAAEBXYgoy3/ve93TFFVeooqJCt912m5KTkyVJbrdbDz74YFwLBAAA6Ipl27btdBH9KRAIyOv1yu/3y+PxOF0OAACIQrSf31GfkXniiSeiHvzuu++OelkAAIBYRX1GZvz48WHPq6qq1NTUpIyMDElSXV2d0tLSNHr0aB05ciTuhcaKMzIAAJgn2s/vqH+1VF5eHno8/PDDuvLKK/XRRx+ptrZWtbW1+uijjzRr1iz9x3/8R1xWAAAAoCcxXSMzceJE/eIXv9DMmTPDpu/du1df/vKXVV5eHrcC+4ozMgAAmCfuZ2Qu9Mknn6itra3T9GAwqFOnTsXSJQAAQK/FFGSuv/56ff3rX9c777wTmrZ3714tX748dE8ZAACA/hZTkPnpT3+q3NxcXXXVVUpOTlZycrJmz56tnJwcPffcc/GuEQAAIKKYbojn8/n0+uuv6y9/+Yv2798vSZoyZYomT54c1+IAAAC6E1OQ6VBQUCDbtjVx4kQlJPSpKwAAgF6L6aulpqYmLVu2TGlpaZo2bZqOHTsmSSouLtbq1avjWiAAAEBXYgoyq1at0p///Gdt27ZNKSkpoemLFy/WCy+8ELfiAAAAuhPT90Evv/yyXnjhBc2dO1eWZYWmT5s2TYcPH45bcQAAAN2J6YxMVVWVRo8e3Wl6Y2NjWLABAADoTzEFmauuukqvvfZa6HlHeHnuuec0b968+FQGAADQg5i+Wvqv//ov3Xjjjdq3b5/a2tr0ox/9SPv27dOuXbu0ffv2eNcIAAAQUUxnZBYsWKA///nPamtr02c+8xm98cYbGj16tEpLS1VUVBTvGgEAACLq9RmZ1tZWff3rX9d3vvMd/eQnP+mPmgAAAKLS6zMyiYmJ+uUvf9kftQAAAPRKTF8t3XTTTXr55ZfjXAoAAEDvxHSxb2Fhob7//e9r586dKioq0ogRI8Lm33333XEpDgAAoDuWbdt2bxuNHz++6w4tS0eOHOlTUfEUCATk9Xrl9/vl8XicLgcAAEQh2s/vmM7IlJeXh/7uyEHcCA8AAAy0mK6RkaR169bpiiuuUEpKilJSUnTFFVfoueeei2dtAAAA3YrpjMx3v/tdPfrooyouLg7dybe0tFT33Xefjh07pu9///txLRIAACCSmK6R8fl8euKJJ7RkyZKw6Zs2bVJxcbGqq6vjVmBfcY0MAADmifbzO6avllpbW3XVVVd1ml5UVKS2trZYugQAAOi1mILMV77yFa1du7bT9B//+Me6/fbb+1wUAABANGK6RkY6f7HvG2+8oblz50qS9uzZo2PHjumOO+7Q/fffH1ru0Ucf7XuVAAAAEcQUZD744APNmjVLknT48GFJUnZ2trKzs/XBBx+EluMn2QAAoD/FFGTefPPNeNcBAADQazHfRwYAAMBpBBkAAGAsggwAADAWQQYAABiLIAMAAIxFkAEAAMYiyAAAAGMRZAAAgLFi/i8K4qGkpES/+tWvtH//fqWmpuqaa67RI488ossuuyy0zLlz5/TAAw9o8+bNam5u1g033KCnn35aOTk5DlYuHalq0Me1TSrIGqHx2SMcrWUouXC72rYd8e+utndPr0nHfLclBW2FluvNa9mX172rtt31eaSqQXvKa2VJmjMhKzT//PQaSZbmXjC9q/Eu3n7RtHdSvN9fvF+BocvRILN9+3atWLFCV199tdra2vStb31LX/jCF7Rv3z6NGHH+YHPffffptdde00svvSSv16uVK1fqlltu0c6dOx2pua6pRXdvKtOOg1WhaQsLfVqzZKa8aYmO1DQURNquXbl4e/f0mnTXd2Zaos40tXbZd3f1Rfu6d9X2P2+apn9/+cOIfdqytfx/3lHpkZqwvq4uyJQl6a2jZ8Kmz5uQpWf+qajH9ZUkb2qi/Gdbu2zvpHi/v3i/AkOfZdu27XQRHaqqqjR69Ght375dCxculN/vl8/n08aNG/XlL39ZkrR//35NnTpVpaWlof+wsjuBQEBer1d+v18ej6fPNd6x7i3tPFSt4AWbzW1Zmj8pWz9fNrvP/Q9XkbZrVy7e3j29Jn3pu7v6on3du2rrSU1Q4GxbxD4lRRXqLrSw0Nfr9Y3U3knxfn/xfgXMFe3n96C6Rsbv90uSRo0aJUnau3evWltbtXjx4tAyU6ZMUX5+vkpLSyP20dzcrEAgEPaIlyNVDdpxsKrTB0TQtrXjYJXKqxvjNtZw0tV27cqF27un12THX6pi7run+qJ53btre6apteu6exliJMW0vhe3d3Ifjvf7i/crMDwMmiDT3t6ue++9V/Pnz9cVV1whSaqsrFRSUpIyMjLCls3JyVFlZWXEfkpKSuT1ekOPvLy8uNX4cW1Tt/OP1nBgjEVP27UrR2sae2z7bsWZbud313eHvrzusa5brGJd3w5O7sPxfn/xfgWGh0ETZFasWKEPPvhAmzdv7lM/q1atkt/vDz0qKiriVKE0blRat/MLsriIMBY9bdeuFGSN6LHtzLzMmPvu0JfXPdZ1i1Ws69vByX043u8v3q/A8DAogszKlSv16quv6s0339TYsWND03Nzc9XS0qK6urqw5U+dOqXc3NyIfSUnJ8vj8YQ94mWCb6QWFvrktqyw6W7L0sJCH7+GiFFX27UrF27vnl6ThZN9MffdU33RvO7dtc1MS+y67kJfVPVeKJb1vbi9k/twvN9fvF+B4cHRIGPbtlauXKktW7boT3/6k8aPHx82v6ioSImJidq6dWto2oEDB3Ts2DHNmzdvoMuVJK1ZMjN0MWaH+ZOytWbJTEfqGSoibdeuXLy9e3pNuus786JfrnT1Wvblde+q7a9XLOiyzzVLZmrehKxOfc0uyNTsglGdps+bkBXV+krnf7XUXXsnxfv9xfsVGPoc/dXSv/zLv2jjxo165ZVXwu4d4/V6lZqaKklavny5Xn/9dW3YsEEej0fFxcWSpF27dkU1Rrx/tdShvLpRR2sauS9FnF24XSVF/Lur7d3Ta9IxP8Flqa3dDi3Xm9eyL697V22767O8ulG7j9R0uo9MeXWj9hypkS11eR+YrrZlxzr31N5J8X5/8X4FzBPt57ejQcbq4vT3+vXr9dWvflXSX2+It2nTprAb4nX11dLF+ivIAACA/mNEkBkIBBkAAMxj5H1kAAAAeoMgAwAAjEWQAQAAxiLIAAAAYxFkAACAsQgyAADAWAQZAABgLIIMAAAwFkEGAAAYiyADAACMRZABAADGIsgAAABjEWQAAICxCDIAAMBYBBkAAGAsggwAADAWQQYAABiLIAMAAIxFkAEAAMYiyAAAAGMRZAAAgLEIMgAAwFgEGQAAYCyCDAAAMBZBBgAAGIsgAwAAjEWQAQAAxiLIAAAAYxFkAACAsQgyAADAWAQZAABgLIIMAAAwFkEGAAAYiyADAACMRZABAADGIsgAAABjEWQAAICxCDIAAMBYBBkAAGAsggwAADAWQQYAABiLIAMAAIxFkAEAAMYiyAAAAGMRZAAAgLEIMgAAwFgEGQAAYCyCDAAAMBZBBgAAGIsgAwAAjEWQAQAAxiLIAAAAYxFkAACAsQgyAADAWAQZAABgLIIMAAAwVoKTg+/YsUM/+MEPtHfvXn3yySfasmWLbrrpptB827b10EMP6Sc/+Ynq6uo0f/58rV27VoWFhc4VPUQcqWrQx7VNKsgaIdu2Q3+Pzx4RNm989og+9d1df92Nc6SqQXvKa2VJmjMhK2IdXY3jtiwFbTuq8feU10iyNPeCMS7s50RdU9j8I1UNevW9kzpa3aQJvhGaPjZDQduW25LKKupU29iqyy/x6KNPAjp+pknjRqWpzbblsixlpiXJl54ctj4d/ZVXNykzLVGZaYnKTk+JWM+F6/PqeydV29iq66eO1qUZqdpTXqPqhhbJtiXLkmyptrFF/rMtkqTAuTZ5UhI1d0KWgna7qupbNDo9WZakjyoDyh6ZrC9OHxM2Zse2cVnSR5/UK3tkkr44fYxs2/50/BZdPzVHl2akhj1vb7f15oGqsOXPb0+prMKvo9WNykhNUGGuR3MnZHXa/yLZfuB0qM+Obd6b/TPSPtDX/RzA4GDZtm07Nfhvf/tb7dy5U0VFRbrllls6BZlHHnlEJSUl+tnPfqbx48frO9/5jt5//33t27dPKSkpUY0RCATk9Xrl9/vl8Xj6aU3MUdfUors3lWnHwaqI8zPTEnWmqTX0fGGhT2uWzJQ3LTGmviP19583TdO/v/xh2HId49iytfx/3lHpkZqwvq+ZmKW1txfJm5YY1Tjdj3+FHvzVe9p1OHyM2QWZcrtcncbukJ7sVn1zsMftEI3ZBaPUGgzq3Qp/l8tcPS5TiQmusDo9KQkKnGuLSw2RzMrPUILL0ltHz/TbGN25eH/7uKZRf/vkTvnPdn5tIy1/sUj7yrwJWbIshW3X3uznAAZGtJ/fjgaZC1mWFRZkbNvWmDFj9MADD+gb3/iGJMnv9ysnJ0cbNmzQP/7jP0bVL0Em3B3r3tLOQ9UKRvmyuy1L8ydl6+fLZselb7dlyZOaoMDZtrDlOsaR1GXIWljo08+Xze71OkQaP1LogfMu3t9mfv+Nbl+rnvbPaPeV3uznAAZGtJ/fg/YamfLyclVWVmrx4sWhaV6vV3PmzFFpaWmX7ZqbmxUIBMIeOO9IVYN2HKzqVQAI2rZ2HKxSeXVjXPoO2rbONLV2Wq5jnK5CjHQ+4Oz4S1Wv1yHS+BicLtzfth843eNr1d3+2Zv9Pdr9HMDgM2iDTGVlpSQpJycnbHpOTk5oXiQlJSXyer2hR15eXr/WaZKPa5tibnu0pvsDfF/67o13K5z5ygMD62hNo8qO1/Vq+YvFsk/2tJ8DGHwGbZCJ1apVq+T3+0OPiooKp0saNMaNSou5bUFW9xdD9qXv3piZlzkg48BZBVkjdOXYjF4tf7FY9sme9nMAg8+gDTK5ubmSpFOnToVNP3XqVGheJMnJyfJ4PGEPnDfBN1ILC31yW1bUbdyWpYWFvh5/1RFt327LUmZaYqflOsZZWOjrsu3CQp8WTvb1eh0ijY/B6cL97XOXje7xtepu/+zN/h7tfg5g8Bm0QWb8+PHKzc3V1q1bQ9MCgYD27NmjefPmOViZ2dYsmRm6qDaSiz845k/K1polM2PuO1J/v16xoNNyHeOsWTJT8yZkder7molZoTqiGaen8a+Z2HmM2QWZEcfukJ7s7nJeb80uGKVZed5ul7l6XGanOj0p/XvHhFn5GZpd4NxZr4v3t1+vWCBvatdhpqf9M9K+Mm9CVqft2pv9HMDg4uivlhoaGnTo0CFJ0syZM/Xoo49q0aJFGjVqlPLz8/XII49o9erVYT+/fu+99/j5dRyUVzfqaE1j6FR6x9/js0eEzYvlX6gXt++qv+7GKa9u1O4jNd3eR6arcRJcltra7ajG33OkRrYUdt+WC/s5ceZs2Pzy6ka99t5JHalq1ATfCM3Iy1Bbu60El6U/V9SpuqFF08acv49MRW2TxmWNULDdlmVJo0YkKXtk+H1kLuwvMy1RmSOSlDUyOWI9F67Pa++dVHVDi66fOlpjM9O050iNqhqaZUmyJVmSzjS26ExTq1yWFDjbppEpCZo3MUvBdltVDc3yjUyW22Xpw5P+TveRuXDbJLgsfXgyELovjKRPx2/W9VNzNDYzLey5JG396HTY8h3bs6yiTuWfruuknHTN/TQ49rS//d+DVaE+O7Z5b/bPSPtAX/dzAP3LiJ9fb9u2TYsWLeo0fenSpdqwYUPohng//vGPVVdXpwULFujpp5/W5MmTox6DIAMAgHmMCDIDgSADAIB5jL+PDAAAQE8IMgAAwFgEGQAAYCyCDAAAMBZBBgAAGIsgAwAAjEWQAQAAxiLIAAAAYxFkAACAsQgyAADAWAQZAABgLIIMAAAwFkEGAAAYiyADAACMRZABAADGIsgAAABjEWQAAICxCDIAAMBYBBkAAGAsggwAADAWQQYAABiLIAMAAIxFkAEAAMYiyAAAAGMRZAAAgLEIMgAAwFgEGQAAYCyCDAAAMBZBBgAAGIsgAwAAjEWQAQAAxiLIAAAAYxFkAACAsQgyAADAWAQZAABgLIIMAAAwFkEGAAAYiyADAACMRZABAADGIsgAAABjEWQAAICxCDIAAMBYBBkAAGAsggwAADAWQQYAABiLIAMAAIxFkAEAAMYiyAAAAGMRZAAAgLEIMgAAwFgEGQAAYCyCDAAAMBZBBgAAGIsgAwAAjEWQAQAAxiLIAAAAYyU4XUA0nnrqKf3gBz9QZWWlZsyYoTVr1mj27NlOlzWgjlQ1aE95jSRLcydkaXz2iE7zP65tUkHWiE7z+jJmNH1eWNulGakK2naXbTr6dFuWTtQ1hdbnWE2j3jxwWpZl6bKcdM2JsI4X12TbdsS/L2wXaR2+/+sPtetItRZM8umzhdkqO16nSzyp8nmS9f/+UqWdh6uVNSJZd107Ue3ttt48UCWXpMLcdF2akaKgLZ0OnNPuIzVqaG5TXkaaMkckKTs9OTTfbVl6c/9pHTgV0CXeVHlSEnXgVL0kWzfPHKuicZl68k+H9NEnATU0t+lsa5sSXC6lJblVlD9KQbtd+z4JqKm5TY0tQU29xKPLRqfrV+8eV93ZNlmSZuVnKCnBpX2fBJSa6FJGWpKa29oVONsqSUp0u3S2LagEy1KC2yXJVlvQVrDd1siURJ2uP6fmNluSlOw+v70SE9xKdFlKTnSrNXi+r3Zbsm2p/dNtaklyWZ8u75YSXC61Bttl25baZctuP7+sJSkl0VLWyGQluV3yn21V/bk2BW1b9kV9SlJKgqXWoK12W0pyW0pPSVB1Y2tofoJLavu0gUvSlNx0na4/pzNNrbJtyeWS3C6ppU2yL6jV0vnnKQmW2iW1ttmd1iVoq5Nkt6XRnmTNHZ+lE3VnVXGmSS1BW/6mFp1r69wgMy1RUy9J19mWoAJnW+U/26q6xlYFL6gjwS0lJ7jV3NZ+fhtaUrsttbfbsnW+fslSstslX3qKkhIsNTS3qbahRU2t7bIkJSdYcrssNbe1a2RygtwuS3VnW8/377JkS8pIS9IYT6pOBs4q0e3Sktn5+psrcvXknw7qo0/qZclWUqJbX7g8V1kjkvRy2Umdaw0qNcmlEckJystMky3peO1ZNba0SraU603VzbMu1aUZqfo/pUe1+0iNUhITVJSfqZqmZjWeCyo9NUGBs61KT0nUzbMu1ckzZz/tu03js0eG2p8/Bkgn6s7JkjRnQpZs29ae8tqw56++94mOVjdqgm+Epo/16kTdWZ0/zqToRN05Vdc3q7apWbZt6fqpo3VpRqr2lNfqL5X1Ol7XpPTkBN08a6w+W+iTJG0/cPr8+92bIl96ityfvvYdx5ALx4/2ONrTsbm7dj0dX/vjuB5vg6FGy7btCG/hweOFF17QHXfcoWeeeUZz5szR448/rpdeekkHDhzQ6NGje2wfCATk9Xrl9/vl8XgGoOL4qmtq0b88/452Ha4Jmz5vQpae+aci2bJ196Yy7ThYFZq3sNCnNUtmypuWGPOY0fTZVW2R2kTqsyfXTMzS2tuLet1+YaFP/3nTNP37yx+GLT8ld6T2VzZEPT6AoWFkkiWXy63Aubao21x4/Imkp2Nzd+16Or5Gewx20kDUGO3n96APMnPmzNHVV1+tJ598UpLU3t6uvLw8FRcX68EHH+yxvelB5o51b3X54b3w039l7DxUreAFL6PbsjR/UrZ+viy2s1Z3rHsrqj67q+3iNpH6jMbCQl+v27stS57UBAXOtvV6PADo0HH8iaSnY3N37Xo6vkZ7DHbSQNQY7ef3oL5GpqWlRXv37tXixYtD01wulxYvXqzS0tKIbZqbmxUIBMIepjpS1dBtUNhxsEo7DlZ1+rAO2rZ2HKxSeXVjzGP21GdPtV3YZsdfTkfsMxqxtA/ats40tRJiAPRJV8fRaI7N3bXr7vga7THYSYOtxkEdZKqrqxUMBpWTkxM2PScnR5WVlRHblJSUyOv1hh55eXkDUWq/+Li2qU/tj9b0fmfqacyOPntT27sVdb2uI57tASBWkY6j0Rz/Yml3tKYx6mOwkwZbjYM6yMRi1apV8vv9oUdFRYXTJcVs3Ki0PrUvyOr9hVc9jdnRZ29qm5mX0es64tkeAGIV6TgazfEvlnYFWSOiPgY7abDVOKiDTHZ2ttxut06dOhU2/dSpU8rNzY3YJjk5WR6PJ+xhqgm+kaHrYCJZWOjTwkKf3JYVNt1tWVpY6IvpCvKOMXvqs6faLmyzcPLoiH1GI5b2bstSZlpiTOMBQIeujqPRHJu7a9fd8TXaY7CTBluNgzrIJCUlqaioSFu3bg1Na29v19atWzVv3jwHKxs4a5bM1DUTszpNnzchS2uWzNSaJTM1f1J22Lz5k7K1ZsnMPo0ZTZ9d1RapTaQ+e3LNxKyY2s+flK1fr1jQafnLc9N7NT6AoWFkkiVPSu/uNnLh8SeSno7N3bXr6fjaH8f1eBtMNQ76Xy298MILWrp0qZ599lnNnj1bjz/+uF588UXt37+/07UzkZj+q6UO5dWN2nOkRrYU8V4F5dWNOlrTGNff8kfb54W1jc1MVVt71/eR6egzwWXpxJmzofU5fqZJf/rotCxLmtzNfWQurElSxL8vbBdpHf7z1X36f4eqtGCST5+7zKd3jp3RGG+qstOTtfNgtXYeqtKoT+8jI0lbPzottyVNykkPrV91fbNKD9eovrlV+Znn7yOTNTI5ND/BZWnb/tP6qDKgMd5UeVMTtb+yXvan95G5qmCUnvzTQe07GVBjc5uaWtuU6HIpNcmtonGjZNu2PjzpV1NzUA0tbZp6iUdTcj365d6KsPvIpCS69eFJv1IS3cpIS1JLW7v8Z1sknb+PzLm2drmt83/bshUM2mprt5WekqhTPdxHpi3YLn8v7yNjy1Z7xPvIuOU/29Kr+8h4UhJU1c19ZKZekq5TgYvvI2Oppc3u5j4yllrb2nt3H5kJ2TpR16SK2ia1Bm3VdXMfmcsvSdfZlvPbzX+2pcf7yFifjt3efr4/16cbtuM+MskJLtU3t6m2obnX95G51JuqE/6/3kfmxs9coqf+dFAfngzIJSkx0aUvXJ4rX3qytrxz4tP7yLiVluxWfub5rw4qOu4jI0s5nhTdPOtSjc1M0/+UHtWuw9Xn7yMzLlO1jS1qONcW+rXgyJQE3TzrUlX6z33ad5sKPr2PzNjMtNAx4PiZs6H7tkjS7iM1Yc9ffe+kyqvO30dmRl5G6JgxNjNVx8+cVU1Ds840tihoS9dPHa2xmWnafaRGh07Vq+JMk0ZedB+Z/3uwKuz9nuCyQseri8eP9jja07G5u3Y9HV/747geb/1Z45D5+bUkPfnkk6Eb4l155ZV64oknNGfOnKjaDpUgAwDAcDKkgkxfEGQAADDPkLiPDAAAQHcIMgAAwFgEGQAAYCyCDAAAMBZBBgAAGIsgAwAAjEWQAQAAxiLIAAAAYxFkAACAsXr3v2gZqOPGxYFAwOFKAABAtDo+t3v6DwiGfJCpr6+XJOXl5TlcCQAA6K36+np5vd4u5w/5/2upvb1dJ0+eVHp6uizLcrqcfhUIBJSXl6eKioph8/9KDcd1lobnerPOw2OdpeG53qxz53W2bVv19fUaM2aMXK6ur4QZ8mdkXC6Xxo4d63QZA8rj8QybN0KH4bjO0vBcb9Z5+BiO6806h+vuTEwHLvYFAADGIsgAAABjEWSGkOTkZD300ENKTk52upQBMxzXWRqe6806Dx/Dcb1Z59gN+Yt9AQDA0MUZGQAAYCyCDAAAMBZBBgAAGIsgAwAAjEWQGaL+9m//Vvn5+UpJSdEll1yir3zlKzp58qTTZfWro0ePatmyZRo/frxSU1M1ceJEPfTQQ2ppaXG6tH718MMP65prrlFaWpoyMjKcLqffPPXUUyooKFBKSormzJmjt956y+mS+tWOHTv0pS99SWPGjJFlWXr55ZedLqlflZSU6Oqrr1Z6erpGjx6tm266SQcOHHC6rH63du1aTZ8+PXRTuHnz5um3v/2t02UNqNWrV8uyLN17770xtSfIDFGLFi3Siy++qAMHDuiXv/ylDh8+rC9/+ctOl9Wv9u/fr/b2dj377LP68MMP9dhjj+mZZ57Rt771LadL61ctLS267bbbtHz5cqdL6TcvvPCC7r//fj300EN65513NGPGDN1www06ffq006X1m8bGRs2YMUNPPfWU06UMiO3bt2vFihXavXu3/vCHP6i1tVVf+MIX1NjY6HRp/Wrs2LFavXq19u7dq//93//Vddddp7/7u7/Thx9+6HRpA+Ltt9/Ws88+q+nTp8feiY1h4ZVXXrEty7JbWlqcLmVA/fd//7c9fvx4p8sYEOvXr7e9Xq/TZfSL2bNn2ytWrAg9DwaD9pgxY+ySkhIHqxo4kuwtW7Y4XcaAOn36tC3J3r59u9OlDLjMzEz7ueeec7qMfldfX28XFhbaf/jDH+zPfe5z9j333BNTP5yRGQZqa2v1/PPP65prrlFiYqLT5Qwov9+vUaNGOV0G+qClpUV79+7V4sWLQ9NcLpcWL16s0tJSBytDf/L7/ZI0rN6/wWBQmzdvVmNjo+bNm+d0Of1uxYoV+uIXvxj23o4FQWYI++Y3v6kRI0YoKytLx44d0yuvvOJ0SQPq0KFDWrNmjb7+9a87XQr6oLq6WsFgUDk5OWHTc3JyVFlZ6VBV6E/t7e269957NX/+fF1xxRVOl9Pv3n//fY0cOVLJycm66667tGXLFl1++eVOl9WvNm/erHfeeUclJSV97osgY5AHH3xQlmV1+9i/f39o+X/913/Vu+++qzfeeENut1t33HGHbANv5Nzb9ZakEydO6G/+5m9022236Z//+Z8dqjx2sawzMFSsWLFCH3zwgTZv3ux0KQPisssuU1lZmfbs2aPly5dr6dKl2rdvn9Nl9ZuKigrdc889ev7555WSktLn/vgvCgxSVVWlmpqabpeZMGGCkpKSOk0/fvy48vLytGvXLuNOWfZ2vU+ePKlrr71Wc+fO1YYNG+RymZfXY3mtN2zYoHvvvVd1dXX9XN3AamlpUVpamn7xi1/opptuCk1funSp6urqhsWZRsuytGXLlrD1H6pWrlypV155RTt27ND48eOdLscRixcv1sSJE/Xss886XUq/ePnll3XzzTfL7XaHpgWDQVmWJZfLpebm5rB5PUnojyLRP3w+n3w+X0xt29vbJUnNzc3xLGlA9Ga9T5w4oUWLFqmoqEjr1683MsRIfXuth5qkpCQVFRVp69atoQ/y9vZ2bd26VStXrnS2OMSNbdsqLi7Wli1btG3btmEbYqTz+7eJx+poXX/99Xr//ffDpt15552aMmWKvvnNb/YqxEgEmSFpz549evvtt7VgwQJlZmbq8OHD+s53vqOJEycadzamN06cOKFrr71W48aN0w9/+ENVVVWF5uXm5jpYWf86duyYamtrdezYMQWDQZWVlUmSJk2apJEjRzpbXJzcf//9Wrp0qa666irNnj1bjz/+uBobG3XnnXc6XVq/aWho0KFDh0LPy8vLVVZWplGjRik/P9/ByvrHihUrtHHjRr3yyitKT08PXf/k9XqVmprqcHX9Z9WqVbrxxhuVn5+v+vp6bdy4Udu2bdPvf/97p0vrN+np6Z2ufeq4njOma6Li90MqDBbvvfeevWjRInvUqFF2cnKyXVBQYN9111328ePHnS6tX61fv96WFPExlC1dujTiOr/55ptOlxZXa9assfPz8+2kpCR79uzZ9u7du50uqV+9+eabEV/XpUuXOl1av+jqvbt+/XqnS+tXX/va1+xx48bZSUlJts/ns6+//nr7jTfecLqsAdeXn19zjQwAADCWmRcQAAAAiCADAAAMRpABAADGIsgAAABjEWQAAICxCDIAAMBYBBkAAGAsggwAADAWQQYAABiLIANg0AoGg7rmmmt0yy23hE33+/3Ky8vTt7/9bUnS3XffraKiIiUnJ+vKK690oFIATiHIABi03G63NmzYoN/97nd6/vnnQ9OLi4s1atQoPfTQQ6FpX/va1/QP//APTpQJwEH879cABrXJkydr9erVKi4u1nXXXae33npLmzdv1ttvv62kpCRJ0hNPPCFJqqqq0nvvvedkuQAGGEEGwKBXXFysLVu26Ctf+Yref/99ffe739WMGTOcLgvAIECQATDoWZaltWvXaurUqfrMZz6jBx980OmSAAwSXCMDwAg//elPlZaWpvLych0/ftzpcgAMEgQZAIPerl279Nhjj+nVV1/V7NmztWzZMtm27XRZAAYBggyAQa2pqUlf/epXtXz5ci1atEjr1q3TW2+9pWeeecbp0gAMAgQZAIPaqlWrZNu2Vq9eLUkqKCjQD3/4Q/3bv/2bjh49Kkk6dOiQysrKVFlZqbNnz6qsrExlZWVqaWlxsHIAA8GyOT8LYJDavn27rr/+em3btk0LFiwIm3fDDTeora1Nf/zjH7Vo0SJt3769U/vy8nIVFBQMULUAnECQAQAAxuKrJQAAYCyCDAAAMBZBBgAAGIsgAwAAjEWQAQAAxiLIAAAAYxFkAACAsQgyAADAWAQZAABgLIIMAAAwFkEGAAAY6/8DpA+A3i6wp7oAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "results_df.plot.scatter(\"X1\", \"preds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "978552e2-24da-478e-9c64-2317cd5fbfb9",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.10.17"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}