src/plot_treatment_model.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import result_processing\n",
    "import os\n",
    "import argparse\n",
    "from reddit.data_cleaning.reddit_posts import load_term_counts\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from sklearn.linear_model import LogisticRegression, LinearRegression\n",
    "from scipy.special import logit\n",
    "from scipy.stats import chi2, ttest_1samp\n",
    "from result_processing.helpers import convert_str_columns_to_float, assign_split, filter_imbalanced_terms\n",
    "from importlib import reload\n",
    "from sklearn.metrics import mean_squared_error as mse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fit_treatment(features, labels, verbose=False, coeff_offset=1):\n",
    "    model = LogisticRegression(solver='liblinear')\n",
    "    model.fit(features, labels)\n",
    "    coeffs = np.array(model.coef_).flatten()[coeff_offset:]\n",
    "    if verbose:\n",
    "        print(\"Model accuracy:\", model.score(features, labels))\n",
    "        print(\"Mean and std. of the word coeffs:\", coeffs.mean(), coeffs.std())     \n",
    "    return coeffs\n",
    "\n",
    "def truncate(df, truncate_level=0.1):\n",
    "    df = df[(df.treatment_probability >= truncate_level) & (df.treatment_probability <= 1.0-truncate_level)]\n",
    "    return df\n",
    "\n",
    "def fit_outcome(features, labels, verbose=False, coeff_offset=1):\n",
    "    model = LinearRegression()\n",
    "    model.fit(features, labels)\n",
    "    coeffs = np.array(model.coef_).flatten()[coeff_offset:]\n",
    "    if verbose:\n",
    "        print(\"Model MSE:\", mse(labels,model.predict(features)))\n",
    "        print(\"Mean and std. of the word coeffs:\", coeffs.mean(), coeffs.std())\n",
    "    return coeffs "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "log_dir = '../logdir/modesimple/beta01.0.beta110.0.gamma1.0'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(93916, 7)\n"
     ]
    }
   ],
   "source": [
    "predictions_file = os.path.join(log_dir, 'predict', 'test_results_all.tsv')\n",
    "predict_df = pd.read_csv(predictions_file, delimiter='\\t')\n",
    "predict_df = convert_str_columns_to_float(predict_df)\n",
    "predict_df = predict_df.rename(columns={'index':'post_index'})\n",
    "predict_df = truncate(predict_df)\n",
    "print(predict_df.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(422206, 529)\n"
     ]
    }
   ],
   "source": [
    "term_counts, author_term_counts, vocab = load_term_counts()\n",
    "print(term_counts.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(422206, 119)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "imbalanced_terms = filter_imbalanced_terms(predict_df, term_counts,imbalance=0.1)\n",
    "term_counts = term_counts[:,imbalanced_terms]\n",
    "term_counts.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model accuracy: 0.6995408944057133\n",
      "Mean and std. of the word coeffs: 0.014636964965330776 0.1787484819598313\n",
      "Model accuracy: 0.5842118687298079\n",
      "Mean and std. of the word coeffs: -0.04361795603189822 0.30104023302986554\n",
      "Model accuracy: 0.6990977187606401\n",
      "Mean and std. of the word coeffs: 0.010864669181085081 0.1690052283496078\n",
      "Model accuracy: 0.5815457950289411\n",
      "Mean and std. of the word coeffs: -0.043216267456407616 0.29639779161934204\n",
      "Model accuracy: 0.6968472109943411\n",
      "Mean and std. of the word coeffs: 0.014937959245977598 0.19088597155457082\n",
      "Model accuracy: 0.5836063481257712\n",
      "Mean and std. of the word coeffs: -0.04322841833286749 0.3109966014721763\n",
      "Model accuracy: 0.7000980141481292\n",
      "Mean and std. of the word coeffs: 0.00850213760582418 0.174732192234564\n",
      "Model accuracy: 0.5835037927213841\n",
      "Mean and std. of the word coeffs: -0.04480548729197767 0.3087422986468522\n",
      "Model accuracy: 0.7000170010200613\n",
      "Mean and std. of the word coeffs: 0.020032151570105128 0.19107010356078805\n",
      "Model accuracy: 0.5848775926555594\n",
      "Mean and std. of the word coeffs: -0.034861372132868176 0.3183590673866064\n",
      "Model accuracy: 0.7026333789329685\n",
      "Mean and std. of the word coeffs: 0.0072631104841584734 0.19627805515658872\n",
      "Model accuracy: 0.5858413132694938\n",
      "Mean and std. of the word coeffs: -0.04978001188782236 0.3224292810308951\n",
      "Model accuracy: 0.697620565115785\n",
      "Mean and std. of the word coeffs: 0.00869816889002352 0.1813066252122628\n",
      "Model accuracy: 0.5841300191204589\n",
      "Mean and std. of the word coeffs: -0.048224683103891906 0.3044934316196951\n",
      "Model accuracy: 0.7019108550122871\n",
      "Mean and std. of the word coeffs: 0.012921389459577607 0.18686123238050908\n",
      "Model accuracy: 0.5866875688500974\n",
      "Mean and std. of the word coeffs: -0.048175683442891686 0.31054683473116257\n",
      "Model accuracy: 0.701049507049719\n",
      "Mean and std. of the word coeffs: 0.014132119938180088 0.18343407160710576\n",
      "Model accuracy: 0.5836319304569066\n",
      "Mean and std. of the word coeffs: -0.03971037268957745 0.3078069855803716\n",
      "Model accuracy: 0.6985635781120706\n",
      "Mean and std. of the word coeffs: 0.01358281299872072 0.17582798078810846\n",
      "Model accuracy: 0.5842651333517218\n",
      "Mean and std. of the word coeffs: -0.03927689058997102 0.30564552567431263\n"
     ]
    }
   ],
   "source": [
    "n_bootstraps = 10\n",
    "n_w = term_counts.shape[1]\n",
    "adjusted = np.zeros((n_bootstraps, n_w))\n",
    "permuted = np.zeros((n_bootstraps, n_w))\n",
    "unadjusted = np.zeros((n_bootstraps, n_w))\n",
    "\n",
    "for i in range(n_bootstraps):\n",
    "    sample = assign_split(predict_df,num_splits=2)\n",
    "    sample = sample[sample.split==1]\n",
    "    indices = sample.post_index.values\n",
    "    labels = sample.treatment.values\n",
    "    words = term_counts[indices, :].toarray()\n",
    "    propensity_score = logit(sample.treatment_probability.values)\n",
    "    all_features = np.column_stack((propensity_score, words))\n",
    "    \n",
    "    adjusted[i,:] = fit_treatment(all_features, labels, verbose=True)\n",
    "    unadjusted[i,:] = fit_treatment(words, labels, coeff_offset=0, verbose=True)\n",
    "    np.random.shuffle(words)\n",
    "    permuted_features = np.column_stack((propensity_score, words))\n",
    "    permuted[i,:] = fit_treatment(permuted_features, labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEKCAYAAAACS67iAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xl8VNX5+PHPmewbgUBYJEgA2QkQNonsKIuCigtVqghSRW3VbzcK1rp9f7bWr627rbUqILVKxbq07iIoKMouIDsSkC2EANmTmcw8vz/uzDCByTqTwGSe94t5zc295577zJ3hmTPn3nuuERGUUkqFPtvZDkAppVRwaEJXSqkmQhO6Uko1EZrQlVKqidCErpRSTYQmdKWUaiI0oSulVBOhCV0ppZoITehKKdVERDbmxlq1aiXp6emNuUmllAp569atOyYiqTWVa9SEnp6eztq1axtzk0opFfKMMftqU067XJRSqonQhK6UUk2EJnSllGoiNKErpVQToQldKaWaCE3oSinVRGhCV0qpJkITulIAe1fA7k/PdhRKBaRRLyxS6pxkL4GFk63pe3MgKvbsxqNUPdXYQjfGvGyMOWqM2eJn2a+NMWKMadUw4SnVCA5tODWdc8bHXKmQUZsulwXAxNNnGmM6AOOA/UGOSanGdfz7U9N5u89eHEoFqMaELiJfAMf9LHoC+A0gwQ5KqUZ10meYjBO1GjJDqXNSvQ6KGmOuAA6KyLe1KDvbGLPWGLM2Nze3PptTqmEV5UBiW4hrAcX6GVWhq84J3RgTD9wL3F+b8iLygogMEpFBqak1jv6oVOMrOQ7xLSG+FZQcO9vRKFVv9TnLpQvQCfjWGAOQBqw3xgwRkSPBDE6pRlGSB/Ep4HJCsSZ0FbrqnNBFZDPQ2vO3MSYbGCQi+j9BhabiY9CmN4gLju0829EoVW+1OW3xNWAV0N0Yc8AY85OGD0upRlSSZ3W5JLTSPnQV0mpsoYvItBqWpwctGqUam8sJpSesZA5Wf7rLCbaIsxuXUvWgl/6r8FZ6EpBTB0URK8ErFYI0oavwVpJnPcelQEySNV1ecPbiUSoAOpaLCm/lhdZzbDI47ZXnKRViNKGr8OZpjcckgtPhnld09uJRKgCa0FV4s7uTd0ySttBVyNOErsKbJ3lHJ4Krwpq2awtdhSZN6Cq8eRJ6TDNwebpc9KCoCk2a0FV48yZ03z507XJRoUkTugpv5YUQEQ2RMdYzRg+KqpCl56Gr8GYvOnX+uTFW14u20FWI0oSuwlt5oXVA1CMmURO6Clma0FV4Ky+0WuUeMUl6UFSFLE3oKryVF57qcgGIigdH6dmLR6kAaEJX4a280Opm8dCErkKYJnQV3nwPigJExYGj5OzFo1QANKGr8Hb6QdFobaGr0KUJXYW38tNb6PHaQlchSxO6Cl8uJziKK5/lEhWnLXQVsmpzT9GXjTFHjTFbfOY9ZozZbozZZIx5yxjTvGHDVKoB+F7276EHRVUIq00LfQEw8bR5nwB9RKQvsBO4J8hxKdXw7MXWs28fuuegqMjZiUmpANSY0EXkC+D4afM+FhH3WKN8DaQ1QGxKNSxvQk84NS8qDsR5aqAupUJIMPrQZwEfBKEepRqXw19Cj3cv0wOjKvQElNCNMfcCFcCr1ZSZbYxZa4xZm5ubG8jmlAquqlrooP3oKiTVO6EbY2YAk4EbRKrucBSRF0RkkIgMSk1Nre/mlAo+vwldW+gqdNVrPHRjzERgLjBKRPSTr0KTJ6FHaQtdNQ21OW3xNWAV0N0Yc8AY8xPgWSAJ+MQYs9EY83wDx6lU8GmXi2piamyhi8g0P7NfaoBYlGpc2uWimhi9UlSFL7v7VnO+CT0y1nrWhK5CkCZ0Fb4cJWCLdN9L1M2T0CvKz05MSgVAE7oKX/Ziq3VuzKl5kTHWsyZ0FYI0oavwZS+qfIYL+LTQyxo/HqUCpAldhS9PC92XdrmoEKYJXYUve4mfhO7ucnFqQlehRxO6Cl/24sojLYJ2uaiQpgldhS97kXXLOV8RkWAitMtFhSRN6Cp8+etDB6vbRVvoKgRpQlfhy+GnDx3cCV1b6Cr0aEJX4cvfaYtg9aNrC12FIE3oKnxV2+WiLXQVejShq/BUYQdXRRUJPVYTugpJmtBVePIOzJV45jJtoasQpQldhSfv0LnxZy7TPnQVojShq/Dkbyx0D22hqxClCV2FJ4cnofvpconQ89BVaNKErsKT936i/rpctIWuQpMmdBWequ1y0T50FZpqc5Pol40xR40xW3zmpRhjPjHG7HI/t2jYMJUKMns1XS562qIKUbVpoS8AJp42bx6wVES6AkvdfysVOmo6KKrD56oQVGNCF5EvgOOnzb4SWOieXghMCXJcSjWsGk9b1ISuQk99+9DbiMhhAPdz6+CFpFQj8B4U1dEWVdPR4AdFjTGzjTFrjTFrc3NzG3pzStWOvQgioiEy+sxlkTHgtIPL1fhxKRWA+ib0HGNMOwD389GqCorICyIySEQGpaam1nNzSgWZo8T/KYugt6FTIau+Cf1dYIZ7egbwTnDCUaqR+Lv9nIfehk6FqNqctvgasArobow5YIz5CfBHYJwxZhcwzv23UqHDXuT/DBc41ULXA6MqxETWVEBEplWx6OIgx6JU47FXcbci8Gmha0JXoUWvFFXhqaqbW4C20FXI0oSuwlO1XS7ah65CkyZ0FZ60ha6aIE3oKjxVd9pihCehawtdhRZN6Co8VXseuh4UVaFJE7oKT44yiIr1vyxSW+gqNGlCV+HH5bSuAo2M879cD4qqEKUJXYUfT6KuqYXutDdOPEoFiSZ0FX4c7oSuLXTVxGhCV+GnotR6rrEPXQ+KqtBS46X/SjUWh8PBgQMHKCtr4Jax0wET/gWRLWHbtjOXi1jLY5r7X65UA4mNjSUtLY2oqKh6ra8JXZ0zDhw4QFJSEunp6RhjGm5DjhLIdUKLThDX/MzlInC4HJLaQlK7hotDKR8iQl5eHgcOHKBTp071qkO7XNQ5o6ysjJYtWzZsMgcrYQOYKj7+xgAGRG9woRqPMYaWLVsG9AtVE7o6pzR4ModTibq6bRnbqcSvVCMJ9POvCV2FH29CP/PjHxERQf/+/ek/7kf0Hz6B7Oxs1q5dy9133w3A8uXL+eqrrxo8xOXLlzN58uR6l9mwYQO33HJLvbZ98uRJ/vKXv9Rr3cbw9ttvs3Xr1jqvl5ho3dAkNzeXiRMnBjusc4ImdBV+qulyiYuLY+PGjWxc+iYbv3if9PR0Bg0axNNPPw00XkIP1B/+8Afuuuuueq3bkAm9oqIi4Drqm9A9UlNTadeuHV9++WXAsZxrNKGr8FObLheMN/F7WsLZ2dk8//zzPPHEE/Tv358VK1ZUWuPBBx9kxowZjB8/nvT0dP7973/zm9/8hoyMDCZOnIjD4QBg6dKlZGZmkpGRwaxZsygvt06P/PDDD+nRowfDhw/n3//+t7fe4uJiZs2axeDBg8nMzOSdd6q/42NhYSGbNm2iX79+3rhmzZrF6NGj6dy5s/fLCeDxxx+nT58+9OnThyeffBKAefPmsWfPHvr378+cOXMq1Z2dnU2PHj2YMWMGffv25dprr6WkpASAdevWMWrUKAYOHMiECRM4fPgwAKNHj+a3v/0to0aN4qmnnmLmzJnccccdjBkzhs6dO/P5558za9YsevbsycyZM73b8rSoAZYsWcLMmTP56quvePfdd5kzZw79+/dnz5497Nmzh4kTJzJw4EBGjBjB9u3bAdi7dy9ZWVkMHjyY++67r9LrmDJlCq+++mq1+zEkiUijPQYOHChKVWXr1q2Ns6GiXJGD60Uc5Wcsstls0q9fP+nXu4dMuWyciIgsW7ZMJk2aJCIiDzzwgDz22GN+q33ggQdk2LBhYrfbZePGjRIXFyfvv/++iIhMmTJF3nrrLSktLZW0tDTZsWOHiIhMnz5dnnjiCe/8nTt3isvlkqlTp3q3ec8998iiRYtEROTEiRPStWtXKSoqqhSXr88++0yuvvrqSnFlZWVJWVmZ5ObmSkpKitjtdlm7dq306dNHioqKpLCwUHr16iXr16+XvXv3Su/evf2+xr179wogK1euFBGRm2++WR577DGx2+2SlZUlR48eFRGR119/XW6++WYRERk1apTccccd3jpmzJgh1113nbhcLnn77bclKSlJNm3aJE6nUwYMGCAbNmwQEZGEhATvOm+88YbMmDHDu/4bb7zhXTZ27FjZuXOniIh8/fXXMmbMGBERufzyy2XhwoUiIvLss89Wqu/AgQPSp08fv6/xbPP3/wBYK7XIsXraojonPfSf79h6qCCodfY6rxkPXN672j50T5cLuTvAVvf/HpdeeilRUVFkZGTgdDq9fbUZGRlkZ2ezY8cOOnXqRLdu3QCYMWMGzz33HKNHj6ZTp0507doVgBtvvJEXXngBgI8//ph3332XP/3pT4B1NtD+/furjOHw4cOkpqZWmjdp0iRiYmKIiYmhdevW5OTksHLlSq666ioSEqxx4a+++mpWrFjBFVdcUe1r7NChA8OGDfPG+fTTTzNx4kS2bNnCuHHjAHA6nbRrd+qUz+uuu65SHZdffjnGGDIyMmjTpg0ZGRkA9O7dm+zsbPr3719tDB5FRUV89dVXTJ061TvP84vnyy+/5M033wRg+vTpzJ0711umdevWHDp0qFbbCCUBJXRjzC+AWwABNgM3i4heL63ObbU+y6Xupy3GxFhXmdpsNqKiorxnLdhsNioqKpBqzpyp6gwHEeHNN9+ke/fulebn5OT4LR8XF3fGqW+euMA68FtTLNU5PU5jDCJC7969WbVqld91PF8ap8djs9kqxebZT6dvp6pT+VwuF82bN7e+hGsRq299cXFVDP0Qwuqd0I0x7YG7gV4iUmqM+RdwPbAgSLGpMPbA5b0brvKazkMHK9m7zkzoSUlJFBTU/5dDjx49yM7OZvfu3VxwwQUsWrSIUaNG0aNHD/bu3cuePXvo0qULr732mnedCRMm8Mwzz/DMM89gjGHDhg1kZmZWuY2ePXvy5z//ucZYRo4cycyZM5k3bx4iwltvvcWiRYtISkqisLCwyvX279/PqlWryMrK4rXXXmP48OF0796d3Nxc73yHw8HOnTvp3bv+72ObNm3Ytm0b3bt356233iIpKQmgUnzNmjWjU6dOvPHGG0ydOhUR8R4/GDZsGK+//jo33njjGf3lO3fupE+fPvWO7VwV6EHRSCDOGBMJxANN7zeManrEBZgaDor6b6FffvnlvPXWW34PitZGbGws8+fPZ+rUqWRkZGCz2bj99tuJjY3lhRdeYNKkSQwfPpyOHTt617nvvvtwOBz07duXPn36nHGA73Q9evQgPz+/2qQMMGDAAGbOnMmQIUO48MILueWWW8jMzKRly5YMGzaMPn36nHFQFKwvjIULF9K3b1+OHz/OHXfcQXR0NEuWLGHu3Ln069eP/v37B3w20B//+EcmT57M2LFjK3XfXH/99Tz22GNkZmayZ88eXn31VV566SX69etH7969vQeNn3rqKZ577jkGDx5Mfn5+pbqXLVvGpEmTAorvXGTq+7MLwBjzP8DvgVLgYxG5wU+Z2cBsgPPPP3/gvn376r091bRt27aNnj17NvyG8n+AkhPQrm/VZY7vtUZbbN0I8TSAJ554gqSkpHqfi16V7OxsJk+ezJYtW4Jab2MbOXIk77zzDi1atDjboZzB3/8DY8w6ERlU07r1bqEbY1oAVwKdgPOABGPMjaeXE5EXRGSQiAw6/UCNUmeFSPXdLVDvPvRzxR133FGpb1qdkpubyy9/+ctzMpkHKpAul0uAvSKSKyIO4N/ARcEJS6kGJK4auluwlofwpf+xsbFMnz496PWmp6eHfOs8NTWVKVOmnO0wGkQgCX0/MNQYE2+sQ8kXAzrWqDr3iavJt9BVeKp3QheRb4AlwHqsUxZtwAtBikuphlOrhB7aLXQVngI6D11EHgAeCFIsSjWO2vShYwNc7rKNMAKkUkGgY7mo8FPbPnSrcIOHo1SwaEJX4aeaLhfP8Ll9LhzN1Nm/oaS4uFFDe/LJJ72DXdVWbYbaDVR2dnaNF+LUpoxqWJrQVfgRF1V99D1juWxZs4Lo6Cief/75WlfrdDoDDq0+CV0pD03oKvyIgK3ms1xGDMlk9+7dAPzjH/9gyJAh9O/fn9tuu82bvBMTE7n//vu58MILWbVqFenp6fz2t78lKyuLQYMGsX79eiZMmECXLl28Xw6nt6jvvPNOFixYwNNPP82hQ4cYM2YMY8aMAayBubKyshgwYABTp06lqKgIqHqoXV8LFixgypQpXH755XTq1Ilnn32Wxx9/nMzMTIYOHcrx48cB2LhxI0OHDqVv375cddVVnDhxArCGw+3Xrx9ZWVk899xz3nqdTidz5sxh8ODB9O3bl7/97W91fQdUA9GErsJPLfrQKyqcfLDsSzL69Gbbtm0sXryYL7/8ko0bNxIREeEdG6S4uJg+ffrwzTffMHz4cMAajXDVqlWMGDGCmTNnsmTJEr7++mvuv//+ard59913c95557Fs2TKWLVvGsWPHePjhh/n0009Zv349gwYN4vHHH6esrIxbb72V//znP6xYsYIjR45UWeeWLVv45z//yerVq7n33nuJj49nw4YNZGVl8corrwBw00038eijj7Jp0yYyMjJ46KGHALj55pt5+umnzxhw66WXXiI5OZk1a9awZs0a/v73v7N3797q97lqFDp8rjo3fTAPjmwObp1tM+DSP1bb5VJaWmoN3epyMmJQb35y8wxemL+IdevWMXjwYG+Z1q1bA1af+zXXXFOpDs/wsxkZGRQVFZGUlERSUhKxsbGcPHmy1uF+/fXXbN261TtUrd1uJysri+3bt1c51O7pxowZ491+cnIyl19+uTe2TZs2kZ+fz8mTJxk1ahRgDec7derUM+ZPnz6dDz74ALB+NWzatIklS5YAkJ+fz65du7xDAquzRxO6Ci8iQNWnLXrHQy/Lh+PfQ1QkIsKMGTN45JFHzigfGxtLREREpXk1DQ0bGRmJy2ckx6qGhhURxo0bV2nkRbC6SGp7M+HTt+8bW3W3gxORaofzfeaZZ5gwYUKl+dnZ2bWKSTUcTejq3HTpHxumXldtbj/HqYQvwsUXX8yVV17JL37xC1q3bs3x48cpLCysNCJiXXTs2JGtW7dSXl5OWVkZS5cu9XbXeIaGbdWqFUOHDuVnP/uZd6jdkpISDhw4UO1Qu3WVnJxMixYtWLFiBSNGjPAO59u8eXOSk5NZuXIlw4cPrzT87IQJE/jrX//K2LFjiYqKYufOnbRv377eMajg0YSuwks1dyuqzHjL9+rVi4cffpjx48fjcrmIioriueeeq3dC79ChAz/60Y/o27cvXbt2rTS2+ezZs7n00ktp164dy5YtY8GCBUybNs17F56HH36Ybt26eYfabdWqFcOHDw9ofJWFCxdy++23U1JSQufOnZk/fz4A8+fPZ9asWcTHx1dqjd9yyy1kZ2czYMAARITU1FTefvvtem9fBU9Aw+fW1aBBg2Tt2rWNtj0VWhpl+NwKOxz9DpI7QEKrqsvZS+DYDmjRGeKSGzYmpXycleFzlQpNtWyhe7tkdIAuFTo0oavw4r39XG370DWhq9ChCV2Fl9r2oXsSvo64qEKIJnQVXmp9UFRb6Cr0aEJX4cXb5aItdNX0aEJX4UXqeB66HhRVIUQTugov3i4U/x99Y4x1L05jAEOFw0Fqamqdh6dNT0/n2LFjAZepytq1a7n77rvrvN5ll11Wp+EHQtktt9zC1q1bAfjDH/5wlqNpHJrQVXipoQ89ISGBLVu2UFpaCsbwyWefn5NXQQ4aNIinn366zuu9//77NG/evAEi8i8YQwrX14svvkivXr2As5/QqxtmIZgCSujGmObGmCXGmO3GmG3GmKxgBaZUg6hFH/qll17Ke++9B8bGa0veZtq0ad5lx48fZ8qUKfTt25ehQ4eyadMmAPLy8hg/fjyZmZncdttt+F6wV9XQu1VJTExk7ty5DBw4kEsuuYTVq1czevRoOnfuzLvvvgtUHoL3888/p3///vTv35/MzEwKCws5fPgwI0eOtG7W0acPK1asAE79KsjOzqZnz57ceuut9O7dm/Hjx1tfYsCaNWvo27cvWVlZzJkzx+9NK5YvX87IkSO56qqr6NWrF7fffrt3fJrThxReunQpmZmZZGRkMGvWLO9Vr+np6cydO5chQ4YwZMgQ71DFubm5XHPNNQwePJjBgwfz5ZdfAvDggw8ya9Ys777wfKEVFxczadIk+vXrR58+fVi8eDEAo0ePZu3atcybN8876NoNN9zAfffdx1NPPeV9Lffee+8ZX45V1blmzRouuugi+vXrx5AhQygsLKSsrIybb76ZjIwMMjMzWbZsGWANXzx16lQuv/xyxo8fD8Bjjz3mHXb4gQca4O6dIlLvB7AQuMU9HQ00r678wIEDRamqbN26teE3UnBE5OB6EWeF38UJCQny7bffyjXXXCOle9dKvz49ZdmyZTJp0iQREbnzzjvlwQcfFBGRpUuXSr9+/URE5K677pKHHnpIRET++9//CiC5ubmydetWmTx5stjtdhERueOOO2ThwoUiItKxY0fJzc09IwZA3n//fRERmTJliowbN07sdrts3LjRuz3fmCZPniwrV64UEZHCwkJxOBzypz/9SR5++GEREamoqJCCgoJK29y7d69ERETIhg0bRERk6tSpsmjRIhER6d27t3z55ZciIjJ37lzp3bv3GTEuW7ZMYmJiZM+ePVJRUSGXXHKJvPHGG974Fy9eLCIipaWlkpaWJjt27BARkenTp8sTTzzhjcUT48KFC72vZ9q0abJixQoREdm3b5/06NFDREQeeOABycrKkrKyMsnNzZWUlBSx2+2yZMkSueWWW7yxnTx5UkRERo0aJWvWrPG+rx579+6VzMxMERFxOp3SuXNnOXbsWKXX56/O8vJy6dSpk6xevVpERPLz8737eubMmSIism3bNunQoYOUlpbK/PnzpX379pKXlyciIh999JHceuut4nK5xOl0yqRJk+Tzzz8/Y9/6+38ArJVa5OR6j+VijGkGjARmur8Y7IA9kC8XpTweXf0o249vD2qdPVJ6MLfHTdYf1bTQ+/btS3Z2Nq+9/SGXjRtTadnKlSt58803ARg7dix5eXnk5+fzxRdfeG80MWnSJFq0aAHA0qVLqxx6tyrR0dFMnDgRsIa5jYmJISoqioyMDL8jGg4bNoxf/vKX3HDDDVx99dWkpaUxePBgZs2ahcPhYMqUKdaQwKfp1KmTd/7AgQPJzs7m5MmTFBYWctFFFwHw4x//mP/+979+4xwyZAidO3cGYNq0aaxcuZJrr7220pDCO3bsoFOnTt6hdWfMmMFzzz3Hz3/+c+96nudf/OIXAHz66afevm+AgoICCgsLvfs2JiaGmJgYWrduTU5ODhkZGfz6179m7ty5TJ48mREjRlS7f9PT02nZsiUbNmwgJyeHzMxMWrZsWamMvzo3b95Mu3btvO9ls2bNAOszcddddwHQo0cPOnbsyM6dOwEYN24cKSkpgDXs8Mcff+wdu6eoqIhdu3YxcuTIauOti0AG5+oM5ALzjTH9gHXA/4hI496EUak6cQGmxrNcrrjiCn79v4+x/N3XyPPp/hQ/pzF6hpn1N9ysVDP0blWioqK8ddVmyNt58+YxadIk3n//fYYOHcqnn37KyJEj+eKLL3jvvfeYPn06c+bM4aabbqq0nu/QuhEREZSWlvp9fVU5/fV6/vYdUrim+nzr8Ey7XC5WrVpFXFzcGeVPj7miooJu3bqxbt063n//fe655x7Gjx9f481EbrnlFhYsWMCRI0eYNWvWGcv91TllypQq3+OqJCQkVCp3zz33cNttt1UbWyACSeiRwADgLhH5xhjzFDAPuM+3kDFmNjAb4Pzzzw9gcyqczB0yt2Eqzj9Qi4uKYNasWSRHlJPRsxvLN//gnT9y5EheffVV7rvvPpYvX06rVq1o1qyZd/7vfvc7PvjgA+9t3II99K4/e/bsISMjg4yMDFatWsX27duJi4ujffv23HrrrRQXF7N+/fozEro/LVq0ICkpia+//pqhQ4fy+uuvV1l29erV7N27l44dO7J48WJmz559RpkePXqQnZ3tHQLYMzyvx+LFi5k3bx6LFy8mK8s6BDd+/HieffZZ5syZA1jjv/v7heFx6NAhUlJSuPHGG0lMTGTBggVnlImKisLhcBAVFQXAVVddxf3334/D4eCf//xnreqcN28ehw4dYs2aNQwePJjCwkLi4uK87/3YsWPZuXMn+/fvp3v37qxfv75SnRMmTOC+++7jhhtuIDExkYMHDxIVFVXjL7a6CCShHwAOiMg37r+XYCX0SkTkBeAFsEZbDGB7SgVOXLVK6GlpafzPbTPOuFL0wQcf5Oabb6Zv377Ex8ezcOFCAB544AGmTZvGgAEDGDVqlLfxEuyhd/158sknWbZsGREREfTq1YtLL72U119/nccee4yoqCgSExO9t5urjZdeeolbb72VhIQERo8eTXKy/9Ems7KymDdvHps3b/YeID1dbGws8+fPZ+rUqVRUVDB48GBuv/127/Ly8nIuvPBCXC6Xd1z3p59+mp/97Gf07duXiooKRo4cWe3Nujdv3sycOXOw2WxERUXx17/+9Ywys2fPpm/fvgwYMIBXX32V6OhoxowZQ/Pmzc+4QUlVdUZHR7N48WLuuusuSktLiYuL49NPP+WnP/0pt99+OxkZGURGRrJgwYJKvyQ8xo8fz7Zt27xfXImJifzjH/8IakIP9KDoCqC7e/pB4LHqyutBUVWdRjkoenyvyJEttSt7bJfI0e0NGs65qLCw0Dv9yCOPyN13331GGd+DsvVV1UHhxuB0OqVfv36yc+fOs7L96gRyUDTQ89DvAl41xmwC+gPhcfa+Cl21bKEDVrkwHMvlvffeq3S64+9+97uzHVJQbd26lQsuuICLL77Ye1/WpkJvcKHOGY1yg4u8PeCqgNTuNZc9kQ32YmjTu2FjUsqH3uBCqdoSV83juHgYmw7OpUJvLc8IAAAgAElEQVSKJnR1TmnwX4za5aLOYYF+/jWhq3NGbGwseXl5DZvUxUWtP/aa0FUjEhHy8vKIjY2tdx2BnLaoVFClpaVx4MABcnNzG24jBYchMhqOltdctizfepzcWvtuGqUCEBsbS1paWr3X14SuzhlRUVF06tSpYTfy56uhy1iY8lzNZVc9Bx/9Fubth1j/52IrdS7RLhcVXipKIaqWP2mj3JeeO0obLh6lgkgTugovjrJTibomUfHudUoaLh6lgkgTugofIlYLPbK2CV1b6Cq0aEJX4aPCfSC01l0unha6JnQVGjShq/Dh6TrxJOqaeFvo2uWiQoMmdBU+PC3tWveha5eLCi2a0FX48Cb0hOrLeXjKaQtdhQhN6Cp8ONw309IWumqiNKGr8FHnLhd3X7td76qoQoMmdBU+PF0n0bXtctEWugotmtBV+NCDoqqJ04Suwoe9jqct2iIgIkYPiqqQoQldhQ/veei1bKF7ymoLXYUITegqfHi7XGrZQveU1Ra6ChEBJ3RjTIQxZoMx5r/BCEipBlPXK0XB3ULXhK5CQzBa6P8DbAtCPUo1LEcJYCAypvbrxCRBeVGDhaRUMAWU0I0xacAk4MXghKNUA3KUWqcs1uXuQzFJUF7YcDEpFUSBttCfBH4D6I0X1bnPUVK3A6IA0Ylg14SuQkO9E7oxZjJwVETW1VButjFmrTFmbYPeK1KpmtjrkdC1y0WFkEBa6MOAK4wx2cDrwFhjzD9OLyQiL4jIIBEZlJqaGsDmlAqQo6T2A3N5xCRql4sKGfVO6CJyj4ikiUg6cD3wmYjcGLTIlAo2R2n9Wuh2baGr0KDnoavw4Sit2ymLANFJUFEGTkfDxKRUEAUloYvIchGZHIy6lGowjuL6tdBBu11USNAWugofjlKIrmMLPSbRetaErkKAJnQVPhwlde9y8bTQtR9dhQBN6Cp81OegaLS20FXo0ISuwoc9gBa6nouuQoAmdBUeXE7roGhMs7qt503oBcGPSakg04SuwoOny8SToGvL0+WifegqBGhCV+GhvgldT1tUIUQTugoPni6T+rbQtQ9dhYDIsx2AUo2ihhb6yRI7O3OKKHM4iY+OoGubJJLjoiAi0jqQqn3oKgRoQlfhwZPQY5O9s0rtThav2c+b6w+y+WD+Gav069Ccawa0Z3pMMqbsZGNFqlS9aUJX4aHMnbDdLfR3Nh7kD+9vI6egnIz2ycyZ0J1e5zWjWWwk+aUOvjtYwAdbjnD/O9+RFRtFwuFDtBPB1OXmGEo1Mk3oKjy4W+hFxPG71zfw9sZD9EtL5unrM7mwc8szio/t0Ya7Lu7K19/nUfZac04cOsj9r6zjz1P7kRwf1djRK1UrelBUhQd3Qp/+j638Z9NhfjmuG2/ecZHfZO5raOeW9Lkgna6Jdj7feZQrn1vJ7qN6xos6N2lCV2HhxMk8XBh2nYQFNw/m7ou7EhlRu4+/iU+hhSnitVuHUlRewdV/+YqNP2ifujr3aEJXTd4Px0v4ZP1Oiolj8W1ZjOhaxztnxaVAyXEGdWzBWz8dRvP4aG588RvWZh9vmICVqidN6KpJO5Jfxo9f/JoEZwExzVLpfV5yzSudLj4FxAll+XRIiWfxbUNpnRTDTS+v5lttqatziCZ01WQVlVcwc/5qThQ7GNneRnRSPe9pG5diPZdaLfJ2yXG8PnsoLROjuXnBGvYeKw5SxEoFRhO6apJcLuHnr29k19Ei/nrjAJJc+RBf/QHQKsW7E3rJCe+s1s1iWXjzEABuevkbjhaWBRqyUgHThK6apD99vINPt+Vw36SeVp95yfH6J/TTWugenVMTmT9zMHlFdma/so4yhzPAqJUKTL0TujGmgzFmmTFmmzHmO2PM/wQzMKXq652NB/nL8j1MG3I+My5Kt2YWHzvV0q4rzxdByZkHQft1aM7jP+rPxh9O8tt/b0ZE6rcNpYIgkBZ6BfArEekJDAV+ZozpFZywlKqfjT+cZM6STQzplMJDV/S2ruy0l0BFKSS0ql+l8f5b6B4T+7TlF5d0498bDvLiir31jFypwNU7oYvIYRFZ754uBLYB7YMVmFJ1dSS/jNmvrKV1UgzP3ziQ6Ej3x7skz3qub5dLbHOwRUJRTpVF7hp7AZf2acsjH2xj+Y6j9duOUgEKSh+6MSYdyAS+8bNstjFmrTFmbW5ubjA2p9QZSu1OZi9aS3F5BS/NGExKQvSphUXuBJtQz7NcbDZIagcFh6opYvjzj/rRvW0z7nptA9/n6nC7qvEFnNCNMYnAm8DPReSMMUZF5AURGSQig1JT6/kfSqlqiAhzlnzL5oP5PHV9Jt3bnjZEbsEB67lZAD8gm7WvNqEDxEdH8sL0gUTaDLMXraOwzFH/7SlVDwEldGNMFFYyf1VE/h2ckJSqm6eW7uK/mw7zmwk9uKRXmzML5B+0npPT6r+RZudB/oEai3VIiee5Gwaw91gxv1i8EZdLD5KqxhPIWS4GeAnYJiKPBy8kpWrvv5sO8eSnu7h6QHtuH9XZf6H8A9ZNKuJa1H9Dye4Wei3OYrmoSyt+N6knn247ypOf7qz/NpWqo0Ba6MOA6cBYY8xG9+OyIMWlVI02HTjJr/71LQM7tuCRqzOqHqu84IDVZRLIWObN0sBZDsW1Ow4086J0rh2YxtOf7ebDLYfrv12l6qDe46GLyEpAR/tXZ8X+vBJ+snAtrRJj+Nv0gcRERlRd+OT+wLpbAFp1tZ5zt0Ni6xqLG2N4eEofdh0t4pf/+pb0Vgn0aNsssBiUqoFeKapCTm5hOdNf/gZ7hYv5Nw+mVWJM1YVFIHcnpHYPbKNtelvPOVtrvUpsVAQvTB9IQkwks19Zx8kSe2AxKFUDTegqpBSUOZjx8mqOFpTz8szBdGvj/6bPXvkHwFEceEJPbGP1wR/9rk6rtWkWy/M3DuRIfhl3/nMDFU5XYHEoVQ29BZ0KGQVlDma+vJqdOYW8OGMQAzvW4iBn7nbrObVHYBs3BtoPxJW9guwTe/ju+FYOFh3kUNEhckpyKHGUUFJRgktcxEXGERcZR8vYlqQlpdEhqQN3TmzB4+8d5dEPt3PvJL2gWjUMTegqJJwssXPTy6vZdriAZ3+cyejuNfdjA3Boo/Xcume9tlviKOG7vO/YeHQjG+PsfJtoJ//dKd7lqXGptIlvQ0J0Ai1iWxBhIiitKKWkooTNxzbz8b6PcYo1aFeLnrH8IzuNnPdHM2f41XRo1qFeMSlVFU3oqtG4xEWhvZDSilLvPJuxER8ZT3xUPDbjvwcwp6CMmfPXsOdoEc/fOJCLe/o517wq+1ZC6961PmXxSPERNh7dyIajG9iYu5Edx3d4E3KXpI5cUmqnX0pP+l76JGlJacREVNN/DzhcDg4VHWLzsc1syNnIf3asZFnuiyx760W6tujKlV2uZFLnSbSKq+c4M0r5MI05OtygQYNk7dq1jbY91fhOlp1k2/Ft7Di+g/2F+73dEnmleRQ5ihD8f94MhsToRJKjk2kd35q2CW1pm9AWcSTz2qoCSoqTeOyqUVzaq0vVpyeezumAP54PmTfCZY+dsdjT+t5ybAubj21mU+4mckqs8VriIuPIaJVBv9R+9G/dn36p/UiOSYZlj8Dnf4Tr/wk9JtV5/5TanVz38nvsKvyGrp13sbdoGxEmghHtRzCtxzSyzsuq/etTYcMYs05EBtVYThO6qi8RYX/hfr45/A1rjqxh87HNHCw66F2eHJNM+8T2pCWmkRqfSrPoZiRFJxEfGW+tj+ASF6UVpRTYCyiyF3Gi/AQ5xTnklORwpDgHp1RU2mZcZBxtE9rSLqEd7RLa0SahDe0S2pEcnUxidCLxUfHERcThEhfOw99S8fbt5I+ey/G2PThRdoJDRYfYW7CX7PxsDhUd8n7BpCWmWQm8tZXAu7foTqTNzw/YCju8eLF1kdFd6yCueZ3328kSOz/62yoOnCjl9z9qTXb557yz+x3yyvLoktyFG3rdwOTOk4mLjKtz3app0oSuGkSFq4LVR1bzcfbHrDi4gqMl1sBXreNbk9k6k14te9EzpSc9U3rSPLbuyQ7AXuHijx9s5+Uv99Crg425k9riMCc4XHyYw8WHOVJ8hMNF1nReWV6d6o6LjCO9WTrpzdLplNyJ3q1606dVH1Ji6zBW+qEN8MJoGP97uOjOur04t6MFZfz4xW84eKKUl2YMYlCnZnyU/RGLti5i2/FtJMck86NuP+LHPX+s3TFKE7oKHqfLydqctXyY/SFL9y3lRPkJ4iPjGdZ+GEPbDWVI2yF0bNYxKF0F+/NKuOv1DXz7w0lmXpTOPZf1qPaiIbvTTk5xDgX2AoodxRQ7iilzlmEzNiJWPE5E0TGaXf0CKbEppMSm0Cy6WXC6NF4aD2UF8LOv611FbmE5N774Ddl5xTwzLZPxvdsiIqzLWceirYtY9sMyIm2RTOo8iZt63UTXFl0Dj1uFJE3oKiBOl5P1R9fzUfZHfLLvE46XHScuMo7RaaOZkD6BYe2HERsZG7TtuVzCK6uyefTDHUTaDP93bV8uzWgXwAuogEfToe9UmPxEsMI85eu/wofz4K710LJLvas5Xmzn5gVr2HTgJPdc2oNbR3T2fuHsK9jHoq2LeGf3O5Q5yxh23jBm9J7B0HZDtZ89zGhCV3XmEhcbj270JvHc0lxiI2IZmTaSCekTGJE2okH6dbOPFfObNzexeu9xRnVL5ZGrMziveYDbObge/j4GrnkJMq4NTqC+TmTDU/1gwh8g62cBVVXmcPKrN77lvU2HuWZAGv9vSm/io0/1358sO8m/dv6Lf277J3lleXRr0Y0ZvWdwafqlREVEBfhCVCjQhK5qRUT4NvdbPsr+iI/3fczRkqPERMQwov0IJqRPYGTaSOKj4htk22UOJ3//4nueW76bqAgb90/uxbUD04LT+vzqGfj4d/DL7dAsgJZ+dZ4ZCC06wY1LAq7K5RKeWrqLpz/bRadWCTx9fSZ92idXKmN32nnv+/d4Zesr7D65m9ZxrZnWcxpTu021zsBRTZYmdFUlEeG7vO/4cO+HfLzvYw4XHybKFsXw9sOZkD6B0R1GkxCV0KAxfLY9h4f+s5V9eSVcltGW+yf3pm1y8LpweG2adZXo3RuCV+fp3p8D6xfBvH0QWf356LX11Z5j/GLxRo4V2bn5onR+Pq4biTGVz7YREb469BULvlvA14e/Ji4yjiu6XMGUC6bQu2Vv7Y5pgjShq0pc4mJT7iY+2fcJn+77lEPFh4i0RTLsvGHeJJ4UXcO4KEHwfW4Rv39vG0u3H6VLagIPXdGH4V2DfBaHywX/1wl6Xg5XPhvcun3t+ABeux5uehc6jwpatSdL7Dz64Q5eW72f1kkx3D6qCz++8Hxio848OLzj+A5e2foKH2V/RLmznC7JXbjigiu4rNNltE1oG7SY1NmlCV3hdDnZcHSDlcT3f8rRkqNE2aLIOi+LcR3HMabDmEb7qX60oIwnl+5i8ZofiI208T+XdGXmRZ1O3cg5mI5shueHw1V/g37XB79+j/JC68Br1p0w7qGgV79+/wke/WA73+w9TmpSDDde2JEfDU6jXfKZxxcK7YV8lP0R7+55lw1HrV8lfVr2Yez5Yxl7/lg6J3fWlnsI04QepvJK8/jy0JesPLCSrw5/RX55PjERMQxvP5xLOl7CqLRRjdIS9zhZYufvK77n5ZXZVLhc/HjI+dw5tiupScHpovDLcwbKz7dA8wYeL2X+ZVZiv31Fg23i6+/zeG7ZblbsOobNwLALWjG+d1vG9Wzjt5tqX8E+Ptn3Ccv2L2PTsU0AnJdwHoPbDubCdhcyuO1gbb2HGE3oYaLIXsTG3I2sy1nHV4e+YmueNV53SmwKw9sPZ0TaCEa2b7gDm1U5nF/KSyv28s/V+ymxO7mi33n8anw3OrZs2L55ABZdDfk/wJ1rGn5bXzwGnz0Mv94NiQ17E/T9eSX8a+0P/GfTIfbllQDQvU0Sg9JbMDg9hUHpLWjfPK5SSzynOIflPyxn1eFVrDmyhgK7dR/3dgnt6NWyl/fRI6UHLWNbaiv+HKUJvQkSEQ4UHmDr8a18m/st63LWsf34dlziIsJEkNEqg+HthzM8bTg9U3pWOdhVQ8a3eu9xFq+xko5L4Ip+53HbqM6Nd7cee4nVDTLkVpjw+4bf3sF18PexcPWL1jnvjUBE2JNbxCdbj7Lq+zzW7ztBUbk1REJKQjS92jWj93nN6HVeM3q1a0bn1EQibAaXuNh5YierD69m87HNbM3byv7C/d56k6KSSE8+dRVth2YdaBtvDbPQKq4VEbZq7gqlGlSjJHRjzETgKSACeFFE/lhdeU3otZdfnk92QTbZ+dnsPLGTbce3sT1vO4WOQgBiImLom9qXAa0HMLDNQPql9mv0VjhYyWVHTiEff5fDWxsOsvdYMUkxkVw9oD23jOhMh5RGjmn7e/D6j2H629BlTMNvz+WEx7pAt4lw1fMNvz0/nC5h+5EC1u07wXcHC9h6uIAdRwqxu2+mERtl44LWiXRtnUTXNu7n1ol0SImnuKKQ7Xnb2Xlip/fztrdgr3dIB49IE0lqfKp3/JyWsS1pGdfy1LN7OiU2heiI6LOxG5q0Bk/oxpgIYCcwDjgArAGmiUiV9+jShG4REQrsBeSU5JBTnMORkiPWc/ER9hfuJzs/mxPlJ7zlYyJi6N6iOz1bWmOk9GzZkwuaX3BW/uOICPvySli77wTr9p3gy93H2H/c+vk/JD2F6wZ34LKMdsRFn6XW3KKrIec7+MUWaKyLbv59G2z/L/x8M8TXYUyYBuRwutiTW+RN8LuOFrErp5DD+WXeMjGRNrqkJtK1TSLd2iTRsWU8HVrEk9YijphoB4eKD3Gk+Ij34RlH50jxEfLK8ioNg+wrKTrJm9yTY5JJjkmmeUxz73Ry9JnzYiNitbunGo2R0LOAB0VkgvvvewBE5JGq1gm1hC4iVEgFFS7r4XQ5K/1dVlFGaUWp91HiKKHAXkxheQmF9mIK7cWcKDtJfvlJ8u35FNnzKaoooLiiAKc4TtuaIYbmxNCaaGmDraI1xpFKRXkr7GUtKHcXN8ZgjHV3bmPAZgxRETbioiKIjY4gLspGbFQEce5HjGc62l0mKoK46FPLY92P6EhDhVOocFkPe4WLEyV2jhfbOVFs51B+Gd/nFrH3WDEldmt88GaxkQxOT+Hinm24pGdrWjcL4nnkdSUC374Gb98BFz8AI37ZeNs+ssU6q6bP1XDlcxB17o6SWFjmYPfRInblFLHraKE70Rdx8GTl5BwfHUGHFvG0bxFH66QYUhKiaZkYQ6vEaFISomkWG4XNZsdOAWWufIorTpLvOE5eaR7Hy04959vzyS/Pp6C8gDJnWRVRQbQtmuSYZJKik0iISiA+Kp7EqERrOjKehKgE73zPdExEDFG2KGIiYoiOiLYetmhrfsSp+REmgggTUe0XhsPpotThpMzhpNzhoszhpMzhM6/CRXmFE3uFy5p2OLE7XZQ7rL8rXILNQITNYDOGCJv1iIowxERGEBtlY2S3VL9nKNVGbRN6IDe4aA/84PP3AeDCAOqr0j0LrmJdxS7vSNq+zwYQQ6VRtv1Nn/61dUZdBioAl/vZacAVYIvBJkIzJyS5hCQnnO8UmrmnmzuF1AqhZYXQ2imkuAxR5gg2jlgJ2524bbgTuG+D1zd4EXAKLjtIkTUkrYg1W7AmrOkz9wGAqWJ8ct9lBrDZDJE2G5Exhsh4Q3SE9bfJA1a6H2cE6EeNDYh6russh7J86DDUOo2wMbXtA2N/B5/9P6vLJ6G1daFRIx/DqI0kINP98EoEV4LgcAoOlwuHU6hwunCUuKgotL7gnZ4PEZXfoRh3nalYnxOM8TY2OG1euYEiGxTaoDDCei5yTxfZyimwHaXEdpRSG+QbOGKDUhuUGiixgTPA/49GBBtgE+tmyr7Txne++L/Zsjnt2XdnGD/lTv9bTv6U68f9NKDXUJNAErq/vXvG/zhjzGxgNsD5559frw21jm1Hev7h0zZqKv1tMFXuSHNaWX/TBkMEVgKNxBAhnmkbEQIRxmATQyQQ6d5atNiIwUacGGKxEW8iiDM24k0k8TZrOtJmwxZpiIwyREQYImw2It3f3p7nyi2Haj601X6gq19PEFwuq7/VKYLTZbUqnC6hwgUuAZs5lbxtBqIjI4iJjCDSZqrYdAPEWt91jQ3aD4CMqRB5FvpwR/4aOlwIOz+EkjyoKKvFl9e5w4aVnKs6mVSwWrH2Cqu1WuG0PjuVP0fibUhYjQrxtjlc7ukIoLn74a24wv82fYsIggOhzAhlxkWZERxGqDDWfO+zzedvI1QguAy4jLV9F1ajTRDvNLgbhZ6HzWoqer+UjOcjaTz/Kn1RyWkZ3vOavbG7X3vP886r4V0InHa5KKXUOa62XS6B/CZcA3Q1xnQyxkQD1wPvBlCfUkqpANS7y0VEKowxdwIfYf2SellEvgtaZEoppeokkD50ROR94P0gxaKUUioA595heKWUUvWiCV0ppZoITehKKdVEaEJXSqkmQhO6Uko1EY06fK4xJhfYV8/VWwHHghhOsGhcdaNx1Y3GVTfnalwQWGwdRaTGAfcbNaEHwhiztjZXSjU2jatuNK660bjq5lyNCxonNu1yUUqpJkITulJKNRGhlNBfONsBVEHjqhuNq240rro5V+OCRogtZPrQlVJKVS+UWuhKKaWqcU4ldGPMVGPMd8YYlzGmyqPBxpiJxpgdxpjdxph5PvM7GWO+McbsMsYsdg/rG4y4Uowxn7jr/cQY08JPmTHGmI0+jzJjzBT3sgXGmL0+y/o3Vlzuck6fbb/rM/9s7q/+xphV7vd7kzHmOp9lQd1fVX1efJbHuF//bvf+SPdZdo97/g5jzIRA4qhHXL80xmx175+lxpiOPsv8vqeNFNdMY0yuz/Zv8Vk2w/2+7zLGzGjkuJ7wiWmnMeakz7KG3F8vG2OOGmO2VLHcGGOedse9yRgzwGdZcPeXiJwzD6An0B1YDgyqokwEsAfoDEQD3wK93Mv+BVzvnn4euCNIcf0fMM89PQ94tIbyKcBxIN799wLg2gbYX7WKCyiqYv5Z219AN6Cre/o84DDQPNj7q7rPi0+ZnwLPu6evBxa7p3u5y8cAndz1RDRiXGN8PkN3eOKq7j1tpLhmAs/6WTcF+N793MI93aKx4jqt/F1YQ3o36P5y1z0SGABsqWL5ZcAHWDc6Ggp801D765xqoYvINhHZUUOxIcBuEfleROzA68CVxhgDjAWWuMstBKYEKbQr3fXVtt5rgQ9EpCRI269KXePyOtv7S0R2isgu9/Qh4CjWrSmDze/npZp4lwAXu/fPlcDrIlIuInuB3e76GiUuEVnm8xn6GkgL0rYDiqsaE4BPROS4iJwAPgEmnqW4pgGvBWnb1RKRL7AacFW5EnhFLF8DzY0x7WiA/XVOJfRa8ndz6vZAS+CkiFScNj8Y2ojIYQD3c+sayl/PmR+m37t/bj1hjKnq1o0NFVesMWatMeZrTzcQ59D+MsYMwWp17fGZHaz9VdXnxW8Z9/7Ix9o/tVm3IePy9ROsVp6Hv/e0MeO6xv3+LDHGdKjjug0ZF+6uqU7AZz6zG2p/1UZVsQd9fwV0g4v6MMZ8CrT1s+heEXmnNlX4mSfVzA84rtrW4a6nHZCBdScnj3uAI1hJ6wVgLvC/jRjX+SJyyBjTGfjMGLMZKPBT7mztr0XADBFxuWfXe3/524Sfeae/zgb5TNWg1nUbY24EBgGjfGaf8Z6KyB5/6zdAXP8BXhORcmPM7Vi/bsbWct2GjMvjemCJiDh95jXU/qqNRvt8NXpCF5FLAqziANDB5+804BDWGAnNjTGR7laWZ37AcRljcowx7UTksDsBHa2mqh8Bb4mIw6fuw+7JcmPMfODXjRmXu0sDEfneGLMcyATe5CzvL2NMM+A94Hfun6Keuuu9v/yo6vPir8wBY0wkkIz1E7o26zZkXBhjLsH6khwlIuWe+VW8p8FIUDXGJSJ5Pn/+HXjUZ93Rp627PAgx1SouH9cDP/Od0YD7qzaqij3o+ysUu1z83pxarKMMy7D6rwFmALVp8dfGu+76alPvGX137qTm6beeAvg9Gt4QcRljWni6LIwxrYBhwNazvb/c791bWH2Lb5y2LJj7qzY3M/eN91rgM/f+eRe43lhnwXQCugKrA4ilTnEZYzKBvwFXiMhRn/l+39NGjKudz59XANvc0x8B493xtQDGU/mXaoPG5Y6tO9YBxlU+8xpyf9XGu8BN7rNdhgL57kZL8PdXQx35rc8DuArrW6scyAE+cs8/D3jfp9xlwE6sb9h7feZ3xvoPtxt4A4gJUlwtgaXALvdzinv+IOBFn3LpwEHAdtr6nwGbsRLTP4DExooLuMi97W/dzz85F/YXcCPgADb6PPo3xP7y93nB6sK5wj0d6379u937o7PPuve619sBXBrkz3tNcX3q/n/g2T/v1vSeNlJcjwDfube/DOjhs+4s937cDdzcmHG5/34Q+ONp6zX0/noN6ywtB1b++glwO3C7e7kBnnPHvRmfM/iCvb/0SlGllGoiQrHLRSmllB+a0JVSqonQhK6UUk2EJnSllGoiNKErpVQToQk9RBlj2hpjXjfG7DHWiHzvG2O61bOuu40x24wxr7rPuf7UWKPSXWeMedEY06uada8wfka+q+V2mxtjflqfdauoryhYdVVR/2vuy91/0ZDbqStjjU55bc0lg77d0caYixp7u6pqjX6lqAqc+4Kbt4CFInK9e15/oA3Webp19VOsc6z3ui98iBIRz5C1i6tbUUTexc8FHrXU3L3tv9Rz/UZjjGkLXCQiHWssfGodz1W4wYwj6HUGYDRQBHx1luNQHsE8wV4fjfPAGjfjiyqWGeAxrItyNgPX+Sybg3XF3SbgIfe85wG7u+xcrAsc8rEuZCfrg9cAAAWRSURBVOmCz1DGWCPBrce6QGOpe95M3EOpYo2W+KZ7G2uAYe75DwIvu+v6HrjbPf91oNS9rcdOex2PAj/1+ftB4FdAItbFSuvdMV/pU6bI/Twa+K/P/GeBme7pgcDnwDqsq/LaueffjXX14CasERZP36+bfGIdAfTHGgFxE9aXawt3ueXAH9zb+NVpdWzG+hIzQB5wk3v+IuASrAuc5rvLbQDG+OzjN7DGUPnMvf6z7njfA97Hz3DDwAVYFyd9695fXaji81HDPssGHvLZ5z2wLqI7gnUhnWefTHXX+y1VfD710cC54WwHoI96vGlW8nmiimXXYA3DGYHVYt8PtMO6rPgF939oG/BfYKR7nWyglXv69P/Yy7Gu8EzFGhmuk3u+5+rPmZxK6P8Ehrunzwe2uacfxGrFxQCt3Mksyp0UqhpDOhP43Ofvre46I4Fm7nmtsL6APBfIVZvQ3dv8Ckh1z78O95jZWGNrxLinm/uJp1KsWIl8lHv6f4EnffbXX6p4Tc8Dk4A+WF94f3fP34X1RfUrYL57Xg/3exfrjv2Azz6/2uc9Pg84if+E/g1wlXs6Foin6s+H333m8/m4yz39U05d7fsg8GufdTYD7avah/po+Id2uTQ9w7FGwnMCOcaYz4HBWIPwj8dq+YGVQLoCX9Sy3qFYra69ACLib/znS4BeVo8QAM2MMUnu6ffEGlyq3BhzFCuZVElENhhjWhtjzsP6MjkhIvuNMVHAH4wxIwEX1nCjbbBaizXpjpVMP3HHGIF1yTZYCfpVY8zbwNvVVWKMScZKWJ+7Zy3EakF7VNVNtQLrfdgH/BWYbYxpDxwXkSJjzHDgGffr326M2Yd1IxBwj5vtnh7Jqff4kDHmM07j3u/tReQtd31l7vlVfT78jb7p69/u53VYXyj+fAksMMb8y6e8akSa0EPTd5waVOt0/obk9Mx/RET+Vs9tGmoe2tMGZIlIaaUVreRZ7jPLSe0+e0uwXmdbrO4ZgBuwEvxAEXEYY7KxWp++Kqh8wN+z3ADfiUiWn21NwkqUVwD3GWN6S/37qourmP8F1iiA52ONEXMV1utb4RNfbeus6b2o7nPgT1X7zMPz/lX53onI7caYC7H25UZjTH+pPDKjamB6lkto+gyIMcbc6plhjBlsjBmFlTSuM8ZEGGNSsZLUaqz+4lnGmER3+fbGmJpu1OFrFTDKPeogxpgUP2U+Bu70iamme4EWAknVLH8da1S9azl1Z6Vk4Kg7mY8B/B2k3If1SyHG3Zq+2D1/B5BqjMlyxxdljOltjLEBHURkGfAbrH7uxKqCEpF84IQxZoR71nSsPvNqicgPWN1EXUXke2Al1tDAnoT+BdYXFu4zls53x3y6/9/eHbpEFkVxHP/+UBD8H0xi8E8wGsxiMIpYjKJB4waLVRaEsWsTDQaTRRARBGFEMG1ZBFEEkQ2uCMdw7rA6zNOksJffJ755777hzp0zb869nHtEVoHsK5UPx3vc65EsB9zZ13ZA0iDN46Opzz7y7vOTNBwRpxHxgyxnPdR4pX0JP6H/hyIiJE0B62XJ4BOZ51wkv7Bj5MRUACsRcQPcSBoFTsoT8x+y4uFHtd3f3vNO0jywWwLgLTDRddoCsCGpTY6tI7LqXFOb95KOlZvrHkTEctfrlyV1cB3/aqRvA/uSzsjJuKse7f4uf/vbZH76vBx/Lsv7fpag1Q+skyuDtsoxkfMTD93tdpkFWiVI/gLmPjm/45RM9UAG8jUysEOu9mkpNyB5IXPYf9+ksDr2yInxi/Lem35MZoBNSatkJcDpcm2v8UGvPvvEPrAjaZLcw3NJ0gjZh4flHvaNXG3RzKwSTrmYmVXCAd3MrBIO6GZmlXBANzOrhAO6mVklHNDNzCrhgG5mVgkHdDOzSrwCElgUCbfax94AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import gaussian_kde\n",
    "density = gaussian_kde(adjusted.mean(axis=0))\n",
    "permutation_density = gaussian_kde(permuted.mean(axis=0))\n",
    "missing_z_density = gaussian_kde(unadjusted.mean(axis=0))\n",
    "xs = np.linspace(-1.,1.,1000)\n",
    "# density.covariance_factor = lambda : .55\n",
    "# density._compute_covariance()\n",
    "plt.plot(xs,density(xs), label='Fit model (not permuted)')\n",
    "plt.plot(xs, permutation_density(xs), label='Permuted model')\n",
    "plt.plot(xs, missing_z_density(xs), label='Model missing propensity score')\n",
    "plt.xlabel('Coefficient values for word counts')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 140,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(11778, 4684)\n"
     ]
    }
   ],
   "source": [
    "blacklist = ['theorem', 'theoretical', 'theoretically', 'theory']\n",
    "from lda_baseline.peerread_fit_topics import load_term_counts as ltc\n",
    "peerread_log_dir = '../logdir/peerread/modesimple/beta00.25.beta11.0.gamma0.0/split0'\n",
    "predictions_file = os.path.join(peerread_log_dir, 'predict', 'test_results.tsv')\n",
    "predict_df = pd.read_csv(predictions_file, delimiter='\\t')\n",
    "# predict_df = predict_df[predict_df.in_test==1.0]\n",
    "term_counts,vocab = ltc(None)\n",
    "term_counts = term_counts[:,~np.isin(vocab,blacklist)]\n",
    "term_counts = term_counts.toarray()\n",
    "print(term_counts.shape)\n",
    "\n",
    "# random_indices = np.random.choice(np.arange(4684), 4000)\n",
    "# term_counts = term_counts[:,random_indices]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 141,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model accuracy: 0.9960447119518486\n",
      "Mean and std. of the word coeffs: -0.012592992339901141 0.3258735530622863\n"
     ]
    }
   ],
   "source": [
    "n_bootstraps = 1\n",
    "n_w = term_counts.shape[1]\n",
    "adjusted = np.zeros((n_bootstraps, n_w))\n",
    "permuted = np.zeros((n_bootstraps, n_w))\n",
    "unadjusted = np.zeros((n_bootstraps, n_w))\n",
    "for i in range(n_bootstraps):\n",
    "    sample = assign_split(predict_df,num_splits=2)\n",
    "    sample = sample[sample.split==0]\n",
    "    indices = sample.index.values\n",
    "    labels = sample.treatment.values\n",
    "    words = term_counts[indices, :]\n",
    "    propensity_score = logit(sample.treatment_probability.values)\n",
    "    all_features = np.column_stack((propensity_score, words))\n",
    "    unadjusted[i,:] = fit_treatment(words, labels, coeff_offset=0, verbose=True)\n",
    "    adjusted[i,:] = fit_treatment(all_features, labels)\n",
    "    np.random.shuffle(words)\n",
    "    permuted_features = np.column_stack((propensity_score, words))\n",
    "    permuted[i,:] = fit_treatment(permuted_features, labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['involves', 'bounded', 'fundamental', 'feedforward', 'achieved',\n",
       "       'multiplication', 'guarantee', 'chosen', 'computable', 'prove'],\n",
       "      dtype='<U17')"
      ]
     },
     "execution_count": 142,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coeffs = unadjusted[0,:]\n",
    "top = np.argsort(abs(coeffs))[-10:]\n",
    "vocab[top]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAEKCAYAAAAyx7/DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xl8VNX5+PHPnez7HiAESIJAIDtJgACyKYsCigvuCFJBrEtrWwq2dfvWX5evfl1QWku1ghaBinXXVkGQXSDshM1AgBDISvZtlvP745KBQEK2STITnvfrNdxk5s69z4TkmTPnnvMcTSmFEEIIx2Ho7ACEEEK0jCRuIYRwMJK4hRDCwUjiFkIIByOJWwghHIwkbiGEcDCSuIUQwsFI4hZCCAcjiVsIIRyMc3N20jQtCygDzIBJKZVytf2Dg4NVREREm4MTQohrRXp6eoFSKqQ5+zYrcV8wVilV0JwdIyIi2LlzZwsOLYQQ1zZN0042d1/pKhFCCAfT3MStgG80TUvXNG1uewYkhBDi6prbVTJCKZWjaVoo8K2maYeVUhsu3eFCQp8L0Lt3bxuHKYQQok6zWtxKqZwL2zzgY2BIA/ssUUqlKKVSQkKa1b8uhBCiFZpM3JqmeWma5lP3NTABONDegQkhhGhYc7pKugEfa5pWt/8HSqn/tGtUQgghGtVk4lZKHQcSOiAWIYQQzSDDAUXXV3oWNi+C/COdHYkQNiGJW3RtZiP883b49hl4ZwKU53d2REK0mSRu0bXtXw15GTDud1BbARte6uyIhGgzSdyia9v/Ifj3get/BYNuhX0rwVjd2VEJ0SaSuEXXVVsBJ77XE7amQcK9UF0CWRs7OzIh2kQSt+i6zuwCiwkirte/jxgBTm6Qua5z4xKijSRxi67r9A/6NvxCFWIXD+iTBsclcQvHJolbdF2nt0PwAPAMvHhf1Fj9YmV5XufFJUQbSeIWXVfOLghPrX9frwtlds7s6vh4hLARSdyia6oogIp8CB1Y//4eCaAZ9KQuhIOSxC26prpZkiHR9e939dLvy9nd8TEJYSOSuEXXVFCXuAdc+VjYYL2rRKmOjUkIG5HELbqm/CPg6g1+4Vc+1iMeKgugPLfj4xLCBiRxi64p/zAE99Mn3lyurvsk/3DHxiSEjUjiFl1T0QkI7NvwY3WJO08St3BMkrhF12OxQGkO+PVs+HHvUPAIkBa3cFiSuEXXU5EPFiP49Wr4cU2DkIGSuIXDksQtup6SbH3r20iLG/TRJnmHZGSJcEiSuEXXU3ohcTfWVQL6xJzqYpn6LhySJG7R9ZSc0be+DQwFrBPcX98WyHJmwvFI4hZdT+kZcHavX1zqckEXRpwUneiYmISwIUncouspydYn3jQ0hruOb08wuEDR8Y6LSwgbkcQtup7SM1e/MAlgcIKACDgvLW7heCRxi66nrsXdlMBIaXELhySJW3QtZiOUnWu6xQ0QGKX3ccuQQOFgJHGLrqXsLKCuPhSwTkAk1JbrtbuFcCCSuEXX0pyhgHUCo/StdJcIByOJW3QtpRcSd3P7uEEuUAqHI4lbdC0lzZg1Wce/N6DB+az2jEgIm5PELbqWkmxw8wM3n6b3dXYD725Qcrr94xLChiRxi66l9EzzWtt1/MIvttKFcBCSuEXXUpLdvKGAdSRxCwfU7MStaZqTpmm7NU37oj0DEqJNSs8078JkHf9eeuKWsdzCgbSkxf0z4FB7BSJEmxmroLKwhV0lvcBULWO5hUNpVuLWNC0cmAy83b7hCNEGpTn6tjljuOvUrZIjFyiFA2lui/s14NeApR1jEaJt6pJvSy9OXvpcIRxAk4lb07QpQJ5SKr2J/eZqmrZT07Sd+fn5NgtQiGazzppsTeKWC5TCcTSnxT0CuEXTtCxgJTBO07R/Xr6TUmqJUipFKZUSEhJi4zCFaIbSViRujwBw9ZbELRxKk4lbKfW0UipcKRUB3AN8p5R6oN0jE6KlSrLBKwRc3Jv/HE3TW93Fp9ovLiFsTMZxi66jOQsoNMSvl7S4hUNpUeJWSq1XSk1pr2CEaJOSFo7hruMXLhcnhUORFrfoOlo6a7KOX099/Lex2vYxCdEOJHGLrqG6BGrLWtfi9umhb8vP2TYmIdqJJG7RNdQNBWzJGO46dYm79Kzt4hGiHUniFl1DaQtWvrlcXeIuk8QtHIMkbtE1tGQBhcv5dNe3kriFg5DELbqG0jOgGcC7e8uf6xEAzu6SuIXDkMQtuoaSbPAJAyfnlj9X0/RWd5lcnBSOQRK36BpKslvXTVLHp4dcnBQOQxK36BpaO2uyjk8P6SoRDkMSt3B8Sum1uNva4i47KyvhCIcgiVs4vspCfRWb1gwFrOPbA4yVUFNqu7iEaCeSuIXja8tQwDoyCUc4EEncwvFZE3cbWtwylls4EEncwvHVJe62dJXI7EnhQCRxC8dXchqcPcAruPXHsHaV5NgmJiHakSRu4fiKT4F/L30iTWu5eoKrD1TIeqnC/kniFo6v5LS+ik1beYdCeW7bjyNEO5PELRxfXYu7rby7QXle248jRDuTxC0cW22FPo7bJi3uEEncwiFI4haOrW5EiX/vth9LWtzCQUjiFo6t+MIiv7ZI3F6hUFMia08KuyeJWzi24pP61lYXJwEqpNUt7JskbuHYSk6DwfnizMe28O6mb8tlSKCwb5K4hWMrPq2XczU4tf1Y3iH6VoYECjsniVs4tpLTtunfhostbukqEXZOErdwbMWnbNO/DeBV1+KWxC3smyRu4bhqK/WiUIFRtjmesxu4+0tXibB7kriF4yo6rm+DbJS44cK0d2lxC/smiVs4rrrEbasWN8gkHOEQJHELx1WUqW9tmrhD5eKksHuSuIXjKjoOnsHg7me7Y3pJV4mwf5K4heMqPA5BfW17TO9QqC3Xi1cJYaeaTNyaprlrmrZd07S9mqYd1DTthY4ITIgmFWVCYDskbpBWt7BrzWlx1wDjlFIJQCIwSdO0Ye0blhBNqK2w7VDAOtZp75K4hf1ybmoHpZQCyi9863LhptozKCGaVPijvrV1V0ndJBy5QCnsmKbn5SZ20jQnIB24DlislFrQwD5zgbkAvXv3Tj558qSNQxVdgdFoJDs7m+rqNpZOrVtAwacHOLnYJjgAixlKz4BHILh52+64Qlzg7u5OeHg4Li71f281TUtXSqU05xhNtrgBlFJmIFHTNH/gY03TYpVSBy7bZwmwBCAlJUVa5KJB2dnZ+Pj4EBERgdaWxX1Lc6DcGXrEgmbDa+zKAmeNerXBupXfhbARpRSFhYVkZ2cTGRnZ6uO06DdeKVUMrAcmtfqM4ppWXV1NUFBQ25I2gLFKn6Juy6QN+vE0JzCbbHtcIQBN0wgKCmrzJ87mjCoJudDSRtM0D+BG4HCbziquaW1O2gCmanB2b/txGuLkDBZj+xxbXPNs8fvfnOZKD2Cdpmn7gB3At0qpL9p8ZiFay2IGcy24tC5xOzk5kZiYaL1lZWWxc+dOnnzySQDWb93Flm3bbRlxg9avX8+UKVNavc/u3bt5+OGHW3Xu4uJi/vKXv7TquR3hk08+ISMjo8XP8/bWr0vk5+czaVLX7RhoMnErpfYppZKUUvFKqVil1P90RGBCNMpUo29b2eL28PBgz5491ltERAQpKSksWrQIgPVbdrDlh3RbRdtu/vCHP/DEE0+06rntmbhNprZ3M7U2cdcJCQmhR48ebN68uc2x2COZOSkcj6lK3zp72OyQdS3brKws3lq6gleXLCMxMZGNGzfW2+/5559n5syZTJgwgYiICP7973/z61//mri4OCZNmoTRqHexrF27lqSkJOLi4pg9ezY1NfqbzX/+8x+io6MZOXIk//73v63HraioYPbs2aSmppKUlMSnn3561XjLysrYt28fCQkJ1rhmz57NmDFjiIqKsr4JAbzyyivExsYSGxvLa6+9BsDChQvJzMwkMTGR+fPn1zt2VlYW0dHRzJw5k/j4eO68804qKysBSE9PZ/To0SQnJzNx4kTOnj0LwJgxY/jNb37D6NGjef3115k1axaPPvooY8eOJSoqiu+//57Zs2czcOBAZs2aZT1XXQsZYPXq1cyaNYstW7bw2WefMX/+fBITE8nMzCQzM5NJkyaRnJzM9ddfz+HDem/tiRMnSEtLIzU1lWeeeabe65g2bRrLly+/6s/RYSmlbH5LTk5WQjQkIyOj7QcpPq1Uzh6lLJZWPd1gMKiEhASVkJCgpk2bppRSat26dWry5MlKKaWeW/hL9dIzP1fKbL7iuc8995waMWKEqq2tVXv27FEeHh7qq6++UkopNW3aNPXxxx+rqqoqFR4ero4cOaKUUmrGjBnq1Vdftd5/9OhRZbFY1PTp063nfPrpp9X777+vlFLq/Pnzql+/fqq8vLxeXJf67rvv1O23314vrrS0NFVdXa3y8/NVYGCgqq2tVTt37lSxsbGqvLxclZWVqUGDBqldu3apEydOqJiYmAZ/PidOnFCA2rRpk1JKqYceeki99NJLqra2VqWlpam8vDyllFIrV65UDz30kFJKqdGjR6tHH33UeoyZM2equ+++W1ksFvXJJ58oHx8ftW/fPmU2m9XgwYPV7t27lVJKeXl5WZ/z4YcfqpkzZ1qf/+GHH1ofGzdunDp69KhSSqlt27apsWPHKqWUmjp1qlq2bJlSSqk333yz3vGys7NVbGxsg6+xszX0dwDsVM3Msc0aDihEe3jh84Nk5JS2/InGCy1ul21XPDQozJfnpsZc9el1XSWN0gyAGSwmMLhe8fBNN92Ei4sLcXFxmM1ma19qXFwcWVlZHDlyhMjISPr37w/AzJkzWbx4MWPGjCEyMpJ+/foB8MADD7BkyRIAvvnmGz777DNefvllQB99c+rUqUZDPHv2LCEhIfXumzx5Mm5ubri5uREaGkpubi6bNm3itttuw8vLC4Dbb7+djRs3csstt1z1Z9SrVy9GjBhhjXPRokVMmjSJAwcOMH78eADMZjM9elwcMnn33XfXO8bUqVPRNI24uDi6detGXFwcADExMWRlZZGYmHjVGOqUl5ezZcsWpk+fbr2v7hPM5s2b+eijjwCYMWMGCxZcnGISGhpKTk5Os87haCRxC8ejzLaddHO5uiGGFiNwZeJ2c3MDwGAw4OLiYh0lYDAYMJlMqKtMamtsRIFSio8++ogBAwbUuz83t+HVeDw8PK4YUlYXF+gXYJuK5Wouj1PTNJRSxMTEsHXr1gafU/fmcHk8BoOhXmx1P6fLz9PYEDmLxYK/v3+jb7aN/Uyrq6vx8LBdd5o9kcQtOk1TLeMGGasg/7C+QLBnkO2DAnx8/SjNPaG3uFshOjqarKwsfvzxR6677jref/99Ro8eTXR0NCdOnCAzM5O+ffuyYsUK63MmTpzIG2+8wRtvvIGmaezevZukpKRGzzFw4ED+7//+r8lYRo0axaxZs1i4cCFKKT7++GPef/99fHx8KCsra/R5p06dYuvWraSlpbFixQpGjhzJgAEDyM/Pt95vNBo5evQoMTGt+H+8oFu3bhw6dIgBAwbw8ccf4+PjA1AvPl9fXyIjI/nwww+ZPn06Silr//6IESNYuXIlDzzwwBX92UePHiU2NrbVsdkzuTgpHIu1m8Sz3U4xdepUPv7POhKHjLzi4mRzuLu78+677zJ9+nTi4uIwGAzMmzcPd3d3lixZwuTJkxk5ciR9+vSxPueZZ57BaDQSHx9PbGzsFRfaLhcdHU1JSclVky/A4MGDmTVrFkOGDGHo0KE8/PDDJCUlERQUxIgRI4iNjb3i4iTobwzLli0jPj6eoqIiHn30UVxdXVm9ejULFiwgISGBxMREtmzZ0uKfz6X+9Kc/MWXKFMaNG1ev2+Wee+7hpZdeIikpiczMTJYvX84777xDQkICMTEx1ou3r7/+OosXLyY1NZWSkpJ6x163bh2TJ09uU3z2qlm1SloqJSVF7dy50+bHFY7v0KFDDBw4sPUHKMmGikLoEQ+2mMjTEIsZzu3Tp7z7dG+fc9jAq6++io+PT6vHcjcmKyuLKVOmcODAgaZ3tmOjRo3i008/JSAgoLNDuUJDfwctqVUiLW7hWIyV+sSb9kraAAYnvZ+7lV0lHeXRRx+t13csLsrPz+cXv/iFXSZtW5DELRyHUnpXSTt2k1gZnO0+cbu7uzNjxgybHzciIsLhW9shISFMmzats8NoN5K4heMw1ejV+1w7InG7gFnqlQj7JIlbOA6jPnsPlw4Y4uUALW5x7ZLELRyHsQrQ2q8q4KWcJHEL+yWJWzgOY6Xe2rZ1De6GGFz0xN0Oo66EaCtJ3MIx2PDCZF1Z19jYWKZPn24toFSP4cLctHZodb/22msNn/MqmlMCtq2ysrKanLDSnH1E+5PELRyDuUaf6m6D/u26WiUHDhzA1dWVt95668qdnBpO3Gazuc3nb03iFuJSkriFY2inGZPXX389P/6orxj/z3/+kyFDhpCYmMgjT/xST9JmI97e3jz77LMMHTqUrVu3EhERwW9+8xvS0tJISUlh165dTJw4kb59+1rfBC5vIT/++OMsXbqURYsWkZOTw9ixYxk7diygF5hKS0tj8ODBTJ8+nfLycqDxErCXWrp0KdOmTWPq1KlERkby5ptv8sorr5CUlMSwYcMoKioCYM+ePQwbNoz4+Hhuu+02zp8/D+hlWhMSEkhLS2Px4sXW45rNZubPn09qairx8fH87W9/s+nPXbSNJG7hGIyVgNbqVW8aYjKZ+Prrr4mLi+PQoUOsWrWKzZs3s2fPHpycXVj+76/BYqKiooLY2Fh++OEHRo4cCejV87Zu3cr111/PrFmzWL16Ndu2bePZZ5+96jmffPJJwsLCWLduHevWraOgoIAXX3yRNWvWsGvXLlJSUnjllVeorq5mzpw5fP7552zcuJFz5841eswDBw7wwQcfsH37dn7729/i6enJ7t27SUtL47333gPgwQcf5M9//jP79u0jLi6OF154AYCHHnqIRYsWXVE46p133sHPz48dO3awY8cO/v73v3PixIm2/LiFDUmRKdF5vl4I5/Y3b19jFaCabnF3j4Ob/nTVXaqqqqwlRa+//np+8pOfsGTJEtLT00lNTbXuE+o1BiwmnJycuOOOO+odo64salxcHOXl5fj4+ODj44O7uzvFxcXNe03Atm3byMjIsJZQra2tJS0tjcOHDzdaAvZyY8eOtZ7fz8+PqVOnWmPbt28fJSUlFBcXM3r0aEAvMzt9+vQr7p8xYwZff/01oH8K2LdvH6tXrwagpKSEY8eOWUvVis4liVs4AKX3bxts8+vaUD1upRQzZ87kj3/8Y90dcHYvWIy4u7vj5ORUb/+mSpY6OztjsVis9zdWslQpxfjx4+tVCgS9a6O5i8pefv5LY7vaMmJKqauWmX3jjTeYOHFivfuzsrKaFZNoX5K4RedpomVsZaqFvIPgFw5eIU3v3wo33HADt956K0899RShoaEUnT9P2Zk8+lzXuloXffr0ISMjg5qaGqqrq1m7dq21m6WuZGlwcDDDhg3jscces5aAraysJDs7+6olYFvKz8+PgIAANm7cyPXXX28tM+vv74+fnx+bNm1i5MiR9cqiTpw4kb/+9a+MGzcOFxcXjh49Ss+ePVsdg7AtSdzC/hltv8bk5QYNGsSLL77IhAkTsFgsuLi4sPjF+fSJat1wwF69enHXXXcRHx9Pv3796tXWnjt3LjfddBM9evRg3bp1LF26lHvvvde6qsuLL75I//79rSVgg4ODGTlyZJvqhyxbtox58+ZRWVlJVFQU7777LgDvvvsus2fPxtPTs17r+uGHHyYrK4vBgwejlCIkJIRPPvmk1ecXtiVlXUWHalVZ17KzUHYOusfrlfs6SmGmvgpOSHTHnVNcE6Ssq+j6aqvAya1jkzbofepmmfYu7I8kbmH/jJUdU8r1cnX1SmTau7AzkriFfTOb9O4K105Y9NXggnVEixB2RBK3sG+m9r8w2ai64YfSXSLsjCRuYd+MF8Y/d0QN7su1Y6EpIdpCErewb6Yq0JxsNvmmRZxc9K1FVsIR9kUSt7BvxmqbLw6saVq9tRpNJhMhISFXlk1tosUdERFBQUHBVc/VnH0as3PnTp588skWP+/mm29u0bR7R/bwww+TkZEBwB/+8IdOjqbjSOIW9kspMFXbfMUbLy8vDhw4QFWV3n/+7bffNjwrsJP7uFNSUli0aFGLn/fVV1/h7+/fDhE1zBalblvr7bffZtCgQUDnJ+6rlRewNUncwn5ZjPqIjna4MHnTTTfx5ZdfArBixQruvfde62NFRUVMmzaN+IQEhk2dyb59+wAoLCxkwoQJJCUl8cgjj3Dp5LV6JWEfeaTJZObt7c2CBQtITk7mxhtvZPv27YwZM4aoqCg+++wzoH5p2O+//57ExEQSExNJSkqirKyMs2fPMmrUKOuiEBs3bgQutvKzsrIYOHAgc+bMISYmhgkTJljfrHbs2EF8fDxpaWnMnz+/wcUR1q9fz6hRo7jtttsYNGgQ8+bNs9ZfubzU7dq1a0lKSiIuLo7Zs2dbZ4FGRESwYMEChgwZwpAhQ6wldPPz87njjjtITU0lNTWVzZs3A/D8888ze/Zs68+i7o2roqKCyZMnk5CQQGxsLKtWrQJgzJgx7Ny5k4ULF1qLh91///0888wzvP7669bX8tvf/vaKN8HGjrljxw6GDx9OQkICQ4YMoaysjOrqah566CHi4uJISkpi3bp1gF5Wd/r06UydOpUJEyYA8NJLL1nL4T733HNX/T1oNaWUzW/JyclKiIZkZGQ0f+eqEqXO7FKqutSmMXh5eam9e/eqO+64Q1VVVamEhAS1bt06NXnyZKWUUo8//rh6/vnnlVJKrV39D5UQO1AppdQTTzyhXnjhBaWUUl988YUCVH5+vsrIyFBTpkxRtbW1SimlHn30UbVs2TKllFJ9+vRR+fn5V8QAqK+++koppdS0adPU+PHjVW1trdqzZ49KSEhQSql6MU2ZMkVt2rRJKaVUWVmZMhqN6uWXX1YvvviiUkopk8mkSktL653zxIkTysnJSe3evVsppdT06dPV+++/r5RSKiYmRm3evFkppdSCBQtUTEzMFTGuW7dOubm5qczMTGUymdSNN96oPvzwQ2v8q1atUkopVVVVpcLDw9WRI0eUUkrNmDFDvfrqq9ZY6mJctmyZ9fXce++9auPGjUoppU6ePKmio6OVUko999xzKi0tTVVXV6v8/HwVGBioamtr1erVq9XDDz9sja24uFgppdTo0aPVjh07rP+vdU6cOKGSkpKUUkqZzWYVFRWlCgoK6r2+ho5ZU1OjIiMj1fbt25VSSpWUlFh/1rNmzVJKKXXo0CHVq1cvVVVVpd59913Vs2dPVVhYqJRS6r///a+aM2eOslgsymw2q8mTJ6vvv//+ip9tQ38HwE7VzBzb5BUfTdN6Ae8B3QELsEQp9frVnyVE0/68/c8cLjrc+A5mo77yjYtXs/u4owOjWTBkQZP7xcfHk5WVxYoVK7j55pvrPbZp0yY++ugjAMaNHklh0W8oKSlhw4YN1gUNJk+eTECAXoBq7dq1V5aEDQ296vldXV2ZNGkSoJdfdXNzw8XFhbi4uAYr8I0YMYJf/OIX3H///dx+++2Eh4eTmprK7NmzMRqNTJs2zVqq9lKRkZHW+5OTk8nKyqK4uJiysjKGDx8OwH333ccXX3zRYJxDhgwhKioKgHvvvZdNmzZx55131it1e+TIESIjI60lX2fOnMnixYv5+c9/bn1e3fapp54CYM2aNda+aYDS0lLKysqsP1s3Nzfc3NwIDQ0lNzeXuLg4fvWrX7FgwQKmTJnC9ddff9Wfb0REBEFBQezevZvc3FySkpIICgqqt09Dx9y/fz89evSw/l/6+voC+u/EE088AUB0dDR9+vTh6NGjAIwfP57AwEBAL4f7zTffWGvTlJeXc+zYMUaNGnXVeFuqOZfqTcAvlVK7NE3zAdI1TftWKZXR1BOFaBNlATSbXpi81C233MKvfvUr1q9fT2Fh4cXTXjpT8sIknLrypw2VQVWXl4RtBhcXF+uxmlOKdeHChUyePJmvvvqKYcOGsWbNGkaNGsWGDRv48ssvmTFjBvPnz+fBBx+s97xLS746OTlRVVVV//U14fLXW/f9paVumzrepceo+9pisbB161Y8PK7sBrs8ZpPJRP/+/UlPT+err77i6aefZsKECU0uWvHwww+zdOlSzp07x+zZs694vKFjTps2rdH/48Z4eXnV2+/pp5/mkUceuWpsbdVk4lZKnQXOXvi6TNO0Q0BPQBK3aJMmW8b5R/QV3YP7tcv5Z8+ejZ+fH3Fxcaxfv956/6hRo1i+fDnPPPMM67fsIDjAH19fX+v9v/vd7/j666+ty39dURK2qIiysjL69Oljs1gzMzOJi4sjLi6OrVu3cvjwYTw8POjZsydz5syhoqKCXbt2XZG4GxIQEICPjw/btm1j2LBhrFy5stF9t2/fzokTJ+jTpw+rVq1i7ty5V+wTHR1NVlaWtTRtXdnYOqtWrWLhwoWsWrWKtLQ0ACZMmMCbb77J/PnzAb3+eEOfGOrk5OQQGBjIAw88gLe3N0uXLr1iHxcXF4xGIy4u+jDO2267jWeffRaj0cgHH3zQrGMuXLiQnJwcduzYQWpqKmVlZXh4eFj/78eNG8fRo0c5deoUAwYMYNeuXfWOOXHiRJ555hnuv/9+vL29OXPmDC4uLk1+AmupFg2O1TQtAkgCfmjgsbnAXIDevXvbIDRxTVMKTDXg0bp62M0RHh7Oz372syvuf/7553nooYeIj4/H082FZa+9ABYzzz33HPfeey+DBw9m9OjR1t/zBkvCLl5s08T92muvsW7dOpycnBg0aBA33XQTK1eu5KWXXsLFxQVvb2/rMmXN8c477zBnzhy8vLwYM2YMfn5+De6XlpbGwoUL2b9/v/VC5eXc3d159913mT59OiaTidTUVObNm2d9vKamhqFDh2KxWKx1xRctWsRjjz1GfHw8JpOJUaNGNbxo8wX79+9n/vz5GAwGXFxc+Otf/3rFPnPnziU+Pp7BgwezfPlyXF1dGTt2LP6jY01jAAAgAElEQVT+/lcshNHYMV1dXVm1ahVPPPEEVVVVeHh4sGbNGn76058yb9484uLicHZ2ZunSpfU+GdSZMGEChw4dsr5BeXt7889//tPmibvZFxwBbyAduL2pfeXipGhMsy9Ommr0C5Nlee0bUFPKC/Q4jNWdG4eNlZWVWb/+4x//qJ588skr9rn04mhrNXZxtiOYzWaVkJCgjh492innv5q2Xpxs1nBATdNcgI+A5UqphpebFsKWTPpwMpyvbNV0KKeuOe39yy+/rDeM8He/+11nh2RTGRkZXHfdddxwww3WdTu7kiYXUtD0nvplQJFS6ufNOagspCAa0+yFFCoKoOQ0hA7q3ORdWwEFRyEgCjwa7k4QoqU6YiGFEcAMYJymaXsu3G5u6klCtImpBtDAybVz4zBIvRJhf5ozqmQT0D7jscQ1SV1ldXErU43e0m6noYDNJhUChY011cvRHDLlXXQod3d3CgsLm/7lNVd3fv82gMGgVyeUxC1sQClFYWEh7u5tq78jq7yLDhUeHk52djb5+fmN76QUlGSDmw/k1nRccI0pLQCnEvAq6+xIRBfg7u5OeHh4m44hiVt0KBcXFyIjI6++U9EJ+Nd0uOUNGNj0hJJ294+n9C6TWQ1PCxeio0lXibA/hZn6Nui6zo2jjlcIlOd1dhRCWEniFvanUC/9aTeJ2zsUKiRxC/shiVvYn6JMcPPVW7r2wCsEqs7r1QqFsAOSuIX9KfwRAqM6fyhgnbo3kIrWLUEmhK3JxUlhfwp/hPAhbT5MRY2J7VlF/JhbzpniKoxmC97uzkQGeTEkMpDIYK+mx5OD3lUCUJ4Lvj3aHJcQbSWJW9gXYzUUn4aE+1r3dLOFtYdyWbH9NFsyCzCa9fHiPm7OuDobKKs2UWvWl9/qG+LFzOERTE/uhYfrldXjrHzC9G3ZWaDx0qNCdBRJ3MK+nM8CVIsvTJrMFlanZ7No7TFySqrp7uvOQyMiGdUvhNievvh76lPnzRbF6aJKNh7LZ/WuMzz76UH+9v1xfjd5IDfFNdKa9r2QuEvPtP51CWFDkriFfbGOKIlq9lPWH8njfz7P4HhBBQm9/Hnh1ljGDgjB2enKSzhOBo2IYC8igr2YkRbBlswCfv/FIR5dvos7BofzP7fG4OV22Z+Fd6g+jrs0py2vTAibkcQt7Etd4g7s2/Su5TX8/osMPtmTQ1SIF0tmJDN+ULfm9VtfMLxvMJ8/PoJF3/3IG98dY192Me8+lEp4gOfFnQxO4NNDErewGzKqRNiXokzwDAYP/6vu9t3hXMa/uoEv95/lyRv68fXPrmdCTPcWJe06zk4GfjG+P//8yVDOlVZz+1+2cDCnpP5OvmHSVSLshiRuYV8Kj0NQ461to9nCH746xOylO+nm684XT1zPL8b3x835KhcXm2nEdcF89OhwnA0a9/39Bw6cuSR5+4ZJi1vYDUncwr4UZTbaTZJbWs1df9vKkg3HeWBYbz7+6XAGdPex6en7d/Nh1SNpeLs588A7P1xsefv21BO3DUpyCtFWkriF/agp14fcNXBh8sCZEm59czNHzpXx5n1JvDgtDneXtreyG9Ir0JMVc4bh6eLEzH9s51Rhpd7iNlZCdXG7nFOIlpDELexH0XF9e1mL+9uMXKa/tRVNg9XzhjMlPqzdQ+kd5Mn7Dw/FZFHMfHc7Za4XJuFId4mwA5K4hf0ourIq4NLNJ5j7/k76d/Pm08dGMCjMt8PC6RvizTszU8gpruL/bbrQZSKJW9gBSdzCftSVcw2MQinForXHeP7zDG4c2I2Vc9MI9W3bqiGtkdwnkNfvSWLDOX0Cj6U4u8NjEOJykriF/Sg6Dt7dUa5e/OGrQ7zy7VFuT+rJX+8ffPUp6e1sUmx35k0ZgVE5sSV9V6fFIUQdSdzCfhRmogKj+M3H+/n7xhM8mNaHl6cnNDgDsqM9OKIvZe7dOX/mKMu2ZHV2OOIa1/l/EUJcoIoy2VEWwIrtp3lsbF9euCUGg8FOSrsCAT37E+t5nhc+P8jaQ7mdHY64hkniFnZBVZegVeTzXZ4PT4y7jvkTo1s1C7I9aQERRBgKiO3px+Mf7K4/QUeIDiSJW3Q6pRRvf7IWgH4DE/jF+P6dHFEjAvqgVRXyzj0DCPRyZfbSHeQUV3V2VOIaJIlbdCqlFH/+zxH270sH4PYbR9ldS9sqIAKAENM53n0olapaM7OX7qCsWpY0Ex1LErfoVK+uOcZb32dySy+95aoFNr+ca4e7kLg5f5L+3Xz46wPJ/JhXzmMf7MZ4YXEGITqCJG7RaRav+5FFa49xV0o4N4SW6fVAXD2bfmJn8e+jb89nATCyXzB/uC2ODUfzefbTgyipYyI6iCRu0Sn+vuE4L/33CLcn9eSPt8ejFR3XFwi2Zx4B4OYH509Y77ortRePje3Liu2n+NuG450YnLiWSOIWHe4fm07w/746xOT4HvzvnfE4aegLKFylnKtd0DS9AFbdYg8X/HL8AKbE9+BPXx/m37tkZqVof7ICjuhQ723N4n++yGBSTHdeuztRn1xTnqdX3QuJ7uzwmhYyEI6vq3eXwaDx8vQEzlfW8qsP9+Jk0Lg1sWcnBSiuBdLiFh1m+Q8nefbTg4wf1I1F9ybhUjcjMv+wvg0Z0HnBNVfIAL30bFX98q7uLk68/WAqQyIDeWrVHj7fK8WoRPtpMnFrmvYPTdPyNE070BEBia5p1Y5T/PbjA4yLDuXN+5Jwdb7kVy//iL51hBZ36EB9WxfzJTxcnXhnZiopfQL5+ao9fLJbljoT7aM5Le6lwKR2jkN0YavTs1n47/2M7h/CX+4ffOUyY/lHwM1XX5DX3tV9Ksg/1ODDXm7O/OOhVFL6BPDzVXt4Z9OJBvcToi2aTNxKqQ1AUQfEIrqgj3dnM3/1Xkb0DeZvM5IbXrUm/zAE99cv/tk7v97g4tlgi7uOt5szy2YPYVJMd37/RQZ/+vqwDBUUNiV93KLdfLY3h1/+ay/DIoP4+4MpjS81ln/EMbpJAAwG/U2mrl++Ee4uTiy+fzD3D+3NW99n8sSK3VTWmjooSNHV2WxUiaZpc4G5AL1797bVYYWD+nLfWZ5atYeUiEDemZXSeD3tyiKoyLPphckKYwXnq89TUlPC+Rp9W2OuodZci9FipNZci6ZpuBpccTG44OrkiruzO/5u/gS6B1q3Lk4uDZ8gdBD8uEZfOPgqnxKcDBovToslPMCT//3vYY7nV7DkwWTCA+x4kpFwCDZL3EqpJcASgJSUFPlceA37z4Fz/GzlbpJ6+fOPWal4ul7l16zgqL5tYYu7rLaMY+ePcez8MU6VneJM+Rnrray2rA3R6zQ0unt1J9wnnHDvcCL8IogOiCY6KJrAnoNh7wdQkg3+va5+HE3j0TF9ie7uw5Mrd3PLm5tZfN9g0voGtTlGce2ScdzCptZk5PLEil3Ehfvx7kOpeLs18StmHQrYeEXA8tpy9hXsY2/+Xg4WHOTo+aOcrThrfdzdyZ0w7zDCvMNICEmgh1cPa8vZ390fPzc/PJw8cHFywcWg3wBr69toMVJprKS4ppjimmLO15wnrzKPM2VnOF12mo1nNvLxjx9bzxfqFkBiSBCpu//CkKQ5RPpFNlkYa2x0KJ88NoI57+3k/re38dSN/fnp2OtwsqN648JxNJm4NU1bAYwBgjVNywaeU0q9096BCcez7nAeP12+i0E9fFk2ewg+7o10NVzq3H5w9dEv+l1QWFXItrPbSM9NZ2/+Xo6dP4ZCoaER5RdFYmgidwXcRf+A/vTz70d3r+7tXlGwpKaEw0WHOVx0mIyCA6RXFPDNqS/g1BcEewQzOnw043qPY2iPobg5uTV4jL4h+oLHv/vkAP/37VG2Hi/ktbsTO2UtTeHYtPa42p2SkqJ27txp8+MK+/X90XzmvKevxr78J8Pw82xG0gZ4ezzVBid2TXyWrWe3sjVnK0fO6yM2vF28iQ+JJzEkkYSQBOJC4vBx9WnHV9F86u0byDYY2Hn9T9mSs4WNZzZSYazAw9mDUeGjuKXvLQwPG46z4cq2kVKKD3dm8+xnB/B2c+aVuxIZ1T+kE16FsCeapqUrpVKata8kbtFWm38sYPbSHUSFeLNizlD8PV2vur9FWThcdJitZ7awdetL7HZ3pxYLzgZnkkKTSOuRRlpYGgMDB+Jk6LxFgq/q6wWw6z1YeBqcnKk117Lj3A7WnlrLtye/pbimmGCPYKZGTWVav2lE+V1ZQOtobhmPf7CLo7nlPDqmL78Y3//ibFJxzZHELTrM1sxCHlq6nYggLz6YM4xAr4aT9tnys9YW9baz2yiu0aeM96utZVjYCNISHiK5WzKeLg4y4uLAv2H1Q/CTNdArtd5DRrORDdkb+CTzEzZlb8KkTAwPG84DAx9gRM8RGLSLybmq1sz/fHGQFdtPM7i3P4vuTZJRJ9coSdyiQ2w/UcTMf2wnPMCDFXOHEex9sW+3wljBjnM72JKzha05W8kqzQIgxCOEtLA0hvUYRlpJIcGf/QzmbYLucZ30Klqpsgj+NwrGPA1jFjS6W2FVIR8d+4iVh1eSX5VPhG8E9w+8n2nXTcPd+WLf9ud7c/jNv/ejafC/d8YzKdYBZpEKm5LELdpd+skiZryzne5+7qycO4xAL2cOFh5ka47eqt6Xvw+TMuHu5E5K9xTSeqQxPGw4ff37XryQ+M0z8MNb8PQZcL5694pdWjIGnN1h9n+a3NVoNvLNyW9Yfmg5+wv2E+QexKyYWdw14C7rp4xThZU8sWIXe7NLmDGsD7+dPLDxSUuiy5HELdrV7lPnmfHOFvwDcrlzRC1HSvayJ28P5cZyNDQGBg20JurE0ERcnRpJykunQE0ZPPJ9x74AW1nzAmx+HRZkgbtvs5+249wO/rbvb/xw9gf83fx5cNCD3BN9Dz6uPtSaLLz8zRGWbDhOdHcf3rxvMNeFerffaxB2QxK3sLlKYyX7Cvbx1bHNfJyxEdxPgaYvktvXry/J3ZJJ7Z7K0B5DCXAPaPqAZiP8sRckz4Kb/tS+wbeXExth2RS46z0YdGuLn74nbw9L9i1h45mN+Lj68JPYn3DfwPvwcPZg3eE8fvnhXr0P/NYY7kwOt99FlIVNSOIWbVZaW8qevD3szN1Jem46GQUZmJQJlIbB1JNbB4xkdO9hJHVLItA9sOUnyE6Ht8fB9KUQc5vN4+8QZhP8X3+IGgN3/qPVhzlYeJDFuxez8cxGQjxCmJcwj9v63UZRuZmfrdzNtuNF3JkczovTYqXrpAtrSeKWmZMC0C+i7crbRXpuOum56RwpOoJC4WxwJi44jsl97uHzH9zxUtexas5YegW2ceTD6W36ttewtgffWZycYeBU2PchGKvAxaNVh4kJiuEvN/6F9Nx0Xt/1Or/f9nuWHVzG40mP8/5PJvDmd5m8vvYYh8+V8tYDUutESIv7mnWu4py1NZ2em86JEr1utLuTOwmhCSR3SyalWwqxwbGcLDBy399/wM3ZwMq5w+gT5NX2AFbNgLN74Of7236szpS5Dt6fBne9D4NuafPhlFJsyN7A67tf59j5Y0QHRvNk0pPUlPbjqX/txdmg8ca9gxnZL9gGwQt7Il0loh6lFKfKTlmTdHpuOmfK9dVZfFx8SOqWRHK3ZJK7JTMocFC9qngZOaXc//Y2XJ0NrJybRmSwDZK2xQIv94PrboDbl7T9eJ3JbIJXBkKvIXDPcpsd1qIsfHXiKxbvXkx2eTYp3VK4q+88Xv2imh/zyvn1pGgeGRUl/d5diHSVXOOMZiMZRRnsyduj3/L3UFBVAECgeyDJ3ZKZMWgGyd2S6effr9HZiQfOlPDAOz/g4eLEijnDiLBF0ga9pV1ZAH1vsM3xOpOTMyTcA1sXQ1ku+HSzyWENmoEpUVOY2Gciq4+t5q29b/HrLQ8zOn4cPc+N509fH2ZfdjH/e2dC04W8RJcj/+NdQHF1MXvy97A7bzd78vZwsPAgNeYaAHp692Roj6HWFnWkb9OV7AD2Z+tJ29vNmRVzhtE7yIb9qj+uATS9xd0VDH4QtizSS72OfMqmh3ZxcuHe6Hu5te+tvJfxHksPLqXKvJ6hqeP4754Uji0u568PJMuQwWuMdJU4mEpjJUfOHyGjMIODBQfZX7DfOivR2eDMoMBBJIYm6reQREI8W168aO/pYma88wM+7i6snDus7RciL/f3cfoiBHPX2fa4nekfN0F5Ljy+U18lp50UVRfx931/Z9WRVSgFlpIRmArH8vIdadwUJ7MtHZn0cXcRlyfpjMIMTpSewKIsgD59PCYohoTQBJJCk4gJiqk3jbo1dmQVMXvpDvw9XVgxZ5jtRzAUnYBFiXDjCzDy57Y9dmfauwo+ngv3fwT9bmz30+WU57B4z2I+z/wcTblTlT+KBwbdz28mJeAshaockiRuB3Rpkq67HS85bk3SwR7BxATFMChokPUW6hlq0xjWZOTy2Ae76OnvwfsPD6Wnf+uGt13Vhpfguxfh5weaXD3GoZhq4bU46DYIZnzc9P42cuz8MV5Nf52NZ77HYvIhTN3CsumPE+YnXSeORhK3nasyVXGk6EJLuvDgFUk6yD2ImGA9Sdcla1sn6cutTs9mwUf7iAnz5d1ZqQR5N7wYQJtYzLAoCfzC4aGvbH/8zrbhZfju9/DoVj2Bd6Bdubt4ZsOfOVWZgWYM4ZH4n/Jo6u31KhEK+yaJ244opThZepLdebvZnbebg4UHySzOxKzMgD7K49KWdExQDKGeoR02zEspxV/WZ/LSf48w4rog/jYjpf1GKRz6HFY90Oop4navsgheGQRxd8Ctizv89Eoplu//Dy/veBWz81kCnHvz62E/5abISfZb11xYSeLuRGaLmUNFh0jPTbcm66LqIgD83fyJCY6xJuqYoBi6eXbrtLG41UYzCz/axyd7crglIYyXpsfj5txOf+AWMywZDdWl8ORu6KqJ5MtfQfpSeGInBER0SghFFdU8+sk77K/4N05uefT06s28xDlMjppsXW9T2B9J3B0styKXLTlb9NrTZ7dSUlMCQLh3OIO7DWZw6GCSuiU1eyheR8grrWbu++nsOV3ML8f35/Fx17VvbDvfhS9+rtf0iL2j/c7T2Upz9O6gmNvhtr92aiif7MnmmW9WoPzXoLnlEOYVxk/ifsItfW9p80VsYXuSuNuZUorDRYf55uQ3rD+9nh+LfwT0C4jDw4aTFpbGkO5D2r1furU2HsvnqVV7qKgx8+rdCe1ftD//qN7a7pkMMz8HO3nzajf//S1s+ws8ugVCB3ZqKGdLqvjVh3vZdm4TAWEbqHE6gb+bP9P7T+ee6Hvs9nf0WiSJux0opThUdIhvsr7h25PfcqrsFE6aEyndUhjRcwTDw4bTP6C/3bSoG2I0W3j126P89ftMrgvx5s37BjOgezsvvlueD+/eBFVFMG8z+F4DY40rCuGNJOgWCzO/aNdx3c2hlOLTPTn8/suDlFiOENl3F7mmdJw0JyZGTmTGwBnEBMd0aoxCErfNKKXIKMrgm6xv+CbrG7LLs3HSnBjSfQgTIiZwQ+8bmld72g7szy5h/uq9HD5Xxj2pvXhuagweru3cz1xZBMtugcIf9SFyfdLa93z2ZNf78NnjcPPLMGROZ0cDQEmVkVe+OcLyH07h4nae2EH7OVW7nkpTBTFBMdze73ZujrwZb1cZStgZJHG3gVKKjMIM/nvyv3yb9a01WQ/tMZQJfSYwrvc4h0nWAOU1Jt5Ye4y/bzxOsLcbv58Wy8SY7u1/4qpieO9WyMuAe1d2nentzaUULL8Tjn8Ps76E3kM7OyKrrIIKXl1zlM/25uDpamRwbCbFTpvIKvsRD2cPxvcZz+39bmdw6GC7/gTZ1UjibiGlFAcLD+ot65PfcKb8DM6aM0PDhjKxz0TG9hqLv7t/Z4fZIiazhVU7T/Pqt0cpKK/lntRePH3zQPw8OmBUQXUpvH8bnN0Ld/8TBkxq/3Pao6rz+vT+qvP6J46wpM6OqJ7D50pZ8v1xPtubg0VZGBpdiXfwLvYVr6PSVEkPrx5MjJjIpIhJDAoaJEm8nUnibgaLsrC/YD/fZn3Ltye/JaciB2fNmWFhw6wtaz83v84Os8VqTRY+2X2Gt77P5HhBBUMiAvnt5IEk9OqgNx6l9LHaR77Wx2sPnNIx57VXRSfgvVv0TyB3v6+vlmNnzpVU897WLP61M5uC8hoCvBWJA05jdN/NweIdmCwmevn0YkKfCYzpNYa44DgZF94OJHE3wmwxsytvF2tOrmHNqTXkVebhbHAmrUcaEyImMLbXWIdM1qAP7/swPZt/bjvJ2ZJqYsJ8efKGfkwY1MHjxLe8Ad/8Dib+AdIe67jz2rOSbFg+HQqOwi1vQOJ9nR1Rg0xmCxuO5fNR+hm+zcil1mwh2NdEdN+T1Ljt4ljpHszKjJ+bHyPCRjAqfBQjwkY43KdReyWJ+xJVpip2nNvBd6e+Y93pdRRVF+Hm5MbwsOGM7zOe0b1G4+va/BW67UlFjYkNR/P5ePcZ1h7Ow2xRpEUFMW9MX0b1C+74j7Ynt+grt0ffrK8IIx+tL6ougX89CMfXw+iFMGahXf98SqqMrD+Sx38PnmP9kXwqa814e9TSLyIHN58jnK7ZTWltMRoa/QL6WVdMSu6WTJBHUGeH75Cu+cR9svQkm85sYuOZjew8t5Macw2ezp6MCh/FDX1uYFTPUXi6ON66fUopsgor2Xa8kLWHctlwrIBak4Vgb1fuSA7n7pReRIV00oiA8jx463pw9YS568HdMT+5tCtTrT4Jac9ySLgXpi4CZ9fOjqpJ1UYzWzILWHMojw1H88k+XwVYCAstoGfYacyumZyuzKDaXA1ApF8kccFx1tnB0YHRMuGnGa65xH22/Kx1/cTt57Zzuuw0ABG+EYzsOZKRPUeS0j0FN6d2KJzUjsqqjRw6W8bBnBJ2nyrmhxOF5JZeWCDB34MJMd2YMKg7qREBnVvK02zS113M3gEPr4XusZ0Xi71TSq+QuO7/QeQo/ZOJh+N0NdQ1HjYey2fD0Xy2ZhZSUWtG00xEhhUTEnoGs+tx8mqOcb5GL/XgpDlxnf91DAgcQF//vlznfx19/fvSw6uHFMG6RJdO3DXmGo4WHeVgob6IwM5zO8mpyAHAx9WH5NBkhvcczsieI+nl4xhlQ6uNZrIKK8jMq+B4fjmHz+nJOquw0rpPqI8bQ6OCGBoZyLCoIPqGeNnPVf41L8CmV+DWv0DS/Z0djWPYuxI+fRyC+sL9H4J/786OqFVqTRZ2nzrP9hNFbM8qYtfJ81TUmgFFz+BawrsV4u6dQ5V2kryaExRWF1if6+nsSV//vvT170sf3z708ulFb5/e9PLpdU2OJe8SiVspRVF1EcdLjnO8+DiHig6RUZjBsfPHMCkTcHH9xLr+tX4B/ez2HVwpRX55jZ6cC8ovbvPLyT5fxaX/Db0DPYkJ82VQD19ievoSE+ZHqI+b/STqSx38GD6cBYNnwi2LOjsax3L8e321e2c3fax7eHJnR9RmJrOFjLOlbD9RRPrJ8xzMKeVU0cUGSLCvme7Bxfj4FKG55VLFGQpqTlFcW1jvOIHugRcTuW8vwr3DCfMOo4dXD0I8Q7pksSyHTNwmi4mVh1eSWZLJ8eLjHC85TnFNsfVxX1dfYoJirNX1YoJi6O7V3e6SWUmVkVOFlZwqqtRb0fnlZObrLemyapN1Pw8XJ6JCvIgK8aZviBd9Q7yJCvEiMtgLT1cHWQr09HZ9kk23WJj1hZ6ARMvkHYYP7tKXPbvlDYi/q7MjsrmSKiMZOaUczCnh8LkysgoqyCqsoKC89uJOWg3unucJ9CvD06sYZ7dCzE75VFjyKDPl1zueQTMQ4hFCmHcY3b2608OrB2FeYfTw7kE3z24EeQQR4BbgcEMWbZ64NU2bBLwOOAFvK6X+dLX9W5O4lVKMWDkCg2agr19fovyjiPKLsn7dmeVPL1VRY+JsSTW5pdWcOV/FyaIKTl5I1KeKKimuNNbbv7uvO31DLyTmYC/6hnoTFeJND193DIbOfz2tdvx7WHk/eAXD7P+ATwfMxuyqKgrhXzPg5GYYNE2fJu/d8rVCHU1ptZGTBfrfzbnSas6VVHG2pJpzJdXWvzGTRYFmxOByHs2lGHf3Uny8y3FzL0VzOU8tRVRYCrEoU71jGzQDge6BBHsEE+QRRLB7sPVrPzc/fF198XH1wcfVx/q1p7Nnp+YYmyZuTdOcgKPAeCAb2AHcq5TKaOw5re0qKa0txcfFp0N/eEopKmvNFFXUUlxppKiylvMVtZyvrKWwvJZzpfov0NmSanJLqimrqf8L4mzQ6BngQe9AT3oHetInyJPegV7690Ge7bcoQXuxWKCmFFy9wamB2CsKYPNrsHUxBF0HMz4Bv54dH2dXYzbBltdh/Z/A4AJDHoah88A3rHnPL8uFrI1wYgNk74SS02CqBhdPCO4H3eOg11D9FhBh10MR61gsioKKGnKKqzlbXMWZ4ipyiqvJKa7ibEkVZ4qrKSivASxozuVozsUYXErw8arG26sKd7dyDM7lmA2lVFuKKTcVY74swV/KoBnwcvHC3ckdNyc33J3dcXVytX7v5uyGq8EVJ4MTTtqF22Vf+7j68Fhi6+Yv2DpxpwHPK6UmXvj+aQCl1B8be05rE3e10Wzt61XoX1z8/sJWKYxmRY3JTI3RQq3ZcmGrf19jtlBZY6a8xkh5jZnyatOFr02UVZuoqDFZv65L1LUmS4PxGDQI9XGnu5873X0vbC/5OszPgzB/d8danFUpfQp22Tn9j7vouD677/wJfVt8EswXPsK6+ugtaq9gcPOFijzIPagfI+l+mPRncLv2LiK1q4Jj8P2fYf9q/fs+IyBqtD5d3jdM/38w10JloV68K3sHZG3SJ/cAuPlBryEQGAUuHlBTBvlH4Nw+/Q0ZwLu7Xjuld5qeyPnmIK4AAA6wSURBVEOi9X0dIJlfrtpo5lxJ9YWkfjGx55RUWe+rNtb9fSswVOHmVkOAtwkvDyOe7rW4utbi4lKDk3M1SqtCaSYUtSiMmKnFoowYVQ1GSy0mixELZizKgkWZ9RsXv/Zx8Wft3f9p1WtpSeJuTnOwJ3D6ku+zgXapmJP5YjKuSk8aGhffUC79+nLOKJwBrwb2rfs1NGigaRpOmv67adA0DJp+v8FLv9+gaTgZLrn/wj4aQC1QcOFW743usrgae+yKN8fWPNbMczX1mLHqYmKu4+oNAZEQGg0DbgLvUKit0BN8RQFUFuhlWb1CYdTN+kIIIQMQ7SC4H9zxNoz9Dez7EDI+hXV/4Ir/xzquPnrVxcT79eGFPRIaXl3IYoa8Q3B6G5zaBqd+0I9dx8kV3P316xSaBpoTaIaLNzvlDkRcuF3BDVQ3sCiF0WzBZFaYLPrX5hqFuVphsSjMFoVZKSwNt99apKyDflTNSdwNvQ1f8VukadpcYC5A796tG9rk2n0gmsV48YT1znzxG4NBw6AZMBguJFiDAae6ZGzQcHYy4Gww4Oyk4WwwYLi0JVGvVdHY/S157PJX0YxztcdjzY3f2U1vcfl0A79eesL2CnbI1laXFhgFYxbot+oS/ZNOea5ewMvZXR/7HdhX7/ZoqEvrcgYnfXx991hIfVi/r+SMnsjPn4TqYr2eitkIynLJzdxA48JxaOgX5ppzmdKC3j2jJ3GFWenfW5Sq9yNQl/x7ucAOmoVtV10lQghxrWpJV0lzGvY7gH6apkVqmuYK3AN81pYAhRBCtF6Tn7GUUiZN0x4H/ov+ieMfSqmD7R6ZEEKIBjVrrJpS6ivgq3aORQghRDPY7+ViIYQQDZLELYQQDkYStxBCOBhJ3EII4WAkcQshhINpl7KumqblAydb+fRg9Mnl9kbiahmJq2UkrpbpinH1UUo1qyxkuyTuttA0bWdzZw91JImrZSSulpG4WuZaj0u6SoQQwsFI4hZCCAdjj4l7SWcH0AiJq2UkrpaRuFrmmo7L7vq4hRBCXJ09triFEEJcRackbk3TpmuadlDTNIumaY1egdU0bZKmaUc0Tfv/7Z15rFTVHcc/X5Y+YlABQVkEBYsLtgmuFRcWpdpqAlKx0FQFtTVIK6l1j22DJtVa/9BYtbjEpWpFRbG4xSKrVcS68RAR2VwQ8BkVLK0+lf76x/lNuQz3zgzvzZ33XnI+yc2c+Z3tN79z5nfvnLnnd1dJujwh7y9psaSVkh7ycLPV0KubpNne7mxJXVPKjJD0RuL4UtKpnnePpLWJvMG10svLbU30PSshb0l7DZa0yMe7XtK4RF5V7ZU1XxL5df75V7k99k3kXeHyFZJOao4eTdDr15LecvvMkbRPIi91TGuk10RJHyf6/1kib4KP+0pJE2qs1w0Jnd6RtCmRl4u9JN0lqUHSmxn5knST61wv6dBEXvVtZWY1P4CDgAOA+cDhGWXaA6uBAcC3gCXAIM97GBjv6WnA+VXS64/A5Z6+HLiuTPluwKfALv7+HmBsDvaqSC9gS4a8xewF7A8M9HRvYAPQpdr2KjVfEmUmA9M8PR54yNODvHwd0N/baV9DvUYk5tD5Bb1KjWmN9JoI3JxStxuwxl+7erprrfQqKn8BIdR03vYaChwKvJmRfzLwDOGhO0cBi/O0VYtccZvZcjNbUabYkcAqM1tjZl8B04HRkgQcD/jTVLkXOLVKqo329iptdyzwjJn9p0r9Z7Gzev2flraXmb1jZis9vR5oACraZLCTpM6XEvrOAE5w+4wGpptZo5mtBVZ5ezXRy8zmJebQS8DeVeq7WXqV4CRgtpl9amafAbOBH7SQXj8BHqxS35mY2ULCRVoWo4G/WOAloIukXuRkq9a8xp32kOI+wB7AJjP7pkheDfYysw0A/rpnmfLj2XHS/N5/Kt0gqa7GenWS9IqklwrLN7Qie0k6knAVtTohrpa9suZLahm3x2aCfSqpm6deSc4lXLkVSBvTWup1mo/PDEl9d7JunnrhS0r9gbkJcV72KkeW3rnYqqIHKTQFSc8BPVOyrjSzv6XId2giRWYl5M3Wq9I2vJ1ewHcJTwYqcAWwkeCcbgcuA66uoV79zGy9pAHAXElLgc9TyrWUve4DJphZ4XnaTbZXWhcpsuLPmcucKkPFbUs6AzgcGJYQ7zCmZrY6rX4Oej0BPGhmjZImEX6tHF9h3Tz1KjAemGFmWxOyvOxVjprOrdwct5mNbGYT64C+ifd7A+sJcQC6SOrgV00FebP1kvSRpF5mtsEdTUOJpn4MzDSzrxNtb/Bko6S7gYtrqZcvRWBmayTNBw4BHqWF7SVpN+Ap4Df+M7LQdpPtlULWfEkrs05SB2B3ws/fSurmqReSRhJOhsPMrLEgzxjTajiisnqZ2SeJt3cA1yXqDi+qO78KOlWkV4LxwC+SghztVY4svXOxVWteKkl9SLGFFf95hPVlgAlAJVfwlTDL26uk3R3W1tx5FdaVTwVS/4HOQy9JXQtLDZK6A8cAb7W0vXzsZhLW/x4pyqumvSp5qHVS37HAXLfPLGC8wl0n/YGBwMvN0GWn9JJ0CHAbMMrMGhLy1DGtoV69Em9HAcs9/SxwouvXFTiR7X955qqX63YA4c++RQlZnvYqxyzgLL+75Chgs1+Y5GOrPP6BLXcAYwhnokbgI+BZl/cGnk6UOxl4h3DGvDIhH0D4Yq0CHgHqqqTXHsAcYKW/dnP54cCdiXL7Ah8C7YrqzwWWEhzQ/UDnWukFHO19L/HXc1uDvYAzgK+BNxLH4DzslTZfCEsvozzdyT//KrfHgETdK73eCuCHVZ7v5fR6zr8HBfvMKjemNdLrWmCZ9z8PODBR9xy34yrg7Frq5e+nAn8oqpebvQgXaRt8Lq8j/BcxCZjk+QJucZ2XkrhbLg9bxZ2TkUgk0sZozUslkUgkEkkhOu5IJBJpY0THHYlEIm2M6LgjkUikjREddyQSibQxouNu5UjqKWm6pNUKEeSelrR/E9uaImm5pAf8nuXnFKKojZN0p6RBJeqOUkqktgr77SJpclPqZrS3pVptZbT/oG/zvjDPfnYWhWiKY8uXrHq/wyUdXet+I9nktnMy0nx8Y8pM4F4zG++ywcBehPtcd5bJhHuU1/omgY5mVgil+lCpimY2i5SNEBXSxfu+tYn1a4aknsDRZrZP2cLb6hR2pVZTj6q32QyGA1uAF1tYj0iBat44H4/qHoS4EAsz8gRcT9i8shQYl8i7hLADrR64ymXTgK+87GWEzQCbCRs+9iMRYpcQvew1wkaGOS6biIf4JET3e9T7+CdwjMunAnd5W2uAKS6fDnzhfV1f9DmuAyYn3k8FLgI6Ezb1vOY6j06U2eKvw4EnE/KbgYmePgxYALxK2KnWy+VTCLvp6gkRAYvtWp/Q9ThgMCFiXz3hJNrVy80HrvE+LipqYynhZCXgE+Asl98HjCRsBLrby70OjEjY+BFCjJC5Xv9m1/cp4GlSwuAC3yZs4lni9tqPjPlRxmbvAlclbH4gYbPZRsKGs4JNTvd2l5AxP+ORs29oaQXiUWJwgpO5ISPvNEKIyPaEK/D3gV6ELbW3+xe3HfAkMNTrvAt093TxF3g+YcdjD0I0s/4uL+yGnMg2x/1X4FhP9wOWe3oq4aqsDujuTqujf/mz4hgfAixIvH/L2+wA7Oay7oQTTWHDWEnH7X2+CPRw+Tg8ZjMhfkSdp7uk6LOdrgSHPczTVwM3Jux1a8ZnmgacAnyHcGK7w+UrCSeki4C7XXagj10n131dwuY/Soxxb2AT6Y57MTDG052AXcieH6k2S8yPCzw9mW27X6cCFyfqLAX6ZNkwHvkfcamk7XIsIXLbVuAjSQuAIwgB308kXMlBcBQDgYUVtnsU4SpqLYCZpcUgHgkMCis5AOwmaVdPP2UhSFKjpAaC08jEzF6XtKek3oSTxmdm9r6kjsA1koYC/yWEwtyLcPVXjgMITnO269iesF0ZgiN+QNLjwOOlGpG0O8ExLXDRvYQr4gJZy0vPE8bhPeDPwHmS+gCfmtkWSccCf/LP/7ak9wgPnACP3ezpoWwb4/WS5lKE272Pmc309r50edb8SIsWmeQxf32VcOJI4wXgHkkPJ8pHakh03K2bZWwLDlVMWrjIgvxaM7utiX2K8mEn2wFDzOyL7SoGJ9mYEG2lsjk2g/A5exKWVQB+SnDkh5nZ15LeJVxNJvmG7f9gL+QLWGZmQ1L6OoXgEEcBv5V0sDV9LfnfGfKFhKh1/QgxUMYQPt/zCf0qbbPcWJSaB2lk2axAYfwyx87MJkn6HsGWb0gabNtHEozkTLyrpHUzF6iT9POCQNIRkoYRnMM4Se0l9SA4o5cJ67nnSOrs5ftIKvdAiCSLgGEeJQ9J3VLK/B34ZUKncs+K/Bewa4n86YQocGPZ9qSe3YEGd9ojgLQ/C98jXPnX+dXxCS5fAfSQNMT16yjpYEntgL5mNg+4lLAO3TlLKTPbDHwm6TgXnUlY0y6JmX1AWN4ZaGZrgH8QQtYWHPdCwokJv0Oon+tczEJC1ML2HqlvREpfnxPC1Baee1onaRey50eWzUqx3fhJ2s/MFpvZ7whhlvtm1ozkQrzibsWYmUkaA9zot+J9SViH/BXhizmE8AeRAZea2UZgo6SDgEV+BbyFEKGvVGzxZJ8fSzoPeMwdXQPw/aJiU4BbJNUT5tBCQqS0rDY/kfSCwoNWnzGzS4ryl/lP/g9tW4zuB4AnJL1C+FPs7ZR2P/Cf6/WE9ePXXf6V3zZ3kzunDsCNhDtx7neZCP8fbCput4gJwDR3hmuAs8uUL7CYsEQDwWFfS3DgEO6umabwoItvCGvMjYmlpwIzCX9QL3Xds04aZwK3SbqaEL3udK+bNj9Is1kZngBmSBpNeMbjhZIGEmw4x/uI1JAYHTASiUTaGHGpJBKJRNoY0XFHIpFIGyM67kgkEmljRMcdiUQibYzouCORSKSNER13JBKJtDGi445EIpE2RnTckUgk0sb4H8S/ni8cW3P5AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import gaussian_kde\n",
    "density = gaussian_kde(adjusted.mean(axis=0))\n",
    "permutation_density = gaussian_kde(permuted.mean(axis=0))\n",
    "missing_z_density = gaussian_kde(unadjusted.mean(axis=0))\n",
    "xs = np.linspace(-1.,1.,1000)\n",
    "plt.plot(xs,density(xs), label='Fit model (not permuted)')\n",
    "plt.plot(xs, permutation_density(xs), label='Permuted model')\n",
    "plt.plot(xs, missing_z_density(xs), label='Model missing propensity score')\n",
    "plt.xlabel('Coefficient values for word counts')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "exog_sim_low = '../logdir/prop_sim/modesimple/beta01.0.beta110.0.exog0.0'\n",
    "exog_sim_high = '../logdir/prop_sim/modesimple/beta01.0.beta110.0.exog0.8'\n",
    "predictions_file = os.path.join(exog_sim_low, 'predict', 'test_results_keto.tsv')\n",
    "predict_df = pd.read_csv(predictions_file, delimiter='\\t')\n",
    "predict_df = predict_df.rename(columns={'index':'post_index'})\n",
    "term_counts,_,_ = load_term_counts()\n",
    "term_counts = term_counts.toarray()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model accuracy: 0.6096851788216013\n",
      "Mean and std. of the word coeffs: -0.016464023555452536 0.16569509106403862\n",
      "Model accuracy: 0.6108376991056782\n",
      "Mean and std. of the word coeffs: -0.015326543983318373 0.1627639900466207\n",
      "Model accuracy: 0.6095973119126371\n",
      "Mean and std. of the word coeffs: -0.012424147449095903 0.15932028828285638\n",
      "Model accuracy: 0.6087036598521325\n",
      "Mean and std. of the word coeffs: -0.01104095630703616 0.16641578254477737\n",
      "Model accuracy: 0.6058602869809292\n",
      "Mean and std. of the word coeffs: -0.012674520240197615 0.164004180674775\n",
      "Model accuracy: 0.6120621334601363\n",
      "Mean and std. of the word coeffs: -0.01667251411468293 0.16358118039087283\n",
      "Model accuracy: 0.6105448895209502\n",
      "Mean and std. of the word coeffs: -0.013781926949811967 0.15966880277178389\n",
      "Model accuracy: 0.6085096268124555\n",
      "Mean and std. of the word coeffs: -0.0162789379375924 0.17024585723303765\n",
      "Model accuracy: 0.6100081947712073\n",
      "Mean and std. of the word coeffs: -0.01488983582065513 0.16324582380766917\n",
      "Model accuracy: 0.6065880739509982\n",
      "Mean and std. of the word coeffs: -0.016613486826960795 0.1642806448206416\n"
     ]
    }
   ],
   "source": [
    "n_bootstraps = 10\n",
    "n_w = term_counts.shape[1]\n",
    "adjusted = np.zeros((n_bootstraps, n_w))\n",
    "permuted = np.zeros((n_bootstraps, n_w))\n",
    "unadjusted = np.zeros((n_bootstraps, n_w))\n",
    "for i in range(n_bootstraps):\n",
    "    sample = assign_split(predict_df,num_splits=2)\n",
    "    sample = sample[sample.split==0]\n",
    "    indices = sample.post_index.values\n",
    "    labels = sample.treatment.values\n",
    "    words = term_counts[indices, :]\n",
    "    propensity_score = logit(sample.treatment_probability.values)\n",
    "    all_features = np.column_stack((propensity_score, words))\n",
    "    unadjusted[i,:] = fit_treatment(words, labels, coeff_offset=0,verbose=True)\n",
    "    adjusted[i,:] = fit_treatment(all_features, labels)\n",
    "    np.random.shuffle(words)\n",
    "    permuted_features = np.column_stack((propensity_score, words))\n",
    "    permuted[i,:] = fit_treatment(permuted_features, labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEKCAYAAAACS67iAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xlc1VX++PHXh31HdlBUwAVQQEBAcDdzyyWtrCy3LE2blqmprGla5vdtZmqaabF1msmlMrO00rJVc9cUEVdAFEEFlB1kh3vv+f3xkZsom3C5l+U8e/C4du9neYP4vuee5X0UIQSSJElS52dm6gAkSZIkw5AJXZIkqYuQCV2SJKmLkAldkiSpi5AJXZIkqYuQCV2SJKmLkAldkiSpi5AJXZIkqYuQCV2SJKmLsDDmzdzd3YWfn58xbylJktTpJSQk5AshPJo7zqgJ3c/Pj0OHDhnzlpIkSZ2eoijnWnKc7HKRJEnqImRClyRJ6iJkQpckSeoijNqHLklNqa2tJTMzk6qqKlOHIkkmYWNjg6+vL5aWlq06XyZ0qcPIzMzE0dERPz8/FEUxdTiSZFRCCAoKCsjMzMTf379V15BdLlKHUVVVhZubm0zmUrekKApubm5t+oQqE7rUochkLnVnbf39l10uUvelrYXKIlDMwdYFzGT7Rurc5G+w1D3VVEBuMlzOgpLzkJ8KWg3m5uaEh4frvzIyMjh06BCPPvooADt27GDfvn3tHt6OHTuYNm1aq49JTEzkgQceaNW9i4uLee+991p1rjF88803JCUl3fB5Dg4OAOTl5TF58mRDh9UhyIQudT86LRSlg5k5eASDawBoqqD4HLa2thw5ckT/5efnR1RUFCtWrACMl9Db6u9//zuPPPJIq85tz4Su0WjafI3WJvQ6Hh4e+Pj4sHfv3jbH0tHIhC51P2WXQFsDPfqCpQ3YOIOTD1RfbvDwupZwRkYGH3zwAW+88Qbh4eHs3r273nEvvfQSCxYsYOLEifj5+fHVV1/x9NNPExoayuTJk6mtrQVg27ZtREREEBoayqJFi6iurgbgxx9/JCgoiJEjR/LVV1/pr1teXs6iRYuIjo4mIiKCTZs2NfntlZaWcuzYMYYMGaKPa9GiRYwdO5aAgAD9mxPA66+/TkhICCEhIbz55psAPPPMM6SlpREeHs5TTz1V79oZGRkEBQWxYMECwsLCuOOOO6ioqAAgISGBMWPGMHToUCZNmsTFixcBGDt2LH/+858ZM2YMb731FgsXLmTZsmWMGzeOgIAAdu7cyaJFiwgODmbhwoX6e9W1qAE2bNjAwoUL2bdvH5s3b+app54iPDyctLQ00tLSmDx5MkOHDmXUqFGkpKQAkJ6eTlxcHNHR0Tz//PP1vo+ZM2eydu3aJn+OnZHsQ5c6pL9+e5Kk7IYTbGsN6unEi1MDoTwfbFzA+veEgb0nlBdQWVlJeHg4AP7+/nz99df6Q/z8/Fi6dCkODg48+eSTDd4jLS2N7du3k5SURFxcHBs3buSf//wns2bNYsuWLUyePJmFCxeybds2Bg4cyPz583n//fdZunQpixcv5tdff6V///7cdddd+mv+7W9/46abbmLlypUUFxcTExPDzTff3Oj3eejQIUJCQuo9l5KSwvbt2yktLSUwMJBly5Zx7NgxVq1axYEDBxBCMGzYMMaMGcMrr7zCiRMnOHLkSIPXP3XqFB999BEjRoxg0aJFvPfeezz22GM88sgjbNq0CQ8PD9avX89zzz3HypUrAbXVv3PnTgAWLlxIUVERv/76K5s3b2b69Ons3buX//3vf0RHR3PkyBH938G1hg8fzowZM5g2bRp33HEHAOPHj+eDDz5gwIABHDhwgIceeohff/2Vxx57jGXLljF//nzefffdeteJioriL3/5S6M/w85KJnSpeykvAKEDR6/6zysKOHhia2PNkd92gY1Tqy4/ZcoULC0tCQ0NRavV6vtqQ0NDycjI4NSpU/j7+zNw4EAAFixYwLvvvsvYsWPx9/dnwIABAMydO5cPP/wQgJ9//pnNmzfzr3/9C1Cnd54/f77RGC5evIiHR/3CfFOnTsXa2hpra2s8PT3Jyclhz549zJo1C3t7ewBuu+02du/ezYwZM5r8Hnv37s2IESP0ca5YsYLJkydz4sQJJkyYAIBWq8XHx0d/ztVvUADTp09HURRCQ0Px8vIiNDQUgMGDB5ORkdFoQr9WWVkZ+/btY/bs2frn6j7x7N27l40bNwIwb948li9frj/G09OT7OzsFt2jM5EJXeqQXpw+2PAXFQLyUsDSHixtr3/dzlV9rChodUK3trYGwMzMDEtLS/00NDMzMzQaDUKIRs9tbMqaEIKNGzcSGBhY7/mcnJwGj7e1tb1uLnNdXADm5ubNxtKUa+NUFAUhBIMHD2b//v0NnlP3pnFtPGZmZvViq/s5XXufxuZm63Q6evTo0einicZ+plVVVdjaNvA70MnJPnSp+6itUAc/6xL3tRQztaVeVQK6hgfvHB0dKS0tbXUIQUFBZGRkcObMGQA++eQTxowZQ1BQEOnp6aSlpQGwbt06/TmTJk3i7bff1ifgxMTEJu8RHBysv35TRo8ezTfffENFRQXl5eV8/fXXjBo1qtnv8fz58/rEvW7dOkaOHElgYCB5eXn652trazl58mSzMTTFy8uL5ORkdDpdva6vq+NzcnLC39+fL7/8ElDf/I4ePQrAiBEj+PzzzwGu6y9PTU29rluqK5AJXeo+qooBBWx7NHGQAgg1qTdg+vTpfP311w0OiraEjY0Nq1atYvbs2YSGhmJmZsbSpUuxsbHhww8/ZOrUqYwcOZK+ffvqz3n++eepra0lLCyMkJCQ6wb4rhUUFERJSUmzbzyRkZEsXLiQmJgYhg0bxgMPPEBERARubm6MGDGCkJCQ6wZFQX3DWLNmDWFhYRQWFrJs2TKsrKzYsGEDy5cvZ8iQIYSHh7d5NtArr7zCtGnTuOmmm+p139x999289tprREREkJaWxtq1a/noo48YMmQIgwcP1g8av/XWW7z77rtER0dTUlL/73P79u1MnTq1TfF1REpzH7sURVkJTANyhRAh17z2JPAa4CGEyG/uZlFRUUJucCE1Jjk5meDg4Pa7QW4ymFmA+4DGjxECck6CpR24BbRfLO3sjTfewNHRsdVz0RuTkZHBtGnTOHHihEGva2yjR49m06ZNuLi4mDqU6zT070BRlAQhRFRz57akhb4auG4WvqIovYEJQOOjM5LUUWiq1e4WG+emj1MUsHWG6lJ1vnontWzZsnp909Lv8vLyeOKJJzpkMm+rZhO6EGIXUNjAS28ATwOtG1mRJGOq60JpLqEDWDsDOjWpd1I2NjbMmzfP4Nf18/Pr9K1zDw8PZs6caeow2kWr+tAVRZkBZAkhjho4HklqH9WlYG4NFi1otVo7qPVdGlloJEkd1Q1PW1QUxQ54DpjYwuOXAEsA+vTpc6O3k6S2EwJqypsZDL2KYgZW9lBd1r5xSZKBtaaF3g/wB44qipIB+AKHFUXxbuhgIcSHQogoIUTUtYsdJMkoaitBaMHKoflj61g7grYaNDXtF5ckGdgNt9CFEMcBz7r/v5LUo1oyy0WSTKLmSkv7RhM6QE0pWLgZPiZJagfNttAVRVkH7AcCFUXJVBTl/vYPS5IMqKYMzK3AwqrZQ+vK54ZERDP7weVUFOUaIcDfvfnmm/piVy3VklK7bZWRkdHsQpyWHCO1r5bMcpkjhPARQlgKIXyFEB9d87qfbJ1LHZYQal94C1vndeVzT5w4gZWNLR/8d6V6jRbQats+zbE1CV2S6siVolLXpqm+0n9u3/yx1xg1chRn0s+BtoZPP/2UmJgYwsPDefDBB/XJ28HBgRdeeIFhw4axf/9+/Pz8+POf/0xcXBxRUVEcPnyYSZMm0a9fPz744APg+hb1ww8/zOrVq1mxYgXZ2dmMGzeOcePGAWphrri4OCIjI5k9ezZlZWr3UWOldq+2evVqZs6cyfTp0/H39+edd97h9ddfJyIigtjYWAoL1dnIR44cITY2lrCwMGbNmkVRURGglsMdMmQIcXFx9aoVarVannrqKaKjowkLC+M///nPDf9spfYhi3NJHdMPz8Cl422/jq5WXVBkaQc+4TDllRadptFo+GHrDiYPDyP52GHWr1/P3r17sbS05KGHHmLt2rXMnz+f8vJyQkJC+H//7//pz+3duzf79+/n8ccfZ+HChezdu5eqqioGDx7M0qVLG73no48+yuuvv8727dtxd3cnPz+fl19+ma1bt2Jvb8+rr77K66+/ztNPP91oqd1rnThxgsTERKqqqujfvz+vvvoqiYmJPP7443z88cf88Y9/ZP78+bz99tuMGTOGF154gb/+9a+8+eab3Hffffrnry4B8NFHH+Hs7Ex8fDzV1dWMGDGCiRMnyv1gOwCZ0KWuTegARZ2K2AJX10MfNXIk999zGx+u/4GEhASio6P1x3h6qvMCzM3Nuf322+tdo678bGhoKGVlZTg6OuLo6IiNjQ3FxcUtDv23334jKSlJX6q2pqaGuLg4UlJSGi21e61x48bp7+/s7Mz06dP1sR07doySkhKKi4sZM2YMoJbznT179nXPz5s3jx9++AFQPzUcO3aMDRs2AFBSUsLp06f1JYEl05EJXeqYWtiSblZeqlpvy71lyaauD10v/wxCU82CBQv4xz/+cd3xNjY2mJub13uuudKwFhYW6HQ6/fONlYYVQjBhwoR6lRdB7SJpaWv42vtfHVtT28EJIZos5/v2228zadKkes9nZGS0KCap/cg+dKnrEkKdg25p1/prWNkzfngEGzZsIDdXnfFSWFjIuXPnWn3Jvn37kpSURHV1NSUlJWzbtk3/2tWlYWNjY9m7d6++FG5FRQWpqalNltq9Uc7Ozri4uOgrR9aV8+3RowfOzs7s2bMHqF9+dtKkSbz//vv6LfVSU1MpLy9vdQyS4cgWutR1aSoBXZsT+qCBAbz84nNMnDgRnU6HpaUl7777br0Stzeid+/e3HnnnYSFhTFgwAAiIiL0ry1ZsoQpU6bg4+PD9u3bWb16NXPmzNHvwvPyyy8zcOBAfaldd3d3Ro4c2ab6KmvWrGHp0qVUVFQQEBDAqlWrAFi1ahWLFi3Czs6uXmv8gQceICMjg8jISIQQeHh48M0337T6/pLhNFs+15Bk+VypKQYvn1ueDyUXwCNY3Qy6NXQadXDW0QccG1wMLUkG1d7lcyWpc6qtUItstaQgV2PMLMDCRq0FI0kdnEzoUtdVW6HuHdrW6XRW9mpCN+KnWUlqDZnQpa5JCKitalv/eR1LO3VxklYW6pI6NpnQpa5JUw2I1vedX83yyu7wtXJJvtSxyYQudU2aK3O7LQyV0BWZ0KUOTyZ0qWsyZEJXzNSWfk1l268lSe1IJnSpa6qtUkvmmpk3f+xVFEWptxenRqPBw8ODaXP/oLbQWzgw6ufnR35+00VIW3JMYw4dOsSjjz56w+fdcsstN1R+oDN74IEHSEpKAuDvf/+7iaMxDpnQpa5JU9Wq1rm9vT0nTpygslJtjf/yyy/06tVLfWPoQAOjUVFRrFix4obP+/777+nRo4Vb8RmAIUoKt9b//vc/Bg0aBJg+oTdVZsGQZEKXuh4hWp3QAaZMmcKWLVsAdVn9nDlz1PnsQGFOFjNnziQsLIzY2FiOHTsGQEFBARMnTiQiIoIHH3yQqxfsNVZ6tzEODg4sX76coUOHcvPNN3Pw4EHGjh1LQEAAmzdvBuqX4N25cyfh4eGEh4cTERFBaWkpFy9eZPTo0epmHSEh+qX9dZ8KMjIyCA4OZvHixQwePJiJEyfq38Ti4+MJCwsjLi6Op556qsFNK3bs2MHo0aOZNWsWgwYNYunSpfr6NNeWFN62bRsRERGEhoayaNEi/apXPz8/li9fTkxMDDExMfoSB3l5edx+++1ER0cTHR3N3r17AXjppZdYtGiR/mdR94ZWXl7O1KlTGTJkCCEhIaxfvx6AsWPHcujQIZ555hl90bV7772X559/nrfeekv/vTz33HPXvTk2ds34+HiGDx/OkCFDiImJobS0lKqqKu677z5CQ0OJiIhg+/btgFq+ePbs2UyfPp2JE9UtmF977TV92eEXX3yxyd+DVhFCGO1r6NChQpIak5SUpP/zKwdeEQt/WNi6r+/ni4Wb7hALt9xb7/lXDrzSbAz29vbi6NGj4vbbbxeVlZViyJAhYvv27WLq1FuEyEoUDy9ZKF566SUhhBDbtm0TQ4YMEUII8cgjj4i//vWvQgghvvvuOwGIvLw8kZSUJKZNmyZqamqEEEIsW7ZMrFmzRgghRN++fUVeXt51MQDi+++/F0IIMXPmTDFhwgRRU1Mjjhw5or+fGtNUIYQQ06ZNE3v27BFCCFFaWipqa2vFv/71L/Hyyy8LIYTQaDTi8uXL9e6Znp4uzM3NRWJiohBCiNmzZ4tPPvlECCHE4MGDxd69e4UQQixfvlwMHjz4uhi3b98urK2tRVpamtBoNOLmm28WX375pT7+9evXCyGEqKysFL6+vuLUqVNCCCHmzZsn3njjDX0sdTGuWbNG//3MmTNH7N69WwghxLlz50RQUJAQQogXX3xRxMXFiaqqKpGXlydcXV1FTU2N2LBhg3jggQf0sRUXFwshhBgzZoyIj4/X/73WSU9PFxEREUIIIbRarQgICBD5+fn1vr+GrlldXS38/f3FwYMHhRBClJSU6H/WCxcuFEIIkZycLHr37i0qKyvFqlWrRK9evURBQYEQQoiffvpJLF68WOh0OqHVasXUqVPFzp07r/vZXv3voA5wSLQgx8oWutT1iCuVDFtYMvdaYWFhZGRksG7dOm655ZYrzypgacuefQf0few33XQTBQUFlJSUsGvXLubOnQvA1KlTcXFxAWDbtm360rvh4eFs27aNs2fPNnl/KysrJk+eDKhlbseMGYOlpSWhoaENVjQcMWIETzzxBCtWrKC4uBgLCwuio6NZtWoVL730EsePH8fR0fG68/z9/fWlgocOHUpGRgbFxcWUlpYyfPhwAO65555G44yJiSEgIABzc3PmzJmjL+R1dUnhU6dO4e/vry+tu2DBAnbt2qW/xpw5c/SP+/fvB2Dr1q08/PDDhIeHM2PGDC5fvqwvWDZ16lSsra1xd3fH09OTnJwcQkND2bp1K8uXL2f37t04Ozs3+fP18/PDzc2NxMREfv75ZyIiInBzq79vbEPXPHXqFD4+Pvoyyk5OTlhYWLBnzx7970RQUBB9+/YlNTUVgAkTJuDq6gqoZYfr7hcZGUlKSgqnT59uMtYb1WxxLkVRVgLTgFwhRMiV514DpgM1QBpwnxCie4y0SEaxPGZ5608uvQSlF8E77IYHRevMmDGDJ598kh07dlBQUKA+aWmL0OmuGxitKzPbULlZIUSjpXcbY2lpqb9WS0rePvPMM0ydOpXvv/+e2NhYtm7dyujRo9m1axdbtmxh3rx5PPXUU8yfP7/eeVeX1jU3N6eysrJeV1Fzrv1+6/7/6pLCzV3v6mvU/Vmn07F//35sbW2vO/7amDUaDQMHDiQhIYHvv/+eZ599lokTJ/LCCy80ed8HHniA1atXc+nSJRYtWnTd6w1dc+bMmY3+HTfG3t6+3nHPPvssDz74YJOxtUVLmjCrgcnXPPcLECKECANSgWcNHJcktV4rZ7hcbdGiRbzwwguEhob+/qSlHaNjI1j76ceA2o/s7u6Ok5MTo0eP1peY/eGHH/TbuI0fP96gpXcbkpaWRmhoKMuXLycqKoqUlBTOnTuHp6cnixcv5v777+fw4cMtupaLiwuOjo789ttvAHz++eeNHnvw4EHS09PR6XSsX7+ekSNHXndMUFAQGRkZ+v7xuvK8der6ptevX09cXBwAEydO5J133tEfU68+fQOys7Oxs7Nj7ty5PPnkkw1+r5aWlvpyvwCzZs3ixx9/JD4+/rq67o1dMygoiOzsbOLj4wEoLS1Fo9HU+7tPTU3l/PnzBAYGXnfNSZMmsXLlSv02gllZWfrfC0NptoUuhNilKIrfNc/9fNX//gbcYdCoJKkt2jAgWsfX15fHHnus/pOWtrz0xIPct/yfhIWFYWdnx5o1awB48cUXmTNnDpGRkYwZM4Y+ffoAMGjQIF5++WWDld5tyJtvvsn27dsxNzdn0KBBTJkyhc8//5zXXnsNS0tLHBwc+Pjjj1t8vY8++ojFixdjb2/P2LFjG+3CiIuL45lnnuH48eP6AdJr2djYsGrVKmbPno1GoyE6OrreNnzV1dUMGzYMnU6nr+u+YsUK/vCHPxAWFqZPmHX7sTbk+PHjPPXUU5iZmWFpacn7779/3TFLliwhLCyMyMhI1q5di5WVFePGjaNHjx7XbVDS2DWtrKxYv349jzzyCJWVldja2rJ161Yeeughli5dSmhoKBYWFqxevbreJ4k6EydOJDk5Wf/G5eDgwKeffqrf/cogWtLRDvgBJxp57VtgbkuuIwdFpaY0NBh0w3Q6IbIShSjObPu1rqXVCJF1WIiSbMNfuwMpLS3V//kf//iHePTRR6875upB2dZqbFDYGLRarRgyZIhITU01yf2bYrJBUUVRngM0wNomjlmiKMohRVEO5eXlteV2ktQ8Q9ZwuZaZOZhbX9k4o+vasmVLvemOf/nLX0wdkkElJSXRv39/xo8fr9+Xtato0QYXV7pcvhNXBkWvPLcAWAqMF0K0qMiF3OBCaopBNrioLIaidHUPUSv75o+/UYXp6opRr8GGv7YkYYINLhRFmQwsB2a0NJlLUku0pIHRJEPWcGmIpa26WlRnnJV/UvfS1t//ZhO6oijrgP1AoKIomYqi3A+8AzgCvyiKckRRlMZHLCSphWxsbCgoKGjbL7Wm7TNcmqQvpVvVPteXui0hBAUFBdjYtL4x0pJZLnMaePqjVt9Rkhrh6+tLZmYmbRprKb2kLtMvTDZcYFfTaeFyLuTWgvX1i3UkqS1sbGzw9fVt9fnNJnRJMhZLS0v8/f1bfwGdFv42BoYtgeiXDRfY1YSA12ZA4C1w6zvNHy9JRiSX/ktdR1EGaKvBI6j97qEo4B0Kl4633z0kqZVkQpe6jtwr3SwebZwp0xyvEPVeWjkwKnUsMqFLXUdeivroMbB97+Mdqn4SKDBsYSVJaiuZ0KWuIy8FnHu3/2Cl95X6LpdOtO99JOkGyYQudR15KeBxfVEkg3MfqE6NvHSs/e8lSTdAJnSpa9BpIf90+w6I1jG3VO+TI1voUsciE7rUNRRlqIuKPNt5QLSO12DISTLOvSSphWRCl7oG/YCoEVroAJ6DoOwSVBQa536S1AIyoUtdg37KohH60EFN6AC5spUudRwyoUtdQ94pcPI13nJ8r7qE3k4lBiSpFWRCl7qGvGTwNFJ3C4CjD9g4Q85J491TkpohE7rU+RlzhksdRQHPwbKFLnUoMqFLnV/dDBdjJnRQZ9TkJqsFuySpA5AJXer88k6pj8aasljHaxBUl8DlLOPeV5IaIRO61PnlXen2cG/nGi7XqpvpIuejSx2ETOhS55ebos5wsXEy7n3rPhHIqYtSB9GSLehWKoqSqyjKiauec1UU5RdFUU5feXRp3zAlqQl5Kcad4VLH1gUce8qELnUYLWmhrwYmX/PcM8A2IcQAYNuV/5ck49NpIT/V+AOidbwGyYQudRjNJnQhxC7g2vXNtwJrrvx5DTDTwHFJUssUnzPNDJc6nsGQlyo3u5A6hNb2oXsJIS4CXHn0NFxIknQDco1cw+VanoPVzS4Kz5rm/pJ0lXYfFFUUZYmiKIcURTnUpt3cJakh+qJcRqrhci39wKhcMSqZXmsTeo6iKD4AVx5zGztQCPGhECJKCBHl4eHRyttJUiPyTDTDpY5HIChmcsWo1CG0NqFvBhZc+fMCYJNhwpGkG5SbbLrWOYClLbgGyJouUofQkmmL64D9QKCiKJmKotwPvAJMUBTlNDDhyv9LknHVzXAx9grRa3kOki10qUOwaO4AIcScRl4ab+BYJOnG6Ge4mLCFDmpCT/4WairAys60sUjdmlwpKnVedTVcPEzdQg8GBOSfMm0cUrcnE7rUeRl7l6LG6HcvSjFtHFK3JxO61HnlnQKnXqab4VLHNQDMreSKUcnkZEKXOq+8ZNMtKLqauQW4B8qBUcnkZEKXOiedTl1y3xESOqjFwWRCl0xMJnSpcyo+B5pK01RZbIhnMFzOhKoSU0cidWMyoUudU11ruG5A0tTq4siTM10k05EJXeqc6gYgTT3DpY7c7ELqAGRClzqnvBRw7gPWjqaOROXcByztZT+6ZFIyoUudU26y6Zf8X83MTP20IFvokgnJhC51PtraKzVcOsiAaB3PQXJxkWRSMqFLnU/hWdDWdJwB0TqewVCeC+X5po5E6qZkQpc6H/0Mlw7U5QJXDYzKfnTJNGRClzqf3GR1Uwn3gaaOpD59TReZ0CXTkAld6nxyk8DFX91coiNx9AYbZzkwKpmMTOhS59PRZrjUURS1lZ4nB0Yl05AJXepcaqvUQdGONiBaxzNYbaELYepIpG5IJnSpcyk4DULb8aYs1vEcpNZzKb1o6kikbqhNCV1RlMcVRTmpKMoJRVHWKYpiY6jAJKlBHa2Gy7Xqqj/KfnTJBFqd0BVF6QU8CkQJIUIAc+BuQwUmSQ3KTQIzC3DtZ+pIGqafuij70SXja2uXiwVgqyiKBWAHZLc9JElqwqUTaivYwsrUkTTM3h3sPeXURckkWp3QhRBZwL+A88BFoEQI8bOhApOkBuWcAO9QU0fRtLqBUUkysrZ0ubgAtwL+QE/AXlGUuQ0ct0RRlEOKohzKy8trfaSSVJanDjZ6hZg6kqbVTV3U6UwdidTNtKXL5WYgXQiRJ4SoBb4Chl97kBDiQyFElBAiysPDow23k7q9nOPqY4dvoQdBbYW6q5IkGVFbEvp5IFZRFDtFURRgPCA7DqX2c6mzJPS63YvkwKhkXG3pQz8AbAAOA8evXOtDA8UlSde7dAKcfMHO1dSRNE1OXZRMxKItJwshXgReNFAsktS0S8c7fuscwMYJnHvLmS6S0cmVolLnUFupbmrRGRJrQYUkAAAgAElEQVQ6qK10mdAlI5MJXeoccpPVJf/eHXyGSx3PYPUNSKsxdSRSNyITutQ55JxQHztLC91zkLqrUuFZU0cidSMyoUudw6XjYOUIPfxMHUnL6EsAyIFRyXhkQpc6h4tH1e4Ws07yK+sRCCiyH10yqk7yr0Pq1rQauHgMekaaOpKWs7QFV3/ZQpeMSiZ0qePLSwZNJfTqRAkd5O5FktHJhC51fFmH1ceeEaaN40Z5BkNBmrrLkiQZQZsWFkmSUWQfVjdfdg1ol8trtDoullRRUlmLuZlCzx62ONtatv3CnsHqVMuC051ndo7UqcmELnV8WYfV1rmiGOySJZW1fH04k59O5pBwrogabf3KiL162DI20IPbh/oS2celdTfxqJvpkiwTumQUMqFLHVttlTqwOPxRg1yupKKWd7af5pPfzlFVqyPI25H5cX0Z4OWAi50VtVrBhaIKDp8r4uvELNYeOM/Qvi48NzX4xhO7W391dyU500UyEpnQpY4t5wToNAbpP990JIsXN5+kpLKWWRG9WDTCn5Bezo0eX1atYWNCJu9uP8Nt7+3jzihfXpg+GAfrFv6zsbACtwEyoUtGIxO61LHVDYi2YYbL5apanv/mBJuOZDO0rwv/d2sIg3o6NXueg7UFC4b7ccdQX97+9Qwf7kojPqOId+6JYHDPxt8I6vEMhqyEVscuSTdCznKROrbsw+oenU69WnX6hcIKbn9vH98du8ifJgxk/ZLYFiXzq9lbW/DMlCA+WxxLZY2W2R/sZ8ep3Jad7DlI3eiiuqwV0UvSjZEJXerYshJaPSB65EIxs97bS87lKj65P4ZHxg/Awrz1v/KxAW5sengEfm723L/mEF8eutD8SZ5XaqPnnWr1fSWppWRClzquikK1YmHvmBs+NT6jkHv/+xt2VhZ8/YcRDO/nbpCQvJxs+GJpHMP7ufH0xmPNJ/W63YtyTxrk/pLUFJnQpY4rM1597D3shk47cLaABSsPqsn3wTj6eTgYNCwHawv+Oz+Kkf3deXrjMb5OzGz8YBd/sHJQSxdIUjuTCV3quC4cAMX8hgZEE88XsXBVPD7ONny+JBZvZ5t2Cc3G0pwP50UxzN+VJ788xq7UvIYPNDMD7zC4eKRd4pCkq7UpoSuK0kNRlA2KoqQoipKsKEqcoQKTJM4fAJ8wsLJv0eFpeWUsWh2Pu6MV65bE4unUPsm8jq2VOf+dH8UATwceWnuYlEuXGz6wZ7i6H6rc7EJqZ21tob8F/CiECAKGAHLCrWQY2lp1QLSF3S05l6uY/9FBzBSFTxYNw9OxfZN5HUcbS1bdF429tTn3rYon53IDdVt8wtXiYvlyYFRqX61O6IqiOAGjgY8AhBA1QohiQwUmdXOXjqtJsAUDohU1Gu5bFU9RRQ2r7ovGz71lLXpD8XG2ZeXCaEoqa3lo7WFqNPXLCNAzXH3Mlt0uUvtqSws9AMgDVimKkqgoyv8URbnuX5KiKEsURTmkKMqhvLxG+hkl6VoXDqiPvWObPEwIwVNfHiP50mXevSeSMN8eRgjueoN7OvPPO8JIOFfEy1uuqYHu1h8s7WU/utTu2pLQLYBI4H0hRARQDjxz7UFCiA+FEFFCiCgPD4823E7qVi4cACdfcG56QdE7v55hy/GLPDM5iHFBnkYKrmHTwnqyZHQAH+8/x4aEq2a+mJmrYwGyhS61s7Yk9EwgUwhxpSnFBtQEL0ltI4Q6INqn6f7zH09c4t+/pDIrohdLRrdPad0b9fSkQIb3c+PPXx/nRFbJ7y/4hKvdSHJgVGpHrU7oQohLwAVFUQKvPDUekPttSW1XlAGl2U12t5y6VMoTXxxhSO8e/OO2UBQDltZtCwtzM96eE4GrnRWPrEukrPpKAu9ZNzCaatoApS6trbNcHgHWKopyDAgH/t72kKRuL2OP+ug/qsGXy6o1LFubgL21BR/OG4qNpbkRg2uem4M1K+ZEcK6gnOe/OYEQQm2hg+xHl9pVmxK6EOLIlf7xMCHETCFEkaECk7qxjD1g5w4eQde9JITgz18dJyO/nBV3R+DVznPNWyvG35U/3jyQrxOz2Hg4C9wHqAOjddUjJakdyJWiUscihJrQ/UY2WJDrs4Pn2Xw0mycmDCSun5sJAmy5P4zrT2yAK89/c4Iz+ZXgOxQyD5o6LKkLkwld6liK0uFypprQr3Eiq4S/fpvE6IEePDS2vwmCuzHmZgpv3R2BrZU5D392GE3PaHXFqCylK7UTmdCljqWu/9yvfv/55apa/vDZYVztrHjjziGYmXWMQdDmeDnZ8O/ZQ0i5VMraiz3VTaPlhhdSO5EJXepY9P3ngfqnhBA8s/EYmUWVvHNPBG4O1iYM8MaNC/Jk8Sh//p10ZWONC7LbRWofMqFLHUcj/ecf7z/H98cv8fSkQKL8XE0YYOs9NSkIP9+enKE3VWf3mjocqYuSCV3qOIrS4XJWvemKRy8U8/KWJMYHebJ4VMdYPNQaVhbq/PTDIhDt+YPUauQCI8nwZEKXOo703epjX3VAtKRC7Tf3dLTh352o37wxfd3sGTB0PPainE++/dnU4UhdkEzoUsdxdjs4+oBHIEIIntxwlJzLVbxzTwQ97KxMHZ1BRIyYBEDaoV/YczrfxNFIXY1M6FLHoNPC2R3Q7yZQFD7ak84vSTk8OyWYiD4upo7OcFwD0Dn2ZKJtCo9/cYS80mpTRyR1ITKhSx3DxSNQWQT9buLw+SJe+SGFSYO9uG+En6kjMyxFwSxgLCMskrlcWc2fvjyKTidMHZXURciELnUMab8CUOIznEc+S8Snhw3/vGNIhym6ZVABY7CoKuL10ebsSs3jv7vPmjoiqYuQCV3qGNK2I3yG8KctWeSWVvHuPZE421qaOqr24T8GgFvsU5kS4s1rP50i8bwsgyS1nUzokulVl8KFAxy1imRrci7P3RJssp2HjMLJB9wHoqTv5JXbwvBysuGRdYmUVNaaOjKpk5MJXTK9jD2g0/CvM72YEuLNguF+po6o/fmPgXP7cLYSrJgTwaWSKh5ff0T2p0ttIhO6ZHJVKb9QiTUXnYbw6h1hXbPf/FoBY6G2Ai4cYGhfF56fNohfU3J5a9tpU0cmdWIyoUsmpdMJio//xAHdIN68dxhONl203/xaAWPA3ApSfwRgflxfbo/05a1tp9malGPi4KTOSiZ0yaTW/bwHb00mdoMmEOrrbOpwjMfaUa1Zk/oTAIqi8LdZIYT0cuLx9Uc4mydL7Eo3rs0JXVEUc0VREhVF+c4QAUndx74z+Zzc/TUA0TfdbuJoTGDgZCg4DQVpANhYmvPB3KFYWpix5JMEOUgq3TBDtNAfA5INcB2pG8kqruThdYlMsz2OzskX5apyud3GwMnq45VuFwBfFzvevSeScwXlPLQ2gRqNzkTBSZ1RmxK6oii+wFTgf4YJR+oOqmq1LP0kAUVTRSzHMQuc3OB2c12eS1/wHASnfqj3dFw/N/5xWxh7zxTwl2+Oq5tMS1ILtLWF/ibwNCCbEVKLCCF4/psTHM8q4b+jqzDTVMKASaYOy3QGToLz+6GisN7Tdwz15dHxA/jiUCbv7UgzUXBSZ9PqhK4oyjQgVwjR5H5aiqIsURTlkKIoh/Ly8lp7O6mLWHvgPF8mZPLoTf2JrD4IFrb16p93O4NuBZ0GUrZc99LjNw9gZnhPXvvpFBsSMk0QnNTZtKWFPgKYoShKBvA5cJOiKJ9ee5AQ4kMhRJQQIsrDw6MNt5M6u0MZhfz125OMDfTgsfED1L7jgDFgaWvq0EzHJxxc/ODk19e9pCgKr94RxqgB7jy94SjfH79o/PikTqXVCV0I8awQwlcI4QfcDfwqhJhrsMikLuV8QQVLPkmgVw9b3rorAvOCVCg+r3Y5dGeKAoNnqaWDr+l2AbC2MOc/84YS2ceFxz5PZHtKrvFjlDoNOQ9daneXq2q5f008Wp1g5cJonO0sf5/Z0Z37z+sMngVCC8nfNviynZUFK++LJtDbkaWfJrAvTW6MITXMIAldCLFDCDHNENeSuhaNVscjnyWSnl/O+/dGEuDhoL5w+mfwCgXnXqYNsCPwDgPXgAa7Xeo42Vjy8aJh9HG1475V8ew4JVvq0vUsTB2A1LW9vCWZnal5vHJbKMP7u6tPVhbB+d9g5OOmDQ7Q6DTkVOSQV5FHXmUe+ZX55FXkUVBVQGlNKWU1ZZRryimvKadCU4FWp0XU/ScE5oo5tpa22FrYYmdhh62lLS7WLrjbuuu/vO296evUFzcbt4br1NR1u+x5E8rywKHhsSZXeys+XxLL/JUHWfzxId6eE8HkEJ92/glJnYlM6FK7+WhPOqv3ZbB4lD93x/T5/YUz29QuBiP1n+uEjqzSLNJK0jh/+TznS8+TWZrJhdILZJdloxGaesebKWa42bjhZOWEvZU99hb2eNl5YWthi4WZBQpqUlYUBZ3QUVlbSaWmkgpNBSVVJWSUZJBXkUeNrqbede0t7enr1Bd/Z3+CXYMJdg0myC0IJysnCLkDdv8bTmyA2GWNfi9uDtZ8tjiW+1Yd5A+fJfLaHVpui/Q1/A9N6pRkQpfaxTeJWfzfd0lMCfHmmSnB9V889T3YuUOvoQa/b35lPqeLTnO66DRnis9wuug0aSVpVGoq9cc4WjrS26k3wW7BTPSbiK+DL552nnjYeeBu646LtQvmZuZtikMIQWltKfkV+WSXZ3Pu8jnOXT7H+cvnSchJYMvZ36cp+jr4EukVSZRvMFFHP8F32NImK04621ryyf3DeGDNIZ744igXS6p4aGy/7lGlUmqSYsxVaFFRUeLQoUNGu59kGjtT87h/dTxRfi6svi8GG8urkqOmGv7ZD0Jugxkr2nSf4qpiThac5GTBSU7kn+BkwUlyK37vW3a1cWVAjwEMcFG/+vXoR1/HvjhbO5s8+RVWFZJSkEJSYRIn80+SkJNAUbW6a5GntSvDeo1gtO9ohvcarrbgG1BVq2X5xmNsOpLN7KG+/G1WKFYWcp5DV6QoSoIQIqq542QLXTKoIxeKWfZpAgO8HPlwflT9ZA5wdifUlELw9Bu6bqWmkpP5Jzmef1yfwLPKsvSv+zn5EeUVxWC3wQx0HciAHgNws3UzxLfULlxtXBneazjDew0H1Bb92UuHObRuFod69WR31m6+Pfst5oo5EZ4RjPEdw7g+4+jr1Fd/DRtLc968K5y+bvas2HaazKJK3p8bSQ87K1N9W5KJyRa6ZDCnc0q58z/7cbCxYOOy4Xg62lx/0KaHIWkTPHUGLKwbvVZ+ZT6JuYkk5iZyJPcIyQXJ+r7unvY9Gew+mMFugwlxDyHYLbjRVmyn8/m9cOEg2j+e4HhRMrsyd7EzcyepRakABLsGM8lvEpP9J9PL4fcZQl8dzuSZjcfxcLTm/bmRXXsLv26opS10mdAlgziTW8bdH/6GosCXD8bh525//UFaDfx7IPS7CW7/vZ6bTuhIK07TJ+/E3EQyy9Sl7lZmVoS4hxDuGU6EZwSh7qEduuXdZsnfwfp74Z4v6g0aXyy7yNbzW/kx/UeO5R8DIMw9jCn+U7gl4BZcbVw5eqGYh9YeJq+0mhemD+LeYX1M3rUkGYZM6JLRpOeXc9d/9qMTgs+XxNLf07GRA3fDmmloZq/ilFcg8Zfiic+JJzEnkdLaUkDtiojwjCDCM4Jwz3AGuQ7C0ryb7GIEoKmB14PAbxTcuabBQzJLM/kp4yd+yviJ5MJkLMwsGNd7HDP7z2RQj2j+9MVxdqbmMTO8J/83MwTH7rILVBcmE7pkFOcKyrnrP79Ro9Xx+ZJYBnpdn8w1Og2nCk8Rv/Ml4vOOctjJlbLacgD8nf0Z6jVUTeAe4fR27C1blT88A4c+gieSwd69yUNPF53mmzPf8G3atxRVF+Fp58mMgBmU5keycsdlfJxtef3OIQwL6MKfaroBmdCldpeRX849//2Nylotny2OJdhH7cfWJ/ArLfDDOYcpq1W3VPPHmuiBtxLtHU2UdxTutk0nrG4pNwXeGwY3/xVG/rFFp9Rqa9mZuZOvz3zNnqw9CCEIdY3h7NlwLl3qy5JR/Xli4kCsLdo2HVMyDZnQpXaVcukyc/93EK1Ox8f3R2Npm8Nv2b9dn8Cd/Yn2iibawpmoH1/E/db/QNidJo6+E1h1C1zOgkcSwezGpiLmlOew8fRGNqRuIK8yDzvFk8JLUfhZjeWft8US0celnYKW2otM6FK7OXy+iIUf/4KlwxmGDcrnZFEChVVqpUA/Jz+ivaPVFrhXFB52V5axf/cEHPlMnd1i7WDC6DuJ4xtg4/0wdyP0v7lVl6jV1bLt/DY+T/mchJwEEJbUloQxuc8d/O2Wydhby1nLnYVM6JJBldeWE38pnq+Sf+XXc3tRrNQFPG42bsT2jCXOJ45Yn1i87L2uP1lbC/8OBP8xMHuVkSPvpDTV8Pog6BMLd69t8+VSi1L55ORnfJv2HVqqMa/uz6LQBTwcOwMzRS5G6uhkQpfaRKPTcLLgJPuz97M/ez/H8o6hERqEzgIb7QAWRkxkgv8oBroMbH4Q8/QvsPYOuHsdBN1inG+gK9j6EuxdAX88brCqlJdrLrPi4Fq+TF2HzrwIW7y4P2wB80Jux87SziD3kAxPJnTphl24fIF92fvYf3E/By8epLS2FAWFINcgbDRB7D3uRoRnBP+bH3tjqxG/ehBSf4AnTze5mEi6RlEGvBUOY56GcX826KUramp49qe1bM3+EjObC1ibOTAn6E7mDrqn4U9ZkknJhC41q6S6hAMXD7D/otoKr1tK72PvQ1zPOOJ84oj0jOb1H7P4PP4CM4b05LXZYTc2U6KmAv414Ertlrfb6Tvpwj67CzIPweMnwbKBlbdtdL6gnOVbNpNYshlLx5OYmZlxi/8U5g2axyC3QQa/n9Q6spaLdJ1abS1H8o7ou1FOFpxEILC3tCfaO5oFgxcQ5xNHX6e+KIpCaVUtD609zO7T+Tw8rj9PTBiImdkNzhFP2gQ1ZRB2V/t8U11d3MOwZhocXQdR9xn88n3c7Fk3fw47U2/mhe92comtfC9+4buz3xHlFcX8QfMZ7Tu6zdUnJeOQLfQuTAhBWnEa+y/uZ1/2PhJyEqjUVGKumBPqHqq2wnvGEeIegqVZ/dWEGfnlLPnkEGl55fx9Vgh3Rfdp5C7NWDkFynLgkQR1IwfpxggB/x0H1aXwh/gbnsJ4I2o0OlbuTWfF9mMIhwM4ex2gXJdPb8fe3Bt8LzP7z8TesoGSDlK7a/cuF0VRegMfA96ADvhQCPFWU+fIhN7+8ivz2Z+9n98u/sZv2b+RW6nORvFz8iPWJ5a4nnFEe0fjaNXI8nzU8rePfHYYMzOFd+ZEMnJAKxf/5J+Gd6Lg5pc6xO5EndaJjbBhkdEGlS+VVPGPH5LZdOQCrp6pePkeJLMyGQdLB24bcBtzgubg6yg31TAmYyR0H8BHCHFYURRHIAGYKYRIauwcmdANr1JTSUJOgtqNcnE/p4tOA9DDugfDfIYR56O2wns69Gz2WkII/rPrLP/8MYWBXo78d34UvV3bMPPhlxdg3zvqEnZHOdDWaloNrIgAJx9Y9JPRPukkni/i/75L4vD5Yvr5FtI34BCJ+bvQoeOm3jcxd9BcIj0jZakGIzD6oKiiKJuAd4QQvzR2jEzobafVaUkuTNa3whNzE6nV1WJpZkmkZ6Q6J7xnHMGuwTc0v7iiRsPTG47x3bGLTA3z4bU7wrCzasMQS20VvDEYeg+DOZ+1/jqSKv4j2PJEmxYatYYQgm+PXeTVH1LIKq5kdLAl/fod45fMbyipLiHYNZh5g+Yx2W9y9yqiZmRGTeiKovgBu4AQIcTla15bAiwB6NOnz9Bz5861+X7dTWZppn4mysFLBympLgEg0CVQPxslwisCWwvbVl3/TG4pf1ibSGpuKU9PCmLpmIC2t7oOfwKbH4b5myBgbNuuJalVGN8Zqm7dt/hXo49HVNVq+WhPOu9tP0ONVse9sT4E+J9i45l1nC05i7utO3cH3s3swNm42rgaNbbuwGgJXVEUB2An8DchxFdNHStb6C1TUl3CwUsH9a3wC6UXAPC089R3oQzzGWaQwlZfHc7kua9PYGdlzht3hTN6YMM7zt8QIeC9ODAzh6V75GCoodS9Sc75HAKnmCSE3NIq/v1TKl8kXKCHrSV/vHkA/r0z+ezUWvZm7cXSzJIJfSdwZ+CdsjvGgIyS0BVFsQS+A34SQrze3PEyoTesRlvD0byj+gR+suAkOqHDzsKOGO8YfTeKv5O/wf6BVNZoeXHzCb44lMkwf1dWzInAy8lA85zPbINPb4OZ70P4PYa5pqT2pb8bc+WNci9YmG6ruZPZJbz8XTL7zxbQz8OepycH0b9nOV+kfsG3ad9SWltKP+d+zA6czfR+07vOjlImYoxBUQVYAxQKIVpU41MmdFWtrpaT+Sc5eOkgBy8e5EjeEaq11fWmE8b6xBLqEXrddEJDSM0p5eHPDnM6t4yHx/XnsfEDsDA34HS4NTMg75S6ZN2ESadLSv0JPrsTJvwfjHjUpKEIIfglKYdXf0whLa+cyD49eGZKMCG+NvyU8RNfpn7J8fzj2JjbMNl/MncOvJMQ9xDZam8FYyT0kcBu4DjqtEWAPwshvm/snO6a0LU6LSmFKRy4dICDlw5yOOcwlZpKQO0Hj/aOZpjPMIZ6DW1yOmFb6XSCVfsyePXHFBytLQzXxXK1jD2weipM/BsMf9iw15ZUa++Ec3vh4UPqzBcT02h1bEjI5I2tqeRcrmZ8kCdPTw4i0NuRpIIkvkz9ki1nt1CpqaR/j/7c2u9WpgZM/b0Sp9QsufTfhHRCx+mi02oL/NJBEi4l6LdYC3AOIMY7hhifGKK8onCxMU5t6uziSp788ij70goYH+TJK7eH4eFo4LoqQqh1vIvS4dFEsGzdIK3UjII0dYyi/3i4+7MOM0ZRWaNl1b503t+RRlm1htsifHli4kB69bClrKaM79O/Z9OZTRzLP4aZYsbwnsO5td+tjO09FhsLw5c16EpkQjeiuh16EnISOJx7mIScBIqriwHo7dhbTeDeMUR7Rxu9VSKEYNORbJ7fdAKtTvDCtEHcFd1O27yd+gHW3Q1T/w3RDxj++tLv9r0DPz8HM96ByHmmjqae4ooa3tuRxup9GQDMi+3Lg2MC8HRUk3Z6STrfpn3Lt2e/5VL5JRwtHZngN4FJfScR7RPdLt2MnZ1M6O2oSlPF8fzjagLPOczRvKNUaCoA6OXQi6FeQxnmM4wY7xi87b1NFmdWcSXPf3OCX1NyiezTg9fvDMfPvZ2WbtdWwrvD1Fb5g7tl33l70+ng4xmQdRge+AW8Bps6outkFVfy5i+pbDyciaW5GXOvSew6oSP+Ujyb0zaz9dxWKjQVOFs7M77PeCb2nUiMT4xM7lfIhG5AJdUlHMk9QkKumsBPFpxEo9OgoNDfpT+RnpEM9RpKpGdkhyg9qtUJ1uzL4F8/n0II+NPEgSwc7mfYgc9r/foy7HoNFnwH/qPa7z7S7y5fhA/HqiWJF28H+465EXR6fjnv/HqGb45kYWGmcO+wviwdE4DnVbOqqrXV7M3ay8/nfmbHhR2U15bjZOXEuN7jGNN7DHE+cThYdd+drmRCbyWNTsOZ4jMcyzumfuUfI70kHQALMwtC3EKI9FIT+BCPIThbO5s44voOny/ir5tPcjSzhDEDPXh5Zkjblu+3xIV4WDkJQmfDbf9p33tJ9WUmwKop4BmsLuKy7WHqiBqVkV/OO9vP8HWimtjnxPTh/pH+1/1+Vmur2Ze1j5/P/czOCzsprS3FwsyCoZ5DGeU7itG+o/Fz8utWs2VkQm+hvIo8feI+lneMkwUn9TNQXG1cCXMPI9QjlAjPCELdQzvs4E12cSWv/pjCpiPZeDpa89zUYGYM6dn+v/RVJfDBKECoi4hsOtYbXLeQ+hN8fi94h8I9X4BDx549UpfYv0nMQicEt4T6sHhUAEN6X/9mVKur5WjuUXZl7WJ35m7OFJ8Bfh+bivaOJsY7psvPmJEJvQEFlQUkFyaTXJBMUkESJwtOcrH8IqC2voNdgwnzCNMncV8H3w7fCigqr+G/u8+ycm86QsCS0QEsHdPPOBsAa2th7WzI2A0Lv4c+w9r/nlLDTv0AX94H9u4wZ52a3Du4iyWVrN6bwWcHzlNarSHG35XFowK4KcgT80bq7meVZbE7czd7s/fWmz3m5+SnT/BDPIbgbe/d4f/t3ohundCFEORW5JJcqCbu5IJkkgqTyK3I1R/Tx7EPg9wGqQncI4wg1yCszTvP9mjFFTX8b3c6q/dlUF6jYVpYT5ZPDsTXxUj7Quq0sOlhOPpZh5xp0S1lJ8K6OVBRCOOehbhHwLzj72FTWlXL+vgLrNyTTnZJFT2dbbg7pg93RfducvWyVqclpSiF+Ivx6vqO3MOU15YD4G7rTqh7KGEeYYS6hxLiHtKpa7l3m4RepakirSSNM0VnOF10mtPFp0kpTKGwqhAAM8UMfyd/gt2CCXYNJtgtmCDXoHZdwNOezhdUsGZ/BuvjL1BWrWFqmA+P3jSAQG8jfj/aWvhmGRz/EsY+C2OfMd69paaV5apVGZO/Bc/BMP4FGDipw8xVb0qtVscvSTl8duA8e87kY26mMD7IkzkxfRg5wB3LZgb166YPH8s/xvG84xzLP8a5y2oxQAWFvk59GegykEDXQPXRJbDTtOS7XELX6DScLz2vJu7i0/rHC6UX0Al1oaqVmRX9evQj0DWQYNdgBrkNYqDLwE6/m7lWJ9h7Jp9PfjvH1uQczBWFW0J9+MO4/sZN5KC2/jY+AGnb1GQx6k/Gvb/UPCEgeTNsfQkKz0KvoTBsGQy6tdNMJ83IL2dd/Hm+PJRJYXkNbvZWTA3z4dbwnkT2cWlxEi6pLuF4/nGO5x0npTCF1KJUMssy9a87WjoywGUA/s7++H2gggIAABV7SURBVDv74+fkR1+nvvRy7NWhpkx2qYT+RsIbfJr0KTW6GkBtdfdx7MMAlwH079Ff/9jHsU+X2ftQCEHyxVK+Tsxk05FsckurcbW34p6YPsyN7Yu3swkGZ7MS4IuFUHYJbvkXDF1g/BikltPWQuIn6iKkwjSw94TQOyDkDugV2Sla7dUaLTtO5bH5aDZbk3Ko1ujwdbFl0mBvxgd5Eu3v2mzL/VplNWWcKT5DalEqpwpPkVqUyrnL5yiqLtIfY6FY4Ovoi5+zH35OfvRy6KX/8nHwaXWp6tbqUgn9x/QfOVlwUp+4A5wDOuxsk7ao0eg4mF7ItpQctiXncr6wAgszhXFBntwW0YtxQZ7YWJrgDau2Cna+CnvfAqeecOcatdUndQ46nfqJ6tAqOPMLaGvAxQ8GTlY3y+g7Aqw6/qfYsmoNP5+8xLdHs9l7poAarQ5HawtGB3owLtCT2ADXNo0hlVSXkF6SzrnL58i4nEFGSQYZlzM4f/m8vjFZx83GTZ/gezr0pKdDTzXZ2/vgbe9t8F6BLpXQu6rKGi2JF4qITy8iPqOQw+eLqKjRYmVhxoh+bowP9uKWUB9c7U34MTljD3z3OOSnQvhcmPQy2Bqn/ozUDiqLIPk7SPpG/bvVVIG5NfQdDgFj1M1IvMPUEr0dWHm1hr1n8vk1JZdtKbnklVYD0KuHLcP8XRkW4EpkHxcCPBwanTHTUjqhI78yn+yybDLLMskuyyarLIussiyyy7K5WHYRjdDUO8fRyhFve2+87bzVR3tvpvhNobdT71bFIBN6B6LVCbKLK0nLKyP5YilJFy+TfPEyZ/PK0An1k2+QtxMxfi6MGuDB8P5ubdv+zRAKz8LPz0PKd/D/2zv34DiKO49/fnp4Zb0syTJ+CGRbPpNgUzkTIIcJBhsIcJcAJnEK3yUEktRRCUeo45VAkaQIVRfC5Y+QuyRnDMUjCQQCF4MhJBwvYxLAPP3mZcs2lm1ZtixZ79Xrd3/8eq3xsruSVqvVatOfqqmZ7enu+U3P7Hd6erp/PakaLrzTnEF5soeeTtj1ivmvr30RGtx0wAVlNtq3ZjHMXgyT52R080x/v/JefSuv72hk3Y5DvL7jEI3tVqOemJ/LvBmlnDijlPlVkzh+agk1U4ooLUhd+3hffx8HOg9Q11pHfUc99e227G/ff+R3c7iZu8+7m9Omn5bUMbygpwlV5XBnDw2tYQ60hmlo7aKhJUxdUye7DnXwUWM7dU2d9PYPlHNV2UROmF7KvOklLKgu4+SZFUyamCEfYNoarGnl9ZWQkw+LroWFV3vPiX8LtO6HHWuhdo0tLe7jYemxVnuf7WrwGT7ht6qy/UAbG3YfZvPew2zec5gte1vo6O47EmdKSYiayiJqphRTU1lEVflEZpRNpKpsIpXFE1Le86Wzt5M8yUt63lUv6IPQ16+0d/fSHralLdzn1gnCugfC2sO9tHT2cLCtm+6+/o/lX1qQx8zJRVRXFFI9uZCZFYXMqizihGmlTCrMEPEO0rwbXv0lvHWftbH+/b/A2d/PCH/bnjFA1d7SIuK+82VrrkGg+jTrMXPChTDp2DE2dGj09ys7GtvZ1tBG7YF2ag+0sf1AG7UH22nu6Dkq7oS8HKqcuFeVTWTqpAKmloaYWlLAMaUhppYWMLlowuj6RooiKwVdVenq6ac13ENrVy9tXSa2rV29tHb10BYeCGtx67auHtrDfbSFe+kICHJnT9/gBwRyBIpCeRSH8ihyS3Eol8IJeZQU5DGlJMQxJQVuHWKKW1L5SjdqdHfAtuesJ8S250By4FPLYdF19prt8UTo74P6jeZmYOtqaNhi4ceeCidcZAJfPnNsbUySwx097GnutKWpgz3Nnext7qKuuZO9zZ0cbAsTLZM5ApXFJu5TS0McU1rA1BIn/KWmBxHhzxlhGz5kmaDf8ef3eGjdR7SFe+nrH9zegvwcikP5lBREhDiXoglHC3IskY7EKQ4NrAvyc8bFwIMh0d0O+zZY98Odf7GaV28XFE+zkZ6f/hqUVY+1lZ7xwMFt8O4TsPUJu6fAej7NWwrzl2bVfdTb18/Btm72t3TZ0hqmocWaVve3drG/xX5H2u2D5OUIlcVWyfvhhfM4dVZFUjaka5LoC4CfA7nAPar6k0TxkxX0J9bv4e1dTRQX5FFSkE9xKO+IWEf/Li7IG3a/1Kwg3Aat+2xp2Qete906KswNwqJ8to0gPP4CmLVoXAwR92Qoh3aYsG9ZBfvWW1jVyTD/Equ5Z5G4J6K7t58DbWH2t3SZ4LeGqT/cdeT72o3nf4ITq5JzXpeOOUVzgQ+AzwF1wBvAP6vq1nhpMqkNfdzQ1wvtDQOC3FoPLXudSLvfrfsg3PLxtKFJUDLN2sFLZtgr8YyTYManM94jn2eccmiHdYncsipQcz8F5iyB6oXWRFNQOng+/f3Qth8O77alObAOt0BvGPp7IL8IQiXWlbas2vrXl8+0CktpFeRkR+UuHYK+ELhVVc93v28GUNXb46VJWtD7euyDTEcjtB+0dUejDUPvPGTrjsaB7b4ewJ1XTq65dC0os/XEciicDIUVtp5YcfTvgkmj3we3r9duysg5Rew/ItiBGnZ7w0CtOkJOHpRMd8s0G+xTMt2tp5l4l0yD0N/uhACeDOBQLWx53Lq+7l0P6r5blUyHijlWqZhQbBN09HRak2BHoxPxPSbYQQrKoOw4W+eFrBdWTzt0tdh/qKXu6P9KbsjEvaLGBL6iBipmQ9EU8xs/sRxCpRndJTPCUAV9JO/ZVcDuwO86YHT8pz59A7x1f+x9E0qgsNwJc4VduMgoUsHEvasFupptct3OJhP+vo+3dx1JFLlZciPLBBPRYaP2Mak3DH1h6O22dbRAB5lYboJcOt2mFYtsR0S6dAYUVmZNzcOTxVTU2Af2RddBuBXq3rDvN4215oqgfrOF93VDfqGNVp1YYbX4eUutB01ZNUw6zoQ8NIjfor4eOFwHTTttovJDO+yh0rTTumP2dHw8jeSY8OdOsGbH3An2n5fcgNAHBP8o8R8kPPpB8YU7YebCxOcwQkYi6LEeax+r7ovIlcCVANXVSbalzf8iTD3RatBFla5G7WrXyTgbUoXutkANvxE6mgZq+b1dVovu67al323HPOVByMkzG3NDgXXIaiZH3gwq7FxKpvn+3p7sJFQCc862ZbTIzbcaeMVsYMnR+1RtjEXTTug4CJ3NrnLXZJWs6P97f+9AuoFMjs4vYXiMlo80uFcYiaDXAcFxrMcCe6MjqepKYCVYk0tSR6o5y5ZUIWI3WKhk3Ha18ng8w0DEBkRl+KCokTKS9/Y3gLkiMltEJgDLgdWpMcvj8Xg8wyXpGrqq9orI1cAzWLfFe1V1S8os83g8Hs+wGFHnY1V9Gng6RbZ4PB6PZwT4rhIej8eTJXhB93g8nizBC7rH4/FkCV7QPR6PJ0vwgu7xeDxZQlrd54rIAWBXkskrgYMpNCdVeLuGh7dreHi7hkem2gUjs22mqg7qUS+tgj4SROTNoTinSTferuHh7Roe3q7hkal2QXps800uHo/HkyV4Qfd4PJ4sYTwJ+sqxNiAO3q7h4e0aHt6u4ZGpdkEabBs3begej8fjScx4qqF7PB6PJwEZJegiUiEiz4rIh25dHiPOAhF5VUS2iMhGEbk0sG+2iKxz6R9xbn3TYpeL92cRaRaRp6LC7xeRHSKy3i0LMsSusS6vy12cD0Xk8kD4GhF5P1Bex4zQngtcfttE5KYY+0Pu/Le58pgV2HezC39fRM4fiR2psktEZolIZ6B8VqTZrjNF5G0R6RWRZVH7Yl7TDLCrL1BeKXXzPQS7rhORrU6vnheRmYF9qS0vVc2YBfhP4Ca3fRNwR4w4xwNz3fYMYB9Q5n7/HljutlcA306XXW7fOcCFwFNR4fcDy8aivAaxa8zKC6gAat263G2Xu31rgFNSZEsusB2oASYAG4B5UXGuAla47eXAI257nosfAma7fHIzwK5ZwOZU30/DsGsW8Cng18H7OtE1HUu73L62MSyvJUCh2/524DqmvLwyqoYOXAw84LYfAJZGR1DVD1T1Q7e9F2gApoiIAGcDjyVKP1p2OXueB1pTdMyhkLRdGVBe5wPPquohVW0CngUuSNHxg3wG2KaqtaraDTzs7Itn72PAOa58LgYeVtWwqu4Atrn8xtqu0WRQu1R1p6puBKInxx3NazoSu0aTodj1oqpGJjR9DZvdDUahvDJN0Keq6j4At074qi0in8GeituByUCzqrrJAKnDJrJOu11x+A/3yvUzEQllgF1jXV6xJhkPHv8+93r8gxGK2GDHOSqOK4/DWPkMJe1Y2AUwW0TeEZGXRGRRimwaql2jkXa08y4QkTdF5DURSVXFJRm7vgn8Kcm0gzKiCS6SQUSeA6bF2HXLMPOZDvwGuFxV++P86YfchSdVdsXhZqAee/isBL4H3DbGdo11eSU6/ldUdY+IlAD/C1yGvUYnw1DOM16cEZXRIIzErn1Atao2isjJwOMiMl9VW9Jk12ikHe28q1V1r4jUAC+IyCZV3Z5Ou0Tkq8ApQGSC5JSXV9oFXVXPjbdPRPaLyHRV3ecEuyFOvFLgj8D3VfU1F3wQKBORPFebiTlp9WjalSDvfW4zLCL3ATdkgF1jXV51wOLA72OxtnNUdY9bt4rIQ9hrbbKCPpTJzCNx6kQkD5gEHBpi2mRJ2i61BtgwgKq+JSLbsW9Lb6bJrkRpF0elXZMCmyJ5J30tXPMsqlorImuAk7A3+7TYJSLnYpWds1Q1HEi7OCrtmpEYk2lNLquByJfey4EnoiOI9cRYBfxaVR+NhLub/EVgWaL0o2VXIpyoRdqtlwKbx9quDCivZ4DzRKRcrBfMecAzIpInIpUAIpIPfIGRlddQJjMP2rsMeMGVz2pguettMhuYC7w+AltSYpeITBGRXABX45yLfVBLl13xiHlNx9ouZ0/IbVcCnwW2pssuETkJuAu4SFWDlZvUl9dofPlNdsHaB58HPnTrChd+CnCP2/4q0AOsDywL3L4a7A+3DXgUCKXLLvf7ZeAA0Ik9fc934S8AmzBh+i1QnCF2jXV5fcMdexvwdRdWBLwFbAS2AD9nhD1LgH8CPsBqZLe4sNuwPxhAgTv/ba48agJpb3Hp3gf+McX3e1J2AV9yZbMBeBu4MM12neruo3agEdiS6JqOtV3A6e7/t8Gtv5lmu54D9jOgV6tHq7z8SFGPx+PJEjKtycXj8Xg8SeIF3ePxeLIEL+gej8eTJXhB93g8nizBC7rH4/FkCV7QxykiMk1EHhaR7c6T29MicnySeV0jIu+KyIOuz/Vzbtj9pSJyj4jMS5D2olge5oZ43DIRuSqZtHHya0tVXnHy/51z33DtaB5nuIh581w2eMyUH3exiJye7uN64pP2kaKekeMGKK0CHlDV5S5sATAV6w87XK7C+ljvEJHTgHxVjbj4fSRRQlVdzdAHnkRT5o79qyTTpw0RmQacrqozB408kCYyCjeVdqQ8zxGwGGgDXhljOzwRUtnB3i/pWTAviWvj7BPgp9ggpk3ApYF9N2Ij2zYCP3JhK4BuF/d72ACHw9gAiDkEXNlinuDexgZoPO/CrgB+4banYL5X3nDLZ134rcC9Lq9a4BoX/jA22Gk98NOo87gDuCrw+1bgeqAYG6z0trP54kCcNrdeTMBVMPAL4Aq3fTLwEjZ46Rlgugu/Bhs9uBHzsBhdrhsDti4CFmCe8zZiD9eg698fu2NcH5XHJuwhJtjAl6+58N8A52IDie5z8d4BlgTK+FHgSWyQmrhz2oq5wHiaGO6Zgb/DBrVEBiDNIc79MUiZ7QR+FCjzT2KuauuBPYEy+bLLdwNx7k+/jLI2jLUBfkniopn4/CzOvi9hbjhzsRr7R8B0bFjxSveHzgGeAs50aXYClW47+o+9BhvhOQXzDDfbhUdGf17BgKA/BJzhtquBd932rVgtLgRUOjHLJ4Ffb8zXxkuB31tdnnlAqQurxB5AkQFyCQXdHfMVYIoLvxS4123vxY2UxfnXj7LnKFsxIT/Lbd8G3Bkor1/FOacVwOeBE7EH3t0u/EPsQXU9cJ8L+6S7dgXO9rpAmX8xcI1nAM3EFvR1wCVuuwAoJP79EbPMAvfHd9z2VQyM2r4VuCGQZhNQFa8M/TL6i29yyT7OAH6nqn3AfhF5CRsSfSYm6u+4eMWYD5C1Q8z3NKzWtQNAVQ/FiHMuMC/g+LLUeUwE+KOaU6KwiDRgYhIXVX1HRI4RkRnYw6RJVT9y/l1+LCJnYn6vq1xe9UM4h09gYvqsszEX81wIJtAPisjjwOOJMhGRSZhgveSCHsBq0BHiNVO9jF2HXcD/AFeKSBXmcKtNRM4A/tud/3sisgtzugXOb7bbPpOBa7xXRF6IYWMJJq6rXH5dLjze/TGYp8Y/uPVb2AMlFn8F7heR3wfie9KIF/TxyRYGnGpFE893uAC3q+pdSR5TGNy1Zw6wUFU7j0po4hkOBPUxtHvvMew8p2HNMwBfwQT+ZFXtEZGdWO0zSC9Hf/CP7BfMv8fCGMf6PCaUFwE/cO5ok22rbo8Tvhb4N+xN4xbgEuz8Xg7YN9Q8B7sWie6DWMQrswiR6xf32qnqt0TkH7CyXC8iC1S1cRA7PSnE93IZn7wAhETkXyMBInKqiJyFicalIpIrIlMwkXoday/+hogUu/hVMry5Ol8FznJeBxGRihhx/g+4OmDTYHOntgIlCfY/jHmvW8bAzEqTgAYn5kuAWB8pd2FvCiFXmz7Hhb+PzW610NmXLyLzRSQHOE5VXwS+i7VzF8czSlUPA00yMLHEZVibeUJUdTfWTDRXVWuBv2CulCOCvhZ7YOF6LFU7m6NZi3mBzHWePJfEOFYL5nZ3qcsvJCKFxL8/4pVZIo66fiIyR1XXqeoPMffMx8VN6RkVfA19HKKqKiKXAHe6LoNdWDvnv2N/2IXYhykFvquq9UC9iJwAvOpqzG2Y58oh+VBX1QMiciXwByeADcDnoqJdA/xSRDZi99Za4FsJ8mwUkb+KyGbgT6p6Y9T+La7pYI8O+JR/EHhSRN7EPsa9FyPf3e61fyPWPv2OC+923fv+y4lWHnAn1jPoty5MsO8TzYMUyeXACieStcDXB4kfYR3W1AMm5Ldjwg7W22eFiGzCasxXqGpYPj53yyrsw/gmZ3u8h8llwF0ichvmofTLLm2s+4NYZTYITwKPicjFwHeAa0VkLlaGz7tjeNKI97bo8Xg8WYJvcvF4PJ4swQu6x+PxZAle0D0ejydL8ILu8Xg8WYIXdI/H48kSvKB7PB5PluAF3ePxeLIEL+gej8eTJfw/uWWmGqK0E+EAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import gaussian_kde\n",
    "density = gaussian_kde(adjusted.mean(axis=0))\n",
    "permutation_density = gaussian_kde(permuted.mean(axis=0))\n",
    "missing_z_density = gaussian_kde(unadjusted.mean(axis=0))\n",
    "xs = np.linspace(-0.2,0.2,1000)\n",
    "plt.plot(xs,density(xs), label='Fit model (not permuted)')\n",
    "plt.plot(xs, permutation_density(xs), label='Permuted model')\n",
    "plt.plot(xs, missing_z_density(xs), label='Model missing propensity score')\n",
    "plt.xlabel('Coefficient values for word counts')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}