{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Metrics in QPyTorch\n", "\n", "In this notebook, we will see how to evaluate QPyTorch models with probabilistic metrics. \n", "\n", "**Note:** It is encouraged to check the Simple QP Regression notebook first if not done already. We'll reuse most of the code from there.\n", "\n", "We'll be modeling the function\n", "\n", "$$\n", "\\begin{align}\n", " y &= \\sin(2\\pi x) + \\epsilon \\\\\n", " \\epsilon &\\sim \\mathcal{N}(0, 0.04) \n", "\\end{align}\n", "$$\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import math\n", "import torch\n", "import qpytorch\n", "from matplotlib import pyplot as plt\n", "\n", "%matplotlib inline\n", "%load_ext autoreload\n", "%autoreload 2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the next cell, we set up the train and test data." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Training data is 100 points in [0,1] inclusive regularly spaced\n", "train_x = torch.linspace(0, 1, 100)\n", "# True function is sin(2*pi*x) with Gaussian noise\n", "train_y = torch.sin(train_x * (2 * math.pi)) + torch.randn(train_x.size()) * math.sqrt(0.04)\n", "\n", "test_x = torch.linspace(0, 1, 51)\n", "test_y = torch.sin(test_x * (2 * math.pi)) + torch.randn(test_x.size()) * math.sqrt(0.04)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the next cell, we define a simple QP regression model." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# We will use the simplest form of QP model, exact inference\n", "POWER = 1.0\n", "class ExactQEPModel(qpytorch.models.ExactQEP):\n", " def __init__(self, train_x, train_y, likelihood):\n", " super(ExactQEPModel, self).__init__(train_x, train_y, likelihood)\n", " self.power = torch.tensor(POWER)\n", " self.mean_module = qpytorch.means.ConstantMean()\n", " self.covar_module = qpytorch.kernels.ScaleKernel(qpytorch.kernels.RBFKernel())\n", " \n", " def forward(self, x):\n", " mean_x = self.mean_module(x)\n", " covar_x = self.covar_module(x)\n", " return qpytorch.distributions.MultivariateQExponential(mean_x, covar_x, power=self.power)\n", "\n", "# initialize likelihood and model\n", "likelihood = qpytorch.likelihoods.QExponentialLikelihood(power=torch.tensor(POWER))\n", "model = ExactQEPModel(train_x, train_y, likelihood)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our model is ready for hyperparameter learning, but, first let us check how it performs on the test data." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhsAAAEYCAYAAAD1WzSOAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAWVBJREFUeJzt3Xl4VNX5wPHvLJk1y2QlCwkQdhBEQBRccEGBWhRt3cpPpXVrhbpgrVKraKtFxVo3tNaNqljcQK07UkBF2RdZQiAhG4Hs+0xmycz9/REyZpkkMyGT9f08zzzJ3Ll37rkh5L5zznveo1IURUEIIYQQIkjU3d0AIYQQQvRtEmwIIYQQIqgk2BBCCCFEUEmwIYQQQoigkmBDCCGEEEElwYYQQgghgkqCDSGEEEIElQQbQgghhAgqCTaEEEIIEVQSbAghhBAiqIIabLz44ouMHz+e8PBwwsPDmTp1Kp9//nkwTymEEEKIHkYVzLVR/vvf/6LRaBg+fDiKovDvf/+bZcuWsWvXLsaOHRus0wohhBCiBwlqsOFLVFQUy5Yt48Ybb+zK0wohhBCim2i76kRut5v33nsPq9XK1KlTfe7jcDhwOBze5x6Ph7KyMqKjo1GpVF3VVCGECBpFUaiuriYxMRG1WtLmRP8Q9GBj7969TJ06FbvdTmhoKGvWrGHMmDE+9126dCkPP/xwsJskhBDdLi8vj4EDB3Z3M4ToEkEfRnE6neTm5lJZWcn777/PK6+8wsaNG30GHM17NiorK0lJSSEvL4/w8PBgNlMIIbpEVVUVycnJVFRUEBER0d3NEaJLdHnOxowZMxg6dCgvvfRSu/tWVVURERFBZWWlBBtCiD5B/q6J/qjLBww9Hk+T3gshhBBC9G1BzdlYvHgxs2fPJiUlherqat5++202bNjAl19+GczTCiGEEKIHCWqwUVRUxPXXX8/x48eJiIhg/PjxfPnll1x00UXBPK0QQgghepCgBhuvvvpqMN9eCCFEGzweD06ns7ubIfqgkJAQNBqN3/t3WZ0NIYQQXcfpdJKVlYXH4+nupog+ymKxEB8f71cdLAk2hBCij1EUhePHj6PRaEhOTpbiYaJTKYqCzWajqKgIgISEhHaPkWBDCCH6mLq6Omw2G4mJiZhMpu5ujuiDjEYjUJ+bGRcX1+6QioS7QgjRx7jdbgB0Ol03t0T0ZQ2BrMvlandfCTaEEKKPkjWlRDAF8vslwYYQQgghgkqCDSGEEL3O4MGDefrpp7u7GZ2mr11PcxJsCCGE6DHy8vL4zW9+Q2JiIjqdjkGDBnHHHXdQWlra3U3rVg899BAqlQqVSoVWqyUmJoZzzz2Xp59+OuAlQDZs2IBKpaKioiI4jfVBgg0hhBCt2r59OxdccAHbt28P+rmOHDnC5MmTOXz4MP/5z3/IyMjgn//8J+vWrWPq1KmUlZUFvQ2tcbvd3V6zZOzYsRw/fpzc3FzWr1/PlVdeydKlS5k2bRrV1dXd2rb2SLAhhBCiVW+88Qbr16/nzTffDPq5FixYgE6n46uvvmL69OmkpKQwe/Zsvv76a/Lz87n//vub7F9dXc21116L2WwmKSmJ5cuXe19TFIWHHnqIlJQU9Ho9iYmJ3H777d7XHQ4Hf/jDH0hKSsJsNnPGGWewYcMG7+srVqzAYrHw8ccfM2bMGPR6Pa+88goGg6FFj8Add9zBBRdc4H3+3Xffcc4552A0GklOTub222/HarV6Xy8qKmLOnDkYjUaGDBnCypUr/fr5aLVa4uPjSUxMZNy4cfz+979n48aN7Nu3j8cff9y735tvvsnkyZMJCwsjPj6eX/3qV96aGNnZ2Zx//vkAREZGolKpmD9/PgBffPEFZ599NhaLhejoaH7+85+TmZnpV9vaI8GGEEKIJnJyctixYwc7d+7knXfeAWDVqlXs3LmTHTt2kJOT0+nnLCsr48svv+S2227z1nBoEB8fz7x583jnnXdQFMW7fdmyZZx66qns2rWL++67jzvuuIO1a9cC8MEHH/CPf/yDl156icOHD/Phhx8ybtw477ELFy7khx9+YNWqVfz4449ceeWVzJo1i8OHD3v3sdlsPP7447zyyivs37+fefPmYbFY+OCDD7z7uN1u3nnnHebNmwdAZmYms2bN4he/+AU//vgj77zzDt999x0LFy70HjN//nzy8vJYv34977//Pi+88II3GAjUqFGjmD17NqtXr/Zuc7lc/PWvf2XPnj18+OGHZGdnewOK5ORkb/vT09M5fvw4zzzzDABWq5VFixaxfft21q1bh1qt5vLLL++cHh2lB6usrFQApbKysrubIoQQnaIr/q7V1tYqBw4cUGprazt0POB9qFSqJl8bHp1t8+bNCqCsWbPG5+tPPfWUAiiFhYWKoijKoEGDlFmzZjXZ5+qrr1Zmz56tKIqi/P3vf1dGjBihOJ3OFu+Vk5OjaDQaJT8/v8n2Cy+8UFm8eLGiKIry+uuvK4Cye/fuJvvccccdygUXXOB9/uWXXyp6vV4pLy9XFEVRbrzxRuWWW25pcsy3336rqNVqpba2VklPT1cAZevWrd7X09LSFED5xz/+0cpPR1GWLFminHrqqT5fu/feexWj0djqsdu2bVMApbq6WlEURVm/fr0CeNvcmuLiYgVQ9u7d6/P1QH7PpGdDCCFEE2+99RZabX2BaeVET0LDV61Wy1tvvRW0cyuNei7aM3Xq1BbP09LSALjyyiupra0lNTWVm2++mTVr1lBXVwfA3r17cbvdjBgxgtDQUO9j48aNTYYNdDod48ePb3KOefPmsWHDBo4dOwbAypUrueSSS7BYLADs2bOHFStWNHnfmTNn4vF4yMrKIi0tDa1Wy6RJk7zvOWrUKO/xHaEoSpOaFzt27GDOnDmkpKQQFhbG9OnTAcjNzW3zfQ4fPsy1115Lamoq4eHhDB482K/j/CHlyoUQQjQxb948Ro8e3eSG2GDLli1MnDix0885bNgwVCoVaWlpXH755S1eT0tLIzIyktjYWL/eLzk5mfT0dL7++mvWrl3LbbfdxrJly9i4cSM1NTVoNBp27NjRosx2aGio93uj0diicNXpp5/O0KFDWbVqFb/73e9Ys2YNK1as8L5eU1PDrbfe2iQ/pEFKSgqHDh3yq/2BSEtLY8iQIUD9UMjMmTOZOXMmK1euJDY2ltzcXGbOnNnuCsBz5sxh0KBBvPzyyyQmJuLxeDjllFM6ZeXgPhtsbM0qo9zmxBiiwaTTYDjx1aTTYtRpMOs0aDXSsSOEEG1Rq9V4PB7v12CJjo7moosu4oUXXuCuu+5qkrdRUFDAypUruf7665vc/Ddv3tzkPTZv3szo0aO9z41GI3PmzGHOnDksWLCAUaNGsXfvXk477TTcbjdFRUWcc845Abd13rx5rFy5koEDB6JWq7nkkku8r02cOJEDBw4wbNgwn8eOGjWKuro6duzYwemnnw7U5050dBrqwYMH+eKLL1i8eLH3eWlpKY899hjJyckALWYSNZSxbyhrD1BaWkp6ejovv/yy92fy3XffdahNvvTZYKOgyk5mUU2b++hD1Jh1Wsx6LWadhlBD/fdh+vqvoQYtoTotarWU/BVC9C9xcXHEx8eTnJzMjTfeyKuvvkpeXh5xcXFBO+fzzz/PtGnTmDlzJo888ghDhgxh//793HPPPSQlJfHoo4822X/Tpk088cQTzJ07l7Vr1/Lee+/x6aefAvWzSdxuN2eccQYmk4m33noLo9HIoEGDiI6OZt68eVx//fX8/e9/57TTTqO4uJh169Yxfvz4JsGDL/PmzeOhhx7i0Ucf5Ze//CV6vd772r333suZZ57JwoULuemmmzCbzRw4cIC1a9fy/PPPM3LkSGbNmsWtt97Kiy++iFar5c4772yRFOtLXV0dBQUFeDweSktL2bBhA4888ggTJkzgnnvuAep7T3Q6Hc899xy//e1v2bdvH3/961+bvM+gQYNQqVR88skn/OxnP8NoNBIZGUl0dDT/+te/SEhIIDc3l/vuu8+vfzd/9Nlgwx8OlweHy0mZtfUuIrVKhVmvIdwQQphBS7ix/muEMYQIYwhhhhA0EowIIfqYgQMHkp2djU6nQ6VSccstt+B0OpvcWDvb8OHD2b59O0uWLOGqq66irKyM+Ph45s6dy5IlS4iKimqy/91338327dt5+OGHCQ8P56mnnmLmzJkAWCwWHnvsMRYtWoTb7WbcuHH897//JTo6GoDXX3+dRx55hLvvvpv8/HxiYmI488wz+fnPf95uO4cNG8aUKVPYunVri6qf48ePZ+PGjdx///2cc845KIrC0KFDufrqq737vP7669x0001Mnz6dAQMG8Mgjj/DAAw+0e979+/eTkJCARqMhIiKCMWPGsHjxYn73u995/11iY2NZsWIFf/rTn3j22WeZOHEiTz75JJdeeqn3fZKSknj44Ye57777+PWvf83111/PihUrWLVqFbfffjunnHIKI0eO5Nlnn+W8885rt13+UCmBZON0saqqKiIiIqisrCQ8PDygYz/ec6zdno3OoFJBqL4++LCYdFhMIViMIUSYQrAYdei0MlQjhPjJyfxd85fdbicrK4shQ4ZgMBiCcg4hAvk969c9G51BUaDaXke1vY6j5bUtXg8zaLGYdESZQ4g06Ygy64g06wjTa2VFRiGEEP2CBBtB1hCI5DWrsqvTqoky64g264gO1RFt1hMdqiPMENI9DRVCCCGCRIKNbuKs81BQaaeg0t5kuyFEQ3SojthQPbFhemJC64OQEJk5I4QQopeSYKOHsbvc5JfXkt9oSEatUhFpDiEuTE9smOHEVz2GEE0b7ySEEEL0DBJs9AIeRaG0xklpjZO04/Ur+6lUEGEMIT7cQFy4gQHheuLCDJKQKoQQoseRYKOXUhSosLmosLk4WPBTABJt1hEfYSQ+3EB8hIGYUJ0kogohhOhWEmz0IYoCJTVOSmqc7MuvBOoTUQeEG0iIMJBoMZIQYZDhFyGEEF1Kgo0+zlnnIa/MRl6ZDajv/Ygy60iIMJJoMZBkMWIx6bq5lUIIIfoyCTb6GUXBm//R0PsRqteSFGkk0WIkyWKUoRchhBCdSoINQY2jjvSCatJP5H4YQjQkRRoZeOIRG6qX4EMIIUSHBXXqwtKlSzn99NMJCwsjLi6OuXPnkp6eHsxTik5gd7nJLKphY3oxKzfn8tI3R/jvnmPszqugtMbR3c0TQvRR8+fPR6VS8dvf/rbFawsWLEClUjF//vyub5g4aUENNjZu3MiCBQvYvHkza9euxeVycfHFF2O1WoN5WtHJap1uMopqWH+wiDd+yOHlb47wxb7j7MuvpMru6u7mCSH6kOTkZFatWkVt7U+1hux2O2+//TYpKSnd2DJxMoIabHzxxRfMnz+fsWPHcuqpp7JixQpyc3PZsWNHME8rgqzGUUfa8WrWHijk1W+zWLEpi/UHi8goqsFR5+7u5gkherGJEyeSnJzM6tWrvdtWr15NSkoKp512mnebx+Nh6dKlDBkyBKPRyKmnnsr777/vfd3tdnPjjTd6Xx85ciTPPPNMk3PNnz+fuXPn8uSTT5KQkEB0dDQLFizA5ZIPUZ2tS3M2KivrExKbLxPcwOFw4HD81E1fVVXVJe0SJ6fc5qLcVsHuvArUKhXxEXoGRZsZFG0iPtwg+R5CdDNFAZute85tMtXPggvEb37zG15//XXmzZsHwGuvvcavf/1rNmzY4N1n6dKlvPXWW/zzn/9k+PDhfPPNN/zf//0fsbGxTJ8+HY/Hw8CBA3nvvfeIjo7m+++/55ZbbiEhIYGrrrrK+z7r168nISGB9evXk5GRwdVXX82ECRO4+eabO+PyxQldtsS8x+Ph0ksvpaKigu+++87nPg899BAPP/xwi+09eYl50TZDiIZB0SZSokwMjjETqpecZNG/dccS81YrhIYG5VTtqqkBs9m/fefPn09FRQUvv/wyycnJ3hy/UaNGkZeXx0033YTFYuGll14iKiqKr7/+mqlTp3qPv+mmm7DZbLz99ts+33/hwoUUFBR4e0Dmz5/Phg0byMzMRKOprz901VVXoVarWbVq1Ulcdf/QI5eYX7BgAfv27Ws10ABYvHgxixYt8j6vqqoiOTm5K5ongsTucjeZ6RIbpmdwtJnBMSYSI4yo1dLrIYRoKjY2lksuuYQVK1agKAqXXHIJMTEx3tczMjKw2WxcdNFFTY5zOp1NhlqWL1/Oa6+9Rm5uLrW1tTidTiZMmNDkmLFjx3oDDYCEhAT27t0bnAvrx7ok2Fi4cCGffPIJ33zzDQMHDmx1P71ej16v74omiW5SXO2guNrBtuwy9CFqBkXVBx5DYsyYdNLrIUQwmEz1PQzdde6O+M1vfsPChQuB+qChsZoTF/Ppp5+SlJTU5LWGe8iqVav4wx/+wN///nemTp1KWFgYy5YtY8uWLU32DwkJafJcpVLh8Xg61mjRqqD+dVcUhd///vesWbOGDRs2MGTIkGCeTvQyDpeHQ4XVHCqsRqWCAeEGhsSYGRJjJi5MansI0VlUKv+HMnqKWbNm4XQ6UalUzJw5s8lrY8aMQa/Xk5uby/Tp030ev2nTJqZNm8Ztt93m3ZaZmRnUNovWBTXYWLBgAW+//TYfffQRYWFhFBQUABAREYHRaAzmqUUvoyhQUGmnoNLOD5mlhOq1DI4xkxprJiXKRIhGVrMVoj/RaDSkpaV5v28sLCyMP/zhD9x11114PB7OPvtsKisr2bRpE+Hh4dxwww0MHz6cN954gy+//JIhQ4bw5ptvsm3bNvnQ202CGmy8+OKLAJx33nlNtr/++utSmEW0qcZRx778SvblVxKiUZEcVT/UkhobKkmmQvQTbSXQ/vWvfyU2NpalS5dy5MgRLBYLEydO5E9/+hMAt956K7t27eLqq69GpVJx7bXXctttt/H55593VfNFI102G6UjTiZrW2aj9E0qFcSFGUiNre/1iAtrOwNaiJ6mO2ajCBEMPXI2ihCdQVGgsMpOYVX9cEu4MYTUWDNDY0IZGCmzW4QQoieSYEP0alW1LnbnVrA7twJDiIYhMSaGxoYyKNqMTit5HkII0RNIsCH6DLvLTdrxatKOV6NVq0iJNpEaE8rQOJlWK4QQ3Un+Aos+qc6jcKTYypFiK+sOQmKEkaFxZobFhhFhCmn/DYQQQnQaCTZEn6cokF9RS35FLd8cKiEmTM/QWDPD4kIlwVQIIbqABBui3ympdlBS7WDLkTIijCEMjQtlWFwoiRGyaJwQQgSDBBuiX6usdbEzp5ydOeWY9RpSY+oDj+QoExqZ2SKEEJ1Cgg0hTrA63OzNr2RvfiX6EDWpMfVDLYOizVLBVAghToIEG0L44HB5vDNbQjQqUqLNDIsNJTXWjCFE0/4bCCGE8JKPa0K0w+VWyCyq4cv9BfzrmyOs3nmUH49WYHXUdXfThBAdpCgKt9xyC1FRUahUKnbv3s15553HnXfe2eZxgwcP5umnn+6SNvYl0rMhRADcHoWcUhs5pTb+d7CIhAgDw+JCZUqt6BX+sfZQl57vrotGdOi4goICHn30UT799FPy8/OJi4tjwoQJ3HnnnVx44YWd0rYvvviCFStWsGHDBlJTU4mJiWH16tUtlpwXnUOCDSE6SFHgWIWdYxV2vjlUQmyYnmFxoQyNDSU2TN/dzROiV8rOzuass87CYrGwbNkyxo0bh8vl4ssvv2TBggUcPHiwU86TmZlJQkIC06ZN826LiorqlPcWLckwihCdpLjawQ+Zpby1OYfXN2Xx7eFijlXU0oPXOhSix7nttttQqVRs3bqVX/ziF4wYMYKxY8eyaNEiNm/eDEBubi6XXXYZoaGhhIeHc9VVV1FYWOh9j4ceeogJEybw5ptvMnjwYCIiIrjmmmuorq4GYP78+fz+978nNzcXlUrF4MGDAVoMoxQVFTFnzhyMRiNDhgxh5cqVLdpbUVHBTTfdRGxsLOHh4VxwwQXs2bPH77YAeDwennjiCYYNG4ZeryclJYVHH33U+3peXh5XXXUVFouFqKgoLrvsMrKzszvjx91lJNgQIggqbC62Z5fzzrY8Xvk2i3VphWSXWHF7JPAQojVlZWV88cUXLFiwALPZ3OJ1i8WCx+Phsssuo6ysjI0bN7J27VqOHDnC1Vdf3WTfzMxMPvzwQz755BM++eQTNm7cyGOPPQbAM888w1/+8hcGDhzI8ePH2bZtm8/2zJ8/n7y8PNavX8/777/PCy+8QFFRUZN9rrzySoqKivj888/ZsWMHEydO5MILL6SsrMyvtgAsXryYxx57jAceeIADBw7w9ttvM2DAAABcLhczZ84kLCyMb7/9lk2bNhEaGsqsWbNwOp0d+0F3AxlGESLIahx1/Hi0kh+P1k+pHRJtZmhcKINlsTghmsjIyEBRFEaNGtXqPuvWrWPv3r1kZWWRnJwMwBtvvMHYsWPZtm0bp59+OlDfW7BixQrCwsIAuO6661i3bh2PPvooERERhIWFodFoiI+P93meQ4cO8fnnn7N161bve7766quMHj3au893333H1q1bKSoqQq+vHzp98skn+fDDD3n//fe55ZZb2m1LdXU1zzzzDM8//zw33HADAEOHDuXss88G4J133sHj8fDKK694iw6+/vrrWCwWNmzYwMUXX9yBn3TXk2BDiC7kcHk4WFDNwYL6xeKSo+pXqU2NNWPWy39H0b/5M+SYlpZGcnKyN9AAGDNmDBaLhbS0NG9gMHjwYO/NHSAhIaFFr0R759FqtUyaNMm7bdSoUVgsFu/zPXv2UFNTQ3R0dJNja2tryczM9D5vqy1paWk4HI5WE1/37NlDRkZGk+MB7HZ7k3P0dPLXTYhuUudRyCqxklViRXUQEiIMpMbWJ5hGmXXd3Twhutzw4cNRqVSdkgTafFaJSqXC4/Gc9Ps2VlNTQ0JCAhs2bGjxWuOgpK22GI3Gds8xadIkn/kisbGxgTe6m0gfbi+Vd2gvL9xzPXmH9nZ3U0QnaJjZ8t3hEv79fTYrNmXxzaFijpbb8Eieh+gnoqKimDlzJsuXL8dqtbZ4vaKigtGjR5OXl0deXp53+4EDB6ioqGDMmDGd1pZRo0ZRV1fHjh07vNvS09OpqKjwPp84cSIFBQVotVqGDRvW5BETE+PXeYYPH47RaGTdunU+X584cSKHDx8mLi6uxTkiIiJO6hq7kgQbvdS2tR+RsWcL27/+SAKPPqjc5mJHTjnvbT/Kv749whf7CjhcWI2jzt3dTRMiqJYvX47b7WbKlCl88MEHHD58mLS0NJ599lmmTp3KjBkzGDduHPPmzWPnzp1s3bqV66+/nunTpzN58uROa8fIkSOZNWsWt956K1u2bGHHjh3cdNNNTXoiZsyYwdSpU5k7dy5fffUV2dnZfP/999x///1s377dr/MYDAbuvfde/vjHP/LGG2+QmZnJ5s2befXVVwGYN28eMTExXHbZZXz77bdkZWWxYcMGbr/9do4ePdpp1xtsEmz0ImWF+eQd2sfRw/vZvfFTAHZt+Iz1771Gxp4tfPfx293cQhEMtU43acer+OTH47y0sb6C6a7cciprXd3dNCE6XWpqKjt37uT888/n7rvv5pRTTuGiiy5i3bp1vPjii6hUKj766CMiIyM599xzmTFjBqmpqbzzzjud3pbXX3+dxMREpk+fzhVXXMEtt9xCXFyc93WVSsVnn33Gueeey69//WtGjBjBNddcQ05Ojnc2iT8eeOAB7r77bh588EFGjx7N1Vdf7c3pMJlMfPPNN6SkpHDFFVcwevRobrzxRux2O+Hh4Z1+zcGiUnpwEYCqqioiIiKorKwM+If68Z5jZBbVBKll3WPRxSPb32f5ahRFwRwRSdSApKC1Je/QXv778jLm3HwPySPGBe08om0xoToGx5gZEmMmMcKIWlaq7fFO5u+av+x2O1lZWQwZMgSDwRCUcwgRyO+ZJIj2IrNuuJ0v/v1sm/s8teCKn77/Kj1obWk8jCPBRvcpqXFSUuNke3Y5+hA1g6PN9Y8YEyad/PcWQvQM8teoF6mpKPdrP7VGw7V/eKz9HRvxp6eirDAfa2V5/aJFJ4Zxtn/9MdkHdnH2ZdcxdPzpQe1NEW1zuDykF1STXlCNSgVxYQYGx5gYHG0mIcLgnaMvhBBdTYKNHs7XDR4Alap+CoMPdz77HgOHjw3oPP70VDxy3QUtttXWVJJ3qJL/LLsXCG5vivCfokBhlZ3CKjtbjpRhCNEwKNpESpSJQdEmwgyy2JQQoutIsNHD+brBAz4DDZVKhaIovPfMg/zy9ofaHd7wFcjs2vAZp190uc+8j3n3LuM/T96Hx+17RsTsG+4g79C+oOeLiMDZXW5vrwfU53qkRJsZFGUiKdJIiEZyxYUQwSMJoj3cjnUft3qDV6nV6I1mYhJTOHP2VWz54n0Kcg7jctg5Z+51XH7bn9t8b38STpv3VBw9vL9JXoi/x4meS6tWkWAxkhJV3/MxIFwvQy5B1JUJooMHD263aJQQHVVbW0t2drYkiPY2vvImJl14KQNShvq8wd/13PvEDxpOZVkRtqoKkkecwr/uvxmXw95mD0WDtnoqOpL3cTLHie5T51HIK7ORV2ZjE2AI0TAw0khylInkSCPRofrubqIIkEajAcDpdEqwIYLGZrMBLSuk+hLUYOObb75h2bJl7Nixg+PHj7NmzRrmzp0bzFP2au3lTTQMkzR8BdDqdDx6fcua+jUVpe3OTGkrkGkt7yPUEk1YZAymMAuFuRl+Hyd6D7vLTUZRDRknegbNeg1JFhPJUUYGRpqklHovoNVqMZlMFBcXExISglotw2Si8yiKgs1mo6ioCIvF4g1u2xLUYMNqtXLqqafym9/8hiuuaL/rvT/yJ2+i4QZviU3gjFm/ZMsX71NRfJxQS/3iP/70ULQ328RXIOOLJTaeB95cz/HsQ/xj4S/8Pk70XlaHm0OF1RwqrM/3MOs1JFqMJFmMJEUaiQ2VYZeeRqVSkZCQQFZWFjk5Od3dHNFHWSyWVlfNbS6owcbs2bOZPXt2ME/R6/lKAPXVK/HAm+vRhISgUqmYesnVuF0utLr6T5j+9FCsXv6Iz16T9gIZX0GKVqcjLDKmzeNE32V1uDlcWMPhwvqeD32ImiSLkYQIIwkRBuIjDJJw2gPodDqGDx+O0+ns7qaIPigkJMSvHo0GPSpnw+Fw4HA4vM+rqqq6sTVdY8TEszi0c5PP1xrnPzQEFlD/qaXx88aa9zQU5tYvQdzWbJPmgUz2gd28/cQfmXPzPa0O7TT0cLQWAIn+w+HycKTYypHi+oWzNGoVsWF6EiIMJEQYiY8wEGGUqbbdQa1WSwVR0SP0qGBj6dKlPPzww93djKBo3ENgjojyDp0cO5LW6jGB5D/46qHIO7SXlY/f02Q/X70mzQOZHz57h4w9W9jw/msc3r0Z8B2kND/ueHa6lDAXuD0KBZV2Cirt7KICqB96GRB+IvgINxAXrscQ4v+nIiFE79ajgo3FixezaNEi7/OqqiqSk5O7sUWdp3EPwbcfvtnp7++rp2HbV2t49+k/+zXbpHHuyPa1a4D6AKOBPwmnUsJctMbqcDfp/QCwmEKICzMwIFzPgHADsWESgAjRV/WoYEOv16PX951pdq0lf86+4Q6+fPN5PB4ftTNUKiIHJOFy2APOf2je05AwZDiJqaM5enhfi32b95q0WjysmbaClMbXmDxiHN+seYOL/28Bp0z1771BFnjrTypsLipsLm/iKUCYQUtsmJ7YMD1xYXpiQw2EG7WSgCpEL9ejgo3ezNdNsrXkz8///Uyr73PX8x+QNGxMp+Q/bFv7kTfQ8DVrpHGb26sO2sCfIKWmopS3n/gjAK8t+V1ABb6kd6R/q7bXUW2va9IDotOqiTTpiA7VEROqI8qsJzpUR5heghAheougBhs1NTVkZPxUiyErK4vdu3cTFRVFSkpKME/d5XzdJNu6gavUahSPx2cQ0FYCaHt89TSoVCrikodyyrQLObBlAzUVpYRaovnfu69423z5bX9udUZLw3v4mtraXpCiN4Vy9PD+NouLBVo2XfQvzjqPd52XxkI0KiwmHVFmHRZTSP1XY/33MhwjRM8S1HLlGzZs4Pzzz2+x/YYbbmDFihXtHt/Ty5U3vkn+6/6bqKkoI9QSzaW33OsdQrDEDPB5A7/5kX+x6u9/ajF19K7nP8AS69+8ZV/8KUF++9Or0IbomrT5lkdfpiA3k7cfv6dFYDHjmltJ3/k9pcdziUtO5fLb7m/S6+BvCXPwnevRkbLpQrRFH6LGYtQRYQzBYgohzKAlzBBC+ImvOm33Tc3tinLlQvQ0sjbKSfDnJrlo+WqeWnBFix6MRctXEz9ouDehU1GUThk6aWstFbVG0+4wCUDyiHGcMeuXbP7sXSpKCli0fDURMQNY/fxf2PTft1usu+JPsNGQ6zHpwksDbnNrxwnRUYYQDWEGLaF6LWa9FrNOg1mvJdSgxaTTYArRYtCp0Ws7v4dEgg3RH0nOxknwZwihurwEU7iFiOgBnDXnV02KX7VVO6OjiZLtFfgqzM1s88Z+1Z2PcPrFl3tntJTk51BdXkJNRSl7vv0CaLvC6YiJ01i36iWf5x44fGzA679I+XMRDHaXG7vLTXG1o839tGoVRp0Go06DXqtBr1VjCKn/qteqOTXZIkM2QvihTwYbXdVX09ZNEsBhq+HlP98CgK2qgmk/v8bv4leNc0CADgUevvJBArmxq1Qqlv5mZov92qpwmp9xgHWrXmq1jHlH1n8RorvUeRRv0qovo+LDJdgQwg99MtiYMwfWrotHZ/CgN3kwGD3ojSe+N3swhXowhroxhdV/NYZ6MIe7CbW4MUfUYQrz0JnrFjWeLtpW8mdriZLV5aVk7NnCdx+/zbV/WNru+dorQd7Anxu7vyvDNlyTr3OXFRylqryEo4f3d3j9FyGEEL1Xn8zZOPdc+Pbbjp9XpVYwhdUHH+GRbsKi6giPqiM8qv77iJg6LCceNZXH+cfCX7Q5hLBo+Wq/hgL8zQHxZ4ZGndPZaj5IRXGBt83+JKe2lpPR2nU1P/fdM0e1e11PfZXeZpuF6Il+c9YQIkyBlWKXnA3RH/XJno2PPoL3fyjkUF4tjlo1dpsaR60ah02N3arGVqOhtkaNrfqnr7YqDTWVGmprNCgeFdZKLdZKLYXtLJgYahlCRHQu4dF1lBUWAgYgG8gCcoBKv9vtT62L9qp4NmgtH6QhZ2L+g88xeMxpAa1r4u8QR/Nzt3ddIyed3WabhRBC9G59MtiIjIQBSW5qQgJf7dBdB9YTgUdNhYaqUi3V5VqqyjRUlWnrHyVayou11DnV1FRoqanQkp8JEAY83fQNVWWsekpLfArEJLmISXQSk+QiNsmJOdzTZNf2ckAaNK/iGYiGnIndGz9jyNiJ9U1s58Z+skMccclDWq1kCpCfmdZuLQ4hhBC9V58MNk6GRgvhUW7Co9qeIqooYKtWU1EcQkWxlooiLWWFIZQe11BepKOsMARrpRaUKI5lwrHMlu8RaqkjLtlJXLKTASe+KooRaL2IFgQ+Q+Nki2ad7AqvjSuZ+uLPuitCCCF6Lwk2OkilAnO4B3O4g6ShvqfPvfv0k2z+bDOjp9xE6rhrKTmmo+RYCIU5UFNh8vaKHNlranTUQKCKEH0mMYnVHDvyHrAP2INKVd6hGRqtlRQP5AYf6BBHa5VMw6LiqC4r8nkdzXtsGk+ThY7NyBFCCNH9JNjoZI1vsvu+/wAoI+/QX5h9QzIjTqvvRdjw/ut899EaJpx3F2PP+B1FeTqK8nQU5OooPqrD4w7DaZ/AsSMA53jfW60pANUetn41mtLjWpKGOYhOcNHe8hD+zijpTL4CHEVRqCotbPWYa+5e2qR4V+NpsoqCrJkihBC9VJ+cjQJdU0HUF39mlIRaopqUCW88lFHnguKjOgqy9RTk6Dh2JISCHAOlx333JBhD3Qwc5iBpmJ3k4Q4GjrATk9gyAAl0RsnJaqsqaPN1YRqcdt4lXHLj3d5g7Z+Lf42tqhJjaAQAtTWVmMIt/Hbpa5LfIXoEmY0ihH8k2Ohkbd1k29PWUIbdquZYlo5jR/Qcy9STn2ngeJaOOlfLgiDGMDcpI+2kjLQzaFT914riH1stmx6sCp2tBTg3P/Iv3n7iXkIt0Zwy7ULWv/cKHrcbc7gFa1VFQOe46/n3ZXhFdBsJNoTwjwyjBKi9MuJtzShRqzV4PB0byjCYPaSeYif1lJ9WvnTXQUGOnrxDeo4eNnD0sJ5jR/TUVmtI324mfbvZu68lNhGt7l1CLYeYfGEMB7f/k8qSo11SNKt5gBMWGYO1qhxrVTmFuT+tChxooDHr+ts7pdKqEEKI4JJgI0DtldturPlN9tp7HmPl4/e02K+j639otJA09ESC6uwqAOpccDxLT+5BA7np9Y/CXD0VxUbgSiqK4Ov/gM5wMykjbPzwmYPUU2oZPKYWvbFzO7namjLbVh5JwzBLewpzMzi8ezPQsUqrQgghuoYEG34IdOpoazdZU5gFCO76H9oQSB7hIHmEg7NOFBSrtarJPWgg+4CB7ANGctIM2G0aMn4MJePHUADUaoWkYQ6GjK0ldVwtqeNshEa0f8NvS1tTZtvqAfrVPY+z0sdS983t2vCZ9/uailJ2b6x/vu2r1Zxz2f9JXocQQvQQEmz4IdCpo63dZGsqy7pl/Q+j2cPISTZGTrIB4PFAYY6O7ANGjuyrf5QXhpB3yEDeIQPfrIkEIGGIg6HjbQwdX8vQcbWEWgLPQ/Fnymzz4MsUFuH9OY076yK+eONZPG7fC2G1Rup2CCFEzyEJon5oK+mzId+i8ZTNtvTU9T/Ki7Rk7TdyZG/9oyBH32Kf+MEOhk+wMfw0G0PH1WIMPbmej7bWaAmNiPL+nFwOB7s2fMq7T/85oMTbQP9thAiUJIgK4R/p2aD94lFtldtOGjqauOQhfp+rp67/ERlXR2RcNRPPrwagpkJD5l4jmT8aydhjqp+Ke+Lx7YeRqNQKycPtDJtQy4iJVoaMtROiCyxu9bcyaYhez5SZV5CYOrLdUu6NdTQXRgghROeSYIP2i0c1LrfdvMs/79C+PlloKtTi5tRzajj1nPreoZoKDZk/Gjm828Th3SaKj+rITTeSm27kf+9EEaL3kHpKLSMm2hgx0UpiqrPdYmPQseCreS7HlIuvYOtXq1t9PRDtzTYSQggRuH4bbDRO+ty5/mMAtn/9sff1LV9+wKGd3zPhvEvYuf6/QP1NLDohhSGnTCJr3w5Kj+eiKEpA64z0VqEWN6eeW8Op59YHH+VFWjL2GDm8y8yhnSaqyrSk7zCTvsMMxBIWWVefJzLZysiJtg7le7RsQ9PE282fvUtFSQHTfn4tadu+ISJmAGfOvuqkcmECmW0khBDCP/02Z8OfSp8d1d8SEhUFCrJ1HNppIn2nmSM/GnE6fio2plIpDBzuYOQkK6NOtzJotB2NpmPnai3nxZ9cmNZ6LRoHnv+6/6ZWq7sK0ZzkbAjhn37bs9FWnYeOCtY6Iz2dSgUJQ5wkDHEy/RcV1DlVZB0wkL7DzMFtJo4dMXhnunz9n2iMoW5GTrIx+vT64CMs0v9/g9aGXRpvP3p4n8+gorVei85YqE4IIUTr+m2w0Vadh/bMu3dZpxbn6mu0OoXhE2oZPqGWn98IVaWa+sBju4n0HWZs1Rp2bwxj98YwAAYOtzPmDCujp1hJHmFH3bICe0AaBxXmiKh2a6R0x0J1QgjRn/TJYGP79u38+bY7mHHDIr/G3f1PKFQBP+0XzOJcfUl4tJvTL67i9Iur8LghN91A2lYzadvMJ8qs1z++eiuaUEsdo6dYGTPFyshJNgxm/6bXtlZ47dsP32yxr69ei9YCTwkghRDi5PXJYOONN95g77bvsSQNbRJsNB+zb5xw2FA8SlE8Pktlz7jmVtJ3fk9F8XFikwZ3S3GuvkCtgcFj7AweY2f2/FKqyjQc3G4mbYuZgztM1FRo2fZVBNu+ikCjVRg6zsaYM62MPbOG6ITWC3u1NhTSdlta9lpIACmEEJ2vzySI5uTkUFJSgkqlYvbs2RQVFbVI8tvw/ut899GbnDP3Oi6/7c9A04RDl8NBfsYBnr3rmhY3m0XLV5M0bExACYkiMHUuyNpn5MBWMwe2hFJ8tOnPM36QgzFnWjllag0po5oOt3Rktd2GFW/zDu1lzfJHKMrPJjo+uUWBMUtsfGddouhjJEFUCP/0mZ6NwYMHe79XnSjw0Ly7PNQSBbS+rkmIXo8lLqHVXovWEhJ7UnGu3kwbAsNPq2X4abVcdmsJRUdDOLA5lP2bzWTtq69qWpCj53/vRBFqqWPsicBj+ETbSeXgbFv7Edlpuzlrzq+4YuGDbRYYE0IIEbg+07OxcuVK5s+fT11dYGtoQMuZBtJr0fNYq9Qc3G5m/w9mDm4zY7f9NHc2RO9h5EQbSUPT+fKtmUDrwycqlYrIAUk4bFau/ePjhEfGyHRX0WHSsyGEf7ok2Fi+fDnLli2joKCAU089leeee44pU6a0e1yg/yl37tzJpEmTWmxXqzV4PCe/ronoGepckPmjif2bzez7PpSK4sZ/7N3ojbtIGZXF4V33AllNjm0YDrt75qh2zyPTXUV7JNgQwj8nOcmwfe+88w6LFi1iyZIl7Ny5k1NPPZWZM2dSVFQUtHOqTwzmNwynXHuP76mLdz77ngQavZA2BEZOsnHFgmIeeCuLu1/M4eL/KyVpqB3Q4KidzOFdVwJHgN3AQ8AE7/EqlYpZN9ze6vurNRrm3bssmJcghBD9StBzNp566iluvvlmfv3rXwPwz3/+k08//ZTXXnuN++67r1PPFRcXR3x8PMnJyUye+Qs+eXclFcXHMYVZAJlp0BepVJA01EHSUAezri+lrEDLvh9C2b1RR/aBUODUE48lqNQ5bPqvhkkzFKrLK1p9T5nuKoQQnSuowYbT6WTHjh0sXrzYu02tVjNjxgx++OGHTj/fwIEDyc7ORqfT8d8fjzPsnLm4XS5qKstkqmo/ERVfx7mXV3Du5VBZcpxDOy3s/T6M9B0mXI5BbPkCtnwBqP4OnAusBv4H2OsjFwlChRCi0wU12CgpKcHtdjNgwIAm2wcMGMDBgwdb7O9wOHA4HN7nVVVVAZ9Tr9d7v2+YJeLvUuaib4mI0XD6xdWcfnE1TruK+y69B5gLXApKDHDjiUc18Bkoawi1bJEgVPgl79Be5jxyK0/9fRmTJ0/u7uYI0aMFPWcjEEuXLiUiIsL7SE5O7rT31up03hwOmara/+gMCvPuPRu15iZgAHAh8DxwFAgDrgZWUVuTyfvPTmTLF+HUVPao/x6ih9m29iO+/WYDb77ZskqtEKKpoP41jYmJQaPRUFhY2GR7YWEh8fEtCyUtXryYyspK7yMvLy+YzRP9zKQLL+XOZ98D6qgfOvk9kAKcwekXHSI2yYm7Ts2BLaG881Q8S64eyvJ7BvLthxbKi7q2JE3eob28cM/15B3a26XnFW0rK8wn79A+jh7e7y2Lv2rVKnbu3MmOHTvIycnp5hYK0TMF9S+oTqdj0qRJrFu3jrlz5wLg8XhYt24dCxcubLG/Xq9vMgwiRLD8lCwMirKVc+bu55o/QGGOjh+/C2Xv96HkZxjI3GMic4+JNS/EkTzSzvizqhl3Vg1xya6gtq+1FWpF9/JVFr+4uLjJlHtJPheipaB/XFu0aBE33HADkydPZsqUKTz99NNYrVbv7BQhulLj9XBaVoiF+MFO4geXcfH/lVF6XMveTfWBR/Z+I3npBvLSDXz6WiwDBjkYf1YN486qIWmYgxMjdC00X4+nLa0tJuer2q3oHr5WCG4ILrRaLStWrOimlgnRs3VJUa/nn3/eW9RrwoQJPPvss5xxxhntHncyxW8+3nOMzKKajjZZ9GEdqRBbVaZh3/f1gcfhXSY87p+ii8gBLm/gMXhMLeqfipuyevkjLdbjac2ii0e223YpNNb9jh7e77Ms/o4dO5g4cWK7x0tRL9EfdclA9MKFC30OmwjRHTqyrk14lJtpP69k2s8rsVWrObDFzN7vQzm4zUx5YQgbV0eycXUkoZY6hk8oYtDoLFJGlgTUQ+HrU3MDXyvUis4RSO9TYw3BqlqtxuNjpWghxE/6zEJsQnQVU5iHyTOqmTyjfkpt+g4TezeFsX+zmZoKLbs2JLJrQyJQRf2MlzXUVHze5NOwrx6KthaTk0JjnatxgBFIfkzDCsGmcAvR8cksWvhb3n5zBXl5ecTFxXVR64XofSTYEOIk6AwK486yYondTHnRU5x67mPs3RTG4V1xQCJw7YmHHfgalepjrlgwod33lWq3wfXtRyvJ2LOFDe+/xuHdmwH/ep+arxD867NTuWPh73A6nZLcLkQb+m2w0dGuU9H/+PO7sm3tR2T++D2JqS/zu8f/TF76Fv7x+8eBK4DLgeHAz1GUn7P6BYVdG2o5ZVoNp0yzEpP408yWthJYxclpnIC7fe0aoD7AaFBTUeqz98lX4u6eb7/kjFm/ZLe5gsED4xk0aFAXXokQvU+/DTZkaqHwV2u/K23NHik8mglsQaXaiqLcC4wF5jIg5W4KcyM5ss/EkX0mPv4XxA92MG5aDadMq2HgcKl2Gyy+pq360pAf0xBkZuzZ0mKfhsDkqRPPpfdJiLb1q2BDphYKf/nzu+Lr5tX80/HA4ac06qF4iVuXXoyiDKyf2bIplCN7jRRk6ynI1rP27WgsMS7GTrUydmoNw8bXotUhgUYnaSsBt7GG/JjVyx8hY88WRk46i8O7N/s8Tqa7CuGfLpn62lGdMfW1cRf4Pxb+st3jZGqhAP+mobZ387p43m3MvP72FlNsG/9ORsWfStpWM/t/CCVtmxmn/aeivnqTm1GTbYw9s4bRU6yYw2XGw8lqbdoq/JQnM+/eZQxIGcq/7r+JmooyQi3RXHbrfax8/J4Wx2z4bgvTz5oSUBtk6qvoj/p8z0bjLnCZWij85c/vSluzRwBqrdU+1+Np/Dt5+W3jvDNbXE4Vh3eb2Pe9mQObQ6kq07LnmzD2fBOGWq0weGwtY8+0MvbM4Fcw7euaJ97OuOZW0nd+T96hvS2CipqK0hbbJHFXiMD0yWAjJyeHjANp5JfXtugCv+YPj/G2j08oMrVQNNZWIJE0dDRxyUPafY/Gwy4ul4OQEH27Q3hjplgZM8WKx1NE3iED+38ws39zKMez9BzZa+LIXhP/fTmW2CQnY86sYeyZVoaMrUXTJ/8nB8afRN7mCbibP3uXipICps35FbN/fRfbvlrDu0//udXeqqj4gVxw1c3exN3YuNhgXpIQfUafHEZRtVY7utk+jacWLlq+WoIN0URDl7uvaagNFUErigv4x8JfUF1ectLna2sIr/S4lgNbQtm/2UzmjybcdT/9jhvMbkZOsjFmipVRp1sJi2w7J6Gv8rdaa3sVZFsbarl60d+YMvOKJsfdcv5IIkwhAbVThlFEf9QnPw+99dZb3HDDfNzuuhavqdRq9EYzsUmDZWqhaFPDp+BQSzRjzjiPfd+voygvE0VRmvRKLPj7W+Sm72VVG8MuZ86+is2fv9vhIbzohDrOmVvBOXMrsFvVpO8wcWCLmQNbzVgrfxpuUakUkkfaGX2ihyRpmAN1UNd27l7tJfJWlhbxzeoVTXo7/K0g23yo5NDOTZwx6xftHieEaKlP9mwA/GPVlyy6dlaL7YuWryZ+0PCA18YQ/VOd08kff97+1Oinvkpv9RNxQ69Ze693pPaLxw25hwyknQg88jMMTV4PtdQxarKNUadbGTmp7yWZ+pPIC/i1Nk2DiuIC/v67uYRaojll2oWsf+8VPG435nALty59rcmw12/OGiI9G0L4oU/2bDTmqwu8I2tjiP6l8Y0/0MTi9qp/Nn/9vWce5Je3P9Sh2i9qDQwebWfwaDuz55dSWaohbauZtG1mDu00UVOhZfvX4Wz/OhyVWiF5hJ1Rk22MnGQlZZQdjab9c/Rkbf3bqNRqdAYTDltNQFPcLbHxWKvKsVaVU5ib4d1urapot+S8EMK3PhtsRERJJUbRcU1njPzZrzVL2qv+6ev149mHyDu0L+Cy2a2pKt3Nzv/VB0nXLR5H9gEjadvMHNxm5niWntyDRnIPGvnqrWgMZjfDT7MxapKNEROtRCe0HHbs6dpK5FU8Hhy2+pWfW6sO2hqZuSZE5+qzwygf7zlG+tEyGS4Rfms8/t+4xsItj75MQW4mbz9+T7uJxe0lH9Y5nVSWFWGrqgDgHwt/0W67AvkE3VaSZEWxlvQdJtJ3mEnfaaK2umm3RnSikxGn2Rgx0cbwCTZMYT1/yCXv0F7ee+Yhjh7e5/d01JGTzubWpa+2u197w16ADKMI4ac+27MBMlwiAhN4RdCWPWXt/c5pdToevf5Cv9rj7ydofyvjWmLrOGNWFWfMqsLjhrxDBg7uMHFoh5mcNAOlx3T8cEzHD59aUKkVBg5zMGxCfeAxZGwtemPP+1yybe1HHD28jxC9gfhBw73/NqXH87BVV/g8Jj8zjaOH9/vdcySL4glx8vp0z0ZmUU2QWib6oh3rPm6z6/yqOx/h9IsvP+mesrbO05i/07H9SZJsr3fEblWT+aORQ7tMHNppojC36QqmGq1Cyqhahk+oZeh4G4NG29HpO/9Phz9Jsr56oMwRUdzy6MsAmMItVJeV8OydV/t1ztZ+Ng3TmpsPi931/AdYYuMB6dkQwl99umdDiEC0Nf7fvOjbyfSUtVd5NNBP0J2RX2Awe06syWIFoKJES8ZuI4d3mzi8y0RFcQhZ+0xk7TMB0WhCPKSMtDN0XC3DTq1l0OjO6fnwJ0nWVw+UtbKsyZDUgys3EhYZQ4jeQHnRMRRPyyGh9n42llhZFE+IziLBhhA+dFXXeWtlswNJZg4kSPKXJabOW0ZdUaD0eAiHd5vI2G0k80cTVWVab/Dx9X9ArVEYOMxO6rhahoytZchYO6GWtntuGnoxpv9iPuFRcW0OA1krywKaHdQ4UMjPONDhn40MxQrROSTYEKKR9maUBOs8zctmd/QTdDCCJJUKYhJdxCRWMvVnlSgKlBwLIfPH+sAjfYeWmgoTuelGctONbHi//ri4gU4Gj61l8JhaBo+2E5fibFJgrKEXo60l3Bucfdl1Ac0Oajwco1Kpg/azEUL4R4INIRoJVtd581yE9s4T6Pm6KkiC+uAjNslFbJKLM2dX8cHzj7Dp428YOek+ouKvIWufkYIcPUVHdRQd1bH1ywigvqx6YmolcQMLSRhSzq71PwCgN4Xistfi8fiulTHrut/z7UdvAj/1eBTkZp5oi+8AovFwzPlX3iTT4IXoZpIgKkQX8HfdDn/5SqRsb9ptZ2prmrCiKKg1sVQUpZKdZiA7zUjeQQNOh6+66UeArSce24DdgP//b5NHjPMGEGUFR7n2j48THhnTok0up5Pw6Fii4wd26s9GEkSF8I8EG0IESXs35EALdjXW2cFLoAKdAeN2w/Ejer758DDbvy4AZTIwysdRHuAQsAPYiUq1G0XZCVQ02cvX7KC7Z/p6v9bb1Bkk2BDCPzKMIkSQ+FO3I5Cbn7/1NNrTkTVYmh8X6AwYjQYGDnfwq3tSOHduNU8tGA1EAJOBKcAUVKrTUZQk6oOQUcA8fvoolAPsAX4E9nDD/bcwduogGhZ4VqlUUvVTiB5Mgg0hgqSzb36dFbx0ZA2W5sf5W8K9LSpVFYqyDpXqfyiKwu//8Q5R8RPJzzRw9JCevMM6cg9qqSozA4NOPC4F4PW/gFbnYUCykwGDnMQPchA/+FrmPzCe1x6aTX0PSeBtEkIEhwQbQgRJZ09JPZngJZBekcY9GOaIqFaPay9JszWtJbNaYuMJj/IQHmVj9Ok2oL6w1lO3/Rpj2LkkDb2Kw7ts1FqHoNacSp1TQ36mgfzMxivdJgE24DBwkPohmTSOZUViiVNjDvd4e0OEEF1Hgg0hukBnTLs8meAlkF6Rxj0Y3374ZrvHtVfCvblAZvxYYuN5cOVHLRJf1ZpMygpDKMjWUZCjpyBbx/EcPUV5IbhdeuCUE496q56s/2owu09M43USk+gickAdUQNcRMW7iIyt43j2jx0aYhJCtE2CDSGCqDOmpHZGzYj2ekXm3HwveYf2tejBmH3DHXz55vM+p6U2T9IMZJpwIMWyWtu3Pmhwcco0q/d1jxtK8hVKC0wUHdVTlKujME9L2XE9FSUh2K0ajh7WcPSwocV5AHSGaJz2RN78m4YxZ8Riia0jIsaFJbaO8Cg3YZF16Aw9NqdeiB5LZqMI0UH+Jlqe7JTUxjNPzr/ypnbX62hNW6uYtlY6vS3+rt3SUzgdKsoKQig5FkJJfgilx3WUF2kpPgoVxQZcTv8+exlMbsJOBB4TRhp47G9qhg/3vx0yG0X0R0Hr2Xj00Uf59NNP2b17NzqdjoqKimCdSohu4W+iZUdKXreVYzH/wee9NSPa6k1oLRjy1SsyYuJZHNq5qc02NT/uvWce5Je3PxS04YaOzpppjU6vED/ISfwgZ5PtP03jjQaGAAObPZKJij+T6nItLocau02D3aah+KiOI3vhkb+cdNOE6POCFmw4nU6uvPJKpk6dyquvvhqs0wjRpTpr+ml7/M2xaCt4aR4MNR/S2fTft6ksKaS6vJRjR9LabE/zJdwLcg6Td2ifXzNaOho0dHTWTKB+GmIqBUqB7d7XGhJvJ10Yj6KAw6amqkxDVbmW6jINoyMGkJysCVrbhOgrgj6MsmLFCu68884O9WycTHfj1wcKyS614vYouBUFRaH+e0+PHTUSvUBnLOfuj/aWu6+/AV7a4rX2ConpTWZikwajUqn8upYG5ogorrn7UWprqjGEhvHO3+/HWulfkbJACpD5Uwit8aJsnRWEtDXE1NZQkRT1EsI/PSpB1OFw4HA4vM+rqqo6/F4zxgzwuV1RFFxuhTqPp/6r20OdR8FZ58FR58FR58ZZ5/E+r3W5sZ941Drd2Os82F1uem6miwimrioc1dGZJ4HMOmnrWpqzVpbx6oO/8/u9O9oD5E/7Gy/K1tk9HrJYmxDB0aOCjaVLl/Lwww8H9RwqlQqdVoUOX+s0+MfjUbC53NgcdVidbqyOOmxONzUOF9X2OqrtddQ46qh1tv9HXPQuwVjOvT2B3ADbCyBGTjrb+31b16JWa3zOQGlN80CrowXI2mp/a4uydcbwVVcuZCdEfxRQsHHffffx+OOPt7lPWloao0a1v0aBL4sXL2bRokXe51VVVSQnJ3fovYJJrVYRqtcSqm/7x1fn9lBtr6PK7qKy1kVV7U/fV9hc2F0SjPRmwf4U3JEbYFsBBEB+ZhpHD+/33qBbu5Zr73mMlY/f0+L4efcu87m9eaDV0R6gttqveDx8/u9nvM9PpvR7c8Fa7VcIUS+gYOPuu+9m/vz5be6Tmpra4cbo9Xr0en2Hj+9ptBo1kWYdkWbff7DsLjflNicVNpf3a6nVSYXVSZ3klvRYXfUpuPkNcODwMfz35WVUlxe3O821Nc1v0A+u3OjzWkxhFqD1gKq9QKszeoCan6O13pbOGr7qyKwhIYR/Ago2YmNjiY2NDVZb+h1DiIaECCMJEcYm2xVFobK2PvAoszoprXFSUuOgzOqUBNceoCs/BTd+z+1ff0zmj1vbzVVoCIZC9AbKi46heDwt9mm4Qbd2LTWVZT6DkNikwa0GWoFMtW1La8HcNXcv5eU/39xif1n3RIieL2g5G7m5uZSVlZGbm4vb7Wb37t0ADBs2jNDQ0GCdtk9QqVRYTDosJh1DG8V2Ho9Cuc1JyYngo7jaQVG1HatDhmO6Wld9Cu5IomXjACI/40C7vQu+rqWtgKq17f9795U2p9qebDnzgpzD3jZKEqcQvUvQgo0HH3yQf//7397np512GgDr16/nvPPOC9Zp+zS1WkV0qJ7oUD0jCfNur3HU1QceVXaKqh0UVtmpttd1Y0tFZ+loomXz4KcjN+jWAqrG28uLjvkMhpJHjOObNW9wxcIHGX/2xd6gIfvAbt5+4o/tTlv1dW5J4hSi9+qz5cr7O6ujjsIqOwVVdoqqHBRU2WV2TC/kT72NuOQhrdadqCgu6HB5c38EWnckkJobvpxs6ffOJnU2hPBPj5r6KjqPWa8lNTaU1NifhqzKrU6OVdZSUGnneKWd0honnp4bawr8S7RcvfyRVutOBDu/pL2ptnpTKAc2r8dWU4UxLPykq65KEqcQvZMEG/1Iw8yYsYkRADjrPBRU2smvqOVYRS0FVXacdS2TCUXP0HwopDA3E6DdG3gwb9DtTbV12Gp45cHfttjefCjoruffl6XdhejDJNjox3RaNSnRJlKiTUB9AmpxjYOj5bXkV9SSX14rtUB6AF+5CnmH9raod9GZdSe6SsNQUFetgyKE6B4SbAgvtVrFgHADA8INTBoUiaIolNQ4OVpu8wYgkvfR9XwNhWz7ag3vPv3noJdN90fjYGjExGmsW/VSi31aKwZ27R8eY0DKUD56aSkQnEXthBDdTxJEhd8Upb7nI6+s1huAyLBL9+no4mHB0JC42TDVtvmQT0OwEeismJ7eOyMJokL4R3o2hN9UKhVxYQbiwup7PjwehYIqO7llNnLLbBRU2qXoWDfoCXUnGvJAWpue6qsYWPHRLJz22qBWBRVC9AzSsyE6jbPOw9FyGzllNvLKbJTWOLu7SX1asKe1dlRr01N9bS/IOdxjemc6Qno2hPCP9GyITqPTqptMt62yu8gpsZFTZiW3zIbDJUMunamnLh7mTzGw5rNiekLvjBAieCTYEEETbghh3MAIxg2MwONROF5lJ6fESlapleJqB3JPOXknM621tbVMupJUBRWif5BgQ3QJtVpFksVIksXItGExWB11ZJdayT7R8yG9Hl2vJ0w37am9M0KIziXBhugWZr2WsYkRjE2s7/XIr6glq8RKdqlVcj2CqCMLuwWbVAUVou+TBFHR41TaXBwpqSGrxMrR8lqZ4dKJAl3LRLRNEkSF8I+6uxsgRHMRphBOS4nkiokDuXV6Kj8fn8CYxHCMOk13N63Xm3fvMtQa3z9HtUbDvHuXdXGLhBD9gQyjiB5Nr9UwfEAYwweEoSgKxyrtHCmu4UixlTKrDLcEyp+F3YQQorNJsCF6DZXqpyTTc4bHUm51knki8DhWWSuzWwIk002FEF1Fgg3Ra0WadUw2RzF5cBQ2Zx1Hiq1kFteQW2qjTvI8WiXTTYUQXU0SREWf46zzkFtmJaPISlaJVVau9aG1Kp8iMJIgKoR/pGdD9Dk6rZphcWEMiwvzTqvNKKohs7iGantddzevR5DppkKIriTBhujT1GoVyVEmkqNMnD8qjsIqO5lFNWQU10g9DyGE6CISbIh+ZUC4gQHhBqYNi6HsRIJpRlENhVV2STAVQoggkWBD9FtRZh1R5ihOHxxFtd1FZrGVjKIa8str8UjkIYQQnUaCDSGAMEMIE5ItTEi2UOt0k1lcIzNbhBCik0iwIUQzRp2GU5IiOCUpAmedh+zS+h6PrBIrzjpZME4IIQIlwYYQbdBp1YwYEMaIAWHUuT3kltnIKKrhSImVWqdMqRVCCH9IsCGEn7QaNamxoaTGhv40pba4hswimVIrhBBtkWBDiA5oMqV2ZP2U2oyi+pktsmaLEEI0JcGGEJ2gYUrtWSem1DYUESuotHd304QQottJsCFEJ4sy65gyJIopQ2RKrRBCAKiD9cbZ2dnceOONDBkyBKPRyNChQ1myZAlOp3Qxi/6jYUrtLycN5NbpqcwcG8/QuFBCNKrubpoQQnSZoPVsHDx4EI/Hw0svvcSwYcPYt28fN998M1arlSeffDJYpxWixzKEaBiTGM6YxHBcbg85pTYyi+un1MrMFiFEX9alq74uW7aMF198kSNHjvi1v6yOKPqDhpktR0qsZBbVUFnr6u4mCT/Jqq9C+KdLczYqKyuJiopq9XWHw4HD4fA+r6qq6opmCdGtGs9smT4ilpIaB5knannImi1CiL6gy4KNjIwMnnvuuTaHUJYuXcrDDz/cVU0SokeKCdUTE6rnjNRorI46jhRbOVJSQ16ZDZdbIg8hRO8T8DDKfffdx+OPP97mPmlpaYwaNcr7PD8/n+nTp3PeeefxyiuvtHqcr56N5ORk6W4UAnCdqGCaVWwlq8RKjUMKiXU3GUYRwj8BBxvFxcWUlpa2uU9qaio6nQ6AY8eOcd5553HmmWeyYsUK1Gr/J8DIf0ohfFMUhaJqB0dOBB5F1TLc0h0k2BDCPwEPo8TGxhIbG+vXvvn5+Zx//vlMmjSJ119/PaBAQwjROpVK5S0kNnVoNDWOOrJL6gOP3DKbLBgnhOhRgpazkZ+fz3nnncegQYN48sknKS4u9r4WHx8frNMK0S+F6rXelWrdHoX88lqySq1kl1ilfLoQotsFLdhYu3YtGRkZZGRkMHDgwCavdeFsWyH6HY1aRUq0iZTo+tktlbUuskusZJdaOVpeK70eQogu16V1NgIlY5tCdC63R+FYRS3ZpVayS22UVDvaP0i0SnI2hPCPrI0iRD+iaVTT45zhUOOoI6fUSm6pjdwyGzapZCqECAIJNoTox0L1WsYmRjA2MQJFUSiucZBXZiOn1Maxilqp6yGE6BQSbAghgPoZLnFhBuLCDEwaFEWd28PxSjt5ZTbyym0UVDpk1VohRIdIsCGE8EmrUXuHXAAcdW7yy2vJK6/laLmN4mqH1PYQQvhFgg0hhF/0Wg2psaGkxoYCYHe5ya+oJa/MRn5FrQQfQohWSbAhhOgQQ4iGobGhDG0UfByvtJNfXkt+hY3CKgduj0QfQggJNoQQncQQomFIjJkhMWagfi2Xgko7xyvtHKuo5VhlLQ6X1PgQoj+SYEMIERQhzXI+FEWh1OrkWEUtxyvtFFTaKbc5ZehFiH5Agg0hRJdQqVTEhOqJCdUz/kRRYbvL7e39KKiqpbDKQa3U+hCiz5FgQwjRbQwhGgbHmBl8YugFoNLmoqDKTmGVnYIqO8XVDimxLkQvJ8GGEKJHiTCFEGEKYWR8GFA//FJuc1FUbaeoykFRtYOiarvkfwjRi0iwIYTo0VQqFVFmHVFmHaMaLRhdWeuiuNpBSY2D4ur6R5XdJTkgQvRAEmwIIXqlCGMIEcYQhsWFerc56zyUWh2U1jgpqan/Wmp1YHVIHogQ3UmCDSFEn6HTqkmIMJIQYWyy3e5yU2p1Um51NvlaLT0hQnQJCTaEEH2eIURDksVIkqVpEFLn9lBR66LC5qTc5qLc6qSi1kWlzYXVWSeBiBCdRIINIUS/pdWovdNxm3O5PVTWuqisdVFhc1Fld1FV66LKXkdVrUtmyAgRAAk2hBDCh5A2AhGoH5rRadRd3CoheicJNoQQogMMIZruboIQvYaE5UIIIYQIKgk2hBBCCBFUEmwIIYQQIqgk2BBCCCFEUEmwIYQQQoigkmBDCCGEEEElwYYQQgghgkqCDSGEEEIElQQbQgghhAiqoAYbl156KSkpKRgMBhISErjuuus4duxYME8phBBCiB4mqMHG+eefz7vvvkt6ejoffPABmZmZ/PKXvwzmKYUQQgjRw6gUpesWUf7444+ZO3cuDoeDkJCQdvevqqoiIiKCyspKwsPDu6CFQggRXPJ3TfRHXbYQW1lZGStXrmTatGmtBhoOhwOHw+F9XllZCdT/5xRCiL6g4e9ZF37OE6LbBT3YuPfee3n++eex2WyceeaZfPLJJ63uu3TpUh5++OEW25OTk4PZRCGE6HLV1dVERER0dzOE6BIBD6Pcd999PP74423uk5aWxqhRowAoKSmhrKyMnJwcHn74YSIiIvjkk09QqVQtjmves+HxeCgrKyM6Otrn/q2pqqoiOTmZvLy8PttN2devUa6v9+vr19jR61MUherqahITE1GrZUKg6B8CDjaKi4spLS1tc5/U1FR0Ol2L7UePHiU5OZnvv/+eqVOnBtbSAPSHMdG+fo1yfb1fX7/Gvn59QnSmgIdRYmNjiY2N7dDJPB4PQJPeCyGEEEL0bUHL2diyZQvbtm3j7LPPJjIykszMTB544AGGDh0a1F4NIYQQQvQsQRswNJlMrF69mgsvvJCRI0dy4403Mn78eDZu3Iherw/WaQHQ6/UsWbIk6OfpTn39GuX6er++fo19/fqE6ExdWmdDCCGEEP2PpEILIYQQIqgk2BBCCCFEUEmwIYQQQoigkmBDCCGEEEHVa4ON5cuXM3jwYAwGA2eccQZbt25tc//33nuPUaNGYTAYGDduHJ999lkXtbTjArnGl19+mXPOOYfIyEgiIyOZMWNGuz+T7hbov2GDVatWoVKpmDt3bnAbeJICvb6KigoWLFhAQkICer2eESNG9Pjf00Cv8emnn2bkyJEYjUaSk5O56667sNvtXdTawHzzzTfMmTOHxMREVCoVH374YbvHbNiwgYkTJ6LX6xk2bBgrVqwIejuF6BWUXmjVqlWKTqdTXnvtNWX//v3KzTffrFgsFqWwsNDn/ps2bVI0Go3yxBNPKAcOHFD+/Oc/KyEhIcrevXu7uOX+C/Qaf/WrXynLly9Xdu3apaSlpSnz589XIiIilKNHj3Zxy/0T6PU1yMrKUpKSkpRzzjlHueyyy7qmsR0Q6PU5HA5l8uTJys9+9jPlu+++U7KyspQNGzYou3fv7uKW+y/Qa1y5cqWi1+uVlStXKllZWcqXX36pJCQkKHfddVcXt9w/n332mXL//fcrq1evVgBlzZo1be5/5MgRxWQyKYsWLVIOHDigPPfcc4pGo1G++OKLrmmwED1Yrww2pkyZoixYsMD73O12K4mJicrSpUt97n/VVVcpl1xySZNtZ5xxhnLrrbcGtZ0nI9BrbK6urk4JCwtT/v3vfweriSelI9dXV1enTJs2TXnllVeUG264oUcHG4Fe34svvqikpqYqTqezq5p40gK9xgULFigXXHBBk22LFi1SzjrrrKC2szP4E2z88Y9/VMaOHdtk29VXX63MnDkziC0TonfodcMoTqeTHTt2MGPGDO82tVrNjBkz+OGHH3we88MPPzTZH2DmzJmt7t/dOnKNzdlsNlwuF1FRUcFqZod19Pr+8pe/EBcXx4033tgVzeywjlzfxx9/zNSpU1mwYAEDBgzglFNO4W9/+xtut7urmh2QjlzjtGnT2LFjh3eo5ciRI3z22Wf87Gc/65I2B1tv+zsjRFcK+hLzna2kpAS3282AAQOabB8wYAAHDx70eUxBQYHP/QsKCoLWzpPRkWts7t577yUxMbHFH7+eoCPX99133/Hqq6+ye/fuLmjhyenI9R05coT//e9/zJs3j88++4yMjAxuu+02XC4XS5Ys6YpmB6Qj1/irX/2KkpISzj77bBRFoa6ujt/+9rf86U9/6oomB11rf2eqqqqora3FaDR2U8uE6H69rmdDtO+xxx5j1apVrFmzBoPB0N3NOWnV1dVcd911vPzyy8TExHR3c4LC4/EQFxfHv/71LyZNmsTVV1/N/fffzz//+c/ublqn2bBhA3/729944YUX2LlzJ6tXr+bTTz/lr3/9a3c3TQgRZL2uZyMmJgaNRkNhYWGT7YWFhcTHx/s8Jj4+PqD9u1tHrrHBk08+yWOPPcbXX3/N+PHjg9nMDgv0+jIzM8nOzmbOnDnebQ0rCGu1WtLT0xk6dGhwGx2Ajvz7JSQkEBISgkaj8W4bPXo0BQUFOJ1OdDpdUNscqI5c4wMPPMB1113HTTfdBMC4ceOwWq3ccsst3H///ajVvfuzT2t/Z8LDw6VXQ/R7ve5/t06nY9KkSaxbt867zePxsG7dulZXk506dWqT/QHWrl3bY1ef7cg1AjzxxBP89a9/5YsvvmDy5Mld0dQOCfT6Ro0axd69e9m9e7f3cemll3L++eeze/dukpOTu7L57erIv99ZZ51FRkaGN4gCOHToEAkJCT0u0ICOXaPNZmsRUDQEV0ofWKKpt/2dEaJLdXeGakesWrVK0ev1yooVK5QDBw4ot9xyi2KxWJSCggJFURTluuuuU+677z7v/ps2bVK0Wq3y5JNPKmlpacqSJUt6xdTXQK7xscceU3Q6nfL+++8rx48f9z6qq6u76xLaFOj1NdfTZ6MEen25ublKWFiYsnDhQiU9PV355JNPlLi4OOWRRx7prktoV6DXuGTJEiUsLEz5z3/+oxw5ckT56quvlKFDhypXXXVVd11Cm6qrq5Vdu3Ypu3btUgDlqaeeUnbt2qXk5OQoiqIo9913n3Ldddd592+Y+nrPPfcoaWlpyvLly2XqqxAn9MpgQ1EU5bnnnlNSUlIUnU6nTJkyRdm8ebP3tenTpys33HBDk/3fffddZcSIEYpOp1PGjh2rfPrpp13c4sAFco2DBg1SgBaPJUuWdH3D/RTov2FjPT3YUJTAr+/7779XzjjjDEWv1yupqanKo48+qtTV1XVxqwMTyDW6XC7loYceUoYOHaoYDAYlOTlZue2225Ty8vKub7gf1q9f7/P/VMM13XDDDcr06dNbHDNhwgRFp9Mpqampyuuvv97l7RaiJ5Il5oUQQggRVL0uZ0MIIYQQvYsEG0IIIYQIKgk2hBBCCBFUEmwIIYQQIqgk2BBCCCFEUEmwIYQQQoigkmBDCCGEEEElwYYQQgghgkqCDSGEEEIElQQbQgghhAgqCTaEEEIIEVQSbAghhBAiqP4frt8/2c3zMSYAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "model.eval()\n", "with torch.no_grad():\n", " untrained_pred_dist = likelihood(model(test_x))\n", " predictive_mean = untrained_pred_dist.mean\n", " lower, upper = untrained_pred_dist.confidence_region()\n", " \n", "f, ax = plt.subplots(1, 1, figsize=(4, 3))\n", "# Plot training data as black stars\n", "ax.plot(train_x, train_y, 'k*')\n", "# Plot predictive means as blue line\n", "ax.plot(test_x, predictive_mean, 'b')\n", "# Shade between the lower and upper confidence bounds\n", "ax.fill_between(test_x, lower, upper, alpha=0.5)\n", "ax.set_ylim([-3, 3])\n", "ax.legend(['Observed Data', 'Mean', 'Confidence'], bbox_to_anchor=(1.6,1));" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Visually, this does not look like a good fit. Now, let us train the model hyperparameters." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Iter 1/80 - Loss: 1.700 lengthscale: 0.693 noise: 0.693\n", "Iter 2/80 - Loss: 1.670 lengthscale: 0.644 noise: 0.644\n", "Iter 3/80 - Loss: 1.636 lengthscale: 0.598 noise: 0.598\n", "Iter 4/80 - Loss: 1.597 lengthscale: 0.555 noise: 0.554\n", "Iter 5/80 - Loss: 1.551 lengthscale: 0.514 noise: 0.513\n", "Iter 6/80 - Loss: 1.496 lengthscale: 0.475 noise: 0.474\n", "Iter 7/80 - Loss: 1.434 lengthscale: 0.439 noise: 0.437\n", "Iter 8/80 - Loss: 1.368 lengthscale: 0.405 noise: 0.402\n", "Iter 9/80 - Loss: 1.304 lengthscale: 0.372 noise: 0.369\n", "Iter 10/80 - Loss: 1.249 lengthscale: 0.342 noise: 0.339\n", "Iter 11/80 - Loss: 1.204 lengthscale: 0.314 noise: 0.310\n", "Iter 12/80 - Loss: 1.170 lengthscale: 0.289 noise: 0.284\n", "Iter 13/80 - Loss: 1.143 lengthscale: 0.267 noise: 0.260\n", "Iter 14/80 - Loss: 1.121 lengthscale: 0.248 noise: 0.237\n", "Iter 15/80 - Loss: 1.101 lengthscale: 0.232 noise: 0.217\n", "Iter 16/80 - Loss: 1.082 lengthscale: 0.219 noise: 0.198\n", "Iter 17/80 - Loss: 1.065 lengthscale: 0.208 noise: 0.181\n", "Iter 18/80 - Loss: 1.048 lengthscale: 0.199 noise: 0.165\n", "Iter 19/80 - Loss: 1.032 lengthscale: 0.192 noise: 0.151\n", "Iter 20/80 - Loss: 1.015 lengthscale: 0.186 noise: 0.138\n", "Iter 21/80 - Loss: 0.998 lengthscale: 0.182 noise: 0.126\n", "Iter 22/80 - Loss: 0.981 lengthscale: 0.180 noise: 0.115\n", "Iter 23/80 - Loss: 0.963 lengthscale: 0.178 noise: 0.105\n", "Iter 24/80 - Loss: 0.945 lengthscale: 0.177 noise: 0.096\n", "Iter 25/80 - Loss: 0.926 lengthscale: 0.178 noise: 0.087\n", "Iter 26/80 - Loss: 0.907 lengthscale: 0.179 noise: 0.080\n", "Iter 27/80 - Loss: 0.888 lengthscale: 0.181 noise: 0.073\n", "Iter 28/80 - Loss: 0.867 lengthscale: 0.184 noise: 0.066\n", "Iter 29/80 - Loss: 0.847 lengthscale: 0.188 noise: 0.060\n", "Iter 30/80 - Loss: 0.826 lengthscale: 0.193 noise: 0.055\n", "Iter 31/80 - Loss: 0.805 lengthscale: 0.198 noise: 0.050\n", "Iter 32/80 - Loss: 0.783 lengthscale: 0.205 noise: 0.046\n", "Iter 33/80 - Loss: 0.761 lengthscale: 0.212 noise: 0.042\n", "Iter 34/80 - Loss: 0.739 lengthscale: 0.220 noise: 0.038\n", "Iter 35/80 - Loss: 0.717 lengthscale: 0.228 noise: 0.035\n", "Iter 36/80 - Loss: 0.694 lengthscale: 0.238 noise: 0.032\n", "Iter 37/80 - Loss: 0.672 lengthscale: 0.248 noise: 0.029\n", "Iter 38/80 - Loss: 0.650 lengthscale: 0.259 noise: 0.026\n", "Iter 39/80 - Loss: 0.627 lengthscale: 0.270 noise: 0.024\n", "Iter 40/80 - Loss: 0.606 lengthscale: 0.283 noise: 0.022\n", "Iter 41/80 - Loss: 0.584 lengthscale: 0.295 noise: 0.020\n", "Iter 42/80 - Loss: 0.563 lengthscale: 0.309 noise: 0.018\n", "Iter 43/80 - Loss: 0.542 lengthscale: 0.322 noise: 0.017\n", "Iter 44/80 - Loss: 0.521 lengthscale: 0.336 noise: 0.015\n", "Iter 45/80 - Loss: 0.501 lengthscale: 0.351 noise: 0.014\n", "Iter 46/80 - Loss: 0.481 lengthscale: 0.365 noise: 0.013\n", "Iter 47/80 - Loss: 0.462 lengthscale: 0.380 noise: 0.012\n", "Iter 48/80 - Loss: 0.444 lengthscale: 0.394 noise: 0.011\n", "Iter 49/80 - Loss: 0.426 lengthscale: 0.408 noise: 0.010\n", "Iter 50/80 - Loss: 0.409 lengthscale: 0.421 noise: 0.009\n", "Iter 51/80 - Loss: 0.393 lengthscale: 0.433 noise: 0.008\n", "Iter 52/80 - Loss: 0.377 lengthscale: 0.444 noise: 0.007\n", "Iter 53/80 - Loss: 0.362 lengthscale: 0.453 noise: 0.007\n", "Iter 54/80 - Loss: 0.346 lengthscale: 0.461 noise: 0.006\n", "Iter 55/80 - Loss: 0.331 lengthscale: 0.467 noise: 0.006\n", "Iter 56/80 - Loss: 0.316 lengthscale: 0.471 noise: 0.005\n", "Iter 57/80 - Loss: 0.302 lengthscale: 0.473 noise: 0.005\n", "Iter 58/80 - Loss: 0.288 lengthscale: 0.474 noise: 0.004\n", "Iter 59/80 - Loss: 0.274 lengthscale: 0.474 noise: 0.004\n", "Iter 60/80 - Loss: 0.260 lengthscale: 0.473 noise: 0.004\n", "Iter 61/80 - Loss: 0.247 lengthscale: 0.471 noise: 0.003\n", "Iter 62/80 - Loss: 0.235 lengthscale: 0.468 noise: 0.003\n", "Iter 63/80 - Loss: 0.224 lengthscale: 0.466 noise: 0.003\n", "Iter 64/80 - Loss: 0.213 lengthscale: 0.463 noise: 0.003\n", "Iter 65/80 - Loss: 0.202 lengthscale: 0.461 noise: 0.002\n", "Iter 66/80 - Loss: 0.192 lengthscale: 0.458 noise: 0.002\n", "Iter 67/80 - Loss: 0.183 lengthscale: 0.457 noise: 0.002\n", "Iter 68/80 - Loss: 0.175 lengthscale: 0.456 noise: 0.002\n", "Iter 69/80 - Loss: 0.166 lengthscale: 0.455 noise: 0.002\n", "Iter 70/80 - Loss: 0.159 lengthscale: 0.455 noise: 0.002\n", "Iter 71/80 - Loss: 0.151 lengthscale: 0.455 noise: 0.002\n", "Iter 72/80 - Loss: 0.145 lengthscale: 0.456 noise: 0.001\n", "Iter 73/80 - Loss: 0.138 lengthscale: 0.458 noise: 0.001\n", "Iter 74/80 - Loss: 0.132 lengthscale: 0.460 noise: 0.001\n", "Iter 75/80 - Loss: 0.127 lengthscale: 0.462 noise: 0.001\n", "Iter 76/80 - Loss: 0.122 lengthscale: 0.464 noise: 0.001\n", "Iter 77/80 - Loss: 0.117 lengthscale: 0.467 noise: 0.001\n", "Iter 78/80 - Loss: 0.113 lengthscale: 0.469 noise: 0.001\n", "Iter 79/80 - Loss: 0.109 lengthscale: 0.472 noise: 0.001\n", "Iter 80/80 - Loss: 0.106 lengthscale: 0.474 noise: 0.001\n" ] } ], "source": [ "# this is for running the notebook in our testing framework\n", "import os\n", "smoke_test = ('CI' in os.environ)\n", "training_iter = 2 if smoke_test else 80\n", "\n", "\n", "# Find optimal model hyperparameters\n", "model.train()\n", "\n", "# Use the adam optimizer\n", "optimizer = torch.optim.Adam(model.parameters(), lr=0.1) # Includes QExponentialLikelihood parameters\n", "\n", "# \"Loss\" for QEPs - the marginal log likelihood\n", "mll = qpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)\n", "\n", "for i in range(training_iter):\n", " # Zero gradients from previous iteration\n", " optimizer.zero_grad()\n", " # Output from model\n", " output = model(train_x)\n", " # Calc loss and backprop gradients\n", " loss = -mll(output, train_y)\n", " loss.backward()\n", " print('Iter %d/%d - Loss: %.3f lengthscale: %.3f noise: %.3f' % (\n", " i + 1, training_iter, loss.item(),\n", " model.covar_module.base_kernel.lengthscale.item(),\n", " model.likelihood.noise.item()\n", " ))\n", " optimizer.step()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the next cell, we reevaluate the model on the test data." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhsAAAEYCAYAAAD1WzSOAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAX+NJREFUeJzt3Xd8U+X+wPFPkjbp3psO9l5SEMGBCAqouDdXQVHwAqKiXvTnHlxQuG4cFxEUUBAE4QKiiBQEZZVVpFAKdFA66V5J25zfH6WhI22T0nR+369XX6Un5+Q8Jy3J9zzP9/k+KkVRFIQQQgghbETd3A0QQgghRNsmwYYQQgghbEqCDSGEEELYlAQbQgghhLApCTaEEEIIYVMSbAghhBDCpiTYEEIIIYRNSbAhhBBCCJuSYEMIIYQQNiXBhhBCCCFsyqbBxueff07//v1xc3PDzc2NYcOG8fPPP9vylEIIIYRoYVS2XBvlf//7HxqNhm7duqEoCt988w3z58/n0KFD9OnTx1anFUIIIUQLYtNgwxwvLy/mz5/P5MmTm/K0QgghhGgmdk11orKyMlavXk1BQQHDhg0zu49er0ev15t+NhqNZGZm4u3tjUqlaqqmCiGEzSiKQl5eHkFBQajVkjYn2gebBxtRUVEMGzaM4uJiXFxcWLduHb179za779y5c3nzzTdt3SQhhGh2iYmJBAcHN3czhGgSNh9GMRgMJCQkkJOTw5o1a/jqq6/YsWOH2YCjes9GTk4OoaGhJCYm4ubmZstmCiFEk8jNzSUkJITs7Gzc3d2buzlCNIkmz9kYPXo0Xbp04csvv6x339zcXNzd3cnJyZFgQwjRJsj7mmiPmnzA0Gg0Vum9EEIIIUTbZtOcjZdeeolx48YRGhpKXl4e3333HREREfzyyy+2PK0QQgghWhCbBhtpaWk88sgjJCcn4+7uTv/+/fnll1+48cYbbXlaIYQQQrQgNg02Fi9ebMunF0IIUQej0YjBYGjuZog2yN7eHo1GY/H+TVZnQwghRNMxGAycPXsWo9HY3E0RbZSHhwcBAQEW1cGSYEMIIdoYRVFITk5Go9EQEhIixcNEo1IUhcLCQtLS0gAIDAys9xgJNoQQoo0pLS2lsLCQoKAgnJycmrs5og1ydHQEynMz/fz86h1SkXBXCCHamLKyMgC0Wm0zt0S0ZRWBbElJSb37SrAhhBBtlKwpJWzJmr8vCTaEEEIIYVMSbAghhGh1OnbsyIcfftjczWg0be16qpNgQwghRIuRmJjIY489RlBQEFqtlrCwMJ5++mkuXLjQ3E1rVm+88QYqlQqVSoWdnR0+Pj5cd911fPjhh1YvARIREYFKpSI7O9s2jTVDgg0hhBC1OnDgADfccAMHDhyw+bnOnDnD4MGDOXXqFN9//z2xsbF88cUXbNu2jWHDhpGZmWnzNtSmrKys2WuW9OnTh+TkZBISEti+fTv33nsvc+fOZfjw4eTl5TVr2+ojwYYQQohaffvtt2zfvp1ly5bZ/FzTp09Hq9Xy66+/MmLECEJDQxk3bhy//fYbSUlJvPzyy1X2z8vL48EHH8TZ2ZkOHTqwcOFC02OKovDGG28QGhqKTqcjKCiImTNnmh7X6/U8//zzdOjQAWdnZ4YOHUpERITp8aVLl+Lh4cGGDRvo3bs3Op2Or776CgcHhxo9Ak8//TQ33HCD6eddu3Zx7bXX4ujoSEhICDNnzqSgoMD0eFpaGuPHj8fR0ZFOnTqxYsUKi14fOzs7AgICCAoKol+/fjz11FPs2LGDY8eO8e6775r2W7ZsGYMHD8bV1ZWAgAAeeughU02MuLg4Ro4cCYCnpycqlYpJkyYBsGXLFq655ho8PDzw9vbm1ltv5fTp0xa1rT4SbAghhKgiPj6eyMhIDh48yKpVqwBYuXIlBw8eJDIykvj4+EY/Z2ZmJr/88gvTpk0z1XCoEBAQwIQJE1i1ahWKopi2z58/nwEDBnDo0CFefPFFnn76abZu3QrAjz/+yAcffMCXX37JqVOn+Omnn+jXr5/p2BkzZvDXX3+xcuVKjh49yr333svYsWM5deqUaZ/CwkLeffddvvrqK/7++28mTJiAh4cHP/74o2mfsrIyVq1axYQJEwA4ffo0Y8eO5e677+bo0aOsWrWKXbt2MWPGDNMxkyZNIjExke3bt7NmzRo+++wzUzBgrZ49ezJu3DjWrl1r2lZSUsLbb7/NkSNH+Omnn4iLizMFFCEhIab2nzx5kuTkZD766CMACgoKmDVrFgcOHGDbtm2o1WruvPPOxunRUVqwnJwcBVBycnKauylCCNEomuJ9raioSDl+/LhSVFTUoOMB05dKparyveKrse3Zs0cBlHXr1pl9/P3331cAJTU1VVEURQkLC1PGjh1bZZ/7779fGTdunKIoivKf//xH6d69u2IwGGo8V3x8vKLRaJSkpKQq20eNGqW89NJLiqIoypIlSxRAOXz4cJV9nn76aeWGG24w/fzLL78oOp1OycrKUhRFUSZPnqxMmTKlyjF//PGHolarlaKiIuXkyZMKoOzbt8/0eHR0tAIoH3zwQS2vjqK8/vrryoABA8w+Nnv2bMXR0bHWY/fv368ASl5enqIoirJ9+3YFMLW5Nunp6QqgREVFmX3cmr8z6dkQQghRxfLly7GzKy8wrVzsSaj4bmdnx/Lly212bqVSz0V9hg0bVuPn6OhoAO69916Kioro3LkzTzzxBOvWraO0tBSAqKgoysrK6N69Oy4uLqavHTt2VBk20Gq19O/fv8o5JkyYQEREBOfPnwdgxYoV3HLLLXh4eABw5MgRli5dWuV5x4wZg9Fo5OzZs0RHR2NnZ0d4eLjpOXv27Gk6viEURalS8yIyMpLx48cTGhqKq6srI0aMACAhIaHO5zl16hQPPvggnTt3xs3NjY4dO1p0nCWkXLkQQogqJkyYQK9evap8IFbYu3cvgwYNavRzdu3aFZVKRXR0NHfeeWeNx6Ojo/H09MTX19ei5wsJCeHkyZP89ttvbN26lWnTpjF//nx27NhBfn4+Go2GyMjIGmW2XVxcTP92dHSsUbhqyJAhdOnShZUrV/LPf/6TdevWsXTpUtPj+fn5TJ06tUp+SIXQ0FBiYmIsar81oqOj6dSpE1A+FDJmzBjGjBnDihUr8PX1JSEhgTFjxtS7AvD48eMJCwtj0aJFBAUFYTQa6du3b6OsHCzBhhBCiFqp1WqMRqPpu614e3tz44038tlnn/Hss89WydtISUlhxYoVPPLII1U+/Pfs2VPlOfbs2UOvXr1MPzs6OjJ+/HjGjx/P9OnT6dmzJ1FRUVxxxRWUlZWRlpbGtddea3VbJ0yYwIoVKwgODkatVnPLLbeYHhs0aBDHjx+na9euZo/t2bMnpaWlREZGMmTIEKA8d6Kh01BPnDjBli1beOmll0w/X7hwgXnz5hESEgJQYyZRRRn7irL2ABcuXODkyZMsWrTI9Jrs2rWrQW0yR4ZRhBBC1ODn50dAQADh4eF88cUXhIeHExAQgJ+fn83O+emnn6LX6xkzZgw7d+4kMTGRLVu2cOONN9KhQwfmzJlTZf/du3fz3nvvERMTw8KFC1m9ejVPP/00UD6bZPHixRw7dowzZ86wfPlyHB0dCQsLo3v37kyYMIFHHnmEtWvXcvbsWfbt28fcuXPZtGlTve2cMGECBw8eZM6cOdxzzz3odDrTY7Nnz+bPP/9kxowZHD58mFOnTrF+/XpTgmiPHj0YO3YsU6dOZe/evURGRvL444/XSIo1p7S0lJSUFM6fP09UVBSffPIJI0aMYODAgbzwwgtAee+JVqvlk08+4cyZM2zYsIG33367yvOEhYWhUqnYuHEj6enp5Ofn4+npibe3N//973+JjY3l999/Z9asWfW2yWL1ZnU0I0kQFUK0Na0hQbRCcXGxYjQaFUVRFKPRqBQXFzdG8+oUFxenTJw4UfH391fs7e2VkJAQ5amnnlIyMjKq7BcWFqa8+eabyr333qs4OTkpAQEBykcffWR6fN26dcrQoUMVNzc3xdnZWbnqqquU3377zfS4wWBQXnvtNaVjx46Kvb29EhgYqNx5553K0aNHFUUpTxB1d3evtZ1XXnmlAii///57jcf27dun3HjjjYqLi4vi7Oys9O/fX5kzZ47p8eTkZOWWW25RdDqdEhoaqnz77bdKWFhYvQmiXEzO1Wg0ipeXl3LNNdcoH3zwQY3fy3fffad07NhR0el0yrBhw5QNGzYogHLo0CHTPm+99ZYSEBCgqFQqZeLEiYqiKMrWrVuVXr16KTqdTunfv78SERFRZ9KuNX9nKkWxIhunieXm5uLu7k5OTg5ubm7N3RwhhLhsTfG+VlxczNmzZ+nUqRMODg42OYcQ1vydyTCKEEIIIWxKgg0hhBBC2JQEG0IIIYSwKQk2hBBCCGFTEmwIIYQQwqYk2BBCCCGETUmwIYQQQgibkmBDCCGEEDYlwYYQQgghbEqCDSGEEELYlE2Djblz5zJkyBBcXV3x8/Pjjjvu4OTJk7Y8pRBCiFZq0qRJqFQqnnzyyRqPTZ8+HZVKxaRJk5q+YeKy2TTY2LFjB9OnT2fPnj1s3bqVkpISbrrpJgoKCmx5WiGEEK1USEgIK1eupKioyLStuLiY7777jtDQ0GZsmbgcNg02tmzZwqRJk+jTpw8DBgxg6dKlJCQkEBkZacvTCiGEaKUGDRpESEgIa9euNW1bu3YtoaGhXHHFFaZtRqORuXPn0qlTJxwdHRkwYABr1qwxPV5WVsbkyZNNj/fo0YOPPvqoyrkmTZrEHXfcwYIFCwgMDMTb25vp06dTUlJi+wttZ+ya8mQ5OTkAeHl5mX1cr9ej1+tNP+fm5jZJu4QQoi1TFCgsbJ5zOzmBSmXdMY899hhLlixhwoQJAHz99dc8+uijREREmPaZO3cuy5cv54svvqBbt27s3LmTf/zjH/j6+jJixAiMRiPBwcGsXr0ab29v/vzzT6ZMmUJgYCD33Xef6Xm2b99OYGAg27dvJzY2lvvvv5+BAwfyxBNPNMbli4uabIl5o9HIbbfdRnZ2Nrt27TK7zxtvvMGbb75ZY7ssMS+EaCuaY4n5ggJwcbHJqeqVnw/OzpbtO2nSJLKzs1m0aBEhISGmHL+ePXuSmJjI448/joeHB19++SVeXl789ttvDBs2zHT8448/TmFhId99953Z558xYwYpKSmmHpBJkyYRERHB6dOn0Wg0ANx3332o1WpWrlx5GVfdPlizxHyT9WxMnz6dY8eO1RpoALz00kvMmjXL9HNubi4hISFN0TwhhBAthK+vL7fccgtLly5FURRuueUWfHx8TI/HxsZSWFjIjTfeWOU4g8FQZahl4cKFfP311yQkJFBUVITBYGDgwIFVjunTp48p0AAIDAwkKirKNhfWjjVJsDFjxgw2btzIzp07CQ4OrnU/nU6HTqdriiYJIUS74eRU3sPQXOduiMcee4wZM2YA5UFDZfkXL2bTpk106NChymMVnyErV67k+eef5z//+Q/Dhg3D1dWV+fPns3fv3ir729vbV/lZpVJhNBob1mhRK5sGG4qi8NRTT7Fu3ToiIiLo1KmTLU8nhBDCDJXK8qGMlmLs2LEYDAZUKhVjxoyp8ljv3r3R6XQkJCQwYsQIs8fv3r2b4cOHM23aNNO206dP27TNonY2DTamT5/Od999x/r163F1dSUlJQUAd3d3HB0dbXlqIYQQrZhGoyE6Otr078pcXV15/vnnefbZZzEajVxzzTXk5OSwe/du3NzcmDhxIt26dePbb7/ll19+oVOnTixbtoz9+/fLTW8zsWmw8fnnnwNw/fXXV9m+ZMkSKcwihBCiTnUl0L799tv4+voyd+5czpw5g4eHB4MGDeL//u//AJg6dSqHDh3i/vvvR6VS8eCDDzJt2jR+/vnnpmq+qKTJZqM0RFNkbQshRFNqjtkoQtiCNX9nsjaKEEIIIWxKgg0hhBBC2JQEG0IIIYSwKQk2hBBCCGFTEmwIIYQQwqYk2BBCCCGETUmwIYQQQgibkmBDCCGEEDYlwYYQQgghbEqCDSGEEO2OoihMmTIFLy8vVCoVhw8f5vrrr+eZZ56p87iOHTvy4YcfNkkb25ImWWJeCCFE8/tga0yTnu/ZG7s36LiUlBTmzJnDpk2bSEpKws/Pj4EDB/LMM88watSoRmnbli1bWLp0KREREXTu3BkfHx/Wrl1bY8l50Tgk2BBCCNFixMXFcfXVV+Ph4cH8+fPp168fJSUl/PLLL0yfPp0TJ040ynlOnz5NYGAgw4cPN23z8vJqlOcWNckwihBCiBZj2rRpqFQq9u3bx91330337t3p06cPs2bNYs+ePQAkJCRw++234+LigpubG/fddx+pqamm53jjjTcYOHAgy5Yto2PHjri7u/PAAw+Ql5cHwKRJk3jqqadISEhApVLRsWNHgBrDKGlpaYwfPx5HR0c6derEihUrarQ3Ozubxx9/HF9fX9zc3Ljhhhs4cuSIxW0BMBqNvPfee3Tt2hWdTkdoaChz5swxPZ6YmMh9992Hh4cHXl5e3H777cTFxTXGy91kJNhoY0rKjOQUlZCaW0xcRgHRybkcSsgiMj6Lo+eyOX4+l1OpeZzNKCApu4i84pLmbrIQQgCQmZnJli1bmD59Os7OzjUe9/DwwGg0cvvtt5OZmcmOHTvYunUrZ86c4f7776+y7+nTp/npp5/YuHEjGzduZMeOHcybNw+Ajz76iLfeeovg4GCSk5PZv3+/2fZMmjSJxMREtm/fzpo1a/jss89IS0urss+9995LWloaP//8M5GRkQwaNIhRo0aRmZlpUVsAXnrpJebNm8err77K8ePH+e677/D39wegpKSEMWPG4Orqyh9//MHu3btxcXFh7NixGAyGhr3QzUCGUVqxAn0paXl60nKLy7/n6cktsj540Nqp8XCyx9NJi6eTFl9XLcGeTjjYa2zQaiGEMC82NhZFUejZs2et+2zbto2oqCjOnj1LSEgIAN9++y19+vRh//79DBkyBCjvLVi6dCmurq4APPzww2zbto05c+bg7u6Oq6srGo2GgIAAs+eJiYnh559/Zt++fabnXLx4Mb169TLts2vXLvbt20daWho6nQ6ABQsW8NNPP7FmzRqmTJlSb1vy8vL46KOP+PTTT5k4cSIAXbp04ZprrgFg1apVGI1GvvrqK1QqFQBLlizBw8ODiIgIbrrppga80k1Pgo1WxGhUSMou4nR6PmfSC8hpQGBhjqHUSFqunrRcvWmbSgV+rg6EeDkS4ulEB09H7DXSESaEsB1FUerdJzo6mpCQEFOgAdC7d288PDyIjo42BQYdO3Y0fbgDBAYG1uiVqO88dnZ2hIeHm7b17NkTDw8P089HjhwhPz8fb2/vKscWFRVx+vRp0891tSU6Ohq9Xl9r4uuRI0eIjY2tcjxAcXFxlXO0dBJstHBlRoW4CwWcTsvnTEYBRYayJjmvokBqbjGpucUciMtCo1YR5u1E70A3Ovk4YyeBhxCikXXr1g2VStUoSaDVZ5WoVCqMRuNlP29l+fn5BAYGEhERUeOxykFJXW1xdHSs9xzh4eFm80V8fX2tb3QzkWCjhSo0lHL0XA5R53LI15fWeDwxJor/LZrP+CdeIKR7P5u3p8yocCa9gDPpBejs1XTzc6VXoCsdPBxNXXtCCHE5vLy8GDNmDAsXLmTmzJk18jays7Pp1asXiYmJJCYmmno3jh8/TnZ2Nr179260tvTs2ZPS0lIiIyNNvSUnT54kOzvbtM+gQYNISUnBzs7OlGRqrW7duuHo6Mi2bdt4/PHHazw+aNAgVq1ahZ+fH25ubg06R0sgt6ctTFpuMb/8ncLiP87y1+kLZgMNgP1b1xN7ZC8HfltPYkwUn73wCIkxUU3SRn2JkWNJOaw+cI4lu+M4nJhNSVnj3jEIIdqnhQsXUlZWxpVXXsmPP/7IqVOniI6O5uOPP2bYsGGMHj2afv36MWHCBA4ePMi+fft45JFHGDFiBIMHD260dvTo0YOxY8cydepU9u7dS2RkJI8//niVnojRo0czbNgw7rjjDn799Vfi4uL4888/efnllzlw4IBF53FwcGD27Nn861//4ttvv+X06dPs2bOHxYsXAzBhwgR8fHy4/fbb+eOPPzh79iwRERHMnDmTc+fONdr12poEGy3E+ewiVh9IZMXeBI6fz6XUWHPsMjM1icSYY5w79TeHd2wC4FDEZrav/prYI3vZteG7pm42OUUlbD+RxuJdZ/nzdAaFBvPBkRBCWKJz584cPHiQkSNH8txzz9G3b19uvPFGtm3bxueff45KpWL9+vV4enpy3XXXMXr0aDp37syqVasavS1LliwhKCiIESNGcNdddzFlyhT8/PxMj6tUKjZv3sx1113Ho48+Svfu3XnggQeIj483zSaxxKuvvspzzz3Ha6+9Rq9evbj//vtNOR1OTk7s3LmT0NBQ7rrrLnr16sXkyZMpLi5uVT0dKsWSjJxmkpubi7u7Ozk5Oa3qRbVGdqGBXbEZnErNr3ffWTf1qH+fhWtRFAVnd0+8/Ds0RhPNqm0Yx06toneQG4PDvHB3kkp8QlTXFO9rxcXFnD17lk6dOuHg4GCTcwhhzd+Z5Gw0kyJDGXvOXiDqXA5lZnoxzBk7cSZbvvm4zn3en37XpX//evKy2liXysM4lYONUqPC0XM5/H0+l/7B7gzt5I2jVqbQCiFEeybBRhNTFIVDidnsOXMBfYl1eQ752VkW7afWaHjw+Xn171iJJQmnmalJFORklS9adHEY58BvG4g7fohrbn+YLv2HmHpTyowKhxKyOZ6cy+AwL64I9ZCps0II0U5JsNGEMgsMbD2ewvnsYsuPMfMBD5QXwqhlBOyZj1cT3K2PVW2rraeisncevqHGtqL8HBJjcvh+/mygZm+KvsTI7tgMjp7L5qrO3vQJcpPZK0II0c5IsNEEjEaFgwlZ/HX6gtnEz7qY+4AHzAYaKpUKRVFY/dFr3DPzjXqnxJoLZA5FbGbIjXeazfuYMHs+3y94EWOZ+Vof4yY+TWLMMbP5InnFpWw9nsqxpBxG9fLH11VXZ9uEEEK0HZIgamMZ+Xq2Hk8lJcfy3ozKIrdtqPUDXqVWo3N0xicolKvG3cfeLWtIiT9Fib6Ya+94mDunvVLnc1uScFq9p+Lcqb+r5IVYelxlapWKK0I9uKqzN1o7GVoR7UtTJoh27Nix3qJRQjRUUVERcXFxFiWIyju9DR1OzOb7vQkWBxrm6mWEj7qNZz5ebXb/Zz9Zw1ur/mTiqx8T0r0v98x8A62DE1DeQ3Hu1N8kxhwjMzXJ7PETZs9HrTGfvKnWaJgwe75F7bb2OKOiEBmfxbI98ZxJr38WjhDCOpqL/69b00JdovUpLCwEalZINcemwyg7d+5k/vz5REZGkpyczLp167jjjjtsecoWwVBqZFt0KidS8urfuZL68iYqhkkqvgPYabXMeaRmTf387Av1zkwJH3Ub/qFdzPZU1Jb34eLhjaunD06uHqQmxFp8nDm5RSWsP3ye7v6u3NDTT2atCNFI7OzscHJyIj09HXt7e9Rqua8UjUdRFAoLC0lLS8PDw8MU3NbFpsFGQUEBAwYM4LHHHuOuu+rvem8LLuTr2Xg0mcwCy+4oLMmbqPiA9/ANZOjYe9i7ZQ3Z6cm4eJQv/lNXLkXFzJT6ZpuYC2TM8fAN4NVl20mOi+GDGXdbfFxdYlLzSMouZHQvfzr7ujToOYQQl6hUKgIDAzl79izx8fHN3RzRRnl4eNS6am51Ng02xo0bx7hx42x5ihYlOjmX30+kYSi1fEqruQRQc70Sry7bjsbeHpVKxbBb7qespAQ7rRYw10OhAbwBH+6e8QV29l3YvHQdsUf82PpdLCPvHYK9TkHrYETn5FtnIGMuSLHTanH19KnzOGsV6MtYf/g8fTu4M6K7r+RyCHGZtFot3bp1k6EUYRP29vYW9WhUaLIEUZVKVe8wil6vR6+/tMx5bm4uISEhLT5B1GhUiIhJ40hijtXHfvHiY8Qc3G32sYpeifBRt5l9vKwUUhO0nDvlwIkDBRzekQr0BLwsPr9KpeDhW4JvcAl+IQZ8OhhAOc7RXW9w+9SZ7N+6nl3rl5lNOC01GEwBkKIoVQKgy+HuaM9NffwJ9nS67OcSoqVpC4nvQlirRU19nTt3Lm+++WZzN8Mq+tIyNkclE5dRWOd+lXsInN29TEMn589E13qMufyH1AQtUbtdOL7XmXOxOkoNlXsAulT6txHIBDIufpUAjoDTxe+OgBuK4kZWmpasNC0xBytWWPQHrmbRyycx6HsBV3Nw+9YaU2IrBxYqlYrkuJONshJtTlEJayLPcWVHL67q7I1aLXU5hBCiNWtRwcZLL73ErFmzTD9X9Gy0VLnF5QmOGXn6evetnPz5x0/LLD6H0QiJJx2I+tOFY7tdSDtXtedA51RGcFc9wd2KCepUQIdupbh5leHoXEbktnX88OErteZyPPDcPHoOuYP0c1rOHi8gJU5NdporsUfKgBDyc/oB/YCXKMjJ5/3pW4EVwEbe/7XmCrOWFAazlKLA3rOZnMsuYlzfAFwdZJ0VIYRorVpUsKHT6dDpWkexp9TcYjYcPl/rEvBQe/LnuIlP88uyTzEazdTOUKnw9O+AoUhL1J/9Wfp2CJkplwIMjb2RbgOL6Hd1Pl0HFOIdWIL5RHMVgZ26EdS5F+dOHavx6KVeEyMu7sV88mzvant0AUYBN1z88gXuBO5Eqyvmhw+KCR+Vh4dvLIV5Na8xpHs/dq77lpv+MZ2+w2opTGZG9RyRpKwiVuxNYEyfADr5ONf/BEIIIVqcFpWzUV1LHduMTctny7FkSsouvXTmEiktKZpV3aOv/UrMocHs3+qGobg8+UbnaKTXlfn0G15ArysLcHC2LAF17cJ32LW+vBel+qyR8tVhjaY2pyWeraM6qAoYADyAi8dM8rMrFwmKAz4HvgTM56xYsyBcRZur54ioVBAe5snVXXxkWEW0ai31fU0IW7Jpz0Z+fj6xsZdqMZw9e5bDhw/j5eVFaGioLU9tM0fPZfP7ibQa1cLNDSHUNSVVpVajGI2VpowOB15myVs3mvYJCNNz7Z3ZhN+Qi9bBspjQXG+KSqXCL6QLfYeP4vjeCPKzL+Di4c3vP3xlavOd016pteYGKKhUR1CUwzz+dnf0RYOJ3ObGkT9c0Bd2BN4FXgb+C3wEnDMdqXNy4dypv+tc9t7SsukH4rI4n13Ezf0CZVhFCCFaEZv2bERERDBy5Mga2ydOnMjSpUvrPb6l3QHsj8tk16kM08+VPyT/+/Lj5Gdn4uLhzW1TZpuGEDx8/M1+gD/xzn9Z+Z//w9l9KBrNfJJO9yp/QKXQZ2gB196ZRbeBRVi7ZpklvSkzP1yJnb22SpunzFlESsJpvnv3hRo1M0Y/MJWTB//kQnICfiGduXPay4R074dBr+JwhCtbVzhxIaXi91MCfA/MB2oO35jr5bC2bLqzTsPN/QJltopolVra+5oQTUHWRrHQrlMZ7I/LrLLNkg/JWQvX8v70u2oMYzw5byNHd13Fnp+9MJapUKkVBo/OYvSDOfh2KGlwO+taS0Wt0dS6iFplId37MXTsPezZ/APZGSnMWrgWdx9/1n76Frv/912NIY7Ek3/zwVOLgBeAysHlcuAVIL7Oabz1tdnccWqVimu7+zAo1LPe6xGiJWlJ72tCNBUJNuqhKAoRJ9M5nJhd47G6PiShfAjhkf97nxXv/Qt3b3+uHv8Qe37+ibTEu4AXMRSXj2L1HprPrZMzCOh4qfhOfRU/61LbYmmzFq4lNeF0nR/s9z3zDkNuutMUFGUkxVNcmF+j92bKnEWmIQ61WsMHM+7GwzeQgI4PsP/XrsB9F59VDyzkn+/2oNsV3Wq9rrraXFf5854Brozu7Y+9RoqAidahJbyvCdHUWtRslJbGaFT49Xgq0cm5Zh+va20RAH1hPotemQJAYW42Id0nsmvDbAzF5TNugrsWM35KOt0GFtU4tnIOCNCgwMNcKXFr1kNRqVTMfWxMjf3qqnCaFHuc/b/eRXkex7vAaGAWX79hYPSDOWSmbbZ6/Ze6nEjJI6PAwPj+gXg4XX5BMSGEEI1Pgo1alBkVNkclE5t2+auSqtRO9Lrydz6aGYrRqMLFvZTbpqYz6Ia8KtNWa0uUzMu6QOyRveza8B0PPj+33vPVt5aKqV0WfLBbsu4KYCrwdencJVw55k92rF1DZvJ09EX92PS1Lyr1a0Cy1eu/1CUjT8/3+xK5tX8gIV6SxyGEEC2NDKOYUWZU2BSVzGkLAo3s9BTTEEL3QcPZtvLLantchZf/r2SmugJwxchc7pyWhot7zemrluaA1DWzo0JdpcQrt7nyB/uzn/6Ih2/NRXWsHeKofu7nxvQCJgD/BiqKtP0APAMkA+W9I5db/lytUnFddx+ukDwO0YLJMIpojyTYqMZ4MdCwpkej4kMyKfZ4pWRQDfAO5UmTaly9Srl3Zip9hxfU+jz15YBUZ039igoVORPjJj1Dx95XWPTBXhFsmKvVYcly8peuywF4g/Igw47yuhwv031QFE/OW2T1tdSmXwd3Rvb0QyP1OEQLJMGGaI8kq64So1Fh8zHrh07stFpUKpVpKCCw00i8g04DswE1/a9N5V//jasz0IDyfIpnPl5d7/nUGg0TZs+3qo0VKnJBDu/YjOrivFqVSlVnD0LFdQV368s9M98kuFtfXD19LF7h1S+kE0GdewEFlAdfg4G9gDvwKbFHPuDwjgQSY46RmZrUoOuqLCophx8PnqPQUHt1VyGEEE1HejYuMhoVfj6WQkxq3mU9z4n99nw3P4T8bDu0DiV4+r/CQy9caXFiZ/VeBHMs7VGoUFs9kMozSuoakoHLW+G1ciXTS9TAVGAu5UGHnvLCYB/w/q+1L05nDTdHe24bEISva+sogS/aB+nZEO2RJIhSHmhs+fvyAg2jEX5f5cXP33ijGFUEdS4msNMrRG77Dwd+e7hKsFHXtNbKiZLVc0AsnaFR3TsP11ybxNyMkrpUX+G1vkCjtkqmrl5+5GWmoShGysucr6O88uh4YAG+wU9zIbkU78CSKq8TWD8jJ7eohB8OJDKubwCdfV0sOkYIIUTja/fBhqKUT289mdLwQKMwT82KdwOI3lf+gdZnWBw33HeEJW8uAWqW3q5rdVQP3wDTNNKcjFT2/fIj7j7+XDXuPqtmaFRm6YySxmQuwFEUhdwLqdW2pgC3AY9hZ/856edCWPCkkdufTCPp9KXXSVFo0IqyhlIjG46c59puvoSHSeKoEEI0h3Y/jPL7iVSOJJpfQMwSF5LtWfRyB9LOabHTGik1PA4sqfMYFw8vi4cyLneGhtZOjYO9hsSYY7w5uWb1TmuHZCxVV7JrzXVhyvUe+hh52f8h8aQHAGq7LRhL/4GjS3nuRVF+Dk5uHjw592uLh38qk8RR0RLIMIpoj9p1sGGuBLk14qId+Pq1IPJz7PDwLWHyW+dJiVtt1YySyhoyu0SlAg9He/zcHPB30+HlrMNZq8FBq8HJXoPdxcqaBw8eJDw8HLVajdFoNH3f9sdfhHbvS3ahgazCErIKDVzIN1BcYn37q6ttyuwT7/yX796bjYuHN32Hj2L76q8wlpXh7OZBQW4u5bNV/g3oKF/U7UFgl9lzPPvpGquGV0K8nLi1fyAO9pqGX5gQl0GCDdEetdthlH1nMxsUaFTkEfQc/AFblg2l1KAmuGsxk99Owt27jA5daq/QqVZrMBovbyhDpQJ/Nwe6+LoQ6O6Ar6vOog9OPz8/AgICCAkJYfLkySxevJjExES6dwwm2K9qPoOiKKTn6UnMKuJcViHnsoowlFq2rL35NledMuvq6UNBbhYFuVmkJlxaFbggN/viv94HtlJei6MnEAG8Rnky6aXYeOwjM62utJqYWcjKfQncPrADns5ScVQIIZpCuww2Didmszs2o/4dzdj363pijwwn9shQQE3vofk8/H/J6BxrdhBV/5B98IV5rHj3hRr7VS8TXp1apaKDpyNd/Vzo4uvcoOXVg4ODiYuLQ3txmu6UKVMwGAzodDVnaqhUKvzcHPBzcyA8zBOjUSE1r5jYtHxOpuSRV2zZlNK6qoLWlUdSPswSRfkU2c+AR4A5wAjgYSANgNSEWE4d3gNYV2k1q7CElful4qgQQjSVdjeM8vf5HLYeT8Waq66YWaEoKhY+d4ESw+MADLjuDNfffQRXL48quQO1Veh84Lm5LHrlCYuLY7k72jMw1INeAW44altGt7+iKCRmFhGdkktsWn69PR515ZzUNswyYfZ8VlRZ6n4SsBBworzi6IPAjjrPa0mlVY1axahefvQJcq/vsoVoNDKMItqjdtWzEZuWx2/H06wKNKBiZoUG+AZ4HDACz3Nk5wcc2Vm+T+V8i8ozSlQqFcNuuZ+ykhLyczItWv8jxMuJK0I96OzjbCq81VKoVCpCvZ0I9Xbihp5GTqfnczA+m9TcYrP7WzJltnrw5eTqbnqd+l19I1u+/Rhj2V7Kh1X6Ar8B/wI+qLWdlkzrLTMq/Pp3KtmFJQzv4t3iXmshhGgr2k3PRmJmIT8dSqLUaP3l7vt1IysXdATuAkqAf1D+wXcp3yJ8VM2ZHubUdqevVqnoGejKoFDPVlmEKjGzkMj4LOIuFFgczNW1RouLu5fpdSrR6zkUsYlVH8xBMS6kfFgFYCUwGSg0+/zW/G66+7tyUx9Zql7YnvRsiPaoXQQbabnFrI48V2uXf13Fowx6FV/MdiXueADlVS7vATaajg3p3pd7Zr5hVe2H6jr7OnNNVx+8XVpfkFFdRr6eg/FZnEjJo8yCwM6aqb2Xhl2mAR8C9kAU5UFgbI39rZ3WG+DuwG0DgnDWtasOP9HEJNgQ7VGbv43LKjCw7lBSnbkFlWc0VP63vkjFV690uBhoFAK3olJdqoYJkBhzzDQTwlr+bg7cEx7M7QM7tIlAA8DHRcdNfQKYOKwjPQNcqW9komJdGbCsMmn5fp8D11Oev9EPjd1h4JZKjzdsOCQlp5i5327kmhHXc+DAgQY9hxBCiJra9C1cXnEJaw8lUWioOeOhcjntg9s3AHDgtw2mx/ds2cqezU9TYugG5AG3oFLtwjswlE59wzl7LJILyQkoilKjQmh9haZcHey4ppsPPfxd22yegLuTPeP6BTIozJNdpzJIyDQ/1GGNqrNbxrB7w4OkJs6jrPQqYCNu3h9z40PZ7PulYZVWAXZsWsvunTv49L9fs3Tw4MtusxBCiDY8jFJcUsbqA4lk5BvMPj7rph51HO1KeRLilUAWMAbYb/G5a0tIVKmgf7A7V3f1QWfXMmaXNJW4jAJ2xWaQnqe/rOepPuyiLyxl05IO7N7gAcDAEXncPysZjcZQo5ektjVpaluo7ovlq+kZ4IaPjw9hYWGX1W4hKsgwimiP2mzPRsTJ9FoDDahrvRAnYBPlgUYGMBo4YtE56yrO5e5oz429/dttXYeOPs6EeTtx5FwOu2MzGlwkrPrsFgdne+6ekUZQZz1rP/Xj8A5Xzp8uxdF1IndNn1glqKhtTZraFqr7x62XtrfgmFw0E6NRwagopiq9Qojatdn/JYayuj/MwkfdxjMfr662VQf8BFwLZAM3YS7QmDB7vtnnfObj1TVmPqhUMDDUg39cFdZuA40KKpWKgSEeTBzekW7+jbsK67Cbc/jnu+dwcS8l7Zwn8dGL+H3V32SmJpEYc4xzp/42rT57KGIz5079TWLMMTJTk5gwez5qjfmeJrVGw1NvfUSRmaE40X6VlBn539HzFOjl70IIS7TJno0DBw7wyrSnGT1xlkWzRMq75O2A1cCNQD4wDjhUfU8ql8uuXh+iOlcHO8b2DSDYs30HGdW56Oy4tX8QZzMK+P1EGrlFJZf9nJmpSdjrsnjgOSe+fiMYo7EvR/6YzZE/ngSWVtk3P/tCjToctZWYr6ju+v2+BG4bGIRPG0nkFQ1XaChl/eHzpOQUc3335m6NEK1Dmww2vv32W6L2/4lHhy5Vgo3qY/YVCYfuPsGUliwhJW4gUASMB/ZUec7RD0zl5ME/yU5PxrdDx3qLc3X0cWJsn8AWU/mzJerk48wjw8L46/QFDiZkWV1srbKqQyHOlAcY91C+Am8v4CXKi7FdYm7Yq7YAMqeohFX7ExnbN4Auvo3bKyNaj+zC8tlt2YWXHyAL0Z60mQTR+Ph4MjIyUKlUjBs3jrS0tBpLuEesWcKu9cu49o6HuXPaKwAYig2s/iiEyG3uaOyM3PLYH2z47/U1PmxmLVxLh669TXUgaqsPoVLBVZ29GdrJq83ONLGFpOwifjmWQk4DezlqLmmvAl6/+AWwHpgAFJiOqajDkRgTxbqF75CWFId3QEiNAmMevgGmY1QquLqrD0M6ejWonaL1SskpZv3hqrPbHru6E+5O1q1VJAmioj1qMz0bHTt2NP274kO+ene5i0f5B0TFVFWjUWH3/4YQuc0dtVrh4f9LIbSHiu2rzfdaVK4DYa4Mt6NWw7i+AYR5OzfBFbctHTwc+cdVYeyMSScqKcfq48NHVV9tVwHeAE4CXwO3U75M/XjKl62/ZP/W9cRFH+bq8Q9x14zXqpSYrz6jRVFg16kMLuTrGdVLKo62F2fS89kclUxJWYu9NxOiRWszwcby5cuZNGkSpaWltc4cyM/OvPi9Igh5gfLqk/DgCyn0vyYfML+uSX3FpgLdHbilf2CDVmQV5bR2akb39qezrzO/Rac2UvLd98BZyhN/BwL7cPOaQlnpbnKzMqokjh754xeGjr3Honop0cl5XCgwcGv/INwd5Xfelu2Py2R3bMZlDfMJ0d41yTDKwoULmT9/PikpKQwYMIBPPvmEK6+8st7jrO1uPHjwIOHh4TW2q9UajMbqH1z/AJYBcNuUNK6/J9uCKzGvu78rY/r4yxS4RlRkKOPX4ymcSS+of+eLKq+10n3QcLat/LLSo6HA/4D+2GmNlBoeoDwhuHa11UupzFGr4ZZ+slR9W1RaZuS36FSik/Nq3UeGUYSwjM0/HVetWsWsWbN4/fXXOXjwIAMGDGDMmDGkpaXZ7JxqdfllVQynPPhC9doXN1LetQ4j7s68rEAjPMyTm/sFSKDRyBy1Gm4f2IERPXzRqC3LfalYbfeZT1Yz4NoxAJVKoScCV9OpbwqlBjXlC+n9n9nnUWs0tU5vrq7IUMbag0lExmdZtL9oHfL1payOPFdnoCGEsJzNPyHff/99nnjiCR599FF69+7NF198gZOTE19//XWjn8vPz4+AgADCw8P55yvzCO7WF1dPH5xcPYCKD54rgB8Be3oMTmT8ExkNOpdapWJkTz+u6+4riaA2NCjUk/sGh1g8VFGx1krFTKPgbn25Z+abF/8WHJgwO4YRd1UEBnMon7VSdYjMXL2UuhgVhZ0x6fwclUxJPfVdRMuXklPM93sTSMkpbu6mCNFm2DRnw2AwEBkZyUsvvWTaplarGT16NH/99Vejny84OJi4uDi0Wi3/O5pM12vvoKykhPycTFw9fXBxH0xm6mr0RS5o7HZy9wwVarW/1eex16gY1y9QpkA2kQB3Bx4aGspv0amcSs236JiKXo7KuTcZSfEU5GQQPiqCPT/vRl80D5gIdALuBFUWlzMwfyIlj4wCA7f0C8TLuf4F5UTL8/f5HH6PTqPUghWLhRCWs2mwkZGRQVlZGf7+VT/Q/f39OXHiRI399Xo9ev2ltTNyc3OtPqdOd6noUsUsEQ/fAGYt3MnCF7qgL9IS1LmYJ+f54OJh/eU7XezeD3B3sPpY0XAO9hpu7R/E4cRsdsakW7R8ffUZQ3MfG1NtjyjKh1OuA/aAciuunpkNWsCtQkaenu/3JTC6lz89Alwb/DyiaZWUGfn9RBrHz1v+npMYE8X4d6by/n/mM1gW7ROiTi0q0WDu3Lm4u7ubvkJCQhrleUtLYPm8MDKStHj6lzBlTlKDAg1nnYZ7woMl0GhGA0M8uCc8GBed9b+/mmXJfwWGUz5bpRtah8M8+MK+KnU1GsJQamRzVDK/n0ilVIZVWrzMAgMr9yVYFWhA+ZTpP3ZGsGzZMhu1TIi2w6bBho+PDxqNhtTU1CrbU1NTCQio+Yb+0ksvkZOTY/pKTEy87DYoCqz52J/TR53QOZXx+FtJuHlbP6XS1cGOe8ND8JZy1c0uyMORh4aG0sHT0arjzK+HcxwYSmCnCxiKHfnq1Y7s+blxZggcSczhhwPnyGlAtckDBw5www03cODAgUZpizAvOjmX7/cl1LloY2Xm1tpZuXIlBw8eJDIykvj4eFs2V4hWy6bDKFqtlvDwcLZt28Ydd9wBgNFoZNu2bcyYMaPG/jqdrsowSGPYvtqTfb+4o1IrPPJyMoGdLHtTqczVwY57woPxcJJx+JbCWWfHPYOC2XEqncMJ2VYfX7UseTr3Pr2LP9aP5NB2N374IIC0RC23Ts5AfZnV5lNzi1mxL55RPa0bVvn222/Zvn07y5Ytky56GygpM7LjpPUF5MytEJyenl5lyn0LLsosRLOx+TDKrFmzWLRoEd988w3R0dH885//pKCggEcffdTWp+boLhc2LfYB4I5/ptNrSKHVz+HuaM+9g0Mk0GiB1GoVI3v4MbZvAPYay2YEmZ+l4oOHryf/eDGFMY+Uz06KWOPFkjeD0Bdd3kyjxJgoPnjmH3z5469sOZaMvrT2XrX4+HgiIyM5ePAgq1atAuSu2RaSc4pYsSe+QZVqza0QXBFc2NnZsXz58kZpoxBtTZMU9fr0009NRb0GDhzIxx9/zNChQ+s97nKK33ywMp0XJ3lj0Ku5+rZs7p5hfV0PDyd77g4Pxk2qgrZ4abnFbDhynrzi0nr3rW1dmwqHtrvy/QJ/SkvUBHUuZvJb5/H0q/95zVm78J0q6/G4Odoztm8AHTxqDgFVnkJd24JwctfccEajwp6zF9h/NgvjZbyO5079bXaF4MjISAYNGlTv8VLUS7RHTZIgOmPGDOLj49Hr9ezdu9eiQONyJCXBOzO9MOjV9BhcwB3/tD7QcHe05x4JNFoNP7fy6bHmPsSrq6jFAVRZ76bCFSPzmLbgHK6epZw/48CHT4USF215UrC5cf1DEZs5d+pv/j5yiP9u2sufsRkYq82oWb58OXZ25SObFUGF3DU3jswCAyv3J7L3TGaNQCMxJorPXniExJgoK59VDQSbiggKIWrXZlZ9raDXw/DhcPAg+IfpmflhIo7O1s0IcNZpuE+GTlqlMqPCtuhU/rZyZoE5WWl2fPVaB5LP6NDYG7n/2VQGj66/ouSsm3rUu8/7v57E382B0b398HO9FMjUVnLf0rtmUZXRqHAoMZu/TmdUWUQtMSaK/y2az/gnXmD/1vU1VoOuTcUKwannMlGM3+KgHUTHsHtISTnK/v37CQ4OrrdN0rMh2qM2F5JrtTBxInj6lM88sTbQ0NmrufMKSQZtrTRqFTf1CWBED1/Ul1nZ1dOvlJkfJNB3eD5lJWq+ey+Q/33lQ41ldqh6d2xuXL9C5VLoqbnFfL83kd2xGTWmyFbcLctdc8Ol5Razcn8iO2PSa6zW+sf6FcQe2UvEmq9r9D4lxhwjMzXJ7HOWrxBcgmLcS3HBCIoKXfj3v7cQFxdnUaAhRHvVZlZ9raBSwcyZ4Dc4jeSC2sfZK9/ZhHTvB5RXBr1jYAd8XWV6a2s3KNQTb2ctm6KS0Zc0vNaFzlFh9AO/kBSbS1bao2z/wYvUeC3/eDEFh0qB7P6t64k9spcDv63nzmmvVFvu/pJnPl5NcLc+pp+NisK+s5nEpuUzure/qeR+SEgIkydPZvHixSQmJuLn59fga2hvDKVG/jpzgcMJ2VWGTDJTkyjIyUKlUnFg6zqgPMCocGk16HIVC/FVPm7/VntgL8UFzrh5FbBg/hl693ZDpwtrmosTopVqc8FGBSdnBepYMLTyh0NI935o1Cpu7R9EkAVj/qJ1CPN25oEhoaw/nER2A2pdVDiwbT1ZacvoMbiI00ef5PheFz5+JoS7njqIziEZlUpV5e54yI13kpJwGqg90bO6zAIDqw8kMiDYg5Oxp3F1ckSlUjFlyhQMBkOjTwlvq86k5/P7iTSzicLmpq2ao9ZoePD5eaYbktgjewENMA9YfHGvreRmPsiUyRcASdwVoj5tNtgwp/IdSpUPh5vu5Oou3qgKtODj3MytFI3Jy1nLg1eG8r8j5zmXVWTxceb+VpJiX+O+Z3qx/ourSIl35LPnewBvUF6JtFz1u+Pgbn0ZOvYe9m5ZQ3Z6cp2l0BUFDidmcyotj+FdfOgT5IZKpZJAwwIX8vXsis3gTHrtdxgTZs/n+wUvYiyru6hfRe/T2oXvEHtkL13638Lpo88BIy/uMQ94GTBiZ2fH0qVLG+kqhGi72lyCaIUNR85zOi2/ynDJBzPuqfe4FvxyiMtgbeJo3UmegcBa4CqgjPKl6t+rsddNE6Yx5pGZNabYmhvCM8ffzYHre/hKb1sdCg2l/HX6AseSci2azlrbtFW41As1YfZ8/EO78N+XHyc/uxcq1WoUJRDIAyZR/rsvF7FrLyOuvtKqNkuCqGiP2nz2WeXhkroS92RqYdtWkTh6bTcfLMkbrTvJM40HnvuLoeNyKO9efxdYBVTtFSsqyDM7xbby32RdUnOLWbU/kZ+jkskrbvgwUFtUWmZk39lMluyO4+i5HKvrZqiq/RGMfmAqwd36ArDi3Rd4f/pd5Gc/CkRcDDSigaFUBBrVjxdC1K1NBhvx8fHEHj9ao86Bf2gXHnh+ntlj9u7dy4QJE5qymaIZDO7oxfgBQWjt6v7TN7+OSrkOXXoR2CmU+59NZfSDhwADcB/wF9DFtF/l2Q1n/o6stfZGXbMfoHzp+m/+jGNHTDoF+oYVF2srSsuMHE7MZumfceyOzcBQeilJ15J6GTUqyHbtg4uHN8PHP8Qzn6zmgefmolJ7AeuABZSPNK8AhgDReAUEV6k86+vna9sLFqKNaJPDKJbcdVR0marVaoxGo9QxaGcsqTha0eVuLsmzoiZDdnoKC56cR2HeYsqHV7KBicAGq9tUMfuhLvYaFf2CPRgc5olzA1a+ba0MpUaikrI5GJ9Nfi0BV/VqrbWpq4JsYoyOr1/3IeeCM6AHnga+BOD+Wf/myjF3VTluysgeuDtZV/hPhlFEe9QmezaWL1+ORmP+jVilVuPo7MrAQYP44osvCA8PJyAgQKYWtjN+bg7cPyQEP7faky8r7oIDOnZn1ANT8QvpYgpkK3ol8rIymPnR09wx7SfKezY8gPWU53CU/w2qNRqG3/qgRbU36lNSpnAwPoslu8+2i56O4pIy9p65wNe7z7IzJqNGoFFXtdbEmGMc++v3Gr0d5irIKgrsXOfBx8+GXAw0zgBXUxFoAMQc3F1n5VkhRO3aZM8GwAcrf2HWg2NrbP/X5+uYdtdIQn3dTHcoMrWw/TKUGtnydwqn0/LNPl5qMPCvW2tP4qzw/q8niY8+zkdPnwFmXdy6G7ifWQs/Ibhbn1qTE2ctXEtwtz4WJ45WplGr6ObnwoAQjzaVSJqWW8zRczmcTM2rMlRSnSXVWoE6ezvysjR8vyCAE/vLc266D8rg3KkhuHra0Xf4KLav/gpjWRnObh5Mnfs1iqLg7O6Jl38HHru6k/RsCGGBNtmzUVnlOxGAYV28CfNzr7JdAo32S2unZnz/QAaFeVbZXjH+nxx30uKKoBo7BXgOuIvy4ZSrgUPEHa/aa1b9b3L1R6+RGBNlceJoZWVGhRMpeazan8iKvfEcS8qhpKzhRcyaU0mZkWNJOXy/L4EVexOISsqpM9CAuhN5VWo1OicXoPb8mJMHnFjwZBgn9jtjZ2/krhmpTJ2bSWFeHKkJsWxb+aVpqmxBbjbvT7+LD2bcbXHNDiFEuTY76Ovu5X1x6fBAU52DgsxUBvfq1NxNEy2MSqViRHdfvJy0/H4iDaOiNKgiaMWwi4dvLH2uWsa2VaMp0fdi7Wc+5FzI5KpxPjX+JpPjYkiMOUbEmq85dXgPcKkwWOU7aEuk5er5+qff2LR4ATNffINbR11DiKcTanXLnTlRZlRIyCzkVGoesen5Vld7DR91W62/G8VoRF9Y3mNVvf7Je5tOsnmJDxFrvAAICNPzj/9LJqiTAai7JkdF0S8hhOXabLDh4x/Eq8u2mxLBHpz4KDf39sXBwfLVO0X74laWTVdVGrvOZDSoIqiHb0CVv7nr74X1n1/gr83ebFvpTcyhITw5bzdax3MU5mYT0r0vH8y423SeCrWVzbbE/q3riTm0h1Xfr8Dg2QkHew1dfJ3p5u9KqJcTmhYQeJSUGYm/UEBsWj5nMgouq5x8YkwUqz96A8Ds78ScsF4T+WhmKEmny98Lrh6fzfgp6Wh1l46tK4ipXnJeCFG/NhtsAKYELjdHe27uF4SDvfnuViEAOnbsWGObtRVBKycNanVw7zMX6D5Izw8f+ZN40oGPZnbFoJ8PfFVveyy9g661Mu7F3pHz7p787d8Bnb2aIHdHgjwc6eDpiL+rDjuN7UdS9aVlJGcXcz67iKTsIlJzi2ssjNZQ+7eu59ypY9jrHAgI62b63VxITqQwL7va3hpgNvHRrwNaHJz1jH8ihmE31z2MamnJeSFE7dp0sAFgp1Yxvn+gBBqiXsuXL2fSpEmUltac4aHWaLjvmXcYctOdqFQqht1yf5Upk3UZcF0+Yb2K+X5+AKcOOwGLgFuBycCFWo+z9A7aXP6Aud4RfYmRsxkFnM0oL+ltp1bh7+aAv7sDHo72eDjZ4+Goxc3RrkFFq8qMCjlFJWQVGsguLCG70EBKbjEZeQaLi25ZkiRrLrjSOjhxz8w3AOgx+BryMjP4+Jn7Kx3VE/gGqKj2+RPFBU+y+sNUht1svufo0rBYoMUl54UQ5rX5YGNkTz/83GToRNRvwoQJ9OrVi/Dw8BqPVf/gt3bqo4dvKVPnnWPHj55sXuJDWentlH/wTQX+V2Vfa++gG5pfUGpUSLrY21CZRq3C3dEeR3sN9nYq7DVq7NRqtHYq7NRqSo1GDKVGDGUKJaVGDGVGCg1l5BWXcLk3/tUXSDTHXHBVkJNpGpICeG3FDlw9fbDTOpGVdh8obwMOQBbwFLDi4mtT+5Tj6sNi1gSYQoiq2nSw0a+DO307uDd3M0QrVFHsreJ7Y5SnVqth5L1ZdB9UyJI3PclMCaS8+NdKYCaQzugHpnLy4J9W3UE3dn5BmVEhs8Bg1TENUdGLMeLuSbh5+dU6DOTs7klBTqapx8OS4MrDN4BHX/+LtQuDyUqtuNnYDDwBnAcse20qBxZSW0OIhmuzwUaAmwPh1aYzClEfPz8/AgICCAkJYfLkySxevJjExETuuro3hy6o652KaYkOXfRMnXuU/zx5EIN+GvAAavVY7B1eYtitDzLu0WcbfAfdHPkFDakPApd6McqXcK+q+jDQNbc/bPHsIJ+gfnz7bzgc0QnQoHM0oC+aBnwNSO6FEM2hzdbZGNLRs0Vk3ovWJTg4mLi4OPbu3cvUqVPZu3cvcXFxDO/fnQeGhOBhZQGnCtXX7fDt4Mc7P97AM5+eI6hzMUajB/rCz1n9YThZafZWBxo11vy4uHZHU+QXWFMfxFzFT52TC2p17bUyxk18ukZ10Mqzgyp/jzkUxLzHO3I4oiugwTf4L6YviMTVcz0h3Rv3tXF1sKt3jR0hRLk2W0FUCFsoLilj09FkEjILrTqurnU7ykph+2pPfl3uTWmJGq3OyKgHMhlxT1aV6ZiVmetNqGvNj8ZWOUmzfCn2TFw8vJkyZ1Gd9UEsrfhpqZDu/Rg69h52rY8kPWkGZaU3AqBSn0ExPomLx0GmzFlEicGAm7cv3gHBjfLa9Ahw5Yaefg1KPJf3NdEeSbAhhJWMRoUdMekcTsyucz9rP5DTEu354QN/zhxzAsDTv4Tbnkin/7X5VE8ZsXTRMVuxJGgwVx8kctuGWvMtyqmoPNShVmswGs3nZtz3zDv0HHIPvy73Yc9md4xGFeUr8M4D5gLFFrXJGg72Gkb29KVnQMPfj+R9TbRHEmwI0UDHknLYfiKNUqP5/0IN+UBWFDgU4crGRT5kZ5QP2XTpX8gd/0zD0eVsg3oTqmtojkXl49ISz9abpBk+6jazz1PbGjFObh54B4RUmWb6wHNzWfTKEzX2feqD9cQeGc7vq7zQF5UPZQR3O0XS6dtQjCesbpMlQr2cuKmPP64ODRtKqyDva6I9arMJokLYWt8O7vi46Nh41PxS9Q2ZkqpSwaCRefS5Kp/tq734/QdPTh914v3pYSjG7cA7lK9IWq4h1UYtmV5a33GWlnCvS/Vk1sff+pKwXgOqTDNNiT9VZV/QAo+w5M3R5GeX9wCFdC/mtinpdOmvcO7Uvxu96qedWsXV3Xy4IsSjUWYlCdEeSbAhxGUIcHfgoaGhbDyaTFJW1XoVlzMlVeeoMPaRC1w5JoeNX/lyeIcr8CjwMOXFqeYAZ03711dttL4qo5V7RSr3YDi7e9V6nCUl3M2prViWh29AjSXcK/Z18+6Ib4e3iNp9BWWlAeRng6dfCTc/lsEV1+ehrpan2VizcnxddYzrG4C3iyzWKMTlkGBDiMvkpLXjnkHB7DiVzuGEbLP7NPTDz8u/lEdeTmbEXVn8ssybEwecKa88+giXgo64eoMXS6uMQtUejD9+WlbvcfWVcK/OmmJZji6BXHfnUXau8yUptvztytWrhBvuzWL4+BzstVVfy8ut+lk50LrzpusY3sVHZrUJ0Qgk2BCiEajVKkb28CPAzYHfT6SZ6nE0RsnrxJgoNn1d/gE48HpvVi4oAsYBjwOTgJ9IiPGhQ1dqJJJWqG9IZ/wTs0mMOVajB2PcxKf5ZdmndSZpNqSEe33FstKT7Nmz2Z29W9wpzCuf8eHpV8IN92dy5ZjcGkFGhcut+lkRaGUf3ca10++u/wAhhEUkQVSIBjpw4AD/+te/eO+99xg8eLBp+4V8PZujksnIL6/CeblTUivPPBl57+N8MONuHF3GoFK9SWpCF9N+/mF6rhmfTfjoXBycav63ri0xc9bCtWa312fWwrWNuvppWSn8vceFPze6E3PQ2bTdJ8jA6AczCR+Vi8YGt0eVh5gWvzqFnMwM/Pz8+Pnnn1EUBR8fH8LCwhrtfPK+Jtojm/VszJkzh02bNnH48GG0Wi3Z2dm2OpUQzeLbb79l+/btLFu2rEqw4e2i44ErQ9kWnUZ0cm6DSl7XlWMx6bVPL9aMKOP8mbPsWu/Gwe2epMbr+PFTfzZ+7cMV1+cR1DmKIztf5LYps6okgpob0uk+6GpiDu6us03Vj1v90WvcM/MNq5JMq1MUOHdKx5E/XDmw1Y3cTLuL51Lo2DuF4oL/cN+sYYT1bPg56lN5iKkiZyQ9Pb3KGjkt+J5MiFbBZsGGwWDg3nvvZdiwYSxevNhWpxGiScXHx5ORkYFKpWLVqlUArFy5kokTJ1a5C7bXqBnbN4BgT0ciTqZZvaS6pTkWQZ1LuO/ZC4yfksWBrW7s3uBB2jktezZ7ANcCP/H9/GPc/qQTvh2qDuns/t935GSkkpd1gfNnoutsT/Ul3FPiT5EYc8yiGS3Vp9oajRAf7cDRP1w5utuFrNRLU0ldPEq5alwOV43LIeLHN9i1fhkHf3/YpsHGs3M+5pPXZ1FaWmoKKiq+29nZsXTpUpudW4j2wubDKEuXLuWZZ55pUM+GdDeKlqby1Mfakj6r/5dKz9Oz5dilYRVL1FX8qq6aERdSkjix346TkSEc+8sVlEu5IQ7Oejr2yqLnlSq69CvmP//sDFi21ouzuxcPPDeHovw8HFxcWfWflynIsazOx4+fzmH3hki6DngON++7OXXYibzMS/c5Wp2RXlcWMGBEHh26xFBckFlnLZHKi7JdTq+K1k7Ndd186RfszsGDB82u9hsZGcmgQYMafA5z5H1NtEctKkFUr9ej1+tNP+fm5jZja4Soafny5UyaNMmqu2BfVx0PXhnKH7EZtc5Wqa6h02bnPFK5R8QOGAHcC9xFcYEvJw4EcOJA+aP22kJKDL8DfwIngVNALJBf43kLcjJZ/No/a2yv3tvy759iuJCs5ezfeaTEqUiJ8+LM3x8CnsQeuXSczqmEvsOK6Hd1Pj0HF6B1KH8NZ900st5zVF6UraHBRoiXEzf28se92lo31Vf7FUI0jhYVbMydO5c333yzuZshRK0mTJhAr169zN4F7927t9a7YDuNmpE9/Ojo7czW4ykU6Gsr112TNdNmq846KQW2XfyaDgzDJ+hhvIMmEHfcAX2hjvJZLeOqPUsycPri9/xqX0WAI+BS7SsQB+eB/N8dLrW0rAD4C9gN/IG+cCcTZh+rp/3VXge1mrEPP8Uf68un49ZWJ6QuWjs113bzoV8H9yq9VLWt9uvn51fvcwoh6mfVMMqLL77Iu+++W+c+0dHR9OzZ0/SzNcMo5no2QkJCpLtRtCgVXe7V74It7XIvNJSy9XgqZ9IL6twvOz2FD2bcXWPa7LOf/oiHb0Ctx9U26wQwDUmUlUJ+bjCnj+iIWHMI6AZ0BS7/w9XZvRRH51Qykn8H5RDwB3CY8uCn4eXM61Nf9dQwbydG9/bHrZZy43q9Hq1WawrqDAYDOl3jF/OSYRTRHlnVs/Hcc88xadKkOvfp3Llzgxuj0+ls8p9biMZ0uXfBTlo7bh/YgWNJOew8lY6+xHx3ffWaEcHdevO/RfPJy0qvM9ioS/UhiddW7CBy2/OmgOavTb+QlebI0LH/x/bVPwKulPdcuALOhPW6mvjo3ZT3VORe/J7HzY9OpEe4Pz5BJTi6lF/PuVMuvD/9gxptaGg587oWZaureqrOvjw3o28H9zrPV/m9R6VSyXuREI3IqmDD19cXX19fW7VFiFYhODiYuLg4013wlClTGnQX3LeDO2HeTvx+Iq3WXo7K02QP/LaB00f31ZurUFFIzF7nQFbaeRQzuQcVH9C1FcHKz8nkwG8/VOpVWUx2ejJ3/vNzFr/+UqXtm8lOT2bw6PvJyzrAkjdrJm42Vjnz2hZlqyt46RngynXdfXHWtagRYyHaHZv9D0xISCAzM5OEhATKyso4fPgwAF27dsXFpbZxXSFah8a6C3Z1sOf2gR04mZJHxMk0Cg1V79ytWdOkQuUAIin2eL1JpubqgNRVibO27b//8FWVxM2GVk+t7dzVF2WrK3jxdLLnhp7+hHo71fMbEEI0BZsFG6+99hrffPON6ecrrrgCgO3bt3P99dfb6rRCtEo9AlwJ83Yi4mQ60cmXZmFZs6ZJZdULhzVkbZbaipFV3p6Vdt5sMBTSvR87133LXTNeo/81N5mChrjjh/nuvX/VO23V3LktCV7s1CqGdPJicJgndhq1uacWQjQDKVcuRAuTmFnIjph00vP0FtXb8AvpVGvdiYYmmVpq1k096t2ncjBUufT6ndNesfp8dZV+7+LnwnXdfPBwsrwUfHOQ9zXRHslAphAtTIiXExOGhhKVlIOj9s56622sXfhOrXUnLndhsvrUNVUVQOfkwvE92ynMz8XR1c3ioaDamOvxCHB34NpuPgR7ypCJEC2VBBtCtEAqlYr+wR5093fl2wtnTdsqD4WkJpwGqPcDvCFrs1iqruJjAPrCfL567cka26sPBT376Rqrq4K6OdpzdVdvevi7VqmZIYRoeSTYEKIFc7DXcPOVPfHz98fdJ4CBo+9iz89rSIyJYsW7L1TZ15JcjpamYiioYml3S6qCOmk1DO7oyYBgD8nLEKKVkP+pQrRwwcHBJMTHczLqEJ+/M5vPV2/hwefnotZozO6v1miYMHt+k7WvInEzpHs/Rj0w1ew+tbXnwefn4R/apUrvzLlTf5MYc4zM1KSq59HZMaKHL49d04nwMC8JNIRoRSRBVIhWKCNfz4pNETzzwNgaj81auNaiolmNqSJxs2KqbfUhnwmz57Pi3ResnhXz/q8ncXO0Z0hHT3oHurWJAEPe10R7JMMoQrRCPi46ru1WXmCvIdNaG1tFHkht01N9O3SssT393FkMxUW1VgWd9vp/GNs3gB7+rqjVkpMhRGsmPRtCtFLnzp1jyJAhhISEMHHSo3z+30WcSzzH05+saZRprQ1V2/RUc9tT4k+ZTS797Y+/GHXNVc3QetuT9zXRHknPhhCtVPWy6dP++SQGg4H0wjKOn8/ldHoBxSWWry7bWCwpBlZ9Vsyl9U/KF7XzbOG1MoQQ1mn9A6BCtGM6nc407bOibHqwpxM39QlgynWdufOKDvQJcsPBvmYyaWJMFJ+98AiJMVFN3Ww0ahUdfZwYO6QH/v7+DB48mC+++ILw8HACAgJkaXch2hjp2RCijSr/QHemo48zRqNCYlYh57KKSMoqIiW32Krppo3By1lLkIcjwZ6OdPJxNgVA8fHxl72onRCiZZNgQ4h2QK1WEebtDPkZOKgz6Oim8NbuLQAc3bmZUbfdR35xCQ6uHhZX86yNSgXOWjs8nOzxd3MgyMORIA8HnLTm325kaXch2j5JEBWiHalcabO2WSxn0/PJKjRQZCjDUGbEUGqkpEyhpMxISZkRO40Ke40ae40arUaNnUaFk1aDu6MWd0d7PJzssW8DU1RtRd7XRHskPRtCtCPLly9n0qRJlJaWmgKMiu92dnYsXbq0fOgF5+ZsphCijZFgQ4h2ZMKECfTq1Yvw8PAaj+3du5dBgwY1Q6uEEG2d9HUK0U6p1eoq34UQwlbkXUaIdsbPz4+AgADCw8NluqkQoklIgqgQ7ZBerzdNN1UURaabNiF5XxPtkeRsCNEOyXRTIURTkmEUIYQQQtiUBBtCCCGEsCkJNoQQQghhUxJsCCGEEMKmJNgQQgghhE1JsCGEEEIIm5JgQwghhBA2JcGGEEIIIWxKgg0hhBBC2JQEG0IIIYSwKZsFG3FxcUyePJlOnTrh6OhIly5deP311zEYDLY6pRBCCCFaIJutjXLixAmMRiNffvklXbt25dixYzzxxBMUFBSwYMECW51WCCGEEC1Mk676On/+fD7//HPOnDlj0f6yOqIQoq2R9zXRHjXpqq85OTl4eXnV+rher0ev15t+zs3NbYpmCSGEEMKGmixBNDY2lk8++YSpU6fWus/cuXNxd3c3fYWEhDRV84QQQghhI1YHGy+++CIqlarOrxMnTlQ5JikpibFjx3LvvffyxBNP1PrcL730Ejk5OaavxMRE669ICCGEEC2K1Tkb6enpXLhwoc59OnfujFarBeD8+fNcf/31XHXVVSxduhS12vL4RsY2hRBtjbyvifbI6pwNX19ffH19Ldo3KSmJkSNHEh4ezpIlS6wKNIQQQgjRNtgsQTQpKYnrr7+esLAwFixYQHp6uumxgIAAW51WCCGEEC2MzYKNrVu3EhsbS2xsLMHBwVUea8LZtkIIIYRoZjYb15g0aRKKopj9EkIIIUT7IUkUQgghhLApCTaEEEIIYVMSbAghhBDCpiTYEEIIIYRNSbAhhBBCCJuSYEMIIYQQNiXBhhBCCCFsSoINIYQQQtiUBBtCCCGEsCkJNoQQQghhUxJsCCGEEMKmJNgQQgghhE1JsCGEEEIIm5JgQwghhBA2JcGGEEIIIWxKgg0hhBBC2JQEG0IIIYSwKQk2hBBCCGFTEmwIIYQQwqYk2BBCCCGETUmwIYQQQgibkmBDCCGEEDYlwYYQQgghbEqCDSGEEELYlAQbQgghhLApCTaEEEIIYVMSbAghhBDCpmwabNx2222Ehobi4OBAYGAgDz/8MOfPn7flKYUQQgjRwtg02Bg5ciQ//PADJ0+e5Mcff+T06dPcc889tjylEEIIIVoYlaIoSlOdbMOGDdxxxx3o9Xrs7e3r3T83Nxd3d3dycnJwc3NrghYKIYRtyfuaaI/smupEmZmZrFixguHDh9caaOj1evR6vennnJwcoPw/pxBCtAUV72dNeJ8nRLOzebAxe/ZsPv30UwoLC7nqqqvYuHFjrfvOnTuXN998s8b2kJAQWzZRCCGaXF5eHu7u7s3dDCGahNXDKC+++CLvvvtunftER0fTs2dPADIyMsjMzCQ+Pp4333wTd3d3Nm7ciEqlqnFc9Z4No9FIZmYm3t7eZvevTW5uLiEhISQmJrbZbsq2fo1yfa1fW7/Ghl6foijk5eURFBSEWi0TAkX7YHWwkZ6ezoULF+rcp3Pnzmi12hrbz507R0hICH/++SfDhg2zrqVWaA9jom39GuX6Wr+2fo1t/fqEaExWD6P4+vri6+vboJMZjUaAKr0XQgghhGjbbJazsXfvXvbv388111yDp6cnp0+f5tVXX6VLly427dUQQgghRMtiswFDJycn1q5dy6hRo+jRoweTJ0+mf//+7NixA51OZ6vTAqDT6Xj99ddtfp7m1NavUa6v9Wvr19jWr0+IxtSkdTaEEEII0f5IKrQQQgghbEqCDSGEEELYlAQbQgghhLApCTaEEEIIYVOtNthYuHAhHTt2xMHBgaFDh7Jv374691+9ejU9e/bEwcGBfv36sXnz5iZqacNZc42LFi3i2muvxdPTE09PT0aPHl3va9LcrP0dVli5ciUqlYo77rjDtg28TNZeX3Z2NtOnTycwMBCdTkf37t1b/N+ptdf44Ycf0qNHDxwdHQkJCeHZZ5+luLi4iVprnZ07dzJ+/HiCgoJQqVT89NNP9R4TERHBoEGD0Ol0dO3alaVLl9q8nUK0CkortHLlSkWr1Spff/218vfffytPPPGE4uHhoaSmpprdf/fu3YpGo1Hee+895fjx48orr7yi2NvbK1FRUU3ccstZe40PPfSQsnDhQuXQoUNKdHS0MmnSJMXd3V05d+5cE7fcMtZeX4WzZ88qHTp0UK699lrl9ttvb5rGNoC116fX65XBgwcrN998s7Jr1y7l7NmzSkREhHL48OEmbrnlrL3GFStWKDqdTlmxYoVy9uxZ5ZdfflECAwOVZ599tolbbpnNmzcrL7/8srJ27VoFUNatW1fn/mfOnFGcnJyUWbNmKcePH1c++eQTRaPRKFu2bGmaBgvRgrXKYOPKK69Upk+fbvq5rKxMCQoKUubOnWt2//vuu0+55ZZbqmwbOnSoMnXqVJu283JYe43VlZaWKq6urso333xjqyZeloZcX2lpqTJ8+HDlq6++UiZOnNiigw1rr+/zzz9XOnfurBgMhqZq4mWz9hqnT5+u3HDDDVW2zZo1S7n66qtt2s7GYEmw8a9//Uvp06dPlW3333+/MmbMGBu2TIjWodUNoxgMBiIjIxk9erRpm1qtZvTo0fz1119mj/nrr7+q7A8wZsyYWvdvbg25xuoKCwspKSnBy8vLVs1ssIZe31tvvYWfnx+TJ09uimY2WEOub8OGDQwbNozp06fj7+9P3759+fe//01ZWVlTNdsqDbnG4cOHExkZaRpqOXPmDJs3b+bmm29ukjbbWmt7nxGiKdl8ifnGlpGRQVlZGf7+/lW2+/v7c+LECbPHpKSkmN0/JSXFZu28HA25xupmz55NUFBQjTe/lqAh17dr1y4WL17M4cOHm6CFl6ch13fmzBl+//13JkyYwObNm4mNjWXatGmUlJTw+uuvN0WzrdKQa3zooYfIyMjgmmuuQVEUSktLefLJJ/m///u/pmiyzdX2PpObm0tRURGOjo7N1DIhml+r69kQ9Zs3bx4rV65k3bp1ODg4NHdzLlteXh4PP/wwixYtwsfHp7mbYxNGoxE/Pz/++9//Eh4ezv3338/LL7/MF1980dxNazQRERH8+9//5rPPPuPgwYOsXbuWTZs28fbbbzd304QQNtbqejZ8fHzQaDSkpqZW2Z6amkpAQIDZYwICAqzav7k15BorLFiwgHnz5vHbb7/Rv39/Wzazway9vtOnTxMXF8f48eNN2ypWELazs+PkyZN06dLFto22QkN+f4GBgdjb26PRaEzbevXqRUpKCgaDAa1Wa9M2W6sh1/jqq6/y8MMP8/jjjwPQr18/CgoKmDJlCi+//DJqdeu+96ntfcbNzU16NUS71+r+d2u1WsLDw9m2bZtpm9FoZNu2bbWuJjts2LAq+wNs3bq1xa4+25BrBHjvvfd4++232bJlC4MHD26KpjaItdfXs2dPoqKiOHz4sOnrtttuY+TIkRw+fJiQkJCmbH69GvL7u/rqq4mNjTUFUQAxMTEEBga2uEADGnaNhYWFNQKKiuBKaQNLNLW29xkhmlRzZ6g2xMqVKxWdTqcsXbpUOX78uDJlyhTFw8NDSUlJURRFUR5++GHlxRdfNO2/e/duxc7OTlmwYIESHR2tvP76661i6qs11zhv3jxFq9Uqa9asUZKTk01feXl5zXUJdbL2+qpr6bNRrL2+hIQExdXVVZkxY4Zy8uRJZePGjYqfn5/yzjvvNNcl1Mvaa3z99dcVV1dX5fvvv1fOnDmj/Prrr0qXLl2U++67r7kuoU55eXnKoUOHlEOHDimA8v777yuHDh1S4uPjFUVRlBdffFF5+OGHTftXTH194YUXlOjoaGXhwoUy9VWIi1plsKEoivLJJ58ooaGhilarVa688kplz549psdGjBihTJw4scr+P/zwg9K9e3dFq9Uqffr0UTZt2tTELbaeNdcYFhamADW+Xn/99aZvuIWs/R1W1tKDDUWx/vr+/PNPZejQoYpOp1M6d+6szJkzRyktLW3iVlvHmmssKSlR3njjDaVLly6Kg4ODEhISokybNk3Jyspq+oZbYPv27Wb/T1Vc08SJE5URI0bUOGbgwIGKVqtVOnfurCxZsqTJ2y1ESyRLzAshhBDCplpdzoYQQgghWhcJNoQQQghhUxJsCCGEEMKmJNgQQgghhE1JsCGEEEIIm5JgQwghhBA2JcGGEEIIIWxKgg0hhBBC2JQEG0IIIYSwKQk2hBBCCGFTEmwIIYQQwqYk2BBCCCGETf0/xpS2GApKSVkAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "model.eval()\n", "with torch.no_grad():\n", " trained_pred_dist = likelihood(model(test_x))\n", " predictive_mean = trained_pred_dist.mean\n", " lower, upper = trained_pred_dist.confidence_region(rescale=True)\n", " \n", "f, ax = plt.subplots(1, 1, figsize=(4, 3))\n", "# Plot training data as black stars\n", "ax.plot(train_x, train_y, 'k*')\n", "# Plot predictive means as blue line\n", "ax.plot(test_x, predictive_mean, 'b')\n", "# Shade between the lower and upper confidence bounds\n", "ax.fill_between(test_x, lower, upper, alpha=0.5)\n", "ax.set_ylim([-3, 3])\n", "ax.legend(['Observed Data', 'Mean', 'Confidence'], bbox_to_anchor=(1.6,1));" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now our model seems to fit well on the data. It is not always possible to visually evaluate the model in high dimensional cases. Thus, now we evaluate the models with help of probabilistic metrics. We have saved predictive distributions from untrained and trained models as `untrained_pred_dist` and `trained_pred_dist` respectively." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Negative Log Predictive Density (NLPD)\n", "\n", "Negative Log Predictive Density (NLPD) is the most standard probabilistic metric for evaluating GP models. In simple terms, it is negative log likelihood of the test data given the predictive distribution. It can be computed as follows:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Untrained model NLPD: 1.45, \n", "Trained model NLPD: -0.29\n" ] } ], "source": [ "init_nlpd = qpytorch.metrics.negative_log_predictive_density(untrained_pred_dist, test_y)\n", "final_nlpd = qpytorch.metrics.negative_log_predictive_density(trained_pred_dist, test_y)\n", "\n", "print(f'Untrained model NLPD: {init_nlpd:.2f}, \\nTrained model NLPD: {final_nlpd:.2f}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Mean Standardized Log Loss (MSLL)\n", "\n", "This metric computes average negative log likelihood of all test points w.r.t their univariate predicitve densities. For more details, see \"page No. 23, Gaussian Processes for Machine Learning, Carl Edward Rasmussen and Christopher K. I. Williams, The MIT Press, 2006. ISBN 0-262-18253-X\"" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Untrained model MSLL: 0.88, \n", "Trained model MSLL: 13.94\n" ] } ], "source": [ "init_msll = qpytorch.metrics.mean_standardized_log_loss(untrained_pred_dist, test_y)\n", "final_msll = qpytorch.metrics.mean_standardized_log_loss(trained_pred_dist, test_y)\n", "\n", "print(f'Untrained model MSLL: {init_msll:.2f}, \\nTrained model MSLL: {final_msll:.2f}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is also possible to calculate the quantile coverage error with `qpytorch.metrics.quantile_coverage_error` function." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Quantile 95% Coverage Error: 0.56\n" ] } ], "source": [ "quantile = 95\n", "qce = qpytorch.metrics.quantile_coverage_error(trained_pred_dist, test_y, quantile=quantile)\n", "print(f'Quantile {quantile}% Coverage Error: {qce:.2f}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Mean Squared Error (MSE)\n", "\n", "Mean Squared Error (MSE) is the mean of the squared difference between the test observations and the predictive mean. It is a well-known conventional metric for evaluating regression models. However, it can not take uncertainty into account unlike NLPD, MLSS and ACE." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Untrained model MSE: 0.19, \n", "Trained model MSE: 0.03\n" ] } ], "source": [ "init_mse = qpytorch.metrics.mean_squared_error(untrained_pred_dist, test_y, squared=True)\n", "final_mse = qpytorch.metrics.mean_squared_error(trained_pred_dist, test_y, squared=True)\n", "\n", "print(f'Untrained model MSE: {init_mse:.2f}, \\nTrained model MSE: {final_mse:.2f}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Mean Absolute Error (MAE)\n", "\n", "Mean Absolute Error (MAE) is the mean of the absolute difference between the test observations and the predictive mean. It is less sensitive to the outliers than MSE." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Untrained model MAE: 0.37, \n", "Trained model MAE: 0.13\n" ] } ], "source": [ "init_mae = qpytorch.metrics.mean_absolute_error(untrained_pred_dist, test_y)\n", "final_mae = qpytorch.metrics.mean_absolute_error(trained_pred_dist, test_y)\n", "\n", "print(f'Untrained model MAE: {init_mae:.2f}, \\nTrained model MAE: {final_mae:.2f}')" ] } ], "metadata": { "interpreter": { "hash": "d4d1e4263499bec80672ea0156c357c1ee493ec2b1c70f0acce89fc37c4a6abe" }, "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.12.12" } }, "nbformat": 4, "nbformat_minor": 4 }