BayesianNetwork.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "BayesianNetwork.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "authorship_tag": "ABX9TyPtADtvrylzNZ+iQCs+RddZ",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/EastbayML/VariationalInference/blob/main/BayesianNetwork.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "CrgFhOUVd7eA",
        "outputId": "f5487d59-2772-4549-f10c-a1def42bf3c8"
      },
      "source": [
        "import sys\n",
        "print(\"python:\",sys.version)\n",
        "import torch\n",
        "print(\"torch:\",torch.__version__)"
      ],
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "python: 3.7.12 (default, Sep 10 2021, 00:21:48) \n",
            "[GCC 7.5.0]\n",
            "torch: 1.10.0+cu111\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "CMiOinONPxSv"
      },
      "source": [
        "#%matplotlib inline\n",
        "import numpy as np\n",
        "import torch\n",
        "from torch import nn\n",
        "from sklearn import datasets\n",
        "from torch import dist\n",
        "import matplotlib.pyplot as plt\n",
        "#from IPython.display import clear_output\n",
        "import copy\n",
        "from dataclasses import dataclass\n"
      ],
      "execution_count": 1,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "s7LwK2qXP3mu"
      },
      "source": [
        "def plotdist(m,X,Y):# draw samples from Q(theta)\n",
        "    plt.figure(figsize=(16, 6))\n",
        "    plt.scatter(X, Y)\n",
        "    if m:\n",
        "        with torch.no_grad():\n",
        "          y_pred = torch.cat([m(X)[0] for _ in range(1000)], dim=1)\n",
        "        \n",
        "        # Get some quantiles\n",
        "        q1, mu, q2 = np.quantile(y_pred, [0.05, 0.5, 0.95], axis=1)\n",
        "\n",
        "\n",
        "        plt.plot(X, mu)\n",
        "        plt.fill_between(X.flatten(), q1, q2, alpha=0.2)"
      ],
      "execution_count": 3,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 374
        },
        "id": "rkDAdFPkPyBn",
        "outputId": "807d23be-3c8f-4dda-e0a2-7d17a7915295"
      },
      "source": [
        "import numpy as np\n",
        "import torch\n",
        "from torch import nn\n",
        "from sklearn import datasets\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "w0 = 0.125\n",
        "b0 = 5.\n",
        "x_range = [-20, 60]\n",
        "\n",
        "np.random.seed(43)\n",
        "\n",
        "def load_dataset(n=150, n_tst=150):\n",
        "    def s(x):\n",
        "        g = (x - x_range[0]) / (x_range[1] - x_range[0])\n",
        "        return 3 * (0.25 + g**2.)\n",
        "\n",
        "    x = (x_range[1] - x_range[0]) * np.random.rand(n) + x_range[0]\n",
        "    eps = np.random.randn(n) * s(x)\n",
        "    y = (w0 * x * (1. + np.sin(x)) + b0) + eps\n",
        "    y = (y - y.mean()) / y.std()\n",
        "    idx = np.argsort(x)\n",
        "    x = x[idx]\n",
        "    y = y[idx]\n",
        "    return y[:, None], x[:, None]\n",
        "\n",
        "y, x = load_dataset()\n",
        "plotdist(None,x,y)"
      ],
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6AAAAFlCAYAAAD1f4WfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3df4wk2V0Y8O+7ubHd/IgH8AW8Yzu3SZx1bBbfhpEDGhThA7IGG1iWgEEkIgmSFSmJAJGNduNIQCS0E62CQVEi5YQJSHEAB68Hiw1ZDHMS4RQbZhmbtX3exMEB3GfwEhh+eWLvrV/+mJnb2dnume6Z6qpXVZ+PZHmnuq/7db/qqvrW+77vSznnAAAAgFl7qOkGAAAA0A8CUAAAAGohAAUAAKAWAlAAAABqIQAFAACgFgJQAAAAavFwE2/6ohe9KD/66KNNvDUAAAAzduPGjT/IOT+yf3sjAeijjz4a6+vrTbw1AAAAM5ZS+u1R26XgAgAAUAsBKAAAALUQgAIAAFALASgAAAC1EIACAABQCwEoAAAAtRCAAgAAUAsBKAAAALUQgAIAAFCLh5tuAAAAQBesbgzjyvVb8czmVpxYGMSFs6fi3JnFpptVFAEoAADAMa1uDOPS1ZuxdeduREQMN7fi0tWbERGC0D2k4AIAABzTleu3ngs+d23duRtXrt9qqEVlEoACAAAc0zObW1Nt7ysBKAAAwDGdWBhMtb2vBKAAAADHdOHsqRjMz923bTA/FxfOnmqoRWVShAgAAOCYdgsNqYJ7MAEoAABABc6dWRRwHkIKLgAAALUQgAIAAFALASgAAAC1EIACAABQCwEoAAAAtRCAAgAAUAsBKAAAALUQgAIAAFALASgAAAC1EIACAABQCwEoAAAAtRCAAgAAUAsBKAAAALUQgAIAAFALASgAAAC1OHYAmlJ6QUrp11JK708pfTCl9INVNAwAAIBuebiC1/hURDyec/6zlNJ8RPxqSukXcs7vqeC1AQAA6IhjB6A55xwRf7bz5/zO//JxXxcAAIBuqWQOaEppLqX0voj4RES8O+f83ipeFwAAgO6oJADNOd/NOT8WES+JiNeklL54/3NSSm9KKa2nlNZv375dxdsCAADQIpVWwc05b0bEkxHxuhGPPZFzXso5Lz3yyCNVvi0AAAAtUEUV3EdSSgs7/x5ExNdExIeP+7oAAAB0SxVVcF8cET+ZUpqL7YD27Tnnn6/gdQEAAOiQKqrg/mZEnKmgLQAAAHRYpXNAAQAAYBwBKAAAALUQgAIAAFALASgAAAC1EIACAABQCwEoAAAAtRCAAgAAUAsBKAAAALUQgAIAAFALASgAAAC1EIACAABQCwEoAAAAtRCAAgAAUAsBKAAAALUQgAIAAFCLh5tuAAAA0F2rG8O4cv1WPLO5FScWBnHh7Kk4d2ax6WbREAEoAAAwE6sbw7h09WZs3bkbERHDza24dPVmRIQgtKek4AIAADNx5fqt54LPXVt37saV67caahFNMwIKAADMxDObW1Nt7zKpyNuMgAIAADNxYmEw1fau2k1FHm5uRY57qcirG8Omm1Y7ASgAADATF86eisH83H3bBvNzceHsqYZa1AypyPdIwQUAAGZiN8W076mnUpHvEYACAAAzc+7MYu8Czv1OLAxiOCLY7FsqcoQUXAAAgJmSinyPEVAAAIAZkop8jwAUAABgxqQib5OCCwAAQC0EoAAAANRCAAoAAEAtBKAAAADUQhEiAAAqtboxVO0TGEkACgBAZVY3hnHp6s3YunM3IiKGm1tx6erNiAhBKCAFFwCA6ly5fuu54HPX1p27ceX6rYZaBJREAAoAQGWe2dyaajvQLwJQAAAqc2JhMNV2oF8EoAAADVrdGMbyylqcvHgtllfWYnVj2HSTjuXC2VMxmJ+7b9tgfi4unD3VUIuAkihCBADQkC4W7Nlttyq4wCgCUACAhhxUsKfNAdu5M4utbj8wO1JwAQAaomAP0DcCUACAhijYA/SNABQAoCEK9gB9Yw4oAEBDFOwB+kYACgDQIAV7gD4RgAIAAI1Y3RjKAOgZASgAAFC7Lq6Dy+EEoAAAwMztH+385Kef7eQ6uBxMAAoAAMzUqNHOcayD222WYQEAAGbqyvVbD4x2jmMd3G4TgAIAADOzujE8cMRzL+vgdp8UXAAAqNi01V27Wg12N/V2nIXBfHz28x/u3OdmPAEoAABUaNrqrl2uBntQ6u1gfi5+4Bte1frPyHSOnYKbUnppSunJlNKHUkofTCl9dxUNAwCANhoVdO1Wd63i+W1yUEGhy+dPCz57qIoR0Gcj4vtyzr+RUvrciLiRUnp3zvlDFbw2AAAd1dW003FBV1Xb2+TEwmDk/M/FhUEn+prpHXsENOf88Zzzb+z8+08j4umIsDcBADDWbtrpcHMrctxLO13dGDbdtGMbV8W1qu1tcuHsqRjMz923TaGhfqu0Cm5K6dGIOBMR7x3x2JtSSusppfXbt29X+bYAALRMl9NOpw26uhyknTuzGJfPn47FhUGk2B75lHrbb5UVIUopfU5EvCMivifn/Cf7H885PxERT0RELC0t5areFwCA9uly2ulucDVpevG0z2+bc2cWO/NZOL5KAtCU0nxsB59vyzlfreI1AQDornFzA7uQdhoxfdAlSKMvqqiCmyLirRHxdM75h4/fJAAAuq7LaafAeFXMAV2OiL8XEY+nlN6387+vq+B1AQDoKHMDoZ+OnYKbc/7ViEgVtAUAgB6Rdgr9U2kVXAAAABhHAAoAAEAtKluGBQAA4KhWN4adXYqGewSgAABAo1Y3hnHp6s3YunM3IiKGm1tx6erNiAhBaMcIQAEAgEZduX7rueBz19adu3Hl+q1WBqBGc8cTgAIAAI16ZnNrqu0lM5p7MEWIAACAI1ndGMbyylqcvHgtllfWYnVjeKTXObEwmGp7yQ4azUUACgAAHMHuSN9wcyty3BvpO0oQeuHsqRjMz923bTA/FxfOnqqotfXp0mjuLAhAAQCoTFUjYpSvypG+c2cW4/L507G4MIgUEYsLg7h8/nQrU1a7NJo7C+aAAgBQCXPf+qXqkb5zZxY7sZ9cOHvqvt9BRHtHc2fBCCgAAJUw961fjPSN1qXR3FkwAgoAQCXMfesXI33jdWU0dxaMgAIAUAkjYv1ipI+jMAIKAEAljIj1T9Ujfasbw7hy/VY8s7kVJxYGceHsKQFtxwhAAQCoxG6gIIDgKBSx6gcBKAAAlTH3jaM6qIiVfao7zAEFAAAap4hVPwhAAQCAxili1Q8CUAAAoHEXzp6KwfzcfdsUseoec0ABAIDGKWLVDwJQAACgCIpYdZ8UXAAAAGohAAUAAKAWAlAAAABqYQ4oAFCE1Y2h4iMAHScABQAat7oxjEtXb8bWnbsRETHc3IpLV29GRAhCocPceOofASgA0Lgr1289F3zu2rpzN65cv+VitIMEHUS48dRX5oACAI17ZnNrqu20127QMdzcihz3go7VjWHTTaNmB914orsEoABA404sDKbaTnsJOtjlxlM/CUABgMZdOHsqBvNz920bzM/FhbOnGmoRsyLoYJcbT/0kAAUAGnfuzGJcPn86FhcGkSJicWEQl8+fNg+sgwQd7HLjqZ8UIQIAinDuzGInAk4Fdg524eyp+wrPRAg6+mr3d+H30i8CUACAiqjqeThBB3t15cYTkxOAAgBUxHIykxF0QH+ZAwoAUBEFdgAOJgAFAKiIAjsABxOAAgBURFVPgIOZAwoAUBEFdgAOJgAFAKiQAjsA40nBBQAAoBYCUAAAAGohAAUAAKAW5oACANRodWPYuyJFffzM9If9ezoCUACAmqxuDOPS1ZuxdeduREQMN7fi0tWbERGdvWDt42emP+zf05OCCwBQkyvXbz13obpr687duHL9VkMtmr0+fmb6w/49PQEoAEBNntncmmp7F/TxM9Mf9u/pCUABAGpyYmEw1fYu6ONnpj/s39MTgAIAtVjdGMbyylqcvHgtllfWYnVj2HSTanfh7KkYzM/dt20wPxcXzp5qqEWz18fPTH/Yv6enCBEAMHMKdWzb/ax9qpjZx89Mf9i/p5dyzrW/6dLSUl5fX6/9fQGAZiyvrMVwxJyoxYVBPHXx8QZaBMAspZRu5JyX9m+vJAU3pfTjKaVPpJQ+UMXrAQDdolAHABHVzQH9iYh4XUWvBQB0jEId3WVuLzCNSgLQnPOvRMQfVvFaAED3KNTRTbtze4ebW5Hj3txeQSgwTm1VcFNKb0opraeU1m/fvl3X2wIABTh3ZjEunz8diwuDSLE99/Py+dMKdbTY6sYwvu/t73+usNSurTt348r1Ww21CihdbVVwc85PRMQTEdtFiOp6XwCgDOfOLAo4O2J35PPumGKW5vYC41gHFACAqVy5fuuBkc+9zO0FxhGAAgAwlYNGOM3tBQ5S1TIsPxUR/yMiTqWUPpZS+q4qXhcAgPKMG+GcS8ncXuBAVVXB/fac84tzzvM555fknN9axesCAFCecVWN/823vlrwCRyotiJEAAC7VjeGceX6rXhmcytOLAziwtlTApcW2e0rfQhMSwAKANRqt4LqbhGb3bUjI0IAU6hxNwz0FzAtASgAUKtRFVR3144U0JTHDYNuk41A3VTBBQBqNa6CqrUjy3TQDQPabffmwnBzK3Lcu7mwujFsuml0mAAUAKjVuAqq1o4skxsG3eXmAk0QgAIAtRpXQdXakWVyw6C73FygCQJQAKBW584sxuXzp2NxYRApIhYXBtaOLJgbBt3l5gJNUIQIAKidCqrtYcmV7rpw9tR9BaYi3Fxg9gSgAAAcyA2De7pUNdbNBZogAAUAgAl0cUkaNxeomwAUAGBCXRr9YnrWsIXjE4ACAEygi6NfTEfVWDg+VXABACZgzURUjYXjE4ACAEzA6BeWpKFpqxvDWF5Zi5MXr8XyylqsbgybbtLUpOACAEzgxMIghiOCTaNf/aFqLE3qyjQAASgAUJwSi/1YM5EIVWNpTleKYAlAAYCilHqX3+gX0KSuTAMQgAIARSn5Lr/RL6ApXZkGoAgRAFCUrtzlB9qvpKI/XSmCZQQUAChKV+7yA+1W2nSArkwDEIACAEVR7AcoQYnTAbowDUAACgAUpSt3+YF2Mx1gNgSgAEBxunCXH2g30wFmQwAKADNW4pqWcFyz2K/9ViiJ6QCzIQAFgBkqrYgFVGEW+7XfCqUxHWA2Us659jddWlrK6+vrtb8vANRteWVtZArX4sIgnrr4eAMtguObxX7ttwIPanNWQErpRs55af92I6AAMEOKWNBGh130zmK/9luB+3U1K0AACgAz1OUiFm2+M894k1z0zmK/Lu23Yv+maSUuA1OFh5puAAB02YWzp2IwP3fftqMUsVjdGMbyylqcvHgtllfWYnVjWGUzp7YbpAw3tyLHvSCl6XZxfAdd9O6qar/eaxaveVT2b0rQ1awAASgAzNC5M4tx+fzpWFwYRIrt+WyXz5+e6u51iRfDkwQppSstqC/FJBe9VezX+83iNY+qC/s37Tdu9L/tGTRScAFgxo67pmWJaVhtvzPf1blVVZg0FXYWa7WWsv5r2/dvuqGry8AIQAGgZtPOLSvxYri0+XrTKi2oL2m+YVcveqfR9v27DiXts13V1WVgBKAAFK1rFzlHGXkr8WK47UFKSUF9aaOxXb3onUbp+3fTx8XS9tkuKyUroEoCUACK1cWLnKOMvJV4Mdz2IKWkoL600diIbl70TqPk/buE42KJ+yztIQAFoFhdvMg5yshbqRfDbQ5SSgrqSxqN5Z5S9+8Sjov2WY5DAApAsbp4kTPtyNv+VLu3vPGxIi+K26akoL6k0VjKV8Jx0T7LcQhAATiSOuYgdfEiZ5qRtxJS7Up23H2wlBGukkZjKV8Jx0X7LMdhHVAAplbXupQlLUxflWnWOrQW4Xglro16VCWtf0n5Sjgu2mc5jpRzrv1Nl5aW8vr6eu3vC0A1llfWRt6BX1wYxFMXH5/qtQ4bxWq62mOTTl68FqPO0ikiPrry+rqbU5Qq90Fomz4fF2mPlNKNnPPS/u1ScAGYWlVzkCZJMS0lTbIJJaTa1W3SC+tR38tB26FL+nxcpP2k4AIwtXEB0LSBkRTTg5WQalenadJq51Ia+RrjtkMXrW4MY3llLU5evBbLK2utTEEvne+4egJQAKZWVWBUQjXHkvVtntU0NyTujplCNG4797ig7oYuzYMule94NqTgAjC1qpawaHOKaV1zsPqUajfNDYnFMfvOYgv2nSaprNwdJawH2nW+49kQgAJwJFUERq99xSPxtvf8zn2FdtqQYuoifjamuSFhGYijcUHdHTJIZs93PBtScAFoxOrGMN5xY3hf8Jki4pu/tPwRP3NXZ2Oa1O6+pSdXxQV1d1Q1F5/xfMezYQQUgEaMCuJyRDz54dvNNGgKLuJnY9rU7j6lJ1elzWnv3E8WwOz5jmdDAArAgWY113GaIK60Ne9cxM+OoHK2XFB3R1Vz8RnPdzwbAlAAxprlXMdJg7gS51u6iKcqdd9ccUHdLW7YzJ7vuHoCUADGmmXBkkmDuBKLpriIpwpN3VxxQQ00SQAKwFiznOs4aRBX6nxLF/EcV4k3VwBmrZIANKX0uoj40YiYi4gfyzmvVPG6ADRr1nMdJwnizLekq0q9uQIwS8dehiWlNBcR/y4ivjYiXhkR355SeuVxXxeA5k2zLEaX2wCzYIkHoI+qWAf0NRHxkZzzb+WcPx0RPx0R31jB6wLQsBLWWiyhDTALbq4AfVRFCu5iRPzunr8/FhF/c/+TUkpviog3RUS87GUvq+BtAZpX2vIgs1DCXMcS2gBVU8wKmtGHc3fJaitClHN+IiKeiIhYWlrKdb0vwKyUuDxI3fp2Eu/b52X23FyBejl3N6+KFNxhRLx0z98v2dkG0GkHVbDsg92T+HBzK3LcO4mvbrTzFLC6MYzllbU4efFaLK+sPfA5uvZ5Afqo7+fuElQRgP56RLw8pXQypfS8iPi2iHhXBa8LULS+V7Ds0kl8kuCyS58XoK/6fu4uwbED0JzzsxHxTyLiekQ8HRFvzzl/8LivC1C6vlew7NJJfJLgskufl4MdNhoOtFffz90lqGIENHLO/zXn/Ndyzn8l5/xDVbwmQOn6XsGySyfxSYLLLn3erqkyYJRqTVu5cTKZvp+7S1BJAArQR31fHqRLJ/FJgssufd5JteGCtuqAUao1beTGyeT6fu4uQW1VcAG6qM8VLOtcQmLW1WcvnD11X1XEiAeDy74tmdGWSpEHBYxHaadUa9qo6t9B1/X53F0CASgAR1bHSbyOQGjS4LJPFy1tuaCtOmA8sTCI4Yj/Vqo1JXPjhDYRgAJwoKbXvqwrEOpTcDmJtlzQVh0wTjIaDqVx44Q2EYAC8Jz9weZrX/FIvOPGsNE0zCYDoaaD7ya15YK26oCxb6nWdIMbJ7SJABSAiBid6vqf3vM7Dzyv7jTMpgKhtsyBnJW2XNDOImA0Gk7buHFCmwhAAYiI0amu49SZhtlUINSWOZCz0qYLWgEj+B3QHgJQACJiuqCyzjTMpgKhtsyBnMRRU4ld0AJQNQEo0Cp9npM3a+NSXUepOw2ziUCoLXMgD9P3VGIAyvJQ0w0AmJSFtmfrwtlTMZifO/R5C4P5iIhYXlmLkxevxfLKWif7YNT3UeIcyMMclEoMAHUzAgq0Rt/n5E3iOCPE+1NdXziYjz//9LNx525+7jmD+bl4w6tf3IsRtTbNgTxIl1KJOZwsEaB0AlCgNSzHcbAqUi33p7qO+tx9uhHQhTmQXUkl5nDSrWE6h53b23DubyMBKNAas7qQnuQE1IaLulkEhqMCsO/9mfeNfG6XR9TafBHSluVUOL4+3RyC4zrs3N6Wc38bmQMKtMYs5uRNMq+0LXPo6hohHhfwd3VEre1zj8+dWYzL50/H4sIgUkQsLgzi8vnTnb6AWt0Ydn6O8ijSrWFyh53b23LubyMjoEBrzGJO3iQjBm25qKsr1bJvI2pdGFXqQirxpPo8aiHdGiZ32Lm9Lef+NjICCrTKuTOL8dTFx+OjK6+Ppy4+fuwLyklOMG0Z8auramvfRtRchLRLn0ctulK5Gepw2Lm9Lef+NhKAAr02yQmmLRd1dQaGVd8IKJmLkHbp8w2Dvt0cguM47NzelnN/G0nBBR7Q5oIr05oknbRNy3H0KdWyLn1LOW67vqehOgbAZA47t7fp3N82Ked8+LMqtrS0lNfX12t/X+Bw++dPRWxfbHf5LnrXA+6uf746+A7bo4/HMIASpZRu5JyXHtguAL2fiwz6bnllbeToweLCIJ66+HgDLSpPm44TLsbpozb9RgG6alwAKgV3jz5XzoNdfZ4/NYm2HSe6UMEVpiUNFaBcihDt0efKebCrioIrXV6Dr23HCTcUAICSGAHdo28XalKUytdEHx234ErbRgin1bbjRN8LsgAAZTECukefSu3vBgnDza3IcS9I6NJIVds11UfHLePfthHCabXtOKGMPABQEiOge/Sp1L55YeVrso+OM3+qbSOE02rbcUIZeQAmJTuOOghA9+jqhdqog0nXg4QuaEMfjdq3pkn5bOOJro3HCQVZ+muS31gbf4dA9bo+hYZyWIal48YtwfCC+Yfijz5554HnW2qjHKUvhzJu3/rmL12Md9wYHrrsR1PLg7jYpi8m+Y1ZpgfYVfp1B+0zbhkWc0A7blwaZ85hXljhSp+7N27fevLDtyeaQ9rEXFFzn+mTSX5jXZ+zDUyuDZlXdIMU3I4bd9D446078ZY3PmYkqGClp3oedKKaJOWziROduc/0ySS/MRecwC5V06mLALTjxh1MXjiYLzaw4Z6S5+4d90TVxInOxTZ9MslvzAUnXWfaxeTaVmSP9pKC2yGrG8NYXlmLkxevxfLKWqxuDEemcc4/lOLPP/2sNMTCjerPkt7zuCnCTaQYt20JFTiOSX5jpaf6U50mzilNM+1iOsddhg0mpQhRRxxUSCLi/jTOT3762dYUIOrLncv9n/O1r3hkokI+Vbdh2mIkx+2fuvtXwRX6RhVcIvp77FNUB5o1rgiRALQjpjnInrx4LUb1eoqIj668fjYNPIK+nDBHfc4UMbKPZnnS7MuJ2sU20Dd9Ob7v15brHeiqcQGoOaAdMc3ctrbM+elLwZhRn3PcbaFZzlUc99rDza04efFaUcHacYLIkufVAsxCX+e/t+V6B/pGADpCG0dIpjnIznKSeZXfXV9OmNN8nr39Oeq7jjh61dxx+1BE3Dd3JqLZBaktlA0wnb4GYorqQJkEoPu09eJ2moPsrJb3qPq768sJc9zn3J+Gu7c/R33XF/7L+yNSxJ27+blt03z/o/ah/faOQDd1o+awdQvbdvMIYNb6GoiVvpxZ3do4wEI3mQO6T5vnSTR9YKn6u+vzHNDB/Fx885cuxpMfvj2yP8d916NM8/3v3YfGHRlSRLzljY811jfj5vTstqHr+wvAUTR9jUCz+nJNRVnMAZ1Qm9M+m57bVvV315c7l0f5nNN8p9M8d+8+NC7IPbEwaHR+7rgR47mUejFnGOAomr5GoFl9qatBOwhA9+lL2ucsvHAwH5tbDy7v8sLB/JFfsy8nzGk/50HzNUc99ygOStn63p9538j/po4bNePaNS51uA03jwBglto8wEL3PNR0A0pjUe6jS2m67RzdqP10/qEU83P3f9nH2XcPWpB6XFBbx42ace1abLBNAFCyJs/bsJ8R0H36kvY5C5uffHD086DtHN24/XTUtuPsu+NGZpssaHHQPKY+FtkAgMP0tRAVZRKAjtCXtM+qSV+u3kHB1rj9tI59t6kbNQdVWo6IeP7DDz332Od91nx8/9e/ym8ZgN4zwEJJBKBUxt21apW+JFATN2rGFVH4gXd9MD717Gfue+z/3flMrW0DgJIZYKEUAtAOqKq0+qSvM+557q5VS8W6B40rljCq+FXfvysAgBIJQPdp2zpZVY2STfo6hz3P3bWjGbXfqVj3oGmq/0b0+7uCiPad0wDoPlVw99gNroabW5HjXnC1ujFsumljHTRKNovXqer9uGfcfrfwWaOXr9mdU7u6MYzllbU4efFaLK+sFb2fVmVclerPO+S7gj5q4zkNgO4TgO7RxuCqqlGySV/HqFz1xu13OcfYJYH6emE5bgmW7//6V1k+CfZp4zkNgO6TgrtHG4OrqirPTvo6bah027aUs3H71x9v3Ym3vPGxkZ9leWWtt/NDD0rzblO/w6y18ZwGQPcJQPdoQ3C1azfIGm5uRYqIvOexo4z8TFrBtvRKt6VXjh3loP1uXLDlwvJB08w/bttNCjiKNp3TAOgPKbh7jJtfVkpwtWtv+mXEdvCZdh7bTUmc9mJ6XGrj/teZ9HlNKSXlbJr5mUfZ78ZdQLqwPFxf05fpn7ac0wDol2ONgKaUviUifiAi/npEvCbnvF5Fo5rSlmVERgVZObaDwacuPn7k1510BKnkSrcljAxOOwp7lP2u9JHoklnehr5oyzkNgH45bgruByLifET8hwraUoSSg6tdJQRZpSoh5ewoAc60+50Ly6Pz+6FPKdhtOKcB0C/HCkBzzk9HRKSUDnsqFSohyCpVCSODdQU4LiyPxu+n39o4TxwAuqS2OaAppTellNZTSuu3b9+u6207ybye8UqYo2p+Ztn8fvqtlHniANBXh46AppR+KSK+aMRDb845/9ykb5RzfiIinoiIWFpayoc8nQNIvzxY0yODJYzCMp7fT79JwQaAZh0agOacv7qOhjCdpoOsqk06J6sNc7cEOOXr2u+HyUnBBoBmWQeUxk06J6tNc7cEOFAmGQoA0KxjzQFNKX1TSuljEfHlEXEtpXS9mmbRJ5POyWrD3K1p1v+EUvRpvy1hnjgA9Nlxq+C+MyLeWVFbWqENKaBtM+mcrNLnbv3L1Zvxtvf8TuxOcC55hBZ29XG/laEAAM2prQpuF+ymgA43tyLHvQu1Lo8W1GHSqrElV5dd3RjedxG/q7QRWtjLfgsA1E0AOoU2pIC20aTLYpS8fMaV67ceuIjfVcoILexnvwUA6qYI0RRKTwFtq0mrxpZcXfagfaCEEVoYxX4LQFuZFtdeAtApKNql3qgAAAlXSURBVN8/O5POySp17ta4fSNFFDFCC6PYbwFoozatjMCDpOBOoeQUUJo1at9IEfEdX/YyB0KKZb8F6Le2VkE3La7djIBOoeQU0NL0LS3CvkEb2W8B+qvNo4imxbVbynlcCYrZWVpayuvr67W/L/XYf0CL2B4pttYeAEAZllfWRk7DWFwYxFMXH2+gRZNrc9v7JKV0I+e8tH+7FFwqJy0CAKBsbR5FNC2u3QSgVK7NBzQAgD4oeX31w5w7sxiXz5+OxYVBpNge+ZRp1x7mgFK5w6oF921+KABAaS6cPTVyylRbRhFLXRmBwxkBpXIHpUXszg8dbm5FjnsT3ttSdQ0AoAuMItIUI6BU7qDKmssra2PnhzrgAQDUxygiTRCAMhPjDmjmhwIAQH9JwaVWbZ7wDgAAHI8AlFopmw0AAP0lBZdaHTQ/FAAA6DYBKLUz4R0AAPpJCi4AAAC1MAJamNWNYWvTU9vcdgAAYPYEoAVZ3RjGpas3n1snc7i5FZeu3oyIKD6Qa3PbAQCAekjBrdDqxjCWV9bi5MVrsbyyFqsbw6n++yvXbz0XwO3aunM3rly/VWUzZ6LNbQcAAOphBLQiVYwAPrO5NdX2krS57QAAQD2MgFakihHAEwuDqbaXpM1tBwAA6iEArUgVI4AXzp6KwfzcfdsG83Nx4eypY7WtDm1uOwB0yXGnBAHMkhTcipxYGMRwRLA5zQjgbqpuGyvJtrntANAVigICpUs559rfdGlpKa+vr9f+vrO0/4AfsT0CePn8aQd8AKAWyytrI2+ILy4M4qmLjzfQIqCvUko3cs5L+7cbAa2IEUAAoGmKAgKlE4BW6NyZRQEnANCYKqYEAcySIkQzpAgAAFAnRQGB0hkBnRFFAACAupkSBJROADojB60L6iQAAMyKKUFAyaTgzogiAAAAAPcTgM7IuMn+igAAAAB9JQCdEUUAAAAA7mcO6IwoAgAAAHA/AegMKQIAAABwjxRcAAAAaiEABQAAoBYCUAAAAGohAAUAAKAWAlAAAABqIQAFAACgFgJQAAAAaiEABQAAoBYCUAAAAGohAAUAAKAWKedc/5umdDsifrv2Nx7vRRHxB003ggfol3LpmzLplzLplzLplzLplzLplzKV3i9/Kef8yP6NjQSgpUkpreecl5puB/fTL+XSN2XSL2XSL2XSL2XSL2XSL2Vqa79IwQUAAKAWAlAAAABqIQDd9kTTDWAk/VIufVMm/VIm/VIm/VIm/VIm/VKmVvaLOaAAAADUwggoAAAAteh1AJpSupJS+nBK6TdTSu9MKS3seexSSukjKaVbKaWzTbazb1JK35JS+mBK6TMppaV9j+mXBqWUXrfz3X8kpXSx6fb0WUrpx1NKn0gpfWDPts9PKb07pfS/dv7/85psY9+klF6aUnoypfShnWPYd+9s1y8NSim9IKX0ayml9+/0yw/ubD+ZUnrvzvHsZ1JKz2u6rX2UUppLKW2klH5+52/9UoCU0v9JKd1MKb0vpbS+s82xrGEppYWU0s/uxC9Pp5S+vI390usANCLeHRFfnHP+koj4nxFxKSIipfTKiPi2iHhVRLwuIv59SmmusVb2zwci4nxE/MrejfqlWTvf9b+LiK+NiFdGxLfv9AnN+InY/h3sdTEifjnn/PKI+OWdv6nPsxHxfTnnV0bEl0XEP975jeiXZn0qIh7POb86Ih6LiNellL4sIv51RLwl5/xXI+KPIuK7Gmxjn313RDy952/9Uo7X5pwf27PMh2NZ8340Iv5bzvkVEfHq2P7ttK5feh2A5px/Mef87M6f74mIl+z8+xsj4qdzzp/KOX80Ij4SEa9poo19lHN+Oud8a8RD+qVZr4mIj+Scfyvn/OmI+OnY7hMakHP+lYj4w32bvzEifnLn3z8ZEedqbVTP5Zw/nnP+jZ1//2lsXxgshn5pVN72Zzt/zu/8L0fE4xHxszvb9UsDUkoviYjXR8SP7fydQr+UzLGsQSmlF0bE34qIt0ZE5Jw/nXPejBb2S68D0H3+YUT8ws6/FyPid/c89rGdbTRLvzTL91++L8w5f3zn378XEV/YZGP6LKX0aESciYj3hn5p3E6a5/si4hOxnf30vyNic89NaMezZvxIRPzziPjMzt9fEPqlFDkifjGldCOl9KadbY5lzToZEbcj4j/upK3/WErps6OF/fJw0w2YtZTSL0XEF4146M0555/bec6bYzt16m11tq3PJukX4OhyzjmlpMx5A1JKnxMR74iI78k5/8n2oM42/dKMnPPdiHhsp9bDOyPiFQ03qfdSSm+IiE/knG+klL6y6fbwgK/IOQ9TSn8xIt6dUvrw3gcdyxrxcET8jYj4pznn96aUfjT2pdu2pV86H4DmnL/6oMdTSn8/It4QEV+V761JM4yIl+552kt2tlGRw/plDP3SLN9/+X4/pfTinPPHU0ovju3RHmqUUpqP7eDzbTnnqzub9Ushcs6bKaUnI+LLI2IhpfTwzmib41n9liPiG1JKXxcRL4iIvxDb89v0SwFyzsOd//9ESumdsT0Nx7GsWR+LiI/lnN+78/fPxnYA2rp+6XUKbkrpdbGd+vENOedP7nnoXRHxbSml56eUTkbEyyPi15poI/fRL8369Yh4+U6FwufFdkGodzXcJu73roj4zp1/f2dEyCao0c78tbdGxNM55x/e85B+aVBK6ZGdkc9IKQ0i4mtie37ukxHxd3aepl9qlnO+lHN+Sc750dg+n6zlnL8j9EvjUkqfnVL63N1/R8Tfju0CkY5lDco5/15E/G5K6dTOpq+KiA9FC/sl3Rv065+U0kci4vkR8X93Nr0n5/yPdh57c2zPC302ttOofmH0q1C1lNI3RcS/jYhHImIzIt6Xcz6785h+adDOneofiYi5iPjxnPMPNdyk3kop/VREfGVEvCgifj8ivj8iViPi7RHxsoj47Yj41pzz/kJFzEhK6Ssi4r9HxM24N6ftX8T2PFD90pCU0pfEdmGOudi+8f72nPO/Sin95dgupvb5EbEREX835/yp5lraXzspuP8s5/wG/dK8nT54586fD0fEf845/1BK6QvCsaxRKaXHYrto1/Mi4rci4h/EznEtWtQvvQ5AAQAAqE+vU3ABAACojwAUAACAWghAAQAAqIUAFAAAgFoIQAEAAKiFABQAAIBaCEABAACohQAUAACAWvx/KgC4o2ajza8AAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 1152x432 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "mhtNrqtBOiHP"
      },
      "source": [
        "def det_loss(y, y_pred, mu, log_var):    \n",
        "    reconstruction_error = (0.5 * (y - y_pred)**2).sum()\n",
        "    kl_divergence = (-0.5 * torch.sum(1 + log_var - mu**2 - log_var.exp()))\n",
        "\n",
        "    return (reconstruction_error + kl_divergence).sum()"
      ],
      "execution_count": 5,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7vaY25CeSEs9"
      },
      "source": [
        "def ll_gaussian(y, mu, log_var):\n",
        "    sigma = torch.exp(0.5 * log_var)\n",
        "    return -0.5 * torch.log(2 * np.pi * sigma**2) - (1 / (2 * sigma**2))* (y-mu)**2"
      ],
      "execution_count": 6,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7HuL6T_ROlgP"
      },
      "source": [
        "class LinearVariational(nn.Module):\n",
        "    \"\"\"\n",
        "    Mean field approximation of nn.Linear\n",
        "    \"\"\"\n",
        "    def __init__(self, in_features, out_features, parent, n_batches, bias=True):\n",
        "        super().__init__()\n",
        "        self.in_features = in_features\n",
        "        self.out_features = out_features\n",
        "        self.include_bias = bias        \n",
        "        self.parent = parent\n",
        "        self.n_batches = n_batches\n",
        "        \n",
        "        if getattr(parent, 'accumulated_kl_div', None) is None:\n",
        "            parent.accumulated_kl_div = 0\n",
        "            \n",
        "        # Initialize the variational parameters.\n",
        "        # 𝑄(𝑤)=N(𝜇_𝜃,𝜎2_𝜃)\n",
        "        # Do some random initialization with 𝜎=0.001\n",
        "        self.w_mu = nn.Parameter(\n",
        "            torch.FloatTensor(in_features, out_features).normal_(mean=0, std=0.001)\n",
        "        )\n",
        "        # proxy for variance\n",
        "        # log(1 + exp(ρ))◦ eps\n",
        "        self.w_p = nn.Parameter(\n",
        "            torch.FloatTensor(in_features, out_features).normal_(mean=-2.5, std=0.001)\n",
        "        )\n",
        "        if self.include_bias:\n",
        "            self.b_mu = nn.Parameter(\n",
        "                torch.zeros(out_features)\n",
        "            )\n",
        "            # proxy for variance\n",
        "            self.b_p = nn.Parameter(\n",
        "                torch.zeros(out_features)\n",
        "            )\n",
        "        \n",
        "    def reparameterize(self, mu, p):\n",
        "        sigma = torch.log(1 + torch.exp(p)) \n",
        "        eps = torch.randn_like(sigma)\n",
        "        return mu + (eps * sigma)\n",
        "    \n",
        "    def kl_divergence(self, z, mu_theta, p_theta, prior_sd=1):\n",
        "        #log_prior = dist.Normal(0, prior_sd).log_prob(z) \n",
        "        #log_p_q = dist.Normal(mu_theta, torch.log(1 + torch.exp(p_theta))).log_prob(z) \n",
        "        log_prior = ll_gaussian(z,torch.zeros(1),torch.log(1 + torch.exp( torch.ones(1)) ))\n",
        "        log_p_q = ll_gaussian(z,mu_theta, torch.log(1 + torch.exp(p_theta)))\n",
        "        return (log_p_q - log_prior).sum() / self.n_batches\n",
        "\n",
        "    def forward(self, x):\n",
        "        w = self.reparameterize(self.w_mu, self.w_p)\n",
        "        \n",
        "        if self.include_bias:\n",
        "            b = self.reparameterize(self.b_mu, self.b_p)\n",
        "        else:\n",
        "            b = 0\n",
        "            \n",
        "        z = x @ w + b\n",
        "        \n",
        "        self.parent.accumulated_kl_div += self.kl_divergence(w, \n",
        "                                                             self.w_mu,\n",
        "                                                             self.w_p, \n",
        "                                                             )\n",
        "        if self.include_bias:\n",
        "            self.parent.accumulated_kl_div += self.kl_divergence(b, \n",
        "                                                                 self.b_mu, \n",
        "                                                                 self.b_p,\n",
        "                                                                 )\n",
        "        return z"
      ],
      "execution_count": 9,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "w6424H66OpD8"
      },
      "source": [
        "@dataclass\n",
        "class KL:\n",
        "    accumulated_kl_div = 0\n",
        "\n",
        "class Model(nn.Module):\n",
        "    def __init__(self, in_size, hidden_size, out_size, n_batches):\n",
        "        super().__init__()\n",
        "        self.kl_loss = KL\n",
        "        \n",
        "        self.layers = nn.Sequential(\n",
        "            LinearVariational(in_size, hidden_size, self.kl_loss, n_batches),\n",
        "            nn.ReLU(),\n",
        "            LinearVariational(hidden_size, hidden_size, self.kl_loss, n_batches),\n",
        "            nn.ReLU(),\n",
        "            LinearVariational(hidden_size, out_size, self.kl_loss, n_batches)\n",
        "        )\n",
        "    \n",
        "    @property\n",
        "    def accumulated_kl_div(self):\n",
        "        return self.kl_loss.accumulated_kl_div\n",
        "    \n",
        "    def reset_kl_div(self):\n",
        "        self.kl_loss.accumulated_kl_div = 0\n",
        "            \n",
        "    def forward(self, x):\n",
        "        return self.layers(x)"
      ],
      "execution_count": 7,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "cA_bqFFdOuLc"
      },
      "source": [
        "epochs = 2000\n",
        "\n",
        "def det_loss(y, y_pred, model):\n",
        "    batch_size = y.shape[0]\n",
        "    reconstruction_error = -ll_gaussian(y,y_pred, torch.ones(1)).sum()\n",
        "    #reconstruction_error = -dist.Normal(y_pred, .1).log_prob(y).sum()\n",
        "    kl = model.accumulated_kl_div\n",
        "    model.reset_kl_div()\n",
        "    return reconstruction_error + kl\n",
        "\n",
        "m = Model(1, 20, 1, n_batches=1)\n",
        "optim = torch.optim.Adam(m.parameters(), lr=0.01)\n",
        "\n",
        "X = torch.tensor(x, dtype=torch.float)\n",
        "Y = torch.tensor(y, dtype=torch.float)\n",
        "\n",
        "for epoch in range(epochs):\n",
        "    optim.zero_grad()\n",
        "    y_pred = m(X)\n",
        "    loss = det_loss(y_pred, Y, m)\n",
        "    loss.backward()\n",
        "    optim.step()"
      ],
      "execution_count": 10,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "gkIgq66oOu-8",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 391
        },
        "outputId": "933a9dd7-ab48-4061-e704-7a0d08af2665"
      },
      "source": [
        "with torch.no_grad():\n",
        "    trace = np.array([m(X).flatten().numpy() for _ in range(1000)]).T\n",
        "q_25, q_95 = np.quantile(trace, [0.05, 0.95], axis=1)\n",
        "plt.figure(figsize=(16, 6))\n",
        "plt.plot(X, trace.mean(1))\n",
        "plt.scatter(X, Y)\n",
        "plt.fill_between(X.flatten(), q_25, q_95, alpha=0.2)"
      ],
      "execution_count": 11,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.collections.PolyCollection at 0x7f4f40e83e50>"
            ]
          },
          "metadata": {},
          "execution_count": 11
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6AAAAFlCAYAAAD1f4WfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzde3Bc6Xnf+d/b5/QNlwbIIXgBCJBDUeJYmpFFGxlrLFn2jGVRcaw1oyS1G2cnJavWUy6vUs5mi7ua8jqOk9qSEpa965KcWitrRc6sdq11ZUx7faNkDx1H9sgOx5RFjTSQZyTPkCCHxIjEjcSlL+/+cfo0uoFuoIE+fS7d309VF4AG2HiBbpL96+d9n8dYawUAAAAAQLelol4AAAAAAKA/EEABAAAAAKEggAIAAAAAQkEABQAAAACEggAKAAAAAAgFARQAAAAAEAo3im964MABe/z48Si+NQAAAACgy1544YU3rLVjm6+PJIAeP35cly9fjuJbAwAAAAC6zBjzarPr2YILAAAAAAgFARQAAAAAEAoCKAAAAAAgFARQAAAAAEAoCKAAAAAAgFAQQAEAAAAAoSCAAgAAAABCQQAFAAAAAISCAAoAAAAACIUb9QIAAAAAoBdcuDKr8xdndGN+ReOjeZ07c0pnT09EvaxYIYACAAAAQIcuXJnV089e1UqxLEmanV/R089elSRCaB224AIAAABAh85fnKmFT99KsazzF2ciWlE8EUABAAAAoEM35ld2dX2/IoACAAAAQIfGR/O7ur5fEUABAAAAoEPnzpxSPu00XJdPOzp35lREK4onmhABAAAAQIf8RkN0wd0eARQAAAAAAnD29ASBcwdswQUAAAAAhIIACgAAAAAIBQEUAAAAABAKAigAAAAAIBQEUAAAAABAKAigAAAAAIBQEEABAAAAAKEggAIAAAAAQkEABQAAAACEggAKAAAAAAgFARQAAAAAEIqOA6gxJmeM+QtjzF8ZY140xvx8EAsDAAAAAPQWN4DbWJP0hLV22RiTlvRFY8zvW2u/FMBtAwAAAAB6RMcB1FprJS1XP0xXL7bT2wUAAAAA9JZAzoAaYxxjzJcl3Zb0BWvtnwdxuwAAAACA3hFIALXWlq2175B0VNKjxpiHN3+NMeYpY8xlY8zlubm5IL4tAAAAACBBAu2Ca62dl3RJ0vubfO5T1tppa+302NhYkN8WAAAAAJAAQXTBHTPGjFbfz0v6IUkvdXq7AAAAAIDeEkQX3COSfs0Y48gLtP+vtfZ3ArhdAAAAAEAPCaIL7lcknQ5gLQAAAACAHhboGVAAAAAAAFohgAIAAAAAQkEABQAAAACEggAKAAAAAAgFARQAAAAAEAoCKAAAAAAgFARQAAAAAEAoCKAAAAAAgFAQQAEAAAAAoSCAAgAAAABCQQAFAAAAAITCjXoBAAAAAHrXhSuzOn9xRjfmVzQ+mte5M6d09vRE1MtCRAigAAAAALriwpVZPf3sVa0Uy5Kk2fkVPf3sVUkihPYpAigAAACArjh/caYWPn0rxbLOX5zpuwBKJdhDAAUAAADQFTfmV3Z1fa+iEryBJkQAAAAAumJ8NL+r63vVdpXgfkMABQAAANAV586cUj7tNFyXTzs6d+ZURCuKBpXgDQRQAAAAAF1x9vSEPvbBRzQxmpeRNDGa18c++EjfbTulEryBM6AAAAAAuubs6Ym+C5ybnTtzquEMqNSflWCJAAoAAAAAXeUHcLrgEkABAAAAoOuoBHs4AwoAAAAACAUBFAAAAAAQCrbgAgAAIFAXrsxy1g1AUwRQAAAABObCldmGbp+z8yt6+tmrkkQIBcAWXAAAAATn/MWZhlETkrRSLOv8xZmIVgQgTgigAAAACMyN+ZVdXQ+gv7AFFwAAIEK9dl5yfDSv2SZhc3w0H8FqAMQNFVAAAICI+OclZ+dXZLVxXvLCldmol7Zn586cUj7tNFyXTzs6d+ZURCsCECcEUAAAgIj04nnJs6cn9LEPPqKJ0byMpInRvD72wUcSXdUFEBy24AIAAESkV89Lnj09QeAE0BQVUAAAgIi0OhfJeUkAvYoACgAAEBHOSwLoN2zBBQAAiIi/TbWXuuACwHYIoAAAABHivCT6Wa+NIcLOCKAAAAAAum5z2Hz8oTH9xxdma52g/TFEkgihPYwzoAAAAAC6qtnM289+6bWeG0OEnRFAAQAAAHRVs5m3tsXXJn0MEbZHAAUAAADQNReuzGp2F6GSMUS9jTOgAAAAQMB221ynV5vx+FtvWzFqrIQyhqj3dVwBNcZMGmMuGWO+Zox50Rjz00EsDAAAAEiiZucdn372qi5cmQ3k65Ok2dZbXz7t6B+9c0oTo3kZSROjeX3sg4/0RPBGa0FUQEuS/kdr7V8aY4YlvWCM+YK19msB3DYAAAB6VK9W/ZqFLr+5TrOfb7dfnyTbneckbPanjiug1tqb1tq/rL6/JOnrkngkAQAAoKVervq1Cl1BXZ8krc5zTozmCZ99KtAmRMaY45JOS/rzJp97yhhz2RhzeW5uLshvCwAAgITZruqXdK1CV1DXJ8m5M6eUTzsN13HOs78FFkCNMUOS/qOkf2qtXdz8eWvtp6y109ba6bGxsaC+LQAAABKol6t+uw1dvRzSzp6e0Mc++AjnPFETSBdcY0xaXvj8rLX22SBuEwAAAL1rfDTfdDRHL1T9/HDV7vnW3X590pw9PdEzPws6Z6xtNQK2zRswxkj6NUl3rLX/tJ0/Mz09bS9fvtzR9wUAAEBy+WdA67fh5tMO1TGgRxhjXrDWTm++PogtuO+S9KSkJ4wxX65efjiA2wUAAECPYmsm0J863oJrrf2ivBmyAAAAQNvYmgn0n0DOgAIAAABAJ3p1LiwaEUABAAAARGrzmWB/LqwkQmiPIYACAAAAiNR2c2GTGEC7Us21ViqueJfciOQkM8olc9UAAAAAekYvzYUNpJpbqUjF+1Jp1XvrB09VJ5hkHiKAAgAAAOgvQVX6emku7K6rueViY8gsrkjltZBWGz4CKAAAAAJDI5n+EeS5zXNnTjWdC3vuzKngFhyS1tXcTSHTr25WSiGvMFoEUAAAAASCRjL9Jchzm/7X98KLF+Ojeb0+v6y81pRTUTmzrpzWNDmckuZeinp5kSOAAgAAIBC91kgG2wv63GYi58JWKl4ls7Raq2r+q3dKv/xH17VW2vi7kHUd/ePHHoxwofFBAAUAAEAgeqmRDHbWS+c2d2RtQ8j03l9telbzibfsk7En9czzr2puaVVjwzk9+dgxPX5qLIKFxw8BFAAAAIHoq0CCrpzbjPwMsR80/YDpv19aU60DbRsePzVG4GyBAAoAAIBA9FIjGews6HOboZ4h3hI0V7yQucugid0jgAIAACAQvdRIBu0J8txmV84QB1TRRHAIoAAAAAhMIhvJIBY6OkNsbbWCuULQjDkCKAAAAIDItXWGuFLxGv/44bLob51dFUEzGQigAAAAACJXf4bYUVlZFTWSLutn33NY+vYrXtBs0nUWyUIABQAAABCdaiXz7ElHA2dG9H/9ydd1d+nexviS4ylpbTHqVSIgBFAAAAAA3edvlW04o7kq2UrtS953Iqf3nTgd4SLRbQRQAAAAAMEpF+vOZq5sjDmpC5roXwRQAAAQC5EPoAewO5Vy86BZKUW9MsQYARQAAEQu1AH0AHanXPKC5eaus5Vixzd9aWZOzzz/quaWVjfOfJ4aC2DRiCsCKAAAiFxXBtAjtqh2x1TD1tm6M5pdqmhempnTJ597WWsl7+/+7aVVffK5lyWJENrDCKAAACByHQ2gR6JQ7Y6BSkUq3q82BFrZeGvLO//ZAD3z/Ku18OlbK5X1zPOvEkB7GAEUAABErq0B9OgJVLtDVlqvns+su8Rklubc0uqurkdvIIACAIDI1Q+g9+XTjs6dORXhqtANVLu7xNqNSqZ/iXlDoLHhnG43CZtjw7kIVoOwEEABAEDk/MpXL5wL5Hzj9qh2B6DWDGi1ca6mbNQr25UnHzvWcAZUkrKuoycfOxbhqtBtBFAAABALZ09PJD6ocb5xZ1S7d6FSqTunWVfZDPmsZrf45zzpgttfCKAAAAAB4Xzjznqp2h2ocslrDFQfNktrSlpVc7cePzVG4OwzBFAAAICAcL6xPb1Q7e5IubgRNv235fWoVwWEggAKAAAQEM43YovS+tawWSlGvSogMgRQAACAgHC+sY9VKtWts6uNMzZ75LwmEBQCKAAAQEA439gniqsbYdN/G5PZmkDcEUABAAAC1PfnG3tJrarpb5+tBk5biXplQGIRQAEAAIDSeuOok+IKVU2gCwigAAAAIbpwZbbvtujG6me2duN8pn8prUiVUjTrQeJdmpnr2izTirW6fa+ia4slXVso69piWdcWSyqlX9T//mPTgXyPsBFAAQAAQnLhymxDk6LZ+RU9/exVSerZEBrpz1xa98JmfeAsrarXZ2siPJdm5vTJ517WWsl7fN9eWtUnn3tZktoKodZa3Vmp6PpSuRowS9WQWdb1xZJmF8ta37Tj+9BgSicOrslaK2NM4D9TtxFAAQAAQnL+4kxDh1xJWimWdf7iTM8G0K7/zOWiVFrzLuW1auBc997nrCa67JnnX62FT99aqaxnnn+1FkCX1iq1QFkfLr23Zd0rNr4gsi+X0mTB0VsPpPW+EzkdLbiaLDiaHHE1Mewo5xpp7CEpgeFTIoACAACE5kaTGaHbXd8LAvmZK+VqyFytC5rVC2NOEKG5pVVVrFFRrtblal1pFeXqtUVXH/j1OV1bLGl+tTFgDqaNFygLrr53Mqujw04tYB4ddjScTUX004SDAAoAAEIRq3OAERkfzWu2SfAaH81HsJpwtP0z+2czG6qZ1UulGNJqga1KFaubSxvnL/3K5bXFkmZ0TOtyGr7eyCqfKmtfLqW3H8prsuDo6LCryREvaO7LpRK5dTYoBFAAANB1/Xj2sZlzZ041/B4kKZ92dO7MqQhX1V2NP7NVRiWNpCv6mR+YlOavSeXqOc3yetRLRZ+qWKs5v9HP4kbQvF59/+ZSWeW6ImbKSEeGvDD56Lirb968q1TFq31mVNKgK/2THzypx089EN0PFWOBBFBjzKcl/Yik29bah4O4TQAA0Dv68exjM/7P2vOV4HKpVs08+yajgTMj+g9/MqP5pWUdHM56XUKnytL9N6JeKfqAtVbzq7YuYFbfLpR1fck7j7m+aSf32IB3DvO7j2Q0ecqpbZk9WnB0ZMhR2tmoYF6acatdcNcC74Lbi4y1nXcBM8a8R9KypP/QTgCdnp62ly9f7vj7AgCAZHjwo7/btO+okfStj/+dsJeDIBVXpPX71Rmaq4w0QSTurVcaqpe1gFmtZC5vavQzmjPettiCU90a64XLyYKjowXXa/QTZ2MPSel4b903xrxgrd0yKyaQCqi19k+MMceDuC0AANB7+vHsY0+pbwLkb5mtfnzppVtdm4EI+FZLVrNLmzrJLmy8f3e1sePxQLXRz9FhR+88mq2FSz9oFnq80U+chXYG1BjzlKSnJGlqaiqsbwsAAGKgH88+Jo4/usR/W17feL9JRfPSzJw+9Z9e0dLaxud2OwMR8JUqVjeXG7fFenMxvYrmrXuNATOTkiaq1cqHD6arAbNa0Sw42p/v70Y/cRZaALXWfkrSpyRvC25Y3xcAAESvb84+xlmlstFZtry+9W3TTdLNXZqZ0yefe3nL/ENp6wxEQPLOYc7dr26TXSjp+lK5GjC9KubN5bJKdRnTb/RztODo+6ayDVtkJwuuDg2llCJgJhJdcAEAQCjOnp4gcIahXKo7j1k3NzPALrPPPP9q0/Dpm1taDex7IRmstVpYqzb6WdjaSfb6Yklrmx4yB6qNfk4fzui/qjt/OVlwdGTYUcYhYPYiAigAAEAS+XMziyvexX8/hJmZOwXMseFc19eA8N1br2ypXPrh8vpiWUvrjVX0QtZosuDqzftdPXHcq2JOjnjnMo8WHOXTnMPsR0GNYfl/JP2ApAPGmOuSfs5a+6tB3DYAAEDfK5ek4v2tgXMX22aDNDac0+0WITTrOnrysWMhrwhBWCtZ3Vj2tshunod5fbGsb680nsPMu6a2LfbR8axXwazrKDtCox80EVQX3H8YxO0AAID+cOHKLOdBm2la1bwfu7EmTz52rOkZ0ELO1U+8502c/4ypcrXRz/XFrVtkry14jX7qX9JIp6SJarXyfSdyXtgc2Wj08wCNfrAHbMEFAAChunBltqEj7uz8ip5+9qok9VcIjVlVczuXZua2jFr5yBMnGb8SM/WNfq7Xh8vqucwbmxr9GEn7clKl6J0RPpaRfvAto3rfqVGv0c9gSk6KgIlgEUABAECozl+caRjHIkkrxbLOX5zp3QBaWvfCZnFl420IZzWDsLnjrT9q5SNPnNSnP7Rlxjy6bGG1suX8pVfB9MaXrJYaX8A4kE/paMHRdx7O6EcKTkMn2b++eUf/xx+/orVK2UujRenLM3f07omTGp/gxQR0BwEUAACE6sb8yq6uT5RKpdp5tq6qWbwv2dYdY+OuWcdbRq10z/1iZUvl0mv84wXNzY1+hjNGkwVHb9rv6geON44rOVpwNLBNo59/eeE17luEjgAKAABCNT6a12yTsDk+mo9gNXtUf1azFjhXvXEnPaZVx1tGrezNetnqxlJdF9mFxkY/b2xq9JNzjY4OO5occTQ9nvE6ydaNKxnJ7b3RD/ctokAABQAAoTp35lTDGVBJyqcdnTtzKsJVbaO0vjFX0z+zWVpTHM9qdkOrjreMWmmuXLG6da9um2y1cnm9WtF8fbmx0Y+bksaHvID53hO5WuXSH1lyoIuNfrhvEQUCKAAACJV/zjNWXXArFa96WVr1Amd5zQuZxZVEb58NQrOOt/08auW5l27r0392XTeXy8rl83ro6AFlMrla4LyxVFZxU6Ofw0MpHS24euxodiNcFrzQeXjQiazRD/ctomCsDf/Vu+npaXv58uXQvy8AAOhTlYpUXq9eil7A9N8vrSWmIVBUmnXB7eUzggtrlerZy5KuL2ycx3xpbk2v36uoosZtr8MZ6cS+dGO4rI4sGR9ylHXj20m23+7bnjH2kJSO97EFY8wL1totncqogAIAgOQrVcNlpVgNmMW6sLkeWMDs1yfrj58a66mfc6Wu0Y+/NfZaXeOfxbWtjX6OFhytr61on9aUUUlplZRRUWmVdCSX1af/62R2BO61+xbxRwAFAADJYW1jd1m/EVAI22RbjSORxBP4mCnWGv1sNPi5Vhc037jf2Ogn66jW1Oe7Duc1OeLWGv9MFlyNZI2MMfrAJ74h26SYSdMeoH0E0GZK61LK8S4AACAa5VLdOJP7Gx1nI2r+wziS+KhYq1vLlS2VS7+T7M3lsip1DxPHSOPD3rbYHzyeqwbLjU6yBwZSSrXR6IemPYiavwvjzxb/Rg+MjkR/fn4PCKDNrNyVlm5I6QEpNyrlRyU3G/WqAADoPf44E7/5j99htrQau+Y/jKwIj7VWd1YqDeHy2mJZ16vnMmcXy1pvLGLq0GBKkwVXj45nNDmyES4nC44ODzlyA2j0Q9MeRKl+F4aVNDu/oqefvSpJiQqhBNDtFO97Fz+M5vd5gdTNRL0yAACSxa9m+uHSf7/F3Mw4nrWk+hWsxbVKw/nL6/UjS5bKul9srHTvz6V0tODorQfSet+b8huNfgquxocd5UJo9OM/BuP22ER/aLYLY6VY1vmLMwTQnuSH0cVZKT0o5Ua8zlPpAcnh1wgAgKSNmZm1oFl9Wym1fRNxPWtJ9Wt3VktW1xcbZ2Beq+sou7Cp0c9Q2ujoiKNjo67eNZWthUt/LuZQJtXiO4WLpj2ISqvdFjfmV0JeSWdITntRvOddfKm0F0bdnPc2Oyw56ejWBwBAN1XKdRXMzdtmKzv/+R3E9awl1a9GxbLVzWW/k2zdWcxqFXNuU6OfTF2jn3cczm8aWeJqNOc1+gHQXKtdGOOj8R7HshkBNAiVorRWlNYWq1cY79zo4JiUGYx0aQAA7Eml4gXK8lpdNXPd+3gX1cy9iPNZy36qflWs1e17lYbK5fWljfdfXy6rvKnRz5Fqo5/Hj+dqlUs/YI4NttfoB4iTOB0HaLYLI592dO7MqUjWs1cE0K6wXiOjlbvedt386EZ1lMooACBO/MY/myuZ5fXIlsRZy3BYa3V3tdroZ6HxLOb1xZKuL5W1vqkP1MHBlCYLjv7WeKYWLP2QeWQ4mEY/QFzE7ThA/S6MlxelidE8XXDRxJbtuq7k5qV0zjs/6ua8Syoe5xoAAD2mUvbCZLlYvWyqZgawZTZonLUMzvJ6pWFbrL9d1p+LeW9To5/RnNFkwdVDB9L6oRO5jU6yI64mQmr0A8RFHI8D1HZhjD3kFbcSiAAatkpJWl/yLjXGG/PiV0nTeS+k0m0XAHrChSuzOn9xRjfmVzTejVesrW1s+lO/bbbL22W7gbOW7VstWc0u1XWP3TQXc361MWAOpk1ta+w7j2Y3OsmOuDo67Gg42/4L4t3Ymhin7Y5AnI8DJBkBNBbqZqCtzm9cnXK9KqkfStMDzCMFgIS5cGVWTz97VStF71X0Pc1ts9arYlbKXqBs2C5bbQQku+PNJEk/nbXcTqlidXOp3NhJtu7t7XubGv2kpKMFbwbmIwfzG51kR7yguS+XCqTRTze2JsZtuyMQ3XEA4+WA2sXZ+rGT3EIVATTOKiWvsVGtuZEk4zQGUjruAkCsnb84UwufvpZz2+q7y9bPyiytqdcCJjwVazXnN/qphcyNgHlzqbHRT8pIR4a8MPn9x7I6OuxvkfWC5sGQGv10Y2tiHLc7or8FcxygjTDZ8LHb80fzCKBJY8vS+rJ38aUHpGxhYxSMm5XoMgcAsbB1PptVRiUtzt+Xluc2zcosRrJGdI+1VvOrtqFyuREyvS2zmxv9jA14jX6++0hGk6ccHR32tshOFhwdGXKUdrr/f/xOW2G7sTWR7Y6Im5bHAR465I1hrFUj/ffTkuPWvZ/2Pt+Brh/hiAABtBcU73uXmk1nSt2sd6bUyfT8KyoAEAuVSu0c5sMjq7qzsKSsKSqjklx5ZzIPDuekxesRL7QznNfz3PMb/WwJmCVdXyhreVOjn5Gs1+jn1H5X730w29BJ9mjBjbzRTztbYbuxNTFu3Y95fPcTU61EVoOks1GNfPzRKT3+zkcbQ2ZIz6cDOcIRQwTQntTiTKnkbeF10huv0Lh5bxtvZiCapQJAUvmNf8prW0eZ1I0w+WfvHNYnn7vVcUfXuD0Z7oXzeu3+Ttf8Rj91DX6uL3hjSq4tlHVntfEcZt41teY+75zINoTLyYKjwi4a/UShna2w3ehUHKfux73w+O57W7a1Vi9+VbIWJquBM4Z2dYQjQeL520b32LJUKkvyX2G8Ky3J+8uXGfLCaG6Ec6UAIFXPZG7qLlsueu+3uV02iI6ucXwynPTzevW/U2ul64slfewPr+k/z1Y0NJCvzsX0tsjeuldpOIGbTkkTw96ZyzMn0xuNfqpBc38+mEY/UWlnK2w3OhXHqftx0h/fPSlVt7V187lJp0nQTPDfQd/WIxzbX58UBFB4KiWvWro6Ly1c88JobkTKjTIOBkDvsrY6G3O9blZm8CNMNnd0vTQzpw9/5nLbT7Lj+GQ4Sef1rLWau19pGFPy2cvzWiwe1LpcFeVKMlJZ+sZXy0qZZR0ZcjQx7OhdU9mGcDlZcHVoKPhGP3GqcLe7FbYbnYrj0v04SY/vqHT0mE253q68lk14mgTMPjQ+mtdsk7A5PprM+Z++/rw3sTO/0dHirJQe9MJofpQxMABCF0gDhnJRKq40r2aG3F12L9XMOD4ZjtN5PWutFtZs4/nLhVJDo5+1TY1+XGWVVkkDWlNay8qoVL0U9Ts/9T3KhNDoxxe3CnectsJGJU6P72aifsHCf8yuliqqyNG1pYr+zXPXVUwP630PT9RVIR3CZAfOnTnVcAZUkvJpR+fOnIpwVZ3jEYCdFe95l6UbXsddvzJKUyMAXbarBgz+rMzNQbO06h0/iIm9VDPj+GQ47JByv1ht9LNQ1vWlUm2LrB86l9YbX0goVBv9vHm/qyeOZ2vnL72zmI7++8/+ZdPf6cHhXKjhU4pfhTtOW2GjEucQ3r0XLLYbF9JYrfz5X7ut14qTqqjueWBR+vkvruh9757qYA2o5/8/Rxdc9De/4+7SzeoVpvEfJpPy3ro5Kb+P7bsAOtKsAcNqsahf+oOv6OxDA1JxtTFoJmBW5l6qmXF8Mhx0SFkve41+/C2yXsDcmI357ZXGRj+5ukY/j45na8HyaMHV5IijkR0a/cTpdxrHCndctsJGJc4hvO0XLEyqcVRIrQFPfais/7j9osKrCyVZbf36pJ9NjKOzpycSHzg3I4CiQ9ZrxNGsGcfSjWoYrY6CaRgJw0MPSLowZpO9Mb+ggorKal15s66c1pVRUWZR0t3DgX6vsOy2mulvtVsrleUYo7K1OhiTJ8O7CSnlitWte5WNTrILG+Hy2mJJry83Nvpxa41+HL3vRM7rJDviVmdiOjqQT+mPv/GGnnn+W/rLV1d1rfo7efhge+uJU8CIY4UbMQ3hxtHsUllF5VSSo5J1VJSjkhy9tuhKB95SFyq7t0utV88mIhykAHRXpSStL3mXen7X3eEjUpr/YIGkCXw2WWldKq3UVTS9y7sKb/TcE/PdVN42b7UrW1v72rg9MbbW6s5K4zxMb1yJ9/7sUlnFuiKmkXRoMKXJEVePHfVHlVS3yY44OjzoyEm13gobxDbEuASMOFVjETKzaZtrQ4Vy8+gQr7Pr8sgbTcPfxGheygyGsuxePZuIcBBAEY1a190FafCANHSYqiiQIHueTebPy/QvxVVdevE1PfNn32paherFJ+a7qbzF7Wzg8nqlaYMf//37xcYt0PtzKU2OOHrbwbTefzJfm4c5WXA1Mewo6+79rGXcfjediFM1Fh2q3/ZaGw2yaRvspkC5W3EIf716NhHh4Bk/Imale3PS/TteNXTgARobAQmw42wy/0zm5jOam5oB7VTF6tUn5u1W3sI+G7hW8s5hXqsLldfrwubd1cZzmINpUwuU3zuZbZiHebTgaHCmMLAAACAASURBVCiz+3/P2+3u2awyvt31cReXaixaMV7zRTfrhUcn4wXLLSGz+89h4hL+evFsIsJBAEU82LK0eN0b+5Id3rikOUsAxIq1UqWkB0dSmltYlquSHFkZVZQ1RU0Np6SbfyXZys63pfaqWP38xDzos4HlitXN5XJD9bJ+bMmte43nMDMpaaLa2Ofhg+mNgDniVTL35VIyAc7D3M22Wv887GZODwyfRxhMi/mTTTrAOmnvEiN++PPP4v8Pn/uyzl+coQrZBWH0O+g3BFDEjJXWFr2L5L2iWB9G3dyetqsA2IVKuW6MSV0Fs1LyLpJ+9p2OPvnc3JatsT/x2Mm2w6cUz+6f7QpjDt9utyBba/XGSqXWQfZ6/XnMhbJuLJdV2nQO88iQV61811R24wxmNWgeGkopFeK/ubvZVtssfG53PTZEPUOyK/wu/O0Eyh6ZRRn4WXxswe+4O5L/tw+9rVKUVu54F8n7DyY94AVRN+eNeXFz3lYYginQHn9eZmlNKhe998tr3vnM8lotZG4nqK2xwzlXi6tbv1/cmwx1bw5fo2a/5783PaWD+0f1B6+sVBv8NI4rWSk1BrAH8ikdLTh6+6G0/s6b87Xq5WTB1fiwE/rMy+3s5gWJgy2qwwdj/tiJWliP3Y5tFx6bXvrv+M6ez+KjbfyOu4MAimSxFWl92bs0qDub4WYlJ+uF0/ozGkA/KRerAdMPmmsb1zUbm7QHnW6NvTQzp/vr5S3XuykT+yZD3WyAs1qytcY+XgUzo+zYCdlsWS8slvTccyVJc7WvH84YHS04enDU1XumstWA6XozMYcdDe7hHGZUdrPluBcbVIUh0uZNtSplelPHV/8sZfVspZPmReU27HgWHx3jd9wdPCtHj7DVJ9hr0lqTT/td6fxzHFve5z88JEilUq1artdVMOverxR3tQ02Ks88/6pKla3bJfNpZ8sT4bhtGexk63CpYnVzqdwwA7O+inn7XuN9l3Gko9Wtse84nG8YVTJZcDWSNYGew4zSbkJlrzao6raubHuvVSHrg2R6U8Cshk4Ehlmc3cfvuDsIoOgPtrIRULeTcqtVU3cjlNb/R+pk+A8U3dcsVDYEzJ23yAapW+Gv1RPe5bXGny+OWwa3q9RZazV3f+u4Ej9g3lwqq1yXu1PGO4c5WXD0/ce2nsMcGwz3HGaUdhsqe6FBVdgvrrRdZTabZlLW/38YwCgRdC4O41h6Hb/j7iCAAvXqmqy0VKumZjZtGaqGVslr4iLrfa2/HRiQvMfGduGyXJQUnyYq3Qx/7T4RjuO8xw9OT+mX/+SalkspFeVqXa5KJqNvl/N66N/e1NqmncUHBlKaLDj6rsMZTb6lGi5HvKB5ZMhROkbnMKPWC6GyXaG+uFL9v+vH3v2Qzn/hFd0rSSWbUkmu0um0fvK9b5cOTBIqEyQu41h6Gb/j7ggkgBpj3i/plyQ5kv5Pa+3Hg7hdIJbarabWM6lq06Tsprd09e0p1u5cvbRbzzzGWTfDX7vbLaPolFt/DvO1hbptstWK5tK6lXSk9vWuKhofSumtB3PVGZgbVcyjBUf5dHLOYSI8nf398seIpDe6utbPo3Tcxq2x1SY97z0kLQ+f0PmLM7pV94T6R3hCnUjM4uw+fsfB6ziAGmMcSb8s6YckXZf0X4wxv22t/Vqntw30DFuRive9y2ZOtnnzJCfbl139Yqlc2qiOV0re/em/74fLAJv7xEk3w1+72y2DnoUpSRVrdfteRdcWS3ptoazX6rbKvladh1kv66g2/3J6PFPbJnu0et1Ilr+r2L3Nf4+sUirK0atLVsqNNAbKZk179ogn1ACiFEQF9FFJL1trvylJxphfl/SjkgigQDu2a56Ucr0tm/Xzzepf1TZOXeipq6Qa4/2Zhkuz65p8Ta+xthoYy95bW95o0GOt93Glen252Bg0/Y9jtCU2bN0If/Xa2W65126nS2uVhhmYr1XDpn9dfQNefx7m5Iij75vKamrE1VTdNtmxgVTPNPpByJpWJb2P1wrHdW1hXSW5KiklK++FjInRvLT/RMQLB4DuCCKATki6VvfxdUnfE8DtAvDPo9qyVC57QbWrRbZqSHUyUmZAcvPV7sD14bT61v/Y2mpgK1fDXN37tdus3q7/vv99au/X3ZYf9ra8L+9j/3pb2Rog/eukuo+72w02bt1ZgxaHURetKqXf9+YDurZQ0muLGxXMWiVzoaw7q433/XDGaGrE1Vv2u/rBB71mP1MjjqZGXI0POcq6BEy0y9Q15anv8rr54+3PUv7U+79LTz97VUUanAChunBllnOdEQqtCZEx5ilJT0nS1NRUWN8WwK5UA11pxbtgW3Hszhq0nbbJdjOAW2s1v2r12mJJy2ZI3/22U7q+6FUy//nzZd34/M2GbrJuSpoYdjRVcPX+k2lNFrxwOTXiXTeS63ybbK+/4AA1jg1pOraregkADU6A8F24MtvQ2XZ2fkVPP3tVkvi7FxJjbWdby4wxj0n6F9baM9WPn5Yka+3HWv2Z6elpe/ny5Y6+b1ct3ZKWbkS9CgAx9+HPXG66PfXgcE6f/tB0BCsK1+YALnnV0Y88cbLtULZasppdqlYuq9tka5XMhbKWi43/Rx3Ip2rzL/1gebTgaGrE0eEhR25q71XMncJlED8vIlY/YsvJ1AXKurFbnL0Hetq7Pv5c09meE6N5/elHn4hgRb3LGPOCtXbLE6IgKqD/RdKbjTEPSpqV9N9I+rEAbhcAYi2K7qxx0k4Hz4q1mrtX8YLlQrm2XfZ69e3rTZr9+Ocvv2ciW6ti+h1lBzPdCQftVLPjOA4GdbbMb86ocYZz++GS7XlA77rRJHxudz2C13EAtdaWjDEfkXRR3hiWT1trX+x4ZQAQc91u0BN3ftAuW6Oi0lqXq3WldXPR1Y//9rd1rTq+ZG1Ts5/DQylNFly9u9rspz5kRtXsp51w2e8vOETLNFYqnczW91NOIN+J7XlIKl44ac/4aL5pBXR8NB/BavpTIGdArbW/J+n3grgtAEiKODToCYO1VndXK/qbea9q+epCWa8ulPRa6qjulb1B9vUcVXT7XkUn97t6otrsxw+ZE8PxbPbTTrjsxxccQjvzmnI3ZiQ72bpxI25bZy6DfOJ9/uJMLXz6Voplnb84w5N5xBYvnLTv3JlTDb8rieZfYQutCREA9Jp251gmQcVavb5c0d/Ml2oh87WFkv6m+nZpfeMspj+y5Eghq7n5Rbm2qIy8y5Br9dNPPKgnHgr2d9DtINROuOyXFxx8gTXZMo62doitq1662Y6ql0E/8WZ7HpKIF07aR/Ov6BFAAaAD7cyxjIv1stX1Ra96+apfyZz33m6ei5lOSUcLjo6NuJo+MqCpEUfHR10dG3F0tOAqV61iXppxux7Aw+g23E647KUXHNqx7bbkhw5vVChT9QEzvVG19GdfdrmpT9BPvNmehyTihZPdOXt6gsAZIQIoAPQQW61kfnO+pG/eLW28vVvS7FJZlbqmsgNpo6kRR2/e7+q9D3rnMeeXlvXFl25oafmeDlZyevKhY5LW9czzr+pCk9AVRgAPo/lPu+EySS847IpJ1Y0f8YLli4t5lTSoohyV5KqslEpy9OJiSjry9qhXXBP0E2+25yGJeOEESUIABYAEWl6v6Ft14fKb82V9825J35ov6X7d6JKBtNGDo67ecTijv/uQV9E8PuqNMDmQb2z4c2lmTs/+1Te9sGe8SuP/9vkZ1fepjWLWaVjNf5qFy0TP/TSprecpt1Qt3ZbdYd3Rcd1qMaogToJ+4s32PCQRL5wgSQigABBT5Yq3ZbZZNfNW3fgSI2+77Il9rh6dGNCJUVcn9rk6Merq8FD7XWWbVRorTb4u7NEjUTX/CWPrb0c2n6UM8FyllJwntN1YJ9vzkDS8cIIkIYACQMTurnhbZl+5W9K36oLmq/MlrdclwJGs0Yl93viSE6Ou3rTP1YP7XB0b2TiT2YndVBTDHD0SVfOfyOZ+tmza426cs3Qy0i7G1eylS2xSntAmZZ1At/HCCZKCAAogUZK6JXKtZPXaQkmvVMPlxvbZsu6ubqTMdEqaGvEqmE88mNObRl09uM/RiVFX+/PdnZHZqtLY6mvDElXzn8C3/vrB0sk0NuqpD5YttsN2opMusUl5QpuUdQIACKAAEiTuWyKttbp9r6JX6qqYftC8ttjYAGhsIKUT+1y9/6QXMk/s8y5HC47c1N5DZicBvVmlsRkjafr4qD78mcuhBcIomv/suPXXOHWzKjefsdx0cdK7qlgGifEM/SXImagA0A0EUACJEdmWSDUGuwNDOf3w6SntHynopTeK+us7G9tnl+saAOVcowdHHT18MK0fPZWvnc08PuqqkA1+NEWnAX1zpXEo62qlWFap0jgD9O1HR/RHX5+L7QsBgUi5+rF3n9K//sK3tFQ0KspVUa6cdFY/ceY7pSPHIguUu8V4hv4R9ExUoNft9IINL+h0BwEUQGJ0qxvqdlXDuysV/fqX5/R/vzCn5cqQVrVfLy5ldOlPSpLuSJLGhxy9ab+rv//WjQZAD+5zdGTIUSrEkBJEQN9caWz2u4nyhYBAORmvS6y7+a3XwOe9h6Xl4Tfp/MUZvT6/IscYlYtW/+YL35RNpRPzJKQfxzP065NGqt1A+3Z6wYYXdLqHAAogMbrRDdWvGq4UK1pTRjOLrv7nL3xbY5dLur1iNHffP5+5X47Kympdo1pSTus6MiD9+ycf6Uo1cy+6EdCbbX39xc/PBP59usNsdIR1s40h08m0ddbSf5KR5CchSelmG5R+ftJItRto304v2PCCTvcQQAEkRlDdUFeKFX3tjZK+cmtdv/L8gu4Wj2hNaXkbTCVTqWjl7pred2pEDz2Q1me/OKOs1uWq3LDrcn1FsQmfUnjjSqIai9KUP+uyvnrph85ddoptJelPQvqtS2zS769O9GO1G9irnV6w4QWd7iGAAkiMvXRD9cPm1dvrunq7qKu3inr5bqnWEMhVRjmtq6B7ymldOa0ro6JSVvqFHzouSXr+K1a3l7Y25okkcG0jrHEloY9FMamNOZfpvOTm6kJmujvfs04vPAnppy6xvXB/7VW/VbuBTuz0gg0v6HQPARTAFnEedbJdN1Q/bH71dlFfub2ur972GgT5YfNAPqVHDqX1/pM5PXIwo0cOpvX0b1zR3PL21byo5lDuVljjSrryfZyslBnYCJr+uBInI6WcgFa+NzwJSZZ+vr/6rdoNdGKnF2x4Qad7jLV2568K2PT0tL18+XLo37cdF67M6lf/4EvS4o3YPfEGwrC5k6rkha2PPHEyVn8XVktWX5srelXNurBZrgubDx9M6+2H0np4LK1HDmZ0eGjrHM12f976UD6UdWWMtLRaSsS/E7F6QcGkJDcvpXONb534vh66+Uyh5D0J+dgHH+GJfQxxfwFoF11wu8sY84K1dnrL9QTQDf5/WkPFb+uQ8bpbxvGJN9BNH/7M5abn+w4O5/TpD235NyQUftj86u2irs6t6+qtxrD5QD6lRw6m9cjBtBc6W4TNVnYT0JIS0H2RrtfJNgbN9IC3bTaBeBKSLNxfABC9VgE0vi85R8BvXDBUd10iRwsAHQiik2onFbfVktXX3yjWttFevV3UX3+7MWw+fDCtHzqR08MHvcrmkV2EzWa229a7WdJGkISy3lTaO5e5uarZRpfZpOinM5S9gPsLAOKLAFqnVYOC+I0WCEastuWhqSjuo047nG6uuN1eWtUnn3tZkras/X6xoq+/UWrYSvuNurC5P+ed2Xzvg7laZbPTsNmpbs0i7ZbA1mucTeNM/Esu8jOaAAAgOQigdVo1Lohbp8sg7CYkIBpR3UedNtxpVXH71T+9Jjc/rBfnivraXEkvzhX1rfmNBkH7c15l84njOT1yyNtOOz7kRBo2m4nVCJI27H69RsoMSpmhxsAZ4zOaAAAgOXhGUcfvdqXixnVx7HS5W82qaEnbRtiPorqPOu1wentxVWtKa1VprSqrVWW0oqxKy66+eME7Wz0x7OitY2n9yJtzettYWm8dS2tiOH5hs5mkdMT1tbVeJytlh6TciJQZ7qmts9gZ5yUB+Pj3AGEggNbx/4J5XXDVE9tSW1XRNgcbX1y3EfajKLd6tnMmcrVk9bkrc/rcl+c0tyLJzSmdzeumHpSVHyStsipqUCsay1b0L3/4hL7jQFr78l7AuTQzp2f+8zcStQ08rFEnQdm83v3DA/pvv+8hPf6dJ7yxJ+lBqps9rJ0Oj/UdY2fnV7wXYiWedAJ9hn8PEBaedWxy9vSEzp58t7R0I+qlBKJVFc0xRuUmHZDjuo2wH8Vlq+fCWkWv3Cnp5bslvXynqJfvlPTynZKuLZblPYL2SbLKlIrKl+/r0bG0bt1ZkFNeVVZFpYz1uq7+wEl97+RGB9SothgHca52N02LImMc73ymm9Xj0+N6/J2PSul8YrvQYvfaeTLpN9+rt1Is6/zFGZ5wAn2Gfw8QFgJoj2tVLStbLxQkZRthP+r2Vs+Ktbq7UtHr9yq6fa+s15fLunWvolvLZd2qfnxjuaz51Y0XKjKOdGKfq7cfSstZu6vS2oqyWldGRaWqRc/BtZx+7r3HqiHPtgx5UWwx7tmzzynXG3GSHqhWNQckJx31qhCxdp5Mtmq+1+p6AL2Lfw8QFgJoj2tVRRvOenf9Wsn7uJBz9RPveVOyn4T3oIxj9nQfLa9XtoTJ+o9vVUNnsbL1z6ZNWY4tadC1esfhAT12bFhv2ufq5H5XkwVHTjVpfuATX5NtcmRzbmm1rQphFFuMe+Lss0lVw2a+GjgHqWqiqXaeTLZqvjc+mu/auhCNfj3b168/917w7wHCQgDtIc22Fjarorkpo5ViWaXKRmVrrbR1Oy6i80cv3dYnnvum7pesKsqorJRW1l09f6Oi62v3tLhW0eK61dJaRQtrFc2vVnR3paL5Vau7qxWtNLk/h9JGB4ccHR5M6dHxjA4NOTo0mNLhIUcHBx1969a8Pvtnr6hYLktGUllauuXoLW87qcdPbA1nnW4RjmKLcdJGqEjG20brn9VM571LApo1IXrtPJn0m+/VV0rzaUfnzpwKZY0IR7+e7evXn3uv+PcAYSGA9ohWWws/8sRJfeSJkw3BdLVY0uJqqeHPx7UKlKRZpdZarZWle+sVLa9bLRetltcrurdutbRe0b2i9a6vfv6e/37RanZhXTcX17VekSpyVJaRNNX4DcrSK1fLkhYkedthC5mUClmj0VxKR4YcfcdYSvtyKR0YSOnQoONdqkFzKLN9Z9NfvviaFz7rbPe46HSLcBTdZONyrrYlJ+sFzMzgxpZaOtJij9p5Mll/FpQKUe/q17N9/fpz7xX/HiAsBNAesd3Wwk9/aLohQHzgE19sehtxqwKFcV5vT6GxuDVA+n+m2ZbWZobSRoMZo6FMSpVyWW8s3VfKVjSoihxVlNr01n/fVUWf/e+mNZxJKecGWwVrdf/fXlrVBz7xxS0vAHTaDbaTP7/XFyZiNULFP7eZ8SubdKNFsNp9Mnn29ARPMHtcv57t69efuxP8e4Aw8Gynid//6k39xqXLiai6+XaztbCbVaAgK5Y7nderWC8cLq5VtLhmtbBW2fK+FxKDCY2DaaOhjNFgJqXhjNFg2mhqxNVQpnp9OlV9P1X9OqPhTKr25/zrB9JGqbotlB/+zGUNalVqI08eHM5pbMCR1Px37f/e9vL7b/W4kCSr5i8AdNoNdi9/vpMXJiIboVI7tzmwsZ3WzXT3ewLiySQ8/Xq2r19/biDuCKCbXLgyq1/4vZc0UvKeiCelS+ZuQmW3qkB7CQa2GiIX1ipaWLNaWPXONC6sVvS1xaxKyqssR2WlqhdH31hM6Tt/5aaW1q0qOxxdHayrNA5Vg+DkiFsLkH4oHMqkql9nNFQLktXrq7eR6tK5u3Yrz/X3UbPf9S/94TckqXa2d7eP3WaPi83qXwCIanv0di9M+J/fbk3dH6FiNhoE+YHTzXFuE0Bk+vVsX7/+3K3QkAlxQQDd5PzFGa0Xyxqpe64Y1/OR9XYTKoOuAllrtVqy+nd/ek0LxZTKSm8ExqKjn790R8/dSHsBs9owZ3GtooVV7+NyyxD5gCQrR+XaNlRXZQ2myzpzakQj2ZRGsikVst45yEL1/ZHq+0MZU+vYGmfbdSrOZ9ym91GzEFZqksZ389jd/LhodbfMLa1GOs5ku63CkawplZayQ15VMzMguXnObQKIlX4929evP3czNGRCnBhrw+9+Oj09bS9fvhz6923Hgx/9XR3QvA6ZOw3XG0n/3z95dzSLalOnFSlrrZbWvSrkvN9ZdbX6sX9Zq2h+xf98davrakXr225jtRrNpTSaTWkktxESR3NegBypfq6QM7VQOZpL6a9e+7b+3X96RevlxlD9kSdOxvrFgN3aHOaknX/OD3ziiy0D4mZ7fex++DOXmwbjg9WqeqvPffpD07v+XkGsyzFG5Sb/ngW+JpOSMkNStiDlCoxAAQDE3rs+/lzT7cgTo3n96UefiGBF6AfGmBestVuehFEB3WR8NK/1+fkt18emS+Y2/K2Flbog+ZVb67pbDY8Lq7YWLBsuaxshs3U10tvOOppLVYOj0akHnFqAHMkaffb5b8nYckPF0lFZaSP99lPv2vXPM/62g8q5JjFdcPdqLxXp7c5rNvvavdiuqv6Ln59p+mfCaGTVal2ttg4HsyYj5Uak/D4veFLhBAAkCA2ZECcE0E3OnTmlX3j2tlQ3pSSyLpnyqpL3ilZ3Vyq6s1rRnZVK7f27K97H9e/fXfUu252NHEobjea9quNozujIsKPRnDe+w69SekHTaF/OD5gpZZztt7P+7vOLTRvptNnjp6nun9eLh93+nK3mu0qNW3E7eexuF4yfef7VSMeZZByjterf0ULO1U+8503dWVNmyAud+X1Sytn77QAAECEaMiFOCKCbnD09oezqQ/qNS0tdqbqtlmxjgKyFx7Lurm6EyG9Xg+bdldbbWx0j7cuntD+X0r58Sm95wNW+fDVI+pes97n6La/pHYLkXh1sUZU7mIDqcVy12lbdKhw2u64bXWqjGmfSbLvyWskL3NPHR/X7V19v2Jq8pzU5WS9wDuxney0AoCfQkAlxwhnQZpZuSUs3dvyyUsV65yRXqoGxLkBuqVRW394rtv59j+ZMLUzuyznan/fef6AaKvfXv82nVMgYGWMC60ba7u20+rq9nGVEa3H/fUbRBbfV+c/hrKv1sm34XRlJf/uRw/qpHzjZxi2bauh8wGsoBABAj6ELLsLGGdA2fe4vXtOvXLys8r27yueyevvUA9o3lN+yzfXOijc2pJXBtNH+fKp2ObnfbQiPftD0w+VILlXbQrkbQXUjbfd22vm6Xj+z2Q3NwtxOc1CjFsX26FbnOZfWSluus5Iu/83W89xb5PdJw0eodgIAehpzgREXBNA6F67M6md/60Wtl7OSDkur0kvfqMg193RgoBoc8ylNFDK1AFl7W1edHM2llHPDGf8RVEhp93Z2+rp+ObMZpFahfqemOlHN4YzSbpovSds0IHKyXgfbgQNSmi3i6F1UPAAAcUMArXP+4ozWyxWNaFkHNF/r4np4KKt//+PdHS2xV62eYO+282e7txPU98OGVqG+1ViRseFcpHM4o9Tq7GnWNVpc3VoFbWxAZKT8aHWb7XAIqwWixdw/AEAcEUDr+K2oMyopb9Zr17+xHL9w5Ve/Wm0C3m3nz1aVpc230+7XRSlplcFW4b1s7ZbxIn5Tnbhvz+2W7ZovtWyK5Oa90JnfJzn8k4f+cf7iTEPDEUlaKZZ1/uIMARQAEBmejdVp1aI6TuFKat6cpt5eOn+229U0qu6n7YpLZXA3IbhVqD9YdxZ08+1EOYczattt8/Z/Vw8M5/UPv//tevxvvU3KDG75OrYloh8w9w8AEEcdBVBjzD+Q9C8kfYekR621MW5tuzO/RbWKG9fFKVz5mlW/fAf3WPFrt4FQ3BsNxaEyuNsQvF2obxW2klCJDtvjp8b0+FuPSkMHpfx+KZVq+nVsS0S/YO4fACCOOq2AflXSByX9SgBriZz/5PNX/+BL0qJiF658rapcRtKnP7T3s6rtNhCKc6OhOJxR3W0I3kuoj3slOnQmJQ0d8i5m+wZgbEtEv1TAmfsHAIijjgKotfbrkmR2eMKXJGdPT+jsyXe3NQc0KlS/WovD72YvIXi3oT7ulehQZQvSyNG2x6iwLbG/9VMF3P95+iFsAwCSI7QzoMaYpyQ9JUlTU1NhfdueRPWrtTj8bsIKwXGuRIcilfaCZ350V3+MbYn9rd8q4Mz9AwDEzY4B1Bjzh5ION/nUz1hrf6vdb2St/ZSkT0nS9PR0q+ataEMvVr/abdqz09fF4XcThxDcu4w3QiW/T8qNtjznuR22JfY3KuAAAERrxwBqrX1vGAvB7vRS9avdpj3tfl3Uv5s4hOCeYlJe6MyNSrkRKeV0dHNsS+xvVMABAIgWY1gQuXab9sShw+1ONldo/9n7TsVmbYliUt7ZzvyolB3ZU6VzO2xLbNQvTXkkKuAAAESt0zEsf1fSJySNSfpdY8yXrbVnAllZTO1mviPa027Tnjh0uN3Ov/3jl/X7V1+Xv788qhmkiWUcKVfwKp3ZQuChE839Lxeu6rNfeq32uO3lpjwSFXAAAKLWaRfc35T0mwGtJfZ2O98R7Wm3aU8cOty2cmlmriF8+uJWoY2dWqVzn7e9toc6aifBhSuzDeHT18tNeSQq4AAARIkSwy5stwUUe/fkY8eUdRvP9TVr2tPu10Xhmedf3fIk3heXCm2sZIalkSnp0MPS/ge9rbaEz9CdvzjT8nFLUx4AQJxduDKrd338OT340d/Vuz7+nC5cmY16SWgTZ0B3Ie5bQJOq3aY9cW7us91jIA4V2lhwc1J+v1ftdDNRrwbaPmTSlAcAEFf9vdNkoQAADypJREFUNNO5FxFAdyHOW0CTrt3OtVF3uG2l1WPDSLGo0EYm5XqBM79fygxEvRps0qojrJFoygMAfSCpTej6baZzr2EL7i7EeQto3FyamdOHP3NZH/jEF/Xhz1zWpZm5qJfUVc0eG0bS337kcCwDc1eZlBc695/wttiOHCV8xtS5M6eUT2993P6jd07xHzgA9Di/ijg7vyKrjSpiErayMtM52aiA7kKct4DGST82a+KxobpmQqN0sE0IOsICQP9KchWRmc7JRgDdpbhuAY2TJMzr7Ia+fGykXGlwTBp4QHLSUa8Ge0BHWADoT0muIjLTOdkIoAjcTs2amKXaA/zgOTgmpZydvx4AAMRKkquI7OBJNgIoArdds6Z+3J7bUwieAAD0hKRXEdnBk1wc1ELgtmvWxCzVhEq50vAR6eBbpeHDhE8AABLu7OkJfeyDj2hiNC8jaWI0r4998BFCHbqOCigCt11Dnl/8/EzTP8Ms1ZhKpasVzwOETgAAegxVRESBAIquaNWQh1mqSWCk7LDXWCg3IhkT9YIAAADQIwigCNWTjx1rOAMqMUs1dClXcjJe11on0/x9AAAAoAsIoAgV8zK7yKS8cGmcaqBsEjBTaWZ0AgAAIDIE0JhJ8oiSdtfel/Myg2ZSUn6fNHDAC5cpl62yAAAAiD0CaIwkeURJkteeKE7Wawg08ABNgQAAAJA4BNAAdVq93G5ESdxDXJLXHnsmJWWGvOCZG4l6NQAAAMCeEUADEkQFsNUokiSMKEny2mMpPSjlClK2IKXzbK8FAABATyCABiSICmCSR5Qkee2xkEp7o0/8C51oAQB7dOHKrM5fnNGN+RWNj+Z17swpZj0CiA3aYQYkiArgk48dU9ZtPNeXlBElSV57NIyUGZaGx6Wxh6TDD0v7jkkD+wmfAIA9u3BlVk8/e1Wz8yuykmbnV/T0s1d14cps1EsDAElUQAMTRAUwySNKkrz20DiZaoWz4L2liRAAIGDnL85opdi4I2ulWNb5izNUQQHEAgE0IE8+dqzhDKi0twpgkkeUJHnt3WG85kHZYe88Zzof9YIAAD3uxvzKrq4HgLARQAPSrAI4fXxUzzz/qn7x8zNUBPtJbkTK76fKCQAI3fhoXrNNwub4KC+CAogHAmiA6iuAzMXsI+kBr9KZGSR0AgAide7MKT397NWGbbj5tKNzZ05FuCoA2EAA7RLmYvYwAicAIKb8c550wQUQVwTQLmEuZg/xA2d2yHtL4AQAxNjZ0xMETgCxRQDtEuZiJhiBEwAAAOgKAmiXBNUVFyEgcAIAAAChIIB2CXMxY4zACQAAAESCANpFzMWMCQInAAAAEAsEUPQmJyMVxqVsgcAJAAAAxAQBFL1n4AGpMEHwBAAAAGKGAIre4WSkkUkpV4h6JQAAAACaIICiN+T3SyNHqXoCAAAAMUYARbKl0tLopJQbiXolAAAAAHZAAEVy5fd5W26pegIAAACJQABFwhgpOywNHqDqCQAAACQMARTxl3K9WZ65US90OjxsAQAAgCTimTwiYqR0XjKOlEp5b01KctLeuU4n7XW1dTLe5wEAAAAkHgEU4UmlvREp2WEpW+DsJgAAANBnOgqgxpjzkj4gaV3SK5J+3Fo7H8TC0AuMlBnyAmeu4FU8AQAAAPStTvc2fkHSw9bat0v6hqSnO18SEs3NSUOHpH3HpcOPSAdOSsOHCJ8AAAAAOquAWms/X/fhlyT9/c6Wg0SqbxCUzkW9GgAAAAAxFeQZ0A9L+lyAt4c4Sw9I+f1e6HQzUa8GAAAAQALsGECNMX8o6XCTT/2Mtfa3ql/zM5JKkj67ze08JekpSZqamtrTYhEDTkYaPiIN7I96JQAAAAASZscAaq1973afN8Z8SNKPSPpBa63d5nY+JelTkjQ9Pd3y6xBTxpGGD0uDY5IxUa8GAAAAQAJ12gX3/ZL+J0nfb629H8ySEC/GC53DhxmbAgAAAKAjnZ4B/aSkrKQvGK8q9iVr7U92vCrEQ36fNDzOGU8AAAAAgei0C+7JoBaCGMkMS4VxKTMQ9UoAAAAA9JAgu+Ai6bIFaeiglB2OeiUAAAAAehABtO8ZKT8qDR2S0vmoFwMAAACghxFA+1HKlTKDUmZIyo1yxhMAAABAKAigfcFIuRFvi21mUErnol4QAAAAgD5EAO1lJiXl93vnOt1s1KsBAAAA0OcIoL3IONLgAW9+p5OOejUAAAAAIIkA2mOMFzqHD0spJ+rFAAAAAEADAmivyI16szvZagsAAAAgpgigSZce9IJndijqlQAAAADAtgigSeXmpeFDUn5f1CsBAAAAgLYQQJOG4AkAAAAgoQigSZEekIYOSfnRqFcCAAAAAHtCAI07gicAAACAHkEAjav0gDdOJTcS9UoAAAAAIBAE0LgwKa+jbXbYu2QGol4RAAAAAASKABql9MBG4EwPSqlU1CsCAAAAgK4hgIbJyVYD55CULUgpJ+oVAQAAAEBoCKDdlhn2RqZkhyU3E/VqAAAAACAyBNBuMCkvdA6OSel81KsBAAAAgFgggAZtcEwaOiw5/GoBAAAAoB4pKSjpAWlkku61AAAAANACAbRTxpGGj0hDY1GvBAAAAABijQDaifw+qTAhOemoVwIAAAAAsUcA3YtUWhqdknKFqFcCAAAAAIlBAN2tbEEaPUaTIQAAAADYJVJU24y33ZazngAAAACwJwTQdrg5r+pJh1sAAAAA2DMC6E4GHpAKR6VUKuqVAAAAAECiEUBbMSlvrufA/qhXAgAAAAA9gQDajJOWHngzW24BAAAAIEAE0GaoegIAAABA4DjYCAAAAAAIBQEUAAAAABAKAigAAAAAIBQEUAAAAABAKAigAAAAAIBQEEABAAAAAKEggAIAAAAAQkEABQAAAACEggAKAAAAAAhFRwHUGPOvjDFfMcZ82RjzeWPMeFALAwAAAAD0lk4roOettW+31r5D0u9I+ucBrAkAAAAA0IM6CqDW2sW6Dwcl2c6WAwAAAADoVW6nN2CM+V8l/WNJC5Ie3+brnpL0lCRNTU11+m0BAAAAAAljrN2+aGmM+UNJh5t86mestb9V93VPS8pZa39up286PT1tL1++vNu1AgAAAAASwBjzgrV2esv1OwXQXXyDKUm/Z619uI2vnZP0aiDfOBgHJL0R9SKwBfdLfHHfxBP3Szxxv8QT90s8cb/EE/dLPMX9fjlmrR3bfGVHW3CNMW+21v519cMflfRSO3+u2UKiZIy53CydI1rcL/HFfRNP3C/xxP0ST9wv8cT9Ek/cL/GU1Pul0zOgHzfGnJJUkVfR/MnOlwQAAAAA6EUdBVBr7d8LaiEAAAAAgN7W6RzQXvGpqBeAprhf4ov7Jp64X+KJ+yWeuF/iifslnrhf4imR90tgTYgAAAAAANgOFVAAAAAAQCj6OoAaY84bY14yxnzFGPObxpjRus89bYx52RgzY4w5E+U6+40x5h8YY140xlSMMdObPsf9EiFjzPurv/uXjTEfjXo9/cwY82ljzG1jzFfrrttvjPmCMeavq2/3RbnGfmOMmTTGXDLGfK36b9hPV6/nfomQMSZnjPkLY8xfVe+Xn69e/6Ax5s+r/559zhiTiXqt/ej/b+9eQq2sAiiO/xfXR2EvMjPRQiNJLNQaiJKEaQ8r0YoIo8AeEEFEQRGZgyBw0KSUqJFWDiwRy5LA0EywiVamYHQrzAoVH1FJRaAYq8G3w6PZg7Kz7+msH8jdj4NszuLse/f59rc/ST2Stkp6q9STSx8g6StJ2yVtk/RhactcVpmksyStLOuXXkmTOzGXrl6AAuuAS22PAz4H5gFIGgvMAS4BZgAvSOqpNsru8zFwC7CxtTG51FXe6+eB64GxwO0lk6jjZZrPQavHgfW2RwPrSz3a5wjwiO2xwCTggfIZSS51HQKm2R4PTABmSJoEPA08a/si4Hvg3opj7GYPAb0t9eTSd1xle0LLYz4yl9W3CHjb9hhgPM1np+Ny6eoFqO21to+U6iZgRCnPBpbbPmT7S2AHMLHGGLuR7V7bn52gK7nUNRHYYXun7cPAcppMogLbG4HvjmueDSwt5aXATW0dVJezvdf2R6X8I80fBsNJLlW58VOp9i//DEwDVpb25FKBpBHAjcDiUhfJpS/LXFaRpDOBK4ElALYP2z5IB+bS1QvQ49wDrCnl4cCulr7dpS3qSi515f3v+4ba3lvK+4ChNQfTzSSNBC4DNpNcqivbPLcBB2h2P30BHGz5EjrzWR0LgcdonicPMJjk0lcYWCtpi6T7SlvmsrpGAd8AL5Vt64slDaIDc/lXzwHtBJLeAc47Qdd822+W18yn2Tq1rJ1j62Z/J5eI+OdsW1KOOa9A0mnAa8DDtn9oLuo0kksdtn8BJpSzHlYBYyoPqetJmgkcsL1F0tTa44nfmWJ7j6RzgXWSPm3tzFxWRT/gcuBB25slLeK47badksv/fgFq++o/65d0FzATmO6jz6TZA5zf8rIRpS1Okr/K5Q8kl7ry/vd9+yUNs71X0jCaqz3RRpL60yw+l9l+vTQnlz7C9kFJG4DJwFmS+pWrbZnP2u8KYJakG4BTgDNo7m9LLn2A7T3l5wFJq2huw8lcVtduYLftzaW+kmYB2nG5dPUWXEkzaLZ+zLL9c0vXamCOpIGSRgGjgfdrjDGOkVzq+gAYXU4oHEBzINTqymOKY60G5pbyXCC7Cdqo3L+2BOi1/UxLV3KpSNKQcuUTSacC19Dcn7sBuLW8LLm0me15tkfYHknz++Rd23eQXKqTNEjS6b+VgWtpDojMXFaR7X3ALkkXl6bpwCd0YC46etGv+0jaAQwEvi1Nm2zfX/rm09wXeoRmG9WaE/8vcbJJuhl4DhgCHAS22b6u9CWXiso31QuBHuBF2wsqD6lrSXoVmAqcA+wHngTeAFYAFwBfA7fZPv6goviPSJoCvAds5+g9bU/Q3AeaXCqRNI7mYI4emi/eV9h+StKFNIepnQ1sBe60fajeSLtX2YL7qO2ZyaW+ksGqUu0HvGJ7gaTBZC6rStIEmkO7BgA7gbsp8xodlEtXL0AjIiIiIiKifbp6C25ERERERES0TxagERERERER0RZZgEZERERERERbZAEaERERERERbZEFaERERERERLRFFqARERERERHRFlmARkRERERERFtkARoRERERERFt8SvmYA2l9mFXTwAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 1152x432 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    }
  ]
}