{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Exercise 9: Multivariate Classification\n", "\n", "This exercise is dedicated to studies of a multivariate data-set, the famous 'Iris flower data-set' used in the year 1936 by \n", "Sir Ronald Fisher as an example of discriminant analysis. The data-set describes the morphological variation of Iris flowers and contains 50 samples from each of the three species Iris-setosa, Iris-virginica and Iris-versicolor.\n", "\n", "The file `iris.data` (available in the exercise repository and cloned into this directory) contains the data. The columns contain four morphological features of the three species: the sepal length and width and the petal length and width (all in $\\mathrm{cm}$), and the species as classified by a botanist (Iris-setosa = 0, Iris-versicolor= 1, Iris-virginica = 2). Fortunately, the biological details are not of too much importance for the studies proposed here, and you may safely name the variables $v_1, \\ldots, v_4$ and call the classes '1', '2' and '3'.\n", "In this exercise you will develop classifiers based on this data-set as a training sample to identify different iris species, for people who are not particularly gifted botanists.\n", "\n", "To guide your own implementation of studies of the data-set, a code basis is provided in this jupyter notebook. Inspect it and use it as a starting point for your own work. " ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "\"\"\"\n", ".. module:: ex9_template\n", " :synopsis example for discriminant analysis\n", "\n", ".. moduleauthor Thomas Keck\n", "\"\"\"\n", "# ------------------------------------------------------------------------\n", "# useful imports\n", "\n", "import bisect\n", "import numpy as np\n", "\n", "import matplotlib\n", "# matplotlib.rcParams['backend'] = 'TkAgg'\n", "\n", "from matplotlib import pyplot as plt\n", "\n", "# ------- load the Data set --------------------------------------------\n", "# Load iris data\n", "data = np.loadtxt('iris.data')\n", "\n", "# Define dictionary with columns names: s=sepal, p=petal, l=lenght, w=width,\n", "columns = {0: 'L(Kelch)', 1: 'W(Kelch)',\n", " 2: 'L(Blatt)', 3: 'W(Blatt)', 4: 'class'}\n", "\n", "# Define boolean arrays corresponding to the three\n", "# classes setosa, versicolor and virginica\n", "setosa = data[:, 4] == 0\n", "versicolor = data[:, 4] == 1\n", "virginica = data[:, 4] == 2\n", "\n", "# Signal is versicolor (can be changed to setosa or virginica)\n", "signal = versicolor\n", "bckgrd = ~signal\n", "# !! note: see indexing of arrays with boolean array in python documentation\n", "#\n", "# exmaples how to access the data:\n", "# data[signal] # All events classified as signal\n", "# data[background] # All events classified as background\n", "# data[setosa] # All events classified as setosa\n", "# data[versicolor | virginica] # All events classified as versicolor or virginica\n", "# data[:, :2] # The first two columns (sepal length and sepal width) of all events\n", "# data[signal, :2] # The first two columns of the signal events\n", "# data[background, 2:4] # The 3. and 4. column (petal length and petal width) of the background events\n", "# data[:, 4] # The label column of all events\n", "# (see the numpy documentation for further examples)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 9.1: (voluntary)\n", "\n", "In multivariate analyses, it is most important to gain an overview and some initial level of understanding before deciding\n", "on strategies for optimal classification. One-dimensional distributions (=histograms) and pair-wise scatter plots of the \n", "available features (=variables) are the best way to visualise the available raw information. Plot the distributions for all features, ideally arranged as two-dimensional array of histograms (on the diagonal) and scatter plots (off-diagonal). You can use the examples in the template class `Plotter` defined in [Exercise 9.2](#exercise92) to make these plots. With appropriate colour coding, all three classes of the data-set can be shown in the same set of histograms. \n", " \n", "Which variables provide the best separation of Iris-setosa from the other two species? How can Iris-versicolor or Iris-virginica be separated from the others?" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA5IAAAKrCAYAAACDc7LDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nOzdbYjl93ke/uu2Nkqo69jF2oDRbiKFruts3YKVQXUJNCp2y0oF7Yv0H7RgUgfhJakVCgkFFRfHKK/c0BQMatwtNbIDsaz4RVnIGkFSGYGJHI2wo1gyChvFqVYx1fqhemNsWfT+vzjH8Wi8o/n9ds+Z85Xn84GF8/D1nIuze2GuOQ+q7g4AAABM9bpNBwAAAOC1xZAEAABgFkMSAACAWQxJAAAAZjEkAQAAmMWQBAAAYJZ9h2RVfayqXqiqL+1xf1XVR6rqYlU9WVW3rD4mkOgjjEQfYQy6CJsx5RXJB5KcepX7b09yYvnnbJLfvfZYwB4eiD7CKB6IPsIIHoguwoHbd0h296NJvvEqR04n+UQvPJbkTVX1llUFBL5PH2Ec+ghj0EXYjCMr+Bk3Jnlux/VLy9u+uvtgVZ3N4jdBef3rX/+zb3vb21bw8DCuJ5544mvdffQAH1IfYQ+j9lEXOYwOuI/+vxH2cC1dXMWQnKy7zyU5lyRbW1u9vb19kA8PB66q/nrTGfaijxw2o/ZRFzmM9BHGcC1dXMW3tj6f5PiO68eWtwEHTx9hHPoIY9BFWINVDMnzSX5p+Y1Y70zyYnf/wFsFgAOhjzAOfYQx6CKswb5vba2qTya5LckNVXUpyW8m+ZEk6e6PJrmQ5I4kF5N8K8kvryssHHb6COPQRxiDLsJm7Dsku/vMPvd3kvevLBGwJ32EcegjjEEXYTNW8dZWAAAADhFDEgAAgFkMSQAAAGYxJAEAAJjFkAQAAGAWQxIAAIBZDEkAAABmMSQBAACYxZAEAABgFkMSAACAWQxJAAAAZjEkAQAAmMWQBAAAYBZDEgAAgFkMSQAAAGYxJAEAAJjFkAQAAGAWQxIAAIBZDEkAAABmmTQkq+pUVT1TVRer6t4r3P+TVfVIVX2hqp6sqjtWHxVI9BFGoYswDn2Eg7fvkKyq65Lcn+T2JCeTnKmqk7uO/cckD3X3O5LcleS/rjoooI8wCl2EcegjbMaUVyRvTXKxu5/t7peSPJjk9K4zneTHl5ffmORvVhcR2EEfYQy6COPQR9iAKUPyxiTP7bh+aXnbTh9K8p6qupTkQpJfu9IPqqqzVbVdVduXL1++irhw6OkjjEEXYRz6CBuwqi/bOZPkge4+luSOJL9XVT/ws7v7XHdvdffW0aNHV/TQwC76CGPQRRiHPsKKTRmSzyc5vuP6seVtO92d5KEk6e4/SfJjSW5YRUDgFfQRxqCLMA59hA2YMiQfT3Kiqm6uquuz+IDy+V1n/neSdyVJVf1MFuX0fgBYPX2EMegijEMfYQP2HZLd/XKSe5I8nOTLWXzj1VNVdV9V3bk89htJ3ldVf5bkk0ne2929rtBwWOkjjEEXYRz6CJtxZMqh7r6QxQeTd972wR2Xn07yc6uNBlyJPsIYdBHGoY9w8Fb1ZTsAAAAcEoYkAAAAsxiSAAAAzGJIAgAAMIshCQAAwCyGJAAAALMYkgAAAMxiSAIAADCLIQkAAMAshiQAAACzGJIAAADMYkgCAAAwiyEJAADALIYkAAAAsxiSAAAAzGJIAgAAMIshCQAAwCyGJAAAALMYkgAAAMwyaUhW1amqeqaqLlbVvXuc+cWqerqqnqqq319tTOB79BHGoIswDn2Eg3dkvwNVdV2S+5P8iySXkjxeVee7++kdZ04k+Q9Jfq67v1lVP7GuwHCY6SOMQRdhHPoImzHlFclbk1zs7me7+6UkDyY5vevM+5Lc393fTJLufmG1MYElfYQx6CKMQx9hA6YMyRuTPLfj+qXlbTu9Nclbq+pzVfVYVZ260g+qqrNVtV1V25cvX766xHC46SOMQRdhHPoIG7CqL9s5kuREktuSnEny36vqTbsPdfe57t7q7q2jR4+u6KGBXfQRxqCLMA59hBWbMiSfT3J8x/Vjy9t2upTkfHd/t7v/KslfZFFWYLX0EcagizAOfYQNmDIkH09yoqpurqrrk9yV5PyuM/8zi9/wpKpuyOLtA8+uMCewoI8wBl2EcegjbMC+Q7K7X05yT5KHk3w5yUPd/VRV3VdVdy6PPZzk61X1dJJHkvz77v76ukLDYaWPMAZdhHHoI2xGdfdGHnhra6u3t7c38thwUKrqie7e2nSO/egjh8FroY+6yGGhjzCGa+niqr5sBwAAgEPCkAQAAGAWQxIAAIBZDEkAAABmMSQBAACYxZAEAABgFkMSAACAWQxJAAAAZjEkAQAAmMWQBAAAYBZDEgAAgFkMSQAAAGYxJAEAAJjFkAQAAGAWQxIAAIBZDEkAAABmMSQBAACYxZAEAABglklDsqpOVdUzVXWxqu59lXO/UFVdVVuriwjspI8wBl2EcegjHLx9h2RVXZfk/iS3JzmZ5ExVnbzCuTck+XdJPr/qkMCCPsIYdBHGoY+wGVNekbw1ycXufra7X0ryYJLTVzj3W0k+nOTbK8wHvJI+whh0Ecahj7ABU4bkjUme23H90vK2v1VVtyQ53t1/+Go/qKrOVtV2VW1fvnx5dlhAH2EQugjj0EfYgGv+sp2qel2S30nyG/ud7e5z3b3V3VtHjx691ocGdtFHGIMuwjj0EdZjypB8PsnxHdePLW/7njckeXuSz1bVV5K8M8l5H2KGtdBHGIMuwjj0ETZgypB8PMmJqrq5qq5PcleS89+7s7tf7O4buvum7r4pyWNJ7uzu7bUkhsNNH2EMugjj0EfYgH2HZHe/nOSeJA8n+XKSh7r7qaq6r6ruXHdA4Pv0EcagizAOfYTNODLlUHdfSHJh120f3OPsbdceC9iLPsIYdBHGoY9w8K75y3YAAAA4XAxJAAAAZjEkAQAAmMWQBAAAYBZDEgAAgFkMSQAAAGYxJAEAAJjFkAQAAGAWQxIAAIBZDEkAAABmMSQBAACYxZAEAABgFkMSAACAWQxJAAAAZjEkAQAAmMWQBAAAYBZDEgAAgFkMSQAAAGYxJAEAAJhl0pCsqlNV9UxVXayqe69w/69X1dNV9WRV/XFV/dTqowKJPsIodBHGoY9w8PYdklV1XZL7k9ye5GSSM1V1ctexLyTZ6u5/nOTTSf7TqoMC+gij0EUYhz7CZkx5RfLWJBe7+9nufinJg0lO7zzQ3Y9097eWVx9Lcmy1MYElfYQx6CKMQx9hA6YMyRuTPLfj+qXlbXu5O8lnrnRHVZ2tqu2q2r58+fL0lMD36COMQRdhHPoIG7DSL9upqvck2Ury21e6v7vPdfdWd28dPXp0lQ8N7KKPMAZdhHHoI6zOkQlnnk9yfMf1Y8vbXqGq3p3kA0l+vru/s5p4wC76CGPQRRiHPsIGTHlF8vEkJ6rq5qq6PsldSc7vPFBV70jy35Lc2d0vrD4msKSPMAZdhHHoI2zAvkOyu19Ock+Sh5N8OclD3f1UVd1XVXcuj/12kr+b5A+q6otVdX6PHwdcA32EMegijEMfYTOmvLU13X0hyYVdt31wx+V3rzgXsAd9hDHoIoxDH+HgrfTLdgAAAPjhZ0gCAAAwiyEJAADALIYkAAAAsxiSAAAAzGJIAgAAMIshCQAAwCyGJAAAALMYkgAAAMxiSAIAADCLIQkAAMAshiQAAACzGJIAAADMYkgCAAAwiyEJAADALIYkAAAAsxiSAAAAzGJIAgAAMMukIVlVp6rqmaq6WFX3XuH+H62qTy3v/3xV3bTqoMCCPsIYdBHGoY9w8PYdklV1XZL7k9ye5GSSM1V1ctexu5N8s7v/fpL/kuTDqw4K6COMQhdhHPoImzHlFclbk1zs7me7+6UkDyY5vevM6SQfX17+dJJ3VVWtLiawpI8wBl2EcegjbMCUIXljkud2XL+0vO2KZ7r75SQvJnnzKgICr6CPMAZdhHHoI2zAkYN8sKo6m+Ts8up3qupLB/n4E9yQ5GubDnEFI+aSaZp/sOkAe9HHqyLTdCPmGrKPr4EuJmP+fco0zYiZEn28WqP+fY6YS6ZprrqLU4bk80mO77h+bHnblc5cqqojSd6Y5Ou7f1B3n0tyLkmqaru7t64m9LqMmCkZM5dM01TV9op/pD5ukEzTjZhrxX08NF1Mxswl0zQjZkr08WqNmCkZM5dM01xLF6e8tfXxJCeq6uaquj7JXUnO7zpzPsm/WV7+10n+V3f31YYC9qSPMAZdhHHoI2zAvq9IdvfLVXVPkoeTXJfkY939VFXdl2S7u88n+R9Jfq+qLib5RhYFBlZMH2EMugjj0EfYjEmfkezuC0ku7LrtgzsufzvJ/zfzsc/NPH8QRsyUjJlLpmlWnkkfN0qm6UbMtdJMh6iLyZi5ZJpmxEyJPl6tETMlY+aSaZqrzlRe1QcAAGCOKZ+RBAAAgL+19iFZVaeq6pmqulhV917h/h+tqk8t7/98Vd00QKZfr6qnq+rJqvrjqvqpTWface4Xqqqrau3f+DQlU1X94vK5eqqqfn/dmabkqqqfrKpHquoLy7/DO9ac52NV9cJeXxFeCx9Z5n2yqm5ZZ55Xo4+rybTj3KHu42hdXD7ma6KPuri6XDvO6aM+XhV9XE2mHecOrItTcx32Pq6ti929tj9ZfOD5L5P8dJLrk/xZkpO7zvzbJB9dXr4ryacGyPTPk/yd5eVfHSHT8twbkjya5LEkW5vOlOREki8k+XvL6z+xzkwzcp1L8qvLyyeTfGXNmf5ZkluSfGmP++9I8pkkleSdST6/7ufpGp47fdTHVWY60C4uH2f4PurianMtz+mjPq7zuTv0fRyxizOeq0Pfx3V1cd2vSN6a5GJ3P9vdLyV5MMnpXWdOJ/n48vKnk7yrqmqTmbr7ke7+1vLqY1n894jWacrzlCS/leTDSb695jxTM70vyf3d/c0k6e4XBsnVSX58efmNSf5mnYG6+9EsvgFuL6eTfKIXHkvypqp6yzoz7UEfV5Rp6bD3cbguJq+ZPuriCnMt6aM+Xi19XFGmpYPs4tRch76P6+riuofkjUme23H90vK2K57p7peTvJjkzRvOtNPdWSz0ddo30/Il5uPd/YdrzjI5U5K3JnlrVX2uqh6rqlOD5PpQkvdU1aUsvsHt1w4g16uZ+29ukzn0UR9XmelDGauLyRh91MXp9HF1mT4UfbzaDPo4Zhcn5Yo+TnFVXZz0n/84rKrqPUm2kvz8hnO8LsnvJHnvJnNcwZEs3i5wWxa/CXu0qv5Rd//fjaZKziR5oLv/c1X90yz+u1Fv7+7/t+FcXAN93NeIfdTFH0KjdHGZRR+n08cfQqP0ceAuJvq4Nut+RfL5JMd3XD+2vO2KZ6rqSBYv7359w5lSVe9O8oEkd3b3d9aYZ0qmNyR5e5LPVtVXsnjv8vk1f4h5yvN0Kcn57v5ud/9Vkr/IoqjrNCXX3UkeSpLu/pMkP5bkhjXnejWT/s0NkkMf9XGVmUbrYjJGH3Vxdbn0cXomfbz6DPo4Zhen5Er0cYqr6+J+H6K8lj9Z/Abg2SQ35/sfNv2Hu868P6/8APNDA2R6RxYfkj2xzixzMu06/9ms/8sEpjxPp5J8fHn5hixeEn/zALk+k+S9y8s/k8X7zmvNuW7K3h9g/ld55QeY//Qg/l1d5XOnj/q4ykwH3sXlYw3dR11cba5d5/VRH9fx3B36Po7YxRnPlT72erp4EP/w7shi+f9lkg8sb7svi9+eJIsF/gdJLib50yQ/PUCmP0ryf5J8cfnn/KYz7Tp7UOXc73mqLN7G8HSSP09y17ozTcx1MsnnlsX9YpJ/ueY8n0zy1STfzeK3Xncn+ZUkv7Ljebp/mffPD+Lv7hqeO32ckGnX2UPbx9G6uHzM10QfdXF1uXad1Ud9XMdzp48TMu06eyBdnPhcHfo+rquLtfwfAwAAwCTr/owkAAAAP2QMSQAAAGYxJAEAAJjFkAQAAGAWQxIAAIBZDEkAAABmMSQBAACYxZAEAABgFkMSAACAWQxJAAAAZjEkAQAAmMWQBAAAYJZ9h2RVfayqXqiqL+1xf1XVR6rqYlU9WVW3rD4mkOgjjEQfYQy6CJsx5RXJB5KcepX7b09yYvnnbJLfvfZYwB4eiD7CKB6IPsIIHoguwoHbd0h296NJvvEqR04n+UQvPJbkTVX1llUFBL5PH2Ec+ghj0EXYjCMr+Bk3Jnlux/VLy9u+uvtgVZ3N4jdBef3rX/+zb3vb21bw8DCuJ5544mvdffQAH1IfYQ+j9lEXOYwOuI/+vxH2cC1dXMWQnKy7zyU5lyRbW1u9vb19kA8PB66q/nrTGfaijxw2o/ZRFzmM9BHGcC1dXMW3tj6f5PiO68eWtwEHTx9hHPoIY9BFWINVDMnzSX5p+Y1Y70zyYnf/wFsFgAOhjzAOfYQx6CKswb5vba2qTya5LckNVXUpyW8m+ZEk6e6PJrmQ5I4kF5N8K8kvryssHHb6COPQRxiDLsJm7Dsku/vMPvd3kvevLBGwJ32EcegjjEEXYTNW8dZWAAAADhFDEgAAgFkMSQAAAGYxJAEAAJjFkAQAAGAWQxIAAIBZDEkAAABmMSQBAACYxZAEAABgFkMSAACAWQxJAAAAZjEkAQAAmMWQBAAAYBZDEgAAgFkMSQAAAGYxJAEAAJjFkAQAAGAWQxIAAIBZDEkAAABmmTQkq+pUVT1TVRer6t4r3P+TVfVIVX2hqp6sqjtWHxVI9BFGoYswDn2Eg7fvkKyq65Lcn+T2JCeTnKmqk7uO/cckD3X3O5LcleS/rjoooI8wCl2EcegjbMaUVyRvTXKxu5/t7peSPJjk9K4zneTHl5ffmORvVhcR2EEfYQy6COPQR9iAKUPyxiTP7bh+aXnbTh9K8p6qupTkQpJfu9IPqqqzVbVdVduXL1++irhw6OkjjEEXYRz6CBuwqi/bOZPkge4+luSOJL9XVT/ws7v7XHdvdffW0aNHV/TQwC76CGPQRRiHPsKKTRmSzyc5vuP6seVtO92d5KEk6e4/SfJjSW5YRUDgFfQRxqCLMA59hA2YMiQfT3Kiqm6uquuz+IDy+V1n/neSdyVJVf1MFuX0fgBYPX2EMegijEMfYQP2HZLd/XKSe5I8nOTLWXzj1VNVdV9V3bk89htJ3ldVf5bkk0ne2929rtBwWOkjjEEXYRz6CJtxZMqh7r6QxQeTd972wR2Xn07yc6uNBlyJPsIYdBHGoY9w8Fb1ZTsAAAAcEoYkAAAAsxiSAAAAzGJIAgAAMIshCQAAwCyGJAAAALMYkgAAAMxiSAIAADCLIQkAAMAshiQAAACzGJIAAADMYkgCAAAwiyEJAADALIYkAAAAsxiSAAAAzGJIAgAAMIshCQAAwCyGJAAAALMYkgAAAMwyaUhW1amqeqaqLlbVvXuc+cWqerqqnqqq319tTOB79BHGoIswDn2Eg3dkvwNVdV2S+5P8iySXkjxeVee7++kdZ04k+Q9Jfq67v1lVP7GuwHCY6SOMQRdhHPoImzHlFclbk1zs7me7+6UkDyY5vevM+5Lc393fTJLufmG1MYElfYQx6CKMQx9hA6YMyRuTPLfj+qXlbTu9Nclbq+pzVfVYVZ260g+qqrNVtV1V25cvX766xHC46SOMQRdhHPoIG7CqL9s5kuREktuSnEny36vqTbsPdfe57t7q7q2jR4+u6KGBXfQRxqCLMA59hBWbMiSfT3J8x/Vjy9t2upTkfHd/t7v/KslfZFFWYLX0EcagizAOfYQNmDIkH09yoqpurqrrk9yV5PyuM/8zi9/wpKpuyOLtA8+uMCewoI8wBl2EcegjbMC+Q7K7X05yT5KHk3w5yUPd/VRV3VdVdy6PPZzk61X1dJJHkvz77v76ukLDYaWPMAZdhHHoI2xGdfdGHnhra6u3t7c38thwUKrqie7e2nSO/egjh8FroY+6yGGhjzCGa+niqr5sBwAAgEPCkAQAAGAWQxIAAIBZDEkAAABmMSQBAACYxZAEAABgFkMSAACAWQxJAAAAZjEkAQAAmMWQBAAAYBZDEgAAgFkMSQAAAGYxJAEAAJjFkAQAAGAWQxIAAIBZDEkAAABmMSQBAACYxZAEAABglklDsqpOVdUzVXWxqu59lXO/UFVdVVuriwjspI8wBl2EcegjHLx9h2RVXZfk/iS3JzmZ5ExVnbzCuTck+XdJPr/qkMCCPsIYdBHGoY+wGVNekbw1ycXufra7X0ryYJLTVzj3W0k+nOTbK8wHvJI+whh0Ecahj7ABU4bkjUme23H90vK2v1VVtyQ53t1/+Go/qKrOVtV2VW1fvnx5dlhAH2EQugjj0EfYgGv+sp2qel2S30nyG/ud7e5z3b3V3VtHjx691ocGdtFHGIMuwjj0EdZjypB8PsnxHdePLW/7njckeXuSz1bVV5K8M8l5H2KGtdBHGIMuwjj0ETZgypB8PMmJqrq5qq5PcleS89+7s7tf7O4buvum7r4pyWNJ7uzu7bUkhsNNH2EMugjj0EfYgH2HZHe/nOSeJA8n+XKSh7r7qaq6r6ruXHdA4Pv0EcagizAOfYTNODLlUHdfSHJh120f3OPsbdceC9iLPsIYdBHGoY9w8K75y3YAAAA4XAxJAAAAZjEkAQAAmMWQBAAAYBZDEgAAgFkMSQAAAGYxJAEAAJjFkAQAAGAWQxIAAIBZDEkAAABmMSQBAACYxZAEAABgFkMSAACAWQxJAAAAZjEkAQAAmMWQBAAAYBZDEgAAgFkMSQAAAGYxJAEAAJhl0pCsqlNV9UxVXayqe69w/69X1dNV9WRV/XFV/dTqowKJPsIodBHGoY9w8PYdklV1XZL7k9ye5GSSM1V1ctexLyTZ6u5/nOTTSf7TqoMC+gij0EUYhz7CZkx5RfLWJBe7+9nufinJg0lO7zzQ3Y9097eWVx9Lcmy1MYElfYQx6CKMQx9hA6YMyRuTPLfj+qXlbXu5O8lnrnRHVZ2tqu2q2r58+fL0lMD36COMQRdhHPoIG7DSL9upqvck2Ury21e6v7vPdfdWd28dPXp0lQ8N7KKPMAZdhHHoI6zOkQlnnk9yfMf1Y8vbXqGq3p3kA0l+vru/s5p4wC76CGPQRRiHPsIGTHlF8vEkJ6rq5qq6PsldSc7vPFBV70jy35Lc2d0vrD4msKSPMAZdhHHoI2zAvkOyu19Ock+Sh5N8OclD3f1UVd1XVXcuj/12kr+b5A+q6otVdX6PHwdcA32EMegijEMfYTOmvLU13X0hyYVdt31wx+V3rzgXsAd9hDHoIoxDH+HgrfTLdgAAAPjhZ0gCAAAwiyEJAADALIYkAAAAsxiSAAAAzGJIAgAAMIshCQAAwCyGJAAAALMYkgAAAMxiSAIAADCLIQkAAMAshiQAAACzGJIAAADMYkgCAAAwiyEJAADALIYkAAAAsxiSAAAAzGJIAgAAMIshCQAAwCyThmRVnaqqZ6rqYlXde4X7f7SqPrW8//NVddOqgwIL+ghj0EUYhz7Cwdt3SFbVdUnuT3J7kpNJzlTVyV3H7k7yze7++0n+S5IPrzoooI8wCl2EcegjbMaUVyRvTXKxu5/t7peSPJjk9K4zp5N8fHn500neVVW1upjAkj7CGHQRxqGPsAFHJpy5MclzO65fSvJP9jrT3S9X1YtJ3pzkazsPVdXZJGeXV79TVV+6mtBrdEN2ZR7EiLlkmuYfrPjn6eNmyTTdiLlW2cfD1MVkzL9PmaYZMVOij1dr1L/PEXPJNM1Vd3HKkFyZ7j6X5FySVNV2d28d5OPvZ8RMyZi5ZJqmqrY3nWEv+jifTNONmGvUPo7exWTMXDJNM2KmRB+v1oiZkjFzyTTNtXRxyltbn09yfMf1Y8vbrnimqo4keWOSr19tKGBP+ghj0EUYhz7CBkwZko8nOVFVN1fV9UnuSnJ+15nzSf7N8vK/TvK/urtXFxNY0kcYgy7COPQRNmDft7Yu30d+T5KHk1yX5GPd/VRV3Zdku7vPJ/kfSX6vqi4m+UYWBd7PuWvIvS4jZkrGzCXTNCvNpI8bJ9N0I+ZaWaZD1sVkzFwyTTNipkQfr9aImZIxc8k0zVVnKr+MAQAAYI4pb20FAACAv2VIAgAAMMvah2RVnaqqZ6rqYlXde4X7f7SqPrW8//NVddMAmX69qp6uqier6o+r6qc2nWnHuV+oqq6qtX918JRMVfWLy+fqqar6/XVnmpKrqn6yqh6pqi8s/w7vWHOej1XVC3v9t6Zq4SPLvE9W1S3rzPNq9HE1mXacO9R9HK2Ly8d8TfRRF1eXa8c5fdTHq6KPq8m049yBdXFqrsPex7V1sbvX9ieLDzz/ZZKfTnJ9kj9LcnLXmX+b5KPLy3cl+dQAmf55kr+zvPyrI2RanntDkkeTPJZka9OZkpxI8oUkf295/SfWmWlGrnNJfnV5+WSSr6w50z9LckuSL+1x/x1JPpOkkrwzyefX/Txdw3Onj/q4ykwH2sXl4wzfR11cba7lOX3Ux3U+d4e+jyN2ccZzdej7uK4urvsVyVuTXOzuZ7v7pSQPJjm968zpJB9fXv50kndVVW0yU3c/0t3fWl59LIv/HtE6TXmekuS3knw4ybfXnGdqpvclub+7v5kk3f3CILk6yY8vL78xyd+sM1B3P5rFN8Dt5XSST/TCY0neVFVvWWemPejjijItHfY+DtfF5DXTR11cYa4lfdTHq6WPK8q0dJBdnJrr0PdxXV1c95C8MclzO65fWt52xTPd/XKSF5O8ecOZdro7i4W+TvtmWr7EfLy7/3DNWSZnSvLWJG+tqs9V1WNVdWqQXB9K8p6qupTkQpJfO4Bcr2buv7lN5tBHfVxlpg9lrC4mY/RRF6fTx9Vl+lD08Woz6OOYXZyUK/o4xVV1cd//juRhVlXvSbKV5Oc3nON1SX4nyXs3meMKjmTxdoHbsvhN2KNV9Y+6+/9uNFVyJskD3f2fq+qfZvHfjXp7d/+/DefiGujjvkbsoy7+EBqli8ss+jidPv4QGqWPA3cx0ce1Wfcrks8nOb7j+rHlbVc8U1VHsnh59+sbzvDvXpgAACAASURBVJSqeneSDyS5s7u/s8Y8UzK9Icnbk3y2qr6SxXuXz6/5Q8xTnqdLSc5393e7+6+S/EUWRV2nKbnuTvJQknT3nyT5sSQ3rDnXq5n0b26QHPqoj6vMNFoXkzH6qIury6WP0zPp49Vn0McxuzglV6KPU1xdF/f7EOW1/MniNwDPJrk53/+w6T/cdeb9eeUHmB8aINM7sviQ7Il1ZpmTadf5z2b9XyYw5Xk6leTjy8s3ZPGS+JsHyPWZJO9dXv6ZLN53XmvOdVP2/gDzv8orP8D8pwfx7+oqnzt91MdVZjrwLi4fa+g+6uJqc+06r4/6uI7n7tD3ccQuzniu9LHX08WD+Id3RxbL/y+TfGB5231Z/PYkWSzwP0hyMcmfJvnpATL9UZL/k+SLyz/nN51p19mDKud+z1Nl8TaGp5P8eZK71p1pYq6TST63LO4Xk/zLNef5ZJKvJvluFr/1ujvJryT5lR3P0/3LvH9+EH931/Dc6eOETLvOHto+jtbF5WO+Jvqoi6vLteusPurjOp47fZyQadfZA+nixOfq0PdxXV2s5f8YAAAAJln3ZyQBAAD4IWNIAgAAMIshCQAAwCyGJAAAALMYkgAAAMxiSAIAADCLIQkAAMAshiQAAACzGJIAAADMYkgCAAAwiyEJAADALPsOyar6WFW9UFVf2uP+qqqPVNXFqnqyqm5ZfUwg0UcYiT7CGHQRNmPKK5IPJDn1KvffnuTE8s/ZJL977bGAPTwQfYRRPBB9hBE8EF2EA7fvkOzuR5N841WOnE7yiV54LMmbquotqwoIfJ8+wjj0Ecagi7AZR1bwM25M8tyO65eWt31198GqOpvFb4Ly+te//mff9ra3reDhYVxPPPHE17r76AE+pD7CHkbtoy5yGB1wH/1/I+zhWrq4iiE5WXefS3IuSba2tnp7e/sgHx4OXFX99aYz7EUfOWxG7aMuchjpI4zhWrq4im9tfT7J8R3Xjy1vAw6ePsI49BHGoIuwBqsYkueT/NLyG7HemeTF7v6BtwoAB0IfYRz6CGPQRViDfd/aWlWfTHJbkhuq6lKS30zyI0nS3R9NciHJHUkuJvlWkl9eV1g47PQRxqGPMAZdhM3Yd0h295l97u8k719ZImBP+gjj0EcYgy7CZqzira0AAAAcIoYkAAAAsxiSAAAAzGJIAgAAMIshCQAAwCyGJAAAALMYkgAAAMxiSAIAADCLIQkAAMAshiQAAACzGJIAAADMYkgCAAAwiyEJAADALIYkAAAAsxiSAAAAzGJIAgAAMIshCQAAwCyGJAAAALMYkgAAAMwyaUhW1amqeqaqLlbVvVe4/yer6pGq+kJVPVlVd6w+KpDoI4xCF2Ec+ggHb98hWVXXJbk/ye1JTiY5U1Undx37j0ke6u53JLkryX9ddVBAH2EUugjj0EfYjCmvSN6a5GJ3P9vdLyV5MMnpXWc6yY8vL78xyd+sLiKwgz7CGHQRxqGPsAFThuSNSZ7bcf3S8radPpTkPVV1KcmFJL92pR9UVWeraruqti9fvnwVceHQ00cYgy7COPQRNmBVX7ZzJskD3X0syR1Jfq+qfuBnd/e57t7q7q2jR4+u6KGBXfQRxqCLMA59hBWbMiSfT3J8x/Vjy9t2ujvJQ0nS3X+S5MeS3LCKgMAr6COMQRdhHPoIGzBlSD6e5ERV3VxV12fxAeXzu8787yTvSpKq+pksyun9ALB6+ghj0EUYhz7CBuw7JLv75ST3JHk4yZez+Marp6rqvqq6c3nsN5K8r6r+LMknk7y3u3tdoeGw0kcYgy7COPQRNuPIlEPdfSGLDybvvO2DOy4/neTnVhsNuBJ9hDHoIoxDH+HgrerLdgAAADgkDEkAAABmMSQBAACYxZAEAABgFkMSAACAWQxJAAAAZjEkAQAAmMWQBAAAYBZDEgAAgFkMSQAAAGYxJAEAAJjFkAQAAGAWQxIAAIBZDEkAAABmMSQBAACYxZAEAABgFkMSAACAWQxJAAAAZjEkAQAAmGXSkKyqU1X1TFVdrKp79zjzi1X1dFU9VVW/v9qYwPfoI4xBF2Ec+ggH78h+B6rquiT3J/kXSS4lebyqznf30zvOnEjyH5L8XHd/s6p+Yl2B4TDTRxiDLsI49BE2Y8orkrcmudjdz3b3S0keTHJ615n3Jbm/u7+ZJN39wmpjAkv6CGPQRRiHPsIGTBmSNyZ5bsf1S8vbdnprkrdW1eeq6rGqOnWlH1RVZ6tqu6q2L1++fHWJ4XDTRxiDLsI49BE2YFVftnMkyYkktyU5k+S/V9Wbdh/q7nPdvdXdW0ePHl3RQwO76COMQRdhHPoIKzZlSD6f5PiO68eWt+10Kcn57v5ud/9Vkr/IoqzAaukjjEEXYRz6CBswZUg+nuREVd1cVdcnuSvJ+V1n/mcWv+FJVd2QxdsHnl1hTmBBH2EMugjj0EfYgH2HZHe/nOSeJA8n+XKSh7r7qaq6r6ruXB57OMnXq+rpJI8k+ffd/fV1hYbDSh9hDLoI49BH2Izq7o088NbWVm9vb2/kseGgVNUT3b216Rz70UcOg9dCH3WRw0IfYQzX0sVVfdkOAAAAh4QhCQAAwCyGJAAAALMYkgAAAMxiSAIAADCLIQkAAMAshiQAAACzGJIAAADMYkgCAAAwiyEJAADALIYkAAAAsxiSAAAAzGJIAgAAMIshCQAAwCyGJAAAALMYkgAAAMxiSAIAADCLIQkAAMAsk4ZkVZ2qqmeq6mJV3fsq536hqrqqtlYXEdhJH2EMugjj0Ec4ePsOyaq6Lsn9SW5PcjLJmao6eYVzb0jy75J8ftUhgQV9hDHoIoxDH2EzprwieWuSi939bHe/lOTBJKevcO63knw4ybdXmA94JX2EMegijEMfYQOmDMkbkzy34/ql5W1/q6puSXK8u//w1X5QVZ2tqu2q2r58+fLssIA+wiB0Ecahj7AB1/xlO1X1uiS/k+Q39jvb3ee6e6u7t44ePXqtDw3soo8wBl2EcegjrMeUIfl8kuM7rh9b3vY9b0jy9iSfraqvJHlnkvM+xAxroY8wBl2EcegjbMCUIfl4khNVdXNVXZ/kriTnv3dnd7/Y3Td0903dfVOSx5Lc2d3ba0kMh5s+whh0Ecahj7AB+w7J7n45yT1JHk7y5SQPdfdTVXVfVd257oDA9+kjjEEXYRz6CJtxZMqh7r6Q5MKu2z64x9nbrj0WsBd9hDHoIoxDH+HgXfOX7QAAAHC4GJIAAADMYkgCAAAwiyEJAADALIYkAAAAsxiSAAAAzGJIAgAAMIshCQAAwCyGJAAAALMYkgAAAMxiSAIAADCLIQkAAMAshiQAAACzGJIAAADMYkgCAAAwiyEJAADALIYkAAAAsxiSAAAAzGJIAgAAMMukIVlVp6rqmaq6WFX3XuH+X6+qp6vqyar646r6qdVHBRJ9hFHoIoxDH+Hg7Tskq+q6JPcnuT3JySRnqurkrmNfSLLV3f84yaeT/KdVBwX0EUahizAOfYTNmPKK5K1JLnb3s939UpIHk5zeeaC7H+nuby2vPpbk2GpjAkv6CGPQRRiHPsIGTBmSNyZ5bsf1S8vb9nJ3ks9c6Y6qOltV21W1ffny5ekpge/RRxiDLsI49BE2YKVftlNV70myleS3r3R/d5/r7q3u3jp69OgqHxrYRR9hDLoI49BHWJ0jE848n+T4juvHlre9QlW9O8kHkvx8d39nNfGAXfQRxqCLMA59hA2Y8ork40lOVNXNVXV9kruSnN95oKrekeS/Jbmzu19YfUxgSR9hDLoI49BH2IB9h2R3v5zkniQPJ/lykoe6+6mquq+q7lwe++0kfzfJH1TVF6vq/B4/DrgG+ghj0EUYhz7CZkx5a2u6+0KSC7tu++COy+9ecS5gD/oIY9BFGIc+wsFb6ZftAAAA8MPPkAQAAGAWQxIAAIBZDEkAAABmMSQBAACYxZAEAABgFkMSAACAWQxJAAAAZjEkAQAAmMWQBAAAYBZDEgAAgFkMSQAAAGYxJAEAAJjFkAQAAGAWQxIAAIBZDEkAAABmMSQBAACYxZAEAABgFkMSAACAWSYNyao6VVXPVNXFqrr3Cvf/aFV9ann/56vqplUHBRb0EcagizAOfYSDt++QrKrrktyf5PYkJ5OcqaqTu47dneSb3f33k/yXJB9edVBAH2EUugjj0EfYjCmvSN6a5GJ3P9vdLyV5MMnpXWdOJ/n48vKnk7yrqmp1MYElfYQx6CKMQx9hA45MOHNjkud2XL+U5J/sdaa7X66qF5O8OcnXdh6qqrNJzi6vfqeqvnQ1odfohuzKPIgRc8k0zT9Y8c/Tx82SaboRc62yj4epi8mYf58yTTNipkQfr9aof58j5pJpmqvu4pQhuTLdfS7JuSSpqu3u3jrIx9/PiJmSMXPJNE1VbW86w170cT6Zphsx16h9HL2LyZi5ZJpmxEyJPl6tETMlY+aSaZpr6eKUt7Y+n+T4juvHlrdd8UxVHUnyxiRfv9pQwJ70EcagizAOfYQNmDIkH09yoqpurqrrk9yV5PyuM+eT/Jvl5X+d5H91d68uJrCkjzAGXYRx6CNswL5vbV2+j/yeJA8nuS7Jx7r7qaq6L8l2d59P8j+S/F5VXUzyjSwKvJ9z15B7XUbMlIyZS6ZpVppJHzdOpulGzLWyTIesi8mYuWSaZsRMiT5erREzJWPmkmmaq85UfhkDAADAHFPe2goAAAB/y5AEAABglrUPyao6VVXPVNXFqrr3Cvf/aFV9ann/56vqpgEy/XpVPV1VT1bVH1fVT206045zv1BVXVVr/+rgKZmq6heXz9VTVfX76840JVdV/WRVPVJVX1j+Hd6x5jwfq6oX9vpvTdXCR5Z5n6yqW9aZ59Xo42oy7Th3qPs4WheXj/ma6KMuri7XjnP6qI9XRR9Xk2nHuQPr4tRch72Pa+tid6/tTxYfeP7LJD+d5Pokf5bk5K4z/zbJR5eX70ryqQEy/fMkf2d5+VdHyLQ894YkjyZ5LMnWpjMlOZHkC0n+3vL6T6wz04xc55L86vLyySRfWXOmf5bkliRf2uP+O5J8JkkleWeSz6/7ebqG504f9XGVmQ60i8vHGb6PurjaXMtz+qiP63zuDn0fR+zijOfq0PdxXV1c9yuStya52N3PdvdLSR5McnrXmdNJPr68/Okk76qq2mSm7n6ku7+1vPpYFv89onWa8jwlyW8l+XCSb685z9RM70tyf3d/M0m6+4VBcnWSH19efmOSv1lnoO5+NItvgNvL6SSf6IXHkrypqt6yzkx70McVZVo67H0crovJa6aPurjCXEv6qI9XSx9XlGnpILs4Ndeh7+O6urjuIXljkud2XL+0vO2KZ7r75SQvJnnzhjPtdHcWC32d9s20fIn5eHf/4ZqzTM6U5K1J3lpVn6uqx6rq1CC5PpTkPVV1KcmFJL92ALlezdx/c5vMoY/6uMpMH8pYXUzG6KMuTqePq8v0oejj1WbQxzG7OClX9HGKq+rivv8dycOsqt6TZCvJz284x+uS/E6S924yxxUcyeLtArdl8ZuwR6vqH3X3/91oquRMkge6+z9X1T/N4r8b9fbu/n8bzsU10Md9jdhHXfwhNEoXl1n0cTp9/CE0Sh8H7mKij2uz7lckn09yfMf1Y8vbrnimqo5k8fLu1zecKVX17iQfSHJnd39njXmmZHpDkrcn+WxVfSWL9y6fX/OHmKc8T5eSnO/u73b3XyX5iyyKuk5Tct2d5KEk6e4/SfJjSW5Yc65XM+nf3CA59FEfV5lptC4mY/RRF1eXSx+nZ9LHq8+gj2N2cUquRB+nuLou7vchymv5k8VvAJ5NcnO+/2HTf7jrzPvzyg8wPzRApndk8SHZE+vMMifTrvOfzfq/TGDK83QqyceXl2/I4iXxNw+Q6zNJ3ru8/DNZvO+81pzrpuz9AeZ/lVd+gPlPD+Lf1VU+d/qoj6vMdOBdXD7W0H3UxdXm2nVeH/VxHc/doe/jiF2c8VzpY6+niwfxD++OLJb/Xyb5wPK2+7L47UmyWOB/kORikj9N8tMDZPqjJP8nyReXf85vOtOuswdVzv2ep8ribQxPJ/nzJHetO9PEXCeTfG5Z3C8m+ZdrzvPJJF9N8t0sfut1d5JfSfIrO56n+5d5//wg/u6u4bnTxwmZdp09tH0crYvLx3xN9FEXV5dr11l91Md1PHf6OCHTrrMH0sWJz9Wh7+O6uljL/zEAAABMsu7PSAIAAPBDxpAEAABgFkMSAACAWQxJAAAAZjEkAQAAmMWQBAAAYBZDEgAAgFkMSQAAAGYxJAEAAJjFkAQAAGAWQxIAAIBZ9h2SVfWxqnqhqr60x/1VVR+pqotV9WRV3bL6mECijzASfYQx6CJsxpRXJB9IcupV7r89yYnln7NJfvfaYwF7eCD6CKN4IPoII3gguggHbt8h2d2PJvnGqxw5neQTvfBYkjdV1VtWFRD4Pn2EcegjjEEXYTNW8RnJG5M8t+P6peVtwMHTRxiHPsIYdBHW4MhBPlhVnc3iLQV5/etf/7Nve9vbDvLh4cA98cQTX+vuo5vOcSX6yGEzah91kcNIH2EM19LFVQzJ55Mc33H92PK2H9Dd55KcS5Ktra3e3t5ewcPDuKrqrw/4IfUR9jBqH3WRw+iA++j/G2EP19LFVby19XySX1p+I9Y7k7zY3V9dwc8F5tNHGIc+whh0EdZg31ckq+qTSW5LckNVXUrym0l+JEm6+6NJLiS5I8nFJN9K8svrCguHnT7COPQRxqCLsBn7DsnuPrPP/Z3k/StLBOxJH2Ec+ghj0EXYjFW8tRUAAIBDxJAEAABgFkMSAACAWQxJAAAAZjEkAQAAmMWQBAAAYBZDkv+/vfsNtfuu7wD+/qyhyljVYSOMJtrK0mHmBsqlcwjToRtpB80Dh7RQmFAsulUGyqDDIaU+cjIHg4ytY6ITtFYfjICVwlxFEFubUa22pRKjrKljzbrOJ6K17LMH91Rv75Le3y85555ve18vKJw/397z5iTvB+97/gQAAGAWQxIAAIBZDEkAAABmMSQBAACYxZAEAABgFkMSAACAWQxJAAAAZjEkAQAAmMWQBAAAYBZDEgAAgFkMSQAAAGYxJAEAAJjFkAQAAGCWSUOyqo5U1aNVdbKqbjnL/a+uqnuq6oGqerCqrll+VCDRRxiFLsI49BF2345DsqouSnIsydVJDie5vqoObzv2F0nu7O43JLkuyd8uOyigjzAKXYRx6COsx5RXJK9KcrK7T3X300nuSHJ025lO8rLF5Zcn+cHyIgJb6COMQRdhHPoIazBlSF6W5LEt108vbtvq1iQ3VNXpJHcled/ZflBV3VRVJ6rqxJkzZ84jLux5+ghj0EUYhz7CGizry3auT/KJ7j6Q5Jokn6qq//ezu/v27t7o7o39+/cv6aGBbfQRxqCLMA59hCWbMiQfT3Jwy/UDi9u2ujHJnUnS3V9L8tIkly4jIPAc+ghj0EUYhz7CGkwZkvcnOVRVV1TVxdn8gPLxbWf+PcnbkqSqXpfNcno/ACyfPsIYdBHGoY+wBjsOye5+JsnNSe5O8kg2v/Hqoaq6raquXRz7QJJ3V9U3k3wmybu6u1cVGvYqfYQx6CKMQx9hPfZNOdTdd2Xzg8lbb/vQlssPJ3nzcqMBZ6OPMAZdhHHoI+y+ZX3ZDgAAAHuEIQkAAMAshiQAAACzGJIAAADMYkgCAAAwiyEJAADALIYkAAAAsxiSAAAAzGJIAgAAMIshCQAAwCyGJAAAALMYkgAAAMxiSAIAADCLIQkAAMAshiQAAACzGJIAAADMYkgCAAAwiyEJAADALIYkAAAAs0waklV1pKoeraqTVXXLOc68s6oerqqHqurTy40JPEsfYQy6COPQR9h9+3Y6UFUXJTmW5PeSnE5yf1Ud7+6Ht5w5lOTPk7y5u5+qqletKjDsZfoIY9BFGIc+wnpMeUXyqiQnu/tUdz+d5I4kR7edeXeSY939VJJ09xPLjQks6COMQRdhHPoIazBlSF6W5LEt108vbtvqyiRXVtVXq+reqjpyth9UVTdV1YmqOnHmzJnzSwx7mz7CGHQRxqGPsAbL+rKdfUkOJXlrkuuT/ENVvWL7oe6+vbs3untj//79S3poYBt9hDHoIoxDH2HJpgzJx5Mc3HL9wOK2rU4nOd7dP+3u7yX5TjbLCiyXPsIYdBHGoY+wBlOG5P1JDlXVFVV1cZLrkhzfduafs/kbnlTVpdl8+8CpJeYENukjjEEXYRz6CGuw45Ds7meS3Jzk7iSPJLmzux+qqtuq6trFsbuTPFlVDye5J8mfdfeTqwoNe5U+whh0Ecahj7Ae1d1reeCNjY0+ceLEWh4bdktV/Vt3b6w7x070kb3ghdBHXWSv0EcYw4V0cVlftgMAAMAeYUgCAAAwiyEJAADALIYkAAAAsxiSAAAAzGJIAgAAMIshCQAAwCyGJAAAALMYkgAAAMxiSAIAADCLIQkAAMAshiQAAACzGJIAAADMYkgCAAAwiyEJAADALIYkAAAAsxiSAAAAzGJIAgAAMMukIVlVR6rq0ao6WVW3PM+5d1RVV9XG8iICW+kjjEEXYRz6CLtvxyFZVRclOZbk6iSHk1xfVYfPcu6SJH+a5L5lhwQ26SOMQRdhHPoI6zHlFcmrkpzs7lPd/XSSO5IcPcu5Dyf5SJIfLzEf8Fz6CGPQRRiHPsIaTBmSlyV5bMv104vbfqaq3pjkYHd/4fl+UFXdVFUnqurEmTNnZocF9BEGoYswDn2ENbjgL9upql9I8rEkH9jpbHff3t0b3b2xf//+C31oYBt9hDHoIoxDH2E1pgzJx5Mc3HL9wOK2Z12S5PVJvlxV30/ypiTHfYgZVkIfYQy6COPQR1iDKUPy/iSHquqKqro4yXVJjj97Z3f/sLsv7e7Lu/vyJPcmuba7T6wkMext+ghj0EUYhz7CGuw4JLv7mSQ3J7k7ySNJ7uzuh6rqtqq6dtUBgZ/TRxiDLsI49BHWY9+UQ919V5K7tt32oXOcfeuFxwLORR9hDLoI49BH2H0X/GU7AAAA7C2GJAAAALMYkgAAAMxiSAIAADCLIQkAAMAshiQAAACzGJIAAADMYkgCAAAwiyEJAADALIYkAAAAsxiSAAAAzGJIAgAAMIshCQAAwCyGJAAAALMYkgAAAMxiSAIAADCLIQkAAMAshiQAAACzGJIAAADMMmlIVtWRqnq0qk5W1S1nuf/9VfVwVT1YVV+qqtcsPyqQ6COMQhdhHPoIu2/HIVlVFyU5luTqJIeTXF9Vh7cdeyDJRnf/ZpLPJ/nLZQcF9BFGoYswDn2E9ZjyiuRVSU5296nufjrJHUmObj3Q3fd0948WV+9NcmC5MYEFfYQx6CKMQx9hDaYMycuSPLbl+unFbedyY5Ivnu2Oqrqpqk5U1YkzZ85MTwk8Sx9hDLoI49BHWIOlftlOVd2QZCPJR892f3ff3t0b3b2xf//+ZT40sI0+whh0Ecahj7A8+yaceTzJwS3XDyxue46qenuSDyZ5S3f/ZDnxgG30EcagizAOfYQ1mPKK5P1JDlXVFVV1cZLrkhzfeqCq3pDk75Nc291PLD8msKCPMAZdhHHoI6zBjkOyu59JcnOSu5M8kuTO7n6oqm6rqmsXxz6a5JeSfK6qvlFVx8/x44ALoI8wBl2EcegjrMeUt7amu+9Kcte22z605fLbl5wLOAd9hDHoIoxDH2H3LfXLdgAAAHjxMyQBAACYxZAEAABgFkMSAACAWQxJAAAAZjEkAQAAmMWQBAAAYBZDEgAAgFkMSQAAAGYxJAEAAJjFkAQAAGAWQxIAAIBZDEkAAABmMSQBAACYxZAEAABgFkMSAACAWQxJAAAAZjEkAQAAmMWQBAAAYJZJQ7KqjlTVo1V1sqpuOcv9L6mqzy7uv6+qLl92UGCTPsIYdBHGoY+w+3YcklV1UZJjSa5OcjjJ9VV1eNuxG5M81d2/muSvk3xk2UEBfYRR6CKMQx9hPaa8InlVkpPdfaq7n05yR5Kj284cTfLJxeXPJ3lbVdXyYgIL+ghj0EUYhz7CGuybcOayJI9tuX46yW+d60x3P1NVP0zyyiT/tfVQVd2U5KbF1Z9U1bfPJ/QKXZptmQcxYi6Zpvm1Jf88fVwvmaYbMdcy+7iXupiM+ecp0zQjZkr08XyN+uc5Yi6ZpjnvLk4ZkkvT3bcnuT1JqupEd2/s5uPvZMRMyZi5ZJqmqk6sO8O56ON8Mk03Yq5R+zh6F5Mxc8k0zYiZEn08XyNmSsbMJdM0F9LFKW9tfTzJwS3XDyxuO+uZqtqX5OVJnjzfUMA56SOMQRdhHPoIazBlSN6f5FBVXVFVFye5LsnxbWeOJ/mjxeU/TPKv3d3Liwks6COMQRdhHPoIa7DjW1sX7yO/OcndSS5K8vHufqiqbktyoruPJ/nHJJ+qqpNJ/jubBd7J7ReQe1VGzJSMmUumaZaaSR/XTqbpRsy1tEx7rIvJmLlkmmbETIk+nq8RMyVj5pJpmvPOVH4ZAwAAwBxT3toKAAAAP2NIAgAAMMvKh2RVHamqR6vqZFXdcpb7X1JVn13cf19VXT5ApvdX1cNV9WBVfamqXrPuTFvOvaOquqpW/tXBUzJV1TsXz9VDVfXpVWeakquqXl1V91TVA4s/w2tWnOfjVfXEuf6tqdr0N4u8D1bVG1eZ5/no43IybTm3p/s4WhcXj/mC6KMuLi/XlnP6qI/nRR+Xk2nLuV3r4tRce72PK+tid6/sv2x+4Pm7SV6b5OIk30xyeNuZP07yd4vL1yX57ACZfjfJLy4uv3eETItzlyT5SpJ7k2ysO1OSQ0keSPLLi+uvWmWmGbluT/LexeXDSb6/4ky/k+SNSb59jvuvSfLFJJXkTUnuW/XzdAHPnT7q4zIz7WoXF48zfB91cbm5Fuf0UR9X+dzt+T6O2MUZz9We7+OqXRq3ZwAAA2VJREFUurjqVySvSnKyu09199NJ7khydNuZo0k+ubj8+SRvq6paZ6buvqe7f7S4em82/z2iVZryPCXJh5N8JMmPV5xnaqZ3JznW3U8lSXc/MUiuTvKyxeWXJ/nBKgN191ey+Q1w53I0yT/1pnuTvKKqfmWVmc5BH5eUaWGv93G4LiYvmD7q4hJzLeijPp4vfVxSpoXd7OLUXHu+j6vq4qqH5GVJHtty/fTitrOe6e5nkvwwySvXnGmrG7O50Fdpx0yLl5gPdvcXVpxlcqYkVya5sqq+WlX3VtWRQXLdmuSGqjqd5K4k79uFXM9n7t+5debQR31cZqZbM1YXkzH6qIvT6ePyMt0afTzfDPo4Zhcn5Yo+TnFeXdzx35Hcy6rqhiQbSd6y5hy/kORjSd61zhxnsS+bbxd4azZ/E/aVqvqN7v6ftaZKrk/yie7+q6r67Wz+u1Gv7+7/XXMuLoA+7mjEPurii9AoXVxk0cfp9PFFaJQ+DtzFRB9XZtWvSD6e5OCW6wcWt531TFXty+bLu0+uOVOq6u1JPpjk2u7+yQrzTMl0SZLXJ/lyVX0/m+9dPr7iDzFPeZ5OJzne3T/t7u8l+U42i7pKU3LdmOTOJOnuryV5aZJLV5zr+Uz6OzdIDn3Ux2VmGq2LyRh91MXl5dLH6Zn08fwz6OOYXZySK9HHKc6vizt9iPJC/svmbwBOJbkiP/+w6a9vO/Mnee4HmO8cINMbsvkh2UOrzDIn07bzX87qv0xgyvN0JMknF5cvzeZL4q8cINcXk7xrcfl12Xzfea041+U59weY/yDP/QDz13fj79V5Pnf6qI/LzLTrXVw81tB91MXl5tp2Xh/1cRXP3Z7v44hdnPFc6WOvpou78Rfvmmwu/+8m+eDittuy+duTZHOBfy7JySRfT/LaATL9S5L/TPKNxX/H151p29ndKudOz1Nl820MDyf5VpLrVp1pYq7DSb66KO43kvz+ivN8Jsl/JPlpNn/rdWOS9yR5z5bn6dgi77d248/uAp47fZyQadvZPdvH0bq4eMwXRB91cXm5tp3VR31cxXOnjxMybTu7K12c+Fzt+T6uqou1+J8BAABgklV/RhIAAIAXGUMSAACAWQxJAAAAZjEkAQAAmMWQBAAAYBZDEgAAgFkMSQAAAGb5P59g0GzfVBs+AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# ------------------------------------------------------------------\n", "# Exercise 1\n", "nvar = 4\n", "f, axarr = plt.subplots(nvar, nvar, figsize=[12.8, 9.6])\n", "plt.tight_layout()\n", "nbins = 20\n", "# add your code here ...\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 9.2: Naive Likelihood-Ratio (obligatory)\n", "In the obligatory parts of this exercise, we only use the first two variables in each line of the data-set, the sepal length and width, to separate one of the species, the 'signal', from the other two, the 'background' in order to simplify the visualization of the results.\n", "\n", "Study the example code, which contains a general python class `Plotter`, which provides methods to display several aspects for the evaluation of the performance of a classifier. The classifier itself is to be implemented as another class containing a method `evaluate()` (an example classifier class is given [below](#cutclass)), which returns the value of the test-statistic under study and is used by the class `Plotter`. The method `plot_contour()` implemented in this class draws scatter plots of the data superimposed with contours representing the result of evaluation of the test-statistic. Also included is a method `plot_test_statistic()` to display the distributions of the test-statistic on background and signal samples, and a method `plot_roc()` to visualise the performance, the so-called 'ROC' curve.\n", "" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# ------- helper functions ----------------------------------\n", "\n", "\n", "class Plotter(object):\n", " \"\"\"\n", " class to display and evaluate the performance of a test-statistic\n", " \"\"\"\n", "\n", " def __init__(self, signal_data, bckgrd_data):\n", " self.signal_data = signal_data\n", " self.bckgrd_data = bckgrd_data\n", " self.data = np.vstack([signal_data, bckgrd_data])\n", "\n", " def plot_contour(self, classifier):\n", " # 1st variable as x-dimension in the plots\n", " xdim = 0\n", " # and 2nd variable as y-dimension\n", " ydim = 1\n", " # Draw the scatter-plots of signal and background\n", " ax = plt.subplot()\n", " ax.scatter(self.signal_data[:, xdim], self.signal_data[:, ydim],\n", " c='r', label='Signal')\n", " ax.scatter(self.bckgrd_data[:, xdim], self.bckgrd_data[:, ydim],\n", " c='b', label='Background')\n", "\n", " # Evaluate the response function on a two-dimensional grid ...\n", " # ... using the mean-values of the data for the remaining dimensions.\n", " xs = np.arange(min(self.data[:, xdim]) - 1, max(self.data[:, xdim]) + 1, 0.1)\n", " ys = np.arange(min(self.data[:, ydim]) - 1, max(self.data[:, ydim]) + 1, 0.1)\n", "\n", " means = np.mean(self.data, axis=0) # calculate mean of each column\n", " responses = np.zeros((len(ys), len(xs)))\n", " for i, x in enumerate(xs):\n", " for j, y in enumerate(ys):\n", " values = np.copy(means)\n", " values[xdim] = x\n", " values[ydim] = y\n", " responses[j, i] = float(classifier.evaluate(values))\n", "\n", " # Draw a contour plot\n", " X, Y = np.meshgrid(xs, ys)\n", " c = ax.contourf(X, Y, responses, alpha=0.5, cmap=plt.cm.coolwarm)\n", " cbar = plt.colorbar(c, orientation='vertical', ax=ax)\n", " # add the direction of the fisher vector (if specified)\n", " if hasattr(classifier, 'fisher_vector'):\n", " vector = classifier.fisher_vector / np.linalg.norm(classifier.fisher_vector)\n", " ax.set_aspect('equal', 'datalim')\n", " ax.plot([-vector[xdim] + means[xdim], vector[xdim] + means[xdim]],\n", " [-vector[ydim] + means[ydim], vector[ydim] + means[ydim]],\n", " 'k-', lw=4, label=\"Fisher Projection\")\n", " plt.title(\n", " \"scatter plot {} vs {} and classifier contour\".format(\n", " columns[xdim], columns[ydim]))\n", " # cbar.draw_all()\n", " plt.show()\n", "\n", " def plot_test_statistic(self, classifier):\n", " # Draw Distribution of the test-statistic\n", " ns, binss, _ = plt.hist(list(map(classifier.evaluate, self.signal_data)),\n", " color='r', alpha=0.5, label='Signal')\n", " nb, binsb, _ = plt.hist(list(map(classifier.evaluate, self.bckgrd_data)),\n", " color='b', alpha=0.5, label='Background')\n", " plt.title(\"test statistic\")\n", " plt.show()\n", "\n", " # calculate efficiencies and plot ROC-curves\n", " def plot_roc(self, classifier):\n", " ns, binss = np.histogram(list(map(classifier.evaluate, self.signal_data)))\n", " nb, binsb = np.histogram(list(map(classifier.evaluate, self.bckgrd_data)))\n", " # enforce common binning for response on bkg and sig\n", " minresp = min([binss[0], binsb[0]])\n", " maxresp = max([binss[len(binss) - 1], binsb[len(binsb) - 1]])\n", " nbins = 100\n", " bins = np.linspace(minresp, maxresp, nbins)\n", " bwid = (maxresp - minresp) / nbins\n", " # calculate cumulative distributions (i.e. bkg and sig efficiencies)\n", " h, b = np.histogram(list(map(classifier.evaluate, self.signal_data)), bins, density=True)\n", " ns = np.cumsum(h) * bwid\n", " h, b = np.histogram(list(map(classifier.evaluate, self.bckgrd_data)), bins, density=True)\n", " nb = np.cumsum(h) * bwid\n", " # finally, draw bkg-eff vs. sig-eff\n", " f2, ax = plt.subplots(1, 1)\n", " ax.plot(1. - ns, nb, 'r-', 1. - ns, nb, 'bo', linewidth=2.0)\n", " ax.set_xlabel(\"signal efficiency\")\n", " ax.set_ylabel(\"background rejection\")\n", " ax.set_title(\"ROC curve\")\n", " plt.show()\n", "\n", "\n", "def find_bin(x, edges):\n", " \"\"\"\n", " returns the bin number (in array of bin edges) corresponding to x\n", " @param x: value for which to find correspoding bin number\n", " @param edges: array of bin edges\n", " \"\"\"\n", " return max(min(bisect.bisect(edges, x) - 1, len(edges) - 2), 0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To guide your implementation of the naive likelihood-ratio classifier an implementation example of a classifier class is given below. The class `CutClassifier` implements a simple, cut-based analysis, where the method `evaluate(x)` simply returns the number of cuts satisfied by the point in variable space given in the array $x$.\n", "" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAEICAYAAACpqsStAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nO2de5hcVZW33193OgmXpBuSAIFcOkzCTRICCRcHBzARhITAzKgjeAH8dFAfGWXUccT5BI3jjDo3cFAhHyLgBVQECRBFJhGBUZAOhLSAQEwCSSDk3kmAdJLu9f2xT6VPnz5Vdaqr+lLd632efrrO3vvsvc5t1am19lpbZobjOI5TndT0tQCO4zhO93El7jiOU8W4Enccx6liXIk7juNUMa7EHcdxqhhX4o7jOFXMoFTikholmaQhvTDWMEnPShpbZj+ZZJZ0lqS1Ber/Q9LHy5Glp5D0r5KurEA/D0n6SIZ2qyW9I0/dNEm/LVeWrEi6TNKj3djvFkn/3BMyRf3vlHRk9Hk/SfdKapH0U0nvl/SrnhrbyUavK3FJX5L0g0RZpoeuL0iTN6VNXmUAXA48bGavRm07PXSS3iLpVUmfrZzUBfl34AuShvbUAJIulvRcouzBPGWfjz6PAS4Bboy2O30ZSRoq6S5J/ytpZE/JnsPMlgPbJM3r6bH6M2Z2oJmtjDbfDRwKjDKz95jZD83snD4Uryjd/XKsJgbEm7ik2r6WoQAfA76fViHpRODXwD+b2b/3hjDRl8kfgQt6cJiHgWMixUz06+EEYL9E2VujtgCXAYvM7M1kZ5KGAXcBDcA5Zra9B2WP80Pgo700VjUwEXjBzPaW21E/f2YzoUDf61AzK/gH/COwDtgBPA/MjsprgS8Af4rqlgLjo7rrgDXA9qj8L6Lyc4HdwB5gJ/A08FWgDdgVlV0ftT0GeBDYEo37NzGZbgG+AywCXgfekSL3Q8C/Ar+P5LgHODiqawQMGBJtHw4sjMZaAfxtPnnznKPVeWSYALyZGycm+z8DpwCbgI/E6mqAz0fndDPwkwIyHwx8D3gF2Ar8PCo/C1gLfAbYALwKfCgh1z8B38tzLL8ArkiUPQ38NSDgv6J+twPNwPF5+vkT8K7o8ymEL6tbE2VvAHXR9hLgA7H9c8exP/Ar4JfAfrH604DfAtsi+c5KXPv4ef1b4DnCffoscFLsun0WWA60AD8Ghsf2OyK6fsPyHOOHYv2uBD6aIn/qdQBGEe657YR79CvAowWew7fFjncNcFn8foo+HwTcB2yM7on7gHGxPi6L5NwBrALeH5VPBn4TnYNNwI9j+1hU/2U6Pwsfjvp7NNa23Gc29Z6OXcMVUd8LgcMTMn4MeDE6P98i3KvHEvRKWyTztqh9PXBbdJ5eAv4vUBPVfQn4QazvRjo/dw8RdNb/RvfG5GI6tKf/iinwo6Mb5vDYAf1Z9PkfCA/x0dEJO4HwMwvgA9FNOoRwE68nejiSJynPQ3dANO6Hoj5OjG6u42I3RAtwOkHxDU+R/SHCl8/xUX8/y42bcmEeBr4NDAemRxd3Vj55U8ZaneemnAs8kyi7haCUtgAfTNR9CngMGAcMI5gWbs8j8/0EpXMQUAecGVMee4H5UfkcgrI8KDbOXwNP5jmWS4D/jW0fR3gwhgHvJHwpN9DxkIzN08/3gOuiz5+N5PnbRNmSWPuNwMmx7bOist8QHtphsbojCF9yc6Lrf3a0PSZ5PwHvie6DkyOZJwMTY9ft94Qv8YMJCvljiePYDkzLc4xzgT+L+j0zOs8nxeTPex2AOwhf0gcQ7tF15FHihDfgHcDFUV+jgOmx+ymnxEcB7yJ88Y0AfkrHl/sB0bEcHW2PBd4Sfb6d8MVeQ3gG3hYb24gUFV0V3GU5manMM5vvnp4V9XUS4T78b4KJMi7jfYT7cgLhvjk3KWOs/W2El7oRhOfqBeDDeY6xka5K/GXgLdFx1nVX+Vbqr5gSn0x4i3hHUljCN+2FmQYJ36onpJ2k5EMXbb8XeCTR5kbgmtgNcVuRMR8CvhbbPo7wJlEbvzDAeMI39YhY238Fbsknb8pYq0lX4u8HHkuU3UJ4mFYBoxN1zxH90ok9aHsiOeMyjwXaiSnm2D5n0fXtfwNwWmz7bGBlnmMZQXhTmhhtfxW4OfYwvUB4C64pck4uA56KPt8TjXlMouyaWPs9wDGJ49gVXbN3Jfr+R+D7ibIHgEuT91NU/qkC1y3+9v8N4IZEm3XAGRnv85/nxip0HaJ7MHm8/0J+JX4VcHeeuluIlHhK3XRga/T5AMKX8buI/aKJ6m4DFhB7a4/VZVXiZT2zRe7p7wLfiG0fGJ2/xpiM8S+enwCfT8oYbddG99RxsbKPAg/lOcZGuirx+Vnuh976K2jPMbMVwJXRgW2QdIekw6Pq8YSfzF2Q9FlJz0Ve7G2Eny+jC42VYCJwqqRtuT+CQjws1mZNhn7ibV4ifLsn5Tgc2GJmOxJtjyhB3nxsJSjFJN8CmoAHJR0UK58I3B075ucIXzCHJvYfH8m8Nc+4m62z3fINwo2fYwThge5CdB7uBy6Kii4m2IYxsyXA9ZH8GyQtKOBkfBiYFh3facDvzOyPwNio7G102MMh/VxtiuS4VdI7Y+UTgfck7o+3ERRBkrz3acT62OfkeYIC50rSeZIek7QlkmEOne+vfNdhDOHLOHl/5qPYMeTk2V/SjZJekrSdcH4bJNWa2esERfsx4FVJ90s6Jtr1c4RfE7+X9Iyk/1NsrBTKfWYL3dOHEzs/ZraT8Msr/owWu445RhP0QPx8l/q8Z9E9vUZRo7yZ/cjM3ka4SAZ8PapaQ/gp2QlJf0G4Kf6G8K3aQPgZpVyXacMkttcAvzGzhtjfgWb28QL7pDE+9nkC4dt7U6LNK8DBkkYk2q4rYZx8LAcmpUwLbAPeR/hZ9kBMEa4Bzksc93AzW5fYf00kc0M35TqWYEfOx+3AxZLeSvh5/etchZl908xmEH7ZHEUwq3XBwoyGVwizc16OHjyA30VlBxJMRzmWR/0l+7mLYIa5U9Lbo+I1hDfx+Hk6wMy+liJK6n2aBUlHAEMJvzqTdcMIJrp/Bw6N7vNFdNznhdhIMLUk7898ZD2GzxDMm6ea2UjgjJy4AGb2gJmdTfiy+yPw/6Ly9Wb2t2Z2OOGt9NuSJmcYLyljOc9soXv6FYL+CQcjHUAwHSWfizSSY24i6IGJsbL48/46wRyVI/4llK/PPqWgEpd0tKRZ0Q27i/DzsD2qvgn4iqQpkZd2mqRRhDeXvYQbdYikq4H429prQGPCq/sacGRs+z7gKEkflFQX/Z0s6dgSj+8Dko6TtD/BNnmnmbXFG5jZGoLD6F8lDZc0jeC0yU0rTJM3jbpo/9zfEDNbS3DGnJJsbGZ7CPbaTcCi6Ma8AfiqpIkQpt1JujBl31cJDshvSzooOj9nJNsV4Mxo/3wsItzk8wlOrvZInpMlnSqpjnCz76LjfkjjEeDT0f8cj0ZlTdZ5JsqiSK4umNntwBXAPZJOJ1ybeZLeKak2Ot9nSRqXsvtNwGclzYju08m585uBMwl2+9aUuqEE++xGYK+k84BM0+2ie/Au4EvR2/NxwKUFdvkh8A5JfyNpiKRRkqantBtBeEa3SToYuCZXIelQSRdG91krwdGXu67viZ27rQQlVei6plHWM1vknr4d+JCk6ZEu+hfgcTNbnaHr14BxiqbURuf+J4TnbER0L3yajud9GXCGpAmS6gmmrH5NMcU0DPgaQdGsBw6h46D+k3AyfkWw8X4X2I9gg/wlwXb6EuFBj//8+Gn0f7OkJ6PP1wHvlrRV0jejn/TnEH5KvxKN/fVInlL4PsEWt57wRvnJPO0uJti+XgHuJtjx/qeAvGksIjxAub8vReU3Ah9M28HMdhOcjLuAewl2yYXAryTtILypnppnvA8S3ij+SLC1ZgqSUQg6Oo5gv00lUlp3EXwhP4pVjSS8vW0lXNvNwL8VGO43hHsmPk/3kajs4UTb24A5kvbLI9OthDfN+wlvkhcSZkdtJNxf/0DK/WxmPyXY9X9EcA7+nODEzML7CV+safLsINxPPyGcj/cRrl1WriD8GllPuEe/l6+hmb1MMNV8huAQX0aYSJDkWsIzuIlw7/wyVldDUFavRH2cCeTekk8GHpe0MzqGT1nH3PBMVOiZTb2no2fxi4RfPq8SfpVclKePJEuAZ4D1knK/wv+O8BKyknBv/gi4ORrrQYJzdTnBiX9fCfL3CYqM9QMOSQ8RHBQ39bEcw4CnCA7LV/tSlkie/wD+ZGbf7mtZkkj6F2CDmV3bD2SZBtxoZm/ta1kcpxCuxB3HcaqYvo82chzHGUBIGi/p1wo5k56R9KmUNpL0TUkrJC2XdFKs7lJJL0Z/hXwlof1AfRN3HMfpCyK/01gzezKa9bYU+EszezbWZg7BNj+H4Pe6zsxOjRzSTcBMgoN5KTCjwHRifxN3HMepJGb2qpk9GX3eQYj3SM5Dv5AQ/GRm9hhhPv9YQlT0g2aWmzP/ICH9R156PBVrdzl41GgbP77Q1FlnILCnDaQaWvdAa2sbdbX+y9ApnVUrnt5kZmO6u/9Zxx9vW3fuLN4QWP7SS88QZpTlWGBmC9LaSmokpCB4PFF1BJ1n7a2NyvKV56XfKvHx4yfwi/8Z0BkkHWB9C9QO2Y+Vr7SzevXrHDZid1+L5FQh75s3plDEa1G27tzJoi/+U6a24z5y+S4zm1msnaQDCdMir7QezLzp5hTHcZwKEwXE/Qz4YRR1nGQdnSN2x0Vl+crz4krccRyngkgSIfjxOTP7zzzNFgKXRLNUTgNaojiSB4BzoqjVgwgBVA8UGq/fmlMcx3GqlNMJ0afNkpZFZV8gyo9jZjcQIrznENJyvEFI4YuZbZH0FeCJaL/5Zral0GCuxB3HcSqImT1KkURoFuZ2fyJP3c1EaQCy4OYUx3GcKsaVuOM4ThXjStxxHKeKcSXuOI5TxbgSdxzHqWJciTuO41QxFVHiklZLapa0TFJTSn3etIuO4zhO96nkPPG3m1lyEeIc5wFTor9Tge+Qf9kxx3EcJyO9ZU7Jl3bRcRzHKYNKKXEjLO67VNLlKfUlp1d0HMdxilMpc8rbzGydpEOAByX90cySq5kXJfoCuBzgiHHji7R2HMdxKvImbmbrov8bgLuBUxJNMqVXNLMFZjbTzGaOGjW6EqI5juMMaMpW4pIOiNaRQ9IBhNSJf0g0y5d20XEcxymDSphTDgXuDil0GQL8yMx+KeljUDjtouM4jlMeZStxM1sJnJBSfkPsc960i47jOE738YhNx3GcKsaVuOM4ThXjStxxHKeK8eXZHMdxKoikm4HzgQ1mdnxK/T8A7482hwDHAmOi9TVXAzuANmCvmc0sNp6/iTuO41SWW4Bz81Wa2b+Z2XQzmw5cBfwmsRjy26P6ogocXIk7juNUlChaveAK9TEuBm4vZzw3pziOM+hpb93F7lXP9+qYkvYnvLFfESvO5aEy4EYzW1CsH1fijuMMeupGHsghs/4iW+Ov/sfoxLoJC7Io2xTmAf+bMKWUnIfKlbjjOE5pbMpqry7CRSRMKfE8VJJyeagKKnG3iTuO4/QykuqBM4F7YmVZ8lB1wd/EHcdxKoik24GzgNGS1gLXAHXQKR3JXwG/MrPXY7um5qEqNp4rccdxnApiZhdnaHMLYSpivCw1D1Ux3JziOI5TxbgSdxzHqWJciTuO41QxrsQdx3GqGFfijuM4VUxFlLikWklPSbovpe4ySRslLYv+PlKJMR3HcZzKTTH8FPAcMDJP/Y/N7Io8dY7jOE43qcRq9+OAucBN5YvjOI7jlEIlzCnXAp8D2gu0eZek5ZLulDQ+XyNJl0tqktS0efOmCojmOI4zsClLiUvKrV6xtECze4FGM5sGPAjcmq+hmS0ws5lmNnPUqNHliOY4jjMoKNcmfjpwgaQ5wHBgpKQfmNkHcg3MbHOs/U3AN8oc06kympth8RJoaYH6epg9C6ZO7WupHGdgUNabuJldZWbjzKyRkFZxSVyBA0gaG9u8gOAAdQYJzc1w773Qsg2w8P/ee0O54zjl0yPzxCXNl3RBtPlJSc9Iehr4JHBZT4zp9E8WL4E9ezqX7dkTyh3HKZ+KZTE0s4eAh6LPV8fKryIsBuoMQlpaSit3HKc0PGLT6VHq60srdxynNFyJOz3K7FlQV9e5rK4ulDuOUz6+KITTo+RmofjsFMfpGVyJOz3O1KmutB2np3BziuM4ThXjStxxHKeCSLpZ0gZJqSvVSzpLUksss+vVsbpzJT0vaYWkz2cZz80pTrfxSEzHSeUW4HrgtgJtHjGz8+MFkmqBbwFnA2uBJyQtNLNnCw3mStzpFrlIzFwgTy4SE1yRO1XI3t20b1pbka7M7GFJjd3Y9RRgRbTqPZLuAC4EXIk7ladQJKYrcafa0H4HMHzazKzNR0tqim0vMLMFJQ751iiK/RXgs2b2DHAEsCbWZi1warGOXIk73cIjMZ1BzCYzy6zxU3gSmGhmO6PkgT8HpnS3M3dsOt3CIzEdp3uY2XYz2xl9XgTUSRoNrAPi6y2Mi8oK4m/ig5zuOidnz+psEwePxHScLEg6DHjNzEzSKYSX6c3ANmCKpEkE5X0R8L5i/bkSH8SU45z0SEzHSUfS7cBZBNv5WuAaoA7AzG4A3g18XNJe4E3gIjMzYK+kK4AHgFrg5shWXhBX4oOYcp2THonpOF0xs4uL1F9PmIKYVrcIWFTKeG4TH8S4c9Jxqh9X4oMYd046TvVTMSUuqVbSU5LuS6kbJunHUSjp492cCO9UmL5ME9vcDNdeBw//Bq67FjZs6PkxHWcgUsk38U+Rf/3MDwNbzWwy8F/A1ys4rtNNpk6FefOgvgFQ+D9vXs/buePrbhqwrQWefwFee61nx3WcgUhFHJuSxgFzga8Cn05pciHwpejzncD1khR5ZJ0+pC+ck2kO1fY2WLUWTpjcu7I4TrVTqTfxa4HPAe156veFk5rZXqAFGJVsJOlySU2SmjZv3lQh0Zz+Rj7HaWtr78rhOAOBspW4pPOBDWa2tNy+zGyBmc00s5mjRo0utzunn5LPcTpsWO/K4TgDgUqYU04HLohyAAwHRkr6gZl9INYmF066VtIQoJ4QoeRUCfcvgqVLwdpBNTBjBsyd072+0qI9a2ph0qTKyOo4g4my38TN7CozG2dmjYQw0SUJBQ6wELg0+vzuqI3bw6uE+xdB0xNBgUP43/REKO8OcYeqgIZ6OPooOPTQionsOIOGHovYlDQfaDKzhcB3ge9LWgFsISh7p0pYmsdQtnRp99/Gcw7V9S3w9tmw8hVYvbrbIjrOoKWiStzMHgIeij5fHSvfBbynkmM5vYflcVfnK3ccp/fwiE2nKMpzl+Qrdxyn9/AEWIOI274Pq1Z2bE86Ei75YPH9ZswINvC08iz4Wpx9z/LlsCR2DWbNgmnTut/O6T/4u9QgIanAIWzf9v3i+86dAzNP7njzVk3YzmIPj0dnYh3pbpubSz4Ep5ssXw733dcxP7+lJWwvX969dk7/wt/EBwlJBV6sPMncOd1zYvpanH3PkjzXYMmSzm/ZWds5/Qt/E3d6FE932/dkvQZ+raoTV+JOj+LpbvuerNfAr1V14uaUQcKkI9NNJ5OO7FqW5ogEX4uzWpk1K9i2k9dg1qzutXP6F67EBwknTodVqwi5X3MolMdJW3fznnvCbu1tHWW+Fmf1kLNnF5t1krWdUxhJNwO5nFLHp9S/H/hHQsDyDuDjZvZ0VLc6KmsD9prZzGLjuRIfJCxeQmcFTthOOhjTHJFtbV3787U4q4tp07Ip46ztnILcQlhD87Y89auAM81sq6TzgAXAqbH6t5tZ5jSursQHCeU6t0rp03GqDWtvY+/OytzQZvZwodXLzOy3sc3HgHHljOdKfJBQXx/N1U4pz9IuX5+OMyCoG46NPyZr69GSmmLbC8xsQTdH/jDwi9i2Ab+SZMCNWfp1Jd6PqHRkY7y//faD2trOppE0B2OaI7K2trNNPN++TvfxSMmqYlMWW3UxJL2doMTfFit+m5mtk3QI8KCkP5rZw4X6cSXeT0hzKGZ1Hmbp7803Qs7u/faHN9/M/yWRzxGZVuZ27sqQi5Tcd+2jSElwRT5QkTQNuAk4z8z2ra1gZuui/xsk3Q2cArgSrwYqHdmYbx3LoUPhc/9QeN98jkhX2j2DR0oOLiRNAO4CPmhmL8TKDwBqzGxH9PkcYH6x/lyJ9xMqHS3n0XfVg1+rgYWk24GzCLbztcA1QB2Amd0AXE1YY/jbkqBjKuGhwN1R2RDgR2b2y2LjuRLvJ2R1PPZVf07PUV+frrD9WlUnZnZxkfqPAB9JKV8JnFDqeGUrcUnDCTabYVF/d5rZNYk2lwH/RlhrE+B6M7up3LEHEpWObJw9C+6+G5KL4G1vgS9/uWOdzAnju2/r9hSzlcEjJZ1yqMSbeCswy8x2SqoDHpX0CzN7LNHux2Z2RQXGG5BUOrLx5TVdFTh0lOXWyXxyKbRHK/SU4kyttCN2MOORkk45lK3EowWPd0abddGfL4LcDSoZ2ZhvXcwkOQWeI6sz1VPMVhaPlHS6S0WyGEqqlbQM2AA8aGaPpzR7l6Tlku6UND5PP5dLapLUtHlz5qhTJ4Vy1r/M4lBzZ5zj9A8qosTNrM3MphPCR0+RlEz6ci/QaGbTgAeBW/P0s8DMZprZzFGjRldCtEFLOetfZnGoedpSx+kfVHq1+22Sfg2cC/whVr451uwm4BuVHHegkNVReP+iYC6x9g4HZXLVnXzrYnZBdDJ+1dSmO1OTYzY2whuv93CK2eXNsGRxzFA8G6Z1PSFp0Y7gNmZncFCJ2SljgD2RAt8POBv4eqLNWDN7Ndq8AHiu3HEHGlkdhfcv6qyccw5KKH35NAlQZweoUtqljblqZchFvmVLD81OWd4M992bCGOMTkhMkadFOy5cGI5pn8PWIyCdAUwl3sTHArdKqiWYZ35iZvdJmg80mdlC4JOSLgD2AluAyyow7oAiq6Mwn8Ny6dLOSjyLY9OMLi7otrbsY65eDVd/sfg43WLJ4jxhjIs7KfG0aMd8qXM9AtIZiFRidspy4MSU8qtjn68Crip3rIFMVkdhPodlsrySjs2sY1aUjCfEU+c6gx1fY7OfkNVRmM9hmSyvpGMz65gVJeMJKcWR6k5XZyDiSryfMHtWcAzGSXMUzpiRvn+yPF+7ODU1wZFZqTEryqzZ6Sdk1uzOzVLOW20tKGEnqq1p73YE5P33w/z5IdJ1/vyw7Tj9BVfi/YSpU2HePKhvABT+z5vX1VE4dw7MPLnjLVg1YTvp1JwwPijpgghOOqlyY1aUaVPh/Hkdr8/19WE7MTtl2jQ4//zOzU6csAnR2TBu7e2w5uWSxbj/fmhqikW6Wth2Re70FzwBVj8ia8Tm3DnFFejiJV2jMZO0t8GLL8KVn6rMmBVn2tTUKYVdmiWiHa+dP5T2xK3dzhCWLG1g2tzSRCjoSC6xL8fpCfxNfICS1Yk3EJ19LTaipPJCpOWfKVTuOL2NK/EBSlYn3kB09tVrR0nlhVDaxPkC5Y7T27gSH6CkOUqTDNR1MmfN2EYduzuV1bGbWTMyrgAdo0+cuo5TAm4T7wWyhtOntXt5TfEQ+zSmTu26b2NjD0ZYlsuG1+CxP8Bv7y0YYp+FaXMnAC+zZGkDLTaCeu3Yp8Cvnb+9c9n4CQXD83N276VLgwlFCtdg/Hi49toMof50P3WAByY5WXAl3sNkDadPa/fzn3d2TpYSYt/cDE8v6wjIsXZYuyZ99kmf09wML2yB1l1hO0+IfSlMmzsh5sQcyfL7t3Ff02HsYWgYwkaysGl/bGk77VbTMWxKeP7cuZ2dmJlD/Re2g61kWntLwePyhZKdcnBzSg9TKJy+WLt8s0uyhNRnHbdfsHhx11j5XIh9hViytGGfAs/RxpB9CrzTsEXOUb5Q/y652dtqWNJ+ZqKw63EVWijZqT4k3Sxpg6Q/5KmXpG9KWhGl5z4pVneppBejv0uzjOdKvIfJGk5fyiyRLOHuVZXvuxeELWVmSrFhSwr1J8VznDF1QL+8Vk4WbiFkcs3HecCU6O9y4DsAkg4mLKp8KnAKcI2kg4oN5kq8h8kaTl/KLJEs4e5Vle+7F4QtZWZKsWFLCvWn+ArIVXWtnKKY2cOERH/5uBC4zQKPAQ2SxgLvJCyqs8XMthLWXij0ZQC4TbzHyboA8uxZcM89na0KUvp85MZGuPa6zg7Kp5aF9LA5xowJ4/Rovu9KMXs2PJe451NC7Mth1oxt3NO0f6cgINGGpE4mlbQFipff9hRLVk2ihXrqaWHKmG08XdfY6dzW1kJbmxFP5ivamVXzG4j/csqTOsAXSu5b2tqNHW/sKd4wMFpSU2x7gZktKGG4I4A1se21UVm+8oK4Eu9hSlkAOamvVQONE0PK1/gMk7VrOjtA77q7684bNwZFvntPP52NEmfqVGh7FdZGJsQyZ6ekMn4CerK9k0KtqREnnlTDiy/mnxWy/LanuG/VWzocojTw9Mb9OWHMal7c3bhvv4OHbmfVxs4mG0OsGXUC03avKjg7xRdK7gfU1tE2YmzW1pvMbGZPilMKrsR7gSzh9IuXhDD4OO1tYUpgPGf3tdd1dYLlW5Z640a45pqSxe0bDjkUTjsQzj6peNtusGQJtLV3tkO1tQcFfuWVBfZbNamLQ3QPQ3lxYwNXxs7t/C8fSNclNcTSjROYe02BASJ8oeRBxTogvs7wuKhsHXBWovyhYp25Tbyf0BMOUKeD7joPUx2TKeWWuiZS/nJnULMQuCSapXIa0BKtfPYAcI6kgyKH5jlRWUH8TbyfUF8fTCNp5VnaOYWpr09X2EWdmLTQQkNqObFyYakKO6TEdUU+mJB0O+GNerSktYQZJ3UAZnYDsAiYA6wA3gA+FNVtkfQVILcY4nwzK+QgBSqzxuZw4GFgWNTfnWZ2TaLNMOA2YAawGXivma0ud+z+SjLycsoUOtld02zTpThA7/55tmmGY8Z0dYB22ybe3NzVaJvW2aJoReX29pALNxfeWGTfDS+28Nhv23h6KQWjKVnTNRIzRGgWJp/zcJXSiG8AABkzSURBVOjQkCc8x6RJcMn0jijLWbXTuadtbieHaA17mTJmG9de27BPtsYxOyObeFxhGzMmbQFGd5Klpxd29oWj+xYzu7hIvQGfyFN3M3BzKeNV4k28FZhlZjsl1QGPSvpFNHUmx4eBrWY2WdJFhIWU31uBsfsdaZGX8UWG80VsZnWAvrwm+7JoW7Z0zHbJN24mkge1LU9nixbBE7GDbW8P2zmlnmff5vtf5vntY2llF7A7bzTlwnvasfbD9ynUFhvJfU3DgZeLKvI05+HQocFvEGfVKuO2VQdwSW5qYNve1PfrJzdNoD3yRbS0wBt1I5k0ZjurNx6IIRQp8LmXdFXgyejMe+4JM5H2XasyIjZ7un+n/1GJNTYN2Blt1kV/SVfbhcCXos93AtdLUrRvOm17qN3xarni9Tq7XoMzMjwYu16D2sbOZdMbYfr/STRMTG+uN5hdRvKltHGL8tqfYNqR6eWNMSVlO2DG0dn7bBwNjGXH8JG0b++8xFAbQ7ou4txeQ9KNs4ehmfOEJ52H8TfwDsQqJu3bWsLsIEuMdmq7yLZnD2zZPZKr9/0GFck3cEiPzkyLzO3uws493b/T/6iITTxa6X4pMBn4lpk9nmiyb/6jme2V1AKMAjYl+rmcEMHE+MMOYeTmP1VCvF5l6CswNEuuacHIzaX3P2JT8TYVH/eVZ9MnrEuwOaaoNq0oQQ7Ruv5gjjgMWjmMacfuYv8RQ4vvl8oo1u8o3e58TN4vQ2M9ZwMwljGMLcGmvb5ITNHYyZB5IluG/nq7f6f/URElbmZtwHRJDcDdko43s9S8AUX6WQAsAJhx/LE25MDqC1l75gXYlmEGSUM9vKMbK+X8bin7fsZ3h26N+8LKPF7XBpgzr2N76dPZbT31DXDGGbD+RSbV7mbl5om8saOOxsY6mn6bsY+IBu3g+KNKX/DhZ7fnrzue3wOwmD9nW54ZKl3kqIfjjyrcZvGibPdH1v56u3+n/1HR2Slmtk3SrwmhonElnpsXuVbSEKCe4ODM39ebr7NreVOhJv2Svzoenn+h65zvOKqBY46GXctL7/8v3wKvvJLWKZ1+4kuhLK5Ta2rh6KO6Me7xU+CFFzqHk9bWwlFTIH6N3jI5o3A1cMwUWPki2t3KzNGt7Nk4BA45gpc2DGPmn9cg2glzO2IRkDIww2ImlRraOPqgN9n5yHOwqxWGD4NJk9jAoaxaBbt2wfDhwWF5CK8RCkO7s88+ka2vD+si7UFsZSdvB2Ae63mekcGEsu9o2kHCrEO23LndubtLd52Y9+6u90e+a3XYobD4kcQxHFK5/rPI6/R/KjE7ZQywJ1Lg+wFnExyXcRYClwK/A94NLCloDwf2vLGbTctXlSter1MDjNoVdEXrLhhSB3v3dDah1gg2tUJNkQcyjYOAzZvglVc78lsfPhZG1neMOSx64KFrWc162LS+Gwe2awisWtNJUbL+dVgfu0YbNsNzKdcsmT9Agt11cMjr2BvbGXUSvHUM1G3aiw6ZwDB2M6lhK7S0sMom0sowhtHKJFYDxiomxcpWcYg2Qp3BSIBWNmx5jVbqOXSkojJo3WLsx2scMnL7vnYTax5h+ZhT2PpGhyI/aP9Wpr3ZsTLyRNYwkVZWcSStDA2yHb4L6utZtRJaW2HYMJh0ZHEFCzDxkPCX3Bc6l40aBevXw6Ej6TiGrbDfIYXHydp/Vnmd/o+K6NLiHUjTgFuBWoIO+4mZzZc0H2gys4XRNMTvAycSEsNcZGYr83YKnNDYaIu++E9lydYfuPa6/JaILAsUVxXXXRdmn2ShoQE+9Sl2r3qe0dMm0b67Ffuz42g9bEoIfy6lrwTXcmVqkE49LVzJtalyFD2GZLseZlDdNxXgiDEHLC0nFP6E6SfZL/7n0V4Zq9JUYnbKcoJyTpZfHfu8C3hPuWNVI4MqzWhJOVo72taMHgeb1jKkrobW7vSV7Dr36pqlPGtIbC9fsH4ihlMFeNh9DzOo0oyWlKO1gvlek7uyPXt51pzAvXzB+okYThXgYfc9TNZIzAHBrFk037OSxW1n7kvbOlu/ZirNnW3iNTU0T/4rFl8Hx9fDumZ479lw+JGd++py4kK+1y7DNtecwOL2t9PCSOrZzhSt4EmdRHt7zPFYYxzcvoX5fBEjOE5n6CnmTnktmFDi4bXLlhXPC5s1grWb5LtvpkzpXhRu1nVenerDlXgPU0oq2mqnmancy1vYE/3Aa6GBe5kH1s7U2GSl5vbjuPfJI9jTDlYfTNDPvwC1R8J+o6JGuRMUV5QHHwwrV3Yd085nTy6Kk3qe1Ekk85WYdQ7iMWposhnQtJS5FkWZbtsWFPj06XTJT1toQdR8EaxlkHbfTJkS1k0ttl5rkqzrvDrViSvxXiBLKtqBwOIlYV3JOHtsCIuZ3UmJL2Y2e9o7R2i2t4WZNMcdGytMnrivfKXrmMxmjyUiKtu7Bufk898vtROZy30xgfcEBV7IiVloUcwKXujk4aelIc6tm1po2ELrrQ6G+3Kg4zZxp2LkdcYlZorkS++6qzW1uIOU+PF8TsysWNoj0N1FNnvY69jtdLruJB3QuBJ3KkZeZ1xincnUdScJ088LUtP1ds3nxMyKSIkO7a7TtYe9jt0d1p2kAxtX4oOd5ubg2Js/P/xvbu52V7NnQV1tZ6VYp73MZnHndizuojyF7QtQyitbY2PXMVlMXU1nZ2dNbVd9rzzpT2boqc4FWRa3nDUrtCt1vzKZnWfYYk7y7u7nVAeuxAczOY/Xtm3BaJxz0HVTkU+lmXncSz3bAKOebcyruZ+pk17v0Ko1Nbw8/KguZgxDtLwce0NPk23NGjjyyE59TT15OPP+spb6BkAhGOakrn7N1HUZJJgw85AQyCOF//PmFTcUT50a2pW6X5nkho0faynilrqfUx24Y3MwU2kH3ZIlTG3bxlSWdZS1AVsb4IsdC4Uu/XJ6gqtXXx/JYcVk27KlU18AU+nqAEzmrknLy2UGi18cz9TuRGL2kbe6u8MOFuf6YMTfxAczlfZ4Zewv1ZlIYj3KMmTrZuCo41QlrsQHM5X2eGXsL9WZSG49ytL66maTbrV1nKxIOlfS85JWSPp8Sv1/SVoW/b0gaVusri1Wt7DYWG5OGcykRUWW4qBLRi1OngxPP921v8mTO0VFzhjxTpp2HNOlu7EHxGaa5JGtecpfs7hIxGJatGNttBhP3MwymJx7HrHZe0SL5HyLkNF1LfCEpIVm9myujZn9faz939E5/9SbZjY963iuxAczaVGRWcPH06IWn34aTjgBVqzIr9i3bWNu7Z3AO1nKDIJxpZ3DeYWjJsQmDKbI1jzlr7l32fiikYf5omTTygaDIvOIzV7nFGBFLlOrpDsIS1Q+m6f9xcA1eeqK4kp8sNNdj1c+x+OKFV1TuybbtbUxl0XMZRG7Gct7+DXtTIZVw+HYt+aVbXEJEYv5DmswKi2P2CzOnjZYn90/MlpSfMWaBdGqZDn2LUcZsRY4Na0jSROBScCSWPHwqP+9wNfM7OeFhHEl7nSPrI7HUjyHrYVDNj3ysHv4eSuOVEPtkP2yNt9UwXziFwF3Rktc5phoZuskHQkskdRsZnkXHHbHptM9sjoeS/EcDiscsumRh93Dz1uvk1uOMse4qCyNi4BOq72a2bro/0rgIVLWa4hTieXZxgO3AYcSfEcLzOy6RJuzgHuA3Npdd5nZ/HLHdjrI5LiqZPrUWbPgnnu6rruZcGIyeTLNT+3tmp625pmu+3YJ2exMKWl9mxetYfHSBlraD6S+ZiezZ2xj6pzxXRtmoNqdgoMqHXL/4AlgiqRJBOV9EfC+ZCNJxxBWXPxdrOwg4A0za5U0Gjgd+EahwSphTtkLfMbMnpQ0Algq6cG4JzbiETM7vwLjOQkyOa56IX0qZvDUUx3Keds2mp/aG6WKDVkLW2jg3poL4cQZTF1xd2jXUA9HHVV00cesaX2bF63h3icOYw8h1rylfQT3PjEcWFOyIh8ITsHBlA65P2BmeyVdATxAWLbyZjN7Jr5kZdT0IuCOxHrDxwI3SmonWEq+lqJLO1GJ5dleBV6NPu+Q9BzBsF9wYKdyZHJc9UB0ZpcFGlKyDC5uO3OfAt83bFsNi1dMCJGSq56HWX8Bm9ZmGjaLH3bx0oZ9CnzfmNSxeGkDU+dkGqajrwHiFPSIzd7FzBYBixJlVye2v5Sy328JQciZqahNXFIjwX7zeEr1WyU9LekXkt6SZ//LJTVJatq8Y0clRRvQZHJc9VZ0ZrJZnrSzPelUa2k/sKTygn25U9Dp51RMiUs6EPgZcKWZJfODPknwuJ4A/DeQOmXGzBaY2UwzmzlqxIhKiTbgyeS46q3ozGSzPGlne9KpVl+zs6Tygn25U9Dp51REiUuqIyjwH5rZXcl6M9tuZjujz4uAusho71SA1BSwte2dHVeVTp86a1Z6vtfazqaT2bW/6SpbDzvVZs/YRh2dbSB17GH2jG159ijQl6dxdfo5lZidIuC7wHNm9p952hwGvGZmJukUwpfH5nLHdgJTaQZWspjYDBB+w1SOZJ95rZzozDRefrmrDdwMJkyArVv3jTF11pFATa861YLzsjKzU9wp6PR3KjE75XTgg0CzpFwO0i8AEwDM7Abg3cDHJe0F3gQuSnhknXLIlwJ2yequ8eiV0j5PPple/tJLRVPF9gZT54yPOTFHRH/d7Mudgk4/phKzUx4lNeV+pzbXA9eXO5aTh77wvqXMRClY7jhOj+ARmwOBvvC+pax3WbDccZwewXOn9BEVjQLMl1I2GT1Zjg08yUknQVNTenl/oJLRqY7Tj3El3gdUPAowzWmZkgK2oqGGEybA0qXBmZlDCuV9TW9EpzpOP8GVeB/QI1GASe9bWgrYciI0kyxZ0lmBQ9iuVP/lUOnoVMfpx7gBsw/oFT9kTw/Sn0MZ+7NsjlNhXIn3Ab3ih+zpQfpzKGN/ls1xKoybU/qAXkkNWu76mX3dPyWkk006MadMgWXLelQ2x+kvuBLvA3olCrDSEZq93H/mdLJpTsxly2D6dHjxRZ+d4gx4XIn3Eb0SBdjTg/Rg/5nTyeZzYr74Yue1Ph1ngOI2cadfkjmdrDsxnUGOK3GnX5I5naw7MZ1Bjitxp1+SOZ1spVPsOk6V4Urc6ZdMnTOeeSevp75mB2DU1+xg3snru85OmToV5s2DhoYQMdrQELbdien0IZLOlfS8pBWSPp9Sf5mkjZKWRX8fidVdKunF6O/SYmO5Y9Ppt2ROJ+u5Yp1+hKRa4FvA2cBa4AlJC1MWPP6xmV2R2Pdg4BpgJmCEhecXmtnWfOP5m7jjOE5lOQVYYWYrzWw3cAdwYcZ93wk8aGZbIsX9IHBuoR38TdxxnEFP6x5Y+UrmXPijJcVTeC4wswWx7SOANbHttcCpKf28S9IZwAvA35vZmjz7HlFImEoszzYeuA04lPD6v8DMrku0EXAdMAd4A7jMzPIsDeNUDE/H6jiZqJFx4NDMSnyTmc0sc8h7gdvNrFXSR4FbgW554ythTtkLfMbMjgNOAz4h6bhEm/OAKdHf5cB3KjCuU4hcJOO2bSG7YC4da3NzX0vmOAOddUDcAz8uKtuHmW02s9Zo8yZgRtZ9k5StxM3s1dxbtZntAJ6j6+v/hcBtFngMaJA0ttyxnQIUSsfqOE5P8gQwRdIkSUOBi4CF8QYJ/XcBQW8CPACcI+kgSQcB50RleamoTVxSI3Ai8HiiKp+d59XE/pcT3tQ54uCDKyna4MMjGR2nTzCzvZKuICjfWuBmM3tG0nygycwWAp+UdAHBkrEFuCzad4ukrxC+CADmm9mWQuNVTIlLOhD4GXClmW3vTh+Rc2ABwAmNjVakuVOI+vpgQkkrdxynRzGzRcCiRNnVsc9XAVfl2fdm4OasY1VkiqGkOoIC/6GZ3ZXSpGQ7j1MmHsnoOIOCspV4NPPku8BzZvafeZotBC5R4DSgxcxezdPWqQQeyeg4g4JKmFNOBz4INEtaFpV9AZgAYGY3EH5WzAFWEKYYfqgC4zrF8EhGxxnwlK3EzexRQEXaGPCJcsdyHMdxOuNh947jOFWMK3HHcZwqxpW44zhOFeNK3HEcp4pxJe44jlPFuBJ3HMepYlyJO47jVDGuxB3HcaoYV+KO4zhVjCtxx3GcKsaVuOM4ThXjStxxHKeKcSXuOI5TxbgSdxzHqWJciTuO41QYSedKel7SCkmfT6n/tKRnJS2XtFjSxFhdm6Rl0d/C5L5JKrpQsuM4zmBHUi3wLeBswqLwT0haaGbPxpo9Bcw0szckfRz4BvDeqO5NM5uedbxKLM92s6QNkv6Qp/4sSS2xb5ar09o5juMMEE4BVpjZSjPbDdwBXBhvYGa/NrM3os3HCOsOd4tKvInfAlwP3FagzSNmdn4FxnIcx6k4b+5q5w8vvFmp7o4A1sS21wKnFmj/YeAXse3hkpqAvcDXzOznhQarxPJsD0tqLLcfx3GcvqKu1jhsxO6szUdHSjbHAjNb0J1xJX0AmAmcGSueaGbrJB0JLJHUbGZ/ytdHb9nE3yrpaeAV4LNm9kwvjes4jlNpNpnZzAL164Dxse1xUVknJL0D+CfgTDNrzZWb2bro/0pJDwEnAnmVeG/MTnmS8M1yAvDfQN6fBpIul9QkqWnzjh29IJrjOE7FeQKYImmSpKHARUCnWSaSTgRuBC4wsw2x8oMkDYs+jwZOB+IO0S70uBI3s+1mtjP6vAioi4RLa7vAzGaa2cxRI0b0tGiO4zgVx8z2AlcADwDPAT8xs2ckzZd0QdTs34ADgZ8mphIeCzRFlotfE2ziBZV4j5tTJB0GvGZmJukUwhfH5p4e13Ecp6+IXlgXJcqujn1+R579fgtMLWWsspW4pNuBswjG/rXANUBdJNANwLuBj0vaC7wJXGRmVu64juM4TmVmp1xcpP56whREx3Ecp8J42L3jOE4V40rccRyninEl7jiOU8W4Enccx6liXIk7juNUMa7EHcdxqhhX4o7jOFWMK3HHcZwqxpW44zhOFeNK3HEcp4pxJe44jlPFuBJ3HMepYlyJO47jVDGuxB3HcaoYV+KO4zhVjCtxx3GcKsaVuOM4ThVTthKXdLOkDZL+kKdekr4paYWk5ZJOKndMx3Gc/oykcyU9H+m9z6fUD5P046j+cUmNsbqrovLnJb2z2FiVeBO/BTi3QP15wJTo73LgOxUY03Ecp18iqRb4FkH3HQdcLOm4RLMPA1vNbDLwX8DXo32PAy4C3kLQq9+O+stL2UrczB4GthRociFwmwUeAxokjS13XMdxnH7KKcAKM1tpZruBOwh6MM6FwK3R5zuB2ZIUld9hZq1mtgpYEfWXl7IXSs7AEcCa2PbaqOzVZENJlxPe1gF2jvvI5c/3vHhFGQ1s6mshyqR/H8NX/6NYi/4tfzb8GHqWieXsvGrF0w+8b96Y0RmbD5fUFNteYGYLYttpOu/URB/72pjZXkktwKio/LHEvkcUEqY3lHhmohOxoGjDXkRSk5nN7Gs5yqHaj6Ha5Qc/hv6OmRUyCfdremN2yjpgfGx7XFTmOI4zEMmi8/a1kTQEqAc2Z9y3E72hxBcCl0SzVE4DWsysiynFcRxngPAEMEXSJElDCY7KhYk2C4FLo8/vBpaYmUXlF0WzVyYRJoT8vtBgZZtTJN0OnAWMlrQWuAaoAzCzG4BFwByCgf4N4EPljtnL9CvzTjep9mOodvnBj2HQENm4rwAeAGqBm83sGUnzgSYzWwh8F/i+pBWEiSEXRfs+I+knwLPAXuATZtZWaDwF5e84juNUIx6x6TiOU8W4Enccx6liXIkXQFKtpKck3dfXsnQHSaslNUtalpjXWjVIapB0p6Q/SnpO0lv7WqZSkHR0dP5zf9slXdnXcpWCpL+X9IykP0i6XdLwvpbJ6cBt4gWQ9GlgJjDSzM7va3lKRdJqYKaZ9dcAjaJIuhV4xMxuijz9+5vZtr6WqztE4dPrgFPN7KW+licLko4AHgWOM7M3I6fbIjO7pW8lc3L4m3geJI0D5gI39bUsgxVJ9cAZBE8+Zra7WhV4xGzgT9WiwGMMAfaL5jPvD7zSx/I4MVyJ5+da4HNAe18LUgYG/ErS0iilQbUxCdgIfC8ya90k6YC+FqoMLgJu72shSsHM1gH/DrxMSJXRYma/6lupnDiuxFOQdD6wwcyW9rUsZfI2MzuJkE3tE5LO6GuBSmQIcBLwHTM7EXgd6JLWsxqITEEXAD/ta1lKQdJBhKRMk4DDgQMkfaBvpXLiuBJP53TggsimfAcwS9IP+lak0oneojCzDcDdFMmG1g9ZC6w1s8ej7TsJSr0aOQ940sxe62tBSuQdwCoz22hme4C7gD/vY5mcGK7EUzCzq8xsnJk1En4CLzGzqnr7kHSApBG5z8A5QOrCHf0VM1sPrJF0dFQ0mxDJVo1cTJWZUiJeBk6TtH+UKnU28Fwfy+TE6FdZDJ2Kcihwd3juGAL8yMx+2bcidYu/A34YmSNWUn1pG3JfomcDH+1rWUrFzB6XdCfwJCEM/Ck8/L5f4VMMHcdxqhg3pziO41QxrsQdx3GqGFfijuM4VYwrccdxnCrGlbjjOE4V40rccRyninEl7jiOU8X8f0sXouqR6oIWAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAEICAYAAABGaK+TAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAS7ElEQVR4nO3df5Bl5V3n8fdHBkJEEiC04xQ/MrCyQXBDiL0IkVISEoNk47C1biTGrcEdazb+qlhRs5hY1kRjmZRVRi2ttWYhm9GNCSwaIdmoGQeoVIwQm4SfIQhMQGAHpgNMgFhiSL77x30Gbnp6pm//uLfnCe9X1a17znOe55wvz718+sw5fW+nqpAk9efbVrsASdLSGOCS1CkDXJI6ZYBLUqcMcEnqlAEuSZ0ywKVlSPKWJJ9c4tg7kpy3wiXpecQA16pIcl+S167Afi5J8ulJjE2yPkklWbO3rao+VFU/PMLYDyZ5z3BbVZ1eVdcvqmhpiAEuSZ0ywDVxSf4UOBH4WJKnkryjtZ+d5DNJ9iS5ZfjyQjtb3pnkySRfapcuvgf4Y+Cctp89+zneyGOTvCHJ55M8keSBJFuGdvWp9rynjTln+Cw+A+9PsruNvy3J9ybZDLwFeEcb97HW/9l/hSQ5JMk7k9zb6rwpyQkrNef6FlVVPnxM/AHcB7x2aP044FHgQgYnFq9r61PAEcATwMta33XA6W35EuDTBzjOosYC5wH/rtXwcuAR4KK2bT1QwJqh/s/uA3g9cBNwFBDge4B1bdsHgffsbw6AXwFuA17Wxp4BvGS1XycfB/fDM3AdLH4S+ERVfaKqvlFV24EZBoEO8A3ge5O8sKp2VdUdi9j3yGOr6vqquq3VcCvwYeCHRjzO14AjgVOBVNWdVbVrxLE/DfxaVd1VA7dU1aMjjtXzlAGug8VLgf/cLp/saZc0zmVwBvtV4MeBtwK7kvzfJKeOstPFjk3y/UmuSzKb5Ctt3LEjHuta4A+BPwJ2J9ma5EWjjAVOAO4dsa8EGOBaPXO/BvMB4E+r6qihxxFV9V6Aqvqbqnodg0sgXwT+5372s++BFjf2z4BrgBOq6sUMrpNnEcf6g6r6PuA04N8yuDQyytgHgH+z0P6lYQa4VssjwMlD6/8beGOS17cbeocnOS/J8UnWJtmQ5AjgaeApBpdF9u7n+CSHzXeQJYw9Enisqv4lyVnATwxtm21jh+sePta/b2fwhwJfBf5lzrHmHddcBvxmklPazdCXJ3nJAfpLBrhWzW8Dv9Yul/xyVT0AbADeySAoH2Bw9vpt7fF24P8BjzG4Jv0zbT/XAncADyf58jzHWezYnwV+I8mTwK8DV+7dUVX9M/BbwN+1us+ec6wXMTi7fxy4n8FN2N9p2y4HTmvj/nKeOn+3HeuTDG66Xg68cJ5+0rNS5R90kKQeeQYuSZ0ywCWpUwa4JHXKAJekTq1ZuMvKOfbYY2v9+vWTPKQkde+mm276clVNzW2faICvX7+emZmZSR5SkrqX5P752r2EIkmdMsAlqVMGuCR1ygCXpE4Z4JLUKQNckjplgEtSpwxwSeqUAS5JnZroJzGlg9WWLc/PY6tvnoFLUqcMcEnqlAEuSZ0ywCWpUwa4JHVqpABPclSSq5J8McmdSc5JckyS7Unubs9Hj7tYSdJzRj0D/33gr6vqVOAM4E7gUmBHVZ0C7GjrkqQJWTDAk7wY+EHgcoCq+teq2gNsALa1btuAi8ZVpCRpX6OcgZ8EzAL/K8nnk1yW5AhgbVXtan0eBtaOq0hJ0r5GCfA1wCuB/1FVZwJfZc7lkqoqoOYbnGRzkpkkM7Ozs8utV5LUjBLgDwIPVtWNbf0qBoH+SJJ1AO1593yDq2prVU1X1fTU1D5/VFmStEQLBnhVPQw8kORlrel84AvANcDG1rYRuHosFUqS5jXql1n9AvChJIcBO4GfYhD+VybZBNwPvGk8JUqS5jNSgFfVzcD0PJvOX9lyJEmj8pOYktQpA1ySOmWAS1KnDHBJ6pQBLkmdMsAlqVMGuCR1ygCXpE4Z4JLUKQNckjplgEtSpwxwSeqUAS5JnTLAJalTBrgkdcoAl6ROGeCS1CkDXJI6ZYBLUqcMcEnqlAEuSZ0ywCWpUwa4JHXKAJekTq0ZpVOS+4Anga8Dz1TVdJJjgCuA9cB9wJuq6vHxlClJmmsxZ+CvrqpXVNV0W78U2FFVpwA72rokaUKWcwllA7CtLW8DLlp+OZKkUY0a4AV8MslNSTa3trVVtastPwysnW9gks1JZpLMzM7OLrNcSdJeI10DB86tqoeSfCewPckXhzdWVSWp+QZW1VZgK8D09PS8fSRJizfSGXhVPdSedwMfBc4CHkmyDqA97x5XkZKkfS0Y4EmOSHLk3mXgh4HbgWuAja3bRuDqcRUpSdrXKJdQ1gIfTbK3/59V1V8n+QfgyiSbgPuBN42vTEnSXAsGeFXtBM6Yp/1R4PxxFCVJWpifxJSkThngktQpA1ySOmWAS1KnDHBJ6pQBLkmdMsAlqVMGuCR1ygCXpE4Z4JLUKQNckjplgEtSpwxwSeqUAS5JnTLAJalTBrgkdcoAl6ROGeCS1CkDXJI6ZYBLUqcMcEnqlAEuSZ0ywCWpUwa4JHVq5ABPckiSzyf5eFs/KcmNSe5JckWSw8ZXpiRprsWcgb8NuHNo/X3A+6vqu4HHgU0rWZgk6cBGCvAkxwNvAC5r6wFeA1zVumwDLhpHgZKk+Y16Bv57wDuAb7T1lwB7quqZtv4gcNx8A5NsTjKTZGZ2dnZZxUqSnrNggCf5D8DuqrppKQeoqq1VNV1V01NTU0vZhSRpHmtG6PMDwI8muRA4HHgR8PvAUUnWtLPw44GHxlemJGmuBc/Aq+pXq+r4qloPXAxcW1VvAa4Dfqx12whcPbYqJUn7WM7vgf934O1J7mFwTfzylSlJkjSKUS6hPKuqrgeub8s7gbNWviRJ0ij8JKYkdcoAl6ROGeCS1CkDXJI6ZYBLUqcMcEnqlAEuSZ0ywCWpUwa4JHXKAJekThngktQpA1ySOmWAS1KnDHBJ6pQBLkmdMsAlqVMGuCR1ygCXpE4Z4JLUKQNckjplgEtSpwxwSeqUAS5JnTLAJalTCwZ4ksOTfDbJLUnuSPLu1n5SkhuT3JPkiiSHjb9cSdJeo5yBPw28pqrOAF4BXJDkbOB9wPur6ruBx4FN4ytTkjTXggFeA0+11UPbo4DXAFe19m3ARWOpUJI0r5GugSc5JMnNwG5gO3AvsKeqnmldHgSO28/YzUlmkszMzs6uRM2SJEYM8Kr6elW9AjgeOAs4ddQDVNXWqpququmpqakllilJmmtRv4VSVXuA64BzgKOSrGmbjgceWuHaJEkHsGahDkmmgK9V1Z4kLwRex+AG5nXAjwEfATYCV4+z0C1bxrn3g++4krSQBQMcWAdsS3IIgzP2K6vq40m+AHwkyXuAzwOXj7FOSdIcCwZ4Vd0KnDlP+04G18MlSavAT2JKUqcMcEnqlAEuSZ0ywCWpUwa4JHXKAJekThngktQpA1ySOmWAS1KnDHBJ6pQBLkmdMsAlqVMGuCR1ygCXpE4Z4JLUKQNckjplgEtSpwxwSeqUAS5JnTLAJalTBrgkdcoAl6ROGeCS1CkDXJI6tWCAJzkhyXVJvpDkjiRva+3HJNme5O72fPT4y5Uk7TXKGfgzwC9V1WnA2cDPJTkNuBTYUVWnADvauiRpQhYM8KraVVWfa8tPAncCxwEbgG2t2zbgonEVKUna15rFdE6yHjgTuBFYW1W72qaHgbX7GbMZ2Axw4oknLrVOSVq+LVu+pY478k3MJN8B/Dnwi1X1xPC2qiqg5htXVVurarqqpqemppZVrCTpOSMFeJJDGYT3h6rqL1rzI0nWte3rgN3jKVGSNJ9RfgslwOXAnVX1u0ObrgE2tuWNwNUrX54kaX9GuQb+A8B/AW5LcnNreyfwXuDKJJuA+4E3jadESdJ8Fgzwqvo0kP1sPn9ly5EkjcpPYkpSpxb1a4TSt6zrr1/Fg5+3isdWzzwDl6ROGeCS1CkDXJI6ZYBLUqcMcEnqlAEuSZ0ywCWpUwa4JHXKAJekThngktQpA1ySOmWAS1KnDHBJ6pQBLkmdMsAlqVMGuCR1ygCXpE4Z4JLUqX7+pNqq/cmr81bpuJJ0YJ6BS1KnDHBJ6pQBLkmdWjDAk3wgye4ktw+1HZNke5K72/PR4y1TkjTXKGfgHwQumNN2KbCjqk4BdrR1SdIELRjgVfUp4LE5zRuAbW15G3DRCtclSVrAUq+Br62qXW35YWDt/jom2ZxkJsnM7OzsEg8nSZpr2Tcxq6qAOsD2rVU1XVXTU1NTyz2cJKlZaoA/kmQdQHvevXIlSZJGsdQAvwbY2JY3AlevTDmSpFGN8muEHwb+HnhZkgeTbALeC7wuyd3Aa9u6JGmCFvwulKp68342nb/CtUiSFsFPYkpSp/r5NkJJWqYt15+3Oscd0349A5ekThngktQpA1ySOmWAS1KnDHBJ6pQBLkmdMsAlqVMGuCR1ygCXpE4Z4JLUKQNckjplgEtSpwxwSeqUAS5JnTLAJalTBrgkdcoAl6ROGeCS1CkDXJI6ZYBLUqcMcEnqlAEuSZ1aVoAnuSDJXUnuSXLpShUlSVrYkgM8ySHAHwE/ApwGvDnJaStVmCTpwJZzBn4WcE9V7ayqfwU+AmxYmbIkSQtZs4yxxwEPDK0/CHz/3E5JNgOb2+pTSe5a4vGOBb68xLFL9u4s2GVV6hqBdS3OqtW1wHvM+Vqcg7Kud2fZdb10vsblBPhIqmorsHW5+0kyU1XTK1DSirKuxbGuxbGuxXm+1bWcSygPAScMrR/f2iRJE7CcAP8H4JQkJyU5DLgYuGZlypIkLWTJl1Cq6pkkPw/8DXAI8IGqumPFKtvXsi/DjIl1LY51LY51Lc7zqq5U1Tj2K0kaMz+JKUmdMsAlqVMHRYAv9JH8JC9IckXbfmOS9UPbfrW135Xk9ROu6+1JvpDk1iQ7krx0aNvXk9zcHit6c3eEui5JMjt0/J8e2rYxyd3tsXHCdb1/qKZ/TLJnaNtY5ivJB5LsTnL7frYnyR+0mm9N8sqhbeOcq4Xqekur57Ykn0lyxtC2+1r7zUlmJlzXeUm+MvRa/frQtrF9tcYIdf3KUE23t/fTMW3bOOfrhCTXtRy4I8nb5ukzvvdYVa3qg8EN0HuBk4HDgFuA0+b0+Vngj9vyxcAVbfm01v8FwEltP4dMsK5XA9/eln9mb11t/alVnK9LgD+cZ+wxwM72fHRbPnpSdc3p/wsMbnyPe75+EHglcPt+tl8I/BUQ4GzgxnHP1Yh1vWrv8Rh8XcWNQ9vuA45dpfk6D/j4cl//la5rTt83AtdOaL7WAa9sy0cC/zjP/49je48dDGfgo3wkfwOwrS1fBZyfJK39I1X1dFV9Cbin7W8idVXVdVX1z231Bga/Cz9uy/kKg9cD26vqsap6HNgOXLBKdb0Z+PAKHXu/qupTwGMH6LIB+JMauAE4Ksk6xjtXC9ZVVZ9px4XJvbdGma/9GetXayyyrom8twCqaldVfa4tPwncyeBT6sPG9h47GAJ8vo/kz52AZ/tU1TPAV4CXjDh2nHUN28Tgp+xehyeZSXJDkotWqKbF1PWf2j/Xrkqy9wNXB8V8tUtNJwHXDjWPa74Wsr+6xzlXizX3vVXAJ5PclMFXVUzaOUluSfJXSU5vbQfFfCX5dgYh+OdDzROZrwwu7Z4J3Dhn09jeY2P/KP3zQZKfBKaBHxpqfmlVPZTkZODaJLdV1b0TKuljwIer6ukk/43Bv15eM6Fjj+Ji4Kqq+vpQ22rO10EryasZBPi5Q83ntrn6TmB7ki+2M9RJ+ByD1+qpJBcCfwmcMqFjj+KNwN9V1fDZ+tjnK8l3MPih8YtV9cRK7vtADoYz8FE+kv9snyRrgBcDj444dpx1keS1wLuAH62qp/e2V9VD7XkncD2Dn8wTqauqHh2q5TLg+0YdO866hlzMnH/ijnG+FrK/ulf9qyKSvJzB67ehqh7d2z40V7uBj7Jylw0XVFVPVNVTbfkTwKFJjuUgmK/mQO+tscxXkkMZhPeHquov5ukyvvfYOC7sL/ImwBoGF+9P4rmbH6fP6fNzfPNNzCvb8ul8803MnazcTcxR6jqTwY2bU+a0Hw28oC0fC9zNCt3QGbGudUPL/xG4oZ67afKlVt/RbfmYSdXV+p3K4KZSJjFfbZ/r2f9NuTfwzTeYPjvuuRqxrhMZ3NN51Zz2I4Ajh5Y/A1wwwbq+a+9rxyAI/6nN3Uiv/7jqattfzOA6+RGTmq/23/4nwO8doM/Y3mMrNrnLnIQLGdy9vRd4V2v7DQZntQCHA/+nvaE/C5w8NPZdbdxdwI9MuK6/BR4Bbm6Pa1r7q4Db2pv4NmDThOv6beCOdvzrgFOHxv7XNo/3AD81ybra+hbgvXPGjW2+GJyN7QK+xuAa4ybgrcBb2/Yw+MMk97ZjT09orhaq6zLg8aH31kxrP7nN0y3tNX7XhOv6+aH31g0M/YCZ7/WfVF2tzyUMfqlheNy45+tcBtfYbx16rS6c1HvMj9JLUqcOhmvgkqQlMMAlqVMGuCR1ygCXpE4Z4JLUKQNckjplgEtSp/4/KspA55hBN4UAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEWCAYAAACT7WsrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nO3deZxVdf3H8dcbEEEERRkVQQUEU0xTG0nNyjILzURcQU39adKipWVpppap5S9zL8rI3MV9wzL3peznwmiiuSEiCrgwioqKC8vn98f3jAzDvTN3YO49d2bez8fjPOacc885930vOp8553vO96uIwMzMrJAueQcwM7Pq5SJhZmZFuUiYmVlRLhJmZlaUi4SZmRXlImFmZkW5SJiZWVEuEtbhSJoh6QNJ70l6TdLFklZtss12ku6R9K6kdyTdIml4k236SDpH0svZsV7IlvtV9hOZ5cdFwjqqb0bEqsAWwJbAcQ0vSNoWuAO4GVgXGAxMAf4taUi2TXfgbmBTYCTQB9gWeBMYUa7QkrqV69hmy8NFwjq0iHgNuJ1ULBqcDlwaEedGxLsRMTciTgAeAk7KtjkQWB8YHRFPR8TiiJgTEadExK2F3kvSppLulDRX0uuSfp6tv1jSqY2220HSrEbLMyQdK+kJ4P1s/romxz5X0nnZ/GqS/irpVUmzJZ0qqesKflVmBblIWIcmaSCwMzAtW14F2A64tsDm1wA7ZfNfBW6LiPdKfJ/ewF3AbaSzk6GkM5FSjQW+AawOXAXskh2TrADsA0zMtr0YWJi9x5bA14Bvt+K9zErmImEd1U2S3gVmAnOAX2br1yD9d/9qgX1eBRraG9Yssk0xuwKvRcSZEfFhdobycCv2Py8iZkbEBxHxEvAYMDp77SvA/Ih4SNLawC7AURHxfkTMAc4GxrTivcxK5iJhHdXuEdEb2AHYmCW//N8CFgP9C+zTH3gjm3+zyDbFrAe8sFxJk5lNlieSzi4A9mPJWcQGwErAq5LelvQ28GdgrRV4b7OiXCSsQ4uI+0mXZ87Ilt8HHgT2LrD5Piy5RHQX8HVJvUp8q5nAkCKvvQ+s0mh5nUJRmyxfC+yQXS4bzZIiMRP4COgXEatnU5+I2LTEnGat4iJhncE5wE6SPpMt/ww4SNIPJfWW1DdrWN4W+FW2zWWkX8jXS9pYUhdJa0r6uaRdCrzH34D+ko6StHJ23M9lrz1OamNYQ9I6wFEtBY6IeuA+4CLgxYh4Jlv/KunOrDOzW3S7SNpQ0peW43sxa5GLhHV42S/cS4FfZMsPAF8H9iC1O7xEagDePiKez7b5iNR4/SxwJzAPeIR02WqZtoaIeJfU6P1N4DXgeeDL2cuXkW6xnUH6BX91idEnZhkmNll/INAdeJp0+ew6WndpzKxk8qBDZmZWjM8kzMysKBcJMzMrykXCzMyKcpEwM7Oi2l1nYv369YtBgwblHcPMrF159NFH34iImtbu1+6KxKBBg6irq8s7hplZuyLppeXZz5ebzMysKBcJMzMrqmxFQtKFkuZI+m+R1yXpPEnTJD0haatyZTEzs+VTzjOJi0kjehWzMzAsm8YBfypjFjMzWw5lKxIR8U9gbjObjCKNDhYR8RCwuiT3P2NmVkXybJMYwNJ96M/K1i1D0jhJdZLq6uvrW/1GV1wBgwZBly7p5xVXLE9cM7POp100XEfEhIiojYjamprW3eZ7xRUwbhy89BJEpJ/jxrlQmJmVIs8iMZs0mleDgdm6NnX88TB//tLr5s+H4496H6ZMgfdKGsLYzKxTyvNhuknAEZKuAj4HvJMNqNKmXn65yPo3esIWW6SFddaBDTeEoUOX/GyYX2ONto5kZtZulK1ISLqSNL5wP0mzSAPRrwQQEecDt5IGdJ8GzAf+pxw51l8/XWJaZn3PN2CDjWH6dHjttTT9+9/Lbti377KFo+HnOuuAVI7YZmZVod0NOlRbWxut6ZajoU2i8SWnVVaBCRNg//2BRYtg9mx44QWYNi1NDfMvvND85ahevVKxaHoWsuGGsN560LXr8n9QM7M2JOnRiKht9X4dvUhAKhTHH58uPa2/Pvz611mBaEkEzJmzbOFoKCZzm7nDd6WVYMiQwpexBg2C7t1b9RnMzFaEi0Qe3norFY1CZyGvNtO80qVLqlaFLmMNGZLOUMzM2pCLRLV5//3U3lGogLz8MixeXHzf/v2XvXzV8LNv38p9BjPrMJa3SLS7rsLbjV69YLPN0tTUxx/DjBmFL2NNn57OQl59Ff71r2X3XWONwndhDR0Ka63lhnQza1MuEnno3h022ihNTS1aBLNmLdv+0TA/dy488kiamlp11cJnH0OHwsCB6TKXmVkr+HJTexIBr79e+C6sadNSG0kx3bun9o5CBWTQoNTQbmYdli83dQZSejZjnXVg++2XfX3u3GULR8PP116DZ59NU1Ndu6aG9EKXsYYMSfcMm1mn5DOJzuK991J7R6HLWC+/nM5Sill33cIPEw4dCqutVrnPYGbLzXc32fL76KMlDelNz0JefBEWLCi+b79+xdtBamrckG5WJXy5yZbfyivDpz6VpqYWLoSZM4tfxnrjjTQ9/PCy+/buXbxPrAED3JBu1g74TMKWX0S6VbfQ0+jTpsE77xTfd+WVlzSkNz0L2WADN6SbtTGfSVjlSam9Yt114QtfWPq1iKUb0puehbz+OjzzTJqa6to13XFV6Cxk8GDo2bMiH8/MfCZheXn33aW7NGn8c+bM5hvSBw4s3jNvnz6V+wxm7Ygbrq3j+PDD1GBe6DLWjBmpnaSYmpriXZr06+eGdOu0fLnJOo4ePWCTTdLU1MKF6ZbdYpex6uvT9OCDy+7bp0/hAjJ0aOovyw3pZsvwmYR1HIsXp4b0Yl27z5tXfN8ePZYeG6RxEVl/fejmv6esffOZhFmXLunW2gED4EtfWvq1CHjzzeIFpL4ennoqTU1165Ya0gudhQwenAqMWQflImGdg5TaJPr1g222Wfb1efOKjw3S0OHitGmFjztwYPF2kN69CzaDtLMTeOvEfLnJrCUffJAa0gudhcyYkXruLUIsApRNS2tn/+tZO1eVl5skjQTOBboCF0TE/zZ5fQPgQqAGmAscEBGzypnJrNV69oThw9PU1IIFqSG9UAF54QX4qHCBMGsvylYkJHUFxgM7AbOAyZImRcTTjTY7A7g0Ii6R9BXgNOBb5cpk1uZWWmlJg3dTixenP4/M2rFy3vM3ApgWEdMj4mPgKmBUk22GA/dk8/cWeN2s/erSheJnEb7WZO1DOYvEAGBmo+VZ2brGpgB7ZPOjgd6S1ixjJrMqEGm67768g5i1KO+nh34CfEnSf4AvAbOBZVoBJY2TVCeprr6+vtIZzZZb4cbpIOgKe+2Vxvgwq2LlLBKzgfUaLQ/M1n0iIl6JiD0iYkvg+Gzd200PFBETIqI2ImpramrKGNms7UU0mRYG7LJLem5j1KjUj5VZlSpnkZgMDJM0WFJ3YAwwqfEGkvpJashwHOlOJ7OOrWtXmDgxdTvy3//CAQekRm6zKlS2IhERC4EjgNuBZ4BrIuIpSSdL2i3bbAfgOUlTgbWBX5crj1lVWW01mDQJ+vZNP088Me9EZgX5YTqzPN11F4wcmR7ImzgRxo7NO5F1UMv7MF3eDddmndtXvwpnn53mDzkEHn003zxmTbhImOXtiCPg299O42iMGpV6sjWrEi4SZnmTYPx42H57mD0bRo9OBcOsCrhImFWD7t3h+uvT2BUPPwzf+Y57ALSq4CJhVi3WWgtuvhlWWQUuvRTOOivvRGYuEmZVZYstUoEAOOYY+Mc/8s1jnZ6LhFm12XNPOOmk9IDdmDHw7LN5J7JOzEXCrBqdeGIqFvPmwW67wVtv5Z3IOikXCbNq1KULXHIJfOYz8PzzsO++sHBh3qmsE3KRMKtWvXqlhuyaGrjzTvjpT/NOZJ2Qi4RZNdtgA7jhhjQC3jnnwIXuA9Mqy0XCrNptvz386U9p/rvfhX//O9881qm4SJi1B4ceCkceCQsWwB57wMsv553IOgkXCbP24owzUoeAc+akPp7efz/vRNYJuEiYtRfdusHVV8PQofD443Dwwe66w8rORcKsPVljjTRIUZ8+cN11cMopeSeyDs5Fwqy92WQTuPLK1HvsL3+Z7n4yKxMXCbP2aJdd4PTT0/y3vgVTpuSbxzqsFouEpM9LulPSVEnTJb0oaXolwplZM44+OhWI+fNTQ/acOXknsg6oWwnb/BX4EfAosKi8ccysZBJMmABTp6YxKPbaK42Z3b173smsAynlctM7EfGPiJgTEW82TKUcXNJISc9JmibpZwVeX1/SvZL+I+kJSbu0+hOYdWY9esCNN8K668K//pWGQvUdT9aGSikS90r6naRtJW3VMLW0k6SuwHhgZ2A4MFbS8CabnQBcExFbAmOAP7Yyv5n175/6eOrRA/7ylzQUqlkbKeVy0+eyn7WN1gXwlRb2GwFMi4jpAJKuAkYBTzc5Tp9sfjXglRLymFlTtbWpX6f99oOjjkp3QO24Y96prANosUhExJeX89gDgJmNlmexpOA0OAm4Q9IPgF7AVwsdSNI4YBzA+uuvv5xxzDq4sWPhySfhtNNg773hkUfSg3dmK6CUu5tWk3SWpLpsOlPSam30/mOBiyNiILALcJmkZTJFxISIqI2I2pqamjZ6a7MO6NRT4ZvfTIMU7bZbGrTIbAWU0iZxIfAusE82zQMuKmG/2cB6jZYHZusaOxS4BiAiHgR6AP1KOLaZFdKlC1x+OWy6KTzzTLr8tMg3JdryK6VIbBgRv4yI6dn0K2BICftNBoZJGiypO6lhelKTbV4GdgSQtAmpSNSXHt/MltGnT+q6Y4014O9/h+OPzzuRtWOlFIkPJG3fsCDp88AHLe0UEQuBI4DbgWdIdzE9JelkSbtlmx0NHCZpCnAlcHCE798zW2FDhqS+nbp2hd/+Np1dmC0HtfQ7WdIWwCWku48EzCX9Ms+lH4Da2tqoq6vL463N2p8//hEOPxxWXhn++U8YMSLvRJYTSY9GRG3LWy6txTOJiHg8Ij4DbA5sFhFb5lUgzKyVvv/9NJrdRx/B7rvD7KbNgmbNK3oLrKQDIuJyST9ush6AiDirzNnMrC2cd15qxL7/fhg9Ov3s2TPvVNZONHcm0Sv72bvAtGqZc5lZW1lppdQ+MWgQTJ4Mhx3mrjusZEXPJCLiz9nsXRGx1MjrWeO1mbUX/fqlO5623RauuAI22wyOPTbvVNYOlHJ30+9LXGdm1WyzzZbc5XTccfC3v+Wbx9qF5toktgW2A2qatEv0AbqWO5iZlcHuu6ensk84IT1o9+CD6cE7syKaO5PoTmp76MbS7RHzgL3KH83MyuLnP4d994V3301dd7xZUs//1kk11yZxP3C/pIsj4qUKZjKzcpJSj7HPPw+PPQb77AO33ZYauM2aKKVN4gJJqzcsSOor6fYyZjKzcltlFbjpJlh7bbjnHvjxj1vexzqlUopEv4h4u2EhIt4C1ipfJDOriPXWS6Pade8Of/hDGgrVrIlSisRiSZ8M4iBpA9JgQWbW3m277ZLicPjhqesOs0ZKKRLHAw9IukzS5cA/gePKG8vMKuagg9LlpoULYc89YcaMvBNZFSml76bbgK2Aq4GrgM9GhNskzDqS00+HkSPhjTdg1Ch47728E1mVKGVkOgEjga0i4m/AKpLclaRZR9K1K1x5JXzqU/DEE3DggbB4cd6prAqUcrnpj8C2pKFGIY1SN75sicwsH6uvnrruWG211KD9q1/lnciqQClF4nMRcTjwIXxyd1P3sqYys3xstBFcfXUaBvXkk+Haa/NOZDkrpUgskNSV7I4mSTWAz0PNOqqvfx3OOCPNH3QQ/Oc/+eaxXJVSJM4DbgTWkvRr4AHgN2VNZWb5Ouoo+J//gQ8+SA3Zr7+edyLLSdFuORpExBWSHgV2JA1funtEPFP2ZGaWHwn+9Cd47jn4v/+DPfZIT2avvHLeyazCip5JSOqT/VwDmANcCUwEXs+65mixJ1hJIyU9J2mapJ8VeP1sSY9n01RJbxc6jpnlYOWV4YYbYODAVCi+9z0PVtQJNXe5aWL281GgrsnPx4DXJBW97JQVkfHAzsBwYKyk4Y23iYgfRcQWEbEFaYyKG5b3g5hZGay9drrjqWdPuOgiOPfcvBNZhRUtEhGxa/ZzcEQMafJzMLAOMLqZY48ApkXE9Ij4mPQg3qhmth9LOlsxs2qy5ZZwySVp/uij4Y478s1jFVXSw3SSDpB0Yra8vqQREbEoIjZpZtcBwMxGy7OydYXeYwNgMHBPkdfHSaqTVFdfX99SZDNra3vvDSeemB6w23dfmDo170RWIa15mG6/bLkcD9ONAa6LiEWFXoyICRFRGxG1NTU1bfzWZlaSk06C0aPh7bfTYEVvuwmxMyjnw3SzgfUaLQ/M1hUyBl9qMqtuXbrApZemsbKfew7GjoVFBf+usw6knA/TTQaGSRosqTupEExqupGkjYG+wIMlpzazfKy6amrI7tcvjWZ37LF5J7IyK9vDdBGxEDgCuB14BrgmIp6SdLKk3RptOga4KsL31pm1C4MGwfXXQ7ducOaZSxq1rUNSc7+bJXUBtgHmsuRhurvzfJiutrY26urq8np7M2swYQJ85ztpZLv77ksDGFnVkvRoRNS2dr9mn7iOiMWSxkfElsCzy53OzDqecePgySfT0KejR0NdXXrwzjqUUi433S1pz2xcCTOzJc46C77yldS306hRMH9+3omsjZVSJL4DXAt8JGmepHclzStzLjNrD1ZaCa65BjbcEB57DA45xF13dDClDF/aOyK6RET3iOiTLfepRDgzawfWXBNuvhl6905jUfzGnUR3JKWcSZiZNW/TTWHixNR77AknpKJhHYKLhJm1jV13hdNOS/MHHJAata3dc5Ews7ZzzDGw337w3nup64433sg7ka2g5saTWKO5qZIhzaydkOCCC2DrrWHGDNhrL1iwIO9UtgKaO5NoPH5EPTAVeD6bf7T80cysXerZE268Efr3h/vvhx/+MO9EtgKaG09icEQMAe4CvhkR/SJiTWBXwB3Km1lxAwbATTel0e3OPz8NhWrtUiltEttExK0NCxHxD2C78kUysw5hxIh06QnS2cS99+abx5ZLKUXiFUknSBqUTccDr5Q7mJl1AAcckBqzFy5M7RPTp+edyFqplCIxFqgh9QR7I7BWts7MrGW/+Q184xswd2664+ndd/NOZK3QbAd/ABExFziyAlnMrCPq2jU9aLfNNvDUU+ns4sYb0yBGVvVKGeN6I0kTJN0h6Z6GqRLhzKyD6NMnDVbUt2/6eeKJeSeyErV4JkHq3O984ALAYxWa2fIZOhSuvRa+/vV0CerTn05DoFpVK+V8b2FE/CkiHomIRxumsiczs45nxx3h7LPT/CGHpDEorKqVUiRukfR9Sf39xLWZrbAjjoDDDoMPP4Tdd4dXX807kTWjlMtNB2U/f9poXQBD2j6OmXV4UhrN7tln4V//SqPa3Xcf9OiRdzIroJTxJAYXmEoqEJJGSnpO0jRJPyuyzT6Snpb0lKSJrf0AZtYOde8O118PG2wADz+chkL1YEVVqcUzCUkHFlofEZe2sF9XYDywEzALmCxpUkQ83WibYcBxwOcj4i1Ja7UmvJm1YzU1adyJ7baDyy6DzTeHn/wk71TWRCltEls3mr4AnATsVsJ+I4BpETE9Ij4GrgJGNdnmMGB8RLwFEBFzSsxtZh3BZz6TCgSkJ7NvvbX57a3iSrnc9ING02HAVsCqJRx7ADCz0fKsbF1jGwEbSfq3pIckjSx0IEnjJNVJqquvry/hrc2s3dhjD/jVr9LlprFjU1uFVY3leeTxfWBwG71/N2AYsAOpq4+/SFq96UYRMSEiaiOitqampo3e2syqxgknpL6d5s1LXXe89VbeiSxTSpvELaS7mQC6ApsA15Rw7NnAeo2WB2brGpsFPBwRC4AXJU0lFY3JJRzfzDqKLl3g4oth2jR4/HHYd9906albKTdgWjmV8i9wRqP5hcBLETGrhP0mA8MkDSYVhzHAfk22uYl0BnGRpH6ky0/uJtKsM+rVKzVk19bCnXemRuxzzsk7VadXSpvE/cCzQG+gL/BxKQeOiIXAEcDtwDPANRHxlKSTJTU0fN8OvCnpaeBe4KcR8WbrP4aZdQjrr586/1tpJTj3XPjrX/NO1OkpWrg3WdI+wO+A+wCR7nD6aURcV/Z0BdTW1kadH+U369guvBAOPTQVi3vuge23zztRuyfp0Yiobe1+pTRcHw9sHREHRcSBpFtb3YWjmZXPIYfAkUfCggXp7qeXX847UadVSpHo0uT5hTdL3M/MbPmdcQbstBPU18OoUfD++3kn6pRK+WV/m6TbJR0s6WDg74CfeDGz8urWDa6+GoYNS3c8HXwwLF6cd6pOp9kiIUnAecCfgc2zaUJEHFuBbGbW2TUMUtSnD1x3HZx6at6JOp1mb4GNiJB0a0RsBtxQoUxmZktsvDFcdVUaJ/uXv4RNN4U998w7VadRyuWmxyRtXfYkZmbF7LwznH56mj/wQJgyJd88nUgpReJzwIOSXpD0hKQnJT1R7mBmZks5+uhUIObPT113zHF/oJVQyhPXXy97CjOzlkjw5z/D1Knw0EPpktPdd6exKaxsSjmTeLfA9Eo5Q5mZFdSjB9xwAwwYAA88AIcf7sGKyqykNgmgHpgKPJ/Nz5D0mKTPljOcmdky+veHm25KBeOCC9JQqFY2pRSJO4FdIqJfRKwJ7Az8Dfg+8MdyhjMzK6i2Fi66KM3/6Edw11355unASikS20TE7Q0LEXEHsG1EPASsXLZkZmbNGTMGfv5zWLQI9tkHnn8+70QdUilF4lVJx0raIJuOAeZkY1j78Uczy88ppywZpGjUKHjnnbwTdTilFIn9SAMG3QTcSBpIaCxpAKJ9yhfNzKwFXbrA5ZenB+yeeQb23z+dWVibKaVI9M7Gt94yIraKiB8AgyLi44iYVu6AZmbN6t07dd2x5prw97+nS1DWZkopEtdLGtCwIOmLwIXli2Rm1kpDhqS+nbp1S09mX3553ok6jFKKxHeAmyStI2kX4PfALuWNZWbWSjvsAOedl+a//W14+OFc43QUpQxfOhn4IXAHcBLw1YiYWeZcZmat973vwXe/Cx99BKNHw+zZeSdq94p2yyHpFqDxo4yrAO8Af5VEROxWeE8zsxyddx48+yzcdx/svjv885/Qs2feqdqt5vpuOmNFDy5pJHAu6U6oCyLif5u8fjBp/OyGcv+HiLhgRd/XzDqxlVaCa6+FESOgri5derr88tT3k7Va0SIREfcDSBoMvBoRH2bLPYG1Wzpw9hzFeGAnYBYwWdKkiHi6yaZXR8QRy5nfzGxZ/fqlO5623RYmToTNN4djPVba8iil4fpaln5oblG2riUjgGkRMT0iPgauAka1PqKZ2XL49KeXnEEcdxzcckveidqlUopEt+yXPADZfCl98w4AGjdwz8rWNbVnNk7FdZLWK3QgSeMk1Umqq6+vL+GtzcxIT2GfemrqKXa//eCpp/JO1O6UUiTqJX3SSC1pFPBGG73/LaQH8zYndSR4SaGNImJCRNRGRG1NTU0bvbWZdQrHHQf77gvvvZe68HjzzbwTtSulFInvAj+X9LKkmcCxpGcnWjKb1IVHg4EsaaAGICLejIiPssULAHc9bmZtS4ILL4TPfhamT4e994YFC/JO1W6U8pzECxGxDTAc2CQitiuxO47JwDBJgyV1B8YAkxpvIKl/o8XdgGdKj25mVqJVVkljUKy9Ntx7b+pe3EpSyvClSPoGsCnQQ9ltZBFxcnP7RMRCSUcAt5Nugb0wIp6SdDJQFxGTgB9ml7IWAnOBg5f3g5iZNWvgwFQovvQlGD8eNtsMvlPKRZHOTdHC0H+Szic9SPdl0iWhvYBHIuLQ8sdbVm1tbdTV1eXx1mbWEVxyCRx8cOrn6e674YtfzDtRRUh6NCJqW7tfKW0S20XEgcBbEfErYFtgo9a+kZlZVTjoIDj6aFi4EPbcE2bMyDtRVSulSHyQ/ZwvaV1gAdC/me3NzKrbb38LI0fCG2+kO57eey/vRFWrlCLxN0mrk7rPeAyYAUwsZygzs7Lq2hWuvBI+9Sl48kn41rdgsQfaLKSUu5tOiYi3I+J6YANg44j4RfmjmZmV0eqrp647Vl89NWifdFLeiapSi0VCUg9JP5Z0A+kM4hBJPcofzcyszDbaCK6+Og2DesopcM01eSeqOqVcbrqUdPvr74E/kJ6XuKycoczMKuZrX4Mzz0zzBx8Mjz2Wa5xqU8pzEp+OiOGNlu+V1LQnVzOz9uvII1PbxIUXpv6e6urSg3dW0pnEY5K2aViQ9DnADyqYWcchwR//CNttB7NmpVHtPvqo5f06gaJFQtKTkp4g9af0f5JmSHoReBBo9QMZZmZVbeWV4YYbYL314MEH0zCoLTxs3Bk0d7lp14qlMDOrBmuvDTffDJ//PFx8cRqsqJP381T0TCIiXmpuqmRIM7OK2XLL1HUHwE9+Arffnm+enJXSJmFm1rnsvTf84hfpAbt994WpU/NOlBsXCTOzQn75y9SA/c478M1vwttv550oFy4SZmaFdOkCl16a2iWmToUxY2DRorxTVZyLhJlZMauumhqy+/VLbRPHHJN3oopzkTAza86gQXD99Wn8ibPOSnc9dSIuEmZmLfniF9PDdpBGs3vwwXzzVJCLhJlZKQ47DI44Aj7+ODVoz5yZd6KKcJEwMyvV2WfDjjvC66/D7rvD/Pl5Jyq7shYJSSMlPSdpmqSfNbPdnpJCkrv7MLPq1a1b6k58ww1Tb7GHHNLhu+4oW5GQ1BUYD+xM6l58rKThBbbrDRwJPFyuLGZmbWaNNdJgRb17p7EofvObvBOVVTnPJEYA0yJiekR8DFwFjCqw3SnAb4EPy5jFzKztDB8OEyem3mNPOCGNbNdBlbNIDAAat+zMytZ9QtJWwHoR8ffmDiRpnKQ6SXX19fVtn9TMrLV23RVOOy3NH3BAGo+iA8qt4VpSF+As4OiWto2ICRFRGxG1NTU15Q9nZlaKY46B/feH96g0TusAAAwuSURBVN+H3XaDN97IO1GbK2eRmA2s12h5YLauQW/g08B9kmYA2wCT3HhtZu2GBH/5C2y9NcyYAXvtlW6R7UDKWSQmA8MkDZbUHRgDTGp4MSLeiYh+ETEoIgYBDwG7RYRHvTOz9qNnz9Qm0b8/3H9/Ggq1AylbkYiIhcARwO3AM8A1EfGUpJMl7Vau9zUzq7h1102FYuWV4fzzlzyd3QEo2tk9vrW1tVFX55MNM6tCV1yRGrG7doU77oCvfCXvRJ+Q9GhEtPpyvp+4NjNrK/vvD8cem7oU33tveOGFvBOtMBcJM7O29Otfwze+AXPnpjue5s3LO9EKcZEwM2tLXbumB+022QSefjpdflq8OO9Uy81FwsysrfXpk7ru6NsXbrklPZXdTrlImJmVw9ChcO216czitNPgyivzTrRcXCTMzMplxx3hnHPS/CGHwOTJ+eZZDi4SZmbldPjhacCiDz9MY1C8+mreiVrFRcLMrJwk+MMf4AtfgFdeSaPafdh+Or12kTAzK7fu3eH662GDDeDhh2HcuHYzWJGLhJlZJdTUpDueevWCyy6DM8/MO1FJXCTMzCpl881TgYDUzfitt+abpwQuEmZmlTR6NJx8crrcNHYsPPNM3oma5SJhZlZpJ5yQ+naaNy913TF3bt6JinKRMDOrNAkuugi22AKmTYN994WFC/NOVZCLhJlZHnr1gptvhrXWgrvugp/8JO9EBblImJnlZf314YYbYKWV4Nxz4a9/zTvRMlwkzMzy9PnPp9HsAL73PXjggXzzNOEiYWaWt0MOgaOOggULYI894KWX8k70CRcJM7Nq8LvfwU47QX09jBoF77+fdyKgzEVC0khJz0maJulnBV7/rqQnJT0u6QFJw8uZx8ysanXrBldfDcOGwZQpcNBBVTFYUdmKhKSuwHhgZ2A4MLZAEZgYEZtFxBbA6cBZ5cpjZlb1+vZNXXestlrq6+mUU/JOVNYziRHAtIiYHhEfA1cBoxpvEBGNB3/tBbSPHq/MzMpl443TAEVdusBJJ6VikaNyFokBwMxGy7OydUuRdLikF0hnEj8sdCBJ4yTVSaqrr68vS1gzs6qx885w+ulp/sAD4fHHc4uSe8N1RIyPiA2BY4GCA8FGxISIqI2I2pqamsoGNDPLw49/nNol5s9HW26OFEh8MlVKOYvEbGC9RssDs3XFXAXsXsY8ZmbthwTnn49YBCibln65EspZJCYDwyQNltQdGANMaryBpGGNFr8BPF/GPGZm7UuPHhQqEJXUrVwHjoiFko4Abge6AhdGxFOSTgbqImIScISkrwILgLeAg8qVx8ysfcqvQEAZiwRARNwK3Npk3S8azR9Zzvc3M7MVk3vDtZmZVS8XCTOzKhZFnh4rtr6tlfVyk5mZrbhKFYRCfCZhZmZFuUiYmVlRLhJmZlaUi4SZmRXlImFmZkUp8mw2Xw6S6oFSx/brB7xRxjhtodozVns+qP6M1Z4Pqj+j8624T0VE79bu1O5ugY2IkruBlVQXEbXlzLOiqj1jteeD6s9Y7fmg+jM634qTVLc8+/lyk5mZFeUiYWZmRXX0IjEh7wAlqPaM1Z4Pqj9jteeD6s/ofCtuuTK2u4ZrMzOrnI5+JmFmZivARcLMzIrqEEVC0khJz0maJulnBV5fWdLV2esPSxpUZfm+KOkxSQsl7VXJbK3I+GNJT0t6QtLdkjaosnzflfSkpMclPSBpeCXzlZKx0XZ7SgpJFb1lsoTv8GBJ9dl3+Likb1cyXykZs232yf5bfErSxGrKJ+nsRt/fVElvVzJfiRnXl3SvpP9k/z/v0uwBI6JdT6ShUV8AhgDdgSnA8CbbfB84P5sfA1xdZfkGAZsDlwJ7Vel3+GVglWz+e1X4HfZpNL8bcFu1fYfZdr2BfwIPAbXVlA84GPhDpf/7a2XGYcB/gL7Z8lrVlK/J9j8gDdtcbd/hBOB72fxwYEZzx+wIZxIjgGkRMT0iPgauAkY12WYUcEk2fx2wo6RKDRzbYr6ImBERTwCLK5SpqVIy3hsR87PFh4CBVZZvXqPFXkCl78go5b9DgFOA3wIfVjIcpefLUykZDwPGR8RbABExp8ryNTYWuLIiyZYoJWMAfbL51YBXmjtgRygSA4CZjZZnZesKbhMRC4F3gDUrkq60fHlrbcZDgX+UNdHSSson6XBJLwCnAz+sULYGLWaUtBWwXkT8vZLBMqX+G++ZXYK4TtJ6lYn2iVIybgRsJOnfkh6SNLJi6Vrx/0l2OXYwcE8FcjVWSsaTgAMkzQJuJZ3xFNURioRVkKQDgFrgd3lnaSoixkfEhsCxwAl552lMUhfgLODovLM04xZgUERsDtzJkrPvatKNdMlpB9Jf6n+RtHquiQobA1wXEYvyDlLAWODiiBgI7AJclv33WVBHKBKzgcZ/8QzM1hXcRlI30inWmxVJV1q+vJWUUdJXgeOB3SLiowplg9Z/h1cBu5c10bJaytgb+DRwn6QZwDbApAo2Xrf4HUbEm43+XS8APluhbA1K+XeeBUyKiAUR8SIwlVQ0qiVfgzFU/lITlJbxUOAagIh4EOhB6qCwsEo2qpSpoaYbMJ10atfQULNpk20OZ+mG62uqKV+jbS8mn4brUr7DLUkNYsOqNN+wRvPfBOqqLWOT7e+jsg3XpXyH/RvNjwYeqrbvEBgJXJLN9yNdWlmzWvJl220MzCB7WLkKv8N/AAdn85uQ2iSKZq3oByjjF7ML6S+KF4Djs3Unk/7ihVQprwWmAY8AQ6os39akv5DeJ53hPFWF3+FdwOvA49k0qcrynQs8lWW7t7lf0HllbLJtRYtEid/hadl3OCX7Djeutu8QEOmy3dPAk8CYasqXLZ8E/G+lv7tWfIfDgX9n/86PA19r7njulsPMzIrqCG0SZmZWJi4SZmZWlIuEmZkV5SJhZmZFuUiYmVlRLhLW7km6oBy9vkp6r5Xb7y3pGUn3ZstXZl1c/EjSydnDiMX2rZV03opmNmtrvgXWrAhJ70XEqq3Y/jbg1Ih4QNI6wAMRMbR8Cc3Kz2cS1m5I6iXp75KmSPqvpH2z9fc1dG8h6dCsH/9HJP1F0h+y9RdLOk/S/0ma3jBuh6RVs/ExHsvGo2ixZ1RJB2THf1zSnyV1lfQLYHvgr5J+B9wBDMi2+UL2/g3vuXWWY0p2nN6SdpD0t0af88Lstf80ZMrGe7hB0m2Snpd0eqNMI7PPMCX7PF2ybWqy17tk4wvUtN2/iHUG3fIOYNYKI4FXIuIbAJJWa/yipHWBE4GtgHdJPXBOabRJf9Iv8o2BSaRu4z8ERkfEPEn9gIckTYoip9iSNgH2BT4fEQsk/RHYPyJOlvQV4CcRUSdpPPC3iNgi2+/Q7Gd34Gpg34iYLKkP8EGTtzkeuCciDsk6r3tE0l3Za1uQukj5CHhO0u+zz/AX4IsR8aKkNSJisaTLgf2Bc4CvAlMior6UL9qsgc8krD15EthJ0m8lfSEi3mny+gjg/oiYGxELSF2xNHZTRCyOiKeBtbN1An4j6QlS1yMDGr1WyI6kju8mS3o8Wx7Sis/wKeDViJgMaRyMSN3XN/Y14GfZ8e8jdSuzfvba3RHxTkR8SOqaYgNSZ4H/jNThHRExN9v2QuDAbP4Q4KJW5DQDfCZh7UhETM3GZNgFOFXS3RFxcisO0bjn2oZBp/YHaoDPZmcGM0i/lIsRqYO541rxvq0lYM+IeG6pldLnWPozLKKZ/4cjYqak17MznBGkz2rWKj6TsHYju5w0PyIuJ41nsVWTTSYDX5LUN+sSfs8SDrsaMCcrEF8m/WXenLuBvSStlWVaQ60b7/s5oL+krbP9e2dZG7sd+IGURk+UtGULx3wI+KKkwQ2ZGr12AXA5cG1U59gGVuV8JmHtyWbA7yQtBhaQxtr+RETMlvQbUk+/c4FnSaMQNucK4BZJTwJ12T5FRcTTkk4A7sgGallA6or+pVI+QER8nDW4/15ST1J7RNNbY08htSM8kb3Hi8CuzRyzXtI44IZs+znATtnLk0iXmXypyZaLb4G1DkXSqhHxXvbX+Y2kgehvzDtXXrK7vs6OiC/kncXaJ19uso7mpKzB97+kv8BvyjlPbiT9DLgeKGf7iXVwPpMwM7OifCZhZmZFuUiYmVlRLhJmZlaUi4SZmRXlImFmZkX9Pzez59XVln6cAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# ------- simple example: Classifier and usage with Plotter class ---------\n", "\n", "# template of Classifier Class\n", "class CutClassifier(object):\n", " \"\"\"\n", " template implementation of a Classifier Class\n", " \"\"\"\n", "\n", " def fit(self, signal_data, bckgrd_data):\n", " \"\"\"\n", " set up classifier (\"training\")\n", " \"\"\"\n", " # some examples of what might be useful:\n", " # 1. signal and background histograms with same binning\n", " _, self.edges = np.histogramdd(np.vstack([signal_data, bckgrd_data]), bins=10)\n", " self.signal_hist, _ = np.histogramdd(signal_data, bins=self.edges)\n", " self.bckgrd_hist, _ = np.histogramdd(bckgrd_data, bins=self.edges)\n", "\n", " # 2. mean and covariance matrix\n", " self.signal_mean = np.mean(signal_data, axis=0)\n", " self.signal_cov = np.cov(signal_data.T)\n", " self.bckgrd_mean = np.mean(bckgrd_data, axis=0)\n", " self.bckgrd_cov = np.cov(bckgrd_data.T)\n", "\n", " def evaluate(self, x):\n", " # simple example of a cut-based classifier\n", " c = 0\n", " for i in range(len(x)):\n", " c += (x[i] < (self.signal_mean[i] + self.bckgrd_mean[i]) / 2.)\n", " return c\n", "\n", "# example how to use a Classifier Class with the Plotter Class\n", "ndim = 2\n", "\n", "# initialise Classifier with training data\n", "cut = CutClassifier()\n", "cut.fit(data[signal, :ndim], data[bckgrd, :ndim])\n", "\n", "# initialize Plotter Class\n", "plotter = Plotter(data[signal, :ndim], data[bckgrd, :ndim])\n", "# and make plots\n", "plotter.plot_contour(cut)\n", "plotter.plot_test_statistic(cut)\n", "plotter.plot_roc(cut)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As a very simple approach, we start by building a two-dimensional histogram of the likelihood ratio of the signal and background distributions. First, choose a 'suitable' binning, and then, in each bin, determine the signal-to-background ratio, which already is the desired likelihood ratio in the bin. \n", "\n", "Implement a class in analogy to the example provided in `CutClassifier`. In particular, implement the method `evaluate(x)` for your new class that returns the likelihood-ratio for each bin corresponding to the test point $x$.\n", "\n", "Use Iris-versicolor as the signal, the other species as background. Visualise the results in the two-dimensional parameter space using the methods of the provided class `Plotter`.\n", "\n", "What are the problems of this simple approach?\n", "What happens if you use Iris-setosa as the signal?\n", "\n", "**Hint**: Use the `numpy` function `histogram2d` (or `histogramdd`) to create a two-dimensional histogram. You may profit from the function `find_bin` defined in the provided code basis to find the bin which corresponds to a given point in variable space. " ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# -----------------------------------------------------------------\n", "# Exercise 2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 9.3: Logarithm of Likelihood-Ratio (obligatory)\n", " \n", "Next, we use parametrised distributions, assuming that the signal and background distributions can be reasonably described by two-dimensional Gaussian distributions, defined by the mean values and covariance matrices of the signal and background data-sets.\n", "\n", "Implement a class with a method `evaluate` returning the logarithm of the likelihood-ratio for each point in the sample.\n", "\n", "Visualise your results, as in the previous exercise, with Iris-versicolor as the signal. \n", "\n", "What is the advantage of this approach compared to the previous case? \n", "What happens if you use Iris-setosa as signal?\n", "\n", "**Hint**: The `numpy` functions `mean`, `cov`, `.dot` and `linalg.inv` are very helpful for this problem. " ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# ------------------------------------------------------------------\n", "# Exercise 3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " \n", "### Exercise 9.4: Fisher's discriminant (obligatory)\n", "\n", "In order to study a truly one-dimensional classifier, implement Fisher's discriminant. Write an appropriate class which calculates the Fisher vector in the `fit()` method, and write a method `evaluate(x)` returning the projection of the feature vector $x$ onto the fisher vector, and study the resulting plots. Use again Iris-versicolor as the signal.\n", " \n", "What is better in this approach, and where are possible limitations? What happens if you use Iris-setosa as signal?" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# ------------------------------------------------------------------\n", "# Exercise 4" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 9.5: (voluntary)\n", "\n", "With the code implemented so far, you are ready now to use the full four-dimensional data, _i.e._ to employ all four features (sepal length, sepal width, petal length, petal width), hopefully with only minor modifications or additions to your code. \n", "\n", "Is it useful to consider the approaches of exercises 9.2 and 9.3 (naive likelihood-ratio in histogram bins, or the parametrisations\n", "of signal an background by four-dimensional Gaussian distributions) in the four-dimensional feature space?\n", "\n", "With your implementation of Fisher's Method you are ready to separate better even Iris-virginica and Iris-versicolor. \n", "Give it a try!\n", "\n", "**Remark**: The basic code skeleton you have developed in this exercise can easily be used for other data-sets, and an extension by other methods, _e.g._ those mentioned in the lecture, is also thinkable." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# ------------------------------------------------------------------\n", "# Exercise 5\n", "# -> Set ndim to 4" ] } ], "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.5" } }, "nbformat": 4, "nbformat_minor": 4 }