gp-bbn-toy.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a11541a8-c5da-4f2b-8a44-aad993c2bc4d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import functools\n",
    "import sys\n",
    "\n",
    "from tqdm.notebook import tqdm\n",
    "import numpy as np\n",
    "import torch\n",
    "import gpytorch\n",
    "from gpytorch.models import ApproximateGP, PyroGP\n",
    "import pyro\n",
    "from pyro.nn import PyroParam, PyroModule, PyroSample\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "sys.path.append(\"../../fermi-gce-gp/notebooks/\")\n",
    "\n",
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "df09d0b4-7e29-47b0-afeb-506c61d65252",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pylab as pylab\n",
    "import warnings\n",
    "import matplotlib.cbook\n",
    "\n",
    "from plot_params import params\n",
    "\n",
    "warnings.filterwarnings(\"ignore\",category=matplotlib.cbook.mplDeprecation)\n",
    "\n",
    "pylab.rcParams.update(params)\n",
    "cols_default = plt.rcParams['axes.prop_cycle'].by_key()['color']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "443756ba-7f77-48f2-9f91-44f82c657e2b",
   "metadata": {},
   "source": [
    "Define an \"observable\" `obs(T_nu, T_gamma)` that represents a mapping from latent temperatures/Hubble to observations (e.g.,  elemental abundances)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ffdd403b-969d-4f48-9fde-3470e3d31f59",
   "metadata": {},
   "outputs": [],
   "source": [
    "def obs(T_nu, T_gamma):\n",
    "    return torch.log(T_nu ** 4 * T_gamma + 0.1)\n",
    "#     return np.log(T_nu.detach() ** 4 * T_gamma + 0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7329a64-29f3-4491-934e-0739ad53cf0b",
   "metadata": {},
   "source": [
    "GP regression model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "4f069eaf-845f-4a56-bc2d-bfc092c6e0b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "class GPRegressionBBN(ApproximateGP):\n",
    "    def __init__(self, obs_sample, name_prefix=\"bbn\", num_inducing=30):\n",
    "        \"\"\" Module for fitting a latent GP to observations\n",
    "            :param obs_sample: Sample of observed values\n",
    "            :param name_prefix: Run tag\n",
    "            :param num_inducing: Number of inducing points for sparse GP\n",
    "        \"\"\"\n",
    "\n",
    "        self.name_prefix = name_prefix\n",
    "        self.num_inducing = num_inducing\n",
    "        \n",
    "        inducing_points = torch.linspace(0.1, 20, num_inducing)\n",
    "        \n",
    "        # Initialize GP variational strategy\n",
    "        variational_strategy = gpytorch.variational.VariationalStrategy(self, inducing_points, gpytorch.variational.CholeskyVariationalDistribution(self.num_inducing), learn_inducing_locations=False)\n",
    "        super().__init__(variational_strategy)\n",
    "        \n",
    "        # Specify GP mean and covariance structure\n",
    "        self.mean_module = gpytorch.means.ZeroMean()\n",
    "        self.covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel())\n",
    "        \n",
    "        self.obs_sample = obs_sample\n",
    "        \n",
    "    def forward(self, x):\n",
    "        \"\"\" Forward samples from the GP\n",
    "        \"\"\"\n",
    "        \n",
    "        mean = self.mean_module(x)\n",
    "        covar = self.covar_module(x)\n",
    "\n",
    "        return gpytorch.distributions.MultivariateNormal(mean, covar)\n",
    "\n",
    "    def guide(self, x):\n",
    "        \"\"\" Forward samples from the GP variational guide distribution\n",
    "        \"\"\"\n",
    "\n",
    "        # Get q(f) - variational (guide) distribution of latent function\n",
    "        function_dist = self.pyro_guide(x)\n",
    "\n",
    "        # Use a plate here to mark conditional independencies\n",
    "        with pyro.plate(self.name_prefix + \".data_plate\", dim=-1):\n",
    "            # Sample from latent function distribution\n",
    "            gp_guide_sample = pyro.sample(self.name_prefix + \".f(x)\", function_dist)\n",
    "\n",
    "    def model(self, x):\n",
    "        \"\"\" GP model\n",
    "        \"\"\"\n",
    "        \n",
    "        pyro.module(self.name_prefix + \".gp\", self)\n",
    "        \n",
    "        with pyro.plate(self.name_prefix + \".data_plate\", dim=-1):\n",
    "        \n",
    "            # Get prior distribution of latent function and sample from it\n",
    "            function_dist = self.pyro_model(x)\n",
    "            function_samples = pyro.sample(self.name_prefix + \".f(x)\", function_dist)\n",
    "                        \n",
    "            # (Exponentiate) to ensure non-zero GP\n",
    "            T_nu = function_samples#.exp()\n",
    "            \n",
    "            # Get corresponding observations\n",
    "            obs_sample = torch.Tensor(obs(T_nu, T_gamma))\n",
    "            \n",
    "#             # If using custom likelihood\n",
    "#             log_likelihood = self.log_likelihood(obs_sample)\n",
    "#             return pyro.factor(self.name_prefix + \".log_likelihood\", log_likelihood)\n",
    "                            \n",
    "            # Sample from observed distribution\n",
    "            pyro.sample(self.name_prefix + \".y\",\n",
    "                pyro.distributions.Normal(loc=obs_sample, scale=1.), \n",
    "                obs=self.obs_sample)\n",
    "\n",
    "#     def log_likelihood(self, obs):\n",
    "#         \"\"\" Custom log-likelihood\n",
    "#         \"\"\"\n",
    "#         delta_obs = 1.\n",
    "#         chi = (obs - self.obs_sample) / delta_obs\n",
    "#         return - chi.square()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1aa9eaee-b021-480e-b997-8697f6ba8e83",
   "metadata": {},
   "source": [
    "Define a toy latent function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "e84e0e77-6906-43d0-bffe-6d0ef031d77e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def T_nu(T_gamma):\n",
    "    \"\"\" The true latent function\n",
    "    \"\"\"\n",
    "    return (T_gamma ** 3 + 2.).log()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "596b478a-32bb-4ad5-8f62-e92d7acf9e8a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '$T_\\\\nu(T_\\\\gamma)$\\\\,[arb.]')"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAFRCAYAAAC1yZDjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAt5ElEQVR4nO3deXwV9b3/8dc3C4GQ5SQhhIgiBCxYbaEQVFBRWRXFKirebmpvK1p7218XK1rbn9Heq4LW9nevrTXaattrLYu7VSugrYALm3VrUSHsBJKQnCQQQnKS7++PcxISsp3MOcmc5f18PHiUzJyZ83E68GbmO9/PGGstIiIiTiS4XYCIiEQvhYiIiDimEBEREccUIiIi4phCREREHEtyu4BwMcbUAclAudu1iIjEkFyg0Vqb2tlKEyuP+BpjGhMSEpLy8vJ6vW1tbS3p6emOv7uuro7U1E6Pb59u6/b2bh63ULd3+7iHcuzcrt3N465zrv+P3YEDB2hubvZZa5M7/YC1NiZ+Afvy8/OtE/7D4NykSZNc2dbt7d08bqFu7/ZxD+XYuV27m8dd55xzTo9dfn6+BfbZLv7u1ZiIiIg4phAJg4ULF7qybSRs7+Z3R/Nxd/O7o/m4hyqa/9vdPnZdiaUxkX35+fn5+/btc7ItsXIc+pOOm3M6ds7ouDnn9NidcMIJlJaWllprT+hsva5ERETEMYWIiIg4phAB7rjjDrdLiEo6bs7p2Dmj4+ZcXx07jYmIiEiXNCYiIiJ9RiEiIhLjmurqaKio6JN9x0zvLBER8bPWUr97N9UbN1GzaTO1//yIrKlTKfjRzWH/rpgKkbq6OgoLCwH/xJxInZwjIhJuTfVHqX3/Pao3bqJ64yYaysvara/Z/C62qQmTmBjU/oqLiykuLqbCfwXTZdMuDayLiESpo6X7qd64keqNG6l9/wOafY0dPpMybBiZkyaRMXEimRO/gEnq3bVDTwPrMXUlIiISy6zPx6EtW6jesBHv+g3U79nd4TMmKYn0008nc9IkMgsnkTJ8OMaYPqtJISIiEsF8hw5Rs2kz3vXrqd64iaa6wx0+MyA7h8zJhWQWFpI+/vMkDhrUb/UpREREIkz9vn1Uv7Me7/r1HPrnv7DNTcd9wpA2diyZZ0wmc3Ihg0aO7NOrje4oREREXGabmjj8ySd431mP9531nd6mShyUSsbEiXgmF5JROInkzEwXKu1IISIi4oLmo0epee89vG+9TfWGjTRWezt8JiUvj8zJk/GceQZpp59OQi8HxftD5FUkIhKjfLW1/kHxt9+mevNmmo8ePe4TgdtUZ07Gc8YZDBwxwrXbVMFSiIiI9KGGigq8b7+D9623qf3www7jGwlJyaRPmIDnrDPxTJ5McnaWS5U6oxAREQmz+tJSvG++RdWbb3H4k487rE8anEbmGZPxnHUmGV/4Qr8+TRVuChERkRBZazmyc2drcBzZuaPDZwbkDPFfbUw5i7TTTovI8Q0nYuO/QkSkn1lrqdtWgnfdOqrefIv6fXs7fGbg8OFkTZmCZ8pZpJ5ySsSPbzihEBERCZK1lrqtW6la9yZV69ZxdP/+Dp9JHTWKrKlT8UydwqARI1yosn8pREREutEaHGvX+YPjwIEOnxn8mbFkTZ1C1tSppOQPc6FK9yhERESOY63lyPbtVL6xpssrjrRx48g65xw8U6aQMjTXhSojQ0yFiFrBi0gojuzeTdUba6h8Y02nYxxp48aRde45ZE2dyoAhQ1yosP+oFbyISBCOlu6ncu0aKv++ptOnqtLGBoLj7NgPjs6oFbyIyHEaK6uoXLuWqjfWcOjjLR3Wp44eQ/a555B17jmkDB3qQoXRQyEiInGh6fBhqt56m8q//53a997H2uZ26weNGEH2tGlkTTuXgfn5LlUZfRQiIhKzmn0+ajZvpvL1v+F9Zz3NjQ3t1qcMzSP7PH9wpI4c6U6RUU4hIiIxxVrL4Y8/pvL1v1G5Zg2+2tp265MyMsk+9xyyzz+PwWPHxuQEwP6kEBGRmHB0/34Ovv43Kl//G/Wl7R+wSUhJwXPWWeRccD4Z48f3+j3j0jUdSRGJWk2HD1O17k0OvvYatR991G6dMQmkT5hAzvTz8Zx5ZlQ3OYxkChERiSq2qYna9z+gYtVqvG+/TXND+3dypI4aRc706WRPmxZ1bdWjUUSHiDGm7WxBj7V2iWvFiIir6ktLObhqNQdXv0bDwYp265Kzssm54HyyLzhfA+T9LGJDxBhzS9vQMMZ4jDGLrbWL3KxLRPpPU309VWvXcXDVqg63qxKSkvFMOYucGdM1zuGiSD7qo9v+YK31GmM8LtUiIv2k5emqildXUbVmDU31R9qtH/yZsQyZMZ2s86aRNHiwS1VKi0gOkUJjzERr7Wa3CxGRvuerqeHg63+j4tVXObJrV7t1yZkecmZMJ2fG9Lhorx5NIjlEFgGbjDGLrLVLAuMjupUlEkOstf5B8r/+laq33sL6fK3rTEIimZMmkjNrFpmTC2PmTYCxJmL/X7HWrjLGTMIfJLcBk6y1XpfLEpEwaKyu5uCq1VT89dUOczpShuUzZPYscmZMZ0B2tksVSrAiNkSMMQXATCALWAxsM8bMstau6mqb2traLmef3nHHHRQVFfVFqSISBGsthz78iPKXX+5w1ZGQlIxn6lRy58wm7XOnaxa5C4qKirjzzju7Wp3e1YqIbQVvjHnYWntDm5+vBB4BRnV2RaJW8CKRyXf4MAdXr6b85Veo37On3bqBJ55I7oUXknPB+SRlZLhToHQrKlvBBwJjZdtl1toVba5OVrhSmIgE7fDWbZS/9BKVf3+j3YRAk5RE1tlnk3vhHNJOO01XHVEuIkOkG3pSSySCNTc0ULV2HWV/eYnDn3zcbl3KsHxyL7qQnBkzSM7UVUesiNQQWYX/1tXxVxyzNNlQJPI0lJdT/vIrlP/1VXw11a3LjUkg84zJ5M6dS8aE8ZiEBBerlL4QkSESmFi4yBhzC+ANLPYA97hWlIi0Y63l0D//SdkLL+J9621sc1PruuRMD0PmzCb3wjkMyM11sUrpaxEZIgDW2hJAvbJEIkxzQwOVb7xB2XMvULdje7t1aWPHkXvJxWSdPZWE5GSXKpT+FLEhIiKRpaGykvKXXqb85Vfa37JKSiJ72jSGzruEwWPGuFihuEEhIiLdqttWwoHnnqNyzZp2czuSs7IZevFchsyZTbLH416B4iqFiIh0YJubqd64kQPPPkftBx+0Wzf4M2PJu3QenrOnqhWJKERE5JjmhgYOvv46B559rt3EQGMS8EydSt5ll5I2bpyLFUqkUYiICL6aGspfepmyF/9CY7W3dXnioFSGzJ7F0HmXkJKX516BErEUIiJx7GhZGQeeeZaKlStpPnpsVvmAIUMYOm8eQ+bM1js7pFsKEZE4dGTHTvY/9RSVb6xpN78jdeQo8q64nKxzztF4hwRFZ4lIHDm0ZQv7ly3Hu2FDu+UZ48cz7IorSJ8wXr2spFdiKkTq6uooLCwEYOHChSxcuNDlikTcZ62l9h/vUbpsGbUffti63JgEss6eSt4V8zW/QzooLi6muLiYiooKgNSuPhexreB7S63gRdqz1lK9fj2lS5dx+NNPW5ebpCRypk9n2JVXMDA/38UKJRpEZSt4EXHONjdTte5NSpcu48jOHa3LE1IGknvRheRdfpneGChhoxARiRG2qYmqtev84bF7V+vyxNTBDJ13CXmXztOLnyTsFCIiUc42NVG5Zi2lS5e2myCYlJ5O3mWXkXvxXD2mK31GISISpWxzM1Vr17HvySfbhUdypoe8yy8jd+5FJA4a5GKFEg8UIiJRxjY3433rbfb96U8c2XXstlVypodhV8xnyEUXkTgwxcUKJZ4oRESihLWW6g0b2ffEE9SVlLQuT8rIZNgV88mdO1fhIf1OISISBWrf/4C9f/gjhz7e0rosKS3dHx6XXEziwIEuVifxTCEiEsEOf/ope//wR2r+8Y/WZYmpqeRddhl5l84jUQPm4jKFiEgEqt+7l71/fIKqdWtblyUMSGHovEsYdsV8ktLTXaxO5BiFiEgEaaispPTJP1Px6srWxogmKYncOXPIX7CA5OwslysUaU8hIhIBmurqOPDMs+x/5pk2LdkNOeefxwlf+TIpw4a5Wp9IVxQiIi6yPh8VK1ex74k/tXsZVOakSQy/7lpSR450rTaRYMRUiKiLr0SLlsd19zz2OPV7drcuHzzmFIZfdy0Z4z/vYnUi6uIrErHqduxgz6O/pea991qXDcgdyonXXkPWuedgEhJcrE6kPXXxFYkQjV4v+/73Cf+guW0G/M0R8xdcxdB5l5AwYIDLFYr0nkJEpI81+3yUv/Ai+/78Z5rq6gAwCYnkXnQh+V/6EsmZ6qwr0UshItKHqjduYvcjj1K/b2/rssxJkzjx37/OoBEjXKxMJDwUIiJ9oL60lD2PPNruXeYDhw/npG98g8zJhS5WJhJeChGRMGqqP8r+FSvY/9RTWJ8P8LcpOeHf/o3ceZeQkKQ/chJbdEaLhIG1lup31rOr+BEaystalw+ZMYPh115DcpZmmktsUoiIhOjo/v3sevgRqjceu3WVOnoMI751A2ljx7pYmUjfU4iIONTs83HgmWco/fMymhv8rUqSBqcx/JqvMWTObExiossVivQ9hYiIA7Uf/ZOdD/6q3WzzITNmMPzr15GcmeliZSL9SyEi0gu+Q4fY89jjVLz6auuyQSeNYMRNN5J++ukuVibiDoWISBCstVStXcfu4kdo9FYBkJA8gPwvXU3e5ZfrqSuJWzrzRXrQcPAgux76Dd533mldljFhAiNu+hYD8/NdrEzEfQoRkS5Ya6l4dSV7fve71nYlSenpnPTNb5J9wfkYY9wtUCQCxFSIqBW8hMvRAwfY+T8Ptuu0mz3tPE5a+E0NnEtcUCt4EQdsczPlL7/Cnscep/loPQADsnMY8e2b8Jwx2eXqRPqfWsGLBOloWRk7/9//UPP+sauPIbNnc+K/f52kwYNdrEwkcilEJO5Zazm4ajW7H3mUpiP+sY8BQ4Yw8rvfJeMLE9wtTiTCKUQkrjVWVbHzwV/hXb++ddmQWbM58Ru6+hAJRrchYoxZGuL+DWCBDdba+0Pcl0hYVb31NjsffBBfTQ0AyVnZjPzud8gsnORyZSLRo6crEWOtXRDql4QhjETCpunIEXY/8lsqVh6bdZ59zjmMuOlbJKWnu1iZSPTpKURKwvQ928O0H5GQHP7kE0ru+zlH95cC/necn3zTt8g+b5rLlYlEp25DxFp7azi+JFz7EXHKNjdz4Oln2Pu/T2Cb/C+LSv/c5xj1/e8xIDfX5epEopcG1iXmNVRWsuOBX7ROHDRJSQz/6lfJu/wyTEKCy9WJRLewhogx5nrgXmttTpj25wEWAl6gEvBaa1eFY98SH6o3bWL7A7/EV1MNQMqwfAoW/YjBY8a4XJlIbAj3lcgywjSOEgiQ5dbaWYGfC4BNgN4zKj1q9vnY979PsP+pp1qX5UyfzogbFpKY2mUHBxHppbCGiLW2Glgdpt3dBjzcZt8lxhg9eyk9aigvp2TxfRz6eAsACSkDOfnb3yLnggtcrkwk9jgOEWPMSGAJMCOwaBWwyFq7I/SyALiF4646rLXhelpMYlT1pk1s//kD+GprAUgdOYqCW29h4PDhLlcmEpschYgxZhSwHFjKsauFicByY8xVoQZJ4NYVQLYxZib+MZGJ1toloexXYpdtamLfk3+mdOky/PNbIffCCznp+m+SMGCAu8WJxDCnj6ZcYa0ttNbeZ61dHfh1n7V2MnBlGOoqwD+24rHWrggMpm82xjzc3Ua1tbUYYzr9VVRUFIayJBL5amr49M67KF26FLAkpKQw6oc/5ORv36QAEQlSUVFRp393lpaWAnQ5C9dRK3hjzAxrbadjH92t68X+ZwIrrbXmuOVVwChrrbeTbdQKPg4d3rqVbXffS0N5GQADTzyJ0bctYtCIES5XJhIb+qoVfGU36w463GdbJfhvYXWmANgchu+QKFexajW7fvVrmn2NAGSdfQ4j/893SBw0yOXKROKH0xCxxpiRx499BAbbs0MtCn9IecKwH4lBzT4fe377O8pefBEAk5DI8Guv8U8e1CtrRfpVjyFijNlAx7/QDTDKGFNy3LIs4PpQi7LWeo0xm40xnuNuXXmstboKiWON1TWULF5C7QfvA5CUkUHBLbeQMf7zLlcmEp+CuRIxwAK6vr10vO5udfXGPfhnqy8BMMZc2fJ7iU9Hduxk68/+k6NlBwBILShg9O0/JmXoUJcrE4lfwYTI9dbad/u8kuNYa1cYYxYaYxYSuBKy1i7q7zokMnjXr2f7fT+nqf4IANnTpnHyd75D4sAUlysTiW89hkhnAWKMyQisq+mLotp8d3Ff7l8in7WWA888y57HHqdl/sfwa77GsCuv1PiHSARwOrC+AngV0NsKpc80+3zseug3VLzqf3lUQspARv3wB2RNOcvlykSkhdMQedha+1RnK4wxE6y1/3Bekgj4Dh+m5J57W9u3D8gZwpg7fkrqqFEuVyYibTmdsb7NGDO9i3VXOy1GBOBoWRkf/2hRa4AMHnMK4x64XwEiEoGcXokswd/X6gu0b/1ugFH4O/CK9FrdthI+vfMuGqv8D/llTZnCyB/8QAPoIhEqlFbw19PxsV8D3BvCPiWOVW/eTMk9i1ufwMr74hc58d+/rrcPikQwpyGyqKvHfo0x94RQj8Spg6+9zo7//p/A+88NJ13/DfIuvdTtskSkB45CpIvHfkfhbwe/LdSiJL7sf/oZ9jz2GAAJScmMuvmHZJ091eWqRCQYYbtPYK3dHnhiqzBc++yturo6CgsLKSwspLhYU0winbWW3b/9XWuAJKYO5pSf3akAEYkAxcXFFBYWUlFRAdDlO6UdtYIHMMb8CJhJ+4aLHmCFtbbfB9bVCj66WJ+PnQ/+iorV/rcGJGdlc8pdRaSOHOlqXSLSXp+0gg8ECMCt+K88SvD3zCrA/5pckS41Nzay/f6fU/XmmwAMzD+BU352Jyl5eS5XJiK95XRg3WutfQTAGOMFbKAt/LvGmAnAP8JRnMSepvqjbLv7bmre9Q+rpY4cxSl3FZGcleVuYSLiiNMxkdYXT1lrt+O/rdUiHO8TkRjUVFfH1qKi1gBJO/VUPnPv3QoQkSjmNESMMWakMebmwM+Fxpjxgd/PCkNdEmN8hw/z6f8tovajjwDImDCBU+66k6TBg12uTERC4ShEAk9hXQWMCSy6FXjdGNNEeF6PKzHEV1vLJ7f/lEMfbwHAM3kyY376ExIHDnS5MhEJleMZ69ba+9r83ou/DcqowO0tESAQID/5KXUl/u44WVOnMupHN5OQFEqzBBGJFN1eiXTTZLFTXQVIb/cjseH4AMk+91wKFCAiMaWn21k3hOl7wrUfiRIdAmTaNEb94PsYBYhITOnpT/QsY8xS/I0Vnc1K9G8rccR3+DCf/PSO9gHy/e8pQERiULd/qq21elxXeqWpro6td9xJ3batgP8WlgJEJHapx7aETVP9Ubb+7D9bn8LKmjJFt7BEYpxCRMKi2edj2933UPvhhwBkFk5m1C0/UoCIxDiFiITMNjWx/f4HqHl3MwAZ48cz+rZFegpLJA7EVIioFXz/s9ay89cPUbVuLQBpY8cx+ie3kzBggMuViUgo+rwVfKRRK3h37Hn89+x/6ikABp08krH33k1SWprLVYlIuPTUCj6mrkSkfx14/vnWAEkZls8pdxUpQETiTK9uWhtjMoCr8TdZLAAy26z24n+vyFJglbW2Jkw1SgSqXLOG3Y/8FoDkrCw+c9edDMjWE+Ei8SaoEAm8P30x/gmHS4FFQKW1trrNZzLxt4GfCCwxxmQB91hr/xHuosVdNe+9z/YHfgFYEgelckrRHaTkD3O7LBFxQY8hYoy5Hv9LpxZ097lAoFQD24GnWrY1xsy01t4fjmLFfUd27GTb3XdjfT5MUhKjb7+N1IICt8sSEZd0GyLGmCuAZW2vOHrDWvuIMSbTGDPfWvu0owolYjRUVvLpXXfRVFcHwKjvfY+M8eN72EpEYllPbU+eCvULAgGkAIlyTfX1bPvZf9FQXg7AiddeS/Z501yuSkTcFrans4wxM9TyPTb5JxP+nMNbPwVgyKzZ5F0x3+WqRCQShG1KsbV2NbR7d0ilBtVjw94//BHvO+8A/tfajrjpRoxRc2YRCWOIABhj5uMfXAcYbYwpBKqALPyP/e4I5/dJ3zu4+jX2P+2/GznopBEU3Kp2JiJyTNj+NjDG/BVY1NXVR+ARYIkih7ZsYceDDwKQlJ7OmJ/+hKTBg12uSkQiSThnrI/u7vaV0ye8xB0N5eVs+697/I/yJiYx+tZbNRdERDroMUSMMQeNMRuMMfcYY24OzFrvzMOB21kS5ZobGth29700eqsAGHHjDaR//nMuVyUikSiY21nbrbWTe/qQtfa+wBNaIzX2Eb2stez69UOtT2INnTuX3AvnuFyViESqYG5nlQS7M2vtajcDRK3gQ1f+0stUrF4NQNqpp3Li9d90uSIRcUPYWsEbY35jrb0xzPWFnVrBh+7Qv/7Fxz++HevzkZyVzam/fEBNFUXiXDhawcfGC0ekW43V1ZTcu+RYT6wf36oAEZEeBRMiVxtjPjXGPGSMubybgXWg3WRDiRK2qYntP3+AhsqDAJz0jW+QNm6cy1WJSDQIJkQ2Au/if4/IU0BVD6FyVbiLlL5Vumw5Ne++C0D2ueeSe/FclysSkWgRzNNZm621t0LrhMFZwMzA/94AWGNMCbASWIX/ZVUSJWree499f3oSgIHDh3Pyf3xbLU1EJGjBhEjrTPPAhMEVgV/Hh8ps4EY0hhI1Gr1ett//AGBJGJBCwa2LSEzt8iEMEZEOggmRLq8sOgkVD/6rEYlwtrmZHb/4ZbsJhakjR7pblIhEnWDGRCYbY4J685C11ot/DEUiXNnzL1C9eTMA2dOmkTNzhssViUg0CiZEFgC399DypK1tIdbUKWPMw32x33h0eOs29vz+9wCk5OUx4qZvaRxERBzp8XaWtXYVgVtUxpgrjDErrbU13Xz+vjDWR+B7r8Q/7iIhaqqvZ/t997c2Vhz1o5vVmVdEHOtVK/hwvC63twLjLBImex97nPp9ewE44StfJm3sWJcrEpFo1u3trHBNHAxxPzPRYH1YVG/cRNlLLwGQftppDJt/ucsViUi062lMZFOgBXwwYyEdGGMyjTE343Cw3RgzEdjsZFtpz1dTw47//m8AEgelMvL738MkJrpclYhEu25vZwUe4b3NGHNvYE7Icmvtaz3t1BgzA//M9a3W2vtDqK/AWrtCt7RCY61l568forHK/zjvSQu/SUpenstViUgsCOrNhoEZ67fif2/6ssBLqv4aaH1yT+B/X21ZDnwB/6tyHQeIMWamtXZFb7apra3FGNPpr6KiIqelRL2qNWupWrcOAM9ZU8iZocd5RaS9oqKiTv/uLC0tBUjvarseW8F3uaExowAPkA1UAl5r7XZHO+u4bw9QGHgyrOXnTdba0d1so1bwnWj0evnopm/jq60lKSOD0379K5Iz9bp7EQlOT63ge/V0VlvhCowuLITWMRGAHCDbGHML/l5eGmgP0q7fFOOrrQX8s9IVICISTo5DpC9Za5e0/dkYUwBcefxy6V7V2nVUrVsLQNaUKWSdc47LFYlIrAlqTCQYxpjfGGPmh2t/EhpfTQ07H/oNAElp6Yz41o2alS4iYRe2EAEWA2cEBtxHhmungdnqi4ECY8ziNre4pBt7fvcYvppqAE66YSHJWVkuVyQisShst7MCYyS3BgbB7w28m/0fYdhva5dgCU7Ne+9TsXo1AJmFhWSfN83likQkVoXzSgTwd/K11t6I/z0j0s+aGxrY9atfA5CQksKIG2/QbSwR6TPhHBPZEJgzMt8Yk2GtvU9jJP2vdNly6kv9jzmf8OUva1KhiPSpHm9nGWMOAiX4+1cdBIq76OK7DH+LklnAjcaYLPzzR54OX7nSnSO7d7P/KX+PzNSCAvIunedyRSIS64IZE9lurZ3c04fatIBf3bLMGPMFp4VJ71hr2fXQw/4W7ybB/670pIh8gltEYkgwt7NKnO7cWvuu022ld6rWrqX2g/cByJ17EYNPOcXlikQkHgQTIpV9XoWEpOnIEfY8+jsAkjIyOeGrX3G5IhGJF8GEiLPmWtJvSv+8lIbKgwCceN21JKWluVyRiMSLYELkamPMp4FOvZf39G6RcL3ISoJzZPceDjz3PABpY8eRM0OHX0T6TzAhshF4F7gaeAqo6iFUrgp3kcGqq6ujsLCQwsJCiouL3SqjX+159FFskw8wnHTjDZiEsE/9EZE4VFxcTGFhIRUVFQCpXX2ux1bwxph7A+8TIfBiqln4X1k7EyjAf7urBFiJ/zHgG6y1c8LxH9Eb8dgKvnrjJj69804AcufM4eT/+LbLFYlIrAlHK/jW3uGBNx22tiE5LlRmAzeiMZR+0ezzsfu3/sH0xNRUDaaLiCuCCZGCrlZ0Eioe/Fcj0scqXvkr9Xt2A5C/YAHJHo+7BYlIXArmBvpkY8z4YHZmrfXiH0ORPuQ7dIh9f/oTACl5eQzVzHQRcUkwIbIAuN0Yc3NPT2YFbAuxJulB6dJlrW8rPPG660hITna5IhGJVz3ezgq8irblXedXGGNWdtE7q+Xz93W1TkJ3tKyMshdfBCDts5/Fc/ZUlysSkXjWq+ZK1tqn+qoQCc6+J/6E9fkAOPHr16nNu4i4SpMKokjdjh0cfO11wP/O9LRx41yuSETinUIkiuz9/R8Bi0lIZPg1X3O7HBERhUi0qP3wQ6o3bgAgZ+YMBp54ossViYgoRKKCtZa9j/8BgIQBKZzw5S+5XJGIiJ9CJArUbNrMoY+3ADB03iUMyMlxuSIRET+FSISz1rLvCf/EwsRBqQy7Qq+tF5HIoRCJcNXvrOfw1k8ByPvipSSlp7tckYjIMTEVIrHWCt42Nx+7CkkdzNDLvuhyRSISL8LWCj5axGIr+Kq169i2eDEAw7/yFfL/7WqXKxKReNNTK/iYuhKJJba5mX1PPglAUlq6miyKSERSiEQo71tvc2TXLgDy5l9OYmqXV5MiIq5RiEQgay2ly5YBgauQi+e6XJGISOcUIhGoZuMm6kpKABh66TxdhYhIxFKIRBhrLaVL/VchiampGgsRkYimEIkwte+/f2x2+ty5JA0e7HJFIiJdU4hEmNKlywF/jyzNCxGRSKcQiSCHtmyh9oP3Aci96EKSMzNdrkhEpHsKkQhy4OlnADBJSeRdfpm7xYiIBEEhEiHq9+6l6q23Acg57zx16hWRqKAQiRAHnn0O8LegyZt/ubvFiIgESSESARqrqji4ajUAnsmTGTRihMsViYgERyESAcpe/AvNvkYA8vS+EBGJIjEVItHYCr7pyBHK//ISAGljx5H22c+6XJGIiFrBR42yF15kVyDwRv/4NrKmTHG5IhGRY9QKPoLZ5mYOPP8CACnD8vGccYbLFYmI9I5CxEXVGzZwdH8pAHmXzsMkJrpckYhI7yhEXHTgOf9VSGJqKjkzprtcjYhI7ylEXFJXUtLa4mTI7Nlq9y4iUUkh4pIDzz0PgDEJDL3kYperERFxRiHigsbKKirfeAMAz9QppOTluVyRiIgzChEXlL/yCtbnAyDv0ktdrkZExLkktwvoijHGAywAPMBowGutXeRmTeHQ7PNR/sorAKSOHsPgU8e5XJGIiHMRGyLAAmtt67RzY8xiY8xKa+0sN4sKlXfdmzRWVQEw9JKLMca4XJGIiHMRGSLGmAL8VyBt3QNUGWM81lpvvxcVJmUvvAhAUkYG2dPOdbkaEZHQRPKYyG1tf2gTHAX9X0p4HN66tfX96blz5pAwYIDLFYmIhCYir0SstSVAVttlgasTrLWbXSkqDFquQkxCIrlzL3K5GhGR0EXylcjxFgFLuvtAbW0txphOfxUVFfVPlV1orK6m6o01AHimnMWAIUNcrUdEpK2ioqJO/+4sLS0FSO9qu6jo4muMmQgs7m5QPdK7+JYuW87eP/4RgLF330365053uSIRkZ7FShff26L5qSzr81H+8ssADDp5JGmnn+ZyRSIi4RHxIWKMWQxc73YdofBu2EiD/8UuDL14rh7rFZGYEdEhYoy5Bbin5cksY4ynZYA9mpT/5S8AJKYOJvv881yuRkQkfCI2RIwxM4EVx80JWQBUulORM0d276HmvfcAyJkxncRBg1yuSEQkfCLyEd/A1cbKwO/brvK2ncUeDcpfeqn190MvnutiJSIi4ReRIRKYJxL1AwdNR45w8LXXAMiYMIGBw4e7XJGISHhF7O2sWFD5+t9oqqsDYOjFemeIiMQehUgfsdZSFriVNSA3l8zJhS5XJCISfgqRPnLogw85snMnALkXXYRJTHS5IhGR8FOI9JGyF/19shKSkhkye7bL1YiI9A2FSB84WlaO9+13AMg+bxrJmRkuVyQi0jcUIn2g/KWXsLYZgNxLLnG5GhGRvqMQCbPmhgYqXn0VgLRx4xg8ZrTLFYmI9J2YCpG6ujoKCwspLCykuNidOYmVb7yBr7YWgKHz5rlSg4hIqIqLiyksLKTC3/cvtavPRUUr+GBEQit4ay3/+v4PqNu2jeSsbD73u0dJSIrI+ZwiIkGJlVbwUaH2/Q+o27YNgNyLLlSAiEjMU4iE0YFnnwUgIXkAuXPVJ0tEYp9CJEyO7NxJ9caNAOTMnKHHekUkLihEwuTAM88GfmfIu+yLbpYiItJvFCJh0FBZycG//x2ArClnMfCETsefRERijkIkDMqefwHr8wGQN/9yl6sREek/CpEQ+Q4fpvzlVwBIO/VU0saNc7kiEZH+oxAJ0YFnnqWp7jAAw+bPd7kaEZH+pRAJQaPXy4FnnwNg8CmnkHnmGS5XJCLSvxQiIShdtpzmo/UADL/ma8e/D15EJOYpRBw6WlZG+csvA5Dx+fFkTJjgbkEiIi5QiDhU+uSfW5/IOuGar7pcjYiIOxQiDtSVlHBw9WsAeM48k7SxY12uSETEHTEVIv3RCr65sZHtD/wCa5sxCYkM/6quQkQk9qgVfB/Z89jj7H/6aQBO+NKXOOHLX+rT7xMRcZNawYdR7Uf/ZP/TzwCQOnoMwxZc5XJFIiLuUogEyXf4MDt++UvAkpCUzKgffl/vCxGRuKcQCUJjdTWf3P4Tju7fD8Dw665l0EknuVyViIj79E/pHjRUVPDJT/8v9Xv2AOA5awpD513iclUiIpFBIdIF29xM9cZN7PrNwzSUlwGQc8EFjPzudzAJuoATEQGFSDtNR47QUFbGoS1bOPDs89Tv2d26bujFF3PSwusVICIibShEgC2LbqN+9y58tbUd1iUOSiX/6gXkzb9cvbFERI6jEAF81d4OATIgZwhDL53HkDmzSRo82KXKREQim0IEyJo6FV9NDQPyhjIgdygpw/JIHTNGj/CKiPRAf0vib+MuIiK9p1FiERFxTCEiIiKOKURERMSxmAoRp63gi4qK+q6oGKbj5pyOnTM6bs719tipFXzvtiVWjkN/0nFzTsfOGR0355weO7WCFxGRPqMQERERxxQiYRDKq3hDfY2v29u7+d3RfNzd/O5oPu6hiub/drePXZestTHxC9iXn59vnfAfBucmTZrkyrZub+/mcQt1e7ePeyjHzu3a3TzuOuecc3rs8vPzLbDPdvF3bywNrDcmJCQk5eXl9Xrb0tJS8vPzHX93RUUFQ4YM6fdt3d7ezeMW6vZuH/dQjp3btbt53HXO9f+xO3DgAM3NzT5rbXJn62MpROqAZKDcwebpQMcWvsFLBepc2Nbt7d08bqFu7/ZxD+XYuV27m8dd55xzTo9dLtBore30Md+YCREREel/GlgXERHHFCIiIuKYQkRERBxTiIiIiGNx/VIqY8xCoDLwY7a1NkJn80QOY8yVwCxgcWDRlcBma+0q96qKTMYYD7AQGG2tvaGT9Tr/OtHdcdP5173AsVsAeIDRgNdau+i4z4T1vIvbK5HAgSyx1q6w1q4AKgPLpHvZQCGwDViJ/yTVH+DjGGMmAjMBbxfrdf51oqfjhs6/niyw1hZba5e0BLAxZmXLyr447+L2EV9jzCZr7aSelkl7xpgrAyefBKHlX86d/Ita5183ujluOv+6YIwpAK601i5ps8wDVAFZ1lpvX5x3cXklEjiwBZ2smhhYJ9JndP5JH7qt7Q/WWm/gtwV9dd7F65hIIcfuCbblxX+QN/drNVHGGDOzzY8T2/7LR4Ki8y8EOv86Z60tAbLaLgtcnWCt3Rw4bmE/7+I1RDx0fs+1Ev89V+laCf770CXg/1e1MWbx8YN30i0POv+c0vnXO4uAlpD10AfnXVzezhLnrLWbW/4AB35eAdziYkkSR3T+BS/wkEJBXwdsvIaIF38qHy+bzi/3pActl80SFC86/8JK51+nbrPWzmrzs5c+OO/iNUQ20vnlmwf/5bJ0whhTYIzZ1skqb3/XEuV0/jmg8y94xpjFwPXHLe6T8y4uQyTwxEJnyVvS5mkG6dzitj8EnurwtL3FIN3T+RcSnX89MMbcAtzTci4Fxo0K+uq8i8sQCXi47SSbwO8Xd/P5uNfFH9TbgA6zsaWdzv71p/OvZ+2Om86/ngWewFpxXCgs4Fh4hP28i9vJhnBs9ib+yzm1nQhCm7YKEHjaQ8eto5aJX8DV+B+fLAZWtp1drfOvo56Om86/rgWOXae3+6y1WW0+F9bzLq5DREREQhPPt7NERCREChEREXFMISIiIo4pRERExDGFiIiIOKYQERERxxQiIiLimEJEREQcU4iIiIhjChEREXFMISIiIo4pRERExDGFiEgfMsZUGWOWB1p0u1XD4kANy92qQWKXuvhK3Av8Bb8c/5vfSvC/KW8iUAgsC3wsG39r8oK2bbWD2Pdya+1VYapzIvAIgLV2koPtw1aLSIsktwsQiQA3ADOstZtbFgT+1V5prW33wiNjzKb+Lq6FtXazMeZ6/IEnEhF0O0sENrQNkICZwKpOPtvZsv7kdfn7RdpRiEhcM8Zcif/teW2XefC/9W1lJ5sc7PuqRKKHQkTi3arj3kcN/qsQ6PyqQ69iFWlDYyIS1zoJEIBZQEln67r4fK8Ern7AP1jvsdYuabNuIrA48OMi/IP5NwBXtf3uwOda9jERWGGtLQm1NpHe0pWISEddjYeEzBhzC7DZWrvCWlvcZhngHzzHHx7ZQLa1dgUdb6sVAF5r7WZr7apACC03xhT0Rc0i3VGIiLQRGA8poPPxkHAYDVzZ5udV+K982vICE621qwCstUuOuwIq6eSqYynHrmBE+o1uZ4m01914SMjaPjIcuHIoxH/Vcbze3praDNwWQmkijuhKRKS9LsdDwsEY4zHGPGyMWRhYpHEMiWoKEZH2uh0PMcbcYozZZoxZ2WZZb1qabAcettYWB25JVbbZTyhjGhNxfw6LxCGFiEhAT+MhgbAosdaOBha1ecoq2P0X4H8aq+3ExgL8c1LAHwTByO4kcK7GPyDfcrWzsONmIuGn3lkS9wJPR43GPz7R8i/6EmB5y+B24HMFbQe0W65A2n6mk32361cV+K4c/EFVGfie24Bt+Ht3Efj5SmAJsPS4diwe/FdL3sD2HR7xbekFdnyPL/XOkr6gEBFxyEmIuCmSapHYodtZIs55aDOmIRKPFCIizmV30rhRJK4oRESc01WIxD2FiIhznU0S7CBS3mwIbHCrBoldGlgXcSDwlFS2mh5KvFOIiIiIY7qdJSIijilERETEMYWIiIg4phARERHHFCIiIuKYQkRERBz7/4dCfncVf+W1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 443.077x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "T_gamma = torch.linspace(0.1, 20, 100)\n",
    "plt.plot(T_gamma, T_nu(T_gamma))\n",
    "plt.xlabel(r\"$T_\\gamma$\\,[arb.]\")\n",
    "plt.ylabel(r\"$T_\\nu(T_\\gamma)$\\,[arb.]\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e03d8360-f04c-42da-84da-009200016d96",
   "metadata": {},
   "source": [
    "Draw samples from latent function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "5da3258b-2159-4155-90a0-f0a792ca5dc9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define mean and variance of observed sample\n",
    "loc = torch.Tensor(obs(T_nu(T_gamma), T_gamma))\n",
    "scale = 1.\n",
    "\n",
    "# Draw sample of observations\n",
    "obs_sample = torch.distributions.Normal(loc=loc, scale=scale).sample()\n",
    "\n",
    "gpr = GPRegressionBBN(obs_sample, num_inducing=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "867f9283-65c6-4240-af52-c4828c17c0f1",
   "metadata": {},
   "source": [
    "Train with SGD."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "bb1d408d-73eb-4591-a6af-6172f87ccb66",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bae78e510516411795c9ba370de24e66",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/1000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "num_particles = 4\n",
    "num_iter = 1000\n",
    "\n",
    "def train():\n",
    "    pyro.clear_param_store()\n",
    "    optimizer = pyro.optim.Adam({\"lr\": 0.1})\n",
    "    elbo = pyro.infer.Trace_ELBO(num_particles=num_particles, vectorize_particles=True, retain_graph=True)\n",
    "    svi = pyro.infer.SVI(gpr.model, gpr.guide, optimizer, elbo)\n",
    "\n",
    "    gpr.train()\n",
    "    iterator = tqdm(range(num_iter))\n",
    "    for i in iterator:\n",
    "        gpr.zero_grad()\n",
    "        loss = svi.step(T_gamma)\n",
    "        iterator.set_postfix(loss=loss, lengthscale=gpr.covar_module.base_kernel.lengthscale.item())\n",
    "\n",
    "train()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f71e9da4-9f9f-4efb-a1c9-209f13a74f16",
   "metadata": {},
   "source": [
    "Draw samples from the conditioned GP."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "0c72e29d-e968-47be-bf84-51476f224dba",
   "metadata": {},
   "outputs": [],
   "source": [
    "gpr.eval()\n",
    "with torch.no_grad():\n",
    "    gpr_samples = gpr(T_gamma)(torch.Size([500])).abs().detach().numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "f689725f-fa42-4e01-85ac-fa294cb548e6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fdcb93b4400>"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAFRCAYAAACmHtAXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABUaUlEQVR4nO3dd3Rc133o++8+Z3oBBoVdpEiQVBclgVAcx3GJRbrEN8U2JafYN9VkbrLeXy9PtG7eWmb+eFch7ZvkvmTFJuyUd53EkUg7znXs2CbllpvEtiTalixZjWARO1EGmHrqfn+cGRAgBn0GgwF+n7W4RMycsnU0mh/23r/920prjRBCCNEIRrMbIIQQYuWSICOEEKJhJMgIIYRoGAkyQgghGkaCjBBCiIYJNbsBS0UpVQTCwPVmt0UIIVaQNYCjtU7UelOtlhRmpZRjGEZo3bp18z43l8uRTqcXfO9isUgiUfP5N/TcZp/fzOe22POb/dwX8+ya3fZmPnf5zC39s7t69Sq+77ta63DNA7TWq+IPcGnDhg16IYLHtHC7d+9uyrnNPr+Zz22x5zf7uS/m2TW77c187vKZW7iFPrsNGzZo4JKe5rtX5mSEEEI0jASZJbB///6mnLsczm/mvVv5uTfz3q383Berlf/dm/3sprOa5mQubdiwYcOlS5cWci6r5TnVkzy3hZNntzDy3BZuoc9u48aNXL58+bLWemOt96UnI4QQomEkyAghhGgYCTJz8NGPfrTZTWhJ8twWTp7dwshzW7hGPbtVNSfT3t6+YceOHUAwSbZcJ8qEEGK56+/vp7+/n+eeew7HcUa11plax62qILPQiX8hhBC1ycS/EEKIppEgI4QQomFWTYFMIYQQk3mej+P5WI5LLBIiGq5/SJAgI4QQq0BQSwxcz6dkOxQtB8/XgMbzNJ3p+OoKMkqpDLAf2K61PlDj/f3AcOXHTq11/xI2TwghlhXf13i+j+drPM/Hdj0cz8d1PXwADQrQSmOgCIdMIiEFgGW7DWvXsgwySqleoAfITvP+fmBAa32y8vM+pdR+CTRCiNXA9XxKlkPZdvF8je/7aA1aAWgUCsNQGEoRDpsYSjWtrcsyyGitTwGnlFL7pjnkgNZ694TjjyulngUkyEzj5MmTPPzww+zfv5+uri4AHn/88fGfh4aGOH78OAcPHpT1Q0I0ieN62I5H0XIwDEXIMAiZBspQaD9YblKyHUq2h6EgZBqYhiJkmqgmBpKZLMsgM5PKMFpPjbd6lVIZrXV2aVvUGrLZLE899RS9vb3jrx08eJADBw7Q0xM8zgMHDnD8+PFmNVGIVUNrjeP6lG0Xx/PwfD0+rKUIgofvaWzHw9caXemdAJiGIhFtna/u1mnpDX3cmIuZKEsQfE4taWtayMQAU9XZ2Tn+956envGAI4Sor2pgKdkOhbKDr30MFGqmYS2zOW2tp1YMMhlqz9UMA501XhdAJpOp63FCiOl5no+v9fgkvOV6lCwHzY1Jd8NYARFkDlbVYsxcLodSquafQ4cONbt5DbVnz545H7t7926OHDkyPkdz/PhxTp06xe7duzl48OD4cQcOHKCjo2PK+UeOHOHkyZMcOXKEgYGBurRfiOXK84Nhr9FCmevZApcGx7g0nOdqNs/gWIGRQgnLcYmGTeKRMNFICMNYnvMnMzl06FDN787Lly8DpKc7b1nXLqtM/O+dmMKslNoDHNVab7/p2BHgoUrSQK1rLbh22fn+T1E6c2be5zVKfNs2tuz/8KKvo5RiZGRkSu+lv7+fEydOcOzYMU6dCh5nb28vx48f5+mnn+bw4cPjx3Z0dDAyMjL+8969ezl27Nj4NXfv3s2zzz676LYKsZz4vsZyXHIlG8f10BoMBWZlIn65TsJPx7Jd0oko6UR03ufOVrusFYfLnqH2sFgGaMivzaUzZ8j96EeNuPSyNHGeZuI8zmxDadWANPG4vr4+Tp48Oa+elBDLjdYa1wvWoJQth4Ll4mufSMgkFmnFr9Gl03JPR2udVUrVmvgfaFRmWXzbtkZcdsGWoj0LSQA4efIkcCPYAAwPD8uQmWhZnu9TKDnkSxbB2vigxxINmyi1OuZUFqsVgkytXsvRiYsvK4szD9c4ri7qMTTVaqpraeZjaGiI3t7eSb2fY8eO1bNZQiwJz/cplh3GSjZoTaTJCxpb2bKc+FdK9SilHgUeA/YopQ5X5mIA0FofqRy3p7pgU1b7L41sNjvtew8++KD0WkTL0jpYl5LNlbgynGe0aBEJGcQiIQkwi7AsezJa6wHgSOXPdMdIUKmD4eHhmnMtQ0NDU17r6elhePjGSOXEYTGAffv28fjjj5PNZsevWQ06sv5GLGeu55PNlyjZLqahKsNhEljqYVkGGdFY/f394xlfBw8e5MEHH+TRRx8FgsBx9OhRBgYG2L59+6QSM9XFmtW5l2qCwIEDBzh69CgATz31FI8//jgPPvjg+Dm1FoEKsVwUyzYjeQulNIlouNnNmZXva2zXC+qWaR9TVTLaDANTBTXLglIzy2OgalmnMNeTbL8shJjI9XzGChYFyyYWXj5rV1xPU7JsSrZHrlhmrGhhuR6+p/F1EGCmUARZCRUaCBkG6XiEVDxCOh4hGg0RC4cIm+Z4IKqSFGYhhKgT39cUyjajBQvDoGm9F601tutTth2KlstY0WI4X6JUtgGFRmMqg2jYxDAVoZDCUAaJaGhOQ3me71fW8lhc8DUaKtXPgn1lDMMgGjIr5Wzgjs3dCwoys5EgI4RYNWzHYyhXxPM10cjSZYxVA9tIocxIrkzRcijZDuhqB0QTNk0iYZNMMlaX+SDTMDAjxrTreHxf42kf1/MZyZXY0Dntov1FkSAjhFjxtNbkSzbZQplIyCQSaewal7LtUrBs8sUbgcXXPoogqSAUMsgkok1NLjAMhYFJ2ISQ2bh2rKogUywW6evrA2D//v2yb4oQq4Dn+QznS+P72Dei9+L7mlzJYjhX4spIgZJlB8X5FcRCIdLxMKaxPCbi66W/v5/+/n4GBwcBEtMdJxP/QogVy3JcBsdKKHRD9q8v2S7XswVevz5GuZL+HIuGiIZaqxrA9dEC92xdx73b1s37XJn4F0KsOjcPj4XM+n3pe74mWyhx4fpYEMCUJhmN0JmO1e0eK4kEGSHEiuJ6PiMNGB5zPZ8rw3nOXM0G1w6H6Eg2d16lFUiQEUKsGIWyTTZfxlCKeKQ+qcmu53NxcIwzV7J42g/WnsTidbn2aiBBRgjR8iaWhanXwkqtNddHi7z0+hCu55FORAitsMn7pSBBRgjRsrTWFC2HbJ3LwowVLE5fHmYoV6qsmF/+5WaWKwkyq9DAwABHjx6lq6uLTCZDZ2cnvb29dHZ28swzz7Bnzx6y2Sz9/f0cPHiQffv2sXfv3vEKzKdPn+bgwYNS9FI0ldaabL5MvlyfsjC+rxnOlzh7JUu2UCIaCtGVlmGxxZIgM0f/9sJ5snmr2c0gk4rypru3LPj848eP88QTT/CpT31qUvXlU6dOceDAAQ4cCHa6zmQyPProozzxxBMcOHBg0s6W2WyWbdu2cebMmVl3y2yk48ePs2/fvqbdXzSP72uGc0XKjrfo3ovna4bGipy+NELBsklEwnSlp132IeZJgswcZfMW3e3N/61mcLS04HMHBgY4ePAgzz777JTg0NvbO+eeSSaTYc+ePTz++OMcPtywveJmJXvXrE6e5zM0VsLxPOKL2PrY9zWDo0VeuzxMyXJIxCLSc2kACTKryIEDBzh48OC0vY8DBw7M+Yu7s7OzqV/yBw4caGovSjSH5bgMjZUAPW1NrrnI5su8fHGIXNEiFY/QuZqDi2MTPT+A547BAhZjzkaCzCpy8uTJ8X1faunt7Z3zF/fJkydn7MWcPHmSgwcPsmfPHvbu3QsEQ3K9vb2Tht4gKE/R09MzPuczcQjs+PHj45ulnThxggMHDnDq1CmGh4cZGBigvz/Yu05KBK1sWmtyRZvR4uIWV5Ysl1cvDXE1m68Mi63C4JLPYZ4bwDg7gHF+AOPiBaKui3PPvbD3p+t+u1UVZFZz7bK57lA52/vVpIEDBw7MOB+yZ88eDhw4wLFjx8aDUTXgVBMNAB5++OFJ80P9/f3jcy1Hjhxh//79U3bZrN736aefXlX/DVcrz/cZzgWLK+ORuZW5v5nraS4MjjJweRjTMOlKxVfHIkqtUUODGGdPY549jXFuADV0veah3pnTaM9DzTGAz7V22aoKMolEgmeeeabZzWiK6i6WE7dGnqsTJ06M9zIymQyPPfbYnK4xMZhUVYfsTpw4Mb7D5sRr7d+/n+3bt48Hkv7+/vFdO/fs2TP+7yFWB9fzGRoLSvMvdHHlcL7Ej88NUnYc2hPRFVeochLfR12+GASUalDJ52ofaxj4Gzfj37qNke6NbHtj35wDDNz4Rb1Su6w43XGrKsisZplMhkwmw8DAwJQv/pMnT3LixAn6+/vp6+vj8OHDk47Zu3fvlCGuhert7R0P9KdOnarZcxoeHiabzbJ//34efvhhHn/8cfbs2cNjjz0madOriO14XB8rYiiIhuc/PGY5HqcvD3NxaIxULEpnagUOjbkuxsXzGGdewzxzGuP8AFjTZMFGY/hbtuJt3Y5/aw/+LbdCJAJAabSA0Z5pSBMlyKwie/bs4eTJk1OCzJ49e9izZw+nTp3i4YcfnvJ+MwwPD9PT0zPei3ryySd5+OGHOXbs2JT2LaR3Jpa3kuUwNFYiHDLmvVe972uujhR45eIgvmZlDY05Nsb5sxhnTmOefQ3j/Blw3ZqH6nQbfiWgeFu3o9dvhCb04iTIrCKf+tSn2LZtG/v27Zu2B9Ho4ahTp06N94p6e3unJCJUh+V6eno4ePAghw8fJpPJsH//fvr6+moGyYGBAXp6eiTQrBCFks1wvrSgBZajhUrWWMEinYgQrmP15aawbYzzZzDPvIZx5jWM82fB92oeqju78bdux9u2HX/rdnRnNyyD4CpBZhXJZDIcO3aMAwcOcPTo0UmBplHpyNV5l6qjR4+OB5Y9e/Zw+PDhST2RJ598cjxRIJvNcvLkyUlDddXjqr0cCIJjNfBks9nxXpBoLVprxooWY0Vr3tWTXc/n9OURXr+eJR4Jt25KsmNjnKsElYFXMV4/N31QWbMOb9sO/G078LbtgLb2JW7s3EiQWWX27NlDX18fBw8eZPv27eNlZQCeeuqpScHmyJEjnDp1alJ22HxVex/A+LUmBoATJ05w5MiRSa9VM8a2b98O3AhUAwMD4+9V062PHz8+qQfT39/P6dOnZ0zVFstPtTx/2Z5/BtlwvsSLZ69jex6drTY05rkY589hDLyCOfAqxrkz0weVdRuDXkrPTryt2yGVXuLGLozsjDlHX/ruq8tmxf973rCz2c2Yk+PHj/P00083tSqAWP5KlsNwroyhIDKPCf6S7TJweYRLw8HEfmwByQFLSWuNZbt4588RPvMK0bOnibx+BuU4NY931m7A3rode0sP1pZt6ESKcMgI1gkZRl0qTVfJzpjLQCYVXVRJl3q2Q4iVIle0yBbKRMPmnFOLXc/n9eujnLmSxTDUsp3Y931N2XZwr1whevZVomdeJXN+ANO2UAoUKpgyqQTHyKZNxO66m+iddxG74y6MtjSer/F9jeW4FMoOowWLsVKZ0ZIFOghcpmGQjIeX7TYEEmTmaDFFKYUQU1UDzHzmX3JFi+fPXqNkL881L47nUbo+TOTMK0TPvkbXudcwc2MopW7MwVd6IOF164nfdRfxe+4hdtddhDIdU64XqnTO4tEQmVSMTd3BEJnWGsfzKZRtrmWLXB0pYLsW0bBJMhpeVkFXgoxoiFOnTnH06FEGBgZ48MEHpVqymGS+Acb3NReHxnj5whCJSHhZrXnxymXsl18mPPAyibOvkRm8ilEjqIQyGWJ330P8nnuI330P4TVrFnxPpRSRkEkkFacjFWfnxk5Gi2XOXxtjcLSIYSiSsdCyyK6TICMaore3dzz7S4gq39eMFcvkSvacJ/jLtstLrw8yOFYkk1wGvRffR125hP/SCxivvETi9TOktcY0KoGlElSMWIz4XXcHQeWeewnfckvDehiGoeioBJxC2eHqSI5LwwVypRIKRSoeblrAWVVBZjXXLhOi2VzPZ3ishF0p0T+XL9xr2QIvnruOUqq5xSxzY5ivvoR65ccYr76EUcijlMI0VWUCXqEMk+iOHcTvvZf4PfcS27EDFVr6r9hkLEzPhk62re+gUHYYzpU4cyVLzrdJxSJEQvUJNnOtXSbZZUKIhiuWbUby1pwzyBzX57VLw1wYHKUtGSWy1L+FV1KLzVdexHjlRdSlC6ABBaYKMruUCuZVErt2Ed+1i/hdd2MkludmZ66nuTaSZ+DKCGXHJRmLTMrGk+wyIURL8nyf0bxFwZr7Fskj+RIvVNa9dKWXLnNMjWYxXv0x5ssvYr72MlhltNZoHQxHmaaBmUgQv/tuErvuI75rF+F19d9/pRFCpmJjd5p1HSmGckVOXx5hKFeaEmwacu+GXl0IsWqVbZehXAkFc9oi2fU0Z6+McObqCKl4lI7Y4rZVnpXnYbx+NggqL7+AujJhlKOSHqyUIr69h7beB4jfdx+xHTubMgRWL6apWJtJ0t2WYKRQ4oWzg1h27dpn9dK6T0sIsWxV64/Ndf1Lrmjxo3PXKVk2nal4XRcaTlIsBEHlpRcwX/0xlG9a+6bBi8Wxd9xBR99uNv7UG4h0ZBrTliYyDEVXOsEb7tjEj89f5+yVbMPuJUFGCFE3WmtGC2VyRZtYdPb0ZM/XXBwc45WLQ8QjITrqnZqsNerq5SCovPQjjNfPQo15aG/jZvJbd8Kd97K5917WdrYRMpfPWpNGiYZNdm1bRyoWIT6H3uZCtHSQUUpNTA/LaK2PNK0xQqxyXqX+WMl2iUdnzx7LFS1efP06+aJd39Rkz8UYeA3zpR9h/vhHqOzw1GOiMbydd+DsvJORzdsx2zvYubGTtR1JzEb1opYpw1BsXtNGOtGYaiItG2SUUo9ODCpKqYxS6rDW+mAz2yXEamQ5LkNjwdDTbPMvrudz5soI566OEo+G6lMxuVQMhsF+/DzmKz8GqzzlEN29Fu+Oe/Buvwtny1bGbB+A7es72NTdNu99a+bK9Xw8z8cn6EEpVOVvesrfq6qLOauVApq+NmgRWjbIANsn/qC1ziqlMk1qixCrktaaXNFmtFgOCjfO8kWdK1r86Ox1SrZDRyq2qLkXNTIcBJUXn8c4+xr4/k0HGPjbtgeB5Y570N1rgsWgJRtt+dy6tp1b1rQvaNfN2fhaYzsevtZEQiapeJRQyBgPHloHx2hfY5rB64ah8H2N6/s4jofr+/i+xvM1juuiCeqUhQxjPAC1glYOMn1KqV6t9almN0SI1cjzfUZyJcrO7OX5fV9zoTL3koiE6UjF5n9DrVFXL2G+8ByhF59HXb4w9ZhYHG/nnXh33Yt3250QT4zfP1ew8HyfW7rb2bKunXik/l9/nu9jOR6GYZCOR4lHQ4Tns/jRhAgm3NQb9Hwfxw1qldmuh/aD1GofTcg06rbAshFaOcgcBJ5VSh3UWh+pzM/IUJkQS8BxPQbHimgN8cjMw2O5osUrF4YYzpfpSM1z7sX3MS6cw3zhOcwXfogaHpxyiG7vCILKnffib9sO5o2vNa01uZKD63lBcFnbTjxa36+9arFKz9OYpqIzHSceCdc1Q840DMyIQSwy+d/NcX3GimVKtouhIBwy57XZ21Jo2SCjtT6plNpNEGgeA3ZrrbNNbpYQK16hbDOSKxMOGYRC0wcMxw3mXs5fHyUaNulum+Pci+dhnD2N+cIPg8CSG5tyiF63EffuXXh37UJv2DRlm2GtNYWyg+W6bOhsY+u6DMk6r7vxfB/b9VFALBIimYoQDZtLt3hUKSJhk+72JLbjUSzbFC0Hn2BOJxIylsWQWssGGaVUD7AH6AAOA6eVUnu11ienOyeXy0370D/60Y9y6NChRjRViBVh0vbIs6zeH86X+NGZ67i+R0dyDnMv1Yyw579P6MXnoFi46QCFv2Ur3l278O7ehe6avoJx0XIoWS7d7Qnu37C+rllT1d6Dp31Mw6AjFSMWCTV9Yj4SNomE47SnYrieT75kky/b897GeiaHDh3iD//wD6d7e9ptOlu2dplS6qjW+sCEn/cBnwK21erRSO0yIRZuUnryDPMvrqc5c2WYc9dGZy9Z4rkYp18h9Pz3MV98HkrFye8rA79nB9499+Peee+se9hbrke+ZNOWjLJzY2fd1txorbFdD98HBcRjYZKxMJHQ0vVaFmJ8O4U5lPOxbJd0IrqggLwia5dVAsqkOvJa6+MTejfHm9IwIVag8fRkrWdMTy5aDs+fvUqhZE/fe/E8jIFXCT13CvPF56YGFsPE23kH3j334d15LySSs7ZPa022aBEyDHZtW0d3W6Ju8yG26+F6mnQsTDwWlMtvWDWCOksnooRMg6GxEiFTzS8BoY5aMsjMQDLNhKiTqenJ039JXcsWeOHcdUKmmtqD8H2MM5WhsB/9YOpQmGHi3XYn3j33B4ElPvceSNl2yVs2m7vb2b6hk/AMc0TzUc0Si0VCdLfFmvYFvVjxaJh1HQaDY0Usx2tIuvZsWjXInCQYGru5x7JXFmMKsXiO6zGSK2O5Mw+PeZ5m4MowZ69mJ5fk1xrj9XOYzz2L+fz3p07eGybebXfg3duLd8c98wosALbnkSvaJGMRdu/cWNedMi3bRStFd1uiYaVWllI4ZLI2k7xRDbuO8zRz0ZJBprLw8qBS6lEgW3k5AzzetEYJsQJorcmXbLIFi5CpZhwey5dsXjh/jXzpRlFLde0KoR88g/nDZ1EjQ5NPMAz8HXfg3vsA3l275h1YAFzfJ1e0CJkh7t6ytq5lYHxfB3utRCO0J6OYDaoA0AymYdCRjhEJmWQLZUxTLdnampYMMgBa6wFAapUJUSeW45LNl7Fdb8bfdn1fc2k4x8uvDxINh+jyLMz//e+EfvgM6vLFm45W+D07cHf14t1z/5zmWKa752jJQqHYsbGLjV3pupWBcT0fx/MwlEFnOk4yFqnLdZcbpRSpRIRoxGQkV6ZoOUvSq2nZICOEqA/X8xkrBEMpkZA5Y++lZLm8/Pp1hgaH6TrzMpEfPosx8CowOUvVv+VWvPt6ce/tnTUrbDbVlOTNa9rZuj5Tl3mF8YwxrQmbJp2pOLE6L6BcrsIhkzWZxHiPNdzgatMSZIRYpXw/GBobK1oYxsyFLX1fc+l6lnP/+zsknz/FpldeANeZdIzuWoN7fx/e/X0zrmOZK8/3GS1aJGMR3nDHprqsdxkv+6IMktEgY2y5pyI3glKKdCJKLBJiJFem5NrTL3RZJAkyQqxCZdtlJF/C8zXRyMylSAZfepWLXzuBeeoZugo5mHhoMo276wG8B34Cf9PmKSvvF2JiKZieDR1sWZPBXORv28HqfA/TMBpS9qVVVXs18VKoYVWoJcgIsYo4rsdoIah1FQmZRCK1h568fJ7hb32Loae+jn/+HJFq1V8FhEJ4d+7CfeBB/J23T6oVtli5ko3teGzoSrNtfWZO2zbPxNcay/EwDUVnKk48Gl51vZbZVHs1jSJBRohVwKtkZeVKNiHTqPnlrT2P0vPPM/rNb5D73vfwbAdUUAcLBf7W7bi9PxFM4Mfqu4Ol7XrkihbdbQm2b++sy5ee7Xh4vqY9GSUZi0jPpUlWVZApFov09fUBsH//fvbv3z/LGUK0PstxGc6V8H1dc82Lc+0auW9+g9y3vok9OITr+6CD33B1RydO7xvwHngQ3dVd97ZVs8ZMZbCrZz1r2hOL7mnYlb1Y4pEQ3cnWXUi53PX399Pf38/g4CBAYrrjWrZ22XxJ7TKx2vi+ZqxYJleyp2woph2HwjNPM/bUU5Re+BFQ2cHR16hIBP+e+3F3vwF/2w5oQPHHifMum9e0s219x6JX61uOi+9rEtEI6UREgssSWZG1y4QQM5uu92JfvEjuG18n961v4eVzQDBv4Xo+zsYt8BNvxNu1e0ELJecqXw7mXdZ3putSgj8YFvMluCxTEmSEWEF8X5MrBeX4qxP72nHIP/09xk6epPTjF8eP1YAbiTF21/3wkz9NeNMtDW2b7XqMFW06UzHu71lcCf4b2xtDMhqW4LKMSZARYgXwPJ+i5TBWskEHvRf32jWGnjpJ7pvfxLupdpjevpNrd/bi3rWLVGrx8yAzcX2fsaJF2DTZtXUdazILr5LsuMF8i6GC7Y0TsXDDUm9FfUiQEaKFVRdU5koWAGEDys/9kCsnvkbxhz+YdKyZSmO+8ae4sOM+im0dtCeixBu42ZbjeeRKNqZhsHMRpWD0eK9FE4uEyKTiREKtU3J/tZMgI0QL0lpTshyyBQutNaFykfy3vsm1Eydwrl+bdGzsjjuJvOVtXNq0k+uFMslYlM4Glnyv7u+iUNy+qZv1nakF9zZs18PzfJKxCKm4DIm1IgkyQrQY2/HIFspYjgsXXyf/ta+S//d/Rzv2+DFGLEb6LW8l8paf4UI4xaWhHOGyQ2cq3tChsaLlULIdNndnFlxnbOIWx9FwiO62hASXFiZBRogWUS1kmS+WcH9wivyJr1F+6ceTjonccgtte99B/KfezMW8xcDlEUyzQGcq1tDgEtQZs0nGwjx42ybak7F5X6Naat9QBvFoqCW2OBazkyAjxDJXnXfJDg5T+vY3KJw8gRssgANAGQaJ3X20v/OdxO66m5F8mafPDVJ2nGBflAbOu0C19+KyfYF1xqqZYkopOtNxYpFQw9sslo4EGSGWKV2pu3Xt9BlyX/kKpX/9Nr5VHn/fTCZJv/0h2va+g/CaNZQslxfOXefycI5UPFrX3SJrKTsehbJNOhHhDdsWViXZcly0hraElH5ZqWYMMkqpJxZ5fUWQjv+01vrji7yWEKuC5/mUbZfrP3yOsS99ifL3n4UJlTkimzbR/q53k3rzWzCiUVxPc/7aKK9dGsY0DLrSjZ13sVyPQskmEYtwX896utLxeQeH6tBYvJItJmnIK9dsPRmltX5ksTepQ7CqC6ldJparaq8lXygx8p3vkv+XL+Gcfm3SMYldu2j/2f9EfNculFK4ns/F66OcuZLF8TzaE40dGgvWu9hEQyb3bl1Hd/v817tUh8ZA0ZWOk1ihu1CuBnWpXaaU+iOt9UcW25h6XWeRbZDaZWLZGV9EOZYn/+1/pfCVL+FevTL+vgqFSP30m8m8+2eJbNkCBAkAl4ZyDFwewdd+kNprNi77yvY88iUbQxns2NjBhs42QvOcdwm2OPaDsvKxCMlYGFN6LyvComqX1SswNDvACLHc2I5HvmyTH85S/ObXKXz1K3ij2fH3zWSStr3voO2d7ySU6QDA8zRXR/K8dnkYx/NoS0QJNbDnUrZdipZDJBzitk1drO9Iz6uIpe0GNcXQinDYoCMVk83CViGZ+BdiiVSHxMaKFqWhYYonvkrhqRP4xeL4MaHubjI/+x7SP/MzGJU9W1xPcz1bYODKCGXHJRUPkzYbN8xUslyKtkMqHuHebevoaktgziMw2I6HpzWxsElHKk7YNKTXsorVNcgopT4M/JHWuque1xWilVUn8sdKFvbgEMV/+RKFb34d376xeDKyeQuZn/t5Um98IyoU/G/puD6Xh3OcuZLF8z0SsQid0fmvP5kLrTVFy6XsuLQlozywZT0dyblP6GutsV0vqPocDZOOR4k0sKqAaB317sk8CQzU+ZpCtCTHDeYyCmUHb/A6xX/5EvlvfQPtuuPHxHbeRuYXfpHEAw+gKkNfluNxcXCUc9dG0UAqFibcoJ5LsK+LjeP5dKXj3L11DZnk3BduVntnvoZkLExaSr+Im9Q1yGitR4Gn6nlNIVqNU0nxzZUd9PVrFL78RfLf/hba88aPid99Dx3vfS+xu+4e/0LPlx0uXh/j4tAYKGiLRxqWLeb5PrmSg+/7rO9Ms2VN27zWudzIEoNUPEIyFpE0ZFHTgoOMUmorcAR4qPLSSeCg1vrs4pslRGtxPR/LdoMNuVwfPXidwj//E/l//fak4JJ44AE63vs+YjtvA4LJ/NFiidevjzE4WsA0DdoT0YZNjttukHCglGJzdxu3dLcTj87+NeB6Pq7nA0HvRSlFOh6VLDExqwUFGaXUNuAY8ARwtPJyL3BMKfWwBBqxWjiuR65oU7AcDAWMDFP44j+R++Y3JgWX5O7ddLxvH9GenmCIqmhxdaTAxaEcru8RCZl0NKi+mNaaQtnBdj1ikTC339LNukxqTplithPs3xINh+hIxwkZBoahMA0lNcXEnCy0J/N+rXXfTa89BXxMKfX7gKzuFyua7QTzLUUr2C8lUsqR/cIXGHvq5KQ5l+Tu3XS8/2Gi27ZhOR4Xro/x+uAYxXJwXjIeJmQ0Zr7F9jwKJQeNZl0mxcbuNJlEbNZeUrXXojXEoyG64gmZxBcLttAg8/0FvidEy/J9jeW4jBUtbNcjZBhEHIvs//onxr76lUnZYokHeunc9zDhrdsYLZZ59dw1ro0UQEEiGqYz3Zi6YtUsMctxiUbC7NzUydpMataS++PZYRrCpkEmFSMWDslQmFi0hQaZ4RneG1rgNYVYljzfp2y5jBYtfO0TNk1i2mP0S/9M9otfxC8Wxo9N3LuLjocfwdy2nWvZAmd/fIGS7RAJmWSS0YYNMWmtyVeGxLrbEtx1a5AlNluvxfeD4AKQjIZJSnaYqLOFBhmtlNp689xLJRmgc7GNahSpXSbmw/P8YFV+ZcgpGjJRvib39acY+fxx3Gx2/NjYbbfR+YFfxuvZyfmhMS68cA60JhGL0NWgXgtMDC4u3W1Jtm/omFOWmOf7WI6HaRi0JaIkojKBL+anLrXLAJRSTwOZm18GtjF5TYwCOoAPa60/v4A2N5TULhNzVU1BzpcdlIJI2EQBhae/x/A/fBbn8uXxYyNbtpB55Jewdt7B+atjjORLmIZBKh5uaLFK2/UolB0A1nYk2dzdNqeNwqqlXqrBRcq8iMVaVO2yCgU8AmTneM+ZhtKEWJZcz6dkOeTLNp6vMYBYJNiVsfzqKwz97Wcov/LK+PHhNWtJ/OL7GLtjF98fKmCfvko8EmpYhhhUssQsB9sJssTmOt9SHRLzNSSiIVLxuOw4KZbMXILMh7XWMpkvViTH9chVVuUbCsIhg0hlTsK5do3hz/49+e/8x/jxZipF/D0/z5W7+ni1YKOuj5GOR0jHw41pn+dRstxgjYqC7vYkW9a00T7HLDHH9TEMFfRaomFZMCmW3KxBplaAUUq1Vd4ba0SjhGi0G8HFxjQU8ciN3+z9YpGRf/oCo1/+0ng6sgqFiLx9L4MPvpXXHE3EculINW4iv7qlcTQcYkNnms62OG2J6HgAnEk1uIRDBt3tCaJh6bWI5lnoxP9x4GvIehjRYhy3UgXZcjAMRTwSGv8C1r5P7tvfYvgf/mG87L4GuG83l9/0DqxUG1FMOlOhhnxp265HwXLQvqYjFeP2zd10pOJzqoDsa43jePhoTEOCi1g+FhpkjmqtP1frDaXU/VrrHyy8SXOnlMoA+wnmi4aBrNb65FLcW7QWz/MZK1rkyzYhwyAenTy8VX71FQb/5m+wBk4DwW7H7i23cu2t78bdspVULEyyARuDub5PvuTg+T6JaJjt6zvoakuQis++QFNrjeP6wRySoUjGI8QjYcIhQ4KLWDYWGmROK6XerrX+eo33PgD8YOFNmptKgDmmtd5b+bkHeJYgw00IoLI4seyQLVgoxaSeC4CbzTL82b8j9+1vV44HN5Vm6K3vxt3VSzpR/yEx39cUrWBNS8g0uXVtG2szKZKx8Jzu5fk+tuujgHg0TDIWlol8sWwtNMgcATqVUg8wNY15G/DYYhs2B49xo24aWusBpdTuJbivaBFl2yWbL+P6HtFQaNJEufY8xk58jeFjT+IXi0FwUQajP/kW7Lc8RDKdquuXtu9r8mUb19cYSrE2k2BdR4qOZBxzDlsZ+1rjuj6eDtKPO1IxYpFQQ9OkhaiHxZT6/zBT05oV8EeLuOZ8PMpNvRattexlI3Bcj9FCmZLtEgmZxCM3DY298jLX//LT2OfPB8HF9yluv4Piu95LYuN6InUMLpbrUSjbKBSbutpY15EkHY/OObDYjodGo1DSaxEtaaFB5uB0ac1KqccX0Z45qQyNQdCb2kMQ7Hq11kcafW+xfLmeT65oVTLGDBI3zbt4uRzDn/17xr7x9SC4eD5OppPiu99L+N5dJOv0xe36PoXxeZYId9yyhjWZxJwywyAIkq4XzLOkZJ5FtLgFBZlp0pq3EZT7P73YRs1BD8EwXUZrfbxyf5RSR7XWB5bg/mIZ8Spf6rmShQJiN827aK3J/+9/Zegzn8EdG8P1fHzDpPS2d2C8/Z1Ewotf4zJx++JIyGTzmjbWZpKk4pE5BYeJlY9jEZNMKirZYWJFmLWszLwvqNRva60/XdeLTr3HHuCE1lrd9PoIsE1rna1xzqVUKrUhn8/XvOZHP/pRDh061IDWikbxfU3JcsgWLdCaSNjEuOlL2blyhet/+WmKzz+P5weZWO6OO/Df+wi6a82i2+D6PrmSjdbQ3RZnU3fbnOdZIOi1OJ5PJGSSikek8rFYtg4dOsQf/uEfTvd2XmudrvXGgoOMUur/AvYwuSBmBjiutW7oxH81k0xr3XHT6yPAQ1rrUzXOkdplK4jteIzkS8FGXOHQlNXv2vMY/fKXGDr2JG7ZxvN9dCqN+3Pvx9/VC4vsIZQdL9gTxjS4dW076zvTxCNzHxiwXQ/XCzYDa0tEic3jXCGWk3rULpuiEmAAPgL0EQxdDRMMYy3FOpVhphbtFKuA72tyJYuxokUkZE6ZdwGwzp/j2ic/Sen0abzKlsH+T7wJ590/D/Fpi8XOauIOk4lYhHu2rqWrLUlojr2Wm2uIdaTiElzEirfQT3hWa/0pAKVUFtCVsv/fV0rdT4PXyWits0qpU0qpzE1DY5lavRixMjiux3CuhOP5U9a7AGjXZeQL/8jQ5z+P6zigQa9Zi/O+X8Hftn3B9x0fEvM13e1Jbl3bTvsc94bx/KDEiwZMqSEmVqGFBpnxjcm01meUUr8NVOdhlmo/mccJVvsfAVBK7av+Xaws1QWVw/kyYVPVHJayzp3j8p//GeWz5/C1RpkG3psfwnno3bCAiX2tNSXbpex4mIZi69p2NnS2EY/O/L/MjaASpB2bhkE6HiUeDclmYGJVWmiQUZUNyvZprT8O9CmlntZa/xDYC9SqBFBXWuvjSqn9Sqn9VIbOtNYHG31fsbRczyebL1GyXWKR0JSJfe15DH7+8wx//nN4jotSoDZswtr3q+hNm+d/v/FeC3Skouzc1EVnOjHjkFh1GExrTThk0p6MEgmFME0liyXFqrfQFObPVeZlqmMQHwEGlFLtwJJ90Wut+5fqXmJp3SgHU0Ypas69lC9c5MKf/b84p4OsecM0cd+2F+ft7wRzfh/t8Yl8w2BrZSK/1j0ncj0fx/MwlDE+eS+9FSEmW/Cso9b6YxP+niVYGLlNa32mHg0Tq5fjeozky9huUOp+Su9Fa6588Utk/+Hv0baNoRR63XrKD38QvWnLvO5VslxKtkMiFuGuLWtZk0nMOF9yYz1L0GvpSieCNsrukkLUNGOQmaEIZk3TBZj5XkesTl5lxX6+khp8czkYgPy1QV7/8z9Hv/gjlALDMHDf9Dacd/ynec29FMoOZcelLRnlvs3r6UzFpw0U1cDia4iGTTLJGNFISCbvhZiD2XoyB6jP/Eq9riNWqJLlMJwLhsZuXrEPQe/l4rf/ndG/6scoFoLeS6YD6+EP4vfsnNM9tNbkyw6269KZinP31jVkktNvl2y7Hp6vCYcMMqkY0bAEFiHma7Ygs1cp9QRB4cuFlgaQcQQxo0LZZjhXqrmoEsApW5z+1F/hf+spTKVAKbxdvdi/8AGIx2e9/o3g4tHdnqRnXYa2ZHTa423Hw9OaWNikKx0jEpZ5FiEWasYgo7VeqnRksUqNFSxGi+WamWMAI2fOcfFP/gTz0oXg/VgU+xcewbv/wTldP1+ysV2PtR1Jtq7NkE7UDi7Vise+hmQ0TCoekeAiRB2squXGxWKRvr4+APbv38/+/fub3KLVS2vNaKFMrmTXXlipNa9/9SS5//nXmI6NMhT+Lbdif+DX0F3ds16/7AQl9rvTcbZv7Jw2uGitsRwPgFQ8QjIWkSExIeagv7+f/v5+BgcHAaYtpVH3ApnLldQuWz5mCzBOuczpv+jH//dvB70XBe6bH8J5x3tmTU22XY9cySYRDXP75i46U/Gacy5+ZetirTWpeIRULCKFKYVYgIbULhNiobTWZPNlClbtADN85jyX/vi/Y16+GASYRBLrkQ/h337XjNd1PI9c0SYaDnH3lrWszSRrVkJ2PR/bDda2pOJhCS5CNJgEGbFkJgeYyenGvq85941/pfjX/ZhWORge27IN+5d+HZ3pmOaKwTXHKqX279y8hnUdqSnBxdcax/HwgbBp0JWOE4uEZW2LEEtAgoxYEr5/I8DcvJLesR1e+evPoE58eTx7zH3Tz+C86+dmHB6zPY+xgsX6jjS33dJF9KaJer8y36JUsMNkIhqWFflCLDEJMqLhfF8znCtSdrwpAaYwMsrAxz5O6JUfB0Nn0Sj2+38F794Hpr2e1prRoo2hYNe29azNJKbshBlM5ivaE1GSsYj0WoRoknkFGaVUG/ABgiKYPUD7hLezBPvKPAGc1FqP1amNooVpHQQYy/GmVE8efGWAKx8/QmhkEGUodPc6rA/+Fnrt+mmvZ7keuaLFpq42tm/snNR70RPSkNOJCKl4RApUCtFkcwoySqltwGGCBZlPEBTBHNZaj044pp2gzH8vcEQp1QE8rrX+Qb0bLVpHNl+mfFOAGZ9/+cujhBw7WFx5573YD38IYrGa19Faky1ahAyTB3ZsoLttcsak7Xp4nk8iFqEtEZU0ZCGWiVmDjFLqwwSbkj0y03GVgDMKnAE+Vz1XKbWnsh2AWGWqdcgmBpiy7fLa3z2B8eUvYAIocN7+bty3vxOm6XWUHY98yWLzmna2b+gkHLpxXDDv4hIJhehuS8icixDLzGwFMt8PPDmxxzIfWutPKaXalVLv01p/fkEtFC2pZDmM5EvEo+Hx+ZJSscwrf/pnRE59F2UoCEewH/kQ3t331bxGde7FNAx279xIZ3pyCZmg96JpT8ZIxSNz2qlSCLG0Zisr87nF3qASoCTArCIly2FwrEQ8Eh4vFVMcGeX0f3ucyNnXgvmX9gzWh/ajN95S8xqu75PNW6zvSHH75i4ioclzLyU72Aaguy0mvRchlrG6ZZcppR4iGFaTasurWKFkM5yfXOwyd+ES5/+f/4fQ4FWUUvibb8X64Ich3VbzGtU9Xu7a0s3GrvSkHkp1MWV7IkY6Ib0XIZa7us2Oaq2f0lp/XSn19sqf++t17Xqp1i7r6+ujv1821ay3sYIVBJjIjQAz+uJLvP5//wFGJcB499yP9dv/x7QBZqxoo4GfuOMWNnW3TQoiZdvF15p1mRRtyagEGCGaqL+/n76+vqWtXaaUeh/B5D9ABugARir/PKm1Plu3m82T1C5rrGo15YmlYoa+812u/o8/RbkOSqmg/tg7f67mBL/WmpFCmUwixj3b1k1KTfZ9TdlxScUitCdjsuZFiGVkyWqXKaW+ChycLmW5kuIsVqBccWqAufa1Ewx9+lMo7aMMA+c/7cN945trnu/5PiN5i1vWtHHbxq5JZWGCdS+arnScRCyyJP8+Qoj6qeeK/+0zrYlZaIaaWN4KZZtsoTy+m6XWmivHjpN98gmUAsJh7F/6dby7dtU8v1oa5vZbutm85sbwmK81ZdslFgnRkYrLuhchWtRc1skMEazkPwkMAf3TrOY/KqnKq0vZdoMdLSsbjmmtufTXf8PYl78UBIt4HOvXDuDf2lPz/JLlUnZc7t++gTXtN4Z0HdfD9TSdqTiJWFjmXoRoYXPpyZzRWs+6DaHW+mNKqYeUUlubOfcilobjegyNFYmGzSDAeB4XPvFJ8t/8RtCjaW/H+o3fRa/bUPP8saKNaSgevG3jpA3FLMfDUIp1HbKwUoiVYC5BZmCuF9NaP7WItogW4fk+g2NFTENhGgbadXn9T/+Uwne+EwSYrm6s3/o9dEfXlHOrE/ztiRj3Tpjgr659iUdCdKTjUnNMiBViLkFmuOGtEC0jKHhZQmuIhE2043D+4x+j+OypIMCs30j5N3+3Zoqy72tGCiU2dbWzc1MXocoEv+f7WI5HWyJKW0JSk4VYSeYSZFbH/sxiTrL58nhFZd+2ef3IEYo/+EGwyHLLrVi//l8gPjVlPsggK9OzvoOeDR3jgaRaNXlNe5JYRHaeEGKlmcuYxAeUUq8qpT6hlHpvpdz/tJRSb69T28QyUy14GQub+JbF64//txsBZmsP1m/87owB5rZNXZMCTNl2MU2DdR0SYIRYqeYSZJ4Bvk+wj8zngJFZgs7D9W6kaL6S5YynKmvb5sLhP6L4/I+CALN9Z9CDicWnnOd4HtmCxZ2b13DrusyNYpmV9OTutoSkJwuxgs3l18dTWuuPwPiCyr3Anso/DwBaKTUAnCBIc66drypalu14DI0F9chwHC597GMUnn8eFPg7bsP60H6ITF0oWXY8imWbe7euY11Hcvz1ouWQjEbIpGT1vhAr3Vx+hRxfqa+1HtVaH9da/47WegdBuZgPAE8B7wCOEwSgZUlql82f5wWZZOGQgfJcrvz3j5P/4Q8B0NtvnzbAFC2Hsu2ye+fG8QATZJA5pGIROtISYIRoZXWrXaaU+qrW+p1zualSKkNQo6xvHm1dElK7bP58XzOUK+K4HhEFV/7kTxh9+ntoX6O334b1a/shEp1yXqHs4Gvo3bGeVDwyfq2y40oGmRArzGy1y+bSk3lQKVV7V6mbaK2zBHM4YgUYLZSxHY+IaXDtk59grBpgtm3H+s+1A0y+bAPQd9uG8QDjej5lx6MzHac9GZMAI8QqMpcg8wjwB0qp358ts6zi9CLbJJaBfNEmX7aJhk0G//qvGP32v+L5Gr15C9avHYBojQBTsjFUsItlIhoGgt0rXU+zNpMgKQUuhVh1Zp3411qfJJjQRyn1fqXUiWlql1WP/1gd2yeaoGy7jBSCmmQjT/wDI1/7Gq7vo9ZtoPwbv1sziyxfsjEMg94dG4hHg4+V5bgopVibkRIxQqxW81qcUI/tmMXy5no+Q2MlomGT3Fe/wtA//iOu50NXN+Xf+j1IJKecUygH+8VMDDBl2yUcMulKxzElRVmIVWvG//vrtbByKRZoKqWONvoeK53va4bGihgGlL77Ha79zd/geD4q3Yb9m78HbVO3BCpZLp6veWBCgJm4BkYCjBCr22zfAM8qpR6f41zMFEqpdqXU79PgZACl1D6Wcep0K9Bak82XcT0f76Ufc+XP/xzX9VCxGNZv/Bd0V/eUc8q2i+169O7cQDIWzMEULYdkJExnOi4pykKImYfLKhuNPaaU+qPKQsxjWuuvz3ZRpdRDBCv/X9Naf7w+TZ32XplGXn+1yJdsCpZN6OplLvzxx3FsG0Ih7A/+NnrjLVOOt1yPkuPSt2Mj6XhkvIpyKhYsspQMMiEEzHFORmv9kUqQeUQp9TvANoLqzANAFsgA2wkWZw4TrP4/uES7Ye6hkpggFsZyXEYLZULFHJePHMbOFwBwHv4g/o7bpxxvex6FskPvjvW0JYMss5Ltko5HJEVZCDHJnCf+KwHjU5U/KKW2EQSXToLAktVan2lAG6ellOoFTi3lPVca1/MZHCthei7XPv5xrOvX8X2N966fw7tv95TjHc9jrGBz//b1dKSCLLPqKn4JMEKImy249O1SB5Rp9Gitj8uQ2cJorRnJl1C+z/An/4LS6ddwPR+/7424b9075XjP9xkr2ty7de34dslFy5EhMiHEtFo29UcptUdrfXw+5+RyOZRSNf8cOnSoQS1dvsaKFmXbpfCFz5F/+nu4ro+34zac9z4CNwUM39eM5MvccUs36ztTQDBEVi10KQFGiJXt0KFDNb87L1++DJCe7rxZa5fNlVLqk8DXtNafr8sFZ75XBuirLBSt/vys1nr7DOdI7bIJyrbL9dEC3rNPc+3P/geO6+N2rcH5vd+H+OTFlr6vGc6X2Lmpi63rMuPnxyNhOtISYIRYzWarXVbPnaIOAweUUgeAA1rrs3W89s32w/icDEAX0KmUepRgawJJBJiB5/nBFsqvn+P60U8EacvRKO6vH5gSYLTWjBTKbF2X4da1wTqZsu0SDZvSgxFCzKpuQaYyR/ORSq/ij5RSn9Ra/6Be17/pXkcm/qyU6gH23fy6mKo6D+NlRxj8kz/GLVt4Gpxf+U1099opx2eLFusyKbZv6EQpheV4RMImnemErIMRQsyq7nMyWuus1vp3CDY1E8tMsexQLFlkP/kX2EODuJ6P+7PvRd9255Rj8yWbRDTMnVu6MYwgwIQMQxZaCiHmrG5BRin1dKU6wPuUUm1a648ppd5Xr+vPcN99BEN1PUqpwxOG0MRNHNdjJF+m9I/HKb74Aq7r496/G++n3zbl2LLj4WvNfdvWEzINLMclZBh0tccxjZbNFxFCLLFZh8uUUkMEiy5PAkNA/zRVmJ8kWLOyF/gdpVR1YWZDEwEqGWbzyjJbjXxfM5wrYT/7PUb/+YtBgFm3Ae99vzIlk8z1fQplm76dG4lHQ9iOh2kYdLVJgBFCzM9c5mTOaK0fnO2gCSX+n6q+ppR6YKENE/WVK1mULlxg+NP9uJ7Gi0bx/vOHp2yd7PtBDbN7t64jk4phux4opNilEGJB5vKtMbDQi2utv7/Qc0X9WI7LaDZH9i/+DLdQxPN93Ec+hO5aM+m4IJOsxI6NnazvTOG4HlrDmvakBBghxILM5ZtjuOGtEA3j+UG6cuGJz2KdO4/r+ThveTv+XbumHDtSKLOpu52t6zJBWrOGNe0JQhJghBALNJdvj/qs1hRNMVawKHznO+SeOhFkkt1yK/47f27KcaMFi85UnNs2duFrjev5rGmTHS2FEIszlyDzAaXUq0qpTyil3jvb3jJLsUGZmJuS5ZB9/QLZv/40nqfxojG8X/1NMCdPxRUth0jY5O6ta1EKLMejqy1BJCwBRgixOHMJMs8A3wc+AHwOGJkl6Dxc70bWS7FYpK+vj76+Pvr7+5vdnIbyPJ+h0QJjRz+BWyjg+j7Owx9Ed3ROOs5yPWzXZ1dPkKpcdly60nFikXoWgxBCrDT9/f309fUxODgIkJjuuFlrlyml/khr/ZHK39sJUpT3VP70EAynDRDsIXOSoKTMO+vxL1FPq6l2mdbBNsqDx44z9o/HsV0P+w1vxv/FRyYd5/k+2bzFAzvX05mKU7JdMskY6US0SS0XQrSaetQuG9/YvbKnzPi6lJuCzjuA30HmcJquWHYY/fHL5P7p88E8TNc6/Pf84qRjgvIyZW6/pZuudIKi5ZCORyTACCHqai5Bpme6N2oEnQyyS2VTOa7H4GCWsf5P4LkenjLwfvnXIDx5PUy2aLGhM83mNW1Yjkc8EqI9GWtSq4UQK9Vc5mQeVErdN5eLaa2zBHM4ogmqq/oLT34W++oVXM/H3vse9KbNk44rWg7xSJg7NnfjeD6moehIx6WishCi7uYSZB4B/kAp9fuzZZZVnF5km8QC5UoW+eeeI//1k8F6mK3b0W95aNIxjudhOR73bg0qLvsautoSUi5GCNEQsw6XVfZmqW4O9n6l1IlpapdVj//YdO+JxrEcl+xQltG//BS+r/FCYfxHPggTgofva0YLFvf1rCcRDVN2PNZmZLGlEKJx5pWnqrX+XKMaIhZufFX/sX/AqZTvd97z8+jO7vFjtNYM58vs2NjJmvYEJdulMx0nGpZUZSFE48ivsCtArmBRfP558l9/Kggw23agf/LNk44JJvpT3Lo2Q9lxaUtEScYi01xRCCHqQ4JMiyvbLqOjY4z+9advDJM9/KuThslyJZtkLJjot12PRCRCm6QqCyGWgASZFub7wVbKpX/6As6160Ev5l2Th8ls18PzNbu2rsfzfcKmSSYVk0wyIcSSkCDTwnIli/LAAGP/8uXx4pcTh8l8XzNWtLhn61pCpkIpRVebbJ0shFg6qyrIrKTaZY7rMZorMfY3f4XnefhK4e37lUnDZCOFMj3rO+hMx/B8LRuPCSHqpm61y1aKlVS7TGvN9dEi2S/9MyN/97c4rofzM+/Am1DCP1eyScUi3NuzDsf1WNOelKKXQoi6m612mfxa24KKlkPx6jVGP3c82Fysaw3e2981/n7ZCXa0vOvWNTiuRyYZkwAjhGgKCTItxvOCysmFJz+LWyqhtcZ77wcgHA7e930KJZtdPcGK/mQ0QiouqcpCiOaQX29bzFjRwnrhR+T/4z9wPR9/1wP4O24HgmG0bMHitlu6SMYiGEpJJpkQoqmkJ9NCLMclly8y+pm/wfV8VCSK+573jb+fLVqsy6TY2JnG8zSdackkE0I0lwSZFlGtsFw++TXKFy+itcbd8y50ewaAkuUSDYWCBZeeR0c6Rjgk2ycLIZpLgkyLyJds7KFhsp//HK7no9euw33T2wBwfZ+S7bBr21p8rUlFI1IyRgixLEiQaQGO6zFWtMh97kmcYgmlwPn5h8EMBfMweYvbN3cTj4ZRCtpSUjJGCLE8SJBpAaOFMu7ZAXLf/BZaa/y7do1P9meLFhu60mzsTGG7Hl1p2RtGCLF8yLfRMle2XYqWw9Bn/ieu76NCIZyf/UUgmIeJhUPcvqkLq7IeJhKWeRghxPIhQWYZqxbAtJ7+HqWXX8JQ4L7preiuNePzMPduXYvWmmg4JOthhBDLzqoKMq1Wu6xQtnHLFoN/97eggVQa523vrMzDlLn9lm6SsQi+1nSk4rIeRgixZKR22U1arXaZ6/lcHcmT/eL/Yvizf4+hFM4vfgD3DW8iW7Dobktwz9Y1lGyXrnSchGSTCSGaQGqXtaCgp1LCy40x8oUvoBTodetx+36Ssu0SMg1u39yF5XikYhEJMEKIZUuCzDJUtBxKtsu1459Dl4oopbDf+fP4yqBg2dy7dS0KCJkG7clYs5srhBDTkiCzzFQLYHrXr1P6+kkMpfC37cC/4+7x/WHaElE8X8rGCCGWPwkyy8xIvoRSmguf+VuU54EC510/T67skEnG2LI2Q9lx6UhJ2RghxPLXslWYlVIZ4BEgA2wHslrrg81s02IVyzYl2yX/2mvoZ7+HoRTevQ9gbdyMW3K469Y12K5HOi7zMEKI1tCyQQZ4RGs9noeslDqslDqhtd7bzEYtlOf7ZAsWpmFw9TN/S1gBpoG99z2MFSzuvnUdplKEw6bMwwghWkZLBhmlVA9BD2aix4ERpVRGa51d8kYtUr5ko7Xm4nefIfTay2Ao3AffRDbezrq2BF1tMVCKTlkPI4RoIa08J/PYxB8mBJaepW/K4gQFMG1s1yN//MkgiITCFN78EKah2LGxE8+H7rYEptnK/8mEEKtNS35jaa0HtNYdE1+r9G7QWp9qTqsWbrRQxgBOn/wWkUvng8n+n3orY6EYd9+6BkD2hxFCtKSWHC6bxkHgSLMbMV9l26XsuFzP5jG/8r+CXkwszmDfT7NlbYZUPEIkZJKIhpvdVCGEmLeW7MncTCnVC/TMll2Wy+VQStX8c+jQoaVp7ASe7zOUK2G7Phe+9hSR61cBKP3U24ik29i2LoPWkJF5GCFEkx06dKjmd+fly5cB0tOdtyJqlymljmmtH57lmGVXu2x4rETBsnlx4DKxP/lvhLJD6GSKiwcepe+ebUTDITrSMdnlUgixbK342mVKqcPAh5vdjvkq2y4Fy2ZwtIj9H/9GKDsEwOhPvZ1bN68lHg0Ri8gwmRCitbV0kFFKPQo8Xs0sU0plqgkAy5nnBcNkjufz6rkrdPz7UwD46Xb8N/w0m9e0o1FSvl8I0fJaNsgopfYAx29aE/MIMNycFs3daMECX/PahSHaf/g91NgoaBh500PcuW09vtZ0t8UlXVkI0fJaMrus0ls5Ufn7xLeyE6sALEeFsk3BshkeK5EdGWPTv30DAKejkzUPvZ1oJER7IkY03JL/aYQQYpKW/CbTWg8ALTeO5LgeI7ky2odXLg7R/dx3oZADrbF+5t1sX9dJLBIinZCJfiHEyiDjMUvE9zUjuTJKwSsXBwk7FtF/fQo02N3r2fbuvRiGIfMwQogVRYLMEsmVLGzPY3CsyFCuRMfT34ZyCa01sV94H/F4hI5UlJDMwwghVpCWHC5rNWXbZbRQRmt45cIQGd8i9G/fAg3u5q1sesubZBtlIcSKJL82N1g1XTkSMnnp9UFCpkHsWyfAtvG1puuXfolYJERbMtrspgohRN1JT6aBtNYM50soNJeHC2QLZbqdAqHv/Ttaa4w77iRz/310puOYhsR7IcTKs6q+2YrFIn19ffT19dHf3/hM53zJxnJcHNfnlYtDZJJRwk99BSrbKq/95V8hFYtIurIQouX09/fT19fH4OAgQGK641ZE7bK5WOraZWXb5fpogUjI5NTpy9i2R2r4KrE/P4Lva+J9D7Lh//x91meSsuhSCNGyVnztsuXIcT2GxkpEwyYXB8fIFSyS0RCRL/8jWmvMUIjuX/plMsmoBBghxIom33B15vk+Q2MlTAOKZYdXLw3TnoxivvQjjNOvgob0nj0kb90sxS+FECueBJk60lozkivha41SBi++fp14JISpfcJf/gK+1oTSKdreu49MMiaLLoUQK54EmToaLZQp2S7RsMnpS8PkSzaJaJjQd/8NBq9hKEX7e99PZm2nbKUshFgVJMjUSaFkkyvZxCMhrgznOX89S0cyBvkcoZNfBg3RDRtI730HqbgsuhRCrA4SZOqgbLsM50vEIiEKZYcXz1+nPRlFKUXki8ehWMQ0DTK/8qt0tCVlTYwQYtWQb7tFCjLJikTDJr6vef7sNSIhk7BpYr7wHOZz3wcFyZ94A+kHHyQRk8l+IcTqIUFmERzX49pokZBpYCjFKxeGKNkOyVgYSkXC//QkvtaEU2naPvifaZfJfiHEKiNBZoGqAcZUEDINLg7luDg0RiYR1CCLfOkfITeKYSg6P/Qh2tatJRaRlf1CiNVFgswCOK7H9UqACYdMxgoWL78+SCYV9FTMHzyN+ex30RoS991H/E1vlgKYQohVaVUFmXrULvM8n8GxIqoSYGzX4/mz14JUZcPAOHuayPG/R2tNKJmk4zd+i/ZkTPaJEUKsKFK77Cb1qF3m+T6Do0V8XxMJm3i+5rmBK4wWrCCbbOg6sU/8MRQKaMNgw0f+K6ld97I2k5S5GCHEiiS1y+rE9zXDuRJeJcAADFwaZihXoj0ZhbFRov9fPxQL+FqT/OCvE7377vEhNCGEWI1kJnqOCmUby/GIVybvLw3lOHstS1c6jvHay0Sf+J9QyKG1xn7zQ6x5x14SkbCU8RdCrGryDThHvtYYlQ7JaKHMi+evk4kahJ/6CuGvfwXQoCF/34Ps+I1fQ6NJJ2RlvxBidZMgMw/a9ylcH+TF//gBa370faI/fg6scvBmKMzou95L8i1vJRmPEgmbUp9MCLHqSZCZg5f/6x9QfP0CTnYUx3bIoCfNs+g16yj90q9TSHeya0Mnru/TnZg22UIIIVYNCTJz4AwO4Q6P4HoeWmuUocAw8W6/C++BB/HuvIds2WXr2nZMUxE2pRcjhBAgQWZOMj/5BpyLV7nqKJJr10BbBq9nBySSANiuR8gw2bI2g+P5dKVjTW6xEEIsDxJk5uCW3/wNrr1yieLlYWLp+JT3cyWbe25di681sbA5nuIshBCrnayTWaSi5dCWjLI2k8T1fNoS0osRQogqCTKLoLWmZLvcvqkL1/NJRMPSixFCiAlWVZCpR+2yiXIlh/WdadqTMTzfpy0hRTCFEKuD1C67yWJrl5165RKvXR6mszIn4/k+Y0WLN965BcOAWDhMZ9vU+RohhFjJpHZZg4yVbLauyxCLBDtiyup+IYSYSoLMAtieR8gw2LI2g+16JGIRWRcjhBA1SJBZgLGCzc6NXZiGwvMhHZdejBBC1CJBZp5yJZs1bXHWdSSxXY9ULCy9GCGEmEZLL8ZUSu0Hhis/dmqtF58yNgPP91Ge4vbNa1AKPB9S0osRQohptWxPphJgBrTWx7XWx4HhymsNky873Lapi3g0hO14pKUXI4QQM2rZIAMc0FqfrP5QCTQHGnnDDR0pNnal0VrjaUjJuhghhJhRSw6XKaUyQE+Nt3qVUhmtdbbe91zfmSKTjqMUlGyXdCxCyGzlGC2EEI3Xqt+SfdyYi5koS+3gs2jJeIRwyKBku7QlorQlpUaZEELMpiV7MkCGIKDcbBjobMQNFYqQadCZjhMNt+pjE0KIpdWqPZkFyeVyKKVq/jl06NCM5ybjYdZmkhJghBCr0qFDh2p+d16+fBkgPd15LVm7TCm1Bziqtd5+0+sjwENa61M1zllU7TIhhBBTrdTaZc9Qe1gsAwwsbVOEEEJMpyWDTCV7rNbE/0AjMsuEEEIsTEsGmYqjExdfVv5+uIntEUIIcZOWDTJa6yMQzM8opfZVXmtIWZnZkgJEbfLcFk6e3cLIc1u4Rj27lpz4X4jFTPwrpVgtz6me5LktnDy7hZHntnALfXYrdeJfCCFEC5AgI4QQomEkyCyB/v6FTxUt5tzlcH4z793Kz72Z927l575Yrfzv3uxnNy2t9ar4A1zasGGDXojgMS3c7t27m3Jus89v5nNb7PnNfu6LeXbNbnszn7t85hZuoc9uw4YNGrikp/nuXU0T/45hGKF169bN+9zLly+zYcOGBd97cHCQ7u7uJT+32ec387kt9vxmP/fFPLtmt72Zz10+c0v/7K5evYrv+67WOlzr/dUUZIpAGLi+gNPTQG4Rt08AxSac2+zzm/ncFnt+s5/7Yp5ds9vezOcun7mFW+izWwM4WutErTdXTZARQgix9GTiXwghRMNIkBFCCNEwEmSEEEI0jAQZIYQQDSPbPM6gUtm5uqVAp25QAc6VpFKsdC83KmLvA05prU82r1XLk1IqA+wHtmutD9R4Xz5/Ncz03OTzN7PKs3uEYO+t7UBWa33wpmPq+rmTnsw0Kg96QGt9XGt9HBieuLWAmFYn0AecBk4QfIjlf/CbKKV6gT1Adpr35fNXw2zPDfn8zeYRrXW/1vpINUArpU5U32zE505SmKehlHpWa717ttfEZEqpfZUPp5iD6m/eNX4jl8/fDGZ4bvL5m4ZSqgfYpyvbpFReywAjQIfWOtuIz530ZGqoPPieGm/1Vt4TomHk8yca6LGJP+gbOwn3NOpzJ3MytfVRe3vnLMF/hFNL2poWo5TaM+HH3om/OYk5kc/fIsjnrzat9QDQMfG1Su8GrfWpynOr++dOgkxtGWqP+Q4TjPmK6Q0QjIMPQPBbuVLq8M2Ti2JGGeTzt1Dy+Zufg0A1CGdowOdOhstEXWmtT1X/B6/8fBx4tIlNEquIfP7mrpJE0dPoACxBprYsQVS/WSe1u5NiFtVuuZiTLPL5qyv5/NX0mNZ674SfszTgcydBprZnqN09zBB0x0UNSqkepdTpGm9ll7otLU4+fwsgn7+5U0odBj5808sN+dxJkKmhknFRK3IPTMjGELUdnvhDJSslM3EIQ8xMPn+LIp+/WSilHgUer36WKvNWPY363EmQmd7RiYuQKn8/PMPxq940/yM/BkxZzS4mqfXbo3z+Zjfpucnnb3aVDLLjNwWNR7gRXOr+uZPFmDOorn4l6C5KWY85mFC2AirZKvLcpqoujAM+QJAe2g+cmLg6XT5/U8323OTzN73Ks6s5nKi17phwXF0/dxJkhBBCNIwMlwkhhGgYCTJCCCEaRoKMEEKIhpEgI4QQomEkyAghhGgYCTJCCCEaRoKMEEKIhpEgI4QQomEkyAghhGgYCTJCCCEaRoKMEEKIhpEgI4QQomEkyAjRREqpEaXUsUoJ9ma14XClDcea1QaxckkVZiFmUQkAxwh2Dhwg2GmxF+gDnqwc1klQer5nYtn0OVz7mNb64Tq1sxf4FIDWevcCzq9bW4SoCjW7AUK0gAPAQ1rrU9UXKr/1D2utJ22IpZR6dqkbV6W1PqWU+jBBQBRiWZDhMiFm9/TEAFOxBzhZ49hary2lbJPvL8QkEmSEmIFSah/B7osTX8sQ7Bp4osYpQ41vlRCtQ4KMEDM7edN+6BD0YqB2r0W2+hViApmTEWIGNQIMwF5goNZ70xw/L5XeEwTJBBmt9ZEJ7/UChys/HiRINjgAPDzx3pXjqtfoBY5rrQcW2zYh5kt6MkLM33TzMYumlHoUOKW1Pq617p/wGhBM7hMEl06gU2t9nKnDdj1AVmt9Smt9shKkjimlehrRZiFmIkFGiHmozMf0UHs+ph62A/sm/HySoOc0URbo1VqfBNBaH7mpBzVQo9fyBDd6QEIsGRkuE2J+ZpqPWbSJKdGVnkcfQa/lZvMd+joFPLaIpgmxINKTEWJ+pp2PqQelVEYpdVQptb/yksyjiJYmQUaI+ZlxPkYp9ahS6rRS6sSE1+ZTMuYMcFRr3V8Z8hqecJ3FzKn00vw1PGIVkiAjxBzNNh9TCSYDWuvtwMEJWWJzvX4PQTbZxIWfPQRrciAIFHPRWSMgfYAgYaDaW9o/9TQh6k9qlwkxi0p213aC+ZFqj2AAOFadfK8c1zNxwr3ag5l4TI1rT6oXVrlXF0EgG67c5zHgNEHtNCo/7wOOAE/cVO4mQ9DbylbOn5LCXK3FdnONNaldJhpBgowQDbKQINNMy6ktYuWQ4TIhGifDhDkVIVYjCTJCNE5njcKaQqwqEmSEaBzpxYhVT4KMEI1TaxHlFMtlZ0zg6Wa1QaxcMvEvRANUsrw6pSilWO0kyAghhGgYGS4TQgjRMBJkhBBCNIwEGSGEEA0jQUYIIUTDSJARQgjRMBJkhBBCNMz/D3A41shKEZDPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 443.077x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.fill_between(T_gamma.detach().numpy(), np.percentile(gpr_samples,[16], axis=0)[0], np.percentile(gpr_samples, [84], axis=0)[0], alpha=0.3, color=cols_default[1], label=\"GP post.\")\n",
    "plt.fill_between(T_gamma.detach().numpy(), np.percentile(gpr_samples,[2.5], axis=0)[0], np.percentile(gpr_samples, [97.5], axis=0)[0], alpha=0.1, color=cols_default[1])\n",
    "plt.plot(T_gamma, T_nu(T_gamma), label=\"True\")\n",
    "plt.xlabel(r\"$T_\\gamma$\\,[arb.]\")\n",
    "plt.ylabel(r\"$T_\\nu(T_\\gamma)$\\,[arb.]\")\n",
    "plt.legend(loc='upper left')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f55affa-c4f9-4fbd-949c-e83848ae1029",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ff531ba-f0ca-4386-8423-60a86e2debc8",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}