{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lab 1 Extra: Uncertainty Propagation\n",
    "### Gaussian Process Summer School 2020\n",
    "\n",
    "This lab is an extension on the work introduced in Lab 1 of the summer school. It is more advanced, and you should make sure you've completed Lab 1 before attempting. It is designed to demonstrate the advantage of using models when we have only a small number of observations of a latent function $f$.\n",
    "\n",
    "Extra labs are for you to explore in your own time, giving details of other uses of Gaussian processes not covered in the summer school. Answers for extra labs will be made available _after_ the summer school."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Support for maths\n",
    "import numpy as np\n",
    "# Plotting tools\n",
    "from matplotlib import pyplot as plt\n",
    "# we use the following for plotting figures in jupyter\n",
    "%matplotlib inline\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "# GPy: Gaussian processes library\n",
    "import GPy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Branin Function\n",
    "\n",
    "For this lab we will focus on the [Branin-Hoo function](https://www.sfu.ca/~ssurjano/branin.html) as our latent function. The Branin-Hoo function in two-dimensions is typically used as an _artificial landscape_, i.e. a test function for optimisation. We will the most common version of the function\n",
    "\n",
    "$$\n",
    "    f(\\mathbf{x}) = a(x_2 - bx_1^2 + cx_1 - r)^2 + s(1-t)\\cos(x_1) + s,\n",
    "$$\n",
    "\n",
    "with paremeters defined:\n",
    "$$\n",
    "    a = 1, \\quad b = \\frac{5.1}{4\\pi^2}, \\quad c=\\frac{5}{\\pi}, \\quad r = 6, \\quad s = 10, \\quad t=\\frac{1}{8\\pi}\n",
    "$$\n",
    "\n",
    "We will define the function over the space $\\Omega = [-5, 10] \\times [0, 15]$.\n",
    "\n",
    "We can create some samples of the function, $N=20$, sampled uniformly over the space:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkUAAAHyCAYAAADlWbqRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdeZhcZZn38e+dpDsJpBM6IYGwKOKwCAwgE3xdUFmEQQbBDQV5URbldQcERQaUQRkUQRTF0YkguDA4gqIsLqATBhkBJ0BkEcQFlBAgEELorN1J3+8fVRUqlaruquqzPM85v8915Up3Lec8tZ5f3/dzzjF3R0RERKTsxuU9ABEREZEQKBSJiIiIoFAkIiIiAigUiYiIiAAKRSIiIiKAQpGIiIgIoFAkUghmdrSZ3ZT3ONphZpPN7HozW2ZmV2e87gfMbN8s1yki8VAoklIxs4lmdpmZ/dXMBszsHjN74yj3edTMVlVv/5yZ/cbM3m9mLT8/ZnaLma02s+XVjf+tZvb3yT+iCne/0t0P6vb+ZuZm9ncNl/2LmX1v7KPbyNuBLYAZ7n5ECssHwMyuMLNz6y9z913d/Za01ikicVMokrKZADwGvB6YBnwK+IGZbTfK/d7k7n3Ai4HPA6cDl41ynw+7+xRgBnAL8N1WNzSzCW2MvSheDDzs7mvzHoiISD2FIikVd1/h7v/i7o+6+7C73wA8AvxDm/df5u7XAe8E3mNmu7Vxn7XA94FdapdVqzDXmNn3zOx54Fgze4WZ3V6tRj1hZpeYWW/dfbxaofqjmS01s6+ZmVWvO9bMbmvntt0ys1eb2f9WK1//a2avrrtuKzO7zsyeNbM/mdn7WizjHODTwDurVbQTGitSZrZddfwTqr/fYmafNbP/qVbrbjKzzetuv0+1evecmT1WfS5OBI4GPlFdz/XV2z5qZm+o/jzRzL5sZouq/75sZhOr1+1rZgvN7FQzW1x9PY4by/MnIuFTKJJSM7MtgB2BBzq5n7v/FlgIvLaNdfRS2UDf0XDV4cA1wGbAlcA64BRgc+BVwAHABxvucyiwN7AH8A7gH0dYdSe3He0xTAduBL5CpfJ1EXCjmc2o3uQqKs/HVlTaY+eZ2QGNy3H3s4HzgP909ynuPlq1reZdwHHALKAXOK06rhcBPwO+CswE9gQWuPtcKs/pF6rreVOTZZ4JvLJ6nz2AVwBn1V2/JZVq4tbACcDXzKy/zfGKSIQUiqS0zKyHyobz2+7+UBeLWARMH+H6r5jZc8By4MPAOQ3X3+7uP65WrFa5+13ufoe7r3X3R4F/p9Lmq/d5d3/O3f8GzKOyQW+lk9sC3F2ttjxXHfcn6677J+CP7v7d6viuAh4C3mRm2wL7AKe7+2p3XwBcChwzyvo6cbm7P+zuq4Af1D2Wo4FfuvtV7j7k7kuq62/H0cBn3H2xuz9N5fWpH/NQ9fohd/8plddxp2QejoiESKFISqk6Sfq7wCCVwFK7/GfVdstyMzt6lMVsDTw7wvUfdffNgElUqjbXmNnuddc/1jCmHc3sBjN7stpSO49K1ajek3U/rwSmjLD+pret7oFVe4z1la693H2z2j8qc6dqtgL+2rD8v1J5DrYCnnX3gSbXJaXV494W+HOXy2x8TH+tXlazpGHe02jPt4hETqFISqc6t+YyKntAvc3dh2rXufsbq+2WKe5+5QjL2JvKRv+2VrepW+awu/8a+BNQv4eYN9z061SqLzu4+1Tgn4ExzQNqMZ5d6x7jr9u82yIqE6TrvQh4vHrddDPra3JdO1YAm9T9vmWb94NKsHxpi+san99GjY/pRdXLRKSkFIqkjL4OvIzKHmWrOrmjmU01s0OpTJz+nrvf1+b9XkVlovVIc5f6gOeB5Wa2M/CBTsaWsp8CO5rZu8xsgpm9k8rjucHdHwN+A3zOzCZVq2EnUGlNtmMB8Doze5GZTQPO6GBcVwJvMLN3VMc1w8xqrbWngO1HuO9VwFlmNrM6cfvTQBqHIBCRSCgUSamY2YuB/0dlTsqTHbTKrjezASqViTOpTDQebW+kS2rLp9KqO8vdfzbC7U+jMqF4APgm8J+jP6JsuPsSKi3AU4ElwCeAQ939mepNjgK2o1JpuRY4291vbnPZN1N5rPcCdwE3dDCuvwGHVMf1LJWAtUf16suAXapzpH7c5O7nAvOr670PuLt6mYiUlLmPVmEWERERKT5VikRERERQKBIREREBFIpEREREAIUiEREREUChSERERASonDE8Gr3jJ/vk8VMrv6wr8Qm2xyf0sk0Y39XdfHznxxPs5j7DXdznhfV1fddE7t90mRn9CZLG2Fsal87eqzY++eWOHzec+DIbTRif/jrytnZdum/kdcPJLd/XJX7sUxge+zIHH1v4jLvPTGA0iXvDfpN8ybPJvo8X3Dv0C3c/ONGFpiSqUDR5/FReveVRAAwvGensCsU2bsZIp9tqn/f3jX6jJob6J3d+n2k9Xa1rzbTutvCDfQl8cXX39IxoKMOTRKydms0Get2Udakst2fKYOLLnNbX0bE6uzZr0+WZrCdLi1ek++ZdNtD598pIhpb3Jrq88cuT+2vjkZNOazxlTjCWPDvMvJ9tkegy+7de2Hi6omCpfRahvANhz9LONyw9y4ZGv1GCegfGXmnoHRj9Np3qyXBbOeH5bD7eSW4s6iW9UYPkN7ytpB0gsrR4xRQFopTe4xKeaENRUtWSMrOlKWz1EzZxWTpViDwpGLVvaHlv4hu4ZQOTMwlHtTARa0DKauwKRBKSaENR2eVdLeqGqkUvyDoYZRGOxi8fr6pRCzGFoyzHqkAkoYlqTpGEo2fpqq7mFnVj4rJ1Xc8tSkLvQDrzi2rBKKt5RhOeH5fJPKPxy8enMs9oaHlv4vOMlg1MzmyeEWzYVgtl3lFeYU2BSEKkUFRytnSg6wnX3ehZNtT1pOtu9A54IpOu0wpGUAlHWQYjSH8SdmzBCLKbhF2TZ0DKs2qVRoUulkCUVTtbuqdQFLHhJc+WZm7VWKpFCkYby6JqVNuwJB2O0ghGkH3VqF6zkJJkUAqldadAJKGLOhSNmzE9yrk1RZFlC60M8ghGEGfVqLYhjL2dNpJOgkwtQIUSfpopc7tMgSgeeqUk873Qup1wPZY90ZKYdF1ZTiKLaSnLCdg1WU3CTkNaE7CznISdhJAncqfxfCoQSVr0akUu70pZN8csyktMwSjrcKRgtLHYglGIyt4uUyCKj14xAcpRLUpS2sEI8glGaX+Jp7XbfprBSOGoc2k9bzEFIolT1HOKRrPfO17FceccwcxtZvD0wiVcfvbVzPvB7XkPq/T2P2AXTnjffsyaNZXFzzzP3O/eyi9vfSiTdSc16bqyrPQmXtdkNc/o0F134tT99mH2tD6eWDbABXf8museTu81SWueURoTsCG/uUYHz96DD+10IFtO2ownVz/H1/5wMz9/4neZj6MTaYXIUALRYTvuzGmv3oet+qayaOB5LvzNbRt8VhSI4lbYULTfO17FyZccz6RNJwKwxYs25+RLjgcoXDDKey+0TiZc73/ALpz68UOYNKnyBbflrGl84kOV8wS2G4zGetyiGIMRpBeODt11J879pwPZpLdyqIStN5vKefsfxLhVxo8fezCdlRLXBGzIftf9g2fvwVm7vZnJEyqPaavJ/Zy125sBggxGaVbUQgpE5x1wEJv0VD4r20ydxnkHHATAdQ8/pEBUANG/gq3CwHHnHLE+ENVM2nQix51zRBbDilIWLbQT3rff+kBUM2lSDyce87rU152WLFppkF477dT99lkfiGo26e3h1P32yaSdloa02mmQ3VyjD+104PpAVDN5Qi8f2unATNbfiTIEIoDTXr3P+kBUs0lPDx9/5WsViAqisK/izG1mdHS5ZGPWrKnNL9+8+eWtjHVuUVKTrl9YXqKLaymNSdizpzUvddUu1zyjjWUx12jLSZt1dHke0n4eQgpEAFv1Nf+eavUZkvgUNhQ9vXBJR5fHLpa90BYvfr755c80vzwmWQUjSDYYPbGs+cDrL89qEnbS0gxGkG4oeHL1cx1dnrW0w1BogQhg0UDz76lWnyGJT2FD0eVnX83qFWs2uGz1ijVcfvbVOY0oDmm30C775jxWr95wvsfq1UPM/e6tHS8rtGpRZZmJL7KlpILRF+fdxsrBDfcGXDk4xBfn3bbRbWMNRjGGo6/94WZWrd3ws7Jq7SBf+8PNia6nUzFWh5J6X134m9tYOdTeZ0XiVNiJ1rXJ1Nr7LCz/9avfA7yw99ni55l7ZXZ7nzVKctL1C8tMf/J1TRKTsG944A8AG+x99sV5t62/vFHaR8KO6bxpjZKcjF2bTB3K3mdZzKUKsTpUr7aX2cdf+dq2PisSH3NP/q/ltEzr3cJfveVRTa/Lu30UiiT2QhvLCWK7Pe3HWE4SO5Y90YDEQ9ELy01lsS1leYqQmjRPEZJGMIJ09kxrJpTThYxVVhPLQw9EkFyl9I+f+thd7j4nkYUl7OV79Pq8n22R6DL7t1444uM1s22B7wBbAsPAXHe/2Mz2BL4BTALWAh9099/W3W9v4A7gne5+TcMyNwGuBl4KrAOud/dPjjbWwrbPRNqVRhutstxUFttS0Y6EHeOeafVqbaZYD/6Y5djLFIikqbXAqe7+MuCVwIfMbBfgC8A57r4n8Onq7wCY2XjgfOAXIyz3QnffGXg58Boze+NoA9GrLBsZy7yibk/70e0RriGZo1wXJRhBsY6EneaeaVmFI4jryNhZhyEFInH3J9z97urPA8CDwNaAA7Vd/qYBi+ru9hHgh8DiFstc6e7zqj8PAncD24w2lsLOKSqrvA/kGLM05hdVlpt9Ky3tAz42M+H5cZpnNIr6sBFSey2PwJZGKFUgip+ZbUelsnMncDLwCzO7kEoR59XV22wNvAXYH9i7jWVuBrwJuHi02+rVlmDkXS1KUx4VI1A7rR1ZVozq5d1ey3P9CkSls7mZza/7d2KzG5nZFCrVn5Pd/XngA8Ap7r4tcApwWfWmXwZOd/dRv/jNbAJwFfAVd//LaLcvTKVo3IzpmmydIFs60PWE605O+xGatKpFlWVnXzGC7KtGqhh1rjGYpFFFCqV9p0AUttVu/HFt4s/nM6NNLDezHiqB6Ep3/1H14vcAJ1V/vhq4tPrzHOD7ZgawOXCIma119x83WfRc4I/u/uV2Bpp6KDKzbwGHAovdfbeG604DLgBmuvszaY+lLMraQhvrOdFq0g5GkF84yjIYQTp7p6UZjCC7vdNGMlqAGSk0hRJ+GqVVkUs6EJU1DOXJKunmMuBBd7+o7qpFwOuBW6i0yv4I4O4vqbvvFcANzQKRmZ1LZS7Se9sdSxaVoiuAS6jsbrdedRe8A4G/ZTAGiUTPsqEx7Z6flDSDUWX5qhqNRW1DWMSqUTtCDT6tKBBVZN3OjshrgGOA+8xsQfWyfwbeB1xcbYGtBpq23eqZ2QJ339PMtgHOBB4C7q5WlS5x90tHun/qocjdb61OnGr0JeATwE/SHoN0J8YWWlLVIihuMILsq0YxttMgjKpR7NQuq1Agas3dbwNafdn+wyj3Pbbh9z2r/y8cYZkt5VInNLPDgMfdPZ9Ds5ZAzPOrxjLhGsKfdF0vrwnYkO1xjdLebT8teU3CLoK0DnugQCRpyjwUVY8yeSaVAzG1c/sTazPWB4dHnnxYxnk0Iev2mEUhSev4RRuuI/9wlBUFo3KIpV0G6QaiPA6oKmOTR6XopcBLgN+Z2aNUDqZ0t5lt2ezG7j7X3ee4+5zecXH10Ysg7RPEthJStSiLYFRZTyaraSrrqlEa0g5GCkejS/N5ijEQSXwyD0Xufp+7z3L37dx9O2AhsJe7P5n1WETaVYZgBApGo1Ewai3N50aBSLKSeigys6uA24GdzGyhmZ2Q9jqlIoR5RWNpoYVULYJsg1HeVaMsxHZqkBpVjTaU9vOhQCRZymLvs+antX/h+u3SHoOMzVj2QiuatPdI23Bd+e6dBunvoRbj8YxqtIda+pWzmAKRwlAxFO4oVZpsXSyhVYuypnba2KTdToNyttSyqJalcQwiBSIZTeFCkWwo9hZaEmJto72wvvzbaVl86ccejMoQjrIKQzEdlFGBqFgUiqQtee2FFqqsg1FlnZmvcgOxByOFo+5l9bhiapeBAlERKRRJJvKccA3ptNHyCkZ5V43SluZGLItgBMUKR1k9DgUiCUEhQ5HmFW0ohBZaCIoSjCrrzWW1QDbttFiPgN0o5nCU5dhjCkQ6IGOxFTIUSTrybKElUS1KS57BKO9wlLYiBCN4IWCEHpDyGGdsgUiKTaFIMpP3hGtIb2+0vIJRZd25rVrBqAshhqO8xqRAJKFRKCqJIrTQkqoWFTUY5RWOsmqnpSGrCdjN5F09ynP9aT3vCkQyVgpF0hHthTayPINRZf35rTvWYAT5VY1q6gNKWiEli3W0I63nWoFIkpD6Ea3zMm7G9EJUR4qmZ+kqhvq7P7Fvz7Ihhqb1jHkcE5etY820dL6cszzqdfP1V/7P42jYPcvTPQr2hOfHpXL0a0j/CNidGC20jHQU7dBac/UUiCR0hQ1FsrHhJc9qz7yM5B2MKmPILxhBeuEo7VODAMGEo1ZCDj6txBSIFIbKS+0z6dhYW2hjnXAd+tyimt4BD6KdludcozQVuZ1WJGnO21IgkqSpUiSllmYbrSaUqhFkXzlSO63cYgpDoEAEMOgT+MvgzISX+teEl5eeQleK1CoqrpCPW9RM3hWjmjyqRrFXjFQ16o4CUXN5n65HRlboUCQbS2ryed4ttCSl3UarCSkYZf3FnPZu+2kGI1A7rROxtcsgm0CU98FWpT0KRRKtJKtFZQtGULyqURbBSOFoZGk+P7EHIomDQpHkJqRqEWQbjEIJR3lVjdKSdjACVY2aSTswKhBJVhSKSiiUFloSkp5blFUwgnJXjdIORqoaZSft5yHWQKR2WZwKH4o02VpCFlowyvJLPPZ5RlDuqlEWwTDmQCRxKnwokrAl0UKLuVoEYbXToHhVo7SVsWqUxeNN6xhECkQyEoUiGZMQWmhpyDoYgapGackiGEE5wlFWjzHGgzKqXVYMCkUlFdJ54UKsFoGCESgYdaOI4SjLxxRrIJJiKEUo0rwiiUmZ22lFCUZQjPlGWYchBSLJWylCkaQrlBZaUapFNaEFo6y+/GPfM61eLVTEFpCyHrMmVEsoFIpKrGgtNChmMAotHGWhCHumNYohHOUxRgUiCYlCkUgb8gxGoKpRGvIIRhBeOMqzmhVjINKE6mJTKJJEJNFCC7laBGEEo9DCURaKGIxgwzCSdSAJoa0XayCSYitNKNJk6+ZCaqHFIO9gBOFVjbJQpHlGraQZkvIMYM0oEEmoJuQ9AJE09CwbYmhaTyrLnrhsHWum5bthqQWjwT7LdRzwwsZisC/d9fQsh6Ep6S1/wvPjWDt1OL0VdKid8LJuyoYhPYTAM5I0w6cCkSQh/z+POjEh7A982YXUQktbCBUjCKullsWGI+3JsyFUjDoRWgVoJApEEoO4vgFEOpDW3KKaUIIRhNNSUzCSZhSIJBal+vRrXlFzoc0riqVaBOEFoxDCURZ752Sxy77CUTIUiGQ0Zratmc0zswfN7AEzO6l6+b+Y2eNmtqD675Dq5Qea2V1mdl/1//1HWPZHzOwP1eV+YbSxaE6RJMqWDuD9KU8u6UCac4tqQphjVK93wIOZa6R5RuWVdqhMKxApDOViLXCqu99tZn3AXWZ2c/W6L7n7hQ23fwZ4k7svMrPdgF8AWzcu1Mz2Aw4Hdnf3NWY2a7SBRPenUEgbXIlD2m00CKtiBGFVjdKWRTtNVaPOKBBJJ9z9CXe/u/rzAPAgTUJO3e3vcfdF1V8fACaZ2cQmN/0A8Hl3X1O93+LRxqJPugBqoSUhtGAEYYSjIgQj0FyjdikQyViY2XbAy4E7qxd92MzuNbNvmVl/k7u8DbinFnwa7Ai81szuNLP/NrO9R1u/2meSuNBaaJBNGw3Ca6XV5N1Sy2K3/bRbaaB22khiDUOgQFRv9XAPf1qzRdKL3dzM5tf9Ptfd5zbeyMymAD8ETnb3583s68BnAa/+/0Xg+Lrb7wqcDxzUYr0TgH7glcDewA/MbHt3b/mXYulC0bgZ04OrikhzPUtXMdQ/ObnlKRgB+R7bKO15RlkFI0DhqI4CUbvryr+lnZNn3H3OSDcwsx4qgehKd/8RgLs/VXf9N4Eb6n7fBrgWeLe7/7nFYhcCP6qGoN+a2TCwOfB0q3FEWQ8OrQoh0ijEVlpN3i21LPZMy4LaadnMt1IgKj4zM+Ay4EF3v6ju8tl1N3sLcH/18s2AG4Ez3P1/Rlj0j4H9q/fZEeilMkm7JX2qZb0kK2hJHMgRkp9blMWk65qJy9YFH47yW3e6y097l/2aMk/CzuJxKxCVxmuAY4D9G3a//0J1t/t7gf2AU6q3/zDwd8Cn6m4/C8DMLjWzWlXqW8D2ZnY/8H3gPSO1zqCE7TORrNpoNaG20yDflloRdtmvKVNLLasQWIRApDDUHne/DWj2JfTTFrc/Fzi3xXXvrft5EPi/nYwl2j9xxtJC00Ec4xLjnmiNQq4YQX5f3lkd6DErRa8aKRB1sh4FohgV+xMsHQuxhZaGLNtoNTEEozzDUZqyDkZFC0dZPiYFIslTsT65UlhpVIvyCkYxhKN81pvu8rMMRlCMcJT1Y4g9EOW9E4OMXdyfWJFIxRCM8vhyL1owgjjDUR5jLkIgkvjF9UltoHlF4Qu5hQb5VItqQg9GkM8XfRGDEcQRjvIaY5pHqVYgkk6E/QmVXIR6cMu0JlznHYxCD0d5VI2Ksst+M7XgEUpAyns8sZ+2Q4GoWLRLvkgAQt5tvybrU4UUaZf9VuqDSJa784cQyNQukxDl/8mQwkuyhVbEalGNqkbN1lfcdlqj+opN0qElzWV3Q4FIQhV9pcj7+7re6Oo8aFIv64M6tqKqUbP1xX/OtE6NFl5aVZZCCD0jUSCSkKX+6TGzb5nZ4uphtmuXXWBmD5nZvWZ2bfU8JhKQpMNiDNUiCKNiBPFUjbJdX7rLD6Vi1K7G6k8oVaCRKBBJ6LL4BF0BHNxw2c3Abu6+O/AwcEYG4xBpSyjBCMIPR3m009KU5wTsolMgkhikHorc/Vbg2YbLbnL3tdVf7wC2SXscUixFOPVHJ0IORpDtBiOLDaCCUbIUiCQWIdRajwd+NpYF6HhFcQj9mEX1QqoW1cRQNcpuXemvQ8Fo7NKuvKX9PtARqssn11BkZmcCa4ErR7jNiWY238zmD65dkd3gJPhJ6GlXi0IMRhB2OMpyI5JVMFI46k7az1sWgUjKJ7dQZGbvAQ4Fjnb3lu8+d5/r7nPcfU7vhE2zG6AI4QYjCLulVqRgBApGnVIgkljlEorM7GDgdOAwd1+ZxxgkH0m30LKYWxR6MAo1HCkYlU8WlTUFIklTFrvkXwXcDuxkZgvN7ATgEqAPuNnMFpjZN8a6Hs0rSkfoLTSpCDUcZdVOy+ocV2qntZbF86JAJGlL/eCN7n5Uk4svS3u9Uh49S1cx1D853XUEcmDH0dSCUWgHf8zqYI9ZnBoEwjzYY54UiKQooj+itcTHlg6MqbKXl1obLZZwpGCULgWj7KpmCkTZGfQJ/HXVjLyHkZsQdskXGbMsj1sU8hyjeiG21Io4z6is7TQFIimiQoUizStKh+YVbSyWYAThhaOiBSMoVzjK8rEWMRCF9FmUjRUqFEk80jiQY9ZHuY4pGEFY4SjLCdhZKnowyvLxKRBJHhSKRMYgtmAEYX0xFzUYFS0cZf2YFIgkLwpFUih5nBMt1mAUypd0EYMRFCMc5fEYihaIQvqsyegKF4o0rygdacwriulcaKOJMRhBOF/YRQ1GEGc4ymvMab5GeZzHLITPlnSmcKFIJI9qEcQbjCCML+8iByOIIxzlOca0A1HWQvhMSecUikQSFHswyvuLvEhHv24lxHCU95iKFIhC+BxJ9xSKJFdptdDyqhZB3MEI8v9SL+Iu+83UgkheYSTv9dcULRBJ3AoZijSvKB2xHa8o72BUhHCUl6Lust9KlgElhCBUo0AkodFpPiR3sZ72ox2xnDOtlbzPpZbFqUGyOi1IuxoDS7enEgkl+DRTpD3MFIaKRaFICi2Lk8WOOobIgxHkey61MgajeiGHm24oEEnICtk+g7G10KS12FpooShKOy2vjUCZWmlFpkAkoStsKJK4pHnMojznFjWKPRhBfhuDMuyZVmQKRBIDhaImNNm6eEILRrGHoyIHo8p6MllNaSgQSSyiCkU+Pt15BSJZKkIwymMDoWAUl6IEorwPVSHZiCoUdUrzitKR1ryitE/7EVK1qEZVo+4oGMWhSIFIyqHQoWgs1EIrphCDEcQfjvL4KzrLYKRw1JksnjMFIklDdKEo792rJV1FOklsN2IORpD9BiTLuSQKRu3J4nlSICoWM9vWzOaZ2YNm9oCZnVS9/AIze8jM7jWza81ss+rlM6q3X25ml4yw3COqyxs2szntjCW6UCQyVqFWi2qKUDXKkoJROBSIpEtrgVPd/WXAK4EPmdkuwM3Abu6+O/AwcEb19quBTwGnjbLc+4G3Are2O5DChyLNK5JmQg9GEHc4UjAql6xajApExeTuT7j73dWfB4AHga3d/SZ3X1u92R3ANtXbrHD326iEo5GW+6C7/6GTsRQ+FI2F5hW1luZBHMveQmsUazjKep5R1sFI4agiq+dBgagczGw74OXAnQ1XHQ/8LO31R3maj6H+yVH8pS9hC+EUIJ2oBaPYThmS5SlCsjglyIbrC/f0IFlQICqeweHxPLZys6QXu7mZza/7fa67z228kZlNAX4InOzuz9ddfiaVFtuVSQ+sUZShSCQpsQUjiDMcZR2MgMzCUS0YlCkcZVklUyAqhGfcfcSJzmbWQyUQXenuP6q7/D3AocAB7p76m6EU7bOxzCtSCy0faqGNLra2WpHnGVXWl+nqcqNAJEkzMwMuAx5094vqLj8YOB04zN1XZjGWaENRbH/dF1FRTg4beys2pnBU5HlGlfUVNxxl/dgUiErlNcAxwP5mtqD67xDgEqAPuLl62TdqdzCzR4GLgGPNbLjkO3YAACAASURBVGF1bzXM7NLa7vdm9hYzWwi8CrjRzH4x2kDUPpNg2dKBzPYejLGN1iimtlqR5xlV1lmcdloeIS+LQKQwFI7qnmTNPqQ/HeE+27W4/L11P18LXNvJWKKtFHVKu+bLaGKvGNXEUjkqcsWoss74q0YKRFI2pQlFY6F5ReVRlGAEcYSjogejynrjC0d5jVmBSPIWdSiKvd1RBGnPK9KE67ELPRyVIRhV1h1+OMpzjGUJRCF/FkVzikQ2UoT5Rc2EPOeo6HOMNlz/Cz+HMO8ohKBWhkCkMBSHqCtFndKu+XHKo1pUpDZao1ArR2WpGNXLtzKjQJSVED9v0lypQpFIJ4ocjCDML+oyBiN4IaCkHVSyWEcnFIgkNNGHoiK2OWJTlOMVNVOGYBTal3aWxzIKKRjVSyK8NAatUIJQTdEDUYifLRmd5hRJFLI8ZlGjos4xqhfifKOs5hnlPcdoNKGFmSSUIRBJnKKvFHVK84qkG0WvGNWE9tdt2StGRVTkQBTa50c6V4hQVPS/4qUi793zyxKMIKy/dLMMRgpH6Sry8xvSZ0a6V4hQJPkr8ryiemULRqF80Zd1AnaRFPlcZqF8TmTsFIo6pBaalCkYQThf+ApG8VIgkliUMhTpPGjxyruFVlPGYBTCl7+CUXyKGohC+UxIsgoTijSvSLJWtmAEYfxVrGAUjyIHIimmwoQiyV9W84pCqRZBeYNR3hsFBaPwKRBJjEobirRrviSlZ+mq0oajPGUdjBSO2pPlc6VAJEkrbSiSuIVULaopazDKc0OR9UZRwWhkWT4/CkSShkKFIs0rkryVMRhBvhsMBaMwKBBJEeg0HyIJqwWjsoX0nmVDuZ0mJKtTgtSEfmqQrBU1EJUxDK1dN47FK6bkPYzcFKpS1CnNK0pelgdxDLGFVq+Mc43ybKflUTEqe9Uo6+dAgUjSVupQJJKFsgUjyG+DksfB+8oajLJ+3ApEkoXUQ5GZfcvMFpvZ/XWXTTezm83sj9X/+5NaX9laFmUXerWopqxVozzkFYzKFI4UiKSosqgUXQEc3HDZJ4FfufsOwK+qv+dCLTTJUtnCUV7ttLzOkl70YJRH+FMgkiylHorc/VagcaLJ4cC3qz9/G3hzW8sar4mNMcj65LCxVIvqlTEcZS3PYFS0cJTXY1IgkqzlNadoC3d/AqD6/6wkF64WmsSiTOGoTMEIilM1yutxKBBJHoKfaG1mJ5rZfDObPzS4IpV1qIUWvxirRfXKEo7KGIxiDUd5jl2BSPKSVyh6ysxmA1T/X9zqhu4+193nuPucnt5NczsOikgWyhCOyhaMIL5wlOdYFYgkT3kdvPE64D3A56v//yTpFQz1T+bAV+7AcR99AzO33Iynn3yOy7/yS2752b1Jr0o6tO/he3Hs6Ycyc6t+nl60lCvOv4FbfnL3mJdrSwfGVPULSX0wyqIdvP8Bu3DC+/Zj1qypLF78PJd9cx7/9avfp7a+PA70mPUBHpuphY12D/x48Jyd+Mhh+7Dl9D6efHaAr153Gz+f/4fUx5eXGAJR1p8VyVbqocjMrgL2BTY3s4XA2VTC0A/M7ATgb8ARSa93/wN24eTTDmHS5F4Attiqn5PPPhwg8WA0bsb0zCcXh254ybNNW4v7Hr4XJ51/JJM2qb4u20znpPOPBEgkGBVR2kfI3v+AXTj144cwaVLlNdlyy2mc+vFDAFIPRkCm4SiEYATthaOD5+zEp991IJMnVp6frWZM5dPvOhAg0WCUdxCqiSUQ5fFZkexksffZUe4+29173H0bd7/M3Ze4+wHuvkP1/8QTxQnv2299IKqZNLmX4z76hubjLEiFIXTHnn7o+kBUM2mTXo49/dBElh/73KKR1FprSbfXTnjffuu/5GsmTerlhPftl+h6Wsm6hZF3K63eSG21jxy2z/pAVDN5Yg8fOWyf1NedtRgCEeT/WZH0BT/Rupl2/rKcNWtq08tnbrlZ0sORDszcqvlxOltdLs0lGY5afVZaXZ6GMgcjeCGg1AeVLac3/0Ot1eWdrisUsQQiCOOzIumKMhS1Y/Hi55te/vSTz6WyPu2F1p6nFy3t6PJuFLla1CiJ6lGrz0qry9NS9mBUr3fAeWpJ8/fxk8+29/5uDFqhhSGIKxBBOJ8VSU9hQ9Fl35zH6tWDG1y2etUgl3/llzmNqHyazbO64vwbWL2y4XVZOcgV59+Q1bAKq9uA1PSzsnqQy745L8nhtUXB6AX/dvWvWbVmw+dj1Zohvv6DX7cMPKGGn7wl9b4K6bMi6chr77MxG5rWM+IbvTbp7b0nvL7tvc+8v69UVYY81CZTp7H3Wb0i7YnWjU4mZ9c+K6HsUZP1nmmhTL5udNPtlcnUHzzitWwxo4+nlgzwb1f/ev3lMYutQlSTxGel6IfciJ25x/MXRd9m2/jLX3/S+t/bfbN38iYcayjSXmgbyrOtWOZQ1CjGo7xnvct+iMGoiGINREmobYt+ce+5d7n7nJyH09SmO8z2nb9yfKLLvPuQ84J9vI0K2z4TUdXvBWntuZYmtdKKp6yBKLbPXpkpFDVQdUGKLqYvaAWj4ihzIJLRmdm3zGyxmd1fd9meZnaHmS2onu7rFdXL+83sWjO718x+a2a7tVjmFWb2SPX+C8xsz9HGEXUoCvGUH9oLbUN5txNVLWotlnCkYBS/MgaiWD5fAbkCOLjhsi8A57j7nsCnq78D/DOwwN13B94NXDzCcj/u7ntW/y0YbRBRh6J2xTifQpKjYDSyGL688whGCkfJKOPzGPrnKUTufivQ+Fe0A7WDQE0DFlV/3gX4VfV+DwHbmdkWSYyjFKFIREYXejgq44lkY5f185d3lSj0z1CETgYuMLPHgAuBM6qX/w54K0C1pfZiYJsWy/jXapvtS2Y2cbQVRh+K0mihjXVekVpo4VG1qH0hf7ErGMWjjIFIWtq8Oieo9u/ENu/3AeAUd98WOAW4rHr554F+M1sAfAS4B1jb5P5nADsDewPTgdNHW2G0xynq1FD/ZL1pS67sxy7qVNonou1W1scxgnCPZRQqBaJ4rRsex7KBxD/zz3S5S/57gNpxeK4GLgVw9+eB4wDMzIBHqv824O5PVH9cY2aXA6eNtsLoK0USvrwnW8vYhPiFr4pRuBSIJEGLgNdXf94f+COAmW1mZrUz874XuLUalDZgZrOr/xvwZuD+xts0Kk2lqFM6unUxqVrUnRCrRqoYhadMgUhhKFlmdhWwL5VW20LgbOB9wMVmNgFYDdTabi8DvmNm64DfAyfULeenwHvdfRFwpZnNBAxYALx/tHEUIhSNdsqPPIybMV0VkkApGHWvZ+kqBSMFo6YUiGQs3P2oFlf9Q5Pb3g7s0GI5h9T9vH+n4yhV+yykL3ORWIU2ETuvVpraaRV5PBcKRJKWUoWiTqmaUFxqjY5dSBuHvDaSZQ9GeTx+BSJJU2FCkY5uHbYQW4kKRmMXUtVIwShbCkRSRIUJRSKSn1A2GHkGozKFIwUiKarShaJO5xWphVZsqhYlJ5QNR54bz6IHo7KFPwjnfS3ZKF0oyppaaOFTMEpOKO00BaPk5fm48ng9Q3kvS7YKFYpCnFckcVAwSlYIG5O8g1FRwlHejyWvQCTlVKhQ1K6sd81XtSgOCkbJCmHDkvfxy2IPRnmPX4FIshZVKBoeb7msV/OKkhHiHmiNFIySFUILIoRglHe46EbeY1YgkjxEFYpEsqBglDxtbOIJRyGMU4FI8lK4UNTuvCK10GQkCkbJy3Ojk3e1qF4IoaOZUMalQCR5ii4U5XXOIbXQykfBKHkKRi8IJYSEMg5QIJL8RReKYqZqUXwUjJKnYLShPEJJbZ2hhCFQIJIwFDIUadf8cMUw2bqRglHyFIw2lkVQCS0I5UmBSJopZChql45uLe1SMEqeglFrSQSk+mWEHoayfj0UiKSVCXkPoBtrpo0P+gM+knEzpkdZLZEXgpHCcXJ6lq7KfKeH9eteNhRFVTnW77p2KRBJSApbKYrhy07ipKpRsvI8llHoFaOiUyCS0BQ2FLUrr79SJW62dEDhKGEKRuWi511CFG0oinnXfO2FVhwKRslSMCoH7WkmoYpyTpHEbXjJs4UKhpprlKy85hnFMscodgpEYfN1xtDy3ryHkZtoK0XtCPXo1qBqURHVWmqqHo2dKkbFpEAkoSt0KEqLKgIyGgWkeCkYpUOBSGIQdSjKa16RSCcUjrqj4xgVhwKRxCLqUNSOkOcIqIVWLvXVIwWl9igYxU+BSGJS+FDULh3dWvKgoDQ6BaN4KRBJbKLf+yzmo1uLNNMqGJU5iOvI1/FRoJQYlaJSFPIXWllbaDrVSeeaVZXKVFlSxSgeeT1fMVSJyvSZjVH0laIkDfVP7uhD5f19eoNL7pq9B4taVVLFKHwKRM1pWxGHQlSKYt8LrazVIklPkatJqhiFS4FoY0X8DBaZKkUiJdD4pRx7JUkVo/AoEG1MYSg+hagUtSOto1sntXFRtUiyVIQqkipG4VAg2lDsn60yK00oEpHmYv4CVzDKnwLRC2L+LElFYUJR7POKQNUiyVesX+h5B6Myh6MyP/ZGMX52ZGOFCUXtSGseQOzzM0Tqxdhay7tqUMZwkOdjzvv1bhTTZ0VGlmsoMrNTzOwBM7vfzK4ys0l5jqcmrwmcoGqRhCW2cJSnMgUjBaIX6PNRLLmFIjPbGvgoMMfddwPGA0eOdB8fpUNWhBZamegAjvGI4Ys/hI1l0YNR3u3CEF7jGv3BUEx5t88mAJPNbAKwCbAo5/F0TS00KboYNgIhbDTzDg5pyfsxhfDa1oT+OZDu5RaK3P1x4ELgb8ATwDJ3vynt9aa1a36S1EKTkIUejkLZeOYdIpJUpMcyViG/92NmZt8ys8Vmdn/dZf9iZo+b2YLqv0Oqlx9dd9kCMxs2sz2bLHMPM7vdzO4zs+vNbOpo48izfdYPHA68BNgK2NTM/m+T251oZvPNbP7aVSsY7LOshyoiTYQcjhSMkhFK1SuU1zPU93tBXAEc3OTyL7n7ntV/PwVw9ytrlwHHAI+6+4Im970U+KS7/z1wLfDx0QaRZ/vsDcAj7v60uw8BPwJe3Xgjd5/r7nPcfc6EyZuOutA85xUl2UJTtUhiEeqGIpQNaSjBolOhjDmU1zHU93lRuPutQDcTTY8Crmpx3U7ArdWfbwbeNtrC8gxFfwNeaWabmJkBBwAPZrHiGFpoIjEJdYMRygYVwgkZ7QhlrCG8fiFXREviw2Z2b7W91t/k+nfSOhTdDxxW/fkIYNvRVpbnnKI7gWuAu4H7qmOZ2859y9JCU7VIYhLqxiOEDWtN6FWjkMYXwusW4vs5UpvXpsFU/53Y5v2+DrwU2JPK3OMv1l9pZv8HWOnu9ze5L8DxwIfM7C6gDxgcbYW5nhDW3c8Gzk56uWumjWfisnVJL7Yt3t+nD5KUmi0d0N6Yo6gFj1BOLBtKEKpRIMrRsDF+eeLTUJ5x9zmd3sndn6r9bGbfBG5ouMmRtK4S4e4PAQdV778j8E+jrTPvXfJzE0sLTdUiiVFoG5QQNrLN5B1GQqoMhSS0929Zmdnsul/fQqUdVrtuHJWW2PdHuP+sutueBXxjtHXmWikSkeIKrWLUs3RV7n/kNNp0Si9HvmN3Zs/uw8zwjP5MteFs1tOVYc919TaczJMz7PDUX5dyzVduY8Wy1Ykss8jM7CpgXyqttoVUukj7Vne1d+BR4P/V3eV1wEJ3/0vDci4FvuHu84GjzOxD1at+BFw+2jiiDUWDfUbvQOsPj1poIvmrfQ5CCUehBaMjj9qdXXd7CRMnbkplf5MX+Phk507aunzDRjtsbc5pbV1y2wx3Z8b0Gbz9o/Dtz/4yseUWlbsf1eTiy0a4/S3AK5tc/t66ny8GLu5kHKVtn4FaaCJZCemPhJBaabNn9zUNRFAJMbV/3ai/vwJRGxIMRABmxsSeyWzx4mY7TEmoSh2KYqJgJLFTMNqYmTUNRBvdriHgtPMvJkULRDVmxrhy7CxdGFGHopB3zQ+lXSASEgWjcjvmmLdx332/a/v2nzzzFH5+U+MOR91b+PhjXH/jtet/v+/+33Huv56V2PIlfm2HIjM70My+WTu/SAfHGchNkke3zruFBqoWSTGEFIxi09c3ie23n8mOO27J9tvPpK9vUt5DGrMsq0SPP/4YN9z44/W///1ue3DWJ8/JbP0Svk4mWn8QOA44y8ymUzmYUvSGpvVol1SRjIWyZ1poE69H0tc3iS23nMa4cZW/ZXt6JrDlltMAGBjobu+mlStXcsop/48nn3yC4eF1fOADJ3PIIYfzta9dxLx5N7NmzWr23HMOn/nMFzAzjjnmbeyyy2488MC9PPvsEs4//yvMnftVHn74Id74xsM4+eTTWbjwMd73vnex++578eCD97Pddttz/vkXM3nyJhus+7bbbuGSr17I4OAg2277Yj537kVsuknrUzndfsdtnH/hZ1m3bi277bon53z6PHp7J3LvfQs47/Nns3LVSnp7e7nisv/kueeW8okzTmLVqpUAfOqfz2Wvl8/hi1/+HH/+y584/G0H8ZbDj+BlO76Mb317Lv9+yeU8t+w5/vnsj/PYwr8xedJkPvPpz7Hzji/jq1//EoueeJyFjz/Goice5z1Hn8C7jz6OlStXcvInPsiTTz3J8Lp1fPDEj3LIwW/q6nWQcHTSPnva3Z9z99OoHAxp75TG1JGytdBULZKiCKViFEsbbebMvvWBqGbcuHHMnNn998yvfz2PWbO24Cc/+SXXXz+P1752PwCOPvo4rrnmZ1x//TzWrFnNvHk3r79PT08P3/vetRx55Lv54AeP41OfOo/rr/8vrr32ByxdWjl11SOP/Jl3vvNorrvuV0yZMoX/+I9vb7DepUuX8I2vX8zl3/w+1179c3bbdXcu/3brExqsWbOaT555Cl+68Otcf+2vWLduLf/xn99lcGiQUz7+Qf75k+dw3Y9u5opLv8+kiZOYMX1zLv/mf3Dt1T/nSxf+G+d+7tMAnHryGczZ6xX85Ic3cezRx2+wjq/+20XssvOuXH/NLzjlIx/n9LM+tv66Rx79M5d9/TtcfeV1fO3fv8zQ0BC//s0tzJq5Bddd/XNu+NHNvPY1r+/6dZBwdBKKbqz94O6fBL6T/HCSV7QWmkiRKBi1b8KE5t9lrS5vx4477sxvfvNrLrzwXObPv5O+vqkA3Hnnb3jHO/6JN71pf+6443/405/+sP4+++9/0Pr77rDDjsyatQW9vRPZdtsX8eSTiwCYPXsr9trrFQAcdtjbuPvu326w3t/ddRd/+vPDHHXMmzn8bQfx459cw6JFC1uO85FH/sI227yIl2y3PQBvOfwI5s+/g0ce+TMzN5/F7n9faVxMmdLHhAkTWLt2iLPO/gRvessBnPSx9/Pnvzy84QKbTKy+657/5fBD3wrAq/7Pa3juuaUMDDwPwOtfuz+9vROZ3j+d6dNnsOTZZ9jx73bmN3fcxgVf+hzz7/7t+udO4jZq+8zMvgyc4u4/qb/c3b+a2qgylmYLLY1jFo2bMZ3hJd2cTFgkPGqltWft2nX09Gz8lb12bfd7Tr3kJS/lhz/8Obfe+l9cdNF5vOY1r+e97/0gn/nMGVxzzc+YPXtrvvrVC1mzZs36+/T0TATAbNz6n2u/18bSuEdd/e+2bhjHec2rXsdFF3ytrXE6zfemc/eme+9d8Z1vsvmMmfzkhzczPDzM7v/w0vo7tVjHxmrL7u3tXX/Z+HHjWbt2LS/Zbnt+9P0b+e9f/xdfvPh8XvOq1/Hh95/U1uORcLVTKVoOXGdmmwCY2UFm9j/pDktGozaaFIkqRqN7+ukBhhuOtjw8PMzTT3f/3D311JNMnjyZww57G8cf/35+//v71geg/v7prFixgptuunGUpWxs0aLHueee+QDceOOP11eNavbcfS/uvud/+evfHgFg1apVPPLoXzZaTs32L3kpjz/+2Prb/+T6H7L3nFey/fZ/x+Knn+Le+xYAsHzFctauXcvA8gFmzpzFuHHj+Mn1P2RdtTK06eRNWLFyRdN17L3XK7iuumfanf97O/2b9TNlSuuw/tTip5g8aRKHH/pWTnjPifz+oVbnJJWYjFopcvezzOxdwH+b2RpgBfDJ1EfWgZCPbi0i7VHFaGS1ydQzZ/YxYcJ41q5dx9NPD3Q9yRrg4Ycf4oILPsu4ccaECT2cffbnmTp1GkcccTSHHXYAW2+9DbvttkfHy33pS3fgxz++mrPPPp0Xv/glHHXUuwGwapVm+vQZfO5fv8THPv5hBgcrIezkj35ifXus0cSJk/jcuRdx0sfev36i9VHvPIbenl6+dMG/ce7nPsXq1auZNGkSl1/6fd515Lv5yMkn8vObbuD/7P1qNpm8Caxbx0477Mz48eM57IiDeethb+dlO++6fh0f/sApnPHp03jT2/+RyZMm8/lzLxr5ufvjQ3zhS+cxbtw4JkyYwL+c+a8dP08SHvMWpcT1NzA7gMqJ1AyYDRzm7n8Y8U4p2WTWtr7TER9ret1IoQhoKxS120Lr5q/JtP4Sjr2NpoqXNAohGGUVij519v7M3urFmawrKwsXPsYHPvBurr9+3gaX53qAxpQOztiOJ556jM+fcPUGl/38sYvv6uas8VmY+KJtfevTTk50mY+cdFqwj7dRO+2zM4FPufu+wNuB/zSz/VMdlYiUVgittJDbaDHK/YjVIm0aNRS5+/7uflv15/uANwLnpj2wToW8az6k99evKi1SRApG8dpmm203qBLlHohyrBJJfDo+zYe7PwEckMJYRl/3GPaub2fX/FhOECtSBgpG8VMgkth0de4zd9c3RUBULZKiKnowcndGm9cpXQogELk7w3p5oxL1CWFjk+YE0hiDUYxjluwVORg98cQAa9asKGQwyr1KlDN3Z83QKp7669K8hyId6OTcZ0EY7IPeFt+RWe6aP9Q/WaV1kRJJY1f97191L0ceBbNn9zU9CGG0ci6P2HD+gWzY4am/LuWar9yW91CkA9GForSlfYLYNI5wXaMjXUtRFfUYRiuWD3LZN+cntrwQ5P3HYgiVRYmX2mcFo5aUFFUoG7u8N/ohy/u5CeU9IvGKslKkFtrIVDGSoipqxagI8v4uVCBKhq2DCc+Xt15S3kc+gnZ3ze9WFl/qqhhJUYWy8cs7BIREz4UUhUKR5EKhTcZCwSgcITwHobwfJH6FDEVZHt065BK6godI+kIIBXkJ4bErEEmSog1Fg2PoQCV5dOtuZTUvQsFIiiqkjWEI4SBrITzmkN4DUgzRhqKQhFwtgvCCUWjjkXiFtFEMISRkpUyPVcqlsKEo9BPEQnbVIpEiUzDKViiPMaTXXYoj6lAUewstS6FUZ0IZhxRLSBvIUEJDGkJ5bCG93lIsUYeikITeQgMFEpGshBIeklTExyTSSKEoZ1m30MbNmJ5bOFIokzSFVj0oUogI6bGE9jpLsUQfikZqoY02r6hsLbR6WQcUBSLJQmgbzJDCRDd6lq4K6jGE9vpK8UQfikLSbQstrwnXCipSRKFtOEMLFu0Kbcyhva7dGF7yrE7BFDiFopLLIhgpfImEFzJGEtNYY6AwFI/Ch6JYWmh57p6fZmhRIJI8hFpViCFshDjGUF/PkdSCkMJQXAoRisaya37SYtgLrZk0wosCkeQp1A1pqO20UMcV6uvYioJQ3AoRiiQZSe6ZpkAkIQh5gxpKAAk1DMVIYSh+UYUi73K0aqF1phaOug02CkQi7ck7kIQehkIOtfVUHSqOCXkPICmDfdAbyOdnqH9y8F827aoFnJE+8ApB4Wv3C7uIr6UtHQjmD45Wat8XWbXfY/h+iikQSXEUJhQVhff3BfllUMSNZZF1+0Xd7H5FeO1jCEawYVhJIyDFEIZAgaiMzOxbwKHAYnffrXrZBcCbgEHgz8Bx7v5c9brdgX8HpgLDwN7uvrphmZ8FDq9evxg41t0XjTSOqNpnAENTurtfLC00kW6ltbdLUVoDsWxoa5JordWWkXebrmiK8pkIzBXAwQ2X3Qzs5u67Aw8DZwCY2QTge8D73X1XYF9gqMkyL3D33d19T+AG4NOjDaJQlaKitNBCrRZJmLL6cq6tpwiVo5iM9D3SrJoUe/gJ/btPYSgd7n6rmW3XcNlNdb/eAby9+vNBwL3u/rvq7Za0WObzdb9uCvho4yhUKBIpk7y+nOvXG1tAiqWN1q7YA1AjBSIZwfHAf1Z/3hFwM/sFMBP4vrt/odmdzOxfgXcDy4D9RluJQlGHhqb10LOsWZUuWaoWSSshfTEPL3lWwUhKIaTPXZpsGHqWJ77Yzc1sft3vc919bttjMjsTWAtcWb1oArAPsDewEviVmd3l7r9qvK+7nwmcaWZnAB8Gzh5pXdHNKYJ85xV1ItYDOUq4Qvxi1vwKSULIfwTq/T1mz7j7nLp/nQSi91CZgH20u9faXwuB/3b3Z9x9JfBTYK9RFvUfwNtGW1+UoWgkIR3dWiRJoX8xhz6+eiFvgMso5Ncjpvd10ZjZwcDpwGHV8FPzC2B3M9ukOun69cDvm9x/h7pfDwMeGm2dhQtFWehkL7SxVItU4peaWL6YY6oahbwhljDE8l4uAjO7Crgd2MnMFprZCcAlQB9ws5ktMLNvALj7UuAi4H+BBcDd7n5jdTmXmtmc6mI/b2b3m9m9VCZnnzTaOKKdUzQ0pbu+52Cf0TvQegL6mmnjmbhs3RhGJpKcWL+UY5lrpPlF+Qs1nMb62YuVux/V5OLLRrj996jslt94+Xvrfh61XdYo10qRmW1mZteY2UNm9qCZvSqJ5RaphaYv7PKK/Us59vFL+hSIJDR5t88uBn7u7jsDewAP5jyeVGjCtXSqKF/KMbTTQt0wSz5Cf79KunILRWY2FXgd1fKYuw/WDt/drrIc3VrVonIp4pdy6I9JwSh7IT7nob9PJX15Voq2B54GLhvDOgAAIABJREFULjeze6qTozZNauGhtdBULZJ2FPlLuciPTTqjQCShyjMUTaByXIGvu/vLgRXAJxtvZGYnmtl8M5u/bsWKrMc4IlWLJEll+FIO+TGGuKGWbIT8vpRs5RmKFgIL3f3O6u/X0OTgS+4+t3bAp/GbblxIiuVAjiISPgWj9Ok5lpDlForc/UngMTPbqXrRATQ5+NJYFK2FpmpRcZXpL9XQH6s22ukJ8bkN/f0o2cr7OEUfAa40s17gL8BxOY+nY1mdC02Kq4xfyrEcx0iSo0AkMch1l3x3X1Btje3u7m+uHqUyM3m00FQtknpl/lIO+bGHuAGXZIX8/pP85H2cokSMNK8otBaaiLwg5A2TglFy9FxKLKIKRR7o3Ocs90IDVYuKIuRAkKWQnwdtzMcuxOcw5Pec5CuqUJSHEFtoEj99KW9Iz4dkRe81GUlhQlFau+aHStUiKZpQN1YhVjpiEdpzF+p7TMIRXShaO3W44/tkMa8o6xaaxEtfzK2F+tyEtnGPgZ4ziVF0oSgPobbQVC0SEWlPqIFbwlKoUJR3C03VIhmNvphHF+pzpMpH+0J7rkJ9T0l4ogxFobbQOqVqkUhzoW7EQtvYh0jPkcQsylCUB50LTcYq1A19qEJ9vrTRby3E5ybU95GEKe/TfARjsM/oHfAxLyeP0354f1+QX0YiInlSIOqcrYPeEm9OClcpiu3o1jpmkUhroW7U9EfIxvScSBFEG4q6mVc0ViG30DS3KGyhbtxjEOpzpxDwghCfi1DfNxK2aENRGvLaCy2papGCkRRVqBu4EMNA1kJ8DkJ9v0j4ChmKut01X0RERMor6lCUx675abXQVC0qLv3VmoxQn8cQKyVZCfGxh/o+kThEHYrSUIQDOSoYSVGFusELMRykrYyPWYqvsLvkD02Bw/9uJ07+x32YvVkfTzw3wJd/cRs/XfCHvIdWevsevhfHnn4oM7fq5+lFS7ni/Bu45Sd35z2sUtvvHa/iuHOOYOY2M3h64RIuP/tq5v3g9ryHFRVbOpD4HyT7vnF3jvvoG5i55WY8/eRzXP6VX3LLz+5NdB1FkkVo1mel2KIPRWunDjPh+Y0LXofuuhOf+acDmdxbqdhs3T+Vz7z1QAB+uuAPYzoOw5pp45m4bF33C2hhqH8yPUtXJbKsUI9dtO/he3HS+UcyaZNeALbYZjonnX8kgIJRTvZ7x6s4+ZLjmbTpRAC2eNHmnHzJ8QBBftkPL3mWcTOm5z2MppIMRvu+cXdOPvtwJk2ufla26ufksw8HyD0YhfjdklUgiumzIp0rbPvs1P32WR+Iaib39nDyP+4z6n2L0EKDMNtox55+6PpAVDNpk16OPf3QnEaUrlDbPfWOO+eI9V/yNZM2nchx5xyR04jillRgOO6jb1gfiGomTe7luI++IZHldyvEQJQVfVaKr7ChaPa05oFg9mbhBYV6RT+Y48yt+ju6XNI3c5sZHV0eghjC5ljN3HKzji7PQqiBKKv3Q4yfFelMXKFoXPPTcDTbC+2JZc0/vE88V7k8q73QVC3a0NOLlnZ0uaTv6YVLOro8FCEHoyTCw9NPPtfR5WkreyCCeD8r0r64QlEHvjjvNlYObngOslWDQ3z5F7e1df+kWmjdKHK16Irzb2D1ysENLlu9cpArzr8hpxHJ5WdfzeoVaza4bPWKNVx+9tU5jagYxhoiLv/KL1m9quGzsmqQy7/yyzEttxuhBqKs6bNSfNFPtG7lhgcqe5mdut8+zJ5W7r3PQpp0XZtMrb3PwlGbIBrjHjUhT7qGsU28rk2mznvvs1C+O5rJuloY82dF2mPuYz8zfFYmbreNb/2xU5pe12wPtJqe5a2XOdJeaL0Doz837e6F1rNsaPQbNd4noT3RakL+ciuqkFs8RRJyMILw2tjtCv07I9bP100rv3uXu8/JexzNbDJrW9/piI8luswF//axYB9vo+jaZ+umNA8haZwgNs8WmogUR+jhIkaxBiIJW3ShKFZ5T7iGeP9ajZW+tLMTw3MdWzCKbbwiSShFKBrpBLGhngsN0plwrWAkkp9Ygkbo44whBEucogxFaqGNjYKRFFEsG8rQA0fo4xNJU5ShKDRpHrOoyLvniyRNwah7tnQgyHE1iuU1ls6Y2Ulmdr+ZPWBmJzdcd5qZuZlt3uR++5nZgrp/q83szd2OozShKM0WWoxULRLJV0gBJKSxSPmY2W7A+4BXAHsAh5rZDtXrtgUOBP7W7L7uPs/d93T3PYH9gZXATd2OJdpQVKYWWlrVIgUjKaKYKgkhhJEQxtCumF5b6cjLgDvcfaW7rwX+G3hL9bovAZ8A2jl+0NuBn7n7ym4HEm0oCk0sp/1opGAkkq8821YxBSIptPuB15nZDDPbBDgE2NbMDgMed/fftbmcI4GrxjKQwh7RupmhKa0P5DjYN/KBHPM21D858YM5ihRV6Ee6bmYsR7/uZl2xUZUoG7auvQMXd2hzM5tf9/tcd59b+8XdHzSz84GbgeXA74C1wJnAQe2swMxmA38P/GIsA426UhRrC03VIhFpJouqUYyBSKL3jLvPqfs3t/EG7n6Zu+/l7q8DngUeBV4C/M7MHgW2Ae42sy1brOMdwLXu3vnpI+pEHYpCk+YxiyDdPdEUjKRoYq4spBGOYtm7rJmYX0tpj5nNqv7/IuCtwHfcfZa7b+fu2wELgb3c/ckWiziKMbbOoIShSHuhtaZgJBKWsYaYWhCKNQxJqfzQzH4PXA98yN2Xtrqhmc0xs0vrft8O2JbKBO0xiX5O0bop6xi/fOMKzdqpwyOeJLYbg32WWK91aFpPVyeJ1dwikfbFOLeoUbNAM9IfMEULQKoSlYO7v3aU67er+3k+8N663x8Ftk5iHNGHotCsmTaeicuaz3WKgff3Fe5LVaRo9BkVSUfp2mcQTgut2wnXaR/lWm00KRJVGuKl106yFlUosvGJ7ybYsbwP5JgVBSMRESmbqEJRK1numt+OtPdCg2zOiaZgJEWhikN89JpJHgoRiroxUgstS6Eds6iRgpGIiJRFdKGoZ8pg6usYbV5ROy20olSLQMFIikGVh3jotZK8RBeKWgmthdaJ0KtFoGAkIiLFV5hQ1I1QWmhjkVW1CBSMJH6qQIRPr5HkKcpQVMQWWgzVIlAwEhGR4ooyFLUScwttLLKsFoGCUbtiP5JyUakSES69NpK3QoWiboRyIMfYKBiJiEjRRBuKsmihjSakFlrW1SJQMJJ4qSIhIs3kHorMbLyZ3WNmN6S5nqK30PKiYCQiSVBQlRDkHoqAk4AHk1pYq3lF3UqihVbkahEoGEmctBEWkUa5hiIz2wb4J+DSPMcxll3zy3IutNEoGIlItxRQJRR5V4q+DHwC6Kq31em8olhaaDFWi0DBSOKjjbGI1MstFJnZocBid79rlNudaGbzzWz+uudXtLXsrFtooZz2o0bBSERioWAqIcmzUvQa4DAzexT4PrC/mX2v8UbuPtfd57j7nPFTN01tMKEd3TqWgzk2o2AkIiIxmpDXit39DOAMADPbFzjN3f9vp8vpmTLI0PLehEeXjjXTxjNxWbJVrFaG+ifTs3RVJutqxvv7sKUDua0/FONmTNdfwoEbXvKsDrQpUjVunWe2nQpR3nOKUpP00a2TaKGVjff3qWokIi3pDwYJTRChyN1vcfdDR7vd+HHpTpQuWgstz7lF9RSMJHTaOIsIBBKKxiqEo1u3K8sJ16BgJCJhUhCVEBUiFLVShBZazBOuGykYSci0kRaR6ELRtL50Jw9n0UIra7UIFIxERCRc0YWiVmJqoXWqSNUi0ARskbJTVU5CVZhQ1EoRWmhJCKlaVFOWYKTdveOhjbVIuUUZisrYQitataimLMFIRCoUPCVkUYai0Kla1BkFIwmJNtoi5VWoUNRqXlHWLbS0JFEtUjASERFpLtpQVMYWWtEpGIkUm6pwErpoQ1Ho0mqhFblaBMUNRppsHRdtvEXKqXChKLYWWl7VIgUjEcmSgqbEIOpQVIQWWjeS2hMt9GCkcCR50kZcpHyiDkWha7eFprlFrSkYiYhIVgoZijo9unXeLbRulKFaVKNgJHlRtSgZeh4lFtGHok5aaK3mFY0k1BZakmIIRiIiImmLKhRNGN9dRSdNSZ32Qy20kRWhWqQ90EREwhZVKOpEVi20PCV56o8YqkVFCEYSH7V+xkbPn7TDzDYzs2vM7CEze9DMXmVmnzWze81sgZndZGZbNbnfnmZ2u5k9UL3tO8cyjsKGolbyaKGpWpQcBSMRkUK6GPi5u+8M7AE8CFzg7ru7+57ADcCnm9xvJfBud98VOBj4splt1u0gChGK0t41fzRFmHANcVSLQMFIRKRIzGwq8DrgMgB3H3T359z9+bqbbQp4433d/WF3/2P150XAYmBmt2OJLhTN2nR527dNsoVWlmqRglG6NK8oTmoBdUfPm1Rtbmbz6/6d2HD99sDTwOVmdo+ZXWpmmwKY2b+a2WPA0TSvFK1nZq8AeoE/dzvQCd3eMWbrpqxj/PLitKqGpvXQs2woueX1T6Znab7Vt3Z4fx+2dCDvYYiIFIat80S3J1XPuPucEa6fAOwFfMTd7zSzi4FPAp9y9zOBM83sDODDwNlNx202G/gu8B5373qScHSVolaK1ELLu1oEqhiJiEhmFgIL3f3O6u/XUAlJ9f4DeFuzO1fbbzcCZ7n7HWMZSJShqOgttG4kObcoNgpGkgW1gjqj50va5e5PAo+Z2U7Viw4Afm9mO9Td7DDgocb7mlkvcC3wHXe/eqxjiTIUJaGbvdBCl3QwiqVaBApGIiKR+whwpZndC+wJnAd83szur152EHASgJnNMbNLq/d7B5VJ2sdWd91fYGZ7djuIQs0pmta3imUD+W3IB/ugN6EpLmumjWfisvyDWyzzi2IybsZ0/RUtIlLH3RcAjfOOmrbL3H0+8N7qz98DvpfUOKKtFKmF1lwabbRYKkaqFknaFGbbo+dJYhVtKEpCGi20ok24jo2CkYiIdKtwoSjvvdDaoWpRuhSMRESkG1GHok5aaJ1Ks4UWq5iCkUha1BoamZ4fiVnUoagTreYVFbWFltYu+rEEo9CrRTqytYhIeAoZitRCqyjzsYsg/GAkIiJhiT4UlaGFFtqE61iqRaBgJOlQi6g5PS8Su6hCUe+4sbW6QmuhxVwtiikYiYiItCOqUNSJGFponQitWgTxBKNQq0WaVyQiEpZChKIytNDGIs25RQpGUlZqFW1Iz4cUQXShaNtNnhvT/WNuoYVYLYJ4gpGIiMhIogtF0p2y74kGqhaJiMjICh2KkppXNJYWmqpFMhLNK4qbWkYixRJlKGrWQkviBLFptNBCkna1KIZgpGqRSPIUDqUoogxFeUh7wnVW1SIFIwUjERFprvChKIRd85M87UcMYghGIiIijaINRWVvoYVcLYLwg1FI1SLNK4qbWkcixRFtKMpDSBOux0rBSESSoFAoRTIh7wFkYVrfKpYNFG8DvWbaeCYui6OyFSLv78OWDuQ9DBGRYNg6p2dp/tNO8hJ1paiILTRVi7IVUhtNRETyFXUoykPaLbROhXrconqhB6MQaF5R3NRCEimG0oSiEPZCa1fRqkUQdjBStUikOwqDUjRRhaJeW7vRZbG10EKsFikYiYiIRBaKQpH2gRwh22pRlkINRiFUi9RCExHJV26hyMy2NbN5ZvagmT1gZielvc6YWmidiqVaBOEGI5GxUCtJJH55VorWAqe6+8uAVwIfMrNdRrvTiycvaWvhebbQ8phwnYSyB6MQqkUiIpKf3EKRuz/h7ndXfx4AHgS27mZZzeYVpW2kFlpSOm2hxbAnmogUgypjUkRBzCkys+2AlwN3pr2uUFpoqhapWtSM5hWJiOQn91BkZlOAHwInu/vzTa4/0czmm9n8FUsrba4itNDalUe1qOzBSEREyinXUGRmPVQC0ZXu/qNmt3H3ue4+x93nbNrf23JZRW2hlUFowSjvapHESy0lkbjlufeZAZcBD7r7RVmuu+gttNiqRRBeMMqTWmgiIvnIs1L0GuAYYH8zW1D9d0i7d1YLLX1lDkaqFomIlM+EvFbs7rcBiW7tt93kOR5buVmSixzV2qnDTHi++2w52Ae9KZyofc208UxcNvajcg9N66Fn2VACI2pzff2TgzlDs/f3YUtTeHFEIqc2oRRV7hOt89JpCy3GalGsu+iHVDESEZHyKEUoSqKFNpKiT7jOuo0WkrzaaJpXJCKSvahDUbN5RXnshTZW7U64zrNaVOb5RSKdUGtJJF5Rh6Kxiq2FVjYKRiIixWdmk8zst2b2u+q5UM+pXv5hM/uTmbmZbT7C/V9kZjdVz6X6++oBobsSVSiaNC6bCb95tNBULWqxzgCCkVpoIiKpWgPs7+57AHsCB5vZK4H/Ad4A/HWU+38HuKB6LtVXAIu7HUhUoaiZdltoncwrSkNRqkVlDUYiIpIOr6htpHuq/9zd73H3R0e6b/VE8hPc/ebqspa7+8puxxJ9KBqrrFpo7YihWlRWOm6RiEh6zGy8mS2gUuW52d3bPRfqjsBzZvYjM7vHzC4ws643eNGFor+b+FQm69FeaK2pWpQdtdBEpAA2r53DtPrvxMYbuPv/b+/eY+SqzzOOfx9srxfb+H6pwTTEKipBhNSIImhUBwqobrFwK4XIlWqZS2VVogGqRJAUKS5pq9Zq1bSIqJJLAshBaYNDGmQlBYc68T+FhkBjSEwa5ApYoNgGx1fWN97+MTPLej3rPTt7zvnNOfN8JMTOYXb2PWZn9/Hv/V1ORsSvAUuAKyRdkvG1JwO/CXwW+HVgKXBzp4Um27wxTx86+x1efW/eKdfabeS4cPohdh9O18c6PgOm5NTFO3aO6DsY4/qcvDZ0hPI3dYT0Gzt6M0czr66rvRMni/g5tzciLs/yxIj4haTvAyuAlzJ8ygDwQkTsApD0b8CVNI4RG7fKjRQVIc8WWlkTrjuVZxvNI0ZmZjZRkhZImt38+Gwak6tfzvjpPwTmSFrQfPxbwE87raWSoaibW2hjyXPCdaoz0VJLGYw8t8jMLHeLgW2SdtAIOVsjYoukOyQN0Gip7ZD0IICky1sfR8RJGq2zpyW9SOP4sH/utJBatM+gOi20LIo6D62l6m20XnPWvLluWZhZbUXEDmBZm+v3A/e3uf4c8EfDHm8FLs2jlkqOFBWhm1po49ENo0W91kbzaJGZWT05FI0hVQutSnOLoPeCkdmZeGTPrJoqG4razSuqykaOeet0tKgOexelCkZljxZ5ab6ZWfEqG4qKMFoLbbTRoom00Oqyw/VwKUaLzMzM8uJQ1MXG00LrltGiXmqjeW6RmVm9VDoUZV2aX2QLrchjP6qql4JRmdxCMzMrVqVDUTvt5hWNx3hbaGeSRwutiqNF4GBkZmbVU6lQ1KcTqUuorTpMuk7BLTQzs/qoVChqZyK7W5fRQstjz6IyRouK4NGi/LmFZmZWnMqHonayLs0fTZ4ttLF00yo0t9E649EiM7N6qFwoWtq3J3UJSVR1tAi8VN/MzKqhcqGonaq30Oo+WgTlB6M6jxa5hWZmVoxahKJ2qtRCy6qs0aK6TLqu+/wiMzPLVyVDUbe20LpttKjbglHd22ieW2QtHs0zq6ZKhqJ2uqGFVoaiD4otWt3baGXxL11Lxd97Vme1CUXtVK2F1gujReBgZGZm3amyoaiOLbSsyhwtqsv8orK4hWZmVl2VDUXt9EoLbby6bYk+eLQoD25jmJnlq1ahqJ1ua6FVaXl+i9to4+PRIjOzapqcuoCJWNq3h13HFqQu4zQnZ5xk0qFi207HzoG+g+N5vug7GB1/vaOzJjF1/+itwU4dnzWFKfuP5/66veKseXN5/513U5dhZnVx8kRP/0yp3UhRHVpo3ThaVKQyR4w8WmRmZqOpVCjqV2cjHVVroWU13gnXE51b5EnXZmZWZ5UKRe1UcRVaFh4tKvBr1Wy0yBOurWz+nrO6qnwoaqeIFppHixqKHC2qYzCy3uPAYFZdlZtofeHkk/z8RLZfzBfP/ATXLFrLrCnz2X98L5sHvsl/vfvM0H8/f9oveP3I7KJKLWXCdSdST7q+bvlFrFuznIXzZ7J77wE2btrO97a/DHjidSrXfOoqbrnvJhYsmceegXd4aP1jbPvGf6Yuq+ddveoybr5nJQvOncOeN/fx8IYtfP/bz6cuq6f5vVJvlQtF7bRbhXbxzE9ww3mfpu+sfgBm9y1kzYduBTglGHWz4zNgSoa53+NdiZaHToPRdcsv4u7bV9Df3xgV+qWFs7j79hUAQ8GoLMfnnM2Ufe1HAPMSc85B+4r/nzORVWjXfOoq7nrgVvqnTwVg0S/P564HGu8V/7BP5+pVl3HnhtX0T+sDYNGSudy5YTWAg1Eifq/UXy3bZwDXLFo7FIhapk6ayu+fd9OYn1u1FlonUm3ouG7N8qFA1NLfP4V1a5YPPa77wbHd5pb7bhr6Id/SP30qt9w39nvFinPzPSuHAlFL/7Q+br5nZaKKzO+V+qtkKLpw8tgjFLOmzG97fW7fvFMej2cVWifKmnCd4qDYTuYXLZw/M9P1soJRGXOLun3C9YIl88Z13cqx4Nw547petl6cO+X3Sv1VMhRlceRE+8nW7x47fXl+SlUfLRpvMNq990Dm63UKRt1sz0D798Ro1210eQaFPW/uG9d1K57fK/VXm1A0cmn+jne+zLH3B0+5duz9Qb71xmOZXi/PFlo3jxaVHYw2btrO4OCpE6kHB4+zcdP2ts+vSyutmzdzfGj9YwwePnrKtcHDR3lofbb3ihXj4Q1bGDxy6s+VwSPHeHjDlkQVmd8r9VfZidZjrUJ77dCTAHxk7h1Dq8+2vf0Ib7/3A+D0FlqRq9DsA63J1KOtPkuljEnXZehkwnVrgqhX1HSX1mRqrz7rHn6v1J8iOl+aXbZlH+uLbd9dNPR4ZChqdw7aK0cXnXbt1fdO7/+OFop2Hz59mGb/wdHbLccP9bW9PtbS/MkHxh60y7ISDTpbiTaRJfotRZyN1lLGMv2iQ1EZq9CAnj63KLVem2fj77XOPHVk048i4vLUdbQza9K8uLL/hlxfs5vvd6TatM+gvN2tR2uhnclEW2hVUPWNHYueW9TNLTSzTvRaCLT6q3QoyrIKrd3u1hM9C60IWSZcd/vcojqow6Rr/6JKw3/uZtVX6VBUhtEmXI+myAnXRUuxGm086jDp2qNFZmbdq3ahqJtbaGPJc3l+in2LWqoejDxaZGbWmyofiurUQssiawsN0rbRqh6MiuTRIqsTB3Crk6ShSNIKST+T9Iqkz6Ws5UzKbKGl3MzRPlCH0SIrj4OB2cSMlQckTZX0r83//qykC4qoI1kokjQJ+DLwO8DFwB9IuvhMnzNpyqWZXrvKLbQsPFpU/dGiMvgXtZlVQcY8cBuwLyJ+BfgSsKGIWlKOFF0BvBIRuyLiGPAvwKpOXqhqLbRun3ANDkZQ7GiRW2hWJw7gNkFZ8sAq4JHmx5uBayXlvmw6ZSg6D3h92OOB5rWulOexH2PJc3k+pJ10DdUORlXnX1bF85+x2YRlyQNDz4mIE8B+Rh5PkYOUx3y0S3inbassaR2wrvnw0Jzz+Fm2l3+148JGmA/szevFupjvs158n2U5UspXSX+f5eiV+/zV1AWM5sD77z751JFN83N+2X5Jzw17vDEiNg57nCUPZMoME5UyFA0A5w97vAR4c+STmn9wG0deL4uk56qyPflE+D7rxfdZL77PehkRELpKRKxI8GWz5IHWcwYkTQZmAbmfM5OyffZD4EJJH5bUB6wGnkhYj5mZmZUvSx54Aljb/PiTwH9EAYe3JhspiogTkv4EeBKYBHw1In6Sqh4zMzMr32h5QNIXgeci4gngK8AmSa/QGCFaXUQtKdtnRMR3gO+krCGDZK27kvk+68X3WS++z3rplfvMrF0eiIgvDPt4ELip6DpUwOiTmZmZWeVU/pgPMzMzszw4FI2DpM9KCkl5L1dMTtLfSnpZ0g5J35I0O3VNearKkTITIel8Sdsk7ZT0E0l3pq6pSJImSXpB0pbUtRRF0mxJm5vvzZ2SrkpdUxEk/Wnze/YlSV+X1J+6pjxI+qqk3ZJeGnZtrqStkn7e/PeclDXaqRyKMpJ0PnA98FrqWgqyFbgkIi4F/gf4fOJ6ctPJkTIVdQL4TER8BLgSuL2m99lyJ7AzdREF+0fg3yPiIuBj1PB+JZ0H3AFcHhGX0JhoW8gk2gQeBkYucf8c8HREXAg83XxsXcKhKLsvAXdTwGZR3SAinmruEgrwDI19IuoityNlullEvBURzzc/PkjjF2jX7hI/EZKWADcAD6aupSiSZgLLaay6ISKORUTx5xClMRk4u7n/zDTa7FlXRRGxndP30hl+XMUjwO+VWpSdkUNRBpJuBN6IiB+nrqUktwLfTV1Ejip1pEwemidILwOeTVtJYf6Bxl9Sxj4Tp7qWAnuAh5ptwgclTU9dVN4i4g3g72iMwr8F7I+Ip9JWVahFEfEWNP4iAyxMXI8N41DUJOl7zX72yH9WAfcCXxjrNbrdGPfYes69NNowj6arNHelbA/fLSTNAL4J3BURB1LXkzdJK4HdEfGj1LUUbDJwGfBPEbEMOEwNWy3NOTWrgA8D5wLTJf1h2qqsVyXdp6ibRMR17a5L+iiNN+uPmwfyLgGel3RFRPxfiSVO2Gj32CJpLbASuLaInUITynSkTB1ImkIjED0aEY+nrqcgHwdulPS7QD8wU9LXIqJuv0gHgIGIaI32baaGoQi4DvjfiNgDIOlx4DeAryWtqjhvS1ocEW9JWgzsTl2QfcAjRWOIiBcjYmFEXBARF9D4QXVZ1QLRWCStAO4BboyIco64LE9PHCmjRmr/CrAzIv4+dT1FiYjPR8SS5vtxNY3t/usWiGj+jHldUuvw0GuBnyYsqSivAVfeo8WqAAABl0lEQVRKmtb8Hr6WGk4oH2b4cRVrgW8nrMVG8EiRtTwATAW2NkfEnomIP05bUj566EiZjwNrgBcl/Xfz2p81d4q1avo08GgzzO8CbklcT+4i4llJm4HnabTuX6AmOz5L+jpwNTBf0gCwHvgb4BuSbqMRCAvfpdmy847WZmZmZrh9ZmZmZgY4FJmZmZkBDkVmZmZmgEORmZmZGeBQZGZmZgY4FJmZmZkBDkVmZmZmgEORWU+TtE3S9c2P/1LS/alrMjNLxTtam/W29cAXJS0ElgE3Jq7HzCwZ72ht1uMk/QCYAVwdEQclLQXuBWZFxCfTVmdmVh63z8x6mKSPAouBoxFxECAidkXEbWkrMzMrn0ORWY+StBh4FFgFHJb024lLMjNLyqHIrAdJmgY8DnwmInYCfwH8edKizMwS85wiMzuFpHnAXwHXAw9GxF8nLsnMrBQORWZmZma4fWZmZmYGOBSZmZmZAQ5FZmZmZoBDkZmZmRngUGRmZmYGOBSZmZmZAQ5FZmZmZoBDkZmZmRngUGRmZmYGwP8Diu5A+RvEK1EAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1008x576 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def f(Xi, Xj): # Define Branin function\n",
    "    a, b, c, r, s, t = 1., 5.1/(4*np.pi**2), 5/np.pi, 6., 10., 1./(8*np.pi)\n",
    "    f_ = lambda x1,x2: \\\n",
    "        a*(x2 - b*x1**2 + c*x1 - r)**2 + s*(1-t)*np.cos(x1) + s\n",
    "    return np.reshape([f_(x1, x2) for (x1, x2) in zip(Xi[:], Xj[:])], Xi.shape)\n",
    "    \n",
    "# Define the function space\n",
    "Xi, Xj = np.meshgrid(np.linspace(-5., 10., 50), np.linspace(0., 15., 50))\n",
    "\n",
    "# We will also define some samples\n",
    "Xsi, Xsj = np.meshgrid(np.linspace(-4.5, 9.5, 5), np.linspace(0.5, 14.5, 4))\n",
    "X = np.vstack((Xsi.ravel(), Xsj.ravel())).T\n",
    "# Our observations\n",
    "y = f(Xsi, Xsj).ravel()[:, None]\n",
    "\n",
    "# Plot the Branin function\n",
    "plt.figure(figsize=(14, 8))\n",
    "\n",
    "plt.contourf(Xi, Xj, f(Xi, Xj), levels=np.linspace(0, 300, 20))\n",
    "plt.plot(X[:,0], X[:,1], 'wo')\n",
    "\n",
    "plt.xlabel('$x_1$'), plt.ylabel('$x_2$'), plt.title(\"2-D Branin-Hoo function\")\n",
    "plt.legend(labels=[\"sample locations\"])\n",
    "plt.axis('square'), plt.colorbar();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "\n",
    "Using `GPy`, fit a GP that is the product of two kernels over the two different inputs. Compare this fit to that of a GP with an isotropic RBF kernel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style type=\"text/css\">\n",
       ".pd{\n",
       "    font-family: \"Courier New\", Courier, monospace !important;\n",
       "    width: 100%;\n",
       "    padding: 3px;\n",
       "}\n",
       "</style>\n",
       "\n",
       "<p class=pd>\n",
       "<b>Model</b>: GP regression<br>\n",
       "<b>Objective</b>: 65.96007391289652<br>\n",
       "<b>Number of Parameters</b>: 5<br>\n",
       "<b>Number of Optimization Parameters</b>: 5<br>\n",
       "<b>Updates</b>: True<br>\n",
       "</p>\n",
       "<style type=\"text/css\">\n",
       ".tg  {font-family:\"Courier New\", Courier, monospace !important;padding:2px 3px;word-break:normal;border-collapse:collapse;border-spacing:0;border-color:#DCDCDC;margin:0px auto;width:100%;}\n",
       ".tg td{font-family:\"Courier New\", Courier, monospace !important;font-weight:bold;color:#444;background-color:#F7FDFA;border-style:solid;border-width:1px;overflow:hidden;word-break:normal;border-color:#DCDCDC;}\n",
       ".tg th{font-family:\"Courier New\", Courier, monospace !important;font-weight:normal;color:#fff;background-color:#26ADE4;border-style:solid;border-width:1px;overflow:hidden;word-break:normal;border-color:#DCDCDC;}\n",
       ".tg .tg-left{font-family:\"Courier New\", Courier, monospace !important;font-weight:normal;text-align:left;}\n",
       ".tg .tg-center{font-family:\"Courier New\", Courier, monospace !important;font-weight:normal;text-align:center;}\n",
       ".tg .tg-right{font-family:\"Courier New\", Courier, monospace !important;font-weight:normal;text-align:right;}\n",
       "</style>\n",
       "<table class=\"tg\"><tr><th><b>  GP_regression.         </b></th><th><b>                 value</b></th><th><b>constraints</b></th><th><b>priors</b></th></tr>\n",
       "<tr><td class=tg-left>  mul.rbf.variance       </td><td class=tg-right>    15308.189387215683</td><td class=tg-center>    +ve    </td><td class=tg-center>      </td></tr>\n",
       "<tr><td class=tg-left>  mul.rbf.lengthscale    </td><td class=tg-right>    23.192541106085837</td><td class=tg-center>    +ve    </td><td class=tg-center>      </td></tr>\n",
       "<tr><td class=tg-left>  mul.rbf_1.variance     </td><td class=tg-right>    15308.190782336731</td><td class=tg-center>    +ve    </td><td class=tg-center>      </td></tr>\n",
       "<tr><td class=tg-left>  mul.rbf_1.lengthscale  </td><td class=tg-right>    146.98147501349206</td><td class=tg-center>    +ve    </td><td class=tg-center>      </td></tr>\n",
       "<tr><td class=tg-left>  Gaussian_noise.variance</td><td class=tg-right>5.562684646268137e-309</td><td class=tg-center>    +ve    </td><td class=tg-center>      </td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<GPy.models.gp_regression.GPRegression at 0x2bed45d9be0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAHSCAYAAAAUvDo/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3de7QcZ3nn+++jvbd8kyzLNmaI7YlhwjWOEYzgOMOEBHOJhzAGchu8CDgDxCtzgBgODOB4ApNJMotghksWZ4Xl2IY4eMyZAE4yXGL7cEi8Mss4EbZsZOQAQwALG4RtIUu2kbak5/zR1ai11b27uqvqvdXvs5aW9u7dXfV2VXXVr5/3rSpzd0RERERytyZ2A0RERETaoFAjIiIiRVCoERERkSIo1IiIiEgRFGpERESkCAo1IiIiUoTFkDNba8f6cWtOCDlLWAj6FgcWF4LNyhcsu2kfanm63vLibmt63uJXhlbf45p2L+NgC91dFmJhzaHOpj3J4kL4eebowMH2vxMfPDT/NP3gDPuVQ/WeawfHP77v3h33u/tj6s8wnBc871h/4MH2t+Gtdy7f4O7ntz7hlgU94h+35gTOPfYXQs6SNaecHHR+Q75xfZD5LG88rtvpb1jqZLr7NrSbRPavbzco7W9p9S2va2c6AAdObG9HdXDdhL31nJbW7W91eqM2rH+0s2mv5rQT9kaZb+p2PtziRr3C7j3z7c+W966d6fkLe6fvfxYfmhywvvY7/9e3ZpphQA88eIgvfO6xrU934+k7Tm19oh1Q91PmlnZ1u8Nf2r3cyXSP2d3uQXXtnjQvIrnU4nFxtZ3srOrs1GexvHftzAeWunbvOW7ug10TOx9e1+kBPDddL4+UAo3kS6FGpuoq2LStzWCzdk9rk+pNsIHZDzCziBFs4PDBvK8BJ8R7j7VuJ2nzcyZhac11xHa1eFScoutqTVfartaAgs2sFGxm07dwE+K9NlmnqtLISgo1hVA3VDcUbObTdbBJJdyUGHBCvbem67GrQKMqTd609joUsloTQi7Bpu3xNX0JNjmNs4E0wg2UEXBCv4em663L7UryplBTkFy7obqgYDOf3Ko2EL9LatRoOEg95MRqZ4z11WaVps3Pr7QvwkVcJGdLu5c7Oc37mN0HWz/Ne+0eb/VU77V72jvVe2lve6d7Lz60ptXTvRf2LrR+yvfy3rWdnvY9PFDGOv17knGBIdap4imErDYCjcbRyGoUajpmu/YEu2YNDKo1XV+7ps/Bpk0KNu3bvee45ILNSiGCTgoBZqUYgWYWqtKUQaGmQCGCTVdSDzZtVmugn8EGur1YX6pVm9WkGELaFCvQqErTPxpTE0BpA4ah22vX9OmMKOjfGBsIM9AzlYHEfZfDOtAZT+XQmixUiEHDuVyUD9IeOAwKNl1SuImjzeWeQpVGXU95UKiRRnI5zRvSDzZtyinYKNyUp83l3HWgUZWmLFqbgcTogsr9FO++BZu2vwnmEmwg7HVHFGy603ZwTOV6NKrS5EOhpnC5d0Mp2DTTRbDJvTsKVLXpQtvLc97tQVWaftMaDajEAcNDOY2vAQWbpkrojoLD4UYBZ35dLL8QgaYuVWnyolDTA+qGCkPBpj0xuh0UbmbXxfIKte5VpSmT1mpgqtbMJ4duqLYp2MQZT6FwM11Xy6jJOleVRkChpjdCVWv6HmxSPtUb8gw2Cjfp6HKZhAw0qtK0y8zONLMvmNl2M7vLzC6pHt9kZl80s61mtsXMnr3idc8ys4Nm9stjpnm8mX3GzO6upvnuOm3Rmo0gVrVGwWY8BZvmur5ya8yzYDTupp8BT1WamRwA3uLuTwXOBV5vZk8D3gP8rrtvAt5Z/Q6AmS0AfwjcsMp03+vuTwGeATzHzP7NtIYo1IigYNOGLs+MgrhVm6G+HdxDvV9VafLm7ve5+23Vz3uA7cDpgAMnVk/bANw78rI3Ap8Edk6Y5iPu/oXq5/3AbcAZ09qitdszqtZMpmDTjpKrNkOlV29ChpnUxtFIM2Z2FoPKyq3Am4DLzewe4L3ApdVzTgdeDny45jRPAv4t8Plpz9UNLSMJfffuUaFueNnV3byhmxtfQvt39U75BpjQ/k0wh7q4GeaoEDfGrGv04J/TTTRXCh3QmobTeQJN3SCfc9fTD9342oFOwt6pZrZl5Pcr3P2K0SeY2ToG1Zc3uftDZvb7wJvd/ZNm9qvAVcALgA8Ab3f3g2ar72/NbBG4Dvgjd//GtEYq1Eincgw2bcsh2ACth5uugw0MDowpBJuhlcEg9ZATq9IUo9qmbqfG7nf3zZP+aGZLDALNte7+qerhi4BLqp//HLiy+nkz8PEq0JwKvNjMDrj7X4yZ9BXA19z9A3UaOXUtm9nVZrbTzLaN+dtbzczN7NQ6M5MjxTy9O/dr1wzlMHB4MM12p9fFN8kcx9lAGmNtJhntpkqlqyp2e9pYV11uUzlXaWKxQTq5Ctju7u8b+dO9wM9WP58HfA3A3R/v7me5+1nAJ4D/c1ygqSo9Gxh0Y9VSp1LzUeBDwDUrZnYm8ELg23VnJmkpoRsKuqnYtN0NNZhm2hUbyLc7CtLqkppkXJDoqpqTSogaFSvQqErTuecArwK+bGZbq8d+G/gN4INVF9IPgYunTcjMtrr7JjM7A7gMuBu4rarqfMjdr1zt9VNDjbvfXA38Wen9wNuAv5w2DZks5tiakLoONl1QsGlXiGAD6XVJTZNi+OhCrGraLIGmTpWm7aprCdz974BJO8t/OeW1v77i903V/ztWmeZEc8VXM7sA+I673zHP6yUdIbuhdEbUcJrtTi+XrigId7ZKyl1SfdTWutDZTjLNzHsuMzueQUnonTWff3F1JcEt+33frLPrhdi3TtD4mtUp2LQrxDibIYWb+GIGGlVp+meevda/AB4P3GFm32RwMZzbzOyfjXuyu1/h7pvdffNaO2b+lkoRur6bd07Bpm05BRsI+61bwSa8NgOlKjRS18x7LHf/srufNjJyeQfwTHf/buuta2jNKSfHbkJtfarWKNgMp9fq5AAFm9WoahNOm8t53m2k7SqN5KHOKd3XAbcATzazHWb22u6bJTEo2ISnYBO2OwoUbrrU9rINEWjqUtdTHuqc/XThlL+f1Vprei6FM6FCneYNOiPq8DTbPSMK8joraijU2VFDOZwCnpNcg6KqNGXRyftSrFy6oQbTbH2SnVVsSumOGlLlprkulp+qNDIPhZrExB5bA+qGqqPPwQbK6o4aUriZXVfLLNT6V5WmPAo1MpaCzXQKNt3uPmKd8aJwM12Xy6jJeleVRooNNTmd+bRSCtWa0BRsRqfZ+iQVbOYwPHAr4BzW9fIIGWhUpSlTsaFGmgt9UT4Fm9Fptj7JrINN7OuU9D3chHj/sdexlEGhJlGpVGtKudpw1xRsut+VpHDQ62P1JsR7bbpuu6rSqOspPwo1MpXG19SjYNPtmVGQRtVmqORwEzK8hQ40UrYit4acx9OMSqVaE1rOwaYLOQUb6E/VZqik6k3o9xEj0KhKU7YiQ01JUgk2Gl9TT1f3iFKwOVpKVZuh0YCTS8iJ1d7U1p2UQaFGalOwqUfBZiBUt0DKB8cUQ04KbWpjnalKI+NMvU1CbkrpehqVwu0ThkLeRiGEY3YfZN+G9g+KXdxOYTDdbm6pAO3fVgG6v7XC0PAgGfI2C/MYFyK6vE1DKkFqVOqBRvJWXKiRsoS4P5SCzUAX94uCwwegUOEm9WCzUt3gMRp+UgwrdaRcVZMyFBVqSqzSDPW5WqNgM266eQUbUNWmqVyDzFBbgabrKk3uXU/7fZFv7H9MB1P+VgfTbJ/G1MhcShtf06WcxthAGeNsIM2BxH0VM9BIvxSzhZRcpRlK5UyoodKCTW7XsBlMt5PJFhNsQOEmttjLvk9VGikk1PQh0AylFmxCU7AZN91OJtt5sIkRbiSctsNkKoODu/ocSzuKCDUST4m3UVCwOazrM0ZUtSlT28s4xHaiKk0Zsg81farSDKVWrSmtGwoUbEaVFmxA4aYrXSzXebcPVWn6KftQI2lQsJmNgs2RYnRHgcJNm7pYjqG2CVVpypF1qOljlWYotWpNDAo2k6bbyWSDXLws1tktCjfNpBZoVKXpr2xDTZ8DzVBqwSbG+BoFm0nT7WSyLO0tsztqSOFmNikur1m3T1VpypJdqFlzyskKNAkrNdjkqMuddandUUMpHqxT0vXy0fVoZF5ZbTkKM0dLrVoDZQabHKs1g2l3Numiu6OGFG6O1vXyCNntVPfzoa6nfGQRalSdWZ2CTRgKNkfrQ7CBw+GmzwEnxPtPYV1L3pLeghRmZBYaX7PatDubdLBgk8oBr2/hJtT7bbp+VaURSDTUKMzMTtWaap4KNqtMu7NJBwk2kNY3+dKrNyHfW+hAI+VKZg8xDDIKM/NTsKnmqWCzyrQ7m3TQYJNSuIFyAk6M9xFjXapKU66oewYFmX5QsJldzsGmj1WbUbkFnJjtbWMdqkojowLfeGVRQaZjKVZrQMFmHrkGG+h31WbUaGBIKeSk0KZYgUbXpSlbunsDmZuCzcg8FWxWmXZnkwbCfoNOOdiMWhlyQoSKGPOcJof1pa6nPC3GboBICY7ZfZB9G7o5WKzd4+xfbx1NG/av72TSwCDYLK/rbvqjhgfKAyceCjPDltQJGQfXjQ/OKQSUWbUVaFSlkXHSj8syF1VrRuYZ6IrDqtiMF3rMQw5VgFmNq7Yo0HRHVZrZmNmZZvYFM9tuZneZ2SXV4//ZzL5jZlurfy+uHn+hmX3JzL5c/X/eKtN+o5n9YzXd90xriyo1BbNde/CNHX4Nn9PSrkdZ3nhc2HnuXmZ5w1LQebYt94oNqGrTZ7EDjao0nToAvMXdbzOz9cCXzOym6m/vd/f3rnj+/cC/dfd7zexs4Abg9JUTNbPnAS8FznH3fWZ22rSGlPeVRo6gis3IPDMfXwN5V2wgTtWmxMpNbmKvAwWabrn7fe5+W/XzHmA7Y0LKyPNvd/d7q1/vAo41s2PGPPU/AO92933V63ZOa4s+7RKNgs18FGxmF/ug2mdtLvsQ2466npoxs7OAZwC3Vg+9wczuNLOrzWzjmJf8EnD7MLis8CTgZ8zsVjP7WzN71rT5q/upB1LthoJyu6K6HDgMeXdFQdgBxEPqkgqr7SCpbqd6fnhoia/ve2wXkz7VzLaM/H6Fu18x+gQzWwd8EniTuz9kZn8M/B7g1f//DXjNyPN/EvhD4EUT5rkIbATOBZ4F/A8ze4K7T0yeCjU9kXKwiUHBZtq0B/+XNM5mSOGme7lWxlSlWdX97r550h/NbIlBoLnW3T8F4O7fG/n7nwCfHvn9DOB64NXu/r8nTHYH8KkqxPy9mR0CTgW+P6kdeW55UpRYd/RWV1Sd6Xc6eSDeFWE13qYbXSxTVWnSZmYGXAVsd/f3jTz+uJGnvRzYVj1+EvAZ4FJ3/1+rTPovgPOq1zwJWMtgkPFE+kT3SKqDhiFesAlBwWa6mJe6V7BpR1chMdS2Ufdz1PXnOVPPAV4FnLfi9O33VKdt3wk8D3hz9fw3AD8B/M7I808DMLMrzWxYEboaeIKZbQM+Dly0WtcTqPupd1Luhip1fA3k3RU1mH6Z42yG1CXVTFfBsEmgUZUmHHf/O2DcDuizE57/+8DvT/jb60Z+3g/82ixt0VeUHlLFZsU8C7g4H5RTsYldtVHlZjYlBBpVacqhT29PKdismKeCTc3pdzr5H4l952WFm+m0jCRF2iIlSQo281Owac/wwK2D92EhlkeKVRrJgz6pPZZytQYUbJooKdikEG5AlYlQ7z+V9b2Sup7y0N9PqADpB5sYFGzqTr/TyR8hpQNd38JNyPfbdD2rSiP9+WTKRCkHm5KvYQMKNrNIKdhA+V1Tod9baut3lKo0+Sjz0ygzU7AZM18Fm5rT72d31KiSAk6M99HGOlWVRkChRjJRerDpWogdeJ+rNqNyDDgx2xwj0Ei58vnUSedSrtZA2cEmRHm7xGCTcriBI8NCaiEnhXalvv5AXU+5mbo1V7cL31ldpnj42OVmdnd1O/Hrq/s4SAEUbOJRsJlPDgfGoZghJ7WA1dZ6m2d7U9dTueps2R8Fzl/x2E3A2e5+DvBV4NKW2yURKdiMmWch42ug3GCTU7gZWhk02godbU+vbTEDzSxUpcnP1Hs/ufvNZnbWisduHPn1i8Avt9ssiS3le0SB7hPVVNf3ihrMo/v7Ra0U8/5RbUstiLQldvhUlaZsbXxqXgN8roXpSGJUsRkzT1VsZpyHqjZyWJvrRYODZZxGocbMLgMOANeu8pyLzWyLmW3Zf6jc8RDSHwo288wnyGyOoGCTlhQCzSzbu7qe8jR3qDGzi4CXAK9094lbirtf4e6b3X3z2jVhuwukOVVrJsxXwWaO+QSZzRFUtUmD1oGEMleoMbPzgbcDF7j7I+02SVKjYDNhvgo2c8wnyGyOonATT9vLXVUaWU2dU7qvA24BnmxmO8zstcCHgPXATWa21cw+3HE7JTIFmwnzVbCZYz5BZjOWgk1YOQYayVuds58uHPPwVR20RRKnM6ImzFdnRc0xn8H/oc+OgsMH2lLOkkqVAqTEMDXUiIxSsJkwXwWbOecVJ9iAwk1Xugozoao0uXc97fdFvvXoKbGbEU2ZF0KQTqkrasJ81RU157yCzWosjbdpT2qBRvpHoUbmomAzYb6FBZs+jLMZUrBpJsVA07cqjSjUSAMKNhPmW1CwgbADiGOHG1VtZqdlJilRqJFGFGziKi3YDOYVbFYT6UBdT9fLSFUamZVCjRSv5NspgIJNl4bhRgHnSCGWSSrbwEohP9syO4UaaSz1ak0sCjZN55XWgU3hJtwyaLreVaXpL4UaaUXqwab08TVQZrAZzC/o7Kbqa7gJ9Z5TW9+jVKVJn0KNtEbBZsJ8FWxamF/Q2dXSl3CT2/vU1YP7TaFGWqVgM2G+CjYtzC/o7GorddxNjPcUYx2r66ksCjXSOgWbCfNVsGlhfumGG8g/4MRsfxvrtcvtUV1PeVCokU4o2EyYr4JNS/MMPsuZ5RJwUmhnrECjKk15FGqkMwo2E+arYNPSPPMIN3BkcEgl5KTSllzWoeRBN7SUTukGmBPmG+gGmBDmJphwONiEuhnm4fnGuynmvMaFia5urJlCcJmkrUDTdZVGXU/5UKiRzinYTJhvgcEGwt7l+/A8B//nFm5G1Qkfw+CTclCpSxUa6ULY7qfFMDtVSY+6oibMt8CuKIh3Wm3pB8pUuoyaanM9qUojo4KPqfGN65P+1i7dUbCZMF8Fm5bnG2W2UlPsQCNlizZQWMGmnxRsJsxXwabl+SrcpEjrRLoW9ewnBZt+UrCZMF8Fmw7mrQNpKtpeD/NuV+p6Klv0U7rVHdVPCjYT5qtg09H8o86+91IJNNINMzvTzL5gZtvN7C4zu6R6/HIzu9vM7jSz683spOrxU6rn7zWzD60y3V+ppnfIzDbXaUv0UDOkYNM/CjYT5ltwsFHVpn9SWua62F5nDgBvcfenAucCrzezpwE3AWe7+znAV4FLq+f/EPgd4K1TprsN+EXg5roNSSbUgKo2faRgM2G+hQYbiP8tO6WDbOm6WNahth91PdXn7ve5+23Vz3uA7cDp7n6jux+onvZF4IzqOQ+7+98xCDerTXe7u//jLG1JKtQMKdj0S+rBJhYFmy7nr3DTtdQCjao0jZ1qZltG/l087klmdhbwDODWFX96DfC5bpuY8MX3fON6Hex6JOUL9MW6OB+Ue4E+iHORvqPbMPg/54v2pSj3wJhzlWb/oQXueeSkLiZ9v7uvOq7FzNYBnwTe5O4PjTx+GYMuqmu7aNioJCs1Q6ke5KQbKYfYWN1QoIpNCKrctKPL5agqTdrMbIlBoLnW3T818vhFwEuAV7p75x/4pEMNaJxN3yjYTJh34cEmpXAj8+ly2aWyfch4ZmbAVcB2d3/fyOPnA28HLnD3R0K0JflQM6Rg0x8KNhPmXXCwgXQOXKrazC7l5TXrtpxz11NEzwFeBZxnZlurfy8GPgSsB26qHvvw8AVm9k3gfcCvm9mO6mwpzOzK4enbZvZyM9sB/DTwGTO7YVpDkh1TM47G2fSHxtjEF3qMDaQxzmZI422mCxFmUgm7Mll1JtO4D+5nV3nNWRMef93Iz9cD18/SlmwqNUOpHuikfSkH2D6c6g39rtgMqXIzXg6BRlWa/sku1ICCTZ8o2IyZb0+CjcJNmkIth9TWv+Qhy1ADCjZ9omAzZr4Rgo2qNgPDg3rfAk5u71lnPPVTtqEGFGz6RMFmzHwjlMoVbI6U24F+HjHeY4x1rq6nMmQdakDBpk8UbMbMt0fBJodwU1rAifF+2ljPXW6jMc+AlOmyDzWgYNMnCjZj5tuTYANpV22Gcg83MQNarECjKk05igg1oGDTJwo2Y+arYJOc3Ko3ObVVZJJiQg0o2PSJgs2Y+SrYJGs04KQUHFJqU+rdTqCupxwUFWpAwaZPFGzGzLdnwSa3cDMUK+SkG67irUd1PZUlqysK16UrD/eHrjw8Zr4B7+w9FOPqw0MpXYV4XqsFjHmuaJxSYJmmrUCjU7gFCg01oGDTJwo2Y+bbw2ADZB9uxskpoMwqdqVtliqNup7yUFz306hUD3TSvpQDrLqiwol9kJT62lxXsbc7SUfRoQYUbPpEwWbMfBVsJEEprCONpSlT8aFG+kXBZsx8IwWbmOEm50HEpWt7vYTYztT1lI9ehBpVa/pFwWbMfCN9K1XVRkblGGgkL70INaBg0zcKNmPm2+Ngo3ATX0rrQF1P5epNqAEFm75JOdjE0tdgA2kdVKW5UNuUup7y0qtQAwo2fZNqsIm5o+x7sFG4CS+lZa4qTdmKvU6NyFCq17GJdQ0biHMdG4h7LZtRJVywLxddBJoUAnKqDhxcw86H18VuRjS9q9SAqjV9pIrNmHn3uGIDqtp0ravlG3L7UddTfnoZakDBpo8UbMbMu+fBBhRuupDq8lTXU/mmhhozu9rMdprZtpHHTjazm8zsa9X/G7ttZjcUbPpHwWbMvCMGG4Wb8nS5DFPaXiRNdSo1HwXOX/HYO4DPu/sTgc9Xv0/lC+rDlvgUbMbMO+I32NQOVAo280t52c26javrKU9TQ4273ww8uOLhlwJ/Wv38p8DL6s4w1sDISVSt6ScFmzHzVrD5EVVtZhNieaW2jUia5h1T81h3vw+g+v+09poUnoJNPynYjJm3gs0RFG6mC7F8mm4bGkvTH50PFDazi81si5lt2b/8MJBetUb6y3btSTLcKNikReFmvFKXibqe8jVvqPmemT0OoPp/56QnuvsV7r7Z3TevXTrhR4+nFmxUrek3BZsV81awGUvhZiDkclCVRmYxb6j5K+Ci6ueLgL+cZyIKNpISBZsV844cbBRu0hP6fae8DUia6pzSfR1wC/BkM9thZq8F3g280My+Bryw+n0uCjaSEgWbFfOO/C039YNan8JN6PfZxrqfZ/tV11Pept4mwd0vnPCn57fViOWNx2lDkmSkeFuFPt5SYSiVWyusZvSAX9rtF/oS2qQMyVxROKWKTWoHNAlPFZsV81bFprZSqjcx34eqNDKvZEINKNhIWhRsVsx797LG2cxgGApyCjgptDmndSzpSSrUQFrBRkTBZsz8VbWZWQphYTUpt20esbdRiSe5UAPpBBtVawQUbMbOX8FmbqMBJ2aQSKENK8Vcr7E/U9KOqQOFRUSDh8fOXwOIWzEuVLQ92Dil4DJJW4EmduCWuJINNamcEeUb1yf5TV3CU7AZM/8Egg1QRLgZlUMIaVPsylsKxxppR5LdT0PqhpLUpBhwY++QU/hmHPugKGkIsS2muA9IgZldbWY7zWzbyGObzOyLZra1ul3Ss6vHN5rZ9WZ2p5n9vZmdPWGaHzWzf6pev9XMNk1rR9KhBtIJNiJDKe7UFGwUbHKl9VaMjwLnr3jsPcDvuvsm4J3V7wC/DWx193OAVwMfXGW6/9HdN1X/tk5rRPKhBtIINqrWyCgFmzHzTyTY6CCZjxTWVezPTSnc/WbgwZUPAydWP28A7q1+fhrw+ep1dwNnmdlj22hHFqFGJEUKNmPmn0CwgTQOlrK6tteRup6S9CbgcjO7B3gvcGn1+B3ALwJUXVI/DpwxYRp/UHVTvd/Mjpk2w2xCjao1kqIUd3IKNgOq2qQrx0AjnFqNixn+u7jGa/4D8GZ3PxN4M3BV9fi7gY1mthV4I3A7cGDM6y8FngI8CzgZePu0GSZ79tM4qZwRJTJKZ0WNmX/ks6JGlXLqdylSCpolHk8OHlrD7j2dfPbvd/fNM77mIuCS6uc/B64EcPeHgH8PYGYG/FP17wjufl/14z4z+wjw1mkzzKZSMxS7YpPawUvSoIrNmPkn9O05pQNpn3WxHlLazuQo9wI/W/18HvA1ADM7yczWVo+/Dri5CjpHMLPHVf8b8DJg28rnrJRVpSYVunaNjKOKzZj5J1axgfKuaZOL1ILlrKFf+/zVmdl1wM8x6KbaAbwL+A3gg2a2CPwQGHZZPRW4xswOAl8BXjsync8Cr3P3e4FrzewxgAFbgd+c1o4sQ426oSRVCjZj5p9QsAF1R5VEVZp0uPuFE/70L8c89xbgiROm8+KRn8+btR3ZdT8NqRtKUpXiN7rYXwJSO/hoEHFYWtYSSrahBuIHG5FJFGzGzD+xYAM62IbQ1TJusj2p66lcWYea2FStkdWkuCNMIdikFm5UtelGl8s1tW1I0pF9qFG1RlKmYDOhDQkelBRu2qPlKLFkH2ogbrBRtUamUbCZ0IYEgw3ogNxU18uv6XajrqeyFRFqRFKnHeN4KQcbhZvZaZlJbMWEGlVrJHWpBZsUqjWQbrABHaRnEWJZha7SSH6ChhpfsE6nr/E1kjoFm/FSDzYKN6srdfmk9nmV6Yqp1MSmao3UldqOUsGmHoWbo4VcJqlvH5KG4KGm6yuLqlojOVCwGS+HA5fCzUBuyyCVbVy6FaVSU2qwUbVGZqFgM14OwQb6G25ivO9ctgmJL1r3U0r3ghGJRcFmvBQv0jdJn8JNjPfZxnYwz3ad2mdT6ok6pqbLYKNqjeQitZ1nKsEG8vqGXnK4Kfm9SWmo2rkAACAASURBVFmKHiis8TWSCwWbyXIKNnA4AJQQAmK/j1hVGslX9FBTYjeUqjUyDwWbyXILNkOxQ8G8Umh3zHWe2mdR6oseaqDMbiiReaS2M1WwaUcu1Zsc2iiymiRCTddiBBtVa2ReCjaT5RxshlILOKm1B9pbzyltuxLGYuwGDC1vWCpihyWHHXrgwVrPW3PKyR23RJpa2vVoMlXP4X6ihK7rlUFi34aFKPNNSezjQGpfKmQ2yYQa6DbYLG88Lnhq943re/sBqRtoVj5XAWfAdu1JrtqXUrCBwcGvhGAzalLYmDfspBxeuqYqTT8lFWpAFZvczRJmpr2+7wFHwWa6EoPNOH0JJ9r3N+cHjeW9a2M3I5pejKkZ0tiabjUNNOOmN/zXVylW+lL7BqwDYRnaXI/zbqMpft5kNkmGGp0NlZcQwaPP4SbFHW2KwUbhRkSSDDVQxiBAaV9fw42CTT0KNnnSepO2JBtquhS6WlNyF1SsgNHHcKNgU48OkHlpe32p66nfkg41qtakLYVQ0bdwk+KOV8FG5qX1JG1LOtRAd8FG1ZqyKNjElWqw0UEzXV2smxS3Qwkr+VDTJQ0aLkufqjYKNvUp2Mg0KX6eZD5ZhJpSuqFKqtakHB5SblubUtwRK9hIHarSSFeyCDVQTjeUhNGXqo2CTX3qjkqD1oF0KZtQU4qSqjU56EO4UbCZjQ6q8XS17Jtsb7N+fkrfn+Quq1Cjao3Mq/QdkYLNbBRswtMylxCyCjVQzvianOUaEHJtd10KNrNRd1QZUt7GJLzsQk1XQlZr1AUVT+ndUQo2s1Ow6V6qyzjFz4s0k2WoUbVGmlKwCSuHYJPqgTd3XS7X0NtVyfuNUjQKNWb2ZjO7y8y2mdl1Znbsas8/tGBNZneELoKNqjX9UnLVRsFmPgo27dLylNDmDjVmdjrwW8Bmdz8bWABeMe11+zYszDvLIDRouH8UbMJRsOmPrpdj020pxc+HNNe0+2kROM7MFoHjgXubN6m+3LuhVK1Jh4JNOEu7Hk0+3Kg7qhktO4ll7lDj7t8B3gt8G7gP2O3uN9Z5rao1kqJSu6NSDDaQT9VGB+j6Qi2vGNtOifuGNpnZ1Wa208y2jTz2n83sO2a2tfr34urxV448ttXMDpnZpjHTfLqZ3WJmXzaz/2lmJ05rR5Pup43AS4HHAz8GnGBmvzbmeReb2RYz23Jg38M/erytYJN7tUbSU+LOS8GmGQWb6XJaRql+HjL3UeD8MY+/3903Vf8+C+Du1w4fA14FfNPdt4557ZXAO9z9p4Drgf84rRFNup9eAPyTu3/f3ZeBTwH/auWT3P0Kd9/s7psXjzmhwewmy3nQsLqg0qRgE05OwSanA3dIIZdLLttL37j7zcA8O84Lgesm/O3JwM3VzzcBvzRtYk1CzbeBc83seDMz4PnA9lkmkHo3lPSbgk04OR2oFG6OpGUhU7zBzO6suqc2jvn7v2NyqNkGXFD9/CvAmdNm1mRMza3AJ4DbgC9X07pi1umk3A2lao2UOs4mRTkFG9DBHMIvgza2kXmCfVb7gEPGwt6F1v8Bpw6HklT/Lq7Rmj8G/gWwicHY2/82+kcz+z+AR9x925jXArwGeL2ZfQlYD+yfNsPFGo2ayN3fBbyryTREcnDogQdZc8rJsZvRCtu1J9kgvbTr0awG6g8P6n0c26dQ1zv3u/vmWV7g7t8b/mxmfwJ8esVTXsHkKg3ufjfwour1TwJ+Ydo8k7iisKo1koOsvq1NkWo3FORXsYF+dUn16b1KM2b2uJFfX86gO2n4tzUMupQ+vsrrTxt57n8CPjxtnkmEGkg72ISQ6jdnOZKCTRg5Bhso+4Af+73luk30hZldB9wCPNnMdpjZa4H3VKdj3wk8D3jzyEueC+xw92+smM6VZjasCF1oZl8F7mZwHbyPTGtHo+6nvljeeJw+UPIjw2BTQneUuqK6UVq3VOyg1tb+t/jxNBG5+4VjHr5qlef/DXDumMdfN/LzB4EPztKOZCo1oGpNqgcXGa+UnZ0qNt2JXd1oKvf2S/8kFWog3dO8c/3GKN1SsOle7sEGDoeDXAJCSm2NWaWR/CQXatqiao2EomDTvRKCzVDKASfVdonUlWSoUbVGcqNg072Sgs1QCiEi6ZAVeZ2X8rnukyRDTVtyrdZInkrZASrYhDcaLILdFDLRINOFlLdpaVeyoSbVQcMhqjXqgsqXgk33Sg02K7UVcsYFphzCTF/Ws7Qr6VO6921Y4JjdB2M3Q2QmpZzyrdO905JDEGmLAo3MK9lKTZtUrZEYSqjapF6x0cFPppl3Gy7h89tHyYeaVAcN91nuFYiQStgxphxsQN/qS6P1KU0kH2qgnWCTY7VGyqBg0z0dCGWc1LdbaV8Woaav1AVVDgWb7inY5E/rUJrKJtSoWiO5KyHYpE4HxXyltO70Wc1X0FDjGh4jPZf7zjL1ag2kdXCUeHLYVqV92VRqoJ/VGnVBlUfBpnsKNnnR+pK2BA81+9dbo9enGGxEZqVg0z0dKEX6J0qlpmmwSY2qNTIPBZvu6Vo26eti/TTZNnP/XPZdVt1PQ6rWxKdr1bQj9x1oDsEGVLVJldaLtC1aqFG1RmRAwSYMHUDTovUhXYh676f96421e3yu17ZxX6jlDUvZ3E/FN67P5uAhszv0wINZV79Svk/UqD7eM6pv+r6ftIOw+FCWnTCt6O87F0mMKjZhqEIQX6rrIPfPoCQQapp0Q6U2tkbfAKWp3HeqOQWbVA+spetyueey/Ul3ooeapvp0w8scyvvSnIJNOAo2YWl5S9eSCDWxBw2rWjOfnMeApE7BJhwdaAXy/8zJQBKhBuJ3Q+VC1Zr+yH0nq2Ajo7pexjltb9KdZEJNbKrWSIpyDzY50Tib7mi5SihJhRpVa0SOlnOwyfHbsw7AIvlKKtTElku1Rl1Q/aNgE5aCTXtCLMum21jOny85UnKhRtWavGiwcDg573gVbPpJy1BCSy7UQNyzoVStkZQp2ISlcTbzC7XcctyupDtJhpomVK2R0inYhKdgMxstL4kl2VATsxtKd/CW1CnYhKeqTT0hl1Eb21LOnyU5WrKhphTqgpKu5LwzzjXYgKoQq9GykdiSDjWq1uRBg4XjUbCJQ1Wbo2l5SAqSDjUQ/xYKbdDF+KRLCjbx6EA+EGM55L7tSDeChhoPPIa3D9UadUEJKNjE1PeqTc7vPefPjYyXfKUGVK0RqSPnHXTuwQbyPrjPK9Z7jrm93PDwNdHmLdMFDzX7AxcWVK0RyUMpwaYv4aYv71PykkWlBlStSZ0GC6ch52pNSUoONyW/N5mfmV1tZjvNbNvIY5eb2d1mdqeZXW9mJ4387Rwzu8XM7jKzL5vZsWOm+XvVa7ea2Y1m9mPT2hEl1KhaI9KdnINNCdWaUaUd/FN4P21tIzl/ThL1UeD8FY/dBJzt7ucAXwUuBTCzReBjwG+6+08CPwcsj5nm5e5+jrtvAj4NvHNaI6JVauYJNiVUa7qiLigZlfMOu8Rgk0IYaKKE9yDdcvebgQdXPHajux+ofv0icEb184uAO939jup5D7j7wTHTfGjk1xMAn9aObLqfmkqlWlNyF5SkRcEmLcNgkFs4SKm9sbcLDRJu5DXA56qfnwS4md1gZreZ2dsmvcjM/sDM7gFeScqVGlC1pjQaV5MeBZs05RJucmijHMkOwdLe9v8Bp5rZlpF/F9duk9llwAHg2uqhReBfMwgq/xp4uZk9f9xr3f0ydz+zeu0bps0ry0rNvMGm9GqNuqBkHAWbdKUYbnKtKM0i589ERPe7++aRf1fUeZGZXQS8BHiluw+7j3YAf+vu97v7I8BngWdOmdR/B35p2vyih5rQg4ZF+ijnnXjpwQbSCBKx5z9NH7aD0pjZ+cDbgQuq8DJ0A3COmR1fDRr+WeArY17/xJFfLwDunjbP6KEGwnZDlV6tEZlEwSYPIQNOCmFKymBm1wG3AE82sx1m9lrgQ8B64KbqtOwPA7j7LuB9wD8AW4Hb3P0z1XSuNLPN1WTfbWbbzOxOBoOLL5nWjsW235jE5RvX9+oAIP1hu/b0rot1NGy09UUoxwCjfVr63P3CMQ9ftcrzP8bgtO6Vj79u5Oep3U0rNQo11YV0rgTOZnCq1Wvc/ZZ5prV/Paydcbvdv95Yu2fqGV5H2bdhgWN2H3X2WG3LG5ZY2j3ulHpZc8rJWVcESnfogQezHtDdx2AzlGMYSdG8+yed+ZSHpt1PHwT+2t2fAjwd2N68Sf2hAcMSQ+6hU9/a+0XrW2Yxd6gxsxOB51KVl9x9v7v/oEljchpbI5IzBRsRKVGTSs0TgO8DHzGz26vBPSe01K7kacCw5E7BRkRK0yTULDI4r/yP3f0ZwMPAO1Y+ycwuHl6s5+AjD0+dqKo17YjZBZXzmI2+UbCRlLW9fnPf3mW6JqFmB7DD3W+tfv8EYy6e4+5XDC/Ws3B8vUJOLteuUbVGSqAdvcjqNEg4H3OHGnf/LnCPmT25euj5jLl4Tiiq1hxNA4alD1StKZPWq8yj6dlPbwSurS6Mswn4r82bNJDLfaHaqtaIxJR7tUYHQBGBhqHG3bdWXUvnuPvLqqsEZieFak1pXVAaV5MfBRtJRRfrMvftW+oJepsEn3FuqtY0py4omUXuO34FG5F+S+LeTylQtUZkQMFGYtL6kyaCh5rldbM9X9UakfAUbKQkTbZnnfmUF1VqRpRarYnVBaVxNRKTgk1+tM6kqSihRtUakfTlXq0BHSRF+iZapWbWYBNKCtUakVQo2EgoWk/Shmy6n3Kp1rShpC4oyZ+CjeSshO1X6osaakJ0Q82jSbVGXVBH0riaMpRwYFCwSZfWjbQlm0rNvFStOUzVGuk7HTxlFjrzKT+LsRuwvA6W9tZ//v71sDbAfmnfhgWO2X1wrtcub1hiafdyyy0SievQAw8WUXmzXXsU8BOioNkuOxjmGJmq4is1kG+1RiQ1JXRDgQ6kfVHK9ir1JRFqNLZmwjQK6YIq4du9HFbKgULBJj6tA2lbEqEmBFVrRNqjYCMiKUom1KhaM2Eauh+USKcUbOLQcpcuJBNqoPsL8qlaM6AuKGlDKdUa0AG2RCVtn1JfUqFmVjlUa0RKVtKBQ8EmnByWtU7nzlNyoaa0ao26oKR0CjYyCy1j6VJyoWZWqtbMR11QIuPpoCuSryRDjao1Y6ahao0krKRqDSjYdCXUci1te5T6kgw1s1K1Zj66qqq0qbQDiYKNSH6ChhqfIROoWjNmGqrWSOIUbGQSLUsJoYhKDYSr1kgzGlcjudHBuH905lO+goeaAyceqv3cFKs1sS/G1zZ1QUnbSqvWgIJNUyGXX4nbn9QXpVIzS7CZRR+qNeqCkhyUeGBRsBFJX/LdT6VVa0RdUH2hYCOgZSZhRQs1XXVDpV6tSXHAsLqgROrTQbo+LSsJLflKTQiq1oh0o8RqDehgLZKqqKFG1Zp0qFojXVGw6acYy6fUbU3qU6Wmklu1JvcBwxpX0y+lHmwUbMbTcukfM7vEzLaZ2V1m9qYVf3urmbmZnTrmdc8zs60j/35oZi+btx3RQ42qNSKSMx3Aj5T78tA1amZnZmcDvwE8G3g68BIze2L1tzOBFwLfHvdad/+Cu29y903AecAjwI3ztiV6qMldSdUadUFJl0qt1sDgQJ77wVykgacCX3T3R9z9APC3wMurv70feBvgNabzy8Dn3P2ReRuSRKhJpVoT+tYJfacuqP4pOdhA/lWKpmK+/9K3rQScamZbRv5dPPK3bcBzzewUMzseeDFwppldAHzH3e+oOY9XANc1aeRikxdLM8sblljavRy7GSLSItu1p5dVz74HulTYQVi7p05RZGb3u/vmcX9w9+1m9ofATcBe4A7gAHAZ8KI6EzezxwE/BdzQpJFJVGog72qNuqBE6uvDN+q+HeD79n7laO5+lbs/092fCzwIfBN4PHCHmX0TOAO4zcz+2YRJ/Cpwvbs3+qafTKiB7m6fkDINGJY+UrARKYuZnVb9/8+BXwSucffT3P0sdz8L2AE8092/O2ESF9Kw6wkSCzWzULXmsJxP79a4GilZHwYQp/D++hCSM/BJM/sK8D+B17v7rklPNLPNZnblyO9nAWcyGGDcSHJjag6ceIjFh7LNWnNJbWyNb1yfxI5KynbogQd7E2pLHWdT2n5Cp3PPz91/Zsrfzxr5eQvwupHfvwmc3kY7sk4PqtaI5K1P37BLCwClvR8pQ9hQs6beiOw+jq1pKucBw335ti5SSndUCe9BypR1pQbKqdZowLD0VZ+qNUM5h4LU2t7H7UcmCx5qDq47WOt5qtbMLucBw9JvfTwwpRYO6sixzdIv2VdqQNWarqgLSqRbOXVH5dJO6bcooUbVGhFZqY/VmqHUA0Pq7RMZKqJSA+lVa+bVtFqT84Bhkb4Hm9TCQ4ptGtXn7UXGixZqSq/W6PTu2akLSkAHqlSCRAptCEnXqClDMZUa6LZakxMNGBbJX6xwk0qoEplH1FCTU7VGA4ZFwul7tWZUqJChMCMlKKpSA6rW5E5dUCLjDUNH28Ej1zCj4CvjRA81qtaMl9qAYZHQdNCarI2Ak2uYEVlN4xtamtkCsAX4jru/ZJ5pHFx3kIW97Q2sXV4HS3vrPXf/elirz/VUusmlxNCnm17OS59LkcPaqNRcAmxvYTpT5VqtmVdfqzU6iImIyDwahRozOwP4BeDKpg2p2w1VV0pja0o5vVsDhiUGdUOJSF1NKzUfAN4GBCuh9K1aIyIKNnKktrcHXaOmHHOHGjN7CbDT3b805XkXm9kWM9tycO/Dq05T1ZqjqQtKRESkniaVmucAF5jZN4GPA+eZ2cdWPsndr3D3ze6+eWHdCQ1md5iqNfGoC0piUbVGRKaZ++wnd78UuBTAzH4OeKu7/1rTBsU8E6pr+zYscMzu2atRyxuWWNq9PPd8lzcex9KuR+d+vYiI5GHNQZ/rOFOKoNepsQVvbVpdVGt0Mb60qAtKVlK1RkRW00qocfe/mfcaNePEHFszi5AX40uJuqAkJgWbftP6l9UEv6Lw0rr9rU2rT9Wavg4YFhERqSv6bRImKblaU4JQ1Rp1Qck4+rYuIuNECTV9r9bo9G4REZH2RavU1Ak2qtaIyCSq1ojISsl2P81C1Zo41AUlsSnY9IvWt0wTNdTEqNZ0JZdqjbqgREQO0y0SylJEpQbqV2tKuHVCH6laI5Po27uIDEUPNarWzC6lAcO6Zo2IiKQieqhpU47VGhFpTtUaEYFEQk1J1ZpZacBwfeqCEukvBVepI4lQ06YuqjWz0IBhkTh00BORZEJNytWaVAcMp1StCUXVGhERmSSZUNMmVWvq0YBhKY2qNSL9llSoUbVGRJpSsBHpr6RCTZv6Uq3p4+nd6oISEZFxgoaahTXTg0afqzUiInI0Vd+krmIrNRC/WjMrDRgWaYcOgiJhmdlJZvYJM7vbzLab2U+b2e+Z2Z1mttXMbjSzHxvzuk1mdouZ3VU99981aUfwULNh/aNTn1NKtUYDhrujLigRkaR8EPhrd38K8HRgO3C5u5/j7puATwPvHPO6R4BXu/tPAucDHzCzk+ZtRNGVGuhPtUZEjqRqjUgYZnYi8FzgKgB33+/uP3D3h0aedgLgK1/r7l91969VP98L7AQeM29booQaVWvapwHDIkdTsBFpzalmtmXk38Ujf3sC8H3gI2Z2u5ldaWYnAJjZH5jZPcArGV+p+REzezawFvjf8zZycd4XNrVh/aPs3hPmqrYHTjzE4kP5FKX2bVjgmN1l3hZCRES6Ywedpd3LXUz6fnffPOFvi8AzgTe6+61m9kHgHcDvuPtlwGVmdinwBuBdY9tt9jjgz4CL3L1eF8sYSR/pQ1dr+n56t0iJVK3JW5fr74aHr+ls2j2zA9jh7rdWv3+CQcgZ9d+BXxr34qr76jPAf3L3LzZpSNRQU6cbqi11x9bUVeLp3eqCEhGRWbn7d4F7zOzJ1UPPB75iZk8cedoFwN0rX2tma4HrgWvc/c+btiXpSg30t1qjAcMiIpKRNwLXmtmdwCbgvwLvNrNt1WMvAi4BMLPNZnZl9bpfZTDI+NerU7+3mtmmeRsRbUzNUMixNW3bvx7W7ondiiMtb1hq1J+6vPE4lnaFq6CJhHDogQdV1RPpkLtvBVaOuRnb3eTuW4DXVT9/DPhYW+1IvlID7VVrdHp3WOqCEhGRkJIINSHH1rStxNO7RUqkAcMi5Usi1NShak04GjAsIilQEJVZJRNqVK1pl6o1IkfTQVKkbMmEmjrqVGvqyK1aE0OO1RoREem3pEJNG9WaWLdOmIVO7+6GuqCkDlVrRMoVNNQsLjS/AF7oak1dKV6MT11QIiLSJ0lVaiC9ak0qXVAaMFyPqjUiIv0VPNScdsLextMooVqjAcMi8agLSqRMyVVqQNWalKhaIyIiuYgSatqo1rQlp2qNBgyLtEfVmrRp/cg8kqzUwPRqTSk3ugxBXVAiItIH0UJNStWatpVWrVEXlIiI5CDZSg2Eq9b04WJ8qtaIHE1dHCJliRpqSq7WdEnVmnpUrRER6ZekKzV1pFqt0endIiIiYUUPNdOqNTnf6LI0bVZrRFKhLiiRckQPNW0IXa2pq7QBw21SF5SIiLRtMXYDYFCt2fnw5H6dDesfZfeedKoEy+tgKdPhQMsblljavRy7GSIi0gE76Czt6m8PRxGVGlC1JhQNGJYSqQtKpAzJhJrcxtbo9G4REZG0BA01a9e0d4XfcVK90eUs+latCUXVGhGR8iVTqYEw1ZoYt06YpQuqb0J1QYlMoy6odGhdyLyCh5ozj/9Bp9Nvq1pTOnVBiYhIaeYONWZ2ppl9wcy2m9ldZnZJGw1KpVqT08X4cu+C0oBhERFpQ5NKzQHgLe7+VOBc4PVm9rQ6L1S1Jg2q1ogcSd0eInmbO9S4+33uflv18x5gO3B6G43KrVpTl6o1k6laIyKh/fwJr47dBGlZK2NqzOws4BnArXVfU1K1Rqd3i4iIxNc41JjZOuCTwJvc/aExf7/YzLaY2ZYf7vph7emqWtO/ak0oqtaIiJSpUagxsyUGgeZad//UuOe4+xXuvtndNx+78dgj/tZ1tSaknKs1qdDp3ZICjasRyVeTs58MuArY7u7va69JhzWt1qR664TU9LELStUaEZHyNKnUPAd4FXCemW2t/r141on0sVqjLqjJVK0REZF5zX2Xbnf/O2C2I+4cmt7Be2ndfpb3rm3cjgMnHmLxoaQuwNyqPt69e80pJ6urQUSkIOUepWeQ+q0TVK0RERGZLolQM60LKsTYGtHYGpEhVfBE8pREqEmBTu9uTtUaERGJKWioWWsHJv6tabVmGl2MTyZRtUZEZH5mdqyZ/b2Z3VHdC/J3q8ffYGZfNzM3s1NXef0/N7Mbq3tJfqW6oO9ciqnU5HoxvtT0sQtKREQa2Qec5+5PBzYB55vZucD/Al4AfGvK668BLq/uJflsYOe8DQkean78uAcm/q2P1Rp1QU0WsgtK1RoRkfn4wPAAvVT9c3e/3d2/udprqxthL7r7TdW09rr7I/O2pZhKDaha0xZVa0REZBZmtmBmWxlUWW5y97r3gnwS8AMz+5SZ3W5ml5vZ3N/Io4SaJtWaplSt6V6u1RqRUToDSuQIpw7v41j9u3j0j+5+0N03AWcAzzazs2tOdxH4GeCtwLOAJwC/Pm8j5774XixNL8ZXx8F1B1nYu3oYKP1ifH2ki/GJpCHkZ/HnT3g1Nzx8TZB5BXHgILZrTxdTvt/dN097krv/wMz+Bjgf2FZjujuA2939GwBm9hfAuQxuwzSzJI/KJd06oa7UqjUpdUGpWiMiki4ze4yZnVT9fByDwcF313z5PwAbzewx1e/nAV+Zty3RQs1qXVDT5HQxvr6e3t1mF1RIGjAsIjKzxwFfMLM7GYSUm9z902b2W2a2g0GX1J1mdiWAmW0e/uzuBxl0PX3ezL7M4PZLfzJvQ5Ltfjrz+B9wzyMnRZu/uqD6eT8oERGZjbvfCTxjzON/BPzRmMe3AK8b+f0m4Jw22hL1iKxqzZFS64JqKtcBw6rWiIjkKekyQ+yxNTq9O62xNSEp2IiI5Cd6qFG15kiq1kymAcMiIrKa6KFmGlVr4lO1RkREcpBEqFG15khdVmtiULVGRERCSCLUTJNDtSYXMbqgcqZqjUj5fv6EV8dugrQkaKg5ds3k04O7rNZM00a1Rl1Q4ahaI1I+faGQeWRRqYFm1Zo2bnTZlhS6oHIfMByadq4iInkIHmp+4pjvTfxbzGrNNBowrGqNiIikLZtKDXRbrenbgGFVa2ajao2ISPqihJquqjVdU7UmLarWiIjIqKwqNdP0/fTuEJp2QalaIyIiXYkWauat1uRwenesak0OXVAiIinSad1lKKpSA3lVa3LV5wHDqtaIhKPPm8wqaqgpuVpThwYMi4iItKe4Sg2kUa0pfcCwqjUiIpKa6KEm12pNW1StyZOCjYhIeqKHmq7oYnzd63O1RkRE0rMYuwEwqNZ8fd9jx/7tx497gG89esrYv515/A+455GT5prnhvWPsnvP5GrB0rr9LO9dO9e0Z7W8DpZqZLD962HtnnrT3L/eWLvHa7dh34YFjtkd9sadyxuPY2lXOrewmNWaU07m0AMPxm6GdETVuDSE/Jz9/Amv5oaHrwkyr84cPNDr/VKxlRpQtaZvYlRrdOATEUlH0FCz1g60Ps1Sbp2Qqz5fjE9ERNKSTKWm5Btd1tHXAcNtU7VGRKS/goeaJ6z9fuvTTL1aU3oXVEoDhkWkPPriIHUlU6mBfG902Za+Vmva7oJStUZE5qXbJeQtSqjpolqzmq4vxqcBw6rWgIJNSbQuRfKUVKUGdDG+FKo1MZRQrRERkbiihZrUqjVNlVStmbcLStUafcMX6Yo+W1JHcpUaiFOtyfH07lmqNTlQtUZE0iPe9gAABnVJREFURJqIGmpUrRmvbhfULHIYMFwKfaPMm9afSL6SrNSAqjW5Su1ifLGqNTowiuRLZ0DlK9lQ05UUqjV1pDBgWNUaEUmJvizINNFDzWpdUKVWa3IZMDwvVWsGtAMWEQkreqiRyVStyZ+CTV60vkTylkSombdas5rVqjW6GF/3VK0REZHQgoaaY81bnV4fbp1QV2mnd5dE3/5F2hPq86TBwnlKolID3VRrVpPCgGGd3r26kqo1Cjbp0zoSyV/wUPPExXbODhrKecBw6XSFYRERCalRqDGz883sH83s62b2jqaN6ePF+OrQgOH2qFojItK+aXnAzI4xs/+n+vutZnZWF+2YO9SY2QLwfwP/BngacKGZPa3Oa+ep1uj07jylNmA4NgWbNGm95EXjatJSMw+8Ftjl7j8BvB/4wy7a0qRS82zg6+7+DXffD3wceGnTBpVWrWmLqjXt0ZlQIiKtqpMHXgr8afXzJ4Dnm9lsB6EamoSa04F7Rn7fUT3WmRSrNdPo9G5Va1ZSVSAtWh8ijdXJAz96jrsfAHYDp7TdkMUGrx2XsI46Z9vMLgYurn7du/H0Hf84fdLfatCso5wK3N/mBBOl91kWvc9QHgkyl/jvM4zi3qfZn417+Mmh21HXQ4cevOHGR/7s1A4mfayZbRn5/Qp3v6L6uU4eqJUZmmoSanYAZ478fgZw78onVW/6ipWPh2JmW9x9c6z5h6L3WRa9z7LofZZlxcE9Ke5+foTZ1skDw+fsMLNFYAPwYNsNadL99A/AE83s8Wa2FngF8FftNEtEREQyUScP/BVwUfXzLwP/n7unU6lx9wNm9gbgBmABuNrd72qtZSIiIpK8SXnAzP4LsMXd/wq4CvgzM/s6gwrNK7poS5PuJ9z9s8BnW2pLV6J1fQWm91kWvc+y6H2WpS/vs7ZxecDd3zny8w+BX+m6HdZB9UdEREQkuGTu/SQiIiLSRK9CjZm91czczLo43S06M7vczO42szvN7HozOyl2m9rS9i05UmRmZ5rZF8xsu5ndZWaXxG5Tl8xswcxuN7NPx25LV8zsJDP7RPW53G5mPx27TV0wszdX2+w2M7vOzI6N3aa2mNnVZrbTzLaNPHaymd1kZl+r/t8Ys41yWG9CjZmdCbwQ+HbstnToJuBsdz8H+CpwaeT2tKLJLTkycwB4i7s/FTgXeH2h73PoEmB77EZ07IPAX7v7U4CnU+D7NbPTgd8CNrv72QwGinYyCDSSjwIrT5N+B/B5d38i8Pnqd0lAb0INg3tNvI0OLvaTCne/sbpSI8AXGVwroASd3JIjNe5+n7vfVv28h8EBsNOrdMdiZmcAvwBcGbstXTGzE4HnMjjrA3ff7+7zX948bYvAcdX1R45nzDXLcuXuN3P09VRGL/n/p8DLgjZKJupFqDGzC4DvuPsdsdsS0GuAz8VuREuC35IjtuoOts8Abo3bks58gMGXjJLvEfIE4PvAR6putivN7ITYjWqbu38HeC+DKvh9wG53vzFuqzr3WHe/DwZfRoDTIrdHKsWEGjP7f6v+3JX/XgpcBrxz2jRyMOV9Dp9zGYOujGvjtbRVQS6vnQozWwd8EniTuz8Uuz1tM7OXADvd/Uux29KxReCZwB+7+zOAhymwm6IaT/JS4PHAjwEnmNmvxW2V9FWj69SkxN1fMO5xM/spBh+2O6obgp4B3GZmz3b37wZsYismvc8hM7sIeAnw/C6u1hhJrVtylMDMlhgEmmvd/VOx29OR5wAXmNmLgWOBE83sY+5e2oFwB7DD3YfVtk9QYKgBXgD8k7t/H8DMPgX8K+BjUVvVre+Z2ePc/T4zexywM3aDZKCYSs0k7v5ldz/N3c9y97MY7GiemWOgmcbMzgfeDlzg7mFu0xdGL27JYYPUfRWw3d3fF7s9XXH3S939jOrz+AoGl0svLdBQ7WPuMbPhzQ+fD3wlYpO68m3gXDM7vtqGn0+BA6JXGL3k/0XAX0Zsi4woplIjAHwIOAa4qapKfdHdfzNuk5rr0S05ngO8CviymW2tHvvt6kqdkqc3AtdWYfwbwL+P3J7WufutZvYJ4DYG3d63U9AVd83sOuDngFPNbAfwLuDdwP8ws9cyCHWdXylX6tEVhUVERKQIxXc/iYiISD8o1IiIiEgRFGpERESkCAo1IiIiUgSFGhERESmCQo2IiIgUQaFGREREiqBQIyIiIkX4/wERfT405S9SVgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1008x576 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "m1 = GPy.models.GPRegression(\n",
    "    X, y,\n",
    "    kernel = GPy.kern.RBF(1, active_dims=[0]) * GPy.kern.RBF(1, active_dims=[1])\n",
    ")\n",
    "\n",
    "m1.optimize()\n",
    "\n",
    "Xnew = np.vstack((Xi.ravel(), Xj.ravel())).T\n",
    "\n",
    "mu, Sig = m1.predict(Xnew, include_likelihood=False)\n",
    "\n",
    "\n",
    "plt.figure(figsize=(14, 8))\n",
    "\n",
    "plt.contourf(Xi, Xj, np.reshape(mu,(50,50)), levels=np.linspace(0, 300, 20))\n",
    "plt.axis('square'), plt.colorbar();\n",
    "\n",
    "display(m1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style type=\"text/css\">\n",
       ".pd{\n",
       "    font-family: \"Courier New\", Courier, monospace !important;\n",
       "    width: 100%;\n",
       "    padding: 3px;\n",
       "}\n",
       "</style>\n",
       "\n",
       "<p class=pd>\n",
       "<b>Model</b>: GP regression<br>\n",
       "<b>Objective</b>: 96.75313670643321<br>\n",
       "<b>Number of Parameters</b>: 3<br>\n",
       "<b>Number of Optimization Parameters</b>: 3<br>\n",
       "<b>Updates</b>: True<br>\n",
       "</p>\n",
       "<style type=\"text/css\">\n",
       ".tg  {font-family:\"Courier New\", Courier, monospace !important;padding:2px 3px;word-break:normal;border-collapse:collapse;border-spacing:0;border-color:#DCDCDC;margin:0px auto;width:100%;}\n",
       ".tg td{font-family:\"Courier New\", Courier, monospace !important;font-weight:bold;color:#444;background-color:#F7FDFA;border-style:solid;border-width:1px;overflow:hidden;word-break:normal;border-color:#DCDCDC;}\n",
       ".tg th{font-family:\"Courier New\", Courier, monospace !important;font-weight:normal;color:#fff;background-color:#26ADE4;border-style:solid;border-width:1px;overflow:hidden;word-break:normal;border-color:#DCDCDC;}\n",
       ".tg .tg-left{font-family:\"Courier New\", Courier, monospace !important;font-weight:normal;text-align:left;}\n",
       ".tg .tg-center{font-family:\"Courier New\", Courier, monospace !important;font-weight:normal;text-align:center;}\n",
       ".tg .tg-right{font-family:\"Courier New\", Courier, monospace !important;font-weight:normal;text-align:right;}\n",
       "</style>\n",
       "<table class=\"tg\"><tr><th><b>  GP_regression.         </b></th><th><b>                 value</b></th><th><b>constraints</b></th><th><b>priors</b></th></tr>\n",
       "<tr><td class=tg-left>  rbf.variance           </td><td class=tg-right>     672398.7240957518</td><td class=tg-center>    +ve    </td><td class=tg-center>      </td></tr>\n",
       "<tr><td class=tg-left>  rbf.lengthscale        </td><td class=tg-right>    12.427228702547053</td><td class=tg-center>    +ve    </td><td class=tg-center>      </td></tr>\n",
       "<tr><td class=tg-left>  Gaussian_noise.variance</td><td class=tg-right>5.562684646268137e-309</td><td class=tg-center>    +ve    </td><td class=tg-center>      </td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<GPy.models.gp_regression.GPRegression at 0x2bed7e8db70>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAHSCAYAAAAUvDo/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3de7RlZXnn++9TtfeGgiqKAsRWoIN2FDUES7v0kPbERLw0MTbG3FqGUdJqGN1HDXq0VZqOdjqXYcR4yfCMOAioMdLYHZQk7SXA8Zgw0gNJSihuFlHbeClBS6AsdgGyd1U95489lywWa+11m3O+z/vO32eMGrX32mvN911zzbXmbz3vO+c0d0dEREQkdxtSd0BERESkDgo1IiIiUgSFGhERESmCQo2IiIgUQaFGREREiqBQIyIiIkVYaLOxJTvSN204us0m12xs9WmuWdjYWlO+0bJc/uGal+s1r/K6luc1fnWo9TluqO90DraxuVNDbNxwuLFlr2dhY5p2Izp4qNnvv4cOz758PzTl58jhye9vh4bf/tCde+5298dM13A7XvC8I/2ee+vfdnfdsnq1u59d+4Jr1ureftOGoznzyJ9vs0kANhx/XOtt+rYtrbW1um1T821sXWxkuQ9trTeJrGypLyit1PgSrm6ub1kHj6nnA+vQ5hGf2DNa3LxS6/L6bd3yYGPLHufEow8kazu1vffXuOGOsH959s+v1QNLUz9m44HJPnMW7hsdtL76W//3N6duuCX33HuYL3zusbUvd9tJe06ofaEN0PBTQ2zfcmttLe5L94EfzdJyfRWDpRpfwsWA+8VJP9wntXpgaaadzCT2L2+aa+c3j733b/7Rv65o6/m2/ZrWvc1LPMWHmhRVmhIt7l9tZLlH7K+3WlC3iMFmvW+Q02riQ76pYAPt7wQHlR5u2nx+876WTW5ndb7HpF165QqRc7Wm7mBTZ7VmbXm1Lq4WCjZplVa9afu5pAg0qtJ0Q4IZtJKrxf2rjc2tqdvSstc6v6Yuiwfqm1+zcN+G2ubXbDywsfY5NqsHlhqbZ9PbKaaca9PTHwZymn+TKpBFDzSq0uRNr16D2pxXA6rWDOrC/JouV2wgRtWmX/QKTur+RXu9pDyq1MhUcqrW1G1pub4jolSxqU+kqk2/weCQoooTKVzVEWhUpZFxFGpkak0FmyP2H6r9EO+6h6EUbGbTdLCBtZ1mtGDTb1TAqCvsRAowg1IFGukehZrCLO57sJXz1jQlh2ATlYJN/GAzTOQwMq/Uw01NVGkinp5BHqZaW8PanlfTlqYO8QYd5h1FjnNsIO05beRhdb4GOtpJJqVQU6CcJww3JfJh3lEnDkNzwUbhpmypA820VKUpR9Ghpssn3msj2ORWrVGwmU1T33jbmiOhYNOuCIFGVZruKjrUSPMUbOpbloJNc1S1aV7d67itbUNVmrIo1LQg1bwaDUPlR8GmWQo3zah7nc6zTahK020KNTI3VWtqXZyCTQsUbOoTaV1Ou33qvDTl0StauLaqNQo2tS6uVjkFG1Vt8tHU+ot4PhoNPeVDoaYlpR7anbPIwabuD9Fcgg2kqdoo3EynqfXV5rCTqjRl0qvaAarWtEfBph4pvq0r3IzX5DqKWKEBVWlyU2yo6fLh3CnlFmzqrtbUrevBRuEmhqbXybyvs6o0aZnZKWb2BTPbbWa3m9kF1e3bzeyLZrbLzHaa2bMHHvcsMztkZr88ZJlHmdlnzOyOapnvmqQvemVblHIIqpQjoXIINpEnDkNewQbSfYPv7ci7HHDaeP5tB5ppqEozsYPAm939qcCZwOvM7GnAu4HfdvftwDuq3wEws43AHwBXr7Pc97j7U4BnAM8xs58b1xGFGqldk9WapijYzGfjgY3FDUf161q4aev5pgg0qtLUz93vcvcbq5+Xgd3ASYADx1R32wrc2fewNwCfBPaOWOYD7v6F6ucV4Ebg5HF90avbIW1Wa3IbhgIFmzqUOBzVr/TqTZvPLfVrOY6qNLMxs1NZq6zcALwRuNjMvg28B7iwus9JwMuAD024zGOBfwN8ftx9dZXultm+ZXzbltTdyF4TV/NuwtIyrNT4ci8egNUaL+pc55W9e5q4wne/Nq72PYnezj+3q4IP03ZIqyPQqEoz3A/d+OrBRj4bTzCznX2/X+Lul/Tfwcw2s1Z9eaO732dmvwu8yd0/aWa/ClwGvAB4P/A2dz9kZus2amYLwBXAH7n718d1UqGmYxb3PcjqtnY+wBb3r7K6dbGx5TcRbJaWnZUt67/Jpl9mvcEmB10JNvDoQJBLyElVbcoh0KhKM9Td7r5j1B/NbJG1QHO5u3+quvk84ILq5z8HLq1+3gF8ogo0JwAvNrOD7v4XQxZ9CfBVd3//JJ0c+yqb2YfNbK+Z3Tbkb28xMzezEyZprC068ml9pQxDNUVHRNWjjQnEEYcwIg9Rpe5bxNdL5mdr6eQyYLe7v7fvT3cCP1P9fBbwVQB3f4K7n+rupwJXAv/XsEBTVXq2sjaMNZFJPs0+Cpw9pLFTgBcC35q0MVmjE/HVR/Nr6pFrsIHYO8r+EJEiTKRse1Bdr5OGnUJ6DvBK4Kzq8O1dZvZi4DeAPzSzm4HfB84ftyAz21X9fzJwEfA04MZqma8d9/ixw0/ufl018WfQ+4C3An85bhkSj4ahxqt7KCqH+TVAdnNs4OEdZpQhqfWMChezDlulDiuTSBlopqWhp+m5+98Boz4s/+WYx/76wO/bq//3rLPMkWaaU2Nm5wDfcfebx03yEYF8g03dogcbyHPycE+kuTbTyiGczCJ1oFGVplumfrXN7CjWSkLvmPD+51dnEty54g9N29zUcplPE2EIqpQT8jWpifk10YeioLnDvdsajoo8JNUlqV+HabdjVWnyN8sn178AngDcbGbfYO1kODea2T8bdmd3v8Tdd7j7jiU7YvaeSiNKmjScy/yatWXWu7xcgg20M4QACjcp1b3u29pmJH9Tf2q5+63ufmLfzOU9wDPd/bu1906Ko2DTHAWb4RRs2lX3+m5r2GnS90/dX0ikXpMc0n0FcD1wmpntMbPXNN+t2eQy9NQTYQgK2h+GUrDpLa/WxQEKNqOoatOOKIFmWhp2KsckRz+dO+bvp9bWG0mmzaOh2tDVicOQz+RheHin1cYkYsjrKKncRAo0TQVxVWni07TwxKJUa9qW40n5II/5NZBXxQbanzOhqk19mqiCtRloVKUpSzGhJrehp4g0DDUZBZtmpAg2Cjfz6dL6U5UmD8WEGsmTgk3/MmtfZCNKCjagcDOLJteZqjQyjyJCTe5VmkhDUCWeu6bLwaapD+2mg43CTUxNr6Ooh27n8oVDCgk1Uq/ShqGapGDTnFQ7OIWb4ZpeJ/O+3qrSCBQQanKv0vREqtZAecGmqWpNUxRs1qT85t4LN10POG2sg7YDjZRLW4KEkWuwyeHEfJB3sEk9LNHFcNPWc07x2k7zXtDQU16yDjWlVGl6ul6taUNOwSaXI6KgnW/KqYMNdKN60+bzq+M1VZVG+mlrkHWVNgwFCjYKNvUoKeCkeC6pAo2qNGVTqJFwFGwGl1n7IhsNNl0YjhqUa8BJ1WdVaKQp2W4VpQ099UQbgoI0w1A6ImpwmbUvstGjP7pWtekXOeD09y1V/1K+bqrSlC/LUFNqoImsxGDT5BFRCjbtBZuo4QYeHSLaDhIRQky/ul4rVWlklLEXtJT22b5lfFvNVzmUoZq88OXSsrOyxWpeZh4XwOxp6kKYgzYe2NjaRTHnNSpczHqRzQhhZRKpA01XqjQrvsDXVx7TwJK/2cAy65ddqFGVJp0UV/Je3L/K6tbFRtvI5YrePQo2w7V9xe+65RJOZpG6mqYT7XVHVjW8LgWaiHNroMxhKMhr4vDacutfZu5DUT2pd6DySHW+Hm1sRzlXaSSjUNOlQBOdgs10FGzWtHFkVE/0uTZdESHQqErTLVmEmq4GmqjVmlR0RNSw5da/zMUDZVVtFG7aV/d6b2ubUZUmf+FDTVcDTU/UYFPi2YYhvyOi1pbbyGKLCTagIak21b2u59lWVKXpntChpuuBJjoNQ01PweZhKYKNwk2zIq3fabddVWnKEDbUKNA8LGq1BhRsZqFg87AU5xtRuKlfU+s04vlocrmAbVeF22I2HH+cAo2MpWAzarmNLLaYCcT9FG7q0dQ6bHPYSVWacoQKNQozo6laM6RdBZsRy21ksY3PT0j1rVzhZjZNrreIFRpQlSYHIbYcVWfyp2DTDaUGG1C4mUaT62nebUBVmm5LGmoUZqYTuVoD5R4R1aTcqjXQTrBRuImp6XXTdqCR8iT55FCYmZ2CzZA2M6/WKNgMl3oIQuHmYW2sixSv9zTvEQ095aHlYyoXFGakEQo2o5bbyGKBbgQbeHiH3sWAk9PzVpVGIMicGpmOqjUj2lWwGbHcRhYLtBdsIoQbyGsnP4+2n2eKYSdVacoU45NCiqNgMxsFm9GiBBsot3qT4jlFel0lf9qaMhW9WgMKNrPKNdh0ZThqUO4BJ2X/63g9VaWRfvE+IWRiCjbrtKtgM2K5jSz2R7o2HDUol4AToZ+pAo2ULeYngxSl5EO9FWwerYvDUcP0B4fUISdSX+oKpbNuZ6rSlC32p4KMlUO1JpU2qjWgYDNMm8EmerjpGQwWTQaMSCGmXy6vlUzHzE4xsy+Y2W4zu93MLqhu/y9m9h0z21X9e3F1+wvN7Etmdmv1/1nrLPsNZvaP1XLfPa4vC/U9LUnF9i3j27ak7sa6Fvc9yOq2Te23u3+V1a2Lrbdbp6VlZ2WLNbBcWGlws1k8AKubm1t+v4X7NnDwmMPtNFazSKGjSXUGmjaqNDKVg8Cb3f1GM9sCfMnMrq3+9j53f8/A/e8G/o2732lmpwNXAycNLtTMnge8FDjD3R8ysxPHdUSxuRA5VGw0v2Z2TVZsSjgyCvKq2nRNjoFGQ0+Tc/e73P3G6udlYDdDQkrf/W9y9zurX28HjjSzI4bc9T8A73L3h6rH7R3XF30CSKsUbGbX5IdsCUdG9SjYxKLXo1vM7FTgGcAN1U2vN7NbzOzDZrZtyEN+CbipF1wGPBn4aTO7wcz+1syeNa59DT8VJIdhKCh7KOqI/Yd4aGtzwwlNDUWtLbus4Sgg2yGpUtQdaFSlGe+Hhxf52kOPbWLRJ5jZzr7fL3H3S/rvYGabgU8Cb3T3+8zsj4HfAbz6/w+BV/fd/yeAPwBeNKLNBWAbcCbwLOB/mNkT3X3kC6RQUxgFmzHtKtiMWXY5wQYUblKKEmikNne7+45RfzSzRdYCzeXu/ikAd/9e39//BPh03+8nA1cBr3L3/z1isXuAT1Uh5u/N7DBwAvD9Uf1QXVCS0VDU7HIdioI0OyfNt2lPE+t6nm2mS1WaVMzMgMuA3e7+3r7bH9d3t5cBt1W3Hwt8BrjQ3f/XOov+C+Cs6jFPBpZYm2Q8kt7lBcph0nAXKNiMlupbt4JNs6KtXx3t1JrnAK8Ezho4fPvd1WHbtwDPA95U3f/1wI8Dv9V3/xMBzOxSM+tVhD4MPNHMbgM+AZy33tATaPipWBqGGtNuS4d6ayhqtF6waXM4CjQk1ZSmAk2bAVhVmtm4+98Bwz6IPjvi/r8L/O6Iv7227+cV4Nem6UusWC21yqViU/IwFDRfsWlSG990U1ZtolUWctTkemxz2EnKoHd04RRsxrRbQLBp+ttlycEGFG7m0eR6a3ubUJWmDHond4CCzZh2FWwmWH6jiwfSH92icDO5ptfVvNuCqjTdpXewhKJgM7tSgo3CTVxtrJsUgWaa907Ow8ldoHduR+RSrQEFm3mUEGwgfbABhZt+ba2LCK+75E3v2A5RsJmg3ZaCTZNKCjYRdnK9HXpXA05Oz1tVGslna5VaKNhM0K5OzjfB8rtVtenpUrhp+7lGep0lX914d8ojKNhM0K6CzYRtNN4EEG+HV2r1JtXzquP1VZVGQKGmsxRsJmhXwWbCNhpvAogzHDWohICTsv+pAo2UaexWXF0ufG91muLebReb2R3V5cSvqq7jIJnJKdikomAzaRuNN/EjEYNNT04BJ0JfU76WqtKUaZKt+aPA2QO3XQuc7u5nAF8BLqy5X9KSXIJNqmpNW9oINqVMIIa4VZt+/aEhQsiJ1p+6Xj9VaaTf2C3b3a8D7h247Rp3P1j9+kXg5Ab6Ji1RsBnTbgGHeveUFGwgfrDpNxgqmg4X0UJMv9SBRmcPLlcdF7R8NfDfa1iOJKQLYI5pt5ALYEKzF8FcW36zF8IclOrCmHWJFjialjrQTEtDT3mZ691kZhcBB4HL17nP+Wa208x2rhwuewhB2qGKzfxKOuS7J6eqTRdFGTJUlaZsM4caMzsPeAnwCncfuZW4+yXuvsPddyxtaP8btkwul2EoULCpQ2kTiCHOjlMeqe7XRFUaGWWmUGNmZwNvA85x9wfq7ZKkpGAzQbsKNlO20XgTj6JgE0ekQKMqTfkmOaT7CuB64DQz22NmrwE+CGwBrjWzXWb2oYb7KS1SsJmgXQWbKdtovIlHUdUmvUjrf9rtXFWaPI2dKOzu5w65+bIG+iKB5DJxGDR5uA5NTx5ea6PdCcQ9uU8kzlUTgUaHb8s4dRz9JIVSsJmgXQWbKdtY+z9VuFGwaUe0QNOlYacVX+CbDx6fuhvJdOtYQpmahqImaFdDUTO000ozj6IhqWY1tX7b3l409JQvhRoZS8FmgnYVbGZop5VmhlK4qV/U9dmlKo0o1MiEFGwmaFfBZoZ2WmlmJIWbejS5DtsedlKVJm8KNTIxBZsJ2lWwmaEdhZtcNb3eUm8Xkh+FGpmKgs0E7SrYzNhWa02NpHAzmTbW07zbg6o03aRQI1NTsJmgXQWbGdtqral1KdyM1sZ6ibIdSH4UamQmtm85m3CjYDO/LgYbULjp19a6qOP1V5WmuxRqZC4KNmPaVbCZsS2Fmyhye+462qnbFGpkbgo2Y9otLNh0tWoDD+/gc9rJzyrF80z1eqtKUw6FGqmFgs2YdgsKNtDd4ah+pYabVM8r1bDTtNp6L8tsFGqkNgo2Y9pVsJmjrfjhJueAk/o5pHxtVaUpi0KN1ErBZky7CjZzttdqc1PrDwc5hJwI/azrNVWVRkChRhqQS7ApXcnBJnq46YkWcKKFrtSBRlWa8rQbahaavcqwxJFDsCm9WgPlBpu1Nltvci6DgaLpYNFmW7PI7vVTlSYLC2036Nu2AHns9GQ+tm/5R693VIv7HmR126b2292/yurWxVbaOmL/IR7a2vwXiqVlZ2WLNd7OI9tc+38l9ma2rknDxurm6e4fWZ2BRlUa6Zds+Cn6zk7qkUN4VcWmPqnOEZLbt/5ZRKy2zCJCoJFyJZ1T49u2KNx0gILNOu0WGmxSDUd1IdzkLMrrM+17QUNP+QgxUVjhpnw5XFZBwaZeqtpIv7pfF1VpZJgQoaZH4aZ8CjYj2lWwqbldhZtIIgUaVWnqZ2anmNkXzGy3md1uZhdUt19sZneY2S1mdpWZHVvdfnx1/wNm9sF1lvsr1fIOm9mOSfoSKtT0KNiUTcFmRLsKNg20naxpqeg16ISDwJvd/anAmcDrzOxpwLXA6e5+BvAV4MLq/j8Efgt4y5jl3gb8InDdpB0JGWpAVZvSKdiMaFfBpoG2tWNNoan1ripNPO5+l7vfWP28DOwGTnL3a9z9YHW3LwInV/e5393/jrVws95yd7v7P07Tl7ChpkfBplwKNiPaVbBpqH2Fm7Y0tZ5Tb0Mdd4KZ7ez7d/6wO5nZqcAzgBsG/vRq4HPNdjHBeWpm4du2hN8Bymyin8tG57GpV2+n1Pb5bB7Zh7zPaxNd1EDTlfPSrBzeyLcfOLaJRd/t7uvOazGzzcAngTe6+319t1/E2hDV5U10rF/4Sk2PhqPKFT2wqmJTv9TfuFW1qV/kdTrL9q2hp+mY2SJrgeZyd/9U3+3nAS8BXuHujb/xswk1PQo2ZVKwGdGugk3DfYi7I85J0+swwrYio5mZAZcBu939vX23nw28DTjH3R9ooy/ZhRpQsCmVgs2IdhVsGqdwM7vogUZVmlY8B3glcJaZ7ar+vRj4ILAFuLa67UO9B5jZN4D3Ar9uZnuqo6Uws0t7h2+b2cvMbA/wU8BnzOzqcR3JYk7NMJpnU6beaxo1uHZljg3QqXk2PSVcS6otbYTAKKFX1lcdyTTsDfzZdR5z6ojbX9v381XAVdP0JctKTU/UHZ/ML3Jg7ULFBrpbtQFVbsbJZd10ZXKwPCzrUAMKNiVTsBnSroJNqxRuHqnN9ZFqW9DQU96yDzWgYFMyBZsh7SrYtK7r4abt51/HNqAqTTcVEWpAwaZkCjZD2u1AsIkcbroScFI815Svu6o0+Ssm1ICCTckUbIa0W3iwgZhVm56Sw02q51bX660qTXcVFWpAwaZkCjZD2lWwSa6k6k3K55E60KhKU4biQg0o2JRMwWZIux0JNtHDDeQZcHLss8goRYYaULApmYLNkHY7EGwgftWmX+SwEK1vqtJIXYoNNaBgUzLbtxw23CjYNCunYNPTHyLSTL6NF2R6cnw9Ja5szyg8KZ15uGxRr/LdhTMPQ7tX+O4X6SzEsxoVLuY9m3G00LKeOgNNW1WaVF9aZDLFhxpQsCmdgs1AuwmCDbR3WYV+S8uedbAZJqdQMo8IgUbKU/TwU7+IOz2pT9TQ2pWhKEg7HKUhjLxEeb1UpSlPZ0INKNiUTsFmoN0OBRuIs6OU9dX9OqlKI/06FWqkfAo2A+12MNgo3MQVKdDoiKcydS7UqFpTPgWbgXY7FmxAVZuIcn9NNPSUh86FGlCw6QIFm4F2Oxpsct+RlqKJ10FVGhmmk6EGFGy6QMFmoN0OBhvIv0KQu2iBZhaq0uSjs6FGuiHqSfoUbNqlqk37oq5zVWnK1onz1Iyi89d0R8Rz2XTlPDaQ9lw2/Uo4aV8OmgwzEUJyZAcPbWDv/ZtTdyOZzldqou3opDkRA2yXKjYQZ4cUsYJQitICjYae8tL5UCPdomDT1+7+1c4OR0Hc4ZGcRV+fGnoq39hQY2YfNrO9ZnZb323Hmdm1ZvbV6v9tzXazWarWdIuCzUDbHQ42oHBThzbW4bzbzCzbuao0+ZmkUvNR4OyB294OfN7dnwR8vvo9awo23aJgM9B2x4MNKNzMqo11Fm1bkbjGhhp3vw64d+DmlwJ/Wv38p8Av1NwvkcYp2Ay0nSjYRNthKdxMLpf1pGGn7ph1Ts1j3f0ugOr/Eyd5kG+MfcSBqjXdo2Az0HbHJxD3y2WHnUKbwS/VtqGhpzw1PlHYzM43s51mtnNl9f4kh7BOQ8GmexRsBtpWsPmR3s5bAWdN2+uijm1CVZpumTXUfM/MHgdQ/b931B3d/RJ33+HuO5YWjwYIH2ykexRsBtpWsHmUroebtp97ym1BVZp8zRpq/go4r/r5POAvp11A5GCjak03RTz7cFeDjcJNHDk/X1VpumeSQ7qvAK4HTjOzPWb2GuBdwAvN7KvAC6vfpxY52Eh3Kdj0tZ1wpxA52EDeO/tJpHx+0V97iWvsZRLc/dwRf3p+HR1Y3bYpZKlPl1DotmiXVUh1SQVIc1mFniP2H0p+aYVx+nf8JVx+IXVQqyvQzBrII+6PZHIhzigctWITaacm7YsWalWxiS/XicVR+p3Tay0xhQg1EDfYSLcp2PS1nTjY5LbDixIU1hOpf3W+vqrSdFeYUAMxg42qNaJg09d24omXuQWbnv6AkzJEROmHSFNChRqIGWxEFGz62lawmdtguGgiYLTRRl0iVGmkDGMnCqcQbfKwJg0LaPLwI9pOOHkYHt4JRp9EPI3IoaNJUUJqpH2OzC5cpaZHFRuJKFq47XLFBuLsEGU2db9+EbZJSStsqIFYwSbSN3RJS8Gmr+0AOxEFmzxFCjSq0szPzD5sZnvN7La+27ab2RfNbFd1uaRnV7dvM7OrzOwWM/t7Mzt9xDI/amb/VD1+l5ltH9eP0KEmGgUb6VGw6Ws7SLBRuMmHXqsifRQ4e+C2dwO/7e7bgXdUvwP8J2CXu58BvAr4wDrL/Y/uvr36t2tcJ8KHmkjVGpF+CjZ9bQcINqCdZQ6aeI1UpUnP3a8D7h28GTim+nkrcGf189OAz1ePuwM41cweW0c/wocaiBVsVK2Rfgo2fW3vXw0RbhRs4tJr0zlvBC42s28D7wEurG6/GfhFgGpI6seAk0cs4/eqYar3mdkR4xrMItRArGAj0k/BZqD9IMFGO9BYmno92q7SRHu/t+iEal5M79/5EzzmPwBvcvdTgDcBl1W3vwvYZma7gDcANwEHhzz+QuApwLOA44C3jWsw5CHdo0Q51FuHeMsgHe490H7iQ757crh2VBdEDDSlOnR4A/uXG3nv3+3uO6Z8zHnABdXPfw5cCuDu9wH/DsDMDPin6t8juPtd1Y8PmdlHgLeMazCbSo1IdNGCbuovAFF2OKrapFXSuo/2Hs/AncDPVD+fBXwVwMyONbOl6vbXAtdVQecRzOxx1f8G/AJw2+B9BmUXaqIMQ0X6Vi5xRPvQU7B5WEk711w0uc7n3bZSvzdKY2ZXANcDp5nZHjN7DfAbwB+a2c3A7wO9IaunAreb2R3Az/FwNQcz+6yZPb769XIzuxW4FTgB+N1x/chq+KknyjCUyDAaihpoP8hQFJR5JuKoFCK7xd3PHfGnfznkvtcDTxqxnBf3/XzWtP3IrlLTE6FiE2nHJbGoYjPQfqCKDWiH27Sm12+KKk2097QMl22oEYku2oeggs0jaa5NM6IHGilb1qFG1RqJTsFmoP2AOySFm3rksh5TvwekWVmHGogRbETWY/uWQ4Wb1B/qUU7SNyiXnXJEba23VNtNpPevrC/7UCOSi0gfjKmDDcSs2oDm20xL60siKSLUpK7WaAhKJqVgM9CHwMFGO+vx2lxHdWwrEbZ5aVYRoQbSBxuRSSnYDPQhaLABhZtR2l4vKbeRSO9XGa+YUJOaqjUyjUgflAo24yncrMl5PUTYzqV5RYUaVWskJwo2A30IHmwg7536vFI97xy2C4mjqFADaYONqjUyLQWbgT5ksgPrUrhJ+Vzr2h5m3bYjvT9lMq2GGt9obTYnkoVIH5xRgiuYuagAACAASURBVE1u4abUgJPyeeWyDUgsrVdq2rgGjKo1MRy+595H/JPRFGweLbedWinhprSgpipNtyS5oOXq1sXGP7B00cs01gsv6/1tw/HHNdGdrES6EGbqi2D+qB+BLoY5qf4wkNOFMyOFmNwCrcSR7CrdbQSbVHzbls6l/HkrMf2P73LAUbAZ0o8Mg01P9IATKcj01Llf0Bfb7kkWatqgak3zmhhW6nrAUbAZ0o+Mg03PYIBIGXIihhmIU6Hp2pfSkiQ9+qnk+TVRdkpNamOeTFfn40T6UI3yxSCnCcSTaGvuSn87Jc2VGSfKdivtSl6pKXkYqmRtB41ee12q3KhiM1wJVZtB44LGelWdUkKK9gNSh+ShBpoPNhqGKkfXwo2CzXAlBpv1lBJcRqn783+ez/tIVVKZXnEn3xslxYdxlJ1R3SIMB3V1WCq1SF8O9M1eRAaFqNSAhqFkNl2o3ESq1kC8ig20Mz9PmhGpSlMCP2SsHlhK3Y1kQlVqmv5gUrVmflGrI1H7VZdoJfFoOw59IcpTtNct2vtMphcq1IC+ccnsSh+SivaBq2Aj82ji9Yq2TUr7woWapkUpm+col8BQcrhRsFmfgk0eIr5O0d5bMpuQoaa0YajShqByoWDTjojBJuJOU9botZEmhQw1oGEoqUepVRsFm/G084ynydck4jYo7Qsbapqmak23KNg0L+JORcEmjsiBJtp7SWYXOtSoWhNHCaGgxKpNtA/jqMFG4UakG0KHGmg22Kha000KNs2KGGxAVZuUIldppCzhQw2oYiP1K61qo2AzGVVt2hd9fUd778h8sgg1TdIh3t2mYNOcqMEG4u9oS9H0eo68jUka2YSaUoahNAQVj4JNcyLvdFS1aVYO6zba+0Xml02oAQ1DSXNKGo6K9kEdOdhAHjvf3LSxTqNvV5JGVqGmSarWjFbKzn4SpTxXBZvpqGpTjy6sx1I+I0qVXahRtUaaVsqHloLN9ErfITepzXVXx7YU7f0h9Zgr1JjZm8zsdjO7zcyuMLMj6+rYepoKNpo0LD0KNs3IJdgo3ExH60uimDnUmNlJwG8CO9z9dGAj8PK6Ola63IaguqiUeTYKNrPRjnoyba+nlNtPCZ8HpZt3+GkB2GRmC8BRwJ3zd2kyqtZIW0r4IFOwmY2qNuvLNdBEez9IfWYONe7+HeA9wLeAu4D97n7Neo85vNFmbW6o3OfXqFqTDwWb+i3ue1DhJmNaH9LPzD5sZnvN7La+2/6LmX3HzHZV/15c3f6Kvtt2mdlhM9s+ZJlPN7PrzexWM/ufZnbMuH7MM/y0DXgp8ATg8cDRZvZrQ+53vpntNLOdBx+6n4e2bpy1ydaoWiPDKNg0I5dgAwo3kG4d5LSddNRHgbOH3P4+d99e/fssgLtf3rsNeCXwDXffNeSxlwJvd/efBK4C/uO4Tswz/PQC4J/c/fvuvgp8CvhXg3dy90vcfYe771g44ug5mhtO1Rppk4JNM3LbYXU13JTwnGfd/kt47zfJ3a8DZllJ5wJXjPjbacB11c/XAr80bmHzhJpvAWea2VFmZsDzgd2TPFDVGslZKROIo8kt2EAZO/lJpA5xOW4b8iOvN7NbquGpbUP+/m8ZHWpuA86pfv4V4JRxjc0zp+YG4ErgRuDWalmXTPr4OoNN7tUayVPOwSZitQby3Hn1dvilBpySnlfU7b5Wh42NBzbW/g84oTeVpPp3/gS9+WPgXwDbWZt7+4f9fzSz/wN4wN1vG/JYgFcDrzOzLwFbgJVxDS5M0KmR3P2dwDtnffxDWzdyxP5D83ThR1a3Ltb+5lvdtqmVD1nftqUbb7YCHb7nXjYcf1zqbszE9i2HHP5c3PdgtpXS3mdQCV+0ooSZCEE35y8wNbrb3XdM8wB3/17vZzP7E+DTA3d5OaOrNLj7HcCLqsc/Gfj5cW1md0ZhaV+uO+225PyBFzVMR9iRzSPnyk2kvue+HXSdmT2u79eXsTac1PvbBtaGlD6xzuNP7LvvfwY+NK7N5KEm+jBUrt8YpV0KNvUrYYcWKSCsp/QhNIi7nZfCzK4ArgdOM7M9ZvYa4N3V4di3AM8D3tT3kOcCe9z96wPLudTMehWhc83sK8AdrJ0H7yPj+jHX8FNdNAylIagSaCiqfjkPRfUb/EyKMjwVOcSUEGq7xN3PHXLzZevc/2+AM4fc/tq+nz8AfGCafiSv1IiURBWb+pW4c0tZGcmhKlP3az7Ptp3ze7qLwoQaDUPpnDWlyPlDUMGmff0ho4mw0eSyRaIJMfzUU+cwlEhKGoqqXylDUZPocvgoOcBK88JUauqWa7VGyqGKTf20w5NpaeipW8KFmhzONtykiN+QZXY5fygq2Ejb9NrKvMKFGqgv2KhaIxHkfFkFBRtpSxOvadTtV5oTMtTUKcqhk9NQtaZMuQabqBb3PahwIyKPEDbURB6G6mK1JtdJr9HkGGyif9tVsMlfxNcwx/eqBA41EHsYSmRWOX5YKthIU5p67aJvs9KM0KEmsqarNRqCKpuCTf0UbEQkfKhRtUZKpWBTPwWbvESt0uT43pQ14UMNdDfYqFpTvhw/PBVspA56naQJWYSaqLo4YVgEFGxkPnp9pCnZhJquVmsi0RFQzcixWgN5BBvtPLsn+nYpzcom1EDMw7w1YVjqoGDTHAWbWKK/Hrm+F2VNq6HGg2QSVWskolw/TBVsZFJNvw45bIvSrNYrNStbbK7Hq1ojJcs12ORAwSYtrX9pQ1bDT3VStUaiyjHY5PINWTtWWU+O7z15pIUUja5sMZaWfebHP7R1I0fsP1Rjj+a3um2TPjClNofvuTe7idm2bzmLqmLvfaqjF9vTxmdjLsG6aXYIFu7rbL0iXaUmwjBUTtWaKDuL3Ha0OcvxW2NOOxZ9CWmH1rO0qbtxrgH65id1U7Bplna4zWpr/ea0zUmzkoYaVWtExlOwaZaCTTNyW685vs/k0ZJXauYNNnWoM9g0Wa2JMgQl7cvxAze3YJPbTjiyNtdlTtuZNC95qJlXxEO8RZqgYNM8BZv5aR1KSiFCTWnDUKVXazRZOB0Fm+Zpp9w9Ob6vZLgQoQZiDEOJ5EAfwM3TcNRs2l5nuQVmaV6YUDOvLlVrRHKT685HwWZyWlcSQahQo2rNZCIMQUlaOVZrcg422mGvL8X6yXV7kmaFCjUwX7BRtUa6RMGmXQo2w+W+XnJ8H8lo4ULNvLpyNFTqao0mC8eQ4wdy7sEm9514nVKti5y3IWlWyFCTehhKJ+STnCjYtK/rwUbhTqIKGWrmFalaoyEoaYOCTfu6ulMv6Xnn+L6R9YUNNarWjJd6CEpiyfEDuoRgU9JOfpwIzzX3bUaaFTbUQPpJw3UptVqjeTVShxJ2Ul0IN6U/PylD6FCTmqo1kpscqzUlKTHcRHpOJQRgaVb4UKNqjch0cgw2pe2sIgWBeZTwHEbJ8X0i44UPNZA22ORQrREZlOMHdmnBBvINBRFDWYnbh9Sv1VDjcQonSTRVrUk5BKV5NXEp2MTQCwjRQsIwufRT4jGzD5vZXjO7re+2i83sDjO7xcyuMrNj+/52hpldb2a3m9mtZnbkkGX+TvXYXWZ2jZk9flw/Wq/UrMy4/1W1RmR6CjaxRA0NUfslWfkocPbAbdcCp7v7GcBXgAsBzGwB+Djw7939J4CfBVaHLPNidz/D3bcDnwbeMa4TWQw/lURza0TGKznYQJzqTYQ+TKLu7SHHsB+du18H3Dtw2zXufrD69YvAydXPLwJucfebq/vd4+6Hhizzvr5fjwZ8XD+ShBpVa+qno6BklFw/wEsPNj0pgkUuYSaiq+//WOou5OrVwOeqn58MuJldbWY3mtlbRz3IzH7PzL4NvIIJKjULtXR1BitbYGmGz6yVLcbS8tiwJi3acPxx2e44u+LwPfdmOf/J9i13JrAPCxl1VXZzDjBdCbd1scOweKCRRZ9gZjv7fr/E3S+ZqE9mFwEHgcurmxaA/xN4FvAA8Hkz+5K7f37wse5+EXCRmV0IvB5453ptJQs1KTy0dSNH7H9UhWtiq1sXWdw/bNhvyuVs29TIh4xv26IPABkp12DTZaM+J8aFnZxDjIR1t7vvmPZBZnYe8BLg+e7eq0jsAf7W3e+u7vNZ4JnAo0JNn/8GfIYxoSbpnJoUw1AikhcF9Ufrn5Mz7F8pmnjtVVVuj5mdDbwNOMfdH+j709XAGWZ2VDVp+GeALw95/JP6fj0HuGNcm8knCs8abGYVZW5NaROGVQHIQ64f6Ao2IrGZ2RXA9cBpZrbHzF4DfBDYAlxbHZb9IQB33we8F/gHYBdwo7t/plrOpWbWqwi9y8xuM7NbWJtcfMG4fmQ7/KS5NcNpCErGyXUYqkvzayRmkNUk4dHc/dwhN1+2zv0/ztph3YO3v7bv51+ath9zVWrM7Fgzu7I6uc5uM/upWZajao1Iu1SxkS7KdbuXyc07/PQB4K/d/SnA04Hd83dpcqVcF0okhVw/4BVsyqfXWGY1c6gxs2OA51KVl9x9xd1/MOvycps0HLlak6pEn+OQhuRJOz0RGWaeSs0Tge8DHzGzm6rJPUfP05nchqFEcpdrtQYUbEql11XmMU+oWWDtuPI/dvdnAPcDbx+8k5mdb2Y7zWznoQfun6O50XKv1jRBEyplUgo20gXzbueaJJyHeULNHmCPu99Q/X4layHnEdz9Enff4e47Nh41vpDTxWpNSROGNQSVp5yDjZRDIVXmNXOocffvAt82s9Oqm57PkJPntEXVGpFu0o5QRHrmPfrpDcDl1YlxtgO/P3+XVK2pi4agZBo5V2sUbPKn11DqMNfJ99x9FzD1tSAmMcsFL1OdkK+ua0KJpJbriflAJ+eT0XIO7DKd5JdJqNusw1Cq1tQn152irMl5B6Bv+yLdFjrUtD0MNQ/NrRGJQcEmP9FfMx35lI9WQ4231JqqNSLzyblaA/F3ktKe3LdlmU7rlZrVzdPdX9Wa+WkISmaR+85AwSYPep2kTqGHn3pmCTY5V2tEpB7aYcam10fqliTUTFutyUkd1ZpSJgxL/nKv1ohIt2RRqQFVa0RSyT3YqBoQUxuvS+7brkwvWahRtWbMMgqYMKx5NeXIfeegYCOz0pFPeUlaqWlj0rCqNQ/TEJR0mYJNHHotpClznVE4hVnONJyr1W2bWNz3YOpuiAB5n224R2cdltLZoe7sI4dJPqemjWGoFNWaqId3ty33naA8Uu7DUKAqQWptrf8StlWZXvJQM4s2h6FKo2+pIgo2IqUKEWoiTxpOXa0pYcKwlKWUb8AKNu3TOpemhQg1EHvScGnartZoCKo8CjbSBTryKT9hQk1kqtaIlEvBph1trudSQrdML1SoUbVGJB8l7TgUbETKECrURFZatUZDUFIHBRuZhNattCVcqFG1RkRS0c63flqn0qZwoQbaCTazSF2tqZsO75Y6lFStAe2Ec1fa9ijTCRlq2pBjtSb3CcMagiqXdiQyjAKitC1sqFG1RkRS0c5YdDh3nloNNR7sGpGq1mgISupTWrVGwWY+Wn+SQuuVmoPHHJ74viVWa7pOQ1BlU7AR0HqTdMIOP7Wl7WqNDu8WyYt20PkoLVTL9JKEGlVrRMpV4o5FwWZyWleSUhaVmqYveJljtSZnGoKSHGlnPZ7WkaSWLNRMU62ZVheqNRqCkshKrNbA2k5bO+7htF4kgqSVmiaHoaalao1IvUoNNqAdeEQlb28yuSyGn2ahas302qzWaAhKcqdg87DS1oXOUTM9M7vAzG4zs9vN7I0Df3uLmbmZnTDkcc8zs119/35oZr8waz+Sh5pIk4ZVrRGpV+nfnkvbmYvMwsxOB34DeDbwdOAlZvak6m+nAC8EvjXsse7+BXff7u7bgbOAB4BrZu1L8lBTgpKqNSIyna4Hm64/fwHgqcAX3f0Bdz8I/C3wsupv7wPeCvgEy/ll4HPu/sCsHQkRalStiUNDUFK30qs10N0de5Tn3YVtLIATzGxn37/z+/52G/BcMzvezI4CXgycYmbnAN9x95snbOPlwBXzdHJhngensroZFg+k7sUjPbR1I0fsP5S6GyIhHb7n3uJDrO1b7tRRhFECjTySHYKl5UmKIlO72913DPuDu+82sz8ArgUOADcDB4GLgBdNsnAzexzwk8DV83QyRKUGYh3inds1oTQEJRJDV3b0XXmeMjl3v8zdn+nuzwXuBb4BPAG42cy+AZwM3Ghm/2zEIn4VuMrdV+fpR5hQA7EO8W6ThqCkC7oyRKAdvnSRmZ1Y/f/PgV8EPubuJ7r7qe5+KrAHeKa7f3fEIs5lzqEnCBZqmtRGtUYThkXWp2CTv2jPrSvbVAY+aWZfBv4n8Dp33zfqjma2w8wu7fv9VOAU1iYYzyXcnJqDxxxm4b7JslbEuTWzWt26yOL+uaputfJtW8J9eInkpPf+KWmeTRc+E3SOmtm4+0+P+fupfT/vBF7b9/s3gJPq6EdnKjWgak1UGoLqlq59sy4lCJTyPKRsIUON5tbEUNI3TJGUcg8EufdfuqPdULOhkcPMplJ6tSZXqtZ0S9eqNZBvMIjc7y5uR7K+1is1hzZPdi4XVWtmfHymQ1AiXZDbVb5z6qsIBB1+6mnq3DWq1kxOQ1DSlC5/y84hLOTQR5FBSULNpNWaaZRUrZlXrtUaDUFJl0Su2kTtl8g4oSs10M1qjSYMS1d0uVrTEylARA5aIpNIFmpUrWmWqjWSCwWbGGEidfvT0nYjwySt1KSeNKxqjYhEkiLcRAhUUejEe/kLP/wkMWgISpqkb92P1FbQUJiR0iQPNarWjNbVw7s1BCWyphdu6g4fqs5IqcJd+0ni0vWgpEmH77lXgXYd/e+9aSunet9KV8wdasxsI7AT+I67v2SWZRzafIiNB8ZXM5q62OXKFlia4j2/ssVYWm7n7MjzXuhyddsmFvc9WGOPRCS1rocUDVfKKHUMP10A7K5hORNp6hDvppVyMr626Bt7N2lnJSLzmCvUmNnJwM8Dl87bkdSHeLcxt6YEmjAsIiJRzVupeT/wVqDV8knXqjVdnTAs3aRqjYjMauZQY2YvAfa6+5fG3O98M9tpZjsPHbh/3WWqWpOHtqo1GoISEZFpzFOpeQ5wjpl9A/gEcJaZfXzwTu5+ibvvcPcdGzcfPXahTRziHYmqNSLjqVojo2jbkPXMfPSTu18IXAhgZj8LvMXdf62mftWqySOhpFkbjj9OH2IiIhPacMg5Yn/9ox65aPXke7ZxssOgc6rWdPVkfJowLE1TmJU26RIJZagl1Lj730x6jprFzSt1NDm1JufWiIiISHrJL5Mwiqo10qMJw92lao2ITCNJqFG1ZnYaghIRERkubKUGyq/WyORUrekuVWukR9uCjJMs1KhaownDIiIidQpdqQFVa0RE39BFZDJJQ42qNemqNXXSGYZFRCSC8JUaULWmCTrDsORG1RoRGSd5qJm0WlP3daFKqNZ0kao1IiIySvJQA/UOQ3WpWqMJw9I1qtZ0l157mUSIUDMpVWtERERklDChRtWa2XSxWqMhqG7TN3YRGSVMqJlUymqNiIiIxBUq1ESv1kQdgop0eHdbVK3pNlVrRGSYUKFmUrlUa7p4eLcmDIuISCqthpqNG8ZXT1St0YRhEZF+qszJpLKs1ICqNYM0YVi6Rjs6qcvV938sdRekJq2Hmq1bHhx7H1VrVK0REZF8mNmxZnalmd1hZrvN7KfM7HfM7BYz22Vm15jZ44c8bruZXW9mt1f3/bfz9CPbSg2oWjNIE4ala1StEQnjA8Bfu/tTgKcDu4GL3f0Md98OfBp4x5DHPQC8yt1/AjgbeL+ZHTtrJ5KEGlVrYspxCEpERNIys2OA5wKXAbj7irv/wN3v67vb0YAPPtbdv+LuX61+vhPYCzxm1r5kXamBcqs1OrxbZDKq1oi04gQz29n37/y+vz0R+D7wETO7ycwuNbOjAczs98zs28ArGF6p+REzezawBPzvWTu5MOsD57V1y4PsX16/MrC4eYXVA0u1tHfwmMMs3FdvhlvZAkvLtS4yudVtm1jcN76SNgnftgXb1/wK2nD8cdqxiYgAdshZ3L/axKLvdvcdI/62ADwTeIO732BmHwDeDvyWu18EXGRmFwKvB945tN9mjwP+DDjP3WceXklaqalrGKruak0UmjAsIl2nLyxZ2APscfcbqt+vZC3k9PtvwC8Ne3A1fPUZ4D+7+xfn6Uj2w0/TmHRuTVMXutSEYZFmaMcnko67fxf4tpmdVt30fODLZvakvrudA9wx+FgzWwKuAj7m7n8+b1+ShxpVa9aXolqT44RhHQUlIpLUG4DLzewWYDvw+8C7zOy26rYXARcAmNkOM7u0etyvsjbJ+NerQ793mdn2WTuRbE5NKpPOrVndDIsH6m9/ZYuxtPyoCeC1W9262NS4qkhIh++5V+FWJBF33wUMzrkZOtzk7juB11Y/fxz4eF39SF6pgfyrNTq8e32q1oiISBtChJq2NTG3Zho6vFtERKR+YUKNqjXx1FmtaYuqNd2mCcMi3RYm1LStK9WaKHSGYRGZlkKqTCtUqFG1pn5dHIJStabbtCMU6a5WQ83CxvqvwTSP1NWaaenwbhERkdFCVWqgW9UanYxPRESkPq2HmhOPbuDkL3PIrVqTQo7VGg1BdZuGoES6KVylBlStWU/uE4ZFRESakiTUqFrTvi4OQalaIyLSLSErNZB/tWYaOVRrchyCkm7TEJRI9yQLNblWayalw7tFRETaFbZSAzGrNVGGoHKv1rRFQ1AieVKlTWaRNNREq9bULeLh3VFoCEraoB2jSLckr9SMCzZtVms0YbhMqtaIiHRD8lBTutIO79aEYRERiSpEqFG1pl2q1kiXaAhKpDtChBp5WNeqNSIiInVZSN2BnhOPPsDe+0eXP7ZueZD9y+vvTBc3r7B6YGnuvhw85jAL943Pe6ubYXGCuc4rW2Bpee5u1Wp16yKL+1dTdwPftgXbF2zliIhkyg45i/vGj26UqnOVmhxOxjctXTphMhqC6i4NQYl0Q6hQU8fcmrqkPBlfDod3a8KwiIhE02qoWdrQfJWkqyfjm4UmDIuISElar9ScctQP1v27qjW9+3ZrwrCqNdI0DUGJlC/U8FNdVK2ZnKo1IiJSiplDjZmdYmZfMLPdZna7mV0w6WNzqtbUrbRqjYhI7v710a9K3QWpyTyVmoPAm939qcCZwOvM7Gn1dGt+OhlfO3IdglK1RiQuDRXKrGYONe5+l7vfWP28DOwGTpr08V2u1kTTxSEo6SbtLEXKVsucGjM7FXgGcEMdy6tL1GpNaUNQuVZrRESkLHOHGjPbDHwSeKO73zfk7+eb2U4z2/nDfT98xN9UrYmji9UaDUGJiJRlrlBjZousBZrL3f1Tw+7j7pe4+w5333HktiPnaW4mbVdrJlVatUZERCS1eY5+MuAyYLe7v3fW5ZRWrdGE4flpwrA0SfNqRMo1T6XmOcArgbPMbFf178U19atWk1RrJlF3tSaaLg5BiYhIOeY5+unv3N3c/Qx33179++wsy4pQrUlxMr7ShqBUrRERkZSKPKPwMHVVa0qnao2IiOQqTKgZV60Zp61qjQ7vFsmf5tWIlClMqBln3BCUxKAhKBERSaXVULNkB9f9e9PVmhIO726DhqBERGRSZnakmf29md1cXQvyt6vbX29mXzMzN7MT1nn8Pzeza6prSX65OqHvTLKp1EB+1ZomDu/OYQhK1RoRkU55CDjL3Z8ObAfONrMzgf8FvAD45pjHfwy4uLqW5LOBvbN2pPVQ82Ob7ln37/NWa8aJWq2JRtUaKZ3m1YjUw9f0qg6L1T9395vc/RvrPba6EPaCu19bLeuAuz8wa1+yqtRAuSfj04ThGFStERGZnpltNLNdrFVZrnX3Sa8F+WTgB2b2KTO7ycwuNrOZd1pJQk0O1RqZT65DUCLSTf/66Fel7kJ0J/Su41j9O7//j+5+yN23AycDzzaz0ydc7gLw08BbgGcBTwR+fdZOLsz6wJROPPoAe+8fXQLZuuVB9i/Pt1M9tPkQGw+sHxYPHnOYhfvqy4UrW2BpedL7GkvLXlvbw6xuXWRx/2qjbUS04fjjNDQhInk6eAjbN+GOZDp3u/uOcXdy9x+Y2d8AZwO3TbDcPcBN7v51ADP7C+BM1i7DNLVkw0/jqjVNa7NaE+F6UJowLCIiTTCzx5jZsdXPm1ibHHzHhA//B2CbmT2m+v0s4Muz9iXsnJpcLp2gCcNl0tyablBFTqQWjwO+YGa3sBZSrnX3T5vZb5rZHtaGpG4xs0sBzGxH72d3P8Ta0NPnzexWwIA/mbUjSYeffmzTPXzzweOTtb+4eYXVA0uttLW6GRYnOCK9ySGoh7Zu5Ij99V3jqm2+bUtTZVUREZmRu98CPGPI7X8E/NGQ23cCr+37/VrgjDr6ErZSA6rWlKDOIai2qVojIpKX5KEm9dyaNkU4vHsWXR2CEhGRvCQPNeM0Xa3p2uHdmjAsIiKlChFqoldrNATV3WqNhqDKp8nCIuUIEWrGmbdaM07Ew7t1huHRVK0REZFhwoSaJqs1mjCcniYMi4hI01oNNUdumP3stCVdOiFCtWYWXR2CEpFu0KUS8hemUgPzVWt0ePd0NGF4OqrWiIjE13qo+fEjvjfzY5uu1sh4qtaIiEhUoSo1kLZa07UhKE0Yno6qNSIisSUJNTlXa0oagkoh5wnDUi4d1h2LvkDIrMJVakDVmkHTVGvaoCEoERGJKFmomadak1rkak0OQ1CaMCwiIk0IWakZp6ST8YmIiEg9koaa9ao10U/GV5dcJwxHGoJStUZERCDTSg1ownDuNGFYRCLSCfjyljzUNFWt0YTh5kWq1rRN1RoRkXiSh5p5qFozWg4Thuuki1yKiMhC6g7AWrXmaw89dujffmzTPXzzweNnWu6JRx9g7/2jyyBbtzzI/uXRSVnDngAAB0xJREFUwyCLm1dYPbA0U9sy3uq2TSzuizO/aVobjj9O5zcpyOF77lUFTvJ36GCnP5eyrtRAHtWaSUSYMDyLSENQqtaIiHRbFqGmybk1bchlwnDu56xJQd/sRUTiaDXULNnBkX9LdTI+TRieX6RqjYiIdFcWlRpYv1qTwxCUJgy3I8UQlKo1IiIxtB5qnrj0/ZF/a6pak9Ph3V2U+xCUiNQv5ZcFnasmX9lUasbJoVozCQ1BzU/VGhGRbkoSamat1qScMFxHtUZDUKOpWiMiIvMqplID81VrcrwelMSiao2ISFrJQk2O1Zpx2p4wrHPWjKZz1oiIdE9RlZp5acLwwzQENRtVa0RE0kkaatar1qwn98O7J6EJw/NTtUZEpFvCVmpKPrxbE4ZHK6FaIyL502HdeQobauaRuloj3aYhqHx1+UKAIiVIHmo0YXg0DUHNT0NQIvnSFwSZVvJQ05QmD+/OecKwhqDaoQ9jEZH2hQg1Kao1udA5a+anao2ISDeECDUp5DIENSmdsyYeVWtERNrVaqg50nzk35qo1mgIqh4aghIRkRzMFWrM7Gwz+0cz+5qZvb2uTrUlQrVmErlOGI4k1RCUqjUi+dJh3ZMblwfM7Agz++/V328ws1Ob6MfMocbMNgL/D/BzwNOAc83saeMe96SF0Tv6WU/Gt57o1ZqSz1mjISgRkfJNmAdeA+xz9x8H3gf8QRN9madS82zga+7+dXdfAT4BvLSebj1a1ycMd1HdQ1Cq1ojkR++fLEySB14K/Gn185XA882s9gme84Sak4Bv9/2+p7otKxqCEhERmcskeeBH93H3g8B+4Pi6O7Iwx2OHJaxHzQQ2s/OB86tfD2w7ac8/rr/Yb87RpaFOAO6ue6EB6XmWRc8zlQcaWWq859mM4p6n2Z8Nu/m0tvsxqfsO33v1NQ/82QkNLPpIM9vZ9/sl7n5J9fMkeWCizDCveULNHuCUvt9PBu4cvFP1pC8ZvL0tZrbT3Xekar8tep5l0fMsi55nWQZ27qG4+9kJmp0kD/Tus8fMFoCtQO3XJZln+OkfgCeZ2RPMbAl4OfBX9XRLREREMjFJHvgr4Lzq518G/j93j1OpcfeDZvZ64GpgI/Bhd7+9tp6JiIhIeKPygJn9V2Cnu/8VcBnwZ2b2NdYqNC9voi/zDD/h7p8FPltTX5qSbOirZXqeZdHzLIueZ1m68jwnNiwPuPs7+n7+IfArTffDGqj+iIiIiLSus9d+EhERkbJ0KtSY2VvMzM2sicPdkjOzi83sDjO7xcyuMrNjU/epLrlfkmMSZnaKmX3BzHab2e1mdkHqPjXJzDaa2U1m9unUfWmKmR1rZldW78vdZvZTqfvUBDN7U7XN3mZmV5jZkan7VBcz+7CZ7TWz2/puO87MrjWzr1b/b0vZR3lYZ0KNmZ0CvBD4Vuq+NOha4HR3PwP4CnBh4v7UYtZLcmToIPBmd38qcCbwukKfZ88FwO7UnWjYB4C/dvenAE+nwOdrZicBvwnscPfTWZso2sgk0EQ+CgweJv124PPu/iTg89XvEkBnQg1r15p4Kw2c7CcKd7+mOlMjwBdZO1dACVq9JEcq7n6Xu99Y/bzM2g4wu7N0T8LMTgZ+Hrg0dV+aYmbHAM9l7agP3H3F3We/GF1sC8Cm6vwjRzHknGW5cvfrePT5VPpP+f+nwC+02ikZqROhxszOAb7j7jen7kuLXg18LnUnalLEJTmmUV3B9hnADWl70pj3s/YlI80VXdvxROD7wEeqYbZLzezo1J2qm7t/B3gPa1Xwu4D97n5N2l417rHufhesfRkBTkzcH6kUE2rM7P+txnMH/70UuAh4x7hl5GDM8+zd5yLWhjIuT9fTWrVyeu0ozGwz8Engje5+X+r+1M3MXgLsdfcvpe5LwxaAZwJ/7O7PAO6nwGGKaj7JS4EnAI8HjjazX0vbK+mquc5TE4m7v2DY7Wb2k6y92W6uLgh6MnCjmT3b3b/bYhdrMep59pjZecBLgOc3cbbGRCa6JEcJzGyRtUBzubt/KnV/GvIc4BwzezFwJHCMmX3c3UvbEe4B9rh7r9p2JQWGGuAFwD+5+/cBzOxTwL8CPp60V836npk9zt3vMrPHAXtTd0jWFFOpGcXdb3X3E939VHc/lbUPmmfmGGjGMbOzgbcB57h7M5flS6MTl+SwtdR9GbDb3d+buj9NcfcL3f3k6v34ctZOl15aoKH6jPm2mfUufvh84MsJu9SUbwFnmtlR1Tb8fAqcED2g/5T/5wF/mbAv0qeYSo0A8EHgCODaqir1RXf/92m7NL8OXZLjOcArgVvNbFd123+qztQpeXoDcHkVxr8O/LvE/amdu99gZlcCN7I27H0TBZ1x18yuAH4WOMHM9gDvBN4F/A8zew1roa7xM+XKZHRGYRERESlC8cNPIiIi0g0KNSIiIlIEhRoREREpgkKNiIiIFEGhRkRERIqgUCMiIiJFUKgRERGRIijUiIiISBH+f4v36koxvUtpAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1008x576 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "m2 = GPy.models.GPRegression(\n",
    "    X, y,\n",
    "    kernel = GPy.kern.RBF(2, ARD=False)\n",
    ")\n",
    "\n",
    "m2.optimize()\n",
    "\n",
    "Xnew = np.vstack((Xi.ravel(), Xj.ravel())).T\n",
    "\n",
    "mu, Sig = m2.predict(Xnew, include_likelihood=False)\n",
    "\n",
    "\n",
    "plt.figure(figsize=(14, 8))\n",
    "\n",
    "plt.contourf(Xi, Xj, np.reshape(mu,(50,50)), levels=np.linspace(0, 300, 20))\n",
    "plt.axis('square'), plt.colorbar();\n",
    "\n",
    "display(m2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis of Uncertainty\n",
    "\n",
    "Let X be a random variable defined over the real numbers, $\\mathbb{R}$, and let $f : \\mathbb{R} \\to \\mathbb{R}$ be a mapping between real numbers. Uncertainty analysis is the study of the distribution of the random variable $f(X)$.\n",
    "\n",
    "We will assume in this example that $f$ is the Branin-Hoo function defined above and that $X$ is a random variable with *uniform* distribution over the input space of $f$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computing $\\mathbb{E}[f(X)]$\n",
    "\n",
    "The expectation of $f(X)$ is given by $\\int_Xf(x)\\mathrm{d}x$. The basic approach to approximating the integral is to compute the mean of our 20 observations and compare it to the true value, $\\mathbb{E}[f(X)] = 54.31$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "67.43511010158977"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(y) # Mean of our samples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, we can do better by instead using Monte Carlo sampling of a Gaussian process. We define $10^5$ Monte Carlo samples over our function space:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the function space\n",
    "Ximc, Xjmc = np.meshgrid(np.linspace(-5., 10., 100), np.linspace(0., 15., 100))\n",
    "Xmc = np.vstack((Ximc.ravel(), Xjmc.ravel())).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2\n",
    "(a) Use your GP fit from the previous exercise to calculate the mean of the samples evaluated over your model of the Branin function. Discuss how it compares to the sample mean and true expectation of the Branin function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "57.91347373791677"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mu, Cov = m1.predict(Xmc, include_likelihood=True)\n",
    "np.mean(mu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(b) Can you utilise the prediction variance to create a confidence interval for the prediction of the mean? Does the true value lie within this bound?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "56.97838025272408\n",
      "58.848567223109455\n"
     ]
    }
   ],
   "source": [
    "lb, ub = mu - 1.96*np.sqrt(np.diag(Cov)), mu + 1.96*np.sqrt(np.diag(Cov))\n",
    "\n",
    "print(np.mean(lb))\n",
    "print(np.mean(ub))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computing $P(f(X) > 200)$\n",
    "\n",
    "In many cases, we might be interested in looking at the probability that $f(x)$ will take the value above a given threshold. For example, in the case that $f$ is the response of some physical model that represents the maximum constraint in a structure, it may be necessary to ascertain the probability that the constraint exceeds some maximum value."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 3\n",
    "(a) Use the original 20 observations to calculate an approximate estimate of the probabilty that $f(X) > 200$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.05"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum(y > 200) / len(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(b) Compute the probability that the best predictor of your GP is greater than the threshold ($200$), with confidence intervals. How does your estimation compare with the true value: $P(f(X) > 200) = 1.23\\times 10^{−2}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.016\n",
      "0.0167\n",
      "0.0173\n"
     ]
    }
   ],
   "source": [
    "print(np.sum(lb > 200) / len(mu))\n",
    "print(np.sum(mu > 200) / len(mu))\n",
    "print(np.sum(ub > 200) / len(mu))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(c) Can you design and implement a procedure that updates the GP model sequentially with new points in order to improve the estimation of $P(f(X) > 200)$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Review Bayesian optimisation notebook -- exercise left to the reader."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Credit\n",
    "\n",
    "This notebook was written by Wil Ward. It is adapted from notebooks by [Rich Wilkinson](https://rich-d-wilkinson.github.io/) and [Neil Lawrence](http://inverseprobability.com/)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}