{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Estimate climate driver influence on rainfall \n", "\n", "* **Compatibility:** Notebook currently compatible the `NCI` environment\n", "* **Products used:** [BOM Rainfall Grids](https://explorer.nci.dea.ga.gov.au/products/bom_rainfall_grids)\n", "\n", "## Background\n", "\n", "There are three main climate drivers for the Australian continent: the El Niño-Southern Oscillation (ENSO), the Southern Annular Mode (SAM), and the Indian Ocean Dipole (IOD). ENSO is the climate driver associated with Pacific Ocean El Niño and La Niña events (ENSO phases). These events together affect climate and rainfall patterns in Australia. Different regions of Australia are affected quite differently. We are interested in the behaviour of waterbodies throughout the continent, but we need a rainfall baseline to compare these to.\n", "\n", "## Description\n", "\n", "This notebook estimates the influence of climate drivers on Australian rainfall.\n", "\n", "***" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Getting started\n", "\n", "Run the first cell, which loads all modules needed for this notebook." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Load modules" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import warnings\n", "\n", "%matplotlib inline\n", "import datacube\n", "import numpy as np\n", "import pandas as pd\n", "import xarray\n", "from matplotlib import pyplot as plt\n", "from odc.ui import with_ui_cbk\n", "from tqdm.notebook import tqdm\n", "import statsmodels.tsa.seasonal\n", "import matplotlib.patches\n", "import mpl_toolkits.basemap\n", "\n", "import sys\n", "sys.path.insert(1, '../Tools/')\n", "from dea_tools.plotting import display_map" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Configuration\n", "\n", "Choose a subset of Australia to examine:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "southwest = (-44.367634, 112.417252)\n", "northeast = (-10.116211, 154.898971)\n", "\n", "# Convert the coordinates into latitude and longitude ranges.\n", "ylim, xlim = zip(southwest, northeast)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Show a map of the the selected region.\n", "display_map(xlim, ylim)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Connect to the datacube" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "dc = datacube.Datacube(app=\"estimate-climate-driver-influence-on-rainfall\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Load gridded rainfall" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "8157cc42d49842e3b315a30bb82f7455", "version_major": 2, "version_minor": 0 }, "text/plain": [ "VBox(children=(HBox(children=(Label(value=''), Label(value='')), layout=Layout(justify_content='space-between'…" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "with warnings.catch_warnings():\n", " warnings.simplefilter(\"ignore\", DeprecationWarning)\n", " # Load all data between 1988 and the time of writing,\n", " # and downsample to 0.5 degree resolution.\n", " rainfall = dc.load(\n", " \"bom_rainfall_grids\",\n", " time=(\"1988-01\", \"2020-10-19\"),\n", " latitude=ylim,\n", " longitude=xlim,\n", " progress_cbk=with_ui_cbk(),\n", " resolution=(0.5, 0.5),\n", " resampling=\"average\",\n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Load Climate Indices\n", "\n", "The Southern Oscillation Index (SOI) tracks ENSO based on pressure differences between Tahiti and Darwin. SAM is tracked by the Antarctic Oscillation Index (AAO), a measure of pressure difference across the Antarctic. Finally, the Dipole Mode Index (DMI) tracks the IOD.\n", "\n", "The United States National Oceanic and Atmospheric Administration (NOAA) has an easily-accessed record of the SOI, AAO, and DMI, which we load here." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "--2020-10-22 14:45:58-- https://stateoftheocean.osmc.noaa.gov/atm/data/soi.nc\n", "Resolving stateoftheocean.osmc.noaa.gov... 161.55.85.40\n", "Connecting to stateoftheocean.osmc.noaa.gov|161.55.85.40|:443... connected.\n", "HTTP request sent, awaiting response... 200 OK\n", "Length: 28604 (28K) [application/x-netcdf]\n", "Saving to: `soi.nc.4'\n", "\n", " 0K .......... .......... ....... 100% 190K=0.1s\n", "\n", "2020-10-22 14:46:02 (190 KB/s) - `soi.nc.4' saved [28604/28604]\n", "\n", "--2020-10-22 14:46:02-- https://stateoftheocean.osmc.noaa.gov/sur/data/dmi.nc\n", "Resolving stateoftheocean.osmc.noaa.gov... 161.55.85.40\n", "Connecting to stateoftheocean.osmc.noaa.gov|161.55.85.40|:443... connected.\n", "HTTP request sent, awaiting response... 200 OK\n", "Length: 33048 (32K) [application/x-netcdf]\n", "Saving to: `dmi.nc.2'\n", "\n", " 0K .......... .......... .......... .. 100% 218K=0.1s\n", "\n", "2020-10-22 14:46:03 (218 KB/s) - `dmi.nc.2' saved [33048/33048]\n", "\n", "--2020-10-22 14:46:03-- https://stateoftheocean.osmc.noaa.gov/atm/data/sam.nc\n", "Resolving stateoftheocean.osmc.noaa.gov... 161.55.85.40\n", "Connecting to stateoftheocean.osmc.noaa.gov|161.55.85.40|:443... connected.\n", "HTTP request sent, awaiting response... 200 OK\n", "Length: 16740 (16K) [application/x-netcdf]\n", "Saving to: `sam.nc.1'\n", "\n", " 0K .......... ...... 100% 112K=0.1s\n", "\n", "2020-10-22 14:46:03 (112 KB/s) - `sam.nc.1' saved [16740/16740]\n", "\n" ] } ], "source": [ "%%bash\n", "# Download the data.\n", "wget https://stateoftheocean.osmc.noaa.gov/atm/data/soi.nc\n", "wget https://stateoftheocean.osmc.noaa.gov/sur/data/dmi.nc\n", "wget https://stateoftheocean.osmc.noaa.gov/atm/data/sam.nc" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# Read the data into a dictionary of pandas DataFrames.\n", "# {'driver name': pd.DataFrame(...)}\n", "driver_indices = {\"enso\": \"soi\", \"iod\": \"dmi\", \"sam\": \"sam\"}\n", "indices = {}\n", "for driver, index in driver_indices.items():\n", " index_ = xarray.open_dataset(f\"{index}.nc\")[index.upper()]\n", " index = (\n", " pd.DataFrame({index.upper(): index_.to_pandas()})\n", " # Resample to daily.\n", " .resample(\"1D\")\n", " .mean()\n", " # Interpolate over missing data.\n", " .interpolate()\n", " # Fill out the first few rows with the first observation.\n", " .bfill()\n", " )\n", " indices[driver] = index" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Deseasonalise the rainfall data\n", "\n", "Climate drivers can be seasonal, i.e. they predominantly occur in a given season. This means that our rainfall influence estimates will be biased by the regular seasonal cycle. For this reason we will deseasonalise the rainfall data." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# Make the rainfall cube (time, height, width) into a matrix of observations (time, height * width)\n", "# We only want one spatial axis, which is simply \"which pixel\", so we can feed this into seasonal_decompose.\n", "ts = rainfall.rainfall.values.reshape((len(rainfall.time), -1))" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "# Perform the seasonal decomposition.\n", "sd = statsmodels.tsa.seasonal.seasonal_decompose(\n", " # Assume multiplicative seasonal effects because there's no such thing as negative rainfall.\n", " # But we *can* have zero rainfall, so add on a tiny quantity to account for that.\n", " ts + 1e-2, model=\"multiplicative\", period=365\n", ")" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "# Convert the deseasonalised rainfall back into an xarray.\n", "deseasonalised_rainfall = xarray.DataArray(\n", " sd.resid.reshape((len(ts),) + rainfall.rainfall.values.shape[1:]),\n", " coords=(\n", " (\"time\", rainfall.time),\n", " (\"latitude\", rainfall.latitude),\n", " (\"longitude\", rainfall.longitude),\n", " ),\n", ")" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "# Get all the times where there is missing data.\n", "null_times = deseasonalised_rainfall.isnull().any(axis=(1, 2))\n", "# ...and toss those values out.\n", "deseasonalised_rainfall = deseasonalised_rainfall[~null_times]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Estimate climate phases\n", "\n", "We first need to estimate when each climate phase was active. For ENSO, this is indicated by the SOI being below -0.8 or above +0.8 for a sustained period (see the [Bureau of Meteorology information page on the SOI](http://www.bom.gov.au/climate/enso/history/ln-2010-12/SOI-what.shtml) for more details). The [DMI uses a threshold of 0.4](http://www.bom.gov.au/climate/model-summary/#tabs=Indian-Ocean®ion=NINO34). There doesn't seem to be a standard threshold for the AAO, but we will use 0.8 (the standard deviation of the data).\n", "\n", "We will say that a climate phase is when the rolling mean climate index is above or below the appropriate threshold. This corresponds to a sustained period of being above or below the threshold. We choose a rolling window of around 4 months here for SOI and DMI, but the AAO needs a shorter rolling window: SAM tends to vary on a fortnightly time scale, compared to many months for ENSO and IOD, so we will use a rolling window of 1 month." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "# Reindex all of the climate indices to match the dates of the rainfall observations.\n", "reindexed_index = {\n", " driver: index.reindex(\n", " pd.to_datetime(deseasonalised_rainfall.time.values, \"D\")\n", " ).interpolate()[driver_indices[driver].upper()]\n", " for driver, index in indices.items()\n", "}\n", "# Choose the rolling window size for each climate driver. This roughly corresponds to the scale of the oscillation.\n", "windows = {\n", " \"enso\": 4 * 28,\n", " \"iod\": 4 * 28,\n", " \"sam\": 28,\n", "}\n", "# Compute a rolling mean of each climate index using the above windows.\n", "rolling_index = {\n", " driver: index.reindex(pd.to_datetime(deseasonalised_rainfall.time.values, \"D\"))\n", " .interpolate()\n", " .rolling(windows[driver])\n", " .mean()[driver_indices[driver].upper()]\n", " for driver, index in indices.items()\n", "}" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "# Find the positive, negative, and neutral phases for each climate driver.\n", "phases = {driver: {} for driver in rolling_index}\n", "# These are the thresholds at which we say there is a non-neutral phase:\n", "thresholds = {\"enso\": 0.8, \"sam\": 0.8, \"iod\": 0.4}\n", "for driver, r_index in rolling_index.items():\n", " phases[driver][\"+ve\"] = r_index > thresholds[driver]\n", " phases[driver][\"-ve\"] = r_index < -thresholds[driver]\n", " phases[driver][\"neutral\"] = (\n", " ~phases[driver][\"+ve\"] & ~phases[driver][\"-ve\"] & pd.notnull(r_index)\n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can plot when the positive and negative phases occurred:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAASeElEQVR4nO3df7BcZX3H8ffHBLSAGIVgMURJpxFMZ/wBMeKMraiDEh2N1XYGdIqiNoMljtqOA7bV1sFaf7YOA5pJLf5oVcYfqKGNUqUqrWIlOBoJEIhoIUKboKL80Mbot3/siS43e+/uvdl7L3l4v2Z2cs7zPOfs8+zZ87lnz56zSVUhSTrwPWC+OyBJGg8DXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQFeTknwvyU+T3NX3uCDJS5NUktdNaL8jycnd9KIkFyX5nyR3JrkhyTl9bZPkdUlu7J7j5iRvTfLAOR6mdC8Gulr23Ko6rO+xriv/IXBOksMnWe7vgcOAxwAPAZ4HfKev/nxgLXAG8GBgNfB04GOzMAZpZAa67o+uA64EXjtJ/ROBj1TVj6rql1V1fVV9AiDJcuBPgBdX1ZVVtaeqtgIvBE5N8vS5GIA0iIGu+6s3AK9N8rABdV8D/ibJmV2A93sGsKOqvt5fWFW3dMudMiu9lUZgoKtln05yR9/jj/dWVNU3gX8Dzhmw3KuADwPrgGuTbE+yuqs7Erhtkue7rauX5oWBrpY9v6oW9T3+YUL9G4FXJvnN/sKq+mlVvaWqTgSOoHdu/OPd0fztwNGTPN/RXb00Lwx03W9V1fXAJcCfT9HmJ8BbgEOBZcC/A0uTrOpvl2QpcBJw+ax1WBrCQNf93ZuAM4FFewuSvCHJE5McnORBwKuBO4BtVXUDsB74cJKTkixI8jvAJ4EvVNUX5mEMEmCgq22XTrgO/VMTG1TVd4F/oncE/qti4P30Tp/cSu+LzudU1V1d/TrgfcA/A3cBnwO+RO9KF2nexP/gQpLa4BG6JDViaKB3t0DvTHLNJPVJcn53adeWJCeMv5uSpGFGOUL/AHDqFPWrgeXdYy3w3v3vliRpuoYGelVdQe+3LyazBvhQ9XwNWJRksut0JUmzZOEY1rEEuKVvfkdXts/ddEnW0juK59BDDz3x+OOPn9ET3vGzn0+r/aIHHbRv4Z339P598CH3Wu/AtsOW3Ts/ir7nm/ayk61jOusZtOwQw17vRQ86aJ82v3odZzi+/d3Go/R5Jm2nZdjYZ7At9uu5p3i+Qa/B0G04Yv+nen33Pseo23vSbTHJeKf7PpqOmbwv9vZnxu8p4Oqrr769qhYPqhtHoGdA2cBLZ6pqA7ABYOXKlbV58+YZPeEl2ya783qwFxw34APDl7vnfurKe613YNthy355GuPoe75pLzvZOqaznkHLDjHs9X7BcUfv0+ZXr+MMx7e/23iUPs+k7bQMG/sMtsV+PfcUzzfoNRi6DUfs/1Sv797nGHV7T7otJhnvdN9H0zGT98Xe/sz4PQUk+e/J6sZxlcsOYGnf/DH0rt2VJM2hcQT6RuCM7mqXk4AfV9Xs/VmUJA009JRLko8CJwNHJtkB/BVwEEBVrQc2Ac8GtgP30LuNWpI0x4YGelWdPqS+gLPH1iNJ0ox4p6gkNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktSIkQI9yalJtiXZnuTcAfUPSXJpkm8l2ZrkzPF3VZI0laGBnmQBcCGwGlgBnJ5kxYRmZwPXVtXjgJOBdyU5eMx9lSRNYZQj9FXA9qq6qap2AxcDaya0KeDBSQIcBvwQ2DPWnkqSpjRKoC8Bbumb39GV9bsAeAxwK/Bt4NVV9cuJK0qyNsnmJJt37do1wy5LkgYZJdAzoKwmzD8L+CbwCODxwAVJDt9noaoNVbWyqlYuXrx42p2VJE1ulEDfASztmz+G3pF4vzOBS6pnO/Bd4PjxdFGSNIpRAv0qYHmSZd0XnacBGye0uRl4BkCShwPHATeNs6OSpKktHNagqvYkWQdcBiwALqqqrUnO6urXA+cBH0jybXqnaM6pqttnsd+SpAmGBjpAVW0CNk0oW983fSvwzPF2TZI0Hd4pKkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRIwV6klOTbEuyPcm5k7Q5Ock3k2xN8uXxdlOSNMzCYQ2SLAAuBE4BdgBXJdlYVdf2tVkEvAc4tapuTnLUbHVYkjTYKEfoq4DtVXVTVe0GLgbWTGjzIuCSqroZoKp2jrebkqRhRgn0JcAtffM7urJ+jwYemuRLSa5OcsagFSVZm2Rzks27du2aWY8lSQONEugZUFYT5hcCJwLPAZ4FvCHJo/dZqGpDVa2sqpWLFy+edmclSZMbeg6d3hH50r75Y4BbB7S5varuBu5OcgXwOOCGsfRSkjTUKEfoVwHLkyxLcjBwGrBxQpvPAL+bZGGSQ4AnAdeNt6uSpKkMPUKvqj1J1gGXAQuAi6pqa5Kzuvr1VXVdks8BW4BfAu+rqmtms+OSpHsb5ZQLVbUJ2DShbP2E+XcA7xhf1yRJ0+GdopLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiNGCvQkpybZlmR7knOnaPfEJL9I8gfj66IkaRRDAz3JAuBCYDWwAjg9yYpJ2r0NuGzcnZQkDTfKEfoqYHtV3VRVu4GLgTUD2r0K+CSwc4z9kySNaJRAXwLc0je/oyv7lSRLgN8H1k+1oiRrk2xOsnnXrl3T7askaQqjBHoGlNWE+XcD51TVL6ZaUVVtqKqVVbVy8eLFo/ZRkjSChSO02QEs7Zs/Brh1QpuVwMVJAI4Enp1kT1V9eiy9lCQNNUqgXwUsT7IM+D5wGvCi/gZVtWzvdJIPAP9imEvS3Boa6FW1J8k6elevLAAuqqqtSc7q6qc8by5JmhujHKFTVZuATRPKBgZ5Vb10/7slSZou7xSVpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1IiRAj3JqUm2Jdme5NwB9S9OsqV7fDXJ48bfVUnSVIYGepIFwIXAamAFcHqSFROafRd4alU9FjgP2DDujkqSpjbKEfoqYHtV3VRVu4GLgTX9Darqq1X1o272a8Ax4+2mJGmYUQJ9CXBL3/yOrmwyLwc+O6giydokm5Ns3rVr1+i9lCQNNUqgZ0BZDWyYPI1eoJ8zqL6qNlTVyqpauXjx4tF7KUkaauEIbXYAS/vmjwFundgoyWOB9wGrq+oH4+meJGlUoxyhXwUsT7IsycHAacDG/gZJHglcAvxRVd0w/m5KkoYZeoReVXuSrAMuAxYAF1XV1iRndfXrgTcCRwDvSQKwp6pWzl63JUkTjXLKharaBGyaULa+b/oVwCvG2zVJ0nR4p6gkNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktSIkQI9yalJtiXZnuTcAfVJcn5XvyXJCePvqiRpKkMDPckC4EJgNbACOD3JignNVgPLu8da4L1j7qckaYhRjtBXAdur6qaq2g1cDKyZ0GYN8KHq+RqwKMnRY+6rJGkKC0doswS4pW9+B/CkEdosAW7rb5RkLb0jeIC7kmwb8HxHAreP0K/7Msdw39DCGKCNcTiG8XnUZBWjBHoGlNUM2lBVG4ANUz5ZsrmqVo7Qr/ssx3Df0MIYoI1xOIa5Mcoplx3A0r75Y4BbZ9BGkjSLRgn0q4DlSZYlORg4Ddg4oc1G4IzuapeTgB9X1W0TVyRJmj1DT7lU1Z4k64DLgAXARVW1NclZXf16YBPwbGA7cA9w5n70acpTMgcIx3Df0MIYoI1xOIY5kKp9TnVLkg5A3ikqSY0w0CWpEbMe6EkuSrIzyTV9ZY9LcmWSbye5NMnhXflBST7YlV+X5PV9y5zelW9J8rkkR85232c4hoOTvL8r/1aSk7vyQ5L8a5Lrk2xN8ta56v+4xtBXtyHJDd1YXjiHY1ia5Ivde2Nrkld35Q9L8vkkN3b/PrRvmdd3P0mxLcmz+spP7Ma3vfvZikGX3t7nx9FXv7F/2x5IY5ivfXu6Y0hyRNf+riQX9K1nXvfte6mqWX0AvwecAFzTV3YV8NRu+mXAed30i4CLu+lDgO8Bx9L78nYncGRX93bgr2e77zMcw9nA+7vpo4Cr6f3hPAR4Wld+MPAfwOoDaQzd/JuAN3fTD9i7TeZoDEcDJ3TTDwZuoPdzFG8Hzu3KzwXe1k2vAL4FPBBYBnwHWNDVfR14Mr17KD47x9tibOPo6l8AfKR/2x4oY5jPfXsGYzgUeApwFnBB33rmdd/uf8z6EXpVXQH8cELxccAV3fTngb1HeQUcmmQh8BvAbuAn9Ha6dHUBDmcOr3Of5hhWAJd3y+0E7gBWVtU9VfXFrnw38A161+vPiXGMoat7GfC3Xd0vq2rO7pyrqtuq6hvd9J3AdfTuSF4DfLBr9kHg+d30GnoHCP9XVd+ldxXWqvR+luLwqrqyenvhh/qWOWDGAZDkMOBPgTfPVf+7fo9rDPO2b093DFV1d1X9J/CzCeuZ132733ydQ78GeF43/Yf8+qakTwB30/vJgJuBd1bVD6vq58ArgW/T29grgH+c0x7va7IxfAtYk2RhkmXAidz7piuSLAKeSxea82haY+j6DXBekm8k+XiSh89tl3uSHAs8Afgv4OHV3ffQ/XtU12yyn6RY0k1PLJ9z+zkOgPOAd9G7XHhe7M8Y7iv79ohjGGU987pvz1egvww4O8nV9D7q7O7KVwG/AB5B72PZnyX5rSQH0dvoT+jqtgCv32etc2uyMVxE7826GXg38FVgz96Fuk8fHwXOr6qb5rTH+5ruGBbSO/L4SlWdAFwJvHOuO90dlX4SeE1V/WSqpgPKaoryObW/40jyeOC3q+pTs9LBEYxhDPO+b09jDMPWM+/79ii/5TJ2VXU98EyAJI8GntNVvQj4XPdXe2eSr9D7qH9Et9x3umU+Ru/c1ryZbAxVtQd47d52Sb4K3Ni36Abgxqp699z1drAZjOEH9I4E9wbIx4GXz2GX6QLgk8CHq+qSrvh/kxxdVbd1p1N2duWT/STFDu79kXjOf6piTON4MnBiku/R25ePSvKlqjr5ABrD42H+9u1pjmGYed+35+UIPclR3b8PAP4SWN9V3Qw8PT2HAicB1wPfB1YkWdy1O4Xe+a55M9kYum+8D+2mTwH2VNW13fybgYcAr5mXTk8w3TF055svBU7uVvEM4No57G/ofRy/rqr+rq9qI/CSbvolwGf6yk9L8sDu1NFy4Ovdx+g7k5zUrfOMvmVm3RjH8d6qekRVHUvvy7ob5jDMxzIG5nHfnsEYplrXfWPfnu1vXel9BLkN+Dm9v9IvB15N7xvlG4C38us7Vg+jd9S3lV5QvK5vPWfR29Bb6IXKEbPd9xmO4VhgW9fXLwCP6sqPofex/jrgm93jFQfSGLq6R9H7InULvfOEj5zDMTylew239L2Gz6b3Ce5yep8iLgce1rfMX9C7omIbfVce0Pvkd01Xd8HesR9o4+irP5a5vcplnNtiXvbtGY7he/QuLrir249WzPe+3f/w1n9JaoR3ikpSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1Ij/B5i3zFZLwERRAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAARPUlEQVR4nO3cf4xlZX3H8ffHXdCAIgKLxV2UTbpAt4koTFETWzAWZbHt1h9tQSMWJYQWGrVJA8a2tsGkUms1BnSz0RVNVaKBlrVdJWpUatHCYmRhwYURqayQsmC1ArHLyrd/3LN4uXtn5s7snZndZ9+vZDLnPM9zznyf++MzZ86cc1NVSJL2f09b7AIkSeNhoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOg6oCS5N8lvd8srknw6ycNJHk1yU5LfGRhfXd8j3bivJvmjxalemp6BrgNSkiOAbwI7gV8HjgI+CHwmyRsGhp9UVc8ETgCuAq5I8p4FLFcaSbxTVAeSJPcC5wOnAa8FXlhVT/T1XwL8KXBcVVWSAlZV1WTfmDcA/wQsr6qHF7J+aToeoetAdQZwTX+Ydz4HPB84fpptrwOWAqfOU23SnBjoOlAdBTwwpP2Bvv6hqupx4CHgiHmoS5ozA10HqoeAY4a0H9PXP1SSg4BlwI/noS5pzgx0Hai+Arw+yeB74A+B+4C7ptl2LbALuGmeapPmxEDXgeqDwGHAx5P8SpJnJDkHeDfwFzXkaoEkRyR5E3AlcLn/ENW+ZuliFyAthqp6OMnLgcuBO4Cnd9/fXFXXDQy/tbvaZSdwK/DOqvrMghYsjcDLFiWpEZ5ykaRGzBjoSTYkeTDJ7VP0J8mHk0wm2ZLk5PGXKUmayShH6FcBZ07TvwZY1X1dAHx078uSJM3WjIFeVTcw/fW2a4FPVc+3gcOTDLu+V5I0j8Zxlctyetft7ra9a9vjLrwkF9A7iufQQw895cQTT5zbT/zZY73vzzpkbttP4yc/f/zJ5cOfcdAebf3tT9ax2+56pmofY21TebK2cRmcS78p5jXnOqf6WXvx+O2uZSyPyzy+7mD4a2+vTVfzHOYz1XM7Xb2zmdeU77X5NMv362yfp2GP2d7M65ZbbnmoqpYN6xtHoGdI29BLZ6pqPbAeYGJiojZv3jy3n/iNbrvTJua2/TSu3fbL30OvO+GYPdr625+sY7fd9UzVPsbapvJkbeMyOJd+U8xrznVO9bP24vHbXctYHpd5fN3B8NfeXpuu5jnMZ6rndrp6ZzOvKd9r82mW79fZPk/DHrO9mVeS/5qqbxxXuWwHju1bXwHcP4b9SpJmYRyBvhE4t7va5aXAT6tq5kM0SdJYzXjKJclngdOBo5JsB94DHARQVeuATcBZwCTwGHDefBUrSZrajIFeVefM0F/ARWOrSJI0J94pKkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNWKkQE9yZpJtSSaTXDqk/9lJvpDk1iRbk5w3/lIlSdOZMdCTLAGuBNYAq4FzkqweGHYRcEdVnQScDnwgycFjrlWSNI1RjtBPBSar6p6q2glcDawdGFPAs5IEeCbwY2DXWCuVJE1rlEBfDtzXt769a+t3BfBrwP3AbcDbq+qJwR0luSDJ5iSbd+zYMceSJUnDjBLoGdJWA+uvBr4LPA94EXBFksP22KhqfVVNVNXEsmXLZl2sJGlqowT6duDYvvUV9I7E+50HXFs9k8APgBPHU6IkaRSjBPrNwKokK7t/dJ4NbBwY80PglQBJngucANwzzkIlSdNbOtOAqtqV5GLgemAJsKGqtia5sOtfB1wGXJXkNnqnaC6pqofmsW5J0oAZAx2gqjYBmwba1vUt3w+8arylSZJmwztFJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDVipEBPcmaSbUkmk1w6xZjTk3w3ydYk3xhvmZKkmSydaUCSJcCVwBnAduDmJBur6o6+MYcDHwHOrKofJjl6vgqWJA03yhH6qcBkVd1TVTuBq4G1A2PeCFxbVT8EqKoHx1umJGkmowT6cuC+vvXtXVu/44HnJPl6kluSnDtsR0kuSLI5yeYdO3bMrWJJ0lCjBHqGtNXA+lLgFOA1wKuBv0py/B4bVa2vqomqmli2bNmsi5UkTW3Gc+j0jsiP7VtfAdw/ZMxDVfUo8GiSG4CTgLvGUqUkaUajHKHfDKxKsjLJwcDZwMaBMdcBv5lkaZJDgJcAd463VEnSdGY8Qq+qXUkuBq4HlgAbqmprkgu7/nVVdWeSLwFbgCeAj1XV7fNZuCTpqUY55UJVbQI2DbStG1h/P/D+8ZUmSZoN7xSVpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaMVKgJzkzybYkk0kunWbcbyT5RZI3jK9ESdIoZgz0JEuAK4E1wGrgnCSrpxh3OXD9uIuUJM1slCP0U4HJqrqnqnYCVwNrh4z7M+Aa4MEx1idJGtEogb4cuK9vfXvX9qQky4HXAuum21GSC5JsTrJ5x44ds61VkjSNUQI9Q9pqYP1DwCVV9YvpdlRV66tqoqomli1bNmqNkqQRLB1hzHbg2L71FcD9A2MmgKuTABwFnJVkV1X9y1iqlCTNaJRAvxlYlWQl8CPgbOCN/QOqauXu5SRXAf9qmEvSwpox0KtqV5KL6V29sgTYUFVbk1zY9U973lyStDBGOUKnqjYBmwbahgZ5Vf3x3pclSZot7xSVpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1IiRAj3JmUm2JZlMcumQ/jcl2dJ93ZjkpPGXKkmazoyBnmQJcCWwBlgNnJNk9cCwHwCnVdULgcuA9eMuVJI0vVGO0E8FJqvqnqraCVwNrO0fUFU3VtX/dKvfBlaMt0xJ0kxGCfTlwH1969u7tqm8DfjisI4kFyTZnGTzjh07Rq9SkjSjUQI9Q9pq6MDkFfQC/ZJh/VW1vqomqmpi2bJlo1cpSZrR0hHGbAeO7VtfAdw/OCjJC4GPAWuq6uHxlCdJGtUoR+g3A6uSrExyMHA2sLF/QJLnA9cCb66qu8ZfpiRpJjMeoVfVriQXA9cDS4ANVbU1yYVd/zrgr4EjgY8kAdhVVRPzV7YkadAop1yoqk3ApoG2dX3L5wPnj7c0SdJseKeoJDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUiJECPcmZSbYlmUxy6ZD+JPlw178lycnjL1WSNJ0ZAz3JEuBKYA2wGjgnyeqBYWuAVd3XBcBHx1ynJGkGoxyhnwpMVtU9VbUTuBpYOzBmLfCp6vk2cHiSY8ZcqyRpGktHGLMcuK9vfTvwkhHGLAce6B+U5AJ6R/AAjyTZNqtq9w1HAQ8tdhFj5Hz2ba3NB9qb00LP5wVTdYwS6BnSVnMYQ1WtB9aP8DP3WUk2V9XEYtcxLs5n39bafKC9Oe1L8xnllMt24Ni+9RXA/XMYI0maR6ME+s3AqiQrkxwMnA1sHBizETi3u9rlpcBPq+qBwR1JkubPjKdcqmpXkouB64ElwIaq2prkwq5/HbAJOAuYBB4Dzpu/khfdfn3KaAjns29rbT7Q3pz2mfmkao9T3ZKk/ZB3ikpSIwx0SWrEAR/oSTYkeTDJ7X1tJyX5VpLbknwhyWFd+8FJPtG135rk9K79kCT/luR7SbYmed8iTWcs8+nrW5/krm5er1+E6ZDk2CRfS3Jn99i+vWs/IsmXk9zdfX9O3zbv6j6GYluSV/e1n9LNdbL7qIphl9vuN/Pp69/Y/3wvtDE/R+d0z9GWJF9KctS+Pp8kR3bjH0lyRd9+Fj4XquqA/gJ+CzgZuL2v7WbgtG75rcBl3fJFwCe65aOBW+j9UjwEeEXXfjDw78Ca/XU+3frfAu/tlp8GHLVI8zkGOLlbfhZwF72PoPh74NKu/VLg8m55NXAr8HRgJfB9YEnXdxPwMnr3TXxxMZ6jcc6n638d8Jn+53t/nRO9izQe3P1a67b/m/1gPocCLwcuBK7o28+C58IBf4ReVTcAPx5oPgG4oVv+MrD76HQ18NVuuweBnwATVfVYVX2ta98JfIfetfgLbhzz6freCvxd1/dEVS3KnX1V9UBVfadb/hlwJ727kNcCn+yGfRL4/W55LXB1Vf1fVf2A3pVXp6b3URSHVdW3qvcO+1TfNgtmXPMBSPJM4M+B9y7cDPY0xjml+zq0++vpMBbhfpbZzqeqHq2qbwI/H9jPgufCAR/oU7gd+L1u+Q/45U1TtwJrkyxNshI4hafeUEWSw4HfpQvKfcSs5tPNAeCyJN9J8vkkz13YkveU5DjgxcB/As+t7l6H7vvR3bCpPoZiebc82L5o9nI+AJcBH6B3qfA+YW/mVFWPA38C3EYvyFcDH1+Qwqcw4nxG2c+C5IKBPtxbgYuS3ELvT66dXfsGei++zcCHgBuBXbs3SrIU+Czw4aq6Z0Ernt5s57OU3pHEf1TVycC3gH9Y6KL7dUej1wDvqKr/nW7okLaapn1R7O18krwI+NWq+ud5KXAOxjCng+gF+ouB5wFbgHeNvdARzWI+M+1nwXJhlM9yOeBU1feAVwEkOR54Tde+C3jn7nFJbgTu7tt0PXB3VX1o4aqd2Rzm8zC9o77dYfF54G0LWPJTdG/0a4BPV9W1XfN/Jzmmqh7oTqc82LVP9TEU23nqn7uL9vEUY5rPy4BTktxL7318dJKvV9XpCzGHQWOa04sAqur73T4/R+9c9YKb5XxmsmC54BH6EEmO7r4/DfhLYF23fkiSQ7vlM4BdVXVHt/5e4NnAOxal6GnMdj7dOeYvAKd3u3glcMdC193VFXp/dt9ZVf/Y17UReEu3/Bbgur72s5M8vTuNtAq4qfsT+WdJXtrt89y+bRbMGOfz0ap6XlUdR+8fcnctYpiPZU7Aj4DVSZZ1486gd/56Qc1hPtPta2FzYT7/47o/fNH7U+gB4HF6Rw5vA95O7z/bdwHv45d31B4HbKP3IvsK8IKufQW9P9/vBL7bfZ2/v86n63sBvX+kbqF33u/5izSfl3eP7Za+x/Ys4Miurru770f0bfNueldObKPvqgJ6//C9veu7YvfjsL/Op6//OBb3KpdxPkcXdq/HLfQOKo7cT+ZzL72LER7p3nerFyMXvPVfkhrhKRdJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhrx/4qcioL4+gCLAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAASzUlEQVR4nO3dfbBcd13H8feHpFVbqAWaMiUptI6BGlGgvRZwUKoO0OBoFHVs6xgep1OHOvgwDmV8ZMqMzw6DrWSiVh5UOipVg0YqMgpKqTZlIDRtU0KBNrTYYOWpCGnK1z/OubDZ7N7dm+7dm/z6fs3c6Z7z+53f/X7v3vPJ2b2721QVkqTj36NWuwBJ0mwY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBroeEZI8N8kNST6X5P4k70/yXQPjJyf5YpKdI479RJKDSU4b2v+hJJXkrJXvQJrMQFfzkpwC/APwh8DjgPXA64CvDEz78X77BUnOGLHMx4GLB9b8DuCbVqpm6WgY6HokeApAVb29qh6qqv+rqn+uqt0Dc14CbAN2Az81Yo23AVuH5r91pQqWjoaBrkeCO4CHkrwlyeYkjx0cTPIk4ALgL/qvrUcuwY3AKUm+Lcka4CeBP1/ZsqXlMdDVvKr6PPBcoIA/Bg4k2ZHkCf2UrcDuqroVeDvw7UmeOWKpxav05wO3A59a8eKlZTDQ9YhQVbdV1UuragPwNOCJwBv64a10V+ZU1T3Ae+meUhn2NuAS4KX4dIuOQQa6HnGq6nbgzcDTknw3sBF4bZJPJ/k08Czg4iRrh477JN0fR18EXDffqqXJDHQ1L8k5SX4xyYZ++0y6V6zcSHcl/m5gE/CM/utpwEnA5hHLvQL4/qp6YB61S8uxdvIU6bj3Bbqr7l9IcirwWbqXMf4S8Elga1V9evCAJG+jC/t3Du6vqo/NpWLpKMT/wYUktcGnXCSpERMDPck1Se5LcsuY8SR5Y5J9SXYnOXf2ZUqSJpnmCv3NwIVLjG+me5XARuBS4E0PvyxJ0nJNDPSqeh9w/xJTtgBvrc6NwKljPgtDkrSCZvEql/XA3QPb+/t99w5PTHIp3VU8J5988nnnnHPO0X3HL3yp++9jTuKzX36QU7/xhMOGP/vlB792e3hseN6p33jCYesdMTZlHUfsf8xJh48v3h41f/i4EX0sWqxpbI/9Gkv9XAbXGNf/YePjDPc3rq9Jx46ob7jWwfHhHobnDc8dPmZiH0N1jTJqfYBTH3xw/DpTrj22xhHHjet/yXVG1TDN78Ck+3j4937MeouG79tRv9uj5h7x/cZtj/m+i+st9bM74twY6GWq4wZrmvLYad18882fqap1o8ZmEegZsW/kS2eqajuwHWBhYaF27dp1dN/xvf1xz1vgur338uKnHv6A4Lq9X/+3ZHhseN6Ln3rGYesdMTZlHUfsf97C4ePvHeh1eP7wcSP6WLRY09ge+zWW+rkMrjGu/8PGxxnub1xfk44dUd9wrYPjwz0MzxueO3zMxD6G6hpl1PoAL/70p8avM+XaY2sccdy4/pdcZ1QN0/wOTLqPh3/vx6y3aPi+HfW7PWruEd9v3PaY77u43lI/uyPOjYFepjpusKYpj51Wkk+OG5vFq1z2A2cObG8A7pnBupKkZZhFoO8Atvavdnk28LmqGn35IklaMROfcknydrqPFj0tyX7g14ETAKpqG7CT7rMt9gFfAl62UsVKksabGOhVdfGE8QJeNbOKJElHxXeKSlIjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjZgq0JNcmGRvkn1Jrhgx/s1J3pnkw0n2JHnZ7EuVJC1lYqAnWQNcDWwGNgEXJ9k0NO1VwK1V9XTgAuD3k5w441olSUuY5gr9fGBfVd1ZVQeBa4EtQ3MKeEySAI8G7gcOzbRSSdKSpgn09cDdA9v7+32DrgK+DbgH+Ajw6qr66vBCSS5NsivJrgMHDhxlyZKkUaYJ9IzYV0PbLwQ+BDwReAZwVZJTjjioantVLVTVwrp165ZdrCRpvGkCfT9w5sD2Bror8UEvA66rzj7g48A5sylRkjSNaQL9JmBjkrP7P3ReBOwYmnMX8AMASZ4APBW4c5aFSpKWtnbShKo6lORy4HpgDXBNVe1Jclk/vg24Enhzko/QPUXzmqr6zArWLUkaMjHQAapqJ7BzaN+2gdv3AC+YbWmSpOXwnaKS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGjFVoCe5MMneJPuSXDFmzgVJPpRkT5L3zrZMSdIkaydNSLIGuBp4PrAfuCnJjqq6dWDOqcAfARdW1V1JTl+pgiVJo01zhX4+sK+q7qyqg8C1wJahOZcA11XVXQBVdd9sy5QkTTJNoK8H7h7Y3t/vG/QU4LFJ/i3JzUm2jlooyaVJdiXZdeDAgaOrWJI00jSBnhH7amh7LXAe8IPAC4FfTfKUIw6q2l5VC1W1sG7dumUXK0kab+Jz6HRX5GcObG8A7hkx5zNV9QDwQJL3AU8H7phJlZKkiaa5Qr8J2Jjk7CQnAhcBO4bm/D3wPUnWJjkJeBZw22xLlSQtZeIVelUdSnI5cD2wBrimqvYkuawf31ZVtyV5F7Ab+CrwJ1V1y0oWLkk63DRPuVBVO4GdQ/u2DW3/LvC7sytNkrQcvlNUkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqxFSBnuTCJHuT7EtyxRLzvivJQ0l+fHYlSpKmMTHQk6wBrgY2A5uAi5NsGjPvt4HrZ12kJGmyaa7Qzwf2VdWdVXUQuBbYMmLezwLvAO6bYX2SpClNE+jrgbsHtvf3+74myXrgR4FtSy2U5NIku5LsOnDgwHJrlSQtYZpAz4h9NbT9BuA1VfXQUgtV1faqWqiqhXXr1k1boyRpCmunmLMfOHNgewNwz9CcBeDaJACnAS9Kcqiq/m4mVUqSJpom0G8CNiY5G/gUcBFwyeCEqjp78XaSNwP/YJhL0nxNDPSqOpTkcrpXr6wBrqmqPUku68eXfN5ckjQf01yhU1U7gZ1D+0YGeVW99OGXJUlaLt8pKkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRUwV6kguT7E2yL8kVI8Z/Ksnu/uuGJE+ffamSpKVMDPQka4Crgc3AJuDiJJuGpn0ceF5VfSdwJbB91oVKkpY2zRX6+cC+qrqzqg4C1wJbBidU1Q1V9b/95o3AhtmWKUmaZJpAXw/cPbC9v983ziuAfxo1kOTSJLuS7Dpw4MD0VUqSJpom0DNiX42cmHwfXaC/ZtR4VW2vqoWqWli3bt30VUqSJlo7xZz9wJkD2xuAe4YnJflO4E+AzVX1P7MpT5I0rWmu0G8CNiY5O8mJwEXAjsEJSZ4EXAf8dFXdMfsyJUmTTLxCr6pDSS4HrgfWANdU1Z4kl/Xj24BfAx4P/FESgENVtbByZUuShk3zlAtVtRPYObRv28DtVwKvnG1pkqTl8J2iktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSI6YK9CQXJtmbZF+SK0aMJ8kb+/HdSc6dfamSpKVMDPQka4Crgc3AJuDiJJuGpm0GNvZflwJvmnGdkqQJprlCPx/YV1V3VtVB4Fpgy9CcLcBbq3MjcGqSM2ZcqyRpCWunmLMeuHtgez/wrCnmrAfuHZyU5FK6K3iALybZO+L7nQZ8Zoq6jnUt9NFCD2Afx5IWeoDV7ePJ4wamCfSM2FdHMYeq2g5sX/KbJbuqamGKuo5pLfTRQg9gH8eSFnqAY7ePaZ5y2Q+cObC9AbjnKOZIklbQNIF+E7AxydlJTgQuAnYMzdkBbO1f7fJs4HNVde/wQpKklTPxKZeqOpTkcuB6YA1wTVXtSXJZP74N2Am8CNgHfAl42cOoacmnZI4jLfTRQg9gH8eSFnqAY7SPVB3xVLck6TjkO0UlqREGuiQ1Yi6BnuSaJPcluWVg39OTfCDJR5K8M8kp/f4Tkryl339bktcOHHNxv393knclOW0e9R9FDycm+bN+/4eTXNDvPynJPya5PcmeJL81r/pn2cfA2PYkd/T9/Ngcezgzyb/2vx97kry63/+4JO9O8tH+v48dOOa1/UdT7E3ywoH95/X97es/vmLUS3CP+T4GxncM3rfHUw+rfH4vq48kj+/nfzHJVQPrrO45XlUr/gV8L3AucMvAvpuA5/W3Xw5c2d++BLi2v30S8AngLLo/4N4HnNaP/Q7wG/Oo/yh6eBXwZ/3t04Gb6f7xPAn4vn7/icC/A5vn1cOs+ui3Xwe8vr/9qMX7ZU49nAGc299+DHAH3cdS/A5wRb//CuC3+9ubgA8D3wCcDXwMWNOP/RfwHLr3UvzTPO+PWfbRj78Y+MvB+/Z46eEYOL+X28fJwHOBy4CrBtZZ1XN8LlfoVfU+4P6h3U8F3tfffjeweIVXwMlJ1gLfBBwEPk93wqUfC3AKc3yt+zJ72AS8pz/uPuCzwEJVfamq/rXffxD4IN1r9udmFn30Yy8HfrMf+2pVze1dc1V1b1V9sL/9BeA2uncmbwHe0k97C/Aj/e0tdBcJX6mqj9O9Guv8dB9PcUpVfaC6M/CtA8ccN30AJHk08AvA6+dVf1/3rHpY7fN7WX1U1QNV9R/Al4fWWdVzfDWfQ78F+OH+9k/w9Tcm/Q3wAN3HBtwF/F5V3V9VDwI/A3yE7o7eBPzpXCs+0rgePgxsSbI2ydnAeRz+xiuSnAr8EH1grrJl9dHXDnBlkg8m+eskT5hvyZ0kZwHPBP4TeEL173/o/3t6P23cR1Os728P75+7h9kHwJXA79O9bHhVPJwejqXze8o+plln7uf4agb6y4FXJbmZ7iHOwX7/+cBDwBPpHpL9YpJvSXIC3R3+zH5sN/DaI1adr3E9XEP3i7oLeANwA3Bo8aD+0cfbgTdW1Z1zrXi05faxlu6q4/1VdS7wAeD35l10f1X6DuDnqurzS00dsa+W2D9XD7ePJM8AvrWq/nZFCpzCDHo4Js7vZfQxaZ1VOcen+SyXFVFVtwMvAEjyFOAH+6FLgHf1/2Lfl+T9dA/zH98f97H+mL+ie05r1YzroaoOAT+/OC/JDcBHBw7dDny0qt4wv2rHO4o+/ofuSnAxQP4aeMUcS6YPgHcAf1FV1/W7/zvJGVV1b/90yn39/nEfTbGfwx8Oz/0jK2bUx3OA85J8gu6cPj3Jv1XVBcdRD8+A1T2/l9nHJKtyjq/aFXqS0/v/Pgr4FWBbP3QX8P3pnAw8G7gd+BSwKcm6ft7z6Z7nWjXjeuj/0n1yf/v5wKGqurXffj3wzcDPrUrRIyy3j/755ncCF/RL/ABw6xzrDd3D8duq6g8GhnYAL+lvvwT4+4H9FyX5hv6po43Af/UPob+Q5Nn9mlsHjllxM+zjTVX1xKo6i+4PdXfMMcxn0gOrfH4fRR9LrbV65/g8/vJK99DjXuBBun+hXwG8mu4vyXcAv8XX37X6aLorvj10IfFLA+tcRncn76YLlMfPo/6j6OEsYG9f678AT+73b6B7SH8b8KH+65Xz6mFWffRjT6b7Q+puuucInzTHHp7b/xx3D/wcX0T3KO49dI8i3gM8buCYX6Z7RcVeBl51QPfo75Z+7KrF3o+3PgbGz2K+r3KZ5X2xmuf30fTxCboXGHyxP5c2rfY57lv/JakRvlNUkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RG/D96/Tdkr07BmwAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "for driver in driver_indices:\n", " la_nina = phases[driver][\"+ve\"]\n", " el_nino = phases[driver][\"-ve\"]\n", " for i, group in la_nina.groupby(np.cumsum(la_nina != la_nina.shift())):\n", " if not group.iloc[0]:\n", " continue\n", "\n", " start = group.index[0]\n", " end = group.index[-1]\n", " plt.axvspan(start, end, facecolor=\"lightblue\")\n", "\n", " for i, group in el_nino.groupby(np.cumsum(el_nino != el_nino.shift())):\n", " if not group.iloc[0]:\n", " continue\n", "\n", " start = group.index[0]\n", " end = group.index[-1]\n", " plt.axvspan(start, end, facecolor=\"pink\")\n", " plt.title(driver.upper())\n", " plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Cumulative distributions of rainfall for each phase\n", "\n", "For the positive, negative, and neutral phase, collect all rainfall observations that occurred during that phase. We can then treat these rainfall observations as samples of a random variable, conditioned on the phase. This allows us to build the cumulative distribution function (CDF), which describes the probability distribution of the rainfall." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "# Figure out what times the phases were active.\n", "times = {}\n", "for driver, ps in phases.items():\n", " times[driver] = {}\n", " for name, ph in ps.items():\n", " times[driver][name] = rolling_index[driver].index[ph]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Estimate the mean and CDF for each phase.\n", "means = {}\n", "cdfs = {}\n", "for driver, times_ in tqdm(times.items()):\n", " means[driver] = {}\n", " cdfs[driver] = {}\n", " for name, time in times_.items():\n", " means[driver][name] = np.mean(deseasonalised_rainfall.sel(time=time), axis=0)\n", " cdfs[driver][name] = np.sort(\n", " np.nan_to_num(deseasonalised_rainfall.sel(time=time)), axis=0\n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We need to line up the CDFs despite having a different number of days. We can interpolate to do this:" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "max_rainfall = deseasonalised_rainfall.max().values\n", "min_rainfall = 0" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "interps = {}\n", "for driver, cdfs_ in tqdm(cdfs.items()):\n", " interps[driver] = {}\n", " for name, cdf in cdfs_.items():\n", " interps[driver][name] = []\n", " for y in range(cdf.shape[1]):\n", " for x in range(cdf.shape[2]):\n", " interps[driver][name].append(\n", " # This interpolation essentially rescales each set of observations to a consistent x scale\n", " # and therefore approximates the CDF. Basically converting the units from \"daily\" to \"quantile\".\n", " # The 200 controls the resolution of the resulting interpolated CDF. Doesn't need to be too high.\n", " np.interp(\n", " np.linspace(min_rainfall, max_rainfall, 200),\n", " cdf[:, y, x],\n", " np.linspace(0, 1, len(cdf)),\n", " left=0,\n", " right=1,\n", " )\n", " )" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "# Convert the interpolated CDFs back into numpy arrays of the appropriate shape for the next step.\n", "for driver, interps_ in interps.items():\n", " for name, interp in interps_.items():\n", " interps[driver][name] = np.array(interp).reshape(cdf.shape[1:] + (-1,))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Difference between active phases compared to neutral\n", "\n", "We can characterise the influence of climate drivers over the rainfall by a) assuming that there is no confounding correlation between climate driver and rainfall, and then b) calculating the distance between the above probability distributions. There are many, many ways of comparing probability distributions. We will employ two:\n", "\n", "1. The Kolmogorov-Smirnov (KS) distance, commonly used to determine whether two distributions are different in a test known as the \"KS test\". Arguably the most common statistical test, besides the Student t-test. While it is parameter-free and works on any distribution, it's also fairly weak and hence tends to under-predict differences in distributions, i.e. if it suggests that two distributions are different they probably are, but if it suggests that they are the same it might not be.\n", "2. The sum-of-squares difference, or the Euclidean distance. This is a measure of how far apart the quantiles are, and is a good choice if you expect your quantiles to have normally-distributed noise." ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "# Kolmogorov-Smirnov.\n", "ks = {}\n", "for driver, cdf in interps.items():\n", " ks[driver] = {}\n", " ks[driver][\"+ve\"] = abs(cdf[\"+ve\"] - cdf[\"neutral\"]).max(axis=-1)\n", " ks[driver][\"-ve\"] = abs(cdf[\"-ve\"] - cdf[\"neutral\"]).max(axis=-1)\n", "\n", "# Euclidean.\n", "euc = {}\n", "for driver, cdf in interps.items():\n", " euc[driver] = {}\n", " euc[driver][\"+ve\"] = np.sqrt(np.mean((cdf[\"+ve\"] - cdf[\"neutral\"]) ** 2, axis=-1))\n", " euc[driver][\"-ve\"] = np.sqrt(np.mean((cdf[\"-ve\"] - cdf[\"neutral\"]) ** 2, axis=-1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then we can plot the climate driver influence on rainfall based on KS or Euclidean distance." ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAABDAAAALICAYAAACJhQBYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzde7hcZXn///fHkAhUENBYkUDhyy/6kwqGGA71VEWRBJWo/QqhVpBqUwTqqfVrsK2CVuTXYrVRJIJSoaUCBQ9RY/HwFVsL2ASMqZxKRA4xqUTAgEUOkfv3x6zEYbOTzCSz98ze+/26rrlm1lrPs9a9NnrNnXue9TypKiRJkiRJkgbZE/odgCRJkiRJ0pZYwJAkSZIkSQPPAoYkSZIkSRp4FjAkSZIkSdLAs4AhSZIkSZIGngUMSZIkSZI08CxgSONcktOS/GO/4xirklyf5CX9jkOSJG1ZkpckWdW2vcnv8aFtJQ0+CxjSGJLktiQvb9uel+TeJL/bz7gGTZI3J7kpyf1Jfprkq0l22ppzVdVvV9WVPQ5RkiSxMbf5ZZJftL0+0avz+z0ujS/b9TsASVsnyfHA3wKvrKqr+h3PtkiyXVWt79G5fhc4A5hdVd9Pshvw6l6ce5hr9SxuSZImsFdX1Tf7HYSkwecIDGkMSjIf+AhwRFVdleQZSRYnuSfJyiR/tIl+eyepJCckubMZvXFikoOSrEjy8/ZfPZI8IclfJLk9yV1JLkzy5LbjxzXH7k7yl+0jRJI8McnHkqxuXh9L8sTm2EuSrEryniT/Dfz9FtrfmORVbdfdLsnPkswc5jYPAq6uqu8DVNU9VXVBVd3f9P1skk8m+VrzK8+/J3l6c717m5EbB7Zdq/2eTktyWZJ/THIf8KYkVyb5YHOe+5N8PclT2/of1Qxf/XnT9tnN/gVJLhvy3+fvkizs5H8DkiSNZ0MfgW3LYbZrtndL8vdNznBvki9u4jzt3+M7NHnAvUluoJUztLd9RpLLk6xN8uMkb2s7dnCSq5vv8zVJPpFkStvxanKqW5rzn50kPf6zSBOeBQxp7Hkr8EHgZVW1rNn3OWAV8AzgfwNnJHnZZs5xCDAdOAb4GPDnwMuB3waOzq8fSXlT83op8L+AJwGfAEiyH/BJ4A3A7sCTgT3arvHnwKHADOC5wMHAX7QdfzqwG/BbwPwttP8ccGxb3yOAn1XVdcPc2/eAI5KcnuQFG4ogQxzdnPupwEPA1cB1zfZltEa2bMrcps0uwEXNvt8HTgCeBkwB/gwgyTOb2N8BTAWWAF9uEp7PAUcm2blpO6mJ6582c21JktTyD8COtHKXpwEf7aDP+4F9m9cRwPEbDiR5AvBl4Ae08pmXAe9IckTT5FfAO2nlCr/THD9pyPlfRaso8lxa3+lHIKmnLGBIY8/hwDXAfwIk2RN4IfCeqnqwqpYDnwbeuJlzfLBp+3Xgf4DPVdVdVfUT4N+ADSMQ3gD8bVXdWlW/AE4F5jW/fvxv4MtV9d2qehh4H1Bt13gD8IHmvGuB04fE9Cjw/qp6qKp+uYX2/wQclWTHZvv32cQ/9Kvq34DXATOBrwJ3J/nbpkCwwReq6tqqehD4AvBgVV1YVb8CLmm7/+FcXVVfrKpHm7gB/r6q/qvZvpRWEQZaBaKvVtU3quoR4CxgB+D5VXU7raLJa5q2hwEPVNU1m7m2JEnj0RebkQ0bXsOOJN0gye7AHODEqrq3qh6pqu90cJ2jgQ81ozPvBNpHPR4ETK2qD1TVw1V1K3AeMA+gyRuuqar1VXUb8Clg6BxkZ1bVz6vqDuDb/DofkNQjFjCksedE4JnAp5uhic8A7tnwiETjdh47GmKon7Z9/uUw209qPj+jOVf7ebcDfrM5dueGA1X1AHB3W9vh+j6jbXttU0DYYvuqWgncCLy6KWIcRVPAyGMn/dqraf+1qno1rREec2mNInnLVtz/cO4cZt9/t31+gE38/arq0ab/hv82/8SvR5ZssigjSdI495qq2qXtdd4W2u9JK/e5t8vrPCZ34bF5x28Bz2gvpADvpZXzkOSZSb6S5L+bx0jPoDUao92m8gFJPWIBQxp77qI1bPFFtB7hWA3slseusrEX8JMeXGs1rS/09vOup/UP/jXAtA0HkuwAPGULfVe3bbeP1uik/YbHSOYCNzRFDarqSW2vO9pP2IyS+Bbwf4HnbP5WOzY07s15zD01Bac9+fV/m38GXpJkGvBaLGBIkrTB/9B6RGSDp7d9vpNW7rNLl+dcQ+t7eIO9hpzzx0MKKTtV1ZHN8XOAm4DpVbUzreKGc1xIo8wChjQGVdVqWo8czAbeBVwFfDjJ9kkOAN7Mr+dn2BafA96ZZJ8kT6L1a8Mlzcobl9EaEfH8Zk6H03nsF/nngL9IMrWZ1PJ9wD+yaVtqfzHwClpzgGzyH/pJ5qa1vOyuaTmY1hDPfjyacSnwyiQvSzIZ+FNac25cBdA8KnMl8Pe0kqYb+xCjJEmDaDnw4iR7pTWB+KkbDlTVGuBrwCeb7/vJSV7cwTkvBU5t+kwD/qTt2H8A96U1wfgOSSYleU6SDRN97gTcB/wiyf9LKx+RNMosYEhjVPPs5mG05qJYCexN6xf/L9CaW+IbPbjM+bQmyfpX4MfAgzRf9lV1ffP5Ylq/aNxPa3TIQ03fvwKWAStozddxXbNvUzbbvklWrgaeT2ueik25F/gj4BZaicY/An9TVb0o6HSlqm4G/gD4OPAzWsu5vrqZM2SDf6I1gaqjLyRJE9WXhzwS+oUmj7mEVl5wLfCVIX3eCDxCa1TEXbQmzN6S02k9NvJj4Ou0chwAmnmwXk1r3oof0/re/jStScqhNUH379PKd85j87mIpBGSqm5GQ0vS8JoRGj+nNbTyx/2OR5IkSdL44ggMSVstyauT7JjkN2itsPGfwG39jUqSJEnSeGQBQ9K2mEvrsZXVwHRgXjmsS5IkSdII8BESSZIkSZI08ByBIUmSJEmSBt52/Q6gn5761KfW3nvv3e8wJEka06699tqfVdXUfsfRb+YVkiT1xqZyiwldwNh7771ZtmxZv8OQJGlMS3J7v2MYBOYVkiT1xqZyCx8hkSRJkiRJA68vBYwks5PcnGRlkgXDHE+Shc3xFUlmDjk+Kcn3k3ylbd9uSb6R5JbmfdfRuBdJkiRJkjTyRr2AkWQScDYwB9gPODbJfkOazaG1JON0YD5wzpDjbwduHLJvAfCtqpoOfKvZliRJkiRJ40A/5sA4GFhZVbcCJLkYmAvc0NZmLnBhtdZ4vSbJLkl2r6o1SaYBrwQ+BLxrSJ+XNJ8vAK4E3jOSNyJJGlseeeQRVq1axYMPPtjvUMak7bffnmnTpjF58uR+hyJJ0kAwt9g23eYW/Shg7AHc2ba9CjikgzZ7AGuAjwH/B9hpSJ/frKo1AE2h42nDXTzJfFqjOthrr7228hYkSWPRqlWr2Gmnndh7771J0u9wxpSq4u6772bVqlXss88+/Q5HkqSBYG6x9bYmt+jHHBjD/VetTtokeRVwV1Vdu7UXr6pzq2pWVc2aOnXCr/gmSRPKgw8+yFOe8hQTjK2QhKc85Sn+wiRJUhtzi623NblFPwoYq4A927anAas7bPMC4KgktwEXA4cl+cemzU+T7A7QvN/V+9AlSWOdCcbW828nSdLj+f249br92/WjgLEUmJ5knyRTgHnA4iFtFgPHNauRHAqsq6o1VXVqVU2rqr2bfv+3qv6grc/xzefjgS+N+J1IkiRJkqRRMepzYFTV+iSnAFcAk4Dzq+r6JCc2xxcBS4AjgZXAA8AJHZz6TODSJG8G7gBePxLxS5LGjy/vu29Pz/fqH/1oi22S8K53vYuPfOQjAJx11ln84he/4LTTTutpLGeccQbvfe97N24///nP56qrrurpNSRJ0mOZW4ysfozAoKqWVNUzq2rfqvpQs29RU7ygWk5uju9fVcuGOceVVfWqtu27q+plVTW9eb9n9O5IkqTOPPGJT+Tzn/88P/vZz0b0OmecccZjti1eSJI0Pk2k3KIvBQxJkiaq7bbbjvnz5/PRj370ccfWrl3L7/3e73HQQQdx0EEH8e///u8b9x9++OHMnDmTP/7jP+a3fuu3NiYpr3nNa3je857Hb//2b3PuuecCsGDBAn75y18yY8YM3vCGNwDwpCc9CYBjjjmGJUuWbLzmm970Ji6//HJ+9atf8e53v5uDDjqIAw44gE996lMj+neQJEm9MZFyCwsYkiSNspNPPpmLLrqIdevWPWb/29/+dt75zneydOlSLr/8ct7ylrcAcPrpp3PYYYdx3XXX8drXvpY77rhjY5/zzz+fa6+9lmXLlrFw4ULuvvtuzjzzTHbYYQeWL1/ORRdd9JhrzJs3j0suuQSAhx9+mG9961sceeSRfOYzn+HJT34yS5cuZenSpZx33nn8+Mc/HuG/hCRJ6oWJkluM+hwYkiRNdDvvvDPHHXccCxcuZIcddti4/5vf/CY33HDDxu377ruP+++/n+9+97t84QtfAGD27NnsuuuuG9ssXLhw47E777yTW265hac85SmbvPacOXN429vexkMPPcS//Mu/8OIXv5gddtiBr3/966xYsYLLLrsMgHXr1nHLLbd0vC67JEnqn4mSW1jAkCSpD97xjncwc+ZMTjjh1/NUP/roo1x99dWPSTwAqmrYc1x55ZV885vf5Oqrr2bHHXfkJS95yRbXUt9+++15yUtewhVXXMEll1zCscceu/EaH//4xzniiCO28c4kSVI/TITcwgLGOHfEB7/a7xA2uuIvX9nvECRpYOy2224cffTRfOYzn+EP//APAXjFK17BJz7xCd797ncDsHz5cmbMmMELX/hCLr30Ut7znvfw9a9/nXvvvRdo/ZKx6667suOOO3LTTTdxzTXXbDz/5MmTeeSRR5g8efLjrj1v3jw+/elPs2zZMj772c8CcMQRR3DOOedw2GGHMXnyZP7rv/6LPfbYg9/4jd8Y4b+EpPGu16sybItOVnSQxqqJkFtYwJAkTVj9TmT/9E//lE984hMbtxcuXMjJJ5/MAQccwPr163nxi1/MokWLeP/738+xxx7LJZdcwu/+7u+y++67s9NOOzF79mwWLVrEAQccwLOe9SwOPfTQjeeaP38+BxxwADNnznzcs6qveMUrOO644zjqqKOYMmUKAG95y1u47bbbmDlzJlXF1KlT+eIXvzg6fwiNGYP0wwj444ikwWNuMbK5RTY1dGQimDVrVi1b9rgVWseVQUo0TDLGhkH6lQT6/yWg8eXGG2/k2c9+dr/D6NpDDz3EpEmT2G677bj66qt561vfyvLly/sSy3B/wyTXVtWsvgQ0QMwrRp+5xdgwSLmFeYV6zdxi23WTWzgCQ5KkAXfHHXdw9NFH8+ijjzJlyhTOO++8fockSZLGsLGaW1jAkCRpwE2fPp3vf//7/Q5DkiSNE2M1t3hCvwOQJEnaFklmJ7k5ycokC4Y5niQLm+MrkszcUt8kM5Jck2R5kmVJDh6t+5EkScOzgCFJksasJJOAs4E5wH7AsUn2G9JsDjC9ec0Hzumg718Dp1fVDOB9zbYkSeojCxiSJGksOxhYWVW3VtXDwMXA3CFt5gIXVss1wC5Jdt9C3wJ2bj4/GVg90jciSZI2zzkwJEnSWLYHcGfb9irgkA7a7LGFvu8ArkhyFq0ffJ4/3MWTzKc1qoO99tpr6+5AkiR1xAKGJGnC6vWSkJ0s6Thp0iT2339/1q9fz7Of/WwuuOACdtxxx46vsXr1at72trdx2WWXsXz5clavXs2RRx4JwOLFi7nhhhtYsOBx00CMZxlm39A14jfVZnN93wq8s6ouT3I08Bng5Y9rXHUucC60llHtNGhJ0vhkbjGyfIREkqRRtMMOO7B8+XJ++MMfMmXKFBYtWtRV/2c84xlcdtllACxfvpwlS5ZsPHbUUUcNTIIxilYBe7ZtT+Pxj3tsqs3m+h4PfL75/M+0HjeRJGngTKTcwgKGJEl98qIXvYiVK1dyzz338JrXvIYDDjiAQw89lBUrVgDwne98hxkzZjBjxgwOPPBA7r//fm677Tae85zn8PDDD/O+972PSy65hBkzZnDJJZfw2c9+llNOOYV169ax99578+ijjwLwwAMPsOeee/LII4/wox/9iNmzZ/O85z2PF73oRdx00039/BP0wlJgepJ9kkwB5gGLh7RZDBzXrEZyKLCuqtZsoe9q4Hebz4cBt4z0jUiStK3Ge25hAUOSpD5Yv349X/va19h///15//vfz4EHHsiKFSs444wzOO644wA466yzOPvss1m+fDn/9m//xg477LCx/5QpU/jABz7AMcccw/LlyznmmGM2Hnvyk5/Mc5/7XL7zne8A8OUvf5kjjjiCyZMnM3/+fD7+8Y9z7bXXctZZZ3HSSSeN7o33WFWtB04BrgBuBC6tquuTnJjkxKbZEuBWYCVwHnDS5vo2ff4I+EiSHwBn0MxzIUnSoJoIuYVzYEiSNIp++ctfMmPGDKD1K8mb3/xmDjnkEC6//HIADjvsMO6++27WrVvHC17wAt71rnfxhje8gde97nVMmzat4+scc8wxXHLJJbz0pS/l4osv5qSTTuIXv/gFV111Fa9//es3tnvooYd6e4N9UFVLaBUp2vctavtcwMmd9m32fxd4Xm8jlSSp9yZSbmEBQ5KkUbThOdV2rX9fP1YSFixYwCtf+UqWLFnCoYceyje/+U223377jq5z1FFHceqpp3LPPfdw7bXXcthhh/E///M/7LLLLo+7viRJGrsmUm7hIySSJPXZi1/8Yi666CIArrzySp761Key884786Mf/Yj999+f97znPcyaNetxz5TutNNO3H///cOe80lPehIHH3wwb3/723nVq17FpEmT2Hnnndlnn33453/+Z6CV3PzgBz8Y2ZuTJEmjbrzmFo7AkCRNWJ0sTTYaTjvtNE444QQOOOAAdtxxRy644AIAPvaxj/Htb3+bSZMmsd9++zFnzhzWrFmzsd9LX/pSzjzzTGbMmMGpp576uPMec8wxvP71r+fKK6/cuO+iiy7irW99K3/1V3/FI488wrx583juc5874vcoSdJEYG4xsrlFhhtaMlHMmjWrli1b1u8wRlSv1yHeFoPyf2Zt3pf33bffITzGq3/0o36HoHHkxhtv5NnPfna/wxjThvsbJrm2qmb1KaSBYV4x+swtxoZByi3MK9Rr5hbbrpvcwhEYkqRxZ1PJ8jM+9Sl+vn79qMayy/77j+r1JEmSxivnwJAkSZIkSQPPAoYkaeJ49NFhZ+VWZ/zbSZL0eH4/br1u/3YWMCRJE8Yjd9zBfY88YqKxFaqKu+++u+Ol1iRJmgi233577r77bnOLrbA1uYVzYEiSJoy7zz4bTj6Zn+21FzxhdGr4O243fr5qt99+e6ZNm9bvMCRJGhjTpk1j1apVrF27tt+hjEnd5hbjJ6uSJGkLHr3vPtZ++MOjek1nvJckafyaPHky++yzT7/DmDB8hESSJEmSJA08CxiSJEmSJGng9aWAkWR2kpuTrEyyYJjjSbKwOb4iycxm//ZJ/iPJD5Jcn+T0tj6nJflJkuXN68jRvCdJkiRJkjRyRn0OjCSTgLOBw4FVwNIki6vqhrZmc4DpzesQ4Jzm/SHgsKr6RZLJwHeTfK2qrmn6fbSqzhqte5EkSZIkSaOjHyMwDgZWVtWtVfUwcDEwd0ibucCF1XINsEuS3ZvtXzRtJjcv16uRJEmSJGmc60cBYw/gzrbtVc2+jtokmZRkOXAX8I2q+l5bu1OaR07OT7LrcBdPMj/JsiTLXOpGkiRJkqSxoR8FjAyzb+goik22qapfVdUMYBpwcJLnNMfPAfYFZgBrgI8Md/GqOreqZlXVrKlTp25N/JIkSZIkaZT1o4CxCtizbXsasLrbNlX1c+BKYHaz/dOmuPEocB6tR1UkSZIkSdI40I8CxlJgepJ9kkwB5gGLh7RZDBzXrEZyKLCuqtYkmZpkF4AkOwAvB25qtndv6/9a4IcjfSOSJEmSJGl0jPoqJFW1PskpwBXAJOD8qro+yYnN8UXAEuBIYCXwAHBC03134IJmJZMnAJdW1VeaY3+dZAatR01uA/54lG5JkiRJkiSNsFEvYABU1RJaRYr2fYvaPhdw8jD9VgAHbuKcb+xxmJIkSZIkaUD04xESSZIkSZKkrvRlBIYkdeqID3613yFsdMVfvrLfIUiSJEkTliMwJEmSJEnSwHMERo99ed99+x3CY71pYb8jkCRJkiRpmzkCQ5IkSZIkDTxHYEiSJEmStAmDNsr+1T/6Ub9D6BtHYEiSJEmSpIFnAUOSJI1pSWYnuTnJyiQLhjmeJAub4yuSzNxS3ySXJFnevG5Lsny07keSJA3PR0gkSdKYlWQScDZwOLAKWJpkcVXd0NZsDjC9eR0CnAMcsrm+VXVM2zU+AqwblRuSJEmb5AgMSZI0lh0MrKyqW6vqYeBiYO6QNnOBC6vlGmCXJLt30jdJgKOBz430jUiSpM2zgCFJksayPYA727ZXNfs6adNJ3xcBP62qW4a7eJL5SZYlWbZ27dqtCF+SJHXKAoYkSRrLMsy+6rBNJ32PZTOjL6rq3KqaVVWzpk6dutlAJUnStnEODEmSNJatAvZs254GrO6wzZTN9U2yHfA64Hk9jFeSJG0lR2BIkqSxbCkwPck+SaYA84DFQ9osBo5rViM5FFhXVWs66Pty4KaqWjXytyFJkrbEERiSJGnMqqr1SU4BrgAmAedX1fVJTmyOLwKWAEcCK4EHgBM217ft9PNw8k5JkgaGBQxJkjSmVdUSWkWK9n2L2j4XcHKnfduOval3UUqSpG3lIySSJEmSJGngWcCQJEmSJEkDzwKGJEmSJEkaeBYwJEmSJEnSwLOAIUmSJEmSBp4FDEmSJEmSNPAsYEiSJEmSpIFnAUOSJEmSJA08CxiSJEmSJGngWcCQJEmSJEkDzwKGJEmSJEkaeBYwJEmSJEnSwNuu3wFIkiRp07687779DuHX3rSw3xFIkiYwR2BIkiRJkqSB15cCRpLZSW5OsjLJgmGOJ8nC5viKJDOb/dsn+Y8kP0hyfZLT2/rsluQbSW5p3ncdzXuSJEmSJEkjZ9QLGEkmAWcDc4D9gGOT7Dek2RxgevOaD5zT7H8IOKyqngvMAGYnObQ5tgD4VlVNB77VbEuSJEmSpHGgHyMwDgZWVtWtVfUwcDEwd0ibucCF1XINsEuS3ZvtXzRtJjevautzQfP5AuA1I3oXkiRJkiRp1PSjgLEHcGfb9qpmX0dtkkxKshy4C/hGVX2vafObVbUGoHl/2nAXTzI/ybIky9auXbvNNyNJkiRJkkZePwoYGWZfddqmqn5VVTOAacDBSZ7TzcWr6tyqmlVVs6ZOndpNV0mSJEmS1Cf9KGCsAvZs254GrO62TVX9HLgSmN3s+mmS3QGa97t6F7IkSZIkSeqnfhQwlgLTk+yTZAowD1g8pM1i4LhmNZJDgXVVtSbJ1CS7ACTZAXg5cFNbn+Obz8cDXxrpG5EkSZIkSaNju9G+YFWtT3IKcAUwCTi/qq5PcmJzfBGwBDgSWAk8AJzQdN8duKBZyeQJwKVV9ZXm2JnApUneDNwBvH607kmSJEmSJI2sUS9gAFTVElpFivZ9i9o+F3DyMP1WAAdu4px3Ay/rbaSSJEmSJGkQ9OMREkmSJEmSpK5YwJAkSZIkSQPPAoYkSZIkSRp4FjAkSZIkSdLA68sknpIkSZI02o744Ff7HcJjXPGXr+x3CNKY4ggMSZI0piWZneTmJCuTLBjmeJIsbI6vSDKzk75J/qQ5dn2Svx6Ne5EkSZvmCAxJkjRmJZkEnA0cDqwCliZZXFU3tDWbA0xvXocA5wCHbK5vkpcCc4EDquqhJE8bvbuSJEnDcQSGJEkayw4GVlbVrVX1MHAxrcJDu7nAhdVyDbBLkt230PetwJlV9RBAVd01GjcjSZI2zQKGJEkay/YA7mzbXtXs66TN5vo+E3hRku8l+U6Sg4a7eJL5SZYlWbZ27dptuA1JkrQlFjAkSdJYlmH2VYdtNtd3O2BX4FDg3cClSR7XvqrOrapZVTVr6tSpnUctSZK65hwYkiRpLFsF7Nm2PQ1Y3WGbKZvpuwr4fFUV8B9JHgWeCjjMQpKkPnEEhiRJGsuWAtOT7JNkCjAPWDykzWLguGY1kkOBdVW1Zgt9vwgcBpDkmbSKHT8b+duRJEmb4ggMSZI0ZlXV+iSnAFcAk4Dzq+r6JCc2xxcBS4AjgZXAA8AJm+vbnPp84PwkPwQeBo5vRmNIkqQ+sYAhSZLGtKpaQqtI0b5vUdvnAk7utG+z/2HgD3obqSRJ2hYWMCRJkiRJGiOO+OBX+x3CRlf85StH9XrOgSFJkiRJkgaeBQxJkiRJkjTwLGBIkiRJkqSBZwFDkiRJkiQNPCfxlCRpBA3SRFsw+pNtSZIk9YojMCRJkiRJ0sDbqgJGkt9IMqnXwUiSpInJ3EKSJG1JRwWMJE9I8vtJvprkLuAmYE2S65P8TZLpIxumJEkaT8wtJElStzodgfFtYF/gVODpVbVnVT0NeBFwDXBmkj8YoRglSdL4Y24hSZK60ukkni+vqkeG7qyqe4DLgcuTTO5pZJIkaTwzt5AkSV3paATGcAnG1rSRJEkCcwtJktS9LRYwkhye5LwkM5rt+SMfliRJGq/MLSRJ0tbo5BGSk4ATgL9IshswY2RDkiRJ45y5hSRJ6lonj5CsraqfV9WfAa8ADhrhmCRJ0vhmbiFJkrrWSQHjqxs+VNUC4MKRC0eSJE0A5haSJKlrWyxgVNWXNnxOsl9Vfbz9eJKXdHvRJLOT3JxkZZIFwxxPkoXN8RVJZjb790zy7SQ3NuvEv72tz2lJfpJkefM6stu4JEnSyBuJ3EKSJI1/Ha1C0ubSJP+nKTDskOTjwIe7OUGSScDZwBxgP+DYJPsNaTYHmN685gPnNPvXA39aVc8GDgVOHtL3o1U1o3kt6fLeJEnS6Nvm3EKSJE0M3RYwDgH2Aq4ClgKrgRd0eY6DgZVVdWtVPQxcDMwd0mYucGG1XAPskmT3qlpTVdcBVNX9wI3AHl1eX5IkDY5e5BaSJL+ONQkAACAASURBVGkC6LaA8QjwS2AHYHvgx1X1aJfn2AO4s217FY8vQmyxTZK9gQOB77XtPqV55OT8JLsOd/Ek85MsS7Js7dq1XYYuSZJ6rBe5hSRJmgC6LWAspZVkHAS8kNbjH5d1eY4Ms6+6aZPkScDlwDuq6r5m9znAvrSWYlsDfGS4i1fVuVU1q6pmTZ06tcvQJUlSj/Uit5AkSRPAdl22f3NVLWs+/zcwN8kbuzzHKmDPtu1ptIaLdtQmyWRaxYuLqurzGxpU1U83fE5yHvCVLuOSJEmjrxe5hSRJmgC6GoHRlmC07/uHLq+5FJieZJ8kU4B5wOIhbRYDxzUTeh0KrKuqNUkCfAa4sar+tr1Dkt3bNl8L/LDLuCRJ0ijrUW4hSZImgK5GYCTZHjiJ1hDPAr4LnFNVD3Z6jqpan+QU4ApgEnB+VV2f5MTm+CJgCXAksBJ4ADih6f4C4I3AfyZZ3ux7b7PiyF8nmdHEdRvwx93cmyRJGn29yC0kSdLE0O0jJBcC9wMb1ms/FvgH4PXdnKQpOCwZsm9R2+cCTh6m33cZfn4MqsrhppIkjT09yS0kSdL4120B41lV9dy27W8n+UEvA5IkSROKuYUkSepIt6uQfL+ZkwKAJIcA/97bkCRJ0gRibiFJkjrS7QiMQ2hNrnlHs70XcGOS/6T15McBPY1OkiSNd+YWkiSpI90WMGaPSBSSJGmi2ubcIsls4O9oTQ7+6ao6c8jxNMePpDU5+Juq6rrN9U1yGvBHwNrmNBsmDZckSX3SUQEjSarl9s216V1YkiRpPOtVbpFkEnA2cDiwCliaZHFV3dDWbA4wvXkdApwDHNJB349W1VlbcXuSJGkEdDoHxreT/EmSvdp3JpmS5LAkFwDH9z48SZI0TvUqtzgYWFlVt1bVw8DFwNwhbeYCFzYFk2uAXZLs3mFfSZI0IDotYMwGfgV8LsnqJDckuRW4hdZyZx+tqs+OUIySJGn86VVusQdwZ9v2qmZfJ2221PeUJCuSnJ9k1+EunmR+kmVJlq1du3a4JpIkqUc6eoSkqh4EPgl8Mslk4KnAL6vq5yMZnCRJGp96mFsM95hJddhmc33PAT7YbH8Q+Ajwh49rXHUucC7ArFmzhl5XkiT1ULeTeFJVjwBrRiAWSZI0AW1jbrEK2LNtexqwusM2UzbVt6p+umFnkvOAr2xlfJIkqUc6fYREkiRpEC0FpifZJ8kUYB6weEibxbSWak2SQ4F1VbVmc32bOTI2eC3ww5G+EUmStHldj8CQJEkaFFW1PskpwBW0lkI9v6quT3Jic3wRsITWEqoraS2jesLm+jan/uskM2g9QnIb8Mejd1eSJGk4XRUwmuXM3gD8r6r6QDNz+NOr6j9GJDpJkjSu9SK3qKoltIoU7fsWtX0u4ORO+zb739jp9SVJ0ujo9hGSTwK/Q2t2cID7aa2fLkmStDXMLSRJUke6fYTkkKqameT7AFV1b/PMqCRJ0tYwt5AkSR3pdgTGI0km0SwxlmQq8GjPo5IkSROFuYUkSepItwWMhcAXgKcl+RDwXeDDPY9KkiRNFOYWkiSpI109QlJVFyW5FngZEOA1VXXjiEQmSZLGPXMLSZLUqa5GYCS5APjvqjq7qj4B/HeS80cmNEmSNN6ZW0iSpE51+wjJAVX18w0bVXUvcGBvQ5IkSROIuYUkSepItwWMJyTZdcNGkt3ofiUTSZKkDcwtJElSR7pNED4CXJXkMlqzhR8NfKjnUUmSpInC3EKSJHWk20k8L0yyDDiM1kRbr6uqG0YkMkmSNO6ZW0iSpE51PUSzSSpMLCRJUk+YW0iSpE50VcBI8kTg94C92/tW1Qd6G5YkSZoIzC0kSVKnuh2B8SVgHXAt8FDvw5EkSROMuYUkSepItwWMaVU1e0QikSRJE5G5hSRJ6ki3y6helWT/EYlEkiRNROYWkiSpI92OwHgh8KYkP6Y1zDNAVdUBPY9MkiRNBOYWkiSpI90WMOaMSBSSJGmiMreQJEkd6aqAUVW3J9kVmA5s33bo9p5GJUmSJgRzC0mS1Kmu5sBI8hbgX4ErgNOb99O6vWiS2UluTrIyyYJhjifJwub4iiQzm/17Jvl2khuTXJ/k7W19dkvyjSS3NO+7dhuXJEkaXb3KLSRJ0vjX7SSebwcOAm6vqpcCBwJruzlBkknA2bSGjO4HHJtkvyHN5tD6JWY6MB84p9m/HvjTqno2cChwclvfBcC3qmo68K1mW5IkDbZtzi0kSdLE0G0B48GqehAgyROr6ibgWV2e42BgZVXdWlUPAxcDc4e0mQtcWC3XALsk2b2q1lTVdQBVdT9wI7BHW58Lms8XAK/pMi5JkjT6epFbSJKkCaDbSTxXJdkF+CLwjST3Aqu7PMcewJ3t5wQO6aDNHsCaDTuS7E3rV5rvNbt+s6rWAFTVmiRP6zIuSZI0+nqRW0iSpAmg20k8X9t8PC3Jt4EnA1/r8poZ7tTdtEnyJOBy4B1VdV9XF0/m03oshb322qubrpIkqcd6lFtIkqQJoNtJPP+/DZ+r6jtVtRj4qy6vuQrYs217Go//pWWTbZJMplW8uKiqPt/W5qdJdm/a7A7cNdzFq+rcqppVVbOmTp3aZeiSJKmXepRbSJKkCaDbOTAOH2Zft+u3LwWmJ9knyRRgHrB4SJvFwHHNaiSHAuuax0ICfAa4sar+dpg+xzefjwe+1GVckiRp9PUit5AkSRNAR4+QJHkrcBKwb5IVG3YDOwH/3s0Fq2p9klNoLZM2CTi/qq5PcmJzfBGwBDgSWAk8AJzQdH8B8EbgP5Msb/a9t6qWAGcClyZ5M3AH8Ppu4pIkSaOnl7mFJEmaGDqdA+OfaD2P+mEeuzzp/VV1T7cXbQoOS4bsW9T2uYCTh+n3XYafH4Oquht4WbexSJKkvuhpbiFJksa/jh4hqap1VXUb8Hngnqq6ndZIiE8nOXAE45MkSeOQuYUkSepWt3Ng/GVV3Z/khcARwAXAoi30kSRJ2pRtzi2SzE5yc5KVSRYMczxJFjbHVySZ2UXfP0tSSZ66FfcmSZJ6qNsCxq+a91cC51TVl4ApvQ1JkiRNINuUWySZBJxNa+LP/YBjk+w3pNkcYHrzmg+c00nfJHvSmmT0ju5vS5Ik9Vq3BYyfJPkUcDSwJMkTt+IckiRJG2xrbnEwsLKqbq2qh4GLgblD2swFLqyWa4BdmiXXt9T3o8D/AWqr7kySJPVUt8WHo2mtHjK7qn4O7Aa8u+dRSZKkiWJbc4s9gDvbtlc1+zpps8m+SY4CflJVP9jcxZPMT7IsybK1a9d2EbYkSepWp6uQAFBVD9CabGvD9hpgTa+DkiRJE0MPcovhVicbOmJiU22G3Z9kR+DPgVds6eJVdS5wLsCsWbMcqSFJ0gjqaARGku827/cnuW/o+8iGKEmSxpse5hargD3btqcBqztss6n9+wL7AD9Icluz/7okT+8iLkmS1GMdjcCoqhc27zuNbDiSJGki6GFusRSYnmQf4CfAPOD3h7RZDJyS5GLgEGBdVa1Jsna4vlV1PfC0DZ2bIsasqvrZNsYqSZK2QUcFjCTv2tzxqvrb3oQjSZImgl7lFlW1PskptObRmAScX1XXJzmxOb4IWAIcCawEHgBO2FzfrbwlSZI0wjqdA2PDryPPAg6i9UsGwKuBf+11UJIkadzrWW5RVUtoFSna9y1q+1zAyZ32HabN3t3EI0mSRkanj5CcDpDk68DMqrq/2T4N+OcRi06SJI1L5haSJKlb3S6juhfwcNv2w8DePYtGkiRNNOYWkiSpI10towr8A/AfSb5Aa/mx1wIX9DwqSZI0UZhbSJKkjnRVwKiqDyX5GvCiZtcJVfX93oclSZImAnMLSZLUqW5HYFBV1wHXjUAskiRpAjK3kCRJneh2DgxJkiRJkqRRZwFDkiRJkiQNPAsYkiRJkiRp4HU1B0aSJwK/R2t5s419q+oDvQ1LkiRNBOYWkiSpU91O4vklYB1wLfBQ78ORJEkTjLmFJEnqSLcFjGlVNXtEIpEkSRORuYUkSepIt3NgXJVk/xGJRJIkTUTmFpIkqSPdjsB4IfCmJD+mNcwzQFXVAT2PTJIkTQTmFpIkqSPdFjDmjEgUkiRpojK3kCRJHemqgFFVtyfZFZgObN926PaeRiVJkiYEcwtJktSpbpdRfQvwdmAasBw4FLgaOKz3oUmSpPHO3EKSJHWq20k83w4cBNxeVS8FDgTW9jwqSZI0UZhbSJKkjnRbwHiwqh4ESPLEqroJeFbvw5IkSROEuYUkSepIt5N4rkqyC/BF4BtJ7gVW9z4sSZI0QZhbSJKkjnQ7iedrm4+nJfk28GTgX3oelSRJmhDMLSRJUqe6eoQkLX+Q5H1V9R1ak23NGJnQJEnSeGduIUmSOtXtHBifBH4HOLbZvh84u9uLJpmd5OYkK5MsGOZ4kixsjq9IMrPt2PlJ7krywyF9TkvykyTLm9eR3cYlSZJGXU9yC0mSNP51W8A4pKpOBh4EqKp7gSndnCDJJFqJyRxgP+DYJPsNaTaH1nrw04H5wDltxz4LzN7E6T9aVTOa15Ju4pIkSX2xzbmFJEmaGLotYDzSFCAKIMlU4NEuz3EwsLKqbq2qh4GLgblD2swFLqyWa4BdkuwOUFX/CtzT5TUlSdJg6kVuIUmSJoBuCxgLgS8Av5nkQ8B3gQ93eY49gDvbtlc1+7ptM5xTmkdOzk+y63ANksxPsizJsrVrXWZekqQ+60VuIUmSJoBuVyG5KMm1wMuaXXOb9dq7keFOvRVthjoH+GDT7oPAR4A/fNxJqs4FzgWYNWvWls4pSZJGUI9yC0mSNAF0VMBIsnjorub9iCRU1VFdXHMVsGfb9jQev957J20eo6p+2hbvecBXuohJkiSNoh7nFpIkaQLodATG79B6pONzwPcYfoREp5YC05PsA/wEmAf8/pA2i2k9DnIxcAiwrqrWbO6kSXZva/Na4Iebay9JkvqqZ7lFktnA3wGTgE9X1ZlDjqc5fiTwAPCmqrpuc32TfJDWnFyPAnc1fTb7Y4okSRpZnc6B8XTgvcBzaH3JHw78rKq+06zZ3rGqWg+cAlwB3AhcWlXXJzkxyYlNsyXArcBK4DzgpA39k3wOuBp4VpJVSd7cHPrrJP+ZZAXwUuCd3cQlSZJGVU9yi21Z3WwLff+mqg6oqhm0RnW+b6vuUpIk9UxHIzCq6lfAvwD/kuSJtNZqvzLJB6rq491etFnidMmQfYvaPhdw8ib6HruJ/W/sNg5JktQfPcwtNq5uBtCM3pwL3NDWZuPqZsA1STasbrb3pvpW1X1t/X+DLc/FJUmSRljHk3g2ycUraSUYe9OaNfzzIxOWJEka73qUWwy3ctkhHbTZY0t9m1VRjgPW0RrdKUmS+qjTSTwvoDXE82vA6VXl/BKSJGmr9TC32JbVzTbbt6r+HPjzJKfSevz1/Y+7eDKf1mMp7LXXXh2GLEmStkanIzDeCPwP8Ezgba25sIDWF39V1c4jEJskSRq/epVbbMvqZlM66AvwT8BXGaaA4fLskiSNnk7nwOh0sk9JkqQt6mFusdWrmyVZu6m+SaZX1S1N/6OAm3oUryRJ2kodz4EhSZI0aKpqfZINq5tNAs7fsLpZc3wRrYnDj6S1utkDwAmb69uc+swkz6K1jOrtwIlIkqS+soAhSZLGtG1c3exxfZv9v9fjMCVJ0jby0RBJkiRJkjTwLGBIkiRJkqSBZwFDkiRJkiQNPAsYkiRJkiRp4FnAkCRJkiRJA88ChiRJkiRJGngWMCRJkiRJ0sCzgCFJkiRJkgaeBQxJkiRJkjTwLGBIkiRJkqSBZwFDkiRJkiQNPAsYkiRJkiRp4FnAkCRJkiRJA88ChiRJkiRJGngWMCRJkiRJ0sCzgCFJkiRJkgaeBQxJkiRJkjTwLGBIkiRJkqSBZwFDkiRJkiQNPAsYkiRJkiRp4FnAkCRJkiRJA88ChiRJkiRJGngWMCRJkiRJ0sCzgCFJkiRJkgZeXwoYSWYnuTnJyiQLhjmeJAub4yuSzGw7dn6Su5L8cEif3ZJ8I8ktzfuuo3EvkiRJkiRp5I16ASPJJOBsYA6wH3Bskv2GNJsDTG9e84Fz2o59Fpg9zKkXAN+qqunAt5ptSZIkSZI0DvRjBMbBwMqqurWqHgYuBuYOaTMXuLBargF2SbI7QFX9K3DPMOedC1zQfL4AeM2IRC9JkiRJkkZdPwoYewB3tm2vavZ122ao36yqNQDN+9OGa5RkfpJlSZatXbu2q8AlSZIkSVJ/9KOAkWH21Va02SpVdW5VzaqqWVOnTu3FKSVJUh9t49xaw/ZN8jdJbmrafyHJLqN1P5IkaXj9KGCsAvZs254GrN6KNkP9dMNjJs37XdsYpyRJGnDbMrfWFvp+A3hOVR0A/Bdw6gjfiiRJ2oJ+FDCWAtOT7JNkCjAPWDykzWLguOYXk0OBdRseD9mMxcDxzefjgS/1MmhJkjSQtmVurU32raqvV9X6pv81tH5MkSRJfTTqBYwmGTgFuAK4Ebi0qq5PcmKSE5tmS4BbgZXAecBJG/on+RxwNfCsJKuSvLk5dCZweJJbgMObbUmSNL5ty9xanc659YfA14a7uHNrSZI0erbrx0WragmtIkX7vkVtnws4eRN9j93E/ruBl/UwTEmSNPi2ZW6tLfZN8ufAeuCi4S5eVecC5wLMmjWrJ/N1SZKk4fWlgCFJktQj2zK31pTN9U1yPPAq4GXNjyuSJKmP+jEHhiRJUq9sy9xam+ybZDbwHuCoqnpgtG5GkiRtmiMwJEnSmFVV65NsmFtrEnD+hrm1muOLaD22eiStubUeAE7YXN/m1J8Angh8IwnANVV1IpIkqW8sYEiSpDFtG+fWelzfZv//0+MwJUnSNvIREkmSJEmSNPAsYEiSJEmSpIFnAUOSJEmSJA08CxiSJEmSJGngWcCQJEmSJEkDzwKGJEmSJEkaeBYwJEmSJEnSwLOAIUmSJEmSBp4FDEmSJEmSNPAsYEiSJEmSpIFnAUOSJEmSJA08CxiSJEmSJGngWcCQJEmSJEkDzwKGJEmSJEkaeBYwJEmSJEnSwLOAIUmSJEmSBp4FDEmSJEmSNPAsYEiSJEmSpIFnAUOSJEmSJA08CxiSJEmSJGngWcCQJEmSJEkDzwKGJEmSJEkaeBYwJEmSJEnSwLOAIUmSJEmSBp4FDEmSJEmSNPAsYEiSJEmSpIHXlwJGktlJbk6yMsmCYY4nycLm+IokM7fUN8lpSX6SZHnzOnK07keSJEmSJI2sUS9gJJkEnA3MAfYDjk2y35Bmc4DpzWs+cE6HfT9aVTOa15KRvRNJkiRJkjRa+jEC42BgZVXdWlUPAxcDc4e0mQtcWC3XALsk2b3DvpIkaQIZoZGdr09yfZJHk8warXuRJEmb1o8Cxh7AnW3bq5p9nbTZUt9TmsTk/CS7DnfxJPOTLEuybO3atVt7D5IkaQCM4MjOHwKvA/51pO9BkiR1ph8FjAyzrzpss7m+58D/z969h1uS1fX9/6667Nu59Onu6RmGYQBBEBEUcJRrDIJA8AaKipefoJgQk3hJTKJEE5OQ+IRHk18UTfQ3GJWoCaCIkIgCcjFBSQIjRINIuAhzYWb6fs7Z19pVtX5/9MGnv/VZM9NMT/fZp/v9eh6eodasqr1qVe1aq6v3+ow90syeYGZ3mtm/SX14jPHmGONNMcabjh07dmEtBgAAq+qS/LIzxviRGONHL99pAACA+7IfLzBuN7Mbz9t+iJl95gLr3OO+Mca7Y4xNjLE1s9fYuUkJAAC4sl3KX3beJ37ZCQDA5bMfLzDeb2aPCiF8XgihZ2bfamZv6dR5i5m9ZG/N6lPMbDvGeOe97bv3Nymf9Q127qefAADgynapftl5QfhlJwAAl09xuT8wxliHEL7PzN5mZrmZ/VKM8cMhhO/d+/e/YGZvNbOvNrOPm9nUzL773vbdO/RPhhCeYOcmHp8ys795+c4KAADsk4v5ZWfvAvYFAAAr4rK/wDAz2/tPnL61U/YL5/3/aGZ/50L33Sv/zge4mQAAYPX95a8zzewOO/frzG/v1HmLnQv6fp2ZPdn2ftkZQjhxAfsCAIAVsS8vMAAAAB4Il+qXnSGEbzCznzWzY2b2OyGED8UYn3d5zw4AAJyPFxgAAOBAu0S/7HyTmb3pgW0pAAC4GPsR4gkAAAAAAPA54QUGAAAAAABYebzAAAAAAAAAK48XGAAAAAAAYOXxAgMAAAAAAKw8XmAAAAAAAICVxwsMAAAAAACw8niBAQAAAAAAVh4vMAAAAAAAwMrjBQYAAAAAAFh5vMAAAAAAAAArjxcYAAAAAABg5fECAwAAAAAArDxeYAAAAAAAgJXHCwwAAAAAALDyeIEBAAAAAABWHi8wAAAAAADAyuMFBgAAAAAAWHm8wAAAAAAAACuPFxgAAAAAAGDl8QIDAAAAAACsPF5gAAAAAACAlccLDAAAAAAAsPJ4gQEAAAAAAFYeLzAAAAAAAMDK4wUGAAAAAABYebzAAAAAAAAAK48XGAAAAAAAYOXtywuMEMJfCyF8NITw8RDCKxL/PoQQXr337/8khPCk+9o3hHAkhPCOEMLH9v55+HKdDwAA2D/MKwAAuDpc9hcYIYTczP6dmT3fzB5rZt8WQnhsp9rzzexRe/97uZn9/AXs+woze2eM8VFm9s69bQAAcAVjXgEAwNVjP36B8eVm9vEY4ydjjJWZvc7MXtCp8wIz+4/xnP9hZlshhOvvY98XmNlr9/7/a83shZf6RAAAwL5jXgEAwFWi2IfPvMHMbjtv+3Yze/IF1LnhPva9LsZ4p5lZjPHOEMK1qQ8PIbzczv3ti5nZOITw0ftzEgfGj3/tNWZ2cr+bYWYWfny/W4ADiXsYB90K3cNml+w+ftglOeqFYV5xOV0d9zOuZNzDuBKs0H18Ce/h5NxiP15ghERZvMA6F7LvvYox3mxmN38u+xxkIYQPxBhv2u92APcX9zAOOu7hS455xWXE/YyDjnsYV4Kr+T7ejyUkt5vZjedtP8TMPnOBde5t37v3fg5qe/88/gC2GQAArCbmFQAAXCX24wXG+83sUSGEzwsh9MzsW83sLZ06bzGzl+ylhj/FzLb3fsZ5b/u+xcxeuvf/X2pmb77UJwIAAPYd8woAAK4Sl30JSYyxDiF8n5m9zcxyM/ulGOOHQwjfu/fvf8HM3mpmX21mHzezqZl9973tu3foV5nZG0II32Nmt5rZN1/G01plV83PWnHF4h7GQcc9fAkxr7jsuJ9x0HEP40pw1d7HIcbPaaknAAAAAADAZbcfS0gAAAAAAAA+J7zAAAAAAAAAK48XGAdICKEJIXzovP+9Yq/8PSGED5xX76YQwnv2/v8ohPDrIYQ/DSH8nxDCe0MI63v/7iEhhDeHED4WQvhECOFn9kLMgMsmhDA+7/9/UQjhXSGE/7t3X/6TEELY+3ffFUI4EUL44N6/e1sI4Wn713LALITwYyGED4cQ/mTvufzkvfIihHAyhPCvOvXfE0K49bP39V7Zb5//PQAuF+YVuBIxr8BBx9zi3vEC42CZxRifcN7/XnXev7s2hPD8xD4/aGZ3xxgfH2N8nJl9j5kt927w3zKz344xPsrMHm1m62b2E5f6JICUEMLQzqX+vyrG+Ggz+xIze5qZ/e3zqr0+xvjEvXv2VWb2WyGEL7z8rQXMQghPNbOvNbMnxRi/2My+ysxu2/vXzzWzj5rZt5w/odhz1syevneMLTO7/vK0GBDMK3DFYl6Bg4i5xX3jBcaV46fM7B8nyq83szs+uxFj/GiMcWFmzzKzeYzxl/fKGzP7e2b2shDC6DK0F+j6djP7wxjj283MYoxTM/s+M3tFqnKM8d12LoH55ZethYB3vZmd3HumWozxZIzxM3v/7tvM7Gfs3H+94imd/V5n5/5znWZm32jn/tAHrBrmFTjomFfgIGJucR94gXGwDDs/9Xzxef/ufWa2CCF8ZWefXzKzHwkhvC+E8C9DCI/aK/8iM7vl/Ioxxh0794X4/Et1AsC9SN2TnzCz9RDC5j3s88dm9phL3TDgHrzdzG7c+2nyvw8h/FWzv/xbv2eb2X81s/9s5yYc53unmX1FCCG3c5ON11/GNgPnY16BKxnzChxEzC3uAy8wDpbuTz27N+a/tM7flsQYP2Rmj7Bzf5NyxMzev/fTuGBmqf+G7j2VA5favd1791Te/fkccNnEGMdm9qV27m/rTpjZ60MI32Xnfvr57r2/7XujmX3D3oTisxoze6+ZvdjMhjHGT13OdgPnYV6BKxnzChw4zC3uW7HfDcADJ8b4rhDCv7DOT4r2vgi/ZefW9bVm9tVm9r/N7EXn19t7G32jmX3i8rQYcD5sZl9xfkEI4RFmNo4x7upSPzMze6KZfeQytA1I2vuZ/HvM7D0hhD81s5ea2dLMnh5C+NRetaNm9pVm9vvn7fo6M3uTmf2zy9VW4HPFvAIHHPMKHEjMLe4dv8C48vyEmf3wZzdCCE8PIRze+/89M3usmX3azv3MaBRCeMnev8vN7N+Y2a/svdkDLrdfN7NnhBC+yuwvfyr3ajP7yVTlvZ/UvdzMXnPZWgicJ4TwBef9fN7M7Al27m9LnmFmD40xPjzG+HAz+zumP/X872b2r+zcz0CBVca8AgcV8wocOMwt7hu/wDhYhiGED523/XsxRhdEFGN8awjhxHlFjzSzn99Lqs3M7HfM7I0xxhhC+AYz+/chhH+y9+/eamY/emlPAUiLMc5CCC8ws58NIfw7M8vN7FfN7OfOq/biEMIzzGxkZn9hZi+KMfI3Jdgv63buft0ys9rMPm5mf2Rmo8+Gb+15s5n9ZAih/9mCGGM0s399ORsLJDCvwBWLeQUOKOYW9yGcO08AAAAAAIDVxRISAAAAAACw8niBAQAAAAAAVh4vMAAAAAAAwMrjlCP/zgAAIABJREFUBQYAAAAAAFh5vMAAAAAAAAArjxcYAAAAAABg5fECAwAAAAAArDxeYAAAAAAAgJXHCwwAAAAAALDyeIEBAAAAAABWHi8wAAAAAADAyiv2uwE4J4QQ97sNAHCQjQbDB+Q40/nsATnOZXAyxnhsvxuB1cXcAgAuzsZo3W2HEC5ovxj943d3On7A2nQpxRgv7AT3ES8wAABXhMc+4tFuO9qFjcHdWh/4sw89QC265D693w0AAOBK9uQv/FK3XZblBe23XC7d9u/f8gcPWJuudiwhAQAAAAAAK48XGAAAAAAAYOXxAgMAAAAAAKy80A0Ywf4gaAsA8Dm6JcZ40343AquLuQUA4HNxEEI8+QUGAAAAAABYebzAAAAAAAAAK48XGAAAAAAAYOUV+90AAFeHm778uVL2gf/19n1oydXjH77sn7ntZa7LGn/6Nf/0MrUGAIAH1t948fdL2Wte/7P70JKrx+/86z9y2+NlI3Ve/I/+yuVqDq5C/AIDAAAAAACsPF5gAAAAAACAlccLDAAAAAAAsPJ4gQEAAAAAAFYeIZ4A7tUXPvHZbrvp6XvPPPoAp4/8r3dJnVRg51Oe+Uy3nYUdqfNH7/7jC2kmUqrjfjub7087AAA4zzOf9YNuu81rqTPo/DHl7e/4GamTCux8yV//YbddhhNS5z+85pcvqJ1Q0+mfu+3FYrFPLcHVil9gAAAAAACAlccLDAAAAAAAsPJ4gQEAAAAAAFYeLzAAAAAAAMDKI8QTwL36yAffeZ91bnjG19yvY1dF5baL5VLqfMlfe5yU5WHmj1P3pU4vv9Zt//HvveeC2vR1L/urbvu//NIfXNB+++07vuPrpayf3eG2R/21y9UcAADu0XvepYGcXU960T+4X8eej32oZNXsSp3v+OYXSFnV7wSJHzkkdUZHH+22X/tPf/KC2vTTr/4Bt/13f+DVF7TffvuPf+d7paysP+K2j25sXK7mAGbGLzAAAAAAAMABwAsMAAAAAACw8niBAQAAAAAAVl6IMe53G2BmIQQuBA6ER3zVt0nZsjzrtodlJXV6zUzKysxnXhRz/bxZ1kpZ3fm6ZG0udUbtutveKIPUOToqpew33/g2bcQB9Te+4+vc9jXrh6XO8OxRKfvx1//bS9YmPKBuiTHetN+NwOpiboGD4vHf/kopK4ptt73en0qdwzMd2wfNxBcsdB4Rsk1txIafhIRyIVXW+z4XY1Brdte112k7f/yf/Zp+3gH1a3//ZW77yMYRqbNoHyxl3/DKH7pkbcIDJ8aoX6oVwy8wAAAAAADAyuMFBgAAAAAAWHm8wAAAAAAAACuPFxgAAAAAAGDlFfvdAACr7eEv+Ra33Z9pYNWa9d12XWhgZ0iEb7U7PlhrOl2XOvVEs4SGPX/82elG6kzy2m1XYVfqvPuPPyRlB9Vf/bbvlrJpz1+r02dukzqjXPsOAIBL6an/8KfcdlHrGJ03PmizseNSp966S8qqyu+3PHGtNqDWEM9+PfSft70jdZa1DxZto9b5uZ9/i37eAfXKv63Bm3Xp52CL5m6p0896l6xNAL/AAAAAAAAAK48XGAAAAAAAYOXxAgMAAAAAAKw8MjAA/KVHfud3SlneybNoe5plMcv9drEWpc4i25CyjbZ028O60kb1NHNjEHzZQ9qx1Lnrzk5ZrKXOlaQpj0rZNPeZF+uF9tNO8/FL1iYAAJ78939OytrwKbe9k2sGRm/T52uVh1qpM6mPSVl/dshtD6P+fW1c6udtbk7cdnHylNTZ/b++zvZsInWuJDONIdN5YNB54enFX1yiFgH8AgMAAAAAABwAvMAAAAAAAAArjxcYAAAAAABg5fEC4yKFELZCCL8ZQvjzEMJHQghPDSEcCSG8I4Twsb1/Ht7vdgIAgIOBuQUAAGmEeF68nzGz34sxflMIoWdmIzP7UTN7Z4zxVSGEV5jZK8zsR/azkTiYnvq8L5eyeeYDLD/4ux98wD6vqWdSljWdoKveXPcbrfl9ekOp08v6UmYTH/w0WNMwrDwOpGww67nt6zqfb2ZWX+sDK3fv0nNL+bxHPcxt/8XHPn1B++23fPuslFXD0277+Ez7qdjKpQzAvmNugUvma/76N0tZXPix9q2/+qsP2OflhYZhLio/ZmW5Bm3Hvv9jSta7QeqEcl3KqpkPCN/ofVLqZOWOlPWXfjy8PvHHpNvWR257u9YAy5Qvesxz3PaH//wdF7Tffmvu1LnTpF647eWgJ3UWaxq4CjxQ+AXGRQghbJrZV5jZfzAzizFWMcazZvYCM3vtXrXXmtkL96eFAADgIGFuAQDAPeMFxsV5hJmdMLNfDiF8MITwiyGENTO7LsZ4p5nZ3j+vTe0cQnh5COEDIYQPXL4mAwCAFcbcAgCAe8ALjItTmNmTzOznY4xPNLOJnftJ5wWJMd4cY7wpxnjTpWogAAA4UJhbAABwD3iBcXFuN7PbY4z/c2/7N+3cpOPuEML1ZmZ7/zy+T+0DAAAHC3MLAADuASGeFyHGeFcI4bYQwhfEGD9qZs82sz/b+99LzexVe/988z42EyvgqV/xxVKWbx3Ssk4Q0mx3IXUOVX77ed/yBKnztjd86HNs4TkPyjUIMlh029VuKXWWoXHb2TDowRvdry79IyhONGQymB6rOerDoc5MxlJn1PeBYL/3h3drmxIOSmhnVz3ZlbJxtem2s0L7qZhquOpzv/Gr3Hac6X3xmPWHuu2f/Y3fuqB2Arh3zC1woZ7//O+Rst6DdWrflj6gczFOBFjWfsx86Xdp0Odrf+U3PtcmmpnZjUP9vCr6MMzltoZ4Vz0/hvWLRuo0PS2rOmHnTalh4EVZSdmi58fIk1LDbLjt2/27v/u+RC11UEI7u6ZTDejcvfuo2y43Ndz9bK5zvh/9+3/dbS8W21LnyLXXu+0f+/FXX1A7cXXhBcbF+34z+/W9lPBPmtl327lftrwhhPA9ZnarmekoAAAAkMbcAgCABF5gXKQY44fMLLXO9NmXuy0AAODgY24BAEAaGRgAAAAAAGDl8QsM4DI4PtA1hIdN118OW79uMxat1Kmjz4To11HqPPebHi5lp3r+fWU2qqWOFbqOcbAcuu1lm8ipCH7daDZ+sNSJMZeyduHPd9lb18+Pui42Nn6/7VJzHPKpzw/5+mc9Vuq85V1/JmUH1cI0AyNrfb/EWtcB55muQa0bX2/dDkudIui1AgBcPtMNnSMMWx3b8zBz2+1wKXUGtT9WP9M/Irzsbz1Fys4M/H7tYZ3vTDc0z2trtOG2q0L/q8DzNZ9dNbAjUqeZ6+dltR8P++VRqXOo1fGwzXzZ6aHOLbItP+d6+Yt+UOrc/MafkbKDqurp3O1UZz7XbzWrbL47lbJ63e+3Obpe6lzX2/pcm4irEL/AAAAAAAAAK48XGAAAAAAAYOXxAgMAAAAAAKw8XmAAAAAAAICVF2LUAEBcfiEELsQV5NFP84GRcXModQZDDSpaG/kgpEGmAZab243bnpjWqXINCD277oOu5lFDtYqQCN8sfJDXZn6N1Am5D8yaZRp8VbZaNl+WbnujOit1NoenpWw49/tNxlLF8vq4294tUu9r9br80Zs+kKi3+h7/nOdLWYj+2vVqDdoK68elbDj0QVuDsCl1bix94NpW3Uidf/P6X0s3Fg+UW2KMqf/UJmBmzC2uNM9//nN8wZaOYWubIynLBz6csknMLbbmnf0SAc9NrfudyH2IZ3VEg77rqGNIf+iDNfvlQ6VO2/dzizposOh6o31Qz/zntTs6tzh6SEPLs9b3wc5M5x/Ztj9WPKltauqZlP3Gb94sZQfBt37d90lZUfg5wuGFzifDMQ0Wt8N+bvqgQxquemx0o9s+c4cG1f7wq/9Jsq14YMQYdbK4YvgFBgAAAAAAWHm8wAAAAAAAACuPFxgAAAAAAGDl6cItABdtM2647XnsSZ1qpGsk86x228NSv6LVxo4/dqXHXvTWpexPf/WDbvtxL/4iqTMol1LWRr/WcWa11CkKn3+wyHTZ9WCqayR7lV9vGpe639nE2tmF+TWRvah9OW6vc9sz0zWw1p6Rosc97wv9Z+1on3zsfR/XY+2zNmgGRWx825vpjtZJrJBfdq75eKghI1u5vy6DQtchAwAeOFul/3vHmOnYV1U6FhSFzyMYDDQnazbymRBVVUqdQV9zsv7rq1/rtr/2R75Z6uSZtimf+nEmHtGxtig7+Rqlzj8Wpu1cND5zK5SHpM6JSudJxWTqtnuTU1LHou+7/qZ+fm93KmUvfMY3+TZOtU9+94/fpJ+3zxZRz2/e+oiEptFr1z+l9+Zg6OcNpzY132LQmQc2dS51AH6BAQAAAAAAVh4vMAAAAAAAwMrjBQYAAAAAAFh5vMAAAAAAAAArjxBPYM8TH/14Kfvg//1Tt/30x3651GkGcykrRj4carauAZZBM4+sLX1Y0nalIZOx8l/bXqnvIUMimfHx3/KlbnuUaejiaENDpSa1P79FrqFWs5kPzGqDtns6e5CUrc1Pu+1Ya0DX0rRNTeMDstpa+2AZ/bH6PQ2LWlaJd7hnOuFma4elyuOe9xVS9n/e9t/0WB1Pe5YGp/7Ruz58n/t1PeVJz5ey6QmtF9f8TRYb7ct2oQFZ5YYPfKsyvVlPj32IV2g1JA0ArnY/9PSvlrL/9w/f6rZ/5LnfJHV2UmNf5oOvT/V2pc6k0XG0mXWe4bk+989O/Dhe5MekzuFCQ8Nf8IN/z22Xmf7RIt88qW3qBJBOo7Z73MmPzksd6BZndRzPl75sOd+SOnWrn2eln1usVTr+lz0/b7lma03qzHa1fzc2OvOyIxoi+q1f9Lek7HW/+vPazo6XfI3W+Y+/o8e6Ly967kulLJvouVRrfp60DBpaavlAig6ZD7gfz05LnTuau932cKDXAOAXGAAAAAAAYOXxAgMAAAAAAKw8XmAAAAAAAICVxwsMAAAAAACw8kJMhP3h8gshcCEusy8+9ky3vT7S4MnQ8yFLzcZI6+SJEK1OsFa7VkqdrNCArsx80NTGYCx16taHUcZuMJSZzUJfykL0bS8S+20NKymrGv+e89TZTamTdcJHm0r7xBoNglprfYjX4JSGhrXtRMrm5oNEU0GfefAhpQObSZ1me0fKitqHTC1K7aeiryFli9z3Ux6C1LnlnbdI2f3xxU94qpSFRCbzfOiD2waJnK3lQO/70dD353Ck57Js/bW6Zv3BUuf6RJjba97wOm0E7q9bYow37XcjsLqYW1x+3/8kH564mSfGnsKPF5OjGjLdJPabmq83LTWMe5FrsHjod8anQj/vTOvH6FDoWF/nG1K22fPjf/tQqWKHN7VNi85c5uT2DVKnO3NqM51bZKd0nOkHP4bF6Y1Sp0kET7ZL387N9qzU2ez5Y/cbnW9NbtNxtVn4r+Kyr8HXI9Pzazb9ddmuNTT0za/7CSm7P174tS+UsqyfmJc1fj45nGxLnd5Ig1uvO+TPZbSuc+plJ7QzDxoA3z+lff5jv/IjUob7J8aok74Vwy8wAAAAAADAyuMFBgAAAAAAWHm8wAAAAAAAACtPF44BV4mi79fi9ftbUmc+9Csws8FC6mz2tKzK/fKxRWIZcrnUTAgb+nWEy1liGVrp12hmlb6HLDNds7iWdz/vhNTJo7Yp7ywT3VomcjJKv44xJh4t06qWstj48+vH01KnyLV/s9yvQT0belJn1PrPa2Z6Ddq+ructOv07mA6kznKqa5Or3Leh3NBr9+SvfIKUzWufU2FLqWJ5p+lN1LW7VdC1pL3W3xv1UPuyV+ta0n7tz6+Y61rdctOvUy02Ev1UaJsA4Eq2OfQhEGtBn7s7dSfHqa/jcW+gY3uv78e+prlG6pSFZmc1wY9rbdAMjM3O/KZXam7FTtSsoyL386RQfErrjPTz2h2fQXXIEjkgw0f5fRodx5tW+6mZ+7nUerxd6oSejuNx6PuuMR374tKf73imc5v+hs6B6tgZVyu9Tr2omV+zXX/N1/qay/HCF2n+w9mzvl3rdpvUGWa+P+NS88SqMpHFVnQmKkcS+XCNZnUsoz9+HnTO2Yx839U9nVuszzSfBVcXfoEBAAAAAABWHi8wAAAAAADAyuMFBgAAAAAAWHm8wAAAAAAAACuPEE9ctfLw+W67zj8hdZq+D7XaSgR2LhcaupSv+5CnQrMirawSX7/WH38xTBy79aGLg0SA1aG5BoItwhG3PexpWmQ91yCmYuKDvOJcQ54WjQ/jKvoatLXsa9DV5LgvO5RrHStLKRq1PjCzaXW/2XLdbYdEWNT6UANJbeGDp6rRjlTZzdalbFDt+s9baPDULEsEki59n7cTDf9cFD7cLMy03f2jGri2DJ1j5XrsXiJcNeS+7cu+BmbFzJ9v3T8qdYpETi0AXMmasR+jm95xqVNu+sDqzfKk1JlrVrPlmZ9MtId0jtDmOj6UnQDJxVADLNdiZ96y1LFh0GjQdlP4oMtee0bbtJMIJN31Y0g51tDH7eZOX1DrfGC+1HlL7Ixr62vamXmmoaFZ7dtQL3SeNK59gHVba3/HXPu3XOsEfU+2tY5Or2wQfT+dmmm755XOuQYTfx9mS5231APfT8uRXvOtmQaZ7q517tdC+7csEyH0h329qtRzmQ0755II2G/La/XYuKrwCwwAAAAAALDyeIEBAAAAAABWHi8wAAAAAADAyuMFBgAAAAAAWHmEeOKq8JWf/ytStj35oNte1rtS54ZOllCsNYRxkkgq7NX+q9VOE+8K1/RYYekDubKggUp5z6c8VQsNWCoT4Vu99U7YV6VpUVkiDHPS96FZ7RkN8cqDP9+6OqJ1Bhpq2c0frXoafDU3DaeMnbCxZdBjW+dU2rEeZ7w5lLL1zAdyzaMeO9vWe6XoBLA2rQaLZZpFZTH6z2uGfamziL5s1NP7qd3RsDHr+3YOEglhodA+aEb+xo9LDRtrc/95O2Otc92GhoZ970tf7rZ/4bU3Sx0AOAh+6Mm/LmXz8YfddhY1VPuGpX/u70SdD5xJzC3yTr2qO9CZ2e5Igz3LcMptD0zHsPVwvf/8eWJus9Cwxix0QqZP61jUrOuYtegEZtr0Lj32WX8uVfkgPU6lY+aw9WPdfKZziyzTkMlg/pyLpY6PvU7Q9vKMHrtOBE/aug8Ez9Z0PK5Oad+VbXde2JM6vahhmDbwfR77GoBaZ8f8LqZ90i713twc+8/rV2P9/JFel3jI92c90ODW1rpzMA1Sn5d6b/7U9/2q2/6HP/ed2iZcMfgFBgAAAAAAWHm8wAAAAAAAACuPJSQXKYTwKTPbNbPGzOoY400hhCNm9noze7iZfcrMviXGqP9hbAAAgA7mFgAApIUY433Xwj3am2TcFGM8eV7ZT5rZ6Rjjq0IIrzCzwzHGH7mP43AhLqFveoyuU/3k+CNuuynfL3XWjvj1nqNEbsV4oOsYe53MgkGjaw+b3ikpqwv/ebOF/kgqH/g1ivVYP7+fJ95NHvFrK/Oo6xNtRz9vXPn1l8uBrpHMOpkU1VBzOYYhkb9Q+8+rZ7qOsih0bee49rkURa3nu+xkV/TO6Fphi7qO0oa+Ddds6nrT6qyuMa7mfh1sLK6ROsF0rWxb+myORatt6pYtlnoNyqneY71ONkl2TNeSjvqaU9EvPs9t531dl1tlnbKYWAecuA0fun6D216f6Rrj177xl3VHpNwSY7xpvxuBBx5zi4PhlV/5Zik7O/0fnRId6/vrfjxsNzVbalJUUhaGflybDA9LnfboHVKW1f7ZH5b63K+H/laZnb1Wjz3Vsac/OOuPncjSWlSanTHrZHxVPR3r2+jH32agGRhh90YpK5a+TdlCx6eQJ/KfzPdv3tOMhuG2rzOaaJ5IvdB5w+SIz8W6pqd92ZtrDsiy+Yzbrqba7mW7IWWxGHXq6Py1Hfh2BtN54TJoftgg+HnStZsfkzrrWzoH2nqwv37lhs4tzhb+2o0rnRONTz5Myh7aPNJ/VqnX4O/+xHdJGVSMUW/OFcMSkkvjBWb22r3//1oze+E+tgUAABx8zC0AAFc9XmBcvGhmbw8h3BJC+Gy8/nUxxjvNzPb+qa+xzSyE8PIQwgdCCB+4TG0FAACrj7kFAAAJZGBcvKfHGD8TQrjWzN4RQvjzC90xxnizmd1sxs88AQDAX2JuAQBAAr/AuEgxxs/s/fO4mb3JzL7czO4O4dx/UHvvn8f3r4UAAOAgYW4BAEAaIZ4XIYSwZmZZjHF37/+/w8xeaWbPNrNT5wVtHYkx/vB9HIsLcT898wm/57b7jQZBLk1DiCbjj3Xq3CZ12kN+fnhoU8Mbm3UNRpq1nfDESi9vr9B2zjshntOpBioNax96FOYajln2tKwa+JDHMEoEM2qTrO0eqtAwriz4vJ9JX9tdrg2krNn1++WaWWZhqeGUMfiyYS8RTmU+wClbapuqu/W+6A98J2wkvpnVQgsnncO3m4mws2xTysa7PrxtLRGuOl7rBLCO9dpljX7etN1y24OBhloNN/Q99kbmy/KhhobF3H9eFRKZT0vdb23h75/h8pjUyU3v33Fzwm2/44/erZ939SHE8wrE3GI1/Iu/+V/cdj3X524xToQXzz7htquF/odi6s7zMgw1YHm2pQPyuOfHh8Uo8dxd10DQEHzby9kRqRNrfy7TxDiTCnSO5sffZqIDeTXXdsbO0NPv6Xjclv72rdZ1HlGferTu1wn/tDrx97XNSIrazM8tyjzRprEv6y31x+z5jrYz71y7crQldUJ2UspmhR/7huNE+Ghfz2VS+HPumYar9o75Oe5m1PnWcpkId9/157IW7pY619yg98rDjj3YbfdHOtbPO3OwcdR7dTk7KmWjHR8a2ksEzg4y/U7NTvn7/u/9yj+XOlebgxDiyRKSi3Odmb0pnBuICjP7TzHG3wshvN/M3hBC+B4zu9XMvnkf2wgAAA4O5hYAANwDXmBchBjjJ83sSxLlp+zc35QAAABcMOYWAADcMzIwAAAAAADAyiMDY0WwTtXsWc/8Tbe93XxU6rRzXY/YmF8PWI3XpE4v0zWhm61ft1lHzbfIc7+u78h6K3VmR3S/nda3M7S6vrXt67HmnXeKk5lmFlxz0q9jvCbblTqncs1RONvvtCHXW67NdY1kN47gmsFhqXKy9X1QVfpuNNvUdZtht9MHja6dbWfazhj9sYpD+nmHS7+usUzke8xPaFl/zd8XRZtLnTqxOnC79ueSr+sP3LK+Hqsa+w6OeltYyHwf5DM9Tm56grWt+zo9bVOT6T12rHPR65C4niOf51HlD5I6/UbXwDa1Xw+et9qmzUzLquCvy0N72u5f/IO3S9kVjgwM3CvmFmY//f1vdtuL03dInSKxNr6bhxB3db1+mOogkjV+7Gmjjquz3GcWFKU+83avOSVlf7Hux8jlUMfMbKBlRe7PpdjWeZLt+j7oV2elyniYmIN15047ei5to2NIXvjBbpjIzqqH/tiNabZENb5R95v5YzW17recJMbo0l/jdk3nd23P3z/9szpoZ5/SvKtizec/HC0/X+rUg9S8zGeqHErkl20n/j560gkwC6XmtQwe5OePR9b0XPKTmq9RnfHn0s1BMzPbOqwZww857OcudSJjZNb6Plj0dW5R9DQDI5se8vvt6Pd1Y6nz5WbX3+d50GySH/zPr5SyK9lByMDgFxgAAAAAAGDl8QIDAAAAAACsPF5gAAAAAACAlccLDAAAAAAAsPL4z6jic/LCJ7zHbf/2h54pdZ7/tP8mZaNOEFMsNLjnrulH3PZifJvUCXPNlWljJ5wp18BOy85I0bT1YUVlrsFIIfPHvn2u7/xuDNdI2ebRTihQqYFOpxca7LnohDNmcSF1BnHqtmOrbdppNDQsBn++o0TQ5rJM5L1t+n7ZXmib6sr3UznWELHliUQY52F/7DyRjpnIc7Ri6cOoKr2dbLcTUnZdowFWw6F+3s5Zfw3aofZTzDUcqir8uQwSYZzFGW3oYuFPcFnotStq359lTwPY2iYRgBZ9m8Juon9H2i87A389Y9DAzGbSdurcqW3Kdb926dvUT4TXTRNJpv3Ml/UO8/4duFL845f/vtv+lzd/ldT5ib/9RimbHfbjoc312Xz89J+67cVd+qzKEmHNvZ4PQbalPr83cg26jLl/Xs+TocS+zm6rwYGlHZOyQ9f6EMLpmo7Hk6BzGeuMa5NG91vvhI3nfe3Ltp1KWc/85/WW+tzvB32m19m2346p8cIHM4ZxKXV6U51blJ1A8qzVca4Oel1C5/5JxJrb9sLP5waz26VO1tO+W1ZH/HbzGamTR507TTYmnToaLNpvExOO2l+rSa1tym3DbWdR5xHDxLw7X/jA13JjQ+oMNnXee6oz55pWGq7a7nSCcQsNyh1Gvcfj3F+tNmrQ507U73ms/f1z3brOzbF6mAECAAAAAICVxwsMAAAAAACw8niBAQAAAAAAVh4vMAAAAAAAwMoLMSZC+3DZhRBW7kJ845f8mpSNZz48qM5OSp35pobkbAxOue1FpoFKs5kPHJrPNcyoGmuAVJ35oKte1NilzO6SsmvLa/1xcg0qKlofejTf1PDEyTENE3p0JwRx+HCpYree0rCk0yd92OhgW8PGhgsfaFRPtU/GYVvKZqVvZ5xou/OhBmSF3Jc1QUOXlt3QzrP6brTtJfq3E+zVrmuIV50IwwqdQMdspp83GPjPO1ppyFXc1ntlXvtgr1lP7+e61Pt32fp7o+gngih3NThtOvLXocn0UdDLO22YrUudddOAuWbp+yAmAjNtqEFbTemD05paw+RC5x4L+d1Sp1fq96We+/t+UGn/rhepd+v+2XP9o/X7s9npu5t/492J41xRbokx3rTfjcDqWsW5xaue9/NS1lv6ca3a1PFisqXPip3Mh2hOFvpMbztD3eYJDelrpomQx8HQbevoaNZPBJL3C/+crQt97sbgz2VejqTOTl/Hw7DuZIJeAAAgAElEQVTmwxNHj9Q6O+FhUnbmjB+j69mfSZ1DSz+f69U6t5gHDUTvF77vhlM9l0HQ61k3vk1nat2v6QSbh3oodUIigLWtfJ+X/cQ41ybmisGfyzRRJ6z5z1uf61yqGus8qdzx41qbuKPakR5r3vNjXVHqPKJMtNOG/lzGiXF1a9MHXQ7H2k8bU/2+DBt/HxYjHevLIw+Ssuaov8Y7A71/q4n/Tsep3nPDRN+Vc9/2otFzqTMNQJ3u+mfPsWu1n9aXh932D/yHfyt1riQxRv3Crhh+gQEAAAAAAFYeLzAAAAAAAMDK4wUGAAAAAABYeYkF0bgaPP/Zv++2ewtdxzk9+edSNg9+rdi4rzkV/UzX8N2x28nAGCTW+ZvPOmgzzRCwgS7LGjUf8wWl5gOEKpGt0Pg1dL1GP29e+nYuQmI95EDfAza5/7zbP3RE6nxxpms5P3rE990dre6XNf7z6oHme4wXp6UsTvy5VFE/f7DQczkcN/znJ67veObPdxp0jWastWxR+fWQYUOvbyLSxEInpyI22u7YyZs4nXjcxZ7eF435D2wKvVfXlnqsynzby0TmRpVoZ9FZalgOEiutO1/PbJE4dq7tzAr/neoPtX+rxJrQrLNoPGv0IsSevy9q00yKWGl+Sb/bB7X25VkbS1m75tuen9FnVhj6Nrzsec+ROr/0tndIGYAHxg88/41ue3Ouz5fFmU9KWRl8veEhfVbNah176tZn8UxyfVbNOllW01KPvbWhc5m8kxsVFpq1sNCoA2s7WQ5FKn+pM7cIps/vZpA4l5P+uXfLB49Lncc/5CFSFh6+5Y/d03yCWefZvLBPa5tK7ScrOv2ypv07nybyJrZ9bsIyNWZ25m51rh0+6iey0SZ+XBkk2p2v6/0UW19vlOmYXUc/nw2l9mXY1DyPduSP1UzPSp1B4h6LnblpbLVNeab7FT3fB71N7bvh1I+Zw4XO09pKs+7aga83SNyrFnUunp3wc9piS3MqlkO/Xz1PjPVBnyt967RpfkrqnM4T89e8s99xzcMZdnL7fv1F3yF1vuONvy5luHT4BQYAAAAAAFh5vMAAAAAAAAArjxcYAAAAAABg5fECAwAAAAAArDxCPK8CT3vSn0hZdeZut10nwrGqiQbZzAsfDlm0GkI0rTSUZx58MFKz1JC+ee7DC7P5IamzNtaAwxB9wM9WT9s0yTS8aJ4fc9uLRFhTVvpjL3INSro2aljTszZ96NLyf7xf6vzK4Aul7MbKt/Oaa9akzvW5D+1am2m7Z62GSn269UFFn/60Bl8NK+2nI00naLOvIU/W90Fqs0LbbVMNzGxzH9aUT/XYRS8RTtkJ+wp9vS7zZSfAstZAyTboI7BY+gCyYU8DyXpzDQQbZv78QuI7ZZn2QTb19fqJgND5wNdp1/R8220NzKpbH3RVN4ngq5EGpw0bH+xV13qPzVvfzmGp7V40u9qm6J8rw6F+f7JG75+w49u+mOl9cbzydY4U2t8AHhgvedbbpCxO/ByhnuozZ5YIDV/0OmHUEx2fJqbP4nHnUXi21DGkyvzzpF7TeUSx1FDrrBP4d22mz7Omp/OdphskWuvzOu/MN1KB0luJse/L+r7ty9s/KnXefecZKXv42Qe77fVHXid1Dq/5ce1Idq3UaXZukLJbx7e67XqwLXUGQc9v0AlcH7V6r8w7w2i7ptdpETRAut/3c5lQ6PifddOxzazJ/DiW9xNjSO7b0B9ru0dndI7QdOYp+UYiIDwxb7ir8v10KnxGmzTS81trfNmRRC7+MPd1lokw7jbo3OLM1B9s0ep341itIZrDyt+H2dnEd+NBfo5QtzofKCaJ+VzTefa0ie95Ys631vr7oF7q9dy5y1+XQ4lj4/LiFxgAAAAAAGDl8QIDAAAAAACsPF5gAAAAAACAlccLDAAAAAAAsPII8TzgvvTxPgSoFzQ8qTf+Cykbmg+62k2E7ZzaqKWsbDoBQ5mGcfYLDfMpZjtuexESIYTRBwW1zXGpE0t959Z0gomWjYb79BPBXru5D/ssSg0vmo98O0eajWkvvUbPdzj21+Vu0/P9s0/+tJQtsqe47RcnghGXjQ+nWsSHSJ0H9zQArWx9sNduIghyMtFgouOd/sxbDZ4aV51gsVYDYEOiD4rC3z/5NBHWNNfkqWXnWG2tQVux9OGQi6D36vpM21l3wtWKqOc7WNN7JUz98ds2dR9KkTXdMKhEYGbWyQNrs0T42Fbicd65NfORhrSGpV6XMvf1yqjnUi5OdBqwI3XaQvvcCn++84HuF5rDUpZ3gtNmiWfIIvqOapfaT8967JOlLI423PYyEWz63v/2LikDrlTf/UwfZBdMn5XN9BNSNp/e4bYXibF3ekwDujd7nXDfROBxYTquHe4815czrbOIvs560Gd62degwrITap0V+lwYDPQZM2l929tF4tlc+udgcUjDGx+/rp/3iFt92R9EfQ7edfwPEh/3CLf9DRs3Sp2HHjritsfhkVLHTK/dIjvqtj8cPiR11oY6ZlonFL231DnneN3PaZueHmcQE+Pxpr8PwvCI1Okl+q6ofRtiqyGtVSe4tZclQksTY/Sy74+9luv42Ev001rtjxWj3iuLSWLe0PjxcNu0f0Ppx/Gy1Ovb5HrsQX292x4O9Hsegn5e3vN/1ljLE/tt+3YvEvdF1U3vNbO68x0O6xroXySued4J4t/u6XxntxPMn/r8H/srL9T9cv+suauvf4h4w9v+PynDfeMXGAAAAAAAYOXxAgMAAAAAAKw8XmAAAAAAAICVRwbGAZeVPt8iT6zXrxNr6updX29a6Vq8ejSSsqazDt0WifVytWZCtJ13ZfNGwwDWC7/ObT7UDAxrtU3Zwh97e6lrYDdDJWWjvj+XaqhrYB+25bdfdkjXpK7nD5ay+rTvz+3E2t3HP/YGKXvexnVu+8Rcr+di7PvlrkrXIe8EfTfZzzttb3WN8azVR0I17Ky3bHRtZ+yENISerg8MS11LWhS+3nIepE4dte+a6LMymoWeb77059I9DTOzYk33axu/32KpO1aZXpd8zdfLovZvVevaymknc6NMHHvWOZe+6XHCXNvZ9nwbau1KO1zod7HfWTsbCv28rHOJd+eJNfK59oEFfx+mMlWyVr+vE+usGc8SOSuFP9ayOCN1tvLEM2Tpz6+MmsUCXE3q2Z1uu9/T58Qy6hiW9/y683Gh3+Vx4rFQDfx4PNzRh9Va1GNV3ayDXmKM7uQMTXN9no1Mn0NN5Y+dZ3rsraBlm53MgHqQyPMY+c97+Lqu83/KQsfRk6c6czfTDIEbH3qNlH3VdQ9z22uZzlv+bOqv+V3FO6TOPJGHELPO3Gmp84jtkY7tG4d9WVzTOqP2WrfdlhtSJ2+1nwahM15Umv80T+Q2Nbt+v940kZkw8+NDv0rcO4l5qHVyKiaFnotFvZ+ywo91azs6hi10WmaT1h+/SOTK3NHJ5TiUOJey0D9D9NY7/TnUY2/VOidZ6/n7tUj8Nfqss9/pWq/T8TJxXQr/DFkL2qYQ9R47E/39etb0u5EP/LOn6d8qdTb6er7Nru+7QWIeivuHX2AAAAAAAICVxwsMAAAAAACw8niBAQAAAAAAVh4vMAAAAAAAwMojxPOAG9Tvd9t10LCdIhEas9v4kMlhIkAqVBMpW2Q+FKhpNBBnYVtS1o++XcOgQZ/z7G63nSUCD9cSgYrHO2FJ64kwrp2lvqtbLu9y28VM97ut83l/OL5W6myunZayuwrfL0945GGp89Lhl0vZp0/6gM56R8/3EzPf7rvLRPBm0ESnbsbjstHrO9KcKwvRB7XWqcdGp39bSwUsJe6Vyt8XbZYI1ewlAhV3fThTTJxvtvRlbaPXdz7WNjWZ/y70TMPrrJ/Yb+br9TPtp9hoqFTbCVyd9RKBmZm/MHPTAMu+drnlja/X6uPBmoWeSyy6gaR6XZpOyG9+SOtkiXYWmf+8NhGqtYj6PKoGPtysjhpe1w1Ta6OGDJ/NT2ibqlNuu3fmC6QOcDWZ5R9y2/VUHx5ZrQPG7a1/Diy39PldZGelbG3pv7tVpvvNaw0vjGNftuxr4F+V+znCLNNn+umz+vzaCJ1nceIBemqiQYxZ7ARILjUkPZT+GffhExpQ+qeJce1TIz93evDjjkmdL1t/jLZp6T/vtoWGpN829G1oE2NYlghEt9IfKyTG+o38EbpfftRt7iaCmdvO563leu3KxJyvmvg2hERwq61rcHuo/PF3d3UOVnbG1X6r41WZmBu3U3+sMnEf1qZzxfm404Yd/bw2Xi9ls+hDLTeG+nmZ+aDrWUz0UyJ0v1z4NkxKvQaTxPmt5f5c2kaDtuPQt2F4ZCh11hP34bQTbj6ZJgLR5zoHm1T+PugnrsGyE1S/neu9Om7ulrK27591663+eQH3D7/AAAAAAAAAK48XGAAAAAAAYOXxAuMBEELIQwgfDCH8173tIyGEd4QQPrb3T/09EgAAQALzCgAA0sjAeGD8oJl9xMw297ZfYWbvjDG+KoTwir3tH7kUH3x25NdzL3d0TV9cnJSy4cDX24yJ7IzdTSmzzprMpq+30CjXxfhZ5o8/CbomdD68wW1H0zbFxFrHteqTnQ/TNXzzRJ5G213emsjAsIlfd//ewWekSj/q+T5k068Nvu6xuv4zu13X+J6Y/Ynb/pNG11+eyfwa2GWbWGuZ+GbnuT+XXJcCWjlspaxpO/0y0/Ot884axaB1ZrUeu7tUNma6trNI3E/1IX+sbKkn3C79wRdLzUwoTNdWBln3rJ/frWJmls38d6rt6fmmrlXROeWF6TXPO+tSE6drdSLsJlgnp2IwlTqLxPrsovN9neb63SiiXxPazQ4xM+vl2tC8s3S1DnojhkqP1W87a9sTeT+x8P1UNbpOtZxdp/vVvlHLVvvpaU98hpQNjvjtd73zvVIHuJ/2bV5hZrZbdsa6sc4tBqbr0NueX/c+2NVnQNMknrujzlrxvuY/DEf6nR/3fdnJMpFlseHH2sFYszTmhbZzvfJt2lhq5sciMR4W1n1ea5vmnZyMceKZV3UnKWY2XPfn+12PeabUWR7Xa/Vb8ze57U8k5mCzTm5D3mi7e6bP67KTa3QkcV8UIx37Jn1/XeaNtrtZ+GNl80QGR6t5Xs3Cj/dNIr+kX+nnVZ3BdZkY10LWHfv083uZzp/Lue/PotJxdZnIaAhT384s1wy5oq/n1+b+OxSWeg0GnSyL2Zp+75Zrmlmz3fkunAg6ZvZbzYfpd+bwzVD3s+CvwSwm+iTxPe+1fg5Ul3q+y12dB+bmvwtF4u/228782RLfjazVPJp55zt85xnNnvnWJ3+dlK2v+WfkL77rDVLnascvMC5SCOEhZvY1ZvaL5xW/wMxeu/f/X2tmL7zc7QIAAAcP8woAAO4ZLzAu3k+b2Q+b2fmv566LMd5pZrb3T/1PVwAAACjmFQAA3ANeYFyEEMLXmtnxGOMt93P/l4cQPhBC+MAD3DQAAHDAXOy8Yu8YzC0AAFcsMjAuztPN7OtDCF9tZgMz2wwh/JqZ3R1CuD7GeGcI4Xoz00VPZhZjvNnMbjYzCyEkAhgAAMBV5KLmFWbMLQAAVzZeYFyEGOM/MrN/ZGYWQnimmf2DGOP/E0L4KTN7qZm9au+fb75UbWiDDxOqitulTj/TwJ+mEwo4b0upM+glAv8WPrBqqFlcFk0Df5aNb0NTJz6v9j8IaoKG7UwbLduoO6FApdapRxratdEJ6Jwt9QdJwXwI0BnNBLRRrv006eRMlbfp+Y6Pa+jSHdG389a1T0qdpvJf22Kh5xYT4Z+LzjS2Zxo+FnuJvpv788sSYWd56YOgYuL69hLhqvnU90G3v83MmkTZMJzy+/UTqZp93wf9sbZ7VuvcPnSCPYuQCNVsNcwtL/xFrxOhpT3TQLC1ThDuIhGYOal8H4SZnm9/Tfsp9HxZNtYwrKbVm3ps/pq3w0TIb+h8XqvfnzwmQksXvu+KQvvScr1/Yud7ViUCfeupv3Z5IugrRL0u3RzestKAu41uiJeZLc76a/eMv/JkqfPe//4/pQy4J6swrzAzmw98+Hd5SL8TqYDDYuG/u2GRCGac6jhTtT7QcdTTZ1xT61xmEO5w21utPq97Sz9GNomg5I2zGl54fe3nMm3U58kkEYDeLHy9Wodom3fGiyoxNkz0MWiDY/6ZMwqPkTqnGx3b78h84vBs7bTUqbuThEqDxueLxBjd6ZZhInQxZnr/hNLfY0e3NPhyVvn+jZU+h6eJ1PKsE77dnjqkdXb0ftow36bYS4z/fd93da7nNi70PxJUdII9E1MLaxOh4fnQX4el6Y2RBT2XY42/8WKt/TS1o257nJhKXX9E/1yRt/47vNjVefAynpGyu0rf9nJd75Wi59sZ6i2pszbXflp2Al/rmPgDSqnzjY2+/+7VeSL0f+7b3VaJ78FCnw/LygeZtqf1exf62nfjziPyG5/x9VLnt977Fim7mrCE5NJ4lZk9J4TwMTN7zt42AADA/cG8AgAA4xcYD5gY43vM7D17//+UmT17P9sDAAAOLuYVAAAofoEBAAAAAABWHi8wAAAAAADAymMJyQH34ff9kNt+0tN/VOpkmnlktvChNWGpQVtl0FCprNhx23Wt+1mp4XptdpvbrsptqZPP/e2YLTekznCpaVh16IQ81RqolO/oucyWfr82EZ64zHywV54I8drJNExopxMwePcJPXYRNGAoq3ygUd7qxZtnPvkpZpoE1SRCrZpOyOIykSDVVNpPzbwTwJo4dqh9oFFrE6lTJ8K/ymVnv0SYbDnT69n2fdkw07580KQT7npEw5u2d/TYVfABTrHQY5eJkLRFJ8ytDXrt1kzLhhv+OizrRPCU+fsiLLXdsdDgqRg716rVgN15rvfBYOD7Kks8Cyzz51smQmFDo/d9zDoBd0Hvp1TAbN0Jhi0SIall5z+4UM302HWmAYKh8kGm2VoiPW+sIX+z1vf5aFOD24CD6A/e9kq3/dyv/wmp0yZCHtu+f3aUOxrSt7ZIfHc7wYSF6X62ps+qZeYD8CZTPfa4EyY8yDS4eHNNgwLbaSe0fKHPTys0xHM3+nnRNNNnx6lOWZZVUicM9bkbRn58+Ezq7yFznSflTaedeSLguBN83TT6+fkyNf77fpolwl2zhQYVzuxO/3l2Suq0jb8v8pAY52b63M0mnXlhImS6F3UsKPNb3XbI9frmjQ+nHOU6Xu2WGmC5E30fxPlRqbOemONWhS+bN3o/xVbLBna92x6P7jtcdbO3I3VCIqi2zPyxZomA/SZ1rJ7vqyIRmBmbTohnqwGlWTd528wy68xbFjqXWmv1zyyz3Le9p4eWMNfpGZ3PVnPtg/7CX5eyp/dhU+t8ebfTLz27Rht1leMXGAAAAAAAYOXxAgMAAAAAAKw8XmAAAAAAAICVRwbGFWbZ02yJ4Tyxnnzk8yXmc11/2c503Vuzc4Pfr9G1pGtbuqatmXeyDhL5D4vKt3MzkaPQJtb+9zprxepGP79NrEHtD/x6+YXpuvd569udh8Qat6DrCsedNYOLTN8Vbma63zWLh7vtYfsJqbNT+nWUbaJNTWKNsUVflie+/YPEK81xZ81erHV9a79T1FR6fetG+7fqrOUsMr12lul9sL70bcpavQ8nPX991yd6DxwptJ9i678Loa8ddUKXZFrbWU+7lujfttR1k7uLzhrfnrbp4Z111nfXup54N5H3EPu+P4tCc2UGPV1XHga+De1S77HYuVZ5qxkjbeLzsujXzTdD/R6EuT7Hqk7mRa/U+6md+37JWr2fatPz7fZ4KnvmRKkXvTS/br4Ouk71OX/j+W77Ha/5XakDrLppvF3KevURKSvCg9z2zlBzcPJKB5qy8t/V9nhinDF9zreFP1ZV6/e7yv28oR81u2O31syC0cA/T8peYk7U02dx2x0zE8+T8cQ/49ZKnYNtDXS8GJl/pp4tdU50w+AGKbt+fsxt37b8tNQJncyC2NOxt5vTZWaWLbs5FXq+iYgmm879vbGca/8OOhlJ+ULHonGjY0/Tme+sDfX6Lg7pPW2tvy5FPKT7jX2OwqxKhCZolIUNev56ZomxvunrfbDsZG4tE1ls1iTyl1p/r/SifhebrDNv0K+YTU/qftVRf83DtXoNzLQsHO7M3Wodx5cz/x3Og96HtqFZFu2Ov1apfI3GTuuxOl2eiIez9e7zKTFXXRbaeb38vjPk5ks9Vrnmb6C4offht7/8x9z2f7pZc4quZPwCAwAAAAAArDxeYAAAAAAAgJXHCwwAAAAAALDyeIEBAAAAAABWHiGeB8jnPempUpav+xC5OL9N6kyXieCp9hFuOzuswUHNTIMY49AHa00bDZaZHdewu5D7enX/pNTpdUKXpolgxjbTNi077+GKRKhmm2kgmIRTNhqYFYIPL+oFDRzqZbpf1nb6qZ8IExpq2cNn17vtL6keKnVuCbe67dNBzy3m2k/rnZBJyzV4KtesJstz35/tRBOO8oU/VkgEog5b7ae80/Y6Eba6SIQlnc39tTtUaxhWVvl2TgaJ0MdW75W8c6/0dzUsKvT1vo+hEyC5pvdvngjfjLnvl16TCK9b89+XjUSI12RXA7OqzhN+1NdrHoIOA3knGC40+nn9TtDnIBFUl8o2q/JOQFfUwM4scR+UW51A3ZneF0XPn0uTCChdJgIEq87jL0uEauW5prLlnVuqbyekTq/elDJglTzz658rZTulv9+b3Vulzqmzx6Vs0ffPvbClz89eXwea7ZOd8WGs38Hh3fp9Lvp+fjPMNeSx6ARk1omwPav0+TnZ6oxrqXnE5o6UTTthifVcx76NTvjn0aEe+1BijNaMR51vbV6jIaWPm3yB275rfofUubUbzJwICM3LxHN37p+7/ahjyihoH1Stb2c90SDKzV0/PuSmfVImynZ7d/k6m9qmUOhYcLr0Y3s507EoK/21WrR6j8dE+HfW6bu61PswxMRcvHOoYqTziH6ufZfXnSDzRAi+LT/qNqcb2ifLRHDqJDvrtg+vaaDvIHFP50d824ux9sF65ss2LNG/3cHXzLbX/LWqEsGx8/5ZKSu686nT2u628n2ZaT65tXPdb9kJmJ3MtX9j1D+jbWb+WVf3dW4xL3UeejXhFxgAAAAAAGDl8QIDAAAAAACsPF5gAAAAAACAlccLDAAAAAAAsPII8TxA/uKP3ydlj3j2C31BIhBvPNHwl7j8z257fXqt1KnrRHhh7W+ZdqnBQXPTYK9s7oOuikpDPNvMh4bF7LTUSQbphWOdD9M6VmrAYKwf7LbX7JNSZ9wJtSpKDeQrKg152ux83E7/lNSZ24Ok7NAxf/wvO/50qbPVCfb89PDDUseqM1J0tPbtfEiroWUfTnTdBzqBhnVMhGFmnRNOBaJGDcyse75N3cOYmeWFBnQNOiGpMWgI0rLwZYcSQVvLVABq48Oopo22+1AnPNfMbBn9dyj2tJ8WE/1uNEv//awTQZtl3QlqSwRmZpmGfxVzX6/p6/mGRBhWHX1CVZkngnk7Ib/NIf38UUiEwLX+umS1Jm0tEiNT0Qk3nYWJVsr9+a4v9R5vgt4HTfDnF2bXSJ1QJu7pof+8aqZ9EHY05A9YJe95y9ul7Bnf7IM9l5kGB857Os4UnVDJtTs0aK5/TeKZM/Df1c/s6nOwWNwtZYeWfgwZFPqsyob+86Z5Ioh6oM/U3R0/FmRDDce0eWoa7ceMk4ln1fqW//vDaxOBnZvaTMu2fTvn65/QJm1+npR98XVf6rabuzRI9f0773fbd9vt+vlLnbsdXvhrcEMiZHo7EaJ9thOKnjf6/Jz3OuGJcz2OZRoy2euM/8tax76dRNl87u/Dtbne95vB1xmMdD6Q97RNjXU+LyZCHxPjTNv68Sir9fMK03DKXuHnnbNc5zI7rZ+HtrVe30Q2pYXO9Lwa6J8zFpsaJt8z/2eGJtf+Xe98P9dDIvy81vs33/LfxaXphHbc6n69zm03yRLz9cL3y2CgfxbYCPqsW85822PQNp0aJf5c02lUv9Lnmk3u0rKrCL/AAAAAAAAAK48XGAAAAAAAYOXxAgMAAAAAAKw8MjAOuE++87fd9o03fYHUqStd89UPfo3Z7JSu+Qqmaw23w/X+OOWnpU6mH2dN3ckoqDSzYNAJQMgGN0idOupa0vmis1Ys0/Xruen60u7a/7HputGi9Z9X1NpPa4n+3er7dfY7jWYfTItHSdmo0+dHH7whddamD3HbT54ekTpVeZuUHT50o9tuTNfzfn5nDayZ2WjXr4N9X1/XCk9b/y50mbgJYtAsgKKz5jVvE+tbE4+ptpPPUlliDWzn/JaJzI9RpusY69Yfq426JnUy0/uwu3q3SGR+xFbfGYfO+uhZYs1t2/r7t4mJ73SinVb5/myneu3m69rOdu6PXw+OSZ3Q+uyKbKKfPwta1uv5siZxffOltqnqZLiEeSLfotMvRdA6hxrtg6Lx/VtHzeWIpgvSq9qfy3yq9++DFolgF2DFvfc3fC7GF32ljlfT5qiUhc73uZ1oLlfV0+9S1vhn8WSgc4S40HG8mfrn/GarmVRrc/90HiTmNkWuz9T2jN9ve6Tf5TrT8XfamW/0E2PBsJPbEGrtk2Ki+2WdZ9NiqjkVnxp8UMoee/jZbvvJ/SdKnc8/6edO2yc0X6uYfErKHtXz1+6aSp+7//Oj/1vK3jj053x2XZ+7k86lanuJcS7X/IVl8PdBv9F7pw06ByobX1aExDUofBvadb2f+qZ9EFo/tldB2xQTf69clb5sPTHfyVrNSFj0PuO2txPZWaX5uVQ+1+9dsaOfV3bG2ta0n2atnl/Wmefn6zoHqwrfd8vZTOo0S+3zqpOX1pSJPLGpXvNs7rO6sqWGfmSdbL1pIuNsMNd2to3vp2VM5LUt9Bk573T5ck33G8409+Rqwi8wAAAAAADAyuMFBgAAADcYGPcAACAASURBVAAAWHm8wAAAAAAAACuPFxgAAAAAAGDlEeJ5hanqRPhMocEys6UPIcoLDU9qEqGHtrjLb9pIqmyGQ1J2NvMBQ22uoUCL6AOdhksN1SoKDYLMmm58orYpa7UPik6QWNZblzpLu9ttn6n0nd/GVL9GdSfEsj2j++30/1zKbhw9xm0PCg32aga+78bXaOjj4LQGmfb6PoB1NNTrezZ8mZQ9c+uw2/7C6pNS50+m/rp8fKxhUWcS4Wohu76znQiUDIkA1IE/1qFESFpWd8JdgwZRZlGDmMaZ789FofdhFhP3WN9/93bGd0md3vxuKVsrttx2MP0O7zQ+wKkou/e8mQW9V2adbqnrxDvrhZ5fve6fD/1lIoyz8/572STCODO95tu5vzfaVsOp8uz/b+/egy29zvrOP89723ufS1+kluTGsrFEhGMjQHY0jmNjm5GJsYdM7JBQMSmmPDOkXCk8A2QyNZj8Q5GqVFFTM1QmU55JKeBEGRgzHnCwxwSMbAHGFwSSkGXZspHxRb7Iakl9OZe993td80cfFXrWs6RuqbvPeXf391Ol6n6X1n73s9/r6vfs9Tu+bYiCtepEQNesiAIETyeC28R/3qqcmeVF7a9Ppxf+dVd1tu3I0cQx/k0bEHbbW29xfe764P2uDRgTTQQAZ7m/5uxmJ81yU/hz4tB85tryKNyvK/01tpz76/zWYK9Np7NEWHOUiHfUZzDLVYW/pmYLe/2oeh8A2Ez9yvSIDdfrcj++ms9tTd/u/TUnbPltl6sN9p6qf91OIly9y79oll9Q+PDRa4/ba2x7KHE/PuGDW48tvtcsH1r3449XVEdc22TLhoR+fNt/3rujMd/jUz8GHKb+OFyLxhLD4UQYeOlDENei++jm3N/7shDVkBjzNrkfW3SdPTYXwY85u+Dv/5vVl8xytfTHYbf096flEN0PE2Owa6OQyar153kcRCkicipKV92e+8+yvuPv41Vr2/rjfjsNa7bOkPkxbigSx/3CnsN1IkRUE79AoGvsPk8Ft+dRyG6euBZlrT/GNAoS18S/ayb+8BXpbeO08e83ae1x+KP/5O+5Ph/4N/8xsfLLA9/AAAAAAAAAo8cDDAAAAAAAMHo8wAAAAAAAAKNHBsZlx2c99Op38+KI7ReCn/eWt37eZhDblg9+rtjpzs9py6M55qWecn3aIprj3j/h+ky6RC5GZ+fCdcXjrk9Z+dflcsj2af38wFDYz9cm5vOenvj5vBu5/SzbvX9WuNacdG3fKh8wy7vLQ67PmcxmQpzZ8us5Vfh5ftfn9rNsrvu6v5nIJpkOdm7jB7rXuD6H121N16vfB1L7fd739jisE3ktRwc/QbDO7HE4axPHXDQHdif440mK613TYmLnTYbWb99ueca1zRtbQ7/m93kI/rPsNPaYGsTvu/hMDLVfT56Y75lH87qbxp+v09bPue2CndPbFn7dRWnPKS383PM2kcmji+gakjg3Nfj56DrYLIks+OvTUNvjJyv8vNyh89tgKO2xkmd+3fmu3+br0TZ4YeGvtfXEzllvar8tgbFbzxJjhNJf05+MzqVi6s+3VnZc21awGQnFeuJ1/lIhQe38+LUuka+1Ya+xXZ/Irdj15+5GlFnQTP29oO/8NWYtel2euF6fGWwNT6i/Xwzuyi9ydZSV0Sz8/WLttN8vj8kfmeWt3N//H20eMcsPbz3q+jy57fMXbhhsxtl1R31+yafUZ0J984XXmuWvhWtdn7Kz+25t6y9cn8XE3y+Gtej+lCXumYk8JFnaXIHUPTNM7X1tt/LnhrR++/bRsVrniXtmYgx/prHb/PDwAtcnn/ux6RBlO0zFf5Y+Ghc2gz83FolcrqG3r8sGvy23T/rz/JHMjpNObfnjPr/Gnoth6s+xUPpt1+V227WZHz9LIjOmzuzrduOMExEJwX6W0Cfyf3r/frvBnp91Ik9svuU/yyQ6Xg9V/nWb0Ti/SGR+Xc74BgYAAAAAABg9HmAAAAAAAIDR4wEGAAAAAAAYPR5gAAAAAACA0SPE8zLz2P3fcm0vueU615ZFOUxD5QNpup1EmOBOFFaUeASWF5u+bdsG/EzaY67PehSK0xaPuD79wtcZqqvMcpEI9+sSh3rXRgFDicCfEIWGqs/akdObp13b9yxskNkNud+Wu8Njru3uxVfM8uNT/1m0tuuW4MN9JPdtTxSfN8vtyZsSNR1xbRulDeT63t6HRXXBhiB9e9OHam1nPjwpzobK80RwbO+DiYbOtg3q328aBUGerPzOG9QHT5VR4OykSYTXhcTrMhv8lC19kJkWiRDPzB6bTe5DtLooU6qo/LbsEuG5WWtP9EwSB3AimDaLwlT7xgeSDXEAayIsWKrKNdXRsTkRv00K9cdBFmwNReG3wVLtcVHrubeJiMikse83BL9Nqqnf55nYY2MzSwSSbdhzP5TcdrF67v70V1zbD7/pRtcWh/LufIe/fvedHyMMUe7z0D/p+mxW3+HbjthzsGv8eRrUXourwd8vuqW/Dj1ZRmHciTDjVvz1+khtrwOHlv4eMult23zdByU+nrh+Zo3dLrMwd31OnfHXxk9EY7cnE6HPbWu3S0iEN/a5H+88cPRh25Bf7fqUh1/i2g7NbPjnyxMBkrvRpnv0iA+w1NZvg2Fit4FOfDjm7iIRID3YN5wGf73Oojq7pd8m2Wk/5tPqGrNcSCKMs0yEaE6ic6HyY8cw+GMsLOw+vjYRZNotbbi5HvbnxhPi26S1x2abGLtNpn6cH6J75ta2v9eeGex9e7bh91192B+bYcOed1kiPD/L/D6vo9DyeuY/b9/b46Ar/flzsvfjq53WHgc7iWtPn6WCxe022Fj3+3d93V7XmiIRSnsZ4xsYAAAAAABg9HiAAQAAAAAARo8HGAAAAAAAYPR4gHEBVHWqqn+qqp9R1c+p6i/utV+lqneq6sN7fx496FoBAMD4MbYAAOCZkSZ2YWoRuS2EsKOqpYh8QlV/V0R+VEQ+FkL4JVV9t4i8W0R+7qCK/Or9PvDnhtfZkKW292F3OvjAn2wShd3ND7s+ee1DiDQKxVn0W66PdDaopxoSYTszH4wkmQ0BKjMfqLRc+HDKbGJDcYo4KVFENAokK1IBlod8wNA3J/bzPt77z3JN5z/LxmC37+HOv24RhSBNEkFb89aHf22XNqT0icWO67PWJj6fVNGyfz/ZsMGITSKIanfij6fjM9tWJ8Iiu0Rw2iTKPCpyH+IVBxpVjX//xfLrrq0YbPjXLBGOuT7491vO7X5ZVI+7Pk2e+HyV3S+l33VSDnZf1YNfT977/RIH05Vrvo8mwjezHRtGpYmArr6xx0q25o/V1A1mrbLBUxoSdSeCU8tpFP7ZJ4IAo4CssveBWbX6DTxEQaJF4to3TQSLzoLtFxI/E9jctsdP6BPXMMBaibHFR37/y67tdW96tVnua39/7Nf9eTJs2PO7TYQnTjN/X6miC+bitH+/7ox9v+3Kr3s+9ed8G4VMavD3gu0tf13Y6ex1Zy1L3MPyKCB86a+Di4m/Hz9c2mvjFzIfFj31l0+ZRuHm2/NEiGdmP1+ZuoI3/nVNFMxcTRN99AnXdjraTjP1YYZZbdsWTeq676+pm8Ful2J9w/XpJr5tWZ8wy/Pe38f7XTveSd2LEoe9TFv7frnfvTLMEwHSC9tWFH6M27Y+QPL08lq7ntwXpZNonJQIKD/UJUK1d6M6j6X+veD357Kz224ncYidGOxnmWV+H6xnvqYiGl8tE/f/du7P4ZDZ46dSHz5az6Lw7y6xgxPnyxCNk4qp37/a+nN/Knbb9cFv36qx22mrTwT6X8b4BsYFCGc99S+Lcu+/ICJvFZE79trvEJG3HUB5AABgxTC2AADgmfEA4wKpaq6q94vICRG5M4Rwt4hcF0J4VERk789rn+G171TVe1T1nv2rGAAAjBljCwAA0niAcYFCCH0I4RYRuV5EXqWqNz+H194eQrg1hHDrpasQAACsEsYWAACkkYFxkYQQTqvqH4rIm0XkMVU9HkJ4VFWPy9mfoIzKsGPnubUhMUcz8xMpu5ntV5fbvk9Ru7aitXPRisLP+Qon7fvNBz/HTVo/py1rr7LryRNz3BLrqqP5rZKYc1sWdk5ZVvr5gXrGz6N8orOv20jM7Wwmibn/wX6+k363SCb2s2St/7xFIg9h83G7Xw73ie1b+nl28XPOZeXnWtaVnQ/YJuZMytLP/dsJNv8hv8rP581P+88yyWydTfC5BhplXiwS84Ib9cfqvLcbvel9nydy33ZVNA83JDJOTg7rrq2MYmTKM36n70Y5Fe2u/zBV7rdTEZ2LXZM4Dgd/3Gtm92cd/D7XKEsinmMtIlKIP19ClC+RJc6NMrWdWruuOO9CRKRSe6wW6rNCMvHvN8vtNtjI/P7dyfz2zaJzQzufK7MQexysJea7As9k1cYWbXQuyfQq1yef+Pnjx9bseaI7/lxqgs8e2o3ylobgs6ymc3vOP75M3AwK//O8Prq3hqm/VvaJzIA6Gktcm/mssCq3dWaJbKk6+PvhMrruhkS+xm4iq6MP9v47V58nMmnibeCv+0Xw189ibl/XJ8YtQyI/LO+iDCH1+6WOVpYncsGazu+7ejvKWiqu8UWpH7tNZlG2U2LMVwx2G+wuE2Mw8ff/dtce931zyvUJ4o/fahEdh90h12fR+RyOJrP30Sd3fOZvG2WvbDZ+LLWZ+W2QV/ZcDMGf04vO7/NTtd0vben3+TIab2wPiVwQ8XWGwR4H7ZA4LpaJMa7aY0w7X9MQjUnmeeK4KH2d62t2XBQS22RI5IBVC1tTIkJOumhdh+sr6zsJV9anvchU9RpVPbL395mI/JCIfEFEPiQi79jr9g4R+eDBVAgAAFYJYwsAAJ4Z38C4MMdF5A5VzeXsw6D3hxA+rKqfFpH3q+pPisgjIvJjB1kkAABYGYwtAAB4BjzAuAAhhAdE5BWJ9idF5I37XxEAAFhljC0AAHhmTCEBAAAAAACjxzcwrlD5xAbLbIhPXao3fQiRRgFd+Y5PlmmDb8uj4KcikUjTTWyf9qTrIpUvSfL5k/Z1vQ846hLhlP3cBvV064lQwMKGPPWJsNPDnQ9w6qJwqu3KB1+dHPy6sih08GTjw32ORG9XJIK2lpn/vEuJQq0S+6Dy+UJSV7aGSe/DmnZrG7hWJgKWrip8nWtbNrCq3vE7OFcf6Nhl9uDogg9rbDO7zReJINchETw5LO1nWYg/Loo4qE5EvrVpnwc3iRDPjc7v80NLGwa10yWC4mp7HNaJMLss9+83i4InU0FbefCBlXkWB/ElQndbu+2Wgw++2u39Z8mX9nVV9QLXZ1Imgnjnts4+EXrXqv18m5IIO/W5eKI7dt91iRBPEX+sHG3t/lwUvu78aruvfuXOP02sG7g8zDbtNX0tu8716cvjrq1et+fu3J9Kspj7a7GIDePOjvhrbHPUnrvdwp/L+anEGOhxe/3Il/66q5V/XTzeOFElgkUL+7p28IGSWeLjllGQehb8zyG7RCB5O9gLX2j8/SLr7UYv4vuAiLSaCC2N7pGauK+WqZDnwV7Dp4lg8SwK0a4Wvqai9WO+aWuPi3DG151NfPBltmH3uWaJAzEK+kzkfMqk9p+lOWWPg6ZLBKmu+89Sr9kA0m7L38SG0h8s09Yeh13jA3XPRAGknfpjdVn58dXRqa09T4wL6/kZ17bW2teFzH/ebNeO3U43/hybJ8bUm9G+qrLE/vVlSlhG4aqJYNwuCknvJBUGmjjG47F4YtyttT+AjkQfb9r7scy0su/3Lz7wEV/TZYxvYAAAAAAAgNHjAQYAAAAAABg9HmAAAAAAAIDR4wEGAAAAAAAYPUI8r1CTKKSm3/DPsrJlIuAwyrbRQz7sZpYIgmyCDRg6nQikmUTBovlxH1QUtv37tVEAadZsuz7a+88XCttW5D6oqFYbAjTpfLhf3Z5ybUNpE3iqVBBU4vlhnAFUNL6mro7WteHDjELrAw4XUehRm3j/3O8WKdso4HAjESYUhXbmEx/iqYmg2Ka0/SbbPvQpqE9zzTZsAFmeOJ76qW0bEqFaXePbrh7s5z2j/jhsZ37bLQd7OdXc75ftpa/zBb0NDd1JbKdhsJ+3zPxxkSfCc+dRvyoR+Jr3iSCxaF194vjNo/NHhsTxnAg3iwNQM/XbRAd/ng2FrbMJPmws6+w2X2aJgDvxnyWb2tcNrS88FL6mtrf7Kg5NExEJidA54HJV2LxB0cwHX5b+kiMSbIilbibGJFPfthvsvS4VZjhEoZpa++ugVP6cL5voWtH6dWeJ63xe2fVrIjl4XtrwwnwnkdhZPOGa2iisuagSP4f0t3/Jo3uPLn04dtPadYdZYkXTRHhhbq+Ntd8kEpZ+XUMUpD5PjD/W1I4RqkRoejn4++FE7fasln5c2CXGFsMJe70eCr9fTk9s2zRxXPZzfx+fZPYYa+Ww67Mrfr/0bTSe7BP7IBH4vh7sibZMBH3P+2NmuQ7+/dc7f6+dV/b4PZIIO5/UiXFgsP9A6NSfi3VurwVZ4T/bmibOuyjQti/8Psgn/nV99O+D5eD/EdP0p81y1/jrxZAYu7XROKnL/bWgKPz+jD9z1vvtlNX+/a4kfAMDAAAAAACMHg8wAAAAAADA6PEAAwAAAAAAjB6Tc69QD33Kzv17xRuuc33yxNExn9h5X4vKz3trEvMBq0ftXK0i83Pjmk07V3Yo/bw3VT//sV7YfpPK5xpo5+eKVSFaV2IufsjtXLi29X2K0rdlwc5XaxJzD7ezRE3R/Lyru0QOiEbz/Do/p2/SXevaBrXbN+uedH26xPzhMth9XiXyJqS1B8siEaaR+10us137Werczyvs80S2wm6cN+HX3UfzETdyfzzVifmQEh0XVemPJ5HE6wbblsqkyIKfD/743G6Y7czXGaJ5wKmKmj6RU5HbbZCY0SyTzs95XRZ2P6TmAQ9x5se6X3ve++1ULu1n6Yod12cIqR1qa8gT588ymufdlv79q0ni+C1s2+7Cb8trGj8vNpRH7Wq6xBz93cR8e+Ayded7HjDL/+U/vcX1GY76c0mjcz4k8qayw75NonNVg88VyNbtupeTxLX5ycRYZt2+37RM5Ebphms7VNh8gDhvS0TE3bYT2VLrbSKzICqzVT9Q6zb964Y1u50OZ5uuj8zjbenvxxr8Pohv91m3cH2aRC5XGWUN9ZnfBkM0dsoSN78skb8Qgt3Adee30ySxffOl3VdN7cdgsw1bhKrPiBgGf68f1mwN0yFRd+J4mkQDnCKRQ9LNfbDMdGHvtcvaf96wbttm06Ouj3a+pqax48etxHExbf2/KybRGLMa/PYtezsmmVb+s2liPNcGu1+GxPYtE/foNsroWybG1FtP2s/XZP7cyGZ+LLWIxoU7iXVvNol/QyztNUSnfh/0iTyYKwnfwAAAAAAAAKPHAwwAAAAAADB6PMAAAAAAAACjxwMMAAAAAAAweoR4QkREyiOJwMxEAF89sWFC2bYPn1mmkgJzG66TJ8Ka8l277mHmQ/rqmQ/OyTdtnWHb190lEiSr3AbgFJ0PE+wWts4w+HXPE4Gka52ts07EJ25UPvCnUtuvEr/uQ50NKTuTCoLq/HaSYJ9XPvrVr/o+z9sps3Tja16U6OP3ZxNtlyERjNQmtm98+FSJZ7HDrt1XxYZ//6uzRBhWFBrWJkJLNxofnrSIcsv6OlF3Itz0VNSta32oVRYdB6Ukzh9NBOMto2ND/fYtNRFMFzX1g/+87TIKJOv8550dSYSktbbObOn3S1clAlej21UYEqFWpV33kKVC0nzAXL+Iwk4TgXO7la/zZGZrGib+dUMiXBW4UkzX/DnYTmeurWvsGKHZ8qGAu5kPRpxt2rC7bNenPLZR+He3ngiLfoE/v9ejIPGNwr8uJPKctbbXmKzz14UsCgRdTP1w/NRp31Yt7RtWE3/d3xR/zSmjMVC83URENuKaTvt1160PJW6i+8Xnf+ePXJ+L5XU/cptrKxPjsqKPjqdEYOdudsSvKzrspgt//B6VKNx9wx9zWSLEu44GLmHi75n5jj8Oh6UdE4Sl37+aCJiN90tY9++3HoVDrs/9MV7N/Xm3G91r8yLxb4jC1zkEe+4PSz8OliG6PiTCZMvcb6d8YY+D1NimGxLHdDQu6ju/P/vB1tmJH0e0uyddW9NEY7fWb5ON3n++Jtp27bo/Dov8yh5b8A0MAAAAAAAwejzAAAAAAAAAo8cDDAAAAAAAMHpkYEBERIrez9frl37OoCzs3LCu9odQmZo/FmVglDP/fu1wtVnOSt8nNfdfDtn5Y0ViPuR08DUtQjSnfnLY9ZF22yzOCj93dzH4bbBz0tautd+WQ+fn4sURCeuJOb8bE7ud7v7Lj7o+B+3Ln/r683rdy258oWur2lRGQpRZEE/2FJHDmd3m09736YdE3kS7a5Z3W3887RS+Ld+0c4OziT8ulon5rZ3YY6X0XSSL+kgik6JMzEEdokwVzRLzZEt/boQ+mnNb+/OujM6fPDHPW7YTGRTxQZ75bTkMfp61ROe1ip8jL5Wte0hsk/nczwMOardBMfH5GsPMn/ut2GNse+Zfd9enP+XrBK4QVUhkzuz4c7DYsv2qPtEnkSfTT7Zsw1GfZTVE9+3JkLjXJ+aY50ub7VQm7jNZ5a/zdXT9GBb+oh5mZ8zyWunX0yXmxrc79rrXnvHrzlLXzx27XY7N/Pjj6g17Pfv1e9/v13PA/vh37nper/vR297m2vKlv/cNu/YeUibufX30ulniZ8HT3I/5FvHNvfb3iza72rUVUZZDnfv72uNh07Wd7u0xNcn9PWwWbLZDNj/t33/wx+FaH+XhrfuxVJc4DIPa83oYfBZL1thtN/HDYGkPJcZXg91O00ROhZaJa8h5ZGBoVMSg/t8n8y7x75PM7qv1IZHdkcgKCVHT9hH/eX/5//uIf78rCN/AAAAAAAAAo8cDDAAAAAAAMHo8wAAAAAAAAKPHAwwAAAAAADB6hHjirM4/y5pkPshmbWn7deKDe+a5D6kZ1s5EDf7QyyY2nKoLqXAq/7oht0FI/ZFTrk8cACgikkeBP8XSBxVNC1tDGXzwVRmFIImItEdsMNFjX/iC67PtWkSeiJa/lugj8rlk6+XgoS9/87z6vfI7v8ssF4lwtemaTUEagt+/s0Rg5hBs0FVIBFjd9+Uvurbveen32NfVPuTJpbSKyKy251CbCNoaogDSJj/j+iwTgaQahX/pxNfU1f48L8W2ZYtEQGe0PeeV3wch+JqG6LzLMx/G2Q4+UK8R+37TVIhnFNAZh/eJiEyCTwRrCps2dvUkEWCcCGXrNuxnueuTBHYCT1clwu6absu1HZ7bi/Fm6c+3ndYHFT5R2fP5UOfDP4cocFAaf30ZEoGObXQvaAu/7ir3SYXtxK5Ly0QIYmE/X+1z/GR6zLeVjb2mPvg7n/adzsMDz+tVq+sDd/32efX7qb/2Ltugfp8v1Y5DNRGkOkuMOdeje1jR+fvjr/35v3Ztb7n5fzDL08TxGwd9iog8ObHHb1f4sWqe2TFJv+E/b7c84tra0o67E1nnIr0/pyZxsHjiXl839pyqq+Ouz1ZizCdr9lrTdv68C5kfY25lNuC2zfx4pyvtdlm2ft/V6sckk8pugyOFf12ZJcZJlb0g/NsP3un6XOn4BgYAAAAAABg9HmAAAAAAAIDR4wEGAAAAAAAYPR5gAAAAAACA0SPEEyIiclh8cM9278OpqoUNK1qb+/Ci4ogP39pe2JCaLhFmWEQhQG3nA/i61iddFb0N4An5UdcnJIL7NAovaqqrXZ+hswFD2vpT5puf+Yprw6V139f+0iz/zRtf7Pose3scZpJImQq+TTdtENT64M8D+YZvyub22NjI/fPhofLHTxdsYNTjbeK4VxsEVeiG67PsE0G1Ys/PbulDaDPx50aIXlds+HCqcm7Pjb7315CdpQ+nWt+wwVqN+OtFW/v36/NF1JIIG4uSxPLWX5/KRNiZRmGj2dSHf7XHfCDZvR+7x7UB+Csbg79WNZ2/pu7W9nzuTyeuZ+v+HDy0sNePsJYIDqzsNb2p/XUpJEK8i+ia2uZ+/NElgj0zsf2K3q87zgzNE6HpD/zbP3FtuLT+jy+9xyy/62Xv8J3qKMw+ce/rM3+fKcWOVa/qTp9XTdPo3rde+FDNzUTOdt4fNsvf7g+5Pm0UJ5/7jyKZ+PBPiQL8886PbTSc8DVl9v2qwv9bIOttSHmz8GOiftff2/PSnq+JfFAZWr9f4s3Sl378Uc+jsNPEuLDa9m1x5WuHfd11Yue9/31/6Npg8Q0MAAAAAAAwejzAAAAAAAAAo8cDDAAAAAAAMHpkYFyhfuQN15nltvNzQvN54vDoo7lhEz93K+v8XM5yEmVgLPycUB12zfKRzM8PnAc/F6+M8gGm0Zw+EZGw8PNpd+PQjTX/ukHs5Li2iefhYwzu/vIjru1VL7vJLGd9Yq5l549DiTISPv4Xnz+vGg5PbL5EN/hsB0lkO4SJnahZHfHnT7Vra9K5f/bciX9dMdh+cSbG2U6+LWzadQ0T36eL5sWGvnV9JolzuAl2XWFIrHvw+2UW7OTcWSKPJo+uBZLYvaXuurastDU0pZ+3P0tNqAVg/PRPfadZXi42XR894+eKT3btWOJ0Iq9HOn/dK3J77m41foxQrds59UfkpOuzSGTj9IW9ph0dfEBA44cN0kzsuqpEjJJGF6d+mQgxwIF7z0N3uLZ//Jp3muXQ+OOiWJ5xbXEsxr/67O+dVw0vEHsfrUs/lqmXfoy72dp71u7aluvTxHkt2/4+12U++yXTKCer8mObbOLvtTK1NdSJDKxpZf+90LZ+3B3UjxuGzrYtZv6c6rpEbl9tkyr6xo+lpLbbZZIYD+SJTMAmuj7Nc5/nUQ+Jax3OiW9gAAAAAACA0eMBBgAAQIzWMQAAFqpJREFUAAAAGD0eYAAAAAAAgNHjAcYFUNUXqeofqOpDqvo5Vf2ZvfarVPVOVX1478+jB10rAAAYP8YWAAA8M0I8L0wnIv8shHCfqm6KyL2qeqeI/Nci8rEQwi+p6rtF5N0i8nMHWKdX2sCddseHAoXOB/dM1mxwTa+JUKtd/1wsiwKr1mf+/XRuA3ce+bNvuz6X1ol9fj9cSrnaY7xIBG19/Et/+bzW/bqbXuTays6eG0MihK4pfKpkHgVtloMPw4pPsyrzgV2aCKrdigIzpxPfJ6wlzv1N269fJoI+1QZrhcx/ti6RollGQabZxCfczdTXtN7atkJ8iNYQhW/1E78tq+DXLaW91g2dvzX+yR/c7V8HXBorO7ZYbtlzsC0SKZd1IpQvm5vlWeaDxXdSGdqH7bmqlb/GNZ29Dtx9x8OJFQHnaWrDKatE5uP//uAHn9eq3/Wa2/zbDVEY9jIRZrvpwzBn0Vi8zf25Ed1WpV3z4/6u9R9wq7bnXVn4+2ohvm0I9rzu4wJEpO9tAGouviZZ+LDetcHWGRKB3aH1Qaa6tNeerH/S9SlqO07JGj++y1sf3FptHjPLzdSPiT78/3zKteHc+AbGBQghPBpCuG/v79si8pCIvFBE3ioiT0UX3yEibzuYCgEAwCphbAEAwDPjGxgXiaq+REReISJ3i8h1IYRHRc4ORFT12md4zTtF5J2p/wcAAK5sjC0AALB4gHERqOqGiPyWiPxsCGFL1X9NKyWEcLuI3L63Dv+9IgAAcEVibAEAgMcUkgukqqWcHWD8egjhA3vNj6nq8b3/f1wIVwAAAOeJsQUAAGl8A+MC6Nkfh/yqiDwUQvjlp/2vD4nIO0Tkl/b+fH5pPpfQ73z09Dn7vOW2Q65tObdhfmu9D/cZgk/a2owCBj9z91+c8/2BC/Hpzz9ill//PX/toq1bJ61ru+vhr9v3++7j/oWt/wlqV9pQK239c+VZZX+IWi988FWWd65tVtn0Tx38+w/dzLdFIZpa+nX3URjneuZvJ2HwAZ1xCW3iNlSpDw0VecyuJ/MhXsVg687mPsRTJv51+VoUSDY9v590A5fCKo8tbv+1b5yzzz/6hze4Nt21QX3Thf/iSF0kQgh7e/267z/84TnfH7gQv3LXr5vld73qH1y0dWvjj/H/9Z7fNcs/9X0/5PqUeeI+vvldZvlwIr+6GexYpu78vXee+ObXtLDBl8vOB2Y2u37coIO9/2ZVIox7sOvK81RAuP8wc7VBplkifHTS+zFB1p20y/Xc9SmX9t81XfAhnjLx2668Jgr/XGdscbHwAOPCvFZE/isR+ayq3r/X9s/l7ODi/ar6kyLyiIj82AHVBwAAVgtjCwAAngEPMC5ACOETIvJMj9PeuJ+1AACA1cfYAgCAZ0YGBgAAAAAAGD0NgYDqMSApHLi83HaD/w2HYc3nL+wu7Nxv9VNJJYrAkCaRN6G5n5MZcrvu+cK/rsv9PNF6ZueXDmt+bufG0s6dHWr/PHwxHPPvV9jXremW6zPt/XaqCrsR8tZfMpvWrrsd/Pxa2fCvq663++reT3zKv26c7g0h3HrQRWC8GFsAl5f//uVvcG1V4h59pn2RWW57n0lV79r79iJsuD5D6ccNtdrXLYO/1y7FD2a6yn6xrJr4da8Xp8xy3vtcvVZ9nU1pM78mm35scWR4wrVlXTR2qhN5Hgub59FMfQZHc3ji2taut2OgD//qR1yfMQohjD6sg29gAAAAAACA0eMBBgAAAAAAGD0eYAAAAAAAgNHjAQYAAAAAABg9QjxHgqAtYJx+4OYXu7ZPPPjI81rXm17yAtc2t9lQEpojrk9W205l5Z8957PWtXVRYOWTuz7Ea0d9W5HNzfJa5t9vIjb8s9PS9WkGH9DV7trwq3XtXJ9JFMYlIlLndVSj6yJ9VOf2jl93cfywa/vcfff4la0GQjzxrBhbAOP0T773ta7t33z2k89rXT/9/T/s2k7XNlRyyK52fZrG3v9D6W+sWeEDQiUK1X6y9vf6rc4HhPf6LbO8kRh/bFR2XZNpYmzT+8vaVtR2tNxxfY72J13bUrbNchP8uGFe2vHNcuo/b3btUdf20d/8qGtbBYR4AgAAAAAAXAQ8wAAAAAAAAKPHAwwAAAAAADB6ZGCMBPNUgSvTK77T5mJUw7rrk4u9PBQ6uD5riSvIMLHzNneCf2a9JWuubaJ27ups6eeESrDzUrN1X3eWuL8My2W0Ht9H4z4istvZOafddMO/rrBzXps1n8sxvcbPA773k59ybSuCDAw8K8YWwJXpx1/2FrNci8+WKssoJ2Pm75lZ51+nau+/86W/zOy2jWvrmzNmeZL4MXpe2HVXhxPjHZm7tiLbNcszOeP6ZNunXNuW2DpPFH4sU0eZH9NjiQyM630Gxu/96p2ubRWQgQEAAAAAAHAR8AADAAAAAACMHg8wAAAAAADA6PEAAwAAAAAAjJ5PIQEA7Js//9q3zfLf+O6bXJ/NaW6Wh9qHai4WPuiqjJrype9zKPNBW/GdoUjcKapNG/6ZL3ddnyzkrq0p7fs1uc+Kyif+DY/O7fP2x7T2r+vs+1XH/Lrrdf86AAAuJ+976HfN8ttvfpPrczi3YZj1wv9ceyfz4ZRS2X7l4O/Z65kfp+SVDQ3Nxd+Pw9TWtNG0rs+s2nJtmtl1hW7h+vSz3rXlra29bvxnWWY23LS7xo9tygl5yfuJb2AAAAAAAIDR4wEGAAAAAAAYPR5gAAAAAACA0eMBBgAAAAAAGD1CPAHgAL3qZTeY5TLRJztlQ6yO5D5Aqgs+oDM7Y5cXM7/uISSeY2d2XVXl191VNiCz6P3tpO58qFU227Trrn2I1+A/nuwWNjRUOh/QuVibmuVZ5j/bpNn2KwcA4DLy91/3OrOsnQ/sznfs/XdDK9dnsvuka2uXNkSzH9ZcH9EN1zSLVp8nxiTtxNZUdEvXJ0x8sGep9vNlwQd2nk4Ei5+WdbM8rPm624n9fMGXJPnVJ3wjLhm+gQEAAAAAAEaPBxgAAAAAAGD0eIABAAAAAABGT0Pwc5Sx/1SVHQFA/vZLr3Nts2i65yT381Tb1j+P3hmitsT8z3I2cW19aS9Hqn4uaTa1uRjLRedrCr7OIbfzTaulf92Z0n+WXTvlVqZzP3l2echOTF0r/brvf/izrm2F3RtCuPWgi8B4MbYAICLyk699hWvb2LX32iJMXZ+2921Nby8r3cRnUmlxxLVlUcpXVvh7tAabp6WZz61q13Z822Brah/3dT8uPryijrIz2n7d9dmt7FimOuQzOP704x91basqhOB36MjwDQwAAAAAADB6PMAAAAAAAACjxwMMAAAAAAAwejzAAAAAAAAAo0eI50gQtAXgfP39l73YtTW1D6zaKm3QZhk2XZ9C566ti4K2Fok4p09+4UGz/KobbnZ9Qmhc2zLbsH182S7oS0Tkgc/f7TuCEE88K8YWAM7XL7zy9a6tCT7o+3QZhXEmMh/LRNtuboPEm2LN9bnjk+8zyz95699zfVRPurYtseGbu4MP+u6m/nUf+eQfuLYrHSGeAAAAAAAAFwEPMAAAAAAAwOjxAAMAAAAAAIweDzAAAAAAAMDoEeI5EgRtATgIt930Utd218NfvCjr/v4b/rpr+8xXvnBR1g0RIcQT58DYAsBB+G9f98Ou7b1//JGLsu6fuOWHXNuv3f/Ri7JuEOIJAAAAAABwUfAAAwAAAAAAjB4PMC6Qqr5XVU+o6oNPa7tKVe9U1Yf3/jx6kDUCAIDVwLgCAIBnRgbGBVLV14vIjoj8hxDCzXtt/7OInAwh/JKqvltEjoYQfu4c62FHAACeCzIwLkMXa1yx9zrGFgCA80YGxhUghPBxETkZNb9VRO7Y+/sdIvK2fS0KAACsJMYVAAA8s+KgC7hMXRdCeFREJITwqKpem+qkqu8UkXfua2UAAGDVnNe4QoSxBQDg8sYDjAMUQrhdRG4X4WueAADgwjG2AABczphCcmk8pqrHRUT2/jxxwPUAAIDVxbgCAADhAcal8iERecfe398hIh88wFoAAMBqY1wBAIDwW0gumKq+T0R+UESOichjIvILIvLbIvJ+EXmxiDwiIj8WQogDueL1sCMAAM8Fv4XkMnSxxhV762JsAQA4b6vwW0h4gDESDDIAAM8RDzDwrBhbAACei1V4gMEUEgAAAAAAMHr8FhIAwMr5vptvdG0PPPjlA6gEAABcDt7+hltc22/80f0HUAmeDd/AAAAAAAAAo8cDDAAAAAAAMHo8wAAAAAAAAKPHAwwAAAAAADB6/BrVkeBXnQEAniN+jSqeFWMLAMBzwa9RBQAAAAAAuAh4gAEAAAAAAEaPBxgAAAAAAGD0eIABAAAAAABGrzjoAgAAq+e1332DWW56n/nU55V/YdOZxXIxcV2yauraQrSuP/nqp8+nTAAAsCL+wX/2FrMcpHN9NF9zbX20PO13XJ9M414ildqf5f+7u+86jypx0PgGBgAAAAAAGD0eYAAAAAAAgNHjAQYAAAAAABg9DSEcdA0QEVVlRwBYWa/9rhtdW135y1o/HDbL2Tx3fcraz3mtsw2znAefuXHPY584Z52XmXtDCLcedBEYL8YWAFbZ217/A66tDqVrK3TTLGtbuz5r9cK/btiIlv3P9v/9Ax8+Z52XkxASA6yR4RsYAAAAAABg9HiAAQAAAAAARo8HGAAAAAAAYPR4gAEAAAAAAEaPEM+RIGgLwOXm1S+9ybWF0gZtDbLp+hQ7Z1xbW15llu95+K4LrO6yQIgnnhVjCwCXmx95422J1sosFbUPCB/mS9e27L/DLN/5mf/rgmq7HBDiCQAAAAAAcBHwAAMAAAAAAIweDzAAAAAAAMDo8QADAAAAAACMXnHQBQAALk8blW+r57VZ3i6nrk+/5tvK+eKi1QUAAFZTmdWurRnsuOFMOOb6tPng2qbaXLzCsG/4BgYAAAAAABg9HmAAAAAAAIDR4wEGAAAAAAAYPTIwAACXRCPBte1mnVnWsOP63Pf5B1zbK2955cUrDAAArKS29BkYfZSvleV+bPHJP/st1/aG17zp4hWGfcM3MAAAAAAAwOjxAAMAAAAAAIweDzAAAAAAAMDo8QADAAAAAACMnobgQ9aw/1SVHQEAeC7uDSHcetBFYLwYWwAAnosQgh50DefCNzAAAAAAAMDo8QADAAAAAACMHg8wLhFVfbOqflFVv6Sq7z7oegAAwGpjbAEAuNLxAOMSUNVcRN4jIm8RkZeLyI+r6ssPtioAALCqGFsAAMADjEvlVSLypRDCl0MIjYj8hoi89YBrAgAAq4uxBQDgiscDjEvjhSLy9actf2OvzVDVd6rqPap6z75VBgAAVhFjCwDAFa846AIuU6lfP+N+lVkI4XYRuV2EX3UGAACeFWMLAMAVj29gXBrfEJEXPW35ehH51gHVAgAAVh9jCwDAFY9vYFwafyYiN6nqDSLyTRF5u4j8o3O85gkR+dre34/tLa8a6t5f1L3/VrV26t5f+1X3d+7De2A8LmRswbm0/1a1dureX9S9/1a19v2oeyXGFTzAuARCCJ2q/nci8hERyUXkvSGEz53jNdc89XdVvSeEcOslLvOio+79Rd37b1Vrp+79tap1Y9wuZGyxqsfkqtYtsrq1U/f+ou79t6q1r2rdlwIPMC6REMJ/EpH/dNB1AACAywNjCwDAlY4MDAAAAAAAMHo8wBin2w+6gOeJuvcXde+/Va2duvfXqtaNy9eqHpOrWrfI6tZO3fuLuvffqta+qnVfdBoCv2ELAAAAAACMG9/AAAAAAAAAo8cDDAAAAAAAMHo8wBgRVX2zqn5RVb+kqu8+6Hqeiaq+V1VPqOqDT2u7SlXvVNWH9/48epA1pqjqi1T1D1T1IVX9nKr+zF77qGtX1amq/qmqfmav7l/cax913U9R1VxV/1xVP7y3vCp1f1VVP6uq96vqPXtto69dVY+o6m+q6hf2jvW/Nfa6VfWle9v5qf+2VPVnx163iIiq/tO98/JBVX3f3vk6+rpx5WBscWkxtjgYjC32F2OL/cXY4tnxAGMkVDUXkfeIyFtE5OUi8uOq+vKDreoZ/XsReXPU9m4R+VgI4SYR+dje8th0IvLPQggvE5FXi8i79rbx2GuvReS2EML3i8gtIvJmVX21jL/up/yMiDz0tOVVqVtE5D8PIdzytN+7vQq1/28i8nshhL8uIt8vZ7f9qOsOIXxxbzvfIiJ/Q0TmIvIfZeR1q+oLReSnReTWEMLNIpKLyNtl5HXjysHYYl8wtjgYjC32F2OLfcLY4jyEEPhvBP+JyN8SkY88bfnnReTnD7quZ6n3JSLy4NOWvygix/f+flxEvnjQNZ7HZ/igiPztVapdRNZE5D4R+ZurULeIXC9nL7K3iciHV+lYEZGvisixqG3UtYvIIRH5iuwFNK9K3VGtbxKRT65C3SLyQhH5uohcJSKFiHx4r/5R181/V85/jC0O5DMwtrj09TK22N+aGVvsb62MLc7xH9/AGI+nDtanfGOvbVVcF0J4VERk789rD7ieZ6WqLxGRV4jI3bICte99VfJ+ETkhIneGEFaibhH5VyLyP4nI8LS2VahbRCSIyO+r6r2q+s69trHXfqOIPC4i/27vq7W/oqrrMv66n+7tIvK+vb+Puu4QwjdF5H8RkUdE5FERORNC+H0Zed24ojC22EeMLfYNY4v9xdhiHzG2ODceYIyHJtr4HbeXgKpuiMhvicjPhhC2Drqe8xFC6MPZr8BdLyKvUtWbD7qmc1HVvyMiJ0II9x50Lc/Ta0MIr5SzX71+l6q+/qALOg+FiLxSRP7PEMIrRGRXVugrhqpaicjfFZH/96BrOR9780/fKiI3iMh3iMi6qv7EwVYFGIwt9glji/3B2OJAMLbYR4wtzo0HGOPxDRF50dOWrxeRbx1QLc/HY6p6XERk788TB1xPkqqWcnaA8eshhA/sNa9E7SIiIYTTIvKHcnae8Njrfq2I/F1V/aqI/IaI3Kaqvybjr1tEREII39r784ScnTP5Khl/7d8QkW/s/RRNROQ35eygY+x1P+UtInJfCOGxveWx1/1DIvKVEMLjIYRWRD4gIq+R8deNKwdji33A2GJfMbbYf4wt9hdji3PgAcZ4/JmI3KSqN+w9KXy7iHzogGt6Lj4kIu/Y+/s75Owc0FFRVRWRXxWRh0IIv/y0/zXq2lX1GlU9svf3mZy9sH1BRl53COHnQwjXhxBeImeP57tCCD8hI69bRERV11V186m/y9m5hw/KyGsPIXxbRL6uqi/da3qjiHxeRl730/y4/NVXPEXGX/cjIvJqVV3bu768Uc4Gm429blw5GFtcYowt9hdji/3H2GLfMbY4Bw2BbxKOhar+F3J2Xl8uIu8NIfzLAy4pSVXfJyI/KCLHROQxEfkFEfltEXm/iLxYzp54PxZCOHlQNaao6g+IyB+LyGflr+ZN/nM5O1d1tLWr6veJyB1y9rjIROT9IYR/oapXy4jrfjpV/UER+R9DCH9nFepW1Rvl7E9GRM5+dfL/DiH8yxWp/RYR+RURqUTkyyLy38jecSPjrntNzs7VvzGEcGavbRW29y+KyD+Us7+J4M9F5B+LyIaMvG5cORhbXFqMLQ4OY4v9w9hifzG2eHY8wAAAAAAAAKPHFBIAAAAAADB6PMAAAAAAAACjxwMMAAAAAAAwejzAAAAAAAAAo8cDDAAAAAAAMHo8wAAAAAAAAKPHAwwAAAAAADB6/z+MNKvhx9ZgGgAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig, axs = plt.subplots(2, 2, figsize=(15, 10))\n", "\n", "for i, driver in enumerate(sorted(ks)):\n", " k = ks[driver]\n", " axs[0, 0].bar(\n", " [i, i + 0.2],\n", " [k[\"-ve\"].mean(), k[\"+ve\"].mean()],\n", " width=0.2,\n", " color=[\"firebrick\", \"steelblue\"],\n", " )\n", "axs[0, 0].legend(\n", " [\n", " matplotlib.patches.Patch(facecolor=\"firebrick\"),\n", " matplotlib.patches.Patch(facecolor=\"steelblue\"),\n", " ],\n", " [\"Negative\", \"Positive\"],\n", ")\n", "axs[0, 0].set_xticks(np.arange(len(ks)) + 0.1)\n", "axs[0, 0].set_xticklabels([k.upper() for k in sorted(ks)])\n", "axs[0, 0].set_title(\"Kolmogorov-Smirnov\")\n", "axs[0, 0].set_ylabel(\"Mean distance (px$^2$)\")\n", "\n", "for i, driver in enumerate(sorted(euc)):\n", " k = euc[driver]\n", " axs[0, 1].bar(\n", " [i, i + 0.2],\n", " [k[\"-ve\"].mean(), k[\"+ve\"].mean()],\n", " width=0.2,\n", " color=[\"firebrick\", \"steelblue\"],\n", " )\n", "axs[0, 1].legend(\n", " [\n", " matplotlib.patches.Patch(facecolor=\"firebrick\"),\n", " matplotlib.patches.Patch(facecolor=\"steelblue\"),\n", " ],\n", " [\"Negative\", \"Positive\"],\n", ")\n", "axs[0, 1].set_xticks(np.arange(len(ks)) + 0.1)\n", "axs[0, 1].set_xticklabels([k.upper() for k in sorted(ks)])\n", "axs[0, 1].set_title(\"Euclidean\")\n", "axs[0, 1].set_ylabel(\"Mean distance (px$^2$)\")\n", "\n", "totals = {}\n", "for driver, k in ks.items():\n", " totals[driver] = np.hypot(k[\"+ve\"], k[\"-ve\"])\n", "\n", "lon, lat = np.meshgrid(rainfall.longitude.values, rainfall.latitude.values)\n", "colours = np.ma.stack(\n", " [\n", " mpl_toolkits.basemap.maskoceans(lon, lat, totals[\"enso\"]),\n", " mpl_toolkits.basemap.maskoceans(lon, lat, totals[\"sam\"]),\n", " mpl_toolkits.basemap.maskoceans(lon, lat, totals[\"iod\"]),\n", " ]\n", ").transpose(1, 2, 0)\n", "max_colour = colours.max()\n", "axs[1, 0].imshow(np.where(colours.mask, np.nan, colours / max_colour), origin=\"lower\")\n", "\n", "totals = {}\n", "for driver, k in euc.items():\n", " totals[driver] = np.hypot(k[\"+ve\"], k[\"-ve\"])\n", "\n", "colours = np.ma.stack(\n", " [\n", " mpl_toolkits.basemap.maskoceans(lon, lat, totals[\"enso\"]),\n", " mpl_toolkits.basemap.maskoceans(lon, lat, totals[\"sam\"]),\n", " mpl_toolkits.basemap.maskoceans(lon, lat, totals[\"iod\"]),\n", " ]\n", ").transpose(1, 2, 0)\n", "max_colour = colours.max()\n", "axs[1, 1].imshow(np.where(colours.mask, np.nan, colours / max_colour), origin=\"lower\")\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "\n", "## Additional information\n", "\n", "**License:** The code in this notebook is licensed under the [Apache License, Version 2.0](https://www.apache.org/licenses/LICENSE-2.0). \n", "Digital Earth Australia data is licensed under the [Creative Commons by Attribution 4.0](https://creativecommons.org/licenses/by/4.0/) license.\n", "\n", "**Contact:** If you need assistance, please post a question on the [Open Data Cube Slack channel](http://slack.opendatacube.org/) or on the [GIS Stack Exchange](https://gis.stackexchange.com/questions/ask?tags=open-data-cube) using the `open-data-cube` tag (you can view previously asked questions [here](https://gis.stackexchange.com/questions/tagged/open-data-cube)).\n", "If you would like to report an issue with this notebook, you can file one on [GitHub](https://github.com/GeoscienceAustralia/dea-notebooks).\n", "\n", "**Last modified:** December 2023" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Tags\n", "." ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "**Tags**: :index:`NCI compatible`, :index:`time series`, :index:`water`, :index:`climate data`" ] } ], "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.8.10" }, "widgets": { "application/vnd.jupyter.widget-state+json": { "state": {}, "version_major": 2, "version_minor": 0 } } }, "nbformat": 4, "nbformat_minor": 4 }