zero_inflated_negative_binomial_regression.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "f0ab5256-b073-44d0-9724-e98d7109d27a",
   "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "af549afc-9003-4e36-813f-74c97479e8a5",
   "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 NegativeBinomialLikelihood(_OneDimensionalLikelihood):\n",
    "    def __init__(\n",
    "        self,\n",
    "        batch_shape: torch.Size = torch.Size([])\n",
    "    ) -> None:\n",
    "        super().__init__()        \n",
    "        \n",
    "        self.raw_r = torch.nn.Parameter(torch.ones(*batch_shape, 1))\n",
    "        r_constraint = Positive()\n",
    "        self.register_constraint(\"raw_r\", r_constraint)\n",
    "\n",
    "        self.raw_mu = torch.nn.Parameter(torch.ones(*batch_shape, 1))\n",
    "        mu_constraint = Positive()\n",
    "        self.register_constraint(\"raw_mu\", mu_constraint)\n",
    "\n",
    "    @property\n",
    "    def mu(self) -> Tensor:\n",
    "        return self.raw_mu_constraint.transform(self.raw_mu)\n",
    "\n",
    "    @mu.setter\n",
    "    def mu(self, value: Tensor) -> None:\n",
    "        self._set_mu(value)\n",
    "        \n",
    "    def _set_mu(self, value: Tensor) -> None:\n",
    "        if not torch.is_tensor(value):\n",
    "            value = torch.as_tensor(value).to(self.raw_mu)\n",
    "        self.initialize(raw_mu=self.raw_mu_constraint.inverse_transform(value))\n",
    "\n",
    "    @property\n",
    "    def r(self) -> Tensor:\n",
    "        return self.raw_r_constraint.transform(self.raw_r)\n",
    "\n",
    "    @r.setter\n",
    "    def r(self, value: Tensor) -> None:\n",
    "        self._set_r(value)\n",
    "        \n",
    "    def _set_r(self, value: Tensor) -> None:\n",
    "        if not torch.is_tensor(value):\n",
    "            value = torch.as_tensor(value).to(self.raw_r)\n",
    "        self.initialize(raw_r=self.raw_r_constraint.inverse_transform(value))\n",
    "        \n",
    "    def forward(self, function_samples: Tensor, *args: Any, **kwargs: Any) -> base_distributions.Poisson:\n",
    "        mu = torch.exp(function_samples) #* self.mu\n",
    "        probs = self.r / (mu + self.r)\n",
    "        return base_distributions.NegativeBinomial(total_count = self.r, probs=probs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "c8d069c7-4032-435f-b5a9-3bbb8c6479b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "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 = 500\n",
    "r = 1\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.NegativeBinomial(total_count = r, probs=probs).sample().double().cuda().flatten()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "a2ae98dc-cf48-4bd3-aa3c-c8528632e304",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100\n",
      "1.5387012438176348\n",
      "200\n",
      "1.5383521260207635\n",
      "300\n",
      "1.538306141805378\n",
      "400\n",
      "1.53826458000918\n",
      "500\n",
      "1.5382279044325708\n",
      "600\n",
      "1.538195906883475\n",
      "700\n",
      "1.5381681523700128\n",
      "800\n",
      "1.5381441471408097\n",
      "900\n",
      "1.5381234096662817\n",
      "1000\n",
      "1.5381054991912098\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "\n",
    "class GPModel(ApproximateGP):\n",
    "    def __init__(self, inducing_points, feature_extractor):\n",
    "        variational_distribution = gpytorch.variational.NaturalVariationalDistribution(inducing_points.size(0))\n",
    "        variational_strategy = gpytorch.variational.VariationalStrategy(\n",
    "            self, inducing_points, variational_distribution, learn_inducing_locations=True\n",
    "        )\n",
    "        super(GPModel, self).__init__(variational_strategy)\n",
    "        self.mean_module = gpytorch.means.ZeroMean()\n",
    "        self.covar_module = gpytorch.kernels.ScaleKernel(\n",
    "            gpytorch.kernels.LinearKernel(),\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",
    "M = 100\n",
    "inducing_points = torch.normal(mean=0, std =1, size=(1,N)).double().cuda().flatten()\n",
    "model = GPModel(inducing_points=inducing_points, feature_extractor=None).double()\n",
    "likelihood = NegativeBinomialLikelihood().double()\n",
    "\n",
    "if torch.cuda.is_available():\n",
    "    model = model.cuda()\n",
    "    likelihood = likelihood.cuda()\n",
    "\n",
    "model.train()\n",
    "#model.covar_module.base_kernel.initialize(lengthscale=5)\n",
    "likelihood.train()\n",
    "likelihood.initialize(r=1.5)\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 = 1000\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": 33,
   "id": "d7a2f16d-26dd-4789-8e96-37d29d7da890",
   "metadata": {},
   "outputs": [],
   "source": [
    "model.eval()\n",
    "likelihood.eval()\n",
    "with torch.no_grad():\n",
    "    preds = likelihood(model(X))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "64c6afe2-5f46-463d-95d8-c71da4b71e2a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Axes: xlabel='X1', ylabel='y'>"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjIAAAGwCAYAAACzXI8XAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAPYlJREFUeJzt3X14FPW9///XkoSQkGRDbiBEwn28L4hUaITmK0pF/B5/pVqvlraKLZfnlCJeantZOaennvbUg7U9l7bUYs91LNTrCFrbYq3fVo+iJC1FsEqKWqEEgiAxkASym2wgCcn+/sBddzd7M7vZ3ZnZfT6ua6+LzM585j2zk+yb2ZnXOrxer1cAAAA2NMrsAgAAABJFIwMAAGyLRgYAANgWjQwAALAtGhkAAGBbNDIAAMC2aGQAAIBt5ZpdQKoNDQ2ptbVVxcXFcjgcZpcDAAAM8Hq96u7uVnV1tUaNinzeJeMbmdbWVtXU1JhdBgAASMDRo0c1adKkiM9nfCNTXFws6dyOKCkpMbkaAABghNvtVk1Njf99PJKMb2R8HyeVlJTQyAAAYDOxLgvhYl8AAGBbNDIAAMC2aGQAAIBt0cgAAADbopEBAAC2RSMDAABsi0YGAADYFo0MAACwLRoZAABgWzQyAADAtkxtZDZs2KBZs2b5vz6grq5Of/jDH/zPX3XVVXI4HEGPr371qyZWDMR2qL1Hr+4/oZYOj9mlAEDGM/W7liZNmqQHH3xQtbW18nq9+sUvfqFPf/rT2rNnjy655BJJ0u23367vfve7/mUKCwvNKheIqqu3X3duaVLjgXb/tPraSq1fPkfOwjwTKwOAzGXqGZkbbrhB119/vWpra3X++efrgQceUFFRkV577TX/PIWFhaqqqvI/+OJHWNWdW5q0o7kjaNqO5g6t2bLHpIoAIPNZ5hqZwcFBPfXUU/J4PKqrq/NPf/LJJ1VRUaFLL71Ua9euVW9vb9Rx+vr65Ha7gx5Aqh1q71HjgXYNer1B0we9XjUeaOdjJgBIEVM/WpKkt956S3V1dTpz5oyKioq0detWXXzxxZKkL3zhC5oyZYqqq6u1d+9effOb39T+/fv1m9/8JuJ469at03e+8510lQ9Ikt47Gb3BPtzp0bSKsWmqBgCyh8PrDfkvZJr19/fryJEjcrlc+tWvfqX//u//VkNDg7+ZCfTKK6/ommuuUXNzs2bMmBF2vL6+PvX19fl/drvdqqmpkcvl4mMppMyh9h5d/Z8NEZ9/9RtX0cgAQBzcbrecTmfM92/TP1oaPXq0Zs6cqblz52rdunWaPXu2fvSjH4Wdd/78+ZKk5ubmiOPl5+f774LyPYBUm15ZpPraSuU4HEHTcxwO1ddW0sQAQIqY3siEGhoaCjqjEqipqUmSNHHixDRWBBizfvkcLZhZETRtwcwKrV8+x6SKACDzmXqNzNq1a7V06VJNnjxZ3d3d2rx5s7Zv364XX3xRBw8e1ObNm3X99dervLxce/fu1d133636+nrNmjXLzLKBsJyFeXpi5Ty1dHh0uNOjqeVjORMDAClmaiNz4sQJ3Xrrrfrggw/kdDo1a9Ysvfjii/rUpz6lo0eP6uWXX9Yjjzwij8ejmpoa3XTTTfrWt75lZslATNMqaGAAIF1Mv9g31YxeLAQAAKzDNhf7AgAAJIpGBgAA2BaNDAAAsC0aGQAAYFs0MgAAwLZoZAAAgG3RyAAAANuikQEAALZlarIvYNSh9h69d7KX2H8AQBAaGVhaV2+/7tzSpMYD7f5p9bWVWr98jpyFeSZWBgCwAj5agqXduaVJO5o7gqbtaO7Qmi17TKoIAGAlNDKwrEPtPWo80K7BkK8DG/R61XigXS0dHpMqAwBYBY0MLOu9k71Rnz/cSSMDANmORgaWNaWsMOrzU8u56BcAsh2NDCxremWR6msrleNwBE3PcThUX1vJ3UsAABoZWNv65XO0YGZF0LQFMyu0fvkckyoCAFgJt1/D0pyFeXpi5Ty1dHh0uNNDjgwAIAiNDGxhWgUNDABgOD5aAgAAtkUjAwAAbItGBgAA2BaNDAAAsC0aGQAAYFs0MgAAwLZoZAAAgG3RyAAAANuikQEAALZFIwMAAGyLRgYAANgWjQwAALAtGhkAAGBbNDIAAMC2aGQAAIBt0cgAAADbopEBAAC2RSMDAABsi0YGAADYlqmNzIYNGzRr1iyVlJSopKREdXV1+sMf/uB//syZM1q9erXKy8tVVFSkm266ScePHzexYgAAYCWmNjKTJk3Sgw8+qDfeeEN/+ctfdPXVV+vTn/603nnnHUnS3Xffrd/97nd65pln1NDQoNbWVt14441mlgwAACzE4fV6vWYXEaisrEw/+MEP9NnPflaVlZXavHmzPvvZz0qS9u3bp4suukg7d+7UJz7xCUPjud1uOZ1OuVwulZSUpLJ0AACQJEbfvy1zjczg4KCeeuopeTwe1dXV6Y033tDAwIAWL17sn+fCCy/U5MmTtXPnzojj9PX1ye12Bz0AAEBmMr2Reeutt1RUVKT8/Hx99atf1datW3XxxRerra1No0ePVmlpadD8EyZMUFtbW8Tx1q1bJ6fT6X/U1NSkeAsAAIBZTG9kLrjgAjU1NWnXrl1atWqVVqxYob/97W8Jj7d27Vq5XC7/4+jRo0msFgAAWEmu2QWMHj1aM2fOlCTNnTtXr7/+un70ox/pc5/7nPr7+9XV1RV0Vub48eOqqqqKOF5+fr7y8/NTXTYAALAA08/IhBoaGlJfX5/mzp2rvLw8bdu2zf/c/v37deTIEdXV1ZlYIQAAsApTz8isXbtWS5cu1eTJk9Xd3a3Nmzdr+/btevHFF+V0OrVy5Urdc889KisrU0lJidasWaO6ujrDdywBAIDMZmojc+LECd1666364IMP5HQ6NWvWLL344ov61Kc+JUl6+OGHNWrUKN10003q6+vTkiVL9NOf/tTMkgEAgIVYLkcm2ciRie5Qe4/eO9mrqeVjNa1irNnlAAAgyfj7t+kX+8IcXb39unNLkxoPtPun1ddWav3yOXIW5plYGQAAxlnuYl+kx51bmrSjuSNo2o7mDq3ZssekigAAiB+NTBY61N6jxgPtGgz5VHHQ61XjgXa1dHhMqgwAgPjQyGSh9072Rn3+cCeNDADAHmhkstCUssKoz08t56JfAIA90MhkoemVRaqvrVSOwxE0PcfhUH1tJXcvAQBsg0YmS61fPkcLZlYETVsws0Lrl88xqSIAAOLH7ddZylmYpydWzlNLh0eHOz3kyAAAbIlGJstNq6CBAQDYFx8tAQAA26KRAQAAtkUjAwAAbItGBgAA2BaNDAAAsC0aGQAAYFs0MgAAwLZoZAAAgG3RyAAAANuikQEAALZFIwMAAGyLRgYAANgWjQwAALAtGhkAAGBbNDIAAMC2aGQAAIBt0cgAAADbopEBAAC2RSMDAABsi0YGAADYFo0MAACwLRoZAABgWzQyAADAtmhkAACAbdHIAAAA26KRAQAAtkUjAwAAbItGBgAA2Japjcy6det0xRVXqLi4WOPHj9eyZcu0f//+oHmuuuoqORyOoMdXv/pVkyoGAABWYmoj09DQoNWrV+u1117TSy+9pIGBAV177bXyeDxB891+++364IMP/I+HHnrIpIoBAICV5Jq58hdeeCHo502bNmn8+PF64403VF9f759eWFioqqqqdJcHAAAszlLXyLhcLklSWVlZ0PQnn3xSFRUVuvTSS7V27Vr19vZGHKOvr09utzvoAQAAMpOpZ2QCDQ0N6a677tKCBQt06aWX+qd/4Qtf0JQpU1RdXa29e/fqm9/8pvbv36/f/OY3YcdZt26dvvOd76SrbAAAYCKH1+v1ml2EJK1atUp/+MMf9Kc//UmTJk2KON8rr7yia665Rs3NzZoxY8aw5/v6+tTX1+f/2e12q6amRi6XSyUlJSmpHQAAJJfb7ZbT6Yz5/m2JMzJ33HGHnn/+eTU2NkZtYiRp/vz5khSxkcnPz1d+fn5K6gQAANZiaiPj9Xq1Zs0abd26Vdu3b9e0adNiLtPU1CRJmjhxYoqrAwAAVmdqI7N69Wpt3rxZv/3tb1VcXKy2tjZJktPpVEFBgQ4ePKjNmzfr+uuvV3l5ufbu3au7775b9fX1mjVrlpmlAwAACzD1GhmHwxF2+saNG3Xbbbfp6NGj+tKXvqS3335bHo9HNTU1+sxnPqNvfetbhq93MfoZGwAAsA5bXCMTq4eqqalRQ0NDmqoBAAB2Y6kcGQAAgHjQyAAAANuikQEAALZFIwMAAGyLRgYAANgWjQwAALAtGhkAAGBbNDIAAMC2aGQAAIBt0cgAAADbopEBAAC2RSMDAABsi0YGAADYFo0MAACwLRoZAABgWzQyAADAtmhkAACAbdHIAAAA28o1uwCY41B7j9472aup5WM1rWKs2eUAAJAQGpks09Xbrzu3NKnxQLt/Wn1tpdYvnyNnYZ6JlQEAED8+Wsoyd25p0o7mjqBpO5o7tGbLHpMqAgAgcTQyWeRQe48aD7Rr0OsNmj7o9arxQLtaOjwmVQYAQGJoZLLIeyd7oz5/uJNGBgBgLzQyWWRKWWHU56eWc9EvAMBeaGSyyPTKItXXVirH4QianuNwqL62kruXAAC2QyOTZdYvn6MFMyuCpi2YWaH1y+eYVBEAAInj9uss4yzM0xMr56mlw6PDnR5yZAAAtsYZGZs71N6jV/efiHrHUbh5plWM1ShJz/31mP4YkCkDAICdcEbGpowE20Wa594l5+uWn+/Wqd4B//RxhXl6bvVC1ZRHvyAYAAArcXi9IaEiGcbtdsvpdMrlcqmkpMTscpLm1sd3a0dzR1AmTI7DoQUzK/TEynlR53E4pLNDw1/2cYV52vPta1NfPAAAMRh9/+ajJRsyEmwXbZ5wTYwkneod4GMmAICt0MjYkJFgu1jzRPLmkVMJLQcAgBm4RsaGjATbJfqJ4eWTxyW0HAAAZuCMjA0ZCbaLNk/uqOBpPuMK8/TJ2sqU1Q0AQLLRyNiUkWC7SPM8t3qBxn14Z5OP764lAADshLuWbM5IsF2kef54oF1vHjmlyyeP40wMAMBSjL5/08gAAADLscXt1+vWrdMVV1yh4uJijR8/XsuWLdP+/fuD5jlz5oxWr16t8vJyFRUV6aabbtLx48dNqjh9jCT2AgCQ7UxtZBoaGrR69Wq99tpreumllzQwMKBrr71WHs9Hb9533323fve73+mZZ55RQ0ODWltbdeONN5pYdWp19fbr1sd36+r/bNCXN76uRT/crlsf3y1XQAovAAA4x1IfLbW3t2v8+PFqaGhQfX29XC6XKisrtXnzZn32s5+VJO3bt08XXXSRdu7cqU984hMxx7TbR0tGEnsBAMh0tvhoKZTL5ZIklZWVSZLeeOMNDQwMaPHixf55LrzwQk2ePFk7d+4MO0ZfX5/cbnfQwy6MJPYCAICPWKaRGRoa0l133aUFCxbo0ksvlSS1tbVp9OjRKi0tDZp3woQJamtrCzvOunXr5HQ6/Y+amppUl540RhJ7AQDARyzTyKxevVpvv/22nnrqqRGNs3btWrlcLv/j6NGjSaow9Ywk9gIAgI9YopG544479Pzzz+vVV1/VpEmT/NOrqqrU39+vrq6uoPmPHz+uqqqqsGPl5+erpKQk6GEXRhJ7AQDAR0xtZLxer+644w5t3bpVr7zyiqZNmxb0/Ny5c5WXl6dt27b5p+3fv19HjhxRXV1dustNCyOJvQAA4BxT71r62te+ps2bN+u3v/2tLrjgAv90p9OpgoICSdKqVav0+9//Xps2bVJJSYnWrFkjSfrzn/9saB12u2vJx0hiLwAAmcoWyb4OR/gvL9y4caNuu+02SecC8b7+9a9ry5Yt6uvr05IlS/TTn/404kdLoezayAAAkM1s0cikA40MAAD2Y8scGQAAgHjQyAAAANuikQEAALZFIwMAAGyLRgYAANgWjQwAALAtGhkAAGBbNDIAAMC2aGQAAIBt0cgAAADbopEBAAC2RSMDAABsi0YGAADYFo0MAACwLRoZAABgW3E3MitWrFBjY2MqagEAAIhL3I2My+XS4sWLVVtbq//4j//QsWPHUlEXAABATHE3Ms8++6yOHTumVatW6emnn9bUqVO1dOlS/epXv9LAwEAqagQAAAgroWtkKisrdc899+ivf/2rdu3apZkzZ+qWW25RdXW17r77bh04cCDZdQK2dKi9R6/uP6GWDo/ZpQBARsodycIffPCBXnrpJb300kvKycnR9ddfr7feeksXX3yxHnroId19993JqhOwla7eft25pUmNB9r90+prK7V++Rw5C/NMrAwAMkvcZ2QGBgb061//Wv/wD/+gKVOm6JlnntFdd92l1tZW/eIXv9DLL7+sX/7yl/rud7+binoBW7hzS5N2NHcETdvR3KE1W/aYVBEAZKa4z8hMnDhRQ0NDWr58uXbv3q3LLrts2DyLFi1SaWlpEsoD7OdQe0/QmRifQa9XjQfa1dLh0bSKsSZUBgCZJ+5G5uGHH9bNN9+sMWPGRJyntLRULS0tIyoMsKv3TvZGff5wJ40MACRL3I3MLbfckoo6gIwxpaww6vNTy2liACBZSPYFkmx6ZZHqayuV43AETc9xOFRfW8nZGABIIhoZIAXWL5+jBTMrgqYtmFmh9cvnmFQRAGSmEd1+DSA8Z2Genlg5Ty0dHh3u9Ghq+VjOxABACtDIAAk41N6j9072xmxQplXQwABAKtHIAHEg6A4ArIVrZIA4EHQHANZCIwMY5Au6G/R6g6YHBt0BANKLRgYwyEjQHQAgvWhkAIMIugMA66GRAQwi6A4ArIdGBogDQXcAYC3cfg3EgaA7ALAWGhkgAQTdAYA1mPrRUmNjo2644QZVV1fL4XDo2WefDXr+tttuk8PhCHpcd9115hQLJNGh9h69uv+EoVu245kXALKNqWdkPB6PZs+era985Su68cYbw85z3XXXaePGjf6f8/Pz01UekHTxJAOTIgwAsZnayCxdulRLly6NOk9+fr6qqqrSVBGQWtGSgZ9YOS/heQEgW1n+rqXt27dr/PjxuuCCC7Rq1Sp1dnZGnb+vr09utzvoAVhBPMnApAgDgDGWbmSuu+46PfHEE9q2bZu+//3vq6GhQUuXLtXg4GDEZdatWyen0+l/1NTUpLFiILJ4koFJEQYAYyx919LnP/95/78/9rGPadasWZoxY4a2b9+ua665Juwya9eu1T333OP/2e1208zAEuJJBiZFGACMsfQZmVDTp09XRUWFmpubI86Tn5+vkpKSoAdgBfEkA5MiDADG2KqRef/999XZ2amJEyeaXQqQkHiSgUkRBoDYTP1oqaenJ+jsSktLi5qamlRWVqaysjJ95zvf0U033aSqqiodPHhQ9957r2bOnKklS5aYWDWQuHiSgUkRBoDYHF5vyG0RabR9+3YtWrRo2PQVK1Zow4YNWrZsmfbs2aOuri5VV1fr2muv1b//+79rwoQJhtfhdrvldDrlcrn4mAkAAJsw+v5taiOTDjQyqXeovUfvnezljAEAIGmMvn9b+q4lWBvJswAAs9nqYl9YS7TkWQAA0oFGBgkheRYAYAU0MkgIybMAACugkUFCSJ4FAFgBjQwSQvIsAMAKaGSQMJJnAQBm4/ZrJIzkWQCA2WhkslQyQ+ymVdDAAADMQSOTZQixAwBkEq6RyTKE2AEAMgmNTBYhxA4AkGloZLIIIXYAgExDI5NFCLEDAGQaGpksQogdACDT0MhkGULsAACZhNuvswwhdgCATEIjk6UIsQMAZAI+WgIAALZFIwMAAGyLRgYAANgWjQwAALAtGhkAAGBbNDIAAMC2aGQAAIBt0cgAAADbopEBAAC2RbIvbOdQe4/eO9k7oq9XSMYYqWDVugDAqmhkYBtdvf26c0uTGg+0+6fV11Zq/fI5chbmpW2MVLBqXQBgdXy0BNu4c0uTdjR3BE3b0dyhNVv2pHWMVLBqXQBgdTQysIVD7T1qPNCuQa83aPqg16vGA+1q6fCkZYxUsGpdAGAHNDKwhfdO9kZ9/nBn7Df7ZIyRClatCwDsgEYGtjClrDDq81PLY18Ym4wxUsGqdQGAHdDIwBamVxapvrZSOQ5H0PQch0P1tZWG7vBJxhipYNW6AMAOaGRgG+uXz9GCmRVB0xbMrND65XPSOkYqWLUuALA6h9cbcoVhhnG73XI6nXK5XCopKTG7HCRBS4dHhzs9I8paScYYqWDVugAg3Yy+f9PIAAAAyzH6/m3qR0uNjY264YYbVF1dLYfDoWeffTboea/Xq29/+9uaOHGiCgoKtHjxYh04cMCcYpERDrX36NX9J7ilGQAyhKmNjMfj0ezZs/Xoo4+Gff6hhx7Sj3/8Yz322GPatWuXxo4dqyVLlujMmTNprhR219Xbr1sf362r/7NBX974uhb9cLtufXy3XL0DZpcGABgBy3y05HA4tHXrVi1btkzSubMx1dXV+vrXv65vfOMbkiSXy6UJEyZo06ZN+vznP29oXD5agiTd+vhu7WjuCAqdy3E4tGBmhZ5YOc/EygAA4djio6VoWlpa1NbWpsWLF/unOZ1OzZ8/Xzt37oy4XF9fn9xud9AD2Y3kXADIXJZtZNra2iRJEyZMCJo+YcIE/3PhrFu3Tk6n0/+oqalJaZ2wPpJzASBzWbaRSdTatWvlcrn8j6NHj5pdEkxGci4AZC7LNjJVVVWSpOPHjwdNP378uP+5cPLz81VSUhL0QHYjORcAMpdlG5lp06apqqpK27Zt809zu93atWuX6urqTKwMdkRyLgBkplwzV97T06Pm5mb/zy0tLWpqalJZWZkmT56su+66S9/73vdUW1uradOm6V//9V9VXV3tv7MJMMpZmKcnVs4jORcAMoypjcxf/vIXLVq0yP/zPffcI0lasWKFNm3apHvvvVcej0f/+I//qK6uLi1cuFAvvPCCxowZY1bJsLlpFTQwAJBJLJMjkyrZliNzqL1H753s5YwDAMDWjL5/m3pGBsnT1duvO7c0qfFAu39afW2l1i+fI2dhnomVAQCQOpa92BfxuXNLk3Y0dwRN29HcoTVb9phUEQAAqUcjkwFIrgUAZCsamQxAci0AIFvRyGQAkmsBANmKRiYDkFwLAMhWNDIZguRaAEA24vbrDEFyLQAgG9HIZJhMSa4NDPbzer2E/AEAwqKRgaWEC/YLRMgfACAQ18jAUsIF+wUi5A8AEIhGBpYRKdgvECF/AIBANDKwjFjBfoEI+QMASDQysJBYwX6BCPkDAEg0MrCQSMF+gQj5AwAEopGBpYQL9gtEyB8AIBC3X8NSwgX7SSLkDwAQFo0MLCk02I8GBgAQDo0MLCcw1TebGphUb3c842frawDAfmhkYBnhUn2zIck31dsdz/jZ+hoAsC8u9oVlhEv1zYYk31RvdzzjZ+trAMC+aGRgCZFSfTM9yTfV2x3P+Nn6GgCwNxoZWEKsVN9MTfJN9XbHM362vgYA7I1GBpYQK9U3U5N8U73d8Yyfra8BAHujkYElREr1zfQk31RvdzzjZ+trAMDeaGRgGeFSfbMhyTfV2x3P+Nn6GgCwL4fXG3JlX4Zxu91yOp1yuVwqKSkxuxwYEJjqm01nAVK93fGMn62vAQDrMPr+TSODlLFCqJrRGgLn83q9I677UHuPdrV0SnLoE9PLh42TjH1jhf1rRaH75dxrcVIOSfPDvBYArMno+zeBeEg6K4SqGa0h3HyB4q27q7dfX3vyTf35YGfQ9Lrp5XrsS3PllXfE+8YK+9eKwu0XZ0GuXKfPBs135Yxybfji3KzeV0Am4YwMku7Wx3drR3NHUB5JjsOhBTMr9MTKeZaqIdx8geKt+9bHd0dtiiSNeN9YYf9aUazXMlB9bWVW7yvADoy+f3OxL5LKCqFqRmuINF+0ZYysN5LGA+0j3jdW2L9WZOS1DJTN+wrINDQySCorhKoZrSHWfOGWGcl607GObA2tS2TfZ+u+AjINjQySygqhakZriDVfuGVGst50rCNbQ+sS2ffZuq+ATEMjg6SyQqia0RoizRdtGSPrjaS+tnLE+8YK+9eKjLyWgbJ5XwGZhkYGSWeFUDWjNYSbL9YysdZ75YzyYdPrppdr/fI5Sdk3Vti/VhRuvzgLht+YeeWM8qzfV0Am4a4lpIwVQtWM1hA4n6QR193S4dGuQ53ySmFzZJKxb6ywf60odL+0dHj02qFOcmQAmyEQ70M0MgAA2E9G3H79b//2b3I4HEGPCy+80OyykAKH2nv06v4TSbslNtnjhY69ZfcRPbX7SFLGj1ZrvNvhm//p14/oR9v+rj9GuR083PiR1pfK/Wm0NiuOb3QMM/YfkC0sn+x7ySWX6OWXX/b/nJtr+ZIRh2Sn1KYy9bart1+r/udN7TwUnNqbaFJstFrjTQCOllA8rjBPz61eqJrywqjzjyvM06negaD1fW/ZpfrWs2+nNUU41cnFyRh/JMnRpDADyWXpMzLSucalqqrK/6ioiHxhJuznzi1N2tHcETRtR3OH1mzZY4nxQscObWIk6c8HOxMaP1qt8W5HuPl9TvUO6P979E8x5w9sYnzr+/Sjf0rZ/owkla9hssY3OkaqtwWADRqZAwcOqLq6WtOnT9cXv/hFHTlyJOr8fX19crvdQQ9YU7JTalOZemsktTee8WPVGs92GEm1PdU74P+YyWgK7qDXq1O9A2lNEU51cnEyxh9pcnS2pzADyWbpRmb+/PnatGmTXnjhBW3YsEEtLS365Cc/qe7u7ojLrFu3Tk6n0/+oqalJY8WIR7JTalOZemskOTae8UeSAhy6HqNjvXnk1IjXHa2OZEh1cnEyxk9WcjTJwkByWLqRWbp0qW6++WbNmjVLS5Ys0e9//3t1dXXpl7/8ZcRl1q5dK5fL5X8cPXo0jRUjHslOqU1l6q2R5Nh4xh9JCnDoeoyOdfnkcSNed7Q6kiHVycXJGD9ZydEkCwPJYelGJlRpaanOP/98NTc3R5wnPz9fJSUlQQ9YU7JTalOZemsktTee8WPVGs92GEm1HVeYp09+WL/RFNwch0PjCvPSmiKc6uTiZIw/0uTobE9hBpLNVo1MT0+PDh48qIkTJ5pdCpIk2Sm1qUy9Xb98juqmD0/tTTQpNlqt8W5HtIRi311LseYfF3IXzYKZFXpu9cK0pwinOrk4nenKpDADqWfpQLxvfOMbuuGGGzRlyhS1trbq/vvvV1NTk/72t7+psjLy/44DEYhnD8lOqU1l6m2yk2Kj1Rrvdvjm7+juU6vrtC6fPM5/JsbIuiOtz4wU4VSvM53pyqQwA/HLiGTfz3/+82psbFRnZ6cqKyu1cOFCPfDAA5oxY4bhMWhkAACwH6Pv35ZOl3vqqafMLgFIqkPtPXrvZG/U/5lHmudQe492tZyUQ1J1aYEGvV7/PIHP+c4SnZvWKckR9vuejNSSSP3RlvF6vRG3Ld5xjY49kvFHyqz1Wq0GIJUs3cgAmcJIwmukeb637BJ989dvhQ3jkyRnQZ5cp4PD7ErG5Mp95mzQtLrp5XrsS3PjTg02Wr+RZQL5tu1bz74Td/KtkbETSUhOFisk+lqhBiAdLP3RUjLw0RKs4NbHd2tHc0dQOFqOw6EFMyv0xMp5UecpKcgdlrqbKN+dV7FqSaR+I8sE8m2b+/TZuMY1OrbvItt4606GRPZXJtYAjERGfGkkkAmMJLxGmydZTYykuFODjdZvdJnQ5RNJDjY6diLbmgxWSPS1Qg1AutDIAClmJOE1WWm7IxUubTaRhNpkbE+k5FsrJxNL1kj0tUINQLpwjQyQYkYSXq3yCW+4tNlEEmqTkR4cKfnWysnEkjUSfa1QA5AunJEBUsxIwmu0eUKD6kYi3tRgo/UbXSZ0+USSg42Onci2JoMVEn2tUAOQLjQyQBoYSXiNNM9zqxeGTRT2cRYMb3RKxgw/2Vo3vTyh1GCj9RtZJnT5RJODjYyd6LYmgxUSfa1QA5AO3LUEpJGRhNdo6bq+ROHzxhXo7JA3KJE3NG24pcOjXYc65ZXC5sgkkjY70mUkJTU52MjYIxl/pKyQ6GuFGoBEZEQgHhCOL+Arx+EICoUbyVjJ/CMfWl+OQxr0yr8Or9frv9jSt87QOiKFxS2fNzlomo/X61WVc0zQdoSOE6hh/wk1vd+lyyePixlMFzjvJ8N8LOF7fmJJgSpL8ofty8D/K4X+vykwtO+80jF67VCnnv9rq056+uQ6fVbTK8fq/86qjrqffOF/0V7HaPsiUUaOHaOvdyqlYtsjCbddBPIh1WhkYBvRQtDiDfpKRVhYrJA26dwXMwbeTl03vVwOh/Tngx+F3UULc7tyRrm8XgWF44WOGW073uv0aNmjO4LmH1eYp//5ynx9/8X9QeuaN3Wc9h/vluv02aB5n1u9UDXlhWHHCqwhXNhdoNKCPHWdHr5sqB/+7991eY1TY0bnDttPiQbqjYTRYyfcfNFebzuH1IXb1nDHaiZsK6yHj5ZgG9FC0OIN+kpFWFiskDajooW5xbN8uO2Y893/Ddt45I5yyOuVoXWNK8zTnm9fG3EsXw3hwu6SaSSBeiNh9NgxejxkQkhdNm0r0odAPGSUWCFo8QR9pSIszEhIm1HRwtziWT50Oxr2n4jYeJwd8hpe16neAT36SnPUoL5IYXfJlGig3kgYPXbiOR7sHlKXTdsKa6KRgS0YDUEzEvSVirAwqwTaBQrdjqb3u5I29p+aI398ZhWpCH0zeuwkcjzYNaQum7YV1kQjA1swGoJmJOgrFWFhyQppS6bQ7bhsUmnSxl44szJpY6VKKkLfjB47iRwPdg2py6ZthTXRyMAWYoWgxRP0lYqwMCMhbUZFC3OLZ/nQ7fg/F4yPGK6XO8pheF3jCvO0+uqZUYP6IoXdJVOigXojYfTYied4sHtIXTZtK6yJRga2ES0ELd6gr1SEhcUKaZM07M2/bnq5rpwRHHYXLcztyhnlw8LxQseMth3PrV44bP5zdyItGLaueVPL5CzIDTPvwohjBdYQLuwuUGmYIL9ILq9xht1PiQbqjYTRYyfcfNFebzszeqxmwrbCerhrCbbjC/jKHeUICoUbyVjJzLgIrS+0znDrjFaHkfnj3Y4/HmjXm0dO+bNhoq0r0ryhY1U7C1RRPDxHJlpoXWBo36RxBXr/1Gl19vTplKdfp3oHhuXIRAsLTHfom9F1xvt621k2bStSz+j7N40MAACwHJJ9MSLxpnEamT9wHq/Xa4m0z49SeIPTd8PN43suGdthZL2pEpim+4np5fJ6vdrV0qn27n6NL87X/A+nRUpoDVfzuTFPqqO7T5XF+aouHRN2u0JTggPHDbd/I/0cWG91aUHYhGdfTYFf2xBpf4R7DcPVGm7/pTtVOtL+CKwnnrGTnbybir8d6RgD9kUjgyDxJt4amT9W4q0ZaZ9GUoLDJeuGpuiGWy7adiQznTheXb39+tqTbwalyhoRLo028LmzQ0N6/fCpiMvX11bq3iXn65af7w7ad6UFeaodX6TX3/to2dD9G/pzrDRgX9rvN3/9VlCirHTumo0NX5wb9bisr63Uvdedr1se3z08/XjlfD3w+3eH7Ye66eV67EtzU54qHW65cPvDWZAnV8C0SGMnO906FX87kr1OZCY+WkKQeBNvjcwfK/XTjLRPIynBUnzJuka2I5npxPG69fHdUb8+IVVyHA45HOdC99KxrpKC3KjNZrTjMlqtvmudYo0bS6Kp0okmR0caO9np1qn425HsdcJeSPZF3OJNvDUyv5HUz3SnfRpNCY43WTfWdiQznThevnWbYdDrTUsT41tXtMThWMdltFqjbUOqU6VHkhwdbuxkp1un4m9HsteJzEUjA794E2+NzB9P6me60j5TncIbaTuSmU4cLysmD5sl3uMynnFjSTRVOhn1Bo6d7HTrVPztSPY6kbm4RgZ+8SbeGpk/nk8u05X2meoU3kjbkcx04nhZMXnYLPEel/GMG0uiqdLJeP0Cx052unUq/nYke53IXJyRgV+8ibdG5jeS+pnutE+jKcHxJuvG2o5kphPHy7duM+Q4HModlbqE39B1RUscjnVcRqs12jakOlV6JMnR4cZOdrp1Kv52JHudyFw0MggSb+KtkfljJd6akfZpJCU43DzR3iSNbEcy04njtX75nGGpskaES6MNfG7e1HFRlz+XwLtg2L4rLcjTFVOClw2XOhy6TOx1LRyWKCudu2sp1nEZqVZf+nG4/VA3vTyu1y3RVOlwy4XbH86QaZHGTna6dSr+diR7nchM3LWEsOJN4zQyf7SEV7MYSQmOlqIrJbYdyUwnjldgmu4nPnzD33WoU+09faosOpcjIw3frmg1t3R49NqhTnX29KmiKF/njSsIu13hUoJjpRSH+zmw3kjr8tUUK0cm0rEbLf04cP+lO1U60v4IrCeesZOdvJuKvx3pGAPWQ7Lvh1LVyNg5gCkdtcdaR7r3n9FwtHjHDNyGwAC180oLIgbK+a7ReH5vq056+nXxRKcqS/KDguZCA+l8oWcdPX2qKBoz7A3UyPZF2udGAv6MbuuRTo+ebWqVQ9JnLj/P/9wJ9xm9+4FbDjkkh+T1enXNRRP8DcLTu49oZ0unzh9frHFj8xQc1ncuaE+STnr65JVUWZyv9u4+dZ8+K6+8Ki0cLa/O5YrMqCwK+lqD0H3j28YchyNskF7oPnt+b6sOd/T6vy7Bt/wJ9xm1uc8o1+FQ+4cNXODzRsIVw60v3OscKYTP9zpMLDkXQGjk+A59vX37JlKIoZFjxOjvs9FxIr02vtfjpKc/6PhBdHZ9v6KR+VCyGxk7BzClo/ZY60j3/uvq7deq/3kzZjhavGOGbkO0jJFogXJGRAqBq5teru/fNEvf/PXeqNsXaZ9/b9kl+taz70QNKgw3T7RtjVfR6FE6fXZIg0NJGS7I5ZNLleNwBAXuRRIuxPH2J/4SNejPiGjhiqG/F5HCCkvG5Mp95mzQtMsmOXW4szdiOGC44ztWMKXRugOFO7bD/T7HWveVM8rl9WrYcewb73vLLtE3nvmrdoe8Hs6CXD1/xydVU87F7OHY+f1KopHxS3YjY+cApnTUHmsd6d5/0ULg4gkxCx0zkWCyVIiVNBxtn5cU5Mp9+mzUoMJY82SKcCGOycjdiRauGPp7keycn9DjO57jNtFQyMBlE113uPGihRyOK8zTnm9fG/e42cDO71cSgXgpYecApnTUHmsdjX8PHzKXqv0XKwQumcFfZokV/tb49xMR9/mp3oGYQYWx5skU4UIckzlu9N+LEykJKww8vuM9bhMNhQxcNtF1hxsv2nF+qndAfzQp7NHK7Px+FS8amTjYOYApHbXHWseeo9FP0yd7/xkJEUt28JfV7DnaZXYJtpKqsLxoUvka+Y5vM47bdK77zSMj+wgwE9n5/SpeBOLFwc4BTOmoPdY65tREv0032fvPSIhYsoO/rGZOTanZJdhKqsLyoknla+Q7vs04btO57ssnR//bko3s/H4VL87IxMHOAUzpqD3WOurPDx8yl6r9FysELpnBX2aJFf5Wf/74iPt8XGFezKDCWPNkinAhjskcN/rvxfiUhBUGHt/xHreJhkIGLpvousONF+04H1eYx91LYdj5/SpeNDJxsnMAUzpqj7WOdO+/9cvnGApHi3fM0G2IlvoaLVDOiEghcHXTyw2Fv0UOflsYM6gw3DzJTOktGj1KOSn6K3T55NJhgXuRhAtxjBX0F8+4Rn4vIh0jJWOGnzifM8kZNRww3PEdK5jSaN2Bwh3b4X6fY41z5YzysMexb7znVi/UvKllw55zFuTqudULo21KVrPz+1U8uGspQXYOYEpH7bHWke79ZzQcLd4xA7chMEBt0rjCiIFyvlO6/29vqzp6+nRJtVMVxflBQXOSwoaedfT0qbwof1iOjJHti7TPjQT8Gd3W90/1auubxySdy5HxPdfR3ae/tbrlOBcjo8GQHJln/nJUOw526PzxxSobOzoorM8XtCdJpzz9GvR6Nb54jE50n1HPmbMa8no1rnD0ued7+zU9JEcmdN/4ttFIGGFLh0f/b2+rDrV7/DkyvuU7uvvU6jqtvFGjdKL7jD9HJtI+jPYaBD4f7nWOFMLnex2qnedCAY0c36Gvt2/fRAoWNHKMGP19NjpOpNfG93p09PSRIxMHu75fcfv1h0j2BQDAfoy+f9viYt9HH31UP/jBD9TW1qbZs2dr/fr1mjfP+vfAw5qipdSG/g/PSBpuaIJsR/cZVRSPCUqmjfU/ZSPJm77x/97mllcOXTyxRJUl+f6E2XD1D1/+o3TYwGTVI52egIRYr9q7+3XK0y9JumhisX+aQ+dSdR0O6d0P3GrtOq2evrO6sKpEt9RN9W/v39u6tb/NLTmkz8yZpJs/XjNs/3R292nHwQ4V5+dp8cUTNOj1fpTg29snV+9ZTa8cq8rifLW5z2hiSYEqS/IDEmk79XrLKXWfGdCUskIVF+bpnWMuNR/vUe2EYt13/UUfpfC6zujdtm65Tw+opCBPF1cXq7J4zIfT3TpwvFue/kHVji/S2SGvPGfO6mM1pZo9yammo11q6fBolMOh4jG5uri6RO+2unX01GlNKS9USUGeP+H4XOrsgHId0iv7Tsjdd1aTxxVqasVY/xkdX1Jv6OvtSzdeMKNCQ0NevfTucTkL8uT1etXm6lN16Rh9fGqZ/zjypdy2dPTKc2ZAQ5KmlI3Vux+4dLCjR+WF+bpyRrnOer0a5XBo3NjR524vP9Gjsfm5umBCkbxy6P1TvSoek6fPXH6ePllbqad3H9GzTcd0wt2n0bkOXXpeqc4ODckhhz5z+XkaGvL6f19aT53WzpZO5Tkcer/rtPr6h1RePFqza0r929qw/4Re3X9CDodD8kpeSRcHHFPji/NVXTpGx7rO+BOGA/8dmjYc+vvq25fvHHOpub1HC2ZU6OaP13y43nZVFI3WrEnOoHF8NY1yOHT+hGJ/4vOulk5/TYG/r7GSs8OlJfuO8/NKC9R09NSwNOJE0nYjbbsvCTkwHTxZZ16sltIejuXPyDz99NO69dZb9dhjj2n+/Pl65JFH9Mwzz2j//v0aP358zOU5IwMfI4m8574ccKGKC3INp+HWTS/X2aEhQymwoYmrRpI3I6UTh+OrPzDpNFpqLNIvNKl3To1Tbx1zx5WWXJKfK3ff2dgzmmyUQ0pSCLQ+PmWcmk/0BKUZJ5IyHc8y86aOk9erYenQvt/jWMnH4ZSMydUFE4qDxoyVtvtep0fLHt0RlKcTaztGmuBrhZT2jPloaf78+briiiv0k5/8RJI0NDSkmpoarVmzRvfdd1/M5Wlk4GM0XXRcYZ4+dl5pQmm4RgQmrhpJ3ow39TU06TQVqbFAtvPdaZaMlO9Yabtzvvu/UUMBExkzFiuktGdEsm9/f7/eeOMNLV682D9t1KhRWrx4sXbu3Bl2mb6+Prnd7qAHEE+66KnegYTTcI0ITZCNlryZSMpsYNJpMlNqAXwk0eTjcKKl7TbsPxF3ExNrzFisltIei6UbmY6ODg0ODmrChAlB0ydMmKC2trawy6xbt05Op9P/qKmpSUepsDirJfIaSZAdScqsL+nUatsNILJwabtN73clfcxYrJbSHoulG5lErF27Vi6Xy/84evSo2SXBAqyWyDu1fKyh5M1E6/YlnVptuwFEFi5t97JJpUkfMxarpbTHYulGpqKiQjk5OTp+/HjQ9OPHj6uqqirsMvn5+SopKQl6APGki44rzEs4DdeI0ATZaMmbiaTMBiadJjOlFsBHEk0+Dida2u7/uWB81GTjRMaMxWop7bFYupEZPXq05s6dq23btvmnDQ0Nadu2baqrqzOxMtiRkURe310/8aTh1k0vN5wCG5q4aiR5M1I6cTi++gNFS41F+oUm9V5e44w7Lbkk3xbJGUpiCLSumDJuWJpxIinT8Swzb2pZ2HRo3+9xPGnJPiVjcoeNGStt97nVC4c1M7G2Y6QJvlZLaY/G8nctPf3001qxYoV+9rOfad68eXrkkUf0y1/+Uvv27Rt27Uw43LWEUNFSakNzWIyk4YYmyHYGpLJKMpQobCR50zd+8/FuDXqlS6pLVFGc70+YjZUjE5oOKykokdeXEDs45FV7T5+6PP0a8koXV5f4pzkkVRTlK2eUQ++0uvRB1xl19w34c2R829t8vFv72rrlldefIxO6f055+vWnAx0qys/Vpy6ZoLND3qAE31O9A5peOVYTSsao1XVa1c4CVRTn+09b7zrUqdcPn5T79ICmlI9VcUGu3m116+/Hu1U7/lyOjG8bO7r79E6rWz1nBlQ0Ji9o373Tei57pqf/rGrHF2twyKueM2f1sRqnLqsp1V+PdulQe49GORwqGpOrS6qd2veB25+dUVyQ6084Ppc626/cUQ69uu+EXGcGNGVcoaaE5MiEe7196cYLZpx7c3jxnTb/G3dr15lhOTKBqcOevnN30k0tH6t3P3DrYHu3ygrzdeXMCp0dGlKOw6HSsaP1XodHB9t7VDg6VxdWFcsr6f2Tp1U0JtefI/PMX45q65vv63h3n0bnnMuRGfzwNt/PXH6eJPl/X9pcZ7TjYIfyRo3Ssa5enfHlyEz6KEfmjwfa9cq7JwJSnc8du75jqrIoX+eNK9D7p077E4YD/x2a8Bv6++rbl++2uvX3E93+HJk/HmjXtndPqKLoXK5N4DiBNflyZHzHlK+mwN/XWMnZ4dKSfcf5pHEF+uvRrmFpxImk7Ubadt+t2EbSquNlZkp7xtx+LUk/+clP/IF4l112mX784x9r/vz5hpalkQEAwH4yqpEZCRoZAADsJyNyZAAAAKKhkQEAALZFIwMAAGyLRgYAANgWjQwAALAtGhkAAGBbNDIAAMC2aGQAAIBt0cgAAADbssc3j42AL7jY7XabXAkAADDK974d6wsIMr6R6e7uliTV1NSYXAkAAIhXd3e3nE5nxOcz/ruWhoaG1NraquLiYjkcSfxO+SzkdrtVU1Ojo0eP8r1VScD+TD72afKxT5OL/Wmc1+tVd3e3qqurNWpU5CthMv6MzKhRozRp0iSzy8goJSUl/AImEfsz+dinycc+TS72pzHRzsT4cLEvAACwLRoZAABgWzQyMCw/P1/333+/8vPzzS4lI7A/k499mnzs0+RifyZfxl/sCwAAMhdnZAAAgG3RyAAAANuikQEAALZFIwMAAGyLRgZxO3z4sFauXKlp06apoKBAM2bM0P3336/+/n6zS7O1Bx54QFdeeaUKCwtVWlpqdjm29Oijj2rq1KkaM2aM5s+fr927d5tdkm01NjbqhhtuUHV1tRwOh5599lmzS7K1devW6YorrlBxcbHGjx+vZcuWaf/+/WaXlRFoZBC3ffv2aWhoSD/72c/0zjvv6OGHH9Zjjz2mf/7nfza7NFvr7+/XzTffrFWrVpldii09/fTTuueee3T//ffrzTff1OzZs7VkyRKdOHHC7NJsyePxaPbs2Xr00UfNLiUjNDQ0aPXq1Xrttdf00ksvaWBgQNdee608Ho/Zpdket18jKX7wgx9ow4YNOnTokNml2N6mTZt01113qaury+xSbGX+/Pm64oor9JOf/ETSue9Zq6mp0Zo1a3TfffeZXJ29ORwObd26VcuWLTO7lIzR3t6u8ePHq6GhQfX19WaXY2uckUFSuFwulZWVmV0GslR/f7/eeOMNLV682D9t1KhRWrx4sXbu3GliZUB4LpdLkvi7mQQ0Mhix5uZmrV+/Xv/0T/9kdinIUh0dHRocHNSECROCpk+YMEFtbW0mVQWENzQ0pLvuuksLFizQpZdeanY5tkcjA7/77rtPDocj6mPfvn1Byxw7dkzXXXedbr75Zt1+++0mVW5diexTAJlt9erVevvtt/XUU0+ZXUpGyDW7AFjH17/+dd12221R55k+fbr/362trVq0aJGuvPJK/dd//VeKq7OnePcpElNRUaGcnBwdP348aPrx48dVVVVlUlXAcHfccYeef/55NTY2atKkSWaXkxFoZOBXWVmpyspKQ/MeO3ZMixYt0ty5c7Vx40aNGsXJvXDi2adI3OjRozV37lxt27bNf0Hq0NCQtm3bpjvuuMPc4gBJXq9Xa9as0datW7V9+3ZNmzbN7JIyBo0M4nbs2DFdddVVmjJlin74wx+qvb3d/xz/+03ckSNHdPLkSR05ckSDg4NqamqSJM2cOVNFRUXmFmcD99xzj1asWKGPf/zjmjdvnh555BF5PB59+ctfNrs0W+rp6VFzc7P/55aWFjU1NamsrEyTJ082sTJ7Wr16tTZv3qzf/va3Ki4u9l+75XQ6VVBQYHJ1NucF4rRx40avpLAPJG7FihVh9+mrr75qdmm2sX79eu/kyZO9o0eP9s6bN8/72muvmV2Sbb366qthj8cVK1aYXZotRfqbuXHjRrNLsz1yZAAAgG1xYQMAALAtGhkAAGBbNDIAAMC2aGQAAIBt0cgAAADbopEBAAC2RSMDAABsi0YGAADYFo0MAACwLRoZAJY1ODioK6+8UjfeeGPQdJfLpZqaGv3Lv/yLJOnOO+/U3LlzlZ+fr8suu8yESgGYhUYGgGXl5ORo06ZNeuGFF/Tkk0/6p69Zs0ZlZWW6//77/dO+8pWv6HOf+5wZZQIwEd9+DcDSzj//fD344INas2aNrr76au3evVtPPfWUXn/9dY0ePVqS9OMf/1iS1N7err1795pZLoA0o5EBYHlr1qzR1q1bdcstt+itt97St7/9bc2ePdvssgBYAI0MAMtzOBzasGGDLrroIn3sYx/TfffdZ3ZJACyCa2QA2MLPf/5zFRYWqqWlRe+//77Z5QCwCBoZAJb35z//WQ8//LCef/55zZs3TytXrpTX6zW7LAAWQCMDwNJ6e3t12223adWqVVq0aJEef/xx7d69W4899pjZpQGwABoZAJa2du1aeb1ePfjgg5KkqVOn6oc//KHuvfdeHT58WJLU3NyspqYmtbW16fTp02pqalJTU5P6+/tNrBxAOji8nJ8FYFENDQ265pprtH37di1cuDDouSVLlujs2bN6+eWXtWjRIjU0NAxbvqWlRVOnTk1TtQDMQCMDAABsi4+WAACAbdHIAAAA26KRAQAAtkUjAwAAbItGBgAA2BaNDAAAsC0aGQAAYFs0MgAAwLZoZAAAgG3RyAAAANuikQEAALb1/wOgrlbYhQ6sXAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "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": 43,
   "id": "89018bd0-5f4d-46c7-b017-9de037f1a788",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Axes: xlabel='X1', ylabel='preds'>"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjIAAAGwCAYAAACzXI8XAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAPKVJREFUeJzt3X14VPWd///XyeSGhCSTuyEhEgiQeAcVAwrGuKmuVEq/2y2t2m62W7Xl6o0/hEvtXq20VqttF6/aq1qji+1qvbm2oNtuoa1ttUqBdCk3FUmRKpiQBIIQmEnI3Qzk9vz+wEwzYZKZhEnOnMnzcV1zyZzzOZ/Pe86cybw8c24M0zRNAQAA2FCc1QUAAACMFUEGAADYFkEGAADYFkEGAADYFkEGAADYFkEGAADYFkEGAADYVrzVBYy3/v5+HT9+XGlpaTIMw+pyAABAGEzTVEdHh/Lz8xUXN/x+l5gPMsePH1dBQYHVZQAAgDFobGzUjBkzhp0f80EmLS1N0rkVkZ6ebnE1AAAgHO3t7SooKPB/jw8n5oPMwM9J6enpBBkAAGwm1GEhHOwLAABsiyADAABsiyADAABsiyADAABsiyADAABsiyADAABsiyADAABsiyADAABsiyADAABsiyADAABsK+ZvUWBXde5OHWnxqTB7qmbnTL3gdgAAxCKCTJRp9XVrzcZqVdW4/dPKi12qrCiRMyVh1O0AAIhllv60tH79el1xxRX+GzqWlpbq97//vX/+9ddfL8MwAh5f+cpXLKx4/K3ZWK0dtZ6AaTtqPVq9cd+Y2gEAEMss3SMzY8YMPfLIIyouLpZpmnrhhRf0iU98Qvv27dO8efMkSV/84hf18MMP+5dJSUmxqtxxV+fuDNjDMqDPNFVV41a9x6vZOVPDbgcAQKyzNMh8/OMfD3j+ve99T+vXr9euXbv8QSYlJUV5eXlWlDfhjrT4Rpzf0HwuoITbDgCAWBc1Zy319fXppZdektfrVWlpqX/6z372M+Xk5Gj+/Plau3atfL6Rv8S7urrU3t4e8LCLWVkj720qzJ46qnYAAMQ6yw/2ffvtt1VaWqqzZ88qNTVVmzZt0uWXXy5J+td//VfNmjVL+fn52r9/v77+9a/r0KFD+uUvfzlsf+vWrdNDDz00UeVH1BxXqsqLXdpR61GfafqnOwxDZUU5/r0s4bYDACDWGaY56JvQAt3d3Tp69Kja2tr0i1/8Qs8884y2b9/uDzOD/fGPf9SNN96o2tpazZ07N2h/XV1d6urq8j9vb29XQUGB2tralJ6ePm6vI1LafD1avXFfyLORwm0HAIAdtbe3y+l0hvz+tjzIDLV06VLNnTtXP/7xj8+b5/V6lZqaqldffVXLli0Lq79wV0S0qfd41dDsDXl9mHDbAQBgJ+F+f1v+09JQ/f39AXtUBquurpYkTZ8+fQIrssbsnPCCSbjtAACIRZYGmbVr12r58uWaOXOmOjo6tGHDBm3btk2vvfaaDh8+rA0bNuhjH/uYsrOztX//ft1zzz0qLy/XFVdcYWXZAAAgSlgaZE6dOqXbbrtNJ06ckNPp1BVXXKHXXntNH/nIR9TY2Kg33nhDjz/+uLxerwoKCnTzzTfr/vvvt7JkAAAQRaLuGJlIs+sxMgAATGbhfn9HzXVkAAAARosgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbCve6gLsqs7dqSMtPhVmT9XsnKlWlzMmsfAaAACTG0FmlFp93VqzsVpVNW7/tPJilyorSuRMSbCwsvDFwmsAAEDip6VRW7OxWjtqPQHTdtR6tHrjPosqGr1YeA0AAEgEmVGpc3eqqsatPtMMmN5nmqqqcave47WosvDFwmsAAGAAQWYUjrT4Rpzf0Bz9ISAWXgMAAAMIMqMwKytlxPmF2dF/wGwsvAYAAAYQZEZhjitV5cUuOQwjYLrDMFRe7LLFmT+x8BoAABhAkBmlyooSlRXlBEwrK8pRZUWJRRWNXiy8BgAAJMkwzSFHfU6g9evXa/369WpoaJAkzZs3Tw888ICWL18uSTp79qy++tWv6qWXXlJXV5eWLVum//zP/1Rubm7YY7S3t8vpdKqtrU3p6ekRq73e41VDs9fW12CJhdcAAIhN4X5/WxpkfvOb38jhcKi4uFimaeqFF17Qo48+qn379mnevHm688479dvf/lbPP/+8nE6n7rrrLsXFxWnHjh1hjzFeQQYAAIwfWwSZYLKysvToo4/qlltukcvl0oYNG3TLLbdIkg4ePKjLLrtMO3fu1DXXXBNWfwQZAADsJ9zv76g5Rqavr08vvfSSvF6vSktLtXfvXvX09Gjp0qX+NpdeeqlmzpypnTt3DttPV1eX2tvbAx4AACA2WR5k3n77baWmpiopKUlf+cpXtGnTJl1++eVqampSYmKiMjIyAtrn5uaqqalp2P7WrVsnp9PpfxQUFIzzKwAAAFaxPMhccsklqq6u1u7du3XnnXfq9ttv1zvvvDPm/tauXau2tjb/o7GxMYLVAgCAaGL5TSMTExNVVFQkSVq0aJH+8pe/6Ec/+pE+85nPqLu7W62trQF7ZU6ePKm8vLxh+0tKSlJSUtJ4lw0AAKKA5Xtkhurv71dXV5cWLVqkhIQEbdmyxT/v0KFDOnr0qEpLSy2sEAAARAtL98isXbtWy5cv18yZM9XR0aENGzZo27Zteu211+R0OrVy5Urde++9ysrKUnp6ulavXq3S0tKwz1gCAACxzdIgc+rUKd122206ceKEnE6nrrjiCr322mv6yEc+Ikl67LHHFBcXp5tvvjnggngAAABSFF5HJtK4jgwAAPZju+vIAAAAjBZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2Fa81QUgsurcnTrS4lNh9lTNzplqdTkAAIwrgkyMaPV1a83GalXVuP3TyotdqqwokTMlwcLKAAAYP/y0FCPWbKzWjlpPwLQdtR6t3rjPoooAABh/BJkYUOfuVFWNW32mGTC9zzRVVeNWvcdrUWUAAIwvgkwMONLiG3F+QzNBBgAQmwgyMWBWVsqI8wuzOegXABCbCDIxYI4rVeXFLjkMI2C6wzBUXuzi7CUAQMwiyMSIyooSlRXlBEwrK8pRZUWJRRUBADD+OP06RjhTEvTiysWq93jV0OzlOjIAgEmBIBNjZucQYAAAkwc/LQEAANsiyAAAANsiyAAAANsiyAAAANuyNMisW7dOV199tdLS0jRt2jStWLFChw4dCmhz/fXXyzCMgMdXvvIViyoGAADRxNIgs337dq1atUq7du3S66+/rp6eHt10003yegMvqf/FL35RJ06c8D++//3vW1QxAACIJpaefv3qq68GPH/++ec1bdo07d27V+Xl5f7pKSkpysvLC6vPrq4udXV1+Z+3t7dHplgAABB1ouoYmba2NklSVlZWwPSf/exnysnJ0fz587V27Vr5fMPfJHHdunVyOp3+R0FBwbjWDAAArGOYpmlaXYQk9ff365//+Z/V2tqq//u///NP/8lPfqJZs2YpPz9f+/fv19e//nUtXrxYv/zlL4P2E2yPTEFBgdra2pSenj7urwMAAFy49vZ2OZ3OkN/fUXNl31WrVunAgQMBIUaSvvSlL/n//aEPfUjTp0/XjTfeqMOHD2vu3Lnn9ZOUlKSkpKRxrxcAAFgvKn5auuuuu/TKK69o69atmjFjxohtlyxZIkmqra2diNIAAEAUs3SPjGmaWr16tTZt2qRt27Zp9uzZIZeprq6WJE2fPn2cq5sc6tydOtLi4yaTAABbsjTIrFq1Shs2bNCvfvUrpaWlqampSZLkdDqVnJysw4cPa8OGDfrYxz6m7Oxs7d+/X/fcc4/Ky8t1xRVXWFm67bX6urVmY7Wqatz+aeXFLlVWlMiZkmBhZQAAhM/Sg30Nwwg6/bnnntMdd9yhxsZG/du//ZsOHDggr9ergoICffKTn9T9998f9oG74R4sNNnc9uwe7aj1qG/Q2+8wDJUV5ejFlYstrAwAAJsc7BsqQxUUFGj79u0TVM3kUefuDNgTM6DPNFVV41a9x8vPTAAAW4iKg30xsY60DH8dHklqaPaOOB8AgGhBkJmEZmWljDi/MJu9MQAAeyDITEJzXKkqL3bJMeQYJYdhqLzYxc9KAADbIMhMUpUVJSorygmYVlaUo8qKEosqAgBg9KLmyr6YWM6UBL24crHqPV41NHu5jgwAwJYIMpPc7BwCDADAvvhpCQAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2JalQWbdunW6+uqrlZaWpmnTpmnFihU6dOhQQJuzZ89q1apVys7OVmpqqm6++WadPHnSoooBAEA0sTTIbN++XatWrdKuXbv0+uuvq6enRzfddJO8Xq+/zT333KPf/OY3+vnPf67t27fr+PHj+tSnPmVh1QAAIFoYpmmaVhcxwO12a9q0adq+fbvKy8vV1tYml8ulDRs26JZbbpEkHTx4UJdddpl27typa665JmSf7e3tcjqdamtrU3p6+ni/hEmpzt2pIy0+FWZP1eycqVaXEyCaaxuLWHs9ADCccL+/4yewppDa2tokSVlZWZKkvXv3qqenR0uXLvW3ufTSSzVz5sxhg0xXV5e6urr8z9vb28e56smr1detNRurVVXj9k8rL3apsqJEzpQECyuL7trGItZeDwBEStQc7Nvf36+7775bZWVlmj9/viSpqalJiYmJysjICGibm5urpqamoP2sW7dOTqfT/ygoKBjv0ietNRurtaPWEzBtR61Hqzfus6iiv4vm2sYi1l4PAERK1ASZVatW6cCBA3rppZcuqJ+1a9eqra3N/2hsbIxQhRiszt2pqhq3+ob8Mtlnmqqqcave4x1myfEXzbWNRay9HgCIpKgIMnfddZdeeeUVbd26VTNmzPBPz8vLU3d3t1pbWwPanzx5Unl5eUH7SkpKUnp6esADkXekxTfi/IZm675co7m2sYi11wMAkWRpkDFNU3fddZc2bdqkP/7xj5o9e3bA/EWLFikhIUFbtmzxTzt06JCOHj2q0tLSiS4Xg8zKShlxfmG2dQeiRnNtYxFrrwcAIsnSILNq1Sr993//tzZs2KC0tDQ1NTWpqalJZ86ckSQ5nU6tXLlS9957r7Zu3aq9e/fq85//vEpLS8M6YwnjZ44rVeXFLjkMI2C6wzBUXuyy9IyaaK5tLGLt9QBAJFkaZNavX6+2tjZdf/31mj59uv/x8ssv+9s89thj+qd/+ifdfPPNKi8vV15enn75y19aWDUGVFaUqKwoJ2BaWVGOKitKLKro76K5trGItdcDAJESVdeRGQ9cR2b81Xu8amj2RuW1TaK5trGItdcDAMMJ9/s7IkGmr69Pb7/9tmbNmqXMzMwL7S6iCDIAANhPuN/fY/pp6e6779azzz4r6VyI+fCHP6yFCxeqoKBA27ZtG1PBAAAAozWmIPOLX/xCCxYskCT95je/UX19vQ4ePKh77rlH3/zmNyNaIAAAwHDGFGQ8Ho//Oi6/+93vdOutt+riiy/WF77wBb399tsRLRAAAGA4Ywoyubm5euedd9TX16dXX31VH/nIRyRJPp9PDocjogUCAAAMZ0w3jfz85z+vT3/605o+fboMw/Df1HH37t269NJLI1ogAADAcMYUZL797W9r/vz5amxs1K233qqkpCRJksPh0H333RfRAgEAAIbDdWQAAEDUCff7O+w9Mk888UTYg69ZsybstgAAAGMV9h6ZoTd0dLvd8vl8ysjIkCS1trYqJSVF06ZNU11dXcQLHSv2yAAAYD8RvyBefX29//G9731PV155pd599121tLSopaVF7777rhYuXKjvfOc7EXkBAAAAoYzpGJm5c+fqF7/4hUpKAm9Yt3fvXt1yyy2qr6+PWIEXij0yAADYT8SPkRnsxIkT6u3tPW96X1+fTp48OZYuMcHq3J060uIbt5sPjnf/AABIYwwyN954o7785S/rmWee0cKFCyWd2xtz5513+q8pg+jU6uvWmo3Vqqpx+6eVF7tUWVEiZ0pC1PcPAMBgY7qy709/+lPl5eXpqquuUlJSkpKSkrR48WLl5ubqmWeeiXSNiKA1G6u1o9YTMG1HrUerN+6zRf8AAAw2pj0yLpdLv/vd7/Tee+/p4MGDkqRLL71UF198cUSLQ2TVuTsD9pQM6DNNVdW4Ve/xXtDPQOPdPwAAQ40pyAwoLCyUaZqaO3eu4uMvqCtMgCMtvhHnNzRfWNAY7/4BABhqTD8t+Xw+rVy5UikpKZo3b56OHj0qSVq9erUeeeSRiBaIyJmVlTLi/MLsCwsZ490/AABDjSnIrF27Vn/961+1bds2TZkyxT996dKlevnllyNWHCJrjitV5cUuOQwjYLrDMFRe7LrgvSXj3T8AAEONKchs3rxZTz75pK677joZg7605s2bp8OHD0esOEReZUWJyopyAqaVFeWosqJkmCWiq38AAAYb04Etbrdb06ZNO2+61+sNCDaIPs6UBL24crHqPV41NHsjfp2X8e4fAIDBxrRH5qqrrtJvf/tb//OB8PLMM8+otLQ0MpVhXM3OmaobLpk2biFjvPsHAEAa4x6Z//iP/9Dy5cv1zjvvqLe3Vz/60Y/0zjvv6M9//rO2b98e6RoBAACCGtMemeuuu05//etf1dvbqw996EP6wx/+oGnTpmnnzp1atGhRpGsEAAAIatR7ZHp6evTlL39Z3/rWt/Rf//Vf41ETAABAWEa9RyYhIUH/+7//Ox61AAAAjMqYflpasWKFNm/eHOFSAAAARmdMB/sWFxfr4Ycf1o4dO7Ro0SJNnRp4ZsqaNWsiUhwAAMBIDNM0zdEuNHv27OE7NAzV1dVdUFGR1N7eLqfTqba2NqWnp1tdDgAACEO4399j2iNTX1/v//dADuJCeAAAYKKN6RgZSXr22Wc1f/58TZkyRVOmTNH8+fP1zDPPRLI2AACAEY1pj8wDDzygH/7wh1q9erX/Sr47d+7UPffco6NHj+rhhx+OaJEAAADBjOkYGZfLpSeeeEIVFRUB0zdu3KjVq1fL4/FErMALxTEyAADYz7geI9PT06OrrrrqvOmLFi1Sb2/vWLoELFPn7tSRFh83uAQAGxpTkPnc5z6n9evX64c//GHA9J/85Cf67Gc/G5HCgPHW6uvWmo3Vqqpx+6eVF7tUWVEiZ0qChZUBAMI1pp+WVq9erRdffFEFBQW65pprJEm7d+/W0aNHddtttykh4e9fAkPDzkTjpyUM57Zn92hHrUd9gz4CDsNQWVGOXly52MLKAADj+tPSgQMHtHDhQknS4cOHJUk5OTnKycnRgQMH/O04JRvRqs7dGbAnZkCfaaqqxq16j5efmQDABsYUZLZu3RqRwauqqvToo49q7969OnHihDZt2qQVK1b4599xxx164YUXApZZtmyZXn311YiMj8nrSItvxPkNzQQZALCDMV9HJhK8Xq8WLFigp556atg2H/3oR3XixAn/Y+PGjRNYIWLVrKyUEecXZhNiAMAOxrRHJlKWL1+u5cuXj9gmKSlJeXl5E1QRJos5rlSVF7uGPUaGvTEAYA+W7pEJx7Zt2zRt2jRdcskluvPOO9Xc3Dxi+66uLrW3twc8gGAqK0pUVpQTMK2sKEeVFSUWVQQAGC1L98iE8tGPflSf+tSnNHv2bB0+fFjf+MY3tHz5cu3cuVMOhyPoMuvWrdNDDz00wZXCjpwpCXpx5WLVe7xqaPZyHRkAsKExnX49HgzDOO9g36Hq6uo0d+5cvfHGG7rxxhuDtunq6lJXV5f/eXt7uwoKCjj9GgAAGwn39Ouo/2lpsDlz5ignJ0e1tbXDtklKSlJ6enrAAwAAxCZbBZljx46publZ06dPt7oUAAAQBSw9RqazszNg70p9fb2qq6uVlZWlrKwsPfTQQ7r55puVl5enw4cP62tf+5qKioq0bNkyC6sGAADRwtIg8+abb+qGG27wP7/33nslSbfffrvWr1+v/fv364UXXlBra6vy8/N100036Tvf+Y6SkpKsKhkAAESRqDnYd7xwryUAAOwnJg/2BQAAGIwgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbCuqbxqJ8VXn7tSRFp8tbpYYyVoH92Wa5riuAzutYwCwI4LMJNTq69aajdWqqnH7p5UXu1RZUSJnSoKFlZ0vkrUG62uwSK4DO61jALAzflqahNZsrNaOWk/AtB21Hq3euM+iioYXyVqD9RWJfsMdK1rXMQDYGUFmkqlzd6qqxq2+IRd07jNNVdW4Ve/xWlTZ+SJZ63B9XWi/oxkrGtcxANgdQWaSOdLiG3F+Q3P0fMlGstZQfY2137GMFU3rGADsjiAzyczKShlxfmF29ByQGslaQ/U11n7HMlY0rWMAsDuCzCQzx5Wq8mKXHIYRMN1hGCovdkXVmTWRrHW4vi6039GMFY3rGADsjiAzCVVWlKisKCdgWllRjiorSiyqaHiRrDVYX5HoN9yxonUdA4CdGaY5wtGPMSDc24BPRvUerxqavba4xkkkax3cl6RxXQd2WscAEE3C/f4myAAAgKgT7vc3Py0BAADbIsgAAADbIsgAAADbIsgAAADbIsgAAADbIsgAAADbIsgAAADbIsgAAADbIsgAAADbIsgAAADbIsgAAADbire6AFijzt2pIy2+YW9mGGr+RLYfru1oxxxvA/U4DEN9pumvK9rqBIBYQpCZZFp93VqzsVpVNW7/tPJilyorSuRMSQg5f7T9XUj74dp+d8V83b/5QNhjjrdgdQ7ITEnQaV+P/7mVdQJALOLu15PMbc/u0Y5aj/oGve0Ow1BZUY5eXLk45PzR9nch7Ydrm54cr/YzvWGPOd6C1TkcK+sEADvh7tc4T527U1U17vO+cPtMU1U1blW95x5xfr3HO6r+LqT9SG1P+3rCHnO8DVfncKyqEwBiFUFmEjnS4htx/r7G0yPOb2gO/PIN1d+FtA/VNtwxx5td6gSAWEWQmURmZaWMOL+kIHPE+YXZgQeqhurvQtqHahvumOPNLnUCQKwiyEwic1ypKi92yWEYAdMdhqHyYpfKL3aNOH/oGTeh+ruQ9iO1zUxJCHvM8TZcncOxqk4AiFUEmUmmsqJEZUU5AdPKinJUWVES1vzR9nch7Ydr++tV141qzPEWrM4BmUPOTrKyTgCIRZy1NEnVe7xqaPYOe22TUPMnsv1wbUc75ngbqCc+zlBv/9+vIxNtdQKAHYT7/U2QAQAAUccWp19XVVXp4x//uPLz82UYhjZv3hww3zRNPfDAA5o+fbqSk5O1dOlS1dTUWFMsAACIOpYGGa/XqwULFuipp54KOv/73/++nnjiCT399NPavXu3pk6dqmXLluns2bMTXCkAAIhGlt6iYPny5Vq+fHnQeaZp6vHHH9f999+vT3ziE5KkF198Ubm5udq8ebP+5V/+JehyXV1d6urq8j9vb2+PfOEAACAqRO1ZS/X19WpqatLSpUv905xOp5YsWaKdO3cOu9y6devkdDr9j4KCgokoFwAAWCBqg0xTU5MkKTc3N2B6bm6uf14wa9euVVtbm//R2Ng4rnUCAADrxNzdr5OSkpSUlGR1GQAAYAJE7R6ZvLw8SdLJkycDpp88edI/DwAATG5RG2Rmz56tvLw8bdmyxT+tvb1du3fvVmlpqYWVAQCAaGHpT0udnZ2qra31P6+vr1d1dbWysrI0c+ZM3X333frud7+r4uJizZ49W9/61reUn5+vFStWWFc0AACIGpYGmTfffFM33HCD//m9994rSbr99tv1/PPP62tf+5q8Xq++9KUvqbW1Vdddd51effVVTZkyxaqSAQBAFOEWBQAAIOqE+/0dc2ctIbQ6d6eOtPjG/SaGI40zdF6o5+Mtmseb6NoAwE4IMpNIq69bazZWq6rG7Z9WXuxSZUWJnCkJEzKOKfO8eZkpCTrt6xn2+XjUGE6tVo830bUBgB3x09Ikctuze7Sj1qO+QW+5wzBUVpSjF1cunpBxJJ03L5TxqDGcWq0eb6JrA4BoYou7X2Pi1Lk7VVXjPi9A9Jmmqmrcqvd4J2ScYPNCiXSNAyZqnYxlvImuDQDsiiAzSRxp8Y04v6E5Ml+Moca5EJGqccBErZOxjDfRtQGAXRFkJolZWSkjzi/MjsxBpKHGuRCRqnHARK2TsYw30bUBgF0RZCaJOa5UlRe75DCMgOkOw1B5sStiZ8OEGifYvFAiXeOAiVonYxlvomsDALsiyEwilRUl/gNuB5QV5aiyomTCxgk2L3PIGThDn49HjeHUavV4E10bANgRZy1NQvUerxqaveN+XZKRxhk6L9Tz8RbN4010bQAQDcL9/ibIAACAqMPp1wAAIOYRZAAAgG0RZAAAgG0RZAAAgG0RZAAAgG0RZAAAgG0RZAAAgG0RZAAAgG0RZAAAgG0RZAAAgG0RZAAAgG3FW10Axl+du1NHWnyW33QwUnVsP3RK1cdatXBmpv6h2BX2WHXuTu2ub5Zk6Jo52ZqdM3XYduHWObTtwHOHIfWZGrbfo83egNcQrLZQY55qO6umjrOanp6sPrM/5LITsR1Ey7YGYPIgyMSwVl+31mysVlWN2z+tvNilyooSOVMSbFfHkWavVjy1Q6d9Pf5pmSkJ+vWq61SQnTLsWKVzstXXb2pPQ0tAf87keLWd6fU/v3ZutkxT2lnXHLLOYONkpiQE1DZ4fMOQ/ny4+bx5khQfZ6i33zxvmaf/bVHAuMHGDGboshOxHUTLtgZg8uHu1zHstmf3aEetR32D3mKHYaisKEcvrlxsuzpKHv5D0KCQmZKgfQ/cNOxYF2K4OiM9TjDlxa6AcUcz5uBlJ2I7iJZtDUDs4O7Xk1ydu1NVNe7zvvT6TFNVNW7Ve7y2qmP7oVNBQ4wknfb16E817mHHuhDB6hyPcYIZPO5oxxxYdiK2g2jZ1gBMTgSZGHWkxTfi/IbmiflyiVQd1cdaR5z/1tHTIce6EIPrHM9xhht3LGM2NHsnZDuIlm0NwOTEMTIxalZWyojzC7Mn5kDMSNVx5YyMEecvnJmpizKSwy1r1AbXGeo1jce4YxmzMHuqQv1yHIntIFq2NQCTE3tkYtQcV6rKi11yGEbAdIdhqLzYNWFnlESqjg9fMk2Zwxw0mpmSoH8odg071oUIVud4jBPM4HFHO+bAshOxHUTLtgZgciLIxLDKihKVFeUETCsrylFlRYkt6/j1quvOCzMDZy2NNFbpnGwtLsw6rz9ncuAOyWvnZqt0TnZYdQYbZ7igVTonW9fOzQ46Tzp31lKwZYaOG2zM4cYbvOxEbAfRsq0BmHw4a2kSqPd41dDstfzaHpGq4081br119PSI15EJNla9x6vddc0yJf/1VoZrF26dQ9sOPB84pXq4fo+d9gW8hmC1hRrT09Gl421nlO9MVl+/GXLZidgOomVbA2B/4X5/E2QAAEDU4fRrAAAQ8wgyAADAtggyAADAtggyAADAtggyAADAtggyAADAtqI6yHz729+WYRgBj0svvdTqsgAAQJSI+nstzZs3T2+88Yb/eXx81JcMAAAmSNSngvj4eOXl5YXdvqurS11dXf7n7e3t41EWAACIAlH905Ik1dTUKD8/X3PmzNFnP/tZHT16dMT269atk9Pp9D8KCgomqFIAADDRovoWBb///e/V2dmpSy65RCdOnNBDDz2k999/XwcOHFBaWlrQZYLtkSkoKOAWBQAA2EhM3muptbVVs2bN0g9/+EOtXLkyrGW415L91Lk7daTFF9EbD4bqc6T5wea9vOeodtY3q2xujm69qiCg3am2s2rqOKuFMzN1UUaydte3yJC05IMbOvrbtZ/VuyfaFWcYykxJUE7aFF0zJ1umaQatpc7d6e8rPyNZfea5m1IGaz+47dBxHYah91vP6L2mDpkydfn0dLnSpwT0da62DuWkJur/XZEf0O/gsYZ77jAMf32jeQ8H9zdQSzh9jcc2A8Ba4X5/R/0xMoNlZGTo4osvVm1trdWlYBy0+rq1ZmO1qmrc/mnlxS5VVpTImZIwLn2ONN+Ued68K2c4deB4u3r7z+X/zfuOa+0v39YVF6Xrrca2kPVkJCeo9UxP2PWXF7v03RXz9PX/fVs765pDtr92bra6e/v15pHTFzTuYD/4w3taWOBUUkJ8QA2ZKQk67esZ9vng1xDqPQz2PgQztK/x2GYA2Iut9sh0dnZq5syZ+va3v601a9aEtQx7ZOzjtmf3aEetR32DNkmHYaisKEcvrlw8Ln2ONF/SefMmmsMwlJ4cHzQg2EU472Gw9yGcvsZjmwEQHWLi7tf//u//ru3bt6uhoUF//vOf9clPflIOh0MVFRVWl4YIq3N3qqrGfd4XWZ9pqqrGrXqPN+J9Vr13auT5QeZNtD7TtHWIkUK/h8O9T6H6Go9tBoD9RHWQOXbsmCoqKnTJJZfo05/+tLKzs7Vr1y65XC6rS0OEHWnxjTi/oXn0X0qh+tzX2DrqPjF2w72Hod6n4foaj20GgP1E9TEyL730ktUlYILMykoZcX5h9ugP4AzVZ0lBxqj7xNgN9x6Gep+G6yvUr+Jj2WYA2E9U75HB5DHHlaryYpcchhEw3WEYKi92jelMlFB9ll88beT5QeZNNMcHZzTZWaj3cLj3KVRf47HNALAfggyiRmVFif8g2wFlRTmqrCgZtz5Hmh9sXskMp+LjAr844+MMLSxwhlVPRvLoQklZUY5+veo6lc7JDqv9tXOzdfWszAsed6iFBRnn1TA0YA0XuMJ5D4Ot63D6Go9tBoC92OqspbHgrCX7qfd41dDsjeg1QUL1OdL8YPN+/majdhz2BFxHZqCdp6NLx9vOaOHMTM3ITNGuuuaA67kMbvfO8XYZhpQ1NVHZqUm65oOwEKyWeo/X39dFmcnq7Tf9P58MbT+47dBx4+MMHTt9RrUnO9RnmpqX71ROWlJAX56OLv3tePt515EZui6Gex4fZ/jrG817OLi/gVrC6Ws8thkA1orJC+KNBUEGAAD7iYnTrwEAAEZCkAEAALZFkAEAALZFkAEAALZFkAEAALZFkAEAALZFkAEAALZFkAEAALZFkAEAALZFkAEAALZFkAEAALYVb3UBgBXq3J060uJTYfZUmabp//fQGw4Objcwb/uhU6o+1qrpzilypU3xzxva5yv7j6vF26MbL5umfyh2DdvfwDSHIb3felaGpPyMKeozJYdhqM805TAMVTeeVoPHq9muVC2YkaE+8+83Utx+6JS2HnIrzpAuzk1TfkbyB8tJfab8Ne2ub5EhyZDU1HHW/xpOtZ/Vuyfa1XG2V6ZpypmcqEvy0rRkTrb21DVrc/VxpSXF63PXztJFGcl6Zf8JNXi8ykhJUKuvW01tZ3Xp9HR9rrTQ/9rffr9NaUkJ+uTCi/QPxS7VuTu17nfv6mBTuy6f7tQ/XjZNb7xzUu7Obk1LTVLqFIcyUhJlSmrzdSsjJVHFuWm66IPXUpg9VXvqmrWzvlnFrjRlpSbqvZMdMk1Tl+enq6ntrN4+1iYZ0oIZGfp/V+TraLP33HuVPkV9pilPZ5dM05AhyZWWpPyMKapubNOuOo9kStcVu5Q5NUGSoThD2lXXos6zPZqRmeJfH0Pf68HbTLD3cmCZkbbBkfoc7fY8mmXr3J3+bWK4OoFox00jMam0+rq1ZmO1qmrcQeeXF7tUWVEiU+Z57RYXZurQyQ61nek9b7nMlASd9vUMO27alHhdlpeuPQ0t/mnXzs2WaUo765ov4BVJ8XFSb/8FdYFRyEhOUOuZv7/X5cUufXfFfN2/+cCw29W1c7O1/rOL5Pwg+A3dtoZuPwPboTMlYcRagvUVzrKtvm7d+d9vnbftDa4TsBp3v/4AQQaD3fbsHu2o9ahvmM3eYRgqK8qRpBHbAQMchqH05Hi1n+kdcXspL3bpxZWLQ26DA32WFeXoxZWLRxw7WF/hLHvbs3tGDPOhxgUmAne/Boaoc3eqqsY94hdIn2mqqsYdsh0woM80ddrXE3J7qapxq+q9U2FtWwPbYb3HO2yb4bbnUMsOLDdSnSONC0QbggwmjSMtPqtLwCS3r7F1VO0bmocPFKG25+GWDedzMNK4QLQhyGDSmJWVYnUJmORKCjJG1b4we/iDb0Ntz8MtG87nYKRxgWhDkMGkMceVqvJilxyGMWwbh2GovNgVsh0wwGEYykxJCLm9lBe7VH7xtLC2rYHtcKSziIbbnkMtO7DcSHVy9hLshCCDSaWyosR/MG8wZUU5qqwoCdpucWGWnMnBr1iQGeIsj/Qp8VpcmBUw7dq52Sqdkx1m5cOL51M8oTKSA9/rsqIc/XrVdSNuV9fOzVZlRYmk4Nvg0O1nYDsMJVhf4SxbWVESdNsbXCdgF5y1hEmp3uNVQ7PXvwt94N9D/090cLuBeX+qceuto6eV70xWTlqSf97QPn+7/7g8nd0B15EJ1t/AtPg4Q8dOn5Eh6aLMZPX2m4qPM/z//Wtjq+rcXs12TdWVBRnq7f/7dWT+VOPWlndPyWHo3LVXhiw/UNOuumYZkhxxho63nfG/Bk9Hl/52vE2dZ3vVb0oZKQm6OPfcdVPebGjRprfeV+oH15GZkZmiV/YfV73bq8yUBLWd6dHx1rO6dHqaPlda6H/tfz3WGnAdmXqPV4/87l29c6JNl093aunluXr9nSa5O7o1LS1JU5McykxJlCSd9nUrMyVRRblpmvHBaynMnqo3G1q047BHF09LU+bURNWe7FCfaWpevlMn289qf2ObTMP0X0fm2Gmf/73q6//gOjI6dx2dnNQkXZSZrOrGVu067JH5wXVksqaeu5ZNfJyhnYeb1dl17joyA+tj6Hs9eJsJ9l4Od32WoX0M1+dot+fRLFvv8fq3Ca4jg2jD6dcfIMgAAGA/nH4NAABiHkEGAADYFkEGAADYFkEGAADYFkEGAADYFkEGAADYFkEGAADYFkEGAADYFkEGAADYFkEGAADYFkEGAADYVvBb+UaZp556So8++qiampq0YMECVVZWavHixVaXBYxZnbtTR1p8Yd3kb2jb4ZYdTZ9D+99d3yJDUn7GFPWZUmH2VJmm6e/PNM1zN4r0+DTXNVWu1CQ1dZzVdOcU9fWfuwmjIWlXfbNOtp9VXnqyrpmTpT5Teq+pXcdOn5EkzchMkQzJNE1dPt0pV3qSHIa09eApHWrq0PSMZKUnJ+hQU7vOdPdrjmuq0pMT1ObrkYxzdxG/fHq63j3RoerG02o706OeXlNpSfG6KDtZvrN9avZ1aX6+U/98Zb76TOlU21ltrn5frb4elc3NUdnFOSrMnqqjzV5tPXRKOalJumKGU32m5DAMvd96Rp6Os8pJm6Jr5mRrT12zXn/3pOIMQ/kZyXIY0t4jrfJ4z8qQlOiIU176FF0yPV2ZUxMl01SLt0fVR0/L09mtxPg4JcYbKshK0dLL8uRKT9KptrNq6jirhTMz/TcU3X7olLYecitOUuNpnxpP+zQza6pSkxyq93hlyFBhzlRdMzdLrrQpchjS5n3HVd14Wj19/UpwxCnfmSxvd68aW3xKTozXNbOz1Gua6uzq1YyMFLWd7dG7x9tkyNBl+elKn3JuXUuGyopydPlF6frb+23aUevRqfYuJcYbmn9Rhlp93fJ29+qTJTPU329qZ32zyubm+P995myvjrT41NPXr8T4OBVkTtXnrp2lfyh26eU9R/X6uyflTE6QKan9TI9KCjLVZ/ar3uNTZkrCB+tNkkzJMNTS2S1Jumx6mlzpUwK2/Sf/WKM6j1c3XZ6n/++GIv/2+5f6ZrWf7dWyeXm69aoC//rMSU30v78D/Zybd0pxhuG/Gahpmtpd3+y/iengm2gO/owEu7nm4M/euX5a/NvQRRnJqm48rRZvt268LNf/fo/l87r90ClVH2v1bzcDfTgMQ32mKYehgNd5oX8bwln2QvqOlKi/aeTLL7+s2267TU8//bSWLFmixx9/XD//+c916NAhTZs2LeTy3DQS0aTV1601G6tVVeP2TysvdqmyokTOlISQbTNTEnTa1xOw7HdXzNf9mw+E1efQ/u/877e0s645Ei8NY5Q+JV6mKXV09VpdSlRLTYxTZ3f/BfczcEf4cCwuzJRpSn85cjpg+rVzs7X+s4tkyjzvMxpK+pR4XZKbFtBnqM/rkWavVjy1I+CzH+p1nPvbME/3b/7bqP82SKH/Vo3mb9lYxczdr5csWaKrr75aTz75pCSpv79fBQUFWr16te67776QyxNkEE1ue3aPdtR61DfoY+cwzv3f8IsrF4dsO5TDMJSeHK/2M71h9Tm0/9H8AQbwd+Uf7FkJ9RkNR6jPa8nDfwgIMeH2Oda/DVLov1Wj+Vs2VjFx9+vu7m7t3btXS5cu9U+Li4vT0qVLtXPnzqDLdHV1qb29PeABRIM6d6eqatzn/dHrM01V1bhV7/GGbDtUn2nqtK8nrD6D1QJgbKpq3GF9RsMx0ud1+6FTow4xA32O5W+DFPpvVdV7wV97OH2Ph6gOMh6PR319fcrNzQ2Ynpubq6ampqDLrFu3Tk6n0/8oKCiYiFKBkI60+Eac39D89w9/qLbhGtznaGoBMPGCfV6rj7VO2FgDQv192Nd4esT5I/U9HqI6yIzF2rVr1dbW5n80NjZaXRIgSZqVlTLi/MLsvx8oF6ptuAb3OZpaAEy8YJ/XK2dkTNhYA0L9fSgpyBxz3+MhqoNMTk6OHA6HTp48GTD95MmTysvLC7pMUlKS0tPTAx5ANJjjSlV5sUsOwwiY7jAMlRe7Ao74H67tUA7DUGZKQlh9BqsFwNiUF7vC+oyGY6TP64cvmabMMRw8O9a/DVLov1XlFwd/7eH0PR6iOsgkJiZq0aJF2rJli39af3+/tmzZotLSUgsrA8amsqJEZUU5AdPKinJUWVESVtuhf9DKinL061XXhd3n0P5L52SPpnyMg/Qp8UpLssWVMCyVmhiZr6v4uPCDx+LCLF096/y9D9fOzVZlRUnQz2go6VPiz+sz1Of116uuO++zH+p1XMjfBin036rR/C0bb1F/1tLLL7+s22+/XT/+8Y+1ePFiPf744/qf//kfHTx48LxjZ4LhrCVEo3qPVw3N3rCuvTC07XDLjqbPof3vqmuWIemizGT19pv+XcMD/UnSb/cfV53bqzmuqcpNn6LjbWeU7zzX3pDkiDO083CzmtrPKC89WaVzs9Xbb6r2ZIcaW85IhqkZmSkydO6gwHn5TuWkJSk+ztC2g6f07okOTc+YoozkBB1s6pCvu09zP7iOTKuvR3GGlDolXvPynXr3RLv2HT2t9jO96u7rP3cdmaxknenqU7O3S/PynfpEyUXq7Tfl6ejS5n3vq8XbrbKiHJUVn7uOzLHTPm1596RyUpO0oCBDvf2m4uMMHTt9Rs2dXcpOTdI1c7L1ZkOLXvtbkxyGoekZyYqPM7T3yGl5Os9KMpToMJTnnKJL89KVMTVRhqQWb7eqj7bK3dl17joyjnPXkfnI5XnKSUuSp6NLx9vOBFxH5k81bm1595QchnTs9BkdbfFqZtZUpU2J12F3p/86MqVzs/3r7Vf73te+o6fV02cqwWFoujNZvu4+HW3xnruOzJxs9fWb6uzq0YzMFLWf6dHfPriOzOX56R9cs6dDpildV5yjy/LT9e7xdv1fjVsnO7qU6Dh3HZm2M93q7Dp3HRlJ2nHYo7K5Of5/n+3u05Fmr7p7TSXFx2lGZor/OjI/f7NRf/hbk9KTE2RIaj3To4UzM9Vvmqpze5WZkuBfb9K5y8m0ervVb0qX56crJy0pYNt/6o81qnV3+q8jM7D97m1oUeuZHv91ZAbWZ05qov/9HejnTzVu/fHdUzIM+a8jI0m765rl7uySKzXwOjKDPyPBriMz+LMnSbvqmv3b0IzMZP21sVWezq6A68iM5fP6pxq33jp62r/dDPQxcCr2wH8j9bchnGUvpO9QYub0a0l68skn/RfEu/LKK/XEE09oyZIlYS1LkAEAwH5iKshcCIIMAAD2ExPXkQEAABgJQQYAANgWQQYAANgWQQYAANgWQQYAANgWQQYAANgWQQYAANgWQQYAANgWQQYAANhWzN+pbODCxe3t7RZXAgAAwjXwvR3qBgQxH2Q6OjokSQUFBRZXAgAARqujo0NOp3PY+TF/r6X+/n4dP35caWlpMozwb9+O87W3t6ugoECNjY3ctyoCWJ+RxzqNPNZpZLE+w2eapjo6OpSfn6+4uOGPhIn5PTJxcXGaMWOG1WXElPT0dD6AEcT6jDzWaeSxTiOL9RmekfbEDOBgXwAAYFsEGQAAYFsEGYQtKSlJDz74oJKSkqwuJSawPiOPdRp5rNPIYn1GXswf7AsAAGIXe2QAAIBtEWQAAIBtEWQAAIBtEWQAAIBtEWQwag0NDVq5cqVmz56t5ORkzZ07Vw8++KC6u7utLs3Wvve97+naa69VSkqKMjIyrC7Hlp566ikVFhZqypQpWrJkifbs2WN1SbZVVVWlj3/848rPz5dhGNq8ebPVJdnaunXrdPXVVystLU3Tpk3TihUrdOjQIavLigkEGYzawYMH1d/frx//+Mf629/+pscee0xPP/20vvGNb1hdmq11d3fr1ltv1Z133ml1Kbb08ssv695779WDDz6ot956SwsWLNCyZct06tQpq0uzJa/XqwULFuipp56yupSYsH37dq1atUq7du3S66+/rp6eHt10003yer1Wl2Z7nH6NiHj00Ue1fv161dXVWV2K7T3//PO6++671draanUptrJkyRJdffXVevLJJyWdu89aQUGBVq9erfvuu8/i6uzNMAxt2rRJK1assLqUmOF2uzVt2jRt375d5eXlVpdja+yRQUS0tbUpKyvL6jIwSXV3d2vv3r1aunSpf1pcXJyWLl2qnTt3WlgZEFxbW5sk8XczAggyuGC1tbWqrKzUl7/8ZatLwSTl8XjU19en3NzcgOm5ublqamqyqCoguP7+ft19990qKyvT/PnzrS7H9ggy8LvvvvtkGMaIj4MHDwYs8/777+ujH/2obr31Vn3xi1+0qPLoNZZ1CiC2rVq1SgcOHNBLL71kdSkxId7qAhA9vvrVr+qOO+4Ysc2cOXP8/z5+/LhuuOEGXXvttfrJT34yztXZ02jXKcYmJydHDodDJ0+eDJh+8uRJ5eXlWVQVcL677rpLr7zyiqqqqjRjxgyry4kJBBn4uVwuuVyusNq+//77uuGGG7Ro0SI999xziotj514wo1mnGLvExEQtWrRIW7Zs8R+Q2t/fry1btuiuu+6ytjhAkmmaWr16tTZt2qRt27Zp9uzZVpcUMwgyGLX3339f119/vWbNmqUf/OAHcrvd/nn83+/YHT16VC0tLTp69Kj6+vpUXV0tSSoqKlJqaqq1xdnAvffeq9tvv11XXXWVFi9erMcff1xer1ef//znrS7Nljo7O1VbW+t/Xl9fr+rqamVlZWnmzJkWVmZPq1at0oYNG/SrX/1KaWlp/mO3nE6nkpOTLa7O5kxglJ577jlTUtAHxu72228Puk63bt1qdWm2UVlZac6cOdNMTEw0Fy9ebO7atcvqkmxr69atQbfH22+/3erSbGm4v5nPPfec1aXZHteRAQAAtsWBDQAAwLYIMgAAwLYIMgAAwLYIMgAAwLYIMgAAwLYIMgAAwLYIMgAAwLYIMgAAwLYIMgAAwLYIMgCiVl9fn6699lp96lOfCpje1tamgoICffOb35QkrVmzRosWLVJSUpKuvPJKCyoFYBWCDICo5XA49Pzzz+vVV1/Vz372M//01atXKysrSw8++KB/2he+8AV95jOfsaJMABbi7tcAotrFF1+sRx55RKtXr9Y//uM/as+ePXrppZf0l7/8RYmJiZKkJ554QpLkdru1f/9+K8sFMMEIMgCi3urVq7Vp0yZ97nOf09tvv60HHnhACxYssLosAFGAIAMg6hmGofXr1+uyyy7Thz70Id13331WlwQgSnCMDABb+OlPf6qUlBTV19fr2LFjVpcDIEoQZABEvT//+c967LHH9Morr2jx4sVauXKlTNO0uiwAUYAgAyCq+Xw+3XHHHbrzzjt1ww036Nlnn9WePXv09NNPW10agChAkAEQ1dauXSvTNPXII49IkgoLC/WDH/xAX/va19TQ0CBJqq2tVXV1tZqamnTmzBlVV1erurpa3d3dFlYOYCIYJvtnAUSp7du368Ybb9S2bdt03XXXBcxbtmyZent79cYbb+iGG27Q9u3bz1u+vr5ehYWFE1QtACsQZAAAgG3x0xIAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALCt/x8CtUgrUmu9cQAAAABJRU5ErkJggg==",
      "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": 51,
   "id": "2a72bfad-e09e-42c0-9273-b87d095f8f1c",
   "metadata": {},
   "outputs": [],
   "source": [
    "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": 52,
   "id": "b70964eb-dbbc-4d90-9035-4efbced2d0e9",
   "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",
    "        self.raw_r = torch.nn.Parameter(torch.ones(*batch_shape, 1))\n",
    "        r_constraint = Positive()\n",
    "        self.register_constraint(\"raw_r\", r_constraint)\n",
    "\n",
    "        self.raw_mu = torch.nn.Parameter(torch.ones(*batch_shape, 1))\n",
    "        mu_constraint = Positive()\n",
    "        self.register_constraint(\"raw_mu\", mu_constraint)\n",
    "\n",
    "        self.raw_pi = torch.nn.Parameter(torch.ones(*batch_shape, 1))\n",
    "        pi_constraint = Interval(0, 1)\n",
    "        self.register_constraint(\"raw_pi\", pi_constraint)\n",
    "\n",
    "    @property\n",
    "    def mu(self) -> Tensor:\n",
    "        return self.raw_mu_constraint.transform(self.raw_mu)\n",
    "\n",
    "    @mu.setter\n",
    "    def mu(self, value: Tensor) -> None:\n",
    "        self._set_mu(value)\n",
    "        \n",
    "    def _set_mu(self, value: Tensor) -> None:\n",
    "        if not torch.is_tensor(value):\n",
    "            value = torch.as_tensor(value).to(self.raw_mu)\n",
    "        self.initialize(raw_mu=self.raw_mu_constraint.inverse_transform(value))\n",
    "\n",
    "    @property\n",
    "    def r(self) -> Tensor:\n",
    "        return self.raw_r_constraint.transform(self.raw_r)\n",
    "\n",
    "    @r.setter\n",
    "    def r(self, value: Tensor) -> None:\n",
    "        self._set_r(value)\n",
    "        \n",
    "    def _set_r(self, value: Tensor) -> None:\n",
    "        if not torch.is_tensor(value):\n",
    "            value = torch.as_tensor(value).to(self.raw_r)\n",
    "        self.initialize(raw_r=self.raw_r_constraint.inverse_transform(value))\n",
    "\n",
    "    @property\n",
    "    def pi(self) -> Tensor:\n",
    "        return self.raw_pi_constraint.transform(self.raw_pi)\n",
    "        \n",
    "    @pi.setter\n",
    "    def pi(self, value: Tensor) -> None:\n",
    "        self._set_pi(value)\n",
    "\n",
    "    def _set_pi(self, value: Tensor) -> None:\n",
    "        if not torch.is_tensor(value):\n",
    "            value = torch.as_tensor(value).to(self.raw_pi)\n",
    "        self.initialize(raw_pi=self.raw_pi_constraint.inverse_transform(value))\n",
    "        \n",
    "    def forward(self, function_samples: Tensor, *args: Any, **kwargs: Any) -> base_distributions.Poisson:\n",
    "        mu = torch.exp(function_samples) #* self.mu\n",
    "        probs = self.r / (mu + self.r)\n",
    "        pi = self.pi\n",
    "\n",
    "        return base_distributions.ZeroInflatedNegativeBinomial(total_count = self.r, probs=probs, gate=pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "d9f6cc96-a52a-4e52-81ec-7e90363351b2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "30\n",
      "1.1875917089786223\n",
      "60\n",
      "1.1871355744182392\n",
      "90\n",
      "1.1871188384897515\n",
      "120\n",
      "1.187111033430796\n",
      "150\n",
      "1.1871037766962085\n",
      "180\n",
      "1.1870963435314605\n",
      "210\n",
      "1.187088922242651\n",
      "240\n",
      "1.1870815668555106\n",
      "270\n",
      "1.1870743130446562\n",
      "300\n",
      "1.1870671893840936\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "\n",
    "class GPModel(ApproximateGP):\n",
    "    def __init__(self, inducing_points, feature_extractor):\n",
    "        variational_distribution = gpytorch.variational.NaturalVariationalDistribution(inducing_points.size(0))\n",
    "        variational_strategy = gpytorch.variational.VariationalStrategy(\n",
    "            self, inducing_points, variational_distribution, learn_inducing_locations=True\n",
    "        )\n",
    "        super(GPModel, self).__init__(variational_strategy)\n",
    "        self.mean_module = gpytorch.means.ZeroMean()\n",
    "        self.covar_module = gpytorch.kernels.ScaleKernel(\n",
    "            gpytorch.kernels.LinearKernel(),\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",
    "M = 300\n",
    "inducing_points = torch.normal(mean=0, std =1, size=(1,N)).double().cuda().flatten()\n",
    "model = GPModel(inducing_points=inducing_points, feature_extractor=None).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",
    "#model.covar_module.base_kernel.initialize(lengthscale=5)\n",
    "likelihood.train()\n",
    "likelihood.initialize(r=r, pi=0.35)\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 = 300\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": 54,
   "id": "dd5f17f1-1890-4657-ba61-c2aab657e0ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "model.eval()\n",
    "likelihood.eval()\n",
    "with torch.no_grad():\n",
    "    preds = likelihood(model(X))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "id": "699539a5-e487-4241-bf90-231a78b5b0ed",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Axes: xlabel='X1', ylabel='y'>"
      ]
     },
     "execution_count": 82,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjIAAAGwCAYAAACzXI8XAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAARBJJREFUeJzt3Xt4lPWd///XndPkPOREQkrCWdQqiFQwYvmKUtHuurKiv5YexMrVXS3iqu3W0m3rdrddXNvfaqWedutC3RVPbdHWrbWVStJFFEXzxWqhBAKEQyAHMpMDmZzu7x9hxkxIJpNkZu77nnk+rmsuzT0zn8/7c98zmTczc79imKZpCgAAwIGSrC4AAABgrGhkAACAY9HIAAAAx6KRAQAAjkUjAwAAHItGBgAAOBaNDAAAcKwUqwuItr6+Ph07dkw5OTkyDMPqcgAAQBhM01Rra6tKS0uVlDT8+y5x38gcO3ZMZWVlVpcBAADGoK6uTpMnTx72+rhvZHJyciT174jc3FyLqwEAAOHwer0qKysLvI4PJ+4bGf/HSbm5uTQyAAA4zEhfC+HLvgAAwLFoZAAAgGPZppG5//77ZRiG7rrrrsC2zs5OrVmzRgUFBcrOztaKFSt04sQJ64oEAAC2YotG5u2339YTTzyhOXPmBG2/++679atf/UovvPCCKisrdezYMd1www0WVQkAAOzG8kamra1Nn//85/Uf//EfysvLC2z3eDx68skn9W//9m+68sorNX/+fG3cuFFvvPGG3nzzzWHH8/l88nq9QRcAABCfLG9k1qxZo7/4i7/Q0qVLg7bv2rVL3d3dQdvPPfdclZeXa8eOHcOOt379ernd7sCFDBkAAOKXpY3Ms88+q3fffVfr168/67r6+nqlpaVpwoQJQduLi4tVX18/7Jjr1q2Tx+MJXOrq6iJdNgAAsAnLcmTq6ur0d3/3d/rd736n9PT0iI3rcrnkcrkiNh4AALAvy96R2bVrl06ePKmLL75YKSkpSklJUWVlpR5++GGlpKSouLhYXV1damlpCbrfiRMnVFJSYk3RAADAVix7R+aqq67S+++/H7TtS1/6ks4991zde++9KisrU2pqqrZu3aoVK1ZIkvbu3avDhw+roqLCipIBAIDNWNbI5OTk6IILLgjalpWVpYKCgsD21atX65577lF+fr5yc3O1du1aVVRU6NJLL7WiZMSJAw1tOtTcoakFWZpWmGV1OQCAcbD131p68MEHlZSUpBUrVsjn82nZsmV69NFHrS4LDtXS0aU7n6lW1b6GwLbFs4q0YeU8uTNTLawMADBWhmmaptVFRJPX65Xb7ZbH4+GPRia4m5/cqe01jeod8JBPNgwtmlmop1YvsLAyAMBg4b5+W54jA8TCgYY2Ve1rCGpiJKnXNFW1r0G1je0WVQYAGA8aGSSEQ80dIa8/2EQjAwBORCODhDAlPzPk9VML+NIvADgRjQwSwvSibC2eVaRkwwjanmwYWjyriLOXAMChaGSQMDasnKdFMwuDti2aWagNK+dZVBEAYLxsffo1EEnuzFQ9tXqBahvbdbCpnRwZAIgDNDJIONMKaWAAIF7w0RIAAHAsGhkAAOBYNDIAAMCxaGQAAIBj0cgAAADHopEBAACORSMDAAAci0YGAAA4Fo0MAABwLBoZAADgWDQyAADAsWhkAACAY9HIAAAAx6KRAQAAjkUjAwAAHItGBgAAOBaNDAAAcCwaGQAA4Fg0MgAAwLFoZAAAgGPRyAAAAMeikQEAAI5FIwMAAByLRgYAADhWitUFIHwHGtp0qLlDUwuyNK0wy+pyAACwnKXvyDz22GOaM2eOcnNzlZubq4qKCr3yyiuB66+44goZhhF0ue222yys2BotHV26+cmduvL/r9SXNr6tJT/cppuf3ClPR7fVpQEAYClLG5nJkyfr/vvv165du/TOO+/oyiuv1PXXX68PPvggcJsvf/nLOn78eODywAMPWFixNe58plrbaxqDtm2vadTaZ96zqCIAAOzB0o+WrrvuuqCfv//97+uxxx7Tm2++qY9//OOSpMzMTJWUlIQ9ps/nk8/nC/zs9XojU6xFDjS0qWpfw1nbe01TVfsaVNvYzsdMAICEZZsv+/b29urZZ59Ve3u7KioqAtuffvppFRYW6oILLtC6devU0dERcpz169fL7XYHLmVlZdEuPaoONYde78Gm9hhVAgCA/Vj+Zd/3339fFRUV6uzsVHZ2trZs2aLzzz9fkvS5z31OU6ZMUWlpqXbv3q17771Xe/fu1S9+8Ythx1u3bp3uueeewM9er9fRzcyU/MyQ108t4N0YAEDisryRmT17tqqrq+XxePSzn/1Mq1atUmVlpc4//3z9zd/8TeB2F154oSZNmqSrrrpK+/fv14wZM4Ycz+VyyeVyxar8qJtelK3Fs4q0vaZRvaYZ2J5sGFo0s5CPlQAACc3yj5bS0tI0c+ZMzZ8/X+vXr9fcuXP1ox/9aMjbLly4UJJUU1MTyxItt2HlPC2aWRi0bdHMQm1YOc+iigAAsAfL35EZrK+vL+jLugNVV1dLkiZNmhTDiqznzkzVU6sXqLaxXQeb2smRAQDgDEsbmXXr1unaa69VeXm5WltbtXnzZm3btk2vvvqq9u/fr82bN+vTn/60CgoKtHv3bt19991avHix5syZY2XZlplWSAMDAMBAljYyJ0+e1M0336zjx4/L7XZrzpw5evXVV/WpT31KdXV1eu211/TQQw+pvb1dZWVlWrFihb71rW9ZWTIAALARwzQHfIM0Dnm9Xrndbnk8HuXm5lpdDgAACEO4r9+Wf9kXAABgrGhkAACAY9HIAAAAx6KRAQAAjkUjAwAAHItGBgAAOBaNDAAAcCwaGQAA4Fg0MgAAwLFoZAAAgGPRyAAAAMeikQEAAI5FIwMAAByLRgYAADhWitUFIPoONLTpUHOHphZkaVphltXlAAAQMTQycaylo0t3PlOtqn0NgW2LZxVpw8p5cmemWlgZAACRwUdLcezOZ6q1vaYxaNv2mkatfeY9iyoCACCyaGTi1IGGNlXta1CvaQZt7zVNVe1rUG1ju0WVAQAQOTQycepQc0fI6w820cgAAJyPRiZOTcnPDHn91AK+9AsAcD4amTg1vShbi2cVKdkwgrYnG4YWzyri7CUAQFygkYljG1bO06KZhUHbFs0s1IaV8yyqCACAyOL06zjmzkzVU6sXqLaxXQeb2smRAQDEHRqZBDCtkAYGABCf+GgJAAA4Fo0MAABwLBoZAADgWDQyAADAsWhkAACAY9HIAAAAx6KRAQAAjkUjAwAAHMvSRuaxxx7TnDlzlJubq9zcXFVUVOiVV14JXN/Z2ak1a9aooKBA2dnZWrFihU6cOGFhxQAAwE4sbWQmT56s+++/X7t27dI777yjK6+8Utdff70++OADSdLdd9+tX/3qV3rhhRdUWVmpY8eO6YYbbrCyZAAAYCOGaZqm1UUMlJ+frx/84Ae68cYbVVRUpM2bN+vGG2+UJO3Zs0fnnXeeduzYoUsvvXTI+/t8Pvl8vsDPXq9XZWVl8ng8ys3NjckaAADA+Hi9Xrnd7hFfv23zHZne3l49++yzam9vV0VFhXbt2qXu7m4tXbo0cJtzzz1X5eXl2rFjx7DjrF+/Xm63O3ApKyuLRfkAAMACljcy77//vrKzs+VyuXTbbbdpy5YtOv/881VfX6+0tDRNmDAh6PbFxcWqr68fdrx169bJ4/EELnV1dVFeAQAAsIrlf/169uzZqq6ulsfj0c9+9jOtWrVKlZWVYx7P5XLJ5XJFsEIAAGBXljcyaWlpmjlzpiRp/vz5evvtt/WjH/1In/nMZ9TV1aWWlpagd2VOnDihkpISi6oFAAB2YvlHS4P19fXJ5/Np/vz5Sk1N1datWwPX7d27V4cPH1ZFRYWFFQIAALuw9B2ZdevW6dprr1V5eblaW1u1efNmbdu2Ta+++qrcbrdWr16te+65R/n5+crNzdXatWtVUVEx7BlLGJ0DDW061NyhqQVZmlaYZXU5AACMmqWNzMmTJ3XzzTfr+PHjcrvdmjNnjl599VV96lOfkiQ9+OCDSkpK0ooVK+Tz+bRs2TI9+uijVpYcF1o6unTnM9Wq2tcQ2LZ4VpE2rJwnd2aqhZUBADA6tsuRibRwz0NPJDc/uVPbaxrVO+DQJxuGFs0s1FOrF1hYGQAA/RyXI4PYONDQpqp9DUFNjCT1mqaq9jWotrHdosoAABg9GpkEc6i5I+T1B5toZAAAzkEjk2Cm5GeGvH5qAV/6BQA4B41MgplelK3Fs4qUbBhB25MNQ4tnFXH2EgDAUWhkEtCGlfO0aGZh0LZFMwu1YeU8iyoCAGBsLE/2Rey5M1P11OoFqm1s18GmdnJkAACORSOTwKYV0sAAAJyNj5YAAIBj0cgAAADHopEBAACORSMDAAAci0YGAAA4Fo0MAABwLBoZAADgWDQyAADAsWhkAACAY9HIAAAAx6KRAQAAjkUjAwAAHItGBgAAOBaNDAAAcCwaGQAA4Fg0MgAAwLFoZAAAgGPRyAAAAMeikQEAAI5FIwMAAByLRgYAADgWjQwAAHAsGhkAAOBYNDIAAMCxaGQAAIBjWdrIrF+/XpdccolycnI0ceJELV++XHv37g26zRVXXCHDMIIut912m0UVAwAAO7G0kamsrNSaNWv05ptv6ne/+526u7t19dVXq729Peh2X/7yl3X8+PHA5YEHHrCoYgAAYCcpVk7+m9/8JujnTZs2aeLEidq1a5cWL14c2J6ZmamSkpJYlwcAAGzOVt+R8Xg8kqT8/Pyg7U8//bQKCwt1wQUXaN26dero6Bh2DJ/PJ6/XG3QBAADxydJ3ZAbq6+vTXXfdpUWLFumCCy4IbP/c5z6nKVOmqLS0VLt379a9996rvXv36he/+MWQ46xfv17f/e53Y1U2AACwkGGapml1EZJ0++2365VXXtH//u//avLkycPe7ve//72uuuoq1dTUaMaMGWdd7/P55PP5Aj97vV6VlZXJ4/EoNzc3KrUDAIDI8nq9crvdI75+2+IdmTvuuEMvv/yyqqqqQjYxkrRw4UJJGraRcblccrlcUakTAADYi6WNjGmaWrt2rbZs2aJt27Zp2rRpI96nurpakjRp0qQoVwcAAOzO0kZmzZo12rx5s1566SXl5OSovr5ekuR2u5WRkaH9+/dr8+bN+vSnP62CggLt3r1bd999txYvXqw5c+ZYWToAALABS78jYxjGkNs3btyoW265RXV1dfrCF76gP/7xj2pvb1dZWZn++q//Wt/61rfC/r5LuJ+xAQAA+3DEd2RG6qHKyspUWVkZo2oAAIDT2CpHBgAAYDRoZAAAgGPRyAAAAMeikQEAAI5FIwMAAByLRgYAADiWLf5EAeLbgYY2HWru0NSCLE0rzLK6HABAHKGRQdS0dHTpzmeqVbWvIbBt8awibVg5T+7MVAsrAwDECz5aQtTc+Uy1ttc0Bm3bXtOotc+8Z1FFAIB4QyODqDjQ0KaqfQ3qHZTe3GuaqtrXoNrGdosqAwDEExoZRMWh5o6Q1x9sopEBAIwfjQyiYkp+ZsjrpxbwpV8AwPjRyCAqphdla/GsIiUP+gvnyYahxbOKOHsJABARNDKImg0r52nRzMKgbYtmFmrDynkWVQQAiDecfo2ocWem6h//6nztrG2WKenS6QW2fifGn3eTbEi9psi9AQAHoJFBVDgpQ2aoWv3sWjMAoB8fLSEqnJQhM1StfnatGQDQj0YGEeekDJnhavWzY80AgI/QyCDinJQhM1KtfnaqGQDwERoZRJyTMmRGqtXPTjUDAD5CI4OIc1KGzHC1+tmxZgDAR2hkEBVOypAZqlY/u9YMAOhnmOYw33KME16vV263Wx6PR7m5uVaXk3BqG9t1sKndEZks/lpTkgz19JmOqBkA4lW4r9/kyCCqphU6pxlwUq0AgH58tAQAAByLRgYAADgWjQwAAHAsGhkAAOBYNDIAAMCxaGQAAIBj0cgAAADHopEBAACOZWkjs379el1yySXKycnRxIkTtXz5cu3duzfoNp2dnVqzZo0KCgqUnZ2tFStW6MSJExZVDCsdaGjT63tPqrYxdn+J2oo5AQDhszTZt7KyUmvWrNEll1yinp4effOb39TVV1+tDz/8UFlZ/Qmrd999t/7nf/5HL7zwgtxut+644w7dcMMN2r59u5WlI4ZaOrp05zPVqtrXENi2eFaRNqycJ3dmatzMCQAYPVv9raWGhgZNnDhRlZWVWrx4sTwej4qKirR582bdeOONkqQ9e/bovPPO044dO3TppZeOOCZ/a8n5bn5yp7bXNKp3wEM12TC0aGahnlq9IG7mBAB8JNzXb1t9R8bj8UiS8vPzJUm7du1Sd3e3li5dGrjNueeeq/Lycu3YsWPIMXw+n7xeb9AFznWgoU1V+xqCGgpJ6jVNVe1riMpHPlbMCQAYG9s0Mn19fbrrrru0aNEiXXDBBZKk+vp6paWlacKECUG3LS4uVn19/ZDjrF+/Xm63O3ApKyuLdumIokPNHSGvP9gU+abCijkBAGNjm0ZmzZo1+uMf/6hnn312XOOsW7dOHo8ncKmrq4tQhbDClPzMkNdPLYj8X6u2Yk4AwNiMupFZtWqVqqqqIlrEHXfcoZdfflmvv/66Jk+eHNheUlKirq4utbS0BN3+xIkTKikpGXIsl8ul3NzcoAuca3pRthbPKlKyYQRtTzYMLZ5VpGmFkW8qrJgTADA2o25kPB6Pli5dqlmzZulf/uVfdPTo0TFPbpqm7rjjDm3ZskW///3vNW3atKDr58+fr9TUVG3dujWwbe/evTp8+LAqKirGPC+cZcPKeVo0szBo26KZhdqwcl5czQkAGL0xnbXU0NCg//qv/9JPf/pTffjhh1q6dKlWr16t66+/Xqmp4Z+a+pWvfEWbN2/WSy+9pNmzZwe2u91uZWRkSJJuv/12/frXv9amTZuUm5urtWvXSpLeeOONsObgrCXrHWho06HmDk0tyBrXuxm1je062NQ+7nHsPicAIPzX73Gffv3uu+9q48aN+slPfqLs7Gx94Qtf0Fe+8hXNmjVrxPsag96699u4caNuueUWSf2BeF/96lf1zDPPyOfzadmyZXr00UeH/WhpMBoZ65DFAgAYq5g0MsePH9dTTz2ljRs36siRI1qxYoWOHj2qyspKPfDAA7r77rvHOnTE0MhYhywWAMBYRS1Hpru7Wz//+c/1l3/5l5oyZYpeeOEF3XXXXTp27Jh++tOf6rXXXtPzzz+vf/qnfxrXAuBsZLEAAGJh1H+iYNKkSerr69PKlSu1c+dOXXTRRWfdZsmSJWdlvyCxhJPFwndOAADjNepG5sEHH9RNN92k9PT0YW8zYcIE1dbWjqswOBtZLACAWBj1R0tf/OIXQzYxgEQWCwAgNmyT7Iv4QxYLACDaRv3REhAud2aqnlq9gCwWAEDU0Mgg6qYV0sAAAKKDj5YAAIBj0cgAAADHopEBAACORSMDAAAci0YGAAA4Fo0MAABwLBoZAADgWOTIJJgDDW061NxBOB0AIC7QyCSIlo4u3flMtar2NQS2LZ5VpA0r58mdmWphZQAAjB0fLSWIO5+p1vaaxqBt22satfaZ9yyqCACA8aORSQAHGtpUta9BvaYZtL3XNFW1r0G1je0WVQYAwPjQyCSAQ80dIa8/2EQjAwBwJhqZBDAlPzPk9VML+NIvAMCZaGQSwPSibC2eVaRkwwjanmwYWjyriLOXAACORSOTIDasnKdFMwuDti2aWagNK+dZVBEAAOPH6dcJwp2ZqqdWL1BtY7sONrWTIwMAiAs0MglmWiENDMaGMEUAdkQjAyAkwhQB2BnfkQEQEmGKAOyMRgbAsAhTBGB3NDIAhkWYIgC7o5EBMCzCFAHYHY0MgGERpgjA7mhkAIREmCIAO+P0awAhEaYIwM4sfUemqqpK1113nUpLS2UYhl588cWg62+55RYZhhF0ueaaa6wpFrZ3oKFNr+89Oa4zaSIxRryaVpilJbMn0sQAsBVL35Fpb2/X3Llzdeutt+qGG24Y8jbXXHONNm7cGPjZ5XLFqjw4RCQC2wh9AwBnsrSRufbaa3XttdeGvI3L5VJJSUmMKoIThQpse2r1gpiNAQCIPdt/2Xfbtm2aOHGiZs+erdtvv11NTU0hb+/z+eT1eoMuiF+RCGwj9A0AnMvWjcw111yjp556Slu3btW//uu/qrKyUtdee616e3uHvc/69evldrsDl7KyshhWjFiLRGAboW8A4Fy2Pmvps5/9bOD/L7zwQs2ZM0czZszQtm3bdNVVVw15n3Xr1umee+4J/Oz1emlm4lgkAtsIfQMA57L1OzKDTZ8+XYWFhaqpqRn2Ni6XS7m5uUEXxK9IBLYR+gYAzuWoRubIkSNqamrSpEmTrC4FNhKJwDZC3wDAmSz9aKmtrS3o3ZXa2lpVV1crPz9f+fn5+u53v6sVK1aopKRE+/fv19e//nXNnDlTy5Yts7Dq+HegoU2HmjscE3wWicA2Qt8AwJkM0xx0qkYMbdu2TUuWLDlr+6pVq/TYY49p+fLleu+999TS0qLS0lJdffXV+ud//mcVFxeHPYfX65Xb7ZbH4+FjphGQpQIAsItwX78tbWRigUYmfDc/uVPbaxqDTkNONgwtmllIlgoAIKbCff121HdkED1kqQAAnIhGBpLIUgEAOBONDCSRpQIAcCYaGUgiSwUA4Ew0MgggSwUA4DS2/hMFiC2yVAAATkMjg7NMK4x+A+O00D0AgD3RyCCmCN0DAEQS35FBTN35TLW21zQGbdte06i1z7xnUUUAACejkUHMELoHAIg0GhnEDKF7AIBIo5FBzBC6BwCINBoZxAyhewCASKORQUwRugcAiCROv0ZMEboHAIgkGpkEYMfwuWiF7kVirZHeX2Mdz47HDQDshkYmjiVS+Fwk1hrp/TXW8RLpuAHAePEdmTiWSOFzkVhrpPfXWMdLpOMGAONFIxOnEil8LhJrjfT+Gut4iXTcACASaGTiVCKFz0VirZHeX2MdL5GOGwBEAo1MnEqk8LlIrDXS+2us4yXScQOASKCRiVOJFD4XibVGen+NdbxEOm4AEAk0MnEskcLnIrHWSO+vsY6XSMcNAMbLMM1B3yqMM16vV263Wx6PR7m5uVaXY4lECp8b71oPNLTprdpmGZIWTi+ISO7LWGtKpOMGAIOF+/pNIwNofNkt5L4AQOSF+/rNR0uAxpfdQu4LAFiHRgYJbzzZLeS+AIC1aGSQ8MaT3ULuCwBYi0YGCW882S3kvgCAtWhkkPDGk91C7gsAWItGBtD4slvIfQEA63D6NTDAeLJbyH0BgMhxxOnXVVVVuu6661RaWirDMPTiiy8GXW+apr7zne9o0qRJysjI0NKlS7Vv3z5rikVCmFaYpSWzJ2paYZYONLTp9b0nwz7zaOB97Wy06wIAO0uxcvL29nbNnTtXt956q2644Yazrn/ggQf08MMP66c//ammTZumb3/721q2bJk+/PBDpaenW1AxEkG8BtzF67oAJDbbfLRkGIa2bNmi5cuXS+p/N6a0tFRf/epX9bWvfU2S5PF4VFxcrE2bNumzn/1sWOPy0RJG6+Ynd2p7TWNQNkyyYWjRzEI9tXqBhZWNT7yuC0B8csRHS6HU1taqvr5eS5cuDWxzu91auHChduzYMez9fD6fvF5v0AUIV7wG3MXrugDAto1MfX29JKm4uDhoe3FxceC6oaxfv15utztwKSsri2qdiC/xGnAXr+sCANs2MmO1bt06eTyewKWurs7qkuAg8RpwF6/rAgDbNjIlJSWSpBMnTgRtP3HiROC6obhcLuXm5gZdgHDFa8BdvK4LAGzbyEybNk0lJSXaunVrYJvX69Vbb72liooKCytDvIvXgLt4XReAxGbp6ddtbW2qqakJ/FxbW6vq6mrl5+ervLxcd911l773ve9p1qxZgdOvS0tLA2c2AaN1oKFNh5o7QobWuTNT9dTqBVENuAunjkjfPxbrAoBYs/T0623btmnJkiVnbV+1apU2bdok0zR133336d///d/V0tKiyy+/XI8++qjOOeecsOfg9GtI9slQGW8ddlkHAERbuK/ftsmRiRYaGUj2yVAZbx12WQcARJvjc2SASLFLhsp467DLOgDATmhkEPfskqEy3jrssg4AsBMaGcQ9u2SojLcOu6wDAOyERgZxzy4ZKuOtwy7rAAA7oZFBQrBLhsp467DLOgDALjhrCQnFLhkq463DLusAgGjh9OszaGRGFm642nhD3GA/oY4pxxuAlcJ9/bY02RfWCjdcjRC2+BPqmJoyOd4AHIN3ZBJYuOFqhLDFn1DHVBLHG4DlCMRDSOGGqxHCFn9GOqYcbwBOQiOToMINVyOELf6MdExD4XgDsBsamQQVbrgaIWzxZ6RjGgrHG4Dd0MgkqHDD1Qhhiz8jHVOONwAnoZFJYOGGqxHCFn9CHVOONwAn4awlhB2uRghb/Al1TDneAKxEIN4ZNDKRZ0VQmn/Ok95O1Xs7dXF5nj45qyhq8yQbUq/Z/5FKr2metdZY7IMDDW16q7ZZhqSF0ws0rTBLz+08rB21TVo0o1A3faJs1PVEom6C8gDEAo3MGTQykWNFMN5Qc/rlZabql2suV1nB2L+8Gs48fotnFel7yy/Qt178Y1T3QUtHl27/73e140BTyNulJBn671sX6tHK/TEJNSQYEUAs0cicQSMTOVYE4w0150B5mal67ztXR30eqX+tuRkp8p7uieo+uPnJnSEbqqHqikWoIcGIAGKJQDxElBXBeMPNOdCpjm79YRQv+mOdR+pf66mO7qjuA38toxGLUEOCEQHYFY0MwmJFMF64wW3vHj4Vk3lGEol9EKlapMiGGhKMCMCuaGQQFiuC8cINbru4PC8m84wkEvsgUrVIkQ01JBgRgF3RyCAsVgTjDTfnQHmZqeM+eymceaT+teZlpkZ1H/hrGY1YhBoSjAjArmhkEDYrgtKGmtPPf9ZStOfxWzSzUL9cc3nU98GGlfNUMb1gxNulJBl67suXxizUkKA8AHbEWUsYtUgGpYWbSeKfs7HVp2Oe0+PKkRk4p2maQfP750lJMtTTZwb+O7i+wftgpHX0Z8I0STJ06ZlMmFD74UBDm3bWNquhzaeibFcgR+aFd+q0fX9jIEfGf7/h6hw87kjHLpzjEYnjnwhZNImwRiCaOP36DBoZe7JbJs1Y5x9pHS0dXfrK0+/qjf3BmTAV0wv0+Bfmy5R51v3zMlN1qqM7ZF3hzDua/Rur45EIWTSJsEYgFmhkzqCRsSc7ZtKMZf6R1hEqE8b/XZhw8mtGmwsz2v0bq+ORCFk0ibBGIBbIkYFt2TWTZrTzj7SOqj83hMyEqdrXEHZ+zWhyYfzzhrt/Y3U8EiGLJhHWCNgNjQxizs6ZNKOZf6Qx36sbX77NYOHmwow07+D1xep4JEIWTSKsEbAbGhnEnJ0zaUYz/0hjzisbX77NYOHmwow07+D1xep4JEIWTSKsEbAbGhnEnF0zaUY7/0jrWHxOUchMmMWzisLOrxlNLox/3nD3b6yORyJk0STCGgG7oZGBJeyWSTPW+Udax4aV83TZjLMzYSqmF2jDynlD3j9v0JktY8mFGe3+jdXxSIQsmkRYI2AnnLUES0Uyk2Ysc0qKyPwjraO2sV1vHWiSKQ2ZIzP4/uHul3DmHc36YnU8rDjusZYIawSiKS5Ov/7Hf/xHffe73w3aNnv2bO3ZsyfsMWhkglkd0hXO/NGqMRqhdYeaO5RsGOo1zw6jC3f+SK/XymMc7tzD7ctY12z18wHA8MJ9/U6JYU1j8vGPf1yvvfZa4OeUFNuXbEtWh3SFM3+0aoxFaF2oeoeb/3vLP65vvfhBxNZr5TEOd+6hbhdOAKBV9QKwP9u/I/Piiy+qurp6zGPwjkw/q0O6wpk/WjXGIrQuVL3DzZ+bkSLv6Z6IrdfKYxzu3COFEsaqZqufDwBGFjeBePv27VNpaammT5+uz3/+8zp8+HDI2/t8Pnm93qBLorM6pCuc+aNV48jhcScjElo3XL2h5j/V0R2x9Vp5jMOdO5xQwljUbPXzAUBk2bqRWbhwoTZt2qTf/OY3euyxx1RbW6tPfvKTam1tHfY+69evl9vtDlzKyspiWLE9WR3SFc780apx5PC4ljGNO5Jww+tGun+4rDzG4c492n0RrZqtfj4AiCxbNzLXXnutbrrpJs2ZM0fLli3Tr3/9a7W0tOj5558f9j7r1q2Tx+MJXOrq6mJYsT1ZHdIVzvzRqnHk8LgJYxp3JOGG1410/3BZeYzDnXu0+yJaNVv9fAAQWbZuZAabMGGCzjnnHNXU1Ax7G5fLpdzc3KBLorM6pCuc+aNV48jhcRMjElo3XL2h5s/LTI3Yeq08xuHOHU4oYSxqtvr5ACCyHNXItLW1af/+/Zo0aZLVpTiO1SFd4cwfrRpjEVoXqt7h5v/lmssjul4rj3G4c481ADDSrH4+AIgcW5+19LWvfU3XXXedpkyZomPHjum+++5TdXW1PvzwQxUVDf+v6IE4aymY1SFd4cxf29iuNw80arg8l/HO7c99STakXlNBIXRjCa1LSTLU0zd8jow/q2Tw7fy5NU1tXSrIdunS6QUyTXPUuSaDs1BGe4wH3n+084917tEGAI6nxlB1j3TsAFgnLgLxPvvZz6qqqkpNTU0qKirS5Zdfru9///uaMWNG2GPQyDhLNPM9hho70nOMNF9/fswF+sYvdp+VW+POSJHndE/YNY13X4XaHyONFasclvHUGO545McA9hQXjUwk0Mg4SzTzPUJlmEQjQyRUfsxQYXqjrWm8+2qkTJdQY8Uqh2U8NYY7HvkxgD3FTY4MEkc08z1GyjCJdIbISPkx4QhV03j3VTiZLsONFasclvHUOJrxyI8BnI1GBrYRzXyPcDNMIpUhMtb8mKEMVdN499Vo6hs8VqxyWMZT41jGIz8GcCYaGdhGNPM9ws0wiVSGyFjzY4YyVE3j3VejqW/wWLHKYRlPjWMZj/wYwJloZGAb0cz3GCnDJNIZIiPlx4QjVE3j3VfhZLoMN1ascljGU+NoxiM/BnA2GhnYSjTzPYYaO9JzjDSfPz9mqNwad0bKWbcNVdN491Wo/THSWLHKYRlPjeGOR34M4GyctQRbimbeTbj5L5Geb/A8Q+XWjGXd491XA+8vaVRjxSqXaDw1jjQe78QA9sTp12fQyIzd4LAzJ9Ux1tojEbw2eO7KvSdVfaRFF5fn6ZOzioKuP9zUHrjuYxMy9FZtsxpbfZKkohyXFo4QCOgfyx/u5w/7G1jzcPuicu9JvVh9VIakS2cUqLfP1OAQwsG1DyXU+kzT1Mu7j6m5vUvnT8pVUW56YPtbtU2jnm/wekMdm0g+fkczVrSfNyONb5fnLTBe4b5+pwx7DRKWXULDxlLHWGuPRPDaUGP43/UZ7udwXDajQI99fn7Q3CPVK/X/iQXDUFDw3uJZRfr6snP0uZ+8JW/nR+F7v3jvWNB9Ly5za39je1BAX15mqn655nKVFfR/afZQU7uWP7I96HTysaxPki6a7NbBpg61nP5orIHzjSbMMJKP39GMFe3nzUjj2+V5C8Qa78jgLHYJDRtLHWOtPRLBayONMR6LZxUFzT3WuZINQ4ahMTUbUn9z8d53rpYkzfun34adiTNW/vlGE2YYycfvaMaK9vNmpPHt8rwFIoVAPIyJXULDxlLHWGuPRPBaOGOMx8C5xzNXr2mOuYmRpFMd3frDvgZV7j0Z9SbGP9/zbx8OO8wwko/f0YwV7efNSONX/bnBFs9bwAo0Mghil9CwsdQx1tojEbwWyQC8keaOxVyhvHv4lKqPtMRsvjcONI18I/Xvn0g+fkczVrSfNyON/17dqajOD9gZjQyC2CU0bCx1jLX2SASvRTIAb6S5YzFXKBeX5+miyRNiNt9l088+VX0oUwuyIvr4Hc1Y0X7ejDT+vLK8qM4P2BmNDILYJTRsLHWMtfZIBK+FM8Z4DJx7PHMlG4ZSksZeY15mqj45q0j/Z/bEsIP9xiMvM1X/3yXlYYcZRvLxO5qxov28GWn8xecU2eJ5C1iBRgZnsUto2FjqGGvtkQheG2qMwU3DWJqIy2YUnDX3SPVK/WctDQ7e6w/kW6Tc9NAnLF5cNuGsgD7/WUR+v1xz+VnNzFibpHmT3ZqQETzWwPlGE2YYycfvaMaK9vNmpPHt8rwFYo2zljAsu4SGxTIkLhLBa4Pn/sO+Br17+FQgG2Xg9UdOdQSum5yXqTcPNKmprT9HpjB75ByZweF+Q4X8Dbcv/rCvQVvePSrJVMWMQvX2mUHhfP7bDKx9KKHWJ0n/s/uYGtt8+nipW4U5rsD2wWGA4cw3mjDDSD5+RzNWtJ83I41vl+ctMF4E4p1BI2OdWARzDTVHrAPBhpqvcu9Jvb73pJIMQ3lZaSoKoykJFZQ3OKRvcCDceNY8+L7P7TysHbVNmlWUo/zsNBlSoPZw5vkotG7ocL7/2nFIdac6dHF5ns4rzT0rfLA/LK/5rHn920onpA8bhuff74YMzS7J0cLpBUMGGw6ucTQBe8mGoaMtHRoY5heJx6EdHstWSIQ1YmwIxINlYhHMNdQcl80okGlKOw4EB8BFKxBsqBoWTM3TnnqvvJ29Q95nrOF2oeRlpgadCh3umkc774SM1KDAunCC6PwqpheovatHu494Atte+9PJsObNTU8JCu8byF9Dy+kuXf/j7UH1DWWox8hwY4azLklyZ6TKM2De0T4Ohxp/uDDDeAq3I8APkcI7Moi4WARzhRsIF81AsLGG0kUq3G444a55vPOGE0QXbf4a3j/aErFcm2isK9QxscNj2QoE+GEkBOLBErEI1BtNIFy0AsHGE0oXqXC74YSz5kjMG04QXbT5a4hkOF801jXcMbHDY9kKdgneRHygkUFExSJQbyyBcJEOBBtvKF0swu1CrTmS84YTROdE0VjX4GNih8eyFewSvIn4wHdkEFGxCNQbSyBcpAPBxhtKF4twu1BrjuS8/i/oxptorGvwMbHDY9kKdgneRHzgHRlEVCwC9UYTCBetQLDxhNJFKtxuOOGsORLzhhNEF23+GiIZzheNdQ13TOzwWLaCXYI3ER9oZBBxsQjmGmqOy2YUqGL62QFw0QoEG6qGBVPzlZuePOx9xhpuF8rgF/Fw1zzaeQcH1oUTROdXMb1Acye7w55roFDhff4afrnm8rPqG8pQj5HhxvQbaT+5B8072sfhUOMPF2YYT+F2BPghUjhrCVETi2CuoeYYvC1UrolfqCyLyr0nA5ku/pA2f15JQbZLF5VNOCuc7Q/7GrT1TyeUbBjqkyRTOm9Sjopy0/XBUY9qGtq0aEah5k/J06HmDp30duqtM6frXjqjQHuOt+pwc4eWfbxEn5iaf1ZI38BAuMNN7Xqiar/qPZ2aV56n6y/6WGCdAzNUBmfT+HNZkgzpw+NeJRmGZhXn6FR7l377Qb28nT1yZ6Sos7tXkqHr5pbq2gsn6c0DTWps7VRhTnpQmN3AY/Kr/3tMBxvbNSEzRbOKc5VkSL/78ISOnTotX0+vunr7lJacpOLcDJ03KVfNHV1q6+xWu69H9d5OGaYhX1+fDJkqz89SZlqy2rt6VJKbofzMVB1s7tBpX69MSZMmuJSbkSpPR4/afD1q8PpkGNLUwiz99cUf07GW09pxoEmLZhTqpk+USZKq/tyg3+85obbOHuWkpyrZkA42d6g8P1NXzJ6o1/ec1JFTHbqobILys11qbO1UbWO7vJ09urg8T719fWpq79JV5xUH9mVjq0+FOS5deia75uXdx3WwsV3Ti7JUlOPSn463ynu6W7kZqbrqvImBfCD/4/LDox7ta2jTORNzAvk6kvTmgUYNzKwZ/Lh9q7Zp2OsHG5jHM1yu0cAxPzYhQ72mqZPeTtV7O4OeA4OfM6PJGBr42Awn3HCo+ydS1k6iIhDvDBqZxBUq/8OfV2HKHDbLouV0l5Y/sj3ojJjc9BSZpqlWX3BOzCVT8vSTVZdENB/Gr2J6gR7/wvyzsjUONbXrLx7+g9p8Q2fWxMrA+g41tZ+1z+xkQkaqphVk6r0BeTZ2NlxGkClTX3n63aCcGWn4x0pLR5du/+93z8rPGZhr1NLRNeSYg7kzUjS7OFc7DzaPWGe4WTwj5cckatZOoqOROYNGJnGFyufw51VIGjbLYrTZJNHMhxk8tiTN+6ff2qZh8Ndnp5ri0cDHbaimYPBj5eYnd454+1C3GWud4WbxjJQfk6hZO4mOZF8kNH9OxXD8eRWjvS4Uf/6F/232SL0oDB5b6v9oy04NQ9W+Bj339mFb1RSPwnlsDn6sjPRYrNrXoKo/n4zo43VgHowZRs0Dbz/Ux2fh1hZqHMQvvuyLuGRVrkk082EGZmtUH2mJ+PjjFSr2H7E18LESzmPxvbqWqNUxmufCUPkxiZq1g/DRyCAuRTOfJZRo5sMMzNa4aPKEiI8/XiOdDYTYGfhYCeexOK9sQtTqGM1zYaj8mETN2kH4aGQQl0bK5/DnVYTKshhtNkk082EGZ2v8n9kTI5qdMl6LZxXpM5eU26qmeDTwcTucwY8V/2Mx1O0XnzMx5G3GWme4WTyh8mMSNWsH4aORQdwKlf/hz6sIlWXxyzWXn/XCnJueohzX2Tkxl0zJi3g+jF/F9LOzZyTpl2suH7KWWBtY31D7zE4mZKRq3hjzbKwwXEbQhpXzzsqZkYZ/rGxYOW/Id8wG5hoNN+Zg7owULZiaH1adA+cP9VwYKT8mUbN2EB7OWkLc8+fKhMqrCJV584d9DXr38KmgDI0/7GvQ7/90UvnZafrLOaUh/wU4eP4/HfPqzydbtWhGYSAjprHVpw+PeZWfnaaLyibo6KnTMqWwskH+sK9B//XGITW09efILDl3YmCdkgLrOnKqI7COyXmZevNAkwxJH8vL0JFTp9XU5lNBdn8OyjsHm/XbD+qVm5Gqts4enWjt1NXnl+grS2aqtrFdbx1oClmff5+VujNUmONSY6tPv/vghBrafJpRlKXcjFQdaT6t7PQUzZyYrZOtnUo2DB05dVqHmztUkJWmjq5etZzu1iR3uopz05WXmaqZxTlKSTL0wTGvjrecVmtnj0onpMud0X/6b15mqvKy0iRJhdkuLZxeELRu//GrbWzXmwea1NTmU2G2S8lJRtD+f33PSdU1d2heeZ7ystLU1ObTqfYu9ZrSVedN1OS8zJDZPpIC4xdkuwI1t3V2Kzs9NWgM//0aW3065jkd2Gf+x2Kox2Y4x2Lw7f3HfbgcmYFjTs7LCKpt8D4cWFc4uVEDbzPwsRnuOyjh5EYhfsTV6dePPPKIfvCDH6i+vl5z587Vhg0btGBBeKfXRauRcWoIUyRCrAaGf31sQrp6TYW9H4ab70BDm17efUy1jR3Ky0zVrOKcQLDYcPU9t/OwdtQGB535w7wa27okUyrK6X8x66/bH/KVrqMtnWps9akox6XSM6Ff/jn8454zMUd5Wan6c32bPKe7Na0oa8imxV/7mweaJNPQ+ZNy1GP2v905qyQn8AJzoKFNP/79PlXXtaj1dLd6+kzlpKcq05WsKfmZKp2Qob0nWnWgoV09vX06tyRXfzm3VPtOtGr7/kYlyZArtf9N1GmF2ZpWmKX3j3h0sq1Tn5iSr+Z2nz487lV5XpYuLHPr9386oT+faFNfn6ksV7LK8zNVmO1Sn6RlHy/RTZ8o04+37tMv3j2i9NRkTc7P1MnWTk3MTld7V7feP+JRT5+pgmyXLpmar7cPNuvImQbLkJRsSEaS1NcnpSUnKcuVrJTkJGWlJaurt09NbT719kmulCR195oyDFOGYSjZSFJyUv8ZJsmGoR7TVE+vqSRDyknv/5d96+luZbhSlJwktXX2KiM1Se6MVLV39aqnt0++nj6d7upVnykZhpSemqSc9FS1dnaru6dPJe4MXTI1Xx8e9+i4p1O+nl75uvt/1WWfeRer9UygXkqSlJmarOQkQ8lJhrp6+pSc1P8xRntXj/r6pExXiiZkpKqxrVOd3aayXcm6bEaBSiZk6tip02rr6lZ2WqrafN2q93bKe7pbHV09SklKkjsjTW2+biUnJ2l2cY5M01Rdc4c8p7vl6+lTkmEoNdlQWkqypuRnqM+UGtt9MiX19JrKy0xVeUGmdh08pdbObhmGlJaSLHdGqrp7+uTt7FZSkqHinHSlpST1jympras/XygzLVmzJuao1detuuYOTcnP0sIZ+apt7FDr6W5Nyc9Sj9mnYy2n1dDqkyQZpqH0tGQlGf21XPixCSrLz1Bze7eKctJUc7JN7Z29mlPm1ps1jfrTyVZNzcvS4nOLNMmdrqKcdCUbUnVdi3Yf8Sg3PVUzirLU0OZTkmH0v4NjGGpu82lvfatkSOeX5ConM1WTcvt/pxiSGlt92r6/UYYhXTajUHMmu3W0pVOGpNIBv3v8z+8/n2iVaRq66ryJ+uSsorPCLD8Kx5R6TQXCCAeGZfp/x4Vq9qIp1O/ksfz+trLe8YqbRua5557TzTffrMcff1wLFy7UQw89pBdeeEF79+7VxIkTR7x/pBuZoYKZnBDCNFTd4w2xGizUfhhuv31v+QX66gvVevvgqbDHP9zcrr9+9A319H300E1J6m8gdh/1hlXrePjD70yZ+pun3tHOMGrPSk1Se3df1GsDMDT/O1+jNTA0MJpCvbYMFdw50u/vaIvFa2HcNDILFy7UJZdcoh//+MeSpL6+PpWVlWnt2rX6xje+MeL9I93IDBXM5IQQpnACpUYbYjXS/UeaP9kwlJuREnb2iH/8N/Y3jukXUiT5vxgZyewNAPY0VMhgpIV6bZHODu4cLNavQ7F4LQz39dvWX/bt6urSrl27tHTp0sC2pKQkLV26VDt27BjyPj6fT16vN+gSKf5gpsEPpoEhTHY0XN2DDVxHuPcZ7v7hzN9rmqMKUPOPb3UTI50JEaOJARJCtH+/j/TaMtrf39Fmt9dCWzcyjY2N6u3tVXFxcdD24uJi1dfXD3mf9evXy+12By5lZWURq2ekYCa7hjCNNlBqtCFWQ91/PPMDgN1E8/d7JH9HxuJ1yG6vhbZuZMZi3bp18ng8gUtdXV3Exh4pmMmuIUyjDZQabYjVUPcfz/wAYDfR/P0eyd+RsXgdsttroa0bmcLCQiUnJ+vEiRNB20+cOKGSkpIh7+NyuZSbmxt0iZThgpnsHsIUbqDUaEOsQt0/nPmT/WcvjHL8lKTIhMyNx0ihZADiR7R/v4/02jLa39/RZrfXQls3MmlpaZo/f762bt0a2NbX16etW7eqoqLCkppCBajZ2VB1jzfEarBQ+2G4/fbLNZdrwdS8UY3/yzWLzmpmUpIMzflYbHKC/OF3G1bOOysYbDhZqbZ+qgFxb6z/ABoYGhhNoV5bxvL7O9rs9Fpo+7OWnnvuOa1atUpPPPGEFixYoIceekjPP/+89uzZc9Z3Z4YSrRwZp4YwRSLEamDQlz8wK9z9MNx8tY3t+p/dx3SgoT0QfHbpmSTS4ep74Z06bd/fGJQj4w/zamjzydBHoWiSgkK+/AFwhdkufWzQGvzjnjMxR/lZafrziVa1dAyfI+Ovfcf+JpmmdH5prnr6TCUbCqzDv68f+f0+vXf4lFpPd6v7TI5Mf8ZLlibnZWhPvVf7B+XI1Jxs0/aaBhkylJ6aLFPS9KKsQI7Midb+HJmWji59cMyj8vwsXTDZrW17TmpPfWtQjkxRTrp6TTOQI/Po6zX6+a46uVKSVV6QqRPeThXlpKvD16PdR1vU09ufI7NgWr7eru3PkenTRzkySWdyZFLDypFRf45MknEmx+Oj/I7uXlNJhqGc9BQZkrynu5V5JkemNZAjk6b2rh719Papq6dPHUE5MsnKSU9Ra2ePunt6VeLO0IJp+frwmFfHPafV2dMn35nT30PlyKQkGfINyJHp6OpR75kcmbyMVDUMzJGZWagSd0Z/MJ+vWzmuMzkynk55BuTITMhMU2tnf47MucU5MiUdbu6Qp6MrOEcmNVlT8jJlylRDW5ek/v2Sl5mqKQVZeudg84AcmRS5M1LU3WvKe7rrTI5MhlwpSers6VWSDLV19chQf47MzIk5avf16FBzu6bkZ+nSGQWqbWyX93S3phRkqbfP1LGW0zrZ2p/RojOPteQkqaGtP0dmSkGmGtu6NDHHpZqTbWrr7NGcMrfeOtCkD+u9mpafpU/OLgoE+qUkGfq/dS3aXdei7PTUoODDvKw0mZJa2ru0p75VpkydP8mtnIwUlbr7n4+GpOb2Lv3vvgYZhqGKGQWaWzZBR06dDgQ5DgwgfOtAk/adaA2EFn5yVtFZYZaDwykH/3fg7zircmRGCkAc7e9vK+sdr7g5/VqSfvzjHwcC8S666CI9/PDDWrhwYVj3JdkXAADniatGZjxoZAAAcJ64yJEBAAAIhUYGAAA4Fo0MAABwLBoZAADgWDQyAADAsWhkAACAY9HIAAAAx6KRAQAAjkUjAwAAHCvF6gKizR9c7PV6La4EAACEy/+6PdIfIIj7Rqa1tVWSVFZWZnElAABgtFpbW+V2u4e9Pu7/1lJfX5+OHTumnJwcGcbY/oy7HXi9XpWVlamuri4h/2ZUIq8/kdcusX7Wz/oTdf2maaq1tVWlpaVKShr+mzBx/45MUlKSJk+ebHUZEZObm5twD+aBEnn9ibx2ifWzftafiOsP9U6MH1/2BQAAjkUjAwAAHItGxiFcLpfuu+8+uVwuq0uxRCKvP5HXLrF+1s/6E3n94Yj7L/sCAID4xTsyAADAsWhkAACAY9HIAAAAx6KRAQAAjkUj40B/9Vd/pfLycqWnp2vSpEn64he/qGPHjlldVkwcPHhQq1ev1rRp05SRkaEZM2bovvvuU1dXl9Wlxcz3v/99XXbZZcrMzNSECROsLifqHnnkEU2dOlXp6elauHChdu7caXVJMVFVVaXrrrtOpaWlMgxDL774otUlxdT69et1ySWXKCcnRxMnTtTy5cu1d+9eq8uKiccee0xz5swJhOBVVFTolVdesbos26KRcaAlS5bo+eef1969e/Xzn/9c+/fv14033mh1WTGxZ88e9fX16YknntAHH3ygBx98UI8//ri++c1vWl1azHR1demmm27S7bffbnUpUffcc8/pnnvu0X333ad3331Xc+fO1bJly3Ty5EmrS4u69vZ2zZ07V4888ojVpViisrJSa9as0Ztvvqnf/e536u7u1tVXX6329narS4u6yZMn6/7779euXbv0zjvv6Morr9T111+vDz74wOrS7MmE47300kumYRhmV1eX1aVY4oEHHjCnTZtmdRkxt3HjRtPtdltdRlQtWLDAXLNmTeDn3t5es7S01Fy/fr2FVcWeJHPLli1Wl2GpkydPmpLMyspKq0uxRF5envmTn/zE6jJsiXdkHK65uVlPP/20LrvsMqWmplpdjiU8Ho/y8/OtLgMR1tXVpV27dmnp0qWBbUlJSVq6dKl27NhhYWWwgsfjkaSEe6739vbq2WefVXt7uyoqKqwux5ZoZBzq3nvvVVZWlgoKCnT48GG99NJLVpdkiZqaGm3YsEF/+7d/a3UpiLDGxkb19vaquLg4aHtxcbHq6+stqgpW6Ovr01133aVFixbpggsusLqcmHj//feVnZ0tl8ul2267TVu2bNH5559vdVm2RCNjE9/4xjdkGEbIy549ewK3//u//3u99957+u1vf6vk5GTdfPPNMh0c0jza9UvS0aNHdc011+imm27Sl7/8ZYsqj4yxrB9IFGvWrNEf//hHPfvss1aXEjOzZ89WdXW13nrrLd1+++1atWqVPvzwQ6vLsiX+RIFNNDQ0qKmpKeRtpk+frrS0tLO2HzlyRGVlZXrjjTcc+9bjaNd/7NgxXXHFFbr00ku1adMmJSU5uycfy/HftGmT7rrrLrW0tES5Omt0dXUpMzNTP/vZz7R8+fLA9lWrVqmlpSWh3oU0DENbtmwJ2g+J4o477tBLL72kqqoqTZs2zepyLLN06VLNmDFDTzzxhNWl2E6K1QWgX1FRkYqKisZ0376+PkmSz+eLZEkxNZr1Hz16VEuWLNH8+fO1ceNGxzcx0viOf7xKS0vT/PnztXXr1sALeF9fn7Zu3ao77rjD2uIQdaZpau3atdqyZYu2bduW0E2M1P/Yd/Lv+GiikXGYt956S2+//bYuv/xy5eXlaf/+/fr2t7+tGTNmOPbdmNE4evSorrjiCk2ZMkU//OEP1dDQELiupKTEwspi5/Dhw2pubtbhw4fV29ur6upqSdLMmTOVnZ1tbXERds8992jVqlX6xCc+oQULFuihhx5Se3u7vvSlL1ldWtS1tbWppqYm8HNtba2qq6uVn5+v8vJyCyuLjTVr1mjz5s166aWXlJOTE/helNvtVkZGhsXVRde6det07bXXqry8XK2trdq8ebO2bdumV1991erS7Mnak6YwWrt37zaXLFli5ufnmy6Xy5w6dap52223mUeOHLG6tJjYuHGjKWnIS6JYtWrVkOt//fXXrS4tKjZs2GCWl5ebaWlp5oIFC8w333zT6pJi4vXXXx/yOK9atcrq0mJiuOf5xo0brS4t6m699VZzypQpZlpamllUVGReddVV5m9/+1ury7ItviMDAAAcy/lfLgAAAAmLRgYAADgWjQwAAHAsGhkAAOBYNDIAAMCxaGQAAIBj0cgAAADHopEBAACORSMDAAAci0YGgG319vbqsssu0w033BC03ePxqKysTP/wD/8gSbrzzjs1f/58uVwuXXTRRRZUCsAqNDIAbCs5OVmbNm3Sb37zGz399NOB7WvXrlV+fr7uu+++wLZbb71Vn/nMZ6woE4CF+OvXAGztnHPO0f3336+1a9fqyiuv1M6dO/Xss8/q7bffVlpamiTp4YcfliQ1NDRo9+7dVpYLIMZoZADY3tq1a7VlyxZ98Ytf1Pvvv6/vfOc7mjt3rtVlAbABGhkAtmcYhh577DGdd955uvDCC/WNb3zD6pIA2ATfkQHgCP/5n/+pzMxM1dbW6siRI1aXA8AmaGQA2N4bb7yhBx98UC+//LIWLFig1atXyzRNq8sCYAM0MgBsraOjQ7fccotuv/12LVmyRE8++aR27typxx9/3OrSANgAjQwAW1u3bp1M09T9998vSZo6dap++MMf6utf/7oOHjwoSaqpqVF1dbXq6+t1+vRpVVdXq7q6Wl1dXRZWDiAWDJP3ZwHYVGVlpa666ipt27ZNl19+edB1y5YtU09Pj1577TUtWbJElZWVZ92/trZWU6dOjVG1AKxAIwMAAByLj5YAAIBj0cgAAADHopEBAACORSMDAAAci0YGAAA4Fo0MAABwLBoZAADgWDQyAADAsWhkAACAY9HIAAAAx6KRAQAAjvX/ANsS/dJoEt6QAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "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": 83,
   "id": "b46e4e07-d32c-4712-a379-ac2512403b9b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Axes: xlabel='X1', ylabel='preds'>"
      ]
     },
     "execution_count": 83,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjIAAAGyCAYAAAD+lC4cAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAOwJJREFUeJzt3X90FPW9//HXZJNsyK+FJEtCaiD8iL+gIKJgwOZKoSK9x1ta9LZcb0HLaSvlR5Vzb5W21tofN962p1opRVss1HMFbG3R21p/VCrkWwRUNEVr4RISIICBTUJ2kw3ZJJv5/oHZspAfm7DZ2Umej3P2yH7mM595z+xm5+XszKxhmqYpAAAAG0qwugAAAID+IsgAAADbIsgAAADbIsgAAADbIsgAAADbIsgAAADbIsgAAADbIsgAAADbIsgAAADbSrS6gIHW0dGhkydPKiMjQ4ZhWF0OAACIgGmaamxsVH5+vhISejjuYsaJ0tJSU5L51a9+NdR29uxZ8ytf+YqZlZVlpqWlmZ/5zGfMmpqaPo1bXV1tSuLBgwcPHjx42PBRXV3d434+Lo7IvPnmm3riiSc0efLksPZ7771XL7zwgn7zm9/I5XJpxYoV+sxnPqNdu3ZFPHZGRoYkqbq6WpmZmVGtGwAADAyfz6eCgoLQfrw7lgeZpqYm3XHHHfrFL36h733ve6F2r9erJ598Ups3b9bHP/5xSdLGjRt11VVXac+ePbrhhhsiGr/z66TMzEyCDAAANtPbaSGWn+y7fPly/fM//7Pmzp0b1r5v3z61tbWFtV955ZUaPXq0du/e3e14gUBAPp8v7AEAAAYnS4/IbN26VW+//bbefPPNi6bV1NQoOTlZw4cPD2vPzc1VTU1Nt2OWlpbqoYceinapAAAgDll2RKa6ulpf/epX9fTTTyslJSVq465Zs0Zerzf0qK6ujtrYAAAgvlgWZPbt26fTp0/r2muvVWJiohITE7Vz50499thjSkxMVG5urlpbW9XQ0BA236lTp5SXl9ftuE6nM3Q+DOfFAAAwuFn21dKcOXP07rvvhrXddddduvLKK3XfffepoKBASUlJ2r59uxYuXChJOnjwoI4dO6bi4mIrSgYAAHHGsiCTkZGhSZMmhbWlpaUpOzs71L506VKtXr1aWVlZyszM1MqVK1VcXBzxFUsAAGBws/zy65488sgjSkhI0MKFCxUIBDRv3jz97Gc/s7osAAAQJwzTNE2rixhIPp9PLpdLXq+X82UAALCJSPfflt9HBgAAoL8IMgAAwLbi+hwZDLxKT5OO1jerMDtNY3PSrC4HAIA+IcgMUQ3NrVq1pVxlhzyhtpIit9YumipXapKFlQEAEDm+WhqiVm0p166K2rC2XRW1WrnlHYsqAgCg7wgyQ1Clp0llhzwKXnDBWtA0VXbIo6pav0WVAQDQNwSZIehofXOP04/UEWQAAPZAkBmCxmSl9ji9MJuTfgEA9kCQGYLGudNVUuSWwzDC2h2GoZIiN1cvAQBsgyAzRK1dNFWzJuSEtc2akKO1i6ZaVBEAAH3H5ddDlCs1SU8tna6qWr+O1Pm5jwwAwJYIMkPc2BwCDADAvvhqCQAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2JalQWb9+vWaPHmyMjMzlZmZqeLiYr344ouh6TfddJMMwwh73H333RZWDAAA4kmilQu/7LLL9PDDD6uoqEimaepXv/qVPvWpT+mdd97RxIkTJUlf/OIX9Z3vfCc0T2pqqlXlAgCAOGNpkLn11lvDnn//+9/X+vXrtWfPnlCQSU1NVV5eXsRjBgIBBQKB0HOfzxedYgEAQNyJm3NkgsGgtm7dKr/fr+Li4lD7008/rZycHE2aNElr1qxRc3Nzj+OUlpbK5XKFHgUFBQNduq1Uepr02sHTqqr1X1IfAADigWGapmllAe+++66Ki4vV0tKi9PR0bd68WZ/85CclST//+c81ZswY5efna//+/brvvvs0ffp0/e53v+t2vK6OyBQUFMjr9SozM3PA1ydeNTS3atWWcpUd8oTaSorcWrtoqlypSRH3AQAgFnw+n1wuV6/7b8uDTGtrq44dOyav16tnn31WGzZs0M6dO3X11Vdf1PfPf/6z5syZo4qKCo0fPz6i8SPdEIPd4iff0K6KWgXPe7kdhqFZE3L01NLpEfcBACAWIt1/W/7VUnJysiZMmKBp06aptLRUU6ZM0U9+8pMu+86YMUOSVFFREcsSba/S06SyQ56wgCJJQdNU2SGPqmr9EfUBACDeWB5kLtTR0RH21dD5ysvLJUmjRo2KYUX2d7S+5/OKjtT5I+oDAEC8sfSqpTVr1mj+/PkaPXq0GhsbtXnzZu3YsUMvv/yyDh8+HDpfJjs7W/v379e9996rkpISTZ482cqybWdMVs+XrBdmp6m3bxgLs9OiWRIAAFFhaZA5ffq0Fi9erA8++EAul0uTJ0/Wyy+/rE984hOqrq7Wq6++qkcffVR+v18FBQVauHChvvnNb1pZsi2Nc6erpMjd7fkvY3POhZRI+gAAEE8sP9l3oHGy7zne5jat3PJOj1ckRdIHAIBYsM1VSwONIBOuqtavI3V+FWandXuUJZI+AAAMpEj335Z+tYTYG5vTeziJpA8AAPEg7q5aAgAAiBRBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2Fai1QXA3io9TTpa36zC7DSNzUkbdMsDAMQ3ggz6paG5Vau2lKvskCfUVlLk1tpFU+VKTbL98gAA9sBXS+iXVVvKtauiNqxtV0WtVm55Z1AsDwBgDwQZ9Fmlp0llhzwKmmZYe9A0VXbIo6pav62XBwCwD4IM+uxofXOP04/URTdYxHp5AAD7IMigz8ZkpfY4vTA7uifhxnp5AAD7IMigz8a501VS5JbDMMLaHYahkiJ31K8mivXyAAD2YWmQWb9+vSZPnqzMzExlZmaquLhYL774Ymh6S0uLli9fruzsbKWnp2vhwoU6deqUhRWj09pFUzVrQk5Y26wJOVq7aOqgWB4AwB4M07zgDMoY+v3vfy+Hw6GioiKZpqlf/epX+uEPf6h33nlHEydO1LJly/TCCy9o06ZNcrlcWrFihRISErRr166Il+Hz+eRyueT1epWZmTmAazM0VdX6daTOH7P7usR6eQAAa0S6/7Y0yHQlKytLP/zhD3XbbbfJ7XZr8+bNuu222yRJBw4c0FVXXaXdu3frhhtuiGg8ggwAAPYT6f47bs6RCQaD2rp1q/x+v4qLi7Vv3z61tbVp7ty5oT5XXnmlRo8erd27d3c7TiAQkM/nC3sAAIDByfIg8+677yo9PV1Op1N33323tm3bpquvvlo1NTVKTk7W8OHDw/rn5uaqpqam2/FKS0vlcrlCj4KCggFeAwAAYBXLg8wVV1yh8vJy7d27V8uWLdOSJUv0/vvv93u8NWvWyOv1hh7V1dVRrBYAAMQTy39rKTk5WRMmTJAkTZs2TW+++aZ+8pOf6LOf/axaW1vV0NAQdlTm1KlTysvL63Y8p9Mpp9M50GUDAIA4YPkRmQt1dHQoEAho2rRpSkpK0vbt20PTDh48qGPHjqm4uNjCCgEAQLyw9IjMmjVrNH/+fI0ePVqNjY3avHmzduzYoZdfflkul0tLly7V6tWrlZWVpczMTK1cuVLFxcURX7EEAAAGN0uDzOnTp7V48WJ98MEHcrlcmjx5sl5++WV94hOfkCQ98sgjSkhI0MKFCxUIBDRv3jz97Gc/s7JkAAAQR+LuPjLRxn1kAACwH9vdRwYAAKCvCDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2CDIAAMC2LA0ypaWluv7665WRkaGRI0dqwYIFOnjwYFifm266SYZhhD3uvvtuiyoGAADxxNIgs3PnTi1fvlx79uzRn/70J7W1tenmm2+W3+8P6/fFL35RH3zwQejxgx/8wKKKAQBAPEm0cuEvvfRS2PNNmzZp5MiR2rdvn0pKSkLtqampysvLi3V5AAAgzsXVOTJer1eSlJWVFdb+9NNPKycnR5MmTdKaNWvU3Nzc7RiBQEA+ny/sAQAABidLj8icr6OjQ/fcc49mzZqlSZMmhdr/7d/+TWPGjFF+fr7279+v++67TwcPHtTvfve7LscpLS3VQw89FKuyAQCAhQzTNE2ri5CkZcuW6cUXX9Rf/vIXXXbZZd32+/Of/6w5c+aooqJC48ePv2h6IBBQIBAIPff5fCooKJDX61VmZuaA1A4AAKLL5/PJ5XL1uv+OiyMyK1as0B/+8AeVlZX1GGIkacaMGZLUbZBxOp1yOp0DUicAAIgvlgYZ0zS1cuVKbdu2TTt27NDYsWN7nae8vFySNGrUqAGuDgAAxDtLg8zy5cu1efNmPf/888rIyFBNTY0kyeVyadiwYTp8+LA2b96sT37yk8rOztb+/ft17733qqSkRJMnT7aydAAAEAcsPUfGMIwu2zdu3Kg777xT1dXV+vd//3e999578vv9Kigo0Kc//Wl985vfjPh8l0i/YwMAAPHDFufI9JahCgoKtHPnzhhVAwAA7Cau7iMDAADQFwQZAABgWwQZAABgWwQZAABgWwQZAABgW3FxZ1/ETqWnSUfrm1WYnaaxOWlR6xuLemIlHmsCAHSNIDNENDS3atWWcpUd8oTaSorcWrtoqlypSf3uG4t6YiUeawIA9IyvloaIVVvKtauiNqxtV0WtVm5555L6xqKeWInHmgAAPSPIDAGVniaVHfIoeMENCIOmqbJDHlXV+vvVNxb1xEo81gQA6B1BZgg4Wt/c4/Qjdf/YSfelbyzqiZV4rAkA0DuCzBAwJiu1x+mF2f84obUvfWNRT6zEY00AgN4RZIaAce50lRS55bjgRzodhqGSInfYlTl96RuLemIlHmsCAPSOIDNErF00VbMm5IS1zZqQo7WLpl5S31jUEyvxWBMAoGeG2dtPUNtcpD8DPlRU1fp1pM4f0T1S+tI3FvXESjzWBABDTaT7b4IMAACIO5Huv/lqCQAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2FZUgkwwGFR5ebnOnDkTjeEAAAAi0q8gc8899+jJJ5+UdC7E/NM//ZOuvfZaFRQUaMeOHdGsDwAAoFv9CjLPPvuspkyZIkn6/e9/r6qqKh04cED33nuvvvGNb0S1QAAAgO70K8jU1tYqLy9PkvTHP/5Rt99+uy6//HJ94Qtf0LvvvhvVAgEAALrTryCTm5ur999/X8FgUC+99JI+8YlPSJKam5vlcDiiWiAAAEB3Evsz01133aV//dd/1ahRo2QYhubOnStJ2rt3r6688sqoFggAANCdfh2R+fa3v60NGzboS1/6knbt2iWn0ylJcjgcuv/++yMep7S0VNdff70yMjI0cuRILViwQAcPHgzr09LSouXLlys7O1vp6elauHChTp061Z+yAQDAIGOYpmlatfBbbrlFn/vc53T99dervb1dX//61/Xee+/p/fffV1pamiRp2bJleuGFF7Rp0ya5XC6tWLFCCQkJ2rVrV0TL8Pl8crlc8nq9yszMHMjVAQAAURLp/jviIPPYY49FvPBVq1ZF3Pd8Ho9HI0eO1M6dO1VSUiKv1yu3263NmzfrtttukyQdOHBAV111lXbv3q0bbrjhojECgYACgUDouc/nU0FBAUEGAAAbiTTIRHyOzCOPPBL23OPxqLm5WcOHD5ckNTQ0KDU1VSNHjux3kPF6vZKkrKwsSdK+ffvU1tYWOgdHkq688kqNHj262yBTWlqqhx56qF/LBwAA9hLxOTJVVVWhx/e//31dc801+vvf/676+nrV19fr73//u6699lp997vf7VchHR0duueeezRr1ixNmjRJklRTU6Pk5ORQWOqUm5urmpqaLsdZs2aNvF5v6FFdXd2vegAAQPzr11VLDzzwgJ599lldccUVobYrrrhCjzzyiG677TbdcccdfR5z+fLleu+99/SXv/ylPyWFOJ3O0MnHAABgcOvXVUsffPCB2tvbL2oPBoP9uqJoxYoV+sMf/qDXXntNl112Wag9Ly9Pra2tamhoCOt/6tSp0A35cOkqPU167eBpVdX6bTk+AGDo6tcRmTlz5ujLX/6yNmzYoGuvvVbSufNZli1bFnY+S29M09TKlSu1bds27dixQ2PHjg2bPm3aNCUlJWn79u1auHChJOngwYM6duyYiouL+1M6ztPQ3KpVW8pVdsgTaispcmvtoqlypSbF/fgAAPTriMwvf/lL5eXl6brrrgt9lTN9+nTl5uZqw4YNEY+zfPly/c///I82b96sjIwM1dTUqKamRmfPnpUkuVwuLV26VKtXr9Zrr72mffv26a677lJxcXGXJ/qib1ZtKdeuitqwtl0VtVq55R1bjA8AQL+OyLjdbv3xj3/U//3f/+nAgQOSzl1NdPnll/dpnPXr10uSbrrpprD2jRs36s4775R07mqphIQELVy4UIFAQPPmzdPPfvaz/pSN81R6msKOlHQKmqbKDnlUVevX2Jy0uB0fAACpn0GmU2FhoUzT1Pjx45WY2PehIrmFTUpKitatW6d169b1p0R042h9c4/Tj9RdWtAY6PEBAJD6+dVSc3Ozli5dqtTUVE2cOFHHjh2TJK1cuVIPP/xwVAvEwBiTldrj9MLsSwsZAz0+AABSP4PMmjVr9Ne//lU7duxQSkpKqH3u3Ll65plnolYcBs44d7pKitxyGEZYu8MwVFLkvuSjJQM9PgAAUj+DzHPPPaef/vSnuvHGG2Wct6OaOHGiDh8+HLXiMLDWLpqqWRNywtpmTcjR2kVTbTE+AAD9Okem8zeRLuT3+8OCDeKbKzVJTy2drqpav47U+VWYnRbVIyUDPT4AAP06InPdddfphRdeCD3vDC8bNmzg/i42NDYnTbOvGDlgIWOgxwcADF39OiLzX//1X5o/f77ef/99tbe36yc/+Ynef/99vf7669q5c2e0awQAAOhSv47I3HjjjfrrX/+q9vZ2ffSjH9Urr7yikSNHavfu3Zo2bVq0awQAAOhSn4/ItLW16ctf/rIeeOAB/eIXvxiImgAAACLS5yMySUlJ+u1vfzsQtQAAAPRJv75aWrBggZ577rkolwIAANA3/TrZt6ioSN/5zne0a9cuTZs2TWlp4VejrFq1KirFAQAA9MQwI/nBowuMHTu2+wENQ5WVlZdUVDT5fD65XC55vV5lZmZaXQ4AAIhApPvvfh2RqaqqCv27MwdxIzwAABBr/TpHRpKefPJJTZo0SSkpKUpJSdGkSZO0YcOGaNYGAADQo34dkfnWt76lH//4x1q5cmXoTr67d+/Wvffeq2PHjuk73/lOVIsEAADoSr/OkXG73Xrssce0aNGisPYtW7Zo5cqVqq2tjVqBl4pzZAAAsJ9I99/9+mqpra1N11133UXt06ZNU3t7e3+GBAAA6LN+BZnPf/7zWr9+/UXtP//5z3XHHXdcclEAAACR6Nc5MtK5k31feeUV3XDDDZKkvXv36tixY1q8eLFWr14d6vfjH//40qsEAADoQr+CzHvvvadrr71WknT48GFJUk5OjnJycvTee++F+nFJNgAAGEj9CjKvvfZatOsAAADos37fRwYAAMBqBBkAAGBbBBkAAGBbBBkAAGBbBBkAAGBbBBkAAGBbBBkAAGBbBBkAAGBb/f6JAgxtlZ4mHa1vVmF2msbmpMVdDZWeJu2tqpchaca4bMtqjLZ42O4AEE8IMuiThuZWrdpSrrJDnlBbSZFbaxdNlSs1yfIaTJla9j9va3dlXdg8M8dna/0d02JWY7TFw3YHgHhkmKZpWl3EQPL5fHK5XPJ6vcrMzLS6HNtb/OQb2lVRq+B5bxuHYWjWhBw9tXS65TVICtvZn6+kyB2zGqMtHrY7AMRSpPtvS8+RKSsr06233qr8/HwZhqHnnnsubPqdd94pwzDCHrfccos1xUKVniaVHfKE7UwlKWiaKjvkUVWt3/IaugsxkmJWY7TFw3YHgHhlaZDx+/2aMmWK1q1b122fW265RR988EHosWXLlhhWiPMdrW/ucfqRuoHfofZWQ29iUWO0xcN2B4B4Zek5MvPnz9f8+fN77ON0OpWXlxejitCTMVmpPU4vzB74k097q6E3sagx2uJhuwNAvIr7y6937NihkSNH6oorrtCyZctUV1fXY/9AICCfzxf2QHSMc6erpMgth2GEtTsMQyVF7phcRdNbDSVF7m7njVWN0RYP2x0A4lVcB5lbbrlFTz31lLZv367//u//1s6dOzV//nwFg8Fu5yktLZXL5Qo9CgoKYljx4Ld20dTQSbWdZk3I0dpFU+OihrWLpqp4XPZF88wcnx3TGqMtHrY7AMSjuLlqyTAMbdu2TQsWLOi2T2VlpcaPH69XX31Vc+bM6bJPIBBQIBAIPff5fCooKOCqpSirqvXrSJ3f0vuZ9FRDVa1feyrrBt19ZOJhuwNALER61ZKt7iMzbtw45eTkqKKiotsg43Q65XQ6Y1zZ0DM2x/odaU81xEN9A2GwrhcA9Fdcf7V0oePHj6uurk6jRo2yuhQAABAHLD0i09TUpIqKitDzqqoqlZeXKysrS1lZWXrooYe0cOFC5eXl6fDhw/ra176mCRMmaN68eRZWDQAA4oWlQeatt97S7NmzQ89Xr14tSVqyZInWr1+v/fv361e/+pUaGhqUn5+vm2++Wd/97nf56ggAAEiKo5N9Bwo/UQAAgP3Y4icKAAAALgVBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2Fai1QXYVaWnSUfrm1WYnaaxOWlWl2N7kW7PaGz3zjEchqGgaaowO02maepofbNO+1pU42vRtaNH6GNF7j4t7x/jSkFTYeNeOH9X43Y1/4XTeL8BQDiCTB81NLdq1ZZylR3yhNpKitxau2iqXKlJFlZmT5Fuz2hs967G6EligqH2DrPX5UU6bkmRW99bMEnffO69sL4zx2fLNKXdlXUXzVM8LluGIb1+uC5sHN5vAHCOYZqm2Xs3+/L5fHK5XPJ6vcrMzLzk8RY/+YZ2VdQqeN5mcxiGZk3I0VNLp1/y+ENNpNszGtu9qzH6orvlRTquwzCUOSxRvrPt/a6hpzoAYDCJdP/NOTJ9UOlpUtkhz0U7oaBpquyQR1W1fosqs6dIt2c0tnt3Y/RFV8vry7hB09SZ5rZLqqG7OgBgqCLI9MHR+uYepx+pY8fSF5Fuz2hs997G6IvzlxfNcS+lDgAYqggyfTAmK7XH6YXZnITZF5Fuz2hs997G6IvzlxfNcS+lDgAYqggyfTDOna6SIrcchhHW7jAMlRS5uZqkjyLdntHY7t2N0RddLa8v4zoMQyNSky6phu7qAIChiiDTR2sXTdWsCTlhbbMm5GjtoqkWVWRvkW7PaGz3rsboSWJCeODobnmRjjtrQo7+d/mNF/WdOT5bxeOyu5yneFy2Zo4Pn8b7DQD+gauW+qmq1q8jdX7u6xElkW7PaGz3zjE6L6/u/IrmSJ1ftY0BnfSeDd1Hpi/L62ncC+fvatyu5r9wGu83AENFpPtvggwAAIg7XH4NAAAGPYIMAACwLYIMAACwLYIMAACwLUuDTFlZmW699Vbl5+fLMAw999xzYdNN09S3vvUtjRo1SsOGDdPcuXN16NAha4oFAABxx9Ig4/f7NWXKFK1bt67L6T/4wQ/02GOP6fHHH9fevXuVlpamefPmqaWlJcaVAgCAeJRo5cLnz5+v+fPndznNNE09+uij+uY3v6lPfepTkqSnnnpKubm5eu655/S5z30ulqUCAIA4FLfnyFRVVammpkZz584NtblcLs2YMUO7d+/udr5AICCfzxf2AAAAg1PcBpmamhpJUm5ublh7bm5uaFpXSktL5XK5Qo+CgoIBrRMAAFgnboNMf61Zs0Zerzf0qK6utrokAAAwQOI2yOTl5UmSTp06FdZ+6tSp0LSuOJ1OZWZmhj0AAMDgFLdBZuzYscrLy9P27dtDbT6fT3v37lVxcbGFlQEAgHhh6VVLTU1NqqioCD2vqqpSeXm5srKyNHr0aN1zzz363ve+p6KiIo0dO1YPPPCA8vPztWDBAuuKBgAAccPSIPPWW29p9uzZoeerV6+WJC1ZskSbNm3S1772Nfn9fn3pS19SQ0ODbrzxRr300ktKSUmxqmQAABBHDNM0TauLGEiR/gw4AACIH5Huv+P2HBkAAIDeEGQAAIBtEWQAAIBtEWQAAIBtEWQAAIBtEWQAAIBtEWQAAIBtEWQAAIBtEWQAAIBtWfoTBbCHSk+TjtY3qzA7TWNz0i563tf5L7Vfb8vaW1Wn2qZWudOdmjEuOzRWJOOf38c0Te2tqtObVWfU2NKmmyfm6fbrCiKuvbPNYRgKmqYchnSioUW1jS0yZWhkhlP5w1MUNBXq01Vt5499rM6v8uMNunb0CH2syC1J2nnw9EVt3a3PpW7faIjG6wwAnfiJAnSroblVq7aUq+yQJ9Q2IjVJZ5rbQs9Litxau2iqXKlJEc3fVf9I+/VW61eefluvH667aNr1hSOU5EgIm3bh+F3V0JXEBEP/u3yWrv6Iq9v5Zo7PlmlKuysvriUSnbWZMnusKTPl3P+H+FraQ20jUpP0v8tvVMawxB7n7ev2jYZovM4Aho5I998EGXRr8ZNvaFdFrYI9vEUchqFZE3L01NLpEc3fVf9I+/VWa28hpKe6I1nXTokJhir+65N9nq+vtUnq19gjUpP00Y8M73Hevm7faIjG6wxg6OC3lnBJKj1NKjvk6XUnGjRNlR3yqKrWH9H8F/aPtF8ktfbF+eNHuq6d2jtM/eat6j7P19fa+jv2mea2Xufty/aNhmi8zgDQFYIMunS0vrlP/Y/Uhe+Iepu/s3+k/XrS11ovHL8/8+86XHtJy40XkWzfaIjG6wwAXeFkX3RpTFZqn/oXZoeftNnb/J39I+3Xk77WeuH4/fl2ddb4nEtabryIZPtGQzReZwDoCkdk0KVx7nSVFLnlMIwe+zkMQyVF7ouuPulu/gv7R9ovklr74vzxI13XTokJhm6/rqDP8/W1tv6OPSI1qdd5+7J9oyEarzMAdIUgg26tXTQ1dNJppxEXXF0ya0KO1i6aGvH8XfWPtF9vtc4cn93ltOmFIy6aduH4XdXQlc6rlnqab+b4bBWP67qWSHTW1ltNmSmJoSuXOnVetdTbvH3dvtEQjdcZAC7EVUvoVVWtX0fq/KH7flz4vK/zX2q/3pa1t7JOnqbARfeRiWT88/tI0t7KOr15pF6+sz3fR6arsTvbEhMMtXeYSkwwdPzMWdU1BWRKcqc79ZERw0LT2ju6vo/M+WMfP9Ost4+dCbtnzP875Lmorbv1udTtGw3ReJ0BDH5cfv0hggwAAPbD5dcAAGDQI8gAAADbIsgAAADbIsgAAADbIsgAAADbIsgAAADbIsgAAADbIsgAAADbIsgAAADbIsgAAADbIsgAAADbIsgAAADbIsgAAADbiusg8+1vf1uGYYQ9rrzySqvLAgAAcSLR6gJ6M3HiRL366quh54mJcV8yAACIkbhPBYmJicrLy7O6DAAAEIfi+qslSTp06JDy8/M1btw43XHHHTp27FiP/QOBgHw+X9gDAAAMTnEdZGbMmKFNmzbppZde0vr161VVVaWPfexjamxs7Hae0tJSuVyu0KOgoCCGFQMAgFgyTNM0rS4iUg0NDRozZox+/OMfa+nSpV32CQQCCgQCoec+n08FBQXyer3KzMyMVakAAOAS+Hw+uVyuXvffcX+OzPmGDx+uyy+/XBUVFd32cTqdcjqdMawKAABYJa6/WrpQU1OTDh8+rFGjRlldCgAAiANxHWT+4z/+Qzt37tSRI0f0+uuv69Of/rQcDocWLVpkdWkAACAOxPVXS8ePH9eiRYtUV1cnt9utG2+8UXv27JHb7ba6NAAAEAfiOshs3brV6hIGhUpPk47WN6swO01jc9IsH2cgVXqatLeqTpKhG8Zla2xOWljdpml2+e/Ofnur6uRpbNXIDKfyhw/TiYazqm1sUU5GSmi8zuVcOG8kzx2G9NqB0yqvbpA7PUVTClzyNAWUne7UlMtcKq/26khtk5pagjJl6uaJeZo2ZkQ3Yxkqrz6jen+r3BlOVZxuUlNLuwqyUlWUm6EEQ/rT+6fUHAjq6lEZykhNUqJhqN00NSpzmNyZTp32tujvNT4lGIaKcjP0keEpKq9uUL2/Te70ZHmaWpUgqSjv3LSgKTkMQ68dOK3qM82aNzFPt1/X9ZWBF27PGeOyZZqm9lbVy5BCzzvXJWiachjSiYaWi7Z5b9u3P++TgX4v2+HvBRgMbHXVUn9EetbzYNTQ3KpVW8pVdsgTaispcmvtoqlypSbFfJyB1NDcqq88/bZeP1wX1u4alijv2fZe5x8+LEkNZ9t67Te9cIQcCQnaXfmP5YxITdKZ5raIn1+KaI4VLQ5D+v2KG3X1R1ySun8t+uvC1/DCbdCX92Is3st2+HsB7CDS/TdBZhBb/OQb2lVRq+B5L7HDMDRrQo6eWjo95uMMpMVPvhG240BsJSYYqvivT0qK/WvRl/diLN7Ldvh7Aewg0v13XJ/si/6r9DSp7JAn7MNUkoKmqbJDHlXV+mM6zkDqrBHWae8w9Zu3qi15LSJ9L8bivWyHvxdgsCHIDFJH65t7nH6kLrIP1GiNM5B6qxGxsetwraWvRW/vxVi8l+3w9wIMNgSZQWpMVmqP0wuzIzv5MFrjDKTeakRszBqfY+lr0dt7MRbvZTv8vQCDDUFmkBrnTldJkVsOwwhrdxiGSorcEV9FEa1xBlJnjbBOYoKh268rsOS1iPS9GIv3sh3+XoDBhiAziK1dNFWzJuSEtc2akKO1i6ZaMs5AWrtoqmaOz76o3TUssjsMDB8W2dUk0wtHqHhc+HJGXHAlSm/PL0U0x4oWhyH97/JZoefdvRb9deFreOE26Mt7MRbvZTv8vQCDCVctDQFVtX4dqfNf8v0sojXOQKqq9WtvZZ1MKXQPkvPrltTlvzv77a2sk6cpIHe6Ux8ZMUzHz5xV3Yf3eTn/PjIXbotInycmGNpx4LTePnZG7vQUXTN6uDyNAWWlJ+uaguEqr25QlcevpkC7TPPcfWSuK8zqdqy/VjeotimgkRkp591HZpgm5GYoMcHQn/52Sk2Bdl2dn6mMYYlKSkhQW0eH8l3DlJPhVG1jQH876ZXDMDQhN0OXjRj24Zjn7v1yujEgh6HQtPYOM7QOR+t7vo/MhdtzxocBcE9lXeg+Mp2vQWKCERq7q23e2/btz/tkoN/Ldvh7AeIZl19/iCADAID9cPk1AAAY9AgyAADAtggyAADAtggyAADAtggyAADAtggyAADAtggyAADAtggyAADAtggyAADAtggyAADAtggyAADAtggyAADAtggyAADAtggyAADAtggyAADAtggyAADAtggyAADAtggyAADAtggyAADAtggyAADAtggyAADAtggyAADAtggyAADAtggyAADAthKtLiAS69at0w9/+EPV1NRoypQpWrt2raZPn251WTFX6WnS0fpmFWanaWxOWr/mM02zy39HOl6lp0l7q+okGbphXHaP83XVd+fB0yo/3qBRmcPkznSGLftc/3r9X41P3rNtGutO15TLhutEw1nVNrZIhpSTnhI21msHPcpJT9Y/T87vcX0qPU36w/6Tqve3yWFIR+ubNSYrVf90xUidaDgrQ1L+8BSdaGjRm1X1qvG2KH94iq4rzNKMcdkyTVN7q+plSJoxLlvH6vx6rvyEDEkj0pJ1tK5Zza3tSk1OVGFWquqaW9XU0q4pBcN1pNav90565U5zqt001dIWlCS1tAXVFGhXW7BD2alOdchUU6Bdza1BtQU7lOxIkDs9RVlpSeqQlJmSqNONAbW0BTUmK00pTofOtgR1qqlFjWfb1XC2TanJDl01KlPpzkSd8rXoA+9ZNbe2K8mRoIn5Lg1Ldsjf0q7JBcMlSb/bd0I1vmYlJiQo1ZmoQFtQ6SmJOtsaVL2/TR2SEhOkjJREyZQazrbL/HCbOh3n/tsalBIMKc2ZqOREQ8Gg+eE6mDKMc/O3B6Xgh/MZUmgMSUpOMNTeYcrUuXE6zHPTHYZkfvjvxAQpwTBkylQwKHVcMIbDkJIcUmv7uWlJhpSd4VR7sEONLe1qC5rq0Ln/c0t3OtTe0aHmtn+M4JBkJJwbx5FgqC1oyjTP1SKda5NMtXV8WLNDys1M0ZWjMuUPtKvG2xJaN39ruxpb2tXyYee0ZIfagh0yJaU7E5XoOLe+ZwNBtbZ3fLhcQ5KhlKQEOZMS1N5hKjHBUHvHuXUfkZqsD3wtOhtoV4cpJSYYcqUmqT1oynu2TYakXNcwpTsdagy0KxjsUGuHqaSEBKUmO1Q0MkONgTZVevxKS07UOHda6IUYk5WmdrNDJxvOyvPh+6upJaikxASlJhpqCLQrNyNFyYkJSk9J1DUFw1VxukmnfC3Kyxymo7VNOnKmWaPSUzRhVIYuH5mhq/Iz5TCk8uoG7T/uVWZKksa701Th8auppU0FWakyTen4mWadONMsU9LEfJcKslM1KjNFQfNcebWNAe06XCvDkGaOz5E73am/1/iUk+7U5MtcCprnXrMTDS2qbQyovjkg0zQ056qR+liRO/R5c+3oEfpYkTv0Wegw9OG8hoKmGfpv5+fi+X/rffmsjYaePucvnNbffUKs6o0VwzRNs/du1nnmmWe0ePFiPf7445oxY4YeffRR/eY3v9HBgwc1cuTIXuf3+XxyuVzyer3KzMyMQcXR19DcqlVbylV2yBNqKylya+2iqXKlJvVpvu70Nl5Dc6u+8vTbev1wXVh78bhsPf7v08Lm665vYoLU3nHx2DPHZ6utvUNvHj3Ta52dHAlSsIuxLlwfU6a++NRbevNI5GMDGHwSPwzMfTVzfLbW3zGtx8/aaOjpc96UedG0EalJOtPcdlHfga4zknqjVUOk+++4DzIzZszQ9ddfr5/+9KeSpI6ODhUUFGjlypW6//77e51/MASZxU++oV0VtQqe91I5DEOzJuToqaXdH5nqar7u9Dbe4iff6DYQlRS5w+brqW+sdK6PJMtrAWBvF37GDYSePucl9fpZHsk+IZr6u1/qi0j333F9jkxra6v27dunuXPnhtoSEhI0d+5c7d69u8t5AoGAfD5f2MPOKj1NKjvkuegNHDRNlR3yqKrW36f5utPTeJ1jdef8+XrrGyud6xMPtQCwt54+a6Oht8/5SD7Le9snRFN/90sDJa6DTG1trYLBoHJzc8Pac3NzVVNT0+U8paWlcrlcoUdBQUEsSh0wR+ube5x+pK7rN0xv8/VlvEjG6pyvv8sFgHjW3WdtNETzc3Mg6+zU3/3SQInrINMfa9askdfrDT2qq6utLumSjMlK7XF6YXbXJ1f1Nl9fxotkrM75+rtcAIhn3X3WRkM0PzcHss5O/d0vDZS4DjI5OTlyOBw6depUWPupU6eUl5fX5TxOp1OZmZlhDzsb505XSZH7w6sa/sFhGCopcnd7lnh383Wnp/E6x+rO+fP11jdWOtcnHmoBYG89fdZGQ2+f85F8lve2T4im/u6XBkpcB5nk5GRNmzZN27dvD7V1dHRo+/btKi4utrCy2Fq7aGrohK9OsybkaO2iqX2erzu9jbd20VTNHJ99UXvxuOyL5uuub2I377aZ47N1/ZgREdXZydHLO7dzfdYumqrphX0bG8Dgk5gQ2f/UXWjm+Is/4wZCT5/zXU0bccGVQZHsE6Kpv/ulgRD3Vy0988wzWrJkiZ544glNnz5djz76qH7961/rwIEDF50705XBcNVSp6pav47U+ft8vf7580nq8t+RjldV69feyjqZUq/3kemq7/875NHbx84o3zVMORnh95GpqvVrT2WdKk416kxzm8a603RNwXAdP3NWdU0BSVJ2ujNsrO1/Px26j0xP61NV69cL+0+qtqlViQmGjtT5NSYrVTddOVLHz5y7j8xHRgzT8TNnte9IvU42hN9HRpL2VNaF7i1x/Eyztr19QpKprDSnjtT51dIaVEqyQ4XZaar3t6qppU2TC4braF2z3j3RIHe6U0HTVEvrubugtLQF1RgIqi0YVHaaUx3mufvInG0NqjVoKtlhyJ2Rouy0ZHXIVIYzSZ7GFp1tC2pMdppSkh1qaQ2qxteiprPtamhp07Akh64elak0Z6JOh+4jE1Siw9DEfJdSkxNDdSUYhn6777g+8J67j0yaM1EtbUFlpCSpubWd+8hEeB+Zs4GgPvCelaFz9Z27j0wwdL+gi+8jc+4+MWdb29Xa1iEjwZDDOLdl/nEfmX/criAxQcpKdeoDX4uaA2093kcmw5moxkCb2oOm2jo6QveRmTAyQ/5Auyo8TUpLTtT4kWmSDJmmqTHZaQp2mDrZcFanG1sU+PB9meRIUGpSghpa2pSXkaLkpHPvkakFI0L3kcnNTNGxOr+q6v0alZGiCXn/uI9MYoKhv1Y3aH91g9JTkjRhZLoqTjepqaVdBVnDZEo6Xn9Wx880y5SpSR8ZrsuyhinfNUztHaYMSfX+Vv3lkEeGYah4fLZyM1P0t5Ne5aQ7NaVgeOh+O52fEWf8rQqaCt1HpvPzpvM+Mp2fhZ2XYl/4387PxfP/1mN9dKGnz/kLp/V3nxCrei/VoLn8WpJ++tOfhm6Id8011+ixxx7TjBkzIpp3MAUZAACGikEVZC4FQQYAAPsZFPeRAQAA6AlBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2Fai1QUMtM4bF/t8PosrAQAAkercb/f2AwSDPsg0NjZKkgoKCiyuBAAA9FVjY6NcLle30wf9by11dHTo5MmTysjIkGH072fc44HP51NBQYGqq6uH5G9GDeX1H8rrLrH+rD/rP1TX3zRNNTY2Kj8/XwkJ3Z8JM+iPyCQkJOiyyy6zuoyoyczMHHJv5vMN5fUfyususf6sP+s/FNe/pyMxnTjZFwAA2BZBBgAA2BZBxiacTqcefPBBOZ1Oq0uxxFBe/6G87hLrz/qz/kN5/SMx6E/2BQAAgxdHZAAAgG0RZAAAgG0RZAAAgG0RZAAAgG0RZGzoX/7lXzR69GilpKRo1KhR+vznP6+TJ09aXVZMHDlyREuXLtXYsWM1bNgwjR8/Xg8++KBaW1utLi1mvv/972vmzJlKTU3V8OHDrS5nwK1bt06FhYVKSUnRjBkz9MYbb1hdUkyUlZXp1ltvVX5+vgzD0HPPPWd1STFVWlqq66+/XhkZGRo5cqQWLFiggwcPWl1WTKxfv16TJ08O3QSvuLhYL774otVlxS2CjA3Nnj1bv/71r3Xw4EH99re/1eHDh3XbbbdZXVZMHDhwQB0dHXriiSf0t7/9TY888ogef/xxff3rX7e6tJhpbW3V7bffrmXLllldyoB75plntHr1aj344IN6++23NWXKFM2bN0+nT5+2urQB5/f7NWXKFK1bt87qUiyxc+dOLV++XHv27NGf/vQntbW16eabb5bf77e6tAF32WWX6eGHH9a+ffv01ltv6eMf/7g+9alP6W9/+5vVpcUnE7b3/PPPm4ZhmK2trVaXYokf/OAH5tixY60uI+Y2btxoulwuq8sYUNOnTzeXL18eeh4MBs38/HyztLTUwqpiT5K5bds2q8uw1OnTp01J5s6dO60uxRIjRowwN2zYYHUZcYkjMjZXX1+vp59+WjNnzlRSUpLV5VjC6/UqKyvL6jIQZa2trdq3b5/mzp0baktISNDcuXO1e/duCyuDFbxeryQNub/1YDCorVu3yu/3q7i42Opy4hJBxqbuu+8+paWlKTs7W8eOHdPzzz9vdUmWqKio0Nq1a/XlL3/Z6lIQZbW1tQoGg8rNzQ1rz83NVU1NjUVVwQodHR265557NGvWLE2aNMnqcmLi3XffVXp6upxOp+6++25t27ZNV199tdVlxSWCTJy4//77ZRhGj48DBw6E+v/nf/6n3nnnHb3yyityOBxavHixTBvfpLmv6y9JJ06c0C233KLbb79dX/ziFy2qPDr6s/7AULF8+XK999572rp1q9WlxMwVV1yh8vJy7d27V8uWLdOSJUv0/vvvW11WXOInCuKEx+NRXV1dj33GjRun5OTki9qPHz+ugoICvf7667Y99NjX9T958qRuuukm3XDDDdq0aZMSEuydyfvz+m/atEn33HOPGhoaBrg6a7S2tio1NVXPPvusFixYEGpfsmSJGhoahtRRSMMwtG3btrDtMFSsWLFCzz//vMrKyjR27Firy7HM3LlzNX78eD3xxBNWlxJ3Eq0uAOe43W653e5+zdvR0SFJCgQC0Swppvqy/idOnNDs2bM1bdo0bdy40fYhRrq013+wSk5O1rRp07R9+/bQDryjo0Pbt2/XihUrrC0OA840Ta1cuVLbtm3Tjh07hnSIkc699+38GT+QCDI2s3fvXr355pu68cYbNWLECB0+fFgPPPCAxo8fb9ujMX1x4sQJ3XTTTRozZox+9KMfyePxhKbl5eVZWFnsHDt2TPX19Tp27JiCwaDKy8slSRMmTFB6erq1xUXZ6tWrtWTJEl133XWaPn26Hn30Ufn9ft11111WlzbgmpqaVFFREXpeVVWl8vJyZWVlafTo0RZWFhvLly/X5s2b9fzzzysjIyN0XpTL5dKwYcMsrm5grVmzRvPnz9fo0aPV2NiozZs3a8eOHXr55ZetLi0+WXvRFPpq//795uzZs82srCzT6XSahYWF5t13320eP37c6tJiYuPGjaakLh9DxZIlS7pc/9dee83q0gbE2rVrzdGjR5vJycnm9OnTzT179lhdUky89tprXb7OS5Yssbq0mOju73zjxo1WlzbgvvCFL5hjxowxk5OTTbfbbc6ZM8d85ZVXrC4rbnGODAAAsC37n1wAAACGLIIMAACwLYIMAACwLYIMAACwLYIMAACwLYIMAACwLYIMAACwLYIMAACwLYIMAACwLYIMgLgVDAY1c+ZMfeYznwlr93q9Kigo0De+8Q1J0qpVqzRt2jQ5nU5dc801FlQKwCoEGQBxy+FwaNOmTXrppZf09NNPh9pXrlyprKwsPfjgg6G2L3zhC/rsZz9rRZkALMSvXwOIa5dffrkefvhhrVy5Uh//+Mf1xhtvaOvWrXrzzTeVnJwsSXrsscckSR6PR/v377eyXAAxRpABEPdWrlypbdu26fOf/7zeffddfetb39KUKVOsLgtAHCDIAIh7hmFo/fr1uuqqq/TRj35U999/v9UlAYgTnCMDwBZ++ctfKjU1VVVVVTp+/LjV5QCIEwQZAHHv9ddf1yOPPKI//OEPmj59upYuXSrTNK0uC0AcIMgAiGvNzc268847tWzZMs2ePVtPPvmk3njjDT3++ONWlwYgDhBkAMS1NWvWyDRNPfzww5KkwsJC/ehHP9LXvvY1HTlyRJJUUVGh8vJy1dTU6OzZsyovL1d5eblaW1strBxALBgmx2cBxKmdO3dqzpw52rFjh2688cawafPmzVN7e7teffVVzZ49Wzt37rxo/qqqKhUWFsaoWgBWIMgAAADb4qslAABgWwQZAABgWwQZAABgWwQZAABgWwQZAABgWwQZAABgWwQZAABgWwQZAABgWwQZAABgWwQZAABgWwQZAABgW/8ftGPZWf4kxpUAAAAASUVORK5CYII=",
      "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": "31c74ca5-2c14-48fa-853f-f7a1a673197e",
   "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
}