{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Implementing linear regression in Python\n",
    "\n",
    "Like many data analysis problems, there are a number of different ways to perform a linear regression in Python. This notebook shows a few different methods. The final method motivates a review of matrix multiplication, which will be helpful in the next section on multivariate regression.\n",
    "\n",
    "## Example: 2007 West Coast Ocean Acidification Cruise "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from scipy import stats\n",
    "import statsmodels.formula.api as smf\n",
    "import pingouin as pg\n",
    "\n",
    "import PyCO2SYS as pyco2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load data\n",
    "\n",
    "##### Load 2007 data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename07 = 'data/wcoa_cruise_2007/32WC20070511.exc.csv'\n",
    "df07 = pd.read_csv(filename07,header=29,na_values=-999,parse_dates=[[6,7]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['DATE_TIME', 'EXPOCODE', 'SECT_ID', 'SAMPNO', 'LINE', 'STNNBR',\n",
       "       'CASTNO', 'LATITUDE', 'LONGITUDE', 'BOT_DEPTH', 'BTLNBR',\n",
       "       'BTLNBR_FLAG_W', 'CTDPRS', 'CTDTMP', 'CTDSAL', 'CTDSAL_FLAG_W',\n",
       "       'SALNTY', 'SALNTY_FLAG_W', 'CTDOXY', 'CTDOXY_FLAG_W', 'OXYGEN',\n",
       "       'OXYGEN_FLAG_W', 'SILCAT', 'SILCAT_FLAG_W', 'NITRAT', 'NITRAT_FLAG_W',\n",
       "       'NITRIT', 'NITRIT_FLAG_W', 'PHSPHT', 'PHSPHT_FLAG_W', 'TCARBN',\n",
       "       'TCARBN_FLAG_W', 'ALKALI', 'ALKALI_FLAG_W'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df07.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear regression: five methods in Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create $x$ and $y$ variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = df07['PHSPHT']\n",
    "y = df07['NITRAT']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'nitrate')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAGwCAYAAABcnuQpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAABy90lEQVR4nO3deXhU5dk/8O+ZSEISQkjYspCNCIQ1JIASQDZfUcAiy9uitErVKpZFqa+FoLWV1pJgq68boKBF+vZCaQ0gP1mKlTVGMUCACAY0JiTCRIiECSRjIpnz+2M4wyznzJzZl3w/10UrZ86ceXIYPTfPcz/3LYiiKIKIiIgoSGn8PQAiIiIidzCYISIioqDGYIaIiIiCGoMZIiIiCmoMZoiIiCioMZghIiKioMZghoiIiILaTf4egLcZDAacP38eMTExEATB38MhIiIiFURRxJUrV5CUlASNxv7cS8gHM+fPn0dKSoq/h0FEREQuqK2tRa9eveyeE/LBTExMDADjzejcubOfR0NERERqNDY2IiUlxfQctyfkgxlpaalz584MZoiIiIKMmhQRJgATERFRUGMwQ0REREGNwQwREREFNQYzREREFNQYzBAREVFQYzBDREREQY3BDBEREQU1BjNEREQU1BjMEBERUVBjMENERERBjcEMERERBTUGM0RERAFAq9OjpLIeWp3e30MJOiHfaJKIiCjQbSqtwbLN5TCIgEYACmYOxuwRqf4eVtDgzAwREZEfaXV6UyADAAYReHrzF5yhcQKDGSIiIj/R6vT48MR5UyAjaRNFVNc3u3Xd9rRkxWUmIiIiPzBfWrIWJghI7xbl9nXby5IVZ2aIiIh8zHppyVyYIGDFzEFIjI10+7rtZcmKMzNEREQ+VlXfJBvIPDu1P6YMSXQpkFG6rrRk5eo1gwFnZoiIiHwso1s0NILlsTBBcCuQAYDyb3U2x9xZsgoWDGaIiIh8LDE2EgUzByNMMEY07iwtSbQ6PVbuqrA5fn9eKt79vAZ//7QqZJebBFEUZSa6QkdjYyNiY2Oh0+nQuXNnfw+HiIjIRKvTo7q+GendolwOZLQ6Parqm/D91RYseveYw/NXzpJPCJauk9EtOiCWpJx5fjNnhoiIyE8SYyPdChysdy6pkV9UjrF9uyMxNtIUwJSf02Hlzoqg3QHFYIaIiCgIye1cUkMEsP2EFm2iiMIdFbB+m7QDKishBk2tbQEzU2MPgxkiIqIgdORsg+oAxtrz27+0+3qbKGL6qhKICI6ZGiYAExERBQDzqr2OKvhuKq3Bwo1lXh2PFCcFQ60azswQERF5kCuJtNbVgAXANCuydHIWBifHmq6n1emRX1TutfFLn20u0GvVMJghIiLyEFdaCchVAzafFSnYYdxuLV3veO1lm2DDk+SuHei1ahjMEBEReYBSKwFp55ASpWrA1gwiFFsgeJMGcLsGjrcxZ4aIiMgD7LUSsCejWzRU7qp2O5B5bFwG/jR9oFPvEdUOzo8YzBAREXmAUosCR8szibGRyJ+cpeoz1NaSUfLG/io8u/WkU+8RgyABmMEMERGRB7jTomDeuEwsm5JlN1iRkoH9MVGiZobJn5gzQ0RE5CGzR6RibN/uii0K7O10mjc2E9Oyk7D9hFa2Dsyr9+YgvlO4V5N/7YkKD9z5DwYzREREHqTUosDeTifzIGfqkESs2PGlRX5MmCAgJT4SJ2S6YvtKc6vBb5/tCIMZIiIiJzlbS8beTqcDZy4iv6gcIow1XgpnDcaMnGQUHT1nen92SixmrC7x+U4mCbdmExERhRBXasko7XQ6erbBFMgAxhovS4vKbfJijtZc9tDoHdMAmJGbjK1l59Emik7l/vgLgxkiIiKVXK0lI+10sl46+r6pRTYHxl95MRoB2DJ/FLJT4vDUnf0Uc38CTeBm8xAREQUYV2vJWO900gjAQ2PS0fjDNW8N1SUFMwcjOyUOgHHMeZldAz6QATgzQ0REpJrcDItGULfTR9rptP6TKqw7UIV1B6u8OFLnPT99YEB3xraHMzNERNSuOepQbc56hgUwBjYzVpdgU2mNqs9762CV35aR7BEDcVAqMZghIqJ2a1NpDUYV7MGcdYcwqmCPw4BEq9MjJT4Kax/ItUjSNaisknu4+pLfdiQ5IgRB2wIlXGYiIqJ2SavT2+wkyi8qV0zmNd/FJAi2SbptoojtJ7SYOiRR9v1v7q9Ewc4Kj/8cAPDzW1OR17srPvqyDh8c07p0jS6R4R4ele9wZoaIiNqlw9WXbAISEcCR6gabc613MSktyTy//UuMLrSd4XnzgPcCGQBYOPFmNLVewzYXAxkAGJYeZ/f147UNWHewEsdrbe+Pv3FmhoiI2iVBYV3F/LBUHO/7qy2ql4est2trdXoU7HAvkOnXsxNOf3dV8fWPv/wOv//gpFu5ONuOn8e8sZmyr/3PP49ZFPGblZuMF3821I1P8ywGM0RE1C4NS4uDAMvlIkEActOMMxTWxfGsz7VH2q6dGBuJqvomp8cmfZZGAO4alIBdX9TZPf9Co/pgS8nKnRUYmRGPptY2i8rGx2sbLAIZACg6eg4P5KWZtnH7G5eZiIioXUqMjUThrMGmTtUaASicOdg0m2JdHA/CjYdmmCBgVm6yxa4mayfOXQZg3M7tLBHAnFtSsXzaQOwsr7MbqAgwtjtwl0EEpq8qwZx1hyyWyj6vviR7/mGZ5Th/4cwMERG1W0pdruWK44ki8PqcHMRHR5jOlarkFn91Eav2VVqc/8LO05iWnYTE2Egsm5Ll9FLTxs/VbfXOn5yFyHDPPM6lH9l8qeyW9HjZc4c7yLHxJc7MEBFRuyZX6VYqjmcuTBCQmxZnca703tY2247S5pWB543NxLLJWV4Z/5BeXWTH64hw/RcA2fdK489OicOs3GSL12blJgfMEhPAmRkiImrn5DpgS8Xxnt78hcNmi1qdHm/JVPPVAEjvFoXjtQ34vPoSRvaOxwcLRuGeVSUeG7vUzToxNhIPj8lQVVVYIwCv3ptj2r1UXd+MqHCNTVdu807ZL/5sKB7IS8Ph6gYMT48LqEAGCKBgpqCgAE8//TSeeOIJvPzyywAAURSxfPlyrF27Fg0NDbj11luxatUqDBw40L+DJSKikGBROwbAI7dl4MExGUiMjVRcgrL2n1PfySYGTxuahCXvn8DBr+pNx27r082t8QowJikbRNgEWA+pCGaE612+785OMh2T3u8oeMtOCbwgRhIQwUxpaSnWrl2LIUOGWBx/4YUX8NJLL+Gdd95B37598fzzz+OOO+7A6dOnERMT46fREhFRKLCpHQNg7cEqvFVchYKZgzF7RCoSYyPtNlrcVFqDZz84Kfva1mPnbY6ZBzauEAG8dm8OunaKsAmwDpy5aLHjSoAxn6ZXXCQuNbUiPjocuWlxij+P2uAtEPk9mLl69Sp+/vOfY926dXj++edNx0VRxMsvv4xnnnkGM2fOBABs2LABPXv2xMaNGzFv3jx/DZmIiEKAXJIvYFsnBrBcigKAI2cb0NDcit9vlQ9k7BmZEYfPqlzfCVRR14in7rTMv5ECM+tt5tOGJjkVlDgK3gKV34OZBQsWYOrUqfiv//ovi2CmqqoKdXV1mDRpkulYREQExo0bh5KSEsVgpqWlBS0tLabfNzY2em/wREQUdKTAJDo8TLF2jHny7t+Kq/B2cZXHeio9MrY3Pqs64vL7V++rxM9HplkEHXKBmUE0VjOO72SZD2RNLmco2Pg1mHnvvfdw9OhRlJaW2rxWV2csENSzZ0+L4z179sTZs2cVr1lQUIDly5d7dqBERBQSrAvh3danGw4oLP0Uf3URc9ZVerTD9azcZNRfbXXrGgYRpkBLCkKk3UzWAc2id8tMxfekpTNz1vdD7pxg4Let2bW1tXjiiSfwj3/8Ax07dlQ8z7rctCiKiiWoAWDZsmXQ6XSmX7W1tR4bMxERBRatTo+SynqH3aqlc60L4RV/rZzDsnqf5wKZu4ck4IMFo/DUnf2wbHO5W9fSCMCJby9jdOEeU4G7A2cuYuldtlu/revGmN8nufuhpvN3IPJbMHPkyBFcuHABw4YNw0033YSbbroJ+/fvx6uvvoqbbrrJNCMjzdBILly4YDNbYy4iIgKdO3e2+EVERKFnU2mNxQPdurmjNaWlGDnOtC5wRBCMuSg9OndUzNNx5lpLJ2dh5a4KmyAkMtz+I9186QyQvx/W5wQLvy0z3X777Sgvt4xOH3zwQWRlZWHp0qXo3bs3EhIS8NFHHyEnJwcA0Nraiv3792PlypX+GDIREQUIpVkFpaTdxNhIxaUYOZ5cWhJFYN3BKqw7WIVlk7NUj0EiAPjjPQNNu5GUghBHy1fmdWMAyN4P63OChd+CmZiYGAwaNMjiWHR0NLp27Wo6vnjxYqxYsQJ9+vRBnz59sGLFCkRFRWHOnDn+GDIREQUIe7MKibGRirkg5rVUnA0qPKFgZwXuHNgT/z75ner3iADio8MxdciN2jByQcjtWT3w2sdfywZicnVjEmMjMSMn2aKJ5PQc53Y/BQq/72ayZ8mSJdDr9Zg/f76paN7u3btZY4aIqJ2zN6tgb9bGvJaKXNVbX3AmkJGIZmNUqk6cnRKHwlmDLYK4pZOzMCS5i2zdGK1Ojy1llt2wt5adx1N39gu6gEYQRdHHf4y+1djYiNjYWOh0OubPEBGFkE2lNTYP9NkjUlFSWY856w7ZnP/uIyORl9nV5hrmgY8/3ZoRj8+rLtnMrAgAXpuTg2HXC96Zby1vbjXYBCpanV5V4Tul+/T6fTkWFYL9xZnnd0DPzBARESlRqlirNGsTFa5BSWW9RT2V2SNS0a1TOB7e4HrdF0+Y0K871j94C7Q6PdYXV2PdwW9MQY0IYOHGMmgEYEZOMraUnbNYPrMO0NQWvlPKIXr8vTI0tV4Lqi3a7JpNRERBS67jtbQME3a9jEeYIGB6ThJmrC6x2PkkbevW/2jb8drX9p2+CK1Oj8TYSDw9tT/yp2TBugiJQQSKjp7z2FZq6T5ZBwLBuEWbMzNERBRyrHNjpq8qsai5kr+5HBCNsx4CPLsV29pf/nsw9p+5iA9P1CmeIwKm5GWtTo+VOytUjcc86dkVs0ekIjriJizcWObR6/oagxkiIgpJ0nLLn7efsgkMzLNFvZ0us6SoHOL1ZaFp2UmyDSgBoLn1R5RU1uNSU6vqHB41W6m1Oj2OnG2AKIoYnh5vE6AMS4sL+i3aDGaIiChkaXV6vF1c5dcxiGbLQtuOn8f0oYnYekxrc96vNhxxaqZIA2DJXfZ3Hm0qrUF+UblFJ+3CWZYtC0JhizZzZoiIKGS5W3HXHXKNdwwi8MExrexrotX/a5Q790AAYACwcleFYuVjuU7aIoBlReU2bQ3ktmgHU84MgxkiIgpJWp0el5rca+rojnuGJskGJCIACDcewErn2AvC7PVckigFcgYg5NoacJmJiIgCmnldlabWNout1Urnlp/TYeVOY/8ibyb32vOBQm4MYFx6en1ODuKjI9wu3qeUrKu09VoDhFxbAwYzREQUsOSK2pm3J3B0LnBjx5KjWMHTQY+9a4UJAnKvF8EDgKV3ZRmDLxc+RynwkLZe528uN+XtCAAKZg2W3cpuXYAwmHJmGMwQEVFAsm5LIFFqKmmvkq8IYEZOEraUyc+W+HL2RiPAIljYVFpj7ILtwrUEq2tZk7aoHz3bAFEEUuIj0dTaZqppY32emsrBgYjBDBERBSR7ybvWSyuOEn3DBAGjMrsqBjPuBjIaADNyk7G17LyxieX1a8pdVxSBy3rjNuzo8DC32ikIIjC2b3eH58VFh6P8Wx0ef6/MpvmmRG3l4EDEYIaIiAKSUs4HYLu0Yu9cABie3gVL3i/30kiB/5nUFzOH9cJTd/YzzW4cOHPRtHRjTgRQsKMCgPszQlIyr1IQorT0Jje7Fcy4m4mIiAKSdVsCiVxOh1QrRcmhqgavLiP9ZfcZjC7cgwNnLpraK8wekYri/An43dT+iu9TO6YJWd1ldz3ZS9R1tPQWbDuW7OHMDBERBSzrtgRSl2gAFk0jtTo9Nh895+BqzhMA5E/OQlPLNby652u758rNdiTGRmLqkESs2PGlW/VuZuX0wooZg7H+kyq8daAKBsgHdebULL0F044lexjMEBFRQJC2VVtvvbbO5bBeOplzSwr6J3X26MyLAOCP9wzEfw3oicTYSHx4QnmbtTm5bdLWu4Ws2cuvMY1HMF7n6SkD8ODoDFWJuo6W6YJtx5I9DGaIiMjvzAMUpa3XAHC8tsGiPD8AbPy81uPjyZ+Shfvz0k2dtVPiImXzWwTBss+T0myH+QzTiW8v44Vdp02JwksnZ+HMd1cs2glYf0ZuWpzp92oTdeW2XC+Z3A9DkrsE5Y4lexjMEBGRX1nndiglp24qrTHWTPHBmIYkd7EJsGbmJmNL2TnLmQ7xRkDjaLZDCkLyMrsCAlB4vajfyl0VistBUmBnfk2lGSw5wb7lWi0GM0RE5Ff2yulb15GRWaXxOA2AqHCNTYC1tew81j0wzNQQEjDO1GiuV/M1L4Jnj1anx8qdFRYNKOU8O7U/pgxJVFxiszeDZS6Yt1yrxd1MRETkV1Juhznr5RpfNowc1687mlrbZAOsqvpmm5khA4D46AjVAYOan0UDYHh6HKrqm0x9l5RmsIKpIaS3cGaGiIj8Sk05fUd1ZOREh4ehf2IMDp+97NR49p6+iAGJnWX7FY1Ij7Pbx0jNEpCan2VISqypX5M0A5MSH+VwBqu94swMERH5nVST5d1HRqI4f4Jp6URKwAUgW3PGnqbWNqcDGcmqfZVYOjnL9HlSgJWdEmcxDg2Ah8ekAzAuAY0u3IM56w5hdOEebCqtkb22Uv0cc8dqdTYzMNHhYQ5nsNorQRR9sQLpP42NjYiNjYVOp0Pnzp39PRwiIlJJLj8kKyEG01eXuJ07M31oIu4YkIBecZGYvqpENqn49ftyMCw9TjZ5VqvTY31xNd4q/sY0PutZkzBBQHH+BLsdvo+ebcCCjWWqxvzuIyNRc6nJZgbLUc6MM5xJLnbmXFc48/zmMhMREQUcpfyQV+4b6pEk4B4xHU0Ju/fdmoqNh2xnUaTaLkoPaimQkcZnTVoCutD4Az6vvoRb0uORnWK5xXrqkEhcbblms91cTlS4xqu7k5xJLnYlEdmbuMxEREQBR2mHU/m3lz1y/bUHq0xLQYsm3ix7Tq845UBBbRLv/31WjXtWleDP2ytwz6oS/M8/j1mco9XpkRIfhbfmDoOjFbTmVmNfbWl7tycDGWeSiwMxEZnBDBERBQwpR0YuPwQA3jxQ5bHPkh7CALBy1mCbz5uxukQ270Wr0+P7qy2y4zM3bWgidpTXWRwrOnoOx2sbAFjm2Dzy9yOYmZOsmEfj7dwYe9vj3TnXV7jMREREAcF66WJoShccrbns1c+UHsKzR6Ta5OPIFe9T6kJtTQCQEi8ffByubkCPzh1l69hsnp+H5laDRZVgX7QekNthpRRAOXOurzCYISIiv5NbuvB2ICOJCjcuUjS1ttnk47SJIo5UNyC+UxOiw8MUAxm5VgdKuT1S/Ri52Y3mVgPyMrsiL7Mrpg1N8lnlXjXb410511cYzBARkd95oyjeT3OT8X7ZOYcJw7WX9GhqbUN0eJhsUPL4e2UwiLZ9mMzZBDIA1uz7BpMHJWDnFzeWmmblJiM7JQ5and7h7Ibayr2e2lXkTHJxoLVJYDBDRER+50pRPHsm9OuOv/xsKIZnxCt2qwaMMyqmYAXynaulMTm7i6pNFPFAXjoeG9cbh6sbMDw9zrSbyVOzG57eVeRM64NAapPAOjNEROQWT80MbCqtMT3clQILRwQA9wxNwtLJWRZ9narrmxEVrsH2E3U3asNc/wxvPQQd1ZkxH5srsxtanR6jC/fYzO44+sxgwTozRETkNFeCEndnBsw/c/aIVFzW/4hCsyaMzhIBbD12HtuOnzeNRfpZquqb8OCYdDw4Jh3V9c34vqkFC1UWrJOECQJu7R2HkspLDs9TM9PiyuyGdM8uNbWyvcF1DGaIiMiloESp3oj57h977/1bcRXeLq4yfeavx2Vi1b5KVeOVNjArxTwGEVhWVI7oiJtw7rKxS7X1z6bV6e3OAGkA9IqPRM2lG/VTuseEOwxkAGDJXf28UkTO/M9JgO3SmL93FfkLgxkionbO1aDEXr0R6/dJswnR4WHYfkKLtQct68UYRKgOZDQCsPSuLJz9vgnvfl6rHNAANjMv1j9b/uQsFOysUHy/eSADAHWNLarG+MKu05g2NMmjMyTHaxuQv7ncNGslwhjMSLlGgbCryF8YzBARtXPOBCXm1NYbUVubRQ2NAMwfn2lcinLxGubbracNTQIEmGZuPMXTyz2bSmtkWx6IAF67NwddO0X4ZVeRt/szqcVghoionXO1CJrSjhwAKKmsR0a3aADwWCADXJ/B2VvpdtKutINJWnb6JH8itp/Q4vntX7p0PW8u90gzZ3I/c5ggYFh6nF8CiUDqz8RghoionXNnm7B1vZEDZy6adthoBODhMRkerx/j6HJqdkKZL6nlF5Vj4cRMDE3p4vL28EfGZuDtg9VeKSKnVINHI8AmePR2UGO+XOhqvpQ3MJghIiK3i6CJEHGh8QebB9zbxa73UhqVGa8q2dacAGBgUmd8cb5R9XtEAK/tMebr5KZ2wfFand26NNavhAkCHhydgQdHZ3iliJzczJkGwJb5o1BRd8UiePTm7Ih18rH1ffDnTioGM0REBMC1bcKOHnDuzMqMvrmb08GMCDgVyFg7WnMZb88dhuO1Ory652ub11+fk4OrLdcUZ7G88SBXmjnr0bkjZqwu8cnsiHWSuNKSl792UjGYISIil6h5wNlrAeDIX/99Brmpnmk2KQVaYYKA6TlJ2FJ2TjHQ2vh5LfZ8ecHmeJggIDfNmJ/i61L+cjNnJZX1PqszY2+pKxB2UjGYISIil6h5wP16fG+8vlfdlmtrIoDjtTr85b8H47fvl8ueozbHZdKAHpg2NBm5aXH4679P233PxzKBjAawmYHx9YPb+jN92b1a6bOkLt/+7s+k8dsnExFRUJMecNZ+PS4Tr9+Xg1fuG4qfj0xDbmoXh9fSCMCdA3raHG8TRdRa1Xoxt2X+KKyakwOZYVj496kLWPRuGf7x6VkUHT3ncDzWXpuT47edOkqk5acwwfjTe3N2ROmzslPikJfZ1e+1bTgzQ0RELkvrGoWq+maLY6v2VZqWdRwFGRKDCPz71Hc2x8MEAd07Ryi+r7nVgKlDkizyWJRma+wV5uuf2An33ZKK339wyuY1jQDkpsWp/El8W3vFl92rA61TtjkGM0RE5JBWp8eRsw0QRRHD0+Nx4MxF2SJuEtHq/1318Jh0DEmOlU0uFgBEhRsXGKQH7ZHqBjQ0t+L3H5x06rMrtFeR3Ut+a7Z500pH/FF7xZdLXoHUKdscu2YTEZFdctVnXe1q7QqNAMzISZZdHjIPGBztrHLk3UdGouZS040ZHhgDmXnjMlW9P9S7WPsau2YTEZGJO8seStVnffm3YIMIbC07jw8WjEL5OR2e3Xpj1kXajpyVEGOzs0ojAHm9u+KTyu8dfoaUOJuX2dU0wwMBGObE8pKrbSHIfQxmiIhCmLvLHko7lnxNSgTW/9gmW6yttLrBZpwG0RiMOApmrBNnD5y56NI98+XuIrLEYIaIKES52g1beq9Utt7VEv+eJOBGPyVrYYKAEelxsoHE7f174LU9X9sEQHNuScXsEb1sthW7c8/caQtB7mEwQ0QUolxd9rCezZmRk4zNR8+5vbQkAHh6Shb+vKPC6fcB8gGVBsCSu/qhR+eOuHdECt79vNZUHE/aOpw/JQsFVp+5qbQWi26/2eY+uLtUFMg7fkIZgxkiohDlyrKH3MzEljL3AxnAmMeiJpAxT969e3AC+id2xl92n5E9zwCgcGcFCnZWWBxfMrkfxvbtjpLKeiR3sQ0o2kQR209oMXVIoscL0QXqjp9Qxt1MREQhbFNpjc2yh738j5LKesxZd8iHI7T19txhOPRNA94q/sbl5S1pNke0+mdrcjkxzt4z8g7uZiIiIgDOL3vIzUz42g8/GtwKZADLwEUKaOR+LrmcGC4VBR+2MyAiCnGJsZGqS84nxkZi6eQsH4xK2fdNLR4PpkQAr96bg99N7W/zmpQTY86Ze0b+x5kZIqIg4IkS+XLXMD8GGBNgk2I7emzcrugaHeFWUT65ZaUwQcCwdGPNmBU7vuT26RDDYIaIKMB5okS+3DUu639E4c4KiC5WzPWGMEFAblocHrktA2sPVql+35xbU9A/oTPqr7ZgYlYPVNRdUdwize3ToYcJwEREAcwTJfLlrhEowcut6XE4VN1g+v2s3GS8+LOh0Or0GFWwR/UYBQDC9ZwYKVizl/ei1emZExPgnHl+M2eGiCiA2at74s41AiGQ0QjA52aBDGBsW6DV6QE4N0YRsCl0B0Ax74U5MaGFwQwRUQCTdheZczbHI6NbNATHp/nc6Ju7ybYmqK5vRlV9k1vXdjbgs6bV6VFSWQ+tTm/xzxSYmDNDRBTAPFUiPxBmYqwVf1Vvc8w8ULO3FCYIgCAai+bJnedMwGedGG3dfRu40bjSlXwl8j4GM0REAc7duifuznI4SyMYt0Gf0+nxws7TaBNFaGAMPMzJBSqzR/QCYAziCmcNRn6RbcduKaCzTmCWrulMwGedGL30riys3FVh0X1b4kyfJvItJgATEQURV7ZoH69twD2rSjzy+WoShxdMyMQvRqaZGlU2txoQFa7BjNUlqurHmM+AaHV6HKlugCAAveIiTY0hAdgkNUtB1LD0OFX3Ri4xWi7osvbuIyORl9nV8Q9CbmEFYCKiEGRvi7a9IKeptc0jnz89OwnhHTT45+Fv7Z539vsmU5BgPk7z5TKNAIiifGBkPQNyd7ZtYFJSWS9bzbdrpwjVQZ5cYrQBxiUspb/msyZNYGIwQ0QUBOQaQEoP/ANnLtqtQ5PRLdruA1qtrcfPqzrvwxN1pn82H+fsEanISohBaXUDRqTH4cMTWqxTqCVj3qlaLlDzRENIpWssmdzPtDwmXJ+Kcnb5inyLwQwRURBQ2qJ9pLpBMcgBgP+c+g4XrvyAO/r3wO5TF3w8astxNrVetMlPsScqXKM4G5UYG4kZOckoOnrOdP70nCSnAg2l5OrZI1IxLTvJlKMEwCM1aTxRxZnkMZghIgoCSrMIkGme2CaKWP9JFdYeUF9B15Pk8moWvVdmmuEAjGMu3FlhNwfn2wbl2SgA2FJ2zuL8rWXn8dSd/ZwKFKxni7JTjC0PEmMjLa7jbvDhiSrOpMyvdWbWrFmDIUOGoHPnzujcuTPy8vKwc+dO0+uiKOK5555DUlISIiMjMX78eJw8edKPIyYi8g9pFiFMMO7b0QB4eEw6IjtoIFgVkdEI8FsgAwDj+3W3OSaXHyNCOZAJEwQYRFGxYKAnigkCxiBjxuoSPL/9S8xYXYJNpTVOvV8NpSVC1q3xHL8GM7169UJhYSEOHz6Mw4cPY+LEibjnnntMAcsLL7yAl156Ca+//jpKS0uRkJCAO+64A1euXPHnsImI/GL2iFQU50/Ao7f1BgRg7cEqPLzhiEUuTJgg4N5bPP83/gUTMlU/MPaevujWZ0nLPcPT4xULBnqimKA7QYYzhfQ8FXiRMr8GMz/5yU8wZcoU9O3bF3379sWf//xndOrUCZ999hlEUcTLL7+MZ555BjNnzsSgQYOwYcMGNDc3Y+PGjYrXbGlpQWNjo8UvIqJQ8lbxN7JbnDUANs/PQ//EGI9+ngBg0oCeWD59oHvXEW48dDSAbFXiObcYAzYpL8Z8Nso8Adfea2q5GmRsKq3B6MI9mLPuEEYX7nE4m+OJwIvsC5icmba2NvzrX/9CU1MT8vLyUFVVhbq6OkyaNMl0TkREBMaNG4eSkhLMmzdP9joFBQVYvny5r4ZNROQxahJE1xdXKdZqMQBobjXA0+XDRMDtOjVyzR+3HTuPgp0VFudtKq3FottvNv3eXsFAudecSbJ1ZUeUvV1lSp/nqSrOpMzvwUx5eTny8vLwww8/oFOnTtiyZQsGDBiAkhLjvzg9e/a0OL9nz544e/as4vWWLVuGJ5980vT7xsZGpKSkeGfwREQeoiZBVKvTY63CVmZJ8dcXsWZfpTeH6pLl0waafh7pIT64V6zNedLOJ/PaMtbJuObMX3M2ydaVHVH2ZnPsvc/dKs5kn98bTfbr1w/Hjh3DZ599hl//+teYO3cuTp06ZXpdsMpsE0XR5pi5iIgIU0Kx9IuIKJCpzd04crZB5t2WVu2tVFVl19fio8NtjsktvwDA4++VOZ2I60r+i1anl90RZe897iwZsVO39/g9mAkPD8fNN9+M4cOHo6CgANnZ2XjllVeQkJAAAKirq7M4/8KFCzazNUREwexw9SW7uRtSsumlphavjcGbXbUFAchNi7M5Li2/WD+IXNnt40r+iyvv8USuDnme35eZrImiiJaWFmRkZCAhIQEfffQRcnJyAACtra3Yv38/Vq5c6edREhF5hrQ0Yk36275cB2dPCxMEPDauN1Y5uTylVCNGrtO00sN+9ohUREfchIUbyyyOq1m6MedK/ourVYStl4wAY3sFFsPzH78GM08//TQmT56MlJQUXLlyBe+99x727duHXbt2QRAELF68GCtWrECfPn3Qp08frFixAlFRUZgzZ44/h01E5BHWSyMSAcCKmYNwofEHi67R3lg9+u2dfTEztxc+PKGuVYE5pfGIAFbNyUF8dISqxNxhaXFutyZwJclW6T2A4+BEytVhMbzA4Ndg5rvvvsP9998PrVaL2NhYDBkyBLt27cIdd9wBAFiyZAn0ej3mz5+PhoYG3Hrrrdi9ezdiYjy77ZCIyB2ulqmXW+YAjMHA51WXkL+53CsBjLnc1Hgkxkaid7doj10zTBCQmxanOjHXU7t9XEmytX7PgTMXZZtkSsz/rAE4vbOJvEMQPb2HL8A400KciMhZ7vzNXKvTmx6c/iAA2LpgFJpa2/D91RYseveY09cIEwQ8Nr431uwzJh6b9zcC5H/GMEFAcf4Emwe+Vqd3GIh4s7+Ro7Fa/1k/PCZDtlHmqjk5iIsO57KTm5x5fgdczgwRUbBwpeaIOWlGQm6pyRdEANNXlUCEa/k40nLY7BGp+MXINNlARK4ujlI+jL0t2ID3+xs5Sgi2/rN+u9g2kBEEYOHGMot8IS47eZ/fdzMREQUrT5Spnz0iFVvmj/LqbiJ73MnHeWvuMIvaMdK2Y2n31fHaBtmZCwFwuvqtL/ob2dt2LfdnLf1eeovmeka0eTNN9mDyDQYzREQu8lSZ+uyUOBTOsmwiGehm5Sbj9v4JNsfNS/1PX13isZwfX/Q3srftWqkmDmCcjXn9vhy8cu9Qm5+XPZh8g8tMREQuciZxVS7Xw/zY7BGpyEqIQWl1AyJu0uDZD076+sdx6KlJffFjmwETs3qgR+eONjt+rGdPlDIyRcCpbdcAUP6tzuaYN/obKSUR21sSNIhA107GnVvu7soi1zCYISJyg5odNHK5HgAsjs3IScaWsnNey52xfsg6a84tKVg4sQ8A488zY3WJTe6K0u4sa650t165q8Lm+JLJ/bySYKuUuyMFnFKekUT6ediDyX8YzBARucle4qpcrkd+UTkg3Ji5MIiw6A/kDQYRGJoSi2O1tjMcaiy63RjI2Et6jg4Pk33v3JFp+L/PzsIAz3W3BoAhyV2c+hk8sRNKWhJUCljYg8k/GMwQEXmB9OD8/mqLzYNYNP2Pb7kayAgADpy5qDj7IuWFiAo/1F2DE/HYhEyXH/CuVuqVaHV6rC+uwrqDVR7ZZeQoYHG0K4s8j8EMEZEDzv6N3npZSansfyB4bGwG3jggv8VYmjkScWP2xVFgofSaqw946d4vvSsLL+w67XRukvmfhcQTxe0YsAQWBjNERHY4W9tEbhlGEABB9G1AozaAkgtkANvk3TZRxJHqBsR3CrcbWHgyZ8T63i+dnIUhyV1U5yYtnZyFlTsrZJeonO39JPFm0T5yHSsAExEpcKZ6raSksh5z1h2yOf789IF4dutJnwU004cmYesx5/stAde3hsskDEsBkqPAwrySLwCXHv7O3nu5881nl6w5+nOUwz5MvuXM8zsYyhkQEfmFK7VNlGrP3N6/JwpnDVasVeJpLgcyAlAwazCW3pVl85p5MbgXdp5WnCGRCuhJfY7mrDuE0YV7sKm0RvU4nL33cucrBTIawOkZI18U7SPXMZghIlLgSlE868JrGgAPj0nHhcYfkBIfheemDfDiiG1pBODxiTfjp7nJjs8FsGX+KMwekYrBvWLtnusoqHP34e/svbdX1E6iEYBHb+uNT5ZNdHpGxRdF+8h1zJkhIlLgat0QqR7Jq3u+xp4vL2DtwSqsvV7W39dtCwwi8Oqerx2eJ/1s2SlxAOR3EFmfby+os/fwd6ZvlXTvNYL9ujKO+lw9O7U/pgxJdDkJ+VJTq00eEgviBQ4GM0REdjhTN0RKDi3/VoeCnbZF3oDA3NX0xO03495bUi1+NutgQgBMtXHkgjrrxFh3t1MDxnt/uflHFF5P4l25swJdIjsozqrYK2rnaiAjtxtKMj0niUnAAYIJwEREHmDvoRfoBACFs+STWa2TeeWCOqXE2E2lNTazWs4s77iSgC2Nx53Ptff55lxJIib1nHl+c2aGiMhF0mxEdHhY0AYygGUdGUcF4OR2LilVBHa3Gq6rS1WeqsLrqD2Dq9u7yfMYzBAROSBXWySYZ2LkuPpgdhRwOFNczvo+K7VHiAp3vHfFE0Xt3M0bIt9hMENEZIfcEsrYvt1DKpABXH8wl5/zTDdrufucEi9/jeZWg9PjdIVN3tD1DGARrvWYIu9hMENEpEBpCeXle7ODKpCRZhc0AO69JQWjMrvh3GW9qvYA9mh1eqyUSXSeNLCnR2q4bJ6f53YSsaPPdVTQz3rJCpDPGyL/YjBDRKRAaQlFp//RPwNy0aIJN0P/owFvFX+DjZ/X4r3SWhTMHIzi/AluPZiVckr+/UUdtDq96msq3efmVoPN1vgld/VDVX0TANv8HWe8ub8ShTsrVDWedJQ3RP7HYIaISIFczoQA4NkPTvptTK4YkhKLR/5+xGbmozh/AvIyu7p83Yxu0bItAwyAU/k39rZx52V2Nc2MnDh32dRrSW2fLLmZlzcPVFpsnfdE40nyL1YAJiJSIFfNF1Aukx+IJg9KQGT4TV6pXpsYG4n8ybZtD6yXgrQ6PUoq6xWr/8rd5yV33SiQlxgbifRuURZNIx1VFN5UWiPbSkGr06NQZmmM1XyDG2dmiKjds5c7YZ4zceTsJfx19xk/jdIx88J2ggAsGJ+Jp+7Mglan90juidx9mjc2E7he0M4A28RYtc0ZZ49IxWW9WYG8XRXoEnWjQJ7SUtSR6gbcna1+u3hVfZNsMKoBuDMpiDGYIaJ2zdHDVnqA/+PTs9jxRZ0fR+qYCEAjAq/PyUFuWpzFzIYrbRnM2csxmTcuE9OGJtnk39gLKuTq1azcWWEKNKzPVdom/fh7ZWhqvWbxZ2Zvu7jSdZZOzuISUxBjMENE7ZZWp0d+UblFN+j8zeWmB2gg1pLJTe2C47U6Y78iGPNTzBkAxEdH2J1hcjbhV02OiVxdF2eK3qmpV1MwczCWFZVb/MxyY7GXgyPX82np5CzjDBMFLQYzRNRuHTnbYNMrSRSBo2cbkJuGgAtkAOB4rQ6b5+ehudWAqHANZqwuUb18ZB1wyC0bWR9zlGNiLyhypj+TmnNnj0hFdMRNWLixzO5YHM1EeapCMAUOBjNE1G4ptaYTRcel7P2lTRRRe0mP+E7h6NG5o8vLR3LLawBkC9e5mmPizPKW2nOHpcWpCpAcBSyeqBBMgYONJomo3dLq9BhVsMdidkYAULJsIk6d1+HhDUf8NTRFAozJvdYViZ2ZZThe24Dpq0ssghTh+v9bd5vePD/PZvYHAJZNzsK8ceqWZo7XNqC0ugEj0uOQnRJn91zzxpZKP4uzjSTVFMejwOPM85vBDBG1a5tKa0x5GBoABde7R7+4uwKv7an09/AsKAUcznRu3lRaY5En5MiqOTm42nLN5RwTtQnWzgYaaoIeNZ9PgcsnwczBgwfx5ptvorKyEu+//z6Sk5Pxf//3f8jIyMCYMWNcGrg3MJghaj889WDcVFqD/M3lAVVPRgDwx+kD8exW24J97z4yUlXxO61Oj9GFe5xaPhMAFM6ynf1Rc6/lPs88+PJ2oGHv8wFwtibAOfP8dilnpqioCPfffz9+/vOfo6ysDC0tLQCAK1euYMWKFdixY4crlyUicpk7D0bz/AlpO3EgBTIAAAEYkhzrVr0YV/KARNhWC1Z7r+3tUAJgdyeZJyh9/vpPqvDWwSrO1oQQlyoAP//883jjjTewbt06dOjQwXR81KhROHr0qMcGR0SkhlI9E6XqsPYEauKvKMLUq0iqlOtsvRhpx5A5ATeWr5SYByBK9/p4bYNNlV+5z5OCL3s7yTxF7vM1ArDuQJXF+JdtLnf4XXFUxZj8y6Vg5vTp0xg7dqzN8c6dO+Py5cvujomIyCmOZgCckdEt2uHD3VOmD01U/VlSEDB7RCqK8yfg3UdGojh/guyMgtKD17ptQJggoHDWYBTOsmwlYD2mMEFAVLgGJZX1OHK2QfZeT19VYtM6QO7zpODL3k4yT5CWwZbelWXx+Q+PybAJogwisP6TKsVrKbVGoMDh0jJTYmIivv76a6Snp1scLy4uRu/evT0xLiIi1ZypZyJRyvnYdvy86uRYd/xyVBqemzYIeZk3duYo0QiwmYERFUbpaAlIacuy+bEDZy5a7BaanpNk2tEkzeTYzKpc/3/rInZKnzc8Pd7mOgKAYen2dzupYX0Plk7OQq8ukTCIIlLjo/B2cZVNQPbWgSo8ODpDtjKx2irG5D8uBTPz5s3DE088gb/97W8QBAHnz5/Hp59+iqeeegq///3vPT1GIiK7nC3Xr1SaX6vTo3CHbYE4TxvfrzuemzYIgDG4yEqIwfRVJTYPdun35nGO+YNaAJBvtkVa7YNXrsaK+THzAMS6MJ94fWxS8CjXGkCuiJ3c5xXOGmyzk8zdAEHuHhTurABEmP68h6fF4fNqy+UspU7fzlQxJv9xKZhZsmQJdDodJkyYgB9++AFjx45FREQEnnrqKSxcuNDTYyQickhtVVel0vxZCTH4uOKCT2ZlDpy5iE2lNaYZk+0ntIozHdI/S2M0f1CLgPFnEYwNH9U8eNXu+JICkJLKeptrigBeuzcHXTtFOF2F2Jw3KvHK3QPzYNAgwiaQAdyrTEz+53IF4D//+c945plncOrUKRgMBgwYMACdOnXy5NiIiJziqKqr0sxLmyjaFJHzJvMAqrZBj7UHlfM1JG2iiNJq23wVwNixelp2ksMHrys7vpSuOSz9RiNLd5pYeroSr1IjSXs0sF3GM/erMRnG3U9wPumafMOlBOCHHnoIV65cQVRUFIYPH45bbrkFnTp1QlNTEx566CFPj5GIyG1anR7vfl6jOPPi663YUgBl3WdIiQBgRHqcbMKwQbyxRKKUcKtmx5dc4rC9a0rUJCX7ivV4NYLj3VqvzcmRHbOU+Lv2YBUgAI/e1tvvPx/Jc6loXlhYGLRaLXr06GFxvL6+HgkJCbh27ZrHBuguFs0jIl90vxYATB6cgJ1f1EEULav1agRjsOTMx9skxwpASf5EbDt23mKZDLCtAixXHbeksh5z1h2y+Ryp4J6aSr3B1JjRfLzbjp3Hyp0VNh3GAeUKyo4K/pH3ea1oXmNjI0RRhCiKuHLlCjp27Gh6ra2tDTt27LAJcIiIvEVtFVpHgYzc7hy1nprUFxndotErLhIzzJaqpCDmtXtzMCw9zmKHkAaQfbCak6vBUl3fbEz2FYxLSwZRfqZEbunG3hKUmsThYGvMKI13U2kNVu66EchMyOqOA6frHS6JMfE3uDgVzHTp0gWCIEAQBPTt29fmdUEQsHz5co8NjohIiTtVaM1NGWScTXHVsLR45GV2lU2UNYhAXeMPAOzvELJm3kxSYp77Mm9sJqZlJzmcKbEO9pRyW+TGHgoPbrlAdm/FRSyckInRN3dHVLgGTa1t0Or0TgV/FHicCmb27t0LURQxceJEFBUVIT4+3vRaeHg40tLSkJSU5PFBEhFJtDo9DldfsplJWCZTCl+r0+P7qy12E0J3uBHImD/clBJPn9/+JVbs+NIUbJknzSr1f3rktt7I7BFtN6nW3kyJVqfH34qrTPVUzIM9ud1DofrgVgpkV++rREzHDsYZG4Vg2Nnt/uRfLuXMnD17FikpKdBoXMof9inmzBCFDke5L3cPTsAzdw+waWLoSdKSVJgg4LFxvdElugNuSY9HdkocNpUqF8CzzrdQ6l6tAfDJsommpF1n81SUruso38N87NKDO9gTXbU6PUYV7JFdQrReWrSXO3P0bAMMoojh6fEMZnzI640m09LSAADNzc2oqalBa2urxetDhgxx5bJERIrU5L58WF6HHV/UYenkLFNOiacJ1/Ngdnyhxap9labjs3KT8eLPhmJs3+7YfkKL57d/afE+82UbUzNLq2tbF45zNk9F6brWny/HGzVf/C0xNhL5k7NsEqbVFPqTHDhz0RQcSh3Egz3IC0UuTa1cvHgRd999N2JiYjBw4EDk5ORY/CIi8jS1DSANIlC4wzuBjHT96u+bsKPccnmq6Og5HK9tQGJsJKYOSVRssKjV6fHhifOy41PaImzOXsNDe/dIzbJRYmwk8jK7hkQgI5k3LhPLpmSZ/jzCBAFLJ2cp/vmY0+r0FrNcIoydvtlsMvC4NDOzePFiNDQ04LPPPsOECROwZcsWfPfdd3j++efx4osvenqMREROFUPzdsmYF3efkT1+uLoB2SnGYnIzcpJRdPSc6bXpOUk4cOai4uxSmCAgN81+XyJHSc9K90iut1N7Ipcw3SWyg8N8mMPVl2QrMx+pbsDd2e3zXgYql4KZPXv24IMPPsCIESOg0WiQlpaGO+64A507d0ZBQQGmTp3q6XESUTtnnZDpT0qf3tBsXHLX6vTYUnbO4rUtZecsghtzapJL1W6fNr9HGgC/Gpsh20CxvbFeslOzrHa5+UfZa/3ny+9wdzY3uwQSl4KZpqYmUz2Z+Ph4XLx4EX379sXgwYNx9OhRjw6QiEgiPYBe+/hrbPy8xt/DsbFm3zf4+cg02eUepRmlZ6f2x5QhiQ6DDbV1T0Ix98VbHOUkdYnqIHv8g2PnsXRyFu9tAHEpZ6Zfv344ffo0AGDo0KF48803ce7cObzxxhtITEz06ACJiKy9G4CBDHAjuJCWexzRAKoCGQCy11TKgwnF3Bd/GJ4eL3tchLGAIQUOl4KZxYsXQ6vVAgD+8Ic/YNeuXUhNTcWrr76KFStWeHSARETmquqbfNLZ2h57gcqJc5dl+wPJ+dVY9cs/anokkWclxkZi2ZQsm+OhUIMn1LhUZ8Zac3MzKioqkJqaim7dunliXB7DOjNEoeF4bQM+r74EiMCfZTpf+8rjE29GdkosHt5wRPZ183olUp0YuYq/GgH4JH+i08FIsPVICgVv7q809XYKlRo8wcCrdWZ+/PFH9OvXDx9++CEGDBgAAIiKikJubq5royUiMiPXb+l//nlMMXnWEwQAT93ZF+cu6/HeoVoYIF+LRABw362pqKpvUryWeR6LeU6GmmqyanpNBVuPpFAwb1wmpg113D6C/MfpYKZDhw5oaWmBIKhYECYicoLc1uOshBivBjLWf9NeNLGP6aH113+ftvjsmbnJpgeZ0jZxpSWIsX2745X7hgIiMCw9zuaBqLbXFPkHg8jA5lLOzKJFi7By5Upcu3bN0+MhonZKaevx1jLvBDJ/umcg3n1kJIrzJ9j05MnL7AoANturt5adNzUlLJg5GNZ/pxMU6rlsKq3B6MI9WLixDI+/V4YDZy5avK70s7M4G5E6Lm3NPnToED7++GPs3r0bgwcPRnR0tMXrmzdv9sjgiKj9UNp67K1Z4LjocNOSgdzyjqOt0GP7drcpOCOIxhkYc2rqw6jddk1E8lwKZrp06YJZs2Z5eixE1I7JVa8VAFxt8ewMsBQaLdxYBo0AzMhJxpaycxbLO2P7dselplbZZoTSEpLcrioDYBOAqAlUQrVrNZGvuBTMrF+/3tPjIKJ2Tq7Crwjgn4e/tTjPOsBQSyMAy6cNxB+2nbSYJTHPiTGIQP7mckCEqbGgeZds8yUkpQAkKlyDksp60yyPmkDF+mfntmsi57gUzEycOBGbN29Gly5dLI43NjZi+vTp2LNnjyfGRkTtjFS99j+nvsOzH5yUPUcKMpwNaAwi0HLN4LC3k3mxChHGIOi1e3Nsknbl+i9lp8SatmCbJ/GqCVRYuZfIdS4FM/v27UNra6vN8R9++AEHDx50e1BE1H4lxkbih2ttHr9umCAgNtL5/+QZRKBrpwjZbdTWCcJHay5bvE/KjVEbqHDHDJFrnPo3+8SJE6Z/PnXqFOrq6ky/b2trw65du5CcnOy50RFRuyIl4vbuFm33PFeWmfonxmBJUbnD8+zlyZiTy4WxplRzhog8y6lgZujQoRAEAYIgYOLEiTavR0ZG4rXXXvPY4Iio/bCus5Kb2sVipsNdX5xvVHXeI7f1xtvFVQ5zV+RyYawxiZfIN5wKZqqqqiCKInr37o3PP/8c3bvf2IIYHh6OHj16ICwszOODJKLQJrd9ucyDgYxaYYKAB8ek48Ex6aqWhKxzYabnJGFr2Xkm8RL5mFPBTFpaGgDAYDB4ZTBE1D7JLdl4u5mkRgDmj8/Emn3fyAYfaoIQuVyYp+7sxyReIh9THcxs27YNkydPRocOHbBt2za7506bNk3VNQsKCrB582ZUVFQgMjISo0aNwsqVK9GvXz/TOaIoYvny5Vi7di0aGhpw6623YtWqVRg4cKDaoRNRgFOzZOMK823VcrMms0ek4ucj09wKPqxzYZgbQ+R7qrtmazQa1NXVoUePHtBolLsgCIKAtjZ1OxHuuusu3HvvvRgxYgSuXbuGZ555BuXl5Th16pSpqvDKlSvx5z//Ge+88w769u2L559/HgcOHMDp06cRExPj8DPYNZsoOLx5oBIFHu6G/ad7BuLmHjEWlX45a0IUHJx5fqsOZnzh4sWL6NGjB/bv34+xY8dCFEUkJSVh8eLFWLp0KQCgpaUFPXv2xMqVKzFv3jyH12QwQxTYpB1M319twaJ3j3n02qvm5GDqkCSPXpPkqen4TeQMZ57fLtWZAYCPP/4YH3/8MS5cuGCRQyMIAt5++22XrqnT6QAA8fHxAIwJx3V1dZg0aZLpnIiICIwbNw4lJSWywUxLSwtaWlpMv29sVLeDgYjcY/4wA6D4YDM/7x+fncWqvZUAjDksniQAMIiiqTEkeQ87fpO/uRTMLF++HH/84x8xfPhwJCYmeqQRnCiKePLJJzFmzBgMGjQIAEx1bHr27Glxbs+ePXH27FnZ6xQUFGD58uVuj4eI1DN/mEn/NZCq55o/2MzPs+bpfBkAWPTusaB4uAbzrIaaRppE3uZSMPPGG2/gnXfewf333++xgSxcuBAnTpxAcXGxzWvWwZJop5PusmXL8OSTT5p+39jYiJSUFI+Nk4gsWT/MzGMS8wcbAOQXlXt9l9LsEb3wz9JvTZ+j5uHqz2Ai2Gc12PGbAoFLwUxraytGjRrlsUEsWrQI27Ztw4EDB9CrVy/T8YSEBADGGZrExETT8QsXLtjM1kgiIiIQERHhsbERkX2OKuFKD7bvm1q8HsgAwG19umNTqWVzSnsPV38GE6Ewq8GO3xQIlLcl2fGrX/0KGzdudPvDRVHEwoULsXnzZuzZswcZGRkWr2dkZCAhIQEfffSR6Vhrayv279/v0WCKiFwnPcyUSJ2kj1Rfcun6zvxHasH4TAxLi5Mdz4lzl22OKQUTWp3epbE6y96sRrCQigeGXZ8tZ7FA8geXZmZ++OEHrF27Fv/5z38wZMgQdOjQweL1l156SdV1FixYgI0bN+KDDz5ATEyMKUcmNjYWkZGREAQBixcvxooVK9CnTx/06dMHK1asQFRUFObMmePK0InIw+S6R0uk+i7TV5W4PCvz6/GZGNOnO17cXYHDZy8rnjehX3f89q4sAMDSu7JQsNNym/cLO09jWnaSxUPW30skoTKrwY7f5G8uBTMnTpzA0KFDAQBffPGFxWvOJAOvWbMGADB+/HiL4+vXr8cvf/lLAMCSJUug1+sxf/58U9G83bt3q6oxQ0TeJ9c9WiMAr96bg5T4SNyzqsSt66/ZVwkIsBvILJiQid/emWX6/eBesTbnyAUpnggm3Mm3kWuJEKyzGiwWSP7kUjCzd+9ej3y4mhI3giDgueeew3PPPeeRzyQiz5Kb3TCIgCAAH57Qun19A2Davm3t2an9MWVIos1DVG2Q4m4w4Yl8G85qELkvoIrmeQOL5hF5l1anx+jCPbJJwFI7AW8QAJQsm6j48Hcm0HClMrDczx0mCCjOn8CAhMgDnHl+u5QATEQkkWY35JJuvfk3pXuG2s7I2Hy+aPn/ShJjI5GX2dWpICQUkneJQgWDGSJy2+wRqXjl3qFOv08AMH1ooqrzrG07rsWm0hrZ86VdSlKsIcLzu5TkdnEFY/IuUShgMENEHjE8Pd7plgQigK3H7OfVCADyp2TZBDT2tlEfrr7k9VkTbkkmChwu92YiIrL2qzEZWHewyqPLS4+MzcC8sZlI7hKJhRvLLF6T26G0qbQGS4vKba7jjVkTJu8SBQYGM0TktjcPVKJwZwXE67uYPBnNdNAYJ5ClYnj2dihpdXrZQEYDeG3WhFuSifyPy0xE5JY391eiYEeF6mRbZ63aV4kPT5wHAIfLOkVHv5W9xv/c2Teo+h0RkXM4M0NELjte22BTaddV9iZ0Fm4sM22v3jw/D6XVDRiRHofslDjTOZtKa/Div8/Ivr/lR4NHxkhEgYnBDBG5ZFNpDfI32y7pOEsA8Md7BmJIr1i7bQ8MorHrtnB9qcm8doz17iVrt/fv4fY4iShwcZmJiBzS6vQoqaw37RwyBQ8eWFISAfxh20l8eEIru2vJ+ly5ppD2OnfPyk22mMEhotDDmRkiskuukm5KfJRi8GDPrNxkPHVnP/zn1Hf4/QcnTTMpBhFYd7DKuA17chZ6xUXiUlMr/rDtpN3PkXYzRYeHyb7+l/8ejJ8OZ64MUajjzAwRmSjNwFjPhuhbr9mdQZEjCMAtGfFIjI1EZo9OsktCIoAXdp1Gbloc7s9Lt0j41cC2eJ4GQFS4BrUN8sXwosL59zWi9oD/phMRAPUzMG2iiIc3HLE4JgjAAyPTsOHTs4rXF6/nvHTrFI4BSbE226zNry/VjrGu43LgzEVTU0jA2IRyxuoSzB6RIvuZDc2tKKmsd6mjNREFD87MEJHiDEx0eJiqqr6iCHSOdPx3IxHAwxuO4K//Pq3Yz8m6dox536TZI1KxeX4eBLP3GURgU2mt7Oc9u/Uk5qw7hNGFexRbHxBR8GMwQ9SOSctKSuX/m1sNFks99ry2p1L15xYdPYeshBh8kj8Rj47NMP2HSE1LgKbWNpvEY4MIi+sIsNzqba/1AREFPy4zEbVT1stKcnVeir+6iNF9umHtA7n41d+PeLQg3uHqBjx8W288PWUAHhydobolgNTg0boS8IOjM0zXqb/6Axa9e8zifXKtD4goNDCYIQpx0tZl87wRuWUlubmXVfsqsWpfpd2Cdq4ann5ju7QzLQGkBo9S7oz1bE5ibCS0Or3D1gdEFDoYzBCFMLmk3tkjUmXrstgLVtwNZIb06owT3zaafu9u7RdHDR4dBTxEFFoYzBCFKKWk3rF9u8su1WgEYyKvmsBFADB1cCJuzYzHs1tP2j03TBDw5v3DcaHxBxyubsBwqzYErnI0m8OO1kTtB4MZohAgt5QkN/si5Y3kZXaVnbm4rP8RK3dWqCqI92G5Fju+0NosQQkCIIjGbdPmMyKJsZE+r8TLjtZE7QODGaIgp7SUpJQoK+WNWM9cbDt2HoU7KyDCOPMyeVACdn1RB7kWjea7hKTgRcSNz+eMiC25gJOIPINbs4mCmNJSklanN+WNSNuq5fJGpBou246fR8H1QAYwBia7TtZh3dxhDiv9mi9NSbudzGvDkDHgHF24hzVviLyEMzNEQczeUpJcBV3z4EKaKYgOD0PhzgqbaxtEYOOhGkwbmogPjmlVjUfEjbwcBjJG9nKXeI+IPIPBDFEQc7SUBMjnjZgvTQkCFOvHfFxx0eaYNFMjLSvZC6bIccBJRO7jMhNREEuMjcSMnGSLY9Nzkuw+JK1nCpwthCed/ujYDGyZP8qmJQHruViSAk5zvEdEnsVghiiAWHetVnN8S9k5i2Nby87bLdt/5GyDqt1K9ogA3j5YjR6dOzrMy2nv1OQuEZF7uMxEFCCUdiUpHQccL2FodXocrr4EQRAwLC0OB85cRH5RuUfGK30O67k4xntE5F0MZogCgFKSaFZCjN3kUaWcmahwDf68/RTWHawyHXfUkmBsn244+HW96mUnATAtlbCei2O8R0Tew2CGKAAozbCUVtsuCbWJIraf0GLqkETTEob5zM30nCRMX1ViE7g4ilEOfFWPP94zAJrryyG/c1DZd/6ETD6ciSggMGeGyI+kXJjo8DCbJFENgIxuUTbHAeD57V9a1CsRzZJ5Nx8953Ivpd9/cArPfnASTa1tWDlrsOJ/ICb0647f3pml+rpanR4fnjiP/3f8nN18HiIiVwii6OxehuDS2NiI2NhY6HQ6dO7c2d/DITKxzoWZkZOMrWXn0Wb2r6R0XClA0QCAzPZoT1g2OQsje8fbzPIIAEqWTVQ9K7OptAb5ReWmawgACmfdyPshIpLjzPObMzNEfiCXI7O17DzWPpALwWwmxnB9pkWJAd4JZABg5c4K1DboZZerquubVV1D+jnNryECWFZUzhkaIvIYBjNEfqCUI1NV32yTgCtCXSdrTzNc/3B3aqTI/ZzStdUGREREjjCYIfIDpUJqI9LjZHNkvCkvMw4LJmTaHA8TBAxLj3O6KJ85uZ8TMP6Hh0XjiMhTGMwQ+YFSIbXslDiL43I8HeuMzjQm8y6bnGX6D4I0HgBOF+UDbiQ2A8a6OOY/jgCgYNZg7oQiIo/h1mwiP1EqpCYdP1LdgMffK7NYptEAWPnfg7GkqNzpNgRKMrpFAwCmDU1CUpeO0AgCctPikBgbiZLKeodF+arqm5DRLdo0frkifyX5E3H0bANEERiWHsdAhog8isEMkZfJPfAlSoXUEmMjcXd2JJpar+HpzV+gTRQRJgiYnpOEpVaBTJgg4LHxvbF6b6VLuTUdO2gsG08CeOS2DDw4JsNuI0u5oGVs3+6yRf6K8ydg6pAkF0ZHROQYt2YTeZG9VgRqaXV6VNc3IypcgxmrS2xmarYsGIWm1jbMWXdI1fUcVQI2Xfv6eAFYBFQrZg7C2L7dMbpwj02Q8/K92Vj07jGba737yEjkZXZVNT4iIsC55zdnZoi8RKlFgdSKQC1p9ub/HT9ns+RjAPBtgx65aXEQBHUdsNX+7cV8VqU4f4LFcpjS8pNGEBRncoiIvIUJwEReYq8JpLM2ldbgifeOyb62cGMZDpy5iCmDElwYpX3m+TF5mV1NQVh0eBisc5TDrufasEM0EfkaZ2aIvMRevokzrGd4rIkwzqCsfSAX28vrXB+wDLnxSktn1nk7UtDCDtFE5GucmSHyEqXt1+YPd2kL8/HaBpRU1stueVYqPGdOKrhnPTuTldBJ1Vhn5Sbj02UT8ehtvU11YZTGax1YaQBsnp9nkQtkPZNDRORNnJkh8iJ7sxTmycESuSTh8m91qj7r+e1fQiMACyZk4iZBQLeYCCR16YiHNxxRfM/AxM6495YU/NeAnkiMjcTTU/vjwTHpirMqcoGVAUBzq0HVGImIvIG7mYj8QKvT2+wGkoQJAorzJ5jquIwq2OPUlmsBgHB9ecvRziXpdbU7reTGbT5eIiJPYaNJIj+Qloysl4rkjttbOjJPEv5bcZXTtWNE3MjTcfRe6XVp55Kjyr5qls6IiHyNy0xEMuwVupOjVE/G+vjDYzLw0PVidPZmTU6cu4z0blF4u7jKoz+XRO6zzXcu2cMEXyIKNAxmqF2yF6w4W+hOqZ5MVkKMzfF1B6vw1sEq5E/Osju+F3aeRnKXSIeJv84SALw+Jwe94iJtCvA5s9NKqXIxEZE/MJihdkcpWNHq9DhytgH5ReU2yy/2Ct0p1ZMprW6QDUZEACt3VthdAmoTRZypu6K6Wq9a+VOyTG0FCmYOtqnsywCFiIIRgxlqV5RmUS7rf8TKnRWywUebKOLo2QbERcvP5CjVkxmRHmdzXKJm78+re772WIdsDYClk7Mwb2ym6Zia5SJnl9uIiPyBwQwFNWcftkqzKIU7KxRbAQiCscqu0q6fxNhIzMhJRtHRc6Zj03OSkJ1irIZrr+CdIyLU91KS8/z0gcjsHqMYrNhbLvJEXykiIl/gbiYKWptKazC6cA/mrDuE0YV7sKm0xuF7pFkUcxoo9zTSXI8k7O360er02FJ2zuJ9W8vOQ6vTY/aIVHySPxFj+3RT/4NZcTWQEQTg9v49XSpepzSD5Wi3ExGRPzCYoaDk6sPWemuxBsCvx2faBjgC8Pp9OXjl3qGKu34kanowFX9d78yP5zYBQOHMwS4vDXmyrxQRkbcxmKGg5M7Ddmzf7ph9Sy8AxtyVNfsrMSMn2aJ2SsHMwbg7OwnD0+NtAh3rXT9ysz3m56hpR+BJggBsXTDKrSUhRz8TEVEgYTBDQcnVh620NLXxUK3pmEE0Lgttnp+Hdx8ZieL8CaZAQE2ROEfn6Fuvuf3zOiN/chayU+LcugaL4xFRMGECMAUl6WHrzNZie92n20QRza0G5GV2tXlNbtePeeIxAKTER2Hz/Dw0txoskm03ldYgf3O5Z35olaZlJ3nkOv4sjsddVETkDAYzFLScfdjaW+5xNKtjvuvHfJePNDlkvtNJCoik4Mmd7mcLxmdi9b5Kp5KA1VTxVcsfxfG4i4qInMVlJgpqibGRqnfryC1NAcZ/CdQuoVjP7ohQ3unkbq7MrNxkdI7s4FQgE+x5LdxFRUSu4MwMtRvWS1MaAfjVmN54cEy66tkHRwFKmyjiSHUD7s6ORHR4mNNj/Mt/D0aj/hqGpxtzXu5ZVWL3/FGZ8Tj0TUPIVPG1l9gdzD8XEXkXgxlqV9zNA5Gr9mvt8ffK0NR6DSnxzs2QTBmcgOS4KIzpE40DZy6qyrUpqbyE8f26Y8zN3TAiPc7txF9/U6qmHMyzTUTkfYIourOiH/gaGxsRGxsLnU6Hzp07+3s4FAI2ldaYZneUqvOGCQJ+Pb43Xt9bqfq60rXM83CcESr5Jeb3V5ptCvafiYic58zzm8EMBS1f7HhR+gytTo/q+mac+PYyCnZWyL5XEJQrC3tLmCCgOH9C0C/JSPfX17uoiChwOPP85jITBSW1O17UBjzW52l1eqwvrsK6g1WyPZmka81Z95niNf3x14RQyS/xxy4qIgpeDGYo6CjteBnbt7vFA9A64Hl4TAYeGpNh85C0Pm9GTjK2lJ2zyNswiMCyonJ06xSOyPCbkNEtGoerL7ncN8lbmF9CRO0RgxkKOmp2vMgFPOsOVuGtg1UonHVjhkXuPPPu1+YMAB7ecASAMegZ36+7x382d2gE9VvMiYhCiV/rzBw4cAA/+clPkJSUBEEQsHXrVovXRVHEc889h6SkJERGRmL8+PE4efKkfwZLAUNNKwOlLdQijDMsH544b1pacqUWjEEE9lRcdP6NXjK2Tzd8kj+RibJE1C75NZhpampCdnY2Xn/9ddnXX3jhBbz00kt4/fXXUVpaioSEBNxxxx24cuWKj0dK/qTV6VFSWW8qnKamb5BSgTzAOMOycGMZRhfuQfk5neJ57hBwY1eSL3zy9fc+/DQiosASMLuZBEHAli1bMH36dADGWZmkpCQsXrwYS5cuBQC0tLSgZ8+eWLlyJebNmyd7nZaWFrS0tJh+39jYiJSUFO5mClL2En0d7Xgxf6+SMEHAkrv64YVdp01bgafnJGFr2XlTYT1XZm6Utmx7gtK1331kpGxvKSKiYBQSu5mqqqpQV1eHSZMmmY5FRERg3LhxKCkpUQxmCgoKsHz5cl8Nk7zIUaKvox0vUoG89Z9U4a0DVTDInNMmihjSqwuK8yeYAiMAmJDVAxCBlPhITF9V4nRg4q1AJkwQsPaBXPzq70csdksx8ZeI2rOA7c1UV1cHAOjZs6fF8Z49e5pek7Ns2TLodDrTr9raWq+Ok7xHKdH36NkGi2UnexJjI/H0lAH4ZNlEvH5fjs2SkgZAc+uPOFx9Cd83tWDb8fMYXbgHCzeW4fH3yvDhCa3fdiw9OCoN04cmmf4llZbTbu+fgEIHy2xERO1JwM7MSATB8ukjiqLNMXMRERGIiIjw9rDIB+RK2wuCMd9FrvaLPYmxkbg7OxJNrddM1WUByx1K1gwi8NbBKg/8JNfHDudmbNaXnIVGAJZOzsKQXl0sltPcbctARBRKAnZmJiEhAQBsZmEuXLhgM1tDoSkxNhIzcpItjomicpdqNWaPSMXm+Xmqk3M9OSsjCMDr9+VgwfhM1e8xiMALu07LBizOdAwnIgplARvMZGRkICEhAR999JHpWGtrK/bv349Ro0b5cWTkSdY7lcyPHa9twJYy+ZovEqm+jDM2Ha51K0hxdZeSQQS6dorAL/LSZK+hdF1XfkYiovbEr8tMV69exddff236fVVVFY4dO4b4+HikpqZi8eLFWLFiBfr06YM+ffpgxYoViIqKwpw5c/w4anKXVN+l/FsdVu6qsNipBMCU9Kumt5Gzia9v7q/ExkPu5VE9MDINf//srNMBkTTWqvom2ffed2sKRmd2My2jWb+PiIjk+TWYOXz4MCZMmGD6/ZNPPgkAmDt3Lt555x0sWbIEer0e8+fPR0NDA2699Vbs3r0bMTEx/hoyuUlpu7TULgBmOTJygYwgAIJozHVxNvFVq9OjUKEppDM2fHbW6fcIsKzOK7fl+71DtVg0sQ8KZw226Rpt3eTS2w02iYiCScDUmfEWds0OHFqdHqML9zhdt0UDy+BFbeLr8doGfF59CbekxyM7JQ4llfWYs+6QWz+DEqk+zeaj5xRnbH6Zl4ZJgxKQ0S0af7vexNKaVCtGqYaO2gabRETBLiTqzFDokGYSvr/aYjeQ0QAWMzOAMUjYPD8Pza0Giwe7oxmJB9d/jr2nb7QbmDI4Ac/ePcDlInj2PD7xZtx3ayq2HTtvd+npnU/P4p1PjTuU7hyYYPO6+XKSXA0dtQ023cWZHyIKNgxmyKusZxKUtidLsy4AbJZYslPinPrMB9+xDGQAYEd5HTK6RqNg5mDkF5V7dJdS354x2Hb8PApULmEZRGDnF7a1kpZM7mc3eFDTYNNdnPkhomDEYIa8Rm4mQRAAjVnOy5LJ/TAk2bKGijv1U47XNmCvQgPI1fsqsXWB53fCVX1/FS/u/srt6wxJ7mL3dbm6O55MDvbVzA8RkacxmCGvkZtJEEXg9Tk5iI+OUAxWHLUpsOfz6kuKr4kA/nXEtW3Zf/nvwTjz3VXZPBdPBDJqghKpwaa95GB3+GLmh4jIGxjMkNcozSTkpsV57eHYu1u03df/8ZnjbdmTByXg31/UWfRyWlJU7lSxO2c4E5R4s/Kvt2d+iIi8JWCL5lHwk2YSvNlDyLzonlanxzf1TW5fc9cXdTazN6IIvL630u1rm9MAWDUnB8X5E5zKS/FW5V9f/HkREXkDt2aT1yltM3aXUs2aYCAFCoGYXOutPy8iImdwazb5nL3tvO7kwNj7vGALZAQByJ+cZZPwHGi88edFRORNDGbIbf7YziuXrOqIlA8SJgh4bFxvnDyvw74z9d4ZoJknJt6MvgkxXs0VIiJqzxjMkFv8tZ1XLlnVkfnjMzH65u44ce4yVu6scDkYcta9t6YyiCEi8iImAJNb7G3n9aYDZy46bEJpbc2+bxAVrnE6kLk13Vi0z9UlrQuNP7j2RiIiUoXBDLlFmiEx5+3tvNJskLOxRZso4l+Hv3U6KDlU3eDkJ1k6bPZ+891XRETkGQxmyC3+2M7rSr6M5B+Hajw7GBWGX5/Z2VRag9GFezBn3SGMLtyDTaW+HwsRUShizgy5zZuF3OS4ki/jbRoBmJGTjKKj5yyOz8pNRnZKHFsFEBF5EYMZkuVs52Rfbue1LusvAJiVk4yisnMebSDpjFfvzcHd2Ul46s5++PjL73ChsQW39+9hapLJVgFERN7DYIYsaHV6rC+uwrqDVRARuJ2TZ49Ixf4zF7Gj3Fit9/2yc8hN7YLjtTpT36LpOUnYWnYebV6uC6kBMOz6UlJibCR+MTLd5hxvtApwNuAkIgpVDGbIRK6ibqAuhxyvbcCO8jqLY0drLuPtucMQFd7BtNz11J398NrHX2Hj5457Mrlq6eQsh/fG000i/VHbh4goUDGYIQD2K+oG4nKIUnfsE9/q8Js7+lkc82YgM3VwAuaNu9GA0t5siadyi5h/Q0RkicEMAbC/QygQOyffkh4ve/yVj79GUpdI0yxFlQcaT9rTeu1Gb201syWeyC1i/g0RkSVuzSYA8vViAOMXJFA6J5vXaMlOicOs3GTZ85ZtLjfVccnoFu3VMU0a2NM0NrnZEm/Uk/FHbR8iokDGYIYA2NaL0QjAo7f1xifLJno1F0NtETm5Gi0v/mwo/njPAJtzDSLw2p6vUFJp7LuUEtfRK2NPjY/ET4ffmAHyVSVkf9T2ISIKZIIoenmrh58500KcjMGFr+rFmC/LCDB2lJbLP4kOD8OM1SU2O4GK8ycAAEYX7lFcIhMEON32wJH+iTF4aHS6KZCRxmo9DmmM3rqPvvyzIiLyNWee38yZIQu+qhdjvSwjAijYWYHGH37Eb+/Msgl0rOMRadYjL7MrHh6TgXUHq2Q/x1OBjEYAfjWmNx4cky57fzy9W0kNX9b2ISIKZAxmyC+UEo5X7a2EAGD1vkqLQEfOtw1NKKkUcfeQRMVgxh0aAEunZGFIchdVsx++roRMRERGDGbILzK6RcvOuADGgEbNhMpv3y8HYFxK8qQnbr8ZI3t3cykg4WwJEZHvMQGY/CIxNhL5k7NkX3N2ZcjTOTETs3ogL7MrgxIioiDBYIb8Zt64TEwfmiT72pxbU27srPLhmKTGkObU7rhyxPw6nromERFxmYn8rH9iZ2w9dt7m+KKJfbBoYh9T/smBMxcVKxR7ws09ovHiT7NtAhlPtQ2wTmgGENC9r4iIgglnZsjrlGYhtDo9CndW2JwvPewTYyORl9kVAJASH4Ut80dhbJ9uHh/fsNQu+M+T42VnZDxRCE9u55YUk3mzuB4RUXvBmRnyKnszG1X1TbL5MSKAo2cbMHVIpGzzS0+akNUdK2YMxocnzqOhuRVdIjtgeHo8EmMjPdY2wF6rCFevSURENzCYIa/Q6vQ4crYB+UXlFrMQy4rKTQ0R7bUaWLixDN9e1qNwR4XTCcFqJHaOwBv3D0NF3RWMKthj8RkCgMJZgzG2b3doBNgUwnO2bYDUfiCYel8REQUTLjORx0mtBxZuLLMJRAwA1hdXm36vtKtaBFDgpUBm7sg0bF4wGrUNeotgy/yzlxUZt317om2AdfsBQbjxc7MVARGR+zgzQx5lnR8i563ib/DgmHTFZSZv+mDBKFTUXbHbAgEwBl3V9c0eK4RnfR1cvz6L6xERuY/BDHmUo/wQwLjc8sKuL3H3kCS7yy/O0gjAq/fm4NxlY2KxaPVawczB6NG5o02fJzkCgO+bWqDV6T1WCM/6OgxiiIg8g8EMqSI1fczoFm33IWyvsq+5LWVabCnTIje1C47X6tDmYuU76bOk5Zq7s411a6YNTUJ1fTOaW39EVX0zRqTHITslDiWV9aqDp4Uby7h1mogoCDCYIYc8VWtFztGayxjVOx4l31xS/R7zpo+A/HJNYmykRW0aadxySb0aAfjjPQMBGKsJ/2HbSZvt2FLSMhERBR4mAJNdztZacSUPxplA5hcjU7Fl/ig8PbW/adlGrvWA0rgB26TegpmD8YuR6fjFyHRk9uikuB2biIgCE2dmyC41tVbMl6DKz+m8Op5/fFaDjYdqHM4O2Ru3vaReuW3U3DpNRBTYGMyQXY4e7nJl+r1NzdJPdHiYTe6O+biVknqlbdRPb/4CbaLIrdNEREGAwQzZZe/hrtVZ1mnx5TZre1VzpQDLOpBRG5R4ajs2ERH5BoMZcsj84R4VrkFTa5upwq+v68RIpFkW611WcnVuNAKweX6eTe8lezy1HZuIiLyPwQypkhgbiW3Hzpvqt2gEYPaIFL+MRZplkdutlBIfZZMrYxCB5laDX8ZKRETex2CmHVNbOwYA3jxQiQKzDtcGEdhUWmv3Pb+9sy9afjTg1T1fuzxGKV8nTBCwZHI/DEnuYsp7Ma/iK+XRbJ6fxwReIqJ2hsFMO+Wodox5oAMAhWaBjMQgAlMHJ2B7eZ3sZ6R3jUZuWhxe3/u1S1V+wwQBm+fnobnVYJO7Ilf8rk0U0dxqYAIvEVE7w2CmHVKqwSLtDrIOdB4ekwGlAr1KgQxgLEBnnUCslgbAipmDFPNc7O2yysvsygReIqJ2hEXz2iF7NVjkAp23i6uc3nYtAEiJNwYRs0ekojh/ArISOsmeO+eWVFMRO40APHpbb3yybCJmj0iFVqdHSWW9TZE+607U1jMwSsX0iIgo9HBmph2KDg+DIMBitkWa1ZALdAyi8T1NrW2qP0MEMGN1CR4ek4GHxmTg1HkdKuquyp676Pabsej2m21mUhwthXELNRERAQxmQpZScq+pBotVILNi5iAAwKWmVtnrORPISAwisO5gFdYdrHJ4rvVWaEdLYUrvIyKi9ofBTAhSmtE4XttgUeQOMK4zbp6fh4q6Kxa7g3xJrvidmjYKREREAIOZkGMdsEgzGpebf0ThrgqbIncGALWXbAvN+YoGkN02zR5JRESkFhOAQ8im0hpMX11iE7C0iaKx2J1MsBImCIBV0GDNWz2XBAAFswbb7ZGklOBLREQk4cxMiJByTJR2P8sdlrY/D0uLs5kF8baHRqXjkXG97QYnTPAlIiI1GMyECLkcE3sEAOvmDsPt/RMAGGvJKCXqejrGEQCHgYyECb5EROQIl5lChLTdWokAy+UiEcAjfz+CTaU1AIBr13zTu0gDoFBhaYmIiMgVnJkJAXLbra2JuBHQWCcHj+3bHQe/rvfqGN+eOwxR4R24XERERB7HmZkgZ12PBTD+oc4bm2FzrgjbJSNpu/PkQQleG+Os3GTc3j+BFXmJiMgrODMT5GQr9gJ484DjQnWSqHCNR/NiZuUm44G8NByubsDw9DjF/kpERESewGAmgJlX8b3Q+ANe2FWBmkvNmD40Gf9zZxYA+XoszlqzvxK7vvjO7fHePzIV/z2slyl4YRBDRES+IIiiE62Mg1BjYyNiY2Oh0+nQuXNnfw9HkXX7AfMqvnI0AjC2Tzfcn5eG+qutpq7U5jkxvhQmCCjOn8BlJCIi8ghnnt+cmQkA5oGLIADzx2dizb5Ku7MtBhHYd6Ye+87UI6NbFJ6bNgBxUeFoaG7Fsx+c9N3gr1tyVz8GMkRE5BcMZjzAfFYFAA5XX4IgCBiWFmd6wMs1ftTq9Hj/SC1e3P2V6VqiCKzaW+nU51fVN+PZD05CADB5sPcSee0Z0quLXz6XiIiIwYybLGZVYLnEI8BYUwWATeNHAFhaVO7RsYgAdpTXefSaarBnEhER+RNzZtyg1ekddpqWy2HxdesAb5J6Js0ekervoRARUQhx5vkdFHVmVq9ejYyMDHTs2BHDhg3DwYMH/T0kAMblJEdBidzLoRLIPDu1P4rzJzCQISIivwr4YGbTpk1YvHgxnnnmGZSVleG2227D5MmTUVNT499xldbgifeO+XUM/qQBMGVIIpN+iYjI7wI+mHnppZfw8MMP41e/+hX69++Pl19+GSkpKVizZo3fxiRXdTfUCVb/XMD+SkREFCACOgG4tbUVR44cQX5+vsXxSZMmoaSkRPY9LS0taGlpMf2+sbHR4+NytkN1sAsTBGyen4dvG/QQRWBYehwDGSIiChgBHczU19ejra0NPXv2tDjes2dP1NXJ79opKCjA8uXLvTouT1TdDRZSgm92CtsSEBFRYAr4ZSYAEATB4veiKNockyxbtgw6nc70q7a21uPjSYyNRMHMwZAfgToKww8YGgF4/b4cJvgSEVHAC+iZmW7duiEsLMxmFubChQs2szWSiIgIREREeH1ss0ekYmzf7lj6/gkc+Kpe9pwwQcBj43sjPioc6d2ioG81QBCA3DTjDEd1fTOiwjVobjVgxfZTKD/v+SUxR2blJuOpO/thfXE13ir+BgbxxmzM3dlJPh8PERGRswK+zsytt96KYcOGYfXq1aZjAwYMwD333IOCggKH7/dFb6bjtQ04XN2A9G5RiArvYApQ0rtFOZVb8vGXddh3+iKG9IpFyzUDDn3zPa61iZjYvwdarhlwobEFt/fvgR6dO+JIdQMEAThcdQnvltbgh2siunS8CTNyklF9qQkXr7Sge6cI/CIvDZeaWrH75HdoaGrFVxeuoHtMR4zt0x335CRZLB1pdXpU1zc7PW4iIiJPc+b5HfDBzKZNm3D//ffjjTfeQF5eHtauXYt169bh5MmTSEtLc/j+YGk0SURERDeEVKPJ2bNn4/vvv8cf//hHaLVaDBo0CDt27FAVyBAREVHoC/iZGXdxZoaIiCj4hFw7AyIiIiIlDGaIiIgoqDGYISIioqDGYIaIiIiCGoMZIiIiCmoMZoiIiCioMZghIiKioMZghoiIiIIagxkiIiIKagHfzsBdUoHjxkbfd6QmIiIi10jPbTWNCkI+mLly5QoAICUlxc8jISIiImdduXIFsbGxds8J+d5MBoMB58+fR0xMDARB8Oi1GxsbkZKSgtraWvZ9Au+HHN4TW7wntnhPLPF+2GqP90QURVy5cgVJSUnQaOxnxYT8zIxGo0GvXr28+hmdO3duN18uNXg/bPGe2OI9scV7Yon3w1Z7uyeOZmQkTAAmIiKioMZghoiIiIIagxk3RERE4A9/+AMiIiL8PZSAwPthi/fEFu+JLd4TS7wftnhP7Av5BGAiIiIKbZyZISIioqDGYIaIiIiCGoMZIiIiCmoMZoiIiCioMZixY/Xq1cjIyEDHjh0xbNgwHDx40O75+/fvx7Bhw9CxY0f07t0bb7zxho9G6jvO3JN9+/ZBEASbXxUVFT4csXcdOHAAP/nJT5CUlARBELB161aH7wnl74mz9yPUvyMFBQUYMWIEYmJi0KNHD0yfPh2nT592+L5Q/o64ck9C/XuyZs0aDBkyxFQQLy8vDzt37rT7nlD+jriCwYyCTZs2YfHixXjmmWdQVlaG2267DZMnT0ZNTY3s+VVVVZgyZQpuu+02lJWV4emnn8bjjz+OoqIiH4/ce5y9J5LTp09Dq9WafvXp08dHI/a+pqYmZGdn4/XXX1d1fqh/T5y9H5JQ/Y7s378fCxYswGeffYaPPvoI165dw6RJk9DU1KT4nlD/jrhyTySh+j3p1asXCgsLcfjwYRw+fBgTJ07EPffcg5MnT8qeH+rfEZeIJOuWW24RH3vsMYtjWVlZYn5+vuz5S5YsEbOysiyOzZs3Txw5cqTXxuhrzt6TvXv3igDEhoYGH4zO/wCIW7ZssXtOe/ieSNTcj/b2Hblw4YIIQNy/f7/iOe3pOyKK6u5Je/ueiKIoxsXFiW+99Zbsa+3tO6IGZ2ZktLa24siRI5g0aZLF8UmTJqGkpET2PZ9++qnN+XfeeScOHz6MH3/80Wtj9RVX7okkJycHiYmJuP3227F3715vDjPghfr3xFXt5Tui0+kAAPHx8YrntLfviJp7ImkP35O2tja89957aGpqQl5enuw57e07ogaDGRn19fVoa2tDz549LY737NkTdXV1su+pq6uTPf/atWuor6/32lh9xZV7kpiYiLVr16KoqAibN29Gv379cPvtt+PAgQO+GHJACvXvibPa03dEFEU8+eSTGDNmDAYNGqR4Xnv6jqi9J+3he1JeXo5OnTohIiICjz32GLZs2YIBAwbIntueviNqhXzXbHcIgmDxe1EUbY45Ol/ueDBz5p7069cP/fr1M/0+Ly8PtbW1+Otf/4qxY8d6dZyBrD18T9RqT9+RhQsX4sSJEyguLnZ4bnv5jqi9J+3he9KvXz8cO3YMly9fRlFREebOnYv9+/crBjTt5TuiFmdmZHTr1g1hYWE2Mw4XLlywiYYlCQkJsuffdNNN6Nq1q9fG6iuu3BM5I0eOxFdffeXp4QWNUP+eeEIofkcWLVqEbdu2Ye/evejVq5fdc9vLd8SZeyIn1L4n4eHhuPnmmzF8+HAUFBQgOzsbr7zyiuy57eU74gwGMzLCw8MxbNgwfPTRRxbHP/roI4waNUr2PXl5eTbn7969G8OHD0eHDh28NlZfceWeyCkrK0NiYqKnhxc0Qv174gmh9B0RRRELFy7E5s2bsWfPHmRkZDh8T6h/R1y5J3JC6XsiRxRFtLS0yL4W6t8Rl/gp8Tjgvffee2KHDh3Et99+Wzx16pS4ePFiMTo6WqyurhZFURTz8/PF+++/33T+N998I0ZFRYm/+c1vxFOnTolvv/222KFDB/H999/314/gcc7ek//93/8Vt2zZIp45c0b84osvxPz8fBGAWFRU5K8fweOuXLkilpWViWVlZSIA8aWXXhLLysrEs2fPiqLY/r4nzt6PUP+O/PrXvxZjY2PFffv2iVqt1vSrubnZdE57+464ck9C/XuybNky8cCBA2JVVZV44sQJ8emnnxY1Go24e/duURTb33fEFQxm7Fi1apWYlpYmhoeHi7m5uRZbB+fOnSuOGzfO4vx9+/aJOTk5Ynh4uJieni6uWbPGxyP2PmfuycqVK8XMzEyxY8eOYlxcnDhmzBhx+/btfhi190hbRq1/zZ07VxTF9vc9cfZ+hPp3RO5eABDXr19vOqe9fUdcuSeh/j156KGHTP9d7d69u3j77bebAhlRbH/fEVcIong9a4iIiIgoCDFnhoiIiIIagxkiIiIKagxmiIiIKKgxmCEiIqKgxmCGiIiIghqDGSIiIgpqDGaIiIgoqDGYISIioqDGYIaIPC49PR0vv/yyX8fwy1/+EtOnT/frGIjINxjMEBEp2LdvHwRBwOXLl/09FCKyg8EMERERBTUGM0TktPHjx2PhwoVYuHAhunTpgq5du+J3v/sdzFu9NTc346GHHkJMTAxSU1Oxdu1ai2uUl5dj4sSJiIyMRNeuXfHoo4/i6tWrptf37duHW265BdHR0ejSpQtGjx6Ns2fPAgCee+45DB06FG+++SZSUlIQFRWFn/70p7IzKH/961+RmJiIrl27YsGCBfjxxx9Nr/3jH//A8OHDERMTg4SEBMyZMwcXLlwAAFRXV2PChAkAgLi4OAiCgF/+8pcAAFEU8cILL6B3796IjIxEdnY23n//fY/cWyJyHoMZInLJhg0bcNNNN+HQoUN49dVX8b//+7946623TK+/+OKLGD58OMrKyjB//nz8+te/RkVFBQBjoHPXXXchLi4OpaWl+Ne//oX//Oc/WLhwIQDg2rVrmD59OsaNG4cTJ07g008/xaOPPgpBEEzX//rrr/HPf/4T/+///T/s2rULx44dw4IFCyzGuHfvXlRWVmLv3r3YsGED3nnnHbzzzjum11tbW/GnP/0Jx48fx9atW1FVVWUKWFJSUlBUVAQAOH36NLRaLV555RUAwO9+9zusX78ea9aswcmTJ/Gb3/wGv/jFL7B//36P32ciUsG/TbuJKBiNGzdO7N+/v2gwGEzHli5dKvbv318URVFMS0sTf/GLX5heMxgMYo8ePcQ1a9aIoiiKa9euFePi4sSrV6+aztm+fbuo0WjEuro68fvvvxcBiPv27ZP9/D/84Q9iWFiYWFtbazq2c+dOUaPRiFqtVhRFUZw7d66YlpYmXrt2zXTOT3/6U3H27NmKP9fnn38uAhCvXLkiiqIo7t27VwQgNjQ0mM65evWq2LFjR7GkpMTivQ8//LB43333KV6biLyHMzNE5JKRI0dazJTk5eXhq6++QltbGwBgyJAhptcEQUBCQoJpCefLL79EdnY2oqOjTeeMHj0aBoMBp0+fRnx8PH75y1/izjvvxE9+8hO88sor0Gq1Fp+fmpqKXr16WXy+9H7JwIEDERYWZvp9YmKiaQwAUFZWhnvuuQdpaWmIiYnB+PHjAQA1NTWKP/epU6fwww8/4I477kCnTp1Mv/7+97+jsrJS1b0jIs9iMENEXtGhQweL3wuCAIPBAMCYc2IeCFmfBwDr16/Hp59+ilGjRmHTpk3o27cvPvvsM8XPk95nfl17Y2hqasKkSZPQqVMn/OMf/0BpaSm2bNkCwLj8pER6//bt23Hs2DHTr1OnTjFvhshPbvL3AIgoOFkHFp999hn69OljMROiZMCAAdiwYQOamppMszOffPIJNBoN+vbtazovJycHOTk5WLZsGfLy8rBx40aMHDkSgHH25Pz580hKSgIAfPrppzbvt6eiogL19fUoLCxESkoKAODw4cMW54SHhwOAabZJGntERARqamowbtw4VZ9FRN7FmRkickltbS2efPJJnD59Gu+++y5ee+01PPHEE6re+/Of/xwdO3bE3Llz8cUXX2Dv3r1YtGgR7r//fvTs2RNVVVVYtmwZPv30U5w9exa7d+/GmTNn0L9/f9M1pPcfP34cBw8exOOPP46f/exnSEhIUDWG1NRUhIeH47XXXsM333yDbdu24U9/+pPFOWlpaRAEAR9++CEuXryIq1evIiYmBk899RR+85vfYMOGDaisrERZWRlWrVqFDRs2qL+BROQxnJkhIpc88MAD0Ov1uOWWWxAWFoZFixbh0UcfVfXeqKgo/Pvf/8YTTzyBESNGICoqCrNmzcJLL71ker2iogIbNmzA999/j8TERCxcuBDz5s0zXePmm2/GzJkzMWXKFFy6dAlTpkzB6tWrVY+/e/fueOedd/D000/j1VdfRW5uLv76179i2rRppnOSk5OxfPly5Ofn48EHH8QDDzyAd955B3/605/Qo0cPFBQU4JtvvkGXLl2Qm5uLp59+WvXnE5HnCKJoVhiCiEiF8ePHY+jQoX5rWfDcc89h69atOHbsmF8+n4gCC5eZiIiIKKgxmCEiIqKgxmUmIiIiCmqcmSEiIqKgxmCGiIiIghqDGSIiIgpqDGaIiIgoqDGYISIioqDGYIaIiIiCGoMZIiIiCmoMZoiIiCio/X8l8zHyEI7hpQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(x,y,'.')\n",
    "plt.xlabel('phosphate')\n",
    "plt.ylabel('nitrate')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "Create a subset where both variables have finite values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "ii = np.isfinite(x+y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0       True\n",
       "1       True\n",
       "2       True\n",
       "3       True\n",
       "4       True\n",
       "        ... \n",
       "2343    True\n",
       "2344    True\n",
       "2345    True\n",
       "2346    True\n",
       "2347    True\n",
       "Length: 2348, dtype: bool"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ii"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Method 1: Scipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = stats.linregress(x[ii],y[ii])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LinregressResult(slope=14.740034517902119, intercept=-3.9325720551998167, rvalue=0.9860645445968044, pvalue=0.0, stderr=0.052923783569700955, intercept_stderr=0.10258209230911229)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "14.740034517902119"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result.slope"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Exercise: Draw the regression line with the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x14b228820>]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGdCAYAAACyzRGfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAABk6klEQVR4nO3df3xU5Zk3/s85o4QKyLcq8uMZiGh47A+LqxBI8AfZqKQVlPpaMOhCQxshirpASCSIu+6uCvkBpLRKIMm0WKnCWLWKuovsMkXXBBNc+yju2gZLwkyVB9r6BEhb3M7c3z/wDDNn7nPmnJkzk8nk83695lVy5syZMydjz5X7vu7rUoQQAkRERERpovb3CRAREdHgwuCDiIiI0orBBxEREaUVgw8iIiJKKwYfRERElFYMPoiIiCitGHwQERFRWjH4ICIiorQ6r79PQC8UCuGTTz7BiBEjoChKf58OERERWSCEwKlTpzBu3DioqvnYRsYFH5988gnGjx/f36dBRERECfD7/XC73ab7ZFzwMWLECABnT/7CCy/s57MhIiIiK06ePInx48eH7+NmMi740KZaLrzwQgYfREREA4yVlAkmnBIREVFaMfggIiKitGLwQURERGnF4IOIiIjSisEHERERpRWDDyIiIkorBh9ERESUVgw+iIiIKK0YfBAREVFaMfggIiKitGLwQURElKBAIACfz4dAINDfpzKgMPggIiJKgMfjQW5uLoqLi5GbmwuPx9PfpzRgKEII0d8nEenkyZMYOXIkent72ViOiIgyUiAQQG5uLkKhUHiby+VCd3d33Hby2crO/ZsjH0RERDYEAgF4vd6owAMAgsEgDh8+nNRxB8sUDoMPIiIii7SpllWrVsU853K5kJeXl9RxB8sUDqddiIiILJBNtWhcLhe2bduG8vJyR447EKdwOO1CRETksK6uLmng0djYiO7u7oQCDwBoa2tzfAon0zH4ICIismDSpElQ1ejbpsvlwrx58xIeofB4PFiwYEHM9mSmcEydOAGcfz5w883OH9sGBh9EREQWuN1uNDc3w+VyATg31ZJo4BEIBLB06VLosx9cLhdqamrQ3NyMpqYm5xJQn30WuPRS4C9/Af793505ZoKY80FERGRDIBDA4cOHkZeXl1DgEQgE0NXVhRMnTqC0tDTm+cmTJ+P999+P2tba2ory8vLwaydNmgQA4X+bnocQwDXXAP/n/5zb9sQTwMMP2z53M3bu3+c5+s5ERERZzu12JzXNsnTpUoRCISiKIt1HH3gAwJIlS/DZZ5/hoYceCo+UKIoCIQRUVUVzc7M852T/fqCoKHrbhx8CX/taQufvFI58EBERpYHZaplkuVwutLe34/Tp0+dGQvTBzf/6X0BPD/DFtJHTOPJBRESUBpHTIPFGQ4xWyzghGAxi+vTpEEJguKLglH5c4a//Gti3LyXvnQgmnBIR0aCXSHXRyMJg48ePR0VFhenrN27c6MSpGhJCoB6IDTz27cuowAPgtAsREQ1ykXkYpvkTEYymUBRFQX19PaZMmRI1GrJ27VqsW7cuZZ8BAKQ381AodvolRezcvxl8EBHRoJVodVGfz4fi4mLTY2uBTElJCcaPH+/YOevNBPAL3bY/A/id35/WCqmscEpERGSBLA/DSnXRSZMmGa5W0YRCIVRUVKCtrS3p8zQiEBt4FAL4YX19RpdmZ/BBRESDllHV0njVRd1uN+rq6uIePxgMQlGUuIFKImTTFgqAAwB+85vfOP5+TmLwQUREg1YyVUurq6vR0NBgGli4XC5cdtllTp0uAGAljAMPzbZt25yrjJoCDD6IiGhQKy8vR3d3N3w+n7RBnGwljLZtwYIFOHr0KJqamqTHrq2txenTp2NKqCdKANik2zYD0YEHcHblS3t7uyPvmQqs80FERIOeUdVS2UoYADHbrrzySulxc3Nz8dFHHyV9fgoAWYWQ9KxjcR5XuxAR0aBgpyCYtr9sJUwoFIoaydCqi2pFviKpqpp0YTGjm7RZ4KEoCo4ePcrVLkRERP0lsiBYbm4uPB5P3NcYrYTRBxjBYBDPPvusdGolFYHH/0Js4DF16tRw4qyqqmhpacno1S4c+SAioqyWaC0P2eucGMmw4goAssW+stEO7bMASKrbbrI48kFERPSFRGt56FfCqKqKG264IWXnqRGwF3hoq3PcbjeKiooyesRDw+CDiIiymqyWBwB0dnbGfW15eTnWr18PRVEQCoWwf//+VJximGwq4nzIA4+ysjLp6pyBgMEHERENKHabwBkVBFuzZk3cYwQCAdTU1Di2VNbIChjX7viLwWsKCgoGxCiHDIMPIiIaMDweDyZMmIDi4mJMmDAhbuKoFqhMmDAh5rl4Uy+BQABerzflOR4CQKNke7xltFOmTEnB2aQHgw8iIhoQAoEAlixZEh6FEEJgyZIlhqMXkStcFixYIN3nv//7v6Wv1167atUq5z6ATnFxseFoh5X6HX19fabPd3Z2YtOmTZaml9KNwQcREQ0IbW1tMdMf+kqe2khHZ2dnuBCYtp/MsmXLYpbeBgKBqNemwikA/75vX8x2O0XD/u3f/s3wucWLF2PatGlYtWoVpk2bhsWLF9s+x1Ri8EFERFkhcqSjoKDAcvAQCoWwdOnS8AiIbHWMVS6XCyUlJdIEV40AMFx/7rBfrXT9+vUIBAIxOTCdnZ14+umno/Z9+umnM2oEhMEHERENCDNmzIhp4qaqKgoLC2NGK2TBg6Io4WWzeqFQCJs3bwZgvDomHkVRsGzZMrzxxhuGwYvRNMs9tt/t3Dnri6e99dZb0v3ffvvtBN4lNRh8EBHRgOB2u9HS0hLVgba5uRlut9twtEILIlwuF1paWtDd3Q2v1ys9fmNjIwKBQLi+h90ARAiBH/7wh9IpHoH4nWjtcrlc2LhxY1TAVVFRgf/9v/+3dP/rrrsuiXdzFoMPIiIaMIw60MpGK1wuFw4cOBC1r9vtxvz581FSUhJz7MjVL+Xl5ejp6cGNN96Y9DnLgg7v0qVwJTC6onG5XFi5cqW01Pvw4cNRVlYWtb2srAz5+fkJv5/Tkgo+tMIrK1asCG8TQuAf//EfMW7cOHzpS19CUVERPvzww2TPk4iIKEy76Wr5DgCiqpFqlT/z8/Njqn4GAgG88cYbMcdUVRV5eXkAzuZNeL1ew1UyVoyFPPDo7OjAndu2obKy0vKxVFWF1+tFR0dHOJhavny5NODKy8vD9u3b0dHRgcbGRnR0dGD79u0Jf46UEAnq6OgQl112mZg8ebJYvnx5eHttba0YMWKEeOGFF8QHH3wgSktLxdixY8XJkyctHbe3t1cAEL29vYmeGhERZanW1lahqqo2iyEURREAhKqqorW1Vfj9fuHz+YTf7zc8RkVFRfj1kY+FCxeKffv2ifnz50uft/MQBo/W1tbwefj9fkvH0j6b0fVwuVwCgHC5XIb7pYOd+3dCjeVOnz6Na6+9Flu2bMHjjz+Ov/qrv8L3v/99CCEwbtw4rFixAqtXrwYAnDlzBqNHj0ZdXR0qKiriHpuN5YiISEbW6C2SvllcIBBAV1cXJk2aFLVtwoQJKa1YKjvyW6+/jonf+EbUCIzH44mqW6IoCurq6pCfn49hw4aFm8UVFhbGbYDXnw3lNHbu3+cl8gb3338/Zs+ejZtvvhmPP/54ePuRI0dw7NgxzJo1K7wtJycHM2fORFtbm6Xgg4iIKJIWRJw4ccJ0CWwwGER7ezsuueQSvPvuu1i9ejVCoRAURcGaNWtw9dVX49ChQykLPBYB+Ilku0tV0aMLPLTVOZHnoigKioqKcPr0aYwdO9Y0R0MfWA20Muu2g4+dO3fiP//zP6XrhY8dOwYAGD16dNT20aNHo6enR3q8M2fO4MyZM+GfT548afeUiIgoS3k8nvASWlVVoSiKafBw5513xmwTQmDdunWpPE3paAfwxWqWUAiHDx+G2+02DaRCoRCmT58OIQRUVUVzc7O0aZz+mhjtl8lsJZz6/X4sX74cO3bswNChQw3306/DFkLEbNOsX78eI0eODD/Gjx9v55SIiGiAsdoYTla7w+he0p/ilUjXkkAji6Dddddd8mN9EVhpy2b110h2TWT7ZTpbwce7776L48ePY8qUKTjvvPNw3nnnYf/+/fjBD36A8847LzzioY2AaI4fPx4zGqJZs2YNent7ww+/35/gRyEiokwXeQPWlzXXk9XuSFcAUlxcjGXLluGJJ54w3OcQ4tfuUBQF27ZtAwDbgZSs8Z3smsRrkJeJbCWcnjp1Kmb65Lvf/S6+8pWvYPXq1fj617+OcePGYeXKlXjooYcAAJ9//jkuvfRSJpwSEQ1ysoRRsyRRAKYJpukim+qR3Tj/HcDNET9XVFTgkUcegdvths/nQ3Fxsa331V8bwNo17C927t+2Rj5GjBiBq666KuoxbNgwXHzxxbjqqqvCNT/WrVuHl156CYcOHcLixYtxwQUX4O67707qQxER0cAW7692/ajInj17omp3JFLy3AlWAg8F0YFHU1MTtm7dGg4IjIqgGdHqlOgDCrfbjUWLFkVtW7hwYb8HHnY5/pt86KGHsGLFCixbtgxTp07Fb3/7W7zxxhsYMWKE029FREQDiNENOC8vzzCXoaSkJFzR9MCBA/0WgAD2SqR/+umnUXkYWsl2fRG01tbWqG0NDQ0x1VsjBQIBPPPMM1HbduzYMeByPhIuMpYqLDJGRJS9jIpi7du3T1pgy+fzGb4+nQ9ZwbAFXxQ4M3pEFj7bt2+f8Pv90iJoVgqjaaxep/6Q8iJjqcScDyKi7CYrimWUy9De3o7Tp0/HFAp79dVXcd9996X8XP8/AJ9Jts+88Ub89Kc/xZ49e6IKhekpigJFURxbFmtUaK2hoQFVVVUJH9cJKcv5ICIiSpbb7Y7ptyKblli4cCEKCgqiVsZoCakXXXRRys9TQB54KADefPNNAGcb0B09ehRNTU3S1StCCEeXxbrdbtTW1sZsr6mpGVBTLwlVOCUiInJaeXk5SkpKcPjwYZw+fRq33357VN0L/QhDvIJjyZAddTSA4xE/t7e3Y/78+XC73bjyyistnYuWYGslQVRWHh4Apk6dmtRxMwFHPoiIKGO43W58/PHHmDt3buwqkzg/O+E2GCeVHpds1wqmDR8+3FIyrKqqGDZsWNz9GhoaMGHCBGk9FLPE3YGCwQcREWUM/aqXdBIAXpFsl61mUVUV3d3d4aXB06dPxy233BI3AAmFQigoKDAtrrZhwwY89NBDhtVOs2G5LRNOiYgoY3i9XpSWlqb9fa0sodWmeVwuF2pra8ON66zQTxEZFQYz69zr8/lQVFSUsYXGUt7VloiIyEwgEEBbWxsAYMaMGXFbwnd1dYU70abTvwG4SbJdNtqxa9cujBo1Cnl5edKCaWb0f+cb5WgYHVdV1fC0ilmxtoEy+sHgg4iIHOXxeKKSQxVFQUtLS9wOrekmG+0IAJC1N3W5XCgsLIy6uauqmvB5G+VoaPkc+uPW1tbGVEvVj3ww54OIiAYlLWcj8i99IQSWLl0at0OrHck2l5MFHi5VlQYeqqpKS51XVlYm9N5GpdOB2CXHqqqioqIiqguuUbXUgTLqATD4ICIiBxlNG4RCIUsdWiOZJW8mmq5oVCL9vC/yOGTvKYTAZ599Bp/Ph0AgEO5Bs2HDBtvvrygK2tvbTQuNlZeXo729HXfeeSdCoRC2bdsWs+KlvLw8XHbeqBR7JuO0CxEROcZo2iAyZyFyX6NaHaqq4lvf+hZee+01x85NFnQcrazEb267Dd1fVFu96KKLYkZjhBCorq4GcG7EJeHgRwj09fWZ7iObiorsdaONcLjd7gE12hGJIx9EROQYbUogcgRBURQ0NzdLO7Tecsst0uOEQiHHAo8vwXg1y8Tvfx8ff/xx+NzKy8vx3HPPGR5LCGEp8FBVFfPnz4/ZHi83w2wqKrID8EDH4IOIiBKiFdjS53KUl5ejp6cHXq8XXq8XR48eRUlJSfhnbf9AIIA33njD8fNSVRVz5swBcDbo+KNkHy1jRFbyfMaMGUl3z33qqafg9XpRX18fPpaV3AyzqaiBllRqhtMuRERkW+TUgKxhmtvtDv/lr1/9AgDXXHMNRo4c6eg5qaqKp556ChdddBF+//vfY/err8bsczmAI7pt+mWq2uhNRUUFgsFgzDEim8UZufjiiwEA1dXVuOuuu2Ia6Rkxm7YaaEmlZlhkjIiIbLFT5CoQCGD8eNkaEufNnz8fL7zwAq4LhfCm5Hmj9TFm53748GF0dnZizZo1CAaDUFUVtbW1KCoqwrRp06THU1UVPT09UV14ZT1ajHg8nnDg43K5sHLlSixfvjzjAw8WGSMiopSxU+Rq8+bNaTuvF154AUGD0QizwMNs2avWgVdRlHBF05qamqilr5G0USDtePFGiGQiG+xZGS0ZiDjyQUREtlgd+TArFZ4KspvZ87t2IfDb30prcjQ2NmLevHlxb+5WP4eqqjhw4ADy8/MNX5cJZdBTxc79mwmnRERki9UiV3ZLkCfqhzCu3VE4Ywbmz58v7QI7b948AJAmzUay+jlqa2sxduzY8PHMRogGOwYfRERkm1GRq8gVMLLW78lSFAVLly4N/ywAPCDZ77yIgEhWNXT9+vXYs2dPuCutvohXJKufY/Xq1ZgwYUL4eAcPHpQGPdmyYiUpIsP09vYKAKK3t7e/T4WIiGxobW0VqqoKAEJVVdHa2irq6+u1oqJJPyoqKoTf7xd+v18oiiIEEPNoamoSPp9P+P3+mPNraGiIOj9FUaKO73K5pK/TPpvL5bJ1vqqqioaGhvDrXC6XaG1tdfy6+/1+sW/fPsNzT3Rfu+zcvznyQURESdMXx9LqZ+Tm5jpyfEVRcMUVV5wdyRg/HiFJuqKCs0tci4qKpCtXtIRR7fyEQadZAOjs7MSmTZvQ2dkJIHqkR5uuiScUCiE3NzelZdC1Uu/xRm/s7ptyjoc+SeLIBxHRwLNv3z7pX/9NTU2OjXxAMtIhAPEPXzynKIrhX/RG5xf50EY+ysrKoraXlZVFHcvv94stW7aIlStXhkdSjB5erzdl19zv98e8v9HojZ19E8WRDyIiSgstx2P48OHS/IadO3c68j4uyJNKVUXBP0f8vGfPHunrhw8fHrcT7vr16/Hpp5/i6aefjtr+9NNPh0dAtNGDZcuWYfPmzVi0aFE4l0RPURQUFhaavmcy7CS0ZlryK4MPIiJKSOQwfkFBAQoKCqKeLykpwf79+5N+HwHgL5LtAb8/KqAQQqCiogKdnZ1RK1g8Hg+mT58etyfL8ePH8dZbb0mfe/vtt6VTSzt27EB7ezt8Pl9UKXVVVdHS0pLSJbWyRFijhFY7+6aFY+MtDuG0CxFR5pMN4+sf+oTORB6yaZb9mzeLffv2iV27dhkmemr/29DQYPk8XC6X2L17t/S5jo4Ow6kbn88XdV2MEl5TITIRNl5Cq519E2Hn/s0iY0REZJvP50NxcbHjx1UUBUIIfAPA+7LngXDvE1kPFD0r+0Ty+XzYvn171NRLWVkZtm/f7ljRMLvl1q0cz2o1VDv72sUiY0REZMioG60dTtfwmDx5Mjo6OtDS0gIB48ADQNS0Rzx2Ag9tGmL79u3o6OhAY2MjOjo6sH37dgDWi6uZScWKE60EvJXzsLNvSjk65uIATrsQEaWOrBZHMsfShvGtTG1861vfMn1eVVXpNMuQL44db5pH9nC5XJZqjdiZhrA7taLV1ujo6Ej5ipP+xGkXIqJBwO7wfbLTBoFAAG1tbQCAGTNmAAB2796N119/Ha+99lrchM7du3dj7ty50tGIxwA8IntPvx8AcPjwYQwbNgwFBQW2RjNmzJiBt99+G4888gieeOIJw/0qKiqwdetWy8e1KrKxnDalpOfz+VBUVOT4e6ebrft3igMh2zjyQUQUXyIjGFYSJiP5/X6xa9cusWvXLlFfX59UAumtt94qdu3aJSoqKmKOIxvtEBGjIZGfLfJzW31897vfjbuPqqqOj0BYScrlyEeG4MgHEZG5REcw7LzO4/FgyZIlcUczrDL6q192dH01DpfLhfb2dpw+fRqTJk3Cp59+ansExAqnRyCqq6uxYcOGmO1aEqyWM+J01dN4nE541TDhlIgoiyVaMMooYRKI7uwaCAQcDTwAxBxL+9NfT1YGLBgMoqCgIJyk+f7776O5udnRhFena14EAgFs3LhR+j4vv/wyNm3ahPb29rQHHplSYp3BBxHRAJNMwSh9N1oAMTejrq4uRwMPPdmR//PrX5cGHprIFS5Lly5FT08PPB5P3KqlZrRrmMiqlXiMruHNN9+MuXPnorKyEgUFBWm5+Wurmzo7O6X9d5JZ9ZSwFE7/JIQ5H0RE8SVTMCre6ouOjo6E8joWLlwo/v7v/950H1luxze+8Y2Ec0lmzJhhu9ts5OdMVUEwWb6HoihpX+0SmSNjlLNjlPNjF3u7EBFlOf0IhtXh+8hh9+nTp0unb959992Ezunaa6/FtGnTpM+ZTbN88MEHCb0fALS1taG1tdU0V8PlcmHevHkxIx35+fkpq3khm+JatWpVWvur6MvBC8lITH+VWGfCKRHRICFLONVzuVx4+OGH8dhjjyX0Hqqq4pprrokKYGQ3mZkA3jQ5jpag6nK5MH369PASX7vnsnPnThQWFsLtdqe0uqeRyPcE4EiFVKuMqtCmKuHVzv37PEfekYiIMp4sURWIvRlNnjw54eAjFArhvffeAwBcBuCIZB+XhZLn4osmceXl5Zg+fXpC51JZWYn58+eHf3a73Wmv7Kl/z+bmZlRUVCAYDKYk1ySSlhukD3ba29vR19eX1iBMjyMfRESDRCAQwPjx42O2P/DAA7j88ssxadIkDBs2DH19fbjtttuSei+jG4uqKHjllVdw++23W0pqveaaa8LBjB12RxRStfzU6L3SNQLj8Xhigp1UrbBhkTEiIhJCnEsu9fv9Yu3atQkndtp5yJJKL9AlOEYmzKqqKm666Sbb7/Pggw9Ki3jZTcB1suR8JkpXp107928GH0REWcjv94uqqqrwCgcn2tvHe1QaBB76/RoaGqLOMZGeLVqQUV9fHxXEVFVV2brJylalZFPV0XSyc/9mzgcRUQZKZhogsp+IRqR4ht3o6LIqHDU1NViwYAEAYNOmTQlXKg0Gg8jPz0d3d3e498vp06dtHcOsYFu/d37NYlxqS0SUYZKpQqlfXpkORktotzY14b777ot5Tru5GyXAWqUtE3W73fj444+jqqBavWbJFGyjxDHhlIgogyTTt6WrqwsnTpxAaWlpOk4VpwAMl2zXRjv0Ky002ucBYpeeqqoKcTYlIPqYioJ58+bhxRdfjEmeTLZbbzqTMrMZl9oSEQ1QiUwDRE6zqKpq2MQtEd/85jexZ88eaW8WvQ4A2qJYRVEMl/WuX78ewNkCYffccw9aW1ujlvoCkE4bvfjii9JloslOnZSXl6OkpCTtNUAGM458EBFlELt/xcv2Nxo9cIrRNItm+vTpeOedd2L3+SIo0gdHiqJg1apVWL58OYCzwcRHH32EZcuWxRzD6/VG1e4AEh8tImexqy0R0QBl1HnW6CYq+6s/FAqlJPDQloPoKQAefPDBcO6ELPAAziW9xoyiCIHGxkY899xz4VyX+++/X9o0bsGCBTH5HHavGfU/jnwQEWUgq4WorJRMd4LsRnEXgJ1f/NuJqR7ZiIhs+sZoVKM/yqfTORz5ICIa4Nxud9ymZ1qS6Zo1a1J2HpfCeLRjZ8TPyQYe2lRRJCEEHnnkkZh9jZqxWblmlBkYfBARpUggEIDP50MgEHD09YFAANXV1eEpinXr1jlxujEEgP8r2S6r3ZEMVVVRV1cnXfI6Z84cLoXNQgw+iIhSIJlaHbLXb9iwIRx0jB8/Hhs2bDBtlZ4s2RFHQR54uFwuLF261Nbx7733XjQ1NeEf/uEfcODAAVRVVUnzNvLz85nPkYWY80FE5LBkV1+kK49D5nsAZGGSFnQ0Njbi7bffxs9+9rPwc2VlZbj//vsxbdo0y+8Tmc+hqiqam5vDNTtkeRvM58h8rPNBRNSPkq07kWzlz0TFK5Hucrlw3XXXYdWqVVHP79ixA7feequ994pYChwKhVBRUYGSkhLDtvdG2yNpOTDDhw/H6dOn09KhlhLDaRciIoclW7J70qRJ0mWmqWSUVBp5FgsXLsSRI0ekgZWiKDGf2Q6jJFIZWS5M5DTVtGnTEp7uovRg8EFE5DAn6k6ka0b8MOSBh+uLSqmRnnnmmXBDuEiqqqKwsDDqM8v2MQtOrAZnslwao3422ohKogm/lDoMPoiIUqC8vBzd3d3wer149tlnUVJSYvm1bW1tKTyzcwSAKyTbn/d68dxzz8UEQEbFy4QQ2LNnT/gz+3w+rF27NhxsuFwuNDc348CBA9IARFVVS8GZPsjQgou2tjbDaSo7IyqUPsz5ICJKkT179kT1XNGSKjVajkIqcxNmzpyJ/fv3x2w3Gu0IhUJQFyxAbW2tYWO4mGMJEZOzkZeXh8mTJ0NRFBQWFsLtdsPn80mPt3PnzpiS6TJGuTTalI9REzsuy81AIsP09vYKAKK3t7e/T4WIKGF+v1+oqqpVJBcAhMvlEn6/XwghRGtra/h5VVVFa2tr1GsVRYl6rVMPYfDQn6uqqqKhoUG4XK7wz/HOyefzmX62eNckmWva2toaPtfI5yKvayK/w3379lk+v8HOzv2bwQcRUQrs27fP8AZtdBPdvXu3WL58uVi+fLm444470hJ4rALEP/zDP0j3LykpCQcciqKIsrKymPOOfHR0dMQNMMrKyqKeKysrs3VdI4MMfXDh9/uFz+cTHR0d4eucKLPgkORSFnxs2bJFfOMb3xAjRowQI0aMEAUFBeL1118PPx8KhcSjjz4qxo4dK4YOHSpmzpwpDh06lLKTJyLKVLKbsKqqYteuXWLXrl0pGdUweowwCDy051esWGHpOIqimI5++Hy+hIIuu0FCR0eH2LRpk+jo6Ejb7y6R8xxs7Ny/bSWcut1u1NbW4uDBgzh48CCKi4sxd+5cfPjhhwCA+vp6bNq0CU8++SQ6OzsxZswY3HLLLTh16pSdtyEiGvD0K160pmmlpaVYsGBB2pbSCgAnJdsj3/373/++tWNF1ObQ03IrzJYZm9U/scrj8aCgoACVlZUoKCiwtZTWarl7J86T4kg20vnyl78sWltbRSgUEmPGjBG1tbXh5/785z+LkSNHiq1bt1o+Hkc+iCib+P1+0dTUJB0xiJzS0D+X7GPFihXS0Y7LUjCyoihK1LSE0dRIKnM+4rEzjcKRj8SkbOQjUjAYxM6dO9HX14fCwkIcOXIEx44dw6xZs8L75OTkYObMmabLxs6cOYOTJ09GPYiIBop4f03v2bMHy5YtMxwx8Hq9ePzxxx09p9kAGiWjGQqAbhvHufXWW8MjN6qk7odGVdWopcSRS267u7vDK3yM6p8ASOmIhNESXaP3c6JOC5mzHXx88MEHGD58OHJycnDvvffipZdewte+9jUcO3YMADB69Oio/UePHh1+Tmb9+vUYOXJk+DF+/Hi7p0RE1C/iNY8LBAJYsmSJYeAhhMCoUaMwZMgQx85JAHhVsj2RSZ7FixeHg4ienh60tLRI63TIAgCj9vb6wASA5QZ8ssqvqqrGXUqbSNBiFECRQ+wOq5w5c0Z0dXWJzs5OUVNTIy655BLx4YcfirffflsAEJ988knU/vfcc48oKSkxPN6f//xn0dvbG374/X5OuxBRxrMyNG8lsfS6664zXUFi5yGbZknmeF6vN+Zzd3R0SM+3oaEhJddQv79+ikpRlLjTIZxGSY+UTrsMGTIEeXl5mDp1KtavX4+rr74amzdvxpgxYwAgZpTj+PHjMaMhkXJycnDhhRdGPYiIMp1TSYlvv/120k3k/h1n76h6Vkc7Fi5cKB1RKCwsjNk3Pz8ftbW1MdtrampslzG3ew27urpiRpGEEHGvOadRMk/S5dWFEDhz5gwmTpyIMWPGYO/eveHnPv/8c+zfvx8zZsxI9m2IiDJGIBDAiRMnYm7Y2qoOLQ/kggsuSNk5KIqCjRs3QgAolj1v41hjxoxBS0tL1M25ubnZ8OY8derUmG2JBF52G/Al07BPP41SUlJiKc+EUsTOkMqaNWvEm2++KY4cOSLef/998fDDDwtVVcUbb7whhBCitrZWjBw5Urz44ovigw8+EHfddZcYO3asOHnyZEqGbYiI0i1y1URk3QttVUd9fX3KqpMiYsqgtbXVsWkWbQpCK9KlTUcYVfh0chrDrGiY1f3tViJlAbHUSFmRse9973siNzdXDBkyRIwaNUrcdNNN4cBDiHNFxsaMGSNycnLEjTfeKD744IOUnTwRUaISKZ0tu+lqD0VRxLx581IadAAQjY2NhiXSkzmuVhpdE+8GbTdoiHdd7VQkjdw/3nnqf8/M/0gdllcnIjKR6F++RtU70/VQFEUadGyyeQz9yIz+5itLKpXdoK0GDanqkWK3f05DQ4PYuHGj9Lp4vV72cUkSgw8iIgPJ/OVrNvKR6scQh0Y77rzzzphRi/r6+vCNt7W11XDaSD86YkUqpzjslnI3ekQ2zeM0TOLs3L8VIQwWoPeTkydPYuTIkejt7eXKFyJynM/nQ3FxbIqmz+dDUVFR3Nd7PB5UVFQgGAym4OzkjP5POpHaHR0dHcjPz0cgEMDhw4dx8OBBrF69GqFQCKqqGpZQV1UVPT09tlaIdHZ2oqCgIGpFi8vlQnd3tyMrTQKBAHJzc6XH7+rqkv6e9VwuF0KhUNRndvIcBxM79++kV7sQEQ0kVldMGFUuLSkpwbPPPguv14v6+vqoCqCpIAs8/gqJBR6zZs3C2LFjAZxdfpqXlxcOPADE3ISjzsPm36kejwfTp09PaY8UsyW0st+zXmNjI5599tmYz8Y+LmmQwhGYhHDahYhSLV6ypNFUgX57WVlZyla2XJfkNMsdd9whfD6f2L17tygtLZV+Hrs5LFanXcymPFKR3GmUexL5ezY6DyagOoc5H0REcRjdsIxuRlu2bElbvocs6LATeADnqpOa3VyNKobKPqedG7JRUNMf+RTa77m+vt4w4HRy5c5gxpwPIqIEGeWEpIvs/5DPA2Anw0RRFBw9ehRut9s0xyUvLw8TJkyImnZQFAXvvPMOvF4vNm3ahFAoFJ7OsNrfRJaLoaoqDhw4gPz8fBuf5Oyxurq6MGnSpKRzMLQ8l7y8vJhjmT1H1jDng4gGvXjdZo32Hz58eMryN8z8FPLAQ4G9wEOzZ88eAOY5Lkblyvv6+tDQ0ICenh7bjdW0YKG2tjamYqpR4KH/XWk/NzQ0WG46Z4VRs7t4z1EKpHYQxj5OuxBRsuwu75TlchjlCqTiYWea5bbbbpNu10+fqKoqOjo6hN/vF1VVVeHPFzmt4HS+g6yuRrw6ILJr72S+SKpqjFAs5nwQ0aCVSKdU2f4dHR3C6/WmvFS61aDD7GFWdTWyfkVVVZVpUqa+5kcqr7vRa+I97NQaYRn19EppV1siokyWSKdU2f59fX2YP38+6uvrU3Ke2t1UL5EltEuWLJFOFYmImh2hUAiNjY0x+0Q2XKutrUVNTU1C0xyJdPmVvcaM1SZywNmpm6VLl0YtI66oqGAjuQzB4IOIsooTnVJVVcV///d/w+v1YsKECY6foyzo2I3EAg9VVfHII4+gubk5bq6KUTBgVPPDzs06kY6zVmpxRB5Lq+FhRSLBEKUPgw8iyipmhafM9leUc7f+UCiEZcuWobS0FAsWLHDs3BQYj3bcnsDxtEROt9uN8vJyHDhwwPRmbhYMJHuzll33lStXxn1NZWWl6fk2NDQklPR64sSJqN+pdjyrIyeUYimfBLKJOR9E5AQrTc+0ZMS1a9emNK8DBrkdieR3aI/GxkbpZ4vM4Yis2aGvX5Gqbq9agqvVXimyWiMul0t4vd6EkkQj8zz0xy0rK7N9PLKOCadERBKRN9zIm1R/BB6zkjxmfX296efUAi9ZEGZWwTXZYluJBDFOFfmKl8DKyqWpxSJjRDTo6YtTeTyeqARERVFs9yux6+sADkm2J5LboZdo87NAIBBTWCzyWHaKbckKgCXauM+JIl9WCsRZbSBI9tm5f5+XpnMiIkqbyEBDVVXU1dVFJVMCSHngYXR0JwIP4Fw+ht0b9ebNmw0bqbnd7vAjHv01bm5uRnl5OYYPHy7df9iwYabHs/q+ZrQEVqMVNMz5yBxMOCWirCJbYqkPPFJNFnhcAOcCDyCxG2kgEMCmTZukz9k5ltky1tOnT0tf09fXZ+tczd7bqHKtPulVVdVw0qnd1TKUWgw+iCiryFZtpCvw+AGMV7P8yeaxIldqKIqCefPmWV7BY8SoroZ+VUgix9FGT4yWLh8/fjypGhuBQADV1dWYMGGCaR2SyLolPT09OHr0qO3VMpR6zPkgoqwia2qWjvwOp6dZ9OfscrnQ3t6Ovr6+hPMiZNdGYycXQnacyLwRj8eDiooKBIPBcGAjhIianjE7tj6PRJ+vI3tP6n9sLEdEg5Zs6D3VjEY7kplmkeVl9PX1JdX8zO12o66uLma7fgonXlM+2TVev359+Ly00Qev1xsVRMUrXObxeGIayemneCKxaNgAlpoFN4njUlsiisdKszC/3y+amprEokWL0rqEVtg8hqwhXH19vWMN32TXKvL4+qWtdvqhNDQ0mO67b98+6WeW9WcxWqK7a9cuLp0dIFjng4iyVrybo3azvfXWW1MWdBgFHr9N4Diqqgqv1ys6Ojqi6nEkW/siXrEvWf0POzU6rOxrVHejoaEh5nhGgYrX65Ueg43iMg+DDyLKSrJqmKqqRt2w01E4TBZ4JHM8o06t+gDBaMRHv93oOsQbKbAzUmF13/r6ekvnYRbM6AMxWXde6n8MPogoKxkNwWuluPWBSTqCjniBR2SJc1VVpaXErdxIzaqSRm6PnAqxE+gI4fzIhxD2Ahqz0R4r5fKpf9m5f7PIGBFlhd27d6d0RYvsyEsAtMZ7nRDYsmUL/vSnP+H666/H+++/H14JYnXJrFFdjc8//xwPPPCA5Zom8WqDaImkVs7P6r6ywl9G51FeXo6SkhK0t7dDCIEZM2ZEvR9XtWSRVEdCdnHkg4iMyEY3FEURfr9f/P3f/31KRjsmJDDaoX/oRyzs/hVvNHpg9JCNANnJkejo6BCbNm0SHR0dln4n8T6LnfwVOwmvlFk47UJEWcvo5tTR0ZER0yzxHoms0JDlTZgdv76+PiZHoqOjI+4KIbPrq7Gy0kjGapdho6mcRN+X0ofBBxENGIncVGQ3soaGhpQHHhc7dGyzvAvZZ7WTRKsoSszoitXRhHh5HKkelTAa4amqquJoyADA4IOIBgSnbmZ2b9Bmj9UOj3YkmmCqsTvlon8Po4BCNhJilhwab6WRE4zOVfa+VqeEOFqSPnbu36xwSkT9wqw5mV1tbW2O9G8RAGol2xOtVOpyuVBXV2epJ4tRVVFZrxRFUUwruAaDQbS3t8Pn80mvTTAYxPTp02N6pMjeS0sObWtri0noDYVCaG9vt3IpLNFXTnW5XFi5cqX0fadPny7t7aKRVUulDJLyUMgmjnwQDQ52lmBqZH/Jtra2OrLENt5oh6IooqmpSdTX11t6v8jEyngJnPFGgGQJm9q0SkdHR8xoQeSSXkVR4p6vfmpFlhxqtszZCZG/W7/fL7xer9i1a5f088nOW38spyrEknWcdiGijGe3poSsWqcTtT1kQYc+8AAg1q5dG3U+TU1N4rbbbjN8fy33IjKwUBQlpkCW1etglrCpDxhkK4K0541u5JFBn1H1U6OVRsnSB19lZWVRP8+bN8/wOidTAI2cxeCDiAYEK0swzap1PvHEE2kJPADE9D+xEvSoqhq3NHi8G6XVvAUtYDAboTAaKUm20FkyksnX4chHZmHwQUQDhtFf9H6/37SpWLIPq0GH/ubV0dHhyDRP5BJSs7Lidm/2Vm68yfSNcbrSaCIJtWbnrQVr+uXGXCGTegw+iGhAS2WPFlnQES/w0B6RSz6tPCJLq8se2uiGUU6HlcZtslERK8FFppQrT2Tko7Gx0XD6SV9qPhM+42DB4IOIBiwnl81aCTzWGlQDlS3vNDsvo+mVyPb18QKJyBtlvOkYK8XABsqNVx8szZs3j0mmAxCX2hJRRjJaThqpq6sr7rJZRbG3+PUinL0b6fn27cMTkn4wO3fuREtLS9SSz8rKStPzeuSRR2K2hUIh5Ofno6enB1VVVaZLbt1uN/Ly8tDV1YVAIGC67NXKMmW3242ioqIB0Q+lvLwc3d3d8Pl8WL9+PV588UWEQiEoioLZs2dHXbfa2trwNYok+94Eg0EcPnw4bZ+DbEhDMGQLRz6IspP+L3WjtuhmSysTechGO8QX/9cn+2s5soBV5OiB2YiMVrTLyjSJ2YoV/UiG0fRJtq7mMLrGa9euFT6fL2oUST/aw5GP/sdpFyLKKGYBRX19fXg/p3M9ZEHHJ21tUecWeYOPDEBkeRIlJSXS96mqqoo5lp0kx3g9TWTLXrPxRmsUVGkBoZ1EWi3ng9KHwQcRZQwry1IbGhoczfW42yDwePjhh8XGjRtjin1ZHbWQnZ++1kUiuRZVVVW2RzKSWbGSqczqtmzatMnSNYosAsc+MOnF4IOIMoLVgEJRFLFly5aUjXYIQMyePTtqv7KysvB5Wqm1sXHjRtNRD7NrYFanw+iGa2UkYyAllVol6+BrZ1orG0eEBgoGH0SUEezUcEhHiXT9IzK3w6h6p9lUULwbm5U6HWadXAerhoaG8HWLHNWJN9pjdC2dKgFP5hh8EFFGSOWy2chHVwKBBwBRWloaPk+jzqlmgYfZkL6dsumypNfB/te6WfE5o9Ge3bt3G/6uOf2SelxqS0QZQetSKuu86hQBIE+y3cpi3Oeffx6BQABdXV3SzqmvvvqqdHltY2Mjuru7UV5ebnhsq0s/ZZ1cm5ubB8QS2VQyWipstoT417/+teHxEu2YTKnB4IOIUqq8vBxLlixJybGFZJsCa4EHcK4lvKymhqqqeOyxx2Je43K5MG/evLjBgVmdDr3IOhfxghoydsMNNxg+x5ofmYXBBxGlVCAQwLZt2xw9pjaWrmev9NhZCxYswJ49e2JGH8TZaemofVVVjSkOZkQ2omH22oFUFCxT5efno6ysTPqcUeBH/UMR+v+6+tnJkycxcuRI9Pb24sILL+zv0yGiBL366qt4/fXXAQBNTU2OHVf2f1g/BPB3BvsrioIHHngAP/zhDw2P6XK50N3dDQA4fPgwjh8/jtLS0pj9vF4v5s+fb+t8A4EADh8+jLy8PAYWadLZ2YkNGzbg+eefhxAiHPhxRCm17Ny/GXwQUVK0nIlJkyaFb67XXXcd2traHH2foQD+JNnuUlWEQiG4XC4Eg8GY5zs6OnD69GkUFxebHt/n86GoqAjA2c+Um5sblbOhBSj6AEL2+SkzMPBLLzv3b067EFHCPB4PcnNzUVxcjNzcXHg8Hrz66quOBx4C8sADQqCnpyecK6Efci8rK0N+fr40/yKSfkje7XajtrY2/BqjKRPZ56fMwamsDJbCVTcJ4VJbooHBaCnpokWLHF1GK1tC+8GOHTGFu+ItbTWq1yFbMitrzW718w/2JbI0eHGpLRGlnNFS0q9+9auOHP9myPM7FpeV4ervfCdqtCEQCMDr9ZoubS0pKZEmkLa3t0flAsg6xtbU1LCLKpGDzuvvEyCigUmbyoi8ASuKgj/96U/4yle+go8++ijhYxsloikA8PTT4Z9DoRCWLl0qXZkCRE+nGNXy6Ovri8rbMAsqIofvZZ+fKyqIrOHIBxElRL+UFACEEHjssceiAo9vfvObto4rCzxcMF5GGwqFDAOPyDwNo7obBw8ejMrbOHjwoKX6HHaX0hLROVztQkRJCQQCePXVV3HfffcldZy9ODvVopdI7Y7GxkZpIbDFixfj6YiRk3nz5uHFF1+MGb2ora1FTU0NgsFg3GWaXFFBdJad+zenXYgoYdp0xSeffGK6n6Io0tEJjdEzV33968CHH9o6tlEF0kAggGeeeSZqmz7wAM5OsUydOhXd3d2Wggq3282gg8gmBh9ElBCPxxOVmGnGbuARHu2IE3gAwKpVq9DY2Bg1SiELBmS5HKFQyDBvg0EFUeow54OIbNOvCEmEEyXSXS4Xli9fbqkvilHOR11dHfM2iNLMVvCxfv165OfnY8SIEbj00kvx7W9/G7/61a+i9hFC4B//8R8xbtw4fOlLX0JRURE+tPDXCxENHLJRBDtkQcduxA88FEWRBgpWikkZJYhWVVWxqRtRmtkKPvbv34/7778fBw4cwN69e/GXv/wFs2bNQl9fX3if+vp6bNq0CU8++SQ6OzsxZswY3HLLLTh16pTjJ09E/SNexVAjLhiPdtyu39flQllZWVSw0NLSklSgYNQ9lpUwidIrqdUuJ06cwKWXXor9+/fjxhtvhBAC48aNw4oVK7B69WoAwJkzZzB69GjU1dWhoqIi7jG52oVoYNiwYQOqq6st729auyPCunXrUFhYGM674GoSooEhbatdent7AQAXXXQRAODIkSM4duwYZs2aFd4nJycHM2fORFtbm6Xgg4jSx25TtEAgEO7bMmHCBMvvIws8/hrALyTbv/zlL4cbvAFcTZIqbIhH/Snh4EMIgcrKSlx//fW46qqrAADHjh0DAIwePTpq39GjR6Onp0d6nDNnzuDMmTPhn0+ePJnoKRGRDZGrVVRVRWVlJZYvXx5zI+rs7MRbb72FP/zhD1i3bl145YqixE8NvRbAu5LtZq9ctmwZzj//fOZepJD+d9/c3MzrTemVaAOZZcuWidzc3KgmSm+//bYAID755JOofe+55x5RUlIiPc6jjz4qbSbFxnJEqSNrigZAKIoS1WStrKzM0YZwwuJr4zVo8/v9MY3l0qm/3z8ZbIhHqZLyxnIPPvggXnnlFfh8vqi/ksaMGQPg3AiI5vjx4zGjIZo1a9agt7c3/PD7/YmcEhHZYLRaRQiBiooKBAIBdHZ2RlUDtUM2zTIU1pfRBoNBtLe3S5/r7zb2/f3+yWJDPMoEtoIPIQQeeOABvPjii9i3bx8mTpwY9fzEiRMxZswY7N27N7zt888/x/79+zFjxgzpMXNycnDhhRdGPYgotcxWq2g3/h//+Me2j7sFxqtZzki2A2enb2TnsmDBgpgbu6zjrBYspUN/v78TjOqdsCEepZOt4OP+++/Hjh078Oyzz2LEiBE4duwYjh07hj/96U8Azv6fyIoVK7Bu3Tq89NJLOHToEBYvXowLLrgAd999d0o+ABHZ53a7sWjRIulzqqqitLQUTU1Nto4pAMi6u5iNdmjLZ5ubm2NuiLIbe3//1d7f7+8ENsSjjGBnPgcG87M//vGPw/uEQiHx6KOPijFjxoicnBxx4403ig8++CAlc0ZElBijnA/ZNisPu7kdkydPFj6fLyrPYNeuXdJ9fT6f6XnbzVdIJl8jm/Il/H5/zO+AKBl27t8JJ5ymCoMPotTbt2+f9EZ/2223pTypdO3atdJzsnpjb21tFS6XK/x8ZIJsPK2treH3UFXV1mudeH+ibGbn/p1UkbFUYJExotQLBALIzc2NmUKI1302kmyvXwK4xuQ1u3fvxpw5cwyf93g8qKioiNvKPpHCY7LP7HK50N3dbXvKgYXPiGLZuX+zsRzRIOR2u1FXVxezPZnAQ4F54AEAc+fOjbs6JDKZ00gi5dCdzNdgOXai5DD4IBqkpkyZYvs12nyIntUltGarQ7SVJFoAJCKW/TqBqzyIMgeDD6JBym5zOFnQcQfkgYfR8llAPtoQCATg9XpTupKEqzyIMgeDD6JByu12o7a2Nm4AchmMRzt+LtuuKGhpaUFPTw+8Xm/c0QaPx4MJEyZg1apVMcdyemTCqKstEaUXgw+iQSgQCKC6uhqrV682za0QAI5ItptNswgh8Ic//AFdXV0oLCw0HW0IBAK45557pLkmqRqZYL4GUf9LqqstEQ08Ho8HS5YsiZtcKnt2JAArrR8feughAAg3LWtvb8d//Md/4Prrr0d+fn54v5qaGunrV65cicrKSgYIRFmKIx9Eg0ggEIgbeKyC8TSLFngoioInnngi7vuFQiEsWbIEBQUFqKysREFBQXi1SyAQwLPPPit93YUXXsjAgyiLMfggykKBQAA+ny9mpcjmzZtNAw8BYINku2ya5Y9//KOlhFUhhLQXSldXl+G5zJ49O+5xiWjgYvBBlGWMuq4GAgFs2rTJ8HVGox2dHR0xiaNCCDzxxBOm+SJGtBUsw4cPlz5/6623Rk3NEFH2Yc4H0QCmjSBMmjQJbrfbsOvq5MmTsXv3bmmwYDQOogAoKytDfn4+Tp8+bRpoqKqKl19+GcOHD8fBgwdRU1ODYDAIVVUhzrZxiNo3Ly8PbW1t0mMtXrzY6scnogGKwQfRAOXxeMKBhpbYefnll0trZUyfPl06xSELPD4FMO6Lf//kJz/BhAkTMG3aNKiqahiAhEIhDB8+HEVFRSgqKsKCBQvC5cf37NkTlWcihMCePXswYsQI6bF+//vfw+fzhQMqIso+7O1CNAAZ9Slpb29HQUGBpekQu5VKZ8yYgXfeeQfBYDDmObMeKWbnqg+KFOXsGQghwgEVa3EQDQzs7UKUpbRE0ra2NukIR3d3d1RdDZlES6S3tbXh5z//OXw+H+rr6y1XCjXqqdLX14eWlpZwLon2v1owYlaKnYgGNo58EGUYfR6HRj/Nos+lAM6OHNTV1SE3Nxd/+MMfcN9990U9L/uP/T4AWy2eW2NjI1asWBE+TyudXeN1k9WOc/z4cZSWlsa83ufzoaioyOIZElF/sXP/ZvBBlEFkeRzl5eXSG7g2UmBliuUSACck2602hNN0dHQktBLF4/GgoqICwWAwPFKin05xsuU9EaUfp12IBiCjlSraSIg+yAiFQnjkkUfiHlfAfuDR2NiIefPmRW3TVr4kwkpPFTZ+Ixo8OPJB1E/00ys+nw/FxcUx+/l8PuTl5SWUYCr7j3sczq5oMRI52tDZ2Ym3334b1113Xdpqb1idziGizMKRD6IMJysEJmtxr3V1lY0K1NbW4siRI1iyZElMgmkZjJNK9YGHqqrh99WPNuTn52PFihVpLfqVKY3fjKrEElHyOPJBlGZmuQ179uwxzY3QRgU6OzuxevXqqITTqqoq5OTk4HGDniuR0yxazQ7tPUpKSjjaEMEo94aIjDHhlCiDmU2vFBUVGU47aNM0w4cPN5xqsbKENrJYmKIoaGlp4Y01AhNfiRLDaReiDGY2vQLIpx0ip2lkgcfHMA48FEUJF+9yuVxRoyVCCNbS0DGqS3L48OF+OiOi7MPggyjN3G43Fi1aFLVt4cKFhn9Vy1bBRBIALpe8Thvx0OqBKIqCe+65J6Y2CG+s0eIFh0SUPAYfREmSJSaaJSsGAgE888wzUdt27NhhOPog+0tcYzjaIdtXCLS2toZHQTS8sUbjkl+i1GNjOaIkyBITAZgmK8Yb1te6vc6YMSP8HnpmnWjNBINBVFVVobGxMSqplTfWaOXl5UzCJUohJpwSJchq1VF9sqJRQmNNTQ2eMFipEkn2H+w/A3jUwjlr5wKAN1YichQTTolSyKy5WygUko5qtLS0hKdVZMP6d9xxR9zAYxjkgcfzXi+u3b0bK1euNH29qqrhUY54tTS0z9jZ2claF0TkPJFhent7BQDR29vb36dCFKO1tVWoqioACFVVhaIoWpNYAUAoihJ+Xv9QFEW0traGj6O9Vn8M2UMYPLTzaGhoEEIIUV9fL339zJkzhd/vt/0ZEfEe2rkTEcnYuX9z2oXIIqNpFhHRXVZRFPzN3/wNXnjhhZhVJdrzTz31FB544AFLDeEA+WjHlQB+rdtWX1+Pu+66CxMmTIh6b0VRcPToUUvTK7LPqGGtCyIyw2kXohQwau4WSQiBn/3sZ9LAQ3t+2bJllgKPb8J4NYs+8ACAmpoatLW1xby3EMLyUlqzlTVckktETuFqFyKLtPoP+pEPqyMYdiSymiUUCoULikUGIKqqWl5KK/uMGi7JJSKncOSDyCKj5m76glTJkgUeKqIDD1mjN5fLhcsuuyz2eBZmVrUEUwBRnzHy2FySS0RO4cgHkQ2y+g8XXXRRuBmcnqqqePzxx3H8+HH84Ac/MB0l2QNglmS7bLTj+uuvx6xZs7Bu3ToIIcLBwenTpw2nXSKX+nZ1dWHSpElwu93SWiXd3d04fPgwhg0bhr6+Pi7JJSJHMeGUSEJ/g7ay/+HDh3Hw4EHU1NSEC3gtXLgQzzzzTNSUiPbvyP/0jP4jVHX7aR5++GHU1taGj7VmzRrcfPPN0qZzkYmi+kCjrq4Oq1evZhM1Ikoau9oSJSHZdupaIDJs2LCYQEBVVezcuRNCCJSWlgKI34l25syZ2L9/v6X3VlUVixYtwo4dO6IqmJaXlxuu1pGNxmgddomIrLJz/+a0C1EEWRO3iooKlJSUWB4J0Ip4+Xw+6eqYQCCA66+/3nJSqdXAQzv+jh070N7eHjNdYrRaRz8Kw8RSIko1JpwSRXCynfq7774r3V5ZWYn8adNitv8I8XuzWBEMBtHX1xdVwTQQCODEiRPSbq11dXVsokZEacXggyiClXbqZh1rI/dZvXp1zPbzYDzNok3s3HfffZbOVd+d1uh8PR4PcnNzUVpaCiFE+PNpgUZ1dTW6u7vh8/nQ3d1ta4qJiCgRDD6IIsRrp67dyIuLi5Gbm4sNGzZIjyPr+yIA/I9k38gQoqysDJdccknc85w/fz6OHj0Kn8+HhoYGw/PVTyNp0yterzcq0IjX64WIyEnM+SDSMWqnLssHqa6uhhAC1dXV4dc3NDTgoYceijqmbLRjGoDOL/6tKAp+9KMfITc3Fz09Pabnd/fdd6Ouri6qQdyCBQukXWqN8jxGjRrFQIOI+g1XuxBJyJba+nw+FBcXx+yrqip6enrgdruxYcOGqEBkKs4FGJFkEyZWqqVqyaFWV+HIVrhwKS0RpQJ7u9CgZCUXw8r+RlMrsnwQ4OxIwuHDh2PyPASsBx7aceLR/lbQVuHE+6zxppGIiPoDgw8aEOIFFvqAwePxmB7PaH+jqZWGhga43W78zd/8jfR4Bw8ejJrikA0nno/EV7PIgh6rq3DKy8uZUEpEGYXTLpTx4hX9sju1YLZ/V1eX4dTKyy+/jNtvv11acVRRFDz++OMYuXYt7pd8hkSDDkVRsGvXLlx22WWmlUuJiPobp10oaxgV/dJGQAKBALxer7Q2R3t7u3S0xKyWh9nUilHgAZydDnnY4cADAOrr6zF//nzk5+dbnj6xO/1ERJR2IsP09vYKAKK3t7e/T4VSxO/3i3379gm/3x9333379gmcncWIevh8PtHa2ipUVZU+r6qqUBQl/O/W1tao99eei9xfO5+GhgbpMc0eQvKwewz9+dTX10uvnc/nM7x2kddE/7mJiFLJzv2b0y6UVnb7phhNkbS3t8dMQ0Q+HwqFYkqGa1MUgUAAEyZMiHpeURQcPXo0PJIwf/58/OxnP4v7eayWSLdq9+7dGD58eEJdZLmyhYj6E6ddKCPFm0KR0a/WUFUVK1euxJEjR6SBR2NjI5599tmY6ZHI5Myuri7DtvPaeb744otxP48s8HgNiQceJSUlmDNnTsLFvpwsDU9ElEosMkZpY3ZzNLvZlpSU4Nlnn8W+ffvQ0tKCDRs2QFVVaUO0efPmAYitmRFZclzL6zB6XnaeevE60dqlqipaW1uTOEL8z0VElCk48kFpY6Vvil5kX5Jt27ZFjZooihLTp0Sr+mmWnOl2u1FbWyt9LQAMHz7c8Hy0pAy9ZAIPl8uF5ubmpKdGWNODiAYK5nxQWnk8HlRUVCAYDIZvjkY5H7IcBj2v14tRo0ZJcyQCgUBUyXGtaum7776L1atXh/NO6urqUFVVFX5ddXW1tGeL7D+UWwH8i6VPLqcoCt555x3k5+cncZRo+s+dLrKqsEQ0eNi5fzP4oLSzenM0KmeusZNMGZnoanYcWTJqHoAuyTGTGe2I5PP5UFRU5NDR+ofdRGIiyj5MOKWMZrWDqlHNDcDelII+0VUvMimzra0tKvAQSG3gkQ05GYkkEhPR4MbggzKWLIehoaHBdplwKwmkBw8ejNkmGxIcgejAY8aMGejo6EBjYyO2bNli6Xw0qqpmRU4GV9kQkV1c7UIZzai9vR2yVSB6NTU1WLBgASZOnIhHADwm2Uc/2rFx40Zcc801GDt2LObNm4fHH3/c9rllA66yISK7mPNBadOfCYmRia76JbqaqqoqNEgSTQHnpln0sqUImJ1EYiLKTkw4pYxjJyHRSpAi2ycQCKCtrQ3A2ekQ2eqX9vZ23HnnndJjOr2E1qpsSDgF+m+VDRFlBjv3b067UMoZJSSWlJTE3KT0QYp+Gaxsn+bmZgDAkiVLwiMaiqKgrq4Oubm5AM4FI7JY+3OcbXevl47AQ1XVrJme0GqsEBHFZbdxzP79+8WcOXPE2LFjBQDx0ksvRT0fCoXEo48+KsaOHSuGDh0qZs6cKQ4dOmT5+Gwsl33MmsNF8vv90kZxDz/8cLgRnWwfVVUNG8xpD0VRRENDg1i+fHnchnCHkmwKZ+fR0NDQP78UIiKH2bl/217t0tfXh6uvvhpPPvmk9Pn6+nps2rQJTz75JDo7OzFmzBjccsstOHXqlN23ogEssq271cqmRqtS1q1bh+LiYuTm5mLz5s0x+4RCofjl0IVAdXU1Nm/efG6bZD8FwFXmH80xiqJgwYIFaXo3IqIMkkyUA93IRygUEmPGjBG1tbXhbX/+85/FyJEjxdatWy0dkyMfA5+srXtra6twuVwCgHC5XNJW70YjH5EPl8slFEWxPfIR+ZCNdog0jnZEPvSjP0REA1VKRz7MHDlyBMeOHcOsWbPC23JycjBz5sxwIqDemTNncPLkyagHDVxm+R3d3d2mNTrcbjfq6upMjx8MBrFq1aqo2h+1tbWorKyEosTP0pCNdtyF1OZ3mBVKy5Z8DyIiOxxNOD127BgAYPTo0VHbR48ejZ6eHulr1q9fj3/6p39y8jSoHxkVnGpvb8cll1wSd5ltVVUVhBCoqakxLIV+55134vLLL8exY8cQDAbD+6qqioqKCgwdOhTf//73o153KYD/K3m/VCeVulwutLe3o6+vDwcPHkRNTU3UclQmaBLRYJSS1S76v0CFEIZ/la5ZswaVlZXhn0+ePInx48en4rQoDWQFp1RVRWlpKYQQlvp+VFdX46677sLhw4djbth33HEHpk2bJn1dKBTCtm3bYla0GK0lNwo8jOqAWKW9XgswtKZxRUVFWLBgAZejEtGg52jwMWbMGABnR0DGjh0b3n78+PGY0RBNTk4OcnJynDwN6kdutxuLFi3C008/Hd4WGYiYLbPVH0frAbNgwQK0t7fj5Zdfxk9/+lPT97cSeFwC4PdxjqEFy4kEIYqiYNeuXSgsLIz5jFyOSkTkcG+XiRMnYsyYMdi7d2942+eff479+/djxowZTr4V9ZPIVSyy7Z2dnXjmmWdMj2G378fOnTtRWloaN/CItAjGq1nMAg+NNkrj9Xrx8MMPW35f4GyANWrUKAYZREQGbAcfp0+fxi9/+Uv88pe/BHA2yfSXv/wljh49CkVRsGLFCqxbtw4vvfQSDh06hMWLF+OCCy7A3Xff7fS5U5pogUVDQwNyc3PDy149Hg+As0W/tO0FBQVxl73aSbRsaGhAdXW1rREIAeAnku128zuCwSBGjRqF++67z3DacOHChTHPMZGUiCgOu0tpfD6fdMlgWVmZEOJckbExY8aInJwcceONN4oPPvjA8vG51DazRC6b1T9cLpfo6OiIu8w1cims0TJbGb/fH7OsNt5DtoR2xYoVMcXFrDxcLpfw+/2GRdK0ferr602XEWvH8Pv9qfgVERFlBDv376TqfKQCg4/MYaXuxqZNmwwDjsibsd/vFz6fz/QG3NHRITZu3Cg6OjqEEMaVUWWP/2OzdofL5RJlZWXhoEH2mD59uti1a1fcAEv7XLLPJ6t5QkSUjezcv9nbhWJoTdtOnDhhOoXicrlw/fXXS9upa8tLI1d1mOVAzJ49G6+//nr455KSEjz22GMxx5aRTcj8AcDFBvs3NTVhzpw5eO655xAMBg2P+84776C0tBSKoqCwsFBaq0abYpElktrpaZOs/uwYTERkWxqCIVs48tG/9H+pG017RE4vWKleamb27NmGowplZWWmUy92Rju0R2Njo2hoaLA9DWN2DWSs9rRJFkdXiCgT2Ll/K0IkUdAgBey05CVnBQIB5ObmxtToAM7+1a5VE506dWpMnYpE26l3dnYa1u3Q3l828mH0pbWSVPrQQw+hvr7e2gma8Hq9mD9/vuHzsuvpcrnQ3d3t2OhEOt6DiMgKO/dvTrtQmKw6aSgUgtfrxahRo0wDi0TrV+zevdv0eauBRwWAZt22v/3bv5Uuz3Ui8HC5XCgsLDTdx+12o7m5GRUVFSmrampUUfbw4cMMPogoYzH4oDBZdVLtJpuKG1kgELDVy+cCAH2S7dVVVci79FLgoYeittupC2KHqqqWg4jy8nKUlJSkrKqp0e+MS32JKJM5WmSMBjbtL/XIpm2p6D8SCASwcOFCTJgwIarFvRkBeeChANi4cSP+9V//1clTNKSqKg4cOGBaHl5Pq9SaigAuXb8zIiInMedjkDJbHZFo/oYVHo8H99xzT9z9Zs6ciTfffBNCCOk0yzgAnzp6ZnKRZda1G7udwCNdUvk7IyKyws79m8HHIOTxeMJLQK00enOKLDnSiKqq2LpwIZb8JLZWaao70WqqqqqwfPlyAOCNnYgoDgYfZKg/V0d4vV6UlpZa2jeZ1SxO6OjoCHejJSKi+Ozcv5nzMciYrY5IJY/HgwULFljaVxZ4KLAXeFx++eWG/VisiOzKTEREzmLwMchoqyMipXp1hFbpM94gWz2MAw+7fvOb39hqRqfX3t4e/rdRJ18iIkoMg49Bpj9WR8hGW/QEgGrdtveRvmkWI5EdeyM7+RIRUeKY8zFIpXN1RLxEU6dGO5yiqip6enoAgNVDiYgsYs7HIGR3aiCVtSdk7xU52qLRmp3opSvwKCkpiTknl8uF5uZmuN3ufsuPISLKdqxwmgUaGhqwevVqCCHSunTWjvLyckyePBnTp083rN2xTFHwx+98B64dO0y7zTpBVVW0trYCOLuMdtiwYTFdeFNRPZTdZ4mIwK62A52sO6vL5RJ+v7+/Ty3Grl27xHkGnWh9Pl/4nOvr65PuOGv2sNP5NdmOvfpjsfssEWUrdrUdJMxyKXw+H4qKitJ/UmYMlr7eW1GBrVu3ArBXiCwR3/zmN9HS0hIedbAyEuFEfgy7zxJRtmPOxyBhtIpEVdXMaywmCTyuwNn8jpaWlnCuipWVMckYOnRo+GZvdSWLE/kxzB8hIjqHwccAJqvZAQC1tbUZ8dd0IBDAOz/+sTTwUAD85ot/h0IhPPHEE/D5fBg+fHhSxcHimTt3bvjctBLz2jlUVFSkrJZHf9RXISLKVAw+BjD9KhJVVdHQ0IDqan3FDOeYraqJfM7j8cA9fjymf+97MfupkuBi69atKC4uRkFBQVLFwczk5uZi8eLFANI/EsHus0RE5zDnIwukq2ZHZEM6RVGwdOlSFBcXY8aMGdizZ0/UcyHJ1yrQ3Q13bi6qq6uxYcOGlJ2nzHe/+1386Ec/Oncu/ZSDwe6zRJSt2FiOHGeWCBrZdr4SwEbJ6xUAK1asQGFhISZOnIhp06al9HxdLhdqa2sxdepUwxu9x+NBRUUFgsFgeCQi05YoExENFAw+yHE+nw/FxcWm+8i+SJ8BuCglZ3TOzTffjNmzZ+O6667D2LFjbY0scCSCiMgZDD7IcZlcIr2jowP5+flpejciIpLhUltynJYwqV+x0QN54CFLKk2FW2+9NSbwcKoLbeRx2NmWiMg5DD7IsvLycvT09GDmzJkAzgYdE3T7HPne9xDw+9HS0hLTN8Vps2fPxmuvvRa1zakutJHHmTBhAiZMmMDOtkREDuG0C8Uwq/oZCASQO348ZJ1XAn5/VOXQn/zkJ/joo4/w8ccfo62tzbHz++pXv4r6+noMGzYs6hydWsESb4qJlUmJiGLZuX+zsRxFiVxOK2tS5zYIPBQAVZs3484778SGDRvg9XpTcn5TpkxBaWkpbr/99nAjvcrKSixfvty0doedQCFeldVEjklEROdw5IPCOjs7w11nNaqq4sCBA2fzKiR5HFcDeD9N5zdx4kTcfffdeOKJJ2KeUxQF9fX1WL16NUc+iIj6ARNOyTaPxyOtLhoKhTB/+nTDEunpCDzGjRsHRVFw5MgRaeABnK0xUlNTg9ra2qSriMoqx2q1TFiZlIgoeZx2oZg+J5EEcLbpvU66ltFOnjwZhw4dslRyPRgM4rzzzkN7ezv6+vqSqt1RXl6OkpKScA0QAKwHQkTkEAYfWcpKq3iNUY6D7HZ/z8KFKCwqguuLyqDJcrlcWL9+PWpqasJ5JmVlZfjDH/6Ab3/72+EVJ1ZVVlaGc1WKioqSOje32x117Rh0EBE5g9MuWcjuctNJkyZFdZL9DoyLhnl27MDo0aPDy22tUBQlfPzIf2tTGNXV1ejp6YHP50NtbS2efvppvPzyyygvL8e7774r7dyrP1akVHeoJSKi5DDhNMsksty0s7Mz3GvF6MuQ6DTLsmXLsHjx4qiy54B8CsPo3NevX481a9aEe7BE9mwBgOeffx6VlZUx7+3z+ZIe/SAiImu41HYQs7LcNHJKRutGC6SmRPqWLVuwdevW2CW7kkDI6Nzz8/PR3d1tmHMxf/58VFVVxQQtWnBCRESZhdMuWWbSpEkx0xSRN2J95c4lS5bgrVAopb1ZrEyDBAIBnDhxImYaRTt3t9uNoqIiadCiX53CFSlERJmNwUeWMbsRBwIBLFmyJDxCIIRASAjM0B3jh3B+NYs2+iKjBUSlpaUAkNCy1vLycnR3d8Pn86G7uztqlIWIiDILcz6yVGSreODslMZHH32EZcuWhfdJZyfayLyTyGkfANI8j+eeew6FhYUcvSAiGiCY85Gl7Cyf1ZaJNjQ0YPXq1RBChEcUnE4qjSdyBENfvr2yslKa5zFq1CgGHkREWYojHwNEvJ4rssBkw4YNqK6ujjqO7Jd9I4C3vvj3o48+CiEE/vmf/zmh81RVFaFQKGZVijbioR/l0PJTki2JTkRE/cvO/ZvBxwAQb/msLDApKSmJes3FAH4nObZ+tKOjowNjx4417W1iROsDY1Rd1OfzSQuGVVVVobGxMbyUdtu2bczZICIaYNjbJcuYLZ/Vl0bXVpa0tbWdSyyFtcADAPr6+sJJq0a0KqSRq2oURUFzczPGjh1rWArdaCXO8uXLmSxKRDSIMPgYAMyWzxoFJmb5HSMgDzxUVQ0nqJaUlBiez1NPPYXt27ejp6cHXq8XXq8XR48eBQDTyqpmK3HMltISEVF2YcJpBjFKKN2zZ09Mm/tt27YBAE6cOBFzHFVVcfVvf2t7NYsQAlu3bsVFF12EvXv3Gu43ZcoUAGeDifnz54fPXTYCU1JSEvVZ9A3bGGwQEQ0+DD4yhFHeRltbG5YsWRIzlfHZZ58Z5mUEQyFg5cqY7fFWswghDFvWR+rr64vZZqWyqkbfsI2IiAYXJpz2s0AggLa2NixYsCBmdEMIYZg/oa0q0Ut17Q5VVdHT0xMTPCTSU4aIiLIHE04HiMjKnvogIxQK2Qo8vEhP0bDm5maWOCcioqRw5KOfyEYKrNC6vNbU1EStZtH79VVX4cpDhxw403Oamppw7733mu4TWVmVgQcR0eDBCqcDgCxHwoyqqnjkkUcwZ84c5Ofn46KLLsLSpUvP5nfoKACUDz908GzPLqWdM2dO3P2Yz0FERPFw2qWfyJbP6mnLZRVFCVcdLSgogMfjQfk99xgGHgAMp2wSoaoqWlpaGFQQEZEjGHz0E/3yWRlFUbBu3bpw8AGczQUpv+eemH1vg/P5Hbt374bP50NPTw8LfxERkWOY89EPEs33GAbgtGS7E0FHQUEBDhw4EP65rKwM27dvd+DIREQ0GDDnox/pC4W9+uqreP3113HrrbeGcybs5nsAqe1EqwUanZ2dePvtt3HdddchPz/fgSMTERHF4shHgmTVSPWFwsaOHYvf/va34dfk5eWhubkZw4cPR0FBgeUARPYLugTA75P8DIqi4JVXXrGUSEpERGSGXW1TzEoXWTOKomDWrFnYu3cvQqFQVE5HpEIAbbLXJ/8RwufNXA4iInICp10siBy5AIC2trO3+RkzZsDtdhv2Wens7MQ9EQmfoVAIS5YsweOPP259JEMI7NmzB4qioKqqCqdOnQr3agnvY/Bap5JKd+7cGe7LQkRElE6DcuQjcuRCP+qgKAq+853v4Jlnnoka2SgvL4fH44kKPJzgcrliqpmmulIpy54TEZHTMqK8+pYtWzBx4kQMHToUU6ZMwVtvvZWqt7JF331VH3sJIfD000/HdGfVj3g4JRgMhs/hIaSnNwvLnhMRUX9KSfCxa9curFixAmvXrsV7772HG264Ad/61rdw9OjRVLydZYFAAF6v1/ZKk2AwiFdffTVFZ3WWAFCn2/YanA08FEXBgQMHmOdBRET9KiXTLtOnT8e1116Lpqam8LavfvWr+Pa3v43169ebvjZV0y6RUy2JMEoKdUI6GsJpfD4fioqKUnR0IiIarPp12uXzzz/Hu+++i1mzZkVtnzVrVjipM930Uy2JSEXg0YbUBR6y0u0ulwt5eXkOHJ2IiChxjgcfv/vd7xAMBjF69Oio7aNHj8axY8di9j9z5gxOnjwZ9XBaIkW9Uk3g7FLaSN+CM4GHNr1SVVUVDkLY4p6IiDJFypbaak3RNEKImG0AsH79evzTP/1Tqk4DwLkmbpkQgJwH4H8k252aZtFW5+Tn5yM/Px/Lly9ni3siIsoojo98XHLJJXC5XDGjHMePH48ZDQGANWvWoLe3N/zw+/1OnxLcbjeam5vhcrni7isLkBRFCY8gqKoq3ccKAWcCjyuuuCLqZ1VV0dDQIG0C53a7UVRUxMCDiIgyhuMjH0OGDMGUKVOwd+9e3HHHHeHte/fuxdy5c2P2z8nJQU5OjtOnEaO8vBwlJSU4fPgwtmzZgueffz78nKIoqK+vx9SpU5GXl4dPP/0Ur732GnJycpCXl4fCwrMTJNoIgvZv/XHMyHI7xgD4vwb76xNchwwZgiVLluCb3/wm5syZg0AggPb2dgBAYWEhgwsiIhowUrLaZdeuXVi0aBG2bt2KwsJCNDc3o6WlBR9++CFyc3NNX5uu8uqdnZ147bXXMGbMGMyZMyfhm7d2nD//+c8YOnQoLrzwQrz//vu48sorcfPNNyPY1YWCv/3bmNd97atfxbx581BcXBwOaF599VUcO3YMs2fPRn5+PrZv346XX34Zc+fOxeLFi5P5uERERCmVEb1dtmzZgvr6enz66ae46qqr0NjYiBtvvDHu6wZCbxfLjKZnMquoLBERUdIyIvhIVNYEH7LAIxgEJEtgiYiIBrqMKK8+aL3yijzwEIKBBxERERh8OEtRAH1S7fe+x2kWIiKiCCmr8zHoGI12EBERURSOfCRryxYGHkRERDZw5CMZsqBj/37AwqoeIiKiwYrBRyKMkkc52kFERBQXp13smjePgQcREVESOPJhh2ya5Xe/Ay6+OP3nQkRENEAx+LDi9GlgxIjY7RztICIiso3TLvHMnBkbeMyaxcCDiIgoQRz5MMMS6URERI7jXVTmN79hiXQiIqIU4Z1U7+KLgSuuiN62cyenWYiIiBzCaZdIrFRKRESUchz5AID332fgQURElCYMPkaNAq6+Onrbhx8y8CAiIkqRwTvtwhLpRERE/WJwjnz4/bGBx4MPMvAgIiJKg8E38tHcDFRURG87eVJewZSIiIgcN7iCj7lzgVdeid7G0Q4iIqK0GjzTLkJEBx4ff8zAg4iIqB8MnuBDUYB//Vfghz8EQiHg8sv7+4yIiIgGpcE17VJScvZBRERE/WbwjHwQERFRRmDwQURERGnF4IOIiIjSisEHERERpRWDDyIiIkorBh9ERESUVgw+iIiIKK0YfBAREVFaMfggIiKitGLwQURERGnF4IOIiIjSisEHERERpRWDDyIiIkqrjOtqK4QAAJw8ebKfz4SIiIis0u7b2n3cTMYFH6dOnQIAjB8/vp/PhIiIiOw6deoURo4cabqPIqyEKGkUCoXwySefYMSIEVAUxdFjnzx5EuPHj4ff78eFF17o6LEHKl6TaLwesXhNYvGaxOI1iTYYr4cQAqdOncK4ceOgquZZHRk38qGqKtxud0rf48ILLxw0XwareE2i8XrE4jWJxWsSi9ck2mC7HvFGPDRMOCUiIqK0YvBBREREaTWogo+cnBw8+uijyMnJ6e9TyRi8JtF4PWLxmsTiNYnFaxKN18NcxiWcEhERUXYbVCMfRERE1P8YfBAREVFaMfggIiKitGLwQURERGmVdcHHli1bMHHiRAwdOhRTpkzBW2+9Zbr//v37MWXKFAwdOhSXX345tm7dmqYzTR871+QXv/gFFEWJeXz00UdpPOPUefPNN3Hbbbdh3LhxUBQFP//5z+O+Jtu/I3avSbZ/R9avX4/8/HyMGDECl156Kb797W/jV7/6VdzXZfP3JJFrks3fk6amJkyePDlcQKywsBD/8i//YvqabP5+JCKrgo9du3ZhxYoVWLt2Ld577z3ccMMN+Na3voWjR49K9z9y5AhuvfVW3HDDDXjvvffw8MMP4+/+7u/wwgsvpPnMU8fuNdH86le/wqeffhp+TJo0KU1nnFp9fX24+uqr8eSTT1rafzB8R+xeE022fkf279+P+++/HwcOHMDevXvxl7/8BbNmzUJfX5/ha7L9e5LINdFk4/fE7XajtrYWBw8exMGDB1FcXIy5c+fiww8/lO6f7d+PhIgsMm3aNHHvvfdGbfvKV74iampqpPs/9NBD4itf+UrUtoqKClFQUJCyc0w3u9fE5/MJAOKzzz5Lw9n1LwDipZdeMt1nMHxHIlm5JoPpOyKEEMePHxcAxP79+w33GWzfEyvXZLB9T7785S+L1tZW6XOD7fthRdaMfHz++ed49913MWvWrKjts2bNQltbm/Q17e3tMfuXlJTg4MGD+J//+Z+UnWu6JHJNNNdccw3Gjh2Lm266CT6fL5WnmdGy/TuSjMHyHent7QUAXHTRRYb7DLbviZVrosn270kwGMTOnTvR19eHwsJC6T6D7fthRdYEH7/73e8QDAYxevToqO2jR4/GsWPHpK85duyYdP+//OUv+N3vfpeyc02XRK7J2LFj0dzcjBdeeAEvvvgirrzyStx00014880303HKGSfbvyOJGEzfESEEKisrcf311+Oqq64y3G8wfU+sXpNs/5588MEHGD58OHJycnDvvffipZdewte+9jXpvoPp+2FVxnW1TZaiKFE/CyFitsXbX7Z9ILNzTa688kpceeWV4Z8LCwvh9/uxYcMG3HjjjSk9z0w1GL4jdgym78gDDzyA999/H//xH/8Rd9/B8j2xek2y/Xty5ZVX4pe//CX+3//7f3jhhRdQVlaG/fv3GwYgg+X7YVXWjHxccsklcLlcMX/RHz9+PCbi1IwZM0a6/3nnnYeLL744ZeeaLolcE5mCggJ0dXU5fXoDQrZ/R5ySjd+RBx98EK+88gp8Ph/cbrfpvoPle2Lnmshk0/dkyJAhyMvLw9SpU7F+/XpcffXV2Lx5s3TfwfL9sCNrgo8hQ4ZgypQp2Lt3b9T2vXv3YsaMGdLXFBYWxuz/xhtvYOrUqTj//PNTdq7pksg1kXnvvfcwduxYp09vQMj274hTsuk7IoTAAw88gBdffBH79u3DxIkT474m278niVwTmWz6nugJIXDmzBnpc9n+/UhIPyW6psTOnTvF+eefLzwej/iv//ovsWLFCjFs2DDR3d0thBCipqZGLFq0KLz/b37zG3HBBReIlStXiv/6r/8SHo9HnH/++eJnP/tZf30Ex9m9Jo2NjeKll14Sv/71r8WhQ4dETU2NACBeeOGF/voIjjp16pR47733xHvvvScAiE2bNon33ntP9PT0CCEG53fE7jXJ9u/IfffdJ0aOHCl+8YtfiE8//TT8+OMf/xjeZ7B9TxK5Jtn8PVmzZo148803xZEjR8T7778vHn74YaGqqnjjjTeEEIPv+5GIrAo+hBDiqaeeErm5uWLIkCHi2muvjVoKVlZWJmbOnBm1/y9+8QtxzTXXiCFDhojLLrtMNDU1pfmMU8/ONamrqxNXXHGFGDp0qPjyl78srr/+evHaa6/1w1mnhrb8T/8oKysTQgzO74jda5Lt3xHZtQAgfvzjH4f3GWzfk0SuSTZ/T773ve+F/z911KhR4qabbgoHHkIMvu9HIhQhvsh6ISIiIkqDrMn5ICIiooGBwQcRERGlFYMPIiIiSisGH0RERJRWDD6IiIgorRh8EBERUVox+CAiIqK0YvBBREREacXgg4iIiNKKwQcRERGlFYMPIiIiSisGH0RERJRW/z94COqMFjHpUwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(x,y,'k.')\n",
    "plt.plot(x,result.slope*x+result.intercept,'r-')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Method 2: statsmodels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ordinary least squares fit using [statsmodels](https://www.statsmodels.org/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "smres = smf.ols('NITRAT ~ PHSPHT',df07).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>         <td>NITRAT</td>      <th>  R-squared:         </th> <td>   0.972</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.972</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>7.757e+04</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Fri, 03 Mar 2023</td> <th>  Prob (F-statistic):</th>  <td>  0.00</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>09:43:25</td>     <th>  Log-Likelihood:    </th> <td> -4993.9</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>  2210</td>      <th>  AIC:               </th> <td>   9992.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>  2208</td>      <th>  BIC:               </th> <td>1.000e+04</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     1</td>      <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>   -3.9326</td> <td>    0.103</td> <td>  -38.336</td> <td> 0.000</td> <td>   -4.134</td> <td>   -3.731</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>PHSPHT</th>    <td>   14.7400</td> <td>    0.053</td> <td>  278.514</td> <td> 0.000</td> <td>   14.636</td> <td>   14.844</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>874.728</td> <th>  Durbin-Watson:     </th> <td>   0.269</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.000</td>  <th>  Jarque-Bera (JB):  </th> <td>5147.310</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td>-1.766</td>  <th>  Prob(JB):          </th> <td>    0.00</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 9.589</td>  <th>  Cond. No.          </th> <td>    4.90</td>\n",
       "</tr>\n",
       "</table><br/><br/>Notes:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                 NITRAT   R-squared:                       0.972\n",
       "Model:                            OLS   Adj. R-squared:                  0.972\n",
       "Method:                 Least Squares   F-statistic:                 7.757e+04\n",
       "Date:                Fri, 03 Mar 2023   Prob (F-statistic):               0.00\n",
       "Time:                        09:43:25   Log-Likelihood:                -4993.9\n",
       "No. Observations:                2210   AIC:                             9992.\n",
       "Df Residuals:                    2208   BIC:                         1.000e+04\n",
       "Df Model:                           1                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept     -3.9326      0.103    -38.336      0.000      -4.134      -3.731\n",
       "PHSPHT        14.7400      0.053    278.514      0.000      14.636      14.844\n",
       "==============================================================================\n",
       "Omnibus:                      874.728   Durbin-Watson:                   0.269\n",
       "Prob(Omnibus):                  0.000   Jarque-Bera (JB):             5147.310\n",
       "Skew:                          -1.766   Prob(JB):                         0.00\n",
       "Kurtosis:                       9.589   Cond. No.                         4.90\n",
       "==============================================================================\n",
       "\n",
       "Notes:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "smres.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Method 3: pingouin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>names</th>\n",
       "      <th>coef</th>\n",
       "      <th>se</th>\n",
       "      <th>T</th>\n",
       "      <th>pval</th>\n",
       "      <th>r2</th>\n",
       "      <th>adj_r2</th>\n",
       "      <th>CI[2.5%]</th>\n",
       "      <th>CI[97.5%]</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Intercept</td>\n",
       "      <td>-3.932572</td>\n",
       "      <td>0.102582</td>\n",
       "      <td>-38.335853</td>\n",
       "      <td>6.540950e-247</td>\n",
       "      <td>0.972323</td>\n",
       "      <td>0.972311</td>\n",
       "      <td>-4.133740</td>\n",
       "      <td>-3.731405</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>PHSPHT</td>\n",
       "      <td>14.740035</td>\n",
       "      <td>0.052924</td>\n",
       "      <td>278.514375</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>0.972323</td>\n",
       "      <td>0.972311</td>\n",
       "      <td>14.636249</td>\n",
       "      <td>14.843820</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       names       coef        se           T           pval        r2  \\\n",
       "0  Intercept  -3.932572  0.102582  -38.335853  6.540950e-247  0.972323   \n",
       "1     PHSPHT  14.740035  0.052924  278.514375   0.000000e+00  0.972323   \n",
       "\n",
       "     adj_r2   CI[2.5%]  CI[97.5%]  \n",
       "0  0.972311  -4.133740  -3.731405  \n",
       "1  0.972311  14.636249  14.843820  "
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pg.linear_regression(x[ii],y[ii])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Method 4: Regression coefficients using numpy's polyfit function\n",
    "We can also use the `polyfit` function from numpy to calculate the coefficients of the line of best fit (minimizing the sum of square errors): "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[14.74003452 -3.93257206]\n"
     ]
    }
   ],
   "source": [
    "coefficients = np.polyfit(x[ii],y[ii],1)\n",
    "print(coefficients)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Method 5: Matrix form\n",
    "\n",
    "In the next section, we'll consider solving for the coefficients of the linear fit using matrices. But first, let's do a quick review of matrix multiplication:\n",
    "\n",
    "#### Review: Matrix Multiplication\n",
    "Suppose we have the following matrices:\n",
    "\n",
    "$$ \\textbf{A}= \\begin{bmatrix} 1 & 2 & 3\\\\ 4 & 5 & 6 \\end{bmatrix} $$\n",
    "\n",
    "$$ \\textbf{B} = \\begin{bmatrix} 7 & 8\\\\ 9 & 10 \\\\ 11 & 12 \\end{bmatrix} $$\n",
    "\n",
    "What will be the matrix product $\\textbf{AB}$?\n",
    "\n",
    "To define matrices in Python, we define 2-d arrays as lists of lists wrapped in numpy's ```array``` function, for example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define matrix A\n",
    "A = np.array([[1, 2, 3], \n",
    "              [4, 5, 6]])\n",
    "\n",
    "# define matrix B \n",
    "B = np.array([[7, 8], \n",
    "              [9, 10], \n",
    "              [11, 12]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check the dimensions of these array's using numpy's ```shape``` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "shape of A:  (2, 3)\n",
      "shape of B:  (3, 2)\n"
     ]
    }
   ],
   "source": [
    "print('shape of A: ', np.shape(A))\n",
    "print('shape of B: ', np.shape(B))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can multiply two arrays using numpy's ```dot``` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 58,  64],\n",
       "       [139, 154]])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.dot(A,B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is important to remember that matrix multiplication is not *commutative*, meaning $\\textbf{AB}$ is generally not the same as $\\textbf{BA}$. In this example, $\\textbf{BA}$ gives us a different size matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 39,  54,  69],\n",
       "       [ 49,  68,  87],\n",
       "       [ 59,  82, 105]])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.dot(B,A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Review: Matrix Transpose\n",
    "\n",
    "The *transpose* of a matrix $\\textbf{A}^T$ has the same values as $\\textbf{A}$, but the rows are converted to columns. One way to do this is with the `np.transpose` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1 2 3]\n",
      " [4 5 6]]\n"
     ]
    }
   ],
   "source": [
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1 4]\n",
      " [2 5]\n",
      " [3 6]]\n"
     ]
    }
   ],
   "source": [
    "print(np.transpose(A))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another way is to use the `.T` method on a Numpy array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1 4]\n",
      " [2 5]\n",
      " [3 6]]\n"
     ]
    }
   ],
   "source": [
    "print(A.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the product $\\textbf{A}^T\\textbf{A}$ is a *square* matrix, which has the same number of rows and columns. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[17, 22, 27],\n",
       "       [22, 29, 36],\n",
       "       [27, 36, 45]])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.dot(A.T, A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Review: Matrix Inverse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The concept of a matrix inverse is similar to the matrix of a single number. If we have a single value $b$, its inverse can be represented as $b^{-1}$. A value times its inverse is equal to 1.\n",
    "\n",
    "$$b^{-1}b = 1$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's say we have a *square* matrix $\\textbf{B}$ where the number of rows and columns are equal. The inverse $\\textbf{B}^{-1}$ is the matrix that gives \n",
    "\n",
    "$$ \\textbf{B}^{-1}\\textbf{B} = \\textbf{I} $$\n",
    "\n",
    "where the *identity matrix* $\\textbf{I}$ is a matrix that has all 0 values, except for 1 values along the diagonal from the upper left to the lower right. For example, a 3x3 identity matrix would be\n",
    "\n",
    "$$ \\textbf{I} = \\begin{bmatrix} 1 & 0 & 0 \\\\ 0 & 1 & 0 \\\\ 0 & 0  & 1  \\end{bmatrix} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is called an *identity matrix* because $\\textbf{B}\\textbf{I} = \\textbf{I}\\textbf{B} = \\textbf{B}$ for square matrices. This is analagous to $1b = b$ for single values.\n",
    "\n",
    "In a linear algebra class, you might calculate $\\textbf{B}^{-1}$ by hand, but in this class we will rely in Numpy to do it for us. Let's set up a $3 \\times 3$ $\\textbf{B}$ matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1 2 1]\n",
      " [3 2 1]\n",
      " [1 2 2]]\n"
     ]
    }
   ],
   "source": [
    "B = np.array([[1, 2, 1],\n",
    "              [3, 2, 1],\n",
    "              [1, 2, 2]])\n",
    "print(B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The inverse $\\textbf{B}^{-1}$ is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.5 ,  0.5 ,  0.  ],\n",
       "       [ 1.25, -0.25, -0.5 ],\n",
       "       [-1.  ,  0.  ,  1.  ]])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.inv(B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The product $\\textbf{B}^{-1}\\textbf{B}$ can be calculated as"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.00000000e+00  3.33066907e-16  1.66533454e-16]\n",
      " [-5.55111512e-17  1.00000000e+00 -2.22044605e-16]\n",
      " [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]\n"
     ]
    }
   ],
   "source": [
    "BinvB = np.dot(np.linalg.inv(B), B)\n",
    "print(BinvB)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we round these values, we can see more clearly that this is nearly identical to the identity matrix $\\textbf{I}$, with some very small round-off error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.  0.  0.]\n",
      " [-0.  1. -0.]\n",
      " [ 0.  0.  1.]]\n"
     ]
    }
   ],
   "source": [
    "print(np.round(BinvB))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For reference, an identity matrix can be created with the `np.eye` function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1. 0. 0.]\n",
      " [0. 1. 0.]\n",
      " [0. 0. 1.]]\n"
     ]
    }
   ],
   "source": [
    "print(np.eye(3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Formulating linear regression in matrix form\n",
    "\n",
    "We can formulate regression in terms of matrices as\n",
    "\n",
    "$$\\begin{bmatrix} y_1\\\\ y_2\\\\ \\vdots \\\\ y_N \\end{bmatrix} = \\begin{bmatrix} 1 & x_1\\\\ 1 & x_2\\\\ \\vdots & \\vdots\\\\ 1 & x_N \\end{bmatrix} \\begin{bmatrix} c_0\\\\ c_1\\end{bmatrix} + \\begin{bmatrix} \\varepsilon_1\\\\ \\varepsilon_2\\\\ \\vdots \\\\ \\varepsilon_N \\end{bmatrix}$$\n",
    "\n",
    "or\n",
    "\n",
    "$$\\vec{y} = \\textbf{X}\\vec{c} + \\vec{\\varepsilon}$$\n",
    "\n",
    "Here, $\\vec{\\varepsilon}$ represents the vector of errors, or differences between the model and data values.\n",
    "\n",
    "$$ \\vec{\\varepsilon} = \\hat{y} - \\vec{y} $$\n",
    "\n",
    "To solve for the parameters c using matrix multiplication, we first need to fomulate the $\\vec{y}$ and $\\textbf{X}$ matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[39.14]\n",
      " [40.36]\n",
      " [42.36]\n",
      " ...\n",
      " [26.46]\n",
      " [25.29]\n",
      " [23.98]]\n",
      "shape of y:  (2210, 1)\n"
     ]
    }
   ],
   "source": [
    "# check to see that the y_subset is only 1-d (and won't work for matrix multiplication)\n",
    "# print(np.shape(y_subset))\n",
    "\n",
    "# define an (N, 1 matrix of the y values)\n",
    "y_matrix = np.reshape(y[ii].ravel(),(len(y[ii]),1))\n",
    "\n",
    "# print the matrix\n",
    "print(y_matrix)\n",
    "\n",
    "# print the shape of the matrix\n",
    "print('shape of y: ', np.shape(y_matrix))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1.   2.73]\n",
      " [1.   2.83]\n",
      " [1.   2.94]\n",
      " ...\n",
      " [1.   2.35]\n",
      " [1.   2.26]\n",
      " [1.   2.1 ]]\n",
      "(2210, 2)\n"
     ]
    }
   ],
   "source": [
    "# define a matrix X with a column of ones and a column of the x values\n",
    "X = np.column_stack([np.ones((len(x[ii]),1)),\n",
    "                     np.ravel(x[ii])])\n",
    "\n",
    "# print the X matrix\n",
    "print(X)\n",
    "\n",
    "# print the shape of the X matrix\n",
    "print(np.shape(X))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have our matrices set up, let's return to our model equation in matrix form.\n",
    "\n",
    "$$\\hat{y} = \\textbf{X}\\vec{c}$$\n",
    "\n",
    "We can multiply each side of the equation by the transpose\n",
    "\n",
    "$$\\textbf{X}^T\\hat{y} = \\textbf{X}^T\\textbf{X}\\vec{c}$$\n",
    "\n",
    "then multiply each side by the inverse of $\\textbf{X}^T\\textbf{X}$\n",
    "\n",
    "$$(\\textbf{X}^T\\textbf{X})^{-1}\\textbf{X}^T\\hat{y} = (\\textbf{X}^T\\textbf{X})^{-1}\\textbf{X}^T\\textbf{X}\\vec{c}$$\n",
    "\n",
    "which reduces to\n",
    "\n",
    "$$(\\textbf{X}^T\\textbf{X})^{-1}\\textbf{X}^T\\hat{y} = \\textbf{I}\\vec{c}.$$\n",
    "\n",
    "$$(\\textbf{X}^T\\textbf{X})^{-1}\\textbf{X}^T\\hat{y} = \\vec{c}$$\n",
    "\n",
    "giving us an expression for the vector of coefficents $\\vec{c}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here $\\hat{y}$ represents the vector of *model* values. With the vector of *data* values $\\vec{y}$ we can solve for the coefficients that minimize the sum of square errors using the same equation:\n",
    "\n",
    "$$\\vec{c} = (\\textbf{X}^T\\textbf{X})^{-1}\\textbf{X}^T\\vec{y}$$\n",
    "\n",
    "Using numpy, we can define the matrix components of the above equation and then run the calculation to find the coefficients:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-3.93257206]\n",
      " [14.74003452]]\n"
     ]
    }
   ],
   "source": [
    "# calculate the transpose of the matrix X\n",
    "XT = np.transpose(X)\n",
    "\n",
    "# calculate the product of XT and X\n",
    "XTX = np.dot(XT,X)\n",
    "\n",
    "# calculate the inverse of XTX\n",
    "XTX_inverse = np.linalg.inv(XTX)\n",
    "\n",
    "# then calculate the products of XTX, XT and using two dot products\n",
    "c = np.dot(XTX_inverse,np.dot(XT,y_matrix))\n",
    "\n",
    "# print the coefficients\n",
    "print(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a sanity check, we can double check that the coefficients are the same as those from numpy's ```polyfit``` function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[14.74003452 -3.93257206]\n"
     ]
    }
   ],
   "source": [
    "coefficients = np.polyfit(x[ii],y[ii],1)\n",
    "print(coefficients)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.8"
  },
  "vscode": {
   "interpreter": {
    "hash": "0ef88d3abb6b62f34a20525ce337090c4512fe8aecf32c74604482b944e1c3bd"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}