LinearRegression_blank.ipynb

{
 "cells": [
  {
   "attachments": {
    "co2_data_mlo_s.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAARUAAADZCAIAAACB0TPOAAAACXBIWXMAAAsTAAALEwEAmpwYAAAAB3RJTUUH5AMTDwYqwrBY+AAAIABJREFUeNrsfWdAFMnzds3mJUnOLDmIICAIKKKIWTGcATNmUTGBZ0DPiFnEcCpwIp7IqaiYQMGAAQxIBhEk5yTsAsuyeXfeD+OPl//dycEd5nk+bZjp6empmq6urnoKQVEUcODA8a9AwIcABw5cf3DgwPUHBw5cf3DgwPUHBw4cuP7gwIHrDw4cuP78G0il0k/aPoqi3b+EVCr9xx0zFEX5fD6fz//T7wKBgM/n/8fbkUgkf2pBIpF8bWMuFov5fL5IJPpsD/FLgbhr166vuX8lJSW///67jY0NhUIBgIaGBiKRSCaTe/ESaWlpeXl5xsbG3dG0uLg4Nputra39sWOam5vDw8MjIiLi4uJ4PJ6FhQWBQBCJRFFRUefPn4+NjS0pKWEwGAoKCt3pW0tLi1AopFKpHb8cP368pKTExsYG+5qYmHjv3j0zMzMajdZbA8JkMqOjo01MTLAx7+nLKCUl5ezZs1evXk1NTaXT6bq6ugiCXLhwgcFg0On0XpeQxsZGqVTaeYg+K9CvGCKRaPv27QAQHx+P/bJly5asrKzevcrBgweXL1/ezYNjY2NTU1O7OMDd3X3UqFE3bty4cOGCvr7+nj17UBT18fExNja+ePHizZs3PTw8bG1t29raunO58PDwGzdudP6FQCBoampWVVVhX+3s7FRVVQsLC3txQHJzc21tbevq6v7FuWlpaSYmJnv27Ll58+aRI0fMzMyePn2Komi/fv1KSko+hZDs27fv7t27X0pEv2r94XK5bm5uM2fOnDlzJoqiGRkZDAZjwoQJT58+vXjx4uHDhwcNGrRnz57ExMQhQ4Z4e3tLJBIURR89ejRixAgXF5fIyEgURXk83tKlS+3s7Pz9/Xk8Xl5e3vbt29euXevg4BAcHIyi6NGjR729vX18fAYNGnTlyhWxWIwpla2t7aJFi9ra2hobG69cubJjx46AgIDbt28nJSWhKBoaGjpw4MBJkyaVlpZ2FndlZWWRSIR9TUpK8vT0TEpKUlZW7pB4FEUNDAyOHDny1/uVSqWPHj0aMmSIs7PzpUuXGhsbBwwYMHjw4Ojo6I5j5OTkBgwYcOnSJRRFq6qqCASCu7t7UVERj8dbvny5ra2tl5cXl8utq6sLCQlhs9koit6+fTsxMfH58+eXL1/29fUdNGjQjRs3pFIpm81euXKlnZ3dnDlzmpqaOi6Rl5fn7OxcX1/f8Uttba2np6eDg8OWLVuwNl++fDls2DBHR8fQ0NDOtzB9+nRfX9+O29m7d+/EiRMxPd+3b5+zs/OcOXMEAgGKohERES4uLmPGjHn58iWKoiwWa+bMmQ4ODuHh4RKJJCsr68iRIytWrFi7du39+/exBjdu3FhRUXHr1q3BgwcPGTLk6tWrhYWFVlZWbm5uMTExbDZ79erVDg4O3t7eLS0tKIpu2rQpJiZmzpw5hw4d4vP5KIpmZWUFBwcLhcIfQn+ePXtmb2+fn5+vqKjIYrFYLJabm1tgYGBZWVm/fv1mzZr16NEjWVlZKyurR48emZqaHjlyJC8vz8zMLDIyMi4uzszMLD4+ftGiRVOmTMnPzx81ahT2OwAcPHjw1q1bSkpKFy9ePHfuHI1GCwoKunDhgrGxcVVVVWhoaN++fV+/fr1w4UIHB4f6+npDQ0N3d/eYmJh169adP3/+4sWLKioqDx488Pf3d3R07Ojw3Llz58yZ86e7CAsLGzlyZOdfNm7c6OHh8df75fP5urq6t2/fjo6OBoDa2tolS5asWbOmuLi44xh5efmAgIAVK1agKLpr166JEyfOnj27urp69erVU6dOffv2raOj45YtW0pLS8eMGdPQ0ICi6IYNG06cOPHbb78BQFBQ0Llz5+Tl5TkcTmBg4MyZM/Py8kaOHDl79uyP6Q+Hw7G2tl67dm1SUtLIkSO3bt0qFov19fWjoqJiY2MBoPOrwcjI6NatWx1fHzx4YGNjI5VKtbS0JkyYkJiYOHny5MGDB2dlZVlaWiYlJV26dMnNzY3JZJqYmPj6+j5//tzAwCA2NvbWrVsAsGTJklOnTpmbm2OvRQ0NjcbGRhsbm/j4+KtXr8rLy+fl5f3000+bNm3Ky8ubO3fu7Nmzk5OTlyxZMmnSJKFQKCcn169fvwsXLlhZWaWlpYlEop9//nnPnj3YK7JX8FX7Dw4fPqysrFxdXU2hUE6cOKGkpKSlpeXg4KCpqclgMPz8/EaMGDFy5EgfH58RI0ZMnjyZyWSmpKSMHTt27ty5Y8eO3bRp04MHDxgMBovFio+PX7Zs2fDhwxEEcXZ23rx58+TJkw8fPvzHH3+IxeKZM2f6+vpOmjTJ2NiYSCSGhIQMGzaspKTEyckpLS2toqLCzMzs6tWrHh4eEomETqefPXv28OHDo0aN2r17t4+PT0eHEQTpbOKjKPr+/fv29vY/WedkMpnL5f71fslk8o4dO/h8fnNzM4VCaW9vNzMz69ev35/WZhYWFtXV1Uwm89y5c4sWLeJwOEKh0NPTc9q0aZmZmXJyci0tLUQikUgkfljjEokIgkil0rFjx/r6+s6dO1dJSQl7GY0bNy4tLU1bW7u2tvZjTyE1NVVBQeHEiRNDhgzZvXv3ixcvhELhL7/8AgBNTU2ysrKVlZWdj+98s2QyGZMzNTW1gIAAV1fX6OjogoKC8vJybW3tO3fuEAiEffv2VVRU1NTUODs7V1VVWVpahoeHk0ikcePGhYWF+fj48Hi8wsLCK1eurFmzRlFR8eeff25qamIymbKyslQqVUdHx8rKSl1dvaqqau/evU5OTseOHWttbS0oKBAKhc+ePfPy8nJ3dw8NDZVIJBkZGSNGjOgYme/Z/1ZXV3fv3j07O7unT5+OGjUqNDS0swOKRqNhK2YajaaoqAgABAIBALhcrqysbIepw2az582bt2XLFiaTefTo0aioKARBOh6wsrIyNpXLy8sDAGZ3EQgENptNo9Ha29tpNNqFCxcQBFFRUREKhR1OJD6fr6SkhMmHubl5R5+trKySkpI6vpaXly9cuFBHRyc7O7uzl+z58+d9+/b96y2/ffs2NDS0uLhYQ0ODQCBgQo+h82Ha2tqKiooXLlygUChDhgwRi8VUKjUkJARbJSooKBAIBOxGEATBOomNm66uLvYBQRAEQf7444+4uDiBQKCjo4Md+X8kg0DocKbJyclhfcCGrqqq6syZM+/evVNVVSUSiR1HAoCmpmZBQUHH18rKSgUFBQRBiEQi9ryIRKKMjIyamtr+/fv19PRiY2MPHDhQWFhIo9GkUmlLS8uMGTNWrlwpEolkZGSwRubMmXPixIl37955eXllZWWFhIRUVFRoaWmRSCTMsyeVSsViMXYV7DlSqVTM+4c9poULF16/fr2oqIjP51tbW/8Q/uszZ864uLgcOnRox44dFy9e5HK58fHxZDJZIBBgA9fZ/4F9IBKJ1tbWT58+bWtrQxDk8uXLgwYNWrZsWUVFRUBAwJQpUwoLCwGgoKAAe2UeO3bM2dkZO7ez73Xo0KGNjY1Lly7V1NTctm2bhoZG5wOkUqmzszOmzwkJCYsXL+74a9WqVZWVlZjPg8/nz5o1S01Nbdq0aYqKinPmzMGOOXnyZGJi4qZNm1paWg4ePHjz5s2O0zMyMmg02tatW+Xl5fl8PoVCwZwonQUUABQVFR0cHDZv3rxkyRKsY2Kx+Pr16+vXr583b15zczPmoqytrRUIBEKh8MqVK5gnrbOSIAjy4sWLhQsXLl26NCsr608uTR6PV1xcXFRUVFxcrKiomJ+fn5ubCwAPHz7U0tLKycnhcrk7duzQ0dFpa2vr7KZbvXr1gQMHsINLS0t37tw5b948BEGYTGZMTAwAXLp0SSgU1tbWBgQErFmzJjAwkMPhaGtrczgcDQ2NFStWJCQkvHnzBnsFYG3OnTs3NDRUU1NTXV09KytLW1t769atysrKzc3NmCTweDxFRUUikZiamoq9nlpaWvT09Dr85nZ2dtra2nPmzBkzZgz2rvz+/W+2tra3b9/u+Lp27dqVK1cGBQXJy8tHRER4e3vn5eWhKLp8+XLMQ7Vnz54DBw6gKLpmzRp9fX0DAwMPD4/29vakpCRVVVVjY2N9ff3MzMz4+HgajWZnZ6ehoTFo0CCJRBISErJlyxYURZuamjw9PYuLi5ubm62trXV0dFRUVOLj4ysrK729vevr6yUSydatWyMiIkQikb29vZaWlpaWVmdzH3NAaWtrq6mpycvLT5w4EdsvamxsdHR0VFRU1NbW1tTUTEhIQFHU39//ypUrU6dO7Vgfs1gsBoPBYDDGjRvn6Og4d+7cmJgYFRWVgICAzguMd+/epaSkGBgY1NfXV1VVzZkzp76+fvPmzUpKSubm5jNmzNDW1s7KyvL29tbX19fW1p41a9aFCxfCw8P9/Pwwr4y9vX1tbe2BAwdUVVX19PTWr19vZGSE+UVQFC0uLlZSUlJQUJCXl5eTkzt8+PDNmzdVVVXNzMysrKwKCgra2tqMjIwYDMbw4cOHDh06evToju4JhcJdu3YpKipaWloqKipu3ryZx+OhKOrs7Ozq6mpgYKCkpJSQkMDlcj08PAwNDQ0MDHx8fFAU/eOPP+Tl5XV1dV1dXdvb22NiYhYtWoS12dLS4u7uHhERIZFIKioq9PX1DQ0Nhw0bNnToUH9//5CQEFVV1f379798+dLAwMDS0lJbWxubinV0dDp8OcHBwXQ6HZOZ79//hnnS/uqeQlEU8/90DUywOp9YVFSEDWVcXJyLiwuKohUVFV03Ulxc3PVCs7y8vL29/W//Kiws/Kv/t7KysqioqEPOKioqsrOzvby8nj9/3nGMWCwuKyvDPmM+7m56ujEvWWtrKyZw2Fi1trb+7Uh2oK6ujsPhYE6CDlH72F5CZzcGdvvYB+yincFisQoKCjr79DCUlZV1vkpFRUVnLx+Xy+2Oj5vP53dcGus8m83G7lckEhUWFv71qbW1tZ0/f97e3r7XBRX50fJPnz59evLkyRs3bnzxnty9ezc/P9/V1dXJyQkPhPmk2L9//5UrV06fPu3q6tq7Lf9w+iMWiwUCgays7JftRmtr65o1a8rLy1VVVTdt2oQtw3B8InC5XJFIhHkycP3BgeNrAekLXvuvEZZ/4x8kEKCH0YdEIrE7UZ7/ZxRIJLFY3LOB+/SnIAhCIBB6FB6Kua1+2Hv/76KCIAiFQun+NPXF9Afb5PrHCEUsgFJdXb37415eXq6hodH9eEoikZidnW1tbd19sSMSiRkZGXZ2dt1/WiiKvnv3rl+/ft0/hcPhtLa26urqdrNjCILU1NQoKCjIycl1X+awiI0/uch7994BIDc318bGpvv6IBAI6urqDA0Nu/9QGhsbSSQStuHTzXsvLCw0MDDovB0iEonk5eVXr17dzdjZL2a/LV++/PDhw9gWRxcC8eLFi7q6Ok9Pz24OPZVK/eWXXxYsWMBgMLrZE1lZ2cmTJ0dFRXX/6crKyg4cODA1NbW9vb37+uPt7f3HH390/5Tc3NwnT574+fl1TgTo+h0fGhrq6OjY/S1CWVnZ6dOnnzt3rrMM/eMp7u7ujx8/7v6NkEik8ePHJyQkdP+U2traCxcu7Nq1SyAQdFOro6KilJSUsFCp7pxCp9PXrVu3bdu2zjtCbW1tFy9e9Pb27maA/Bebf7BQl38MO6dSqWQyuUfB+WvWrFFVVe1RQDuCID2N/79w4QImTN08XiqVYu/47p9ibW2tr69PoVC6n0eAjVWPvCMHDx78dwvrHl3lzJkzPTqFTqcTiUQSidR9xfbw8CCRSB1RC92cguh0eudeYcEZ38b6pzvo06dPT01tHR2dnl6lf//+PT3F0tKyp++Lnp4iIyPTI2nA7r37xhsGExOTnt77gAEDenpK5yinbtoRhoaGPTpFRUWlp70yMzP7q372SH++mP3m7e198uTJf5wlsGi33k2Y+yt4PN6nSO36/FcRiUR/Ckj7FOByuT1V7J4C2wz9Fwl8PfVgUanUzgrDZrPPnj27bNmyr91+6/5qtRejZbuwFj7DvXyGq3zqF03HxPgZzPtPrTwA8N+TdnH+EBw/JI4dg3fvQCIBofC/NEPCRxLHj4iICFBTg8JCUFSEoUPx+QcHjp6gogLa2s48f17030w4XH9w/HhgMsHCIu/hw6inT1W6wbuE6w8OHAAoCiIRAEhDQ/MGD579+PFaDkeJRPpW9YdEwldfOD4jWlvh2DEAuH7z5rqsLL8RI6bNmoX8X+9uT2Xyi+mPWCx+/PhxY2Mj/lhxfCZUVkJw8Lt79/YymQbq6gvk5IDHg04Kw+FwHj582NLS0v0t1C+mP1gC49/S0ODA0Zt48QKyswEAUlLA3PzE0aNOFhZnz50DFgukUui01ywQCEpLSwUCwTegP2QyefHixfr6+vjzxfEJwedDWBg8eAAAkJ2dYmxMys4+aGgINBqoq4OaGnSy31RUVLy9vTU0NLofWv4l1z89DWzDgaMHqKsDiQQkEmhrQ+vqQCwuT0o60dg41d1dxdUVEATGjgUHB/i/65+eyuQ37H/7pjNnv8LOf1eZyPX1sHMn8HggEoGCAiIvH3Ptmmt2Nikry7WsDPr1AwCYMQPc3P7jdb5e/UlPT+/bt6+mpqafnx+Px/vrAQiCeHp6NjU19ahZkUi0a9eu6urqjx3AYrFev35dUVERFBRUX1/fmQ8RQ01NTUVFhUgk6qB06xFmzJixcOHCDgvb19dXU1NTR0engyDy8+PevXu//vprr3MDfElUVcHNm8DjgUAgBpj1+nXgunVPAwN/3rCBVFwM3c6x+0d8vR5kDw+PR48e9evX76effsrJyXFwcKitrRWJROrq6iQSiclkysnJzZo1S1ZWlsVi8Xg8oVCoqqoqIyPT0tLC4XAAQENDg0Qi1dfXi8ViRUVFjKYU47LicrksFgujldDQ0MDYBrFTEhMTg4KCbt++7e7u/ubNGz8/v7S0ND6fj2VZsVisuLi4+vr69evXz549GwAaGhra29vpdLqWlhbG2MThcPr06aOkpIRJJI/Hw2psqKio1NTUVFdXY7TRAODj41NZWVlbW9vY2Ojk5GRoaDh69Oi6ujqMEFBZWZnJZLa2tmKNt7W1MZlMANDU1BSLxTwej8Ph0Gg0TU1NPp/f0NAAAKqqqhh/L4vFIpPJcnJyMjIyTCaTRCLJyso2NjaiKKqqqkqhUFpaWvh8PoIgGL9xXl5eQ0ODVCpVV1d///69lpaWUChsaWnpfubv1wUuFxAEGhurystDrl51HzRouZcXuLsDgQCKitB7FIpfr/707dv3xIkTHh4e/v7+/fv3v3379pkzZ5SVleXl5WfNmrV06dKAgIDAwEAHB4cNGzZ0sN3u2bNnw4YNKioqPB6vb9++lpaWt27d0tLSAoDDhw9j4ckEAqGhoWHx4sX9+vWrra0dM2YMiUSKjY1lMBjm5ub5+fl1dXW3b9+ur6+vqKhobm6+ePFiWVnZhg0bSCTS1KlThw4dmpSUNHbs2KCgIF1d3V27dqmpqTGZzG3btu3cuVNRURErR3P06FEVFRWRSHTy5MnMzEwKhUIikRwdHauqqm7durVkyRKJRBIbG5uTk0MgEDQ0NHbv3n3q1CktLa1NmzZpampWV1eHhITs3buXRqNVVVXt2LHj8ePHRUVFbW1t9vb2ysrKd+/eVVFRqaysPHr06B9//FFZWYkgiKGhoa+vr5eXl7q6ekNDg4mJycyZM2fMmLFp0yYWi1VdXY2iqK6u7rBhw1atWjVy5MjCwkKMvjAhIQEAXr9+HR0dPXXq1JSUlCtXriQlJZ09e/ab1B8OB8zM2hsb5x49arpgwVaRCDIyYOdO4HBgwQLovSj1r9d+u3HjxtChQx8+fLhy5crc3NyTJ08eOnTo6tWrixYtwhh0vby8Oniud+zYceXKlaqqqvT0dCMjo9OnT58/fz4rK6ukpKSqqkpeXt7FxeVPeRBCofDUqVMBAQFsNpvNZldVVampqdna2i5YsMDS0nLUqFFMJnPZsmX6+vpTpkzhcDjYZEKj0VxdXRctWmRnZ8fn8y9fvjx//vywsDBvb++7d++yWKzly5eHh4eTyWQWiwUAtbW1hYWFx48fj4iIIBKJKioqNjY2S5YsAYD29nYikdinTx+sP1paWlwu98GDB3PmzDl//vyOHTsyMzMzMjIMDAy4XO7Zs2cdHR1NTExEIlFycrJEIhk5cuS5c+cmTJhw//793Nzc33///bfffsNI5SkUyu+//75x40Y6nY6iqLm5ua+v79ChQ/v27YsgyLt379ra2vT19YODg/ft2xcTEyMWi4cMGRIcHMxgMJqamoYOHRobG3vx4sV169Z9o+absL29iEbbuWnTqIKCE7/8IisvD0VFIC8PGhqweTP0XiLJV6o/lZWVM2fOnDdv3q+//urp6fnq1SuxWIylNLW2thKJxM40EVQqFcvDU1BQ4HK52GEUCoVMJg8cOPDSpUsMBuPMmTMVFRWdL6GgoIBVByCRSC4uLnfv3pWXlw8ODm5sbJRIJFgeL0aWKZFIMF4bPp+PkcKIxWKpVIogiEQi6dOnj1QqlZeXx37U09MTCoUduSsY1T12uoyMDJ/PxwhBsQ7Q6fSXL19iRz548MDCwkIsFmPZNUQiUSgUOjk5TZs2bf/+/WPHjj137pylpeXPP/+MXUhdXR0bE6FQKCMjgzFJyMrK8ng8TCfl5eXJZDKCIKqqqgAQEBDAYDCWLl2qq6srlUoxaxY7AEEQbDzpdLpQKJw1a9bJkydVVVV7mjT6VXirW1tBIskqKhqakHCPRNpuZyenpQXy8qCqCgBAIACdDr230vtK7TcGg4Gi6OrVq93c3O7du3fgwAGJRLJjxw47O7snT55s3rwZUxgZGRkEQUpKSkJDQ/X19evr64cNGxYcHIx9lpGRSUlJiYuLc3d3V1dX78h1xQQay2YjEAiY2qiqqjo4OGhpadFoNCaTmZOTIycnRyaTW1tb37x509jYGBsbGx8fr6qqigm9q6srjUYbPXr06dOnBQLBlStXZs+enZqaiokjlUrFkkC1tbUJBMLvv/+ura2dnZ29Zs2a69evd9xmYGCgl5fX9u3b379/HxMT8+rVq9evX589e1YgEGzbtu327dtXrlzJzMwMDg7GVnpSqfTXX38Vi8VsNhszSjHK/LS0tEOHDikoKDCZzJ9//nnYsGGXLl16/PgxgUAgEAhYlhh2p3FxcRkZGQMHDkxISLh27drdu3eHDBlCpVIxhafRaEQi0cDAoLy8fMKECZ8hc7GXEREBiYlw8WJ+cfH2tWvny8sD9nrS14dPw1D5xeqfxsbGjhs3rotwozlz5rx//76iosLb23vQoEHOzs5YmYp9+/YpKiqamJhoamoaGRkZGRmlpqZaWVlJJJKgoKDW1lYejycnJ6eoqLhjx45hw4Y1NzfX1dUtXLgQq6JDIBB0dXVNTEwsLCwYDAaNRtPQ0Fi4cGFjYyOHw1mwYIGdnZ2cnBxWaMjCwkJTU1NLS2vEiBF5eXkzZsxwdHQcMGAAj8czNDS0s7MbMWKErq7u27dvR44cOW7cOAaDYWZmhhWlYTAY2JpnyJAhlZWVbW1tO3bs0NbWNjAw6OAGMjU1HT58+Js3b8hk8tGjR9XU1IyMjJSVlSsqKgIDA01NTS0sLLKzs2fOnDl9+nQNDY3c3NzFixcbGBhYWFhYW1srKCgoKSkxGAwvL6+ioiIAWLNmjbKyMp/PZ7PZWKmiSZMmmZmZ6erqOjs7p6Wl2drauru7czgcHo9nampqbm4+e/ZsZWVlU1NTDQ0NHR0dExMTzMUyc+bMf8Ek8YVx5Qq8eMFetuzM5cszZ85k/P47WFnBhAlgZgbjxkE3YtsEAkFGRoa9vX13+We+FEP88uXLsZJ6/xECgeDQoUOYZwkjbg8LC+uaCv37hkgkWrt27eLFixcsWPDixYu/PSY1NdXf3/9jLWzfvt3f3x+jZv82wOGg2BNftQo1NLwWEzNuyJCWtjZUSQndv79HLbW2tgYGBv6VEf8L1F9ob2/vIP+XSCR/6lNv6U9HXYa/fv5hIZFIOBwOl8v92GhIJJIuaoDy+XysROk3A3d3NCQERVHUzw+dPn2Yri5r0yYURVFjY7RT6dhPoT+97D/YsmXLli1bACA0NNTBwUFXV/fBgwclJSXDhw93dHScMmVK53lP/D/8bVObNm169uwZAPD5fMzIbGho2Ldv36JFi1atWpWTk/O/jbIqzDmLfa2urt6xY0fnTdW6urply5b5+fn5+fmtXLlS+JF897q6OolEcu3atb9uyBYVFfn5+c2fP//XX3/tIjKKzWZ3PxwWRdFLly6tXr3az89v9+7dWGGvXsH79+/3798vKytLp9M7tkTfvHkTHh7eEWFAIBA6mEbCw8Pr6+s7t9CxHPoYuth9/jIoLIT0dOBy2wmECS0trlZWSmpqAABnzoCHRzfbwNxCWO2TL+N/e/ToUXJyslAoFAgEsbGxCQkJcXFxq1atevbs2ezZszG3aXx8PHYwVjLt6tWrWK3cvyIqKmrx4sVNTU0SieTKlSsAsH37dhaLtXnzZjc3N29v7/z8fADA9joxTUNRNDk5+cCBA69eveocT3D16tXJkyePGDFi9OjR2ObPnTt30tPTAaC+vv7NmzdpaWlDhgxJTk7OyMjIyMi4e/duRz1QsVi8d+/evn377t27F/Pq1tbWYkyleXl5TCYzISHh8ePHWEnQ3bt383i82tra27dvY0JWVlaWlZX17NkzPp//4MGD4uLijo7Fx8ebmJhMmDABQZDdu3cLhcKysrKYmJiysjLM8X3nzp38/HzscT548ODFixfYddPT0+/cufP+/XsAyMnJyc3NFQgE7969i4mJqaysbG1tvX79em5u7vPnzzuISysqKh49etSh/A0NDXfv3n3x4gUQuiTXAAAgAElEQVQAPH78uLS09O7duyUlJdi/iYmJ8fHxbW1tAJCVlfX27Vs+n5+RkREbG1tZWdnU1GRtbf3kyRMASE5Ojo2NZbFYYrG4srIyJSUlLS2ttbUVAIqLiz95ckpNDQBAYyM4OsL799X5+V43bowRiQIGDABMf0aPhu6R+PD5/EePHl29evXGjRuNjY3dJwDrHf3BJr6oqKgVK1YoKCgQiUSRSKSsrGxjY1NdXV1XV4dtY48aNQpb5mKvQCsrq/79+1tZWX2sza1bt/r4+BAIBCqVyufzMzMzt23bZmFh4enpOXLkSOzxR0VFHTlyBCtqi+UU7d69OyIionOYD51Op9FoVCpVT08PRdGlS5fW1NRs3bo1IyPj3Llza9asyczMxOrPMZnMyMjIkpKStWvXstnsjkaio6OLi4vDw8NtbW33799fX1/f0NCwdu3akJCQ7OzshIQEX19fAoGAFQlftGhRTU3N4sWLa2pqfvnll4iIiN9++238+PEZGRl79uzpeHkTiURjY2N7e3tzc3M6nV5YWLh8+fKamppdu3aVlZVNnz6dxWL5+/uXl5evXr368ePHUVFRUVFR8fHxe/bsKSoq8vX1xbwsZ86cefr0qbe3d1lZGVZ5u6Sk5Pr168HBwXfv3sV0hkgkdkwpzc3Nq1atqq2tvX79+q1bt1RUVPbu3VtSUrJp06bKysqAgIDw8PDk5OQNGzYAwOzZs8+cORMREbF///7a2trAwMCcnBwikdjQ0BAWFnbx4sWSkpKFCxfm5eUtWLDg9OnTCQkJFy9exMKUukm9+++xciWw2VBeLunbNzEjY46Xl8mYMasXLoS4ONDW7pkbmkQyMTGxtra2srLC6ud+Vv81giABAQFZWVnNzc3FxcU2NjYUCoXP52O1XalUKvYUW1tbO7iGMcdrF3aCSCTy9PS8f//+jRs3sF0XbAsFm2rV1dURBImMjLSysqLRaPn5+U1NTQiCPHv27Pz58ydPnqyqqsIqYAIAl8u9c+cOl8vV19d3cHDYtm1bTU1NXV1dcXExgUBYuXLlzJkzQ0JC3NzcYmNjN27caGRk9OLFCxaLpaCgQCKRTp8+HRUVtXv37qysrDt37iAIUlpampKSMnDgwIaGBiaT6eXlNWTIkFOnThkYGLx8+VIsFuvq6lKp1Fu3bsnIyCxatKi0tDQqKmrz5s2rVq1qbm7G6vjS6fRt27apqqoWFBRs3779xo0b+vr6ampqLBYLqzdaXl6+c+dOVVXVyMjI8+fP5+fnX7lyJTw8fO3atenp6SkpKdg2V2Bg4IYNG3x8fDw9PTFmUBMTk127dkVERDQ3N/9VFGRkZHbv3l1QUMBkMrOzswUCwfr160ePHs3hcN68eXPkyBEmk0kmk+3t7Zubm7H62CKRyNzcPCcnp7i4uF+/fnQ6feDAgUFBQcuXL7e2tq6vr79165aamtqpU6cqKioOHTo0cOBAOTk57DY/FaRSePwYKiuhqenwkyeHW1vjPDyc588HqRQKC+F/j777+oO5Z9va2jCny+e237y9vY8cOeLk5IS5VlVVVdPS0qKjo0eMGGFsbJyZmdnY2JiQkNC5ANg/9lJGRsbPz+/q1avNzc00Go1Op5eWlmKhA7m5uWpqaiEhIbq6unfu3CkoKCguLg4NDVVSUoqMjDQzMwsKCupoR0VFZd++fceOHVu/fn1DQ8O0adOqq6tVVVUJBEJHXWhMpVVUVGRkZLC1AdY9kUi0efNmT0/PxMTEXbt2nTp1Ctvy//3331esWLF58+aRI0ceP3585syZdDqdTCbz+fympqa3b9/q6ur27dsXRVHMi41VxUCQ/0/42t7eHhoa+vjx45iYmOzs7JaWltra2oKCAl1dXWVl5aSkJCqVOnXq1MTERLFY/PbtWzabPWrUqDNnzmzfvl1BQQF7lWCMzwKBAPtqYGCAbZtiy5gOOwTzrIhEIrFYXFhYuHTp0ry8PCMjI8zBgJ1LoVDEYnHH0gir7UEkEhUVFYODgw8ePCgWi7HtLBRFsa1nMplMIBAUFBR4PJ68vDyVStXW1lZVVd22bdvq1as/7eSTlQVKSlBa+jwj4zqCvDI3d7a0BFlZsLQEGRn4t573ngah99r+qampqampKbaBYGJisnjx4gsXLnC53JCQkD59+qSlpW3cuHHGjBndZ5rGnpOTk5OZmVlycjKRSJw2bdrOnTtdXFyw7RSM6Dk8PBwL9jl37lxaWtrjx481NDSYTKa1tfWhQ4ew+Y3NZu/bt4/P54tEolmzZmF1f4VC4ZMnT0xNTTEhI5PJUVFRHavqDsnD5Gn9+vXW1tYvX76cPn36yJEj/f39VVVV9fX1161bp6WlNWHChPT0dFlZ2ejo6G3btqmrq2tqar569QpjZMYW8ViDHTXlMdHncrlYAASXyx0/fnx+fr6WllZhYSGCIIsWLfLz83NxcdHS0ho5cmTHyp7JZDo4OEgkEizok0wmi8ViT0/PY8eOFRcXJyYmLl68GNv37BxPTSKRysrK9u3bh4VEMBgMAwODhw8fYuNz/Pjx7OzspKSkCRMmzJ8/f/Xq1VpaWgYGBioqKli3GxsbHRwcqFRqbm5uZWWlnJzcnTt3LCwsQkJCTExMXr58uX79+j/++AMAsP2o8PDwW7dufVr9uXtXtH79qdOnr+fl7Q0Ls7h1C3JzwcMDlJTg2DHo3SLbXVheny7rg8vlYq9ezOLCopI7z1dd81+/fv3a0dERQRDMJhw4cKBQKExJSamoqFBSUnJycsIK5WL05+3t7Xl5eWKxeNCgQdjpr169srW1pdPpXC43PT0dW0lLJJLBgwenpaVxOBw9Pb3m5mYzMzMKhaKiopKTk9Pe3q6rq6umpkahULDKMNjU1NLS8vz589bWViMjI6z9oUOH7tq1y93dvba2Nj09nUAgODk5EQiElJQUNze3oqKid+/eGRoaOjg4FBYW6unp8Xi8pqYmMzOzkpISTU1NjPC/qKhITU1NUVGRy+WWlZWZm5unp6dXVlaam5v3798/JyensLBQS0vL0dGRzWa/ePGCRCINHTqUy+W+evVKVVVVKpUaGRkxmUwrKysCgfDq1auampq+ffsaGBgUFRXZ2tpiQd+YodvS0pKfn4+VwtbW1m5ubm5qatLR0cEKTwgEAqyutY2NjUAgePbsmUAgcHR01NLSysjIGDBgAIvFev78uYyMjJycnIGBQUNDA5vNHjBgwOvXr9va2rCLlpeXm5qaIggSExPz7NmzgwcPfhJ+mMJCiIuDdetg+PAAF5f3d+54Dx5sGhRE3bsX0tMhOhp6SJ//Vydqj/ivv7390y+yw4Mt5LDPW7ZswapHddGff93JLvay/nGbq0cX7aLxv23qY1f8U68yMzPl5OSSkpI+1ZM4eRL18EBRdIyGhpWeXo2bG3r8OCoUomfPojNm/Pfmv/D+z+eYMb9Emhe2dME+79+/PyMjo+Pr3/bnX3ey84l/aqSLv/7FRbto/G+b+tgV/9QrW1vbtra2IUOGfKonUV4OQmFydLTQ2Djw8GFtoRDk5QFBYOhQmDXr8wsGzp/4bSjwjw4uFzBv+Pv3hQYGa6dPD/3ppzFTp0JbG1CpQCSCmRlMnYrrDw4cf4FQCEeOYJHU/DdvTtPpG11dTfv3B5EINDVBQQG+3BsNZwDF8bVCIAAuFxQVgcuF+/dBS4sPMDknZ4C6+kShEPr2BSoVDhzo6VYPPv/g+DGQlATTpwOfD21tUgLhaWbmFH9/nUmTdsyfTyssBGVlIJHA3h6+KEMDrj84vlYwmZCVBSgKBQU5FhYbysq28njh58/T+/QBBQXoNsUhrj84fkhIJNCnj1Qi2bR/v8+DB5uUlIYiCCgpgYYGmJvDJ6vx2qMd0Y+uf6RSaWFhYUNDA7YDjXm7Bw8e3IsVNj9PsU4c3yS4XGhvh6FDF48bl5mTM3XVqmGtrcDnAwA4OsKNG/BparCSyeQeZa1/VH/EYnFMTMzz588pFIpQKHz79m1JSUlDQ0NvEYLxeLwZM2Zs3rzZ0dERlxYc/x+bN4OGBqxcGRsfH/T6td64cdk0Ghw8CD/9BFitKwQBWdlPceXy8vKVK1dqamp2P3Lio5MghULZuHHj7du3T58+bWNj4+rq+uTJE4zJpVdAp9Ojo6Nx5cHxJ7MHysrgzh0+mZzap8/WFSsu2NujWNiXlRV8YmkxMDC4e/eura1t91mw/8GIDAsLmzp1qpGR0ZkzZ9zc3Ai9anR+8vwQHN8QVq+GpiaoqAB7e6irC7lwoYFKHTFiBJw/j2Bqs349rFz5qXvB4XA6kg7/k/6IRKJNmzb9/PPP48aNEwqFkZGRoaGhf8tDjQNHL+D6dXj+HJhM0NPzkZU97ONzQFcX6d8fMjPB2hoAQEXl87iqexRfQuqiFScnJxkZGSxWF8uBw/KHceDoHYjF0NwMamrQ3g4Iwk1NLeDx/LZvN6itTT12TAlFQVYWdHSg2zkvnx8f1R8SiTRt2jQKhYLF0vbp02fRokVy/y04HAeO/4PUVNiyBZ49g5QUmD07KDg4y9p64pIlfo2NUFgI7u4AAA8fwifNY/1vIHx8IScNCws7e/asq6vrxIkTURT18vLCqCFw4Og1/amuhtLSgrt3r6ip3ZNIDo8Y4TdrFlhbw7VrYGoKAGBi8jXfwUfnHx6Pl5mZuX379oEDBwKAm5tbfn5+Xl5eR4IaDhz/EikpwGLB2LFQUQGjRpU8euSVnNzY0vLOzIyipAQUClhZAYfTUw6Qr2v+wfKKO2fhycnJfVclynB8KRw/jhUklebk5Dk4bFm7dkhNzcL58ykKClBfD336gIMDzJsH3cwA/TrnHxqNpqOjExgYuHDhQiUlpbS0tPz8fIyjBAeO/4SaGlBRARRNyspaVls718Zm55o1MG8eJCdDQ8MH6oJTp76JW+nKf7B8+fLffvvt2LFjUqlUSUnpwIED32o1MhxfA4RCIJGAQICmJtDTS710aSGXu2nSpAXPnoGqKkil4O0Nycnf1j11FadAo9G8vb2XLl0KABjVEy4DOP4l6upg+nQ4eRJUVGrI5NjMzGObNz+8dMlk8GAIDgYZGSAQYMwY+B/T0DevP0KhcOvWrdHR0ViJJays55s3b9QwZtS/Q2tra2trq5qaGp1OF4lEDQ0NZDJZQ0MD+6utra1HkUU4vrfJJzcXJJL6mJhB2dkTdHXfTJpEHjwYiERQUPhQWQRB4CsQj96Jv6ZQKEFBQZ1ZCMPCwrqgm6qrq9u6dSvGE3vs2LHLly8nJydzOJzp06c7OTnt379fJBJhRmDHKd0tsYLj20V9PQgEwGAAj4eamxcVF689e3aWj8/h0aPB3x9UVEAohCVL4H8Fkb44ZGRkevSK78GhN2/enDBhwsd4sV68eGFsbPzLL7+oqqquXr06Li5u27ZtFRUVu3bt8vPz09bWXrFixfTp0589ezZs2DAA4PP5ixcvplAoNBrt+PHjuKR9n4iNhTdv4OhRtKDgmoKCz44dJ+rqZufkQFgYsNlAJIKiImzc+ImSEbqPtra24ODg0tJSoVDYoxSGHugPZsh97N+pU6eWlpZOnz6dwWBYWlo6OTktX768urr60KFDRUVF1tbWSkpKrq6uOTk5mP5gFarxKeh7BooCkwlJSQBQkpi4PSHhzOLFU9PSEACQSMDEBAgEIBC+uPIAgJycnK+vL4qibDb7/Pnz3Y9T+6j+iESikJAQjPYS+6WkpKSLqa2pqUlRUXH37t1Tpky5d+/evXv3fvnll4aGhq1bt65cuRJTPARBOiK4EQShUCiU7pWXwPHteQvk5YFAABIJiEQBn3+usHC+j88EPp+orw8AwGDA7Nnw1SyGEQTBsjkpFEqPkgy68l9PmjTJxcUFE30URe/fv9+FuJ87d66lpWXnzp3a2totLS18Pl9PT09NTa21tdXAwCAjI8PFxeX27du//fYbLl3fPxYuhIkTwcurvbW1mUY74u/PLis7d/kyYeLED8HU48Z9HzdK6EIj9fX1BwwYYGdnZ2dnN2DAgJSUlC7qq3l5eTU2Ni5atGjw4MFTpkxZsWLFgQMHDh06FB0dPWrUqPb29jVr1syePdvBwQGXru8cXC6gKKSnA4Vy9o8/7F69Ai73pI0NQU4O2Gzo2/d7uldSl+Yr2vlz10ahjo7O2bNnpVIptlM0Z86cWbNmAQCRSEQQ5MiRIyiK4jtIPwSqq8HJqe3x4wnTpw+cMOF0W9vEoUPppaUAAPr6MGLED6E/Eonk9evX5eXlmDmIIEhDQ0PXqUVEIrFDQwgEQmc7Et/2+SEQHw9jx0J5ea5UOvPVq+3z5s0aNQpaWiAyEpYsAQC4ceM7u2NSF5NPS0vL+/fvO9Rg4cKFCt9CSB+OL2a2jR8PFRXw/v3+a9fCdXWdBg0CNTVQVISnT+HYse/ypruaFthsdn19PTZ1oCgqEonw/FMcH0VSEpiZvb9xIyAtzXDOnIEvXsDLlzB9OrS1gVAIRkY/lv4gCCInJ9enT58Ok0wikRAION8ijr9MO6WlYGUFz5/Dzz+vWLbM1sVlx5kzUFcHDx+CmhqQybB2LfyvSOYPpD+6uroEAoFEIonFYqw6Eu4AwPFnnD4NkZGQnQ3FxacVFMDQcN2wYUChwJAhUFgIAKCoCCdOfK933xV/YlhY2OnTp7W1td3c3BgMRnt7u5ubGx4xgAMAQCr9wKB7/z4wGNL6+mgW69zlyxmmpqCmBlQqjB//IQH7u8ZH7TEymbxr166XL19u2rRJLBZXVFRMmzZN5isItcDxVWD9ekhMBABoa4MBA07s3PlGTy8iMhJycoDDAQBQUQFn5x9XfxAEUVVVdXZ2XrduXVhYGIPBcHNza2pqwiXnx0V7O2Rnfyh8EB0Nr15Bfn67vX14be27mze3jR1r1bcvaGrCj2ShdMV/UF1dff369aVLl7q4uIhEooKCgl7MP0VRtKmpSSgU4mL5zeDdO/D0hLo64HBAReV9a+ulq1cHv3xZwOXuGjaMqq4OCAKTJ4Om5jd6f2KxGJPJ7lMoflR/sPw5T0/Pd+/eTZo0SUZGJjg4mINNzb0BkUh069at6upqXCy/dohEgKKAotDQABIJMJmQnNwyYMDR16+TAwOTHz485OOjlZsLWGzk9u3fboRBS0vLtWvXmExm9/Wnq/y5PXv2rFu3DvO8AYBEIqH1nheSTCZPnTpVRUUFl8+vHdeuAYcDy5ZBamrNuHHqEsnPhw7dTUy0nTgx0tmZpqAAdDowmR94P4jEb4J36m+hpKTk6ekZFRXV/RTUrvzXBgYGBgYGn6ivCIIoKyvj+QtfOyQSyM5GMzOR5cvjo6I2kEj1V69O4fG2BAW5ZWZ+eJvKy4OpKXz7sSlEIlFFRYVMJveC/vytxUUikfDq7T8WeDxQVS1oaYk5eDCqvX2YsbGMsnKgggL4+IC394fa8draEBb27U47/wVd6U9paemZM2fIZLKrq+uLFy8qKipOnTqlqKiIC9UPAWyHh8fjt7dvravTffr0pY8PZdw4+O23DxmjenqAIIAgICPznWUldB9d1S+JjIyUkZGhUCjLli2ztLScO3cuvv/zo4DNhmnTQCAAgSAqJmYcj3d86FAKgwFUKiQkgJUVAMCKFbBiBfzYIV1d5S+0t7dPnjzZxsbm6tWrc+fORVEUN95+FDQ0QFISFBZy6PRSiWT6hAmEZ89g3DhQUIDWVrCxAQDovWKE3+H8g01BVVVVZWVlRCKxtra2oKCg+3XtcHyrauPlBQDw/j2YmUFGRsCpU1I3N0sXF0hIgH79QEYGjI0B4zDA0cX8QyQSjY2NHz16lJSUZGVldfDgQRRF9+/fL4+5KXF8lzh/Hu7fBxYLKithypSKp0/TX7+Ou3GDWFYGRkZAoQCJBGFhgKeB/aP+EAiE/v37jxo1SklJSSwWy8vLX7lyBR+v7xPx8UAkwqhR8OYNLFgAb940VFSsOHkyg8U6ZmREVlcHXV24cwcTCzAzwwesW/6D2NjY1tbWjIyMQ4cOycnJPX78uL29HR+y7xCBgXDyJLS3S7W1U1RURixfbr5/P2nUqJeDB08dNAgUFEBODiws8HHqwfwD/yukiuf8fP9obARDQ3ZFxbibN6saGzO9vVuUlBj+/mRbW5CXB5y74t/5D8RiMZPJZLFYXC5XIBB0QV7Vedbq+i+8AtdXhAsXoKkJeDwwNGxQVp7j46NnYXHE31+lqMiYSMTYBAF/e/47/UEQRCKR7N2799dff01NTf3pp5/y8/O7oNFhs9keHh6DBg1as2YNAKSkpIwYMcLe3v7y5cutra3z5s1zdXXdunVrZw84ng3+JSGRwJo1kJsLmZm5FMrY1NQBDQ1nfX1nrlkDKSkfgkF1dEBL68fShx7K5Ef1gUqlBgUF1dTUYDVIsBDSLlq/ceOGubl5bGysjo7OqlWrgoODt2/f7uDgYG1tjSBIv379IiMjHR0dp02bZm9vDwASiSQ7O5tMJhMIBGuMkxLHZwCHA9XVYGEBRUWgrw+VlW0czszExKPu7iYEgrysLMjKgkDwQW2GD/9BAgvEYnFZWRmXy+VwOFwutxfyF1AUjYiIcHFxsbW1PXLkyJ/43P6KCRMmrFu37vr16wCgoqKSm5t769YtPz+/uLi4wsJCc3NzABg5cmTy/wqMSSSSBw8ePHr06OnTp7hUfxKgKHSY3GFhkJwMKAp5ebBsGXA4kJoKI0Y8iI4eERBwMT5+rJOTSVsb9OkDAODoCFjc8Nq1MGrUjzBUQqEwNTUVk8ampqbuz0IfnX8EAkFWVtbFixdNTEycnJw2btzYdUMUCoXP57e0tOjp6eXl5TU1Nbm6uiII4ubmtmbNGkyhZWRkOhiwKBTK1q1bcSH/hKiqghMn4MgRaG+HW7fg/XsYOBBycoDLBbEYkpLu6+lFpKXdMTHRtLGB16+hoeHDxk5o6I9mtsnIyMyZMwdTpNOnT3efp60rPaPRaCoqKtrdi6sNDg5OSEhYunSpoqJifn6+ra0tRoTNYrG0tLSwPLmEhIS+newBgUCAC/knRGEhhIcDjwccDpBIwOUCgQD5+eDoCHz+pJs349+987W01Bw4ECQSsLYGGg3k5AAA9PR+WJ8bn8/Hliq9oD/V1dWDBw/W0tKqr6/X0dFRVlZubGz82MEzZ87cuXOnsbFxS0vLsmXLJkyYMG7cOB0dnePHj48ePfry5cumpqaysrKjfgx74MuKALS1AQDU1YGSEtTVAZuN0mgoiYSiKNTUVBsbjx89eoyV1bGzZx2YTKDTAUHAzg5CQoBOx8evd/wHNBrt3LlzmNMZq3+KTXMfO97Q0LC8vLy+vl5LSwsAli5dOm3aNCqVip2SlJTU2Nio+c1mxn9LiIyE+/fh2jVoawM7O2hoEBGJMSxWfnq6ekhI0p07T169OjV37mQ+H4hEYLFAWRkAgE7Hd0h7U3+wJUqP8kMRBNHqZDcrKSl1fCYSibjyfFq8eAFmZqCmBmlp8PYtAJSzWFdzc73fvAnIyOAoKdUkJtKrq+WHD7/Vt689iwUDBgCJBMrKwGAAvpHwKfQHx7cELy/w94f584FOBzm55ubmvU+eKAwcODssbAaNNjsqSuzrSzQ2JgwfThUIYOdOmD8fiETYvx8w3hwc/wpdvXgqKyuXLVvW3t6+b9++Pn367Ny5E1/xf6WQSoFIhNRUNC8PTExqjIyG2ttXoqifp+e9ESMW6uvTdHTkSCR6bi7VwQFUVKCg4IO1Zm4OncwEHL2mPwKB4OjRo/r6+iKRKCAg4N69e4mJiZWVlfiQfUWQSGDvXuDxoLgY3N0lNTXXi4qsfv3VIyXlyKxZDyZO1B0+HEpLEWwXSFERsrNBWRlUVWHQoA9bPTg+kf5gzIZr1669ffv2mDFjXFxcrK2ta2tr8SH7ilBSAocOQX4+68mT5wD+5eWPAwPXzJt3aPTosW1toKgI8vLw6hUwGAAAqqrQrx8AgJ4eXLiAB7Z92vUPlUqVl5ffsmVLVFRUZGRkRkbG06dP/fz88CH7ipCeDuPGNcbE/BQbm5meftTGZv7IkbKenlBbC/Pnw9SpAAAkEhgaAgDMnw9z5wIA0Omgo4MP3qedfygUyvr167W0tAIDA0eNGnX9+vVdu3YxsDcZji+LsjKIjQUASEpKd3d3O3lyMoHw4NWrFYMGyfL5QKWCoSG0toKuLgDAuHEfil2bmOCpb59v/hEIBL///runpycWurZ//358sL4WV0FMDDx4AB4eJQ8fXpeTO+jhMVFfH5ycICoKystBVhYUFMDZGfT0AACCgz+ciDvZPqf+YPzxSUlJHA6nI2nHxsaGTCb31rXxUkI9AIpCbS0oKwOdDqWlQKej9fULRKIlFMrE2lrAYnA0NaGiAhQUgEKBmBjAh7fnkJGR6VHCaFf7P1wu9/r16y9fvsT0RyQSYXSKvdJRHo83b948X19fLJ0Bxz+gtRXWroUNG2DQIBAKW3V05g0dutDdfZGfH7i4fAj3tLQEIvFD6g6uPD1HZWWlr6+voqJi98vFd3WciorKokWLBg0a9Cn6SqfTT58+jde1/2drTSwGCgXa2uDNG6irA4mkPiNjmbLyhuHDh/XvDyQSqKqCri4QiTBmDIwfj4/ZvwaDwTh//vy5c+e6z9P2Uf8BgUAwMzMTiUTNzc0A0NzcXFlZ2bv1t/Fq3v+MBw9g/PgPmTwaGvXl5Q9fvlybmTlfInFjMhEnJyAQYOnSDx42MhmPxPnMYvnR1z+ZTDY1NT158qRIJLKysiosLMQ8Cjj/9edAfT3IyoK8PBQXfyj5JhJxlJS8Tp0icTiqs2d7Dh4Mu3aBpSXIyMCCBfiA9SJ6RLL7Uf0RiUT37t2bMmWKurr6rFmz0tPTJRIJTp74mXD4MOjqgp/fh90bPj/34cNjSUl9HByCpNI+e/dCUhLQaPgi54uD8HHDWyqRSNLwF94AACAASURBVPT19Z2dnfX09PT19Q0NDXEuq0+LDnKi9HR49AhEokYOp9HWdquvr42fn+XatdeGDNF7/15BVxcUFWHAAHzAvl79AYD29vb79+9fv36dxWLduHEjMjKSz+fjQ/YJsW4dVFWBQACWllBRUcVieUZETEhOVlZQkEyevGH3bsjJ+VA7xM0Nrl3DY3C+Xv0hkUgjRowQCoUFBQXLly/Pzc19/fo1Hn/9CSEWQ0gIPHwIxcUwYMBdMnn10qWWRkaR69f/TCZ/yKymUD4Q4mBJozi+Zv356aefDh48eODAgaVLl/J4vLi4uC7oEXH8S7S2wuLFAAAlJeDsDKWlpfn5m2/dOl1UxKitPT1jhtno0RAd/UFtliyBKVPwMft60NX2i1AofPbsWWBgYGJioqur68uXL1Xxki+9jrg4ePIEeDyorJQOHnwkOvpJbOygqVPv9e0rkZUFZWXQ0YGSEujfHwBg5Eh8wL6N+UcoFPr4+MybN2/q1KlPnz7V19fX0NDAx6vXUFkJ5eUAAGlpYGoK9fWQlvZbRUVyc/PpQYN2rlkDWlrER49ASwsQBBwdwdgYH7Nvaf4hkUgrVqwwMDDIysrCWBFbW1v74ElXvQKpFDZvBgMDOHAAkpNh9WpgMq89fnwNQRLs7D5shmprQ0XFhxjq2FjAK2d+W/pDIBDs7e3t7e1bWlry8/OvXbs2cOBA3IT7T+DzITcX+vYFOh0yM8HODgCkjY2bHzyoPXasvrZ2W0QEHDsGDQ1AowGDASj6oUaiigo+eN+Y/QYATCbzzJkzZ8+etbW1DQoKysjIUMa4jv4OAoHg4MGDHh4e0dHR2C9isfjgwYNtbW0CgSAoKGjy5MkPHz78oQe7vh4WL4bSUuDxQEYGlJXjw8NN6+tRNtvP0PCht7f78OGgoABSKdDp4OAA8fG4gH6r849IJDpx4kReXh6dTp8yZcr9+/flMBfqR5CUlJScnHzkyBFbW9va2loVFZVjx475+/t7eXk9f/48Kyvrl19+mTx58osXLwwNDQFAIpFkZmZi/PF2dnY/xGBjdRHJZHF2dqqWVvXTpxcvXcq/c4diZQWenjB3Lkgk0L//B0IPKhVwWv3PBZFIVFJSwuFw2tvb29vbux/C01X9bYFAsG3bNjs7O61usCHr6elt27bt9evXxsbGsrKyT548AQAzMzOpVJqXlzdp0qSBAwfOnTv30aNHy5YtAwCpVJqdnU2lUkkk0neuPyIR1NWBnh60tYGt7dN792LfvCknkQxFohMDBlD69wc2G0pLQV39QzBo72VY4egmJBJJcXFxY2Mjj8djs9m9oD8AIBQKX716xWQyRSJRRkYGm80ePHjwxxgVzc3NeTxeTU0Nm83Oyso6efLkxo0bw8LCysvLRSIRRmivrKzcUQGSTCYvXbr0hwgIOncOYmPh9m1gsSLI5Ic5OU7FxduDg2VfviSmp4OeHhQWgoICyMgAgsDHLWQcnw40Gs3DwwMA2tvbQ0NDu0+B3VX+Ap1O/+233/z9/Y2MjLy9vVevXt2GESv/HbZv3x4TEzNlyhQPD4/09HSpVHrixIn6+vqgoCAdHR02mw0AWVlZBlhhjP8tkL7nZ1JdDU1NIBBAczMUF4sBdty6tTgpaXb//j7v3/extiY1NiLYVEOlfigNj+MrmIh6VCLxo/MPhULZv39/92kPRo4cuWPHjtLS0pSUlL179/r4+ADA+PHjIyMj8/Pzjx49WldXV1dXN3HixB/lUYwcCbNng68vyMq2yMn9furUw/v3c2/csLh//wOhlEgEBgaAIKCvD9evQ0+oknF8Jei1dKthw4adOHHC3Nw8Li6uw8d95swZCoVib2+/detWCwuLmJiY799gCwoCAHj7FhwcoLQUUDSvpmZKUxO7pCTJ3d10/Hh49AiGDAEA0NeHKVMARQFBgErF49m+RfRm+rStra2trW3nXzqsNSsrKysrq+9/OEtLYeNGWLwYysrgp5/e7dr1/O7d9WfOzBEKt9nbE0tLP0w7AwcCAHh6AomEqw2uPzjgwzSSkgKGhpCdLWEyEyoqpr59a7l69e1Dh9ypVCQ6GmbMAADIzf2Qd9DlfgAOXH9+GOTnQ0kJeHhAWhps2VIZE7O+qKiotjbdwMD00CGCigpQqbBqFQQEYCtLfMDw9Q8OgIYGYLMBAPbsgdOnAQCqqwNev14dEWHW2Hg3Otrcxobw9Cno6Hzgx7GxwccM1x8cAAAgkcCWLXD3LgBAZiaYmADAvYaGlyUlm/r3PzhzJkNdHbS04PVrMDAAVVV48AAfM9x++4G1BUGAQID0dPj1VzhxAmRkICEB7OyguVmkp0cyNS2+cmVdbu6b8nKarS3IygKNBhYW8OYNUKlApX7wuf0Y4AoljW1CrkhKJBCkUqkEAEGlJCJJKkXlaUR1eTKZ2GsvbhQFFJVKJGhlZRO/XVhRw5IhUVXV6WIpoqRMRRAyg6GM68+XRmAg9OkD3t7Q0ADPn0N7OxCJQKOBvPzrR49OiERqr18zX750s7SkycgAhwNCIQBA//4fPvwY4AmlIomU2S5mcYUMJZpQLBWIxGIpaqunkFHZCkKBrhKdTCa+qeGoylP60EgyFOL/a+/K46qqtv/aZ7oz916GywwyowgIojnPA4ZDgxOaJWYa1cvSVy/NMssyLSccKjOHBn2ZpJkaiqkooTnhAIjMk3DhcoE7j+fs3x+HfL0yH/jLsrzfD38An3POvvvc/d1r77XX+i6avHP3o9PJWS22yoqmytKmgcOi27RGb2+5r03OYifnBABnU50dA9tU2+zuI/f1VVI0SdOkiz9/ODgOamqguhqefho0GhCLweFwHDxYPXjwR3v2NNbXT5ox41BW1iypdNgbb4DTCbGxwKvkDRkCQ4a0u+b+7rjRZqcpVK0xMyRqM9oECIQ0KaEprdHWZnJQGIkYSiqgmo12wFxjm9UhYcQCCgB85XfgUMEmk0PXZiopqOsS7hUW7Xf+fFVZqbah/oreYHU47LSAlsmEKi+3uLgA/wCFrs2sbTJiAickhNAM+Tfhz19GP76tDXx9644fF2u17pWV0Levsb6+78KFWK9fPHHiwJoa/zlzxsbF0fPmgY8P0DR89hlIJP+5/e9OHg5DTavVamftDqeT5YQ0GeolySlvOVjaWKI3m1nWiTGFkJKhe3jJxkX6JPi5tVnsxQ3GuACZRERXaS1dPESdarH8eh0tFHl4Snv0Dvtm39ndX+a36YxWi411ciwiAGOCw5gAEtNiOeUul06e0nvaY/0unL5+7nRxZLcATy/5bcbk76OfeLdht9uXL18+ceLEiIiIe3dofPYZzJgBra16p3Om1Vo/cKB/a+tZnS7y/PmX7fakrKyuublw5gxIpbRYDA0NwMt530+J7higpNGMEJhtTpuT85Axe642LD9TbsROIEkADD9FkyGwnGvRf1xYq2IE742IjlfJrtbpFRI60tetrs3mLxd0ZNxyHK6uaq69YfD3Jusduvkv7Kos1zAMQojgMCJIkuAQJjEQCAHGCJuMdqup9d0VB3ftOvP+2sdoCovEQpPRJpHeYu5uamrasGGDTqfreJTMn+Z/I0ly1KhRKpXqXhwUfJhsSQnMnAlqNej1b27evJDjtmdkKGNjKzMyzi1cOKNHj65JSVBWBj4+fH/Azw9EIrifgDFUNlukAkKrt9S22VgCz9x7+dVjhUbEAUEA/g95eKYBBiCJJtb+xDf5q85WuUlojOFyTZvZzqn1HZJG0zTpy0vUEZHe18o1j0xYW13dJGAIwIjjOIQRYEAYI8CAOUCAECYAscAKabLuRsv0yRuuXVMXFdZfza9UN+h//XCZTDZy5EgPD4/fIf76D+BPfHz8PSSocFMasq6uXVCqqQmioyE39+UdO0QTJ/adMqV3QcH6IUPce/aE3bvbD3MslnZlHC8vGD0ahML7ij9VWounlKnUmAiE9A77jK8vFegMw6N9iFuaEnxTYBWBkPrsau3c7wr0Nqe7mG42WDgMav3/8LVUFDc0Nur8A5Qnjhe/+drXFEWSBIE5LBBQCoWIxU43udBdJQGMZQqRQEQ6HaxUzggZmsNAUQQG/ObrX9fWNovFwutX64yGXzJWJBIlJiZKJJKOh2D/mfsflmXvofolTzwBDz0Eqalw9Srs2QN6PVRXw+zZB5cutRoMyzMzyZISSE/3LiwEloWTJ9triaanA0kCxhAcDEuW3FchOXqLkySI+lZrQrD8+zLtU99d3T0x8Z3cUg8hHe0mLNJZAAFgANROHk8R5SVkrreZOf4/ArJIb3zuaNHHyd1ZFjfqrXKxwOrghPSt53Rdm0kko3UGs2+Q57HjhXaHg6FJjJDObN62aUbh1TqhgG7UGGLj/A58e3n0iFg7y+bmXu/XJ1JvtHy67ZTJbCcQ7hLmFdM9IDzc5/vv8gvyax7oH4b+25Pe2Zog9/f5KcfBgQPQ0gIAcO5c+2FoZiY88YTt1Kl9WVmTP/98YU3NqokTSR8fCAsDmQz8/cHdHSgK+vQBAEhKgoQEQAhIEtzc7p/yIRzGZgenNdqa9NYrdYZXjl2P83N7KPPCoED3KxrD7IQgDzFD0KRIQANFAkUCTbqJBSoRgwRk+38oEmiqwmrbeKUWAGx2zlNCaYzO315T2yqq9eHRfl7e8k+2znnuheEsYAAMDiyTMVlHrkRGe9+o01aVN7q7y3wDFT6+bu7uUplC6OEuIikCMI6ND9q565moKF+SRH0GRTmxvVVnumVbCCGj0Th9+vRx48bt2rXrNqS6v/ljs8GcOXDuHOh0EBnJSiTlpaUNjY3ZQUELlyyZ/+23jSrVJ9270w4H+PtD795w5QoAgFQKP/4I97canoPFNgd3XW0M8hSdUeuKW02Ipp5JCDpQ3fxcUnCBxtDfT0EShJ+EAYQAEYAIo5OL8ZZxHAKC+M8PRX5d0cwJUMENQ8ENA0WA1XGLvYfFbNPpTAEqqUzWvsOcM2f4+o0z3NwEHEZGgw1hQq+3SmRCmVJoNFjOnqkoLlYD4GPZZTV1bcDhgCDlJ1tnMz85r2tsSCoWNTUYbrlSYxhm9+7d48aNW7duXXFxsdlsvhfXb38aWlpAp4OQELhxAxQKMBohL29+eXk0QRybNs1ks5nKy0fbbBdeeMHtjTfI7t3bDRRAu7IHgKv+u97qtDtZiYBUiuiFeRUgF55tNtWY7Y908bSw3GPx/m+drnJQpIdE6KAoo511ArgLqUazY3y46pvKZqCI9o0QYCDwqjM1m0ZG1bZZhQzh4LDwVrtlsVjoH+Tu5H66FWDAgKiXF45749XMjRu+b2zUf7jphFBAajTmosKa0uuNiCIcdqfdzhUVM5pm47KVU6ifDk+3XNF8Xtz8RX+fkqLaiGjvX3vb7Hb7I488MnPmzPLy8ri4uNsctNxP/DEYQCwGkoSLF2HTJvj8czh7tm7o0K2HDllZ9qynp6WuLiMkpOjBB1UREdEvvkjwWaLh4S75tVvsRizOmhaLhCG/r2yxMyQwFACoMWyt0bnTBra4aXFiYMNV3IohxkOyfljk/vLmKZEqikAmB3tkp8GCAJwccBhIEgREtsYgFJDNBntZoyncSwKCX3r5mjWG1ma92Esxdv+1N/p3SenSvs8cMyb+6tXaf392BgFUGjWIAOzEFE3eMLQCABCIQJy2GULCvHvEB/G3bC/UrD7XwBGohiX8/JQtLWYvL9mv7c/nn38+a9as8vJyjUbjdDp/S/bj775+wxh4zXu7HVJT4eJFcDqhrAyqqsDprDl79vm8vM9zctp27z567Ni6Z55RNTYOGTGiW+/eRFlZuzPg0UdhwQIXYX7p++Ggm68kMUj+VanW013S3UMSoRT5ygQKMY1pkmTo1y/d+FdS0Khg95d6BYcqRE/G+vnJBCoJE6IQTY5SxXnIPhoedWZq4op+IVKGEsoEF9XGKG+JgCT0Vgf3yzUVNposnl7yS1pjdavtmeyyrytbbi7yXn55bNce/gghkgISIZIiOIQRSZAIEQQHCDkd3ISHe4hEFAbIqm5b/mOjGQFwuExv9/SSWQ3WW3q2EhMTP/zww/Dw8Nzc3HvU//ZHoKQEli6FL74ArRYaGqCgAHr2BJ3OERhYf/XqlMzM2TNmfB0WBgcOgEgEYjFYLCCRgEAAnp7g5wcAMGOGiy2/htPJFbdaI70kTic7QC6mEBoZrMhvNrVZnCzGYXJhNy/xDaN9qLvUk6YQgNvPoma2DI2giHYH93cVLYOVEhag1WSPkgt1Focc0b/y8uCWJn1IsMfVZhtJAnCw8NgN0VB6TGi73chYO2Pu3K2lhY0EA4hAwGKCAExgjEnAXGiY1+iRcSRJHqpqXXC8zuxgSQQsiZodbFVpGyEiAsHrV7MuVqlUgwcPLi4ufvTRR29Tc/7vzp+zZ+HoUaiuBpsNwsOhqSk/P1+Unb3dYmlavvyV+PiHXnkF5s5tPwNVqSA2FtzcAADGjXPJF97OfwDYXymSieiRQUqJjCYBhgYrZ8R4ay0OAoGUIt0EFO9V++Zyo7eYUbn9Z/1zkzwAEKMUhrkJaIogrTjAXSg0EL/2YbIsxwgYqVSsr9AQBEkgZHU6nzle9ZxOtSDBGwA8PKQZGTMWL95z9nSFgCGAJDgAxCEEHCD88svJ/l083v2xfktBs83JEohEBMdhsDs5FrPeHreo50uS5I8//tjc3JyYmHj69OlRo0b9FoX+FvyxWNoP/o8cgQceALkcWluBrzlVWwvJyXDxIvj5ndHpVOfOTVy9OkCp/GdUVN+ePT1JEmQyaGuDpCQAgCFDoHfv9ke98047kVy4FWgCRXpLAGBn7o2uAW42B5vxXfWgMLnNgc/W6AeEKQwW54UbhkGhihK1SSGmitSmS/XG0mbz0HBlq8lxucE4KFzZ2Ga7UG8YEKq80WoeEu4uEVASAVXbavvVaCY4lrNa7Uo3AXBOTNIIEYhlN55vKmmzrhgYpKCQn59y1fup77z97cED+QKGRhhZbHaFuyhj3WOqbqGPHag402AEzNGIxMAhDhEEEhLYx9e9vKwxIPiX9sfhcKSkpOTl5ZEkOW7cOOa3U4Z/T/5gjH9+JMqyLEKIV07kOA5j/HuK7/B+R4TAYICRIyErCxQKmD8f0tPh2WchPR3q6+HkSU1dncbHx3bsWMrevTaTCRjmykMP+c+aBdevQ2YmpKUBQhATAw88AAAgFP4ngMBFntuCD50GgFhvyb4CjYAkolXiHefUAoqI9ZZ8dk4tE5ARnsJ/5zfKBOTa76vSBwd9V6hxY8iPTtfLhWS0p2jr6XovCR3sIdx7uUlEEk/2DQQAJ4cdLPuzY9f29ZtAIqAY4gFfqZXFDIkwYA4wwuyh4tZTNYYl/QPHh7jJlJL3VqU+PrPf51+ctphtQwZFjx3Xc191ywsHi+oNDgASAXAEB5jAmGUQ5SsUsKw9Itrn170TCARbtmzp2rVrSUlJXV1dYmLiXbc/Dofju+++U6vVMTEx/fv3r6ioyMvLczqdo0aN8vT0zM7ObmhoiI+P78VLz/C+i//PaePly6DXw6BBUFEBZWVQVwdCIdjtUFoKANDaCgEBFcePb6mqyjx3LtjhWBMaOnzfvsp58/z1eoiLA4MBcnPhvfcAAN5/38WHTvsPnO1Hikmhyso2G5CousU6MERhc7IVzZa+wTIWQ3GjuXeQ3OpwljdbdpxXT0vy3X1RPSRMaXWyVVrLwDC5xc6VNpn7h8irm83dfCQAYHVwQor4OXkAgKIIDw9p9Y3W7pEBFIEQcLg9mRETJDLZnM8eLl0mFU+M8ojzFEYH+Kx4dyp/Y4XeHuAmmZsYsLNAc73FzGEEHEaAMRAkgYMESNNojPAK+C9/GkHwg3nYsGGffvppSkrKgQMHbjPv/27+t6ysrLVr13p5eaWmptbW1r7++us1NTWVlZWzZs3Ky8vbuXMnwzBz5sxRq9U/TSpcVVVVVVVVdXV1R9uw2WDbNmhtBacTduyAXbuAZaGkBKKicFsbHD0Kjz8OFotDrTYNGPAeQQycONHg43No794jQ4ZM6dXLMyioF0KQnw9iMQQEAMYQGeliwp1BJqR1FicAjI31ajTY9BanmCQaDfYWs0NMEc1GR6vJIWNIrdGmtzhlFKkx2A4UaMbHqhgCawx2/mK9xSmjiXqdLchd6K8UAoCD5WQimvjv6DmEECOg6mt1UgpGdJFZWRYBQoCdACQAAJLStJ11fF+lOV2v01vtFg4frTW+dKo27buKqd+Uzj9aWdhs5jnJRxQRQIQqBIE0GPRWiZjh10oNDQ1VVVWVlZU2m43juOjo6Hfeead///7Lly8X/nZY4+9mfyIjIzMyMiiK4oV/R4wY8dhjj6nV6g0bNhQUFIwdOzY1NTUvLy8rK2vmzJk8xT/77DOGYQQCwUsvvXS7R//4I2zaBDt2QFsbbNgASUkQEQFaLcjlYLNBYWFDeLivzbYoI8MzLm7bzp3U6dOtzc1gNH4yduygnj3FAwbAlCnw6KMAAAEB0NwMABAVBfv3u7TX7hhyEdVicshFlIeUGRCm3HymPsRdVKgx0STh78YUa60EgQLdmJIbZoYi/GRM0Q2jkCRqWiyvJoe9m13ZaHFoNHYBQ/lKqIJa/bZpMTRJYAxmBye/VWq3SCQICvHgEJqf6H+44hpNYgIIBAQLLAIMQIop4uFu3r4S+vPi1vycmhYLZ7JzGDgSIRlNcYAxBgSIA0CAdXbrvF5BXh4iUhDEp9NZrdbDhw/fuHGD14/nrdC33367bt26gQMHvvLKK791hPq78Sc8PLysrOzVV1/ljeDMmTOPHDkya9asjRs3FhYW+vn58dfcVNAWCARLliy53Y6I49pdMdeuwTffAABoNGA0cuXlEBLiNBgqgoKosrIfdux4rbq696VLXgCePj7vJydbHY761FTq2LFkpxNkMt4qA19CYvBg4NONCAIGDnTR4I5BkYjF2ObkBBSxcHRoi9Fe3Wbr28XHznI1LZYHgtycGCo1pqQ4LxuLa5otqfEqs51VG+x7Lze9My5897kGKwYry1VrLf8aGjyqqweBwObkOA5LBbcYEhIJ0y0mAGMcqRSNDVN+X21ABMLAEkBiAARcq5VbnVfLYUyTBACBABMIABAHiAIOAeKD5RAgloMB/ooRAXIA8BExPz1fwk/rP9ePP3/+/MGDB2+f5fm78efAgQO+vr579ux59NFH9+7da7PZcnJy8vLygoKCli9fztPm4sWL48eP/9kZgvN2/Dl2DLy8ID4e6uogIAA0GtBoID5+/xdf1NXUHL1woer6dfWKFYO6d3/72WfJq1enCYXwz3/C669DeTnMmQPDh8PDD8PChQAA/fq1S07zVuj+SKi+u/whkIQhjTZWQBGeUmZApMfpY9WlWi2BQCygyrStiEBihqwsaf+l4noLRSABTV7XWIoajM8OCdpwslZtdHhK6LggNw8pAwBak/OW5LmJ4uIbHgrJK338C7WlGjMmAGHgCAAOAAEiCUQDwkBicAIQGIAAIABYIAAwAQgBdnDYS0x9PCrkWmGdr59SrpD8YiDw+vEURZ05c8ZkMj399NPx8fHPPPPMb7ngfjf+kCSZnp7et2/fysrKpKSkUaNG9erVa9GiRUqlMi0tbfHixWfOnKmvr09JSbnttpSFzZshPR0AYOlS8PGBr74Cmy3P17dfXV1lZeWbeXm5bW3JVVWTnnkm6NIlaVRUwsSJkJgIP/4IsbHAcWAygc0GDAORkSAStYfeLF/+X5k5LvL8HlCK6TKN2U1E0QSa3S9ATKLl2VUcAQ4WyxgSCORkORlDYYSdTqwUkiyHnRynFJI6B7cup2ZImPJsZdvKh6P6hisBQK23sxynFN8ugYqiCKPJpmDxG30C/3G8isWYAAoDhwAwYH55hoElf1qnYQAWMAWYA8AALEZyIWwYHaqr02i1elpAyxWS33KG9enTJzc398UXX9y7d6/dbr/r/BkzZkxYWJhWq33llVf8/f1Pnz7NG0GKooKDg9evX69Wq6Ojo2Uy2e2e0tQE8+ZBXBz06wc+Pri6+tyJE2UHDnyg1Xq8/vq1goJJkyfvb2oKJQjB3Lnw0UewezcMGAAOB1y5Ao89BhwHSiXYbO21dL76qn3ZFhrqGu6/OxCCUC/R5TpTYqAUIUh9wB8Q+vKimiKRQkBRJNKYHe5iBjBusTg8xbSDxVaWk1AEh6DV5FDrbKsmde3dRQ4AJjtnsDojVOL/tUfwLSioAgwpET5BMtHUA0VtDidDEACIA0wA5gARABhI3iYhICjgWMAEIBZzCiFsfzBKrNYWXquL6xkaEqq6/USakJCwYMGCYcOG/RHnPyRJdu3a9eafof89ZMPCwsI6UoG9pAT69rXm56PGRjR6dPabbz7+6KNeEsnJf/0r//jxIX36CFavhkcegZoaUCohKAhqaiA4GNRqkEohMhIIAlJSgOPaBabvk7KQfx5IhBICpGerDUlBMgLB9Af8BoQrn/ny2mdFjSzGMhFlMDkJAklFlM5op0jEMKTVzgpIIq2X79sTIhRiGgAsDq621RrpLe4AY1FsbEiL1pBzoqj3A2GnUuNmHS693OxAmON93hRwCCgnOAggebvEAQZAGLh4T/H2B0N0da1NBmPvfhF+Af87Jvj06dPz5s2LjIz8y8TvcBxXM2HC5LfeCg4M9Ojd+3Bzc/aLLyb27QuenqPPnWuPDKiqgp49AQBiY2H4cAAAb2+YMwc8PIAgXJz5461QvJ+kTGP1ldMyARnsITr4TOKRQk2zyVHbYm2zOO0sa3JgGY0ENOkvFyKE+4S59whoX4a0mJ06syPEQ9jxgxShgPFSSI8eLoyIUP57dNjRBuNXpdryVpvG4jTagcNO3qNAIRDSWCWmIhTCRyPckwNl+ZdrtRqDj6+8I+QBgD59+nz88cddunR5+eWXf8uFfW/xp9DX1hIgVAAAESNJREFUd+o//tFDLn88OJh68MFnNZpYoxF8fKBLF8jKgsWL22mTlgYAEB3dfgCKEMyd6xrKfxYENBHszrSYnY16e4iHiCRgVIwXAFRozG0WR4POJqBJh4NVSBgRTfQIbA/sMDu4+jarTEgHKoVUZ1QUxVJBVHd/L39lXXVjbk5JTBdl3yRVC8nUmpxlWlODxWGzEwKG9ZOIwhTCAAklt1gdNtv3Rws4gJ69wj1Vsg42lJWV5evrW19f73A4/hr8iQoNXbls2dA9e8QqFQwfDnV1sHEjLFoEnp7Q1AQxMQAA69a1B9fQdLtXzYU/nUIU4SNjLEK2usWCMQS5C2kSNeptAprwdmPiA9yOFjZLGEJjsAEGi4OrabUKKSJAIfwttYP/sW6kSC8vmYeHtKG+raSoQaazGnUWlbc0QMh4BXk0txgFAhFrs9eW31B08ykqUctlgriEEHcvWaf0R2fMmJGXl+fn5/cHxb/9/8EwTEpKCvz732AyAUIwdCh88kl7Sd1du9qXbQqFa7zemws5MUOGeorsLHdDZzPZWF+FSEQTbRan3oaj/GQCipCJmEK1US6kI1Qi4v/tBSUI5B+g9A9Q2iy2mprmitJGT5X0yqWapKQws86mM1iUMsZLpegrlzAM3anWhELhtm3bVq1a5ePj8wAfG/mX4E/799C9O4jFIBBARASsW9dez33KFNcY/UuAIYku7kIA4DBuMzukQspsZ4U0KWWQm4i+K9ZPJIiI8o+I8rfZ7D17hTtZjqRIp4PjrY1A0OlGbTbbzJkzZTLZ999/36dPn9uoRN17+acUBenpMH06HyJ4Izy8pKTkbrd5kFfeuZvAGB+5+yVMrl69qtFo7nYre/fu7ZB9QMhdwvjLmUAF4ydnOkUeo9H4448/dppIAgYQQVEUAtSRpdqpU6dstlvoNtI0nZ2dvWnTpvr6+qKior+a/o6Hx011m8qSksuXL9/tBj/99NM/gD9ffvnl3W7lzJkzDQ0Nd7uVjz/++M7eQMcv1ul02dnZd7sj+/fvt1pvkb/tdDpjY2PffPPNpUuXTp069V60PxjjjrxQjuOcTmennrx58+aGhoZOfVv2zlcZefLJJzvb2c62UlpaumvXrk7d4nQ6O6sAuHjxYgcvEXE3X9e8efM6pct+B997dnZ2Z00W38TPhwr/TXEc5+vrO2DAgB49egQHB98m0eZP2/9IpdKcnByBQHCbgY4QKioqam1tzcnJ6aAkMU3TV65cUSqVXl5eHaEQQkgoFNrt9pycnA5SGiEkEAjKysrOnDljtVo7eAvHcRaL5fTp0zabrYPcrqmpKSws9Pf37yAlSJKsra0ViUQ6na7jHSkuLj558iRJkh1/XRzHdbwjfA7l9evXO/W6tFptU1NTbm5uB7lNEMT58+elUqnFYung62UYprW1NTc396ZgL0LIZDI5HI6Oy+KiTs3TvyNKS0sLCwv/ZwodSZIIoU5NRQzDOByOjvcLYywWiy0WS6fsiVgsNpvNnZpThUKhxWLp+C0EQZAk2SnjQNM0y7Idlz/HGItEoluuYX7fvotEok7dghCiabpTho6iKD4DujP7JYHdbv/FUAkODo7vcKXaP40/LrjwN4CrfrALLrj444ILLv644IKLPy644OKPCy644OKPCy7cPZBvvPHGvfNp1q9fz4e7nj17dvXq1WVlZb179waAzMzMzZs3A0B4ePiVK1fee+89jUYTGxuL7iiG99NPP+Ud/Nu2bfv000/d3d09PDy2bduWn59/+fLl0tLSmJiYPXv2bNu2LSIiQnmz5k9ncOzYMZFIJJFINBrNypUrT5482a9fP4TQd999t3HjRqvVGh0drdFo1qxZk5eX17t37ztTZl22bNmgQYMAIDs7e+PGjU1NTfHx8S0tLe+9996xY8ciIiLc3Nx++OGHjIwMAAgNDe3s62publ62bNnFixdjY2MFAsHu3bt37NgRFhamVCpv3Ljx9ttvq9Xq+Ph4p9O5Y8eO3bt3JyQkiMXiTjXBcdzRo0c3bNhgs9mio6ObmppWr1599uzZPn368GeDH3zwQc+ePRFCZWVlq1atKi4uTkpK6mxHLBbLqlWrDh065OPjo1KpcnJy1q9fz3FcRESE2Wx+++23T548GR0dLZVK+bA3T09PX1/fjrRyr9ifpqam+fPnL1q0CACuX7/+7rvvjhs3rrKycsmSJbt37/7yyy+nTZu2YcOGffv2zZw5c8yYMcePH7+DcLLs7Ozk5GQ+WnT9+vX79u2bMGHC0qVLq6urtVqtXq8/evRoZmbm2bNn169fn5iYOGnSpM42wXFcRkbG448/3tzczHHcokWLvLy8pFLpuHHjKioqNm/ePGXKlBUrVhw8eHDZsmXe3t52u30WX664MygpKZk1a9a7774LADk5OWvXrp0yZcqJEyd27Ngxbdo0X1/f6OjopUuXFhQUpKWljRgxYuPGjfn5+Z1tZcCAAVFRUW1tbatXr87KytqwYUOPHj2ee+45AEhLS+vevfuuXbsOHjx48ODBY8eOhYSEjBkzprNNlJWVvfPOOykpKW+99db+/ftXrFjh6emp0Wjmzp1rNptffvnlf/zjHxhjg8Hw3HPP9erVq6KiYtmyZZ1tJT09nWXZnj17vvrqqzk5OfPnzx8zZszWrVuvXr06ZswYkUikUCjefPPNCxcuPP/884MHD54/f35ra+tfaf1G0/SYMWP4LL+ampqAgIDevXs///zzBw4cyMzMTEtL4zhuzZo1arV62LBhw4cPnzRp0vnz5zsbiKVUKh966CH+uD0rK2vBggVDhgzx9vauqKh45ZVX0tLSEELr1q3bunXrE088kZqaSpJkQUFBZ/sSExPTrVs3XgmpuLh4woQJCxYsOHLkSFVVlaenZ2xsrFwur6qqcjqdY8eOXbx48ZUrV27K4nUQAoEgNTWVP2u/evVqcnJyz54909LS9u7dq9VqZ82aNWzYMJIkP/zww8mTJycnJw8cODA/P79TZ/NtbW1xcXFPPPHEiBEjLBbL9u3b582bN23aNI1Gk5+fX1VVNWPGjKeeeurLL788efJkamrqnDlzeAnPjjeBMVar1VOnTh0xYsTQoUNzcnLsdvv48eNXrlx56NAhs9k8fvx4hmEIgrhy5UqXLl2Sk5NffPHF/Pz82xRU/DWsVqunp+f8+fOHDBkSHBy8bdu26dOn8y2eOnVKJpP961//GjZsGE3TfM2s8ePH+/n5dXC6uVfyf5RK5ciRI0UiEb/SaGpq2r59+7Vr1ziOMxgMK1as6NGjR1FRUVxcnEqlAgCxWMzH9dwmN/DXSEpKksvl+/fvB4AJEyasXbu2vLz8hx9+SE5OZll23759ERERPj4+ZrM5ICAAAFQqVWVlZffu3TuxoSSI4cOH79mzB2OMEEpISFi/fr1MJqNpmqKoysrKyZMnl5aWhoeHX7hwgRem8PX11el0/0OZ6L8RHBwcHBx8s1PvvfeeVCr9/vvvlUolQRArVqwQCARnz56Njo7u0aMH/7qaeeHVDkOhUOzevfvw4cObN2+eM2fOO++8w7fo7u5eUFDAC2IGBgbqdDqCIBQKBb+61mq1Xbp06WATCKFBgwYlJCRs2rQpJydn06ZNGRkZ/BeqUCiEQuGAAQN4+UKTyeTt3V5oEWNsMpk6vlAUCoXvv/9+fn7+zp07RSKRw+HgP6FYLK6trT1w4EBubu6WLVsGDhx4+PDhhIQEAPDz8+vgjHZv+Q/4YCKFQjFp0iSGYXr27Onm5kbT9EcffbR27dphw4bV1tY2NjYCAJ+2cRthlN/CzXCymJiYKVOmUBTl6+sbGBjIcdzevXvHjh0LADKZrL6+HgDUanXQHaWI8xFoCKFRo0Z17dq1W7duCKHa2to+ffocOnQoPT09IyODJwxBEPX19W53VO6BX6B7e3tPnTqVZdmYmBixWPz444+HhoZKJJKuXbt6e3vzBsFoNCoUis5uGxYuXHjkyJHXXntt9OjRPj4+dXV1vF2KiYnh309dXZ1CoVAqlXq9HgBKS0v52a3j0Gq1zz//PMdxR44ciYiI4GOf+VakUunNy3j+0zSNELLb7RKJpLM73k2bNg0fPvzdd98NCgqqqanhN0W+vr6rVq364osvZs+e/eSTT3bp0oXP/lCr1R3c995b+ae8XdZoNFu2bElNTf3www+fffbZlpaWRYsWjR8//ptvvlmxYkVaWlpoaOjBgwcfe+yxO+APx3E8906fPp2dnZ2SkiISiUJDQ9Vq9fXr1/nZevbs2ZMnT66treU4Li4u7g46wmuQ8zk/ISEhBQUF6enpERERK1euDA4O3rhx45IlSy5evLhmzZrGxsa4uLg74w8f81pVVfXBBx/MmTNn+/btr7322urVqwcNGqTVamUy2UsvvRQbG+vt7X3o0KFVq1Z1quBFSUnJmjVrZs+evXXr1ri4uKeeemru3LlFRUVSqTQuLs7Ly+u11177+uuv161bZ7fbV65ceeLECW9v78DAwE51Yd++fcePH5fL5UuWLBk0aJC3t/eqVavKysqSk5P5C0wmE8aY37qsXr36woULffr06ZSXoqWlZcOGDVFRUYcPH87Pzx88ePDkyZPFYnFmZuYLL7ywYMGCqVOnfvXVV+fPn09LSxs6dCjDMLwPqUNT2D0VP5qXl9evXz8AuHz58p49ewYPHjxixAje/3bp0qVp06Z17dq1rq7uk08+SUhI+LkUcMdhMpnKysp4/9uuXbvKy8ufeOIJfh1SUVHBm2+McU5OTk5OTnp6emcnVB7FxcUBAQFSqbS1tXXz5s0BAQHTp0/nOC43NzcrKys5OXnQoEEGg2Hr1q0Mw6TzequdR05OzuDBg/lfjhw5kpKS0q9fP4PB8OGHH8pksrS0NIFAUF5evn379uTk5P79+3fq4a2trZcuXXI6nRzHeXl5JSYmnjhxIicn5+mnn1apVK2trR988EFSUtLo0aNZlj1w4EBRUdG8efM663+rqakpLi5GCGGMw8PDvb29t2zZIhAInn76af6CU6dODRw4kP88GzduDA4OntHJippWq/XKlSs6nY63Y7169SovL//8889TUlIeeOCBH374gU/EUCgUvXv3vnz5cmZm5vTp06Oiov56/HHBhb8WXOenLrjg4o8LLrj444ILLv644ML9Asr1Cu59YIx/+OGH8PBwHx8fADh//jxCqLGxkdfKUCgU/fv3F4lE165dEwqFISEhAHD48GGDwUBRVI8ePTp+oOlCZ3FvxY+68FvIysrauXPngw8+WFtbO2XKFB8fn23btlkslpqaGrPZ3KNHD6FQOHLkyJqamuTkZIIgJkyYgDGura395z//OWDAAH9/f9c7dNmf+xQIoQkTJuzdu7esrGzu3LkvvfQSSZL9+/d/7bXXCILgT0WvXbvW1tbGS1iFhIS4ubm9//77YrFYIBBkZWXxYewuuOzPfQo+mOWtt94iCGL16tXnz5/Pysqqqqo6ceKEwWCIiopasGDBww8/rNfrJRJJZGTkqlWraJo+d+5cRkbG66+/7rI/Lvtzv2P69OmLFy/et28fAHAcFxoayheTValUBoNh3759FouloqKisbGRj8wQi8U+Pj7u7u6uI3IXf1wAhJBIJJLL5bxHQaVSxcfHI4QoilqzZs20adPef//95ubmkSNH6nQ6sVg8cuTIwMDAzMzM/Px81/rNxR8XICYmhg/v9/f3/+STTzIzMzHGERERFotl48aNUqlUKpVOmjTpyJEjcXFxvAhtUlJSU1OT69XdrUnNZdz/TuCTjn795y/+74KLPy648OfDFX/gggsu/rjggos/Lrjg4o8LLrj444ILLrj444ILLv644MK9iP8Dlmuf+QluggkAAAAASUVORK5CYII="
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise on Linear Regression\n",
    "\n",
    "This notebook will show you how to apply a linear regression to a time series of monthly averaged CO2 air concentration measured in Mauna Loa (data extracted from here: https://www.esrl.noaa.gov/gmd/ccgg/trends/mlo.html)\n",
    "\n",
    "![co2_data_mlo_s.png](attachment:co2_data_mlo_s.png)\n",
    "\n",
    "The linear regression will allow us to calculate the rate of increase of CO2 in the atmosphere, and how this trend has changed over the last decades."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using PyPlot\n",
    "using DelimitedFiles\n",
    "using Statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#This cell contains the script that calculates the linear regression.\n",
    "#As you can see it follows the equations seen in the slides in class\n",
    "#You need also to execute this cell (Shift + Enter as always)\n",
    "#Once you execute it you can call it later on in your code\n",
    "\n",
    "\"\"\"\n",
    "    a,b,r2,sigma_a,sigma_b,sigma_e = mylinreg(x,y)\n",
    "\n",
    "Linear regression\n",
    "\n",
    "y ≈ a + b *  x\n",
    "\n",
    "\"\"\"\n",
    "function mylinreg(x,y)\n",
    "\n",
    "    # remove NaNs\n",
    "    ind = .!isnan.(x) .& .!isnan.(y);\n",
    "    x = x[ind]\n",
    "    y = y[ind]\n",
    "\n",
    "    xm = mean(x); xa = x[:] .- xm;\n",
    "    ym = mean(y); ya = y[:] .- ym;\n",
    "\n",
    "    ssxy = xa'*ya;\n",
    "    ssxx = xa'*xa;\n",
    "    ssyy = ya'*ya;\n",
    "\n",
    "    b = ssxy/ssxx;\n",
    "    a = ym-b*xm;\n",
    "    r2 = ssxy^2/(ssxx*ssyy);\n",
    "\n",
    "    n = length(x);\n",
    "\n",
    "    sigma_e = sqrt((ssyy - b^2 * ssxx)/(n-2))\n",
    "    sigma_a = sigma_e * sqrt( sum(x.^2)/( n* ssxx ))\n",
    "    sigma_b =  sigma_e / sqrt(ssxx);\n",
    "\n",
    "\n",
    "    return a,b,r2,sigma_a,sigma_b,sigma_e\n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#We download the data if they're not already in your disk \n",
    "#(if it's already in your disk a message will appear over a red background)\n",
    "\n",
    "filename = \"co2_mm_mlo.txt\";\n",
    "if !isfile(filename)\n",
    "    @info(\"downloading $filename\")\n",
    "    cp(download(\"https://dox.uliege.be/index.php/s/lDbmzx8SQae6sdt/download\"),filename)\n",
    "else\n",
    "    @info(\"$filename is already downloaded, nothing to be done\")\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = readdlm(\"co2_mm_mlo.txt\", comments=true);\n",
    "@show size(data)\n",
    "\n",
    "#If you look at the text file (with a text editor), you will see it has a long header and then the data organised in 7 columns:\n",
    "#{Year} {Month} {Decimal Year} {Average CO2} {Interpolated CO2} {Trend (seasonal correlation removed)} {#days for daily average}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#We will creat two time series: one with the decimal year and another with the interpolated monthly CO2 averages:\n",
    "time = data[:,3];\n",
    "co2 = data[:,5];\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#We visualise the data with a simple plot, which gives the already famous curve\n",
    "plot(time,co2);\n",
    "ylabel(\"[CO2] (ppm)\");\n",
    "xlabel(\"time (years)\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now, we simply call the regression function we wrote above (mylinreg).\n",
    "#You provide variable x (time) and then the independent variable y (CO2)\n",
    "\n",
    "a,b,r2,sigma_a,sigma_b,sigma_e = mylinreg(time,co2);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Parameter a gives us the y-axis zero crossing\n",
    "#i.e. the time in which, if the trend would be constant in time, the line would cross 0 ppm of CO2\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#with our 2 regression parameters we can then vuild the regression line with equation y' =  a + bx\n",
    "yp= a .+ b*time;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#We re-do the plot above, adding a new line in red with the regression\n",
    "plot(time,co2,label =\"CO2\")\n",
    "plot(time,yp,\"r\",label = \"Linear fit\");\n",
    "ylabel(\"[CO2] (ppm)\");\n",
    "xlabel(\"time (years)\");\n",
    "legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Your interpretation of the results\n",
    "\n",
    "\n",
    "$\\rightarrow$ What is the trend provided by the regression? (in ppm per day) \n",
    "\n",
    "$\\rightarrow$ What do you think of the regression? Is it good/bad? \n",
    "\n",
    "$\\rightarrow$ What would you do differently?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----------------------------------------------------------------------------------------------------------------\n",
    "\n",
    "Before proceeding with the exercise, let's comment on something:\n",
    "When doing a linear regression with data that has a seasonal variability in it, as the CO2 concentration in Mauna Loa, the choice of teh initial and final points will have an effect on the trend.\n",
    "\n",
    "Let's do a kind of an extreme example. If we had only a (very) limited amount of CO2 data and wanted to know if \n",
    "there is a trend of increasing concentration:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(co2[708:715,1]);\n",
    "ylabel(\"[CO2] (ppm)\");\n",
    "xlabel(\"time (years)\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Applying the regression over such short series:\n",
    "\n",
    "aS,bS,r2S,sigma_aS,sigma_bS,sigma_eS = mylinreg(time[708:715],co2[708:715]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#We plot our reduced times series and the regression line\n",
    "\n",
    "plot(co2[708:715,1],label=\"CO2\");\n",
    "plot(aS .+ bS * time[708:715],\"r\",label=\"Linear fit\");\n",
    "ylabel(\"[CO2] (ppm)\");\n",
    "xlabel(\"time (years)\");\n",
    "legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see in the above plot, the trend of CO2 concentration over this short period is decreasing with time!!\n",
    "Just an example of how you need to be careful on how you handle your data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's return to the full time series. We saw that the linear fit might not be perfect for the data we have.\n",
    "An alternative answer would be to apply another type of fit, like exponential fit.\n",
    "But let's keep it simple :) \n",
    "\n",
    "### Exercise: \n",
    "Apply the linear regression separately to different periods of time (early period, mid-period late period) and \n",
    "look at how the trend has changed with time.\n",
    "\n",
    "$\\rightarrow$ What is the CO2 increasing trend in the 1960s?\n",
    "\n",
    "$\\rightarrow$ What is the CO2 increasing trend in the 1990s?\n",
    "\n",
    "$\\rightarrow$ What is the CO2 increasing trend in the 2010s?\n",
    "\n",
    "Calculate your regression and make a plot including all three lines. You can use different colour for each period.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#You can use these cells for your calculations\n",
    "#If you want additional cells, you can do so in the \"Insert\" menu above\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#When you have the regressions, add them to this plot, remember to give a label to each line\n",
    "#add also the trend ([CO2/day] for each period to the figure legend)\n",
    "#you can save the figure by changing YOURNAME by... your name :) in the command below\n",
    "#We'll visualise them all in a slack channel I've set up: https://join.slack.com/t/ocot/shared_invite/zt-cxh2uirp-dfQe6HTezFA42xuJ~lnBQg\n",
    "\n",
    "#You'll need to log-in/create an account to use that. When you have your figure, upload it there.\n",
    "\n",
    "\n",
    "plot(time,co2,label=\"CO2\");\n",
    "ylabel(\"[CO2] (ppm)\");\n",
    "xlabel(\"time (years)\");\n",
    "legend();\n",
    "savefig(\"regression_YOURNAME.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Julia 1.5.2",
   "language": "julia",
   "name": "julia-1.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}