gravity_plot.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from scipy.optimize import curve_fit\n",
    "from itertools import product\n",
    "from tqdm import tqdm\n",
    "import json\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['font.size'] = 18"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# reading in byte-saved variables created by preferences.py that are the same for each run\n",
    "with open('adjacency.pickle', 'rb') as f:\n",
    "    # A is the adjacency matrix\n",
    "    # indexmap maps user_ids to the row indices of A\n",
    "    # indexmap_back maps row indices of A to user_ids\n",
    "    # seed is a boolean vector of the initial infected nodes\n",
    "    A, indexmap, indexmap_back, seed = pickle.load(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>id</th>\n",
       "      <th>comm_mem</th>\n",
       "      <th>cityid</th>\n",
       "      <th>month</th>\n",
       "      <th>node_id</th>\n",
       "      <th>ABM1</th>\n",
       "      <th>ABM2</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>new_id</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>24550</th>\n",
       "      <td>72367576</td>\n",
       "      <td>55</td>\n",
       "      <td>78567</td>\n",
       "      <td>78.0</td>\n",
       "      <td>72367576</td>\n",
       "      <td>96.0</td>\n",
       "      <td>55.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>68</th>\n",
       "      <td>64248926</td>\n",
       "      <td>55</td>\n",
       "      <td>379</td>\n",
       "      <td>120.0</td>\n",
       "      <td>64248926</td>\n",
       "      <td>77.0</td>\n",
       "      <td>81.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>23282</th>\n",
       "      <td>70127822</td>\n",
       "      <td>218</td>\n",
       "      <td>46758</td>\n",
       "      <td>94.0</td>\n",
       "      <td>70127822</td>\n",
       "      <td>69.0</td>\n",
       "      <td>65.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>39036</th>\n",
       "      <td>71295318</td>\n",
       "      <td>167</td>\n",
       "      <td>393</td>\n",
       "      <td>84.0</td>\n",
       "      <td>71295318</td>\n",
       "      <td>53.0</td>\n",
       "      <td>71.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22556</th>\n",
       "      <td>71889448</td>\n",
       "      <td>55</td>\n",
       "      <td>13943</td>\n",
       "      <td>80.0</td>\n",
       "      <td>71889448</td>\n",
       "      <td>49.0</td>\n",
       "      <td>72.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              id  comm_mem  cityid  month   node_id  ABM1  ABM2\n",
       "new_id                                                         \n",
       "24550   72367576        55   78567   78.0  72367576  96.0  55.0\n",
       "68      64248926        55     379  120.0  64248926  77.0  81.0\n",
       "23282   70127822       218   46758   94.0  70127822  69.0  65.0\n",
       "39036   71295318       167     393   84.0  71295318  53.0  71.0\n",
       "22556   71889448        55   13943   80.0  71889448  49.0  72.0"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# adoption times of nodes, index is the adjacency matrix index\n",
    "original_timeline = pd.read_csv('./vertices_fs.csv')\n",
    "abm_timelines = pd.read_csv('adoption_times.csv')\n",
    "timelines = pd.merge(original_timeline,abm_timelines,left_on=\"id\",right_on=\"node_id\")\n",
    "timelines[\"new_id\"] = timelines[\"id\"].map(indexmap)\n",
    "timelines.set_index(\"new_id\",inplace=True)\n",
    "timelines.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculating the adoption category time bins\n",
    "percentile_dict = {}\n",
    "for col in [\"month\",\"ABM1\",\"ABM2\"]:\n",
    "    percentile_dict[col] = np.percentile(timelines[col].dropna(),[0,2.5,16,50,84,100])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# getting the adjacency matrix ids of innovators and early adopters\n",
    "categories = {}\n",
    "category_names = {\n",
    "    0: \"innovators\",\n",
    "    1: \"early_adopters\",\n",
    "    2: \"early_majority\",\n",
    "    3: \"late_majority\",\n",
    "    4: \"laggards\"\n",
    "}\n",
    "\n",
    "\n",
    "for i in range(5):\n",
    "    category = category_names[i]\n",
    "    for col in [\"month\",\"ABM1\",\"ABM2\"]:\n",
    "        if category in categories:\n",
    "            categories[category][col] = np.array(list(map(lambda x: indexmap[x],list(timelines[\"id\"][(timelines[col]>=percentile_dict[col][i])&(timelines[col]<percentile_dict[col][i+1])]))))\n",
    "        else:\n",
    "            categories[category] = {}\n",
    "            categories[category][col] = np.array(list(map(lambda x: indexmap[x],list(timelines[\"id\"][(timelines[col]>=percentile_dict[col][i])&(timelines[col]<percentile_dict[col][i+1])]))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "innovators month 6428\n",
      "innovators ABM1 6776\n",
      "innovators ABM2 6555\n",
      "\n",
      "early_adopters month 32821\n",
      "early_adopters ABM1 36279\n",
      "early_adopters ABM2 34595\n",
      "\n",
      "early_majority month 94511\n",
      "early_majority ABM1 91575\n",
      "early_majority ABM2 90289\n",
      "\n",
      "late_majority month 93013\n",
      "late_majority ABM1 91004\n",
      "late_majority ABM2 92999\n",
      "\n",
      "laggards month 45101\n",
      "laggards ABM1 46079\n",
      "laggards ABM2 47010\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# checking how many users have been sorted into the different categories\n",
    "for k in categories:\n",
    "    for k2 in categories[k]:\n",
    "        print(k,k2,len(categories[k][k2]))\n",
    "    print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# city-city distance matrix\n",
    "citydist = pd.read_csv(\"city_dist.csv\",sep=\";\")\n",
    "citydist[\"tup\"] = citydist.apply(lambda row: (int(row[\"cityid1\"]),int(row[\"cityid2\"])),axis=1)\n",
    "cd = citydist.fillna(0).set_index(\"tup\")[\"distkm\"].to_dict()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# creating dicts for faster access by user_id\n",
    "cities = timelines[\"cityid\"].to_dict()\n",
    "\n",
    "# creating set from edgelist for fast containment checks\n",
    "elist_set = set(list(zip(*A.nonzero())))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# dumping output\n",
    "output = \"gravity_res_big.json\"\n",
    "# clear file (it should be more sophisticated!)\n",
    "open(output,\"w\").close()\n",
    "\n",
    "for col in [\"month\",\"ABM1\",\"ABM2\"]:\n",
    "    for cat in [category_names[i] for i in range(5)]:\n",
    "        print(col,cat)\n",
    "        gravity_dict={}\n",
    "        def attr(e,b=10):\n",
    "            global gravity_dict\n",
    "            if e[0]<e[1]:\n",
    "                edge=int(e in elist_set)\n",
    "                d=cd[(cities[e[0]],cities[e[1]])]\n",
    "                d_trunc = int(d/b)*b\n",
    "                if d_trunc in gravity_dict:\n",
    "                    gravity_dict[d_trunc][edge]+=1\n",
    "                else:\n",
    "                    gravity_dict[d_trunc] = {edge: 1, int(not edge):0}\n",
    "        total_pairs = len(categories[cat][col])**2\n",
    "        for i,e in enumerate(product(categories[cat][col],categories[cat][col])):\n",
    "            attr(e)\n",
    "            if i%1000000==0:\n",
    "                print(str(round(100*i/total_pairs,2))+'%',end=\", \")\n",
    "        \n",
    "        # writing calculated output\n",
    "        print(\"Dumping gravity data row.\")\n",
    "        f = open(output,\"a\")\n",
    "        f.write(json.dumps({\"data_source\":col,\"category\":cat,\"gravity_dict\":gravity_dict}))\n",
    "        f.write(\"\\n\")\n",
    "        f.close()\n",
    "        print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "# exporting gravity plot data for Balazs\n",
    "d = []\n",
    "for col in pd.DataFrame(gravity_res):\n",
    "    d.append(pd.DataFrame(gravity_res)[col].map(lambda d: d[1]/(d[1]+d[0])))\n",
    "d = pd.DataFrame(d).T\n",
    "d.index = pd.DataFrame(gravity_res).index\n",
    "d.reset_index().rename({\"index\":\"d\",\"month\":\"data\"},axis=1).to_csv(\"gravity_plot.csv\",index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/ubuntu/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:10: RuntimeWarning: divide by zero encountered in log\n",
      "  # Remove the CWD from sys.path while we load stuff.\n",
      "/home/ubuntu/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:16: RuntimeWarning: divide by zero encountered in log\n",
      "  app.launch_new_instance()\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAokAAAG/CAYAAAA0HE0MAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3XlcVOX+wPHPA4IsAmoIbgmYWqJeyVzKMs3Ubm6llqZQ5lpdyyzLFu1q3tQyzSWtLipu4NL9aZlXr1o3qNTKJUgEb5oh7ltuKILCPL8/RkaHGWCAGWaA7/v1mpfynOec8z0a9PV7nkVprRFCCCGEEOJWbs4OQAghhBBCuB5JEoUQQgghhAVJEoUQQgghhAVJEoUQQgghhAVJEoUQQgghhAVJEoUQQgghhAVJEoUQQgghhAVJEoUQQgghhAVJEoUQQgghhIUqzg6gIggMDNShoaEOv8+VK1fw9fV1+H2EEKKk5OeUEK5v9+7dZ7XWtYrqJ0miHYSGhrJr1y6H3ychIYFOnTo5/D5CCFFS8nNKCNenlEq3pZ+8bhZCCCGEEBYkSRRCCCGEEBYkSRRCCCGEEBYkSRRCCCGEEBYkSRRCCCGEEBYkSRRCCCGEEBYkSRRCCCGEEBZknUQhhBB2l52dzblz58jIyCA3N9fZ4QhR4Xl6ehIYGEhAQIDdrilJ4i2UUp8AvYAAIAP4FzBOa33NqYEJIUQ5kp2dzeHDh6lRowahoaF4eHiglHJ2WEJUWFprrl69ytGjR6latSpeXl52ua68bjY3D7hLa+0PRAAtgTecG5IQQpQv586do0aNGgQGBuLp6SkJohAOppTCx8eHwMBAzpw5Y7frSpJ4C611qtb6yq1NQCNnxSOEEOVRRkYG/v7+zg5DiErHz8+PrKwsu13PpZJEpdRbSql/KaX+UEpppdShQvq6KaVeUUr9TymVpZQ6opSaqZQq1c7ySqk3lVIZwGngL8DHpbmeEEJUNrm5uXh4eDg7DCEqnSpVqpCTk2O367lUkghMBToDB4HzRfSdBXwEpAIvYRw/OBpYr5Qyey6lVFWlVLVCPu55fbXW72ut/YBwIBo4abenE0KISkJeMQtR9uz9fedqE1fu0Fr/AaCU2gtUs9ZJKdUMY2K4Vmvd75b2NGAu8BSw4pZTFgGRhdz3ISDh1gat9T6lVBKwBOhS3AcRQgghhCjPXKqSmJcg2mAgoIDZ+doXAJlAVL7rRmmtVSGfhALu4w40LsYjCCGEEEJUCC6VJBZDG8AA7Li1UWudBSTdOF4sN147D1FKVVdGLYB3gM32CFgIIYQQojwpr0liXeCs1jrbyrFjQKBSyrOY19TAIOAPjGskrgM2AGNKE6gQQgghRHnkamMSbeUDWEsQAbJu6WPzItg3lr7pamt/pdRIYCRAcHAwCQkJtp5aYpcvXy6T+wghREldvnyZgIAAMjIynB2KS2jevDkNGjRg48aNzg6l2LZt28b06dPZtWsX165do0mTJrz55pv06tXL2aGJQmRlZdktVyivSWImEFTAMa9b+jiM1joa4+xnWrdurTt16uTI2wGQkJBAWdxHCCFKKiEhAS8vL/z8/JwdiktQSuHu7l6iP4+kpCS+/PJLnn32WUJDQ+0fXCE2bdpEz549ufPOO5k4cSIeHh7MmTOHwYMH89tvv9GwYcMyjceeDAYDc+bM4Z///CeHDh2iVq1a9O/fn8mTJ+Pra/sqeufOnWPq1Kl8+eWXHD16FD8/P5o3b87kyZPp0KGDWd9Tp04xceJENmzYwKlTp6hduzZ9+vTh3XffpXr16nZ9Pi8vL+6++267XKu8JonHgXClVFUrr5zrYXwVLVvpCSGEKLeSkpJ499136dSpU5kmiQaDgRdeeIFGjRqxa9cuvL29AQgNDeWxxx7jxx9/LNdJ4iuvvMLcuXPp06cPY8eOZd++fcydO5fExES++eYb3NyKHomXnp5Op06duHz5MsOGDaNJkyZcvHiRPXv2cOzYMbO+p0+fpl27dhw/fpznnnuO5s2bs3fvXj799FO+//57tm3bho+Pj6Met1TKa5K4E+gGtAV+yGtUSnlh3E7veyfFJYQQQpRr27Zt49ChQ8yfP9+UIAKmBdJvbStvUlJS+Pjjj+nbty9r1qwxtYeFhTF69GhWrVrFoEGDirxOVFQUOTk57Nmzhzp16hTad+rUqaSnp7NixQoGDhxoam/fvj2DBg3io48+YsKECSV/KAcqrxNXVmOcaJJ/UskIjGMR48o8IiGEEJXWkSNH6N+/PwEBAfj7+9OrVy8OHjxo0S8jI4MJEybQrl07AgMDqVq1Ko0aNeLNN98kM/PmKKlJkyYxZMgQAB566CGUUiilePbZZ4t1nZLYvn07AA8//LBZ+7fffgtAq1atSnV9Z1q5ciVaa8aMMU8fRowYgY+PD7GxsUVe4/vvv2fr1q2MGzeOOnXqcP369UL/zOPj4/H29uapp54yax8wYABeXl4sXry4ZA9TBlyqkqiUehoIufFlLcBTKZWXXqdrrZcDaK2TlVLzgReVUmuBjUBTjDuufIf5QtpCCCGEw1y4cIEHH3yQI0eO8PzzzxMeHs53333HQw89xNWrV836Hjt2jIULF9KvXz8GDRpElSpV+O6775g+fTqJiYls3mxcda1v376cOHGC6Oho3n77bZo2bQrAHXfcUazrlERiYiLVqlWjceObywQfPXqUmJgY2rRpU+Crb4PBwLlz52y+T82aNW16tWtPO3fuxM3NjbZt25q1e3l5ERERwc6dO4u8Rt4kpAYNGtCrVy/+85//kJubS+PGjfn73/9OVJTZUs1kZ2fj5eVlsRuKm5sb3t7e/PHHH5w9e5bAwMBSPp39uVSSCAwDOuZr+8eNX78Dlt/SPgY4hHGGcQ/gLMZ9lv+utTY4NkwjpVQvoFejRo3K4nZCCFGufTtmDKeTkpwdhpmgiAg6z86/L0PxTJ8+nUOHDhETE2Oq/v3tb39jzJgxzJkzx6xvw4YNOXLkiNne1qNGjeKdd97hvffeY8eOHbRt25a//OUv3HfffURHR9O1a1eLSYu2XqckkpKSiIiIwM3Njf379/Pzzz8zceJEMjIymDVrVoHnHT58mLCwMJvvk5aWVuYTco4fP26qvOZXr149tm/fzrVr1/D0LHgVvd9++w0wVh8bN27M0qVLyc7O5qOPPuLpp5/m+vXrpv8OAJo1a8Zvv/1m+nPNk5SUxPnzxh2IDx8+LEliUbTWnYrRNxeYeePjFFrr9cD61q1bj3BWDEIIIZzryy+/JDg4mGeeecas/Y033rBIEm9NPnJycsjIyCA3N5cuXbrw3nvv8fPPP9uU3NnrOvllZmZy4MABHnnkEdLT07nzzjtNx15++WXuu+++As+tXbs2X3/9tc33ql27ts19L1y4wOxiJPOjR4+mZs2aFu2ZmZlWE0QwVhPz+hSWJOYt7+Tn50d8fLypb58+fWjYsCFvv/02gwcPNlVJx4wZw5dffkn//v2ZPXs2zZs3JyUlhTFjxuDh4VHk62pncqkkUQghRMVV2oqdq/rjjz9o06YN7u7uZu116tSxurzJJ598wmeffUZKSgoGg/mLr7zKki3sdZ1b/frrrxgMBu655x78/f1Zv349hw8fZt26dcyZM4fz58+zdOlSq+d6eXnRpUuXEt23KBcuXODdd9+1uX9UVJTVJNHHx4fTp09bPScrK8vUpzB5E3cGDhxolkzWqFGD3r17s2zZMn777TfTEIEOHTqwatUqRo8eTY8ePQBwd3dn+PDhNGvWjC+++AJ/f3+bn60sSZIohBBClFL+8WZ5tNZmX3/00UeMHTuWbt26MXr0aOrWrYunpyfHjh3j2WeftUj2CmKv6+SXdGM4QKtWrahRowY9e/YEjK/Pu3btSmxsLAsXLjR7zZ0nNzeXM2fO2HyvWrVqWSTWBQkNDbX4syyJunXrkpqaSnZ2tkVF8dixYwQGBhZaRQSoX78+YL0SmjfTOX+S/uSTT9K3b1+Sk5PJyMjgzjvvJCgoiLZt21KlShVcddiaJIlCCCFEKTRs2JD9+/eTm5trlvScOHGCixcvmvVdvnw5oaGh/Oc//zGbtLFp0yaL6xaUeBb3OsWRlJSEt7e3qQp2K4PBQFBQkNUEEYwzvF19TGKbNm3YsmULO3bsMFvwOisri6SkJB588MEir9G2bVs+++wzjh49anEsry0oyHK/D3d3d7MxiSdPniQxMZGOHTvKOolCCCFERfTYY4/x/vvvs2zZMrMJCx988IFFX3d3d5RSZlWxnJwc3n//fYu+1apVA7A6Y7g41ymOxMREDAYDaWlpZtWt7du3k5CQwNixYws815FjEu1lwIABTJ06ldmzZ5sliQsWLCAzM5PIyEiz/tevX+fgwYP4+PjQoEEDAB5//HFefvllYmNjmTBhgunv6cSJE3z55Zc0bty4yMqgwWBg9OjR5ObmMn78eDs/pf1IkiiEEEKUwrhx41ixYgUjRoxg9+7dNGvWjISEBH788UeLGatPPPEEb731Fo8++ih9+/bl0qVLrFixwmp1rk2bNri5uTFlyhTOnz+Pr68vYWFhtGvXrljXAePr2vT09EJf2ebm5rJ3716ys7Pp0KEDo0aNIjg4mOTkZGJiYmjZsiXvvPNOgec7ckyivbRo0YJRo0Yxb948+vbtS/fu3U07rnTs2NFiIe1jx47RtGlTOnbsaNoPuUaNGsyYMYPnnnuOe++9l6FDh3Lt2jU+/fRTrl27xrx588yucfnyZdq2bUufPn0ICwvj4sWLrFy5kt27dzNlyhQeeuihsnr84tNay6eUn3vuuUeXhfj4+DK5jxBClFR8fLxOTU11dhhlLj09Xffr10/7+fnpatWq6Z49e+rff/9dh4SE6I4dO5r65eTk6KlTp+o77rhDe3p66gYNGujXX39dp6amakBPnDjR7LpLlizRTZs21R4eHhrQgwcPLtF1brvtNl23bt1Cn2Hv3r0a0AMHDtQRERHa09NTV69eXUdEROhp06bpixcv2uFPyvlycnL0jBkzdJMmTbSnp6euW7eufuWVV3RGRoZF37S0NA2Y/R3mWbNmjW7Xrp328fHR1apV0127dtVbt2616Jedna0HDBigQ0NDddWqVXWNGjV0t27d9KZNmxzxeDZ9/wG7tA35jdJ2GAhaWd2yTuKIAwcOOPx+CQkJFmtlCSGEK0lISCA4ONjqmDbhHHv27KFly5Zm6zhaExcXR1RUFNu2baN9+/ZlGKGwp3379hX5/aeU2q21bl3UtcrrtnwuQWu9Xms9MiAgwNmhCCGEEFZt3ryZli1bMnjw4EL7JSUloZSiRYsWZRSZcHWSJAohhBAV2Ouvv05SUlKRW+AlJiYSEhKCn59fGUUmXJ0kiUIIIYTg119/lSqiMCOzm4UQQghRrIWwReUglUQhhBBCCGFBkkQhhBBCCGFBkkQhhBBCCGFBksRSUEr1UkpF59+bUwghhBCivJMksRRknUQhhBBCVFSSJAohhBBCCAuSJAohhBBCCAuSJAohhBBCCAuSJAohhBBCCAuSJAohhBBCCAuSJAohhBBCCAuSJAohhBBCCAuSJAohhBCi1L7//nu6du2Kv78/Xl5eRERE8MUXXzg7LFEKkiSWguy4IoQQ4lbnz5/Hy8sLpRSxsbFW+yQkJKCUMvt4eXnRsGFDhgwZwr59+wrs/+KLL1q95unTp/H09EQpRadOnSyOT5s2jSeffJKGDRuilCI0NLS0j2pm06ZNdO7cmePHjzNx4kSmT59ORkYG/fv3548//rDrvcqSwWBg1qxZ3HXXXXh5eXH77bczduxYrly5YvM1Ll++zNSpU2nRogV+fn4EBgbSvn17lixZgtbaor+j/66Ko4rT7lwBaK3XA+tbt249wtmxCCGEcL64uDiuXbtGWFgYixYtIioqqsC+AwcOpHv37gBcvXqVPXv2sHDhQtasWUNycjIhISFm/b28vFixYgUzZ86katWqZseWL1+O1poqVaz/b/3tt9+mZs2atGrVigsXLpTyKc0ZDAZeeOEFGjVqxK5du/D29gYgNDSUxx57jB9//JGGDRva9Z5l5ZVXXmHu3Ln06dOHsWPHsm/fPubOnUtiYiLffPMNbm6F19oMBgOPPvoo27dvZ/Dgwbz00ktkZmaycuVK0z8IPvjgA7NzHPl3VWxaa/mU8nPPPffoshAfH18m9xFCiJKKj4/Xqampzg7DaSIiInTnzp31nDlztFJK//777xZ94uPjNaA//PBDi2Nz5szRgP7oo48s+g8cOFADevXq1RbnNWvWTPfu3Vv7+vrqjh07Whw/ePCgWd+QkJCSPaAV33//vQb0/Pnzzdo3btyoAb1mzRq73ass7d27VyuldN++fc3a586dqwEdFxdX5DW2b9+uAT1mzBiz9uzsbB0WFqYDAgIszint35Ut33/ALm1DfiOvm4UQQgg7+OWXX0hKSmLw4MFERkbi4eHB4sWLi3WNunXrAuDp6WlxrFWrVrRs2dLimjt27CAlJYUhQ4YUeF1HVvK2b98OwMMPP2zW/u233wLGuMujlStXorVmzJgxZu0jRozAx8enwOEEt7p06RJw8+81j6enJ4GBgfj6+lqc40pVV3ndLIQQQtjBokWL8PX1pV+/fvj6+tKjRw+WLl3K5MmTrb6WzMzM5OzZs4DxdfPevXsZP348gYGB9OvXz+o9hgwZwquvvsrRo0epX78+ADExMQQFBdGzZ0/HPVwhEhMTqVatGo0bNza1HT16lJiYGNq0aVPgmDqDwcC5c+dsvk/NmjWLfL1rTzt37sTNzY22bduatedNytm5c2eR12jbti3Vq1dn+vTphIaG0q5dO65evcqSJUvYvXs3n332maPCtwtJEoUQQpSN3WPgfJKzozBXIwLumV3qy2RlZbFy5UqeeOIJU3Vo8ODBfPHFF2zevJlHH33U4pyJEycyceJEs7bw8HB++OEHateubfU+UVFRjBs3jmXLlvH2229z9epVVq1axfDhwwscj+hoSUlJRERE4Obmxv79+/n555+ZOHEiGRkZzJo1q8DzDh8+TFhYmM33SUtLK9NJHMePHycwMNBi/CdAvXr12L59O9euXbNa9c1To0YNvvrqK4YPH07//v1N7X5+fqxZs4bHH3/cIbHbiySJQgghRCmtXbuW8+fPM3jwYFNbjx49CAoKIiYmxmqSOHLkSJ588knAmGSmpqYyc+ZMunfvTnx8vMXEFYDbbruN3r17s2TJEt5++23Wrl3LxYsXGTp0qOMerhCZmZkcOHCARx55hPT0dO68807TsZdffpn77ruvwHNr167N119/bfO9Ckqc87tw4QKzZ9ue+I8ePZqaNWtatGdmZlpNEMFYTczrU1iSCFCtWjWaN29O7969ad++PefOnWP+/PkMGjSIdevW0bVrV5tjLWuSJAohhCgbdqjYuapFixZRq1Yt6tevz++//25q79q1K//61784e/YsgYGBZuc0btyYLl26mL7u2bMnHTt25N577+WNN95g1apVVu81ZMgQevTowdatW4mJiaFt27aEh4c75sGK8Ouvv2IwGLjnnnvw9/dn/fr1HD58mHXr1jFnzhzOnz/P0qVLrZ7r5eVl9vz2cuHCBd59912b+0dFRVlNEn18fDh9+rTVc7Kyskx9CpOcnEz79u2ZNWsWzz//vKl94MCBNG/enBEjRnDw4EHc3d1tjrcsSZIohBBClEJaWhrx8fForWnSpInVPrGxsRYTIKxp164dAQEBpkkf1jzyyCPUq1ePd999l/j4eD799NMSx15aSUnG4QOtWrWiRo0apnGRf/vb3+jatSuxsbEsXLgQDw8Pi3Nzc3M5c+aMzfeqVauWTclUaGio1fUHi6tu3bqkpqaSnZ1tUVE8duwYgYGBRVYRZ82aRVZWlqlinMfHx4cePXowb948Dh06xB133FHqeB1BksRyIDUujq3jx3Pp8GH2N2jAA1OmEB4Z6eywhBBCAIsXL0ZrzYIFC6hevbrF8QkTJrBo0SKbkkSAnJwcsrOzCzzu7u7OM888w7Rp0/D29uapp54qceyllZSUhLe3N02bNrU4ZjAYCAoKspogAhw5csSlxyS2adOGLVu2sGPHDjp06GBqz8rKIikpiQcffLDIaxw7dgwwJsT55eTkmP3qiiRJdHGpcXFsGTmSnMxMAC6lp7Nl5EgASRSFEMLJDAYDS5YsoUWLFgwfPtxqn5SUFCZNmsTOnTtp06ZNodf7+uuvuXLlCvfff3+h/Z5//nk8PT1p2LAhAQEBJY6/tBITEzEYDKSlpdGoUSNT+/bt20lISGDs2LEFnuuoMYn2MmDAAKZOncrs2bPNksQFCxaQmZlJZL7/B1+/fp2DBw/i4+NDgwYNAONEpC1btrBkyRLGjRtn6nvhwgXWrVtHjRo1XLaKCJIkuryt48ebEsQ8OZmZbB0/XpJEIYRwsi1btnDkyBGGDRtWYJ9+/foxadIkFi1aZJYk/vLLL6a19rKzs0lJSWHBggV4eHjw3nvvFXrfBg0aMGnSJJvjXL58Oenp6QCcOXOGa9eume4REhLC008/beobGhpKenp6ka9sc3Nz2bt3L9nZ2XTo0IFRo0YRHBxMcnIyMTExtGzZknfeeafA8x01JtFeWrRowahRo5g3bx59+/ale/fuph1XOnbsyKBBg8z6Hzt2jKZNm9KxY0cSEhIAGDNmDMuWLePNN98kOTmZ+++/n3PnzrFgwQJOnDjB/PnzLWalF+fvyuFsWXFbPtY/QC8gulGjRkWubl5SHyqlPwTLj1IOu6cQQpRUZdtx5YknntCA3rNnT6H9mjRpogMCAnRmZqZpB5VbP25ubrpWrVq6T58+eseOHWbnFrZDS34F7bjSsWNHi3vmffL3v+2223TdunWLvNfevXtNO8FERERoT09PXb16dR0REaGnTZumL168WOQ1XF1OTo6eMWOGbtKkifb09NR169bVr7zyis7IyLDom5aWZvXP8/fff9fPPPOMrlevnq5SpYr28/PTHTp0KHAnmuL8XVljzx1XlLbD4M7KrnXr1nrXrl0OuXZ0aCiXbvyL4lb+ISGMPHTIIfcUQoiSSkhIIDg42OoYNeH69uzZQ8uWLYmJiSl0Bxcw7lMdFRXFtm3baN++fRlFKIqyb9++Ir//lFK7tdati7qWbMvn4h6YMoUq+abYV/Hx4YEpU5wUkRBCiIpq8+bNtGzZ0my9x4IkJSWhlKJFixZlEJlwBhmT6OLCvccSviAzX2smeI0FZEyiEEII+3n99dd5/fXXbeqbmJhISEgIfn5+Do5KOItUEl1d1qnitQshhBBl4Ndff5UqYgUnlUQhhBBCFFtxFsIW5ZNUEsuJlKPhNB+XTMpR52y9VBmkxsURHRrKDDc3okNDSY2Lc3ZIQgghhNNIklgOXMnyofv0DaQeb0qPDzdwJavwvSJF8eUtWn4pPR20Ni1aLomiEEKIykqSxHJgaHQMpy8FobU7py4GM2zBImeHVOEUtmi5EEIIURlJkujiYraNZkNSD7KuG6uHWde9WZ/Yi5hto50cWcVy6fDhYrULIYQQFZ0kiS7urf+bw5XsamZtmdm+vPV/c5wUUcXkf2OfTVvbhRBCiIpOkkQXN20a+Pqat/n4wPvvOyeeikoWLRdCCCHMSZLo4oYOhR49wMvL+LWXF/TqBUXsliSKKTwykm7R0fiHhIBS+IeE0C06mvBIWbBcCCFE5STrJJYDMTEQHg5HjmiCgxWLZN6KQ4RHRkpSKIQQQtwglcRSUEr1UkpFX7x40aH38fWFjRshJCSTDRssXz8LIYQQQtibJImloLVer7UeGRAQ4PB7NWsGixfvpFkzh99KCCGEEEKSRCGEEEIIYUmSRCGEEEIIYUGSRCGEEEKU2vfff0/Xrl3x9/fHy8uLiIgIvvjiC2eHJUpBkkQhhBDCTs6fP4+XlxdKKWJjY632SUhIQCll9vHy8qJhw4YMGTKEffv2Fdj/xRdftHrN06dP4+npiVKKTp06mR3bv38/f//737n33nupVasWfn5+REREMGXKFK5cuWKX5960aROdO3fm+PHjTJw4kenTp5ORkUH//v35448/7HIPZzAYDMyaNYu77roLLy8vbr/9dsaOHVusP7fLly8zdepUWrRogZ+fH4GBgbRv354lS5agtTbr+9tvvxEZGUnTpk0JCAjAx8eHu+66i1dffZUTJ07Y+/GKJEvgCFFBpcbFsXX8eC4dPox/gwY8MGWKLPEjhIPFxcVx7do1wsLCWLRoEVFRUQX2HThwIN27dwfg6tWr7Nmzh4ULF7JmzRqSk5MJCQkx6+/l5cWKFSuYOXMmVatWNTu2fPlytNZUqWL5v/WYmBjmz59P7969iYyMxMPDg/j4eCZMmMDnn3/OTz/9hLe3d4mf2WAw8MILL9CoUSN27dplulZoaCiPPfYYP/74Iw0bNizx9Z3plVdeYe7cufTp04exY8eyb98+5s6dS2JiIt988w1uboXX2gwGA48++ijbt29n8ODBvPTSS2RmZrJy5UrTPwg++OADU/+jR49y4sQJ+vTpQ/369alSpQrJyclER0ezatUqkpKSCAoKcvRj36S1lk8pP/fcc48uC/Hx8WVyH1H+pcTG6lk+PvpDMH1m+fjolNhYZ4cmKrj4+Hidmprq7DCcJiIiQnfu3FnPmTNHK6X077//btEnPj5eA/rDDz+0ODZnzhwN6I8++sii/8CBAzWgV69ebXFes2bNdO/evbWvr6/u2LGj2bGdO3fqCxcuWJwzfvx4DeiPP/64BE960/fff68BPX/+fLP2jRs3akCvWbOmVNd3lr1792qllO7bt69Z+9y5czWg4+LiirzG9u3bNaDHjBlj1p6dna3DwsJ0QECATbF8/vnnGtAffPBBkX1t+f4Ddmkb8ht53SxEBbR1/HhyMjPN2nIyM9k6fryTIhKi4vvll19ISkpi8ODBpord4sWLi3WNunXrAuDp6WlxrFWrVrRs2dLimjt27CAlJYUhBWzF1bp1a6wt1TZgwAAA9u7dW6wY89u+fTsADz/8sFn7t99+a4q7PFq5ciVaa8aMGWPWPmLECHx8fAocTnDbIrKSAAAgAElEQVSrS5cuATf/XvN4enoSGBiIr40LH+dVlc+fP29Tf3uR181CVECXDh8uVrsQovQWLVqEr68v/fr1w9fXlx49erB06VImT55s9bVkZmYmZ8+eBYyvm/fu3cv48eMJDAykX79+Vu8xZMgQXn31VY4ePUr9+vUB4+vkoKAgevbsWax4jx49CkBwcHCxzssvMTGRatWq0bhxY7Nrx8TE0KZNG0JDQ62eZzAYOHfunM33qVmzZpGvd+1p586duLm50bZtW7P2vEk5O3fuLPIabdu2pXr16kyfPp3Q0FDatWvH1atXWbJkCbt37+azzz6zel5WVhaXL18mKyuL1NRU3njjDQDT8ISyIkmiEBWQf4MGXEpPt9ouhLOMGQNJSc6OwlxEBMyeXfrrZGVlsXLlSp544glTdWjw4MF88cUXbN68mUcffdTinIkTJzJx4kSztvDwcH744Qdq165t9T5RUVGMGzeOZcuW8fbbb3P16lVWrVrF8OHDrY5HLEhubi6TJ0+mSpUqDBo0qBhPaikpKYmIiAjc3NzYv38/P//8MxMnTiQjI4NZs2YVeN7hw4cJCwuz+T5paWkFJpyOcPz4cQIDAy3GfwLUq1eP7du3c+3aNatV3zw1atTgq6++Yvjw4fTv39/U7ufnx5o1a3j88cetnrdw4UJeeukl09ehoaHExsbSoUOHUjxR8UmSKEQF9MCUKWwZOdLslXMVHx8emDLFiVEJUXGtXbuW8+fPM3jwYFNbjx49CAoKIiYmxmqSOHLkSJ588kkAU8Vo5syZdO/enfj4eIuJKwC33XYbvXv3ZsmSJbz99tusXbuWixcvMnTo0GLFO2bMGH766SemTp3KnXfeWcynvSkzM5MDBw7wyCOPkJ6ebnatl19+mfvuu6/Ac2vXrs3XX39t870KSpzzu3DhArOLkfmPHj2amjVrWrRnZmZaTRDBWE3M61NYkghQrVo1mjdvTu/evWnfvj3nzp1j/vz5DBo0iHXr1tG1a1eLcx5//HHuuusuLl++TGJiIl999RVnzpyx+ZnsxpaBi/KRiSui/EmJjdX/DAnRHyql/xkSIpNWRJmorBNXOnfurGvVqqX379+vDxw4YPpERkZqT09PfebMGVPfwiau/PTTTxrQAwYMKLD/hg0bNKB/+OEH3blzZ922bVtTX2sTV/KbMGGCBvTIkSNL+dQ3J2YsXbpUnzt3Tq9fv17Pnz9fd+vWTQP6mWeeKfU9iistLU0DNn8OHDhg9TrNmzfXQUFBVo89+eSTGtDZ2dmFxrJnzx7t5eWlP/30U7P2K1eu6LCwMB0SEqJzcnKKfKZff/1Ve3p66qlTpxbZ154TV5yeYFWEjySJQthGEteKrzImiX/88YdWShWaiMyaNcvUv7AkUWutAwICdK1atQrsn5OTo+vVq6e7dOmilVL6s88+M/UtKkmcOHGiBvSQIUO0wWAo5ZNr/cknn2hAJycnWxzr0qWLdnNz09euXbN6bk5Ojj5x4oTNH1uSKXvq1q2bdnNz01lZWRbH2rdvrwMDA4u8xpAhQzSgz549a3HsxRdf1IDVGfDWtGvXTtevX7/IfvZMEuV1sxCiTKTGxZm9Ar+Uns6WkSMBZP1GUa4tXrwYrTULFiygevXqFscnTJjAokWLLGbJFiQnJ4fs7OwCj7u7u/PMM88wbdo0vL29eeqpp2y67rvvvsu7777LM888w8KFC1FK2XReYZKSkvD29qZp06YWxwwGA0FBQXh4eFg998iRIy49JrFNmzZs2bKFHTt2mI0FzMrKIikpiQcffLDIaxw7dgwwjgHNLycnx+zXoly9erVYE33sQZJEIUSZKGxZHkkSRXllMBhYsmQJLVq0YPjw4Vb7pKSkMGnSJHbu3EmbNm0Kvd7XX3/NlStXuP/++wvt9/zzz+Pp6UnDhg2tLm+T3+TJk5k0aRJPP/00ixcvttss4cTERAwGA2lpaTRq1MjUvn37dhISEhg7dmyB5zpqTKK9DBgwgKlTpzJ79myzJHHBggVkZmYSme/n1vXr1zl48CA+Pj40uDFJMDw8nC1btrBkyRLGjRtn6nvhwgXWrVtHjRo1uOOOO0ztJ0+etPqc8fHx7N2712I3HUeTJFEIUSYKXJYnPZ1Patem1+efc7sN/zIXwpVs2bKFI0eOMGzYsAL79OvXj0mTJrFo0SKzJPGXX34xrbWXnZ1NSkoKCxYswMPDg/fee6/Q+zZo0IBJkybZFOP8+fOZOHEiDRo0oEuXLqxYscLseHBwsNnkidDQUNLT041j0gqRm5vL3r17yc7OpkOHDowaNYrg4GCSk5OJiYmhZcuWvPPOOwWe7+XlRZcuXWx6Bmdo0aIFo0aNYt68efTt25fu3bubdlzp2LGjxazwY8eO0bRpUzp27EhCQgJgnCC0bNky3nzzTZKTk7n//vs5d+4cCxYs4MSJE8yfP99sVvoLL7zAiRMn6Ny5MyEhIWRlZbF7925WrVqFn58fM2fOLMs/AkkShRBlo6BleQAyT51idceOAPjWqUOvzz+n/gMPlGV4QpTIokWLAOjbt2+BfZo3b06TJk1YtWqV2ZIwK1euZOXKlQC4ublx22230bVrV956660iK47Fkbee3+HDh81mX+fp2LGjWZJ4+fJli8Wfrfnf//7H1atXGThwIPv27eMf//gHPj4+hIaGMmHCBP72t7/h5+dnt+dwhtmzZxMaGkp0dDQbNmwgMDCQl156qcC1L/MLCQlhx44dTJ48mf/+97+sWrUKb29vIiIimDlzpsV/NwMHDmTp0qUsX76cM2fOoJQiJCSE5557jtdff91UoSwrqqh/KYiCKaV6Ab0aNWo04sCBAw6/X0JCQpmXmoWwl/xjEsG4LE+36Gh8g4NZ378/Wfl2E/CrX5+eq1dTr337sg5XlFBCQgLBwcFWx6gJ17dnzx5atmxJTExMgTu45ImLiyMqKopt27bRXr5HXca+ffuK/P5TSu3WWrcu6lqyLV8paK3Xa61H2jIeRIjKLjwykm7R0fiHhIBS+IeE0C06mvDISEK6dOHFc+d4TWue2LyZqje+pzKOHmXl/fczQymiQ0M5/tNPTn4KISq2zZs307JlS6sVx/ySkpJQStGiRYsyiEw4g1QS7aB169Z6165dDr+PVBIdaG1tyDpFytFwBsxdzerRA2hWPxW8gqHvSWdHV2mlbdrE+v79uZaRYdYeEBZGz1WrqJNvuyzhfFJJrDy6dOnCwYMHSUtLc3Yo4hZSSRTC3rJOcSXLh+7TN5B6vCk9PtzAlSwfyDrl7MgqtbC//pXRly7xmtb0+fe/qeLjA8DFtDTi2rVjhlIsbNSIk2XwjzQhhLlff/1VqogVnCSJQtwwNDqG05eC0NqdUxeDGbZgkbNDEre4o0cPxly5YkwYv/oK9xvbZV04eJDYNm2YoRSL7ryTk7t3OzlSISqHM2fO8NVXXzk7DOFAkiQKAcQkPMuGpB5kXTdWqrKue7M+sRcxCc86NzBh1R29evFKVhavac1jX3yB243Fes/v309s69bMUIqY8HBOJyU5OVIhhCi/JEkUAnhr9ftcya5m1paZ7ctbq993UkTCVo0ff5xXr13jNa3pvWYN6sayFOf27WPZ3XczQykWN2/OqV9+cXKkQghRvkiSKAQwbcCb+Fa9bNbmU/UK7z/1hpMiEiXRpG9fxubm8prW9Pr8c1P7nykpLL/nHmYoxQyl2HdjbTohhBAFkyRRCGDoX/9Dj4gNeHkY1/Dz8rhKr7vXM+SRTU6OTJTUnU8+yWta8+jy5XjddpvZsQ2DBpkSxt/+9S8nRSiEEK5NkkQhAPqeJOa/Awiq44NSEFzXm0XfPCXL35RzqXFxfP3cc2T9+aepTeXbJWF9//6mhHH/2rVlHWKFJcurCVH27P19J0miEDf4+sLGjRAeDhs2GL8W5dvW8ePNdngB0AYD/iEhjDUYCH3kEbNjX/XrZ0oYf5dZmyXm7u7O9evXnR2GEJVOTk6O2V7QpSV7Nwtxi2bNYO9eZ0ch7OXS4cMFtiuleGKTcTiB1prPH36YI/Hxpj5fPvaY6fd9/v1v7ujRw7HBViB+fn5cunSJwMBAZ4ciRKWSkZGBl5eX3a4nlUQhRIXl36CBTe1KKQZ8+y2vac1Yg4F6DzxgdvyLnj1NFca0TTJOtSg1a9bk/PnznD17lmvXrsmrZyEcTGtNZmYmZ8+epVatWna7rlQShRAV1gNTprBl5EizV85VfHx4YMqUAs9RSjHwhx8A46vpFe3bc+Lnn03H1zz6qOn3T379NSFdujgg8vKtatWqNGjQgHPnznHo0CFyc3OdHZIQFV7VqlUJDg62ayVRkkQhRIUVHhkJGMcmXjp8GP8GDXhgyhRTe1GUmxuRP/0EGBPG2LZtOXXLji7/6trV9Pten3/OnU8+acfoy7eqVatSp04d6tSp4+xQhBAlJEmiEKJCC4+MtDkpLIxyc+PpG3tEG3JzWXb33ZxNTjYdX9+/P+tv/L778uWER0WV+p5CCOFMMiZRCCGKyc3dnWf37OE1rXk1J8fi+MannzaNYfzfqlVOiFAIIUpPkkQhhCgFN3d3XtOa17Tmlexsi+P/Hjjw5jqMa9Y4IUIhhCgZed0shBB24u7pyWs3ZvLmZGczO98A8q+eeML0+8e//JJGtyyzI4QQrkYqiUII4QBVqlY1VRjHXL1qcfzLxx83VRgPbthg13unxsURHRrKDDc3okNDSY2Ls+v1hRCVg1QShRDCwap4eZkqjNczM5mTbzufL3r2NP2+36ZNhOXbCaY4UuPizJb9uZSezpaRIwHsMoFHCFF5SCVRCCHKkIePj6nCOPryZYvja/76V1OF8ZPg4GJXAa1tRZiTmcnW8eNLFbcQovKRJFEIIZzE09f3ZsJ46ZLF8czTp9kYFcUMpTjy3Xc2XbOwrQiFEKI4JEkUQggX4Onnh39ISIHHV3fqZKowHtu2rcB+tm5FKIQQRZEkUQghXISt1b6VDzxgShiP39gRJs8DU6ZQxcfHrK2orQiFEMIaSRKFEKIU7DmTuMAqYEgIr2nNqD//tDi24r77TAnjyV27CI+MpFt0tLEqqRT+ISF0i46WSStCiGKT2c2loJTqBfRq1KiRs0MRQjiBvWcSPzBlitn1wLwK6F2zpmmWdObZs3xSq5bZ+bFt2ph+/0xiIkEREcWOQQgh8kglsRS01uu11iMDAgKcHYoQwgnsPZO4OFVAn8BA06SXF06etDi+7O67TRXGEzt3ligeIUTlJpVEIYQoIUfMJA6PjCx2FdI3ONhUYbx84gSf1a1rdjyubVvT75/+5ReC7767xPEJISoPqSQKIUQJueJM4mp16pgqjCPS0iyOL2/VylRhPJOcbPUasmOLEAIkSRSiYlpbG1YoUqY3o3n9vaRMbwYrlLFd2I2rzyQOCA01JYzDf//d4vjSv/zFlDCeTU0Fbo6zvJSeDlqbxllKoihE5SNJohAVUdYprmT50H36BlKPN6XHhxu4kuUDWaecHVmFUp5mEle/4w5Twjj0t98sji9p1owZSrExKkp2bBFCADImUYgKa2h0DKcvBaG1O6cuBjNswSJWvTTQ2WFVOCUZQ+hsNZs0MY1h/HPfPhaHhxd5juzYIkTlI5VEISqgmIRn2ZDUg6zrxlehWde9WZ/Yi5iEZ50bmHA5tzVtaqowDi5gjCIAWnPx0KEyi0sI4XySJApRAb21+n2uZFcza8vM9uWt1e87KSJRHtRq3pzusbEW4yzzLAgLM41htEdlUSbICOHaJEkUogKaNuBNfKteNmvzqXqF9596w0kRifLC2jjL9pMnW/SLDgkxJYwZx44V+z4yQUYI16f0jXEpouRat26td+3a5fD7JCQk0KlTJ4ffR1QAa2szYPocvvqlF1nXffDyuMpj96xj1etjoK/lwstC2OrEzz8Td++9hfbxDwnhgSlTCh2rGR0aakwQrZw7Ul5rC+FQSqndWuvWRfWTSqIQFVHfk8T8dwBBdXxQCoLrerPom6ckQRSlVqddO9MYxoFbt1rtcyk9nY1RUSR+8kmB13HEQuRCCPuSJFGICsrXFzZuhPBw2LDB+LUQ9lTv/vt5TWvjq2kr/jtqlOmVdObZs2bHXHEhciGEOUkShajAmjWDvXuNvwrhKLZU/z6pVcuUMGadP+/yC5ELISRJFEIIUUoFVgVDQuj3n/9YtM+rWdO0aLff7bfbtBC5rTOhZca0EPYjSaIQQohSKawqGPbXv5rGMPb56iuLczOOHDHNbm7YvbvV69s6E1pmTAthX5IkCiGEKBVbtye8o1cvU8L42Nq1FteZV7Om6ZV09sWLpvat48fbtFWgrf2EELaRbfmEEEKUWt72hLYu1dW4Tx/T1oC/zJ3Lty+/bHb84+rVi7xG/rGQMmNaCPuSSqIQQginajV6tKnC+OD06Tafl38spMyYFsK+JEkUQgjhMtq+/ropYXzgvfcK7Xvf3/9u9rXMmBbCviRJFEII4ZLuHT/elDA26tPH4vjmYcNMYxhzsrJsHhsJMgtaCFvImEQhRJlKSYEBA2D1alm/Udju8Vsmuvx74ED+t2qV2fHZ3t6m34/JyqJK1aoFXitvFnTeJJe8WdBAoVsJClHZSCVRCFFmrlyB7t0hNRV69DB+LURx9Vy58maF8bHHLI7P9vIyVRhzr1+3OC6zoIWwjSSJQogyM3QonD4NWsOpUzBsmLMjEuXd419+aUoYQx95xOL4LE9PU8JoyMkBZBa0ELaSJFEIUSZiYox7SGdlGb/OyoL1643tQtjDE5s2mRLG+g8+aHH8Iw8PZihl/FeKFTILWghzkiQKIcrEW29Zvl7OzDS2C2FvT333nSlhrN2mTZH9ZRa0EJYkSRRClIlp08DX17zNxwfef9858YjKI2rHDlPCGNi8udU+OZmZbIyKQhdQZbRGZkiLik5mNwshysTQobB5M3z1lfFVs5cX9OoFQ4Y4OzJRmTybnGz6/coOHTi2davZ8ZluN2snYw0GlFJWryMzpEVlIJVEIUTZWFubmEd8CfI9hMJAcLVDLOrmC2trOzsyUUkN/OEHU4Xx9ocesjg+083NNOklf4XRlhnSUmkU5Z0kiUKIspF1Cl+vTDaO60F4vVQ2vN4DX69MyDrl7MiEYMC335oSxrr33Wdx/NaEEYqeIZ1XabyUng5amyqNtiSKklwKVyFJohCiTDWrn8re6S1oVj/V2aEIYdWg7dtNCaNH/oG0YNMM6aIqjQUlgqVJLoWwN0kShRBCiAK8fPmyKWEsyq0zpAurNBaWCJY0uRTCESRJFEIIIWyQlywWlDDmzZCeoVSBay76N2hQaCJY0uRSCEeQJFEIIYQopqISxkvp6RZteZXGwhLBkiaXQjiCJIlCiLLhFVy8diHKiaISxjw5mZmER0YWmgg+MGUKVXx8zNptSS6FcARZJzEfpZQ3kAwEaq2rOzseISqMviedHYEQDndrojjDyhqL1trgZiKYt8Zi3qvnvMQxPDLS2GalQinbCQpHkUqipcnAUWcHIYQQwrlKO0nE1gojQLfoaFOCGB4ZychDh3jNYGDkoUOm9sKqjEI4giSJt1BKtQK6A9OcHYsQQgjnsfckkaISxrwJLwVVGsGYPHaLjsY/JASUwj8kxCy5FMLeXOZ1s1LqLaAVcA8QBqRrrUML6OsGvAw8B4QCZ4DPgb9rra+U8P5VgAXAqJKcL4QQouIobJJIaZOy4rySzp9UhkdGSlIoyozLJInAVOAc8AtQ1FjAWcBo4AtgJtD0xtd3K6W6aK0NeR2VUlUBj0KudVVrnQuMBfZorROUUp1K/BRCCCHKvbKaJFKahFEIR3OlJPEOrfUfAEqpvUA1a52UUs2Al4C1Wut+t7SnAXOBp4AVt5yyCCjsn10PKaWOYKwg3l2qJxBCCFEh+DdoUOaTRMoyYcxbuDv/5BghbuUyYxLzEkQbDAQUMDtf+wIgE4jKd90orbUq5JMAdABqASlKqZPAWsBfKXVSKdW+NM8lhBCi/HH2JJGixjDmjV/MSxyLM8lGFuUWtnKlSqKt2gAGYMetjVrrLKVU0o3jxbUa2HTL1+2BxUAExlfgQgghKpHClqIpa8VdVicv6QOsxuvI8ZaiYimPSWJd4KzWOtvKsWNAe6WUp9b6mq0X1FpfBa7mfa2UOmds1rKwmxB2lpMD334L998Pvr7OjkaIgrniJJGiEsY8hSV9sii3sFV5TBJ9AGsJIkDWLX1sThLzu/EKutDJM0qpkcBIgODgYBISEkp6O5tdvny5TO4jhCPt31+N555rbfp61qxEIiIuOjEiYU/yc6rstI6PB2BX585g5bX0pfR0YyLp5kbr//7X1O4ZFMS1U6cs+nsGBVn9u/vzm284tnAh106fxjMoiHrDh3Nbly72exDhspR2wdlSeRNXrC2Bo5RKBoK01hZ7eSmlPgeeBKoWp5JYWq1bt9a7du1y+H0SEhLo1KmTw+8jhKO98vI1Zs/1NGsbMQLmzQNPzwJOEuWC/Jwqe9GhoVYn2Vjj6edHl08/ZcvIkWavnKv4+FhdczFv/KItfUX5oZTarbVuXVQ/l5m4UgzHgcAbS9vkVw/jq+gySxCFEMV0bjez2lVFxyl+mPWSqXnBAqha1Zgk7t7txPiEKGcKmmTTPTbWou+1jAw2RkWRk5mJcnMrclHuwsYvioqvPCaJOzHG3fbWRqWUF8aJJo4v6QkhSq7G3eB/FwAPBM1DxykyF3sT2ft3AK5fh9atQSkYNw5yc50ZrBCur7CdWLrHxpra89MGg2l28/9WrbJ6bRm/WLmVxyRxNaCBMfnaR2Aciyhz+IVwZcoNeu6DgbnQ0riciLdnFrEDGqPjFP+ZfHPTow8/hCpVICgI9u1zVsBCuD5r+z3nX+oGCq4w/vHvf5uW1Fn3xBOm9oLWhXTkepHCdZR4TKJSyh3j9nlXtNYnSh2IUk8DITe+fAnwxLibChi36Ft+S9+PgRcx7riykZs7rmwDOt+640pZkDGJQpTS5TRI6A6X/mdqupjpT2RsIhu+a2jWdcoUeOstq4UR4QLk55TrKGison9ICCMPHQJAa81Mt4LrRV6BgeRkZlqMSWw2eDBpGzc6fXkgUTIOHZOolHoQ43Izu4ADN9pClFIWk0mKYRjwjxufIIyzi/O+Hpav7xjgNaAZMB/jLisfAz3LMkFUSvVSSkVfvCgzM4UolWphN6qLBmg9D4AAn0v8e+Qd6DjF6tdeNHUdPx7c3KBJE7jx/7kipaRA8+bGX4WoLGx5VayUMi3a/ejy5RZ9s86eNUsQ/UNCaDZ4MClLl8pi3JVASV83zwPGaq2rA9dvtNXBuAB1iWitOxWyK0qnfH1ztdYztdZ3aq2raq3raa1f1VpfLun9Sxjzeq31yICAgLK8rRAVl1LQZBQM0tDnONzWDoD+d89HxylOf1qL+yOMlZEDByAszHjKJ59YXQEEgCtXoHt3SE2FHj2MXwtRGRT3VfG2CROKvOal9HR+/fRTmcxSSZQ0Sbxda533T4a8H81JQKvShySEEIB3HXjkJ2PCeJ+xwlHL/yxbXw9FxykWPP+yqeuoUcbqYps2cDLfEvhDh8Lp08Yk8tQpGJb/vYQQFVRxtxYs7WQUmcxS8ZQ0STyglMpLCBUYt8XDOHFECCHsKyzKmCz2+xPqPALA8A5z0XGKI3Pr0/yO4wDs2gV16hiri7GxEBMDGzZA1o1l9rOyYP16Y7sQFV1hs56tKbDyGBJieiU9trDlBrTm13/+0x6hCxdR0iRxGvAvpVQHblQSlVKdMK5hKIQQjlG1Jjy0yZgwPrgOgPq3HSN5cj0MsYoZka+Zuj79tLFqmP/1cmamceKLEJWBtVnPBbGl8qjc3OgeG2vRL8/Xzz9vmiV9+Ntv7fMQwmlKlCRqrb8ApgP/BgKUUt8C6zAmj0II4Xj1exuTxSczoEF/lIKx3Wei4xT7Zzbm9uA/rZ7mU/UK7z85uoyDFcL12Vp5tNbv0WXLLK73+cMPmxLGCwcPltVjCDsq1bZ8Sil/4CHgNmCX1nqPvQIrT2QJHCFcxMlv4duHTV8aDIoWbyaTeqyZWbd+bf6PpfFP4Otb1gFWfPJzShhyc1n6l7/wZ2qqxbGQrl15ZOFCWWfRyRy9BM4LSqn2QK7Wep3WOqayJohCCBdSu7OxujggCxqNxM1Ns2NyWxrcdggwoJRxhaw1O5+gWjXj2MX4eKdGLESFkRoXR3RoKB95eHD9yhW6x8by6vXr3DVwoKlP+tdfEx0Swsc1alCaIpUoGyUdk9gT+Ba4qJTar5T6XCn1tlKqi1LKw47xuTRZJ1EIF+VeFdoaB9D7emWycVwPmtVLJfn9Flxb6sFL3eaaunbubEwWhw+H7GxnBSxE+ZZ/d5e8tRP/t3o1PVes4DWteTUnh4dmzQKg7n33oWRFfJdXotfNSqnpQACwBOOr5r5AFJAG1AYmaq1n2y9M1yavm4VwUSsK/5/Q9v33cf+7283a3N3hxx+Ny+mI4pOfU5WTLbu72Co1Lo6t48fLbi4O5NDXzcBw4AWt9Y9a639rrYdi3AVlIdABGKaUer6E1xZCiDLRvsmP6DhF5mJvnn7AOPA+NxfatjVWF197DXJynBykEOWALbu72KKgiqTs5uIcJU0SzwKN87UtBEbfGJsYhXFvZSGEcB6vAnYK9Qo2jl18KgfcvfH2zGLZC4PRcYrNb3QzdZs5Ezw8IDDQuGOLEMK64u7uUpCt48e79G4ueeMuZ7i5ER0aWuGT15ImibOAtUqpe29pawH4AmitfwVk6pIQwrn6njQmg/k/fW9sy+LmDgMyjW2d/wtAt798jY5TXFzgT69WXwHw55/QrJmxuviPf4ChzHaIF6J8KGiNxbDu3YuVVNmrIukIlbHKWdJ1Ej8FYoAtSqn/KaV+BH4APgVQSjUCynQfZSGEKJW8mdEDDeDfFFP3tJYAACAASURBVH+fDL4a+xg6TvF/L/czdfv7343jFhs2hLQ0J8YrhAuxtnZis8GDSVm6tFhJlb0qko7g6lVORyhpJRGt9UygPvAGxlfND2mt8/6kagLvlD48IYQoY0pBz1RjwvjX3QD0a7sWHac481kgD971HWBMEBs2NHafO9e4N7QQlVn+3V3SNm4sdlJV3P2mrXHUK2FXrnI6is1JolLqolLqe6XUXKXUMKXUPUD2jXUSF2mtf87rq7XeobVe5JCIhRCirNRsdfMVdZ1HCPT7k+/e6YSOU8SMHGLq9vLL4OZmTBgTE50YrxAupCRJVXH3m87Pka+EXbnK6SjFqSR2AD7D+Br5cWAbkKGU2quUilVKva6U6uqIIF2VrJMoRCWSt2d0rwMADOm4BB2nOPpxPbNurVoZk8WuleqnoRCWSppU5a9IAjZXBh35StgeVc7yxuYkUWu9R2u9Qmv9NpACLALuB8YDVYD3gX86JEoXpbVer7UeGRAQ4OxQhBBlxa/RzepiwyHUq3kcHacwxCpCAg+Zun3zjTFZVAp+/9154QrhLPZ6dVycyqAjXwmXtspZHlUp4XkjgUCtdd4cv3VKqe3AXfYJSwghyoF7Y4yfjIOo9Y04NCcMgC17uvLIB1tM3RrfWDDs2Wdh8WInxCmEE+QlT6VZGLuwyqC16/g3aGB9UW87vRIOj4ys0ElhfiWduHIaaJav7Z/Ao6ULRwghypG1tY27uqxvZNact4xOznJ3s/YlS25WF0+eLMM4hXCS0rw6huJXBivjK2FHKmmSOANYpZS6deOqCKB66UMSQohyIutUwcf6nsLdzYCOU+g4RdzfBpkdrlPHmCxOmODgGIVwESWZVFLccY2V8ZWwI5V0ncSFQDTwjVLqoFLqJ+A7KtmYRCGEKJBX0M2xi2GDGXT/SnSc4upiL7NuU6bcrC5euuSkWIUoAyWZVFKSymD+6qUkiCVXmnUS5wC3A68ACzCuk/imvQITQogK474lxmSx/2W8PLNN1cU5T4826xYQYEwWP/7YOWEK4UjOWBKnpCrb9nsFKXGSqJRyB4KAnTfWSfzRfmEJIUT5kXI0nObjkkk5Gl54xyq+N6uLzScy+q8fo+MUFxaYr5AwevTN6uKVKw4MXIgyZK8lccoiQaxs2+8VpERJolLqQeAYsAs4cKMtRCkVbMfYhBDC5V3J8qH79A2kHm9Kjw83cCXLp+iTAP4yyZgsPnWdAJ9LpuriW72nmnWrVs2YLL4je1iJcq68TCqpjNvvFaSklcR5wFitdXXg+o22OoAs7iCEqDy8ghkaHcPpS0Fo7c6pi8EMW7AIvIrx72W3Kjeri+0WMnXAeHScYv/Mxmbd3nvvZnUxJ8fOzyFEGSgvk0oq4/Z7BSlpkni71jqv7pq3Y2kS0Kr0IZUfsuOKEJVbzIWTbNg7gKzrxupI1nVv1u95ipgLJVzf5o5hxmRxoIHGtX83VRf9vc1/xnh4GJPF6OjSPoEQZcvaq2NXG/9n62txV4vbEUqaJB5QSuUlhApAa50F2PiepWKQHVeEqNzeestyzGBmprG9VJQyqy5eXFgdHaf4eXJbs27PPXezuqh1AdcSwoW54vg/W16Lu2LcjlDSJHEa8C+lVAduVBKVUp2A43aKSwghXN60aeDra97m4wPvv1+666akQPPmxl//v717j7OyLPc//rk4yDhDYKUMoQJubZsM5gk7bUtKtzuZJt3un0KiCUNOVp46kWQHf7YN0iwzrRxkkhRl0h8l41C4K7F0W4lnhmrvVBBCwSPIDIvj9fvjWWtmnQ+zjrPm+3691ms5z3M/z3MvDsuL677v6+7NLp7rvOfwR3uzi/GGDAmCxY6O/J4tUkqlnv+XTfYvm2HxbPsd/7zl37mv7+/2ANCvbfnc/RdmNga4DxhpZr8DTgAuTX+liEj1aG6GlSth+XIIhaCmBpqaYPbs/t+zuxumTYMNG6CxMfifSW8gem44XbhpZW+geN/jjTRdf1/v9R//eN+9lF2USlfK+X+R7F8kuItk/4CEeZGZtt/Lpt/xz3tl/ctcMe/dvME+GhuHxP7drlD51Em8haBO4lnAHcAH3X1xoTomIjIQtLXBmDFBFq++HhYtyuNmy8bSfEo7W17swR02b9rBnFOXBtv/RRv3b73ZxY8d34kvMfbdkZhdjAxFr1iRR59Eiqi/ZXH6o5BZy2z6Hf+8dtp40w/EfQibN8OcOTk/tuQyBolmVm9ml5jZ9Wb2f83sTDOrBXD3be5+r7u3ufvTxe+uiEhlqasLgrBJk6CzM7/MQNuvT6fzycbYhTBPNNH269NTXxQOFm3aE71D0a1zLoxp0tjYFzCKVJJSlsUpZNYym35H3/fPzOKvNLInvHQjFAqmhrS15fzokkobJIbnHP4duIFgZ5WvA8uA9WamoWUREaChAdasCd7zMa99Ad07R8Yc69lZx7z2LCY5vvXY3oDxwo/cii8xdi0entAsEiw+/nh+fRUphFKWxSlk1jKbfkff91csYBdxf7cLscityDJlEr8bfm8GxgMTgRnAs8ANZnZn8bomIjK4zJ9+BXUjtsccqx3RzYIZX8n+JlFD08OH7enNLv7HiffENDvhBGUXpTKUakeVQmctM/U7+nmncwX7Efd3uwCL3IotU5A4Gfi+uy92943u/oK7/9zd3wdcCEw3s4uL300RkerXPPU2Go/tpGZ4MI+pZvgOmo7rYPbJOUz3Dm0GErcKvOfys/ElxtaFoxIuiQSL69fn/xlEKlWpi3lHP+89tpija1cxYnhQCb8Qi9xKwTzN8jczewn4mrvfmuL8EuAYd59cpP4NCFOmTPHVq1cX/TmrVq1i6tSpRX+OiJTJnUZ3qJZJc7vY8Op4xh/4Al3faaCupqdvZXMu93jtUMa/fUPsPXa/CXcHgeKkL3fxl02J+01PnAjPP9+/j6DvKZHkuruDucsbNsD48ZR1dbOZPebuUzK1y5RJfACYluH84bl0TEREUqipp66mhxVzG5l08Fo6v9wYBHe5bPMHybcKjBj+lt65i2t/eTO+xNj4w4Njrl+3ri+7+Ic/FOBziUhBF7mVSqY6ia3AXWZ2mbv/IMn5iaiAtohIYZwVbOfXAKyZC5B7xd22VbOSr5BeNYvmc+Man3gznHgzB+/bjb8tmJxoM2Mzlh/6UPB+yCFBBkRE+i+yyG2gyJRJ/C3wNuB7ZrbKzM4zs6PM7J1m9ingMuD6oveyQmnvZhGpNP1aIT1keG920buu45z3trPf0FBMk40b+7KLoVCK+4hIVckUJH4LWAFsAD4E/AxYA/wVuAX4G/CGmU02s37t3jKQae9mEak088/7dvIV0udnt4Kz7eaNdD7ZyK69NSnb7L9/ECxG7+4iUi7ZbLUn/ZN24UpMQ7MDgOOAY8Ov44B3EQxZO7CbIGh82t3PL0pvK5QWrohIJZk+PXarwDPOgKVLs7u2fvRmtmxLnAM5ZtRmbpp1MefceHfS6/buhb/8BT72sW7uu68u75qRItmI3/oOgrI2xVy1XA0KtXCll7u/4e4PuPv33f0Cd383MBKYQlAOZyHwJtDU306LiEj+8tkqMF2txrN/cHdv3cV4Q4fC5Mmwbl0tjY3BSk6RYivkVnuSqN97NwO4+y53fzy8Ld8l7n6Sux9QqM6JiEju8llFmbFWY2Tu4gu/xJcY35kxN+4Oxvr1MHIkiXtOixRYPlvtaZg6s7yCRBERqUz5bBXY1tLMmFFbMPZRP3oziy6ck9jo0DPhXGfuXdey6MLZ1O63PaGJ/cdLmFX+/rQycPV3q73IMPW29evBnW3r13N/S4sCxTgKEkVEJEbSWo1pzGtfQM+ukSnPz5mjLQClOPq71Z6GqbOjIFFERPqEC3c3HLKWNdceTcMha2OOJ5NqHmPrpxIzkJFg8be/LVyXZXCJHiZ+6Morabjggpy32stnmHowUZAoIiJ9znqpd95hzCtc6DuZVPMYL/xwG/7GX/ElxvETH4u55tRTlV2U3CUbJu5avJiTrrmGL+3bR8u6dTEBYqp5h/0dph5sFCSKiEjeUs5jHP0uONd57PkT8CXGGwsT68pGgsXHHks4JRIjl2HidPMOcxmmHswLXBQkiohIfrLdc/pcZ/SntuI7XklaRmfKFGUXJb1chonTBZSTZs7ktNbWjMPUg32Bi4JEERHJT3iIumFuFzfd8QoNc7vSD1HXHBiU0XFY850Gjqj/n4QmkWBx48Yi910GlFyGiTMFlJNmzqRl3bqkw9QR/VngUk2ZRwWJIiJSFt3dMO3mLp7d8s9MmLCP7YsSCzoeeqiyi9Inl2HiQsw7zHWBS7VlHhUkiohIWTQ3w5Yt4A6bNw9hzv3duIMvP5IHv/ahhPaRYPHNN8vQWakI2Q4TQ//L40TLNdCsttI6ChLzYGZNZta6devWcndFRGRAaWsLdoMJhYKfQyHo6AgX3m76Gx/61u/xfZ507uKoUUGwODpxDYwMAtkME0faZRtQppJroFltpXXM3cvdhwFvypQpvnr16qI/Z9WqVUydOrXozxER6a9sv6fq64MsYrwxY2Dz5iQXPPwJbmkbxUVttyS93549wf7RJbNsLISSdLSmPm25IBl41i5ZwkNXXsm2F15g1PjxnHTNNSkDzdaJE4Oh5jijJkygZd26Ivc0e2b2mLtPydROmUQRESm5+fMT95SurYUFC1Jc8C938elFtwTD0Umyi8OGBdnFj34Uurpg8uTgvWiSBYjpjsuAlW3mEgozxF1JFCSKiEjJNR8wlsbJ7bEFuN+9lNmjx2a+OLwy2p/+FpecdmPMqZUr+wLExsZgcYxIqRRiiLuSKEgUEZHSC21OXoA7l0zc0V/nxpWXpswurl8PI0fCvHkF7HeO+pvVLEk2VIoil8xjpVOQKCIiZZG0AHd/nessWgR1++9JOLVgQfIyOsUOxLq7Ydo0WLs2t6xmf68TKTQFiSIiUjYNh6xlzbVH03DI2rzvNe8LL9O9Y1jaNpFgceHC4gdisSV+YM6cLC5aNpbmU9rZ8mJPcN2mHcw5dWmwUGbZWLjTEl/LshiiF+kHBYkiIlIV5p8zl7oR22OO1Y7opq1lVkLblhZ44YUcA7ho8VsOxh1PW+InjbZfn07nk42EdgeLH0K796fjiSbafn26FsuUioLxXgoSRUSkKjRPvY3GYztjF8Mc18HskxcH8xYdXm89IOG6UAja2+Eb38jhYeGtCBNe4fI38+YlZid7ejLPj5zXvoDunSNjr9tZx7z2VMu+peAUjPdSkCgiIqWXIRMXI4fMTtLFMFEOuPANxox+Oemjv/Wt3LcATDWvMecSP5Hrpl+RNBu6YMZXsu+USIEoSBQRkdLLkImLEc7gdG2cxOS5z9C1cVLM8WjZLIZJNiwdLxIsbtiQuk26BSbNzcGxmprg55oaaGqC2bPTPjZtNjQi4ddBym7tkiW0TpzId4cMoXXixAG7V3M8BYkiIlLxukO1TLu2k7WbjqLxuk66Q7Up22ZaDJMsEJv+vqVJy+iMHx8Ei0eMS9xFI+3ClGVjafu3OsbUrQuymiPXsei0uqwWoKTLhuby6yClsXbJEu5vaQl2WnFn2/r13N/SUhWBooJEERGpeM2tbWzZNgb3oWzeWs+chYsSG+UwhJ0qEIvMXfzz1SfGtH/2xQm92cXuN/dmXpgS2pw8qxnanH7OW0198utq6qGmPvmvQ6rPLSXx0JVXsqcnNmO9p6eHh668skw9Kpz0tQJERETKrG3VrOQrflfNovncqIY57JkcCcSm39hO+6XTE4alT/z6o/jXAd+HDYnNp4wclXyT6MjClObmvmORrGbWwp+hAVgzF6BvsmNbG3SugdDu4OfQ7v3peHoGbb99hOZQkomU2ke6f2rqU+/LncS2F17I6fhAokyiiIhUtGKt+M2qRqMN6c0u/mb+59LeL5uFKdFynVuYcsX0HV9NfsEgXI1bELnMlwVGjR+f0/GBREGiiIhUtPnnfTv5it/zr+nfDXMYlo5evXzKFTf3BozJ9PRA6A+fyaoL/ZlbmHLFtFY+l9VJ11zDsNpaXmIS3+UZXmISw2prOemafv75rCDmqf60S9amTJniq1evLvpzVq1axdSpU4v+HBGR/irW99T06bB8eTD3r6YGzjgDli4t+GNidHfDpEnBCufx44NAMTpI6+6GCQdv59WtI5Nen2whTMT0G5ey/PEmQrtrqRm+gzNOuJell3wiyFilkfTX4eNpavakuV9XV3C/9nZoaEj72MFr2djUQ89RmcXVi5byb58+idf3voO3Dn2Rlbc8xJQ5M0rY0dyY2WPuPiVTO2USRUSk4rW1wZgxwcKR+npYlGTdSqFl2lavrg4efHgkDQ2wZk3i9TbTsZnO8sea+g7W1NP28KXJ51g+fGnGPhXq1yFd+Z5i72k9oGRZWPu6+2fQM/wQnKH0DD+E7/5X5QaIuVCQmAczazKz1q1bt5a7KyIiVa2uDlasCDJ7nZ2Jw66Flu22epEAsaGhb2V06/zYkaUzvre8N2Bk6grm3fOD5HMs7/lBxn4V6tchVQCcLniU5Pq7BeNAoCAxD+7e4e4to0ePLndXRESqXnRAVmz93VYP4MIrpuAOe/cmnrO3H8+WLfTWaIzIZdFLwq9DLrvXkD6oyZQ9jVZtGcf+FsTO589KpVOQKCIiEqe/2+pFGzKkL7t45Re3xJyLDDUD1IzYm9VuLCnluBp33hdeThrUXP65bVlnxKot45iqIHY2CvFnpVIpSBQREYnT3231UvnP747BHXbsSDwX2jmU9nbYdPO4/nc4B8m2Jawd0Q141hmxXDKOA0GqgtjZKPSflUqiIFFERCSJYiyWqanpyy6eftrOmHMHX7wJM/jUh2+FZ4s3oS3V/tA3nH95VhmxapyDl6rwdfcbKS6IG8ovx8KqUlCQKCIikkSxF8usWDkCd9i2Lfb4olWfwo5oxgx23LZ/YR8almxbwuapt2WVEavGOXipCl8vuXZCVkP5pV5YVSoKEkVEROItGwt3Gg1PGWu+GrxzpwXHC+wtb+nLLl78uT0x52pn78AM7vvSx+CRC9LeJ5eFJEn3hya7jFg1zsGLFMSOlmtB7FIurCoVBYkiIiLxsqyPV2g/vGlY7zy/aE3X34d9YDHHH/Y4e24fBr4v5nx/FpIk25aw7l5jxecamDRuDZ2fbaDu3sTAuBrn4E2aOZPTWlsZNWECmDFqwgROa21l0syZ5e5aWSlIFBERqTBjxvRlF2/6YV9A+MS64xn+yT3YkCE89M2T4MEzgBwXkqQqmROWEDwmCYyrcQ7epJkzaVm3ji/t20fLunWDPkAEBYkiIiIV7XMXD8E9CAKjffDqh7Cp93L0oU9z3y+3Z7+QJFXJnBxU6xw8iaUgUUREZAA46KC+7OJVV/UdX7Px3fTsitvBpdALSe602NeysVU5B09iKUgUEREZYL75zSBYXLcu+fkRw3cw/4xZicFdoRbfFHluplQGBYkiIiLxctzqrlwmTAiCxXPOgaFDI0ednbv35+b/uoSFv/sU20NxY8EK8CRLChJFRETi5bjVXbm1tcHBBwcLSQ491Pj+NRvYvXc4LYsW8vZPv4rNdK5YOj/9TdIEwF0bJzF57jN0bZxU4J5LJVOQKCIiMsBFLyT51a/g8q8eylPzj+G/r3o/H3jnfwPwnY4rsJmOzXRevPkdwdBztGSBMdAdqmXatZ2s3XQUjdd10h2qjX+8VCkFiSIiIlUgfiGJGbz/nX/kga99hD9d/Z6YtuMufhGb6Vw/84tBsPjqoynv29zaxpZtY3Afyuat9cxZWAX1biQrChJFRESq3HsOfxRfYuy5fSgfPPL3vce/dOf1QXbxwBN5Y+EBCdnFtocvpfPJRkK7g+xhaPf+dDzRRNvDl5a0//0S3jWnKAt3BgkFiSIiItUoyRzDoUP28ftvn4M73H9/7Lm3tryBzXQWX3RBEExtXM68e35A98648jo765h3zw+K2fPCKNOuOdVEQaKIiEg1yrD45l//NVgZHQrBO9/Zd9msWxZjM53l37+V/2yaQ92I7TG3Hej7NEv2FCSKiIgMYiNGwP/8TxAw/vznfcfP+N5yrv7FVUw88HlGDNsBQM2IvQN+n2bJnoJEERERAeDss4Ngcfdu+MUvoOFdPXT942h27qkBnFEjXmbhqSOD4WjPbSs/GXgUJIqIiEiMYcPgzDPh1w8fyd//Ds2zYeiQvWzZNpZjr3yKazu+zMu3jAmCxU0ry91dKRIFiSIiIpLS4YfDojaju2cYd90Fh/7TaL6y9FoOuWQj5960hN//+Bp8icGdQ2DfnnJ3t88A2TWnkg0rdwdERESk8o0YATNmwIwZB/KXv8Att4xg8aIm7nrkXI4at5aLTvkJ53cfxFvr3oD33wGHzSxvh9PtjrNsbPJVzjX1FburTjkokygiIiI5OeoouOEG+Mfmt/DTn8KocYdx2e03cvDF/6C5dRF/uv3GcHbRYM+Ocnc3kcrjZEVBooiISBXr6oLJk4P3gogqUl37S2PWfsYfP1/L49edyidn1/LzP5/L+775J46/8nFu+W0Lby4Oz118+ZECdUBKRUFiHsysycxat27dWu6uiIiIJOjuhmnTYO1aaGwMfs5bimzbceN+y09+AptequHHPwYfNZmL2m5h3MWb+Ezbj3jq1ouCYPGPsyszuygJFCTmwd073L1l9OjR5e6KiIhIguZm2LIlqFazeTPMmVP8Z44aBRddBE88NZxHHoH/M+Mt3PZwC8d+9Sne/83/ZvFt+9hxx9tUQmcAUJAoIiJShdraoLMz2FEFgveOjuB4KZjB+94HP/0pbNo0lBtugDeGvZdZtyzm4MtfZ9OLlvkmUlYKEkVERKrQvHmJw8s9PcHxUnvrW+Gyy2Dt2iE88AB8+jM1jBtX+n70UnmcrKgEjoiISBWaPx8uvTQ2UCz3vstmMHVq8CorlbnJioJEERGRKtTcDCtXwvLlwVBzTQ2F2Xe5pj51jcFclbJeYRGe1dUF06dDezs0NOTZvwqkIFFERKRKtbXBpEmwYQPU18OiRQW4aSGDt1LWKyzwsyIrxzds2Evj1A10faeBupqevgZVUJhbcxJFRESqVF0drFgRBIqdncHPlahr4yQmz32Gro2Tyt2VrPWtHB/K5q31zFkYF4FXQWFuBYkiIiJVrKEB1qwpw3BoVNHt6FfXDSfHFPfuDtUy7dpO1m46isbrOukO1Za4o7lLWDm+e386nmiibdWssvar0BQkioiISOGFM2nRWcLuUC3T/nNxX3HvUC3NrW1s2TYmdUauAiVdOb6zjnntZVwVVASakygiIiJFEckSbnjtUBqv6+SEw1aHA8KguPeHr/kda//RQGh3kD2Mzsg1E1VHscLm9yVdOT6imwUzvlK+ThWBgkQREREpiugs4T9eG8fG185k774g9AiF4NHn3gPEFtWOZOSap97WdzDV/L5cViwXcFV2wsrx4TtoOq6D2ScvzvlelUxBooiIiBRc26pZdD7Z2Jsl3LNvvyStDNhH9Oy3ZBm5ro2TmD45SamZdCuW7yxuJrJv5fg+6kdvZtGFcXseVkFhbgWJIiIikp8kGb157S/RvXNk2stqR3TTMG4Nz7z43pQZub4h62AeY1dXP1ZpF2GlcWTl+PTpQ2hvn0hdQ3fmiwYYLVwRERGR/CQJwuZPv4K6Edtjjg0dsofhQ3cBfQHhA1/7CGPGBLuxJMvI9Q1ZB/MY58Ql7MqpbCvHS0RBooiIiBRc89TbaDy2k5qa4Oea4Ts484Rf8I4DNmH0DdHW1fT01XL8cmNMQer4IetQCDo6gqFeKT4FiSIiIlIUbS3NMVnCxRfNYsXcRiYdvDYmIOzNyB3xasz189oXJAxZ9/QEJWjSGYjFuSuR5iSKiIhIUUSyhNOnQ/sng6Cw4ZC1rLn26OQXxC0umR9KUmqmFhZEyhEmWbEcX3YnYbu8ZEq5h/QAokyiiIiIFE2qLGGvNKuAm5uDxSq9Q9Y10NQEs2eHG5z1Epzrfa+a+uTFuTOtNC7lHtIDiDKJIiIikp9sahD2MyPXV2oG6uthUZoNWdreeInONRDaHfwc2r0/HU/PoO2NGTT36+mDm4JEERERyU8Rh2T7Ss0EdRLTlb9Jul1eD7S0wHvfW72rkItFQaKIiIhUtMiQNcvGwlOp5w4m2y7PDPbty6PG4iCmOYkiIiJSEF1dMHly8F4UGeYOxs9hHDIkeFVijcWBQEGiiIiI5K27G6ZNg7Vrg0Atfti3VNraYMyY4L/dYe/e4L/T1lhMtbClCrbWy4eGm0VERCRvzc2wZUts1m7p0tL3IzKH8Zhj+gLEiEiNxeb4VSyDuMxNOsokioiISF7a2qCzM8jWQfl3RmlogNbWxPmHMTUWJSMFiSIiIpKXVKuKM+2MUkwZayxKRgoSRUREJC/z55coa5fj3MHI/ESzzDUWJZHmJIqIiEhempth5UpYvjwYai5a1i7HuYO51FiURMokioiISN4qNWvXuy2gCmnnTEGiiIiI5C2StZs0KVjEUk1Zu6LXf6xQGm4WERGRgujdGaWKROo/btgw+HZtUSZRREREJIVk9R8HCwWJIiIiIklUWv3HUlOQKCIiIhWldw7gDSfDnZb4Wja2JP2oxPqPpaQgMYqZ3WZmu8xse9TrxHL3S0REZLCI2QP6msV0h2oTG4U2l6QvJav/WKEUJCb6kbuPjHo9Wu4OiYiIDBYxcwC31jNnYflq6Qz2XVsUJIqIiEhFSJgDuHt/Op5oom3VrLL2qRLrP5ZCRQWJZjbPzO42s+fMzM1sXZq2Q8zs82b2VzMLmdkGM7vezPJdmP5JM3vNzLrM7ItmVlG/RiIiItUq6RzAnXXMay/f+G4113/MpNICoG8DHwGeBV7P0Pb7wPeAtcAlwN3ApUBHfGBnZiPMbGSa19Bw0xuBI4GDgE8Bl4Vfe82w4AAAEeFJREFUIiIiUmRJ5wCO6GbBjK+Up0Nhg3XXlkoLEg9397e7+78Cm1I1MrMGgsBwmbuf5e4L3f0LwBeADwMz4i5ZBLyZ5vVBAHd/3N1fdve97v4IsACYXtBPKCIiIkklzAEcvoOm4zqYffLi2IY19aXv3CBUUUGiuz+XZdNPAAbcEHd8IdADnBd33/Pc3dK8VqV4zr5c+i8iIiL5iZkDOG5/Fv1mBpzrsa+zXirIswbrdnvZqqggMQcnEgRwf44+6O4h4Mnw+ZyZ2TlmNsoCU4ArgP+Xb2dFREQkO6WaAxhTaqcxcS6kDNwgcRzwirvvTHLuH8CBZrZfP+57MfACwRD0EuBHwPX97qWIiIjkrBRzAAfzdnvZGlbuDvRTLZAsQAQIRbXZlctN3f1D2bY1sxagBaC+vp5Vq1bl8qh+2b59e0meIyLSX/qekoFgxYqxLF/+TkKhYN1qKAT33ruXuXP/l2nTgqHs55+v5eqrG/jGN7o47LCecna3bMzdy92HpMxsDTDS3ScmOfcMMMbdE2aumtnPgbOBEe6eU5DYX1OmTPHVq1cX/TmrVq1i6tSpRX+OiEh/6XtKKl1XFxxzDOzdm3huzJggq9jdHQx3b9gA48cH11RT6Rsze8zdp2RqN1CHmzcRDCmPSHLuYIKh6JIEiCIiIjIwROYh7t0bLIyJFr3dnoaiAwM1SHyUoO/viT5oZjXAsUDx03oiIiIyoESCPwiCxKHhKsnR2+0l7PoSgo6O4PhgM1CDxHbAgcvjjl9IMBdxScl7JCIiIhUrPvjbty94Qd92e11d0NKSZNeXnmA3mMGmohaumNn5wITwjwcB+5nZ18I/r3f32wHc/Rkzuxm42MyWASuAowh2XHkQuLO0PRcREZFKlmzLP/cgm9jZGfwcPRQdvWSjthYWnH0p3PnDxBvX1BesbmOlqaggEZgDnBx37Fvh9weB26OOXw6sI1hh3Ai8AvwQ+Ia7l6QItpk1AU1HHHFEKR4nIiIi/TR/Plx6aWygWFsLN90UlNqZPj12KHrIkCBg7B2K/kCSABEgtLn4nS+TihpudvepaXZFmRrXdq+7X+/uR7r7CHc/2N2/4O7bS9jfDndvGT16dKkeKSIiMuj1Z6eUhC3/0sxDTDYUPRhVVJAoIiIikk4+O6XEbPkXFfxlGoqupvI3uVCQKCIiIgNGPuVpUm35N39+YiBYWwsLFxZ315dKpyBRREREBoRClKdJtuVfuqHowUxBooiIiAwIyYaFC1WeJtVQdK+ahE3e0h+vApW2ullEREQkqVQrlCM7peQjMhQ9fTq0tyeZh1ilZW7SUSZRREREBoRiDwsnG4oezBQk5sHMmsysdevWreXuioiIyKCQcVhYCkZBYh5UJ1FERKS0kq1Q7k/dRMlMcxJFRERkQIkMC0Nf3cQNG4Kh6K6uwVvXsNCUSRQREZEBK5+6iZKegkQREREZkApRN1FSU5AoIiIiA1Ix6yaKgkQREREZoFJtp1eIuomiIFFEREQGKG2nV1wKEvOgOokiIiLlpbqJxaMgMQ+qkygiIlJeyeomSmGoTqKIiIgMaNF1E6VwFCSKiIiIFMKysRDanHi8ph7Oeqn0/cmThptFRERECiFZgJjueIVTkCgiIiIiCRQkioiIiEgCBYkiIiIikkBBooiIiIgkUJAoIiIiUgg19bkdr3AqgZMHM2sCmo444ohyd0VERETKbQCWuUlHmcQ8aMcVERERqVYKEkVEREQkgYJEEREREUmgIFFEREREEihIFBEREclSVxdMnhy8VzsFiSIiIiJZ6O6GadNg7VpobAx+rmYKEkVERESy0NwMW7aAO2zeDHPmlLtHxaUgUURERCSDtjbo7IRQKPg5FIKOjuB4tVKQKCIiIpLBvHmJw8s9PcHxaqUgUURERCSD+fOhri72WG0tLFhQnv6UgoJEERERkQyam4PFKjU1wc81NdDUBLNnl7dfxaQgMQ9m1mRmrVu3bi13V0RERKTI2tpgzBgwg/p6WLSo3D0qLgWJedDezSIiIoNHXR2sWAGTJgWLWOKHn6vNsHJ3QERERGSgaGiANWvK3YvSUCZRRERERBIoSBQRERGRBBpuFhEREcnFsrEQ2px4vKYeznqp9P0pEmUSRURERHKRLEBMd3yAUpAoIiIiIgkUJIqIiIhIAgWJIiIiIpJAQaKIiIiIJFCQKCIiIpKLmvrcjg9QKoEjIiIikosqKnOTjjKJIiIiIpJAQWIezKzJzFq3bt1a7q6IiIiIFJSCxDy4e4e7t4wePbrcXREREREpKAWJIiIiIpJAQaKIiIiIJFCQKCIiIiIJFCSKiIiISAIFiSIiIiKSQEGiiIiIiCRQkCgiIiIiCRQkioiIiEgCBYkiIiIikkBBooiIiIgkMHcvdx8GPDN7GVifwyWjgf5s+Hwg8Eo/rpPc9Pf3p1JV4ucpR5+K/cxC378Q98v3Hv25Xt9TpVGJf6/zVYmfqVq/qw5w94MyNVSQWAZm1uruLf24brW7TylGn6RPf39/KlUlfp5y9KnYzyz0/Qtxv3zv0Z/r9T1VGpX49zpflfiZBvt3lYaby6Oj3B2QtKrt96cSP085+lTsZxb6/oW4X773qMQ/OxKoxt+bSvxMg/q7SpnEAUT/QheRSqfvKZHqoUziwNJa7g6IiGSg7ymRKqFMooiIiIgkUCaxypjZ58zsz2YWMrNV5e6PiEg0MxthZgvN7Dkz225m/2tml5e7XyKSaFi5OyAF9yKwADgReH+Z+yIiEm8Y8BJwGvAc8G5gpZm96O7tZe2ZiMTQcHOVCv/L/Ex3n1ruvoiIpGNmbcB2d7+03H0RkT4abi4TM5tnZneHh1zczNalaTvEzD5vZn8NDyNvMLPrzayuhF0WkUGmFN9TZjYMOAl4usDdF5E8KUgsn28DHwGeBV7P0Pb7wPeAtcAlwN3ApUCHmen3UESKpRTfUzcS7C7xs7x7KyIFpTmJ5XO4uz8HYGZrgJHJGplZA8EX7jJ3/4+o488TfLnOAO4sfndFZBAq6veUmV1PkEX8iLvvKnz3RSQfykKVSeSLNwufAAy4Ie74QqAHOK+Q/RIRiSjm95SZ3UCweOUUd9dezyIVSJnEyncisA/4c/RBdw+Z2ZPh873C83siryFmVgPs07/SRaSIcv2eupFgGPvD7v5yyXopIjlRJrHyjQNecfedSc79AzjQzPaLOvY1YAdwHfDB8H/fX/ReishglvX3lJlNIBiaPgJ4PlwrcbuZ/ap03RWRbCiTWPlqgWRfvAChqDa7ANz9KuCqovdKRKRP1t9T7r6eYGhaRCqcMomVrwcYkeJcTVQbEZFy0feUSBVSkFj5NhEM1ST7Aj6YYIhH8w1FpJz0PSVShRQkVr5HCX6f3hN9MLwg5VhgdTk6JSISRd9TIlVIQWLlawccuDzu+IUEc3yWlLxHIiKx9D0lUoW0cKVMzOx8YEL4x4OA/czsa+Gf17v77QDu/oyZ3QxcbGbLgBXAUQQ7GTyICmmLSJHoe0pkcDN3L3cfBiUzWwWcnOL0g+4+NartUIJ/obcAE4FXCP7l/g13317UjorIoKXvKZHBTUGiiIiIiCTQnEQRERERSaAgUUREREQSKEgUERERkQQKEkVEREQkgYJEEREREUmgIFFEREREEihIFBEREZEEChJFREREJIGCRBGpWmY2y8zczKaWuy+Vysy+Y2bPm9nwqGOrzGxdiZ5vZva4mf20FM8TkewpSBSRimdmU8PBXuS118xeN7M1ZrbYzD5qZlbgZ15lZmcW8p6VxswOAy4Drnb33eXogwfbfl0FfNLMji1HH0QkOW3LJyIVL5wJfAC4C1gBGPAW4EjgTGA88BvgbHd/I+q6ocBwYJe778vxmQ4sdvdZBfgIFcnMbgH+HTg4OkgM79k80d0nlrAvzwKPu/vZpXqmiKQ3rNwdEBHJwePufkf0ATP7AnAt8AWCIPL0yDl33wvsLWkPBwgzGwWcC7SVK4sY5w5gnpmNdfeXyt0ZEdFws4gMcO6+192/CDwEfNTMToqcSzYn0cxqwkPJfzOzHjN7w8yeMbPrwucnhrOIABdED3NH3WO6mS03sxfMbKeZvWJmvzSzd8f3z8zWhef4vcvMOs3sTTPbamb3mNnYJO1Hmdk1ZvYXMwuZ2atm9pCZzYhr9w4z+3G4D7vMbJOZtZrZmCx/6aYBIwkysxmZ2dvN7JFw308JH4v8+p5iZt8ws/VmtsPM/mRm7wu3OTnc/24ze9HMvp7iEb8iyPpW9RC/yECiTKKIVItFwElAI0HAmMrNQDPwM+D7wFDgncBHwudfBs4Hbgf+ALQmucfFwGvhcy8BhwMtwMNmdry7/29c+4OBVcAvgC8DxwCfBkYBp0UamdkB4b43APcAPw737zjgY8DScLvxwCPAfuHP/SxwBPAZ4MNmNsXdt6b5NQA4Ofz+aIZ2mNlEYCXBEP/J7v5kXJMF4X7+INynLwIrzeyCcP9agSXAOcDVZvZ8fEYYeBzYCUwFfpKpTyJSfAoSRaRaPB1+/+cM7f4d+JW7X5DspLt3A3eY2e3Ac0mCGYCPhtv1MrOfAU8Cnwc+G9f+CGC6u/88qv0+4LNm9i53/2v48LcJAsRPu3tMcGpm0SM/PyTIuh3n7huj2twN/DHch6uSfb4ok4DX3f21dI3Ci0l+BWwFPuDu65I0Gwq8z913ha9ZC9xLEOi+390fDR9fBKwHPkcwvNzL3XeZ2UaCzy8iFUDDzSJSLbaF30dlaLcVaDCzyf19UCRADJdvGWVmBxJkIP8GvDfJJZuiA8Sw34XfjwjfawgwA/gLsDDJM/eF240myCouB0JmdmDkBawD/k5UdjKNgwiyoSmZ2anA78P3/ZcUASLAjyMBYtgfwu9/jASI4c+wC/gzQeY2mVeBbIfLRaTIFCSKSLWIBIfb0raCy4G3As+Y2bNmdquZnRGXqUvLzI4zs/uANwmCzpfDr6PD9473XJJjr4bf3x5+PzB87ZOevuzEkQTf3XOinhv9OhKoz+JjOMEq8VTqCeYrPgec4u6vpmkb8/nc/fXwfz6fpO3r9H3meBbul4hUAA03i0i1iCwa+Vu6Ru5+b3iO3TSCeXmnEgRcfzCzU+MyYgnC8wF/TxCMfiv8vG6C4OYGgsUg8dKtsLa490xBUqTdHcDiFG12ZLgHBAHlMWnOv0YwT7ARmEmS7GaUVJ8v15Xlbwv3S0QqgIJEEakWc8LvnZkahufh3UEw99AIFl7MBc4A7s5w+b8TBIIfd/cHok+Y2dsJFl/0x8sEWbZMBaX/ThBI7ufuv+nnswDWACeb2YHu/kqS87uBs4B24BYzG+7uP8rjeWmZ2QjgUGBZsZ4hIrnRcLOIDGhmNtTMvkuwsnmFuz+coe0B0cfCQ7tPhH98W9Sp7XE/R0SyYzFDtWZ2IZBQ0iZb4TmHdwGTzGxO/PlwMEt42HcFcFakzEx8OzM7KItHrgq/J9wjqk+7CVYk3wPcbGaXZXHf/jqOYGX0g0V8hojkQJlEERlIjjez88L/Hb3jygTgfoLi0Om8BXjRzJYTBIZbgMMISse8DnREtf0jcKqZfQV4gSCeXEqw0rcHuN3Mbgpf9y8Ew9fPkt/36tcISvHcamanEZTDMYIAahhBaR7C/X0I+H14VfUTBP/o/yeCbOjPyLy6+dcEcyqnAfelauTue8zsEwSZxRvMbJi7X9+vT5deY/gZvyzCvUWkHxQkishA8onwax9Bpm8jQebpLnf/dRbX9xDMGzyFYC7iSOBFgpXC8919U1TbzxLUVLySILgEWOruz5rZ6QTlar5KkFl8mGB+403AxP5+OHd/3czeH77vWQRD228CawnK3kTabTCzE4CvEASF5wEhYANBoBu/kjrZs7ab2R3AdDO7PN1cTHffa2bnEwRx3zWz/dx9fn8/ZwozgXu124pI5dDezSIig1R4Ac9fgYvd/dYy9uMMgrmIU9z9iUztRaQ0FCSKiAxiZraAoD7jP2da2V3EPjwOPOXus8vxfBFJTkGiiIiIiCTQ6mYRERERSaAgUUREREQSKEgUERERkQQKEkVEREQkgYJEEREREUmgIFFEREREEihIFBEREZEEChJFREREJMH/BxhffSSqJq/CAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# gravity plots of selected categories\n",
    "plt.figure(figsize=(10,7))\n",
    "\n",
    "for col,marker,line,color in zip(gravity_res.keys(),['o','s','d'],['-','-','-'],['darkred','orange','blue']):\n",
    "    res = {k: gravity_res[col][k][1] / (gravity_res[col][k][0] + gravity_res[col][k][1]) for k in gravity_res[col]}\n",
    "    x = np.array([k for k in res])\n",
    "    y = np.array([res[k] for k in res])\n",
    "    def f(x,a,b):\n",
    "        return a*x+b\n",
    "    m = (x>30)&(x<300)\n",
    "    p,pcov = curve_fit(f,np.log(x)[m],np.log(y)[m])\n",
    "    plt.plot(x,y,marker,color=color)\n",
    "    if col==\"month\":\n",
    "        label=\"data\"\n",
    "    else:\n",
    "        label=col\n",
    "    plt.plot(x[m],np.exp(f(np.log(x)[m],*p)),line,color=color,label=label+r\", $\\beta=%.2f$\" % p[0])\n",
    "plt.loglog()\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.xlabel(\"Distance (km)\")\n",
    "plt.ylabel(r\"$p_{edge}$\")\n",
    "plt.savefig(\"gravity_plots.png\",density=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data export"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = []\n",
    "with open(\"gravity_res_big_formatted.json\") as f:\n",
    "    for line in f:\n",
    "        d.append(json.loads(line.strip(\"\\n\")))\n",
    "\n",
    "df = pd.DataFrame.from_dict(d)\n",
    "\n",
    "def data_series(d,t):\n",
    "    temp = []\n",
    "    for k in d:\n",
    "        temp.append((float(k),d[k]['1']/d[k]['0']))\n",
    "    return pd.DataFrame(temp).rename({1:t},axis=1).set_index(0)\n",
    "\n",
    "for key in [\"month\",\"ABM1\",\"ABM2\"]:\n",
    "    pd.concat(df[df.data_source == key]\\\n",
    "        .apply(lambda row: data_series(row[\"gravity_dict\"],row[\"category\"]),axis=1)\\\n",
    "        .to_list(),axis=1)\\\n",
    "        .fillna(0)\\\n",
    "        .reset_index()\\\n",
    "        .rename({0:\"d_km\"},axis=1)\\\n",
    "        .to_csv(\"gravity_plot_\"+key+\".csv\")"
   ]
  }
 ],
 "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}