6a.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1><center>Colour-Colour Redshift Plot</center></h1>\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "# Complete the following to make the plot\n",
    "if __name__ == \"__main__\":\n",
    "    data = pd.read_csv('opticaldatafinals_SDSS.csv')\n",
    "    # Get a colour map\n",
    "    cmap = plt.get_cmap('YlOrRd')\n",
    "\n",
    "    # Define our colour indexes u-g and r-i\n",
    "    u_g = data['u'] - data['g']\n",
    "    r_i = data['r'] - data['i']\n",
    "\n",
    "    # Make a redshift array\n",
    "    redshift = data['redshift']\n",
    "\n",
    "    # Create the plot with plt.scatter\n",
    "    plot = plt.scatter(u_g, r_i, s=0.5, lw=0, c=redshift, cmap=cmap)\n",
    "\n",
    "    cb = plt.colorbar(plot)\n",
    "    cb.set_label('Redshift')\n",
    "\n",
    "    # Define your axis labels and plot title\n",
    "    plt.xlabel('Colour index  u-g')\n",
    "    plt.ylabel('Colour index  r-i')\n",
    "    plt.title('Redshift (colour) u-g versus r-i')\n",
    "\n",
    "    # Set any axis limits\n",
    "    plt.xlim(-0.5, 2.5)\n",
    "    plt.ylim(-0.5, 1)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1><center>Overfitting Trees</center></h1>\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Depth with lowest median difference : 13\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEKCAYAAAA4t9PUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xd4VGX2wPHvSYHQEiCEGpAixRCSEEJREKwURbCgwtoLiGLZZdVFf7q6rK7oWrCgCDZ0UWRBFFFsC0pRgdA7BAgSOgEhdJKc3x/3BoeQMpCZTMr5PM88M3PnvXfOjDgn7/vee15RVYwxxpizFRToAIwxxpRulkiMMcYUiSUSY4wxRWKJxBhjTJFYIjHGGFMklkiMMcYUiV8TiYj0FJG1IpIiIsPyeL2iiHzqvj5PRBq72y8XkYUisty9v8Rjnx/dYy5xb7X9+RmMMcYULMRfBxaRYGAUcDmQBiwQkamqusqj2V3APlU9V0T6A88DNwJ7gKtUdZuIxALfAg089rtJVZP9Fbsxxhjv+bNH0gFIUdWNqnocmAD0zdWmLzDOfTwJuFRERFUXq+o2d/tKIExEKvoxVmOMMWfJbz0SnB7EFo/naUDH/NqoaqaI7AcicXokOa4DFqvqMY9t74tIFjAZeEYLuTy/Vq1a2rhx47P6EMYYU14tXLhwj6pGFdbOn4lE8tiW+we/wDYi0hpnuKu7x+s3qepWEamGk0huAT487c1FBgGDABo1akRyso2EGWPMmRCRzd608+fQVhrQ0ON5NLAtvzYiEgJEAHvd59HAFOBWVd2Qs4OqbnXvM4CPcYbQTqOqY1Q1SVWToqIKTajGGGPOkj8TyQKguYg0EZEKQH9gaq42U4Hb3Mf9gBmqqiJSHfgKeExV5+Y0FpEQEanlPg4FegMr/PgZjDHGFMJviURVM4H7cc64Wg1MVNWVIjJcRPq4zd4FIkUkBRgK5JwifD9wLvBkrtN8KwLfisgyYAmwFRjrr89gjDGmcFIeysgnJSWpzZEYEzgnTpwgLS2No0ePBjoUk4ewsDCio6MJDQ09ZbuILFTVpML29+dkuzHGAJCWlka1atVo3LgxInmdY2MCRVVJT08nLS2NJk2anNUxrESKMcbvjh49SmRkpCWREkhEiIyMLFJv0RKJMaZYWBIpuYr638YSSX5UYdFHsHpaoCMxxpgSzRJJfrKzYME7MO0vcGRfoKMxxhRBeno6CQkJJCQkULduXRo0aHDy+fHjx706xh133MHatWsLbDNq1CjGjx/vi5DPyIwZM/j111+L/X1z2GR7foJDoM9rMOZi+P4p57ExplSKjIxkyZIlADz99NNUrVqVhx9++JQ2qoqqEhSU99/X77//fqHvM2TIkKIHexZmzJhBrVq16NSpU0De33okBakXD+ffB4vGQercwtsbY0qVlJQUYmNjGTx4MImJiWzfvp1BgwaRlJRE69atGT58+Mm2Xbp0YcmSJWRmZlK9enWGDRtGfHw8559/Prt27QLgiSeeYOTIkSfbDxs2jA4dOtCyZUt+/vlnAA4dOsR1111HfHw8AwYMICkp6WSS8/TII48QExNDXFwcf/vb3wDYuXMn1157LUlJSXTo0IFff/2VDRs28M477/Dvf/+bhISEk+9TnKxHUpiLHoNVX8CXD8G9cyHEihAbUxT/+HIlq7Yd8OkxY+qH89RVrc9q31WrVvH+++8zevRoAEaMGEHNmjXJzMzk4osvpl+/fsTExJyyz/79++nWrRsjRoxg6NChvPfeewwbdtqSS6gq8+fPZ+rUqQwfPpxvvvmG119/nbp16zJ58mSWLl1KYmLiafvt3LmTr7/+mpUrVyIi/P777wA8+OCDPProo3Tq1InU1FR69+7NihUruPvuu6lVqxZ//vOfz+o7KCrrkRSmQhXo/Qqkr4fZLwc6GmOMjzVr1oz27duffP7JJ5+QmJhIYmIiq1evZtWqVaftU6lSJXr16gVAu3btSE1NzfPY11577Wlt5syZQ//+/QGIj4+ndevTE2DNmjUJCgpi4MCBTJkyhSpVqgDwww8/MHjwYBISErj66qvZt28fR44cOevP7ivWI8lHdrby3txN1A4Po0/8ZdDmBpj9ErS+Bmq3CnR4xpRaZ9tz8JecH2mA9evX8+qrrzJ//nyqV6/OzTffnOf1FRUqVDj5ODg4mMzMzDyPXbFixdPaeFNNJDQ0lOTkZL7//nsmTJjAW2+9xXfffXeyh+P5/iWB9Ujyka3Ktyt38Phny9mcfgh6/AsqVnWGuLKzAx2eMcYPDhw4QLVq1QgPD2f79u18++23Pn+PLl26MHHiRACWL1+eZ48nIyODAwcO0Lt3b1555RUWL14MwGWXXcaoUaNOtsuZW6lWrRoZGRk+j9VblkjyERIcxMj+bQkOEh74ZDHHwyKh+7Ow5VdY9EGgwzPG+EFiYiIxMTHExsYycOBAOnfu7PP3eOCBB9i6dStxcXG89NJLxMbGEhERcUqb/fv3c+WVVxIfH88ll1zCyy87w+qjRo1i7ty5xMXFERMTw9ixTs3avn37MnHiRNq2bRuQyXYr2liIb1bsYPB/FjKoa1Me79UKPuwD25bAkPkQXs/HkRpTNq1evZrzzjsv0GGUCJmZmWRmZhIWFsb69evp3r0769evJyQksDMNef038rZoo/VICtEzti43d2rEmFkb+XHdbug9ErKOw/RHAx2aMaYUOnjwIJ07dyY+Pp7rrruOt99+O+BJpKhKd/TF5IkrY0hO3cfD/13K1w9dSO1uj8L/hsOar6DVlYEOzxhTilSvXp2FCxcGOgyfsh6JF8JCg3l9QFsOHstk6KdLye70ANSOga8ehqO+PR/eGGNKG0skXmpepxpPXdWaOSl7GD33N+jzOmRshxnPBDo0Y4wJKEskZ6B/+4ZcGVePl75bx6LsZtBhEMwfA2m2+qIxpvyyRHIGRITnrm1DvYgwHvxkMfsvGAbh9WHqg5B1ItDhGWNMQFgiOUPhYaG8NqAt2/cf5fGvU9Er/g27VsLPVh3YmJLqoosuOu3iwpEjR3LfffcVuF/VqlUB2LZtG/369cv32IVdXjBy5EgOHz588vkVV1xxsn5WcUlNTeXjjz/2y7EtkZyFxEY1+Gv3Fny1bDufHmgD5/WBH5+H9A2BDs0Yk4cBAwYwYcKEU7ZNmDCBAQMGeLV//fr1mTRp0lm/f+5E8vXXX1O9evWzPt7ZsERSAg3u2owu59bi6S9XsrH93yEkDKb92VlZ0RhTovTr149p06Zx7NgxwPlR3bZtG126dOHgwYNceumlJCYm0qZNG7744ovT9k9NTSU2NhaAI0eO0L9/f+Li4rjxxhtPKZp47733nixB/9RTTwHw2muvsW3bNi6++GIuvvhiABo3bsyePXsAePnll4mNjSU2NvZkCfrU1FTOO+88Bg4cSOvWrenevXuexRn/+9//EhsbS3x8PF27dgUgKyuLRx55hPbt2xMXF8fbb78NwLBhw5g9ezYJCQm88sorPvlec9h1JGcpKEh4+YZ4er06m/umbufLS54idPpQWPIxtL0p0OEZU3JNHwY7lvv2mHXbQK8R+b4cGRlJhw4d+Oabb+jbty8TJkzgxhtvREQICwtjypQphIeHs2fPHjp16kSfPn3yXcf8rbfeonLlyixbtoxly5adUgb+2WefpWbNmmRlZXHppZeybNkyHnzwQV5++WVmzpxJrVq1TjnWwoULef/995k3bx6qSseOHenWrRs1atRg/fr1fPLJJ4wdO5YbbriByZMnc/PNN5+y//Dhw/n2229p0KDByaGyd999l4iICBYsWMCxY8fo3Lkz3bt3Z8SIEbz44otMm+b75cOtR1IEtcPDeOmGeNbsyGD4tiRo2Am++z84uDvQoRljcvEc3vIc1lJVHn/8ceLi4rjsssvYunUrO3fuzPc4s2bNOvmDHhcXR1xc3MnXJk6cSGJiIm3btmXlypV5FmT0NGfOHK655hqqVKlC1apVufbaa5k9ezYATZo0ISEhAci/VH3nzp25/fbbGTt2LFlZWQB89913fPjhhyQkJNCxY0fS09NZv369l9/S2bEeSRFd1LI2Ay9swtjZm+je5wku/OEa+PZxuG5soEMzpmQqoOfgT1dffTVDhw5l0aJFHDly5GRPYvz48ezevZuFCxcSGhpK48aN8ywd7ymv3sqmTZt48cUXWbBgATVq1OD2228v9DgF1TrMKUEPThn6vIa2Ro8ezbx58/jqq69ISEhgyZIlqCqvv/46PXr0OKXtjz/+WGAsRWE9Eh94pEcr4qIjGPLdYQ60fwCWT4SUHwIdljHGQ9WqVbnooou48847T5lk379/P7Vr1yY0NJSZM2eyefPmAo/TtWtXxo8fD8CKFStYtmwZ4JSgr1KlChEREezcuZPp06ef3Ce/Mu9du3bl888/5/Dhwxw6dIgpU6Zw4YUXev2ZNmzYQMeOHRk+fDi1atViy5Yt9OjRg7feeosTJ5xLEtatW8ehQ4f8WmreEokPVAgJ4rX+bclWGLSpGxrZHKb9BY4fCnRoxhgPAwYMYOnSpSdXKAS46aabSE5OJikpifHjx9OqVcEL1917770cPHiQuLg4XnjhBTp06AA4qx22bduW1q1bc+edd55Sgn7QoEH06tXr5GR7jsTERG6//XY6dOhAx44dufvuu2nbtq3Xn+eRRx6hTZs2xMbG0rVrV+Lj47n77ruJiYkhMTGR2NhY7rnnHjIzM4mLiyMkJIT4+HifT7ZbGXkf+mLJVh6asITn2x3gxpWD4YIHoLuVUDHGysiXfFZGvoTom9CAfu2iGbYonJ3n9odfRjlrlxhjTBlmicTH/tGnNU1qVeGmzb3IrhTpLM2blfd6zsYYUxZYIvGxKhVDeH1AW347XJG3q9wD25fAvNGBDsuYgCsPw+ilVVH/21gi8YPW9SN47IpWPL/lPH6r1RVmPgv7Cj4TxJiyLCwsjPT0dEsmJZCqkp6eTlhY2Fkfw64j8ZPbL2jM3JQ93LLuBmZUSiZ46gNwy+cQZLnblD/R0dGkpaWxe7ddrFsShYWFER0dfdb7+zWRiEhP4FUgGHhHVUfker0i8CHQDkgHblTVVBG5HBgBVACOA4+o6gx3n3bAB0Al4GvgIS2Bf+aICC/0i6fXq/t5Re7g4U2jYP7b0OneQIdmTLELDQ2lSZMmgQ7D+Inf/jwWkWBgFNALiAEGiEhMrmZ3AftU9VzgFeB5d/se4CpVbQPcBnzksc9bwCCguXvr6a/PUFQ1q1Tg1f5tefPABaysej58/xTsWhPosIwxxqf8Oc7SAUhR1Y2qehyYAPTN1aYvMM59PAm4VEREVRer6jZ3+0ogTEQqikg9IFxVf3F7IR8CV/vxMxRZp6aRPHBJC27bcyvHgivDZwMh83igwzLGGJ/xZyJpAGzxeJ7mbsuzjapmAvuByFxtrgMWq+oxt31aIccEQEQGiUiyiCQHelz2wUubc27TJvz1yF2wYxn8FJhaQ8YY4w/+TCR51WDOPZdRYBsRaY0z3HXPGRzT2ag6RlWTVDUpKirKi3D9JzhIeLV/W36u0IlvK1yOznkFfvs1oDEZY4yv+DORpAENPZ5HA9vyayMiIUAEsNd9Hg1MAW5V1Q0e7T1PLcjrmCVSHbfk/NAD/dkXWgem3APH/FNAzRhjipM/E8kCoLmINBGRCkB/YGquNlNxJtMB+gEzVFVFpDrwFfCYqs7Naayq24EMEekkTh3nW4HTlzMroS5uWZubu7VmYMYgdN9v8M1jgQ7JGGOKzG+JxJ3zuB/4FlgNTFTVlSIyXET6uM3eBSJFJAUYCgxzt98PnAs8KSJL3Ftt97V7gXeAFGAD8Eet5lLg4e4tyW7YkXe0Dyz+CNZ8FeiQjDGmSKz6bwCk7TtMn1dnMin4CZpUPIDc9wtUrV34jsYYU4ys+m8JFl2jMv/q1457Dg8m68gBmPoglIOEbowpmyyRBEjP2Lp0Pr8L/zp+I6ybDos+DHRIxhhzViyRBNBjV7RiQZ3r+ZU2ZH8zDPZuDHRIxhhzxiyRBFDFkGBe+1MSj2ffy+HMIPSzQbZ2iTGm1LFEEmBNalXhoWsv4vFjtyNpC2Cub9dSNsYYf7NEUgL0TWhA5XY3MjXrfLJnjoBtiwMdkjHGeM0SSQnx1FWtGVf9fnZrOJmTBsKJI4EOyRhjvFJoIhGRKiIS5D5uISJ9RCTU/6GVL5UqBDPi5m4My7qXkL3ryf7+qUCHZIwxXvGmRzILp4x7A+B/wB04C0sZH2tepxq9+gzgvcyeBM1/GzbMCHRIxhhTKG8SiajqYeBa4HVVvQZnoSrjB9cnRbM65i+sz27A8UmD4fDeQIdkjDEF8iqRiMj5wE04hRTB1nr3GxHhqeuS+HeVvxJ0ZA/Hpg4NdEjGGFMgbxLJn4HHgClu0cWmwEz/hlW+Va0YwoM39+O1rH5UXDOF7KUTAx2SMcbkq9BEoqo/qWof4A33+UZVfdDvkZVzsQ0iiOr5KMnZLTjx5V9gf1rhOxljTAB4c9bW+SKyCqcUPCISLyJv+j0yw80XNOPzxk9y4kQmGRMGQnZ2oEMyxpjTeDO0NRLoAaQDqOpSoKs/gzIOEeGRAb14PfQuqm3/mSNzRgU6JGOMOY1XFySq6pZcm7L8EIvJQ0SlUHrc8jA/ZLcjZOY/ULvq3RhTwnhz9tUWEbkAUHfJ3Adxh7lM8Ug8pyYfdn2ePbNvIPKdnoT2fRWJ7x/osIwx/pKdDScOwbGDcPwgHMtw7714fvJxBhw/BH9dC8H+PdHWm6MPBl4FGgBpwHfAEH8GZU538yVJPLPrXXqseYyOU+4hc/M8Qq4YASEVAx2aMcZTdrb7Q34Aju6HowfcxwfgmOfz3K953B/LALxc7K5CVedWMee+GlRv9Me27Ey/JxJbarcUUVXenrkOZvyTwSFfcqJuAqH9P3L+0RhjfCvrBBz5HY7s87jtzfXcvR3O2f67kwgKSwJBoRAWDhXDPe4jTn3umRhyJ4qc56FVIMh/JRO9XWq30EQiIuOAh1T1d/d5DeAlVb3TJ5EWg7KSSHJ8s2IH0z4dy3PBb1KpYgVC+r0DzS8PdFjGlB6H90J6CqRvcO73bYLD6R6JYR8cz8h/fwmCsOpQqYZzq1zTuQ+r7iSEwpJEaCUQKb7Pe5a8TSTe9HficpIIgKruE5G2RYrOFEnP2LpE13iQ2z5owr+OvkDL8dcjXR+Bi4ZBUHCgwzOmZDiW4SSKvRvchOEmjb0bnGSRQ4KhekOoUhuq1oWo805NDpVqQKWcpOFuqxju155AaeNNIgkSkRqqug9ARGp6uZ/xo9gGEbx5fz+GjKvHjbtf5YZZL6BpC5Dr3oEqtQIdnjHF48RRpzfhmSRyksbBHae2DY+GyGbQ+hqIPBdqNnPuqzeCkAqBib+M8CYhvAT8LCKT3OfXA8/6LyTjrboRYfxn8EUMnVid5NWf8OymcYSMvhC54UNo2D7Q4RnjW0d+hx3LYccy2L4Mti+FPetAPa5GqBLlJIdzL3OSRqSbLGo0gQqVAxd7GefVZLuItAYuBgT4n6qu8ndgvlTW5khyy85WXvp+LT/++D/eq/w6tXUP0uNZ6DCoVIzDGnOajJ1uwljiJI0dy2Bf6h+vV6sHdeOgXhxEtYKaTZ2kERYRsJDLIl/OkQCsAfbltBeRRqr6WxHiMz4UFCQ80qMVzaKqcuXk2rwa9jadpz8Kv/0KfV53zu4wpiRSdRKEZy9jxzI4uPOPNjWbQr0ESLwV6sY7yaNq7YCFbE5XaCIRkQeAp4CdOFe0C865bXH+Dc2cqWsTo2lY81IGfxjOrXzOg6smIDtXwA0fQe1WgQ7PGOfaic2/wOY5sG2JkzSO7ndek2Cnd9HsEqgX7/Q46sZaL6MU8Ob03xSgo6qmF09IvlfWh7Zy+y39MHeNW0Cd9Pm8U+VNwvQY9HkN2vQLdGimvDl20OkZp86CTbOdoSrNhuAKfwxN1Y1zEkftGAgNC3TExoMvh7a2APuLHpIpLo0iKzP5vgu4/+NKdF1Xj0m1xtBo8l3O/9A9nrWr4Y3/HD8MW351kkbqHNi2yLmyOigUopPgwoehyYUQ3d65lsKUCd4kko3AjyLyFXAsZ6Oqvuy3qEyRhYeF8t5tSfxzWmUu+eWvjKr9JT0WjHX+x75+nHPevDFFdeIIpC1wE8dsSEuG7BMQFAL1E6HzQ9C4CzTsCBWqBDpa4yfeJJLf3FsF92ZKiZDgIP7RN5Zmtaty35eh3F69Gf+3+w2C3u4KvZ6HNtfbWV3mzGQec5JF6mwneaQtgKxjzpXe9RLg/PugcVdo1MlO8ihHvK61JSJVVPWQn+Pxi/I2R5KXn9bt5v7xi2gesoOParxDlT3LnL8Sez0P9a1QgSnAsQxY/z2s/hLWf+cUJESc+Y3GF0ITN3HYpHiZ48taW+cD7wJVVbWRiMQD96jqfb4J1f8skTjW78zgznEL2HXgCO/Hr+P81FHIoT3Q9ia45O9QrU6gQzQlxeG9sPZrWD0NNsxweh1VoqDVldC8O5xzgVMqxJRpvkwk84B+wFRVbetuW6GqsV4E0ROnBH0w8I6qjsj1ekXgQ6AdzgqMN6pqqohEApOA9sAHqnq/xz4/AvWAI+6m7qq6q6A4LJH8If3gMf786RJmr9/DJY3DGFn/e8KXvAMhYdDtUeg42MpFlFcHtsGar2D1VEid61wxHtEQzrvKuTXsaLXcyhmfXpCoqlvk1LH0QldIFJFgYBRwOc46JgtEZGquq+LvAvap6rki0h94HrgROAo8CcS6t9xuUlXLDGchsmpFPryzAxMWbOHZr1Zz/raL+Fe3XvTZMQr5/klY+AH0fM75q9PmT8q+9A2wZpozbJW2wNlWqyV0+YuTPOrF278DUyh/rpDYAUhR1Y0AIjIB6At4JpK+wNPu40nAGyIi7lzMHBE517uPYc6EiDCgQyO6tojib5OW8dD3e5h47lBe7XsLteb8Az6+walV1OM5iGoR6HCNL6nCrlVO4lj9Jexc4WyvlwCXPOkkj6iWgY3RlDr+XCGxAc41KDnSgI75tVHVTBHZD0QCewo59vsikgVMBp7R8rA6lx80qF6Jj+7qwMfzf+NfX62m22+h/F+vCQzQb5Cfnoe3znfqdXX7m1NG25ROqrB1oTNktfpL2LsREGh0vvPHQqsrocY5gY7SlGIFJhJ3eOoWVb3pLI6dV3849w++N21yu0lVt4pINZxEcgvOPMupBxYZBAwCaNTIVhDMj4hwU8dz6No8imGfLePxL9YyvXl7XrhtLvWSX4Rf34Jlnzp/rSbeamPkpUV2NmxNhpWfw6ov4ECac21Hk25wwYNO8rB6VcZHClyZRVWzcIafzkYa4HnVWzSwLb82IhICRAB7C4lpq3ufAXyMM4SWV7sxqpqkqklRUVFn9QHKk4Y1K/OfuzryzNWxLNy8j8vfXs0ndR9GB/3ojJlP+zOM6eZMwpqSKTvbqV4wfRiMjIV3L4cFY516VVePhkdS4JbPIOkOSyLGp7wZ2porIm8AnwInryNR1UWF7LcAaC4iTYCtQH/gT7naTAVuA37BOTNsRkHDVG6yqa6qe0QkFOgN/ODFZzBeEBFu7nQO3VpE8eikZTz22XKmt4hixDWTqb91Onz3d/jgCmdhoMv/aVfHlwTZ2U5JklVfwKqpkLHNqWN17mVw6VPQsqdd32H8zpvTf2fmsVlV9ZJCDy5yBTAS5/Tf91T1WREZDiSr6lQRCQM+Atri9ET6e0zOpwLhOFfT/w50BzYDs4BQ95g/AEPdnlO+7PTfM5edrYyft5nnpq8hWIQnep/HDfGRyNzXYO5Ip1HnhyDpTqhWN7DBljfZWU7PY9XnTvI4uAOCK0LzyyGmL7To6awNbkwR+ew6krLAEsnZ+y39MI9OXsqvG/fSrUUUI65rQz3dA98/CSunAAINOzhn+7TqDTWbBDrksik7Czb/7CSP1V8663WEhLnJ42po0QMqVgt0lKaM8eUFiXWAfwH1VbWXiMQA56vqu74J1f8skRRNdrby0a+bGTF9DSHBwpO9Y7i+XTSyZ53zF/Hqqc66EgB128B5fdzTSFvZNQhFcWiPU9dq/bdO8ji0G0IqOcmj9dXQvIfVszJ+5ctEMh14H/g/VY135ykWq2ob34Tqf5ZIfGNz+iEembSM+Zv2cnHLKJ67No66Ee76EftSnXIaq7+ELfMAhcjmf1wVXb+tJZWCZJ1wrulIS3YuDExb4J6mC4RWdi4QbX21c29VdE0x8WUiWaCq7UVksUeJlCWqmuCjWP3OEonvZGcr435J5flv1hAaHMQdnZtw6/nnUKuqxxonGTv+KLWxabZTaiM8+o+k0qiTnUZ8YPsfCSMtGbYthky36k/VOs56HTm3+m2hQuXAxmvKJV8mkh+B64DvVTVRRDoBz6tqN59EWgwskfhe6p5DPPPVKn5YvYsKIUFclxjN3Rc2oVlUrqGWw3th3TdOTyXlf07xv8q1nOsYYvo4JcfLem2vE0edoT/PxLHfvVY3uIJThiS6vbPwU3QHiIi23pspEXyZSBKB13FqXq0AooB+qrrMF4EWB0sk/pOy6yDvztnE5EVpHM/M5rLzajPwwqZ0aFITyf1jeOwgpLjlyNd965QjrxgBzS9z5lZqtYRaLaBGYwj2qgxcyaLqzGukp8DeDbBjhZM4diyDrONOm4hGbsJwexv14mzFSlNiFTmRiMj1qvpf9zqQLUBLnCvR16rqCZ9G62eWSPxvz8FjfPjLZj76JZV9h08QFx3BwAub0iu2LiHBeVz3euIobPrJmazf8D/I2P7Ha0GhULOpU+erlueteck4M+noASdRpG9wkkZ6ivt4AxzzWJU6pBI0SDw1cdip0qYU8UUiWeQOZS1S1USfR1iMLJEUnyPHs5i8KI1352xi055DNKheibu6NOGG9g2pWrGAXsbR/bBnPexZ59x2u/d7NzpzLDmq1c8jwbRwfqB9ORx04ijs23R6okhPgUOeqxaIU2o9spl7O9e9NXN6H6WxZ2WMyxeJ5Aeci/4SgNm5X1fVPkUNsrhYIil+2dnKD6t3Mnb2Rhak7qNaWAg3dTyH2y9o/MeZXt7IPO7v1b44AAAc6klEQVT8oOdOMHvWw/GMP9pVqOYUHpSc3o/Hv+tT/ol7bte8tx876M5heGyrUvv0ZFGzmXPdTGgl7z+PMaWILxJJBSAR58rzu3O/rqo/FTXI4mKJJLAW/7aPd2ZvYvqK7QQHCVfF12fghU05r14Rrr5WdYbDcpLK7rWwP+3UNqf0UCSf7Xm0D6nkJIyazf5IHlZmxJRDvkgkH6nqLSLyqKq+4PMIi5ElkpLht/TDvDd3ExOTt3D4eBYXNq/FwAubcmHzWqdPzBtjAs4XiWQV0AunsOJF5Cr5rqoFVuktSSyRlCy/Hz7O+Hm/8cHPqezOOEarutUYeGFTroqvT4WQAgtSG2OKkS8SyYPAvUBTnOq9nolEVbWpLwItDpZISqZjmVlMXbKNsbM3sm7nQeqGh3FH58YM6NiI8LDQQIdnTLnny+tI3lLVe30WWQBYIinZVJWf1u1mzKyN/LwhnaoVQxjQoSF3dG5C/eo2kW1MoPiiRxKuqgdEpGZer9vQlvGHFVv3M2bWRr5avh0BesfVY2DXprSub5PdxhQ3XySSaaraW0Q25fGyDW0Zv0rbd5j356YyYf5vHDqeRZdzazGwa1O62sS8McXG1iPxYImk9Np/+AQfz/+N9+duYpdNzBtTrHySSNyS8b2AVu6mVcC3qprpkyiLiSWS0i/3xHyd8Irc0bkJf7KJeWP8xhdDW/WBmcB2YDHOWVttgbrAxaq6zXfh+pclkrJDVflx3W7GekzM92/fkDu72MS8Mb7mi0TyAbBEVUfm2v4g0E5Vb/NFoMXBEknZtDxtP2Nn/zEx3ye+Pn/t0ZIGllCM8QlfJJI1qtoqn9fWqmrLIsZYbCyRlG1p+w7z3pxUxs/bTJAIQy5uxt0XNiUstJwvnmVMEXmbSAqarTxSwGuHzzwkY/wjukZl/n5VDD8M7Ua3FlG8+N06eoycxYw1OwMdmjHlQkE1riNE5No8tgtQhGp7xvhHw5qVGX1LO2av383TU1dy5wfJXNKqNn/vHUPjWrbOuTH+UtDQ1vsF7aiqd/glIj+woa3y53hmNuN+TmXkD+s4kaUM7NqEIRefS+UKtj6IMd6y60g8WCIpv3YdOMpz09cwZfFW6kWE8cSVMVzRpq5d1GiMF3wxR2JMqVc7PIxXbkzgv4PPp0blCgz5eBF/GjuPdTszCt/ZGOMVSySmXGjfuCZfPtCFf14dy6rtB+j16mz+OW0VB46eCHRoxpR6+SYSEbnevW9SfOEY4z/BQcItnc5h5sMXcUNSQ96bu4lLXvyJSQvTyM4u+0O8xvhLQT2Sx9z7ycURiDHFpWaVCjx3bRumDulCw5qVePi/S+k3+mdWbN0f6NCMKZUKOmvre5zTgxOA2blfV9U+/g3Nd2yy3eQnO1v5bPFWRkxfTfqh4wzo0IhHurekRpUKgQ7NmIDzdrK9oHMhrwQSgY+Al3wVmDElSVCQ0K9dNN1b1+HVH9bzwc+pfLNiB2/8qS0XNKsV6PCMKRW8WSExSlV3i0g1nHVIDhZPaL5jPRLjrbU7Mhjy8SI27TnEE1eex+0XNLZThU255cvTf+uIyGJgBbBKRBaKSGyRIzSmBGpZtxpT7ruAi1vW5h9fruKRScs4eiIr0GEZU6J5k0jGAENV9RxVbQT81d1WKBHpKSJrRSRFRIbl8XpFEfnUfX2eiDR2t0eKyEwROSgib+Tap52ILHf3eU3sz0XjY9XCQhlzSzseurQ5kxamceOYX9mx/2igwzKmxPImkVRR1Zk5T1T1R6DQwkUiEgyMwlkYKwYYICIxuZrdBexT1XOBV4Dn3e1HgSeBh/M49FvAIKC5e+vpxWcw5owEBQl/ubwFb9/SjpSdGfR+fQ4LN+8NdFjGlEjeJJKNIvKkiDR2b08Aea3jnlsHIEVVN6rqcWAC0DdXm77AOPfxJOBSERFVPaSqc3ASykkiUg8IV9Vf1Jnc+RC42otYjDkrPVrXZcqQzlStGEz/Mb/yyfzfAh2SMSWON4nkTiAK+My91QK8KdjYANji8TzN3ZZnG3f53v1AZCHHTCvkmMb4VIs61fhiSBcuaFaLxz5bzhOfL+d4ZnagwzKmxCi0FKqq7gMePItj5zV3kfsUMW/anFV7ERmEMwRGo0aNCjikMYWLqBzKe7e359/frmX0TxtYuyODN29qR1S1ioEOzZiA82etrTSgocfzaCD3Ou8n24hICBABFDQQneYep6BjAqCqY1Q1SVWToqKizjB0Y04XHCQM69WK1wa0ZfnW/fR5Yw7L0n4PdFjGBJw/E8kCoLmINBGRCkB/YGquNlOBnLXf+wEztIALW1R1O5AhIp3cs7VuBb7wfejG5K9PfH0m33sBQSL0G/0Lny1KK3wnY8owvyUSd87jfuBbYDUwUVVXishwEckpr/IuECkiKcBQ4OQpwiKSCrwM3C4iaR5nfN0LvAOkABuA6f76DMbkp3X9CKbe35nERtUZOnEp/5y2iswsmzcx5ZNXV7YDA4HGeMypqOqdfo3Mh+zKduMvJ7Ky+dfXq3l/biqdz43kjQGJVqfLlBm+vLL9C5y5ix+ArzxuxpR7ocFBPHVVa/7dL44Fqfu46o05rN5+INBhGVOsvFnAurKq/s3vkRhTil2f1JDmdaox+KOFXPvmz7x4fTxXxtULdFjGFAtveiTTROQKv0diTCmX0LA6Ux/oTEz9cIZ8vIiXv1tLYUPHxpQF3iSSh3CSyREROSAiGSJifXdj8lC7WhifDOzEDUnRvDYjhddnpAQ6JGP8zpsLEqsVRyDGlBUVQoIYcW0c2Qovf7+O8LAQbu9sK1abssubORJEpAZOgcSwnG2qOstfQRlT2gUFCSOubcOBIyd4+stVRFQO5Zq20YXvaEwpVOjQlojcDczCuR7kH+790/4Ny5jSLyQ4iNcGtOWCZpE8/N9l/LBqZ6BDMsYvvJ0jaQ9sVtWLgbbAbr9GZUwZERYazJhbk4itH859Hy/ilw3pgQ7JGJ/zJpEcVdWj4CxEpaprgJb+DcuYsqNqxRA+uKMD59SszMAPk60+lylzvEkkaSJSHfgc+F5EviCfQonGmLzVqFKBj+7qSPXKodz23nxSdmUEOiRjfKbQRKKq16jq76r6NM6qhe9ii0kZc8bqRoTxn7s6EhwUxC3vzidt3+FAh2SMT+SbSEQk3L2vmXMDlgNzgKrFFJ8xZUrjWlX46K4OHDqWyS3vzmd3xrFAh2RMkRXUI/nYvV8IJLv3Cz2eG2POwnn1wnn/jvZs33+E296bz/4jJwIdkjFFkm8iUdXe7n0TVW3q3ufcmhZfiMaUPe3OqcnbtySxflcGd49bwJHjWYEOyZizVtDQVmJBt+IM0piyqFuLKF65MYHkzfu4d/xCWwfelFoFXdn+knsfBiQBS3HWTI8D5gFd/BuaMWVf77j6ZBzN5LHPljN04hJe7d+W4CAJdFjGnJF8E4l78SEiMgEYpKrL3eexwMPFE54xZd+ADo3Yf+QEI6avIaJSKM9cHYuzkrQxpYM3tbZa5SQRAFVdISIJfozJmHJncLdm/H74BKN/2kBEpVAe7dkq0CEZ4zVvEslqEXkH+A+gwM04a7AbY3zobz1bsv/ICd780Ukm93RrFuiQjPGKN4nkDuBenJpb4BRwfMtvERlTTokIz1wdy4GjJ3jOHebq36FRoMMyplDerEdyVERGA1+r6tpiiMmYcis4SHjlhgQOHs3k8SnLCa8UyhVtbMleU7J5U0a+D7AE+MZ9niAiU/0dmDHlVYWQIN66OZG2jWrw0ITFzFpnxbZNyeZN0cangA7A7wCqugRo7MeYjCn3KlcI4b3b2tMsqir3jV9Eyq6DgQ7JmHx5k0gyVXW/3yMxxpwionIo793enoohQQz6MJkDR62UiimZvEkkK0TkT0CwiDQXkdeBn/0clzEGqF+9Em/elMhvew/z5wlLyM7WQIdkzGm8SSQPAK2BY8AnwAHgz/4Myhjzh45NI/n7VTHMWLOLV35YF+hwjDmNN2dtHQb+z70ZYwLglk7nsGLrfl6fkULr+uH0jLUzuUzJkW8iKezMLFXt4/twjDF5ERGG941l3c6DDJ24lCa1qtKybrVAh2UMUHCP5HxgC85w1jycgo3GmAAJCw1m9M3tuOqNOQz6KJmpQ7oQUTk00GEZU+AcSV3gcSAWeBW4HNijqj+p6k/FEZwx5lR1I8IYfXMi234/wgMTFpNlk++mBChoYassVf1GVW8DOgEpwI8i8kCxRWeMOU27c2ryjz6xzFq3m39/a8UmTOAVONkuIhWBK4EBOBchvgZ85v+wjDEF+VPHRqzYtp/RP20gtkE4vePqBzokU44VNNk+DmdYazrwD1VdUWxRGWMK9fRVrVm3I4NH/ruMprWqElM/PNAhmXKqoDmSW4AWOFV/fxaRA+4tQ0QOeHNwEekpImtFJEVEhuXxekUR+dR9fZ6INPZ47TF3+1oR6eGxPVVElovIEhFJ9vaDGlPWVAgJ4s2bEwmvFMKgj5LZd+h4oEMy5VRBcyRBqlrNvYV73KqpaqF/+ohIMDAK6AXEAANEJCZXs7uAfap6LvAK8Ly7bwzQH+dCyJ7Am+7xclysqgmqmnQGn9WYMqd2tTBG39yOXQeOcf8ni8jMsnXfTfHz5sr2s9UBSFHVjap6HJgA9M3Vpi8wzn08CbhUnDVG+wITVPWYqm7Cmejv4MdYjSm12jaqwTPXxDI3JZ3nv1kT6HBMOeTPRNIA5zqUHGnutjzbqGomsB+ILGRfBb4TkYUiMii/NxeRQSKSLCLJu3dbGW5Ttt2Q1JDbzj+HsbM38fnirYEOx5Qz/kwkeV3AmPuk9/zaFLRvZ1VNxBkyGyIiXfN6c1Udo6pJqpoUFRXlbczGlFpP9I6hQ5Oa/G3yMlZstYLdpvj4M5GkAQ09nkcD2/JrIyIhQASwt6B9VTXnfhcwBRvyMgaA0OAg3rwpkcgqFRj0YTJ7Dh4LdEimnPBnIlkANBeRJiJSAWfyPHf9rqnAbe7jfsAMVVV3e3/3rK4mQHNgvohUEZFqACJSBegO2GnJxrhqVa3I27ckkX7oOEPGL+KETb6bYuC3ROLOedwPfAusBiaq6koRGe4u3wvwLhApIinAUGCYu+9KYCKwCmeJ3yGqmgXUAeaIyFJgPvCVqn7jr89gTGnUJjqCEde1Yd6mvTz71epAh2PKAXE6AGVbUlKSJifbJSemfPnntFW8O2cT/+4Xx/VJDQvfwZhcRGShN5dZ+HNoyxgTQI/1asUFzSL5v89XsGTL74EOx5RhlkiMKaNCgoN440+J1K5WkcEfLWRXxtFAh2TKKEskxpRhNatUYMwtSfx+5Dj3/WcRxzKzAh2SKYMskRhTxsXUD+fF6+NJ3ryPJz9fQXmYFzXFq9A1240xpV/vuPqs3ZHB6zNSOK9eOHd0bhLokEwZYj0SY8qJv1zWgu4xdfjntFXMXm9lg4zvWCIxppwIChJeuTGBFnWqMWT8IjbtORTokEwZYYnEmHKkSsUQxt6aRHCQcPe4BRw4eiLQIZkywBKJMeVMw5qVeevmdmxOP8xDnywmK9sm303RWCIxphzq1DSSf/Rtzcy1u3nB1jAxRWRnbRlTTt3U8RzWbM/g7VkbaVm3GtcmRgc6JFNKWY/EmHLs71fF0KlpTYZ9tpzFv+0LdDimlLJEYkw55qxh0o464RW556OF7NhvZVTMmbNEYkw5V7NKBd65tT2HjmUy6KNkjp6wMirmzFgiMcbQsm41RvZvy/Kt+/nb5GVWRsWcEUskxhgALo+pw8PdW/LFkm2M/mljoMMxpYidtWWMOem+i5qxZkcGL3y7hua1q3JZTJ1Ah2RKAeuRGGNOEhFeuC6O2PoRPDRhMet2ZgQ6JFMKWCIxxpyiUoVgxtzajkoVQrh7XDL7Dh0PdEimhLNEYow5Tb2ISrx9Szt27D/KkI8XcSIrO9AhmRLMEokxJk/tzqnBv65tw88b0nlm2qpAh2NKMJtsN8bkq1+7aNbuOMDY2ZtoWTecP3VsFOiQTAlkPRJjTIGG9TqPbi2i+PsXK5i3MT3Q4ZgSyBKJMaZAwUHCawPa0iiyMveOX8SWvYcDHZIpYSyRGGMKFVEplHduTSIzK5ver89hYvIWu/rdnGSJxBjjlaZRVZkypDMt61Tj0UnLuOmdeWxOt+V6jSUSY8wZaBZVlQmDOvHsNbEsT9tP91dmMfqnDWTa6cHlmiUSY8wZCQoSbup4Dt8P7UbXFlGMmL6GvqPmsmLr/kCHZgLEEokx5qzUjQhjzC3teOumRHZlHKPvqLk89/Vqjhy3MvTljSUSY8xZExF6tanHD0O7cUNSNG/P2kiPkbOYs35PoEMzxcgSiTGmyCIqhfLctXF8MrATwUHCze/O4+H/LrU6XeWEJRJjjM+c3yyS6Q9dyJCLm/H54q1c/spPTF26zU4VLuP8mkhEpKeIrBWRFBEZlsfrFUXkU/f1eSLS2OO1x9zta0Wkh7fHNMYEVlhoMI/0aMXU+7vQoHolHvxkMXeNS2br70cCHZrxE78lEhEJBkYBvYAYYICIxORqdhewT1XPBV4Bnnf3jQH6A62BnsCbIhLs5TGNMSVATP1wPruvM0/2juGXDel0f/knxv2cSla29U7KGn/2SDoAKaq6UVWPAxOAvrna9AXGuY8nAZeKiLjbJ6jqMVXdBKS4x/PmmMaYEiI4SLirSxO++0tX2jWuyVNTV3L96J9twawyxp/VfxsAWzyepwEd82ujqpkish+IdLf/mmvfBu7jwo5pjClhGtaszLg72vP5kq0M/3IVV742m8aRVQIdVrkw7cEuVAwJ9ut7+DORSB7bcvdp82uT3/a8elB59pNFZBAwCKBRIyt9bUygiQjXtI2ma/Mo3piZws4DRwMdUrkgef6c+pY/E0ka0NDjeTSwLZ82aSISAkQAewvZt7BjAqCqY4AxAElJSTYoa0wJEVm1Ik9d1TrQYRgf8uccyQKguYg0EZEKOJPnU3O1mQrc5j7uB8xQ5zzBqUB/96yuJkBzYL6XxzTGGFOM/NYjcec87ge+BYKB91R1pYgMB5JVdSrwLvCRiKTg9ET6u/uuFJGJwCogExiiqlkAeR3TX5/BGGNM4aQ8XCiUlJSkycnJgQ7DGGNKFRFZqKpJhbWzK9uNMcYUiSUSY4wxRWKJxBhjTJFYIjHGGFMklkiMMcYUSbk4a0tEdgOb83ipFlDaVuApbTGXtnjBYi4uFrP/FTXec1Q1qrBG5SKR5EdEkr05ta0kKW0xl7Z4wWIuLhaz/xVXvDa0ZYwxpkgskRhjjCmS8p5IxgQ6gLNQ2mIubfGCxVxcLGb/K5Z4y/UciTHGmKIr7z0SY4wxRVQuE4mI9BSRtSKSIiLDAh2PN0QkVUSWi8gSESmRFShF5D0R2SUiKzy21RSR70VkvXtfI5Ax5pZPzE+LyFb3u14iIlcEMkZPItJQRGaKyGoRWSkiD7nbS+z3XEDMJfl7DhOR+SKy1I35H+72JiIyz/2eP3WXsygRCoj5AxHZ5PE9J/j8vcvb0JaIBAPrgMtxFtBaAAxQ1VUBDawQIpIKJKlqiT2HXUS6AgeBD1U11t32ArBXVUe4SbuGqv4tkHF6yifmp4GDqvpiIGPLi4jUA+qp6iIRqQYsBK4GbqeEfs8FxHwDJfd7FqCKqh4UkVBgDvAQMBT4TFUniMhoYKmqvhXIWHMUEPNgYJqqTvLXe5fHHkkHIEVVN6rqcWAC0DfAMZUJqjoLZ10ZT32Bce7jcTg/ICVGPjGXWKq6XVUXuY8zgNVAA0rw91xAzCWWOg66T0PdmwKXADk/yCXte84vZr8rj4mkAbDF43kaJfwftUuB70RkobsefWlRR1W3g/ODAtQOcDzeul9ElrlDXyVmmMiTiDQG2gLzKCXfc66YoQR/zyISLCJLgF3A98AG4HdVzXSblLjfjtwxq2rO9/ys+z2/IiIVff2+5TGRSB7bSsP4XmdVTQR6AUPcIRnjH28BzYAEYDvwUmDDOZ2IVAUmA39W1QOBjscbecRcor9nVc1S1QQgGmck47y8mhVvVAXLHbOIxAKPAa2A9kBNwOdDnuUxkaQBDT2eRwPbAhSL11R1m3u/C5iC8w+7NNjpjpHnjJXvCnA8hVLVne7/kNnAWErYd+2Of08GxqvqZ+7mEv095xVzSf+ec6jq78CPQCeguojkLFFeYn87PGLu6Q4tqqoeA97HD99zeUwkC4Dm7tkXFXDWiZ8a4JgKJCJV3ElKRKQK0B1YUfBeJcZU4Db38W3AFwGMxSs5P8iuayhB37U7ofousFpVX/Z4qcR+z/nFXMK/5ygRqe4+rgRchjO3MxPo5zYrad9zXjGv8fgDQ3DmdHz+PZe7s7YA3NMMRwLBwHuq+myAQyqQiDTF6YUAhAAfl8SYReQT4CKciqM7gaeAz4GJQCPgN+B6VS0xk9v5xHwRznCLAqnAPTnzD4EmIl2A2cByINvd/DjOnEOJ/J4LiHkAJfd7jsOZTA/G+YN7oqoOd/9fnIAzRLQYuNn9Sz/gCoh5BhCFM6y/BBjsMSnvm/cuj4nEGGOM75THoS1jjDE+ZInEGGNMkVgiMcYYUySWSIwxxhSJJRJjjDFFYonElBoioiLykcfzEBHZLSLTzvJ4fSRA1Z9FZIpbiTVFRPZ7VGa9wMfvc5l7/MUisk5EfipKlV0RaSoi/T2e3y0iI30TrSmtQgpvYkyJcQiIFZFKqnoEp4Lz1rM9mKpOJUAXo6rqNQAichHwsKr2zqudiIR41HY6WzNV9Wr3eInAFBG5VVV/OotjNcW5iHdCEWMyZYj1SExpMx240n08APgk5wUR6SAiP7t/ff8sIi3d7UNF5D33cRsRWSEilUXkdhF5w93+gYi8Jc66GRtFpJtbSHC1iHzg8R4HPR73y3nN2/29ISJpIvKkiMwFrhGR5iLyrVuwc5aItHDb1RGRz0QkWZx1KDoVdmy3Cu+zwP0FHUNEnhGRce7nWS8id7qHGAFc7PaeHnS3RbvxrReR587ks5oyQlXtZrdSccNZNyQOp4x3GM5VuhfhrLUAEA6EuI8vAya7j4OAWThlOJJxCmCCs4bHG+7jD3D+yhackuwHgDbuvguBhJwYPOLpB3xwJvvn8ZlOxu+xLQ0Y6vF8JtDMfdwZ+M59/CnQyX3cGFiRx/EvAz7PtS0JWF7QMYBngEXu91zbjalO7uMBdwPrgWpAJZzK2vUD/W/FbsV7s6EtU6qo6jJxSpEPAL7O9XIEME5EmuOU3Qh198kWkduBZcDbqjo3n8N/qaoqIsuBnaq6HEBEVuL8yC4pJLyi7u/pU3ff6jjFAic7pZKAP4akLwNaemyv4THsVxDPCth5HsN9/LmqHgWOisgsnOqxR/M43g/qrDOCiKzBKdNSIosZGv+wRGJKo6nAizh/zUd6bP8nznzANW6y+dHjteY4PZr6BRw3p2ZStsfjnOc5/6941hQKO4v9vXXIvRdgjzqlwXMToIM6C7SdibY4BQjzPYabWHLXT8qvnpLnZ83CflfKHZsjMaXRe8DwnL/4PUTwx+T77TkbRSQCeBXoCkSKSD/O3k4ROU9EgnCGyvxKVfcB20UkZ3I+SETi3Zd/AIbktBUv1uJ22zwOjPLiGFeLSEURqQVciDMsmIEzjGXMSZZITKmjqmmq+moeL70APOdOUgd7bH8FeFNV1wF3ASNE5GxXEBwGTANm4CzGVBz6A4NFZCmwEsg5w2sI0Fmcle9WAQPz2f9i9wSEtcBrwH36xxlbBR1jAc7JDb8AT6nqTpyKt8EistRjst2Uc1b91xhzGhF5BmdIza4RMYWyHokxxpgisR6JMcaYIrEeiTHGmCKxRGKMMaZILJEYY4wpEkskxhhjisQSiTHGmCKxRGKMMaZI/h9qokE2s9wVTgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from matplotlib import pyplot as plt\n",
    "from sklearn.tree import DecisionTreeRegressor\n",
    "\n",
    "# paste your get_features_targets function here\n",
    "def get_features_targets(data):\n",
    "  features = np.zeros((data.shape[0], 4))\n",
    "  features[:, 0] = data['u'] - data['g']\n",
    "  features[:, 1] = data['g'] - data['r']\n",
    "  features[:, 2] = data['r'] - data['i']\n",
    "  features[:, 3] = data['i'] - data['z']\n",
    "  targets = data['redshift']\n",
    "  return features, targets\n",
    "\n",
    "# paste your median_diff function here\n",
    "def median_diff(predicted, actual):\n",
    "  return np.median(np.abs(predicted - actual))\n",
    "\n",
    "# Complete the following function\n",
    "def accuracy_by_treedepth(features, targets, depths):\n",
    "  # split the data into testing and training sets\n",
    "  split = features.shape[0]//2\n",
    "  train_features, test_features = features[:split], features[split:]\n",
    "  train_targets, test_targets = targets[:split], targets[split:]\n",
    "\n",
    "  # Initialise arrays or lists to store the accuracies for the below loop\n",
    "  train_diffs = []\n",
    "  test_diffs = []\n",
    "\n",
    "  # Loop through depths\n",
    "  for depth in depths:\n",
    "    # initialize model with the maximum depth. \n",
    "    dtr = DecisionTreeRegressor(max_depth=depth)\n",
    "\n",
    "    # train the model using the training set\n",
    "    dtr.fit(train_features, train_targets)\n",
    "\n",
    "    # Get the predictions for the training set and calculate their med_diff\n",
    "    predictions = dtr.predict(train_features)\n",
    "    train_diffs.append(median_diff(train_targets, predictions))\n",
    "\n",
    "    # Get the predictions for the testing set and calculate their med_diff\n",
    "    predictions = dtr.predict(test_features)\n",
    "    test_diffs.append(median_diff(test_targets, predictions))\n",
    "        \n",
    "  # Return the accuracies for the training and testing sets\n",
    "  return train_diffs, test_diffs    \n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    data = pd.read_csv('opticaldatafinals_SDSS.csv')\n",
    "    features, targets = get_features_targets(data)\n",
    "\n",
    "  # Generate several depths to test\n",
    "    tree_depths = [i for i in range(1, 36, 2)]\n",
    "\n",
    "  # Call the function\n",
    "    train_med_diffs, test_med_diffs = accuracy_by_treedepth(features, targets, tree_depths)\n",
    "    print(\"Depth with lowest median difference : {}\".format(tree_depths[test_med_diffs.index(min(test_med_diffs))]))\n",
    "\n",
    "  # Plot the results\n",
    "    train_plot = plt.plot(tree_depths, train_med_diffs, label='Training set')\n",
    "    test_plot = plt.plot(tree_depths, test_med_diffs, label='Validation set')\n",
    "    plt.xlabel(\"Maximum Tree Depth\")\n",
    "    plt.ylabel(\"Median of Differences\")\n",
    "    plt.legend()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1><center> KFold Cross Validation</center></h1>\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Differences: 0.014, 0.014, 0.014, 0.014, 0.014, 0.014, 0.014, 0.014, 0.014, 0.014\n",
      "Mean difference: 0.014\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from sklearn.model_selection import KFold\n",
    "from sklearn.tree import DecisionTreeRegressor\n",
    "\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from sklearn.tree import DecisionTreeRegressor\n",
    "\n",
    "def get_features_targets(data):\n",
    "  features = np.zeros((data.shape[0], 4)) #n lines, 4 columns\n",
    "  features[:,0] = data['u'] - data['g']\n",
    "  features[:,1] = data['g'] - data['r']\n",
    "  features[:,2] = data['r'] - data['i']\n",
    "  features[:,3] = data['i'] - data['z']\n",
    "  targets = data['redshift']\n",
    "  return (features, targets)\n",
    "\n",
    "def median_diff(predicted, actual):\n",
    "  diff = np.median(np.absolute(predicted - actual))\n",
    "  return diff\n",
    "\n",
    "# complete this function\n",
    "def cross_validate_model(model, features, targets, k):\n",
    "  kf = KFold(n_splits=k, shuffle=True)\n",
    "\n",
    "  # initialise a list to collect median_diffs for each iteration of the loop below\n",
    "  mediandiffs = []\n",
    "  \n",
    "  for train_indices, test_indices in kf.split(features):\n",
    "    train_features, test_features = features[train_indices], features[test_indices]\n",
    "    train_targets, test_targets = targets[train_indices], targets[test_indices]\n",
    "    \n",
    "    # fit the model for the current set\n",
    "    model.fit(train_features, train_targets)\n",
    "    \n",
    "    # predict using the model\n",
    "    predictions = model.predict(test_features)\n",
    "\n",
    "    # calculate the median_diff from predicted values and append to results array\n",
    "    mediandiffs.append(median_diff(test_targets, predictions))\n",
    " \n",
    "  # return the list with your median difference values\n",
    "  return mediandiffs\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "  data = pd.read_csv('opticaldatafinals_SDSS.csv')\n",
    "  features, targets = get_features_targets(data)\n",
    "\n",
    "  # initialize model with a maximum depth of 19\n",
    "  dtr = DecisionTreeRegressor(max_depth=19)\n",
    "\n",
    "  # call your cross validation function\n",
    "  diffs = cross_validate_model(dtr, features, targets, 10)\n",
    "\n",
    "  # Print the values\n",
    "  print('Differences: {}'.format(', '.join(['{:.3f}'.format(val) for val in diffs])))\n",
    "  print('Mean difference: {:.3f}'.format(np.mean(diffs)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " <h2><center>KFold Cross Validated Predictions</h2></center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Median difference: 0.014\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEKCAYAAAAB0GKPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzsnXd4XNWZ8H9nRpbciyxhyZJtWVZxrxIYGwPBBkwS2m72SwgugAuQBEiyu/QEkuxuNt9uCDbNDYxtWMIXNja92bjbAsm9YFULY8tFllyxsayZ9/tj5g5Xoyl3pBnNSD6/57mPNLed915pznvO244SETQajUajCYYt2gJoNBqNpnWgFYZGo9FoLKEVhkaj0WgsoRWGRqPRaCyhFYZGo9FoLKEVhkaj0WgsoRWGRqPRaCyhFYZGo9FoLKEVhkaj0WgsERdtAcJJUlKSJKWmU3n8GwRI6dqe5C4J0RZLo9FoYpYtW7YcF5FkK+e2KYWRkZHB518UsqGkmr1HTjPjqkzi4/QkSqPRaPyhlPrK6rltrjctqKjhqpxk7r82C7tNsbHsOA6nrpel0Wg0zaVNKYyzF+qZvbSIgooawKU8zJ81Go1G03TalMJAYN6UMeRnJLKx7Dj5GYksmJbH2Mye0ZZMo9FoWj1tSmF8VXsOm01RWFnL7KVFFFbWMj4rCbtNRVs0jUajafW0Kad3v54dPbMJPbPQaDSa8NKmZhidE+Kw2xR2m2o1MwuHU5rsmG/OtRqNRhMqbUphtEaa45jXTn2NRtOSaIURZkId9Y/N7Nlk85mVa/UsRKPRhAutMMJMqKP+5pjPrFyrZyEajSZcKJG2M/LMy8uToqKiqMrgcAoFFTWMzewZEz6UWJNHo9HEFkqpLSKSZ+XcNhUlFQsYo/5YIdbk0Wg0rRdtktJoNBqNJbTC0Gg0Go0ltMK4RNHRUxqNJlS0wrhE0dFTGo0mVC4JhaGzqRvTnPwPjUZzaXJJKAydTd2Y1lQ+RaPRxAaXhMLIz0jkwYnZ5GckhnxtuEfibXXGotFo2j6XhMIorKxl7qpSCitrQ7423CPxtjpj0Wg0bZ+IZXorpV4BfggcE5GhPo7/K3Cn+2McMAhIFpFapVQlcAZwAPVWsxD9ZXrHUrZzLMmi0Wg0oWR6R1JhXA2cBZb6Uhhe594M/EpErnN/rgTyROR4KG3GQmkQjUajaU2EojAiZpISkXWAVRvQHcAbkZJFo9FoNM0n6j4MpVRHYDLwv6bdAnyilNqilJod5PrZSqkipVRRdXV1JEXVaDSaS5qoKwzgZmCjiJhnI+NFZDRwE/Bzt3nLJyKyQETyRCQvOTk50rJqNBrNJUssKIyf4GWOEpEq989jwHLg8ijIpdFoNBoTUVUYSqluwDXA26Z9nZRSXYzfgRuA3dGRUKPRaDQGEVsPQyn1BnAtkKSUOgg8BbQDEJF57tNuBz4RkW9Ml/YCliulDPn+R0Q+ipScGo1Go7FGxBSGiNxh4ZxXgVe99lUAIyIjlUaj0WiaSiz4MDQajUbTCtAKQ6PRaDSW0ApDo9FoNJbQCkOj0Wg0ltAKQ6PRaDSWaFMK4+yF+kt6nQm91oZGo4kkbUphfFVz7pJeZ0KvtaHRaCJJxMqbR4OBw0bKnh3bLtl1JvRaGxqNJlRCKW8escS9aNA5Ie6S7iiN1QE1Go0mErQpk5RGo9FoIodWGBqNRqOxhFYYGo1Go7GEVhgajUajsYRWGBqNRqOxhFYYmjaHTmDUaCKDVhiaNodOYNRoIoNWGJo2x9jMniyYlsfYzJ7RFkWjaVO0qcQ9jQZ0AqNGEykiNsNQSr2ilDqmlNrt5/i1SqlTSqnt7u23pmOTlVLFSqkypdSjkZJRo9FoNNaJpEnqVWBykHPWi8hI9/Z7AKWUHXgBuAkYDNyhlBocQTk1Go1GY4GIKQwRWQfUNuHSy4EyEakQkTrgr8CtVi40ypvrKBmNRqMJP9F2el+plNqhlPpQKTXEvS8N+Np0zkH3vqAY5c11lIxGo9GEn2g6vbcC/UTkrFLq+8AKIBvwVW7W71RBKTUbmA3Qq3e6JzJGR8loNBpNeInaDENETovIWffvHwDtlFJJuGYUfUynpgNVAe6zQETyRCQvPbUXdpvyRMmEo9S5Nm9pNBqNi6gpDKVUilJKuX+/3C1LDVAIZCul+iul4oGfAO9ES05t3tJoNBoXETNJKaXeAK4FkpRSB4GngHYAIjIP+BFwv1KqHjgP/ERcy//VK6V+AXwM2IFXRGRPpOQMhk4C02g0GhdtaonWfrnDpHTPDuLjou3Lb4ivpVP1cqrRQ797jeY7QlmiNbZ61mZy5PS3LN64P9piNMKXWUubuqKHfvcaTdPQM4wWQM8wYgv97jWa7whlhtGmakkld0mIOWUBvmsb6XpH0UO/e42macRe7xpGdEhs9NDvXqNpe7QphWGUBjHQturood+9RtP2aFM+jA69c2Tluk0ec4O2VUcP/e41mtZBKD6MNqUwBg4bKXt2bNMdlEaj0VgkrGG1Sqk/WdkXC3ROiNPKQqPRaCKEFR/G9T723RRuQTQajUYT2/gNq1VK3Q/8DBiglNppOtQF2BhpwTQajUYTWwTKwygAPgT+CJiXST0jIk1ZGEmj0Wg0rZhACmORiIxRSl0mIl+1mEQajUajiUkCKQybUuopIEcp9WvvgyLyTOTE0mg0Gk2sEcjp/RPgW1xKpYuPTaPRaDSXEH5nGCJSDPxJKbVTRD5sQZk0bRSdzKfRtG4CRUlNEZHXgMFKqUHex7VJShMqRrmQBdPydPE/jaYVEsiH0cn9s3NLCKJp++jVCzWa1k0gk9R898/ftZw4mraMLiuu0bRugq6HoZRKBmYBGebzReSeyIml0Wg0mljDSmmQt4FuwErgfdMWEKXUK0qpY0qp3X6O36mU2uneNimlRpiOVSqldimltiuliqw9Chw6cZ7zdQ6rp2suAfS6HBpN+LCiMDqKyCMi8v9E5H+NzcJ1rwKTAxzfD1wjIsOBPwALvI5/T0RGWq2iCFB7ro7fv7vH6umaSwC9LodGEz6sKIz3lFLfD/XGIrIO8FtCREQ2icgJ98cCID3UNrxJ7BjPb28e0tzbaNoQ2tGu0YSPQGG1ZwABFPC4UuoCcNH9WUSkaxjlmIGrbpWBAJ8opQSYLyLesw+fpPXoQId4exjF0rR2tKNdowkfgaKkWiSbWyn1PVwK4yrT7vEiUqWUugz4VCm1zz1j8XX9bGA2QN++fSMur0aj0VyqWFlAabxSqpP79ylKqWeUUmHpmZVSw4FFwK0i4jEyi0iV++cxYDlwub97iMgCEckTkbzk5ORwiKXRaDQaH1jxYbwEnHNHMT0MfAUsa27DbqXzd2CqiJSY9ndSSnUxfgduAHxGWmk0Go2m5QiahwHUi4gopW4F5ojIy0qp6cEuUkq9AVwLJCmlDgJPAe0ARGQe8FugJ/CiUspoJw/oBSx374sD/kdEPgr5yWIAXTtJo9G0JawojDNKqceAKcDVSik77o4/ECJyR5DjM4GZPvZXACMaX9H60LWTYgOtuDWa8GDFJPVj4AIwQ0SOAGnAf0VUqjaCDumMDXQuhkYTHpRI28mAzcvLk6Iiy4nhmksEPcPQaPyjlNpiNUHaSh6GT8Kch6HRRAydi6HRhAe/JikR6eJWCs8Cj+IyRaUDjwD/1jLiNQ/vOkLRrisU7fY1Go2mOVjxYdwoIi+KyBkROS0iLwH/GGnBwoG37Tratuxot6/RaDTNwYrCcLgry9qVUjal1J1AqygJ6+10jrYTOtrtazQaTXOwojB+Cvwf4Kh7+yf3vpjHsF0bjk7vz9GWJ9Zpiya0tvhMGk1LEVRhiEiliNwqIkkikiwit4lIZQvIFjJnL9SHrSPQHUvbNKFZfSb999doGmOlllSOUmqVsRCSUmq4UurJyIsWOl/VnAtb59YWO8tQaYsmNKvPpP/+Gk1jguZhKKXWAv+Kq8z4KPe+3SIytAXkC4mBw0bKnh3bwmLy0bH7lzb676+5VAglD8PqintfeO2rD12syNM5IS5sX25z7L42TVx6tDZ/k0bTElhRGMeVUgNwJ/EppX4EHI6oVDHEprLjzHi1kE1lx6Mtikaj0UQVK8UHf45rve2BSqlDuNbinhJRqWIIl5YMkPKu0Wg0lwhBFYa7euwk99oUNhE5E3mxYofxWUm8PD2/TTl+NRqNpikENEm5k/WSAETkG+CCUmqWUurLFpEuBtC2bI1Go3HhV2EopX4C1AI7lVJr3WtvVwDfB+5sIfkiio6112g0GusEmmE8CYwRkd7Ar4CPgAdE5HYR2doi0jUDh1NYX1LNupJqvwpBx9prNBqNdQIpjDoRKQNwK4j9IrK8ZcRqPgUVNcxaWsTsZf4VQltMTNNoNJpIEcjpfZlS6temz53Nn0XkmciJ1TwcTsHhFOZPGYOyKb8KIdg6CTp5qzGReiexfl/9v6DRBJ5hLAS6mDbvz0FRSr2ilDpmlBXxcVwppeYqpcqUUjuVUqNNx6YrpUrd23Qr7Rm1pAoqarj/tS3Exdm4Oie5yV9wbbJqTKTeSSTu63AKi9ZXNPu+xn1mLSnU/wuaSxsRidgGXA2MBnb7Of594ENAAWOBz937E3E52BOBHu7fewRrr31qtmworZZ6h9PzszmE6z5tiUi9k0jcd0NptQx68gOZt6asWffdUFotg3/zYbPvo9HEIkCRWOzTI76mt1IqA3hPfNSeUkrNB9aIyBvuz8XAtcYmIvf6Os8fqQOGSOmeHew4eNKn6UCbFS4ttDlKowlOuGtJRZI04GvT54Puff72N0IpNVspVaSUKqo+Xs0f3t/DjCWFrC+pbhQya5g95q8t56U1ZdTVO8P/RDFAtMKFw9luOO4VrhyacN1Hh3FrWjvRVhi+voESYH/jnSILRCRPRPKSkpLp3b094hT2Vp1mxquFLFxX7vmCjs3syYMTs/nLp8X86aNiFm/c77lPOL7MsdIhRMv30tx2ze+vLfqP2uIzaS4t/EZJeUVINULCEyV1EOhj+pwOVLn3X+u1f02wm508X8dLqyv49Q255PbqgkOEv3xawrD07p4R4swJmWRf1pmP9x5h2pUZns7J4RTuf20LC6blBYycCoTRITTnHuEgWuHCzW3X/P5a4hla2tSkw7g1rR2/Pgyl1FPuX3OBfOAd9+ebgXUiMtNSA4F9GD8AfoHL+X0FMFdELldKJQJbcDnMAbbiSiKsDdRWv9xh8vC8vzNrwgDsNsXGsuMoYJyXOWFj2XFPxwQwe2kR86aMweYOwW1q56Ft3c2jpd+f+f8gmgpeo4kmofgwrCyg9Anwj+IuOqiU6gL8TUQmWxDkDVwzhSRc64E/BbQDEJF5SikFPA9MBs4Bd4tIkfvae4DH3bf6dxFZHKy99r2zJeOeubx8V36DDsC7IzJ/djiFxRv3c/f4/sTH2fxeo2l76L+xRhN+p3dfoM70uQ7IsHJzEblDRFJFpJ2IpIvIyyIyT0TmuY+LiPxcRAaIyDBDWbiPvSIiWe4tqLIA6NWlPfOmjGk05fe2HZudmAUVNfzl0xJe3lDh00Futjd7+yhixckbK76T1oYuLKnRhIYVhbEM+EIp9bTbTPU5sDSyYjWNY2cu0C7O5ukAjI40PyORBdPyGN23B/PXlnO+zsGafceYt6aMi3UOLjqcPPtpSQPlMLpvD344IpXzF+pZV1LN+ToHTy7f1SB5K5xOzGD3CqQUWqsztTUrutYsu0bTVCzlYbgzsCe4P64TkW0RlaqJ9MsdJqV7dnhMS4aN+sGJ2cyckMmi9RX88cN93JHfh7e2HuSiQ/hxfjrLtx7k9tHp/OHWYZ5r568t548f7sNug3Z2G7ePTOONwq+5I78P/3b7sAamrfyMRAorayPq/whkb2+tppXW7ENozbJrNGbC6sNw3/AqIFtEFiulkoHOIrI/2HUtTfvUbPl07SYm5CQDro50wbpynv20hAVT8xAFXx4+zdSxGSzZvB8liruv6s/SzZXMXVXa4MtfV+/klQ0V5FzWhbg4G/kZiSzdXMnd4/t7TFlGB90SnYdZKQCtUkF401oVHbRu2TUaM2H1YbjNUI8Aj7l3tQNea7p4EcRrKVW7TTGkdzeUTbH3yGl+/vpWhqd3Z8fBk7y4upwRfbvTId7OzAmZjcId4+Ns3HdtFtcN7sXVOcl0iLdz7zUDiI+zNTIBtUS4pLffpTWaoLxpzT6E1iy7RtNUrERJbQdGAVtFZJR7304RGd4C8oVEv9xhsm/XdrYeONEoImp03x4s2bSfrOQuFB87zaCUroCi+OhpZlyV2SBCKhhNHV3qUhVtC/130LQFwh0lVecuUCXum3dqjnCR5NiZCyzdXOkZfZu/0FsPnODZlaXc9/oW/uvjEj7ec4TZywobZXxboamjS38zg1AdqHp0GxuEa6anHeia1oKVGca/ANnA9cAfgXuAN0RkbuTFC42Bw0bKzm1bPQ7ogooaZi0p5OaRvblxUIprxiFQfPQ0ub268rPXt3DrqN4NnN2RHDX6u7d2oLZOmvu/Es4qAxpNU4mE0/t64AZcNZ4+FpFPmydiZBg4bKTs2eEK4DKil556ezdvFH5NOxv8y40DmTkhs1HyXlM770AJgdEwVbU12vp7Mf7XwlFlQKNpKuHO9P6TiDwSbF8sYERJ2WzK0+nnZyTy8oYKnCK8uLo8qCIIpZPyVi7rSqqZvayIBVPzuNodqaVpOm195tXWFaKmdRBuH8b1PvbdFJpILceuQ6cYkd6dW0b0ZkR6dwora5l99QDuuybLUiRTKP4B7+goBSCun9ou3XwiHX0W7b+R9kVpWhuBig/eD/wMGACUmQ51ATaJyJ2RFy80ElKzJeOeOfzD6HTeKPya7+Um83lFDQun57fICNU8YoyVyrUthdH5+ir2GGsYfyenU7hP+w40lzjhmmH8D67KtG+7fxrbmFhUFga3jUrjiR8MdimL/bU8NCkn5BFqXb2T+WvLLS+wZHSWgGfE2FKlrKM9SjYoqKhh9rIiZrWC/BBDmQvocuMaTQj4VRgickpEKoE5QK2IfCUiXwEXlVJXtJSAoVJXf5HF6ytYV1zNL67NYuaETMBlD6+rd1rqXBdv3M8fP9xnOdzWHF7pS3k0l0jVkQqnshmb2ZMFU/NYaLEDDrXtpshqvsbhFNaXVLOupNpTW2x8VpJnNhhthRvJwpYaTbiw4sN4CThr+vyNe19Msnz7UZ79rBQHcODEN8B3nerijfsbdey+vpB3j+/PYzcN5O7x/S21mZ+RyIMTs8nPSGxyB95UpdCcmUw4M8btNsXVOclMyEm2pCRDbbspspqvKaioYdbSImYvK6Kwsjbmsua95YgVuTSaBohIwA3Y7mPfzmDXRWOLT8mSOxduko93HZaH/7ZdBj7xvmworZZ6h1M2lFbLhYuOBp8H/+ZDWVt8rNExb4zzzcfM+9YVH5PcJz6QNV8elTVfHpUXV5fKuQv1fu/nC0OeDaXVltoPB4HuG6k2m3r/pshjvsb4O63ed1TWFR/z3CfSz9kUWQPJFSvyatoOQJFY7GOtzDAqlFIPKqXaubeHgIoI6rBmsbmslrJjZ7lxcAoPTcohPyPRE40SH2dr5GNQ0Gj24Y2v0Z55nysFHvYeOc29r21hzqpSlmzaz4xXC9nkNk8Fw99MIdyhl+aZTKAonUiNcJtqsgsWUeRrhma+xm5TTMhJpp3dxn2vbWmwNkosmKW8n8/f8+qZhyaqBNMowGXAX4FjuFbN+x/gMqsaqSW3hJQsGfDYe5L9+PuS9dj7kvvkBz5H7GaMEdu5C/Uyb02ZXLjo8HuO9wxjbfExWVd8zDM7uXDRIau/PCovrS6VVXuPSO6TH8ja4mMB2w9GoJlHJO8XqZFsuJ8n1Pv6eq5IyRQJrM5EWqJtTduAEGYYUe/kw7n1yhwkz36yT55bVSKf7Tkia92mhwsXHUGVwdriYz47DeNaXyYmXx3NuuJjkvP4+/LCZyWe9ptDU7+ksWrSCGf73ianpt63JTvhcN+7JZVda1KsGuuEojD8mqSUUg+7fz6nlJrrvVmZvSilJiulipVSZUqpR30c/4tSart7K1FKnTQdc5iOvWOlvRPnLjJ3dRnPry4jPt7O1W4HbKCoJ2OKr/AdYmlc+/Q7uxuZmHyZkRwiOESYu6rMYwrxhdUomHAXOox2spi/9oO9D++Ip41lx9nkzgQvqKhp1nN5XxtJs0+4791S4dst3ZYmRvGnSYCb3T+n+9qCaSLADpQDmUA8sAMYHOD8B4BXTJ/PWtV6xtYrc7B8uvuwrDWZiazOMLxHl8b1hqlq1d6jlkxMxgzjpdWlAUeRkXZyR3sm4U2wGY+/GZ6B+X2ZAxbWFR9r8kwuGrOwWPu7aDTEgkkKuBJXoULj82PAYwHO3wRcb/ocssKIT8mSP7y/U344d608+HqRZDzynny8s0rmrSmTM+cvyourS2X1vqMuX8O+o/L8ZyWez97RU/PWlDXowKx+0cNxXihT/2hEUDUFf89k7F/njlbz154/81NzzCTaxKKxQltX8qEojEClQd6l4QJ23jOTWwLNXJRSPwImi8hM9+epwBUi8gsf5/YDCoB0EXG499UD24F64D9FZEWg9sBVGiR1+rMN9vVP6sj+4+f4Xm4yq4uribcr/vmGXP78STF1DsGu4J/y0nlnexULp+d7ImbyMxL53B0B5W3qMEcuge/lUpsT3RTKtZEq0Bfu+/p7puZGgbXUe24uLdmWJry09SKY4SoN8t/An4H9wHlgoXs7C+y2IoePff4U0E+Atwxl4aav+yF+CjyrlBrgsxGlZiulipRSRY5zp/jhkGSGpXXh1uEpABysOccd+X149sej+HFeOs/fMZp6p5MX7xzDP+WlYVOKFduqPCVEzPbsvYdPc78pBNPAOyHMl03a1/5I+C0CheM2JzO6OfbqYCGuZprrUwmn7yKS6HDY1ov23ZgINgUB1lnZ5+McyyYpYBswLsC9XgV+FKzNXpmDZc7KYlnzpcvMZLZvG+aHh/+2Xfo98p686PYx+ErmEvnOF/G8KdrJVwJgKHbwpppAmmJrN7dldUq9obRaBj35gcxbU9Ysu36g+1yqtHWzhqb1QpgT95KVUpnGB6VUf8DKYg+FQLZSqr9SKh7XLKJRtJNSKhfoAWw27euhlEpw/54EjAf2BmvwxLk6nvm0lHuWFPJ5RQ0TcpI9kVJjM3syb8oY+iR2op0NBqV2bZTMtaGk2lN0UFxtc6D2HPe/toVF6ys8UTnm0hKhYIxU8jMSg47+zaP0UGYx3m2ZK+cGG92OzezJQ5NymLuqNKS2gt2nKbOdtka0o9M0mrAQTKMAk4EDwBr3VgncaEUbAd8HSnBFSz3h3vd74BbTOU/j8lGYrxsH7MIVWbULmGGlvfiULBnx9EeS/fh7DaKZvCNxvEe+xvEXV5dKv0fe8xyft6ZMBj7xvjz61g4Z5I6Q8ldixNs5vs4r6idUR63Z+e6vbEm4HfGBzm1OKY/W4lz2NYPUaNo6hDtKCkgARri3BKs3b+ktPiVLHnz9c7+moGCRON7ht/46EKN21DqTucs4Zg759LXfqonIUFiDLGSrtyRNMa2EyxwTbrOOv7+dd4ScRtOWCavCADoCTwIL3Z+zgR9abaAlt/iULJmzsqTRCwl3R7O2+JjfnAx/4Z+R6GhbMkPZIJqzhea0bcWn5D1A0DMNzaVAuBXGm8DDwG735w74qGAbC1tCSpas2ns0fG/Si1BNFpHuXL3vb6W95iqVaDpvm9O2r3cT7H6txZSm0TSHcCuMIvfPbaZ9O6w20JJbQkqWPLeypEFBQH8Z3E2x0RuRUw//bbvPrHFvIt25NmWGcal2gtE0pWk0sUwoCsNKlFSdUqoD7hwKdz7EBQvXRYVnV5Vy9+Iv+O3bu5i1pNATpbNofQWzlhQ2KmNuRPBsNNUl8ocADoE3iw5aWo3POzIm3NFCVktim7lUY8qbEqWkI5s0moZYURhPAR8BfZRSrwOrcJmoYhKnUxBokIxXUFHD3FWlPDQph7vH92/QYQYrPmhmfFYSi6bn8cjkXMur8ZmJheSt1t4JWlW6/s4Ldb32SKNDjjWtCb+lQQCUUgpIB84BY3FlbxeIiLVVgVqYfrnDZPHyldhsCps796Kwspb8jEQKK2t9lmVoyZINwdoyHwffJUcudayWafA+z3i3Ow+e5E8fFfPYTQO59xqfxQNalLZedkIT+4SrNAhu+9YKEakRkfdF5L1YVRYAx85coLT6LOPcI+jP3SN6f4l2RidirMXd1FGrVYKN7q2UHLnUsWpS8z7PeJ8DL+vCHfl9mHZlRgtIG5xL1USoaZ1YMUkVKKXyIy5JGEiIs/HMx/s8fgohsJnJ6EQCLc/q6/xIdeLmzkN3JL6xalLzPs94n/Y4G+/sqGLrgRMtIW5QWruJUHNpEdAkBaCU2gvk4srw/gaXWUpEZHjEpQuRhNRs+el/vM78qXlsKj/Ol4dPM3VsBq9//hV3j+9PfFxD/WieYfgzWZmpq3fy8oYKBqV2ZUJ2sv6St0Ja0gR5KdEa3mtrkDEahM0k5eYmXIsgXQfcDPzQ/TPmSOwYz4t3jiE+zkbxkTP86aNiHnhjq2e1PYdTWFdSzfqSahxO8Yzu4uNslqKZPq+o4ZlPSrhv2RZdI6mVYtQVs2KC1FinNZhQW4OMsU6gJVrbK6V+CfwrrnpSh0TkK2NrMQlDoFvHdgDMX1vOnVf047GbBvLcHaN57KaBTLsywxVau7SQWUH+afz9Ywlgs8GDE7NxupXF7KVFbHKH5Zo7oFiIxvGl0FqTkouUrLrjCD+twYTaGmSMefwlaODK8H4NuBdYAcyxmtwRrS0hJUumvbxZ+j3ynry0urRBcoqRsPbi6tJGpcy9CVaAzygsaNSm8rW86Lw1ZZ5Chs0h3EX/rCbuxULSmj9Zm1t0MRaeTdMY/XeJDoQj0xvYZfo9Dthq9abR2hJSsmTAo+9J5qPvyQurSuXCRYdn3edw1gWykmF97kK9PPrWDjl3ob5ZbYWamR2syKHV7PBQ221u3axg9/QlW7BikpdqVntrRf+9okO4FMbWQJ9jccsdMkLakVyWAAAgAElEQVRWf3lUXlxd6qk6mvvkB5L7RMtXfA3XP3+o5T9C7aybOor3V+l1Q2l1xL/43uXq/bWjR6ytC/33ig7hUhgO4LR7O4NrbW3j99NWG2jJbcyYMSLSsG7U6n0uBXLuQn2TzRNN6aQj9c8fzs643uH0zMBCldNfpVd/M4xIvI9w3jOWOqtYkkXT9glFYfh1eouIXUS6urcuIhJn+r1ruH0p4aD6zAXOflvPovUV5GckYrcpio+c4fnPyli6udJnDan1pdXMWFLI2uJq5q0pY82+Y9TVO1lXUs2afcdY746qmvFqIZvcDtiNZcc5X+fwOLV9OVEjFV8fzrW2CypquP+1LdhsKqCcvpzPhhyDUrpy83Pree6zUs879/XskXA0h/MdhypfJIMHtFNeE6vYn3766WjLEDYe/Y+/PH2w1ziWFRygc0IcdQ4nTyzf5akhNaZfD8Zm9sSmFJvLXV/KzOROfF5RS5xdMX/dft7fdZiO8XaeWLGLt3dU8d7Ow2Rd1pmC/bXcPKI3VSfPM3tpEUdOfctLayvonBDHLSPTGNX3u3tHEptS9E3siE0pHE5hc3kNad07WGrXeOZRfXvQN7Ejad07WJLb+zqzHD9ZsJmdh05TUFFL54Q48jISfd7DalvRIlT5fL2TaMmi0TSH3/3ud4effvrpBVbODZq415rolztM9uzY5knUs9uU30Qd76S90X17sGTTfgamdGVcVhIFFTU4neLKUlRgV4px7lo/BRU1jO7bg6WbK30mBLYUodYhamriUqDrTp27yJSXC5g4qBc/uzYrau8iGOFO2tJJYJq2QiiJe21KYeTl5cmcv37k+RJbzeQOdF44i8PpTit66CJ/Gk1jHE6hXccuJc5vz+ZaOT+iw0Gl1GSlVLFSqkwp9aiP43cppaqVUtvd20zTselKqVL3Nt1Ke5XHv2HG4s9ZtL7CY2OesaSQVzZUBLQJB6oplZ+RyIMTs8n3YWoJtdT2JgtrboRCMBt+uOzs0Uz2C1fb4Uraak2JjxpNMAoqamjXvVeW1fMjpjCUUnbgBVylRQYDdyilBvs49U0RGeneFrmvTcS1DscVwOXAU0qpHsHaPHOhntzUbsxdVcqi9RU4nQICA1O6BuwsjM7Ee60McL3Qv6ws8dnJW3VOGucFK4YYbsLlPI2mE9Zf26F23L6Ua1M6f+2Q1rQlxmb25OLJo2WWL7AaThXqBlwJfGz6/BjwmNc5dwHP+7j2DmC+6fN84I5gbfbsN1BOflMn89aUyaAnPwia2GWFdcXHJPcJ1728aW7GcaRpapa4v+Vtm5Jt3lz8yRSO8OKm3CMWQl5jQQZN24EwL9HaVNKAr02fD7r3efOPSqmdSqm3lFJ9Qry2AWcv1PNv7+0l+7LOzJ+axxU+RvKh1ngal5XEy3flexzeZppSarulTBrB/Bu+5DBGzy9vqPCEEXvLbwUro3Cr78Fou6CihhlLCtnolikcJqam3MOQB2ggf0uaqvQsRxMtIqkwfPUu3t+md4EMcZVKXwksCeFa14lKzVZKFSmliurPneJvWw9yrzu3oLCyttEXa/HG/Z7qtVYI97rcoXzZm9NWsHZ8HffkVqR2BdX09q10xKF2ek6n4HCKy8wYJpqTx+Etf0t24rqIniZqWJ2KhLphwSTldb4dOCXNMEn1zRkqz60skVV7j/itH3XhokPmrSmTCxcdTZq+NdUU4su0Eqw4XrDSF1baa0oJkWi37wtv02C06w6FWrJFo4lVCEdpkOZuuAoWVgD9gXhgBzDE65xU0++341ovHCAR2A/0cG/7gcRgbRqlQYzOZG0YfBjeNLVjCKVyrNXiepGk3uGUdcXHZG0TyoaIhL9D1x20RhMZYkJhuOTg+0AJUA484d73e+AW9+9/BPa4lclqYKDp2nuAMvd2t5X2coeOkHqH0zOLWPPlUUuF9ZrT+ZhnDoE62FDqK4WjvlVzaW6HH86Ktc1RXBqNJjChKIw2lbjXPjVb/vw/HzAktRv3vbaFByZmMTytO+NMDmdjpbX7X9viSeIKJanL25m8vqSaWUuLeGhSNnNWloKCl6fnRzw5LFyJaP6c41b2AxFPHNxYdpwZrxa22HvVaC41wr1Ea6tBgGc+LWFv1WmcIsxdWdqgsJ7hmFQ0zIcIxYno7dwUAAWDUruycFoeC6b6v09TnciBiv/lZyRGxAnvzyFsPr8lHL1jM3s2eK86cU6jiSJWpyKtYevWJ1f6P/KefLbniKzae0Qefmt7gwWMwmHesGrO8rXfbOYxO9+DmW8CmYesmo78mcRCLW8e6kJJTTFN+cu9EGn+KnwajaYhxEgeRotz/qLDNaOIs1F67CxvFh7k9+/u8YxGjdGysQ53MHyNZr1H3lZG4gbmmYw5vDfYSD3QDMjq7MhXG1bLm5sxP6+VkOOmzEIClWrx97w6N0GjaQGsapbWsPXNGSqr9x31LJz08N+2S+7jrnW1jRHxi6tLJfeJD2T1vqNBR6TNcfwGG/H6m2GEwwHfXKd7U/H1vpqSDR9ohmHl+kuJS/W5Nf4J9X+CS3WG0SHezoTsZAoqarjvtS3cNCyVX90wkLmrSj029+dWlfKr63NQQoNsZgPzKNnfaNaKHT1Y7aL4OBv3XjPAUw7c4RQ2NbNAYaBRti95gmWgh+ovMBdqNK4FLCXHmWU35IqPs1lOrIvUglVmYtF/omdWGm8i+T/RphRGZc03bCw77koTF9caFjMnZHo6/bGZPVk4PZ+ZEzJRNgWqcfq4r47L2+TS1E69oKKGWUsKPdV0G7S5rIhZzSxQaFZwoXZu/kxWoTxnYWUtc1eVUlhZG/K1oQQeRKvjDmfJk3Chs7413kT0f8LqVKQ1bPEpWbJy12HLzlhfDnDDVGReA7ze4fSYuAY+8b7fhEArZqiH39ouuU98IGtNxQy9nc/eJhl/phnzM/grzudLVqvmqVCntubzfb3HcBGtLG8r7yPaGegaTahwqZqkAD7ed9SSecJuU9hsivtf28LGsuOeUaExSjavAW6MLN8sOshto9IZn5XE2MyerjXBS6obrJHta91w43hhZS0rth3C4XSyp+pUA2f81TnJnpX+jBwLw+lrdv6aiycWVNQwa2kRs5c1dhAbowwFlmYOTTFZ+Xqnxvm+3mO4aIlRtZWAh2jJ1lRi0aSmaWVY1SytYYtPyZIPt1f5DH31NZswzltnqpnky+FqzDBeXF3qcVLPW1MmuU98ILlPfOAZTXq3a4w2Daf7hYsOeXF1qbywqrTBKNRbjhdXl8rqL4/K6n1HZV3xMTl3od7jIJ+3pkz6PfKdI9/fDMP7GZvr7A515NwUx3Us0RZnCm3xmTTNh0s10zshNVt+OfctUnokMGdlObcPTyErtasn8xsFC6bmYbepBtnJdfVOXt5QweCUrlyVk4zDKSzeuN/vet0by44za0khD0zMZmjvblyR2ZPCytpGS7w6nMKi9RXMWVnC/Kl57D18mjkrS3jpzjHsO3rGc3/jfjeP6I1ThBXbDvEPY9J5d3sVC6fnA3iyuvMzEgPKZuAIkKm9qew4gjVndLD7tVVi8XmbK1MsPpMm+lyymd5dEuK4YVAv5q4sB2D5ziP818cl7D18il9OymbenWMQpzSKjiqsrGXOqlLufW0LBRU1nhyJ37y9izX7jrFm3zHWl1RTV+9kY9lxRvftwUOTcphxVSbjspJYvHE/s5YUNjIL2W0up/vC6fkIMHdVKbeOTEMUHucwuMwYN4/szZtFB/n7tkM4Bf6+9RC3jOxNfkaiJ/podN8eFFTUMDi1a0hRR977DTNWKGailohCijZmk00oz9tSpp5gTvdgcvhby0OjsUqbm2G8v3I94hDmrC3je7nJtG9nZ0jvbvz89a0smJaHwynMXlrELydlc89VmRRW1jK6bw9e3bS/wQzjN2/v4u9bD7nCqJRLs/7D6HTe2VHFgxOzmbuqlAXTXEp59tIiHpyYzd3j+zeYYZjxnm0IIE5Bub/EG0qqmbGkkNtHpdE3qSNzV5ZhU/DK3Zd72njguiz+/EkxIrD47suZkJPs912Ee4ZhhdY+gm1qfa5w1fUKRrD3a1WO5srb2v/OmoaEMsOIut8hnFt8Spb8y5vbfUYdGZ/PXaiXR9/aIQOfeN+TxPfi6sY+hRdXl0r2Y+/JC5+51td4+G/bZdCTH3h8Cd7RQFbW1zD7DjaUVnt8IGb/hBHZ5CtqavWXRyXr8fcl67H3G0RZxQrNSdyLBZoqa6w8Y1OSJJuC9oW0LYiV8uYtvcWnZMlzK0s8nXHWY+83WnDH6JTnrSmT1fuOSu6TjbO+vZ3VG0qrPcrCOxu7qV+eeodT1nzpCtXNfvw9j0Pel6PYvD/U2k/m6yPdoQWrnxXJdjQth37/bYtLVmF0Sc+RmjMXPDOEnMffl4ff2t6g/IaV2YGvKCtzxNS8NWUy6MkPGhUR9IevvIp1xcc8M5ysx1yzHaO9NV8elazH3pfP9hzxRE/5mt1YJZojwkh0LnqEq9GEj1AURptyel+odzJtcQEL1pUzdWwG47KSeLPwIPcuK8LhFE+EUX5GIuBa33vOyhKP89kX3vbagooa5q4q5aFJOYx1R0eZHdgOp7CupJq1+46xzp2j4V1M76m3dzNzSSF//qSYm4f35pYRKfzpo2LueqWAT3YdZsH6Ci46hfd2HWa2O/v7gYnZ/OXTEl7ZUNHA8WnOy/Dn9IxkboBVR2s4bd3NeZ7mOKh1HoPmUsf+9NNPR1uGsPFvf37u6fOZ17GxrIaaby7w2ZfHGNy7KwX7a+mcEMfWAyf4z4+K6Rhv56JDeHLFbh6alMMtI9Owqe86tM3lrg5+VN8eHDxxntlLi2hntzGqbw/Se3RkVN8e3DIyDYADNef4wfDejBuQhE0pNpe7opDe3n6Id3ZU0c5mo2eneG4amsKRE+fBJnyy+wi5KZ05dqaOvYdPs/fIWde9as/z7q7DfH3iPADnLtTxX/9nJFcOSGLbgRMUlNeQ0q09D03MYXxWEiLw5PJdzFtXQcd4OxcuOpm5pIh2dsXovj08z2RTir6JHT2fHU5hc3kNad07NHjupmB+V30TOzbrXlbxfh6w/kzNkTcaz6rRRJrf/e53h59++ukFVs5tUzMMA7sN+nTvyIt3juGhiTk8MjmXu8f3Z1BqV9rZwOGE+nonL00Zw93j+3syqH0VHczPSOSWEb159tPiRvWljCKHRqlvcI9+p+bxy+tzQOCZlcXc9/oWPt17lD+vKmVtSQ2D0rqyq+osToGu7f3/Cb46cYEvD59m8cb9zF1ZwpVZPXmz6CB7q057MsJXbDtEnE0hTnCI4BTh2U9LGoVemkfH4SxOZmW039yRuZXrN5Yd555Xv2DBuvIGf8umyOuPWM7i1mhagjanMBQwISuJF1aXUnz0DA+8sZWhvbtRWFnLuAFJ/MuNA3n+s1JmLyviy8On2VR2nBlLCnllQwWzlhSyYF05G8uOezqFxRv3886OKn5+XQ7vbD/E+TqHpy1fxf4Ars5JZtaEAfzzDTksmOJaLe63Nw/h4Rty+eWkbDonfPfaT33rBGBwaifGZXRv8CyDUzqRndyFuatKuW1UOlMu70c7u0KUKzRYnMIvr8/h15NyeH51KXalWDQ9n0XT8wOuFxHOjs+Kyam5CsrK9QpwOmHOylKf62gEkzfU0iehoE1ZmrZCRPMwlFKTgTmAHVgkIv/pdfzXwEygHqgG7hGRr9zHHMAu96kHROSWYO0lZQyS3yxcwUtrK3hwYhYDe3Wl+NgZsi/rws9ed63hPSE7mfWl1Xy4+zDv7jjMLSNSWb71EC9NGcOne4+yfNshz/rRALOWFPLQpBwqqs/yZtFBfpyXzi0j0xrFoK8pPsaspUX8w6g0/nDbMAoraz2x7mMze1JQUUN+RiKzln7B2pLGHZnClfKRkdiBytrznv0P35DL1yfO8fb2Q8yfmofNpsjPSKSgooa9Vad47rMyXpoyxpO97i9LPVDsfF2901L2eFMJZ4Yy+F5H3JxfYviWzOdYzWF4cGI2MydkhsXnYrTp8FpDXqOJJWIi01spZQdeAG4CBgN3KKUGe522DcgTkeHAW8D/NR07LyIj3VtQZQFw5kI9B0+ex+EU6h3CzKVF/N+Pinl1YwV1DuHDXYcBaGe3sWLLQbq1b+cqwzE6HZtSvL39EA9NymbelDE4neLJ6J52ZQbXD07h2pwkJg7sxYxXC3li+c4GjuY9Vae46BDeLDrIb1bsYnTfHjw4MZsR6d15fPlO7ln8BY//706fygK+K7P+lUlZTMjqyaCULqzYdpBbR6ZxhbtEe0FFDXuqTjF3VSkPTsxmfFaSZ+RrXsnPTKAsX3/XGDR3hNxcx7e3GdDfGuQTcpK5OifZ5zoagWYpDqfgcAq/uC6LOSsbm/N8nR/ofRjHDSXkvYa8RtNaiaRJ6nKgTEQqRKQO+Ctwq/kEEVktIufcHwuA9OY0aFOQ1CnBUy3Viasj3lBeS2ZSR97ZUeUZZQ7s3Y3Dp7+ld/eOvL39EHabYuH0fGZfPYB2dhv3vbaFp9/dzTOfFrNgbTmzlxWxsew4JcdO8229k78WHmTBujLmrynnrsVfsL/6G67KckVfvVl0kJc+K+PPnxRz37Ii3iw8yAWHsHz7oaDPIMDQ1C788w3ZbD1wEnucjV9en8s7O6oorKxlfWk197xayDOfFPPQpByPD8bo9HJ7deFRt8/GF+aO0+jYpl2ZwWM3DQx4ja91PKJBU81pga4zlqod0rsbC32Y83ydH8hEZhw3FMU4k0LXaFozETNJKaV+BEwWkZnuz1OBK0TkF37Ofx44IiL/5v5cD2zHZa76TxFZ4ee62cBsAHvX5DHp9y/m2uwk1pQ2XrP7e7nJzJ+aR3ycjdqzdfzjSxsZmd6FHw5LR9lt7Dl8EhyAHRQ25qwsxumE5C7xHDlTR6/O8YzLSmL59ioAure3c+pbh2d2cGN2Nz4uPQVAgg0uOBu2P6hXJ748+o2l9/fwjbkMS+vWyMSyYF05f/qomGtzknhpSh5LN1c2KlMSyPRhNs0YHVswU4nDKSxYV86claUsnJYXsCRJoPbMZdIVMC7Mnai3ac2qKSxUk1mw8309sy6loYlVQjFJxUVSDh/7fGonpdQUIA+4xrS7r4hUKaUygc+UUrtEpLzRDUUWAAvAVUuqHbC+rLGySOkaz6byGhatq2BYejd2V51if8059teco/z4OXZXncE8eL5tRC9AMTC1I3sOf0OCDY6erWOFW1kAnPzW0aANQ1lAY2UBBFQW4zN7YLMpkjsn8M6OwwxJ7YrNphp16HeN688X+2spKD/O0s2VzFlZwgMTswMuKWvg3XFZjXAqqKhhUGpXnysUBsNbKRmrCyLw8l35YbXpG6Y1gHuvGWBZIZrNdQaBOnlf5wc6blUOjSbWieQM40rgaRG50f35MQAR+aPXeZOA54BrROSYn3u9CrwnIm8FajMhNVtSpz/LhKye1Jz9lr1Hvuug/3FUb97ZUYUAcTYbL/10NK9sqGB9RcOkPbsCh8BlXeKp+aYOh4+OP1wo4LaRqWT16uIpkGh2YENjB69RCt0wRxVW1jZyqvrr7HwVnQs2GjaumTdlDDavsvBWiOUZRqDj4SwoqGcYmlgmJpzeQCGQrZTqr5SKB34CvGM+QSk1CpgP3GJWFkqpHkqpBPfvScB4YK+VRgdd1pENZTUeZZHUKY5xmd1Yvq2Kn1+bxb3X9qd3t/acr3OwsaJxhrfDrT+PnYmssgDo0t7G8u2HGdCzM3sOnmLykF68u/0QI9K7N/BLbDI5WM3rkhvO3fFZSQ1mCv5W/svPSGw0o/A+1/uzMQvxZ4cPNdPbWF1wQk5y2DtPu00xNK2bZTNQoBUSrfpKIhmOq9HEGpEOq/0+8CyusNpXROTflVK/x1W75B2l1EpgGHDYfckBEblFKTUOlyJx4lJqz4rIy8HaM2YY/ugcb+NsnUsLdIiD8/XNeLgwclnnOI6d/U6Ya3OTKNx/ggcnZvOXlSUhm2+8R9qBRsvNtbdbGYn7umckRt1mWeC7kGhfYbLmMFyjM/d+Fisymq8x/EJ6JqFpTcTKDAMR+UBEckRkgIj8u3vfb0XkHffvk0Skl3f4rIhsEpFhIjLC/TOosjDoEgft4777ssabnvBsnZP2bq9NrCgLgGNn67l9VCpXDUhEAcPTunHziFTuvKIfC6bmsdA00vU1ovXeV1hZy5yVJSzeuN9TQ+vBidmeGlrm67w7OO/RcLARtJWRuK+oIiuLAa0vqfbU47KC+TnHZvbkoUk5zF1V6rMNv1n6AWZqwZ4/nBn0Gk0s0uYyvc/Uw7f133UwdV5mpW9jSFGYOXHuItu+PsU1OUk8/1k5fy08yOuff9XAfGNEK814tdCTVQ6+zUjmztK7QKLBprLjjVYf9CZYJ2jF3OJLYQVTNE1ZGdD8nMZqh/7a8NW+97NYUYbma2KldIjOLNdEikhGSWksMqhXJ0b27sZlXdqzfOtBlIJ/Gp3eKC+ioKKGZ1eW4HQ2DEHz7qiMznJoWjfPPl8dmUDQyKdwdIJGRz7YHfllzGgCmdjGZvZk4bQ8T1ixGX+mIl/vwZ+ZKFj7Vs+JRXRUliZStLkZRmvj6pyelFV/w7Ory/nfLQcZm5mIQ6Bvzw6eCCj4Lht5/hSXicq8b1PZcS4GKLjnbxYwPiuJl6c39I14j06DzSDM5wcrry5g2WRjztz2bjtQtrehIAwZgmV4h3MkHsiJ3pKEouT1bCQ2idW/i1YYUaRXl3g2ldVyy4jeXJ3dE4dASrcOxNsVdmVrFL10/2tbsNkUxUfPcK/bVGOYbmYta2i+sWJP96UMQrXDm8/3Z+Iy2vGO5moKhmKaN2WMz5nHovWuIpLeUV7+MrzD6XMI5gNpqU4glKgs7XeJTWL17xLRKKmWJliUVKxiU/DK9HyKj55h2pUZbD1wgvyMxAYF9AwzjNMp3PfaFk+RPMATdmszdRShRCF5F/fzlSfh737m/RvLjjN7WRELpuZxdQjZ4IHkCZQf4W1uMo49cF0Wg3u7zHGfV9Q0iISy2k445PXev76kmllLi0LOlo8kOkckNmnJv0vMRElpfHPz8BQeuDaToalduH1kKoumujqQe68ZQId4O+OzkoiPszUwrxijxnHuUboRKmqYbq4deJnHfOMvjNXf6NY8mjHued9rWxrUm9pQUh1w9mCwcGrz7eb+RlcOp1Bf7+QX12V5Kvb6yhkZ0rsb9y37bmVDY+ZlXp3QkN3bhNUceY2Cg/7MeVZ8Ri2NzhGJTWL176IVRhToGB/HS+v2s/vwGeLb2UmItzc47l3t1GzS8JU17d3ZhRrG6m1K8RUquvfIaQRhd9Upv0rHMJl5y+YrRDaQAvNlRjLMTbOXFfGsO2TYnIhofjfjspJ4aFIOb++o4rbRaR7zla+qvOGY+htRYE6nBLyXL59RsHcRSWLVTm4Fb+WvaRna3BKtXUZOjrYYARma1pV/HJ7GqQsXyc/owU1DU/nZ61sZkd6dgyfOk9a9g6v8x9IifjAsldtGudbeMJZ/nb20yHPuVzXnuG/ZlkZLhqZ178Covj081/nbZ+C95Kn5s3HdD4f3pkM7O899VuZziVLv+3svcztzSRHv7qxiTD/XtYGWO/W1BOvm8hqeXLGbW0b15sYhqTz/WRlj+vVgfFZSo/b69ezEqL49iLfbWLqpkttGpdOvZyeGp3enc0Icd4/v7xm5BXovVimocMs2ordnrRRf9/L1XMazRWPp19a85Oyi9RX88cN9dE6II88rv0gTGqEs0ap9GFHAWCzpx/np/Mftwz0mEaMelNMpjWzdDqewoaSaPVWnESW8uLqcF346muKjZ0KuztpUQvWLGL6QK0y+BKPyrrePxkrbi9ZXMGdlCS/eOcbz3MYaGb7uZ9W3EMpzQ+P6Xuas8dF9e/Dv7+/ltzcPoYN75ticmlbhIFzvIZaI9KJflxLahxHDjEjryoPXDSDOBjcOTvHY0cUpHtv8uKwkXr4rn3FeFU/vfW0Lz64qYe5K18JJdptqkJAX6cgKX3ZVf6YBsy+ksLLWEyJrrERYWFkbUka5kVsyf2oe+46cZs7KEgoraz3P7H0/f/JCaO/JYwpb+l1Umnk9kfUl1awvrWZ31Snuf20LD7yxlTcKv+b37+5p1N4mLx9HMDnDga/IsZZoN9LEx9m495oBWlm0MHqGEWE6x8PZOtesol9ie6pOfsvCafnExdk8o2KnU5i5pBABFrlt3BvLjiNOQbkVyueGY1aBXSmPMjGPEJsy6mruKPOlNWX86aNiHpmcy/3XZgW9dyiVdH3JOn9tOc98WsyvJuYwrE93z4wl0Ojdu2aUIYOVWY65OrARlWbIX1BRw4xXC3GKoIBf35DLnVf08zvDMM8iWyqhLhJLz2raFqHMMLTCaCEMM5RNwaJp+Vw36DLPl/mlKWPYU3WKOStLeWBiFjalPBnddrviV+4yH76+9GYlYV5H3GqHFGoZb+8OeNXeo9z3+hbm3TmGG4amNDrXbMrx7rj9nRuo875r8RdcdAg/zkvnvZ2HfZZqNysCo1MH+NX13xUiNEbexuJToRRO9H4XDhGPEg/UIUfDBNRShR81rZdQFIZ2ercANmBwSheqz9YhQLs4xcRBvUjt1oGEOBs3j+hNfkYiCXE2nv20hM/31/LQxBzGDkjknvH9uWVEGglxNuasLGF0vx6kde/A5vIa0rp3YNH6Cv7zo2KOnfmWKWMzGNWnO/VO4avj33Cg9hzpPRo7Wc2E6vTdXF7TwIH9/q5D7Kk6Q0I7xQ1DUhud68/xbX4Gm1INHMIOpzQ4Bq5O7kDNOb4/rDcDU7vwy0m55LmLDBrnGO21syse//suhqd1QwE/GJZKZnJnnl1Zwsg+Lqf45vIanli+i4cm5XDLyDS/z244/jeVHeermobv0zUr/l0AABb0SURBVKYU/Xp2onf3Dny85wjnLzjok9gRERrJb5zvy+ltBV/vxAr+Aghaq7NbE35CcXprhdEC2BTUflPH0LSuHD1zgX2Hz3Dk1Hk+23eUpZu/4sipb7kqO5k6h5Obh7sibXp1bc9vVuzCble0j7Nzy8g0xrg7yIIK1xrbcXZF/UUnm/fXsu/wafIzEl2r9C0r4u3tVby/8zCHTpyj5OgZRvbp4XM0GawT8+6o0rp3YFSf7tw8ojfjBiRRdvQsmypq+P7QVC73CoM9UHOOHwx3nZfeoyPD07qRmdSZnp3iOXTiPPeaIrzM7RhRYiP7dKdfz06Aq5O7d9kWbh+dxj+MTqed3dZI7sSO8ew8eJKrs5P5cM8Rsi7rzG/f3sNto9K4rEsC7+44zE1DU/hk71EmD00lv38it4xMazDy9tUpe5TkjioS4myM6tujgSJ7/O87WbB+P+9uryK/fyIHT5z3dMhmxehPkQR632YZwtXJhyMyTNN2CEVh6OKDLYBD4IGJWQxP78G+I6eoqP6GN4sOAjAivRtvFH7NkdPf8nlFDQvdPgyHU7htVDpvFH7NW1sO8updl3sipoxqtM98UoxDXKsE2mwKhwg4Yd6UMSiBD/Yc5s1CVztxNpeTsIFcQcwtvuzuRqKgwexrBtAuzuazUOJ97uuMe7eLszH3s1JQsGBqXoMcCrN5SAFI4AKLhpNcnIK4/Trv76pidXE1l3VNYOHUPBzSsITIy3fls/PgSf70UTFAo/fhr2ifUQjRZTYsYWhaN8/xgooalm89BMCPxqQ3KvZo3HPelDHsPexy1i/0kYthRYZwVcJtrUUVNdFH+zAiRPcOdk6ed9CjYztOnLtIRs+OHDp5nl9NysEpTv77k1Kuzkokv19PDp46z7s7D3PriN787tahHod1Xb2T3769ixXbqlgwrWG5jbp6J79ZsYvlWw/y0PU5DEvrjgJPJz0+K4m6eieL1lcgCLMmDGgUehuolIf30qyGXyA/I7FBuQ2HU3w62n2FopqvNTurN5YdZ+aSL7h9dDp/uHWYJ1Q2mD9jxpJCHA4ngkth/Oy6AWw7cJI5Px7NXwsP+PRPBFrG1ZAzlPBXI9x575HTzLgq0/MOjHNH9+3B0s2V5KZ04eevbw3qfG6Kf0H7JDTNQfswokyndoozF5wMSe3EgRPfAnDy/EWcYnRISfRN7MCgXl3586pSfjAshSsHJLFkUyVjMhI9Jge7TfG9gb0Y068H4wYkNTAfFFTU8OdPSvjV9bnMvnoA/ZM6kd6jYwNTg92myO+fyJh+iXyxv5aUru15ZcN+nlyxmxHp3dlbdYrN5TXcOjLNY/oxSOnanoQ4Gz8c0Zv+SZ08I9+EOBuPL9/lMc9sO3CC//yomM4JcYzq24PN5TWkdG3PF/u/iz4yzClj+vXgquxk+vX87n4j0ruzp+oUm8pq2Hf4NO3b2RnZpwdVJ8838BcYym94enfsNpdpbHh6N+LsipKj3zA2M5E3Cw/xVc05OibYeP6zMp/+CbtNkZeR6OlYzcmQh06eD9jp1juET/YcYUSf7p5zbErRJ7EjFx3SwERm3Ld9OztzV5V6kvrMJjBfNMXP0ZI+iab6Ulpbm7FEpJ9f+zCiTKc4+NYB1Wcv+jy+uaKGL6tOY7Mrvq49z/Sx/Ujt3oEfDEtFgNRuHSioaNzxmv9xDOVg7oB8dTaGueeJ5btIaGdnzsoSfn5dFoVf1bKs4AC/uj6HW0e5OlVzp1xYWcuTK3Z7MrMNu/fkoam0b2fjiv49eW5VKdPH9WdQalfuHt+fL/bXNnA6j3L7IAJlnjsFnnL7Gb48fIbNFTUkxNl4YsWuBj4Mc2bvqL49KKioQZzCf39cTG5KZ7YeOIknveGik9/dPozbRqW53nd5DcmdE3hlw36PwnE4hY2lx/n2ohObDb6tc/LUO3toZ2/oozDjL7vYV4fd8H3ZPYrXih8jVFrSJxENh/ml7qSP9POHojAimvWilJqslCpWSpUppR71cTxBKfWm+/jnSqkM07HH3PuLlVI3RlLOcHOyruHn8f0TGZ/pMhkpd/E5B7C+zGWeWfb5V8xaUsi7Ow9z9+IvmLe2jHsWf8FvVuxi1pJCT12ptcXV3L34C+avKQcIWjjPnB39wMRssi/rwgMTs1Eo3iw8yOX9e5Dbq4vnfHOtJe9V8hxOYfehU3xeUcPzn5UxLK0bC6fnc4XJpj42syfzpoxBBESEPe66U4bN3MijOF/n8JhQjJLnf7htGL++IQcbIEgjH8a0KzO4I78Pd17Rj5fWlDHtlc8prKzhgkPYeegMDqdr4XiAzQdO8uXh08B3PoHfv7uHP364j6fe3u0x4cxcUsi9y7bwt6JDzFlVws0jUnnmk2IWriv31MD6bO9RHvnbDj7be5RpV2bw2E0DG/lrzP4Fw7cCLpPd1gMnmLOyhKfe3s35OoffJLrmEK4iilaIxqqC4WwzVutnhVpbLVpETGEopezAC8BNwGDgDqXUYK/TZgAnRCQL+AvwJ/e1g4GfAEOAycCL7vu1Sjbtr2VjxUkE8PV/uqm8hnqnsHzbQRwCn+w5Qr0I/6/oIDeP6I24i9p9vOcwdQ7hmZXFjbKOfVFQUcPcVaX8/Losvqr5hp+9voU5K0sZktqVO/L7sKn0OPe9vsVz/d3j+3s6RO9lXQ1lsu/IaRZMy2Oce32LJZtc+3/79i7A5Xx/YXUZ/zA6nTkrS1i0vsLzJTDu8ft39zSojmtU55199QBeuftyZk0YwMt3uZSR8SUqrKxl+bZD/Px/ivjzJyU4nPDW1oOeZ7UplxI2mLuq1FP2fd6UMfz25iHckd+Hd3ZUeZTVL6/PIc4G/5SXxryfjsHhhIsO4S+flrCp7DgL1pUzc2kRb245yKzXiiiqrGVoWrdGJiVzxrTHab2unHUl1eRnJHLryDRP9vfcVaU8NCkn7F/+llo/IVLZ4YE6zHC2Ge73FC4FFEiucDy/Wc7mFG6M5AzjcqBMRCpEpA74K3Cr1zm3Akvcv78FTFRKKff+v4rIBRHZD5S579cq6da+oa7LTWpPp3YKBfxoVG8euDaL8VmJzP3JaK7NSaL48BnXKB34+7ZDiIJfXJfF9YNTePjGHBZMyeP8hXre2X6IOT8ehdM9Gvb+5zVGJnal+GvhQW4d1Zv5U8agbIrf3jyEh67P4Rffy6a+3onDKQ3KLZhnGA6nkNurC49OzuWeqzI9TuSNZcfJuawLdhus2PZdR/zSlDFMHpLCg6Z1xevqnVysd/LjvHSe+MHggCMm4wtiJCIWVNSggItOJ2tLXErm2twk3v3F1YxI74od+FFeuuf6vPTOnhX+jAq6HeLt/Nvtwzzt2m2KGVdl8usbcvn324ZTVvMNf9viWh73l9fnIMCzn5ZgtyuuzU5i4ZSGKwb66yjGZvbkwYnZ/OWTYu55tZBN5cf53a1DeXRyLtcP6sULd45mcGrXsPxfebcbyig01kbaRqi4eYARClafJ9yj9XApoEjPIsxy+qrabJVIhtWmAV+bPh8ErvB3jojUK6VOAT3d+wu8rk2LnKiR5eS3DpI62jl+zjUGLj7ucoQP6d2Fo2cv8Na2KgBsysbIPt2pOXuBPVVnABibmcjOr08yd1UpNpti0fR8FDB98Rc4BY6e+ZaCiloWTHUFOcxeWsQvJ2Vzz1WZniik0X17sL/mG64fmMLeI6d5/rMyHpyYzXOrynCKIAKL3Mu+GtE+xgxjaFo3AH72+hYemJjtUQqLN+5n7qpS/vJPIxjSuyvXDbyMfLcz2W5TzFxayC0j03jhztGM7tuDe5cVsbq4mnZ2xS0j0zwmqvUl1Z6IK++QUrPSstsU8386hmVfVDJ1bH++N/AyHE5h8tBUfj0plyuzkkjtnMCc1eXMvDrHc3/vNb6N/RvLjuN0iucZ7x7fH+f/b+/MY6yq7jj++fqAAWYQmKKgtSogLmjc21hxQetWm2rrErHW4BZrXdLG2CWhMYbGVrtoakxrXIhxCaK2NVQxFivGio44VRYXlEXbTFkGxUFmcMaZN7/+cc7DO895M3dg3lZ/n+SFe8/6vb97mHPPPef+jhlTJ+zKcXHF2L0zv/o5lyLJ5bI5lyGXTpvYY+XYFcdPYu2mVuY1NvH2+o+ZfsDuTN1zNFc+2MiPvjGFO59bPejuQQa6VLbQ8t1CFHslVm6peO5+9OYCvq/6017Pziwp7k3DYP2hH2xd+SR15l4zXzptIp92dZOpqx+ftq6iLauVdD5wupldEc8vBr5mZtcl0rwZ0zTF8zWEkcRs4GUzeyiG3wcsMLM/91LPlcCVAGSGHDVst32Lcj0Dxuj5Ej6eW3e2q7u9dUhmxGjADEmYmXV9uk1Dho1EcdbSzLo72z/eZdiI0SRuUmfLhtUAQ8dMmNzd2b61e9uWTUNGj5/U2bJhTQzfD7Bs20frMrVj9uxs2bg6Ea5cXLZ1c7NqakftMnT4yEztmD1C+rFfRlK2dXNTtnXzRtXUjrKOthrgg0xd/fhM7dg9AXJlZ9ta1mnYiFE5jZ0tG1ZZR9tWgCGjx++9y4hRu2VbNzcBZOrq9+r+9JMt3W0tzd0dbR8DqKZ2VNRFZ8uG1dbRtjXWufWz+PH7dbZsXG0dbVszdfXjM3X1e+X05Z1nM3X1mWR8oVuTKHcNYLn6BkrOJtF2eyavA1Cmrn73bOvmZsBUU7vr0DETJne1bFjTne3cla5Pm/ouvfgkbV2AccAHubTJe1FqTSnqH6ea2o5iaxsEG2y3aSXoytTVj7ds117ZbVtSPQUUc4TRBHwlcb4XsK5AmiZJQ4DRwOaUeQEws7uBuwEkNXasX5VqPXE5kdTYtaW54nVC0Jp2jXY5kdTYtfXDitcJ1WXTatAJQWt3e2vFa60mm/ZGMecwXgWmSJooaRhhEnt+Xpr5wMx4fB7wnIWn6fnAjLiKaiIwBVhSRK2O4zhOPxRthBHnJK4FniGseJxjZm9Kmg00mtl84D7gQUmrCSOLGTHvm5IeBd4CuoBrzCzba0WO4zhOSSiqLykzWwAsyAu7MXHcDpxfIO/NwM0DrDLVxycVQLXohOrRWi06oXq0VotOqB6t1aKzV/6vfEk5juM4xcP3N3Qcx3FSURUdRjW5GEmh9XpJb0laLukfkvZJxGUlLY2//AUCpdZ5iaRNCT1XJOJmSloVfzPz85ZB6+0Jne9KaknEldKmcyQ1S3qjQLwk3RGvY7mkIxNxJbNpCp0XRX3LJb0k6bBE3PuSVkR7NhZTZ0qt0yVtSdzjGxNxfbabEuv8SULjG7Fd1se4ktp0pzCziv4RJszXAJOAYcAyYGpemquBu+LxDGBePJ4a09cAE2M5mTJrPQkYGY9/mNMaz1sryKaXAHf2krceWBv/HRuPx5ZTa1766wgLLEpq01jXCcCRwBsF4s8EniZ8D3MM8EqZbNqfzmNz9RNc+7ySiHsfGFdBNp0OPLmz7abYOvPSfpuwIrQsNt2ZXzWMMKrJxUi/Ws1skZlti6cNhG9MSk0amxbidGChmW02s4+AhQR/X8VioFovBOYWUU9BzOwFwmq/QpwNPGCBBmCMpD0osU3702lmL0UdUL42mtPSn00LsTNtfMAMUGfZ2ujOUg0dRm8uRvLdhPRwMQIkXYz0l3cwGWh9lxOeOHMMl9QoqUHSd4ohMJJW57nxtcTjknIfUlasTePrvYnAc4ngUtk0DYWupdQ2HQj5bdSAv0v6l4KXhUrg65KWSXpa0sExrCJtKmkk4WEg6bWiEm3aK9WwRWtvn6znL+0qlCZN3sEkdX2Svg8cDZyYCN7bzNZJmgQ8J2mFma0pk86/AXPNrEPSVYQR3Mkp8w4mA6lvBvC49fxmp1Q2TUOltNNUSDqJ0GEclwieFu25O7BQ0sr4dF0uXgP2MbNWSWcCTxA+9K1ImxJeRy02s+RopNJsWpBqGGEMxMUI2kEXI4NEqvoknQLMAs4ys45cuJmti/+uBZ4HjiiXTjP7MKHtHuCotHkHmYHUN4O8oX4JbZqGQtdSapv2i6RDgXuBs81suyvWhD2bgb9SZi/SZvaxmbXG4wXAUEnjqECbRvpqoxVh0z4p9yRKfz/CKGgt4VVDbvLq4Lw019Bz0vvReHwwPSe911LcSe80Wo8gTMZNyQsfC9TE43HAKoo0SZdS5x6J4+8CDfG4Hngv6h0bj+vLadOY7gDC5KHKYdNEnftSeIL2W/Sc9F5SDpum0Lk3Yb7v2LzwWmBU4vgl4Ixi6kyhdULunhP+0P4n2jdVuymVzhife5CtLbdNd/gayy0g5Y04E3g3/qGdFcNmE57QAYYDj8VGvgSYlMg7K+Z7B/hmBWh9FtgILI2/+TH8WGBFbNgrgMvLrPPXwJtRzyLgwETey6KtVwOXltum8fwm4Ja8fKW26VxgPdBJeMK9HLgKuCrGi7Cp2Jqo5+hy2DSFznuBjxJttDGGT4q2XBbbxqwS3Pv+tF6baKcNJDq53tpNuXTGNJcQFuEk85Xcpjvz8y+9HcdxnFRUwxyG4ziOUwF4h+E4juOkwjsMx3EcJxXeYTiO4zip8A7DcRzHSYV3GE5VIckkPZg4HxK96j5ZTl39Iel5SZ/byzmGvxNdW7wq6fABlju9r2vvK17SvZKmxuPzJb0taZGkw+NX047TA+8wnGqjDThE0oh4firw33IIiV4FBoOLzOww4I/AbwepzH4xsyvM7K14ejlwtZmdBBxO+IbBcXrgHYZTjTxN+Goa8jx/SqqNexO8Kul1SWfH8H0l/VPSa/F3bAzfQ9ILiX0Kjo/hrYkyz5N0fzy+X9JtkhYBt/ZR3whJj0TnjfOAXAfXFy+TcJAn6TRJL0e9j0mqi+FnSFop6UXgnET6ExN7LrwuaVSMqosOJFdKejh6ct4+6lHYQ+I44C5JtxM+irwglnNBulvifBGoBueDjpPPI8CN8VXLocAc4PgYN4uw18BlksYASyQ9CzQDp5pZu6QphE7maOB7wDNmdrOkDDAyRf37A6eYWVbSrwrU9wNgm5kdGv0yvZai3DMIzvOI/pB+Eetpk/Qz4HpJvyH49jqZ8FX4vET+G4BrzGxx7FzaY/gRBDc564DFwDTgxVwmM5st6WTgBjNrlLSM8BX6tSk0O18gvMNwqg4zW66wq+KFwIK86NOAsyTdEM+HE3wjrQPujHMEWcIffYBXgTmShgJPmNnSFBIes8884haq7wTgjoTe5X2U97CkWsKmP7ld+I4hbAC2OA4IhhFGIAcC75nZKgBJDwE5l9iLgdskPQz8xcyaYt4lZtYU0y8l+Dza3mE4Tlq8w3CqlfnA7wg7rn0pES7gXDN7J5lY0k0EH16HEV7FtkPY+EbSCYRXXA9K+q2ZPUBPV9jD8+puS1EfpHenfRHBl9AtBF9T58RyF5rZhXnlHl6oXDO7RdJThPmHhugVGaAjkSyL/793dhCfw3CqlTnAbDNbkRf+DHBd4j19zp35aGC9mXUDFxOe5nObLjWb2T3AfXz2hL9R0kGSdiF46y1EofpeIHQESDqE8OqsIGbWSXgFdYykgwiO9KZJ2i+WMVLS/sBKYKKkyTHr9g5F0mQzW2FmtwKNhNHIjrAVGNVvKucLh3cYTlViZk1m9odeon4JDAWWS3ojnkNYgTRTUgPhdVRulDAdWCrpdeBcIFfmz4EnCbv3re9DSqH6/kSYbF4O/JTgRbm/a/oE+D1hLmETwbvp3FhGA8FjcDvhFdRTcdL734kifhwn7pcBn9Bzp7yBsAiY6pPeTj7urdZxHMdJhY8wHMdxnFR4h+E4juOkwjsMx3EcJxXeYTiO4zip8A7DcRzHSYV3GI7jOE4qvMNwHMdxUuEdhuM4jpOK/wGTUAU0UWIGdAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from matplotlib import pyplot as plt\n",
    "from sklearn.model_selection import KFold\n",
    "from sklearn.tree import DecisionTreeRegressor\n",
    "\n",
    "# paste your get_features_targets function here\n",
    "def get_features_targets(data):\n",
    "  features = np.zeros((data.shape[0], 4))\n",
    "  features[:, 0] = data['u'] - data['g']\n",
    "  features[:, 1] = data['g'] - data['r']\n",
    "  features[:, 2] = data['r'] - data['i']\n",
    "  features[:, 3] = data['i'] - data['z']\n",
    "  targets = data['redshift']\n",
    "  return features, targets\n",
    "\n",
    "# paste your median_diff function here\n",
    "def median_diff(predicted, actual):\n",
    "  return np.median(np.abs(predicted - actual))\n",
    "\n",
    "# complete this function\n",
    "def cross_validate_predictions(model, features, targets, k):\n",
    "  kf = KFold(n_splits=k, shuffle=True)\n",
    "\n",
    "  # declare an array for predicted redshifts from each iteration\n",
    "  all_predictions = np.zeros_like(targets)\n",
    "\n",
    "  for train_indices, test_indices in kf.split(features):\n",
    "    # split the data into training and testing\n",
    "    train_features, test_features = features[train_indices], features[test_indices]\n",
    "    train_targets, test_targets = targets[train_indices], targets[test_indices]\n",
    "    \n",
    "    # fit the model for the current set\n",
    "    model.fit(train_features, train_targets)\n",
    "        \n",
    "    # predict using the model\n",
    "    predictions = model.predict(test_features)\n",
    "        \n",
    "    # put the predicted values in the all_predictions array defined above\n",
    "    all_predictions[test_indices] = predictions\n",
    "\n",
    "  # return the predictions\n",
    "  return all_predictions\n",
    "\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "  data = pd.read_csv('opticaldatafinals_SDSS.csv')\n",
    "  features, targets = get_features_targets(data)\n",
    "\n",
    "  # initialize model\n",
    "  dtr = DecisionTreeRegressor(max_depth=19)\n",
    "\n",
    "  # call your cross validation function\n",
    "  predictions = cross_validate_predictions(dtr, features, targets, 10)\n",
    "\n",
    "  # calculate and print the rmsd as a sanity check\n",
    "  diffs = median_diff(predictions, targets)\n",
    "  print('Median difference: {:.3f}'.format(diffs))\n",
    "\n",
    "  # plot the results to see how well our model looks\n",
    "  plt.scatter(targets, predictions, s=0.4)\n",
    "  plt.xlim((0, targets.max()))\n",
    "  plt.ylim((0, predictions.max()))\n",
    "  plt.xlabel('Measured Redshift')\n",
    "  plt.ylabel('Predicted Redshift')\n",
    "  plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1><center>Splitting the train and test sets</center></h1>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number data galaxies: 370082\n",
      "Train fraction: 0.7\n",
      "Number of galaxies in training set: 111024\n",
      "Number of galaxies in testing set: 259058\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "data = pd.read_csv('opticaldatafinals_SDSS.csv',delimiter=\",\")\n",
    "\n",
    "\n",
    "fraction_training = 0.7\n",
    "  # split the data using your function\n",
    "training,testing = train_test_split(data, test_size = 0.7, random_state = 0)\n",
    "\n",
    "  # print the key values\n",
    "print('Number data galaxies:', len(data))\n",
    "print('Train fraction:', fraction_training)\n",
    "print('Number of galaxies in training set:', len(training))\n",
    "print('Number of galaxies in testing set:', len(testing))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1><center>Generating features and targets</center></h1>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Features shape: (370082, 13)\n",
      "Targets shape: (370082,)\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "def generate_features_targets(data):\n",
    "  # complete the function by calculating the concentrations\n",
    "\n",
    "  targets = data['class']\n",
    "\n",
    "  features = np.empty(shape=(len(data), 13))\n",
    "  features[:, 0] = data['u-g']\n",
    "  features[:, 1] = data['g-r']\n",
    "  features[:, 2] = data['r-i']\n",
    "  features[:, 3] = data['i-z']\n",
    "  features[:, 4] = data['ellipticity']\n",
    "  features[:, 5] = data['mCr4_u']\n",
    "  features[:, 6] = data['mCr4_g']\n",
    "  features[:, 7] = data['mCr4_r']\n",
    "  features[:, 8] = data['mCr4_i']\n",
    "  features[:, 9] = data['mCr4_z']\n",
    "\n",
    "  # fill the remaining 3 columns with concentrations in the u, r and z filters\n",
    "  # concentration in u filter\n",
    "  features[:, 10] = data['petroR50_u']/data['petroR90_u']\n",
    "  # concentration in r filter\n",
    "  features[:, 11] = data['petroR50_r']/data['petroR90_r']\n",
    "  # concentration in z filter\n",
    "  features[:, 12] = data['petroR50_z']/data['petroR90_z']\n",
    "\n",
    "  return features, targets\n",
    "\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "  data = pd.read_csv('opticaldatafinals_SDSS.csv')\n",
    "\n",
    "  features, targets = generate_features_targets(data)\n",
    "\n",
    "  # Print the shape of each array to check the arrays are the correct dimensions. \n",
    "  print(\"Features shape:\", features.shape)\n",
    "  print(\"Targets shape:\", targets.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1><center>Train the decision tree classifier</center></h1>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Some initial results...\n",
      "   predicted,  actual\n",
      "0. E2, E4\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import math\n",
    "from sklearn.tree import DecisionTreeClassifier\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "def generate_features_targets(data):\n",
    "    targets = data['class1']\n",
    "    \n",
    "    features = np.empty(shape=(len(data), 13))\n",
    "    features[:, 0] = data['u-g']\n",
    "    features[:, 1] = data['g-r']\n",
    "    features[:, 2] = data['r-i']\n",
    "    features[:, 3] = data['i-z']\n",
    "    features[:, 4] = data['ellipticity'] \n",
    "    features[:, 5] = data['mCr4_u']\n",
    "    features[:, 6] = data['mCr4_g']\n",
    "    features[:, 7] = data['mCr4_r']\n",
    "    features[:, 8] = data['mCr4_i']\n",
    "    features[:, 9] = data['mCr4_z']\n",
    "    features[:, 10] = data['petroR50_u']/data['petroR90_u']\n",
    "    features[:, 11] = data['petroR50_r']/data['petroR90_r']\n",
    "    features[:, 12] = data['petroR50_z']/data['petroR90_z']\n",
    "    return features, targets\n",
    "\n",
    "# complete this function by splitting the data set and training a decision tree classifier\n",
    "def dtc_predict_actual(data):\n",
    "    training,testing = train_test_split(data, test_size = 0.7)\n",
    "    features_training, targets_training = generate_features_targets(training)\n",
    "    features_testing, targets_testing = generate_features_targets(testing)\n",
    "    dtc = DecisionTreeClassifier()\n",
    "    dtc.fit(features_training, targets_training)\n",
    "    predictions = dtc.predict(features_testing)\n",
    "    return predictions, targets_testing\n",
    "\n",
    "\n",
    "if __name__ == '__main__':\n",
    "    data = pd.read_csv('opticaldatafinals_SDSS.csv')  \n",
    "    predicted_class, actual_class = dtc_predict_actual(data)\n",
    "\n",
    "  # Print some of the initial results\n",
    "    print(\"Some initial results...\\n   predicted,  actual\")\n",
    "    for i in range(1):\n",
    "        print(\"{}. {}, {}\".format(i, predicted_class[i], actual_class[i]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center><h1> Accuracy in classification</h1></center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Our accuracy score: 0.9999891915845678\n",
      "Confusion matrix, without normalization\n",
      "[[  5943      0      0      0      0      0      0      0      0      0]\n",
      " [     0  33985      0      0      0      0      0      0      0      0]\n",
      " [     0      0    902      0      0      0      0      0      0      0]\n",
      " [     0      0      0  82826      0      0      0      0      0      0]\n",
      " [     0      0      0      1 100339      0      0      0      1      0]\n",
      " [     0      0      0      0      0    124      0      0      0      0]\n",
      " [     0      0      0      0      0      1     73      0      0      0]\n",
      " [     0      0      0      0      0      0      0  51765      1      0]\n",
      " [     0      0      0      0      0      0      0      0  74823      0]\n",
      " [     0      0      0      0      0      0      0      0      0  19298]]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVUAAAEYCAYAAADsymWcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzsnXeYVEXWh98zM0RBcg5KRrJkURAlKBjQNYFhTajryvqZF3UVMYc1YcbsGjCsChJERHGNRCUqAoJIkpzTMHO+P6pmphlmpnt6eqbTeXnqoW/duqfO7e45XVW3qn6iqhiGYRiRISXaDhiGYSQSFlQNwzAiiAVVwzCMCGJB1TAMI4JYUDUMw4ggFlQNwzAiiAVV4yBEpJyIfCIi20Tk/SLYuUBEPoukb9FCRHqKyOJo+2HEB2LzVOMTETkfuAFoCewAfgLuU9Vvimj3IuAfQA9VPVBkR2McEVGgmaoujbYvRmJgLdU4RERuAJ4A7gdqAQ2BZ4FBETB/BPBrMgTUUBCRtGj7YMQZqmopjhJQCdgJnFNAmTK4oLvGpyeAMv5cb2AVcCOwHlgLXOrPjQT2A+m+jsuBu4A3A2wfCSiQ5o8vAX7DtZaXAxcE5H8TcF0PYCawzf/fI+DcNOAe4Ftv5zOgej73luX/LQH+nwEMBH4FNgO3BZTvCnwPbPVlnwZK+3P/8/eyy9/veQH2/wmsA/6TleevaeLr6OiP6wIbgd7R/m5Yio1kLdX44xigLPBRAWVuB7oDHYD2uMDyr4DztXHBuR4ucD4jIlVUdQSu9fuuqlZQ1ZcLckREDgNGAQNUtSIucP6UR7mqwARfthrwGDBBRKoFFDsfuBSoCZQGbiqg6tq496AecCfwInAh0AnoCdwpIo192QzgeqA67r3rA/wdQFV7+TLt/f2+G2C/Kq7VfmVgxaq6DBdw3xKR8sCrwGuqOq0Af40kwoJq/FEN2KgFd88vAO5W1fWqugHXAr0o4Hy6P5+uqhNxrbQWYfqTCbQRkXKqulZVF+ZR5hRgiar+R1UPqOo7wC/AaQFlXlXVX1V1D/Ae7gchP9Jx48fpwBhcwHxSVXf4+hcC7QBUdbaq/uDrXQG8ABwfwj2NUNV93p+DUNUXgSXAdKAO7kfMMAALqvHIJqB6kLG+usDvAce/+7xsG7mC8m6gQmEdUdVduC7z34C1IjJBRFqG4E+WT/UCjtcVwp9NqprhX2cFvT8Dzu/Jul5EmovIeBFZJyLbcS3x6gXYBtigqnuDlHkRaAM8par7gpQ1kggLqvHH98Be3DhifqzBdV2zaOjzwmEXUD7guHbgSVWdrKr9cC22X3DBJpg/WT6tDtOnwvAczq9mqno4cBsgQa4pcEqMiFTAjVO/DNzlhzcMA7CgGneo6jbcOOIzInKGiJQXkVIiMkBEHvbF3gH+JSI1RKS6L/9mmFX+BPQSkYYiUgm4NeuEiNQSkdP92Oo+3DBCRh42JgLNReR8EUkTkfOAVsD4MH0qDBWB7cBO34q+Otf5P4HGh1xVME8Cs1V1KG6s+Pkie2kkDBZU4xBVfQw3R/VfwAbgD2AY8LEvci8wC5gHzAfm+Lxw6poCvOttzebgQJiCm0WwBvdE/Hj8Q6BcNjYBp/qym3BP7k9V1Y3h+FRIbsI9BNuBa0W/m+v8XcDrIrJVRM4NZkxEBgEn44Y8wH0OHUXkgoh5bMQ1NvnfMAwjglhL1TAMI4JYUDUMw4ggFlQNw4h5ROQVEVkvIgsC8qqKyBQRWeL/r+LzRURGichSEZknIh0DrrnYl18iIhcH5HcSkfn+mlEiIgXVURAWVA3DiAdewz0gDGQ4MFVVmwFT/THAAKCZT1fiptVlrewbAXTDrTIcERAkn/Nls647OUgd+ZJUD6qqVa+uDRrmni4ZPqkSbLqjYcQPv/++go0bN0b0S516+BGqBw5ZlHYIumfDZFXNHTQPQkSOBMaraht/vBi358JaEakDTFPVFiLygn/9TmC5rKSqV/n8F3D7TkwDvlTVlj5/SFa5/OooyM+k2oGnQcMj+Px/0yNmr0LZpHr7jATn2G6dI25TD+yhTIugM9XY+9MzLUVkVkDWaFUdHeSyWqq6FsAHvZo+vx5ummEWq3xeQfmr8sgvqI58sahgGEbxIQIpqaGU3KiqkYrqebW2NYz8sEiaMVURWbH450X07tGJvr26AbBg/lwGnHgcvbp14IJzzmDH9u0HXbPqj5UcUbsyzzz5GAB79+6lf+9j6H1MR47r0p57Ro4osM7PJn9Ku9YtaN2yKY88/GCR7yHS9orDptlLbHthISnBU3j86bvk+P/X+/xVQIOAcvVxC1QKyq+fR35BdeRPtPceLKkErGjdtp1u2JGenTp07KQfT5qqG3ak6xPPjNYbbrntoPOnnH6mnnbGWXrXvQ/phh3pun77fl2+dotu2JGuazbv1s5duuq0r7/XPel6SNq594A2atxYFy1eptt27dO2bdvpnLkL8ywbSoq0vXjw0eyVrL2OHTtppP/upHxNLdv5+qAJmBXC3/CRwIKA40eA4f71cOBh//oUYBKuBdodmOHzq+L2/K3i03Kgqj8305cVf+3AguooKCVNSzUvli75lR7H9gSg94l9GT82Z4vSiZ+M5cgjG9HyqFbZeSJChQpu86T09HQOpKcj+TysmjljBk2aNKVR48aULl2ac84bzPhPxobta6TtxYOPZi+27IWHRKSlKiLv4DYTaiEiq0TkcuBBoJ+ILAH6+WNwe038BizFLU3O2j93M24z9Jk+3e3zwO0J8ZK/ZhkusFJAHfmSTEFVf1u6hD49u/LGK24jpaOOas2nEz4BYNxHH7B6tRvD3rVrF089/gg33XrHIUYyMjLo3aMTRzWuy4l9+9G1W7c8K1uzZjX16+f0NOrVq8/q1eFvyhRpe/Hgo9mLLXthIbgx1WApCKo6RFXrqGopVa2vqi+r6iZV7aOqzfz/m31ZVdVrVLWJqrZV1VkBdl5R1aY+vRqQP0tV2/hrhqlvmuZXR0HEfFAVkQwR+SkgDff5IiL3icivIvKziFwbxNSxzVsexZgPx/PKi8/x3Tdf8+SzL/LKi8/Rp2dXdu7cSelSpQF4+L6RXDXs/7JbpYGkpqYy7bvZzPtlBbNmzmDhggWHlAGyuie576WQd1989orDptlLbHvhIe5hVbCUQMTD0/89qprXLvCX4AadW6pqZrCpDqq6pkPHTtSoUZOBp53Bj7Nncs3/3cD7Y10rf9mSX5kyeSIAs2fN4JOxH3L3HbeybdtWUlJSKFO2DEOvuibbXqXKlel1fG8+++xTWrdpc0h99erVZ9WqnNkbq1evom7duoeUC5VI24sHH81ebNkLm/AfRMUl8Xy3V+PGRDIBVDXfp3IicpiIVATXtZ82dQotW7VmwwZ3SWZmJo89cj8XX+bkiMZ/No05C5cyZ+FSrvr7tVx343CGXnUNGzdsYNvWrQDs2bOHL6Z+TosWeW10D527dGHp0iWsWL6c/fv38/67Yzjl1NPDvtlI24sHH81ebNkLG2upxhzlRCRQTO4BdQJtTYDzRORM3J6i16rqktwXi8iVuL1Gmy5aMJ+TevfgL+cOpk+/k3jh2VG8MtrtL3zK6Wdw/kWXFOjIn3+uZdhVl5GZkUFmpnLOuecy8JRT8yyblpbG408+zWmnnERGRgYXX3IZrVq3DuP2i8dePPho9mLLXliEPk81YYj5ZaoislNVDxncFJGdOHG2R0XkL8D1qtqzIFsdOnZSW1FlGHlzbLfOzJ49K6LNxpSKdbXM0VcGLbf365GzNXKT/6NKPHf/VwH/9a8/wqtnGoYRS0RmSlU8Ec938zFwon99PPBrFH0xDCMvBEhNDZ4SiHjov+YeU/1UVYfjJuG+JSLX4wTnhkbFO8MwCibBHkQFI+aDqqrm+TOmqltxy9EMw4hZJOG698GI+aBqGEacYy1VwzCMCJGEU6osqBqGUbxY998wDCOCWPc/cUkVieiE/U8WrAleqBCc1iYK67INo1ix7r9hGEbkEKz7bxiGETmSb0pVct1tEArS89m7dy/HHdOVrh3b07F962x9qufuupGbz+3LTef25dGbrmDv7l0AbFiziruvOpebzu3LXUPPZtOfOUMFbz5xLzecdQLX/+V4Xnnojux9L/v36U271i3o1qkD3Tp1YP36gzfeMo0qsxdte2GRZLtURV07qiRTx46dwtbz2b0/Uzds2aF70lW3796vnbt01XtfH6evff2Lvvfjan3vx9V6ygVX6Pn/uFXf+3G1du97iv595OP63o+r9c4X3tWep/xF3/txtd7z6sfaon1nHTNrpY6ZtVKbte2oI158X9/7cbX27HW8fvP9TNOoMnuJo1FVqaGWHfRC0EQIGlXxkqyl6gmm55OfPlX5ChUB9+O0f9/e7F/dVb8toW234wBo3eVYZk37LNvO/v37OJC+n/T9+8k4cIBKVWsU2b/iuGezZ/aKjNiGKklLKHo+GRkZdOvUgYZ1a3Ji3340a9sRgGdHXM+VfTuwZsVSBgy+DIAjmrdi+lSnJDDji0ns2bWTHVs307x9Z1p37sGV/TpyZf+jad/jeOo3bpZdx1VDL6Vbpw48cN89qOZsy2gaVWYv2vbCJsm6/zEfVPPTqAo4/5TfW7VIBAawANsHHaempjJ99k8sXbGKWTNnsHLpLwD8feTjvPDZHOo1asZ3n40D4KLr72DR7B+4ZXB/Fs3+gao1a5Oamsa6lctZvXwJz0+exQuTZ7Ngxrcsmv0DAK++8RazfprP59O+5ttvvubtN/9TKP+K457NntkrKiISNCUSMR9U8RpVASl7tF1EOgOVI1FJYfR8Knt9qp++m5adl5KaSo/+pzN96gQAqtaszU2PvsTDYz5jyLB/AlC+4uHM+PJTmrXtSNnyh1G2/GEcfeyJLJk/x/tQD4CKFSty3uDzmTlzRlj+Fcc9mz2zFw6u9y9BUyIRD0E1T0QkFXgEuCUS9oLp+WzYsIGtufSp6h7RmHUrlwOuVTDrf1Ooe2RTALZv2UxmZiYAH73yFCcMGgxA9dp1+Xn2D2QcOMCB9HQWzfmeeo2aknHgABs3bgTcmO3EieNp3bpNyP4Vxz2bPbNXdIK3UhOtpRoP81Tz06gaBoxT1bWR+FCC6fmsW7uWKy67mIyMDDI1k7POPpfWPfsy4rIz2b1rJ6hyRPNWDL3tAQAWzfqOt596ABHhqI7dufzW+wDo3vdUFsz8lpvO7QMIHXr0pvPx/dm7ZzenDzyT9PR0MjIzOOHEvlw29IqQ/SuOezZ7Zi8SJFrQDEZcalSJSF3gPaC3qh7IT8fKl70SuBKgQcOGnX5d9nvEfLNlqkYiURwaValVG2mFk+4OWm77mL+aRlWUORpoCiwVkRVAeRFZmldBVR2tqp1VtXON6sGnLhmGEUEkxJRAxEP3/xBUdQJQO+vYt1SbRtElwzDyQEi8MdNgxENQzU+jyjCMOMCCaoyh+WhU5SqT53iqYRjRJyUlXkcZwyPmg6phGHFMAo6ZBsOCqmEYxUqydf+Tq11uGEaJIhGc/C8i14vIQhFZICLviEhZEWkkItNFZImIvCsipX3ZMv54qT9/ZICdW33+YhE5KSD/ZJ+3NPdy+MJgQdUwjGIlEstURaQecC3QWVXbAKnAYOAh4HFVbQZsAS73l1wObPGzgh735RCRVv661sDJwLMikupXaD4DDABaAUN82UJj3f8iEOnJ+pmZkV2IkZJga6qNOEQi2v1Pw80GSgfKA2uBE4Hz/fnXgbuA54BB/jXAB8DT4hwZBIxR1X3Acj+/vasvt1RVfwMQkTG+7KLCOmktVcMwipUQu//VRWRWQLoy0Iaqrgb+DazEBdNtwGxgq6oe8MVWAfX863rAH/7aA758tcD8XNfkl19orKVqGEaxIUioU6o2FrRMVUSq4FqOjYCtwPu4rnpusrp7eTWPtYD8vJwMq+toLdUAYkUf6JmnnqTz0W3p3KENT496AoDNmzdz6oD+tGvVnFMH9GfLli0AjHnnLbp2ak/XTu058fhjmTdvbon4aPaS015YRGaZal9guapuUNV04EOgB1BZRLIah/WBrA05VgENAPz5SsDmwPxc1+SXX3iiredSkqkoGlUloTe0a1+mzpgzT49q1Vo3bNmp23bt194n9NG5CxbrdTfcpCPvvV937cvUkffer9ffeLPu2pepU6d9o6vWbdJd+zL1w7ETtHOXrrprX6bu2pcZ85pIZi+27BWHRlWpGk207t8+DJoIolEFdAMW4sZSBTd++g9ci3WwL/M88Hf/+hrgef96MPCef90amAuUwbV6f8M99ErzrxsBpX2Z1uHcs7VUPbGiD7T4l5/p2q0b5cuXJy0tjZ69ejFu7EdM+GQcF1x4MQAXXHgx48c5W92P6UGVKlUA6NqtO6tXryp2H81ectoLl0hMqVLV6bgHTnOA+bhe9mjgn8AN/oFTNeBlf8nLQDWffwMw3NtZiNvhbhHwKXCNqmb4cddhwGTgZ1wQXhjO/VpQ9cSKPlCrVm349uuv2bRpE7t372byp5NYveoP1q//kzp16gBQp04dNmxYf8i1r7/6Mv1POrnYfTR7yWkvXCK187+qjlDVlqraRlUvUtV9qvqbqnZV1aaqeo5/qo+q7vXHTf353wLs3KeqTVS1hapOCsifqKrN/bn7wr3fmH9QJSIZuF+mLMao6oMi8jVQ0efVBGao6hnh1qMxog/U8qijuOGmWzhtYH8qVKhA27btSE0L/jF9Ne1L3njtFaZ8+XWx+2j2ktNeuCTbiqqYD6p4jarcmaraM+u1iPwXKFK/Jpb0gS6+9HIuvtTNYR5xx23Uq1efmjVrsXbtWurUqcPatWupUaNmdvn58+dxzd+u4KNxE6lWrVqJ+Gj2ks9eOCSiXEow4r77LyIVcROAPy6KnVjSB1q/3nXt/1i5knEff8Q55w1h4Kmn8dabrwPw1puvc8ppp2eXOf/cs3jp1Tdo1rx5iflo9pLPXrikpKQETYlEPLRU89OoyuJMYKqqbi9KJbGkD3TB4LPZvGkTaaVK8diTT1OlShVuvHk4F51/Hm+8+gr1GzTkzXfeA+CB++9m8+ZNXHftNdn1fvP9zGL30ewln72wSa6GanxqVOU6Pwl4SVX/m8/5YtOoijS2TNWIJsWhUVWmVjOtd8GTQcstf/wU06iKBUSkGm7d7oT8yqhpVBlG9JDITKmKJ+Kh+18Q5wDjVXVvtB0xDONQ3DLVxAqawYiHoFqQRtVgIEpr7wzDCIUEa4gGJeaDqhagUaWqvUvQFcMwwiDRuvfBiPmgahhGHCPWUjUMw4gYAqSmJldUtaBqGEaxYt1/wzCMSGHdfyOaRHrqSd3L3o6ovTWvnB+8kGEEUIid/xMGC6qGYRQr1lI1DMOIIMk2pppc7fIgxLo+UDB7qQKlU10qlQJlSqXwwt96MP2hU/n2/oE8NbQbaf5JbMVypXj7+uP5370D+O7+gZzfszEAbRpWZvKd/fnu/oF8fe8AzuzW8KA60lJy6sjroW68v4dmL8L4MdVgKZGwoOrJyMjgumuvYewnk/hx3iLeH/MOPy8qtOR3VO2lpsD+DJcA/tLtCN7/bgXd/jmeY2+bSNlSqVx0fBMAhvZtxuLV2+j1r0mc9sBU7hlyNKVSU9izP4OrX/ieHrdN5Jx/T+O+CzpxePlSzr64KTJZdWTk2v8lEd5DsxdZBPesIFhKJCyoemJdHygUe7m/muu27uHzeTmCkHN+20TdquUBUIUK5dzoz2Fl0tiyaz8HMjNZtm4Hv/25I/v6jdv3Ur1iWcAF7QOZsXXPZq/k7IVLsm2oYkHVE+v6QKHYO5AJZVJdAvhywbrsc2mpwrnHNmLqvLUAvPT5rzSvU4lFo87km/sHcuubs8m9C2THxtUonZbC8vUuyAqQIjnDC7n/FBLhPTR7kce6/zGGiGSIyE8BabjPbyQi00VkiYi8KyKli1JPrOsDhWIvRWBfhksA5/Q4Mvvcvy/uwveL1/PDrxsAOLFtHRas3EKraz/i+H9N4uG/dqZi2ZznlrUqleW5q45h2Is/HBJss7r+pXLtypAI76HZizBJuPVfzAdVvEZVQMoabX8IeFxVmwFbgMuLUkms6wMFs5ciEPgnlKHQtVl1AG45ow3VKpbh9rfnZJ8/v2djPpnl7C1fv5PfN+ykWd1KAFQsm8aYG3tz/wdzmbVsU/Y1CmTto52ph7ZU4/09NHuRJ2vrPxtTjXHE/bSdiNMBB3gdCFtJFWJfHyiYPVUXWLNIFfh1zXYuOr4JJ7atwxXPfndQi3PVpt0c37o2ADUOL0vT2oezYv1OSqWm8Mb/9eLdb5czdmbOHyS4QJpVR+4gHo17Nnslay9ckq37Hw/zVA/RqAKmAltV9YDPWwXUK0olsa4PFMxeViuytO+Sq8LrXy5l1Yvn8sfGXUy+sz8A42f9wSNjF/DvsQt45orufHPfQERg5Hs/sXnnPs7pcSQ9WtSkaoUyDDnOTbO65sXvWbByKwcy3Vhqqv8pTs+I7j2bvZK1Fy6J1r0PRlxqVIlIDeB7VW3qjxsAE1W1bR7Xx41GVaSxZapGYSgOjaqKDVpqh+teClrum5t6mkZVlNkIVBaRrJZ2fWBNXgVNo8owoos9qIoD1DWvvwTO9lkXAyU/Ac8wjKAk25hqPATVcrmmVGU9/f8ncIOILAWqAS9Hz0XDMPIj2VqqMf+gKj+NKlX9DSdPbRhGjCKSeFOmghHzQdUwjPgmwRqiQQna/ReR7iJS3r8eIiIP+6fthmEYQUkRCZoSiVDGVEcDe0SkHXAb8CfwZrF6ZRhGwhCpB1UiUllEPhCRX0TkZxE5RkSqisgUv1x9iohU8WVFREaJyFIRmSciHQPsXOzLLxGRiwPyO4nIfH/NKAlzsDeUoHrAP20fBDypqo8CFcOpzDCM5EIEUlMkaAqRJ4FPVbUl0B74GRgOTPXL1af6Y4ABQDOfrgSec/5IVWAE0A33TGZEViD2Za4MuO7kcO45lDHVXSJyM3Ah0FtEUoBS4VRmlCyRnqxfpcuwiNrbMvPpiNozYpNIPN0XkcOBXsAlAKq6H9gvIoOA3r7Y68A03MygQcAbvkH4g2/l1vFlp6jqZm93CnCyiEwDDlfV733+G7il75MK62soLdXzcHtn/E1V1+Im2j9W2IoMw0hOQuz+VxeRWQHpylxmGgMbgFdF5EcReUlEDgNq+biE/7+mL18PCNy8Imspe0H5q/LILzShtFS3AP9W1UwRaQK0AP4TTmWGYSQXAqSG1lLdGGSZahrQEfiHqk4XkSfJ6ernV3VuNIz8QhNKS/VroKxvOn8FXA28Ek5lhmEkGSFM/A9xeGAVsEpVp/vjD3BB9k8fm/D/rw8oHzhLKWspe0H59fPILzShBNUUVd0NnAU8raqn4QaJE45YF12LtL2rhl5Gw7o16dShDeBE/cqk5ux0lUUpL/Y3/rlhVK5YLjv/0VvOZsHYEcx491Y6tHTfx4Z1qvDtW7fww5jhzP7gdoaefVx2+bFP/53p77r8UbcPzs4XDhYsLM57jmV7uT+PSBB14T8i8/RfVdcBf4hIC5/VB1gEjMMtU4eDl6uPA/7qZwF0B7b54YHJQH8RqeIfUPUHJvtzO/wUUgH+SphL30MKqiLSBTgfGF+I6+KKWBddKw4Rt4suvoSx4z/NqSMzRzQwi7QUt6Xg/gyYNmMxN13qthA86bhWNGlYgzaDRjLs3ncYdZsLkms3bOeESx6j++AH6XXRI9x0aT/q1HCbX1/4z1fodt6DdDr7PmpUqZC9N2upVLeN4P4MV1daSvHcc6zby/15FJWYEf6L3DzVfwBvicg8oANwP/Ag0E9ElgD9/DHAROA3YCnwIvB3AP+A6h5gpk93Zz20wvXCX/LXLCOMh1QQWnC8ARgJTFDVBSLSGDckkFDEuuhacYi4HdezF1WrVs0+zmsAKUVyVFPf/GQ6p53QDoBTj2/H2+NnADBj/goqVSxH7eqHk34gg/3pbpvbMqVLHfQHs2PXXgDS0lIolZbTHJaAujMCNsKO9fcw0vZyfx5FJVaE/yK187+q/uR3nGunqmeo6hZV3aSqfVS1mf9/sy+rqnqNqjZR1baqOivAziuq2tSnVwPyZ6lqG3/NMA1zX9SgQVVVv1DVgap6nz/+TVX/Hk5lsUysi65FS8Qt8Ou+buN2alR1U5Tr1qzMqnVbss+t/nMrdWtWBqB+rcrMePdWlky6h0df+5y1G7Zllxv3zDWsnPogO3fvy5ZmUXICaZYMNsT+exgrwnr5EQv+hdL1T7AFVSEtU60uIg+IyDgR+SwrlYRzvv78hP+G+ZUPKiLVi1pPrIuuxYSI20F1H5qX5eOqP7fS9bwHaDNoJBee1pWaVXPWipx+zTM06ncbZUqnZQfS9AwXTHOP5cb6exhrn0luYsU/W6Z6KG8CK4DmOLG9dcBPBV0QYfIT/vsW6AtEZCv/WBddi5aIW+CfZe3qh7Nhs5OrXv3nVurXrpJ9rl6tyge1SAHWbtjGomXrOLZjk4Py9+0/wPiv5mcHVQXSM3PGVLPqjPX3MBaE9QoiVvyzoHooNVT1BWC/qk7FPWGL+pZ7qvqjqq6IlL1YF12LlohbprpWJMCFp3Vj/LR5AEz4aj7nn+q+Bl3bHsn2nXtYt3E79WpWpmwZt+CucsVyHNOhMb+uWM9h5UpTu/rhAKSmpnDysa0Okb4G95AqI9O9jvX3MFaE9fIjFvxzD6qCp0QilMn/6f7/dSJyEm7uVknuUnWI8J+qvhvqxXKwRlW+5WJddK04RNz+euEQvv5qGhs3bqTJkfWZNu0rmjZxrcoyqXAgk4PE/k7s1pILbnF7gX/6zUJOOq41C8eNYPfedK66y+2x06JRbR684UwURRCeeGMqC5euoWbVinzwxFWULpVGamoKX838NfsBWKrkiAlmas6DsVh/DyNtL/fnccedI7nksvCV12NC+C8BN6EORlDhPxE5HTfp/wjgGeBwYKSqflj87uUt/Jfr/Aqgs6puDGarU6fO+u30WcGKGflga/8Tm+IQ/qvWuLUOvCe4AOWbF3ZIGOG/oC1VVR3nX84DehavO4ZhJBIChdmFKiEKzgI6AAAgAElEQVTIN6iKyOMUsPZVVW8oFo8Mw0gokq37X1BLdUGJeVEwucdUP1XV4SJyLXALUBuYJyITVXVodFw0DCM/kiukFhxU3wQqqOqmwEwRqQbsLFavAihA+G8UMKqk/DAMo/CIkHBTpoJR0JSqJ4ET88g/BdtP1TCMEInUMtV4oaCg2ktV388j/z/k7LRtGIZRIMm2TLWg7n+et6qqGq4glmEYyYWQeCumglFQUN0oIp1UdXZgplcl3JzPNUYCY/NKjUIjJFz3PhgFBdWbgf+KyEtAVmDtDFyG21vVMAwjKAm3+XIQ8g2qqvqD3zH7H8DffPZCoEeW0JZhGEZBCDZP9SC8hMHtJeSLYRgJSJL1/pOuZV4gsaxfVBz2isNmOPby0ma69Z83075NS7oc3Y5zzz6TrVu3HnTNypUrqV65Ao8/9u9i98/shY+IW6YaLCUSFlQ9sa5fVBx6Q7HiY17aTH369mP2TwuY+eM8mjVrziMPPXDQ+Vtuup7+Jw8oEf/MXtFItq3/Qg6qIlKmOB2JNrGuX1QcekOx4mNe2kx9+/UnLc2NTnXt1p3Vq1Zlnxs39mMaNWpMq1aF28YuVu43Xu2FS7LNUw1FTqWriMwHlvjj9iLyVLF7VsLEun5RcegNxYOPAG+89gon+Vbprl27ePSRh7j9jhFR9y/Z7IVDhNVU44JQWqqjgFOBTQCqOhc4IdhFeWhLHRlw7kkRWS0iKQF5l4jIIRMhReQyEZkvIvNEZIGIDPL5r4nI8gD734VwL/kS6/pFxaE3FA8+PvTAfaSmpTH4/AsAuGfkCP7xf9dToUK+W+yWmH/JZi9cUiV4SiRC2fk/RVV/z/VhZORXOIA9qtohd6YPpGcCfwC9gGn5GRCR+rjZBx1VdZuIVABqBBS5WVU/CMGXoMS6flFx6A3Fuo9vvvE6EyeMZ9JnU7ODwcwZ0/noww+4/dZb2LZ1KykpKZQtU5arrwm+gXas32+s2wsHScCWaDBCaan+ISJdARWRVBG5Dvi1CHWegNtW8DlgSJCyNYEd+F2xVHWnqi4vQt35Euv6RcWhNxTLPn42+VMe/fdDfPDROMqXL5+dP3Xa1yxeuoLFS1cw7NrruHn4bSEF1Ej7l4z2wiXZxlRDaalejRsCaAj8CXzu84IRuA/qclU9078eArwDjAXuF5FSqpqepwWY6+tcLiJTgQ9V9ZOA84+IyL/864WqekFuA6ZRFfs+5qXN9MjDD7Bv3z5OPbkf4B5WPfXs82H7VhT/zF74CJCWaI/3gxBUoypsw3loS4lIaZzcdQtV3SEiHwIvq+oEEbkEpzU1LNc1AnQB+gCXA2+q6l0i8howvjDdf9OoMoz8KQ6NqnrN2+rfnv0oaLk7+zVLHo0qEXmRPGRVVPXKMOo7GagEzPdjZOWB3cCE/C5QF/VnADNEZArwKnBXGHUbhlHSJOA81GCE0v3/POB1WXIeMoXDEGCoqr4DICKH4br25fMqLCJ1gdqqOsdndQB+D7NuwzCigCSZoEooaqrvBh6LyH+AKYWtyAfOk4CrAmzvEpFvgNN81iUickbAZccC//bBdS+wgZzNXeDgMVWArqq6v7C+GYZRPLgx1Wh7UbKE0lLNTSPgiGCFco+nqupuoGoe5f4ScPhaHqbyknRBVS8J5oNhGNHHdqnKhYhsIWdMNQW3QfXw4nTKMIzEwK2oirYXJUuBDXP/5L09bsJ9DaCKqjZW1fdKwjnDMOKcCO5S5efJ/ygi4/1xIxGZLiJLRORdP7sIESnjj5f680cG2LjV5y8WkZMC8k/2eUtFpEiNxgKDqn/y/pGqZvhUPPOvDMNISLJaqhHaper/gJ8Djh8CHlfVZsAW3JRL/P9bVLUp8Lgvh4i0AgYDrXEzkZ71gToVeAYYALQChviyYRHKEPIMr0tlGIZRaCKxosovWT8FeMkfC+55S9Y89deBrIfcg/wx/nwfX34QMEZV9/mVmUuBrj4tVdXf/IPuMb5sWOQ7pioiaap6ADgOuEJElgG7cD8+qqoWaA3DCIKQEtqUquoiErgyZ7Sqjg44fgK4Bajoj6sBW32MAlgF1POv6+GnfarqARHZ5svXA34IsBl4zR+58ruF4nReFPSgagbQkZzobxiGUSjczv8hFd2Y34oqETkVWK+qs0Wkd1Z2HkU1yLn88vPyMOyhzoKCqgCo6rJwjRuGYURgl6pjgdNFZCBuAdLhuJZr5YAedX1gjS+/CmgArBKRNNwqzs0B+VkEXpNffqEp6DekhojckF8Kt8JYJtb1gRJVoyo/8tKuKiy/Ll5Mt04dslPNqofz1JNPMHLEHXQ5uh3dOnXg1AH9WbMmvL+hWH7/isNeYXFqqkUbU1XVW1W1vqoeiXvQ9IXfPOlL4Gxf7GLcJk0A4/wx/vwX/iH7OGCwnx3QCGiG65HPBJr52QSlfR3jwr5pVc0zAWuBO4EReaX8rovl1LFjJ92TrnmmnXsPaKPGjXXR4mW6bdc+bdu2nc6ZuzDf8sFSrNuLBx+nfPGVfjd9trZq3bpI9xnoX61atfSXpSv0z03bsvP//diTOvSKq6J+v9G217FjJ43039wRLdvqS9N/D5qAWaHYA3rjNlICaOyD4lLgfaCMzy/rj5f6840Drr8dWAYsBgYE5A/EbWm6DLi9KPdcUPd/rareHXa0jjMC9XyAbD2fo1qFN7Mi1u3Fg4/H9ezF7ytWhHVtXnz5xVQaNW7CEUccvCBw9+5dYa36ifX3rzi+M4VFiOzO/qo6Db+xvar+hntyn7vMXuCcfK6/D7gvj/yJwMRI+FhQ9z+p1kHEuj5QMmtURYr33x3Duefl7Is+4o7badqoAWPeeYs77ip8+yHW37+Y+DzELVMNlhKJgoJqnxLzogDkUK2r4T6/j4jM8XnfiEjTotSjMa4PFGl7xWGzOHyMFPv372fC+HH85eycBszIe+5j6fI/GDzkAp5/9hB5tKDE+vsXK5+HhJASiXyDqqpuLklHCmCPqnYISFmj7c8BF6jTwXob+Ff+JoIT6/pAyahRFUkmfzqJDkd3pFatWoecO3fw+Xz80X8LbTPW379Y+Dxc91+CpkQinjflUtzUCnBTJsKeAgGxrw+UbBpVkea9d985qOu/dMmS7NcTPhlH8xYtC20z1t+/WPk8TKMq9gjUugJ4QN0er0OBiSKyB9gOdM/rYtOoil8f89KuuuSyy4NfmIvdu3fzxedTePrZF7Lz/nX7cJb8upgUSaHhEUcw6pnC61/F+vsXCxpVkHhjpsEoNo2qSJGX1pXP/xB4SFWni8jNON2roQXZMo0qw8if4tCoatKqvd7/VvCH6oM71k8ejapYRERqAO1VdbrPehf4NIouGYaRDxFYURVXxOuY6hagkog098f9OHhLMMMwYoEknFIVDy3V3GOqn6rqcBG5AviviGTiguxl0XHPMIz8EOK35RYuMR9UVTU1n/yPgOCC4oZhRJVEa4kGI+aDqmEY8U2yaVRZUDUMo9hw3f/kiqoWVA3DKFaSrPdvQdUwjOJEkm5KlQVVwyiAaz9aEFF7o84Mf8PteMS6/4ZhGJEkAdf2B8OCqmEYxUqyBdVkm5dbILGuD5RsGlUlba9MKpQOSACd6h/OXf2b8vzZrTmiStnssl0bVuKOfk2y0/Nnt6Z+JXc+VYQLO9XlnpObcfdJzehYz22mdswRlQ+qI/eO+JHQ5CrM/ZYEtvVfEpORkcF1117D2E8m8eO8Rbw/5h1+XrQoYe3Fg4/RsLc/IycBrN62j+e+W8mSDbsPKjdj5TbumbKMe6Ys45Xpq9i0K51V2/YCMPCoGuzYe4A7Pl3CiMlL+HXDrhwfNMd+Rq69jC66+BLGjo/cFhbF8Z0JBwnhXyJhQdUTqOdTunTpbD2fRLUXDz7Ggr11O/bx5879BZbp0rAyM//Ymn18bKMqTPplA+A2/d2ZFaGDcFzPXlStWjWksqFQHN+ZcEi2/VQtqHpiXR8oGTWqStqekn/XvCC6NKjEjJXbAChXyv1JDWpTi3/1bcJV3RtQsUzOSutUcfZLlcBfXixoVFn3PwYpQKPqZRGZKyLzROQDETlkz9XCEOv6QMmoUVXS9gK7/qkpoWknNapajv0ZmazZvg9wAaRq+VIs27ibez9fxrJNuzmnfR0A5q3dwT5vP1OLP7DGhkZVKJ3/xAqq8fD0f4/XocrN9aq6HUBEHgOGAWGPxMe6PlAyalRF016mhrZmPbCVCq6rv+9AJj+u3g7A7FXbOa5RFQB2BQwDZCikFXNQjQWNqmScUhXzLdX8CAioApTD9d7CJtb1gZJRoyqa9lIk+BdKgE71Kx00ngowb812mtc4DICWNQ9jrW/FViqb04YJxX5RiRmNqhBSIhEPLdX8NKoQkVeBgcAi4Ma8LjaNqvj1sSTtCVAqYJPJjEzXWu1QtyJDjq5LhTKp/OO4I/lj6x6e/Pp3AJrVOIwte9LZuCv9oHr+O/9PLutan/NKpbJj3wFen+nGMU9sWi17qpYqpOd6fhUpTa5Q7rekyBpTTSbiVqMq4Hwq8BQwU1VfLciWaVQZhSWZlqkWh0bVUW2P1lc//jJouWOaVkkYjaq47f5noaoZOI2qs6Lti2EYh5JsD6riMqiKo2nWa+A04JfoemUYRl4k2zzVeBxT/RS4DXhdRA7HDdvMBa6OhnOGYRRMogXNYMR8UM1Powo4tkQdMQyj0Lin+8kVVWM+qBqGEcckYPc+GHE5pmoYRvwQiTFVEWkgIl+KyM8islBE/s/nVxWRKSKyxP9fxeeLiIwSkaV+1WXHAFsX+/JLROTigPxOIjLfXzNKwlx+ZkHVMIxiJGLLVA8AN6rqUUB34BoRaQUMB6aqajNgqj8GGAA08+lK4DlwQRgYAXQDugIjsgKxL3NlwHUnh3PHFlQNwyhWItFSVdW1qjrHv94B/AzUAwYBr/tirwNn+NeDgDfU8QNQWUTqACcBU1R1s6puAaYAJ/tzh6vq9+om778RYKtQ2JiqYRRApCfrt/7nxIjaW/jQwIjaizSFWIZaXUQCV+aMVtXRedoUORI4GpgO1FLVteACr4jU9MXqAX8EXLbK5xWUvyqP/EJjQdUwjGIlxKHJjaGsqPK70f0XuE5VtxdgO68TGkZ+obHuv2EYxUqkJv+LSClcQH1LVT/02X/6rjv+//U+fxXQIODy+sCaIPn188gvNBZUA4hnfaVYsWn2HMLBeldlcm18PbR3I5Y9OpAqh5UCoELZNEZf1onxNx7HpJt7clYX9/d9VN2KvP+PY5h0c08m3Hgcp3Sok23jgXPbZtvPb2/WaGtUQWR2qfJP4l8GflbVxwJOjQOynuBfDIwNyP+rnwXQHdjmhwkmA/1FpIp/QNUfmOzP7RCR7r6uvwbYKhyqmjSpY8dOuidd80w79x7QRo0b66LFy3Tbrn3atm07nTN3Yb7lg6VYtxcPPiaSvcxM93/jGybosXdP1a9+Wa+rNu3WTnd8po1vmKCPTPhFn5+6VBvfMEE73zFFt+zapy1unqgn3j9NT7z/S218wwTtftfn+ue2Pdr+tsna+IYJ2u7Wydn20zNU9x8o2v127NhJI/0316rd0Tp/1Y6gCZhVkB3gOFx3fB7wk08DgWq4p/5L/P9VfXkBngGWAfOBzgG2LgOW+nRpQH5nYIG/5mn8hlOFTdZS9cSCHlJJ2osHHxPFXu69U28//Sge+uQXNCBXFQ4r4x5xlC+Tyrbd6RzIVFZs3MWKjU50cP32fWzauZ9qFUoDsHPfgRK933AQIEUkaAqGqn6jqqKq7VS1g08TVXWTqvZR1Wb+/82+vKrqNaraRFXbquqsAFuvqGpTn14NyJ+lqm38NcNUw9vCz4KqJ971lWLBptnL216KuP1ZAfq0rsmf2/byy9odB5X5z7craFqrAt+POJGJN/Xk7o8XkftPul2DSpRKTeH3TTnKrmkpbmhBOFSdNRY0qiD5NqmO+aBagEbVWyKyWEQWiMgrfhA7bPL6UYonfaVYsGn28raXKjkB7+99mvL45CWHlOnZogaLVm/nmJFfcNqj33DXma2pUCZnck6NimV49Pz2/HPMvIOC7YFM2JfhWsK5xQqL4zsTFkkWVWM+qOI1qgJS1mj7W0BLoC1OTmVoUSpJJH2lRPExEeyliFMQABc7GlQtx4Qbj+Or23tTu1JZxl1/HNUrlubsLvWZPH8dAL9v2s2qzbtpXNNJslQok8ZLQzvz2KRf+WnlwdItWWRkHqqpFRMaVUSm+x9PxENQzRM/nqJ+3GMGB0+HKDSJpK+UKD4mgr3UgKCqQNe7pnL8fdM4/r5prNu2l9Mf/4aNO/azZuseejSrDkC1CqVpVLMCf2zeTalU4blLO/LRrNVMmrfuINtHVCufU0/KoZMqTaMqOsTD5P98Naoge+7aRcD/5XWxaVTFr4+JYC9FID0zuK2npyzl4cHtmHhTTwR4ePwvbNmVzqCOdenSuCqVy5fOnmZ1y5i5/LJ2B48MaXew5lWuemJBowpIvKgZhETQqHoR2KWq1wWzZRpVRrSJ5WWqxaFR1bZ9R/3ws2+Dlmteu3zCaFTFQ0s1X0RkBFADuCravhiGkQdy6FhvohO3QVVEhuJ2nOmjqiF0sAzDiAoWVGOOQzSqVHU48DzwO/C9nybyoareHQ0HDcPIj8RTSw1GzAdVzUejSlVj3nfDSHbciqpoe1GyWGAyDKN4saBqGIYROaz7bxiGEUESbMFUUCyoGoZRfNiUKsMwipNIa0r9mmu3q6KwN5SlX2GRXFHVgqphGMWGYN1/wzCMiJJkMTV+d6kqDuJVDymWbJq94rN31dDLaFi3Jp065MhmL140n4vO6MNZ/brzj0vPZeeO7QB8/78vGDywF2f1687ggb2Y/u1X2dd8Ou6/nN3/GM7s05XH77sjOz99/35E5EsR+VFE5olIRMYqkm3rv6jrRpVkMo2q+PLR7B1sb8oXX+l302drq9ats/NatztaX35vos5duV3veuQZveLam3Xuyu06ZuLXOmXGYp27crt+MOUHrVGrjs5duV2/mrtca9etr1/++JvOXbldTztriI5+e5zOXbldK1eppsDVfpOlVsCKov7NtevQUddu2x80EUSjKp6StVQ9iaKHlEg+mr2D7R3XsxdVq1Y96JoVvy2lU7djATim5wlMnTgOgKPatKdmbae82rT5Uezft5f9+/axauUKjmjUlKrV3N6t3Y7rzeeTfB2uwXi4N12JMCWac5Ns+6laUPUkih5SIvlo9oLba9riKKZNcdsJfjbhY9atPbT85xPH0rJ1e0qXKUPDIxqzfNmvrP7jdw4cOMCXn01g3Rp3TY2adQAuFJFVwETgH2HfjEck+br/MR9UC9Coek1ElgfkdyhKPZrHvrLxqIcUTZtmr+TtjXzkWca8PprBA3uxe+cOSpU6WKpt6eKfeeKBO7njgScAOLxyFW6/73FuueYSLj37JOrWb0hqmntevX3rFoDXVLU+Tv75PyJS9BiRZE3VeHj6v0dV8wuYN6vqB5GoJBH0kKJt0+yVvL1GTZvzwluu+77ityX874vJ2ef+XLua6688n3sfH02DIxtn5/fuN4De/QYA8MFbr5Ka4vYs2rplE8B7AKr6vYiUBaoD68O+KRIuZgYl5luqJUUi6CFF26bZK3l7mzZuACAzM5MXRz3CORdeDsD2bVsZdsk5/N8/7+LoLt3zvGb71i2895+XOHPIXwFIc63cPgAichRQFtgQ9g15RIKnRCIeWqoFaVTdJyJ3AlOB4aq6L/fFplEVvz6avYPt/fXCIXz91TQ2btxIkyPrc8edI1m+ZiNj3ngRgD4nn84Z514IwJjXR7NyxW+MHvUwo0c9DMBzb35Mteo1ePiuW/h10QIArrzunxzZuBkAtevUY9mSX64QketxOoKXaF5jEoVASLwx02DErUaViNQB1gGlgdHAMg2ySbVpVBmJRiSXqQ455XgWzpsT0Qh4dMfO+sU304OWq3pYWsJoVMVt919V16pjH/Aq0DXaPhmGcSjJ1v2P26DqW6qIezx6BrAguh4ZhpEXEsK/RCIex1SzNKreEpEauIeLPwF/i4p3hmHki9jWf7GH5q9RdWJJ+2IYRhhYUDUMw4gcida9D0bcjqkahhEfpEjwFAoicrKILBaRpVkrK2MRC6qGYRQvEVimKiKpwDPAANwOWkNEpFXxOFw0LKgahlGsROjpf1dgqar+pqr7gTHAoGJ1PEySakx1zpzZG8uVkt9DKFod2BjBqs2e2YsHe0dEsE4Afpwze3L50lI9hKJlRSRwZc5oVR0dcFwP+CPgeBXQLRI+RpqkCqqqWiOUciIyK5KrO8ye2Uske4VBVU+OkKm8mrMxuRzUuv+GYcQDq4AGAcf1idAm2pHGgqphGPHATKCZiDQSkdLAYGBclH3Kk6Tq/heC0cGLmD2zl7T2ShxVPSAiw4DJQCrwiqoujLJbeRLzu1QZhmHEE9b9NwzDiCAWVA3DMCKIBdU88Ks3Im0zYgug/UB9pGwl18LsGMc+j/jHgmouRKQncEGkApeItAZQVY3EH4yI9AGeF5FKRXbOUdPbLRWsYCwhIgNF5NIwrmsqInVEpEqE/KgrIikiUi4S9oBq3m5EPg8R6SQijYOXNCKFBdUARGQg8BawGjgQkB9WMBSRZsB8EXkDih5YReQk4BWgMVC1KL75a08BPhKR0cBIETkyXFve3gAROa8oNkKspz/wAAevsAnlupOAj4G7gSdF5PAi+nEy8AHwFHCLiBxWRHsDgDEi8gpwfVEDv7/fd4FD5IiMYkRVkz7hVmtUAN4B+vu8isBhQGV/nBKG3XrAWGA5MLaIPp4KzACOA0YBHxTRXhPgN+B4oCdwOzAdaBamvR5AOrALGFKMn9VJwEKghz+uDwwM4bo2wC/ACUBTXCBMA0qH6Uc/YJ63Nwh4CShVhPvqBywGegHnAk9m3WOY9noDS4AT/HE5/39acX02llyyliquBamqO3FyvNtEpAEwCffF/p+ItFbVzMK2ClV1NfAt0BY4XETeFJHWIpK/rGseiEga8FfgVlX9BtfSKicifQtjJxebgc9V9SvgG+B+4L/AGyISzhrwBsApOInj+0XkAu97JMeSqwG3AuNU9TsRqYXTqa8VwuWlgS9V9UtcL+QMXGAd43sUhfEjBRdMb/T25uF+7EaKyA3++1NYewOAB1X1f6r6Hi7gh7VhiH/PBwBzgB/8922UiIwC7i2sf0bhSPp5qiLSCSjj/0jvA/YB24B0VX3Wy/UOA7qpatBNKby9Uqr6gz9+AvheVd8VkZ+BFsDxqvq1iIgG+QBEpIv3Z5mqZvgAWwa4B1irqo+EYie3f7jW3ue4Fu8j/lwKrsWqOCnwjBDsdcYFqUVAJVXd4Md9XwRGqOp/fLmyqro3FB/zqaeKqm4RkfNxrbkFwCXAS6r6fJD7FWAl8DPwPnAO8BDwEXA2rpV4hqpuD/F+dwOL/edRCdfF/h8wH+iEaz0PC+V+A+ytBKoAq/0P+JlAH1Ud5suVUtX0QthbB9wBlMX9ALyAW9bZCjeOfp3mIeluRIBoN5WjmXC/5rOB9v64Ia4L9iPQIqDcy0DtQthrF5B3CnApUBvX3V5AiEMBuf3Lda4HbuehkLuIAfY6+OPmuGD494AyJwHPRsC/vv5+T8IFrgeA1DA/p1OAt7M+E+BCYApuVU1gucuBq/Pw72h/XBcX9J4NKFMNp8ZbMdz7zfV5Hwu8GM7nketcf2C8fz0Y98Ne4BBUHvdbHTdUFPj59grVP0vhpag7ELUbh5NxXfOsMdQaQDkfWJfjuti1/B/wPKBmIe1VBergWi0/A1sDzn0ONCikvSpAI/86q4dxI641EnRcMJ/7rQp0wY29Dff5lwBTcWPKEqZ/qf7/1kAGroXUJszP6VScsOOAXPlnAW8A5wUcz8kKcPncbwNcK/0z4BqffwFurLp6GJ9H0zzKXYQbOiowSBf0/vnjnrgf80HAXKB5Ie1V959v6cDvh/dvIu4ZQr6fr6XwU9QdiMpNuy9bJq7LB+6hzTe47hbAkf4L/QjwNdA6DHtfk/OQYCjQu4j+fY0bNggsdzHwIVA2THu9/XEzYJq/55+BthHyrzdunLpVmJ9TRdzT+qyHUuVxXdfWuDHHk3CtzNE+8BwV4ufRBTem/J6/33A/3+MDypQCrse1FCNhrxGww38vw7XXO1e5YbgfngLtWSpairoDUbtx16WcA7TDdSVv9Pkp/v+s1lblotjLo1xIswgKYa/AFnQI95t1n2V9KrDFVhj/cGOXLUKxl08d5XCt+tO9b//2x7NwPwJVcC3UseQK3MH8w7XkOgD1InG/uJkibxNiizwEe5VxAT9S9soCr4Vqz1L4KeoORPXmXZcpk5yub1aAGUBOqzXkLlIe9tIC8k+MoH8nA/0iaG9ghP0bQBGmA3kbWT9uf/FB9A9cq3Qwriv/DPA3X6ZCIf07BegV4ffvaAo57S6If60o5Bh0CP5Zd78EUtQdiHbCPfn9hZz5qJfgxtgamb0i2Wscpr2KuY4F153tmyv/EeDmBLjfmLZnKYzPINoOxELCtazmA1cTwhiq2Ssee7hxxOeBLgF5h7TWcJPjZxLi0EKs3m+82LNUyPc/2g7ESsI9Zd4fqS+g2QvLxlHAXbhJ+R0D8rOGAmoA1/mWWKHGBmPxfuPJnqXQU9JP/g9ERMqr6m6zFz17Xst9EG7600uqOifg3LG46Vmqqiui4V8y2zNCw4KqEVVEpDduetT/1Om5IyJtcS2tI3AT1WeLyBXAY7i9CdZFy1/DCIYFVSNq+CWeE3BjqU8AGar6mD/XCDcxvyJOk+hU3ET/uVFy1zBCwjZUMaKGqm4DxuPG/pYAA0XkNREZhFuC+zRuitA5wGALqEY8YC1Vo8QRkdrAn6qqfjPmJ4H/qupUEXkb6I7bRGYEsAX4WUPYzMYwYgFrqRolipAuAtIAAASVSURBVN8YexxQ3e+KJcB64Gi/Q3133AY0T+BWSy2ygGrEE9ZSNUoMv1P+7cB9qvqpiJRW1f1+f8/ZuE0+zlXV8b68Pb024g4LqkaJICJVceOkf1HVj0WkCXAnblXUehG5ErfD1LCsYBtVhw0jTKz7b5QIqroZOA24U0Ta4XaW+lFV1/sic4E+ItLcAqoRz6RF2wEjeVDVCSKSgdsf9TZVfUJEUlU1Q1Wn+4dUhhHXWPffKHFEpB9uKWo3Vd0mImXUpD2MBMGCqhEVvBzzE8AxfmjAMBIC6/4bUUFVJ4lIaeBzL1anar/wRgJgLVUjqohIBXXy4IaREFhQNQzDiCA2pcowDCOCWFA1DMOIIBZUDcMwIogFVcMwjAhiQTVJEJEMEflJRBaIyPsiUr4ItnqLSNamJ6eLyPACylYWkb+HUcddInJTPuf+6u9joYgsyirn92I9u7B1GUYksaCaPOxR1Q6q2ga3KfTfAk+Ko9DfB1Udp6oPFlCkMlDooJofftHAdUB/VW0NdMTtvWoYMYEF1eTka6CpiBwpIj+LyLPAHKCBiPQXke9FZI5v0VYAt22fiPwiIt8Af8kyJCKXiMjT/nUtEflIROb61AN4EGjiW8mP+HI3i8hMEZknIiMDbN0uIotF5HOgRT6+3wrcpKprAFR1r6q+mLuQiNzp61ggIqNFRHz+tb51O09Exvi8471/P4nIjyJSsYjvr5HEWFBNMkQkjRxdeHDB6w1VPRrYBfwL6KuqHYFZwA0iUhZ4EbfLVE+gdj7mRwFfqWp7XAtyITAcWOZbyTeLSH+gGdAV6AB0EpFeItIJGAwcjQvaXfKpow1u79VgPK2qXXzLvBxO4wrvz9Gq2o6c1vpNwDWq2sHf354Q7BtGnlhQTR7KichPuEC5EnjZ5/+uqj/4192BVsC3vuzFOEXTlsByVV3il5K+mU8dJwLPAfidp/Lqlvf36Udc67glLsj2BD5S1d2quh2nDlAUThCR6SIy3/vV2ufPA94SkQuBAz7vW+AxEbkWqKyqBw41ZxihYWv/k4c9viWWje8R7wrMAqao6pBc5ToAkVp6J8ADqvpCrjquC7GOhUAn4It8K3At62eBzqr6h4jcBZT1p08BegGnA3eISGtVfVBEJgADgR9EpK+q/lLI+zIMwFqqxsH8ABwrIk3ByZmISHPgF6CR360fYEg+108FrvbXporI4cAOnMx0FpOBywLGauuJSE3gf8CZIlLOj2melk8dDwAPe/FARKSMb2EGkhVAN/p6zvZlU4AGqvolcAvuIVoFEWmiqvNV9SFcS75lQW+SYRSEtVSNbFR1g4hcArwjImV89r9U9VcvdzJBRDYC3+DGNnPzf8BoEbkcyACuVtXvReRbEVkATPLjqkcB3/uW8k7gQlWdIyLv4jaw/h33MC0vHyeKSC3c7laCa92+kqvMVhF5ETduvAKY6U+lAm+KSCVci/lxX/YeETnB+7wImFS4d84wcrANVQzDMCKIdf8NwzAiiAVVwzCMCGJB1TAMI4L8fzt1LAAAAAAwyN96EjsLIqkCjKQKMJIqwEiqAKMAGNtfHuJ/EH0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from matplotlib import pyplot as plt\n",
    "from sklearn.metrics import confusion_matrix\n",
    "from sklearn.model_selection import cross_val_predict\n",
    "from sklearn.tree import DecisionTreeClassifier\n",
    "import itertools\n",
    "\n",
    "\n",
    "def generate_features_targets(data):\n",
    "    output_targets = np.empty(shape=(len(data)), dtype='<U20')\n",
    "    output_targets[:] = data['class1']\n",
    "\n",
    "    input_features = np.empty(shape=(len(data), 13))\n",
    "    input_features[:, 0] = data['u-g']\n",
    "    input_features[:, 1] = data['g-r']\n",
    "    input_features[:, 2] = data['r-i']\n",
    "    input_features[:, 3] = data['i-z']\n",
    "    input_features[:, 4] = data['ellipticity']\n",
    "    input_features[:, 5] = data['mCr4_u']\n",
    "    input_features[:, 6] = data['mCr4_g']\n",
    "    input_features[:, 7] = data['mCr4_r']\n",
    "    input_features[:, 8] = data['mCr4_i']\n",
    "    input_features[:, 9] = data['mCr4_z']\n",
    "    input_features[:, 10] = data['petroR50_u'] / data['petroR90_u']\n",
    "    input_features[:, 11] = data['petroR50_r'] / data['petroR90_r']\n",
    "    input_features[:, 12] = data['petroR50_z'] / data['petroR90_z']\n",
    "\n",
    "    return input_features, output_targets\n",
    "\n",
    "\n",
    "\n",
    "def plot_confusion_matrix(cm, classes,normalize=False,title='Confusion matrix',cmap=plt.cm.Blues):\n",
    "    \"\"\" This function prints and plots the confusion matrix.Normalization can be applied by setting `normalize=True` \"\"\"\n",
    "    plt.imshow(cm, interpolation='nearest', cmap=cmap)\n",
    "    plt.title(title)\n",
    "    plt.colorbar()\n",
    "    tick_marks = np.arange(len(classes))\n",
    "    plt.xticks(tick_marks, classes, rotation=45)\n",
    "    plt.yticks(tick_marks, classes)\n",
    "\n",
    "    if normalize:\n",
    "        cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]\n",
    "        print(\"Normalized confusion matrix\")\n",
    "    else:\n",
    "        print('Confusion matrix, without normalization')\n",
    "\n",
    "    print(cm)\n",
    "\n",
    "    thresh = cm.max() / 2.\n",
    "    for i, j in itertools.product(range(cm.shape[0]), range(cm.shape[1])):\n",
    "        plt.text(j, i, \"{}\".format(cm[i, j]),horizontalalignment=\"center\",color=\"white\" if cm[i, j] > thresh else \"black\")\n",
    "        plt.tight_layout()\n",
    "        plt.ylabel('True Class')\n",
    "        plt.xlabel('Predicted Class')\n",
    "\n",
    "\n",
    "def calculate_accuracy(predicted, actual):\n",
    "    correct = 0\n",
    "  \n",
    "  # iterate over the two lists simultaneously\n",
    "    for p,a in zip(predicted,actual):\n",
    "        if p == a:\n",
    "            correct += 1\n",
    "  \n",
    "    accuracy = correct / len(actual)\n",
    "    return accuracy\n",
    "\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    data = pd.read_csv('opticaldatafinals_SDSS.csv') \n",
    "\n",
    "  # split the data\n",
    "    features, targets = generate_features_targets(data)\n",
    "\n",
    "      # train the model to get predicted and actual classes\n",
    "    dtc = DecisionTreeClassifier()\n",
    "    predicted = cross_val_predict(dtc, features, targets, cv=10)\n",
    "\n",
    "    # calculate the model score using your function\n",
    "    model_score = calculate_accuracy(predicted, targets)\n",
    "    print(\"Our accuracy score:\", model_score)\n",
    "\n",
    "    # calculate the models confusion matrix using sklearns confusion_matrix function\n",
    "    class_labels = list(set(targets))\n",
    "    model_cm = confusion_matrix(y_true=targets, y_pred=predicted, labels=class_labels)\n",
    "\n",
    "    # Plot the confusion matrix using the provided functions.\n",
    "    plt.figure()\n",
    "    plot_confusion_matrix(model_cm, classes=class_labels, normalize=False)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}