notebooks/GWL-plot_Python.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3fd786b1-fd60-4f0f-aab4-383980aebb59",
   "metadata": {},
   "source": [
    "# Atlas notebooks\n",
    "\n",
    "***\n",
    "\n",
    "> This notebook reproduces and extends parts of the figures and products of the AR6-WGI Atlas. It is part of a notebook collection available at https://github.com/IPCC-WG1/Atlas for reproducibility and reusability purposes. This work is licensed under a [Creative Commons Attribution 4.0 International License](http://creativecommons.org/licenses/by/4.0).\n",
    ">\n",
    "> ![Creative Commons License >](https://i.creativecommons.org/l/by/4.0/88x31.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7b4f5e7-524e-4090-be33-f70c53c4a8c2",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "## Creating GWL plots\n",
    "\n",
    "11/02/2021\n",
    "\n",
    "*J. Fernández* - [Santander Meteorology Group](http://www.meteo.unican.es). Instituto de Física de Cantabria, CSIC-UC, Santander, Spain.\n",
    "\n",
    "> This notebook is a companion of the article _Implementation of FAIR principles in the IPCC: The WGI AR6 Atlas repository_ by Itubide et al., submitted to Scientific Data. It shows how to reuse the repository resources to create GWL plots as shown in the panels of Figure 3.\n",
    "\n",
    "Global Warming Level (GWL) plots summarize the regional climate sensitivity of changes in different variables as a function of global mean surface air temperature changes. Decadal mean changes for the different members of a multimodel ensemble provide the sample to study the relationship between regional climate sensitivity and GWL. These plots are generated by the IPCC WGI AR6 Interactive Atlas (http://interactive-atlas.ipcc.ch) and are reproduced in this notebook using the data available in the Atlas repository. As in the interactive Atlas, the notebook allows to set different variables (tas or pr), masks (land, sea, landsea), scenarios, reference regions, etc."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b56d3b0e-55a3-4cc3-af04-af88111e9f51",
   "metadata": {},
   "source": [
    "### Import packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "4a60058d-67f4-4612-8610-430367976c87",
   "metadata": {},
   "outputs": [],
   "source": [
    "import glob\n",
    "import matplotlib.pyplot as plt\n",
    "import os.path\n",
    "import pandas as pd\n",
    "from matplotlib import cm\n",
    "from sklearn import linear_model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49b14a70-3a89-49e3-95d0-ebc0c86c8d78",
   "metadata": {},
   "source": [
    "### Set parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e46d8cbf-f5c0-4beb-8d6b-1a3692cb3822",
   "metadata": {},
   "outputs": [],
   "source": [
    "var = 'tas'\n",
    "mask = 'landsea'\n",
    "project = 'CMIP6' # CMIP5, CMIP6 or CORDEX\n",
    "project_gsat = 'CMIP6' # CMIP5 or CMIP6\n",
    "scenario = 'ssp585' # e.g. ssp585 or rcp45\n",
    "region = 'MED'\n",
    "season = 'Annual'\n",
    "period_hist = slice('1850','1900')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc7dd615-714f-49e8-ae5f-786b7e7ba307",
   "metadata": {},
   "source": [
    "The following, in principle, require no modification. It is just information for labels and season selection. The abscissas in this plot always represent the GWL so the variable is always `tas` and no mask is applied (`landsea`), regardless of the selection above for the variable in the Y axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "7bea44c0-4efd-49d8-a654-93e637588d59",
   "metadata": {},
   "outputs": [],
   "source": [
    "basedir = f'../datasets-aggregated-regionally/data/{project}/{project}_{var}_{mask}'\n",
    "basedir_gsat = f'../datasets-aggregated-regionally/data/{project_gsat}/{project_gsat}_tas_landsea'\n",
    "longname = dict(pr = 'Precipitation', tas = 'Near surface temperature')\n",
    "units = dict(pr='%', tas='K')\n",
    "months = dict(DJF=[1,2,12], MAM=[3,4,5], JJA=[6,7,8], SON=[9,10,11], Annual=range(1,13))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f4926ac-14d5-4032-8e6a-fa9ea0dc62dd",
   "metadata": {},
   "source": [
    "### Compute decadal anomalies\n",
    "\n",
    "We define first some commodity functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "9a40c044-08e9-4732-879a-fb3b06dc6383",
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_files(csvlist):\n",
    "  '''Concatenates a list of CSV files as a dataframe'''\n",
    "  csvdata = [pd.read_csv(csv, comment = '#') for csv in csvlist]\n",
    "  csvdata = pd.concat(csvdata)\n",
    "  csvdata['date'] = pd.to_datetime(csvdata['date'], format='%Y-%m')\n",
    "  return(csvdata.set_index('date'))\n",
    "\n",
    "def get_average(df, period, region = 'world', season = 'Annual'):\n",
    "  '''Computes the seasonal average for a given region and period'''\n",
    "  rval = df.loc[period,region]\n",
    "  return(rval[rval.index.month.isin(months[season])].mean())\n",
    "\n",
    "def get_run(filepath):\n",
    "  '''Extracts member identifier from filename'''\n",
    "  return( filepath.split('_')[-1].split('.')[0] )\n",
    "\n",
    "def get_model(filepath):\n",
    "  '''Extracts model name from filename'''\n",
    "  return( filepath.split('_')[-3].split('.')[0] )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36959718-9413-4f70-8531-607e0358a7f8",
   "metadata": {},
   "source": [
    "We can now compute decadal anomalies for the selected variable and for global mean surface air temperature as reference:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ddb691f0-5839-4154-92b9-b868386fe668",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def decadal_anomalies(basedir, project, region = 'world', season = 'Annual', relative = False):\n",
    "  data = pd.DataFrame(columns = ['model', 'run', 'decade', region])\n",
    "  files_scen = glob.glob(f'{basedir}/{project}_*_{scenario}*.csv')\n",
    "  for scenfile in files_scen:\n",
    "    histfile = scenfile.replace(scenario, 'historical')\n",
    "    if not os.path.exists(histfile):\n",
    "      print(f'Missing historical file for {scenfile}')\n",
    "      continue\n",
    "    member_data = read_files([histfile, scenfile])\n",
    "    reference_value = get_average(member_data, period_hist, region, season)\n",
    "    # get decadal anomalies w.r.t. reference\n",
    "    for decade_start in range(2010,2099,10):\n",
    "      decade = slice(str(decade_start), str(decade_start+9))\n",
    "      dfrow = dict(\n",
    "        model = get_model(scenfile),\n",
    "        run = get_run(scenfile),\n",
    "        decade = f'{decade.start}-{decade.stop}'\n",
    "      )\n",
    "      dfrow[region] = get_average(member_data, decade, region, season) - reference_value\n",
    "      if relative:\n",
    "        dfrow[region] = 100. * dfrow[region] / reference_value\n",
    "      data = data.append(dfrow, ignore_index=True)\n",
    "  return(data.set_index(['model', 'run', 'decade']))\n",
    "\n",
    "ydata = decadal_anomalies(basedir, project, region, season, relative = var=='pr')\n",
    "xdata = decadal_anomalies(basedir_gsat, project_gsat)\n",
    "data = pd.concat([xdata, ydata], axis=1, join='inner').dropna()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ac4a4c1-0129-4f32-bf46-37ab0c782593",
   "metadata": {},
   "source": [
    "We have now all the information we need to construct the GWL Plot in a single dataframe object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "5cbbf211-a1bb-4f73-8af4-c1e6d4996edf",
   "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></th>\n",
       "      <th></th>\n",
       "      <th>world</th>\n",
       "      <th>MED</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>model</th>\n",
       "      <th>run</th>\n",
       "      <th>decade</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"5\" valign=\"top\">CAMS-CSM1-0</th>\n",
       "      <th rowspan=\"5\" valign=\"top\">r2i1p1f1</th>\n",
       "      <th>2010-2019</th>\n",
       "      <td>0.770013</td>\n",
       "      <td>0.856081</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2020-2029</th>\n",
       "      <td>0.898854</td>\n",
       "      <td>1.085415</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2030-2039</th>\n",
       "      <td>1.024613</td>\n",
       "      <td>1.087681</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2040-2049</th>\n",
       "      <td>1.429946</td>\n",
       "      <td>1.660690</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2050-2059</th>\n",
       "      <td>1.745204</td>\n",
       "      <td>2.003315</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <th>...</th>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"5\" valign=\"top\">INM-CM4-8</th>\n",
       "      <th rowspan=\"5\" valign=\"top\">r1i1p1f1</th>\n",
       "      <th>2050-2059</th>\n",
       "      <td>2.341992</td>\n",
       "      <td>2.526830</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2060-2069</th>\n",
       "      <td>2.824592</td>\n",
       "      <td>2.974939</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2070-2079</th>\n",
       "      <td>3.203709</td>\n",
       "      <td>3.524455</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2080-2089</th>\n",
       "      <td>3.508034</td>\n",
       "      <td>4.190764</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2090-2099</th>\n",
       "      <td>3.930476</td>\n",
       "      <td>4.518714</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>306 rows × 2 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                                   world       MED\n",
       "model       run      decade                       \n",
       "CAMS-CSM1-0 r2i1p1f1 2010-2019  0.770013  0.856081\n",
       "                     2020-2029  0.898854  1.085415\n",
       "                     2030-2039  1.024613  1.087681\n",
       "                     2040-2049  1.429946  1.660690\n",
       "                     2050-2059  1.745204  2.003315\n",
       "...                                  ...       ...\n",
       "INM-CM4-8   r1i1p1f1 2050-2059  2.341992  2.526830\n",
       "                     2060-2069  2.824592  2.974939\n",
       "                     2070-2079  3.203709  3.524455\n",
       "                     2080-2089  3.508034  4.190764\n",
       "                     2090-2099  3.930476  4.518714\n",
       "\n",
       "[306 rows x 2 columns]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "915e909a-404c-48f4-92ec-a7fe4b23c51e",
   "metadata": {},
   "source": [
    "### GWL Plot\n",
    "\n",
    "The GWL Plot is now just a scatterplot of the region vs. the world columns in this dataframe. We might colour the individual points according to their decade, and add decadal averages on top, as in the Interactive Atlas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "ae1ce6a6-c557-4cb9-8ad8-fd2e7e6c6b7b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWgAAAEGCAYAAABIGw//AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAB0BElEQVR4nO2dd3iURdeH79nNpvcCgUAKJUBCFRGkSJEqooJiAcUOqCigYuGzoK/6qqiIvjYUVBQURVQQQWro0nvoJYGQSkhvm935/tjskk1jQ7JpzH1de7HPPM80dnMyOTPnd4SUEoVCoVDUPTS1PQCFQqFQlI0y0AqFQlFHUQZaoVAo6ijKQCsUCkUdRRlohUKhqKM41PYAiuPt7S1btWpV28OwO9nZ2bi5udX2MOzOtTJPqP9z3b17d4qUMqC2x6Gwpk4Z6MaNG7Nr167aHobdiYqKol+/frU9DLtzrcwT6v9chRAxtT0GRWmUi0OhUCjqKMpAKxQKRR1FGWiFQqGooygDrVAoFHUUZaAVCoWijqIMtEKhUNRRlIFWKBSKOooy0ApFA2TNgo24OXgghGBMyETWLtxU20NSXAV1KlBFoVBUnX/mr2fogwMs1+dizzNr/JcA3DymT20NS3EV1AsD/fHHH5OWllbbw6g2EhISiIqKqu1h2BVvb286d+5c28O45sjKyrIyzkGE4SxcyM8pYN70hcpA1zPqhYFOS0tjxowZtT2MauPs2bOEhobW9jDsSkP6vOoLSUlJNG7c2HLdgkhaiHaW6+RzF2tjWIoqYFcftBBiqhDisBDikBDiJyGEsz37UyiuVU6dOmVlnNvR1co4AwQ096vpYSmqiN0MtBAiCHgGuF5K2R7QAvfaqz+F4lpl165dFFeBfOv592jh1sbqGSdXRx55Z0xND01RRex9isMBcBFCOACuwAU796dQXFP8888/dOvWzXK9bds2/m/mC0ydM5FGwf4IIWgU7M/UOROV/7keYjcftJQyTgjxARAL5AKrpJSrSj4nhBgPjAcICAgoc/MsISGBs2fP2muoNU5BQUGDmk9ZJCQkkJWV1eA3Q83UxlxXrVrFf//7X8v1999/T15eHlFRUWibwuPf3VPsacM181k0JOxmoIUQPsDtQBiQBvwqhLhfSvlj8eeklHOAOQBt2rSRZWnqRkVFNahNtWthkzAwMBB3d/d6rZFcGWpaD/r999+3Ms4XLlygSZMmNda/omawp4tjIHBGSpkspdQDS4CeduxPobgmePrpp3nxxRct1+np6co4N1DsecwuFughhHDF5OK4GWj46VIUCjsyYsQI/vrrL8t1Xl4eTk5OtTgihT2x2wpaSrkdWAzsAQ4W9TXHXv0pFLXF2oWb+PTe7xmsvZuxoU9Ue1j12oWbGBv6BO7Cy8o4GwwGZZwbOHY9xSGlfF1K2VZK2V5K+YCUMt9efR08eJBbbrkFPz8//Pz8GDlyJElJSdXS9smTJ3n88ceJjIxEo9EghMDBwfY/PkrWDwsLK7d+Vfp66KGHEEIghLDaENq1axeenp4IIfD29mbfvn1W9ZYvX44Qgl9//bVK7VyLrF24iVnjvyQjKQspJUmxKcwa/2W1Gem1Czfx0eNf8FPMV2STAYAOJ9Ys2IhGo6R0GjoN4hOOi4ujf//+CCFYsGABX3zxBRs3bmTWrFnV0v6hQ4f45ptviI6ORkpp1/pV7askR44cYdiwYWRmZuLi4sKyZctKhWD/9ddf6HQ6hgwZUqV2rkXmTV9Ifk6BVZk5rLo6+OblH/k75yckpu+CJ770FSOqrX1F3aZBGOhPPvkET09PlixZwtChQ7n77rvp1q0bmZmZ1dJ+UFAQ06dPZ9myZdxwww12rV/VvooTGxvL4MGDSUlJwcHBgV9++YU+fUqfhf3777/p06cPnp6eVWrnWqS88OnqCKvOzc3l59jLXsHGNOcGMaDa2lfUfWw20EIINyGE1p6DuVqWL1/OyJEjLf64S5cusXnzZqsD/FWhW7duvP3229x66624uLjYtX5V+zKTlJTEoEGDOH/+PEIIvvvuO2699dZSz+3fv5/Y2Ngy71WmnWuV8sKnqxpWnZqaiqurq+U6hDZ0EN2rrX0z2dnZbNu2rVraUlQ/5RpoIYRGCDFGCLFcCJEEHAXii7Q1ZgohWtfcMMsnOzubI0eOWFbMmzZtYujQoTRr1ox77rmn1PNRUVGlfKsNkQkTJnD8+HEAZs+ezdixY8t8zrzpVJ7RtbWda5VH3hmDk6ujVVlVw6rPnTuHn99lAxzh2JXWokO1tW8mNjaW3r17M3To0Cq3pbAPFe0+rQfWAC8Dh6SURgAhhC/QH3hXCPF7ycCTmubAgQMYjUa6dOlCSEgIly5dwsnJiaioKJydS2szCSHQarUIIWphtDWHWZ514sSJPP300+U+99dffxEeHk7r1mX/vrW1nWsVc/j058/OIzM5m4DmfjzyzpirDqs+dOgQHTpcNsa//PILvvpA5k1fSPK5i1Vu38y2bdsYOXKkyY3y88/ccsstVWpPYR8qMtADiwJMrJBSpgohlkkpfxNC6Ow4NpvYt28f7u7utGjRgsWLF3PixAk+/vhjhg8fzuHDhwkMDLR6vm/fvhQWFtbSaGuexYsXM2XKFNq0aVPqXkpKCjt27GDKlClVauda5+YxfdA2NVQ5knDjxo307dvXch0VFWW5rm4djYMHD+Lm5sb69etp167dlSsoaoWKfNAvlVUohPACVgGUZcBrmr1799KpUyd0Oh0DBgxgwoQJLF26lNTUVDZs2FDbw6s1evfuDZiM8ODBgzl//nypZ5YvX47RaKzQp2xLO4qq89tvv1kZ5/3791tdVwdGo5GDBw8CMH78eA4cOKCMcx2nIgPdRwjxdvECIUQgsBFYZ9dRVYJ9+/bRpUsXqzKza6O4Pu61xptvvslDDz0EmHyNQ4YMITU11eqZv/76Cy8vL4sRvtp2FFXjs88+46677rJcx8TE0LFjx2rtIzMzk5EjR9KjRw9iY2MBcHNzq9Y+FNVPRQb6NqCTEOIjgKJNwc3A51LKN2ticFfCYDBw6NAhGjVqZFW+YMECfH19KzQ8lUGv15OSkkJKSgp6/eU/Gsxl+fmX429mzJhhCfIwK9ZVpn5lnq0IIQTffPMNt912GwDR0dHccsstZGdnW/pZtWoVQ4YMQacr31N1pXYUVePll19m0qRJluuUlBSCg4OrtY+YmBh69erF8uXLee+992jevHm1tq+wI1LKcl+ADvgN+BmIAUZW9HxVX+Hh4bIsXn/99TLLDx06JAHp7u4uP/vsM7lu3Tr55ptvSicnJzlv3rwy60RFRUmtViujoqLKvF8W69evl0C5r2+//dZqrObyM2fOVLp+ZZ4tyYMPPmh5bv369VJKKXNzc+VNN91kKR8yZIgsKCiQa9askYCcP39+ldopj9dff91S91rgSnNds2CjHBMyUQ7SjJZjQibKNQs2yrFjx1p9tjk5OdU+rs2bN8uAgADp5eUlV61aVe5zwC5px59t9bq6V7mbhEKIZ4ve7gBeADYBYeZyKeVHV/k7odrYt28fOp2OkSNH8tJLLyGEoF27dixYsIA777yzzDpSSgwGg/kXUIPH2dmZZcuW0bdvX/bt28c///zDuHHjCAwMRKPRMGzYsCq1s2DBAhVyfAXM4eDmiMOk2BTuemAUacYUyzN6vb5S8gG28ssvv+Dt7c2yZcvUBm89RJRnqIQQr1dUUUr5RnUPpk2bNvLYsWOlymfMmFFmEtJp06axatUq9u/fX91DsSt1QQ86PDycgIAAtmzZYpf2Z8yYQb9+/ZQeNDA29AmSYi8b4w1yGXpMriohNKxeuoubb+1SZt2rwWAwEB8fT7NmzdDr9WRnZ+Pt7V1hHSHEbinl9dU2CEW1UO6vbHsY4Opm7969pTYIFbZhDj5R2B9zWLaUkrX8Zil3xYM+kU/zyRt/IoRgwPDOVe4rIyODsWPHcuDAAQ4dOoSHh8cVjbOi7lJRJOErRVlRyrs/QAhRqzG/+/fvV4I9ijpPQHM/jNJoZZz9CKSn43AA8vP0fDd7dZX7OX36ND179mTFihW8+OKLeHh4VLlNRe1SkdPrIPCXECIPk6ZzMuAMtAY6Y4oyfMfeA6yI5OTk2uxeobCJe165ndsfv+zrb0ZL2mquQzTyt5QlJ6RVqY8NGzZw5513YjQa+eeff7j55pur1J6iblDuClpK+aeUshcwETgMaIEM4EfgBinlVCmlspCKBoNZGL86hfcTEhKsjHNL2tPWsRuiaWM03pfVAwMCva+6Dykl77zzDv7+/mzfvl0Z5wbEFbeNpZQngBOVbVgI0QZYVKyoBfCalPLjyralUNibsk5azBr/JXD1YdbHjx+3Ojnx3Xff0dy/E7Nn/EF+3uUz7k7OOh6aPKjS7RcWFpKdnY2Xlxc//fQTGo1G+ZsbGPZMeXVMStlZStkZ6ArkAL/bqz+FoipUt/D+jh07rIzzihUrePDBBxkwvDOTZ9xBoybeCAHuFOIQe453b3unUqv29PR0RowYwYgRIzAYDPj6+irj3ACxZ9LY4twMnJJSxtRQfwpFpahO4f2///6b4cOHW6537NhhpU0+YHhnBgzvfNWr9pMnTzJixAhOnjzJZ599hlZrH5n23bt3N3JwcPgGaE8DSe5RxzAChwoLCx/r2rVrmfn5aspA3wv8VNYNIcR4YDxAQEBAmTrNCQkJlrDpq+XChQtMmzaNxMRENBoN/fv3twS31DQFBQVVnk9dJyEhgaysrHqju+0R4EZGUlaZ5bbMwTzXFStW8P7771vKf/zxR7Kzs8ts4/Nnvy9z1f75s/PQNjWU2c/evXstMQEzZ84kPDzcbv/HDg4O3wQGBrYLCAi4pNForo3IrhrEaDSK5OTkiISEhG8wSWuU4ooGWggRDnwBNJZSthdCdARuk1K+ZcsghBCORZ2/XNZ9KeUcirJ9t2nTRpZ12D8qKqrKgR1OTk7Mnj2b66+/noKCAgYNGsSePXvKjTi0J3UhUMXeBAYG4u7uXm8CVQwfaa1Ws2ASxn/yo0fo1+/KPuioqCi2bt1qZZwTEhIqFOx6J/nzMsszk7PL/H8rLCxkwoQJBAUFsWzZMlq2bHnFcVWR9so42w+NRiMDAgLSExIS2pf7jA3tfI3JuOoBpJQHMK2IbWUYsEdKmViJOpVi6tSpFoEiIQSNGzdm3LhxXLx4+c/TJk2acP31pkApR0dHOnbsyLlz56p1HLGxsdx///0EBATg4uJCp06dmDt3rs31zbkB3dzc8PDwYMCAAWzevLnc53/44Qd69uyJh4cHbm5utG3blhdeeKE6pnLNcfOYPkydM5FGwf4IIWgU7M/UORNt3iCcNWsW//d//2e5zsjIuKKaoq3psgoLCy2h4H/99Rfbtm2rCeMMoFHG2b4U/f+Wa4dtcXG4Sil3lHAFVEbx/j7KcW9UFwcPHqRHjx7MmjULo9HI5s2befnllykoKODnn38u9fzFixf5448/WLVqVbWNISEhgZ49exIXF2cpO3DgAI899hgJCQlWP7xlMX/+fB566CErjZD169czYMAAVqxYUero1KRJk/jss8+syo4dO8bChQutVnEK27l5TJ+rOrExbNgwVq5cabnOz8/H0dGxghomHnlnTJmr9uLprC5dusQ999xDaGgoc+bMKTfzjaJhYssKOkUI0RKT4hZCiLuAeFsaF0K4AoOAJVc9Qhs4ePAgXbt2pUePHvTs2ZMXXniB4cOHs3p16eis/Px87rrrLqZMmVKtYuWvv/66xTj/8MMPxMfH0727KcnnG2+8YdHgLYvCwkKee+45pJS0aNGCM2fOcOrUKUJCQtDr9TzxxBNWhvuvv/6yGOfRo0dz/PhxcnJyOHToENOnT6+2OSmuTHh4uJVxNhgMNhlnuPKq/fjx4/To0YOoqCh69Ohhl/HXdU6ePKnr3r17eIsWLSJbtWoV+Z///KcRQGJiorZnz56tQ0JC2vfs2bN1cnKyFiAhIUHbvXv3cFdX1y7jxo2z0m3dtGmTa3h4eERwcHD7hx56qLnRaCzVX2ZmpqZfv36twsLCIlu1ahX55JNPBpnv5ebmiuHDh7cIDg5u37Fjx7bHjh2zfNB9+vRp7eHh0bl///6tire3dOlSj4iIiHatW7eOHDVqVGhxCWFbsMVAPwV8BbQVQsQBUzAFr1wRKWWOlNJPSpleqVFVgqSkJJKSkkoZ20aNGpXa3TYYDIwdO5YuXbrw3HPPVdsYjEYjixaZjny3a9eO+++/n8DAQJ591iQIqNfrWbx4cbn1Dx06REqKSUznzjvvJDQ0lBYtWjBq1CgATpw4wc6dOy3Pf/LJJwCEhoby448/0rp1a1xcXIiMjOTJJ5+stnkpykdKiRCCEydMIQKODm48MOh9olYcqFQ7N4/pw4KzX7DK8AsLzn5hMc5r1qyhe/fupKamsnbtWh555JFqn0N1Y8xe6GtM6tXBmNCmqzGpVwdj9kLfqrap0+n48MMPz58+ffrwzp07j8ydO7fR7t27nV9//fUm/fr1y4yJiTnUr1+/zNdeey0QwNXVVb755psXZsyYUSr1z5NPPhny+eefx5w9e/bQ6dOnnRcvXuxZukd47rnnEs+cOXP40KFD0du3b3f/5ZdfPAFmz57t7+XlVRgbG3to0qRJic8++2wzc53nn38+4auvvjpTvB2DwcD48ePDfv7559MnTpw4HBwcXPC///3Pv2R/FVGhgRZCaIEnpJQDgQCgrZSyd106LnfggOkHom3btpYyo9HItm3bGDFihNWzEyZMwMPDgw8//LDCNiub+fvUqVOkp5t+B4WHh1vKi7/fu3dvufVzc3Mt78s7VWKubzAYLAp0zZo1Y+TIkXh7e+Pr68vYsWOJj7fpjxtFJSgZYbjqh/VWEqvers3o33YySfFpzJ7xB+uW76tSf+np6YwePZpmzZqxY8cO+vSp3nyE9sCYvdCXzP+GYEx2BAnGZEcy/xtSVSMdEhKi7927dw6Aj4+PsWXLlrmxsbGOK1eu9J4wYcJFgAkTJlxcsWKFD4Cnp6dxyJAhWc7OzlbL45iYGF1WVpZm4MCB2RqNhrFjx178448/SmkNeXh4GEeMGJEJ4OzsLDt27Jhz7tw5R4C//vrL+5FHHrkI8PDDD1/aunWrh3kVfvvtt2d6enpa9ZmYmOjg6Oho7NixYz7A0KFDM/744w/vysy/QgMtpTRgCjJBSpktpcysTOM1gdlAt2rVisLCQmJjY5kwYQKenp5WvtgtW7Ywd+5cdu3aRZcuXejcubNlJVqSymb+Nq9+ATw9Pct8n5RU5jFHwPTLxZzVZPHixcTExHDmzBmWLLnsGTJveF68eJGcnBwANm/ezN9//016ejqXLl1i4cKF9OvXz3JfUXXMZ5WTYlOQUhIfk8CQcQMs95t6t6d7i3GW66oIHxkMpqN1Xl5eLF++nK1btxIWFla1CdQU2Z8FQX4Je5KvMZVXD8eOHXOMjo527du3b9bFixcdQkJC9GAy4qmpqRXup8XExOiaNGli8S+EhIQUxMfHV5j0OiUlRbt69WrvYcOGZQAkJiY6hoWFFYBpZe/u7m5ITEwst9/AwMDCwsJCsXHjRleARYsW+cTHx9vm/yrClk3CvUKIpcCvgCXPkZTSrn5lWzEnwSx+bK1Zs2bs3LkTP7/Lu+G9evWyWaS/spm/y2u3eHlFxt7Hx4dJkyYxa9YsTp8+XeYRPLMBL+nD+v777xk5ciTTpk3jq6++4vjx4yxYsIDHH3/c5vErLrN24SbmTV9I8rmLBDT3Iy8737KJVyDz2cgyy7MtAnrSunG/Um1cjfBRamoqo0ePZtSoUTz11FP07NnzaqdQOxhTyjY85ZVXkvT0dM2oUaNavvvuu+d8fX1LO4+vQFk/oxX9TOr1ekaNGtVi/PjxiREREQUVtFGuUdFoNMyfP//01KlTmxcUFGj69++fXtmgIlt80L7ARWAAMKLoVasyo8Uxn+DYuXMnW7du5Y033uD8+fO89ZZNx7SrhYCAAMv7jIwMy/vMzMwynymLDz74gHfeeYfQ0FCcnJxo374948ePt9w355Hz8fGxfLF8fHwYN24cHh4eVr7n+pbAoK5QcrWcFJtCxkXTZ5grs62Mc1vRhV6d7yizncoKHx09epTu3buzefPm+isRqvEvqFR5JcjPzxfDhw9vOXr06NQHH3wwDcDPz68wJiZGB6bVsa+vb4UrqtDQUH3xFXNMTIxjYGCgvrCwkLZt20a0bds2YsqUKU3N98eMGRPaokWLvNdee83yp29gYGDBmTNnHMFkwLOysrSNGjUqO6KoiIEDB2bv3r372MGDB4/069cvKywsLK8yc7+igZZSPlzGq07sWBiNRqKjo+nZsyfXX389N954I6+99hpDhw7l119/paxdWnvQsmVLvLy8AGsh/OLvr5RYQKPR8PLLL3PmzBny8vI4ePAgPj4+lnvmBLiurq4WjYfyVgAuLi5XP5kGyLrl+xg3eCbDOv4f4wbPLNdHXJYeB0CmTGMLKyzXHbmR64K7061TE4wnzmA4fBzD8dMY0zIqLXz0zz//0KNHDzIyMli/fj3jxo27cqW6iNtTceBU4gfOyWgqv3qMRiP33ntvSHh4eN6MGTMssRRDhgxJ++qrr/wAvvrqK7+hQ4emVdROSEiI3s3Nzbh27Vo3o9HIggUL/G6//fY0BwcHjh49Gn306NHojz/++ALAM8880zQjI0M7d+5cq0CJ4cOHp82bN88P4Ntvv/W58cYbM6+U7i0uLs4BTCdAZs6cGThx4sRKKYBe0UALIb4VQswr+apMJ/bixIkT5ObmljJ+Y8aMISkpiR07dtTIODQaDffccw8AR44cYcGCBSQmJvLRR6a0jTqdjtGjR1ue//jjj0tl/l6zZg0bNmwgIyODixcv8tVXXzFr1iwARo0aRbNmlg1j7r3XFCeUmprK/PnzycrK4vPPL0el9e3b167zrU+sW76P2TP+ICk+DSkhKT6ND57+jjsbP1pKVrQs3Y1UmcR21liuu9KX5m5hdB/elX++WIEsKHI56QuR8Um0bezE3Ke+skmy9MyZM9x6662EhISwY8eO+ufWKIbGbUwqHi/HoAkoAAGagAI8Xo7RuI1JrUq7q1evdv/jjz/8Nm/e7GFe6S5atMjrjTfeiF+/fr1nSEhI+/Xr13u+8cYblt3xoKCgDq+++mrzxYsX+zVu3Ljj7t27nQE+//zzmIkTJ4aGhIS0Dw0NzR89enSp02WnTp3Sffrpp01OnDjhHBkZGdG2bduIjz76yB9g8uTJKZcuXXIIDg5u/+mnnwZ+8MEHlpMiXbt2bfPAAw+02LZtm2fjxo07/vbbb54Ab775ZmCLFi0i27VrFzls2LC02267rVL7eLb4oP8q9t4ZGAlcqEwn9sLsfy6ZVWXYsGFoNBpWrlx5VedHN2zYwM0338zatWttNnZvvPEGy5cvJy4ujvvvv9/q3uuvv37FVPdRUVG8/fbbpcpbtWrFp59+alX23HPPsWjRIo4cOcKDDz7Igw8+aLk3aNAgK6Gea4F1y/fx3ezVJCekERDozUOTB1nSR303e7WVtKcxLQPDhUT0Rf7E4gJFAc39rHIHJshzHGK75fpGMZiw4BY88s6YslfbRiN7VuwpihgoX/zIfEQvLCyMn3/+mSFDhuDu7l5t/x+1hcZtTCpVNMglGTJkSJaUcndZ97Zt21Zm3ra4uLiDZZXfdNNNOSdOnDhcUX8tW7bUl9efq6urXLFixemy7u3evbt0MlXgq6++Og+UOvJnK7a4OH4r9loA3I1J3arWOXjwIM7OzlZH7AD8/f3p3r27VQBBZbiazN+BgYFs3bqVMWPG4Ofnh5OTEx07duSbb765YhQhQPfu3enRowc+Pj44OjoSFhbGs88+y/bt2wkMDLR61t3dnQ0bNjB+/HgaN26MTqejRYsWvPLKKyxbtqxWBKBqi7JWyMWPupXcsJNJKVDiczXLij7yzhhwMP1IxMoTVsZ56Jh32Wr8x2Kcixty6w6sL0tKlqakpDB48GDWrVsHmM69NwTjrLAP5Wb1LreCSYh/uZSy1RUfriSVzepdX7kWxJJqKqv3uMEzSTgSazK8+kLQOSAa+RPYLpj5q6YxbvBMkuLTLM8bDpedLFcIwSrDL/S94SXO7l1KbOERy73uw99A5+TKjDu7lwrNtgVz24cPH+a2224jLi6O+fPnc/fdd1/VnO1BWVm99+/ff7ZTp07l/CZSVBf79+/379SpU2hZ92zxQWcKITLML2AZ8GI1j1GhuCoSjsYiLySajDOYfMEXEkk4agqtf2jyIJycix131ZXt1Qto7se65fs4fnKZlXG+8bZ3cHB0pZG/Z7mbiFcioLkfy5cv58YbbyQnJ4cNGzbUKeOsqLvY4uLwkFJ6FnuFSyl/u1I9haImECmppVwWSGkqh1IZTLzCg3FwsjbSTq6O9BjdizvvHk7CpWhL+aD2L+GYbcTJyYEJY3tflXi/k6sjPR7uxIgRI2jVqhU7duywaLQoFFfCFj3otVLKm69UplDUBsb8sle0xcvNGUzMlAxGefjt+xj+0FAKCk0h9xrhwMCIaQgELmn5TH1tJINviuD7EpuIV6JRsD+PvDOGAff1RvoU8Oijj+Lm5nZ1E1Vck5RroIUQzoAr4C+E8AHMO0+eQNPy6ikUNUmj5v5lGs1GzcvXpCkuK2o0Gq1EtdydAujVulgUZoGBwTdFAGXLg5aHngJyuiQQ1iMIIQTPPPOMrVNSKCxU5OKYAOwG2hb9a379CXxWQT2FosboPrxrpcqLk5+fb2WcAzxaWxtnrKMCi8uDApeXLCXIkuns0qxj5cqVREdHl/2QwiYqKzf6+++/e0ZGRrYLDw+PiIyMbLd06VJLaKa95Ea3bt3q0rlz57atWrWKDA8Pj/j6668tIkx2kxuVUs6WUoYBz0spW0gpw4penaSU/6tULwqFndi+vMwjq+WWm0lPT8fZ2dlyPXLEfdzYZozVM2VFBZrlQVcbf+WlH57Bw9f6iFyyvMBO1uPs4cTGjRu59dY6o4pgd06lL/ZddnZIh19Pde267OyQDqfSF9e43GijRo30y5cvP3n8+PHo77777sxjjz1mUZuyl9you7u78Ycffjhz8uTJw6tWrToxffr05ikpKVq7y40CSCk/FUK0F0LcLYQYZ37Z0rgQwlsIsVgIcVQIcUQIcWNlBqdQXImrycYdHx+Pt7e35fr9999nydKFVpuJjZp4M3nGHVa+65LcPKYPS1K+5aUfn8GzkTvJxLOfrbRs0Yr9h/ZbZfJu6JxKX+y77+KHIXkGkzhSniHFcd/FD0OqaqQrKzfaq1ev3NDQUD1A165d8woKCjS5ubnCnnKjHTt2zO/QoUM+mDQ/fH19C+Pj4x2qQ27Ulk3C14F+QATwN6Ycg5uB+Ta0PxtYKaW8qyh5rGtlBqdQXImS0X/Fy8vi6NGjVskdfvzxR8aOHQuU3kw0a3gUj1CU6ZlWG4yPvDOGm8f0QdvUwA033MA777zD9OnTcXW9tr7q0Ze+DjLKAqsFn1EWaKIvfR3U0uuuaokurKzc6Pfff+8TERGR4+LiIqsiNzpt2rREKF9utEmTJhahpvXr17vq9XoRERGRL4TALDd600035dhLbvQuoBOwV0r5sBCiMfDNlSoJITyBm4CHAKSUBcBVKVt5e3s3qECVhISEUtGBDY3iK9TqpGRYd4/RvfjnixUV5vUzs23bNiu9i1WrVjFoUNnCRuYIRXOYuFnDQ8YnUphv+nlMik3h3cdn8cGcd3n6hadwdXWtURXFuoR55WxreWWprNzorl27nF977bWglStXnoCakRuNiYnRPfzwwy3mzp17xry3UVW5UVsMdK6U0iiEKCwyuklACxvqtQCSgW+FEJ0wbTBOllJmF39ICDEeGA8mSc6yspiU1Nqo72RlZV0T4b1ZWVk2Z6WxhcM741nxUzSFetPPZ1J8Gv+kZNB+REdOrD9KRnIWngHu9H/sRrRNDVZ9b9261Srk/ssvv0Sn05U7vi/f22Sl4QGgjysWEINJ5W5/zlb0G/LpO7TPNbdqLo6z1r+gLGPsrLWv3GhISIi+pNzoqVOndHfddVeruXPnnomMjLS4HsqTG23fvn0EwNChQ9PMinYVyY22bNlSX1JuNDU1VTNs2LBWr732WtzNN99ssXFmuVGAJUuWeJ48efLyxocN2GKgdwkhvIGvMRnZLMAWmTgH4DrgaSnldiHEbOAl4NXiD0kp5wBzwBTqbe/Q4LpAVFSU3UOg6wLVOc91y/ex/Mc1GI3Wq5hCvZELFzX8lvBtuXVL6qH88M0yVi06wh8Jq0uJK5l5L62MrCjFjHOSjOMwO3DAketFf3r06HFNfKblEeHzeNy+ix+GFHdzaISjMcLncbvKjb7zzjsJxeVGU1JStLfcckvrGTNmnB88eLDFUBaXG+3fv3/2ggUL/J566qkks9xo8T7NcqM///zz2eLlZrnRgQMHZheXG83LyxPDhw9vde+991585JFHLhWvExcX5xAUFFRolht9+eWXK5WT7ko5CQXwXyllmpTyS0wZuh+UUj5sQ9vngfNSSrPizGJMBluhqBRmd0NJ42zGIoi0YAGEhoJGY/p3wQLefPNNq+wyi39cy69f7ipXXMlMmaL7RWHiF+RZDrANNzy5gQG0DK52WZp6R0uvu1I7+z0XY14xO2v9Czr7PRdTVf9zZeVG33///UaxsbFO7777blPz82ZNZnvJjc6bN89n586d7gsXLvQ397l161YXqLrc6BXFkopEVK58qLTsupuAx6SUx4QQMwA3KeW08p4vTyypoaFW0JWjpOBRSRo18Wb+g01h/Hgolo/xMa2WuYbLCS8yMzN5ctQXZbbVqIk381dd/mqW9EEDaLKzkfGJZOVlcJZjtKYjrm4uTJ0zEW1TQ73+TJVYUu1RJbEk4F8hxNWeF3oaWCCEOAB0Bt65ynYU1zAV5fiznFX+v/+zMs4Dwco4FxQU4O7uXm5bJctLani4+2rwvSGbZ7+eSHBIMO0019EkpDFT50y00npWKKoTW3zQ/YGJQoizmJLGCkBKKTteqaKUch9w/ZWeUygqIiDQu8xVr0YjLp9VHhFrKQ8FYoo9N/LxLxlw72wa+Xvi7u1K5qXSWc/LcmmYj93t2bOH22+/ndRdqbz0ymQWnP2iqlNSKGzClhX0MEwnMsxJY28t+lehqBFKSYZicje4JMTz7m3vMDb0CXL9GiExrR7Mxrk5EO/qQ1JKJlJCYnIGl9y0kJmF4fhpSy5BTXZ2uXkEf/vtN3r37o0Qgi1bttC+fZ3IVaG4RrAlkjAG03d9QNH7HFvqKRQVsWpjNHdOmEOfOz/gzglzWLWxfM2KUu4GCpHxiWQkZ1iyb89OD7P6UvYDjjs48mXEUKu29ClpGOISrPWj4xOR6aX3bmbNmsVdd91F586d2blzZ4M77qmo+9gi2P86JoH+l4uKdMCP9hyUomGzamM0732xisTkDMvK9r0vVl3RSM9fNY0VB97GKTPDEiwCYJCFvKz/13L9MLA+JIT3Ot/JmubWB4ccos+BwTrOoTC/0CotlZlu3brxyCOPsG7dOho3bnyVs1Uorh5bVsIjgdsw+Z+RUl4APCqsoVBUwFcLNpNfzMAC5OcX8tWCzTbVL66zUSDzWc8fluvXXnuNeVLC2bMcuK4fDhn5uJ5Jx+3EJVzPpCNyy46bMLcZFxfH119/DUDv3r2ZO3eulaiSQlGT2GKgC6TpLJ4EEEIoxXFFlUhKySizPDE5wyaXh1lnI1dms5FllvII3Q3s+L2QcYNnsm75Pm5qHYRTUg6aQiMCTP9WkPJq165ddOvWjeeee474+ErFEyjsRGXlRtevX+9qPovcpk2biPnz53ub27KX3Ojx48cdIyMj27Vt2zaiVatWke+//36AuY7d5EaL8YsQ4ivAWwjxOLAGU1ShQlFp1i7chMvq/Tj/sR2nf/aiPWd9zNYWl8cj74whzymLLaywlHUUPQlqfINV8MnmP/YiShzzF438oaQGg0ZD2IjG9OnTBycnJ7Zs2UKTJk2qZb7XEmsTo3yf3vNchwe2P9b16T3PdVibGFXjcqPXX3993sGDB6OPHj0avWrVqhNTpkwJMRtFe8mNBgcH63ft2nX06NGj0bt37z4ye/bswLNnz+pqSm70A0xRgL8B4cBrUspPK9OJQgEm4zxr/JfIrDzTija3AN2+M6WMNFTs8pCN8tmct9Jy3c1xEI2DrkPjffnnLT9PT2Z6bqm6Gm9PRNPGl5PH6hw44xbD2/97g+uvv54dO3bQoUOHqk30GmRtYpTvgphFIWn6dEeANH2644KYRSFVNdKVlRv18PAw6nSmEz+5ubnCLIhkT7lRZ2dn6eLiIs19mlfm1SE3autpjIPAJmBj0XuFotKUlRVbGIw4Hoix+IgdMvIt98pyhUx7e6aVAt0Xvy7Dr003K+N8JTTenmjDW6CNDEcb3oKAZs15+OGHWbNmDQEBAVduQFGKP+KWBeml3sqe6KVe80fcsqDy6lQWW+VG161b59aqVavI6667LnLWrFkxOp2OqsiNDhs2LAPKlxsFkysmPDw8IiwsrOMzzzyTEBoaqg8MDCw0y40CXI3cqC2nOB7DJI40CpP06L9CiEcq04ni2mDtwk2MDX2Cwdq7GRv6BIfWHLe6X66Ivr7Q4iN2SsyxGOlG/tZG99GpL/DBKy9Yrts8/hpzd57B3bfsbREhBFqHsr/iefoMUjJP4+Ss4+33/4+5c+fi5ORk40wVJTGvnG0tryyVkRsdMGBA9smTJw9v3rz5yMyZM5vk5OSUKWlRnXKjrVq10h8/fjz6yJEjhxYuXOh/7tw5B41GY5Eb7dChQzsPDw9DZeVGbVlBTwO6SCkfklI+CHTFdOxOobBgdl8kxaZYziYv/3A9axdusjxTMj2UhWIbdwJwTMzBycmBCWN7W8qfffZZ5n0803Ld7om30bl7k6cvJL2VS6lAFjD9QJX1Q5iWE8e2U99yMG4ZE14awoDhnSv8YVVcGW+dV5nHY8orrwwVyY2CyX1RXG7UzHXXXZfn6upq2LVrl0tFcqPmTcUpU6ZYkmFXJDcKJgNeXG7UTGhoqL5Nmza5a9as8YDLcqMHDx480q9fv6ywsLC8yszdFgN9Hih+ij8TOFeZThQNn7LcF8XPF69duInMtOzSFUXRxl0xNMBtXVtbsmmPGjWKWbNmWe5HPP0+WmcXy3WCJ0yecQcaTWkjW6g3WJVfSDvMzjM/ohU6hnSbwC13qixs1cEdQSPidEJntbLVCZ3xjqARdpUbBSguN3r06FFH86bg8ePHHc+cOePcunXrguJyo0ajkQULFvjdfvvtaWa50aNHj0abtaDNcqNz5861snNmuVGA4nKjp06d0mVlZQmA5ORk7a5du9wjIyPzwCQ3Cibf9MyZMwMnTpyYXJn526LFEQdsF0L8iemo3e3ADiHEswBSyo8q06GiYVJRbsC1Czfx3rhPkWXJhQpRpv9496poeH44Xbp0Yd++fZbyyMkfIDTW64pAHw8GDO/MzJd/LXMMRqPE0cmBw7FrOJ28FR/X5tzQ+l6m/F/prCuKq+Pmxv1SweSLTtOnO3rrvAruCBoRZy6/Wsxyo61bt85t27ZtBMAbb7wR98Ybb8SPHDmyZUhIiH/Tpk0L/vjjj1MAa9eudb/11lubODg4SI1GIz/88MNYc0qqzz//PObRRx8Ny8vLE/3798+oSG40LCwsLzIyMgJg/PjxSc8++2zK5MmTU+68886w4ODg9l5eXoZFixadAjhw4IDLiy++2EwIgZSSSZMmJdxwww25YJIbXb16tZfRaBSPPPJIkj3kRl+v6L6U8o3KdFgRSm60/jI29IkycwN6+LpTkFdQanVdHG1keBmlkk2nZ5NTpFCndXCi+63/Aa0g11egdzcZaWedAzPuHsjwru0Y3futMk9ueHi58MTLt/LkE5PIzsrmpuvu5tGpwypMCFtZ6vtnquRGa4+K5EavuIKuTgOsaLg88s4YZo3/spQhzkzNqriizoF+aYd4ODGKAH0GyTpP5jbqy//FXQ5AcfMOonP/KaYLI7helORgxL+5F5Nv6cXwrqYksGX5kXML0pE5GQwY3pnosxsRQih/s6LeYEtW7+uB/wNCij9vi9yo4trBrIk8b/pC00paUBR7Wj5Sq2FgQAZT4lbhLE17PF76DCvj7Ne0A227jytREcKEB7+9+phVcWa6tYzopZzz7ItZjKODG0bjbDQapfGlqF/Y4oNegOkkx0Hgitl0i1OkIZ0JGIDCkn9CKRoWN4/pw81j+pTr7iiOFCC7tWTy6e8sxjkNKB450LTVTYR1KFvZtqwz0sV1o+MuHeTwhb9x1nky4Lpxyjgr6iW2GOhkKeXSKvTRX0qp/FgNlOW7jzD77y0kXMok0MeDIcLdJuPs0q89T78+Gpd+7wKmnehmxZ75CNg8YCyJyWXrdjTy92TVxmi+WrCZpJQMGvl7ctPgCNb+vJ2DMas5m/Ivvm4h3ND6biarzUBFPcWWZcXrQohvhBD3CSFGmV92H5mizrN89xFm/LKG+EuZSODSttOsePP3K9YTQjDl0QGmY3TBwURjbZx/AqaGhDBhbG+02tL+Yp2Dhp5dw0pJli7dfYJ+o7uSZ0yhuW8XhvWYyLS3xlTrZqBCUZPYsoJ+GGiLSQfa7OKQwBIb6kpgVVG0zVdSyjklHxBCjAfGAwQEBBAVFWVDs/WbrKysBjHP91YfIk9/OT7AfUssotAGL5hR8vmz89A2NXC+b18emD/fcmsN0M/Jiej778fRmMTIgaEsj4olN98UD+DqrOWWvsGs3nLESrI0LzuVgnxH1kbH8P1Pn+HoaA5gS6uR/+uG8pkq6ha2GOhOUsqrVY/pJaW8IIRoBKwWQhyVUm4s/kCR0Z4DpmN29fmokq3U9yNZZtKW7rW61mbml/NkaTKTs0lLS7MyznuBziEh8PbbRIwdSwTQrx9MmVC6/m+rPrC8z0g5w5Ht3+PhE4xTr0cYPHhwJWdSdRrKZ1rXOHnypG7s2LFhycnJOo1Gw4MPPpj86quvJiUmJmpHjhzZIi4uzikoKCj/zz//PB0QEGCJ6jtx4oRjp06dIp9//vkLb775ZiKY5EYfffTR0Ly8PM2AAQPS582bd67k3kRmZqZmxIgRLWJiYpy0Wi2DBw9O+/zzz+PAFGxy1113hR08eNDV29u78Ndffz3dpk0by7Gl1NRUTdu2bdsPHTo0bf78+bFgkht96aWXmun1ek2HDh2yFy1adNYs5mQLtmb1jrC5xWIUifsjpUwCfgduuJp2FHUTL1drIXuDh+1aFmneCYwcOdJy3Xvqfxk39SMGP/oKy9teV0FNE2adjsSYnRza/BUOjq6EdRhRSr9DUXP8dHqXb5+/P+rQbsmbXfv8/VGHn07vqnG5UTOTJk1q3rdvX6tAFHvJjRarF9S9e3dLIEqNyI0CvYF9QohjQogDQoiDQogDV6okhHATQniY3wODgUOVGZyibiNLnKPL6hWMsRxxouKcdTjKztTLGh2dJr1DGk5IIP5SJjN+WcPy3UcqbOOxe28kNno5J/f8gldACzr1nYS3fxMr/Q5FzfHT6V2+7x5YFZKcl+UogeS8LMd3D6wKqaqRrqzcKMAPP/zgHRoamt+uXTuL7oU95UbBtDpPTk7WDRo0yLKrXVNyo0OB1pgMbGWyejcGNgsh9mNSw1supVx5hTqKekRGjrVLI69tABkDW2LwcEIIgWcjd1768Rle+vEZGgX7I4TglOsBTuov/54e8Mr/MOisV+J5+kJm/72lwsSyN3ZpTmHWKVpG9iey56MEBQXy4hODLfodiprl86Mbg/KNhVb2JN9YqPn86MYalRvNyMjQfPjhh4Hvv//+heJ17Sk3ajAYeO6555p//PHHVtod1SE3akskYYwQojfQWkr5rRAiAChHlsyq3mmgU2UGo6gbrFu+j+9mryY5IY2AQG8emjyIAcM7lyoPDHYgvsQfiXltA/C5sQVPje7FW+vWMD55J008PJi28hk+mjCRMxsvS5D2uuNdco/movO5HLqtyzLitj8Zw6FdvP/6P0gXRzQRzUmUmE5tJFzgnjtuwtvbm+hDB/D29q7B/xlFeaTkZZVpeMorryy2yo0+//zzTSdNmpTo5eVl9Yw95Ubfe++9gMGDB6e1atXKKp9VcbnRgoICTf/+/dMrKzdqSyTh68D1QBvgWy5n9e5VqZ4U9YJ1y/cxe8Yf5OeZvmvm9FHRe2NY/edeq3JNgsC3kQt6Hyfyioyss86BXp3DmL52NbmFplMWFzIzuevGGylMu+wS7HnH+yahJAO4XJRo8wy4xOfjdDwRzieYcl8BoijrCkCSczqPjXuF3RMf4+OPP1bGuQ7h7+xekFyGMfZ3drer3GhISIi+uNzo7t273ZYvX+7z+uuvN8vIyNBqNBqcnZ2NY8eOTStPbrR9+/YRAEOHDk0zK9pVJDfasmVLfXG50X///dd9586d7t9++22jnJwcjV6v17i7uxs+//zzOLPcKMCSJUs8T548WakMxLac4hgJdAH2gGnjz+xbVjQ8vpu92mKEzeTn6Vn+y85SKwjjpXQcj51Cpy/ExcUReVNLnnvrPv67f7PFOEspOTv1eUsdNw9/Og+0lhMXElwv5OOclIMxIdlinC33DUYSD6znuH4Xzu7+TJo0qTqnrKgGnmx7U9y7B1aFFHdzOGkcjE+2vcmucqPvvPNOQnG5UbMxBHj22Weburu7G6ZPn54MYJYb7d+/f/aCBQv8nnrqqSSz3GjxPs1yoz///PPZ4uVmudGBAwdmF5cbXbp06RnzM5988onfrl273MwnP+Li4hyCgoIKzXKjL7/8cqWyEdtioAuklNKcOUBl9W7YJCeklVleyjinZSAvJIKUCEwrXdYcw/mBFOIzTRvZ0mDg7HOXM6C4tAmnS+TjJe0vAE4X80wJXvXWuutGaeQEBzhXcBLvxm3oM+wJWrVqVZUpKuzAfS2uTwWTLzolL8vR39m94Mm2N8WZy6+WysqNVoS95EYroibkRp/HtEk4CPgv8Ajwk5Tyk8p0ZAtKbrT2GTd4pkXPoiIMx0+XMqYAQiMwGiX53oLNly7rM3vc2IOOjz1C4+2yzPBttxOXEGW0myOz2M4amjq0pOWox3n5qWF1ciOwLn+mtqDkRmuPiuRGK5vVuw2mrN7VbpwVdYOHJg8q957Vr/IyjDOANEr0Ms/KOHvfMpTmY+9jWs/eTBjbGyen0n+4yaLjeaKRPwhBvsxDSomrcKeHdiidbrqvzhpnhcJe2JI09j0p5Wop5TQp5fNSytVCiPdqYnCKmmfA8M54erte+UFd2d6xHJnFRv6yXLdx7k6zNoPxSdRxbHscXy3YTH5+oSUNVeMAT0YO7YQx0A0pTBm3L/nq2cYqznMaF2833pj/AsvWvq6Ms+Kaw5Zz0GUtqYZV90AUdYeJLw0vlYS1pCPMvNItTp7MYSuXj7p3phfN8psDkJOQw59L91ncG0ajtCSGfW78IF54+TacWvkRk7aXPSl/4urqzSvfvMnS1O8sWtMKxbVGuQZaCPGEEOIg0KYogtD8OgNcMZJQUX8ZMLwzk2fcgYuro8UwC6w1+DXenoimjemvOc+P8m8+k4uJZYWljW70x180sYR/O1+Spk3AYuTnF/LVgs2mPnuGE9z0NEfPr2DYLcM4F3+C+x+91a7zVCjqOhWd4lgIrMC0MfhSsfJMKWWVdmYVdYc/jx5h5tbNxGdmmgJKevbm9rbtGDC8Mx9MX0zJo/zFjfQAYpli3MUxChkCJCLpCgymCzuFHxLID/M21TNQJmbh/W3btvH555/z3HPP8d5771HZA/0KRUOkXAMtpUwH0oH7am44iprkz6NHSgWUvP7jSj49s4q8XD1uRlnKQAOWsocTo9gtCxmO6YtyMyZFrByOcT8tEYDTmTQyMYn0l1xBA7g4mVrr06cPBw4coH379tU9TYWi3qLyAF3DzNy6mYK0QlzOg2sMeB4x4hptIC9Xf+XKwE59BoMwGedRwHLAAwjgcm5AbWY+LvF5uJ9Ox+3EJVzPpOOQYdLwSEs6zoYlM1izZg2AMs6KUpw8eVLXvXv38BYtWkS2atUq8j//+U8jgMTERG3Pnj1bh4SEtO/Zs2fr5ORkLZj0Opydna9r27ZtRNu2bSPGjBkTbG5r06ZNruHh4RHBwcHtH3rooeZmoaPiZGZmavr169cqLCwsslWrVpFPPvmkRUskNzdXDB8+vEVwcHD7jh07tj127JglclKr1XY19zlgwADLQf2lS5d6REREtGvdunXkqFGjQvV62362zCgDfQ2TFJ+JUypoDAKBwPWi9RfC4KwpN+9rfNphRgK5wKPAIsAsNprM5VMgws0J7+gsNIVGBKApNOKYmE3i4Q0c3joXByd3FXjSQFhwYL9vj2++6tBy9kdde3zzVYcFB/bXitxo8+bN848ePRp99OjR6IULF8aay+0pN+rk5GQ097lu3bqTUHNyowghQoQQA4veu6hQ74aBc4YGIS87MYr7iR0y8tHmG8t0ccRe3M2B839SCExF8DWXfWV5aJmLaSVsdNDgFRqEpthCxSgNHL2wkpPH/8K3cVtuHvUyoaGh1TwzRU2z4MB+37c2RoUkZWc7SiApO9vxrY1RIVU10lcjN1oW9pYbLYsakRsVQjyOKVDlq6KiZsAflelEUbNUJNNZHKm/vD52ic/D7cxlN4Rjcm4pn7GUkpNJmzgS/w8ArbSd8KcbScIVIxDv7MkHzt1ZJ4Ip9HDCcEs42cJ6myMh/QjnUvcS5n8jnfo+yqSHB1br3BW1w6c7/g3KNxis5UYNBs2nO/6tUblRgPPnzzu2a9cuolu3bm1WrlzpDvaVGwUoKCjQtG/fvl2nTp3a/vDDD95QQ3KjwFOYMqFsB5BSnihKYaWog6zaGM17X6yy5OtLTM7gvS9WAZQK9PCXjuSfz0V3ySRUZDbIotBYyrUhpeRowmpiL+4CBO3EdQQZw1gvYD3BGHUaMvq3JK9tAADOOgfu7hDEwfmnSIpPwygNaISWJl6ROOs88fMNY+qTQ1XwSQMhOTu7TMNTXnllsVVuNDg4WH/mzJkDgYGBhk2bNrmOHj26VXR09CF7yo0CnDx58kBoaKg+OjracdCgQW2uu+663MjIyHy7y40C+VLKAvNkhBAOlI5bKBchhBbYBcRJKdXBVjtjjtQrjvm8cW4jcflIXboTmvN5aAzFhIqKUfyra5QGDsct50LaITRCS0+/AThf9LJ6XlNoxOufE3itPIHRw4nIh3rTrZkvXSe3YPrUmRw4+zddQ+/FzcmXJn4tmTzjDgYo49xgCHBzK0gqwxgHuLnVqNyoi4uLdHFxMQD06dMnJzg4OP/QoUPOoaGhenvJjQKEhobqASIiIgp69OiRuWPHDtfIyMj8qsqN2uKD3iCEmA64CCEGAb8CyyrRx2Sg4vxFimrDfK64JInJGUxfu5oLmZlIoPBILkZDkeZyOZm4JWAw6tkXu4QLaYfQanS8+5//4ZLqXebzQpoMuzYznyNfrGfON1t44n//ZfvJn3DQOaERWho18TYZ5+GdqzxXRd3h6Rt6xDlptVZfJCet1vj0DT3sKjcKUFxu9MKFCw6FRcdGo6OjHc+ePevUpk2b/JCQEL1ZbtRoNLJgwQK/22+/Pc0sN3r06NFos3E2y43OnTvXKkOKWW4UoLjcaHJysjY3N1cAxMfHO+zatcu9Y8eOuWCSGwXTCZCZM2cGTpw4Mbky87dlBf0Spo36g8AE4G/gG1saF0I0A4YDbwPPVmZgiqujkb9nmWpxuAjLeWcAbbFsVcaMTEhINgkg6RwQjfzReHtSIPLYF7uYtKxYHLQutBrwID9k5BHs70Z+claF45D6Qtb/8i0J+lN4tupA0JAx5Li5cc/dAxnQtV11TVdRRxjbsVMqmHzRydnZjgFubgVP39Ajzlx+tVRWbnTVqlXub731VpBWq5VarVZ+/PHHMY0bNzaA/eRG9+3b5/zUU0+FCCGQUjJlypSErl275kHNyI26AXlSSkPRtRZwklLmVFjR9OxiTJGIHsDzZbk4hBDjgfEAAQEBXX/55ZfKjL9ekpWVhbv7FbOGVZp9xy6yPCqW3HzrsD2dg4aUZkaMWg3CAFILnrFGtHrQnktBt/c0wljseyAEBYEe7MlbSVbqBRzcPAkdNRFnf9NJJqejyfiuO40sKCc8EDgtj3CawwTcMJBGPYcihOmPNR8XHW8Oanjnne31mdYU/fv3V3KjtURFcqO2rKDXAgMB85LJBVgF9KyokhDiViBJSrlbCNGvvOeklHOAOWDSg67Pmrq2Yg/t4FUbo1kWtY/8EsZZCOjQsRnrL8WhKbolDJDnBa4XwSH6nLVxBnKNWeyNX0mOzMTFyYfW/R/D6H85q31+2wAuCWi+N4mMhHQoFiUopUQIQQjhuLr4Q69brNpOy9XXa93k8qjvetCKuoktPmhnKaXl79mi9zboUdILuE0IcRb4GRgghPjxqkapuCJlbQ6CKXvUnr2xOCcb8DhnwPOs6V+BINdPmDKhFCNLZrCLKHJkJh7Ojeke9gCNY3W4xOdZPZfbxp9Hlz3H1J1voh/RDqODhhQZzy7WUyj1CJ0O536lf4cH+qgj9AqFrdhioLOFENeZL4QQXTEFkFWIlPJlKWUzKWUocC+wTkp5/1WPVFEh5W0OAiDBMQs0BtMmnsYALikS54sS6XJ54z1dprKbKPLJxVv40y1sLE4O7miM4HnS2qMlEMz+ewvDu7Yj6vc3cB2gYx9bMGLE4Kkl4on+GCKtT2M66xyYfIvKNaxQ2IotLo4pwK9CiAtF102Ae+w2IsVVUe7mYBFlqdINit3DxMLlNJIZLMGJsejRY8SfJnRschs67eUTQdo8I7oso0k2tMiPnZKVTvtnZnJp8x9c2LuV0aNH89133+HqavoDa+aPS1h9JpWES5kE+ngw+ZZeDFcbhAqFzVzRQEspdwoh2mJKdyWAo1LKSil+SCmjgKirGaDCNiaM7W0VoHIlBp7bw4t7F+Ni0LMEGEs+BUAn4Y9/0ztw8LGOgjU6anC5eFnTWRjA5aLk3JHfSTmyjaa9hjHuhRkW4wzQrZkv0+4fVU0zVCiuPWwVS2oDRABdgPuEEOPsNyTF1TD4pghefGIwnh62nYOfeHgFLgY9c4HRQAHwNLBCm4+2pHHWQL6fc+lgFgnBYQNofutD+N4wiE9WbK2OqSgUiiJs0eJ4Hfi06NUfeB+4zc7jUlwFg2+K4O/vJ/HalFtoHOCJpPyQz0a5acwEHgOMwAxgNtC4MNPiDpGAwQHSItwxuDtZ6l6MP8yxnQuQ0oizkzderTsCkHCpUkc8FYorUlm5UYDt27e7dO7cuW2rVq0iw8PDI3JycgTYV270xIkTjr169WrdokWLyJYtW0aa79WE3OhdmLTYE6SUDwOduKwsqagj/Hn0CL3nfU3L2R/xzrJ/EMdTcD9xCdeYdLSZ+VbPSimZ7ODEC0XXnwCvY/JfJesuKzAKTNm2c5s4Y3Qw1Tt/bB1H//2evKwUCvV5yGLSAuqExrXNoi37ffvPmNOh47OzuvafMafDoi01Lzeq1+t54IEHwr744ouYkydPHt64ceMxR0dHCfaVGx07dmzY888/n3j69OnDe/bsOdK0adPCmpIbzZVSGoFCIYQnkAS0qEwnCvtizoxyITMT95N5OK5OISsl2xR2XWDEOSkHTVa+KcRbY+Toqd/4X2E+DsCPmFwbAHnCgW8b97NqW5tnpOnqFFxOX+TkvwuJiV6Bf7NOtL/pSRycXMnzMa231QmNa5tFW/b7vv/nhpCUDJPcaEpGtuP7f24IqaqRrqzc6JIlS7zatWuXe+ONN+YCBAYGGhwcHOwqN7p7925ng8HAyJEjMwC8vLyMHh4exhqRGwV2CSG8ga+B3cAeYEdlOlFUP8tiDzJgxWzaLXmT59eusIRx++/KtQSkmBFGcLyUR2qIkUNHF5B6cDsajQOjW93EzS7eGIFEnRcfB91ClLd1lJ85Wezh40tISthHm3ZDadNtDD4+njgGu1DorqGJjwcz7h6oTmhcw3y5entQQaG13GhBoUHz5ertNSo3euzYMSchBL17924dERHR7pVXXmkM9pUbjY6Odvb09DQMHjy4Zbt27SImTJjQrLCw0L5yo0KIXlLKLcBUKWU+8KUQYiXgKaVUWb1rkWWxB3l171/kGUzfN0OhxHyQziGrHOGj7DzOrviFvGPH0eqcadfjYWL9W3BXhxGmehn5OBWTHC1JWEAPmni3p2OzHsxfMq3a56So31zMKFtWtLzyymKr3GhhYaHYuXOn+65du464u7sb+/TpE96tW7ccb2/vUroE1SU3WlhYKHbt2uW+ffv26NatWxfceuutLT/99FP/qVOnplRVbrSiFfQnRf9uMxdIKc8q41z7zDq8jjyDHpHkgHa7G64xApfzoM2CQvfSH2lBYQ67YheSd+w4Gnd32veeiJe/tZeq0NOJvEau6N01ls3FxPSjnEraDIC3azMCvdqSFJ9m/wkq6h1+nmXLipZXXhkqkhsF0+rYLDfarFmzgh49emQ2adKk0MPDwzho0KD0Xbt2uVYkN2rOJThlypSm5vsVyY2CyYCb5UaDg4ML2rVrlxsREVGg0+m47bbbLu3Zs8cVwCw3evDgwSP9+vXLCgsLsw7JvQIVGWi9EOJboJkQ4pOSr8p0oqhe4nPTTcb5hDMavQaBQGMQOKVCansXjMV+SefpM9lx9kcysi/g4OND08mTcG5c9l+d+QFOxN7vw/FHfTl5cQv7zi0hOfMURuPls9XCqXIrAMW1wcRB3eMcHazlRh0dtMaJg7rXqNzoyJEjM44cOeKSmZmp0ev1bNmyxSMyMjLPnnKjffv2zU5PT9deuHDBAWD9+vWeERERdpcbvRWTSNIATL5nRR3BJ82DjOMSUSI+UEhBgbczyTcJWh6UJF44x66Yn8jLT8cxMJDGE8cTENiYTGMOXsexzhWoAYd2LrzdqyePP/YYGfF7aOIdSWTT4Wg0pq+JFJDnow7wKEpzTy+TrOiXq7cHXczIdvTzdCuYOKh7nLn8aqms3GhAQIBh0qRJiV26dGknhODmm29Ov/fee9PBfnKjDg4OvPvuu+f79esXDtChQ4ecqVOnpoCd5UaLpEUnSyk/qkyjV0ubNm3ksWPHaqKrWqUqymfLdx/h1UWr0Jcrsi/JDYFbclz56o1p6POzcPMLpvndj6Jxc7eEaWsMRtziTbrQBifIa6nl9TGDmTXxCTZt2kRIp2GE+PXEKTXflALLQUOBnzN+LQP47avxdp9nfaO+z1UIoeRGa4mrlhuVUhqEECOAGjHQCmuW7z7C7L+3EH8pE40z6INz0cW6YKwgmltqQe46xRc/z8NQmId3o9a07f4gmkQncv0kevciTWih4WJ7MLiDVgg+GDyU29u2I3v8eJ599llcfcN574tV5Hhdjkx0cnJgwtjeNTBzhUIBtoklbRVC/A9YBGSbC6WUe+w2KgXLdx9hxi9ryNObrLExDzQnnDHI0q4NMxLJpYRoEn76HmkoxK9pR8Kvvw+N1gGkScHOJcWA1EKej0BqBLnukLF3H3mNm0Lbdtx/v7Xg4FcLNpOUkkEjf08mjO2tkrwqFDWILQbaLOr7ZrEyick3rbATs//eYjHOZoQUyBLB27cc2c3kLX8TmHmJz51deSY/FykljUO707LzKEsmE7isaGcWOjJKIxcOrSFt5Sq+PhfH3XffbXX0aPBNEcogKxS1iC1qdv1rYiAKayrStZCYVtG3HNnNjDW/4FKo51PgmTyTZnPPpu2g850VnvM06gs4t/IX0s7vp/8dd7D8558rfF6hUNQ8VzTQQojXyiqXUr5ZVrmi8iyLPcisw+uIz02niYsXUyMH4OXqTFpO6SOTUgsUCe9P3vI3zoV6ZgBvFN2fCdyXFc9orUCUc5zfYNBzcNMXZKfF8dC055n33vvKOCsUdRBbXBzZxd47Yzp+d8Q+w7n2WBZ7kOl7lqI3mgKdLuSmM33FcjS5pY+zSSQF3uB00XTdKPMSzwD/w3Sg/WvgEcCYkUZ6OHicuZy9u7j51Wp1+DRuQ8fut/Pt+zPtNjeFQlE1rqjFIaX8sNjrbaAfcMX4eiGEsxBihxBivxDisBDijSvVuRZ5+8BKi3E2I8/qMJZx+lFiOnUhtSANBu52cOR/gCPwKybjDJDg4YMu1ZzIVZqeB5LP7yPr0nkAWnYazozpT9prWgpFtVBZudEvvvjC1xwZ2LZt2wiNRtN169atLmA/udFly5Z5FO/Tycnpuh9++MEbakZutCSu2KZmlw8MkFJ2AjoDQ4UQPa6ivwbLn0ePkHRWUHDOjYILrhiyHTBkO0B+2e4GDQJtFuS55ROzbB6/FxbgDvwNmPOW5Dro+N/1Q3FLBm2BKXegKDQSG72S4zsXEHciCk9vF6Y/NURtACqqlT/+2ed7+6NfdOg96oOutz/6RYc//tlX43KjTzzxRKo5MnD+/PlnmjZtWtCzZ89csJ/c6IgRIzLNfW7YsOGYs7Oz8Y477sioEblRIcRBIcSBotdh4BgmbfcKkSbM2cB1Ra/yo2KuMcwSoRg0gACDBkOqE4ZUJyuN5ZI0yjSQtuR/ZJ05gtbFlTtuHEJbDx+MwAUPH2YMvJsov+ssokeGwgKO7fiB88fWEtKmF3Gnd/D3vKeUcVZUK3/8s8/3k3nrQy5eMokjXbyU7fjJvPUhVTXSlZUbLc78+fN9R44cmQrYVW60OD/88INP375906tLbtQWH/Stxd4XAolSSpsS3xVFIu4GWgGfSSm3l/HMeGA8QEBAAFFRUbY0Xa/JysrirXUHLRKhlzGtnE1+ZoljlrRK0prplMWeVV+Tl3wBB19Pmr0+jh2yGUNuHIyQl1fdnmdNLhN9fhaHt3xNdno8oR1G0LRlH7Zurbm0VFlZWdfE5wnX1lzL4ttftgUV6EvIjeoNmm9/2RZ0x5DOVQr3NmOL3Ghx/vzzT58lS5achKrJjU6bNi0RypcbbdKkieUHefHixb6TJ09OBKzkRm+66aacapUbLfHMeSllvhCiH3CnEGK+lDLtShWllAagc5Ge9O9CiPZSykMlnpkDzAFTqHd9Dpe1laioKFIr8EUZ3MEh1ojLRSwr4fyMVE5u+Zq87BR0jQJoPPFxhMYHh3PCyjiDyZgLAzjoXHB29yckYhg+gW1pHOBZo+HI9T38uTJcS3MtC/PK2dbyymKr3KiZdevWubm4uBi7deuWB+VKhZZbvzJyo+b3MTExumPHjrmMGjUqA0Cj0WBPuVEzvwEGIUQrYC4QBiysTCdFxjwKGFqp0TVgmnhUnB7KJf2ycc7OSODgxs/Iy07B1TuIJk8/hc7XDwwahL70lyw+dT/5+RkIjZa2NzyAT2BbFaatsCt+PuXIjZZTXhkqIzdqZsGCBb6jRo2yrNztKTdqvj9//nyfoUOHpjk5OVmMtj3lRs0Yi1wao4CPpZRTgSZXqiSECChaOSOEcMGkjHe0MoNriGxJ+ZefnX7HrdERtJryFwKaoq91ZmoMhzZ+TkFeBp5+LWjfewIeqR4W/efi/mopjSRs/ouz637kbOx6DM4SiaRxgCcvPjFY+Z0VduPhu2+Mc9SVkBvVaY0P331jjcqNAhgMBv766y+fcePGWQy0PeVGzSxevNh3zJgxVu4ce8qNmtELIe4DxgEjisoq9N0U0QT4vsgPrQF+kVL+VZnBNTS2pPzLvDPzKdAU0CjAVBZ7rjH5BQ5QQl/D4ARZ545z5N/vMRoK8A2MIPyG+xGOOtPJDAM4pUr0rqDLkRjz8zm/ciGZpw7h07EHfo8P4lLTXJq6ePHbMNvU5xSKq8XsZ/72l21BFy9lO/r5uBU8fPeNcVX1P1dWbhRgxYoVHoGBgQVm14QZe8mNgsk/Hh8f73jLLbdYhQDbVW4UQAgRAUwEtkkpfxJChAH3SCnfrUxHttDQ5Uan7H2BiwWlv69b/m1HSQOt37CPuD8WIqWBgOZdaX3daNBqyfUT6ItlTTFqJdmaSyR8O5f8lAQCb74Nr1HdoLEBZ62O/3S5lRHBHew9tTK5lvyy9X2uSm609rhquVEAKWU08Eyx6zNAtRvna4GSxjkl25W4TG/QSjBcNtAZ2/7l4h+LQUpatb+ZRq0GIx005PlYG2cwbQZqWoLw0tBs4v20vrGrVch4bRlnhUJRdWxxcSiukj+PHmHm1s1cyMxE6wB4huHimU+QRxqZ+Y6k5HoAAq1XAcNXH+K5ZSv4MS2N6UX1/e4ciPa+3mRsKdujlHHqEK6twnFspCP43ccJcvVm3bDJNTU9hUJhZ64mklBhA+ZAlAuZJpeToRAMl5zJzXDmbJqfxTgD3LZ7N2///CufFBlnAczSaLjfJYCCeDeMDtZuKGk0krBxKbFL53HxbBQALg6OTI1UCrAKRUOiQgMthNAKIZSazlUwc+vmUoEo2kyB4wFntJs9cNjpjkgy/QHzzA9reUZfyAeY/qRZAEwxGnl++UowaCjwksiiM3eG/Dxil84jZXcUPjf2wOeeG9AKUau+ZoVCYR9sSXnVVQgh5JV2ExVWmFfOZrRZ4JTK5aCSfIH2pDP5WTlMSU3nD8AF06HzYUV1ml66BIDBXZAPcDaFc7/PIz81iYDbRuL7YOda3whUKBT2wxYf9F7gTyHEr1invFpit1HVc7ak/IuTk578/Mu+Y8c0SkX8GfPyufDp95wEvIG/gF7F7l/wuSwVYHCHAv9CDLKAwCcexyW8NU1dHNRGoELRgLHFB+0LXMSU4mpE0evWCmtc4/x6bgnBzRIprg0lrBVFKczN4uziL8iJPYmrszMrHRysjHOOTscHw01r6dxTp5BS4hjYmGb/9xIu4eE09fBk3bDJyjgrGjSVlRvNz88Xo0aNCg0PD49o0aJF5MsvvxxobstecqMAEydObNaqVavIFi1aRBZv2+5yo1LKh8t4PXKletcCxpylGJP6YUxoY/o3ZylgOk7XKCADB4fLVrl4xJ8+M40zv3xGbuI5dF6++D4/lZ/uvZvzPiZVuvM+Pky/ZzR/XteFi7//ScKnn5O9bz8AQqvFxcGBaT1V2LaibrF80XbfMf3f7TC04/91HdP/3Q7LF22vcbnRb7/91qegoEBz/Pjx6P379x+ZP39+gNmQ2ktudPXq1W47duxwP3r06OHjx48f3rdvn9vff//tUR1yo7akvHIGHgUiMWVUAeBaN9LGnKWQ8QpQFFpvvAAZr2AE/Bx9uViQSlhIAidPN0VKjUmhLlVScDGZs0u+Qp95CSf/QJo/8xAGLz+W+vuz9PquRa1LjHm5JM35htyjx/C8qQ9uHU0rZR9nZ17r25/b27arhVkrFGWzfNF236/e/zukoKBQA5Cakun41ft/hwAMv6f7VUcThoSE6M2qdSXlRjds2HAMTHKjffv2bQPECSHIycnR6PV6srOzhU6nk97e3obicqOARW707rvvzije35XkRmfMmHEBTHKjL774YrDRaEQIQX5+vsjLyxNSSlFYWCiaNm2qL0tu9L///W/g1KlTbQ7+scXF8QMQCAwBNgDNgEqFKzZIsj4C8ojVO/B3jhu/Zbvzd46W2NR3Gd18FI4aRxoFZNCqxQUcdIUY3CXp+ec4/eun6DMv4RoUTOjLj6KNdEHrk48pgaAEJPqLyVz45BNyT5zA79478Rt1O0He3swaMozdE55UxllR51jw5fogs3E2U1BQqFnw5forZl+yFVvkRh966KFLrq6uxkaNGnUKCwvrOGnSpITGjRsbqiI3OmzYsAwoX2504MCB2b169cps0qRJp6ZNm3bs379/xnXXXZdXXG4UwF5yo62klKOFELdLKb8XQiwE/qlMJw0SYzzfxbdgzvnOpOld8NLlMigwmlzf87RwXM8jYeP49dwSCDC5OxL3xBP17WoMOXqa9AjCd8p9FOqKf1YC87lofeJFDBk5NJ32MB7dgvF2NLL51sdrZZoKhS2kXsws0/CUV15ZbJUb3bBhg6tGo5EJCQkHUlJStL169Wp7yy23ZNhTbvTQoUNOx48fdz5//vwBgL59+4avWLHCfdiwYVk1ITdq/q2TJoRoD3gBoZXqpQHyV0pnZp29gTS9KyBI17uy+FxXXt1/O09uT+d46glGNzclojq3MYa1k/9Bn6PHo1cknhMfJSfZ/3KqqzRHkIKC+HgAXCPa0fyV6TgFtQUgvSC3tqapUNiEr59HmbKi5ZVXhsrIjf7www9+Q4YMSXdycpJBQUGF3bp1y9q6daubPeVGFy1a5N2tW7dsLy8vo5eXl3HgwIHpW7ZscYOakRudI4TwAV4FlgLRwPuV6aQhMiu2M3pZ8g8Q0yo4Xe/Ke4dOM37TSvb/HsPGl9Zi1BvwGnw9jR6/l4J0T6tUV1Jv5OJvvxP3/ofknY0BQOPsbNHnaOLiVaNzUygqy9iJ/eMcHR2s5UYdHYxjJ/avUbnR4ODggvXr13sajUYyMjI0e/bscevQoUOePeVGg4ODC7Zs2eKh1+vJz88XW7Zs8YiIiMiDGpAblVJ+U/R2A7Yli70mSMgr/7jM8I37eXbBWn5JSefZojLfO2/C/74B6OOdodh5aENODknf/0DeseN49uuLU3Dzyw1pJc5anQrhVtR5zBuBC75cH5R6MdPR18+jYOzE/nFV2SCEysuNvvDCC0n33ntvaHh4eKSUkjFjxqR07949F+wnN/rwww9fWr9+vWebNm0ihRD0798/fcyYMelQM3KjjYF3gKZSymFF8qM3SinnVqYjW6hPcqMDVszmQm6pz5fhG/fzny+W8t+CQt4qKntXqyFp0h0sv6kTBefcsPiak5JJ/Hou+tRU/EePwqNHsaTnQuIVYOTNXkPr7Vnn+i7BWRnq+1yV3GjtUZHcqC0uju8wbQqa/TPHgSlXqiSEaC6EWC+EOCKEOCyEaFAya1MjB+BQxibDlB/WMq3IOGuBb4EXDUaeXbDW9ID28i/EnOhoDDk5BD45Ae9e1+PkWABInBwLaNsynv8N6lNvjbNCoag6tpzi8JdS/iKEeBlASlkoRMm4uDIpBJ6TUu4RQngAu4UQq4v0pes9vq7ZNPdKJi7TmwKDaWe2MF0wLTWdRYATsAi4vej5JhdNq22tVwH5Z3Jx8PLBs+9NuHe9DkcvN1qEXqBRgNWRTH49t4Re/j1QKBTXJrYY6GwhhB9FcctCiB5A6b/tSyCljAfii95nCiGOAEGYNhnrLVtS/uXXc0u4WJCKvxv4u+UAkHhRx+ZXo1gEeAB/Av2L1Yv380IWGkhZ8BcZmw7S9MVn0Xn54uDtWqZxhtIC/wqF4trCFgP9LKbTGy2FEFuAAOCuynQihAgFugDby7g3HhgPEBAQQFRUVGWatgsntWfY5bCPbJGNm3Tj+sLOtDKEkeS6mn+Ip0Bae4byM/LZ99IycqOTcHFyYoXBQK9iUqM5Oh1v39mH82/9SM7B0/je0QvXdjqENhuQhPkZLqtQFcPN6FYn/j+ulqysrHo9/spwLc1VUXPYcopjjxCiL9AG0+7WMSmlzYofQgh3TCqaU6SUpZaJUso5wBwwbRLW9kbLlpR/2XZmJwVG0/HNbJHNNuedaL3OE3UpAWMJt31OSg7rJq8k7dQlHHy88J04gV/PnSdk+QqaXrrEBR8fpt/Uk5//3EphSjqBk0bi1b+zpb6fkwsPtL7PlEzWePnIqKPGkQda3levXRz1feOsMlxLc1XUHOUaaCHETeXculEIgZRy45UaF0LoMBnnBfVFnvTXc0usDCVAgbGAdZfOUjKxa2ZcBmufWUlWXCaeIV50eu1OLmQHsLRx48u6GkJy8Y9FGHPzaf7GQ7i0DbbUd9RoeKnDMHr5d7D0fbEgFT9HX0Y3H1WvjbNCoag6Fa2gp5VRJoFOmPQ4KoxZFKY4yrnAESnlR1c9whqmfL+vtXG+dDKVdVP+ITclB9+2/gyYNRhnH3BJvsCZ2ED0BRqMBbk4NtHSePwQDBl90PlfDjhpWiKpay//HsogKxQlOHnypG7s2LFhycnJOo1Gw4MPPpj86quvJiUmJmpHjhzZIi4uzikoKCj/zz//PB0QEGDIy8sT999/f8iBAwdchRB8+OGH52699dZMMMmNPvroo6F5eXmaAQMGpM+bN++cRmP9F3FmZqZmxIgRLWJiYpy0Wi2DBw9O+/zzz+PAFGxy1113hR08eNDV29u78Ndffz3dpk2bAoAnnngiaM2aNd4AL7zwwoXHH3/8EpjkRl966aVmer1e06FDh+xFixad1ekqlACxotxjdlLKEcVfwHuADtPG3x02tN0LeAAYIITYV/S6xeaR1RJ+jldWSEw+mMjqJ5aTm5JD465NGPjZMLIc/difEERMgRdOjdPJXfUDCZ98AiITjaPOYpydtToedmujtJwVDY5lX67yvSdofIdB2tFd7wka32HZl6tqXG501qxZ/gDHjx+PXrdu3fEXX3yxmcFgOnRmL7nRn3/+2Wv//v2u0dHRh3fv3n1k9uzZgampqZrqkBu94jloIcTNQogo4D/AR1LKHlLKZVeqJ6XcLKUUUsqOUsrORa+/KzO42sCsn1GSlGxXmi66yNDBP7Hp8b8oyCygYxs/Bnw0mAy8iUn3pcDggCEzl1OvLyRxxQGa3RRCmyY5+Dk5IzCtmv/T5Va6OzWu2UkpFHZm2ZerfL989ruQ1PhLjkhIjb/k+OWz34VU1UiHhIToe/funQOl5UYnTJhwEUxyoytWrPABiI6OdhkwYEAGQFBQUKGnp6dh48aNrsXlRjUajUVutGR/V5IbfeSRRy6CKXpw69atHkajkcOHDzv37t07S6fT4enpaYyIiMhZsmSJV1lyo3/88Yd3ZeZfroEWQgwXQmwFngf+T0rZX0q5ujKN10d6+fdAlHBnpGS70n5ZHO0/XcbDGTnkAA8BW86k0Xv9WeIyvTFKDfnnkoh5cQ55x8/TZPKduN8zjJc73ceW4S9wZNRratWsaLD8+J/FQQV5emu50Ty95sf/LK5RudFOnTrlLFu2zFuv13P06FHHQ4cOucbExDjaU260S5cuuWvWrPHKzMzUxMfHO2zdutXz3LlzjvaWG10GnMeU7urFktJ8UsrbKtNRfWJAo76sTYqyXMdletP36695sNCAxHTucCagKTBw95e7+OQz04nnpLkrMObraf7mQ7iENyffIJVfWXFNkJpwqWy50XLKK4utcqOTJ09OOXLkiEuHDh0igoKC8q+77rosBweH8qRCy+2vMnKjo0aNyti+fbtrt27d2vr6+uqL+pQajYaqyo1WZKD7V3CvwWLMWUo/3eeccvInJj8Ao5TE/7KFl7NNKoHvAC9xecvQNyEbB2M+hRpnmjwzEmmUFn+zj2OVlRYVinqBb6BPQWp8aWPsG+hjV7nRkJAQfXG5UZ1OR3EVui5durRt165dnr+/v6E8udH27dtHAAwdOjTNrGhXkdxoy5Yt9cXlRgHee++9hPfeey8BYMSIEWHh4eH5cFluFGDJkiWeJ0+etGSlsoWKNgk3VPSqTCf1BXMaq0MFuXT0iGO4717Sv15OysJ1COBL4GUuG+cC4AEXB5I+WIgwFuLg62kxzjphYHx429qZiEJRw9z/6l1xjs46a7lRZ53x/lfvqlG50czMTE1GRoYG4Pfff/fUarWya9eudpUbLSwsJCEhQQuwfft2l6NHj7qOGjUqHWpAbvSaoiiNVa50x1Ao+XJ6DBuWpCIcNAwfdD3j1u2FfJMbKxkYKWBLbiEdIn1o6n2JC9k+FBi0uDvkM7TJIRo57SIm058Qj2G1Oi2Fwt6MmDg4FUy+6NSES46+gT4F9796V5y5/GqprNzohQsXHIYMGRKu0WhkYGCgfuHChWfMbdlLbrSgoED06tWrLYC7u7vh+++/P20+Smd3udGapLblRo0JbYjVa9mS6cjHk8+ya006Ti4aHvy4I+db9aL9X3E888N6UlLSGaERJGrgvv92xvvmgVwsSMXLwYXWLidp6mT5qwitcKZrwCtWRvpaiTq7VuYJ9X+uSm609qhIbvSaXUHHZK7gUOr/yClMRCc8EEJQYHQnJ7OQmRNPcXh7Fm5eWl7+uiXh12k4n3eUI/c15YW7RrL8viXI3DxmfBlGq04aRrc0JZhZHjOcnMIkq34MMo9Dqf9Tq2iFQlFpKjTQQogHgcmYdDgAjgCfSCnn23tg9iQmcwW7k9/CIE0bf3qZARIyLup559FTnD6Ug08jHf/3XSuCw10ACHJKpaljKhqNoNOnTfHyc8A30BFXh8tnmnMKE8vsr7xyhUKhqIiKtDjGYRLmfxbYg2lv7DpgZpEWR7010odS/2cxzmZSLhTw1kMnuHA6n8bBTrzyXSsaBzsBoM838tUrsXgH6Lj/hSDCIl0t9dr7TrK8d3VoTE5hQqn+ihtxhaIeYTQajUKj0dQdP2gDw2g0CqDcY4MVRRI+CYyUUq6XUqZLKdOklOuAO4vu1TtiMlcUuSGsjajh6wTe7neIC6fzae8Acx8JoHGwE+fzfPjzRAsm33uBjb+nkq9ztToL2cLjLivXRXvfSWiF9SkarXC2MuIKRT3iUHJysleREVFUM0ajUSQnJ3sBh8p7piIXh6eU8mzJQinlWSFEmTHsdZmSbg0zWZ9d4IVZCaRgEg9ZVgge78axytmT9UFurJu2grxLefR5ewCBNwcTl3+OcHcn2vtOKuVXNl+bfduuDo3LfE6hqA8UFhY+lpCQ8E1CQkJ7bEuPp6gcRuBQYWHhY+U9UJGBzr3Ke3WSstwah7Zl8tHHCWQBw4DFgCtAniTio7Osyj+F1lHL4K9uxa+tPwbgSE5TpkaW790J8RimDLKiQdC1a9ckoMFGDNcHKjLQ7YQQB8ooF0ALO43HbpTcqNuxOo3Zk8+gl3Afpsy4xcOgQpJzuPGDQfi28cfV/7LPOddou1SgQqFQVIUKDXSNjcKOGHOWQtZHuAgDuUWpqqJ+u8gXL8cgjfCwqwPf5BSiAfKAx4GhwNDGbjTrFVyqPS+tzclkFAqFokpUZKBdpJRHAYQQTlLKfPONosSxMfYeXFUxh25DHoEaJ84YdPw1N4n5/zVFn3Z4pDNdm3uif28Ll/IM3IEpaWKkA2yfEoIWA4ZieQl0wsC9TZT4kUKhqBkqcvwvLPZ+W4l7n1+pYSHEPCFEkhCi3B1Ku1MUug0Qb9Dy04cXLMb5+qnd6TS+K/8Oa82MBzvRVSM4CCzw0dLv/RCy73Kmk/s5XDQFgMRLm8cjvnH0DqqXB1gUCkU9pKIVtCjnfVnXZfEd8D+g9s5LG+MBMBgkn7x6njU/XUSjhR6v3ESLYa0ByE7MYub3B3AKcKXfzIHoep7jHCCQNHNOpbXLRdrr8gnWOYDnW7U2FYVCce1RkYGW5bwv67p0ZSk3CiFCr2ZQ1YcXBQWXGPd0AmuWZqFzFEz9NIxLXULILToa7tbYna5P30DzfiH4Blz+75DAnW4OQDpomoD7s2hc1Ya2QqGoOSoy0M2EEJ9gWi2b31N0XW1ZEoQQ44HxAAEBAURFRVVLuy0bLcDLKZW7Ho1n9YYcPDw0TPuyBW27e3I67RxfzzhH6zsj8Y8MIPzOdmgx0M71srqgg8GNjUf/U6LV6hlbVlZWtc2zLnOtzBOurbkqag5bs3rvKnGv5PVVI6WcA8wBk5pddSiCGXOWkhqznqH3xrFtVx4Bflr+XtgU/3Y6NpzP5+sJOzh9KIeg9r4Q6Y+LpoB2rvE0c74EmKL/rmsynZDwqo+lLOq78pmtXCvzhGtrroqao1wDLaX8viYHUp1cOPlfht1zjkNHCwgOcuCfRUGEt3Rk174spj98gfQMI7/Na8LtQ9OATcTqHTikdyJXClw03nTwn6aCTRQKRa1TkVjS0ooq1tWchCdPnmTwiN2cidXTrrUjK39uSrOmOvYcyKPvyPM09teyeVlzOkY4WeoE6woJdmqk/MwKhaJOUZGL40bgHPATpuPBlRJMEUL8BPQD/IUQ54HXpZRzr3KcNnHgwAEGDx5MYqKebp2d+OvHIPz9TOeYO0U6MeVxbyaP96aRf/FpO4PnW8owKxSKOkdFBjoQGIQpEnoMsBz4SUp52JaGpZT3VX14trNlyxaGDx9Oeno6N/fryG9f56DVSCa9nMT0yb40DXTg7en+1pU0TdWqWaFQ1FkqShprkFKulFI+CPQATgJRQoina2x0NrJixQoGDRpEeno6o0aNYvnKHaTn30rfO87z5ffprN+SU7qSpimaRlHKOCsUijrLlTKqOAHDMa2iQ4FPgCX2H5bt/PTTT4wbN47CwkIeffRRvvzyS3bv3s0dd/xMdraGP+e3ZvjAkrWcwf3Z2hiuQqFQ2Ey5K2ghxPfAVkxZVN6QUnaTUv5HSlmlNOrVyeeff87YsWMpLCxk2rRpfP3112zYsIG+ffvi4uLCtm07GXH/cfD8wOTOQJj+VT5nhUJRD6hoBf0AkA2EA88IYdkjFICUUtaaaL+UkrfeeovXXnsNgHfffZcXX3wRgK5duzJ27Fjee+89/P1NPmeN622gDLJCoahnVOSD1kgpPYpensVeHrVpnI1GI1OnTuW1115Do9EwZ84cJk2axCuvvEJeXh7e3t7MnTvXYpwVCoWivlKhD7quodfrefTRR/nhhx/Q6XQsXLiQG264gV69enHw4EF69+7N0KFDa3uYCoVCUS3UGwOdm5vLPffcw7Jly3Bzc+P333/H3d2dG264gdzcXJYvX66Ms0KhaFDUi0SQ6enpDB06lGXLluHj48PatWvJysqiX79+uLu78++//yrjrFAoGhx13kAnJSXRv39/Nm7cSNOmTdm0aRPdu3cnIiKC4cOHs337dtq1axDZuRQKhcKKOm2gY2Ji6N27N3v37qVVq1b8888/rFu3Diklbdq0YcmSJfj5+dX2MBUKhcIu1FkDHR0dTa9evThx4gSdO3dm4cKFjBkzhqlTp7J3797aHp5CoVDYnTq5Sbhjxw5uueUWLl68SO/evZk+fTrDhw+noKCAFStWcN1119X2EBUKhcLu1DkDvXbtWm6//Xays7MZPnw4I0eO5Pbbbyc0NJRly5bRpk2b2h6iQqFQ1Ah1ysWRlZXFLbfcQnZ2NmPHjuX3338nKCiIAQMGsH37dmWcFQrFNUWdWkFfuHABgAkTJjB48GB0Oh1Dhw5lyJAhFAs1VygUimsCu66ghRBDhRDHhBAnhRAv2VLnmWeeYfPmzdx3332cP3/e3I49h6lQKBR1ErutoIUQWuAzTKL/54GdQoilUsro8up4e3uzYMECjEYjK1asoFmzZvYankKhUNR57LmCvgE4KaU8LaUsAH4Gbq+oQlpaGgEBAezYsYMBAwbYcWgKhUJR97GnDzoIU05DM+eB7hVVcHV1Zdu2bXh7e9txWAqFQlE/sKeBLstxLEs9JMR4YDxAQEAA+/bts+OQ6gZZWVlERUXV9jDszrUyT7i25qqoOexpoM8DzYtdNwMulHxISjkHmAPQpk0b2a9fPzsOqW4QFRWFmmfD4lqaq6LmsKcPeifQWggRJoRwBO4FltqxP4VCoWhQ2G0FLaUsFEJMAv4BtMA8KeVhe/WnUCgUDQ27BqpIKf8G/rZnHwqFQtFQqVOh3gqFQqG4jDLQCoVCUUcRUpY6+VZrCCEygWO1PY4awB9Iqe1B1ADXyjyh/s81REoZUNuDUFhTp8SSgGNSyutrexD2RgixS82zYXEtzVVRcygXh0KhUNRRlIFWKBSKOkpdM9BzansANYSaZ8PjWpqrooaoU5uECoVCobhMXVtBKxQKhaIIZaAVCoWijlInDPTVpMaqjwghmgsh1gshjgghDgshJtf2mOyJEEIrhNgrhPirtsdiL4QQ3kKIxUKIo0Wf6421PSZFw6HWfdBFqbGOUyw1FnBfRamx6itCiCZAEynlHiGEB7AbuKMhzhVACPEscD3gKaW8tbbHYw+EEN8Dm6SU3xSpNrpKKdNqeViKBkJdWEFXOjVWfUVKGS+l3FP0PhM4ginzTINDCNEMGA58U9tjsRdCCE/gJmAugJSyQBlnRXVSFwx0WamxGqTRKo4QIhToAmyv5aHYi4+BFwBjLY/DnrQAkoFvi1w53wgh3Gp7UIqGQ10w0DalxmpICCHcgd+AKVLKjNoeT3UjhLgVSJJS7q7tsdgZB+A64AspZRcgG2iweyiKmqcuGGibUmM1FIQQOkzGeYGUckltj8dO9AJuE0KcxeSyGiCE+LF2h2QXzgPnpZTmv4IWYzLYCkW1UBcM9DWTGksIITD5K49IKT+q7fHYCynly1LKZlLKUEyf5zop5f21PKxqR0qZAJwTQrQpKroZaJAbvoraodbV7K6x1Fi9gAeAg0KIfUVl04syzyjqJ08DC4oWF6eBh2t5PIoGRK0fs1MoFApF2dQFF4dCoVAoykAZaIVCoaijKAOtUCgUdRRloBUKhaKOogy0QqFQ1FGUga6nCCEChRA/CyFOCSGihRB/CyHChRChQggphPhPsWf9hRB6IcT/iq5nCCGeL3r/nRDijBBinxBiT3E1NiHE00Uqg4eFEO+XM47WQoi/isaxu0it7yZhIkUI4VP0XJOicfUuVjdZCOFXfDwVzFcIIdYV6V8ghMgqdu8WIcQJIUSwEGKSEEIddVM0CJSBrocUBbz8DkRJKVtKKSOA6UDjokdOA8XV40YDFZ0tnyal7IwpTPmroj76YxKt6iiljAQ+KGMczsByYE7ROLpiOhfcQprOb24HzAa/J7C36F+KgjtSpJQXbZz2LcD+kqHxQoibgU+BoVLKWGAe8IyNbSoUdRploOsn/QG9lPJLc4GUcp+UclPRZS5wRAhxfdH1PcAvNrS7EWhV9P4J4F0pZX5R+0llPD8W2CaltER+SikPSSm/K7rcQpFBLvr3I6wN9lYbxlS8rz+LFwgh+gBfA8OllKeK+s8BzgohbqhE2wpFnUQZ6PpJe0xa0hXxM3BvkeynAdv0TUYAB4vehwN9hBDbhRAbhBDdyng+EthTQXtbuWygbwD+4LLuSk9MBtxWemE9ZydMBvsOKeXREs/uAvpUom2Fok6iDHTDZSWmJAj3AYuu8OzMotDz8cCjRWUOgA/QA5gG/FLkWikXIcTvQohDQgizCNQOoEuRBKdOSpkFnBZCtKLyK2jfIg1tM/qi+o+W8WwS0LQSbSsUdRJloOsnh4GuFT1QlPxgN/AcJvW8ipgmpewspRwkpTxUVHYeWCJN7MCk6+xfxjgs6m1SypHAQ4Bv0XUOcBJ4hMsr7X8x+ZMbAceuMK7iFAohin9fjcDdQDchxPQSzzpjcvMoFPUaZaDrJ+sAJyHE4+YCIUQ3IUTfEs99CLxYiY244vwBDChqOxxwBFJKPLMQ6CWEuK1YmWuJZ7YAU4BtRdfbgMnAv7JyQjDHMAnkWyj6BXArMFYIUXwlHQ4cQqGo5ygDXQ8pMmwjgUFFx9sOAzMo4WeWUh6WUn5/ld3MA1oIIQ5h8mc/WNKgSilzMRnIiUKI00KIbcArwFvFHtuCybCaDfQeTJrfJd0brwghzptfZYxnOdCvZKGUMhUYWlTfnCqtF7DG5pkqFHUUpWanqBcUJdydL6UcdIXnugDPSikfqJmRKRT2Q62gFfUCKWU88LU5UKUC/IFXa2BICoXdUStohUKhqKOoFbRCoVDUUZSBVigUijqKMtAKhUJRR1EGWqFQKOooykArFApFHeX/Acagiu/TFhhkAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "decades = sorted(list(set(data.index.get_level_values(2))))\n",
    "colors = cm.get_cmap('viridis_r', len(decades))\n",
    "fig, ax = plt.subplots(figsize=(4,4))\n",
    "h = []\n",
    "for k,decade in enumerate(decades):\n",
    "  decdata = data.xs(decade, level=2)\n",
    "  h.append(ax.scatter(x=decdata['world'], y=decdata[region], color=colors(k)))\n",
    "# Set legend outside the plot\n",
    "ax.legend(h, decades, loc='center left', bbox_to_anchor=(1, 0.5), ncol=1)\n",
    "# Add decadal means\n",
    "decadal_means = data.groupby('decade').mean()\n",
    "ax.plot(decadal_means['world'], decadal_means[region], 'or')\n",
    "# Add linear fit\n",
    "lm = linear_model.LinearRegression()\n",
    "X = data['world'].values.reshape(-1, 1)\n",
    "y = data[region]\n",
    "model = lm.fit(X, y)\n",
    "ax.set_xlim(left=0)\n",
    "#ax.set_xlim(0,8)\n",
    "#ax.set_ylim(0,16)\n",
    "ax.text(0.05, 0.92, '$\\\\beta:$ %.2f %s/K\\n$R^2:$ %4.2f' % (model.coef_, units[var], lm.score(X, y)), color='k',\n",
    "        verticalalignment='top', horizontalalignment='left', transform=ax.transAxes, \n",
    "        fontsize=15, fontweight = 'bold', bbox={'facecolor': 'w', 'alpha': 0.5, 'pad': 5})\n",
    "X[0] = 0\n",
    "ax.plot(X, model.predict(X), color = 'k', linewidth = 2)\n",
    "# Grid and labels\n",
    "ax.set_axisbelow(True)\n",
    "ax.grid()\n",
    "ax.axline((1, 1), slope=1, ls=\"--\", c=\"k\")\n",
    "plt.xlabel(f'{project_gsat} GWL (K)')\n",
    "plt.ylabel(f'{region} {longname[var]} ({units[var]})')\n",
    "plt.savefig(f'GWL-plot-{project}_{scenario}_{region}_{var}_{mask}.png', dpi=150, bbox_inches = 'tight', facecolor = 'w')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}