Linear regression for GWAS created by GPT4

gpt4_linreg.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "59aef73f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style type=\"text/css\">\n",
       "#T_1ea9a_row0_col1, #T_1ea9a_row0_col2 {\n",
       "  background-color: #6fa7ce;\n",
       "  color: #f1f1f1;\n",
       "}\n",
       "#T_1ea9a_row0_col3, #T_1ea9a_row0_col4, #T_1ea9a_row2_col1, #T_1ea9a_row2_col2 {\n",
       "  background-color: #023858;\n",
       "  color: #f1f1f1;\n",
       "}\n",
       "#T_1ea9a_row1_col1, #T_1ea9a_row1_col2, #T_1ea9a_row1_col3, #T_1ea9a_row1_col4 {\n",
       "  background-color: #167bb6;\n",
       "  color: #f1f1f1;\n",
       "}\n",
       "#T_1ea9a_row2_col3, #T_1ea9a_row2_col4 {\n",
       "  background-color: #ede8f3;\n",
       "  color: #000000;\n",
       "}\n",
       "#T_1ea9a_row3_col1, #T_1ea9a_row3_col2, #T_1ea9a_row3_col3, #T_1ea9a_row3_col4 {\n",
       "  background-color: #fff7fb;\n",
       "  color: #000000;\n",
       "}\n",
       "#T_1ea9a_row4_col1, #T_1ea9a_row4_col2 {\n",
       "  background-color: #a1bbda;\n",
       "  color: #000000;\n",
       "}\n",
       "#T_1ea9a_row4_col3, #T_1ea9a_row4_col4 {\n",
       "  background-color: #2081b9;\n",
       "  color: #f1f1f1;\n",
       "}\n",
       "#T_1ea9a_row5_col1, #T_1ea9a_row5_col2 {\n",
       "  background-color: #056ba9;\n",
       "  color: #f1f1f1;\n",
       "}\n",
       "#T_1ea9a_row5_col3, #T_1ea9a_row5_col4 {\n",
       "  background-color: #62a2cb;\n",
       "  color: #f1f1f1;\n",
       "}\n",
       "#T_1ea9a_row6_col1, #T_1ea9a_row6_col2, #T_1ea9a_row6_col3, #T_1ea9a_row6_col4 {\n",
       "  background-color: #fdf5fa;\n",
       "  color: #000000;\n",
       "}\n",
       "#T_1ea9a_row7_col1, #T_1ea9a_row7_col2 {\n",
       "  background-color: #03466e;\n",
       "  color: #f1f1f1;\n",
       "}\n",
       "#T_1ea9a_row7_col3, #T_1ea9a_row7_col4 {\n",
       "  background-color: #dddbec;\n",
       "  color: #000000;\n",
       "}\n",
       "</style>\n",
       "<table id=\"T_1ea9a\">\n",
       "  <thead>\n",
       "    <tr>\n",
       "      <th class=\"blank level0\" >&nbsp;</th>\n",
       "      <th id=\"T_1ea9a_level0_col0\" class=\"col_heading level0 col0\" >feature</th>\n",
       "      <th id=\"T_1ea9a_level0_col1\" class=\"col_heading level0 col1\" >t_stat_custom</th>\n",
       "      <th id=\"T_1ea9a_level0_col2\" class=\"col_heading level0 col2\" >t_stat_statsmodels</th>\n",
       "      <th id=\"T_1ea9a_level0_col3\" class=\"col_heading level0 col3\" >p_value_custom</th>\n",
       "      <th id=\"T_1ea9a_level0_col4\" class=\"col_heading level0 col4\" >p_value_statsmodels</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th id=\"T_1ea9a_level0_row0\" class=\"row_heading level0 row0\" >0</th>\n",
       "      <td id=\"T_1ea9a_row0_col0\" class=\"data row0 col0\" >feature_0</td>\n",
       "      <td id=\"T_1ea9a_row0_col1\" class=\"data row0 col1\" >-0.112678</td>\n",
       "      <td id=\"T_1ea9a_row0_col2\" class=\"data row0 col2\" >-0.112678</td>\n",
       "      <td id=\"T_1ea9a_row0_col3\" class=\"data row0 col3\" >0.910309</td>\n",
       "      <td id=\"T_1ea9a_row0_col4\" class=\"data row0 col4\" >0.910309</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th id=\"T_1ea9a_level0_row1\" class=\"row_heading level0 row1\" >1</th>\n",
       "      <td id=\"T_1ea9a_row1_col0\" class=\"data row1 col0\" >feature_1</td>\n",
       "      <td id=\"T_1ea9a_row1_col1\" class=\"data row1 col1\" >0.409890</td>\n",
       "      <td id=\"T_1ea9a_row1_col2\" class=\"data row1 col2\" >0.409890</td>\n",
       "      <td id=\"T_1ea9a_row1_col3\" class=\"data row1 col3\" >0.681975</td>\n",
       "      <td id=\"T_1ea9a_row1_col4\" class=\"data row1 col4\" >0.681975</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th id=\"T_1ea9a_level0_row2\" class=\"row_heading level0 row2\" >2</th>\n",
       "      <td id=\"T_1ea9a_row2_col0\" class=\"data row2 col0\" >feature_2</td>\n",
       "      <td id=\"T_1ea9a_row2_col1\" class=\"data row2 col1\" >1.194775</td>\n",
       "      <td id=\"T_1ea9a_row2_col2\" class=\"data row2 col2\" >1.194775</td>\n",
       "      <td id=\"T_1ea9a_row2_col3\" class=\"data row2 col3\" >0.232462</td>\n",
       "      <td id=\"T_1ea9a_row2_col4\" class=\"data row2 col4\" >0.232462</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th id=\"T_1ea9a_level0_row3\" class=\"row_heading level0 row3\" >3</th>\n",
       "      <td id=\"T_1ea9a_row3_col0\" class=\"data row3 col0\" >feature_3</td>\n",
       "      <td id=\"T_1ea9a_row3_col1\" class=\"data row3 col1\" >-1.474552</td>\n",
       "      <td id=\"T_1ea9a_row3_col2\" class=\"data row3 col2\" >-1.474552</td>\n",
       "      <td id=\"T_1ea9a_row3_col3\" class=\"data row3 col3\" >0.140652</td>\n",
       "      <td id=\"T_1ea9a_row3_col4\" class=\"data row3 col4\" >0.140652</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th id=\"T_1ea9a_level0_row4\" class=\"row_heading level0 row4\" >4</th>\n",
       "      <td id=\"T_1ea9a_row4_col0\" class=\"data row4 col0\" >feature_4</td>\n",
       "      <td id=\"T_1ea9a_row4_col1\" class=\"data row4 col1\" >-0.432544</td>\n",
       "      <td id=\"T_1ea9a_row4_col2\" class=\"data row4 col2\" >-0.432544</td>\n",
       "      <td id=\"T_1ea9a_row4_col3\" class=\"data row4 col3\" >0.665441</td>\n",
       "      <td id=\"T_1ea9a_row4_col4\" class=\"data row4 col4\" >0.665441</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th id=\"T_1ea9a_level0_row5\" class=\"row_heading level0 row5\" >5</th>\n",
       "      <td id=\"T_1ea9a_row5_col0\" class=\"data row5 col0\" >feature_5</td>\n",
       "      <td id=\"T_1ea9a_row5_col1\" class=\"data row5 col1\" >0.593910</td>\n",
       "      <td id=\"T_1ea9a_row5_col2\" class=\"data row5 col2\" >0.593910</td>\n",
       "      <td id=\"T_1ea9a_row5_col3\" class=\"data row5 col3\" >0.552708</td>\n",
       "      <td id=\"T_1ea9a_row5_col4\" class=\"data row5 col4\" >0.552708</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th id=\"T_1ea9a_level0_row6\" class=\"row_heading level0 row6\" >6</th>\n",
       "      <td id=\"T_1ea9a_row6_col0\" class=\"data row6 col0\" >feature_6</td>\n",
       "      <td id=\"T_1ea9a_row6_col1\" class=\"data row6 col1\" >-1.434138</td>\n",
       "      <td id=\"T_1ea9a_row6_col2\" class=\"data row6 col2\" >-1.434138</td>\n",
       "      <td id=\"T_1ea9a_row6_col3\" class=\"data row6 col3\" >0.151849</td>\n",
       "      <td id=\"T_1ea9a_row6_col4\" class=\"data row6 col4\" >0.151849</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th id=\"T_1ea9a_level0_row7\" class=\"row_heading level0 row7\" >7</th>\n",
       "      <td id=\"T_1ea9a_row7_col0\" class=\"data row7 col0\" >feature_7</td>\n",
       "      <td id=\"T_1ea9a_row7_col1\" class=\"data row7 col1\" >1.058532</td>\n",
       "      <td id=\"T_1ea9a_row7_col2\" class=\"data row7 col2\" >1.058532</td>\n",
       "      <td id=\"T_1ea9a_row7_col3\" class=\"data row7 col3\" >0.290072</td>\n",
       "      <td id=\"T_1ea9a_row7_col4\" class=\"data row7 col4\" >0.290072</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n"
      ],
      "text/plain": [
       "<pandas.io.formats.style.Styler at 0x10ab0f340>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import statsmodels.api as sm\n",
    "import pandas as pd\n",
    "import scipy.stats as stats\n",
    "from scipy import linalg\n",
    "import matplotlib.pyplot as plt \n",
    "\n",
    "def compute_t_stat_p_value(X, y, covariates):\n",
    "    \"\"\"\n",
    "    Compute t-statistics and p-values for each feature in an OLS regression.\n",
    "\n",
    "    Parameters:\n",
    "    X (numpy.ndarray): The feature matrix of shape (n_samples, n_features)\n",
    "    y (numpy.ndarray): The outcome vector of shape (n_samples,)\n",
    "    covariates (numpy.ndarray): The shared covariate matrix of shape (n_samples, n_covariates)\n",
    "\n",
    "    Returns:\n",
    "    t_statistics (numpy.ndarray): t-statistics for each feature of shape (n_features,)\n",
    "    p_values (numpy.ndarray): p-values for each feature of shape (n_features,)\n",
    "    \"\"\"\n",
    "    n_samples, n_features = X.shape\n",
    "    n_covariates = covariates.shape[1]\n",
    "    t_statistics = np.zeros(n_features)\n",
    "    p_values = np.zeros(n_features)\n",
    "\n",
    "    # Adding the shared covariates to the feature matrix\n",
    "    X_with_covariates = np.hstack((X, covariates))\n",
    "\n",
    "    for i in range(n_features):\n",
    "        # Perform OLS regression on feature i with shared covariates\n",
    "        X_i = X_with_covariates[:, [i] + list(range(n_features, n_features + n_covariates))]\n",
    "        X_i = np.hstack((np.ones((n_samples, 1)), X_i))  # Add intercept\n",
    "\n",
    "        beta_i = linalg.lstsq(X_i, y)[0]\n",
    "        residuals = y - X_i @ beta_i\n",
    "\n",
    "        # Compute the standard error of the beta_i\n",
    "        residual_variance = np.sum(residuals ** 2) / (n_samples - n_covariates - 2)\n",
    "        XTX_inv = linalg.inv(X_i.T @ X_i)\n",
    "        se_beta_i = np.sqrt(residual_variance * XTX_inv[1, 1])\n",
    "\n",
    "        # Compute the t-statistic and p-value\n",
    "        t_stat = beta_i[1] / se_beta_i\n",
    "        p_value = 2 * (1 - stats.t.cdf(np.abs(t_stat), df=n_samples - n_covariates - 2))\n",
    "\n",
    "        t_statistics[i] = t_stat\n",
    "        p_values[i] = p_value\n",
    "\n",
    "    return t_statistics, p_values\n",
    "\n",
    "\n",
    "\n",
    "def compute_t_stat_p_value_short(X, y, covariates):\n",
    "    n_samples, n_features = X.shape\n",
    "    t_statistics = np.zeros(n_features)\n",
    "    p_values = np.zeros(n_features)\n",
    "\n",
    "    for i in range(n_features):\n",
    "        X_i = np.hstack((np.ones((n_samples, 1)), X[:, [i]], covariates))\n",
    "        beta_i, _, _, _ = linalg.lstsq(X_i, y)\n",
    "        residuals = y - X_i @ beta_i\n",
    "        residual_variance = np.sum(residuals ** 2) / (n_samples - covariates.shape[1] - 2)\n",
    "        se_beta_i = np.sqrt(residual_variance * linalg.inv(X_i.T @ X_i)[1, 1])\n",
    "        t_stat = beta_i[1] / se_beta_i\n",
    "        p_value = 2 * (1 - stats.t.cdf(np.abs(t_stat), df=n_samples - covariates.shape[1] - 2))\n",
    "        t_statistics[i], p_values[i] = t_stat, p_value\n",
    "\n",
    "    return t_statistics, p_values\n",
    "\n",
    "# def compute_t_stat_p_value_short_qr(X, y, covariates):\n",
    "#     n_samples, n_features = X.shape\n",
    "#     t_statistics = np.zeros(n_features)\n",
    "#     p_values = np.zeros(n_features)\n",
    "\n",
    "#     # Project the covariates out of the outcome and the features using QR decomposition\n",
    "#     Q, _ = linalg.qr(covariates)\n",
    "#     Q = Q[:, :covariates.shape[1]]  # Take only the first n_covariates columns of Q\n",
    "#     X_proj = X - Q @ Q.T @ X\n",
    "# #     y -= y.mean()\n",
    "#     y_proj = y - Q @ Q.T @ y\n",
    "# #     y_proj = y_proj / np.sqrt(np.sum(y_proj**2) / (n_samples - Q.shape[1]))\n",
    "\n",
    "#     dof = n_samples - covariates.shape[1] - 2\n",
    "#     for i in range(n_features):\n",
    "#         X_i = np.vstack((np.ones(n_samples), X_proj[:, i])).T\n",
    "#         beta_i, _, _, _ = linalg.lstsq(X_i, y_proj)\n",
    "#         residuals = y_proj - X_i @ beta_i\n",
    "#         residual_variance = np.sum(residuals ** 2) / dof\n",
    "#         se_beta_i = np.sqrt(residual_variance * linalg.inv(X_i.T @ X_i)[1, 1])\n",
    "#         t_stat = beta_i[1] / se_beta_i\n",
    "#         p_value = 2 * (1 - stats.t.cdf(np.abs(t_stat), df=n_samples - covariates.shape[1] - 2))\n",
    "#         t_statistics[i], p_values[i] = t_stat, p_value\n",
    "\n",
    "#     return t_statistics, p_values\n",
    "\n",
    "def compute_t_stat_p_value_short_qr(X, y, covariates):\n",
    "    n_samples, n_features = X.shape\n",
    "    t_statistics = np.zeros(n_features)\n",
    "    p_values = np.zeros(n_features)\n",
    "\n",
    "    # Demean the outcome\n",
    "    y_mean = np.mean(y)\n",
    "    y_demeaned = y - y_mean\n",
    "\n",
    "    # Project the covariates out of the outcome and the features using QR decomposition\n",
    "    Q, _ = linalg.qr(covariates)\n",
    "    Q = Q[:, :covariates.shape[1]]  # Take only the first n_covariates columns of Q\n",
    "    X_proj = X - Q @ Q.T @ X\n",
    "    y_proj = y_demeaned - Q @ Q.T @ y_demeaned\n",
    "\n",
    "    # Scale the projected outcome by its standard deviation\n",
    "    y_proj_std = np.std(y_proj)\n",
    "    y_proj_scaled = y_proj / y_proj_std\n",
    "\n",
    "    for i in range(n_features):\n",
    "        X_i = np.vstack((np.ones(n_samples), X_proj[:, i])).T\n",
    "        beta_i, _, _, _ = linalg.lstsq(X_i, y_proj_scaled)\n",
    "        residuals = y_proj_scaled - X_i @ beta_i\n",
    "        residual_variance = np.sum(residuals ** 2) / (n_samples - covariates.shape[1] - 2)\n",
    "        se_beta_i = np.sqrt(residual_variance * linalg.inv(X_i.T @ X_i)[1, 1])\n",
    "        t_stat = beta_i[1] / se_beta_i\n",
    "        p_value = 2 * (1 - stats.t.cdf(np.abs(t_stat), df=n_samples - covariates.shape[1] - 2))\n",
    "        t_statistics[i], p_values[i] = t_stat, p_value\n",
    "\n",
    "    return t_statistics, p_values\n",
    "\n",
    "\n",
    "def compare_results(X, y, covariates, t_statistics_custom, p_values_custom, feature_names=None):\n",
    "    \"\"\"\n",
    "    Compare t-statistics and p-values from the custom function to those from the statsmodels library.\n",
    "\n",
    "    Parameters:\n",
    "    X (numpy.ndarray): The feature matrix of shape (n_samples, n_features)\n",
    "    y (numpy.ndarray): The outcome vector of shape (n_samples,)\n",
    "    covariates (numpy.ndarray): The shared covariate matrix of shape (n_samples, n_covariates)\n",
    "    t_statistics_custom (numpy.ndarray): t-statistics from custom function of shape (n_features,)\n",
    "    p_values_custom (numpy.ndarray): p-values from custom function of shape (n_features,)\n",
    "    feature_names (list, optional): List of feature names for prettier output\n",
    "\n",
    "    Returns:\n",
    "    comparison_df (pandas.DataFrame): DataFrame comparing the custom and statsmodels results\n",
    "    \"\"\"\n",
    "    n_features = X.shape[1]\n",
    "    if feature_names is None:\n",
    "        feature_names = [f'Feature_{i}' for i in range(n_features)]\n",
    "\n",
    "    t_statistics_statsmodels = np.zeros(n_features)\n",
    "    p_values_statsmodels = np.zeros(n_features)\n",
    "\n",
    "    for i in range(n_features):\n",
    "        # Perform OLS regression on feature i with shared covariates using statsmodels\n",
    "        X_i = np.hstack((X[:, [i]], covariates))\n",
    "        X_i = sm.add_constant(X_i)  # Add intercept\n",
    "\n",
    "        model = sm.OLS(y, X_i)\n",
    "        results = model.fit()\n",
    "\n",
    "        t_statistics_statsmodels[i] = results.tvalues[1]\n",
    "        p_values_statsmodels[i] = results.pvalues[1]\n",
    "\n",
    "    # Create a comparison DataFrame\n",
    "    comparison_df = pd.DataFrame({\n",
    "        'feature': feature_names,\n",
    "        't_stat_custom': t_statistics_custom,\n",
    "        't_stat_statsmodels': t_statistics_statsmodels,\n",
    "        'p_value_custom': p_values_custom,\n",
    "        'p_value_statsmodels': p_values_statsmodels\n",
    "    })\n",
    "\n",
    "    return comparison_df\n",
    "\n",
    "\n",
    "# Generate random data for the small dataset\n",
    "np.random.seed(42)\n",
    "n_samples = 1000\n",
    "n_features = 50\n",
    "n_covariates = 10\n",
    "\n",
    "X = np.random.randn(n_samples, n_features)\n",
    "y = np.random.randn(n_samples)\n",
    "covariates = np.random.randn(n_samples, n_covariates)\n",
    "\n",
    "# Calculate t-statistics and p-values using the custom function\n",
    "t_statistics_custom, p_values_custom = compute_t_stat_p_value_short(X, y, covariates)\n",
    "\n",
    "# Compare the custom function results to those from the statsmodels library\n",
    "feature_names = [f'feature_{i}' for i in range(n_features)]\n",
    "comparison_df = compare_results(X, y, covariates, t_statistics_custom, p_values_custom, feature_names)\n",
    "\n",
    "np.testing.assert_allclose(comparison_df['t_stat_custom'], comparison_df['t_stat_statsmodels'])\n",
    "np.testing.assert_allclose(comparison_df['p_value_custom'], comparison_df['p_value_statsmodels'])\n",
    "\n",
    "comparison_df.head(8).style.background_gradient()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b661853e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABR0AAAFzCAYAAACzTa7TAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAACMjElEQVR4nOzdeVjVZf7/8ddh3xHZ3FBAUDDN3FNzTU1tLJupbLGsrMYWySxL6/dtsmmyTVssLWfMsoVxWqeZXCt3rdTALAVENkVUQNllO+fz+4NkRBBBDxyW5+O6vK7Ofe778MZP4u37Xt4mwzAMAQAAAAAAAICV2Nk6AAAAAAAAAAAtC0lHAAAAAAAAAFZF0hEAAAAAAACAVZF0BAAAAAAAAGBVJB0BAAAAAAAAWBVJRwAAAAAAAABWRdIRAAAAAAAAgFWRdAQAAAAAAABgVQ62DqAxWSwWHT16VJ6enjKZTLYOBwAAoF4Mw1B+fr46dOggOzvWjpsj5qMAAKC5q+uctFUlHY8ePaqgoCBbhwEAAHBJDh8+rE6dOtk6DFwE5qMAAKCluNCctFUlHT09PSVV/KZ4eXnZOBoAAICamS2G/rP3qBZ/f1DH80okSd3beWrG4Pa6ccQVlXMaND/MRwEAQHMRl5GnRd8maEditiTJw9le04eF6PoebRUeGnzBOWmrSjqeOcLi5eXFJA8AADQ5hmFoc0KmXlwTp7hj+ZLs1SmgrR6/prsmX9FRBQX5ksSx3GaM+SgAAGjq0nNOa+H6eH0Zky7DkJzd3DX1yi6aOTpcbd2dlJeXJ+nCc9JWlXQEAABoqn5Nz9WCNQe0/feVZE8XBz08KkzThgTLxdHextEBAACgpcs9XaYlmxK1YnuKSsstkqQ/XN5ec67pri6+7vX+PJKOAAAANnT4ZJEWro/XV7FHJUlO9na6c3AXPTQqTD7uTjaODgAAAC1dSblZH+5M1VsbE5VTVCZJGhTSVk9NjFTvoDYX/bkkHQEAAGwgp6hUb29M1Ac7UlVqrlhJnnxFBz02rruC2rrZODoAAAC0dBaLof/8clSvrIvXkVOnJUnhAR6aOyFCoyMCLvlKH5KOAAAAjai4zKyVO1P01veJyisulyQN6eqreRMi1auTt42jAwAAQGuwIzFLC9bEaV96riQpwNNZj43rpj/17SQHezurfA2SjgAAAI3AYjH0773penVdgtJzKlaSI9p5au6ECI3o5k9xGAAAADS4uGN5enFNnDbFZ0qSPJwdNGNEqO65KkRuTtZNE5J0BAAAaGBbD2Zqweo47c+oqPTXzstFj43rpj/27SR7O5KNAAAAaFgZuae1aH2CPvv5iAxDcrAz6fZBnTXz6nD5eTg3yNck6QgAANBA9h/N04I1B7T1YJYkydPZQQ+M6qp7hoZQkRoAAAANLq+4TO9sOqTl25JV8ntF6mt7tdfj13RXiF/9K1LXB0lHAAAAK0jKLFDqySIF+7rLycFOC9fH68uYdBmG5Ghv0h1XBuvh0WFqS0VqAAAAWFlSZoF+TD4pk6RBob7q2MZVH/+Yqje/O6hTv1ekHhDso3kTI9W3s0+jxETSEQAA4BLkFJUqKjpWWw5mVraZTJJhVPz3pN4dNGdcd3X2pSI1AAAArCunqFQPfPSzdiZlV2l3drCr3NnY1d9dcydEakzkpVekrg+SjgAAAJcgKjpW285KOEoVCUdvV0etvGegege1sU1gAAAAaPGiomOrJRwlqaTcIid7Oz173WW6ub/1KlLXB0lHAACAi5R4Ir/KDsez5Z4uk5erYyNHBAAAgNYiKbPgvHNRSSo1WzS4q69NEo6SZJuvCgAA0MztSMzS9A9219onJbuwkaIBAABAaxN7JOeCfWw5H2WnIwAAQD3EHcvTi2vitCn+/KvKZwT7NmxFQAAAALQ++cVlWrYlScu2JF2wry3noyQdAQAA6iAj97QWrU/QZz8fkWFIDnYmTb2yi+Iy8rUr5aTMZyrHSLI3mTQ0zE8hfiQdAQAAcPGSMguUerJIwb7u6tjGVdE/penN7w4qu7BUkuTl4qC84vIaxw4P97fpfJSkIwAAwDnOntz5ejhp6aZDem9bcmUFwGt7tdeca7or2M9duUVlmhkdU+U+naFhflp8ax9bhQ8AAIBmLqeoVFHRsVXmmC6O9iouM0uSQv3c9cT4CF0Z0lYPfFy9evWQrr42n4+SdAQAAFBFonH/0Tx9sCNFu1JPVbY72JlUbqnYxTgwuK3mTYxQn84+le97uzlq5fSBSs4qVEp2oYJ93dnhCAAAgEsSFR2r7YlZVdqKy8xytDfpmUmX6ZYBQXL8vUBM9P1XKjmrUD8kZcskaVCob5OYj5J0BAAArVpNq8hnK7cYcnOy1xu39NGYyACZTKYa+4X4kWwEAADApautKnWZ2dBVYX6VCcczmuJclOrVAACgVatpFflcRaVmhQV4nDfhCAAAAFjDibxiPfff/bX2sWVF6vpgpyMAAGi1altFPldKdmGTWz0GAABA83X2PeL+ns5atiVJf9+SpNO/39t4PrasSF0fJB0BAECrFXs4p859m8vkDgAAAE1bTdf7ONqbVGauuEf8iqA2MgxDv6bnyWwYlX3sTSYNDfNrNgvhJB0BAECrk19cpnc3J+nvW5Mu2Le5Te4AAADQtEVFx2rbOadtysyGXBzttejm3prQs53yTpdrZnRMlcTk0DA/m1ekrg+SjgAAoNUoLbco+qc0vfHdQZ0sLJUkebk4qKC4XJbzjGlukzsAAAA0XbVd71NcZlZkey+ZTCZ5uzlq5fSBSs4qVEp2oYJ9m16hmAsh6QgAAFo8wzC05tdjenltnFKyiyRJoX7uenJChAYFt1XUP6sebxkQ7KNpQ4J1WQfvZje5AwAAQNOUlFmgOZ/9Umufc+8Rb4pVqeuKpCMAAGgxzr6M+8zk7Kfkk3ph9YHK+xv9PJw0a0w3TRkQJEd7O0lq9qvIAAAAaLoy80v0xncJiv7psMwWo9a+LekecZKOAACg2cspKtW9H+zW7tRTlW39u/jIw9lBmxIqdjC6OdnrvmGhum94qDycq0+BmvMqMgAAAJqewpJy/WNrspZtOaTC0oqK1FdHBCinqEyxh3OadZGYuiDpCAAAmrWcolKNenWTThWVVWk/k4C0tzNpyoAgzbo6XAFeLrYIEQAAAK1Iudmif+0+ote+TVBmfokkqXcnb82dEKnBXX2VW1TW7IvE1AVJRwAA0Kzdt3J3tYTj2d6b1l8jugc0YkQAAABojQzD0LcHTujFNQd0KLNQktS5rZueGN9d1/ZqL5PJJEktokhMXZB0BAAAzVZSZoF2pZyqtc/5qlIDAAAA1hKTdkoLVsfpp5STkiQfN0dFXR2u2wd1kZODXY1jWvr1PiQdAQBAs3F2oZhgXzd9EZN+wTEt6TJuAAAANC0pWYV6ZV28vtmXIUlydrDT9KtCNGNkV3m5ONo4Otsi6QgAAJq8nKJSRUXHVrn3xsvFQXnF5bWOG9DFp0WvHgMAAMA2sgtK9OZ3B/Xxj2kqtxgymaQb+3bS7HHd1N7b1dbhNQkkHQEAQJMXFR2r7YlZVdryistlZ5I6+bjpyMmiaseofdwc9Y9pAxovSAAAALR4p0vNWr4tSe9sTlJBScUC+Mju/po7IUIR7bxsHF3TQtIRAAA0aUmZBVV2OJ7NYkhv3dZHr65LqNJnQLCP/nHnAHm7te4jLQAAALAOs8XQZ3sOa9GGBB3Pq6hI3bOjl56aEKkhYX42jq5pIukIAACarMKScr35/cFa+2QXlraK6n8AAABofIZhaGP8Cb24Jk4JxwskSZ18XDXnmu6adHkH2dmZbBxh00XSEQAANDnlZov+tfuIXvs2QZn5JbX2PVMopqVX/wMAAEDj2ns4RwvWHNAPSRUVqb1dHTVzdJjuGNxFzg72No6u6SPpCAAAmgzDMLRh/3G9tDZOhzILJUmd27rJ3cle8cfyq9zbaG8yaWiYH4lGAAAAXJKkzAKlniyqPDGTml1Rkfq/v1RUpHZysNPdQ4P14Igwru+pB5KOAACgSfg57ZQWrD6gXSmnJFUUgnnk6nDdNqiLTpeaNTM6psq9jUPD/LT41j62ChcAAADNXE5RqaKiY6vMMTu2cdXxvOLKitQ39Omox8Z1V8c2VKSuL5KOAADAppKzCvXKujit3ndMkuTsYKd7h4XozyO6ysulYiXZycGOexsBAABgVVHRsdqemFWlLT3ntCRpeDd/zR0foR4dqEh9sUg6AgCARrE5/oRij+Sob2cfDQv3V1ZBiRZ/d1Af/5hWuZJ8U79OenRsN7X3rnklmXsbAQAAcCnOHKW2N6nKDsdzzb/uMuadl4ikIwAAaFCp2YWa/PZ2nSoqq2xzcbSTncmkolKzJGlUd389OSFCEe1YSQYAAID11XSUujYp2YUkHS8RSUcAANCgzk04SlJxWUVJmF4dvTVvQoSGhPnZIjQAAAC0EjUdpa5NsC8Jx0tF0hEAADSYzfEnqiUcz/b4Nd1IOAIAAKBBJWUW1HmHo73JpKFhfuxytAI7WwdQH1u2bNGkSZPUoUMHmUwmffXVV7YOCQAA1GLtb8dqfT/2cE7jBAIAAIBWKaeoVC+tjatz/6Fhflp8a58GjKj1aFY7HQsLC9W7d2/dfffd+tOf/mTrcAAAwHmkZhfqlXXx+u8vGbX269vZp5EiAgAAQGtSXGbWBztS9PbGROUVl9fa98PpA1VuMRTsS9FCa2pWSccJEyZowoQJtg4DAACcx8nCUi3+/qA++iFVZeaKitSOdiaVmo1qfX3cHDUs3N8GUQIAAKClslgMfRmTroXr43U0t1iSFNHOUw52dtp/NFeWs/qeOUrNnLRhNKukY32VlJSopKSk8nVeXp4NowEAoOU6XWrWe9uT9c6mQ8ovqVhJHt7NX3PHR8jT2UHXvb2tyt2OPm6O+vqhq2wVLgAAAFqgLQmZWrAmTgcyKvI/7b1d9Ni47rqhT0cVFJdrZnRMlbsdOUrdsFp00nHBggWaP3++rcMAAKDFMlsMff7zES1an6BjeRUryT3ae2nexIgqK8Yxz4zT1oOZ+jntlPp29mE1GQAAAFbz29FcvbgmTlsPVlSn9nRx0IMjw3T30GC5ONpLkrzdHLVy+kAlZxUqJbuQo9SNoEUnHefNm6fZs2dXvs7Ly1NQUJANIwIAoGUwDEObEjL14uo4xR/PlyR1bOOqOdd013W9O8jOzlRtzLBwf5KNAAAAsJojp4q0aH2CvoxNl2FIjvYm3Tk4WA+PCpOPu1ONY0L8SDY2lhaddHR2dpazs7OtwwAAoEXZdyRXC9Yc0I5D2ZIkLxcHzRwdrjsGd6lcSQYAAAAaSm5RmZZsStSKHSkqLa+4pfG63h30+Lju6uzrZuPocEaLTjoCAADrOXyySK+uj9e/Y49Kkpzs7XTX0GA9OLKr2rjVvJIMAAAAWEtxmVkf7kzVWxsTlXu64r7wK0Pb6qmJkbq8UxvbBodqmlXSsaCgQImJiZWvk5OTFRsbq7Zt26pz5842jAwAgJbrVGGp3tqYqA93pqrUXLGSfEOfjpo9tpuC2rKSDAAAgIZlsRj6eu9RvbIuXuk5pyVJ3QI9NG9CpEZ295fJVP1qH9iena0DqI/du3erT58+6tOnorLQ7Nmz1adPHz3zzDM2jgwAgJanuMysdzYf0vBXNmr5tmSVmi0aGuar/868Sq9NuYKEI1qEJUuWKCQkRC4uLurXr5+2bt1aa/+PP/5YvXv3lpubm9q3b6+7775b2dnZjRQtAACtz/bELE16a5tmrYpVes5pBXo56+UbL9eaR4ZrVEQACccmrFntdBw5cqQMw7B1GAAAtGhmi6GvYtK1cH28juZWVKSOaOepeRMjNTzcj4kdWoxVq1Zp1qxZWrJkiYYOHap3331XEyZM0P79+2s8RbNt2zbdeeedeu211zRp0iSlp6drxowZuvfee/Xll1/a4DsAAKDlOpCRpxfXxGlzQqYkycPZQQ+M7Kp7hobI1Yl7xJuDZpV0BAAADWtLQqYWrInTgYw8SVIHbxc9Nq67JvfpKPsaKlIDzdmiRYs0ffp03XvvvZKk119/XevWrdPSpUu1YMGCav1/+OEHBQcHKyoqSpIUEhKiP//5z3r55ZcbNW4AAFqKpMwCpZ4sUrDv/ypKH805rUUbEvT5z0dkGJKDnUlTr+yimaPD5OtBseDmhKQjAADQr+m5emltnLYezJIkebo46KFRYbprSDAVqdEilZaWas+ePZo7d26V9nHjxmnHjh01jhkyZIiefvpprV69WhMmTNCJEyf02Wef6dprrz3v1ykpKVFJSUnl67y8POt8AwAANGM5RaWKio7VloOZlW1DQn3VvZ2nPvkpTSW/V6S+9vL2mjOuu4J/T0iieSHpCABAK3bkVJEWrk/QV7HpMgzJ0d6kOwcH6+FRYfJxpyI1Wq6srCyZzWYFBgZWaQ8MDNSxY8dqHDNkyBB9/PHHmjJlioqLi1VeXq7rrrtOixcvPu/XWbBggebPn2/V2AEAaO6iomO1PTGrStuOpGztSKq4J3lgSFvNmxChPp19bBEerKRZFZIBAADWkVtUphdWH9DoVzfry5iKhON1vTvo+8dG6v/+0IOEI1qNc+8oNQzjvPeW7t+/X1FRUXrmmWe0Z88erV27VsnJyZoxY8Z5P3/evHnKzc2t/HX48GGrxg8AQHOTlFmgLQczZT5PzY6/3dBTq+6/koRjC8BORwAAWpHiMrM+3JmqtzYmKvd0mSRpcKivnpoYqV6dvG0cHdB4/Pz8ZG9vX21X44kTJ6rtfjxjwYIFGjp0qObMmSNJuvzyy+Xu7q5hw4bp+eefV/v27auNcXZ2lrMz908BAHBG6smiWt/v0MaVwoUtBElHAABaAYvF0Nd7j+qVdfFKzzktSeoe6Km5EyM0sps/Ezu0Ok5OTurXr582bNigG264obJ9w4YNuv7662scU1RUJAeHqtNne/uKO0+N8+zWAAAA/xN/LF/vbj5Ua59gX+5vbClIOgIA0MJtT8zSC6sP6LejFQUs2nm5aPa4bvpT305UpEarNnv2bN1xxx3q37+/Bg8erGXLliktLa3yuPS8efOUnp6ulStXSpImTZqk++67T0uXLtU111yjjIwMzZo1SwMHDlSHDh1s+a0AANCkHcst1msbEvTpnsOyGNKZGejZS3b2JpOGhvlVVrFG80fSEQCAFupARp5eXBOnzQkVVQE9nR00Y2RX3TM0RK5OVKQGpkyZouzsbD333HPKyMhQz549tXr1anXp0kWSlJGRobS0tMr+d911l/Lz8/XWW2/pscceU5s2bTR69Gi99NJLtvoWAABo0vKKy/Tu5kNavi1ZxWUVFakn9GynGSO6auH6hCrVq4eG+WnxrX1sFSoagMloRWdB8vLy5O3trdzcXHl5edk6HAAArCIps0CpJ4sU7OuuED93Hc05rYXrE/RFzJHKitRTr+yimaPD1ZYCMc0ac5nmj2cIAGgNSsst+uTHVL35faJOFpZKkvp38dG8iZHq1+V/BWKSswqVkl1YOY9F81DX+Qw7HQEAaKZyikoVFR1bZYU4yMdVJ/JLVFJesZL8h8vba8413dWFu3EAAADQwAzD0Op9x/TyujilZlcUjAn1d9fc8REa2yOw2j3iIX4kG1syko4AADRTUdGx2p6YVaXt8KmKIjEDQ9rqqYmRuiKojQ0iAwAAQEt27kkbSfoxKVsvrInT3sM5kiQ/D2c9OjZcU/oHycHezobRwlZIOgIA0AwlZRZU2eF4rhf/2Euh/h6NGBEAAABauppO2vTv4iMPZwdt+v0ecTcne90/PFT3DQuVuzNpp9aMpw8AQDO09rdjtb6ferKIpCMAAACsqqaTNrtTT0mS7O1MumVAkB4ZE64ATxdbhIcmhqQjAABN2LlHV+KP5evFNQe0Mf78uxwlKZg7HAEAAGBFFzpp8960/hrRPaARI0JTR9IRAIAmqKajK4FeLsrML5bFkBzsTArwdNax3GJZzhpnbzJpaJgfF3IDAADAqg5lFdT6vqXWd9EacZMnAABNUE1HV47nVSQcJ/Zqpw2zR2jNI8N1Vbh/lT5Dw/y0+NY+jRkqAAAAWjDDMLRmX4bmf72/1n6ctMG52OkIAEATsDn+hGKP5KhvZx91bONa69GVOddEVO5kXDl9oJKzCpWSXVileiAAAABwqXannNQLqw/o57QcSZKjvUnlZkPGWX04aYPzIekIAIANpWYXavLb23WqqKyyzdnBVOuYlOzCKpO6ED+SjQAAALCexBMFenltnNbvPy5JcnW0133DQnTLgM6a+8W+KgvknLTB+ZB0BADAhs5NOEpSSblxnt4VOLoCAACAhnAiv1hvfHtQ/9x1WGaLITuTNGVAkGaN6aZAr4qK1Jy0QV2RdAQAwEY2x5+olnCsDUdXAAAA0BAKS8r1961JWrYlSUWlZknSmMhAPTm+u8IDPav156QN6oKkIwAANrLtnEIx5+ri66rU7NOVrzm6AgAAAGsqN1u0avdhvbbhoLIKSiRJvYPa6KkJERoU6mvj6NDckXQEAKCBJWUWKPVkUeXxk/ziMi3bkqQPdqbWOu75yb3UyceNoysAAACwKsMwtH7/cb20Nk5JmYWSpC6+bnrimghN7NVOJlPtd4wDdUHSEQCABpJTVKr7Vu7WrpRTlW1h/h7KLiypPFZtb2eS2VL9DkcfN0cNC/eXJJKNAAAAuCjnLn5L0p7UU1qw+oB2p1bMUdu6OylqdJhuG9RFTg52tgwXLQxJRwAAGkBs2ilNXf6jCkrMVdoTMwskSaF+7npifHf1aOel65dULSbj4+aorx+6qlHjBQAAQMuRU1SqqOjYKlWm+3fxURs3R3174IQkycXRTvdeFao/jwiVp4ujrUJFC0bSEQAAK6ppgleTd+/oV3kpd8wz47T1YKZ+Tjulvp19Knc4AgAAABcjKjpW28+5P/zMzkY7k3RTvyA9Orab2nm72CI8tBIkHQEAsKKo6Fhtu0DCUZKO5JyuUglwWLg/yUYAAABckqTMAv2YfLLWBfB/3NlfoyMDGzEqtFYkHQEAsJKkzIIL7nA8I9iXexoBAABgHXU9bSNJJjuKxKBxkHQEAOAinHspd2FJud747mCdxg4I9qE4DAAAAC7ZmTnpku8TtSf11IUHiMVvNB6SjgAA1ENNq8hh/h46VVSq7MLSC473cXPUP+4c0JAhAgAAoAVLyizQ/qN5+mBHinbVMdEoSfYmk4aG+bH4jUZD0hEAgHq4b+XuaqvIZypSB/u6ydXRQfHH8mSpYeyALj76x7QB8najOiAAAADqpz5HqGsyNMxPi2/tY+WogPMj6QgAwAWcWU1+Z/Mh/Xo077z9lt3ZX4GeLpoZHVNlMtizg5deuKGXLg9q0wjRAgAAoCWqqSJ1Xbz4x14aFOrLDkc0OpKOAACcR2zaKf2/f/+qX9PPn2g8W3rOaXUL9NTK6QOVnFWolOzCyjsfAQAAgItVn4KFZ5w5Tn3LwM4NFBVQO5KOAACc42KPrpx9KXeIH8lGAAAAWEdyVmG9x3CcGrZG0hEAgHNERcdqWz0SjnYm6aowf5KMAAAAuCRnqlGfOS1jGIY2xWfq+W8OXHCsvcmkvl3a6MFRYZy2QZNA0hEAgLNczNGVfl18WEUGAADARavppE2foDZysDdpV0pFEUMHO5PMFkPGeT7jzM5GihaiqSDpCADA7wzD0L9jj9a5v50qEo6fzhjScEEBAACgxaupSEzM4RxJkpODne4eEqypg7ro6a9+rZKYHBDso2lDgnVZB292NqLJIekIAICkmLRTWrAmTj8ln6zzmKvC/dnhCAAAgEtyoZM2H9w9UIO7+koSBQvRrJB0BAC0ailZhXplXby+2ZchSXJ2sFOAp7PST52WpYb+rCYDAADgUpx7b+PBEwW19i8uN1d5TcFCNBckHQEArVJ2QYkWf5+oj35IVbnFkMkk3di3kx4d203uTg6aGR1TZcW5ZwcvvXBDL10e1MZ2QQMAAKDZ+uaXo3p5XZxSs09XtnUL9FR2YUmt44J9STCieSLpCABoVU6XmvXe9mQt3XRIBSXlkqSR3f315PgIRbb3quzH0RUAAABYQ2p2oSa/vV2nisqqvZdwPF9SxWmb0nJLlSIx9iaThob5MQ9Fs0XSEQDQKpgthj7fc0QLN8TreF7FanLPjl6aNyFSQ8P8ahzD0RUAAABcqvMlHM/28fRBevP7xConbc5UowaaK5KOAIAWzTAMbYrP1II1B5RwvOK+nI5tXPXE+O6adHkH2dmZbBwhAAAAWqrN8ScumHCUpPzSck7aoMUh6QgAaLF+OZKjF1Yf0A9JFRWpvV0dNXN0mO4Y3EXODvY2jg4AAAAt3Q/J2XXqd+beRk7aoCUh6QgAaHHSsov0yvp4/WfvUUmSk4Od7h4arAdHhMnbzdHG0QEAAKClObcidXGZWSt3puiDHam1jjNJGhbuT6IRLRJJRwBAs3Xu5O5UYakWf5+oD39IUZm5oiL1DX066rFx3dWxjautwwUAAEALk1NUqqjo2Cp3MXYP9FRecZkycoslSfYmyWzUPH5YuD/3NqLFIukIAGh2aprchfi6K6ugRPm/V6QeFu6nuRMidFkHb1uFCQAAgBYuKjpW2xOzqrTF/16Rup2Xix4b100DurTVDUurFpPxcLbXsqn9NSS85oKGQEtA0hEA0OzUNLlLzi6UJEW299K8CREa3s3fFqEBAACglUjKLKiyCH6u9+8eoIj2XpKkmGfGaevBTP2cdkp9O/toWDhzVbR8JB0BAM3KhSZ3b93WR139PRoxIgCNzTAMffbZZ9q4caNOnDghi8VS5f0vvvjCRpEBAFqTn9NO1fp+Rl5xZdJRqjhKTbIRrQlJRwBAs7I54fwJR0lKO1lE0hFo4R555BEtW7ZMo0aNUmBgoEwmk61DAgC0Irmny7RkU6Le25Zca78zFamB1oqkIwCgWTh8skgL18frq9ijtfZjcge0fB999JG++OILTZw40dahAABakZJysz7cmaq3NiYq5/f7Gb1dHZV/ukxn77m3N5k0NMyPitRo9Ug6AgCatJyiUr29MVEf7EhVqbliOhfg6ays/BImd0Ar5e3trdDQUFuHAQBoJSwWQ//55aheWRevI6dOS5LCAzw0b2KE+gb5KOqfVQscDg3zoyI1IJKOAIAmqrjMrJU7U/TW94nKK66oSD2kq6+emhipIB83zYyOYXIHtFLPPvus5s+fr/fee0+urq62DgcA0ILtSMzSC2sO6Nf0PElSoJezZo/tpj/17SQHeztJ0srpA5WcVaiU7EIF+7qzCA78jqQjAMDmkjILlHqySMG+7urS1k1fxaZr4foEpedUrCRHtPPU3AkRGtHNv/LuNiZ3QOt10003KTo6WgEBAQoODpajo2OV93/++WcbRQYAaI6SMgv031+O6mRhma6ODNCwcH/FHcvTi2vitCm+YpHbw9lBM0aE6p6rQuTmVD2VEuLHfBQ41yUlHQsKCqpVC/Ty8jpPbwAAqsopKlVUdNXjKO7ODiosqdjZ2N7bRY+N664b+nSUvV31QhFM7oDW6a677tKePXs0depUCskAAC5aTlGp7lu5W7tS/leF+v0dKXK0N6ncbMiQ5GBn0tQru2jm6DD5ejjbLligGap30jE5OVkPP/ywNm3apOLi4sp2wzBkMplkNputGiAAoOWKio7V9sSsKm2FJeWytzPp8XHddffQYLk42tsoOgBN1TfffKN169bpqquusnUoAIBmLCo6tkrC8YwysyFJurZXe825pruCWeQGLkq9k4633367JOm9995jZRkAcNGSMguq7HA8m9liaHzPdiQcAdQoKCiI0zUAgIuWlFmgH5NPnncuesYtA4NIOAKXoN5Jx19++UV79uxR9+7dGyKeC1qyZIleeeUVZWRk6LLLLtPrr7+uYcOG2SQWAMDFyS0q0yvr4mvtk5JdyNFpADVauHChnnjiCb3zzjsKDg62dTgAgGaipqt9avNz2ikNC/dv4KiAlsuuvgMGDBigw4cPN0QsF7Rq1SrNmjVLTz/9tGJiYjRs2DBNmDBBaWlpNokHAFA/JeVm/WNrkoa/slFrfj1Wa99gXxKOAGo2depUbdy4UV27dpWnp6fatm1b5Vd9LFmyRCEhIXJxcVG/fv20devWWvuXlJTo6aefVpcuXeTs7KyuXbvqvffeu5RvBwDQCDbHn9D1b23T1jomHCWpb2efBowIaPnqvdPxH//4h2bMmKH09HT17NmzWrXAyy+/3GrBnWvRokWaPn267r33XknS66+/rnXr1mnp0qVasGBBg31dAMClsVgM/eeXo3plXbyOnKqoSN0t0ENO9nbafzRPZ5ckszeZNDTMj12OAM7r9ddft8rnnFnQXrJkiYYOHap3331XEyZM0P79+9W5c+cax9x88806fvy4li9frrCwMJ04cULl5eVWiQcAYH2p2YWa/PZ2nSoqq9c4HzdHdjkCl6jeScfMzEwdOnRId999d2WbyWRq8EIypaWl2rNnj+bOnVulfdy4cdqxY0eNY0pKSlRSUlL5Oi8vr0FiAwCc347ELL2w5oB+Ta/4GRzo5azHxnbXn/p1UkFxuWZGx1Q54jI0zE+Lb+1jq3ABNAPTpk2zyufUd0F77dq12rx5s5KSkip3VHK8GwCatotJOHq7OurrhyhWBlyqeicd77nnHvXp00fR0dGNWkgmKytLZrNZgYGBVdoDAwN17FjNR/QWLFig+fPnN0Z4ANCqJWUWKPVkkYJ93St3KMYdy9OLa+K0Kb4ioejh7KAHRnbVPUND5OpUUSDG281RK6cPVHJWoVKyC6uMB4DamM1mffXVVzpw4IBMJpN69Oih6667Tvb2dStAdTEL2l9//bX69++vl19+WR9++KHc3d113XXX6a9//atcXV1rHMMiOADYztp9GXVOOE4fGiyzIV0dGcAOR8BK6p10TE1N1ddff62wsLCGiOeCzk1yntlhWZN58+Zp9uzZla/z8vIUFBTUoPEBQGtS02Xcg4LbKtDbRf/55agMQ3KwM2nqlV00c3SYfD2ca/ycED+SjQDqLjExURMnTlR6erq6d+8uwzCUkJCgoKAgffPNN+ratesFP+NiFrSTkpK0bds2ubi46Msvv1RWVpYefPBBnTx58rz3OrIIDgCNr7Tcouif0vTimrgL9j1ztc//TbqsESIDWpd6Jx1Hjx6tvXv3NnrS0c/PT/b29tUmgSdOnKg2WTzD2dlZzs41/wMXAHBpkjILFPXPGO0/WnXXzo8pJyv/+9rL22vOuO4KJqEIwIqioqLUtWtX/fDDD5XHnLOzszV16lRFRUXpm2++qfNn1WdB22KxyGQy6eOPP5a3t7ekiiPaN954o95+++0adzuyCA4ADefc0zaGYWjNr8f08to4pWQX1ekzuNoHaDj1TjpOmjRJjz76qPbt26devXpVKyRz3XXXWS24szk5Oalfv37asGGDbrjhhsr2DRs26Prrr2+QrwkAqK6m3Y01efv2vrq2V/tGigpAa7J58+YqCUdJ8vX11YsvvqihQ4fW6TMuZkG7ffv26tixY2XCUZIiIyNlGIaOHDmi8PDwamNYBAcA66tpPtq7k7cMQ/olPVeS5OfhpFljuunVdXHKOV294JeHs73+M3MYp22ABlTvpOOMGTMkSc8991y19xqykIwkzZ49W3fccYf69++vwYMHa9myZUpLS6uMCQDQsM63u7Embk51u1cNAOrL2dlZ+fn51doLCgrk5ORUp8+4mAXtoUOH6tNPP1VBQYE8PDwkSQkJCbKzs1OnTp0u4jsBAFyMqOhYbU/MqtK290hFstHV0V73Dw/VfcND5eHsoBHh/rru7W1V7nb0casoFBPk69aocQOtTb2TjhaLpSHiqJMpU6YoOztbzz33nDIyMtSzZ0+tXr1aXbp0sVlMANAa1HV349mCfVk1BtAw/vCHP+j+++/X8uXLNXDgQEnSjz/+qBkzZtTr1M2FFrTnzZun9PR0rVy5UpJ022236a9//avuvvtuzZ8/X1lZWZozZ47uueee8xaSAQBYV1JmQa1z0g+nD1T/4P/thA/ydVPMM+O09WCmfk47pb6dfSgUAzSSeicdbe3BBx/Ugw8+aOswAKBVqWk1+XzOXMbNURUADeXNN9/UtGnTNHjw4MqrfsrLy3Xdddfp9ddfr/PnXGhBOyMjQ2lpaZX9PTw8tGHDBs2cOVP9+/eXr6+vbr75Zj3//PNW/f4AAOcXd7z6Tvez5ZdUP0otScPC/Uk2Ao3MZBiGUd9Bmzdv1quvvqoDBw7IZDIpMjJSc+bM0bBhwxoiRqvJy8uTt7e3cnNz5eXlZetwAKBZSMos0OiFm+vcf3i4vxbf2kfebo4X7gygXpjLVJWYmKgDBw7IMAz16NGj0QsdXgyeIQBcnDKzRf/cdViL1sdXOSp9ro2Pj2TxG2hgdZ3P2NX3gz/66CONGTNGbm5uioqK0sMPPyxXV1ddffXV+uSTTy4paABA05JfXKZFGxIu2M9OUs8OXtr4+EitnD6QhCOABvXcc8+pqKhIYWFhmjRpkq677jqFhYXp9OnTNd47DgBovgzD0Npfj+ma17bo/776VaeKyuTiaF8tmWFvMml4uD8JR6AJqfdOx8jISN1///169NFHq7QvWrRIf//733XgwAGrBmhNrCwDQM2SMguUerJIwb7uCvFzV2m5RdE/pemN7w7qZGHpBcezuxFoHMxlKtjb2ysjI0MBAQFV2rOzsxUQENCghQ0vFc8QAOpuT+pJvbA6TntST0mSfN2d9MiYcE3s2V6z/7W3yt2OzEeBxlPX+Uy973RMSkrSpEmTqrVfd911euqpp+r7cQAAG6qpQExkOy8VlpYr7WSRJCnU310uDvaKy8jT2aXE7CT16OClxbf1ZUUZQKMyDEMmk6la+969e9W2bdsaRgAAmpOkzAK9vDZea387JklycbTTfcNCdf/wUHm6VCQVV04fqOSsQqVkF1YunANoWuqddAwKCtJ3331X7c6c7777TkFBQVYLDADQ8GoqEHPgWJ4kyc/DWY+ODdeU/kEqLDFrZnRMleTkVawmA2hkPj4+MplMMplM6tatW5XEo9lsVkFBQWXlaQBA05aUWaD//pKhU4WlGh0ZoGHh/srML9Eb3yUo+qfDMlsM2Zmkm/sH6dGx3RTo5VLtM0L8SDYCTVm9k46PPfaYoqKiFBsbqyFDhshkMmnbtm16//339cYbbzREjACABpCUWVAliXiulfcMUI8O3pIkbzc7VpMB2Nzrr78uwzB0zz33aP78+fL29q58z8nJScHBwRo8eLANIwQAXEhOUanu/WC3dv9+ZFqSVuxIkYuDnUwmk06XVVyRMSYyQE+Oj1B4oKetQgVwieqddHzggQfUrl07LVy4UP/6178kVdzzuGrVKl1//fVWDxAA0DD2Hsmp9f3j+SXqcU4bq8kAbGnatGmSpJCQEA0dOlQODvWeygIAbOyBj36uknA8o7i84iKf3p28NW9ipK4M9W3s0ABY2UXN1G644QbdcMMN1o4FANBAzi4U4+/prGVbkrRs86FaxwT7klwE0DR5enrqwIED6tWrlyTp3//+t1asWKEePXro2WeflZOTk40jBADUJCmzQDuTsmvt89i4biQcgRai3knH0NBQ7dq1S76+VX8I5OTkqG/fvkpKSrJacACAS1NToRhHe5PKzIYkydPFQYXF5VUKxNibTBoa5seORgBN1p///GfNnTtXvXr1UlJSkqZMmaI//vGP+vTTT1VUVKTXX3/d1iECAGrwY/LJC/aJOZyj4d0CGiEaAA3Nrr4DUlJSZDabq7WXlJQoPT3dKkEBAKwjKjpW2865t7HMbMjV0V5Lb++rrXNG6apw/yrvDw3z0+Jb+zRmmABQLwkJCbriiiskSZ9++qlGjBihTz75RO+//74+//xz2wYHAJBUsatxY/wJJWcVSpKSswr18Y+pFxzXt7NPQ4cGoJHUeafj119/Xfnf69atq3Jxt9ls1nfffafg4GCrBgcAqJ+zj1EbhnHeQjGny8yKaO+lNu5OFIgB0OwYhiGLpWKP9rfffqs//OEPkqSgoCBlZWXZMjQAaPVqOmnTwdtFJ/JLVG4xah3r5eKgYecsiANovuqcdJw8ebIkyWQyVV7ifYajo6OCg4O1cOFCqwYHAKibmiZ3Xi61/4hPyS6sTDBSIAZAc9K/f389//zzGjNmjDZv3qylS5dKkpKTkxUYGGjj6ACgdYuKjtX2xKoLQEdziyVJo7r762RhqfYeya02zsHOpG9mDmuUGAE0jjonHc+sJoeEhGjXrl3y8/NrsKAAAPVT0+Qur7i81jEUigHQXL3++uu6/fbb9dVXX+npp59WWFiYJOmzzz7TkCFDbBwdALReSZkF5z1pI0nPTLpMbd2cNDM6pkq/Xh289NG9V8rbzbExwgTQSOpdSCY5OblaW05Ojtq0aWONeAAAdVDXY9SSZGeSzj7JQqEYAM3d5Zdfrn379lVrf+WVV2Rvb2+DiACg9Tl7PhriVzEn/Xrv0VrHnDlpw/U+QOtQ76TjSy+9pODgYE2ZMkWSdNNNN+nzzz9X+/bttXr1avXu3dvqQQIAKpzvjpza9OjgpV/T8ypfUygGQEvl4lL7z0MAwKWraT7aJ6iN7O1M2p16qtaxZ5+04XofoOWrd9Lx3Xff1UcffSRJ2rBhg7799lutXbtW//rXvzRnzhytX7/e6kECACpWk6P+GaP9R/OqtJ+5I+d8Ft/aV5JYSQbQYtjZ2clkMp33fbPZ3IjRAEDrUtO1PjGHcyRJTg52CvB01tFTp2U5631O2gCtU72TjhkZGQoKCpIk/fe//9XNN9+scePGKTg4WIMGDbJ6gADQ2tW0mlwTk6Sz6wGeO7ljkgegpfjyyy+rvC4rK1NMTIw++OADzZ8/30ZRAUDLd6E7G1feM1CR7byq3dnISRugdap30tHHx0eHDx9WUFCQ1q5dq+eff16SZBgGq8oA0ABqWk2uyWUdOUYNoHW4/vrrq7XdeOONuuyyy7Rq1SpNnz7dBlEBQMtz7r2NB08U1Nr/dJlZ3m6O3NkIQNJFJB3/+Mc/6rbbblN4eLiys7M1YcIESVJsbGxl5UAAwMXbHH9CsUdy1Lezjzq2cb3gDsczOEYNoLUbNGiQ7rvvPluHAQDNXk0nbboFeiq7sKTWcdzZCOBs9U46vvbaawoODtbhw4f18ssvy8PDQ1LFsesHH3zQ6gECQGuRml2oyW9v16misso2J/vz31l2BseoAUA6ffq0Fi9erE6dOtk6FABo9mo6aZNwPF+S5Oxgp9JyS63X+gCAdBFJR0dHRz3++OPV2mfNmmWNeACgVUrKLNB1b21TQUnVaypKzcZ5RvwPx6gBtDY+Pj5VCskYhqH8/Hy5ublVFjwEAFycC93b+PH0QXrz+0TubARwQfVOOq5cubLW9++8886LDgYAWpuLLRJjJ6lHBy8tvq0vK8oAWp3XXnutStLRzs5O/v7+GjRokHx8fGwYGQA0f7tTT9X6fn5pOXc2AqiTeicdH3nkkSqvy8rKVFRUJCcnJ7m5uZF0BIB6qGuRmM6+bkrNLqp8fVW4vxbf2kfebo4NGR4ANEl33XWXrUMAgBYnp6hUb32fqA92ptTa78y9jdzZCOBC6p10PHWq+qrHwYMH9cADD2jOnDlWCQoAWqqkzAL9mHxSJkkd2rjUuUjM85N7qpOPG6vJACBpxYoV8vDw0E033VSl/dNPP1VRUZGmTZtmo8gAoPkpLjPrgx0pentjovKKyyVJbVwdlXe6TJaz+nFvI4D6qnfSsSbh4eF68cUXNXXqVMXFxVnjIwGgRdl7+JSe/Hyf4o7l13usj5ujhoX7S6JIDABI0osvvqh33nmnWntAQIDuv/9+ko4AUAcWi6EvY9K1cH28juYWS5Ii2nlq3sRI9e7orah/xnJvI4BLYpWkoyTZ29vr6NGj1vo4AGgR6npn4/n4uDnq64eusnJUANC8paamKiQkpFp7ly5dlJaWZoOIAKB52ZKQqQVr4nQgI0+S1MHbRbPHddcNfTrK3q7izlzubQRwqeqddPz666+rvDYMQxkZGXrrrbc0dOhQqwUGAM1ZUmaBUk8WacnGRO25wGXc57IzSUFt3fT85J6VOxwBAP8TEBCgX375RcHBwVXa9+7dK19fX9sEBQDNwG9Hc/XimjhtPVhxp7ini4MeGhWmu4YEy8XRvlp/7m0EcCnqnXScPHlyldcmk0n+/v4aPXq0Fi5caK24AKBZutSdjZJ0VRhFYgCgNrfccouioqLk6emp4cOHS5I2b96sRx55RLfccouNowMA2zr7DvFBob4K8XPXkVNFWrQ+QV/GpsswJEd7k+4cHKyHR4XJx93J1iEDaKHqnXS0WCwX7gQArVRdq1Gf68PpA1VuMTi6AgB18Pzzzys1NVVXX321HBwqprMWi0V33nmnXnjhBRtHBwC2kVNUqgc++lk7k7KrtHdo46Ks/FKVmiv+LX9d7w6ac013BbV1s0WYAFoRq93pCACt1eb4E4o9kqP2Xq4XtcNxeLg/x6gBoB6cnJy0atUqPf/884qNjZWrq6t69eqlLl262Do0ALCZqOjYaglHSTqaU1EkZnCor+ZNjNDlndo0cmQAWqt6Jx1vvPFG9e/fX3Pnzq3S/sorr+inn37Sp59+arXgAKApS80u1OS3t+tUUdlFf8aQrr5UAQSAixQeHq7w8PDzvu/l5aXY2FiFhoY2YlQA0PiSMgsuuPj9txt6KtTfo5EiAoCLSDpu3rxZf/nLX6q1jx8/Xq+++qpVggKApi4ps0DXvbVNBSXmeo8ND3TXbQO7aGT3AI5SA0ADMgzD1iEAQKNIPVlUpz4kHQE0pnonHQsKCuTkVP2iWUdHR+Xl5VklKABoqi6mUIy9yaS+ndvowdFh3NkIAAAAqzqQkaelGxMv2C/YlzkogMZlV98BPXv21KpVq6q1//Of/1SPHj2sEhQANFUXUyhmaJif/jFtgEaxsxEAAABWcjTntB7/dK8mvrlVP6WckqmWvsPD/ZmHAmh09d7p+H//93/605/+pEOHDmn06NGSpO+++07R0dHc5wigxUnKLFDqySIF+7rLMIw673B85cbL5efpzM5GAAAAWFXu6TIt3XRIK7Ynq6S8oiL1tZe31wPDu+pvqw9UKybDHeIAbKXeScfrrrtOX331lV544QV99tlncnV11eWXX65vv/1WI0aMaIgYAaDR1XSM2tvVsU5jfdwcdVP/oIYKDQBQRyZTbft+AKB5KSk366Mf0rT4+4PK+b2Q4cCQtnpqYqSuCGojSYq+/0olZxXqh6RsmSQNCvVlARyAzdQ76ShJ1157ra699tpa+0RHR+u6666Tuzs/4AA0PzUdo849feEq1T5ujvr6oasaKiwAQD1QSAZAS2CxGPrvvgy9si5Oh0+eliSFBXho3oQIjY4IqLbAEuLHSRsATcNFJR3r4s9//rMGDRqk0NDQhvoSAGA19TlGbSfJctZrk6TObV31/A29NCzcv6FDBQD8rrS0VMnJyeratascHKpPa9esWaOOHTvaIDIAsI4dh7L04po4/XIkV5IU4Oms2WO76cZ+neRgX+8SDQDQqBos6cjKMoDmoKZj1J18XGsd06ODl349mlf5eli4vxbf2kfebnU7fg0AuDRFRUWaOXOmPvjgA0lSQkKCQkNDFRUVpQ4dOmju3LmSpKuuYuc5gOYp/li+Xlobp+/jTkiS3J3sNWNEV00fFiI3pwb7ZzwAWBU/rQC0ajUdoz5y6nStYxbf1leSlJJdSKEYALCBefPmae/evdq0aZPGjx9f2T5mzBj95S9/qUw6AkBTd/ZpmxA/dx3LLdaiDfH6bM8RWQzJwc6k2wZ1VtTV4fLzcLZ1uABQLyQdAbRaSZkFtR6jNkk6e8+2vcmkoWF+lUlGko0AYBtfffWVVq1apSuvvLLKXWY9evTQoUOHbBgZANRNTadtgnzclJlfrOLfK1JP7NVOc66JYM4JoNki6Qig1fo57VSt7192zjHqoWF+Wnxrn4YOCwBwAZmZmQoICKjWXlhYSMVqAM1CTadtDp8qkiQNCPbR3AmR6tfFxxahAYDVkHQE0CqcfXSlrbuTlm46pOXbkmodwzFqAGiaBgwYoG+++UYzZ86UpMpE49///ncNHjzYlqEBwAVd6LTNS3+6XKH+Ho0YEQA0jAZLOnbp0kWOjhRVAGBbNR1dcbAzqdxScXDa29VR+afLqlSj5hg1ADRtCxYs0Pjx47V//36Vl5frjTfe0G+//aadO3dq8+bNtg4PAGq1fv/xWt9PPVlE0hFAi2BX3wGhoaHKzs6u1p6Tk6PQ0NDK17/++quCgoIuLToAuAhJmQXaGH9CyVmFioqO1bZzVpLLLYbcnOy1fFp/bX58pK4K96/yPseoAaBpGzJkiLZv366ioiJ17dpV69evV2BgoHbu3Kl+/frZOjwAqNHB4/m694NdenFNXK39gn1Z8AbQMtR7p2NKSorMZnO19pKSEqWnp1slKAC4GDXtajyfolKzQv091MbdSSunD1RyViHHqAGgGenVq5c++OADW4cBABd0PK9Yr3+boFW7DstiSPZ2JgV4Out4bnGtp20AoLmrc9Lx66+/rvzvdevWydvbu/K12WzWd999p+DgYKsGBwD1UdOF3LVJyS6scoSaCR4ANA9paWm1vt+5c+dGigQA/ufsO8RD/NxVUFKuZZsP6e9bk3W6rGLjzjWXBeqJ8RHyc3fWzOiYKovlnLYB0NLUOek4efJkSRUXdU+bNq3Ke46OjgoODtbChQutGhwA1NWFLuSuCUdXAKB5Cg4OrrVKdU2ncgCgodR02ibM30PZhSU6VVQmSerbuY2emhip/sFtK/tw2gZAS1fnpKPFUrHxOyQkRLt27ZKfn1+DBQUA9XXgWF6d+3J0BQCat5iYmCqvy8rKFBMTo0WLFulvf/ubjaIC0FrVdNomMbNAUsVpmifHd9c1l7WrcbGE0zYAWrJ6F5JJTk4m4QigySgtt2jF9mTN+3xfncdwdAUAmrfevXtX+dW/f3/dd999evXVV/Xmm2/W67OWLFmikJAQubi4qF+/ftq6dWudxm3fvl0ODg664oorLuI7ANBSnDltYzaMGt9fdkc/je/Zvtbd2QDQUtW7kIwkFRYWavPmzUpLS1NpaWmV96KioqwSGADUxjAMfbMvQy+vjVfaySJJkqujvUrKzDVeyD3/+ss4ugIALVy3bt20a9euOvdftWqVZs2apSVLlmjo0KF69913NWHCBO3fv7/WeyFzc3N155136uqrr9bx48etETqAZuqHpOxa3z+Sc1rhgZ6NFA0ANC0mwzjPksx5xMTEaOLEiSoqKlJhYaHatm2rrKwsubm5KSAgQElJSQ0V6yXLy8uTt7e3cnNz5eXlZetwAFykH5KytWBNnPYezpEk+Xs669Ex3XTNZYF6dNXeKvfpDA/31+Jb+8jbzdFG0QKA9TCXqZCXV/VKDcMwlJGRoWeffVZxcXGKjY2t0+cMGjRIffv21dKlSyvbIiMjNXnyZC1YsOC842655RaFh4fL3t5eX331VZ2/3pnYeYZA83civ1hvfHtQ0T+lyVLLv6g3Pj6SBW8ALU5d5zP13un46KOPatKkSVq6dKnatGmjH374QY6Ojpo6daoeeeSRSwoaAGpz8Hi+Xlobp28PnJAkuTvZ6/7hXXXvsBC5O1f8OONCbgBo+dq0aVPtqKJhGAoKCtI///nPOn1GaWmp9uzZo7lz51ZpHzdunHbs2HHecStWrNChQ4f00Ucf6fnnn7/g1ykpKVFJSUnl63MTpgCal8KScv19a5KWbUlSUWlF0SpfdyedKiyt8bQNc1EArVm9k46xsbF69913ZW9vL3t7e5WUlCg0NFQvv/yypk2bpj/+8Y8NESeAVux4XrFe25Cgf+0+LIsh2duZdNvAzoq6Olz+ns7V+nMhNwC0bBs3bqzy2s7OTv7+/goLC5ODQ92mt1lZWTKbzQoMDKzSHhgYqGPHjtU45uDBg5o7d662bt1a56+zYMECzZ8/v059ATRd5WaLVu0+rNc2HFRWQcVCQu+gNnpqQoQi2nlpZnRMldM23CEOABeRdHR0dKxcWQ4MDFRaWpoiIyPl7e2ttLQ0qwcIoPXKLy7Tsi1J+vvWJBWXVawdj7+sneaM766u/h42jg4AYCsjRoyw2mfVtGOypoIPZrNZt912m+bPn69u3brV+fPnzZun2bNnV77Oy8tTUFDQxQcMoFEZhqH1+4/rpbVxSsoslCR18XXTE9dEaGKv/1Wk5rQNAFRX76Rjnz59tHv3bnXr1k2jRo3SM888o6ysLH344Yfq1atXQ8QIoJUpM1sU/VOa3vj2oLILK4pV9evio3kTItQ/uK2NowMA2MLXX39d577XXXfdBfv4+fnJ3t6+2q7GEydOVNv9KEn5+fnavXu3YmJi9PDDD0uSLBaLDMOQg4OD1q9fr9GjR1cb5+zsLGfn6rvyATR9e1JPacHqA9qdekqS1NbdSY9cHa5bB3aWk4Ndtf6ctgGAquqddHzhhReUn58vSfrrX/+qadOm6YEHHlBYWJjee+89qwd4xt/+9jd98803io2NlZOTk3JychrsawGwDcMwtPbXY3p5XbySsypWkkP93PXE+Ahdc1lgjTtPAACtw+TJk+vUz2QyyWw2X7Cfk5OT+vXrpw0bNuiGG26obN+wYYOuv/76av29vLy0b9++Km1LlizR999/r88++0whISF1ig9A05eUWaBX1sVrza8VixIujna696pQ/XlEqDxdKE4IAHVV76Rj//79K//b399fq1evtmpA51NaWqqbbrpJgwcP1vLlyxvlawJoPLtSTuqF1QcUk5YjSfLzcNIjY7rplgFBcrSvvpIMAGhdLBbLhTvV0+zZs3XHHXeof//+Gjx4sJYtW6a0tDTNmDFDUsXR6PT0dK1cuVJ2dnbq2bNnlfEBAQFycXGp1g6gecoqKNGb3x3UJz+mqdxiyM4k3dQvSI+O7aZ23i62Dg8Amp16Jx1Hjx6tL774Qm3atKnSnpeXp8mTJ+v777+3VmxVnLmA+/3332+QzwdgG4knCvTy2jit339ckuTqaK/7hofq/uGh8nCu948oAADqbMqUKcrOztZzzz2njIwM9ezZU6tXr1aXLl0kSRkZGdxZDrQCRaXl+sfWZL27+ZAKf69IPToiQE+Oj1D3dp42jg4Ami+TYRhGfQbY2dnp2LFjCggIqNJ+4sQJdezYUWVlZVYN8Fzvv/++Zs2aVafj1SUlJSopKal8febi7tzcXHl5eTVglAAu5ER+sd749qD+ueuwzBZD9nYmTRkQpFlXhyvAi5VkAKhJXl6evL29mctIKiws1ObNm5WWlqbS0tIq70VFRdkoqgvjGQJNR7nZok/3HNFrGxJ0Ir/i342Xd/LW3AkRGtLVz8bRAUDTVdf5TJ23Ef3yyy+V/71///4ql26bzWatXbtWHTt2vMhwG8aCBQsqd0gCaBoKS8orK1IX/b6SPLZHoJ4c311hAawkAwAuLCYmRhMnTlRRUZEKCwvVtm1bZWVlyc3NTQEBAU066QjA9gzD0HcHTujFtXFKPFEgSQpq66o510ToD73ay86Oe8QBwBrqnHS84oorZDKZZDKZaqzM5+rqqsWLF9friz/77LMXTAru2rWryj2S9TFv3jzNnj278vWZnY4AGl+Z2aJVuw7r9W8PKqugYiX5iqA2empipAaGUJEaAFB3jz76qCZNmqSlS5eqTZs2+uGHH+To6KipU6fqkUcesXV4AJqw2MM5emH1Af2UfFKS1MbNUVGjw3X7lZ3l7GBv4+gAoGWpc9IxOTlZhmEoNDRUP/30k/z9/Svfc3JyUkBAgOzt6/dD+uGHH9Ytt9xSa5/g4OB6febZnJ2d5ezsfNHjAVw6wzC0fv9xvbQ2TkmZFRWpg33d9MT4CE3o2Y6K1ACAeouNjdW7774re3t72dvbq6SkRKGhoXr55Zc1bdo0/fGPf7R1iACamNTsQr28Ll7f/JIhSXJ2sNM9V4Voxoiu8nalIjUANIQ6Jx3PXKhtzcqBfn5+8vPjrgygpdqTekoLVh/Q7tRTkiRfdydFXR2uWwd2lpMDFakBABfH0dGxctEqMDBQaWlpioyMlLe3N4VfAFSRXVCixd8n6uMfU1VmNmQySX/q20mzx3ZThzautg4PAFq0epeG/eCDD+Tn56drr71WkvTEE09o2bJl6tGjh6KjoyuTk9aWlpamkydPKi0tTWazWbGxsZKksLAweXh4NMjXBHBxkjIL9PLaeK39reLuVxdHO903rKIitacLK8kAgEvTp08f7d69W926ddOoUaP0zDPPKCsrSx9++KF69epl6/AANAGnS816b3uylm46pIKScknSyO7+enJ8hCLbU8QJABpDvatXd+/eXUuXLtXo0aO1c+dOXX311Xr99df13//+Vw4ODvriiy8aJNC77rpLH3zwQbX2jRs3auTIkXX6DKoFAg0rM79Eb353UJ/8lCazxZCdSbq5f5Bmjemmdt5UpAaAS8VcpsLu3buVn5+vUaNGKTMzU9OmTdO2bdsUFhamFStWqHfv3rYO8bx4hkDDMlsMfb7niBZuiNfxvIp7xHt29NK8CZEaGsYpOwCwhrrOZ+qddHRzc1NcXJw6d+6sJ598UhkZGVq5cqV+++03jRw5UpmZmZccfENhkgc0jKLScv1ja7Le3XxIhb9XpL46IkBPTohQt0AqUgOAtTCXaf54hkDDMAxDm+Iz9eKaOMUfz5ckdWzjqifGd9ekyztQkRoArKiu85l6H6/28PBQdna2OnfurPXr1+vRRx+VJLm4uOj06dMXHzGAZqfcbNGne45o0YYEZeZXrCT37uStuRMiNbirr42jAwC0VPPnz9fUqVPVtWtXW4cCoAn45UiOFqyO086kbEmSt6ujZo4O0x2Du1CRGgBsqN5Jx7Fjx+ree+9Vnz59lJCQUHm342+//XZJlaYBNB+GYejbAyf00to4JZ4okCR1buumJ8Z317W92lORGgDQoD7//HM999xzGjBggKZOnaopU6bI39/f1mEBaGRp2UV6ZX28/rP3qCTJycFOdw8J1oMjw+Ttxj3iAGBr9U46vv322/p//+//6fDhw/r888/l61uxm2nPnj269dZbrR4ggKYlJu2UFqyO008pJyVJPm6Oiro6XLcP6kJFagBAo/jll1/022+/6eOPP9aiRYs0e/ZsjRkzRlOnTtXkyZPl5uZm6xABNKBThaVa/H2iPvwhpbIi9Q19Ouqxcd3VkYrUANBk1PtOx7p68MEH9dxzz8nPr+lc1ssdOsDFS8kq1Cvr4vXNvgxJkrODnaZfFaIZI7vKi4rUANAomMvUbPv27frkk0/06aefqri4WHl5ebYO6bx4hsDFKy4za8X2FC3ZlKj84oqK1MPC/TR3QoQu6+Bt4+gAoPVosDsd6+qjjz7S448/3qSSjgAuLCmzQKknixTs664QP3dlF5Ro8feJ+uiHVJVbKlaSb+zbSbPHdVN7b1aSAQC25+7uLldXVzk5OSk/P9/W4QCwMrPF0Jcx6Vq4Pl4ZucWSpMj2Xpo3IULDu3G1AgA0VQ2WdGygDZQAGkhOUamiomO15eD/KtAH+7opM7+ksiL1yO7+mjshQhHt2JkBALCt5ORkffLJJ/r444+VkJCg4cOH69lnn9VNN91k69AAWIlhGNqcUFGROu5YxYJCB28XPX5Nd02+oiMVqQGgiWuwpCOA5iUqOlbbE7OqtKVkF0mSenb00lMTIjUkjJ3LAADbGzx4sH766Sf16tVLd999t2677TZ17NjR1mEBsKJf03O1YM0BbU+sqEjt6eKgh0eFadqQYLk4UpEaAJoDko4AlJRZUGWH47neuKWPuvp7NGJEAACc36hRo/SPf/xDl112ma1DAWBlh08WaeH6eH0V+3tFans7TRvSRQ+NClMbNycbRwcAqA+SjgC0KeH8CUdJSjtZRNIRANBkvPDCC3Xq5+XlpdjYWIWGhjZwRAAuVU5Rqd7emKgPdqSq1GyRJE2+ooMeG9ddQW2pSA8AzVG9k45paWkKCgqSyVT1/gzDMHT48GF17tzZasEBsJ5zC8RIUlp2kV5ZH6//7D1a69hgX/fGCBEAAKvijnGg6SsuM2vlzhS99X2i8n6vSD2kq6/mTYhUr05UpAaA5qzeSceQkBBlZGQoICCgSvvJkycVEhIis7mi4MTUqVNrLZsNoHHUVCBmcKivuvq7a9XuwyozV1Sk9vdwVlZ+iSxnjbU3mTQ0zK8ySQkAAABYg8Vi6N970/XqugSl55yWJEW089TcCREa0c2/2iYXAEDzU++ko2EYNf4FUFBQIBcXl8rXS5cuvbTIAFhFTQVidiZla2dSxaXcw7v5a+74CHVs46qZ0TFVkpNDw/y0+NY+jRovAAAAWrZtB7P0wuoD2p+RJ0lq7+2i2WO76Y99O8meitQA0GLUOek4e/ZsSZLJZNL//d//yc3tf/dqmM1m/fjjj7riiiusHiCAi3ehAjGv3Hi5buofVPl65fSBSs4qVEp2YZVj2AAAAMCl2n80TwvWHNDWgxUL4p7ODnpwVJjuHkpFagBoieqcdIyJiZFUsdNx3759cnL6X+UwJycn9e7dW48//rj1IwRw0VKyC2t938/TuVpbiB/JRgBAy8DxTKBpSM85rYXr4/VlTLoMQ3K0N+mOK4P18OgwtXWnIjUAtFR1Tjpu3LhRknT33XfrjTfe4L5GoIn7NT1Xb32fWGsfCsQAAJq7M8ViakowUkgGaFznFi7MPV2mJZsStWJ7ikrLK24On9S7g+aM667OvlSkBoCWrt53Oq5YsaIh4gBgJYdPFmnh+nh9FVtRkdpkkmRIZ/+ziwIxAIDmbvny5Xrttdd08OBBSVJ4eLhmzZqle++9t7LPmjVr1LFjR1uFCLQaOUWlum/lbu1KOVXZFurnruzCUuWeLpMkDQppq6cmRqp3UBsbRQkAaGz1TjoCaJpyikr11veJWrkzVaXmipXkG/p01H3DQvXimjgKxAAAWoz/+7//02uvvaaZM2dq8ODBkqSdO3fq0UcfVUpKip5//nlJ0lVXXWXLMIFWIaeoVKNe3aRTRWVV2pOyKq756RboobkTIjSqewBXHgBAK2MyWtG5k7y8PHl7eys3N5fj4WgxisvM+mBHit7emKi84nJJ0tAwX82bEKmeHb0r+1EgBgCaP+YyFfz8/LR48WLdeuutVdqjo6M1c+ZMZWVl2SiyC+MZoqW5cekO7U49dd73v509QmEBHo0YEQCgodV1PsNOR6CZslgMfRWbroXrE5Sec1qSFNHOU/MmRmp4uF+1lWQKxAAAWgqz2az+/ftXa+/Xr5/Ky8ttEBHQOiVlFtSacJSkw6eKSDoCQCtF0hFoBs69lHvrwUwtWB2n/Rl5kqT23i56bFx33dCno+ztOLYCAGjZpk6dqqVLl2rRokVV2pctW6bbb7/dRlEBrUtG7mk98/VvF+xH4UIAaL1IOgJNWE5RqaKiY6vcx+jj5lh5Z46ni4MeHBmmu4cGy8XR3lZhAgDQ6JYvX67169fryiuvlCT98MMPOnz4sO68807Nnj27st+5iUkAlyavuEzvbDqk5duSVfJ7RerzGRDsw0kbAGjFSDoCTczZuxr/8u/ftD2x6r1Up4rKZJJ0z1UhenhUmHzcnWwTKAAANvLrr7+qb9++kqRDhw5Jkvz9/eXv769ff/21sh9FKwDrKS236OMfU/XmdwcrF8AHBrdVqdmiXw7n6Nz0o4+bo/5x54DGDxQA0GSQdASaiJp2NZ6PIWnqlV1IOAIAWqWNGzfaOgSg1TAMQ//9JUOvrItX2skiSVJYgIeeHB+hMZEByjtdrpnRMVXmsAO6+Ogf0wbI283RVmEDAJoAko5AExEVHVttV2NtUrILOa4CAACABvNDUrYWrD6gvUdyJUn+ns6aPbabburXSQ72dpIkbzdHrZw+UMlZhUrJLqy8gxwAAJKOQBOQlFlQpx2OZ+NSbgAAADSEhOP5emlNnL6LOyFJcney159HdNW9w0Lk5lTzPyFD/Eg2AgCqIukINAGpvx9VqQt7k0lDw/yY1AEAAMCqjucV67UNCfrX7sOyGJK9nUm3DeysqKvD5e/pbOvwAADNDElHwMbijuVp6aZDde4/NMxPi2/t04ARAQAAoDXJLy7Tu5uT9I9tSSouqygJM/6ydpozvru6+nvYODoAQHNF0hGwkYzc01q4PkGf/3xEhiGdqa9pnNXnzK7G+ddfxh05AAAAsKrScouif0rTG98d1MnCUklSvy4+empihPp1aWvj6AAAzR1JR6CR5RWXaemmQ3pvW7JKyitWkq+9vL1mDO+qV9bFV7nb8cyuRm83R5KNAAAAsArDMLTm12N6eW2cUrIrrvkJ9XPXkxMiNK5HoEwm0wU+AQCACyPpCDSS0nKLPvohVYu/P6hTRWWSpIEhbTVvQoT6dPaRJCr/AQAAoEH9lHxSL6w+oNjDOZIkPw8nzRrTTVMGBMnx94rUAABYA0lHoIFZLIa+2ZehV9bFK+33gjFhAR6aOz5CV0cGVFtJpvIfAAAArC3xRL5eWhuvDfuPS5LcnOx137BQ3Tc8VB7O/LMQAGB9/O0CNKCdh7L14poD2nskV5IU4Oms2WO76cZ+neTASjIAAAAa2Im8Yr327UGt2pVWWZF6yoAgzbo6XAFeLrYODwDQgpF0BBpAwvF8vbgmTt/HnZAkuTvZa8aIrpo+LERuTvyxAwAAgHUlZRYo9WRR5RU9BSXlWrYlSX/fkqTTZWZJ0rgegXpifITCAqhIDQBoeGQ/ACs6llus1zYk6NM9h2UxJAc7k24b1FlRV4fLz8PZ1uEBAACghckpKlVUdGyVYoRh/h46WVRaWZG6T+c2empipAYEU5EaANB4SDoCVpBfXKZ3NyfpH9uSVFxWUZF6Qs92mnNNd4X6s5IMAACAhhEVHavtiVlV2hIzCyRJwb5uenJ8hMb3bEdFagBAoyPpCFyC0nKLon9K0xvfHaxcSe7fxUfzJkaqXxcfG0cHAACAliops0A/JmdX2eF4rr/f2V/hgZ6NGBUAAP9D0hG4CIZhaPW+Y3plXZxSsisqUof6u2vu+AiN7RHISjIAAAAaRE3Hqc/nSM5pko4AAJsh6QjU00/JJ/XC6gOKPZwjSfLzcNajY8M1pX8QFakBAADQoKKiY7WtDglHSQr2dW/gaAAAOD+SjkAdJZ7I14tr4vXtgeOSJDcne90/PFT3DQuVuzN/lAAAAGB9Z1elLiopr9MOR3uTSUPD/BTiR9IRAGA7ZEqACziRV6zXvj2oVbvSZDEkezuTbhkQpEfGhCvA08XW4QEAAKCFScos0G8ZeVq5I0W7Uk5VtjvY1e0Kn6Fhflp8a5+GCg8AgDoh6QicR0FJuZZtSdLftyTpdJlZkjSuR6CeGB+hsAAqUgMAAMC6LnRfY7nFqHX8gj/20pWhvuxwBAA0CVxAB5yjzGzRhztTNPKVjXrzu4M6XWZW385t9NmMwVp2Z38SjgAAtCBLlixRSEiIXFxc1K9fP23duvW8fb/44guNHTtW/v7+8vLy0uDBg7Vu3bpGjBYtXVR0rLYnZl2w37n/iLM3mTQ83F+3DuxMwhEA0GSQdAR+ZxiG1v6aoWte26L/+/dvyiooVYifu96Z2lefPzBE/YPb2jpEAABgRatWrdKsWbP09NNPKyYmRsOGDdOECROUlpZWY/8tW7Zo7NixWr16tfbs2aNRo0Zp0qRJiomJaeTI0RIlZRZoy8FMmY3adzNKUo8OXlVec5waANAUmQyjDn+rtRB5eXny9vZWbm6uvLy8LjwArcbulIqK1D+n5UiSfN2dNGtMuG4Z2FmOVKQGADQRzGWsa9CgQerbt6+WLl1a2RYZGanJkydrwYIFdfqMyy67TFOmTNEzzzxTp/48Q5xxdoGYED93/Ts2XY/8M7ZOYzc+PlKSlJJdWDkeAIDGUtf5DHc6olU7lFmgl9fGad1vFRWpXR3tdd+wEN0/oqs8qEgNAECLVVpaqj179mju3LlV2seNG6cdO3bU6TMsFovy8/PVti2nIVB3Nd3bGOzrpuN5JRcce25VapKNAICmjKwKWqXM/BK98V2Con86LLPFkJ1JmjIgSLPGdFOgFxWpAQBo6bKysmQ2mxUYGFilPTAwUMeOHavTZyxcuFCFhYW6+eabz9unpKREJSX/Sybl5eVdXMBoMWq6tzElu0iS5OHsoKKSclnOM5Zj1ACA5oSkI1qVwpJy/X1rkpZtSVJRaUVF6jGRgXpyfHeFB3raODoAANDYTCZTldeGYVRrq0l0dLSeffZZ/fvf/1ZAQMB5+y1YsEDz58+/5DjRMpy5t/F8Prp3oBatP1ilz4AuPrprSLB6dPRmZyMAoFkh6YhWodxs0ardh/XahoPKKqjYbdA7qI2emhChQaG+No4OAAA0Nj8/P9nb21fb1XjixIlqux/PtWrVKk2fPl2ffvqpxowZU2vfefPmafbs2ZWv8/LyFBQUdPGBo1nbGH+i1vdPFZVp5fSBSs4q5L5GAECzR9IRLZphGNqw/7heWhunQ5mFkqQuvm564poITezVrk47GQAAQMvj5OSkfv36acOGDbrhhhsq2zds2KDrr7/+vOOio6N1zz33KDo6Wtdee+0Fv46zs7OcnZ2tEjOar9TsQr2yLl7//SWj1n7Bvv+7q5FkIwCguSPpiBbr57RTWrD6gHalnJIktXV3UtToMN02qIucHKhIDQBAazd79mzdcccd6t+/vwYPHqxly5YpLS1NM2bMkFSxSzE9PV0rV66UVJFwvPPOO/XGG2/oyiuvrNwl6erqKm9vb5t9H2i6ThaWavH3B/XRD6kqMxsymaQAD2dl5pdUubfx3AIxAAC0BCQd0eIkZxXqlXVxWr2v4h8CLo52uveqUP15RKg8XRxtHB0AAGgqpkyZouzsbD333HPKyMhQz549tXr1anXp0kWSlJGRobS0tMr+7777rsrLy/XQQw/poYceqmyfNm2a3n///cYOH01YcZlZ721P1tKNh5RfUi5JGt7NX3PHR6hjG1fNjI6pcm8jBWIAAC2RyTAMw9ZBNJa8vDx5e3srNzdXXl5etg4HVpZVUKLF3x3Uxz+mqfz3itQ39QvSo2O7qZ03FakBAM0fc5nmj2fYsiRlFij1ZFHl3Ytmi6HPfz6i1zYkKCO3WJLUo72X5k2M0LBw/ypjubcRANBc1XU+w05HNHtFpeVavjVZ72w+pMLfK1KPjgjQk+Mj1L0dFakBAABgXTlFpYqKjq2yW7FnBy8Vl1uUeKJAktSxjavmXNNd1/XuIDu76veIc28jAKClI+mIZqvcbNFne45o0YYEncivqEh9eSdvzZ0QoSFd/WwcHQAAAFqipMwCRf0zRvuP5lVp//X3114uDpo5Olx3DO4iF0d7W4QIAECTQNIRzcLZR1eCfd30fdwJvbgmTgd/X0kOauuqOddE6A+92te4kgwAAABcipp2N9bkw+mD1DuoTeMEBQBAE0bSEU3a3sOn9PSXv1auHEuSt6ujck+XSZLauDlq5uhwTb2ys5wdWEkGAABAw4iKjtX2xKwL9jtZVNoI0QAA0PSRdESTVNtKcu7pMtmZpD+P6KoZI7rK25WK1AAAALC+M6dt7E2mC+5wPCPYl3saAQCQmknSMSUlRX/961/1/fff69ixY+rQoYOmTp2qp59+Wk5OTrYODw0gKjpW2xLPP7GzGNLN/YNIOAIAAMDq6nqU+mz2JpOGhvlRHAYAgN81i6RjXFycLBaL3n33XYWFhenXX3/Vfffdp8LCQr366qu2Dg9WlpRZUKcJXkp2IZM6AAAAWF1dj1KfbWiYnxbf2qeBIgIAoPlpFknH8ePHa/z48ZWvQ0NDFR8fr6VLl5J0bGHMFkMf/Zhap74cXQEAAIA1nF200DCMOu9wtJPUo4OXFt/Wl8VwAADO0SySjjXJzc1V27Zta+1TUlKikpKSytd5eXm19IYtGYahTQmZenF1nOKP59fa107SVeH+TOwAAABwSWo6Rt2mHtf3XBXur8W39pG3G1f+AABwrmaZdDx06JAWL16shQsX1tpvwYIFmj9/fiNFhYu170iuFqw5oB2HsiVVVKf2dXdSSlahLDX0PzO5AwAAAC5FTceoc06X1Trmw+kDVW4xFOzrziI4AAC1sGnS8dlnn71gUnDXrl3q379/5eujR49q/Pjxuummm3TvvffWOnbevHmaPXt25eu8vDwFBQVdWtCwmsMni/TKunh9vfeoJMnJwU53DwnWgyPDJEkzo2OqrDr37OilF27opcs7tbFFuAAAAGgB/leRWrUeo7YzVRQvPONMoZhh4f6NECUAAM2fyTAM48LdGkZWVpaysmq/oDk4OFguLi6SKhKOo0aN0qBBg/T+++/Lzs6uXl8vLy9P3t7eys3NlZeX10XHjUtzqrBUb21M1Ic7U1Vqtshkkm64oqNmj+umTj5uVfomZxUqJbuQlWQAAMRcpiXgGdpOfStS9+zopV/T/3c903COUgMAIKnu8xmb7nT08/OTn59fnfqmp6dr1KhR6tevn1asWFHvhCNsr7jMrPd3pOjtjYnKLy6XJA0L99PcCRG6rIN3jWNC/Eg2AgAA4NLVtyL14lv7ShIL4AAAXKRmcafj0aNHNXLkSHXu3FmvvvqqMjP/tzrZrl07G0aGujBbDH0Zk65F6+N1NLdYkhTZ3kvzJkRoeDeOpwAAAKBhJWUW1HmH45lj1GeSjCQbAQC4OM0i6bh+/XolJiYqMTFRnTp1qvKeDU+How42J2RqweoDijtWUZG6g7eLHr+muyZf0VF2diYbRwcAAICW7mjOaT3z79/q3H9omB9FCwEAsIJmkXS86667dNddd9k6DNTDr+m5enFNnLb9foTF08VBD48K07QhwXJxtLdxdAAAAGjpck+XaemmQ1qxPVkl5ZZa+1KRGgAA62sWSUc0H0dOFWnh+gR9GZMuSXKyt9Odg7vooVFh8nF3snF0AAAAaOlKys366Ic0Lf7+oHKKyiRJA0PaqrTcon1HcmU+66QUFakBAGg4JB1hFblFZXp7U6Le356iUnPFSvLkKzrosXHdFdTW7QKjAQAAgEtjsRj6774MvbIuTodPnpYkhQd4aO6ECI2OCFDe6XLNjI6pcrcjR6kBAGg4JB1xSYrLzPpwZ6re2pio3NMVK8lDuvpq3oRI9epUc0VqAAAAwJp2HMrSi2vi9MuRXElSgKezZo/tphv7dZKDvZ0kydvNUSunD1RyViEVqQEAaAQkHXFRLBZD/96brlfXJSg9p2IlOaKdp+ZOiNCIbv4ymSgSAwAAgIYVfyxfL62N0/dxJyRJHs4O+vPwUE0fFiI3p5r/qRPiR7IRAIDGQNIR9bbtYJYWrDmg347mSZLaebnosXHd9Me+nWRPRWoAAAA0sGO5xVq0IV6f7TkiiyE52Jl0+6DOmnl1uPw8nG0dHgAAEElH1MP+o3l6cW2ctiRU3IPj6eygB0Z11T1DQ6hIDQAAgAaXV1ymdzcf0vJtySouq7hHfGKvdppzTQS7FwEAaGJIOuKCjuac1sL1Cfoi5ogMQ3K0N+mOK4P18OgwtaUiNQAAABpYablFn/yYqje/T9TJwlJJ0oBgH82bGKm+nX1sHB0AAKgJSUecV+7pMi3ZlKgV21NUWl6xkjypdwfNGdddnX2pSA0AAICGZRiGVu87ppfXxSk1u0iSFOrvrrnjIzS2RyD3iAMA0ISRdEQ1JeX/q0idU1RRkXpQSFs9NTFSvYPa2DY4AAAAtAo/JmXrhTVx2ns4R5Lk51FRkfrm/v+rSA0AAJouko6oZLEY+s8vR/XKungdOVVRkTo8wEPzJkZoVPcAVpIBAADQ4A4er6hI/e2BiorUbk72+vPwrrp3WIjcnfnnCwAAzQV/a0OStONQlhasjtO+9FxJUqBXxUryn/qykgwAAICGdzyvWK9/m6BVuw7LYkj2dibdOjBIj1zdTf6eVKQGAKC5IenYysUdy9NLa+K0Mb6iIrWHs4NmjAjVPVeFyM2J/z0AAADQsApKyrVs8yH9fWuyTpeZJUnXXBaoJ8ZHqKu/h42jAwAAF4usUiuVkXtai9Yn6LOfKypSO9iZNPXKLpo5Oky+HqwkAwAAoGGVmS2K/ilNb3x7UNm/V6Tu27mNnpoYqf7BbW0cHQAAuFQkHVuZvOIyvbv5kJZvS1ZxWUVF6mt7tdeca7or2M/dxtEBAACgpTMMQ2t/PaaX18UrOatQkhTq564nxnfXNZe14x5xAABaCJKOrURpuUUf/5iqN787qFO/V6QeGNxWcydGqG9nHxtHBwAAgNZgV8pJLVh9QD+n5UiS/Dyc9MiYbrplQJAcuUccAIAWhaRjC2cYhr7Zl6GX18Yr7WSRJKmrv7vmTojUmEgqUgMAAKDhJZ4o0Mtr47R+/3FJkqujve4bHqr7h4fKg4rUAAC0SPwN34L9kJStBWvitPdwjiTJ39NZj47pppv7U5EaAAAADe9EfrHe+Pag/rnrsMwWQ3YmacqAznp0TLgCvFxsHR4AAGhAJB1boIPH8/XS2jh9e+CEJMndyV5/HtFV9w6jIjUAAAAaXmFJuf6+NUnLtiSpqLSiIvWYyEDNndBdYQGeNo4OAAA0BjJQLcjxvGK9tiFB/9p9WBZDsrcz6baBnRV1dbj8PalIDQAAAOtKyizQj8nZkky6MtRXnXxctWrXYb3+7UFlFZRIknoHtdFTEyI0KNTXtsECAIBGRdKxBcgvLtOyLUn6+9akyorU4y9rpznju6urv4eNowMAAEBLk1NUqgc//lk7DmVXaXdxtKucj3bxddOT4yM0oScVqQEAaI1IOjZjZWaLon9K0xvfHlR2YakkqV8XHz01MUL9urS1cXQAAABoqaKiY6slHCWpuMwiR3uT/t+1PXTrwM5ycuAecQAAWiuSjs2QYRha++sxvbwuXslZhZKkUD93PTkhQuN6BLKSDAAAgAaTlFmgLQczz/t+mdnQ8G7+JBwBAGjlSDo2M7tSTuqF1QcUk5YjSfLzcNKsMd00ZUCQHKlIDQAAgAb2S3ruBfukZBcqxM+9EaIBAABNFUnHZiLxRIFeWhunDfuPS5JcHe11//BQ3Tc8VB7OPEYAAAA0rKLScv1ja7KWbjp0wb7BviQcAQBo7chWNXEn8ov1+rcHtWrXYZkthuztTJoyIEizrg5XgJeLrcMDAABAC1dutujTPUf02oYEncivqEjt4eyggpLyGvsPD/dnlyMAACDp2FQVlpRXVqQuKjVLksb2CNST47srLMDTxtEBAACgpTMMQ98dOKEX18Yp8USBJCmoraueuCZCV4X56aFPqlevHhzqq8W39rFFuAAAoIkh6djElJktWrXrsF7/9qCyCipWkvt0bqOnJkZqQDAVqQEAANDwYg/n6IXVB/RT8klJko+bo2aODtftV3aWs4O9JOmT+65UclahfkzKliHpylBfdjgCAIBKJB2bCMMwtH7/cb20Nk5JmRUVqYN93fTk+AiN79mOitQAAABocKnZhXp5Xby++SVDkuTsYKd7rgrRjBFd5e3qWK1/iJ87iUYAAFAjko5NwJ7Uk3phdZz2pJ6SJPm6O+mRMeG6dWBnKlIDAACgwWUXlGjx94n6+MdUlZkNmUzSjX076dGx3dShjautwwMAAM0QSUcbSsos0Mtr47X2t2OSJBdHO903LFT3Dw+Vp0v1lWQAAADAmk6XmvXe9oqK1GcKw4zs7q8nx0cosr2XjaMDAADNGUlHG8jML9Gb3x3UJz+lyWwxZGeSbu4fpEfHdlMgFakBAABgZUmZBUo9WaRg34rj0GaLoc/3HNHCDfE6nldxj3jPjl6aNyFSQ8P8bBwtAABoCUg6NqKi0nL9Y2uy3t18SIW/V6QeExmgJ8dHKDyQitQAAACwrpyiUkVFx2rLwczKtp4dvFRcbqmsSN2xjaueGN9dky7vIDs77hEHAADWQdKxEZSbLfrX7iN67dsEZeZXrCT37uSteRMjdWWor42jAwAAQEsVFR2r7YlZVdp+PZonSfJ2ddTM0WG6Y3CXyorUAAAA1kKVEitLyizQxvgTSs4qlGEY2rD/uK55fYue+nKfMvNL1Lmtm966rY++emgoCUcAAAAbW7JkiUJCQuTi4qJ+/fpp69attfbfvHmz+vXrJxcXF4WGhuqdd95ppEjrLymzQFsOZspsGDW+/+H0gbp3WCgJRwAA0CDY6WglNR1d8XJxVF5xmSTJx81RUVeH6/ZBXeTkQK4XAADA1latWqVZs2ZpyZIlGjp0qN59911NmDBB+/fvV+fOnav1T05O1sSJE3Xffffpo48+0vbt2/Xggw/K399ff/rTn2zwHdTut993NJ5PdmFpI0UCAABaI5NhnGfpswXKy8uTt7e3cnNz5eVl3Wp8dy7/SdsTs6qtJNuZpBkjumrGyK7yoiI1AAC4BA05l2mNBg0apL59+2rp0qWVbZGRkZo8ebIWLFhQrf+TTz6pr7/+WgcOHKhsmzFjhvbu3audO3fW6Ws2xjMsLjNrxfYUvbXxoApLzOftt/HxkQrxc2+QGAAAQMtV1/kMOx2t4MzRlZpYDOmm/kEkHAEAAJqQ0tJS7dmzR3Pnzq3SPm7cOO3YsaPGMTt37tS4ceOqtF1zzTVavny5ysrK5OhYfb5XUlKikpKSytd5ebXvPrwUZouhL2PStXB9vDJyiyVJ7k72Ol1qluWsfvYmk4aG+ZFwBAAADYpzvpfodKlZb29MrLVPSnZhI0UDAACAusjKypLZbFZgYGCV9sDAQB07dqzGMceOHauxf3l5ubKysmocs2DBAnl7e1f+CgoKss43cBbDMLQp/oSufXOrHv90rzJyi9XB20WLbu6tbU+O1lXh/lX6Dw3z0+Jb+1g9DgAAgLOx0/EimS2GPttzWIs2JOh4XkmtfYN9WUUGAABoikwmU5XXhmFUa7tQ/5raz5g3b55mz55d+TovL8+qicdf03O1YM0BbU/MliR5uTjo4dFhunNwsFwcKwrErJw+UMlZhUrJLlSwrzs7HAEAQKMg6VhPhmFoY/wJvbgmTgnHCyRJnXxc5eHsoIRj+RxdAQAAaAb8/Pxkb29fbVfjiRMnqu1mPKNdu3Y19ndwcJCvr2+NY5ydneXs7GydoM9y+GSRFq6P11exRyVJTvZ2mjakix4aFaY2bk7V+of4kWwEAACNi6RjPew9nKMFaw7oh6STkqQ2bo56eFSY7hjcRcWlFs2MjqlytyNHVwAAAJomJycn9evXTxs2bNANN9xQ2b5hwwZdf/31NY4ZPHiw/vOf/1RpW79+vfr371/jfY4NIaeoVG9vTNQHO1JVaq5Y7p58RQc9Nq67gtq6NUoMAAAAdUHSsQ7Ssov0yvp4/Wfv7yvJDna6Z2iIHhjZVd6uFRNMZwd7jq4AAAA0I7Nnz9Ydd9yh/v37a/DgwVq2bJnS0tI0Y8YMSRVHo9PT07Vy5UpJFZWq33rrLc2ePVv33Xefdu7cqeXLlys6OrrBYy0uM2vlzhS99X2i8orLJUlDuvrqqYmR6tnRu8G/PgAAQH2RdKzFqcJSLf4+UR/+kKIysyGTSfpjn06aPa6bOrZxrXEMR1cAAACahylTpig7O1vPPfecMjIy1LNnT61evVpdunSRJGVkZCgtLa2yf0hIiFavXq1HH31Ub7/9tjp06KA333xTf/rTnxosRovF0L/3puvVdQlKzzktSYpo56m5EyI0opt/rfdPAgAA2JLJOHP7dSuQl5cnb29v5ebmysvL67z9isvMem97spZuOqT831eSh3fz19zxEerR4fzjAAAAGlJd5zJouurzDLcezNSC1XHan5EnSWrv7aLHxnXXDX06yt6OZCMAALCNus5n2Ol4FrPF0Bc/H9GiDQnKyC2WJPVo76V5EyM0LNzfxtEBAACgNdh/NE8L1hzQ1oNZkiRPZwc9OCpMdw/9X0VqAACApo6koyoqUm9OyNSLa+IUdyxfktSxjavmXNNd1/XuIDtWkgEAANDA0nNOa+H6eH0Zky7DkBztTbrjymA9PDpMbd2rV6QGAABoylp90vHX9FwtWHNA2xOzJUleLg6aOTpcdwzuwkoyAAAAGlzu6TIt2ZSoFdtTVFpeUZF6Uu8OmjOuuzr7UpEaAAA0T6026Xj4ZJEWro/XV7G/V6S2t9NdQ4P14MiuauPGSjIAAAAaVkm5WR/uTNVbGxOVU1QmSboytK2emhipyzu1sW1wAAAAl6hVJh1fXhunf+3NVqm5YiX5hj4dNXtsNwW1ZSUZAAAADe+bX47q7e1HdeRURUXqboEemjchUiO7U5EaAAC0DK0y6bhyZ6rsnN00NMxX8yZEqmdHb1uHBAAAgFbkyc/3yc7ZTYFeznpsbHf9qV8nKlIDAIAWpVUlHQ3DkCSFetvpiesiNbSrr0wmk/Ly8mwcGQAAwIWdmbOcmdOg+Tnz7FxUovuvCtMdVwbL1clehQX5No4MAACgbuo6JzUZrWjWeuTIEQUFBdk6DAAAgEty+PBhderUydZh4CIwHwUAAC3FheakrSrpaLFYdPToUXl6enJXThOVl5enoKAgHT58WF5eXrYOBw2M59168KxbD551wzIMQ/n5+erQoYPs7OxsHQ4uQmPNR/mz2LzwvJoXnlfzwvNqXnhezUNd56St6ni1nZ0duwKaCS8vL37AtCI879aDZ9168Kwbjrc3d1E3Z409H+XPYvPC82peeF7NC8+reeF5NX11mZOyRA4AAAAAAADAqkg6AgAAAAAAALAqko5oUpydnfWXv/xFzs7Otg4FjYDn3XrwrFsPnjXQNPBnsXnheTUvPK/mhefVvPC8WpZWVUgGAAAAAAAAQMNjpyMAAAAAAAAAqyLpCAAAAAAAAMCqSDoCAAAAAAAAsCqSjgAAAAAAAACsiqQjmqyUlBRNnz5dISEhcnV1VdeuXfWXv/xFpaWltg4NDeBvf/ubhgwZIjc3N7Vp08bW4cCKlixZopCQELm4uKhfv37aunWrrUNCA9iyZYsmTZqkDh06yGQy6auvvrJ1SECLV9+fr5s3b1a/fv3k4uKi0NBQvfPOO40UKaT6Pa8vvvhCY8eOlb+/v7y8vDR48GCtW7euEaPFxc5ftm/fLgcHB11xxRUNGyCqqO/zKikp0dNPP60uXbrI2dlZXbt21XvvvddI0aK+z+vjjz9W79695ebmpvbt2+vuu+9WdnZ2I0WLS0HSEU1WXFycLBaL3n33Xf3222967bXX9M477+ipp56ydWhoAKWlpbrpppv0wAMP2DoUWNGqVas0a9YsPf3004qJidGwYcM0YcIEpaWl2To0WFlhYaF69+6tt956y9ahAK1CfX++Jicna+LEiRo2bJhiYmL01FNPKSoqSp9//nkjR9461fd5bdmyRWPHjtXq1au1Z88ejRo1SpMmTVJMTEwjR946Xez8JTc3V3feeaeuvvrqRooU0sU9r5tvvlnfffedli9frvj4eEVHRysiIqIRo2696vu8tm3bpjvvvFPTp0/Xb7/9pk8//VS7du3Svffe28iR42KYDMMwbB0EUFevvPKKli5dqqSkJFuHggby/vvva9asWcrJybF1KLCCQYMGqW/fvlq6dGllW2RkpCZPnqwFCxbYMDI0JJPJpC+//FKTJ0+2dShAi1Xfn69PPvmkvv76ax04cKCybcaMGdq7d6927tzZKDG3Ztb4+/Cyyy7TlClT9MwzzzRUmPjdxT6vW265ReHh4bK3t9dXX32l2NjYRogW9X1ea9eu1S233KKkpCS1bdu2MUOF6v+8Xn31VS1dulSHDh2qbFu8eLFefvllHT58uFFixsVjpyOaldzcXP5iAJqJ0tJS7dmzR+PGjavSPm7cOO3YscNGUQFA83cxP1937txZrf8111yj3bt3q6ysrMFihXX+PrRYLMrPz2ce3Agu9nmtWLFChw4d0l/+8peGDhFnuZjn9fXXX6t///56+eWX1bFjR3Xr1k2PP/64Tp8+3Rght2oX87yGDBmiI0eOaPXq1TIMQ8ePH9dnn32ma6+9tjFCxiVysHUAQF0dOnRIixcv1sKFC20dCoA6yMrKktlsVmBgYJX2wMBAHTt2zEZRAUDzdzE/X48dO1Zj//LycmVlZal9+/YNFm9rZ42/DxcuXKjCwkLdfPPNDREiznIxz+vgwYOaO3eutm7dKgcH/ondmC7meSUlJWnbtm1ycXHRl19+qaysLD344IM6efIk9zo2sIt5XkOGDNHHH3+sKVOmqLi4WOXl5bruuuu0ePHixggZl4idjmh0zz77rEwmU62/du/eXWXM0aNHNX78eN10003c3dCMXMyzRstjMpmqvDYMo1obAKD+6vvztab+NbWjYVzs34fR0dF69tlntWrVKgUEBDRUeDhHXZ+X2WzWbbfdpvnz56tbt26NFR7OUZ8/XxaLRSaTSR9//LEGDhyoiRMnatGiRXr//ffZ7dhI6vO89u/fr6ioKD3zzDPas2eP1q5dq+TkZM2YMaMxQsUlYhkGje7hhx/WLbfcUmuf4ODgyv8+evSoRo0apcGDB2vZsmUNHB2sqb7PGi2Ln5+f7O3tq61anjhxotrqJgCg7i7m52u7du1q7O/g4CBfX98GixWX9vfhqlWrNH36dH366acaM2ZMQ4aJ39X3eeXn52v37t2KiYnRww8/LKkiqWUYhhwcHLR+/XqNHj26UWJvjS7mz1f79u3VsWNHeXt7V7ZFRkbKMAwdOXJE4eHhDRpza3Yxz2vBggUaOnSo5syZI0m6/PLL5e7urmHDhun5559np34TR9IRjc7Pz09+fn516puenq5Ro0apX79+WrFihezs2JzbnNTnWaPlcXJyUr9+/bRhwwbdcMMNle0bNmzQ9ddfb8PIAKB5u5ifr4MHD9Z//vOfKm3r169X//795ejo2KDxtnYX+/dhdHS07rnnHkVHR3N3WSOq7/Py8vLSvn37qrQtWbJE33//vT777DOFhIQ0eMyt2cX8+Ro6dKg+/fRTFRQUyMPDQ5KUkJAgOzs7derUqVHibq0u5nkVFRVVu7bA3t5e0v927KMJM4AmKj093QgLCzNGjx5tHDlyxMjIyKj8hZYnNTXViImJMebPn294eHgYMTExRkxMjJGfn2/r0HAJ/vnPfxqOjo7G8uXLjf379xuzZs0y3N3djZSUFFuHBivLz8+v/HMryVi0aJERExNjpKam2jo0oEW60M/XuXPnGnfccUdl/6SkJMPNzc149NFHjf379xvLly83HB0djc8++8xW30KrUt/n9cknnxgODg7G22+/XWUOnJOTY6tvoVWp7/M611/+8hejd+/ejRQt6vu88vPzjU6dOhk33nij8dtvvxmbN282wsPDjXvvvddW30KrUt/ntWLFCsPBwcFYsmSJcejQIWPbtm1G//79jYEDB9rqW0A9kHREk7VixQpDUo2/0PJMmzatxme9ceNGW4eGS/T2228bXbp0MZycnIy+ffsamzdvtnVIaAAbN26s8c/wtGnTbB0a0GLV9vN12rRpxogRI6r037Rpk9GnTx/DycnJCA4ONpYuXdrIEbdu9XleI0aM4GeqjdX3z9fZSDo2vvo+rwMHDhhjxowxXF1djU6dOhmzZ882ioqKGjnq1qu+z+vNN980evToYbi6uhrt27c3br/9duPIkSONHDUuhskw2I8KAAAAAAAAwHq4IA8AAAAAAACAVZF0BAAAAAAAAGBVJB0BAAAAAAAAWBVJRwAAAAAAAABWRdIRAAAAAAAAgFWRdAQAAAAAAABgVSQdAQAAAAAAAFgVSUcAaKXuuusuTZ48uc79N23aJJPJpJycnAaLCQAAoLkwmUz66quvbB1Gi1bf3+P6zm8BNCySjgCsbuTIkZo1a9Ylf05KSopMJpNiY2PrPdaak8D3339fbdq0qfc4knQAAABoDBc7/yZJB6AhkXQEAAAAAAAAYFUkHQFY1V133aXNmzfrjTfekMlkkslkUkpKynn7nzp1Srfffrv8/f3l6uqq8PBwrVixQpIUEhIiSerTp49MJpNGjhwpSdq1a5fGjh0rPz8/eXt7a8SIEfr5558rPzM4OFiSdMMNN8hkMlW+rs3evXs1atQoeXp6ysvLS/369dPu3bu1adMm3X333crNza38fp599llJ0kcffaT+/fvL09NT7dq102233aYTJ05IqtilOWrUKEmSj4+PTCaT7rrrLknSZ599pl69esnV1VW+vr4aM2aMCgsLK3//Jk+erBdeeEGBgYFq06aN5s+fr/Lycs2ZM0dt27ZVp06d9N5771WJf9++fRo9enTlZ95///0qKCiofN9sNmv27Nlq06aNfH199cQTT8gwjCqfYRiGXn75ZYWGhsrV1VW9e/fWZ599dt7fs9TUVE2aNEk+Pj5yd3fXZZddptWrV1/w9xoAAMAaRo4cqYcfflgPP/xw5Rzn//2//1dtjnOuefPm6corr6zWfvnll+svf/mLpAvPN89V0wmX2NjYanPhHTt2aPjw4XJ1dVVQUJCioqIq54EXsmTJEoWHh8vFxUWBgYG68cYbJZ1//m02mzV9+nSFhITI1dVV3bt31xtvvFH5ec8++6w++OAD/fvf/64ct2nTJpWWlurhhx9W+/bt5eLiouDgYC1YsKBynMlk0rvvvqs//OEPcnNzU2RkpHbu3KnExESNHDlS7u7uGjx4sA4dOlQl/qVLl6pr165ycnJS9+7d9eGHH1Z5/+DBgxo+fLhcXFzUo0cPbdiwodrvQXp6uqZMmSIfHx/5+vrq+uuvr/XfGrXNuwE0AgMArCgnJ8cYPHiwcd999xkZGRlGRkaGUV5eft7+Dz30kHHFFVcYu3btMpKTk40NGzYYX3/9tWEYhvHTTz8Zkoxvv/3WyMjIMLKzsw3DMIzvvvvO+PDDD439+/cb+/fvN6ZPn24EBgYaeXl5hmEYxokTJwxJxooVK4yMjAzjxIkTF4z7sssuM6ZOnWocOHDASEhIMP71r38ZsbGxRklJifH6668bXl5eld9Pfn6+YRiGsXz5cmP16tXGoUOHjJ07dxpXXnmlMWHCBMMwDKO8vNz4/PPPDUlGfHy8kZGRYeTk5BhHjx41HBwcjEWLFhnJycnGL7/8Yrz99tuVnzlt2jTD09PTeOihh4y4uDhj+fLlhiTjmmuuMf72t78ZCQkJxl//+lfD0dHRSEtLMwzDMAoLC40OHToYf/zjH419+/YZ3333nRESEmJMmzat8vt76aWXDG9vb+Ozzz6r/D3z9PQ0rr/++so+Tz31lBEREWGsXbvWOHTokLFixQrD2dnZ2LRpk2EYhrFx40ZDknHq1CnDMAzj2muvNcaOHWv88ssvxqFDh4z//Oc/xubNmy/4ew0AAGANI0aMMDw8PIxHHnnEiIuLMz766CPDzc3NWLZsWa3j9u3bZ0gyEhMTK9t+/fXXynmbYVx4vmkYhiHJ+PLLLw3DqD5PMgzDiImJMSQZycnJhmEYxi+//GJ4eHgYr732mpGQkGBs377d6NOnj3HXXXdd8HvdtWuXYW9vb3zyySdGSkqK8fPPPxtvvPGGYRjnn3+XlpYazzzzjPHTTz8ZSUlJlb8/q1atMgzDMPLz842bb77ZGD9+fOW4kpIS45VXXjGCgoKMLVu2GCkpKcbWrVuNTz75pMr33bFjR2PVqlVGfHy8MXnyZCM4ONgYPfr/t3f/MVXVfRzA3/y8EFwQEdAGwvDKDYgfGm0qGRYYZWtsNcsgEE0LSUDI0lapaTXWZlrM5mAuNQmKiYt0Ck3MXcUIgSsFCJgIlBiynA4FA+/n+cNxHo5cuNgD+jzP3q/tbpzz/Z7v+ZzzB3vve+/5niflyJEj0tDQIHPmzJGnn35aOaa4uFjs7Oxkx44d0tTUJFu3bhUbGxspLy8XEZFbt27Jww8/LAsWLJDa2lo5fvy4zJo1S3WPr1+/LjNnzpTly5dLXV2dNDQ0SHx8vOj1erl586aI3M7Sg/nWUu4moonHSUciGndRUVGSkZExpr7PPfecLFu2zGxba2urAJDa2tpRxxgYGBCtVivff/+9sm9oQBkLrVYru3fvNtv25Zdfiqurq8UxBidJB4OMufBZXV0tAOTChQtmx1i6dKn4+vrKrVu3lH16vV7mz5+vbA8MDIiTk5MUFBSIiEhubq64ublJT0+P0ufQoUNibW0tly5dEhGRadOmSXZ2ttLe398v3t7eSijr6ekRBwcHqaioUNXz6quvyssvv2z2ekJCQmTTpk0W7wsRERHRRIiKipLAwEAxmUzKvnXr1klgYKDFY0NDQ2Xz5s3K9jvvvCOPPvroiP0t5c2xTDomJibKa6+9phrXYDCItbW19Pb2jlrv/v37xcXFRTXpOdRY83dqaqq88MILyvbQSbpBaWlp8uSTT6ru61AA5L333lO2T506JQBk165dyr6CggJxcHBQtufNmycrV65UjbN48WJZtGiRiIiUlpaKjY2NdHR0KO2HDx9W3eNdu3aJXq9X1XXz5k1xdHSU0tLSYddjKXcT0cTj49VEdF+tWrUKhYWFCA8Px9tvv42KigqLx3R1dSElJQUBAQFwdXWFq6srenp60N7e/o/ryMrKwooVKxATE4Ps7Oxhj4OYU1tbi7i4OPj6+kKr1SqPf49WR1hYGKKjoxESEoLFixcjLy8PV65cUfUJDg6GtfW//z17eXkhJCRE2baxsYG7u7vyKHdjYyPCwsLg5OSk9ImMjITJZEJTUxOuXr2Kzs5OzJ07V2m3tbVFRESEst3Q0IC+vj4sXLgQzs7Oymfv3r0j3ov09HR8+OGHiIyMxMaNG1FXV2fxnhERERGNpzlz5sDKykrZnjt3LlpaWnDr1q1Rj0tISEB+fj6A20vMFBQUICEhQWmfiLxZXV2N3bt3q7JWbGwsTCYTWltbRz124cKF8PX1hb+/PxITE5Gfn48bN25YPOfOnTsREREBDw8PODs7Iy8vz+I1JCcnw2g0Qq/XIz09HWVlZcP6hIaGKn97eXkBgCqvenl5oa+vD9euXQNwO69GRkaqxoiMjERjY6PSPn36dHh7eyvtQ7MrcPv+nTt3DlqtVrl/kydPRl9fn9m8OpbcTUQTi5OORHRfPfPMM2hra8OaNWtw8eJFREdHY+3ataMek5ycjOrqamzfvh0VFRUwGo1wd3fH33///Y/r2LRpE+rr6/Hss8+ivLwcQUFBOHDgwIj9r1+/jqeeegrOzs7Yt28fqqqqlP6j1WFjY4MffvgBhw8fRlBQEHJycqDX61VB087OTnWMlZWV2X0mkwnA7aA8NGzf2W8sBsc6dOgQjEaj8mloaBhxXccVK1bg/PnzSExMxC+//IKIiAjk5OSM6XxERERE91N8fDyam5tRU1ODiooKdHR0YMmSJUr73ebNwS+MZch6kv39/ao+JpMJr7/+uiprnTlzBi0tLZgxY8ao9Wq1WtTU1KCgoADTpk3Dhg0bEBYWplpD8k7ffvstMjMzsXz5cpSVlcFoNGLZsmUWM/Ps2bPR2tqKLVu2oLe3Fy+++KKyfuSgodl0MG+a2zeYMYfuGzQ0w4qZdTjv7G8ymfDII4+o7p/RaERzczPi4+OHHT+W3E1EE4uTjkQ07uzt7S1+uzyUh4cHkpOTsW/fPmzfvh25ubnKOACGjWUwGJCeno5FixYhODgYGo0G3d3dqj52dnZ3VQMABAQEIDMzE2VlZXj++eeVF9qYu56zZ8+iu7sb2dnZmD9/Ph566CHll4eDRqrfysoKkZGR+OCDD1BbWwt7e/tRJzgtCQoKgtFoVC2KffLkSVhbWyvfzk+bNg0//fST0j4wMIDq6mrVGBqNBu3t7dDpdKqPj4/PiOf28fFBSkoKiouL8eabbyIvL+8fXwcRERHR3Rqabwa3Z86cCRsbm1GP8/b2xuOPP478/Hzk5+cjJiZG+cUeMLa8OZSHhwcAoLOzU9lnNBpVfWbPno36+vphWUun0ym5cTS2traIiYnBJ598grq6Oly4cAHl5eUAzOdVg8GAefPmITU1FbNmzYJOpxv2i8CRcruLiwteeukl5OXl4ZtvvsH+/fvx119/WaxxJIGBgThx4oRqX0VFBQIDAwHczqLt7e24ePGi0n7q1ClV/9mzZ6OlpQWenp7D7p+rq6vZ84537iaiu8NJRyIad35+fqisrMSFCxfQ3d2t+obzThs2bMB3332Hc+fOob6+HgcPHlTCh6enJxwdHXHkyBH8+eefuHr1KgBAp9Phq6++QmNjIyorK5GQkABHR8dhNRw9ehSXLl2y+BhFb28vVq9ejR9//BFtbW04efIkqqqqlDr8/PzQ09ODo0ePoru7Gzdu3MD06dNhb2+PnJwcnD9/HiUlJdiyZYtqXF9fX1hZWeHgwYO4fPkyenp6UFlZiY8//hinT59Ge3s7iouLcfnyZeVc/0RCQgIcHBywdOlS/Prrrzh27BjS0tKQmJiohOeMjAxkZ2fjwIEDOHv2LFJTU1XfjGu1WqxduxaZmZnYs2cPfvvtN9TW1mLHjh3Ys2eP2fOuWbMGpaWlaG1tRU1NDcrLy/+j6yAiIiK6Wx0dHcjKykJTUxMKCgqQk5ODjIyMMR2bkJCAwsJCFBUV4ZVXXlG1jSVv3tnfx8cHmzZtQnNzMw4dOoStW7eq+qxbtw6nTp3CG2+8AaPRiJaWFpSUlCAtLc1irQcPHsTnn38Oo9GItrY27N27FyaTCXq9HoD5/K3T6XD69GmUlpaiubkZ77//PqqqqlTj+vn5oa6uDk1NTeju7kZ/fz+2bduGwsJCnD17Fs3NzSgqKsLUqVMxadKkMd1Xc9566y3s3r0bO3fuREtLCz799FMUFxcrTzjFxMRAr9cjKSkJZ86cgcFgwLvvvqsaIyEhAVOmTEFcXBwMBgNaW1tx/PhxZGRk4Pfffx92zonI3UR0l+7ripJE9H+pqalJ5syZI46OjqrFs83ZsmWLBAYGiqOjo0yePFni4uLk/PnzSnteXp74+PiItbW1REVFiYhITU2NREREiEajkZkzZ0pRUZH4+vrKtm3blONKSkpEp9OJra2t+Pr6jlrvzZs3ZcmSJeLj4yP29vby4IMPyurVq1ULeqekpIi7u7sAkI0bN4qIyNdffy1+fn6i0Whk7ty5UlJSMuzFN5s3b5apU6eKlZWVLF26VBoaGiQ2NlY8PDxEo9FIQECA5OTkKP3NLeZtbmHwO6+3rq5OnnjiCXFwcJDJkyfLypUrVW/m6+/vl4yMDHFxcZFJkyZJVlaWJCUlqc5lMpnks88+E71eL3Z2duLh4SGxsbHKG6nvXCB99erVMmPGDNFoNOLh4SGJiYnS3d096r0mIiIiGi9RUVGSmpoqKSkp4uLiIm5ubrJ+/foRX4BypytXrohGo5EHHnhg2BuNx5I3cceLC0+cOCEhISHi4OAg8+fPl6KiomFZ+Oeff5aFCxeKs7OzODk5SWhoqHz00UcWazUYDBIVFSVubm7i6OgooaGhyluoRczn776+PklOThZXV1eZNGmSrFq1StavXy9hYWHKcV1dXUo9AOTYsWOSm5sr4eHh4uTkJC4uLhIdHS01NTUjXre5lz+ae7HOF198If7+/mJnZycBAQGyd+9e1TU2NTXJY489Jvb29hIQECBHjhwZdq7Ozk5JSkqSKVOmiEajEX9/f1m5cqVcvXpVRNRZ2lLuJqKJZyViZvEEIiIiIiIiov9iCxYsQHh4OLZv336/SyEiIjP4eDURERERERERERGNK046EtGESklJgbOzs9lPSkrKPasjODh4xDry8/PvWR1ERERENLEMBsOIuc/Z2fl+l2fW/2LNRESW8PFqIppQXV1duHbtmtk2FxcXeHp63pM62tra0N/fb7bNy8sLWq32ntRBRERERBOrt7cXf/zxx4jtOp3uHlYzNv+LNRMRWcJJRyIiIiIiIiIiIhpXfLyaiIiIiIiIiIiIxhUnHYmIiIiIiIiIiGhccdKRiIiIiIiIiIiIxhUnHYmIiIiIiIiIiGhccdKRiIiIiIiIiIiIxhUnHYmIiIiIiIiIiGhccdKRiIiIiIiIiIiIxhUnHYmIiIiIiIiIiGhc/QtYacvuHE59UQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1600x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "t_statistics_custom, p_values_custom = compute_t_stat_p_value_short_qr(X, y, covariates)\n",
    "\n",
    "# Compare the custom function results to those from the statsmodels library\n",
    "feature_names = [f'feature_{i}' for i in range(n_features)]\n",
    "comparison_df = compare_results(X, y, covariates, t_statistics_custom, p_values_custom, feature_names)\n",
    "\n",
    "np.testing.assert_allclose(comparison_df['t_stat_custom'], comparison_df['t_stat_statsmodels'], rtol=1e-1)\n",
    "np.testing.assert_allclose(comparison_df['p_value_custom'], comparison_df['p_value_statsmodels'], rtol=1e-1)\n",
    "\n",
    "fig, axs = plt.subplots(1, 2)\n",
    "comparison_df.plot.scatter(x='t_stat_statsmodels', y='t_stat_custom', ax=axs[0])\n",
    "comparison_df.plot.scatter(x='p_value_statsmodels', y='p_value_custom', ax=axs[1])\n",
    "axs[0].axline((0, 0), slope=1)\n",
    "axs[1].axline((0, 0), slope=1)\n",
    "fig.set_size_inches(16, 4)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}