OSSEN23-MF-tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd \n",
    "import glob\n",
    "\n",
    "import networkx as nx\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# DATA FOLDER\n",
    "DATA_FOLDER = \"./data/\" # Change to \"/content/\" if working with Colab"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The tutorial will focus on the paper _Abhijit Banerjee et al. ,The Diffusion of Microfinance. Science 341,1236498 (2013)_. In the paper, the authors explore the diffusion of information on micro-credit in several villages in India. In this tutorial, we will work with the original data from the paper, test some of the authors assumptions, and reproduce some of their results. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "    1. Load data on villages\n",
    "\"\"\"\n",
    "\n",
    "# Household info\n",
    "info = ...\n",
    "# Leaders-adopters\n",
    "leadapt = ...\n",
    "\n",
    "# Load network data\n",
    "# Adjacency_matrix\n",
    "adjmats = dict()\n",
    " \n",
    "for file in sorted(glob.glob(DATA_FOLDER + \"adjmat*.csv\")):\n",
    "    # find key adjmat 1 -> 0\n",
    "    key = int(file.split(\"/\")[-1].split(\".\")[0].strip(\"adjmat\")) -1\n",
    "    # load data as numpy array\n",
    "    village_adj = np.loadtxt(file, delimiter=\",\")\n",
    "    adjmats[key] = village_adj\n",
    "\n",
    "graphs = ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "    1.1 Has everything worked? Plot a graph as an example.\n",
    "\"\"\"\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 10), dpi=200)\n",
    "\n",
    "# Compute nodes positions using spring layout\n",
    "pos = ...\n",
    "\n",
    "...\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "    2. Are the villages similar? \n",
    "    Compute and plot a few basic statistics of our networks: number of nodes, number of edges, average degree, density, and number of components.\n",
    "\"\"\"\n",
    "\n",
    "n_nodes = np.zeros(len(graphs))\n",
    "n_edges = np.zeros(len(graphs))\n",
    "avg_degree = np.zeros(len(graphs))\n",
    "density = np.zeros(len(graphs))\n",
    "diameter = np.zeros(len(graphs))\n",
    "\n",
    "...\n",
    "\n",
    "# Plot distribution of number of nodes\n",
    "...\n",
    "\n",
    "# Plot distribution of density\n",
    "..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "    3. BSS, the microfinance institution mentioned in the paper, targeted \"leaders\" in the villages, assuming that they have a central role in the social network.\n",
    "    Let's see if this is true.\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "    Plot results\n",
    "\"\"\"\n",
    "# Make scatterplot of leaders vs non-leaders centrality\n",
    "fig, ax = plt.subplots(figsize=(10, 10))\n",
    "\n",
    "# Scatter plot of degree centrality; blue circles\n",
    "ax.scatter(\n",
    "    avg_dc_leaders,\n",
    "    avg_dc_nonleaders,\n",
    "    color=niceblue,\n",
    "    s=100,\n",
    "    marker=\"o\",\n",
    "    facecolor=\"white\",\n",
    "    linewidths=2,\n",
    "    label=\"Degree centrality\",\n",
    ")\n",
    "# Scatter plot of betweenness centrality; red squares\n",
    "ax.scatter(\n",
    "    avg_bc_leaders,\n",
    "    avg_bc_nonleaders,\n",
    "    color=nicered,\n",
    "    s=100,\n",
    "    marker=\"s\",\n",
    "    facecolor=\"white\",\n",
    "    linewidths=2,\n",
    "    label=\"Betweenness centrality\",\n",
    ")\n",
    "# Scatter plot of closeness centrality; green triangles\n",
    "ax.scatter(\n",
    "    avg_cc_leaders,\n",
    "    avg_cc_nonleaders,\n",
    "    color=nicegreen,\n",
    "    s=100,\n",
    "    marker=\"^\",\n",
    "    facecolor=\"white\",\n",
    "    linewidths=2,\n",
    "    label=\"Closeness centrality\",\n",
    ")\n",
    "\n",
    "# Plot diagonal line\n",
    "ax.plot([0, 1], [0, 1], transform=ax.transAxes, ls=\"--\", c=\"red\")\n",
    "\n",
    "# Add legend\n",
    "ax.legend(loc=\"upper left\", fontsize=20)\n",
    "\n",
    "# Add labels\n",
    "ax.set_xlabel(\"Average centrality of leaders\", fontsize=20)\n",
    "ax.set_ylabel(\"Average centrality of non-leaders\", fontsize=20)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "    4. Assortative/disassortative mixing: do adopters and non-adopters connect more between themselves?\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now let's compute some random benchmarks!\n",
    "..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 162,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABScAAALpCAYAAABBiZOaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAACTsklEQVR4nOz9f3zedX0v/j/SHyRNm5bSInNYToC1JlPGmuBQ89lHj+O2CXibRY+fKc3ZFFqdHwU/YqersAEHHBwPh6G73fYTmHMt+vk4YepE9/ni2PkckXkkQWa3xAYl0tvx0CEtNDS50lKu7x81IWmupLny60py3e+3W2/J9f5xvV+B9sp1Pd6v1/NZUywWiwEAAAAAmGNLKj0AAAAAAKA6CScBAAAAgIoQTgIAAAAAFSGcBAAAAAAqQjgJAAAAAFSEcBIAAAAAqAjhJAAAAABQEcLJExSLxRw6dCjFYrHSQwEAAACARU04eYK+vr6sWbMmfX19lR4KAAAAACxqwkkAAAAAoCKEkwAAAABARQgnAQAAAICKEE4CAAAAABUhnAQAAAAAKkI4CQAAAABUhHASAAAAAKgI4SQAAAAAUBHCSQAAAACgIoSTAAAAAEBFCCcBAAAAgIoQTgIAAAAAFSGcBAAAAAAqQjgJAAAAAFSEcBIAAAAAqAjhJAAAAABQEcJJAAAAAKAihJMAAAAAQEUIJwEAAACAihBOAgAAAAAVIZwEAAAAACpCOAkAAAAAVIRwEgAAAACoCOEkAAAAAFARwkkAAAAAoCKEkwAAAABARQgnAQAAAICKEE4CAAAAABUhnAQAAAAAKkI4CQAAAABUxLJKDwAAmFv9/f3p7u4efjwwMJDe3t40NjZmxYoVw9ubmppSX19fiSECAABVQjgJAFWmu7s7ra2tJz2uo6MjLS0tczAiAACgWgknAaDKNDU1paOjY/hxV1dX2tvbs2vXrjQ3N486DgAAYDYJJwGgytTX15ecEdnc3GymJAAAMKc0xAEAAAAAKkI4CQAAAABUhHASAAAAAKgI4SQAAAAAUBHCSQAAAACgIoSTAAAAAEBFCCcBAAAAgIoQTgIAAAAAFSGcBAAAAAAqQjgJAAAAAFSEcBIAAAAAqAjhJAAAAABQEcJJAAAAAKAihJMAAAAAQEUIJwEAAACAihBOAgAAAAAVIZwEAAAAACpCOAkAAAAAVIRwEgAAAACoCOEkAAAAAFAR0won77///lx00UU57bTTsnLlyrS0tOSP/uiP8uKLL5b1PDfccENqamom/NPd3T3u+V1dXdm6dWte/vKXp66uLueee2527NiRZ599djo/HgAAAAAwi5ZN9cRbb701O3fuTJKcc845WbVqVR577LFcffXVeeCBB3LfffdlyZLyss8NGzbkrLPOKrmvvr6+5PYHH3wwl156aQYGBnL66afnVa96Vbq7u/Nf/+t/zX333ZdvfetbOeOMM8r74QAAAACAWTelmZMPP/xwPv7xj2fJkiW555578oMf/CCPPfZYOjs7c8YZZ+TLX/5ybr/99rKf94orrsg3v/nNkn9KhZZ9fX35jd/4jQwMDOTqq6/O//yf/zMdHR158skn09bWlh/+8Ie58sorp/IjAgAAAACzbErh5M0335xisZht27blXe961/D2888/fziUvPXWW3P06NGZGeU4/vRP/zRPP/10mpubc/vtt2f58uVJknXr1uWee+7JsmXL8tWvfjWdnZ2zOg4AAAAAoHxlh5OHDh3KAw88kCQlZyW+4x3vyOrVq/PMM8/kwQcfnP4IJ3DvvfcmSd797ndn6dKlo/adddZZueiii5Ikf/M3fzOr4wAAAAAAyld2OPnoo4/myJEjqaurS0tLy5j9y5cvz2te85okybe//e2ynvvBBx/MO97xjrzpTW/Kf/gP/yGf/OQn89RTT5U89oUXXkhHR0eSpK2treQxQ9vLHQcAAAAAMPvKbojT09OT5PjMxGXLSp9+zjnn5Bvf+MbwsZP1//1//9+ox1/84hdzww035I//+I/z7ne/e9S+3t7e4WXj55xzzrjjGDlmAAAAAGD+KHvm5MGDB5Mka9euHfeYoX1Dx57My1/+8nz84x/Pd77znTzzzDPp7+/PQw89lIsvvjgDAwO54oor8pWvfKXkOCYay2TGMTg4mEOHDo36AwAAAADMvrLDyUKhkCQ55ZRTxj2mtrY2STIwMDCp53zf+96XT3ziE7ngggty2mmnZcWKFXn961+fr371q7nssstSLBbz4Q9/OMViccw4JhrLZMZxyy23ZM2aNcN/NmzYMKkxAwAAAADTU3Y4WVdXlyQ5cuTIuMcMDg4mSVasWDHFYR1XU1OTW2+9NUnygx/8IP/8z/88ZhwTjWUy49i5c2eee+654T/79u2b1pgBAAAAgMkpO5yczFLpySz9nqxNmzbltNNOS5I8/vjjY8Yx0VgmM47a2tqsXr161B8AAAAAYPaVHU5u3LgxSfLkk0/mhRdeKHnMD3/4w1HHTtfy5cuTZNT1Ghsbh7cPXW+2xwEAAAAAzJyyw8nNmzdn+fLlKRQK6ezsHLP/6NGj+c53vpMkufDCC6c9wJ/85Cf5t3/7tyTJK17xiuHty5YtS0tLS5LkoYceKnnu0PaZGAcAAAAAMLPKDidXr16diy66KEly1113jdn/hS98IYcOHcq6devyxje+cdoDvP3221MsFrNmzZq85jWvGbXvbW97W5LkM5/5TI4dOzZq35NPPpkHHnggSfL2t7992uMAAAAAAGZW2eFkklx77bWpqanJnXfemc997nPD2x977LFcc801SZKPfvSjo7po33HHHWlsbMw73/nOUc/1L//yL/k//8//M//yL/8yanuhUMgf/MEf5D//5/+cJPnYxz42piv3b//2b2f9+vXp6urKNddck6NHjyZJnnnmmVx++eV54YUXcvHFF6e1tXUqPyYAAAAAMIumFE62tbXlpptuyosvvpjLL7885557bs4///y0tLRk//79ufTSS/ORj3xk1DnPPvtsfvSjH+Wpp54atf3o0aP5kz/5k7z61a/Oy172slxwwQW54IILsm7dulx77bV58cUXc+WVV+Z3f/d3x4xj9erV+fznP5+6urp8+tOfzplnnpkLLrggZ511Vh566KE0Njbm7rvvnsqPCAAAAADMsimFk8nx2ZNf+cpX8qY3vSnPPPNMHn/88Zx33nm544478qUvfSlLly6d1PM0NjbmpptuysUXX5xVq1bl+9//fr73ve/ltNNOy3/4D/8hX//613PnnXempqam5Pm/8iu/kkceeSTvfOc7U1NTk+9973s544wzcs0116SzszM/8zM/M9UfEQAAAACYRTXFYrFY6UHMJ4cOHcqaNWvy3HPPZfXq1ZUeDgDMus7OzrS2tqajo2O42RwAAMBcmPLMSQAAAACA6RBOAgAAAAAVIZwEAAAAACpCOAkAAAAAVIRwEgCqWLFQyPM7dhz/fnCwwqMBAACqjXASAAAAAKgI4SQAAAAAUBHCSQAAAACgIoSTAAAAAEBFLKv0AACAuVEsFEpse6kJTvHIkZLH1NTVzeq4AACA6iWcBIAq0Xf51gn39994Y/pOXTtm++p7vzhbQwIAAKqcZd0AAAAAQEWYOQkAVaLhnt1jthULg8nbLkuS1F9/fRouvHCuhwUAAFQx4SQAVInH9+1LX1/fqG3FwcHs/em27z/xRJY0NIw5r6GhIRs3bpyTMQIAANWlplgsFis9iPnk0KFDWbNmTZ577rmsXr260sMBgBnR09OTTZs2Tfn8vXv3CigBAIAZZ+YkAFSBoRmTu3btSnNz8/D24uBgnvnd382T/f1p/uQns2LNmlHndXV1pb29fcyMSwAAgJkgnASAKtLc3JyWlpbhx8VCIX3r1ue165KG170uNXV1FRwdAABQbXTrBgAAAAAqwsxJAKhiNXV1WX3vFys9DAAAoEqZOQkAAAAAVIRwEgAAAACoCOEkAAAAAFARwkkAAAAAoCKEkwAAAABARQgnAQAAAICKEE4CAAAAABUhnAQAAAAAKkI4CQAAAABUhHASAAAAAKgI4SQAAAAAUBHCSQAAAACgIoSTAAAAAEBFCCcBAAAAgIoQTgIAAAAAFbGs0gMAABa//v7+dHd3Dz8eGBhIb29vGhsbs2LFiuHtTU1Nqa+vr8QQAQCAChBOAgCzrru7O62trSc9rqOjIy0tLXMwIgAAFoNioZC+y7cmSRru2Z2auroKj4hyCScBgFnX1NSUjo6O4cddXV1pb2/Prl270tzcPOo4gMVg5Ixxs8UBYHzCSQBg1tXX15ecEdnc3GymJLAoTWbGuNniACCcBAAAmHEjZ4ybLQ4A4xNOAgAAzLBSM8bNFgeAsYSTAAAAAMx7xUKhxLbBkt+PpEnO/CacBAAAAGDeG+rKPZ7nr7ii5PbV935xNobDDBFOAgDMkJHdeZPSHXp15wUAgJcIJwEAZojuvAAAs6fhnt1jthULg8MzJlfdfXdq6mrnelhMk3ASAGCGjOzOm5Tu0Ks7LwDA1JysdmRNXa36kguQcBIAYIaU6s6b6NALAADjWVLpAQAAAAAA1Uk4CQAAAABUhHASAAAAAKgINScBAAAAWJBq6uqy+t4vVnoYTIOZkwAAALOkWCjk+R07jn8/OFjh0QDA/COcBAAAAAAqQjgJAAAAAFSEcBIAAAAAqAgNcQBghP7+/nR3dydJBgYG0tvbm8bGxqxYsWL4mKamptTX11dqiAAAAIuGcBIARuju7k5ra+uEx3R0dKSlpWWORjQzBo4eS9ZuSNf+w8kTT0/6vK79h5O1G46fD8CEioVCiW0vNcEpHjlS8piaurpZHRfAYmRSweIhnASAEZqamtLR0ZEk6erqSnt7e3bt2pXm5uZRxyw0vQcKySU70757T5I95Z18yc70HiikbVZGBrB49F2+dcL9/TfemL5T147ZvvreL87WkAAWrcU6qaAaCScBYIT6+voxb2Cam5sX/JuaM+qKyf235Lrf+70x4ergYCE//vH/ys/+7MtTWzt69s4TTzyR37vuujRe5YMzAADzx2KdVFCNhJMAUAV6f7A3ObgvN1/z3imd/7LTTp3ZAQEsQg337B6zrVgYTN52WZKk/vrr03DhhXM9LIBFabFOKqhGwkkAqAJbtmxJcvzucU1NTXp7e4f3PfHEE/m93/u93HTTTTn77LOHtw/V7GloaMjGjRvneMQAC8/JakfWnHKK+pIAcALhJABUgfXr12fbtm1Jks7OzrS3t4855vd+7/dGPVajBwAAmG3CSQCoMiPr8yQTdzcEAACYTcJJAKgyperztLXpxQ0AAMy9JZUeAAAAAABQnYSTAAAAs6Smri6rbrvt+Pe1tRUeDQDMP8JJAAAAAKAihJMAAAAAQEUIJwEAAACAihBOAgAAAAAVIZwEAAAAACpiWaUHAADj6e/vT3d39/DjgYGB9Pb2prGxMStWrEiSNDU1pb6+vlJDBAAAYBqEkwDMW93d3WltbZ3wmI6OjrS0tMzRiAAAAJhJwkkA5q2mpqZ0dHQMP+7q6kp7e3t27dqV5ubm4WOYv3p6etLX1zdme1dX16ivJ2poaMjGjRtndWyzrVgo5PkdO45/PzhY4dEAAMD8JJwEYN6qr68vOSuyubnZbMkFoKenJ5s2bZrwmPb29nH37d27d8EHlAAAwMSEkwDArBiaMTlypuuQUvVDhwzNkC014xIAAFhchJMAwKwab6ZrW1tbBUYDAADMJ8JJAIBpGK+uZnFwMHt/ur37+99PTW3tmGMWQ21N4CXVXGcXAKZKOAkAMEWTqauZJO3vec+4+9TWhMVBnV0AmBrhJADAFA3NkPqLCy7IpoaGMfsLx47lyf7+nFVfn7qlS0ft29vXl+2PPKK2JiwS6uwCwNQIJwEApmlTQ0N+8dS1Jfe9dt0cDwaoKHV2AaA8wkkAgGla+YlPpGHz5lHbioXBPH/FFUmSVXffnZq60TUnVz76aPL618/ZGAEAYD4STgIATFNNbW1q6urG3183dn+pBjkAAFBthJMAC1B/f3+6u7uHH5eqZdXU1JT6+vpKDREAAABOSjgJsAB1d3entbV1wmM6OjpK1rwCAACA+UI4CbAANTU1paOjY/jxUKfPkR1Cm5qaKjU8AAAAmBThJMACVF9fX3JW5HgdQgEAAGA+Ek4CTNNk6j8makBCtampq8vqe79Y6WEAAMC8JpwEmKbJ1H9M1IAEAACYjp6envT19Y3Z3tXVNerriRoaGrJx48ZZHRtTJ5wEmKbJ1H8cOg6qycDRY8naDenafzh54ulJn9e1/3CydsPx8wEAIMeDyU2bNk14THt7+7j79u7dK6Ccp4STANOk/iOU1nugkFyyM+279yTZU97Jl+xM74FC2mZlZAAALDRDMyZPnASSjF9aK3lp8kipGZfMD8JJ5h31+wAWh8bT6pL7b8mu3bvTXMbM4a7u7rRv3ZrGq9RrBBYOs8UB5sZ4k0Da2tzWXqiEk8w76vcBLA4rli9NDu5L8xkr03L26ZM/8eC+5OC+4+cDLBBmiwPA1AgnmXfU7wOA8hULhfRdvjVJ0nDP7tTU1VV4RFBdzBYHgKkRTjLvqN8HAMBCY7Y4AEzNkkoPAAAAAACoTsJJAAAAAKAiLOsGAJiiqXbnTXToBQCARDgJADBl0+rOm0y5Q2+xUCixbbDk9yNpkgMAwHwjnAQAmKKpdudNptehd6gr93iev+KKkttX36sbMAAA84twEgBgiqbcnTfRoRcAoAxTLaejlM78J5wEAGZFf39/kqSzs3PMvoGBgfT29qaxsTErVqwYta+rq2tOxreQNdyze8y2YmFweMbkqrvvTk1d7VwPCwBg1kyrnM4US+kwN4STAMCs6O7uTpJs3759Suc3NDTM5HAWlZPVjqypq1VfEgBYVKZaTmc6pXSYG8JJAGBWbNmyJUnS1NSU+vr6Ufu6urrS3t6eXbt2pbm5ecy5DQ0N2bhx41wMEwCABWDK5XSU0pn3hJMAwKxYv359tm3bNuExzc3NaWlpmaMRAQAA882SSg8AAAAAAKhOZk4CMCn9/f3DNQST0g1NSi3fBYCRRv4+Ga851kL8faIJGABMjXASgEnp7u5Oa2vrhMd0dHRYostJFQuFPL9jx/HvBwcrPJrFo6auLqvvVeid+W+x/j7RBAwApkY4CcCkNDU1paOjY/hxqYYmTWV0zQOgOo38fTJec6yF+PtEEzAAmBrhJACTUl9fX3IWi4YmVLOJlnEmlnJCKaV+nyyG3yWagAHA1AgnmdcW89I/9fsAFr7pLuNMLOUEAFgIRn6GLw4O5pmP/W6eHOhP8yc/mRVr1gwf53N8+YSTUCGLtd4SQDWZaBlnYiknAMBMqXTjsXE/w7/pTaMe+hxfPuEkVIj6fUA1KBYKJba9NBO+eORIyWNq6upmdVwzZTLLOBNLOQEApqvSjcdGfoYvDg7m0auvzvZHHsmuv/zLNP/CL4w6jvIIJ6FC1O8DqkHf5Vsn3N9/443pO3XtmO26TgMAMFKlG4+N/AxfLBRy+KdhZ9MrX+kz/DQJJwGoej09Penr6xuzfWgJyERLQSzLBQCA2afx2OIlnGTeWOxL/4D5qaenJ5s2bZrwmPb29gn37927V0A5joZ7do/ZViwMJm+7LElSf/31abjwwrke1qIysjj7ePWWFGZfeIqFwvDM44Z7dnu/AwAVIquYfcJJ5g1L/4BKGJoxWWoJyESFtZOXlo+UmnXJcSd7U1ZzyineuE2TBmsAALNHVjH7hJMAkPGXgLS1tVVgNDB5I4uzj1dvSWF2AADmK+Ek84alfwBQvlIN1tRbAljYlOyA+UNWMfuEk8wblv4BAAAo2QHziaxi9gknAQCggk5aaH/E9yP5IASLl5IdQDURTgIAc6qmri6rbrstaW1NTW1tpYcDFXeyQvvPX3FFye0K7cPipWQHUE2Ek5yUeidMh78/AAAAwHimFU7ef//9uf3229PZ2ZnBwcG88pWvzHve85584AMfyJIlS6Y1sDvvvDPbt29Pklx55ZW58847xxzzmc98Ju95z3smfJ6vfe1refOb3zytsVQ79U6YDn9/AGBi4xXaH5oxueruu1NTZ5YxwEwqFgrDM9cb7tmtVAZU0JTDyVtvvTU7d+5MkpxzzjlZtWpVHnvssVx99dV54IEHct999005oHz66afzsY99bNLHv+xlL8vGjRtL7lu7du2UxsBLKlnvxNK/hU+9HACY2EkL7dfV+tAMAPOIrGJmTSmcfPjhh/Pxj388S5Ysya5du/Kud70rSfLYY4/l137t1/LlL385t99+e3bs2DGlQX34wx/Os88+m0svvTRf/epXT3r8xRdfnM985jNTutZCU4klsuqdMB3+/gAAsBiYaQcwO6YUTt58880pFovZvn37cDCZJOeff35uv/32bN26Nbfeems+9KEPZfny5WU99wMPPJDdu3fn/e9/f172spdNKpysJpbIAtWqWCjk+Z/e9CoOlu5cCwAAwMJSdjh56NChPPDAA0mO14I80Tve8Y68//3vzzPPPJMHH3wwv/qrvzrp5y4UCsOh5B/8wR/kjjvuKHd4i54lsrNj5IzUpPSsVE1bAAAAAGZW2eHko48+miNHjqSurq7k7Lzly5fnNa95Tb7xjW/k29/+dlnh5M0335zHH388f/VXf5VTTz110uc99thjufzyy/PUU09l9erV2bx5c9rb23PuuedO+jkWCktkZ4cZqQCz68SbQF1dXaO+DnEjCACYacVCocS2wZLfj2TpPsyNssPJnp6eJMlZZ52VZctKn37OOefkG9/4xvCxk9HV1ZX/8l/+S375l385v/mbv1nWmL773e/mu9/97vDjL33pS7npppty44035tprr53w3MHBwQyOWB546NChsq7N4jByRmpSelaqGakAUzfeTaD29vZRj90IguNq6uqy+t4vVnoYAIvCUK3Q8Tx/xRUlt3sdhrlRdjh58ODBJBN3wR7aN3TsyRSLxbzvfe/Liy++mD/+4z+e9FhOPfXUXHXVVXnnO9+Zn/u5n8uaNWvS1dWV22+/PX/913+d6667LmvWrMkHP/jBcZ/jlltuyY033jjpa7I4lZqRmpiVCjBTTrwJNFFTNwCoNDPtAOZO2eFk4acv0qeccsq4x9T+tI36wMDApJ7zrrvuyn//7/89O3bsyKtf/epJj2XLli3ZsmXLqG2/+Iu/mM9+9rNZt25d7rjjjlx33XX5rd/6rTQ0NJR8jp07d+aaa64Zfnzo0KFs2LBh0mMAYGb19PTk0E9+MmZ78ciR7O3rS5J07dlT8tzV69dn48aNMzIOHTlnVqmbQG1tbRUaDQBMzEy7xaXhnt1Jflpm5vvfT3L8veWB3/v9PDnQn6abb86K1auTJE2vfKUSMzDHyg4n63764ezIkSPjHjO0THrkTIjxPP300/nYxz6WV7ziFbn++uvLHc64brzxxvzJn/xJnnvuufzDP/xD3vrWt5Y8rra2djhMBaCyenp6smnTppMe9x/f+95x9+3du3fGAsq5MrIe40QzCr1RBlg4Rr62q7MLlTV0k/n7//qvueD1rx97wJvfPPytEjMw98oOJyezZHsyS7+HfPSjH82BAwfyZ3/2Z1m1alW5wxnX6tWr86pXvSqdnZ15/PHHZ+x5mT09PT3p++msqJHGezM3pKGhYcEFEUBpQ68Bf3HBBdlUYsZ74dixPNnfn7Pq61O3dOmofXv7+rL9kUdKvo5MZODosWTthnTtP5w88fTw9uKRwRyuOT6GlT/6SWpOGXsjq2v/4WTthuPPMQ2acgHJ2MZRpW5WCLMWjlKv7ersLhxDM+1GKhYGh2dMrrr77tTUmeSy0IwsM1McHMyjV1+d7Y88kl1/+Zdp/oVfGD4G5otqeW9Qdjg5FAI9+eSTeeGFF0o2xfnhD3846tiJPProo0mSD37wg2NqQz7//PNJknvuuSd/93d/lyR56qmnJj3W5cuXJ0leeOGFSZ9DZUxmttSJb+ZGWogzpYDxbf70p9OyefOobUMfCF67rvQHgpWPPpqUuhN+Er0HCsklO9O+e0+SE5aLr7jw+NebvjL+E1yyM70HCpnOAuWRb5RLNeQaOgZY3NyoWFxGvrars7vwnKycS01drZIvC9DIMjPFQiGHf3ozvOmVr/TayrxULe8Nyg4nN2/enOXLl6dQKKSzszO/9Eu/NGr/0aNH853vfCdJcuGFF076effv3z/uvoGBgUnXrxxy7NixfP+ntSRe8YpXlHUuc29optOJH8aT8d/MJS99iC93phQwv9XUTvyGv9QHgpopluhoPK0uuf+W7Nq9O80jPiQWjwzm8MevS5Ks/IObS8+c7O5O+9atabxqevWlStVj1JBr+qrlTjOLx4mNo0rdrBBmLRwnvrarswtAuarlvUHZ4eTq1atz0UUX5Wtf+1ruuuuuMeHkF77whRw6dCjr1q3LG9/4xpM+33e/+91x991www258cYbc+WVV+bOO+8sa5x33XVXnn322SxdunRS42B+GO/DuDdzwGxZsXxpcnBfmk5dls0vf2kpebFwSp4vHr/xseqMhtJLtw6uTA7uO/4czDvVcqeZxaPUjYrEzQoAqFbV8t6g7HAySa699tp8/etfz5133pk3vvGNede73pUkeeyxx4Y7X3/0ox8d1dH7jjvuyB133JHXvva1+fznPz/tgR86dCjve9/78uEPf3hUQHrs2LHcfffd+dCHPpQkufLKK3PmmWdO+3oALG6Hr702faeWrpU8XkfOXHftLI6I6aqWO80AQPlq6uqy6rbbktbWKa/AAWbGlMLJtra23HTTTbnuuuty+eWX57rrrsuqVauyZ8+evPjii7n00kvzkY98ZNQ5zz77bH70ox+lsbFxJsadF198MZ///Ofz+c9/PqeeemrOPvvsLFu2LD09PXn22WeTJBdffHE+9alPzcj1AICFpVruNAMAwEI2pXAyOT578vzzz88f/uEfpqOjI0899VTOO++8vOc978kHP/jBLF06u0vcVq5cmU9+8pP51re+lT179uQHP/hBBgYGsm7dulx66aX5zd/8zbzjHe9ITU3NrI4DgMVh5Sc+kYYRTXgm1ZHzX/91roYHwDSMrEE7UXMaNWiZSE1dXVbfO7060wCMNeVwMkne8pa35C1vecukjr3hhhtyww03lPX8E52zfPny/M7v/E5ZzwfAwjXbHwgmasKjIyfAwqYGLQDMX9MKJwEAAOa7kTVoS9WfHToGAJh7wkkAAGBRK1WDVv1ZAJgfhJMAADNkZF275PgMrZFfE3XtAABgJOEkAMAMGa+uXXt7+/D36toBAEzfyJvCpW4IJ24KLxTCSQAoQUdOpmJkXbukdFdgde0AAKav1E3hkTeEEzeFFwrhJADADClV166tra1CowEAmBsjZzGWujmbzPwsxpE3hSe6JvOfcBIAAACAKRuvtM1IMz2L8cSbwm4IL1zCSQAAAJiHioVCnt+x4/j3g4MVHg2Mb+Qsxq6urrS3t2fXrl1pbm4edQyUIpwEAAAAYMpKlbZpbm5W75FJWVLpASwWxUIhh9729hx629tTLBQqPRwAAAAAmPfMnARg3hg4eixZuyFd+w8nTzxd1rld+w8nazccfw4AAAAWBOHkNIzsRlUcHMwzP/lJnhzoT/PDD2fFmjVJZr4bFcBi1nugkFyyM+279yTZU/4TXLIzvQcKUQobAABgYRBOTsO43aje9Kbhb2e6GxXAYtZ4Wl1y/y3ZtXt3msssmN3V3Z32rVvTeNUXZ2l0AACzZ++ePenr6xu1rXjkSPb+dFvXntI3bmtqa9PQ0JCNGzfO+hihGvX09Iz5t5kcb/wz8uuJ/LucPOHkNIzsRlUcHMyjV1+d7Y88kl1/+Zdp/oVfGD4GgMlZsXxpcnBfms9YmZazTy/v5IP7koP7jj9HGfr7+5MknZ2dY/YNDAykt7c3jY2NWbFixZj9470RAQAoR09PT1553nkTHvMf3/veCffv3btXEAIzrKenJ5s2bZrwmPb29nH3+Xc5OcLJKRhqeLNiyZJs/vmf/+m2wRxuaEiSvPLss4e3jzy+pq5ujkcKwMkMlefYvn37lJ+j4aev/9NVLBTy/I4dx78fHJyR5wQA5r+hWVl/ccEF2XTC+4rCsWN5sr8/Z9XXp27p2Juw+668Mu3t7SVndgHTM/TvateuXWlubh61b6KJDF1dXf5dlkE4OQV9l2+dcH//jTem79S1Y7avvtdSQ4D5ZsuWLUlK1wgeelNR6s3IEMs1AICZsvnTn07L5s3Dj4uFwTx/xRV57bpk1d13p6audsw5j/7rv87lEKEqNTc3lyzZ19am2v1MEE4CUNXWr1+fbdu2TXjMeG9GgMVrZOPDZPzZEZofAjOpprZ23BV3NXXj7wNYyISTU9Bwz+4x24qFweRtlyVJ6q+/Pg0XXjjXwwIAYIaM2/jwBJofAgBMj3ByCk52t6rmlFPc0QIAZs1Uu0YmShFM1sjGh8n4ZR40PwQAmB7hJABUwFCztNHbXmqCUzxypOQxbn4x3a6Ric6Rk1FfX19yRqQyDwAwPg0emQrhJPPCwNFjydoN6dp/OHni6Umf17X/cLJ2w/HzoUp5A7Awaa7GVE21a2SicyTAVIysQav+LMDME04yL/QeKCSX7Ez77j1J9pR38iU703ugED2yAKgmukYCzI3J1KCdrfqzNXV1bkwCi55wknmh8bS65P5bsmv37jSXUbupq7s77Vu3pvEqv7BhLpmtOX2aq8HC4TUPqtvIGrTqzwLMPOHkDKmpq8uq225LWltTU1tb6eEsOCuWL00O7kvzGSvTcvbpkz/x4L7k4L7j5wOzolTjjeLgYPb+dFv3979f8nVP042Jaa4GAAtDqRq06s8CzBzhJADjmlTjjfe8Z9x9mm4AAMDipMEjM0U4CbCA7N2zZ+wsxiNHhmcxdu0pXbN19fr1UwoJx2u8URwczDO/+7t5sr8/zZ/8ZFasWTPqPE03AABgcdPgkZkinAQWncXaUbGnpyevPO+8CY/5j+9977j7pjOL8cSlS8VCIX3r1ue165KG173O3U8AAACmRDgJLDqV7Kg4m4ZmIf7FBRdkU0PDqH2FY8fyZH9/zqqvT93S0TVY9/b1Zfsjj0xrFmNxcHDUkoxRyzUKY5tDaBgBLFSWqM2OUvWLk+Mz7Ud+LUUNY2C6pvoatNBff4qFwvDsxoZ7ds/47yoNHpkpwklg0VnsHRU3f/rTadm8efhxsTCY56+4Iq9dl6y6++7U1I1uTrPy0UeT179+Wtc8fO21JZdkJMnzV1wx9vhnD07retVKczWoPEvUZt6k6he3t0+4Xw1jYKqm+xrk9Wd8GjwyU4STwKKz2Dsq1tTWjvtLvqZu7D4hFwCVNF794mT88itD1DAGpmuqr0Fef2DuCCcBOKmVn/hEGkrM1kxmdrZmf39/kqSzs7Pk/pO9gQSYLkvUZs94Nwrb2toqMBqg2ngNgvlLOAnASc3VbM2hRkbbt2+f0vnJ8dpAAFNliRpQCQNHjyVrN6Rr/+HkiafLOrdr/+Fk7YbjzwF5qUFocXAwz/zkJ3lyoD/NDz+cFWvWDB+zEBuEsngJJwGYN7Zs2ZJk/DdL49UQHbLQi5YDANWp90AhuWRn2nfvSbKn/Ce4ZGd6DxRiDmD1Gtmsreuxx3LBiauY3vSmUQ8f+da30rJ5sxtuJzHVGwduGpRHOAnAvLF+/fps27btpMctphqiAACNp9Ul99+SXbt3p7nMxo1d3d1p37o1jVdpzFXNRjZ0O/OFF/Lf/v2/T5Ls7evL9kceyV9ccEE2jVhhdOYtt6Zv2bIZb+i22Bo8TuvGgZsGkyacBFjgaurqdIkFAFjAVixfmhzcl+YzVqbl7NPLO/ngvuTgvuPPAUnqly3LL566dtS2TQ0NY7ZxclO9ceCmQXmEkwAAAACzrFgoDM9wbLhn94wuqdbQbXZM+caBmwZlEU5SUk9PT/r6+sZsH+qGO15XXPXeoDqYrQlUk8W2RA2AxUdDNxYy4SRj9PT0ZNOmTRMe097ePu6+vXv3CigBYJbo6AoAwGIinGSMoRmTpbrhDgwMpLe3N42NjVmxYsWofUNddEvNuOS48WakJmalAjA5OroCALCYCCcZ13jdcNvaZv7jTH9/f5Kks7NzzL6TBaILxWRmpCZmpQIwMR1dAQBYTISTTNpsFu/t7u5Okmzfvn1K5zc0NMzYWGbLRDNSE7NSAZgcHV0BYP4rFgoltg2W/H4kdSGpRsJJ5oUtW7YkSZqamlJfXz9q31AwN16ot9CWO5eakToU/J6XmQ9+AQAAmFtDE3vG8/wVV5TcPpNNJzV0Y6EQTjIvrF+/Ptu2bZvwmPGWmQMAAACwMAknAaACxmuQpTkWAMDMGTh6LFm7IV37DydPPD3p87r2H07Wbjh+/hQ03LN7zLZiYXB4xuSqu+9OTZ3ZjJAIJwGYQKXezC12k2mQpTkWwMwrFgp5fseO498Plq73BiwuvQcKySU70757T5I95Z18yc70HihkKi1hT1aqq6auVjkv+CnhJOMqDg6OKuJ7suK93uDB4lOpN3OL3UQNsjTHAgCYOY2n1SX335Jdu3enualp0ud1dXenfevWNF41czUgWXj6+/uTJJ2dnWP2nex9O5MnnGRch6+9Nn2nri25r1Tx3sPPHpztIS0KJ4a+ieCX+cubudk1Xi3dtjaRLgDATFixfGlycF+az1iZlrNPn/yJB/clB/cdP5+q1d3dnSTZvn37lM5vaGiYyeEsWsJJmGMThb6J4Jf5xZs5AACgWm3ZsiVJ0tTUlPr6+lH7hlY1lVoNlagVXw7hJONa+YlPpGHz5uHHJyveu/LRR5PXv35Ox8j8NNVGH4kXcAAAAOaH9evXZ9u2bRMeM95qKCZPOMm4amrHL9BbqnhvTa1OY5NxYuibLK7gd7qNPhLNPgCA6TuxjM7xbSNK6Rw5UvIYDSqA2VJTV5fV9yp7BCcSTsIcmyj0TRZ+8DvVRh+JZh8AwMzpu3zrhPv7b7yxZKkdwQEAzC3hJGMMHD2WrN2Qrv2HkyeeHt5ePDKYwzXHi7mu/NFPUnPK6MCsa//hZO2G4+dT9TT6AGCm9Pf3DxekT8a/2VWqHhQAMDuU82KmCCcXqGKhkOd37Dj+/Qx3cu49UEgu2Zn23XuS7Bm9c8WFx7/e9JXSJ1+yM70HChE/wcwb78bBybhxACx03d3daW1tPelxHR0d0675NNUPWj5kzT8N9+wes61YGEzedlmSpP7669Nw4YVzPSwoqb+/P0nS2dk5Zt9kVh/BXFPOa3ZV2/sR4SRjNJ5Wl9x/S3bt3p3mpqbh7cUjgzn88euSJCv/4OaxMye7u9O+dWsar7IUBmbDhDcOTsaNA1g0qvEDbFNTUzo6OoYfj9cds2nE+5apmO4HLR+y5peT1Y6sOeUU9SWZN4Zmh2/fvn3Kz9HQ0DBTw4GTUs5r9lTj+xHhJGOsWL40ObgvzWesTMvZp4/e+cXPjH/iwX3JwX3Hz4d5YDZnGFfCeDcOTsaNA1hc5sMH2JHLrOdiiXV9fX3JGZEz3R1zqh+0fMgCpmvLli1JSr92jndDZqSFOluKhU85r5lXje9HhJMwD+jatjDN9VT7CW8cTMSNA1hU5sMH2Mkss56JJdaV4oMWMNfWr1+fbdu2TXjMTN+QAea3ano/IpwEmIJqnGoPzA/z4QPsyGXWs7XEGgCA6iCcXACKhUKJbS8tUS0eOVLyGDV0YPZU41R7gCGlllmb0QMAwFQIJxeA7/76W9P3wtEx2/f+NNz47o4d6S9RO2rVbbepPQKzrJqm2jNzdF4HqIyaurqsuu22pLU1NbW1Jz8BWPCm2khuoTaRg4VIODnP9fT0pOX/9/9OeMz2Rx4pveOntaAsH6UamGHMQqLzOlBNpnpDJpneTZm5rg0NzE/TbSSnCzrMPuHkPDf0huqv7757VHfe4pEjOfD7v58n+/vTdPPNWbF69Zhzu3/4Q8tHqRp9l2+dcH//jTem79S1Y7ZrREQl6LwOVJNp3ZBJpnRTRm1oYMh0Gsm5WQFzQzi5QPz8+eePWjpaLBTSt259XrsuaXjDG0rO/rJUBWB+0nkdpq5YKOT5HTuOfz84eJKjmQ+mekMmmfpNGbWhgSHzoZEcMDHhJLAoNNyze8y2YmEwedtlSZL6669Pw4UXzvWwFjw1egCYrinfkEmmfVNGbWgAFotioTC8YrDhnt2LqkSZcBJYFE72wlxzyimL6sV7rqjRAwAAwGwSTgIwLjV6AIBqNNWGSon3QADlEk4CMC41egCAajPdhkqJpkoA5RBOLlA1dXW6DC8wA0ePJWs3pGv/4eSJp8s6t2v/4WTthuPPAcyq/v7+4eXs482OKDWTFJh5xUKhxLaXmuAUjxwpeYwyHrB4jfw9nZSugT3d39NTbaiUaKpE9fD5lpkknIQ50nugkFyyM+279yTZU/4TXLIzvQcKUb4dZld3d3daW1tHbTtxdkRHR4fZojAHhoq+j6f/xhvTd+raMdvdwIXFq9Tv6RPN1O9pDZVgfD7fzr7i4OCom7CjbtCO+H7k8QuVcBLmSONpdcn9t2TX7t1pbmoq69yu7u60b92axqt82CpHTV1dVt12W9Lampra2koPhwWiqakpHR0dScafHdFU5r9hAGBmjPw9nZSuge33NMw+n29n3+Frry15EzZJnr/iirHHP3twtoc0a4STMEdWLF+aHNyX5jNWpuXs08s7+eC+5OC+488BzKr6+vpRsyTMjoDKabhn95htxcJg8rbLkiT111+fhgsvnOthwbxXLBSGZx433LN7UZU6OPH39BA1sGFu+XzLTBJOAgCLXrFQyPM7dhz/fgEveak2JwtUak45ZVGFLgAAQ1Z+4hNp2Lx5+HGxMDg8Y3LV3Xenpm706sCVjz6avP71czrGmSKcBAAAAIB5pKa2dtybsDV1Y/ct5FJmwkkAAAAAIP39/enu7h5+XKoOf1NTU+rr62fsmsJJAAAAACDd3d1pbW2d8JiOjo4ZrfMrnGSM/v7+JElnZ+eYfeN1rk2Od8oDAADmTrFQKLFtsOT3I6nXCkApTU1N6ejoGH7c1dWV9vb27Nq1K83NzcPHzCThJGMMTd/dvn37lM5vaGiYyeEAVezEJQVDN0FG3gyZ6SUFLHwn/aB+5EjJY3xQXxhq6uqy6rbbktbWBV1baaSenp709fWN2V7qNW+khoaGbNy4cVbHxvw31JW7/4UXsvf5l/4eFV44licH+nPWr/966pa91BF306qG1C9bltX3fnHOxwrA/FdfX19yVmRzc/OMzpYcSTjJGFu2bElS+gN/qcR8JG+SgZk03pKC9vb24e9nekkBC9/QB/Xx9N94Y/pOXTtmuw/qVEJPT082bdo04TEjX/NOtHfvXu+9SJLsfb4vb3jwwZMe99/+/b/PL5Z4DQRgfqupq1u071eFk2Ua7852snjubq9fvz7btm2b8JjZTMwBhpy4pGC8YswAM2Hg6LFk7YZ07T+cPPH0pM/r2n84Wbvh+PllGnpfWerG78nK6bS3t4/7vpTq0XDP7iRJa39/Hvn+95McnyH+3R07sv2RR/LXf/7naX71q4ePb3rlK604AGBeEU6WYTJ3thN3t6luU/1gl0zvwx2LU6klBW1tbRUaDQvFU//pxjGBTfHIkezdsSNJ8uS73pWVIz6oD1nd0+N3dJXrPVBILtmZ9t17kuwp7+RLdqb3QCFTfYUa78av17ypqUTQXClDJSlW1tWl9XWvS3K8vEX/T0stNb/61cPbAWA+Ek6WYaI724m72zNlMjXmEnXm5qtpfbBLpv3hDqhuPT09eeV55014zH9873vH3ecmYnVrPK0uuf+W7Nq9O81lzMru6u5O+9atabxqcS61WogqGTQDAOURTk7BREua3d2evsnUmEvUmZuvpvrBLpnehzvNBIBk/BuJxcHBPPO7v5sn+/vT/MlPZsWaNaPOcxORJFmxfGlycF+az1iZlrNPn/yJB/clB/cdP595QdAMwEJVTbP/hwgnmXcmU2Nu6Djmnyl/sEum/OGuWpoJ9Pf3J0k6OzvH7DvZzG2oNifeSCwWCulbtz6vXZc0vO51OnPDIidoBmChqsbZ/8JJ5h015ihXtTQTGCp3sH379imd3/DT2lMAAADMT9U4+184CSwai72ZwJYtW5KUrrc6FLSOVxPX8nUAKmWqM/8Ts/+nqqauLqtuuy1pbU1NbW2lhwMsQl7bZ081zv4XTgIsEOvXr8+2bdsmPGaimrhQzWrq6rL63oV3FxkWg+nO/E/M/qf6jGwSqkEo85HXdmaScBIAAJg105n5n5j9T3Uq1SRUg1Dmk2p8bR950yBx42AmCScBAIBZY+Y/lG9kk1ANQpmPqvG1vdRNg8SNg5kgnIQ5MlFNjkS3ZagmOq8DABM5sUnoYqmhDgvZyJsGiRsHM0k4Oc8NHD2WrN2Qrv2HkyeeLuvcrv2Hk7Ubjj8HFacmBzBE53UAmL98BgNKOfGmQeLGwUwRTs5zvQcKySU70757T5I95T/BJTvTe6AQ/1wqb6KaHIluy1BNdF4HgPnLZzCAuSWcnOcaT6tL7r8lu3bvTnOZU4O7urvTvnVrGq/SnXQ+mExNjmTx1eUAxqrGGj0AsFD4DAYwt4ST89yK5UuTg/vSfMbKtJx9enknH9yXHNx3/DkAAAA4KZ/BAObWkkoPAAAAAACoTmZOAkzBVAulK5IOAAAALxFOAkzBtAqlK5IOAAAASYSTwBwpFgrpu3xrkqThnt2pqaur8IimZ6qF0hVJBwAAgJcIJwGmYMqF0hVJB2ZYf39/uru7kyRdXV2jvg5pampKfX39nI8NAABORjgJALCAdXd3p7W1ddS29vb2UY87OjrS0tIyl8MCgHln5A29pPRNPTf0YO4JJwEAFrCmpqZ0dHQkSQYGBtLb25vGxsasWLFi1DEAMF/NVQmoUjf0ktE39dzQg7knnAQAWMDq6+tHfYhqa9NuayEZOHosWbshXfsPJ088PenzuvYfTtZuOH7+AlUsFPL8jh3Hvx8crPBogGow8oZeUvqmnht6MPeEk8CsKA4OplgovPS4MFjyewCoZr0HCsklO9O+e0+SPeWdfMnO9B4oRBw9Vn9/f5Kks7NzzL7xZhgnY+u1AovLiTf0Ejf1mH+q8XeYcBKYUUMvpN/6wAdy+NRTR+0rHDuWJ/v7c9Zbfz11S8c2hNl35ZVzMUQAmDcaT6tL7r8lu3bvTnMZs3W6urvTvnVrGq/64iyObuEaqim3ffv2KZ3f0NAwk8MBgEmrxt9hwklgRg29kF796KPln/zgg0kW5ovpkLmqlwPA4rBi+dLk4L40n7EyLWefPvkTD+5LDu47fj5jbNmyJUnpxhZdXV1pb2/Prl270tzcPObchoaGbNy4cS6GCQBjVOPvMOEkMKO2bNmS4tGjadq0adQLafHIkXx3x45sf+SR/PWf/3maX/3qMefW1NYu2BdTABa+alxGtVitX78+27Ztm/CY5uZmTS+gAkaWfnpp28lLQLnpT7Woxt9hwklgRq1fvz7b3//+MduLhUL6fzojsvnVr07r614310MDgAlV4zKquXLSMOLIkZLHCCNg8RlaZTSe56+4ouT21fcqYwGLlXASAABSncuo5srJwoj+G29M36lrx2wXRgDA4iecpKr19PSkr69vzPah5VnjLdPyAQQAFp9qXEYFMNca7tk9ZluxMDg8Y3LV3Xenpq52rocFVJBwkqrV09OTTZs2TXhMe3v7uPv27t0roCTFwcFRy9BOVi+nOFi6hg4ALGbjhRF522VJkvrrr0/DhRfO9bCACjhZuYaaulolHaDKCCepWkMzJkstzzpZ0fv29vaSMy6pPoevvbbkMrSkdL2cw88enO0hAcC8c9Iw4pRThBEAUKWEk2UYOHosWbshXfsPJ088Xda5XfsPJ2s3HH8O5pXxlme1tbVVYDSLV01dXVbddlvS2pqaWss0gNkx1d/Vfk8DAEBlCCfL0HugkFyyM+279yTZU/4TXLIzvQcKEXnB4rHyE59Iw+bNw49PVi9n5aOPJq9//ZyOEarJtH5X+z0NAABzTjhZhsbT6pL7b8mu3bvT3NRU1rld3d1p37o1jVfpOEjl9Pf3p7u7e/hxqeXrpTqUMr6a2vFr4pSql2PWKMyuqf6u9nsagCH9/f1Jks7OzjH7Jir/lIzfUBOA8Qkny7Bi+dLk4L40n7EyLWefXt7JB/clB/cdfw6okO7u7rS2tk54TEdHhy6kwII15d/Vfk8D8FNDN/O3b98+5edoaGiYqeFUhZq6uqy+1w1CqFbCSagiTU1N6ejoGH481NxnZFOgpjJnBQMATIV61MxXW7ZsSVJ6RVGp988namhoyMaNG2d7mACLhnASqkh9fX3JWZHjNQUCAIBqs379+mzbtm3CY7x/Bpg5wklgwdOdd26cWLN0qKbSyNpKapYCAABQDuEksODNp+68i7lezng1S9vb24e/V7MUAACAcggngQVPd965cWLN0vG6vQMAAMBkCSeBBU933rlRqmZpW9tMzTkFAACgGi2p9AAAAAAAgOpk5uQ819/fnyTp7Owcs6/UksqRRjapAAAAAID5Zlrh5P3335/bb789nZ2dGRwczCtf+cq85z3vyQc+8IEsWTK9SZl33nlntm/fniS58sorc+edd457bFdXV26++eb8wz/8Qw4ePJgzzzwzl112Wa677rqceuqp0xpHpQ11xh36bzEVDQ0NMzUcAABgHunp6UlfX9+Y7UMTFcabsNDQ0JCNGzfO6tiAuTPV14LE6wGVN+Vw8tZbb83OnTuTJOecc05WrVqVxx57LFdffXUeeOCB3HfffVMOKJ9++ul87GMfm9SxDz74YC699NIMDAzk9NNPz6te9ap0d3fnv/7X/5r77rsv3/rWt3LGGWdMaRzzwZYtW5IcbzJRX18/al9XV1fa29uza9euNDc3lzzfiwwAACxOPT092bRp04THtLe3j7tv7969PivAIjDd14LE6wGVNaVw8uGHH87HP/7xLFmyJLt27cq73vWuJMljjz2WX/u1X8uXv/zl3H777dmxY8eUBvXhD384zz77bC699NJ89atfHfe4vr6+/MZv/EYGBgZy9dVX57bbbsvy5cvzzDPP5K1vfWseeuihXHnllfm7v/u7KY1jPli/fn22bds24THNzc1jmlQw//X39w/PjE1K39EqFUoDAECS4VlSpSYrTFQCamiSQ6lZVsDCM9XXgmTmXg9Gfr4db7amz7eMZ0rh5M0335xisZjt27cPB5NJcv755+f222/P1q1bc+utt+ZDH/pQli9fXtZzP/DAA9m9e3fe//7352Uve9mE4eSf/umf5umnn05zc3Nuv/32LF16vOPuunXrcs899+Tcc8/NV7/61XR2dgrvmHe6u7vT2to6ZvvIO1odHR3+7gIAMKHxJiu0tbVVYDRApVTytaDU59sTZ2v6fMt4yg4nDx06lAceeCDJ8VqQJ3rHO96R97///XnmmWfy4IMP5ld/9Vcn/dyFQmE4lPyDP/iD3HHHHRMef++99yZJ3v3udw8Hk0POOuusXHTRRfn617+ev/mbv/EPgHmnqakpHR0dw49L3dFqamqq1PAAAABgUkZ+vh1vtqbPt4yn7HDy0UcfzZEjR1JXV1cy8Fu+fHle85rX5Bvf+Ea+/e1vlxVO3nzzzXn88cfzV3/1VydtZPPCCy8M/8Uf7y5AW1tbvv71r+fb3/72pMcAc6W+vn7MvyF3twEAAFhoTvx867Mt5Si7Y01PT0+S4zMTly0rnW2ec845o46djK6urvyX//Jf8su//Mv5zd/8zZMe39vbm6NHj4663lTGMTg4mEOHDo36AwAAAADMvrJnTh48eDBJsnbt2nGPGdo3dOzJFIvFvO9978uLL76YP/7jPy5rHBONZTLjuOWWW3LjjTdO6poAALCYnNigb7wyMxoYAEyN11k4ubLDyUKhkCQ55ZRTxj2mtrY2yfF/dJNx11135b//9/+eHTt25NWvfnVZ45hoLJMZx86dO3PNNdcMPz506FA2bNgwqTEAAMBCNl6DvpE0MACYOq+zcHJlh5N1dXVJkiNHjox7zODgYJKUbFN/oqeffjof+9jH8opXvCLXX3992eMYGsvIx+WMo7a2djjEBACAanJig76urq60t7dn165daW5uHj4GgKnxOgsnV3Y4OZml0pNZ+j3kox/9aA4cOJA/+7M/y6pVq8oex9D1Xv7yl09rHIxv5DT0rq6uUV+HmIZOtenv70+SdHZ2jtk3Xne6ZOy/HQCqW6V/n5Rq0Jckzc3NZvEAzACvsyxkxUIhz+/Ycfz7n04AnA1lh5MbN25Mkjz55JN54YUXSjbF+eEPfzjq2Ik8+uijSZIPfvCD+eAHPzhq3/PPP58kueeee/J3f/d3SZKnnnoqSdLY2Jjly5fn6NGj+eEPf1gynCxnHIyv1DT09vb2UY9NQ6faDAX227dvn9L5DQ0NMzkcABYov08AgGpXdji5efPmLF++PIVCIZ2dnfmlX/qlUfuPHj2a73znO0mSCy+8cNLPu3///nH3DQwMjKkbuWzZsrS0tOTb3/52HnrooZJt6h966KGyx8FYI6ehj3cH3zR0qs2WLVuSlJ41XGqpxkgNDQ1umgCQxO8TAICyw8nVq1fnoosuyte+9rXcddddY8LJL3zhCzl06FDWrVuXN77xjSd9vu9+97vj7rvhhhty44035sorr8ydd945Zv/b3va2fPvb385nPvOZfOQjH8nSpUuH9z355JN54IEHkiRvf/vbJ/fDUdKJ09BLBcFQbdavX59t27ZNeIylGgCcjN8nAEC1WzKVk6699trU1NTkzjvvzOc+97nh7Y899thw5+uPfvSjo7po33HHHWlsbMw73/nOaQ75Jb/927+d9evXp6urK9dcc02OHj2aJHnmmWdy+eWX54UXXsjFF1980s5YAAAAlKe/vz+dnZ3Dfx566KHs3r07Dz300PC2obqqADCesmdOJsdnzt1000257rrrcvnll+e6667LqlWrsmfPnrz44ou59NJL85GPfGTUOc8++2x+9KMfpbGxcSbGneT4LM7Pf/7zectb3pJPf/rT+dznPpezzjorXV1d6e/vT2NjY+6+++4Zux4AAADHlapNfyK16QEWhmKhUGLbS01wikeOlDympq5u2teeUjiZHJ89ef755+cP//AP09HRkaeeeirnnXde3vOe9+SDH/zgqCXWs+lXfuVX8sgjj+Tmm2/OP/zDP+R73/tezjzzzFx22WW57rrrZrRT90TdFBMdegEAgOoxsjZ9UrpOqtr0AAtD3+VbJ9zff+ON6Tt1bMa2+t4vTvvaUw4nk+Qtb3lL3vKWt0zq2BtuuCE33HBDWc8/2XNe9apXjVpePlum200x0VERAACYGQNHjyVrN6Rr/+HkiacnfV7X/sPJ2g3Hz5+GE2vTD1EnFYByTCucrDYTdVNMdFQEAADmTu+BQnLJzrTv3pNkT3knX7IzvQcK0eoSgCRpuGf3mG3FwmDytsuSJPXXX5+GCy+clWsLJ8swmW6KiTuFAADA7Gs8rS65/5bs2r07zWUsn+7q7k771q1pvGr6S/EAWBxOVjuy5pRTZqS+ZCnCSWDBm6gerFqwAMBitWL50uTgvjSfsTItZ58++RMP7ksO7jt+PgBUmHASWPCmWw92IdaC7e/vH/65k5eC1pGB63glKGAxc7MCAAAWFuEksOBNVA92sdaC7e7uTmtr65jt7e3tw993dHQoMUHVqcabFQAAsJAJJ4EFbzL1YBdbLdimpqZ0dHQMPy41I6ypjNpTsFhU480KAABYyISTAAtQfX39mLC1rU2/TajGmxVzpaenJ319fWO2lyorMZLQFwBgYaqpq8uq225LWltTU1s7a9cRTgIAMKGenp5s2rRpwmNGlpU40d69ewWUAACUJJwEgAob2eBovFloGhxRSUMzJk9cEl8cHMwzv/u7ebK/P82f/GRWrFkz6ryhpfSlZlwCAEAinISqVSwU8vyOHce/Hxys8GigupVqcHTiLDQNjpgPTlwSXywU0rdufV67Lml43etSU1dXwdEBQHUaOHosWbshXfsPJ088Xda5XfsPJ2s3HH8OqBDhJABU2MgGR6WaGw0dA7AYmC0OMLN6DxSSS3amffeeJHvKf4JLdqb3QCEq2Fev8WqLJ3NTX1w4CQAVdmKDI82NgMXMbHGAmdV4Wl1y/y3ZtXt3msu8od3V3Z32rVvTeNUXZ2l0zHeTqS2ezG59ceEkAAAwZ8wWB5hZK5YvTQ7uS/MZK9Ny9unlnXxwX3Jw3/HnoCqNV1t8yHi/q5OZqy8unAQAYFKKg4MpFgovPS4Mlvx+5PFwIrPFAWD+ObG2+Eiz/btaOAlVYOQHyZe2jfhAeeRIyWM0NgBgpMPXXpu+U9eW3Pf8FVeMPf7Zg9O63nj1j+ai9hGLizqXADB/CSehCvRdvnXC/f033ljyw+bqe9UdAaAyJlP/aDZrH7G4LNY6l/39/UmSzs7OMftOtgwPAOYL4SQAAJOy8hOfSMPmzcOPi4XB4RmTq+6+OzV1taOPf/TR5PWvn9K1Jqp/NBe1j1hcFmudy6HZoNu3b5/S+Q0NDTM5HACYEuEkVIGGe3aP2VYsDCZvuyxJUn/99Wm48MK5HhYAC0xNbe24JT9q6sbuq6mtLXlsOcarf6ROIeVYrHUut2zZkqT0kvShoH68BgfKHwAwXwgnoQqcrHZkzSmnqC8JALDArF+/Ptu2bZvwmIkaHADAfCCcpGoNHD2WrN2Qrv2HkyeenvR5XfsPJ2s3HD8fAAAAgCkTTlK1eg8Ukkt2pn33niR7yjv5kp3pPVDI4lgQBAAA81dPT8+4dWTH674+xPJ1gPlPOEnVajytLrn/luzavTvNZRRA7+ruTvvWrWm8Sifryejv7x8u1j7em8dSdZIAAKCnpyebNm066XEndl8fae/evQJKgHlMOEnVWrF8aXJwX5rPWJmWs0+f/IkH9yUH9x0/n5Pq7u5Oa2vrqG0nvnns6OhQCwlgAaqpq8vqe92sA2bP0IzJ8Rr7jNd9PXmpKdB4sy4BmB+Ek1Claurqsuq225LW1hnppjqepqamdHR0JBn/zWNTGTNXAQCoPhM19lks3dcBTjRyJWJSejXiYliJKJwEZlV9ff2oN5LePAIAAMDJlVqJmIxejbgYViIKJwEAAIA5deKMsIlWWS30WWEwVSNXIial/50shpWIwkkAAABgTo03I+xEi2FWGEzViSsRk8W5GlE4CQAAAMypE2eEDTUwOrH50WKYFQZMTDgJAAAAzKlSM8KSiZsfAYvTkkoPAAAAAMpVLBRy6G1vz6G3vT3FQqHSwwFgioSTAAAAAEBFCCcBAAAAgIpQcxIAAADG0d/fn+7u7iTHm7aM/Dqkqakp9fX1cz42gMVAOAkAACcYGUYkAgmoZt3d3WltbR21rb29fdTjjo4OTVwApkg4CQAAJygVRiQCCaiU4uDgmKY3xcJgye9HnjMTmpqa0tHRkSQZGBhIb29vGhsbs2LFilHHADA1wkkAADjByDAiEUhApR2+9tr0nbp23P3PX3HF2HOePTgj166vrx91E6KtrW1GnhdgPhg4eixZuyFd+w8nTzxd1rld+w8nazccf45pEE4CAMAJTgwjEoEEALD49B4oJJfsTPvuPUn2lP8El+xM74FCpvMuSTgJAADAvLbyE59Iw+bNo7YVC4PDMyZX3X13aupqR5/z6KPJ618/Z2MEWIgaT6tL7r8lu3bvTnOZK0K6urvTvnVrGq/64rTGIJwEAABgXquprU1NXd34++vG7q+prR3naACGrFi+NDm4L81nrEzL2aeXd/LBfcnBfcefYxqWTOtsAAAAql6xUMjzO3Yc/36GGtEAUB3MnFxg+vv7093dnSTp6uoa9XVIU1NT6uvr53xsAAAAAFAO4eQC093dndbW1lHb2tvbRz3u6OgYU8AdAAAAWHz6+/uTJJ2dnWP2DQwMpLe3N42NjVmxYsWY/SdOdoJKEE4uME1NTeno6Egy/otMU5kFTAEAJjJw9FiydkO69h9Onnh60ud17T+crN1w/HwAYFYMra7cvn37lJ+joaFhpoYDZRNOLjD19fWjZkW2tU2nWTsAwMn1Higkl+xM++49SfaUd/IlO9N7oBDvWICZVlNXl9X3Tq9DLC8ZWUJsookwSojNP1u2bElS+v9PV1dX2tvbs2vXrjQ3N5c8v6GhIRs3bpztYcK4hJMAAEyo8bS65P5bsmv37jSXsUKjq7s77Vu3pvEq4QEsJsVCocS2l5rgFI8cKXnMRN22qbxSJcROpITY/LR+/fps27ZtwmOam5v9v2PeEk4CADChFcuXJgf3pfmMlWk5+/TJn3hwX3Jw3/HzgUWj7/KtE+7vv/HG9J26dsx2sxznt5ElxMabbaeEGDAbhJMAAABQ5U4sIZaYbQfMDeEkAAAAk9Zwz+4x24qFweRtlyVJ6q+/Pg0XXjjXwwJggRJOAgAAMGknqx1Zc8op6ksCMGnCSQAAAOalgaPHkrUb0rX/cPLE02Wd27X/cLJ2w/HnAGDeEk4CAAAwL/UeKCSX7Ez77j1J9pT/BJfsTO+BQtpmfGQAzBThJLDo9Pf3p7u7O8nxToMjvw5pampKfX39nI8NAIDJazytLrn/luzavTvNZXaK7uruTvvWrWm8SpdwgPlMOAksOt3d3WltbR21rb29fdTjjo4OnQcB5rGpLuW0jBMWlxXLlyYH96X5jJVpOfv08k4+uC85uO/4cwAwbwkngUWnqakpHR0dSZKBgYH09vamsbExK1asGHUMAPPXtJZyWsYJc66mri6rbrstaW1NTW1tpYcDwAIinAQWnfr6+lGzItvafDwFWGimupTTMk4AgIVFOAkAwLwz5aWclnECACwoSyo9AAAAAACgOgknAQAAAICKEE4CAAAAABUhnAQAAAAAKkI4CQAAAABUhHASAAAAAKgI4SQAAAAAUBHCSQAAAACgIpZVegAAAABMT39/f7q7u4cfd3V1jfo6pKmpKfX19XM6NgCYiHByGibzBsAvfwAAYLZ1d3entbV1zPb29vZRjzs6OtLS0jJXwwKAkxJOTsNk3gD45Q8AAMy2pqamdHR0DD8eGBhIb29vGhsbs2LFilHHAcCQ/v7+JElnZ2fJ/eP9PknGzs6fKuHkNEzmDYBf/gAAwGyrr68fMymira2tQqMBYKEYWhG8ffv2KT9HQ0PDtMYgnJwGbwAAAIBq09PTk76+vjHbx6tzOaShoSEbN26c1bEBUJ4tW7YkGb8sYVdXV9rb27Nr1640NzeP2T8Tr+3CSQAAACalp6cnmzZtmvCYE+tcjrR3714BJcA8sn79+mzbtu2kxzU3N89a2ULhJAAAzAOaLbIQDM2YLDWD5mR1ydrb20vOuASgugknAYBFSdDDQqPZIgvJeDNolLkCoFzCSQBgURL0sNBotggAVCPhJACwKAl6WGg0WwQAqpFwEgBYlAQ9AEAljNfRPtHVHkoRTgIAAADMgMl0tE90tYeRhJMAAAAAM2CijvaJrvZQinASAACqiOWGALNvvI72iTIzcCLhJAAAVAnLDQGA+UY4CQAAVcJyQwBgvhFOAgBAlbHcEACYL4STsMipKwXAdPX39ydJOjs7x+w72Uw7AACYiHASFjF1pQCYCd3d3UmS7du3T+n8hoaGmRwOAACLiHCSqlUNs0DUlQJgJmzZsiVJ0tTUlPr6+lH7hn5fjPe7xix8AAAmIpykalXTLBB1pQCYjvXr12fbtm0THjPR7xoAABiPcJKqZRYIAAAAQGUJJ6laZoEAwPxVDeVXAAAQTgIAMA9VU/kVAIBqJpwEAGDeUX4FAKA6CCcBAJh3lF8BAKgOwkkAAADmpYnqzyZq0AIsBsJJgGnq7+8fro2WvPRG+MQ3xKWWJgIAML7p1p9N1KAFmO+EkwDT1N3dndbW1jHb29vbRz3u6Oiw/BAAoAwT1Z9N1KAFWAyEkwDT1NTUlI6OjuHH4y0vampqqsTwAGBSioVC+i7fmiRpuGd3aurqKjwimFz92UQNWoCFTDgJME319fVj3gy3tbVVaDQAAACwcCyp9AAAAAAAgOoknAQAAAAAKkI4CQAAAABUhJqTAABQZYqDgykWCqO3FQZLfj/yHACAmSacBACAKnP42mvTd+racfc/f8UVY8959uBsDgmoEj09Penr6xuzvaura9TXEzU0NGTjxo2zOjagMoSTAAAAwKzr6enJpk2bJjymvb193H179+4VUMIiJJwEAIAqs/ITn0jD5s2jthULg8MzJlfdfXdq6mpHn/Poo8nrXz9nYwQWn6EZk7t27Upzc/OofQMDA+nt7U1jY2NWrFgxal9XV1fa29tLzrgEFj7hJAAAVJma2trU1NWNv79u7P6a2tpxjgYoT3Nzc1paWsZsb2trq8BogEoTTgIAAADMgIGjx5K1G9K1/3DyxNNlndu1/3CydsPx54AqIpwEAAAAmAG9BwrJJTvTvntPkj3lP8ElO9N7oBBzSKkmwkkAAACAGdB4Wl1y/y3ZtXt3mpuayjq3q7s77Vu3pvGqL87S6GB+Ek4CAAAwKVNdsmq5KtVixfKlycF9aT5jZVrOPr28kw/uSw7uO/4cUEWEkwAAQGrq6rL6XrN1mNi0lqxargpACcJJAAAAJmWqS1YtVwVgPMJJAAAAJmXKS1YtVwVgHEsqPQAAAAAAoDoJJwEAAACAihBOAgAAAAAVIZwEAAAAACpCOAkAAAAAVIRwEgAAAACoiGWVHgAAADA3Bo4eS9ZuSNf+w8kTT5d1btf+w8naDcefAwBghggnAQCgSvQeKCSX7Ez77j1J9pT/BJfsTO+BQtpmfGQAQLUSTgIAQJVoPK0uuf+W7Nq9O81NTWWd29XdnfatW9N41RdnaXQAQDUSTgIAQJVYsXxpcnBfms9YmZazTy/v5IP7koP7jj8HAMAM0RAHAAAAAKiIaYWT999/fy666KKcdtppWblyZVpaWvJHf/RHefHFF8t6ngcffDBXX311Xve61+XMM89MbW1tGhoa0tramptuuil9fX0lz/vMZz6TmpqaCf98/etfn86PCAAAAADMkikv67711luzc+fOJMk555yTVatW5bHHHsvVV1+dBx54IPfdd1+WLJlc9nnXXXdl9+7dWbZsWX72Z382v/ALv5Cnn346jz76aDo7O/OXf/mX+cd//MecddZZJc9/2ctelo0bN5bct3bt2qn9gAAAAADArJpSOPnwww/n4x//eJYsWZJdu3blXe96V5Lksccey6/92q/ly1/+cm6//fbs2LFjUs932WWXpb29PW94wxuyYsWK4e3/+q//mne9613553/+57z//e/PV7/61ZLnX3zxxfnMZz4zlR8FAAAAAKiQKS3rvvnmm1MsFrNt27bhYDJJzj///Nx+++1Jjs+sPHr06KSe7+1vf3ve/OY3jwomk+Tnf/7nc+eddyZJ/v7v/z6FQmEqwwUAAAAA5qGyZ04eOnQoDzzwQJLkyiuvHLP/He94R97//vfnmWeeyYMPPphf/dVfndYAm5qakiTHjh3L4OBg6urqpvV8AAAAwNwbOHosWbshXfsPJ088PenzuvYfTtZuOH4+J9Xf35/u7u4kSVdX16ivQ5qamlJfXz/nY4NSyg4nH3300Rw5ciR1dXVpaWkZs3/58uV5zWtek2984xv59re/Pe1w8uGHH05yvK7lmjVrSh7z2GOP5fLLL89TTz2V1atXZ/PmzWlvb8+55547rWsDAAAAM6P3QCG5ZGfad+9Jsqe8ky/Zmd4DhbTNysgWl+7u7rS2to7a1t7ePupxR0dHyUwHKqHscLKnpydJctZZZ2XZstKnn3POOfnGN74xfGy5isVi9u/fn2984xv5nd/5nSxbtmx4uXgp3/3ud/Pd7353+PGXvvSl3HTTTbnxxhtz7bXXTmkMAAAAwMxpPK0uuf+W7Nq9O80/XSU5GV3d3WnfujWNV31xFke3eDQ1NaWjoyNJMjAwkN7e3jQ2No4qpddUxn9/mG1lh5MHDx5MMnEX7KF9Q8dO1t/+7d/msssuG7XtDW94Q77whS+krW3s/ZFTTz01V111Vd75znfm537u57JmzZp0dXXl9ttvz1//9V/nuuuuy5o1a/LBD35w3GsODg5mcHBw+PGhQ4fKGjMAAABwciuWL00O7kvzGSvTcvbpkz/x4L7k4L7j53NS9fX1o2ZFlspTYD4pO5wcakpzyimnjHtMbW1tkuMJfTnWrVuXtra2HDt2LPv27cuPf/zj/I//8T/y2c9+Ni0tLWMa5mzZsiVbtmwZte0Xf/EX89nPfjbr1q3LHXfckeuuuy6/9Vu/lYaGhpLXvOWWW3LjjTeWNU5YKKZa0yVR1wUAAACYfWWHk0MNaY4cOTLuMUMzEU8ME0/ml3/5l/PNb35z+HFXV1c+8IEP5M///M/z5JNP5mtf+9qkn+vGG2/Mn/zJn+S5557LP/zDP+Stb31ryeN27tyZa665ZvjxoUOHsmHDhrLGDfPVtGq6JOq6AAAAALOq7HByMku2J7P0ezKam5vzla98Jeeee26+/vWv55vf/Gb+t//tf5vUuatXr86rXvWqdHZ25vHHHx/3uNra2uGZnrDYTLWmS6KuCwAAADD7yg4nN27cmCR58skn88ILL5RsivPDH/5w1LHTsXLlyrzxjW/M//1//9/p7OycdDiZHO8cniQvvPDCtMcBC9GUa7ok6roAAAAAs25JuSds3rw5y5cvT6FQSGdn55j9R48ezXe+850kyYUXXjj9EealcLGckPHYsWP5/ve/nyR5xSteMSPjAAAAAABmTtkzJ1evXp2LLrooX/va13LXXXfll37pl0bt/8IXvpBDhw5l3bp1eeMb3zjtAT733HN58MEHkxxvdjNZd911V5599tksXbp0RsYBAABQ7fr7+5Ok5ESVgYGB9Pb2prGxcUz/ga6urjkZHwALT9nhZJJce+21+frXv54777wzb3zjG/Oud70rSfLYY48NN5f56Ec/Oqqj9x133JE77rgjr33ta/P5z39+ePuPf/zjfPKTn8z27dvzqle9atR1/umf/ikf/vCHc+DAgZx33nl5wxveMLzv0KFDed/73pcPf/jDowLSY8eO5e67786HPvShJMmVV16ZM888cyo/JgAAACN0d3cnSbZv3z6l8xsaGmZyOAAsAlMKJ9va2nLTTTfluuuuy+WXX57rrrsuq1atyp49e/Liiy/m0ksvzUc+8pFR5zz77LP50Y9+lMbGxlHbjxw5kk996lP51Kc+ldNOOy2NjY0pFovZt29ffvKTnyRJzj333Nx3331ZuvSl2ncvvvhiPv/5z+fzn/98Tj311Jx99tlZtmxZenp68uyzzyZJLr744nzqU5+ayo8IAADACbZs2ZIkaWpqSn19/ah9XV1daW9vz65du9Lc3Dzm3IaGhhnpSwDA4jKlcDI5Pnvy/PPPzx/+4R+mo6MjTz31VM4777y85z3vyQc/+MFRQeJEfuZnfiZ/9md/lm984xv57ne/mx/84Ac5fPhw1q5dmze96U3ZsmVLtm3bNmZZwMqVK/PJT34y3/rWt7Jnz5784Ac/yMDAQNatW5dLL700v/mbv5l3vOMdqampmeqPCAAAwAjr16/Ptm3bJjymubk5LS0tczQiABa6KYeTSfKWt7wlb3nLWyZ17A033JAbbrhhzPa6urq8973vzXvf+96yrr18+fL8zu/8TlnnAAAAAADzR9ndumGxKxYKeX7HjuPfDw5WeDQAAAAAi5dwEgAAAACoCOEkAAAAAFAR06o5CQAAAAAsDv39/enu7h5+3NXVNeprkjQ1NaW+vn7GrimcBAAAAJgB/f39SZLOzs6S+wcGBtLb25vGxsasWLFi1L6R4Q9USnd3d1pbW8dsb29vH/6+o6MjLS0tM3ZN4SRVrVgolNj2UhOc4pEjJY+pqaub1XEBAABjVWJGD5Rj6O/n9u3bp/wcDQ0NMzUcKFtTU1M6OjqGH5cK1Juammb0msJJqlrf5Vsn3N9/443pO3XtmO2r7/3ibA0JAAAYRyVm9EA5tmzZkmT8kLyrqyvt7e3ZtWtXmpubx+xvaGjIxo0bZ3uYMK76+voxr6FtbW2zek3hJAAAAAtCJWb0QDnWr1+fbdu2nfS45uZmITr8lHCSqtZwz+4x24qFweRtlyVJ6q+/Pg0XXjjXwwIAAEqoxIweAGaXcJKqdrLakTWnnKK+JAAAsCj19PSkr69vzPZStTxHsvQYmEnCSQAAAKgyPT092bRp04THjKzleaK9e/cKKIEZIZwEAACAKjM0Y7JUY5ZStTyHDDV0KTXjEmAqhJMAAABQpcZrzKKWJzBXhJMAAFAl+vv7kySdnZ0l959sthSw+BULhfRdvjXJ8QaiavADs004CSeoqavLqttuS1pbU1NbW+nhAADMmO7u7iTJ9u3bp/wcDQ0NMzUcAADhJAAAVIstW7YkSZqamlJfXz9m/1AtuVI16BIdegGAmSecBACgLP39/cMz8IaW+p645He88IvKWr9+fbZt23bS48arQQcAMNOEkwAAlKW7uzutra2jtrW3t4963NHRIdwCAOCkhJMAAJSlqakpHR0dScZvoNLU1FSp4QFQhuLgYIqFwkuPC4Mlvx95PMBMEk4CAFCW+vr6UbMi29raKjgaAKbj8LXXpu/UtSX3PX/FFWOPf/bgbA8JqDJLKj0AAAAAAKA6mTkJAAAAVWrlJz6Rhs2bhx8XC4PDMyZX3X13aupqRx//6KPJ618/p2MEFjfhJAAAAFSpmtra1NTVld5XN3ZfTW1tyWMBpko4CQAAAMy6/v7+JElnZ+eYfeM1WEuSrq6uORkfUBnCSQAAAGDWdXd3J0m2b98+pfMbGhpmcjjAPCGcBAAAAGbdli1bkiRNTU2pr68fta+rqyvt7e3ZtWtXmpubx5zb0NCQjRs3zsUwgTkmnAQAAABm3fr167Nt27YJj2lubk5LS8scjQiYD4STAAAAQJKkpq4uq+/9YqWHAVSRJZUeAAAAAABQnYSTAAAAAEBFCCcBAAAAgIoQTgIAAAAAFSGcBAAAAAAqQjgJAAAAAFSEcBIAAAAAqIhllR4AMHv6+/uTJJ2dnSX3DwwMpLe3N42NjVmxYsWofV1dXbM+PgAAAKC6CSdhEevu7k6SbN++fcrP0dDQMFPDAQAAABhFOAmL2JYtW5IkTU1Nqa+vH7O/q6sr7e3t2bVrV5qbm8fsb2hoyMaNG2d7mAAAAECVEk7CIrZ+/fps27btpMc1NzenpaVlDkYEAAAA8BINcQAAAACAihBOAgAAAAAVIZwEAAAAACpCOAkAAAAAVISGOAAAAFBlBo4eS9ZuSNf+w8kTT0/6vK79h5O1G46fDzADhJMAAABQZXoPFJJLdqZ9954ke8o7+ZKd6T1QSNusjAyoNsJJAAAAqDKNp9Ul99+SXbt3p7mpadLndXV3p33r1jRe9cVZHB1QTYSTAAAAUGVWLF+aHNyX5jNWpuXs0yd/4sF9ycF9x88HmAEa4gAAAAAAFSGcBAAAAAAqQjgJAAAAAFSEcBIAAAAAqAjhJAAAAABQEcJJAAAAAKAihJMAAAAAQEUsq/QAAABgIv39/enu7h5+3NXVNerrkKamptTX18/p2AAAmB7hJAAA81p3d3daW1vHbG9vbx/1uKOjIy0tLXM1LAAAZoBwEgCAea2pqSkdHR3DjwcGBtLb25vGxsasWLFi1HEAACwswknI5JaLWSoGAJVRX18/ZkZkW1tbhUYDAMBMEk5CJrdczFIxAAAAgJklnIRMbrmYpWIAAAAAM0s4CbFcDAAAAKASllR6AAAAAABAdRJOAgAAAAAVIZwEAAAAACpCOAkAAAAAVIRwEgAAAACoCOEkAAAAAFARwkkAAAAAoCKEkwAAAABARQgnAQAAAICKEE4CAAAAABUhnAQAAAAAKkI4CQAAAABUhHASAAAAAKgI4SQAAAAAUBHCSQAAAACgIoSTAAAAAEBFCCcBAAAAgIoQTgIAAAAAFSGcBAAAAAAqQjgJAAAAAFTEskoPAAAAgIWnv78/3d3dw4+7urpGfU2Spqam1NfXz/nYAFg4hJMAAACUrbu7O62trWO2t7e3D3/f0dGRlpaWuRwWAAuMcBIAAICyNTU1paOjY/jxwMBAent709jYmBUrVgwfw/zU39+fJOns7Byzr9T/yyEjZ8YCzAThJAAAAGWrr68fMyuyra2tQqOhXENL8rdv3z6l8xsaGmZyOEAVE04CAABAldmyZUuS0nVBu7q60t7enl27dqW5uXnMuQ0NDdm4ceNcDBOoAsJJAAAAqDLr16/Ptm3bJjymublZzVBg1i2p9AAAAAAAgOpk5iQAAADALOjv7x+u75m81FBoZGOhUkvroZoIJwEAAABmQXd3d1pbW8dsb29vH/6+o6PD8nmqmnASAAAAYBY0NTWlo6Nj+PHAwEB6e3vT2NiYFStWDB8D1Uw4CQAAADAL6uvrx8yKbGtrq9BoYH7SEAcAAAAAqAjhJAAAAABQEcJJAAAAAKAihJMAAAAAQEUIJwEAAACAihBOAgAAAAAVIZwEAAAAACpCOAkAAAAAVIRwEgAAAACoCOEkAAAAAFARwkkAAAAAoCKEkwAAAABARQgnAQAAAICKEE4CAAAAABUhnAQAAAAAKkI4CQAAAABUhHASAAAAAKgI4SQAAAAAUBHCSQAAAACgIoSTAAAAAEBFCCcBAAAAgIoQTgIAAAAAFSGcBAAAAAAqQjgJAAAAAFSEcBIAAAAAqIhllR4AAAAAUF36+/vT3d09/Lirq2vU1yFNTU2pr6+f07EBc0s4CQAAAMyp7u7utLa2jtne3t4+6nFHR0daWlrmalhABQgnAQAAgDnV1NSUjo6O4ccDAwPp7e1NY2NjVqxYMeo4YHETTgIAQJWazLJKSyqB2VBfXz9mRmRbW1uFRgNUknASAACq1GSWVVpSCQDMpml1677//vtz0UUX5bTTTsvKlSvT0tKSP/qjP8qLL75Y1vM8+OCDufrqq/O6170uZ555Zmpra9PQ0JDW1tbcdNNN6evrm/D8rq6ubN26NS9/+ctTV1eXc889Nzt27Mizzz47jZ8OAAAWt6FllUN/vvnNb2bXrl355je/ObzNkkoAYDbVFIvF4lROvPXWW7Nz584kyTnnnJNVq1Zlz549efHFF/Prv/7rue+++7JkyeSyz/b29uzevTvLli3Lz/7sz+ZlL3tZnn766Tz55JMpFos5++yz84//+I8566yzxpz74IMP5tJLL83AwEBOP/30bNiwId3d3env788555yTb33rWznjjDMm/XMdOnQoa9asyXPPPZfVq1dP+jxYiDo7O9Pa2mpGBAAAMMznBGAuTWnm5MMPP5yPf/zjWbJkSe6555784Ac/yGOPPZbOzs6cccYZ+fKXv5zbb7990s932WWX5Wtf+1oOHTqUH/3oR/nOd76T3t7e7NmzJ7/wC7+QJ554Iu9///vHnNfX15ff+I3fyMDAQK6++ur8z//5P9PR0ZEnn3wybW1t+eEPf5grr7xyKj8iAAAAADDLphRO3nzzzSkWi9m2bVve9a53DW8///zzh0PJW2+9NUePHp3U87397W/Pm9/85lEduZLk53/+53PnnXcmSf7+7/8+hUJh1P4//dM/zdNPP53m5ubcfvvtWb58eZJk3bp1ueeee7Js2bJ89atfTWdn51R+TAAAAABgFpUdTh46dCgPPPBAkpSclfiOd7wjq1evzjPPPJMHH3xw2gMcqnFz7NixDA4Ojtp37733Jkne/e53Z+nSpaP2nXXWWbnooouSJH/zN38z7XEAAAAAADOr7HDy0UcfzZEjR1JXV1ey9sTy5cvzmte8Jkny7W9/e9oDfPjhh5Mcr2u5Zs2a4e0vvPBCOjo6kiRtbW0lzx3aPhPjAAAAAABmVtnhZE9PT5LjMxOXLVtW8phzzjln1LHlKhaLeeqpp7J79+68+93vzrJly8bUsOzt7R1eNj50vamMY3BwMIcOHRr1BwAAAACYfWWHkwcPHkySrF27dtxjhvYNHTtZf/u3f5uamposWbIkL3/5y9Pe3p5NmzblH//xH/PWt7615DgmGstkxnHLLbdkzZo1w382bNhQ1pgBAAAAgKkpO5wcakpzyimnjHtMbW1tkmRgYKCs5163bl3a2try2te+NmeeeWZqamryP/7H/8hnP/vZMc81sjnOeGOZzDh27tyZ5557bvjPvn37yhozAAAAADA1pddlT6Curi5JcuTIkXGPGWpcc2L37ZP55V/+5Xzzm98cftzV1ZUPfOAD+fM///M8+eST+drXvjZmHENjGfm4nHHU1tYOh5gAAAAAwNwpe+bkZJZKT2bp92Q0NzfnK1/5Ss4444x8/etfHxVcjnzu8cYyU+MAAAAAAGZe2eHkxo0bkyRPPvlkXnjhhZLH/PCHPxx17HSsXLkyb3zjG5MknZ2dw9sbGxuzfPnyUdebzXEAAAAAADOr7HBy8+bNWb58eQqFwqiwcMjRo0fzne98J0ly4YUXTn+EyXAIOjIMXbZsWVpaWpIkDz30UMnzhrbP1DgAAAAAgJlTdji5evXqXHTRRUmSu+66a8z+L3zhCzl06FDWrVs3PONxOp577rk8+OCDSZJf/MVfHLXvbW97W5LkM5/5TI4dOzZq35NPPpkHHnggSfL2t7992uMAAAAAAGZW2eFkklx77bWpqanJnXfemc997nPD2x977LFcc801SZKPfvSjo7po33HHHWlsbMw73/nOUc/14x//OP/X//V/5V/+5V/GXOef/umf8uY3vzkHDhzIeeedlze84Q2j9v/2b/921q9fn66urlxzzTU5evRokuSZZ57J5ZdfnhdeeCEXX3xxWltbp/JjAgAAAACzaErhZFtbW2666aa8+OKLufzyy3Puuefm/PPPT0tLS/bv359LL700H/nIR0ad8+yzz+ZHP/pRnnrqqVHbjxw5kk996lN59atfnXXr1qW1tTUtLS05/fTT87rXvS7/9E//lHPPPTf33Xdfli5dOurc1atX5/Of/3zq6ury6U9/OmeeeWYuuOCCnHXWWXnooYfS2NiYu+++eyo/IgAAAAAwy6YUTibHZ09+5StfyZve9KY888wzefzxx3PeeefljjvuyJe+9KUxQeJ4fuZnfiZ/9md/lv/j//g/sn79+vzgBz/I9773vdTU1ORNb3pTPv3pT+d73/tezj333JLn/8qv/EoeeeSRvPOd70xNTU2+973v5Ywzzsg111yTzs7O/MzP/MxUf0QAAAAAYBbVFIvFYqUHMZ8cOnQoa9asyXPPPZfVq1dXejgwqzo7O9Pa2pqOjo7hBlMAAEB18zkBmEtTnjkJAAAAADAdwkkAAAAAoCKEkwAAAABARQgnAQAAAICKEE4CAAAAABUhnAQAAAAAKkI4CQAAAABUhHASAAAAAKgI4SQAAAAAUBHCSQAAAACgIoSTAAAAAEBFCCcBAAAAgIoQTgIAAAAAFSGcBAAAAAAqQjgJAAAAAFSEcBIAAAAAqAjhJAAAAABQEcsqPQAAAACgsvr7+9Pd3Z0k6erqGvV1SFNTU+rr6+d8bMDiJpwEAACAKtfd3Z3W1tZR29rb20c97ujoSEtLy1wOC6gCwkkAAACock1NTeno6EiSDAwMpLe3N42NjVmxYsWoYwBmmnASAAAAqlx9ff2oWZFtbW0VHA1QTTTEAQAAAAAqQjgJAAAAAFSEcBIAAAAAqAjhJAAAAABQEcJJAAAAAKAihJMAAAAAQEUIJwEAAACAihBOAgAAAAAVIZwEAAAAACpiWaUHAMyd/v7+dHd3Dz/u6uoa9TVJmpqaUl9fP+djAwAAAKpPTbFYLFZ6EPPJoUOHsmbNmjz33HNZvXp1pYcDM6qzszOtra0THtPR0ZGWlpY5GhEAAABQzcychCrS1NSUjo6O4ccDAwPp7e1NY2NjVqxYMXwMAAAAwFwwc/IEZk4CAAAAwNzQEAcAAAAAqAjhJAAAAABQEcJJAAAAAKAihJMAAAAAQEUIJwEAAACAihBOAgAAAAAVIZwEAAAAACpCOAkAAAAAVIRwEgAAAACoCOEkAAAAAFARwkkAAAAAoCKEkwAAAABARQgnAQAAAICKEE4CAAAAABUhnAQAAAAAKkI4CQAAAABUhHASAAAAAKgI4SQAAAAAUBHCSQAAAACgIoSTAAAAAEBFCCcBAAAAgIoQTgIAAAAAFSGcBAAAAAAqQjgJAAAAAFSEcBIAAAAAqAjhJAAAAABQEcJJAAAAAKAillV6APNNsVhMkhw6dKjCIwEAAACAhauhoSE1NTUTHiOcPEFfX1+SZMOGDRUeCQAAAAAsXM8991xWr1494TE1xaGpgiRJXnzxxfz4xz+eVLJ7okOHDmXDhg3Zt2/fSf/Dz4S5vp5ruuZCu55ruuZCu55ruuZCu55ruuZCu55ruuZCu55ruuZCu55ruuaJzJycgiVLluQVr3jFtJ5j9erVc/aXpBLXc03XXGjXc03XXGjXc03XXGjXc03XXGjXc03XXGjXc03XXGjXc03XLIeGOAAAAABARQgnAQAAAICKEE7OoNra2lx//fWpra1dlNdzTddcaNdzTddcaNdzTddcaNdzTddcaNdzTddcaNdzTddcaNdzTdecCg1xAAAAAICKMHMSAAAAAKgI4SQAAAAAUBHCSQAAAACgIoSTAAAAAEBFCCcBAAAAgIoQTk7DE088kb/4i7/I9u3bc/7552fZsmWpqanJzTffPCvXKxaL+eY3v5nf+Z3fyWtf+9qceuqpOeWUU/KzP/uzefvb354HH3xwVq77t3/7t3nf+96X1tbWvPzlL88pp5ySU089Na9//evzqU99KkeOHJmV64503XXXpaamZtb++7773e8efv7x/hQKhRm/bpIcO3Ysf/EXf5E3vOENWb9+ferq6vLv/t2/y5YtW/KlL31pxq7T29t70p9x6M9/+2//bcaumyT/9m//lh07duRVr3pV6uvrU1dXl3PPPTfvfe978/jjj8/otYY89dRT+fCHP5yNGzemrq4u69evz5vf/Ob8/d///ZSebzr/3h9++OG89a1vzemnn54VK1bk53/+53PTTTed9O/UVK751FNP5bOf/Ww++MEP5pd+6ZdSW1ubmpqabNu2bdZ+zkcffTS///u/P/x3ePny5XnZy16Wiy++OPfdd9+sXPPBBx/M1Vdfnde97nU588wzU1tbm4aGhrS2tuamm25KX1/fjF+zlDvvvHP4383J/htP5Zo33HDDSf+9dnd3z8rP+P/8P/9P3vzmN+eMM85IbW1tzjzzzLz5zW/O3XffPaM/42Rfl/7qr/5qRn/Ovr6+/Kf/9J+yefPmrFq1KqecckrOOuusbN26NZ2dnROeO9VrPvfcc/n93//9vPrVr059fX1OPfXU/O//+/+ez33ucxOeN93f/1N5DZrqNafzGjTVa071NWiq15vO689Mvpeb7OvPVK85ndefmfg5y30Nmuo1p/MaNJ2fcyqvQdO53lRff5LpfRaY6nugqVxzuu+BpnLN6b4Hmso1p/seaKY+25XzHmgq15zOa9B0f8apvAeayjWn+x5oqj/ndN4DTfWa03kNOtFks4Gpvv5M5ZrTff2ZyjWn+/ozlWtO9/VnQkWm7EMf+lAxyZg/N91006xc74EHHhi+xpIlS4qbNm0qbt68ubhq1arh7dddd92MX7etra2YpFhbW1s8++yzixdccEHxzDPPHL5ma2tr8eDBgzN+3SH/+q//WjzllFNm9b/vb/3WbxWTFDdu3Fhsa2sr+WdwcHDGr3vgwIHia1/72mKSYk1NTfGVr3xlsbW1tfjyl7+8mKT49re/fcau9b/+1/8a92dra2srnnPOOcUkxbq6uuKzzz47Y9ft7u4uvuxlLysmKS5fvrz4yle+svjqV7+6WFdXV0xSrK+vL/7jP/7jjF2vWCwW//mf/7l4xhlnDP+9bW1tLf7cz/3c8N+hW265peznnOq/9127dhWXLl1aTFI888wzi5s3by4uX768mKT4mte8pnj48OEZveYf/uEfljznyiuvnJWf8/HHHx913Nlnn11sbW0trl27dnjbb/3WbxWPHTs2oz/n1q1bi0mKy5YtK5511lnFCy64oPjv/t2/K9bU1AyP40c/+tGMXvNE//Zv/1Y87bTTJv3feCrXvP7664tJihs2bBj33+54P+dUf8ZCoVD89V//9eHjzznnnOJrXvOa4oYNG4pLliwptra2zujPONHr0v+/vTMPiuJM//h3kFNQlGBEQEFF0THigaArnrhmvaPrUW48S7dcMWsZtyqbuO5Gy2ujRmPWKk1iLI2rJsYjGy2vqMGoeCtoPDAqsngTxQPkkOP5/UFN/4CZYbrf920m2Xo+VVTM9Ex/5u3p+fbzvt3zttVq1daRnp6uzPnw4UNq2bKldjxt3rw5tWvXTjue1qpVizZv3qy0nXfu3KEWLVpo62/Xrh1ZrVZtn506darT18oc/0UzSNQpk0EiTpkMEm2jTP6oquWM5I+oUyZ/ZNopmkGiTpkMEnWKZpCoTyZ/bNsIMN4XkKmBRJyyNZBRp4oaSKSdsjWQir6d0RpIxCmTQaJtlKmBRJyyNZCIU7YGEnHKZlBF9I4NyOSPiFM2f4w6VeSPSDtl86c6eHBSgvnz59OgQYNo3rx5tHfvXho+fHi1XxBZDhw4QFFRUbRq1SrKycnRHi8qKqJZs2ZpO9GuXbuUetetW0fJycn08uXLSo+fOHGCwsPDCQBNmzZNqdNGWVkZde/enfz9/SkxMdG07WsbnFy3bp3ydTujtLSUunXrRgDo97//Pd2+fbvS8tu3b9MPP/xQY+/HFjSjRo1Sut4+ffoQAEpISKjUxkePHmkH/6ZNm1JZWZkSX3FxsXbA7dWrF2VnZ2vLDh06RHXq1CGLxWJ424p832/dukU+Pj4EgJYsWaK1MTMzk6KjowkAvfXWW0qda9eupb59+9Ls2bPp22+/penTpxs6MBp1Xr9+nRo1akSLFy+me/fuaY+XlpbSypUrtQPVypUrlbZz27ZttHfvXsrPz6/0+OXLlykmJoYA0IABA5Q6qzJmzBjy8PCggQMH6trGIk5bYT5nzhzd70vGR0T0hz/8gQBQjx497Irh7Oxs2r9/v3KnM2bPnk0AKD4+Xqlz8uTJBICio6Pp6tWr2uN5eXk0ZcoUAkB169alZ8+eKXP27t2bAFCbNm3o1q1b2uNpaWkUGhpKAGjDhg0OXyt6/JfJIFGnTAaJOGUySLSNMvmjqpYzkj+iTpn8kWmnaAaZUSe7yiBRp2gGifpk8odIrC8gWwOJOGVrIKNOFTWQSDtlayAVfTujNZCIUyaDRNsoUwOp7jPrqYFEnLI1kIhTNoNs6B0bkM0fEads/hh1qsgfkXbK5k918OCkQmwDXGYNTj579oyKi4udLu/fvz8BoCFDhpjid8TXX39NACg0NNSU9a9Zs4YA0OLFi03dvu4YnFy9ejUBoN69exs6o2EGubm55O/vr3xw+8WLF+Th4UEA6OLFi3bLc3JytOC8cuWKEud//vMfAsrP5mVmZtot/+CDDwgAJSYmSnn07I/Tpk0jAPT666/bLUtJSSGg/GrSBw8eKHNWxVbYiZ61c+UsKCio9szj1KlTCQDFxMQoc7ri9OnTBJSfmS0oKDDFeeDAAQJASUlJwttYj1OmMBfx7d27lwBQq1at7IoOs5zOKCsro8jISEOFlV5nSEgIAaCdO3faLSsuLqbg4GACQHv27FHiTEtL0wYNTpw4Ybf8q6++IqD8Cg1HiB7/ZTJIVc1h5Psh4pTJIDPqKlf5o8JpNH9EnTL5I+qUySDVn6eeDBJ1imaQiE82f1zhrC+gugbS46yKbA3kymlGDeTK6QqRGsioU0UNpMepsgbS41NdA+lxOkO0BtLjVF0DuXKqzCC9YwMq80d0PELmu6HHqTp/VIy7yOYPzzn5K6Ju3brw9PR0urxv374AgJ9++qmm3hJatWoFAMjPz1e+7p9//hnvvvsurFYrZs6cqXz97ubjjz8GAMyfPx8eHu79Ku7YsQMvXrxAgwYN0K9fP2XrffnyJcrKygAAzZo1s1tev359BAUFAQBKSkqUOFNSUgAAcXFxiIiIsFs+fPhwAMDhw4eRnZ2txOkIItLm+pg8ebLd8q5du6JVq1YoLi5WOrdoTePr64vatWs7Xf76668DcE8ulZaWoqioSPn6CwsLkZSUhFdffRWLFi1Svn53smLFCgDl8834+fm59b0cPXoUmZmZ8PLywujRo5Wuu6CgAIDjXPL09NSyQ3UuhYeHo0uXLnbLhw0bBg8PD2RkZODcuXN2y0WO/7IZ5I6aQ8Qpk0FmtNFV/sg6RfLn1/JZAnIZpLqdejJI1CmaQSI+2fxxhaO+gNk1kJn9DyNOs2sgkXbK1kCunGbUQDX9eTrzmVkDGW2jihrImdPMGsiRU1UG6R0bUJk/7hiP0OtUmT+q2imbPzw4+T+EbWLXmuxQnjhxAgDQsWNH5eueOXMmcnJysGrVKnh5eSlfvyO2bduGoUOHIjExEaNHj8bKlSvx7Nkz5Z7r168jPT0dQUFB6Nq1K7799luMHTsWffr0wejRo/H555+bMqjijI0bNwIARo8eXW3Ra5R69eqhcePGAIDjx4/bLb927RoeP36MevXqoUWLFkqcT548AQCEhYU5XG57vKysDGfOnFHidERWVhbu378PAEhISHD4HNvjp06dMu19uBt35lKzZs0QGBiofP0LFizAjRs3sHTpUtSrV0/5+h2RnJyMkSNHIjExESNGjMCSJUvw4MEDpY6CggIcOnQIFosFAwcOxOHDhzF58mT06dMHw4cPx4oVK+QmuTaILZf69euH4OBgpeuOiYkB4DiXcnJykJ6eDk9PT7Rv316Jz1UueXt7a208efKk4fU7+p6ZnUHu+G6LOGXep8hrZfPHldOM/HHlNCN/HDnNziCjn6eKDHLmNCuDHPnMzh9HfQGz88fM/odKp2xOijhlM8iV04wMcuVUnUGOfGbnj9HPUkX+OHOaWQM5cqrKIL1jAyrzxx3jEaqcRvJHlVM2f9SNQjBuhYiwdetWAM6/hKooLS3F/fv3sXPnTrz33nvw9/fHP//5T6WOQ4cOYdOmTRg7dix69uypdN3VsXv37kr/v2XLFsyZMwebN29WekWh7axQq1atMG7cOGzatMnOu2zZMuzbt8/h1X8quX//Pg4dOgQAGDdunPL1L1iwABMmTMCkSZOwYsUK9OrVC56enjh58iTefvttWCwWLFmyBL6+vkp8tiC8e/euw+UVH7927RoGDhyoxFuV69evAwB8fHwQGhrq8Dm2M5a25/4v8vXXXwMwP5eICA8fPsShQ4fwzjvvwNPTE8uXL1fuuXr1KpYuXYru3btj/PjxytfvjCNHjlT6/+3bt2Pu3LlYtWoVJk6cqMRx4cIFlJSUICwsDIsXL8YHH3xQafmOHTuwdOlS7N69W9mgnTOKioq0Y5oZuTR37lz069dP21cGDBiAgIAApKWl4Z133sGLFy/w97//XTu5IourXHr58iUePXoEoDyXjODs+G9mBtVkzSHrFM0gIz5V+ePKaUb+6Gmn6vxx5jQzg4zuPyoyqDqnGRnkzGdG/rjqC5iRPzXR/1DtFMkfEadsBul1qswgI+1UkUGufGbkj+j+I5M/epyq88eVU0UGGRkbUJU/7hiPUOnUmz+yTqV9MMM/BGecYvack9Xx6aefEgDy9vamGzdumOJwdAeqoUOH0o8//qjUU1BQQFFRURQYGFhpHggzt++8efNo0aJFdOHCBXr+/Dnl5ubSd999R507d9bmLzxz5owy3/Lly7W5LgDQH//4R8rMzKTCwkI6ePCgdufs2NhY0+ejXLp0KQHlkyKbxfbt26l9+/Z2+09MTAzt3btXqcs2b4mvry9lZWXZLV+8eLHml7m7vav90TbnSsOGDZ2uY9WqVQSAXnvtNSVOR5g952R17N+/X9vWRm5AZMT5zTff2O1XPXv2pGPHjhl6r3qctomiPT09K+WemXNOfvLJJ/S3v/2Nzpw5Q48fP6b8/HxKSUnR5hWzWCwO5wwS8e3YsaNSLg0ePJjS09OpqKiITp8+TR07diQA1KRJE8rNzVXWRkds3bqVAFBgYKDhOWv0Og8fPkzdu3e3238iIyNp48aNSp0nT57U1n/q1Cm75Vu2bNGWjx071pDb2fHfjAxy5XSEqjnfROoc0QzS61OVP3qcqvNHj1Nl/uhxmpFBetrpCJkM0utUmUHV+VTmj96+gMr8kel/iH4/VPR5jOaPiFM2g4w4VWWQEaeKDNLrU5k/svuPSP4YdarIH71O2QwyOjagIn9UjEcY/W6oHAPRmz8yTtU1EBHPOfk/wfnz5zFjxgwA5VepNW/e3BRPWFgYEhISEB8fj4YNGwIov8z+yy+/RGlpqTKP7ecCCxcu1Dxm849//AOzZs1CTEwM6tSpg4CAAPTt2xdHjhxBfHw8ioqK8O677yrzvXjxAgBQXFyM7t27Y82aNYiIiICPjw/69OmDHTt2wGKx4Ny5c3ZXc6rG9rMBM65OAgAiQkZGBh4/foxatWohKioKVqsV3t7euHTpEj777DPk5OQo873xxhsIDQ1FYWEh3nzzTe2yfqD8ytiFCxdq/2+bc8UMbJfSe3t7O32Oj4+P6e/DXWRlZWHMmDEAgGnTpqFHjx6meF555RUkJCSgS5cuCAsLg8ViwenTp7Fhwwbl23Xt2rU4evQo3n77bbz22mtK1+2MP/3pT1i4cCE6deqEoKAg+Pn5oWvXrti9ezeGDRsGIsLMmTNBRNKuirnUrFkzbN++HdHR0fD29kZcXBx2796N2rVrIysrC+vWrZP2VYctl0aOHKnsquqq3Lp1C9nZ2bBYLIiIiEDbtm3h5+eHzMxMfP7558jMzFTm6ty5M2JjYwEAEydOrDT/z6lTpyrN72Nkv63u+G9WBtVUzSHrlMkgvT6V+ePKaUb+uHKakT/VOc3KIJH9RzaD9DhVZlB1PpX5o7cvoDJ/aqr/odIpkj8iTtkMMuJUlUFGnCoySK9PZf7I7j8i+WPUqSJ/9DplM8jo2ICK/HHHeIQqp5H8kXGa0gcTHtZk7HDHlZMZGRnUqFEjAkBvvvkmlZWV1Zj75MmT1K5dOwJAU6dOVbLOK1eukLe3N3Xs2NHuikF3XZlqO/Pg4eFBOTk5StZpu1oRAG3bts3hcxITEwkATZs2TYnTERcvXtTOPDq6s7UKpkyZQgCoa9eudOvWLe3xhw8f0qBBgwgov4KypKREmfPo0aNUp04dAkCenp7Upk0bCg0N1c549ujRgwC5u//xlZPOefz4MVmtVgJAvXr1osLCQtOdNq5cuUK9e/cmANSvXz9lzuzsbAoKCqLw8HC7M+ZmXjlZHdeuXdNyJC0tTdpnO1MPgD788EOHz5k0aRIBoAEDBuh6jyJtfPTokXblwpEjR3S/zohz0aJFBJTfkfPChQva47m5uTR58mQCyu80+fTpU2XO9PR07Q6ZHh4e1LJlS+1OnPXq1aPBgwcTAJowYYIup6vjvxkZJFJzyGaQiFMmg2TqKtH8ceU0I39k60ej+aPHaUYGibRTNoP0OFVmkB6f6vyx4awvYOaV20b6H6qu3DbilK2BRJw2RDPIldOMDHLldIVIBlXnMyN/XDkdIZs/epyqayA9TtEMEhkbkM0fVeMRRr4bqpxG8kf1uIts/hDxlZO/ah48eIC+ffvi/v37GDhwINavXw+LxVJj/s6dO2PPnj3w8fHBZ599hv/+97/S65w2bRpKSkqwevVqt9/B2sZvfvMbAOU3UMnIyFCyzvr162v/tt3VqiqtW7cGAKVX71Tl3//+NwCgR48epsxteeHCBaxZswZeXl746quvEBkZqS179dVXsWnTJgQHB+PixYvavBgq6NatG86fP49JkyYhJCREO0M3depUnD17VjubFxISosxZFdtn/PTpU6dndG0TRFfcH37t5OXlYcCAAbhy5QpiY2Oxc+dO7exkTdC6dWvs2rULDRs2xL59+3Ds2DEl6/3rX/+KnJwcfPTRRwgICFCyTllatmyp3e3+xo0b0uv7peTSli1bUFxcjMjISHTr1k35+rOzszFv3jwAwPr167WJ4QEgICAAn3zyCaxWK+7du4dVq1Yp80ZHRyM1NRUzZsxAZGQkMjMz8eLFC4wZMwbnz59H3bp1AejLJT3Hf9UZ5I6aQ8Qpk0GybRTJHz1O1fmj4rM0mj9G9llATQaJtlMmg/Q4VWaQ3jaqzJ+KOOsLmFkDmdH/UOVUWQOJtFO2BnLmNLMGEv08RWsgV/ssoL4GMtJGVTWQM6eZNVB17RTNIJGxAdn8ccd4hAqn0fxR3U4VfbBfxugPY5icnBz07dsXN2/eRM+ePbF169Yau4NURUJDQ9G+fXuUlZXhwoUL0utLTU2FxWLBkCFDEBISUulvy5YtAIDFixcjJCQEcXFx0j49VNyuJSUlStYZHR2t/dtZaNgeN+snK2VlZfjyyy8BmPeT7pSUFBARWrZs6XBS5bp16yI+Ph4AcPbsWaXuqKgorF27Frdv38bLly9x9+5drF69GvXr19f2VdtPDMzAdvfxoqIi3Lt3z+FzbIPdqu5U7m6Kiorwxhtv4NSpU7Bardi3bx/q1KlT4+/D398fvXr1AlD+UzcVpKamAgD+/Oc/22XThx9+CADYvHmz9lhNYcsnFdn0S8gl4P9/zjR27FhTBr/Onj2LwsJCBAQEaPlTEU9PT23/UZ1LISEhWLFiBW7evImioiJkZ2dj48aNaNq0qeZylUt6j/8qM8gdNYeIUyaDVLXRSP7odarMH5Wfpd780etUmUEy7RTNIL1OVRlktI0q8scRjvoCZtdAqvsfKpxm1EAi7ZStgRw5za6BRD9P0RrIkc/sGkhvG1XWQI6cZtdA1bVTJINExgZk88cd4xGyTpH8MaOdsvnDd+v+FWIbFb906RLi4uKwa9cuXbeINwvbAUHVwF1paSkePnzodHleXh7y8vJMm4esKpcvX9b+HR4ermSdHTp0gK+vLwoLC5GRkYGoqCi759hCMywsTImzKsnJybhz5w58fX0xYsQIUxy5ubkun2M7o2WbH8Rs9u/fj7y8PISGhqJjx46meZo0aYKQkBA8ePAAKSkpGDVqlN1zUlJSAJSfafy1U1JSglGjRuH7779Hs2bNcODAAQQHB7v1/VT8ryqqy6aCgoIanT/00aNHyM7OBqAmm8LDw9G4cWPcvn3b6VXiZufSzZs3ceLECQDlhbkZ/BJz6fLly7h27Rp8fX3x29/+1unzjBz/VWWQO2oOEadMBqluo578EXHK5o/KdurNHyNOVRkk007RDDLiVJFBKj9LvflTHVX3+Zqogcw6zos4zayBRNopu22cvd7MGsjoe5atgar6aqIGctVGM2qgqs6aqIGMfpauMsjo2ICK/HHHeISoUyZ/zGinTP7wlZO/MiqOirdp08ZtVybZyMzM1M6KtGvXTnp9tsuvHf1NmDABADB//nwQkak/K6zIsmXLAJRf4q+qQ+7v748BAwYAAL744gu75Q8ePMD+/fsBAImJiUqcVbH9pHvIkCEIDAw0xWE7G/XTTz/h9u3bdsufP3+OM2fOACj/eYbZvHz5Eu+//z4AICkpCbVq1TLNZbFYMGzYMADlk4hX5fjx40hPT4eXlxeGDBli2vuoCYgIEydOxM6dOxEaGoqDBw8iNDTUbe/n2bNnSE5OBgC0b99eyTrT0tKcZtOcOXMAAJMnT9YeqwmWL18OIkJgYKCyM7cjR44EAGzYsMFuWWFhoXYm1excio+Pr3QVg0psuZSXl4fTp0/bLS8pKcEPP/wAoGZyiYgwa9YsAMCYMWOc/sTR6PFfRQa5o+YQccpkkOo26skfo04V+aO6nXryR8Qpm0Gy7RTJIKNO2QxS+VnqzZ/qcNQXMLsGUt3/kHGaWQOJtFO2BnLkNLsGEmmnTA3kzGdmDaSnjaprIEdOs2sgo5+lqwwSGRuQzR93jEeIOmXyx4x2SvfBRCaqZBxj9g1bSkpKaOjQoQSAmjdvTvfu3TPFU5GzZ8/S+++/Tzdv3rRbtnfvXmrVqpXQxMAimLV9v/vuO3rvvfcoIyOj0uNPnz6l6dOna5Mjb968Wak3LS2NatWqRR4eHrR+/Xrt8SdPntDvfvc7AkDNmjWjoqIipV4iovz8fO2GMbt27VK+fhu5ubkUHBxMcHFDHF9fX7pz544y7+7du+nkyZOVHsvKyqL+/fsTALJarcITlNvQsz9mZGSQt7c3AaAlS5Zok9RnZmZSdHQ0AaCkpCSlzqrUxA1xbN+T4OBgunLlipDHiPPu3bs0Y8YMunTpkt2yEydOUJcuXQgAtW3bVveNlmTyxawb4ly6dImSkpLs2llQUEALFy4kDw8PAkCLFi1S4iMiun//PgUEBBAAWrBggTZBdn5+vvb6+vXrU3Z2tjJnRaKioggArVy5UtfzRZxlZWXaZOFVJ4N//vy5Nhk8ADp79qwSJ1H5jboOHjxY6WYVjx490l7bsGFD+vnnnx2+VvT4L5NBqmoOI98PUadoBon4ZPNHdS2nZ/uKOGXzR7SdMhmkYtsazSARp0wGibZRJn9E+wIy+aOq/2Ekf0SdMjWQiFM2g1T37fRsYxGnTAaJtlEmf1RsV6P5I+KUrYFE2ymTQc6orvZS3QfT46yKqhtyVedU3Qdz5TSjD1YRHpyU4NixY/TKK69ofz4+PgSAateuXenxrKwsJb7NmzdrYdGiRQtKSEhw+DdixAglPiKi5ORkzRkSEkKdOnWimJgYqlevnvZ4XFyc4TARwazByW+++UZrS1hYGMXFxVH79u21QLNYLFJ3da6O1atXk8ViIaD8LtKdOnWi2rVrayGTmppqite2LzVo0ICKi4tNcdjYs2cP+fr6EgCqVasWtWjRgqxWq7Z9PT09Kw3OqmDGjBlaAdGhQwdq3bq1tp2tVqvQQKjo9/2LL77QCqiwsDDq0KGDdhe+2NhYysvLU+rMysqqtMzPz48AkI+PT6XHjx07psR5/Phx7fvTuHFjp7mUkJCgrJ23bt3SnEFBQdSxY0fq0KGDNhBu67TduHFD+efpCL3Fh1Fnamqq1p4GDRpQbGwsxcbGahlhczq7+6xoG3fu3Kl9Pxs2bEhxcXEUGBiovXb//v2mbFfbvuTl5WXomCLiPHfuHNWvX1/L+MjISIqJidG+L7aOiUrnRx99RACoTp06FBMTQ23btiVPT08tG3788UenPpnjv2gGiTplMkjEKZNBIj7Z/FFdy+nJHxGnbP7ItFM0g2S3rUgGiTpFM0jUJ5M/Mn0B0fwRdcrkj4hTtgYSccpmkOq+nZ4MEnHKZJBMG0XzR3a7iuSPqFOmBhJ1ymSQM1yNDcj0wUScsn0wo04VfTCjThV9sOrgwUkJKn45q/ureLWYDOvWrdPli4iIUOIjIsrJyaGPP/6YhgwZQs2bN6eAgADy9vamRo0aUf/+/WndunWmD27ZMGtwMisri2bPnk2JiYnUpEkT8vPzI19fX2ratCmNHz/e7go81Rw5coQGDx5MwcHB5O3tTZGRkfTWW28pvZKwKrYrCKdPn26aoyLXrl2jKVOmUFRUFPn4+JC3tzdFRETQuHHj6Ny5c8p9ycnJNGrUKIqMjCRfX18KDAyk+Ph4WrZsmfAVkzLf95SUFBo0aBAFBQWRj48PRUdH09y5c6mgoEC5s+JBo7q/5ORkJU69zwecH26MOgsKCujTTz+lUaNGUcuWLSkwMJA8PT2pQYMGlJiYSP/6178oPz9f+bZ1ht7BSaPOJ0+e0Pz586l///7UtGlTLX/Dw8NpxIgRtG/fPtPaePHiRRo9ejSFhISQl5cXhYaG0vjx4+nq1aumOZOSkggADR48uFqHKufdu3fpL3/5C1mtVvLz89PaOXz4cPr++++VO1NTU2n8+PHUokUL8vf3p4CAAIqJiaG5c+fSs2fPqvXJHv9FMkjUKZNBIk6ZDBLxyeaP6lpOT/6IOGXzR7adIhkk6xTJIBmnSAaJ+mTyR7YvIJI/ok6Z/BFxytZAIk7ZDFLdt9OTQSJOmQySbaNI/sg6RfJHxilaA4k6ZTLIGXrGBkT7YCJO2T6YUaeKPphRp4o+WHVYiIjAMAzDMAzDMAzDMAzDMAxTw/ANcRiGYRiGYRiGYRiGYRiGcQs8OMkwDMMwDMMwDMMwDMMwjFvgwUmGYRiGYRiGYRiGYRiGYdwCD04yDMMwDMMwDMMwDMMwDOMWeHCSYRiGYRiGYRiGYRiGYRi3wIOTDMMwDMMwDMMwDMMwDMO4BR6cZBiGYRiGYRiGYRiGYRjGLfDgJMMwDMMwDMMwDMMwDMMwboEHJxmGYRiGYRiGYRiGYRiGcQs8OMkwDMMwDMMwDMMwDMMwjFvgwUmGYRiGYRiGYRiGYRiGYdwCD04yDMMwDMMwDMMwDMMwDOMWeHCSYRiGYRiGYRiGYRiGYRi38H8ymNiUVyk9wQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1600x900 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot results\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(16, 9))\n",
    "# boxplot plot for the benchmarks. Make the median line red\n",
    "bp = ax.boxplot(\n",
    "    assortativity_random.T,\n",
    "    positions=np.arange(len(graphs)) + 1,\n",
    "    showfliers=False,\n",
    ")\n",
    "# Scatter plot of assortativity\n",
    "ax.scatter(\n",
    "    np.arange(len(graphs)) + 1,\n",
    "    s=100,\n",
    "    marker=\"+\",\n",
    ")\n",
    "plt.show()\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "    5. BSS targeted villages leaders as seed nodes for their diffusion process. Is this a good strategy?\n",
    "    Simulate the diffusion process using the BSS seed nodes and compare it with a random seed nodes strategy.\n",
    "    Can you identify better seed nodes?\n",
    "\n",
    "    You can use the function below to simulate the diffusion process.\n",
    "\"\"\"\n",
    "\n",
    "qN = 0.05\n",
    "qP = 0.35\n",
    "\n",
    "\n",
    "def diffusion_InformationModel(village_num, qN, qP, A, seeds, T, leadapt=leadapt):\n",
    "    \"\"\"\n",
    "    Diffusion information model with endorsement. Simpler version of the model in the paper.\n",
    "    qN: probability of non-infected to transmit the innovation\n",
    "    qP: probability of infected to transmit the innovation\n",
    "    A: adjacency matrix\n",
    "    leaders: boolean array with leaders\n",
    "    T: number of periods\n",
    "    \"\"\"\n",
    "    N = A.shape[0]  # Number of individuals.\n",
    "    \n",
    "    infected = np.zeros(N, dtype=bool)  # Nobody has been infected yet.\n",
    "    infectedbefore = np.zeros(N, dtype=bool)  # Nobody has been infected yet.\n",
    "    contagiousbefore = np.zeros(N, dtype=bool)  # People who were contagious before\n",
    "    contagious = np.array(seeds).astype(bool)  # Newly informed/contagious.\n",
    "    dynamicInfection = []  # Will be a vector that tracks the infection rate for the number of periods it takes place\n",
    "\n",
    "    contagious_at_t = np.zeros((T, N), dtype=bool)\n",
    "    infected_at_t = np.zeros((T, N), dtype=bool)\n",
    "    \n",
    "    p = leadapt.loc[leadapt.village == village_num, \"adopter\"].mean()\n",
    "    beta = np.log(p / (1 - p))\n",
    "    \n",
    "    for t in range(T):\n",
    "\n",
    "        contagious_at_t[t, :] = contagious.copy()\n",
    "        infected_at_t[t, :] = infected.copy()\n",
    "\n",
    "        # Step 1: Take-up decision based on Newly informed, p. \n",
    "        LOGITprob = 1. / (1 + np.exp(-beta))\n",
    "        infected = ((~contagiousbefore & contagious & (np.random.rand(N) < LOGITprob)) | infected)\n",
    "        infectedbefore = (infectedbefore | infected)\n",
    "        contagiousbefore = (contagious | contagiousbefore)\n",
    "        \n",
    "        # Step 2: Information flows\n",
    "        transmitPROB = (contagious & infected)*qP + (contagious & ~infected)*qN  # Probability of transmission: Pr(.|infected)==qN, Pr(.|~infected)==qP. Individual (node) specific.\n",
    "        contagionlikelihood = transmitPROB.reshape(-1, 1)*A\n",
    "        newcontagious = ((contagionlikelihood > (1. -np.random.random(size=contagionlikelihood.shape))).astype(bool)).any(axis=0)\n",
    "        \n",
    "        # Step 3:\n",
    "        contagious = newcontagious | contagiousbefore\n",
    "        dynamicInfection.append(np.sum(infectedbefore)/N)\n",
    "    \n",
    "    return contagious_at_t, infected_at_t"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "yfinance",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.0"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}