Revision 1db48f3a735eb0fba06a7d503f080a7ead512604 authored by Artem Artemev on 11 July 2018, 12:50:44 UTC, committed by GitHub on 11 July 2018, 12:50:44 UTC
1 parent 707b195
coreg_demo.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import gpflow\n",
"import numpy as np\n",
"import matplotlib\n",
"%matplotlib inline\n",
"matplotlib.rcParams['figure.figsize'] = (12, 6)\n",
"plt = matplotlib.pyplot\n",
"from gpflow.test_util import notebook_niter"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A Simple Demonstration of Coregionalization\n",
"--\n",
"\n",
"*James Hensman 2016, 2017*\n",
"\n",
"In this notebook I'll demonstrate how to fit a simple model with two correlated outputs. For a little added complexity, the noise on the observations will have a heavy tail, so that there are outliers in the data. \n",
"\n",
"In GPflow, multiple output models are specified by adding an extra _input_ dimension. We'll augment the training data X with an extra column contining 1 or 0 to indicate which output an observation is associated with. This also works at prediction time: you have to specify two columns containing the location and the output of the prediction. "
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x7fb22490c5f8>]"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X18lOWZ6PHflZeJyQwgZCKliEKqqCjbFXM00ZXa9QWMnqJbTcuyQN0KpaX7OcXderqnW92627Nbt8V2VykFbQWWpY3u+rKKsG2twqmJGrFbFOWlAUWqkAkpMpMwCcl9/njmeeYlM8kkmfe5vp9PPk5mnsw8E+L13HPd133dYoxBKaVUcSnJ9gkopZTKPA3+SilVhDT4K6VUEdLgr5RSRUiDv1JKFSEN/kopVYQ0+CulVBHS4K+UUkVIg79SShWhsmyfQCJer9dMnz4926ehlFJ55bXXXvMZY2qGOy5ng//06dNpa2vL9mkopVReEZF3kjlO0z5KKVWENPgrpVQR0uCvlFJFSIO/UkoVIQ3+sV5ZD/6O8Pf+Dus+pZQqIDlb7ZMVr6yHrX8Frz4MS5+x7ttwM3S8bd2+fFn2zk0ppVJIg3+kWbdYgb/jbVhTb93X7YOaC63HlFKqQGjaJ5KnxhrxV3mtoN/ts24vfcZ6TCmlCoQGf5U8nQ9RqmBo8I/k77By/PaI3/4EsOHm6KBXjOz5EPt3Yf+utv6VXgCUykMa/CPtedLK99dcCF9qtb5qLrTu2/Nkts8uu2bdEv5drKm3vuzflc6HKJV3dMI3kl3NM+uWcI5/6TNW4C/2Sh97PmRNvfVpCHQ+RKk8psE/VmyQ99Ro4FdKFRxN+6jk6HyIUgVFg79Kjs6HKFVQNO2jkqPzIUoVFB35p1oh18Jfvix6clfnQ5TKWzryTyXtDaSUyhMa/FNJewMppfKEpn1SaSS9gQo5PaSUynk68s8GTQ8ppbJMR/6plGwtvLZKUEplmQb/VEq2Fl5bRyulskzTPqmktfBKqTyhI/9US6YWXlslKKWyTIN/NmirBKVUlmnaJxs0PaSUyjIN/tmiraOVUlmkaR+llCpCGvxVlI0th/D5g873Pn+QjS2HsnU6Sqk00eCfT9LcEmJjyyHueepNFq5rxecP4vMHWbiulXueelMvAEoVGA3++cJuCWGXg9rlolv/CnauDh8XuiCMZgTfOHsK55/lYf8xP/Me2MG8B3aw/5if88/y0Dh7SlrellIqO8QYk+1ziKuurs60tbVl+zRyhx3sO9621gVAeCP16vPgjm3W7dAx3+j7HNuq/ieb77yCSW4XC9e1sv+Yn5tmf4SHFl2W8GV8/iDzHthBZ6DXemq3i+2r5uL1VKTz3SmlUkREXjPG1A13nFb75Is9T8JtG2DDTeGgX14F4z8KnQeiWkifrr6Al0/MpeNkkBu/t5PxlWV0dffhKi3h2d0fcEXLIZY0TM/WO1FK5YCUpH1E5EcickxE3kjwuIjIP4vIARH5jYjMScXrFg075dO8CMxA+P6+bmutQEyPoLI7nuXBZTdQKkK/MXR19yFAb//A4BROxDyCzx/ki2u30XjqGardLqrdLjoDvdz4/Z3sO3rS+RGdBFYq/6Uq5/8oMH+Ix28Ezg99LQd+kKLXLQ6zbrFSO50HoOc4IOHH3ng8+oIQwVUWPs4AE6vK2bK8PpzCiZlHeL7tDb714V/zd+WP8uI1+9m+ai414yroOBlk0fpW1rxwgH1HTzqTwGteOOBcBLRKSKn8krKcv4hMB54xxlwS57EfAi8YY7aEvt8LXGOMeT/R82nOP8bO1fCLb4a/r5wEZ4yHrkMA+Mx4SkuEieYEwYkz+cSxv+SD/nFRT1EqwnNfuZqZk0P3R8wjDFR6KRGg20dw4kye+vhamq65jH1HT7JofSsd/l7nOfqNodbrxmA46Ovm4o+O483fneT8szxsWV7P8UCv8zP3LbhYU0xKZVCyOf9MVftMBQ5HfP9e6L78ke2dty5dDC53+HspgUtuA+B09QV8cdyDXNfzj/yWs6no2sf1tDifDyZWlTtBe9H61vAI3VND88Vr8JnxlPRYaaOBSi9/1vc33L3tAza2HGLm5HE895W5TKwqB6DfGATo6u7loK+bGo+LN393EldpCfuP+bl+9Yvc+L2ddPh7qfG4tEpIqRyVU6WeIrJcRNpEpK2jI4e6Ww5VZpmJC4D9er2B6C6gbz8D195L2R3P8oMV8znlmkTTqa/zjb7P8RPmYYCacRVsu2ovd111JjUeFx3+Xp5ve8M572svmkxpSTg99PueXtp9Ac4/y4M/eNq5UJRIdAqpq7uPareLzcvqOf8sD739A6GLQh/9xlAqwuZl9VolpFSOylS1zxFgWsT3Z4fui2KMWQesAyvtk5lTS0I6NmZ/ZX10Yzd/R7ijZ+z9W78a7gK69Bl4fRP8erN1X4WVwtnz1HcJ9H6cABPY1H8DJWKYWFnG/O7/ZPLOR7l+YCpVf/QoADe8eicE2tm1/13+4Ph2yswJjjOeAQNe+ZDmM77Fjkt+zDe37eWxtsMIQmeglxKBgYh/lQFjmOR2sWV5PdevfpGu7j7nsfGVZUxyuxK+/Y0th2icPcW5OPj8Qbbufl9TREplSKZG/k8DS0JVP/XAiaHy/Tkn1TtvDfVJIt79e56AWbeGO3/auf9r74VZt3D6xzcxd/+3WVz6X85LDBjo6jnN1v4r2DcwlZklR1i867Pc0nIbZwbaed91Lv+15yhlnXsJTpzJbbKaecFvs29gKh/jPRpLX+b8szwc9HXT7gtEBX77g0JXdx8L17Xy2w4/H/acdl5bIh6LnAS26UpipbIvJRO+IrIFuAbwAkeBe4FyAGPMWhER4EGsiqBu4A5jzJCzuTk34evvsEb9do19ldfqwz+a4J9owVb1edbtzgPR99sjfk9Nwp/dNzCVr3n+LwdPVUWNwAGqOcHPzvjfTOJDALpkAseXvMiKJ96lvvM/2D5QzzEzHoCz5EPmlbTSWv0nPLRoDgvXtToLvgBqvW7WLr6Mn791lCd2HWH/MT+eijL8wdOUikStKejtH4g74WsH+/3H/FSHPh10BnqdCWNNFSk1ehmd8DXGLDTGTDHGlBtjzjbGPGKMWWuMWRt63BhjVhpjPmaMmT1c4E+b0U7apmLnrcjX9tRYC7Zc7uhPEndss76G+ISxcXeAztv+nR7XROeYHtdEfvVHj9Lw8QujcvM2ITpdM+GMcj5WYwXa/yi90Qn8E6vK+dev3Exr9Z+w/5ifn791NOp53K5S1i6+jJmTx/Gla85jy/J6bpr9EfzB09SMq+C5r1zNz+76hDMHcN1FZzmBP7L00+upYMvyemcdQWegl+pQ+kgDv1KZUTwrfO1Uy6sPW8EUwiNoGLqXfuTOW7E/m8wGLPZr7/gnWPw0uKutBVu9gaRPf2PLIfzB09y/bS//6jnFv/X1UxmK84FgP5tb3+VAdyUQLscEqOEEm11/j1c+tKp6BCb1WBcuue3fOaO8lEBvP2BN6to5/Oa2wzyx64gTmMEana/cvMsJ0l5PBQ8tuowrWqLz97fOmcr92/byTme3k/axR/qA5vWVygE5Ve2TEolG97NuCW+VuKbe+rIDerxJ28jnuXwZzLwRxoeqU+2R+8zG5M5p1i3gmQz+o7D2Sniwzkrt2OxPAD+eb33FfML4/Q9u4PtPvcRjbYeZM6mPB3u/4QRznxmPVz5kzel7qOYE1e5ypw7/nEmVzC99mZklR9g3MJVbzXe4/tS3OVQyDTre5t9+/M9OcLdH4QvXtVpvsaLMaeq2fdVctq+a6zR927o7erpmScP0qBF7U920IRvE2WmfeK8db45AKZV6hTXyH250v/SZwXn7eJO2sc/z+ibY95z12A8aYMkz1si98wDs22rdP9To31NjjfjXXgmmH079PvxY9XnQtBkeXxo+z5hPGGd2vM3nzvw13/VN4IbKHU4wX9j7NwBscf09M0uO8MDsQ8xa8Gm27n6fxtlTaG47zP3bbsDrdjG+7naeuPLjLFzXyqeP/R/unvYW3z1c7+TZITw6j6y6iRzRb1len1RFjp3WiW0QZ39i2NhyiP3H/NR63TSvaACgaW3LoNdWSqVPYQX/VJVkxj5PZPuEQAf8oB6r2h0reCfz3O5qOGNCqD1DyBlnWjn+UDXRy88+wsUfHY/n0tvBU4PPH+T5i9fQVLWLhbOW8OgDO/hh4Fq6S/vZ2n8FnUwAYGHv33B7ZRt3LvhLvJ4KJ3h+6Zrz8FSU0Tj7ukEB/DMNf0qw5dCQwT02CEc+dyoYTNzbSqn0K7yWzomqciB+hU1kJc1Qz1M5yboIRI7aKyfBylfif3KIrNU/+jZs+pSV9kFwLhxSCitegskXOuWP8Ubi9y24mMbZUwbV0tvsHH+uVMsMV80DaLWPUmmSa+0dsi9y0vZLrdaXPQew58nBcwWBTujrGfnrxKvhtwN/SSlgrIuGlFopoE2fAn/HkBup1NdWs3BdqxP47YKe8hLh3ElV9BtDjccVNx+fDVt3vz/kfIFW+yiVfYWV9oktyQTr9tqrrJx743fCKZo9T4YXTUF0jj/QCeuuhv5eKHdDWUV0usYevfcctyZo7dQNJE49ucZB78nwJ41AZ/iisOdJvJcvS5gnjwymt86ZynUXTWbRwy/TcTLIZy6fFkrtTMmZfPlY5guUUplRWGkfe9QdOWG69iorwHomw4pfWffZ6Z/G71gTtbELp073WGWYpS5YvtOa1LVX1bproid8Ifw8tkSppz1PDmrd8PKzj/Cxm1bh9VTg8wejUjuRu2ilvB1CovYSw5WtpoAu8lIqfZJN+xRW8Ieh8+1D5fpjA7bLDZ9/HiZfaH2/czUcfgU+9S/hlbavb7J668QGzNjnKnfDneHnan7hNRpLX+Y/yhq556k3meGt4qbZH2Xr7vdp91m1/26XVX+floAY7yIZe0FMo6HmN26a/RG+ueAS7fmj1CgVb/CPJ5nWDKlq3xD7KaKvB/rCnyKa3wrwh88vZmbJEY5e/fd89vVLOOjrjnqKGd4qfri4jpWbdzkTvikNfonaSySa/E6DeJ9k7n3qDZ7d/QE1Hqtb6Ej2HlZKWXTCdyRS0b7BFjuxfOfzVuDv74VH/pjbXm5y6vQX/WoKRz+MXtQ0obKM2+umMXPyOLYsrx8+8I+mZUWyjerSuIdB7MIwr6eCby64xGk7feP3dnL96hfZf8zv7D2sTd+USp3CD/7JBPbhKoFG4vJlVurEDqSTL7TmDVxu6A1Q0uPDZ8bzZ33f4EB3Jd2h1grO6Z7q5/5te9nYcmj4uvp07jOQhT0MvJ4KNi+rj9p7GKL3HtbtIZVKjcKq9oknmb48do47cq7ArgQaTf479DNOasNdDWWVTi+f0pJw751IgrVTlqu0hPra6uFfZ7SL2hJVRW24OXzRSsceBkmY5HY5nUFtJQIPLZoDaI8gpVKlOHL+WahssSc1/4f3NFsqvkVZ5166ZAL9AwavfOi0Z7BX6drsUW/Sef7RzFUkO+E71IK5NPz+IquAIpbCATChspyyEtGqIKWGkWzOv/BH/jA4SHlq0l7R0jh7Cpta3uHCzv+grHwvv+Vsmnq+TrXbxYN99zCz5AiNpS+zqf+GqJ8bX1nGsrm16R3VJvqks/Wr0aP6QCcEPwx/bwas+yLLXFP4e7z3qTecHH9v/wATKsv4sOc0BjjREy5/1cCv1NgVfs4/C+yc9Jbl9Ww942a+0fc5mk59Hdw1PLh8Hl8uv49v9H2OZypuojS0XFew+ul3dffxxK4jyXW3HMtE9eXLoj8d7HnS2jHM/tmjb4cXutl6jlvN6ezAHzw59GuMYMJ4Y8shnt39AZ6KMnr7B5jhraLaXZGw44/m/pUaGw3+KRa5ReHxQC8DxrCp/wY6mcCpvn5+/tZR9gUqaa3+E5bNrXVy/AZYNrc2YdvkuFI5UR3b8vrhP7YCf6kLJpwTPs6EJqirz4NLFyd+vhFOGNvtLfzB07hdpfy+u492XwBXafhPtFSstE/T2hZuX/uSbvuo1BgUR84/gyLz1pGbqkQ2X7t1zlSa6qY5K3fra6tpbe9kScP0kS9qSuV8RqLFaRDdyRSBL7aGF8Aleq4RriXw+YNR7S3shW61XjcGw0Ffd9TvVHP/Sg2mdf5Z4vVUcOucqVSVlzhBakJlGV+85mPOqN5TUeYErCUN05k5eZwT7EfcNjk2fZPMfEay6Zhya2cwmhcRPf1qrPuGSi2lYNP7M8pLuXv+BTSvaOCxFVdS7XY5v1PN/Ss1Nhr8U2xjyyHu37aX4OlwsPSf6ufBXx7g1jlTU79ad6SGSsesvSr+/IGd46+cZH2Bdd/rm1J2Wol293pi15GUvYZSKqw4qn0yqL622qlWsbdSt/P61100mZmTx2X1/BLW79vbTMYr/wQrx3/HNuv2j+dbwb9iiPeSzFqCCJGdS2P7/STaT3jhulYd/Ss1Shr8U6y1vZPe/oFB+f7e/gFa2zuzH/ztdExs/f4VK6zbly4ON677g8+EA3zkvMId24afVxjhpvdDtYEGEl4YtOmbUqOjE75psOaFA6zf0e6sUp1YVc6yubU01U3LjWAVr4NpbyD1XT5TOBmd8pbWShUo7eqZJZHVPhCu8omsWMlq3j9RFY7dfC5LXT6VUqmh1T5ZYueua71uZnitLRZLRWj3BTjo63YalGVNorUB9q5l8SpzhqsOSmP3T6VUemjOP8Uic9dA3G0ZszpBmai1w+uboOUha++BSHZ1kL3FJURPBMPQj2dgZzCl1Mhp8E8D+wKQVIuGbIgXkH/z0/iVObdtiF75az8W2d0zC90/lVJjo2mfNElUt75wXWvuXRSGahPxzs6hF2ulYDGXUirzdOSfJkPVredclcpw+xmMdDczpVTO02qfJIy2zLAgyhOH69EDWd8PWCkVptU+KRLZpdPnDzrpnGQ6SsbbpzavAj9Ep4QaVsLSZ8Mpodc3WXsApKqzaBpsbDkUlWbTVtBKWTTtMwx7U5b9x/zMe2AHgLObVFZLNjPFTgkFT8IvvmlNDN+2AfZttW53vA2zboXGf0rNFpgpZF+4N7W8Myj1BroNpCpuGvyHYXfpXL+j3SnZnFhVzq1zphZPTxk7728H+w03Wffb6Z3IwA8Z2SktGUV/4VZqCClJ+4jIfBHZKyIHRORrcR7/nIh0iMivQ193puJ1M8Hu0vlhz2nnvg97TnP/tr3FlT4YS1VPlhaBeT0VbFle71Ra2ZVX9qeAovr3UyrGmEf+IlIKPARcD7wHvCoiTxtj9sQc+lNjzJfH+nqZNlSXzvra6qyeW14YbpFYBj4hnOrrd24PGMPxQC8rN+/S9I8qaqkY+V8OHDDGtBtjeoGfAAtS8Lw5IbJLp8Ha0iSyS2fRGO1+wbHbQ66pD08Q24vA0vTJwJ6cD/T2O3sld3X3ceP3djpluJr+UcUqFcF/KnA44vv3QvfF+rSI/EZEHheRafGeSESWi0ibiLR1dORGbfmShuncPf8CxleGPySNryzj7vkXFNeIcbT7BQ+XLhrhXr8jEbnW4rmvXM3EqnLA+uTmdpVmv9WGUlmUqQnf/wS2GGOCIvIFYAPwx7EHGWPWAevAqvPP0LkNyecP8sSuI3R190VtJPLEriPOPrxFYbiFYEN5fROYgfD3ZsC67+q7Em8uk4L2ELF9lkpEnMfOKC8d03Mrle9SMfI/AkSO5M8O3ecwxnQaY+xi64eBy1LwuhkROXrcvmou21fNdfbitTcaKRqj2S9452qrRLTneHgbyJ7j1n07V6e9PYR9AcibVhtKZUgqRv6vAueLyAysoP9Z4E8jDxCRKcYYO1J+CngrBa+bEUPtMFVUaZ88lletNpTKkJS0dxCRRuB7QCnwI2PMt0TkPqDNGPO0iPwDVtA/DRwHvmiMeTvxM+ZWewc1RjtXw0v/Yo34wRr9X/kXVtpnuPYRKWoPkajVBpD/LTiUipBse4eU5PyNMVuBrTH33RNx+6+Bv07Fa6k8dOlia68Am5RY98GI9/odrdhgbgd7XQGsipWu8FXpFVsiCuES0aXPjG0ieYx0BbAqZtrYTaVXMiWio5lIToGhVgAXTRWXKlo68lfplcWRvVIqMR35q/TL0sh+OHm125pSKaYjf6xKEH/wtLNoy+cP0tx2GE9FmU76FbDIEtBb50wlEDzNc7s/YP8xP81t1qJ1/RtQharog7/d8x3g8bb3WLv4Mr6wqY2Dvm7nGP2fP4+8sj46xeTvSJhisv9d/UGrSyvAuZOq+PInz+OxtsP6N6AKWtGnfRpnT2GGtwqAdl+AeQ/scP6nr/W6teojn4yiT9CShuk01U1z/gbeOd7NQ788oH8DquAVffD3eip4bMWVTtMve8nbxKpymlc0aNVHPkmmg2gc+jegilHRB39VQNLcJ0ipQlL0wd/nD3L72pfo6u4DcDZs6eruo2lti1Z9FIGh/gY+veYl529AN39XhaTog//W3e9H5Xe3r5obNQdQdJ0789koN5yJ9zdQHUoBvXO8m+a2w05Z6D1PvakXAFUQir7aJ7Liwy71fGzFlVrqmY9G2Sco3t/Ali80OGsAHtl5kEd2HtTWD6qgpKSrZzpoV081KiMo9RyOzx9k3gM76Az0AlDtdrF91VydAFY5LdmunkWf9lEFJgOriTX3rwpB0ad9lIonsvUDQKkInYFemta2YDDOHIGmBVW+0uCvVBx264dar9sJ9qUitPsCAJr7V3lP0z4qP72yPrqCx9+RcBUvWG08Ist2h0vdLGmYzn0LLqZ5RQOPrbiSareL/tD8mN32eevu90f0nErlEh35q/xjt3F49eHBVT0wKMdv928a6Y5d9v2xaz0GjKG57TD3b9vLj391kNvrptFUN013AVN5Rat9VP4Z4b6/dv5+/zE/1W4XEN6xa7iNWyJ/1u0qJdDbD1gN4AyGd4/3ADiPJfOcSqWTVvuowjXCNg5j2bErsu3zEyuvotbrBqzFXz5/r3NcoLdfdwFTeUXTPkoNwU7fNM6egtdTQfOKBq5f/SJd3X109/YjhBvBKZVPdOSv8o+/A9ZeNbiNw9qr4rZxGOuOXUsapkeN5ktEnNsGq/un7gKm8o0Gf5V/tn4V/Eeh1AVLn7W+Skqt+7Z+1TomovonMnWzfdVctq+ay/lnedh/zD+i3k2RFxG3q9S5f2KVle4ZzXMqlS2a9lH5p/Gf4N2XrGC/4SY43QMD/dYF4BNfi54QBpY0WNU/duoGcEo1R1KVE3kR2bK8nua2wzze9h7tvgCt7Z2jek6lskWrfVR+8ndYm7V0+6zvpRRMf1LVP2OxseVQ1EXE5w9qwFc5JdlqHx35q8JwxgTrv/bFIE2buMQGea+nQgO/yksFlfMf6SpOlad2roYfzw9P+FZOgp7jcOpEts9MqbxRMMHfXsVpV1vo5hsF6pX18ItvQucBqD7Pmuw9Y7z1mOmHcnfSm7goVcwKJvg3zp5CzbgK9h/zM++BHcx7YAf7j/mpGVehDbgKib1JO0D3cWvCt+uQ9X3lJPhf/21dENze8CYuSqlBCib4b939Ph0ng5QIzirOEoGOk0EtvctnsQ3cAP7gM9bovud4OMcP1ieAQCc0L4KAD6VUYgUz4esPngZgIKJ4yb5tP6byzFAN3Mrd4ePsSp+uQ/CDepw1t9XnWZ8UlFKDFMzI/7qLJlNeKoPuLy8VrrtochbOSI2ZneLpeNsq61xTb90udUFfILy61/RbFwDACfyVk+CObSmv9lGqUKQk+IvIfBHZKyIHRORrcR6vEJGfhh5/WUSmp+J1I7W2d9LXP3jNQl+/obW9M9UvpzIhXgM3lxv6e62Lwpdara/q86wLgFIqaWMO/iJSCjwE3AjMAhaKyKyYwz4PdBljzgMeAL491teN1Th7CtXu8kH3V7vLdcK3kJRVwrX3RtfwD0QG/tCnv57jVjmoVvsoFVcqRv6XAweMMe3GmF7gJ8CCmGMWABtCtx8HrhWRwTmaMWhuO0xnoA+wGm1NrLIuBJ2BPprbDqfypVSm2G0aYhu4/ean4WP2PAldB63b1efBF0OfBMAqB81wtY+uNVH5IhUTvlOByOj6HnBFomOMMadF5ARQDaSsJMNTYb2VWq+b5hUNADStbaHdF3AeU3lmz5NWjt9u0wDhCd89T1o7dtm7dgVPwqWLrU8Dd2yD1zdBxbhBu3ql02h3DFMqG3IqKorIcmA5wDnnnDOin43tuw7QvKJB+67kMztwz7olnOJZ+kw48MceZ/PUwNV3ZeYcIzTOnsKmlnectSYQ3jFMU48q14y5sZuINAB/a4yZF/r+rwGMMf8Qccz20DEtIlIGfADUmCFeXBu7qXzk8weZ98AOOgPWLl/VbhfbV83V3b1UxmRyG8dXgfNFZIaIuIDPAk/HHPM0sDR0+zbg+aECv1KFRucCVK4Zc9onlMP/MrAdKAV+ZIx5U0TuA9qMMU8DjwCbROQAcBzrAqFUQYndMQystM+N399Jx8mgzgWonJKSnL8xZiuwNea+eyJunwJuT8VrKZWrYjd7gXCQr/G4dC5A5ZScmvBVKp/FKzqwd/dqnD1l0FzAluX1OhegsqZg2jsolQtiN3vXzV5UrtLgr1Saxc4FVLtddAZ6nb0nlMoGDf5KpVnkXMD2VXPZvmou55/lYf8xv7YbV1mjOX+l0myouQBNCals0eCvVAboxu8q12jaRxWP2F3B/B3WfUoVIR35q+Iw1K5gkNEGcErlAh35q+KQaFewmgtzZqtHbQGhMkmDvyoO8XYFq/JGbwqTRXY7aLv80y4PveepN/UCoNJC0z5K5QBtB60yTUf+qjgk2hVsw805sdWj11PBluX1zgIwe0GYtoBQ6aLBXxWm2Mqe1zeFc/z2xu/2HECGt3pUKhdo8FeFx67ssUf1/o7wvr9/8Bkrx2/PATR+JycqfbQFhMo0Df6q8AxV2XPp4vBxnpqcCPygLSBU5umEryo89qh+Tb2V14ecquyJR1tAqEzT4K9UjtAWECqTNO2jCk+OV/aMhi4AU6mmwV8Vnj1PFlRljy4AU+mgaR9VeOxJ3Fm3hHP8S5+xAn+OTPCOhC4AU+kgxphsn0NcdXVDd+SDAAASNUlEQVR1pq2tLdunoVRO8PmDg/YA3r5qri4AU4OIyGvGmLrhjtO0j1JKFSEN/krlOF0AptJBg79SOU4XgKl00AlfpXKcLgBT6aDBX6k8oAvAVKpp2kcppYqQBn+llMoBmV7FrWkfpZTKMnsV96aWd9iyvB6Aheta2X/MDwxO+6WCBn+llMqybKzi1rSPUmmQjo/wiZ5Tm77lv2xs46kjf6VSLB0f4Yd7zkymC1Rh0N4+SqWYvSJ3/zE/1W4XEP4IP9qRXKLnrPW6MRgO+rpT9loq81L5N5OR3j4iMklEfiYi+0P/nZjguH4R+XXo6+mxvKZSuS4dH+ETPWfzigYeW3FlRtMFKjXsdN3GlkM0tx1m/zE/NeOsf+fPXz0j7au4x5r2+RrwC2PMP4rI10Lf/+84x/UYY/5wjK+llFIFwU7j/csv9tPhtz7BVbvL6TgZZNHDL9NxMsjd8y/AU1GWttTdmNI+IrIXuMYY876ITAFeMMZcEOc4vzHGM5Ln1rSPylea9lHDifz3LBWhPxSH7dtj+ffLVEvnycYY+zPJB8DkBMedISJtItIqIrckejIRWR46rq2jIz+321MqHY3YEj1nuy/AQV+3Nn3LM5FpvP6IAXi/MRlL2w2b9hGRnwMfifPQ1yO/McYYEUn0MeJcY8wREakFnheR3caY38YeZIxZB6wDa+Q/7NkrlUqvrI/e/cvfMardv9LRiG2o54y9/9Y5U6N+xucPahM4NUhG0j4xP/Mo8Iwx5vGhjtO0j8qoV9bD1r+y9vpd+ox134abrX1/G7+TN9s/2rlkO20A4dLP+xZcrBeAHFEIaZ+ngaWh20uBp+KcyEQRqQjd9gJXAXvG+LpKpdasW8KbvK+pt77sTeBnJcxU5pzG2VOctM+8B3Yw74EdTrpI9/vNHXYar7qqnH5jqPW6OWdSJf3GUDOuIiNpu7EG/38ErheR/cB1oe8RkToReTh0zEVAm4j8N/BL4B+NMRr8VW7x1Fgj/iovdPusryqvdZ+dBsoD2VgpqoYXuwrbHzwNwIQqF3fPv4C1iy+jrMQKx3dcNT0jn9LGVOppjOkEro1zfxtwZ+j2S8DssbyOUkrlq3irsx9rOwxAuy/AIzsP8sjOg06VVlPdtIxcqLW3j1JgTe5uuDk84rc/AWy42XosT+h+v7knXiruoK+bGd4qJlaVZ+0TmgZ/pcCq6rFz/F9qtb7sOYA9T2b77JKm+/3mnkSpuB8urqNEJGvnpY3dlIJwNU9kqefSZ0ZV6plNiUpC733qjagJXy3/zK4BY1ix6TXnQgA4n9AyNfrXkb9StsuXRU/uemryKvDbljRMjwoeW3e/z7O7P+DG7+1g39GTTmronqfeZOXm17J4poUpdnJ339GT3Pi9HVGpuK7uPtp9AWq97qx9QtPgr1SBa5w9hRqPiw5/Lzd+byfXr36R/cf8uEpLeHb3B9r7P4Xsyd2F61pZ88IB9h09afXq8fdSXVXO56+e4QR6gNvqzsbrqXBSQ5lci6EtnZUqAtboc6ezmEgAA9oDKEU2thxy0mqR+ynYi7aq3S48FWW8c7yb+xZcTOPsKWlLu2VqkZdSKg9McrsYXxme4jPAxKpyDfwpYI/2m9a2APDQojmUhOZx+41BsHL87xzvdhbbeT0VWZ9v0eCvVIGzc/xd3X1E1pZ82HOa44HejJ9PoW07aS/YavcFuH71i9y+toWBiISKAbq6+3JusZ0Gf6UKnF3+6SotcUb8djpi0cMvZ7T+PzIn7vMHoyaf8/UC0FQ3jRneKsAK8id6+gArtRZ5sR3IsRS7Bn+lCtyShuncNPsj9PYPcP5ZHn521yd47itXUzOugo6TwYzW/xdi7yGvp4IfLq4jtmJ/6pmVGKy8P1gXhlxabKcTvkoVCXtS0k47ZKvW3+cPMu8Bq/QRoNrtYvuquTmTDhkpnz9I09oW2n2BQY/N8Fbxw8V1/Pytozyx60hGuqsmO+Gri7yUKhKRAce+EGjP/7FrbjvsBP6JVeWANcoHuGn2R5k5eRwzJ4+jqW5aTv2ONe2jVJHJZt69EHsPvXnkBADlJcJPv9DAT7/QQHmo3Oegz+8clwsVPpF05K9UkWmcPYVNLe84eXcI7/ub7rx7ZO+h2M1mcmlUHE+itNnFUyfQevA4nX7rIjZgDH0DVl/+by64JLsnPQQN/koVGXs1aWzePRNliOnY4jIT4rVljlzMNcNbxcCAcX6fYPXlz+V5DA3+SqmMig3yuZYOiSfy09L1q1+kRITOQC/nTqrCH+zjoK87qtqnVITrLpqctfNNhub8lSoyw+XdC20RVip4PRXcOmcqYE3mdgZ6rfUSJUJnoM9plwFWbX+/MazcvCun5zE0+CtVZIbq+X/vU28U3CKssYi8EF530WSnbQPA70OdOas9LiIL5s+sKqfW6875PRQ07aNUkRkq7944ewr7jvqzMhmcayLz/A8tmsOKTa8NattQIuB2ldKJVeZpp4POrCrn7vkX5HQ6S0f+ShWh2J7/dt49mQ3giyUtFLka+ZYH/59Tyx/dsgHePd7jrJy2P0Ud9HXjqcjtsbUGf6VU0jK5RiCTF5l4r7V19/tOnr+7bwAgKrdfWW6Fz5pxFTy0aE7W+vKPlgZ/pZRjuMngTPXmyfRFJtFrBYKnnd48EA78tV43T335jzj/LA8dJ4O0tnc6x+RD9RJobx+lVAQ7EMZbhGWPZjPRm8cOwPuP+aP2uE3H5jM+f5Abv7+TjpPBqNeq9rgYf0aZU8ZpR8oSgW1fmcvMyeNysi2G9vZRSo1YrizCytRCtI0th/AHT9NxMkhpaLIWrDr9Tn8vnf5eZnir+H13n9OvZ8DAz986yszJ4/JmlB+Ppn2UUlESTQZDdFrI7SplYlW5kxbad/RkXk382p9yHm97j3MmVTpbXIJVpz/DW8XKT34MQZzNWOzGbU/sOpLTNfzJ0OCvlEqavUagZlwFgd7+qJr2RetTl5NPRwO42End+tpqasZV0O4L8F5Xz6DjBcFdUUa7L+CsifjZXZ9w5jxyuYY/GZr2UUolzf4EUF9bzcrNu9h/zM/EqnLcrlI6/KlbD5DqBnCxvXma2w7zWNthOk4GcZWW0Ns/EHV8qQjtvgCeijJnw/V86kWUDJ3wVUqNSqKJX3ux2Fg3jYnsormx5RD1tdW0tneypGE6a144AMCXrjkvqdeInEB2u0oJ9PYDcM6kSrq6+zh56rRz7MSqcrq6+6jxuHjuK/m3yYxO+CqlMq657TD3b9vLo786RPOKBoCoXa7s4JzMrmKRx0ZWIK154QD3b9vrHNdUNy2qw2a8C0C8CWSA97p6olbtAkyscnFmVTkHfd0FMcJPRIO/UmrEYnPyYJVHPrKjHYB2X4DrV78IhHe18get0fXKza/x7O4PnBTM8UAvix5+mY6T4Xx85IWhvraaGo/LWVsQuRH6+h3tPLLz4KhbUNiBv7xU+Nc7r+BvnniD/cf83D3/AjwVZQUb+EGDv1JqFIbKyduTs3bQB6vffVPdNDa2HOLZ3R9QXipOe+QTPX0MGGulrD94mvu37Y3qm79y8y46/FZ1kT1qj90ucbgy0MiLFViVLpFZ/r5+w2vvdBVMPj8ZGvyVUiM21HqA+tpq5j2wI6rF8Q8X1+H1VDij/75+g0DUBaKp7mya6qY5G53H9s3vOHkq4flEfhqIl0KyL1a1XjenBwZ49/jg6p7H296jqW5aUQR+GGOpp4jcLiJvisiAiCScYBCR+SKyV0QOiMjXxvKaSqncEG89QOPsKazY9FpUi2MDrNj0Gj5/kKa6adR63c79thKBBX84NW7f/AmV5fzu9z109w04awu6Qouu7P46Xd193L72JfYdPRm3DcSShunct+Bimlc08NnLz4l6HxND5artvkDel2+OxFjr/N8A/gTYkegAESkFHgJuBGYBC0Vk1hhfVymVg5rbDjuTuxOryp30TLsvQHPbYbyeCv7h07MH/dyAwdn85LqLJkf10znR00ffgMFVWsITK6/isnMnOo/9+R/N4NxJVQAc9HVz60O/SthryL5YNdVNw+0qde4vEWHt4svyohlbKo0p+Btj3jLG7B3msMuBA8aYdmNML/ATYMFYXlcplZvsNsYzvFX87K5P8LO7PsEMb5XzmM8f5M9/9MqgnysV2H/MT3PbYVZu3kW/MVGtk0Vg052X09reyc/fOka128XKT36MO66aQUlEFAv09g+Z/7dz//Zx9vzEys27im6/gkzk/KcChyO+fw+4IgOvq5TKsHhzAY+tuNLJwd+54RUCofbI484oo0SEEz199Bu4YLIHT0WZk5vv6g5PGhsDr73TRVPdNGcv3Z+8cpifvHKYzkAvpSJR7RkSSfXisXw27MhfRH4uIm/E+Ur56F1ElotIm4i0dXR0pPrplVIZMFRvoDnnTnLuLyuRqG0RF1w6lSUN07l7/gUYTNx+OsCgzWbswJ9MGwg7929/Msin/vupNuzI3xhz3Rhf4wgwLeL7s0P3xXutdcA6sFb4jvF1lVI5xl6Ru35HuzOqn1hVzrK5tc5jngqrjXKi0XlseqbfGGq9bmdR2XAj+dj78rkz51hkIu3zKnC+iMzACvqfBf40A6+rlMpBTXXTeGTnQef7EhGa6sLjw+H2GI63uMyEaofskXyxpXBGY6ylnreKyHtAA/CsiGwP3f9REdkKYIw5DXwZ2A68BTQbY94c22krpfJRst06E6WOInP221fNjdoz1y7TLNaR/EiNaeRvjHkCeCLO/b8DGiO+3wpsHctrKaXy31gnXHNls5lCoF09lVIZlUxTNzV62tVTKZWTdMI1N+hOXkopVYQ0+CulVBHS4K+UUkVIg79SShUhDf5KKVWEcrbUU0Q6gHeSPNwL+NJ4OrlO37++f33/xSneez/XGFMz3A/mbPAfCRFpS6autVDp+9f3r++/ON//WN67pn2UUqoIafBXSqkiVCjBf122TyDL9P0XN33/xWvU770gcv5KKaVGplBG/koppUYgr4K/iMwXkb0ickBEvhbn8QoR+Wno8ZdFZHrmzzJ9knj/d4nIHhH5jYj8QkTOzcZ5pstw7z/iuE+LiBGRgqkASea9i0hT6N//TRH5t0yfYzol8bd/joj8UkReD/39N8Z7nnwkIj8SkWMi8kaCx0VE/jn0u/mNiMxJ6omNMXnxBZQCvwVqARfw38CsmGO+BKwN3f4s8NNsn3eG3/8ngarQ7S8W2/sPHTcO2AG0AnXZPu8M/tufD7wOTAx9f1a2zzvD738d8MXQ7VnAoWyfdwrf/1xgDvBGgscbgecAAeqBl5N53nwa+V8OHDDGtBtjeoGfALGbyC8ANoRuPw5cKyJCYRj2/RtjfmmM6Q5924q1X3KhSObfH+DvgG8DpzJ5cmmWzHtfBjxkjOkCMMYcy/A5plMy798A40O3JwC/y+D5pZUxZgdwfIhDFgAbjaUVOFNEpgxxPJBfaZ+pwOGI798L3Rf3GGNtH3kCqM7I2aVfMu8/0uexRgOFYtj3H/q4O80Y82wmTywDkvm3nwnMFJFfiUiriMzP2NmlXzLv/2+BPwttK7sV+IvMnFpOGGlsAHQzl4IkIn8G1AGfyPa5ZIqIlACrgc9l+VSypQwr9XMN1ie+HSIy2xjz+6yeVeYsBB41xnxXRBqATSJyiTFmINsnlqvyaeR/BJgW8f3ZofviHiMiZVgf/zozcnbpl8z7R0SuA74OfMoYE4x9PI8N9/7HAZcAL4jIIazc59MFMumbzL/9e8DTxpg+Y8xBYB/WxaAQJPP+Pw80AxhjWoAzsPreFIOkYkOsfAr+rwLni8gMEXFhTeg+HXPM08DS0O3bgOdNaEakAAz7/kXkUuCHWIG/kHK+MMz7N8acMMZ4jTHTjTHTseY8PmWMKYSNoJP5238Sa9SPiHix0kDtmTzJNErm/b8LXAsgIhdhBf+OjJ5l9jwNLAlV/dQDJ4wx7w/3Q3mT9jHGnBaRLwPbsWb/f2SMeVNE7gPajDFPA49gfdw7gDVB8tnsnXFqJfn+/wnwAI+F5rnfNcZ8KmsnnUJJvv+ClOR73w7cICJ7gH7gq8aYgvjUm+T7/0tgvYiswpr8/VyhDPxEZAvWhd0bmtO4FygHMMasxZrjaAQOAN3AHUk9b4H8fpRSSo1APqV9lFJKpYgGf6WUKkIa/JVSqghp8FdKqSKkwV8ppYqQBn+llCpCGvyVUqoIafBXSqki9P8BFBOT8bvTsuMAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# make a dataset with two outputs, correlated, heavy-tail noise. One has more noise than the other.\n",
"X1 = np.random.rand(100, 1)\n",
"X2 = np.random.rand(50, 1) * 0.5\n",
"Y1 = np.sin(6*X1) + np.random.standard_t(3, X1.shape)*0.03\n",
"Y2 = np.sin(6*X2+ 0.7) + np.random.standard_t(3, X2.shape)*0.1\n",
"\n",
"plt.plot(X1, Y1, 'x', mew=2)\n",
"plt.plot(X2, Y2, 'x', mew=2)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To build a GP with correlated multiple outputs, we'll use a kernel of the form\n",
"\n",
"$$\\textrm{cov}(f_i(x), f_j(y)) = k_1(x, y) \\times B[i, j]$$\n",
"\n",
"The covariance of the i'th function at x and the j'th function at y is a kernel applied at x and y, times the i, j'th entry of a positive definite matrix B. This is known as the _intrinsic model of coregionalization_.\n",
"\n",
"To make sure that B is positive-definite, we parameterize is as\n",
"\n",
"$$B = W W^\\top + \\textrm{diag}(\\kappa)$$. These parameters will be formed by the Coregion kernel below."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# a Coregionalization kernel. The base kernel is Matern, and acts on the first ([0]) data dimension.\n",
"# the 'Coregion' kernel indexes the outputs, and actos on the second ([1]) data dimension\n",
"k1 = gpflow.kernels.Matern32(1, active_dims=[0])\n",
"coreg = gpflow.kernels.Coregion(1, output_dim=2, rank=1, active_dims=[1])\n",
"kern = k1 * coreg\n",
"\n",
"# build a variational model. This likelihood switches between Student-T noise with different variances:\n",
"lik = gpflow.likelihoods.SwitchedLikelihood([gpflow.likelihoods.StudentT(), gpflow.likelihoods.StudentT()])\n",
"\n",
"# Augment the time data with ones or zeros to indicate the required output dimension\n",
"X_augmented = np.vstack((np.hstack((X1, np.zeros_like(X1))), np.hstack((X2, np.ones_like(X2)))))\n",
"\n",
"# Augment the Y data to indicate which likeloihood we should use\n",
"Y_augmented = np.vstack((np.hstack((Y1, np.zeros_like(X1))), np.hstack((Y2, np.ones_like(X2)))))\n",
"\n",
"# now buld the GP model as normal\n",
"m = gpflow.models.VGP(X_augmented, Y_augmented, kern=kern, likelihood=lik, num_latent=1)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"lib/python3.5/site-packages/tensorflow/python/ops/gradients_impl.py:100: UserWarning: Converting sparse IndexedSlices to a dense Tensor of unknown shape. This may consume a large amount of memory.\n",
" \"Converting sparse IndexedSlices to a dense Tensor of unknown shape. \"\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"INFO:tensorflow:Optimization terminated with:\n",
" Message: b'STOP: TOTAL NO. of ITERATIONS EXCEEDS LIMIT'\n",
" Objective function value: -146.909142\n",
" Number of iterations: 1001\n",
" Number of functions evaluations: 1082\n"
]
}
],
"source": [
"# fit the covariance function parameters\n",
"gpflow.train.ScipyOptimizer().minimize(m, maxiter=notebook_niter(1000))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"That's it: the model has trained. Let's plot the model fit to see what's happened."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXd4VOX2tu89NWUSUum9S1HpRYoKigIKNhRRxAL27tHf+WxHPcXjsSuI2EFEsYII2CnSqzSpobf0MilT9/fHykwKSUifZPLe1zVXZvbs2ftNe/ba613vszRd11EoFApF8GAI9AAUCoVCUb0oYVcoFIogQwm7QqFQBBlK2BUKhSLIUMKuUCgUQYYSdoVCoQgylLArFApFkKGEXaFQKIIMJewKhUIRZJgCcdK4uDi9bdu2gTi1QqFQ1Fs2bdqUrOt6/Nn2C4iwt23blo0bNwbi1AqFQlFv0TTtcHn2U6kYhUKhCDKUsCsUCkWQoYRdoVAoggwl7AqFQhFkNFxhX/8e2JMKXtuTZJtCoVDUcwJSFRNw1r8Hix+DDe/DLYtk2ydjIWm3PO8/NXBjUygUiirSMCP2buMhvqsI+YyB8kjaLdu6jQ/06AKDuoNRKIKGhinstniJ1MPiICdZHmFxss121tr/4MN3B/PJWBF0e5I8X/yYEneFoh7SMFMxiqJ0Gy9pKd8dDMjFriHfwSgU9ZiGGbH7IlJfpO6L3H0Ra0ND3cEoFEFFwxT2Xd8V5NTvWSsPX85913dn/7zKRysUijpMw0zF+Kpeuo0viEhvWSSifraKmGCsqCl+BwMFdzAqalco6h0NM2IHEeDCgmWLL58oV6Sipr5E9lW9g1EoFHWKhhmxVwVfPnrGQIlqoeR8dH2K7KtyB6NQKOocSthrivpWaVJcwMt7B6NQKOocDTcVU1nKW1GjKk0UCkWAUMJeUVQ+WqFQ1HGUsFeU/lNh9MsFkbcvMh/9ctHURWVq5Wt4snX2mkMk2x3+18l2B7PXHKq24ysUirqBEvbKUJ6KmopG9iUt6595QdFl/flCXxmBnr3mEM8s2MnEWWtJtjtItjuYOGstzyzYqcRdoQgy1ORpTVHRShNHFsR2LJhsdeXIw2iBNkML7gCSdrPPNYWJa67mqt4tGHlOE+6du5l9iXbsDjc2q4nJg9qecfjRPZsxZ81h9iXaGfXaCgBSsp10amxjdM9mNfRDUCgUgUDTdb3qB9G0D4GxQKKu6z3Otn/fvn111cy6EL5oPbYj5KRCbmqRt3Mt0ViNRgy5ybhju3Bd7t/ZkmoBwKhpeHSd9nHh6OgcTM7h+XHdRdzXv1fkwpJy+hjvzXyVmbkjAIgNt3D70HZM6NuKxdtP+gXe93zx9pOAXBTibFZA7g4Wbz9Z4sVDoVDULJqmbdJ1ve/Z9quuiP1j4G1gdjUdr3RO74SY9mAOrf5jFxNC7EnVW8udcgDSDol4O+3gzAajGbpfXVAaWQQDuaZwQp1pAHhD48iY8A0Zc/YD2QB4dB0NSMtxkpbjKojA8y8W6ctn4Jn8PbHhVhrNv5r/0/eQbfQwx3MpeS4PLy3dw0erDpGU5eCjVQfR0EhIzuatX/eRZHcCMGfNYeZNG0hqtpNJ7631b1firlDUTapF2HVdX6FpWtvqOFaZeL0w9zrIy4RuV0KPa6DtEDBZq37syiwocuZAdiJEt5XXB1fC0XVgPw1H1oDBBF4P3LVSLhJf3gKnthc9RmQL0Axw7SfwzoDi3zAh7iz/q6zcXO6c+SMJuU0xaODV4WbjTyz2DCAlpxFhFiPPjWzMgR9eY7Z1CFdprWiXnUDKjEGkahBDJnu9LVjsGYDJoJHt9ACQlCX5+oPJOYDcBfjE22SAfYl2Lnl1OZm5bjy6TrzNUmL6ZvaaQyq6VyjqAPUvxz7+Hdg2H3YtgK1zwRwOlzwnwuv1gNddOaEvvKBo+gBAl5RIZEvoMlr22fYlbP4EMk+IeDvtsv3J02AOgd0/wLp3wBQC7jwwh0HnyyHrNMy+suAiEd0OJswBSzjMu0EuKJ1HFxuQBmhoePGiYUCnEdl8pT/CHktLlnj7Y9Zd3Gv+npuNPzPR+RQ4ofHXj9JRO8Yi1xSu1/8fS8xPEKtlApCsRzLR+RTphig83tJTcB5dx2zQcHl13F7ZlpbjAmS2fe7UgX7x9uGbnPVF9wATZ61lX6L8jJS4KxS1R60Ju6Zp04BpAK1bt67cQQwGaD9cHmNehoRlsO9niO0g75/8Ez64FMLjoXFXiaRtTaHjSDixGTqOgP2/gscFHgc47DJpOfh+iGoFA+6GHx6G3JSCc2YeE5EG+epxQrNz5bi2xvIAifj7T4ORz8pxP7oMUvbDweUwY5AcM6q1XHzSDsJHlwO6XBzCG5NzaB1hgNMQhtmbh4YX0NHzRf2AtylWnLQ0pGLFyQPGbzFoOqm6jc6G4/xofQKAOK0gKkfOcAZliboPVyn7GA2Q7XAz/ff9HEnJYXDHWFKyHVx5Xgs1OatQlMbJP2H9LLjgYYjrWOOnq5bJU4D8VMyiWps8LSkfvuF9SYUk/A5Gq+Th89KhUSvIOAq9JsOWQtMAmkH2u/5T6DQSFj4AWz4FXVIUmMMl5dP0XBh0T8nj8Hpg1evw6/MS3Tc7Ty4sO74GZ1bJn6kgObqFp1y3sdnbkaHGHczxXEocGVxi3MgRb2PesEwnLj8qT9PDuc74JgAz3M/Q2XCcZD0SENHfT0vuNz/PX1khZz2vRskXBoAmEVbSclw4PRLSD2gXQ48Wjfhm8zF/dB8dZmbqsPbcc2HN/yErFHUOjwv+WgjrZsHRtXIHP2469Li60ocs7+Rp/RR2Xz48vuuZ+fARz8K2L+R5WBzoukTL8V3hxq8kZWI0Sxnhlk9hyd/kva5jYeXLBecwhIA3r+D1OVfKpG12kqRhEnfJLy4nBXRv+cbdqBUpXhsx2fvRvCJ+utHKvtiL8Cbv55grkkgth6ak0EpLQtOKftyra+zU2/C793x+9PRnp96GWDL50fqEX9gB9tj6s4/WjLV/xV5vC0nTAPND/kUHjvGyaSqxF93Li0t243AXjL1NTBhGg0yehpgM5BV6z5eaKQmjpqGj49ULqnR82z26XlClo1A0FDxueLuv3J1Ht5O7+fNvhNCoKh22VoVd07R5wIVAHHAaeFbX9Q9K27/Kwl6opruIf3hhoS/uvnj3GjBbYfXb0Ox8ScWkJsCatyE3TaL38gp0cUKjpcLF4yxlBwPgJSOsDcl2Bx0Mp/CGyrgNucns9bZgmuEfHMoLJ5YM5ln+SWfDcVJ1GyY8RGq55OgWzHgwax7/UY97YwnVHMRodlL0SAyal2jseDBgxMtObxsmO58gXYvGo+v0jnEx3L2K1zIvJN5mIcnuJNxixGIykJbjol1cGO/e3Je1CSnYHW5eWrqH+AgrMyb1Jtxi5MHPt/pz5uXFYjSw6IEhdG4SUbmfrUJRXzixFfb+CBdKWpS170gw2PESSSNXA7UesVeEaknF2JNkktOXD7c2gkv/CQd+k6vi1s9EvAHQJEovVXjLQDNAp1GSgw9vDBFNwBgCv/8T0g9DaCxomlxENGNBGiefXFMkp6/+jrgFN2JznAIktXJIa4lJd6HpHjoZTrDIM5BIsonQcuhlOMBpbxSPu6ayS2/LXMu/6Ww4ztvmW2lzTm8yNn/HKONG4rUM/3n+9LYns9NVtE74gjb6MZZ5zuVwWA+yBz3KqDZG/pj3Is9mXsnjl3Vh5/EMfth+ik6NbWdMdBaOrkuqcvliwxFmLU8gI8/tP7fFaCAi1ESKveSf70Vd4vjfdeerahlF8OFxw54fRMSPrJH07X0boFGLGjldcAr7L/+QCcms05B1UvLmFcEcBl6XpFCMFhFtdx5lZpNDY+De9QW5fK8H8jIgcTfMux4cmUV292JA073+NEqaHsbl7te4np+4J2QpGR4Lx7yx2PVQnJh41nULFxu3oqFztfEP4rV0mpKKUdPJ0830cHxAjJbDfcav2GEbwuN33sri3ekMbBtFwuZf6b33deKy92NwS6mibgpDc+eQEdoa95QlxIZb4b2LIOMoyZHdiZvyKcS0r1Rpos+GYF+indhwC7quk5rjolGoiSUPDuWdZQf4bN0RPIV+lL6fbJjFyHf3XkBMuMV/jDE9mzJ9Up+K/Q4VirrCqe0wb6LoUFQbGHAn9LoJQhrV2CmDU9hnDJLcdnF8aRRzuETMPrG2hEvVicEsgj76ZZlwLZym8UXZIVEi2CUJfGxHuPDvsPhvkrYpdUoRvKZQDO5cdB2/uO/1tuBG51MkU/YvPNxiJNvpIT7ciCX7FLeZFjPPOoH92aF8H/IsPdmHVzNh6DwKuo2DLpfJH5EzG3Z+K7PuJ/8sOGDhCeTQGLEo8Lrljy+yJfSZUqHFWL6SxpIi/ccv68K3m4+XmaoxAI3CzKTluPw5++fHdfevclURvKLOk3JA5tjaDJb/u69ug143Q5fLwWCs8dMHp7D/9b1E28c2wNoZENsJbl0i7xWePF3xMrhkZaZ4oP8glTItestnf3uhoIQRRLjbDoNNHxZsM4Xk18XLJCdDHxXhD42ROvkNH0DWCVL0CIwaRFFQAePWDZg0L0neCDKx0cFwkqddU5jjubTEbyvOZmFi/9bcMritfzn/rgWvMGzff3HHduGbnu8womMk1jljJZ0T0kjG0nUs3DBXDuJxyYKoo+vhj9dg75KCExjMcP1caH4erHxFxq57JLU0+XsIjy34+XW7CiZ8XOqvoLRIH+CZBTvpGB/OY6O68Nn6I6zYm1zmr9OX0/d53ahJVkWdRNfh0EpYMwP2LoUm3eHuVQEZSnAKe2FKW/5fPCK3RsL9myVafePcgs9rRrnCepwi7JO+gU0fieCHx8Gwx+SYW+aANaJoJJtfleOI7sxNrqdISM7m65DnactJ/y7JeiSjHP8FYLRxXRFRNxng6bHdaRkdyn8W/8X+pOwzRa2sCeLJCyH9iIy/RW9IPwozh0gZVe/JENECpveHvLSiP7M2Q2DYo+AF5t8IrlzZbgwBT56kpzzOs4p7aT/72Z5Lioj+xsOpPDhvC8fT80o8jAZEhpjIyHPTqbGN6ZN6szYhRYm7ou6w/xf45Tk4tQ3CYqHv7dDvDplrCwDBL+wlkZIgC3/sp4puj+8qkemy/0hUHtcFpvwg7/nEs7ifehnMXnOIfYteY3f0xfzr5ou5/t01xOYk8J31WWyaiJhP2FNKSL+0iQnj63sGE2ezlp3btiedWd1zz9ozOzClJsCyF2U1rjtP7ijcDpnYdeVCfv7dj29lbEkUTluV9PMoq9S0hM/sPZ3F2Df/8Ne7l0So2cCC+4aoyF1RN8hNk8AvJFLWoyz7r6xjOff6mvGoqgDlFfbg8mP/ZIyIutEC/e+EG+aJiCfthr8WwBWvifhM+aHsJhlnYXTPZiwNHcuGZBPXv7sGQ04yMyxvYtPySNYjSdYjidMymWf5J7FkFPmsUdM4nJrjT1/E2axVF7GY9nD1LHh0D3S/SkQd4JbvYervoOUvMDaaAa1A1LUSfv1el9yhdBhR8rm6jS/wlZ8xUB4+3/kSermuTUgpU9QBcl1eJsxczb5EOx3iw9VKVUVgSD0Iix+HV7vD+ndlW7erJJjqMyXgol4R6q+w67pYCsy/pSClcOVbMPgBeGgnjH4Juo4WES8s3OVpknEWFm8/SZLdiUETD5XRxnV0Nhxnr7cFoxz/5TLnfzlASzobjjPauM7/Od+CndJMtIpQmQ5MoVHQ5gJ5HtVa7A7CYyXnAZKHR5eIvbS6faMVYjoUGJtlHC/6fgV7uU4e1JaR5xTdrp2xF6TnurEYNbIdHlbtTyYQd5KKBsqxjTB/MrzVGzZ+mF+YMEbeMxiqrQa9Nql/JmAgv4ifnoYjq8WzJXmf+Ld0HCmPwlRCuM/GsM7xRIWaSc+ViVVf/nyxZ4CkXnSYkPck14dtYo7nYtrFhfntcH0Lg85aBVK4A1PxlEdZ1SuFzcx8DTu87qL7+CN2TS6QhfE44MoZ8seclwnvDIImPWH449BuGGcshz0LyXYHh1PkwhsdZsbp9pLt9GAyQNdmkRxMspPtlAuM06OTmefiwc+38vHqQzw9thu9W0dX6HwKRbkoXLa2/L9iRXLBg3KnH1n/7xjrV47dmQ0L7pXSvvB4GP6ETBZWh21vOflp5yn+37fbSS5hMY5Bg0ahUs43pmdTnhvXo9QGFuVKv1TWH754br4s/JF7oVp+zQDjZohvvCsHNn4g6wZa9IEmPeSfoLRVv8Wi9pJKJC9/YyVJWQ7/Ra5RiJHMPI+/iNSogdVsJMfp4aGRnXhoZOezfx8KRXlw5cKf82DtTJg0X+5M04/K6nGrLdCjOyvBOXmq6/Dp1dCynzgyWmtvmXqy3cE/Fu5k0baTmI0aLk/pPzefiBW3tq01ShR2sTUo8jyiBWQdFyG/6TtZDr1uetFjNWoFty6F7V/Kaltf9F/OyVMouUTy2QU7/KtfP71jAB/+cZB3VyQU+VyE1cTUYe14YERnvF4dg6FidwsKhZ+cVLmTXfeu/F807wVjX5Ov9YjgFHYoegtVjZxtJeZ9n21myY5TaIDbq/sbXRQ3y4oOM/PzI8MDK+ollUkWtjsIjZELY6+bpQTUlV81YwoFdy5YIsBoyl+MhVTKWMJloVNkc+g3VT57aIWUUGqGSnWaKv4zX30gmUfm/8mpjDwiQkxk5dsWXNWrBfE2C6cyHTw/rjtRYZaq/pQUDQlnNrx6jqz96HSpzMO1HVIjOlLT1HZrvNqjhkS9tCYRuU4P43u14FRGnt/H3GTQcHt1f+319e+u8VvVOt2VNBKrLorn5rNTYNZQqU83h8vMfk6yOGD2uhke3Aav9xRBd+cCmizGMpql0cjKV6VSJi9dJlZv+FwWOuWkwoL75IIx/HEYeHeFh1o8HTW4QxwrH7+ImcsO8MXGo9x2QTtmLj/At1uOE2I24HR7WXcwhVHdm/LAiE7Ke0ZROie2SK+G4Y9LUHLpv6BlX2h8TqBHVivUv4i9BijsgRIdZsagaaRkOwkxG4iwmnB7ddJyXP7ovFuzSEac05grzmvur70uPDEY8FRM4dy8r+7c1gTuyl8tVzh14siCX58r+KzBJAsx7Kfh4qfFDbPwQidTKAx9RCL+la9Jvv3gMmjcHS55XmxKq2GyOs/lIcRs5EBiFjd9sJ6TGTLhazFqOD06UaFmljw0FLPRUKKBmaIBoutiArjqDWlwY40UQ66IpoEeWbURvBF7DRBns3JV7xa8tHRPQQs4DfJcXvJcMkkaajaS5/Jw1/AOPHJJZywmgyxUSrSX6J0S0AiysLD6nheehL1lUcEq3Y8uK/RBTXLo9tPiIb12uoh6SJS8l5cmUf3v/xJLh9w0uTMY/44YtM29Fv8EbP+pVWoOHmIW340cp5esPDdmo4ZB0/z+8em5Li55dTlWk5GUbCft4sKwO9xlHVIRzCTtga9uh9PbIaKZBBl9ptSoIVddRkXs+ew9ncXlr6/0N4koTJ820RxLy+H163sxqENskffqdQNnXzQf21FSK7mpBe91Hg17F0v+XDNC31tl3UDyXvnHySqwTyAkWlJkuakFKaBNH8kFoAKTrKVxMiOX5xbuYunOU36jtMJEhJiIt1lJSC7BmkERvDizIeMYxHeR0txPr5EquXMn1GqlXG0SvJOnNUDhVExxXrr2XK7u1YLMPDcx4UE4abf+PWgzVFbt+qpoLOHwwJ8SXbcZCmvegq2fyi3tOeNg1L9hw3vw89P5i57yCYmGCx6Axt3gi0lSQpadVK6yyPLw++5Enl24k1yXh6Qsxxnvd4gP54s7BwUuBaaoHbJTxMl0/SxJMd6zpl5OhFaGhmkpUEl+2HaCfYl2mjY6UxBOZuRiMhqCU9RBUiVf3VJ0haszW6LrbuOhSVcYPx1u/1ncIM+9XipmDKaiog6Sqvn1OUjeA+dcIaJuMJdrhWp5uKhrYz6bOoBwS8n2qB6vzuw1h0i2F4h+st3B7DWHKnU+RR0j7bBYZ7/WHZa/CK0HwRVvNBhRrwgNXti9Xp2DyTm0iQ3jVEaBIPj+VBb9ebKIUAQdhato7lkrD58XzK7vCvZr1R/uXAEt8xtjnNgqKZoSMcK1H0lk73WVsk/l+G13IodScmgXF0aIqeif76GUHN78dT9j3/yDZLvDfyf2zIKdStzrM76swpE1sPEj6HGNNL+Z+Bm0HhDYsdVRGvTkqcvj5dH5f7LwzxNFtneID+eV687jb19tC/xEaE1T1uRq8Ty4LzLyuKTrerE2gH5+fhJO/Qknt+Z/ziDdq3xeN1WI2icPaovd4earjcfIc3tLzLmfyszj0leXo+VXN3VqbFPGYvUNXYeDK2DV69D+Qlnu3+MasbWIbB7o0dV5GpSwz15zCLvDzYS+rQi3mLjjkw2sOpDifz/OZqFjYxvv39IPm9XEvGkDg1vUfRQX8LP56xjNcO8G+PzGog09QqKhzSDY/ytsny/botvBzd+CxVYwebrzWxgwrdLDtVlNJCRn+6uRTqbnMuHdNeS6CtYQpOZXNzUKMfkrlmYs24/Nagr+32d9xuuBvxbCH69LYBDeWNJ6IH93StTLRYMRdt8iJIAvNxwlxmZl0+GC+uyrerXg31f1xGTUMBvlFr9aLHWDlZwUqVkHybd73ZJjT02Am76GHx6RCpqsU3Bis0RbY9+QSdqUffIPXN5WYsXKJif3DKfLke10GPMwcTYrcTYrKx6/iBcW7aJv2xieXbDT7zuTkedm0bYTfLL6EAeTC3zp1e+1jrLgPvjzM3EYveINOPcGMIcEelT1jgaTYx/dsxnt4sIAOJiSU0TUTQaNv43qQqjF6Bd1xVkonJt/aAfc/F1Bbv7wKpi2TCZa3bnSF/LHJ6Xb08C7pZrhyykFvvFl4SvJ9NkV51smDPjrP8Ttmu3fLT4ihGeu6M6cNYfP6Ej7j4W7/KLePk75vdcpclJh+f+kbBGg720wYbYsLOozRYl6JWkwEXuczcrL153Hze+vI8dVdNn/QyM70Tyq/pjo1wmK5+Yjm0nu/KenxIbA1hjGz4QWfeHHv8Oat6Wr+4RPIKK5bPv8RpgwByxhpZ+nuA0xFJRNFmvssXj7SfYl2ukYH06ftjF8seFokfeNGkyf1Lt+rjkINtKPSA/RzbOlP7EtXoS8Vb9AjywoaDDCvu1YOrd9vPEMUQ+zGLmhf+sAjaqeU1JuftS/RXgXPQyHV8PY16FpT2lkcHA5vD8SbpwvtfLfPyiNDQbfV/o5fI09ircILGEC1ifQvoj8h20nsDsKJlY9Otz+yQYW3jcEoMjaBSXutYSuw3d3w7b5Mhnfc4LYUzTpFuiRBRUNYoHSH/uSmTp7I7mukqs42seFM/8utbCl2vB64Y9X4Pd/Q2wnuGGudG2adwOc3iELl67/VMolW/U/e669vL1f80m2O7hu5uoiOfXCmAxgNRn9vj6qiXYNo+syz9Iiv1R20cNSJTXwbmjUMrBjq2eoBUr5LNp2gikfrS/Sd7NZIyv/vaanP+eekJzt70GqqAYMBhj2N6mGyUmRSdSoVnDbUuh8mXjMzB4veVWDURodfDkFctPPPFYpLQLT37mUlNPH/LsVXoi0ePvJIjn1JQ8OpXlUwUXb7YVspwerSePNib24d+5mVeteE7id8OfnMHMIvHexpOJAfNBH/UuJeg0S9KmYrUfSMWgaTo+XMIuRyYPacMfQ9sTZrIw4pwnzNx5VJXA1RfsL4cGtBQ1REndLpP7T07DuHfjmDsg8Dk26w1+LpOXhzQsgroPsb0+SlYbFWgSmv3MpUdkJvPXRm0y87wWg5LSKr7Q1zmZl4X1D+fCPBGYuS/C3G3G4dca+tRKPV5qjnM7MY8ay/dxzYUdA5eArjcMuXbfWzoSsE/K7Gzcd4lQnrNoiKIVd13VOZeaRmu3k2y3HcHq8fndA3z86yISq759YUUP4RD15nzhJdhgB17wnEfyPT8Ivz0rjjvMnweaPYeZgsRcOaVRQ997tKhj9P3/qxTP5e9766E1eSR/Ox6+tADhjIVJJYvzzrkSKNQHE45XXA9rHMP33A/59J/RtpXLwFcXtBJNFSl+XvyRVUFe+Kb/zetgQuj4TdDl2Xdd5cclu5qw9DECO04PZqGGzmphz+wB6tGiYNp4BR9clilvyhPSZvOEzSNwF39wpDbS7jIZDq8CRIbn3kEbiFunzkS+WT0+2Oxj12gpSssVWOTbcwo8PDyt1nqRw79UXxvfgxvfW4i30px9qNvgXOBX25A+4t3594OgGqXpKPwxTf5dJ0axTQeWDXldokDl2t8fLE19v490VCeS5PH5Rjwm38OVdg5WoBxJNg353wC3fS4uy90ZI046bv5GGCHsWi6ijiVVBbqp8xn66wLPGniR17ZVg8qC2PD+uO9Mn9ebp73bg1WX9go9clxerScNiMpCW4yIl20lsuEWJeml43LKC+P1L4IORcOB3aJvfqQuUqAeYoBH2XKeHO+dsYv5GmVDz6jCqexPaxYXz1V2D6di47ncgbxC0GSyLl2I7wKlt0ntyyg8F1r6FlxfpurhDthlaMIm6+DHsK2cwcdZav/jGhltIyXYycdbaMg3bJg9qy9qEFH9zlEUPDCHCWlCR43DrgW9tWF/Y+a1MeGcnweX/g0d2waUvBK0Pen2jXubYS2pu8diXf7JsT5J/n0cv6cx9F3fEq4NRdbevWzRqCbf9CMZ8K2SDEXpcBxvfk/xsYbwu+PBS8QnJX5i02DOAfYmnKtW5yvfewPax3Dt3M1kOD9FhZjJzXXiKZSV9FwsVtQNJe2H9u9IztN8d0O1KsMyDzqPKbw2hqDXqnbD7cqVv/LqPeVMHEhFi4sZZazmYIuVtvomx81tHoWkaRqXpdY/C3i+OLFm05Cq55hzIT9Egtc/nXs+EoX3Isxa9uFfEsG3yoLZntDV8ZsEOFm8/5d/HajLgcHuD392zLLweaQi9/l3pJWq0ymIikMi86+jAjk9RKvVO2H19LVPsTka9tgJNo8gkmA7cNLA1g9rHlnwARWDxeb9seL+gZV5IZNnC7sOVI408rBFMHlR01WtFDdsKr1KNs1lHW9QaAAAgAElEQVSZMakPD3+xhYVbT2I0SG/VELOBif1aNUxRB/juHtj2ubRCvPgp6HMrhMed/XOKgFMtVTGapl0GvAEYgfd1XX+xrP2rUhWTbHcw5o0VnM5ylvj+o5d25r6LOqKprip1E1+uPGl30ZZ5IdHiDnk2NJNUyTTpWiPDW30gmamfbMTt1XG4vYSajbwy4TyW70ni/hEdaRldhq9NfefEFrngXvj/oFELsYTIOiW2uUZzoEenoBZ7nmqaZgT2ApcAx4ANwERd13eV9pnKCnu2w80j87fy487TJb7/6oTzuLq3Ws1W5ynJIuCWH+C9C8UNEqRixvfcT36i7fKXoP+0GmuJtuN4BqsPJLN4+0m2HpU0kNVkIMxi5OreLXh6bPcaOW9AcGbDjm/Es+fEZkl3XfOBSrPUUWqz3LE/sF/X9QRd153A58C4ajjuGYRZjNjz3ISYzxy2Bqqcsb6ie2H+JBHy0BgIjS1B1JHFTiOeBXsifHmL5OdrgB4tGhFiNrL1aAYRIZKtdLi9ZOS6+OCPQ9w/b3ONnLfWcWbDaz1g4X2S5rrsv/DobiXqQUB1CHsLoLA/6rH8bdVOlsPNqcw88lxnlqTpcNZyN0UdoCTvl9xUSNkPsR2ll+WUxVLmWBxHJmz6RBp7/PW91MIn76uRYY7u2YyO8eFk5RVU6Xh1MGjw/Z8neXHJX/7t9aZhdl6G9Axd8oS8toTDsMfg1qViqjbwLlkYpqj31Fodu6Zp0zRN26hp2sakpKSzf6AYGbkuxr+9igNJ2f5tZqPGp3f0p3WMeKmnZDuVmVddp7Tm2SC2ArZ4OLyyaBPsEc9JRQZA+iGpyLj5O7k4zLpIjKaqeQV1nM3K53cOIjqs6AXGN1H/3Zbj9aNhttcrvUO/uRNe7gKLHoKDK8GVf0c06F5pZ6jmpIKK6qiKOQ60KvS6Zf62Iui6PguYBZJjr+hJGoWaGdghFq+uczIjj1CLkXlTB3JOs0i+uecCZeZVXyhP82zf1+wU2LsYlv0LxrwKvz0vaZjtX0lThmnL4ZupsPB+aD0IottU+3ANxQQvwmrKv3N0MPyl3/B4Ic/trXsNs3VdxHrTR9Km0BoJ50+EXjdB895KyIOc6pg8NSGTpyMQQd8A3Kjr+s7SPlPZydMDSVlc+84aQs1G5t81KLgrFBRCTip8PBbSDkmDjh8eFhvg5r1h8gJJJ5zYAi3z55OObYKWfap8Wl8kvi/RTmy4BV3XSc1xEW4xckHHWH7alVhk/zuGtOPvo88J7GK4zJOw42tpJN5/moh4djIkLIOuY8CsuoTVd2pt8lTXdTdwH/Aj8BcwvyxRrwptYsIZ36sFc6cOVKLeUAiLkebYbQZDdGsR86g2UsExazhkHC8Q9R3fwvsXw+eTCnpoVhJfm71OjW38+PAwfnpkOJ0a28h2ejivVRShhSbwNeD9Pw4y6f21JGbmVem8FUbXZW3AR6Ph1XPgpydlRNZIeT88Dnpeq0S9gRF07o6KIMfrhdQD0rjBkQnmcJlwNVnh4zGQvEcmXo0WuOABySH7rIMrSEnWFfM3HuXbzcf9dr6+lc4aYDJqtI4J46eHh9ds5J5+VDpRdblcXr87DFx50ONq6H41xCvf82ClvBF7vVt5qmjAeD1S5hgWK2mZj8dII+S3eoElAnJTZCJ2/Ez441VY9h/YtRDuXlWpnHLx+Zo4mxWb1cS+RDsd4sPp1DiCpTtPYdSkn6rLo9OjRSPmrD3EmJ7NiAm3YjRoVW/Y4fVIumnPEti7VETdFApPHJRIfPICCIlSeXOFHyXsivqDwSglkX+8Bi37wU3fwJzxYhWbmyL1774m19fPgeObxfZX06QJxMpXJO8c1ers5yqFwlYEseEWZq85zPOLdhFi0nC4vCzYeoIFW0/w5i/76JDv/X7/Z1sq3rAj/aikUcyh8v3+9gJoBpkkvuSFojnz0OhKfz+K4ESlYhT1A59xWFgMzB4HxzZIZczWeQUmYaYQeHA7RDQ+8/OHVkn9PBp0uAh6XCviGBJZ5aFtOpzKm7/u55ZBbXjsq22kZhfYXWiapMHLbNih65CaAEfXwZE1Uo6YdlDuSjqPguT9cHIrdLhYvn9FvSTP5QEgxFx5N8xasxSoDErYFRXCZxzm63uanSQNknX5R8EaKfl2kOj1nvUli3v6EVmgs+MreW60wn0bpEzSkSX5+iq2cDuYbOeGWWs5nVl0odxDIzvx4IhOaCATu6d3SjOK5ufD6V3wziDZMaQRtBkC7YaKR4tq+BwU5LmkX4RBgw+n9Ku0l5XKsSuCh27jxZwqabd4zECBqDdqJe3YDvwO306F3DT47k5J0xT/54lqDSOfhRHPSOPsvUtkG8DSv0vziKbnQrPzoHFXSfu0HVKhof51Mov0zCxaaelo6BzRm6Ch03L5o+TuyyYsbS84860Q+twKzV8nOawd+7o9w6ALR0NcF9UfNMjIc3m469NNLN+bxH+v6VkrBoVK2BV1H1u8ROrFjcP6T4O+t8n7502Q1aoL7xPv8BX/g+GPl3w8TYNW/eTho8vlkrM+sRU2fSxeNXGdJaIHmHud2BdYI6QCx2iRC8Bl/5H3P7seTm3nstx0RofI6uifPX2Y6noUHY1O2jHcxLO7yWjmJIST3agzT10wAewOJr63nn2JXXm+VSiTGytRDzacHi/pOS5evLon1/drXSvnVMKuqL/0vU1yziv+J1F9r0lgDoGv74Df/yWTrUMfLd+xuo6RB0gVSuZx8Vbx0XqgpHycdpmsdTvBXSjdEt0OwmLZnaaxeL8Dd3hTTllaQzIYNY1xzn9iOWbg5oGt+Xz/ITypsHK6mIn5mmbXqZWriiqT5/Kg6xAZYuaruwZhMtbeRVvl2BV1n9I83OO7wrUfwcejZdHS7T+DySLeMd/eBehw8dNidFWL+OrfbVYTLy7ZzeieTXlh0S62H5d5gOgwM2k5BV44seEWfnx4mGq/F0TkuTxMm7MJXdf55Nb+GAwadoebGb/vZ8oFbWkcEVKp49amba9CUbOUZhyWtBsOr4Ir35aqkd//KfufdwOMnwFoUia44n+1OtzJg9oSZ7MSYjbyjyu7079dLN/dcwH92kajAWk5LgpnWQMRXClqjjyXh6mzN7JyXxJXnNscTYMFW48z4pVlzFh2gBeX7K7xMahUjKLuUx7jsD63wqo3oONIaDcMzr9R6r6/vQt++6ekV4Y/EbBFPCcy8jiYnO3vpepbraoDqTku1TS7nlJ8dfLhlGxu+Wg9h1Ny+O8159KrVRQ3zFrLuoOpAJzXKoopg9vW+LhUKkYRHDhz4N2h4HXDfZvAmB+z/PkFfHeXNPMYdB9c+s/AiXt6LtfOXM2J9DwM+b16wy1GGoWaOJHh4Lkru3HL4HYBGZui4sxec4hnFuz0r1EAGP7S72Q7PYw/vzmtYsKYufwALo9OTLiFJy7rwnV9WmGogt2EqmNXNDxObBEBb1HM3XHndzKh6nXJoqYxr8rEajkpyTOmshYBeS4PI19ZxrH0PJo2CuFURh5Gg0b7uDByXV6m39ib81pFVfk8ipqnuAMoyER4TLiZcIuJo2nieT+xf2ueuKwLUWGWKp9T1bErGh7NexU8z06B8Fh53n282Pt+cZOUMuZlwlUzpWzxLPiisjlrDvujMt8/M1TAIiCf+RuPciw9j3ibhX9c0Y1V+1OYs/Yw+xKlRPLqGav55LZ+dG0WWaXzKGqeOJuVedMGMvKV5aTkrza2mgykZrtIzXbRuYmNf1/Vk75ta3+1sIrYFcHHmuniC3P3alnd6ePQH/DZDbJAqO1QuGHuWVvBlRaVlWkRUInjFcfX0KOy51HUDLPXHMLucDOhbyvibFYWbD3Og59vLbKPyaDxwIhO3DW8A26vl+NpuXRqIg6jv+9J5KIuJayKLieqKkbRcOl4iTRqXvhA0ZZ5bYfAbUvA1hQOrRQP88wTZR7KF5XFhltIyXaSku0kNtxSabEt6Xim/Jyr1VTw75jlcBMZYlKiXofw3b29tHQPE2au4fVf9p4h6gBxNgs5Tjc3fbCO8577iXvmFjQ/L/w7rkmUsCuCj/jOMPIfsO9HWDujYLs9CY6shdt/gthOYn/7/kg4uS1QIwUgMsTExP6tcLiLNmn3eAte15uG2UHM6J7NaBcnDX4SkrN5/ZeijdRvH9KORqEmTmU6eHd5AjlON7dd0I6nx3bz7zO4Q1ytjFXl2BXBiZYfs/z0FLQeLGZavkVOo18Wcf/8RnFT/PAyuOa9gpWnhfClTnyROkjqpLLliaUdb11CKvERFpKyCtIy2U4vl7y6nBHnNGbLkXR/I3eVb68dik+aA3SIt5Fid5CZ5ymy79s39mLsuc3p1zaahVtP8Ny4HsRHBO5OS0XsiuCk+9UQ00GqZD65QnxmfIucfPa/kxfAeROlWcfnk8T3vNicU/EWeT8+PIxOjW3sS7SzePvJCg+rtOMlJGeTlOWkTUwY/3d5FxqFSsyVluPiq03HOZCUTbu4MGU7UEv40i4TZ60l2e4g2e7g4peX8ctfiWeIOuD/W7isRzNm3NQnoKIOavJUEczYk2D6AGnCAWJHcM/agkVOIEL+x6vw6/PyuvtVspLVavPvUp3ljmUdD/BvT7Y7uODF34qkZ2JtFt6b3JferVVjjZqm8CR3dJiZbIcbp6d0rWwTE8bX9wyu8fkQNXmqUEDRxUhed8nvD30UbvhM2uvt/BY+uARSDvh38VkE+IizWauUDinteMW3m4otZEnLdnLtO6t5eP5Wku1iQJZsdzBj2X6Vf68GZq85xIxl+0m2O4izWZk+qTdWo0ZajqtUUff9ig6n5lTqDq6mUMKuCE58xmE5yRAaA2iQly59Uu1JZ+7fdQxM/U0mVRN3wawLZWFTAPBFi9lOD5EhJv+1yavL49vNx7nundXsPZ3FdTNX89LSPTyzYKcS9woye80h/wWyeMXL3tNZ3PLBOhyFBF0DrjyvGaHmAtlsFGrm3os68Py47nVq7kMJuyI4KWwcdu96uPIt2Z68R94rifjOIu7nXCEdmb68BRb/rag9by1QOA//22MX8tNDw4qICcDBlBwufW0FB5NzAGgfF15r+ffCggj1s2LHJ+QTZq4h2e5gYPtYzEa5giYkZ3Ppays4WawLlg4s3XGaXJeX2HALseEW0nJc/LTzdJ2b+1BVMYrgpLhxWO+bYd9PsHvRmZYDhQmJhAlzYN27UlGzfpaUSF7zgQh/LVC4YXaczUqczcqyv13EG7/spVfraP72VdHyzBCzgblTB9RKvXt1r8QNFHaHpOUSkrO55NXlALhKSbdYTAYeGtmJD1cmkJztIj7CypIHhwIF33tds35Qk6eKhkNuGswYLBOjd68Go7ns/Y9vgi9vhfTDYAqFS1+AfncEzEQMJDr2GU0VpkN8OE+N6cbRtJwaFZjqXokbKJLtDka9upyUQr74JWEyaHx37wX0aNGIZLuDZxfs4LlxPaptIr2iKBMwhaIkDq6QValdLi/f/nmZsOQJ+PMzed3xErjyTYhsXnNjLIVku4PrZq72p1+KYzJouL16jed7k+0ORr22wm+FUJcbhZRWgbQuIYUftp8q1zHax4Uz/65BdeL7U1UxCkVJtBtWIOquvLPvHxIJV70D130MIVGw/2eYPhC2fHpGzXtNs3j7Sb+ot44OJTKkaCbV7dWJs1m4oGMsn68/gsfbsBt4FK5Fn7FsP3tPZzFx1lqeWbCzTFFvHmnlp4eHFVllWpcqXsqDitgVDZPNs2HFy3DncggtZ1145klY9BDsXSqvO4yAsa9CdNsaG2ZxCptQASWmZZpGhnAqM4+OjW08OfocLuwST0q2s1pSBvUlFTN7zSEGto/l3rmb/fl/TZNrsc9grTTuu6gjj43qQrLdwfyNR7FZTXUmf65sexWKsmjaUxpW//CoTIyWJ28e2Qwmfg7b5sOSx+HArxK9D38cBt9/9px9NVBYYJLtDkLMRr+wR1hNWEwGTmXKnciBJDu3fryB3q2jSMxycCzfH7wqIlW4Yqf45GldmUD0RerxNgtvT+rNtNkbych1+2+wyhJ1gMaRcnGKs1m558KONTzamkGlYhQNk+a94ML/gx1fw/Yvy/85TYPzrof7NkDP68CdC78+BzOHQMKyGhtucYp7zsSGW8hyuMnKc2HMXzXjE7LNR9I5lpZLp8Y2f1leZUsWJw9qy/Pjuvujc59bZUl5/dooiyzpHHaHm3ibhSS7ePpk5J4p5MUv477X8RHWOle6WBmUsCsaLkMegVYDJWpPO1yxz9oawzXvw83fQkx7qZmfPU48Z9IO1chwC1Oa54zTo9MhPhw4U7yGd47nVHoeF728jGcW7OSGd9f4fVB8uefyivvZVuKW5LVSkXOcDd8qUV8t+t7TWcxYtp8JM9fw0tI9jOohPvylTTP4Nhs1jegwMzpgMRpIynLUu3x6Sagcu6Jhk3YI3h0Ol70I50+s3DFcebDmbVj5qhiKGa0wYJpYFZQ3f18JSqv4uHlgG77cdIz/+3rbGcJmNRkwGjRy8tM3YRYjISYDqTku2seFc23fltWSfqjJXLzvohFrs5Bil8ocXw/ZimA0aHw2dQAd4m3+sY7p2ZTpk8pY5xBgVFWMQlEeotvCg39WXtQBzCEw7DG4fyOcez14HLD6LXjjPFj1Brhyq224hSktck7JdvLeigS8OtisRmLCJPdvMWo43F6/qAPkOD2k5riICjWho/PS0j3cO3dTlVMoJTUUiQ4zM2/aQBZvP8ne01n+Y57t+MWX/ndpGkF8hJUUu9N/V1JeUZ81uQ+PX9aF+AgrHq/OnlNZRdJJdVnUK4KaPFUoQqV5NPt/EV+ZFr0rd5zI5nD1LBh4N/z8jNTM//yMtOob+ij0vkUuAjVMSROcN7y7hv35fu5RoWYy8lxFqjXTc92k57oJMRv4Yfspdh7P5Kt7BgMwYeYaEpKLesFXxvEyM9fNR6sOMv33A1iMBpweL3aHm283Hy915Wrhla5tYsP45a9EzAYNl1fHZIBivUmKYNQ0PLpOZIhUweg6PPntDpY8OJQJfVsVGW9Vjd3qGioVo1CA+MG81ReMJrhzBVgjqnY8XYf9v8Jvz8PJP2VbRHO44EHoPRksYVUfcxmUJLwLth4nxGzk/ZUHOZgv1MUxGjR//bvJoGE2auS6RD0fv6wLNquJFXuT+OWvRDo1tjF9Um8WbD3O/I3HSMpyMKZnU54b1wMoqJYJsxhxuLx4dF3cEHXwIsJrMcnxC1+EfIKb43Sz+1QWT3y1zS/85cGgQbu4cA4kZfuPm5rtZNL760jKctQ5w66KUCsrTzVNuw74B3AO0F/X9XKptRJ2RZ3k8GpxfzxvIoyfcfb9y4Ouw57F8Pt/4PR22RYWCwPuhn63S8OPWuaT1Qd5duEuwq1GNB3szjMbRxSndUwol3ZvyvsrDwJgNIDHKxO0PgUJt0jpZWy4hSvObc7Haw75V8MO6xTHH/uSKSnANgBLHx7GuoQU/r14N7kuD6FmI7kuD5oGH97Sl3s/21IkhXQ2fBeh6vTRrwvUlrCfg1x83wUeU8KuqPf89i9Y8ZLUtve8tvqO6/XKwqaVr8Dx/L99U6jk9gfcXWsGYz5mLNvPlxuPcjA5h9hwab7si8zNRq1EQ6xGoaYSSwd9XNunBV9tOn7GdqMGHr3ga0mEWYy4PV6cHh2TQSM63ILT7SEzz33WBb5NI62YjQaOpuUSG24mJdtVr6PysqiVBUq6rv+Vf7KqHEahqDsMf0Lq0Rc9DG0ukEVJ1YHBAF1Hi53BoZUyqbr/F9j4oTzaXwh9b4cuoyUdVMPYrCYOJucUSYGU5UMDlCnqQImiDgViXkYDoiLRuNurk5Qlk6UGTe4Oysqlh1lMzLy5D2sTUhjds1m9j8qrg2rJsWuatgwVsSuChbTDIrp9b6tZJ8fE3bDuHfjzc3Dn+9ZENIPzJ8H5N0Jsh5o7NyXn4f/v62388lciAJEhJowG6SAERdMupWE1GYq084OCUkRD/pL+so6hAaEWIzlOD6N7NqVPm2heWPQXYRYjuq777yp8+CZhgzVCL061pWI0TfsFaFrCW0/qur4gf59lnEXYNU2bBkwDaN26dZ/Dhyu4IEShCAS56QVVMzV2jjTYOk8i95R9BdtbD4LzboBzrqy1XLyvCqVdXBhf3iVVMb5I/pFLOzF//TGOpZ9ZvukT75LSOEZNo3VsKBaTiPO+RLs/veKrXDFqYDDIZyNDTNwyuC1DOsYxoH2s3x/npaV7/DYBmw6n+atp6nrteXVSq7a9KmJXBCXHNsKcq+G6D6HjyJo/n67D4VWwZa50eXLlp0UMZjl/96ug86gav9CU1Wz7mQU7gZKj9+7NI9h5IssfRUeHmcnMdePRdeJtFpY8NAygiOUvFNjiQkElTWkWBcE2GVpRlLArFFXFlQvvXQz2RLjrj+rLt5cHRxbsWgg7vpKcv56fgjCYof1w6dHaaRQ0alF7Y0ImXT9YedDvUePVddJyCroKPbtgBz9sP1W0zPC9tSTZnTx+WRd/lF14Narv7iDOZm2QYl0Raqsq5irgLSAeSAe26ro+6myfU8KuqDck7ZHG1s17weSFtTKxeQb2RGms/ddCiej1QnnmJj2h0yXQ4SJoNQBMNWub60vVlOTu6Iuyzxbxl/VZRdmoDkoKRXXx5+fw7Z0w9DEY8XRgx5KdLHXxe3+EA7+LN40PUyi0HghtL4A2Q2QFbQ0IfVVSIg0qneL1QOYJSDsIqQfzvyZIH94eV1fqkErYFYrqZMF9YDDC2NcD2vO0CG4HHPoDDvwm6ZrTO4q+b7RC8/OhZT9o2VfuOqLa1J3x13d0HXJSpCdu2mFIPyKmcoVfe0voqTrgbrj8xUqdUgm7QlGdeD0i7HWZrNNw+A9ZQXtoFST9deY+IVHQ7Dxo0gOadIcm3SCuM1jCa3+8dR1njkTcmcflkXGs0OOofHWVXvcPgK2JGM1FtxN755h20Oz8Si9IU8KuUNQEp3bI4qJx08FkCfRoyiY3DY5tgmPr4fhmOLEFcpJL3rdRK4jrlC8+7UWIolrJ9pou96xNdB0cmWBPguxEmb+wn4asU/Kwn5IWiFknIS/97MezNoLo1nInFJX/NbpNwddqvmCq1ngKRU2Qsh+2zxef9dEvBXo0ZRMaDZ1GygNE1DKPw8ltkLgTTu+ExL8g5UB+BHpU0jrFsUbKwqnIZvI1PL7gERYj5wmNlv2sEWAOrdl0j65LxZLTLtVDTjvkZUJehoh2brpc1HLTIDdV0iU5+V+zk8DjPPs5AIwW+X4btRTnzkYtIbKFXOx8F72QyJr7PquAEnaFoiJ0Hw9H74G1M6BV/+r1k6lpNE3EqVFLsTfw4XFLXjh5r0zypSbIRF96vtg7MuWRvKec5zGCxSYCbwkDU4iIpMkq5ZoGozw0A6AVdJnWvaB7JO3lcYkAe1yyKtftkDaEzpz89EcVMg3mcLDFQ3hj6YRlawIRTQu+RjQTIQ+NESuIeogSdoWiolzyvKQ2Fj4AjbtJnro+YzSJfUFJFga6LpFvVn56IvOkRL3ZyZLK8EfGaRI9O7JEiB0Z8qixMVvBapMLiDUCQhrJHUNIpNw9hERJCiksNv+uIgbC4yAsrsYtk+sCStgViopiNMN1H8Os4bD6TbhqZqBHVHNomghjWIxMtpYHt0Oial907XZIVym3U6pEvPlRue6lSOStGeV8BqNE+EYLGEwS+ZusEvmbw+QRiPUE9Qj101EoKkNkM5iyWCbMFEUxWeVRg/1eFWVTPxNICkVdIK6jVMZkp8C2+YEejULhRwm7QlFVVr0G30yF3T8EeiQKBaCEXaGoOhc9Cc17wzfT4PSuQI9GoVDCrlBUGXMo3DBXFqN8PlFqphWKAKKEXaGoDiKbw/WfyhL0n54K9GgUDRxVFaNQVBet+ou4N+8d6JEoGjgqYlcoqpPOo2RVo8cFB1cEejSKBooSdoWiJlj5CsweJ02xFYpaRgm7QlETDLpX7Aa+vFWMthSKWkQJu0JRE1gjYOLnUjEz9zrxSlcoagkl7ApFTRHVCm78Quxiv75dDLUUilpAVcUoFDVJ814wYY7Yw6qWdIpaQkXsCkVN02kkNDtXnicsV5G7osZRwq5Q1BZ7f4TZV0rFjEJRgyhhVyhqi46XQM8J8NsLsGVuoEejCGJUjl2hqC0MBmmCnZ0IC++Xjj6dRwV6VIogREXsCkVtYrLIZGrTHlLjbk8K9IgUQYiK2BWK2iYkEiZ9DcfWi/2AQlHNqIhdoQgEtnjoOkaeH/gdUhMCOx5FUKGEXaEIJM4c+PYu+GQcZBwL9GgUQYISdoUikFjC4MbPIS8dPrlSWQ8oqgUl7ApFoGneCyZ9CVknYc541YFJUWWUsCsUdYHWA2HiPEg5AJs/CfRoFPUcVRWjUNQV2l8I05ZB43MCOw5FvUdF7ApFXaJJNzELS02AL26CvIxAj0hRD1HCrlDURZL3w54lMOcqyE0L9GgU9YwqCbumaf/TNG23pmnbNE37VtO0qOoamELRoOl8qaxQPblNWuypCVVFBahqxP4z0EPX9XOBvcDfqz4khUIBQNfRcMNnkLgbPrkCslMCPSJFPaFKwq7r+k+6rrvzX64FWlZ9SAqFwk/nS6XOPTRafGYUinJQnVUxtwFfVOPxFAoFQIeLof1FMqnqzJZWe1GtAz0qRR3mrBG7pmm/aJq2o4THuEL7PAm4gVJNpjVNm6Zp2kZN0zYmJSlHO4WiQvja6i24Fz64VNIzCkUpaHoV23RpmjYFuBMYoet6Tnk+07dvX33jxo1VOq9C0SA5vUtWp3qccOOX0KpfoEekqEU0Tduk63rfs+1X1aqYy4DHgSvLK+oKhaIKNOkGt/0IIVHSZm/fz4EekaIOUtWqmLeBCOBnTdO2apo2sxrGpFAoyiKmHdz+E8R2hO8fAldeoEekqMlP+LAAAAquSURBVGNUafJU1/WO1TUQhUJRAWyNYcoPkHUKzCHg9Uoe3peLVzRolFeMQlFfCYmUB8DPT4PTDqNfAaP6t27oKEsBhaK+o+tgtMCmj2HeDZCXGegRKQKMEnaFor6jaTDyWRj7Ohz4DT4cBelHAj0qRQBRwq5QBAt9b4WbvoaM4/DRGHA7Aj0iRYBQyTiFIpjocBHc8Quk7AOTNdCjUQQIFbErFMFGfGfoOkaeb/0Mlv4/8LjL/owiqFDCrlAEM4l/wdrp8OnVyh2yAaGEXaEIZi59AcZNhyNrYdZwOL4p0CNS1AJK2BWKYKfXTXDbUnn+4eWqYqYBoCZPFYqGQIveMG057PmhwPLX6wWDiu2CEfVbVSgaCuGx0HuyPD+2Ed4dKm6RiqBDCbtC0RBxO8CeCO9dDBs/lNWriqBBCbtC0RBpewHc9Qe0GQSLHob5N6uG2UGEEnaFoqES0QQmfQ2XvAB7lkrNuyIoUJOnCkVDxmCACx6ATpdAXGfZdmoHxHYAc2hgx6aoNCpiVygU0PgcMBjBmQNzroJ3Vc17fUYJu0KhKMASBlfNBEcWvD8SfvmH6tBUD1HCrlAoitJxBNy7Fs6/Ef54Dd4dpiZW6xkqx65QKM4kpJFYEXS7ShY1hUbLdrWoqV6gfkMKhaJ0Oo2Esa9JM4/UBJjeD/76XtW913GUsCsUivLhygVTCHxxE3x2PaQeDPSIFKWghF2hUJSPJt3Fb2bUv+HwKpgxEJb/L9CjUpSAEnaFQlF+jCYYdC/ctwG6XA7ZiYEekaIElLArFIqKE9kcrvsYLntRXh9eI31Wj28O6LAUghJ2hUJReQxG+ZqTDEm74b2L4JtpkH40sONq4ChhVygUVeecK+CBLTDkYdj5HbzVB5b9N9CjarAoYVcoFNVDSCSM/Afcvwl6XgeWcNnu9UBeZiBH1uBQwq5QKKqXqFYwfjoMvk9eb/8K3jgXVr4CDntgx9ZAUMKuUChqlibdoGV/+PV5eL0nrHgZ8jICPaqgRgm7QqGoWZr2hEnz4Y5foWU/+O0FmHtdoEcV1CivGIVCUTu07CsCf/JPcY8Eyb0v+w/0nwox7QM7viBCRewKhaJ2aXYetB0iz4+thw3vw5u94fNJcHi18qGpBpSwKxSKwNFxJDy4DYY+IjYFH10Osy5Uk6xVRAm7QqEILJHNYMQz8PAuGPOq5OStNnlv25eQciCw46uHVCnHrmnaC8A4wAskAlN0XT9RHQNTKBQNDEsY9Lu94HVeJiy8H9y50G449LkFuo4FkzVwY6wnVDVi/5+u6+fqun4+sAh4phrGpFAoFLLg6YEtcPFTYhH81W3wShfY+2OgR1bnqZKw67peeDlZOKBmPRQKRfUR2QyG/Q0e3Ao3fQMdLoa4TvJewnJY+arypSmBKpc7apr2L2AykAFcVOURKRQKRXEMRunF2nFEwbaDy2U166/PQauB0OMa6DYOIpoEbpx1BE0/S2mRpmm/AE1LeOtJXdcXFNrv70CIruvPlnKcacA0gNatW/c5fPhwpQetUCgUgKRotn8FO7+BxF0Q2wnu3yjv5WVI79YgQtO0Tbqu9z3rfmcT9gqcsDWwWNf1Hmfbt2/fvvrGjRur5bwKhUIBQOJuyDoJHS4Cjwte7gRRbaDrGGkK0qSH9G6tx5RX2KuUY9c0rVOhl+OA3VU5nkKhUFSaxl1F1AE8Thj8ABgt8Pu/YeYQeK077FpQ9jGChKrm2F/UNK0LUu54GLir6kNSKBSKKmIJl0VPQx+BrNOw70fY9xPY8vPvB1eK4He4GNpfCM17Sdu/IKHaUjEVQaViFApFQNn/C/z6gvjWoIPFBm0Gw/h3IDwu0KMrlfKmYoLnEqVQKBTlpeNIeWSnwKGV8ji2EUKi5P2fn4UTm8VuuNUAMTALiwnsmCuAEnaFQtFwCY+F7uPlUWR7nKx8/eM10D2yre1QmLJIniftgcgWBdYHdQwl7AqFQlGcwffLw5kNxzdJNG8oJJefXgMZxyC2IzTtIf42bYdBq36BG3MhlLArFApFaVjCod0wefjQdRj9MpzcCie3ifDv/Bb63SHC7nHDJ1dAbHuI6wLxXeQCENUajOZaGbYSdoVCoagImgZdLpOHj7wMcOXlP08HzSCeNls+LdhnxDMw9NFaGaISdoVCofj/7d1diFR1HMbx71NiYpi9bIHkyyYotBiUSNlNLxghe+FeFGEgvbAEGnTTVeBN1FUXdREI5UX0ApXVRQxYN5WyIK0laGpCtatSW5La1t5IZfTr4hxoWNad4+55mXPm+cDAOXP+7Dy/OTO/PW8zM1+Llv7/Kder++DJvcn0hUk4/wNMjsOy20uL48ZuZlaUxdfDyruSW4n8QxtmZg3jxm5m1jBu7GZmDePGbmbWMG7sZmYN48ZuZtYwbuxmZg3jxm5m1jCVfB+7pHMkP8wxF33A+Rzj1IFr7g2uuTfMp+ZVEXFjp0GVNPb5kHQoyxfNN4lr7g2uuTeUUbMPxZiZNYwbu5lZw9Sxse+uOkAFXHNvcM29ofCaa3eM3czMZlfHLXYzM5tF1zZ2SZslfSdpTNJzMyy/StKedPlBSf3lp8xXhpqflXRC0lFJn0taVUXOPHWquW3cQ5JCUq2voMhSr6RH0vX8raR3y86Ytwyv65WS9kk6nL62B6vImSdJb0g6K+n4JZZL0qvpc3JU0vpcA0RE192AK4FxYDWwEPgGGJg25mngtXR6K7Cn6twl1Hw/sDid3tELNafjlgAjwCiwoercBa/jNcBh4Lp0/qaqc5dQ825gRzo9AJyuOncOdd8DrAeOX2L5IPApIGAjcDDPx+/WLfY7gbGIOBkRfwPvA0PTxgwBb6XTHwGbJKnEjHnrWHNE7IuIC+nsKLC85Ix5y7KeAV4EXgL+LDNcAbLU+xSwKyJ+B4iIsyVnzFuWmgO4Jp1eCvxSYr5CRMQIMDnLkCHg7UiMAtdKWpbX43drY78Z+KltfiK9b8YxEfEPMAXcUEq6YmSpud0wyX/8OutYc7qLuiIi9pYZrCBZ1vFaYK2kA5JGJW2m3rLU/DywTdIE8AnwTDnRKnW57/fL4t88rSFJ24ANwL1VZymSpCuAV4AnKo5SpgUkh2PuI9kjG5F0W0T8UWmqYj0KvBkRL0u6G3hH0rqI+LfqYHXVrVvsPwMr2uaXp/fNOEbSApJduN9KSVeMLDUj6QFgJ7AlIv4qKVtROtW8BFgH7Jd0muRYZKvGJ1CzrOMJoBURFyPiFPA9SaOvqyw1DwMfAETEl8Aiku9TabJM7/e56tbG/jWwRtItkhaSnBxtTRvTAh5Ppx8Gvoj0rERNdaxZ0h3A6yRNve7HXqFDzRExFRF9EdEfEf0k5xW2RMShauLOW5bX9cckW+tI6iM5NHOyzJA5y1Lzj8AmAEm3kjT2c6WmLF8LeCy9OmYjMBURZ3L761WfPZ7lrPIgydbKOLAzve8Fkjc2JCv/Q2AM+ApYXXXmEmr+DPgVOJLeWlVnLrrmaWP3U+OrYjKuY5EcfjoBHAO2Vp25hJoHgAMkV8wcAR6sOnMONb8HnAEukuyFDQPbge1t63lX+pwcy/t17U+empk1TLceijEzszlyYzczaxg3djOzhnFjNzNrGDd2M7OGcWM3M2sYN3Yzs4ZxYzcza5j/AGxxeuhsMNRWAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"def plot_gp(x, mu, var, color='k'):\n",
" plt.plot(x, mu, color=color, lw=2)\n",
" plt.plot(x, mu + 2*np.sqrt(var), '--', color=color)\n",
" plt.plot(x, mu - 2*np.sqrt(var), '--', color=color)\n",
"\n",
"def plot(m):\n",
" xtest = np.linspace(0, 1, 100)[:,None]\n",
" line, = plt.plot(X1, Y1, 'x', mew=2)\n",
" mu, var = m.predict_f(np.hstack((xtest, np.zeros_like(xtest))))\n",
" plot_gp(xtest, mu, var, line.get_color())\n",
"\n",
" line, = plt.plot(X2, Y2, 'x', mew=2)\n",
" mu, var = m.predict_f(np.hstack((xtest, np.ones_like(xtest))))\n",
" plot_gp(xtest, mu, var, line.get_color())\n",
"\n",
"plot(m)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"From the plots we see:\n",
"\n",
" - The first function (blue) has low posterior variance everywhere because there are so many observations, and the noise variance is small. \n",
" - The second function (orange) has higher posterior variance near the data, because the data are more noisy, and very high posterior variance where there are no observations (x > 0.5). \n",
" - The model has done a reasonable job of estimating the noise variance and lengthscales.\n",
" - The model has done a poor job of estimating the correlation between functions: at x>0.5, the upturn of the blue function is not reflected in the orange function.\n",
" \n",
"To see why this is the case, we'll inspect the parameters of the kernel:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>class</th>\n",
" <th>prior</th>\n",
" <th>transform</th>\n",
" <th>trainable</th>\n",
" <th>shape</th>\n",
" <th>fixed_shape</th>\n",
" <th>value</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>VGP/kern/kernels/0/lengthscales</th>\n",
" <td>Parameter</td>\n",
" <td>None</td>\n",
" <td>+ve</td>\n",
" <td>True</td>\n",
" <td>()</td>\n",
" <td>True</td>\n",
" <td>0.6668020738578904</td>\n",
" </tr>\n",
" <tr>\n",
" <th>VGP/kern/kernels/0/variance</th>\n",
" <td>Parameter</td>\n",
" <td>None</td>\n",
" <td>+ve</td>\n",
" <td>True</td>\n",
" <td>()</td>\n",
" <td>True</td>\n",
" <td>1.3274634849562423</td>\n",
" </tr>\n",
" <tr>\n",
" <th>VGP/kern/kernels/1/kappa</th>\n",
" <td>Parameter</td>\n",
" <td>None</td>\n",
" <td>+ve</td>\n",
" <td>True</td>\n",
" <td>(2,)</td>\n",
" <td>True</td>\n",
" <td>[0.8848631418180702, 1.3974074371796048]</td>\n",
" </tr>\n",
" <tr>\n",
" <th>VGP/kern/kernels/1/W</th>\n",
" <td>Parameter</td>\n",
" <td>None</td>\n",
" <td>(none)</td>\n",
" <td>True</td>\n",
" <td>(2, 1)</td>\n",
" <td>True</td>\n",
" <td>[[0.0], [0.0]]</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" class prior transform trainable shape \\\n",
"VGP/kern/kernels/0/lengthscales Parameter None +ve True () \n",
"VGP/kern/kernels/0/variance Parameter None +ve True () \n",
"VGP/kern/kernels/1/kappa Parameter None +ve True (2,) \n",
"VGP/kern/kernels/1/W Parameter None (none) True (2, 1) \n",
"\n",
" fixed_shape \\\n",
"VGP/kern/kernels/0/lengthscales True \n",
"VGP/kern/kernels/0/variance True \n",
"VGP/kern/kernels/1/kappa True \n",
"VGP/kern/kernels/1/W True \n",
"\n",
" value \n",
"VGP/kern/kernels/0/lengthscales 0.6668020738578904 \n",
"VGP/kern/kernels/0/variance 1.3274634849562423 \n",
"VGP/kern/kernels/1/kappa [0.8848631418180702, 1.3974074371796048] \n",
"VGP/kern/kernels/1/W [[0.0], [0.0]] "
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m.kern.as_pandas_table()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see that the W matrix has entries of zero. This is a caused by a saddle point in the objective: re-initializing the matrix to random entries should give a better result."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"m.kern.kernels[1].W = np.random.randn(2, 1)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"lib/python3.5/site-packages/tensorflow/python/ops/gradients_impl.py:100: UserWarning: Converting sparse IndexedSlices to a dense Tensor of unknown shape. This may consume a large amount of memory.\n",
" \"Converting sparse IndexedSlices to a dense Tensor of unknown shape. \"\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"INFO:tensorflow:Optimization terminated with:\n",
" Message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'\n",
" Objective function value: -148.744827\n",
" Number of iterations: 1495\n",
" Number of functions evaluations: 1626\n"
]
}
],
"source": [
"gpflow.train.ScipyOptimizer().minimize(m, maxiter=notebook_niter(2000))"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXd4FNX6xz+zNWXTCx1C6CAoRZpUQUBAigqKVEVRwXqt91qv7f7sioIoCgIiiA1FUJogvSMiNfQikE6yKVvn98fJ7qZsIGXTz+d59kkyMztz0r5z5j3v+30VVVWRSCQSSfVBU9EDkEgkEolvkcIukUgk1Qwp7BKJRFLNkMIukUgk1Qwp7BKJRFLNkMIukUgk1Qwp7BKJRFLNkMIukUgk1Qwp7BKJRFLN0FXERSMjI9WYmJiKuLREIpFUWXbv3p2oqmrU1Y6rEGGPiYlh165dFXFpiUQiqbIoinK6KMfJUIxEIpFUM6SwSyQSSTVDCrtEIpFUM6SwSyQSSTWj5gr7jtlgTvB8bU4Q2yQSiaSKUyFZMRXOjtmw4knY+TlM/EVsmzcUEg6LzzvfV3Fjk0gkklJSM2fsrUdAVEsh5DO7ilfCYbGt9YiKHl3FIJ9gJJJqQ80UdlOUmKkHREJmongFRIptpqvm/lc/XE8w84YKQTcniM9XPCnFXSKpgtTMUIwkL61HiLCU6wkGxM2uJj/BSCRVmJo5Y3fNSF0zddfM3TVjrWnIJxiJpFpRM4X94FJPTH3qNvFyxdwPLq3o0UkkEkmpqJnC3vk+GPyOZ0bqmrEOfqdoGTHVbaFRPsFIJNWKminsIAQ8d5jBFFV0US/qQmNVuQHIJxiJpFohF0+LS1EXGqtSrrxrLK1HeG52E38Rol6ZximRSIqEoqpquV+0U6dOapW27TUnCFHPTBRfB0SKWW7uJwDXTD7hsNgPnhuAXJSUSCQlQFGU3aqqdrracTU3FFPWyEwTiURSQUhhLy5yoVEikVRypLAXl6IuNMobgEQiqSCksBeXoqZKliTTpIyzaOZvPUWi2eL+OtFsYf7WUz47v0QiqRzIrJiSkD9TxFuqZHEzTbxl0cy6AcyXPOczJ8DBpcx33MTgtnWINBkBIdAr9l9gQreYQoc8f+spXvzpAAu2nmbRFJHNM+azbcTFmwGu+F6JRFK1kMJelhTlBuDCkg4RTT1plLZM8dIaoFHPPFk2cbZJjNl6KyM71KN/q1pMW7iHuHgzZosdk1HnVaQHt63Dgq2niYs3M/D9DQAkZVhpFm1icNs6Pv7GJRJJReKTUIyiKHMURYlXFOVvX5yvxrFjNqz9r/jcP1zE4m2Z4muHlawvBuGcIayF7REt+DukD3HxZt767Qg3f7CRuHgzsZGBfLvrLC/+dMATXskV2ok0GVk8tgkP+K8lKcNKUoaViEADIzvUAzxhGld4xvVRhm8kkqqHT/LYFUXpBZiB+aqqXnO14yttHvuO2XlDJzmhjyIX6TjsYM8Gpx1UJzgdYrt/GGh1YM0U+zVa0PmDVg+KkjfnHQXI+Z0oWrK0Afjb0wFw+keSMukPRi04xonEDPdlFSA0QE9Kpo1m0SYWTelK5MH5sOJJUgNjcUxYRkSgEfvcIeiSjvCCbRILHAMINGjJsDqICjKSkG6hcWQACgonEjOIMhlIMFsB3OdMzrAydvY2EsxWXhneRoZvJJJypqh57D4JxaiqukFRlBhfnKvCuFKlqMMO3R6E83tg7wJI+wcyEiE7FbIvw4SfoVZr2D1XnCM/j+wFQxAsexSOLPdsV7RgCISHdsLt8+CTboDTs1914Gf3CHhqlpU7Pt3GiUx/NAo4VRivXcUKRxeSMkOICDQwa2RDln/+Cms03XhV04CYjBMkzexGuk5DkCOVo856rHB0ASDD6kABEtItKMDJRPGUoIBb1HUaiIs3c9N7f5CWZcehqkSZDDJ8I5FUYmSM3UXrEbD9UyHkH7QVs26nTewLqi0+pp2HA0shuB4ERkJoA/ALAWOQ2N+gC9z0CpzbBYd+BlMt6PKgmIC7bhK12kGrIWLWbssSN4i/fxAz+dyinoOCEytGrKpCuJLGp/bnGM1LJKkhjNeu4lX9l4zXrmaM9XkM2aDMf4yJ6jmO2cyMdv6HXw3PEKGkgQMS1WDGWJ8nRQlxPxS4ntdyP7fl/tyeM6SUTPGz0CoKn03o5F64zc38raeKvagrkUh8T7kJu6IoU4ApAA0bNizbi10tpGK3QtIxuPgXXNgHTftD034w+G1YMALsWeI4rRG6Pwx1rxVftxwKrW4p/Joth0KnydBiiDh3yinY/KF4ZaeIXPZLf4nzT/hZhGLmDYU988RNALApfujVbPdpnTp/LmsiibKexaLqaKK5yHLjf/jG0ZftjpYcddalueY8K43PABCppuWZlXsLtDlLEX3TaVRu+2QLbeuFMPCa2vRoGsmfZ1MBZNaNRFJJ8JlXTE4o5pcyj7GnnhXhi/3fQZuRBcUbRDgkohkM+wjM8bD6BUg9DQNehz3zIfm4mJGDiHX3ewG6TYM/3oYt08GSJvb5hUGne6D5AEg+CWnnIP2ieGUmQVZOKCYrOWfGXUIULagOrNpAllo60E/3FxFcBiBDNRCoWDnqrMf7tlsZodtCqJLB9coRNIrKGkd7rtMcJ1IRY05UgxmlvEuAUcv7WS/QXHOeRDUYgEhFiP5d1udJJOTqw1KgsD+PsJyYvotJ3Rux+VgScfFmIgINgCfrZtGUrl5n+BKJpHiUa4y9XFk4ChIOic9/exZM0dDoBrj0d47hVoTYlxQHcwd53hfVEtrdIWbSrYZCVCuo3RZOrIdrboMN78K618SxilaoWnYKbHpXvEqNAgHhWKw2DPZ0lNxzaVUsshocGYzWbczzrkDFik3VssXZhizFj6dsD5BGIBFc5ibtbuyqlus0x93Ha3Hgp9cyQN1Gc815jjrrMcb6PAqw2Pg6zTXnuN1/F5YOk1m47TRWh2ccUSYDVofK5SwbOo2CPWdqrwH8DFoyrQ73sblFXatR6N+qFtHBfsxaf5ykDBGfDwvQM7JDPSnqEkk546usmEVAHyASuAS8pKrqF4UdX6oZ+6FlEH8Itn0iZspaI2h0YMsQ4t2kPzitsHeh2AZi5v3QzrzmWw67SDHcMl08ATjtYLd4vyZAnesgtg8E1REx9zPboOUQCG8MBhOknIU5N3nCOLnJmZGnBsZye9L9fOv/BmGqmJGnEsS9lsfx16mEOVKI0SYSo55jqHYrBsVR4FROVeFPtQmrHJ3Y6WzB//Sf01xznstqAMFKJgqQhZF3jNOwmpNY4ehKqhKKQ1XpEG6jt30z76f1cWe9BBq0GHQaUjJtNI4M4NPxndh2Igmzxc5bvx0h0mTg6/u6Eh5ooPdb68iwFhwT5MnlcaNBrBrIDBqJxDcUdcZedW17C7POPbhUFPTMG+LZZwiEh/+ErCRY9z+wposFTlfI5Wq4rHYDIyErBbbOgI3vQGQLmLQczu+Gn6aK8Ex+/MOh+2Ow+V3Ivkw2BvywkoWBTNVIhJLOKWc0rzsmEOesg0XVM8/wZk4IJQgNKuGKmWTVxCVtbZqrp9Cq9jyXSFBDuNPyHHf3as7gnRMJz7lpHNM0Zn33ufRq29RdxPT0oBYcOH+Z5fsvusMk4ImH5xZhb4uhX2w6yfwtp9wCb9Aq9GtVi1UHLuLw8qek0ygsf6QHLWoHF+1nLZFICqVmCnu3aWIWrjWAwwp+oSIGjgqKRuSW5ya0kRDqKwm8fzj0fkaIefoFT6aMi4BIsJrzxNhVxAwWIEsfQlbrO/E79D0BznSwW3CqClZ03GCZziLDazTXnHe/16mCRoEM1Ugvy/uAwreG/xKruciHhimMfeBZ/t7wI+0zt+J3YhVGuxi7U9GhaTWEy81uI2XLPOLDO9LOcB6/22aCOYHMde/wffi9jL+hGVCyDJZEs8V9A4gINOBwqqRmidz5lnWCWLbvgtf3tW8QwuyJ18tsGYmklFRvYS+siUV4LGQkgeWy9/cZTELgLWkiPKPR5twYXEEDLyhaGPwenNkMIfVE9kpglAj//PIv8RTgPlZHml9tgrPO5RH3OLU+d1qeIzy6HoMyf+Zrc0eSCAYUIrjM3Ya1HNK1wpidQH0lkT7avZxXo9ittmJvUG9eyHidazUn2OtoSt1a0dQf8TLU7SDCR0d/g71fwbHVeW9cQbVh8lrQGWHOIEg+Jm50fZ+DDuPFe4pZjOXym/E207+3Z2NWH7zE6aRMr+8ND9Dzv9va0bZeCBPn7CAu3syQtrWZMbZjodeTSCR5qd7C7iomimwhXBVPb4YtH3li6rnR6KHrVFH5eeMLIm98RhePIOfEvwtgCAJrBuAUYv7A5oIdknI/MeSKMttVDTrFSYIziDRMNNFccFd7esO1UNks2uT2f1n66cs87fyc5IBYDsTeQyflMPoDS9A5xcIkrYbBHQvE56oqniZ2zRUFVlnJOUPSgM5P2BME1xcLzf/sETc4q1lkDt39qzjWdaO8SkPvwmb6g9vWcYu8UafBYvd+o3QVVuk1CjanyivD2zC4bR05g5dIikD1zYpJOS1moXXaw+UzMD9fXnl0a5Ga6FrE9AsRueimKEg4KqpDs1M8x6sOCG8Ct34m0h3jfhPbe/4Lmg+GBcOEw2Lu2Wwur/VENRitRiFMvYyqaFFUBzrFSaIazCDrmwAM1m4vIOoK0K9VNIcvpnMuJYshbWvz3+HXuAVz8v2Pkzr/d8IzTtDzxAfiTU6rEOMu94swEkBmMszuK1I/rxsLPZ+A3fNg5b/F9+bynGnaD3o8Lqpm178BpzaJzKEP2wnxz0oWJmTpF8WNsxBxzy++kSYjE7rFMH/rKeLize7Z/IXULB74ajfnU7Pz2BO4cuhtTpUGYX50jY2Q+e4SiY+pOjN2pxM+6e5JdXQRVBea9IXYvlCrLXw3sWCIJrKFCE2c/MPzPoNJxOJdgnb3b0L8N74n9vf8l/joLUSR88RgCWvOONvznEjM4Hu/V4jBE2NOIpgB2W+S5CVfPMhPx8J7u9CufuiV481F6a2adBx++zccWyOEvHEvuGYUrH4JspPznk/RQOvh0H4ibP4ATq737NMaxYz+8lnx9ZVm7oUUgHmzE37pp79Zvv+i9/MAwX460rLtNIs2MWNsB7adSJLiLpEUQvWbsWs0QkjS/4GYntC4N8T2hsjmIucchOC4mlvcNgeOrID930LiERF/bzkUDv/iyXIBTwjCJd4uQXfhxWp3vuMm4myTOGy7kdfH38gDs37D4QA0okBIo0AEaSwyvMYY6/MFxD3KZKRuqD/gmfGWmIgmMHYJpF8SPja75sCyh8W+gMgc47FMIeqKAgd+FK/8OCweUQ+NEcLtjSt46kwY/A6YPD+rSJORR/s3Z/XBeKwO76GZtGw7of46nr25JVO/2s2xBBFOk+IukZScqiPsALd9IQqQNFrv+zvfBxkJYhY7uw84bEJ8Eg55xDn/bPNKzS8KYXDbOty8digJiVbu+HQrQ62baKK/4C4GAtzZLvnDMBoFTiRmXD2mnL+1Hnha63lriB1UC3o9KXxrfn3a8xSy9AE4tlYsrOr8xTqxt1z73KT/I1I4WwwquK/1CCHqLt9417iiWnq9GWw7kVSoqLtIzbIzeZ54gmsY5i8NxiSSUlK1hN0UXfi+zGRY9TzsWyxCLB0miEXTiCbQuIfnuOI0vyiEFfsvkGC2oiAqMBcghHuFowtJhKAAY6zPM9K4k2WGQZBldy+tOlWK5o6Yu7VeYU8X3lByWeyrKoTk8uWxuTxw/MCR7X3hWGsUlbgNOouvE46Ioiy/nDx0VyvA/CEibzcbxMx7w9F41hzytPzzVszk4lxqFh//fozHb2pOiL++kKMkEsmVqPo9T505s0G9P5zeAl0fhMf+giHvClEvAwa2qU1sZGAecVrgGEAywjpXBVo1jSW87zRSs0T8eOXjvXh6UAuaRZtIMFtZsd97zrebovZWzU/rEeJmkHRM2ADv/4YCMurIqbA1evGLcVjgxO+iMldV4Ycp8FEH2P2lx1++GCSaLZxOEjeUsAA9gQYtKqDXKtx9QyOMOiXP8QadhnlbT3HTe3+Qnm0reEKJRHJVqs7iaX5s2bB9Fhz4ASavFpkyDptwTCxDDl1I45FFe0kyW0jOzCs8Wo1C3RAjZ1Oy3RWcFWJlWyAV0wsarUeogxvAtWNEdaxrW1iMWJNo0h82vAlntoobRsPucGZLwQVq15NFvlm7t9z3mz/cSEK6xZ0tY9AqeTxrAFrVCWL+PV2ICjKSabUTYKhaD5cSSVlQvfPYj/8Oy5+A5BPQfBAM+9hrGMCXOJwqn288wTurjqCquA2yvFHhjoZehT1XEZYrBGOqldMsWxGhq1YjYcnYvPUABhOM/xni/4bl//K4YnoLERXyNOHt5ubKlnH9rDbFJfLEkn04VBWtouBQVUL89Uzq3oj5W0/z8rA2DL+unq9/UhJJlaJ6CrslHX5+RMzSw5vAkHegyY2+H2A+Es0WHv56L1tPJBGQz+Uwf7w4LEDP6n/1rlhR91aVmzue7h8ucvvbj4cvboKUk3nPEdEMYnqI8Ivru9MHiJz4wCjofD90nCRm/QHhxW8hmEN+wT8Wn84TS/bx3JBWfLzuOBuOirh8aICe1Ewbt3WozyvD2xBolLN3Sc2kqMJetWLs+kDhr973OZi61aeifqWmzcv2/cPOU8loFci0OqgVLIQoOsjItw90IyzAE/6xFlJxWW7kXnSdug0mLheLyapD/PwCIkXu/l/fiOMnrxbbXegDYfRXEN0K7lsnFlNBiLqigSHvQe+nhKjP7Cba/en9iy3qIBZWc98Am0YH8dNDPejcOIIvJ3Xi2vohGHUaUjNt+Os1/LDnHMM+3sThi0U0b5NIaihVa+qj0cDEZeKjD3HFgT9aG8fCHItaVzXk5rgEwONNHmDQ8vGYDhy6mEbX2AimLdxDSqaNsAA9VruTDKuDMZ9tq7hQjEtgXSmdB5cKQzSXLQLkza6xpOcNveiMIhyTdAwCo3PWLnJueKoTlowXXaL6/gfajRZWDifWi7qB+r7zfcm0OQk06rDYnYT660nNEusZ51OzeOu3I7x1eztpKiaRFELVCsWUEYlmCzd/sIEEsxWtohDkpyU1K681rmuBr2XtIGaM7UCTKNMVTbEqlQd5Ya0CW48QzUiSjuUcmCuwZKolbIiddjFr1wfktWIIaQgtB0NML/jtGWFF0OdZMAZDlyk+GbbTqbJg22n+9+shsa7hUHHk/L2G+ut48/ZrubZ+COO/kKZikppB9YyxlyFHL6Vz8wcb3cKRG3+9hiybkzGdG/DSLW3w03sKpKp0A2dXFWlEU1EHkJXLgqD5YDi6QhSEZaeB0STEPe28aOadlmM1HBAhFlFXPQdntwvjNNci6tV6zxaRk4kZPPXtPo4nmAnx13Mql4OkVgMOJxi0GqwOZ+W6oUokPqb6WQqUMeGBBq8z9SZRgZxNyeKt29sxulODAu8rzBSrSuASWFdjEheGQBg2HQ7eKIQ5I0HE0rtOFQuxHSfB6pdh28diVv9pT9CbhKibakGDrqKhyR//59V6IM+1i0DjyEC+ub8bZ5MzqR3ix7+W/MmKHP8ZV1Gr1eGkWbRJVq1KJEhhBzwNJPKLugJMH9Meo05D0+igihlcWdN6xNWtC0xRMHmVx5Nny0ce10gQ4RpLKqAR6ZPf3yMWuUMaFNl64GpoNQoxkWKRt32DMLew5+auLg1lf1WJhKqWFVNGrNh/gbh4M0q+7Sowae5OQgMMFTGs8iF/Fs3UbeJz1+KqC5eo263w1xLYPUdYIufBCf5hcMtHYuaefkHkwWcmem4chVgPFJVEs4Vvdp7xuu/r7Wf4cvOJQrObJJKaQo0X9jNJmaw7HE9UkCFPPvqNLaOJNBlISLdcvfy/KlNc6wKdAe5dI8Ix2V46VWWlwLpXYdx30LCraOjhQ1bsv8CxhAwMWvGnq8l1N46LN/PyskPcOnMzJxPN7iexF386IMVdUqOo0YunxxPMjPlsG8kZVuxO1d3d57H+zXi0XzOSMqxVZyG0vHEtvOYhV3VreBORHplyUhRH+YWKrlWFWA8Uh2kLd7urVufefT1vrzzC6oOX8hSOaRQI9tO7e7JWaCWwROIjZFaMF+ZvPYXZYmd0pwYkma3c+dlWLmfZcKrgr9eSbXPw+si23NWl4VXPVeMxJ8CsG0RM3RAk8t0zE8E/AgIjIPGoOC64AYz/QYRovhwivPGvZmRWBLxlI83dfJL5W06TbvGslZiMWtY/1ReAJbvOYjLq5I1aUmWRWTH5cOWcAyzcdpq0bDvp2R4BmNS9EV1iI+jT4grWwBIPB5cKUc/tGfNZb5EG2fl++Gc3xK2EjIsQf0AsmIY3FsVl7ccV71pe0iYnaFcXaOrx8I3N+P1wPIcupLu3my0OZv1xnLWHLnEy0bPgK8VdUp2pMTH2wW3r0DgyAIDzqdl5RL1eqB+Te8ZKUS8O3mLzNzwmrAb2L4abXoHOU4Tj5rd3w/ZPodkAiD8MC0eJitei4Ar5zBsqnhJcXjgrnhT7cmG22LF5sXT4fONJt6jHRgbKlEhJtafGCHukycgbI9sS7Few+9Ibt7aV8deS0Pm+vLHyLlNg0gpR0DRnALQcAv1eAlRRnZp6BkZ+KnzzF4yErNSrX8PlL+9Km5zZ1ZPFky9t0rWw2izaxHODWxU4lZ9O4d3R1+YJ38hFVUl1pMYI+/oj8Uycs5O07LzNIoKMWtrU9dJwQlIyGnaF+34XTca/uk3YKo+YBRqdaKB9ehOM+hL++RO+GSeaeVwJV5ZOQORV0yYndIvhleFtWDSlKyM71CPYL2+kMduuMmrWVvaeSZEZM5JqTY2Isa86cJGpC/d49VBPtzgYPWsrSx7oJmftviKskSho+vt74RJZq7Ww910yAfbMFzP10fMhMNKTH+8jJnSLIdFsYdSsLaRl2wvstztVRs7cgsmoxWxx0CzaRNfYCOZvPSXj7pJqQ7WfsS/b9w8PLtyTJ0ddr1WYPuY6d8zd1Vxa4kP8gqHT3UK4L/4tGmqP/U604zv0M+z4FGq1Ecfungdp/3g/T/6m3q6Zuyvm7oUV+y/kiamverwXQflCcGaLgwC9hhljOzBt4R45c5dUK6p1umNKhpUb3vwdu1PFandSJ9jI7Z0aMK5rI2oF+5FotsgUuPJgy0ei0XizAdDrKVg8FjLioUEX0f1q9o1i9n7dGOh4T17TsBVPwcEfvXZsMvf7H6aeU4GC5mu5U1sjTUYSzRYe/2YvG+OS8gwtyKgh3SJ8Zga0qUWgUcfUPk29nlMiqWhkuiOQlm0jQK8lMcOKQavh5rZ1eWJAC/f+SJPR/U8sKUO6PyyMxZY/ISpTxyyGbycKN8jv74FRc0W8fd0bsG8x3LNKvM9lGtZ6JAx+2y34S9rMZP/q+Wzb2YZF7YV9gMsuGUQ4xpsYX7xsKbAt3eIkyKihV/NIZqw77t4+ulODAueUSEqFaxLt4/CjN6rljH3JzrOcT81i0Y4zxKdbMOg0GLUaFtzbhesahJbZdSVX4eDP8P29olH2yE/h+8mQfFy04uv9DPw4RVSr+oWKxdZCmmS7Fj7j4s1EBAofn6QM6xUrTHN75z/SrxlPLPkzTwPt3C0OwwL0aBTlqueUSIqM3SImNqZa0O+FEp+merbGKwJfbz/D09//xcz1x4hPt2DUafDTafhKinrF03qYqEKNaiEWVO/5DWpdA0lx8Pur0G4MoEB2qhB1fQC0HJo3NLNjNpEmI4umdCUi0EBShpWkDCsRgYYrCrArY2bG2A5MXxuH1aFiMnri7i5RDzRoScm0FemcEkmRUFX4ejTsXSBm6+Uwma5WoZg5m07yyi8H0WkUbA4VP50Go17LV5O70La+TGmsFMT0EC8Qs/KB/4M1L8E/eyD1NOT22LRlwsZ3RBin/fi8fu6tJxT70hO6xTB/6yni4s00izbx8V3tGTx9o9vTHSDD6ij0/RJJiVAU8ffbYQJcc1v5XLK6hGI+WX+cN387jE4j+pP2axnNbR3qExsdSMvawT69lsRHfHs3HF0Jt80Wi6Surkz5yReaSbr9e+5ceLxYoZjczN96yt2vNi7ejFYBR75/gxB/HZez7DIUIyk5h5eL5jPtRvvslNU6FDNt4W6OXvKUpB+9lM7X20+jUUSecvsGoXwyriOD29WRol6ZGfQ/EW9fMr5wUQdPaMY/HNrdwfITdvese+XjvVj5eC+aRZuIizcXKW11QrcYtp1Icp9j7RN9qBfqn+eYx/s3L9Y5JRI3qgqb3hfZX7vmgLOgzUVZU+Vm7C7LVq1G4fMJHUnLtvPEkn3u4iOtRiE2MpAVj/ZEr62S963qT25Tr8xk+HKoMAq7KjlLnIPfYb7jplL3ms3tEJlhsXPzhxswWxwkZ1gBGNimFjtPpbDk/m40jTYV+9uU1EDsVvjlMfhzIbS5FUbMBL3/1d9XRKqtbe/RS+nc/OFGHF6qSAHqhPix7OEe8tG5suIy9cqdlz73ZrGAWhQULTywBWq19PnQsqwOtBqFhdtP899lBwFh52zQafhsfEe6xEb4/JqSaoTDJjyQTm0UWV59/u3z1MZyDcUoijJIUZQjiqIcUxTlWV+cszAahgcw7Nq6Xvf56TUyHlrZ8WbqlRTn6bd6NVQHnN5YJkPzNwgRH3FdPWqH+KFRIMvmwOlUGffFdn7eV0h1rEQCoNVDo+5w62zo+59yyVcvjFILu6IoWmAGcDPQGhijKErr0p7XG2nZNgZ+sIEf9xaMx+o0Cmv+1ZuYiMCyuLTEV3gz9fIPB/+cVFS/MNAWcmP2CxXWwBf+LLrtbwkINOq4rn4oTlX8XaVb7AQYdDyyaC9v/Xa4zK4rqaKc2QbndovP+/7Hp4ulJcUXM/bOwDFVVU+oqmoFFgPDfXDeAgT76bk+JtzrPrtT5XxqVllcVlLWOCyQdEzM5B/YBA27eT8uO1U00v5zEcy5GS6fK5PhGHQausSKvzO7U6TNXs6yodMozFx/XHrKSDz8/T3MGwarniuX/PSi4gthrweczfX1uZxtPifJbGHNwUt5tuUuMpn85a7iCKehAAAgAElEQVQ8HeollRBvpl7WDFGRd/s8CK0Pw2eA1pDzBgU0es/70/+BTveInPfZ/eCfvWUyzFuurUvTKPH0l213oiBEXqdR2HUqmfFfbOdypg2Qvu41ElWFje/Bd/dAvQ5w59cVGnrJT7mljSiKMkVRlF2KouxKSPDuynclLmfZuPWTLaRm2XAlu/jrtax8vBff3N8Vk1GH2WKXqWmVnYNLPY0ypm4Tr6iWos2eK3Z+9FdwWIWga7Qw4hPQ54TYdH7QbRrcs1KI/5yb4dQmnw8z0mRk8f3dCA8QNxUVkXFld6r8vO8CG+MSueXjTew6nSx93WsaDhssewTW/heuuR3GLxW21JUIX1Sengca5Pq6fs62PKiq+hnwGYismOJeJMRfT7+WtfjjaDz/pGYTEqRn8ZSu1AsNoF5oAOuf6iOd+KoCribWuXuYTvxFCL5rn+tjbF9Y+qD4p5m6DRYMh+QT4tH3nl/hvrWw5mWo3a7MhqvkmoUFG7WY/PScTREhvzPJmYz6ZCsq0CzaJFvu1RQUDZjjoeeT0Pc50ce3klHqdEdFUXTAUaAfQtB3AnepqlpoYnJJ0x3PJGUy8pPN+Om0/DC1O7WC/Uo6bElVQVU9j7iZSbBwNJzfBRFN4e7fPDcHayas/5+wBfYrfVFaYUZjAEFGHemWvE08Zo3rwKBrpLBXa1LPigro4DrgsIO2/B1Zyi3dUVVVO/AQsBI4BCy5kqiXhsggA72aRfHl3ddLUa8puER95+dipj5qbo5x2LGcvqkpYv+ZLbB1BszuKxpml5IV+y94rW4FsOWrJNQoyAKm6s75PfB5P+FAChUi6sWhyhUoSWooJ9bDgluhyY0wbDrMu0WIe932IpRjNMGh5fDj/cL6d/jHcM2tpbpk7spUwN2Y5ZsdZzmdnJnnWINW4ZNxHfl6+xn+O7wN9cMCSnVtSSXi0DL4Pqdx+13fQrTvi+OKSrX2ipHUQGL7wJB34Nhq2Dwd2t0htv+zV1iipp6D318Ba7roxvTd3WIGXwomdIvJU+wWaTJiMuo4nZxJbGQg3ZtE8M6odhi0GqwOlQe/2sOW40kMeH8DKw9cdL9PZs1UUVyZL9+ME20c711boaJeHCr384REkptO90DCUdj+CQx6C8JiIeUEnN4MH7UXmTRRLWHcD6KnavNBPh+Ca3E+90y+U0w4Yz7bxoXL2ZDj+nv/gt08M6gFt3esz12zt8tOTFURWyb89Y3IfBn+sU89X8oaGYqRVC0cdlh4GzQbCG1HwcedROESiIrVx/6GoGjP8aoqTJmaDYSWg8tkSHGX0hn4wQbCAw0kmq159gUatGRYHdL+typhThChPb2/MKnzD6s0OeoyFCOpXuyYLf7htDoY96MQ9b0LRJaCC4cFNryV933Zl0W4ZvEYWP4k2HxfndysVhBfTLoei81JkJ8OvdYjAhlWh+zEVJW48Bd81ke0sQORausjUf9q22k+33jCJ+e6GlLYJZUflyPkvKFC3DOTYO4gUSCSmQj+EWDMSXHcORvWvuZ5r38oTF4N3R4S+2b3g/hDPh9i3xbR/DitO+GBBpGhmWufze7gPz/s58jFsvO3kfiAAz/CnIGACp2n+Pz0dUP9+PNsKs5CnGl9iQzFSCo/LhuChMMeF8jMRM/+7o9C94fh017CcgBg6AfQ6e6854lbA0sfENWrj+wVbnw+5niCmaHTN5Flc7i7eYFIifTTa5k5tgNnkjNL7SUv8SFOh+i5u+l9aNAF7vgKTNFXf18RScmwEhZouPqBRUCGYiTVB2+OkAGRcOOLooJ1y3SIPwj3bxDNDQB+eVzMwHLTrD88uBVunyNE3emA9IsFr1cKNh9LJMsmYuobnu5Ll8ai1NypCkuCe77cyYs/HWDMZ9tINFvchVDSkqACSb8Iu+ZCx7th4jKfivqeMyn0emtduVudyKwYSdWlwwSxwBV/EL6/VzhDjpoL0a1h3Wsi99gvROS+uzBFeapVt0wXs7Sb3xZWqz6IpebOmlGAv85dJtIkFlXTs+1EBxmJT7cQF29m4PsbAE+/VmlJUM4kHoOIJhBSD6ZuhWDvfR5KyqELaUyas4Nwk4FOMWEAXM60ERLg+yfF/MgZu6Ty480RMjNRbLNlwah5wp9939fi+F5PQtep4LTB4nFwrpCwX6thIj3yxykiV9kc75PhuvLfI0xGPhnXgWybkxB/PQadQny6hdrBfgT76UjKsJKUYZWLq+WNqooZ+ifdxfoN+FzUTyZmMP6LHQQadXw1uQtRJiPf7jpLj7d+Z3U+h9qyQAq7pPJTmCNkwmGxr1ZreGCjaMIBYuY94HVodyfYMmDh7d5tBiKawN2/wk2vQtwqmNEFjq706dD7tIjm+we7YzLqUFAI8tNxMS2bDIvDfUxFrHPVWCzpojr5l8cgpgdcc5vPL5GWbWPc59txqioLJndBq1G4+8udPPXdX6Rn2wtYj5cFcvFUUjXI3QAbxCw+tyOki8RjYEkTHtkOm5iJH/0NgurCPb9BWCPv548/LKxYB74B9a+6NlVs4tOzuW/+bppGBbJs3z9YHXn/72IjA1nyQDc5ay9LLuyDb++GlJPQ+1lhGFdKZ0ZvthMr9l8g0+rghiYRHLqYzqvLDpJusRPir+elW1ozsn29PK6hxaHaNrOWSApFVWFWD8hOEzN4/1ARqvnqNlGdGh4rfNwLWxzL7SS55mUIi4EOE32Wx5xtc7B45xle/vkgRp0Gi92ZY0cgTMWeHtiCqX2b+uRaEi+c3ChsoEfOErP1UjJ/6yle/OmAu/gsMd3ClPm7OJOSxdMDW7DnTCprDonZef9WtXhj5DVEl9K8UGbFSGoeigK3fChSHn95TAi13h/GLBKe7ckn4KtbISu18PeDmOmf3w3LHoUvh/LjmvV5OnOV1PvFT69lUvfGDGgdjb9BS1SQEatDdGdSFPhi00n+Pn+51NeR5CLtAuz9SnzeuCc8vNsnog5igbxZtIm4eDM3vfsHQ6Zv5ExKFrWDjXy28QRrDl0iyKjjvdHXMntCx1KLenGQwi6pXtTvJJofHPjR8w/tFyL8YyKawsX9sOhO4d9eGFo9TPgZbpmO9fw+Bm+8jZXTHyIxJaXU6Ynzt55i1cF4DFoNyRlWQvz1qIh7UFKGlfvm75RpkL7iwFL4pBv8+owI3QHofBfqijQZWTSlK6H+elKybDhU0GsVLqZZSM200aNpJCsf78WtHeqXOPRSUmQoRlL9cDph/jDhof3gJhGCAUg9A3MGQdp5aNof7lwEuisXjiRdOsOezx/hBusWRmk/4KIS5U5PLEkmS+4GHsF+OtKz7eT/DwwPEGKfkmmTHjMlISMJfn1KNJqu2wFunQ2Rvglxzd96CrPFzuhODYg0GVl76BL3ztuV53do0Gp45uaW3N09htQsGxvjEthxMpldp1L45v6uhAaUvFhJxtglNZvL52DbJ2L2bsjljZ5wVNgRZCZBm5Fw2xeir+oVSDRbGPveUo5kmgCV1/wXM3TMg4Q2v6FEQ0s0Wxj4/gZ3RyatAq3rhhAXn062zdPEI9hPx+9P9pGiXhzsFvioE6RfgN5PQ4/HfVZh7Iqpg1js/t9tbbnz020FbswNw/15bURbPl53jF2nknGqYDLq6NgojP8Oa0NMZGCJxyBj7JKaTUh9sfhpzRAVpiAex0/+IcIyxmARrvnlcREHuQqJiqggrUUKg9SNhH49GL4Z7xPfmWB/PR/d1Z5Z4zrm2Z6WbWfr8SRxfRlvvzKZyeL3qDNCvxdhynoh7D60jRjctg6NI8Uk4URiBnfkE/WuseHUC/XjTHIWu08ncznTxrS+Tflp2g3se2kA8+7pXCpRLw5yxi6pnriMw8Iai+bDA9+ANS+J3PfB74jGCQtGgj1b+Mzc9KrX7BdvvU+zMtJ4NngV4/gFjTVDdGoa9H9FKkUvrJdqbGQgTtXJqaSC7pNPD2zBj3vPExdv5pXhbaSnTG6cTtg9V2QxDfsI2owo08vNWBfH9LVxWOx5dbN+mD/nUrIY2KYWNzSNLLPfUVFn7NJSQFI9aT1C9ElNOCyEffEY0TIvqqUnH370ArF9y0digbXXUwVOk7v36aIpXQHEomb8MIw3P8gdtp/EzN+Q0/M0+7I4VyEUdj5XI45awUaSM6zYcuW5v7XyCIC0HcjPuV3i5v3PXmjcG2q39enp8+eov7vqCB/9fszrsWlZNt66vR3Drq2Ln/7Kob3yQM7YJdUXcwLM7OpxgtQa4fEDniIngL9/gO8nC9G/+S3ocn+B0xRWhOKelbk61jvsovFHeCz0eAxienp9CijsfCAe98+nZDF53s48TTsC9Bo2PHOjjLe7WP0SbP4ATLVhwKvCn9+HmSf5c9T3n7vMvfN24ihELhuG+fPDtBvK/PcjF08lkvzCDjBsBnQYl/e4PfPh54fF5yM+gevuKtn1bNmwbYZYtM1IEI22b3hUeNJcZYE2P4cupDH0o004cnl3D21Xh2PxZmaO7UBslHhCqFGWv1kpoPMHvZ/IeLmwTzxlGYN8fqncIbPwAD2pWTa82agr4I6zl0eYTAq7pGaT38NdVSErSYRMHvkz76wdROPrlf8RYZvb55YuVmvLgn2LRIgn+QSM+QZaFL3/am5RCfXXY3M4ybB6vGU0Csye0IlrG4S6jxvStjYzxna8wlmrMFkpsP0zcdPs9ZRYEykD8qcybjmeyD1ztpPtyHucRsEt8mEBeu7q0pBawX7lcnOVMXZJzSa3cdjEX8S2uYMg6Zh3j5lu08BihvVvCAtgQyA0u6lk19b7i8bbHSYKczHXebZ8LPqzdp5yxYVWb3H4nm/+TlZOKqRThcnzdrn7qRq0Gpbvv0iXraeq18w9/RJsmwk7vwBrOjS/GWL7+uz00xbu5tH+zWleKyhPKuN3u84xrmsjXl1+sEDClMmoxWxx5Fn4XnXgkvv3VFmQM3ZJ9aUw47DWI8RjfLP+eY9XVVj1PGz9WHRZGvudKEP3FT9Ng70LQWuAa+8U7fqimns9NH8cfumf53n6u7+wO5xeQwLlWch01TUHX7FgJBxfJ+oNej4Bta/x2amnLdzN8v0XMWg1/PJID1Iyrdzx6bZCj88dcnEZtoFn4bu8spVkKEYiKYzv74XDy+H+jQUrElVV+Mzs/hL0gTDhJ2hwve+unRgnwj77FolUyxtfEP7xReB4gplJc3dwNjlvSmSov47fHutF7RB/342zEPIvKoKPxM2SDvu/EzYQdywQ/ugX94M+QNgr+5ijl9IZOn2T26sHKFBo5GJCt0aM69qIaQv3EBdv5ulBLZjaR/zdlPcahxR2iaQw0i6IRdWIJsLtMX8Ri9MpeqP+9Q0YQ2DSMqhzrW/HkJEo0jGb9BM3jqTjcHa78AcvxM8k0Wxh9KytnEjMKLCvdrAfY7o0IMRfz6TujX071nxj8JaHX6InBqcDTm2Efd/AwZ+Ed35UKxg+A+qX/XrB2ysPM3Pd8UIFHSDEX8/aJ3oTaTJWioVqKewSyZX4+wf47m7o82/o82zB/Q47fDcJDi0D/3CY9Isoaior1v8frP8fBEZDx0miEXe+rj6u2bJeo2C7Qqf7h29syhMDWpTZUPNbIkQEGlj5eK+iibqqinUG/zARQ3+vpVjQbj1cfN/1Ovo8bdFb2Gj7iSSW7y9av9vK5JUvLQUkkitxza3QdjT88ZYwC8uPVge3zYFmAyArGeYN84l9QKH0fgbGLxUpkhvehvevgaVT8xwyoVsMQ9rWxuZUiY0MdJe350YBPll/nHdXHcHucBbYX2FcOghrX4Xp7eHbSWJbUC3hovnkURj+sXDmLINc9DGfbWPm+mMcvZTudsy8kqiH+OtY+VjPPPYB5d2MurTIGbuk5pKVIkrRb3wBAiO9H2PLhsV3wfG1EBgFk5ZDVNnNhgFIPgm7vhCz24Gvi4+7v4TmAyG4rnsWCtD/3T9IzbLleXuAQUuov56fHrqBqCDhAe6rMEKxQzF7F4qm4a4K4Ma9oN0dJa8VKCLzt56ia2yEOy4O4p6hqmAyaDFbHYW+t3+raD6feD2JZgtLdp3FZNRVmmwjme4okVwN/zDRmAPydk/Kjd4P7lwIi8bAiXUw7xYxy4xuWXbjCm8MA17zfJ0YJxZ0UaBxTya0HQ3aW0h0+KPV5B2zv15DptVBptXBje/+wUdj2vPZhhOcS8nkTM6ia2lE6kqWCCv++ocJTbPh0C+igtcvWDzt+IcLf57Ww4vkp1NaXDP1KJOBj8d2YMr8XVzOsrtTFwsTdZNRS9fYCD6fKBbLI01G9yJpVUPO2CWStAsiPND33xDbx/sx1kzRoOPkHxAQIbJlfOxNckWSjsNfS2D/Ekg+garR8aTfK3yfHENUgBanoiUpw4pGgUHX1ObX/RcLLAoGGXUsmNyZ6xqGlSplMc977VYuH93AuR0/0yZtEyQfFweN/Y75ic0YfE1tIn381FDoWHKusWTXWeZuOkmC2ZonTfFKuI6r7P73cvFUIikq1kz4tKcIu0zdUriJly1LNMc+tgb8QmH8j6JpdnmiqnB+D/vXfsWdh26gbnQUS9ttxW//1/yW2YI1GU3o07c/DZpdy+2f7siT867XKqgqdGwUxvaTySVLWbRmCtMtYxDUaSf87WdcDxq9yPlvORRaDGb+AUvZpEXmwjUzz51X7soaur9XYz7bcPKqoq5VFIL9daRk2tz9Zyuzg6YUdomkOJzbDV/cJOK/Iz8p/Di7RXS6P7JceLqPWQwxJWu4UVrcs9Wzq+HPr3Ge2ojGkgaAwz+CZqkf4lQ1DNLsIEjJJFENwWIII7JWPc5maNmbqCHET4dOqyEpw0rTCD/ubB/OvZ1rAQoE5zhJbnxPFHTFHxRPDqoDrr1L/JxUVdzoGnYDo8k9Np+mRRbyvV9Ky2bGOvGEEOKvw2p3uqtzXc3Cr4RWo/D1fV1oEmWqMtYMUtglkuLy++uw4S244ytodUvhxzls8MN9wq5X5wej54uFzYrG6YDEONJP7mTB73t563I//HQavtK8SCfN0TyHHqQxg7NfB+AXw3+4RnMq77li+4hwE8BHHYX7ZXRrYdFQvxM06AIB4VccTqnSInORP9wyc/0x3vrtCLGRgdidTvfagQujTingl+5i1rgOnEjMYO7mUySkW9yz88qQo14U5OKpRFJcej0FcSuFp0vLoYWn3mn1oqWeMRj2zBNZMyM/hba3l+9486PRQnRLfjzux1uXw92z49Px7RkwbzUB1mT6NYC05HjOZHr+9b929CPamYKqKihGE21j69Lp2mtxBaQSJ25gxYGEUoueM2cS6cpY2XYi6aqi6gq3LNh6mkVTurJk11kW7zgDiDREg65gxnZ+UQ80aMm0OVBVeOGnA/z6aE9Gd2qQ55qRJmOlF/XiIGfsEkluUs+KxVFDwRzxAqiqSJfc/AGgwM1vevVzrwgKW1Sc2qcpF9Oy6PXmOqxezMV1GgV7Tp78G7e2pU6IH3fP3cmJxAz37Da/C2LutEDAnYrpCm+4YtcAESYDSWare9vTg1rw457Cu0PlDunoNOCKroT467icZS/0+zcZdZgtdppGBbL4/m4kZ1gZ+/n2PLP0qogMxUgkpcGaCZf+hgadr37spg9E2z2AHv8SPTd9WGjja+ZuPsl/lx28asZI/v23dajH5Swbaw7FA56GzlMX7CYp05NLHxVk5I5ODfh43TG3gAcYtGTmpBm6zqtVFAw6hSybk9jIQP7vtraEBxrYcjyJHk0jmbP5JGeSszgeb+Z8asGWgfkJMGjItArlf2pgc4L89OVjVlaOlIuwK4oyCngZaAV0VlW1SGothV1S6Vk6TfiXPLgZwhpd/fg/v4afHhILi+3HwdAPRfVqJWXm+mMs3nGGM8lZaDWKu6GHQafBepVFR50C3kLYwX7i+03LtqNB3NscqkfIc/uYF0a0yUC82cqI6+qy6uAlAg06VNQ83aQKI3fMvSrPyq9EeVkK/A3cCmwo5XkkkspF76fFx6VThSnY1bjuLhizSHT42fuV6KVqLWjWVVkwGXWcSc6iWbSJ7f/px6L7uhAbGYjV7uT6mLArvreQdUnSs+3UCfVHAZzgbiNXVFEHiM8R8KV//kOm1UGC2eIWdW/PQAatyGxpFm3iRGIGd3ZuWG1FvTj4JBSjKMp64Ek5Y5dUK/Z+JTzUB7wO3R8q2nvO7oSvR4uKy7od4K4lBbs1VRIKi8Mv2n6GsykFQx+9m0fy17nLpGTaCuzzJf56DcH+evRaDamZNswWOwEGLXMmdeKu2dsL3CCaRZuYMbaDezG2OiNNwCSS0nLdWGgxBNa+ApcOFO09Da6HyasgtCH8swfmDBDt8SohE7rF5Ek9dH3uEvWwAD2h/p5wUpfYCGaNLzzHO9RfT3igvtD9tYKNDGhd64pjCg/Q0ygikEtpFix2J23rhdCnRRSvj7iGh7/ei1MVs/8AvZAug1ZDXLy5Roh6cbhqEFBRlDVAbS+7nlNV9aeiXkhRlCnAFICGDRsWeYASSYWhKMJLZtkjot1dUYlsBpNXw8LbRbOIz2+Csd+Wf5VqCXBltjSODODbB7oDMGrWFk4mZgLw0Nd7C32vy4zMlVnjr9eQbXOiAkatwoLJXQgPNLD+SLw7I8cVfzdqFUx+epIyrAT761n/VB9iIgLd556/9RQJZitRJgML7+tKeKAhT1GRFPW8yFCMRFJWWNLhm/HCPEwfCKPnlbyPajlSmI8MwIs/HSDCZCDDYifb5ll76NUskv6ta7l9zl059MkZVsbO3kaC2ZontTEi0EC2zeHu2frLIz3yiLW3OHm5teSrxJRruqMUdkm1JysFfn4Yrr8PYnsX/X12q3jfX4tB0cKw6SJrpooyc/0xvt11lpOJmUQEGnCqKimZtjxWAVe7MeT2j7n5gw0kmK1VrgK0oiiXylNFUUYCHwFRwHJFUf5UVbUS1FZLJD5GaxCNNn58QKRAXqWc3o3OACNniW5Im94Ti7GXz4usm0qc614YJqOOk4mZ3m17cwQ5vyjnr+rMLfq/PtarWleAVhSyQEkiKSr/7BXx8uYDhZ9McYV5x2z49Wnhu9JhIgx5r1LnuheGDIlUHDIrRiLxNXXbi6rSw7+IjkbFpfN9MHqBMA5zecxU4lz3wvCWTSNFvXIhhV0iKQ7dHoImN4qwiv3q1ZAFaDVUdGDyDxeGY18OAXO878cpqdFIYZdIioNGI5wc710r4ucloWGXnFz3RiK888VNkHjMt+OU1GiksEskxcUULV5OBxxdWbJzRDaDe9dAnesg5RR80R/ObPPpMCU1FynsEklJ2TVH2AccLHKdXl5M0TBpOTQfJNIp5w2Dv3/w7RglNRIp7BJJSekwEep1hJ8eFrPukmA0wR0LodNkcFjgu7th0/vC610iKSFS2CWSkqIziE5KqPDd5JItpoJIeRzyLtz0ivh6zcuiqMlRtmZbknLm/G6YP1ysq5QxUtglktIQ3hiGfQTnd8Hvr5T8PIoCNzwKo+aJdMi9C+Cr2yAr1XdjlVQMicdgyQSYfSOcWA8b3inzS1a96giJpLLRZgQkPg9Nb/TNuUIawKI74eQf8Hl/uOsbiGhS+nNLyhdzPKz/P1HzoDrEDbvL/XDDY2V+aVl5KpH4Gls26P1Kd47UM/D1nRB/APxCYdSX0KSvT4YnKWOsGbB1Bmz+EKxmUDTCH6j3sxBSr1SnlpWnEklFsPYVUXRU0ni7i9CGMHml8IPPThVhme2fykXVyozTKVokftQR1r0uRL35zfDgVhGuK6WoFwcp7BKJL6lzrYi3u5pblwZjkPCk6fmEeJT/9WlhQmbNLP25Jb7l1GaY3QeWPgjpF0R9wqTlcNdiiG5Z7sORwi6R+JLWw6Hz/bBtJhxaVvrzaTTCn+b2ucLT/a/FoitTSdMrJb4l5bRYGP1yMFzYB0F1RWXyfesgpkeFDUsKu0Tiawa8KvqdLp0GySd9c85rbhWVquGxoivTp718c+OQlAyLGda+Ch9fLwrUdP7Q59/w8G649k5xQ65ApLBLJL5GZxSLnYZASDjsu/PWai1mgi2GQPZl+GYcrHga7BbfXUNyZXLH0Te+I4rK2o6Ch3dBn2fBEFDRIwRkVoxEUnb4IjvGG6oK22fBqhfAaYPabeHW2RDdyvfXkng4tRlWPecpMKrbAW5+Exp0LrchyKwYiaSi0fsJEd77FRz40XfnVRTo+qBwiAyLyQnN9IYtH4kZpcS3JB6DxWNFHP2fvRBUx+PwWY6iXhxkgZJEUpaoTtgzHy4dgFptIbKp785drwM8sAlWPicad6x6Hg6vgFs+gKgWvrtOTSX9EvzxpqfASB8giou6PyTCbJUYOWOXSMoSjRZunyN6pi6ZALYs357fGCQaZI/5BgKj4cwW+OQG+P11EQqSFJ+sVPj9NZh+HezK8QJqPx4e2Qt9nqn0og5S2CWSsiekPtz6magi/fXpsrlGi0EwbbtwnHTaYMNbMLOryJyRRU1Fw5IOf7wNH7aDDW+DLRNaDoWp22D4xxBUu6JHWGTk4qlEUl6s+a9oqXffOhFGKStOb4VfHoeEQ+LrRj1g4GuiZ6ukIFmpotH4thnCFx+gcS/o+7zodlWJKOriqRR2iaS8cNjh9GaI7V0O17KJ2PC6NyArWWxrNUzkWtdqXfbXrwqkX4Idn4mXJU1sa9AVbnxOCHslRAq7RFKZOb9HLHCWdbw2K1XkW++YDfZsQBEOkj0eF/YHNZFLB8Xs/K8l4Mjx9InpCb2fFh8VpWLHdwWksEsklZW0C/DhtdBmJIycVT5CknZBhIF2f+kRs9g+0P0RaHJjpRYzn2C3iPWGXXPEUxMACrQcIn4GlSzkUhhS2CWSysz6N2H9G3DLdOg4sfyue/kcbJ0pBN6WIbZFNIPrJ8O1Y8A/tPzGUh5c+Av2LRKz88xEsc1gEt9r1wernM+9FHaJpDLjdAgr3tNbhAdMnQRrSc4AAA4pSURBVHble/2sFNj5Bez8XLgRgvA7aTkE2t0hvN+1+vIdk69IOg4Hl4rG4Jf+9myPbgPX3yO+P2NQxY2vFEhhl0gqOxmJMKun8Ja5fwP4BZf/GBw2OPKryNc+sd6zPSACmg+C5gMhtm/FjK2oOJ1w4U84uhKOLBeVuC78QoWXy3VjhAVAFQ85SWGXSKoCp7eKFni9nhLFTBVJymnY/y389Q0kHvVs1+iEv3iDLiIWXbsthMZUnIOh0wlJx+D0JuHfcmojmC959huCoOVgaD0CmvYTN85qghR2iaSqUVamYcVFVSH+oJgBx62Cs9uFNUJudP4Q1RxCG4kCrOC6ovLVPwwCwsEYLJwODYGiFF9rKN5s2WET4SLzJXHDST0jxPzS3yKrxZqe9/jg+uLpovlAaNy7cvwcywAp7BJJVeLCPtHjdNRcaNi1okeTl+zLcG6XEPhzOyH+kCcuX1QUjWjm7H4ZQKMX2xWNuHE4LKKloDUDLJevfD5TbfFziukhXlEtq3yYpSgUVdilCZhEUhkIixEhg28nwf0bwRRV0SPy4BciQhpN+3m2ZaVAYpyYSaedh7R/xJpBVoooiLKYhUBbzaI032kXH21FbeuniNl/YJTo/xrWSPyMarURZmqV6edTCZHCLpFUBvxC4I4F8Hl/+P4eGPcjaCvxv6d/mLCsLaptrcMO9iyRT263iGIpp13M1FUnoIgbm84oQjd+IRW/5lCFqcR/ORJJDaN2WxjyHvw0Fda9Bv1frugR+Q6tDrRBVTbNsKohhV0iqUy0HwvndoiFQqdDzlolJUIKu0RS2Rj8jkgxrAGLgZKyQfqxSySVDa1eiHrKafj+XrEQKZEUAynsEkllJfkE/P09/PyQbJYhKRalEnZFUd5WFOWwoih/KYryo6Io1cxBSCKpQJr0hX4vikbYWz6q6NFIqhClnbGvBq5RVbUdcBT4d+mHJJFI3NzwGLQeDmteguPrKno0kipCqYRdVdVVqqrac77cBtQv/ZAkEokbRYHhMyGyBWx8V4ZkJEXCl1kx9wDfFLZTUZQpwBSAhg0b+vCyEkk1x2iCcd8Jp0KZKSMpAledsSuKskZRlL+9vIbnOuY5wA4sLOw8qqp+pqpqJ1VVO0VFyXJgiaRYhNQXAm8xw+YPhcOhRFIIV52xq6ra/0r7FUWZBAwF+qkV4SgmkdQkDi+H1S8KT5b+L1f0aCSVlFKFYhRFGQQ8DfRWVbWo7j4SiaSktBsNZ7bApvdFS7v2Yyt6RJJKSGlj7B8DRmC1ImJ/21RVfaDUo5JIJN5RFFGZmnIKlj0qnA8b96zoUUkqGaXNimmqqmoDVVWvy3lJUZdIyhqtHkbNg/DG8MvjwlNGIsmF9IqRSKoi/qFw1xIxg5dGYZJ8SEsBiaSqEt5YNJ9QVdj+mfSUkbiRwi6RVHUu7IPfnoUlE0SvUEmNRwq7RFLVqXsdDH0fjq+Fn6bJHHeJjLFLJNWCjhMhIx5+fw0MJhjyrqxSrcHIGbtEUl3o+STc8CjsmQ+XDlT0aCQViJyxSyTVBUWB/v+FdndCrdYVPRpJBSJn7BJJdUJRPKK+/ztY/6Z0hKyByBm7RFJdOfmHCMs4bdD3ORlzr0FIYZdIqitDPxSz9Q1vg9MO/V6S4l5DkMIukVRXNBq4ZTpodMI0LDtNZsvUEKSwSyTVGY1G5Lj7BQuBl6JeI5DCLpFUdxQFbnrFs4gafwhCGojGHZJqicyKkUhqCooC1kyYPwK+HALplyp6RJIyQgq7RFKTMATALR9C4lH4vD8kHKnoEUnKACnsEklNo8UgmLQc7NnwxU1wfF1Fj0jiY6SwSyQ1kXod4N41EFwPdn1R0aOR+Bi5eCqR1FTCGsHkVZ6vzQkie0ZnrLgxSXyCFHaJpCZjDBIfHXb46laREjl6nuilKqmyyFCMRCIBrQ56PQVJx+DTXhC3uqJHJCkFUtglEomg9TCYsl7E3RfeDiufA7ulokclKQFS2CUSiYeIJmJR9fp74fjv0hmyiiJj7BKJJC96f+EpY80EvZ/wmNk9F7o8IBdWqwhyxi6RSLxjCBAfD/0Mq1+EWT3h9NaKHZOkSEhhl0gkV6b9OBj7HdiyYO4g+GEKXD5X0aOSXAEp7BKJ5Oo0uwmmbYMe/4IDS+G3Zyt6RJIrIGPsEomkaBgCof9L0HESkLOomhgHh5ZB5ynSLbISIWfsEomkeIQ1grAY8fmRX+H/27v/2KruMo7j74fSDloB5woZG4UKwhwwDE2ZbBjHBjEEF1hkURYZLpItzP0y8x8Nf0wlRkxgiRoSJYthOjdwm8HGzajARg3aQpGfw40xYJQf0s7xYwJjhT7+8T21zQK9l95zz+m9/bySE87p+ZbzPD23T8/9nu/9ng0/gJ9OgvrlcP5kqqFJoMIuIj037XFYtB5umAwbl8LTE2DD0rSj6vNU2EUkN1VTYMHLsHgz3Hx3mDUSwhj45q0aCw9hyoYEqY9dROJx/UT4yqrOQv7u32H1bBg2HqYsgklf65ybpthdvADNW+DgJjhYD2eOwbd3J/ZoQl2xi0i8OorXjTUw5+dhYrFXvgMrPgt/fBLOvZ9ufPng3vkHrXEVLBsFz94Nf1sB7ZfglnsTnZ5BV+wikh+lA6FmIUy+H45ug63PwL4/w6wfh/3/3h1uwhbqVfyFD+DAJnj7L7B/A8x/LtxrGHpTyHv0dKieBgOGJB6aCruI5JcZjKgNy6U2KCmF9nZY83U42wrj58Ln7oPqL0C/krSjzezkIVj3CDQ3QPtFKBsEY6YD0TuV0XeEJUUq7CKSnJLS8K8ZzHsGtj8He34PO1+AiqEw4ymouT/dGDu0X4ITb8DhBjhUDyOmwLQnoGIYXDwPtz8Gn5kJVZ/vzKuXUGEXkeSZQdWtYZm1LHRn7F0HFZVhf8ubYXz8mLtCm2ETwpzx+eIOH56CgdeG7Ze+Geakv3AmbA+pghtqwnpZOTy4MX+xxECFXUTSVVYOE+4JS4cPjsGJPfDWq2G7tCL0X9+zMvTLnzke+rg/MSz0YWcabeIObefCp2ch9I0fbQrdKq37oOVfMHg4PNIY9g/8FEycB6Nuh5FTC+6JUjkVdjNbCswF2oEW4AF3PxZHYCLSh425C57YBacOw5Gt0NwIR/8J1wwO+7ethk3LwnpJWSjYpRXw6JawXr8cdq0NE5d9dDYs/UpgyfHwPTueh11roLwSKseFUSvDJ3Ue/8vLE003buY5fHjAzAa7+5lo/XFgvLsvzvR9tbW13tTU1OPjikgf17oPju+E/56Asy2hcLedj4ZXloS++7f/CqXlodCXlUP5dXDbo2H/2feg/4CCm9/GzLa5e22mdjldsXcU9UgF/58ZSEQkj4aOC8uVTF4Qlivp6MsvUjn3sZvZj4CFwGngzm7aPQQ8BDByZGH1V4mIFJKMXTFmth64/jK7lrj7H7q0+x4wwN2fynRQdcWIiFy92Lpi3H1mlsf8LfAqkLGwi4hI/uQ0V4yZje2yORd4M7dwREQkV7n2sS8zs5sIwx3fBTKOiBERkfzKdVTMvLgCERGReGjaXhGRIqPCLiJSZHL65GmPD2rWSuiT74lK4L0YwykEyrlvUM59Qy45j3L3oZkapVLYc2FmTdmM4ywmyrlvUM59QxI5qytGRKTIqLCLiBSZQizsq9IOIAXKuW9Qzn1D3nMuuD52ERHpXiFesYuISDd6bWE3s1lm9paZ7Tez715m/zVmtjba32hm1clHGa8scn7SzPaa2S4z22Bmo9KIM06Zcu7Sbp6ZuZkV9AiKbPI1s69G5/kNM3s+6RjjlsXreqSZvWZm26PX9uw04oyTmf3KzFrMbM8V9puZ/Sz6mewys5pYA3D3XrcAJcA7wGigDNhJeDpT1zbfAn4Rrc8H1qYddwI53wmUR+sP94Wco3aDgHqgAahNO+48n+OxwHbg2mh7WNpxJ5DzKuDhaH08cCjtuGPI+4tADbDnCvtnA38CDJgKNMZ5/N56xX4rsN/dD7j7R8AawuyRXc0Fno3WXwJmmGV6om2vljFnd3/N3c9Fmw3AiIRjjFs25xlgKfAT4MMkg8uDbPJ9EFjp7icB3L0l4Rjjlk3ODkQPM2UIUPDPTXb3euD9bprMBX7tQQPwSTMbHtfxe2thvxFo7rJ9JPraZdu4+0XCE5yuSyS6/Mgm564WEf7iF7KMOUdvUavc/ZUkA8uTbM7xOGCcmW02swYzm5VYdPmRTc7fBxaY2RHCMx0eSya0VF3t7/tVyfnReJI8M1sA1AJ3pB1LPplZP+Bp4IGUQ0lSf0J3zHTCO7J6M7vF3U+lGlV+3QesdvcVZnYb8Bszm+ju7WkHVqh66xX7UaCqy/aI6GuXbWNm/Qlv4f6TSHT5kU3OmNlMYAkwx90vJBRbvmTKeRAwEXjdzA4R+iLrCvgGajbn+AhQ5+5t7n4Q2Eco9IUqm5wXAb8DcPd/AAMI86kUs6x+33uqtxb2rcBYM/u0mZURbo7WfaxNHfCNaP1eYKNHdyUKVMaczWwy8EtCUS/0vlfIkLO7n3b3Snevdvdqwn2FOe5eqA/MzeZ1vY5wtY6ZVRK6Zg4kGWTMssn5MDADwMxuJhT21kSjTF4dsDAaHTMVOO3ux2P739O+e9zNXeXZhKuVdwgPzgb4IeEXG8LJfxHYD2wBRqcdcwI5rwdOADuipS7tmPOd88favk4Bj4rJ8hwboftpL7AbmJ92zAnkPB7YTBgxswP4Utoxx5DzC8BxoI3wLmwR4Qlzi7uc55XRz2R33K9rffJURKTI9NauGBER6SEVdhGRIqPCLiJSZFTYRUSKjAq7iEiRUWEXESkyKuwiIkVGhV1EpMj8D+p/ufo7/7PLAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plot(m)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The plot is much more satisfying. The model now recognises the correlation between the two functions and is able to suggest (with uncertainty) that as x > 0.5 the orange curve should follow the blue curve (which we know to be the truth from the data generating procedure above). "
]
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.4"
}
},
"nbformat": 4,
"nbformat_minor": 2
}

Computing file changes ...