{ "cells": [ { "cell_type": "markdown", "id": "ea863e05", "metadata": {}, "source": [ "# HW#1: Extreme Rainfall Deficit in Singapore \n", "\n", "```{admonition} Objectives\n", ":class: tip\n", "\n", "This homework will help you gain a better understanding in terms of the ways how to:\n", "* Fit Generalized Extreme Value (GEV) distribution \n", "* Estimate the return level of extreme rainfall deficit\n", "\n", "Happy coding!\n", "```" ] }, { "cell_type": "markdown", "id": "af48d568", "metadata": {}, "source": [ "```{admonition} Submission Guide\n", "\n", "Deadline: **Sunday 11:59 pm, 3rd November 2024** \n", "(Note: Late submissions will not be accepted). \n", "\n", "Please upload your solutions to [Canvas](https://canvas.nus.edu.sg/courses/61921/assignments) in a Jupyter Notebook format with the name \"Homework1_StudentID.ipynb\". Make sure to write down your student ID and full name in the cell below. \n", "\n", "For any questions, feel free to contact Prof. Xiaogang HE ([hexg@nus.edu.sg](mailto:hexg@nus.edu.sg)), Haoling CHEN ([h.chen@u.nus.edu](mailto:h.chen@u.nus.edu)) or Meilian LI ([limeilian@u.nus.edu](mailto:limeilian@u.nus.edu)).\n", "\n", "```" ] }, { "cell_type": "code", "execution_count": 1, "id": "c8ae81f2", "metadata": {}, "outputs": [], "source": [ "## Fill your student ID and full name below.\n", "\n", "# Student ID:\n", "# Full name:" ] }, { "cell_type": "markdown", "id": "5abc96cb", "metadata": {}, "source": [ "**Data**:\n", "You will need to use the historical (1981-2020) daily total rainfall at Singapore's Changi station for this homework. \n", "You can create a DataFrame using Pandas by reading file \"../../assets/data/Changi_daily_rainfall.csv\"." ] }, { "cell_type": "markdown", "id": "89af3af8", "metadata": {}, "source": [ "## Q1: Calculate daily rainfall statistics (10 marks)\n", "\n", "Calculate the following statistics for daily rainfall during DJF (**D**ecember-**J**anuary-**F**ebruary): (i) mean, (ii) variance, (iii) skewness, and (iv) kurtosis.\n", "\n", "Hint: \n", "- You can filter the daily rainfall time series for DJF using Pandas' boolean filtering method. Details on filtering values can be found in the [Pandas tutorial](https://xiaoganghe.github.io/python-climate-visuals/chapters/data-analytics/pandas-basic.html).\n", "- DJF spans across two calendar years. Make sure you only include complete DJF seasons. For the period 1891 to 2020, this results in 39 complete DJF seasons, from DJF 1981-1982 to DJF 2019-2020." ] }, { "cell_type": "code", "execution_count": 2, "id": "92d51020", "metadata": {}, "outputs": [], "source": [ "# Your solutions go here.\n", "# using the + icon in the toolbar to add a cell." ] }, { "cell_type": "code", "execution_count": 1, "id": "d6f7205c", "metadata": {}, "outputs": [], "source": [ "import math\n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "\n", "from scipy.optimize import fsolve\n", "from scipy.stats import genextreme as gev" ] }, { "cell_type": "code", "execution_count": 9, "id": "36f1211e", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Daily Rainfall Total (mm)
Date
1981-01-010.0
1981-01-020.0
1981-01-030.0
1981-01-040.0
1981-01-050.0
1981-01-060.3
1981-01-0751.4
1981-01-080.0
1981-01-094.4
1981-01-100.0
\n", "
" ], "text/plain": [ " Daily Rainfall Total (mm)\n", "Date \n", "1981-01-01 0.0\n", "1981-01-02 0.0\n", "1981-01-03 0.0\n", "1981-01-04 0.0\n", "1981-01-05 0.0\n", "1981-01-06 0.3\n", "1981-01-07 51.4\n", "1981-01-08 0.0\n", "1981-01-09 4.4\n", "1981-01-10 0.0" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "rain = pd.read_csv('../../assets/data/Changi_daily_rainfall.csv', index_col = 0, header = 0, parse_dates = True)\n", "display(rain.head(10))" ] }, { "cell_type": "code", "execution_count": 3, "id": "f3718d42", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Mean 7.203409\n", "Variance 340.698022\n", "Skewness 4.909512\n", "Kurtosis 33.412606\n", "dtype: float64" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Filter daily rainfall for December-January-February (DJF)\n", "DJF_rain = rain.loc[(rain.index.month < 3) | (rain.index.month > 11)]\n", "# Select the complete DJF periods\n", "DJF_rain = DJF_rain.loc['19811201':'20200301']\n", "# Calculate the following statistics for daily rainfall during DJF: (i) mean; (ii) variance; (iii) skewness; and (iv) kurtosis.\n", "DJF_rain_stats = pd.Series([\n", " DJF_rain.mean().values[0], \n", " DJF_rain.var().values[0], \n", " DJF_rain.skew().values[0], \n", " DJF_rain.kurtosis().values[0]],\n", " index=['Mean', 'Variance', 'Skewness', 'Kurtosis']\n", " )\n", "display(DJF_rain_stats)" ] }, { "cell_type": "markdown", "id": "c03f010a", "metadata": {}, "source": [ "## Q2: Fit the GEV distribution (40 marks)\n", "\n", "Find the seasonal maximum rainfall deficit for DJF, based on the 30-day centered moving average rainfall deficit. This will result in a data set of 39 values, one value for each year. Fit a GEV distribution to the time series of seasonal maximum rainfall deficits. To do this, estimate the GEV parameters using (i) Maximum Likelihood and (ii) L-Moments, respectively. (Details on fitting a GEV distribution can be found in the [Scipy tutorial](https://xiaoganghe.github.io/python-climate-visuals/chapters/data-analytics/scipy-basic.html))\n", "\n", "Hint: The rainfall deficit is calculated by subtracting the 30-day moving average rainfall from the mean rainfall calculated in Q1.\n" ] }, { "cell_type": "code", "execution_count": 6, "id": "1fc8bc93", "metadata": {}, "outputs": [], "source": [ "# Your solutions go here.\n", "# using the + icon in the toolbar to add a cell." ] }, { "cell_type": "code", "execution_count": 11, "id": "714372c8", "metadata": {}, "outputs": [], "source": [ "# Add 'Year' column to denote the year of the days\n", "DJF_rain['Year'] = np.where(DJF_rain.index.month==12, DJF_rain.index.year.values, DJF_rain.index.year.values-1)\n", "# Calculate the 30-day moving average for each year\n", "DJF_rain_ma = DJF_rain.groupby(['Year']).rolling(30, center=True).mean()\n", "DJF_rain_ma = DJF_rain_ma.dropna(subset='Daily Rainfall Total (mm)')\n", "DJF_rain_ma.rename(columns={'Daily Rainfall Total (mm)': '30-day Moving Average Rainfall (mm)'}, inplace=True)\n" ] }, { "cell_type": "code", "execution_count": 65, "id": "fda7c3f2", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
30-day Moving Average Rainfall Deficit (mm)
YearDate
19811981-12-16-5.286591
1981-12-17-5.099924
1981-12-18-5.616591
1981-12-19-6.663258
1981-12-20-4.976591
.........
20192020-02-113.170076
2020-02-124.436742
2020-02-134.436742
2020-02-145.036742
2020-02-155.036742
\n", "

2389 rows × 1 columns

\n", "
" ], "text/plain": [ " 30-day Moving Average Rainfall Deficit (mm)\n", "Year Date \n", "1981 1981-12-16 -5.286591\n", " 1981-12-17 -5.099924\n", " 1981-12-18 -5.616591\n", " 1981-12-19 -6.663258\n", " 1981-12-20 -4.976591\n", "... ...\n", "2019 2020-02-11 3.170076\n", " 2020-02-12 4.436742\n", " 2020-02-13 4.436742\n", " 2020-02-14 5.036742\n", " 2020-02-15 5.036742\n", "\n", "[2389 rows x 1 columns]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Calculate the 30-day moving average rainfall deficit\n", "DJF_rain_def = DJF_rain['Daily Rainfall Total (mm)'].mean() - DJF_rain_ma \n", "DJF_rain_def.columns = ['30-day Moving Average Rainfall Deficit (mm)'] \n", "display(DJF_rain_def)" ] }, { "cell_type": "code", "execution_count": 64, "id": "181db785", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
DJF Maximum Rainfall Deficit (mm)
Year
19817.040076
19827.016742
19832.786742
19845.956742
19856.306742
\n", "
" ], "text/plain": [ " DJF Maximum Rainfall Deficit (mm)\n", "Year \n", "1981 7.040076\n", "1982 7.016742\n", "1983 2.786742\n", "1984 5.956742\n", "1985 6.306742" ] }, "execution_count": 64, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Calculate the maximum deficit for each year\n", "DJF_max_def = pd.DataFrame(DJF_rain_def['30-day Moving Average Rainfall Deficit (mm)'].groupby('Year').max())\n", "DJF_max_def.columns = ['DJF Maximum Rainfall Deficit (mm)']\n", "DJF_max_def.head()" ] }, { "cell_type": "code", "execution_count": 66, "id": "a8c4b8fe", "metadata": {}, "outputs": [], "source": [ "# Fit the distribution using MLE\n", "cMLE, locMLE, scaleMLE = gev.fit(DJF_max_def, method = \"MLE\")\n", "MLEGEV = gev(cMLE, loc = locMLE, scale = scaleMLE) " ] }, { "cell_type": "code", "execution_count": 67, "id": "af85744e", "metadata": {}, "outputs": [], "source": [ "# Calculate L-moments based on samples\n", "def samlmom3(sample):\n", " \"\"\"\n", " samlmom3 returns the first three L-moments of samples\n", " sample is the 1-d array\n", " n is the total number of the samples, j is the j_th sample\n", " \"\"\"\n", " n = len(sample)\n", " sample = np.sort(sample.reshape(n))[::-1]\n", " b0 = np.mean(sample)\n", " b1 = np.array([(n - j - 1) * sample[j] / n / (n - 1)\n", " for j in range(n)]).sum()\n", " b2 = np.array([(n - j - 1) * (n - j - 2) * sample[j] / n / (n - 1) / (n - 2)\n", " for j in range(n - 1)]).sum()\n", " lmom1 = b0\n", " lmom2 = 2 * b1 - b0\n", " lmom3 = 6 * (b2 - b1) + b0\n", "\n", " return lmom1, lmom2, lmom3\n", "\n", "def pargev_fsolve(lmom):\n", " \"\"\"\n", " pargev_fsolve estimates the parameters of the Generalized Extreme Value \n", " distribution given the L-moments of samples\n", " \"\"\"\n", " lmom_ratios = [lmom[0], lmom[1], lmom[2] / lmom[1]]\n", " f = lambda x, t: 2 * (1 - 3**(-x)) / (1 - 2**(-x)) - 3 - t\n", " G = fsolve(f, 0.01, lmom_ratios[2])[0]\n", " para3 = G\n", " GAM = math.gamma(1 + G)\n", " para2 = lmom_ratios[1] * G / (GAM * (1 - 2**-G))\n", " para1 = lmom_ratios[0] - para2 * (1 - GAM) / G\n", " return para1, para2, para3" ] }, { "cell_type": "code", "execution_count": 68, "id": "d495faf9", "metadata": {}, "outputs": [], "source": [ "# Fit the distribution using LMM\n", "LMM = samlmom3(DJF_max_def.values)\n", "locLMM, scaleLMM, cLMM = pargev_fsolve(LMM)\n", "LMMGEV = gev(cLMM, loc=locLMM, scale=scaleLMM) " ] }, { "cell_type": "code", "execution_count": 69, "id": "df8166ac", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Parameters fitted by Maximum Likelihood: shape=1.082, loc=5.252, scale=2.104\n", "\n", "Parameters fitted by L-Moments: shape=0.945, loc=5.164, scale=2.052\n", "\n" ] } ], "source": [ "print('Parameters fitted by Maximum Likelihood: shape={:.3f}, loc={:.3f}, scale={:.3f}\\n'\n", " .format(cMLE, locMLE, scaleMLE,))\n", "print('Parameters fitted by L-Moments: shape={:.3f}, loc={:.3f}, scale={:.3f}\\n'\n", " .format(cLMM, locLMM, scaleLMM))" ] }, { "cell_type": "markdown", "id": "3ddc466e", "metadata": {}, "source": [ "## Q3: Estimate the return level of the extreme events (20 marks)\n", "\n", "Using the GEV parameters estimated with L-Moments in Q2, estimate the rainfall deficit for events with return periods of 50 years, 100 years, and 1000 years." ] }, { "cell_type": "code", "execution_count": 14, "id": "43c2f39c", "metadata": {}, "outputs": [], "source": [ "# Your solutions go here.\n", "# using the + icon in the toolbar to add a cell." ] }, { "cell_type": "code", "execution_count": 70, "id": "af770e8b", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
DJF Maximum Rainfall Deficit (mm)
Return period
507.281076
1007.307344
10007.332279
\n", "
" ], "text/plain": [ " DJF Maximum Rainfall Deficit (mm)\n", "Return period \n", "50 7.281076\n", "100 7.307344\n", "1000 7.332279" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "T = np.arange(2, 1001)\n", "rainDef_lmm = LMMGEV.ppf(1 - 1.0 / T)\n", "rainDef_lmm = pd.DataFrame(rainDef_lmm, index = T, \n", " columns = [\"DJF Maximum Rainfall Deficit (mm)\"])\n", "rainDef_lmm.index.name = \"Return period\"\n", "\n", "display(rainDef_lmm.loc[[50, 100, 1000]])" ] }, { "cell_type": "markdown", "id": "0006273e", "metadata": {}, "source": [ "## Q4: Test the goodness-of-fit (30 marks)\n", "\n", "In this task, you will compare how different distributions fit the same dataset and interpret the results using both statistical analyses. \n", "- Repeat the distribution fitting as in Q2, but this time using a [normal distribution](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.norm.html) and the Maximum Likelihood method. (5 marks)\n", "- Use the Kolmogorov-Smirnov (KS) test to evaluate the goodness-of-fit for both the normal distribution and the GEV distribution you obtained in Q2. (Details on the KS test can be found in the [Scipy tutorial](https://xiaoganghe.github.io/python-climate-visuals/chapters/data-analytics/scipy-basic.html)) (10 marks)\n", "- Based on the KS test results, discuss how well each distribution (Normal and GEV) fits the data. (15 marks)\n", "\n", "**Bonus (10 marks):**\n", "- Plot the CDF (Cumulative Distribution Function) to visually compare the fitted normal distribution, the GEV distribution from Q2, and the empirical distribution derived from the data. Compare the behavior of the two distributions at different return periods. Are the KS statistic results consistent with your observations from the CDF plot?\n", "\n", "\n", "Hint: You can reuse the empirical distribution estimation and CDF plotting code from the [Scipy tutorial](https://xiaoganghe.github.io/python-climate-visuals/chapters/data-analytics/scipy-basic.html)." ] }, { "cell_type": "code", "execution_count": 16, "id": "80ee34d7", "metadata": {}, "outputs": [], "source": [ "# Your solutions go here.\n", "# using the + icon in the toolbar to add a cell." ] }, { "cell_type": "code", "execution_count": 17, "id": "587e1357", "metadata": {}, "outputs": [], "source": [ "# Fit the normal distribution using MLE\n", "from scipy.stats import norm\n", "locMLE, scaleMLE = norm.fit(DJF_max_def, method = \"MLE\")\n", "MLEnorm = norm(loc = locMLE, scale = scaleMLE) " ] }, { "cell_type": "code", "execution_count": 18, "id": "f7dd7bfa", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Goodness-of-fit for GEV distribution: statisic is 0.077, pvalue is 0.961\n", "Goodness-of-fit for Normal dirstribution: statisic is 0.164, pvalue is 0.217\n" ] } ], "source": [ "from scipy.stats import kstest\n", "GEVKS = kstest(DJF_max_def[\"DJF Maximum Rainfall Deficit (mm)\"], MLEGEV.cdf)\n", "print('Goodness-of-fit for GEV distribution: statisic is {:.3f}, pvalue is {:.3f}'.format(GEVKS.statistic, GEVKS.pvalue))\n", "NORMKS = kstest(DJF_max_def[\"DJF Maximum Rainfall Deficit (mm)\"], MLEnorm.cdf)\n", "print('Goodness-of-fit for Normal dirstribution: statisic is {:.3f}, pvalue is {:.3f}'.format(NORMKS.statistic, NORMKS.pvalue))\n" ] }, { "cell_type": "code", "execution_count": 19, "id": "9b35c747", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGwCAYAAAB7MGXBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACW00lEQVR4nOzdd3xN5x/A8c/NzZ4kkUUIIvbee6+0FLXVrkpRNYri11LULopatYrWqE3V3pvUHrEFSUSQvXPP74/LTW8zJCQyfN+v13055znPec5zbiL3e5/zDJWiKApCCCGEEFnEIKsrIIQQQogPmwQjQgghhMhSEowIIYQQIktJMCKEEEKILCXBiBBCCCGylAQjQgghhMhSEowIIYQQIksZZnUF0kKj0eDn54eVlRUqlSqrqyOEEEKINFAUhbCwMFxcXDAwSLn9I0cEI35+fri6umZ1NYQQQgjxFh49ekSBAgVSPJ4jghErKytAezPW1tZZXBshhBBCpEVoaCiurq66z/GU5Ihg5PWjGWtrawlGhBBCiBzmTV0spAOrEEIIIbKUBCNCCCGEyFISjAghhBAiS0kwIoQQQogsJcGIEEIIIbKUBCNCCCGEyFISjAghhBAiS0kwIoQQQogsJcGIEEIIIbKUBCNCCCGypSdPnjB69GiKFCmCnZ0dNWrUYNmyZcTFxWV11UQGS3cwcvToUVq1aoWLiwsqlYqtW7e+8ZwjR45QuXJlTE1NKVKkCIsWLXqbugohhPhAXL58mQoVKrBw4UJatGjBN998g52dHf369eOjjz4iOjo6q6soMlC6g5GIiAjKly/P/Pnz05T//v37eHp6UrduXS5cuMCYMWMYPHgwmzZtSndlRcbo1asXbdq0eWO+tAabaeXm5sacOXMyrDwhRO6k0Who3749rq6u3LlzhwULFjB69Gj++usvDhw4wLFjx5g4cWJWV1NkJOUdAMqWLVtSzTNy5EilRIkSemn9+/dXatSokebrhISEKIASEhLyNtV8oydPnij/+9//lNKlSyuFCxdWWrVqpezcuVPRaDSZcj1FUZSePXsqQJJX8+bNM+2arwUHBysvX758Yz5/f38lOjo6w65bqFAhZfbs2anmCQkJUcaMGaMUL15cMTExURwdHZXGjRsrmzZt0v086tevr3u/jI2NFRcXF+Xjjz9WNm3alKS85N7j2rVrZ9g9CSEy3q5duxRAOX36dLLHhwwZotjb22fo3yeROdL6+Z3pfUZOnTpFs2bN9NKaN2/O+fPnU3zuFxMTQ2hoqN4rs5w7d44yZcrw888/U716dTp27MiTJ0/4+OOP8fLyQlGUTLt2ixYt8Pf313utXbs20673mo2NDXny5EnxeGxsLABOTk6YmJhken1eCw4OplatWqxatYrRo0fzzz//cPToUTp16sTIkSMJCQnR5e3Xrx/+/v7cuXOHTZs2UapUKTp37swXX3yRpNwVK1bovcfbt29/b/ckhEi/06dP4+zsTLVq1ZI93rZtW4KCgrhz5857rpnILJkejAQEBODo6KiX5ujoSHx8PEFBQcmeM2XKFGxsbHQvV1fXTKlbdHQ0n3zyCcWLF+fBgwcsW7aMqVOncv78eZYvX86SJUtYunRpplwbwMTEBCcnJ71X3rx5dcdVKhWLFy/m448/xtzcnJIlS3Lq1Cnu3LlDgwYNsLCwoGbNmty9e1d3zvjx46lQoQKLFy/G1dUVc3NzOnToQHBwsC7Pfx/TNGjQgEGDBjFs2DDs7e1p2rSp7vr/fkzz+PFjOnfujK2tLRYWFlSpUoUzZ84AcPfuXT755BMcHR2xtLSkatWq7N+/P13vx5gxY3jw4AFnzpyhZ8+elCpVCg8PD/r168fFixextLTU5TU3N8fJyQlXV1dq1KjBtGnTWLx4Mb/++muS6+bJk0fvPba1tU1XvYQQ75darSY2NjbFL4MxMTEAGBoavs9qiUz0XkbTqFQqvf3Xv2D/TX9t9OjRhISE6F6PHj3KlHpt3LgRf39/Vq5cqfcBpVKp6N27N+3atePnn3/O1NaRN5k4cSI9evTg4sWLlChRgq5du9K/f39Gjx7N+fPnARg0aJDeOXfu3GHDhg3s2LGD3bt3c/HiRQYOHJjqdX777TcMDQ05ceIEixcvTnI8PDyc+vXr4+fnx/bt27l06RIjR45Eo9Hojnt6erJ//34uXLhA8+bNadWqFb6+vmm6T41Gw7p16+jWrRsuLi5JjltaWr7xD0/Pnj3JmzcvmzdvTtM1hRDZU9OmTXn+/Dm7d+/meXgMdwLDiE/Q6I6vWbMGNzc33N3ds7CWuYhGA8GPIORxllUh08NKJycnAgIC9NICAwMxNDTEzs4u2XNMTEzey+OBY8eOUa5cOYoXL57s8Q4dOtClSxeCg4P1Wiwyys6dO/W+7QOMGjWK7777Trffu3dvOnbsqDtWs2ZNvvvuO5o3bw7A119/Te/evfXKiI6O5rfffqNAgQIAzJs3j48++oiffvoJJyenZOvi7u7O9OnTU6zrH3/8wbNnzzh37pwucPv3H4Ly5ctTvnx53f6kSZPYsmUL27dvTxIsJScoKIiXL19SokSJN+ZNiYGBAR4eHjx48EAvvUuXLqjVat3+mjVr0tSBVwiRNWrUqEHNmjXp378/A2evZ9H5lxgbGvBT+3L4ndnJqlWrmDNnjt7/a5EGEc8hyAeCbkHQbXh+B17ch5cPICEGqvSBj2dnSdUyPRipWbMmO3bs0Evbu3cvVapUwcjIKLMvnyqVSqX7Zp+c18dSasF5Vw0bNmThwoV6af99hFCuXDnd9uvHXWXLltVLi46OJjQ0FGtrawAKFiyoC0RA+zPQaDT4+PikGIxUqVIl1bpevHiRihUrpviIIyIigh9++IGdO3fi5+dHfHw8UVFRaW4ZeVNrWVopipKkjNmzZ9OkSRPdvrOz8ztdQwiRuVQqFRs3bqRp06bM+PUPrCq2JDZew5D+vbh3dj9eXl4MHjw4q6uZfSXEwbOb4H8ZAq5A4HXtK+JZ6ue9uPd+6peMdAcj4eHhep2G7t+/z8WLF7G1taVgwYKMHj2aJ0+esGrVKgC8vLyYP38+w4YNo1+/fpw6dYply5a9l46ab9KwYUMWL17M1atXKVOmTJLja9eupUKFCql29nwXFhYWb2xm/HfA9vpDNrm01IKq13lS+6C3sLBItR5mZmapHh8xYgR79uxh5syZuLu7Y2ZmRvv27XWdYd8kX7585M2blxs3bqQpf3ISEhK4ffs2VatW1Ut3cnKS5lwhchhzc3OqV6/O7jg3XdrTWxf46quvmDNnTqZ9ScxxFEUbRDw+D0/OwxNvbQCSkLa/vRiaQl43sC0CBVL/UpqZ0h2MnD9/noYNG+r2hw0bBmif169cuRJ/f3+9b8OFCxdm165dDB06lF9++QUXFxfmzp3Lp59+mgHVfzdt27alYMGC9OjRg127dulaDRRFYf78+ezcuZOVK1dmbSXfgq+vL35+frq+F6dOndI9wnhb5cqVY+nSpbx48SLZ1pFjx47Rq1cv2rZtC2iD1v8+LkmNgYEBnTp1YvXq1YwbNy5Jv5GIiAhMTExS7Tfy22+/8fLly2zxuyWEeHthYWE0atSI+/fvY9+/DXGAvZmKRm1bM2/ePNRqNbNnZ83jhCyn0cCzG3D/KDw8Cb6nISLwzeeZ24NjKbAvDvmKg30xsHMHKxcwyPrJ2NMdjDRo0CDVDp3JfXjXr1+ff/75J72XynTGxsbs2LGDZs2a4ebmRuvWrXFwcGDfvn3cunWLoUOH0qNHj0y7fkxMTJL+NIaGhtjb279TuaampvTs2ZOZM2cSGhrK4MGD6dixY4qPaNKiS5cuTJ48mTZt2jBlyhScnZ25cOECLi4u1KxZE3d3dzZv3kyrVq1QqVR89913qbbWJGfy5MkcPnyY6tWr8+OPP+oe5R07dowpU6Zw7tw5XStVZGQkAQEBxMfH8+TJEzZv3szs2bP58ssv9YJlIUTOM3/+fG7cuMHOgyfpu80PgPKF8rFs3HLKlSvH0KFD6d27t95j7Fwt1B/uHoA7B7RBSGTyI1F17NzBuQI4lwOnsuBYFizzvZeqvq0PflxUuXLluHbtGsuWLWPLli34+PhQuXJlfv31V+rVq5ep1969e3eS/gvFixfn5s2b71Suu7s77dq1w9PTkxcvXuDp6cmCBQveqUxjY2P27t3L8OHD8fT0JD4+nlKlSvHLL78A2n4Zffr0oVatWtjb2zNq1Kh0zw+TN29eTp8+zdSpU5k0aRIPHz4kb968lC1blhkzZmBjY6PL++uvv/Lrr79ibGyMnZ0dlStXZv369bqWGSFEzrV06VK6dOkCeVwAbTBS3MkKgIEDBzJjxgyWLl3K3Llzs7CWmUhRwO8C+OwCn7/h6dWU85pYg2s1cK2ufcziUhHMMn7ARWZTKVk5bjWNQkNDsbGxISQkRNdJUyRv/PjxbN26lYsXL2Z1VYQQIt0URUGtVrNw4ULiPRozY48PAD93rsAnFfID0Lp1azQaDTt37szKqmYsjQYen4VrW+DGDgh9knw+E2soVBsK1wO3OuBYGgyy76iitH5+f/AtI0IIIbIPlUqFra0tDx48INw6TJdewkn7QaYoCvfv36dSpUpZVcWM9fQaXFoLVzenHIC4VAT3puDeBPJXBnXu++jOfXckhBAiR+vSpQvLly+nzAjtUiJGahVF8mlH/B04cICrV68yc+bMrKziu4l8AZfXw8XftSNf/kttDEUaQHFPKN4SrN6+v19OIY9phBBCZCv379+nSrXqWPVeCgZqSjhZsdWrGhs2bODrr7+mXLlyHDp0CINsMAokzRQFHp2F88u1j2ISYvSPGxhC0UZQuh2U8ARTm+TLyWHkMY0QQogcqXDhwizdsJOhe7STdF04tBOrEY2Jj4+nTZs2rFy5MucEIvEx2kcwZxaC/6Wkx/NXhvJdtEGIRfKzkn8IJBgRQgiR7Ry+eAt4NSokxB+VSoVKpaJAgQJvnKQxW4gOgXNL4fSipPOAmOaBCl2hUk9wePslMHITCUaEEEJkK3/++Sertu3HpkYHANYt+okq+RewZMkSRo0ahaGhYfad9Cz8GZxeoA1EYv4zvYFzBajuBaXbgFHqs1p/aCQYEUIIkW0oisK4cePI3/Arwl+lFXeywsrKjOHDhxMbG8u4ceP49ttvdet1ZQuRL+DEz3B2CcRFJqarDKDEx1BzoHYuEJnGPlk55KGbEEKID8G1a9e4ceMGRvYFAbAyNcTZxlR33MvLC0VR2Lp1axbV8D9iwuHQFJhTDk7MSQxE1MbaxzCDzkOn1VCwhgQiqZCWESGEENlGSEgIBiYWvHw12KSkk7Xeonh58+bF0tKSkJCQLKrhKwnxcGE1HJqs3ydEbQyVe0OdIWDtkuLpQp+0jOQiK1euzLQVhjNLTqrzd999xxdffJGldWjfvj2zZs3K0joIkZmKFCmCsUNh3f7raeBfu3LlCsHBwRQrVux9Vy3RvcOwqDbsHJIYiBgYaoOQwRfAc7oEIukkwUgW6dWrl653+L9fLVq0SNP5bm5uzJkzRy+tU6dO3Lp1KxNqqy8rAohDhw7h6emJnZ0d5ubmlCpViuHDh/PkiXbGwsOHD+veQwMDA2xsbKhYsSIjR47E399fr6zx48cn+97v378/xes/ffqUn3/+mTFjxuilBwQE8PXXX+Pu7o6pqSmOjo7UqVOHRYsWERmZ+NzYzc0t2WtOnToVb29vVCoVx48fT/bazZs3p3Xr1gB8//33/Pjjj+le90eInMLZ2ZlKjT7W7f87GImPj2fs2LE4OTnx8ccfJ3d65gp+BBt6wqpP4Nm/1hAr2RoGnoVWc8CmwPuvVy4gj2myUIsWLVixYoVemomJyVuXZ2ZmhplZ7uuhvXjxYgYMGEDPnj3ZtGkTbm5u+Pr6smrVKn766Se9lgIfHx+sra0JDQ3ln3/+Yfr06SxbtozDhw9TtmxZXb7SpUsnCT5sbW1TrMOyZcuoWbMmbm5uurR79+5Ru3Zt8uTJw+TJkylbtizx8fHcunWL5cuX4+LiogsiACZMmEC/fv30yrWyssLCwoLy5cuzYsUK6tSpo3f80aNH7N+/n82bNwPahR3d3Nz4/fff+fLLL9P+JgqRg1Rs8BFPbmofw1w++jenlPLcu3ePefPm4e3tzebNmzEyMnp/FUqI146QOTxFv3Nq/srQfLK2P4h4N0oOEBISogBKSEhIVlclw/Ts2VP55JNPUs0zbtw4xdXVVTE2NlacnZ2Vr776SlEURalfv74C6L0URVFWrFih2NjY6J1fvnx5ZdmyZYqrq6tiYWGheHl5KfHx8cq0adMUR0dHJV++fMqkSZP0rvvTTz8pZcqUUczNzZUCBQooX375pRIWFqYoiqIcOnQoybXHjRunKIqixMTEKCNGjFBcXFwUc3NzpVq1asqhQ4f0yl6xYoXi6uqqmJmZKW3atFFmzpypV+f/evTokWJsbKwMGTIk2eMvX77Uq9fr/dciIyOV4sWLK7Vr107yvqRH2bJllfnz5+ulNW/eXClQoIASHh6e7DkajUa3XahQIWX27Nkplj937lzF0tIySVkTJkxQHB0dlbi4OF3a+PHjlbp166ar/kJkZxqNRjl48KDSuXNnpVy5ckoxr4VKoVE7lUKjdipG5ta6vzX169dXDh8+/H4r53dRURbVVZRx1omvaUUU5Z/VipKQ8H7rkgOl9fM7V7aMtJp3nGdhMW/OmMHyWZmw46s6b86YBhs3bmT27NmsW7eO0qVLExAQwKVL2tn7Nm/eTPny5fniiy+SfNP+r7t37/L333+ze/du7t69S/v27bl//z4eHh4cOXKEkydP0qdPHxo3bkyNGtro3sDAgLlz5+Lm5sb9+/cZMGAAI0eOZMGCBdSqVYs5c+bw/fff4+OjXU3T0tISgN69e/PgwQPWrVuHi4sLW7ZsoUWLFly5coVixYpx5swZ+vTpw+TJk2nXrh27d+9m3Lhxqdb/zz//JDY2lpEjRyZ7/E2Pi8zMzPDy8mLo0KEEBgbi4OCQav7kvHz5kqtXr1KlShVd2vPnz9m7dy+TJ09OcQImVTp6znfr1o0RI0bw559/0qtXL0A7xHHlypX07NkTQ8PE/6rVqlVjypQpxMTEvFNLmhDZgaIoDBo0iAULFlCyZEkaNGzIHnPtkF0l/DlHD+zBwsICW1tb8ufP//4qFh8LR6bC8TmgJLxKVEG1ftBwLJjleX91+QDkymDkWVgMAaHRWV2NN9q5c6fug/y1UaNG8d133+Hr64uTkxNNmjTByMiIggULUq1aNUD7OEGtVmNlZYWTU+oLKGk0GpYvX46VlRWlSpWiYcOG+Pj4sGvXLgwMDChevDjTpk3j8OHDumBkyJAhuvMLFy7MxIkT+fLLL1mwYAHGxsbY2NigUqn0rn337l3Wrl3L48ePcXHRdtz65ptv2L17NytWrGDy5Mn8/PPPNG/enG+//RYADw8PTp48ye7du1Os/+3bt7G2tsbZ2Tntb+x/lCihneHwwYMHumDkypUreu99qVKlOHv2bLLnP3z4EEVRdPcFcOfOHRRFoXjx4np57e3tiY7W/u4NHDiQadOm6Y6NGjWK//3vf3r5d+7cSYMGDbC1taVNmzasWLFCF4wcPnyYe/fu0adPH71z8ufPT0xMDAEBARQqVCg9b4UQ2c6SJUtYsGABixYt4osvvuBJcBS7ph0CwDDiKT16TOLGjRuo1er3V6mAq7DFC57+axG7fCWh9Txwrfr+6vEByZXBSD6rrPm2mN7rNmzYkIULF+qlve630KFDB+bMmUORIkVo0aIFnp6etGrVSu8bclq4ublhZZXYAczR0RG1Wq23roOjoyOBgYlD0w4dOsTkyZO5fv06oaGhxMfHEx0dTURERIqtAP/88w+KouDh4aGXHhMTg52ddr2FGzdu0LZtW73jNWvWTDUYURQlXS0MKZUB+i0VxYsXZ/v27br91FoYoqKiADA1NU1y7L91O3v2LBqNhm7duhETo986N2LECF2g8dq/v+n17duXZs2acefOHdzd3Vm+fDm1a9dOEvC87hf07w6yQuREiqIwZ84cOnToQP/+/QG45pfYObt1vSrMWTmSv/76S6//VabRaODUPDg4CRJitWkGRlBvBNQZCobGmV+HD1SuDEYy6lFJZrOwsMDd3T3ZY66urvj4+LBv3z7279/PgAEDmDFjBkeOHElXx63/5lWpVMmmaTQaQNsK4OnpiZeXFxMnTsTW1pbjx4/Tt29f4uLiUryORqNBrVbj7e2d5BvM6xYI5S0WiPbw8CAkJAR/f/+3bh25ceMGgF7nU2Nj4xTf+/+yt7cHtI9r8uXLB4C7uzsqlYqbN2/q5S1SpAhAsh2J7e3tU71mkyZNKFSoECtXrmTkyJFs3ryZ+fPnJ8n34sULAF1dhMipAgICuHnzJj/++KMu7crjxPlDmlctyY6iRTl06FDmByPhgbD5C7h3KDHNoRS0XQzO5TL32kKG9mZnZmZmtG7dmrlz53L48GFOnTrFlSvaZkNjY2MSEhLeUEL6nT9/nvj4eH766Sdq1KiBh4cHfn5+enmSu3bFihVJSEggMDAQd3d3vdfrxzmlSpXi9OnTeuf9d/+/2rdvj7GxMdOnT0/2eHBwcKrnR0VFsWTJEurVq/fWH95FixbF2tqa69ev69Ls7Oxo2rQp8+fPJyIi4q3K/S+VSkXv3r357bff+OOPPzAwMKBjx45J8l29epUCBQrogiQhcqrXX1D+/QXm8pPEYKRcARvUavVbfZFJl7sHYWHtfwUiKqj9NXxxWAKR90SCkSz0+rn/v19BQUGAdi6PZcuWcfXqVe7du8fq1asxMzPT9RFwc3Pj6NGjPHnyRHdORihatCjx8fHMmzdPd91Fixbp5XFzcyM8PJwDBw4QFBREZGQkHh4edOvWjR49erB582bu37/PuXPnmDZtGrt27QJg8ODB7N69m+nTp3Pr1i3mz5+f6iMa0LYQzZ49m59//pm+ffty5MgRHj58yIkTJ+jfvz8TJ07Uyx8YGEhAQAC3b99m3bp11K5dm6CgoCSPw9LDwMCAJk2aJJkHZMGCBcTHx1OlShXWr1/PjRs38PHxYc2aNdy8eTNJC1FYWFiSn/d/5wvp3bs3fn5+jBkzhs6dOyf7WOzYsWM0a9bsre9HiOzCycmJokWL8ueffwLa4OTy42AA7C1NCHzgw61bt6hbt27mVECjgcPTYHW7xMnLLJ2gxzZoOgEMpYP4e5OJI3oyTG4d2st/hsgCSvHixRVFUZQtW7Yo1atXV6ytrRULCwulRo0ayv79+3Xnnzp1SilXrpxiYmLyxqG9/73uf4cU169fX/n66691+7NmzVKcnZ0VMzMzpXnz5sqqVauSDJv18vJS7Ozs9Ib2xsbGKt9//73i5uamGBkZKU5OTkrbtm2Vy5cv685btmyZUqBAAcXMzExp1arVG4f2vrZv3z6lefPmSt68eRVTU1OlRIkSyjfffKP4+fkpiqI/5FilUilWVlZK+fLllREjRij+/v56Zb3N0N7du3cr+fPnVxL+M5TPz89PGTRokFK4cGHFyMhIsbS0VKpVq6bMmDFDiYiI0OUrVKhQsj/v/v37J7lWs2bNFEA5efJkkmNRUVGKtbW1curUqXTVX4jsas6cOYqBgYGydu1axfd5hG5Ib7fFx5XKlSsrBQsWVGJjYzP+wpEvFGVNB/0hu6s/VZTwZxl/rQ9YWj+/VYqS2e1f7y40NBQbGxtCQkKwtrbO6uqID5CiKNSoUYMhQ4bQpUuXLKvHL7/8wrZt29i7d2+W1UGIjJSQkECvXr1Ys2YNpVt8Rnj5zgBEnv0T5cpf7Nu3jwoVKmTsRZ9eh3Vd4OUD7b7KABp9B7WHgIE8MMhIaf38lnddiDRQqVQsWbKE+Pj4LK2HkZER8+bNy9I6CJGRNBoNJUqUwM7OjidRiWMqyhbIw5UrVzI+EPHZDcuaJgYi5nbw2WaoO0wCkSwkLSNCCCGyRHx8PG3btmXPnj307NmTh0U+4VaIdrj84/ndGTNsEBMmTMiYiykKnJwH+75H+5QUcK4AndZAHteMuYZIQlpGhBBCZGurVq3ir7/+Yvv27SxevAT/GO20A07WpkwY8w0TJ07UjSB8JwnxsONr2PcdukCkdFvo/bcEItmEBCNCCCGyxOLFi/H09KRFixY8fBFJWLT2MWjZAjaMGDECZ2dnFi9e/G4XiY2AdV3hn98S0xqMhvYrwNj83coWGUaCESGEEFni+vXrNGzYEEA3pBegXH4bjIyMqFevnt78PukW/gxWfgy392j3DYzg02XQ4Ft4x5mdRcbKlTOwCiGEyP6srKx4+vQpAJf/NfNqOdc8gHaG1rfuJxjsC6vawIu72n0TG+i8BgrXe4cai8wiLSNCCCGyRLt27Vi1ahXh4eF608CXzW/D9evXOXLkCO3bt09/wUG3YXnLxEDEOj/02S2BSDYmwYgQQogsMWTIECIiIvi4VWsuP34JQIG8Zlz75wwff/wxHh4eyS6JkCr/y7C8BYQ+1u7buUPfveBYKoNrLzKSBCNCCCGyhLu7O3///TfXHz8nOl47ysX34nHq1auHtbU1e/fuTXa17BQ98YbfPobIV0tkOJWF3rvBpkAm1F5kJAlGcpGVK1eSJ0+erK5GuuSkOn/33Xd88cUXWV2Nd3b48GFUKpVukcGdO3dSsWJF3crNQrwvL1++ZNasWUSaJC76GHL/MuXKlWPbtm26tbjS5LE3rGoL0a8e97hWh547wVJWt84JJBjJIr169UKlUiV5tWjRIk3nu7m5MWfOHL20Tp06cevWrUyorb6sCCAOHTqEp6cndnZ2mJubU6pUKYYPH86TJ0+AxA9YlUqFgYEBNjY2VKxYkZEjR+Lv769X1vjx45N97/fv35/i9Z8+fcrPP//MmDFjdGmvf4ZTp07Vy7t161ZUOain/scff4xKpeKPP/7I6qqID0hsbCwtW7bkyJEjtOzWX5f+7RddCQkJoVGjRrx8+TJthT0+D6vbQMyrQMStLnTfAmZ5MrzeInNIMJKFWrRogb+/v95r7dq1b12emZkZDg4OGVjD7GHx4sU0adIEJycnNm3axPXr11m0aBEhISH89NNPenl9fHzw8/Pj3LlzjBo1iv3791OmTJkkEyeVLl06yXtfr17KnduWLVtGzZo1cXNz00s3NTVl2rRpaf+jmUaxsbEZWt6b9O7dW6aZF+/V5s2bOXPmDH/99RfxNvl16QM6t+LQoUM8efKEJUuWvLmgJ//A6rYQ82oFbLe60HU9GCdd8VpkXxKMZCETExOcnJz0Xnnz5tUdHz9+PAULFsTExAQXFxcGDx4MQIMGDXj48CFDhw7VfauHpC0W48ePp0KFCixfvpyCBQtiaWnJl19+SUJCAtOnT8fJyQkHBwd+/PFHvXrNmjWLsmXLYmFhgaurKwMGDCA8PBzQtkD07t2bkJAQ3bXHjx8PaD9AR44cSf78+bGwsKB69eocPnxYr+yVK1dSsGBBzM3Nadu2Lc+fP0/1PXr8+DGDBw9m8ODBLF++nAYNGuDm5ka9evVYunQp33//vV5+BwcHnJyc8PDwoHPnzpw4cYJ8+fLx5Zdf6uUzNDRM8t4bGxunWI9169bRunXrJOmvg6QpU6akeh+bNm2idOnSmJiY4ObmliSIcnNzY9KkSfTq1QsbGxv69eun+3nu3LmT4sWLY25uTvv27YmIiOC3337Dzc2NvHnz8tVXX5GQkKAra82aNVSpUgUrKyucnJzo2rUrgYGBqdavdevWnD17lnv37qWaT4iMsmbNGurWrUvFKlW59kQbSBS2t8DG3IjChQvToUMHVq9enXohT6/BmnaJgUjhetB1gwQiOVDunGdkcX0IT/2Pb6awdID+RzKkqI0bNzJ79mzWrVtH6dKlCQgI4NKlS4D2G0X58uX54osv6NevX6rl3L17l7///pvdu3dz9+5d2rdvz/379/Hw8ODIkSOcPHmSPn360LhxY2rUqAGAgYEBc+fOxc3Njfv37zNgwABGjhzJggULqFWrFnPmzOH777/Hx8dHe9uWloD22/WDBw9Yt24dLi4ubNmyhRYtWnDlyhWKFSvGmTNn6NOnD5MnT6Zdu3bs3r2bcePGpVr/P//8UxfkJOdNj4vMzMzw8vJi6NChBAYGvlXL0cuXL7l69SpVqlRJckytVjN58mS6du3K4MGDKVAgaUc5b29vOnbsyPjx4+nUqRMnT55kwIAB2NnZ0atXL12+GTNm8N133/G///0PgOPHjxMZGcncuXNZt24dYWFhtGvXjnbt2pEnTx527drFvXv3+PTTT6lTpw6dOnUCtEHhxIkTKV68OIGBgQwdOpRevXqxa9euFO+xUKFCODg4cOzYMYoUKZLu90iI9Hr27BllypThml8osQna/koVC+bRHS9evHjqq1MH3dHOIxL1qlWyUB3osl5mVc2hcmcwEh4IYX5ZXYs32rlzp+6D/LVRo0bx3Xff4evri5OTE02aNMHIyIiCBQtSrVo1AGxtbVGr1bpvvqnRaDQsX74cKysrSpUqRcOGDfHx8WHXrl0YGBhQvHhxpk2bxuHDh3XByJAhQ3TnFy5cmIkTJ/Lll1+yYMECjI2NsbGxQaVS6V377t27rF27lsePH+Pi4gLAN998w+7du1mxYgWTJ0/m559/pnnz5nz77bcAeHh4cPLkSXbv3p1i/W/fvo21tTXOzs5pf2P/o0SJEgA8ePBAF4xcuXJF770vVaoUZ8+eTfb8hw8foiiK7r7+q23btlSoUIFx48axbNmyJMdnzZpF48aN+e677wDtfV+/fp0ZM2boBSONGjXim2++0e0fP36cuLg4Fi5cSNGiRQFo3749q1ev5unTp1haWup+pocOHdIFI3369NGVUaRIEebOnUu1atUIDw9P8vv2b/nz5+fBgwcpHhciI7m6unLhwgX+eZj4iLNSwcSWYW9vbwoWLJj8ycG+sOoTiHj1pTN/Fei6TgKRHCx3BiOWWdRvIp3XbdiwIQsXLtRLs7W1BaBDhw7MmTOHIkWK0KJFCzw9PWnVqhWGhun7kbm5uWFlZaXbd3R0RK1WY/CvpbIdHR31mvEPHTrE5MmTuX79OqGhocTHxxMdHU1ERAQWFsk3f/7zzz8oioKHh4deekxMDHZ2dgDcuHGDtm3b6h2vWbNmqsGIoijv3Bn09cLU/y6nePHibN++XbdvYmKS4vlRUVEAqQ4xnDZtGo0aNWL48OFJjt24cYNPPvlEL6127drMmTOHhIQE1Go1QLItL+bm5rpABLQ/Kzc3N72g4r8/vwsXLjB+/HguXrzIixcvdKNkfH19KVUq5bkWzMzMiIyMTPG4EBmpT58+fPTRR+w8dY3XPQYqF9IGI97e3mzfvp358+cnPTHiOaxulziPiGNZ+GwjmFglzStyjNwZjGTQo5LMZmFhgbu7e7LHXF1d8fHxYd++fezfv58BAwYwY8YMjhw5gpGRUZqv8d+8KpUq2bTXH1gPHz7E09MTLy8vJk6ciK2tLcePH6dv377ExcWleB2NRoNarcbb21v34fra6w/O10FBenh4eBASEoK/v/9bt47cuHEDQK/zqbGxcYrv/X/Z22uHHb58+ZJ8+ZIfJlivXj2aN2/OmDFj9Fo7IPmAKrn3IrlAL70/v4iICJo1a0azZs1Ys2YN+fLlw9fXl+bNm7+xU+yLFy9SvD8hMlqLFi1o3bo13g9forayw8xIhUHYU6atWsiUKVOoUqVKkv9LxEbAHx3g+W3tvp37q1EzeZOUL3KW3BmM5BJmZma0bt2a1q1bM3DgQEqUKMGVK1eoVKkSxsbGep0WM8r58+eJj4/np59+0rWebNiwQS9PcteuWLEiCQkJBAYGUrdu3WTLLlWqFKdPn9ZL++/+f7Vv355vv/2W6dOnM3v27CTHg4ODU+03EhUVxZIlS6hXr95bf9AWLVoUa2trrl+/nqTl59+mTp1KhQoVkuQpVaoUx48f10s7efIkHh4eSQK3d3Xz5k2CgoKYOnUqrq7apdHPnz//xvOio6O5e/cuFStWzND6CPFaWFgYv//+O3v27CE+Pp5q1arxv8kz6bBaOx3Bi9sXKFH8I0xMTOjWrRuzZ8/GzMwssYCEONjQUzuxGYClE3y2WeYRySUkGMlCMTExBAQE6KUZGhpib2/PypUrSUhIoHr16pibm7N69WrMzMx0kwC5ublx9OhROnfujImJie7b+7sqWrQo8fHxzJs3j1atWnHixAkWLVqkl8fNzY3w8HAOHDhA+fLlMTc3x8PDg27dutGjRw9++uknKlasSFBQEAcPHqRs2bJ4enoyePBgatWqxfTp02nTpg179+5N9RENaFuIZs+ezaBBgwgNDaVHjx64ubnx+PFjVq1ahaWlpd7IlMDAQKKjowkLC8Pb25vp06cTFBTE5s2b3/o9MTAwoEmTJhw/fpw2bdqkmK9s2bJ069YtyRDZ4cOHU7VqVSZOnEinTp04deoU8+fPZ8GCBW9dp5QULFgQY2Nj5s2bh5eXF1evXmXixIlvPO/06dOYmJhQs2bNDK+TEBcuXMDT05PAwEDq16+Pubk5U6dO5ac/D5HHU/tos0OjKnzSdy+VKlXSPdrVURTY8TXc2afdN7GBzzZB3nRMiiayNRnam4V2796Ns7Oz3qtOnTqAdpTIr7/+Su3atSlXrhwHDhxgx44duv+kEyZM4MGDBxQtWjRDm9YrVKjArFmzmDZtGmXKlOH3339PMmy1Vq1aeHl50alTJ/Lly8f06dMBWLFiBT169GD48OEUL16c1q1bc+bMGd039Bo1arB06VLmzZtHhQoV2Lt3r27kSGoGDBjA3r17efLkCW3btqVEiRJ8/vnnWFtb63X4BG1fEBcXFypXrszUqVNp0qQJV69eTbWvRFp88cUXrFu37o2zlE6cODHJI5hKlSqxYcMG1q1bR5kyZfj++++ZMGFC0iboDJAvXz5WrlzJn3/+SalSpZg6dSozZ85843lr166lW7dumJtLB0CRsUJCQmjZsiX58+fn3r17HDx4kJ07d+Ln50fJuh/p8n1avyJNmzZNGogAHJ0JF3/XbqtNoMtacCrznu5AvA8q5W0e5L9noaGh2NjYEBIS8vbLSQvxDhRFoUaNGgwZMoQuXbpkdXUy1LNnzyhRogTnz5+ncOHCWV0dkcvMmzePoUOHcv/+fd0Xk9dazz/G5cfaOUIufd8MG/Nk+sNd/hM2f5643+E3KN0mE2ssMlJaP7+lZUSINFCpVCxZsoT4+PisrkqGu3//PgsWLJBARGSKv//+myZNmiQJRKLjErjuFwaAJtgv+UDk4UnYNiBxv+kECURyKekzIkQalS9fnvLly2d1NTJctWrVdHPYCJHRYmNj9WaWfu3KkxDiNdqG+Tj/ZNbUenEf1nWDhFejwCr3glqDM7GmIitJy4gQQohMU7lyZfbv359kDpt/T3aW3/Q/w85jwmBdV4h6od0v2hg8f4IctAClSB8JRoQQQmSa/v37ExYWxrBhw/Q6gP/jmxiM9G5VP/EEjQY294fA69p9u2LQYQWopSE/N5OfrhBCiExTpEgRFi9eTL9+/Th+/DjdunXD1NSM/Y9cwNgCQyUOry5tEk849CP4/KXdNrWBLuu0/4pcTYIRIYQQmeqzzz7j0qVLrFixgjFjxmBo40h+L+06TjU9nFGrXzXSX98Ox14NRVcZQPsVYJ+2mZJFzibBiBBCiEwTExPDRx99xJEjR/j0009p0qQJJx7HcChae9w69rl249kt2Ppl4onNJoF74/dfYZElJBgRQgiRaWbMmMHx48fZt28fDRo0AODh1qtw+iEAGxZMYVqnSliu7wax4dqTynaAGgNSKFHkRtKBVQghRKZISEhg0aJF9OzZUxeIAJy9rx0lY6CCl7cvEPhrBwh6NbzXoTS0+llGznxgpGVECCFEpnj27BlPnjzB09NTl/YiIhafp9rJzsrkt2HAR/koEnVZe9DEBjqtBuOkK1iL3E2CESGEEJnCxMQE0E4J/tq5By90223y+dO9dOI+7RaDXdH3Vj+RfchjGiGEEJkib9681KpVixUrVugWkHz9iMaGcNrf+Raj159CtYdA8ZZZU1GR5SQYEUIIkWlGjhzJoUOHGDFiBBEREa+CEYWZRouwjgsCQHGtAY2+y9qKiiz1VsHI60W1TE1NqVy5MseOHUs1/++//0758uUxNzfH2dmZ3r178/z587eqsBBCiJzjk08+Yfbs2cyePRunAm5cefySvupdNFX/A0CCSR5U7ZfLDKsfuHQHI+vXr2fIkCGMHTuWCxcuULduXVq2bImvr2+y+Y8fP06PHj3o27cv165d488//+TcuXN8/vnnyeYXQgiRe0RGRuLj44OhoSHxeQtS2uAhowzX6Y6r2y8Fm/xZWEORHaiU1w/y0qh69epUqlSJhQsX6tJKlixJmzZtmDJlSpL8M2fOZOHChdy9e1eXNm/ePKZPn86jR4/SdM3Q0FBsbGwICQnB2to6PdUVQgiRReLj42nRogWnTp1i7NixxLjVosdNL4oa+AOw5r4t7RffwNTUNItrKjJLWj+/09UyEhsbi7e3N82aNdNLb9asGSdPnkz2nFq1avH48WN27dqFoig8ffqUjRs38tFHH6V4nZiYGEJDQ/VeQgghcpatW7dy4MABduzYwZgxY6jru0AXiIRaFaPv7w9Zs2ZNFtdSZAfpCkaCgoJISEjA0dFRL93R0ZGAgIBkz6lVqxa///47nTp1wtjYGCcnJ/LkycO8efNSvM6UKVOwsbHRvVxdXdNTTSGEENnAihUrqF27No0aNSLm0iaaRO8BIApTrHuuo2kLT1asWJHFtRTZwVt1YFX9Z2Y8RVGSpL12/fp1Bg8ezPfff4+3tze7d+/m/v37eHl5pVj+6NGjCQkJ0b3S+jhHCCFE9vH48WMqVKgAoX6o/xqiS9+ZfwjYu1OhQgUeP36cVdUT2Ui6ui/b29ujVquTtIIEBgYmaS15bcqUKdSuXZsRI0YAUK5cOSwsLKhbty6TJk3C2dk5yTkmJia6yXKEEELkTI6Ojty8cR22DsAwVvu4fWdCDdSVPgPg5s2bKX52iA9LulpGjI2NqVy5Mvv27dNL37dvH7Vq1Ur2nMjISAwM9C+jVqsBSGffWSGEEDlI9+7dKRF2HO4dAiBAycvYuD5UK2LHzZs32bZtG927d8/iWorsIN2PaYYNG8bSpUtZvnw5N27cYOjQofj6+uoeu4wePZoePXro8rdq1YrNmzezcOFC7t27x4kTJxg8eDDVqlXDxcUl4+5ECCFEttKpcUVmNjfX7Y+I64+5jR0Hd2ykYcOGeHh40Lt37yysocgu0j3LTKdOnXj+/DkTJkzA39+fMmXKsGvXLgoVKgSAv7+/3pwjvXr1IiwsjPnz5zN8+HDy5MlDo0aNmDZtWsbdhRBCiGwlOjKc8N+6YK/WtoCvjG/GMU05wk/upufOn/D09GT58uVYWlpmcU1FdpDueUaygswzIoQQOcfJkyc5MrEVo6vHA3A32pqP+JloTKhv9oj/dW1MsWLFsriW4n3IlHlGhBBCiNQ8ePCA4Z+15JuqCdoElQFL808kGu2ghENrF0mnVZGEBCNCCCEyzIL5c1nYAowMtI3ucTUGszHAAQAXa2P87lyVic5EEhKMCCGEyDBOt/+ggsOrnXwlOF2wH3EJ2sCkfglHWrRowebNm7OugiJbkmBECCFExnh6jUHlorTbKjW0WcDxB2G6w7WK2uPg4EBkZGQWVVBkVxKMCCGEeHeaBNg2CGP1q/06QyB/ZU7eea7LUr1wXg4ePEipUqWypIoi+5JgRAghxLs7vRD8/gHgxrMENgW5ExwZy1W/EABKOFnx2+L5PHz4MNXlQMSHKd3zjAghhBB6XtyDg5MAUFCxProuE7t2p0mvSyj56gHw7OoxRi77H2PGjKFKlSpZWVuRDUnLiBBCiHSLiopi4cKFVK1ahWOjKkO8tq9IbIWefLdkJ3PmzOF2aOJHTOwj7SiaSZMmZVWVRTYmLSNCCCHSJTg4mKZNm3LhwgXm9KxM3fwaAB6GKLQZu4u/qo3iq6++YmvsYe4+i0CtglPbV2NpIh85InnSMiKEECJdvv76a+7evcuF43sZVCxxFXejtvMJCo2md+/eBIREc/dZBADlXfNIICJSJb8dQggh0iwwMJC1a9cydepUyvqth+hg7YGyHXCp14Np0wzp1q0bG49f0Z1T290+ayorcgxpGRFCCJFm58+fJy4ujq418sPlddpEUxtoPhmAdu3aAbDv0kPdObWKSjAiUifBiBBCiDRTq9WYGoLtqYmJiU0ngqV22tX4eO3iePcjjQEwNTKgUqE877uaIoeRYEQIIUSaVa9ene8bWmIc9kibULAmVOyuO/7HH39gnK8QofHa2c+qF7bDxFCdXFFC6EgwIoQQIs3yxAcxoqa2u6FGpYaP54CB9qPk1KlTfPvtt1T7pJcuf32PfFlQS5HTSAdWIYQQaaMosGs4hirtUN6pRyP5c2sXqlatys2bNzl27Bg1atTApXITHj3QzrxaT4IRkQbSMiKEECJtrm2Be4cBeBis4cdjMVy6dIlNmzZhZGTEunXr+HvfQS480i6Olz+PGUXzWWRhhUVOIcGIEEKIN4sOJWLzYN3u08rDOXrqPDNmzMDExIS7d+9Sr149LjwOIzZB23JSzyMfKpUqq2oschB5TCOEEOKNonePx0KjbfFQPJpTrcs4UKmoXLkyHTt2pEqVKowfPx7njxIDFukvItJKWkaEEEKkLvAGxhdXAKCoTVC1nA7/avFwdXVl4MCBrFmzhsM+TwEwNFBRy90uS6orch4JRoQQQqRMUWDXCAzQPnpR1R0Oed2SZKtVqxaxRlY8eK5dMK9SwbxYmxq9z5qKHEyCESGEECm7vhUeHAPgQQjEVPVKNtvDhw8xLVxJt1+/uDyiEWknwYgQQojkxUbAnrG63cG7Ilm99s8k2eLi4vjll18oVL2FLq1eMQlGRNpJMCKEECJ5x2ZB6BPttnsTrKt2ZNCgQfzyyy9ERGhX5L1x4wbt2rXj6vWbxNsXBcDOwpjSLtZZVWuRA8loGiGEEEm9uA8n5wKQgJoFd1ypUcMdjUbhq6++YsSIEdjY2BAQEICDgwPTlm3g51cL9dbzyIeBgQzpFWknwYgQQoik9v4PEmIBmHUyitlXN/Ds2TMMDQ0ZNmwYLi4uhIeHU6JECT755BNmH7wP3AWgnoes0ivSRx7TCCGE0HfvCNzcCUCIxoweS6/g5+fHkydP+Oqrr/jpp59QqVR8//33dOzYERMTEw7c0A7pVamkv4hIPwlGhBBCJEqIR/P3t7pdm7YzcSzoDoCDgwPTp0/nyy+/ZNKkSURFaYfx+j6P5HZgOKAd0mtnafL+6y1yNAlGhBBCJPrnNwyeXQcg0sYDyndNkuXrr7/mxYsX7N+/H4CDN5/qjjUq4fB+6ilyFQlGhBBCaEUFw8FJut34ZlPAIOnHRKFChQAIDg4G4MDNQN2xJiUdM7WKIneSYEQIIYTW0RkQ9QKAP67EcfR+dLLZTp48CUDRokUJj4nn9L3ngHaVXg9Hy/dTV5GrSDAihBACnt+FM4sBUAxN+f1pMSZMmEBkZKRetujoaMaNG0epUqWoWbMmx249Iy5BAaBxSQdZpVe8FRnaK4QQAvaPA00cAKpaXzGuWTMaNmxI9erVGT58OGXLluX69evMmjWLmzdvsnfvXlQqld4jmsbyiEa8JQlGhBDiQ/fgONzYod22dITaQ7B7/JRPP/2U7du307t3b13Wxo0bc+TIEapVq0aCRuHQq2DE3FhN9cK2WVF7kQtIMCKEEB8yjQb2jEncb/Qd67bspEePHlhaWtK6dWvCwsLYt28f8fHxDBkyhGrVqgFw6XEwzyO0E6PVLWaPqZE6K+5A5AISjAghxIfs8nrwv6TddizLNePydO9emS5durB48WLMzMwACA8Pp0ePHnTo0IEbN27g5uamm+gMoHEJeUQj3p50YBVCiA9VXBQcnJi433wS8+YvwMHBgWXLlukCEQBLS0tWr16NiYkJixYtAuDAjcT+Ig1lfhHxDiQYEUKID9XpBYmr8hZrDkUacODAATp06ICRkVGS7BYWFrRp04YDBw7w+GUkNwPCACjvmod8VjLrqnh7EowIIcSHKCIIjs3WbqsMoOkEABISEjA0TPkJvpGRERqNhr3XEh/RNJFWEfGOJBgRQogP0ZFpEKtt2aBid3AoAUDt2rXZsmULGo0mySmxsbFs376dWrVqsftagC69RRmn91JlkXtJMCKEEB+aoDtwfrl228gCGiaOpvnqq6+4d+8e3377rV5AkpCQwODBgwkKCqJr7/6cf6CdqbWIvQXuDjLrqng3MppGCCE+NAfGgyZeu117MFgltmxUq1aN2bNnM3ToULZs2UKLFi0wMTFh06ZN+Pr6snTpUh7GW6PRTrpK8zJOMuuqeGcSjAghxIfk0Tn9Cc5qDkqS5eOPP2bv3r3s2bOH+fPnA9p1aLZs2ULr1q3pveKsLm/z0vKIRrw7CUaEEOJDoSiw7/vE/Qbfgon+I5arV69Sv359LCwsmDRpEmXKlOHKlSssXLiQPn36sPvAEU7c0S6M52RtSrn8Nu/zDkQupVIURcnqSrxJaGgoNjY2hISEYG1tndXVEUKInMnnb1jbWbtt5w4DToNafwhvzZo1CQ8P5+jRo+TNm1eX/uLFC+rWrYupR22eF/8EgJ41C/HDJ2XeW/VFzpPWz2/pwCqEEB8CTQLsH5+433hckkDk0qVLnD59mkmTJukFIgC2trb88MMP+GoS0+URjcgo8phGCCE+BBf/gGc3tdsFqkLJVkmyXLt2DdAuhpecug0aYXZC25iex9yIarIwnsgg0jIihBC5XVwUHJqcuN90AiQzAsbCwgKAwMDAJMcA9lx8gIGJOQBNSjpiqJaPEJEx5DdJCCFyu7O/QpifdtujBRSqlWy2Ro0aYWVlxeLFi5M9vmyvt25bHtGIjCTBiBBC5GZRwXDsp1c7Kmj8fYpZraysGDp0KDNnzmTOnDlERUVpi4iKYsZPs7gXox15Y26spm4x+0yuuPiQSJ8RIYTIzU7Ohehg7Xb5zuBYOtXs48aNIyAggKFDhzJ27FgKFixIQEAAMbZFcejwAwCNSzpiaqTO5IqLD4m0jAghRG4VFgCnF2q31cbQYHSq2TUaDRMmTOCPP/4AIDIykps3b2JoaIhn///p8n1U1jnTqiw+TBKMCCFEbnV0BsRFarer9IW8hVLNPmrUKCZMmMDAgQPx9fUlPj6eQ4cOUax4Cc76xwJgaWJIg+L5Mrvm4gMjwYgQQuRGL+6B90rttrEl1B2eavbHjx8za9YsfvzxR6ZOnYqrqytqtZoGDRrwv1/+wMBEO9KmaSl5RCMyngQjQgiRGx2akrgYXs1BYJl6a8batWsxNTVl0KCka9Xsu/lct12vsKzQKzKeBCNCCJHbPL0GV/7UbpvZQs2Bbz7l6VMKFCiAlZWVXnp0XAL7rj8FQBMdTjGr+AyvrhASjAghRG5zcBLwatmxusPA9M1rehUoUICHDx/y4sULvfQjt54RHqMNQKLunKWAs8wvIjLeWwUjCxYsoHDhwpiamlK5cmWOHTuWav6YmBjGjh1LoUKFMDExoWjRoixfvvytKiyEECIVj86Bzy7ttpULVP08Tad17doVjUbDjBkz9NL/uuyv265kryFPnjwZVVMhdNI9z8j69esZMmQICxYsoHbt2ixevJiWLVty/fp1ChYsmOw5HTt25OnTpyxbtgx3d3cCAwOJj5emPiGEyFCKAgd+SNyvPxKMzNJ0qoODA+PHj2fs2LE8ffqUQYMGYZvPid1XngAqlOhwZo3qnzn1Fh88laIoSnpOqF69OpUqVWLhwoW6tJIlS9KmTRumTJmSJP/u3bvp3Lkz9+7dw9b27RZVSusSxEII8UG7ewhWt9Fu2xaBgWeTrMybkgsXLrBs2TKOHj3K3bt3iYyMxLx4bfK10c5N0szdkiWf18+kiovcKq2f3+l6TBMbG4u3tzfNmjXTS2/WrBknT55M9pzt27dTpUoVpk+fTv78+fHw8OCbb77RTTOcnJiYGEJDQ/VeQgghUqEocHBi4n7DsWkKRBISEvjiiy+oVKkSW7dupUCBAjg7ayc1K1CrjS5f9/qlMrrGQuik6zFNUFAQCQkJODo66qU7OjoSEBCQ7Dn37t3j+PHjmJqasmXLFoKCghgwYAAvXrxIsd/IlClT+OGHH5I9JoQQIhk+f8OTVwvZOZSG0u3SdNqPP/7IsmXLWLRoEX379sXQ0BBFUfh941bGnlWhAvJZmVCziF3m1V188N6qA6vqP0tPK4qSJO01jUaDSqXi999/p1q1anh6ejJr1ixWrlyZYuvI6NGjCQkJ0b0ePXr0NtUUQogPg0YDh35M3G80Fgze/Oc9KiqKn3/+mcGDB9O/f38MDbXfT1UqFQYFK6N61bLS2N0GQ7UMvhSZJ12/Xfb29qjV6iStIIGBgUlaS15zdnYmf/782NjY6NJKliyJoig8fvw42XNMTEywtrbWewkhhEjB9S3w9Kp226USFPdM02lnzpzhxYsX9O7dO8mxrRef6LatX/pkSDWFSEm6ghFjY2MqV67Mvn379NL37dtHrVq1kj2ndu3a+Pn5ER4erku7desWBgYGFChQ4C2qLIQQQichHg5NTtxv9D9IoaX6v2JiYgD0viwCPH4Zydn72vlG4oJ8sSUiY+oqRArS3e42bNgwli5dyvLly7lx4wZDhw7F19cXLy8vQPuIpUePHrr8Xbt2xc7Ojt69e3P9+nWOHj3KiBEj6NOnD2ZmaRtyJoQQIgWX18HzO9rtQrWhaKM0n1q2bFnUajV//fWXXvq2i3667fBrh6hcuVKGVFWIlKR7npFOnTrx/PlzJkyYgL+/P2XKlGHXrl0UKqRdDdLf3x9fX19dfktLS/bt28dXX31FlSpVsLOzo2PHjkyaNCnj7kIIIT5E8bFwZFrifqPv0twqAuDi4kLbtm2ZNGkSLVq0oEiRIiiKwpYLiY9oihq+pEqVKhlZayGSSPc8I1lB5hkRQohknFsGfw3TbhdtDN03p7uIgIAA6tWrh7+/P927d8eldHWWPrIHICHAh72jWlKyZMmMrLX4gGTKPCNCCCGyibhoODozcb/h2LcqJigoiMqVKxMXF8fChQuZtSlxeY9RnRpJICLei3Q/phFCCJENeK+AsFd9O4p7QoHK6S5i3759tG7dGgcHB0aNGkU+B0fm3LcnHkCTQMeaxTK0ykKkRFpGhBAip4mNhGOzEvcbjkl3EREREXTq1IkGDRrg4+PDDz/8QJkm7Yk3tAAg6t55li+al1E1FiJVEowIIUROc+5XiAjUbpf6BJzKpruIdevWERwczMKFCzE1NQXgT+/EuZ+q2WtYsGABGo0mQ6osRGokGBFCiJwkJgyOz3m1o4IG6W8VATh//jxly5bFzc0NgODIWPZdewqAvaUxvZpV4dGjRzx79uzd6yzEG0gwIoQQOcmZRRClnZCMsu3BocRbFWNsbEx4eDivB1Ruu+hHbIK2FaRNhfxERoTp8gmR2SQYEUKInCI6BE6+6sehMoD63751UZ6enrqFTAE2nE9cA6x95QKsXLmSGjVqkDdv3neqshBpIaNphBAipzi1QBuQAJTvAvbub11U06ZNqVChAt27d2fWig1c8wsFoKyLFSt/nszBgwfZtGlTRtRaiDeSlhEhhMgJIl/A6QXabQNDqDfinYozMDBg+/btmJub8/mPy3Tp5zbMZfr06UyfPp127dq90zWESCtpGRFCiJzg1HyI0bZeUKEb2BZ+5yLPnTtHvAIWpeoDoImLweDRBdavX0+HDh3euXwh0kpaRoQQIruLeA6nF2m31cbv3CoCsHbtWj799FMcKjRGba5dtbeKkyH5HWzp2bMnZ86ceedrCJFWEowIIUR2d/JniIvQblfqAXlc36m4mJgYhgwZQqdOnXBv2l2XPqxNLY4ePUqpUqUYPnz4O11DiPSQYEQIIbKz8EA4+6t2W20CdYa9c5G7du0iMDCQ/sNGc+S2dh6R/HnMqFXUDjMzM7799ltOnDiBj4/PO19LiLSQYEQIIbKzEz9DXKR2u0pvsMn/zkU+evQIU1NTvF+a8nrd9s5VXTEwUAFQqVIlXT4h3gcJRoQQIrsKC4BzS7XbhqZQZ2iGFJsvXz6iY+P448wDANQGKjpVTXz0c+vWLV0+Id4HCUaEECK7Oj4H4qO121X6gpVThhTbqlUrbMvW50VkPABNSzriYK1dn0aj0fDTTz9RtmxZypUrlyHXE+JNJBgRQojsKNQfzi/XbhuaQZ0hGVa0paUlpVv31+03djNBURQuXrzIp59+yoEDB5g6dSoqlSrDrilEamSeESGEyI6Oz4KEGO121b5g6fDORcbHx7Nz5052HTmLr0ltABJCntKxXivUagMSEhLInz8/GzduxNPT852vJ0RaSTAihBDZTcgT8F6p3TYyh9pD3rnIa9eu0bp1a+7du0fhNkOhuDbd4P5Jpk+fhoWFBYUKFaJ58+YYGspHg3i/5DdOCCGym2M/QUKsdrvaF2D5bh1Jnz9/TpMmTXBwcODMOW8G7H5OUHgshgZgGezDzJk7uHLlCg4O7976IsTbkD4jQgiRnQQ/gn9WabeNLaHW4HcucunSpbx8+ZLdu3cTaOxEULg20Glexpk92zYSFhbGr7/++s7XEeJtSTAihBDZybGfQBOn3a7eHyzs3rnITZs20bZtW5ydnfnt5ANderfqBXF0dKRdu3ayQq/IUhKMCCFEdhHsCxfWaLeNraDmoAwpNjw8HEdHR675hXDuwUsAPBwtqVlEG+g4OzsTHh6eIdcS4m1InxEhhMgujs5IbBWp4QXmthlSbMmSJTl06BDGdXrr0nrUdEOlUqEoCvv376dkyZIZci0h3oa0jAghRHbw8gFc/EO7bWINNQdmWNFeXl5cvf2Azd7a6d2tTA1pW1E7rfzq1au5ePEi/fv3T60IITKVtIwIIUR2cHQGaLQzolJjAJjlzbCimzRpQqPP/8dt5dXaM3mi2btrBxs2bGD9+vX06dOHli1bZtj1hEgvaRkRQois9uIeXFyr3TaxgRpfpuk0b29vevXqRdGiRXF3d+eLL77g8uXLSfIlaBTCnbWL36Fo+GN8f9q1a4e3tzdz587l119/ldlWRZaSlhEhhMhqR2eCkqDdrjkQzPK88ZSFCxcycOBAChUqRLt27dBoNGzcuJHly5ezfPlyevToocu7/8ZT/EO0a9w0KunEjzcuoCgKzs7OEoSIbEGCESGEyErP78KlV60ipjbajqtv4O3tzcCBA/nqq6+YNWsWarUagBkzZuDl5UWfPn2oVq0aJUqUAGDFiQe6c3vVcsPZWVbjFdmLPKYRQoisdGQ6KBrtds2vtAHJG8yfP59ChQrpBSIAhoaG/PLLL9jZ2bFw4UIArvmFcOb+CwCK5rOgbjH7jL8HId6RBCNCCJFVgm7DlQ3abbO82knO0uD48eO0a9dOLxB5zcTEhNatW3P06FEAlh67rzvWq5abPJYR2ZIEI0IIkVX+3SpS6yswtU7TaSqVioSEhBSPazQaVCoV/iFR7LjkB0BecyPaV3Z95yoLkRkkGBFCiKzwzAeu/KndNrPVLoiXRg0bNmTjxo3ExcUlORYVFcXWrVtp1KgRK088IF6jAPBZjUKYGSdtSREiO5BgRAghssKRaYA2UKD2YDCxSvOpgwYNIiAggP79+xMTE6NLj4qKolevXoSHh9Oj7xf8ccYXAGO1AT1qumVg5YXIWDKaRggh3rfAG3B1s3bb3B6q9kvX6WXLlmXlypX07t2bnTt3UrduXVQqFQcPHiQiIoK1a9dyNsiQsBjtJGptK+Ynn5VJRt+FEBlGWkaEEOJ902sV+RpMLNNdRJs2bfjss88ICQlh8+bNulV3x48fT6vWn+gN5/28buEMqLQQmUdaRoQQ4n16eg2ubdFuW+SDqn3TXURUVBTNmjXjypUrDB06lI8//piwsDBWrFjBmDFjuBpqwhOlOAANi+ejmGPaHwEJkRUkGBFCiPfp8JTE7TpDwdgi3UUsWLCA8+fPc/z4capVq6ZLb9myJZMnT2HeLTUmTtq0fvWKvGuNhch08phGCCHeF//LcGOHdtvSEar0eatilixZQseOHfUCkdeqfdITEyd3AMrkt6ZmEbu3rq4Q74sEI0II8b4cnpq4XXc4GJmluwhFUbh9+zZ16tRJ9viS4w912wMbuMskZyJHkGBECCHeB78L4POXdtvKBSr1fKtiVCoVNjY2PHr0KMmxcw9ecPbV1O9mcSE0L+301tUV4n2SYEQIId6HQ//qK1J3GBiZvnVRnTp1YsWKFYSFhemlLzh0R7f9aUkrDAykVUTkDBKMCCFEZnt8Hm7v0W5bF4BKPd6puOHDhxMeHk7Lli25cOECAJcfveCQzzMADKKC+V+PFu90DSHeJxlNI4QQme3Qj4nb9YaDYfonIAsMDGTLli28ePGCwoULs2PHDrp3706lSpXInz8/Ss0+GBWtDsCIjytgamyUUbUXItNJMCKEEJnp4Sm4e1C7nacgVPgsXadrNBrGjBnD7Nmz0Wg0WFtb8+LFC+zs7Jg7dy7m5uYcPH+N7fHlAbC3NKF3gxIZfRdCZCp5TCOEEJnp360i9UeBoXG6Th8zZgzTp09n7NixBAQE8Pz5c3x8fGjcuDGfffYZRkZGqMt+BGj7h3xetzCmRrIgnshZVIqiKFldiTcJDQ3FxsaGkJAQrK3TtsS2EEJkuftH4bdW2m3bIjDwHKjT3iAdGBiIq6srY8eO5fvvv9c7ptFoaNq0KUExasLqfo1GgbzmRhwb1QhLE2n0FtlDWj+/pWVECCEyg6LAocmJ+/W/TVcgArBlyxY0Gg0DBw5McszAwICvv/6aJ3nLoXn1lbJfvSISiIgcSYIRIYTIDHcPgu8p7bZdMSjbPt1FvHjxAmtra+zskp9F1dC2ABYl6wHaVpGeNd3etrZCZCkJRoQQIqMpin5fkQbfgkH6+3EULlyYFy9ecOvWrWSP/3L4LqpX5X5RrygW0ioicigJRoQQIqPd2g1PvLXbDqWhdLu3KqZNmzbY2dnx3XffodFo9I6dv/2Ea2HaidNsLYzpUbPQO1VZiKwkYbQQQmQkjUa/VaThGDB4u+99pqam/PLLL3Tp0oWgoCC+/vprChcuzOnTp5lyOABcKwHQv14RaRUROZr89gohREa6uQMCrmi3nStAiY/eqbgOHTpw8+ZN5syZwyeffAKAsUNhnHvNBcDe0pju0ioicjgJRoQQIqNoEvRH0DT6H7zDqrkxMTF06NCBHTt2UK5cORo0aICPjw9BpTroyvWqXxRzY/lTLnI2+Q0WQoiMcnUzPLup3S5QDdybvFNxo0aNYu/evWzbto1WrVqhUqk49+AFHRZpR+kYx0fwWQ1pFRE5n3RgFUKIjJAQB4f/3Soy9p1aRUJCQvj111/59ttvad26NSqVCkVRmPb3TV2egP3LCHrq/y61FiJbkGBECCEywsU/4MU97bZbXSjS4J2KO3XqFJGRkXz2WeJaNgdvBnL+4UsACtuZEXZ5P4cOHXqn6wiRHbxVMLJgwQIKFy6MqakplStX5tixY2k678SJExgaGlKhQoW3uawQQmRP8TFwZHrifuPvU86bRgkJCQCYmGhX+E3QKEzf7aM7PrSxOyga4uPj3/laQmS1dAcj69evZ8iQIYwdO5YLFy5Qt25dWrZsia+vb6rnhYSE0KNHDxo3bvzWlRVCiGzJeyWEPtZuF2sOrtXeucjKlStjaGjI1q1bAdh+6Qk+T8MAKF/AhvCbJwCoWbPmO19LiKyW7oXyqlevTqVKlVi4cKEurWTJkrRp04YpU6akeF7nzp0pVqwYarWarVu3cvHixTRfUxbKE0JkW7ER8HMFiAjU7vc/Cs7lM6Torl27snfvXvbsP8jg3UE8CY4C4KePCzG0qyfFixdn3759GXItITJDpiyUFxsbi7e3N82aNdNLb9asGSdPnkzxvBUrVnD37l3GjRuXpuvExMQQGhqq9xJCiGzp7JLEQKTUJ+8ciERHRxMcHIyiKMyfP58CBQrQZOCPukAkb0wAnzWtirGxMStWrHjX2guRLaQrGAkKCiIhIQFHR0e9dEdHRwICApI95/bt23z77bf8/vvvGBqmbSTxlClTsLGx0b1cXV3TU00hhHg/ooLh+BzttsoAGo5966KOHj2Kp6cn5ubm5M2bl4IFC7Jw4ULWb9tFvnraTqyKRkPUyd8ZN24c586do0CBAu9+D0JkA281z4jqP8PVFEVJkgbaDlhdu3blhx9+wMPDI83ljx49mmHDhun2Q0NDJSARQmQ/J+dBdLB2u1xnyFf8rYr5448/6N69O+XLl2f+/PnY29uzb98+Jk6cyNrbCrFOFQHoWsONKdMPZlDlhcg+0hWM2Nvbo1ark7SCBAYGJmktAQgLC+P8+fNcuHCBQYMGAaDRaFAUBUNDQ/bu3UujRo2SnGdiYqLrQS6EENlSeCCcftV3zsBIuzLvW3j+/Dmff/453bp1Y8WKFajV2lV4O3bsSNP2PRhx4CUqwMJYzbCmaf9SJ0ROkq7HNMbGxlSuXDlJh6l9+/ZRq1atJPmtra25cuUKFy9e1L28vLwoXrw4Fy9epHr16u9WeyGEyCrHfoK4CO12ld6Q9+1mQl21ahUJCQnMnDlTF4i8tsvPBJWBNs2rQVHyWcmXNJE7pfsxzbBhw+jevTtVqlShZs2aLFmyBF9fX7y8vADtI5YnT56watUqDAwMKFOmjN75Dg4OmJqaJkkXQogcI9gXzi/XbhuaQd1v3rqoK1euULFiRRwcHPTSj9x6xmGfZwDEhwbSvox8eRO5V7qDkU6dOvH8+XMmTJiAv78/ZcqUYdeuXRQqpP1W4O/v/8Y5R4QQIkc7PA0SYrXbNbzAKulj6rQyNzcnKChIr+9dbLyGH3Zc0+UJPrIKG8tO71RlIbKzdM8zkhVknhEhRLbxzAcW1ABFAyY2MOQSmOV96+L27dtHs2bNOHDggK4P3eIjd5nyag0aw2Bfij7Yxp7duzOk+kK8T5kyz4gQQnzwDkzQBiIAdb5+p0AEoHHjxlSrVo1u3bpx6NAhAkKimHvgtvagovB42yy+HTXqHSstRPYmwYgQQqTV4/Nwc6d229IRqnu98ZSwsDC2bNnC6tWr8fb25r+N0QYGBmzfvp2CBQvSqFEjanpNIyJWuy5N5JW9LJ0xjoYNG2b4rQiRnbzVPCNCCPHBURTYPz5xv/4oMLZIMbtGo+GHH35g9uzZhIWF6dIrV67M0qVL9RYMdXR05NSpUyzatI/p3tqF70xVCRxYMoaiBZwy+k6EyHakZUQIIdLi7gF48GqFctsiUKlHqtmHDRvGxIkT+fLLL7l//z6RkZH89ddfaDQaGjRowM2bN/XyaxTY6Weu2x/7SXkJRMQHQ4IRIYR4E40G9v+QuN/of6A2SjH7vXv3mDt3LtOnT2fatGm4ublhZmaGp6cnhw8fxtbWlgkTJuid89uph1z3167DVcrZmq7VCmbKrQiRHUkwIoQQb3JtMwRc1m47lYNSbVPN/vvvv2NlZcWAAQOSHLO2tmbgwIFs3LiRiAjtpGl+wVH8tNcHAJUKJrYpg9og6RIbQuRWEowIIURq4mO0I2heazIeDFL/0xkQEEDhwoUxNzdP9nipUqWIi4vjxYsXAIzbfo3IV51Wu1YrSOVC7zZCR4icRoIRIYRIzfkVEPxQu12kAbg3fuMpLi4u3Lt3j/Dw8GSPX7lyBWNjY+zs7NhzLYB9158CYG9pwsgWJTKq5kLkGBKMCCFESqJD4ej0xP0m49N02meffUZERARz585Ncuzly5f88ssvdOrUCY3amPHbE2da/b5VKWzMUu6LIkRuJcGIEEKk5ORciHyu3S7THlwqpum0QoUKMWLECMaOHcvgwYO5fv06L168YOPGjdSpU4fw8HC+//57Zu29hX9INAD1PPLRqpxzZt2JENmaBCNCCJGcsAA49Yt228BIO4Imje7cuQNA2bJlWbJkCaVLl8bOzo4OHTpgZ2fH0aNHCTGyY8XJ+wCYGBow6ZMyurVphPjQSDAihBDJOfQjxEVqt6v2BdvCbzxFURTGjRtHsWLF+PXXX3FwcKBAgQIAlClThjNnznD06FGKepRg5MZLvJ6MdXgzDwraJd/ZVYgPgQQjQgjxX0+vw4U12m0Ta6g3Mk2nLV26lAkTJjBhwgQeP37M/v37uX37Nn/99Re+vr5MnToVgLkHbnP3mXZYb3nXPPStUyRTbkOInEJW7RVCiP9a0x7u7NNuNxkPdYa+8RSNRoOHhwdVqlRh3bp1SY6vWrWKnj17sv3YRYbsekKCRsFIreKvwXXxcLTK4BsQInuQVXuFEOJt3DucGIhYF0jTYngAN2/e5O7du/Tt2zfZ4506dcLc0orxf98hQaP9Dji4UTEJRIRAghEhhEik0cDef3VUbfwdGJml6dToaO2omLx5k5+wzMTEhDy1u/A8wRTQTvnu1aDou9VXiFxCghEhhHjt8noIuKLddioHZTum+VR3d3fMzc35+++/kz3+++6TGJb7CABDAxXT25fDSC1/goUACUaEEEIrNkJ/2vdmk9447fu/WVtb061bN2bPns3169f1jgW+CGHc33dRGagB+LpxMcrkt8mQaguRGxhmdQWEECJbODkfwvy02x4toUj9dBcxdepUzpw5Q7Vq1ejevTu1atXi4cOH/HohDFWxegBULJiHL+XxjBB6pGVECCFC/eHEHO22gSE0m6g75OPjg5eXFw4ODlhYWFC1alWWLVtGfHx8kmJsbW05evQow4cPZ/v27fTo0YOZv+/SBSJmRmpmdayAoTyeEUKP/I8QQoiDExMnOKvSF+yLAXDo0CEqVarEjh076Nu3L5MmTcLR0ZF+/frRrl074uLikhRlY2PDDz/8wJMnT/B7HkqxzxIf/Yz5qCSF7S3eyy0JkZPIYxohxIfN7yJc/EO7bWoDDb4FIDIyko4dO1KrVi22bduGubl2htShQ4fy999/07p1a2bPns3IkclPiKYoCmO23eBZWAwA9T3y8Vn1gpl+O0LkRNIyIoT4cCkK7BkLvJr7sd5IMLcFYMOGDTx//pzFixfrApHXWrZsyWeffcaCBQvQaDTJFr3ixAMO+TwDwM7CmBkdysnaM0KkQIIRIcSH6/o2eHhcu523MFTrpzt07tw5ypQpQ5EiyU/V/sknn/Dw4UOePXuW5NjVJyFM/fumbv+njuVxsDLN2LoLkYtIMCKE+DDFRcHe7xL3m08GQxPdrqGhIZGRkaS0YkZEhHZtGSMjI/30mHgGr71AbIK2xaRf3cI0KO6QwZUXIneRYEQI8WE6NR9CfLXbRRpC8ZZ6h1u0aMHdu3c5depUsqevWrWKqlWrYmtrq0tTFIXvtl3lXpA2UCmb34YRzUtkTv2FyEUkGBFCfHhC/eDYLO22Sg0tpsB/+nM0b96cMmXK0L17d65du6ZLj42NZdy4cezdu5dvvvlG75x15x6x+Z8nAJgbq5nbpSLGhvJnVog3kdE0QogPz/7xiUN5q/YlwqIg/xw7hqIoVKhQAWtrawwMDNixYwfNmjWjTJky1K1bF0dHR44ePUpgYCCTJk2iY8fE6eKvPglh3PbEoGXqp+VkGK8QaSQhuxDiw+J7WrsGDaCY5uX7Q9G4uLhQr1496tevj4uLC4MGDSIyMhI3NzcuX77MqlWrsLW1JTg4mM6dO3Pt2jXGjh2rKzIkMo4vf/cmNl7bT6RHzUK0Lu+SJbcnRE6kUlLqnZWNhIaGYmNjQ0hICNbW1lldHSFETqVJgCX1dYvhLfR1Z8jvVxk2bBhdunRBrVbz559/MnPmTKpWrcqePXswNjZOtUhFUei3ypv9N54CUN41Dxv618DEUJ3ptyNEdpfWz295TCOE+HB4r9QFIiFmBRm08h+279jJRx99pMtSunRpGjduTL169Vi7di09e/ZMtcgFh+/qApE85kb80rWiBCJCpJM8phFCfBgiX2infX9l8kVbqlStpheIvFa3bl2aNWvGsmXLUi3y4M2nzNzrA2j7v87uVIECec1TPUcIkZQEI0KID8PBiRD1UrtdtiO7rwdTtWrVFLNXrVqVhw8fpnj8TmA4X6+9yOsH3cObetBQ5hMR4q1IMCKEyP38LsD5FdptY0toOgF7e3vu3LmT4im3b9/G3t4+2WOh0XF8sfo8YTHalXs9yzoxsKF7hldbiA+FBCNCiNxNkwA7h6Fbf6b+SLB2plu3buzdu5fLly8nOeXu3bts2bKFbt26JTmWoFEYuu4i955pJzYr4WTFjPblZd0ZId6BBCNCiNzNeyX4/aPdzlcCagwAoEuXLpQtW5ZmzZqxdu1aYmNjiY+PZ/PmzTRq1IhChQrRt2/fJMVN2XWDAzcDAW2H1SXdq2BhImMBhHgX8j9ICJF7hT+DAxMS9z/6ifCoGNavX8X169dp0qQJp0+fpmvXrhgaGqJSqYiLi6Nu3br88ccf2NjY6BX3xxlflh6/D4DaQMUvXStR0E46rArxriQYEULkXvvHQXSwdrtcZzZ5P6Vv3wKEhYVRrFgxgoKCeP78OfXr16dVq1aYmJhQp04dKlSokKSo47eD+G7bVd3+xE/KUNs9+T4lQoj0kWBECJE7PTwJF3/XbpvYcMrak04tPqVdu3bMmDGDQoUKERcXx6ZNm+jfvz/m5ubs2rUr2aLuBIbx5e/eJGi0/U761S1M1+oF39edCJHrSTAihMh94mNgx9eJ+42/Y/x3Cyhfvjx//PEHhobaP31GRkZ07twZQ0NDOnTowLlz55IM930aGk3P5ecIi9aOnGlS0pFvW5Z8b7cixIdAOrAKIXKfE3Mh6JZ2O39lQop9yt69e/Hy8tIFIv/Wtm1bnJyc2Lhxo156WHQcvVac40lwFAClnK35uXMF1AYyckaIjCTBiBAid3l+F47O0G6r1PDxHCKiogFwcUl+8Tq1Wo2TkxNhYWG6tNh4DV5rvLnhHwpAgbxmrOxdVUbOCJEJJBgRQuQeigI7h0BCjHa/5gBwLke+fPnImzcvR44cSfa0wMBArl27RvHixQHQaBRGbrzEiTvPAe0Q3t/6VMPB2vR93IUQHxwJRoQQuceldXD/qHbbpiA0GA1o+4b06dOHxYsXc+PGDb1TNBoNY8aMQa1W0717dxRFYcLO62y96AeAqZEBy3pWpWg+y/d6K0J8SKS9UQiRO4Q/gz2jE/c/mgnGFrrdsWPHsmfPHmrVqsWXX35J48aNefr0KYsWLeLYsWMsX74cW1tbZu27xcqTDwAwUMG8LpWoXCjve74ZIT4sEowIIXKHv0cmLoRX5lPwaA6Aoijs27ePxYsXExUVhaGhIbNmzWLKlCkA1K5dm7/++gtPT0+WHrvH3AO3dUVOb1+epqUc3/utCPGhkcc0Qoic7+ZfcG2zdtvMFlpMA7SByMCBA2nevDl37tyhTZs2eHp6olaryZs3L3v27OH48eN4enqy7qwvk/5KfIQzrlUp2lcukBV3I8QHR6UorxfAzr5CQ0OxsbEhJCQEa2vrrK6OECI7iQ6BX6pDmL92v+0SKN8JgGXLlvH555+zZMkSPv/8c91idkFBQXh6euLn58fdu3fZefUZ32y8xOu/hsOaejC4cbGsuBshcpW0fn5LMCKEyNm2D4Z/ftNuuzeFbn+CSoWiKJQrVw53d3e2bNmS5LQbN25QqlQpRvyygT8fmesCkc/rFGbsRyVlFV4hMkBaP7/lMY0QIue6cyAxEDG2hI9nw6sgIjg4mKtXr9KxY8dkTy1ZsiQlmnVjw0MzXSDSs2YhCUSEyAISjAghcqboEG2ryGtNJ0AeV93umwKKHZf8iCrfURe8fFajIONbl5ZARIgsIMGIECJn2jMWQh9rt4s0gCp99A7b2NhQrlw51q9fn+TUTd6P+XrdBTBQA9ClWkEmtC4jgYgQWUSCESFEznN7P1xYrd02toLW83QtHK+pVCqGDBnCtm3bWLRoEa+7x/1xxpfhf17i1QK8dKicnx/blMFA1psRIsvIPCNCiBxBo9Hw4MEDVNEhuP09CF3o0HwS5CmY7Dm9evXin3/+4csvv2T+/PkUadGHy4bFdcc93S2Y9ml5CUSEyGLSMiKEyNY0Gg1z586lWLFiFC1alFPf1Ub1ahivpkhDqNQzxXNVKhVz585l//792NTqrBeIfFbFiV/61pdARIhsQIIRIUS2pSgK/fv3Z8iQIdSuXZsLq8bStawRAC+jFLx2xZOg0byhDDga7sgT24q6tMGNizHx00rSR0SIbEKCESFEtnXo0CGWLl3KsmXLWDVvChWerNIdu1tyIL+u/4vNmzeneH5svIYh6y/y26mHurSxniUZ1tRDAhEhspG3CkYWLFhA4cKFMTU1pXLlyhw7dizFvJs3b6Zp06bky5cPa2tratasyZ49e966wkKID8eSJUsoXbo0vXr2gG0DtMN5Acp8SpVeU6hbty5LlixJ9tzwmHg+X3We7Ze0q++qDVTM7FCefvWKvK/qCyHSKN3ByPr16xkyZAhjx47lwoUL1K1bl5YtW+Lr65ts/qNHj9K0aVN27dqFt7c3DRs2pFWrVly4cOGdKy+EyN18fHyoV68eqtML4N5hbaKVC3jOBKBevXr4+PgkOe9paDQdF53i6K1nAJgYGrD4s8qy1owQ2VS6p4OvXr06lSpVYuHChbq0kiVL0qZNG90qmG9SunRpOnXqxPfff5/s8ZiYGGJiYnT7oaGhuLq6ynTwQnxg6tevTyVnNbNLXQZNnDax+1Yo2hCAPn36cPbsWa5evao7xycgjN4rzuIXEg2AjZkRv/aoQrXCtu+7+kJ88DJlOvjY2Fi8vb1p1qyZXnqzZs04efJkmsrQaDSEhYVha5vyH4YpU6ZgY2Oje7m6uqaYVwiRe3Vt3xqvfOcSA5Fag3WBSGBgIBs2bNCb7v347SDaLzypC0QK5DVj05e1JBARIptLVzASFBREQkICjo6OeumOjo4EBASkqYyffvqJiIiIFNeLABg9ejQhISG616NHj9JTTSFELtHbyYfidto/U2FW7iiN/gfAuXPnaNasGZaWlvTv3x+AVace0HPFWcJi4gEoV8CGLQNq4+5gmTWVF0Kk2VtNevbfXuiKoqSpZ/ratWsZP34827Ztw8HBIcV8JiYmmJiYvE3VhBDZREJCAnv37uXGjRtYWFjQunVrnJ2d017AlY0YX9NO5R4Zr6LS5ItEzS6KWq3G19cXd3d39u/fj619Pv639QprTif2W2tS0pG5XSpgbizzOgqRE6Trf6q9vT1qtTpJK0hgYGCS1pL/Wr9+PX379uXPP/+kSZMm6a+pECLHOHLkCL169eLBgwdYWFgQHR3NoEGD6NevH3PmzMHY2Dj1Ap7fhR1f63bNPp3PwgZOHDx4EEVRqFOnDi1atCAkOoGey89y8u5zXV6v+kUZ0bw4apnMTIgcI13BiLGxMZUrV2bfvn20bdtWl75v3z4++eSTFM9bu3Ytffr0Ye3atXz00UdvX1shRLb3zz//0KJFC2rUqMGGDRuoWrUqwcHBLF26lLFjxxIVFcWKFStSLiAuCjb0hNhw7X65zqgqdKOJSqX3RebqkxD6r/bmSXAUAMZqA6a0K8unMmJGiBwn3W2Yw4YNo3v37lSpUoWaNWuyZMkSfH198fLyArT9PZ48ecKqVdrJidauXUuPHj34+eefqVGjhq5VxczMDBsbmwy8FSFEdjBhwgQKFy7M33//jampKQB58uThm2++wdramv79+zNy5EhKliyZfAG7R8PTK9ptew/46Kcki+Bt8n7MmC1XiInXzr5qb2nM4u6VqVxIOqoKkROle56RTp06MWfOHCZMmECFChU4evQou3btolChQgD4+/vrzTmyePFi4uPjGThwIM7OzrrX119/ndIlhBA5VFhYGDt27GDAgAG6QOTfevbsia2tLWvXrk2+gMt/gverVhNDM+jwG5gkdkCNiU9g3LarDP/zki4QqeCah51f1ZVARIgc7K16dw0YMIABAwYke2zlypV6+4cPH36bSwghcqDg4GA0Gg3u7u7JHjcxMcHV1ZXnz58nPfj0GuwYnLj/0UxwLKXbffQikoF//MPlxyG6tC7VCjK+dSlMDNUZdg9CiPdP1qYRQmQYe3t7zM3NOXv2bLLHg4ODuXXrFm5ubvoHooJh/WcQF6ndr/AZVOimO7z3WgAfzT2mC0SMDbX9Q6a0KyuBiBC5gAQjQogMY2ZmRrdu3ViwYAF+fn5Jjk+ZMoW4uDi6d++emKjRwBYveHFPu+9cXtsqolIRHZfA+O3X+GK1N6HR2vlD3OzM2fxlLbpUK/g+bkkI8R5IMCKEyFDjxo3D0NCQGjVqsGDBAnx8fDhy5AidO3dm+vTpTJo0CScnp8QTjs2EW39rt83yQsfVYGTGradhtPnlBCtPPtBl/aisMzu+qkOZ/NL5XYjcRGYEEkKkS3h4OH/88Qfbtm0jKiqKcuXK0b9/f93omPz583Py5EmGDh3K4MGDSUhIAKBw4cIsW7aMPn36JBZ2Ywcc+vHVjgo+XYaSpyBrTj1g0l83dJ1UjQ0N+N9HJeleo1CaJlgUQuQs6V4oLyukdaEdIUTmunnzJs2bN+fx48c0atQIW1tbDh06xLNnz5g5cybDhw/Xy+/v78/t27extLSkQoUKGBj8qzH26TVY2hTiIrT7jccRUG4AIzdd1q22C1Dc0Yq5XSpS3MnqfdyiECIDpfXzW1pGhBBpEhsby0cffYSlpSW3bt2iaNGigHaV7e+//55vvvmGEiVK6E1s+HoofxIRz2FtZ10gopTtwHbLjnw3+4iubwhAz5qFGO1ZElMj6aQqRG4mLSNCiDRZv349nTt35vLly5QtW1bvmKIo1K1bF2NjYw4ePJh6QfGxsLotPDwOQJxjBb6xnMK2ay91WRysTJj2aTkalkh5DSshRPYnLSNCiAy1Z88eKlSokCQQAe3imd27d8fLy4uYmJiUF7pUFNj+lS4QiTaxp3WgF7ceJgYircq7MPGT0uQxf8P6NUKIXEOCESFEmsTHx2NmZpbicXNzcwBdh9VkHZkGl9cBEKsyplPoYG4p2m9LthbG/NC6NK3Ku2RcpYUQOYIM7RVCpEm1atU4d+5csvOHAGzZsoXSpUvrgpIkLq2Hw1MA0CgqvooZwCVFO1Nrq/Iu7BtaTwIRIT5QEowIIdKke/fumJub88UXXxAdHa13bMOGDWzdupVBgwYlf/K9w2i2DdTtTo7vyh5NNRysTPi1RxXmdamInWUKj3aEELmePKYRQqSJjY0N69ato127dri7u9O9e3fs7Oz4+++/OXjwIN26deOLL75Ict7Lu+cw+70Lppo4ANbEN2aZxpMeNQvxTfPiWJsave9bEUJkMxKMCCHeKC4ujl9++YWFCxcSHR3NkydPmDVrFkZGRlSsWJE1a9bQpUsXvXlE4hM0bDt4nAYnPsMU7Zoz+xIqs85+MFs+rUAF1zxZdDdCiOxGghEhRKri4uJo27Yte/bsoUOHDvzvf/8jODiYFStWcPHiRT7//HO6deumd87RW89YsOMkU0NGYGegXdzuH6U4fk3ms7VOCQzV8oRYCJFI5hkRQqRq3rx5DB06lL/++ovmzZvr0jUaDf369WP16tX4+vri5OTEncAwJv11g4s+91hnPIkSBo8ACDApjGHf3dg7OKV0GSFELpTWz28JRoQQqSpZsiTlypVj/fr1SY69fPmS/PnzM3T0OJQynqw/9wgzTQS/G0+mvIF2Fd4Yi/yYfLEPbPK/76oLIbJYWj+/pa1UCJGi2NhYbt68SYsWLZI9rjazxKP9N6wNK8EfZ3wx1kSx3HiGLhBRLB0x6bNDAhEhRKqkz4gQIkWGhoYYGhry8uVLvfTI2Hh+O/mQxUfvEuxSHQAzollp8hPVVD7aTGa2qHpsA7ui77vaQogcRlpGhBApMjAw4OOPP2bFihUkJCQQHZfA8uP3qTf9ENN23yQ4Ujtc11IVxS77+VRXXdOeaGID3beAQ8ksrL0QIqeQlhEhRKpGjBhB/UZNaDpwMqH5qxMUEZd4UNFg6XeWc3XPYub3jzbtdSDiUiFL6iuEyHkkGBFCpCg8Jp4rcY4U/2Y99+IN4F+BSMSNo7gFneLYZwaY+l3QJprYQI8tkL9yFtVYCJETSTAihEjiWVgMK0/eZ/Wph4RGx/PvJ7oRPicIObGO/KogdgxwwvTpU+0BUxvovhXyV8qSOgshci4JRoQQOveehbP8xH3+PP+YmHiNLl1RNJgEXGV4yzJ0GfMtT283x2rLZzgprwIRczv4bLM8mhFCvBUJRoT4wCmKwpn7L1h67D4Hbj7l3zMPGalVWL/wIfTMJryP/I2VlRU888H66EAwDAXgcagG446rcZBARAjxliQYEeIDFR2XwLaLT/jt5EOu+4fqHTM3VtO1WkE6lLOjdJHWzJkzRxuIPDwFaztDdDAACTaFaLLoDv3yn2Z4qdpZcBdCiNxAghEhPjC+zyP5/exD1p97pBua+5qzjSm9arnRuWpBbMyNuH37NgkJCZQpUwaub4NN/SAhRpvZqRzqbhtRFtXDz88vC+5ECJFbSDAixAcgPkHD/huB/HHWl6O3niU5Xr6ADX3qFMazrDNG/1rEzt7eHrVajfG5BRCxF3j1DKdoI+i4ipBoDY8ePcLZ2fk93YkQIjeSYESIXOx+UAQbzj9ik/djAsNi9I4Zqw34uJwzPWq5UcE1T7Ln57UyZ9+AotSK2JOYWL4rtJ4LaiN+njaBuLg4unbtmol3IYTI7SQYESKXiYiJ5++rAfx5/hFn7r9IctzV1owu1QrSsYor9pYmKRcU9hQ2dKehbYAuySd/ewo0ncGTu/f55ZdfmDt3LmPGjMHFxSUzbkUI8YGQYESIXCBBo3Dq7nM2//OYv68GEBWXoHfc0EBFoxIOdKtRiLru9hgYqFIv0PcMbOgB4dpAJEZjQP/d8Nu55fDFcgDy5s3LlClTGDVqVKbckxDiwyHBiBA5lKIoXHkSwraLfuy45JfkMQxAkXwWdKriSrtKBchnlUorSGKhcPZX2DMaNPGAdujuiH/yc9/AGpXqLGZmZowZM4ahQ4dibm6e0bclhPgAqRTl37MKZE+hoaHY2NgQEhKCtbV1VldHiCyjKAo3/MPYdcWfv674cz8oIkkea1NDWpV3oV2lAlQqmAeV6g2tIK9Fh8COIXBtsy7p7DMTrPtsokTlugA8evSIfv36cezYMS5cuICHh0dG3JYQIpdK6+e3BCNCZHOKonDNL5S/r/rz99UA7j1LGoAYqw1oUDwfn1TIT+OSDpgaqdN3kcfesLE3BD/UJc08GUvHJdco6FZEL2tkZCRFihShffv2zJ8//63uSQjxYUjr57c8phEiG4pP0HD+4Uv2XX/K7qsBPAmOSpLHQAXVC9vxSQUXWpZxxsbcKP0X0iTAyblwcJLusQwmNky45MAZI0e++U8gAmBubk63bt3YsGGDBCNCiAwhwYgQ2URodBzHbgVx4MZTDvoEJpmQDEClgmputnxczpnmZZxwsDJ9+wu+uAdbvoRHpxPTClSDT5eyr/VnuLnZpniqnZ0dUVFJAyQhhHgbEowIkUUUReF2YDiHbgZy8GYg3g9fEq9J+tTU0EBFzaJ2tCjjRLNSTmnriJoajQa8l8Pe7yHu9SMfFdQZAg3HgtqIcuXKsXXrVuLi4jAyStrismfPHsqWLftu9RBCiFckGBHiPXoZEcvJu885eusZR249IyA0Otl8liaG1C+ej6YlHWlY3OHtHsEk59kt2DEYfE8lpuV1gzaLoFBNXZKXlxcLFixg0qRJjB8/Xq8T7Nq1azl69CgbNmzImDoJIT54EowIkYmiYhPwfviS43eCOHEniKt+IaTUZbyQnTkNizvQqIQD1YvYYmKYzk6oqYmPgRM/w9EZkBCbmF65NzSbBCaWetnLli3L5MmTGTNmDEeOHOGzzz7DxMSELVu2sHXrVnr06MGnn36acfUTQnzQJBgRGerkyZP8/PPPHDx4EEVRqFOnDoMHD6ZRo0ZZXbX3Iio2gQuPXnL63gtO333OhUcviUtIPvowMTSgRhE76nnko0HxfBSxt0j7MNz0uLUXdo/S9hF5LW9haDUHijRI8bTRo0dTsmRJZs2aRb9+/QAoU6YMixYt4vPPP8fAwCDFc4UQIj1kaK/IMAsXLmTAgAEUL16czp07o1ar+fPPP7ly5QqTJ09m9OjRWV3FDBcSGYe37wvOPXjJ2fsvuPw4OMXgA6CkszV13O2oUywf1Qvbpn8Ibno8vwt7xsKtvxPTVGqoNQjqfwvGaZ+wLDo6moSEBCwsLDKhokKI3ErmGRHv1dWrVylXrhxfffUVs2fP1n1rVhSF8ePHM2HCBI4fP07t2rWzuKZvT6NRuBcUwT++L7ng+xLvhy+59TQ81XPc7MypWdSOGkXsqO1un/paMBkl4jkcmQbnlyUO1wUoWAs8Z4BTmcyvgxBCIMGIeM8GDBjA1q1befjwYZLRFxqNhpIlS1KpUiXWrl2bRTVMv6DwGC49CubSo2AuPg7h0qNgQqKSDrf9tyL2FlQrbEu1wrbUKGKHSx6z91RbICYMziyCE3MhJjQx3dJJ2y+kbHvt2GAhhHhPZNIz8V6dOXOGjz/+ONlhoAYGBrRp0yZbj74ICo/h6pMQrvmFcvlxMFceh+AXkvxIl9fUBipKu1hTuVBeqrrZUqVQXhys32Hej7cVGwnnfoXjcyDqX6v0GplDrcFQ66skHVSFECI7kWBEZAi1Wk10dMof3lFRURgaZv2vm0aj4Psikuv+odzwD+W6XyjX/EJTHGL7b/aWxlRwzUulQnmoVDAv5QrYYG6chfcUFawNQk4vhMjniekqNVTsBg3GgLVzllVPCCHSKus/HUSu0KxZM+bOnUtYWBhWVlZ6x2JiYtiwYcN7HQqqKApB4bHcfhqGz9MwfALCuBkQxq2nYUTGJrzxfAtjNWXy21CugA0VXPNS3tWG/HnMMme0S3oFP4Kzi+H8SogN+9cBFZTtAA2+BbuiWVU7IYRIN+kzIjLEo0ePKFmyJPXq1eP3338nb968AISFhfH555+zdetWLl26RIkSJTL0uoqi8DQ0hjuB4dwJDON2YLj29TSMl8lMp54cK1NDSrtYU9rFhjL5rSmbPw9F7C0wMMgGgcdrigKPzmhbQW7sAOVfAZXKAEq3g3rfgEPJrKujEEL8h/QZEe+Vq6srW7ZsoV27duTPn59q1aqhVqs5e/YscXFxrFu37p0CkfCYeB4ERXAvKIL7zyK4FxTOvWcR3HsWTkQaWjpeK2hrTgknK0o6W1PS2ZpSzta42maTFo/kRAXD5Q3gvQICr+sfUxtDha5Q+2uwTbqgnRBC5BQSjIgMU7NmTdq3b8/vv//OkSNHAO0Kr1999RWtW7dO9VxFUXgREYvvi0jt63kkD55H8vB5BA+eRxIUHpOuuuSzMqG4oxXFHC3xcLSihJMVHo5WWJjkgF95TQLcOwyX1mlbQeL/syCdhQNU7QtV+oClQ5ZUUQghMlIO+MsscoLo6GiaN2/OlStXGDt2LG3atCE6Opo1a9Ywa9YsngYGMueXJfiFRPP4ZSSPX0a9ekXy6IX23/S0cAAYqKBAXnOK5LOgaD5LijlY4v7qlcfcOJPuNJMoCjz5B65thqubIcwvaR7X6trp28u0A8P3MF+JEEK8JxKMiAyx4NfleN9+wvzfd2HtVIhDz6LwD1YRXOEzKo5twuHweCpN2v9WZTtYmVDIzpzC9hYUtrd89a8FbvbmGbt+y/umSdD2A7n5F9zYDsG+SfOY5tF2Sq3SGxxLv/cqCiHE+yAdWEWqImPjeRYWQ2BYDE9DowkM1W4HhkYTEBrN09BonobGEB4T/+bCUmCkVpE/jxmutuYUtDWnkJ05BW0tdNs54tFKWkW+gHuH4PZ+uL1Hf0juawaG4N4UKnQBjxbSCiKEyLGkA6tIlqIohEbF8zwihhcRsQSFx/I8IoagsFiCwmN0r2dh2ld6H50kx9TIAFVkMEp4EG2a1Sd/HjMK2JqRP485rrZmOFqZZq+RKxkpLkrb+nH/qPb1xBsUTdJ8KjUUqa8dFVPiIzC3ff91FUKILCLBSA6m0SiExcQTHBlLcGQcL//178vIOF5GxPIiMlb7779e8ZqMawwzM1LjZGPK/esXcMljxifNGuBiY4qTjRnONqY425hia2FMnTp1yJcvH1M//SrDrp0thT+Dx+fA95Q2CPG7AAmxyec1sgD3xtrgo1gzCUCEEB8sCUayWHRcAmHR8YRFxxH6+t+oeEKj4wiNiiPkP6//pmVgXKHHysSQfNYm5LM0IZ+VCQ5WpjhYm+DwatvJxgQHa1OsTAxRqVR8990B5syZw+aRt3FyctIr68yZM5w8eZL169dnTmWzSsRzCLgE/pfB/yI89oaQZPp9/Fu+ktoAxL0JFKolj2CEEALpM5JuiqIQE68hIiaeyNgEImLjiYiJJyImgYiYeMJfvbTb+mlh0XHa7ej4VwFIPLEJyTTZZwJjtQG2FsbYWhhjZ/nqXwsT7CyNsbfUbttbmWBvaYy9pUm6l7Z/+vQplSpVwsbGhjlz5tCkSRPi4+PZuHEjX3/9NW5ubpw4cQJj4xw2ygUg6iUE3YagWxB449XrOoT5v/lc26LaoKNIA3CrC1aOmV5dIYTILqTPyBvEJWhYfeohUXEJRMUmEBWXQGRsAtFxCUTGagONqFht2uv919uZ1RqRFioVWJsaYWNmRB7z1/8ak+fVfl5zY/JaaNPymhtjZ2FMXgtjLIzVmTqxl6OjI4cOHaJTp040b94cCwsLEhISiI6OpmXLlqxevTr7ByLxMdp5PV7cS3w9vwuRQWk738gcnCtA/kraYbgFa8g8IEIIkQYfbDACMGHn9TdnymBqAxWWJoZYmRpiZWqE1attazMjrF+nmRpiY2b0Kk0bcGj3tcfV2bSzp4eHB//88w8nTpzgzJkzGBoa0qRJE0qXziFDUlUGsPkL/anWU2KaRzvU1qkcOJfT/puvBKg/6P9SQgjxVt7qL+eCBQuYMWMG/v7+lC5dmjlz5lC3bt0U8x85coRhw4Zx7do1XFxcGDlyJF5eXm9d6YxgpDbASK0iLiH1Zg5jtQFmxmosjNXaf00MMTdWY2FsiLmJIZYmasyNDbF4dczCxBDLV/9amKixMjHCwkT9KgAx0o4sya5Tj2cAlUpFnTp1qFOnTlZXJf3URpDHFV4+SEyzdAS7YmBfDOw9IJ8HOJQGKydtM5UQQoh3lu5gZP369QwZMoQFCxZQu3ZtFi9eTMuWLbl+/ToFCxZMkv/+/ft4enrSr18/1qxZw4kTJxgwYAD58uV7r6u4JmdOp4oYqVWYGxtiZmyAmZEhZsZqzF8FHuZGagzVBllaR/GeNR6nDUpsi0BeNzC2yOoaCSFErpfuDqzVq1enUqVKLFy4UJdWsmRJ2rRpw5QpU5LkHzVqFNu3b+fGjRu6NC8vLy5dusSpU6eSvUZMTAwxMYlrkYSGhuLq6potOrAKIYQQIm3S2oE1XV/7Y2Nj8fb2plmzZnrpzZo14+TJk8mec+rUqST5mzdvzvnz54mLS36J9ylTpmBjY6N7ubq6pqeaQgghhMhB0hWMBAUFkZCQgKOj/vBER0dHAgICkj0nICAg2fzx8fEEBSU/SmH06NGEhIToXo8ePUpPNYUQQgiRg7xVB9b/dsBUFCXVTpnJ5U8u/TUTExNMTGQyKCGEEOJDkK6WEXt7e9RqdZJWkMDAwCStH685OTklm9/Q0BA7O7t0VlcIIYQQuU26ghFjY2MqV67Mvn379NL37dtHrVq1kj2nZs2aSfLv3buXKlWqYGRklM7qCiGEECK3Sfe41WHDhrF06VKWL1/OjRs3GDp0KL6+vrp5Q0aPHk2PHj10+b28vHj48CHDhg3jxo0bLF++nGXLlvHNN99k3F0IIYQQIsdKd5+RTp068fz5cyZMmIC/vz9lypRh165dFCpUCAB/f398fRMXCytcuDC7du1i6NCh/PLLL7i4uDB37twsn2NECCGEENmDLJQnhBBCiEyRKfOMCCGEEEJkNAlGhBBCCJGlJBgRQgghRJaSYEQIIYQQWUqCESGEEEJkqbeaDv59ez3gJzQ0NItrIoQQQoi0ev25/aaBuzkiGAkLCwOQ1XuFEEKIHCgsLAwbG5sUj+eIeUY0Gg1+fn5YWVmluiBfThYaGoqrqyuPHj2SuVT+Q96b1Mn7kzJ5b1In70/K5L1JXVrfH0VRCAsLw8XFBQODlHuG5IiWEQMDAwoUKJDV1XgvrK2t5Rc/BfLepE7en5TJe5M6eX9SJu9N6tLy/qTWIvKadGAVQgghRJaSYEQIIYQQWUqCkWzCxMSEcePGYWJiktVVyXbkvUmdvD8pk/cmdfL+pEzem9Rl9PuTIzqwCiGEECL3kpYRIYQQQmQpCUaEEEIIkaUkGBFCCCFElpJgRAghhBBZSoKRbObBgwf07duXwoULY2ZmRtGiRRk3bhyxsbFZXbUss2DBAgoXLoypqSmVK1fm2LFjWV2lLDdlyv/bu/eoJs70D+DfgAYwpliIELqUAEKVRcQI1hsu0FWgqAdrV6wtN9GeUi0XPd6OWo22hkoFvOCKaEXU1dVzrEV3l1osCrZaQW6KsNIFQYWwKHqsguIS3t8f/pgl3KPIhOX5nJNznHfey5PJOHnyzgwTjfHjx0MsFsPMzAyzZ8/GjRs3+A5LJ0VHR0MgECAqKorvUHRGVVUVAgICYGpqiiFDhmDs2LHIzc3lOyyd0NTUhHXr1nHHYFtbW2zatAnNzc18h9bnsrKyMGvWLLzxxhsQCAT47rvvNNYzxqBQKPDGG2/AyMgIHh4euH79+guNRcmIjvnnP/+J5uZm7NmzB9evX0d8fDwSExOxZs0avkPjxbFjxxAVFYW1a9ciPz8fU6dOxbvvvotbt27xHRqvMjMzsWTJEvzyyy9IT09HU1MTvLy8UF9fz3doOiUnJwdJSUkYM2YM36HojAcPHmDKlCkYPHgw0tLSUFxcjNjYWAwbNozv0HTCli1bkJiYiISEBJSUlCAmJgZff/01du7cyXdofa6+vh7Ozs5ISEjocH1MTAzi4uKQkJCAnJwcSKVSTJ8+nXuenFYY0XkxMTHMxsaG7zB48fbbb7OwsDCNslGjRrHVq1fzFJFuqq2tZQBYZmYm36HojEePHjF7e3uWnp7O3N3dWWRkJN8h6YRVq1YxNzc3vsPQWTNmzGChoaEaZXPmzGEBAQE8RaQbALCTJ09yy83NzUwqlbKvvvqKK3v69CkzNjZmiYmJWvdPMyP9wMOHD2FiYsJ3GH3u2bNnyM3NhZeXl0a5l5cXLl68yFNUuunhw4cAMCD3k84sWbIEM2bMwLRp0/gORaecOnUKrq6umDt3LszMzCCXy7F3716+w9IZbm5u+PHHH1FaWgoAKCwsxE8//QRfX1+eI9MtN2/eRE1Njcbx2cDAAO7u7i90fO4XD8obyMrKyrBz507ExsbyHUqfu3fvHtRqNczNzTXKzc3NUVNTw1NUuocxhmXLlsHNzQ2jR4/mOxyd8Ne//hV5eXnIycnhOxSdU15ejt27d2PZsmVYs2YNsrOzERERAQMDAwQFBfEdHu9WrVqFhw8fYtSoUdDX14darcbmzZsxf/58vkPTKS3H4I6Oz5WVlVr3RzMjfUShUEAgEHT5unLlikab6upq+Pj4YO7cuVi0aBFPkfNPIBBoLDPG2pUNZJ999hmuXr2Ko0eP8h2KTrh9+zYiIyNx+PBhGBoa8h2Ozmlubsa4ceOgVCohl8vxySef4OOPP8bu3bv5Dk0nHDt2DIcPH8aRI0eQl5eHlJQUbN26FSkpKXyHppN66/hMMyN95LPPPsMHH3zQZR1ra2vu39XV1fD09MSkSZOQlJT0iqPTTRKJBPr6+u1mQWpra9tl4wNVeHg4Tp06haysLFhaWvIdjk7Izc1FbW0tXFxcuDK1Wo2srCwkJCSgsbER+vr6PEbILwsLC/z+97/XKHNwcMCJEyd4iki3rFixAqtXr+aO105OTqisrER0dDSCg4N5jk53SKVSAM9nSCwsLLjyFz0+UzLSRyQSCSQSSY/qVlVVwdPTEy4uLkhOToae3sCcwBIKhXBxcUF6ejree+89rjw9PR1+fn48RsY/xhjCw8Nx8uRJnD9/HjY2NnyHpDP++Mc/4tq1axplCxYswKhRo7Bq1aoBnYgAwJQpU9rdBl5aWgqZTMZTRLqloaGh3TFXX19/QN7a2xUbGxtIpVKkp6dDLpcDeH6dX2ZmJrZs2aJ1f5SM6Jjq6mp4eHjAysoKW7duxd27d7l1LZnoQLJs2TIEBgbC1dWVmyW6desWwsLC+A6NV0uWLMGRI0eQmpoKsVjMzR4ZGxvDyMiI5+j4JRaL2107IxKJYGpqStfUAFi6dCkmT54MpVIJf39/ZGdnIykpacDOwLY1a9YsbN68GVZWVnB0dER+fj7i4uIQGhrKd2h97vHjx/jXv/7FLd+8eRMFBQUwMTGBlZUVoqKioFQqYW9vD3t7eyiVSgwZMgQffvih9oO97O0+pHclJyczAB2+Bqpdu3YxmUzGhEIhGzduHN2+ylin+0hycjLfoekkurVX0+nTp9no0aOZgYEBGzVqFEtKSuI7JJ3x22+/scjISGZlZcUMDQ2Zra0tW7t2LWtsbOQ7tD537ty5Do8zwcHBjLHnt/du2LCBSaVSZmBgwP7whz+wa9euvdBYAsYYe/G8iRBCCCHk5QzMixEIIYQQojMoGSGEEEIIrygZIYQQQgivKBkhhBBCCK8oGSGEEEIIrygZIYQQQgivKBkhhBBCCK8oGSGEEEIIrygZIYRnHh4eiIqK4juMPhcSEoLZs2dr1aampgbTp0+HSCTCsGHDetTmwIEDGnUVCgXGjh2r1bgvo6OYBQIBvvvuux61f5l46+rqYGZmhoqKihdq/zKuXbsGS0tL1NfX9/nYpP+hZITovJCQEAgEAggEAgwePBjm5uaYPn069u/f3+7hVdbW1ti2bZvGckvblldXT7dVKBQQCATw8fFpty4mJgYCgQAeHh699dYAAN9++y2++OKLXu3zVWr9eQwaNAhWVlb49NNP8eDBA6362b59Ow4cOKBVm/j4eKhUKhQUFKC0tFSrtj1VUVGhsb+IxWI4OjpiyZIl+PXXX7Xur6OYVSoV3n333R61X758OX788UduWZskLjo6GrNmzdJ4InhfcXJywttvv434+Pg+H5v0P5SMkH7Bx8cHKpUKFRUVSEtLg6enJyIjIzFz5kw0NTV12XbTpk1QqVTcKz8/v8v6FhYWOHfuHO7cuaNRnpycDCsrq5d+L22ZmJhALBb3er+vUuvPY9++fTh9+jQWL16sVR/GxsY9nt1oUVZWBhcXF9jb28PMzEyrtto6e/YsVCoVCgsLoVQqUVJSAmdnZ43EoCc6ilkqlcLAwKBH7YcOHQpTU1Ot43/y5Am++eYbLFq0SOu2vWXBggXYvXs31Go1bzGQ/oGSEdIvGBgYQCqV4ne/+x3GjRuHNWvWIDU1FWlpad3+uhaLxZBKpdxr+PDhXdY3MzODl5cXUlJSuLKLFy/i3r17mDFjhkbdnJwcTJ8+HRKJBMbGxnB3d0deXh63/vz58xAKhbhw4QJXFhsbC4lEApVKBaD9aRpra2t8+eWXCAoKwtChQyGTyZCamoq7d+/Cz88PQ4cOhZOTE65cucK16Wgqf9u2bRq/iFt+USuVSpibm2PYsGHYuHEjmpqasGLFCpiYmMDS0hL79+/vcvsA//08LC0t4eXlhXnz5uGHH37g1qvVaixcuBA2NjYwMjLCyJEjsX37do0+2v7C9/DwQEREBFauXAkTExNIpVIoFAqN7XLixAkcPHgQAoEAISEhAIC4uDg4OTlBJBLhzTffxOLFi/H48eNu30N3TE1NIZVKYWtrCz8/P5w9exYTJkzAwoULNb5cT58+DRcXFxgaGsLW1pbbpl3F3PY0zZ07d/DBBx/AxMQEIpEIrq6uuHz5MgDNz1ahUCAlJQWpqanczM358+c7jD8tLQ2DBg3CpEmTuLLz589DIBDgzJkzkMvlMDIywjvvvIPa2lqkpaXBwcEBr732GubPn4+GhgaunYeHB8LDwxEVFYXXX38d5ubmSEpKQn19PRYsWACxWIwRI0YgLS1NIwZvb2/U1dUhMzPzRT8GMkBQMkL6rXfeeQfOzs749ttve73v0NBQjSRn//79+OijjyAUCjXqPXr0CMHBwbhw4QJ++eUX2Nvbw9fXF48ePQLw30QjMDAQDx8+RGFhIdauXYu9e/fCwsKi0/Hj4+MxZcoU5OfnY8aMGQgMDERQUBACAgKQl5cHOzs7BAUFQdvnXGZkZKC6uhpZWVmIi4uDQqHAzJkz8frrr+Py5csICwtDWFgYbt++3eM+y8vL8f3332Pw4MFcWXNzMywtLXH8+HEUFxdj/fr1WLNmDY4fP95lXykpKRCJRLh8+TJiYmKwadMmpKenA3ie+Pn4+MDf3x8qlYpLbvT09LBjxw4UFRUhJSUFGRkZWLlypVbbpSf09PQQGRmJyspK5ObmAgDOnDmDgIAAREREoLi4GHv27MGBAwewefPmLmNu7fHjx3B3d0d1dTVOnTqFwsJCrFy5st0pSOD5KRt/f39uZkqlUmHy5MkdxpuVlQVXV9cO1ykUCiQkJODixYu4ffs2/P39sW3bNhw5cgR///vfkZ6ejp07d2q0SUlJgUQiQXZ2NsLDw/Hpp59i7ty5mDx5MvLy8uDt7Y3AwECNJEYoFMLZ2VkjGSekQ731qGFCXpXg4GDm5+fX4bp58+YxBwcHblkmk7H4+HiNZaFQyEQiEffavn17p2Nt2LCBOTs7s2fPnjEzMzOWmZnJHj9+zMRiMSssLGSRkZHM3d290/ZNTU1MLBaz06dPc2WNjY1MLpczf39/5ujoyBYtWqTRpu3j7WUyGQsICOCWVSoVA8A+//xzruzSpUsMAFOpVBpxtxYfH89kMhm3HBwczGQyGVOr1VzZyJEj2dSpUzXiF4lE7OjRo52+x+DgYKavr89EIhEzNDTkHiseFxfXaRvGGFu8eDF7//33Nfpp/bm6u7szNzc3jTbjx49nq1at4pb9/Py4x5d35vjx48zU1JRbTk5OZsbGxtxyR9uqtZs3bzIALD8/v926kpISBoAdO3aMMcbY1KlTmVKp1Khz6NAhZmFh0WXMANjJkycZY4zt2bOHicViVldX12E8bePt6v9Da35+fiw0NFSjrOWR8GfPnuXKoqOjGQBWVlbGlX3yySfM29ubW2772bTsJ4GBgVxZy3566dIljTHfe+89FhIS0m28ZGAbxFcSREhvYIxBIBB0WWfFihXc9DgASCSSbvsdPHgwAgICkJycjPLycrz11lsYM2ZMu3q1tbVYv349MjIy8O9//xtqtRoNDQ24desWV0coFOLw4cMYM2YMZDKZxgW2nWk9lrm5OYDnFwS2LautrYVUKu22vxaOjo7Q0/vvhKi5uTlGjx7NLevr68PU1BS1tbVd9uPp6Yndu3ejoaEB+/btQ2lpKcLDwzXqJCYmYt++faisrMSTJ0/w7Nmzbu8KabuNLSwsuo3l3LlzUCqVKC4uxm+//YampiY8ffoU9fX1EIlEXbbVFvv/maiWfS43Nxc5OTncTAjw/BTV06dP0dDQgCFDhnTbZ0FBAeRyOUxMTHo11idPnsDQ0LDDdW33ryFDhsDW1lajLDs7u9M2LftJZ/tka0ZGRhqzJYR0hJIR0q+VlJTAxsamyzoSiQR2dnZa9x0aGooJEyagqKgIoaGhHdYJCQnB3bt3sW3bNshkMhgYGGDSpEl49uyZRr2LFy8CAO7fv4/79+93+yXZ+pRHyxdfR2UtU/l6enrtTtn85z//6bLfln46KuvoFEFrIpGI26Y7duyAp6cnNm7cyN0VdPz4cSxduhSxsbGYNGkSxGIxvv76a+46iM5oG0tlZSV8fX0RFhaGL774AiYmJvjpp5+wcOHCDt//yyopKQEAbp9rbm7Gxo0bMWfOnHZ1O0sE2jIyMuq9AFuRSCSd3uHUdl/qyXbvbt9pu0+2uH//PkaMGKH9GyADCl0zQvqtjIwMXLt2De+///4r6d/R0RGOjo4oKirChx9+2GGdCxcuICIiAr6+vnB0dISBgQHu3bunUaesrAxLly7F3r17MXHiRAQFBXX7Za+t4cOHo6amRiMhKSgo6NUxurJhwwZs3boV1dXVAJ5vl8mTJ2Px4sWQy+Wws7NDWVlZr4975coVNDU1ITY2FhMnTsRbb73FxdDbmpubsWPHDtjY2EAulwMAxo0bhxs3bsDOzq7dq/UMVFfGjBmDgoIC3L9/v0f1hUJhj+5OkcvlKC4u7lGfr1JRURG3vQjpDCUjpF9obGxETU0NqqqqkJeXB6VSCT8/P8ycORNBQUGvbNyMjAyoVKpOb0G1s7PDoUOHUFJSgsuXL+Ojjz7S+KWrVqsRGBgILy8vLFiwAMnJySgqKkJsbGyvxunh4YG7d+8iJiYGZWVl2LVrV7s7G14lDw8PODo6QqlUAni+Xa5cuYIzZ86gtLQUn3/+OXJycnp93BEjRqCpqQk7d+5EeXk5Dh06hMTExF7pu66uDjU1NSgvL8epU6cwbdo0ZGdn45tvvoG+vj4AYP369Th48CAUCgWuX7+OkpISHDt2DOvWrevxOPPnz4dUKsXs2bPx888/o7y8HCdOnMClS5c6rG9tbY2rV6/ixo0buHfvXqczQN7e3rh+/brWf/+lN1VUVKCqqgrTpk3jLQbSP1AyQvqF77//HhYWFrC2toaPjw/OnTuHHTt2IDU1lftiAJ7/eh00qPfOPnb3lz7379+PBw8eQC6XIzAwEBERERp//2Lz5s2oqKhAUlISgOd/X2Lfvn1Yt25dr85cODg44M9//jN27doFZ2dnZGdnY/ny5b3Wf08sW7YMe/fuxe3btxEWFoY5c+Zg3rx5mDBhAurq6rT+OyQ9MXbsWMTFxWHLli0YPXo0/vKXvyA6OrpX+p42bRosLCzg5OSE1atXw8HBAVevXoWnpydXx9vbG3/729+Qnp6O8ePHY+LEiYiLi4NMJuvxOEKhED/88APMzMzg6+sLJycnfPXVVxr7dWsff/wxRo4cCVdXVwwfPhw///xzh/WcnJzg6ura7R1Mr9LRo0fh5eWl1fYgA5OAtT3RTEg/pVar8dprryElJQV/+tOf+A6HEN794x//wPLly1FUVNTj00a9pbGxEfb29jh69CimTJnSp2OT/ocuYCX/E+7cuYODBw9CrVbDzc2N73AI0Qm+vr749ddfUVVVhTfffLNPx66srMTatWspESE9QjMj5H+CRCKBqakpvvzyS8ydO5fvcAghhGiBkhFCCCGE8IouYCWEEEIIrygZIYQQQgivKBkhhBBCCK8oGSGEEEIIrygZIYQQQgivKBkhhBBCCK8oGSGEEEIIrygZIYQQQgiv/g8m2pHKL7JYtgAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Plot empirical CDF\n", "plt.scatter(DJF_max_def.sort_values(by = [\"DJF Maximum Rainfall Deficit (mm)\"]),\n", " np.arange(1, DJF_max_def.size+1, dtype=float)/DJF_max_def.size,\n", " color = 'k', facecolors='none', label='Empirical CDF')\n", "\n", "# Plot estimated PDF based on maximum likelihood method\n", "bins = np.linspace(MLEGEV.ppf(0.01), MLEGEV.ppf(0.99), 100)\n", "plt.plot(bins, MLEGEV.cdf(bins), \n", " lw=2, label='Estimated CDF (GEV)')\n", "# Plot estimated PDF based on the method of L-moments\n", "bins = np.linspace(MLEnorm.ppf(0.01), MLEnorm.ppf(0.99), 100)\n", "plt.plot(bins, MLEnorm.cdf(bins), \n", " lw=2, label='Estimated CDF (Normal)')\n", "\n", "plt.xlabel(\"DJF Maximum Rainfall Deficit (mm)\")\n", "plt.legend(loc='best', frameon=False)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 71, "id": "3a789edb", "metadata": {}, "outputs": [], "source": [ "# GEV distribution describes the data better than the normal distribution. \n", "# The Kolmogorov-Smirnov test confirms that the GEV distribution fits the data better than the normal distribution, \n", "# as indicated by a lower KS statistic and a higher p-value.\n", "# The CDF plot provides a visual confirmation of these findings.\n", "# Especially, the normal distribution overestimates the probability of occurrence for moderate rainfall deficits (around 4-6 mm).\n", "# When we move towards longer return periods (more extreme events), the two curves diverge, with GEV distribution providing a better fit." ] } ], "metadata": { "kernelspec": { "display_name": "base", "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.11.0" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 5 }