{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Putting it all together - a case study\n",
"> A Summary of lecture \"Statistical Thinking in Python (Part 2)\", via datacamp\n",
"\n",
"- toc: true \n",
"- badges: true\n",
"- comments: true\n",
"- author: Chanseok Kang\n",
"- categories: [Python, Datacamp, Data_Science, Statistics]\n",
"- image: images/bootstrap-plot.png"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt\n",
"import seaborn as sns\n",
"\n",
"sns.set()"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def ecdf(data):\n",
" \"\"\"Compute ECDF for a one-dimensional array of measurements.\"\"\"\n",
" # Number of data points: n\n",
" n = len(data)\n",
"\n",
" # x-data for the ECDF: x\n",
" x = np.sort(data)\n",
"\n",
" # y-data for the ECDF: y\n",
" y = np.arange(1, n + 1) / n\n",
"\n",
" return x, y"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Finch beaks and the need for statistics\n",
"- Data source\n",
" - Peter and Rosemary Grant\n",
" - 40 Years of Evolution: Darwins's Finches on Dphne Major Island\n",
" - Princeton University Press, 2014\n",
" - Data acquired from Dryad Digital Repository [link](http://dx.doi.org/10.5061/dryad.g6g3h)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### EDA of beak depths of Darwin's finches\n",
"For your first foray into the Darwin finch data, you will study how the beak depth (the distance, top to bottom, of a closed beak) of the finch species Geospiza scandens has changed over time. The Grants have noticed some changes of beak geometry depending on the types of seeds available on the island, and they also noticed that there was some interbreeding with another major species on Daphne Major, Geospiza fortis. These effects can lead to changes in the species over time.\n",
"\n",
"In the next few problems, you will look at the beak depth of G. scandens on Daphne Major in 1975 and in 2012. To start with, let's plot all of the beak depth measurements in 1975 and 2012 in a bee swarm plot."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Create bee swarm plot\n",
"_ = sns.swarmplot(x='year', y='beak_depth', data=df)\n",
"\n",
"# Label the axes\n",
"_ = plt.xlabel('year')\n",
"_ = plt.ylabel('beak depth (mm)')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### ECDFs of beak depths\n",
"While bee swarm plots are useful, we found that ECDFs are often even better when doing EDA. Plot the ECDFs for the 1975 and 2012 beak depth measurements on the same plot."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"bd_1975 = np.array(df[(df['year'] == 1975) & (df['species'] == 'scandens')]['beak_depth'])\n",
"bd_2012 = np.array(df[(df['year'] == 2012) & (df['species'] == 'scandens')]['beak_depth'])"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEJCAYAAACOr7BbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dfVxUZd4/8M888CBJQjCAYJrFCmb4kLX5UFqugQ9Qae62ldl9u1FmrWWv9SG1zbVMc93sp3ZXulttm7aabRburmmZ5fpQ2u2uZQL3pIlgDAiIIAhzzrl+f4wMzjAzzDBzZs7A5/169YrDOTNzXZziy7m+1/W9dEIIASIiIif6UDeAiIi0iQGCiIhcYoAgIiKXGCCIiMglBggiInKJAYKIiFxSPUDU19cjNzcXpaWlbc4dO3YMkydPRk5ODhYuXAhJktRuDhERecmo5pv/5z//waJFi/DDDz+4PD9nzhw8//zzGDx4MBYsWIDNmzfjvvvu8/r9a2rOQ1G0sYwjIaE7qqrqQ90Mv4R7H8K9/QD7oBWdtQ96vQ7x8Zd5/R6qBojNmzfj2Wefxdy5c9ucKysrw4ULFzB48GAAwOTJk7F69WqfAoSiCM0ECACaaktHhXsfwr39APugFnNZLYpKapDROx7paT3avV6LffCWbDGjpvAErD36wpCc3uH3UTVALF261O25iooKmEwm+7HJZILFYlGzOUTURZnLavH7dw9DkhUYDXrMuXeIV0EiHMkWMxq2rUCDIgF6I2Jy53Y4SKgaIDxRFAU6nc5+LIRwOPZGQkL3QDfLLyZTbKib4Ldw70O4tx9gH9Sw+8iPkGQFQgCSrKC0qgHDB/fy+Jpg9uFCaREaTx5Ftz4DEN0rw6/3qik+YQsOQgEUCVG1JxB/3ZAOvVfIAkRKSgoqKyvtx2fOnEFSUpJP71FVVa+Zx0CTKRaVlXWhboZfwr0P4d5+gH1QjWILDgBs/1YUj230pg+yxQzpdCGMqZl+DeO0/MUPxYoafYRff/EDgNyjL6A3AhefIJp69LX3Ra/X+fSHdcgCRFpaGqKiovD1119j6NCh+PDDDzFq1KhQNYeIOrHzjVb71zqn44649Jd6s5+/1KXThYBitUUuRYJ0utCvAGFITkdM7lxE1Z5Ak585iKCvg8jPz8c333wDAFi5ciWWLVuGcePGoaGhAdOmTQt2c4hIw8xltfj7/h9gLqv1630yescjwqiHXgcYjXpk9I736/2k04WAfPGXumy1HXeQMTUT0EcAOj2gN9qO/WRITkf8yMl+BQcgSE8Qu3btsn+9fv16+9eZmZnYsmVLMJpARGEmkInl9LQemHPvEJ9mMXmii+4OoGV4W1w87piWv/gDMVwVaCEbYiKirsubKadFJTX2xLIsKygqqfHrF3t6Wg+vXi9bzKgpPgHZw/CMuHDp+gKd07HvDMnpmgoMLRggiCiovH0yyOgdD6NBD1lWYDD4PyzkDdliRkPB8tYponnzXf7iNqZmotkQaU8EB2JYSIsYIIgoqLx9Mgj0sJA3rMV7bb/0AUCRYC3e6zJAaHlYKJAYIIgoqDJ6x8Og10GSBfR6nccng77GSlwZXQijMROA+gHCedK8p0n0Wh0WCiQGCCIKOuH0b1cCOZXUW5H9RkIq2gMoMqA3ILLfSFU/T+sYIIioXb7WMfKkqKTGvsBVKMLtEJN9Kilap5KqHSAMyemIyZsfkDUEnQEDBBF5FOg6Rt4mnwM5ldQXhuR0xF83RHurwUOAAYKIPFJjuqk3yWflzEmPx6Q+BgiiLsyboSM1ppt6sybBl4QxqYMBgqiL8nboKBTTTQHAkNgHktMxBRcDBFEXVVRSA0lSIGArge1p6MjbVciB5Lw62d/VyuQ7BgiiLuqybhGtKWBhO/aXN2UqvOWclA5WkppaMUAQdRHO+Qa1SmAHYiczIPD1jsh3DBBEXYCrfENLCexAJZ8DvW6hq9Q70jIGCKIuwNVU1YnDr9JsCWyg69Q70jIGCKJO6tIhJXdTVQOZfFYjqdwV6h1pGQMEUSfkakhJ7amqTCp3PgwQRJ2QqymsE4dfpepUVa587nyCvic1EalPjSms7eHK586HAYKoEwr0FFZvOK905srn8McAQdQJmMtq8ff9P8BcVgsA9imseh1gNPo/hVW2mNF0eBtki9ntNVy30PkwB0EU5tzVVApUUrpln2YoEpq5T3OXwgBBFGacV0S7K8cdqCmsvu7TzM12Og8GCKIw4m5FdEfLccsWc7sL0Xzdp5mb7XQeDBBEYcTd9NWODCd5u+cz92nuuhggiMKIu+mrHRlO8rZ2Uss+zSx50fUwQBCFkUBOX/WldhJLXnRNnOZKFALO01K9Fcjpq1z5TO3hEwRRkHm71acrgZy+ypXP1B4GCKIg82WrT1cCNX2Vez5TezjERBRkoaiT5Ar3fKb2MEAQBVko6iS5wvLc1B4GCKIgsy1s0wEADAad33WSOopJamoPAwRRCAinf4eyDe6OiVQNEAUFBZgwYQKys7OxYcOGNuePHj2Ku+++G3fccQceeeQRnDt3Ts3mEGlCUUkNFNn261hRBIpKakLSDpbnpvaoFiAsFgtWrVqFjRs3YuvWrdi0aRPMZsdSwUuXLsWsWbPw0UcfoW/fvvjTn/6kVnOINEMrSWoOMVF7VAsQ+/btw7BhwxAXF4eYmBjk5ORg+/btDtcoioLz588DABobGxEdHa1Wc4g0I1hJ6vb2cOAQE7VHtXUQFRUVMJlM9uOkpCQcOXLE4Zr58+dj+vTpeOGFF9CtWzds3rzZp89ISNDWrAuTKTbUTfBbuPchHNqfYoqFQa+DIgQijHoMG5jm0O5A9OFCaRF+/PsKCFmC1WBEz/sXI7pXhuM1P70dp4v3ALIMGAxI+untiA7Qzy8c7kN72AcVA4SiKNDpdPZjIYTD8YULF7Bw4UK89dZbGDhwIN58803MmzcP69at8/ozqqrqoSja+LvHZIoN+xLH4d6HcGi/uawW67Z+A1kR0Ot1+OXPfoKEyyLs7Q5UH5q++18IyVaIT0hWnPnufxEVlep4UVQqYnJbi/DVRaWiLgCfHQ73oT2dtQ96vc6nP6xVG2JKSUlBZWWl/biyshJJSUn24+LiYkRFRWHgwIEAgHvuuQdfffWVWs0hUoWvNZVaNvcBAAjRoeElb7b/9LYQnyE5HVFDclmIj1xSLUCMGDEC+/fvR3V1NRobG7Fjxw6MGjXKfr5Pnz4oLy/H8ePHAQCffvopsrKy1GoOUcC11FT62xfH8ft3D3sVJFo299Hr4PPmPkDrHg7Nh95Hw7YV7vML3B+aAkC1Iabk5GTMnj0b06ZNg9VqxZQpUzBw4EDk5+dj1qxZyMrKwrJly/Dkk09CCIGEhAS88MILajWHKOA6UlPJ32J73u7hwP2hKRBULdaXl5eHvLw8h++tX7/e/vXo0aMxevRoNZtApJqOTlf1p9ieL0NHMblzuckP+YXVXIk6KBQ1lXwpsMdNfshfLLVB1EGhqKnEAnsUTAwQRH4Idk0lrn6mYGKAIHLBm+mrRSU19nU4Ikg1lbj6mYKJOQgiJ95uCdoyZVWWlQ5NWe0I7gJHwcQAQeSkZTGbEIDsYfpqIPeHBmxrHGqKT0Du0ddtcpnrGyiYGCCInGT0jodBr4Mk28pheHoyCNT+0LLFjIaC5Wi4uG4hJm8+1zdQyDFAELkQ7OSztXiv7Zc+ACgSrMV7XQYIrm+gYGKAIHLiakOfQDwleOJL8pnrGyhYOIuJyEkoNvTh7m6kRQwQRE60vkKaKFgYIIiccIU0kQ0DBJELXCFNxABBXQxXSBN5j7OYqMvgCmki3zBAUJfh7QY/gV4h7Q0mqUmLGCCoy/Bl+mqgVkh7i0lq0iLmIKjLCMX0VW8xSU1axABBXUYopq96i0lq0iIGCOpSgj191VtcSU1axABBXYarGktawSQ1aREDBHUZoaix5C0mqUmLGCCoy2CSmsg3DBDUZYQqSS1bzGg6vA2yxez2GiapSYu4DoK6lGAnqWWLGQ3bVgCKFc36CMTkznW5l0Nkv5GQivYAigzoDYjsNzJILSRyjwGCuoxQbAQknS4EZCsAAchWSKcL3e8UlzcfUbUn0ORhT2qiYGKAoC4jFElqW7K59bnFU/LZkJyO+OuGoLKyTvV2EXmDOQjqMkKRpGbymcIZAwR1GaFIUjP5TOGMAYK6lGAnqblCmsIZAwR1GaFYSc0V0hTOGCCoywhdktr9MZGWMUBQl+GclGaSmsgzVQNEQUEBJkyYgOzsbGzYsKHN+ePHj+OBBx7AHXfcgV/96leorXW/TzCRv5yfGILxBMEkNYUz1QKExWLBqlWrsHHjRmzduhWbNm2C2dxaakAIgUcffRT5+fn46KOP0L9/f6xbt06t5hCFZJork9QUzlQLEPv27cOwYcMQFxeHmJgY5OTkYPv27fbzR48eRUxMDEaNGgUAmDFjBu6//361mkNdkLmsFn/f/wPMZbYn04ze8Ygw6qHXAUajPjjTXJmkpjCm2krqiooKmEwm+3FSUhKOHDliPy4pKUFiYiIWLFiAY8eO4eqrr8Yzzzzj02ckJGgr4WcyxYa6CX4L9z60tL/wh2qs/OthSJICo1GPpTNGYvjgXnghLgbffH8GWdckIvOqK/z+vAulRWg8eRTd+gxAdK+MNufPJSbizCXHlycm4vJ2fsbhfg8A9kEr/O2DagFCURTodDr7sRDC4ViSJHz11Vd45513kJWVhZdffhnLly/H8uXLvf6Mqqp6KIo2RnVNptiwL5EQ7n24tP0HjpTBalUgAFglBQeOlCHhsggkXBaBWwf2BAC/+3ppIb4aN4X4ms5cGh50OHfmDJo8fG643wOAfdAKV33Q63U+/WGt2hBTSkoKKisr7ceVlZVISkqyH5tMJvTp0wdZWVkAgNzcXIcnDCJ/BGNKq70Qn2gtxOfMmJoJGCIBnR4wRNiOicKEagFixIgR2L9/P6qrq9HY2IgdO3bY8w0AMGTIEFRXV6Ow0PY/1a5duzBgwAC1mkNdTDAS0t4U4jMkpyMmdy4ib5jsttQ3kVapNsSUnJyM2bNnY9q0abBarZgyZQoGDhyI/Px8zJo1C1lZWXjllVewaNEiNDY2IiUlBStWrFCrOdTFXNYtAga9DooQMBrUSUh7m4A2JKczMFBYUrXcd15eHvLy8hy+t379evvXgwYNwpYtW9RsAnVB5rJavPvJ/0FWBPR6He4d+xNV9n3gKmnq7LgfBHUK5rJa7D7yI3olxKCopAaSrNhOCKHaegeukqbOzmMOYtKkSfavv/jiC9UbQ9QR5rJa/P7dw3jnn8fw+3cP47JuETAabOsdDCoNLwFcJU2dn8cnCCFa/5NftWqVQ5KZSCuKSmogSbYprUJWcL7Rijn3DkFRSQ0yesertq2oIbEPJKdjos7EY4BwXsdApEWuprSmp/VQfb9prpKmzs7raa6XBgsiLQlFjSWASWrq/Dw+QZw7dw47d+6EEAJ1dXXYsWOHw/ns7GxVG0fkjZatRGVZBG0rUYBJaur8PAaI1NRUvP322wCAnj174i9/+Yv9nE6nY4AgzQj2VqKuPouDsNTZeAwQlwYEIq1q2UpUoHUrUbXzDwCT1NT5tbsO4vz589i2bRuKi4sRHR2NjIwMjBs3DpGRkcFoH1G7QrGVKMAkNXV+HpPUJ0+exMSJE7Fjxw5ERUUBALZs2YJx48ahrKwsKA0kao8aSWrZYkbT4W2QLWa31zBJTZ2dxyeI1atXY/bs2bjzzjsdvv/ee+9h5cqVWLVqlaqNI/JGoJPUl5bxbnZTxhtwfmLQ8QmCOh2PTxDFxcVtggMA/PznP8eJEydUaxSRrwKZpPamjDfAUt7U+Xl8gjAYDG7PcV0EaUWgk9TelPEGWkt5S6cLYUzNZMVW6nS8XklNpFWBTlL7knxmKW/qzDwGiPLycjz//PMuz1ksFlUaROSrQCepmXwmsvEYIO6//3635+67776AN4bIFXNZrcfCexm94xFh1EOWlYBUb+UKaSIbjwHi8ccfb/O95uZmroGgoGkp5S3JCowGPebcO6RNkEhP64E59w5BaVUDeiXE+Jx/kC1mhzwCV0gT2XicxdTc3Ix58+Zh586d9u/9+te/xtNPPw1Jkjy8kigw7KW8BSDJCopKalxel57WAz//Wb8OBYeGbSvQfOh9NGxbAdliRmS/kYDeCEAH6I22Y6IuyGOAWL16Nerr63H99dfbv7dkyRLU1tZizZo1qjeOKNAJaOcFcNLpQkC5OKVVkSCdLrTNTsqbj8gb70ZM3nwmoanL8jjEtHv3bmzZsgXR0dH27yUnJ2PFihW45557MHv2bNUbSF2bc8LZnwS0bDGjoWA5oEho1hsRkzcfxtRMNOsjAEUC9Eb7WgbOTiJqJ0BEREQ4BIcW3bt3Zx6CgsL5icGfJwhr8V5bIAAARYK1eC+ib3mQaxmI3PAYIPR6Perr69G9u+M0v/r6euYgKCgCOYXVXfKZTwtErnnMQeTm5mLRokVoaGiwf6+hoQGLFi3iXhAUFC11lgD4XWfJuRw3y3MTeeYxQDz44IOIjY3FyJEj8Ytf/AJTpkzByJEjcfnll+Oxxx4LVhupiwtUnSWubyDyTbtDTM899xxmzJiBo0ePQq/XY+DAgUhKSgpW+6iLa6mzBPhfZ4nrG4h84zFAnD59GqmpqUhLS0NaWprDuS+++AKjRo1StXFEgZzmyh3giHzjcYjp0mGkX//61w7nuBcEBUNAk9TcAY7IJx6fIIRofQg/deqU23NEHdFejSXA+zpLssWMmuITkHv0dTsjiUX4iHzjdblv59LfLAVO/vCmxhLQWmfJUyBpKZfRcHGxG3eAIwoMr58giALJXmMJrTWW3D1FpKf18JiYtu8Ah9Yd4FwFCGNqJpoNkW1WTRORax4DhKIoqK2thRACsizbvwYAWZaD0kDqnAKZfOYOcETq8BggiouLMWzYMHtQuOmmm+znOMRE/ghV8pmrpom85zFAFBa63qydyF8tK6QlWfi9QprJZyJ1eJzm6q+CggJMmDAB2dnZ2LBhg9vrdu/ejTFjxqjZFNIgrpAm0jbVAoTFYsGqVauwceNGbN26FZs2bYLZbG5z3ZkzZ/Diiy+q1QzSKFcrpDuKK6SJ1KFagNi3bx+GDRuGuLg4xMTEICcnB9u3b29z3aJFi1xubUqdW6BXSHs6JqKO8ZiD8EdFRQVMJpP9OCkpCUeOHHG45u2338a1116LQYMGdegzEhK0NdZsMsWGugl+C1of9Po2xx397JpiK5ouOY4xWBEfxveC/x1pA/ugYoBQFMVhppMQwuG4uLgYO3bswFtvvYXy8vIOfUZVVT0URRsDCiZTLCor60LdDL8EtQ+K0ua4o5/dLDs+fTTIEZDC9F7wvyNt6Kx90Ot1Pv1hrdoQU0pKCiorK+3HlZWVDlVgt2/fjsrKStx99914+OGHUVFRgfvuu0+t5pDGBHKaK5PUROpQLUCMGDEC+/fvR3V1NRobG7Fjxw6H6q+zZs3Cxx9/jA8//BDr1q1DUlISNm7cqFZzKMTMZbX4+/4fYC6rBdBaY0mvA4xG9zWW3JEtZjQd3gbZYmaSmkglqg0xJScnY/bs2Zg2bRqsViumTJmCgQMHIj8/H7NmzUJWVpZaH00a467uUns1ltxpqb0ExYpmfQSiRtwHSW8EFBnQGxDZb6SKvSHqOlQLEACQl5eHvLw8h++tX7++zXW9evXCrl271GwKBZFzldaikhpIsgIhAPmSukvt1VhqIVvMDuUxpNOFgGK1TX9SJIgL9YjJm4+o2hNo8lDNlYh8o2qAoK7H1dOCbdV0+yW7XZEtZjQULAcUCc16I2Ly5tuK7ukjHIruGZLTEX/dkLBPLBJpCQMEBZSrp4WJw6/q8HCStXivLRAAgCLBWrwX0bc8yKJ7REHAAEEBldE7Hga9rcaSXt9aY8nb4SRn7hLQLLpHpD5VazFR1xSoGksAV0kThRIDBAVUIGssAdxHmiiUGCAooAJZYwlgKW+iUGKAIL84L4AL5AppgKukiUKJSWrqMHdTWiOMHZvS6gpXSROFDgMEdVigp7S6YkjsA8npmIiCgwGCvOa8QtrdAriOTml1hUlqotBhgCCvBLqekreYpCYKHQYI8kpRSQ0kSYEAIHWgnlJHMUlNFDqcxUReCfT0VW8xSU0UOgwQ5JVAT1/1FldSE4UOAwR5xZaQtm0ZazDo/J6+6i0mqYlChwGCvBbIGkveYpKaKHQYIMgrga6x5C0mqYlChwGCvBKsJPWle00DTFIThRKnuZJXgpGkdt5rOiZ3LiL7jYRUtIf7TROFAAMEeaUlSS3JImBJ6vb2mpZOFyJqSC5i8uZz9ziiEGCAIK8FMknt7V7TAHePIwoVBgjyiqsktT8rqLnXNJH2MUCQVwKdpOZe00Tax1lMBMBWjO+9T4vtG/84C3SSmiukibSPTxBkr9TaUra7pVLrpQK+ERBXSBNpHgMEOWz8g0sqtV4q0KW9uUKaSPsYIAgZveNh0OsgywJ6vfsprIEs7c0V0kTaxxwEAQj8FNZLV0N7+jx3x0QUegwQZJ/CKuB/naWW1dDNh95Hw7YVboNEZL+RgN4IQAfojVwhTaRBHGKigE5hlU4XArIVgABkK6TThS6nrRqS07lCmkjjGCAooFNYbcnm1gErT8lnrnkg0jYOMZG9zpIO/m8GxOQzUefBAEEAApekZvKZqPNQNUAUFBRgwoQJyM7OxoYNG9qc/+STT3DnnXfijjvuwMyZM1Fb63oVL6krkElqrpAm6jxUCxAWiwWrVq3Cxo0bsXXrVmzatAlmc+uMlvr6eixevBjr1q3DRx99hIyMDKxZs0at5pAHgUxSc4U0UeehWoDYt28fhg0bhri4OMTExCAnJwfbt2+3n7darXj22WeRnJwMAMjIyMCPP/6oVnPIA2+T1N6sb+AKaaLOQ7UAUVFRAZPJZD9OSkqCxWKxH8fHx+P2228HAFy4cAHr1q3D2LFj1WpOl2Uuq8Xf9//gtggf0FpnSa8DjEbXdZa8Xd/g+MSg4xMEURhTbZqroijQ6XT2YyGEw3GLuro6PPbYY8jMzMSkSZN8+oyEBG39dWoyxYa6CQ4Kf6jGyr8ehiQpMBr1WDpjJDKvuqLNdSZTLF6Ii8E3359B1jWJLq+pKT6Bhpb1DYoVUbUnEH/dkDbXXbj2evx4+CMIWYLOYETitdcjOog/F63dg45gH7SBfVAxQKSkpODQoUP248rKSiQlJTlcU1FRgV/96lcYNmwYFixY4PNnVFXVQ1G0MU/GZIpFZWVdqJvh4MCRMlitCgQAq6TgwJEyJFzmOr+QcFkEfv6zfqisrHPZj2Y5AvY5SUKgQY6A5Kq/UanoNrF105+6qFTUBennosV74Cv2QRs6ax/0ep1Pf1irNsQ0YsQI7N+/H9XV1WhsbMSOHTswatQo+3lZljFjxgyMHz8eCxcudPl0Qf7xJfksW8yo2fu3gAwdGZLTETUkl4vgiMKcak8QycnJmD17NqZNmwar1YopU6Zg4MCByM/Px6xZs1BeXo7vvvsOsizj448/BgBcd911WLp0qVpN6nKck82eks8NBcvRcHEv6Ji8+W1+uRtTM9FsiGyzXzQRdV6qltrIy8tDXl6ew/fWr18PAMjKykJhYaGaH9/lOT8xuHuCcLU/tHOAMCSnc79oCjuNjedRX38Wsiz59LqKCj0URVGpVeozGIyIjOwJfweJWIupE/N2+qq3q59ZO4nCSWPjedTV1SAuzoSIiEifhrGNRj0kKTwDhBACVmszfvyxHDExPdCt22Udfi+W2ujEWmosAZ5rLHH1M3VG9fVnERdnQmRkVJfKcep0OkRGRiEuLhH19Wf9ei8GiE7OmxpLXP1MnZEsS4iIiAx1M0ImMjLK56E1ZwwQnYjzoriikhr01lVgbPQ36K2rcFtjiaufqbPqSk8OzgLRd+YgOglzWS1+/+5hSLICo0GPOfcOwXWxZzEsdgcMkCHDgPrYDJev5epnouA5f74eM2ZMx4oVL6Nnz1T84x8F2Ljxbej1elx//Y14/PEnUVd3DrNnP+7wmrNna7Bz5x4cPvw1Fi6ci6QkW5mifv0ysGDBs6q0lQGikygqqYEkKxACkGUFRSU1GBtdiiadAh0APRSkWEsBDG7zWk5hJQqOo0e/xYoVz+PUqRIAQEnJD1i//lWsX/82EhMTsXLlcmzZ8lf88pdT8dZbGwHYqlI88cSjyM+fCQAoLDyGe++digce+G/V28shpjDlPJxkS0jb6ikZDLZ6Srro7tBdzD7oPOzu1jKFNX70vYjJncuZStTl/V/p2XZrmHVEQcEHeOqpeUhMtNWpM5vNGDAgC4mJiQCAkSNvxp49nzu85h//+AjR0dHIzh4HACgsPIqvvjqABx/8JebNmw2LpTygbbwUnyDCkKvhpPS0Hphz7xAUldQgo3c80tN6oKnC++SzITkd8dcNCfvyAkT+MpfVYuW7h2F1+v8rEObPf8bhOD39J1i7dhUslnIkJprw2Wefoqqqyn5elmX8+c9vYvnyP9i/1717LMaMuR2jR4/B1q1bsHjxArz66hsBaZ8zBogwVFRSA0my1ViSLg4npaf1sP/TgslnIt+5Gq4NVIBw1rt3H8yY8Tjmz38KUVHRGDNmLI4dO2o//+WX+3HllVfimmtan+rnzGmtW3fXXVPw2mtrUV9fj+7dA///N4eYwpC3NZa4PzSR71wN16qlqakJ/fsPwJtvbsRrr72BxMQkpKX1sp/fs2c3fvazbPuxoij485//BFmWHd7HYDCo0j4GiDAU6BXSRNQqPa0H5k29HpNGXR3Q4SVXLlxoxBNPPIqGhvOwWq14//1NGDOmNSB8++0RDBrUWlZfr9fjiy92Y/fuXQCAf/5zG6699jp069ZNlfZxiCkMtayQlmTR7gppyemYiNr3k15x6Jtyueqf06NHHKZPz8fDD/83JEnC7bfn2JPRAHD6dBlMJsdtEhYuXIwVK5bizTfXIz4+HosW/U619jFAhCmukCYKX1u2FNi/zs29C7m5d7m87tNP97b53tVXX4PXXlMnKe2MQ0xhiCukibQCC3EAABKuSURBVCgY+AQRhpKkHzHzkhXSx6XeAK5qcx1XSBORP/gEoTHOC+Bciao2wwgZBh1ghIyoate7wBlTMwFDJKDTA4YIrpAmIp/wCUJDzGW12Lz5Y/TV/4jNB3riF7/IcTmDIiEpAbpy2xRX3cVjV7jJDxH5gwFCQ04fO4JHYj6+OHR0BEeP9UR62i1trkuMVtCkswUHAR0So91vbMJNfoioozjEpCHpERaHoaP0CIvL64ypmdBdHDrSceiIiFTCJwgNEZExtqeCi0NHIjLG5XUcOiKiYGCACCLZYvb4S72qogqXATDoAFnYjtPcvBeHjojCzxtvrMOuXZ8AAEaMGImZM5/AwYNfYu3aVWhqasKYMbfj4YdnOrzmued+i6FDb8SECXkAgCNH/o01a16C1SqhR48eePrp3yIlpacq7eUQU5DIFjMaCpaj+eAWNBQsh2xpO/Oo6Yp0yDBAFjrIMKDpCgYAos7i4MEvcfDgAbz55ga89dZGFBUVYufO7Vi2bAmWLfsD3nnnPRQWfof9+22L486cqcTcubOxe/enDu+zZMkzmDfvGbz11kbcfvs4vPzy71VrMwNEkFiL99o25AEARbIdO6kw9sQrddn4R+Ng/E9dNiqM6vxVQESeSeX/h6bD21z+IddRCQmJeOyx2YiIiIDRaESfPlfh1KkSXHllb6SmpsFoNCI7ezw++8z2hLFjxz9xyy2jMWbM7fb3aG5uRn7+o0hP/wkAW7lw7gehcbLFjJriE5B79HU77CMu/qO75GtnGb3j8dHeZJQ0mWAw6HGvilUkicg12WJGw99XALIVzfqIgG2idfXV19i/PnWqBLt2fYIpU+5BQkKi/fsJCYmorKwAANx33zQAtiGlFpGRkcjJmQDAVtn1jTfW4ZZbbvW7be4wQPhJtpjRsG0FGi5u1+nuP6aK+EGIFZ/DAAUy9KiJHwTn0nmuNv0houCSThcCsmSbLaJIkE4XBjTfd/z495g790k89tgTMBgMOHXq0jL8Ajpd+wM7VqsVzz//LCRJxrRp0wPWNmcMEH6SThdCyFboICBkq9v/mL6ti8P/1uXgGmM5vpdScH1dXJsAAaDNpj9EFFy2PdqNtiAR4D3ajxz5NxYtmodZs57C2LE5OHz4a5w507qDXFVVlX37UXcaGhowf/5TuPzyHli+/A8wGtX7Nc4A4aczF/SIFeLikJHAmQt6lzOPWoaPTnL4iEjTDMnpiL1jHppOHQvoNHKLpRwLFvwGv/vdMgwdeiMA4Nprr8OpUydRWnoKPXumYufOjzFx4h0e3+e5555BWtqVmDPnaej16qaRGSD8VF96HLEAdDrbE2l96XGX13H4iCh8GFN+AiRe0/6FPnj33XfQ1NSMNWtW2b93112TsWDBs1i4cC6am5swfPhI3Hbbz9y+R3FxIfbs+RxXXXU1pk+fCgBITEzEypWrA9rWFjohRNhuNFZVVQ9FUa/57a1bAIDjH/wPEiu+sgeIM0k/xdWTZrq8VutMplhUVtaFuhkdFu7tB9iHQCovP4mUlI5tkmU06iFJ7kvYhAOjUY/S0hMOPwO9XoeEBO/L/vMJwo2W5DMUzzMZSmOzEFdxCAZhSz6Xxmbh6hC0l4go0LpkgPDmyUA6XQgo1nZnMqT2H4hXj+Tgan05jispuKf/QLWbT0QUFJ0uQLT3y9/bJwNjaiaa9RG2xW3tzGT4QU7C8WYT9AZdQPtCRBRKnSpAePPL39snA0NyOqpvmona779Bj2uyEOvmSaOopAaKYpvFJBSBopIaJqCJNEEHIRSv1hV0RkIosC3N7ThVf3IFBQWYMGECsrOzsWHDhjbnjx07hsmTJyMnJwcLFy6EJEl+fZ5tgcvFX/4X1yQ4M6ZmAvoI2y5rHp4MzGW1eGH7Wfy/wjS8sP2s2x3eMnrHw2jQQ68DDAY9Mjh9lUgTIiOjcfbsGUiSFWE8F8dnQghIkhXV1WcQGRnt13up9gRhsViwatUq/O1vf0NkZCR++ctf4qabbkJ6eutf4nPmzMHzzz+PwYMHY8GCBdi8eTPuu+++Dn+mLro7WotYiIvHjrwtlV1UUgNJVmyxRlbcPhm0TF8trWpAr4QYPj0QaUR8vAn19bWorrZAUWSfXqvX66Eo4TuLSa83wGRKAKDRALFv3z4MGzYMcXFxAICcnBxs374djz/+OACgrKwMFy5cwODBgwEAkydPxurVq/0KEOJC/SVHOqfjVt6Uym55MpBlpd0ng/S0Hhg+uJcmpvYRkY1Op0NsbBxiY+N8fq1Wpur6IxB9UC1AVFRUwGQy2Y+TkpJw5MgRt+dNJhMsFtc7qHnLtkQ+0qvEcnu4sI2IujrVAoSiKNDpWhMkQgiH4/bOe6PNgg/TEFyIW4zGk0fRrc8ARPfK6FjjW97OFIvhg3v5dH24C/c+hHv7AfZBK9gHFQNESkoKDh06ZD+urKxEUlKSw/nKykr78ZkzZxzOe6Om5nzbldTd0oDMNJwHcL7K9RCTGhISuqMqiJ+nhnDvQ7i3H2AftKKz9kGv1yE+/jKv30O1ADFixAisWbMG1dXV6NatG3bs2IHnnnvOfj4tLQ1RUVH4+uuvMXToUHz44YcYNWqUT5/hS0eDwZcl7FoV7n0I9/YD7INWsA8q12IqKCjA66+/DqvViilTpiA/Px/5+fmYNWsWsrKyUFhYiEWLFqG+vh4DBgzAsmXLEBkZqVZziIjIB2FdrI+IiNTTNZcYEhFRuxggiIjIJQYIIiJyiQGCiIhcYoAgIiKXGCCIiMglBggiInKJAcJHH374ISZOnIiJEyfixRdfbHP+9OnTuP/++zFu3Dg8+uijOH/+fAha6V577f/ggw9w8803484778Sdd96JVatWhaCVnq1btw45OTnIy8vDq6++2ua81u8B0H4ftHof6uvrkZubi9LSUgC2qs15eXnIzs5220at3Y+O9EFr98O5DwBgtVrx4IMP4ssvv3T5mg7dB0Fea2hoEDfeeKOoqqoSVqtVTJkyRezdu9fhmocfflhs27ZNCCHE2rVrxYoVK0LRVJe8af+SJUtEQUFBiFrYvr1794rc3FxRV1cnJEkSjzzyiPj4448drtHyPRDCuz5o8T78+9//Frm5uWLAgAHi1KlTorGxUYwePVqUlJQIq9Uqpk+fLnbv3t3mdVq6Hx3tg5buh3MfhBDi+++/F/fcc4/IysoSBw4ccPm6jtwHPkH4QJZlKIqCxsZGSJIESZIQFRVlP2+1WnHw4EHk5OQAsO1xsX379lA1t4322g8A33zzDT744APk5eXhN7/5DWprXe+kFyrfffcdbr75ZnTv3h0GgwG33HILPvnkE/t5rd8DoP0+ANq8D5s3b8azzz5rL6p55MgR9OnTB1deeSWMRiPy8vLa/Ky1dj860gdAW/fDuQ8AsGXLFjz00EMYNGiQy9d09D4wQPige/fueOKJJzB+/HiMHj0aaWlpuP766+3na2pq0L17dxiNthqIgdjjIpDaaz9ga/PMmTPx0UcfoWfPnliyZEmIWuvagAED8K9//Qtnz55FU1MTdu3ahTNnztjPa/0eAO33AdDmfVi6dCluuOEG+7GrPV+cf9Zaux8d6QOgrfvh3AcAmDt3LsaOHev2NR29DwwQPigsLMT777+Pzz77DHv27IFer8ef/vQn+3nhYk8LX/e4UFN77QeAV155BUOHDoVOp8NDDz2EPXv2hKi1rg0fPhyTJ0/GAw88gIceeghDhw5FRESE/bzW7wHQfh8A7d8HwLs9XbR+P7zdlyYc7ocnHb0PDBA++Ne//oXhw4cjISEBkZGRmDx5Mr766iv7+SuuuAJ1dXWQZdv+t857YIRae+2vq6vDW2+9ZT8WQsBgMISgpe7V19cjOzsbBQUF+Mtf/oLIyEhceeWV9vNavwdA+30Ih/sAtN3TxdXPWuv3w5s+hMv98KSj94EBwgeZmZnYt28fGhoaIITArl27kJWVZT8fERGBG264Af/4xz8AAFu3bvV5jws1tdf+mJgY/PGPf8R//vMfAMA777yD22+/PVTNdam0tBQzZ86EJEmoq6vDli1bMH78ePt5rd8DoP0+hMN9AIBBgwbhxIkTOHnyJGRZxrZt29r8rLV+P7zpQ7jcD086fB/8TKh3Oa+//rrIyckRubm54umnnxYXLlwQCxYsEJ988okQQojS0lIxdepUMX78eDF9+nRx9uzZELfYUXvtP3jwoLjrrrvEuHHjxIwZM8S5c+dC3OK21q5dK8aPHy+ys7PFxo0bhRAirO6BEO33Qcv34bbbbrPPntm3b5/Iy8sT2dnZYunSpUJRFCGE9u+Hr33Q4v24tA8tpk6d6jCLyd/7wP0giIjIJQ4xERGRSwwQRETkEgMEERG5xABBREQuMUAQEZFLDBCkaV9++SVyc3ND/n5DhgxxqJzpq/feew8bNmwAAKxZs8brUg2yLOORRx5pU4ojEL799ls888wzAX9f6jwYIIiC4Ouvv8aFCxd8ft0bb7yBn/70p0hMTAx4m6677jpIkoTPPvss4O9NnYMx1A0gak9DQwNmzZqFkydP4vLLL8eSJUvQt29fNDc3Y+XKlTh48CBkWca1116LRYsWoXv37vjss8/w+uuvo7m5GdXV1bjrrrvw5JNPOrzvoUOH8Jvf/AYvvfRSm6KFhw4dwnPPPQedToesrCwoimI/t2vXLrz66quwWq2Ijo7GvHnzMGTIEKxZswYnT55EeXk5KisrkZmZiaVLl2L//v3YtWsX9u7di+joaADA8ePH8cADD6CyshKJiYl46aWX2pQ+aGxsxJ///GcUFBQAsD15lJSUwGKxoLKyEgMGDMBNN92ErVu3orS0FHPmzEFubq7X1wHAPffcg8WLF+O2224L+H2jTiCQK/uIAu3AgQMiMzNTfP3110IIIf7617+KKVOmCCGEWLNmjVi+fLl95esf/vAH8eyzzwpFUcTUqVPFiRMnhBBClJeXi/79+4uqqipx4MABMXHiRLF//34xduxYcezYsTaf2dTUJEaMGCH27dsnhBCioKBA9OvXT5w6dUqcOHFC5ObmiurqaiGEEMXFxWLkyJHi/PnzYvXq1WLUqFGisrJSyLIsnnrqKbF8+XIhhBDz5s0Tf/zjH4UQQqxevVqMGTNGVFVVCSGEePTRR8XatWvbtGPXrl1i6tSp9uPVq1eL2267TZw7d040NjaKG2+8USxbtkwIIcTOnTtFdna2T9e1GDJkiCgpKfHpvlDXwCcI0ryMjAz7X/iTJk3C4sWLUVdXh927d6Ourg779u0DYKt5n5CQAJ1Oh9deew27d+/Gtm3b8P3330MIgcbGRgBAeXk5ZsyYgXvvvReZmZltPq+4uBhGoxHDhw8HAOTm5uK3v/0tAGDv3r2oqKjAf/3Xf9mv1+l0KCkpAQCMGzfOPhw0ZcoUvPDCC5g3b16bzxg5ciSuuOIKALYaWdXV1W2uOX78OHr37u3wvREjRiA2NhaArTT1LbfcAgDo3bs3zp496/N1ANCrVy+cOHHCoWAgEcAhJgoDer1jqkyn08FoNEJRFCxYsACjR48GAJw/fx5NTU1oaGjApEmTMHbsWNxwww24++678cknn0BcrCpjMBiwbt06zJw5E+PGjXO5yYpwqkDTUkdfURQMHz4cL7/8sv3cjz/+iKSkJOzcudOhyqeiKG3a7vx+Lf1x/ryW7186tAUAkZGRbt+nI9e1nAu36qQUHExSk+YVFRXh2LFjAIBNmzZh6NCh6NatG26++WZs2LABzc3NUBQFzzzzDF566SWcPHkS9fX1ePLJJzFmzBh8+eWX9msA22Yp119/PebNm4e5c+fanyxaZGRkQAiBzz//HADw6aef2ncQGz58OPbu3Yvvv/8eAPD555/jjjvusCegP/30U9TV1UFRFGzevNk+tm8wGCBJkk/97tu3L06dOtXBn5p3hBA4ffo0+vbtq+rnUHjiEwRp3tVXX421a9fi1KlTSEhIwPLlywEAM2fOxIsvvohJkyZBlmX0798f8+fPR0xMDG699VaMHz8ekZGR6NevH9LT03Hy5EmHv6wnTZqEjz/+GMuXL8fvfvc7+/cjIiLwyiuvYPHixXjppZfQv39/JCQkAADS09OxZMkSPPXUUxBCwGg04tVXX8Vll10GAEhMTER+fj5qampw4403YsaMGQCAUaNG2dvtrREjRmDhwoU4d+4cLr/8cr9+hu5888036N27N1JTU1V5fwpvrOZKFCBr1qxBTU2NPV8RCK+99hoMBgPy8/MD9p6Xmj9/PsaNG4dbb71Vlfen8MYhJiINmz59Og4cOOCw61mgfPvtt9DpdAwO5BafIIiIyCU+QRARkUsMEERE5BIDBBERucQAQURELjFAEBGRSwwQRETk0v8H0lAR1rRvzugAAAAASUVORK5CYII=\n",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Compute ECDFs\n",
"x_1975, y_1975 = ecdf(bd_1975)\n",
"x_2012, y_2012 = ecdf(bd_2012)\n",
"\n",
"# Plot the ECDFs\n",
"_ = plt.plot(x_1975, y_1975, marker='.', linestyle='none')\n",
"_ = plt.plot(x_2012, y_2012, marker='.', linestyle='none')\n",
"\n",
"# Set margins\n",
"plt.margins(0.02)\n",
"\n",
"# Add axis labels and legend\n",
"_ = plt.xlabel('beak depth (mm)')\n",
"_ = plt.ylabel('ECDF')\n",
"_ = plt.legend(('1975', '2012'), loc='lower right')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Parameter estimates of beak depths\n",
"Estimate the difference of the mean beak depth of the G. scandens samples from 1975 and 2012 and report a 95% confidence interval."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"def bootstrap_replicate_1d(data, func):\n",
" \"\"\"Generate bootstrap replicate of 1D data.\"\"\"\n",
" bs_sample = np.random.choice(data, len(data))\n",
" return func(bs_sample)\n",
"\n",
"def draw_bs_reps(data, func, size=1):\n",
" \"\"\"Draw bootstrap replicates.\"\"\"\n",
" \n",
" # Initialize array of replicates: bs_replicates\n",
" bs_replicates = np.empty(size)\n",
" \n",
" # Generate replicates\n",
" for i in range(size):\n",
" bs_replicates[i] = bootstrap_replicate_1d(data, func)\n",
" \n",
" return bs_replicates"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"difference of means = 0.22622047244094645 mm\n",
"95% confidence interval = [0.05979819 0.39109804] mm\n"
]
}
],
"source": [
"# Compute the difference of the sample means: mean_diff\n",
"mean_diff = np.mean(bd_2012) - np.mean(bd_1975)\n",
"\n",
"# Get bootstrap replicates of means\n",
"bs_replicates_1975 = draw_bs_reps(bd_1975, np.mean, 10000)\n",
"bs_replicates_2012 = draw_bs_reps(bd_2012, np.mean, 10000)\n",
"\n",
"# Compute samples of difference of means: bs_diff_replicates\n",
"bs_diff_replicates = bs_replicates_2012 - bs_replicates_1975\n",
"\n",
"# Compute 95% confidence interval: conf_int\n",
"conf_int = np.percentile(bs_diff_replicates, [2.5, 97.5])\n",
"\n",
"# Print the results\n",
"print('difference of means =', mean_diff, 'mm')\n",
"print('95% confidence interval =', conf_int, 'mm')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Hypothesis test: Are beaks deeper in 2012?\n",
"Your plot of the ECDF and determination of the confidence interval make it pretty clear that the beaks of G. scandens on Daphne Major have gotten deeper. But is it possible that this effect is just due to random chance? In other words, what is the probability that we would get the observed difference in mean beak depth if the means were the same?\n",
"\n",
"Be careful! The hypothesis we are testing is not that the beak depths come from the same distribution. For that we could use a permutation test. The hypothesis is that the means are equal. To perform this hypothesis test, we need to shift the two data sets so that they have the same mean and then use bootstrap sampling to compute the difference of means."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"p = 0.0034\n"
]
}
],
"source": [
"# Compute mean of combined data set: combined_mean\n",
"combined_mean = np.mean(np.concatenate((bd_1975, bd_2012)))\n",
"\n",
"# Shift the samples\n",
"bd_1975_shifted = bd_1975 - np.mean(bd_1975) + combined_mean\n",
"bd_2012_shifted = bd_2012 - np.mean(bd_2012) + combined_mean\n",
"\n",
"# Get bootstrap replicates of shifted data sets\n",
"bs_replicates_1975 = draw_bs_reps(bd_1975_shifted, np.mean, 10000)\n",
"bs_replicates_2012 = draw_bs_reps(bd_2012_shifted, np.mean, 10000)\n",
"\n",
"# Compute replicates of difference of means: bs_diff_replicates\n",
"bs_diff_replicates = bs_replicates_2012 - bs_replicates_1975\n",
"\n",
"# Compute the p-value\n",
"p = np.sum(bs_diff_replicates >= mean_diff) / len(bs_diff_replicates)\n",
"\n",
"# Print p-value\n",
"print('p =', p)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We get a p-value of 0.0033, which suggests that there is a statistically significant difference. But remember: it is very important to know how different they are! In the previous exercise, you got a difference of 0.2 mm between the means. You should combine this with the statistical significance. Changing by 0.2 mm in 37 years is substantial by evolutionary standards. If it kept changing at that rate, the beak depth would double in only 400 years."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Variation in beak shapes\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### EDA of beak length and depth\n",
"Make scatter plots of beak depth (y-axis) versus beak length (x-axis) for the 1975 and 2012 specimens."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"bl_1975 = np.array(df[(df['year'] == 1975) & (df['species'] == 'scandens')]['beak_length'])\n",
"bl_2012 = np.array(df[(df['year'] == 2012) & (df['species'] == 'scandens')]['beak_length'])"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEJCAYAAACKWmBmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXwT5do//k+SSaGlQaRNBdJCwcruw1IQUDbLVmQTOMcfIggqBxRUwKdA2SygyCLgAohwQOQnIKI9wEERj6yylVMqi4CVB6EQUqALW0u3TDLfP2KSbpkkk5nJdr1fL1+Sbeaau0nuzH1fc90KjuM4EEIIIQ4ovR0AIYQQ30YdBSGEEF7UURBCCOFFHQUhhBBe1FEQQgjhRR0FIYQQXtRREEII4cV4OwAp3L37EGZzYF4eEhERjvz8Qm+H4ROoLeyoLeyoLexcbQulUoFHH63l8PGA7CjMZi5gOwoAAX1s7qK2sKO2sKO2sBOjLWjoiRBCCC/qKAghhPAKyKGn6nAch7t3c1FWVgLAf09Lc3KUMJvNbrxCgZCQmnj0US0UCoVkcRFCApekHUVhYSFGjBiBzz//HNHR0QAAo9GIcePGYeLEiejUqVOV12RnZ2PatGnIz89H48aNsWzZMtSq5XiSxfVY7kOhUOCxx6KhUPjviRTDKMGyrncUHGfGvXt5KCy8D42mjoSREUIClWTfmGfPnsWLL76IrKws231XrlzB6NGjcfr0aYevmz9/PkaOHIm9e/eidevW+Oyzz0SJp7i4EBpNHb/uJIRQKJTQaB5FcTFlgRBChJHsW3P79u1ISUlBVFSU7b7vvvsO48aNQ5s2bap9jdFoRHp6Ovr16wcAGDZsGPbu3StKPGazCSpV0Iy0VaBSMTCbTd4OgzihNNyAOu04lIYb3g6FkAok++ZcuHBhlfumT58OANi0aVO1r7l79y7Cw8PBMJawtFotbt++LVpMwTpGH6zH7U+UhhsIW7YYCtYIjlGjKCkZZl20t8MiBICPTWZzHFflS03Il1xERHiV+3JylGAY3xp2eviwEP/4xytYtuwTNGjQAN9//29s3rwJKpUK8fEd8fbbU/HgwQNMnjzR9prCwkLcu3cXBw8ew6+/ZiA5OQmPPfYYAKBp02aYO3d+tftSKpXQajWyHJecAuaYMvMARgHENQWyshBakAdoW7i1iYBpCxFQW9iJ0RY+1VHUrVsXBQUFMJlMUKlUyM3NrTB05ar8/MIqF5mYzWa3JoGlduHCeSxd+j6uX78Gk8mMK1euYO3az/DPf/7/iIyMxLJli7Ft21aMGDEKGzduBWA5hilTJuIf/5gIljXjwoULePHFURg9+hXbdh0do9lsRm5ugSzHJhetVhMwx6TURCKM5aDIvGQ5o9BEwuzGsQVSW3iK2sLO1bZQKhXV/sC2PS5mUJ5Sq9Xo0KED9uzZAwDYuXMnunfv7tWYDAYF0tJUMBjEHb7ZvXsH3nlnBiIjtQCAy5cvo1WrJxEZGQkAeOaZrjhy5HCF1+zZ82/UrFkTffsmAgAyMy/gv/9Nw5gxIzBjxlTcvn1L1BiJfMy6aBQlJaPkpTE07ER8jk90FLNnz8b+/fsBACkpKdi+fTuee+45nDp1ClOmTPFaXAaDAsuWhWDLFgbLloWI2lkkJ89FmzbtbLfj4p7AxYvncfv2LZhMJhw8uB/5+fm2x00mEzZt2oiJE9+y3RcersHf/vb/YdOmbejS5RnMmzdLtPiI/My6aBg7P02dBPE5kg89HThwoMp9X331VYXb5Se+dTpdlce9Ra9XgmWBmBgOer0Cer0SOp002UMNGzbC66+/ieTkd1CjRk0kJPTG779fsD1+8uQJxMTEIC7uCdvw0rRp9o7h+ef/hs8/X4XCwkKEhzs+hSSEEHf5xBmFr4qJMYNhAL1eAYax3JZKaWkpWrRohY0bt+Lzz79AZGQUdOV+WR45cgi9evW13Tabzdi0aQNMpoodl0qlkixGQkhwoo6Ch07HISmpDC+9xCIpqQw6nXSlP0pKijF58hsoKnoIo9GI1NRvkJBg7xjOnz9XYahKqVTil18O4dAhyxnbjz9+j5YtWyM0NFSyGAkhwcmnsp58kU7HSTbcVN4jj9TBq6/+A+PHvwKWZdGnTz/bpDUAZGcboNVWzACbPXseli5diI0b/4lHH30Uc+ZUnxpLCCGeUHAc578V8hyoLj321q1rqFevkZciEo+7tZ6sAuX4y6M0SDtqCztqC7uATI8lhBDie2joiRASMJSGG0BmHpSaSEozFhF1FISQgGCtlwVGgTCWowsXRURDT4SQgKDSX4eCNQKxsVCwRqj0170dUsCgjoIQEhBMMQ3BMWogKwsco4YppqG3QwoYNPRECAkI1npZoQV5lqKKNOwkGuooCCEBw6yLBrQt3Kq8S5yjjsJLvvhiHQ4c2AcAePrpZzBx4mSkp5/EqlUfobS0FAkJfTB+/MQKr3nvvXfRseNTSEwcCAA4d+4MVq5cAaORxSOPPIKZM99FvXr1ZT8WQkhgozkKJ6RYnjI9/STS09OwceMWfPnlVvzxRyZ+/nkvFi1agEWLlmPz5m+RmXkRJ04cAwDk5eVi+vSpOHRof4XtLFgwFzNmzMWXX25Fnz6J+PjjD0WLkRBCrKij4GFNt6u5ZRPCli0WrbOIiIjEpElToVarwTAMGjWKhV5/HTExDdGggQ4Mw6Bv3/44eNByxvGf//yIbt16ICGhj20bZWVl+Mc/3kBc3BMALGXKaT0KQogUqKPgYU23M8U0EjXdrkmTx9G69ZMAAL3+Og4c2AelUomIiEjbcyIiIpGbmwMAGDnyZQwa9HyFbYSEhKBfv+cAWCrJfvHFOnTr1lOU+AghpDzqKHhY0+1U+muSpNtdufInpk6dhEmTJqNBAx0qLg/OQaFw/ucxGo2YP38OWNaEl19+VdT4CCEEoMlsXtZ0O5X+OkwxDUVNtzt37gzmzJmBt99+B71798Pp0xnIy7OvaJefn29bFtWRoqIiJCe/g9q1H8HixcvBMPTnJISIj75ZnDDrokXPx759+xZmzUrC/PmLEB/fEQDQsmVr6PXXcOOGHvXrN8DPP/+EAQMG827nvffmQqeLwbRpM6FU0smhv1Aaboj644PqGxGpUUfhBV9/vRmlpWVYufIj233PPz8Ms2alYPbs6SgrK0WXLs/g2Wd7OdzGpUuZOHLkMGJjm+DVV0cBACIjI7Fs2aeSx0+EsyZIKFgjOEbtcT0iqm9E5EAdhRdMmZKEKVOSqn1s06avHb5u9ux5tvUomjZtjqNHT0kVIpFI+QQJlf4aVPrrHn2x2+obxTWFIvOSx9sjpDo0XkGIjMROkKD6RkQOkp9RFBYWYsSIEfj8888RHR2N48ePY9GiRSgtLUX//v0xderUKq/ZsWMHli9fjoiICABAz549q30eIf5G7AQJqm9E5CBpR3H27FnMmTMHWVlZAICSkhLMmjULX331FerXr48JEybg8OHD6NGjR4XXnT9/HsnJyRg4cKCo8XAcB0XFHNSgEICr3fo1sRMkqL4RkZqkQ0/bt29HSkoKoqKiAADnzp1Do0aNEBMTA4ZhMGjQIOzdu7fK63777Tfs2LEDgwYNQlJSEu7fv+9xLEqlCiYT6/F2/JHJxEKpVHk7jKAipPSLFOViCBGDpB3FwoUL0aFDB9vtnJwcaLVa2+2oqCjcvn27yuu0Wi0mTpyIf//736hfvz4WLFjgcSyhoeEoKLgHjjN7vC1/wnFmFBTcRWio44XTibiElH6RqlwMIWKQNevJbDZXGPpxNBS0evVq27/HjRuHPn36VHkOn4iIql+KERG1oNfrkZdnQDCNxCgUQK1atRATowvIay20Wo23Q6gqMw9gFEBcUyArC6EFeYC2hfivqcQn28JLqC3sxGgLWTuKevXqITc313Y7NzfXNixlVVBQgNTUVIwdOxaApTNRqdwbNsnPL4TZXLU3CAuri7Cwuu4H7kO0Wg1yBYxF5+c/lCAa7xLaFlJTaiIRxnJQZF6yXCuhiXQ6fyDkNeX5alt4A7WFnattoVQqqv2BbSVrR9GmTRtcvXoV165dQ3R0NL7//nsMHz68wnPCwsKwfv16tGvXDm3atMHmzZvdPqMgxJuEZDZJWS6GEE/J2lHUqFEDixcvxltvvYXS0lL06NEDiYmJAIDZs2cjISEBvXr1wscff4x58+ahpKQEsbGxWLp0qZxhEuIxIZlNUpSLIUQMCi4AcycdDT0FAjqttqO2sKO2sHPUFgaDAnq9EjExZuh0vvv9IGacfjn0RAgh3mAwKLBsWQhYFmAYICmpzCc7C1+NM/DSYAghpBK9XgmWBWJiOLCs5bYv8tU4fSMKQgiRUEyMGQwD6PUKMIzlti/y1Thp6IkQEvB0Og5JSWU+P0fhq3FSR0EICQo6HQedzuTtMJzyxThp6IkQQggv6igIIYTwoo6CEEIIL+ooCCGE8KKOghBCCC/qKAghhPCijoJ4lZyruom9LzljZzLSUXPDOjAZ6aJsz2BQIC1NBYMh+JYGJu6j6yiI11hXdVOwRssaDEnJklVPFXtfcsbOZKRD8+YEWAsAFaxaCza+o+Dt+Wo9IeK76IyCeI1Kfx0K1ghTTCMoWCNU+ut+sy85Y2fOnAZYFuYGOoBlLbc94Kv1hIjvoncI8RpTTENwjBoq/TVwjBqmmIZ+sy85Y2fbtgMYBspsA8Awltse8NV6QsR30XoUfibQ1h1QGm4IXtXN3bbwZF9ybI8Pk5EO5sxpsG3bVTvs5G5b+MvaDEIE2mfEE2KtR0EdhZ+hD4EdtYUdtYUdtYWdWB0FDT0R4gRlCBGx+Ot7ibKeCOFBGUJELP78XnLpjMJgMODo0aM4ceIEbt26JXVMhPgMyhAiYvHn9xLvGcWhQ4fwySef4Pr164iOjgbDMMjOzkajRo3wxhtvoEePHnLFSYhXUIYQEYs/v5ccTma/++67MJlMeOGFF9CmTZsKj507dw5ff/01FAoFPvjgA4cbLywsxIgRI/D5558jOjoax48fx6JFi1BaWor+/ftj6tSpVV6TnZ2NadOmIT8/H40bN8ayZctQq1Yttw6KJrODg1xt4Q8ZQvS+sPPltpD7vSR51tOFCxfQqlUr3o3zPefs2bOYM2cOrl69ir179yIyMhKJiYn46quvUL9+fUyYMAEvv/xylbOSCRMmYPDgwRgwYABWr16NoqIiTJs2zdlxVkAdRXCgtrCjtrCjtrCTPOvJWSfh7Dnbt29HSkoKoqKiAFjOQho1aoSYmBgwDINBgwZh7969FV5jNBqRnp6Ofv36AQCGDRtW5TmEAJZrGHDkiCx1lvyV0PpQctawEsLX4wtETrOe9uzZg08++QQPHjwAAHAcB4VCgRMnTvC+buHChRVu5+TkQKvV2m5HRUXh9u3bFZ5z9+5dhIeHg2EsYWm12irPIcRaZwmMAmEsJ2mdJX8ltD6UnDWshPD1+AKV047iww8/xJw5c9CwoWclCsxmMxQKe+6wtcMpr7r7Kt92Bd8pVCDQajXeDsG7MvMARgHExiI0KwuhBXmAtoW3o/K6Cu+LyxcBzgzENgL0ejx6+SKQmOB8I9a2jWsK+GLbuhhf0H9GyhGjLZx2FDqdDr169fJ4R/Xq1UNubq7tdm5urm1Yyqpu3booKCiAyWSCSqWq9jmuoDmKwKbURCKM5RCalYVilkORJhLmIG+Tyu8LJq4lNAolkHXNckYR1xKsK2PVf7WtIvOS5Re7j7WtK/HRZ8ROrDkKpx3F888/jyVLlqB79+62ISEA6NjRvTLHbdq0wdWrV3Ht2jVER0fj+++/x/Dhwys8R61Wo0OHDtizZw8GDRqEnTt3onv37m7thwQ+sy4aRUnJCC3Is3xR0NBDFWx8RxSsWstbH6o61raVq4aVu3w9vkDltKM4efIkfvnlFxw9erTC/bt373ZrRzVq1MDixYvx1ltvobS0FD169EBiYiIAYPbs2UhISECvXr2QkpKC5ORkrFmzBvXr18eKFSvc2g8JDmZdNKBt4VO/dn0NG99R0LoVZl20T38B+3p8gchpUcABAwbgX//6F2rUqCFXTB6joafgUF1byFnRVU7OjoveFxYGgwIFBeHQaAp99poXOck29BQZGQmWZf2qoyDBKVAzYgL1uMRmraXEMADLhvhVLSVf57SjeOyxxzBkyBA8/fTTCAkJsd0/Z84cSQMjxF3lV51T6a9Bpb8eEF+ogXpcYrPWUoqLAzIzLbd1OpO3wwoITjuKhg0bepwaS4gc5Fx1Tk6Belxis9ZSysqC39VS8nW0cJGfobFoO5qjsKP3hQXNUVQk2xzFnj178Omnn+L+/fsV7nd2ZTYh3hCoGTGBelxi0+k4aLVAbi51EmKS7cpsQsSkNNwAMvOgrHQdhRRnFM7WqxaL0NiZjHTg8kUwcS0rxOdPZ1dyVFX1hyrAvkq2K7MJEYujWk9SZAcJrZnkLqGx2+LjzNAolLb4/ClTSo6V3/x5dTlf4HSJJeuV2SdOnEB6errtP0K8xZoFhNhYKFgjVPrrFe43xTSqcL8nmDOnAZaFuYEOYFnLbQkIjd0aH2JiKsQnRVtIRY6V3/x5dTlfINuV2YSIxZoFhKysCllAUmQHsW3bAQwDZbYBYBjLbQkIjd0aH/T6CvH5U6aUHCu/eWt1uUAZ7qIrs/0MZbdYKA03EFGQh3yaowCTkY5HL1/EXZqjAOD4MyL3l7YvDHfRldkkqDmq9SRFdpDQmknuEho7G98RSEyoUh3WnzKldDpO8ovj5NhHeeWHu/R6hV9fAEhXZhPiI/zpDEBqfL/+czIMuHdGjzptYxAVr/NShM55a7hLCnRlNiE+wJ+ylKTGN2STk2FAzptLoWCNyGHUwKrpPttZ6HQckpLKAmKOwmFHcefOHdStWxdvvvmmwxfn5+cjIiJCksAICSZUz8mOb8jm3hk9FKwRZQ0aIST7Gu6d0ftsRwHIP9wlFYc5YrNmzcLGjRurXJENAIWFhVi/fj2Sk5MlDY6QYOGrWUoGgwJpaSoYDO4vSSwU35BNnbYx4Bg1QrIt7VSnbYxscQUzh1lPZrMZX3zxBTZt2oTGjRujUaNGMJvNuH79Oq5evYqXX34Zr7zyCtRqtdwxO0VZT8Eh0NrCkzkKKdrCm1k7nsxRBNr7whNiZT05TY8tKSlBWloarly5AoVCgcaNG1eZ2PY11FEEB2oLOynaIi1NhS1bGNsQ0Esvsejc2feHUeh9YSdbemzNmjXRs2dP9OzZ060ACQkUQn7py/UaALiz9xTuZ1yAKb4V6iZ2cPl1zgRS1g4QOBe/eYPTjoKQYCYkG0mu1wCWTkIzaQKUnAnmL1S4s3qtaJ1FIGXt+MLFb/6MCp4QwkNIzSS5XgMAZUdPQ2kyoSQiBkqTCWVHxa1FpdNx6NzZ5PdfqlTryTPUWoTwEJKNJNdrACCkazuYVSrUzNfDrFIhpKs0taj8XaANo8nN6WR2Xl4evv32W+Tn51e435Mrs9etW4fU1FSEhITgueeewxtvvFHh8R07dmD58uW2azR69uyJqVOnurx9mswODnK1hT/MUagkmKPwV75S68kXyDaZPXXqVNSuXRvNmzeHQuF5LvXx48exe/dupKamIjQ0FJMmTcJ//vMf9O3b1/ac8+fPIzk5GQMHDvR4f4R4SkjNJLleAwB1EztAO/pZ+gHhRKBc/OYNTjuKvLw8fPXVV6Lt8OLFi+jatSvCwy29V7du3bBv374KHcVvv/2GrKwsrF27Fs2aNcPcuXPxyCOPiBYDIe5wVD1W6BmAeu8PCDl6BGVdu8GYOMCl7Tl7TK7V/oTwlTiIcE47Cq1Wi3v37qFOnTqi7LBVq1b44IMPMGHCBISGhuLAgQOoPPql1Wrx6quvon379lixYgUWLFiA5cuXi7J/QtzhaIU7oVlK6r0/oPak8eBMZtT8+is8WL0OxsQBvNtz5TE5VvsTwlfiIJ5x2FG8//77AACVSoW///3vePrppytchS10jqJLly4YNmwYRo8ejTp16qBLly44e/ZsheesXr3a9u9x48ahT58+bu2Db6wtEGi1Gm+H4DMkb4vLFwHODMQ2AvR6PHr5IpCYAGTmAYwCiGsKZGUhtCAP0LZwvr2MNIDjAG0kkJ+POhlpwOgR/Ntz5bHYWISWf0xofGLzUhz0GbEToy0cdhTWM4j4+HjEx8d7vCOrwsJC9O3bF6+88goAYP369YiJsddrKSgoQGpqKsaOHQsA4DgOKpXKrX3QZHZwkKMtmLiW0CiUQNY1yxlFXEuwuQVQaiIRxnJQZF6y/FLWRFZZG6M66vjOqP3FF+By86BQKfEgvjOMTrbnymOhWVkoZjlceajFld1FaKLWoomA+MRinTiWOo7qJqit74tgm7zmawtnPC7hsXXrVowcObLCfevWrcP48eNdib2KzMxMzJgxA6mpqSguLsbf//53LFy40NYZmUwm9OjRA6tXr0abNm2watUq5OTkYMGCBS7vgzqK4CBXW/jDHEVEQR4uPNRi0eYmtovKZo66Ap3xmuxzA5UvbpMqDkcX0Wm1Gpw5UxhUF9jxtYWkWU9ff/01SkpK8OWXX6K0tNR2v9FoxLZt2wR3FM2bN0ffvn0xePBgmEwmjB07FvHx8Zg9ezYSEhLQq1cvfPzxx5g3bx5KSkoQGxuLpUuXCtoXIWJwtMKd0CwlY+KACh2EK9tz9hi0LXBld1GF8txXjA0R1Vn+EtyVy4RLFQdfOfJAWl3OFVIfr8OOgmEYXLp0CSUlJbh06ZLtfpVK5XF58UmTJmHSpEkV7lu4cKHt3x06dMCOHTs82gch7hCyLrbQtbQdnR14mh1U+aKyJurrUKe5/0te7DikuriNbz/BdoGd1MfrdOhp37596N27t6g7lRoNPQUHsdrCUWaT2K8BHGcBeZodVHlcvon6OppsXuT29sTKUpJrfoDmKOyknKNwWsKjc+fOWLp0KQYPHoyhQ4di5cqVKCsrcyN8Qnwbc+Y0wLIwN9ABLGu5LcFrAMc1nYTWeqrMVpvJeE3Q9kSPQ+IvaL79BEqdKldJebxOO4q5c+fi9u3bmDlzJqZNm4Y///zTljpLSCBg27YDGAbKbAPAMJbbErwGcFzTSewV7oRuz1dX2iPe5XToqV+/fvjpp59st81mMwYMGIAff/xR8uCEoqGn4CBmW/j7HEV1bSF0e/5+JTV9Ruxkq/UUFRWFO3fuoG7dugCAoqIiPProo26ESojvc5TZJPZrAMcZTEKzqNzdj1SvI4HLaUdRr149DB8+HImJiVCpVNi/fz8iIyNtw0+eVJEl/seff206un5BKL7JUrEnUvnOXgwGBTIzAY1GUWFfQq/LEPI8EticdhSNGjVCo0aNbLcHDPD8A0b8kz/X7XFUY0kovhXTxF5NjS/DyrovhgFYNsS2L6G1o8rz5783EZfTjuLNN99ESUkJrl27hieeeAKlpaUIDQ2VIzbiY8pnxKj016DSX/ebL46Qo0fAmcxA3brg7txByNEjHnUUcl7sVT7DSpltsJxZ/NVRWPcVFwdkZsK2L76/lat/R3/+exNxOc16Onv2LHr37o0JEyYgJycHPXv2xK+//ipHbMTH+HNGTFnXblColMCdO1ColCjr2s2j7cl5sRdfhpV1X1lZqLAvvr+Vq39Hf/57E3E5zXoaOXIkFixYgKSkJOzcuROHDx/Gp59+itTUVLlidBtlPUnHl8as3W2LQJ6jKCgIh0ZTSHMU8P5nxJfIlvVUUlKCuLg42+0ePXrgo48+cjFMEmj8OSPGUY0lofhWTBN7NTW+DCtL8TcgN7dihyS0dpSQ55HA5rSjYBgG9+/fty2DeuXKFcmDIsQZR6u6SbUvb/+q5ovBUdaTFIScKYndfkLP1oKtpIeYnHYUb7zxBkaNGoW8vDy88847OHbsmFslvwkRm6NV3aTclzczf/hicJT1JAUh2Vxit5/QjDKxM9GCjdPJ7GeffRarVq3CW2+9hfbt22Pr1q3o16+fHLERUi1rNg5iYz2qR+TOvjytfSRVDNasp9hYS/asXu/0Iy1Y+WwuV/cldvsJicGT1xELh2cU2dnZtn+r1Wr07NmzwmMNGjSQNDBCHLFm4yAry61sHCFDIL6Q+cMXg6OsJykIyeYSu/2EZpQFW9lxsTnMemrXrh0UCgU4jkNJSQlq1aoFlUqFBw8eICIiAkePHpU7VpdR1lPgs67qlu/iHIUnQyD+MEdRXdaTFPxhjsLRZyQY5ygkz3o6fdpSNvndd99Fp06dbFdk79+/H/v27XM3XkJEZV3VzdX1lz25eMwXMn/4YnCU9SQFIdlcYref0IwysTPRgonTgbrz589XKNvRq1cvZGZmShoUIYDll6g67bglw6max3DkSJXHHL1GiiEkvvjEfp3BoEBamgoGg6Lax44cQbWPOZKRocSGDWpkZMg3Vs93DK48biV27DkZBlzakIacDIMo2wtETrOezGYzTp48iU6dOgEAfvnlF1uqLCFScaVWUeWsJ77XmHXRKEpKFm0IROhQlpDXuVJXyp2sp4wMJd58s+Zf21Nj1aoSxMdLO2bvLOvI1awksWPPyTAg582lULBG5DBqYNV0RMXLv864r3PaJc+ZMwf/+7//i2effRY9e/bE3LlzKT2WSI4vW8ZR1pOzDBuzLhrGzk+LMgwiNJtHyOv4MnaEZD2dOaMCywINGli2d+aMyqXYPeEs68jVrCSxY793Rg8Fa0RZA8vf494ZvUfbC1ROzyg6dOiAgwcP4tKlSwCAZs2agWGcvowQj7hSq6hy1pOcGUpyriDnSl0pa9aTWs0hLU3FW1ZEp7O8Jjvbsr22baUft3eWdeRqVlLbtiYwjNqj2MtPatdpG4McRo2QbMvfo07bGEHHF+ic1nqSwrp165CamoqQkBA899xzeOONNyo8np2djWnTpiE/Px+NGzfGsmXLUKtWLZe3T1lPgcFZraLqsp7kzFCScwU5Z3WlCgrC8fDhQ2zerHap9Hm/fiwMBiXatjVJPuzkyjG48rhVRoYSZ86oHMbO9xmpbohLfesG7p3Ro07bmIAbdhIr60n2q06OHz+O3bt3IzU1FTt37sTZs2fxny+EkAgAAB3YSURBVP/8p8Jz5s+fj5EjR2Lv3r1o3bo1PvvsM7nDDBh8E6dMRjpqblgHJiNdlO3JyayLBrp1q/JFK+bwkpWjSVY9YnAE3aCHe79ChcSo03Ho3NlU7ReoTsehWzfAaFQ4HaKyPlanDvDaa0bUq8e5NIFs5eqEc3nWyWL1rRtVjqH8+4nvGMurV49Dq1Zm1Kvn/o/B6oa4ouJ1aPpa54DrJMQk+xjSxYsX0bVrV4SHW3qvbt26Yd++fejbty8AwGg0Ij09HatXrwYADBs2DKNGjcK0adPkDtXv8U2c8i2GI2R7csYuJ0eTrL5YEsLd0ufuHoOQY+abLBZ7Yt/TNiKOOe0ozpw5g7Zt29pul5aWYvHixUhJSRG0w1atWuGDDz7AhAkTEBoaigMHDqD86Nfdu3cRHh5umwfRarW4ffu2W/vgO4UKBFqtxrUnZuYBjAKIawpkZSG0IA/QtrA8dvkiwJmB2EaAXo9HL18EEhOEb09sLu7L5bYQGkam5QspLs4yD1BQEAKt1vH93tS2bTgWLbLEExsLxMSE2B7TalHlsSNH3DsGIcecfTkXDGeGKfZxMPossJdzoU1s/tcG3X8/uRqDo/dFde0Q6MT4jDjtKKZOnYp169bhiSeewLlz5zBt2jQ8/vjjgnfYpUsXDBs2DKNHj0adOnXQpUsXnD171vY4x3FV0m/dTcelOQoLpSYSYSwHReYlyy82TaTtAjUmriU0CiWQdc1yRhHXEqyT7So1kah1vxCKI0fB1a6Dh+W2Jza+2K3kmK/RaBRg2RDbF5RGU4bcXM7h/d5ibYuaNYHmf30P5+ZWfE7lx9w9BiHHzMRpwSqUUGT9CZZRg4nT2v5mrvyNhcTg7H3B10aBRqw5CqeT2adOncKsWbPQu3dv7NixAzNmzMDzzz/vfsR/KSwsxP3796HTWU4/169fj1u3bmHOnDkALENPnTp1Qnp6OlQqFW7evIlRo0Zh//79Lu+DOgo7volTvsVwHG2r1vy5UDy4D672I3iY8p6kw0HOJn3lmth3NMnqSyUhhLaFu8cg5JhzMgwOJ4vFntgHgivhwxnZFi7q0KED3nvvPbz++uvYuHFjhWEoIW7cuIEZM2YgNTUVxcXF+O6777Bw4ULb42q1Gh06dMCePXswaNAg7Ny5E927d/don8GMr3wC32I41VHprwM1QsC27yDLGsq+UDoDcFz6IRBKQrh7DEKOOSpe53CiWMjfOBDa3d847CgGDRpU8YkMgzfeeAORkZEAgN27dwvaYfPmzdG3b18MHjwYJpMJY8eORXx8PGbPno2EhAT06tULKSkpSE5Oxpo1a1C/fn2sWLFC0L6IuEwxDYHSMjC/ngJX+xHRymAIOeO5s/cU7mdcgCm+FeomdvA4Dj58v4h9YT9yLlwk5AxA7DMKIj+HQ0///e9/eV/41FNPSRKQGGjoSRr2oad7ljkKD4eehGZl3dl7CppJE6DkTDArVChYvVayzqJ81g7HqBElUYkHofuxl/AIAcuWSZp9JSRLie+4pMpso6EnO8mvo3jqqads/zVt2hQxMTGIjo5G/fr1YTQahUVN/Jp96KkjUCPE40Vo+MpZMGdOAywLcwMdwLKW238pO3oaSpMJJRExUJpMKDt6urrNi0KuEg9C9yPnwkVCyo/wHZcvLApFXON0juKTTz7BunXrAAAqlQpGoxFxcXGCh56I/xK7RAbf9ti27QCGgTLbADAMTDod1GnHYYppiJCu7WD+WoWa+XqYVSqEdG3n6aE5JFeJB6H7kXPhIiF/f77j8oVFoYhrnGY9JSQkYNu2bVi8eDGmT5+OtLQ0HD58GB999JFcMbqNhp6kI/aYsitzFCadDjV++rHCEEXeb7egojkKAPIuXOQPcxTe/oz4EtmynurWrYuoqCg0adIEmZmZeP755/HPf/7TvWhJwBA7E8mVrCx12vEqiw7VTXwa2tHPyvKFwJe14wv7kXPhIiF/f7Gznoj8nA5oMgyD69evo0mTJjh16hRYlkVpaakcsQUsX6mX5AtcaQt3hijkbFu+fQmpiSTUli0M/v53y//FQO9PUpnTd9aECRMwd+5crFmzBp988gl27tyJnj17yhBaYPKVGka+wNW2cHXRIV+pRSVnHagtWxhMm1YDALBzp+X/L73ECt4evT9JdZyeUTz77LPYtGkTwsLCsHPnTqxfvx7vvfeeHLEFJMr0sHOnLVypuCpn2/Lty9VFeMRw4IAKHAfUqgVwnOW2J+j9Sarj9B388OFDzJ8/H2PGjEFpaSm2bt2KoqIiOWILSJTpYSdnFpXY+PYlZ4XShAQTFArg4UNAobDc9oQ7bSjn8BrxLqdZTzNnzkRUVBT279+Pb7/9FrNnz4ZCocDy5cvlitFtvp715EmmR6BldIjdFr6ycJGcdaC2bGFw7Fgonnmm2KNhJytX2tAXy6xbBdpnxBOyZT39/vvvWLRoEQ4fPozQ0FAsW7YMAwcOdC9aUgFletjJmUUlNr59yVmP6KWXWEyZAuTmet5JAK61YfnhNb3eusQq1V8KVE47CqWy4uiUyWSqch+Rh9JwA8jMg7LS8p/ETs4zCrGvD+Bb4tPXaiLRAkDBxWlH0bFjR3z44YcoKSnBkSNHsGXLFnTq1EmO2Eg51mwUMAqEsRxlo1RDzowdsVduy8hQ4s03a/41lKPGqlUlts7CFzORdDoOSUllPlNmnUjL6alBUlISwsLCoNFo8NFHH6FZs2aYPn26HLGRcqzZKIiNpWwUB+TM2BG7htGZMyqwLNCggSVT6swZe/aSr2YiubrGNfF/Ts8o1Go1Jk2ahDFjxkCtVqNGjRpyxEUqsWajICsr6LOlHJEz60nsGkZt25rAMGpkZ1uGctq2tY/3V96eQd0IV9JU9EueyMZp1lNWVhamT5+OCxcuQKFQoH379liyZAnq168vV4xu8/WsJ6GUhhuIKMhDPs1RAPB+1pM35igM6kZYtLlJlWwjyvSxo7awk20p1JdffhkDBw7E0KFDwXEcvvnmGxw6dAgbNmxwP2qZBGpHAdCHoLxgbIu0NBW2bGFs2UYvvcSic2dTULaFI9QWdpKvR2H14MEDvPDCC1Cr1QgJCcHo0aORl5fnXrQkYIhdB4jvoi2+x3IyDDi76ghyMgwuv8bXuRJ7ddlGBoMCR47AL4+Z+AencxQNGzbE2bNn0aZNGwBAZmYmGjak8fFgJHb2Dd9FW3yPWTOOGM4MVqG0ZRz58kVgzrgae+VsIwB/rXAHsGyIXx0z8R9O18x++PAhRo4ciWbNmkGpVCIzMxOPP/64bAES31E++8Za7tuTjoLvoi2+x6wZR6bYx6HI+hP3zugRFa/z64vA3Im9/MV8aWmWbKm4OCAzE351zMR/OOwo5s6dK2ccxA+InVXEd9EW32PWjCNGnwW2XMaRP18EJjR2OVe4I8HL6WS2P6LJbOmInVXEVxOJ77GcDAPYy7lg4rQVMo7krLEkNqGxy7nCnT/w9mfEl8iW9SSFXbt22dbh7t69O2bMmFHh8R07dmD58uWIiIgAAPTs2RNTp051efvUUQQHags7ags7ags72YoCiq24uBgLFy7E3r17Ubt2bbz44os4fvw4nn76adtzzp8/j+TkZCo+6EVyXY8g9Fe0o7pXcl5H4StnL0JqgPFdsyE2X2knIpzsHYXJZILZbEZxcTHCwsLAsmyVq71/++03ZGVlYe3atWjWrBnmzp2LRx55RO5Qg5ZctYWEZik5qnslZ00kX8mwElIDjK+ulNh8pZ2IZ2TvKMLDwzF58mT0798foaGh6NixI9q3b1/hOVqtFq+++irat2+PFStWYMGCBW6tf8F3ChUItFqNtDvIzAMYBRDXFMjKQmhBHqBtIf5uMi1fHnFxlsnYgoIQaLVuxBcbi9Dy8ckUt0exix6Ig7bgcfmyZTW82FhArwcuX66FxESJwvNSO0n+GfEjYrSF7B1FZmYmUlNTcfDgQWg0GiQlJWHDhg0YN26c7TmrV6+2/XvcuHHo06ePW/ugOQrPKDWRCGM5KDIvWX6ZayJhlmCfGo0CLBti+zLRaMqQm+vCGcVf8YVmZaGY5WzxyRW3J7GLzVFb8ImLU0KhqGnLlIqLK0FurjRnFN5oJ5qjsPPbyez169cjPz/fNoF96NAhbN261Ta5XVBQgNTUVIwdOxYAcO/ePfTv3x8nTpxweR/UUXjOH+Yoqqt7FaxzFO7WAAvkOQrqKOz8djK7efPm+PDDD1FUVITQ0FAcOHAATz75pO3xsLAwrF+/Hu3atUObNm2wefNmt88oiOfkWilO6EpwZl00oG1R5deznCvcybmKHR9HbcEnPt4seQdh5SvtRISTvaPo2rUrLl68iGHDhkGtVuPJJ5/E+PHjMXv2bCQkJKBXr174+OOPMW/ePJSUlCA2NhZLly6VO8ygwPfrW8gvczmvsWAy0oHLF8HEtQQb39Er8cnJ2fUmmZmWYR6aKBaXr5w1ehtdcOdnxDqt5ssQEpI9JGcdKCYjHZo3J9hqPRWsWuu0s/DFVeJc5UpNLIYJAcuWUVYRxPuMBELGlmzVY0lg4ls1TciKamKvwla+9hHLWm5bMWdO468HAZa13JY5PjnxtYX1sdhYVHmMeIav3YNN8B55kOOr2ySkppOcdaDYtu3w14MAw1huuxmfQd3Ib8qRu1ITi2o9ic+fa4eJjYae/IyYGR3+Pkfx6OWLuCtgjsLRKnG+zNkcBdV6shPzM+LvcxR+mx4rB+oogoPQtnC0Spw/o/eFHbWFnd+mx5Lg5EsZR64OKZQ/A7libCjar0ohv1L9/ZetUMF63L6GOgoiOV/LOKq8Slx1X0DWmIvus8j5vQb2tpyN+7VjPB6mEpJJEwjZN0IE63H7IprMJpLzxYwjnY5D584mh1881pjzwmOhYI1oVStLlMwXIZk0wZp9E6zH7Yuo5YnkxM6IkoM15sjCLHCMGhcexoqS+SIkkyZYs2+C9bh9EU1m+xl/naiTYo5C6rbwpzkKf31fOCOkrQK1LYSgyWziV+SswSQWa8xRAKIgXlaUkNpHwVovKViP29fQ0FMQUxpuQJ123LJCmp/ty2BQ4MgRyHLBnMGgqPbiPEf3ExJo6IwiSMmZiSRVHSiGAVg2RNJsGEeZN5SRQ4IJnVEEKTkzkaSqAyVHfSNHmTeUkUOCCb27g5ScmUhS1YGSo76Ro8wbysghwYSynvyMXLWexCZFHSi56hs5yrzxpauGKdPHjtrCjrKeiMfkzEQSe186HQetFrKsU+0o84YyckiwoKEn4rOEZErJmYkkZ9YYId5EZxTEJwnJlJIzE8nX6lcRIiU6oyA+SUimlJyZSL5Yv4oQqVBHQXySkEwpOTOR/LF+FSFCUdaTnwmmjA5nmVLVtYWcmUi+tMZGML0vnKG2sPPrrKddu3Zh3bp1AIDu3btjxowZFR7Pzs7GtGnTkJ+fj8aNG2PZsmWoVauWN0IlXiQkU0rOTCR/rF9FiBCyDz0VFxdj4cKF+Oqrr7Br1y6cOnUKx48fr/Cc+fPnY+TIkdi7dy9at26Nzz77TO4wfZLScAM4cqRKlo0/Z98IzVJy1BZyEho71Ygi/kb2MwqTyQSz2Yzi4mKEhYWBZVnUqFHD9rjRaER6ejpWr14NABg2bBhGjRqFadOmyR2qT7Fm2YBRIIzlbFk2/px9IzRLyVFbyElo7FQjivgj2TuK8PBwTJ48Gf3790doaCg6duyI9u3b2x6/e/cuwsPDwTCW0LRaLW7fvu3WPvjG2vxWZh7AKIDYWIRmZSG0IA/QtrDfH9cUKH+/H8jMtHxZxsVZynEUFIRAq3XlhQ7aQkZCYxd8zC7QajXibCgAUFvYidEWsncUmZmZSE1NxcGDB6HRaJCUlIQNGzZg3LhxAACO46BQVDwlr3zbmUCczFZqIhHGcgjNykIxy6FIEwlzboHtfkXmJcsZxV/3+wONRgGWDbF9eWo0ZS5dae2oLeQkNHahr3OGJnDtqC3sxJrMlj3raf369cjPz7dNYB86dAhbt261TW4bjUZ06tQJ6enpUKlUuHnzJkaNGoX9+/e7vI9A7CgAy5BLREEe8jWRFYZafCn7xl1Cs5QctYWchMYuRWYWfTnaUVvY+W3WU/PmzfHhhx+iqKgIoaGhOHDgAJ588knb42q1Gh06dMCePXswaNAg7Ny5E927d5c7TJ9k1kUD2hZVfj37c/aN0CwlR20hJ6GxU40o4m9kz3rq2rUrBgwYgGHDhmHw4MFgWRbjx4/H7NmzbWcNKSkp2L59O5577jmcOnUKU6ZMkTtM4gP4soN8IeuJkGBBF9z5mWA5rebLDrJmPYUyCsschR9lekklWN4XrqC2sBNr6IlKeBCfxFe3yVpnCbGxVGeJEBlQR0F8El/dJmudJWRlUZ0lQmRAZcaJT9LpOCQllVWbHWTWRaMoKRmhBXmW1NggH3YiRGrUURCfxZcd5AtZT4QECxp6+os/10sihBAp0RkFaLUyQgjhQ2cUoNXKCCGED3UUoNXKCCGEDw09wZ5F46/1kgghRErUUfzFn+slEUKIlKijIAGFyUgHc+Y02LbtwMZ39HY4hAQE6ihIwGAy0qF5cwKsBaIKVq2lzoIQEdBkNgkYzJnTAMvC3EAHsKzlNiHEY9RRkIDBtm0HMAyU2QaAYSy3CSEeo6EnEjDY+I4oWLWW5igIERl1FCSgsPEdqYMgRGQ09EQIIYQXdRSEEEJ4UUdBCCGEF3UUhBBCeFFHQQghhFdAZj0plQpvhyCpQD8+d1Bb2FFb2FFb2LnSFs6eo+A4juN9BiGEkKBGQ0+EEEJ4UUdBCCGEF3UUhBBCeFFHQQghhBd1FIQQQnhRR0EIIYQXdRSEEEJ4UUdBCCGEF3UUhBBCeFFH4eMKCwsxcOBA3LhxAwDwzTffYODAgRg0aBBmzpyJsrIyL0con8ptsXXrVgwYMADPPfcclixZgmAqMlC5Law2b96M0aNHeykq76jcFjNnzkTfvn0xZMgQDBkyBD///LOXI5RP5bY4ffo0XnjhBQwYMADvvPOO4O8L6ih82NmzZ/Hiiy8iKysLAHD16lVs2LAB27Ztw7///W+YzWZs3brVu0HKpHJb6PV6fPnll/j222+xe/dunD59GseOHfNukDKp3BZWly9fxrp167wTlJdU1xbnz5/H5s2bsWvXLuzatQt9+vTxXoAyqtwWhYWFeOutt7BgwQL88MMPAIDvvvtO0Lapo/Bh27dvR0pKCqKiogAAISEhSElJQXh4OBQKBZo2bYrs7GwvRymPym0RExODH374AWFhYXjw4AEKCwtRu3ZtL0cpj8ptAQBlZWV499138fbbb3sxMvlVbovi4mJkZ2dj1qxZGDRoED799FOYzWYvRymPym1x7NgxtG3bFs2bNwcAzJkzR3CnGZDVYwPFwoULK9zW6XTQ6XQAgDt37mDLli1YtGiRN0KTXeW2AAC1Wo3t27djyZIl+J//+R/bByLQVdcWy5cvx/DhwxEdHe2FiLynclvk5eWhc+fOSElJgUajwYQJE/Ddd9/hhRde8FKE8qncFteuXUNYWBimTp2KK1euoH379khOTha0bTqj8EO3b9/GmDFjMHz4cHTq1Mnb4XjVCy+8gJMnTyIyMhKrVq3ydjhecezYMdy8eRPDhw/3diheFxMTg9WrVyMqKgqhoaEYPXo0Dh8+7O2wvMJkMuHo0aN455138K9//QvFxcWChyapo/Azf/75J0aMGIGhQ4di0qRJ3g7Ha27evImMjAwAAMMwGDBgAP744w8vR+Ud33//Pf7v//4PQ4YMwZw5c3D+/HlMmTLF22F5xR9//IGffvrJdpvjODBMcA6cREZGok2bNoiJiYFKpUL//v1x7tw5QduijsKPFBYW4rXXXsPkyZPx6quvejscryooKMC0adPw4MEDcByHn376CfHx8d4OyysWLVqEH3/8Ebt27cL777+P1q1b4+OPP/Z2WF7BcRw++OAD3L9/H0ajEd98803QTGZX1rVrV1y4cAE3b94EABw8eBCtWrUStK3g7Gr91HfffYe8vDxs3LgRGzduBAAkJCRg8uTJXo5Mfk2bNsX48eMxYsQIqFQqdOjQAa+88oq3wyJe1rx5c4wfPx4vvvgiWJZF3759MXDgQG+H5RX169fHggUL8Prrr6O0tBQtWrTAjBkzBG2LVrgjhBDCi4aeCCGE8KKOghBCCC/qKAghhPCijoIQQggv6igIIYTwoo6C+J2TJ0+KmvLo6vYSEhLw22+/ibbf8l599VXcuXPH7f38/vvvmDlzpiQxLV68GCdPnpRk28S/UEdBiA8QUvnWbDZj9uzZkl2FPWnSJLz//vsoKSmRZPvEf9AFd8QvFRUV4e2338a1a9dQu3ZtLFiwAI0bN0ZZWRmWLVuG9PR0mEwmtGzZEnPmzEF4eDgOHjyItWvXoqysDHfu3MHzzz9f5Uv21KlTSEpKwooVK9C+fXuH+z9w4ADWrFkDo9GImjVrYsaMGWjXrh1WrlwJg8GA3NxcGAwGPPbYY/jwww8RFRWFc+fOYd68eTAajWjYsCGys7ORnJyMnTt3AgDGjBljq8XzzTffICUlBXfu3MGQIUMwderUKjH8+OOPiI6OxmOPPQbAciYycOBApKWl4f79+xg3bhx+/fVXXLhwAQzDYM2aNXjsscdcfp5Go0G7du3wzTffYMyYMWL96Yg/4gjxM2lpaVzz5s25jIwMjuM4btu2bdzf/vY3juM4buXKldzixYs5s9nMcRzHLV++nEtJSeHMZjM3atQo7urVqxzHcdytW7e4Fi1acPn5+VxaWho3YMAA7sSJE1zv3r2533//vdr9Pvvss9y5c+e4q1evcgMHDuTu3LnDcRzHXbp0iXvmmWe4hw8fcp9++inXq1cvrqCggOM4jpswYQL3ySefcEajkevevTt36NAhjuM47sSJE1yzZs24tLQ0juM4rmnTplx+fr5tPwsWLOA4juNycnK41q1bc9nZ2VXieeutt7jU1NQK8X3wwQccx3HcDz/8wDVv3tx2LBMnTuTWrFnj1vM4juN++ukn7qWXXnLlz0ICGJ1REL/UrFkz2y/+oUOHYt68eSgoKMChQ4dQUFCA48ePAwCMRiMiIiKgUCjw+eef49ChQ/j+++/x559/guM4FBcXAwBu3bqF119/HS+++KLTcuXHjh1DTk4Oxo4da7tPoVDg+vXrAICnnnoK4eHhAICWLVvi/v37uHTpEgCgR48eAIDOnTvjiSeecLgP65yJVqtFZGQk8vPzUb9+/QrPuXLlCl5++eUK9/Xt2xeApYpqZGSk7VgaNmyI+/fvu/286OhoXL16lbc9SOCjjoL4JaWy4vSaQqEAwzAwm82YNWuW7Qv54cOHKC0tRVFREYYOHYrevXujQ4cOGD58OPbt22dbPlWlUmHdunWYOHEiEhMT0aZNG4f7NpvN6NKlS4XCezdv3kRUVBR+/vln1KxZs0JcHMdBpVJVWapVpVI53Ef5iqfWbVRW3f0hISG2f6vVaofbd/V5DMNUaWsSfOgdQPzSH3/8gd9//x2AZTw/Pj4eoaGh6Nq1K7Zs2YKysjKYzWbMnTsXK1aswLVr11BYWIgpU6YgISEBJ0+etD0HsPxyb9++PWbMmIHp06fbzjSq06VLFxw7dgx//vknAODw4cMYPHgw76Tv448/jpCQEPzyyy8AgHPnzuHSpUtQKBQALJ0Gy7JutUHjxo1tZzFSuXHjBpo0aSLpPojvozMK4peaNGmCVatWQa/XIyIiAosXLwYATJw4EUuWLMHQoUNhMpnQokULJCcnIywsDD179kT//v0REhKCpk2bIi4uDteuXavw63ro0KH46aefsHjxYsyfP7/afcfFxWHBggV45513bOsdrFmzBrVq1XIYL8MwWLlyJVJSUrBixQrExsYiMjLSdvaRmJiI0aNHY+XKlS63Qb9+/fDzzz9LumDRkSNHkJiYKNn2iX+g6rGEyGTJkiV47bXXEBkZiZs3b2LIkCHYt2+f4LW+TSYThg0bhnXr1tkyn8RUWFiIESNGIDU1FTVq1BB9+8R/UEdBiEw2b96Mbdu2gWEYcByHSZMm2SaVhTp37hy2bNmCJUuWiBSl3aJFi9CjRw88/fTTom+b+BfqKAghhPCiyWxCCCG8qKMghBDCizoKQgghvKijIIQQwos6CkIIIbyooyCEEMLr/wHnmj5To5tAMwAAAABJRU5ErkJggg==\n",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Make scatter plot of 1975 data\n",
"_ = plt.plot(bl_1975, bd_1975, marker='.', linestyle='None', color='blue', alpha=0.5)\n",
"\n",
"# Make scatter plot of 2012 data\n",
"_ = plt.plot(bl_2012, bd_2012, marker='.', linestyle='None', color='red', alpha=0.5)\n",
"\n",
"# Label axes and make legend\n",
"_ = plt.xlabel('beak length (mm)')\n",
"_ = plt.ylabel('beak depth (mm)')\n",
"_ = plt.legend(('1975', '2012'), loc='upper left')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Linear regressions\n",
"Perform a linear regression for both the 1975 and 2012 data. Then, perform pairs bootstrap estimates for the regression parameters. Report 95% confidence intervals on the slope and intercept of the regression line."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"def draw_bs_pairs_linreg(x, y, size=1):\n",
" \"\"\"Perform pairs bootstrap for linear regression.\"\"\"\n",
"\n",
" # Set up array of indices to sample from: inds\n",
" inds = np.arange(len(x))\n",
"\n",
" # Initialize replicates: bs_slope_reps, bs_intercept_reps\n",
" bs_slope_reps = np.empty(size)\n",
" bs_intercept_reps = np.empty(size)\n",
"\n",
" # Generate replicates\n",
" for i in range(size):\n",
" bs_inds = np.random.choice(inds, size=len(inds))\n",
" bs_x, bs_y = x[bs_inds], y[bs_inds]\n",
" bs_slope_reps[i], bs_intercept_reps[i] = np.polyfit(bs_x, bs_y, deg=1)\n",
"\n",
" return bs_slope_reps, bs_intercept_reps"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"def bootstrap_replicate_1d(data, func):\n",
" \"\"\"Generate bootstrap replicate of 1D data.\"\"\"\n",
" bs_sample = np.random.choice(data, len(data))\n",
" return func(bs_sample)\n",
"\n",
"def draw_bs_reps(data, func, size=1):\n",
" \"\"\"Draw bootstrap replicates.\"\"\"\n",
" \n",
" # Initialize array of replicates: bs_replicates\n",
" bs_replicates = np.empty(size)\n",
" \n",
" # Generate replicates\n",
" for i in range(size):\n",
" bs_replicates[i] = bootstrap_replicate_1d(data, func)\n",
" \n",
" return bs_replicates"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1975: slope = 0.4652051691605937 conf int = [0.32492037 0.59321466]\n",
"1975: intercept = 2.3908752365842276 conf int = [0.57667688 4.3381558 ]\n",
"2012: slope = 0.46263035883531306 conf int = [0.32689249 0.59623919]\n",
"2012: intercept = 2.9772474982360198 conf int = [1.14976518 4.79237901]\n"
]
}
],
"source": [
"# Compute the linear regressions\n",
"slope_1975, intercept_1975 = np.polyfit(bl_1975, bd_1975, deg=1)\n",
"slope_2012, intercept_2012 = np.polyfit(bl_2012, bd_2012, deg=1)\n",
"\n",
"# Perform pairs bootstrap for the linear regressions\n",
"bs_slope_reps_1975, bs_intercept_reps_1975 = draw_bs_pairs_linreg(bl_1975, bd_1975, 1000)\n",
"bs_slope_reps_2012, bs_intercept_reps_2012 = draw_bs_pairs_linreg(bl_2012, bd_2012, 1000)\n",
"\n",
"# Compute confidence intervals of slopes\n",
"slope_conf_int_1975 = np.percentile(bs_slope_reps_1975, [2.5, 97.5])\n",
"slope_conf_int_2012 = np.percentile(bs_slope_reps_2012, [2.5, 97.5])\n",
"intercept_conf_int_1975 = np.percentile(bs_intercept_reps_1975, [2.5, 97.5])\n",
"intercept_conf_int_2012 = np.percentile(bs_intercept_reps_2012, [2.5, 97.5])\n",
"\n",
"# Print the results\n",
"print('1975: slope =', slope_1975, 'conf int =', slope_conf_int_1975)\n",
"print('1975: intercept =', intercept_1975, 'conf int =', intercept_conf_int_1975)\n",
"print('2012: slope =', slope_2012, 'conf int =', slope_conf_int_2012)\n",
"print('2012: intercept =', intercept_2012, 'conf int =', intercept_conf_int_2012)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Displaying the linear regression results\n",
"Now, you will display your linear regression results on the scatter plot."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEJCAYAAACdePCvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOy9eZQl113n+bk34u37npkvt1pUKlXJLi0ILMu2wAZjsww+A8yYweCGA9gY02CO2jZGbbkbjA2oD3iaw9AwB9ttMGDM9Aymbdmydll77fuSlfueb9/fi7h3/ghVaUGlzKrKrKpUxeccHdV7+SLiF/Ey7y/it3x/QmutcXFxcXG57pBX2wAXFxcXl6uD6wBcXFxcrlNcB+Di4uJyneI6ABcXF5frFNcBuLi4uFynuA7AxcXF5TrFdQAuLi4u1ynm1TbgYimVGih18a0LqVSYQqG+ARZtDJvJ3s1kK2wuezeTrbC57N1MtsKl2yulIJEIvebPNp0DUEpfkgM4t+1mYjPZu5lshc1l72ayFTaXvZvJVlh/e90QkIuLi8t1iusAXFxcXK5TNl0I6LVotRrU62Vs27rgZ5aWJEqpK2jV5bEWew3DJByOEwi8dnzPxcXF5fXY9A6g1WpQq5WIxzN4PF6EEK/5OdOUWNbmcQCr2au1ptfrUi4vA7hOwMXF5aLZ9CGger1MPJ7B6/VdcPF/IyKEwOv1EY9nqNfLV9scFxeXTcimdwC2beHxeK+2GVcNj8f7uqEvFxcXlwux6R0AcF3d+b+a6/ncXVyuCywLUdmYp/xNnwO4Vmk06nz4w7/MH//xn9HfP8A3v/kNvvrV/46Ukttuu4OPfvS3qdWqfOxjH33FNuVyiQcffIJ9+/byyU/eQzabA2DHjhv51Kfuu1qn4+LicqXRGrGygrAtVDqzIYdwHcAGcPToEf74j/+A6ekpAKamJvjrv/6/+Ou//u+k02nuv//zfP3r/8D73/8BvvSlrwKglOK3fuvX+dVf/QgAx48f4+d+7gP8wi/80lU7DxcXl6uDKBURrRYqmUL7/Rt2nDdECOhSmZ0VPPOMwezs+oZRvvGN/8Hv/M4nSL/otc+cOcPu3W8inU4DcNddb+OJJx57xTbf/Oa/4Pf7efe73wPA8eNHee65Z/jgB9/PJz7xMRYXF9bVRhcXl2uQeh05O4P2eFEDedjAxR+uYwcwOyu4/34vf/d3Jvff711XJ/DJT/5H9uy59fzr7dtv4NixIywuLmDbNo888hCFQuH8z23b5stf/iIf/vBvnn8vEonwMz/zv/PlL/8Dd955F5/5zKfWzT4XF5drjHYbOTeLsHqo/CCEw1fksNetA5iellgWDA1pLMt5vVEMD4/w4Q9/lE9+8nf4jd/4VbZvvwGP56Xo27PPPs3Q0BDbtm0//94nPvF73H33OwF43/t+hvHxs9Trm0e4ysXFZQ1YFnJhHlGrofoH0PHEFT38desAhoYUpgnT0wLTdF5vFJ1Oh5tu2s0Xv/hV/vIv/4Z0Oks+P3j+50888Sjvete7z79WSvHFL/7f2Lb9iv0YhrFhNrq4uFxBlEIsLSELK6hMFp3JwFWo6LtuHUA+r7nnni4///MW99zTJZ/fOFXAdrvFb/3Wr9NsNuj1evzzP/8j73znSwv+kSOHXhEyklLy2GOP8OijDwPwrW/9K7t23UwgENgwG11cXK4MolhALsyjYzFUrg9WubETlTLiZSHj9eS6rgLK5zX5vL36By+TWCzOL//yr/Jrv/ZLWJbFj/zIj55P9gLMzc2SyWRfsc2nP/2f+cM//H2++MW/JpFIcO+9/2nD7XRxcdk4RK2KqFZR8QQ6mVp9g2YTWSqiYzF0LL4xNmmtN5UgdqFQf4Um9sLCJH19I6tu90bTAno5a70GG0UmE2F5uXbVjn+xbCZ7N5OtsLnsvWK2tlrIYgEdDq9tIe/1kCvLaJ8PnUieDw1dqr1SClKp104qX9dPAC4uLi4bxrmF/FxJ52oxfqUQKyugNSqbc0JDWiNKRajXIbN73U10HYCLi4vLeqIUYtlR6T2/kK+CKBUR7TYqlQavo20mKmXk2bOI6Sl0Lge3uQ7AxcXF5dpEa0SxiOh2XrGQvy71OrJSdvICiSQAol5Dnj2LnJpAeTwwMIDKD22Iya4DcHFxcblMRLXi1PInU+jUGhK8nQ6ysIIOBp3GL4BmE2PiLMbEOEpIdDaHyvWj+/sRzcaG2O06ABcXF5dL5VylTjT60kL+eti2kxcQEtXXD1I6XcCTExhjZxBaoRJJ7P48enAQMT+HcewoKhYDhtfdfNcBuLi4uFws3a5zB+/zrS3Bey481Os64SGPx9nH5ARy/CzCstDhMHZ+EDU8gliYxzxyCOXzO06mUoTv37Pup7GhjWD1ep2f+ImfYGZm5vx7vV6PD37wgzz77LMbeWgXFxeX9ce2kYsLiHIZletz6vlXWfxFrYqcm3XCPX39IARy/CzmY49gnDmN8JjY23dg3XkX2u/HPHwQUa2imy2Mg/sRS4uoLdtf9xiXyoY9ARw8eJB7772XiYmJ8++dPXuWT33qUxw7dmyjDntN8Dd/81c8/PB3AXjrW+/iIx/5LZ5//ln+/M//lE6nwzvf+SP82q995BXb/P7vf5rbb7+DH/uxnwTg4MED/Nmf3U+vZxGLxfjd3/00fX39V/xcXFxccO7gCwVHrO3cHfxqtNvn6/9VfhBsGzEx4Sz63Q54vajRUeztOxDLS5iHD6L8AXSrjXH8KDoUxr7jLWhTItrtDTmtDXsC+NrXvsZ9991HNvtSh+vXv/51fuVXfoU9e9b/UeZSkLMzeJ55Cjk7s/qH18jzzz/L888/wxe/+Hd86Utf5eTJEzz44AN87nP/mc997r/wt3/7T5w4cYynn/4eACsry3z84x/j0UcfesV+7rvv9/jEJ/4jX/rSV/mRH3kPf/Znf7JuNrpcfTbid89lYxDlEnJ+zlnI+/pXX/zPCbzV647AWziCmJnGfPwxjONHQNmo4RF67/hBdDzh3PGXS+hmG2Pf84iVJey33IXaswc5PYn57DOIF2eLrDcb9gTw2c9+9t+89/GPfxyAL3/5yxt12DUjZ2cI3v95hNVDmx6a93xybUmcVUil0vzGb3wMz4u/JCMjo0xPTzE0NMzAQB6Ad7/7vTzyyHe58867+M53vsXb3343sVjs/D663S4f+tBH2L79BsCRk/7nf/7Hy7bN5dpgo373XNaZRgNZLjmaPS/+7b4u5yZ4KduZ4GUYiMUFjBPHod1GI2Agj71zF6JUxDx8EB3woxtNjCOHIRrHvvPtgEKeHUOUCgjLRkeiqOGN6fTfdEngV7c0Ly1JTHNtDzIv/5w5N420LdTIKHJqEu/cNNbI5WfZd+y44fy/p6ameOSR7/KzP/t+0un0+eNns1lWVpYxTckv/uK/A+Dw4YNIKTBNiWn6ee97fxw4pwz6V9x99w9d8DyllGQykcu2/XK42se/WK6qvSdWwBSwfQdMTBCorUDmpgt+3L22G8dr2trpwPIyRIMweuHv5RWUy0637tYBZ4jL0hIc2gvNJniB4VF485udz02fdT5jKjh5BBIJ+Mn3gG3D6dNQKoFlQV8aduxw7CmXLmzvZbDpHMCrtYCUUmvSzHm1to4aGMI0TMTkBMr00B0YQq2jVtDZs2N8/OO/zUc+8lsYhsHk5MT54zsyz+IV9mitUUqff880Ja1Whz/4g/vo9Ww+8IFfuuB5KqWuqv7KZtJ/gatvr4ykCVoaceKU8wQQSaMuYM/VtvVi2Uz2/htbLcup7JGGU8tvS1jtXF4u2BaOIaYWMY4ehmYLbAXJBNaumxGtBsYTz6G8PsTyEnLyLCRS2G/+fui2Mb73PLrRQPR6kE5h79iFqJYxvv0wulJBD4+QAFcLaL1Q+UGa93wSY3oKe2h4XR/BDx06wL33foJ//+9/hx/+4R9l//69rKy8JOdaKBTOj4e8EM1mk3vu+W2i0Rif//x/wTSv26/qDcdG/u65XAJKOQnec8PX1/K39jLBNjWQh1oV86knoN5A2Aodi2Ht3o1otzHGzqA8XnSphDl5Fp3JYr/tB53GryMH0a0OWB3I5LC2bkMuzjv7qlbRgQAi17e2ENQlcF2vKio/uO5/fIuLC3zqU/fwn/7T57j99jsA2LXrZqanJ5mZmaa/f4AHH/w2P/7j/8vr7uczn7mXfH6I//Affhcpr9uxDW9YNuJ3z+XiOT98PZVG+3yrb/BqnZ9WC/PZp6FaAxTCH6D3pjcjuj2no9fjhZUVPNMT6EwO6+0/BPUaxqH96E7HCfXk+rFGRzAmJ/A+/CC63UYHwjA0iPYH0baNKKxsyPlvuAN4+OGH/817X/nKVzb6sFeNv//7v6XT6fJf/+ufnn/vfe/7X/nUp+7j937v43S7He688y5+6IfedcF9nDp1gscff5TR0a388i9/AIB0Os399/+fG26/i8v1gKjXoF1Ge7znNXhW3aZYQHRe1PmxbYwXnkdUymgE0mNi7boZAGN2Bm0YiMVlPDMT6Fyfs/BXyhgH9qF7HejZ6KEh7PwgxqmTeB/6DrrXQ4XCkMmCz4dutRAzU8jCMqpvY54A3HkA1yjuPICNYzPZu5lshU1g77na/GCQ9A3Da7P1ZYJtmCbGkUPOKEdhIA2JdeNNYBhOSMgwnO7euVlUNofaeRMUixiTE9BtAxI1mEdnchgnjjnlpRpUOIKIRtEeA1muostFjHIBLIUKh9GjW4h/5NfcHICLi8u1iZydgRMryEj62gtv9XpOgtf0oPoH1jZ/9+WCbbk+5JFDGEtLKNMAaaC2bUcFgshSAa1BTE1gzM+h+gfove1uZGEZY/9e6LRBSNToFnQ4gnH8KOLgIYTHRIWj6EgYoRWyWIRSAVGrI0yBHQxDIIiOR7FHt23IZXEdgIuLy2VzrrcBUxC09LXT2/DikBWhFSqTXZM2/ysE2zJZ5MnjmPNzKGmiDYkaGkGFI8haBd3pIM6cwVicR/UP0nvbDzpJ3AMvOMldr8fp9DVNjGNHoVwEbwCScZTfj7Y1cmkBCiVoNRD+ACoShkgUnUyhM2l0ZOPKat8QDkBrjViLR38DsskieC5vUIzpKYTVg+07ECdOYUxPXV0H8OIkrXMxe70Wbf6XC7YlksixMcwDe7E9XhAC3d+PisWRrSa6UkacPomxtIjqz9N7290Ys7OY+/ciWk10KITetRt6XcxjRxH1Ojrgh6Rji261ELNzGKUi0rJQgQCk0qhoDJXJIhJJbENgzM5iTE3T27VrQy7TpncAhmHS63XxeteQwX8D0ut1MYxN/zW6bHLsoWG06YGJCbTpwR5af+nitXJ++Hoiubbh6y/fJplCzBTwHDqI7fWCYSKSKexUEno9RK2GOH4Uo7iMyg/Te8vbMGemMffthV4bgmGsm9+ErFaRB/c6Gj7BICrjlH2LahU5X4NKCSEFKhhGBYPYsQTkcuhIFFoN5OQ4Uki0Voips04T6I/98Lpfq02/coTDccrlZeLxDB6P97p5EtBa0+t1KZeXiUQSV9scl+ucc70NgdqK09h2Ne7+zw1fj0TWfvx224nzRyLQ7eJ5/BFsrx/t8SJCYexMBo1A1JvIY4eQK0XU8BDWHW/FmJnCPLgPel0IhbG3bUcWC44zsHvoQAiVDiM6XcTyCqJSRrTbaJ8XImF0NIZKZCCdRIfCiJVlzFMnUH4/FAoYlTJ6dAud/+PfQTC4IZds0zuAQCAEQKWygm1bF/yclBKlNk8V0FrsNQyTSCRx/hq4uFxNVH4QMjddsKt5wzjXlOVdozY/OPX3c3OIegdsG/Ppp5xmrWAQ4fGiMlm014doNjAOHUCWSqihEaw7tmHOTMPhg4hex1nEt2xDLi1g7t8LykKHoyCC0GxjLC+hW01kz3aaxiIxSCawczlkNIaSJsbiPHJqCi0NxMIshumle/v3wbYbkPPzyOkpGqEMvHn9L92mdwDgOIHVFsFrvjztVWw2e11crji2jVhxGqTWOnz95YJtxPyY+15AS4mKRECBSifRkRiiUcd45nvIYhE1Moo1uhVjbhp55DDYFsST2OlR5OIcxuEDaFuhQwFQLw5zL5XQXQthAKYHHQuhcjl0JguhCKLTQUxOYlhdRKfjHKevn+6Pvw9MD7KwBFMTlLw5ml4/icCFb24vhzeEA3BxcbmOePnw9XRmbdr8OLLOotlEaY1x/BhE/ahYDNHtQTSGSmcQlTLmk48hKiXU4CjW0DDG3Dxy8RhC2+h4AjueQC7OI48cAq3R4ZCzoM8vOHkEw0D5vAjDQCdi2AN5SMTRviCiVMSYOYqyFKJec6qTbryJ7rvejaxWnFLQgJ9qfJD6XJNYYYqkXUF5XDVQF5dNi5yduWjtn2u6rv4qISplR2d/rcPX4SXBNqWQp05i9HqoTMZR47QN7NGtjpbPYw87k76GhrAHBjGWFjCPz4NpomNhp1lraQW5MAeADkedpO6ZMUSvjfb60IEgwupCIo090I8IR9GGRC7MI1dOgVZQryESSXo/8FZIJpHlInJ5CZVK0/bGqJ4tEqodYiAh0fEYLWOEVsvH2vqVLw7XAbi4bDCXov9/zdbVXy0udvg6vJQb6HYdff1WC7u/H1mtOdLLt74ZNT6L+fB3EbUqajCPyuaQy8uYi4tgmOhkAu0PIApF5MKC0xsQCiOXlzDGxtA+DyIQAPwgJXpgENWXg4Af2l3E2TGMVgPRswCNGhzGvvNtgI1st6FRR2WyWCsVKs+NY6oefcMh9OAAvWaX+mKLjs+PObQxvQCuA3Bx2WDO1cjbQyMY05NrqpG/5urqrxYXO3wdXhJsazSQUxOIeh17aAhZriCbTexduxDVCnzrW3imF1H5AVQ6jSysIBaXnCRwPA4eL6JexpibR5smyudHTk/hKRax43FUMoXRrGN7fNhbhxHJFNo0EeUi5qmTTl5A2xAMYN+4E53NIVotRLuJjsdRpg+xuEjp0BwqFCK5NYYMeLHLVWrzbdoyhBWNY2EgpOsAXFw2Jedq5I3pyTXXyF9LdfVXhZd34+b6YI2KuKJYQJTLyOkpRLWCPTyKNA1kuYx90y5EvYH53QcRlTLs3IYaGXFyAysFhMdAhSNor8epxZ+fB48XZRgYp08gO12svn6svhxGpYJOJLB27oRQyHE6c7NO8lZIlOFBp5Oo0W3gNRG2gl4PFY0iSiXEqZNUGz7awQSJ3Uk8dKHZplzVtO0IPV+Ejj9ETwQoLtlEihvT8Ok6ABeXDeDVMf9X6/+vlhO4JurqrwaXMnwdHMG25SXk1CSiXMIeGUWYBsbKMtbumxHNFuZDD0KlhB4YREcjUK0iVioI00QHg2ivF9FqIpYWnIld3Q7GscPg96MGBlDVBqJeh5ERrF1vQvh9zlPG8ePIdgu8XgiF0amMozckJQLQPj+qZ2MszkKnS8WI0fBtIZZRxFQDunUqNYO6SGCF47RCSUSzRWGqga5XyBglIApsX/fL7ToAF5d15kIx/3OL+FpzAletrv4qcb5KJ5lC+/1r26jTQS7MY5wdg3IZe3QU4fdhLC9h7dqNtmzMRx5yFv5sDvrzyGoVrB7EQo7Ymt8PzQZieQHt8yEKBYypKexsDnvLNozCClRr9LbfiEjGQQjE4jJyaR4tDJTPD8kUKp9HRCIgDbSQ4PUiigXk3Cx4PdQDWSpmmIi3Sb8soztQVUFqIkYrn8fqaMzyMpXTE9jVFqlQh2WZ5v87+2bi2xO8ZQOuuesAXFzWmdVi/peSE3hDc05uORZHr3XylW0jF+YRp08gCwXs0W2ISARjaQlr1y5UbgDzyScQxaJTLZTJIWsNwEYr4dzhJ5MwPuvc8Xt8joxzuYK1dRvWjTsc5U/bprfndmdhb7WQY2edKqSAH6IxdCgMuT6nacw0nLBdp40sFEArdDRCu2+ElZqPSLvCgH8BDC8NX4qCp4+uL4I9v4J37AzlFYVV6RDOhymEhvmnM/0Ipbh5dJF0oATcsu6X3nUALi7rzGox/0vJCbwhOSe3HAis3QFqjVheRp446lTlDI9ib0tiLC5gbb0BlclhPvM0YmUJHY2h00lkvYEWAiwLHQigE0nE/DyMjUGlgXn6JNoGa9eNEI1hFArYW7Zg3fUOMA1nhu/UJNK20AE/Op2GeAw7nnSeHqSTmJbFIqLXQxsGur8PSxmsFAX+Yon+GMi+BI3IKEs6g7VcRk0t4BcTNIo2yx0/4aEk1dFBHnzOi6dZZU/wEL2VIo2xHpEbw7gOwMVlE7DazN/rZSbw7KxgeloyNKTI51+WxLSsF4enmGvX5scJEcmDB5ELc+ihYaxdN2POzkB0mNn4zTS//hw5PUc0H4JECtFqobs9dK+HDoVRff0Yk5OI5SV0rQ5zk5j+EL0bd2EsLWLOLWLduAN9Sw7RbCJmppy7fZ8PHQpCwI9OZZyOX5/fGfrebCBrdTAkKh5DBELYGiozdXSvR7o/iNy5jVYkTXG2Q/voAnZtgUDQot40OdtOELghR48QDzzUwrN8lu8PT9LtQbtoomNReqEI1URsQ74j1wG4uGwAq838faPPBJ6dFdx/vxfLcmas33NPl3y//ZI2/1qHrwM0Ghj79yGnJtCDg1hvugVzZgoRDmPtuJHad/ex7x/nackAypvmve9sEos3Eb0uxGLY/XnMsdOIuWnEchFzaQE7PwC33AKnxxHlEr1du5DBMLqwjHHqJKCdcZHJJDoQQieS6FAQIaWjCrrsVPsI00QPDKAMAxpNauMrdEWA8I4+PAMZepagdGKRxtwJuniI+rtU/UGOd/oJbkkhS2We/ad5fJUVbs+W6USDlEljBQxqkSSG3ySb1Xgzbhmoi4vLJmF6WmJZMDSkmZ4WLBwrMyTqax++DtDtYhzcj3H8GCo/iHXb7Rizs4h2E3vbDYjDB/BMT1Gb9dMwYiQTFtVCg/JCl9hQCjsSxTx+DHnqFHJlBVGvY23dgpXYgVxeAcvCuuMHEHYPY7mA6C2gvQbaa4IvgI7H0NG4U4lkWchKBXAEGnUiCV4v2lLIUoFW20MtkCHyfTtIhg1UqULp6dMUGx66lodYwkfPE2dvY5BkqEZi8jSHvnUMrRS7RxXNTJz5RpieGaRpRAhEPWSzMLYQ48xsmO8bsDfke3IdgIuLy7ozNKQwTVg+WyerygxsCTuNXGtBKeThg5iH9qOyffTueAvG0jyiUkWNbkEcO4L5xOPOdK5QmGgWfDRpLtvUIoMEb5bIs4cwalXkygoo6G3dhlGvYRSL2Fu30xvdCrqLmJ0GIZ3kbSDgDG2JJ1HBIEIaiI7Tras9XpTXB/GoE/qp15HNJi1vjHLiZgJRk2zQgvkpKidtVqo+WlaUWJ+PXqyfg6cVmeVjbFnZxwsTCTAM8jv8NDtexu0IVs9DJxQhkvISj5qcXYgyPi0YGVFkvBp7Y9Z/1wG4uLisP/lUi9/9hVmmihH6b8rRn19DI5PWzvjFZ59GJZP0vu8tGMUCRqmAyvUjzp7BfOpJkKD9IafrttUk6Re84xe2UZhpM1TcR+hwFaNUQgeDWFu3IwvLyHIJe+s2hNeAWh1jbhZSMYTHg5YGRCOoRAohhdPUVa8ifE6CV6cyaKEdXf9SGR0M0ekbotAM4QtKcpSQ5TrVacFiyUNDJIhuT+FTkrGHzpCoHmLQb3JwJodljtI33GFZpJns+ujioxeMEMt4MIMBzk4ayBoMDytAsLIi8PsFO3a4jWAuLi7XOr0ecn4O7fGSvbWfrBDAKouX1ojJSTxPPILyB+jdcSdGvYKxsoJKJpBTU3iefgIQqHAIpAfR7iCUhb19J3JpjuzxJ8hWy8huDzsSwtq2FWN5GdFsYm+/AaE1olZxwjlSojw+MAxULosOhkCBbLXQHhPt8bwY+jGhZzs9Ah4vKpXC8gQoVT3QkvTJRWRL0Sj3mGunKMkUqa0G2fETTPzDMbwBg0QmwuGFbdjLHcJZLyUjy7Q0qBJDBAPEc156lmBsShAIwOiootuVFAoQ9NsMRyvUixa1KQNY/zyA6wBcXFwun3P6O1Zo7dr8gJibw/Pod9EIrFtvQ1g2xvIiKpZAVqp4Hvge9GxUJIowDWfEYsDE3rUL8+RxPA99B12rItDY8QTa60VWKyitsW/Ygej1kJUaOuRD+/xO810ghM5kIOSBYh0aTWcIu2k6s4OtHgLAstHxmKPnLw1WOhHsSoNMoITZqtM0Isw00xS6YQbaZ9i+sJdjT6XRqQShGwc5ftCitujDnw7STeZY9oUo96LIoJdUGmo1wYmTglhMs3WrotmUFBcsEmaVIX+PUtXLghEit8VLpn9jhlm5DsDFxeXSebk2fyoN/SlYQ+eyWFnB/O53oNPG2nMLeLzIpSV0JIRWNp5v/0/odtHxBAQkst1BB2JYo1vwHjyA8Y3/gW51kF4TO51GW12MTgcrGkP19yOaLehZKNODQCLaPXQ6jQ6HodNFNmqQ7EfHDOcBxe45d/xao2NxNAJpmhAJU6p76BaaJPQU3pBB2xNnrJqldnSBvvZeMmHNvsp2vL43ERmsc/CESeFMlEjGC5EUtUCEYt2HzwvpPk2hIDl1CqJRzY4bFI2VNs2pOtGgjSfqZb4Zo65Nsls1xSIcOCApl10H4OLyulyK5v567GM9jgtg7n0e88B+rFtuxbr9jjVtc8Fa+yuAqFYQtZozfH2N2vyiXHKkGYpF1JvejEqkMOamwetDC4H54Heh3UQnEhAKOY1V8TS9/BDeQwfwPv8coMHjQ+VyiHod0WmjUhnw+5FWF5RCBfyIZgsBqL5+pG0j2i2nJyCVRNvaeWrpdtDBEDoeR6ORlo0OhSAcplbo0jheJuZrk+yP0SLD5KFlimNnSCcU8a39nBi/EaPWJOmr8tyZLCVrmHAuTCCbpG14WVmR+HuabA4WFwVjY4JIyGZ7rkGz1KEzq4klffSSKRYrBiGhSedgfl6wf78gGtXksopkuLsh36HrAFzeEFyK5v567GM9jgvO4h/56Ic4Vzhf+/P/tqoTeM1a+yvhBM5p81/M8PVaFfOJx5AzM9g7b0LtuQ1jegqjNYsWBp6HH0K3GuhUGkJBp2s3kcLyePAd2o9nbn/MnswAACAASURBVA7l9zudvMEg1KoYraaj7y9fvHP3+9E9A1mvQzKBGhpCNOrIZt1RFDVSiGYDUalCLAp9fShvA9FuIgwDO51CC0lneonaiSKBVJDciA+1XGP+sVPMFkOEdw0wcneUM/sb9A528IV67J3Ps9yNEx2JEY0H6HZhsSAIBjXptGZ5GWqFHmlvlf6IptE2aRg+wiMRWi1YqAsihiYWU8zOSmZmIBXrkA82icgO2wdt+m8Mb8hX6ToAlzcE66Gvczm6/Zer62Me2A+WhRrII+dmnSeBVRzAq2vtp6cl+fwG1QvCS9r8FzN8vVbDfP5pjNNj2Nu20vvRH3Ou08QYyvTgfeJRqNZQqRSEM05JZiYN3R7G3r14SgV0KIidzYHqIRtNbJ8fNTwKgPB6sP0BRKWKrJSw+/rQ8TiyXkVJgRoeQdYbyEoFFQ6jc30oBLLddEJMyQQqPIwsrmCPTVGtSoyAj0y4gVhcYP6Yh1NiJ5mRIXZsW+TU0RJzJ314UlH2FfuYrWbIjXrpD0OrJZif13i9kEppCtMN7Nk6yQjImJ+WiFHzm0QymmZTsLwsiMWcctnZKQgaDbIBC2HYJCKSgduDCG+McAQSg24IyMXlgqyHvs7l6PZfrq6PdcutYJqOcqRpOq9X4Vyt/fS0wDSd1xvCy4evr1Wbv9HAfOF5jJMnsPsH6L73xzFmpjBOHkf5vHifeAJRLaEyOXRfn9NUFY1Du4nxvScQ9SbE4lh9WYxGA1Eto9JprNwAQmsIR9CAKCwjexb2wACyWkV2LdRQP6SSyGIJWS6hMznsvn5ksw7NJiKdxs4PQthEHBuDqVlWaiay3SYty6BNpjzbONNKMhBb4Hb/FCdOenmKLJ50kH2LIyxMB+jrgy1hqNcFk5PgMRVZX5XKXJP6nCaWjaAyGTrCwOfTJIKacllQKIDfp9HNFotzHcIRyEcBv59UPkQ67WjV5XKaoSFFoQBLS7Bz5/p/tUJrfWUDh5dJoVBHqYs3OZOJsLyJZHU3k73Xiq1rjcW/nr1r2cer4+4bmQNY7dpuaA7gnDZ/r7vm4euZkEHpwccwjh5xyiZvuwO5tIgsFrA9XnxPPIool7DSWQgEwOdFx+JQKuE5fABt2+hkGgyJKJbB78dOpxDBACBRqSTVqRrNuQqBkQyxfBBZr6Az/c6w9lIB0euhcn3oaAzR7UDPgkgIlcyAx3SmftUqxE3N+Hgdq61IpgVieJDpVoLZg1VS3hqD8RqHl7MUdBo9kOXweJJCQdPXB5EIVCqCwpJNVNaIexsUCpK2J0IgHcLrcx6QfD6B1+ss/KLbwder0qpretpDKOXDCPvweCWZjCYeh3jcKQX1emFsTHL2rKBe1/T1KX7+5y/t70xKQSr12iGkNTmA2dlZxsfHMQyDLVu20NfXd9FGrBeuA7j22Ey2wuXZe6Xj7lfr2opKGdFooBJJZ6FejXYb8+A+ElNjVLSH3u3fhyiXMZYWUaaB98nHEYUiVjqDiITRPi86lkTMzeA9fBDbH8Duy2HUmsh6GTscR6VTCNOEQAA7lsCYn6Vetvifh4YxrB6GKbjr5/pJRXvIag0ViznlnUIjul0wfahkAkJhaNSRiwuIUgnabWoNAyM9ANkQvkyU5eMrTJzoEfW2Gdzh43ihj2kGIRbh5EmTQgEGBhShkKC01KGx2Cbo7RIMaOYacbrCRzzuXCopNYGAQKse9cUOnl4DraFl+9DBEOGEicdpPCYa1WQyMDioGRjQLC3BiROS+XmBUs6TQD6vSSQU73zn+juA1w0BPfroo3zhC19gamqKwcFBTNNkbm6OkZERfv3Xf5277777oo1xcdnMXPG4+5Wm0XDCJrHY2qQbul3MIweRR46A1wfv+kGseWcmrrJszOeeRhSWsdNZ9PAwIhRCRaLIU6fwPvY4djJBd3QLxtIy5vw8KpGmN7IN4TGdhLBlIWdnMJsN7HSWQtmihxdjKEd9oUZlqkbyrTlUJgMIBNrpEs70oaVEzM8hjx1FVipoj4dqIEc5tp3onjhDRoXxA3McfWgJXzbKjXclOVYf4nAzjYgoTp0yKB6FgX7FDfkW5YUO01WBLyIJ5MLML5tQg0wGfD6NYQgCoomuN6nOKISEji9Cw5shEJDEwhrDcG4cMhnN4KBmyxbnvZMnJc89J2m3nSeHTAaCQcXSkuT0ablh4b0LOoBPf/rT2LbNZz7zGfbs2fOKnx06dIi///u/59vf/jZ/+Id/eMGd1+t13v/+9/OXf/mXDA4O8tRTT/G5z32OTqfDe9/7Xj72sY+t35m4uFwBrljc/UpzTpvf719bGKvbxTh6GPPoUZQhsG69HQwJY2OI+RXMF55DLC2hcjns4a2IRAzt8yH3H8A3O43V309v21aMpUWM+UV0KuWUY4b8qEwaY34BY/wsdjbrzO1VCh2OEtpmYh7tUCiAFR0kfIsJYdDSQHv9EAlDoYA8fBhzYR7tNVH9eaoDN1Ay0kQaiwz2Fqg9f4anGwM0YyOM/m8pTlX6eHDOoNfTjI0J6mXNSLJIKqlZnjeY7/jxRuIYGcH8isBsQzYLQV+XiF2DTpdW12Sx4aPrSeCNGQjh3OVHIgohNKYJ+Tzs2KHI5TTz8/C97zldv7YtiERgcFBRqwmWliSGIchmHUeRy23M137BENDRo0fZvXv36278ep85ePAg9957L+Pj4zzwwAOk02ne85738JWvfIX+/n4+9KEP8Yu/+IsX/RThhoCuPTaTrXJ2hlRthcJlzNm9krX3L7+2a+kTuOh8hGU5C7+QzqCT1RK8vR7GsSPOtK1SEeuO70cPjyKmpqFaInH0AI2peexsCp1MQSaL7ik8zzyJrNXoDQ0jAbGyjAoGIRJDh8PoeALt82FMjiOkQW+gDyklGCYqFkf0LEcDKJFkuZdkvugjl1PkkhYqHod2Gzk9jTE5gZYSNTSIGhyiF0pQPFshUFsiIWtUO16O1EeRA/3c9qMj7D/UYnJS0mppJk91oVIlm1b0pI/5eohm10sgoFEKR5fHp8mn2sRFBb/HpmH7KPailBpevF4IBBSG4ej3BIMapZz/796t2LLF+V05elQyOSnpdDQ+H6RSYNuKSsWg24V0WpFMOk8VuRyEw5pYDPr6rmAIaLXFf7XPfO1rX+O+++7j4x//OOA8NYyMjDA0NATAT/7kT/LAAw+4YSSXK8a5mn1MQdDSl1yzn8/rKx72WUufwEX1JGjtaPPb1tq0+S0LeeIY5uHDUKvgeeoJhBB4Dh6g8467nXGKK8sw0Ie16yZ0th+KK3i/8wAoRW9kFCMUxFxZwY445Zg6nkBlsojSCnJ6wtHuH92K7HWQ/gAEQuhuB2k6g2NUMoWwLDL1GukhP1qasLCMZ98L0OuihkbovvVtiEQCVaxQObmEbE6RS0AtkefJ6h5UPsO2nQazs4JHHpPMnmqyPNHC7lqk+/20oinOFiSdjiAc1uiOZnmmS9DscFPOIhJo0/FGqYksZ1eckE00CsmkBgSxmMDj0fR6gmQSbr3VJpfTzMwIHn7YoFx2QjyJhBPXbzQki4sQiwnyeYXPx8sSwk54aHkZpqdhI1Kvq5aBfvOb3+QLX/gC1WoVAK01Qgiefvrp193us5/97CteLy0tkclkzr/OZrMsLi5etMEX8mRrIbNBQxU2is1k76aw9cQKmAJGRwlMTBCorUDmpqtt1apkMhE4cwy0gtERmJ4mceYYvOedr/zgufPbvgNe7/xKJWg0YOuAk2V8PWwbjh6FffscJ3HLLmeUorbBMGFmEv+D34Ldu+Htd0E+T+zgQfjW/+vse+cO53jlgvN6yxDkchCLOfuZOgODg5Dvg24XkgnQ2lklczkYGXG2W1wEuwPRIOg2nDzknEMuBz/2bme+b62Gnp6htHcBS3sZubmPZjLPgYV+tNfHbe+A+Rmbk89UKc61mZmr0zUCJHak6dmS2TnHv0aCFr5mlfK4Ipow2HmzxIymsIVJowkLC6CUE6d/UVuOZJLzks033wx79jin8fzzzqXr9RxHMTrqbFuvO9tt2eKEieJxZ4FPJJzXKyswM+PsM5N5qQR0vf/OVnUAf/Inf8K9997L8PDlzS1VSiFe1jhyzpFcLG4I6Npjs9gqI2mCliYwMUHL0jQjadQ1bve5a2tu30VESJiYdJ4Atu/CepXt585PnDjlPAG8+vxeNnydcAxqPee/18K2kadOYBw8AF4vqq8PgmHEchV57BS+U6eh1UL7/fTedjf2yFaM7z1O9JvfpBGOYW/ZhigUEWPj4POjszlUug86LYzT4+hAALu/D4FECx8o0+kr64G9fQcqnUYuLSGPn0V7DegpxMwURqmEikWxdu1GZ/sRrSbi9CSycohqy0MtPEDs5hE6sQwPTISwSoLRfJ3a9ByP/p3N3LKHmUoEyxMjPxykXmhx9Hgb0e6Q8NSoL/dY6EZIDfoYuc2P1wvVFpSneqys9PD7oa9PYVlOqMzvt7Esg05Hc8cdisFBzfi44KtfNSiXNcEgBINOc1ilIimXNZmMJp12FvpcTpNKacJhqFTg1ClJrweJhCKbhfl5J0E8P6/4gR+4wlVAAPl8nne9610XfdBX09fXx/Ly8vnXy8vLZLPZy96vi8taOTeLN1BbcRbH/OC61fCvB69ni3X7HdT+/L+9bg7ggrOG221ksYAOBlc/R6UQZ05i7j8IHqexTQQC0GxhPPkY8tABzEIR69bb0Plh7FgEc/8+PI8+SjU5RDFzG8HGDKHxcYhEsbdtR0WjyJUVjKlxVDyOvXUbUki034dAI7we1PAI1g03IGpVp+pnYQEtBZTLGIUCOuDD3nYj9juGoNtFTM9gTD+PRtAIZSkO7CK6JY7PH+LImKQz3WI4sUijaXPkSS9nCimWih48HohlNc1Kj6njFXr1JhFvl4VKgAWdIT8Itw6BUpp2WzA9ralUBNksbN+uaLUkPp8gnda0WpDNSm6/3QY0Bw5IHn/ccQyRiBPGMQzodASxmBPb9/udLuFMxonrt1pOZVm77cT6R0cVKytQKklqNYjHFX19ei0tGJfEqg7gfe97H3/0R3/EO97xDsyXxQnvuGNtYlXn2LNnD+Pj40xOTjI4OMi//uu/8tM//dMXb7GLy2Wg8oOQuQm1XFs3HZ/1YC22WLffsao8xCtmDV/M8HWtEadPYu7fB34famgQ4fFAr4t8/DE8e19A1mvYW7fR/uH3IGZm8D7+KKbVobf9RqqxYcYfn8Gnj1H2Ztj5w3uIhm1ksYio1dCpFHZ+yFHeNAyU14B0BuvmPSivF+PsWcznn0UjEJYF1QpSgxrM07vtdtDONZL79oLQqGiSxk23UZBZ/NkwATQTJ+rUVwr05WxaYT/7ZzKMjZusrAjCAZu+YIXCkmJmWmL4JDKYYKruRUoYGdHcMqTpdKBcFszMCISA/n6nTl8p8WJFjxPa2bFDs2OHYmxM8MADBo2GJhrVhEJgGE7y1++HWEy9WO8vzt/tWxbMzAgmJwWBgCafV5TLUCg4jV+RiFMl1Os5zWa9niMKtxGs6gCeffZZHn/8cZ588slXvP+Nb3zjog7k8/n4/Oc/z2/+5m/S6XS4++67ec973nNx1rq4rCPrpeNzzdmi1EvD1zPZ19fm1xoxdhpz/37wepyRi0JAt4d86EE8B/Yhu13Uzp20b/1+5P7n8f/NX0EwRHf3LoRtY46fRRcVVc8I3nwfrcUynYlF5EgAe3AQ4QugpESgIZHAvmkXqj+PnJnG2L8Xo9UGj+mEhdsdyKaxtm9HRKLI6SmMw4cchc9EHLX7zXTjaVasBBJFsF1h8dAixapJcihIaGuMQ+OC06cktaUGSU+RAVOwWPSz0ApgBj3UpaC0BNGohxtv7JHLaep1ydgYziCWIGzZouh0BB6PJJNR2DYEg4Jbb7XxejV79xq88IKJz6fx+xXptNO4ZZpOHN/vd8I6qRQMDGiE0MzNCY4elXg8mv5+TSCgWVmRjI8LQiFn0e92oVoVlEpOmXE06nQFr3WM8sWyqgM4duwYjz/+OL5LtODhhx8+/+8777yTf/mXf7mk/bi4rDfrpeNzLdkiigVEu7368HWtEWfHMF94Hvxe1Mioo7HT6WA+9B3ksSOARO2+mc4NO/A+9iD+v/oL7HSazl3vwFhawnP4MAT8WIPDiJwf45k29tIythnBu7sPO2EgtEYF/eiRbVg33oiu1/GcOIY88AJaG5BIIAyJDoXQyTR2Ko2xvIB59gxo0JEY9s7d6HQaO5ZgZcFGLlfxsUCx7mOsHSOaMQkmBSdO95g8WqVV7hAIaGKRMLO1DM2qdOQYmoLKvJNwve02zdAQTE0JDh8WNBpO9c2uXZpqVeLxCPr7odVypB927rQZH5c88ohBp+MszKkUCOFU/0SjTtjnRSVpBgedip6lJcHJkxLD0GSzGq/XWfSnppxFP5NxnjoqFUG16iz6kYizrZQQCkEopNeku3dJvy+rSUF88IMf5C/+4i8IhUIbY8FF4iaBrz02k62v7gO4FA3+S9nmtex4daz+1e+Ze58nceYYpe27sG6/43VzBKJWRVSrqLijo/96iHEn3ILXi8oPOd2zpTLeB7+FMXYG2+9H7X4zK2YK/3ceJGSVMG7chh4ZQY6NIRbnIRbDGhhEKgXdNvjDlAnTCOcJejtEEoKiZ4CJ8C6yowHypaPImWlEu4saHkKjnRIafwA7k8GoVRHz846sczSO7u9HpTPOXIBul8psk05bQ8hPRcUoVkz8XoVs1Zk52WbsrKDaMjATMbTphH06HQDN0pJ4MV4Pu3Zp/H7N4qJkYSGAZbUYHlaYpqDTEQwN2QQCkm5XsH27IhRSHD5ssrAA4bDC75cv9gVogkFBNKoIBgWhkJMcTiSgWBQUCgKtnfJQrZ0hMFo7fQKhkKDddkJNhuEs9NGo0xNgGM6/X6s461L/zi4rCZzL5fipn/op3vrWt+L1es+/f++99160IS4uV5NX9wG0P/BB/H/7ZYTVwzxyiGZf/6phl0vR7b+QHa+O9788fn/+OFoREZL6fb+P79vf+rc5glbLSfCGw6vaLibHMZ9/DkwTtWWbkxOYmcF85DsYMzOoSITerXcAPexHnqAxLTgSuxnlCfD2xWMEZ6dQiSRq23YnzNTrop1SFlQ8SdTqMpiPUUrmmPGO8v/86QK58jMsWUVCP5snvmUIkAhDoqIxpN1DzM5jFoqO9MS27ehU2mn8snqITofadJmqCqMjfXR8gvJ8G7NVIqhtJqZNzsxHKDSjBEMSArBUdEouez0n5GLbgoEBuPNO5059elqwsODcse/e7ZRbmqZkZEShFHi9ki1bFHNzgmeekdi2JJnUDA46sX/DcJ4gIhHnzjwW0/T1aep1QbHoSDzH45p4XFEqSZaXBT6f4wiaTacSqNl09IJSKWfR93qdRf+1Er1aO9Wu58pB15tVHcDw8PBll4C6uFwLnIuzs30H4sQpzAP7Lzrufim6/Rey4/WOe+44jI7AxCTeJ5945TbjZ0FKtMe7qja/mJrA89yzKI+JGh4BIRGnTuB56knk4gIqFsPa9SYorOB54mFUJMrs4O3MlJps756kW7JohmP4tow4A9N9PjBN9OAw2mMipIGKJVDbt0PEi3x6L2Lvc2Qr/aiRLRwv72TEahEPGk6cf3kJo1RBx6OorVvRyRQ6EgZLIZSNaLdpmFEKVholepidKqWpJZQWKNPP6WKKs5MmlYrG5xP4/I5csm1DraYpFgVaC7ZsgZ07FY2G4NgxQbN5rqbemb8bjTp37c2mJB4Hn8/mxAmTkycNYjEnTi+EkwCORjXRqCIScdoSBgYUliUol8WLiVtNOKypVp27f6/XCQk1m05Mv9kE0zwXBnJKQyOR1268bjYdiWmleFFOwvnsRrCqA/joRz+6MUd2cbnCnIuzMzGBNj1Yt9yKeeTQRcXdL0W3/0J2vN5xzx2H6WkwTbpvezu+b38LY2oCbdmoSGTV4etiZhrP009hez3Yg0OgBfLwIcznnkGWS6hgCHvLFsT8PJ7nnsLODdC94y2Y83P0T75AtyqoeiK0IwP0jQZQEQnRGGpgACEMQKEHhxz55qlJzIf/f/beO0qu6zr3/J1z760cu7o6R6CRCAaQYJRIM4qiLdnSk6zx2KaD1uhZ0iPXjGXLVvCyNc9PsmzZy37jkeRn2iPzSVwSJYsjj5VIiaQCQQowkUhkoNFAo7s6h+quXHXvOfPH6W4EdqMBkgCD6lurFlnxnrrV2Pvcvb/9fT+EpgZUugPr3iuYHqvizSiCokSLbwox7qHjDajeNehozHRbhVi4SXQqSbkmmR4sUJ3J4vfNMlsKUJRRhvJ+jhyRlEomL4ZCZuc8NWV2/NmsZnZW4DiCdetME3d6WrBjh1xi8zQ1gWUJurs1fr9eaK5q8nnNzp0Sn0/Q2mpcvMplCAQEiYQikTABvbERfD5NPg+ZjCQcNjV9c3xTvw8GNZ5ngn6pZBhBzc2mpBMOmx7Bubm6XDYBv1YzzwUCmoYGvepw9quBVXsA3/ve9/j7v/975ubmznp8tUngS4V6D+D1hzfSWs/tAbycOYBXowdwIZ9xVg/guuuxDuzHGjyBu/lKVM+alT88M4zvZ8+hLAvd3AKVMvLoUZzdz8PcPDrgh0AAMTGOrJRxOzohFEaOjSFmJhEaVKKBYjDFDEniEY9qUweTOk0qUiKRtpipRCgeHibl5Ah1xNGd3Xhd3aTCNrMj06hykdxgjqkZi9iaBI0bGo32vz9gKC22bRq/kSiUy6hsjqkxj1zBwp8KkdMR5nOSfN7o4n/jGxb5vMmJN9/sUauZnfLsLMzPm4SwYYMp1QwOCkZHzWOplEIISSplOPbVqoUQGs9TTE1FmJws0d6ul5qxIGloUDQ2mquLUEgTi5kh5UJB4PdrtDZBu1YTSGnKOKUSlMsmCTiOkXrw+5ev51erZhK4UjFUU5/PMIbOqLAvi0vRA1g1Adx5553LTgKvW7fuohfyaqCeAF5/eCOtFV779V7M/EE6HWXq+LAxX29InV+bP5PBt+NZw/tvaEAUi8jjx7Ff3A3ZeRM9LYmcnUZLG7ejC+m5iNERY/CuFKqtHbelDRkIgOPg9vSSnbd49rt5ZpwWKv4ov7RmP/a3v4snbHJOA6lP/i6pzjDMzpL0SmSrGp1sMHLOUpirrlgM7Q+gYzFwHNO0rlRQCiZzASYrcSIJI4Y2NSWYnTW74vFx2LlTLJVlJieNZo6UZieeSBh1zcZGQX+/eV88bm6OA2vXekSjgmpVorXH0JDFzIyhZ27ZEiSbLVIomLp9c7MiEjFaPuGwxrYXd/GmvFSrLQZtszuvVMyu3bI0Pp8gmTSTv+fW813XBPxy2TSGHccMfa2mwrHc38IbdhK4jjrqMLhgzn+xCENZEPL8VyejI/i2P4tGoCJRRLlMfts+3F0vEKtMY4eM9IqVzeEGAngd3Vi5eaz+I8hcHhXwobu7cDt6kJ4LjWm8VBqZz2EVCkzPp5mxQ3SGphjNVclMBEj4mlGpRgLjA6gf/hR5zxVGk3/9RtT4LNq2oaEBHY2ho1FEtYIoFBDZWZOgIlGm3AZGJyTxuCbgmOGoyUlwXcHEBAwMiAVxNNi3z0hwu67Z9be2msZuMCjo75ecOmWCfmen4eGvX+/heZJyWTI0pMlkwOez2LhRsXGjmbQtlyGdVmzaZJqtfr8JzLWaqfsXCibou65AKdMPMIlAYNsmgKfTZgAsGtVL9XylzFVJsSiWGseRiCaRuPiNq9bmz+A1awK/WpPAddRRh8GqPYBF83W/H9Z3o6fyy3/Q2Ai+nz2HrrmoUBBRKiMGTlD52W4mnp/E05I5UaW3tYLdnEA1NSFns4gD+xHlEjqeoHbzLZBIIBSohgaU42Dl5lFzWVRTE1QqNKoa5UCaHfYVyGaLWzbsIvgfUzAySdmJIO65HdUdXqLIqEQz2vGZgF+tIGZrRoYi3QRSMjMDQ8cksZimoUFx8qRkbMyUQ2Zm4NAhid9vAuupU2YoKhbzsG0zTXvttUZtc3jYJIRYzGjC9fQo0mmoVIzC5sCAoVv29cEdd3jMzUmqVYlte1x7rUdDg9HfCYfNxdFisK1UjD6/UkYG2txM0A8GoaXFsH8W+fmLTJ1CwTRuDX/f1P5fDn//zCYwmL5C5OVrYJ4Xl20SuI463iy4kL7B+V6j2jso3/87Sz2Apec9z0g3CIlqbmHmB7vJ/sN38LZupuG+609/wGgG3/bt6EIeFYkhVBUxnMHeuxMmJylMKSzXI+KHcdVAzq7SODWFnp9HKIHqaad2xTUIKcCy0ZGw2WLOZ6GpGbdnDcIznVbV0UkskeTd7dOUdv+MhoRHZHMHkz1/SHXfCeJ9zcRu6MVtbjGDZCELMZUDv8dwIcnQmN/4JkQ18/Nw4oRpnjY1KY4dk2QyEscxgf/AAfNcIqEYGIDpaVOC8TwT4DduNA3bEydMcEwkzOPr13uAJJ+X7NplbBUbGix+4Rc8ajVBNiuZn5esX6/w+03SWLyyEMLU743YscB1zRWAlItB3wx2GflmM0OgtXnP5ORppk4oZPoIq9kpLIdSyQT8RTXRy9kEvuSTwHXU8WbChdTvV3uNzAyfPX/Q3IIOBBFuDZVqBMdh5vGdRB/4IFJ7qC9ZzHzhH2m4pg3nuW2IuTlUIomQEnnqBPLQQazRMXStinRrhPx+ynYYqgXaSv0kJisQ9uNdfTV6TR+iWEIKjYrGITuLKNl4vT0QDCIAEQyimlqgVEQeOICslGlqbsZ997VYhTwCSK1tQb39RqPbozxEuYKORKCtGeXkjHfy/2W8k7WGd73LpatL0dqqOHjQBP5AwAT0F1+URKMmyB47Zhg1pRKA2fGvW6eYnRUcPgy2LWhshLVrNa2tikpF0t8vOHnSXEHceKPiF34BRkctpqclXV2aDRs8ymVwHLEwcKXJz3YYAQAAIABJREFU5Ux/YXZWUK1qqlWB45idv9+/2EA2gdhxTNM3lzM9hkVq5ssN0ueyfvx+Ux66VIJv58Oqy29sbMR13XoCqKMOLqx+v9prznr++DHsvbup3vN29Bldweq2PUjPo5zuJDg2QOBLD+G75zpUUxqERA4cR54cwMoMIQpFM1kbCKC1RaA4zxWBcVxPIDsb0NduphpPYuXy6GoFJS2k9tDRKHrTFQitEDUX1dSEsh1j/HL4EKqhEW/DJrRbM/z8gN9o+2D0+oVjoxMJ9DJU1EWFy1BIMz5ugl0mY4arIhHD19+xQxIOm6Gq/n5T269WTVDs7DR2mzMzgkOHzKTt2rWwaZPCccww1w9/KCmXBV1d8N73eszPC3I5i0JBcd11HtWqaYAumqzMzZmAv6B4QTBodt9asyTb3NSkSSZPN26npsQrDtLVqkke1aq5v1jeWo31czlQnwSuo46LwIVw+Fd7jdfZZUQxjxxEB8PGT/ccSojv1mvRX/GIj+xHA9VNfehaFav/GIwMYw8cR8xmF1TCNKJag/l5BIDPh7VhLbKry6hvKg9RLIDPhnAMb906vFjcsIECAVQ4hBw6hbNzByocwevswuvqgpoLfh+qqxMcB+3zG/vGVTaD5bLR0CkWoVo1ge/oUUFHB+RysHu3IBQyO+4jR0z5Z2Gwmu5uaG1VTE5Kjh4VJJNw442ari5NoSDYvdtM24bDcNttimRSMD4umZgw/P9AwMPzBFIKOjoU+fzpnXulYhq6tm00fBIJU25pazPCbIWCYeqMjxumTjSqF5y+Lg612mmaJ5xm/bwe99D1SeA66rgIqPYOKm//RXzbnqF6623LzhIsp8vvPP5d856bbsbbegOV9/4vWMePn9UDWPwc5Xq0TB2n/PY7cQcHoTlFPKbg0AGcF/ZCtYIOhUzULJdAClO0jsfxmtvQiSTCAoolRDiCSqVwe3sRza2IuTmEz0FbFlpK7IP7EfPzaNvG6+lFh6PgD6A6Oo1JeyRyemBrFVQqsH8/jI9LNmxQvPOdLvv2GSmGalXxzDOSYNDsuA8cOL3jD4UEfX1GdnlkBPr75ZJ8QyBgjMO++12JUrB5s+KeexTZrKRQsGhsVFxzjYdSpyUXFgezTpw4XUqS0iSYRMI0b1tajCnLqVOmeVupvHymzpk0T2Chb2B6Bq93rDgHMDMzQ0NDw3nfPD09TSqVuiQLW/mY9TmA1xveSGuFV7bec7WAVtToOQPO498l9l/+M9rzEFKS+/PP4tu986z34HmE/vTjiPFxNODecCPOtp8SsCSlXN4wbEZHEDXXGKWkm8GxTbBOxFHJBgiFEa5CBQKIdAq3dw26pQ1RqyI8D88fRJZLyJERqFXRzU1Q8/B/8+sIodGJBKWPfgJv/cbzS0ifg0rFqGqWSoI1a8Ls3FngxAmx5Gc7MmLKKLmcCfBTUybwx+OCnh5jjJLJmByzbh0LNX/Jiwuuj62tmptuMg3c6WnDHursVAvzZGKBgqkXkgIUixohJFob1k4kYlg5yaTG88SS42RPT4RiMXfRTB3XNesqlcxnLQb8841ovBq4rHMAn/zkJ7npppt4z3veQzweP+u5fD7Po48+yo4dO/inf/qni15QHXW8UXGuFtBLNHrO7Qkohe+HT6BdD1Ip9Ows/h8/DaGgec9AP85TP0QODiAyw+hYEjE9gTV4Auaypp4yl8e1HJxaDR0KIlwXLQS6rQMVCYNl9HiIxXDXbUS3tSFcD6E9457l1hCTE/iKJbxEEtXcDJZguhRl7uQUPZFmAleuwRoaRM7N4V1g8C+XYe9eyfHjku5uD63hiSfM5G4mA9WqGZAqFEzgn5w0gdfw9c3ufmzMBNQbb1TE47B/v+Tf/12STMKVVyq6uwXz84bamUxqOjs9fD5BJGLKOPPzRuCtWBRYlkAIo7YZDBqd/mj0tOKm328SwSJTJxJh4Qrh/PA8FhKLyRRSmoAfi106mebLhRUTwBe/+EW+9KUv8c53vpPe3l66u7tRSnHq1ClOnDjBb//2b/PFL37xcq61jjpec5yrBbSk0bNMvV/MziBKJap33EXg3x5Dz84iLEnljrvwPfsM1t5d6PkcnDqJ1qCVRsxMGrrm/gMwNkYZSdm1yBMl4nn4lQfhIKqtDR0MQTiCt/EKo9ApBKJSBrcK1RoiO4ssFVGhCCoaNwyjWBzV2cmQfx1/83dB4vND/Oqxn3G1fYpQ/ML0kBYpkMeOCb78ZWdhN2xx880ePh9861sW5bLhz8diHoWCuaBobDTG55ZlmrHpNNxxh2n0Pv+8icp9fZp3v9sod1YqYsl9KxYzO+xFQ/WhIUEuZ7R/wOj12LZprC4yd8zuXy9dzGQyxpCls1PR3r5yFUGp0wF/8WrhlfD6X89YMQFIKfnABz7A/fffz/bt2xkYGEAIwdve9raXNITrqOPnBct583pXXXNWvV/kc4aqGU+gkw2otnbmHQffj39E7eqr8TZuojaXxT7ej7dhMzqRRJUrWJEwHD+OUBrtc9BNLdQKLnNlH/l4J7MVj5ZkgWhnA94116A2bkZHouZ4pRJagJieRpaKRnahIYVOpdGJBKqjC9XdzSKNZWi7ZUTVNnTwTf4E31sG2Hhfx3knjotFw83PZg3/fWDADFwZgTRBfz8LksfgeXppmKmry/D1FzVwOjuNbs/hw4Knnzb1/jvu8IjHJbmcRmtD9ezqMh66tZqp009OGr6+bRsNnnjcyCpHo9DYqAmFzMDUckydTEbwN3/jW2o2f/Sj1aXJWqXMuguFswP+y+X1v5GwahM4EAhwxx13cMcdd1yG5dRRx+sP5wq3nevNK8dGsQ7sRy9Ej7PM15Uy8gepRtz1GxEzs1g/egpRyKNsG6Ym8D33jNEzVi44DkoLhFboUBDaW5jrr+LVJIVQE003pijfcCUiEoVKGZGdRczOIvJ5CIZRAQfL83Cbm/E2X2Xkn5cpTnd2qgWxUYEd6yR8XzNqhV2xEV1bZNIYSufhwybgV6vmuXLZ1PYXg7QQgmjUaO4LYXbUPT2Gb3/yJIyNCTZvVtx9t1hIFMaF64YbFA0NRo1zbs68bn7elHEW+feLHPymJojHL4xdMzQkcV1Tejp1SnD0qKSlxTh2SfnKBrneyLgMs2Z11PHGxWoGMPau54k+8HumE+o45L74T7jX3wjlMnIuawxb9u0j8I//NyKfR0sbb/MVWM9sQ8xlTbkm1WiSh+uBNM1YkUwhqhVijREaN65n2m7m6he/h79/CN2/h+pNtyIl6EQcr7sHEY0jZmfwf/vfEMEA1vAQxZtuWVE8rr1d89GPVhkaWrkkUigYq8LZWb1U/963z5iZx2Ia29akUopaTVMqwcmTiyJrZqRVShO0W1o0o6OCo0cFXV2ad75TA3KBGy/YsMGjr88cY3JScOiQSTpanzZjb2gwMs3t7Rcvi6A1NDYqajVDR7UsQzVNp42Q288z6gmgjjrOg/MawCiF88xPoVJBdXaZ57c/h2pphWIRMZ9FHu/H2f4cYmYGHQwixkexxkcRuXkTIV0PPT8H4Yh5n8+HsB285ia8q7cQ3LCW8Eye+J5d2PlJdDiCKJcRiTi1t7wVkc+hY3Hc3jU4/ccQsdgFG9y0t2va272XPJ7Pm8CfzZ4O/Hv2mGnbhgZNtarZvVsyPq4ZHzfSCMmkEWgLBs0VQCBgguv8vGBiQnDddYq2NuOaVSyaRHDnnR4gyGQEu3dLymUjySyEJBo1NfeODo+WFi564nZRXsF1TRJKpzWf/GSVkZHTCe8iiE5vWtQTQB11nAfLGsBobRq8lQq1m28h8OgjyKFTICVeewdMTmL1H0MceBGZzUIuh57LIsZGTUE7GjNXDFqjfT50Vw86HEb4Q6juTmrXXAO2D2anYc8e7PkiOhBAt7QiHAfteujubnRbG15755KtlFepvCJj+Xze+NkuctprNdi1SzI4aGrsSsHzzwsyGZieNs3clhaz61/U518UVSuXobcXbr9dU6tpymWJ58Htt7skk3DqlGTPHot8XlMum8GwRMLw/FtbXRKJ5d2yVsKivILrmvuBgKF9npk44nFNV9dLE97PM1ZNAFNTU/zrv/4r09PTZz1enwSu41Li5Ri1vBpYGti69TZq973jJU1fb/0G5EjGmK/bNqq1jaEPfJLyC/34r+olMZ/D//A/U5ytkSsKEvlhQtlRCASNHaLS4LnoljYI+NGxGKqpBa9vHdKtIQeO45RK6HgDWX8jQ523IW9uoU1MGnrn1BTu3W/Dvfkty56nxQE07ThYQ6cAVj1/udxp4bVKRVAo6KXA39xsmrLbthmTlbk5QSCg6ekx2vfFImSzLEzXmpmxLVvMhG4+b6aA163T9PW5jIxIjhyxmJpioe4uaGtTrF1ravqx2IUH/UrltJ4OGHmFePy10dN5I2PVBPCRj3yEWCzGxo0bEW82DlQdr0tcjGHKqwnn8e8Se+D30J4i8LWvMP+Fh5aSgHvFlciZaajW0JaNNTwEUjAyH+VrP17PxqEhOp59msjmKlXpZ+q5EwQqk8xJiS9cxKkU0LZjtHwiUViwdFTpZnD82NMT+L//HbMQf4CT/8dn+PJPNxA7miP3kxK/9sFWeg//qxGQ+/pXKS7s7lc6Txdy/ubnDZ2zWjW19qkpY404PAxNTYaK+fTTJvCXSoa/v2GD8cItl03w9/tN2WftWk1np1HoLBSMxv+117rk84LhYYuDBw1lM5HQXHWVoqeHBZ/dCwv6y7loRaOvDz2dNzIu6ArgK1/5yuVYSx11ABdhmPIqw7ftGbSnjJHJzAy+bc9Qu/te5NTkojMIIps1kgvBEPbB/chv7+Ge/ceJhFyKHqhjwzj5GcIVTYgilCp4jg87noB4AtWQQETjhh7amEaHAghPIUeqEAyhmpuRY2PI53YxGrgPeV0nRw5XmT36E9acc06AZc/Taudvbg7Gx09THoeHBTt3SiYnTbPVdQVPPGFMWZQyJZ7OTtMPOM2aMfTO7m5NMikWtHOgowNKJcXMjOSnP7UIBASNjZobb1R0dhp9nQsJ+mfq6Wht7BIjkZenzVPHylg1AaTTabLZLIlE4nKsp446Lkhw7VKgeuttBL72FfTMDEIKaldcgTxy2OjcB/yoQBAxP4f1/H9gHz6AdeQILbNlZnSQ4PgQLdVpwnEXlCJQ9XABZQcJpJJ47Y2IWBza2qmuXYtVKKBtC2/tBlRPDyR+jLPjOeTUFEQi1H791yg90czJkwsaNls60ftfek6WO08rnb/ZWRgdNcHatjVHj0q2b5fkcoYGmc8LnnvOGLBIqWluNs3cbNa8b1FDv7fXXCGkUhAIKLQ2xutSGiaQ5wnicc3NNxvJhgsJ+m9kPZ03MlbUAvr0pz8NwPHjxxkeHuYtb3kLzhkFtteqB1DXAnr94VKs9VL1AM41hT8Xzve/g+/x7+P29OLefic6HIZKBVkqIjLDyP0vYj//PKJSxIvEscZH8WZzVLSDz/KwVdXo7AiLQrQZ0dVBoKsZd8NGCASR+RyqqYnalVdhlStGeiGVQm28AufHT53Vf9i1S9LfH6avr8DWrWrZc7LSeVp83O3o4sXZTl580aK3V9HRodizR7Jjh0W5rLEszeixEoeOSuYrAQIh09h1XaPSKYSRTOjshPZ2cwuFNNPTpmHc0QGRiDFeDwYVV18dJhrNk0gYBlAmszzN1PNMwC+VTMC3LDN8FQq9aj/1qngj/RuDy6wFtLjj37p1K1u3br3og9ZRxyvBoqrmq4nF3gK2IOTqs2vjnoc8dgwch/Jv/BbYNkIIdC6PHBrEeeanWMf70V4Vcfw4oljAUgqam7H8NqFCDkpVhN+P17sWtfEKAo2NeE0tKOUhggFUdw9uKo2YmcEuFPDWrsNLNYIQyMzwkqic/4nvk0lv4ZFH1mDbsH27Q0tLlfZlzslK58lr62DE6mDfPotHH7XRGmZmLLq6PIRQeJ5ieFhw4qiiNKOIiBy94QxzdieZjA/HMdINRpffmLX4fCw4bMG2bTaeB4cPw/vfr9m61SOZ1KxZA9PTmkxG8Ld/e3ry9iMfqZJI6CUBtVdqm1jHq4MVE8CDDz4IwFe/+lV+4zd+46znHnrooUu7qjrquARYrI3Ttx5x+KipjccTWMePIYeHUD4feB72ieO4iSQiM4x/53ZkZgSkhVIucnAQkZtHOwFEKQcTE2ghIRLB3XojatMmhLRQqUaEz4F4glpXN9KrgadQ8Tj6is0v0Ss4t26f3TuE666hr88E2aEhuSxn/1xobZQ1JycliYQmnzf3XVczPS2o1QQgGRwU1GqamFWjwT/BlEoxlA8QdDw6OjXd3YLubkUsJgCBz2car0IIBgYEgYDR6JmdFbS2milfOF3mGRyUlEqmVDQyYjR4brvNo6mpHvBfT1gxAXzta1+jXC7z8MMPU6lUlh6v1Wo8+uij/N7v/d5lWWAddbxaWKyNc+IEulyBsTGcUych3oDq6oFcDv9//2vE+DjOXBZSKVQkinZryKlRRLUCloWuVhHVKhoglaJ22x2IhhTacaAxjYpFDU002YCo1RDRKF57OzoaW1FX/9y6fWJLJ/Z+lnoAnZ3qvN9NaxgcNPX7REKRSGieeUbys58ZFs+ihEOpZBq70agmHIbcbIC5UpqIU6E7Ps+m29O0bWBJeM3v13ie0d9JpzXr13vccgv8j/9hUSgIQiHo7lZLxuieZ+QVQiHDwZ+ZMZIQV1+tLpmxeR0vHysmANu2OXr0KOVymaNHjy49blkWH//4xy/L4uqo45XgJUYt8QTVG28i+MIu3A0b0MkE2nKwhk+Z245nEQf2G1G1ahXt9yGzc6YgLjQUilCrQiSKbkzjXnsdurUVq1xGATQ1o3rXQNiYqKjmFnQqBY6DzAxjHzywYk9DtXcwcP8nyO4dIrGlk6at7dxPjf5+H319tRXVK3fskDz9tEVPj+aWW1wCAdi2TfL885KxMUOfDAaNIqeUpqkrhGECaW2ones2OVzfNY2MJxmvRtFaI6Wp63ueJhwWXHedx8aNp3fvb3+7y49/bHHttR5Sng76TU1mAri1VfOJT5wtNTGxK3PW96vjtceKCeB973sf73vf+3jyySe55557Luea6qjjFWNplqBaQZfKVO+6G6ouvu9/B/w29tgkqrkF3yNfQcxOo/N5dEOD2d27LtqtgbOgMlYugeeZhNFzJaqjE5IpdDKJDgSQP3gCy7bQx/vxtlyLt+kKdCT60rWch5efyQj+5pE1uO4a7P1wPzUeecQ5pwdwOglUq/Dtb1v8t//mX6qpv+99moEBY2noeYsDXoZS6fcLPE8zP396gveqq0wZJxAIUCx28+ST9pJX7oc+VGPtWo+HHvLhefD88xYPPlglFoM9eyR/9md+lIKf/MRm3boyW7eaK5QzmT5nSk1M7Mow8eDnEG6NCduBz/9xPQm8DrAqDfTmm2/mc5/7HNu2bcOyLO666y4++MEPviI56IceeojHHnsMn8/HL/3SL/HhD3/4ZX9WHXUsB2v/i4jBk+hozMgluy4IEJaxThQHDuL78sOIoUG05UOUS5DLQzgEpbKpu5RKJlqm0+ieHlQqjW5uQqUazQBXLIo1No6Ix/C6e7GmJhCed1bwhwubazhTrXJoSLB3r5Fr7uuD3bvh8cdt7rvPJZnUnDghGR4WPPWUtWCpaOrsX/+6pLXV1P0nJkwSEEJgWZpKxUzJrl0LW7Zo1q41XzEc1mSzZgcfCMANN2hmZgTptBFmKxROi7mdPCm54w6P4WET5Ts6zHH37rWWEsByyGQE+781TGvFnAPfiOlx1BPAa49VE8Cf/umfIqXkE5/4BFprvvGNb/DpT3+aP//zP39ZB3zuuef49re/zWOPPUYwGOSBBx7gBz/4Affee+/L+rw66lhCqYT1wh4zJJXPGeG0BYllbdkwO40eGYHjx8yOP5mGmosoFo1do7+MVhqtlLFbbG/B6+kxU7uNzajmNCLVhGpuRrW2QSSCWpvDeWEP1swU2ud/WSbxcI48s20UNffvlxw4AAcPmkD/k5/4+PVfr+E48PDDNsePC6amTvcUPE9z9Kip3QcCJvDXaibwb9wouOkmTVOTqf3n85pSSRIMCn7xF11sW/C5z1mMjooF9ywjrRAKGX/daBQ2bVJLa7Nth5GR02tdCYs6/M5YL/fN+WhkEO03PY46XnusmgAOHjzIE088sXT/5ptv5h3veMfLPuDBgwe59dZbiSx0hG677TaefPLJegJ4E+PcWvy5+voX8p4VUa0ijh/DPnQYqhW89g7cm9+CLpcR01M4+/fhdfeiR0ewTgxAwA/aQ9s+RLUEgQAaAT7jyEUgCC3NuL1rELYPlYyjO7tRPb3oljZ0QwMjs0FOTjTSGRK0r9Nn1+7bX7qrXc4k/lwsJ8/c3FzliSd8jI0ZiuXEhGBgQLJ/P+zcCaAWqJkKUMzPQzhsEoDrGqOUG28UrFnjIaUkEFhU2xRcc43HlVcqSiXByM5xcgdH+F9va8VraWfNGkUnQ1hDp/iT3+5moNZ1Fpd/61bFpz5VYds2i1tv9c67+1+8smm9toPH+SR3rz3Btf+po777f51g1QTQ1NR0lkF8sVgkmUy+7ANu3ryZv/iLv+CDH/wgwWCQp59+mhVm0ep4E+Dc+nfl7b9I5L/+6Yr6+su95yU182oVMTqKffQQTEyiEwm8a65Bh8LIwwdx/t9/hdw89o6fIebnsfbsxkrEwedHVyqnqTBGBAcCPvAF0J1duH19iEoNEYvgbroatWE9OtUIkTAqGiOTjfA3XzrNb7///hqPnFG7/+g5tfpFXMhcw2LNXCkj1TA5KVi71ujtHD1qpBmGhsySjfrm4u5fABaVCvh8Hq2tcO21gmuv9ZidlXz/+w5CGKmGP/zDKr29JgmUSlA6mkH8n39FfOFcN33+j2lBL53/NbZDyznnP5MRPPGEjevCE0/YXHXVyhaLZ13ZtHRy5YebaTqPHWMdlxerJoCWlhbe+973ct9992FZFk899RSNjY1Lk8IXOxF8yy238J73vIff+q3fIpFIcMstt/DCCy9c8PtXmmi7EKTT0dVf9DrCG2m9K6718BTYAvrWw8mTBHdtB62gpxuGhkj2H4T77nrpeypFM4KazxPMTUF8relo9vfDwICJvps2wR1vNUT5Xc+ZrmdPD9y0FR5+GAaOm8+rVoEF4ZuFW6haNR3LaBTWrzfF9tlZiEfgjjtg3TpDkYlEjNj9gnj84WfMofv6DEWzv9931v1czrdkNXixqNVgfNzw9mdmzO3IEbBtH8WF01GrGectpUyzVmvzmOOYRHHXXfDWt5rnFxuyTU2wYYNx6wqHfVx99eljvvDkBLZWeD1rsYdO4vZPkgpz9m+Wm4L0ptM/z2HO+53P/FtIp+GznzWv6+mBzs7Xl3rbG+nfGLz66101AXR3d9Pd3b10/5WUfwDy+Tz33nsv73//+wH453/+Zzo7L7weWJeCeP3hfGuV0UZCrkYcPmquALbeTOQ734WTg+YKoO8K3HPeaxdqRPfsZVHrt7D3APpwvzFeaWlHb7gaMTqC9ZNnwXXxWtuhrRt58BDO1x+DzBCyWIByxbCAlIJQGISEapWgz6GUTJlJ3M4OnOwsngvVX34PorUdHY8bHn8wCBqYKS6tLRoVuK7vjCBYY/t2Z+l+NFplcvLi/j5LJWN92N8vGB83nrtTUzA6qpmcDDE8XKZQ0EsuWT6fWZplCfx+TbFokUwaimdTk4dSivZ2w/Tx+zVDQz4mJ836GhvPXp/dl8YVEnHyOK7tYPelmY7qs36zYrQRdcZvdO45OPM7L/e3EAjAxo3m/ycnL+rUXFK8kf6NwWWWgljEgw8+SLlcZnBwkHXr1lGpVAiuYDN3IRgeHuZjH/sYjz32GKVSiW9+85t85jOfedmfV8frG8vVv3PpppV7ANUqcmLCSCQUi6Bc5MgwtdvuRFTKyEwGRoZQXT24116PdeQQ1vbnsA4fQo6PIyplNMpMJUmBlhLCUajVED4fpVCcoh3B7VlLuDkJoRDlu9+Gbu9Ep1Im8J9Htay9XfOBt59geNsQHbd2snlrGy0t57dWXA5jz2fof2aCYrqTrL+JiQmJ1pq5OcHIiLnQGRqSlMumnm8asybgBgLmlk4bXX7b9rAsuO46zebNZg3J5OmZs0/cP7Bij6JpazvzD/4280/vIXbXtTRtbUcB5ft/Z+k3Ord0dSF2knW8MbBqAnjhhRd44IEHsG2bRx99lHe96138wz/8A9ddd93LOuDGjRu59957+ZVf+RU8z+N3f/d361pDb3KcW/8+11SdahUxN4fIziCHM4iDB0zD1vPMdtd2kJkhVCKJuvIqEAL72WeQL+7BzoygSwVEqQKDJxC1qol8iSTO3Jz5/FKJ6pVXU4g2Yv3sP7CFwDpxnJn//Y8J/+a7UE0LspcXgIldGfR//Rwdbg39lMNE+o9p39q+qkzD4qRsPi8YeG6CgU99k0w1zbweIfae26iEkgwOCg4cMNOztZqmUDCJSEqLpiaPUOi0PPO6ddDaqrnhBiOzHAoZvZ5zp21lZpg1jyz0U/Y7FFvOrudP7MpQ/vyX8bs1ygf3MbGxhZYWTeCR/2m8B/a/SLGlddkkcCHSFHW8vrFqAvirv/orHn74YT760Y/S0tLC5z73OT7zmc/w2GOPveyDPvDAAzzwwAMv+/11vAmwFPRnEWNjiJkZZCmPiiUQAT9uby/a5zc7+lgMb8t16NlpnH//N5y9u9DlKrhVtOvBbBYxP2OGtUIRRCEPU8bBToXDyJqLCASoFRU+aVNItiLzeeYKfoJd3ass9Gxk9w4h3BrVtvPz2bU+7UvreeZ+pQKjo/Dk/1ehUFyLjEeYmfDY8aTFtGdkmcFILyzW8CMR8z6ljDjbli1w1VWKa67RS1cCiYReMX+tNoOw3Pdp36zS35gBAAAgAElEQVReEz+GOi4/Vk0A5XKZvr6+pfu33347f/d3f3dJF1XHmxSLQX8ui8hkkNkZqNXQgYBRy0w1oMtVVCqFjkYRjg/dmEb5A/j/9i+Rg4MgbbRXQ5ZKMJ9D1Kpo20Y1pJG5E4j5LDoQQl9xBRw+iNQgPRcxMUY4kqLkjyMqFTxfAN9tF3/lmdjSyYTt4Bs5rdmziEVf2lrNXIQEg8a1amJCLJiqS0CjEgmOVrsZGU4ySwJvMo5wjLm6ZRnvWhBMLvjQJJPw3vdq3vIWzYYNCtsWWBakUmd73i6H1WYQlvs+Xot+TfwY6rj8WDUB2LbN3Nzckh3kwMDAJV9UHa8PvCqa/ItBf3baKGmWS6A02u9HBwIQT6Aty0TPmotuaaJ60424vWvxf+NryJMD+L/2CJ4/iHYVspRFlIwmj7Z8EIvCTBapy+jeXghFUGvXUrn9TvSeXcg9uxFa4d1+F77xUbL/6SaqeQtv62Ya7rv+or9j09Z2+Pwfk907RHhzJ3ZvO+N7xrBGMji9rcTXteE4cOiQ4Hvfs8lmjYKmUopTpzTHjwtmZxPMhq6hpDSesAgFDYsnEDCNXdc1w1rhsKJahVtvVXz4wx6uK/D7BcmkXrFN8RL9o1VmEM78PosaPQpWnVuo482BVRPAhz/8Ye6//36mpqb4gz/4A5599tmXPQVcxxsHr8iXt1ZDzM4ixseQw6cQWqOlhYpEsTwXLBvl2EhPIzwPFU/gbbkOlML66U8IfOV/Yg30IzwXz3KQuoY1PQ7lGrpSBr8PnACiXIBKGbV+A7qlCfvQIbAE8thRnNZ2Kr/6a/D+/0zwH7+ANT6Kth3Cv3ovXVs2MTmZu+jvuOhapdo7iLV3GBvEuWGSX/4swq1Rk34O/dqf8OxAG1/+sk25bBio6bSH60rm581UbbGoEcIiGDcsncXmbrUqiMfhyis11Sp8/et+hIDBQR8bN2o+8AH3Zf1mq80gNG1tf0kZ61L4MdTx+sOqCeDOO+9kzZo1PPvssyileOCBB1i7du3lWFsdryEu2pe3VkNMTCBHM4ixUYTjoIMhVEsbcnIMUXORhQIqHEHUqshIFK9vHdofQB7Yh/9f/gnrwH6YmUIg8LTAcl2smRkTeas1sCTYFqJcBr8fddXVeGvWoRvTyIkJ0AfRsTjCkqibbsK7/kYQYsXd7Grf0fNM43Zx4MqyjE1hInGa9WIdGWQm7+NY6HqOHdOc+maV3TOCsTEQwvjoFotGjrlSMaYqySRIaSQafD5wXRP4r71Ws2WLZssWxac+5VvqARQK8LOfWasmgNW+TyYj6sydOs7CiglgZGRk6f8dx+GOO+4467m2trZLurA6XltckC9vpWJomYfnsftPIYIBvIYUqrsXa3IKMTuLlcuhYjHwK3QohOrsRsdiiMww9g+eQO57AevUSUSpjPI5CMuC+RxWvgDKg2oFtEJbDsLTiEgEd00fuq8PlUiZHX+ljGpvQ7d3InwOOhLFvfLqJR7kiq5Z53zHWnsX+TwUCuKMJuxLXau0XlTaFAyc2sDgqbcykkuQKSUZ163kPdMALpVMP0BrQSBgpBmEMDv+UMhUvRIJo9GzZYtm0yaF328Szb33evzoRzbFovkad921OuPmfL/ZoibP4gTzRz+6/MRyHT9fWDEBvOMd7zCWeFpTLpcJh8NYlsX8/DypVIpt27ZdznXWcZmxYu14fh6ZySCz04Z9Uy6DX6AjSUBgnTwJUqKSDWYM1HFQ6WZ0QwMim8U6tB9x8AD2kYOI6RmwHLQvgKhUsMbHoFJDKxdRrYKr0EE/wvFDOITX24vXswYdiSIFCMtCdbRTa+uEZAL33l/EGh+74Lq119bB9K++n/Lzh6hsugrt6ySkNen08jX2fB5yOWOw0t8vOHxYMjbWzEjl7ZTmStipMJ4OMD4OwvOI2FUCMYkVtBeUNw2Xv1w2SaGtDe66S3P33R4+H4SyI6RHTkBPF7/zOx3YNjz7bJC3vrXCb/7m+Xf/5/3NeKna6IU6jNXx5saKCWDPnj0A/Nmf/Rk33XTT0gTwU089xZNPPnl5VlfHawrV3oFqajY7+X0vInLzRhs/FEE3NqGVJvTl/wfcGkFXU3n//4bX1o5wLIhEcZMphOtijWQQO55DHjqINTKMLpUhGkE1NGBlxpDZQXA9tFtFVGsIDToWB78PHY6gOtvxunrQoTAW4EWieO2tqMZmiMZQjY0sEuBV75rzfqdSyUglj40J5Ngoqa9+nVbmkYPbKV75cVS04yWvn58XzM7C8LBg927ByZNGJrlQgNJMFY6cJOslyE4GEDFNNOrhlGexLA9/qUa4tQkZ8FGpCBoaoLsbvvMdh4kJM/C1fn2ZO/tOEf7S2fX73/zNDn7/92FycvXgf9ZvtkzyO1dtdDWHsTp+PrBqD2D//v1nNX3vvvtuPv/5z1/SRdXxGkJrRCGPyGbNrZBH11wIBtDpJpRtIU8NYe3bhxw4hlAa+voQp4aQhSJuWxvadhDZGew9OxGHDmGfGEDk8oBCpxpRkSr20aOIuTk0Gl2pIZQZ+tLpZgj40aEguqkV1dlhuPy2g04kqLW0QyKGDsfQsSi6IXXer1OpGGpmdWE+zO/XtLaaUowzOEBAzuF1diPOqJlXqzA7a4L++Ljg0CE4eFBSKJip3FptwSRMwPiIS8HrxEYQpkRNh7BdTSwwRzAiUCWXarlGa6+f227TrF+v+eEPLaQ0evrj44Ljxy3eFjhdvy8fPsXhx4cJ39f5snWFzkV9ereO5bBqAlBKsWPHDm666SYAfvrTny5RQut4k6BcRuRyhl6ZyyGqZah5YFuGk29ZWENDyJMnkX4/XnMLbiqJFQ5iv/gCTE+jI1G83l7kcAbZfxjr6GHIziOrFZTjR3W2w+wczt7dWLl5c1ylTJkxFETFG8BvQyCEbmw0HrqJJMK20ckkbmsbBEOoSASCIXQ6vaxkQ61mSjWVivkbdRxTx/f7T79m0Y/9zJp5TfoZCfQycUCQyZjb4KDR5NH6dENYKbF0FeG6YIVsfDMVQGKhSLdr/CHQAwpdge7oLG/99RRd14UIBo18w9ve5vH97zuMj5/W0/dazFrKh0/x4iE/3xRrmNvv47OfveAh5VVRn96t41wIvYoW886dO/n93/99HMdBa43Wmi984QtceeWVl2uNZ6EuBvcqwHWNS9YiT7FUMhIK1RpaCIjFzJXA2BjW+Ag6GMbr6IRICGbnEFovlWfk9BTJ0UGyswXk2CgiOw2VCtp1kU4AL5mA0VF8B/ZBqQDCNurFjoMKhSAcRgfDEPSj4wlUazuqqRkrO4twXdzNm9Gt7eho1DB/Uo2GOrMAzzMBv1QSjI4aMbV16xRr1+plWS8TuzK4/ZPYfWlSW9o59OQoh7fNUmlsQzWkyOcFk5NGk6daNc3eQsEkjelpyGbNzj8cNo9rDX5ZxV8tYEdDODE/lgVrmnK8tXuY1iuS2O1NtLWZZvLi4Nbjj1tLevr33WeCsswMc/jxYb763BpCGzo4fFjw9rf7uPXW/Btmx/5G+nf2RlorvEZicNdffz0/+tGPlozhN2zYgL3a+GEdry8sCNHIQh6UMv4LlarZ8VcqaCnQ0RjYVeTkBOLkCQgG8XrWUOtdg5yeRBbyRks/EjaWh5UyYugU1r4XoZDFmsshqjWjWRCKoEMhOHYM3/ZnzTEcGxGMoB0bHQwgwlG03w+OhY41oDo7UW1tCNdFInCe34FQHvaBfZQ/9ADe2j5jtK6guAxTx3Xhy192ztHpd85ivThjw4z+l79motbAcK2J0q/+Ot/f2YPn9aAUXHmlR7Vq+Pj5vEku1aqRZS6VDHMnFluc9tU0NJjEoLUPIXwEg9DRAbfc4pFOh/H7N9DdrWhuVmexiFbS01ftHYTv62Ruv49Th+HQIYtAALZv99VZO3VcElxQJHcch82bN1/qtdTxaqJSMWWdBXE07fjQrofMziJqLtqWJpB7HnJqAjE0hPY5qN61qM1XIacnscbH0VPj6FAE1dYOUiDGxrF37ECOZ6BcAU+BpRGuh5dIguviHDqINXDcPBcJQmMa7VhmACyeMKUe4aAakuj2LlRrs7mqSCTx2jtx9uwyXP6uLsTEBMW8ZjYXRs2ZgB8Ov5Sps2/f8p66HR3GKP3HP5a4ewrk569CNjchxscYPzhFsdgEaLJZwZEjhqNfqRimzuysOV4waNg7xQVV6FRKY1mGrrn4fGen4IYbXBIJQTBorkJSK7QnzsfIWazVP/64jRCwebPF4cPUWTt1XBLUt/JvFnieKeuUSqA12udHB4NQrSCnpxCeQjsOOhyCmVnk6ARCDaMCAVR3DzrZAKUycmgQe/AEOhTGa21Dh0KImRms/9iBzAyhazUj2VAuI2plI8eQSqDzYzjPP48YOgkCvGQKEV1IMJaFTjWh0eCz0U3NqJY2vKZmBALiUdyOLkStBuEw+Z4r8NhJbVSj7XZkSyeNjXrRk2VZLLJcTp0yW+2ODsVPfiLJZARCQCCg0X2NlGyBOzZFf62T6Won09OGkgmmsjQxcXq3H4+bnsKZ3jGOY57zPPPfjg7Bdde5RCKCSESwaZMikTj/T7UaI6e9XXPffS779/s4eZI6a6eOS4ZVewCvN9R7AAvQGopFRD6PUB5aWuhIBBwHMTmJnJkCzUKZxUFMTCDmZhGWjYrF0M0t6GAI7fNhDQ0ih4fRto1ubzd19kIB++BBxPApRKmE1gpRKiOrNRSmpKG1RuTmiB49QikzivD7cVtawOdDlsuIQBC3sQFZqaEjJqGIhhRucysSgY5FUV1d1HJVihMFahMz6JZWrC2biRXG8I+urkWzcBooFAS7dkmef17S2qpZs0aTy2kmJyWxmCYS0UxOCvb/cJzJ40WqsQSB9gZUPs9Mpsa8ihCIO4RCZmdfLJpqllICn8+YowcC5uogGoWuLti82SMSEVSrgnAYNmy4cHbNhUzlntmvWMlDd9cuyd69Flu2nN+b95Wu5ULxRvp39kZaK7xGPYC9e/eyZcuWpfuVSoW//Mu/5FOf+tRFL6SOV4hq1ZR1qhVT1gkG0Y2NJhCPjmIdPwYIdCiEDoWQmRFkPocOBvEa04iODrAsVDSGnJrE3r0TXSmj2zpwt14Pnod15DD2ju3IfB4tLaiWEYUiWitQCoULHuipcZxDB6GQh0QCb/16pDYywvgcvHTaBP5oFHdNKyIawWttR0iBG06QT7RTK7joaQenXCL99X/E5wMdCFLc9HFUdwe17uUD/2LAX6zRF4tw7JjkkUfMn/PRo5BOuziOkUrOZGBwUDKTKSIPj5LyVfAmRxhyryJ3KkuAMj3Bk5SaN1HVgSXphsVSk2WxZLh+xRXQ26uIxyEeF6TTmi98wbnoCdvVGDmLOv5BW1Darl+i4w8m+D/4YGDh2A6f/3z5opNAfUL45xurJoCPfOQjPPTQQ6xbt44XX3yRP/qjP6prAV0uKGU4+Qt0E+340NEo2p8yTJ6REazBk4AwjweCyOEhZKmEjkTwOrogGjVBORBACYn14h58E5Ooxkbca7aA7WD1H8N+8geQy6GlhRDKXF1ojRYaalWEVggNYmAA69QA1FxUqgnaOkC7CE+jojHwhxC1CqoxjUo2IEJh3OZ2KhUoWI14qWZsXEJ+iHamkJUK1pFxfCF7RQ2bM3X1wVx5VKvGP7dUEvh8ZkLXlGk0IyOCXbuMcmYmY4K33w/tvlnmZZETXjeFeZeUb5bm8BhZu5FiwU9usko1ECSR0KxbpygUTDN40VM3nVak09DQoLnqKnNVsH27dUkmbBd1fehbjzh8dFktpsU+R1ub+c5791oXnQDqE8I/31g1Afz1X/81DzzwAPfccw/f+ta3+NjHPsa73/3uy7G2n0+USohcDqrzyNkiKhxBN7cY7mGthhzJILKzIC1UNAKODzF8CnnqJDqZxNu4CSIRRKGAkBIVCiGOHcU5fgzl96OuvJra1VuQJwawn9tmhrGEREhh6vULPEeFhz0/j7J96GIB5/hRE3EdP6q5De1zEJ6H9lxobUa7CuFp3LY2iIaoODHyDe0IJVDBJnx9KZKOixAuOtmAcF00wjSX4SwNm2JzN/PTYtESGJ9PY1l6geoJtZpc2t3HYobuady2zOBWqWTM1C3LqFF0d5vyz87hNFahwvroMMQVkw1XMDFQZN6TCBGmuc2hc71iYkIwMSFobzc7/oYGQ+FMp2HTJkU4fPrnulQTtoszCpw8uaIW05YtHrbtMDJyep7gYlGfEP75xgX1AHbs2MGHPvQh/uVf/uWsctBrgTddD8B1Ebl5w8kHU9aJREm3Js16azVkZhiZzRrjk0gEWcgjhocQ1Sq6MY3XuxadSBjdfeWZxu30FNaL+8Ct4q1bh+rtQw6dQg6eQExNYbbzFlpIrPERWGgSy/l5dK0KSMTsNPbxfsjn0OEoqrEB6XrgeSi/Dx1PghZEYkGmIk1UrCAqEMFtaiUYljjdLYh43MwUKGX8dpVCjo4gC3m8nt6lydvys3tw9xzC3bwZ+4ZrCAZNvX1mRjA9LfBGJwhMjxHqS+PrbCGTgfFxuSTRPDNj7BQnJkypZu1aTaUiOHjQPN/cDH19mvKRU2SHysz6GygE0ziqwtrULE3dQU5l45RmKvQkpmjuDhHpTNDerhfeq4hGl/8JX60a+rmf8/PQA3gtFUpftzFhBVyKHsCKCeCXf/mXz7o/NjaGbds0NjYC8O1v///svXmsZdd15vfb+5w7z/N9w33zq1cjq1jF4qTBktsRIkQBZMMwYHQ7BgJHCewGnD8EqYPYMgLIgAzZcIQYMFoNBUbQEmDEdlvo2K04tmhSokiKxeJQA2t483jfncdzp3POzh+7ipNIk6JIVpXqfv+8qod7z9nv3aq11lnrW9/3n3/qg3wQuOcTwG2pBb1WijI9engbCLz+muGQTL9BfX0PZZo44QiyUUcW9xH2CDc3gTs7j0qldPKwLF0tOg7y2hWMchknl8c5eQpRKWNubUGtqpONx0RJE9otjHIZZRgow0C2mih7hBgNkKUy6mAfY2ijwkGcZBqjZ6FcpQfH0Si2DSNlMsxNE0wlGHoNjOkMHo8He3JKM4BcF2HbOvCbpjaDGQzw/O/foNMz6Ikgvf/p3+L1KLLf/CMMZ0jNiXH43/3PVDwTuK6Wbgi1iqT/rz+j1A2w2c/T+MyvYEeTdLuC7W29rdvr6eFsNqvo9wUXL8LqqkEopKmdc3MO/YMG/Zdu0idAxOgy80vzxJYybG/rp4bFVJP8lSeIGF1mozViX/hV5h/NEI1++P8s3tqLv73HYJpebHt4z/Tmf5r/Z3d6/nDXxIT3iI90CPz7v//7P/WNxngH3JZasEcgpZZFzmTfLGUwGCB3d5HtFq5pwmweF5A723gcF3dyAuehR1CpFPR6yFYTUTpEKRAHexhb2/q6R48ymp3H2NrEfOpJRL+nOfUeL8I0Efv7OvmEQ7geE6PRQNk2ot/TS2CVKkI5OJEYTiSiE4LVZRiIMgglUI6LUCbMTOKJ+vEmosSncjSGYE8VcIIBzekfDiEaxQkGcett2maQnpHCvHqBsOUhOJ8lebBGo71Jo22y0UjQSs0higeweUjsoTyhELRaivUXGlzaX0YlU8hBmZ2rHWqRNKWSpmMmkzA/r9jfh2ee0Qti+Ty3DFkUBwe6nRS1BmRki3OZPcLNXV7dTlEmy/Ky3i2IViocj62RnI2x1L1BZ3CTq1fzH0l1+tZe/O3+/tISP7d7AOP5w53HOyaAhx9++LU/NxoNer0eSikcx2F7e/sjOdw9C8d5va2jbtkfxuOo2yI0t9Hv66Ftp4NrGKhQCKVAbm9Bs4QIJnAefhSVTOprNhrI/T2wR4iuhdzZQtk27lQB+8yDGHu7yMuXdIIQEldKVCyGqtYw99ZQpoETjmMMB3Cwj3BdlGUh6g3tz4vSpuyxCO5giKq16UfSOLkkXjXCF5DIiSncWExLOHj9ejX25HGczlALw41GOP4QLW+G0WEH5XMgVSAcVsQCCnE6h/oviupqnTV3nu5wEfyKsPQSqW4RDo9oLWS5uS0olyWO4xIKZ6kPI1Q2TZrOItVaVuvszILHo7h+XXDpkiCRgM9+1sWyNCX08BBGIz0cLhTgaNQl8k+vsllfYE8eYfaEH8+sIhaDI0dcksMQRxv7pO2rdPHwf/znJUo+8yOpTt/aiz9zxuHyZflzvQcwnj/cebzrDOAb3/gG3/zmNwEwDIPRaMTS0tK4BfRGvFVqwTB1WycY/MnX9no66He7uEKiggFkt4s42AcEKp/HnZ0lfWSWcrmtl7vabbB64NjISlnTWuJx3FgcWatBp6U3f20HpMANBhCOg7F6A9Fq48aSqFgUubWle/y4KHuELFeQ7TauUOAPMYiksa0+ZrOOnZ/EyKfwDy2UEKipSa3x7w+ANCCRwJ5fQAhBPBpkd6eF5QRwYgk8gy7BoMJXyOinD1f34et1QacjMCtFzMN91OQEw+QE4bDCU95n+8UKm4MZBpEkkYii3RZUKpqJ09pv0zocEEwHmDkWpF6Hay8P6dWHzMxLHvlXQba34fJlSbutDVxMUw+C5+b09u7OjsAp1VjxbGNMpggeLXDsmEMkIigUXNLp1z11ny/P861/mHutOv3X/9rm0Ucd5N4utYvbbDJP5uzkB8r7f2s/f29P0G6HiUTuDi2g9/IzjGcAHx7uyB7Ad7/7XZ544gm+9rWv8aUvfYlnn32WJ5988qc+xM8d3iq18HZtnduwLB30LUsnr0gYJUEe7AEClc1iP/Y4xOKa7TMYwMEBxvqefn/PgkZDm4uHIwhAdjqasTMaIkZD3GhMyzDvbuO5/ArK52M0s4Dhq2NurKJ20O0ne4hxWER1e4wwGPni2LEUctglWD7AnJ9FHJ/H12ig1AC3MIubyYDH0E8zyTSjpRV6HZd+2dZKmUpv96a9I+SwiTuTwRYeag3tgdvrCXw+zeTx+6GfzuObzeH3Kw4PFT/+saTbLRBNTaEG0Gtqts9wKCiVdCtnYiJC4ViIGzckTz0FfvqcqT7JcnCfl67P87eVxxgKHx6PIpvV+kCZjMDnc9neluztwcQE+OcTTM9OMTnZIRLRQSebfT3w3NbTz+wJzO+/uTqVe7vwv32Ng5ddJB7+z9P/K//9H+TeNXC9l1733p54Tbvo8mVJPq9fk8lAuXx3BP8Po18/Vii9s3jXBJBMJslmsywsLHDt2jU+//nP8x/+w3/4KM52d+HtpBaiUdQblCnfhE5Hs3e6XZTj6BaQFBiVElTLuKkU9qO3gj7ogWm1gqxWtf5AOgLNBmIwRJmm5uGPHGStinIdaHcgHEalMlCt4nn5ZUS3gzM5iX36QYzr1/BeeBa8PhyPH1mr4BarOEObARKCMVRhGl+rTLi2xWhxGefkMma5BJaFMzuLmpwCV6Ecm14gQ3P6ONg2cr9HMOYhmTEQ8SjpiSTVtV36RpK6TNNY18E7ElEEAlqzx7IEgYBm6DiO4pVXBPW6QSDgYhgKw9ASzP2+oNnUrB4pdZtnMICbN+HaNUkqBR/7GGRrO1zcjnCh82lUb0Aw2SE348PnU6RSOuHs7gpcV5DN6hl7On17WxcCAUU+/84B7O30841nt2m1HErBeWbVJsnWFjs7E+8awN5Lr/tu74ff7ecb4/3hXROAaZpsb2+zsLDAhQsX+PjHP85gMPgoznbnYVm6yn+D1ILKxeCd/BDabW2X2GmjXBeVSuH6vMhSCdls4CaTjM4/8nrQv/Ue42BPe+L6ffrJotcDBohWEzG0tUGKEDoBmSZufgLl9SNXryP29yESYXTkCKLVxvPqZYyb13FicWxfBHdnB9lq4yiFNA08uThqYgpRKmHurjFcXmJw7DjmYRHRbuAsLOFMTjNoDxhWRowyE4yOnSToGZJWLWTMB8KDCoZR4TC9vRo7A4f1ZgHVEMRiimRS0evp5Szb5hZvXrG6Ktjfl/h8Cp9PEQjoBa9+X8ssN5vaQD0a1bTL/X3Jyy9rjZ5cDiYnQUqXixclT+7OEBqOSIoq8XAHVZgmMaV3BkolieNAPK4Xw7JZLc7m9wvyecWZM1Auv/vH/9bq1CnMEIoaZDc2aeGhFp19T33r99Lrvtv74Xf7+cZ4f3jXGcATTzzBX/zFX/Dnf/7n/PIv/zLdbpdPfepTfPWrX/2ozvgmfKgzgNFID28Hb5BaCEf4F1XIWk0tudBqopSLyuZQwwHG4SFiNMJNJHCWV3gTl3AwQO7sIDttXH8AJQWy2wWre0vqoU98IkN1KDCaNRiMcKenUIaJcfMmcnMNbBdndgYnk8H76quIg31sM0jfjCCqJcxaGXNoYUoXYQjcRAqVzyEP9pBti9HCAiRTyOIByuvBmliglZjBbLcxnAHG7BSesycx3aFu4Hu9KCFRXh8db4LaRot218TIxllaitDttmm3BY2GwOPRBuiOA1tb2lhFSvB4tAHMqFhneFilKnM0VIzRSMs5xOOa9rmzoxNBMqkDeCol6HQUV64IGg39qwyHIT4oIstVwnMJ/At5qlWD0cjB75fEYi5TU5KZGedW4HeZnob9/Z+try73dln9xx0uVOZZ/NTka736d+tjv9/X3E196g9jBnAncS+dFT7iPYC3Q6/XY2tri5WVlTvmCvaBJoB3kFp4k33UW6EUotlA7u8j6nX9ZJDPoRwH42AfMRjgxuM4R1YgGnv9fY6D3NtB1uq4pokKhV4P+rU6otNChUK4iQSya5FQA+reME5+AnNtFXHzut7azeRwZmaRnRZceZVRd8ggEEUIhb+0h8dqIO0hEoErBSqXw42nMYr7CKuLs7iIisRx9/boKx/dwjFG2QmiwypBYwhzMzjHTsDIRrRa4PNiO4L20EfFTdGrdAmYI+ILCQIRg3pdIESYZrPz2qJUsQjr67ploLd4wXU1o6exUaf+d1B4gi0AACAASURBVBcw3AFd10fv7ONEJsIMBnqL1+PR27umqfD5dDvo5k2dIGIxvf0bDgucZpv4pacJ06FBjO4DjyKjERIJl0JBMDXl4vMJJiZcCgX90Ha7j/2zcOvfia//YXHZ75cgdSdwL50V7tAQuNvt8sd//Mesr6/zjW98g+985zt8+ctfJvTGffh7CbekFoRj66FoKIR7W2rhnXA76O9saxkGjw83n8cNBpH7+xhrq7jxOPbpB99c6TsOolzCKB2iACeZQgSDyE5b0zAbdTBN3IkpVDSGLBd1u2ZxCdwe8smnMX/4FCoQxJmewZ5ZwlndRPzjDxkFwxihOD4axEubMOghBn1AoDwenMkCTiSMubeDsb1Bf3aJgQgiSiVUq4laehBfIUWm00TZRdTynJaRsB1krcYQD41ugFrRRy8QJ+obkPOW8T4Qp9aLc1AWiIqu0HM53VK5cEFiWTro66GvwO/X7J/bEg7+wybWSGKHcsSsMv1Wmx07QjKpXtve9XgUrZbglVf0OCSZ1O2YQEC7dPl8EPeV6TmKcnAa1WqTFnUKp8NMTOhBcyajRdve+LHe7mP/LNz6d+Lrj3vjY9yLeNcE8NWvfpVsNku1WsXn89HpdPjKV77Cn/zJn3wU5/vgUCwiS02U349KJlHv5mrmupp3v72pA3UgqPXx02nMjU2M1Zu4sTjOyVM4bw369TqydAj2CDeWwInEkM0GxtpNZKWMcMFOpRELS4jDoh7QzswxymTwvPA83pf/bwj7ccMJOkunsattzMsbyGEfTyKEOZsieFBEbqzi2AphD1CuC+Ew7sQkTjCEWN9ieNimXjiOYbgEaiWCiRHikROQySAaDVSjgjO/gLtyFFzFsFin1pQ0+lEcVxKeDDNZgJB1SMWOsWVN4mwJolHd52824dIlLcHQagmCQYXfLwgGBb2ew8GBSaulh8FSKup1gU9kiHiuUGl6qckwyUKIREwPiQcDRbsN6+s6amezinhcVzC2rbX402lFtyuouGlwN8mO9pnNNUn94knMSfVa4H87MtbtPvbPwq1/J77+uDc+xr2Id20Bff7zn+dv//ZvX/vqui6f+9zn+Pu///uP6oxvwoc6A7gdvLfWEc0WKhTGLRRQhoG5vo7odnCjMd3eeaPrh+Mgmg1ErYboWahIFOX1IKtVrc9TKSEGQ5xoDJXNIpstRLernyLyExg3rmNcfRmsIf1YBssTIxT00T/Yx9+t4o34EIaJsgcYm5vQH6AcG3MwwhGgEnFGuWl6the5tobtCyAWCgS8Lr5aGScUQi0sQDINjRYoG3d2AffIETpdQX2zRactIegjEXVITAfxxII01moc1APYoRiRqB7Gdjo6+LVaguFQEY2CzxdiNLKQ0mV726BW41YygN1dLd3s97uYpqTRgITZJK2qND0pukYMv187ct0OzBMTEIsplNIbvIGAIJVyb1FKXRxHkkwq5sJVJtjDmMyRPJplcdF907jmNqf/jZ4Cb+XWvx8e+ttp9jRe2iF+pvCOmj3vF/dLm+JO4F46K9yhFpB8SynlOM5PfO+ehm0jKhWM7Q3odFCRGO7sHMxJjJs3Ma9ewY3GsI+f0E3o27gd9DsdHRWDQVQojLC1eBv1mtbRDwSw85NIBKJW1RTL+Xlot/E8+wyqUqfrT2BHC4iAi1eMiPeKxF1o+BQEYohqCbm/r5fMXBfZ7+OYXjqJSbrZaaTVx3dzE38qgnF+BWk7iHoFRRD75AO6H1KrQ72OPbtAM7dMtenFeraF3xwRn/QxWbAxIgZ1Uqy/2mI4aBOYzBCbNej3YXNT0O0K2m1dlZvm7WpfM35u3BCADsyBgB62uq524vL5dPsnHndJpwWlUox9J0Y4DJ2i3uT1+eDIEUUwKG65cAnCYcjl9FNBqaQYjQSJhGBmRpHLuXg8SZLJJEtLLqb55spb7u0S/OOvIeyRVhj9otbTfyO3/v1y29/IDrqt2y/sEeqy5211+8cY427FuyaA8+fP8/Wvf51+v88PfvADvv3tb/PII4/8TDf97ne/+9p28Sc/+Um+/OUv/0zX+6kxGiGKBxg7W9Af4MaT2Msrent2bRXj0ksQjWGvHNVKYrdh27qiH/TB6oJhattFvx9ZLGpv3X4PTC9uLIrK5hDNJsagj8rmcCYmUM88j/vEc4w8ftr5KXz5KEGfg7TbCKsHroPyBcB0kKtrGM06rhC4Q5uRZTPwBXGzMziZLKHWIZndKzi5HOrhU9rkvVKHgE8PctMZqNewKw1qiWXKsQX6yk9kp0ImVCN0OoSBSWtgcqObZ7DdIUiV6HScAT4ODgT1NYFlKUIhkFIRiwkCAYVSLltbJv2+YHlZD2d3d/XiFmj3rH4fvF7NwLEsQbUqAUUqpbh5U/Dqq9pU5eRJhderdwD6fc3wyef1HKBe13OBWEwwNaWYmHAxTZ0Ilpdd3qqucRu39fTfyWMAPhhu+3u5zxhj3K141wTwxS9+kW9+85tEIhH+9E//lE984hP89m//9vu+Ya/X4w//8A/53ve+RzQa5dd//df50Y9+xOOPP/6+r/meMBggNtYx9nZg5OhFrJOnEcMBxvoa8uIFVDSGvbyip463cTvoDweo0QgB4PGgvF7k4SHy8BAxGILXo+cL8SS4jta79/vpJ3MMbuzi+U//n+ZFTk/hO7aIz2ohPD2UPQJrgECCx0Q02xg3buLg0B96cZsKoUa4qTTmQpZQIoA8OEBsvoQ7O8do5qyWkqhWwOvDPXFct4RKbSo7NtXwKeyJaSJZLwVRJkQRFYvRs4KsbXrp+JP4GJBw9hmmYhx286xdF/R6moXj8+mK2O8XCOGysyOp1wXpNOTzLpubkqtXtTqFlLpN5DgQjSqmp3W7ZHtbEg4rkkmX69cl168LslnFI49oFY1eT/f/ta+uloBoNhW2rU1dFhd14DcMQTwuOHLE5Z32727jtp7+bY+Bt9PT/yC47e/lPmOMcbfiXROAx+Phd37nd/jN3/xNPB4Pvn+JIvke4DgOruvS6/UIBoPYtv0zX/M94dlnEcKPfeacVr5cX8e88BwqGsNZXEIlU68zgd4Y9IXQ5rAAhoEoHSIO9hFdS0s6eL24wRDC5wN7iO3x0zEiDLo2wR9cwNeoEJ5II88u6y1f20KpgI6WowFKGEjHQZWLDLeLDEcGrsdHbNTCJ3twZBoSSZTHg7m1jjro4S4uwdwitJva+1cauCtH6IbyVLcs6iU/MjtD/ESG+bkIgXYZ2h0G/hg3immsHYFMxohkFJl6kUbfz+XhNNaBboWEQgql9J9NU9HtwtWrEr9f8+odR7K9LRmNFDSbeAYlvEEfZjxGJqNwazX2XxlQi0ZIzYfxehXXrkm6XcH0tOLoUYVt69mAELriDwbVrRbT7Q9MMD+vbgVl+VrF7/e/80f85t78NNYX/91PzADk3i5cqyAjaaampn9i2/enhTs1zfq/+V9enwFMvf0M4E5q3owxxjvhXRPA5uYmX/rSl7hy5QpCCM6ePcsf/dEfMTEx8b5uGA6H+d3f/V0++9nPEggEOH/+PGfPnn1f1/qpcOwY4oXLmD9+FhWN4i4uvjnoj0Y66I+GKMMEeev7rouoVl934orHcP0B8PjA68HGoK/89EYBHMNDcGOd6MHTmNEgzvwEsh+CWg3RbmmDduVqeQclGJbbjHZK0KiD4cUbMEgM6ri+KJFzJ2m6EoZDjO11XAT2whIiGtVrs+UDQNKYOsGhyGHt9PBGJckj06wsZzFyKWS1gr29wXY3Rs3KIV2H8FSYaNbHcLfC4R4cuhMYHkk4fNt0TOA4CtNUbG9LBgPJ/LxiZsZhbc3gwgUDIRSDAahmi+ALT5PyDfAbHQ4f/m/ZeMUk+cpzLHmqbFoZXtj7JEMjyNyc7p33+5oxZJp6jh4IuHQ6kn5fU0a7XZiYUMzPK4QQhMOClRXnbXX13oi37+dPv6kdc3sugCkI2grri/+Oqanpn4m2ubcn+OP/uIBtL2Behi/m317nZ+y7O8bdiHdNAF/5ylf41V/9Vb797W+jlOIv//Iv+b3f+z2+9a1vva8bXrt2jb/+67/miSeeIBKJ8MUvfpFvfetb/NZv/dZ7ev87TbPfFS+ukjx3UgvC3B5ij0Y6mA6HmmNYyOiBrm1Dswn7O9piKhiEqQyko7gu9ESAru3T7/F4SLktAjde0rOBfB5mT6PVzLqQiEHIixqOsGwP/aGLrJTgYI9wv4M/GUB6TWg1ID8N50/pKNFsEtve0CI2589pt5N6HbdWodH3Ucx9gp43TthusDwtiX6sgJzMQy6HfVhl/0qJ+jCE7c8SjVhEZ+JYIoS1W6NyrY2bLhBZ8DCvi3g8Hl2JF4uwv69/jHPntMH6jRuae3+bOevx6OOkjS1kuMGOZ4HDWo8Jt8tUbsS1YZCrwyOYgy7Hc01mP5aiXod2W/8qEwn9tV5/nfVTr2vv3WPH9D3CYW3HGH6PH/e1a/paS0ua5tlue8lk3vqiCpgC5uYIbG4SaFcgc+z9/Xv6Ke77ns72LyCTeQcrsrsU99J576Wzwgd/3ndNAK1Wi1/7tV977e+/8Ru/wV/91V+97xv+8Ic/5LHHHiOVSgHwK7/yK3znO995zwngfdNAH3xQU6iKDUSj8Ro7REUiiJ6N6HdQnR2ty3NrkukWChBMMmgPsJpdnGAUKcDnVwT8fYzd6xjFIsNYjEEsgTAs2NiDVArXDDBoj+hXOjhK4q2X8XcrBEq7CFycaBTV69Mr9bAXl2HlJGpkI7aLGJUi4YkcjZUTEA0zqlm0r2xStUJYM2cJT3qYMMpMTDgQT2Kn05QzeYo3OjT+4SrC68GTjuCnR29gUnHjVJ7pEBzu452IIwsxquURtfaIQEDR6ykODw1MUzEzo3AcxeXLBqORwutVWJakVNJLX6GQ3gPo9+FSNUm0k2Qhsk47IHileJ5yzSDUr3AidJVcpsVB/HOsrlr4fAK/X93KbbrnHwopDg8lySScPu0Aegi8suISjWqXr17vvX2+kYjAtr2vBdtIZPgTKpoykiZoKwKbm/RshRVJ4/6MNMD3ct/38pp3wv1CVbwTuJfOCneIBjozM8PLL7/M6dOnAV3Bz8y8/0HX0aNH+frXv45lWQQCAb7//e9z6tSp932994xyGXlQ00E/HkeNRtpV62AfcXCALB4gej3c2Tms3ByDUhtnu48TTxJIQtQcYcqRVvlc30QhULEYznQBeVgEx6UfjNMNTONWB0i/SXBkkeiWkdUSZvEQN+DDiUcwq2UMq8PogdMQCCJ6A9jcwGw3cXN57IceoZ/MsLvWo/tyCSUkoeMrFCY8BIYVLQYXWsSOxih6CpQ2e8iLewSDglAhybA5YGi5rHanMd0hSWefmakgO+1p6m0I2Aq/36VYNCiVdLW/uGizuWnw7LOSQEAH6nYbLEvi9ytmZ8EwXBoNyeqqZpaefjzIfvxRnr3pUO6H8Es/5z4OmU/OUVrNsh2K44lFCZuKQMDFsiRSKiYnXXZ3JVIKHnvMQUqBaQqOHnXftF7x0+Dt1Dvf2nd3p6bp/5vfJLB6lf7S8Q+ErfN2930/rxljjDuBd/UE7na7HB4esrKygpSSa9eusbi4yHe/+933fdNvfvOb/M3f/A0ej4dTp07xB3/wB+95EPy+nwCSQcqlFqJeR/QsROkAuV+EVoNBepJOIItbayFGQ4x0jGAygMe2UMMhwjCgWMRoNlHBIE4ohNHpMKp16AaTDGUAObLx3Vp+8pb2Ud0OcnMLs9vCicTAYyCrNZxMCvfIcUBotc/9XWS/jztVoJ0sUHSydDuKbGkNj2gTfXAeM+xH9Pu40wVcj5+yHefALGB3+oS7hwR9Lh1/GjkaUKyY9IJpggGXKU+Jct3gYJTB69Oc/WpVYFmCYFCzeBoNwdqaZDBQBIOaa1+t6uWrQEBz/kcjqFQEgwFMT7skk3D5spZv7vUgkwmwvGwRCCiKRYkQun/v9So8nteloeNxxfa2wDAEZ87owC+lrvhvPRB+YHi7vnuBHYJ//DUCptBPAF+8+zn790uVeidwL50Vfo48gb/whS/whS984UO7/tvi4AB5ZQ2juIdbqtINJLGSBYQ5ha/fIhyyMI4kMYZDaLdR1hAx6GOUSyhXgT9AP5ZmtFdlaLqocBQjFSMke8SiAppdZKkMnRbm5jrYLnYqiWvGwBqg4hlGD5wDZ4goFZEHRRwX6vmjlDwFut44oVGfyZ0fE/b2if/iQ9QHCmGPGOWnaXQ9FCtRrPQUXndAvL3BaKBwUyl2S4peycY7kWXyAehulTnchldCWeIpg8BAG7O023oEkk677OwInn1Wt30CAcVwCIeHklBIyy57vXrBa3tbU0BnZ/Vg+IUXtOLnaKT7848/7hKJwM2bgm73dSqo16sYDnUSmZpy2djQ73voIc3ddxxd8afTH87H/XYc/zk0Z5+lI4hrN8ac/THue7wnT+CfB3T//kk6tT6DiVnkzBRh1SHr6+Lm4wgniWg2oFIGvx9labct2zVomAnodBD1HjJu4J/JEBYOyiNQhsTYq8PhLpRKeHZ3cD1e7Gweo95EWkOcQgE1NwOdLmJ9FSpVKm6cw8zHGAXihLIB8uEOsSv/qBeKzj0IHhNlmFRTcxTLHnobHsTMFKGYTexgGx8DyjJDwzIJ9IfkT0TpOX52r7Z59bpFcDqBb9aLVde+uJGIYnpa0zlXVwXdrrxV3WtHrmZTEg67zM9r85VuV1CtCmIxOH3apdlUPP+8xLLEazr7Dz/s0O1KSiXJaKQVs19X/rxd8btsbgoqFcnZs5rJMxwKFhddcrkPtw3ydhx/B83ZZ3NzzNkfYwx+SjnouwHvtwUU3tigVy4hBn3tbatcZLOl7RbDEZRycQ+r9CttBp4ISAPvoI0/amDG9OOTUK62ZGw2kcUiatBHbq5jNuo4sQQqHEE2Gyi/D2dhETebw2g2GV1fp1JSVGQWNTtHLG2QKfjxeV3MZ59GDIbYpx4A00Or52M3eQKP9NNpdfEuF3BdQai5T6/WoyKyKK+PQryBzCRZ249Q3esRc+pEp0LU7SidjsDTqpAeFFG5DPt2lsNDzduPxcB1FQcHulWTTmuVzW5X0etJlFIkk7o9tLUlePVV/T7bFmQy8OijNuWySauldwT6fUE2G8CyLOAWOyit2NmBSsXg3Dmb2C2xt5UVl4mJf1l49YPE23Hv5d4uqXaFaiR9T1T/90ub4k7gXjor3AV+AHcD3vcMoLxDtdhAtptaVz8cYuQN0St3cfeKKAEk4oRcC68cQjgEholAobxelMeLsbML7Sa0mhhr64hBH2diUs8IrC4EgzhHVnDjMQZbNcqXD2l1TIxoiMQDE6TTYCRiKAnm0z9EWBb2iVO0VZjDTpj65ApGf0DY0yd46gG6gz6j3UOsgw6DSIrUTIC02WC9GmezFkO6IwrBKo4wKdkpHFerUQc6JfiPf8l+O8qhncD89CeJLiSo1xWViiQY1IlASs16dRytw59K6ar9yhXJ1pZAKd3GKRR0xb+zY9Dvq9e2d0MhnQSi0SBgkUqpWzRSyQMPOORy0GpJlpZcCgX1kQX+d8O99B//Xjor3FvnvZfOCneIBfRzg40NVNWiEcoxNIOYa9t4B7t4MzECMxH9JCC7qEAIXC94vKhwBNVpYWxtokY2qniId30NJ+DFmS4g+gOE7aCCAeyTp2jKJJVLJXqbawRkn/S0n6lPzyCDPghFdMvo6aeg3aa+dJbqIEqzHGY4v0wo3CYxquEsL9KXfuydBu2tPUJTceY+XaBys8n6dS8XmGUi77AYPaRRF+z0MwTDgkRUa/W0WrD5nIVVmSGV8zLb2uVmqU3JTRKLKSYnecsWriIWU/j9iosXJaWSgVJ68Hv0KJw65bC5abCxIRmNNJMnFuM1r99gEBYXYXVVv//YMZeTJx0aDUkq5XL2rPO20sxjjDHGncd9kwAO48t0izcIrV0mEfHiTsYQA4lwhrgygIhEEK5CRcIojwe5u4vc2gJ7hNhYxVMsYifT2CsriG5Xyzsn0lQyRyn1IjjPbRCpP8dUsEXoVAZ35SSYHvD5cQM+zCefpFe12Jk8TcMTQ3X9uEtLBJwuCadIf3qRmhuhv14h7W6x8vAEWxMzbF7s8MqNEbH5KdJHIFqq09wcUEwmCOa95ExQSlfeu7sGhqGYnIsSe6nF9n6SopwiPBkhntXVfq+nq/ZcTtsvui5cvKhlGhxHq16cPOmyvKxYXzfY3JRYlsKyBKmUlokQQu+nZTKKeh2ee05TST/zGYd6XRKLKR58cBz4xxjjbsf90wL6x/+H+hCt09/poDxe8PkRoz4qEMKNxpBWF7m1iRoOEB0L49XLyGGf0fQMRKLIVpuh8FAKzVOOzOMiSa+9QN5axRcycefnsJeOYChw/X5cr4/hPz1NpzyglFxh5AvjjXhxjh4nqCy8VpNmeo6KnSQ6KDNtFjHiEa5UcnQOXEZ2n9yJBMKQHG72MVoNfJkwvlQE1wXbVqyvSyoVQTqtyGY1zbJWE4ScOiGrBpk4PU8SwxDEYg4ej4Hf79JswuXLeoA7GimkFJw96zA9Ldje1r2ael0PbbNZ9aYt4HxeSzTfuCGZmnJ57LEQq6sWhYJOHP+ShfLdgDc+SpsvPI/50ovYZx7EPnf+Dp/sJ3G/tCnuBO6ls8K4BfSzIR5HHNZAoY3ebVtrAvlzyJ0tjO1NcBXqYB/PzRvaVnFhAUcaDLojiuUwjdgJRDZN0qhz6sb/i6ddxYlE4ewJhktHMIYDhD9E2/Fgff8i3bJFO1FA5gL4kyHM0w/gGVoYtT32fXPIiUWyRpk5+xIbzRBPNI7AjS5TuTanPlXg6o0hu1sjkk6ZTMKLPTmBx6M4OBCsr+se/dycIpVyuXnTYH9fG6nkcgKlYlhWgmAQJuIujiPwerWu/tqaLs21jYHg0592CIcFxaJkdxeqVR34p6cVHo/CthUejyCX09TQF1+UpNOKT33KptEw8Png0592eDeTtbsN5gvPE/m3/yO3lwXaf/bv78okMMYYHxbusf+yPwO8Xj319PlwE0nEoI9cW0P0LZQL8uZVPLtF7Ewa+4HTdNsuxUqEFmG8hSyZuSBT3V28l59CNFs4yQTOI49iLyxhWB26Iz+dXpDeE5cZ1No4sTjeQpb4VIr+0dOYtsVodYtGdgb/8Yc4Ea7TvHKV6xtefiCWyKeGnJyuU5dJijU/9oZBtFckFRWM4lkcJVi9CYeHBqmU4sQJl/19uHZNYhhaTTOZFFiWy2AAyaRm+AwG+sfe34f9fQPHUa8ZrH/uc9rcp1aTNJuKSgVsWxuuaB6/fjKYmYF+X/HKK5oa+slP2jSbBj6f4FOfcpic1J7A9xrMl14E28adnELu7+kngXECGOM+wv2TANJpbG8EY3cb88WLuLjQ62G+/DKi12MwO0/l1C9wWITBlpdQ3CB7Ls7sZAy5sYr59CVEr4+bz+OcPoczN0f/sEW96qPX9TK6cAnRbSIjQWJLGYZTc6jjx2lV+pg3NuhNTpD53EOM6h02nrzO9ZLAmJ1l9gGXrNVko57k1UaCaMRlIVYhIg2auSTFqpe1C5qKOT+vmJx0uHrV4Mc/FoRC2pvX49F6/IMBLCzoAD8cCkYjxdaW1u8fDhW9niCfF3zmMw6uK2k2Je22olwWKCVZWNAVf6+ncBzB9DTYtsvlyxK/Hx57zKHTMZBS8IlPOHwUKt4fJuwzD4JpIvf3wDT138cY4z7C/ZMAVlcxN3e1i9fhAd5rr+J4vGxPnqUyiENrQGw4ZH5J4pvPQiSMfPll5IVVGNowNY198gSt2ATWfpfGboBeTeK7/hJeq0kiamAv5ZBLSzTzy/Q7NuFr62SXU3gePserL/S48RcbmGpA7uwkR895KN7ocmMtRDA9QX5OYVgtPN0mXX+C/U6WV5/vEQ4rTpxwqNXg5k0D19UiarkcuK6LZRmEQi6nT7sUiyaNhsR1FVtbWvLBsjRXf3ZWcP68w3AoaDYlrZaWbZBSsLSk8PkU7baWcCgUwHG0eYsQgtOntQcACB5/3CEQuNMf5gcD+9x52n/27+/qGcAYY3yYuH8SQL+Pcf0azs4hu/55alP/DWI0Iul0OJ4qYWbjuBPzYHgwXryI2N5ESQUz8zSPn6dFjN5Bi2YjgGraRDeeISfaCNNELufpzRyhkVlCuYrJxjqhQoLV3DkuXnIwntxlJtXm5CezlPsxdq93YM8mPZ9mMSRQvT7eWoWSFeL6YYFuV3DmjObeX79u8MwzBoGAIhrVw9h+X/PwFxYgEHA4ODA4ONDfPzgQ9Hpa9sFxtATx6dMuwyHUapJKRVGpaA/eU6e0vWK9rtk/09MKpRSrqxLH0c5bt4Xazp51CYXu9If4wcM+d34c+Me4b3HfJIDVv3qJXd8S5tJ5JmSJk6FdDJ8Hlc7g5o7i9nuYLzyP3NnBCQTpLJykvvgQPdukvdbF8nsIWIL8/lPgunh9I1Q8SX3mAaz8IqGI4Eh/nfIgyo8G52hdcJlQu5zPtOjPZlmrTrL7UpNEssPiwzFGjsQdOtg7ZQ4rJjcbBYIhWFlRDAYuq6tQrUricUUmo0XT2m1tlH7mjE2jIalWDUzTodOBWk1r8TQaAsOA48ddTpzQbaFqVXBwAPW6Xt46c8bFMLQVgpQwNeUiJayvS4ZDmJ/XbSXTFJw8+f4VOscYY4y7G/dNApj41AqZ8i5KWgi/D5WZw0mmEPU65tNPYezt0Q2lqS1/it7SSRo1hb3WYRQKkPUNye48jaGG+EWXfihOaeljMDtHvmDgXF9j7UqIZ9SD+PySo/E9wp4GxWGaC51jBGpN5qcriIUUfduD1XMZ7JXZ33LZt3NMzhg8dlw7bj33nO6353KacjkYKDodzch56CGH1VWDmzcNPB5Fq6Xodg1aLe2j6/PBww+7LC5qcbeDA63Y2WwKEgnFKSGUaQAAHj1JREFUww+7qFsGMK4rbgV+zfcfjVymp8Hv1zz/kyfdN1kjjzHGGD9/uG8SQMgzpBFPoLI5VDgExQPMf/4nBntVqrEZrOOfpZo7Sr8+gJsdwlMxEvE+wYNX8ToWwmrT8qaonf4MmbMFMl4fpafXee7HPmqpMxQW4DGjTH+7wtZOEityjOlog0cmirT9aXpOEqsBnb0W5fUOw2iaxdMeclJLK6+uGqTTkM262La4ZYYiWFlRxOM216+bXLyoq/JWSzEcSsplF8sSRKOCxx/Xgd+ybhuxCyxL6/Y/9pgO/PW6ln2YmnLx+XTgHwy0MJu2XFQcO+aSzd7pT2uMMcb4KHDfJABOnEAVa4idHYZPrtMq9xmkJ2ic/gVq4Vlk1yLerBGZiBFuWwR3r0Lfwjrs0o3GMX/pvya8PEu9HODSf1mnM/TgO3ac2QdNZqpVDl+pccUJ4p0/yrEjLYz+PlXS7JGkWxVUdnr0d+uEpyIc+6/y7O1JXnhBApquGQ5rP9xeT2/S/sIvwP6+y9qaZGeHWy0g3ecvl/XXbFby8MMuc3OvK33u7Gjd/slJePBBvelbrb4e+INBxeqqrvgTCUU+ryv+20JtY4wxxv2D+yYB9K+sUXn6Ok7Lop2ZoX7mAaxwjoSsM+Wt4pmKEGm08a5fplkZ0G1bePMxwr/2GfY982xVggz+elPz7R9ZZCIRoHSjxc2/O8QT9pN/aI5jqQGtrQOqzTgdUcCyBPtbQ4KdChNzJr5fzHPlqsm1f9AyCtmsg21LWi19xpkZxdGjDpubkh//GCxLIYS65cylqFYlgwFMTAiOH1cUCi6djuDqVV3xCwFTU3qYCzrwj0Y68IdCuuIfDvUwORjUWj5LS+quEmobY4wxPjrcNwlg55lNDswpmg+ewZsKkDbbFPwVvPEAgVqJzvNX6bZ72HaH5Gyc0uOfYF3Msb8TJlrfJhM9IPRLs3QJsXrJwlNcZWIaHvm1afo9RXO9wlo1SN9XoFaF4j5MmCUemh9x4GR58YaX0brW5V9YUDQagnpdEgwKHnzQJZvVfPunntL9/UgEul1JowGNhjY3mZvTw91sVtFqCV56SbK3J/B4tGVjPq9ZPfoJQVAouEQimtUD3PIA0K+Zm4P5eXcc+McY4z7GfZMAmr/0q3hHFkuyiz/k4A+btK7sYl2qMhhaJGnQy+TZST/GrmeBfjdCbrTL6eAe6uECa5U49Yt9ssNVHl62Cf+rAnuHXnZermK5PvqBSYqHgl5PMhutsrLc4eJOnr+76CORUMzMuPR6glJJG5Qkk3D0qMIwtPzy1ataUsE0tQ9vtQrFoq7qZ2cVDzygSCR04njuOUmxKIhEXJaXIZO5HfgF/b6u+GMxxY0bklJJ9/elVPh8UCjA0pI7FmobY4wx7p8EsFwYMuj1aLRMGld38bTrxI0m4UGPQ3+Bl1KPUYnM4UnFmFa7xNU2e8YMF5oLGNeGHA9d59MPDKhE5tivhxg8UwHZp2HkOaxIPB7BkakW1OpcWE/zEkmmp+FoyqFWk2xvCwIBWF52WVzU5uovv6w3fP1+rdDZaknqda3RE43C8rLi9GmXYFBTOJ96SlKrCRIJh5UVrcfj80G5LOh2BdPTLvG44uZNSbWqF8aUEvj9islJxZEj6p7T6xljjDE+PNw34WD/UNB/dZ+EbBIetqhXh2yGp6hljzDITOGfTnHU3cet7vBqu0DbO89EasBnF6/hdftsuHP8uBtBbdcRToWSm6LW8pJOKx56oMfGhQbPPBXEk5lh5phiMNCbtrYticfh/HlFPu+yuSn5wQ9MXNclHFY4jjZpaTbBsrTM8vHjeghsWS7lsuBHP9LLXem0w4kTr+vxVyqa8TMzoxlAN25I6nUd+B1HP2lMTipWVly83jv9CYwxxhh3G+6bBDAzWGWjXqRetajG5+jPz8LEJNHlNJH2AXuXXuUiE/imZ3jg4wMm7ZtUt7qsNaapqiTBYQOnXeJgkGQo/MzOOkxMDrj6gya7L0JiKcfRR7Xa5s2bEr9fkc3qFk0o5PLqq7rNI4QO0JalBdo6HUG/D5GI4OxZxfHj2jS9XIZ//mcD24apKYdIRCeSQEBRLuvAXygo5uY0U6he18PdXk+glDZDX1lx8fvv9G9+jDHGuFtx3ySAzeeK7HlmcU/k8U8lyR6ZpHGjzM73VnESCaYfPs7njrmMNjYpvtTkx8YkvfAivlEbo7LPei+FP5lkZcVhY93hlSdbhM0+06eTDF0PxaLk8Kq2Wjx+XFsgKgWXLklaLQPThGhU0WwKtrYEo5Hm+mvjdXVLdgHW1gSXL0tCIZiedvH7BfG4VvusVgW7uzrwz887rK7qJ4doVC+LDYeaAbSy8vMp2zDGGGN8sLhvEsDo8U+TlC69yVk2Lzbgu6vEZ8N87H9YIhQWHLx4wOp/qlIxcsjJWdxuD3trn0MnTm5xiuMJl+eeE+xc6zMXb3LyfIxKJ8b1NW2unkgojh7VbZ56XfD885J+H0IhCIeh2RQcHgpsW9Hva7nmhx7SCp+OA9euCa5d008OS0suiQSA3g+oVHTgn5pyOXfOZXVV0mrp7d52W9DpaJOWo0ddotE7/ZseY4wx7hXcNwmglD3FxvM7pG6ucvahALlfX6Blmey+XKF7s8QolkAVTkFvSOVSEUJBph+YwC4JLl6UBI0hJyeKuPkgu51JNq/p4D456TI3B7GYy+6upnGORjo4e70upZJBvw+uq1U502nFsWOK2VmXfl/w4ouCjQ1JJKK3cEEQjytmZvRi1/Xrerh77pzL+rrkxg19jUZD0GgIcjk4cmQs2zDGGGP89LhvEsB8/wpnfsnFyc+zV/Hz/BN1husHBPNB7JWjdDvQeqVKMiuZPp/l0mUPN74P0xNDPnmkSLXp4UZ7kl5ZEInoQe3UFHg8ips3BRcvah/EVAp6PZftbYnjGK/57GoTF8X0tEuzKfjhDw0ODiSZjMPx4y5K6cAfiWhht1de0e2hhx5yWFuT3LypHbnqdXHLAlJv72Yyd/gXO8YYY9yzuG8SQPrRI7x4ZUT5H1tEauv443544Aib+wbmjRr5nEIsJ3j5qhexrzi6YrOcqFIrKS40ckiPSTyuOHpMD3cHA92yqVYFpqlIpRStlq7apRSMRvq+uZyu+CcmNKPn+983qFQk+bzDqVMOo5EgGtWtm0ZDcOWKFn77+Mfh+efh+nXJxIS+9v6+NmY/c2Ys2zDGGGP87LhvEsCzPxzh3b1M3m9SLSyxVvaS2a2zEB1wrZ/hqWteslnFI484tPbadFZ7rDtJ/FEfuZTL7KxLLOZSqwmef96g09G9/XRa6+vfuKH59koJHEfdGsYqcjmXrS3J975n0OkIJiZcTpxwsG092E2lFLWa4OpVST7v8sgjWgriyhV97U4H9va0PtD58y6Fwp3+TY4xxhg/L7hvEsCM3OXFwAJtJ0S+X2PWU+PlvTRXimmOHnUpLLiUdgZsPVOjKyN4U3kWczqAh8OKnR3B1au6nx+L6ZZOsSjY3ZX4/S6mqQ1YZmcVy8su6bR+Gnj+eYPhEKamtObOcKgDfyaj1TkvXTLI513On7fZ3pZcuyaZnNQWxhsbkmBQceaMy/z8WK9njDHG+GBx3ySA3eAK0fAB9uY+F2sJorkkDz7q0u06VA8dKpdruB4vZj7PkYIgldI6OZubkmIRbFuSSrkEArC7KxmNtCevxwNCCP7/9u49KKrzjOP49+wuIMqdBURFQQPiJTGaTCSJl9ZEahRLqY0xbaytycSMNiY6rZfUBuukEdpq25iOrf847WgnpNppZ8ykjqnRpnipsRqMQSRAUOQqIBcF9nLe/vHGNSIajegB9vn8x23Pb5ed9zl7zvs+b1KSvhk7YACcPGnj4EFQSi/ECgoCj0fP4U9M1AP/Rx/ZcTrVNQN/cLC+1ON06umkI0Yo7HarXz0hRF/kNwWg6ng1Z87bGDomnkfGK6qr4XQR2C/UE2BXBA+NZsgwOxERJpcuQXGx7q1jt4PTqVflfvqpXl0bHKwv9SilGD4cUlP1HP6PP7ZRWQmBgXrKZkAAvks9cXF64/Zjx+zExCgeekgP/EVFNuLiICREUVNjwzBgxAjF5MnQ2KisftmEEH2Y3xSAMdMGEl7RSm2tvtEaeKmRENVByPAInIODCAtT1NQoTpyw0dho0K8fxMfrdguFhTaCgvTWiC6Xgc2m++qMGqVoa1McP65vBvfvDwkJJgEBCrdbD/zx8XrrxWPH7ERFKdLSrh74w8P1PQSXy0Ziop7Lf7kpnBBC3El+M8ycPg11dTYC3G04XY2EJYYSmRBOQIDuullcrBdUDRigN0mprDQ4d85GSIgiKkrhduuz/2HD9KrdxkY4eFAvyAoN1d0+7XbweAyCgkySkvQZ/9GjdqKj9Rn/uXO2z2/2QkSE8nXvHDQIxozxfr4rlxBC3B13vQD89a9/Zdu2bb6vKyoqyMzM5NVXX72jxw00XES1VxGTEEjI0Di8XoPKSht1dQqXyyA8XN+4LS/XO3CFhSliYvRNW4dDb5ySlGRy7pzuytnRoW8GDxlifr53r8GAASbx8Xp3rqNH7YSHKyZO1AP/qVP6jD8s7Moirvh4eOQRk9DQO/rUhRCiS3e9ADz55JM8+eSTABQXF7NkyRJ+9KMf3fHjpsY30zI0hgvNNkpLDRoa9MboEREKj0evsgU9uyc0FNxuPRMnOVn39Skrs7F3rx2bTd+ojYyEwECT9nYbAQEmw4frgf/4cTshIXrgr6qyUVioB/6ICN0HqL7eICoKpk3Tl5SEEMIqll4CWrt2LcuWLSPqLvQxaApw8umnF2ls1HMpo6IUzc0mRUU2AgIUgwaZeL16P13dktkkOlrvprV3r43gYN3oLSBAD/wdHTYCAw2SkvRN42PH7AQHKyZM8FBba+PUKTvR0Xre/sWLBlVV+nLS5Mn6k4YQQljNUEpZMtXkwIEDbNiwgZ07d96V4/35z3qeflQUVFVBdTWEhkJMjP5+e7vehGXkSN3j55NPdEvmkBAICtKbsvfrBy6X/vnQodDWBidO6J+PH6/bO1+4oNtBDBig+/tXV+uf338/sohLCNGjWFYAli5dSnp6OhkZGbf0d/X1rZjmrUc+ciSUDz+8xMWLBmFh+saux6N78UdEwIgRJoahPm+xbBAefnn3LIN+/fSlnv799WKujg4oLNR7944e7aWuTt8MDg/X6wBcLoPz5xVgY9Qok2HDbn0RV0xMKHV1Lbf8PK3Qm7JC78rbm7JC78rbm7LCV89rsxlER4d0+TNLLgG5XC6OHDlCTk7OXTvmyZP6rDw21sTj0b34IyIUY8fqAf3UKT0VMzxct3Y2DIPgYN3zx243SE01cbvho49s2Gxw330eGhpslJTotQODBys8Hqit1TeEk5MVycle2XtXCNFjWVIAioqKSExMpP9dnPc4YoS+RKObrymGDdObqHzyiQ2PB6KiTDweUEpvtt7WZmAYBikpeuA/ccKGacLYsV4uXLBRVmYnPNxk0CC98Ut9PbS22khMNBk92iurd4UQPZ4lBeDs2bMMHDjwrh7T69Vn/IMG6WmYBQV2HA6F06kHe7fbds3A7/Xq1b1utx74m5ttlJfrM/64OH1Zp7kZGhttDBpkMmmSl4CAu/q0hBDiK7OkAMycOZOZM2fe1WOmpkJRERQU2D5f5eulpcWgrc1GWJiX1lY7hqGnfSoFhYUG7e02xo710Nxs4+xZfcYfEwN2u74BXFene/ZMn+4lOPiuPh0hhLhtfrMS+MQJPStn2DCTCxcMWlttREV5aW62o5ThG/hPnTK4dMnGmDEe2tqgokIP/FFRel2Ay6WoqtKbwkydKlswCiF6L78pAMnJUFAATU164L9wwY7Xqwd+w4CiIoOWFj3wt7fDuXN2IiPNzxd86aZulZUGDodBWpqX6Girn5EQQtwevykA9fUQE2PS2KjP+FNSdAfP4mJobLQzerQHj0dRVaUH/qgoCAoyME2orTXweg3GjfPKTlxCiD7DbwqAaYJp2khONgkMhJISqKnRA39srJfa2suXehQBAQY2m14I1t5ukJrqJSnJ6mcghBDdy28KwMiRcOGCoqQEqqvtjBzpITLSy/nzVwZ+m80gMFBv0djcbDBihElysilz+YUQfZLfFICzZ6GgwM7IkV4iI73U1+uBPzJSAQYhIXrgr6y0kZCgSEvzSk9+IUSf5jfntpGRuud+Y6Ne0BURcWXg93j0/r2BgTBtmpd77zVl8BdC9Hl+M8zV1+uZPJGRCtM0iIxUtLQYlJYaREbClCkmIV23yxBCiD7JbwpAVJSexhkVpbh0CUpLDQYMgIkT9VRPIYTwN35VAJqbFZ99ZhAQAGPHmjKlUwjh1/ymAJSV6dYNycm6pfOttmcWQoi+xm8KwNChkJoqUzqFEOIyvxkOk5KQwV8IIb5AhkQhhPBTUgCEEMJPSQEQQgg/JQVACCH8lBQAIYTwU1IAhBDCT0kBEEIIP9XrFoLZbF99Ce/t/K0VelPe3pQVelfe3pQVelfe3pQVvlreG/2NoZRStxNICCFE7ySXgIQQwk9JARBCCD8lBUAIIfyUFAAhhPBTUgCEEMJPSQEQQgg/JQVACCH8lBQAIYTwU1IAhBDCT/XZAtDa2kpGRgYVFRUAHDhwgNmzZ5Oens5vfvMbi9Ndq3NeALfbzYIFCzh8+LCFya7VOWteXh4ZGRnMnj2b1atX43K5LE54tc55//KXvzBr1ixmzpxJbm4uPWkxfFfvA4Bt27Yxf/58i1JdX+e8q1evJj09nczMTDIzM9mzZ4/FCa/onPXYsWPMnTuXWbNmsXz58h79vt2/f7/vNc3MzCQtLY1Fixbd/kFUH3T8+HGVkZGhxowZo86ePava2trU1KlT1ZkzZ5Tb7VYLFy5U+/btszqmT+e8SilVUlKinnrqKXXvvfeqQ4cOWZzwis5ZS0tL1fTp01VLS4syTVOtWLFCbd261eqYPp3znjlzRk2fPl1dvHhReTwe9dRTT6kPPvjA6phKqa7fB0opVVxcrCZPnqyeeeYZC9Ndq6u8GRkZqqamxuJk1+qctaWlRT366KOqsLBQKaXUsmXL1Pbt2y1OecX13gtKKVVbW6see+wxVVZWdtvH6ZOfAN5++22ys7OJjY0FoKCggGHDhpGQkIDD4WD27Nn885//tDjlFZ3zAuzYsYPnnnuOcePGWZjsWp2zBgYGkp2dTUhICIZhkJKSQmVlpcUpr+icNyEhgXfeeYf+/fvT3NxMa2srYWFhFqfUunofuFwuXn31VZYuXWphsq51ztvW1kZlZSWvvPIKs2fP5o033sA0TYtTap2z5ufnc//995OamgrAmjVrmD59upURr9LVe+GyX/7yl8ybN4/ExMTbPk6v6wZ6M37xi19c9XVtbS0xMTG+r2NjY6mpqbnbsa6rc16AFStWAPCnP/3pbse5oc5ZBw8ezODBgwFoaGhg+/btrF+/3opoXerqtQ0ICODtt98mNzeX++67zzcIWK2rrBs2bGDOnDkMGTLEgkQ31jnv+fPnSUtLIzs7m9DQUBYtWsSOHTuYO3euRQmv6Jy1vLyc/v37s2zZMkpLS5kwYQKrVq2yKN21unovAHz22Wf897//ve7Pb1Wf/ATQmWmaGMaVlqhKqau+FrevpqaGBQsWMGfOHCZOnGh1nC81d+5cDh8+jNPp5M0337Q6Tpfy8/Opqqpizpw5Vke5KQkJCfz+978nNjaW4OBg5s+fz/79+62O1SWv18t//vMfli9fzt/+9jfa2trYsmWL1bG+VF5eHt/97ncJDAzslsfziwIwcOBA6urqfF/X1dV1+dFKfDUlJSXMmzePrKwslixZYnWcG6qqquLo0aMAOBwOZs2aRVFRkcWpurZr1y6Ki4vJzMxkzZo1fPzxx7z88stWx7quoqIidu/e7ftaKYXD0TMvMjidTsaNG0dCQgJ2u50nnniCgoICq2N9qX/961/MnDmz2x7PLwrAuHHjKCsro7y8HK/Xy65du5gyZYrVsfqE1tZWnn32WV566SUWLlxodZwv1dLSwk9+8hOam5tRSrF7924eeOABq2N1af369bz77rv84x//4LXXXmPs2LH89re/tTrWdSmleP3112lqasLtdpOXl9ejrqt/0aRJkzh58iRVVVUAvP/++4wZM8biVDfW0NBAe3s7CQkJ3faYPbM8d7OgoCBycnJ48cUX6ejoYOrUqcyYMcPqWH3Cjh07OH/+PFu3bmXr1q0ATJs2jZdeesniZF1LSUnh+eefZ968edjtdh588EF++MMfWh2rT0hNTeX555/n6aefxuPxkJ6eTkZGhtWxuhQfH8+6det44YUX6OjoYNSoUaxcudLqWDdUUVHBwIEDu/UxZUcwIYTwU35xCUgIIcS1pAAIIYSfkgIghBB+SgqAEEL4KSkAQgjhp6QAiB7l8OHD3Tp18GYfb9q0aZw4caLbjvtFCxcupKGh4ZaPU1hYyOrVq+9IppycnB7XZVbcfVIAhLjD8vPzb/lvTNPkpz/96R1b+btkyRJee+012tvb78jji97BLxaCid7l0qVLLF26lPLycsLCwli3bh1JSUm4XC5+/etfc+TIEbxeL6NHj2bNmjWEhITw/vvv88c//hGXy0VDQwPf+ta3rhk8P/zwQ3784x+zceNGJkyYcN3j7927l82bN+N2u+nXrx8rV65k/PjxbNq0iXPnzlFXV8e5c+eIi4vjV7/6FbGxsRQUFLB27VrcbjdDhw6lsrKSVatW8fe//x2ABQsW+HrN5OXlkZ2dTUNDA5mZmSxbtuyaDO+++y5DhgwhLi4O0J8cMjIyOHToEE1NTTz33HP873//4+TJkzgcDjZv3kxcXNxN/15oaCjjx48nLy+PBQsWdNe/TvQ2t91QWohudOjQIZWamqqOHj2qlFLqrbfeUt/5zneUUkpt2rRJ5eTkKNM0lVJKbdiwQWVnZyvTNNUzzzzj649eXV2tRo0aperr69WhQ4fUrFmz1MGDB9Xjjz/u6//e2de//nVVUFCgysrKVEZGhmpoaFBKKXX69Gn16KOPqosXL6o33nhDPfbYY6qlpUUppdSiRYvU7373O+V2u9WUKVN8e0wcPHhQjRw50rePQ0pKiqqvr/cdZ926dUop3dd97NixqrKy8po8L774otq5c+dV+V5//XWllFLvvPOOSk1N9T2XxYsXq82bN9/S7yml1O7du9X3vve9m/m3iD5KPgGIHmfkyJG+M/SsrCzWrl1LS0sL+/bto6WlhQMHDgB6x7To6GgMw+APf/gD+/btY9euXZSUlKCUoq2tDYDq6mpeeOEFnn766S9t/Zyfn09tbS0/+MEPfN8zDIMzZ84A8NBDDxESEgLA6NGjaWpq4vTp0wBMnToVgLS0NJKTk697jMv3JGJiYnA6ndTX1xMfH3/V75SWlvL973//qu+lp6cDuuum0+n0PZehQ4fS1NR0y783ZMgQysrKbvh6iL5NCoDocWy2q29NGYaBw+HANE1eeeUV30B78eJFOjo6uHTpEllZWTz++OM8+OCDzJkzh/fee8+31aPdbmfLli0sXryYGTNm3HCTHdM0efjhh69qulZVVUVsbCx79uyhX79+V+VSSmG326/ZVtJut1/3GF/skHn5MTrr6vtfbAEcEBBw3ce/2d9zOBzXvNbCv8h/X/Q4RUVFFBYWAvp6+QMPPEBwcDCTJk1i+/btuFwuTNPkZz/7GRs3bqS8vJzW1lZefvllpk2bxuHDh32/A/pMe8KECaxcuZIVK1b4Phl05eGHHyY/P5+SkhIA9u/fzze/+c0b3iwdMWIEgYGB/Pvf/wb0DnSnT5/27Tlht9vxeDy39BokJSX5PnXcKRUVFQwfPvyOHkP0bPIJQPQ4w4cP58033+Ts2bNER0eTk5MDwOLFi8nNzSUrKwuv18uoUaNYtWoV/fv352tf+xpPPPEEgYGBpKSkcM8991BeXn7V2XBWVha7d+8mJyeHn//8510e+5577mHdunUsX77c189+8+bNDBgw4Lp5HQ4HmzZtIjs7m40bN5KYmIjT6fR9WpgxYwbz589n06ZNN/0afOMb32DPnj13dDOYDz74QLri+jnpBipEN8jNzeXZZ5/F6XRSVVVFZmYm77333lfeb9jr9fLtb3+bLVu2+GYCdafW1lbmzZvHzp07CQoK6vbHF72DFAAhusG2bdt46623cDgcKKVYsmSJ72bsV1VQUMD27dvJzc3tppRXrF+/nqlTp/LII490+2OL3kMKgBBC+Cm5CSyEEH5KCoAQQvgpKQBCCOGnpAAIIYSfkgIghBB+SgqAEEL4qf8DHdnLKYvhADcAAAAASUVORK5CYII=\n",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Make scatter plot of 1975 data\n",
"_ = plt.plot(bl_1975, bd_1975, marker='.', linestyle='none', color='blue', alpha=0.5)\n",
"\n",
"# Make scatter plot of 2012 data\n",
"_ = plt.plot(bl_2012, bd_2012, marker='.', linestyle='none', color='red', alpha=0.5)\n",
"\n",
"# Label axes and make legend\n",
"_ = plt.xlabel('beak length (mm)')\n",
"_ = plt.ylabel('beak depth (mm)')\n",
"_ = plt.legend(('1975', '2012'), loc='upper left')\n",
"\n",
"# Generate x-values for bootstrap lines: x\n",
"x = np.array([10, 17])\n",
"\n",
"# Plot the bootstrap lines\n",
"for i in range(100):\n",
" plt.plot(x, bs_slope_reps_1975[i] * x + bs_intercept_reps_1975[i], \n",
" linewidth=0.5, alpha=0.2, color='blue')\n",
" plt.plot(x, bs_slope_reps_2012[i] * x + bs_intercept_reps_2012[i],\n",
" linewidth=0.5, alpha=0.2, color='red')\n",
"plt.savefig('../images/bootstrap-plot.png')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Beak length to depth ratio\n",
"The linear regressions showed interesting information about the beak geometry. The slope was the same in 1975 and 2012, suggesting that for every millimeter gained in beak length, the birds gained about half a millimeter in depth in both years. However, if we are interested in the shape of the beak, we want to compare the ratio of beak length to beak depth. Let's make that comparison."
]
},
{
"cell_type": "code",
"execution_count": 55,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1975: mean ratio = 1.5788823771858533 conf int = [1.55709396 1.6006049 ]\n",
"2012: mean ratio = 1.4658342276847767 conf int = [1.44422244 1.48826439]\n"
]
}
],
"source": [
"# Compute length-to-depth ratios\n",
"ratio_1975 = bl_1975 / bd_1975\n",
"ratio_2012 = bl_2012 / bd_2012\n",
"\n",
"# Compute means\n",
"mean_ratio_1975 = np.mean(ratio_1975)\n",
"mean_ratio_2012 = np.mean(ratio_2012)\n",
"\n",
"# Generate bootstrap replicates of the means\n",
"bs_replicates_1975 = draw_bs_reps(ratio_1975, np.mean, 10000)\n",
"bs_replicates_2012 = draw_bs_reps(ratio_2012, np.mean, 10000)\n",
"\n",
"# Compute the 99% confidence intervals\n",
"conf_int_1975 = np.percentile(bs_replicates_1975, [0.5, 99.5])\n",
"conf_int_2012 = np.percentile(bs_replicates_2012, [0.5, 99.5])\n",
"\n",
"# Print the results\n",
"print('1975: mean ratio =', mean_ratio_1975, 'conf int =', conf_int_1975)\n",
"print('2012: mean ratio =', mean_ratio_2012, 'conf int =', conf_int_2012)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The mean beak length-to-depth ratio decreased by about 0.1, or 7%, from 1975 to 2012. The 99% confidence intervals are not even close to overlapping, so this is a real change. The beak shape changed."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Calculation of heritability\n",
"- Heredity\n",
" - The tendency for parental traits to be inherited by offspring"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### EDA of heritability\n",
"The array ```bd_parent_scandens``` contains the average beak depth (in mm) of two parents of the species G. scandens. The array ```bd_offspring_scandens``` contains the average beak depth of the offspring of the respective parents. The arrays ```bd_parent_fortis``` and ```bd_offspring_fortis``` contain the same information about measurements from G. fortis birds.\n",
"\n",
"Make a scatter plot of the average offspring beak depth (y-axis) versus average parental beak depth (x-axis) for both species. Use the alpha=0.5 keyword argument to help you see overlapping points."
]
},
{
"cell_type": "code",
"execution_count": 69,
"metadata": {},
"outputs": [],
"source": [
"df_fortis = pd.read_csv('./dataset/fortis_beak_depth_heredity.csv')\n",
"df_scandens = pd.read_csv('./dataset/scandens_beak_depth_heredity.csv')"
]
},
{
"cell_type": "code",
"execution_count": 70,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"