quantA3-rep.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d4462900-2034-46fb-9996-0f1af3459203",
   "metadata": {},
   "source": [
    "# Quantitiative Methods Reproducible Analysis for Assessment 3\n",
    "Student ID: 21197822\n",
    "\n",
    "References for code are linked for come from class practicals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "94967fb8-04b4-47ae-aa3d-139a89914ab3",
   "metadata": {},
   "outputs": [],
   "source": [
    "#load packages\n",
    "\n",
    "import os\n",
    "import re\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from matplotlib import pyplot as plt\n",
    "import seaborn as sns\n",
    "from sklearn.linear_model import LinearRegression\n",
    "import statsmodels.api as sm\n",
    "import sklearn\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b49736d-a527-4f5a-aeb2-cb16dc39d58b",
   "metadata": {},
   "source": [
    "## Loading and Cleaning Data"
   ]
  },
  {
   "cell_type": "raw",
   "id": "7417dd1e-9d56-49ea-aff5-c0c326515493",
   "metadata": {},
   "source": [
    "# census data was not used\n",
    "\n",
    "# https://data.london.gov.uk/download/ward-profiles-and-atlas/5a367ec8-3bdb-4216-9f75-8368a00fe6a2/ward-atlas-data.xls\n",
    "wardatlas = pd.read_excel(\n",
    "    'https://data.london.gov.uk/download/ward-profiles-and-atlas/5a367ec8-3bdb-4216-9f75-8368a00fe6a2/ward-atlas-data.xls',\n",
    "    #sheet_name='Summary - Mean Scores',\n",
    "    header=[0,1,2])            # Where are the column names... there's three of them!\n"
   ]
  },
  {
   "cell_type": "raw",
   "id": "c1fcef81-2066-4489-b959-0a0bb3279284",
   "metadata": {},
   "source": [
    "#wardatlas.columns.values"
   ]
  },
  {
   "cell_type": "raw",
   "id": "d6ad6c38-c907-4d43-93e7-00d1e3f4b086",
   "metadata": {},
   "source": [
    "atlas0 = wardatlas[[('Unnamed: 1_level_0', 'Unnamed: 1_level_1', 'New Code'),\n",
    "                    ('Unnamed: 2_level_0', 'Unnamed: 2_level_1', 'Borough'),\n",
    "                    ('Unnamed: 3_level_0', 'Unnamed: 3_level_1', 'Names'),\n",
    "                    ('Population and Age', 'Aged 16-64', 2013),\n",
    "                    ('Population and Age', 'Aged 65+', 2013)]]\n",
    "atlas = atlas0.droplevel(0,axis=1)\n",
    "atlas = atlas.droplevel(0,axis=1)\n",
    "atlas.columns = ['New Code','Borough','Names','Aged 16-64','Aged 65+']\n",
    "atlas['Aged 16+'] = atlas['Aged 16-64'] + atlas['Aged 65+']\n",
    "borosum = atlas.groupby('Borough').agg({'Aged 16+' :'sum'})\n",
    "atlas = pd.merge(atlas, borosum, left_on = 'Borough', right_on = 'Borough', how = 'left')\n",
    "atlas.rename(columns = {'Aged 16+_x' : 'Aged 16+', 'Aged 16+_y':'Boro 16+'}, inplace=True)\n",
    "atlas['% 16+'] = atlas['Aged 16+'] / atlas['Boro 16+']\n",
    "atlas = atlas[['New Code','% 16+']]"
   ]
  },
  {
   "cell_type": "raw",
   "id": "fe6e884c-e354-413f-9a48-3934b20d9335",
   "metadata": {},
   "source": [
    "atlas0.isnull().sum(axis=0).sort_values(ascending=False)\n",
    "atlas0.tail()\n",
    "atlas.tail()\n",
    "atlas.dropna(inplace=True)\n",
    "atlas.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d37ffbb7-a862-418e-a6aa-b6120af966cc",
   "metadata": {},
   "source": [
    "### Personal Wellbeing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e2fa0293-a059-41f3-9f11-6e58094c311f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load data\n",
    "# https://data.london.gov.uk/dataset/subjective-personal-well-being-borough\n",
    "wellb0 = pd.read_excel(\n",
    "    'https://data.london.gov.uk/download/subjective-personal-well-being-borough/d1bae8e8-f07b-4c2f-aec7-deed8006b52a/personal-well-being-borough.xlsx', \n",
    "    sheet_name='Summary - Mean Scores',\n",
    "    header=[0,1])            # Where are the column names... there's three of them!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ba63023a-eeed-4e53-91d8-135891119208",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Clean data\n",
    "wellb = wellb0[[(' ', 'Code'), \n",
    "               (' ', 'Area'),\n",
    "               ('Life Satisfaction', '2018/19'),\n",
    "               ('Worthwhile', '2018/19'),\n",
    "               ('Happiness', '2018/19'),\n",
    "               ('Anxiety', '2018/19')]]\n",
    "wellb = wellb.droplevel(1, axis=1)\n",
    "wellb.columns = ['Code', 'Area', 'Life Satisfaction', 'Worthwhile', 'Happiness', 'Anxiety']\n",
    "wellb.dropna(inplace = True)\n",
    "wellb = wellb.iloc[:32]\n",
    "wellb['Happiness'] = pd.to_numeric(wellb['Happiness'])\n",
    "wellb['Anxiety'] = pd.to_numeric(wellb['Anxiety'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7532cfa-69bf-4391-a555-bc2af219039b",
   "metadata": {},
   "source": [
    "### GreenBlue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "2bcb1242-35ee-4ecd-98d8-95f7492296f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load data\n",
    "# https://data.london.gov.uk/dataset/green-and-blue-cover\n",
    "greenblue0 = pd.read_excel(\n",
    "    'https://data.london.gov.uk/download/green-and-blue-cover/7056efeb-d34a-4394-9d62-2f335edeacec/green_cover_ward_summary_0.05.xlsx',\n",
    "    sheet_name='ward_green_cover')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a6a6295e-afce-415d-99ca-2ea51e332b94",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/conda/envs/sds2021/lib/python3.9/site-packages/pandas/util/_decorators.py:311: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  return func(*args, **kwargs)\n"
     ]
    }
   ],
   "source": [
    "# Clean data\n",
    "greenblue = greenblue0[['lb_name', 'lb_code', 'ward_code', 'ward_name',\n",
    "       'percent_green', 'percent_blue', 'percent_green+blue']]\n",
    "#greenblue['percent_blue'] = np.log10(greenblue['percent_blue'])\n",
    "#greenblue['percent_blue'].replace([np.inf, -np.inf], np.nan, inplace=True)\n",
    "greenblue.dropna(inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "523bfb10-7ef6-4fac-b987-f865c138f3ad",
   "metadata": {},
   "source": [
    "### Load Access data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "4d7666fc-4b33-4fd1-8011-24c6ef5cd2af",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load data\n",
    "# https://data.london.gov.uk/dataset/access-public-open-space-and-nature-ward\n",
    "access_h0 = pd.read_excel(\n",
    "    'https://data.london.gov.uk/download/access-public-open-space-and-nature-ward/18fc0bef-68af-42c1-b77c-d6babbb66465/access-public-open-space-ward.xls',\n",
    "    sheet_name='Access to open space 2013 wards',\n",
    "    header=[0,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "73a3e6ce-e8dd-4f63-8404-540ffd346c33",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Clean data\n",
    "\n",
    "# Units: Percentage of households with access to: each type of park\n",
    "access_h= access_h0.droplevel(0, axis=1)\n",
    "access_h = access_h.iloc[:625]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "13439975-add2-4602-80d2-f41b5477c8a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# https://data.london.gov.uk/dataset/access-public-open-space-and-nature-ward\n",
    "access_o0 = pd.read_excel(\n",
    "    'https://data.london.gov.uk/download/access-public-open-space-and-nature-ward/18fc0bef-68af-42c1-b77c-d6babbb66465/access-public-open-space-ward.xls',\n",
    "    sheet_name='Open space 2013 wards',\n",
    "    header=[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "8cd94ee1-21fa-45af-89a4-da04ddaba982",
   "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>Ward_GSS_CODE</th>\n",
       "      <th>Ward_NAME</th>\n",
       "      <th>BOROUGH_LB_GSS_CD</th>\n",
       "      <th>BOROUGH_Name</th>\n",
       "      <th>% of open space that has access</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>619</th>\n",
       "      <td>E05000645</td>\n",
       "      <td>Tachbrook</td>\n",
       "      <td>E09000033</td>\n",
       "      <td>Westminster</td>\n",
       "      <td>64.713326</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>620</th>\n",
       "      <td>E05000646</td>\n",
       "      <td>Vincent Square</td>\n",
       "      <td>E09000033</td>\n",
       "      <td>Westminster</td>\n",
       "      <td>58.893270</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>621</th>\n",
       "      <td>E05000647</td>\n",
       "      <td>Warwick</td>\n",
       "      <td>E09000033</td>\n",
       "      <td>Westminster</td>\n",
       "      <td>20.429739</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>622</th>\n",
       "      <td>E05000648</td>\n",
       "      <td>Westbourne</td>\n",
       "      <td>E09000033</td>\n",
       "      <td>Westminster</td>\n",
       "      <td>41.628560</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>623</th>\n",
       "      <td>E05000649</td>\n",
       "      <td>West End</td>\n",
       "      <td>E09000033</td>\n",
       "      <td>Westminster</td>\n",
       "      <td>38.690776</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    Ward_GSS_CODE       Ward_NAME BOROUGH_LB_GSS_CD BOROUGH_Name  \\\n",
       "619     E05000645       Tachbrook         E09000033  Westminster   \n",
       "620     E05000646  Vincent Square         E09000033  Westminster   \n",
       "621     E05000647         Warwick         E09000033  Westminster   \n",
       "622     E05000648      Westbourne         E09000033  Westminster   \n",
       "623     E05000649        West End         E09000033  Westminster   \n",
       "\n",
       "     % of open space that has access  \n",
       "619                        64.713326  \n",
       "620                        58.893270  \n",
       "621                        20.429739  \n",
       "622                        41.628560  \n",
       "623                        38.690776  "
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "access_o = access_o0[['Ward_GSS_CODE','Ward_NAME','BOROUGH_LB_GSS_CD','BOROUGH_Name',\n",
    "                      '% of open space that has access']]\n",
    "access_o = access_o.iloc[:624]\n",
    "access_o.tail()\n",
    "# no city of London, 'Ward_GSS_CODE'=='E09000001'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "931827ba-f7bb-4bd8-a75d-b7102cb52172",
   "metadata": {},
   "source": [
    "## Merging Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "e991ded5-6171-469f-bfdb-64717b128e13",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(625, 8)\n",
      "(624, 5)\n",
      "(657, 7)\n",
      "(32, 6)\n",
      "\n",
      "(33,)\n",
      "(32,)\n",
      "(33,)\n",
      "(32,)\n",
      "\n",
      "(604,)\n",
      "(624,)\n",
      "(638,)\n"
     ]
    }
   ],
   "source": [
    "#Checking data before merging\n",
    "\n",
    "#Greater London is composed of London boroughs, the 32 local authority districts\n",
    "#City of London has 25 wards\n",
    "print(access_h.shape)\n",
    "print(access_o.shape)\n",
    "print(greenblue.shape)\n",
    "print(wellb.shape)\n",
    "#print(atlas.shape)\n",
    "print('')\n",
    "print(access_h['Borough name'].unique().shape)\n",
    "print(access_o['BOROUGH_Name'].unique().shape)\n",
    "print(greenblue['lb_name'].unique().shape)\n",
    "print(wellb['Area'].unique().shape)\n",
    "print('')\n",
    "print(access_h['Ward name'].unique().shape)\n",
    "print(access_o['Ward_GSS_CODE'].unique().shape)\n",
    "print(greenblue['ward_name'].unique().shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "624d687c-77a4-48fa-a9b9-e676fb51c0ca",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(624, 13)\n",
      "(600, 20)\n",
      "(600, 26)\n",
      "(600, 16)\n"
     ]
    }
   ],
   "source": [
    "#Merging\n",
    "\n",
    "df = pd.merge(access_h, access_o, left_on='WD13CD', right_on='Ward_GSS_CODE', how='inner')\n",
    "print(df.shape)\n",
    "\n",
    "#df = pd.merge(df, greenblue, left_on='WD13CD', right_on='ward_code', how='inner')\n",
    "df = pd.merge(df, greenblue, left_on='Ward_NAME', right_on='ward_name', how='inner')\n",
    "print(df.shape)\n",
    "\n",
    "df = pd.merge(df, wellb, left_on='lb_code', right_on='Code', how='left')\n",
    "print(df.shape)\n",
    "\n",
    "df = df[['WD13CD', 'Ward name', 'Borough name', \n",
    "         \n",
    "         # Mean per ward = borough mean * (pop ratio of ward /borough)\n",
    "         'Life Satisfaction', 'Worthwhile', 'Happiness', 'Anxiety',\n",
    "         \n",
    "         # Percentage of households in ward with access to\n",
    "         'Open Space', 'Local Parks', 'District Parks', 'Metropolitan Parks', 'Regional Parks',\n",
    "         \n",
    "         # Per ward\n",
    "         '% of open space that has access',\n",
    "         \n",
    "         # Coverage per ward\n",
    "         'percent_green', 'percent_blue', 'percent_green+blue']]\n",
    "             \n",
    "print(df.shape)\n",
    "\n",
    "#df.isnull().sum(axis=0).sort_values(ascending=False)\n",
    "df.dropna(inplace=True)\n",
    "#df.dtypes\n",
    "#df.isnull().sum(axis=0).sort_values(ascending=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "41dd05dc-d09f-409e-ad01-dfe8ffdfe32c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Some subsets for manipulation\n",
    "\n",
    "dfplot = df.drop(columns=['WD13CD', 'Ward name', 'Borough name'])\n",
    "\n",
    "df0 = df[['Life Satisfaction', 'Worthwhile', 'Happiness', 'Anxiety']]\n",
    "\n",
    "df1 = df[['Open Space', 'Local Parks', 'District Parks', 'Metropolitan Parks', 'Regional Parks',\n",
    "                     '% of open space that has access']]\n",
    "\n",
    "df2 = df[['percent_green', 'percent_blue', 'percent_green+blue']]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5c61b3e-6119-4541-a35b-2ab9a110edbb",
   "metadata": {},
   "source": [
    "## Plotting Distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "0d84a212-88f5-41b8-a0c8-65f33dfc6f3b",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/conda/envs/sds2021/lib/python3.9/site-packages/seaborn/_decorators.py:36: FutureWarning: Pass the following variable as a keyword arg: x. From version 0.12, the only valid positional argument will be `data`, and passing other arguments without an explicit keyword will result in an error or misinterpretation.\n",
      "  warnings.warn(\n",
      "/opt/conda/envs/sds2021/lib/python3.9/site-packages/seaborn/_decorators.py:36: FutureWarning: Pass the following variable as a keyword arg: x. From version 0.12, the only valid positional argument will be `data`, and passing other arguments without an explicit keyword will result in an error or misinterpretation.\n",
      "  warnings.warn(\n",
      "/opt/conda/envs/sds2021/lib/python3.9/site-packages/seaborn/_decorators.py:36: FutureWarning: Pass the following variable as a keyword arg: x. From version 0.12, the only valid positional argument will be `data`, and passing other arguments without an explicit keyword will result in an error or misinterpretation.\n",
      "  warnings.warn(\n",
      "/opt/conda/envs/sds2021/lib/python3.9/site-packages/seaborn/_decorators.py:36: FutureWarning: Pass the following variable as a keyword arg: x. From version 0.12, the only valid positional argument will be `data`, and passing other arguments without an explicit keyword will result in an error or misinterpretation.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA/0AAAEcCAYAAACVuNRuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAB3H0lEQVR4nO3deVhU1f8H8PcMm+DCooDgWn4TzTJhWCwgFXdlczeXSir3JSuV3HdFK7dwS6Ov2WoKKmqWW6aVSqFmLinuyiKLCsgyzNzfH/y4X0ZgmIHZGN6v5+F5mLnb59y559x77j3nXIkgCAKIiIiIiIiIyOxIjR0AEREREREREekHK/1EREREREREZoqVfiIiIiIiIiIzxUo/ERERERERkZlipZ+IiIiIiIjITLHST0RERERERGSmam2l39PTE3fu3DHoNtPT09G7d28UFhYadLs10dy5cxEdHS1+/vrrr/HKK6/A09MTWVlZRoxM82NHXzEnJCSgZ8+eOlufOdNXPo+MjMSqVavKnVZYWIhevXohIyND59slzXh4eODWrVvVWse6devwwQcfVDi9b9++OHXqlEbzkmkz9PUArwVMy8aNGzFr1ixjh0FGZqhyYNeuXXjttdcA8HrBkO7fvw9PT08oFApjh2I0Zl/pDwoKQvv27eHp6Sn+paamIjExEc2aNdN6fYWFhZg8eTKCgoLg4eEhXvRpolGjRvDz88N3331X4TyRkZF44YUXxFiDg4Px8ccfIzs7W+tYtTV37lxMnz69zPeXL1/GCy+8gIcPH+p0e6dOncKrr75a7rSFCxdiwoQJAAC5XI7ly5fj888/R2JiIhwdHbXaztWrVxEREQEfHx94e3ujf//++OWXXzRaduTIkdixY4fKd5ocO9WNubSnKzDe3t44ePBglddnjkwpn1tbW2PAgAH47LPPKpxn3bp1aNeunRhrz549sXDhQqSlpWkdq7Y2bdqE4cOHl/k+MzMTL7zwAv7991+dbu+tt95S2Repqanw8PAo97sHDx5ovf7y8qgh7Nu3D35+fgbfLlWdqZQTvBZQVdG1gKHy9tixY7FkyRK9b4dMg6mUAwCvFyry5MkTeHp64p133tHZOt3d3ZGYmAgLC4tK5y19Y8acmH2lHyi+i5uYmCj+ubq6Vmt9Xl5eWLFiBZydnbVeNiQkRO2JHii+SE5MTMQff/yBpUuX4uzZs3jttdfw5MmTqoaskX79+uHnn38us53du3ejS5cucHBw0HhdRUVFOosrIyMDBQUF+M9//lOl5ceOHYtXXnkFJ0+exG+//YZZs2ahbt26OouvPNWNmbRnavk8NjZW7ZO83r17IzExEadPn8ann36K9PR09O/fX+8n8rCwMCQmJpZ5orF//360bt0arVu31nhdgiBAqVSqncfb2xtnzpwRP585cwbPPvtsme9atmyp1b7WZNtETzOVcoLXAkTGYyrlAMDrhfIcPHgQ1tbWOHnypEFubtQWtaLSX57ST0+zsrIwduxYeHl5YcCAAVi1alWFd3isra3x5ptvwtvbG1Jp2d137NgxhIeHw8vLC506dcK6detUpr/00ku4c+cO7t27V2mMNjY2aN++PTZs2ICHDx9i165d4rQffvgBvXv3ho+PD9566y2V9Xl4eGDbtm3o2rUr/Pz8EBUVpVFG8/T0hIuLC3766SfxO4VCgb179yI8PBy3b9/G66+/Dj8/P/j5+eH999/H48ePxXmDgoKwefNmhISEoEOHDtU62Zc0n75x4wZ69eoFAPDx8cHrr78OAEhKSsKoUaPg6+uLnj17Yv/+/eWuJzMzE3fv3sXgwYNhbW0Na2tryGQyeHt7AwAePXqEMWPGoGPHjvDx8cGYMWOQkpICAFi1ahUSEhKwcOFCeHp6YuHChQBUj51ffvkFffr0gaenJwIDA7F169YKY168eDE6deoELy8v9O/fHwkJCSr7eePGjejWrRs8PT3Rv39/JCcni3dYw8LC4Onpif3795d5KpKUlISRI0fC29sbffv2xeHDh1X244IFCzB69Gh4enpi0KBBuH37dpV/l5pGX/m8ZH2jRo2Cp6cnRowYoZIHGzduDHt7e5w9e7bSGK2srPDcc89h1apVcHJyQkxMjDjt6NGjCAsLg7e3N4YOHYrLly+L04KCgrBp0yb06dMHPj4++PDDD1FQUFDp9ho3boyOHTti9+7dKt/HxcUhPDxcbZ4Aip++rVq1CkOHDhXLM3V8fHzw119/iWVQQkIC3njjDVy4cEHlu5I8+ddff2HAgAGQyWQYMGAA/vrrrwq3PW3atHLzKAD89ttv6NGjB3x8fLBgwQIIggAA6NKlCy5cuACguBLj4eGBa9euAQB27NiB8ePHi+uQy+WYPn06PD090bdvX/z9998q+/+3334rN81nz57F0KFD4e3tjdDQUJUnPyNHjsTq1asxdOhQeHp6IiIiApmZmRotu2vXLnTt2hWenp4ICgrCnj17AAC3bt3CiBEjIJPJ4Ofnh3fffVftb0Kq9FVO5OfnY/ny5ejSpQtkMhlee+015OfnA+C1gLY0KZc+/vhjDBw4EDKZDOPGjRNbJNy9exceHh747rvvEBAQgICAAHz++efisqW755TMGxsbi86dO8PPzw8bNmwQ51Uqldi8eTO6desGPz8/TJkyRdxOQUEBPvjgA/j5+cHb2xsDBgxAeno6gIrzLpkOfZUDpdc1cODAMtdgvF4oKzY2FkOHDoWHhwf27t2rMi0oKAhbt25FSEgIZDIZ3n33XTEtmzdvxuDBg8Uy5+uvv0bfvn1RUFAg5u2SadnZ2Zg5cyYCAgIQGBiIVatWQaFQICkpCfPmzcPZs2fh6ekJb29vnD9/Hq+88opKWXbw4EGEhYVplB5TUWsr/aUtXLgQtra2OHnyJKKiohAXF1flddna2iIqKgoJCQnYtGkTvvnmGxw6dEicbmlpiebNm6tkxsrUq1cPr7zyilhJPHToEDZt2oRPP/0Uv//+O2QyGd5//32VZX7++Wfs3LkTsbGxOHLkCHbu3KnRtsLDw1XS/9tvv6GoqAivvvoqBEHAmDFj8Ouvv+LAgQNISUkpc1Nj37592Lx5MxISEmBpaYmQkJAyGVYbzzzzDOLj4wEUPw3ctm0bnjx5goiICAQHB+O3337DJ598ggULFuDq1atllnd0dESLFi0wbdo0HDp0SDwBl1Aqlejfvz+OHj2Ko0ePwsbGRqw4TJ06Fd7e3pg7dy4SExMxd+7cMuufNWsWFi5ciMTERMTHx6Njx47lxgwAL774IuLi4nD69GkEBwdjypQpYkEVExMj7ru//voLS5cuRZ06dfDVV18BKK6cJCYmok+fPirbl8vlGDt2LPz9/fHbb79h9uzZ+OCDD3D9+nVxnn379mHixIk4c+YMmjdvXmFfdHOny3wOAHv37sX48eNx6tQptGnTpkyf7meffVarfG5hYYGuXbuK+fyff/7BzJkzsXDhQpw6dQpDhgzB+PHjVZ4G7N27F1u3bsXPP/+MGzduYP369RptKzw8XOWi8/r167h8+TKCg4PV5okSu3fvxqJFi/DXX3/B3d0dY8aMwebNm8vdVvv27VFYWCjui4SEBPj7+6uUgwkJCfDx8cHDhw8xZswYjBw5EqdOncKoUaMwZswYlTExSm97+fLlFebRY8eO4YcffsDu3btx4MAB/PrrrwCKb0KcPn1a3G6zZs3Ez2fOnIGvr6+4jiNHjqBv375ISEhAUFAQFi1aVOm+TU1NxZgxYzBu3DicPn0aM2bMwOTJk1Uq9vHx8Vi2bBl+//13yOVysQKibtknT55g8eLF+Oyzz5CYmIhvv/0Wbdu2BQCsWbMG/v7+OHPmDI4fP44RI0ZUGieVT5flRFRUFP755x98++23OH36NKZNmyZWCngtoB1NyqW4uDgsXboUv/76KywtLbF48WKV6adOncJPP/2ErVu3YvPmzRXetAOAP//8Ez/++CP++9//Ijo6GklJSQCAbdu24dChQ9i+fTt+/fVX2Nvbi3HExsYiJycHx44dw6lTp7BgwQLUqVNHbd4l06TLcmDhwoWwsbHBiRMnsHTp0nLzIK8X/uf+/fs4ffo0QkJCEBISUu6+P3DgALZs2YLDhw/jypUr4o3Qt99+G1ZWVtiwYQNu3ryJVatWYeXKlbCxsSmzjhkzZsDS0hI//fQT4uLicPLkSezYsQOtWrXCggUL0KFDByQmJiIhIQHt27eHg4MDTp48KS6/Z88eVvpN0YQJE+Dt7Q1vb2+VpzhA8d3rn376CZMmTYKtrS3+85//IDw8vMrb8vPzg4eHB6RSKdq0aYO+ffuKF5Ql6tatq3W/PBcXFzx69AgA8O2332L06NFo1aoVLC0tMXbsWFy6dEnlDv8777wDBwcHuLu74/XXXxcroZUJCwvDmTNnxDt1cXFxCA4OhpWVFVq0aAF/f39YW1vDyckJo0aNUmmiCxTf1XNzc0OdOnUAFBcyISEhWqW1MseOHUOTJk0wYMAAWFpaol27dujZs2e5/dwlEgm2bduGJk2aYPny5QgICMDw4cNx8+ZNAMU3BXr27AlbW1vUq1cP48aNK5MmdSwtLXHt2jXk5OTA3t4e7dq1q3DesLAwODo6wtLSEhERESgsLMSNGzcAFD9hnDJlCp599llIJBK0adNGo3EAzp07hydPnmD06NGwtrbGyy+/jC5dumDfvn3iPN27d0f79u1haWmJ0NBQXLp0SeP01SSGzOcA0LlzZ/j4+MDa2hpTp07F2bNnkZycLE6vW7euytMvTZTO599//z2GDBmCl156CRYWFujXrx+srKxUngYMHz4cbm5ucHBwwLhx41R+d3W6d++O9PR08Sn67t27ERgYCCcnJ43yRL9+/fDcc8/B0tISVlZW2LRpE0aPHl3utqytrfHSSy/hzJkzePjwIR4/foxmzZqJzf4fPnyIa9euwcfHB8eOHUOLFi0QHh4OS0tLBAcH49lnn8XRo0cr3HZF3nnnHTRo0ADu7u7w8/MTL6iervSPGTNGTN+ZM2fg4+MjrkMmk6FTp06wsLBAWFiYRhdlu3fvxquvvopOnTpBKpXC398fL7zwgso4Iv3798czzzyDOnXqoFevXmKerGxZqVSKq1evIj8/Hy4uLnjuuecAFJdD9+/fR1paGmxsbMRWE1SWocoJpVKJnTt3YtasWXB1dYWFhQW8vLxgbW0tzsNrgf9JS0sTf5eSvz///FOcrkm5FBYWhtatW8POzg5TpkzBjz/+qDJo14QJE2BnZwcPDw/0799f7b6YOHEi6tSpgzZt2qBNmzZi3v/uu+8wdepUNG7cGNbW1pg4cSIOHjyIoqIiWFpa4uHDh7h16xYsLCzwwgsvoF69egAqzrtkHIYqB0rWNXnyZNjZ2aF169bo169fmfl4vfA/cXFx8PDwwH/+8x/07dsX165dw8WLF1XmGTlyJFxdXeHg4IAuXbqI51CpVIqoqCh8+eWXGDduHN5++208//zzZbaRnp6O48ePY+bMmbCzs0PDhg3x5ptvqt0npW9+PHz4ECdOnEBwcLBG+9BUWBo7AEOIjo7GK6+8Uu60zMxMFBUVwc3NTfyu9P/aOnfuHD766CNcvXoVcrlcHJmztNzcXNSvX1+r9aampsLe3h5A8V2wpUuXIioqSpwuCAJSU1PRpEmTMmlo0qSJxn1i3N3d4e3tjT179mD48OE4fPiw+LQ5IyMDixcvRkJCAnJzcyEIAho0aKCyfHX2nabu3buH8+fPq1zYKhQKhIaGljt/48aNxSeAycnJmDNnDmbMmIHvvvsOeXl5WLZsGX799Vex8MzNzYVCodBosI+1a9diw4YN+Pjjj+Hh4YH3338fnp6e5c77+eefY8eOHUhLS4NEIkFOTo74BDMlJQXNmzfXaj8AxRdKjRs3VmlS5u7ujtTUVPFzo0aNxP9LnjqYI0Pmc6D4uCpRt25d2NvbIy0tTVxvbm5umfxRmafzeVxcHLZv3y5Ol8vlKnm5dBrc3d01zue2trbo1asX4uLi4Onpib179yIyMhIANMoT2u47b29vJCQkoEmTJpDJZACKK9S7du1CkyZN4ObmhiZNmmDfvn1wd3dXWfbp41nTbZfuW2lra4vc3FwAgK+vL1asWIEHDx5AqVSid+/e+PTTT3H37l1kZ2erPIF7Ou8UFBSIF/cVuX//Pn788UeVGxVFRUUqg/49HVtJnlS3rJ2dHVatWoXPP/8cs2bNgpeXF2bMmIFWrVph2rRpWLNmDQYOHAh7e3uMGjUKAwcO1Gg/1TaGKieysrJQUFCgdmAwXgv8j4uLC44fP67y3ciRI8X/tS2X3N3dIZfLVVoJPb0v1A1CVjrvP51HJ0yYoHLOlUqlyMjIQFhYGFJSUvDee+/h8ePHCA0NxdSpU9XmXTIOQ5UD5a3r6XMcwOuF0nbv3o1BgwYBAFxdXeHj44PY2FiVyvvT59DSaWnatCn8/Pzwyy+/lDsIIVC8v4qKihAQECB+p1Qq1cYaFhaG3r17Izc3FwcOHIC3tzdcXFy0Spux1YpKvzpOTk6wtLRESkoKnnnmGQBQeVqnrffffx8jRozAli1bYGNjgyVLlqicdIqKinD79m20adNG43Xm5ubi999/x9ixYwEUZ6CxY8dWWMktSUPJneT79+9rdWD269cPmzdvhrOzM5o2bSo+vf74448hkUiwZ88eODo64tChQ2Wa8UgkEo23U1Vubm7w8fFR6cukzbLDhw/He++9B6C4In7jxg18//33cHZ2xqVLlxAeHi72/61MST9LuVyOr776Cu+++265bwZISEjAZ599hi+++ALPPfccpFIpfHx8xO00btwYt2/f1mpQFKD4QiklJQVKpVK8CElOTkbLli21Wo+503U+B6DSby03NxePHj1SyWfXr19HRESExutTKpU4evSoeCFSks/HjRtX4TKl01CVfD5hwgT06NEDubm56Ny5MwDN8oS2+dzHxwfffvstmjRpIt6s8/LywuzZs1W+c3Fxwf3798ukMTAwsMrbflqLFi1Qp04dfPnll/D29ka9evXQqFEjfP/995DJZBWO4aApNzc3hIWFlWlarItlAwMDERgYiPz8fKxevRpz5szB119/DWdnZ3GZhIQEjBo1Cj4+PmjRokW10lLb6LKccHR0hI2NDe7cuVPu+Z7XAtrRpFwq/VslJyfDysoKjo6O4vfJycliRVvbfVGicePGWLp0qXjz8mkTJ07ExIkTcffuXYwePRrPPPMMBg0aVGHeJdOjy3KgZF2lj73y1sXrhWJ//fUXbt68ic2bN4vX+Lm5ubh27ZrYHL8yv/zyCxITE/Hyyy9jxYoVZcomAGJLnT/++KPcdZYXs6urKzw9PfHzzz9j9+7dNXJ0/1rRvF8dCwsLdO/eHZ9++iny8vKQlJRUZsCKpxUWFop9seVyOQoKCsQDPDc3F/b29rCxscH58+fLNB87f/48mjRpIt6Fr2w7Fy5cwIQJE9CgQQP0798fADB06FBs3rxZ7MOenZ2NAwcOqCy7detWPHr0CMnJydi2bVuZvuDq9OjRA8nJyVi3bp1Kk6bc3FzY2dmhQYMGSE1NxZYtWzRepzoFBQUqf5VVuDt37oybN28iLi4Ocrkccrkc58+fF/vclfbo0SOsXbsWt27dglKpRGZmJnbu3IkOHTqIabKxsUGDBg3w8OFDfPrppyrLN2rUqMKBRwoLC7Fnzx5kZ2fDysoKdevWrbB1QG5uLiwsLODk5ISioiJ8+umnyMnJEacPGjQIa9aswc2bNyEIAi5fvizeLFIXQ/v27WFra4stW7ZALpfj1KlTOHLkiFa/d22g63wOFJ9YEhISUFhYiDVr1uCll14S7xKnpqbi0aNH4nGmjlwuR1JSEt577z2kp6fjzTffBFB8THz77bc4d+4cBEHAkydPcOzYMZXj5uuvv0ZKSgoePnwoDtKjKW9vb9SvXx9z585Fnz59xGbHleWJqvD09ER2djb27NkjXizb29vDyckJe/bsEZvUd+rUCTdv3sTevXtRVFSE/fv349q1a+IFRnnU5Y+K+Pr6Yvv27eJ2n/5cHaGhoTh69Ch+/fVXKBQKFBQU4NSpUyo3iaqybHp6Og4fPownT57A2toadnZ2YnlT0q8aKN6vEomk2jcvaiNdlhNSqRQDBgzAsmXLkJqaCoVCgcTERLGPLa8FtKNJubRnzx5cu3YNeXl5WLNmDXr27KlyTl6/fj3y8vJw9epV7Nq1q0rnyddeew2rV68Wu1BkZmaK4zb98ccfuHLlChQKBerVqwdLS0tYWFiozbtkenRZDjy9rmvXriE2NlZlWV4v/E9cXBz8/f2xb98+xMXFIS4uDnv37kVeXl6ZlkDlyczMxKxZs7BkyRIsX74cR44cKfdBnIuLC/z9/bF8+XLk5ORAqVTi9u3bYte/hg0bIjU1tcwbFcLCwrB161b8+++/6N69e7XSagy8KkDxO2mzs7Ph7++P6dOno2/fvir97vr27asyiEWvXr3Qvn17pKam4q233kL79u3FE8C8efOwdu1aeHp6Ijo6Gr1791bZ1t69ezF06FC18WzduhWenp7w9fXFjBkz0K5dO3z77bews7MDUNy/5u2338Z7770HLy8vBAcHl8kMXbt2Rf/+/REeHo7OnTuLTT0TEhIqbH5ews7ODj179kRycrJKH7yJEyfi4sWL8Pb2xujRo9GjRw+16wHK7runpaamon379ip/lY0uX69ePWzduhX79+9HYGAgAgIC8NFHH5X7uhMrKyvcu3cPo0aNgkwmQ0hICKytrbF8+XIAwBtvvIGCggJ07NgRQ4YMUXmiCACvv/46Dh48CB8fn3Kfvu3evRtBQUHw8vLCt99+ixUrVpQbc0BAAF599VX07NkTQUFBsLGxUWlGNGrUKPTu3RsRERHw8vLCrFmzxBPIxIkTERkZCW9v7zJvKbC2tsaGDRtw/PhxdOzYEQsWLMCKFSvYbLAcusznABAcHIzo6Gj4+fnhn3/+wcqVK8VpJaNcl17/0w4cOABPT0/4+Phg3LhxcHBwwK5du8RXB7344otYtGgRFi5cCB8fH/To0UNl1O6SGCIiItCtWzc0a9ZMvMt///59eHp6lnlqXppEIkF4eDju3bunckFfWZ4oz9tvv42NGzdWON3W1hbt2rVDYWGhSmsWmUyGjIwM8Um/o6MjNm7ciJiYGPj5+WHLli3YuHEjnJycKlx3ZXm0PD4+PsjNzVWp9Jf+XB1ubm5Yv349Nm3ahJdffhmdOnXC1q1bNRo1Xd2ySqUSMTExCAwMhK+vL86cOYN58+YBAP7++28MGjQInp6eGDduHGbNmlWl902TbsuJGTNmoHXr1hg4cCB8fX3x0UcficcBrwW0o0m5FBYWhsjISPj7+6OwsBCzZs1Sme7r64vu3bvjzTffREREhErTXk29/vrrCAoKQkREBDw9PTF48GCcP38eQHE/4cmTJ0Mmk6FPnz7w9fVFaGio2rxLpkmX5cDcuXPx5MkT+Pv7IzIyUrxpV4LXC8UKCgpw4MABjBgxAs7OzuJfs2bNEBYWptFginPnzkVQUBA6deoER0dHLFmyBLNmzVJpcV1ixYoVkMvl4hsNJk+ejAcPHgAAOnbsiP/85z8ICAhQ6ZrXvXt33Lt3D927dxfL4ZpEImjajrkWWblyJdLT01X6yelCRkYGRowYgbi4uHJHktQVDw8P/PTTT2zaSaSGvvJ5YWEhQkND8dVXX6Fhw4Y6XXdpQUFBWLx4cYX9Eomo+vRRTvBaQPdGjhyJ0NBQsS9waXfv3kXXrl3xzz//aNQ8mOhpvF6gEt26dcPChQtr5L7kk34Uv+f88uXLEAQB58+fxw8//KCXZhsNGzbEgQMH9HqSJ6LyGSqfW1tb48cff9TrCZyI9MMQ5QSvBYhMG68XqDwHDx6ERCJBx44djR1KlfCWJ4r7pLz//vtIS0tDw4YNERERga5duxo7LCLSIeZzIqoMywkiYjlATxs5ciSuXbuGFStW1Ngxc9i8n4iIiIiIiMhM1cxbFURERERERERUKVb6iYiIiIiIiMwUK/1EREREREREZsrsB/LLysqFUimgYcN6yMjIMXY4esG01UzGSptUKoGjY12Db1dfSvK4uTDnY74ytTXt+kg387lp47Fe+zCfV84Y+dxUjklTiMMUYmAcZWPQFbOv9CuVgliAmNMFw9MMlbavv96GO3duab3co0cPAQD29g5aLWdlZQG5XKH19nShWbMWGDbsdb1uw5yPSUMpncfNhT7TU9U8rA1D5XdD5FFDMbdjWNdqYj5Xl9eMeW4rUdV8Wh1Pp9uc8rAmatoxbGjGyuem8rsY+1peF+WSLsoVdXEYuswwlWNDF8y+0k+6defOLVy5eg0WdRy0Wk6R/xAA8OBxke6D0oOSeInMTVXzsDYMkd+ZR8nUGSKvVYexz8vMw0TGoc+ySZ/lCsuM6mGln7RmUccBdi20e1/pk1uHAUDr5YylJF4ic1SVPKwNQ+R35lGqCfSd16rD2Odl5mEi49FX2aTPcoVlRvVwID8iIiIiIiIiM8VKPxEREREREZGZYqWfiIiIiIiIyEzV6kr/yZPHcfLkcWOHQaQzPKZVcX+QOeJxrYr7g8wRj+uahb8XacvQx0ytHsjvxIlfAAD+/q8aORIi3eAxrYr7g8wRj2tV3B9kjnhc1yz8vUhbhj5mavWTfiIiIiIiIiJzVquf9BMRERERkWkbP/4t5Ofnwda2LqKjP9N6OlFtxyf9RERERERksvLz8wAAeXm5VZpOVNux0k9ERERERCZp/Pi3VD5PmPCOVtOJqJY373/06CEePXqEqKhFxg6lWqysLCCXKwyyrdu3b0GpsDDItoxJWZSP27dv6fXY0Mfvdvv2Ldjb2+t0nTVZTczj+s7P5pKHDZFHDUXb35z5XJWp5nNzyWv6Yk55WBPM51VX8hS/xNNP8yubbgjVLYd4LV85Q5cZhrgeM2Qe55N+IiIiIiIiIjNVq5/029s7wN7eATNmzDF2KNXi7FwfDx5kG2RbUVGLcO1OukG2ZUxSyzpo3qyRXo8NffxuteWJiaZqYh7Xd342lzxsiDxqKNr+5sznqkw1n5tLXtMXc8rDmmA+N2/VLYd4LV85Q5cZhrgeMyQ+6SciIiIiIpNUp46tymdb27paTSciVvqJiIiIiMhErV+/VeXz06/kq2w6EbHST0REREREJqzkaX5FT/Erm05U29XqPv1ERERERGTann6ar+10otquVlf6AwI6GTsEIp3iMa2K+4PMEY9rVdwfZI54XNcs/L1IW4Y+Zmp1pd/f/1Vjh0CkUzymVXF/kDnica2K+4PMEY/rmoW/F2nL0McM+/QTERERERERmSlW+omIiIiIiIjMFCv9RERERERERGaqVvfpp6pR5D/Ek1uHtV4GgNbLGUtxvI2MHQaRXlQlD2u7fkC/+Z15lGoCfee16jD2eZl5mMh49FU26bNcYZlRPaz0k1aaNWtRpeUePSo+1OztHbRazsrKAnK5okrbrJ5GVU4rkSkzxHFtmPzOPEqmTd3xabxz2/9UNZ9Wh2q6mYeJjKGifKeLckkX5UrFcbDMqA5W+kkrw4a9btDtOTvXx4MH2QbdJpE5M3Qe1gbzO5kTdXmtth7rtTXdRKakorLJVPKnqcRhbsy+0i+VSsr939wwbTWTMdJmbvvT3NIDmGeaNFVb067rdJvbfjS39ADmmSZN1NZ0A8znlTFWekxlP5pCHKYQA8A49EEiCIJg7CCIiIiIiIiISPc4ej8RERERERGRmWKln4iIiIiIiMhMsdJPREREREREZKZY6SciIiIiIiIyU6z0ExEREREREZkpVvqJiIiIiIiIzBQr/URERERERERmipV+IiIiIiIiIjNlaewAdOXu3buYMGGC+Dk7Oxs5OTk4ffq0ynzR0dHYv38/LCwsYGlpialTpyIwMNDQ4WpF07SVuH79Ovr164dhw4ZhxowZhgqzSrRJ2/79+7FhwwYIggCJRIKYmBg0atTIkOFqRdO0ZWRk4MMPP0RycjLkcjk6duyI2bNnw9LSbLInlWLOx3xlanueOHr0KNasWQNBEKBUKjFp0iT06NFDZR6FQoHFixfj119/hUQiwejRozFo0CAjRUxVZc7n7crU1jKutpdvpiQrKwvTp0/H7du3YW1tjRYtWmDhwoVwcnJSmW/dunX4+uuv4eLiAgDw8vLCvHnzAFS/LNY0hunTp+PKlSvi5ytXriA6Ohpdu3ZVG582xo8fj7t370IqlcLOzg5z5sxB27ZtVeZRl15dnZc0iUNdPU0X+0OTGPR5XGgTh86PDcFMLV68WFiwYEGZ748fPy48efJEEARBuHTpkiCTyYS8vDxDh1ctFaVNEAShqKhIGDFihPDee+8Jy5cvN3Bk1VdR2s6fPy/07t1bSEtLEwRBEB4/fizk5+cbOrxqqShtixcvFn+rwsJCYeDAgcK+ffsMHR4ZiTkf85WpTXlCqVQK3t7ewpUrVwRBKD7/dOjQQVAoFCrzxcbGChEREYJCoRAyMjKEwMBA4c6dO8YImXTInM/blamtZVxtKt9MTVZWlvDHH3+In5cvXy58+OGHZeZbu3ZthXmuumWxpjGUdunSJcHX11coKCioND5tPH78WPz/559/FsLDw8vMoy69ujovaRKHunqaLvaHJjHo87jQJo7SdHFsmGXz/sLCQuzduxcDBgwoMy0wMBC2trYAAA8PDwiCgIcPHxo4wqpTlzYA2Lx5Mzp37oyWLVsaNjAdUJe2L774AhEREXB2dgYA1K9fHzY2NoYOscrUpU0ikSA3NxdKpRKFhYWQy+VwdXU1QpRkaOZ8zFemNuYJqVSK7OxsAMVPAV1cXCCVqp6G9+/fj0GDBkEqlcLJyQndunXDjz/+aIxwSUfM+bxdmdpaxtXG8s2UODg4wM/PT/zcoUMH3L9/X6t1VLcsrkoMP/zwA0JCQmBtba1VrJWpX7+++H9OTg4kEkmZedSlV1fnJU3i0Hc9TZMY1DHkvihNF8eGWbYlOnLkCFxdXdGuXTu188XFxaF58+Zo3LixgSKrPnVpu3z5Mk6cOIFt27Zh/fr1RoiuetSlLSkpCU2bNsXw4cPx5MkTdO/eHePGjdM6sxqLurSNHz8ekyZNQkBAAPLy8jB8+HDIZDIjREmGZs7HfGVqW56QSCRYvXo1xo8fDzs7O+Tm5mLTpk1l5ktOToa7u7v42c3NDSkpKYYMlXTMnM/blamtZVxtK99MmVKpxDfffIOgoKByp+/btw8nTpyAs7MzJk2aBE9PTwC6LYsriwH4342iL774QqP4tDVr1iycPHkSgiBgy5YtZaarS68u90VlcZRWXj1NF/tDkxgMcVxoui90dWyY5ZP+nTt3VnhHvcTp06exZs0afPzxxwaKSjcqSptcLsecOXOwYMECWFhYGCGy6lP3uykUCly5cgUxMTH48ssvcfz4cezevdvAEVadurT9+OOP8PDwwIkTJ3D8+HEkJCTwyV4tYc7HfGVqW54oKirCpk2bsH79ehw9ehQbNmzA1KlTkZuba+zQSM/M+bxdmdpaxtW28s2ULVq0CHZ2dhgxYkSZaUOHDsXhw4exd+9evPXWWxg/fjyysrIMGkOJQ4cOwd3dXaVfty7jW7JkCY4dO4apU6dixYoVVVqHLmgaR3n1NF3tj8piMNRxoem+0NWxYXaV/tTUVJw5cwYhISEVzpOYmIhp06YhOjoazz77rAGjqx51aXvw4AFu376N0aNHIygoCP/973/x/fffY86cOUaIVHuV/W7u7u7o1asXrK2tUa9ePXTt2hXnz583cJRVU1natm/fjtDQUEilUtSvXx9BQUE4deqUgaMkQzPnY74ytTFPXLp0CWlpaeITPZlMBltbWyQlJanM5+bmptIENDk5uUa1RiNV5nzerkxtLeNqY/lmqqKionDr1i2sXr26TFcqAHB2doaVlRUAwN/fH25ubrh69SoA3ZXFlcVQorwbReriq6rw8HCcOnWqTAVRXXr1cV6qKA6g4nqarvdHRTEY4rjQJI4Sujo2zK7SHxsbi06dOsHR0bHc6efPn8fUqVOxdu3aSpv/mxp1aXN3d8epU6dw5MgRHDlyBG+88QYGDx6MRYsWGSFS7VX2uwUHB+PEiRMQBAFyuRx//PEH2rRpY+Aoq6aytDVt2hTHjx8HUNyE5/fff8dzzz1nyBDJCMz5mK9MbcwTjRs3RkpKCq5fvw6guGlzeno6mjdvrjJfr169sGPHDiiVSmRmZuLQoUPo2bOnMUImHTDn83ZlamsZVxvLN1O0atUqXLhwAdHR0RX2g05NTRX/v3TpEu7du4dnnnkGgG7KYk1iAICUlBT8+eefCA4O1jg+TeXm5iI5OVn8fOTIEdjb28PBwUFlPnXp1cW+0DQOdfW06u4PTWPQ93GhaRyAbo8Ns+vTHxsbi1mzZql8984772Dy5Ml48cUXsWDBAuTn52Pu3Lni9BUrVsDDw8PQoWqtsrTVZJWlrW/fvrhw4QL69OkDqVSKgIAADBw40EjRaqeytM2cORPz5s1DSEgIFAoF/Pz8MHjwYCNFS4Zizsd8ZWpjnnB2dsb8+fMxZcoUsc/ysmXL4ODgoJL2sLAwnDt3TnyV34QJE9CsWTNjhk7VYM7n7crU1jKuNpZvpubq1avYuHEjWrZsiaFDhwIovtkSHR2t8lt88skn+OeffyCVSmFlZYUVK1aIA0tWtyzWNAag+Jjp0qVLmUqfuvg0lZeXhylTpiAvLw9SqRT29vbYuHEjJBKJxuceXZyXNI1DXT2tuvtD0xj0eVxoEweg22NDIgiCoFWkRERERERERFQjmF3zfiIiIiIiIiIqxko/ERERERERkZkyuz79T8vKyoVSKaBhw3rIyMgxdjh6wbTVTMZKm1QqgaNjXYNvV19K8ri5MOdjvjK1Ne36SDfzuWnjsV77MJ9Xzhj53FSOSVOIwxRiYBxlY9AVs6/0K5WCWICY0wXD05i2msmc02YopfO4uTC39Gijtqa9tqZbU8zn5qO2phuo3WnXhLHyuan8LqYQhynEADAOfWDzfiIiIiIiIiIzxUo/ERERERERkZlipZ+IiIiIiIjITLHST0RERERERGSmWOknIiIiIiIiMlNmP3o/6UcDe1vYWGt3+BQUFuHxozw9RURENU155Yizc/1y52X5QaZO2/Niecc6j3Miqum0KQtZ5hmOQSr9UVFROHjwIO7du4e9e/eidevWyMrKwvTp03H79m1YW1ujRYsWWLhwIZycnAAAN27cQGRkJB4+fAgHBwdERUWhZcuWhgiXNGBjbYlpa37RapmVUzrpKRoiqomeLkesrCwhlxeVOy/LDzJ12pwXKzrWeZwTUU2nTVnIMs9wDNK8v2vXrvjqq6/QpEkT8TuJRIK3334bBw8exN69e9GsWTN89NFH4vR58+Zh2LBhOHjwIIYNG4a5c+caIlQiIiIiIiIis2GQSr+3tzfc3NxUvnNwcICfn5/4uUOHDrh//z4AICMjAxcvXkRwcDAAIDg4GBcvXkRmZqYhwiUiIiIiIiIyCybRp1+pVOKbb75BUFAQACA5ORmurq6wsLAAAFhYWMDFxQXJycli839NNWxYT/y/or6i5sAYabOy0v7wqUqc/N2IiIiIiIiqxiQq/YsWLYKdnR1GjBih83VnZORAqRTg7FwfDx5k63z9psAYaXN2rl9h31t1tI2Tv5vuSaUSlZthRERERERkvoz+yr6oqCjcunULq1evhlRaHI6bmxtSU1OhUCgAAAqFAmlpaWW6CBAREZFpOHr0KMLDwxEWFoaQkBD89NNPAIoH5h0yZAh69uyJIUOG4ObNm8YNlIiIqJYxaqV/1apVuHDhAqKjo2FtbS1+37BhQ7Rt2xbx8fEAgPj4eLRt21brpv1ERESkf4IgYPr06VixYgV2796NlStXYsaMGVAqlRyYl4iIyMgM0rx/8eLF+Omnn5Ceno5Ro0bBwcEBq1evxsaNG9GyZUsMHToUANC0aVNER0cDAObPn4/IyEisX78eDRo0QFRUlCFCJSIioiqQSqXIzi7uspSdnQ0XFxdkZWXh4sWLiImJAVA8MO+iRYuQmZnJG/lmqrJ3dGsylg3f3U1EpFsGqfTPnj0bs2fPLvP9lStXKlymVatW2LFjhz7DIiIiIh2QSCRYvXo1xo8fDzs7O+Tm5mLTpk06HZiXagZ17+i2srLUaDwgvrubiEi3TGIgPyIiIqq5ioqKsGnTJqxfvx4ymQx//vknpk6dihUrVuhk/TVl8FFt3mpT0bzm8FYXdftB031kDvvhaeaYJiKqGVjpJyIiomq5dOkS0tLSIJPJAAAymQy2trawsbERB+a1sLCo8sC8JW/iMWXavNVG3RPvmv7GGnX7QdMn/UDN3w9P08cbe/g2HiLSlNFH7yciIqKarXHjxkhJScH169cBAElJSUhPT0eLFi04MC8REZGR8Uk/ERERVYuzszPmz5+PKVOmQCKRAACWLVsGBwcHDsxLRERkZKz0ExERUbWFhoYiNDS0zPccmJeIiMi42LyfiIiIiIiIyEyx0k9EBnH06FGEh4cjLCwMISEh+OmnnwAAN27cwJAhQ9CzZ08MGTIEN2/eNG6gRERERERmhM37iUjvBEHA9OnT8dVXX6F169a4fPkyXnvtNXTr1g3z5s3DsGHDEBYWht27d2Pu3LnYtm2bsUMmIiIiMmsN7G1hY/2/6qC610oWFBbh8aM8Q4RFesBKPxEZhFQqRXZ28euKsrOz4eLigqysLFy8eBExMTEAgODgYCxatAiZmZkc3ZuIiIhIj2ysLTFtzS8AKn+l5sopnQwVFukBK/1EpHcSiQSrV6/G+PHjYWdnh9zcXGzatAnJyclwdXWFhYUFAMDCwgIuLi5ITk5mpZ+IiIiISAdY6ScivSsqKsKmTZuwfv16yGQy/Pnnn5g6dSpWrFhR7XU3bFhPBxGaFnXN68yNlZWl2s+lmfN+Mee0EZH5OHr0KNasWQNBEKBUKjFp0iT06NEDN27cQGRkJB4+fAgHBwdERUWhZcuWxg6XiP4fK/1EpHeXLl1CWloaZDIZAEAmk8HW1hY2NjZITU2FQqGAhYUFFAoF0tLS4ObmpvG6MzJyoFQK+grd4Jyd6+PBg2y9rf/p/nuV0WcfPmfn+ipNCStrWqjP/WJM+vjNpVKJWd4QIyLj4fg8RDUXK/1EpHeNGzdGSkoKrl+/jmeffRZJSUlIT09HixYt0LZtW8THxyMsLAzx8fFo27Ytm/brUen+e5pgHz4iIirB8XmIaiZW+olI75ydnTF//nxMmTIFEokEALBs2TI4ODhg/vz5iIyMxPr169GgQQNERUUZOVoiIiJ6GsfnIaq5WOknIoMIDQ1FaGhome9btWqFHTt2GCEiIiIi0pQ+x+cBjDdGj6mMqWKsOEqPpaNuXB1A8xgrW09l66ztv4k+sNJPRERERAC0H/eDag99js8DGGeMHn2Po2PqcZQeW6eycXUAzcbWeXq8nso8vc7a/ps8HYOusFQnIiIiIgDaj/vxNI4DYr44Pg8ZkzY3JPU5CHFNZZBKf1RUFA4ePIh79+5h7969aN26NQCofb0HX/1BRERERGQaOD4PGZM2NyR587Esg1T6u3btitdffx3Dhw9X+V7d6z346g8iIiIiItPB8XmIaiaDVPq9vb3LfJeRkVHh6z0EQeCrP4iIiEgjuuiHzuagRERkrozWp1/d6z0EQdDZqz9KjwRqTiMwPs0YadNmZM4SVYmTvxsREalT3X7oAJuDEhGR+TL7gfxKRgI1hREY9cUYadN2ZM4S2sbJ3033pFKJ0V6LQ0REREREhmW0Sr+bm1uFr/cQBEEnr/4gIiIiIiIiqs2MVulv2LCh2td78NUfVBVV6dfJfpxERERERGSuDFLpX7x4MX766Sekp6dj1KhRcHBwwL59+9S+3oOv/qCqqEq/TvbjJCIiIiIic2WQSv/s2bMxe/bsMt+re70HX/1BREREREREVD1SYwdARERERERERPrBSj8RERERERGRmWKln4iIiIiIiMhMsdJPREREREREZKZY6SciIiIiIiIyUwYZvZ+IiIiIiIhqJnmREs7O9Y0dBlWRRpX+c+fO4aWXXirz/fnz59G+fXudB0VERERERESmwcpSimlrfql0vpVTOhkgGtKWRs37R40aVe73b7/9tk6DISIiIiIiIiLdUfukX6lUQhAElb8St2/fhoWFhd4DJCIi49GmOV9BYREeP8rTc0RkqgoKCrB06VL8/vvvsLGxQYcOHbBo0SLcuHEDkZGRePjwIRwcHBAVFYWWLVsaO1yTVN3ms8yDRERUHrWV/ueffx4SiUT8vzSpVIqxY8fqLzIiIjI6TZvzAWzSV9utXLkSNjY2OHjwICQSCdLT0wEA8+bNw7BhwxAWFobdu3dj7ty52LZtm5GjNU3a5LfyLJ0QyD63RERUhtpK/+HDhyEIAkaOHInt27eL30skEjg5OaFOnTp6D5CIiIhMW25uLuLi4vDLL7+IDwsaNWqEjIwMXLx4ETExMQCA4OBgLFq0CJmZmXBycjJmyGapujcNAN68IyIyR2or/U2aNAEAHD161CDBEBERUc1z584dODg44NNPP8WpU6dQt25dTJkyBXXq1IGrq6vYHdDCwgIuLi5ITk7WqtLfsGG9Suexsqr+C4mq+5Rcmxgqmre66dDFftBnDJqu2xxbLJhjmoioZtCo5H348CE+//xzXLp0CU+ePFGZ9tVXX+klMCIiqln4Op/aq6ioCHfu3MHzzz+PGTNm4Ny5cxg7dizWrFmjk/VnZORAqRQqnO7sXB9yeVG1t/PgQXaVl9UmBisrywrnrW46dLEf9BWDunQ/rTq/hSlydq6v8zRJpRKNbogREWlU6X///fdRWFiI3r17w9bWVt8xERFRDaRt02I2IzYf7u7usLS0RHBwMADgpZdegqOjI+rUqYPU1FQoFApYWFhAoVAgLS0Nbm5uRo6YiIio9tCo0p+YmIg//vgD1tbW+o6HiIiIahgnJyf4+fnh5MmTCAgIwI0bN5CRkYGWLVuibdu2iI+PR1hYGOLj49G2bVv25yciIjIgjSr9Hh4eSElJQfPmzfUdDxGZKb7Oi8i8LViwADNnzkRUVBQsLS2xYsUKNGjQAPPnz0dkZCTWr1+PBg0aICoqytihEhER1SoaVfo7duyIt99+G/3790ejRo1Upg0cOFAvgRGReeHrvIjMW7NmzfDll1+W+b5Vq1bYsWOHESIiIiIiQMNKf0JCAlxdXXHy5EmV7yUSiU4q/UePHsWaNWsgCAKUSiUmTZqEHj168AkgkZng67yIiIiIiIxDo0p/eXfudUUQBEyfPh1fffUVWrdujcuXL+O1115Dt27d+ASQyEzo+3VeREREpH/sqkdUM2lU6VcqlRVOk0ql1Q5CKpUiO7v4NSbZ2dlwcXFBVlYWnwASmQl9vs7LHF9XpO/X3mn7Dm5dvHtc0/nVLW/OrwM057QRkflgVz2imkmjq7Pnn39ebJL7tEuXLlUrAIlEgtWrV2P8+PGws7NDbm4uNm3ahOTkZJ08ASxdITDniypjpE3bi3uganFqu4yh4tIFcz4mS9Pn67wqe393TaOPdzk/vX5t38GtzfzVWXdl7/A2t/d2l+D7u4moJmBXPaKaS6Pa0eHDh1U+P3jwAJs3b0aXLl2qHUBRURE2bdqE9evXQyaT4c8//8TUqVOxYsWKaq8b+F+FQN8X0sZkjLRVpeIAaH/Rrm3aDBWXLhjrmDRGZYCv8yIiIqrZ2FWPqObSqNLfpEmTMp+joqIwcOBADBo0qFoBXLp0CWlpaZDJZAAAmUwGW1tb2NjYVPsJIBGZDr7Oi4iIqObSZ1c9wHjd9Uyl1aWx4ijdQray1rKatqbVptVteemuaF9Ud73aMpVjQxe0bwf9/3JycpCZmVntABo3boyUlBRcv34dzz77LJKSkpCeno4WLVrwCSCRGeHrvIiIiGoufXbVA4zTXc9UWgIbK47SLWQr62IHaN6FT5tWt0+nu6J9oW1r3uruT1M4NnR500GjSv+0adNU+vTn5+fjzJkzCA0NrXYAzs7OmD9/PqZMmSJuY9myZXBwcOATQCIiIiIiE8CuekQ1l0aV/hYtWqh8trW1xdChQ/HKK6/oJIjQ0NBybyDwCSARERERkWlgVz2imkmjSv/EiRP1HQcREREREZkwdtUjqpk07tO/c+dO7N69G6mpqXB1dUVYWBgGDBigz9iIiIiIiIiIqBo0qvRv2LABcXFxiIiIgLu7O+7fv48tW7YgLS0N48aN03eMRERERERERFQFGlX6d+zYgS+//FLl1X0BAQEYMWIEK/1EREREREREJkqqyUx5eXllRuB0cHBAfn6+XoIiIiIiIiIiourTqNIfGBiIDz74ANevX0d+fj6SkpIQGRmJgIAAfcdHRERERERERFWkUfP+uXPnYuHChQgLC0NRUREsLS3Rq1cvzJkzR9/xkQE0sLeFjbXGYzoSERERERFRDaFRTa9evXpYsWIFli9fjqysLDg6OkIq1aiRANUANtaWmLbmF62WWTmlk56iISIiIiIiIl1RW3O/evUqPvvss//NLJWiYcOGkEql+Oyzz5CUlKT3AImIiIiIiIioatRW+qOjo+Hm5lbutCZNmiA6OlovQRERERERERFR9alt3n/27FlERUWVO61bt24VTiMiItIleZESzs71NZq3oLAIjx/l6TkiIiIioppBbaX/0aNHFfbdl0gkePz4sV6CIiIiKs3KUqrx2CMcc4SIiIjof9Q272/atCkSExPLnZaYmIgmTZroJSgiIiIiIiIiqj61lf5BgwZh9uzZuHDhgsr3//zzD+bMmYMhQ4boNTgiIiIiIiIiqjq1zftff/113L59G4MHD0bjxo3h4uKCtLQ0pKam4rXXXsPIkSMNFSdRrdTA3hY21hq9WVPE/sxERERERFSi0trE7NmzMXLkSPz+++94+PAhHBwc8PLLL6NFixaGiI+oVrOxttS4H3MJ9mcmIiIyrvJu2ms6GGkJ3sQnIl3R6BFiixYtWMknIiIiItLA0zftrawsIZcXabUO3sQnIl3Rrt2wnhQUFGDp0qX4/fffYWNjgw4dOmDRokW4ceMGIiMjxRYGUVFRaNmypbHDJSIiIiI90eYVneXhE3LzpGmXR/7+RGWZRKV/5cqVsLGxwcGDByGRSJCeng4AmDdvHoYNG4awsDDs3r0bc+fOxbZt24wcLRERERHpizav6CwPn5CbJ027PPL3JypL7ej9hpCbm4u4uDhMmTIFEokEANCoUSNkZGTg4sWLCA4OBgAEBwfj4sWLyMzMNGa4REREetHA3hbOzvU1+mtgb2vscCv06aefwsPDA//++y8A4MaNGxgyZAh69uyJIUOG4ObNm8YNkIiIqJYx+pP+O3fuwMHBAZ9++ilOnTqFunXrYsqUKahTpw5cXV1hYWEBALCwsICLiwuSk5Ph5OSk8fobNqwn/l+dpmKmrrpps7LS/lCoyjJViVPbZbSNqyrNCOVFSlhZVv+emSbbNdR+JiLj0mbgTlN9kvXPP//g7NmzcHd3F79jqz0iIiLjqrA20alTJ/HJuzrHjh2rVgBFRUW4c+cOnn/+ecyYMQPnzp3D2LFjsWbNmmqtt0RGRg6USgHOzvXx4EG2TtZpaqqbNmfn+loPLgOgSstoG6e2aatKWqrSjHDllE7VPp40SVtVfxt165VKJSo3w4iIdKGwsBALFy7ERx99hDfeeAMAxFZ7MTExAIpb7S1atAiZmZla3cAnIiKiqquw0r9y5Urx/7///htxcXEYOXIk3N3dcf/+fWzfvh3h4eHVDsDd3R2WlpZiM/6XXnoJjo6OqFOnDlJTU6FQKGBhYQGFQoG0tDS4ublVe5tERESkW2vWrEFoaCiaNWsmfpecnKyTVntERERUdRVW+n19fcX/Fy5ciK1bt8LV1VX87tVXX8Xbb7+NiIiIagXg5OQEPz8/nDx5EgEBAbhx4wYyMjLQsmVLtG3bFvHx8QgLC0N8fDzatm3LiwQiIqoRNB1pGqj5XXISExPx999/44MPPtDL+jVpnVSVrlBPM2RXuYrmrW46dLEf9BmDpuuubgymkKeeTgO76xGRsWhU+qSlpcHOzk7lOzs7O6SmpuokiAULFmDmzJmIioqCpaUlVqxYgQYNGmD+/PmIjIzE+vXr0aBBA0RFRelke0RkHJ9++inWrVuHvXv3onXr1nwtJ5k1Tfvol7y/21T76WvizJkzuH79Orp27QoASElJwVtvvYUPP/xQJ632SrrqVaSqXaGeZqiucure2V7ddOhiP+grBm3eVV/dGIzdpfPp40GbtJfG7npEpAsaVfqDgoIwbtw4jBs3Do0bN0ZycjI2bdqEoKAgnQTRrFkzfPnll2W+b9WqFXbs2KGTbRCRcXGALyLzNXr0aIwePVr8HBQUhI0bN6J169b45ptvakSrveq+G56IiMhUaVTpX7BgAdatW4d58+YhLS0NLi4u6NWrFyZOnKjv+IjIDHCAL6Laq6a02uO74YmIyFxVWulXKBSYO3cuFi1apLe+ekRk3jjAF1HtcuTIEfF/ttojMj/srmd4mo4TU1BYhMeP8gwQEdUklR45FhYWOHnypEav7yMiepopDPBV0+i7ibG2g0npYnAyTefXxQBggOkMfqXtoGU1MY1EVLuwu55xaDpODFsdUXk0urp44403sG7dOkycOBHW1tb6jomIzIixB/iqaZyd6+t1AKqqDHimzfzVWXdlA11ps25jD+IFaL6vS6dbV2nkAF9EpA/srkdUM2lU6d++fTvS09MRExMDJycnlaf+x44d01dsRGQGzGGALyIiItJvdz1dvppTm5ZQptJqSpM49JH+0uusbP36eOVmebFWFL+m69V0YFZ5kRJWllKtYqupNNpzK1eu1HccRFQL1ZQBvnRB23e2s08eERGZEn1319Plqzk1be2l79Z1mtIkDn2lv2SdmrxWUh+v3Hw61or2hVavRdVwYNaVUzpVuK9M4djQ5U0Hja5AfX19dbZBIqrdausAX7Xpne1ERGR+9N1dj4j0R+O2F5cuXUJCQgKysrIgCP+7CzdlyhS9BEZERERERKaB3fWIai6NKv3fffcdli1bBn9/fxw/fhyvvvoqTp48Kd7pIyIiIiKi2qk2ddcjqok0qvRv2bIFW7Zsgbe3N3x8fBAdHY1ffvkF+/fv13d8RERERERkYmprdz2imqji4QpLycjIgLe3d/ECUimUSiU6deqEo0eP6jU4IiIiIiIiIqo6jZ70N27cGHfv3kXTpk3RsmVLHD58GI6OjrCystJ3fERERERERERURRpV+t9++20kJSWhadOmGD9+PKZMmQK5XI5Zs2bpOz4iIiIiIiIiqiK1lf79+/fDx8cH/fv3F7/r1KkTTp8+Dblcjrp16+o9QCIiIiIiIqqcvEip0/e7k3lQW+lfs2YNbt++jebNm4uD+Pn4+KBJkyawtrY2VIxkJqpSCMmLlHqKhoiIiIjIvFhZSjFtzS8azbtySic9R1M7NbC3hY115Q3qCwqL8PhRngEiqqTSf/DgQaSnp+PMmTNISEhATEwMZs6cCVdXV3h7e8PX1xeDBg0ySKBU82lTCJVgYURERERERDWFjbWlRnUeQ9ZzKr0F0ahRI/Tu3Ru9e/cGADx+/BjfffcdvvjiC8THx7PST0RERERERGSiKq30C4KAS5cuiU/7ExMT4eLigt69e0Mmk+k0mE8//RTr1q3D3r170bp1a9y4cQORkZF4+PAhHBwcEBUVhZYtW+p0m0RERERERETmSm2lf8yYMfjnn3/wzDPPQCaTYfDgwVi2bBnq1aun80D++ecfnD17Fu7u7uJ38+bNw7BhwxAWFobdu3dj7ty52LZtm863TURERERERGSOpOom3rhxA9bW1mjatCmaN2+OFi1a6KXCX1hYiIULF2LevHmQSCQAgIyMDFy8eBHBwcEAgODgYFy8eBGZmZk63z4RERERERGROVL7pP+nn35SGcjvv//9L7KysuDl5QVvb2/IZDK0bdu22kGsWbMGoaGhaNasmfhdcnIyXF1dYWFhAQCwsLCAi4sLkpOT4eTkpPG6Gzb8300Kc359RXXTZmVV+QiTxlpG27QZIq6qvonAylL1Ppsm6zDEPiMiIiIiIvNU5YH8NmzYgMzMTFy6dKlaASQmJuLvv//GBx98UK31VCQjIwdKpQBn5/p48CBbL9swtuqmzdm5PuTyIq2XM9Qy2qTNUGmp6psISqdFk9+tqulRt16pVKJyM4yIiMic6OI95YZ8lRYRkb5pPZDfn3/+icePH+OFF17AgAEDqh3AmTNncP36dXTt2hUAkJKSgrfeegsffvghUlNToVAoYGFhAYVCgbS0NLi5uVV7m0RERERknqpyY/5pfGUwEZkTtZX+0aNHIzExEXK5HO3bt4evry+GDx8OT09P2NjY6CSA0aNHY/To0eLnoKAgbNy4Ea1bt8Y333yD+Ph4hIWFIT4+Hm3bttWqaT8RERERERFRbaa20u/t7Y2xY8fixRdfhJWVlaFiEs2fPx+RkZFYv349GjRogKioKIPHQERERERERFRTVfqk39COHDki/t+qVSvs2LHD4DEQERERERERmQO1r+wjIiIiIiIiopqLlX4iIiIiIiIiM6X9C8DNTAN7W9hYa7cb+BoXIiIiIiIiqglqfaXfxtqySu9bJyIiIiIiIjJ1bN5PREREREREZKZq/ZN+Q2E3AiIiMldZWVmYPn06bt++DWtra7Ro0QILFy6Ek5MTbty4gcjISDx8+BAODg6IiopCy5YtjR0yEZHRla4fODvXN3I0ZM5Y6TcQdiMgIiJzJZFI8Pbbb8PPzw8AEBUVhY8++ghLly7FvHnzMGzYMISFhWH37t2YO3cutm3bZuSIiYiMr6R+YGVlCbm8SO28rBdQdbB5PxHpXVZWFt555x307NkTISEhmDhxIjIzMwEAN27cwJAhQ9CzZ08MGTIEN2/eNG6wRKQ1BwcHscIPAB06dMD9+/eRkZGBixcvIjg4GAAQHByMixcvivmfiIiI9I+VfiLSu5KngAcPHsTevXvRrFkzfPTRRwAgPgU8ePAghg0bhrlz5xo52pqngb0tnJ3ra/RHpG9KpRLffPMNgoKCkJycDFdXV1hYWAAALCws4OLiguTkZCNHSUREVHuweT8R6V15TwG/+eYb8SlgTEwMgOKngIsWLUJmZiacnJyMFa5JkBcptaqka9p9iM0DSd8WLVoEOzs7jBgxAhcvXtTJOhs2rFfpPFZW1b+kqe46tFm+onkNGYO+1qFueU3XbQr7obo3Sp+OoSoxmdLNWo7dQVRzsdJPRAal6VPA2l7pt7KUsiJPNU5UVBRu3bqFjRs3QiqVws3NDampqVAoFLCwsIBCoUBaWhrc3Ny0Wm9GRg6USqHC6c7O9SvtD6uJ6q5D0+XV9d81VAz6XEdFy2vSb1nfMWjjwYPsKi/79DGpTdo1jUEqlWh0Q0xXOHYHUc3FSj8RGZSunwIa8oKnurR9wqWLp4bVndfQ69bFE0JtW0nIi5SwstRPbzd9/uam9AQQAFatWoULFy5g8+bNsLa2BgA0bNgQbdu2RXx8PMLCwhAfH4+2bdvW+pt6RDURW+0R1Vys9BORwejjKWBlTwBNhaZPIks/DdLmqZC+5jXkuit7EqbxU1QtWkkAxS0lqvNEryL6/s1N6Qng1atXsXHjRrRs2RJDhw4FADRt2hTR0dGYP38+IiMjsX79ejRo0ABRUVEGi4uI9EMfrfZ02Y1Hm5uims6r6Q1ibW8ka3PTVx/dY0rPW9ly+th+efu/ot9EHw8h1P3+1bm5ro9jtTpY6Scig+BTQCLz9dxzz+HKlSvlTmvVqhV27Nhh4IiISJ/0MXaHLrvxaHoj19m5vlbzanJDWZsbySVp0rT7hz66x5TMq0kM+tj+0/uqot9E225c1T1WtDk2yltWF8eqLm8IsNJPRHrHp4BERETmQV9jdxCR/rDST0R6x6eAZMq0GQOgoLAIjx/l6TkiIiLTxFZ7RDUTK/1ERFSraTMGwNIJgSY3gB4RkSGw1R5RzWX0Sj/f+UlERDUFX6VIRLUVW+2RrlXU0o4313XP6JV+vvOzYto0OS2Zj01PiYiIiIjI1JV3I72iAQV5I716jF7p5zs/K6bpE6XSmYMZgoiIiIiIiEoYvdJfmr7f+amLdz5Wti51qrKdqrwP09Riq+4y2qanJqVFk7QZ6vgkIiIiIiLzY1KVfn2+81NX73wsoe17G6u6HU2WeboZjCnFpotltElPTUqLJu//1MfxKZVKVG6GERERERGR+TKZSj/f+UlEREREpkCbcZWIiEydSVT6+c5PItKXBva2sLHWvKjjYJhERKTNmzrKwzGWiMiUGL3Sz3d+EpE+2VhbanXhxgs1IiIiIjInRq/0852fRERERERERPohNXYARERERERERKQfrPQTERERERERmSlW+omIiIiIiIjMlNH79BMRmRJtXtPEkf6JiIiIzJ+2b4MyNTU3ciIiPdDmNU0c6Z+IiIjI/GnzNihTvD5k834iIiIiIiIiM8VKPxEREREREZGZYqWfiIiIiIiIyEyx0k9ERERERERkpljpJyIiIiIiIjJTrPQTERERERERmSm+so/IQMp7/7um74MnIiIiIiKqClb6iQzk6fe/W1lZQi4vUruMKb7nk4iIiIiIag427yciIiIiIiIyU6z0ExEREREREZkpVvqJiIiIiIiIzBQr/URERERERERmyuQr/Tdu3MCQIUPQs2dPDBkyBDdv3jR2SESkY8znROaNeZzI/DGfE5kuk6/0z5s3D8OGDcPBgwcxbNgwzJ0719ghEZGOMZ8TmTfmcSLzx3xOZLpM+pV9GRkZuHjxImJiYgAAwcHBWLRoETIzM+Hk5KTROqRSSbn/l+ZY30br2CpalzpV2Y4my1haWaJIbiF+NqXYdLGMtukx5bSUXubp302X21G3z6pyfOhTdfO5JunRdh9qM782+1Pb/KxNHPqa15DrrixP1MQ06vM3ryn5XNfn8opUpazU9To0XV7dsW6oGPS5joqW1/S8p88YDLmOqpzzn8Z8rkrT30SbfaPr83hV1qnra8KqnKc0iUGf2y+hq/JRF7/V09MMvX1dkgiCIBhkS1Vw4cIFzJgxA/v27RO/69OnD1auXIl27doZMTIi0hXmcyLzxjxOZP6Yz4lMm8k37yciIiIiIiKiqjHpSr+bmxtSU1OhUCgAAAqFAmlpaXBzczNyZESkK8znROaNeZzI/DGfE5k2k670N2zYEG3btkV8fDwAID4+Hm3bttW4bxARmT7mcyLzxjxOZP6Yz4lMm0n36QeApKQkREZG4vHjx2jQoAGioqLw7LPPGjssItIh5nMi88Y8TmT+mM+JTJfJV/qJiIiIiIiIqGpMunk/EREREREREVUdK/1EREREREREZoqVfiIiIiIiIiIzxUo/ERERERERkZmyNHYAunL37l1MmDBB/JydnY2cnBycPn1aZb7o6Gjs378fFhYWsLS0xNSpUxEYGGjocLWiadpKXL9+Hf369cOwYcMwY8YMQ4VZJdqkbf/+/diwYQMEQYBEIkFMTAwaNWpkyHC1omnaMjIy8OGHHyI5ORlyuRwdO3bE7NmzYWlpNtmTSjHnY74ytT1PHD16FGvWrIEgCFAqlZg0aRJ69OihMo9CocDixYvx66+/QiKRYPTo0Rg0aJCRIqaqMufzdmVqaxlX28s3U5KVlYXp06fj9u3bsLa2RosWLbBw4cIyrw9ct24dvv76a7i4uAAAvLy8MG/ePADVL4s1jWH69Om4cuWK+PnKlSuIjo5G165d1canjfHjx+Pu3buQSqWws7PDnDlz0LZtW5V51KVXV+clTeJQV0/Txf7QJAZ9HhfaxKHzY0MwU4sXLxYWLFhQ5vvjx48LT548EQRBEC5duiTIZDIhLy/P0OFVS0VpEwRBKCoqEkaMGCG89957wvLlyw0cWfVVlLbz588LvXv3FtLS0gRBEITHjx8L+fn5hg6vWipK2+LFi8XfqrCwUBg4cKCwb98+Q4dHRmLOx3xlalOeUCqVgre3t3DlyhVBEIrPPx06dBAUCoXKfLGxsUJERISgUCiEjIwMITAwULhz544xQiYdMufzdmVqaxlXm8o3U5OVlSX88ccf4ufly5cLH374YZn51q5dW2Geq25ZrGkMpV26dEnw9fUVCgoKKo1PG48fPxb///nnn4Xw8PAy86hLr67OS5rEoa6epov9oUkM+jwutImjNF0cG2bZvL+wsBB79+7FgAEDykwLDAyEra0tAMDDwwOCIODhw4cGjrDq1KUNADZv3ozOnTujZcuWhg1MB9Sl7YsvvkBERAScnZ0BAPXr14eNjY2hQ6wydWmTSCTIzc2FUqlEYWEh5HI5XF1djRAlGZo5H/OVqY15QiqVIjs7G0DxU0AXFxdIpaqn4f3792PQoEGQSqVwcnJCt27d8OOPPxojXNIRcz5vV6a2lnG1sXwzJQ4ODvDz8xM/d+jQAffv39dqHdUti6sSww8//ICQkBBYW1trFWtl6tevL/6fk5MDiURSZh516dXVeUmTOPRdT9MkBnUMuS9K08WxYZZtiY4cOQJXV1e0a9dO7XxxcXFo3rw5GjdubKDIqk9d2i5fvowTJ05g27ZtWL9+vRGiqx51aUtKSkLTpk0xfPhwPHnyBN27d8e4ceO0zqzGoi5t48ePx6RJkxAQEIC8vDwMHz4cMpnMCFGSoZnzMV+Z2pYnJBIJVq9ejfHjx8POzg65ubnYtGlTmfmSk5Ph7u4ufnZzc0NKSoohQyUdM+fzdmVqaxlX28o3U6ZUKvHNN98gKCio3On79u3DiRMn4OzsjEmTJsHT0xOAbsviymIA/nej6IsvvtAoPm3NmjULJ0+ehCAI2LJlS5np6tKry31RWRyllVdP08X+0CQGQxwXmu4LXR0bZvmkf+fOnRXeUS9x+vRprFmzBh9//LGBotKNitIml8sxZ84cLFiwABYWFkaIrPrU/W4KhQJXrlxBTEwMvvzySxw/fhy7d+82cIRVpy5tP/74Izw8PHDixAkcP34cCQkJfLJXS5jzMV+Z2pYnioqKsGnTJqxfvx5Hjx7Fhg0bMHXqVOTm5ho7NNIzcz5vV6a2lnG1rXwzZYsWLYKdnR1GjBhRZtrQoUNx+PBh7N27F2+99RbGjx+PrKwsg8ZQ4tChQ3B3d1fp163L+JYsWYJjx45h6tSpWLFiRZXWoQuaxlFePU1X+6OyGAx1XGi6L3R1bJhdpT81NRVnzpxBSEhIhfMkJiZi2rRpiI6OxrPPPmvA6KpHXdoePHiA27dvY/To0QgKCsJ///tffP/995gzZ44RItVeZb+bu7s7evXqBWtra9SrVw9du3bF+fPnDRxl1VSWtu3btyM0NBRSqRT169dHUFAQTp06ZeAoydDM+ZivTG3ME5cuXUJaWpr4RE8mk8HW1hZJSUkq87m5uak0AU1OTq5RrdFIlTmftytTW8u42li+maqoqCjcunULq1evLtOVCgCcnZ1hZWUFAPD394ebmxuuXr0KQHdlcWUxlCjvRpG6+KoqPDwcp06dKlNBVJdefZyXKooDqLiepuv9UVEMhjguNImjhK6ODbOr9MfGxqJTp05wdHQsd/r58+cxdepUrF27ttLm/6ZGXdrc3d1x6tQpHDlyBEeOHMEbb7yBwYMHY9GiRUaIVHuV/W7BwcE4ceIEBEGAXC7HH3/8gTZt2hg4yqqpLG1NmzbF8ePHARQ34fn999/x3HPPGTJEMgJzPuYrUxvzROPGjZGSkoLr168DKG7anJ6ejubNm6vM16tXL+zYsQNKpRKZmZk4dOgQevbsaYyQSQfM+bxdmdpaxtXG8s0UrVq1ChcuXEB0dHSF/aBTU1PF/y9duoR79+7hmWeeAaCbsliTGAAgJSUFf/75J4KDgzWOT1O5ublITk4WPx85cgT29vZwcHBQmU9denWxLzSNQ109rbr7Q9MY9H1caBoHoNtjw+z69MfGxmLWrFkq373zzjuYPHkyXnzxRSxYsAD5+fmYO3euOH3FihXw8PAwdKhaqyxtNVllaevbty8uXLiAPn36QCqVIiAgAAMHDjRStNqpLG0zZ87EvHnzEBISAoVCAT8/PwwePNhI0ZKhmPMxX5namCecnZ0xf/58TJkyReyzvGzZMjg4OKikPSwsDOfOnRNf5TdhwgQ0a9bMmKFTNZjzebsytbWMq43lm6m5evUqNm7ciJYtW2Lo0KEAim+2REdHq/wWn3zyCf755x9IpVJYWVlhxYoV4sCS1S2LNY0BKD5munTpUqbSpy4+TeXl5WHKlCnIy8uDVCqFvb09Nm7cCIlEovG5RxfnJU3jUFdPq+7+0DQGfR4X2sQB6PbYkAiCIGgVKRERERERERHVCGbXvJ+IiIiIiIiIirHST0RERERERGSmWOknIiIiIiIiMlOs9BMRERERERGZKVb6iYiIiIiIiMwUK/0m6v79+/D09IRCoah03lWrVsHPzw/+/v46jWHPnj2IiIjQ6TqJ6H9GjhyJHTt2VGsdp06dwquvvlrh9Llz5yI6OlqjeYnItGhzLUBENVffvn1x6tQpY4dBZoyVfh1KSEjA0KFDIZPJ4Ovri6FDh+L8+fMaLRsUFITffvtN/Ozu7o7ExERYWFioXS45ORkxMTHYv38/Tp48WeXY7969Cw8PDxQVFYnfhYaG4vPPP6/yOity6tQpeHh4YOLEiSrfX758GR4eHhg5cqTOt1mZlJQUTJo0CX5+fpDJZAgJCcGuXbsMHgcZz6ZNm/DOO++ofNejR49yv9u3b5/W61+3bh0++OCDasVYFQsXLsSECRMMvt3KsBygmuLp8zMA7Nq1C6+99pret63ptYAhlFwn9OvXT+X7zMxMvPDCCwgKCjJ4TI8fP8aHH34If39/eHp6omfPnti8ebPB46DaZeTIkfDx8UFhYaHO1rlv3z74+flVOl951+va2LVrFzw8PLBs2TKV7w8dOgQPDw9ERkZWab3VcfXqVURERMDHxwfe3t7o378/fvnlF4PHYe4sjR2AucjJycHYsWMxf/589O7dG3K5HAkJCbC2ttbrdu/duwcHBwc0bNhQr9vRNScnJyQmJiIrKwuOjo4AgNjYWLRs2dIo8UybNg1t2rTB0aNHYW1tjX///RcPHjzQ6TaKiopgacksZ6q8vb2xefNmKBQKWFhY4MGDBygqKsLFixdVvrt16xa8vb21WndVT87mjuUAUc3z5MkT/Pvvv2jdujUAID4+Hk2aNIFcLjd4LMuWLcOTJ0+wf/9+1K9fHzdu3MDVq1d1ug3mWSrt7t27SEhIQP369XH48GH07t3b2CFprXnz5ti/fz+mTZsmHttxcXFGO/eOHTsWr732GjZu3AgA+PvvvyEIgk63wXzMJ/06c+PGDQBAcHAwLCwsUKdOHQQEBKBNmzYAgNu3b+P111+Hn58f/Pz88P777+Px48cAii8079+/j7Fjx8LT0xOfffZZmTt5u3btQteuXeHp6YmgoCDs2bMHv/32GyIiIpCWlgZPT0/x7tzkyZPh7+8PmUyG4cOHq5wA8/PzsXz5cnTp0gUymQyvvfYa8vPzMWLECACAj48PPD09kZiYWOZJxl9//YUBAwZAJpNhwIAB+Ouvv8RpI0eOxOrVqzF06FB4enoiIiICmZmZFe4vKysrdO3aFfv37wcAKBQKHDhwACEhISrzJSUlYdSoUfD19UXPnj3F+QHg2LFjCA8Ph5eXFzp16oR169aJ00r2X2xsLDp37gw/Pz9s2LChwnguXLiA/v37w87ODpaWlnj++efRqVMncXpJKw5vb2906tRJfPqXnZ2N6dOno2PHjujSpQvWr18PpVIp/mZDhw7F0qVL4evri3Xr1qGwsBBRUVHo3LkzXnnlFcydOxf5+fkVxkWG8+KLL6KoqAiXLl0CAJw5cwZ+fn545plnVL5r3rw5XF1dkZqairFjx8LX1xfdu3fH999/L65r3bp1mDx5Mj744AN4eXnh22+/xaZNm3DgwAF4enoiNDRUnPfevXvl5psZM2aILW1SU1Ph4eGBr776CgBw69Yt+Pr6qpwUP//8c7z88ssICAjAzp07xe8jIyOxatWqctMcFBSErVu3IiQkBDKZDO+++y4KCgrE6UePHkVYWBi8vb0xdOhQXL58WZy2efNmBAYGik/Xfv/9dwDA+fPn0b9/f3h5eeGVV14p8zShNJYDLAfMxebNm9GtWzd4enqiT58++Pnnn8VpJcfAokWLIJPJ0KtXLzG/AMXnz48//hgDBw6ETCbDuHHj8PDhQwBln+pVdq49e/aseIyGhoaqNBcu7zoCKC5PRowYAZlMBj8/P7z77rtq0xoWFobY2Fjxc1xcHMLDw1XmSU1NxaRJk9CxY0cEBQVh27Zt4rTz589jyJAh8Pb2RkBAABYuXKjyxNTDwwPffPMNevToAR8fHyxYsKDCCsDff/+NkJAQ2NvbQyqVolWrVujVq5c4/erVq2LZ8corr4iVisLCQixZsgQBAQEICAjAkiVLxBhKukFt3rwZ/v7++PDDD6FUKsXf2M/PD1OmTBF/I6pd4uLi8NJLL6Ffv36Ii4tTmRYZGYkFCxZg9OjR8PT0xKBBg3D79m0AxdfQfn5+SE5OBlDcqs3b2xtJSUkAVFsUqTvenr5eP336NHx9fXHlyhUxjoyMDLRv377C6/BGjRqhdevWOHHiBADg4cOHSExMLNNaR115snPnTvTu3Ruenp7o2rUrvv32W3FaSR6q6LqktMzMTNy9exeDBw+GtbU1rK2tIZPJVB6uHDp0CGFhYfDy8kK3bt1w/PhxANDqOiw2NhbZ2dmYOXMmAgICEBgYiFWrVtWurlMC6UR2drbg6+srTJ8+XTh27Jjw8OFDlek3b94UTpw4IRQUFAgZGRnCsGHDhMWLF4vTu3TpIpw8eVL8fOfOHaF169aCXC4XcnNzBU9PTyEpKUkQBEFITU0V/v33X0EQBOGPP/4QAgMDVba1Y8cOITs7WygoKBAWL14shIaGitPmz58vjBgxQkhJSRGKioqEP//8UygoKFDZXomdO3cKQ4cOFQRBELKysgRvb28hNjZWkMvlwt69ewVvb28hMzNTEARBGDFihNC1a1fh+vXrQl5enjBixAhh5cqV5e6rkpj//PNPYeDAgYIgCMKxY8eEiIgI4fvvvxdGjBghCIIg5ObmCq+++qrwww8/CHK5XLhw4YLg6+urkvbLly8LCoVCuHTpkvDyyy8LP//8s8r+mzVrlpCXlydcunRJaNeunXDt2rVyY3rjjTeEIUOGCPHx8cK9e/dUpt27d0/o0KGDsHfvXqGwsFDIzMwULl68KAiCIEybNk0YO3askJ2dLdy5c0fo0aOH8P3334v7r23btsK2bdsEuVwu5OXlCYsXLxbGjBkjZGVlCdnZ2cKYMWOEjz76qNyYyPBGjBghxMTECIIgCAsWLBB27NghfPLJJyrfRUZGCoIgCMOHDxfmzZsn5OfnCxcvXhT8/PyE3377TRAEQVi7dq3w/PPPCz///LOgUCiEvLw8Ye3atcL7779fZnsV5ZsdO3YIY8aMEQRBEPbs2SN07dpVmDJlijht7NixgiAU54O2bdsKq1evFgoLC4Vjx44J7du3F8ugGTNmCJ988ok4b+nyokuXLsKAAQOElJQUISsrS+jVq5fw9ddfC4IgCBcuXBA6duwonD17VigqKhJ27doldOnSRSgoKBCSkpKEV199VUhJSREEoTi/3bp1SxAEQRg8eLAQGxsrCIIg5OTkCImJieXua5YDLAdqiqfPz4Kgen4UBEHYv3+/kJKSIigUCmHfvn3CSy+9JKSmporztm3bVoiJiREKCwuFffv2CV5eXkJWVpYgCMXlQEBAgHDlyhUhNzdXmDhxolhWPH1uVldmpKSkCL6+vsKxY8cEhUIhnDhxQvD19RUyMjLUXkdMnTpVWL9+vaBQKIT8/HzhzJkz5e6Hklju3LkjvPrqq0JRUZFw7do1oUePHsLJkyeFLl26CIIgCAqFQujXr5+wbt06oaCgQLh9+7YQFBQkHD9+XBAEQfj777+FxMREQS6XC3fu3BF69eollrGCIAitW7cWRo8eLTx69Ei4d++e4OfnJ/zyyy/lxjRz5kyhT58+wg8//CDcuHFDZVp2drbg7+8vbN26VcjPzxeys7OFs2fPCoIgCKtXrxYGDRokpKenCxkZGcKQIUOEVatWCYLwvzJ1xYoVQkFBgZCXlyfExMQIgwYNEpKTk4WCggJhzpw5wtSpU8uNicxbt27dhO3btwt///238PzzzwsPHjwQp82YMUPw8fERzp07J8jlcuG9994T3n33XXH6J598IowcOVLIy8sTgoODhS+//FKcVrqcUXe8lXe9Pm/ePGHFihXi5y+++EK8fnhaSdm1Z88e8Zpi+/btwpw5c4RPPvlEmDFjhiAI6ssTQRCEo0ePCrdu3RKUSqVw6tQpoX379sKFCxcEQaj8uqQ0pVIpdO/eXRg9erTw888/q+xPQRCEc+fOCV5eXsKJEycEhUIhpKSkiOdwba/Dxo0bJ8yZM0fIzc0V0tPThQEDBgjffPNNhb+1ueGTfh2pV68evv76a0gkEsyZMwcvv/wyxo4di/T0dABAixYt4O/vD2trazg5OWHUqFE4c+aMxuuXSqW4evUq8vPz4eLigueee67CeQcOHIh69erB2toakyZNwuXLl5GdnQ2lUomdO3di1qxZcHV1hYWFBby8vDTqgnDs2DG0aNEC4eHhsLS0RHBwMJ599lkcPXpUnKd///545plnUKdOHfTq1Ut8OloRLy8vPHr0CNevX0dcXBzCwsLKbLNJkyYYMGAALC0t0a5dO/Ts2RMHDx4EAPj5+cHDwwNSqRRt2rRB3759cfr0aZV1TJw4EXXq1EGbNm3Qpk0blSeVpa1Zswbe3t5Yv349unbtirCwMHE8hr179+KVV15BcHAwrKys4OjoiLZt20KhUGD//v14//33Ua9ePTRt2hSjRo0Sn54AgIuLC0aOHAlLS0vY2Nhgx44dmDlzJhwcHFCvXj2MGTOmSv3DST98fX3FfJmQkABvb2/IZDKV73x9fZGcnIw///wTH3zwAWxsbNC2bVsMGjQIu3fvFtfVoUMHdOvWDVKpFHXq1KlwmxXlG19fXyQkJECpVOLMmTN4++23xdY1Z86cga+vr7gOS0tLTJgwAVZWVujUqRPs7OzE1keVGTlyJFxdXeHg4IAuXbqI2//+++8xZMgQvPTSS7CwsEC/fv1gZWWFs2fPwsLCAoWFhUhKSoJcLkfTpk3RvHlzMZbbt28jMzMTdevWRYcOHdRun+UAy4GaYMKECfD29hb/FixYoDK9d+/ecHV1hVQqRZ8+fdCiRQuVMX2cnJzwxhtvwMrKCn369MEzzzyDY8eOidPDwsLQunVr2NnZYcqUKfjxxx8rfAJVUZmxe/duvPrqq+jUqROkUin8/f3xwgsviH1jK7qOsLS0xP3795GWlgYbG5tKuy81btwYzzzzDH777TfExsaWecr/999/IzMzExMnToS1tTWaNWuGwYMHiy10XnjhBXTo0AGWlpZo2rQphgwZUuZ66J133kGDBg3g7u4OPz+/CvPsnDlzEBISgq+++gp9+/ZF9+7dxfQeO3YMjRo1QkREBGxsbFCvXj289NJLAIrz84QJE9CwYUM4OTlhwoQJKnlWKpVi8uTJsLa2Rp06dfDdd99h6tSpaNy4MaytrTFx4kQcPHiQXbdqmYSEBNy/fx+9e/fGCy+8gGbNmiE+Pl5lnu7du6N9+/awtLREaGioyrXwxIkTkZOTg0GDBsHFxQXDhw8vdzvaHm/9+vVDfHy82MJs9+7dKi0Ky9O9e3ecPn0a2dnZ2L17d5lzb2XlSefOndG8eXNIJBL4+vrC398fCQkJ4vKaXpdIJBJs27YNTZo0wfLlyxEQEIDhw4fj5s2bAIAffvgBAwYMgL+/P6RSKVxdXdGqVSutr8NycnJw/PhxzJw5E3Z2dmjYsCHefPPNWnXurd2dG3SsVatWWL58OYDi5qjTpk3D0qVL8cknnyAjIwOLFy9GQkICcnNzIQgCGjRooNF67ezssGrVKnz++eeYNWsWvLy8MGPGDLRq1arMvAqFAqtWrcKPP/6IzMxMSKXF93WysrJQWFiIgoICNGvWTOu0paWlwd3dXeU7d3d3pKamip+dnZ3F/21tbfHkyZNK1xsaGoqvvvoKp06dwtKlS1UKz3v37uH8+fMqFyAKhUIsyM6dO4ePPvoIV69ehVwuR2FhoUqzPqC4CZMmMdnb2+ODDz7ABx98gMzMTKxYsQITJkzA8ePHkZycLFZoSsvKyoJcLlfZL0/vk8aNG4v/Z2ZmIi8vD/379xe/EwRBLKTJ+Ly9vfHVV1/h0aNHyMzMRMuWLdGoUSNERkbi0aNHuHr1Kry9vZGWlgZ7e3vUq1dPXNbd3R0XLlwQP5f+7dWpKN80b94cdnZ2uHTpEv78809MmDABP/zwA65fv44zZ86oDHTn4OCg0ldN0/xX3vbT0tIAFI8aHhcXh+3bt4vT5XI50tLS4Ovri5kzZ2LdunW4du0aAgICEBkZCVdXVyxZsgRr165F79690bRpU0ycOBFdunRRGwPLAZYDpi46OhqvvPKK+HnXrl0qb96Ii4tDTEwM7t27B6C433tWVpY43dXVFRKJRPzs7u4u5jUAcHNzU5kml8tVli+tojLj/v37+PHHH1VuxhcVFcHPz0/tdcS0adOwZs0aDBw4EPb29hg1ahQGDhyodn+Eh4cjNjYWiYmJ2L59O27duiVOu3fvHtLS0srk2ZLPN27cwPLly3HhwgXk5eVBoVCgXbt2atOYm5tbbhx16tTB2LFjMXbsWOTk5GDz5s149913cfTo0QrzLFD2mubp38PR0RE2Njbi5/v372PChAniNRVQfGMgIyMDrq6uavcVmY+4uDj4+/vDyckJQHGX3tjYWLz55pviPKXPN3Xq1FE531hZWaFfv35YvHgxIiMjVcqE0tQdb+V56aWXYGtri9OnT8PZ2Rm3b99G165d1aalTp066NSpE9avX4+srCzIZDKx2XxJDBWVJwDwyy+/IDo6Gjdv3oRSqUR+fr44zgeg3XVJ48aNMXfuXADFA5TPmTMHM2bMwHfffYfk5GSVbnYltL0Ou3//PoqKihAQECB+p1QqVcpec8dKv560atUK/fv3x3fffQcA+PjjjyGRSLBnzx44Ojri0KFDWLhwocbrCwwMRGBgIPLz87F69WrMmTMHX3/9dZn59u7di8OHDyMmJgZNmzZFdnY2fHx8IAiCeBK7c+eOONZAiYoKnhIuLi64f/++ynfJyckIDAzUOA3lCQsLQ48ePRAeHg5bW1uVaW5ubvDx8UFMTEy5y77//vsYMWIEtmzZAhsbGyxZsqTCiyRtODk5ISIiArGxsXj48CHc3NzKfQuDo6MjrKyscP/+ffznP/8BULxPSl8AlN6vjo6OqFOnDvbt28eLBBPl6emJnJwcfPfdd/Dy8gJQ3IrHxcUF3333HVxcXNCsWTNYWlri0aNHyMnJEU846n778j5rwsfHBwcPHoRcLoerqyt8fHywe/duPHr0CG3btq1GSivn5uaGsWPHYty4ceVODwkJQUhICHJycjB37lx89NFHWLlyJVq2bIlPPvkESqUSP/30EyZPnoxTp07Bzs6uwm2xHKCa7N69e5g9eza++OILeHp6wsLCosxTs9TUVAiCIB4LycnJKv1nS/r5lvxf0pqk9PeVcXNzQ1hYGBYvXlzu9IquI5ydncVlEhISMGrUKPj4+KBFixYVbqtHjx5YuHAh2rVrhyZNmqhU+t3c3NC0aVP89NNP5S47f/58PP/88/j4449Rr149fPHFF2LLneooaTWzadMm3L17F25ubhU+xSu5pilp7ZCcnAwXFxdx+tPldePGjbF06VLIZLJqx0k1U35+Pg4cOAClUim+IruwsBCPHz/G5cuXy1xXlyc1NRWffvop+vfvj+XLl2Pnzp3ltrZVd7yV3Fh8Wr9+/bBnzx44OzujZ8+eKjetKhIeHo433nijzFt0APXlSWFhISZPnoyoqCh07doVVlZWGD9+vE4G33Nzc8Pw4cPx3nvviZ9LxkUozcXFRavrsJJWE3/88UetHdCPzft1JCkpCZ9//jlSUlIAFB948fHxYlOy3Nxc2NnZoUGDBkhNTcWWLVtUlm/UqBHu3LlT7rrT09Nx+PBhPHnyBNbW1rCzs6vw9T25ubmwtraGo6Mj8vLy8Mknn4jTpFIpBgwYgGXLliE1NRUKhQKJiYkoLCyEk5MTpFJphTF06tQJN2/exN69e1FUVIT9+/fj2rVr6Ny5s7a7SkWzZs3w5ZdfljtwUOfOnXHz5k3ExcVBLpdDLpfj/Pnz4qAnubm5sLe3h42NDc6fP1+miZU2Vq5ciX///RdFRUXIycnBN998gxYtWsDR0REhISH47bffsH//fhQVFSErKwuXLl2ChYUFevXqhVWrViEnJwf37t1DTExMhU2qpFIpBg0ahKVLl4p3a1NTU/Hrr79WOW7SrTp16uCFF17AF198ofKUSiaTqXzn5uYGT09PfPLJJygoKMDly5fxww8/lBmArrSGDRvi3r17Wj3R9fX1xfbt28Xt+vn54csvv4RMJtP7K7wGDRqEb7/9FufOnYMgCHjy5AmOHTuGnJwcXL9+Hb///jsKCwthbW0NGxsbMZ7du3eLrYxKWjNVFivLAZYDNVleXh4kEon49G/nzp1lRpDPzMzEtm3bIJfLceDAASQlJak8vdqzZw+uXbuGvLw8rFmzBj179tQ6j4eGhuLo0aP49ddfoVAoUFBQgFOnTiElJUXtdcSBAwfEaxd7e3tIJBKVJ4zlsbOzw3//+18sWbKkzLT27dujXr162Lx5M/Lz86FQKPDvv/+KN81yc3NRt25d1K1bF0lJSfjmm2+0Smdp0dHROH/+vNiScdu2bWjQoAGeeeYZdO7cGenp6fjiiy9QWFiInJwcnDt3DkDxO9E3bNiAzMxMZGZmIjo6Wm35/dprr2H16tVihSszMxOHDh2qctxU8xw6dAgWFhbYt28f4uLiEBcXh/3798Pb27vMgH7lEQQBkZGRGDhwIJYuXQoXFxesXr263HnVHW8VXa+HhYXh0KFD2LNnT5kuNxXx9fVFTEyMODhgaerKk8LCQrHuYGlpiV9++aXKrw1/9OgR1q5di1u3bkGpVCIzMxM7d+4UuwYOHDgQu3btwu+//w6lUonU1FQkJSVpfR3m4uICf39/LF++HDk5OVAqlbh9+3aZ7oDmjJV+HalXrx7OnTuHQYMGoUOHDhg8eDBat24tjqg/ceJEXLx4Ed7e3hg9ejR69Oihsvzo0aOxYcMGeHt7Y+vWrSrTlEolYmJiEBgYKPY5njdvXrlxhIeHw93dHYGBgejbt2+Z/rQzZsxA69atMXDgQPj6+uKjjz6CUqmEra2t+MoMb29vnD17VmU5R0dHbNy4ETExMfDz88OWLVuwceNG8SKnOry9vct94lWvXj1s3boV+/fvR2BgIAICAvDRRx+JI+zOmzcPa9euhaenJ6Kjo6v12pT8/HxMnDgRPj4+6NatG+7fvy+O8u3u7o7PPvsMMTEx8PX1RXh4uNi/cM6cObC1tUW3bt0wbNgwBAcHY8CAARVuZ9q0aWjRogUGDx4MLy8vvPnmmxr3vSbD8PHxQUZGhsoddplMhoyMDPj4+IjfffLJJ7h37x4CAwMxceJETJo0Sbz7X56SJud+fn5l3nOtLpbc3FxxuzKZDPn5+Vq/MrAqXnzxRSxatAgLFy6Ej48PevToIY5WX1hYiI8//hh+fn4ICAhAZmYmpk6dCgD49ddf0bdvX3h6emLJkiVYtWqVRk8cWA6wHKip/vOf/yAiIgJDhw7FK6+8gn///VdsKVSiffv2uHXrFjp27IjVq1dj7dq14msqgeIL9sjISPj7+6OwsBCzZs3SOg43NzesX78emzZtwssvv4xOnTph69atUCqVaq8j/v77bwwaNAienp4YN24cZs2apVE3wBdffLHc5vMWFhbYsGEDLl++jK5du6Jjx46YPXs2cnJyABRfh8THx8PLywtz5sxBnz59tE5rCYlEgpkzZ6Jjx44IDAzEb7/9hk2bNqFu3bqoV68ePv/8cxw9ehT+/v7o2bOnOPr4+PHj8cILLyA0NBShoaFo164dxo8fX+F2Xn/9dQQFBSEiIgKenp4YPHhwuS1/yHzFxsaif//+cHd3h7Ozs/g3fPhw8YGYOtu2bUN6ejqmTJkCiUSCpUuXYteuXSr94EuoO94qul5v3Lgxnn/+eUgkEo2vESQSCV5++WU4ODiUmaauPKlXrx5mz56Nd999Fz4+PoiPjy8z8r+mrKyscO/ePYwaNQoymQwhISGwtrYWu0u3b98ey5YtE1s+jBgxQmx5rO112IoVKyCXy9GnTx/4+Phg8uTJOn8trymTCLpoi0FEREREZZT0/6/oifbIkSMRGhqKQYMGGTgyIjInH374IVxcXMSb8ESl1c5ODURERERERGbg7t27+PnnnxEbG2vsUMhEsXk/ERERERFRDbR69WqEhITgrbfeqtIbuqh2YPN+IiIiIiIiIjPFJ/1EREREREREZoqVfiIiIiIiIiIzxUo/ERERERERkZlipZ+IiIiIiIjITLHST0RERERERGSmWOknIiIiIiIiMlP/B6+8Od9gRS4AAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1224x288 with 8 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.set(style=\"darkgrid\")\n",
    "fig = plt.figure(figsize=(17, 4))\n",
    "gs = fig.add_gridspec(1, 4, wspace=0.25, hspace=0)\n",
    "\n",
    "gssub0 = gs[0].subgridspec(2, 1, height_ratios=(.15, .85)) \n",
    "ax0h = fig.add_subplot(gssub0[1])\n",
    "ax0b = fig.add_subplot(gssub0[0], sharex = ax0h)\n",
    "\n",
    "gssub1 = gs[1].subgridspec(2, 1, height_ratios=(.15, .85)) \n",
    "ax1h = fig.add_subplot(gssub1[1])\n",
    "ax1b = fig.add_subplot(gssub1[0], sharex = ax1h)\n",
    "\n",
    "gssub2 = gs[2].subgridspec(2, 1, height_ratios=(.15, .85)) \n",
    "ax2h = fig.add_subplot(gssub2[1])\n",
    "ax2b = fig.add_subplot(gssub2[0], sharex = ax2h)\n",
    "\n",
    "gssub3 = gs[3].subgridspec(2, 1, height_ratios=(.15, .85)) \n",
    "ax3h = fig.add_subplot(gssub3[1])\n",
    "ax3b = fig.add_subplot(gssub3[0], sharex = ax3h)\n",
    "\n",
    "sns.boxplot(df['Life Satisfaction'], ax=ax0b)\n",
    "sns.boxplot(df['Worthwhile'], ax=ax1b)\n",
    "sns.boxplot(df['Happiness'], ax=ax2b)\n",
    "sns.boxplot(df['Anxiety'], ax=ax3b)\n",
    "\n",
    "sns.histplot(data=df, x=\"Life Satisfaction\", ax=ax0h)\n",
    "sns.histplot(data=df, x='Worthwhile', ax=ax1h)\n",
    "sns.histplot(data=df, x='Happiness', ax=ax2h)\n",
    "sns.histplot(data=df, x='Anxiety', ax=ax3h)\n",
    "\n",
    "ax0b.set(xlabel='', title='Fig.1a) Dep. Var: Life Satisfaction')\n",
    "ax1b.set(xlabel='', title='Fig.1b) Dep. Var: Worthwhileness')\n",
    "ax2b.set(xlabel='', title='Fig.1c) Dep. Var: Happiness')\n",
    "ax3b.set(xlabel='', title='Fig.1d) Dep. Var: Anxiety')\n",
    "\n",
    "ax0h.set(ylabel='Ward Count')\n",
    "ax1h.set(ylabel='')\n",
    "ax2h.set(ylabel='')\n",
    "ax3h.set(ylabel='')\n",
    "\n",
    "ax0h.set(xlabel='Satisfaction Mean Score')\n",
    "ax1h.set(xlabel='Worthwhilness Mean Score')\n",
    "ax2h.set(xlabel='Happiness Mean Score')\n",
    "ax3h.set(xlabel='Anxiety Mean Score')\n",
    "\n",
    "plt.show()\n",
    "\n",
    "#https://www.python-graph-gallery.com/24-histogram-with-a-boxplot-on-top-seaborn\n",
    "#https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.gridspec.GridSpecFromSubplotSpec.html#matplotlib.gridspec.GridSpecFromSubplotSpec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1292dc13-66e3-4343-9013-7a15a0b3b168",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/conda/envs/sds2021/lib/python3.9/site-packages/seaborn/_decorators.py:36: FutureWarning: Pass the following variable as a keyword arg: x. From version 0.12, the only valid positional argument will be `data`, and passing other arguments without an explicit keyword will result in an error or misinterpretation.\n",
      "  warnings.warn(\n",
      "/opt/conda/envs/sds2021/lib/python3.9/site-packages/seaborn/_decorators.py:36: FutureWarning: Pass the following variable as a keyword arg: x. From version 0.12, the only valid positional argument will be `data`, and passing other arguments without an explicit keyword will result in an error or misinterpretation.\n",
      "  warnings.warn(\n",
      "/opt/conda/envs/sds2021/lib/python3.9/site-packages/seaborn/_decorators.py:36: FutureWarning: Pass the following variable as a keyword arg: x. From version 0.12, the only valid positional argument will be `data`, and passing other arguments without an explicit keyword will result in an error or misinterpretation.\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "sns.set(style=\"darkgrid\")\n",
    "fig = plt.figure(figsize=(17, 4))\n",
    "gs = fig.add_gridspec(1, 4, wspace=0.25, hspace=0)\n",
    "\n",
    "gssub0 = gs[0].subgridspec(2, 1, height_ratios=(.15, .85)) \n",
    "ax0h = fig.add_subplot(gssub0[1])\n",
    "ax0b = fig.add_subplot(gssub0[0], sharex = ax0h)\n",
    "\n",
    "gssub1 = gs[1].subgridspec(2, 1, height_ratios=(.15, .85)) \n",
    "ax1h = fig.add_subplot(gssub1[1])\n",
    "ax1b = fig.add_subplot(gssub1[0], sharex = ax1h)\n",
    "\n",
    "gssub2 = gs[2].subgridspec(2, 1, height_ratios=(.15, .85)) \n",
    "ax2h = fig.add_subplot(gssub2[1])\n",
    "ax2b = fig.add_subplot(gssub2[0], sharex = ax2h)\n",
    "\n",
    "sns.boxplot(df['Open Space'], ax=ax0b)\n",
    "sns.boxplot(df['Local Parks'], ax=ax1b)\n",
    "sns.boxplot(df['District Parks'], ax=ax2b)\n",
    "\n",
    "sns.histplot(data=df, x=\"Open Space\", ax=ax0h)\n",
    "sns.histplot(data=df, x='Local Parks', ax=ax1h)\n",
    "sns.histplot(data=df, x='District Parks', ax=ax2h)\n",
    "\n",
    "ax0b.set(xlabel='', title='Fig.3a) Ind. Var: Open Space')\n",
    "ax1b.set(xlabel='', title='Fig.3b) Ind. Var: Local Parks')\n",
    "ax2b.set(xlabel='', title='Fig.3c) Ind. Var: District Parks')\n",
    "\n",
    "ax0h.set(ylabel='Ward Count')\n",
    "ax1h.set(ylabel='')\n",
    "ax2h.set(ylabel='')\n",
    "\n",
    "ax0h.set(xlabel='% Households w/access to Open Space')\n",
    "ax1h.set(xlabel='% Households w/access to Local Parks')\n",
    "ax2h.set(xlabel='% Households w/access to District Parks')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa6221f3-8715-4c2a-9c20-b0cfcc05e635",
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set(style=\"darkgrid\")\n",
    "fig = plt.figure(figsize=(17, 4))\n",
    "gs = fig.add_gridspec(1, 4, wspace=0.25, hspace=0)\n",
    "\n",
    "gssub0 = gs[0].subgridspec(2, 1, height_ratios=(.15, .85)) \n",
    "ax0h = fig.add_subplot(gssub0[1])\n",
    "ax0b = fig.add_subplot(gssub0[0], sharex = ax0h)\n",
    "\n",
    "gssub1 = gs[1].subgridspec(2, 1, height_ratios=(.15, .85)) \n",
    "ax1h = fig.add_subplot(gssub1[1])\n",
    "ax1b = fig.add_subplot(gssub1[0], sharex = ax1h)\n",
    "\n",
    "gssub2 = gs[2].subgridspec(2, 1, height_ratios=(.15, .85)) \n",
    "ax2h = fig.add_subplot(gssub2[1])\n",
    "ax2b = fig.add_subplot(gssub2[0], sharex = ax2h)\n",
    "\n",
    "sns.boxplot(df['Metropolitan Parks'], ax=ax0b)\n",
    "sns.boxplot(df['Regional Parks'], ax=ax1b)\n",
    "sns.boxplot(df['% of open space that has access'], ax=ax2b)\n",
    "\n",
    "sns.histplot(data=df, x=\"Metropolitan Parks\", ax=ax0h)\n",
    "sns.histplot(data=df, x='Regional Parks', ax=ax1h)\n",
    "sns.histplot(data=df, x='% of open space that has access', ax=ax2h)\n",
    "\n",
    "ax0b.set(xlabel='', title='Fig.3d) Ind. Var: Metropolitan Parks')\n",
    "ax1b.set(xlabel='', title='Fig.3e) Ind. Var: Regional Parks')\n",
    "ax2b.set(xlabel='', title='Fig.3f) Ind. Var: % Open Space w/Access')\n",
    "\n",
    "ax0h.set(ylabel='Ward Count')\n",
    "ax1h.set(ylabel='')\n",
    "ax2h.set(ylabel='')\n",
    "\n",
    "ax0h.set(xlabel='% Hh w/access to Metrop. Parks')\n",
    "ax1h.set(xlabel='% Households w/access to Regional Parks')\n",
    "ax2h.set(xlabel='% Open Space with Access')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9448f853-dd4c-400f-b03c-9f2ddd298eec",
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set(style=\"darkgrid\")\n",
    "fig = plt.figure(figsize=(17, 4))\n",
    "gs = fig.add_gridspec(1, 4, wspace=0.25, hspace=0)\n",
    "\n",
    "gssub0 = gs[0].subgridspec(2, 1, height_ratios=(.15, .85)) \n",
    "ax0h = fig.add_subplot(gssub0[1])\n",
    "ax0b = fig.add_subplot(gssub0[0], sharex = ax0h)\n",
    "\n",
    "gssub1 = gs[1].subgridspec(2, 1, height_ratios=(.15, .85)) \n",
    "ax1h = fig.add_subplot(gssub1[1])\n",
    "ax1b = fig.add_subplot(gssub1[0], sharex = ax1h)\n",
    "\n",
    "gssub2 = gs[2].subgridspec(2, 1, height_ratios=(.15, .85)) \n",
    "ax2h = fig.add_subplot(gssub2[1])\n",
    "ax2b = fig.add_subplot(gssub2[0], sharex = ax2h)\n",
    "\n",
    "sns.boxplot(df['percent_green'], ax=ax0b)\n",
    "sns.boxplot(df['percent_blue'], ax=ax1b)\n",
    "sns.boxplot(df['percent_green+blue'], ax=ax2b)\n",
    "\n",
    "sns.histplot(data=df, x=\"percent_green\", ax=ax0h)\n",
    "sns.histplot(data=df, x='percent_blue', ax=ax1h)\n",
    "sns.histplot(data=df, x='percent_green+blue', ax=ax2h)\n",
    "\n",
    "ax0b.set(xlabel='', title='Fig.2a) Ind. Var: % Green Cover')\n",
    "ax1b.set(xlabel='', title='Fig.2b) Ind. Var: % Blue Cover')\n",
    "ax2b.set(xlabel='', title='Fig.2c) Ind. Var: % Green + Blue Cover')\n",
    "\n",
    "ax0h.set(ylabel='Ward Count')\n",
    "ax1h.set(ylabel='')\n",
    "ax2h.set(ylabel='')\n",
    "\n",
    "ax0h.set(xlabel='% Green Cover')\n",
    "ax1h.set(xlabel='% Blue Cover')\n",
    "ax2h.set(xlabel='% Green + Blue Cover')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88449a36-4ebe-4822-89ba-5fe9450edc5c",
   "metadata": {},
   "outputs": [],
   "source": [
    "df1 = df[['Open Space', 'Local Parks', 'District Parks', 'Metropolitan Parks', 'Regional Parks',\n",
    "                     '% of open space that has access']]\n",
    "\n",
    "df2 = df[['percent_green', 'percent_blue', 'percent_green+blue']]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7fc8978-445b-467b-94e7-1d337ea6fb58",
   "metadata": {},
   "source": [
    "## Multicollinearity\n",
    "(Multi)Collinearity occurs when two (or more) variables are highly correlated.\n",
    "- Visualising correlation with correlation matrices.\n",
    "- Detecting collinearity with Variance Inflation Factor (VIF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2b7c97c-eca1-43e9-b612-6adc47ad0877",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Correlation Matrix\n",
    "\n",
    "# Input data must be in columns of numbers without column labels, saved in a csv.\n",
    "# Each column is a data series, X0, X1, X2, etc.\n",
    "# resaving replaces existing file\n",
    "dfplot.to_csv('prematrix.csv', header=False, index=False)\n",
    "\n",
    "# read csv\n",
    "prematrix = np.genfromtxt('prematrix.csv',delimiter = ',')\n",
    "\n",
    "# generate correlation matrix\n",
    "correlation_matrix = np.corrcoef(prematrix.T)\n",
    "\n",
    "# correlation_matrix[i,j] is the correlation between the series Xi and Xj.\n",
    "print(correlation_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f97bf89e-92be-40e8-8411-bac1fe39c42e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Correlation Matrix Visualisation\n",
    "\n",
    "plt.rcParams[\"axes.grid\"] = False\n",
    "f = plt.figure(figsize=(15, 13))\n",
    "\n",
    "plt.matshow(dfplot.corr(), fignum=f.number)\n",
    "plt.xticks(range(dfplot.shape[1]), dfplot.columns, fontsize=14, rotation=80)\n",
    "plt.yticks(range(dfplot.shape[1]), dfplot.columns, fontsize=14)\n",
    "plt.title('Fig.4) Correlation Matrix', fontsize=16)\n",
    "\n",
    "cb = plt.colorbar()\n",
    "cb.ax.tick_params(labelsize=14)\n",
    "\n",
    "f.savefig(os.path.join('corr_matrix'), dpi=200)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c10399b2-16ab-4e62-badb-cbeb1861fd4f",
   "metadata": {},
   "source": [
    "I predict ... to be collinear"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09986bda-ac5f-49d7-8bcc-a2a9cbe50be6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculating VIF to drop collinear variables\n",
    "# This function is adjusted from: https://stackoverflow.com/a/51329496/4667568\n",
    "from statsmodels.stats.outliers_influence import variance_inflation_factor \n",
    "from statsmodels.tools.tools import add_constant\n",
    "\n",
    "def drop_column_using_vif_(df, thresh=5):\n",
    "    '''\n",
    "    Calculates VIF each feature in a pandas dataframe, and repeatedly drop the columns with the highest VIF\n",
    "    A constant must be added to variance_inflation_factor or the results will be incorrect\n",
    "\n",
    "    :param df: the pandas dataframe containing only the predictor features, not the response variable\n",
    "    :param thresh: (default 5) the threshould VIF value. If the VIF of a variable is greater than thresh, it should be removed from the dataframe\n",
    "    :return: dataframe with multicollinear features removed\n",
    "    '''\n",
    "    while True:\n",
    "        # adding a constatnt item to the data\n",
    "        df_with_const = add_constant(df)\n",
    "\n",
    "        vif_df = pd.Series([variance_inflation_factor(df_with_const.values, i) \n",
    "               for i in range(df_with_const.shape[1])], name= \"VIF\",\n",
    "              index=df_with_const.columns).to_frame()\n",
    "\n",
    "        # drop the const\n",
    "        vif_df = vif_df.drop('const')\n",
    "        \n",
    "        # if the largest VIF is above the thresh, remove a variable with the largest VIF\n",
    "        # If there are multiple variabels with VIF>thresh, only one of them is removed. This is because we want to keep as many variables as possible\n",
    "        if vif_df.VIF.max() > thresh:\n",
    "            # If there are multiple variables with the maximum VIF, choose the first one\n",
    "            index_to_drop = vif_df.index[vif_df.VIF == vif_df.VIF.max()].tolist()[0]\n",
    "            print('Dropping: {}'.format(index_to_drop))\n",
    "            df = df.drop(columns = index_to_drop)\n",
    "        else:\n",
    "            # No VIF is above threshold. Exit the loop\n",
    "            break\n",
    "\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "399fa284-c604-4bfc-b95b-1ea1ac975fd1",
   "metadata": {},
   "outputs": [],
   "source": [
    "noncoll = drop_column_using_vif_(\n",
    "    dfplot.drop('Life Satisfaction',axis=1)\n",
    "            .drop('Worthwhile',axis=1)\n",
    "            .drop('Happiness',axis=1)\n",
    "            .drop('Anxiety',axis=1)\n",
    "#            .drop('Open Space',axis=1)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e774114d-bd61-4490-819d-c3bec39499c3",
   "metadata": {},
   "source": [
    "## Multiple Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86de7e25-9e5b-4a6a-9f27-1ebb0ba64dba",
   "metadata": {},
   "outputs": [],
   "source": [
    "import statsmodels.api as sm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5bf7f1eb-88c5-463c-894c-564f6f18c746",
   "metadata": {},
   "outputs": [],
   "source": [
    "satis = sm.OLS(endog = dfplot[['Life Satisfaction']], \n",
    "              exog = sm.add_constant(noncoll)).fit()\n",
    "worth = sm.OLS(endog = dfplot[['Worthwhile']], \n",
    "              exog = sm.add_constant(noncoll)).fit()\n",
    "happy = sm.OLS(endog = dfplot[['Happiness']], \n",
    "              exog = sm.add_constant(noncoll)).fit()\n",
    "anxi = sm.OLS(endog = dfplot[['Anxiety']], \n",
    "              exog = sm.add_constant(noncoll)).fit()\n",
    "\n",
    "#satis.summary()\n",
    "worth.summary()\n",
    "#happy.summary()\n",
    "#anxi.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3cb91364-c4be-4b4c-b9fb-b0de4138200b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualise Residuals \n",
    "\n",
    "fig1, ax = plt.subplots()#figsize=(8,6))\n",
    "ax.scatter(x=satis.fittedvalues , y=satis.resid)\n",
    "ax.axhline(y=0, color = 'black', linestyle='-')\n",
    "plt.title('Fig.5a) Life Satisfaction: Residuals vs. Fits')\n",
    "plt.xlabel('Fitted Values')\n",
    "plt.ylabel('Residuals')\n",
    "plt.show()\n",
    "\n",
    "fig2 = plt.hist(satis.resid)\n",
    "plt.title('Fig.5b) Life Satisfaction: Residuals')\n",
    "plt.xlabel('Residuals')\n",
    "plt.ylabel('Counts')\n",
    "plt.show()\n",
    "\n",
    "\n",
    "\n",
    "fig3, ax = plt.subplots()#figsize=(8,6))\n",
    "ax.scatter(x=worth.fittedvalues , y=worth.resid)\n",
    "ax.axhline(y=0, color = 'black', linestyle='-')\n",
    "plt.title('Fig.6a) Worthwhile: Residuals vs. Fits')\n",
    "plt.xlabel('Fitted Values')\n",
    "plt.ylabel('Residuals')\n",
    "plt.show()\n",
    "\n",
    "fig4 = plt.hist(worth.resid)\n",
    "plt.title('Fig.6b) Worthwhile: Residuals')\n",
    "plt.xlabel('Residuals')\n",
    "plt.ylabel('Counts')\n",
    "plt.show()\n",
    "\n",
    "\n",
    "\n",
    "fig5, ax = plt.subplots()#figsize=(8,6))\n",
    "ax.scatter(x=happy.fittedvalues , y=happy.resid)\n",
    "ax.axhline(y=0, color = 'black', linestyle='-')\n",
    "plt.title('Fig.7a) Happiness: Residuals vs. Fits')\n",
    "plt.xlabel('Fitted Values')\n",
    "plt.ylabel('Residuals')\n",
    "plt.show()\n",
    "\n",
    "fig6 = plt.hist(happy.resid)\n",
    "plt.title('Fig.7b) Happiness: Residuals')\n",
    "plt.xlabel('Residuals')\n",
    "plt.ylabel('Counts')\n",
    "plt.show()\n",
    "\n",
    "\n",
    "fig7, ax = plt.subplots()#figsize=(8,6))\n",
    "ax.scatter(x=anxi.fittedvalues , y=anxi.resid)\n",
    "ax.axhline(y=0, color = 'black', linestyle='-')\n",
    "plt.title('Fig.8a) Anxiety: Residuals vs. Fits')\n",
    "plt.xlabel('Fitted Values')\n",
    "plt.ylabel('Residuals')\n",
    "plt.show()\n",
    "\n",
    "fig8 = plt.hist(anxi.resid)\n",
    "plt.title('Fig.8b) Anxiety: Residuals')\n",
    "plt.xlabel('Residuals')\n",
    "plt.ylabel('Counts')\n",
    "plt.show()"
   ]
  }
 ],
 "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.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}