{ "cells": [ { "cell_type": "markdown", "metadata": { "nbpresent": { "id": "11fab20f-a5cb-42ff-8112-19cb19457f46" } }, "source": [ "# Introduction to AI: assignment 2 - Regression" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this assignment you are going to apply linear regression methods to analyze a genomic dataset1: *genome_sizes.csv*.\n", "\n", "### Motivation\n", "The genome of an organism encodes all the information for producing proteins, which in turn are the actors of all cellular functions. The number of different coded proteins varies widely between species, and is a measure of the complexity of the chemical pathways of an organism. Moreover, not all genetic material gets translated into proteins: *coding* and *non-coding* parts are present, and their extension and number widely changes between different organisms.\n", "**Is there a relationship between genome length, number of coded proteins, and extension of coding/noncoding parts?**\n", "\n", "### Aim\n", "We will base our enquiry on the data provided by Hou & Lin in a study on genome sizes of a large variety of species1. We will first read, clean and visualize the dataset (part A), then we will analyze it by applying linear regression methods of increasing complexity (parts B, C, D).\n", "\n", "1The data is extracted from: *Hou Y, Lin S. Distinct gene number-genome size relationships for eukaryotes and non-eukaryotes: gene content estimation for dinoflagellate genomes. PLoS One. 2009;4(9):e6978*\n", "[doi:10.1371/journal.pone.0006978](https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0006978)" ] }, { "cell_type": "markdown", "metadata": { "nbpresent": { "id": "5e44b8fa-6165-4862-9c1c-3411063b4d93" } }, "source": [ "## A. Reading and visualizing the Hou-Lin data set" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Always have a look at the data first!**\n", "\n", "Open the genome_sizes.csv file with a text editor (or print it on a terminal with your favorite command: cat, less, more, head,...). \n", "1. How is the file formatted?\n", "2. Are there missing data?\n", "3. Are there header lines?" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To read the csv file genome_sizes.csv, you can use a library called [**pandas**](https://pandas.pydata.org/docs/getting_started/intro_tutorials/index.html). \n", "\n", "* Import the pandas library here, and use the function [*read_csv*](https://pandas.pydata.org/docs/reference/api/pandas.read_csv.html) to read the dataset into a DataFrame.\n", "\n", "* Read the documentation in order to import this file correctly: remember for example to consider the function arguments for parsing header lines.\n", "\n", "* Give these names to the columns in the file:\n", " 1. \"ScientificName\"\n", " 2. \"Domain\"\n", " 3. \"GenomeLength\"\n", " 4. \"ProteinCodingGenes\"2\n", " 5. \"Genes\"\n", " \n", "2Genes usually code for proteins, but can also code for types of RNA that do not get translated into polypeptides. For example, *ribosomal* and *transfer RNA* (rRNA and tRNA) genes code for special RNA machineries that do not have any corresponding protein." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import pandas as pd\n", "# Use pd.read_csv here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Show what is inside the data frame" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "How many different values does the \"Domain\" column take? Enumerate them using a DataFrame method" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you might have noticed, some lines contain NaN (Not a Number) values, which means that the table was not complete. Keeping NaN values in your table can soon generate issues: it is thus recommended to remove them, when this does not compromise the integrity of the dataset.\n", "\n", "Look for all lines containing NaN. How many do you find? How many lines does the file contain? Will the dataset be compromised if you eliminate all NaN-containing lines?\n", "If not, drop all the lines containing NaN values" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Enumerate NaN values in each column. Are they few enough not to change statistical significance of the dataset?\n", "\n", "# If so, delete the lines with NaN values\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We consider the number of protein coding genes as the output and the genome size and percentage of coding genome as features. Import numpy, generate feature matrix $\\mathbf{X}$ and generate output array $\\mathbf{y}$ " ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\n", "\n", "# Feature matrix X\n", "\n", "# Output array y" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Print first row to check" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Scatter plots for $\\mathbf{y}$ against each considered feature.\n", "Given the diversity of the organisms, quantities might need to be plotted in logarithmic scale!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Protein coding genes vs. genome length" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import matplotlib.pyplot as plt # This is the graphic library of Python\n", "# Command to make figures appear in Jupyter\n", "%matplotlib inline \n", "\n", "\n", "## plt.scatter(...) # Produces a scatter plot. It needs two arguments (Xs and Ys)\n", "## plt.title(\"Scatter plot of protein coding genes vs. genome length\") # Set the title\n", "## plt.xlabel(\"Genome length\") # Label the x axis\n", "## plt.ylabel(\"PCG\") # Label the y axis" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Protein coding genes vs. percentage of coding genome" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [], "source": [ "# Same..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now plot the two scatter plots at once in a single figure with the command **plt.subplots()**. You will have to specify that the Y axis is the same for the two plots." ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# There are two ways of using matplotlib: by calling functions (ex: plt.function())\n", "# or by defining figure and axis objects and using their methods (ex: ax.set_function())\n", "# Here, you will have to use the second way, which gives you more options\n", "\n", "## fig, axs = plt.subplots() \n", "## axs[0].scatter(...)\n", "## axs[1].scatter(...)\n", "## etc..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Are both the two features correlated to the output? What could be the cause?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Make the same two graphs, but this time color values belonging to each domain differently (one color for Eukaryotes, one for Bacteria, etc.). What do you notice?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## B. Simple linear regression" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To start, we will work only with data from **Eukaryotes** only. Extract the Eukaryotes data from the DataFrame and define the usual feature matrix and output vector.\n", "Import linear regression functions from Scikit-learn." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from sklearn import linear_model as lm\n", "from sklearn import metrics as met\n", "\n", "# Define X and y again" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Make a linear regression for the Genome length feature. Be careful: what are the variables that linearly correlate?" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Transform X and y so that they linearly correlate\n", "\n", "# Define the linear regressor object\n", "\n", "## linreg_1 = lm.LinearRegression()\n", "\n", "# Train the linear regressor\n", "\n", "## linreg_1.fit(...)\n", "\n", "# Calculate the MSE, RMSE, R^2\n", "\n", "## MSE_1 = met.mean_squared_error(...)\n", "## RMSE_1 = np.sqrt(...)\n", "## R_2_1 = met.r2_score(...)\n", "\n", "# Print the quantities you calculated\n", "\n", "## print('beta_0 = {}, beta_1 = {}'.format(linreg_1.intercept_[0],\n", "## linreg_1.coef_[0,0]))\n", "## print('MSE = ',MSE_1)\n", "## print('RMSE = ',RMSE_1)\n", "## print('R^2 = ',R_2_1)\n", "\n", "# Use the regressor to predict the regression segment. \n", "# What are the x-extrema of the segment?\n", "\n", "# Scatter plot the data and the prediction line\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Make a linear regression for the Coding genome percentage feature." ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Same..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## C. Multiple linear regression with 2 explanatory variables" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Linear regression for both features" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [], "source": [ "# Generalize what you did in B to account for 2 input dimensions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot predicted surface and data points for genome length and percentage of coding genome" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWQAAADuCAYAAAAOR30qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsfXlwJGd5/jO3jtF9r6TVfaz2kFerXe8GMAZSBhuy4Y84\nMVWx2ThAcOEDigpxoOKY8DM4xAbjA6cSG0wC2DiGKlPEcYFNmcP2HnYteO3VSqPRSKNjNKPR3GfP\ndPfvj62v/fWoZ9Qz06NpzfZTtRUiyTM9Rz/99vu+z/PoeJ6HBg0aNGgoP/TlPgANGjRo0HAZGiFr\n0KBBg0qgEbIGDRo0qAQaIWvQoEGDSqARsgYNGjSoBBoha9CgQYNKoBGyBg0aNKgEGiFr0KBBg0qg\nEbIGDRo0qATGPP9ek/Vp0KBBQ/7QyfkjrULWoEGDBpVAI2QNGjRoUAk0QtagQYMGlUAjZA0aNGhQ\nCTRC1qBBgwaVQCNkDRo0aFAJNELWoEGDBpVAI2QNGjRoUAk0QtagQYMGlUAjZA0aNGhQCTRC1qBB\ngwaVQCNkDRo0aFAJ8jUX0qAhJ3ieB8uyAACDwQCdTpanigYNGqARsgaFwHEcGIYBz/NIJpPCzzc2\nNtDV1QWDwQCDwQC9Xg+9Xg+dTqeRtQYNGdBaFhqKAiHiZDKJ1157DQCQSCSg0+mg1+vhdDrBsiyS\nySRisRgikQhCoRBCoRDC4TBisRgYhkE6nQbHceB5zeFVw5ULXZ4ngHa2aADP8+B5HqlUChzHCT/7\n3e9+B7PZDJPJJFTLiUQCPT09qKmpQW1tLWpqamAwGITHyPz+ESI3GAwwGo1aRa2hUiDry6u1LDTI\nBs/z4DhOqGaByxXy2toalpeXwXEcpqamBPLkeR5nz55FfX09otEofD4fYrEYOI5DVVWVQNC1tbWo\nra0ViJrjOLAsC4ZhRM9P2h50+0Mjag2VBI2QNWyLTCLW6XRgWRbLy8tYW1tDV1cXjh07hnPnzsFi\nsQhESqrd1tZWtLa2ih4vmUwiGo0iGo1idXUVsVgMLMvCYrFsIWqj0Sgiap7nRSRME3Rmn1qDht0E\njZA1ZAXZmEin0wIJplIpOJ1OeDwe9PT04MSJEzAYDHk9rk6nQ1VVFaqqqtDS0iJ6PoZhBKJ2uVyI\nRqNgWRZms1kgaELYJpMpK1F7vV40NTWhqqpKcqCoQYMaoRGyhi2QImKGYeBwOOD3+7F3716cOHEC\ner2yM2GdTgeLxQKLxYLm5mbR8aRSKYGo19fXEY1GkU6nYTKZRCRdW1sLs9kMj8cDq9UKg8EAhmFE\nJEwqaY2oNagNGiFrEMDzPNLptKjajMfjcDgcCIfDGBgYwPj4+LbEpTSx6XQ6mM1mmM1mNDU1iX7H\nMAxisRii0Sg2NjawuLiIVCqFZDKJxcVFNDQ0iIiaIJ1OI5VKiR5LI2oN5YZGyBoEIk6n0wAuE2Ak\nEsHCwgIYhsHAwAD279+vSmIiRN3Y2Cj6+R//+Ee0tbUhnU5jc3MTTqcTDMPAYDCIWh9yiFqn023Z\n+tBELxpKAY2Qr2BwHIdYLIZwOIzGxkbodDoEAgEsLCwAAAYHB7dUpLsFer0edXV1qK6uFv08nU4L\nFbXP58Py8jKSyST0ev2WHnVVVZXw32W2cHiez1lRa2StoRBohHwFgmxMsCyLcDiMtbU1sCwLh8MB\ns9mMkZER1NfXF/z4mVsQaoLRaER9ff2W18eyrEDUgUAAq6urSCQS0Ov1oo2PmpoaEckToqaxsbGB\nzs5OGI1Gjag15AWNkK8QSIk5ACAQCMDj8QAA9u/fj9ra2qKfZzfCYDCgrq4OdXV1op/TRB0KheBy\nuZBIJABARNS1tbWoqqqCTqeD0+lEe3u7SEIOaKIXDdtDI+QKh5SYg+d5uFwuOJ1OWK1WNDU14cCB\nA0U/F7mVJ89BSEbNFfN2yEbUpN1D5OButxvxeBzAZen40tKSQNTV1dXQ6/Wa6EXDttAIuUIhJebg\nOA4rKytYXV1Fe3s7pqenkUqlYLPZFHlOQjo0MZP/XWnEotfrYbVaYbVaRT/nOA7nzp2D1WoVNj9i\nsRh4nkd1dfWWPnUmUWcTvdCtD23zo3KhEXKFQWqHOJ1OY3l5GS6XC93d3bj66qthNF7+6AkJKAGa\niHP9rJJBCLOtrQ1tbW3Cz3meRzweF3apNzc3t8jIaaLOlJFnErW2oleZ0Ai5QpBNVbe4uAiv14ve\n3l5JVR2pnJWARsjZodPpUFNTg5qami1EnUgkEI1GEYvF4Pf7EY1GwXGcICOniZrIyMk8IFP0kkgk\nYLFYUF1drRH1LoRGyLscUmKORCKBxcVFBAIB9Pf3Y2RkJKuqTq/XK07IV/rJn88FSKfTobq6est6\nXi6/j0wZOe33sba2hubm5i2fgVZR7w5ohLxLISXmiMViWFhYQCwWw8DAAPbt27ftCVdqQtYq5MIg\n1+8jU0aeSqUEIUttba3g9wFkVyfq9XrJFT0NOw+NkHcZiCE8cVnT6XQIhUJYWFhAOp3G4OCgZIWU\nDUoSMnksuhq/Egm5lGSWze8DuCwjn5mZAQB4PB5Eo1GkUikYjcYtFbXJZBL+O030oh5ohLxLQIs5\nWJaF0+lEVVUVFhYWYDAYMDg4uEU+LAdKEqbWQy4vSDhAZ2enaPuDNmai/T62k5FLiV7ILrUmeikN\nNEJWMbKJOXw+HyKRCFZXVzE+Pr5lRzYfKHkSaYRcfmFM5h0KAJhMJjQ2Nm65YKfTadHWB5GRGwyG\nLaIXi8Uieo5sohee52GxWLRd6gKhEbIKkU3M4Xa7sbS0hPr6elRXV+PQoUNlPlIx6B5y5h7ylYJy\nDzWlCDkbjEYjGhoa0NDQIPo57ffh9/uzyshpdSL5zr7zzjvo7+9HTU2N8Hia6EU+NEJWEbKJOUhE\nUmtrK6ampmCxWIRAUTVBp9PB7XbD4/EIZj3xeBzBYBBGo1HUt6xU7CZCzgY5fh/BYBBra2tCoC0h\n6kQiAZZlBdKld6kJyHukiV62QiNkFUBqh5hlWUFVRyKS1Epo5KLh8XjQ1taGyclJ8DyPWCyGUCiE\nzc1NrK2tiQzlrVaraGWrUlAJhJwN28nIo9EoGIaB0+kUWhpk95qWkWuil+yonDNhFyKbqm5paQke\njwfd3d04fvy4agmL4zisrq7C6XSira0N7e3t6OnpgdlsBsdxaGxsRG1tLfbu3SsMmcjKViQSEUU0\nZYogSOjpbkO52zOlJORsoGXkTqcTBw8eFO7uaHWix+MR/D6qq6u3uOiRHnQ20cuVQNTqPNMrHFJi\nDoZhsLi4CJ/PV7KIJKXAsixWV1exvLyM9vZ2oXqfmZnZdqgnlfxBdmsjkYgggiBqNVpWbLVahRNX\nzSgnQZS7Qqefn/aYpsFxnKBOjEaj8Hq9Ir+PTKImFTVQ+UkvGiHvIKTEHIlEAgsLCwiHw+jv78fY\n2JjsL9JOn3wkaXplZUWyjVLolgW9W5spgqBPXJ/Ph2g0KjLqIa0P4qhWbpAd3nJC7UREhoO5ZOTk\n8870+6DJWoqoySCyp6cHer0e3/nOd3D33Xertt2XCY2QdwAcxyEcDgtfLDoiKZFIYHBwMO+IJIPB\nAJZlFWtn5CJ3Yk60urqKPXv2ZG2jKL32RsuKW1tbRcdKboUjkciWW2Gr1YpkMol4PC683zuFcleo\n5UYxrz3X551NRk5aXYSoGYZBKpUSLoo//elP8ZWvfKXo17VT0Ai5hKDFHC6XSxhoLSwsgOd5QVVX\nCEqhsMvs2abTaTidTqytrcnqZ5MeII1SrL1lM+qhe5bk2Ofn53OuaykNpT4TDe9Crozc5XIhEAgg\nnU7jj3/8I375y1+C53m8/vrrmJiYyCuO7NZbb8UvfvELtLe34+23397ye57ncdddd+GFF15ATU0N\nnnrqKUxNTRX9WjVCVhjZxByJRALLy8uor68vOiIJKC0hp1IpOJ1OuFwu9PT0SLrESYEMcsrlZUH3\nLN1uN4aHh1FdXZ11XSszR89qtcJsNhdN1FdqhZz52ZcaUjJyEke2b98+WCwWnD9/Hj/+8Y/xzjvv\n4AMf+AD++Z//WdZjnzp1CrfffjtuueUWyd//3//9H2w2G2w2G86cOYPbbrsNZ86cKfo1aYSsEKTE\nHMBlT4HFxUXo9Xq0t7djfHxckecjLQslQAg5lUphaWkJbrc7q11nLqhVqZcrnolUV36/HysrK4JS\nLXM1j5YU58KV3LIg+8flRDqdhtlsRktLCz784Q/jgQcewGOPPZb341xzzTVYXFzM+vvnn38et9xy\nC3Q6HY4fP45AIACXy4Wurq4ijl4j5KIhJebgeR7r6+tYWlpCU1MTJicnEQ6HEQwGFXteJStkAFhY\nWCh6w0OthJwNBoNBUgBBS4pp7wej0Sgi6UyTHkAj5HKvaKbTaUElGAwGt6gQlcLq6ip6e3uF/7+n\np0fQDBQDjZALhNQOMdnLXVlZQVtbG6anp4XKigwhlIISFTLDMHA4HAgEAmhubi561S5Xpt5uQjZJ\ncSqVElbz3G630Kum/Yml+ug7hXL3r9VSIZOLQigUKshwSw6kPmMlLsQaIeeJbKq65eVlrK2tYc+e\nPaKIJAIlWwxAcRVyMpmEw+GAz+dDX18fWltb0d7eXvTJJJU+UknVoslkQlNT05YdapqoNzc3EQgE\ncPbsWZHYhexQl1LsUu6VO5Zlyy7moQk5EAiUrELu6enB8vKy8P+vrKxgz549RT+uRsgyISXmID3X\njY2NbYdfShNyIY+XSCTgcDjg9/tFO88+n0+R6mqntizUBJ1OB7PZjObmZjQ3N6OhoQEmkwljY2Oi\nVa3l5eUtO7Wk/aGU2KUcKj0aaiPkUrYsTp48iUcffRQ33XQTzpw5g4aGhqLbFYBGyNtCSsxBKsxA\nIIC+vj4MDw9veyKUs0KOx+NwOBwIBoMYGBjA+Pj4FkmqEoS823rIpQC5WOda1SLih0gkIlKpZa7m\n5St20Qh5KyEX2rL4xCc+gVdeeQVerxc9PT346le/KigEP/vZz+KGG27ACy+8gOHhYdTU1OD73/++\nIsevEXIWEOMTmoij0SgcDgei0ajsiCQCo9G4xey7GMgh+FgsBofDgVAohMHBwazHq1TfM1vL4kok\n5GzIJn4gcmLS+qDFLlJELfUc5SZkqV32nQZ9USiGkJ9++umcv9fpdAVtb2wHjZAzQDYmVlZWYDAY\n0NHRIYpIGhgYQEtLS9690Z2skGOxGOx2O6LRKAYHBzExMZHzeJWukK/kTL1CtyxoOTEN2kktHA5j\nfX0d8XhcJHYhrY9yD9XUUCHT738gEBBtQuwGaIQMaTEHy7IIhUJYWVmBXq8vOCKJQOk1NSmCj0aj\nsNvtiMfjGBwcRGtrqyxykKpsC4HWslB+7Y12UqNBi10CgQBWVlYQj8eRSqUwMzOzZYd6J4araiBk\nGqXcsigVrmhCzibm8Hq9WFpagl6vx6FDh4qKSCJQ+oTQ6/VCTysSicButyORSGBoaCjvCl6pi4XU\n42iEXBpIiV2ICrG7uxuRSASbm5twOp1gGEYy6FSu2EUuWJYtq4lP5nsfCoVKNtQrFa5IQiarayzL\niuSebrcbi4uLqK+vx8DAAGKxmCJkXAoYDAbEYjH84Q9/QCqVyjttmoY21FMO5RSGcByXNe0jlUoh\nFoshEomIxC7EX4VufRQq7ih3hZz5/MX0kMuFK4qQpXaIeZ4XvH2bm5tx+PBhVFVVwe/3IxQKlfuQ\nJREKheBwOJBMJnHw4MGCDYoISk3I5RYs7CTKTcjZesgmk0lS7EKb86yvr0uKXQhRb0e25SZkesMC\nuEzI+RgKqQFXBCFnU9WRiKSOjg6Rqg5QfiuCPpZCT9hgMIj5+XkAQFdXF2KxWNFkDCi3ZQFc7mO/\n8cYbYBgGVVVVAC6/l42NjbvCXL5YqJWQsyFXYEA2u0u6P02LXTRCLh4VTchSO8TEltHtdmdV1QGl\nIWRyIcj3S+v3+2G326HX6zE8PIyGhgYEAgFEIhFFj6sY+Hw+2Gw28DyPQ4cOwWQyIZVKCYKIpaUl\nxGIxAO96FpMTe6c9i0uJ3UbIUpByUQPEvsSRSERkIF9dXY1YLCZcdMtx8c0k5Gg0uiWtRO2oSELO\nJuZYWlqC1+uVZaBTCkI2Go15VRE+nw92ux1GoxGjo6OivqDSbm+ZsThy4ff7MT8/D7PZjP7+foRC\nIVitVqFCrq+vh9lsRl9fH4B3PYsjkYjICjPTYc1qte6alAcalUDI2ZBL7BKPx/HOO++AYZgtF1+6\n9ZFth1oJ0IRM7vh22x1ZRRGylJiDqNRCoRD6+vowMjIi60MyGAyKE7IcEuV5XiBis9mM8fFxycFi\nKfyQ80EgEMD8/DyMRqNwjH6/f4ujXWZfmfYg7ujoEH5OHNbI0MnhcCCVSgl+EKSiLrUfRLGoZELO\nBhIYYDQa0dfXJ1xIc4WcZu5QK3GXlFkhE8XkbkJFEDLLskgmk8KbnxmRNDAwsK04IhOl+GLnImSe\n57G5uQm73Y6qqipMTExs2T2V+1j5Ih9CJn1svV6PsbEx0cWimC0LKYc1upcZiUREfhCZbY9SVl75\noNyEXE77y8y7v1whp3ICA2pra2GxWGS/nzQh77RZvlLY1YRMTthUKoU333wTV199NYLBIBYWFsBx\nHAYHB9HU1KSaD0aKRHmeh9frxcLCAqqrq3HgwAFZfS+lK+TtSDMUCmF+fh48zwt97EyUIlMvWy+T\ntD0y1WvkRGYYBgzDoLq6uqDnLhTlJuRyftflPj8tdqHvkojYJRKJSAYG0K0Pk8m05bloL+RQKKTa\nldVc2NWETMeNMwyDN954A0ajEUNDQ6pcCKfbIDzPY2NjAwsLC7BarTh48OAW2Wwu7FTLIhwOY35+\nHhzHYWhoKOdeZzky9drb24Wf0wkgDMNgfn5eECvQbQ85K1yFopwWmGrwkijmgpAt2YUODMgmdiHB\ntmTOEggEdt0OMrDLCRkANjY2YLfbkUqlMDU1lfM2vxAoWfEQQl5fX4fD4UB9fT0OHTqUFxHTj6VU\ny0JqyyISiWB+fh6pVArDw8Oy1ofKLQKhE0C8Xq+QqUe3PVZXVxGNRgUbTHqImK+7mhSu5Aq5VMgV\nGED3pzc2NuDz+fDyyy9jdnYWgUAAp0+fxsTERF4Zli+++CLuuususCyLT33qU7j77rtFv3c6nfjk\nJz+JQCAAlmVx//3344YbblDmtSryKGUEx3GYnJzE+fPnFSdjsmmhxLSf3Ga7XC60tbXhqquuKup2\nWmm/BELIRIbNMAyGhoby2nNWq9tbtl1b2l1tY2NDGDjRVVe+waflJORyG9Tv9Os2mUxobGwUKmGG\nYTA4OIiBgQE8++yzeOmll/DUU0/h4sWL+PrXv473vve92z4my7L43Oc+h1/96lfo6enB0aNHcfLk\nSUxMTAh/8//+3//DX/7lX+K2227DxYsXccMNN+TM38sHu56Qu7q6hMpA6SmzEoTM8zxcLhcWFxdh\nNBqxd+9eDAwMKHaMSoC0fN566y3E43EMDw+L1prkgiZfesBabkKWAm2D2dbWJvyc4zjJ4FNye0xX\n1FIDtHJXyOUiZDV8xmSo19nZiYGBAbzvfe/Dv/zLv+T1GGfPnsXw8DAGBwcBADfddBOef/55ESHr\ndDpBxRsMBhVJCiHY9YRMvvxEiGCxWBR77GJ2kTmOg8vlwtLSEpqbm3HkyBFsbGwoasGpBIhVp8/n\nw+TkZEHWogSV4GWh1+sl+5jk9jgSicDtdsNut29Rrlmt1rK2DcpJyOW2/gQuEzLthVzIHEkqvPTM\nmTOiv7n33ntx3XXX4ZFHHkE0GsVLL71U3IFT2PWETEDIU0lCLqRPy3Ec1tbWsLS0hNbWVpEk22Aw\ngGEYxY6vGMTjcdjtdkQiEfT09ACAyDC9EFQCIWdD5u0x8K5yjc7T8/v9MBgM2NjYUHzPdjuUk5DL\nvXJHjoEQciAQEG1wyIWc8NKnn34ap06dwhe/+EW8/vrruPnmm/H2228r8t7vekLOrJCVRD4VMvHG\nWF5eRnt7O44ePbrF3lBpk/pCkEgksLCwgGAwiKGhIezfvx+JRAIej6fox6Z70XSUUSUQshRo5Rq5\nmC0sLKCurg7V1dVb9mxLrUYsJyGn0+myV8g0QqEQRkdH8/7v5ISXPvnkk3jxxRcBACdOnEAikYDX\n6xVt/BSKXU/IBKWSOm/3mCzLYmVlBSsrK+jo6MCxY8eynmSlUP/J7VmSgNNAILAlzqkUiSFq7yGX\nCmSwJrVnK6VGJM5qNEkXqkYsd4Vc7pU7GoUaCx09ehQ2mw0OhwPd3d145pln8OMf/1j0N3v37sXL\nL7+MU6dOYWZmBolEQjSHKAa7npDLVSGzLIvl5WWsrKygq6srq0lR5uOVIsYp14lAAll9Ph8GBwe3\nBJzSj1MsKrlloQRyqRFJ26MYNWK5e8jlJOTM116oF7LRaMSjjz6KD3/4w2BZFrfeeiv279+Pe+65\nB9PT0zh58iQefPBBfPrTn8a3v/1t6HQ6PPXUU4q1o3Y9IRMYjcaSEHLmY6bTaSwvL2N1dRV79uzB\n8ePHZffOlG5ZkMeTOhEYhoHD4cDm5ib6+/sxNjaW9UujEbJyyHfLglYjZhr2ZMvSy2x7kNbYlUzI\n9EAPKM5684YbbtiyV0xva0xMTODVV18t7EC3QcUQsslkEvZIlYLRaBQek9h2koicfIiYYCeCThmG\nweLiIrxeL/r7+zE6OrotQZSCkK+EHrIUlFp70+l0wj60lBqRRDQtLS0JyR+xWAxutxt1dXUlVSNK\nQQ2ETJ+PuzG+CagAQqZbFkonfBgMBqRSKdjtdqyvr6OnpwcnTpwo+ItXqgoZuLyWtbi4CI/Hg76+\nPhw/flx2taTU7dZu2kMuFUq9h0yrEWkwDIM333wTLMtKqhHptkcpqmi1EXI4HM5LnacW7HpCJlB6\nqJdKpeByueDxeDAyMoLjx48X/YUrRYXMMAxcLhfcbrcsn+dSQoqINELeGZjNZkF4RB8LUSNGIhGR\n/aVU26OY41YbIfM8r6oho1zsekJWeqhHbvk3NjbQ0dGBlpYW0Ze8GChpCJROpxGJRHDhwgX09/eX\nlYgzQdboSGCmRsjlQTY1YjZXNVqNmG/gKcuyiqdY5wMpc/rdiF1PyATFVsjJZFLovfb19eHEiRNg\nGAbvvPOOYseoxIlKhopra2tCkkixgg6lwDAMEokEzp8/j9bWViFmKhaL4e23367Y2CYaaiLkbMjm\nqkarEdfX1xGJRCTViFLxTGqqkIlaUu2fgxQqgpB1Ol3BFTK9Ftbf3y9KFClV0GkhoNfsenp6cPz4\nccGfuNxIp9NYXFyE2+2GXq/H8ePHkUqloNPpkEqlcOHCBQwMDEjGNhGSJv92420mjd1AyNmwnRqR\nDBJJPBOd+kFbX5YDtOdMJBLZlV7IQIUQMpB/cjIRSvj9fgwMDEiuhalBWUcLTzLX7AwGg2ItkEKP\nzel0YnV1Vehfnz59WvQ+kv8tFdtEV2QulwuRSGTL/q3Vat1V1fRuJmQpSKkRga2pHz6fD36/H2az\nuSzZiOl0WnBPDAaDu3KgB1QIIeczOCIZe8FgEAMDA5JCCfpxSwE5Jy0txe7s7JQUnijZk84HHMdh\ndXUVTqcTXV1dBW+eZKvI6DQQl8slWU3n09/cSVQaIWdDphoxmUxi7969qKqqKqkaMRtYlhW+D4Ua\nC6kB6vtGF4lsJ0QsFsPCwgLC4fAW6fBOYjt1HU12HR0dORWASlfw25EJbSXa1taWUyZOkO+WRbY0\nECI7JiIJ0t+k17pKnWosB+Ui5HKb05Mecr5qxMyw00I/v1QqJZwnuzUtBKgQQs70ZKDJLhqNYmFh\nAdFoFIODg9i/f39Zv7jZ1HW0S1x7e7ssslOyQqZ9KDJB4qbsdjsaGxtx5MgR2a56Sq29ZTvR6bUu\nt9stqNkYhsH6+jqampqyeheXAlc6IUtBSTViruenK2SNkFUAMtgzGAxC6nQ8Hsfg4CBaW1sL+sIq\nbXyfWdXSVWdra6ukS1yux0omk4ocFyH3zNe5ubmJ+fl51NbWFpRyUso95GxrXel0Gm+99RZ0Op3I\nu3gnqulyEXK500IK2bIoRI1IkzStRqQrZI2QVQKj0YhgMIjZ2Vkkk0khgqiYE4RsWii1Y0kImed5\nIVuvpaVF5JssF6UMOg0EArDZbDCZTLKTsKVQDnIyGo0wmUzo6OgQLiBSIolYLCaqxsi/YhNiylUh\nl5OQlS5asqkRSduDViNWV1cjFovB5/MhEonA5/MpYoVZDlQEIet0OoTDYQQCAUSjUYyPj+eVBZcL\nShOyXq+H2+2G2+0WkkQKNdVXsodMtlTolOmxsbFdO63ORC6RhNQQiuze0tW0HMK5UgkZKP3F12w2\no7m5WXRukyHwH/7wB8RiMXzjG9/AmTNnYDQa8bvf/Q4HDx7EF77wBdnV+3YBpwDw7LPP4t5774VO\np8Pk5OQWe85iUBGE7Ha74XA40NTUhM7OTsXIGFCO9HieF5JxGxsbMTU1haqqqqIeU8kKmeM4XLx4\nEalUCiMjIwU7ZQG7a9NAqhrL3L0lAajkFjtXNX0lE3I5QIbAJpMJAwMDeOKJJ/DlL38ZH/nIR9DT\n04NLly7JJmM5Aac2mw3f+MY38Oqrr6KpqUmRYAcaFUHIROJMqhslUaw4hB6INTQ0oKurCy0tLUWT\nMaDMxSKRSMButyMUCmFsbAzd3d1FEUoleFdk272lq2mv1yuqpglRp9Ppsrz+K5WQgcvnGP2eB4NB\ntLe3Y3JyEpOTk7IfR07A6X/+53/ic5/7nFCwKN0aqQhCVkuMEw2e5+H1emG322G1WoWB2MLCgqJt\nhkIrZNoveXBwEOl0Gg0NDUVXd4SQd1OVLBfZqmnS24xEIoJ0XKqaLqXXw5WcOJ25WVWoF7KcgNO5\nuTkAwHve8x6wLIt7770XH/nIRwo88q2oKEI2Go2KbR0Q5EvIPM9jc3MTdrsdNTU1OHToEGpqaoTf\nK933zfexaJkz7Zfs9XpLZlJfychc6fJ4PJiengbP85KbAkQgQf5J+UIUgnITcjkvvplOb4VuWcgJ\nOE2n07CuO8QQAAAgAElEQVTZbHjllVewsrKC973vfXj77bcV2+qoCEImKHeFTIi4qqoq62aCkoSc\nj3RaSuZMn8D5Ss+zoVzqQTWBGNtIGfjQKdVLS0siXwhiLG+1WvMe9F7paSGZXsiFKPXkBJwSHxnS\nsx4bG4PNZsPRo0cLfwEUKoqQSxV0ul3V7fP5YLfbYTabMTExAavVmvVvldwdlkPucmXOZN+6WNAC\nk3JXTmqElECC+EIQO8zl5WUwDCPs3dJy8WykqxHyu1TGcVxBQiA5Aacf//jH8fTTT+PUqVPwer2Y\nm5sTes5KoCIIudQ95Gg0Kvk7v9+P+fl5mEwmjI+Py3KYUrplkY1E85U5KxnjRP7vldS6KAa0LwQN\nuje9vLwsfA9ramq29KY1Qn7XC7nQ752cgNMPf/jD+OUvf4mJiQkYDAb827/9m+jiWiwqgpCB4iw4\nc8FgMGypuoPBIObn56HX62UTMYGSydNSJFqozLmUQacaCoPU3m1mNU3M5TmOg9lsFqnZdook1UbI\nQOE70dsFnOp0OnzrW9/Ct771rQKPNjcqhpCB0rUsCIGGQiHYbDYAwPDwcEF9qlKIOQiKkTlrhLw7\nkK2aXlpaAsMwWzL1Mqtpi8WieBtJTYQci8UKVpWqARVDyKVKCCDJ0+fPnwfHcQUTMUEpPJaDwSDm\n5uaKkjkrScgcxyEajcJkMpU11udKAiHqrq4u4WccxwlWpsFgEKurq0JUU2ZvuhhCVQMhk73+QCCw\na603gQoi5FIgHA5jbm4OoVAIU1NTRanXCJQk5HA4jFgshvn5+aJlzkpVtizL4p133hGIOZ1OIx6P\nY2FhQTUWmZUIqR4y8emQCgYgvenMaprenZYbDKAGQq4EYyGgQgm52Ol+JBKB3W4HwzAYGBhAKpVS\nhIwBZQiZkHAikYDZbMaRI0eKPq5iK+R4PA6bzYZQKITR0VG0tbUJBH/27FlYrVZEo1HBIpMYztfV\n1e14z7MSkc9Qz2QyoampSfSdlhsMIBWzpTZC1ipkFYAQMCG8QtZeotEo7HY7EokEhoaG0NLSAp7n\nBXWOEpAaEsoFkTmHw2EMDw+jpaUFr7/+uiLHVYjIBLi8CbCwsAC/34/h4WHodDpRf1On00Gv16Ot\nrW2L4Tw5+UmVxvO8qOdZV1entTxkotgti1zBAKSadrlciEajYFlWFLOVSCQUsQIoFBohqxhk0yIf\nQo7FYrDb7YjFYgIR0+tbSqIQAUamzHliYkJ0XErs++r1+rw2VFiWxdLSElwuF/r7+4VMQo/Hs+X1\nSe0kG43GLfFNmRsETqdTULeRSpqo27SWhxilWnuT+pwyq2mPxwOO4+ByuYS7HSJy2YlgAK1loUJk\n7iLL2TCIx+Ow2+2IRCIYGhoq2MS+kOOUg2wyZxrbRULJhdyWBS006e7uxvHjx0XPnSkMoX8m5xgy\nNwhor4hwOIyNjQ3Bx5iupNWasbdT2Mk95Mxqmud5NDY2oqGhQZCL72TMlkbIKoac1Tfa4WxoaEhW\nrNNOqs5YlsXy8jJWVlYkZc40lCTkXKRJ7EPtdjtaW1uzCk2kiL2YgWG2+B+WZYVbafrkr66uRl1d\nHZLJJJLJ5K5KrC4GahCG5BuzldmbLvSCyrKs8NqDwaDIIGi3oWIIWY5aL5FIYGFhAcFgUPLWPxuU\nIr3tUEiaM+mZFxu1nks67ff7MTc3h9ra2m19nKXItxS7yQaDQfLkj8fjCIfDSKfTWFhYQDqdFkmQ\n6+rqFDP0URPKmam3XZ5etmAAOr2FfFZ0NV1bWyu7PUX+RquQVQapCjmZTAqDp0ISp8ljKknIdMVd\nSJozgVL7w1KPEw6HYbPZoNPpsH///pweHQQ7RcjZnpvcSns8HgwPD6O6ulokQSaGPnLM5vNFOQUx\n5czUK2TLItsFlQ4G8Hg8soIB6HM5FApphKwm0BUyPQwbGBjA+Ph4QVUEIeRCo5YyQQeKFprmTKDU\nXjNNyPF4HPPz84jH43mnh5STkLNBSoLMsixisZjQlyZm86RCI0PE3dLyUEPLoljICQbIjNmqra1F\nKpUShEihUEjbslAD6JYFEXR4vd6sw7B8oLQk22AwwOv1YnFxseA0ZwIlK+R0Oo3Z2Vlsbm5ieHgY\nbW1teb9vUr3ochOyFAwGwxZ7TLrfSe/i0so20vKQIqByuttVAiFnQ66YrWAwiI2NDVy8eBGf/exn\nEY1Gce+99+L48eM4ceIErr76alnPISdLDwCee+453HjjjTh37hymp6cVeX00KoaQgcsKpPX1dWxs\nbGBsbAzHjx9X5EtazO5wJoLBIMLhMFZWVopKc6aPrVhCZlkWa2tr8Hg8GB8fL+oCJtWLViMhSyFb\nv5NWtq2srCASiQB418OYFkxcqYS8089Nqmme52G1WnHw4EG8+eabeP/734/Pf/7zuHTpEmZnZ2UR\nspwsPeByC+/hhx+WTfKFoGIImeM4vPHGG2htbUVrayt6enoUe2wlHNroNOeGhgaMjo4qYoJSqKAD\nuFxlrK6uYmlpCS0tLYq8b1I7x7uFkLNBStlGdqbD4bCQCMIwDJLJpBDbVVdXt2My8XJn6pXruaW8\nkK+++mocP35c9mPIydIDgH/6p3/Cl770JTzwwAPKHLwEKoaQDQYDTpw4gXg8jpmZGUUfu5iWBS1z\nJv3Yd955p6y5esSic35+Hi0tLTh27BhSqRRmZ2eLPh419pBLAamdaYZhcOHCBTQ0NIiGUnLkx8Xi\nSg0DkLLezBdysvTOnz+P5eVlfOxjH9MIWS5K5YlcCCFLyZwz5d1KIN/H8vv9sNlsqKmpEa2wsSxb\nMvvNSiTkbDAYDMJdGkE6nUY0GhX60pFIRGSNSdoeZrO56MTvKw2ZFTKQ//uwXZYex3H4whe+gKee\neqqgY8wHFUXIQGnsLYkFpxxsJ3NW+hjltiwikQhsNht4nse+ffu2mOorbb+Z+bMrAdmq1GyCCSIT\nDwQCotimTJl4pe1MKwmakBOJREHD8e2y9MLhMN5++21ce+21AID19XWcPHkSP//5zxUf7FUUIZfS\nE3m7ClmOzJlgJ4NOE4kE5ufnEY1GMTo6mnWFTcltDfLaMnetKx35tA3Ibm2mNSbDMAiHw0JaNb0z\nTRO1WmTi5f5caUIOBAIFWdBul6XX0NAAr9cr/P/XXnstHnjgAW3LolzIRchE5kz6ULlkzgRKV8hS\nLZpUKgWHwwGv1ytLHq5U6rRUe0Kpx1Y7lOjjms1mtLS0bJGJkz1ct9sNu92+xXGN47iy9JHVYL1J\n2m6FqvTkZOntFCqKkGmHNiWnzlKEnClzzjTZyQUl1+gMBgMSiYTouJaWlrC2toa+vj7Zq39Kp05v\n97NKRKkIMdseLu24lkgkcO7cOWFnmvaZLmXLQw2ErIT15nZZejReeeWVgp5DDiqKkAkIgSrlpUtX\ntETm7HA40N7enpfMmT6+ZDKpyLGRVgPP81hbW8Pi4mLeFwhAuT4vId9AIIBkMimQiEbIyoKWibe1\ntWFzcxNHjx4VdqbD4bAoqZqWHtfV1RUtEydQEyEHAoFdLZsGNEKW/XipVEpwO2tsbMT09HTBUmql\nWxbhcBinT59GU1MTjh49WlZT92QyidXVVQQCAVRVVWF1dRXhcBjBYBDNzc1C5VaJMU7lWj2jnzfb\nzjRpeWxubmJxcVGQHtN96UI+EzUQMnn+3W4sBFQYIctxfCsEoVAIwWAQHo+nKJkzgVKEHAgEMDc3\nB5ZlMT09XfRxFYNkMon5+Xn4/X40NjZibGwMDMPAYDDAbrejvr4eRqMR4XBYsGDc6dvrUqNchLxd\ne06v10vKxImRD/2Z0DvTxGc6F+GqgZBJta8RskphMpkU6dHSac5VVVU4cOCAAkdXPCFHo1HMzc2B\n4zgMDAzA6/WWjYxZlsXi4iLW19cxODiI1tZWrKysIBAIoKamBuFwGIFAAK2trairq0N9fb0gMU6l\nUoLazel0btkoIES9W7L21ErIUshm5ENHNsmJ1io3IdPPHwqF0NnZWbZjUQIVRcjkZCAthkJBy5xJ\nmvNrr72m1GEWPNQjYpNIJIKRkRE0NzcLirCdBt2zJskhwOXtjng8jrW1NWxuboLjODQ2NgonttVq\nBc/zAomQynjPnj3CyhzxNM4UUdAkrcasvXK2LJQcYEtFa5HPxO/3CzvTZrMZer0eer0e0Wi0LNFa\n9GvXKmSVotCWhZTMuRTIt0JOp9NwOBzY2NjYIjYphRBmO2xubmJubk7oWdO70DqdDgzDIBqNYv/+\n/WhubhaqYJ/Ph6WlJaTTaSHZg25VkNdBDH46OjqEISFxYaN7oLRVZl1dHSwWS1n70uUi5FKb++j1\nemFnmkYymRTubBwOh+BdrFQSSL7QCFlloCvkfLYYcsmcaSh1wsklUY7j4HQ6sbq6ir1790qusCkl\n6JCDSCSC2dlZGAwGHDp0CFVVVeA4Tnh+l8sFp9OJnp4eHDt2TDhWQphE/URUauFwGKFQCKurq0gm\nk6IhU21tLSwWi/DYFosFZrNZlHtIZ+2tra0hmUwKSrdEIoFYLLajfsblbFmU43ktFguqq6uFOxzg\n3Z1p0pfO3Jkmn28pLp4aIasUJpNJsEjMBTkyZwJCokpc7bdzj6NX6zo7O3H11VdnfV6lK2QpUiED\nu0gkgtHRUdTX1wttB51OB5/Ph/n5eTQ3N2N6enrblSpapUZ6fmTIFA6HRUMmOn6JRPoQkqZvr/V6\nPXQ6HdLptODAtrKyArvdnvegSsn3bidQ7rQQ+vPebmc6GAwKF+Bio7VYlhW938FgsGR3tTuFiiTk\n7aTO+cicMx9TCULOJZTwer2w2WxobGyUtcKmZIVMFHXkvSADO7fbLURfkYpYp9MhFovBZrMJFXMx\ng0V6yJTpRUwq6ZWVFUSjUZHTGiFp+gJRX18Ps9ksRDil02nBN4IeVNFy5Lq6uqI/29001FMKcoZ6\nmSnVBHS0ltPpFO1M0+t42S7wmQXSbo9vAiqMkLdbe6PTnOXKnAmUTA2ROmnJRofFYsHk5CRqampk\nPZaSsmRC7jqdTvBJ7u7uFgy5SUXCMAwWFhYQi8UwMjJS0sgck8kkGb9EV9Jk8EfaHH6/H2azGSaT\nSThmQtydnZ3C6ySDqo2NDSwsLIhSqwlJ57NrrhFyfpCK1iI703KitTILJIZhcgbw7gZUFCED71pw\n0uQpJXPOtxpSOsaJIBaLYW5uDul0WtjoKBd0Oh28Xi8cDgeam5u3DOw4jsPy8jLcbjcGBgbyDotV\nCgaDYcsmQDKZxNzcHDweD+rq6sAwDP74xz+K/B5Iq4K0eEhF3t7eLrwO0jYJBoNYWVkRtgnoDY9s\nAgqNkItHtp1pMtSNRCJCtBb53czMDDY3Nwt+D7aLb/rWt76FJ554AkajEW1tbfje976Hvr6+wl9k\nDlQcIQPvrr3xPI/19XU4HI6805ylHlNJQuY4DhcvXkQoFMLIyIjITKYcIBXn2tqa5MBufX0dS0tL\n2LNnj2hgV26Q1JPl5WX09fXhwIEDIoe5fIeHROlGD3ZJX3o7UUu55OGVRMhSyBat5fF44PF4sLm5\niccffxxOpxPT09MYHx/Hpz/9acEuMxfkxDcdPnwYb7zxBmpqavD444/jS1/6En7yk5+U4qVWHiHr\ndDro9Xokk0mcPn264DTnTChFyKR/HYvFMDQ0VLYqkyCZTMJmsyEWi8FqtWJ0dFQgKJ1OB7/fj/n5\neTQ0NMga2O0kyDCRrN9JGZXnGh6GQiFZw0O9Xi/4GZPhIcuyWzwjyC007Wm8E6KJchNyOS/OVqsV\nBw4cwPT0ND760Y/i1VdfxezsrOx4NDnxTR/4wAeE/338+HH88Ic/VPZFUKg4Qvb5fJidnUU6ncax\nY8cUU7AVu81AbvdJ/7q+vr6gVGelQC4MHo8HQ0NDaGtrw8WLFzEzM4PGxkaYzWZ4PB4YDAYcOHBA\ndk97JxCPxzE3NwcAeR9btuEh8SEmAbRSw0NSBROiptsgOp0OLpcL0WgUOp0O6+vriEQiYFm25KKW\nchNyubyZpZzeTCZTXopaOfFNNJ588klcf/31hR/0Nqg4Qk4mk9i/fz/eeustReXEhVbIpG2ysLCA\njo4OYYXN4/Eo+mWW27+kg017enpEA7vR0VGEQiE4HA5EIhGYzWbwPI+5uTlB9lxOAQYRyPh8PkGp\nqBSy+RATkl5fX0c0GhWGh2R9jpB0Op1GMBhEbW0t2tvbRaIWUpFnE7WQIVWh72k5Cbmcz00biBW6\ng7xdfBONH/7wh3jjjTfwm9/8Ju/nkYuKI+Q9e/YIk3WlPZHztczc3NyEzWZDfX39Fnc4JfeHpZKe\npUBW6qQGdjzPY3l5GS6XC/39/ejs7BSp5IjB0vLystCHJQRdX19fUgEG2cteWlpCb28vjh07tiMX\nBKnhIdkCCIVCAsEmk0mkUinhrocQBXlvCdnLEbXQykO5UmSO48pWpZabkMndUSgUKmggvl18E8FL\nL72E++67D7/5zW+Kbn/mQsURMgFZfVPqzcunQg6FQoIp0aFDhyRvqUsRdJrtxAiHw5idnYXJZMLk\n5KTQIyaE4fF4sLi4iI6ODhw7dkzU96QHKiRqiFR9oVBIIJREIgGz2SwiaSUsNgOBgOiiVu4eNr0F\nEI/HMTs7C7PZjN7eXqHtQQiWVMGktUHuOIB3c/aampqE6DEyPIxEIvB6vYjFYsLz5RK1lJMUgfJl\nJirhhbxdfBNwOXH67/7u7/Diiy+K9qhLgYoj5Mxd5J0kZCKUYBgGo6OjOfdzlfZElhKHkDy9WCyG\n0dFR1NXViQQUwWAQNpsNVqsVU1NTsnubdB+W/oJmDstisZgw4CJEXVtbK+sETiQSsNlsSKfTmJiY\nkD2k2QkQwYzX693SOpE7PCTvRXV1tcgHhKjcyPCQ4zhJ9zVaPJFOp1Wz9bKTyOwhlyq+6e///u8R\niURw4403AgD27t2Ln//854q+FuF4SvKoKoDSa2q5Ho9hGNjtdgQCAYyMjIjsDLOhFBUyAW1GNDw8\njNbWVpHCjhA1y7LYt28frFarIsdhsVhgsVhErz+VSgmV9MbGBmKxGAwGg1Bl1tfXi3yQWZbF0tIS\nPB6PcOxqAc/zgoikq6sLR48ezUqEcoaHZDtju+Eh2fro6uoSSJoWtWxsbMDtdmN9fb1gUctuhBKE\nDGwf3/TSSy8VfpB5ouIIuVQm9VIEmk6nsbS0hPX1dQwMDGB8fFz27VspKuRcAztyvIuLiwgEAhge\nHlZ0KJYNJpNpy7AsnU4LJL20tIRIJAK9Xg+DwYBoNIqOjg5MT0+rJlkZeNeD2mQy4fDhwwWTndTw\nkHgQk+FhJBIRVcFSJJ15h9Lc3Cz4TxciaikE5Y7lyiRksrq2m6Geb7zCKGWFzHEcVlZWsLy8jJ6e\nnrwk2KU4PoPBgM3NTVy4cCHrwG51dRWrq6vo6+vDyMhIWXefjUajSDIbCoUEF7ne3l7EYjG8+eab\nAC6vlpF2R11d3Y6bobMsK2x2jI6OlsQrIZsHcebwMFPaXVtbi0gkAp/Ph87OTkFmXoiopZCWR7l7\n17QopRKc3oAKJmSlK2Si/iMrbG1tbTld2LaDUhVyOByG1+tFPB6XHNh5vV4sLCwIgaxqSt8gLnKJ\nRALj4+MiuSwAkfhibW1N8KwgZEKIuhSVNM/z8Hg8WFhYEBRcO3kRyyYhJsrDzc1NzMzMgGVZ1NfX\nw+v1Sg4PiailsbFRGB5KiVoKSWopd1oI8O4dcSUYCwEVSMh0yyIWiyn2uH6/H7FYDD6fTxHln8Fg\nKOqCQYZe8XgcTU1N6OrqElatdDodQqEQbDYbampqirrFLgWIzzOJfcomkDEYDIJCjv5vSeXodrth\ns9mERBF6w6OYbYxIJIK5uTlUVVXhyJEjqkknIa5pPp8PgUAA+/btQ1tbm2jjhR4e0hUwPTwkyS1W\nq1WkPCR9abmiFjUQMgERhux2VBwhExQb40QQDocxNzcHg8GAqqoqkaSyGBgMBsEgJR9IDeycTqdg\n2VldXQ2fzyfET2VWneUEPRTr6OgQWiv5IFvlSEia3BGQHVWapLcj1nQ6jYWFBQQCAYyNjanuBCet\nncbGRtHdjtTGCz08JKkeZHhIdpwz+9LV1dWoqqpCR0eHMJegRS1LS0tIpVLCOh+56JXDVCmzf10J\nXshABRKyUkO9eDwOm82GZDIprLApnauXT8uCdqzr7e0VDey6u7vR0tKC+fl5eDweVFdXg2VZXLp0\nSURIclfOSgFSdVosFsUrdjo2iIC+vSfRUQzDCD1YWnUIXDZPWlxcRG9vb9l77JlIp9NCoo3crZhc\nw8NQKCRUwQCEKthqtQo78yzLguf5nKIWr9eLcDiMc+fOFSxqKRSZ1blGyCqGTqcreGhGvH79fr9Q\ngZIvlpLqP7mEzPO8oLBrbW0VnNbo28+1tTXBI+PAgQPC8dHm7sS/mCQ67BRJk/czHA5vu5utJLIZ\nC2WqDuPxuEDUfX19qjqp6T723r17ZQUp5EK24SHpJxOCTafTqK2tFa3imUymLUktwOXv8fDwsGhT\nJB9RS6HI9EJOJBJlS15XEhVJyED+FTIdZz8wMICxsbEtX/7tFHH5QE7yNLlFtVgsuOqqq7YM7DY3\nN2G329Ha2irpdiZl7p6NpAlBK0XSZBNldXUV/f39ku/nToNWHTY3N8Nutws+1DzPC33pTNXhTlR8\nmYjFYrh06RIsFktJ+9h6vV4yconcXZALFzF/J+SaTqfhdDrR398vWBWQ7xAZHuYStdBbHoX0+2lC\nJu2Lcn+/lEDFErLcCpluBXR3d+dcYSOPqYR8N1eFTAZ2iUQCo6OjsFqtIoUdfft/1VVX5ZWSsBMk\nvbm5ifn5eaGiV8vgBxD7YvT19YkuFLSAIzPfrxjVYT7gOA6Li4vY2NjA6OhoWSr2XLalfr9f8O8w\nmUxYW1tDOByWHB4C24ta6KSWzATxXJCKU9MIWaUgnsi5FtfJ7SCpMOWY1yu5OywVdLqdwo4oAhOJ\nBEZGRhRLF1GKpIl4QomMvVKA3HHI8cVQQnWYL3w+H+bm5tDZ2ZlTBVguBINBLC4uisynChkeVlVV\nwWKxoL29XThP801qoQlZqaxLNaAyXkWe8Pv9mJubQ21tLaampmRXmEqLOQgh00KTzIEd+f3S0hI2\nNjYwNDQk6muXCvmQdG1tLRKJBJLJJMbGxnZEAZgPUqmU4Okhte8sF9lUh+Q9oVWHtKBlu51eEhKQ\nTqcxOTmpugtZIpHApUuXYDQat7RP8hkeZvalgXe/4/mKWjiOE8RPwWCwrNFnSkKXp/yxvFpJmUin\n02BZFq+99hpOnDghfMBkhU2v12NkZCRvD4e5uTk0NTWJbm0LBcdxOH36NIaHhzE/P4+2tjb09/dv\nMQpyuVxwOp3o6elBd3e3qqomnufhdDrhdDqF3mGpetKFHh+JdxoYGBA8iksN2kc5FAoJhJSpOtTr\n9VhZWcHKygqGhoZK7iSWL3ieF+YAxcaM0cND8o/sOWcbHpJjILvStKhlZWUFyWQSv/zlL/Gzn/0M\nPM/jzjvvxOHDh3HkyBFZRdZ2WXrJZBK33HIL3nzzTbS0tOAnP/kJ+vv7C30LZH3xKrpCJlUoqZDi\n8XhR8lcl/ScikQii0SjcbjcOHz4siDpIe4LEEzU3N6vCdjITfr8fNpsNTU1NOHHihOiWka6kHQ4H\notHojpN0MBjE7Oxs1ninUkLKRzlTdRgMBhGPx1FVVYWuri4hmFctt96RSAQzMzNoaGgoaF88E7mG\nh6FQCIFAQHJ4WFNTI8o8BCCs57W1teHuu+/Gtddei8ceewxmsxk/+tGP0NTUtK1eQE6W3pNPPomm\npibMz8/jmWeewT/8wz+ULEuPQB2fvsIgJ7per8fs7CxCoZAQU1Ts2lCxLQt6v9lisWBiYkI0sItG\no7DZbKrtw5Lj5zgua3zSdu2OUpI0wzDC+7t//37V2HYS1WFtbS2i0ShMJhMmJiag1+uF7Q7iwKek\n6jBfcBwHh8OBzc1NjI+Pl7QVQA8Pu7q6ALy7mkiqaJIwbTabRf4dfr8fnZ2dYFkWv/71r5FIJPCZ\nz3xG9nPLydJ7/vnnce+99wIA/uIv/gK33357yUUwFUnIxBAmGAyiubkZExMTiryJxaj/iAqMeOi2\ntLTg1VdfhcPhEBI3lpaWEIvFMDIyojqVGHlPNzc3MTw8nPfta6EkXVNTI6tNQ6/ZDQ4Oor29XVVT\n91w7xVKqQ7LP63A4kEqlBPEGeV9KsQYXCARw6dIldHZ2Ynp6uiztMXo1MVN5uLGxAYfDIbR67rzz\nThiNRhiNRvzjP/6jUOTIgZwsPfpvSKDA5uZmSS1hK5KQfT4fAKCrqwvNzc2KnZhGoxHxeDyv/4Ye\n2O3du1c0sJucnITf74fdbkckEoHJZEJjYyMCgQA4jiuZcU4+IJmAi4uLwq2dUieqFEkTa858SJoM\naVtaWlS3Zgdc3ikmySLb7RTTqkO6aozH4wiFQvD7/YLqsKqqSnhPisk6pFt62RJuygkifnK73Th0\n6BDq6urw05/+FCaTCTfeeCPq6urwwgsvIBKJ4NSpU7IfMxOZ710+eXtKoSIJub29HU1NTbDZbIo7\nvsltWRDfBjKwy1TYARCW7vfs2YPp6WkAl0/eYDCI9fV1oTVAhkHkxNspwgkGg0LA6U71sTOtOQGx\nfzJN0jU1NYJT2YEDB1TTniBQaqeYmArV1NTkVB0WknVIVj/7+vry8vPeKZBeNrGV3djYwG233Yba\n2lq89NJLBVercrL0yN/09PQIIbal3iCqSEImIIMSpSCXkAmRVVVVSQ7s/H4/5ufn0dDQsIXoSHXU\n3d0N4N3pdCgUEpbwyd/V19ejoaFBcO1SCiRRhGEYRRNFCkUmSROiW1tbQ1NTEziOw4ULFwpud5QC\nZKeYmCgpfRy5sg5JG8jlciEejwv7vISoa2pqkEwmBQ9qNTnaEdCrnvv27UNtbS2ee+45PPjgg/ja\n1zDIxCYAACAASURBVL6GP//zPy/q4iEnS+/kyZP4wQ9+gBMnTuC5557DBz/4Qa1CLgTkTVPK8Y1g\nO7kzGXiRTL26ujoREZPMPVLRybk1lJpOk4l9KBTC8vIyIpGISLpaqECBZVk4nU643e4d23fOF4To\n2tvbceLECdHdQrZKOlNdV0qSZhgGc3NzSKVSO75TnCs2irbnDAaDSKfTwkWOYRiYTCbVfNbhcBgz\nMzNobW3F9PQ0PB4PPvOZz6C+vh6vvPJKUet3BHKy9P72b/8WN998s5Cu88wzzyjw6nKjIveQeZ4H\nwzBwu90Ih8MYHh5W5HETiQTeeecdHDlyRPTzVCqFhYUF+Hw+QWFH3LJ0Op3we3IspZDDkt1X0n+N\nRCKytxjogVNXVxf27t2rqn1n4PJ7Pzc3B57nMTo6KpvoaOFGKBQqGUnTO89q3CkGLispZ2ZmUFdX\nh76+PmHlLBQKCQq7TOOpnfwekDsfr9crVMXPPvssvv3tb+O+++7Dn/3Zn6nmolEAtD3kUgedchyH\n5eVlrKysCNFIPM8Lf0OEEy6Xq+QGO1K7r/SALNPtjfxjWRY2mw3V1dWqvnV1u90FhZ4ajUY0NTWJ\nLoI0SS8uLgpBo3SfPh8yCofDuHTpkrCzW+5BbCZooqN9nquqqrb06sm6WaGqw0JBquK2tjZMT0/D\n7XbjU5/6FJqbm/Gb3/xGderPUqEiK2QAQorC0tISDh48qMhj8jyP119/HSdOnBCGIe3t7ejr69sy\nsPN4PHA4HOjs7MTevXtVM/knfgx+vx8ulwupVAq1tbVobm4WCGm7QdBOwev1Yn5+Hh0dHcJ7XCpI\nVdLbkTTxKQ6FQkVJskuJQCCA2dlZ0fc0H+RSHdJ3GYV+vzmOE+xu9+3bh5qaGjz99NN45JFH8PWv\nfx0f/ehHVfFdVACyXkTFEjLDMIhGo5idncXhw4cVe9zf/e53sFgsqKmpwdDQkDCwAy738ILBIGw2\nG6xWq/B7NYFU9Wtra8K+LiFp8i8ejwvTevKv0JWqQhCPxzE7Owu9Xo/R0dG83OyURDaSJsS7ubmJ\n/v5+dHd3q4400uk05ufnEY1GBaJTCvSgmZA0nXVIiHq7O4VQKISZmRnhgutyuXDXXXeho6MDDz74\noKq8qRWARsgMw+D8+fM4duxY0Y8Xj8cxNzeHjY0NHDt2TDA4IX3ieDwuKK0K8ckoNYjRPV3V56pq\nyEoV+ZdMJkV7r4SklQTxpCbiGTXeppJba57nUVVVhXg8XvbeaybIumVfX59ge1lqkKxDcgGjvSoy\nVYccx8FutyMYDGLfvn2orq7Gj370Izz22GO4//77cf3116vuAqcAruwecjGpITTogd3IyAgikYgQ\nkUTMThwOBwKBgDCNVRuIf7LZbJbtn5yZ00bvvQYCATidTjAMI5xw5KQr5I6Aztrr6upSpfUk6WV7\nPJ4tO8XZHN92mqTJKhsATE1N7WiwLZ11SPZ5pVSHyWQSDMPAarXC4XBAp9Ph61//Orq7u/Hb3/62\nIpKji0HFVsipVAocx+G1117Dn/zJn+T932cO7Pbs2QOe52G32+HxeGA0GoXKuLe3V5WbCalUSshi\nGxkZUfzLTivIyD/Sk6ZJOpeghKjYTCYTRkZGVJWOTUDvFMvtw9IkHQ6HBZIudjVRCkTJ5nQ6MTw8\nrIgbodJgWVaoigcHBxEIBPDlL38Z58+fh8lkwoEDB3DjjTfilltuKfehlgpXdsuiUEKmjeuzDezI\nsIkMwMLhMBKJhOiWvqGhoWz9Y5KCQi4mO3XbCogdvMg/lmVFJE12qh0OB3w+X1EOfKUEvVM8Pj5e\n9E5xKUg6Go3i0qVLqK2txfDwsOo2PIB3PTL27NmD3t5erK6u4o477kB/fz+++c1vor6+HsvLywiF\nQjhw4IAiz3nrrbfiF7/4Bdrb2/H2229v+T3P87jrrrvwwgsvoKamBk899RSmpqYUee4suLIJOZsn\nci4EAgHMzc2hpqYGw8PDIm9WnU6HUCgEm80mDPToao6opILBoEBE5Ja+oaFBOOFKfcKQ+KSWlhb0\n9/er4gQl/UXyvvh8PiQSCcGvYacl4dshc6e4WJfAXCiUpOkWyvj4uOrMqIDLVfH8/DwikQgmJiZg\nsVjwgx/8AP/xH/+BBx54AH/6p39asvf1t7/9LaxWK2655RZJQn7hhRfwyCOP4IUXXsCZM2dw1113\nbTEXUhhXdg+ZQE4OXiwWw9zcHFiWxfj4uDCwIwo7IiVOpVIYGxuTXG+iVVK0lJVUi6TqZll2i+xZ\nCSIir0Gn06nOtpMQjE6ng9vtRlNTEwYHBwUFGS0Jp4lIaUm4HOz0TnGuPenMfWDy3uj1eiwtLZVM\nlq0E/H4/Zmdn0dPTg9HRUSwvL+OOO+7A0NAQfv/735d8RfCaa67B4uJi1t8///zzuOWWW6DT6XD8\n+HEEAgG4XC7B0KlcqFhCJlde4mchRcikx+r3+wVLTJZlhYEd8Yb1+XyClDjfY8j0e6VXhlZXVxEO\nh7fInq1Wq+zKgeTw+f1+1Q4VifVoIBDYIkzIJQkPh8Ml67tKHaNadoqzkXQgEIDD4RACVzc2NpBI\nJEr+3uQDsm4Xi8UwOTkJi8WC733ve3jiiSfw4IMP4kMf+pAqNiik7DdXV1c1Qi41iJ8FXTFyHAen\n04nV1VX09fVhdHR0i8KO9GB7e3sVrUKyeVOQgEeiHCOKOtLuoAMeyTGSQU5vby+Gh4dV8UWnQVt3\n9vb2YmRkJOcxEhN3+vabloRLqeqKNbanZeO9vb0in2I1gRhS9fb2CnvP21XSO03SZPjZ29uLsbEx\nOJ1O3H777RgfH8err76qqlXQclhrykHFEjJdIRODIZ7n4Xa7sbCwgI6ODlx99dXC6hrB5uamkES9\nUzJYKdkzLdZYX18XiTUMBgPcbrdq452Ay7f+s7OzqK2tLeoYs703dNhqNBqF0WgUDVQzL2BSiMfj\nuHTpkiyf4nKBYRjMzs6C5/ktq2xSlTR9Adspkk6n07DZbEgkErjqqqtgNpvxxBNP4Pvf/z6+/e1v\n49prr1UF2dGQY79ZDlQsIROQHjKRkFqtVkxNTQkDOyLsILu6FotF9q5uKSGVcEyGjgzDwGKxYHNz\nE/F4XFQtlpuc6VW7sbGxkkQAZUsfIRcwj8eDWCwGs9ksem+IJJz2xxgdHVVlm4fnebhcLiwtLeVl\nVpQtz69UJL25uQmbzYa9e/difHwci4uLuOOOO7B//368+uqrqvOoJjh58iQeffRR3HTTTThz5gwa\nGhrK3q4AKnjLguM4pFIp2Gw2eL1emM1mjI6Oora2VkTEyWQSdrsdiUQCIyMjqowTJwq2jY0NkcEO\n2QOmNzvIihlpdezU9gJNIDu9apcNxM+E/EskEtDr9UgkEmhubsbQ0JCsSnqnEYvFMDMzU9JVNil3\nwHxImpxbDMNgfHwcJpMJTzzxBP7rv/4LDz30EK655pqyvq+f+MQn8Morr8Dr9aKjowNf/epXhTvl\nz372s+B5HrfffjtefPFF1NTU4Pvf/74QElEiXNlrb6lUChcvXoTH40FLSwsmJiZElph0moNavX/p\nHmx3dzd6enq2rWLoFbNgMLhle4EEbSp5yxoKhTA7O4v6+noMDg6WvUqXAtkpTiQS6OzsFFSHmZLw\ncu+PEz/qsbGxHd/NlkvSPp8PNpsN/f396OzshMPhwB133IHJyUncd999qq2Ky4wrm5DT6TRWVlag\n0+ng9XoxPDwsVIoulwtOpxM9PT3o7u4u+2RaCnR80uDgYFEkkc0rmb6dr6mpyfuCRLLYYrGYYMiv\nNmy3U0xLwjP3x3eyFRQKhXDp0iW0traiv79fNd9J+rsTCATg8/kEX5Tl5WUEAgG89NJLePjhh3HN\nNdeU+3DVjCubkEmFWFtbK5jDp9NppFIpoZJraGhQXVWcTCYxPz+PRCJRUpIjgzHS7qB7rqTdka2P\nTpPcwMAAOjo6VPc+AuKd4sHBQdm3/nIk4UqJfIh4IhwOC6bsagQxLBoYGEBraytefvllPPTQQ/D7\n/TAajTCbzbj//vtx7bXXlvtQ1Yorm5DPnj2LL37xiwgGg+jq6oLb7cb111+PU6dOgeM4wU6RxHsT\nIirXMI8eNA0ODpZUHZYNdM81GAwimUyiurpadDtPvCeampowMDCgCiVgJsjeczAYVGynmBjlkPeH\nuJllknQ+/XoyECN3amq8qNFbHuPj4zAYDHj88cfxzDPP4OGHH8Z73/teAJcl3DzPK7ra9uKLL+Ku\nu+4Cy7L41Kc+hbvvvlv0e6fTiU9+8pMIBAJgWRb3338/brjhBsWeX2Fc2YRMcOedd+LVV1/Fdddd\nB4/Hgz/84Q/Q6/U4fPgwpqamMDU1hfb2dqFaJCREk3QpSYd2OiNm9mq5XSW388FgEH6/Hx6PByzL\noqmpSWRorybJ88bGBux2u2hft1TIlISHw2FwHLdF5JP5/pB+NsuyGBsbK/tGTzaQFdHBwUF0dHRg\nbm4Od955J44dO4avfe1rJVWDsiyL0dFR/OpXv0JPTw+OHj2Kp59+GhMTE8LffOYzn8Hhw4dx2223\n4eLFi7jhhhtyqvPKDE06DVz+0L7zne8IJybP84hEInjzzTdx+vRp/Ou//ivm5uaEQMUjR47gqquu\ngslkgtfrxcLCgiB3JqIFpSS94XBYlE6tNqcznU4Hi8UChmEQCAQwPj6OtrY2QQ6+vr4Om80mIqFS\npGDLAdkpNplMO7ZTTFtOSqWEEyUmAOHvSNbj8PCwKnP3gMsXjEuXLkGn0+HIkSPQ6/V4+OGH8T//\n8z945JFHCnJPzBdnz57F8PAwBgcHAQA33XQTnn/+eREhE38Z4PLMRQ17xMWi4itkOSArW2fPnsXp\n06dx7tw5eDweDA8P48iRIzhy5AjGxsaEXVci6aWr6HzWpxiGgd1uRzQaxejoqCpX7YDL6jCbzYbm\n5mYMDAxkrYQ5jhMNDcn7o5SaLhd2w04xy7KC6RPP89Dr9ZIBtGq4MyJVMdl9vnTpEu6880685z3v\nwVe/+tUdq+afe+45vPjii3jiiScAAP/93/+NM2fO4NFHHxX+xuVy4brrroPf70c0GsVLL720JYBY\nRdAqZLnQ6XTYs2cPPv7xj+PjH/84gMsn0ezsLM6cOYPnn38e58+fRyqVwqFDh3DkyBFMTU2hqakJ\n0WhUECJYLBYRSWdWaRzHYWVlBaurqxgYGMD4+Lgq+4bJZBI2mw2pVAoHDhzYNv6HXJxoyXNmwGo0\nGoXJZBI53xWb3ef3+zE3N4f29nYcO3ZMFYSWCeKrvb6+jn379gmqOqmgVbkp4aVAMpnEpUuXYDAY\nMD09DZ1Oh4ceegg/+9nP8Nhjj+Hqq6/ekeMgkCNtfvrpp3Hq1Cl88YtfxOuvv46bb74Zb7/9tiq/\nB3KhEXIWGAwGTExMYGJiAn/zN38D4PLC/vnz53H27Fk8+uijuHjxIurq6nDkyBFMT0/jqquuQm1t\nrShRg4g0SBVOyEMtfVcahDxcLpewIlYojEbjFjUdcXcLBoNYW1sr2EOaYRhBlKA2ZzsaZJWtpaVl\nywVDSvYsJQk3mUyiOw2lhSz0rjsxt5+ZmcEdd9yB97///fj9739flh63HGnzk08+iRdffBEAcOLE\nCSQSCXi9XtW2guRAa1kUAZ7nsbm5ibNnz+LMmTM4e/YslpeXsXfvXhw9ehRTU1OoqqrCzMwMJiYm\nhJSRuro6oaLcySooF8iyf1tb27Z5e0oh2w5wtvUyet2OBLSq4b3LBMuygrvdvn37ito8YBhG1A7K\nJQnPF8lkEjMzMzCZTIKp0ne+8x38/Oc/x3e/+10cPXq04OMuFul0GqOjo3j55ZfR3d2No0eP4sc/\n/jH2798v/M3111+Pv/qrv8KpU6cwMzODD33oQ1hdXVXldwLalkV5QGLNf/3rX+Pxxx/H2toahoeH\nsXfvXmFoODAwIAzGpFbvdjLhOZFIYG5uDjzPY3R0tOzVJvGQpuXgHMfBYrEgGo2ivr4e4+PjqjQC\nAt5dZSPKylJ8jlKS8MyU8FxVLS1zHxkZQWtrKy5evIg77rgDH/zgB3HPPfeoYsD8wgsv4POf/zxY\nlsWtt96Kr3zlK7jnnnswPT2NkydP4uLFi/j0pz+NSCQCnU6Hb37zm7juuuvKfdjZoBFyOfHv//7v\nqK6uxs0334x0Oo233noLZ86cwZkzZ3DhwgWYzWYcPnwY09PTOHz4MFpaWoRKKJFIiFbvSqEUo4dh\ntD+G2kB8in0+H9ra2pBMJoUTMFMOXs7KiKyypdNpjI+P7+htPkmryZYSTr5HZrMZiUQCMzMzsFgs\ngu3sQw89hP/93//Fd7/73VL7OVzJ0AhZreB5HqFQCOfOnRNaHXa7HV1dXUI/+uDBgzAajSLTILJ6\nR0yDCh1ekEzAfEI7dxrb7RSzLCsiIDIUoy9iO2EcRPdg1dRGkWoHRaNRsCwLs9kMj8eDtrY23HPP\nPbjuuuvwla98RRVVcQVDI+TdBJ7nsbKygtOnT+Ps2bM4e/asEABKWh0jIyPCYIxeLSP96O0IKB6P\nY3Z2Fnq9HqOjo6oVJJDjNBqNGB0dld2eIGuJpN1Be0gTolaSdMjus8ViwcjIiCpNlYDLxzkzM4Pq\n6mp0d3fj0qVLuO+++3DhwgU0NDTg8OHD+Ou//mt87GMfK/ehVjI0Qt7tSKfTmJmZEXajz58/D57n\nMTk5ienpaUxNTaGzs1MQImRbvSP2nV6vFyMjI6rc1QVKs1NMV4nBYFAR4yCe5+F0OuFyuVS7+wyI\nk2/IcV64cAF33nknrr/+enz5y1+GwWDA3NwcOI4TDcyUwHbSZwB49tlnce+990Kn02FycvL/t3eu\nMVFeXRu+HkBEsBysREQUAZFjrSJ8osZUG5XWA2r6qW2irV9T0xhNlUZEkyrEWrUpJe2roaaKqW9t\nAfWPpFVsrMeqnBRUEAUUrGO1iOLIIYAM+/uB83RQhKcIw4D7SogM7LCX6KzZs/a97sXPP//cqTFY\nEDIh9zaMfgrnz59XT9FGXwlT6V2/fv3UBFRbW6u2Ow8bNgwnJyeLlNyZaoq7soxiOnj26XKQMUG3\n5SFtNCwyenlY4u8Smk/FV65coX///owYMQKDwUB8fDxHjx5lx44djB49ukv319L6XFxczIIFCzh2\n7BguLi6Ul5f3aMlaO8iE/DJgnAlnvDDMzs7mzp07uLq6cv/+faZMmcKyZcuwsrJSSx1CCIuR3plq\niv38/NptQukKjJ4UxlJHax7SdnZ2lJWVqS3klmg1Cv+Uvm7fvo2/vz/Ozs5cvHiRlStXMmvWLNau\nXWsWhcq5c+eIi4vjyJEjAGzZsgWAdevWqWvWrFnDyJEj+eijj7o8HgtAduq9DCiKwqBBg4iMjCQy\nMhKArVu3kpKSwuzZs6msrGTJkiXU1dURHBysnqSNnhSlpaXdIr0zHdLa3Zdhpp4URkx9gIuKinj4\n8CF9+/Zl4MCBqnF7RzykuxLjpJFXXnmFsLAwGhsb2bRpE8ePH2fXrl2MGjXKbLG0NtU5MzOzxZqi\noiIAJk6ciMFgIC4ujrfeestsMVoiPS4ht1eXqq+v5/333+f8+fO8+uqrpKamMnz48O4JtpuYNWsW\n0dHRLd5O19fXk5eXR0ZGBjt27CA/Px97e3tCQkJU6Z2TkxNVVVVqF11XSe+MA1CNicMSLTytra1x\ncHDg9u3bWFlZMWHChBaqF9O5fU+3g5sbIQS3bt3ir7/+IiAgACcnJ/Ly8li5ciVz587l1KlTZr9w\n1NL6bByOeuLECXQ6HZMmTSI/P9/sk1IsCct7JrSBwWBg+fLlLepSkZGRLepSSUlJuLi4UFJSQkpK\nCjExMaSmpnZj1OYnODj4ma/17duXcePGqZ4EQggqKyvJzs4mIyODAwcOUFZWhoeHB6GhoYSGhuLp\n6alOXHna9a4j0juDwcD169c71ae4KzBOJy8tLX3GgP/pwbNG/a9er0en07XwkDb+nroyGdbU1FBY\nWIiTk5N6Kt64cSOnT59m9+7dvPbaa122d1toaX328PAgPDycPn364OXlhZ+fH8XFxd3aIdjd9Kga\nspa6VEREBHFxcYwfP57Gxkbc3Ny4d++eRb21tFSMKgfTenRVVRX+/v5qqcPHx0dVLvwb6V15ebnZ\nfIpfBGPjhHEobkcUGKbTRvR6/Qsb2T9vH6Mixd/fHycnJy5cuMCqVat45513WL16dbfK8LS0Pqen\np5OcnMyePXuoqKhgzJgx5OXltXjB60X0vhqylrqU6RpjXfT+/fsW24lmSVhZWeHl5YWXlxfvvvsu\n0Kztzc/PJyMjgx9//JFLly5hbW2tGvyPHTuWgQMHUlVV1arrXd++fblx4wY2NjaEhIRYbPOB6dt+\nX1/fDicFRVGwt7fH3t4eNzc39WcbLw2f9pA2/p7+jYd0dXU1hYWFuLi4EBYWRkNDA7GxsZw7d449\ne/Z0unytI9jY2LB9+3YiIiLU1uegoKAWrc8RERH89ttvBAYGYm1tzVdffdVbk7FmetQJef/+/Rw5\ncqSFR2pWVhbbtm1T1wQFBXHkyBE8PDwA8PHxISsr66X/h+4shBBUVVWpBv/Z2dmqKZGp9E5RFPLy\n8rC3t39mVp8lTRmBfxKcs7Mz3t7eZonNYDCo+nG9Xq9eFJo63z19aWh8B3Pv3j38/f1xdHQkJyeH\nqKgoFi5cyKeffmqR9XgJ0BtPyFrrUrdu3cLDw4PGxkb0er3FCvd7Ioqi4OjoyJQpU5gyZQrwj2LC\naPD/5ZdfcvPmTUJCQpgyZQpjxozBzc2Nx48fqydES5DeGQwGSktLefDgAQEBAWataRvbvJ2cnNR3\ndKYe0sZp3kYPaVtbW+7cuaNOtmloaGDDhg1kZmayd+9eAgICzBa7pOvoUQk5LCyM4uJiSktLGTJk\nCCkpKc909kRGRrJnzx7Gjx/PgQMHePPNNzvlid6euiMhIYFdu3ZhY2ODq6sru3fvxtPT84X37Qko\nisKQIUOYN28e/fr1Iz8/n9TUVBoaGlSD/9jYWAwGg2rwP3bsWAYNGkR1dTWlpaVUV1e3MLDvauld\nZWUl165dY/DgwYSGhlqEn0drHtJ1dXUUFxdz584d+vXrx7fffsvly5cpLy9n6tSp7N+/Xy2NSHo+\nPapkAe1b8tXV1bF48WJyc3MZMGAAKSkp6lyujqKl6+j48eOMGzcOe3t7vvvuO06cOPHSqTug+W21\noiitJtLa2louXLigdhkWFhbi6OioljrGjBmDg4NDl7rePX78mOLiYurq6ggICOh2u9G2qKqqorCw\nUPWorq+v5/PPP+fSpUvMnj2bu3fvkpOTw+bNmzt9ooeWtmdoHrU0f/58srOzpVNc28hOvc5Ci7rD\nlNzcXFasWMGZM2fMFmNPRAhBRUVFC4N/nU6Hp6enaqgUHBysuuO9qOudcV7c8OHDcXNzs1ilR1NT\nU4tSSv/+/cnIyGD16tUsWrSIlStXdmmdW8sBBJpfMGbOnElDQwPbt2+XCbltel8NubvQou4wJSkp\nibffftscofVoFEXB1dWVmTNnMnPmTKA5GV2/fp3MzEyOHj3K1q1bqa2tJTAwUE3SgwcPpra2Fp1O\np0l6V1dXx7Vr17C2tjbbROqO8ujRIwoLCxk0aBBjx46lrq6OdevWkZeXR0pKCiNHjuzyGLRMfAZY\nv349a9asIT4+vstjelmQCVkDWrqOjOzdu5ecnBxOnjzZ1WH1SqysrPD19cXX15dFixYBzX4XFy9e\nJDMzk127dpGfn0/fvn1Vg/+QkBDV4N9Ueufo6Mjjx4+prKzEz8/PopU2xkkzlZWVBAcH4+DgwNmz\nZ4mOjuaDDz4gISHBbMoULQeQ3Nxcbt26xaxZs2RC7kRkQtaAFnUHwNGjR/niiy84efKkxepteyK2\ntraEhYURFhbGihUrEEKg1+tVg//Y2Fhu3LiBu7s7ISEhhIWFYWdnx+nTpwkPD8fW1paioiJ14Kyl\nSe/0ej1Xr17Fzc2N0NBQamtriYmJIT8/n3379uHr62vWeNo7gDQ1NREVFcUPP/xgxqheDmRC1oAW\ndUdubi4ff/wx6enpvdlC0CJQFAVnZ2emTZvGtGnTgH88iv/44w++/vprSkpK8PHx4fz582qpw93d\nnfr6eouR3hmHoer1eoKDg7G3t+fMmTPExMTw4Ycf8s0333TLi0Z7B5Cqqiry8/OZPHkyAHfv3iUy\nMpK0tDRZR35B5KWeRtpTd0ydOpXLly8zePBgAIYNG0ZaWlqn7C1vvLWzb98+bt68SVRUFAAFBQVq\nG3hubq5qhG5ajzaqOswpvXv48CFXr17F3d2doUOHUltbS1xcHFevXuX777/Hx8en0/fUipa2Z1Mm\nT55MfHz8S/t/TiNSZdEbkDfenYexhTknJ4esrCyys7O5du0aAwYMaCG9Mxr8d4X0zmAwUFJSQnV1\ntSq7O336NGvXrmXp0qWqd3V3094BxBSZkDUhE3JvQKvkbtWqVUydOpX4+Hj55PgXGJ3dTA2V7t69\ni7e3t5qkAwMDWwxV7aj0ztiMMmTIEDw8PKipqWHDhg2UlJSwc+dOvLy8zPA3lnQTUvbWG5A33l2L\noii4ubkxZ84c5syZAzRfWhUVFZGRkcEvv/zCxo0baWhoaGHw7+7uTk1NTavSO0dHxxY+FAaDgeLi\nYmpra3n99dexs7Pj1KlTrF27lmXLlpGYmGgRp2JJ9yMTsoUjb7zNj5WVFf7+/vj7+7NkyRKgWcts\nNPhPTEykoKAABwcH1eA/JCRENfg3ld7Z2try8OFD3N3d8fPzo7q6mqioKMrKyjh48OBLNzxB0jYy\nIVs48sbbMrCzsyM8PJzw8HCg+YXywYMHqsG/8TJx6NChhIaGEhQUxOHDh5k3bx5Dhw7lp59+IiUl\nhfr6eiZOnMhnn30m1TiSZ5A1ZAtH3nj3HJqamigrK2PXrl3s3LkTf39/9Ho9fn5+lJeXY2dnAnMI\nNAAABOlJREFUx5o1a9DpdGRlZREREfHMBVlnII2wLBJtUh0hxL/5kHQDv/76q/D19RXe3t5i06ZN\nQggh1q9fLw4ePPjM2jfeeENkZ2d36v6HDx8WI0eOFD4+PmLLli2trklNTRUBAQEiMDBQvPfee526\nf0+isbFRfPLJJ0Kn0wkhhGhoaBA5OTkiNjZWGAwGs+zv7e0trl+/Lurr68WoUaNEQUFBizXHjh0T\nNTU1QgghEhMTxYIFC7o8Lom2HCsTsqRNtDzBi4qKxOjRo8WDBw+EEEL8/fff3RGqRAhx9uxZMX36\ndPXx5s2bxebNm5+7/sKFC2LChAnmCO1lR1OOlVe7kjYxNZqxtbVVjWZM2blzJ8uXL8fFxQVA1ka7\nkdZUObdv337uemmEZVnIhCxpEy1P8KKiIoqKipg4cSLh4eGkp6ebO0zJE0QHjLCio6O7OiyJRqTK\nQtImWp7gjY2NFBcXc+LECXQ6HZMmTSI/Px9nZ2dzhSl5gjTC6tnIE7KkTbTOMZwzZw59+vTBy8sL\nPz8/iouLzR2qhJZGWA0NDaSkpDyj5DAaYaWlpcnykoUhE7KkTbQ8wefOncvx48cBqKiooKio6IXH\nZpmSnp6On58fI0aMYOvWrc98/88//1SHqY4aNYpDhw512t49DRsbG7Zv305ERAQBAQEsWLCAoKAg\nNmzYoJpdRUdHU11dzfz58xk9enSXSO8kHUTr7Z+QKouXlvZkd01NTSIqKkoEBASI4OBgkZyc3Gl7\na1F5LF26VCQmJgohhCgoKBCenp6dtr9E0kloyrGyhixplxkzZjBjxowWX9u4caP6uaIoJCQkkJCQ\n0Ol7axknpCgKjx49AprN3lurmUokPQGZkCUWjRZzpbi4OKZPn862bduoqanh6NGj5g5TIukUZA1Z\nYtEIDSqP5ORklixZgk6n49ChQyxevJimpiZzhSiRdBoyIUssGi0qj6SkJBYsWADA+PHjqauro6Ki\nwqxxdhbtXWDW19ezcOFCRowYwbhx4ygrKzN/kJIuQyZkiUWjReUxbNgwfv/9dwAKCwupq6vD1dW1\nO8J9IQwGA8uXL+fw4cNcuXKF5ORkrly50mJNUlISLi4ulJSUEBUVRUxMTDdFK+kK/q3bm0RidhRF\nmQF8A1gDu4UQXyiKshHIEUKkKYoSCOwE+tPsSLhGCPFbJ+29G5gFlAshglv5vgJ8C8wAaoElQogL\nHdxrPBAnhIh48ngdgBBii8maI0/WnFMUxQa4C7gK+UTuFchLPYnFI4Q4BBx66msbTD6/Akzsou1/\nALYD/33O998GfJ98jAO+e/JnRxgC3DJ5rGvlZ6lrhBCNiqLogVeBnlmjkbRAliwkkjYQQpwCHrSx\nZA7w3yda0wzAWVGUwR3crjXTiadPvlrWSHooMiFLJC9Ga6faIR38WTpgqMljD+Cv5615UrJwou0X\nDEkPQiZkieTF6MwTazbgqyiKl6IotsC7QNpTa9KAD558/r/AMVk/7j3IGrJE8mJoOdVq4klNeAVw\nhH8uMAtMLzCBJOBHRVFKaD4Zv/tC0UssCqmykEjaQVGU4cAvz1FZzARW0KyyGAf8RwjxP2YNUNJr\nkCdkiaQNFEVJBiYDAxVF0QGxQB8AIcQOmtUfM4ASmmVv/9c9kUp6A/KELJFIJBaCvNSTSCQSC0Em\nZIlEIrEQ/h+BpMUpVXf34QAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from mpl_toolkits.mplot3d import Axes3D # For plotting 3D graphs\n", "fig = plt.figure()\n", "ax = fig.gca(projection='3d')\n", "\n", "# This time you will have to predict a plane, not a line. A plane embedded in a 3D space \n", "# has 2 input dimensions and 1 output dimension...\n", "\n", "# Plot the predicted plane\n", "# ax.plot_surface(x1, x2, y, alpha = 0.2) # alpha is the color channel for transparency\n", "# etc..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## D. Multiple linear regression with interaction" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's work with a more difficult part of the dataset: **Mitochondrial DNA**.\n", "Try to apply the simple regression model of part B to the two features. What do you notice?" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Define X and y for Mitochondrial records (don't forget the usual transformations)\n", "\n", "# Plot the two features side by side\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Is there a better explanatory feature than the two we have worked with so far? Provided the genes are all about of the same length, the length of the genome multiplied by the percentage of coding genome should give us a better estimate of how many coding genes there are. We will thus try to use an interaction term that relates the two variables.\n", "\n", "Define again the matrix and the vector, and add to the matrix a new column given by the product of the two features. Linear regression for genome length and percentage of coding genome spendings with additional interaction term. Plot this new interaction feature against the output and make a simple regression model. Is it any better?" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Define an interaction feature and plot it against the PCG" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Obtain linear regression coefficients and evaluate model adequacy indicators" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [], "source": [ "# Train a linear regressor on the data you plotted, \n", "# and plot the usual statistics and the resulting graph\n" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "Do a residual/leverage analysis to check if there are outliers and outliers with high leverage. To obtain the leverage scores you will need to evaluate a matrix inverse and to retrieve the diagonal elements of a matrix, these operations can be done with the numpy commands linalg.inv and diag respectively." ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Calculate h\n", "\n", "# Calculate the residuals\n", "\n", "# Scatter plot h against the residuals\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Make a multivariate analysis like you did in Section C, but this time add to the feacture matrix a column containing the interaction term you calculated in this section. How does the fit look?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Confirm your analysis by repeating the residual/leverage analysis for the multivariate+interaction case" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] } ], "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.1" }, "toc": { "colors": { "hover_highlight": "#DAA520", "running_highlight": "#FF0000", "selected_highlight": "#FFD700" }, "moveMenuLeft": true, "nav_menu": { "height": "63px", "width": "252px" }, "navigate_menu": true, "number_sections": true, "sideBar": true, "threshold": 4, "toc_cell": false, "toc_section_display": "block", "toc_window_display": false, "widenNotebook": false } }, "nbformat": 4, "nbformat_minor": 1 }