GaussianProcess.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/bezhvin/Optimization-with-Tensorflow/blob/main/GaussianProcess.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "8911bb01",
      "metadata": {
        "id": "8911bb01"
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import tensorflow as tf\n",
        "import scipy"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "ELGiUpDAAo-U",
      "metadata": {
        "id": "ELGiUpDAAo-U"
      },
      "outputs": [],
      "source": [
        "#!pip install git+https://github.com/google-research/tensorflow_constrained_optimization\n",
        "import tensorflow_constrained_optimization as tfco"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b27e8c13",
      "metadata": {
        "id": "b27e8c13"
      },
      "source": [
        "To fully understand everything you may want to learn about Gaussian process first. Here the focus is on how to implement its optimization, so only a brief description is provided; feel free to check the literature.. \n",
        "\n",
        "\n",
        "Gaussian process defines a prior probability over functions. Given an arbitrary set of points $x_1,...,x_N$, the assumption is that the function values $y(x_1),...,y(x_N)$ jointly have a Gaussian distribution. The goal is to fit the parameters of this join distribution. The main trick is that the mean and covariance are defined as functions of $x$:\n",
        "$$ y(x) \\sim GP(m(x),k(x,x'))$$\n",
        "\n",
        "$k(x,x')$ is a kernel function which estimates the covariance and can be defined in different ways. The one we used in this example is:\n",
        "\n",
        "$$k(x,x') = \\theta_1 exp(-\\frac{||x-x'||^2}{2\\theta_2})$$\n",
        "\n",
        "And the optimization is to find out proper parameters $\\theta_1$ and $\\theta_2$; both of the parameters are positive, so we have two constraints.\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "3891bb45",
      "metadata": {
        "id": "3891bb45"
      },
      "source": [
        "### Dataset generation"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "06d300c5",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 296
        },
        "id": "06d300c5",
        "outputId": "5673055a-8116-4157-bbfc-b46d7f928394"
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEGCAYAAACHGfl5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAzqklEQVR4nO3deXyU1dn/8c+VfScJJEBIMgkg+xIhKqKyCOL6gNpWbWO1Lk1bly5urQX71LaorVtta9tfam2tRmsr4q4oiigqaNhk30lICJAQCAnZM+f3RwKPQICEzMy5Z+Z6v155Se5M5v6OgVxzn/uc64gxBqWUUsEtxHYApZRS9mkxUEoppcVAKaWUFgOllFJoMVBKKQWE2Q5wqnr16mWysrJsx1BKKb+ydOnSSmNMytHH/bYYZGVlUVRUZDuGUkr5FREp7ui4DhMppZTSYqCUUkqLgVJKKfz4nkFHmpubKS0tpaGhwXYUBURFRZGenk54eLjtKEqpk3BMMRCRnwA3AwZYBdxgjOnSb/XS0lLi4+PJyspCRLwRU3WSMYa9e/dSWlpKdna27ThKqZNwxDCRiPQDfgjkGmNGAKHANV19noaGBnr27KmFwAFEhJ49e+pVmvKawsJCsrKyCAkJISsri8LCQtuR/JpjrgxoyxItIs1ADLDzVJ5EC4Fz6M9CeUthYSH5+fnU1dUBUFxcTH5+PgB5eXk2o/ktR1wZGGPKgEeAEqAcqDbGvHv040QkX0SKRKSooqLC1zGVUg4xc+bMw4XgkLq6OmbOnGkpkf9zRDEQkSRgBpANpAGxInLt0Y8zxhQYY3KNMbkpKccsoFNKBbjlJftYvHUvJSUlHX79eMfVyTllmGgqsM0YUwEgIi8D44HnrKZSSllnjGHR5kr+vGALn23dy5lZyWRmZlJcfOxC2ogeKWzeU8vA1DgLSf2bI64MaBseGiciMdI20DwFWGc5U7fceeedjB49mu9+97tMnDiR1tbW4z62qamJCRMm0NLS4pNs48ePP+lj6uvrHZdbBZ+VO/Yz/U+f8O2/f87WylpmXTqUf9xwBrNnzyYmJuaIx0ZERTP+6tvISI62lNa/OaIYGGOWAC8By2ibVhoCFFgN1Q1btmzhk08+YeXKleTk5HDllVcSGhp63MdHREQwZcoUXnzxRZ/k+/TTT0/6mKefftpxuVXwWb2zmqqDTTx05Ug+umcyN5/Xn9jIMPLy8igoKMDlciEiuFwunn7qb3zw1/uIDAtlR1UdD89bj27r23nir/+zcnNzzdGN6tatW8fQoUMPf371//vsmO+7bFRfvn12FvVNrXznH58f8/Wvj03nG7kZVB1s4gfPLT3iay9+7+yT5tqwYQNTp06lpaWFPn36ADB37lwOdVidPHkyP//5z7nggguYNWsW1dXV/PGPf2TlypXce++9vPXWWyc9R0c6et6HHnqIq666itLSUlpbW7nvvvu4+uqriYuLo7a2lu3bt3PxxRdz7rnn8umnn9KvXz9effVVoqOjGT9+PM8//3y3cx/9M1GqqxpbWokMO/6bko78fdE2fv3GWn54/kDumDbYS8n8k4gsNcbkHn3cKfcMAsbgwYO5/vrrycrK4rrrriMzM5Ovttq+//77+cUvfsGePXtYvnw5r732GgAjRozgiy++OOb5zjvvPGpqao45/sgjjzB16tQTPu8rr7xCWloab775JgDV1dXHPM+mTZt44YUX+Nvf/sZVV13FnDlzuOqqq9i6dWu3cit1qtxuw90vfcllo/oyeUhqlwsBwI3nZLFxVw1/+GAzGckxfCM3wwtJA4wxxi8/xo4da462du3aY47ZMH36dPP555+bsrIyM3jw4GO+PmHCBDNmzBhz4MCBI46npaUdc6wrjn7eDRs2GJfLZe655x7z0UcfHX5cbGysMcaYbdu2mYEDBx4+/tBDD5lf//rXHs3tlJ+J8h+PzltvXD99w/ztoy3dep6mllaT97fFZsC9b5pPNlV4KJ3/A4pMB79THXHPINCsWbOGESNGEB0dfcwK3FWrVlFeXk5ERATx8fFHfK2xsZGoqKgjjp133nnk5OQc8zF//vyTPu+gQYNYtmwZI0eOZNasWfzqV786JmtkZOThP4eGhtLS0uKR3Eqdije+3MkfPtjMN8amc9O53WtjEh4awp+vHUN2r1j++MFmvX9wEloMPKympobw8HCio6NJSkqitbX18C/W8vJy8vLyePXVV4mLi+Odd945/H179+6lV69exzR1+/jjj1mxYsUxH18dIjre8+7cuZOYmBiuvfZa7r77bpYtW9ap1+CJ3Ep11eqyau7670rGupL4zRUjPLKCPSEqnGdvOounrs/VFfEnocXAw1avXs2IESMOfz5t2jQWLVpEXV0dV155JY8++ihDhw7lvvvu4/777z/8uAULFnDppZd2+Xwnet5Vq1Zx5plnkpOTw/3338+sWbM6/bzezq3U0V7/cifJMRH89dqxp3Sf4Hj69IgiNjKMfQebeGV5mceeN+B0NHbkDx9OvmfwVUuXLjXXXnvtSR93xRVXmA0bNvggUed4KrcTfybKmdxut9lzoMFrz//4exuM66dvmDVl1V47hz9A7xnYMWbMGCZPnnzSxVuXX345gwYN8mGyE/PX3Mr/rCqtZtPuGkSElPjIk3/DKbphfDbxUWH8fv5Gr53Dn2kx8IEbb7zxpIu3rrvuOh8m6hx/za38h9tt+NnLX/LdfxXhdnv3Bm+PmHBuOjebd9fuZnXZsdOsg50WA6WUNW+sKmfNzgP8aOpphIR4/wbvjedmk6BXBx3SYqCUsqKpxc2j725gSJ94Zozu55NzJkSFkz+hPyEiNLW4fXJOf6ErkJVSVrz4RQnFe+v4x3fO8MlVwSG3Th6o00w7oFcGSikrahpbmDAohUmDfbs3yaFCsLWils17an16bifTYqCUsuKWSQN55oYzrLxLb251c9X/W8yv31jr83M7lRYDpZRPVR1s4qONFRhjrA3XhIeGcNO52SzcWMHS4n1WMjiNFgOllE89uWAz3/nH55Tuq7ea47qzXSTHRujMonZaDLygM7uEdZbuJqYCSem+Op79rJivj00nIznm5N/gRbGRYXx/Yn8+3lTJF9urrGZxgqAuBoWFhWRlZRESEkJWVhaFhYUeed7O7BLWWbqbmAokj7+3CQR+PNUZq9avHeeiX2I063cdu2dIsAnaYlBYWEh+fj7FxcUYYyguLiY/P98jBaGwsJAZM2YAbTuEvffeewDMmjWL22+/vcPvWb169RF7Ey9btowpU6YAcPnll3usUCllS9n+euYuL+W6cS7SEp2xT3FMRBgf3DWRb49z2Y5iX0cNi/zho7uN6lwulwGO+XC5XJ1+jo40Njaa3r17H/584cKFZuLEiea5554zl1xyiWlpaenw+1pbW03v3r0Pf33ixIlm6dKlxhhjWlpaTK9evbqVyxZtVKcO+WRzhRn/4PtmR9VB21E6tL+uyXYEn+A4jeqCdtFZSUlJl453VmVlJYmJiYc/nzBhAsYYHnvsMT788MPjDh2FhIQwfPhw1qxZw6ZNm3C5XIwZMwZo23QmIiKCmpqaYzaWUcpfjB/Qi4/vmezTBWaddc9LK1mxYz/zfjwhaBekOWaYSEQSReQlEVkvIutE5OS7z3dDZmZml4531tG7hJ1oh7CjjRs3jk8++YRf/vKXPPDAA0d8TXcTU/6svLqe5la3IwsBwOiMRDburmXNzgO2o1jjmGIAPAG8Y4wZAowG1nnzZLNnzyYm5sjZDDExMcyePbtbz/vVXcJOtEPYlClTKCs7cqONcePGMWvWLK644gr69fu/Xi26m5jydz98YTl5Ty2xHeO4LhuZRkRoCHODePMbRxQDEekBTAD+DmCMaTLG7PfmOfPy8igoKMDlciEiuFwuCgoKyMvL6/ZzT5s2jXffffe4O4S53W42b95McnLyEd83ZMgQIiMj+elPf3rEcd1NTPmzdeUH+GL7Pi4Y2tt2lOPqERPO+UNSeXXFTlpag7OBnSOKAZANVAD/EJHlIvKUiMR6+6R5eXls374dt9vN9u3bPVIIAG699Vb++9//8tlnn3HBBRcAbfcOPvvsMwDWrl3L1772NaKjj5xR8cQTT/Dggw8SG3vkS3/++ef53ve+55FsSvnac4uLiQwL4Ru56bajnNAVY/pRWdvIx5srbUexwinFIAwYA/zFGHM6cBD42dEPEpF8ESkSkaKKigpfZ+y0k+0SNmLECB577LHDn2/ZsoUhQ4ZQX1/P9ddff8RjdTcx5c9qGpqZu7yM6aPTSIyJsB3nhCYPTuWxq0ZzRlbyyR8cgJwym6gUKDXGHBpUfIkOioExpgAoAMjNzfXutkjddOONN3b6sQMGDGD9+vUdfk13E1P+7O3Vu6hrauXbZzt/Hn9EWAhXjnH21Ys3OaIYGGN2icgOERlsjNkATAG0naBSfu7rY9Lp3yuWUemJtqN0SnOrm2c+3U5Wz1imDnPuPQ5vcEQxaHc7UCgiEcBW4AbLeZRS3RQSIuT60bBLWIhQuKSE3gmRQVcMnHLPAGPMCmNMrjFmlDHmcmPMKfWVbVtgp5xAfxbB7d6XV/GH9zfZjtElIsLlOf1YvLWKsv12u6r6mmOKgSdERUWxd+9e/SXkAMYY9u7dqwvlgtTuAw38p2gHtY3+1233itPb1vi8EmRrDpw0TNRt6enplJaW4uSZRsEkKiqK9PTgvSEXzF74vAS3MeSd1b0V/TZk9ozhjKwk5i4v45ZJA4KmPUVAFYPw8HCys7Ntx1AqqDW3unl+SQkTB6Xg6un15UJe8fWx6by6Yic1jS0kRAXHyv+AKgZKKfsWbqhgT00js89y/nTS47n6jEyuPsP/rmq6I6DuGSil7OuXFM31Z7uYNDjFdpRu23ewiVZ3cNyD1GKglPKooX0TuH/GCMJD/fvXy8ebKhj7m/dYXnJKExv9jn//tJRSjrJyx36+LN0fEDP6RqUnIiIs2LDHdhSf0GKglPKYR9/byC2Fy2zH8Ige0eGMdSXxwfrgmJ2oxUAp5REVNY0s2lTBjJy0gJmOef6QVNaVH2BXdcPJH+zntBgopTzi9ZU7cRu4PKffyR/sJ84fkgoQFENFWgyUUh7x6ooyhvVN4LTegbNP92mpcTxwxUgmDPL/mVEno8VAKdVtVQeb2FJxkMtPT7MdxaNEhG+dlUm/xOiTP9jP6aIzpVS3JcdGUDRrakDOyW9obuWtVeUM6ZPAsLQE23G8Rq8MlFLdYozBGENUeCixkYH3/tKYtg6sLy0ttR3Fq7QYKKW65cvSaqY9/hFrdx6wHcUroiNCOXtAz4C/iazFQCnVLa+sKKN4bx39kgJ3XP38IalsqzzItsqDtqN4jRYDpdQpa2l18/rKcs4fkkqP6MDt7jl5cPsU0/WBe3WgxUApdco+3bKXytrGgJtFdLSM5BgGpsaxYVeN7SheE3h3e5RSPvPKijLio8KY1P7OOZDNvWU88QG8t4EWA6XUKZs2rDej+vUgKjzUdhSvC+RCAFoMlFLdcNGIvrYj+NTP5nxJfFQYMy8dZjuKxznqnoGIhIrIchF5w3YWpdSJzV+7m9J9dbZj+NSBhmZeX1keEC26j+aoYgD8CFhnO4RS6sTqm1q5/YXl/PnDLbaj+NSkwansOtDA2vLAW1PhmGIgIunApcBTtrMopU5s4cYK6ptbuSTIhokObeX54YbA2+PAMcUA+D1wD+A+3gNEJF9EikSkqKIi8H4YSvmLd1aXkxgTzln9k21H8anU+ChG9uvBBwG43sARxUBELgP2GGOWnuhxxpgCY0yuMSY3JSXwW8oq5USNLa28v24P04b19vt9jk/FNWdmcFZ2csDdN3DKbKJzgOkicgkQBSSIyHPGmGst51JKHWV12QEONrVwcZANER2Sd5bLdgSvEKdVNxGZBNxljLnsRI/Lzc01RUVFPsmklDpSZW0jCVHhRIQF35UBtF0d7a5uJLNnjO0oXSYiS40xuUcfD86fpFKqW3rFRQZtIQC4tXAZNz3zhe0YHuW4n6Yx5sOTXRUopez4dEsleU8tZkdVcK0vONpYVzKb9tRSWdtoO4rHOK4YKKWc680vy1lWvJ9ecZG2o1g1rn0W1eKtey0n8RwtBkqpTml1G+at2c3kISlERwR+L6ITGdmvB7ERoVoMlFLBp2h7FZW1jUHXj6gjYaEhnJGdzOKtVbajeIxTppYqpRzu7dW7iAgL4fwhgd+uujNuP38grcddIut/tBgopTplWN8Ebjwnm7gA3PT+VIx1Bdbqa/2pKqU65aozMmxHcJxFmyqpbWwOiKEzvWeglDqpNTurqWloth3DcZ5atJVH3t1oO4ZHaDFQSp2QMYYfPLeMH76w3HYUxxnXvyeb99RSUeP/6w20GCilTmht+QFKquqYNryP7SiOM65/TyAw1htoMVBKndA7q3cRIm37HasjjUhLIC4yTIuBUirwzVuzizOykukZ5KuOOxIWGsIZWUkBsfOZziZSSh1Xyd46Nu6u5ReXBd4G8J7y+NU5JESF247RbVoMlFLHlZEczfw7JpAcq1cFx5MYE2E7gkdoMVBKHZeIMDA13nYMx3vgrXXERYbxwymn2Y5yyvSegVKqQxU1jfzkxRVs2l1jO4rjbdpdwysrymzH6BYtBkqpDs1ft5u5y8toddhuiE509oCebK04yJ4DDbajnDItBkqpDr27ZheZyTEM7q3DRCdzeL3BNv/tYqrFQCl1jJqGZj7ZvJdpw3ojIrbjON6wvgnE+/l6Ay0GSqljLNxYQVOrW1cdd1JYaAiXje5Lsh/PLNLZREqpY7S0Gkan92CsK8l2FL/x4JWjbEfoFi0GSqljXH56Py4/vZ/tGH7HGENji5uocP/bFtQRw0QikiEiC0RkrYisEZEf2c6kVLCqrm+mJZC28PIRYwznP7qQ2W+usx3llDiiGAAtwJ3GmGHAOOBWEdH170pZ8Lt31jPx4Q9xu3VKaVeICP0So1lavM92lFPiiGJgjCk3xixr/3MNsA7Qa1QPKCwsJCsri5CQELKysigsLLQdSTmY2214b+1uRqX3ICREZxF11RhXEut3HaC2scV2lC5zRDH4KhHJAk4HlnTwtXwRKRKRooqKCp9n8zeFhYXk5+dTXFyMMYbi4mLy8/O1IKjjWlG6nz01jUwbru2qT0WuKwm3gZU79tuO0mWOKgYiEgfMAX5sjDmmJ6wxpsAYk2uMyU1JSfF9QD8zc+ZM6urqjjhWV1fHzJkzLSVSTvfumt2EhQjnD9ZicCpyMhMRwS+Hihwzm0hEwmkrBIXGmJdt5/FHDc2th2cx3P7CcoqLSzp8XElJCd/9VxFxkWEMSInltN7xTByU4pczIJRnvbt2F+P696RHjP+3ZLYhISqcuy8cTK4r2XaULnNEMZC2JY5/B9YZYx6zncefVNQ08vbqct5YWc7++ibm/XgCIkJSTDjxvfpQU1l+zPdkZGTQ2OJmTdle5i5va67VKy6S+y4byowcvVUTrIwx/GbGCEL1XkG33DJpoO0Ip6TLxUBEYoEGY0yrB3OcA3wbWCUiK9qP/dwY85YHzxFQPli/m6c+3sbirXtxGzgtNY7LRqXR4jaEhwq/mjGCwbUPk5+ff8RQUUxMDA888AB5eWcCcLCxhaLifTy9aBuxEW1/HaoONlHb0EJmzxgrr03ZISKMH9jLdgy/19zqZs3OA6QnRdPLj3aHO2kxEJEQ4BogDzgDaAQiRaQSeBP4f8aYzd0JYYxZBOjbkZNodRta3G4iw0JpaHZTXt3ArZMHctmoNAb3ObaZWF5eHtB276CkpITMzExmz559+DhAbGQYEwelMHHQ/92D+fuirfzlwy1cNKIPd184hOxesd5/ccq6pxdt46z+yQxP62E7il/bub+ey5/8hAeuGMm3zsq0HafTxJykPa2ILATmA68Cq40x7vbjycBk4FvAXGPMc17OeoTc3FxTVFTky1Na43Yb3l69i8fnb2T66DR+OOU0Dv3cvNFEbPeBBp75dDvPLi7GGPjt10Zx6ai+Hj+Pco5d1Q2Me/B97r5wMLdO9s9hDqcwxnDG7PlMHJTKo1eNth3nGCKy1BiTe/TxzgwTTTXGNB990BhTRdsN3zntN3+VhxljeH/dHh59byPryg8wMDWOoX0TAO8UgUN6J0Rxz0VDuHaci1ufX8atzy9jT80wbjgn22vnVHa9u3YXABfqlNJuExHGZCaxrMS/ZhSddGrpoUIgIk/IcX4DdVQsVPf9bt4Gbv5XEfVNLfz+6hzm/XgCFwzz3T/WtMRoXsw/m1snD+BC7V4Z0Oat2UX/lFjd4tJDxrqS2FZ5kMraRttROq0r6wxqgNfabyAjIheKyCfeiRW8mlvdh1cvTh+dxqxLh/LeHRO5/PR+VmZ5RISFcPeFQ0hLjMbtNvzkxRW8v263z3Mo79lf18TirVVa8D3oULfXZX603qDTs4mMMbNE5FvAhyLSBNQCP/NasiC0vGQfP5+7mqF94nns6hyG9k04PCzkBNX1zWzcXcNNz5Rx5wWDuN2PN/9W/2fj7lpiwkO1GHjQyPQePH/zWeRkJtqO0mmdLgYiMgX4LnAQ6AvcaIzZ4K1gwaSmoZmH523g2cXF9I6PcuyGIkmxEcz5wXjufXkVj763kYiwEL43cYDtWKqbzsxOZul9FxAeqhP6PCUyLNTvpul2ZZ3BTOA+Y8wiERkJvCgidxhjPvBStqCwqrSa7z+3lJ3V9Vx/dhZ3ThtEfJRz78dHhYfyyDdG09zq5sG31xMfFe5X0+fUkYwxiAgRYY7qTBMQNu+p4bUVO7n1/IFEhjl/dX+n/wYYY85vXw+AMWYVcDHwG28FCxapCZGkJkTy0vfH88vpwx1dCA4JDREevzqHK0/vx7A05wxjqa6bt2Y3Ux79kB1VdSd/sOqSzXtq+cMHm1mz85g2a4500mJwghlE5cCUEz1GdWzx1r3c+Z+VuN2G3glRvPyD8X63vWB4aAiPXZ1DTkYiACV79ZeJP3p3zS4qa5vo0yPKdpSAMyaz7d/00u3+cRO5M1cGH4jI7SJyxFiAiEQAZ4vIM8D1XkkXYOqbWvnla2u4pmAxRcVV7Klpm3bm77X0P0U7mPrYQj7epG3F/Ulzq5v563YzZWgq4aE6TORpqQlRZCT7z2Y3nfkbsAloBeaKyM72rSm3th//JvB7Y8w/vZgxIHy2ZS8XPfER//x0O98Zn8XbPzovYN6NTRvWm/4psXz3X0V8sb3KdhzVSUu2VnGgoUVnEXlRriuZpSX7OFmnByfoTDE4wxjzZ9p6B2XSNjQ0xhjjMsZ81xiz3KsJA0Bzq5t75qwE4N/54/jl9OHERDiiYaxHJMZE8OxNZ5HWI5qbnymidJ8OGfmDeWt2ERUewoTTdG8QbxnjSqK+qZUKP1h81pneRI/Q1lU0E7gPWElbjyKrr84fehN9srmSsa4kosJD2bynhn6JMURHOH9WwanaXnmQy/64iIGpccz5wXhthexw89fuZmtlLfkTdHqwtzQ0txIeGuKofwun3JvIGHOXiAwAFgDZwHRgePvCs9XGmKs9ntbP7a9r4tdvrGPOslJ+fskQ8icMCIpl/lm9YnnkG6NpbGl11F9+1bGpw3oD2ovIm/xpw6hOjVUYY7aIyFRjzMZDx9q3qBzhtWR+qNVt+PcXJTz67kaq65u5/fyBXD8+y3Ysn7poxP+NP9c3tQb0lZA/+3xbFX0SonTPCh945tPtfLG9ij99a4ztKCfUlXUGG4/6vNYYs9jzkfzXrFdWM3PuagamxvH6bedy57TBfrHYxBveWb2LCQ8v0CmnDmSM4adzvmTmK6tsRwkKVQebeGtVOTUNzu7nqfPJuql0Xx17ahoAuO5sF09+awwv5o8L+sVYw9MSaGxu5dbnl9HY4slN8VR3bdpTy7bKg45texJoxrqScBtYuaPadpQT0mJwiiprG3l43nqmPLqQR+a1tWga2jeBS0f19ft1A56QkRzDw98Yzaqyah54c53tOOor5q1u27tgmg/boQeznMxERHD8eoPAmd/oI1sqannq423MWVZKU4ub/xmdxo+nDrIdy5EuHN6Hm8/N5qlF2zgzu6fuluYQ89bu4vTMRHonBMY6F6dLiApncO94ljp8sxstBp3w1S0m//nJduYsK+VrY9K56dxsBqbGWU7nbD+9eAjLSvaxaU8Nbc1ulU2VtY2sL6/h7gsH244SVC4Y1pvK2ibbMU7opOsMfEVELgKeAEKBp4wxD53o8aeyzqCwsPCEm8N/lTGGdeU1LNiwh3lrdjHr0mGcmZ3MnpoGQkToFRfZpXMHs6YWt3bFdJCqg02EitAjxvlNEZXndWcPZK8TkVDgSeACoBT4QkReM8as9dQ5CgsLyc/Pp66ubXZLcXEx+fn5AEcUhAMNzcx+Yx0fbtzD7gNt6+pG9Eugrqlt97HUeL207qpDheDzbVVU1zf7dOtOdazk2AjbEYKWk98YOaIYAGcCm40xWwFE5N/ADMBjxWDmzJmHC8EhdXV1fO9Hd/GffS5yMhK5f8YIYiPCWLJtL7muZCYOTmHSoBRSdWy124wx/O6d9WypqOW9OybqlZUFO6rquOu/K7nvsmGM6NfDdpygYozhgsc/4uz+Pfn15c5cnuWUEtUP2PGVz0vbjx1BRPJFpEhEiioqutYhs6SkpMPjB/fuJiE6nN7tTeNCQ4QFd03iybwxXJWboYXAQ0SEB68cycHGVu5/3WM1XnXBm6vKWbKtih7ROjzkayJCanwkK3bstx3luJxSDDrFGFNgjMk1xuSmpHStuVZmZse7cblcmTx701ncMmng4WM6NdQ7Tusdz23nD+T1lTuZv3a37ThB580vyxmdkUhGsq46tiEnI5F15QdoaHbmuhunFIMyIOMrn6e3H/OY2bNnExNz5D+CmJgYZs+e7cnTqJP4/sQBDOkTz6xXVnPA4SsyA0nJ3jpWlVVz6UhdaGZLTkYiLW7D6jJnLj5zSjH4AjhNRLLbN825BnjNkyfIy8ujoKAAl8uFiOByuSgoKDjubCLlHRFhIfz2a6O4+bxsYgOojbfTvbmqHIBLRur0XltyMhMBHDtU5Ih/jcaYFhG5DZhH29TSp40xazx9nry8PP3l7wCjMxIZ3b5d5qEN2ZV3pSdF880zM0lP0iEiW1Ljo7ht8kBGpSfajtIhx6wz6Cp/2M9Andj8tbv584ebKbx5nHY3VcpHjrfOwCnDRCoIxUSEsqxkP39ZuMV2lIC2YVcN1fV6f8YJWlrdrNlZ7cgOploMlDXjB/Zi+ug0/rpwCzuqtNW1t9zxnxXc9M8vbMdQwMrS/Vz6h0V8tmWv7SjH0GKgrLr3kiGEivCbN3XtgTdsqzzImp0Hjth0SNkzPK0HYSHCcgfeRNZioKzq2yOa284fyLw1u1npwH8g/u6t9llEF+ssIkeICg9lWFoCK0r2245yDEfMJlLB7ebzshmelsCodG2R4GlvflnO6ZmJ9EuMth1FtcvJSGTO0lJa3cZRe4XrlYGyLjIslEmDUxERx67O9Ec7qupYW36AS/WqwFFyMhI52NTa3tbdObQYKMd4dUUZ5/52ARU1jbajBISM5Bjev3MiV45Jtx1FfcWEQSn868YzyXRYWxAtBsoxRvbrQXV9E797Z73tKAFjQEqctqx2mF5xkUwYlEKMw1bgazFQjtE/JY4bz83mv0tLWe7wLQKdbmtFLbc9v4xtlQdtR1EdWF1WzbOLi23HOIIWA+Uot59/GqnxkfzytTW43f65Ot4JXl5WxluryonRld2O9N7a3fzi1dXUNrbYjnKYFgPlKHGRYdx7yRBWlVWzTK8OTonbbZi7vIxzT0vRTe8dKiczEWPgy9L9tqMc5qxBK6WAy3P6MSo9kQEpcbaj+KXF2/ZStr+eey7STe+dKqe9Wd3ykv2MH9DLbph2emWgHEdEDheCPTUNltP4nzlLy4iPDOPC4brq2KmSYiPI7hXrqHbWWgyUY/3rs+1M/N2H7NxfbzuKX+mfEst3zskiKlzvFzjZ6RmJrN91wHaMw7SFtXKsHVV1THlsIZeN7MtjV+fYjqOUR+072ERcVBjhob59T64trJXfyUiO4aZzs3l5eZmjbrQ52dLiKlpa3bZjqE5Iio3weSE4EeckUaoDt0waQM/YCH7zxjr89SrWFwoLC0nPyCQ3qxepaRkUFhbajqQ64cG31/HPT7bZjgFoMVAOFx8Vzh3TBrG2/ADb9+qeBx0pLCwkPz+fstIdgGHfnp3k5+drQfADRdv38fqX5bZjAFoMlB+4OjeDBXdNIrtXrO0ojjRz5kzq6o4slHV1dcycOdNSItVZORmJrC6rptkBQ3taDJTjhYWGkBIfiTFG2yt0oKSkpEvHlXPkZCTS2OJmwy77HUy1GCi/8eDb65n+p0VUHWyyHcVRMjMzu3RcOUdORiKAI3Y+s14MRORhEVkvIl+KyFwRSbSdSTnT18emc7CxhSfmb7QdxVFmz55NTMyR7ZBjYmKYPXu2pUSqs9KTohmelkCrDhMB8B4wwhgzCtgI3Gs5j3KoQb3j+eaZmTy3pITNe2ptx3GMvLw8CgoKcLlciAgul4uCggLy8vJsR1MnISK8+cPz+M452baj2C8Gxph3jTGHWvctBnQnDnVcP7lgEDHhoTzw1jrbURzjrv+uZH/aWWzfvh2328327du1EPgh21OnrReDo9wIvH28L4pIvogUiUhRRUWFD2Mpp+gVF8lt5w9kw64avXdAW++mucvLHNUKWXXN+l0HOPe3H7Boc6XVHD4pBiIyX0RWd/Ax4yuPmQm0AMedHG2MKTDG5BpjclNSUnwRXTnQDedk8/6dE3UHL9r2LWh1G76mW1v6rb49oindV8+Kkv1Wc/ikhbUxZuqJvi4i3wEuA6YY29dKyvEiwtrewzQ0t7J5Ty0j+vWwnMiOllY3z35WzFnZyQxM1Xbf/qpHdDj9U+x3MLU+TCQiFwH3ANONMbrEVHXaPS99yfVPf86BhmbbUayYv243ZfvrucEBNx9V9+RkJLKydL/V+wbWiwHwJyAeeE9EVojIX20HUv7hu+f1p6quiSc/2Gw7ihVZvWK5/mwXFwzrbTuK6qbTMxKprG2idJ+9du3Wdzozxgy0nUH5p5HpPfjamHT+8cl28s5ykdkz5uTfFECG9Eng/hkjbMdQHjCuf0++eabdRYJOuDJQ6pTdfeFgwkKFB98Orqmm/yna4aiNUVT3nNY7ngevHElGsr03NFoMlF/rnRDF9ycOYPeBBuqbWm3H8YmKmkZmzV3N80u091AgcbsNO6rs3TbVYqD83g8mDWDOD8YTHREc2zw+v6SEplY314/Psh1FedBv31nP1McWWutgqsVA+b3w0BBEhIqaRj6xvHDH2xpbWnluSTGTBqcwIEWnkwaSEf160NjiZn25nQ6mWgxUwJg5dxU/eG4plbWNtqN4zZtfllNR06jTSQPQoQ6mKyxt8arFQAWMey4aQn1zKw+9vd52FK+pOtjE6PQeTDitl+0oysPSk6LpFRdhbSWyFgMVMAamxnHzef3557+epW96BiEhIWRlZQXU9o83n9efubecg4jYjqI8TETIyUhkxY59Vs5vfZ2BUp6UuvsL9s37E+7mtqGi4uJi8vPzAfy+k+em3TUMTI0jJEQLQaC6+bz+1DXZaTqoVwYqoNz/v/cdLgSHBMJ+wGX767noiY/528dbbUdRXjSuf0/OH2JnRbkWAxVQAnU/4H99uh1jDJeM7Gs7ivKyJVv3srS4yufn1WKgAkog7ge8q7qBZz7bzv+MTiM9KbhabgSjn89dxZ8XbPH5ebUYqIDS0X7AkVHRfr0f8OPvbaTVbbhr2mDbUZQP5GQksWKH7zuYajFQAeXo/YAjE1MZcMUdfP3qa2xHOyX1Ta0s2lzJt8dlWe1bo3wnJzORvQd938FUi4EKOHl5eYf3A3578SoOZp7N7+dvsh3rlERHhDL/jon85ILTbEdRPnL6ocVnPt7sRouBCmiTB6dyzRkZ/HXhFpaV2Jm/farK9tfT2NJKdEQo8VHhtuMoHxncJ57IsBCfFwNdZ6AC3qzLhtEvMZphfRNsR+k0Ywy3PLeUyPBQ/vO9s23HUT4UHhrC67efi8vH+3PolYEKeHGRYdw+5TSiwkNpbPGPNtdvripnZWk13xirG90Ho0G944kM820XXi0GKmhsrzzI1McW8t7a3bajnFBTi5uH521gSJ94rhyjxSAYlVfXM/vNtWza7bsOploMVNBIS4wmPjKcn8350tGdTZ9fUkzx3jp+evEQQrX1RFByG/jbx9t82pJdi4EKGhFhITx+dQ41jS38bM4qn8/j7qwFGyoYP6Ankwal2I6iLEnrEUWfhCiW+bCDqWOKgYjcKSJGRLQ3r/KawX3iuefCwcxft5v/FpXajtOhf3znDJ781hjtTBrERIQxrkSfzoBzRDEQkQxgGuDfDWSUX7jxnGzO7t+TeWt2OerqYF35AfYcaCAkREiKjbAdR1k2JjOJ0n317Klp8Mn5nDK19HHgHuBV20FU4AsJEf567Vjio8Ic8+5738Embn6miJT4SObeMt4xuZQ9p2cmkRwbwY6qelLjo7x+PutXBiIyAygzxqy0nUUFjx4x4YSEtO2bPOuVVTQ025ty6nYb7vjPCipqGrl/+nAtBApoW4m8dNZUxrqSfHI+n1wZiMh8oE8HX5oJ/Jy2IaLOPE8+kA/+3YVSOceKHfspXFJCZU0TT+aNsTJ75y8Lt7BgQwW/njGc0e2tCJTy9SZGPrkyMMZMNcaMOPoD2ApkAytFZDuQDiwTkY4KB8aYAmNMrjEmNyVFZ1qo7rtgWG9mXjKUd9bs4jdvrvX5+T/fVsWj727gf0ance04l8/Pr5zt1RVlXPT7j2hqcXv9XFbvGRhjVgGphz5vLwi5xhjfTa5VQe+mc7Mp3VfPPz7ZTr/EaG4+r7/Pzj24TzzXnZ3FXRcO1uEhdYzQEGH9rhrW7zrAqPREr57L+j0DpWwTEe67bBgXDe/DC5+X+OT+QUurm8aWVnpEh/PL6cOJi3TKXA7lJGMy2+4XLCv2/hRTRxUDY0yWXhUoG0JDhN9fk8NL3x9PVHio16ecPvbeRr7+l8+sbX6u/ENaYrTPFp85qhgoZVNUeChJsREYY/jes0t5csFm3G7PFoWmFjf3v76GP3+4heFpCcRE6BWBOjFfLT7Tv4lKHaWxxU1EWAgPz9vAkm1VPH7VaHrGRXb7eXfur+fW55exvGQ/N5yTxb0XD/VAWhXopg7tTWxEGM2tbsJDvff+XZy0ArMrcnNzTVFRke0YKkAZY3j+8xLuf30tSTHh/OGa0zmrf89uPee3/76E5SX7+e3XRnHpqL4eSqpU14jIUmNM7jHHtRgodXxrdlZz2/PLqW9q5cO7JxEV3rUe861uQ1OLm+iIUEr21tHidtM/Jc5LaVWgMsZQXd9MYkz325QcrxjoMJFSJzA8rQev334uJXvriAoPpbK2kX99VszkwSmMSk887iK1ytpGlhXv49nFxcRGhPGXa8eQ6eOdq1TguO7pz2lobuW/3x/vtXNoMVDqJOIiwxiW1rZl5sod+/nTB5v4w/ubSIoJZ+KgFCYPSW0b140M428fbeW59v0IAKLCQ/jf/xluM74KAIN6x/Pc4mKa2u9neYMWA6W6YMrQ3iy77wI+2lTJh+v38OHGCl5ZsZNFP51MbGQYYaHCkD7xfOvMTMa6khjRr0eXh5aUOtqYzCT+vmgb68oPeK1liRYDpbooMSaC6aPTmD46DbfbsH5XDelJbUNAN5yTzQ3nZFtOqALN6ZmJACwr2ee1YqDrDJTqhpAQOTyEpJS3+GLxmV4ZKKWUH7jnosGkxHd/vcvxaDFQSik/cOWYdK8+vw4TKaWUH2hpdbO0uIptlQe98vxaDJRSyg+0GsM3C5bw78+9s1W8FgOllPIDkWGhJJQt5jffnkxISAhZWVkUFhZ67Pn1noFSSvmBwsJCvvz372hubACguLiY/Px8APLy8rr9/HploJRSfmDmzJmHC8EhdXV1zJw50yPPr8VAKaX8QElJx/cKjne8q7QYKKWUH8jMzOzS8a7SYqCUUn5g9uzZxMQc2fk2JiaG2bNne+T5tRgopZQfyMvLo6CgAJfLhYjgcrkoKCjwyM1j0M1tlFIqqBxvcxtHXBmIyO0isl5E1ojI72znUUqpYGN9nYGITAZmAKONMY0ikmo7k1JKBRsnXBn8AHjIGNMIYIzZYzmPUkoFHScUg0HAeSKyREQWisgZx3ugiOSLSJGIFFVUVPgwolJKBTafDBOJyHygTwdfmtmeIRkYB5wB/EdE+psO7mwbYwqAAmi7gey9xEopFVyszyYSkXeA3xpjFrR/vgUYZ4w54Vt/EakAik/xtL2AylP8Xn+lrzk46GsODt15zS5jTMrRB63fQAZeASYDC0RkEBBBJ15kRy+ms0SkqKOpVYFMX3Nw0NccHLzxmp1QDJ4GnhaR1UATcH1HQ0RKKaW8x3oxMMY0AdfazqGUUsHMCbOJbCiwHcACfc3BQV9zcPD4a7Z+A1kppZR9wXploJRS6iu0GCillAq+YiAiF4nIBhHZLCI/s53H20QkQ0QWiMja9kaAP7KdyRdEJFRElovIG7az+IKIJIrIS+0NH9eJyNm2M3mbiPyk/e/0ahF5QUSibGfyNBF5WkT2tM+2PHQsWUTeE5FN7f9N8sS5gqoYiEgo8CRwMTAM+KaIDLObyutagDuNMcNoW+V9axC8ZoAfAetsh/ChJ4B3jDFDgNEE+GsXkX7AD4FcY8wIIBS4xm4qr/gncNFRx34GvG+MOQ14v/3zbguqYgCcCWw2xmxtn9L6b9o6pgYsY0y5MWZZ+59raPsl0c9uKu8SkXTgUuAp21l8QUR6ABOAv0PbdG1jzH6roXwjDIgWkTAgBthpOY/HGWM+AqqOOjwDeKb9z88Al3viXMFWDPoBO77yeSkB/ovxq0QkCzgdWGI5irf9HrgHcFvO4SvZQAXwj/ahsadEJNZ2KG8yxpQBjwAlQDlQbYx5124qn+ltjClv//MuoLcnnjTYikHQEpE4YA7wY2PMAdt5vEVELgP2GGOW2s7iQ2HAGOAvxpjTgYN4aOjAqdrHyWfQVgjTgFgRCbrFq+3dGjyyPiDYikEZkPGVz9PbjwU0EQmnrRAUGmNetp3Hy84BpovIdtqGAc8XkefsRvK6UqDUGHPoiu8l2opDIJsKbDPGVBhjmoGXgfGWM/nKbhHpC9D+X4/sARNsxeAL4DQRyRaRCNpuOL1mOZNXiYjQNpa8zhjzmO083maMudcYk26MyaLt5/uBMSag3zEaY3YBO0RkcPuhKcBai5F8oQQYJyIx7X/HpxDgN82/4jXg+vY/Xw+86okntd6byJeMMS0ichswj7bZB08bY9ZYjuVt5wDfBlaJyIr2Yz83xrxlL5LygtuBwvY3OVuBGyzn8SpjzBIReQlYRtuMueUEYFsKEXkBmAT0EpFS4H+Bh2jb9+Um2tr4X+WRc2k7CqWUUsE2TKSUUqoDWgyUUkppMVBKKaXFQCmlFFoMlFJKocVAKaUUWgyUUkqhxUApj2nfN+KC9j//RkT+aDuTUp0VVCuQlfKy/wV+JSKptHWHnW45j1KdpiuQlfIgEVkIxAGT2vePUMov6DCRUh4iIiOBvkCTFgLlb7QYKOUB7a2EC2nrsV8rIkdvVaiUo2kxUKqbRCSGtn76dxpj1gG/pu3+gVJ+Q+8ZKKWU0isDpZRSWgyUUkqhxUAppRRaDJRSSqHFQCmlFFoMlFJKocVAKaUU8P8BSjxE/hw1uHUAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "# Lets assume the generative process is the follwoing:\n",
        "y_process = lambda x: x * np.sin(x) ## the same function as Sikit learn example on GP\n",
        "x_full_range = np.linspace(start=0, stop=10, num=50).reshape(-1, 1)\n",
        "\n",
        "# and the samples drwan from this process are:\n",
        "X = np.linspace(start=0, stop=10, num=6).reshape(-1, 1)\n",
        "y = np.squeeze(y_process(X)).reshape(-1, 1)\n",
        "\n",
        "plt.plot(x_full_range, y_process(x_full_range), linestyle='--', label=r\"$f(x) =x\\sin(x)$\",)\n",
        "plt.plot(X, y, 'ko', label=r\"$(x,y)$\")\n",
        "plt.legend()\n",
        "plt.xlabel(\"$x$\")\n",
        "_ = plt.ylabel(\"$f(x)$\")\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "aa6d601d",
      "metadata": {
        "id": "aa6d601d"
      },
      "source": [
        "### Loss function"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "8649098d",
      "metadata": {
        "id": "8649098d"
      },
      "source": [
        "The loss function is negative log likelihood, so we would like to find the parameters which fit a multivariate Gaussian distribution best:\n",
        "$$ -ln P(y|X,\\theta_1, \\theta_2) = \\frac{1}{2}(y-\\mu_y)^T \\Sigma_{\\theta}^{-1}(y-\\mu_y) + \\frac{1}{2}ln|\\Sigma_\\theta|$$\n",
        "\n",
        "The loss function only shows the terms which are dependent on the trainable parameters. The first term is called 'data fit', and the second term is called 'complexity penalty'. The implementation is:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "76542615",
      "metadata": {
        "id": "76542615"
      },
      "outputs": [],
      "source": [
        "def objective():\n",
        "    mean_y = lambda : np.mean(y_tr)\n",
        "    y_tr_norm = y_tr - mean_y()\n",
        "    Covariance = Kernel(X_tr,X_tr)\n",
        "    Kernel_inv = tf.linalg.inv(Covariance)\n",
        "    data_fit  = 0.5*tf.tensordot(tf.matmul(y_tr_norm, Kernel_inv, transpose_a=True), y_tr_norm, axes=1)\n",
        "    complexity_penalty = 0.5*tf.math.log(tf.linalg.det(Covariance))\n",
        "    return complexity_penalty + data_fit[0]"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b0a69621",
      "metadata": {
        "id": "b0a69621"
      },
      "source": [
        "Let's make the Kernel function:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "e2963382",
      "metadata": {
        "id": "e2963382"
      },
      "outputs": [],
      "source": [
        "def Kernel(x1,x2):\n",
        "    kernel = theta1 * tf.math.exp(-tf.math.squared_difference(x1,tf.transpose(x2))/2*theta2)\n",
        "    return kernel"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ba2d5d74",
      "metadata": {
        "id": "ba2d5d74"
      },
      "source": [
        "### Model fitting"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "399742c5",
      "metadata": {
        "id": "399742c5"
      },
      "outputs": [],
      "source": [
        "# setting the data in TF constant form\n",
        "X_tr = tf.constant(X,dtype=tf.float32)\n",
        "y_tr = tf.constant(y, dtype=tf.float32)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "3d605a13",
      "metadata": {
        "id": "3d605a13"
      },
      "source": [
        "we can make the class now"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "dc7d8431",
      "metadata": {
        "id": "dc7d8431"
      },
      "outputs": [],
      "source": [
        "class FitGP(tfco.ConstrainedMinimizationProblem):\n",
        "    def __init__(self, theta1, theta2):\n",
        "        self.theta1 = theta1\n",
        "        self.theta2 = theta2\n",
        "        \n",
        "    def Kernel(self,x1,x2):\n",
        "        kernel_1 = self.theta1 * tf.math.exp(-tf.math.squared_difference(x1,tf.transpose(x2))/2*self.theta2)\n",
        "        return kernel_1 \n",
        "        \n",
        "    def objective(self):\n",
        "        mean_y = lambda : np.mean(y_tr)\n",
        "        y_tr_norm = y_tr - mean_y()\n",
        "        Covariance = self.Kernel(X_tr,X_tr)\n",
        "        Kernel_inv = tf.linalg.inv(Covariance)\n",
        "        data_fit  = 0.5*tf.tensordot(tf.matmul(y_tr_norm, Kernel_inv, transpose_a=True), y_tr_norm, axes=1)\n",
        "        complexity_penalty = 0.5*tf.math.log(tf.linalg.det(Covariance))\n",
        "        return complexity_penalty + data_fit[0]\n",
        "    @property\n",
        "    def num_constraints(self):\n",
        "        return 2\n",
        "    \n",
        "    def constraints(self):\n",
        "        constraint1 = -self.theta1\n",
        "        constraint2 = -self.theta2\n",
        "        return tf.stack([constraint1, constraint2])\n",
        "    "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "980d0e8f",
      "metadata": {
        "id": "980d0e8f"
      },
      "outputs": [],
      "source": [
        "# Initializing the parameters\n",
        "theta1= tf.Variable([1], dtype=tf.float32)\n",
        "theta2=tf.Variable([.1],dtype=tf.float32)\n",
        "## Instantiating the problem\n",
        "problem = FitGP(theta1, theta2)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "b1dadab0",
      "metadata": {
        "id": "b1dadab0"
      },
      "outputs": [],
      "source": [
        "optim = tf.optimizers.Adagrad(1.)\n",
        "optimizer = tfco.LagrangianOptimizer(optim,\n",
        "                                    num_constraints=problem.num_constraints)\n",
        "var_list = list(problem.trainable_variables) + optimizer.trainable_variables()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "60a2dd60",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "60a2dd60",
        "outputId": "ce38d61a-22f8-4b1b-d702-16770c7d88b3",
        "scrolled": true
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "iteration 0, the objective is [32.585915] theta1 =[2.], and theta2=[1.1]\n",
            "iteration 100, the objective is [27.130434] theta1 =[2.4936635], and theta2=[1.110772]\n",
            "iteration 200, the objective is [24.518663] theta1 =[2.8359442], and theta2=[1.1196165]\n",
            "iteration 300, the objective is [22.89058] theta1 =[3.1058428], and theta2=[1.1273657]\n",
            "iteration 400, the objective is [21.744608] theta1 =[3.3318968], and theta2=[1.1343697]\n",
            "iteration 500, the objective is [20.878738] theta1 =[3.528057], and theta2=[1.1408195]\n",
            "iteration 600, the objective is [20.19324] theta1 =[3.7023127], and theta2=[1.1468331]\n",
            "iteration 700, the objective is [19.632254] theta1 =[3.859712], and theta2=[1.1524903]\n",
            "iteration 800, the objective is [19.161654] theta1 =[4.00367], and theta2=[1.1578485]\n",
            "iteration 900, the objective is [18.759207] theta1 =[4.1366186], and theta2=[1.1629506]\n"
          ]
        }
      ],
      "source": [
        "for ii in range(1000):\n",
        "    optimizer.minimize(problem, var_list)\n",
        "    if ii % 100 == 0:\n",
        "        print(f'iteration {ii}, the objective is {problem.objective()}', end=\"\")\n",
        "        print(f' theta1 ={problem.theta1.numpy()}, and theta2={problem.theta2.numpy()}')"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "84b43c31",
      "metadata": {
        "id": "84b43c31"
      },
      "source": [
        "### Prediction"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "fc141dd7",
      "metadata": {
        "id": "fc141dd7"
      },
      "source": [
        "Gaussian process is a regression model and we can predict the mean and standard deviation for any new given point $x_{test}$. The prediction function is provided below. For more information on how it is done, please refer to literature (e.g. https://peterroelants.github.io/posts/gaussian-process-tutorial/)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "3b56204f",
      "metadata": {
        "id": "3b56204f"
      },
      "outputs": [],
      "source": [
        "def predict(X_tr,y_tr, X_ts, mu_tr,mu_tst, kernel_func):\n",
        "    Sigma11 = kernel_func(X_tr, X_tr)\n",
        "    Sigma22 = kernel_func(X_tst, X_tst)\n",
        "    Sigma12 = kernel_func(X_tr, X_tst)\n",
        "    solved = tf.transpose(tf.linalg.solve(Sigma11, Sigma12))\n",
        "    mu21 = mu_tst + tf.tensordot(solved, (y_tr - mu_tr), axes=1)\n",
        "    Sigma2 = Sigma22 - tf.tensordot(solved, Sigma12, axes=1)\n",
        "    return mu21.numpy().ravel(), np.diag(Sigma2)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "52152e68",
      "metadata": {
        "id": "52152e68"
      },
      "outputs": [],
      "source": [
        "X_tst = np.linspace(start=0, stop=12, num=100).reshape(-1, 1)\n",
        "X_tst = tf.constant(X_tst, dtype=tf.float32)\n",
        "mu1 = mu2 = np.mean(y_tr)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "d7e64a23",
      "metadata": {
        "id": "d7e64a23"
      },
      "outputs": [],
      "source": [
        "mu21, Sigma2 = predict(X_tr,y_tr, X_tst, mu1,mu2, problem.Kernel)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "2657beaa",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 408
        },
        "id": "2657beaa",
        "outputId": "ded39d30-d735-4f93-ad22-113c84452873"
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<matplotlib.legend.Legend at 0x2420daaeee0>"
            ]
          },
          "execution_count": 76,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAF2CAYAAABJSW4pAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAACEmElEQVR4nO2dd3hb1d2A3yPJe8Yje+8dZxBIgCQQEhIKIWGUEVYZAcpuaaFNy0cpKbRQWiirYbeYvQJlhbBHyB5kT2fHifeQba3z/XEkW7ZlWZI17fM+jx9J917de6517/2d3xZSSjQajUajaSuGSA9Ao9FoNO0DLVA0Go1GExS0QNFoNBpNUNACRaPRaDRBQQsUjUaj0QQFLVA0Go1GExRMkR5AJMnJyZF9+/aN9DA0Go0mplizZk2RlDK36fIOLVD69u3L6tWrIz0MjUajiSmEEPs8LdcmL41Go9EEBS1QNBqNRhMUtEDRaDQaTVDQAkWj0Wg0QUELFI1Go9EEBS1QNBqNRhMUtEDRaDQaTVDQAkWj0Wg0QUELFI1Go9EEBS1QNBqNRhMUtEDRaDQaTVDQAkWj0Wg0QaFDF4fUaDQ+UF0NFgvEx6u/uLhIj0gTpWiBotFoPCMlHD0K+/aBEGqZENCpE/TtqwWLphlaoGg0muZYrbBnD5SWQmYmGNys4+XlsGkTDB4MKSkRG6Im+tA+FI1G05w9e6CiArKyGgsTgPR0MBph40Y4fjwy49NEJVpD0Wg0jamsVJpJVlbL2yQmKpPXrl3KNNa5c/jGp4latIai0Wgac+CAEhitYTQqc9iuXVBYGPJhaaIfLVA0Gk0DFRXKR5Kc7Nv2RqNy0u/erYWKRgsUjUbjREoV0eWvo91dqBw+HJqxaWICLVA0Go2ivFzlnPhi7mqK0ah8Lvv2KZOZlOTn59O3b18MBgN9+/YlPz8/+GPWRBXaKa/RaBRHjvhu6vKEwaCEysGD5L/1Fgv+8AfMZjMA+/btY8GCBQDMnz8/GKPVRCFaQ9FoNCrvpKIiMO3EnaNH4b77WHjX3Vji0kgaeCKm7J4AmM1mFi5cGITBaqIVraFoNBqoqmr7PnbupHbB9dx80tU4bvoFPRIafDEVq96j9Itn2b9/f9uPo4latEDRaDRQVAQJCYF/f/VquPFGEisq6B9nw7D7R4r3b8FyfB+JPYdjKVKCpHfPXlSWVZGWntw8YVIT82iBotF0dOx2KClROSUBIDdv4Z///oQZSbmMPPFEfn/nefT58EMWLPsGS00NlsPbAEgGRs+9i8v+/T3/PTmd9IxUSE2FtDRlaktJaagZpgkttbVqAhHk/3dMTBGEEEOEEOvd/iqEELc32WaaEKLcbZt7IjRcjSa2qK5Wr4E8XKTkxf8u49FJF/Hx2VfBo4+C3c78m25i8TPP0KdPH4QQ9ElJYTFw8zdvs6XcxhUrzFRYpcrI371b1QbbtElFmkkZzLPTNKWmBrZuBWfARDCJCYEipdwupcyTUuYB4wEz8K6HTb91bSelvC+sg9RoYpWSkoArB6/58FsWDZjBGQVr+PUNsxrMWBkZzJ8/n4KCAhwOBwUlJcwfM4aZG77gCctGNhXVctP3xciUFJXDkpWlBMnWrbBliyqXrwk+Nhvs2BEcn5kHYkKgNGE6sFtKuS/SA9FoYh4plf8kKcnvrxZX1HDz3kS6VRbx98FgSE9Xs96srOYCKj4enBFeM5/7G/eMzeDbQ2aW7Kls2CYxUX23thZ27lSmOE3wcDiUNmi1ts1f5oVYFCgXA6+2sG6SEGKDEOJjIcSIcA5Ko4lJqqvVg9to9Purz721nOKEFJ5a/V8yLr5ALbRYWi4Ued55MGoUlJUxf9USzhuYRs9UD5pRWpoa19692vwVTA4dgrIyVS06RMSUQBFCxANzgDc9rF4L9JFSjgH+BbzXwj4WCCFWCyFWH9eltzUdnfLywKKtSkv51dO/441X7mbk9fPBZFIzYINBOdo9YTTC7bert88/zyMT0pnQpQXNKDNTaU4HD/o/Nk1zqquVQAkw8MJXYkqgALOBtVLKZlXopJQVUsoq5/uPgDghRI6H7RZLKSdIKSfk5uaGfsQaTTRTXOy3uavW5qD8lTcwlZeR1y8HpkxRK6qrlXbiTds55xwYMUL5bV57jWqrgz/9eJzVhTXNt+3USQmUENn7OwwuU1dSUsij6GJNoFxCC+YuIURXIdR/SwgxEXVuxWEcm0YTW9hsyl8RH+/X157bVMpUw0SOJ2fCtdc2PKSsVsjO9v7l1FS45hr1/vnnwWZj6b4q7v3xOLKpeUsI5VfRBSfbxvHjKrIrAD+Zv8SMQBFCpAAzgHfclt0ghLjB+fECYJMQYgPwGHCxbHaFajSaempr/fZRFNXYeGp9ESfs+4ncTikwaZJa4XL0tlapODERTj0V+vWD48dJWbeaW8dmsam4jh+OeNBSUlKUNlPjYZ2mderqVMHOEPpN3IkZgSKlrJZSZkspy92WPS2lfNr5/nEp5Qgp5Rgp5UlSyh8iN1qNJgYIIA/h0XUl1Njh7q9fhAsuaPC/1NRAbm7rJhUhVCTXGWeozx99xLn908hJMrL4p1LP34mLUzXCNP6zb5/ybwUQdBEIMSNQNBpNkCkv96sYZEGFhVe2l3Pphk8ZUHZERW25sNtVdJYvdOoEp5+u3n/2GYnSzi+GZ/L1ITPbS+uab5+aCseOqdm2xneqq5V211KQRAjQAkWj6YhIqQSKH/kIy/ZXY5QObvkuH6ZOha5dG1a6/B2+kJysTF6DBqkw1uXLmT80gwsGpZFg9KDhCKE0IR2V6R+HDvntH2srWqBoNB2RurqGMF8fuXZEJl998H90ri5T5i4XrjwWX4VTQoJ60M2cqT5//DGZCUYePrUrfdNbeACmpSnnvM3m83g7NGZz2LUT0AJFo+mY+OmQt9glrFlD9y3rlK9k2rSGlXV1kJHh3/E7dYLTTlPvP/usvtTKluI6/ueePe/CYFDjrfSwTtOcI0fCrp2AFigaTcekosLn+l1Wh2TGO/tY/PVutWDuXOXodRGIQMnMhF69YMgQJSS+/x6AJzaUsPCHY9TZHc2/Ex+v8mY03qmpUebB1iLuQoAWKBpNR6SszOcZ7CcFVeyrtNJ/zbdqwaxZzTfyt3Wwa/uzzlKvH38MwIWD0ym3OPjqoIcItKQkVZ3Y4UHYaBo4elRNFiLQCkALFI2mo2G3q1msjwLluc1l9E1wcPqaZdCtm8p0b4q/rYPj45WW4/KjLFsGFgundE8mJ9HIkt0tmL0cjpCUXW83WCwqIi7MvhMXWqBoNB2NmhqfZ6/rjtWy/ngtvyjdjAGp8kfcv2uxKNNKIHkOqalKQA0erEJcN2zAZBCc0z+NZQeqKa/zUG3YYFDRaRrPuP43EWpUpgWKRtPR8CPr/I2d5SSZBOd/+h+1YPr0xhsE4j9xkZqqMuxPOkl9/kHlIs8bmEaCQbC91ENPlORkVTRS45nDhyPiO3GhBYpG09EoK/M5xHfByE78fTCk7tiqBMcJJzTewJ+ExqYkJ6vvu8q3LF8OwKjsBFZd0o+JXT3UnoqLUxFqtbWBHbM9U10dUG22YKIFikbT0aio8Pmh0y8jnrM2fqE+nHZa4+guUKG8gRYddAm1iROVyWzjRqiqQghBosmAlBKro4XQZl2BuDnFxWErsdISWqBoNB0Jq9XnhlqPry/hh8Nm5TCHhvpbLhwOtZ9AZ8Su76WmwujRalyrVgFgtjqY8e5+nvFU3ysxUYcPN8Vuh8LCiJq7QAsUjaZj4WOv9qPVNh5ZV8wPu47Dhg1Kmzj55MYbufwngTqAjUYlHNz9KE6zV3Kcgcx4A+95ivZKTFRmO90iuIHKyoA7bwYTLVA0mo5EXZ1PGfLv7q7AIeH8A2vUglNOaZ5rUlfX9rLoLsd8Ez8KwNn909hRZmFveRMh6BJg1dVtO3Z74ujRsPQ7aQ0tUDSajkRlZasZ8lJK3t5VyfjOifT77jO10FUmpSltfYilpSmtaexYpXns2FEfxTW9lzLffH7Ag+AwGHQ+iguLRYULa4Gi0WjCSlVVqz6PjUV17CqzcEH/FPjxR7Xw1FObbyhl2yOKEhMb9jNhglrmPGavtDiGdIr3LFASEnQ+ioso0tS0QNFoOgpSqodPKxpKca2dgZnxnFW5R2kBAwY0LlXv2he0XaAkJDSYsJr4UQBuHpPFZUM95LnEx6toNd2UVWl0frQhCCWm1jfRaDTtAotFPYBbcaKf3iuF03ulwKPvqQVNnfGg/B5JSX6Vv/dIfLwak5SN/SjOcZ7Tv4UcF1f14bo6/8u+tCfsdlWmPjMz0iMBtIai0XQcfOh4WGmxN+R+OCsAtyhQAk1odEcIFepqtcKwYSpq7NAhOHiwfpOCCgtfeDJ7ge7i6DJ3RajUSlO0QNFoOgo+1PB6fEMpk1/fS11JGfz0kzKPNc2OB6XtBKsAYWqq2p/RqJzzoEKVnfxrfQl3fH0UW9MkR6NR90cpKfG5DUE40AJFo+koVFZ69XlIKflobyXDsxJIWLVCJS6OHes5Wc6flr+t4QodBsjLU6/r19evPqN3CuUWB6sLm9Qg6+iOeSlVgmcURHe50AJFo+kotOKQ31xcx4EqGz/rl+rd3AXBifBy4S6YPAiUU3ukEGfwED4cF6ei1jpqfxSzWbVEjnAyoztaoGg0HQG7XRUO9CJQPiqowihgRq8U+O47tdCTQHE59oMlUOLjG0xxo0ap91u31heATI0zcFLXZD7f30SguL7TUf0oFRVtD4oIMtE1Go1GExrq6rz6T6SUfFRQxeRuyXQ6dkg5xjMzYfjw5htbrSprPliOYJdAcTiU+WvQIDXz3ry5fpPTe6VQUGml0Gxr/v2OWnn42LGoMneBFigaTcfAFTLshQdP7syteVkN2smkSZ7NKRZLcCK83ElOVkIEGhzzbmav8welse7S/nRJbpLpEBenZuodjbq6iJeq94QWKBpNR6C62qutXQjBSd2SOaFrUkN2fEv+E6s1+FVtk5MbHPNjxqhXt0iv9HgjGQkexh8f3zEd8340SQsnWqBoNB2BViK8nt1UyqaiWmV2WrlSLXRlrjclmBFeLly5KNDgmF+3rpFW9dn+Km784gjSXdOKj1cP145WebisLKrChV1ogaLRdASqq1sUKEerbdy/sohvDplh5071sOrWDXr29LyvYEZ4uXAvHdKvn6pifOyYqqLrpKTWzscFVewq81CCv6P5UUpLo7JCgBYoGk17x2ZTM/gWIoJcWehn9E5p0E4mTvTsdHc41H6CPTt235/BoBpugdJSnEzupsrnf3/Eg7mnIwmUujrlx2raPTMKiBmBIoQoEEL8JIRYL4RY7WG9EEI8JoTYJYTYKIQYF4lxajRRRytNtT4/UE2vVBODMuMbCxRPuPwnwS71ERfXOGjAQ8Z8r7Q4eqWaVBfJpt+Nooq7ISdK/ScQQwLFyWlSyjwp5QQP62YDg5x/C4CnwjoyjSZa8RLhVWNz8N1hM9N7pyCkbBAoJ57oeV/BquHVlLi4hoKP0OCYd4v0Aji5ezI/Hq3B7l6GxZXg2FGIUv8JxJ5A8ca5wH+k4kcgUwjRLdKD0mgiTk1Ni+au3eUW4o1CNbNy+U+6dm3ZfxKKCC8XSUnNI702b26kYU3tmcyI7ARK6tyc8C4NpaOUsi8ri0r/CcSWQJHAUiHEGiHEAg/rewAH3D4fdC5rhBBigRBitRBi9fHjx0M0VI0miqiqanFGOzI7kbWX9lf+idb8JxDcDPmmuIcOp6erPixWq8qadzK7bxqvzu5JbpKb/8Cl2bi+256pq1N/Ueg/gdgSKKdIKcehTFs3CSGmBLITKeViKeUEKeWE3Nzc4I5Qo4lGzGaPAsUVfhtnEBgNonX/iYtQmVtSUhqSGwFGjlSvW7Y027TW1qR+l5St+oraBT5UjI4kMSNQpJSHnK/HgHeBplf9IaCX2+eezmUaTcfF4Wixhtfm4jpOf7uAjU3zT1ryn7gaYYVKQ0lIaFzo0VX2xa0EC6icmQmv7qXO3kSodASBUlYWtdoJxIhAEUKkCCHSXO+BmcCmJpu9D1zhjPY6CSiXUh4J81BDg8XSMW4WTfDxct18fqCaveVWuqeYYNeuBv9Jr16ev+Dq0hiqGXJToTdihHptoqH0SoujyurgpyK3opAdJdIriv0nEDstgLsA7wp1IZuAV6SUnwghbgCQUj4NfAScBewCzMAvIjTW4GGzqeQuV/e6AQMgOzuyY9LEFhZLiwLg60NmRuckkJNk8s1/YrUq30aocK86DKqDI8COHeo8nJrR+M7qgbq6sIYJXZzFETtCpJfFovwnycmRHkmLxIRAkVLuAcZ4WP6023sJ3BTOcYWU8nJ1I0mpwjSlVJ+7dIHevaNa7dVEES2Udi+vs7P+eC03j8lSC3zxn9hswevS6AnXNe0qj5+aCn37QkGB0qCcJrCcJBP90uNYfcwtmTEurv13b4yB5M2YMHl1OOx22LNHmRcyM1VRP5NJaSdFRbB7d6RHqIkVqqs9Tj6+O2zGIWFqj2T1AF+1Sq3w1O7XhZShNbcI0Th0GFr0o0zoksSawpqGul4Gg/K/tOdIr6qqqOt/0pToHl1HpaiokYrfiMxM1Ue6I1ZY1fiP2ezxOuqabOLCQemMyU2EvXvVNZWbC336eN9fqBPq3MvYQ4NAaeJHmTcgjZvzsrC6++WFaN++xoqKxjXPohBtN4k2rFbYv997NnJqqjIDjBoV9TMWTYSprvbo9xjfJYnxLv/Damclo/HjvTvcQxnh5SI5WQk3Fy2EDk/unszk7k18CVIqE1+oEi8jiZRKoGRkRHokXtFPo2jjiDMwzZuPJCFB2VPdbzyNpilWa4M/wo2iGht7yi0N5iKXQJngqaKRE7tdXZOh9t0lJjbOeHc55rdta6y5AMfMNjYcd/MrmExKI2uP1NV5/C2jDS1QoomaGjh82LdImrQ0paXYPLRE1WigRfPPu7sqOf3tfRSaneVLXALFm//EZgvPzL+pBpSZCT16qAnU3r2NVv3+h2Pc/vXRxt9tr5FeMeCQBy1QoovCQmWj9mUWYjIpJ6RbvwiNphEtFIX8+lA1gzLj6ZpiUhrxoUNqgjJoUMv7CmUNL3eaVh2GhnyUpo75zonsrbBSVGNr+G57zUWprIyJyE4tUKIFux2OH/fvpk1LUw+EjtatTuMbZnOzh5DZ6mDl0VoV3QUN2sm4cV5bBGO1hif/weX0dxcqLTjmXTkoqwuds3eDQd0L7THSq7w86h3yoAVK9FBV1dC8yFeMRvWdiorQjUsTu3gIGV5xtAaLQzKlZxOB4s1/4iLUDnloaC/sKdKriYYyKieBeKNgzTG3/iDtsaaXw+G142Y0oQVKtFBUFFhIZlJSgyNfo3GnurrZNfX1ITMJRsHEphFe0SJQoHnosHsJFrdaXwlGA6OzE1hV2MS/0N4ESoz4T0CHDUcHdrsSKJmZ/n83MVH1l66pUcJFowF1TVkszUyot+VlMatPCokmg7pudu1SphRXeK43wiVQkpIaa905OdC5sypDtG+f6jnv5N6TckmPdzPVGY3qXujUKTxjDQcxJFC0hhINuEpGBBoSaDBAcXHwxqOJfaxWj9dTp0QjJzl7s7NmjXodM8a7sLBa1cQlXCGrSUnNoxddZi+33igAI3MS6Z3upoW1R8d8eXnUdmhsihYo0cDx421zuKWmqmgvh6P1bTUdAw8RXt8fNrP4p9KGXiK+mrvCFeHlwlOk49Ch6nX79kaLHVLy8rYyvj7oFCImU1T3XA+IGMiQd6EFSqSx2ZTpoS3mKqNRmTi0c17jwkNRyLd3VvDvn0qJNzof1tEsUJoyZIh6bSJQDELw1IZS3tzpvPZNJmUiai/tgG22FvvZRCNaoESaiorgZMAmJuqcFE0DTbo0Sin57rCZyd2SMAihzEJbtqjJSF6e9305HOH1z/khUADG5CY2ZMwL0b7aAdfWRn12vDtaoESa4uLgqLNJScrW2t4iXDSB0STCa1eZhWM1dk5x5Z9s2KC02mHDWtc+QtlH3hMmU4NgcNGnj7pPDh9uVhh1TG4iB6psFLsSHIVoPwKlheTUaEULlEjicLTd3NUUXYVYA8qP4CZQvjus/Aondw8g/wTCnwORkNBYKJhMDZn8O3Y02jQvV5XU3+jewbG9TKwqK2PG3AVaoEQWszm4Bd+Sk7XZS6Ps7jZboyTZQrONgZnx9Ex1Ppx8FSiuZNtwl/1ISmpeAcJl9tq2rdHikdkJGAXsqXAKESFiKtTWK5WVMZHQ6ELnoUSSysrglp9PSFAaT21t5PpO22wNvcc1kcFD29+7T8jhzvHZDes3bFDvx4/3vq9wlVxpSnJy8yCTFiK9UuIMbLhsAKlxznspLq59VB2WUpkuA8lPixBaoESS4uLgP3iFUGavSAiUykqVKFdbq0KZu3dv6DipCR8t+A9MBqeQ2bJF/Ub9+0NWVuv7CmUf+ZbwlIvixTFfL0yg/ZSxd0Xqaae8plWs1tDU54mE2cvhUM7STZsaWhUD7Nyp7N0x5FRsFzQx9zyxoYSLPjqIzeFH/xMXNltkNBSXY94dl0DZsaOZOWxTUS1Xf3aIA5VWpaHU1MT+dech9Dva0QIlUpjNoZl5xMermymcM7RDh1SXyU6dGiLWEhLU7Le8XDcCCzdNIry+PlhNjc3RoKG4MuRbM3e5iERSnacy9pmZ0LWrEpj79zdaZRCCLw6YWXustiFCLNZ7BdXUxFxH1tgabXuirCx0jk6jUflSwkFVFRw8qISJp4tfNwILP245KGarg3XHa5nsKrficMDateq9LxpKONr+eqKlyKYWHPODO8WTZBJsKHLLR4n10OEYc8iDFiiRo6QkdH6OlBTVrCvUKr/DAbt3q+O1pG3FxSnzhK6IHB6kbCRQVhXWYHXA5O5OX93u3Woy07Wr6oToC5F4qJlMaoLS9BpuwTFvMghGZSc0JDi2hzL2WqBofKK2Vl3sodJQTCY1Owt1O9SjR32LKEtPV2ax9lZjKRqxuSX3AT8cqSHOACd4KlffmsnVZlPmrkiZXZKSmmsZrWTMbyquw+qQaswx6IOoxxX6HWMBLVqgRIJw9L2Oi1Ml8UOF2azs2BkZrW9rMKiZ1oEDoRuPRtEks3pYVjy/GJ5Jksl5q7sEii/+k3DX8GpK00Zb0KKGAqqD45icBEpq7bEf6RWjwlCHDUeC0tLQOzqTk1UV4969QzPL2b/fv9lrSoo6b6s1pjJ/Y44mOShzB6Qzd4Dzg5QNDnlf/CeRykFxkZyszHPu9OmjJieHDimTUFpa/aoz+6RyZp9U9cHiiH2BEoNRalpDCTdShidPxGV/dvVaCSaVlepG92f26nrI6YrIoaW2tl7IF5ptlNW5hdcePKh8WZmZMHBg6/uy2yMrUBITm2fLu5dg8aClgCqEGfMaSoyVXHGhBUq4qalRN0k47NIJCco5H0ykVKarQARiYqLSmjSho6qq3jf35IYSTn59b0P+yapV6nXCBN+uv3AXhWyKp74o4NWPsvCHY/z8o4MNE6pYjS6MQYc8aIESfsI5a0pKUppEMKNdKiuVhhXIzFVXRA49bhFePxypYXyXpIb8E5dAOeEE3/cXyVlyS8cePFi9NikSCZBiEmwocjrmITavNVfJFa2hBB8hRC8hxJdCiC1CiM1CiNs8bDNNCFEuhFjv/LsnEmP1iXD4T1y4ZnfBqkAspfKdtMVRK0RozHAa9fs4mzEdM9vYWWZhcje30j7+CBSX/T6aBcrOnc1WjchOxGKX7Cy1xG5flBgsueIi6gUKYAN+LaUcDpwE3CSEGO5hu2+llHnOv/vCO0QfCZf/xJ2UFFUWJRhUVCiTSlvGn5gIx44FZzyaxrjNxpcfUZrwZFe5+qNHlakyNbUhUsobLod8JB9qRqMyXTVtbe2uoTRxXI/KUZO1TcXOjPlY1FA8FPeMFaJeoEgpj0gp1zrfVwJbAR8zsqKM2trw+U9cuEqxtFUrkFJlvLc1jDQxUQmmWLzRox23B9EPR2pIjzcwIsupDbu0k/HjfYv6i3SElwtPocM5OSqwoLKymY+wb3ocqXEGNhfXxa5jPoZL70e9QHFHCNEXGAus8LB6khBigxDiYyHECC/7WCCEWC2EWH083A7iSF3cCQltz1QvKVGCKVjmOh3tFXys1voZ+42jO/Ho1K4Ym/pPfG2oZbUqbSbSeKo6LESDluKhx/zVIzIZ2zmxoUhkrFFVFZP+E4ghgSKESAXeBm6XUjZ9Gq0F+kgpxwD/At5raT9SysVSyglSygm5ubkhG69Hwuk/cSclpUEgBILDoXwnwXrAJCXpaK9Q4ObI7Zsez2m93LTJlSvV68SJvu1Lyshcq01JTvYcqeXFMf+rcdnMHZCuNDEtUMJKTAgUIUQcSpjkSynfabpeSlkhpaxyvv8IiBNC5IR5mN6RUkVcRarxlckUuO+iqEiZU4IVxugyezXNMdC0jepqMJlYebSGt3ZWNEQ6FRXB3r1KkI9oUXlvTjSErSYkeL5OvDjmAYpqbFRJg7pum/pgohkpm7VvjiWiXqAIIQTwHLBVSvlIC9t0dW6HEGIi6ryKwzdKH4iE/8SdtDRlb/bXd2GzKWeuW0Zy0IhF+3Y049RQXttezgOrijC5/Louc9fYsf49qKJBoAQQOryrzMKEV/eydF9V7FUddvnBtFM+ZJwMXA6c7hYWfJYQ4gYhxA3ObS4ANgkhNgCPARdLGWV1CyKtert6RBT7KWePHVNCJdiFLA2G8NQ06yjYbGC3Iw0Glh+t4aSuSQgRYP6Jzaa0yGgoTNhScqMrW3737mYmsX7pcSQaBT8VOcNvYym5McaDVaK+lpeU8jvAq7iWUj4OPB6eEQVISUnkZ3yuqr85Ob7NVKurle8kFD2tExOVT6lbt+DvuyPinIUXVFg5Um1j0pgA809APdQiWRTSnZau09RUVX7/0CHYtw8GDKhfZTQIhmcnsLm4FoiLLQ0lhiO8IDY0lNgnEvknnnDNOPfubb3wnN2u+sMnJ4fGTBcfr8I+tR8lODirDC8/ojThSa6GWiUlyiyUkACjR/u2L6s1NCbOQDAalYbiyQ/i0lI8+FFGZSewubgOB8RW5V6nHyxW0QIlHNTWRk9vg/R0ZfY6fpz8/Hz69u2LwWCgb9++5OfnN2x34IC6EUMlBF1mjEibAtsLtSqRb0eZhS7JRvqnO2f2P/6oXseN8z1qS0rlwI8WPIUOg1c/yojsBKptkr1mYstXV1UVeUtGG4hdURhLRPKCtliU2cr1Z7VC587kP/MMC157DbPzgb5v3z4WLFgAwPzZs1XeSlZWaMcmhLqBoiHfIdZxOuTvPSmTX43LavCfuATKSSf5vq9Itf1tiaQkz8USvQiUU7on89CpXchONcWOGcnVbdOXHkNRihYo4SDc+Sd2O6xYAR98AJ9+qh427quFgYU5vZE980iJS0DEJ4HDTvW2b1l4113MHzBAmTxCHWmSlKT+N127hvY4HQG3HJT0eDdNeIUzB3jSJP/2F20CpbS0+XIvJq/uqXFcOMjZfjpWtGCXrydGI7xAC5TQ48o/Cccs3G6H99+Hf/1LOStd9OpF8YChZPToginOxNOOHtBtEp2bfN2880f2Hz7M28VGNu0o5/yBaYzITmiY7Qab+Hj1v7Hbo8McGKs4cxfePCb4/EAR/5jaVXVoPHKkoVyOr/knrra/0fR7eOqLAtC/v/I37NunhEYTM93ucguHqqxMSbYqH0ykQvZ9JZZ8PS2gBUqoCVf/k2++gYcealD/e/WCOXNYfsrP+G9ZMp/tr+Kp07txRu9Uziyz8NAdf+TI7p9wWGqQlhowxuGoqaBPYiJ7jlaQv7uOF7aUMbRTPHefkMO0niGI+nH3o2izV+A4S64sO2Bmc3FdQ7tfl7lr4kTfHb2RbvvriZZCh+PjoW9fFTyyezeMHNlo9dMbS/nyYDWrZ2ep84qGzH9vNGnfHItEuchuB4Taf1JRAb/9LVx3nRIm3bvDX//KhvwlXDb4PC5Z7+DHozVcMSyTQZnKjDEQM4vOOwVT5RFspYexV5dirzhGMrCotpbf3HkBK0dUsmhyZ2rtkquWHua+FSEqlSJEM5Ocxk8sFhxS8uORmoboLgjMfxItNbzc8Rbi7sWPMiwrnqIaO8dqHbEROhzDJVdcaIESakpLQxcp9f33cM45sGSJOsZvfgOffIJ9zrnc/PUxtpRY+OOJOfzw87788cRc+qTHq4s2OZn5t93G4sWL6dOnD0II+vTuzeLbb2f+iSdCWRmZC65mfslmPp3Xm7snZHO6U0MJer5oYqIKbdUEjsXClnI75RZHQ/8TKQMTKHZ7dEV4gdKuWrruvJRgGe6stLylzBobAiVGm2q5o01eocTlPwl2TL/dDo89Bk8/rT7n5SEfeID3ZS5nGuNINAient6dPs5S3o2+Z7Mp04DJxPz585k/f37D+iNHYM8eZTpbsgQWLCDhn//khjPOqN/kgVVFWByShRNziTMEwbeSkKD+R7Fg445WzGaWFysfQ72Gsm+f6oHSqVPDQ9dXoskhD0qguFr6NjV9uRzzHjUUp0CpsDMtFiK9qqujJ/8nQPQdHEpqaoL/oCwthQULlDAxGOC22zj6zEtctj2e274+yju7VN+TEdkJjYUJKO2ke/eWHxi5ueoBf//9cNllalZ3663w9deA0k7sEl7cUs7lnxyivC4ISYmuB0Q7cEhGDLOZzCQTs/um0jXFOUd0aScnnujf9SdEdPoakpI8axleTF4ZCUZ6pJrYUh4DkV6u1gMxPqmK7dFHO9XVwQ0B3LYNzj8fvvtOzTyff55lZ1/B7A8OsfZ4LYsmd+biIektf9/hgOzsltebTMqZbzbDH/6g/DJ2O9xxB+zahRCCP56YyyNTurDmWA1XLj1MpSVIQiUWZpDRSnU1Fw7J5KnT3crYLF+uXv0JF7bblcklmiK8XHhqtAXQs6cSNseOKU23Cf+e3o17T8yJfoHSTiZUWqCEkpKS4M32vvoKLrlEhQOPGgXvvsu/U4dy7bIjdEuJ43/n9mb+0AwMLQmwmhqlTrdmH8/JUQ8Uux1+/Ws46ywlGG+4oT4X4LyB6Tx+Wjc2FdWy4PMjbfermEy64Vag2O1UVddiEW5CwOFoyD/xx38STTW8mtJStrzBAAMHqvce/CgjsxPJSU2I/glLO4jwAi1QQofdHpz+J1LCiy/CjTcqzeHssyE/H7p147SeyVw3MpN3z+nJgIxW7N61tcrc1RpGoyrY6NKu/vIXlcNw4IAyfznNDmf2SeXRaV25anhm2/NU4uNVrTON/1gsvLSnlrz83VRZnfWutmxRwr97d+jTx/d9RWOEl4ukpJbrvnnxoxyvsfHET2XsLquL7rpxNTUxXcPLhRYoocIVCtuWh63NBn/+MzzwgJp13nILq37zZx7YWImUksGdElg4MZcEYys/oytxMN2LOcydjIyGYnxJSfDkk8q/snKlCgZw8rN+aZzZRz2AtpXUBa6puPreR/MNH61YLCwvstI7zS0A49tv1eupp/p3/Tkc0Rfh5aKlXBSAIUPUqwcNpdYmeWhNMT8et0R3pFeMF4V0oQVKqCgpaVsIYHU13Hyz0kbi4uDvf+eVMy7nko8P8em+KsotfnShq65W5U18dfglJzeYvUB999FH1feffRbWrm20+bpjtZy1ZD/PbS7zfUyeiHazRBRSV1HFqiJr4/wTd4HiD9HS9tcT3kKHvZRg6ZlqIi3eoBzz0SxQzOaYDxkGLVBCg6uRVaCzvWPH4PLL4csvITMT6wsvck/2RH7/wzFO7p7MknN6kZngh+PUZlO+EV8RQm3v7sgcPx6uvVbNYu+6q1Ey4pjcBGb2TuGBVUV8f7gNiZzR7jiNQtbtLaLOAZO7O6+1igpYv149gP3xn7iIVoHi7WHrbvJqInSEEAzPSmBruS16BYrdrpzyWkPReMRsDrxc/fbt8POfw+bN0Ls3vPYaN5V35z9by1kwMpPnZ3Qnwx9h4nK0+uvLycpq7gS95RZlXti/X+WqODEIwcNTujIgI56bvjzCgcoAblztRwmI5fvKMQiY2MUpUJYvVw+osWP9y2mwWNQEKBojvKBxLkpTcnNVE7iKCjUZa8KwrAS2lduxm6N0wmK1xnRBSHe0QAkFFRWBXSDffKMiuY4cUQ+E11+Hfv24eEgGfz+1C7+fmIvR32TCmhro3LQMpA+kpKib1/0Gjo+Hv/1NzRZffVWFLztJjTOw+Ixu2B3w62+OYnf46U9JSNCRXv5it3NmFyN/ntS5YZIRqLnLYvHdxxYpEhM9axlCeM1HGZ4Vj01KjhRVhniAAdJOIrxAC5TQcOyY8kP4ipTKV3LDDcqUdNZZfL7oaV4+pn6e03ulcP6gAG92hyOw7FujUWkpTf0aQ4eqaC+AP/2pUfx83/R47puUyxm9A4gUMpnUjRWtZolopK6O4Rkm5g919s+QMnCBYrVGv0BpKXQYvEZ6zemfxuaLe9PTFKW95dtJDgpogRJ86urUQ9jX8hUWC9xzD9x3H9jtyBtu4Ikr/8C1Xx/nrZ2V2Pyd6bvjiu4K1JeTleX5Yv/FL1Ts//798NxzjVbNG5jOglGd/NekXGjHvM/sOlzGl0ctWOzOa2T3blVuJSdHCX5/iVb/iQtfBIoHx3yiyUBcfFz0+ujaQQ0vF1qgBJuqKt+3LS6Gq66CN96AhASq/vp3bhx9EQ+tLWFO/zRend0DU1vqZdXUqMz4QO2zaWmeVfG4OCUEQZWAOXiw2Saf76/isk8OUmf3IxpNiOhu1yolFBaqOlm7d6u6Zy094MLAG2sPcf2KCuyu38ilnZx8sv8lPKK15Io73nJRvJi8AJ7aVMbffqqMztD06uroq58WIFqgBJuiIt8c4KtXw9y5sGYNdO6M9T8vc4EcxdL91fxhYg7/nNqloa9FoFitbWvjGx+vTHcWS/N1J56okizr6lTyYxMk8N3hGp7Y4KHTXkskJES3Y95VPLO4WLWkLSpSgsXhh9AMIt/vK2dsTnzDdRKoucvVVCvao4x8ifTatcuj0NhaUseSg3XRaVI1m6P/f+8jWqAEE4ul9ex4hwMWL4YrrlC+lvHj4e23icsbzYWD0/nPzB5cO7JT8LoktrWURtPwYXd++1u1/88/ry8g6eKM3qnM6Z/K0xtL2VvuQSB5IiFBPaij0UFZUaE0k06dVDZ5crKKLCorU10Rwzzm0moLW4rrOLm78/c1m2HVKqVpnHyyfzurq4uNKrfeBEp6usqXqqtTVR2aMLRTAofMDioqo8zsZY2RbpI+0j7OIlooKlI3dEvCoLBQFVz8+9/Bbqf2uuv5/S8f4WuLeihcM6ITp/Tww5nvDddDoq0zn9TUlmfgXbqoUGKAv/61mfnnDxNzSTAK7vnxuG9Z9AaDml160ogiSV2dMqWkpTW/8TMz1e965EhYh7R85zEkcHIPp0D57jv1fxszxn+tNBYc8tC6n8FLbxRXKfvtR6JMA7ZY2k3IMGiBEjzsdlW40dNMT0rVX+Tss9WNn5nJun8+z1n9zuOVHZVsLg5BlEdNjX/JjC3RmkN//nxVoXj3bnjnnUarOieb+PX4bL49ZObHoz7ODKWMrqgXKZWZy2DwbOcWQmkt+/eH1f+zak8RKSbB6FynNvzZZ+rVrXeNz0gZuiZwwcRoVBOk1mp6bd/ebNXQLPXbbTscZaHp0TZ5aiNaoAQLV5OopolhBw6owo6//S1UVFB3+hk89NCbnH+0M7U2Sf6sHvxyTBv8HC0hZXDMGCaTeti0dOHHx8OvfqXeP/ZYs4fq5UMzeGlmd07q6mOkmcEQXY75ykrl1/FWNNElbMKopfxxWm8+Pj1TNTmzWlU1aghMoAgRGwIFvEd6eXHMd002MTTDhCPaHuC1te3G3AVaoAQHKVWkk7u/orZWPWDPOkuVUElJgb/8hWW/eZAndtVxwcB0PpnXm5O7B8nE5Y7d3iAIgkGnTt5nUrNmqZL6x4+ryshuGA2CqT1TEEJgtvrgvI6Pj64Ex6NHfYt+SklR5x+msGeD2UzvDOe4Vq9W/7MBA6BfP/92ZLMpU1KsOIVbarQFDaHSHjQUIQSfnNODK/tHmeCsqoqd/70PaIESDKqqlIkpPl49ePPzYeZMeOIJsFjYfuGVLH3xAzj/fGb3S+Wtn/Xkb6d2IT0+RGUuamqUHT1Yttn0dO/hsQaD6mcP8MwzKgqqCZ8UVHHS63tbL8vicsxHA3V1qsinL4ENQqgHQ2FhyIf16eaj/H7pXqpxXj8uc9eMGf7vzGKJDYe8i+Tklq/F/v2VhWD/fs+BJCZT9OU5tZOikC5iRqAIIWYJIbYLIXYJIe72sD5BCPG6c/0KIUTfsA3uyBE1a3r5ZXVT33cfFBaydvKZXHv/u5zZ/0Lu3WrB6pAYhGBClwATDX3FalXO4mDhS2LkiSfCtGnqBnnyyWarx+QmYLFL/ra6yPt+DAb1wIgG00RxsX+1rVJTlUYT4rF//NMRlu6rIjk5QWnHy5apFdOn+78zi0W1K4gVEhJajqiLj1camsOhwoeb8N2RWqZ9XMT+41EyYXE4lIDTAiW8CCGMwBPAbGA4cIkQYniTza4BSqWUA4F/AH8NxVjy8/Pp27cvBoOBvn36kP/738Pddyuzz5//DEeP8sPkszj3D29z3qm3sNqSwO1js/jw3N7K3h0OhPCv9EtrJCQ0LmffEr/+tTr266/D4cONVnVLieO6kZ34YG8V6461MkuMhpbADoc6B3/Crg0GNfbjx0M2LCklP+wuZlKuCWE0wqZNSivq0gVGjgxkh7HjP4HWH76u3ige/Cip8QYKqu1sOehHblQoaUdFIV3EhEABJgK7pJR7pJQW4DXg3CbbnAu85Hz/FjBdBC2ZQ5Gfn8+CBQvYt28fxCVxqLyOX77wJnfuKOOeyVewcsYF8NhjxP3lfmzp6fzxxBy+/3k/bh+bTafEMFVxddnEg531nJnZ+kN+8GD42c/UjfLUU81WXz+6EzmJRv6yqpUw4mjImK+oaPBF+UNamhJEIcqg33msimOVdZyS64w4c2knZ5wRuHO3PQkUl2Pegx9lcGY8giiK9GpHRSFdxIo3qAfgnq10EDixpW2klDYhRDmQDbRiY/GdhQsXYjabMaZm0fOm/9QvfwtINsDwSZ2ZOCSDE4APz41Qb+7a2rZlx7dERobKs2ltxn7TTfDRRyqEeMECFVLsJDXOwB3jsvnDD8fYVFzHqJwWHmQux3zXrkE8AT85fDiwB63RqLSb8nJV9ibIfLNDaT+ndPEgUPzFYlGabCyZXNqgoSTHGeiXamRbYZSYvKLBrBtkYkWgBA0hxAJgAUDv3r39+u7+/fsBsJsrKPniWRw1lThqK7GXH8P8zSeBF0QMJjZbcP0nLnw1ofXvr0rKvPOOCkp48MFGqy8anM6ELokM6eRFg4q0Y762Vh2/U6fAvp+UpMxQIRAoBiGY3C2JHplJKvdn1y4VNHHCCf7vrLZW9Z2PJVwhzjabZ+3Ri4YCMDTDxJbCao/rwk47c8hD7Ji8DgG93D73dC7zuI0QwgRkAM3CjaSUi6WUE6SUE3Jzc/0aRL0ActioXPUe1Zs+p2bXSrqZ6jDWRUn0iJTB9Z+4SExUN7MvKvovf6lu9iVL1EPPDZNB1AuTFsOII+2Yb6swS0xU+wiBH+jqU/rxyikZSui+/75aOGNGYA8muz22IrxceAsd7t5dBUeUlCiNugmn9UhicrckHG2p4h0s2lGVYRexIlBWAYOEEP2EEPHAxcD7TbZ5H7jS+f4C4AvpU70P31m0aBHJTR7WycnJLLr/fhWmGKEigfVYrepBE4rKpQaDMnv5ksXeqxecf776fzzxhMdNHllbzNlL9mNt6caOpGPe3342nhACSoPr/K212pG1tUoQCAEffKBWnNvUnegjwQ7eCBfeqg67N9vyoKVcODiDv0zIwBAN1oTq6naVgwIxIlCklDbgZuBTYCvwhpRysxDiPiHEHOdmzwHZQohdwK+AZqHFbWX+/PksXryYPn36IISgT58+LF68mPlXX63a9Ua6Um5dXWj8Jy58ccy7uPFGNfv66KNmWgrA6JwE9lRYeWeXFwdpJPpXWCwqr6itQQ0pKcrsFcQ5zaOf7+SUf3ynhPDatarUT9eugZm76urUGKO15a83kpO9Vw324kfBZMJRU4PZEuFmW3Z7Q7+idkRMCBQAKeVHUsrBUsoBUspFzmX3SCnfd76vlVJeKKUcKKWcKKXcE4pxzJ8/n4KCAhwOBwUFBcyfP1+t6NpV3aCRjE6yWkObU5CU5PsDsls3uOACtb2HvJTpvVIYk5PAY+tKPPdMiVRL4Oog2dfj4pTwDdb+gG93HqdHiok4k7HB3HXOOYFFd7l65cQicXHew229aChSCE78XxF/+3hriAbnI+0wwgtiSKBEPQaDKn1RVxdZ01eg3Rl9weVH8ZUFC9TN/+GHzbQUIQS/GpfNoWobb+zwIDgiVYKlqCh4IddxcR6rBgRCcVUdmw9XcGrnePUbfPKJWjFnjvcvesNbfbJoJi7O+8PYi4YihKBnsoFtRyIcOtwOI7xAC5TgkpSkNBV/ujYGC6tVPfBD2fktPt57tdemdO+ufCktaClTeiQzoXMiL2wua56XYjSG3zFvtytnbrCEckqK8scEoUvg97uLkRJOzZTw44/KvDp0aMNs3B+kVH+hnHyEEl81lJ07PeYDDc0wsa2wyreWCqGinZWtd6EFSrDp3DkybUZrawMPc/WH9HT/ystff32LvhQhBA+e0oXXf9az5YZi4XTMV1d772fjLwaD0laDMMH4budxMhJNjOpkanDGB6qd1NUp02is2u9NJu8Rh2lp0KOHemjv29ds9ZB0I2U1No5VRrBNQjsrCulCC5Rgk5SkbtZw+1JstvDUZMrI8E9r6N4dzjtPPVg9ZM8PzIwnN8mElBJH0weEwRBUH0SrlJQE/yZPSFBaShuZM6YHd0/tg7GyEr74Qj1Qzz47sJ2Fa/IRSryFDoNXP8qQTGdvlKMRzHUymyPXR76mRvVlCgFaoISC7t0jE/IajhBQfxzzLlxayocfqmZVTSiusXHBhwd5d1eTGzycPealVP6TYJuBkpJU+HAbe5mfMiiHS/onKd+J1QqTJqn6XYEgZez6T1x464sCDX4UDwJleHYCt43LpmenCJr8IpXUWFEBV18Nd97ZrCFeMNACJRSkpXlvShVsLBZ1g4XjAg2kHEmPHl61lKxEI2abgyc2lGB3z0txRXqFw9ZtNocmjNNlPmtDgMGafaX8dLBcCaY331QLXdGF/mKzqZlxLOafuONr6LAHgZKRksAdw1MZkBshoerqIx9uH0pREVx+uQo579w5sGKiraAFSigQAnr2DJ+5prY2fCGgJpN60Ptb/HDBAvXd//0P9u5ttEoIwc1jsthTYeXDgir3FeFrCVxREbobPDGxTX1S/r50O3e+uV41atu/X2nAp50W2M6qqlTgSKw7hJOSvEdTDhumXrd6CA82mSivMCshHQki4ZA/eBAuvRS2bYO+fWHx4sACOlpBC5RQkZmpLppwhBCHu4SGv455UALWi5Yyu28qAzPjeXJDSfPom3CYD4uKQjdrd5ViCUAwVtXZWFVQwrT+mfDWW2rhJZcErkk5HKFNfg0X8fHeH8q9e6vf8+hR5Rtzx2Dg8S1VnP/0D9g85UCFmjaaP/2mqAiuvFIFKAwfDq+8ErLCq1qghAqTSf1o4dBSwl1CIz09MHPe9der/8sHHzTTUgxCcOOoTmwrtfDVQbeABpMp9GHYVmtDx81QUlbm91e+31WE1S6ZJkth+XI1xgsuCOz4ZrNyxge7tUEkaC0XxWj0avYakmHCYnNQUByBROSamvBpKGazuu8OHlQmrv/8J6TWDC1QQkl2dsj6YtTjKkEezhDEQB3XPXvCvHktailzBqTxl8mdOamb2/4TEgJ6EPtFdXXo/TQpKWq27CdfbT9OaoKJCW89p8Z49tmBaxh1dZFtCRBM4uIawrJbwtVjfsuW5qvSlYa3PRKRXuEqCmmzwR13qCZsPXvCv/8dcktGqwJFCPGQECJCzT1inORk9UAMpYobiRIabWnI5K6lNMlLiTMILh2aQZLJ7bKMi1M3YChNh2Vlob/B4+PVb+VnOPmKPcWc3DuNuCXvqQWXXRbY8W02NWuPxerCLeGtvzx49aMMTDNiELD9aAQy5sMV4bVoEXz1lTK/P/ss5OSE/JC+aCgDgZ1CiCtb3VLTGCFCb/ZyOMIfAmo0Ki0lELOXD5WI39tdwb0/OnM3XKaBUPlRpFTlUcKRNR4X53dOyoe3nsq9B75SZr+xY2HEiMCOXV2tnPmBdnWMRpKTvV+DLoGybVuzVYkJcfTNSAh/LoqUamIRaovChx8qX0l8vLIG9OsX2uM5afXqklLOA64A7hRC/CiECKC0aQcmMzN05hTXfiMRAupvgqM77pWIPdRbKqiw8uKWcnaUujmxQyVQamoaZu+hxlWKxQ+NNclSQ7fHHlYfbrghsONK2X6c8e6kpnrXUAYPVgJ0z57m14/JxJ9OyuHW6YNCO8amuCK8QulD2b8f/vhH9f53v4Nx40J3rCb4NF2RUi4D8oBXgI+FEM8LIQLMqupgJCYqM0MoHogWi3KQR6KERlpa4Ka8bt3goovUg+5f/2q2+sphmSSZBE//5OwnEhcXukKR4czENxjUOfvYJ+X/lmzilYdeVtuPGQNTpwZ23NJSlQsUS73jfaG14ILERNVB1G5Xdb3cMZk4tZNgZI8wVJdwJ9QRXhaL8ptUV8OZZ6qIwDDis/4rpbRLKR8DzgSmANuFEL8WQsRoQaAw0qVLaHp7mM2RK0Hub+Xhplx/vXogLF0Kmzc3WtUp0cjFgzN4f3clh6usoc2YD0V2vDdSU1Ufk1a01qo6G6+s2M++FRvUgttuC+z/XVmptMlYa/XrC61FekGDY76pH8VkoqKqhvfXH+JQWRj77oS6bP0jjygnfI8ecP/9Yc938cUpP0MIsVAIsUQIcQRYAdQAbwCzgHVCiDDrjTFGRobv7XP9JVJO1sTEhqq1gdC5c0O292OPNVt9zchMJPD85rKG3iLBjpiz25XmE84w2rg4FW3Visb19fbjWB2S0zd/q0wWkyerh5Hrz5eZrsskOWBA+/KduPAlzHv4cPXaVKAIQUmd5NbX1vPdzuPBH1tLmM2h85+sXQsvvqgsFo88oqwXYcaXq+xNlEayDtViN0tKOUpKuUBKOQN4HHg9hGOMfUwmZb8OpnnFVUIjUmYMg0H5BNqiwl93nfL/fPUVrFrVaFXP1DhuG5vFxK5u2kOwtTzX7xHurOWkJDhyxOsmn63fT6eaSsYf2go336yS86RU/3eDQf3+paXqr7a2sWB3OJRmUlUFgwZFrghhqDEalYD2Vt3bS6RX7xQDSXGG8DrmQxUybLHAPfeo6+DaayEvL/jH8AFfBEonKeWZUsp7pZRLpZSNplZSysVA/9AMrx3RuXNw7ac1NZCbG9kSGr72mG+JrCy45hr1/m9/a6bt3JqXzcw+zgg2g0E9JINJWVlkSognJaljewghzs/Pp2+//ryztoDSXSt5vXcvlZA2dCiMGqVm3MOHK5/KmDHKRxAXp0yCLgFTWalMoaNGRWSWGlZaq+nlMnlt394s9NwgBINzksObixKqkOHnn1d+ot694Ze/DP7+fcSXKC9fbBozgjCW9k1qqnp4BctsE65y9d5oLcrGF37xCxUfv3FjQxdCN0pr7Ty1sYS6+MTmJTTaQjjDhT0RHw8FBY0ecvn5+SxYsIADx0qpPbiFY9u/Z8HRo+Rv2eK53HxioppUDBum+sqPGaOEzbhxql5TSgdIH0tJ8X4NZmUpH6bZrKKfmjAkOzF8AsVuV8Iv2EE0BQUNIfj33RfR4IugGFallKta36qDYzAELydFSqWZRPqB0VbHPKhzuOUW9f4f/2gWirylpI6/ri7mvQO16qEQLIFcW6uOFakmR6mpyo/iZvpauHAhZrMZh7mM42/fR83ulZjr6lj4f//X+v6EaIgojNXGWYHgi9m1JT+KycSQdBPF1RaOh6PZVigc8lLCvfeqfc+dq9oaRJB26KmLYrKygtPNsaZG7SvSjlbXTKitN8kFF6jEq3374I03Gq2a3C2JEdkJ/PunUhwOR/D8KJWVka+4m5kJBw7Um/L2O2fQxpTG2sh+DzNrjRNfzEctRXrFxTGvdyLf3306Oalh8DOFop3F55+rGm+ZmXDXXcHfv59ogRJOkpLUjKqt5djr6iIXLuyOEGqm3dYbxWRSDX8AHn+8UYiwEILrR3ViT7mVZYW24PlRjh+PfE91V2DDzp1QXU3vrl2Jy+lDz5v/S/KQk+s36927dwQHGeW0VnUYGjSUJuHpmExkSQs9MpNabkEdTGprgzuJsVrhoYfU+1tuiYrEVS1Qwk23bsGZZUfa3OWirY55F9OnKz9AaSk8+mijVWf1TaVnqonFu2qD40exWlUEVDRU3U1IUA+ZlStZZLWSMehEAOoOqoKGycnJLFq0KJIjjG58yUVxNZLatKnxtkYjWCy88mMBb6w6ELoxugh2hNcbbyj/Sd++KlE4CtACJdy4HOmBmoksFiVMoiUUNDU1OGY8IVS5CKMRXn21kXnCZBBcO7ITGYkmzBXVbfejhDM73heSk+HBB5lfVMTgoadA0V4c5jL69OnD4sWLmR9od8aOgMGgTK/e/CjduimNvqxMlXF3Rwg+2HCY/BX7QjpMILgCpbKyocrEnXdGpp2wB7RACTdxcW3LSamuDryXeChwzbCDwZAhKtnR4YA//alRBNSVwzJ4bkZ3ko20XcMrKYmaGxCAp5+GpUsp7NyLQ53785uLZ+BwOCgoKNDCxBdaqzoshAqhBvjpp2arh+Qks6OwCocjyA5zd1xFIYN13T3zjNLmx42DM84Izj6DgBYokaBr18D8Dq4+1JmZQR9SwLgESrDKy996qwojXrcOliypX+yyce83Sw4fLg58/1IqgRItPdU//VSZ+ITgk9vuA2DmqHZYJiWU+BLp1ZJAkZKh2YnUWO0cKA1hsy3X+IIx+SosVBnxAHffHfngEje0QIkEqakq4czfmXZ1tco7iKbZtRCBd3D0RFoa/OY36v3f/tbIZ1JtdTDr8xL+9W0bzBNmszLRRTpCDlS15TvuUELuttv42SlD+MecIQzq0o56loSD5OTWJzQuP0pTgWI0MiRDhY6HNB8lmBFezzyj/JYzZ6rcoygiCu6qDkrPnv4LFKtVZdxHG4H0mPfGuefCxIlKmPzpT/X+ppQ4A3MHpPH2niqOlQZoMiwtjQ5h8u678OtfK+F2/fVw9dXkpCcxb9KASI8s9khIaN0n6dJQNm9u7PMzmRicJDEZBEcrQtQiAYKXg1JYCK87K13ddFPb9xdkouDO6qCkpfnXpKquTmk20RLd5U5KSnATtoSAv/xFzTw/+UQ1C3KyYFQnbA54/uudXnbQAg6HasMb7oZk7kgJzz2n+lQ4HCrc8447WLqrlFeL4gilGb/d4ku0Xna2qrhsNqv+KC5MJlKstWy+70yumNQ3ZEMMmkP+2WfVM2PmzIb8migiqgWKs/3wNiHERiHEu0KIzBa2KxBC/CSEWC+EWB3mYQaGEKp7YVWVb9ubzaokdTSSlBT8DOBevZR9GFQ5icJCAPqmxzO7VxL5a45QUetnbbSqqvA10/JEWZlqLuaqW/brX6vCj8DinTW89FMRBkP02MNjBoPBt8mZJz9KXBzU1JBgDPGjsLq67VUZjh2Lau0EolygAJ8BI6WUo4EdwO+8bHualDJPSjkhPEMLApmZKvy3NYeiy+YfrYX+4uLULDHY5eV//nOYMkUlOv7hD/VC68a8bOrsDtbs9rPs+PHjkQu3XrkS5s2DL79Uv+MTT8CCBQAcOV7B6hIbZ4/RzviASUtrXaB48qM420p8ufkIlz+3gjpbEELgPREMDeXZZ5WlYsaMqNROIMoFirO6sesp9SPQM5LjCToGA/Tpo2o6eZvhV1Yq7SSaazRlZAS/K6UQqklQRgZ8840KrwVG5iaxYlYWp3X1IzHRZlPFIMNtMjx6FH71K7j8cjh8GEaPVv4Tt1DPj/Ypf9BZo7qFd2ztCV86iLo0lE2bGi8XgipzHd/uLGL3sRDkKNlsbdeMi4rgtdfU+yjVTiDKBUoTrgY+bmGdBJYKIdYIIRZ424kQYoEQYrUQYvXx42FsrNMS2dkq8aqljoRVVepm6do1vOPyl8zM0LQ37dJFmYiEgH/+U9UuAjplpkBhIeVmH49ZXt5QVDMcuDL+Z81SPqCEBOUvyc9XARkupOTDQ3UM75ZG/9wI+nZiHV8q7Lo0lK1bG2szUjIkU2muOwpDEOnl6iPfFv77X6WdnH56Q4+XKCTiAkUIsUwIscnD37lu2ywEbEB+C7s5RUo5DpgN3CSEmNLS8aSUi6WUE6SUE3Jzc4N6LgHTu7cSGk3rVNXUqAtx4MDoiEzyRjAqD7fEtGkqvBZUVvDOnZCYyP99d4RzH/8Wuy+e7MLC8NTuKiyEv/5V3fhPPql+wzPPhI8/Vv6SJia32iozdqOJn43W5q424YtjPi1NFSG1WmHHjoblRiP9EiVxRhGaZlttDRmuqoJXXlHvr7uu7eMJIRGq3d2AlNJrmqcQ4irgbGB6S71ZpJSHnK/HhBDvAhOBb4I81NBhMKg2rZs2qZmtyaQezlarUtOjKe+kJVwzRIcjNMJvwQLVJOnDD1UDoVdeYVJOHC/tqeR/Gw9zbp6XgIW6OiWsPfUUCQZSwurV8PLLsGxZgy/p1FPhhhtgQstuvUS7lSW/GIsjM0Rj6yjExan7xm73bloaNQr27lV+FJfGEhdHXF0NA3JT2X7Ue2vmgKira9tk6623lFl83Dj1F8VE9bRXCDEL+C0wR0rpMY1VCJEihEhzvQdmAps8bRvVxMcrVbZ/f2XeyspSjrdoyehuDSGCVyiypf0vWgQjRqhGSZddxsx0C4PSTTzx5S7vZTOOHQuN9lRSojKWzz4bLrtMhThLqTSSt99WTlQvwkRKSYXFDqmpOrorGKSmtm529eSYj4sDs5mT+meTmRyCoI3KysAnhVZrQ1b8tdcGbUihIuIaSis8DiQAnzlLb/wopbxBCNEdeFZKeRbQBXjXud4EvCKlbN76LxZISop8SfW2kJmpepqE6hySktRD+he/gG3bMFx9NTf/8V/cttvG0i2FzBrpwc9kNsOhQ8HTTmw2+P57eOcd5c9xPcByclRU2kUX+ezvWnu4ikuXlfJCtwomD8gJzvg6Mi6zsTd/iiuzfP36hmUmE1RWcu/ZJ4RGu66qCjy68KOPVBO2/v3htNOCO64QENUCRUo5sIXlh4GznO/3ANFVf6CjEuwER09kZcFLL6le9Js2cfa9v+S++f/ghofz2f+f39C7d28WLVqkiipKqcwbwfDv7N6torOWLFEaD6iHz9SpcOGFys/j5yz0rR1lGIRgdM/Mto1No0hJab0Ey/Dh6nrYvVtpmO49RCyW4LfPtduV1h6IpUFKNYECdb1Hux+VKBcomhgjXNpVZqYyA1x3Ha+tW8fOF++iqrwQKSX79u1jgTO/Y/6MGWrGGmjjoZIS5Ux/7z3V895Fnz5w3nmq5WqA0Xe1Ngf/O1DH7OFdSE3Qt2FQ8MUxHx+vQrdXroS1axvCt4WgvLKGC55cwS9O7selJwapqVlb/CfffaeCB3JzYc6c4IwnxOgrWRM8jEY1E7NYQp9AmJYGL77IwsmTqTxe4FwoAInZbGbh737H/EGD/E8GLS9XyYcffqhMW666TykpMHu2Sk4cP77NGs+neyuotEkumNinTfvRuBEfr2b1rYWHjx+vBMqaNQ0CRUrSsVFYUcvmwy2E8AdCW2p4vfSSer388ujpf9QKWqBogkunTip0Nhw3QHw8+80qVsPUqTu55/2B0i+epXbvWvYfcHbga63chd2uZoErVypBsmpVQ5SW0ahMWuecox48QdTA3tpRTs+0OE7qHwWtnNsLBoOa0Fit3q+/8ePV65o1Dcvi4hBmM0O7pge36rDZHFhC4+7d8O23ygT3858HbzwhRgsUTXBJS1NO8FBjs0FZGb27d2ffoUPYKo5hiEsg89TLOLp3Lb1BCYPRo2HsWCXoEhOVn6OoSGWtHzigTFnu9dSMRjjpJFV8b/bskPXp/r9RyRzN7amju4JNWpoyVXoTKGPHKuGzebPKE0pKqo/0GtI1jffWH0JKGZw+85WVgU2uXNrJueeGLtw9BGiBogku4QhztliUEBgyhEV//SsLFizAbDZT9v1r5Jx1G5lDJrOoeIcSHCtXqj9v9OypZq2TJyvnehgamA1Mj2PgiCgt9hnLpKXVFxJtkdRUFZK/ZQts2KAmEHFxUF3NkK45VNbaOFxeS4/MIGikVVX+3xOlpQ3N5a68su1jCCNaoGiCS1ycumHr6nxzkvqLw6Fu0uHDIT29vkXuwoUL2b/5CzjlYgb//DdcMquz2nbjRjUTrapSY6qrU+VuevRQJW9GjAhrS+U6u4PffXOUa4anM6Kt1Wc1zfHVLDl+vBIoa9YogWIwgM3G2B5pnDeuB3Z7EKIVrdbWEy098cYbqi7eqaeqhOcYQl/RmuCTk6OSD0MhUMrLVZSVm7N9/vz59YJlyfpD3Pbaej5zZHFmXLnyfURRz+0P9lTxzt5qzs/ThSBDQlKSEg6tVWwYP17Vx1q7tmGZEIzITuSRn+cFZyyBlFyxWlXFBYCrrgrOOMKIFiia4JOW1rgrXrBw2bu9aBTnjO5OvNHAGUNzYeMG5WuJEk1ASskLW8oYnGZk8nBduyskCKFMlmazd23F5Zhfu7bxNWK1qgoGtTYyktpY8qiuzv8Ir08/VXlOAwfCySe37fgRIPozZTSxR1JSQ12lYCGlEigDBnideRoMgtmjumGMM+Ho0aN5wc0Isqqwls3FdVw1IAkRKyV1YpHMzNZbKXTurIqyms2qRpyL2lpufmUdP396edvHUVXlf8mV//xHvV5+efgqYwcRLVA0wUcIFR1VUxO8fZaVKb+Hj/1MPtl0lJn526h0iNBoSwHw4pYyMuINzB2eE929bWIdX3veNA0fdjrme2cns/t4FRZbK1n3reFvyZUNG9RfRoaK7opBtEDRhIasrOD1R7FY1M3e3XczUffMRHYdr+bJQ4ao0FKklAzLSuD6QUkkd9F1u0KKq9ROa+Yml0BZ7ewa7hQoQ7umYXNI9hT52J7bE1L6r6G4tJMLL4zZmn5aoGhCQzA7I1ZVKUe8H7P60T0zOW9sD55bf5wD1faIaylCCG7Jy+KXg5NUFJwmdBgMvlW+PuEE9bpihbo+XLkondXvs+1IGyYiLoe8r2arwkJVrdpoBGeASSyiBYomNLiHD7cFV2G9ABIMfzNrCAYBD+62N05eDDOHq6z8b08lDrtdPexidPYZU/jiR+nTR+UglZWpXkTOh3//dBMmQxubbflbcuXVV1VwwBln+KWJRxtaoGhCR05O2/0o1dXQt29ADspuGUlcP2UAH+4qZ1d5G7vmtYF/bSjhjm+OcrTUrB50MVA1NuZJTW39gS4ETHE2d/322/rF8VYLd88eytTBbejoWlvr+zVbVwevv67eX3FF4MeMAvSVrQkdaWmtlxP3RnW1egD7W+DRjRumDiD/2hMZ2CU9dM2/vLC/wsqbOyq4ZEgG3Y22mCqjEdO4ytC3JlROPVW9ugSKyQTV1Vx7an8mDWhDnTV/mmp98IEqFzN8eINfJ0bRAkUTOpKTlTAIREuRUgmAXr3aNISkeCMnD8yBbt0oKQ2/2evR9cUYDYKbxjhNdsH0LWlaxmhU115ryYUnnqge/Bs2qJIn8fFQUUGt1c66/aWYLbbAjl9R4Vtir5QNdbuuvDImQ4Xd0QJFE1p69Gjdlu2JykqVKxCkB/DHB2qY/GkJ20rCp6XsKrPw7u5KLh+aQZdEg3rIBbuBk6ZlOnVqfTKTkqK0Ainhhx+UQKmqYsWeYuY9+QMbDwZQyt5qVZMhXxJqf/ihoefJWWf5f6woQwsUTWhJS1M3qT8hxDabMpX17Bm0YZw0uAspcUbu+vYIdm/954NIucXOgIx4bhjdSQnVzMyYn4HGFL6aXN39KM7fZ1iW0i62Hqnw/7g1Nb7/zq5+8ZddFjM9T7yhBYomtBgMKmrFnyirykrliA/iDdYpJZ57ZvRnQ7GVpzaWBm2/3hjfOYlP5/UmJ8mkZqxhqGKscSM5WWmErZm93P0oDgdISW6cJCslPrDQ4epq3wTKrl3wzTdqjBdd5P9xohAtUDShJytL3WC+zBZra9UNlhP85L85E/sxp2cCj6wrZvkRc9D376K4xsbj60uoszswuB4sQmj/SbgRQlWUNrfyWw8apOrDFRXBtm1gNDqbbaWx7WgAGkpZmW+mTZfvZO7cdhOsoQWKJvTExSl/SHW19+2kVNv07x+S0FoRF8dfpvdmYJqJXWWhCyO+b0URj64v5kCl06Hrqnyr/SfhJyOj9YmMEI21lPh4qKxkWLd0thdW+mcidTh8a6pVUgLvvafex1jPE29ogaIJD126KN+IzUvUTEmJ8pukpYVsGKk9uvK/0zK4fFhmSPb/xYFqluyp5KYxWQzMdD5UtP8kciQmquuptZBxlx/l66/rI70uOaEXL/5ion/Hc1UYbu23zs9XprjTTlMTqHaCFiia8JCUpEwL5eWey6CUliqhE0RHvEdSUog3GkBKvjhQzV9WHkf6W2K8BXaW1nHb10cZ2imeG0e7mTDq6tqNSSMm6dq19WivyZOVIFm7VpWPl5KBmfGc1D8boz9tmn0Jka+qUr1YAK65xvd9xwBaoGjCR1aWEiplZQ1CxW5Xnzt1Cjgj3i+MRjUOs5nVhTUs3lTGfSuK2ixUHFJy05dHSTQKnpvRnQRjk1tLl6uPHOnpSmvw9hunpSltQUr43//Usro6PttSyA+7i3w/VkVF6wmNr72mJlbjxzfUE2snREfnIU3HISdHCZG9e5XwiItTLXn79g1fSZKcHCgu5jfjs6mxqaZXiSbBb8dnIwIUaAYheGRqF6SEHqluDxRXC1jtP4kccXFqElFV5T0wYs4c1eDq/ffhggugpoa/fbKNPtkpTB7gY5BIWZn3hMbaWnjhBfX++ut9PoVYQQsUTfjp0kVpJEZjZPqCpKaCEAjgnhNzsDgkT20spdLi4L5JuQ2RWT5wvMbGJwVVXD4sk5HZHoSG9p9EB507qygubwJlyhT1W23fDvv2QadODO2Wzrr9PoaZW63KL+JNG337bTWOESMa/DbtCG3y0kSG+PjINZkymZQZpLYWIQR/npTLgpGZOKT0S5isOFrDz97bz6KVRRyqaiFx02IJqFKyJsikpyth4S3SMD4eZs1S7z/5BCoqGNo1lYOlNVTU+pCYW1vr3axmtcKzz6r311/fLicZWqBoOia5ufUlYQxC8PuJuSya3BmATUW1PLCqiN0tVCjeVWZh0crjXPrxQVLjDLx3Tq/GZq6maP9J5BFCmVXr6ryHEbs6JX74IdTVMTxLaZ3bfSll31pC4/vvw+HDKqprxgzfxx5DaJOXpmPioby5y3/yYUEVz24q5d8/lTKucyKjcxLISjRya142dofk0o8PUlxr59z+afxpUi5p8S1oWjab0oZ8KRKoCT1JSSqK8NChlqPuxo5V2xw8COvXM3SmKk667UgFJ/RtRdMsLW3ZV1ZbC489pt7fcEO7bWEQ9WclhLhXCHFICLHe+eexgpoQYpYQYrsQYpcQ4u5wj1MTYyQkKHu6h7Icd03IYflF/fj9CTmYrQ7e2VXJR3tV6RijQfCv07qx/KJ+PDK1a8vCBNRDJDu7XZo2YpZu3ZSTvqVyLEIo5zzAZ5/R1VzKl3dOY/6Jfbzv12JREV4tCZSXXoKjR2HYMDjnnMDHH+VEvUBx8g8pZZ7z76OmK4UQRuAJYDYwHLhECDE83IPUxBhesvc7J5tYMKoTn8zrw8bLBvDJvIYHyoldk+ic7INyb7Xq+l3RhtGoTE6VlS0n2boEytKliCNH6JdmwtBaLkqFlxItJSXw73+r93fd1W61E4gdgdIaE4FdUso9UkoL8BpwboTHpIl2MjJCfwxdvyv6yMyEwYOVEPCUQd+vn4rAMpvhlVdYte0wC9/9yXsJlsLCln1ljz+uJi5Tp8KkSUE5hWglVgTKzUKIjUKI54UQnoyfPYADbp8POpc1QwixQAixWgix+vjx46EYqyZWSEz0rRptILjCR33t2qcJL9nZMHKk+p2qqppHZ91yi3p9800KthaQv2I/e4taiBCrq1P78OQr27NHJTIaDPCb3wT3HKKQqBAoQohlQohNHv7OBZ4CBgB5wBHg7205lpRysZRygpRyQm5uG3pGa9oHnTu3ve+9J2pq1ENLE72kpiqhkpqqHOruZYFGj1aZ8zU1DP/0bQC2tNQbpbyFJlwOB9x7r9rnBReoKhHtnKgQKFLKM6SUIz38LZFSFkop7VJKB/AMyrzVlEOAe6/Yns5lGo13fKlGGwgOh8p90EQ3CQnK/DV2rOouWl2tfB5mM9x0EwCDXn+ROANsOdyCQDl61LNp87nnYMUKNbG4/fbQnUMUERUCxRtCiG5uH+cBmzxstgoYJIToJ4SIBy4G3g/H+DQxTlKSSmjzVgXZX1zVZnX+SeyQkKAawY0bpyKxkpPV56lTia+pZpClzLOGUlurtNGm5eo3bYJHH1XvH3igw2irUS9QgL8JIX4SQmwETgPuABBCdBdCfAQgpbQBNwOfAluBN6SUmyM1YE0MIYRKcmytCZM/1NYqzacdR/O0WwwG9dsNHgxjxtRrKYlbV/PNd8sxGAz07duX/Px8tX1ZWfOwcLMZ7rxTRfldfrlyxncQoj6xUUp5eQvLDwNnuX3+CGgWUqzRtEqnTirZLVjU1irziSa2SU6G888nf+pUPv3oUcwox/2+fftYsGABWCzMHzmysbnLYoG771bFTwcNUoKlA6GnUBpNcrLKaPfUpyVQdLhw+8BgYOHevfXCxIXZbGbh73+vTGWuSL6qKrjuOlWxOCUFHn64w1WZ1gJFowmm2ctiUQ+RpKS270sTFew/cACEgdzz/khq3uyG5YWFDQJj/36YPx9+/FG1R3j5ZRg6NEIjjhxaoGg0oMxewXDM19SoUGRNu6F3794gHcTl9iGxz+iG5SYT/PnPcNZZqtjjtm2qAOVrr8HwjlmoQwsUjQZULkJcXNuFisOhy620MxYtWkRycjLWwj3Ed1b935OBRVar0kR271bXz9lnw6uvQq9e3nfYjol6p7xGExaEUL3HDx0KvCSLNne1S+bPnw/Awle+Q3Y6iT79BrLo0ouYn5ICR47AiSdCXp6uioAWKBpNA1lZyhYeKGYz9O4dvPFooob58+fTZewZXPuf1bz95QrGW0tUhnxaWqSHFlVok5dG4yIxsb6TY0BIGZ6Ck5qIMKJHOhP7ZmF3oHqm2GyhqbIQw2iBotG4061bYLW9tLmr3dMtI4k3bpjExH5Z6rfu0cN72foOiBYoGo076enKn+LvzNNsVj4YTbunvox9164qsz6YZXtiHC1QNBp3jEavjbdaRJu7OgTPfbeXvD8txWJzqGTYHj1UQqMG0AJFo2lOTo6qw+QrVVUqj6WDZUV3RHLTEqiss7GjsFItyM5Wk4mm/VQ6KFqgaDRNSU1VAsKXmaeUyn/Ss2fox6WJOKN7KC100yFnD5T4eDUB8VejbadogaLReKJ3byUoWpt5VlerB4ouVd8h6JOdTFqiiY2H3Jpqde0amq6fMYgWKBqNJ5KSlH28pW58oISN1aorC3cghBCM7J7RoKGAKgSZmuq5P30HQwsUjaYlunXzHsVTWakc+DpUuEPx8xN6MmdM98YLXd0eOzg6U16jaQmTSRX727lTZdG7N1Kqq1Ohxd26tfh1Tftk3lgP/rKMDHW92GzqtYOiNRSNxhvZ2aoVbGmpmoFKqbr02e2qoqyO7OqQFFXVcaTcLQHWYFDXSTA7f4aI93ZX8N7+AKtBtIIWKBqNN4SAPn1g9GglPEpKVO+UkSOV3VzT4ZBScvrDX/H4F7sar+jUKbhN2kLE85vLeGNfaARKx9XNWsBqtXLw4EFqA63npGnfZGSoWl87d0Z6JGEhMTGRnj17Eqcr6dYjhGBkjwx+OtQkYCMpSUX7WSwqnDgKsTkk20stXN4vNJq1FihNOHjwIGlpafTt2xfhbjPXaDoYUkqKi4s5ePAg/fr1i/RwoopRPTJ44fsCLDYH8SY3Q0/XrrBnT9QKlL3lFurskmEZoXn0a5NXE2pra8nOztbCRNPhEUKQnZ2ttXUPjOyRgcXuaMiYd5GRocykUZo5v6VE5csMz9QCJWxoYaLRKPS94JnRPZtkzLuIj1dCJZCK1WHgQKWVeINgQJoxJPvXAkWj0Wj8pHdWMn+/cAynDs5tvrJLl6hNcrw5L4u18/sTbwjNREH7UDQajcZPhBCcP76F+m1paQ0tEAzRN2dPjQvdmKLvbDUajSYGKKyo5Y3VB6i1NgkVdrVAiLKclKIaGwuWHWbtsdCZ47SG0hpbtqgSG8EiLU0lxHnhwIEDXHHFFRQWFiKEYMGCBdx2223BG4MfPProozzzzDNIKbnuuuu4/fbb69f17duXtLQ0jEYjJpOJ1atXc/z4cebNm0dZWRn3338/c+fOBeDcc8/lqaeeonv37p4PFCQee+wxnnrqKcaNG8fevXv54Ycfmm1z7733kpqayp133hnSsXhi8uTJHsfkoqysjFdeeYVf/vKXIR9LamoqVbqXR8Cs21/Gb9/ayMDOqYzr3anxyuxsKCyMzMBaYHNxHUv3V3PViMyQHUMLlNaorFQJS8GitLTVTUwmE3//+98ZN24clZWVjB8/nhkzZjC8FUEUbDZt2sQzzzzDypUriY+PZ9asWZx99tkMHDiwfpsvv/ySnJyc+s+vvvoqN9xwA+eddx5nnXUWc+fO5YMPPmDs2LEhFyYATz75JMuWLaNnlJaT9yZMQAmUJ5980i+BIqVESokhCs0r7ZmxvTMBWL+/rLlASUlRmordrl6jgK0lyq8zPCsBakJTHVlfgVFIt27dGDduHABpaWkMGzaMQ4cOAbBhwwamTJnC8OHDMRgMCCG45557QjKOrVu3cuKJJ5KcnIzJZGLq1Km88847Xr8TFxeH2Wymrq4Oo9GIzWbjn//8J7/97W+9fu8///kPo0ePZsyYMVx++eUAPPLII4wcOZKRI0fyz3/+E4CCggKGDRvGddddx4gRI5g5cyY1zoiaG264gT179jB79mz+8Y9/kOqWyb5o0SIGDx7MKaecwvbt2+uXv/zyy0ycOJG8vDyuv/567Ha712O0NFZP+/FEamqq1/3ffffd7N69m7y8PH7zm994HeOQIUO44oorGDlyJNdccw1PPPFE/XHuvfdeHn74YQDmzp3L+PHjGTFiBIsXL/b6O2h8p0t6It0zEll3oKz5SoMh6sxeW0rq6J5iIjMhdAIuqgWKEOJ1IcR651+BEGJ9C9sVCCF+cm63OszDDCkFBQWsW7eOE088kdraWi666CIefvhhtmzZwsKFC7nzzjv505/+FPD+zzrrLA4fPuxx3ciRI/n2228pLi7GbDbz0UcfceDAgfr1QghmzpzJ+PHj6x9Ul156KUuWLGHGjBn8/ve/58knn+Tyyy8n2Uu/kM2bN3P//ffzxRdfsGHDBh599FHWrFnDCy+8wIoVK/jxxx955plnWLduHQA7d+7kpptuYvPmzWRmZvL2228D8PTTT9O9e3e+/PJL7rjjjvr9r1mzhtdee43169fz0UcfsWrVKkAJzNdff53vv/+e9evXYzQayc/P93oMT2P1tp+WaGn/Dz74IAMGDGD9+vU89NBDrY7xl7/8JZs3b+bWW2/ljTfeqN//G2+8wUUXXQTA888/z5o1a1i9ejWPPfYYxcXFXsem8Z2xvTuxbn8LVoesrKjqN7+1xMKwrISQHiOqTV5Syotc74UQfwe8NKfgNCllUehHFT6qqqo4//zz+ec//0l6ejr/+9//GDduHBMnTgRg9OjRfPLJJyxZsoQPP/yQiooKrrnmGmbOnElBQQEvvvgiV1xxBYsWLaK8vJy33nqr2TE++uijFo8/bNgw7rrrLmbOnElKSgp5eXkY3dT37777jh49enDs2DFmzJjB0KFDmTJlCh9++CEApaWlPPjgg7z77rtcd911lJaW8utf/5pJkyY1Os4XX3zBhRdeWG86y8rK4r///S/z5s0jJSUFgPPOO49vv/2WOXPm0K9fP/Ly8gAYP348BQUFXv+P3377LfPmzasXanPmzAHg888/Z82aNZxwwgkA1NTU0LlzZ6ZMmdLiMTyN9ZVXXvG4H2/4eg7extinTx9OOukkAMaOHcuxY8c4fPgwx48fp1OnTvTq1QtQfqV3330XUP65nTt3kp2d7XV8Gt/I65XJhz8dobiqjuzUJg/rlJSoqUDskJLMBAPjOoe2mGlUCxQXQmVX/Rw4PdJjCRdWq5Xzzz+f+fPnc9555wHKpzFq1Kj6bdauXcu4ceOYO3cuc+fOpbS0lDvvvJOZM2fWb9O/f3+ee+45LrjggoDGcc0113DNNdcA8Pvf/76Rb6KHs7FU586dmTdvHitXrmTKlCn16//85z+zcOFCXn31VU455RQuuOACzjvvPD799NOAxuIiIaHhxjUajY3MUf4gpeTKK6/kgQceaLS8oKDAr2O0tB9v+Lp/b2N0CVsXF154IW+99RZHjx6t106++uorli1bxvLly0lOTmbatGk68z2IXDihJ/PG9WguTECFDnfpAkePQnp6+AfnhkEI3vxZr9AfJ+RHCA6nAoVSypYq8klgqRBijRBiQRjHFRKklFxzzTUMGzaMX/3qV/XLs7Oz2bhxIwA7duzgnXfe4eKLL65ff//993PTTTcFdSzHjh0DYP/+/bzzzjtceumlAFRXV1PpjH6rrq5m6dKljBw5sv57O3fu5ODBg0ybNg2z2Vzv7/H04Dz99NN58803600xJSUlnHrqqbz33nuYzWaqq6t59913OfXUUwM6hylTpvDee+9RU1NDZWUlH3zwAQDTp0/nrbfeqj/HkpIS9u3b53VfnsYayH5aIi0trf7/6u8YL7roIl577TXeeustLrzwQgDKy8vp1KkTycnJbNu2jR9//DGgcWk8k5kcT44nYeKiU6eoMnuFmohrKEKIZUBXD6sWSimXON9fArzqZTenSCkPCSE6A58JIbZJKb9p4XgLgAUAvXv3bn2AaWk+RWb5TFpaq5t8//33/Pe//2XUqFH1ZpG//OUvXHLJJbz//vuMHDmSnJwcXn31VbKzs5FScvfddzN79ux6Z76vnHXWWTz77LMtRmCdf/75FBcXExcXxxNPPEFmZiYAhYWFzJs3DwCbzcall17KrFmz6r+3cOFCFi1aBMAll1zC3LlzefDBB7nvvvuaHWPEiBEsXLiQqVOnYjQaGTt2LC+++CJXXXVVvXnv2muvZezYsa2atzwxbtw4LrroIsaMGUPnzp3rzUfDhw/n/vvvZ+bMmTgcjvpz7NrV0+Xofaye9tOnTx+/x5qdnc3JJ5/MyJEjmT17Ng899JDPYxwxYgSVlZX06NGDbs7GX7NmzeLpp59m2LBhDBkypN5Epgke7647yN7j1fxq5pDmK5OTISFBtYqOYMXmB1YVsamolvzZoY1+FDJKi5i5EEKYgEPAeCnlQR+2vxeoklI+3Nq2EyZMkKtXN/bhb926lWHDhgU42sjw2GOP8dJLL3HCCSeQl5fHDTfcUO9DueWWW1i4cCGfffYZ1157Lb/73e8iPVxNjBGL90Q4+b8lm3hrzUE23nsmRk8lTQ4fhkOHVI2vCHHu+/tJMhl47SynQCkrgxEjlJ8nAIQQa6SUE5ouj7iG4gNnANtaEiZCiBTAIKWsdL6fCTSfBrdjbr31Vm699VaP67Kzs3n66afDPCKNpuOQ1zuTl5bvY+exSoZ29eAr6dQJ3KIjw02d3cHWEgu/CGFCo4tY8KFcTBNzlxCiuxDCFZ7UBfhOCLEBWAl8KKX8JMxj1Gg0HZSxvVRS47r9ZZ43SEpS3T4toUkmbI1tJRYsDklebmhDhiEGNBQp5VUelh0GznK+3wOMCfOwop6+ffty7733RnoYGk27p092Mp2S41i3v5RLJrbgl+3aFQoKItJ4a8NxFdU3Jie0IcMQGxqKRqPRRC1CCCb2y6LO5mh5I2cwSyTokWrivIFpdEsJvf4Q9RqKRqPRRDtPXzbeezOyhARITVV9UhJCb3pyZ3rvVKb3Tm19wyCgNRSNRqNpIz51tuzSJey1versDkpqPdeVCwVaoGg0Gk0bkVJy2bMrePjT7S1v5MqWD2OqxurCWsa9soflR8IjyLRA0Wg0mjYihKDOZufbXV7KCcbHK19KGPvNr3c65IeHuCikCy1QNBqNJgic2C+bTYfKqa7zUmqlSxcIYy21Dcdr6Z8eR0YIS9a7owWKRqPRBIGJ/bKwOyRrWypnD6r0kqvxVhjYUFTLmNzQhwu70AJFo9FogsC4Pp0wGgQr9pS0vJHRqHJSqqtDPp6j1TYKzXZGhyH/xIUWKFHKtGnT2LZtGwDFxcWNKvlqNJroIzXBxBWT+jCoSyshutnZYdFQkkyCv57SmdN6tdzcLtjoPJQoZdeuXQwePBiAjRs3NuqDotFoopP/O2dE6xslJ6uijCHOSclIMHLR4PAWpNQaijeECM1fK+zbt48ePXpgMKifZ+PGjYwePZr33nuP6667josuuoilS5dSUFBQX15lz549XHPNNQE30tJoNMGhzGyhtLqVul3duoU8J2XpvioOVFpDeoymaIEShWzYsIHRo0fXf16zZg2jR49m7ty5PPPMMzz99NO8/vrrjb7j6syo0WgiR2WtlfH3L+PlH1tpsOYqZe/wUq6lDdTYHNz05RFe3uata3rw0QLFG1KG5q8V1q9fX9+mdefOnSxZsqSRySsUnRk1Gk3bSUuMY1DnVFYWeHHMg+ox36VLyJzza4/VYnXASV2TQrL/ltACJQrZsGEDDoeDMWPGcN999zF8+HBeeuklpJTcddddAXVm1Gg04eGk/tms2VeK1d6K9tG5s+rkGAJWHK3BIGBCl/BFeIF2ykclGzduZO3ataQ1aRf82GOPsWzZMsrLy9m1a1ejlrvFxcUsXLiQdevW8cADD+jOjBpNhDipfzYv/lDAuv1lTOyX1fKGycmQlaW0lAA7J7bEj0fMjMxOIC0+PAmNLrRAiTIqKysRQjQTJtC8M6N7f3XdmVGjiQ4mD8zGZBB8veOYd4EC0L07/PRTUAVKnd3B+qI6rhoW/pbDWqBEGWlpaezYsSPSw9BoNAGSnhjHY5eMZVQPHx7oqamqvpfZrDSWIJBgNPDDz/viCF8Nynq0QIlhdFdGjSY6OWtUN9837tEDtmwJmkAByEmKzKNdO+U1Go0myFjtDt5Ze5CVe1uJ9gJV3ys1NWhFIx9aXcSS3ZVB2Ze/aIGi0Wg0QcYoBH/5aCv5K1rJRwGV7NyrF1RVtfm4tTYHz2wu46ciL8Kprk6lLxiC//jXAkWj0WiCjMEgmDIol292HMfuizMjI0OFEVe2TbNYf7wWi11yUjcP+ScOB5SWqjpiw4dDUvBzVLRA0Wg0mhAwdUgupWYrmw75mK3eu7d66LchN2XF0RoEcEKXJsKipgbKypS/ZtQoZWYLAVqgaDQaTQg4ZWAOQsA3O4779oX4eOjfHyoqAj7m94fNDM9OaGioJaXSSkAJkp49VQn9EKGjvDQajSYEZKcmMLpHBruO++EbycpS5e0rK/3WImwOCUJwhqtcfW2tCkfu0UPlu4RQkLjQAkWj0WhCxCvXnURKgh+PWSGgb1/YvFk56VNb6a3ihskgeOOsnki7XZm3EhJg5Ei/9tFWtMlLo9FoQoRLmEgfisLWEx+vnOZGo1+RX3V2B1RVIcrLlVYSZmECWqBoNBpNSPnLR1v5Zf5a/76UkADDhqmqxOXlrVYpr62sZuIre3h+nxXGjFECJQwmrqZogRKl1NTUMHXqVOxBaBVqsViYMmUKNpstCCPTaDT+YBCCpVsKKWmt6VZT4uOVUMnJUUKlrAwsFhUFZrEoH0lpKZSWsrzUQblV0n94v5CEA/uKFihtJD8/n759+2IwGOjbty/5+flB2e/zzz/PeeedhzEIs4z4+HimT5/erCmXRqMJPWeP7obdIflk01H/vxwXB/36wdix0KeP0lRck8y4OBgwAMaOZVllPCnxRiYNyA7u4P0kKgSKEOJCIcRmIYRDCDGhybrfCSF2CSG2CyHObOH7/YQQK5zbvS6EiA/HuPPz81mwYAH79u1DSsm+fftYsGBBUIRKfn4+5557LgCnnXYan332GQB/+MMfuOWWWzx+Z9OmTUyePLn+89q1a5k+fToAc+fODZqw02g0vjOiezr9clL438bDge8kLk415Bo1CkaPVv4Rp/Yi4+JYtrWQKYNzSTCF38zlTlQIFGATcB7wjftCIcRw4GJgBDALeFII4ek/9lfgH1LKgUApcE1oh6tYuHAh5iZ9oc1mMwsXLmzTfi0WC3v27KFv374A/OlPf2LRokXk5+ezbt06/vnPf3r83vDhw9mzZ0+9mexXv/oVDz30EAAjR45k1apVbRqXRqPxHyEEZ4/uxo97ijleWRf0/f90qJzCijrOGNYl6Pv2l6gIG5ZSbgX1j2/CucBrUso6YK8QYhcwEVju2kCoL50OXOpc9BJwL/BUaEcN+/fv92u5rxQVFZGZmVn/ecqUKUgpeeSRR/jqq69aNIMZDAZGjBjB5s2b2blzJ3369Knv7Gg0GomPj6eystJjrxWNRhM6zs3rjtUukQS/pnyX9ER+O2sIpw3tHPR9+0tUCBQv9AB+dPt80LnMnWygTEpp87JNSOjduzf79jUv/ta7d+827TcpKam+pzzATz/9xJEjR8jOzm5VGJx00kl8//33PPnkk3zyySeN1tXV1ZGYGN6WoBqNBgZ2TuPu2UNDsu8u6Yn8ctrAkOzbX8Jm8hJCLBNCbPLwd264xuAcxwIhxGohxOrjx30sidACixYtIrlJD4Pk5GQWLVrUpv126tQJu91ObW0tR44cYf78+SxZsoTU1NRGQmL69OkcOnSo0XdPOukk/vCHPzBv3jx69GiQq8XFxeTk5BAXF9emsWk0msCw2h18s+M4xyqDU6YeYPPhcj7ceKT1/vVhImwCRUp5hpRypIe/JV6+dgjo5fa5p3OZO8VAphDC5GUb93EsllJOkFJOyM3NDeRU6pk/fz6LFy+mT58+CCHo06cPixcvZv78+W3aL8DMmTNZunQp5513Hn//+98ZNmwYf/zjH/nTn/4EgMPhYNeuXWRlNW4xOnToUBISErjrrrsaLf/yyy/52c9+1uZxaTSawDhSVsuVL6zkpR8KgrbPxd/s4e53NmKzR6A9oweixSnfEu8DFwshEoQQ/YBBwEr3DaRKQf0SuMC56ErAm5AKKvPnz6egoACHw0FBQUFQhAnATTfdxJtvvsny5cuZMWMGoHwpy5cr99GWLVs4//zzSWoSc/7oo4/ywAMPkNKkR/Urr7zC9ddfH5SxaTQa/+mdncyMYV3IX7GfGkvb88tKqi18/NNRzh/Xk6T4yEZ3uYgKgSKEmCeEOAhMAj4UQnwKIKXcDLwBbAE+AW6SUtqd3/lICNHduYu7gF85nfbZwHPhPodgM27cOE477bQWExtHjhzJI488Uv959+7dDB06lJqaGq688spG21osFubOncvgwYNDOmaNRuOda07pR5nZyjvrDrZ5X2+vOYjF7uDSE9vmsw0mwq8aM+2MCRMmyNWrVzdatnXrVoYNGxahEWk00Ye+J4KHlJJzHv+OGoudz+6YisHQLLLV5/1M//vXZKXE89aNk1v/QpARQqyRUk5oujwqNBSNRqPpCAghuPaU/tRY7Bwurwl4P8cq60AQVdoJRH/YsEaj0bQrzh7djZ+N7kacMfD5fJf0RD7/1VR86S4cTrSGotFoNGHEZDQQZzRQY7FzsNTc+heasK+4mopaK0IIjAGazEKFFigajUYTZqSUXPbcCm58eS02P3JI7A7JLa+u45LFP/rXYyVMaIGi0Wg0YUYIwdUn9+OnQ+U88+1en7/32qr9bDxYzoIp/T2Vqoo4WqBoNBpNBDhrVFfOHNGFfyzbwW4f+s6XVlt46NPtnNgvizljure6fSTQAkXTjHvvvZeHH3642fL33nuPLVu2+L2/goICXnnllfrPL774IjfffHObxhhump77Pffcw7JlyyI4Ik2sI4Tgz+eOJNFk4O63N+Lw4mGXUvKXj7ZSWWvjvnNHRqV2AlqgxCyR6L7oTaB4G09TgRKteOuO2fTc77vvPs4444xwDEt32mzHdE5P5I9nD0dKsHkVKFBqtnLtqf0Y0jV6q4VrgRKF/PnPf2bIkCGccsopXHLJJfXawrRp07j99tuZMGECjz76KJ9//jljx45l1KhRXH311dTVqV4Lffv2paioCIDVq1czbdo0QGkeV199NdOmTaN///489thj9cdctGgRgwcP5pRTTmH79u3NxvTDDz/w/vvv85vf/Ia8vDx2797dbDxXXXUVb731Vv13UlNTAbj77rv59ttvycvL4x//+AcAhw8fZtasWQwaNIjf/va3Hv8PH330EUOHDmX8+PHceuutnH322QBUV1dz9dVXM3HiRMaOHcuSJarSzosvvsh5553ncb9Lly5l0qRJjBs3jgsvvJCqqqr6/9Vdd93FuHHjePPNN3nmmWc44YQTGDNmDOeffz5ms9njubvO9ZNPPuHCCy+sP85XX31VP86WjunOtGnTuO2228jLy2PkyJGsXLmy/re6/PLLOfnkk7n88sspKCjg9NNPZ/To0UyfPr2+RUJhYSHz5s1jzJgxjBkzhh9++AGAl19+mYkTJ5KXl8f111+P3W7Hbrdz1VVXMXLkSEaNGlX/Wzz22GMMHz6c0aNHc/HFF3v8LTSh44LxPXnuqhOINxkoqbaw61hl/bpys5XCiloMBsHjl47lrjNDU7E4WOg8lFa46N/Lmy07e3Q3Lp/UlxqLnateWNls/QXje3LhhF6UVFu48eU1jda9fv0kr8dbtWoVb7/9Nhs2bMBqtTJu3DjGjx9fv95isbB69Wpqa2sZNGgQn3/+OYMHD+aKK67gqaee4vbbb/e6/23btvHll19SWVnJkCFDuPHGG9m4cSOvvfYa69evx2azNTsmwOTJk5kzZw5nn302F1xwQf1y13gArrrqKo/HfPDBB3n44Yf53//+B6gH//r161m3bh0JCQkMGTKEW265hV69GuqA1tbWcv311/PNN9/Qr18/Lrnkkvp1ixYt4vTTT+f555+nrKyMiRMn1msLnvablJTE/fffz7Jly0hJSeGvf/0rjzzyCPfccw8A2dnZrF27FlBVma+77jpAdcd87rnnuOWWWzyeO8AZZ5zBggULqK6uJiUlhddff52LL76YoqIir8d0x2w2s379er755huuvvpqNm3aBKh6bd999x1JSUmcc845XHnllVx55ZU8//zz3Hrrrbz33nvceuutTJ06lXfffRe73U5VVRVbt27l9ddf5/vvvycuLo5f/vKX5OfnM2LECA4dOlS//7KysvrfZ+/evSQkJNQv04QPIQQZSaoK+H0fbObTzYX0y0mh1mqnqKqO/rmpvPvLySTGRUe9Lm9oDSXK+P777zn33HNJTEwkLS2Nc845p9H6iy66CIDt27fTr1+/+vpcV155Jd98802z/TXlZz/7GQkJCeTk5NC5c2cKCwv59ttvmTdvHsnJyaSnpzNnzhyfx+saj79Mnz6djIwMEhMTGT58eLO+Mtu2baN///7069cPoJFAWbp0KQ8++CB5eXlMmzaN2tra+hm7p/3++OOPbNmyhZNPPpm8vDxeeumlRsdzP4dNmzZx6qmnMmrUKPLz89m8ebPX8zCZTMyaNYsPPvgAm83Ghx9+yLnnntvqMd1xnduUKVOoqKiof6jPmTOnvvjn8uXLufRS1UPu8ssv57vvvgPgiy++4MYbbwRUE7WMjAw+//xz1qxZwwknnEBeXh6ff/45e/bsoX///uzZs4dbbrmFTz75hPT0dABGjx7N/PnzefnllzGZ9Bwzktw9exizR3Wle2Yiw7qnc+aIrvzhZ8Oi1mfSFH31tII3jSIp3uh1fVZKfKsaib80rSLsCZPJhMOhYtvdG3UBJCQk1L83Go1tts+7j8f9uA6HA4vF0uL32jIOKSVvv/02Q4YMabR8xYoVHvcrpWTGjBm8+uqrrZ7DVVddxXvvvceYMWN48cUX+eqrr1odz8UXX8zjjz9OVlYWEyZMIC0trdVjutP0YeH67Mtv7QkpJVdeeSUPPPBAs3UbNmzg008/5emnn+aNN97g+eef58MPP+Sbb77hgw8+YNGiRfz0009asESIrhmJPPLzvEgPI2C0hhJlnHzyyXzwwQfU1tZSVVVVbyZqypAhQygoKGDXrl0A/Pe//2Xq1KmA8gusWaNMbW+//Xarx5wyZQrvvfceNTU1VFZW8sEHH3jcLi0tjcrKSo/rmh73/fffx2q1+vQ9TwwZMoQ9e/ZQUFAAwOuvv16/7swzz+Rf//pXfWLXunXrvO7L1cXS9b+qrq5mx44dHretrKykW7duWK1W8vPz65d7O4epU6eydu1annnmmXofhD/HdJ3bd999R0ZGBhkZGc22mTx5Mq+99hoA+fn5nHrqqYDSyJ56SnW7ttvtlJeXM336dN566y2OHTsGQElJCfv27aOoqAiHw8H555/P/fffz9q1a3E4HBw4cIDTTjuNv/71r5SXl3v09Wg0vqAFSpRxwgknMGfOHEaPHs3s2bMZNWqUxwdMYmIiL7zwAhdeeCGjRo3CYDBwww03APB///d/3HbbbUyYMKHF/vPujBs3josuuogxY8Ywe/ZsTjjhBI/bXXzxxTz00EOMHTuW3bt3N1t/3XXX8fXXXzNmzBiWL19eP8MePXo0RqORMWPG1DuCWyMpKYknn3ySWbNmMX78eNLS0ur/D3/84x+xWq2MHj2aESNG8Mc//tHrvnJzc3nxxRe55JJLGD16NJMmTWLbtm0et/3zn//MiSeeyMknn8zQoQ0OUG/nbjQaOfvss/n444/rHfL+HDMxMZGxY8dyww038Nxznjsv/Otf/+KFF15g9OjR/Pe//+XRRx8FVP+bL7/8klGjRjF+/Hi2bNnC8OHDuf/++5k5cyajR49mxowZHDlyhEOHDjFt2jTy8vK47LLLeOCBB7Db7Vx22WWMGjWKsWPHcuutt5KZmen1/6nRtIQuXx+F5eurqqpITU3FbDYzZcoUFi9ezLhx4yI6pkjg+j9IKbnpppsYNGgQd9xxR6SHFVSmTZvGww8/zIQJzSqBRw3RcE9oogtdvj6GWLBgAXl5eYwbN47zzz+/QwoTgGeeeYa8vDxGjBhBeXm57jip0UQ5WkOJQg1Fo4km9D2haYrWUDQajUYTUrRA8UBH1to0Gnf0vaDxBy1QmpCYmEhxcbG+kTQdHiklxcXFJCYmRnoomhhBZy81oWfPnhw8eJDjx49HeigaTcRJTEykZ8+ekR6GJkbQAqUJcXFx9eU+NBqNRuM72uSl0Wg0mqCgBYpGo9FogoIWKBqNRqMJCh06sVEIcRzwXFO8dXKAoiAOJ5K0l3NpL+cB+lyilfZyLm09jz5SytymCzu0QGkLQojVnjJFY5H2ci7t5TxAn0u00l7OJVTnoU1eGo1GowkKWqBoNBqNJihogRI4iyM9gCDSXs6lvZwH6HOJVtrLuYTkPLQPRaPRaDRBQWsoGo1GowkKWqD4iRBilhBiuxBilxDi7kiPJ1CEEL2EEF8KIbYIITYLIW6L9JjaihDCKIRYJ4T4X6TH0haEEJlCiLeEENuEEFuFEJMiPaZAEELc4by2NgkhXhVCxEyVSSHE80KIY0KITW7LsoQQnwkhdjpfO0VyjL7Swrk85Ly+Ngoh3hVCZAbjWFqg+IEQwgg8AcwGhgOXCCGGR3ZUAWMDfi2lHA6cBNwUw+fi4jZga6QHEQQeBT6RUg4FxhCD5ySE6AHcCkyQUo4EjMDFkR2VX7wIzGqy7G7gcynlIOBz5+dY4EWan8tnwEgp5WhgB/C7YBxICxT/mAjsklLukVJagNeAcyM8poCQUh6RUq51vq9EPbR6RHZUgSOE6An8DHg20mNpC0KIDGAK8ByAlNIipSyL6KACxwQkCSFMQDJwOMLj8Rkp5TdASZPF5wIvOd+/BMwN55gCxdO5SCmXSiltzo8/AkEpKa0Fin/0AA64fT5IDD+EXQgh+gJjgRURHkpb+CfwW8AR4XG0lX7AceAFp/nuWSFESqQH5S9SykPAw8B+4AhQLqVcGtlRtZkuUsojzvdHgS6RHEwQuRr4OBg70gKlgyOESAXeBm6XUlZEejyBIIQ4GzgmpVwT6bEEARMwDnhKSjkWqCZ2TCv1OP0L56IEZHcgRQhxWWRHFTykCo+N+RBZIcRClPk7Pxj70wLFPw4Bvdw+93Qui0mEEHEoYZIvpXwn0uNpAycDc4QQBSgz5OlCiJcjO6SAOQgclFK6tMW3UAIm1jgD2CulPC6ltALvAJMjPKa2UiiE6AbgfD0W4fG0CSHEVcDZwHwZpPwRLVD8YxUwSAjRTwgRj3Iyvh/hMQWEEEKg7PRbpZSPRHo8bUFK+TspZU8pZV/Ub/KFlDImZ8NSyqPAASHEEOei6cCWCA4pUPYDJwkhkp3X2nRiMLigCe8DVzrfXwksieBY2oQQYhbKRDxHSmkO1n61QPEDpxPrZuBT1M3xhpRyc2RHFTAnA5ejZvPrnX9nRXpQGgBuAfKFEBuBPOAvkR2O/zg1rLeAtcBPqGdNzGSZCyFeBZYDQ4QQB4UQ1wAPAjOEEDtRGtiDkRyjr7RwLo8DacBnznv/6aAcS2fKazQajSYYaA1Fo9FoNEFBCxSNRqPRBAUtUDQajUYTFLRA0Wg0Gk1Q0AJFo9FoNEFBCxSNRqPRBAUtUDQajUYTFLRA0Wg0Gk1Q0AJFo4kShBA3CyH2CSESnJ+HCSGOCiEujPTYNBpf0JnyGk2U4Owbsgn4N/Ae8A1wj5TyhUiOS6PxFS1QNJooQghxDvACUA48KqV8LMJD0mh8Rpu8NJroYhOQguoMqoWJJqbQAkWjiRKcfdg/R5m8pgohhkZ4SBqNX2iBotFEAUKIXOAz4EUp5e3Ay6gWuhpNzKB9KBpNhBFCZABfohqD3elc1gPYCZwjpfw8kuPTaHxFCxSNRqPRBAVt8tJoNBpNUNACRaPRaDRBQQsUjUaj0QQFLVA0Go1GExS0QNFoNBpNUNACRaPRaDRBQQsUjUaj0QQFLVA0Go1GExS0QNFoNBpNUPh/jC3R+MU/smwAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x432 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "fig, ax = plt.subplots(\n",
        "    nrows=1, ncols=1, figsize=(6, 6))\n",
        "ax.fill_between(X_tst.numpy().ravel().flat, mu21-2*Sigma2, mu21+2*Sigma2, color='red', \n",
        "                  alpha=0.15, label='$2 \\sigma_{2|1}:$ 95% confidence interval')\n",
        "ax.plot(X_tst, mu21, 'r-', lw=2, label='$\\mu_{2|1}$')\n",
        "ax.plot(X_tr, y_tr, 'ko', linewidth=2, label='$(x, y)$')\n",
        "ax.plot(X_tst, y_process(X_tst), linestyle='--', label='ground truth generative process')\n",
        "ax.set_xlabel('$x$', fontsize=13)\n",
        "ax.set_ylabel('$y$', fontsize=13)\n",
        "ax.legend()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "04bf6d46",
      "metadata": {
        "id": "04bf6d46"
      },
      "source": [
        "### Note on optimization and the Gaussian process:\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "7cfd6b90",
      "metadata": {
        "id": "7cfd6b90"
      },
      "source": [
        "Fitting a Gaussian Process is not as easy as the simple examples like normal distribution fitting as provided in the previous tutorial. This means that the results could be different if you change the initial values of the parameters, for instance. Give it a try!\n",
        "\n",
        "Therefore, we may need to do take extra considerations into account; for example, rerunning the optimization multiple times and use the results of the one with the lowest loss function value..."
      ]
    }
  ],
  "metadata": {
    "colab": {
      "provenance": [],
      "include_colab_link": true
    },
    "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.9.12"
    },
    "toc": {
      "base_numbering": 1,
      "nav_menu": {},
      "number_sections": true,
      "sideBar": true,
      "skip_h1_title": false,
      "title_cell": "Table of Contents",
      "title_sidebar": "Contents",
      "toc_cell": false,
      "toc_position": {},
      "toc_section_display": true,
      "toc_window_display": false
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}