apoorvalal/event_study_inference_sims.ipynb Secret

## event_study_inference_sims.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-08-10T00:12:19.923053Z",
     "iopub.status.busy": "2024-08-10T00:12:19.922443Z",
     "iopub.status.idle": "2024-08-10T00:12:20.809113Z",
     "shell.execute_reply": "2024-08-10T00:12:20.808749Z",
     "shell.execute_reply.started": "2024-08-10T00:12:19.923024Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "            <div id=\"bQxrnZ\"></div>\n",
       "            <script type=\"text/javascript\" data-lets-plot-script=\"library\">\n",
       "                if(!window.letsPlotCallQueue) {\n",
       "                    window.letsPlotCallQueue = [];\n",
       "                }; \n",
       "                window.letsPlotCall = function(f) {\n",
       "                    window.letsPlotCallQueue.push(f);\n",
       "                };\n",
       "                (function() {\n",
       "                    var script = document.createElement(\"script\");\n",
       "                    script.type = \"text/javascript\";\n",
       "                    script.src = \"https://cdn.jsdelivr.net/gh/JetBrains/lets-plot@v4.3.3/js-package/distr/lets-plot.min.js\";\n",
       "                    script.onload = function() {\n",
       "                        window.letsPlotCall = function(f) {f();};\n",
       "                        window.letsPlotCallQueue.forEach(function(f) {f();});\n",
       "                        window.letsPlotCallQueue = [];\n",
       "                        \n",
       "                    };\n",
       "                    script.onerror = function(event) {\n",
       "                        window.letsPlotCall = function(f) {};    // noop\n",
       "                        window.letsPlotCallQueue = [];\n",
       "                        var div = document.createElement(\"div\");\n",
       "                        div.style.color = 'darkred';\n",
       "                        div.textContent = 'Error loading Lets-Plot JS';\n",
       "                        document.getElementById(\"bQxrnZ\").appendChild(div);\n",
       "                    };\n",
       "                    var e = document.getElementById(\"bQxrnZ\");\n",
       "                    e.appendChild(script);\n",
       "                })()\n",
       "            </script>\n",
       "            "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "\n",
       "            <div id=\"mSK3gL\"></div>\n",
       "            <script type=\"text/javascript\" data-lets-plot-script=\"library\">\n",
       "                if(!window.letsPlotCallQueue) {\n",
       "                    window.letsPlotCallQueue = [];\n",
       "                }; \n",
       "                window.letsPlotCall = function(f) {\n",
       "                    window.letsPlotCallQueue.push(f);\n",
       "                };\n",
       "                (function() {\n",
       "                    var script = document.createElement(\"script\");\n",
       "                    script.type = \"text/javascript\";\n",
       "                    script.src = \"https://cdn.jsdelivr.net/gh/JetBrains/lets-plot@v4.3.3/js-package/distr/lets-plot.min.js\";\n",
       "                    script.onload = function() {\n",
       "                        window.letsPlotCall = function(f) {f();};\n",
       "                        window.letsPlotCallQueue.forEach(function(f) {f();});\n",
       "                        window.letsPlotCallQueue = [];\n",
       "                        \n",
       "                    };\n",
       "                    script.onerror = function(event) {\n",
       "                        window.letsPlotCall = function(f) {};    // noop\n",
       "                        window.letsPlotCallQueue = [];\n",
       "                        var div = document.createElement(\"div\");\n",
       "                        div.style.color = 'darkred';\n",
       "                        div.textContent = 'Error loading Lets-Plot JS';\n",
       "                        document.getElementById(\"mSK3gL\").appendChild(div);\n",
       "                    };\n",
       "                    var e = document.getElementById(\"mSK3gL\");\n",
       "                    e.appendChild(script);\n",
       "                })()\n",
       "            </script>\n",
       "            "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import pyfixest as pf\n",
    "from joblib import Parallel, delayed\n",
    "from tqdm import tqdm\n",
    "\n",
    "np.random.seed(42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-08-10T00:12:22.931694Z",
     "iopub.status.busy": "2024-08-10T00:12:22.931308Z",
     "iopub.status.idle": "2024-08-10T00:12:23.001930Z",
     "shell.execute_reply": "2024-08-10T00:12:23.001507Z",
     "shell.execute_reply.started": "2024-08-10T00:12:22.931682Z"
    }
   },
   "outputs": [],
   "source": [
    "def panel_dgp(\n",
    "    num_units=100,\n",
    "    num_periods=30,\n",
    "    num_treated=50,\n",
    "    treatment_start=15,\n",
    "    hetfx=True,\n",
    "    base_treatment_effect=0.1 * np.log(np.arange(1, 30 - 15 + 1)),\n",
    "    return_dataframe=True,\n",
    "    sigma_unit=1,\n",
    "    sigma_time=0.5,\n",
    "    sigma_epsilon=0.5,  # Standard deviation of epsilon_it\n",
    "    ar_coef=0.8,\n",
    "    het_AR=False,\n",
    "):\n",
    "    unit_intercepts = np.random.normal(0, sigma_unit, num_units)\n",
    "\n",
    "    # Generate day-of-the-week pattern\n",
    "    day_effects = np.array(\n",
    "        [-0.1, 0.1, 0, 0, 0.1, 0.5, 0.5]\n",
    "    )  # Stronger effects on weekends\n",
    "    day_pattern = np.tile(day_effects, num_periods // 7 + 1)[:num_periods]\n",
    "\n",
    "    # Generate autoregressive structure\n",
    "    ar_coef_time = 0.2\n",
    "    ar_noise_time = np.random.normal(0, sigma_time, num_periods)\n",
    "    time_intercepts = np.zeros(num_periods)\n",
    "    time_intercepts[0] = ar_noise_time[0]\n",
    "    for t in range(1, num_periods):\n",
    "        time_intercepts[t] = ar_coef_time * time_intercepts[t - 1] + ar_noise_time[t]\n",
    "    # Combine day-of-the-week pattern and autoregressive structure\n",
    "    time_intercepts = day_pattern + time_intercepts - np.mean(time_intercepts)\n",
    "    # Generate autoregressive noise for each unit\n",
    "    ar_noise = np.random.normal(0, sigma_epsilon, (num_units, num_periods))\n",
    "    if het_AR:\n",
    "        ar_coef = np.random.normal(ar_coef, 0.1, num_units)\n",
    "    noise = np.zeros((num_units, num_periods))\n",
    "    noise[:, 0] = ar_noise[:, 0]\n",
    "    for t in range(1, num_periods):\n",
    "        noise[:, t] = ar_coef * noise[:, t - 1] + ar_noise[:, t]\n",
    "    # N X T matrix of potential outcomes under control\n",
    "    Y0 = unit_intercepts[:, np.newaxis] + time_intercepts[np.newaxis, :] + noise\n",
    "    # Generate the base treatment effect (concave structure)\n",
    "    # Generate heterogeneous multipliers for each unit\n",
    "    if hetfx:\n",
    "        heterogeneous_multipliers = np.random.uniform(0.5, 1.5, num_units)\n",
    "    else:\n",
    "        heterogeneous_multipliers = np.ones(num_units)\n",
    "\n",
    "    # Create a 2D array to store the heterogeneous treatment effects\n",
    "    treatment_effect = np.zeros((num_units, num_periods - treatment_start))\n",
    "    for i in range(num_units):\n",
    "        treatment_effect[i, :] = heterogeneous_multipliers[i] * base_treatment_effect\n",
    "\n",
    "    # random assignment\n",
    "    treated_units = np.random.choice(num_units, num_treated, replace=False)\n",
    "    treatment_status = np.zeros((num_units, num_periods), dtype=bool)\n",
    "    treatment_status[treated_units, treatment_start:] = True\n",
    "\n",
    "    # Apply the heterogeneous treatment effect to the treated units\n",
    "    Y1 = Y0.copy()\n",
    "    for t in range(treatment_start, num_periods):\n",
    "        Y1[:, t][treatment_status[:, t]] += treatment_effect[:, t - treatment_start][\n",
    "            treatment_status[:, t]\n",
    "        ]\n",
    "\n",
    "    result = {\n",
    "        \"Y1\": Y1,\n",
    "        \"Y0\": Y0,\n",
    "        \"W\": treatment_status,\n",
    "        \"unit_intercepts\": unit_intercepts,\n",
    "        \"time_intercepts\": time_intercepts,\n",
    "    }\n",
    "\n",
    "    if return_dataframe:\n",
    "        # Create a DataFrame\n",
    "        unit_ids = np.repeat(np.arange(num_units), num_periods)\n",
    "        time_ids = np.tile(np.arange(num_periods), num_units)\n",
    "        W_it = treatment_status.flatten()\n",
    "        Y_it = np.where(W_it, Y1.flatten(), Y0.flatten())\n",
    "        unit_intercepts_flat = np.repeat(unit_intercepts, num_periods)\n",
    "        time_intercepts_flat = np.tile(time_intercepts, num_units)\n",
    "        df = pd.DataFrame(\n",
    "            {\n",
    "                \"unit_id\": unit_ids,\n",
    "                \"time_id\": time_ids,\n",
    "                \"W_it\": W_it,\n",
    "                \"Y_it\": Y_it,\n",
    "                \"unit_intercept\": unit_intercepts_flat,\n",
    "                \"time_intercept\": time_intercepts_flat,\n",
    "            }\n",
    "        )\n",
    "        # assign units to ever treated if the max of W_it is 1\n",
    "        df[\"ever_treated\"] = df.groupby(\"unit_id\")[\"W_it\"].transform(\"max\")\n",
    "        result[\"dataframe\"] = df\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_periods = 30\n",
    "treat_start_period = 15\n",
    "treat_effect_vector = 0.2 * np.log(\n",
    "    2 * np.arange(1, num_periods - treat_start_period + 1)\n",
    ")\n",
    "treat_effect_vector[8:] = 0  # switch off effects after a week\n",
    "\n",
    "sigma_i, sigma_t = 2, 1\n",
    "\n",
    "df = panel_dgp(\n",
    "    num_units=1000,\n",
    "    num_treated=200,\n",
    "    num_periods=num_periods,\n",
    "    treatment_start=15,\n",
    "    hetfx=False,\n",
    "    base_treatment_effect=treat_effect_vector,\n",
    "    sigma_unit=sigma_i,\n",
    "    sigma_time=sigma_t,\n",
    ")\n",
    "\n",
    "# interact treatment indicator with time dummies omitting the last period before treatment\n",
    "m2 = pf.feols(\n",
    "    \"\"\"\n",
    "              Y_it ~ i(time_id, ever_treated,  ref = 14) |\n",
    "                unit_id + time_id\n",
    "              \"\"\",\n",
    "    df,\n",
    "    vcov={\"CRV1\": \"unit_id\"},\n",
    ")\n",
    "# confidence regions\n",
    "restable = m2.tidy()\n",
    "restable0 = m2.vcov(\"HC1\").tidy()\n",
    "timestamps = restable.index.str.extract(r\"\\[T\\.(\\d+)\\]\", expand=False).astype(int)\n",
    "# create confidence envelopes for all pre-treat and all post-treat coefs\n",
    "kv = [f\"C(time_id, contr.treatment(base=14))[T.{t}]:ever_treated\" for t in range(15)]\n",
    "uniform_confint_pre = m2.confint(keep=kv, joint=True, exact_match=True)\n",
    "kv = [\n",
    "    f\"C(time_id, contr.treatment(base=14))[T.{t}]:ever_treated\" for t in range(15, 30)\n",
    "]\n",
    "uniform_confint_post = m2.confint(keep=kv, joint=True, exact_match=True)\n",
    "uniform_res = pd.concat([uniform_confint_pre, uniform_confint_post])\n",
    "uniform_res.columns = [\"unif_25\", \"unif_975\"]\n",
    "\n",
    "uniform_full = m2.confint(joint=True)\n",
    "uniform_full.columns = [\"unif_full_25\", \"unif_full_975\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7567f08c3200>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABMkAAAH5CAYAAACbEmmeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAACb2klEQVR4nOzdeXgTdf4H8PfkbNOTnim10EILFLnBIipSXFiQs6uw7qrLsS66CisuAgoqBTxQVARXFH8egOu1LmLlEldZi6AICLLcUMpRhB6cvdI0aWZ+f7SdNm0S2tJ0kvT9ep48z+Sb78x8csxk8sn3ECRJkkBERERERERERNSKqZQOgIiIiIiIiIiISGlMkhERERERERERUavHJBkREREREREREbV6TJIREREREREREVGrxyQZERERERERERG1ekySERERERERERFRq8ckGRERERERERERtXoapQNobqIo4vz58wgKCoIgCEqHQ0RERERERERECpIkCcXFxWjbti1UKuftxXwuSXb+/HnExcUpHQYREREREREREXmQs2fP4oYbbnD6uM8lyYKCggBUPvHg4GCFoyEiIiJSltVqxcqVKwEAkydPhlarVTgiIiIiopZVVFSEuLg4OWfkjCBJktRCMbWIoqIihISEoLCwkEkyIiIiavVKS0sRGBgIACgpKUFAQIDCERERERG1rIbmijhwPxERERERERERtXpMkhERERERERERUavHJBkREREREREREbV6PjdwPxERERERERG5j81mg9VqVToMIplWq4Varb7u7TBJRkRERERERETXJEkS8vLycPXqVaVDIaonNDQURqMRgiA0eRtMkhERERERERHRNVUnyKKiomAwGK4rGUHUXCRJgslkQkFBAQAgJiamydtikoyIiIjIh+n1emzYsEFeJiIiagqbzSYnyMLDw5UOh8iOv78/AKCgoABRUVFN7nrJJBkRERGRD9NoNBg5cqTSYRARkZerHoPMYDAoHAmRY9WfTavV2uQkGWe3JCIiIiIiIqIGYRdL8lTN8dlkSzIiIiIiH2a1WvHRRx8BAO677z5otVqFIyIiIiLyTEySEREREfkwi8WCyZMnAwDGjx/PJBkRERGRE0ySEREREREREVGLsYkSdp26jIJiM6KC/JCSEAa1it04SXkck4yIiIiIiIiIWsTmg7m47aX/4o/v/ITpn+7DH9/5Cbe99F9sPpjrlv0JguDyNn/+fLfs15VPPvkEarUaU6dOlctSU1NdxpmamgoAiI+Pd/j4iy++2OLPwxexJRkRERERERERud3mg7l4+MO9kOqU5xWa8fCHe/HW/X0wvFtMs+4zN7cm+favf/0L8+bNw7Fjx+SywMBAeVmSJNhsNmg07k2VvPfee5g9ezbefvttvPrqq/Dz88PatWthsVgAAGfPnkVKSgq+/fZb3HjjjQAAnU4nr79w4UJMmTLFbptBQUFujbm1YEsyIiIiIiIiImo0SZJgslQ06FZstiJ93aF6CTIActn8dYdRbLZec1uS5GgrjhmNRvkWEhICQRDk+0ePHkVQUBC++uor9O3bF3q9Htu3b8ekSZOQlpZmt53HHntMbs0FAKIoYtGiRUhISIC/vz969uyJNWvWXDOeU6dO4ccff8STTz6JTp06Ye3atQCAsLAwOa7IyEgAQHh4uFwWFhYmbyMoKMjueRmNRgQEBDT4NSHn2JKMiIiIiIiIiBqtzGpD13lfN8u2JAB5RWZ0n/+fa9Y9vHAYDLrmS2c8+eSTeOWVV9ChQwe0adOmQessWrQIH374IVasWIGkpCR8//33uP/++xEZGYlBgwY5XW/lypUYOXIkQkJCcP/99+O9997Dvffe21xPha4TW5IRERERERERUau1cOFCDB06FB07drRrseVMeXk5XnjhBbz//vsYNmwYOnTogEmTJuH+++/H22+/7XQ9URSxatUq3H///QCAP/zhD9i+fTtOnTrVqHifeOIJBAYG2t22bdvWqG2QY2xJRkREROTD9Ho9PvvsM3mZiIioufhr1Ti8cFiD6u46dRmTVu6+Zr1Vk29CSoLrRJW/Vt2gfTZUv379GlX/xIkTMJlMGDp0qF25xWJB7969na73zTffoLS0FCNGjAAAREREYOjQoXj//ffx7LPPNnj/s2bNwqRJk+zKYmNjG/4EyCkmyYiIiIh8mEajwfjx45UOgzyAaDLhWJ++AIDOe/dAZTAoHBEReTtBEBrc7XFgUiRiQvyQV2h2OC6ZAMAY4oeBSZFQq4RmjfNa6o7npVKp6o17ZrVa5eWSkhIAwMaNG+slp1z9IfXee+/h8uXL8Pf3l8tEUcT+/fuxYMECqFQN6+wXERGBxMTEBtWlxmGSjIiIiIiIiIjcSq0SkD66Kx7+cC8EwC5RVp0SSx/dtcUTZI5ERkbi4MGDdmX79u2DVqsFAHTt2hV6vR45OTkuxx+r7dKlS/jyyy/x6aefyjNWAoDNZsNtt92G//znPxg+fHjzPQlqEibJiIiIiHxYRUUFvvjiCwDA7373O7dPa09EROTM8G4xeOv+Pliw/jByC81yuTHED+mju2J4txgFo6txxx134OWXX8YHH3yAAQMG4MMPP8TBgwflrpRBQUGYOXMm/v73v0MURdx2220oLCzEDz/8gODgYEycOLHeNv/5z38iPDwcv//97yEI9onAESNG4L333mtwkqy4uBh5eXl2ZQaDAcHBwU18xlSNV0lEREREPqy8vBy///3vAVR2D2GSjIiIlDS8WwyGdjVi16nLKCg2IyrIDykJYR7RgqzasGHD8Mwzz2D27Nkwm83485//jAkTJuDAgQNynWeffRaRkZFYtGgRTp48idDQUPTp0wdz5851uM33338fv/vd7+olyADg7rvvxp/+9CdcvHgRERER14xv3rx5mDdvnl3ZQw89hBUrVjTymVJdglS3o62XKyoqQkhICAoLC5lFJSIiolavtLQUgYGBACqTZHXHXaHWg2OSEdH1MJvNOHXqFBISEuDn56d0OET1uPqMNjRX1LBR4YiIiIiIiIiIiHwYk2RERERERETkkmgy4UiXZBzpkgzRZFI6HCIit2CSjIiIiIiIiIiIWj0myYiIiIiIiIiIqNVjkoyIiIiIiIiIiFo9zgFORERE5MN0Oh1WrlwpLxORd+BspERELY9JMiIiIiIfptVqMWnSJKXDICIiIvJ47G5JRERERERERC3HUgrMD6m8WUqVjoZIxiQZERERkQ+rqKjAxo0bsXHjRlRUVCgdDhERkVdatWoVQkNDW3y/8fHxWLp0qVv3sWXLFiQnJ8Nms7l1P66sWLECo0ePVmz/1ZgkIyIiIvJh5eXlGDVqFEaNGoXy8nKlwyGSiSYTjnRJxpEuyRBNJqXDISJy6Z577sHx48cbtU5qaioee+yx69rv7t278eCDD17XNq5l9uzZePrpp6FWqwG4TggKgoCMjAy7su+++w4jRoxAeHg4DAYDunbtiscffxznzp0DAJjNZkyaNAndu3eHRqNBWlpave3++c9/xt69e7Ft27bmfGqNxiQZEREREREREZEL/v7+iIqKavH9RkZGwuDGiTu2b9+O7Oxs3H333U1a/+2338aQIUNgNBrx+eef4/Dhw1ixYgUKCwvx6quvAgBsNhv8/f3x6KOPYsiQIQ63o9PpcO+99+L1119v8nNpDkySEREREREREZEyis67fRepqamYNm0apk2bhpCQEEREROCZZ56BJElynStXrmDChAlo06YNDAYD7rzzTmRlZcmP121dNX/+fPTq1Qv//Oc/ER8fj5CQEPzhD39AcXExAGDSpEnYunUrli1bBkEQIAgCTp8+jX79+uGVV16Rt5OWlgatVouSkhIAwK+//gpBEHDixAkA9t0tJUnC/Pnz0a5dO+j1erRt2xaPPvqovK3y8nLMnDkTsbGxCAgIQP/+/ZGZmenytfn0008xdOhQ+Pn5Nfp1/fXXX/Hoo4/i0Ucfxfvvv4/U1FTEx8fj9ttvx7vvvot58+YBAAICAvDWW29hypQpMBqNTrc3evRorFu3DmVlZY2OpbkwSUZERERERERELWffxzXLy1OAvR+4fZerV6+GRqPBrl27sGzZMixZsgTvvvuu/PikSZPw888/Y926ddixYwckScKIESNgtVqdbjM7OxsZGRnYsGEDNmzYgK1bt+LFF18EACxbtgwDBgzAlClTkJubi9zcXMTFxWHQoEFy4kqSJGzbtg2hoaHYvn07AGDr1q2IjY1FYmJivf19/vnneO211/D2228jKysLGRkZ6N69u/z4tGnTsGPHDnz66afYv38/xo8fj+HDh9sl++ratm0b+vXr16jXstq///1vWCwWzJ492+HjjR3DrV+/fqioqMDOnTubFE9z0Ci2ZyIiIiIiIiJqXQrPAV/VSqpIIrD+MaDjb4CQWLftNi4uDq+99hoEQUDnzp1x4MABvPbaa5gyZQqysrKwbt06/PDDD7jlllsAAB999BHi4uKQkZGB8ePHO9ymKIpYtWoVgoKCAAB/+tOfsGXLFjz//PMICQmBTqeDwWCwaz2VmpqK9957DzabDQcPHoROp8M999yDzMxMDB8+HJmZmRg0aJDD/eXk5MBoNGLIkCHQarVo164dUlJS5MdWrlyJnJwctG3bFgAwc+ZMbN68GStXrsQLL7zgcJtnzpyR69dWWFiIwMBAl69pVlYWgoODERMT47JeQxkMBoSEhODMmTPNsr2mYEsyIiIiIiIiImoZl7MrE2O1STbg8km37vbmm2+GIAjy/QEDBiArKws2mw1HjhyBRqNB//795cfDw8PRuXNnHDlyxOk24+Pj5QQZAMTExKCgoMBlHAMHDkRxcTF++eUXbN26FYMGDUJqaqrcumzr1q1ITU11uO748eNRVlaGDh06YMqUKfjiiy/kmasPHDgAm82GTp06ITAwUL5t3boV2dnZTuMpKytz2NUyKCgI+/btq3erTZIku9e0Ofj7+8Ok4GQubElGRERERERERC0jrCMgqOwTZYIaCOugXExNpNVq7e4LggBRFJ3UrhQaGoqePXsiMzMTO3bswNChQ3H77bfLs2dmZWU5bUkWFxeHY8eO4dtvv8U333yDRx55BC+//DK2bt2KkpISqNVq7NmzR56lspqrFmERERG4cuVKvXKVSuWwy2dtnTp1QmFhIXJzc5utNdnly5cRGRnZLNtqCrYkIyIiIvJhOp0Ob7zxBt544w3odDqlwyEiotYuJBa4c3HNfUENjF7q1q6WAOqNc/XTTz8hKSkJarUaycnJ9cbCunTpEo4dO4auXbs2eZ86nQ42m61e+aBBg/Ddd9/h+++/R2pqKsLCwpCcnIznn38eMTEx6NSpk9Nt+vv7Y/To0Xj99dflRNuBAwfQu3dv2Gw2FBQUIDEx0e7marD83r174/Dhw016fuPGjYNOp8PixYsdPn716tVGbS87Oxtmsxm9e/duUjzNgS3JiIiIiHyYVqvF1KlTlQ6DiIioRq97gU0zK5en7gQikty+y5ycHMyYMQMPPfQQ9u7di3/84x949dVXAQBJSUkYO3YspkyZgrfffhtBQUF48sknERsbi7FjxzZ5n/Hx8di5cydOnz6NwMBAhIWFQaVSITU1Ff/4xz8QGRmJLl26AKgcq+yNN95wOv4ZUDnDps1mQ//+/WEwGPDhhx/C398f7du3R3h4OO677z5MmDABr776Knr37o0LFy5gy5Yt6NGjB0aOHOlwm8OGDcPq1aub9Pyqx3mbNm0aioqKMGHCBMTHx+PXX3/FBx98gMDAQPk1Pnz4MCwWCy5fvozi4mK562avXr3k7W3btg0dOnRAx44dmxRPc2BLMiIiIiIiavVEkwlHuiTjSJdkiAqOh0PU6gTXHzTeHSZMmICysjKkpKRg6tSpmD59Oh588EH58ZUrV6Jv374YNWoUBgwYAEmSsGnTpnpdKhtj5syZUKvV6Nq1KyIjI5GTkwOgclwyURTtulWmpqbCZrM5HY8MqOyq+c477+DWW29Fjx498O2332L9+vUIDw+Xn8OECRPw+OOPo3PnzkhLS8Pu3bvRrl07p9u87777cOjQIRw7dqxJz/GRRx7Bf/7zH5w7dw6/+93v0KVLF/zlL39BcHAwZs6cKdcbMWIEevfujfXr1yMzMxO9e/eu12Lsk08+wZQpU5oUR3MRJEmSFI2gmRUVFSEkJASFhYUIDg5WOhwiIiIiRdlsNmzbtg1A5UV53XFKqPUQTSYc69MXANB57x6oDAbGw3i8Nh5qeWazGadOnUJCQoLDgd4bxVIKvFCVHJt7HtAFXH+ALqSmpqJXr15YunSpW/fjrWbNmoWioiK8/fbbisVw6NAh3HHHHTh+/DhCQkKatA1Xn9GG5orYkoyIiIjIh5nNZgwePBiDBw+G2WxWOhwiIiLyME899RTat29/zUkH3Ck3NxcffPBBkxNkzYVjkhEREREReYMWbnlBROQ2ugBgfqHSUVCV0NBQzJ07V9EYhgwZouj+q7m9Jdny5csRHx8PPz8/9O/fH7t27XJZ/+rVq5g6dSpiYmKg1+vRqVMnbNq0yd1hEhEREREREZEPyszMZFdLahC3tiT717/+hRkzZmDFihXo378/li5dimHDhuHYsWOIioqqV99isWDo0KGIiorCmjVrEBsbizNnziA0NNSdYRIRERERERERUSvn1iTZkiVLMGXKFEyePBkAsGLFCmzcuBHvv/8+nnzyyXr133//fVy+fBk//vijPINEfHy8y32Ul5ejvLxcvl9UVNR8T4CIiIiIiIiIiFoFt3W3tFgs2LNnj12/UpVKhSFDhmDHjh0O11m3bh0GDBiAqVOnIjo6Gt26dcMLL7wAm83mdD+LFi1CSEiIfIuLi2v250JERERERORrTFYTuq/uju6ru8NkNSkdDhGR4tyWJLt48SJsNhuio6PtyqOjo5GXl+dwnZMnT2LNmjWw2WzYtGkTnnnmGbz66qt47rnnnO5nzpw5KCwslG9nz55t1udBRERERERERES+z6NmtxRFEVFRUfi///s/qNVq9O3bF+fOncPLL7+M9PR0h+vo9Xro9foWjpSIiIjIO2i1WixevFheJt9isprQ/+P+AICd9+6EQWtQOCLyWJwdlYjomtyWJIuIiIBarUZ+fr5deX5+PoxGo8N1YmJioNVqoVar5bLk5GTk5eXBYrFAp9O5K1wiIiIin6TT6TBr1iylwyAiIpIxwU+eym3dLXU6Hfr27YstW7bIZaIoYsuWLRgwYIDDdW699VacOHECoijKZcePH0dMTAwTZERERERERESkiFWrViE0NLTF9xsfH4+lS5e6dR9btmxBcnKyy/HglbZixQqMHj3a7ftxW5IMAGbMmIF33nkHq1evxpEjR/Dwww+jtLRUnu1ywoQJmDNnjlz/4YcfxuXLlzF9+nQcP34cGzduxAsvvICpU6e6M0wiIiIin2Wz2bB7927s3r3boy9+iYiIPNk999yD48ePN2qd1NRUPPbYY9e13927d+PBBx+8rm1cy+zZs/H000/b9eqzWCxYvHgxevbsCYPBgIiICNx6661YuXIlrFYrAGDSpEkQBAGCIECr1SIhIQGzZ8+G2WwGAHTv3h1//etfHe7zn//8J/R6PS5evAiz2YxJkyahe/fu0Gg0SEtLq1f/z3/+M/bu3Ytt27Y1/wtQi1vHJLvnnntw4cIFzJs3D3l5eejVqxc2b94sD+afk5MDlaomTxcXF4evv/4af//739GjRw/ExsZi+vTpeOKJJ9wZJhEREZHPMpvNSElJAQCUlJQgIIDjEBG1CAdjgLGLGZH38vf3h7+/f4vvNzIy0q3b3759O7Kzs3H33XfLZRaLBcOGDcP//vc/PPvss7j11lsRHByMn376Ca+88gp69+6NXr16AQCGDx8uJ8727NmDiRMnQhAEvPTSS3jggQcwf/58vPbaa/Veu5UrV2LMmDGIiIhAaWkp/P398eijj+Lzzz93GKdOp8O9996L119/HQMHDnTb6+HWlmQAMG3aNJw5cwbl5eXYuXMn+vfvLz+WmZmJVatW2dUfMGAAfvrpJ5jNZmRnZ2Pu3Ll22UwiIiIiIiIi8g35pvxrV7pOqampmDZtGqZNm4aQkBBERETgmWeegSRJcp0rV65gwoQJaNOmDQwGA+68805kZWXJj9ftbjl//nz06tUL//znPxEfH4+QkBD84Q9/QHFxMYDKVlZbt27FsmXL5NZWp0+fRr9+/fDKK6/I20lLS4NWq0VJSQkA4Ndff4UgCDhx4gQA++6WkiRh/vz5aNeuHfR6Pdq2bYtHH31U3lZ5eTlmzpyJ2NhYBAQEoH///sjMzHT52nz66acYOnQo/Pz85LKlS5fi+++/x5YtWzB16lT06tULHTp0wL333oudO3ciKSlJrqvX62E0GhEXF4e0tDQMGTIE33zzDQDg/vvvR1lZWb3E16lTp5CZmYkHHngAABAQEIC33noLU6ZMcTqGPQCMHj0a69atQ1lZmcvndD3cniQjIiIiIiIiIqq2LnudvJyWkYa1WWvdvs/Vq1dDo9Fg165dWLZsGZYsWYJ3331XfnzSpEn4+eefsW7dOuzYsQOSJGHEiBFy10JHsrOzkZGRgQ0bNmDDhg3YunUrXnzxRQDAsmXLMGDAAEyZMgW5ubnIzc1FXFwcBg0aJCeuJEnCtm3bEBoaiu3btwMAtm7ditjYWCQmJtbb3+eff47XXnsNb7/9NrKyspCRkYHu3bvLj0+bNg07duzAp59+iv3792P8+PEYPny4XbKvrm3btqFfv352ZR999BGGDBmC3r1716uv1Wqdtko/ePAgfvzxR3lM+YiICIwdOxbvv/++Xb1Vq1bhhhtuwG9/+1uncTnSr18/VFRUYOfOnY1arzGYJCMiIiIiosazlALzQypvllKloyEiL5FXmodFOxfJ90WIWLBjAfJK89y637i4OLz22mvo3Lkz7rvvPvztb3/Da6+9BgDIysrCunXr8O6772LgwIHo2bMnPvroI5w7dw4ZGRlOtymKIlatWoVu3bph4MCB+NOf/iRPXhgSEgKdTgeDwQCj0Qij0Qi1Wo3U1FRs374dNpsN+/fvh06nw3333ScnzjIzMzFo0CCH+8vJyYHRaMSQIUPQrl07pKSkYMqUKfJjK1euxL///W8MHDgQHTt2xMyZM3Hbbbdh5cqVTp/DmTNn0LZtW7uyrKwsdOnSpUGv64YNGxAYGAg/Pz90794dBQUFdrNqP/DAA8jMzMSpU6cAVCYGV69ejYkTJ9oNv9UQBoMBISEhOHPmTKPWawwmyYiIiIiIiIioReQU5UCEaFcmSiLOFp91635vvvlmCIIg3x8wYACysrJgs9lw5MgRaDQau+GhwsPD0blzZxw5csTpNuPj4xEUFCTfj4mJQUFBgcs4Bg4ciOLiYvzyyy/YunUrBg0ahNTUVDlJtnXrVqSmpjpcd/z48SgrK0OHDh0wZcoUfPHFF6ioqAAAHDhwADabDZ06dUJgYKB827p1K7Kzs53GU1ZWZtfVEoBdN9RrGTx4MPbt24edO3di4sSJmDx5st34ZkOHDsUNN9wgJ+q2bNmCnJwceULHxvL394fJZGrSug3BJBkRERERETULk9WE7qu7o/vq7jBZ3fcjhoi8V7vgdlDVSUWoBBXiguIUiqjptFqt3X1BECCKopPalUJDQ9GzZ09kZmbKCbHbb78dv/zyC44fP46srCynLcni4uJw7NgxvPnmm/D398cjjzyC22+/HVarFSUlJVCr1dizZw/27dsn344cOYJly5Y5jSciIgJXrlyxK+vUqROOHj3aoNcgICAAiYmJ6NmzJ95//33s3LkT7733nvy4SqXCpEmTsHr1aoiiiJUrV2Lw4MHo0KFDg7Zf1+XLl906mQGTZERERERERETUIowBRszpP0e+rxJUSB+QDmOA8wHbm0Pdcax++uknJCUlQa1WIzk5ud5YV5cuXcKxY8fQtWvXJu9Tp9PBZrPVKx80aBC+++47fP/990hNTUVYWBiSk5Px/PPPIyYmBp06dXK6TX9/f4wePRqvv/46MjMzsWPHDhw4cAC9e/eGzWZDQUEBEhMT7W6uBsPv3bs3Dh8+bFd277334ttvv8Uvv/xSr77VakVpqeMu9iqVCnPnzsXTTz9tN7j+5MmTcfbsWaxduxZffPGFPGB/Y2VnZ8NsNjscK625MElGRERE5MO0Wi3S09ORnp5e7x9vci/RZMKRLsk40iUZohu7hhA1J5soYUf2JXy57xx2ZF+CTWx4tyuihhrTcYy8nDE2A3cl3eX2febk5GDGjBk4duwYPvnkE/zjH//A9OnTAQBJSUkYO3YspkyZgu3bt+N///sf7r//fsTGxmLs2LFN3md8fDx27tyJ06dP4+LFi3Irs9TUVHz99dfQaDTy2F+pqan46KOPnLYiAyoHvH/vvfdw8OBBnDx5Eh9++CH8/f3Rvn17dOrUCffddx8mTJiAtWvX4tSpU9i1axcWLVqEjRs3Ot3msGHD5EkDqj322GO49dZb8Zvf/AbLly/H//73P5w8eRKfffYZbr75ZpcTAYwfPx5qtRrLly+XyxISEnDHHXfgwQcfhF6vx1131X+/Dx8+jH379uHy5csoLCyUW8LVtm3bNnTo0AEdO3Z0uv/rpXHblomIiIhIcTqdDvPnz1c6DCLyApsP5mLB+sPILTTLZTEhfkgf3RW/7RCiYGTky6IN0S2ynwkTJqCsrAwpKSlQq9WYPn06HnzwQfnxlStXYvr06Rg1ahQsFgtuv/12bNq06br+YJo5cyYmTpyIrl27oqysDKdOnUJ8fDwGDhwIURTtEmKpqalYtmyZ0/HIgMqumi+++CJmzJgBm82G7t27Y/369QgPD5efw3PPPYfHH38c586dQ0REBG6++WaMGjXK6Tbvu+8+zJ49G8eOHUPnzp0BAHq9Ht988408k+bMmTNhMBiQnJyMRx99FN26dXO6PY1Gg2nTpmHx4sV4+OGH5ZkwH3jgAWzZsgWPPPJIvTHQAGDEiBF2A/JXtxarPT7aJ598Ik9U4C6C1JgR2bxAUVERQkJCUFhYiODgYKXDISIiIqJWSjSZcKxPXwBA5717oDIYrm+DllLghaoZyOaeB3QBMFlN6P9x5UDTO+/dCYPW+T58Pp7r5POvzzXieb7Xl5j+yRHU/XFYPcz5inHJaHf/yOaLh7yO2WzGqVOnkJCQ4DDJ0RiNORaaQ2pqKnr16oWlS5e6dT/eatasWSgqKsLbb7+tdChOHTp0CHfccQeOHz+OkBDHSXtXn9GG5orY3ZKIiIjIh4miiEOHDuHQoUPXHEyYiFqvF746Wi9BBkAue2FTwwbxJiLv89RTT6F9+/YefZ2Qm5uLDz74wGmCrLmwuyURERGRDysrK5O7RZSUlMjdHohag5ZureLN8grLAegcPiYByCsyO3yMqCkMWgMOTDygdBhUJTQ0FHPnzlU6DJeGDBnSIvthkoyIiIiIiLyfg+6EREQAkJmZqXQI5CWYJCMiIiIiImpFDp0vwn+zz+PbY2cB5hKJiGRMkhEREREREfmw0vIK7Diaj+rOSuNW/Igy+AGCBUFdKssC9WqUOOlRKQAwBl/fQO1ERN6ASTIiIiIiIiIfc/piKf57tADfHSvAzpOXobaZcKQqz2XQqXFbYjRu6xSMV45Vlr3wu26Y/skRALAbwL96dsu5I7oAH7dY+EREimCSjIiIiIiIyMtZKkTsPn25MjF2tAAnL5baPd6pjQEoq1z+8ck7oDcEw2Q1yUmyIV2NeOt+fyxYfxi5hTVNyqKD9Zg/5kYM7RCCYy31ZIiIFMIkGRERERERkRcqKDYj8+gF/PdoAbafuIiS8gr5MY1KwE3xYbijSxQGd4lCxxAAiyof02vUDrc3vFsMhnY1Ytepy3jwg59RXF6B5ff1Rd/2bSCaTC3wjIiIlMUkGREREZEP02q1mDlzprxMRJ7FJkrYdeoyCorNiAryQ0pCGNQqwWFdUZRw4FwhtlS1FjtwrtDu8YhAHVI7R+GOLlG4LSkCwX61jnlLKRpCrRIwoGM4ut8Qgh+zL+HkhRL0bd+myc+PyBHRZMKxPn0BAJ337oHKYFA4IqJKTJIRERER+TCdToeXX35Z6TCIyIHNB3PrdW+MCfFD+uiu+G2HELns60N52HK6GJnHLuBiSbndNnrcEILBVYmx7rEhUDlJsDVWUlQgfsy+hBMFJc2yPSJvkZmZicGDB+PKlSsIDQ1t9u1funQJycnJ2LVrF+Lj45t9+57qD3/4A2666SY8/vjjSofiEpNkRERERERELezbw3mY/skRu0HyASCv0Iy/frgX93QLx6Sqssf+tQ/lGj0AIFCvwcCkCAzuEoXUzpGICnLPrJOJ0UEAgCwmycgHTJo0CVevXkVGRsY1695yyy3Izc1FSEjINes2ZfvPP/88xo4d26IJsq1bt+L+++/H2bNn3b4vZ0nGp59+Grfffjv+8pe/NOq1bWlMkhERERH5MFEUkZOTAwBo164dVCqVwhEREQC88NXRegkyoGZmyYx95+UkWUJEAG7rFoc7ukShX3wYdBr3H8dJUYEAgKyCYrfvi8iT6HQ6GI1Gt2zbZDLhvffew9dff31d27HZbBAEocHf6V9++SVGjx59Xfu8Xt26dUPHjh3x4YcfYurUqYrG4gqvkoiIiIh8WFlZGRISEpCQkICysjKlwyGiKnmF5deuVGXTowPx9KiuuCUxokUSZEBNkuzXK2UwWSquUZuo6az5+S26v/Lycjz66KOIioqCn58fbrvtNuzevVt+PDMzE4Ig4OrVqwCAVatWITQ0FF9//TWSk5MRGBiI4cOHIzc3FwAwf/58rF69Gl9++SUEQYAgCMjMzHS4702bNkGv1+Pmm2+ut7+NGzeiR48e8PPzw80334yDBw/KdapjWLduHbp27Qq9Xo+cnByUl5dj5syZiI2NRUBAAPr37+9w3+vWrcOYMWMAAKmpqZg2bRqmTZuGkJAQRERE4JlnnoEk1aTtr1y5ggkTJqBNmzYwGAy48847kZWVJT9+5swZjB49Gm3atEFAQABuvPFGbNq0CadPn8bgwYMBAG3atIEgCJg0aZK83ujRo/Hpp5827I1SCJNkREREROQTRJMJR7ok40iXZM7ERx7JJopKh9Bg4YF6hAXoIEnAyQsNG/SfqKGu1uqWeHLkKFxds6bF9j179mx8/vnnWL16Nfbu3YvExEQMGzYMly9fdrqOyWTCK6+8gn/+85/4/vvvkZOTI0+KM3PmTPz+97+XE2e5ubm45ZZbHG5n27Zt6Nu3r8PHZs2ahVdffRW7d+9GZGQkRo8eDavVahfDSy+9hHfffReHDh1CVFQUpk2bhh07duDTTz/F/v37MX78eAwfPtwuoXXo0CEUFBTgjjvukMtWr14NjUaDXbt2YdmyZViyZAneffdd+fFJkybh559/xrp167Bjxw5IkoQRI0bI8UydOhXl5eX4/vvvceDAAbz00ksIDAxEXFwcPv/8cwDAsWPHkJubi2XLlsnbTUlJwa5du1Be3vA/CVoau1sSERERERG50c6Tl7H+6El8degMcIPS0TRcYlQgdp26jKyCYnRtE6Z0OOQjrHl5yH/u+ZoCUUTuvHQE3HYbtG7q5littLQUb731FlatWoU777wTAPDOO+/gm2++wXvvvYdZs2Y5jtlqxYoVK9CxY0cAwLRp07Bw4UIAQGBgIPz9/VFeXn7NbppnzpxB27ZtHT6Wnp6OoUOHAqhMYt1www344osv8Pvf/16O4c0330TPnj0BADk5OVi5ciVycnLkbc6cORObN2/GypUr8cILLwCo7Go5bNgw6HQ6eV9xcXF47bXXIAgCOnfujAMHDuC1117DlClTkJWVhXXr1uGHH36Qk30fffQR4uLikJGRgfHjxyMnJwd33303unfvDgDo0KGDvO2wsMpzRVRUVL2JD9q2bQuLxYK8vDy0b9/e5WulFLYkIyIiIiIiakaiKGHPmSvy/UmrduHDn3JwqaSmVUibAC2czUMpADAGu2dA/saQxyXL5+D91Hwsp88AdVtViiIsZ3Lcvu/s7GxYrVbceuutcplWq0VKSgqOHDnidD2DwSAnyAAgJiYGBQUFjd5/WVkZ/PwcH9sDBgyQl8PCwtC5c2e7mHQ6HXr06CHfP3DgAGw2Gzp16oTAwED5tnXrVmRnZ8v1vvzyS7mrZbWbb74ZglBzBhowYACysrJgs9lw5MgRaDQa9O/fX348PDzcLp5HH30Uzz33HG699Vakp6dj//79DXr+/v7+ACpbxXkqtiQjIiIiIiK6TpIk4ZezV7Hhf7nYdCAXhUVXcaTqt3CInxaju92AoTe2wYxdlWULRnfF9E+OQADsBvCv/tk6d0QX4OMWfAIO1AzezyQZNR9dfHtApbJPlKlU0LVvp1xQ16DVau3uC4JgN4ZXQ0VERODKlSvXruiAv7+/XWKrpKQEarUae/bsgVqttqsbGFh57Obm5uKXX37ByJEjm7RPZ/7yl79g2LBh2LhxI/7zn/9g0aJFePXVV/G3v/3N5XrVXVojIyObNZ7mxJZkRERERERETSBJEvb/ehUvbDqC2176Dne9+SPe/+EU8orMCNLXtEf4fvZgLB7XE7cmRshlQ7oa8db9fWAMsW9VYgzxw1v398HQru7tdtYQSdFBAIATTJJRM9IajYh++qmaApUKMQsXuL2rJQB07NgROp0OP/zwg1xmtVqxe/dudO3atcnb1el0sNls16zXu3dvHD582OFjP/30k7x85coVHD9+HMnJyS63ZbPZUFBQgMTERLtbdbfP9evX45ZbbpG7QFbbuXNnvX0nJSVBrVYjOTkZFRUVdnUuXbqEY8eO2b1GcXFx+Otf/4q1a9fi8ccfxzvvvCO/FgAcvh4HDx7EDTfcgIiIiHqPeQq2JCMiIiIiImogSZJwOLcIG/bnYuP+XORcruk2FKBTY0jXaIzq0Ra3x/sDiyvLnc1IObxbDIZ2NWLXqcsoKDYjKsgPKQlhUKsEj5h8orol2ZlLpTBbr50AIGqo0LQ05C98FgDQYeMG6BMSWmS/AQEBePjhhzFr1iyEhYWhXbt2WLx4MUwmEx544IEmbzc+Ph5ff/01jh07hvDwcISEhNRrfQYAw4YNw5w5c3DlyhW0adPG7rGFCxciPDwc0dHReOqppxAREYG0tDSn++zUqRPuu+8+TJgwAa+++ip69+6NCxcuYMuWLejRowdGjhxpN6tlbTk5OZgxYwYeeugh7N27F//4xz/w6quvAgCSkpIwduxYTJkyBW+//TaCgoLw5JNPIjY2FmPHjgUAPPbYY7jzzjvRqVMnXLlyBd99952c0Gvfvj0EQcCGDRswYsQI+Pv7yy3btm3bht/+9rdNeo1bCpNkRERERD5Mo9HgkUcekZeJqPEkScKx/GJs3J+LDftzcepizWyP/lo17kiOwugeMUjtHAU/bVW3J0vDZoRUqwQM6BjujrCvW2SQHsF+GhSZK3DmkvJJO/JN2uhot+9DFEX5O/DFF1+EKIr405/+hOLiYvTr1w9ff/11vaRVY0yZMgWZmZno168fSkpK8N133yE1NbVeve7du6NPnz747LPP8NBDD9k99uKLL2L69OnIyspCr169sH79ervB9h1ZuXIlnnvuOTz++OM4d+4cIiIicPPNN2PUqFEoLS3Fli1bsHTp0nrrTZgwAWVlZUhJSYFarcb06dPx4IMP2m13+vTpGDVqFCwWC26//XZs2rRJTvzZbDZMnToVv/76K4KDgzF8+HC89tprAIDY2FgsWLAATz75JCZPnowJEyZg1apVMJvNyMjIwObNmxv56rYsXikRERER+TC9Xo/ly5crHQaRR7GJUr3WW44G0T9RUIwNVYmx2l0O9RoVBneOwqieMbijSxQMOt/8WSUIApKig7DnzBWcuFCCRKUDImqi6i6JAODn54fXX38dr7/+usO6qampduONTZo0CZMmTbKrk5aWZlcnMjIS//nPfxoUy7x58zBr1ixMmTIFKlVNK9PbbrsNBw8edLiOoxiAyrHSFixYgAULFtR7bO3atUhISJCfd931li5dirfeesvh/tq0aYMPPvjA6XP4xz/+4fQxAHjmmWfwzDPP2JWtXLkSKSkpuPnmm12uqzTfPJsTERERERE5sPlgLhasP4zcQrNcFhPih/lDE1A9bPhbmdlYf/wKjuYVy3V0ahUGdY7EqB4x+E1yNAL1reOnVFJUIPacuYKTBUySkfe5cuUKfvjhB2RmZuKvf/2r0uEAAEaOHImsrCycO3cOcXFxbttPYGAgXnrpJbdtv7G0Wu01k2ueoHWc2YmIiIhaKUmScPHiRQCVs2rVnhmLqLXZfDAXD3+4F3XnpMstNOPRT/cho+r+6//NQrlGD61awMCkSIzsHoOhN0Yj2K/+GEO+LrHWDJeePZIQeROVwYDko0fcvp8///nP2L17Nx5//HF5PC1P8Nhjj7l9H5429tdf/vIXpUNoECbJiIiIiHyYyWRCVFQUgMrp4gMCAhSOyItYSoEX2lYuzz0P6PjaeTNRlLBg/eF6CTJHbkuMxLC+8Rh2oxEhhtaXGKutU9UMl9kXOMMleZ8vvvhC6RCuqW73TnfLzMxssX15IybJiIiIiIjI5/185opdF0tX3pnQFyqDwc0ReYek6KoZLi9z4H4i8n2O5yImIiIiIiLyIat2nFY6BK9kDPZDoF4Dm9hyLV3Is7VkqyeixmiOzyaTZERERERE5HMsFSK+Ppgn3888ekHBaLyXIAjyuGTUumm1lV2PTSa2KiTPVP3ZrP6sNgW7WxIRERERkU955T/H8Mkvl3DJVIKgLpVltyaG4/A5M66arA7HJeOUFs4lRQXiyKl8pcMghanVaoSGhqKgoAAAYDAYOBkMeQRJkmAymVBQUIDQ0FCo1eomb4tJMiIiIiIi8mqWChH/PZiH4VX339t+CmXwQ2SwDtWjkL0zoR++P1aIhz/cCwGwS5TxZ75r1eOSERmNRgCQE2VEniQ0NFT+jDYVk2REREREROSVTl8sxae7z2LNnrMoLSnCcL/K8tsSIzBuQGcMSAzErZ/W1B/eLQZv3d8HC9YfthvE3xjih/lDk4ENLfwEvERSVJDSIZCHEAQBMTExiIqKgtVqVTocIplWq72uFmTVmCQjIiIi8mEajQYTJ06Ul4m8naVCxDeH8/HJrhxsP3FRLo8L1AMVlcvvTOgH6AJgstYfO2l4txgM7WrErlOXUVBsRlSQH1ISwiCYy3CspZ6El+GYZFSXWq1uloQEkafhlRIRERGRD9Pr9Vi1apXSYRBdtzOXSvHJrspWYxdLLAAAQQBuT4rEH1Pa4TcdA4CXGrYttUrAgI7hdmVicwfsQ2JD/eGvZUKEiHwfk2RERERERA1ksprQ/+P+AICd9+6EQWtQOCLf5qzVWGSQHvf0i8M9N8UhLqzqPbCUKhSl71OpBHSIZGsyIvJ9TJIRERER+bDqGZ8AzkTmU4rOAxFJdkX5pnwkhCQ0aHVrfj70CQ2rq4RrthpLjoJWrVI4ytYlMTJA6RCIiNyOSTIiIiIiH2YymRAYWNkCpKSkBAEB/KHrtfZ9XLO8PAUYvQzrDDq5KC0jDem3pOOupLscrn41I0NePjlyFGIWLkDouHHuitaOTZTqjwFWp46lQsS3R/Lx8c4GtBqjFteB45IRUSvAJBkRERERkacrPAd8NbvmviQib9MMLIprKxeJELFgxwLc0vYWGAOMdqtb8/KQ/9zzNQWiiNx56Qi47TZojfZ1m9vmg7n1ZpOMCfHD/KEJaFd1f8k3x/HpgQtsNebBEtndkohaASbJiIiIiIg83eVsQLIfWj5Ho4IIya5MlEScLT5bL0lmOX0GEOsMTS+KsJzJcWuSbPPBXDz84d46UQJ5hWY8+uk+ZFTdf2fbSZRr9Gw15sE6RgaiOs1pEyUwbUlEvohJMiIiIiIiTxfWERBUdomydhUiVBDsEmUqQYW4oLh6q+vi2wMqlX2iTKWCrn27enWbiyhKWLD+cL0EGYB6ZQMTIzHutiS2GvNgN7Txx4mq5XNXyxAfyK7bROR7+A1EREREROTpQmKBOxfX3BfUMI5Ygjn958pFKkGF9AHp9VqRAYDWaET000/VFKhUiFm4wK2tyH4+c8Wui6Ur/zehL4Z3MzJB5sHUqppR5E4UlCgYCRGR+/BbiIiIiIg8g6UUmB9SebOUKh2N5+l1b83y1J1AnwkY03GMXJQxNsPpoP0AEJqWJi932LjB7YP2Xygud+v2STknLjBJRkS+iUkyIiIiIiJvE9y2XlG0IbrBq2ujG163qd787sS1K5Eyis7XK8o35Td49bwTOc0ZDRGRx2CSjIiIiMiHqdVqjBs3DuPGjYNarVY6HPJhoijhm0N58v1Tl0wu6wsuH6Vmt+/jmuXlKcDeD7Aue51clJaRhrVZa52ufjUjQ14e/8ZMXF2zxh1REhEpigP3ExEREfkwPz8//Pvf/1Y6DPJxGw/kYtn353C84DKCulSWTbk9AZ0iwzF7zX4A9oP1M0HWwgrPAV/NrrkvicjbNAOL4mpaJIoQsWDHAtzS9pZ649pZ8/KQ/9zz8n2VJCF3XjoCbrvNrePaERG1NLYkIyIiIiKiRquw1cyUOfPf/8Px/BIE+dW0Vvz7kE4Y3y8Ob93fB8YQP7t1jSF+WPaHXi0VKl3OtpsZFQByNCq7mVEBQJREnC0+W291y+kz9jOjAoAownKG3S6JyLewJRkRERERETWYpULEF7/8ive/O4Svq8pC/LV45LZOGJ8SjaF1euwN7xaDoV2N2HXqMgqKzYgK8kNKQhgEcxmOtXj0Cik6D0Qk2RXlm/KREJLQoNWt+fnQJzSsrkNhHQFBZZcoa1chQgXBLlGmElSIC4qrt7ouvj2gUtklyiSVCrr27ZoeExGRB2JLMiIiIiIfVlpaCkEQIAgCSks5YyQ1XXmFDR/+dAaDX8nEE58fQM7lmjHHtkzphL/9JgnBfjX/wdceCF6tEjCgYzjG9orFgI7hUKvsO1xa8xs+aLzXaMYxwE6OHHV9Y4CFxAJ3Lq65L6hhHLEEc/rPlYtUggrpA9LrdbUEAK3RiOinn5Lv2wQBWfdPY1dLIvI5TJIREREREZFTZqsNq344hUGLM/F0xkGcu1qGyCA93u2ZJdcJeGeAskkgT+NsDLCdL8hF1WOA5ZXm1Vu97hhgEEXkzkuHNa9+3QbrdW/N8tSdQJ8JGNNxjFyUMTYDdyXd5XT10LQ0efnBO2Zhe9ItTY+FiMhDsbslERERERHVY7JU4OOdOXj7+5O4UFwOADAG++GvgzrgD13U8HtjfE3l6xwIvjoJ5DMDwTdyDLC6r4+rMcCa5fUJbluvKNoQ3eDVL/mHIKug5PrjICLyMEySERERERGRrKS8Av/ccQbvbDuJy6UWAEBsqD8eGdwR4/reAL1GDZz63rOTQEpzwxhg8LAxwE4UlECSJAgC5yolIt/B7pZERERERITCMite35KF2176L17afBSXSy1oF2bA4rt7IHNWKu7r374yQQbUJIFqqU4C1XbNJJBdZc9KAl2XZh4DDCoVYhYu8JgEololoKS8AnlFZqVDISJqVkySERERERG1YldNFiz5zzHc9tJ/seSb47hqsqJDZACW/L4n/vv4IPz+pjho1XV+Nvh4EqhZNOMYYB02bkDouHHuiLJJ2ocZAABZ+exySUS+hd0tiYiIiIhaoUsl5Xh3+yl88ONplFpsAIBO0YGYdkcSRnaPqTcDZT297gU2zaxcnroTiEjCGKsJz++sHGssY2wGEkISnK4empaG/IXPAqhMAukTnNf1etc5Bpg2uuF1W0LHyEAcuXoVWQUluL1TpNLhEBE1GybJiIiIiHyYWq3GiBEj5GVqHWyihF2nLqOg2IyoID+kJITZdYRcvPkoVv9SgDJrZXIsOSYYj96RiGE3GqG6VnLMER9LApFriVGBQNZVnCgoVjoUIqJmxSQZERERkQ/z8/PDxo0blQ7jupisJvT/uD8AYOe9O2HQGhSOyLNtPpiLBesPI7ewZryomBA/PHZbHHpU3V/542mUa/ToeUMI/nZHEn6THOVbA7AXnQcikuyK8k35Llu21WbNz/ftlm3XqWNUIAB2tyQi38MxyYiIiIiIfMTmg7l4+MO9dgkyAMgtNGPelwfl+73iQrFq8k3ImHorhnSN9o0E2b6Pa5aXpwB7P8C67HVyUVpGGtZmrXW6+tWMDHn55MhRuLpmjTui9AmJkZVJsuP5xZAk6Rq1iYi8B5NkREREREQ+QBQlLFh/GA1JWXz8l/5I7exDrccKzwFfza65L4nI2zQDi3a+IBeJELFgxwLklebVW92al4f8556vKRBF5M5LhzWvfl0C2ocboBKAInMFLhSXKx0OEVGzYZKMiIiIyIeVlpYiICAAAQEBKC0tVToccqOfz1yp14LMGZ9JjlW7nA1Iol1RjkYFsU7KUJREnC0+W291y+kzgGi/PkQRljM5zR6qL/DTqtE+PAAAkFXALpdE5DuYJCMiIiLycSaTCSaTSekwyM1adYuesI6AYP/Tpl2FCBXsk4EqQYW4oLh6q+vi2wOqOj+NVCro2rdr9lB9RaI8LhkH7yci38EkGRERERGRD/hi3zmlQ1BOSCxw5+Ka+4IaxhFLMKf/XLlIJaiQPiAdxgBjvdW1RiOin36qpkClQszCBdAa69elSknVSTK2JCMiH9IiSbLly5cjPj4efn5+6N+/P3bt2tWg9T799FMIgoC0tDT3BkhERERE5IVqD5r+44lLLuv6WAfL+nrdW7M8dSfQZwLGdBwjF2WMzcBdSXc5XT201m+ODhs3IHTcOHdE6TOSopkkIyLf4/Yk2b/+9S/MmDED6enp2Lt3L3r27Ilhw4ahoKDA5XqnT5/GzJkzMXDgQHeHSERERETkdU5fLMVfPvhZvt8pOhCzh3WGgPoJMZ9PkNUV3LZeUbQhusGra6MbXre1SooKAgCcYJKMiHyI25NkS5YswZQpUzB58mR07doVK1asgMFgwPvvv+90HZvNhvvuuw8LFixAhw4d3B0iEREREZFXefv7kxi29HvsyL4sl615eAAeGZyIt+7vA2OIn119Y4gflv2hVwtHSb6sY2QgBAG4XGrBpZJWPB4eEfkUjTs3brFYsGfPHsyZM0cuU6lUGDJkCHbs2OF0vYULFyIqKgoPPPAAtm3b5nIf5eXlKC+vOSkXFRVdf+BERERERB5s6bfHUQ4/3JIUjgNVZZqqgeeHd4vB0K5G7Dp1GQXFZkQF+SElIQyCuQzHlAuZfIy/To0b2vjj7OUyZBWUIDxQr3RIRETXza1JsosXL8JmsyG6TnPl6OhoHD161OE627dvx3vvvYd9+/Y1aB+LFi3CggULrjdUIiIiIp+kUqkwaNAgeZm8U7HZiqWbDuOZqvttDDq8MLonfntjG9z8Sf36apWAAR3D7cpE94dJrUxSVJCcJLu5Q/i1VyAi8nAedaVUXFyMP/3pT3jnnXcQERHRoHXmzJmDwsJC+Xb27Fk3R0lERETkPfz9/ZGZmYnMzEz4+/srHQ41wdeH8jB0yff4eFeOXLbhb7fhd71vgCC0utHGyINUz3B5Ir9Y4UiIiJqHW1uSRUREQK1WIz8/3648Pz8fRgfTKWdnZ+P06dMYPXq0XCaKlf95aTQaHDt2DB07drRbR6/XQ69n014iIiIi8i15hWakrzuIrw9VXkt3CTMApsrHwgJ0CkZGVCkxijNcEpFvcWtLMp1Oh759+2LLli1ymSiK2LJlCwYMGFCvfpcuXXDgwAHs27dPvo0ZMwaDBw/Gvn37EBcX585wiYiIiIgUZxMlfLDjNIYs2YqvD+VDoxIwbXAiMqbeqnRoRHaSoitnuGSSjIh8hVtbkgHAjBkzMHHiRPTr1w8pKSlYunQpSktLMXnyZADAhAkTEBsbi0WLFsHPzw/dunWzWz80NBQA6pUTERER0bWVlpYiPj4eAHD69GkEBAQoGxC5dDSvCE9+fgD7zl4FAPRpF4pFd/VAZ2MQYClVNjiiOqpbkl0oLsdVkwWhBrZwJCLv5vYk2T333IMLFy5g3rx5yMvLQ69evbB582Z5MP+cnBwOIktERETkRhcvXlQ6BLoGs9WG17dk4f++P4kKUUKQXoPZd3bBfSntoFJx3DHyTIF6DWJD/XHuahlOFJSgX3yY0iEREV0XtyfJAGDatGmYNm2aw8cyMzNdrrtq1armD4iIiIiIyENsz7qIpzIO4MylygHHht9oxPwxN8IY4qdwZETXlhgViHNXK2e4ZJKMiLxdiyTJiIiIiIjI3qWScjy/8QjW/nIOAGAM9sPCsTfitzfWn+CKyFMlRQVi6/ELyMrnuGRE5P2YJCMiIiIiakGSJGHt3nN4buNhXDFZIQjAxAHxePy3nRDkp1U6PKJGSYqunuGyWOFIiIiuH5NkREREREQt5PTFUjyVcQA/nLgEAOhiDMKiu7qjd7s2CkdG1DSJUVUzXLIlGRH5ACbJiIiIiIjczGoT8e53J/D6liyUV4jQa1T4+9BOeOC2BGjVnMSKvFf1DJd5RWYUma0IZmtIIvJiTJIRERER+TCVSoV+/frJy6SMcW/twIHLFgDAwKQIPJfWDe3DAxSOiuj6hfhrER2sR35ROU4UlKAPW0USkRdjkoyIiIjIh/n7+2P37t1Kh+HzbKKEXacuo6DYjKggP6QkhKGk3Co/frygGGEhQXhmVDLSesVCEAQFoyVqXklRQZVJsnwmyYjIuzFJRkRERER0LUXngYgku6J8Uz4SQhKw+WAuFqw/jNxCs/xYqEELndWM1VX37+voj79NGIQ2AboWDJqoZSRGBWL7iYscvJ+IvB6TZEREREREjuz7uGZ5eQowehnWGWqSXGkZabi7/XSs3BwNqc6qV01WjDz5k3x//D9mQohcAIwb5+agiVpezQyXHLyfiLwbB6YgIiIi8mEmkwnx8fGIj4+HyWRSOhzvUXgO+Gp2zX1JRN6mGVi08wW5SISIf59eCmgK660eUXYVj+zPqCkQReTOS4c1L899MRMpJIkzXBKRj2CSjIiIiMiHSZKEM2fO4MyZM5Ckuu2dyKnL2YAk2hXlaFQQ67YZEySodBfrrd625CJUdeuKIixncpo7UiLFJVXNcHnuahlKyysUjoaIqOmYJCMiIiIiqiusIyDYXyq3qxChgv2A+5IkQLRE1Fv9fGAExDp1oVJB175ds4dKpLQ2ATpEBFZ2Rc6+wNZkROS9mCQjIiIiIqorJBa4c3HNfUEN44glmNN/bk0RVCjPvQtSRUi91S/6h+LNHmk1BSoVYhYugNZodGPQRMpJrGpNxi6XROTNmCQjIiIiInKk1701y1N3An0mYEzHMXJRyMU5sBbe5HT1Le36ycsdNm5AKAftJx8mj0vGwfuJyIsxSUZEREREdC3BbesVnb2oR1jVbJd1OlbWu6+NjnZPXEQeolPVDJcnCooVjoSIqOmYJCMiIiIiaqBtWRfk5S7GIGx+bCBW3N8HxhA/u3rGED8s+0OvFo6OSDmJbElGRD5Ao3QAREREROQ+giCga9eu8jI13Zo9v+LJtb/Av1Pl/Q8eSEFUoB+Gd4vB0K5G7Dp1GQXFZkQF+SElIQyCuQzHlA2ZqMUkVbUky7lsgtlqg59WrXBERESNxyQZERERkQ8zGAw4dOiQW7Ytmkw41qcvAKDz3j1QGQxu2Y8nePv7k1j0bQ4gSHJZoL7mUlqtEjCgY7jdOmKLRUekvPAAHdoYtLhisiL7QglubFt/QgsiIk/H7pZERERERNew9NvjAIC/DExQOBJ7Bq0BByYewIGJB2DQ+m6SkjyfIAg1g/dzhksi8lJMkhEREREROVBmscnLggAsGHMjZgztpGBERJ4tsarLZRYH7yciL8UkGREREZEPM5lMuPHGG3HjjTfCZDIpHY7XuFJqwZ9X75bvL/19L0y8JV65gIi8QFJUVZKMLcmIyEtxTDIiIiIiHyZJEg4fPiwv07WdvWzCxJW7kHvhKlA1aeVvbzQqGhORN6jubnmCM1wSkZdiSzIiIiIioioHzxXirrd+xMkLpYgJ9lM6HCKvUj3D5elLpSivsF2jNhGR52GSjIiIiIg8T9F5h8X5pvwGrW7Nb1i92rZlXcA9b+/AheJydDEG4ZMHb270Nohas6ggPYL8NBAl4NTFUqXDISJqNCbJiIiIiMgz7Pu4Znl5CrD3AwDAuux1cnFaRhrWZq11uPrVjAx5+eTIUbi6Zk2Dd71276+YvHI3Si023NIxHJ/9dQCiPa0lmS4AmF9YedMFKB0NUT2VM1xyXDIi8l5MkhERERGR8grPAV/NrrkvicD6x5CXtw+Ldi6Si0WIWLBjAfJK8+xWt+blIf+552sKRBG589JhzbOvV5ckSXgz8wRmfPY/VIgSxvRsi1WTUxDsp22Wp0XU2lSPS5bFccmIyAsxSUZEREREyrucXZkYq02yISf3Z4iwLxclEWeLz9qVWU6fAcQ664siLGdynO7SJkpIX3cIizcfAwA8dHsHLL2nF3QaXiITNVX1uGQnCooVjoSIqPE4uyURERGRDxMEAe3bt5eXPVZYR0BQ2SfKBDXaxfSDar/KLlGmElSIC4qzW10X3x5QqewTZSoVdO3bOdyd2WrD9E9/wdeH8iEIwDMju+LPtyU061Miao0S2d2SiLwY/yYjIiIi8mEGgwGnT5/G6dOnYTAYlA7HuZBY4M7FNfcFNTB6KYzGXpjTf45crBJUSB+QDmOA0W51rdGI6KefqilQqRCzcAG0Rvt6AHDVZMH97+7E14fyoVOr8MYf+zBBRtRMkqIru1ueulgKq028Rm0iIs/CJBkREREReYZe99YsT90J9JkAABjTcYxcnDE2A3cl3eVw9dC0NHm5w8YNCB03rl6dX6+YcPdbP+LnM1cQ7KfBBw+kYGSPmOaJn4jQNsQPATo1KkQJZy5xhksi8i5MkhERERGR5wlu67A42hDdoNW10fXrHT5fhLve/BHZF0oRE+KHNQ/fgps7hF9XmERkTxAEJFa1JmOXSyLyNkySEREREfmwsrIy3HTTTbjppptQVlamdDiK+eHERfz+7R0oKC5H5+ggrH3kFnSq+iFPRM0rqXpcMs5wSURehgP3ExEREfkwURTx888/y8utUcYv5zBrzf9gtUm4uUMY3v5TP4T4a5UOyycZtAYcmHhA6TDIEV0AML+wRXbFJBkReSsmyYiIiIjIJ0mShP/7/iQWfXUUADCqRwxe/X1P6DVqhSMj8m1J0dUzXBYrHAkRUeMwSUZEREREPscmSnh+/WGs+vE0AOAvtyVg7ohkqFSCsoH5khZsmUTu4a6Wf0lRlV2ZT14oRYVNhEbNUX6IyDswSUZEREREPufxf+/DumNXAQBPj0zGXwZ2UDYgolYkNtQffloVzFYROZdN6BAZqHRIREQNwiQZEREREfmcrw/lQ6f3x6u/74nRPR3PlEnkVg5a2rWWMdtUKgGJUYE4eK4IWQUlTJIRkddgu1ciIiIi8gnnC83ycpCfBqv/nMIEGZFCqrtcnuDg/UTkRdiSjIiIiMjHRUREKB1Cs7OJEnaduoyCYjOigvwQ5KfBw+/uwP9VPf7PB/ojOSFc0RiJWrPEKA7eT0Teh0kyIiIiIh8WEBCACxcuKB1Gs9p8MBcL1h9Gbq2WYwIAXUW5fL9zdJACkRFRtaTqJBlbkhGRF2GSjIiIiIi8xreH8zD9kyOQ6pTXvU9EykqKruluaRMlqDmzLBF5AY5JRkRERERe44WvjjpNiPEnOJHniGvjD51GhfIKEeeulCkdDhFRg7AlGREREZEPKysrw5133gkA+Oqrr+Dv71/zoKUUeKFqYPu55ytn4/NweYXlAHQOH2NrslauFc8m6Yk0ahU6RATgaF4xsgqK0S7coHRIRETXxJZkRERERD5MFEVs3boVW7duhSiKSodDRK1IdZdLjktGRN6CSTIiIiIiIiJqdp3kGS6ZJCMi78AkGRERERF5jbAArdPHOCYZkWdJiq5Mkp0oKFY4EiKihmGSjIiIiIi8hsbJDHlMkBF5nsSomu6WksRRA4nI8zFJRkREREReo6DYgmA/DaKD9XblxhA/LPtDL2WCIiKH2ocboFULMFlsOF9oVjocIqJr4uyWRERERORVlvy+FwZ3icKuU5dRUGxGVJAfUhLCIJjLcEzp4IhIplWrkBARgOP5JcjKL0ZsqP+1VyIiUhCTZEREREQ+zmAwKB3CdbHYamblHN0jBkO6RgMABnQMt6vHuTuJPE9SVBCO55fgREEJUjtHKR0OEZFLTJIRERER+bCAgACUlpYqHcZ1eTszW16eMyJZwUi8m8pgQPLRI0qHQa1MIme4JCIvwjHJiIiIiMhjHTpfiHe2nZLvhxqcz25JRJ6neobL45zhkoi8AJNkREREROSRrDYRs9fsR4XIWfGIvFVS1QyXJ/I5wyUReT4myYiIiIh8mNlsxsiRIzFy5EiYzd41u9zbW7Nx6HwRQgwcIYTIW8VHGKBWCSgur0B+UbnS4RARucQkGREREZEPs9ls2LRpEzZt2gSbzaZ0OA2WVVCC17ecAAA8NaKrwtEQUVPpNWq0D6+cPCSLXS6JyMMxSUZEREREHuepjIOw2ET8pksURnY3Kh0OEV2HJA7eT0RegkkyIiIiIvI4B369iiA/DZ7/XXcIgqB0OER0HarHJcsqYJKMiDwbB3ggIiIiIo/0zMiuMIb4wWQ1KR2KzKA14MDEA0qHQeRVqme4PMHulkTk4ZgkIyIiIiKPIIqS3M3h1o7hGN/vBkXjIaIaKoMByUePNGndxKrulserZrhk61Ai8lRMkhERERGRR/hkVw7uq1peMLYbf0j7uOtJupB36RgZCJUAFJZZcbHEgsggvdIhERE5xDHJiIiIiEhxZy+bsOTb4/L92FB/BaMhoubkp1WjXRhnuCQiz8ckGREREZEPCwgIgCRJkCQJAQEBSofjkCRJeOLz/TBZbEqHQkRuklg1eP8JDt5PRB6M3S2JiIiISFGf7j6LH7MvoY2W/98SVfO17qhJ0YH49kg+svKZJCMiz8UrESIiIiJSzPmrZXh+Y2UiYPpvOikcDRG5S1LV4P3sbklEnoxJMiIiIiIfZjabMX78eIwfPx5ms1npcOxIkoS5XxxASXkF+rQLxZ9ubq90SETkJknsbklEXoBJMiIiIiIfZrPZsGbNGqxZswY2m2eN+bV27zlkHrsAnUaFxeN6Qq3ibJZEvqpjVOWYiBdLLLhcalE4GiIix5gkIyIiIiI7JqsJ3Vd3R/fV3WGymtyyj4IiMxasPwQAeGxIEhKrumIRKaV6DLDko0egMhiUDsfnGHQa3NCmctZatiYjIk/FgfuJyGuIJhOO9ekLAOi8dw8vYImIvJQkSXg64yCKzBXoHhuCBwd2cFnfoDXgwMQDLRRdLboAYH5hy++XyEclRQXi1ytlOJ5fjJSEMKXDISKqhy3JiIiIiKhFbdifi/8czodWLeDl8T2gUfOSlKg1SIrmuGRE5Nla5Ipk+fLliI+Ph5+fH/r3749du3Y5rfvOO+9g4MCBaNOmDdq0aYMhQ4a4rE9ERERE3uNSSTnS11V2s5w6OBFdjMEKR0RELSWRM1wSkYdze5LsX//6F2bMmIH09HTs3bsXPXv2xLBhw1BQUOCwfmZmJv74xz/iu+++w44dOxAXF4ff/va3OHfunLtDJSIiIiI3m7/+MC6XWtDFGIRHUhOVDoeIWlBSdZIsny3JiMgzuT1JtmTJEkyZMgWTJ09G165dsWLFChgMBrz//vsO63/00Ud45JFH0KtXL3Tp0gXvvvsuRFHEli1b3B0qEREREbnR14fysP5/56FWCXh5XE/oNOxmSdSaVLckKyguR6HJqnA0RET1ufXKxGKxYM+ePRgyZEjNDlUqDBkyBDt27GjQNkwmE6xWK8LCHA/sWF5ejqKiIrsbEREREVUyGAwoKSlBSUkJDApOeHLVZMHTGQcBAA/d3gHdbwhRLBYiUkaQnxYxIX4AgBMX2OWSiDyPW5NkFy9ehM1mQ3R0tF15dHQ08vLyGrSNJ554Am3btrVLtNW2aNEihISEyLe4uLjrjpuIKokmE450ScaRLskQTSalwyEioiYQBAEBAQEICAiAIAiKxfHshiO4UFyOjpEBePQ3SYrFQUTKSmSXSyLyYB7dxv3FF1/Ep59+ii+++AJ+fn4O68yZMweFhYXy7ezZsy0cJRERERG58t2xAny+91cIArB4XE/4adVKh0RECulUNcNlFme4JCIPpHHnxiMiIqBWq5Gfn29Xnp+fD6PR6HLdV155BS+++CK+/fZb9OjRw2k9vV4PvV7fLPESERER+Zry8nI89NBDAIC33367xa+bisxWzF17AADw51sT0Ld9mxbdPxF5FnnwfibJiMgDubUlmU6nQ9++fe0G3a8ehH/AgAFO11u8eDGeffZZbN68Gf369XNniEREREQ+raKiAqtXr8bq1atRUVHR4vtftOkocgvNaB9uwMzfdm7x/RORZ0mKrkySncjnmGRE5Hnc2pIMAGbMmIGJEyeiX79+SElJwdKlS1FaWorJkycDACZMmIDY2FgsWrQIAPDSSy9h3rx5+PjjjxEfHy+PXRYYGIjAwEB3h0tEREREzeSHExfxya4cAMBLd/eAv47dLIlau8TIyu6W5wvNKDZbEeSnVTgiIqIabk+S3XPPPbhw4QLmzZuHvLw89OrVC5s3b5YH88/JyYFKVdOg7a233oLFYsG4cePstpOeno758+e7O1wiIiIiagal5RV44vP9AIA/3dweN3cIVzgiIvIEIQYtooL0KCguR/aFUvSKC1U6JCIimduTZAAwbdo0TJs2zeFjmZmZdvdPnz7t/oCIiIiIyK1e/voYfr1ShthQfzxxZxelwyEiD5IUHYiC4nJk5RczSUZEHsWjZ7ckIiIiIu+z69RlrPrxNADgxbu7I1DfIv/LEpGXSIqq7HJ5goP3E5GH4RULERERETWbMotN7mZ5T784DEyKbPjKugBgfqGbIiMiT5HIGS6JyEOxJRkRERERNZvXvj2OUxdLER2sx9yRyUqHQ0QeKElOknGGSyLyLGxJRkTURKLJhGN9+gIAOu/dA5XBoHBERET1GQwGFBQUyMvu9EvOFby77SQA4IXfdUeIP2etI6L6kqIru1uevVwGk6UCBh1/lhKRZ2BLMiIiIiIXRJMJR7ok40iXZIgmk9LhNJogCIiMjERkZCQEQXDbfsorbJi9Zj9ECfhd71j8JjnabfsiIu8WFqBDeIAOAJBdUKpwNERENZgkIyIiIqLrtiLzJLIKShARqMO8UV2VDoeIPFwiu1wSkQdikoyIFGWymtB9dXd0X90dJqv3tdAgIvJ05eXlmDp1KqZOnYry8nK37eed7ZXdLJ8d2w1tqlqIEBE5kxTNwfuJyPMwSUZE5AO8vTsYEblPRUUF3nzzTbz55puoqKhw235sooSR3WNwZ/cYt+2DiHxHUlTluGRZ+UySEZHnYJKMiIiaHZN2RK1PG4MW88fcqHQYROQlqme4PMHulkTkQZgkIyIiIo/CJKv3OF6rBcjcEcmIDNIrGA0ReZPqGS5zLptgttoUjoaIqBKTZERERETUaBU2EU9nHJDvj2Q3SyJqhIhAHUINWogScPICZ7gkIs/AJBkRtQxLKTA/pPJm4YUQEZG3e2/7KRw4VyjfFwRBwWiIyNsIgiB3ueQMl0TkKTRKB0BEREREns8mSth16jIKis2oECW88p9j/LeViK5LYlQQdp++ghOc4ZKIPASTZNSiRJMJx/r0BQB03rsHKoNB4YhaH5PVhP4f9wcA7Lx3JwxavgdEROTa5oO5WLD+MHILzXblPY1BCkXkHVQGA5KPHlE6DCKPJbck4wyXROQhmCQjIiIi8mH+/v44deqUvNxY3x7Ow/RPjkBy8NjRPHaRIqKmS4pmd0si8ixsJU9ERETkw1QqFeLj4xEfHw+VysWlX9F5h8XPff2TwwQZANQehcyan9/kGImodUqKqmyNevqSCZYKUeFoiIiYJCMiIiJqvfZ9XLO8PAXY+wEAYF32OrnYFPUitCG7Ha7+m5yf5eWTI0fh6po17omTiHxSdLAeQXoNbKKE05c4sRMRKY9JMiJfxhkliYhaPYvFglmzZmHWrFmwWCw1DxSeA76aXXNfEoH1jyEvbx8W7VwkFwuCBH3MWgiaQtQWUXYVj+zPqCkQReTOS4c1L89Nz4SIfI0gCEiM5rhkROQ5mCQjIiIi8mFWqxWvvPIKXnnlFVit1poHLmdXJsZqk2zIyf0ZIuzLBUGCSnfRrqxtyUWo6nbEFEVYzuQ0Z/hE5OPkwfs5LhkReQAmyYiIiIhao7COgFDnUlBQo11MP6jqXCJKkgDREmFXdj4wAqLdqGQAVCro2rdzR7RE5KOqxyVjSzIi8gRMkhERERG1RiGxwJ2La+4LamD0UhiNvTCn/5yaYqhQnnsXpIoQu9Uv+YfizR5pNQUqFWIWLoDWaHRz4ETkSxI5wyUReRAmyYiIiIhaq1731ixP3Qn0mQAAGNNxjFwsnZ0Ja+FNCNJr7FY1hvhhzOy/yPc7bNyA0HHj3BsvEfmc6u6Wpy6Wwmpr+AyXosmEI12ScaRLMkSTyV3hEVEro7l2FSLfJZpMONanLwCg8949UBkMCkdESjNZTej/cX8AwM57d8Kg5WeCiFqJ4LYOi4tLA9HjhhD8+6EB2JtzFQXFZkQF+SElIQyCuQzHquppo6NbLlYi8hltQ/xh0Klhsthw5pIJiVVJMyIiJTBJRkREREROadUCFo/rAb1WjQEdw+0ea3ibDyIix1QqAYlRgdj/ayFOFBQzSUZEimJ3SyIPwmbjRETkCS6VWuTlh27viC7GYAWjISJfx8H7ichTsCUZERFRK8eu577N398fBw8elJcb4vmNR+TlKbcnuCUuIqJqSfLg/UySEZGymCQjIvJgHCONiK6XSqXCjTfe2OD6mw/mYvPBPAR1qbyvVbPjARG5V/Xg/UySEZHSmCQjIqIGY9KOvJKlFHihalD6uecBXYCy8XiwqyYLns44pHQYTaYyGJB89Mi1KxKRR6nubpl9oQQ2UYJaJSgcERG1VvxrkIhaJ0spMD+k8mYpVToaIiK3sVgsmD9/PubPnw+LxeKy7sL1h3GxpBwdI5lIJKKWE9vGH35aFSwVIs5e5ri8RKQcJsmIiIiIfJjVasWCBQuwYMECWK1Wp/Uyjxdg7S/noBKA537XrQUjJKLWTq0S0DGSXS6JSHlMkhERERER5n95GADwwG0J6HlDqLLBEFGrUzMuWbHCkRBRa8YkGREREREhv9iM+HADZgztrHQoRNQKJUVXjkt2Ip8tyYiUJJpMONIlGUe6JEM0tb7uz0ySERF5Ak8bI83T4iGiFvHS3T3gr1MrHQYRtUKJnOGSiDwAk2R03UxWE7qv7o7uq7vDZG19mWYiIqqDSVavUVpeIS/fm9IO/TuEKxgNEbVm1d0tTxSUQBQlhaMhotaKSTIiIiKiVuq1b47LyzOGdlIwEiJq7dqFGaBTq1BmteHc1TKlwyGiVopJMiIiIqJWaNepy/hoV458P0CvUTAaImrtNGoVOkQGAKhsTUZEpAQmyYiIiIh8mJ+fH3bt2oVdu3bBz88PAFBmsWH2mv8pHBkRkb3qccmO53OGSyJSBv8yJCIiIvJharUaN910k13Zkm+O4fQlE+KD/ACrQoEREdWRFBUEINfrBu8XTSYc69MXANB57x6oDAaFIyKipmJLMmocDsZMRETk1fbmXMF7208BAOaP7apwNERENZKiOcMlESmLLcmIiIiIfJjFYsGyZcsAAH+dOg2z1+yHKAF39Y5FaqcohaMjIqrRqSpJdiK/GJIkQRAEhSMicj+2RPQsTJKRzzFZTej/cX8AwM57d8Kg5UmGiIhaL6vVitmzZwMALJ1+gxMFJYgI1GPe6K5gX0si8iTtwwOgUQkotdiQW2hG21B/pUMiolaG3S19nGgy4UiXZBzpkgzRZFI6nFbJZDWh++ru6L66O0xWvgdERKScd7edBgA8l9YNoQadssEQEdWhVauQEFE5wyW7XBKREpgkIyIiImolbKKEkT1iMLybUelQiIgckscl4wyXRKQAJsmIiMjnsVUtUaVQgxYLxtyodBhERA6ZrCZsq5iEoOQncTT/otLhEFErxCQZERERkQ87ll8kLz81MhkRgXoFoyEiapjsC6VKh+C1+OcgUdMxSUZERERux/EZlVFhE/H0Fwfl+yO7xygYDRFRw50oqJzhkoioJTFJRkREROSj3t1+CgfP1bQkEwRBwWiIiBqu2GzDheJypcMgH8SWduSKRukAiK6LpRR4oW3l8tzzgC5A2XiIiIg8RPaFEiz55jgEjRbpK/6F1M5R8PPzUzosIqIGX8NnFZQgKpjnLSJqOUySETUnJu2IiMgD2EQJs9fsh6VCxKAuRqRPvomtyIjI62TlF+PWxAilwyCiVoRJMiIiIm/HBD3VsfrH09hz5goC9Rosuqt7oxNkBq0BByYecFN0REQNk1VQonQIRNTKMElGRERE5EPOXCrF4q+PAgDmjOiCqAANli9fDgB48MEHodVqlQyPiKjBmCQjopbGgfuJiIio1fHV2TZFUcKTnx+A2SpiQIdw/PGmdrBYLJg2bRqmTZsGi8WidIhERA12gkkyombhq9c97sAkGREREZGP+HhXDnacvAR/rRov3d0DKhXHISMi7yQIwOVSCy6WcIZLImo5TJIRERER+YBzV8uwaNMRAMDs4Z3RLtygcERERE13Q6g/ACArn63JvJ1oMuFIl2Qc6ZIM0cRWTOTZmCQjIiIi8nKSJGHO2gMotdjQr30bTBwQr3RIRETXJTE6EABwoqBY4UiIqDVhkoyIiIjIy63Z8yu+P34Beo0KL41jN0si8n4dIyuTZBy8n4haEpNkREREPogDtLYe+UVmPLvhMABgxtBO8g9LIiJvJifJFOxuye9SotaHSTIiIiIiLyVJEp764gCKzBXoeUMIHrgtQemQiIiaRceoAABsSUbeiQlW76VROgAiIiIiapp1/zuPb48UQKsWsHhcT2jU9f//1Ov12LBhg7xMROQNOkRUtiS7WFKOK6UWtAnQKRwRkZewlAIvtK1cnnse0AUoG4+XYZKMiIiIyAtdLCnH/HWHAAB/uyMJnY1BDutpNBqMHDmyJUMjIrpuBp0asaH+OHe1DCculOCmgDClQyKiVoDdLb0Qm24SERFR+peHcMVkRdeYYDyc2lHpcIiIml1StPLjkhFR68IkGREREZGX+epALjYeyIVaJWDxuB7QOuhmWc1qtWLVqlVYtWoVrFZrC0ZJRHR9kqKqZ7gsVjgSImot2N2SiIiIyItcKbXgmS8ru1k+PKgjusWGuKxvsVgwefJkAMD48eOh1WrdHiMRUXNIiqrsRn7C3YP3cwwnIqrClmREREREXuTZDYdxsaQcSVGB+NtvEpUOh4jIbRLZ3ZKIWhiTZJ7OUgrMD6m8WUqVjoaIyGNwfEZqjf57NB9rfzkHlQAsHtcDeo1a6ZCIiNwmsaq7ZV6RGUVmdhfntY9zfG2ouTBJRkREROQFCsusmLP2AADgLwM7oHe7NgpHRETkXsF+WhiD/QC0QJdLoqZiwxafwiQZERFRCxNNJhzpkowjXZIhmvhvJzXMok1HkF9UjoSIAMwY2knpcIiIWkT1DJcn2OWSiFoAB+73ITZRwq5Tl1FQbEZUkB9SEsIgMB6PjMdRLGqVcq9Oi8dTdB6ISKpXnG/KR0JIQoPisebnQ5+Q4L4Ya/GG98vnj61GDKjbKl8fL43HG44tT4knv8iMT3efhSAAL93dA35adrMkotahQ2QAtmVdxFcHcxEXZvCoc7PS36OMh/E0dyxKHlueokWSZMuXL8fLL7+MvLw89OzZE//4xz+QkpLitP6///1vPPPMMzh9+jSSkpLw0ksvYcSIES0RqtfafDAXC9YfRm6hWS6LCfHD/KEJaMd4PCoeZ7Gkj+6K33ZwPUOZV8ez7+Oa5eUpwOhlQJ8JWJe9Ti5Oy0jD3e2n46sd7R3Gk7L/O7ns5MhRiFm4AKHjxjVfjA54y/vFY4vxeFs83nJseVI8AJDaKRIpCWEtHg8RkRI2H8zFF3vPAwC+O3YB3x274HHnZn6vMx5vjMfTrns8idu7W/7rX//CjBkzkJ6ejr1796Jnz54YNmwYCgoKHNb/8ccf8cc//hEPPPAAfvnlF6SlpSEtLQ0HDx50d6hea/PBXDz84d56F9J5hWZM/3Qf4/GgeL49nOc0loc/3ItvDuf5ZjyF54CvZtfcl0Rg/WPIy9uHRTsXycUiRHx2einySvPrxfP0O/9F3rPP1RSKInLnpcOa577XzJveLyWOLU+Lx5OOdcbjmjcdW0rE4+y9AoDMYxew+WBuo7an1+vx2Wef4bPPPoNer2+uMImI3Kr63Fx3wH5P/K7gdZhnXWcwHtc87brH07i9JdmSJUswZcoUTJ48GQCwYsUKbNy4Ee+//z6efPLJevWXLVuG4cOHY9asWQCAZ599Ft988w3eeOMNrFixosH7LS0thVpdvyuCWq2Gn5+fXT1nVCoV/P39m1TXZDJBkiSHdQVBgMFgaFhdqwmGWvfLyspQWl4KsVwEABQXl2De53tgs5RXxqGreW42azlsFWaYxMq6X+zOhqSvidGvVgyWcjNEm+j0+VXXVVeYMLhCgk0E/rM7Gza1ARbRLMezdmc2Xt5yHtXPRqqwQhJtNTHZyu3i0YWGQRAqm3RaLeWwVdTUrUvn5weVqjKva7VaYLNWQG0z4beWyr3VjkcSJQgqAZv252Lhl0dhs9T/oVEdj02SoBYEfHUwFxa1GjZrhdMYtHq9/LmqsFpRYbX/0q4dz3f7zgJ+wbDYzJAqJEg2CU+t2QObpX5uWtBoIajUmLfuEN6VJFgkqd77VU2j1UKj1VbGX1EBq8XiNF69YMXoquUN//sVZpva7v2qHY+gVkNQV25XFG2QKqx4as0evO3g86PWaqDV6qrqirCYHb++1XWNV3/BAKnO50uyYf2OzRBhXy4IElS6i6iwBkOylsvlMVfPQ6h7nIgitn63F2Vde6K8rMxpDCq1Cjp9zbFhNpkcfnYAoAI17+kz6w7B5uL1ffrzvfigavmrg7koEx0fxwAgqATo/WreT3NZGVDr+djHcxLaoHAAgMXm+P2q2TCg0vrJx5xZFJ1+dgD7477cXAbJScxqmwm/q/3cbCqINtHp56f2uUessMBmdX7u0fv7N/i4r66rrjDhjgoJFa7OPd+eqzn32KyQbM7PPdqQNvXOJ87UPu4bWlcUJTzzxf9cnnsqJAmaqnOPVa2pdz6x265OB7Wm8ivb0XFf+/Pz319yIPiHVJ57bBKkCglrd2bjpc2nYLPY70NQawC1BhIAmySh3MW5p/Zxb7PZYC0vr1enmk6wYEzV8sb951BWoWrYuUcSIVktzs89GjW0usokjyRJLo/72nVV1lIMcvBdIcfz+S+QUOszXOd9qx3Pl3tOQxsSKj9mdjGmXL3jvgF1RVHCvHWHYLOaAQeHpwBg3ud7MbTrCLlLREOuOcaPHw+g8jqivPZ7ZykFql4blJYioFZX5rKyMoii82uDgICauuWiiNLSUqicxFG7rtlshs3m/Lg3GAzyOaK8vBwVFc6PucbU9ff3l497i8UCq4tjrjF1/fz85HNEY+parVZYXH2H6/XQVB33jalbUVFh/x7XodPpoK26jmhMXZvNBrOL73utVgudTtfouqIooszFsdyYuhqNRk4ES5IEk4tjrjF1G/P7wet+a9Sp25jjvjF16x33Ls49okUEJODZL/fBZqkft6DVQxAEzFt3CP8niqhA/WuNag295lDbTBgrSVBVfS+bJQ1s1grn311aHQShclm0WV1e9zTlOkJdYcKQL/rAYgP+M3K7fJ1aO56nv9gPCZXHp2SrgGSr2W7d6x5NUIjL64jaGvNbo7quKEp4JmO/y+seiyRBV/X6Vmh1Lq8jGvtbQ6vVVV73iBIkq4vrHpUa0GghofK4L3N13dOIaw6tVI6xVctfHcxFqQX1rlOr4xFUKggaXc1rYzHXe7+q42nKdQRQ+Zt9kFWCJDm57lmzFxJqtivWuuYQADybsQ/LquIxmUwIVOIccR11nV0buDrH2pHcqLy8XFKr1dIXX3xhVz5hwgRpzJgxDteJi4uTXnvtNbuyefPmST169HBY32w2S4WFhfLt7NmzEirfYoe3ESNG2K1vMBic1h00aJBd3YiICKd1+/XrZ1e3ffv2Tut27drVrm7Xrl2d1m3fvp0kpQdX3spLpH79+jmtq/IPlto/sUG+6eO6Oa0raPV2df07ON8uALlelyfWSOO6alzWjfv7Grl+QLffuKx7w98+kusG9h7psm7sX9+T6wan3OWybuLziVK3Vd2k9k+ulUJu/aPLuv9q11463LmL1Onxz6XQ1Mku60b/8QU5hrChf3X93MbNqaz75Fop9oFYl3Ujxj4ptX9ig9Tp8c+lJW3buqwbPuIxOYbIcemu4x36gPz5ifvjfJd1Q1Mny9s1Tljism7IrX+U68b8ebnLusEpd0n9n1gtVcwLqfkspwdL1meCpBvitNKN798odVvVTb7duLK7FP/UR9INf/vI/rloNNKBTp2lw527yLcDXZKlvo9+KMX9fY3LGAydb7X7vLuq69+xjxxL+yfXSoJW77zuDTfKsXR6/HNJ5R/stK7OmGQXgzo4ynnd8Btq6j65VtK3dR6DOjhK/uwc7txF6ubn57RuY88R1e9VlyfWNPgc0f6JDZKh860u6zblHNHliTXSI/20Lus25hwR8+flct1rnSOME5bIdZvzHPFW7A3yZyd8xGMu61afI9o/sUGKGPuk63hHTJU/O+0fc/5dBEAKG/pX+fOzKi7OZd3GnCPCbx0vf34S/vyay7rBKXfJ243963su6wb2HinXrXuOqHsL6PYbuW6nv3/osq6h84AGnyMCEvrY1XV1jtDHdbOr21znCG14O+nHExcbeB3R3u6aw9V1REREuF3dQYMGOX/NDAZJkiTJVloqHe7cRbo9IMDl61bbuHHjXNYtKSmR606cONFl3YKCArnuI4884rLuqVOn5LozZ850WffgwYNy3fR019+1u3btkusuXrzYZd3vvvtOrvvGG2+4rLthwwa57sqVK13W/eyzz+S6n332mcu6K1eulOtu2OD68/7GG2/Idb/77juXdRcvXizX3bVrl8u66enpct2DBw+6rDtz5ky57qlTp1zWfeSRR+S6BQUFLutOnDhRrltSUuKy7rhx4+w+w67qet9vjfZ2dV2fIyLs6jbkHFFtxIgRLl+3aqWWUim4n/NzJVBzHdHp8c+ltGDXdRvzW+PU9ED5ukfp64guT6yR3rjT+TUdACly3Fx5u815HdGY3xrV1xHtn9ggRf/xBZd1H4+MlK97mvu3RvV1T6eXO7msW30d0enxz6XtHRNd1q19HXGt3xpBnW+2u252Vde/Qz+3X0d0eWKN1D5EcFpXGxFnt11teDunddu3a+dR5whJuv7riMLCQskVt7Yku3jxImw2G6Kjo+3Ko6OjcfToUYfr5OXlOayf56RL1aJFi7BgwYLmCbgVUQkCbu5QM6bJ9wYtnOfGIdfVi2Y05zxsfduHwi+ocgr7PcF6nHBRt3dcKAIiKuPYt9cPRQ3YfsdIA/Y0MJaUhDAcOGXAVRd1usYEI6rqtcg6E4DLLup2igpCZIcwiCjHjxcaGEQDdIwMQEJVDOdLA+Fq0+3DajLuA2MkVI8KpmmjgT5aj/L8clRccf5vFgBEazRor9XhjNWC/Kos/A1t/NGtKoZCv6tw1fEnJsQP8R2SsKrkEfy58A0IgoAKUcJDG8z49awVpavOo+3EthDUAiQRKM+7C1JFCIBCu+3kV1QgPT8PC6KNUAsCbAA2DpmIxBs7IL68DGddxBAWoLP7vJ9xUTfEXwvA+T9mzqQkhOGMWgVn/4kF6NV2MazXqpweS/66mroiyrHOSb3G0qpVdjH811/r9POjFmoG7UyJD0PhNc4RjZGS0Aaaqn/Idgbp4eo/nepzhF50/g9iU/SMC0FI28rX4uAh/zqfNnvdYoMRHl9Z92hWw84RF0vK8UsDY0lJCMPxiwG45KJOp+hAxFW9d2evBOKii7odIgIQV3Xu2XWkgUE0QPswA7pUxXBJFQxXDfFjQ2v+neweG4xTzRSDMViPvlUxmIsF/OqibmSQHv2r6qrKBBxvphhCDFq74+jXqvORI8H+9nXzGnCOuFhSjhMF1/6ns6C44cdERUUFvvjiCwBw2pqEiMgXNea3RrWU+DDsDGnYb42GaMp1hF40A/ubKQA07jqiMb819NFfIih5O6IuvYJDOQ2LJSUhDHm6Ky6vIxrzW+OGYA26V133/K+53jTYX0dUNOC3RrWU+DA4znQ0XlOuI4DK3+xbmymG1kiQ3Hi1dP78ecTGxuLHH3/EgAED5PLZs2dj69at2LlzZ711dDodVq9ejT/+8Y9y2ZtvvokFCxYgPz+/Xv3y8nK7ZuJFRUWIi4vD+fPnERwcXK++1zWBtppgWJpYeWfueZTZVCgpL0HqZ6kAgFdSvsBDHxyoiaN2lydrOfQVZnz6VWUSsdMP26Gqtd8mNW+0lMK8MAY2EcCsE4AuACaryS6ev35yXG7eWLe7pd5WbhdPYETE9XWTsJQCL1e9PrXiGfzFYAgqAf8Y8BUmvvOzXber2vS2cnyxeSHUgoDOe/egQqO5vm4SteLRP30WGkMITFYTUj5IgWSTUHzsaUDS1d2s3N1SX1GONevnwiJJ9d6vao3pJqHf/09ovn4CACAJKlh+uxifGnR4+ZeXIQgCJElA2a+jYb3cx67Lk1TV3XLUmZ8w7eD6yvdIpULYU3MRmJbWtG4SllKIz8egvAKQHtoGKbzydZI/PwJQevLvgBRTGYMk2XW3BCrfr399tQACgMT16xDUubNct9HdJBx8dgCgXCzHoLWDAADFRxfKzZIdvr42C9Z9vRAA0HnvHpS5OJ1e8xxRKx5hdjYMoZGV9awm3LTypsp4HH1+qrpb6ivKkbHhKZhFEYlOPjtAI5pAW0oRsKxqBtK552EW1bDZbHbHe+14ap97pAoLdNYyp+eeJnWlspSifGEMKlyee47JXR/qdrese+4JCA93a1eqHdmX8IcV2+26PtSmt5Vj7eaF0FSde2xa7fV1par1+dE9lQNtQGjlueefKZAqpHrfFdUEtQaCWgN9RTk+Xz8X5S7OPY3qSiVZoHslHgAgPvkrymwqp58du3NPVXfLuu9XdTxN7kpVXgLTgqrZUR18fkqy5kFQB8rr1u1uWTue5J92wBBWc8Ha3NcRO7Iv4Y/v/GTX9aEeAfjXI6kY0LGyW/a1rjkkSUJgYOXzu3Dhgl1MdTW2m4RoMuFYn74oF0V0bOC5h90t2d2S3S0recRvDSW7W9a6DgtoE1UZq9WEm1bdBEhOrntQ091SX1GOz9bNQQXqX2tUa/A5wlIK/6WJUAkCMPc8LNDCarU6/+6q1d1Sslmhs5icXvc06RxhKYX12RhYbLC7TrX77speAEHwr4rBvrtl3e9R/zZtmv0cYbKaMPjzwRA0Vb+53t0LqcLxc9PbyrFm80Loqq57JL3++s4RtT4/2rlnoAtsU3nd82EKJKuL6x6VGoJGC31FOb5YPxdlLq57GnWOsJnht6RD5Z2551FqRb3r1Op46na3FC1mp9c9TT5HWEphWhBTObKLg+ue4uPPQKUJqomhzjVH7Xg6//gDAiMi5Me8ubtlUVER2rZti8LCQoe5ompubUkWEREBtVpdL7mVn58Po9HocB2j0dio+nq93uEgtAEBAXYvpDMNqdOUugYnF4mNrlvnesjf3x+SRoJKX3lSHpgci9jIs8grNNe7llZp9VALgKHqIi8gIMDpxWvtL/Nr8dNUtTAJCAB0ARCsgl08bUNr4hE0WghVfeUBQF0h2MUj1Gqt4uy9dESn01WeOLUAdPXjEarGaenXvg3ahgU5fH2q46ndYkbebgNotVr5C6SmsFY8mprDS9AIEDQC2kaGIP+q5DAWAYAx2A8aQYBGEFy+X9U0Go38hVdP4TngP3Nqti+JuLJlDpbEtZVfd0GQ4H/DeojmG6tab1XVVakRaSvG1EMba94jUcTl519A2G9+A22t41GlUjX42FAJAvy1AIyJ8pd97c+PMaSt/PoIggChVuJFANDW4A//6s9PXFzNY1WvV0MFBAQ4/OxUx1PNGKJ3+X61Ndj/0Lyuc0TteOo8Vv36XOvzAwB+Ve/HtT47AFz+UEadj3b1OaL2++UsHkGjgxpSg849jTnu9RoB+soNXvvco9bKiReg/rmn+scv0LjjvqF1UxLC0DYs0OW5R1Pr3OPwfOKEw+O+9uen1nYEtQBBLbj8rgAqPz9qQYChgecetVrt+vNe67tLpVIhwC+gYZ8dQQWVzg9tDf7X/Pw05rgXBAEBDr4rquOJCQ+yi0fl4NxTHU/d78vmvo5ISQhDTIgf8god58gEAMYQP7sZLq91zVH7otrf37/BMbs8R9Shb8S5pzHXHI06RzTlOkLBuo057q/7HNEMda953DexbqOuIxpRt9HnCDf9JvD43xp1NOa4b0zdesd93euwWlS6hl33GIP9oFOpoIPra41qLs8RWgAOfhM06LtLrYVa69eg655GnSPUArRq1LtOlb+72hhqrpur/vSqVu+6p9Zx3lznCMFa+fsGqPrN1SbA5XWPrtbre93niNqfn1qvp6ASIOgbdt0jNOK655rniDq/2QMC6l+nOv3NrvOr/3414A+nazFonV/3tI0Itr/u0Tq/7ql7/mixc0Qz1a193LtKrNXm1tktdTod+vbtiy1btshloihiy5Ytdi3LahswYIBdfQD45ptvnNZv7VQqAemjuwKo+ZFcre79lsB4XJt7ZxeH+66+P3dEl+bb2eVsoM5g+TkaFcQ6p2ZBkKDW2Te6FgDEllyEqu5pXBRhOdPAttRN0KKvjw/G09I8KR5PO9bVHhaPp70+3nZstWQ8DfnspI/uKg/aT0Tkqzzp3NyQeFqaJ8XjadcZjMc1Tzu2PI1bk2QAMGPGDLzzzjtYvXo1jhw5gocffhilpaXybJcTJkzAnDk1rV2mT5+OzZs349VXX8XRo0cxf/58/Pzzz5g2bZq7Q/Vaw7vF4K37+8AYYp9RNYb4YdkfejEeD4pnSFej01jeur8PhnZ13GKyScI6AoL9Id6uQoSqzulQgArh+th68UybMBhQ1TlFqFTQtW/XfDHW0aKvj4/Fo8Sx5WnxeNKxznhc86ZjS4l4XL1Xb93fB8O7xbRoPERESvC0c7OnXfd4WjyedJ3BeFzztGPL07i1uyUA3HPPPbhw4QLmzZuHvLw89OrVC5s3b5YH58/JybHr9nLLLbfg448/xtNPP425c+ciKSkJGRkZ6Natm7tD9WrDu8VgaFcjdp26jIJiM6KCKrtiCOYyHGM8HhWPs1jUKgGii77ujRYSC9y5GNg0s/K+oIZxxBLMMejw/M7nAQAqQYX0AekY2/F3DuO5/PRTyF/4bOX6KhViFi6w62rpDi32+vhYPDy2GI+3xeMtx5YnxkNE1Fp4y7mZ3+uMx9vi8bRjy5O4PUkGANOmTXPaEiwzM7Ne2fjx4zF+/Hg3R+V71CpBHsS3mvNh8tyP8TQuFrfodW9NkmzqTiAiCWOsJjlJljE2AwkhCQDgMJ7QtDQ5SdZh4wboExLcHzNa8PVpIG+Ih8dWDcbjmifF4w3HlpI8LR4iIiV42rnQk75HAcZzLYyncbFQC3S3JCIPEdy2XlG0IbrBq2ujG16XiIiIiIiIyNu0SEsyIiIiIo9QdB6ISLIryjfly61qr8Wan99irWqbi06nw8qVK+VlIiKiRrGUAi9U/eE+97w82yaRL2JLMiIiIvJt+z6uWV6eAuz9AOuy18lFaRlpWJu11unqVzMy5OWTI0fh6po17ojSbbRaLSZNmoRJkyZBq9UqHQ4RERGRx2KSjIiIiHxX4Tngq9k19yUReZtmYNHOF+QiESIW7FiAvNK8eqtb8/KQ/9zzNQWiiNx56bDm1a9LRERERN6NSTIiIiLyXZezAcl+SNwcjQoiJLsyURJxtvhsvdUtp88AYp0hdUURljM5zR6qu1RUVGDjxo3YuHEjKioqlA6HiIiIyGMxSUZEVHTeYXG+Kb/Bm7DmN7wuEbWgsI6AYH+5065ChAqCXZlKUCEuKK7e6rr49oCqzuWSSgVd+3bNHqq7lJeXY9SoURg1ahTKy8uVDoeIiIjIYzFJRkStk4MxigC0qnGKiFqFkFjgzsU19wU1jCOWYE7/uXKRSlAhfUA6jAHGeqtrjUZEP/1UTYFKhZiFC6A11q9LRERERN6NSTIian0cjFGE9Y8hL28fFu1cJBdznCKiZuCgpWaLt9LsdW/N8tSdQJ8JGNNxjFyUMTYDdyXd5XT10LQ0ebnDxg0IHTfu+mMiIiIiIo/DJBkRtT4OxiiCZENO7s8QYV/uy+MUNdh1dkdlV9RWyJNnkwxuW68o2hDd4NW10Q2vS0RERETehUkyImp9HIxRBEGNdjH9oKpzWlRknCJPSEpdZ3dUdkVtxTibJBERERF5KSbJiKj1cTBGEUYvhdHYC3P6z5GLW3ScIk9KSl1nd9RWkeSwlALzQypvllKlo/EsnE2SiIiIiLwUk2RE1Do5GKMIgDLjFHlaUuo6u6MyydHKcTZJIiIiIvJSTJIRETkYowhowXGKPC0pdZ3dUZnkaOU4m6TH0el0eOONN/DGG29Ap9MpHQ4RERGRx2KSjJrOE2YsI/IFnpaUus7uqK01yWGymtB9dXd0X90dJqtJ6XCUxdkkPYpWq8XUqVMxdepUaLVapcMhIiIi8lhMklHjePKMZUTeyhOTUtfZHZVJDpJxNkkiIiIi8hJMklHDccYyIvfx5KTUdXZHZZKDSFk2mw2ZmZnIzMyEzWZTOhwiIiIij6VROgDyIo2csaxuaxdX4yb5ejcsokZhUoqImpHZbMbgwYMBACUlJQgICFA4IiIiIiLPxJZk1HCcsYyIiIiIiIiIfBSTZNRwnLGMiIiIiIiIiHwUk2TUOJyxjIiIiKqoDAYkHz2C5KNHoDIYlA6HiIiI6LowSUZNxxnLiIiIiIiIiMhHMElGvqPofL2ifFN+g1e35je8LhERERERERH5FibJyLvt+7hmeXkKsPcDrMteJxelZaRhbdZap6tfzciQl0+OHIWra9a4I0oiIiIiInLEwR/dAP/sJiJlMElG3qvwHPDV7Jr7koi8TTOwaOcLcpEIEQt2LEBeaV691a15ech/7vmaAlFE7rx0WPPq120StmwjIiXw3EN1aLVaLF68GIsXL4ZWq1U6HCIih390A+Cf3USkOCbJyHtdzgYk0a4oR6OCCMmuTJREnC0+W291y+kzgGi/PkQRljM5TY+JLduISAk895ALOp0Os2bNwqxZs6DT6ZQOh4haOwd/dGP9Y8jL24dFOxfJxYr+2U1ErRaTZOS9wjoCgv1HuF2FCBUEuzKVoEJcUFy91XXx7QFVnUNApYKufbumxePpLduIyDfx3ENERN7EwR/dkGzIyf0ZIuzLW/TPbiIiMElG3iwkFrhzcc19QQ3jiCWY03+uXKQSVEgfkA5jgLHe6lqjEdFPP1VToFIhZuECaI316zaIJ7ZsIyLfx3MPXYPNZsPu3buxe/du2Gw2pcMhotbOwR/dENRoF9MPqjo/T1vsz24ioipMkpF363VvzfLUnUCfCRjTcYxclDE2A3cl3eV09dC0NHm5w8YNCB03rumxeFrLNiJqHXjuoWswm81ISUlBSkoKzGaz0uEQUWvn4I9ujF4Ko7EX5vSfIxe36J/dRERVmCQj3xHctl5RtCG6watroxte1yFPa9lGRK0Dzz1ERORtHPzRDUC5P7uJiKowSUbUnDypZRsRtR489xARkbdy8Ec30MJ/dhMRVWGSjMhdlG7ZRkStE889RERERERNwiQZERERERERERG1ekySERERERERERFRq8ckGRERERERERERtXoapQMgIiIiIvfRarVIT0+Xl4mIiIjIMSbJiIiIiHyYTqfD/PnzlQ6DiMg7FJ0HIpLqFeeb8pEQknDN1a35+dAnXLseEXkmdrf0JkXn6xXlm/IbvLo1v+F1iYiIfJlBa8CBiQdwYOIBGLQGpcMhIiIl7fu4Znl5CrD3AwDAuux1cnFaRhrWZq11uPrVjAx5+eTIUbi6Zo1bwiRqNOYQGo1JMk/n4ITd0JM1wBM2ERFRayeKIg4dOoRDhw5BFEWlwyEi8iyF54CvZtfcl0Rg/WPIy9uHRTsXycUiRCzYsQB5pXl2q1vz8pD/3PM1BaKI3HnpsObZ1yNqMcwhXBcmyTyZgxN23qYZWLTzBbnI2cka4AmbiIiIgLKyMnTr1g3dunVDWVmZ0uEQEXmWy9mVibHaJBtycn+GCPtyURJxtvisXZnl9Bmg7h8QogjLmRx3REvkGnMI141JMk/m4ISdo1FBhGRX5uhkDfCETURERERE5FJYR0Co87NYUKNdTD+o6vxcVgkqxAXF2ZXp4tsDqjrrq1TQtW/njmiJXGMO4boxSebJHJyw21WIUEGwK3N0sgZ4wiYiIiIiInIpJBa4c3HNfUENjF4Ko7EX5vSfIxerBBXSB6TDGGC0W11rNCL66adqClQqxCxcAK3Rvh5Ri2AO4boxSebJHJywjSOWYE7/uXKRs5M1wBM2ERERERHRNfW6t2Z56k6gzwQAwJiOY+TijLEZuCvpLoerh6alycsdNm5A6LhxbgmTvIDSA+Uzh3DdmCTzdA5O2A09WQM8YRORD1L64oOIiIh8V3Bbh8XRhugGra6Nblg98iGeNlA+cwjXhUkyb+LghN3QkzXAEzYReTFPu/ggIiIiIvL0gfKZQ2g0JsmIiMizefrFBxEREZG7OGhBDzS8FT1b0LsZB8r3OUySERGRZ+PFB9F10Wq1mDlzJmbOnAmtVqt0OEREdC0OWtADaHAreragb0EcKN/nMElGRESejRcfRNdFp9Ph5ZdfxssvvwydTqd0OERE5IqDFvRY/xjy8vZh0c5FcrGzVvRsQd/COFC+z2GSjIiIPBsvPoiIiKi1cNCCHpINObk/Q4R9uaNW9GxBrwAOlO9TmCQjIiLPx4sPoiYTRRGnT5/G6dOnIdb94URERJ7FQQt6CGq0i+kHVZ2f745a0be6FvSeNus5B8r3ekySEVHL0AUA8wsrb7oApaMhb8aLD6JGKSsrQ0JCAhISElBWVqZ0OERE5IqDFvQYvRRGYy/M6T9HLnbWir5VtKDnrOfkRkySEREREREREXkKBy3oATS4Fb1Pt6DnrOfkZkySEREREREREXkiBy3ogYa3ove5FvSc9ZzcjEkyIiIiIiIiIvJ8nPWc3IxJMiIiInI7g9aAAxMP4MDEAzBoDUqHQ0RERN6Is56Tm2mUDoCI3Kh6sHwPVv3DmYiIiIiI6Jp63Qtsmlm5PHUnEJGEMVYTnt9ZOdZYxtgMJIQkOF09NC0N+QufBVA5Zps+wXldan3YkoyIiIiIiIiIvA9nPadmxpZkRERERD5Mo9HgkUcekZeJiIiarOg8EJFkV5RvynfZcqs2a34+W26RR2NLMiIiIiIfptfrsXz5cixfvhx6vV7pcIiIyNvs+7hmeXkKsPcDrMteJxelZaRhbdZap6tfzciQl0+OHIWra9a4I0qiZsEkGRERERERERH9f3t3H1tlff8N/HMKtFCg7SotpZNnNsApLjdg1236M4PIg/HnAyQ6WW7cGGYbuPi0DXbfDry33SzO/cZ0ZmbJMrJE9uA2mHPTzYhgtiAihogMiBIdTKyIhBZaQORc9x/+ON4dBVqwXqfnvF7JlVznOtfFeRO/fks/5/twsubXIh77+nuvk2w0/fn2WLbh/+YuZSMbd6+/O5pam056/FhTU7zxne++dyGbjde/tSSONZ18L+QDY+6hyFgoH6C4JEkS+/bti4iIQYMGRSaTSTkRAD3G/p0RSbbdpV29SyIbSbtr2SQbuw/uPmlHybdf/WdEtv3zkc3G2//cZUdJ8pIiGQBQ2HrATr/dqa2tLWprayMi4tChQ9G/f/+UEwHQY1SPjsiUtCuUDXsnGyWRaVcoK8mUxNCBQ096vHTE8IiSkvaFspKSKB0+rFtjw9ky3RIAAAA4WeWHI2bc897rTK+om/lfsbjhm7lLJZmSWNK45KRRZBERferqYvD//l/vXSgpiSH/526jyMhbimQAAABAxz5+43vnCzZE/I//Gf85+j9zl1ZfvTqu+8h1p3y86pprcuej/vRoVM2e3R0p4X2hSAYA56Jlz0mX3mh7o9OPH3uj8/cCAKSqov6kS4PLB3f68T6DO38vpEGRDAC6ylboAABQcBTJAKArbIUOAAAFye6WQHEq8t3uOAe2Qj8z/38BANADKZIBQFfYCp0epnfv3jF37tzcOQAAHeu26Zb79++POXPmREVFRVRVVcW8efPi0KFDp73/lltuibFjx0a/fv1i2LBh8dWvfjWam30TDfCBOzESaGnzu+e8x1bo9DBlZWWxYsWKWLFiRZSVlaUdBwAgb3VbkWzOnDmxdevWeOKJJ+LRRx+Np59+Om6++eZT3r9nz57Ys2dP3HvvvfHiiy/GihUr4vHHH4958+Z1V0QAODu2QgcAgILTLWPut23bFo8//nhs3LgxJk2aFBER999/f8ycOTPuvffeqK8/edvYCy+8MH73u9/lXo8ePTq++93vxuc+97l45513Tjk94OjRo3H06NHc65aWlvf5bwPwAbCGU89lK3TyXJIk0dbWFhER5eXlkclkUk4EAJCfumUk2fr166OqqipXIIuImDp1apSUlMSGDRs6/ec0NzdHRUXFadfPWLZsWVRWVuaOoUNPXv8FAKBYtbW1xYABA2LAgAG5YhkAACfrliJZU1NT1NbWtrvWu3fvqK6ujqZObnG/b9+++Pa3v33aKZoREYsXL47m5ubcsXv37rPODVDepzy2zN0SW+ZuifI+5WnHAQAA4APSpSLZokWLIpPJnPbYvn37OYdqaWmJK6+8Mi644IJYunTpae8tKyuLioqKdgfdyGLeAAAAQAHq0ppkd9xxR9x0002nvWfUqFFRV1cXe/fubXf9nXfeif3790fdGXbvOnjwYEyfPj0GDhwYq1atij59+nQlIlDASsrLY/z2bWnHKGonRtoBAAAUmi4VyWpqaqKmpuaM9zU2NsaBAwdi06ZNMXHixIiIWLNmTWSz2WhoaDjlcy0tLTFt2rQoKyuLRx55JPr27duVeADAf1PQBACArumWNcnGjx8f06dPj/nz58ezzz4bf//732PhwoVxww035Ha2fO2112LcuHHx7LPPRsS7BbIrrrgiWltb42c/+1m0tLREU1NTNDU1xfHjx7sjJkDes0YaAADAB6NLI8m64qGHHoqFCxfGlClToqSkJGbNmhX33Xdf7v1jx47Fjh07crssPf/887mdL8eMGdPuz3rllVdixIgR3RWVc2S0AgAAANDTdVuRrLq6OlauXHnK90eMGBFJkuReX3755e1eAwBw7nr16hWzZ8/OnQMA0LFuK5IBPZ+F8gF6vr59+8bDDz+cdgwAgLzXLWuSAQAAAEBPYiQZPVtp/4ilzWmnAAAAAHo4RTIA6Ol8YcBptLa2xoABAyIi4tChQ9G/f/+UEwFFz88tIE8pksH7yQ98gIJjfUYAgOJgTTIAAAAAip4iGQAAAABFz3RLAAAAUlfepzy2zN2SdgygiCmSUXD8cAUAAAC6ynRLAAAAAIqekWTQzYxsAyBNvXr1ipkzZ+bOAQDomCIZAJBXSsrLY/z2bWnHKBh9+/aNP/3pT2nHAADIe4pkAGfJL/IAAACFw5pkAAAAABQ9I8nyXWn/iKXNaacA6BKj7CB/tLa2Rm1tbURE7N27N/r3759yIgCA/KRIBgBQ4Nra2tKOAJC/DEwA/psiGUABMHILAADg3FiTDAAAAICiZyRZD1Tepzy2zN2SdoyzYrQLAAAAkI8UyQAAAAAKgTX2zonplgAAAAAUPSPJAAAKWElJSfzHf/xH7hyAzunJy9wAZ0eRDACggPXr1y/Wrl2bdgwAgLynSAZA/rO2AgAA+ci/UwuKMfcAAAAAFD0jyQAAClhra2uMGDEiIiJeffXV6N+/f7qBADg9I5MgNYpkFLWS8vIYv31b2jGAs2AxXei8ffv2pR0BACDvKZIBQJHzhQEAQOHy5XLnKZJBHvGLKgAAAKRDkQwAAIAexZfLRcyabXQju1sCAAAAUPSMJAMAAIA8Z12pnsN/q55LkQwAoICVlJTEpEmTcucAAHRMkQwAoID169cvNm7cmHYMAIC85+tEAAAAAIqeIhkAAAAARU+RDACggLW1tcWIESNixIgR0dbWlnYcAIC8ZU2yAldSXh7jt29LOwYAkJIkSeKf//xn7hwACo3dJHm/KJIBAAAAHSvtH7G0Oe0U8IEw3RIAAACAomckGQAAAJwlS9xA4VAkA4D3gbUwAACgZzPdEgAAAICiZyQZAHSVBWzpQTKZTFxwwQW5cwAAOqZIBgBQwMrLy2Pr1q1pxwAAyHumWwIAAABQ9IwkAwAAADrNhkWFq9h3a1UkA6DgFfsPe4pbW1tbTJ48OSIiNm7cGOXl5SknAgDIT4pkAAAFLEmS+Mc//pE7BwCgY4pkAAAAACkw4yG/KJIBAABAgVB0gbNnd0sAAAAAip4iGQAAAABFT5EMAAAAgKJnTTIAgAKWyWRi+PDhuXMAADqmSAYAUMDKy8vj1VdfTTsGAOQFGxtwOqZbAgAAAFD0FMkAAAAAKHqKZAAABezw4cMxefLkmDx5chw+fDjtOAAAecuaZAAABSybzcZzzz2XOwcAoGNGkgEAAABQ9BTJAAAAACh6imQAAAAAFD1FMgAAAACKnoX7AQAAgG5RUl4e47dvSzsGdIoiGQBAgRs0aFDaEQAA8p4iGQBAAevfv3+8+eabaccAAMh73bYm2f79+2POnDlRUVERVVVVMW/evDh06FCnnk2SJGbMmBGZTCZWr17dXREBIBUnph2M374tSsrL044DAABENxbJ5syZE1u3bo0nnngiHn300Xj66afj5ptv7tSzy5cvj0wm013RAAAAAKCdbpluuW3btnj88cdj48aNMWnSpIiIuP/++2PmzJlx7733Rn19/Smf3bx5c/zgBz+I5557LoYMGdId8QAAisbhw4djxowZERHx2GOPRb9+/VJOBACQn7qlSLZ+/fqoqqrKFcgiIqZOnRolJSWxYcOGuPbaazt8rq2tLW688cZ44IEHoq6urlOfdfTo0Th69GjudUtLy7mFBwAoINlsNtatW5c7BwCgY90y3bKpqSlqa2vbXevdu3dUV1dHU1PTKZ+77bbb4pOf/GRcffXVnf6sZcuWRWVlZe4YOnToWecGAAAAoDh1qUi2aNGiyGQypz22b99+VkEeeeSRWLNmTSxfvrxLzy1evDiam5tzx+7du8/q8wEAAAAoXl2abnnHHXfETTfddNp7Ro0aFXV1dbF379521995553Yv3//KadRrlmzJnbu3BlVVVXtrs+aNSsuvfTSWLt2bYfPlZWVRVlZWWf/CgAAAABwki4VyWpqaqKmpuaM9zU2NsaBAwdi06ZNMXHixIh4twiWzWajoaGhw2cWLVoUX/ziF9tdu+iii+KHP/xhXHXVVV2JCQAAAABd0i0L948fPz6mT58e8+fPjwcffDCOHTsWCxcujBtuuCG3s+Vrr70WU6ZMiV/84hdxySWXRF1dXYejzIYNGxYjR47sjpgAAAAAEBHdtHB/RMRDDz0U48aNiylTpsTMmTPj05/+dPz0pz/NvX/s2LHYsWNHtLW1dVcEAAAiory8PMrLy9OOAQCQ1zJJkiRph3g/tbS0RGVlZTQ3N0dFRUXacQAAAABIUWdrRd02kgwAAAAAegpFMgAAAACKniIZAEABO3LkSFx55ZVx5ZVXxpEjR9KOAwCQt7pld0sAAPLD8ePH489//nPuHACAjhlJBgAAAEDRUyQDAAAAoOgpkgEAAABQ9BTJAAAAACh6imQAAAAAFL2C290ySZKIiGhpaUk5CQBA+lpbW3PnLS0tdrgEAIrOiRrRiZrRqRRckezgwYMRETF06NCUkwAA5Jf6+vq0IwAApObgwYNRWVl5yvczyZnKaD1MNpuNPXv2xMCBAyOTyaQd5wPX0tISQ4cOjd27d0dFRUXaceCsacsUEu2ZQqEtU0i0ZwqFtkyh6M62nCRJHDx4MOrr66Ok5NQrjxXcSLKSkpI4//zz046RuoqKCh0kBUFbppBozxQKbZlCoj1TKLRlCkV3teXTjSA7wcL9AAAAABQ9RTIAAAAAip4iWYEpKyuLJUuWRFlZWdpR4JxoyxQS7ZlCoS1TSLRnCoW2TKHIh7ZccAv3AwAAAEBXGUkGAAAAQNFTJAMAAACg6CmSAQAAAFD0FMkAAAAAKHqKZAAAAAAUPUWyAvLAAw/EiBEjom/fvtHQ0BDPPvts2pGgy5YuXRqZTKbdMW7cuLRjwRk9/fTTcdVVV0V9fX1kMplYvXp1u/eTJIlvfetbMWTIkOjXr19MnTo1XnrppXTCwhmcqT3fdNNNJ/XV06dPTycsnMayZcti8uTJMXDgwKitrY1rrrkmduzY0e6eI0eOxIIFC+K8886LAQMGxKxZs+KNN95IKTF0rDNt+fLLLz+pb/7Sl76UUmI4tZ/85CcxYcKEqKioiIqKimhsbIzHHnss936a/bIiWYH49a9/HbfffnssWbIknn/++bj44otj2rRpsXfv3rSjQZd97GMfi9dffz13/O1vf0s7EpxRa2trXHzxxfHAAw90+P4999wT9913Xzz44IOxYcOG6N+/f0ybNi2OHDnyASeFMztTe46ImD59eru++pe//OUHmBA6Z926dbFgwYJ45pln4oknnohjx47FFVdcEa2trbl7brvttvjjH/8YDz/8cKxbty727NkT1113XYqp4WSdacsREfPnz2/XN99zzz0pJYZTO//88+N73/tebNq0KZ577rn4zGc+E1dffXVs3bo1ItLtlzNJkiQfyCfRrRoaGmLy5Mnx4x//OCIistlsDB06NG655ZZYtGhRyumg85YuXRqrV6+OzZs3px0Fzlomk4lVq1bFNddcExHvjiKrr6+PO+64I+68886IiGhubo7BgwfHihUr4oYbbkgxLZzev7fniHdHkh04cOCkEWaQ7958882ora2NdevWxWWXXRbNzc1RU1MTK1eujNmzZ0dExPbt22P8+PGxfv36+MQnPpFyYujYv7fliHdHkn384x+P5cuXpxsOzkJ1dXV8//vfj9mzZ6faLxtJVgDefvvt2LRpU0ydOjV3raSkJKZOnRrr169PMRmcnZdeeinq6+tj1KhRMWfOnNi1a1fakeCcvPLKK9HU1NSun66srIyGhgb9ND3W2rVro7a2NsaOHRtf/vKX46233ko7EpxRc3NzRLz7y1hExKZNm+LYsWPt+udx48bFsGHD9M/ktX9vyyc89NBDMWjQoLjwwgtj8eLF0dbWlkY86LTjx4/Hr371q2htbY3GxsbU++Xe3f4JdLt9+/bF8ePHY/Dgwe2uDx48OLZv355SKjg7DQ0NsWLFihg7dmy8/vrrcffdd8ell14aL774YgwcODDteHBWmpqaIiI67KdPvAc9yfTp0+O6666LkSNHxs6dO+Ob3/xmzJgxI9avXx+9evVKOx50KJvNxq233hqf+tSn4sILL4yId/vn0tLSqKqqanev/pl81lFbjoi48cYbY/jw4VFfXx8vvPBCfOMb34gdO3bE73//+xTTQse2bNkSjY2NceTIkRgwYECsWrUqLrjggti8eXOq/bIiGZBXZsyYkTufMGFCNDQ0xPDhw+M3v/lNzJs3L8VkAJzw/08Rvuiii2LChAkxevToWLt2bUyZMiXFZHBqCxYsiBdffNFap/R4p2rLN998c+78oosuiiFDhsSUKVNi586dMXr06A86JpzW2LFjY/PmzdHc3By//e1vY+7cubFu3bq0Y5luWQgGDRoUvXr1Omm3hzfeeCPq6upSSgXvj6qqqvjoRz8aL7/8ctpR4Kyd6Iv10xSqUaNGxaBBg/TV5K2FCxfGo48+Gk899VScf/75uet1dXXx9ttvx4EDB9rdr38mX52qLXekoaEhIkLfTF4qLS2NMWPGxMSJE2PZsmVx8cUXx49+9KPU+2VFsgJQWloaEydOjCeffDJ3LZvNxpNPPhmNjY0pJoNzd+jQodi5c2cMGTIk7Shw1kaOHBl1dXXt+umWlpbYsGGDfpqC8K9//SveeustfTV5J0mSWLhwYaxatSrWrFkTI0eObPf+xIkTo0+fPu365x07dsSuXbv0z+SVM7XljpzYCEvfTE+QzWbj6NGjqffLplsWiNtvvz3mzp0bkyZNiksuuSSWL18era2t8fnPfz7taNAld955Z1x11VUxfPjw2LNnTyxZsiR69eoVn/3sZ9OOBqd16NChdt/UvvLKK7F58+aorq6OYcOGxa233hrf+c534iMf+UiMHDky7rrrrqivr2+3YyDki9O15+rq6rj77rtj1qxZUVdXFzt37oyvf/3rMWbMmJg2bVqKqeFkCxYsiJUrV8Yf/vCHGDhwYG49m8rKyujXr19UVlbGvHnz4vbbb4/q6uqoqKiIW265JRobG+1sSV45U1veuXNnrFy5MmbOnBnnnXdevPDCC3HbbbfFZZddFhMmTEg5PbS3ePHimDFjRgwbNiwOHjwYK1eujLVr18Zf/vKX9PvlhIJx//33J8OGDUtKS0uTSy65JHnmmWfSjgRddv311ydDhgxJSktLkw9/+MPJ9ddfn7z88stpx4Izeuqpp5KIOOmYO3dukiRJks1mk7vuuisZPHhwUlZWlkyZMiXZsWNHuqHhFE7Xntva2pIrrrgiqampSfr06ZMMHz48mT9/ftLU1JR2bDhJR+04IpKf//znuXsOHz6cfOUrX0k+9KEPJeXl5cm1116bvP766+mFhg6cqS3v2rUrueyyy5Lq6uqkrKwsGTNmTPK1r30taW5uTjc4dOALX/hCMnz48KS0tDSpqalJpkyZkvz1r3/NvZ9mv5xJkiTp/lIcAAAAAOQva5IBAAAAUPQUyQAAAAAoeopkAAAAABQ9RTIAAAAAip4iGQAAAABFT5EMAAAAgKKnSAYAAABA0VMkAwAAAKDoKZIBAAAAUPQUyQAAAAAoeopkAAAAABS9/wd1Jzqc6t8CvgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1500x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "f, ax = plt.subplots(1, 1, figsize=(15, 6))\n",
    "ax.plot(\n",
    "    np.arange(30),\n",
    "    np.r_[np.repeat(0, 30 - 15), treat_effect_vector],\n",
    "    marker=\"o\",\n",
    "    label=\"True ATE\",\n",
    ")\n",
    "\n",
    "ax.errorbar(\n",
    "    timestamps -.3,\n",
    "    restable0[\"Estimate\"],\n",
    "    yerr=restable0[\"Std. Error\"] * 1.96,\n",
    "    fmt=\".\",\n",
    "    label=\"pointwise (HC1)\",\n",
    ")\n",
    "\n",
    "ax.errorbar(\n",
    "    timestamps -.15,\n",
    "    restable[\"Estimate\"],\n",
    "    yerr=restable[\"Std. Error\"] * 1.96,\n",
    "    fmt=\".\",\n",
    "    label=\"pointwise (CRV1)\",\n",
    ")\n",
    "ax.errorbar(\n",
    "    timestamps + 0.15,\n",
    "    restable[\"Estimate\"],\n",
    "    yerr=[\n",
    "        restable[\"Estimate\"] - uniform_res[\"unif_25\"],\n",
    "        uniform_res[\"unif_975\"] - restable[\"Estimate\"],\n",
    "    ],\n",
    "    fmt=\".\",\n",
    "    label=\"Joint (pre/post)\",\n",
    ")\n",
    "\n",
    "ax.axhline(y=0, color=\"k\", linestyle=\"--\")\n",
    "ax.axvline(x=14.5, color=\"k\", linestyle=\"--\")\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulation study"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def onerep(seed=1, **dgp_params):\n",
    "    np.random.seed(seed)\n",
    "    df = panel_dgp(**dgp_params)\n",
    "    # Estimate model\n",
    "    m2 = pf.feols(\"Y_it ~ i(time_id, ever_treated, ref = 14) | unit_id + time_id\", df)\n",
    "    restable0 = m2.vcov(\"HC1\").tidy()\n",
    "    restable = m2.vcov({\"CRV1\": \"unit_id\"}).tidy()\n",
    "    # Pointwise confidence intervals\n",
    "    pointwise_lower_0 = restable0[\"Estimate\"] - 1.96 * restable0[\"Std. Error\"]\n",
    "    pointwise_upper_0 = restable0[\"Estimate\"] + 1.96 * restable0[\"Std. Error\"]\n",
    "    pointwise_lower = restable[\"Estimate\"] - 1.96 * restable[\"Std. Error\"]\n",
    "    pointwise_upper = restable[\"Estimate\"] + 1.96 * restable[\"Std. Error\"]\n",
    "    # uniform 1: pre-treatment as one chunk, post-treatment as another\n",
    "    kv_pre = [\n",
    "        f\"C(time_id, contr.treatment(base=14))[T.{t}]:ever_treated\" for t in range(14)\n",
    "    ]\n",
    "    kv_post = [\n",
    "        f\"C(time_id, contr.treatment(base=14))[T.{t}]:ever_treated\"\n",
    "        for t in range(15, 30)\n",
    "    ]\n",
    "    uniform_confint_pre = m2.confint(keep=kv_pre, joint=True, exact_match=True)\n",
    "    uniform_confint_post = m2.confint(keep=kv_post, joint=True, exact_match=True)\n",
    "    uniform_res = pd.concat([uniform_confint_pre, uniform_confint_post])\n",
    "    uniform_res.columns = [\"unif_25\", \"unif_975\"]\n",
    "    true_effects = np.r_[np.repeat(0, treat_start_period - 1), treat_effect_vector]\n",
    "    pointwise_coverage_0 = (true_effects >= pointwise_lower_0) & (\n",
    "        true_effects <= pointwise_upper_0\n",
    "    )\n",
    "    pointwise_coverage = (true_effects >= pointwise_lower) & (\n",
    "        true_effects <= pointwise_upper\n",
    "    )\n",
    "    uniform_coverage_1 = (\n",
    "        np.r_[\n",
    "            (\n",
    "                np.repeat(0, treat_start_period - 1)\n",
    "                >= uniform_res.loc[kv_pre, \"unif_25\"]\n",
    "            ),\n",
    "            (treat_effect_vector >= uniform_res.loc[kv_post, \"unif_25\"]),\n",
    "        ]\n",
    "        & np.r_[\n",
    "            (\n",
    "                np.repeat(0, treat_start_period - 1)\n",
    "                <= uniform_res.loc[kv_pre, \"unif_975\"]\n",
    "            ),\n",
    "            (treat_effect_vector <= uniform_res.loc[kv_post, \"unif_975\"]),\n",
    "        ]\n",
    "    )\n",
    "\n",
    "    uniform_coverage_1 = pd.Series(uniform_coverage_1, index=uniform_res.index)\n",
    "    return (\n",
    "        pointwise_coverage_0,\n",
    "        pointwise_coverage,\n",
    "        uniform_coverage_1,\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def simulate_coverage(num_simulations=1000, **dgp_params):\n",
    "    num_periods = dgp_params.get(\"num_periods\", 30)\n",
    "    treat_start_period = dgp_params.get(\"treatment_start\", 15)\n",
    "    treat_effect_vector = dgp_params.get(\n",
    "        \"base_treatment_effect\",\n",
    "        0.3 * np.log(2 * np.arange(1, num_periods - treat_start_period + 1)),\n",
    "    )\n",
    "    treat_effect_vector[8:] = 0\n",
    "\n",
    "    results = Parallel(n_jobs=12)(\n",
    "        delayed(onerep)(**dgp_params) for _ in tqdm(range(num_simulations))\n",
    "    )\n",
    "    (\n",
    "        pointwise_coverage0,\n",
    "        pointwise_coverage,\n",
    "        uniform_coverage,\n",
    "    ) = zip(*results)\n",
    "    return (\n",
    "        np.array(pointwise_coverage0),\n",
    "        np.array(pointwise_coverage),\n",
    "        np.array(uniform_coverage),\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1000/1000 [00:25<00:00, 39.63it/s]\n"
     ]
    }
   ],
   "source": [
    "# Run simulation\n",
    "num_periods = 30\n",
    "treat_start_period = 15\n",
    "treat_effect_vector = 0.2 * np.log(\n",
    "    2 * np.arange(1, num_periods - treat_start_period + 1)\n",
    ")\n",
    "treat_effect_vector[8:] = 0  # switch off effects after a week\n",
    "\n",
    "sigma_i, sigma_t = 2, 2\n",
    "\n",
    "dgp_params = {\n",
    "    \"num_units\": 1_000,\n",
    "    \"num_treated\": 5_00,\n",
    "    \"num_periods\": num_periods,\n",
    "    \"treatment_start\": treat_start_period,\n",
    "    \"hetfx\": False,\n",
    "    \"base_treatment_effect\": treat_effect_vector,\n",
    "    \"sigma_unit\": sigma_i,\n",
    "    \"sigma_time\": sigma_t,\n",
    "    \"ar_coef\": 0.8,\n",
    "}\n",
    "pointwise_coverage0, pointwise_coverage, uniform_coverage = simulate_coverage(\n",
    "    num_simulations=1000, **dgp_params\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(np.float64(0.896551724137931),\n",
       " np.float64(0.896551724137931),\n",
       " np.float64(0.9655172413793104))"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pointwise_coverage0.mean(), pointwise_coverage.mean(), uniform_coverage.mean()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "py311",
   "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.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
No results found