doc/examples/py_double_ml_did_pretest.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Python: Difference-in-Differences Pre-Testing\n",
    "\n",
    "This example illustrates how to use the Difference-in-Differences implmentation `DoubleMLDID` of the [DoubleML](https://docs.doubleml.org/stable/index.html) package can be used to pre-test the parallel trends assumptions.\n",
    "The example is based on the great implmentation of the [did-package](https://cran.r-project.org/web/packages/did/vignettes/did-basics.html) in `R`. \n",
    "You can find further references and a detailed guide on pre-testing with the `did`-package at the [did-package pre-testing documentation](https://cran.r-project.org/web/packages/did/vignettes/pre-testing.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from doubleml import DoubleMLData, DoubleMLDID\n",
    "from lightgbm import LGBMClassifier, LGBMRegressor"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At first, we will generate some data to evaluate the corresponding effect estimates. The DGP is slightly adjusted from [did-package pre-testing documentation](https://cran.r-project.org/web/packages/did/vignettes/pre-testing.html) to consider conditional non-linear parallel trends.\n",
    "\n",
    "We will consider $n_{obs}=1000$ units for $2n_t$ time periods $t\\in\\{-n_t, .., n_t - 1\\}$, where the treatment is applied just before $t=0$.\n",
    "The data generating process will take the following form\n",
    "\n",
    "$$\n",
    "Y_{it} = \\theta_t + \\eta_i + D_{i}\\mu_{it} + g(X_{it}) + \\epsilon_{it}.\n",
    "$$\n",
    "Here $\\theta_t = t$ denotes a time fixed effect, $\\eta_i\\sim \\mathcal{N}(0,1)$ is a unit fixed effect, $X_{it}\\sim\\mathcal{N}(0, I_4)\\in \\mathbb{R}^4$ denote 'pre-treatment' covariates for the time period $t$ and $\\epsilon_{it}\\sim \\mathcal{N}(0,1)$ are independent error terms.\n",
    "Further, $\\mu_{it} \\sim \\mathcal{N}(\\tau,1)$ denote time varying individual effects with $\\tau = (\\underbrace{0,\\dots,0}_{n_t\\ times}, n_t, \\dots, 1)^T$ beeing the average effects (coinciding with the ATTE since $\\mu_{it}$ is independent of $X$.) The treatment indicator $D_{i}$ denotes whether unit $i$ was treated between $t=-1$ and $t=0$ and is generated as\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "f_{ps}(X) &= 0.75\\cdot(-X_1 + 0.5X_2 - 0.25X_3 - 0.1X_4X_3 + \\cos(5X_2))\\\\\n",
    "D_{i} &= \\frac{\\exp(f_{ps}(X_{i0}))}{1+\\exp(f_{ps}(X_{i0}))}\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "Finally, the direct effect of the covariates is given via\n",
    "\n",
    "$$\n",
    "g(X) = \\exp(X_2) - \\sin(5X_3) + 2X_4.\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(42)\n",
    "n_obs = 1000\n",
    "n_time_periods = 5\n",
    "\n",
    "time_periods = (np.arange(2*n_time_periods) - n_time_periods)\n",
    "\n",
    "# fixed effects\n",
    "theta = np.zeros(shape=(n_obs, 2*n_time_periods)) + time_periods \n",
    "eta = np.random.normal(loc=0, scale=1, size=(n_obs,1))\n",
    "\n",
    "# covariates\n",
    "X = np.random.normal(loc=0, scale=1, size=(n_obs, 4, 2*n_time_periods))\n",
    "\n",
    "# treatment effects\n",
    "mu_means = np.concatenate((np.zeros(n_time_periods) , np.arange(n_time_periods, 0, -1)))\n",
    "mu = np.random.normal(loc=0, scale=1, size=(n_obs, 2*n_time_periods)) + mu_means\n",
    "\n",
    "# treatment\n",
    "f_ps = 0.75*(-X[:, 0, 0] + 0.5*X[:, 1, 0] - 0.25*X[:, 2, 0] - 0.1*X[:, 3, 0]*X[:, 2, 0] + np.cos(5*X[:, 1, 0]))\n",
    "ps = (np.exp(f_ps) / (1 + np.exp(f_ps))).reshape(-1,1)\n",
    "u = np.random.uniform(low=0, high=1, size=(n_obs,1))\n",
    "treatment = np.ones(shape=(n_obs, 2*n_time_periods)) * (ps >= u)\n",
    "\n",
    "\n",
    "# outcome\n",
    "g_X = (np.exp(X[:, 1, :]) - np.sin(5*X[:, 2, :]) + 2*X[:, 3, :])\n",
    "epsilon = np.random.normal(loc=0, scale=1, size=(n_obs, 2*n_time_periods))\n",
    "Y = theta + eta + treatment*mu + g_X + epsilon"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate a corresponding dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "              y    d       X_0       X_1       X_2       X_3    t      i\n",
      "0     -2.354188  1.0  1.399355  1.317394  0.522835 -0.257377 -5.0    0.0\n",
      "1     -5.045624  1.0  0.924634  0.197600 -0.573700 -1.668584 -4.0    0.0\n",
      "2      5.014432  1.0  0.059630  2.075261 -0.024355  0.399223 -3.0    0.0\n",
      "3     -0.002290  1.0 -0.646937 -0.689188  2.142270  0.647196 -2.0    0.0\n",
      "4      4.024926  1.0  0.698223  1.735964  1.727543 -0.483186 -1.0    0.0\n",
      "...         ...  ...       ...       ...       ...       ...  ...    ...\n",
      "9995   6.412477  1.0  1.152926 -0.405203 -0.654755  0.590320  0.0  999.0\n",
      "9996   9.881465  1.0 -0.410681  0.698244 -0.132454  0.760915  1.0  999.0\n",
      "9997   1.591080  1.0 -1.082804 -0.533489 -0.344834 -1.482790  2.0  999.0\n",
      "9998  13.709026  1.0  0.619454  1.666307  0.304201  1.264086  3.0  999.0\n",
      "9999   7.230956  1.0  0.336612  1.687854 -1.245720 -0.644799  4.0  999.0\n",
      "\n",
      "[10000 rows x 8 columns]\n"
     ]
    }
   ],
   "source": [
    "# generate dataframe\n",
    "Y_df = np.hstack(Y)\n",
    "treatment_df = np.hstack(treatment)\n",
    "X_df = np.transpose(np.hstack(X))\n",
    "time_df = np.hstack(np.ones(shape=(n_obs, 2*n_time_periods))*time_periods)\n",
    "individual_df = np.hstack(np.ones(shape=(n_obs, 2*n_time_periods))*np.arange(n_obs).reshape(-1,1))\n",
    "\n",
    "data = pd.DataFrame(np.column_stack((Y_df,treatment_df, X_df, time_df, individual_df)),\n",
    "                    columns = ['y'] + ['d'] + [f'X_{i}' for i in range(4)] + ['t'] + ['i'])\n",
    "print(data)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will just rely on [LightGBM](https://lightgbm.readthedocs.io/en/latest/index.html) with a small number of trees."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_estimators = 50\n",
    "ml_g = LGBMRegressor(n_estimators=n_estimators, num_leaves=5, verbose=-1)\n",
    "ml_m = LGBMClassifier(n_estimators=n_estimators, num_leaves=5, verbose=-1)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, construct the corresponding `DoubleMLData` for each time period and fit a `DoubleMLDID` model.\n",
    "Remark thatwe have to include the covariates from both time periods (and they are either pre-treatment or independent of previous outcomes)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.DataFrame(np.nan, \n",
    "                  index=range(2*n_time_periods-1), \n",
    "                  columns=['lower', 'effect', 'upper'])\n",
    "df['time'] = time_periods[1:]\n",
    "df['true_effect'] = mu_means[1:]\n",
    "\n",
    "np.random.seed(42)\n",
    "for t_idx, t in enumerate(time_periods[1:]):\n",
    "    if t <= 0:\n",
    "        t_diff = t-1\n",
    "    else:\n",
    "        # compare to outcome before treatment\n",
    "        t_diff = -1\n",
    "    # outcome as the difference for each model\n",
    "    y_diff = data[data['t'] == t]['y'].values - data[data['t'] == t_diff]['y'].values\n",
    "    covariates = np.column_stack((data[data['t'] == t][[f'X_{i}' for i in range(4)]].values, data[data['t'] == t_diff][[f'X_{i}' for i in range(4)]].values))\n",
    "    dml_data = DoubleMLData.from_arrays(x=covariates,\n",
    "                                        y=y_diff,\n",
    "                                        d=data[data['t'] == t]['d'].values)\n",
    "    dml_did = DoubleMLDID(dml_data,\n",
    "                          ml_g=ml_g,\n",
    "                          ml_m=ml_m)\n",
    "    dml_did.fit()\n",
    "\n",
    "    df.at[t_idx, 'effect'] = dml_did.coef\n",
    "    confint = dml_did.confint(level=0.95)\n",
    "    df.at[t_idx, 'lower'] = confint['2.5 %'].iloc[0]\n",
    "    df.at[t_idx, 'upper'] = confint['97.5 %'].iloc[0]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Take a look at the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0wAAAKACAYAAACvwPjMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACJrklEQVR4nOzdeVzUdeLH8fd3hvsaQEFAUPAOD7zPSu20Wq22a7tMa93OX9udtVtqu6W127G1ZW2Hdpfdt5WVmll540GaGgooiooMIHLNfH9/kCQiBgp853g9H495xHy/35l5+w0G3vP9fD9fwzRNUwAAAACAemxWBwAAAAAAT0VhAgAAAIAGUJgAAAAAoAEUJgAAAABoAIUJAAAAABpAYQIAAACABlCYAAAAAKABAVYHaE1ut1vbt29XZGSkDMOwOg4AAAAAi5imqZKSEiUlJclma/g4kl8Vpu3btyslJcXqGAAAAAA8RG5urpKTkxtc71eFKTIyUlLNTomKirI4DQAAAACrFBcXKyUlpbYjNMSvCtOBYXhRUVEUJgAAAAC/e6oOkz4AAAAAQAMoTAAAAADQAAoTAAAAADTAr85hAgAAaE0ul0tVVVVWxwD8UmBgoOx2+zE/D4UJAACgmZmmqR07dqioqMjqKIBfi46OVkJCwjFdg5XCBAAA0MwOlKX4+HiFhYUd0x9rAJrONE2VlZWpoKBAkpSYmHjUz0VhAgAAaEYul6u2LLVp08bqOIDfCg0NlSQVFBQoPj7+qIfnMekDAABAMzpwzlJYWJjFSQAc+Dk8lnMJKUwAAAAtgGF4gPWa4+eQwgQAAAAADaAwAQAAAEADKEwAAADwKuvXr9fQoUMVEhKivn37NrgMaA4UJgAAAE/kdknZ30pr3q75r9vVYi9lGMYRb1OnTm2x1z7UqFGjDpvhmmuuqd1mypQpCg8P14YNG/TVV181uOxYGYah999/v1meC96LacUBAAA8TdaH0tw7peLtvy2LSpLGPCilj2v2l8vPz6/9+s0339S9996rDRs21C6LiIio/do0TblcLgUEtNyfkZMmTdJ9991XZ9nBsw5u3rxZZ511ljp27HjEZUBz4AgTAACAJ8n6UJozvm5ZkqTi/JrlWR82+0smJCTU3hwOhwzDqL2/fv16RUZG6rPPPtOAAQMUHBysRYsWacKECTrnnHPqPM9NN92kUaNG1d53u92aPn260tLSFBoaqoyMDL399tu/mycsLKxOpoSEBEVFRUmqOeqzfPly3XfffbVHvw63TJJyc3N14YUXKjo6WrGxsTr77LO1ZcuWOq/1wgsvqGfPngoODlZiYqJuuOEGSVJqaqok6dxzz5VhGLX34X8oTAAAAJ7C7ao5siTzMCt/XTZ3cosOz2vI5MmTNWPGDP3000/q06dPox4zffp0vfTSS3r66ae1bt063Xzzzbrsssu0YMGCo86Rn5+vnj176tZbb1V+fr5uu+22wy6rqqrS6aefrsjISH377bf67rvvFBERoTFjxqiyslKSNHPmTF1//fX6y1/+ojVr1ujDDz9Uly5dJElLly6VJM2aNUv5+fm19+F/GJIHAADgKbYurn9kqQ5TKt5Ws13aCa0WS5Luu+8+nXrqqY3evqKiQg888IDmzZunYcOGSZI6deqkRYsW6ZlnntHIkSMbfOxTTz2l5557rs6yZ555RpdeeqkSEhIUEBCgiIgIJSQkSKoZMnjosldeeUVut1vPPfdc7bV4Zs2apejoaM2fP1+nnXaa/vnPf+rWW2/VX//619rXGTRokCQpLi5OkhQdHV37nPBPFCYAAABPUbqzebdrRgMHDmzS9ps2bVJZWVm9klVZWal+/fod8bGXXnqp/va3v9VZ1q5duya9fmZmpjZt2qTIyMg6y8vLy7V582YVFBRo+/btOvnkk5v0vPA/FCYAAABPEdHIUtDY7ZpReHh4nfs2m02mWXfoYFVVVe3XpaWlkqRPPvlE7du3r7NdcHDwEV/L4XDUDo07WqWlpRowYIBeffXVeuvi4uJks3FmChqHwgQAAOApOg6vmQ2vOF+HP4/JqFnfcXhrJ6snLi5Oa9eurbNs1apVCgwMlCSlp6crODhYOTk5Rxx+11L69++vN998U/Hx8bUTRhwqNTVVX331lUaPHn3Y9YGBgXK5Wv98MXgWqjUAAICnsNlrpg6XJBmHrPz1/pgZNdtZ7KSTTtKyZcv00ksvaePGjZoyZUqdAhUZGanbbrtNN998s1588UVt3rxZK1as0BNPPKEXX3zxiM9dVlamHTt21Lnt3bu3SfkuvfRStW3bVmeffba+/fZbZWdna/78+brxxhuVl5cnSZo6daoefvhhPf7449q4cWNtvgMOFKqjeX34DgoTAACAJ0kfJ134khSVWHd5VFLN8ha4DtPROP3003XPPffojjvu0KBBg1RSUqLx48fX2eYf//iH7rnnHk2fPl3HHXecxowZo08++URpaWlHfO5nn31WiYmJdW4XX3xxk/KFhYVp4cKF6tChg/74xz/quOOO01VXXaXy8vLaI05XXHGFHnvsMT311FPq2bOn/vCHP2jjxo21z/Hwww/ryy+/VEpKyu+edwXfZZiHDj71YcXFxXI4HHI6nQ0emgUAADgW5eXlys7OVlpamkJCQo7+idyumtnwSnfWnLPUcbhHHFkCvMmRfh4b2w04hwkA4DfKKquVfu/nkqSs+05XWBC/BuHBbPZWnzocQH0MyQMAAACABlCYAAAAAKABFCYAAAAAaACFCQAAAAAaQGECAAAAgAZQmAAAAACgARQmAAAAD1RWWa3UyZ8odfInKqustjoO4LcoTAAAAADQAAoTAACAB3K5zdqvl2QX1rnvLWbPnq3o6GirYzRJc2V+//331aVLF9ntdt10000NLmtNl19+uR544IGjeuzUqVPVt2/fI26zZcsWGYahVatWHdVrNEVWVpaSk5O1b9++Fn8tChMAAICHmbs2X6c8sqD2/oRZS3X8g19r7tr8FnvNCRMmyDCMercxY8Y06vGpqal67LHH6iy76KKL9PPPP7dA2rpau5gdbj8ZhqE33nijdpurr75a559/vnJzc/WPf/yjwWXHYv78+TIMQ0VFRb+7bWZmpj799FPdeOONR/Vat912m7766qva+xMmTNA555xzVM91OCtXrtQFF1ygdu3aKSQkRF27dtWkSZNqv38OLWPp6ekaOnSoHnnkkWbL0BAKEwAAgAeZuzZf176yQjuLK+os3+Es17WvrGjR0jRmzBjl5+fXub3++utH/XyhoaGKj49vxoSeY9asWfX21YECUVpaqoKCAp1++ulKSkpSZGTkYZe1pieeeEIXXHCBIiIijurxERERatOmTTOnqvHxxx9r6NChqqio0KuvvqqffvpJr7zyihwOh+65554GHzdx4kTNnDlT1dUte44fhQkAAMBDuNympn2UpcMNvjuwbNpHWS02PC84OFgJCQl1bjExMTWvb5qaOnWqOnTooODgYCUlJdUerRg1apS2bt2qm2++ufZoi1T/yM+BYV0vvPCCOnTooIiICF133XVyuVx66KGHlJCQoPj4eN1///11cj3yyCPq3bu3wsPDlZKSouuuu06lpaWSao6yTJw4UU6ns/a1p06dKkmqqKjQbbfdpvbt2ys8PFxDhgzR/Pnz6zz37Nmz1aFDB4WFhencc8/Vnj17GrWvoqOj6+2rkJAQzZ8/v7YMnXTSSTIMo8FlkrRo0SKdcMIJCg0NVUpKim688cY6w8wqKip05513KiUlRcHBwerSpYuef/55bdmyRaNHj5YkxcTEyDAMTZgw4bBZXS6X3n77bY0dO7Z22X//+1/16tWr9v77778vwzD09NNP1y475ZRT9Pe//11S3SF5U6dO1YsvvqgPPvigdp8fvF9/+eUXjR49WmFhYcrIyND333/f4H4sKyvTxIkTdeaZZ+rDDz/UKaecorS0NA0ZMkT//ve/9cwzzzT42FNPPVWFhYVasGBBg9s0BwoTAACAh1iSXah8Z3mD601J+c5yLckubL1Qv3rnnXf06KOP6plnntHGjRv1/vvvq3fv3pKkd999V8nJybrvvvtqj7Y0ZPPmzfrss880d+5cvf7663r++ed11llnKS8vTwsWLNCDDz6ov//97/rxxx9rH2Oz2fT4449r3bp1evHFF/X111/rjjvukCQNHz5cjz32mKKiompf+7bbbpMk3XDDDfr+++/1xhtvaPXq1brgggs0ZswYbdy4UZL0448/6qqrrtINN9ygVatWafTo0frnP/95TPtp+PDh2rBhQ+0+y8/Pb3DZ5s2bNWbMGJ133nlavXq13nzzTS1atEg33HBD7fONHz9er7/+uh5//HH99NNPeuaZZxQREaGUlBS98847kqQNGzYoPz9f//nPfw6bafXq1XI6nRo4cGDtspEjRyorK0u7du2SJC1YsEBt27atLT5VVVX6/vvvNWrUqHrPd9ttt+nCCy+sc0Ry+PDhtev/9re/6bbbbtOqVavUrVs3XXzxxQ0eBfr888+1e/fu2v+fhzrSUMugoCD17dtX3377bYPbNIeAFn12AAAANFpBScNl6Wi2a6qPP/643pCtu+++W3fffbdycnKUkJCgU045RYGBgerQoYMGDx4sSYqNjZXdbldkZKQSEhKO+Bput1svvPCCIiMjlZ6ertGjR2vDhg369NNPZbPZ1L17dz344IP65ptvNGTIEEmqM0FCamqq/vnPf+qaa67RU089paCgIDkcDhmGUee1c3JyNGvWLOXk5CgpKUlSzR/6c+fO1axZs/TAAw/oP//5j8aMGVP7x3q3bt20ePFizZ0793f31cUXXyy73V5nWVZWljp06FA7DDE2NrY20+GWTZ8+XZdeemntv69r1656/PHHNXLkSM2cOVM5OTmaM2eOvvzyS51yyimSpE6dOtW+XmxsbO1zH6lYbN26VXa7vc7wyF69eik2NlYLFizQ+eefr/nz5+vWW2+tLV1LlixRVVVVnSJ0QEREhEJDQ1VRUXHY/9+33XabzjrrLEnStGnT1LNnT23atEk9evSot+2B8nq4dY2RlJSkrVu3HtVjG4vCBAAA4CHiI0OadbumGj16tGbOnFln2YE/yi+44AI99thj6tSpk8aMGaMzzzxTY8eOVUBA0/6cTE1NrXP+Trt27WS322Wz2eosKygoqL0/b948TZ8+XevXr1dxcbGqq6tVXl6usrIyhYWFHfZ11qxZI5fLpW7dutVZXlFRUXsuzk8//aRzzz23zvphw4Y1qjA9+uijtSXmgAPFrLEyMzO1evVqvfrqq7XLTNOU2+1Wdna21qxZI7vdrpEjRzbpeQ+1f/9+BQcH1w6VlGomrjjxxBM1f/58nXLKKcrKytJ1112nhx56SOvXr9eCBQs0aNCgBvfvkfTp06f268TERElSQUHBYUuRaR7b8NLQ0FCVlZUd03P8HgoTAACAhxicFqtER4h2OMsPex6TISnBEaLBabEt8vrh4eHq0qXLYdelpKRow4YNmjdvnr788ktdd911+te//qUFCxYoMDCw0a9x6LaGYRx2mdvtllQzO9of/vAHXXvttbr//vsVGxurRYsW6aqrrlJlZWWDf9CXlpbKbrdr+fLl9Y4EHe3EBwdLSEhocF81Vmlpqa6++urDzlzXoUMHbdq06Zie/4C2bduqrKxMlZWVCgoKql0+atQo/e9//9O3336rfv36KSoqqrZELViw4KiL2sH/Pw+UtAP/Pw91oNCuX79ew4YNa/JrFRYWqnPnzkeRsvE4hwkAAMBD2G2GpoxNl1RTjg524P6Usemy2w5d2zpCQ0M1duxYPf7445o/f76+//57rVmzRlLN+SQul6vZX3P58uVyu916+OGHNXToUHXr1k3bt2+vs83hXrtfv35yuVwqKChQly5d6twODCM77rjj6pwrJUk//PBDs/8bGtK/f39lZWXVy9elSxcFBQWpd+/ecrvdDU5qcKD8/N5+PzBZQ1ZWVp3lB85jeuutt2rPVRo1apTmzZun77777rDnLx382s3x//u0005T27Zt9dBDDx12/e9Nmb527Vr169fvmHMcCYUJAADAg4zplaiZl/VXfFRwneUJjhDNvKy/xvRKbLHXrqio0I4dO+rcdu/eLalmNrnnn39ea9eu1S+//KJXXnlFoaGh6tixo6SaoXYLFy7Utm3bah/THLp06aKqqio98cQT+uWXX/Tyyy/XmcntwGuXlpbqq6++0u7du1VWVqZu3brp0ksv1fjx4/Xuu+8qOztbS5Ys0fTp0/XJJ59Ikm688UbNnTtX//73v7Vx40b997//bdRwPKnmD/lD91VTL6J65513avHixbWTTmzcuFEffPBB7aQPqampuuKKK3TllVfq/fffV3Z2tubPn685c+ZIkjp27CjDMPTxxx9r165dtTMHHiouLk79+/fXokWL6izv06ePYmJi9Nprr9UpTO+//74qKio0YsSIBrOnpqZq9erV2rBhg3bv3q2qqqom/dsPCA8P13PPPadPPvlE48aN07x587RlyxYtW7ZMd9xxh6655poGH7tlyxZt27at3tDI5kZhAgAA8DBjeiVq3i2/DYeaPXGQFt15UouWJUmaO3euEhMT69yOP/54STWzlT377LMaMWKE+vTpo3nz5umjjz6qPR/ovvvu05YtW9S5c2fFxcU1W6aMjAw98sgjevDBB9WrVy+9+uqrmj59ep1thg8frmuuuUYXXXSR4uLiao9WzJo1S+PHj9ett96q7t2765xzztHSpUvVoUMHSdLQoUP17LPP6j//+Y8yMjL0xRdf1E6j/XsmTpxYb1898cQTTfq39enTRwsWLNDPP/+sE044Qf369dO9995b51yomTNn6vzzz9d1112nHj16aNKkSbXFrH379po2bZomT56sdu3a1Zld71B//vOf65wrJdUMlzvhhBNkGEbt/+c+ffooKipKAwcOVHh4eIPPN2nSJHXv3l0DBw5UXFycvvvuuyb92w929tlna/HixQoMDNQll1yiHj166OKLL5bT6TzirIWvv/66TjvttNrS3lIM81jPtPIixcXFcjgccjqdioqKsjoOAKCVlVVWK/3ezyVJWfedrrAgTuVF8ysvL1d2drbS0tIUEnL0kzPw/YrmtH//fnXv3l1vvvnmUZ0r5GkqKyvVtWtXvfbaa0c8Enakn8fGdgN+8gAAADxQWFCAtsw4y+oY8BGhoaF66aWXmnW4pJVycnJ09913H7EsNRcKEwAAAOAHjjSJg7c5MDlGa+AcJgAAAABoAIUJAACgBfjRaeKAx2qOn0MKEwAAQDM6cNHOsrIyi5MAOPBz2JSLKx+Kc5gAAACakd1uV3R0tAoKCiRJYWFhMgxrLjQL+CvTNFVWVqaCggJFR0fLbrcf9XNRmAAAAJpZQkKCJNWWJgDWiI6Orv15PFoUJgAAgGZmGIYSExMVHx+vqqoqq+MAfikwMPCYjiwdQGECAABoIXa7vVn+YANgHSZ9AAD4DVd1de3XS75fWOc+AACH4zWFaerUqTIMo86tR48eVscCAHiJuZ99qFP++U7t/Qmf7dfx976puZ99aGEqAICn86oheT179tS8efNq7wcEeFV8AIBF5n72oa5dYJMpR53lO9xRunaBoZn6UGPOGGdROgCAJ/OqxhEQEHDMs1wAAPyLq7pa074t+bUs1Z3a2ZRNhtya9m2xTj21WnY+iAMAHMJrhuRJ0saNG5WUlKROnTrp0ksvVU5OzhG3r6ioUHFxcZ0bAMC/LPnxW+W7o3VoWTrAlE357hgt+fHbVs0FAPAOXlOYhgwZotmzZ2vu3LmaOXOmsrOzdcIJJ6ikpKTBx0yfPl0Oh6P2lpKS0oqJAQCeoKCwqFm3AwD4F8M0TdPqEEejqKhIHTt21COPPKKrrrrqsNtUVFSooqKi9n5xcbFSUlLkdDoVFRXVWlEBABb6/rtvdPFHZb+73et/CNWw409qhUQAAE9QXFwsh8Pxu93Aa44wHSo6OlrdunXTpk2bGtwmODhYUVFRdW4AAP8yeMgJSrQVyZC7gS1qPjf838/hynfub71gAACv4LWFqbS0VJs3b1ZiYqLVUQAAHsweEKApJ0RKMuqVJkNuGZKu7eXSuvwSnfbIQr2+JEdeOvgCANACvKYw3XbbbVqwYIG2bNmixYsX69xzz5XdbtfFF19sdTQAgIcbc8Y4zRzpVryt7nmvCTanZo50687LxunLW0bqjN4J+tt7a7R+R8PnxwIA/IvXzJ+al5eniy++WHv27FFcXJyOP/54/fDDD4qLi7M6GgDAC4w5Y5w+2L1Un60rkCTNPiNUJ4w4vXYqcUdooB46P0PXjOysTnERcrlNfZS5XWMzkmS3HX6GPQCA7/OawvTGG29YHQEA4OXWbv/tyNHgYSce9rpLneIiJEk/Zu/RTW+u0ss/bNWD5/VRl/iIVssJAPAcXjMkDwCAY1G4r1K5exs/qcPwzm015+phKtxXqTMf/1ZPfrNJ1a6GJo4AAPgqChMAwC+szitq8mMGp8Xqs7+eoIkjUvXwFxs0d92O5g8GAPBoXjMkDwCAY7Emz6mo0AAV769u0uNCAu2664zjdF7/ZHX9dVje5+t2aFT3OAUH2FsiKgDAg3CECQDgF64b3UXvXDv8qB/frV2kDMNQ3t4y3fDaCv3h8UVambO3GRMCADwRhQkA4BfsNkPto0OP+XmSY8L00f8dr9Agu86buVj//DhL+ytdzZAQAOCJKEwAAJ+Xt7dM581crE0Fpc3yfD0SovTutcN1x5geeumHrXp6weZmeV4AgOfhHCYAgM9blVuk5Vv3KiYssNmeM8Bu0zUjO+u09HZqFxUiSfp+8x71TnYoIphfrwDgK3hHBwD4vMzcIrWPDlWbiOBmf+4D120qr3Lp/15fqSC7oQf+2Fujusc3+2sBAFofQ/IAAD4vM8+pvinRLfoaIYF2vXfdcHWOj9CEWUt165xMFZVVtuhrAgBaHoUJAODTql1urclzqk+yo8VfKyU2TC9dOVgPnddHX2Tt0N/eW9virwkAaFkMyQMA+LxnLh+gtLbhrfJahmHowkEpGtk9TlUutyRp3Xan4iNDFBfZ/EMCAQAti8IEAPBpAXabTuwWJ0kqq2zaRWuPxYGJIEzT1N/eW6ste/Zp6tieOrtvkgzDaLUcAIBjw5A8AIBPe+3HHL29PM+y1zcMQ89dMVDHd2mrm95cpateXKZ8537L8rSksspqpU7+RKmTP2nVcgoALYnCBADwaa/+uFVLsvdYmqFtRLD+e0l/PXP5AK3Z5tSEF5bKNE1LMwEAGocheQAAn1Ve5dL6HSW6ZEgHq6NIkk7vmaChaW2Uu7dMhmFoe9F+VbtMdWgTZnU0AEADOMIEAPBZ67Y75XKbykiOtjpKLUdYoHq1r5mx75Evf9bpjy3UC4uy5XJzxAkAPBGFCQDgs1blOhUcYFP3hEiroxzW1HE9dcHAZN33cZYufOZ7bSootToSAOAQFCYAgM/qlRSlm0/tpkC7Z/66iwgO0H1n99Kbfxmqwn2V+uNT36m0gskSAMCTcA4TAMBnDenURkM6tbE6xu8a0qmNPvvrCVqZU6SI4ADtr3Rpa+E+9UiIsjoaAPg9z/zIDQCAY1RcXqX3VubJub/K6iiNEhJo17DONeXupe+36A+PL9IjX2xQRbXL4mQA4N8oTAAAn7Qyp0g3v5mpvfsqrY7SZBNHpOmGk7po5oLNGvvEIq3KLbI6EgD4LQoTAMAnZeYWyREaqI5eOGV3UIBNN53STR/93/EKCbTrj099x4QQAGARzmECAPikzNwi9Ul2yDAMq6MctR4JUXr32uGav2GXusRHyDRN/ZRfovQkzm0CgNbCESYAgM8xTVOZeUXqmxJtdZRjFmC36ZT0dpKkbzYU6MzHv9Xf31/DbHoA0EooTAAAn7O/yqX+HWI01AtmyGuKUd3iNW1cT727YptOf3ShFvy8y+pIAODzKEwAAJ8TFhSg/40fqBFd2lodpVnZbIauGJ6qz286UZ3iwnXFC0s0L2un1bEAwKdxDhMAwOfk7ClTTHigIkMCrY7SIlJiw/TSlYP12dodGtk9TpK0qaBEXeIjLU4GAL6HwgQA8Dm3v52p2PAgzbxsQJ3lYUEB2jLjLItSNS/DMHRm70RJ0i+7SnX6Y99qTM8ETR3XU3GRwRanAwDfwZA8AIBPcblNrdnmVIYPTPjQWGltw/XIhRn6/pc9OvXRBXp/5TaZpml1LADwCRQmAIBP2VRQqrJKlzKSo62O0moMw9DZfdvry5tP1Ald43TTm6v07Le/WB0LAHwCQ/IAAD4lM7dIhiH1TnZYHaXVtYkI1hMX99O4jCT17xAtSdq6Z586xIZ59fWoAMBKHGECAPiUXaUV6pkUpYhg//1M8NT0dmoTEax9FdU696nFuvS5H5Wzp8zqWADglShMAACfcv3oLvrw+uOtjuERwoMD9NhFfbV1T5lOf2yhXliULZebc5sAoCkoTAAAn+F2m3K7TdlsDD874MRucfr85hN1wcBk3fdxlv723hqrIwGAV6EwAQB8xsrcImXc94Wyd++zOopHiQgO0H1n99Kcq4dp/LBUSdLO4nJVu9zWBgMAL+C/A7wBAD4nM7dIFdVuJceEWh3FIw1Oi5UkmaapSS8tk9s09dB5GUpPirI4GQB4Lo4wAQB8RmZekXomRSnQzq+3IzEMQ/84u5eqqk2N++8iPfLFBlVUu6yOBQAeid8oAACfkZlb5FfXXzoWGSnR+uj/jtf1o7voqfmbNf75JVzsFgAOgyF5AACf4Cyr0pY9ZeqbEm11FK8RFGDTzad20xm9E5TvLJdhGHLur1KQ3abQILvV8QDAI1CYAAA+wREWqB/vPllh/KHfZD0SotQjoeY8pvs+ytLyrYWacV4fDe3UxuJkAGA9huQBAHxGu6gQRYYEWh3Dq103urPaRgTrT//7QX9/f41KK6qtjgQAlqIwAQB8wpQP1up/CzdbHcPrdY6L0Jyrh2nauJ56d8U2nfGfhdpfyYQQAPwXQ/IAAF7PNE19smaHLh6cYnUUn2CzGbpieKpO6hGv73/Zo9Agu6pdbu2rcMkRxhE8AP6FI0wAAK+33Vmu3aUV6sMMec0qJTZMFw6sKaEvfJetUx5doLlrd1icCgBaF4UJAOD1MnOLJEkZyQ5rg/iws/u2V0ayQ9e8slzXv7ZCu0srrI4EAK2CwgQA8HqZuUVKcoQoPirE6ig+q11UiJ4dP1D/+VNffb95j059ZIFyC8vqbONy/3YdpyXZhXXuA4C34hwmAIDXu2BgioZ1ZgrslmYYhs7u217Hd2mrt5bnKTkmVJLk3F+l7zfv1pQP19VuO2HWUiU6QjRlbLrG9Eq0KjIAHDPD9KPLehcXF8vhcMjpdCoqKsrqOAAAeL35Gwp0zcvLVV7trrfO+PW/My/rT2kC4HEa2w0YkgcA8Go5e8o047P1KtxXaXUUv9SnfbQMwzjsugOfyE77KIvheQC8FoUJAODVfsjeo2cWblZQAL/SrLBhZ4n2VzV8nSZTUr6zXEuyC1svFAA0I367AAC8WmZukbrGRygimNNyrVBQUt6s2wGAp6EwAQC8WmZeEddfslB8ZONmJmzsdgDgaShMAACvVV7l0vr8EmWkRFsdxW8NTotVoiNEhz+LqUaiI0SD02JbLRMANCcKEwDAa1VUu3XtqM4azpTilrHbDE0Zmy5JDZamG0/uKrvtSJUKADwXhQkA4LUcoYG69bTu6hwXYXUUvzamV6JmXtZf8VHBdZa3iwpWdFigXvlhq/ZVVFuUDgCODYUJAOC1Fvy8S+t3FFsdA6opTfNuGVl7f/bEQVo8+WS9PmmoisqqlL17n4XpAODoUZgAAF5r6ofr9MaSXKtj4FcHD7sbnBYru83QcYlR+ua2UerV3iE312IC4IUoTAAAr1RUVqns3fvUlwkfPF5QgE1VLrf+8vIyvfzDVqvjAECTUJgAAF5pdZ5TktQn2WFxEjRGoN2m5JgwTf1wneZvKLA6DgA0GoUJAOCVMnOLFBUSoNQ24VZHQSPd84d0jewWpxteW8m5ZwC8BoUJAOCVokIDNTYjSTamq/Yadpuhxy/up5TYMF01e5mKyiqtjgQAvyvA6gAAAByNK4anWh0BRyEiOEAvTBiojzPz5QgNtDoOAPwujjABALzOvopqbSvaL9Nk1jVvlOgI1aQTO8kwDK3I2cvseQA8GoUJAOB1Fv68SyNmfK1dpRVWR8ExyC0s04VPf68HP19vdRQAaBCFCQDgdVblFSnJEaL4yBCro+AYpMSGafIZPfTMgl/0xpIcq+MAwGFxDhMAwOtk5hYpg+sv+YSrjk9T9u59+vv7a5USG6YRXdpaHQkA6uAIEwDAq7jcptbkOdUnOdrqKGgGhmFo2rieGt6lrR7+YgPnpQHwOBxhAgB4lR3F5QoKsCkjhQvW+ooAu01PXtJPbndNgQIAT8IRJgCAV2kfHaoV95yqoWltrI6CZhQZEihHWKByC8t021uZKq9yWR0JACRRmAAAXsbtNmUYBhes9VGF+yr18ertuu2tTKYbB+ARKEwAAK9y3tOL9cgXG6yOgRaSkRKtRy/sq49X5+vReT9bHQcAKEwAAO9RXuXSmjyn4qKYTtyXndE7UZPP6KEnvt6kt5fnWR0HgJ9j0gcAgNfIyi9WtdtUX2bI83lXn9hJO5zligi2Wx0FgJ+jMAEAvEZmbpGC7DZ1T4i0OgoOx33QRA1bFktdRki2oys8hmFo6riekiTTNLW3rEqx4UHNkRIAmoQheQAAr5GZW6T0pCgFBfDry+NkfSg9Ofi3+6+eLz3Wq2b5MXp03kad+9R3KtxXeczPBQBN5bW/cWbMmCHDMHTTTTdZHQUA0Eoe+GNvPXFxP6tj4FBZH0pzxkvF+XWXF+fXLD/G0nR+/2SVllfr6peXqaKa6cYBtC6vLExLly7VM888oz59+lgdBQDQisKCApQSG2Z1DBzM7ZLm3inpcFOA/7ps7uS6w/WaqEObMP1v/EBl5jk1+Z01Mk2mGwfQeryuMJWWlurSSy/Vs88+q5iYmCNuW1FRoeLi4jo3AIB3WpJdqCtnL5Vzf5XVUXCwrYul4u1H2MCUirfVbHcMBnSM0cMXZOi9ldv0waojvR4ANC+vK0zXX3+9zjrrLJ1yyim/u+306dPlcDhqbykpKa2QEADQEpZk79GyLYWKDGa+Io9SurN5tzuCsRlJevHKwfpDn8Rjfi4AaCyvKkxvvPGGVqxYoenTpzdq+7vuuktOp7P2lpub28IJAQAtZVWuUxkp0bLZDKuj4GAR7Zp3u98xslucAuw2/fhLTYEGgJbmNR/T5ebm6q9//au+/PJLhYQ07oKFwcHBCg4ObuFkAICWZpqmVuUW6aJByVZHwaE6DpeikupP+FDLqFnfcXizvaRpmvrvN5u0dptT7103Qqltw5vtuQHgUF5zhGn58uUqKChQ//79FRAQoICAAC1YsECPP/64AgIC5HIxaw4A+Kp8Z7l2l1YogwvWeh6bXRrz4K93Dj369+v9MTOO+npMh2MYhp64uJ9iwoJqzmsr47w2AC3HawrTySefrDVr1mjVqlW1t4EDB+rSSy/VqlWrZLdzJXAA8FVRoYF64uJ+GpQaa3UUHE76OOnCl6SohLrLo5JqlqePa/aXjA4L0gsTBqmwrFJXv7JMldXuZn8NAJC8aEheZGSkevXqVWdZeHi42rRpU285AMC3RAQHaGxGktUxcCTp46ROp0tT59Xcv/RtqcuIZj2ydKjUtuH63+UDdfvbmdrhLFeHNkw5D6D5ec0RJgCA//rfws1atHG31THwew4uR6nDW7QsHTA4LVbzbhmpDm3C5HZzfSYAzc+rC9P8+fP12GOPWR0DANCCXG5Tj3+1SWu2Oa2OAg8VaLeptKJaFz7zvT5ezTWaADQvry5MAADf98uuUpVWVCsjxWF1FHiw8CC72seE6pY5mVqRs9fqOAB8CIUJAODRVuUWyTCk3u0pTGiYYRh68Lw+6tPeoUkvLlNuYZnVkQD4CAoTAMCjZeYVqXNchCJDAq2OAg8XEmjX/8YPVHhwgK6cvVTlVVxyBMCx85pZ8gAA/mlwWhv1SIiyOgYaISwoQFtmnGVphtjwIM2aOEjLt+xVSCCXHAFw7ChMAACPNo7pxNFEneMi1DkuQpK0JLtQg1JjZBiHXlQXABqHIXkAAI+1rWi/vl6/k4uS4qiszivShc98r2e//cXqKAC8GIUJAOCxvly3Q9e8vMLqGPBSfZKjdd2ozpr+2XrNXbvD6jgAvBSFCQDgsTLznEpPilJQAL+ucHRuO627zuyVqJveXKnVeUVWxwHghfgNBADwWJm5ReqbEm11DHgxm83Qwxdm6LjEKD35zSar4wDwQkz6AADwSM6yKv2ye59uOKmL1VHg5UIC7XrhikEKDWLWPABNxxEmAIBHKiyr1OC0WI4woVnEhAcpJNCurO3FuvPt1ap2MZEIgMbhCBMAwCOltQ3XnKuHWR0DPmbPvgq9vSJPwYE2TRvXk+nGAfwujjABADxSbmGZqjgKgGZ2Qtc4/ePsXnrp+62a9d0Wq+MA8AIUJgCAxzFNU3+cuVj/mbfR6ijwQZcM6aC/nNhJ//gkS/OydlodB4CHY0geAMDj7Cgu166SCvVJdlgdBT5q8pgecpZVKTKEP4UAHBnvEgAAj5OZWyRJTPiAFmOzGXrw/D6SpGqXWyXl1YoJD7I4FQBPxJA8AIDHWZXrVEJUiOKjQqyOAj9wzwdrddnzP2pfRbXVUQB4IAoTAMDj7HDuV0YKw/HQOsYPS9XWPWX66xsr5XKbVscB4GEM0zT95p2huLhYDodDTqdTUVFRVscBABxBtcutADuf66F1fLOhQFfNXqqJI9J0zx/SrY4DoBU0thvwmwgA4FEOfI5HWUJrGt09XlPH9dTzi7K1aONuq+MA8CD8NgIAeJR3VmzTCQ99rYpql9VR4GfGD0vVS1cO1ogubayOAsCDUJgAAB4lM7dIwQF2BQfYrY4CP3RitzgZhqG5a3do/Y5iq+MA8AAUJgCAR8nMK1JGcrTVMeDHXG5T//1mo66ctVQFJeVWxwFgMQoTAMBjlFe59FN+MTPkwVJ2m6Fnxw+UyzQ16cVl2l/J8FDAn1GYAAAe46f8YlW5TI4wwXKJjlA9f8Ug/byzVDe/uUpuphsH/BaFCQDgMfokR+uLm0/UcYlc+gHW69Xeoccv7qfs3fu0t6zS6jgALBJgdQAAAA6w2wx1axdpdQyg1qnp7TS6e5wC7Da53KbsNsPqSABaGUeYAAAe44bXVuizNflWxwDqCLDbVFBSrrMe/5ZrNAF+iMIEAPAIzv1V+nh1vvZXcYI9PE9sWJDaRYXo2leXa+POEqvjAGhFFCYAgEdYk+eUJGWkRFsbBDiMALtN/72kn9pHh2ri7KXaXVphdSQArYTCBADwCJl5RYoMDlBam3CrowCHFRkSqOcnDFJFtVtXv7ycmfMAP8GkDwAAj7Aqt0h9UhyycVI9PFj76FA9f8VA5TvL+V4F/ASFCQDgEf5yYidVudxWxwB+V5/kaPVJlkzT1Pe/7NHwzm2tjgSgBTEkDwDgEQalxvKHJ7zK1+sLdMmzP+qtZblWRwHQgihMAADLrcjZq5nzN6uaI0zwIif1iNefBqXo7vfW6PvNe6yOA6CFUJgAAJb7MmunXly8RQF2fi3BexiGoX+c00uD02J1zSvLtXlXqdWRALQAfjMBACyXmVukjBSH1TGAJgu02/TUpQMUFxmsFxdvsToOgBbApA8AAEu53abW5Dl17ejOVkcBjoojNFBv/GWoYsKCWvR1yiqrlX7v55KkrPtOV1gQf8YBrYEjTAAAS/2ye59KKqqVkRxtdRTgqLWNCJbdZuiHX/bo3g/WyjS5RhPgKyhMAABLBdoNjR/WUb2TGZIH77erpEIvfb9V//lqo9VRADQTjuUCACzVsU247ju7l9UxgGYxNiNJW/fs07+/+FmpbcJ1Tr/2VkcCcIw4wgQAsNTX63cq37nf6hhAs7l+dBed1z9Zd7y9WkuyC62OA+AYUZgAAJapqHbp6peX68usnVZHAZqNYRia/sfeOqtPoqJCGcwDeDt+igEAlvkpv0RVLlN9mPABPiYowKZHL+orSSqvcqmi2i1HaKC1oQAcFY4wAQAsk5lbpEC7oeMSI62OArSY615doatfXqbKarfVUQAcBQoTAMAymXlFSk+MUnCA3eooQIu5ZmRnrdhapLvfW8N044AXojABACzjCA3UyO7xVscAWtTgtFg9eH5vvb08T0/N32x1HABNxDlMAADLTBnb0+oIQKs4t1+ytuwu078+36CTj4tXj4QoqyMBaCQKEwDAEs6yKhk2KSqEE+HhH246pauGpMVSlgAvw5A8AIAlXvlxq46f8TXndMBvGIah4V3aSpLmLM1VbmGZxYkANAaFCQBgiczcIvVq75BhGFZHAVrV/kqXnpy/SVfOXirn/iqr4wD4HRQmAIAlMvOKlJESbXUMoNWFBtn1/BWDtLO4XNe/ukJVLqYbBzwZhQkA0Op2OMu1s7hCGVywFn6qS3yEnr58gH74ZY/u/WAtQ1MBD0ZhAgC0uuzd+xRktykjxWF1FMAywzu31QN/7K2NO0u1v8pldRwADWCWPABAqxvWuY3WTjtdgXbOX4J/u3Bgis7rnyy7zVC1y60AO59lA56Gn0oAQKszTVNBATYmfAAk2W2Gsnfv06mPLlRmbpHVcQAcgsIEAGhVbrep4x/8Ru+v3GZ1FMBjJDpCFB0WqD+/tEzbivZbHQfAQShMAIBWlb1nn7YV7VfbiGCrowAeIyTQrv9dPlDBATZdNXupSsqZbhzwFBQmAECrOjDkqHcyEz4AB4uLDNasCYO0be9+3Ton0+o4AH7FpA8AgFaVmVukTnHhcoQGWh0F8Dhd20Xq6csHKCiAz7QBT8FPIwCgVa3Kc6ov118CGjSiS1sNSo1VtcutxZt2Wx0H8HsUJgBAq3p2/ADdfGo3q2MAHu+dFXm69PkfNS9rp9VRAL9GYQIAtKr4yBClxIZZHQPweBcMSNFp6e104xsrtXabUy63WbtuSXZhnfsAWo5hmqbf/LQVFxfL4XDI6XQqKirK6jgA4Hc+zNyuRRt36cHz+nANJqAR9le6dNH/vtfWPfsUZLdpV2ll7bpER4imjE3XmF6JFiYEvFdjuwFHmAAArWbBhl1av6OEsgQ0UmiQXZcO7iDn/uo6ZUmSdjjLde0rKzR3bb5F6QD/QGECALSazLwiZTDhA9BoLrepx77aeNh1B4YITfsoi+F5QAuiMAEAWkVxeZU27ypVRkq01VEAr7Eku1D5zvIG15uS8p3lWpJd2HqhAD9DYQIAtIq1eU6ZptQ3hQvWAo1VUNJwWTqa7QA0HYUJANAqOsdH6IFze6tT2wirowBeIz4ypFm3A9B0FCYAQKtoFxWiS4Z0kM3GhA9AYw1Oi1WiI0RH+qlJdIRocFpsq2UC/A2FCQDQKh7/aqPW7yi2OgbgVew2Q1PGpktSg6Vpyth02fkgAmgxFCYAQIvbWVyuR778WVt277M6CuB1xvRK1MzL+is+KrjO8kRHiJ6+rL+GpLXRE19tlJuZ8oAWEWB1AACA78vMLZIkZsgDjtKYXoka0aWtek/9QpI0e+IgndA1TnaboS/W7dAj837W1sIyPXheH442Ac2MI0wAgBaXmVek+MhgJURxYjpwtOxy13492Fhfe/+0ngl67KK+endFnm6Zs0rVLndDTwHgKHhNYZo5c6b69OmjqKgoRUVFadiwYfrss8+sjgUAaITMXKcyUqJlGHzyDRyVrA+lJwf/dv/V86XHetUsl3R23/Z64uL++nh1vm6ekynTZHge0FwaPSTv8ccfb9R2N95441GHOZLk5GTNmDFDXbt2lWmaevHFF3X22Wdr5cqV6tmzZ4u8JgCgeZzQta2SokOtjgF4p6wPpTnjJTOo7vLi/JrlF74kpY/TWX0SFWA3tHdfJR9OAM3IMBv5EURaWtrvP5lh6JdffjnmUI0VGxurf/3rX7rqqqsatX1xcbEcDoecTqeioqJaOB0AAMAxcrtqjiQVb1eZGaz0ilmSpKzgiQozKiQZUlSSdNMayWav89CPMrfrtJ7tFBxgP8wTA2hsN2j0Eabs7OxmCdYcXC6X3nrrLe3bt0/Dhg1rcLuKigpVVFTU3i8uZjpbAGhtmwpKVFxerf4dYqyOAnifrYul4u1H2MCUirfVbJd2Qu3S3MIy3fZWpoYub6NnLh+gkEBKE3C0vOYcJklas2aNIiIiFBwcrGuuuUbvvfee0tPTG9x++vTpcjgctbeUlJRWTAsAkKRXfsjRbXMyrY4BeKfSnUe1XUpsmGZNGKQl2YW66sWl2l/paoFwgH9odGH6+uuvlZ6eftijNE6nUz179tTChQubNdyhunfvrlWrVunHH3/UtddeqyuuuEJZWVkNbn/XXXfJ6XTW3nJzc1s0HwCgvlW5RUwnDhytiHZHvd3wLm01e+Igrcwp0oRZS7SvorqZwwH+odGF6bHHHtOkSZMOO77P4XDo6quv1qOPPtqs4Q4VFBSkLl26aMCAAZo+fboyMjL0n//8p8Htg4ODa2fVO3ADALSeymq3srYXKyPZYXUUwDt1HF5zjpIamsTBkKLa12x3GEM6tdHLVw1WZEiAmAcCODqNLkyZmZkaM2ZMg+tPO+00LV++vFlCNZbb7a5zjhIAwLOs31GsSpdbfTjCBBwdm10a8+Cvdw5tPL/eHzOj3oQPBxvQMVbPXTFIYUEB2ryrVM79VS0SFfBVjZ70YefOnQoMDGz4iQICtGvXrmYJdTh33XWXzjjjDHXo0EElJSV67bXXNH/+fH3++ect9poAgGNTWlGtnklRSk/kCD9w1NLH1Uwd/tk90sF/akUl1ZSl9HGNehq329Q1Ly9XSKBdL181WNFhQb//IACNP8LUvn17rV27tsH1q1evVmJiYrOEOpyCggKNHz9e3bt318knn6ylS5fq888/16mnntpirwkAODbDO7fVJzeewAxdwLFKHyddv+S3+5e+XTOVeCPLkiTZbIb+86d+2la0Xxc/+6P2lDJKB2iMRhemM888U/fcc4/Ky8vrrdu/f7+mTJmiP/zhD80a7mDPP/+8tmzZooqKChUUFGjevHmUJQDwcDuc5Wrk5f4A/J6Dh92lDj/iMLyGpCdF6Y2/DNWukgpd/OwP2lVCaQJ+T6ML09///ncVFhaqW7dueuihh/TBBx/ogw8+0IMPPqju3bursLBQf/vb31oyKwDAi5SUV2nYjK/0/qptVkcBcJBu7SL1xl+GqrLarZzCfVbHATxeo89hateunRYvXqxrr71Wd911V+0nhoZh6PTTT9eTTz6pdu0aOfUlAMDnrdnmlGlKvZKYIQ/wNF3iIzTvlpEKsNtU7XJrb1mV4iKDrY4FeKRGFyZJ6tixoz799FPt3btXmzZtkmma6tq1q2JiuHo7AKCuzFynIoID1CkuwuooAA4jwF4z0GjGZ+v1RdZOvTZpiJJjwixOBXieRg/JO1hMTIwGDRqkzZs3KyiIGVYAAPVl5hapd3uH7DYu/gJ4sgkjUiVJFz3zg3L2lFkbBvBAR1WYDrj66qu1c+fO5soCAPAhWwvL1CeF4XiAp0uOCdObVw9VUIBNFz7zvbJ3c14TcLBjKkzMfAQAaMinNx6vm0/pZnUMAI2Q6AjVm38ZqoiQAL25NNfqOIBHadI5TAAANJZhGFx/CfAi8VEheuea4YoMqfnzsLzKxc8woGM8wvTZZ5+pffv2zZUFAOAjHv3yZ41/YcnvbwjAozjCAmWzGfrhlz0a9a/5WrfdaXUkwHLHVJh69OihefPm6cMPP1R+fn5zZQIAeLkVOXsVZD+mXzEALNQjIVJxkcG65NkftSaP0gT/dtS/zd555x116dJF06ZN05QpU9S5c2fNmjWrObMBALyQ220qM7dIfZnwAfBa0WFBeuXPQ5TWNlyXPPeDVubstToSYJlGF6bS0tI696dNm6YlS5ZoyZIlWrlypd566y397W9/a/aAAADvsmXPPhWXVysjJdrqKACOgSM0UC9fNVg9EiJ15zur5XIz2Rf8U6ML04ABA/TBBx/U3g8ICFBBQUHt/Z07d3JNJgCAMvOKJEl92kdbmgPAsYsMCdTsiYP1/BWDZLcZzJAMv9ToWfI+//xzXX/99Zo9e7aefPJJ/ec//9FFF10kl8ul6upq2Ww2zZ49uwWjAgC8wZieiUq7PkKOsECrowBoBuHBAQoPDlBJeZUmvbRMN4zuquO7trU6FtBqGl2YUlNT9cknn+j111/XyJEjdeONN2rTpk3atGmTXC6XevTooZCQkJbMCgDwAqFBdvVlOB7Q7MKCArRlxlmWvX6g3abgALuufHGp/nf5AI3qHm9ZFqA1NXnSh4svvlhLly5VZmamRo0aJbfbrb59+1KWAACqrHbrzy8u5QRxwAeFBNr1v/EDdGLXtvrLS8v11U87rY4EtIomFaZPP/1UDz/8sJYtW6bnnntODz30kC699FLdfvvt2r9/f0tlBAB4iQ07SjTvpwJxlgPgm4ID7Hrq0gE6qUe8rnlluXILy6yOBLS4RhemW2+9VRMnTtTSpUt19dVX6x//+IdGjhypFStWKCQkRP369dNnn33WklkBAB5uVV6RAmyG0hOjrI4CoIUEBdj0xCX99MzlA5QSG2Z1HKDFGWYjpztp06aNvvjiCw0YMECFhYUaOnSofv7559r1WVlZuvrqq/Xtt9+2WNhjVVxcLIfDIafTqagofpkDQHO77a1MbdhRoo/+73irowBoJS8sylZMeKDO7ZdsdRSgSRrbDRp9hCk8PFzZ2dmSpNzc3HrnLKWnp3t0WQIAtLzM3CJlcMFawG+Ypqmf8ot1y5xMvbUs1+o4QIto9Cx506dP1/jx43XjjTeqrKxML774YkvmAgB4oZtP7ab20aFWxwDQSgzD0IPn9VGA3abb316tarepiwd3sDoW0KwaPSRPkvbs2aNffvlFXbt2VXR0dAvGahkMyQMAAGh+pmlq6ofr9OL3W/XUpf11Zu9EqyMBv6ux3aDRR5ikmvOY2rRpc8zhAAC+Z/6GAhXuq9Qf+3MeA+BvDMPQ1HE9ldo2nIvawuc0+hymgoKCOvdXrVqlK664QiNGjND555+v+fPnN3c2AIAXeXNpruZwDgPgtwzD0MQRaYoKCdSW3fv00vdbrI4ENItGF6bExMTa0rR48WINHjxYW7du1YgRI1RcXKxTTz1VCxcubLGgAADPVjPhQ7TVMQB4gHk/7dS9H6zTE19ttDoKcMwaPSTv4FOdpk6dqssvv1zPP/987bKbbrpJ06ZN01dffdW8CQEAHq+guFzbneXqmxxtdRQAHuDPJ3RSWaVLD3/5s6pcbt18ajcZhmF1LOCoNOkcpgPWrl2r++67r86ySZMmadSoUc2RCQDgZTLznJLEESYAtW48uasC7TY9OHe9JOmW07pbnAg4Ok0qTCUlJQoJCVFISIiCg4PrrAsJCVFZWVmzhgMAeIfY8EBdPLiDEh0hv78xAL9x7ajOCg6wqXN8hNVRgKPWpMLUrVs3STXD85YtW6Z+/frVrlu3bp2SkpKaNx0AwCsM6BirAR1jrY4BwANdeXyaJMntNvXR6u0al5HE8Dx4lUYXpm+++abO/cTEuvPrZ2dn6y9/+UvzpAIAeA2329T8nws0MDVWUSGBVscB4KGWbinUX99YpR+zC/XPs3vJZqM0wTs06cK13o4L1wJA8/tlV6lOeniBXrpysE7sFmd1HAAebM6yXN35zmqd3z9ZM87rIzulCRZqkQvXHpCTk6P8/HzZbDZ16tSJi9kCgB9b/euED32SHRYnAeDpLhyYokC7oVvnZKrabepf5/dRgL3RV7kBLNGk79CnnnpKHTt2VFpamoYPH66hQ4cqPj5exx9/vJYvX95SGQEAHmxVbpHS2oYrOizI6igAvMC5/ZL1nz/1k2GIc5ngFRpdmP7973/r/vvv1+23365nnnlG3bt319SpU/XJJ5+oU6dOOvHEE7Vs2bKWzAoA8ECZeUXK4OgSgCYYm5GkRy7sK7vN0LrtTlVWu62OBDSo0UPynnzyST333HM644wzJEknnniihg8frh07dmjMmDGKiYnR3XffrS+++KLFwgIAPItpmooKCdSQTgzNBtB0JeVVuuTZHzUoNVZPXtpPwQF2qyMB9TR60ofw8HCtW7dOqampkmp+SQYFBSknJ0eJiYnKzMzU8ccfr5KSkpbMe0yY9AEAAMCzfLOhQFe/vFzDO7fR05cNUEggpQmto7HdoNFD8rp166Yvv/yy9v4333yjoKAgJSQkSKq5cC3jUAHAvxSVVarKxVAaAEdvdPd4vXDFIP3wyx79+cVl2l/psjoSUEejh+TddddduuyyyzRv3jyFhITo3Xff1Y033lhbkubPn69evXq1WFAAgOf55yc/afOuUr133QirowDwYsd3bavZEwfrzndWq6CkXB3bhFsdCajVpOswffbZZ3rllVdUUVGh008/XZMmTapdt2fPHkny6CnGGZIHAM3rtEcXaFBqrO4/t7fVUQD4gCqXW4F2m/ZVVMuUFBF8VFfAOSplldVKv/dzSVLWfacrLKj1XhvWaJHrMJ1xxhm1kz4cypOLEgCg+ZVWVGtjQan+fEInq6MA8BGBv16T6aY3V2l3aYVmTxwsR2igxang77hSGADgqKzJc8o0pb4p0VZHAeBj/u+kLvpl1z5d/vyPKiqrtDoO/ByFCQBwVDbvKlVEcIA6x0VYHQWAj+mTHK3XJg1RbmGZLnn2RxXuozTBOhQmAMBRuWxoRy3528my25ghFUDz65nk0Bt/GaaCknJ9sW6H1XHgxzibDQBw1DgpGkBL6p4QqS9uHqnY8CBJUnmVi+s0odVxhAkA0GQFJeUaMeNrLd+61+ooAHzcgbL03so8nfn4t8p37rc4EfxNoz4a/OMf/9joJ3z33XePOgwAwDusznVqW9F+JThCrI4CwE/07xCjiiq3LnrmB702aYiSY8KsjgQ/0agjTA6Ho/YWFRWlr776SsuWLatdv3z5cn311VdyOBwtFhQA4Dky84rUNiJYSRQmAK2kY5twvfGXoTJl6qJnflBuYZnVkeAnGnWEadasWbVf33nnnbrwwgv19NNPy26vGUPqcrl03XXXcTFYAPATq3KL1DfFIcNgwgcArSclNkxv/mWYLnn2B037aJ2eu2KQ1ZHgB5p8tu4LL7ygRYsW1ZYlSbLb7brllls0fPhw/etf/2rWgAAAz2KapjJzizSJC9YCsEBSdKjevHpY7UVuTdPkwxu0qCZP+lBdXa3169fXW75+/Xq53e5mCQUA8GzvXjdcFwxMsToGAD/VLipEseFBynfu13kzF+vnnSVWR4IPa/IRpokTJ+qqq67S5s2bNXjwYEnSjz/+qBkzZmjixInNHhAA4FkMw1CX+EirYwCAguw2lVW69Kf//aBXrhqi9CROD0Hza3Jh+ve//62EhAQ9/PDDys/PlyQlJibq9ttv16233trsAQEAnmXWd9kqq3Tp+tFdrI4CwM+1iQjW65OG6vIXftQlz9WUpl7tmYQMzavJQ/JsNpvuuOMObdu2TUVFRSoqKtK2bdt0xx131DmvCQDgmz7K3M7wFwAeIyY8SK/+eag6tgnXFS8s0b6Kaqsjwccc0yXamRUPAPxLlcuttduL9Yc+SVZHAYBajtBAvXLVYGXmOhUefEx/3gL1NPkI086dO3X55ZcrKSlJAQEBstvtdW4AAN+1YUeJKqvdykiJtjoKANQRGRKo47u2lWmamv7ZT/rhlz1WR4KPaHIFnzBhgnJycnTPPfcoMTGRaRwBwI+syi1SgM1QT06sBuChKl1urd3m1IuLt+j5KwZpRJe2VkeCl2tyYVq0aJG+/fZb9e3btwXiAAA82eC0WP3znF4KCWREAQDPFBxg1/NXDNJfXl6uK2cv1f/GD9TIbnFWx4IXa/KQvJSUFJmm2RJZAAAerlu7SP1pcAerYwDAEYUE2vW/ywdoRJe2mvTiMq3KLbI6ErxYkwvTY489psmTJ2vLli0tEAcA4Kn2VVTrqfmbtMNZbnUUAPhdIYF2PX3ZAN16WjeGEeOYNHlI3kUXXaSysjJ17txZYWFhCgwMrLO+sLCw2cIBADzHmm1OPTR3g07u0U4JjhCr4wDA7woKsOnqkZ0lSSty9mqHs1xn9k60OBW8TZML02OPPdYCMQAAni4zt0hhQXZ1iY+wOgoANNk7y/P0+pIcPXpRX53dt73VceBFmlyYrrjiipbIAQDwcJl5Rerd3iG7jdlRAXif+87upfIqt25+c5WqXKbOH5BsdSR4iWO6sld5ebkqKyvrLONitgDgmzJznfpDH4ayAPBOdpuhf53fR4F2Q7e/nSlJlCY0SpML0759+3TnnXdqzpw52rOn/gXBXC5XswQDAHiOapdbpxwXz9S8ALyazWbogXN7yxEWqOMSI62OAy/R5Fny7rjjDn399deaOXOmgoOD9dxzz2natGlKSkrSSy+91BIZAQAWC7DbNO3sXhrOBSABeDmbzdBdZxynnkkOlVe5NHdtviTJ5f7tsjlLsgvr3Id/M8wmXlSpQ4cOeumllzRq1ChFRUVpxYoV6tKli15++WW9/vrr+vTTT1sq6zErLi6Ww+GQ0+lk6CAANMG67U6FBQUorW241VEAoNnMWZarO95erT/2b6/vNu3WzuKK2nWJjhBNGZuuMb0YiuyrGtsNmnyEqbCwUJ06dZJUc77SgWnEjz/+eC1cuPAo4wIAPNmMz9brnx9nWR0DAJrVBQOSdUavBL27YludsiRJO5zluvaVFbVHoOC/mlyYOnXqpOzsbElSjx49NGfOHEnSRx99pOjo6GYNBwCwnmmaWp3nVEZKtNVRAKBZuU1pVW7RYdcdGII17aMshuf5uSYXpokTJyozs2ZmkcmTJ+vJJ59USEiIbr75Zt1+++3NHhAAYK2te8rk3F9FYQLgc5ZkFyrfWd7gelNSvrNcS7ILWy8UPE6TZ8m7+eaba78+5ZRTtH79ei1fvlxdunRRnz59mjUcAMB6mXlFkqSMZIe1QQCgmRWUNFyWjmY7+KZjug6TJHXs2FEdO3ZsjiwAAA9UUeXWkLRYRYcFWR0FAJpVfGRIs24H33TMhQkA4NsuHJSiCwelWB0DAJrd4LRYJTpCtMNZrobOUkp0hGhwWmyr5oJnafI5TAAA/+Fym3Lur7I6BgC0CLvN0JSx6ZIko4FtBnaMkd3W0Fr4AwoTAKBBP+UXK2PaF1r963lMAOBrxvRK1MzL+is+KrjO8kRHiMb1SVRWfrHKKqstSgdPwJA8AECDMvOKZLcZ6tYu0uooANBixvRK1IhOMep931eSpNlnhOqEESfKHhCgsspqhQUFqLzKpZBAu8VJYYUmH2Gy2+0qKCiot3zPnj2y2/kmAgBfkplbpB4JkfyRAMC3ZX0o+8whtXcHf/0n2R/vLWV9qLCgAJVWVOsPTyzSswt/sTAkrNLkwmSahz8lrqKiQkFBLTeD0vTp0zVo0CBFRkYqPj5e55xzjjZs2NBirwcAkDJzuWAtAB+X9aE0Z7xUnF93eXF+zfKsDxUeZNfpPdvp/k9/0vOLsq3JCcs0ekje448/LkkyDEPPPfecIiIiate5XC4tXLhQPXr0aP6Ev1qwYIGuv/56DRo0SNXV1br77rt12mmnKSsrS+Hh4S32ugDgr8qrXMrdW6arjk+zOgoAtAy3S5p7p3TYOfJMSYY0d7KMHmfpttO6q9pt6h8fZynAZuiK4amtmxWWaXRhevTRRyXVHGF6+umn6wy/CwoKUmpqqp5++unmT/iruXPn1rk/e/ZsxcfHa/ny5TrxxBMP+5iKigpVVFTU3i8uLm6xfADga0IC7Vp172lyNzCyAAC83tbFUvH2I2xgSsXbpK2LZaSdoMljesjlMnXfx1k6sVuc0tryob0/aHRhys6uOfw4evRovfvuu4qJiWmxUI3hdDolSbGxDc+LP336dE2bNq21IgGAzwkKYDJVAD6sdGeTtjMMQ3876ziN65tEWfIjTf5N+M0331heltxut2666SaNGDFCvXr1anC7u+66S06ns/aWm5vbiikBwLvd8Xam7vsoy+oYANByIto1eTvDMNQnOVqmaerhLzZozlL+vvR1TS5M5513nh588MF6yx966CFdcMEFzRLq91x//fVau3at3njjjSNuFxwcrKioqDo3AEDjfLdpjwLsXKwRgA/rOFyKSlLDl601pKj2Ndsdxt6ySt357mq9vTyvxSLCek0uTAsXLtSZZ55Zb/kZZ5yhhQsXNkuoI7nhhhv08ccf65tvvlFycnKLvx4A+KNdJRXaVrRffZIdVkcBgJZjs0tjDhwIOLQ0/Xp/zIya7Q5daxi6b1wv/WlQim5/O1PvrqA0+aomF6bS0tLDTh8eGBjYopMqmKapG264Qe+9956+/vprpaUxaxMAtJTVeUWSpIzkaEtzAECLSx8nXfiSFJVQd3lUUs3y9HENPtRmM3T/Ob114YAU3fZWpr7duKuFw8IKTS5MvXv31ptvvllv+RtvvKH09PRmCXU4119/vV555RW99tprioyM1I4dO7Rjxw7t37+/xV4TAPxVZp5TbcKDlBwTanUUAGh56eOk65f8dv/St6Wb1hyxLB1gsxma/sfe+vtZ6RqU2vBkZPBejZ4l74B77rlHf/zjH7V582addNJJkqSvvvpKr7/+ut56661mD3jAzJkzJUmjRo2qs3zWrFmaMGFCi70uAPijK4Z11OjucTIMzmEC4CcOHnaXOvyww/AafKjN0JW/XrMua3uxcgr3aUyvxOZOCIs0uTCNHTtW77//vh544AG9/fbbCg0NVZ8+fTRv3jyNHDmyJTJKqhmSBwBoHW0igtUmItjqGADgdV5bslWvL8nVk5eI0uQjmlyYJOmss87SWWed1dxZAAAeILewTA9/sUG3j+mh9tEMyQOAppg6tqf2llXphtdW6qlLDZ3WM+H3HwSPdlRXJCwqKtJzzz2nu+++W4WFhZKkFStWaNu2bc0aDgDQ+lbk7NX7q7YrLLDxw1EAADUC7DY9dlFfnZreTte/tkJf/dTIi+PCYzX5CNPq1at1yimnyOFwaMuWLfrzn/+s2NhYvfvuu8rJydFLL73UEjkBAK0kM9epDrFhigmvPyMqAOD3Bdptevzifpr8zhpFh/Fe6u2afITplltu0YQJE7Rx40aFhITULj/zzDNb5TpMAICWlZlXpIyUaKtjAIBXC7Tb9PCFGRrQMUZVLrfW5DmtjoSj1OTCtHTpUl199dX1lrdv3147duxollAAAGtUudxau82pDC5YCwDN5rlvs3X+04u1aONuq6PgKDS5MAUHBx/2ArU///yz4uLimiUUAMAaLrepKWN7anSPeKujAIDPuPL4VA3r3EZXvbhUizdRmrxNkwvTuHHjdN9996mqqkqSZBiGcnJydOedd+q8885r9oAAgNYTEmjXJUM6qHNchNVRAMBnBAfY9fRlAzQ4LVZXvrhUP/yyx+pIaIImF6aHH35YpaWlio+P1/79+zVy5Eh16dJFkZGRuv/++1siIwCglXyyOp8hIwDQAkIC7Xp2/EAN7Birb9YXWB0HTdDkWfIcDoe+/PJLfffdd8rMzFRpaan69++vU045pSXyAQBa0RNfb1S/DjE6vmtbq6MAgM8JCbTr+QkDFWSvOWZRUl6lyJBAi1Ph9zTqCFNsbKx27675xPHKK69USUmJRowYoeuuu0533HEHZQkAfEBZZbV+3lmivilM+AAALSU4wC7DMLR4826d8NA3Wr51r9WR8DsaVZgqKytrJ3p48cUXVV5e3qKhAACtb+22YrlNqU9ytNVRAMDnZSRHq2t8hCa8sESrcousjoMjaNSQvGHDhumcc87RgAEDZJqmbrzxRoWGhh522xdeeKFZAwIAWkdmbpFCA+3qGs+EDwDQ0sKDAzRr4mBd8cISXf78j3rtz0PVm0s6eKRGHWF65ZVXdOaZZ6q0tFSS5HQ6tXfv3sPeAADeqX1MqC4f1lEB9ibPBwQAOAoRwQGaPXGQOsdF6O731sg0Tasj4TAMs4n/Z9LS0rRs2TK1adOmpTK1mOLiYjkcDjmdTkVFRVkdBwAAAB6irLJa6fd+LknKuu90hQU1eW60o1ZcXqXS8molRR9+BBdaRmO7QZMnfRg9erSCgoKaJyUAwCOUlFdp2ZZCVVS7rI4CAJYICwrQlhlnacuMs1q1LElSVEigkqJDtXdfpS5//ket31Hcqq+PI2PSBwCAlm4p1PlPf6+C4gqrowCA3zIMqXBfpS599kf9vLPE6jj4FZM+AAC0Ktep2PAgJccwHAQArBIdFqRXrhqiS577UZc8+4NenzRUXdtFWh3L7zV50gfDMJj0AQB8zOq8ImUkO2QYhtVRAMCvxYQH6dU/D1HbiGBd+tyPKimvsjqS32vUEaZ27dppxowZkmomfXj55Ze9ctIHAEB9pmkqM7dIVwxPtToKAEBS7K+lafHmPYoMCbQ6jt9r8tyx2dnZlCUA8CHO/VVqFxWifh1irI4CAPhVm4hgjc1IkiS9uHiLtuzeZ3Ei/9XownTmmWfK6XTW3p8xY4aKiopq7+/Zs0fp6enNGg4A0PKiw4I096YTNbJbnNVRAACHKKus1ovfb9HFz/6gnD1lVsfxS40uTJ9//rkqKn6bPemBBx5QYWFh7f3q6mpt2LChedMBAFpccXkVF0sEAA8VFhSg1ycNVUigXRc/+4NyCylNra3RhenQX6b8cgUA33DlrKW6853VVscAADSgXVSIXp80VAF2Qxc/+4O2F+23OpJfafI5TAAA31Htcmvtdqe6MW0tAHi0BEdNaRqcFitHKBNBtKZGX8bYMIx6080y/SwAeLefd5aqvMqtjJRoq6MAAH5HUnSoHrmwryRp865ShQcFKMERYm0oP9DowmSapiZMmKDg4GBJUnl5ua655hqFh4dLUp3zmwAA3iEzr0h2m6GeSVFWRwEANJJpmrp1Tqac+6v0xl+Gql0UpaklNXpI3hVXXKH4+Hg5HA45HA5ddtllSkpKqr0fHx+v8ePHt2RWAEAz27CjRN3aRSosqNGfnwEALGYYhh7/Uz+VV7l08bM/qKCk3OpIPs0w/Wj2huLiYjkcDjmdTkVF8WkqAJimqeL91XKEMR4eALzNlt379Kf//aCIkJqZ9OIig62O5FUa2w2Y9AEA/JhhGJQlAPBSqW3D9dqkIap2ubVlDxe2bSkUJgDwU8u2FGrMYwu1w8lQDgDwVp3iIvTlLSM1KDVWLrcp5/4qqyP5HAoTAPiplTlF2rJnn9pGBFkdBQBwDALtNX/S/+PjLP3pfz+oqKzS4kS+hcIEAH4qM69Ivds7FGDnVwEA+IKLB3fQzuJyXfb8j3KWcaSpufBbEgD8VGZekTKSo62OAQBoJt0TIvXqn4do2979uvyFHxme10woTADgh/aUVii3cD8XrAUAH3NcYpRe+fMQ5RSW6ZUftlodxycwrTgA+KHKarfWbCtS57gIRYdxDhMA+JrcwjK1jw6VzWbINE0ZhmF1JI/DtOIAgAYFBdg0oGMsZQkAfFRKbJhsNkMLf96ly59fotKKaqsjeS0KEwD4oUe+2KD3V26zOgYAoIVFhQYqM7dIV85aqn2UpqNCYQIAP2Oapl7+Yauyd3ORQwDwdX1TovXiVYOVlV+sK2cvVVklpampKEwA4GdyC/drb1mV+jLhAwD4hf4dYjR74iCt2ebU3e+usTqO1wmwOgAAoHVl5hVJkvokO6wNAgBoNQNTY/XSlYPVNiLY6ihehyNMAOBnMnOLlBIbqjb80gQAvzIwNVapbcNVUl6l6Z/9pPIql9WRvAJHmADAz5x8XDv1bM+lFQDAX/28s1Szv9uin3eU6OnLByg4wG51JI/GESYA8DPDOrfRuf2SrY4BALDIgI4xeu6Kgfpu8x5d/+oKVVa7rY7k0ShMAOBHcgvL9OqPW5laFgD83Ald4/Ts+IFa+PNuXf/aCrndptWRPBaFCQD8yKJNu3XP+2vFBd8BACO7xemZywdoWKc2stla/hdDWWW1Uid/otTJn3jV9OacwwQAfiQzt0jd2kUqLIi3fwCANLpHvEb/+vX8DQU6vktbBdg5pnIw9gYA+JHMPCfXXwIA1JNbWKZJLy3TzXMyVe3inKaDUZgAwE+UVVbr550lyqAwAQAOkRIbpicu7qdP1+Tr1rcy5eKcploUJgDwE6UV1RrbJ1GDUmOsjgIA8EBjeiXq8T/108er83U7pakWg9gBwE/ER4bosT/1szoGAMCDndUnUW7T1AvfZausslqRIYFWR7IchQkA/MTabU4lOELUNiLY6igAAA82NiNJZ/ZOlN1maFdJhdqEB7XKLHqeiiF5AOAnrn11uWbO32x1DACAF7DbDJVXuXTuU9/p7x+s9evrNFGYAMAP7CmtUG7hfiZ8AAA0WkigXTee3FWvL8nRlA/XyTT9szQxJA8A/MDqbU5JUt/kaGuDAAC8yoUDU2Sapu58Z43sNkNTxqbL8LOrn3OECQD8QGZukWLCApUSG2p1FACAl7loUAc9cG5vvbYkRz/vLLU6TqvjCBMA+IGKareGd2nrd58KAgCaxyVDOmh0jzglOkJrh+b5y+8UChMA+IE7x/SwOgIAwMsdKEt/f3+tIkMCdeeY7n5RmhiSBwA+rqLapcpqt9UxAAA+wDAMdY6L0NMLNuvfX2zwi4kgKEwA4OO+WLdTvaZ+Luf+KqujAAB8wJXHp+nvZx2nJ7/ZrEfnbbQ6TotjSB4A+LjM3CLFRwbLEcrV2gEAzePPJ3RStdvUjM/Wq0dCpM7snWh1pBZDYQIAH7c6z8n1lwAAze6akZ2V1jZcJ/eItzpKi2JIHgD4sGqXW2u2Obn+EgCgRZzeM0EBdpuWZBfquW9/sTpOi+AIEwD4sC179ml/lYsjTACAFrUke4/+/cXPcpum/nJiZ6vjNCsKEwD4sC7xkcqccppCA+1WRwEA+LDrR3dReZVbD3y6XnabTVcdn2Z1pGZDYQIAH8dkDwCAlmYYhm49rZuq3ab+8XGWQgPtumRIB6tjNQsKEwD4sAmzlujkHvG6fFiq1VEAAD7OMAzdOaa7woLsGtAxxuo4zYZJHwDAR5VVVuvbjbtlt/FWDwBoHYZh6MaTu6p7QqT2V7r0ZdZOqyMdM36LAoCPWre9WC63qYwUh9VRAAB+6O3luZr00jK9sSRHkuRym7XrlmQX1rnvyRiSBwA+KjO3SCGBNnVrF2l1FACAH7psaEf9vLNUd723Ruvyi/XFuh216ybMWqpER4imjE3XmF6efdFbjjABgI/KzHOqV5JDgXbe6gEArc8wDE0b11PHd2mrl7/fqp3FFXXW73CW69pXVmju2nyLEjYOR5gAwEfdemo3FZdXWR0DAODHTEkbd5Y2uM6QNO2jLJ2aniC7zWjNaI3Gx44A4KNS24arT3K01TEAAH5sSXahdhSXN7jelJTvLNeS7MLWC9VEFCYA8EHLt+7VlA/Wan+ly+ooAAA/VlDScFk6mu2sQGECAB/07cZd+iBzu0ICeZsHAFgnPjKkWbezglf9Jl24cKHGjh2rpKQkGYah999/3+pIAOCRMnOLlJEcLcPwzPHgAAD/MDgtVomOEBk6/BTihkwlOkI0OC22lZM1nlcVpn379ikjI0NPPvmk1VEAwGOZpqnVeU5lpERbHQUA4OfsNkNT+u6TJBly11l34P6Uvvs8dsIHyctmyTvjjDN0xhlnWB0DADxa3t792rOvUn25YC0AwGpul8b8dKdmBrbXlKortFO/HUlKUKGmBL6sMT9tl04/S7LZLQzaMK8qTE1VUVGhiorf5nsvLi62MA0AtI6QQLvuHNND/VJirI4CAPB3WxdLxds1xr5dI4w16l35giRpduAMnWBbI7thSsW/bpd2grVZG+BVQ/Kaavr06XI4HLW3lJQUqyMBQIuLiwzWtaM6KyY8yOooAAB/V7qz9ku78dt5TINtG+rcP3g7T+PThemuu+6S0+msveXm5lodCQBa3DvL8/RTPkfUAQAeIKJd825nAZ8uTMHBwYqKiqpzAwBfVu1y6+/vr9W3G3dZHQUAAKnjcCkqSVJDkzoYUlT7mu08lE8XJgDwNxsLSrW/yqWM5GirowAAUDORw5gHf71zaGn69f6YGR474YPkZYWptLRUq1at0qpVqyRJ2dnZWrVqlXJycqwNBgAeIjO3SDZD6tWeGfIAAB4ifZx04UtSVELd5VFJNcvTx1mTq5G8apa8ZcuWafTo0bX3b7nlFknSFVdcodmzZ1uUCgA8R2aeU93aRSo82Kve3gEAvi59nNTpdGnqvJr7l74tdRnh0UeWDvCq36ijRo2SaR7+KsEAAKlT23ClxIZaHQMAgPoOLkepw72iLEleVpgAAEc26cROVkcAAMCneNU5TACAhu0qqdDGnSVyuzkSDwBAc6EwAYCP+GDVNv3hiUVyMXQZAIBmQ2ECAB+RmedUr/YOBdp5awcAoLnwWxUAfERmbhHXXwIAoJlRmADAB+zdV6mcwjJlpHD9JQAAmhOFCQB8wLai/UpyhHCECQCAZsa04gDgA3q1d2jxXSdbHQMAAJ/DESYA8AHlVS6rIwAA4JMoTADg5UzT1PEPfq3nvv3F6igAAPgcChMAeLm8vfu1u7RSaW3DrY4CAIDPoTABgJfLzCuSJPVhwgcAAJodhQkAvFxmbpHaR4cqLjLY6igAAPgcChMAeLms/GL1TYm2OgYAAD6JacUBwMvNnjhYJeXVVscAAMAnUZgAwMsF2m2KDQ+yOgYAAEcUFhSgLTPOsjpGkzEkDwC82Psrt+mSZ3+Qy21aHQUAAJ9EYQIAL/ZjdqH2lFbKbjOsjgIAgE+iMAGAF8vMLVJGisPqGAAA+CwKEwB4qf2VLm3YWaIMZsgDAKDFUJgAwEtl5TvlcpvK4IK1AAC0GAoTAHipbu0i9cKEgeqeEGl1FAAAfBbTigOAl4oMCdRJPdpZHQMAAJ/GESYA8FJTP1ynH37ZY3UMAAB8GoUJALzQ3n2Vmr14i3YWl1sdBQAAn0ZhAgAvlJlXJElM+AAAQAujMAGAF8rMdcoRGqiObcKsjgIAgE+jMAGAF1qdV6SMlGgZhmF1FAAAfBqz5AGAFxrXN0mhgXarYwAA4PMoTADghc7u297qCAAA+AWG5AGAl8naXqyPV2+XaZpWRwEAwOdRmADAy7y/apumf7qe85cAAGgFFCYA8DKZuUXKSHFYHQMAAL9AYQIAL+Jym1qzzcn1lwAAaCUUJgDwIpsKSlVW6VJGSrTVUQAA8AsUJgDwItVut0Z3j1Ov9gzJAwCgNTCtOAB4kZ5JDs2aONjqGAAA+A2OMAGAF1m7zan9lS6rYwAA4DcoTPA5ZZXVSp38iVInf6Kyymqr4wDNprzKpbOf/E7vrMizOgoAAH6DwgQAXmLddqdcblN9mfABAIBWQ2ECAC+RmetUUIBN3RMirY4CAIDfoDABgJfIzCtSz6QoBdp56wYAoLXwWxcAPMiRzsHbX+nSgA4xFiUDAMA/Ma04AHiJ/40fKNM0rY4BAIBf4QgTAHiBimqXTNOUYRhWRwEAwK9QmADACzyz4Bed/PACq2MAAOB3KEwA4AUyc4vUPibU6hgAAPgdChMAeDjTNJWZ5+T6SwAAWIDCBAAebruzXLtLK9QnOdrqKAAA+B0KEwB4uDV5TklSRrLD4iQAAPgfChMAeLjTe7bTd5NPUnxUiNVRAADwOxQmAPBwhmGofTQTPgAAYAUKEwB4MJfb1DlPfqdvNhRYHQUAAL9EYQIAD7Z5V6lW5RYpJMBudRQAAPwShQkAPNiq3CIZhtSbCR8AALAEhQkAPFhmbpG6xEUoIjjA6igAAPglChMAeLDVeU5lcMFaAAAsw0eWAODB/nFOLwUH8NkWAABWoTABgAfry9ElAAAsxceWAOChvl5foP9+vdHqGAAA+DUKEwB4qLlrd2jeT1x/CQAAK1GYAMBDrdnmZEgeAAAWozABgIfauqdMGSlcfwkAACtRmAA0SVlltVInf6LUyZ+orLLa6jg+LyM52uoIAAD4NQoTAHgQl9us/Xpsn0SlxIRZmAYAAFCYAMBDzF2br1MeWVB7/6PV+TrxX99o7tp8C1MBAODfKEzwOQd/Qr8ku7DOfcBTzV2br2tfWaGdxRV1lu9wluvaV1ZQmgAAsAiFyQKcA9JyDv2EfsKspTr+wa/5YxMezeU2Ne2jLB2u2h9YNu2jLMo/AAAWoDDBZ/AJPbzVkuxC5TvLG1xvSsp3lmtJdmHrhQIAAJIoTPARfEIPb1VWWa31+cWN2ragpOFSBQAAWkaA1QGA5tCUT+iHdW7TesGAQzj3V+mHX/Zo2ZZCLdmyV+u2OdW1XUSjHhsfGdLC6QAAwKEoTPAJO5z7G7XdW8tzFR0WqOMSo1o4ESCZpqm8vfu1dEuhEh2hGta5jZZmF+rql5cryRGiQWmxOn9AsgZ1jNHE2Uu1w1l+2KOkhqQER4gGp8W29j8BAAC/R2GC1/tsTb6mfZTVqG0/XLVddsPQvy7I0M7icj29YLP6pkSrf4cYJceEyjCMFk4Lf7Bo4269sTRHS7cU1p5TN+mENA3r3EYjurTVd5NPUvvo0DqPmTI2Xde+skKGVKc0GQett9v4/gQAoLVRmKzgdv329ZbFUpcRks1uXR4vU1xepY8z89UuKlgnH9dO7Rwh+kNGoj5bu0OFpZVH/IT+q1tHqrzKLalmiN7X6ws067stkqQ24UEa2T1Oj1zYV1LNuSVhQfyI1MP3b63yKpdW5zm1dEuhlm4p1OVDO+rk49op37lf24v265y+7TUoNVYDOsYoJjxIkhQaZFf7oNB6zzWmV6JmXtZfUz5cV2fikgRHiKaMTdeYXomt9u8CAAC/4a/B1pb1ofTZPZIeqLn/6vmSo4005kEpfZyl0TyZ221q8eY9emt5ruau3aEql1uTTuykk49rp/4dYtS/Q4yO79JW176y/NdP6H/7JN74tUJNGZuusKAAhdX83aq+KdFacPtoFe6rVGZukVbm7NWBOSHKq1zqO+1LdWwTpr4p0erXIUZ9U6LVPSHSvz/l9/PvX+f+KoUF2RVot+mfH2fppe+3qtLlVkRwgPp3jFFQQM08OhcMTNEFA1Oa/PxjbEs1Iuhe9db9kqTZgTN0QtAe2W0zJPn+/gUAwBN5XWF68skn9a9//Us7duxQRkaGnnjiCQ0ePNjqWI2T9aE0Z7xkBtVdXpxfs/zCl/zij86mqHa5FWC36bO1O3T9ayvUKS5cfz2lq/7YL1kJjronwI+xLdXMwMc1peoK7dRv53okaI+mBL6sMbYbdbg/OmPDgzS6R7xG94ivXWaa0v3n9tKq3CKtzCnSuyu3yTRNrZ12usKCAvTWslxFhgSqX4dotYvykxPx/fD7N9+5X0u37NXS7JojSBt2lmjO1cM0KDVWA1NrhnEOTI3VcYlRx16kf92/9oP272DbBtlLKn12/wIA4A28qjC9+eabuuWWW/T0009ryJAheuyxx3T66adrw4YNio+P//0nsJLbJc29U2pw4mtDmjtZ6nGW3w5vOmBfRbU+WZOvt5fnKSEqRI9f3E8nHxevd64drv4dog9/ntGv+3eMfbtGGGvUu/IFSb9+Qm9bI7uhJu3f0CB7naMEZZXV2lRQWjtE78Xvt2jttpqpoBMdIeqbEq1bT+uuLvERMk3T986F8oPvX9M0tXlXqVbkFOmCAckyDENXzl6mn/KL1altuAalxuqq49PUqW24JDXvEDk/2L8AAHgrrypMjzzyiCZNmqSJEydKkp5++ml98skneuGFFzR58uR621dUVKii4rdzAYqLf73WyapVUsRB0/jGxEhpaVJ5uZR1mMkD+vev+e+GDdK+fXXXpaZKsbHSrl1Sbm7ddZGRUteuksslfTpL2lCz3jCr1LNykzbEpdZsV+iWKkwpP0f6ZJbUvr/Uvr3Urp20d6+UnV33eUNDpeOOq/l65cqawyEHO+64mm22bpX27Km7rl27mucuKZE2bqy7LjBQ6t275us1a6Sqqrrru3at+Tdt2ybt3Fl3XZs2UseO0v790k8/1V1nGFK/fjVf//RTzTYHS0uTYmK0fUO23nz/Ry3atFsV1S5lJEfr5KFdJUkhcmvAnmzpkH+OMjIku12a/2bt/g0wK9WzcpPyo+I0OGaD7OVuqchdd/+Gh0vdu9c8x4oVqic9XQoJqdn3e/cqTFIfSSqQlJioj//vBO3MLdDm71dq/Y4ibfgpVxHty6Vd0ZqyydSKnL063VWg9LgwdUuIVHtHqGw2Q+rWreZ7Ly9PKiio+5pt20odOkhlZdL69XXX2WxS3741X2dl1XyvHqxTJyk6WtqxQ9q+ve666Oia9ZWV0tq19f+tffvWPP/PP0ulpXXXdehQk2vlp/W+f8uCQqVE1Xz/7XDV3b+S1KuXFBQk/fKLVFRU93mTkqSEhJrlv/xSd11ISM3+l2p+Vt3uuut79JDCwqScHGn37rrr4uOl5OSaf8fPP9ddFxAg9elT8/W6dVJFhapcbn2YuU1Z20v0eXm48txBSijbq1H7Oik+MkRP9jAVM6CNYtq3O+g9Yl39fXis7xG/fFtv/0qSks2ad+nCw+xfP3yP0M6dNc99MIdD6ty5JsuaNarnwHvExo01/6aDpaRIcXFSYaG0ZUvddU18j6gjMbHmVlwsbdpUd11wsNSzZ83Xq1dL1dV113vre8Tu3TU/kweLiKj597jdNT/Lh/KC94g6unSRoqKk/Pya28Fa4++IzMz6z9u7d83P5ebNktNZdx3vETV4j6jBe8RvDn6POPT/a0NML1FRUWHa7Xbzvffeq7N8/Pjx5rhx4w77mClTppiq+Xi2zs1Z89bw2+3SS2sesHFj3eUHbgcMHVp/3csv16z773/rrzvttJp1Tudhn7ff/71q7rs3zjS7BdRf//DDNY+dM6f+un79fssUFFR//dq1Neuuuqr+usmTa9Z98039de3b//a87dvXX//NNzXrJk+uv+6qq2rWrV1bf11Q0G/P269fvfVZ/3nONE3TLL5/Rv3Hjh1b87iCgsP/v3E6a9YPz6i37u+nXlOzf88Nqf+4oUN/y3S45924sWbdpZfWXzdlSs26uXPrr+vc2fx09Xbz1jmrzL3hjnrr8z6ZZ379005z/w031n/sddfVPO/y5fXXRUb+ljc9vf76Dz6oWffAA/XXnX9+zbrc3MP/W8vLa9aPHFl/3bPP/vrDdE29dd+n9KrZv3+LPPzz5ubWPPb88+uve+CBmnUffFB/XXr6b//WyMM89/LlNeuuu67+uptvrlm3eHG9de62bc0FGwrMhz9fb+5ql1Jv/dsPPGd++/Mus+Jv99R/3pZ+j1j84mGfd9+tbUxzSpTfv0eYc+bUrHv44frrGvsecdpp9df99781615+uf66FnyPqNW2bf31ixfXrLv55vrrPPk94tln668bObJmXXn54Z/Xw94jzLZtf3vezp3rr587t2bdlCn111n0d4RZUFCzfuzY+ut4j6i58R5Rc+M94rfbQe8RTv3aDQ58HzTAqPn/7Pm2b9+u9u3ba/HixRo2bFjt8jvuuEMLFizQjz/+WO8xhzvClJKSIueCBYqy4gjTxzdJkvabgTq/cqo2xKVqddgkhe3dX3OESZL+8JjPH2EqLy7V95v36KufdmpVXpHM1DR9fO9YGQUFMvPy6g5na+wnQ1+/Jr35lzr7Nz8qTotiblJYeXnNEaaD928rfTJUUlqujTtL1D4mRO2iQvVMvk3Tv92mhOLd6h1Yrh6JURrWuY2Gd27r2Z8MLf9I+t/FdfZvWVCoPkn8u8JULu04ZP9Kln96XJn1k4IC7NpZvF8PfPqTNu4pV1ZcmmLDg3R2wF79/dTOstsMud1mzdE/Kz893jRf+tcf6uxfSXor+Z8KC6j87Sj0wfvXh98j+PSYT485wvQrjjDV4D3iN7xH1Gim94jiLVvkGDlSTqdTUVENX6PTpwvToYqLi+VwOH53p7QIt0t6rJdUnK8yM0jpFbMkSVnBExVmVEgypKgk6aY1Pn2OQr5zv057ZKFKKqo1ODVW5w9M1pm9ExURfIyjQ71k/5qmqZzCstrJJFbmFumELm112+ndtX5Hse5+d03tjHx9U6I959pQHr5/TdPUlj1lNdN7Zxdq2da9So4J1ctXDVFFtUt3v7tWA1NjNCg1Vp3jwj1jnx7Mw/cvAAC+qLHdwGvOYWrbtq3sdrt2HvKJxM6dO5WQkGBRqiaw2WumXp4zXtKhf6z9en/MDJ/7Y2hncbneXbFNa7c79eQl/ZUQFaKbTu2mk3vEK/XXk+ebhZfsX8Mw1LFNuDq2CdfZfdvXW98+Jkyfr9uh5xfVfBrYv0O03r1uhCRp6ZZCHZcYdezl8mh42P6tdrm1fkeJAu02dU+I1JdZO/WXl5fLMKQeCVE6oWtbjejSVpIUHGDXwxdmtEquo+Zh+xcAAPzGawpTUFCQBgwYoK+++krnnHOOJMntduurr77SDTfcYG24xkofVzM18Gf3SLsOWh6VVPPHkI9MGVzlcuuLdTv11vJcLfx5lwLtNo3plaDyKpdCAu266vi0lnlhL9+/PRKi9MTFNUMOdpVUKDO3SOXVrtr7Fzz9vQxD6hYf+eu1oaJ1bv/2Cg5opT+iLd6/P+8s0edrd2jJlkKtzClSaUW1zuufrIcvzNCQtDaaNWGQ+neMkSM0sEVztBgv//4FAMBXec2QPKlmWvErrrhCzzzzjAYPHqzHHntMc+bM0fr169WuXbvffbylQ/IOUlZeofSp8yRJWVdGKqzLCK//5PjAULOObcJVWe3WsOlfKSU2TBcMTNYf+iS16h+xvrh/3e6aKa9rhvHt1cqcIm0v2q+V954mu83Q1A/XKTTIrn4p0erbIVrxkS13bajW2L9791Vq2da9WralUKekt9Og1Fi9/P0WPfT5Bg3sGKNBabEalBqr3u0dCgn07v+3h/LF718AADyRzw3Jk6SLLrpIu3bt0r333qsdO3aob9++mjt3bqPKkkc5+I+f1OFe/cfQ7tIKvb9ym95enqeNBaX64a6TFRcZrHm3jFRMeNDvP0FL8KH9e4DNZqhru0h1bRepCwfVXBuqotolu82QaZoq3FepH9bs0cz5myVJ7aND9fyEgeqREKWC4nJFhQY2X7Fo5v174DMbwzD0wqJsvb4kRxsLak4aTYgK0XGJURqUKl04KEWXDOl47BeI9XQ++P0LAIA386rCJEk33HCD9wzB82GmaeqG11fq87U7ZDMMnZIerzvH9FBMWM2RJMvKkh85MBTPMAw9fnE/maapfGe5VuYUaVXuXiU6QiVJ936wTvN+2qn0pKjaySRGdGmrdlFHdxTK5f7toPSS7EKd0DWuSSXG7Tb1c0GJlmYXasmWmqNIT13aX/06xCgowKaBqbG6bnRnDewYW2fSi1YbeggAAHAQrytMsM76HcV6f+V23XRKV4UE2tUxNkx/P+s4nd23PQXJAxiGoaToUCVFh+qsPom1y286tatGdGmjlTlFWrRxt176fqv+fUGGzh+QrMWbd2tp9l717RCtvsnRcoQdeejk3LX5mvLhbxdunTBrqRIdIZoyNl1jeiUe9jEV1S6t216sfinRMgxD5z29WCtzihRgM9Q72aGxGUmKDqv5/rlsaMdm2BMAAADNh8KEIyoqq9QHq7brreW5WrutWLHhQRqbkaieSQ7dMaaH1fHQCD0SotQjIUqX/zobf1FZpQLsNknS5oJSvfBdtpzzaq6V0SkuXOOHdtSEEWmqdrllSgr8ddu5a/N17SsrdOhJjzuc5br2lRWaeVl/jemVKJfb1MKNu7RsS6GWZu/VqrwiVVa7tejO0UqOCdP/ndRFIYF29UuJUWgQR40AAIBnozChnmqXW3abIcMwNOmlZVqRU6TR3eP1fyd11eju8QoKsFkdEcfgwNEcSbp8WKouG9pRW/aUaWXOXq3KLVJkSM1Rph+zC3XVi0vVu71DfZIdemfFtnplSVLtstvfXq1T0xNkM6Rb3lwlu83QoNRYTR7TQ4NSY5Xw6xDAk3p42TmHAADAr1GYUGtTQaneXp6nd1fk6enLB6h/hxhNG9dLcZHBiosMtjoeWohhGEprG660tuH6Y//k2uWpbcN122ndtTK3SB+s3K6isqojPItUUl6t7zfv1vFd4/T5zScqLiLY8y4QCwAA0EQUJujj1dv1/KJsrcwpkiM0UGf3TVLsr0ch0pOsm34d1mofHao/n9BJkvTBqm366xurfvcxe/ZVSlKLTmsOAADQmihMfsjtNrV48x51axeh+KgQ/byjRI7QQP33kn465bh2PnddGxy7xhYgihIAAPA1FCY/snXPPr29PE/vLM/Tdme5poxN18QRabr51G4MncIRDU6LVaIjRDuc5Yc9j8mQlOAI0eC02NaOBgAA0KIoTH7iv19v1L+/+FkRwQEam5Go8wekqH+HaEmiLOF32W2GpoxN17WvrJAh1SlNB757poxN9/2LygIAAL9DYfJBpmlqSXah3lqepxO7xWlcRpJGdotXckyYTu+ZwFTOOCpjeiVq5mX9NeXDddpZXFG7POF3rsMEAADgzShMPiTfuV9vLcvT28vzlFNYppTYUI3o0kaS1DvZod7JDosTwtuN6ZWoEV3aqvfULyRJsycO0gld4ziyBAAAfBaFycuVV7lUXF6l+MgQLd60RzPnb9aZvRP10Pl9NDg1Vjb+kEUzO7gcDU6LpSwBAACfRmHyQqZpamVukd5alqePM7fr1PR2euSivjqrT6JO75WgiGD+twLeKiwoQFtmnGV1DAAA8Cv+svYyy7fu1R1vZ2rzrn1KcoRowohUnffrxUaZDhwAAABoXhQmC7jcv80xtiS78IjngFRUuzQvq0AV1S79sX+yEh0h6pnk0NRxPTW8c1uGQwEAAAAtiMLUyuauzdeUD9fV3p8wa6kSD5llzDRNrd1WrLeX5+qDzO0qKqvS6T3b6Y/9k5UUHarHL+5nVXwAAADAr1CYWtHctfm69pUV9S78ucNZrmtfWaEnL+mvM/skat32Yo397yLFRQbrokEpumBAsrrER1qSGQAAAPBnFKZW4nKbmvZRVr2yJP12EdC/vrlSp/dKUM+kKL1y1RAN7RSrALutNWMCAAAAOAh/jbeSJdmFyneWH3GbKpepH3/ZI8MwdHzXtpQlAAAAwGIcYWolBSVHLksH7CqtaOEkvo9pmQEAANBcOITRSuIjQ5p1OwAAAAAtj8LUSganxSrREaKGJgE3JCU6QjQ4LbY1YwEAAAA4AgpTK7HbDE0Zmy5J9UrTgftTxqZzXSUAAADAg1CYWtGYXomaeVl/xUcF11me4AjRzMv6116HCQAAAIBnYNKHVjamV6JGdGmr3lO/kCTNnjhIJ3SN48gSAAAA4IE4wmSBg8vR4LRYyhIAAADgoShMAAAAANAAChMAAAAANIDCBAAAAAANoDABAAAAQAMoTAAAAADQAAoTAAAAADSA6zABaJKwoABtmXGW1TEAAABaBUeYAAAAAKABFCYAAAAAaACFCQAAAAAaQGECAAAAgAZQmAAAAACgARQmAAAAAGgAhQkAAAAAGkBhAgAAAIAGUJgAAAAAoAEUJgAAAABoAIUJAAAAABpAYQIAAACABlCYAAAAAKABFCYAAAAAaACFCQAAAAAaQGECAAAAgAZQmAAAAACgARQmAAAAAGgAhQkAAAAAGkBhAgAAAIAGUJgAAAAAoAEUJgAAAABoAIUJAAAAABoQYHUAfxQWFKAtM86yOgYAAACA38ERJgAAAABoAIUJAAAAABpAYQIAAACABlCYAAAAAKABFCYAAAAAaACFCQAAAAAaQGECAAAAgAZQmAAAAACgARQmAAAAAGgAhQkAAAAAGkBhAgAAAIAGUJgAAAAAoAEUJgAAAABoAIUJAAAAABpAYQIAAACABlCYAAAAAKABFCYAAAAAaACFCQAAAAAa4DWF6f7779fw4cMVFham6Ohoq+MAAAAA8ANeU5gqKyt1wQUX6Nprr7U6CgAAAAA/EWB1gMaaNm2aJGn27NnWBgEAAADgN7ymMB2NiooKVVRU1N4vLi62MA0AAAAAb+M1Q/KOxvTp0+VwOGpvKSkpVkcCAAAA4EUsLUyTJ0+WYRhHvK1fv/6on/+uu+6S0+msveXm5jZjegAAAAC+ztIhebfeeqsmTJhwxG06dep01M8fHBys4ODgo348AAAAAP9maWGKi4tTXFxcq72eaZqSOJcJAAAA8HcHOsGBjtAQr5n0IScnR4WFhcrJyZHL5dKqVaskSV26dFFERESjnqOkpESSOJcJAAAAgKSajuBwOBpcb5i/V6k8xIQJE/Tiiy/WW/7NN99o1KhRjXoOt9ut7du3KzIyUoZhNHPCpikuLlZKSopyc3MVFRVlaRZfxP5tWezflsX+bVns35bF/m1Z7N+Wxf5teZ60j03TVElJiZKSkmSzNTy1g9ccYZo9e/YxX4PJZrMpOTm5eQI1k6ioKMu/WXwZ+7dlsX9bFvu3ZbF/Wxb7t2Wxf1sW+7fleco+PtKRpQN8elpxAAAAADgWFCYAAAAAaACFySLBwcGaMmUK0563EPZvy2L/tiz2b8ti/7Ys9m/LYv+2LPZvy/PGfew1kz4AAAAAQGvjCBMAAAAANIDCBAAAAAANoDABAAAAQAMoTAAAAADQAAqTB6moqFDfvn1lGIZWrVpldRyfMW7cOHXo0EEhISFKTEzU5Zdfru3bt1sdyyds2bJFV111ldLS0hQaGqrOnTtrypQpqqystDqaz7j//vs1fPhwhYWFKTo62uo4PuHJJ59UamqqQkJCNGTIEC1ZssTqSD5h4cKFGjt2rJKSkmQYht5//32rI/mU6dOna9CgQYqMjFR8fLzOOeccbdiwwepYPmPmzJnq06dP7cVUhw0bps8++8zqWD5rxowZMgxDN910k9VRGoXC5EHuuOMOJSUlWR3D54wePVpz5szRhg0b9M4772jz5s06//zzrY7lE9avXy+3261nnnlG69at06OPPqqnn35ad999t9XRfEZlZaUuuOACXXvttVZH8QlvvvmmbrnlFk2ZMkUrVqxQRkaGTj/9dBUUFFgdzevt27dPGRkZevLJJ62O4pMWLFig66+/Xj/88IO+/PJLVVVV6bTTTtO+ffusjuYTkpOTNWPGDC1fvlzLli3TSSedpLPPPlvr1q2zOprPWbp0qZ555hn16dPH6iiNZ8IjfPrpp2aPHj3MdevWmZLMlStXWh3JZ33wwQemYRhmZWWl1VF80kMPPWSmpaVZHcPnzJo1y3Q4HFbH8HqDBw82r7/++tr7LpfLTEpKMqdPn25hKt8jyXzvvfesjuHTCgoKTEnmggULrI7is2JiYsznnnvO6hg+paSkxOzatav55ZdfmiNHjjT/+te/Wh2pUTjC5AF27typSZMm6eWXX1ZYWJjVcXxaYWGhXn31VQ0fPlyBgYFWx/FJTqdTsbGxVscA6qmsrNTy5ct1yimn1C6z2Ww65ZRT9P3331uYDGg6p9MpSbzftgCXy6U33nhD+/bt07Bhw6yO41Ouv/56nXXWWXXeh70BhclipmlqwoQJuuaaazRw4ECr4/isO++8U+Hh4WrTpo1ycnL0wQcfWB3JJ23atElPPPGErr76aqujAPXs3r1bLpdL7dq1q7O8Xbt22rFjh0WpgKZzu9266aabNGLECPXq1cvqOD5jzZo1ioiIUHBwsK655hq99957Sk9PtzqWz3jjjTe0YsUKTZ8+3eooTUZhaiGTJ0+WYRhHvK1fv15PPPGESkpKdNddd1kd2as0dv8ecPvtt2vlypX64osvZLfbNX78eJmmaeG/wLM1df9K0rZt/9/e/cdUVf9xHH9eoIuBRqi3ZIVMQRGbDKWJNNcuSePa8kdzVFvKhU1c7mIQ/dEiMf+wtZqE1ZYz2e7VBeamiFYkOCczQ1E0EKdzwgwKJUiiCXNQl9sfrbvvjW6KYQf4vh7b+eN8zudy3+ezC5wX5/M5tGOz2UhPTyc7O9ugyseGuxlfEZE/ORwOLly4wGeffWZ0KeNKbGwsDQ0N1NXVsX79eux2OxcvXjS6rHHh+++/Jzc3l9LSUiZMmGB0OcNm8uiq8Z7o6urixo0b/9hn5syZPP/883z++eeYTCZvu9vtJjAwkJdeeoldu3bd61LHpDsdX7PZPKT9hx9+IDIyktraWt1q92O443vt2jWsViuLFi3C5XIREKC/xfyTu/n8ulwu8vLy6OnpucfVjV8DAwOEhISwb98+Vq5c6W232+309PTozvMIMplMHDhwwGecZWTk5ORw8OBBjh8/zowZM4wuZ1xLTU0lOjqaHTt2GF3KmFdRUcFzzz1HYGCgt83tdmMymQgICKC/v9/n2GgTZHQB45XFYsFisdy234cffsiWLVu8+9euXSMtLY29e/eSlJR0L0sc0+50fP/O4OAg8Mdj3OXvDWd829vbSUlJITExEafTqbB0B/7N51funtlsJjExkaNHj3ov5AcHBzl69Cg5OTnGFidyGx6Phw0bNnDgwAFqamoUlv4Dg4ODulYYIUuWLKGpqcmnLSsrizlz5vD666+P6rAECkyGmz59us/+xIkTAYiOjubRRx81oqRxpa6ujjNnzrB48WLCw8NpaWmhsLCQ6Oho3V0aAe3t7VitVqKioti6dStdXV3eY9OmTTOwsvGjra2N7u5u2tracLvd3v/RFhMT4/15IXcuPz8fu93O448/zsKFC9m2bRt9fX1kZWUZXdqY19vbS3Nzs3f/6tWrNDQ0MHny5CG/62T4HA4HZWVlHDx4kEmTJnnX3YWFhXH//fcbXN3Y98Ybb7B06VKmT5/OzZs3KSsro6amhqqqKqNLGxcmTZo0ZL3dn2vLx8I6PAUmGddCQkIoLy/nrbfeoq+vj4iICGw2Gxs3biQ4ONjo8sa8I0eO0NzcTHNz85CAr9m+I2PTpk0+U3Pnz58PwLFjx7BarQZVNXa98MILdHV1sWnTJjo6OkhISODw4cNDHgQhw1dfX09KSop3Pz8/H/hjyqPL5TKoqvFj+/btAEO+751OJ5mZmf99QeNMZ2cnGRkZXL9+nbCwMOLj46mqquLpp582ujQZBbSGSURERERExA8tNhAREREREfFDgUlERERERMQPBSYRERERERE/FJhERERERET8UGASERERERHxQ4FJRERERETEDwUmERERERERPxSYRERERERE/FBgEhERERER8UOBSUREDJOZmYnJZBqy2Ww2o0sTEREBIMjoAkRE5P+bzWbD6XT6tAUHBxtUzb83MDCA2Ww2ugwRERkhusMkIiKGCg4OZtq0aT5beHg4NTU1mM1mvv76a2/f9957j4ceeogff/wRAKvVSk5ODjk5OYSFhTF16lQKCwvxeDze1/z8889kZGQQHh5OSEgIS5cu5cqVK97jra2tLFu2jPDwcEJDQ3nssceorKwEwOVy8eCDD/rUW1FRgclk8u5v3ryZhIQESkpKmDFjBhMmTACgp6eHtWvXYrFYeOCBB3jqqadobGwc8fETEZF7S4FJRERGJavVSl5eHmvWrOGXX37h22+/pbCwkJKSEh5++GFvv127dhEUFMTp06f54IMPeP/99ykpKfEez8zMpL6+nkOHDnHy5Ek8Hg/PPPMMv/76KwAOh4P+/n6OHz9OU1MT7777LhMnThxWrc3Nzezfv5/y8nIaGhoASE9Pp7Ozk6+++oqzZ8+yYMEClixZQnd3978fHBER+c9oSp6IiBjqiy++GBJQCgoKKCgoYMuWLRw5coR169Zx4cIF7HY7y5cv9+kbGRlJcXExJpOJ2NhYmpqaKC4uJjs7mytXrnDo0CG++eYbnnjiCQBKS0uJjIykoqKC9PR02traWLVqFfPmzQNg5syZwz6HgYEBdu/ejcViAeDEiROcPn2azs5O7/TCrVu3UlFRwb59+1i3bt2w30NERIyhwCQiIoZKSUlh+/btPm2TJ08GwGw2U1paSnx8PFFRURQXFw95/aJFi3ymyCUnJ1NUVITb7ebSpUsEBQWRlJTkPT5lyhRiY2O5dOkSAK+88grr16+nurqa1NRUVq1aRXx8/LDOISoqyhuWABobG+nt7WXKlCk+/W7dukVLS8uwvraIiBhLgUlERAwVGhpKTEyM3+O1tbUAdHd3093dTWho6Ii+/9q1a0lLS+PLL7+kurqad955h6KiIjZs2EBAQIDPeijAO5Xvr+fwv3p7e4mIiKCmpmZI37+uiRIRkdFNa5hERGTUamlp4dVXX2Xnzp0kJSVht9sZHBz06VNXV+ezf+rUKWbNmkVgYCBxcXH89ttvPn1u3LjB5cuXmTt3rrctMjKSl19+mfLycl577TV27twJgMVi4ebNm/T19Xn7/rlG6Z8sWLCAjo4OgoKCiImJ8dmmTp16N0MhIiIGUWASERFD9ff309HR4bP99NNPuN1uVq9eTVpaGllZWTidTs6fP09RUZHP69va2sjPz+fy5cvs2bOHjz76iNzcXABmzZrFihUryM7O5sSJEzQ2NrJ69WoeeeQRVqxYAUBeXh5VVVVcvXqVc+fOcezYMeLi4gBISkoiJCSEgoICWlpaKCsrw+Vy3facUlNTSU5OZuXKlVRXV/Pdd99RW1vLm2++SX19/cgOoIiI3FMKTCIiYqjDhw8TERHhsy1evJi3336b1tZWduzYAUBERASffPIJGzdu9Hk8d0ZGBrdu3WLhwoU4HA5yc3N9HqrgdDpJTEzk2WefJTk5GY/HQ2VlJffddx8Abrcbh8NBXFwcNpuN2bNn8/HHHwN/rKX69NNPqaysZN68eezZs4fNmzff9pxMJhOVlZU8+eSTZGVlMXv2bF588UVaW1t9nvAnIiKjn8nz18nZIiIiY4TVaiUhIYFt27YZXYqIiIxTusMkIiIiIiLihwKTiIiIiIiIH5qSJyIiIiIi4ofuMImIiIiIiPihwCQiIiIiIuKHApOIiIiIiIgfCkwiIiIiIiJ+KDCJiIiIiIj4ocAkIiIiIiLihwKTiIiIiIiIHwpMIiIiIiIifvwOMACOPsUs4+gAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1000x750 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['figure.figsize'] = 10., 7.5\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "errors = np.full((2, 2*n_time_periods - 1), np.nan)\n",
    "errors[0, :] = df['effect'] -   df['lower']\n",
    "errors[1, :] = df['upper'] - df['effect']\n",
    "\n",
    "plt.errorbar(df['time'], df['effect'], fmt='o', yerr=errors, color='#1F77B4',\n",
    "             ecolor='#1F77B4', label='Estimated Effect (with CI)')\n",
    "ax.plot(time_periods[1:], df['effect'], linestyle='--', color='#1F77B4', linewidth=1)\n",
    "\n",
    "# add horizontal line\n",
    "ax.axhline(y=0, color='r', linestyle='--', linewidth=1)\n",
    "\n",
    "# add true effect\n",
    "ax.scatter(x=df['time'], y=df['true_effect'], c='#FF7F0E', label='True Effect')\n",
    "\n",
    "plt.xlabel('Exposure')\n",
    "plt.legend()\n",
    "_ = plt.ylabel('Effect and 95%-CI')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "dml_dev",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.2"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}