{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# NTDS'18 milestone 1: network collection and properties\n",
    "[Effrosyni Simou](https://lts4.epfl.ch/simou), [EPFL LTS4](https://lts4.epfl.ch)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Students\n",
    "\n",
    "* Team: `42`\n",
    "* Students: `Alexandre Poussard, Robin Leurent, Vincent Coriou, Pierre Fouché`\n",
    "* Dataset: [`Flight routes`](https://openflights.org/data.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rules\n",
    "\n",
    "* Milestones have to be completed by teams. No collaboration between teams is allowed.\n",
    "* Textual answers shall be short. Typically one to three sentences.\n",
    "* Code has to be clean.\n",
    "* You cannot import any other library than we imported.\n",
    "* When submitting, the notebook is executed and the results are stored. I.e., if you open the notebook again it should show numerical results and plots. We won't be able to execute your notebooks.\n",
    "* The notebook is re-executed from a blank state before submission. That is to be sure it is reproducible. You can click \"Kernel\" then \"Restart & Run All\" in Jupyter."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Objective "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The purpose of this milestone is to start getting acquainted to the network that you will use for this class. In the first part of the milestone you will import your data using [Pandas](http://pandas.pydata.org) and you will create the adjacency matrix using [Numpy](http://www.numpy.org). This part is project specific. In the second part you will have to compute some basic properties of your network. **For the computation of the properties you are only allowed to use the packages that have been imported in the cell below.** You are not allowed to use any graph-specific toolboxes for this milestone (such as networkx and PyGSP). Furthermore, the aim is not to blindly compute the network properties, but to also start to think about what kind of network you will be working with this semester. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 1 - Import your data and manipulate them. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  A. Load your data in a Panda dataframe."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, you should define and understand what are your nodes, what features you have and what are your labels. Please provide below a Panda dataframe where each row corresponds to a node with its features and labels. For example, in the the case of the Free Music Archive (FMA) Project, each row of the dataframe would be of the following form:\n",
    "\n",
    "\n",
    "| Track   |  Feature 1  | Feature 2 | . . . | Feature 518|  Label 1 |  Label 2 |. . .|Label 16|\n",
    "|:-------:|:-----------:|:---------:|:-----:|:----------:|:--------:|:--------:|:---:|:------:|\n",
    "|         |             |           |       |            |          |          |     |        |\n",
    "\n",
    "It is possible that in some of the projects either the features or the labels are not available. This is OK, in that case just make sure that you create a dataframe where each of the rows corresponds to a node and its associated features or labels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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>node_idx</th>\n",
       "      <th>AirportID</th>\n",
       "      <th>Airline</th>\n",
       "      <th>Stops</th>\n",
       "      <th>DestAirportID</th>\n",
       "      <th>SourceAirportID</th>\n",
       "      <th>DestSourceRatio</th>\n",
       "      <th>Country</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>2.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>Papua New Guinea</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>2.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>8.0</td>\n",
       "      <td>8.0</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>Papua New Guinea</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "      <td>2.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>12.0</td>\n",
       "      <td>0.833333</td>\n",
       "      <td>Papua New Guinea</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3</td>\n",
       "      <td>4</td>\n",
       "      <td>2.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>11.0</td>\n",
       "      <td>11.0</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>Papua New Guinea</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4</td>\n",
       "      <td>5</td>\n",
       "      <td>4.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>51.0</td>\n",
       "      <td>49.0</td>\n",
       "      <td>1.040816</td>\n",
       "      <td>Papua New Guinea</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   node_idx  AirportID  Airline  Stops  DestAirportID  SourceAirportID  \\\n",
       "0         0          1      2.0    0.0            5.0              5.0   \n",
       "1         1          2      2.0    0.0            8.0              8.0   \n",
       "2         2          3      2.0    0.0           10.0             12.0   \n",
       "3         3          4      2.0    0.0           11.0             11.0   \n",
       "4         4          5      4.0    0.0           51.0             49.0   \n",
       "\n",
       "   DestSourceRatio           Country  \n",
       "0         1.000000  Papua New Guinea  \n",
       "1         1.000000  Papua New Guinea  \n",
       "2         0.833333  Papua New Guinea  \n",
       "3         1.000000  Papua New Guinea  \n",
       "4         1.040816  Papua New Guinea  "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "routesCol = ['Airline','AirlineID','SourceAirport','SourceAirportID',\n",
    "             'DestAirport','DestAirportID','CodeShare','Stops','Equipment']\n",
    "\n",
    "airportsCol = ['AirportID', 'Name', 'City', 'Country', 'IATA', 'ICAO', 'Latitude', 'Longitude',\n",
    "               'Altitude', 'Timezone', 'DST', 'TzDatabaseTimezone', 'Type', 'Source']\n",
    "\n",
    "routes = pd.read_csv(\"data/routes.dat\",header = None,names = routesCol,encoding = 'utf-8', na_values='\\\\N')\n",
    "airports= pd.read_csv(\"data/airports.dat\",header = None, names = airportsCol, encoding = 'utf-8')\n",
    "# We drop nan value for source and destination airport ID because they are not in the airports file\n",
    "# and it does not make sense to keep routes that go nowhere.\n",
    "routes.dropna(subset=[\"SourceAirportID\", \"DestAirportID\"], inplace=True)\n",
    "\n",
    "# Get unique airlines for a source airport.\n",
    "routesAirline = routes[['SourceAirportID','Airline']].groupby('SourceAirportID').nunique().drop(['SourceAirportID'], axis=1)\n",
    "# Get average number of stops of outbound flights (from Source)\n",
    "routesStop = routes[['SourceAirportID', 'Stops']].groupby('SourceAirportID').mean()\n",
    "# Get number of routes that leave the airport\n",
    "routesSource = routes[['SourceAirportID', 'DestAirportID']].groupby('SourceAirportID').count()\n",
    "# Get number of routes that arrive to the airport\n",
    "routesDest = routes[['SourceAirportID', 'DestAirportID']].groupby('DestAirportID').count()\n",
    "\n",
    "# Concatenate everything and fill nan with 0 (if no departure = 0 airline, 0 stop and ratio of 0)\n",
    "features = pd.concat([routesAirline, routesStop, routesSource, routesDest], axis=1, sort=True).fillna(0)\n",
    "features.index = features.index.astype('int64')\n",
    "# Create Ratio\n",
    "features['DestSourceRatio'] = features['DestAirportID']/features['SourceAirportID']\n",
    "# Add countries (as labels)\n",
    "airports_countries = airports[['AirportID', 'Country']].set_index(['AirportID'])\n",
    "features = features.join(airports_countries).sort_index()\n",
    "\n",
    "features.reset_index(level=0, inplace=True)\n",
    "features = features.rename(columns={'index':'AirportID'})\n",
    "features.reset_index(level=0, inplace=True)\n",
    "features = features.rename(columns={'index':'node_idx'})\n",
    "\n",
    "features.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For now our features are:\n",
    "* Airline : the number our unique airline leaving an airport.\n",
    "* Stops : the average number of stops of an airport's outgoing flights.\n",
    "* DestAirportID: number of the airport incoming flights.\n",
    "* SourceAirportID: number of the airport outgoing flights.\n",
    "* DestSourceRatio: Ratio of the incoming fligths over the outgoing flights. (can be infinite)\n",
    "\n",
    "The only label we have for now is the country.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### B. Create the adjacency matrix of your network."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember that there are edges connecting the attributed nodes that you organized in the dataframe above. The connectivity of the network is captured by the adjacency matrix $W$. If $N$ is the number of nodes, the adjacency matrix is an $N \\times N$ matrix where the value of $W(i,j)$ is the weight of the edge connecting node $i$ to node $j$.  \n",
    "\n",
    "There are two possible scenarios for your adjacency matrix construction, as you already learned in the tutorial by Benjamin:\n",
    "\n",
    "1) The edges are given to you explicitly. In this case you should simply load the file containing the edge information and parse it in order to create your adjacency matrix. See how to do that in the  [graph from edge list]() demo.\n",
    "\n",
    "2) The edges are not given to you. In that case you will have to create a feature graph. In order to do that you will have to chose a distance that will quantify how similar two nodes are based on the values in their corresponding feature vectors. In the [graph from features]() demo Benjamin showed you how to build feature graphs when using Euclidean distances between feature vectors. Be curious and explore other distances as well! For instance, in the case of high-dimensional feature vectors, you might want to consider using the cosine distance. Once you compute the distances between your nodes you will have a fully connected network. Do not forget to sparsify by keeping the most important edges in your network.\n",
    "\n",
    "Follow the appropriate steps for the construction of the adjacency matrix of your network and provide it in the Numpy array ``adjacency`` below: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The number of nodes in the network is 3330\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[0, 1, 1, ..., 0, 0, 0],\n",
       "       [1, 0, 1, ..., 0, 0, 0],\n",
       "       [1, 1, 0, ..., 0, 0, 0],\n",
       "       ...,\n",
       "       [0, 0, 0, ..., 0, 0, 0],\n",
       "       [0, 0, 0, ..., 0, 0, 0],\n",
       "       [0, 0, 0, ..., 0, 0, 0]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "edges = routes[['SourceAirportID', 'DestAirportID']]\n",
    "edges = edges.astype('int64')\n",
    "airportID2idx = features[['node_idx', 'AirportID']]\n",
    "airportID2idx = airportID2idx.set_index('AirportID')\n",
    "edges = edges.join(airportID2idx, on='SourceAirportID')\n",
    "edges = edges.join(airportID2idx, on='DestAirportID', rsuffix='_dest')\n",
    "edges = edges.drop(columns=['SourceAirportID','DestAirportID'])\n",
    "\n",
    "n_nodes = len(features)\n",
    "adjacency = np.zeros((n_nodes, n_nodes), dtype=int)\n",
    "\n",
    "# The weights of the adjacency matrix are the sum of the outgoing flights\n",
    "for idx, row in edges.iterrows():\n",
    "    i, j = int(row.node_idx), int(row.node_idx_dest)\n",
    "    adjacency[i, j] += 1\n",
    "\n",
    "print(\"The number of nodes in the network is {}\".format(n_nodes))\n",
    "adjacency"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Execute the cell below to plot the (weighted) adjacency matrix of your network."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1.05,'adjacency matrix')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQsAAAEGCAYAAAByy7CMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzsnXl8lNW9/99nJpPJZLIvhJAACSSRRRABZa0roqC2tsW2v9taim21i73ayr1VW69yvS7XqxVvbb21Fkrt4kKvrQoo4gKyCxjWQMhGFpKQZLJMZiYzk5nn98ez5JnJbFnA2JvP6zWvmXmWc85znuf5nnO+y+crJEliFKMYxSiiwfBpN2AUoxjFZwOjwmIUoxhFTBgVFqMYxShiwqiwGMUoRhETRoXFKEYxipgwKixGMYpRxIRRYfEZgxCiQAghCSHilP9bhBArP+12fRYhhDguhLjq027HZwVi1M/iswUhRAFQDZgkSer9dFszMiGE+D1QL0nSzz/ttvwjYXRmMYr/c1BnZaMYGEaFxQiAEOI+IUSlEMIuhDghhPiibp9RCPGUEKJVCFEF3Bh07odCiO8ovycLId4XQrQpx/9JCJGmO3a8EOJ/hRAtyjHP6fbdLoQoE0K0CyHeEUJM1O2ThBDfE0KcVvb/SgghdPu/q5yrtn+2EOJfhBB/DWrrL4UQa8P0QY1yzhEhhEMI8TshRI6yzLILIbYJIdJ1x78mhGgSQnQKIXYIIaYr2+8Avg78qxCiWwjxpq78nwohjgAOIUScsm2Jsn+zEOJpXfmvCCHWxXYH/49AkqTRz6f8AW4FxiEL768CDiBX2fc94CQwHsgAPgAkIE7Z/yHwHeV3EXAdYAaygR3AWmWfETgMPANYgQRgsbLvFqACmArEAT8HduvaJwFvAWnABKAFuEHX9gbgMkAobZgI5CrXkaYcFwecA+aE6YMaYC+QA+Qpxx4CLlWu533gId3xtwPJyr61QKlu3++B/whRfqnSjxbdtiXK77FKndcgC5sqIPnTfjZG0udTb8DoJ8RNkR/qLyi/3we+p9u3NJywCFHOLcAnyu8FykseF+K4LcC3df8NgBOYqPyXVMGi/H8VuE/5/Q5wd5j6twDfVX7fBJyIcM01wNd1//8KPK/7/yPgb2HOTVPamKr8Dycsbg+xbYnu/5eAOqBVf72jH/kzugwZARBCfFMIUSqE6BBCdAAXA1nK7nHID7CKMxHKGSOEeFkI0SCE6AL+qCtnPHBGCq0UnQg8q6vfhjxLyNMd06T77QSSdOVWhmnSBuAbyu9vAC+Fa7uCZt1vV4j/SaAtzZ5Qlm5dyC899F1rONRF2f8W8gzslCRJO6Mc+38Oo8LiU4aiG/gtcBeQKUlSGnAM+WUFaER+IVVMiFDc48gj7ExJklKQX1C1nDpgQhjlXh1wpyRJabqPRZKk3TFcQh0wOcy+vwEzhRAXI88s/hRDebHgn4AvAEuAVKBA2a5eazgTXzTT36NAGZArhPh/Q2zjPxxGhcWnDyvyQ9wCIIRYhTyzUPEq8M9CiHxFwXdfhLKSgW6gQwiRB/yLbt9+ZMHzhBDCKoRIEEIsUvb9D3C/TkmYKoS4Ncb2vwisFkLMETKKVOWoJEk9wEbgz8B+SZJqYywzGpIBN9AGJAKPBe1vBiYNpEAhxBXAKuCbyueXSh+OQsGosPiUIUnSCeBpYA/yQz4D2KU75LfIeoHDyAq//41Q3BpgNtAJbNIfK0mSD7gZWQFZC9QjK1ORJOl14D+Bl5Vp/TFgWYztfw15RP4zYEeeTWToDtmgXFO0JchA8Afk5VgDcAJZMarH74BpyrLqb9EKE0KkKGXeJUlSg7IE+R2wXm/1+T+PT1tpEu4D3ACcQtbS33ee66oBjiIrFg8o2zKAd4HTyne6sl0A/6206wgwexD1rUPWvB/TbRtwfcBK5LV8E7ByGNrwMPILWKp8luv23a+04RRwfaz3CXnZ5ARSgraPR7bslAHHUZSkQ+iH08on5n6I0IZh74cIbUhAnvUdVtqwRtleCOxTrukVIF7Zblb+Vyj7C6K1bdjek/P5Eg66UbKSqRJ5KhmvdOS081hfDZAVtO1J+jT+9wH/qfxejqzlF8B8YN8g6rsCeQZwbLD1KS9VtfJQL0M29aUPsQ0PA6tDHDtNuQdm5SGuVO5RxPuEPHNdC6wLUWau+sIjLyvKlXoG0w9Vynf6QPohQhuGtR+itEEAScpvE7IAmI+8/Pyasv1/gO8rv38A/I/y+2vAK5HaNpzvyUhdhlwOVEiSVCVJkgd4GVmhdSHxBeQpNMr3Lbrtf5Bk7AXShBC5AylYkqQdyBaHodR3K7K14kPgbeRR+IYhtiEcvgC8LEmSW5KkauTR63Ii3CchhBXoQvb7eChE/Y2SJB1SftuRR/e8QfTD9cC7kiTZJElqH0g/RGjDsPVDDG2QJEnqVv6alI+E7O+xUdke3A9q/2wErlWWSuHaNmwYqcIij0AzVz2Rb+JQIQFbhRAHFQ9AgBxJkhpBfqiAMee5bQOtLxn4d0mSvi7JQ8twteMuxYtync5jMlwbwvaFJEkOSZKSJEmaLklSRJOlEu9yKfKoOtB+GJb7EdQGGKZ+iLFuoxCiFHlZ+C7yrKBD6jNz68vT6lL2dwKZQ21DLBipwiKUUul8RrwtkiRpNvJ0/oeKZjwcLnTbwtV3PtrxPLIZdBay5UR1fz5vbRBCJCE7YN0jSVJXpEMvYBsuaD9IkuSTJGkWkI88G5gaobwL+TwEYKQKi3oCfQvygbPnqzJJks4q3+eA15FvWLO6vFC+z53ntg20vmFvhyRJzcqD60e2wqjT2PPSBiGECfkl/ZMkSarl5oL2Q6g2XOh+UCFJUgfysnI+8jJL9YnRl6fVpexPRV5Onv93ZjgVIMP1QY4jqEJW1KgKo+nnqS4rSgyA8ns38pr3vwhUtD2p/L6RQEXb/kHWW0CgcnFA9dGn4ExXPtVAxhDbkKv7/WPkNTDAdAKVZ1XISr1B3yflev6AErvyafRDhDZcyH7Ipi9+xgJ8hOzA9hqBCs4fKL9/SKCC89VIbRvWd+V8vIDD0jBZ+12OvH772XmsZ5LSyarp6mfK9kzgPWTT1XvqA6g8YL9S2nUUmDuIOv+CPL31Io8I3x5MfcjBVBXKZ9UwtOElpY4jwBtBL83PlDacApYN9T4Bi5GnyUfQmSgvZD9EaMOF7IeZwCdKXceAf9M9l/uVa3oNMCvbE5T/Fcr+SdHaNlyfUfKbUYxiFDFhpOosRjGKUYwwjAqLUYxiFDFhVFiMYhSjiAmjwmIUoxhFTLjgwkIIcYMQ4pQQokIIESncehSjGMUIwgUVFkIII7L5axly4Mv/E0JMi3D8HeH2XSiMtuHTr3+0DSOj/gs9sxhowM2n/oAw2oaRUD+MtuFTr/9CC4sLHSA2ilGMYphwQZ2yFKq26yVJUvNc3AZcLknSj3TH3IEiQQ2WlDlxqXLQYZxBkJwQR4fT2y86ZkyymfTEeJo6e3B4ehmTbMbl9eFw+/D5JYwGOcbG4/MTZxBMyEjEao7D0+un3ekhPTGe+LjQcrOlpYXs7Gz5/KDjPb1+WrvdAGQlmUOWEUsdzV09nLO7GZNsJiclIWIb9McCEc8bDkRqm7oPwGQ04PX5tXvR7vSQZI6j292rfev7QO2XJHMcnS4vENiHatkGIfBLEkYkpuWlhWyPp9dPfbsLh0cO0rSYjBSNSdLaebShs991ZVrjSbWY+rUrEvT3Ifj6z+c9CFe/HrE8Z+Fw8ODBVkmSQhesw4XOzBQ12EWSpBeAFwDMucVS7srAnDRjwxS86tpinn3vNIlAktmI5PYR7tbNLMzglTsX8My75Tz73mlWXVvMj68ridr44OPV/2r9ocoIVUedzcm6ndUA3L64EICNB+tZMSef8RmJEdtQZ3Nqxw7kvMFCX19wHXU2J89uK6fiXDffuaKQ080OVszJZ+PBerlfTAKjVyLBYsLl8vLlhQXcvriQjQfr6XJ5Wb+7hovyUymtl19mS3I8syemMzbFwrIZY9lV0UaGNY4XdlRz//Ip3DQzL2R71D7ON8fR5e7lypIsNtw+T2vnynV72V7e1u/afEaB0Sfx5YUFPPT56YPqG/19HI57EKm/I2HNG8dZv7sm5LVEK1MIEZYxXo8LLSw+BoqFEIXIDE9fQ2ZqDgnzACTkijn5lNa2c+BMO9/+XAH1th7ePtZEt8dHcryR2QXp1NtcxJsM3Lu0RDtH/x1LHcHfXS4v9h4vXS4vdTZnv5uhHruoKJM1bxzXtq/fXaP9fujz02MSVgDjMxIDjo31vMEiuL7gfU99ZVa/7YuKMlm3swq72wdAuzJzsPd4NUEyKz+VVQsLaOpyUVrfSbwRmu0ethzrY/9XH/qVC/u4d0O1R+3jrSea6Gq00+n0BuzfcPt87XedzcnSZ7bj8vrx+CStXYPB+IzEQQmZSFD7p8vlJcViGpaBQBPeDO15uaDCQpKkXiHEXcgEtEZkurXj4Y4vyLQSieBAj/EZiUzKTmL76VY6nT7lIS5l46EGrr94bMiHOtKLEK6O4ONTLCZAfvlTLKZ++9VzVMkPsGphAbN0I+o/GrYcbdIERTC6XF6mjk3Wrj3cgDCQF1jt40VFmTy9tVwbDMIhwxpPQ0cPAjmKLDnBFHNd5xuq4Gtod7J+dw1dLm9MAmlOQRpbTyQwpyAtYHudzUmXy8uqhQUxD4rhcMH9LCRJ2ixJUokkSZMlSXo00rFen5+4IXArqw/B+XoY9BL77muLY7oZs/JTuX1xIb/8p9ncfW2xtgz5R0SCThDMyk8lOcHE+t01moAore/kiKJP8IaWLQPC5YWZvHLnAi4vzAx7zLqd1TR09ACyoLDGG/q9YJ8mVME30Gf2pT21NHT08NKewGwLGw/Ws353DScaYx12w2NEe3A2dvbQG6P+dX91G+8cayTOAKmJRupsTk40dmIygN3t4au/2cP+6v7r1kioszl55t1y6mzOgHqWrd3Bsmd30NDhwGSAT2ptdClT7TqbkzVvHGfNG8fZX93G6ldLueW5naQmGslJjic10cTqV0u59ukPWLeziq//dg+F923ill99RJ3N2a9O/f/gstXf+vYNBPur26L2S6g+CIW3jjSw6In32LC7iqYuF0nxBnp6/dr++nYn1a3dTM6y0tLtJtFkwGQAr3JM8G3+6HRLTNemb18s19PU6Qr47/D4eeB/j0WsJ9Y+GA6o/ZiaaGReYQbLZoTT0gXi3qUlzCvM6DerWlSUSV5aAvuqbWw8WB+wT70ujKb4WOoY0annE0yxy7Knt5ZztkvWTP92RzWdTh/7qtsB2HhQ1qE+tqmMv921OOYyQ631nt5aTlmTHYBTTXb8EpTWd1Fa3xWwJAE40djFvmqZE7eypRu720ezvVUr3+PzadP10rou7Wbq69S3QV/2J7Xt2lQ+1PInFjy9tZx91Tae3lrOK3cuiLkPQuHxzSdp6OjhqXfKQy5BWh1ePixvDXFmaDTbPWGXduHat7eqLeL17K9u490Tzf22d/X0sm5nddjp/nCt+WOB2o+/+6gGu9vHroq2iDMlFeqsKhi7Ktpo6OhhXmFGv5mvel3GxNRoaR8BMD788MOxXcWngN+88MLD/ouujenY+5dfxNZjTfiB6bnJJJnjONFoDzjGHGdgydQculxe1u2sJj/dQqol/HQvP90CErR2u9lT2cbk7CRyUs3sr7KRm5bANVOyOd1s56KcJPLTLHx9/gQmZCRy6Ew7Pp+Ew+3FACTEGRmTbKbT6SUx3kBSvDFg1FVx24IJXDMlh5QEE4uKMnn1QD2LijLp7vGy6Ugjp5u7EAjcPj+XTkjjc0XZXDo+ja/Pn6hdh6qhV6+tzuZk7bvlbD/VwuTsJG3bup3VpFiMHK7r5GynCwk/8ydlsb+6jZ+8epiJmYnkpSdiNECtzclXLxtPXnp4RVucAQ7WtLNibh72nl7anf11DglxEOKyw2Lq2GS+uXBiQL3B15efbiElQVYE5qSaOVzXwdcuz+eDky397u93NhzgXLcnZF1F2VaWTg89isfaB8OBnBT5Gr6xYAKSJGKuc8PuKm57cR9Ws4FZ49O17Wr/3LOkpJ+iVN339u+fqXn43x58PlodI3pm0dbtISXGY+/6c6n2297jC7nmO2NzhRy9g6GamjKscWw8WKeNlDvKW2ju6qHb46PT5aWyxQFAWZPM5L6roo0ul5fKVnm7MtFRfss+AA6Pn6RkE9A/P/FLe2q56c48fnxdCatflZWz5U1dfFh+DpdXol03gx6bYgk5EgaPguqaFfpmIOoxyWYjEiBJ8OsPqrhnyRQe21RGaX2nNgvbVSGP1tFGuJpWF3a3j0/OdGr9EoyeUCmZI6Csyc6Wo03sqmjTrALB16dXOm88WE9DRw+bjzZrM7qA+xtB/xVpiRFrHwwVdTYnp5sdvHzHAjYerB9QneqM7ql3yqNajwaLES0sclLMGOMNdHsGMBwBXr+f2xcX8srHtTi9fuKNkJIQz2WF6QFTsUVFmax+rZQdp1qwOT3ccuk4KpodtNjdNHT29CtXFQJyHf2VKSvm5Gt293AoyrZyz3XFHKzp4MTZTj6uaccPZCWZAtabFedkAXSotgOXV67LJGDJ9BzGpljCKkbDmXdD7cuwxvFfb5/C5fXzg6vlB6woJ4nS+k6KcpJCnhMNRTlJXDohnRONnRyobkddkBgFLC7OZH+1DWt8HB6fRLLZiNcv0e7wEB9nwGqOw9HTi9cvcemENKblpgKBgj3DGkey2YgQfp55tzzA32RRkfxSLSrK1ASMHms+P53HNpWRm5qA1Sw/+q9/0oBPQlO0xtKn5wuqINxb1RZg3o/F92L19SU89U45q6+PXTAMdBkyooWF1RzHgX8PTLmp77jPPflByPMc7l7GZyQSZzSA14/fD60OD9b4OK2zVXPmxoMN2nmqbiMSDEAo0ZViNgJQ1dKtmeSCYTLA+lWXMz4jkZtmRvZyv2pKFicau5g9MZUtx2SC64R4I89/Y27E84JHklC+APpj9KMQwN3XlpCXlqi9GLGOTLcvLgzpFxDqQQ/eppqVvzJ3XEiHohSLvCx75t1y/rRXXsv/dkc1Dt0gEjxTDDUaX16Y2U9n1drt5sPyVuYVyulZ3zrSwOObT2pOYBcSK+bks/3UOfZV29hytEnrC7V/IplRVy6c1O9eRkKdzUl5Uxc5yWbqet2u6GeMcGHR6fKy6In3uH/5FMYkJ/D01nKWz8hh+6lzlDd3hX1x2xxeCu7bpP1XLSobDzWw8VCfcFDebw3hytMj3P4uty+s8FLh9cN1v/iQP3x7HluONgFQkGUJ8FBUsWHXGTw+SRMUAHa3T3PsUj0cg19O9WG/dW4eeyrbuXdpCbmplrAvbIY1jl++V4HPL9HZ04tRyO381funKciyctMlY/nLvjoKsqy02j3Emwys+fx0Li/MZH91W4Bfw/ZT5yita+dcl5uCrERsDg8Hz7TT64eX99fy4M1TOVjTwY7yFipbHfxp7xm+sWA8f94vOxD+9WA9n9TaONXcTXZyPHW2Hk3ovrijAodXIidZVty7PH4mZljIsMbxhz21ZFnjKc6xUmdz8q11+6hsdZKeGMdvbpsbcRr/g6uLlJlVEW8daeCf/1KKX5IVjer9WLezekA+D8MBe483YOYUK9a8eZQNu2tZuXACD908Q9se7Pm7+rXD2lLNaE2PKaPeiBYWZztceDt6eOTNExRkJbGv2saxhg4cHv+wODQFK+0vRJRMT6+k6QUALCaByyvxyJsnAkeyMMm7Vf3D3qo2yprs/R7gR948QbPdw6/er8Trly0e8ydl9ht51SmoNd4QMEKrq6teCSpaHDz3XiW9kmydUKFaG/T6DXtPb8AyTbUYqWi2u3n478dpdfQpPlsdHn65rVJbqnS5eymtl/0Bam2By0CHshRT40j8yDqoF3b0+U386v1Klk4fS2WrrH9od/by0N+Ps+We8Dmj1u+sZl+1jfU7q/m4xqZd//3Lp4Q9BwgQlMOlx9h4sJ7S+k7mFWaQnGDS7pl+1hYN63fVat+3L5qsCQi9rqfL5WVftU02XfsBYYjJ7Dii/SwsJnnon5KbotmRLyvI0PbHG+UXagh+Wxry0hL41qIJ5CTHM68wPfoJQ4CqDwAYkyxHsMyeGFjnI7dMxxrf//ZMHZsMyEFxoaCWM3tiumZ3XzEnv5/TmLpN35/BSE80kZYoK4otJsHkLCtTc5O1mURuWoL2bXOGtjLokZYYz6qFBQHOWonB07so8AQ53ty/fAoW1cQuQugVojwcH9fI5vW9VTYU72+Ksq0Bgvv2xYX9HOj0Zmc9huKTsWJOPqsWFjAtN4VlM8Zq90xdCg7U7Xvdzmqefe8063ZWs6gok3mFGZpeByBJNQJI/piUgiN6ZqFGz03KStLsyPrgnTkFaby0p5Z7l5awav3+gBFyoNh1n2yiVadu+mXMuJR4znZFfxliQboljruvLSHZLN8o/XJCj5tm5nFJfjrPbivnw1PnaHN4+daiCSAZKGuyc+mENG6aOa7feQ8sn0ZJTkq/5Uk4N3R1eqr35QAwCHjjrsXc8dIBWh1eCrKS2HJ34Ag9NsWifc+flBEQ15FljSfVEkdlq5MJ6QlkWM08cONULi/MZE5BGo+8WcbsiWncODOXB18/TrvLS7xRkJuawBlb+CX0ZYXpJJgMbC9v48qSTG6amceY5AQe21SmCeEsq4lWhxejkJWakXBZYTpbjjWTmRRPRYuDvLQEHvvSjIBjQult7l1aEtK1fCg+GeMzEkmxyDOKwfrOrJgzjo0Hz7JizriA7XqLjjpTKc6x8tKeWl7raq0NU1wARrSfxdrnnn844eKlfK4kk/mTZIVtl8vLjvIW9lXZ+PDkOZo63bR2u5mYmUjFudAmu1ggSQTY5ddu63t51JHL6xv6QsXT62f2xDRmjU/jjcNnSU+M44NTLczITw2wp++vbmPNmyf47hWTmJRtZU+ljfLmbpbPyOGS/HRuX1TI9ReP7ecnkmoxsWByZkT/ERVvHWngjg0HOVzfwZQcK01dHiwmgclo4Oc3TuGKkjGkJcZxuK6DlQsn9vNdmJydREqCia/Pn0jRmCTF5ySZ9EQTT6yYic3h4USjneum5bBu1eX4JXm0m5qbQqolnjuvnMzYFAvn7D1IEnyuJIspuUkcb+giPy2BziBba15aAgVZibxf1kqcQbBy4URmjU8nLz2Rxs4e1u+qISXBRIJJUHHOwaSsRG5bUBCxL8YkmznW0IlPkrh2yhieunUW08alRu27vPREbp3b3wdC7/cRyz3Qo87mZG9VWz/fmYHgT3trOWNzkp1k5t7rL9LuT1OXk4+rbRgNgmum5HD9xWOxmOLo6unlnZeea4zFz2JEJxmyjCuWcr65FotJcOOMcVSc6yY3NYEtx/t74Q0HVszOIy89MaKlZTiQkxxPnNFAQ0cPyWYjdrePvLQEXr5jgTYbWLZ2B2VNdtItJi1qE2Sry6rFk1hUlKkpSfXh0QNZSy964j1tvR+MvFQzM/PT2Ftlo93l1dqxShcCra/r1Y/r2HiogSuLs+h0eUlNNLGnshWPD5LiZX6JVIuJ7adbMccJ3L0SU8cm4+n1B+g6BgKLSVD2yHIgvJVM395Q+Or/7GafshSZlZ86IA/f4YZq9YjW5ki46OebcfdKmOMEp/5jubZ93qPvanqnuxW6BLW++udvP9fb2ZwTreyRvQwxyiO6NT5Os2LUt8dk5RkQVFPn8bNdbDzUwPZT56KdMiRMyU1he3krFpNgTLIZh9tJQ0cPGw/Wa1NPVSehFxQA3W4fz753mpf319KsEK/op6wP/f04ZU12Hvr7cV745lzN/0C/1FGXcXdcUciaN8pCWngaOt00dPYJ5eB2BNflUVwzt5/u79Ld7fFRWt+JWYkKdCt6h2AlaCyIN4JH0YimKGvu/dVtPPR32UqkX5ODHGMSijpAxcEz7bp/kmaBOF/8IJGgRtoONmQe5OVfQ6ebLGtguEduagLNdg/WeAPFOVaeebdci5MRBkNMiqMRreBUp/09Xj8rZucxKz+Vy0IoH60DiCEJBQlZ2k7Pk/1Fz3fo+KQseW3t8kpUtjrxA8lmY4D+4dIJciRkemLgVFR9sZvtbo0TIkBvIfq+1fXz01tlcpiNB+s1j871u2uwOXo1FrFQKMq2asrI9EQTqxTyGhWqQPP4/DFpmccOkUnq7muLmZBh1f5nJMlsYWq8TlmTvZ/CsaLF0S+ASo85BX3PU1mTXeunzyoWFGUFfKuoapEVrg6Pnxd3yIrP8ubuAZU9omcW+gWSykdRZ3MyNsVCU5eLD0620NPrJzfNQmWLY0imT1XZpyoe9eQ0w4l4Yx87lr3HS4vdzckmOw/ePDVgNFOdoxYVZXLbi/twK4IzyxpPq8PDrPxUfvlPs/uNgGs+P11bGuSmygrI4JmF3qNTCD+//qCKsalmam09JJoMOL1+po5N5oVvzqWx0xV2WXPphDQqWxyaYKtscRBvFBqpjB5Txyaz5gvT2XK0STP7AqQmxPXTTYSDyllx/1+P0uny8sOrJwOyslGd3YzPsPDJGVB13QbgD3tqKM6xhnSyeurWWVzx5AdIgN8vxUw1cD4wHJQKwU51GhRhHm8UmpfuQM2II1pnkTL+Iinj67+gKNvKtnuv6rdfv2b+l9cOR9SiR4LJAKcfuzFgm94akmw2cvviSf0sBoNBuGuJhLeONPDTjUdwePysmJOnPQzDOVW+5bmdlNZ3YjEZcHn9Ma2bQ1H8yQFwdXxypoOmLpdmodKXF0zH9+DfjtPu9GI2CsYq1pCJGZaA+ymA6ifke6TS6N0dRGWorsGvLMlie1CEa15agmbxCob+Xtc8cWPIYy4EBkupFwtWrtvH9vJWrizJ4j9umaHE0TjZeLDhH0NnkWmNJ8ls5LYFE0Luz021MH9SJrmpFq6Zmq05pAwUoWbiV5Zksr28jXgjPP7lGdw0M2/IwkIAN10SGz+BHpfkp/OVuXIfDBfXYzCKcqyU1nfi8vpjJoTRmxRVv4LcVAtP3do3CwzmGlXPC2Yue+TNE8yemM6qRYXsqmhlTLecAAAgAElEQVSjOMcaEBw4KavvmqPFakzKSuKj8lZtyZZljQ/rZLV228mo1xkrhvqyD2fQVzASFZ+lI3WdNHa6+PF1JeyvbqPO5qK2p7s9yunACNdZdLp6sbt9bD4a2vqhrsk3Hqxn44GGkMfEglAMcLPGy85K37+qeNhiBCTgj3sGLtBUPUOKxXTeFG9vlDZqvx0efz/GpWjQ3wsValzKQ5+fHrHdp5sdGv/mroo2fnxdCaebAy0kZ3WBfeGclPTOU1eUyGv2q0qyOPDgdWHv4a8/qNJ+pyYMbewM1QcjBTVtilery6vpdVTfC4M5Mabg7hEtLDKT4slLSwg7s9B7Jq6YO7QX+lvr9gWwTxXnWJmVn0p5c9eQ2KiCoUY7DgTFOVby0hIozrFGP3iQ+Pr88QggJcFIljU+bJ+HQ7CHYJ3NyepXS1m2dgerXyuN2H+LijKZOjaZqWOTNU19cY6VeJ2OviCzP7N4MFOYfmTOSjIzKz+VH1xdFLHdP7h6EoojMJ09vaxctzem6w3lqRnKU3agUMvdX902rOxcqvI+yxqvOZKpz5X0jxBI1uP10dDR02+UUaH3Qtx8pGlIdX1Y3hrA5LR+t/ytWkbeOdYYcLxB9MVRDASXFYZ3r9ZDP6V9ZutpGjp6eGbr6aiznMFOhVMSzEhAV48P8PHSnlouyZctBZFSAah6CjVOZMvRJi4vzJRne4q5u6zJTl5a4JJFX+aWo02awvPFHdWU1ncyKz9VM5EC2BzegHMb2p1a+cHejvq6o/FB3LNkCvcsmaLpLUKlDAiFUJ6aw7GMUMvdfuocpfWdwxa89pW546mzuQIU1S8qcTXGxLSo+goY4cIiPTGeVVEk9dptJ3l2W2VES0ic6Is8TTQJnN6Bv+VndUw2Ali5cAKbDjdxrlv2ery8IJPqVget3W6cXr9Wr9nUF6g1OdvK3dfG9jDpH0Y17iKW+IvBuhurvBdNnS4+rmlnX7WNdTurZfr+Qw1sPd7EC9+cGyAw1GjMrcf7XvamThfPvFvOoqJMGtqdHD/bxfS8lIB7GK6Ns/JTA/g0isYksb/ahtPr48GbpwLw89ePsv10K2OSZD+CzEQTv/6gAiH83LNkSsC1qL9jgWqViXUpcr44LjR27w7nsJrwQxH4qH0t+byhPfOCMKKFRXycIeoD/+sPqvoJihVz8njq1lmahn9ilpWsJHM/859eCx4rcpLjabZ7SEkw0668vC6vRHaymQ3flhPbqFr52xbIfgk/+vMhSus7uaI4O+bRXv8wlta1s728lZn50d2QB/sQq/qFZ94tD/CQVUl4yprsrH7tME/dekn/a9ApiBs7e7TzQ6VfCNVGfVRlY6eLOpuLr8wdH3JGsLtSnv2p9HhtCn2fyvSlv5aB4EtzxrFhdy1fCoqpCIfzpYzUz5ZDmkAHieA+V90EVi0s4OGn2qITuTDCdRax4AdXTyJOyAFaKQlxFGVb+cpcOenZd64oJC8tgfx0i0YoMhDkpfYnPU5KMGlrc70RRX2poE/RNqcgjdWvHeY7V/SPWlQRLkpRr8TLVpyP1O9IaOx0sbeqjcbOwZmR1chH1QHrgRunMis/lam5yf0YotXrXPP56ZrTnHqtkR7yYAWl/r9+BAyFUK4BBoHG9DVY/GFXLX5J/h4JGGykaazl6ZXm+LwxRUmO6JlFLPjy7Al0On3a9LmixaGtm083O2jo6An7kkUju0kwmQC5H9MtMhVcfrqF7eWt7Kpo49KJ6eyrbic53sgDN07VzlNvjBrf0Wp3c9MloUesWJYNDpW/0x3deSkWxu5oUJmpVL3C3+5a3M+nAgJH110VbWw81MDpZseQRtxoM6PxGYlUtjo1t+8sazw3XzKOL88emEI2GD7ddyT38H8U6Pv5JzGeMyRhIYSoAezIfdwrSdJcIUQG8ApQANQAX5EkqV0IIYBngeWAE/iWJEmHhlI/9EnIUNCnDgwVBh4toF0f4GTv6aVXgnqbSxs5y5tkopbFxVmhlWjKMNjS7ebZ907T0O7sNzWPZdnQqJgNG0PwggYjXOh0rNDzQOpJb6NNu4drDR+tnlplBub1yc5yk8dYtfs/FEXg1Nwkyhrl2WGktAD/KBifkUhxjpUVz+8mLn1cTNOyIXlwKsJiriRJrbptTwI2SZKeEELcB6RLkvRTIcRy4EfIwmIe8KwkSfNClati7ty50oEDB8LuVz0BNx9t1JSK0WAxGbjqomweWD5tQJGl0WYh8wrTeOXORWGT8MLgoxrPBytTsEVCrWP5jBw2H23mtgUTtETHn9Yoq3fqUnk/fvXe6QBe9JxkM812NwIwGQU/u3HKgLgo9XXd/N8f0dHTy+SsRN5bffXwXMQIRZ3NyfJnd2B3+2jccA/uxtNRnb/PxzLkC8BVyu8NwIfAT5Xtf5Bk6bRXCJEmhMiVJKkxZClhoH/I1+2sDuDUjAUur58tx5opyUkhJcGomApljEuNp7XbGzK2IcNqwubwhhUY+6o7qLM5I5revnNF7KkK9S9KQZaF+nYn5+zRZxaxmk6Dlz/q8qWmtZtmu4dpuSkXbHQN12b9rFFN2KSPJclLS+COKwp5+I0yJMDjk3hscxnXTBk7YPPx+IxEOpRyVVq+f2Q8u01OHWAU4Pe4YjK7DFVYSMBWIYQE/EaSpBeAHFUASJLUKIQYoxybB9Tpzq1Xtg1IWOgf8sGG8hZlWwOsDCqaOz2ES7mp546M1LZI+PBkS8zeoI9tOqFZFVTT730bj0Q9P1bTafCyQV2+ZFhNbDnWzCe17dTZnBEDyfRQZyaxzkiChb5KiHv74kJt+4o5+TS0O6k4181tCyYwf1ImGdY4fv1Blcay9fjmkwHWsIS4uAuaQWykIdZZ6H5liemTQOr1usMeqMNQhcUiSZLOKgLhXSFEJEf7UNOcfkO4EOIO4A6ACRP6K630D/ljm04MuMETMywaHf+5rsA+Gjwpn7xMWVSUybqPqrB7QoucT2o7Yi5vW1mf+VL1EYklf0qsuoNg3YBKW7h220neOdZMaX0n63ZWayN6sMI0eDagcluUne2ky+0LqZ/RI5TQt/d4NcHR0O7E4e5lW1kzXj+8uL0Kc3wcrfYemu1uDp1p53hDV0A29DgBj37pYs2ZbCD6k2Br1IbdVdpy5nwGeA03Qim4Q7XfrlOWC2NcTGGuQzKdSpJ0Vvk+B7wOXA40CyFyAZRvlUmmHhivOz0f6GfflSTpBUmS5kqSNDc7O7tfnXoTkOrvHgp6igs97+01U3K0DitQgpMS4gRxAm64eAxJ8bERyE7OtjIxw6L9X7mwgF0VbWEFBUCHK3YezxhVMP0wVJPb7z6qCRCa4RLu9ouDUIYCl5IOPZpg1LtGhwrNrjjXzZbjzVo/HG6QhZa6RGi2e2jploW9Stx8y+w8bpqZN6g+CJ4VPvRGmfb72W1ylOuz28qDTztviCXJcyjctmCCFiKhlqG2X3+N8yf1eRILY1xMRCODnlkIIayAQZIku/J7KfDvwBvASuAJ5fvvyilvAHcJIV5GVnB2DlRfEYzpeSmUNdm5qiSrX9Jd/cumH5Bf2lPDn/efQSDISJQvv0cZumtaXXQrL/uy6Tk0dvaE9aLr7vHyuZIszthcxBlgR/m5qCxebQ4vbx1puODJa2LB6tc+CUiyNC4lXotwDWWCDZ7BqDwaTZ0uzthcxEchJNLPbIKp7lMsJho6HAF9n51sorXbG+Bif9+yi3hhRzV3XFGIzdEb4HA00JmA6vUZyrKm+tDofWnONwZrAj9Y00FDRw8Hazq05NldruR+vi+rFhXyUXkr3R4fkt8fE6HIUGYWOcBOIcRhYD+wSZKkt5GFxHVCiNPAdcp/gM1AFVAB/Bb4QbQKPL3+fg5LeiemGXkpWEwGqgfA4dgrybRuPb3+/ozdAqzKzGLL8eaI7rbNdo/2cvX6ZaWYOwZC3wdfPxZTO9V0BMETnWgjTSgnr1jo6YOzsZ3t8kR80cZnJGr5KOpsTm0Z81+3XsK8wgzWfH56TEFRqiJXdc9WhcjfDgW255w9UFAYhJyF6+U7FrC30sZvd1Ry1X99wNptJwMo8AcK9YVQZysAD9w4lXmFGQG+NOcT+6vb6OrxMnVs8qBN4NCXcmJ6XoqSqvOw9vzsqmjTBkYRZzq/MwtJkqqAS0JsbwP6sYwoVpAfDqSOdqenn6JKrwzbfPQsLq9/0KQ3KgQwLi2Br12Wz94qGzvKW4aUViASHJ5e7QVp6nRR3twdkOVLxb1LL+LpreUK7X3frCnaSPPstnItoXLJ2BSKc6w8+Pox2l29lNa2Myk7iaZOFzVtTqbnpXD3taGn60khcpYEQ693UAXHijn5PHXrJWw8WE95Uxdbjjfz1uGzVLY6NAWmnuNCvZ8An9S288CNU9lytAlTnKA3RAxPHGAwCs1jc+PB+gD39F9uq+SLs+WZ20AU4HrLS7LZyONfltMB1Nmc7KpoC+3mfp7w9NZyyhrtzCvMGLCpPHiWprqMqxnI1OdnxZx8/vZJPWdsLqReT0zmnxHtwRkqkEwlGW3qdFGQZQ3IlDVYpCTEBWTfXrWw4LzR6vn8oR3JgoWA6vY8SxcPkpNsjjrS6BMqbzneTF5aAu0ueZZ54Ex7AKGuPhpUZRnX2hmD+41+KRKce+TZ906Tkyx7znbqXtrgrO56lNZ3atPvcEhONPHJvy0NaMOmI2epUDK3J5qNg6Kn0wef6QmGLqRlRV0+qfQAg5lVBCuu1d/BznrjMxK5ZkoO63fXxGwNGdF5Q3734m8fXv+f9wfkT3jy7ZN09vTKFPLn7DE91OFgEmCJN7Jibh5ur59udy/XTBnD968qYv2umoBjVy0swO7y0O4MP1rFQml495LJrJgzHiTITTEjSTA2NYEHlk8NyEGh5p/4+vwJdPd4sff08uDNU7l6SuRo4knZVmptTn50bRHFY5L5wqXjONVox2gQ/OS6Ygozk8hNMWMyys5pty8qJNViwuPr5WBNO9lJ8bg8Pr6xYDxXXRS5Ln2OEn2+jGnjUkhJMPHF2eNo6/Zw9xK5LV+fP5Fp41JAQsuNcemEdOwuL5IEV5Vk890rJpFkjsPW7e7HzZlsNvKTpcXMGp8e0Iarp4yhqcNFh9PLgsmZXJyXwunmbr5yWT4lOZF5XfZXt/GTVw8zIz+V1MQ4NuyuoTgnSTtvKHlABopH3pQDEBNNRn678rJ+OUmGAr8EXT29zCvM7Jf35e31v6iKJW/IiJ5ZqDoLVVFVZ3PS1CXPLDpc3pAMVwOBVwKvxxdAx1drc1DX3n9Wtre6jTNRnHVikVshlaAhTmzsdGkJoD86LSuifvxyKZuONDI2xRIyMXKdzcmrH9fh9voYk5zATTPzqLM5uaJEtipdM0Wm9Fu3sxqrOY66NidLfvEhWUlm7K5evP6+nKavfdwAkkyv9+KOaopykrjqomwtA1xwYuRz9h42HqwjwxpHTasLe4+XvZWtHK7vpPKcHYQI2Ke2Y9mMsRw/20VZkx2Hu5fjZ7soyEqkI0goxwk5MXRNqysgOfSWo03Ye7xYE+LIT7ew5VgzeyvbaHf18uDrx7kkPz3i8kGfr/VMm0M575imhB5sdOlglKyqBWkgJvZYMRzJnUe0sGhVYirUte7q1w5rhCjJZiNdMbJCDwQeH+yr7k9JWNY48BwXobDx4Fny0qxRlyFPby2ntL4zQMnq9aOlCFT9HyAw2bHq0aqWF2raH1x3qERD3R6fzFVxIoGGDtkqtKeyjYaOHq1svca+vl3OffLUO+UByxmAHiWnSKh9Jxq7NC6MWkWQhswnIsG8wgyNXwPQtP16zCvMoKxRjtlpd3kDcrGEgp4/44yqgA2TlHogGMzyRbUgRbMkDRfUNhoTU7OiHz3ClyFP//evHzZOu45Lx8uU83891ECWNR6n10d+RiKt3UPTV0xIN+PxSSSZjfj8fiwmA7mpCUwZmxw2U9dQMTU3iZ/eMDVwGZKSwAM3Bi5DJmYmcrpZzmmamRhPa7ebsclmFhdn8bmibL65cCLj0iwB0+P8dAt2l5d4o4GfLptCXnoi+emWgGm/ugwoyraSbI6Ty01JwO/z49WZG6zxBr4xbyJfuSyfs+0uFhdncduCibR1e7h3aQl56YlMzEyk1ubk3qUlzJ6YxuG6Dn50bRGFmUla+W3dHtIscSSa4/ixsgwqyrYyLTeFeYWZfHPhRCRJQpJgem4KifFGLitI75eK0hwnqLG5yE4yk55ooqnLzeLiLOYVZlLX5sDl9ZNoMvDkipm0OdxUnHNQlG3lgeVTIy4fpubKS6bbFxUybVwyh+s6+Lebp0VdvkTDYJYvJTlJ1Nqc/ZakwwF9qkn985KSYOLt3z9T85lPXzhz1mxp1X+93I9qfldFW0Ak6fiMRArv2zTgvCGWOIGrV+LK4iyNuEaF6ncggC/PGUedrYeK5i7anIOfzVjjDfznipkjzs+izuZk2dodmikNYFZ+Cn+763OfWpv0xEQGIDHeSLfHR5Y1ntd/uChgij/5vk2am/68gnRuWziRxzef5P7lU0ZcX19I6JdCodz21f0/WXbxUanXMzNaeSN6GRLMlKX+Vi82N9WidcZgRJ5LccbaXdkaNm2dBHxypmNYgoscHj+PvVWm6RKCUwvq69YHkqUmGvntjmouK8jgB1cXhTxexWAiVNftrA4QFCAvD9R2ROPgDCXMw20L1eZQ5ejhpy9dQ356Qr/2pCtZ0wH217TT2NlDQ0cPD/z1WFSdxfnCSHAR1y+FVOrDH/zxEL/+xmx2VbRpPKZGa/pnn4MzFPQvw/qd1Ww53symIzGxgoWF30+/9aXeSWk4oxBVLs9oxKx6/wMVH5a30tzlpqzJHjb2Qp9/dMs9V4Rth/5hDuWP4PHJ9IA7yluobHWErE+vNAMC+Di3nzqnJbHW+1k8pDhrqYJwTkEaD//9BK0OD28dPkt8nIGyJrsW76GiS9F1VLU4Ka0PvFepifGasBCgsYR1uXs/NW6KkRDMpjdtv3ZAVuK3Ojzc/9ejVLY6GJcim7ZjjQ35zAkLvfZac8NtGVwWbhVqrtMLmbZusMSsTV2yLuX42a7QB4ig7zBQH+Yul5fdFf2TGUOgIjSahl71f3F45GVaaX0nZ5TYHTWGQ1+3Wvabh/tmBXqyoXAzRY/Px7LpOby8v1ZLSViSk0Sl8gzEGWRfFhV7q9ouGPOVXgCfL0LfgbRB72Xr0sU/qH4vbQqHbKyxIZ85Ds7ctATt2xUhaGsg8CNL/8P17cx79F1Wrts3LOVGgmqSu/vakpD8nLcvLmTVwoKAbNjWeIPm56HmgdBjf3UbLYowqTxnZ8nTH4Z1D1dzRjR1ugKYy/Wwxhu1xMyhNPT6pD4qi5fVHMey6TnkJJsxKqdYTAbmFWawbIZsul0xJ19zNgvOo2KJYgnwS7IAaLa7Ndd5a3xfGdPGpQTQDJQ12WNO+hPOLf6tIw0seuI93joSmTtFH1w3mGC2WNzyoyE4wE/9P39SBjnJZpZdnMPDn5/GvMIMFkySl6mxenCOaGERKjZEG3KkvnXsUA1davjF45tP0mz39MuTOZxIih9Ya/PT+4R+akKcRoyrkhLr8fTWcm2U9vjkGVdwVnEVL+2ppaGjh8bOHq4sDm05+9cbSrhnSRHJZiNfuyzyCPnF2eNINhu5fvoYDtW202x3a22x9/Syr9rGqx/X8cy75TR2ushNSyAn2czYVPn6DMj0AV+7PB+zEgUcCpOzrTiUQSK0Q54gJ1kWsFlWU/8s8xEQLqPYI2+eoKGjh0fejEyJMNQkQ8OR0Uyf7KnO5qTL5WXVwgIeWD6NfT9bwvPfmMtNM/N45c4FXD0lm2SzEb/XHdPUfESbTh/7xS8f3p8wm5QE2VMQ4MWPqmjqcpNkjsPl9dHT68dkGFzCHxUmgxwMtrg4k7KzXVycl3LeTKceH9yzpG8NqwY96a+xzubklud2srvKRpNu1Le7fWQnmdlX006SOY6rLhoTUPbEzEQ2HT2L1ydhFPKL9eBN00Ka4VSz56rFBRw60x7yel1uH28fa6LD1cvh2k6au3qYnJ2kmd5Uj0O7y8sHJ1to6OzhWH2n5l6uwidJSEC80cDfD5+l1ubkw/JWHB4fbq8fp9eHhJyust7mxO72heUW8Xj9GoVifloCT759kla7B7u7l5xkM49/eQYbdp8BwOn1093jZd6kjJhMkaHMnWvePMrOCtmfxWSA710lZzhT9S756Rbt2FNNXbx2oJ4Z+amDMn2q9WdY4/jhnw5R1tjJva8exuPrZU+lLaCucPiPt07wYXkrh860Y3N4lETR2aRbTfzk1cP8aW819/3vMT6ptfHW4UY6XL04jn+Q+OC/3P1ItPaN6JlFemJ8P0mt0vt/54pCnAp3wmC5H1T0KgpOm6OXfT+7jvmTYvJRGTTUmdL+6ja2nzrHsuk5dLm82vaNB+u1fBjBOHAmfA7bywsz+aqSQPmbCwrYdu9VYS0iapTowZoO9tWELvNQbYc2O2hXwrf1o15A6LbGZ6EkWDKgjfB+Ce2erVpYQIbVRLry0E/PS9F+Q/SYlF5/382utclZ2hs6e5hXmMHG7y/k8sLMAH1HpNlVMEItHTbs7vPu7dBRMIaaBeid1AYDtf4XlExhGw+exe728esPqmKecRxSdEuVrQ7sPV5tlqG2rbRednjbXt5GWqKyxPX7Y1rPj2gFp950umF3FU++XY413si5bg8fnmwhP81CZauTFLMRh8c36DiRywrTmT8pq59SalFRJluONtHU6WJ/tS3gBQ7W1gfzeUZCMFGwquBs6nTx/G1zWTEnn63HGilr7sYIAWvwf72hJIC7IRhzCtLYeiJBy4JeZ3Py2KYTHKrt4MGbp3LTzDz2V7fx2KYy7D3eiJaeH14ziXpbj+aCPTbFElDvAzdODQhOenprOQsmp/P8h5W4e6FNF+TX0NHDT187jIQIIFfeV9Wm8YkAtDu9GAQkxBlCkjAbhNAylyUnGPnjnlpmT0wPGxValG0dUpj3yoUTAsIB9lfLGb1CKTCHyqyu4v7lU3h880kWTM7gnWPNfPtzBUiSIablzYM3T+WRN8uYPTGN5AQT+6ob2FXRprXN7fVSWm/nypJMfnh1MU9vLee1rpYzsbRrRDtl6dm9Zzz0doCr8NSxydTanDg8PqzxRv5zxQzu+nOptn9WfrImRVUYkAOtPL0+att7WHbxGEpyUmO2hQcTxOgZv6fmJlOQkcg7x5sHTc+Xk2xm38+W9Nv+1pGGmJ2MvvqbPeyrtjGvMINX7lzAM++Waya8vLQEdt13rXZMKJgMAq9f0rK6DQbzHt1Gsz2mQMawSE80MTEjkRa7m4bOHqzxchrIomwr2+69CkC7truvLQ4wTy57djtljd1MzU1iy91XDqkdAN//4wHNzV7t188CYvX1EEIclCRpbrTyRvTMQo/V15fw5NvlWEwGWh1evD4/X7ksjz/treO7VxRw08w8XjtQx/byNq4syWTD7fPZsLuKp94p5/qLc9hTaRuyR99Tt17KjLxUnny7nLkF6UzKTuSlPbXkJJs1Pgq9D4Fq4dh4sJ63DtdT2Rqed8MUIavWTTPzYm73bQsmUN/u1MKcFxVlsvFALc12N7cqmeZvWzCBynPdeHw+bTY0dayVpdPHhc2xMhA8ePPUAME9GEzMsPC3uxb3S1GgH7XDmSdfuO2yfgmRhoIHlk8DSZ7iDzS7/KeJ4U6x+JkQFnU2JzZHL2/fc4VG7lLR4qChw4XHJyFJsuqlRyFLUb9tjl7sbh91tp6I2dgHApujF4fHp4VJ9/phxdwJmm4gVJ5N9Ybp+R6C4ZXksqMh2mihZmFTr3VXRRsNnfIov6eynXuWyMe0OuQlQvCoDAwpN0mdzTks/VzfLitcVV6P+ZMy+43o4V6G4X5JxmckUjI2hS3Hm4fl2j6rGNEKThWqxeCOPxwIcA5yef0YgAyrLPOWz8gh2Wxk+Ywc6mxOGtqdzMpPZcHkdPLSEijOsQ65Lap/QnGONeC3ijqbkzVvHGfNG8cDTL4NHZEfsnSLKab2BSvW9Lb5/dVtvH6oHnOcYG9VK2veOI4QfYuixg4H8x59l9cPyRkZLCbBW4cbKH5gM2u3nexX3mAQTIKjx0DoIFSTcag+/jSQYY0j2WzUnrU6m5PVr5Zyy3M7B0yqe6EQq39IrPhMCAsVZU32AC8/kHUGj26SmZh/91ENdreP331UoyUgKq3v5I+KT8GLOyJzMsbyojyz9TQNHT08s/U0Lypa60feLNPOeWzTCdbvrmH9brkNapnBHJfBaHd5+dX7lVH7INiWrxceT28tp7bdhbtXYl91O+t31/Bb3TXXtrtptnuoVUZtl1eistWJ1y9pdQ+Fv1JtX15aaIdA1wDSvJxpk5dsah//6v3KITssDRZ1NiePbz6J3e1j7buyIFTpAFR2r+FC8As+FOH9+OaTNHT08PjmSBk6YsdnQljcvrhQI9INBdWRU3UrDnYvVk1EKoFpOMTiFKO6ynb2eLXymu1u7ZxDQW7RkUbaYJwKxeMQhGDznl543Lu0hMlZVtITTVxZksWqhQVcVpDRr4xQ/k6zJ6aH2DpwjM9I5JmvziI1YWgr3BSLfL52zwQxmw8HS6MfDhsP1msWG9XitmJOvuYgN1Trhx7BL/hQHLXuXz6FvLQE7l8+ZVjaNqKdsp779f883DP5aqaNS2FGfgqH6zqYlptEQ4dsyahtc+CTZLPl968qIjnByMGadv71hov40uzxGrfDD66ZTPGYZI1CLhxi4SAYl5bA4boOHrxpKl+YlR/AFZFqMZGbZuZgTTtXFGfxz0tKNIo5t9dLk8ImbjbK3qfBjmQS8Nz7p/nFu6c51dQZoNSsszlZ+24520+1BDhG6ant8tITWbmwgO9dOZkvXprPVReNYc7EDD6psXHO7ubLc8Zx3bSxLJuRw7H6TtITTQrXahIAACAASURBVOSkJHDNlDH8283TSbWY+vEehHI+0kNt1/MfVvDwG8dJMhuYmZ9Oj9fP8fqOQVuGspPMXD1lDPMnZZKSYOKGi3No6/bQ2t3Dj185TH27g6XTc0O257t/OEBZk51am5NbdZ6u0a4lHPLTLbyy/wxeP7h7/bxz7Cw/uLqYpdPH8rXLJwwL94TatsXFmZxu7ub+5VMoyUnBaJCTQX/1svER61HpAS3xBjYdaSI/3cJlBZl8e/EkLKa4iNe9Zs2axocffviFaG0c0QrOYHbvYIvA6tdK2XiwgaXT5ZiDlQsnBSTFzUtPZOOhBk43O2JSeMWiGAu2TAQrM0NZLn58XQl7q/pGOYPBgDuMJ5lCKsWWY+cCtgczXsWqwBufkcjfftSflyJc8uDgPogWPRlMPvzUO+XYHL2s313TbwYTZ+i7vmioaHFoLFc/vq6EZ94NJPLdePAsT916acj2NHT0kJeWEDYxUrhrCYfxGYk4dUzjZU3Dr+RU23b3tcXsuq+PHF9V8O6qaIuoeFadrlTGMghkUBuOCNgRPbP4/e9efPh7d94ZdqTXsxyF2n8hyVajYWJmIofOtGMUgnuXKsS5qWZcbl9AKjk9fvPhae66Rr65KuNVboqZYw2dTMq2xjSiDXY0VRFpZHvrSAOPbT4REOV56YRUUhLiKcq2IiFpbGZ5aQk8+sWLSTQZcXv9ZFrjKR6ThLvXz2UT05mWm4wkyUzraRYT41ItfHPhRPySrEdZVJTJuDRLgMBYODkzoE11Nid7q9q4dHwaj35xBtPG9TGjq3042OfhVFNnAHvX2m2nw85uBoNwbYu1zar7/g+VWXQwg1qkMmKdWUR1yhJCrANuAs5JknSxsi0DeAUoAGqAr0iS1C6EEMCzwHLACXxLkqRDyjkrgZ8rxf6HJEkbojVO75T1j4xghzM9ap64MeB/sNNVNIRzXIoVkc5f9MR7/WJKVM/WeYUZ2os9kLr1TmR3X1sMEFC/nhEtuA+Geq3RsHbbSdZuC1RCB9+fzyKG0ynr98BzwB902+4D3pMk6QkhxH3K/58Cy4Bi5TMPeB6YpwiXh4C5yM/SQSHEG5IkhQ90GASCWaKisVHFWubqVw9zzu7mkvxUpo1L1fJKxOIhp5rYDtV28MNrJnHPkv7KptXXl/DY5pO4e/sLbj29HMC8wjQt92gsrFix8CoEM1U9u628z0QtJATwu48qefa905rDG8gKtJ+8WopHNzGalJ1IVpKca3NabkrUuoOxqCiT7afOUTQmSYlpOIXJgGYCfvjzU3lscxkJcXEsn5ETwHCmv9bv//Fjthw7x7KLx/D8Ny6Luf5IuGfJFPZUtrGvukOpZ9ywlBsKg2E8U6HezwxrHC/sqA7rjKh6BhssyTFpt2Ny9xZCFABv6WYWp4CrJElqVJIffyhJ0kVCiN8ov/+iP079SJJ0p7I94LhwiDazUGMcinKSmJGXwr+/UYYPmJxl5fe3X8631u2nstXB5Cwrla2OfiNOqJc9mLdw5bp9uIIyY80rzOCpWy/RHMRWzM7TWKT2V7fx0N9lqvo1X5jOr96v0BL7CKA6ykgULBxCQR3Nlq3dQVmTnaljkyOyYkWDfkSGyM5jav1vHWnggb8ew+HpDYjJGZMUz7luD8um5/DAjdPCClP9ywBov3dVtPHse6fJS0sg2xpPaYNM8hNvhPJH5Xrv/+tR7G6fNnsJvq9vHWkI8CAdyOgfi2v9ht1VPLHlFNnJZp669ZIhObGFQp3NyRd/tZNWh3fA91bfPxaTwOWVSDEbteXami/0Zb6b+dDbdLl9nF3/z72e5sqo67LBKjhz1KTGisBQY6XzgDrdcfXKtnDbhwQ9Xf47x5q0gKvOHpkCXvXJ6OzxhuQZCKX40W/bW9XWT1AA7Ku2sfFgfciEuU9vLdeo7J/eWs6xhj4WrFiicOYVpmkjV1TEyIoVDcGzjz/tPaN5eIbD45tP0qXoWkwGQarFRKvDg01hXzpU2xFRsRYcoan+llMf1slcG7oljlAuUvV3SDYbNeESfF/1fgURLO5hr6uho4efvHKYMckJIQXBo5vK8PjkqNe7/nSI/T+/bmCVRMFjm05o0b4Dvbdq/5gMaM+uxycFPJPq0s2jSHkhREwuFMPtZxHq0qQI2/sXIMQdQogDQogDLS0t/fbrnVTGZ8jOPxPSE+jx9q35Oxweunrc2sV1KlaVH/3lIHU2p1ZGQ4cDg4Dy5k6++ps9vHWkgfLmLnKSzRTnWJk2LrRfhkHI0+LcNJnDUP0GAs75uNqmEbWoiGb7f+XORTz3T5EDuG55bidvHWlgTLIZa7whgJhmoF57+6vbWPH8Lp597zSfe/IDPvfkB1EFxYbdVSyd3selkWQ20q6coxIS+fx+ypu7mDo2mfLmLs2jVfVwdXt7MRlgweR0po1LQgA1rd2s21VJi90dEKQHIEkSMx56m5QE+e0fk2Lm1QO1/PL906zbFahHWDC5z7fE42NATk3qdXl8kjZDDHaMErrH+ZwuHcVw+Xd8rKMMuH76mAhH9sf9y6eQYjYG0DbMm5SOxWRgTFI8bd1uVr9aSp3NyZSxslfs+c512iyEyNUtQ1Q7Xz2gp3DKB84q268K2v5hqIIlSXoBeAHkZUjwfv1o9briFal6JKroleAPu2q1h03tuNK6Ln7050NcOiGd9btrNF8H1UypNzu9tKeWE2dDc2P6JfjN9ipNcuvNnH/e1xfOHMpKGC2xMUTPtF5a36mNgACbjzZrplC9U08swWdPby0fcL7Yp94pD8hTqie7UWVjq8OrRWqqo1qoJEevHWigy+VFQs6Gpg8H18PjB4/bp5ktK1ucVLbIz/iG3bU8dPMM7dh3jjUHnDsQs+HW4zqTtSITgmdID9w4hYfeKOt3rn62NJTI1MsK07W+++Oe2pB6rnC4aWYe9//1aMC23RVteP0QZ/BxrttDRYuDvPRELSpbmMyRvRUVDHZm8QawUvm9Evi7bvs3hYz5QKeyXHkHWCqESBdCpANLlW0Dht5jMdEsjzIWk4FEU5+0F6DtC4bKHXH3tcWsXDiBeKNgxZxxzCvM4P7lUwK88lZfXxK2g8YkJ5CTInuGqt8ACXF98jcUNVwob79+8SQxZMRyenpJt8QhCJzNDNRr796lJRpJTaxYfX1J1PJNBlgxO0/7qPR2K+bks2phAVeWZJFijsPp8XH9xTlYTAYmZlgIR8FpiZOz0M0rTMOAnNdkxZxxGITMORHcPj0GQnV3//IpZFnjmZxlZY3iQxPsYh/OR+XepSWa8nkoeGD5NI37VCOoGQBWX1+CNd5IllUuY/bEdOYVZrD6+hJm5afK92QQEbmxmE7/gjwryAKaka0afwNeBSYAtcCtkiTZFNPpc8ANyKbTVZIkHVDKuR14QCn2UUmS1kdrnD7JkF5Bpio23T4fFc3d/PCaydyzZAprt53k1x9U8fX542nqdPPByXO4eyVuuHgMNa0uvD4/s8ancfeS2IlUVT8Fe4+X5AQTcwrStHyfZY2dPPVOOauvL9EeoLXbTvKr9yspyknCbDTyxdnjeOrtcuweH/MK03nlzoX9FKt6c+GK2XmU1nUMmLH8o3+9etAM1nqCnP83L5/XD52lvt2FXwq9VkxPjOOWS8fx+121EfUwE9It2Ht6Kfn/5H15fFvlme5ztO+25EXedztx9pB9IQlJCiRQaKdh6RJoGAq3BQpT6C3QaROGlnY6oZCBlpJpSWkpbSEdGiAJhKQkgSwOceIsthPvsSXvlizJ2pdz/zj6Pp8jHVnykrnh3uf34ydyLJ31O+/3fu/yPDk6+IMR3L+qFLUdw3D5gnD7Q9hfz82eUgBVOXqUZGpQ3+1Apy2e4k+vkMKgkaMiS4cjzYPIMyhh9waRl6ZCr9OPb60aJYgpNGkoPcHKygx82jwkeEZTgVk/3o+RQAQKKZClV42L/oAfBC/JVOPD+n7cNDMbNndIVEd2ogHUZNk6wvvR/bvveAIDV5J26l3T5DeFVbNY6T/9e1y0O5a8heTbSd5fr5QKouVmvYK62lOZgxereYg9t/x0FV0ubFlegq23zYyrB+BzYBD9TrLPLpsH658/DH8SGrDJXBffWCmkDA18jQWxcvVk4N+LiUKrkMbFgQDOkwlG4u/DzB/vhzsQgVYhQf2/bZjUsfkgLyIJBpv1CtT8MLVA51jkQ5PNbI0Hz7zLcahaXrmvP+ToSyo09Lms4CQ6oNNy9DCq5VTX02xQUr3NuQVG3LWoEHnpaihkErT0u1Gdq8cPbp4+ZdWcfK1PkpqKPbdH11fCrFdhfmE67lvJVZrGVtSlqeVYMy0ba6ZloyxLK9hnmlqO0iwNPqrvHfPl3H7H3AlfV4FRjd5hL1y+EO5bWQxfIILiDA1CYRYAi2DUeEgZztOQRl3/c13ckk6WgDC5yKiGhGEwtzANWTolvRcVWVoUmzTodngRDLPI1MpRnKHFohIjvIEwXCKC1woJkJOmwnejWqpKGQObOwC9UgpfKIJFJUZsWlAYN1aONg3AOuzDdUXpgj6RyYL05Pz1s07YPUHkGJS4d3lp8h+CGyMXrQ5k6ZRYVJKOtgE30jUy+IIRZOmV+MbS4ik7z7FA+oA+2PXLtlS0Tq/p3pBY+UKC3DQ1Vk/LjnOvsvUqFBg1qM5No0VZALBlRSmqzIYpl5LLTVNjaVkGctPUdNvi0gz8/eGVgu/Fuqdj9aAQIl0+5hYYkZumpirjU41Ck4bWRHzlugI8tn46nTmtw27aXl+SoUH7kAf3LC/CfSvK0dbvwekrdvzvm6vwq3+0oH9E2IO+rpqbrDbMzsGxliHMLTAK7oWYqx1bI0GQk6bC0R+si9v+/nkrnn2vESq5lCqj8WGKrvlNE1j7pwLCYToe1brFpRnUezjVPgSbO4Rl5Ua8fdqKh24ovyrnKQZSyPb9FFXUr2ljQXRDYl9yvmwev5ErNhqd6HtThaspUce9BA1c6ziLpIbi5hcOo/7ZjRM+Xuy1kHvHDxO3DXrAgsv4tPV7aLHZv++/JGi0IiBZD7K0Ivsm2H+hFzXtNuy/0EuNRSLOkU67jyqL8WNTB+r70efyo8/lx5GmQbxZ04HLPx0twiLiR+TzauLe107S6tZUEdsAVtsxjOY+95RMbMliFl02D554+xyk2vSUGlyuaWMR23WaDFPFrpwqJitRN9bDJIJH+y/2waxXJtjDKNxBFiVP7oUUACMBFhQbsf2OeThnsXNsz0XpePqWGSg0aURn9ETXoonGfwBAHo1nqGQynOoYXXOLFa6lqWRwRJcTFWYd9QSTocKsSyjnSO7Vfx5sRQSgAVazXkkJgmNbbGIZyKca+elKWIe5Yx9pGn99BRmzm5cVobmP04PdcagZVrsHXXbvpAKcsRNA7HjbXWtBTbsNkaBvJMmuAHxOYxaxnAsE+UYN7lg42h1Jvre2OhtvnbZMuPMyEfhcEhOBmMAQgdmgRG2HDauqMnH/qlJ0270wauQY9gQwzayDViGFwxu/tmfBxQ+swz4YVHLsPNqGHqcfLQNuepzvvXUONe02Ad9D7LWQe3fnogJ0271YWZGJTQvz0dw3gq23zYDTG0THkAcMgG+uKEKuQYVhTxCLSo1YOy0b37uxCjqlDPML0/HtNRW4aVZO3H0Se47VuQa4vFy2RC6RIBThJB5kDPCLTXPx2qftONvFVbjOKzTgrkVF2HbbTPyl5gpCEUAtY/Dw2kp6jNgxMdW4cUYO/tHYB4c3hNVVGfjy/PFNHOT8qswGLCvPoPfkbKcdtZ3DcZwc40FsbCx2vJG/7//df1xOJWZxTWdDpqrr9Gp1IyZz85Klv8ZD8Z/o+LEaJGJQSIF107lejR6HF1vfrQdYYZ/ARI6dCs08EH+d4/ntrB9/gJFAGDqFFBf/7WYawTfrlXjpa/MFsY7J3MvJYCqPTe5NpVlLU/SpPqMXD17CSwdboZBLkJemxs++Mhu5aWqaaSOxo9j7nmrX6eeCVi9VxJblkn9XmrVUmWk8v0/2HTG+Sn6B1RPRGfypv10Q3fev/tFK+SUncl6prmkDYSAnGoR9/kATGntcuNzrwq/+0TJhTstUhX+7bB58/+1zsA77sG0PV3T2yJtnBPeNXN+p9qG461RFq7TI530rSzGvIA19Lj/2X+il3yParX88IV4Bmgomyne5bU89rMM+fP/tc5PmCCVLh+Y+N/764LJxGfOXD7UiDI7IunWQU2Ij5ES7jndg/4Ve+swmcq3XdMwiGWJbqx+I0qlZ7R5sv3MeDdJV5+jR2OsSBNJif19o0sQFRMVmQNJparWP3mSXbzQKT/YBjFZwXhly03NzeoO4b2Updtda4I72dgfCY9NH8Y9JulsJlDIgAXeOALuOd6BtYASmaFVfBMCR5kE88uYZvPS166Y0S0Rwqn0Ij7x5hsY0lDIJHnnzDI1JnGwbokZ31/EOWrNw5HI/PSdCDEQ+C00a5KapUGdxoNcxGvSdkafDZx02zMjTjctz4WOiAetwNG/sDbKU3WuimEwcLEOnEJTub15WhGy9Cv91tBWeYARnr/CZxrhrPdk2BEjlKaWKPteeBZ/MdHethfYg1Hc7hV9MUD2djAxV7O/8TlO9invxyGccoiu8EDvaH8HfL0GBUR37SwHEulsJlLLU7f3pK3bac0BQZ3FMSrV7LMT2nUilEkHwsrHXJTg2aWDjn9O8wnTBJyCe4XjrMwsiLPc5UZLbiaqg83tjknmviUCa0Hoc3pQ8NjHEki4397lxrGWIykDW94yOwU0LCmjRovT/hdRpMsRa4ffPdaN10I2STO5G37eyFAa1PKHKVuzvyfdjt8cKM/9s3yXcv6oUcwuMgu+TfRDsPt0FVyAMhQS4bV4+9Cq54O8HGjg32pFABJlgrIj+iooM7L/YBxnDGaVEyNQq8Mi6cnQMevH3OivsniCMGjm+NG9ifQKp4PEbq/Dc3ka0Dbrh9IUgl0qwZXkJ9cT0KuG9E9s+IzcNNe12zMgdpcjjPwOCIpMGjb0jKIohwRkPJipOxNejTcaVmQhT0YSmVYy+zktKjfT63znThU67D8srRrtxC02aKB2ABd/b7hhMZf+fa8+Cj0KTBquqsgAAOQY13fYvX+ACRGLWOnbdnezfAFDbMUzz4YnOY+ttM7H1tplYWcUZ7HUzzHh0fZWgU9PpDSJbr4RZrxQMerG1JCnUEhuEt8zJhVmvQLpm7IzMoDuADJ0SC0rSEQxFIGGAKrOOsn5NFonO++8Pr6SzLTmeXiWnhpMc26CW49H1Vdh+5zxsvW0m3b5hdg6WlJqwYXYOPQ6JT/DVwWxR/gebO4hzFjt213bhnGVKidgSQheVPVBKmQkb3qloQiMe7pblJfjrg8tRaNKg0KSBPToZnbkyOmYFS/hwMKW242s6dbpz585tDzzwQMK/v/hRE35ztA1ggTXTsqFRSNFp8+Ce5cVxRK6TIa3l45WPW2Ad9kEpZWCxewXHj8XcgnQYVHI8uLocz+1twBs1negd9uJMpx1vnurClSFO2PlgQx+WlJmQb9SMmU4Vw/9++zxaBz2iiuOx2HehF/su9HKSj+DSqw1WB1482ASzQYkqs2EitwSAeBq4y+bBd/98BocaOV6SDI0CFrsXu453oK5rGCdbh7C6Kgs7DjZh1/EOuLxBVOcaBM/qJ+834HDTIEZ8Qdw6Nw+vfdqOv52xYkmpCY+tr6LPMzedK/X/4a3VePa9BvQ4/ajtsOFbq65+RaRGIcEnlwcRZAGXLzAhEl9+inei43XEH8S5rmHcuahA8Cy1SglqO+xc12lUdpM8r4ZuJ+r3/2EoldTp/zOeBSCkTedjMkItsTgfZb46bxUvHOKD75kQ8aEzncNxsYdAmKWMUeNdNxMRnoosLVZXpbT0FOBM5/CUqFaJnffuWguONA3SztSzXcKZnsQsSIypvtsZ96z4941/nO13zBV4RLfOycexJ9fh1jn5dO0uJpw0WXlGMdjcIboMSaY8lwomOl6JeltsFey9y8tw4ZmbBV23/9/FLFKJMYy1fSKYU8CtoecUpMUdfyz86IvVNBefrVfhub2NgmAfcT/Hu25+dF0V8tM1ca3ZpCV72g/3JexYlUuA2+fn4sOLfXhgVWpNUIkgdt6bFhRg9+lOKsyslEpx38pSuHxB1Hc7MTOf69dp6nOisdeFksz4eAP/viU6Tiye3jiD9gLF4mqU6G9aUIDfHGmBP8TCqJn8KzXR8UqqX5Mp7wETi1l8ro1F7MCZavVs0WMaNahpt6PQqBnX8WLFh/7+8EpUPr0XwQj30k60OIp/Dl02DzoGvdi0oBBrp+egy+bB2unZlDciFsEIcPaKAy5/GCdbbVPK90DO7caZuTSVfF2JkRukMelfEmPKMajj7uncAq6bdG5BcgJqUhw1luGbyomDoNCkwbzCdNS021Fl1k9qXxNN+wKAXiWFhOE+UwG5199LMWbxuTYWqd7Yq9VQxuehSBYo5Fdzkqo6tVyKoD8MnVJGm6Qmcg7kHvBrPD5pHsD1lVkJDQWBZZhzx2M1WqcCp9qHcLRpgHJkZOmUovdsw+wcNPQ4aRCTj/E8O0IpuP3DJtrPEmvMr9aE0tQ7IvicKCbj+fzpZBciLPfJpxmMxUQN0uc6ZhG7tpvK9WiifbmjxUFuf4gO5F3HOwTrS/Lb989bKYHrc3sbUdNuw3N7G0dfgOiAtntD2HFwYkrc/CpSfnFYy4AbLl8wIU0dAdEq0Smnfig8f6AJrYNuSqbDL8AiKvOAsPs0FqTwil+AlejZEErBJ26qmlC9xGSgUUgEnxPFRGs9AOA7N5RBIWXwnRvG9hAnGhP53HkWfKtYadYiP12FSjPHCMavdHx0fRV2HGxCS/8I7l9VCoN6fA8gkYXnFwRpldzty9QqBMU45LeEGer5A01QRt/a1v4R5Kap4o4XV0g2AcQWh+lVcuhUcpo60/M6SGPROjj1XBmblxXhvGUY3mim5myXPU5SEBg1BG0DI3GUBGIFWImeDVnqES9uRUXGValMFYNKIQPgj36OH3zPc6Kez2Prp6dE7jvRpdjnLnXKT5f+/tN2dDv8qLtixz9fX4Z/338Jdk8Q/mAEoQiLXcc70Ov046LFgfpuByrNupTTg4n0IQmT1Q82TEd9txMNPS54gmGwLEsFmslvb5+fh6GRAB6/sQovHWwBCy7zwbLAPctK0D7ggjvAvUhZeiVWV2WNO7XL79ycX2SkHZDFJjWe/dJszMjT4+jlAYQjLG6dm4ueYR98CdSJH1s/Ne75++et2LLrM3xyuV9Q3RhmAX8wDF8ghECYRV3XMN46dQVnLZyhvGLzoKbdhl3HOmDUyDCv0Ci43yQdnujZEM/ljRNXJt2xyQdRKC/O0CTsXn3nTBd6nX4EgmEsKjWNu8v1qztPoqHHhZq2IWyJFu6R+zjZtHYsYjuMU9U6/VwvQwiHAfkMRiL0c9OCAuQZOB6I7mEvrMM+PPteQ8r7TtQoxS+Q4s/kselQpzeI2o5hbL9jLhp7HBDM5ww3Iw7wmKUae1wJ3cJn3ruAsqf24pn3LsT9jX+ehSYN3rh/KR5dV4k37l+KQpMGt87JR5pGjjDLUeTbRdikCEqe3Cv478WDlwRaGIka9bpsHpxqH8KXXv4UX3r5Ezz8Zh2swz50O+PjZo29LowERo2V2HcA4Nn3G3GqfQjHWoaw/Y652HWsDSVP7sW33/gMfzvTiVcOt+BvZ4RNY8TjUMolkEuAixY7Sp7ci8qn96asoyKGWEEkMVzu456/OxjBltdOjXtJPDDiF3wCQlmHawGfO8+CP5Ne7HbAOuzD4lIjlpZl4LdH2xAIs/CHwvis3QbrsBfBCEvFemflG2CxeyGVICm/hVhhTJfNg2ffq8evP25BWZYWoUgYNW02ZOsUqDDrcPhyP945Y8G29y6itnMYdV3D2He+G+9fEAYZ11Vn48aZOTh8uQ+9Tj+qc3X4+pJiOlPyj32514kfvlMPFkBdlwNGjQybf1sDrVKCDK0y7hzFODb43KR9Dj8GRwKQMMkV0mo77DgfVXy7aHXAH4oICq/4hVg/ea8BHTYvep2jgz0/XQmXb9RMZmrlqDLroZAylBgnz6AQXRpFWODvZy043mbD2Q4bPm3lmqBa+t042WZHmAVOttkF3hDxOP7R2A93MEL1YiIscPhSHx7i8VyMB2Jcq7E43WFDxxBnGIIRFi5vELuOd6RUXNdl86Cpz4UBlw8/2DCNFk4RTpN0jRw/23cJOqWE/m0sEI9EJgE+vjRAx8ep9iHc//pp/OnkFVSZdVSh/sM3X02pKOtzzWfBj1/ENmfxQdStDEoZnP4QZYKel2+AUiET5QwgvAmEkRsAVv37ISpoJJdwzFGeIEv3lwo2XZdPpQgSRaX5/Bsn24YoE7QUo8xVeqUU960sGzdPB58v4cnd5wWzfCzyDAqkaZVo7HGhOlePnZsXCs6Xf/6xvBpSBjj8/Ruw42ATral4dJ14g9QTb9dhd60Vmxbko77biUZew1MyiEkgLHz2wKj8Hw9XU/E8EbcIX0g6EcjzJhq6/Osh4xBASgzlXTYPNu44Cpc/TMf96spM/OTLs3H3zhOUXX1JqQlLyzKw41AzLL/+Zk/IOZBU5flzF+Dkg58G27SgIKmgL9HmJC/2easTESRXCSMvBV/5LBjhZhAASNJhLtyX3RvXexILfgBqRUUG7nz1JHcccAIypOhq7fQcwfdTAf+Yt87Jx7ffOB3XiUrQ7QzgzsWFaO0fEZXR4++r2KTCFZsPKhngC3EygGI1FWLgF5ads9jxxFvnE8ZVYrHjYJPgGO+ft4oaiquNQpMGRo0Mdo+QLyAVqr1NCwropDBWi3tRCsHa3bUWqnVKxvnx1kG89ilX3WlUy5GTpqIpfAD/fzSS8XU9exyJo/mZOjnUcgnMegXy05RQyxlIGeD6qgzolVLMC8KDpwAAIABJREFUyNPFrcsXlKRjXkEaep1e3P3qcew41CzodNcpJMjWc66+QiauPiaGjbM5xuuxCF8KTRpUmrXY9Mox/OrjFsHvM3RKGNRyZOiUKRPQjIWeJETALx5sRSDM4sWDrbj+Fx9jx6FmPLeXi/3wYymksIqw+B9rHqLP58WDl+I0QPlr+h6HFyfbhtDj8GLvuR74QpGU9YBPxehvPLn7fIq/FMdk0u/3Li+ekEZ1oUmDzcuKBJk9gg2zc+jYGkhBZnJFRQaWlJqwfsaoDEg4MtrRu646G/sfW4XFpRn0vjNTxWfBMMxrDMP0MwxzkbdtG8MwVoZh6qL/beT97SmGYVoYhrnMMMxNvO03R7e1MAzzZConlwz8ANBYwaeBkSC8wQj6XAFYHX54gyzCLLf+dfnD2H3aSgNYZB3+26PtqLM4sP9iHy1X5i/YAmEWtmiA0hcauz2cj33R+AVZNj1/oEk05/2jdy6izxXAkSah0X/2vYZxB2vHQv04XH6CM53D6LJ5sOtYJyIssOtYZ1yvjC8Uoc/n1x+3xQUI+bl+fgCRiAKnujj2BIXxDm+q68EEmEwf0dunraLnnYrhScT0daxliI6tRaXJ4xWkP4oYbwISjK/vHhWp3rqnHjXtNsgMmUVi+4pFKsuQ34OTJPxDzPYXWJbdzt/AMMwMAHcDmAkgD8BBhmGIT/UrAF8AJ5L8GcMw77IsO6kR/9TG6YJ+i7t3nqRiN8UmNYxaOeq6nJhbaABYBrlpKlobwZci3DjbjH0X+vD4jVW0MKjCrMP8IiNcviC67B5ctDowK9+AmvZhSMAVwJxotaGm3Y5MrRzpGgWc3iAGRwICQeQ8g4JG/DO0ctoDQpYOibg2EumdkgKn64qNUyJx99DaMrx4cGxav7hzCIUFRWQKKZCXphJoZ8zI1eP+VWX42b5LuGNhPk602gXt17FLLXIdW/dcpCQ4MgAhcGt1mUQChy8EvVIKBgw8wRA0chm23TZDcG6PrCvHrz9uw7rqLHxU309ftA2zUlMjn0w5+FMbp2PbngbY3MIxELtUEkMiZvpNCwqoHgqfCyUR+Off0ONATbsdi0qNuG9lKRp6nKhpt1EiplGGthTEdZGCsWBZ9ijDMCWp7AzA7QD+wrKsH0A7wzAtABZH/9bCsmwbd2rMX6LfHdNYJNINITGEFRUZgr6BRSVGnLlix+3z89DS5wZYFlqFBPOL0tHQPYItK0vjaPWay9zI0I1S7fObw/i0+bu2LI57ISuyrbDYvQKiVk6LoQ61HXYsKDFi87JivHCgGQ5vENtun0H3wV/v56apBfSAz+1twEjUbSSZHAJSYHW+axjffO0UPMEI7n/9NPZ+9/oJLUeWl2fh8KVB5KapkJOmRq/TI1CFF4PdG8LRpgH670A4XmSnY9CDved6EAqzqMjWxxUL8a+/x+Hlyrn31AvKpRkJ8PLd8ziKfJ8ffzrZhU0L83Gy1Y7GXhfSNXK8cKAZv/q4FQ/dUB7V2yjCV64rwmuftmPzspJx83VMphyc3/9T8uReul2s4I7owuiUMnTavHhobRmVq3zmXU4HlZz7eFoU+Oe/eVkxLHYvNi8rpo1jfFpDPREPZyPilXoxmEyA82GGYe4BcBrA4yzL2gHkAzjJ+44lug0AumK2L0l2gES6IXz+QL5mZE0758L+95lugZzen0520TZwfiBTrNLyrw8uExxrLAYjvutIBgmnxWCn52Oxe2kE+tn3GkTZn2MzOfx+jkTuOFkaAYDTF5ow9+PzB5pQF02Pjgc2z9jr52HfqPDxz/ZdGpP1+vkDTQLaQYJghLtnfa4AtAoJAmEWu09baaqVL7z0yJt1gntFMggGtfyqNxeKoSJLg5YBzoA6vIG43h+iC0NoB18+1IrH1k+nBLtTce78lvVb5+Sj0KTBS1+7jk5MN77ALXkYmWJsXscoJhrgfAVAOYB5AHoAPB/dLubOsGNsjwPDMA8wDHOaYZjTrG9EtE6e1M8/fuNoD8CmBQWYV2CABMBNM7MxryAN1Tk6yCUMbpuXK8pCRPbzwKpSmPUKmLTyuPXlWAxGJCi1ednokm/TggKsrsyEWi7BmqpMPLVxOozRmgcxfgX+eZDr2DAzqUYtgNGHl6aSTYqhaV5BGjbMNGPL8pKU3fWH145NKiOTABtmmmHWK2l7+VjnUJ2jR3VOfMcmuWeLSky07yNTx8XjpLxRxYKbKck93LK8BFuWl1y1/pBkgdDNy0b1Srsd/rgYyFMbp8OsV0AVjV6aoy0AU3nuuekqwWcswTW5f2wknFL6aEKeBcuydOpjGOa/ALwf/acFAL++tgAAYQNJtD123zsB7AS4Ogsxy8p3tYRLAwYRAD3Dfvz94ZX40sufIhhh0dLnjtMf5e/nhY+aqPpXldkgsOZi2qMEzX3uOHo3AJhXZMRPvjybziTZehWeP9CELSvE15yxru8rmxcK3Fg+zHoFSjK18AcjWDM9E2+ftuKpjdMnnBEh1HdkIKVp9EmXIQBgjNEOJTl9gll5eryyOakUBT0Hov05/V/3gq+LrFVKKW9naaYOa6fnIEOnxM/2XcKcAgM+auiHTiEDI2GwtGyUY5K/lLwaSNYdupNHQGPWK+Je/Ng+Fj6fyVR1RhNOTvLJ753KN2pojQ0jkSWXvMMEjQXDMLksy/ZE//llACRT8i6ANxmG+SW4AGclgFPgPItKhmFKAVjBBUG/NpFjj4W0KA8l+UyVDIQfRBqPNY9tZAPEBxGfwSvVQGR5plq0uavfFUAgFIFBLcfOI+3wBCN4+zPLhMVtiJFwRisOU00Bb323UfDv2BxEnWX8WRYAiBVQ313bTZeJAJcCzDdqYB32wekNIhThljzZOjk19gCumgYtQbJA6FMbp1OB5z5XAE+8fRZ/fXBF3PfExLWnCrHs86Ql4ZPmQbz0tfn41T+aowHgSEoppKTGgmGYPwNYAyCTYRgLgK0A1jAMMw+c99cB4EEAYFm2nmGYt8AFLkMAHmJZNhzdz8MAPgRXiPgay7L1qZzgeEBYssknv9hnLEzUmovFLMQG0UQi7IeeWEv/f8XPD9GXhQUXYOQ3aJ2+MnFiWmLctiwvwaPrKrG7toseazKYSL0BgDiWcik4rlC5hOG8xP4RPBot8TZpZdRo9Y8E45asV7NFPVkg9NY5+QI1+Jr2YVHOkqsprs0P1p9qH4I/HIFRLaMCTcUZGrQOesCGw/7ke/ucl3vHYipSiePB/9TxyHF6HR5csXEvsl4hhUzKwO4NYV6hAS99dcGkyXP42R+nL4DGnokTuWxZUTQmAUsiTP/X/YLqzfw0JQpMWpreXlZupEuvuQVGbHrlGPpcAcwr0OPvD6+a8Pkmw0QIYza8eBiNvaNLVLGy/ET7nWo5xrtePYGadhsVkd6yvARvnLyCYIRFz+8fg7+3Oal9/1xXcMZiLMr8ySBRMOtqHQ/gBsuSn36Eb79xGv0uzkCUZIwup1ZWZkIq4R5fXZdzwmTEsVWg5JqeuW0WqnP00Cm49NqSUiPy01Ja2gIATrbYkn9JBJEYj7jApMHSsgysnZ6Dvz64DH+u6aJFaa992o4+VwBblpdcVUMBTKxYa/9jawT/FvN0ElXhTkXHKX/cbl5WBLNegepcPbYs51LK31gWDSOm6AZe070hDm8QK35+KGXryudgtLlDVFx242wz3jnTjQqzTtDMxK/XEBOMJeC7iqRpbdOCAvQ4vHGeRSxtHABKwvP0LdWihoX8xuUL4mhTP/p5rev7L/bho/q+uArRjy/3w8fbaNJO7FGKeRZPvHVOkJYERtPSqaKxbyTpbCz299i+ts/a7ahpt2PHoWaoZIBOya2/ZRIGDd1curdtYERAVAwAv/igCQuLRwPNk+G2BIRFdGK1P4mglktoVSm/0Yx4XonOi19wGItUPVo+JaHLF4ymagdRmqnjpCG9XAqaDQVSYj66po1Fj8OHYNS6pmIsYjkYSWCsoZsjpa2zOJCfPrrWFKvXSNbYxTcc5Hf8GozYPDkA7D7DcSkkaljjc2eKQayU3BezcefR9gkR7saumZ8/0BRnKCaKZPyZqfBr8m2HLwT4QpwhtTr8GI7GbS71unCm0w6XP4xffHAZvkAEYXBarqT+ZLKxAX7mbDz7ef2+xbQRkI9dxzpRkqGFzR0ak/VLDBNRLxOTvhwvO9s1bSyydArolNI4tmbiQfDLiHPT1MjSKdA97EO2QYlypRwVZi3erevBdcXpONVuR7ZeiRUVGaMZAJ8fCikDk1aG/HQVTFpZ3Kxxqn0IW/dwsdgVFRlYUZGBk21DWFGRAZNWhoZuB20OI995/1w3bO4AmvqcWFpmgkEpQyAcEXyPD0IrF1utCXA1FI7YFAG4OgZ+Y2aiWoZkM+qKigy8c8aClw814+V/NKMqppFpMnj7NFeHt+8ClyXfMDuHltNvmJ1Dq0DbBkZo1aJOIRG0zYvdE7LNHeBmxj6XHwopoJYzmJWfJvCCmvocqPrhPnx9aSEeXVcJk1Y2Lm81FrHB6mT3d6yZn59R2nGoGTsONWPDrGw8vXEmXvu0HW2DIzjf5YBUwl3rkxum4d7lZdi8rAgWu0dQ38PHM+9dwOvHOmHUyLg4AxPB/OI0Wnj3+vEONPQM8wvhUgpcXtPGIhRh4fKHYXMLXxZ+gxKpzFxaloG6KD1ba7RyrrnfhUCYxbFmrhnnis1LBYh2HGqmrNMf1vcjwnKzM8kEECvPry4kxyFp0COX++Hyh/HOmW46qx9rGULrIBfU2n+xD+ctDtoav+9Cn+jsT/glxZ6YmKHg7o3w34kGfrIZ9VjL0KgnwUIQkJssRujLHMCu4x20NwEAGnqc9D5d6nXhSLN4l7TYPRHbFggDS4qMmJFrEBgLUjOy61gnOn5+C+Y/8yHs3hB+9M7FCRmL2CxIKt6RhIGgongs7L/YD7CMKCv7c/sace/ysoT1PQS7jnGVmUPRdvldxzqh5REJR8BlZwgYmSKlddk1bSyMGgW2iFRwkvVcrGdxsnUQtZ12pKvlGHQHUZmtQ+uAG0qZBIPuIIwxgkAME8HvPunApoX5aOgeweZlRdH+gtHvcM1N9fT/ST5804ICWIc9cXUcmxYUwGr3UMKXNdOyuN4QXzDhTPD0LdV46m8XYPMEoFdK0WX3IUMrx6JSE3IMahy61ItO28TSmcnStisqMvDrfzQjGB3MfJHfyaLIqEa/y4+5hWmYkZsW51m89VkXJVQm2rFjLceSgTyfsfZB9JYS6C5NCWK9je+uK8d/HmyNq0VJhFhFPQIZw73wE0nFZ2gVcAcmlxK/po2FQiYRnQ2J8Mzy8gywrAS5aZw4jdMXQigCBKLreZcvDJc/DKNWAbiDyIiWCfNnh8fWT6cPd26BMW62yU1TY2aeAWc7h7H13Xo8c9tM+ts7Fxaiy+bFnTGksHqVHEvLMmiAUyGTYHAkgMOXBgQNZ2RALS7NwKEn1gi2M0wErx5pg0mjQL8rpTS4KJLVAxxrGaKGAkjdUMwr0CctvPKHOE2UZeUmnGi1oyRTjYYeJzYvK8KxliGsmZ6Fs53D+NHf61FsVKPHOfHrBEDHwViYX5SOI02DmF+ULvp3sWVFoqVGl80Dlz/IlcvzNE9ivblY1u1E1bkEoWhGSCERBnxz09VUVqIiO7nqGB9W++RrZ65pY0G6TldUZOCt011o6RvBl6/Lw0uHWjDoDuJ3n7RhJBBGU68Tr2xeCHeAc7si0doRQuBLtD5aBtyChqvY6kUg3lXfXWuhAUoAeG5vI1ZPy8amBQUCvQuyNuUHOE+2cTMEWcb89xkrwLA41NAPp5fTx3z7s068+NX5NBtDBhpZIvEbxqYaXTYP6jrtNP5hNihg0irQ2DMypussYUaXemOBNEmRFniyBKlpj0+rkm7ayWDHwSaAEcY5lFIG/vCorGBZpg5HmgZRlhn/snEdw+fo+W1aUIB/fecCjrcOUVa02IbG3bXc2OBX5/LjWpNBIDIq3yABkG9U45E3z6LP5eeC9UbxiYAE9kkJ/rxCAwZcAUGxHf/5XtXekP8pDI74seNQM945Y6Hran53JFkTH2kaQJfNg+7od9yBMJaUGHEqSqQyxKNZO9dlp5V0pFZ+w0xzQmEXsqw42zkMhVwCpVyCHYea0dTrhC3Ku3Cgvpe2E6+oyKBGLLaTMoJ44dxup59Gy49c6sP9q8uQn66ifQ8ZWgXc/hBGApG46kYCfvOXWFotUeEVWBY10XukkTOoefoLtHU6N12FnmEfSjI1qGkfRnmmBnqlDPU9Ljy0tgxGjYIG6DK18v8rVHaxONs5TOMgarkEcwrScKnHBY0CePZLswDE6+PywXUM21Cdo8f757rx2idtcEXHmBTxbr9JK4NWIcXCYiP926n2IfzLXzl280Tl/aurMsak21PKpbhrURFcviDcgRBOtg7B7g1RIiSFlIvRWIdHYxb8Z/zCXfPwnTdqMegOojpHj5e+ugA9Di+e+tsFSpWQrVdxy2sG6Bzua0vl/l7TxoKAT48uBk8wgt21FvoisQDOdg2LBsION3F8hFtvm0lTRx1DnoRNT7E8kkt+ehAAxxbl9JE0no96LMdahqgRGy/qrE5aQg5ws/3di4tx5HI/6iwOKOVShKL7JuXPAAQBYLG02nN7G7C/vg91nXbMKzLiQH0vGntdgoo8T3Qt8scTnbR1+tF1lXjtU24c9bv8OPTEDYLzvXd5mSCV+H8bNk+AvkhZOgUuWh1Ul6W2Y5i2aSdalpEX/sjl/jhDr1FI45Y4O4+2wx0Io2VghP7tub2NsA77RJvHCGIJfEmdTUO3Axe7nXhsfYXg3hpjGOjJ8Npd242W/hG89NUFcUufAqMGg24HlDIJPbdb5+YJllKkeY95zJNSqe41XcGZppZjSakJt8yJ18DkQy7hXD8i1ScB8PWlhdBENyRqjpqZbxB8poKvLimAQsrgq0sKkBdtKzaqR1vEV1RkID9NybVoz8qGTiEd100mLfFPbZxOvZ0vX5cHvVKKDbPN0CmkkDLA+hlZyNQqUJ6ppR2Lrx9vQ23UhVbJGdz16gm8fryNRtaPNA9ix6FmtA9yY0Ms4LZxthlahRSrqzKxaUEBnripCnqllBY7xYK0VBcZuXtRnTN1qVeAc5eTIS9NCRnDLWXIi9Rp91FDkSqIIblfRFh5hsgYIXKJ/LQ1aQe/rtiYcvEXWbr6QxG4A2F0RBsIyb0NxTBC85Xa67qceO3Tdji9QUFb+/2rSpGfrqLXIlaBSio8kSIH5zXtWYz4Q6hpt8FiH5X7E5PgC0eigTqiEwHAoFLCqFXAM+yjL4VWIcWdCwtp4DHVRjM+TrTaEQizONFqpxoZgTBLB8axliEaZ7C5Q+P2MsRa4vdd6IPLH8aJVhvd33mLE4PuAL6+tJi6uts/bALxMYibS6oc+Ygt6Io9ljsQhi8YQaFJg3uXl41Z7EWa8Bp6nOi0++D0TVUuhUMqKcdAiKVeJfEsANA+CLNeKSq6nAgkJclvu/eLcHtm61UoMGqQrR8dn3xF+FRBxl9dF7ckJHU3hSYNDGo5XQoRmkGwwKYFefhbbTfKsjT4W60FTn8IG2aa6Tgk6dW953poVTOZfGJjdVJNWmYq53lNiwz9/ne/3fa/HnwQt8/PQ7fdi5WVmbj/+lKcarMJiFqn5+jwgw3VaLA60Of0YVGJEY+tr4JJK0dthx1fnJsLhzeIZ26fiQdWlY8pyJMMfMGZaTk61HbY8f2bq6j4S4FRDbDA/MJ03LO8GDqlDBXZWgy7/Ql1Rgmqc3X4xtKShMd8aG05NHIpFFIJHl1ficpsvUDCT6uU4Fi0XuErC/KgV8lx7/Ji1HUOwxeKYHVVBtZOM0MpYzDg8oFhR1+GJaXpuGNhUUqCOmPdl7sXF+Dw5ZSY5ZPCqJFjbkEabCN+QV2JTiHFuuosXBl0Y1GJESatAtZhH4pMahi1Stg9XJp8Zr4Bgy4/7N4g8tLVScV+CMgzZFkWvU4/8tNV+OmXZ8fdj++9xQVD+TKJfBGsVMcVGYcvH2qB3RPE0EgAt8zJhdMbxMm2Ibi8Qdg9QfqsfKEIvjAjB0vLjdh/sR/+qOfh8oXwrVVl9BoMKjk+auhDj9OP5r4R/Ne9i5CmllNS6vmF6VhdlYUPfv9CRyoiQ9e0Z0FSp102D+YXjRaRZOmVSFPL0eP0whOIwO0Po9CkweM3TaPBvUKTBidbbXD5w+ga8sbNABNFLP/AjLw0VOeOiv3GtrsTfs0H15TjPz5oEvU0DEopghEWdy8S1+Ukx5xbYKREOtl6VVya997lZZQYZs20bNw6Jx9dNg91a0kQlt+/UpKpxs6jHF8lMDbZjxj4/TVLyzJwwTJ5gWcCuycIlUwKb4wnNBII48OL/YgAqLc6YdJyXrRcIoEjSvc3EgjSgKBOIRmX99jj8GLfhW70uQIw6xV44a55ooFKMZLdyXB4KqLLZrs3KJj5N12Xj0y9kmZptNHrWff8YcHvf/TF6rjzqDRr43pM+HUahSYNvhcOJtcYwDVuLAj46ch9FxQ0JUdgcwdwqn0Im393Ev4Q8OAfavHqPQvwj8tc9d6pDjtYAFv31NOgDjCxtmN+xd7ZTjvqLA48t7dRwMTF3y/5fqJMBgA4ox7H8x82ibr8qRyTHPf7b5+DN8hi25563DonX9B34vIFsf3OeaK9KE/uPj+hisbY/pp5BfEq6ZPBsRZxL4XMsiOBMCIst+xz+YM0K8NXCBhLdU0Mj/35LB1jfa4A7vndKWy/c07c/RmvYU2Gm2Zmo7nXhQydAtZhD9zR6l23P4RCkxqDI1rOay3makTmF6XRSszH1pcLzo/fWHfsyXWC40zUoF3TAU5+ncXqykxoFVIoZdwpm/UKyHnBr+f2NiJaTgG7Nxj9N/d2kou0DnsFbebJ2o4J0zLRWYhFLMchAXHzHnnzDNoGuGBiKroiAV5ZIf/Y/O7KRMck1+ONZjWI5AEfYzUOkRdqvCI7hD9087IiLCk1iQYHJwrpGAaWDxJPsLsTT5DjaS3viymCIxooVxtvn7YixHIGanetlbYBdAx5sLvWitYBN/QqGXbXWrG71oKzndy4UEiZOPZ00ky5/cPEejrjDXBe08aCsHsfaxlCy8AI3IEwLNFKtMGRAE2NpqllgpJro0Yu+Lcq6t4RFmwCPlGuGIhHs+t4B/3dfSu5QNF9K0uTBrPqLA5cEmGtTgSdSor3z1vxwkdNdPbfdbwDF6Mv+afNgzh6mWu+OnPFjg0vHsUTb9fRF3vTggJKeruo1ETPl2wjWZ/7VpZiy/ISAUEusbsT4W1weoN44UAzatpt2HuuJ/kPUoROKUORMXmgUBodxQw4UiACVXRiqcjSjmsZsrBESKwslyRu1IvFZBTNnto4nfZwzCtIw9O3VGPL8hJabAgALn+IZj2+c0MZFFIG37kh3htNlsUCRp91qgHOa3oZolPKMKfUhEqzFlk6LoillDNUUYwgGGFx58JCnGgZRJ/Lj3uXF2F5eRYONfTD7g0iEjUrRSa1YNAkc8fEuDmJ5NuKigxsmJ2Dhh6nINLeZfOg1+FFplaBmXkGgAHCERYahQSddp9oFyXB4EiQNskRhmeAiyv8dO8lBMIs3NFZlNRCNPa6aNt9oUmDnfeMiheTaxTbtvW2mYKlW1mmRnCdqb5c/CUiwNWfTBUcvlDCRjo+Ck0cPVwgAhjkErgCYWiiPBJ6pRSPfaFyXBwWj984Dff//jRtAJRKJKLxLjE2Kz4pbio6r7HQKeVYWKxHWZYOuWlquHxBdNpGKQNaB9xweILYMDsnroycj2RZLGC00rTT70kp0HRNZ0Oe++VL23rzVqLb7qUdpWJ6uYFQGOctDrQNeRBhgdoOO85bHOga9gp+EwxFcKrdBp1Kir3nezHk9nF9CRka0cj/5V4n3j3XjZtnm/Fp8xAKjGo8/tY51FkcuGh1YMDlx0eN/bhodWB1VRbS1HI8+1493j3fA08wjH6XDy0DbniC4ZQGvVwC/OTLs1CZrcfa6my09I/A5Qui3urE4EiAc8ujuUStQoKKLB3WTMvCfStKaeT9cq8Tb5+2YHZBGvKNXLXm42+dQ+vgCIoy1Nh7vhdSCfDbT9rxywOX4fKTkniO/DYYCWPfhV7B77/31rmE94hkDtJUMviCYXx1SQFOtk2cE3Qi8AXD9L54osaUFKwFwiyGRgI0W5EK3jptwVFeF2wowuJoUz9unpUryHBs/m0Nepx+1HbY8K1VnDTCrz9uQa/TD4VUgrsXp6QKKNhfnyuAwREfTnXYYVDJcbZzNEVP4AmG8XatBb892gatUoI9Z7tx5PIAyrN0KWVgSID7ZOsQPmrsh6tuf+THT31/e7LfXdOeBUGaRg6tQpKwyEYhlQgq7oIRiIrWeIIR1FkcVLiGlCmLBQuB0WpIi91DqyqJ5FsgHKExgMZeF63g5JOMeIMRKGUMjZ0kQ1GGls5QL3zUNGb35KqqLLzyjfiq0+f2NgoCoPwWe+K1xIozcdcDUcGlrXvq0djrigsOExAv5Zl363GkeRAfpiAjMJWQAJhTwHFYVOfo4fYHBWr3AGCxe3CqPXVm9UqzFgalDJ5ACGC4yYaLIwiFnK4rNmL/xT6BHszTt1SLyhCmgpJMLfpcAVRm67B6mhmbFhSg0qzFj96px4gvKGj4A7jgLolNAMCBht6EmRs+yBK3PJMroGMkEumYP4jimo5ZZOqUeHRdJfqd/jGr8XQqWZxATaZ2NGbDj/WZ9Ur6cNOj2heJpALEqilJt+L8onSURF33ct6a+OlbqgXHLhxHrUIBb31Oqvf411VsUtMHTLQgYkGuhXxSAZ9cPb2Ox2+swpblJcjUjs5CCilHKPvUxulCUSUSzEiVrnuitN4JICY8xEcEAFjuuZZkatDniu+uJMYvVbxwoBlOfwghNtpgp1di03X5cUuzpzfOwKPrKvGvllmZAAAgAElEQVT0xhkJ9jQ+jBZ+MbSp8PClAdi9QRRlaOMyTTqFFCsrMyCXMNAppILrTCV2QlK1bOTqyxdedZA6i6Y+p6inQOD2h/HS12bi/t9/Bqc/jDyDAlU5BhyO5tlJlkSjkOD6qkzcubAQVWYDTFoZdh5tx5ppWaJp1MWlGdh+x1xakwAI29J3RbdXZevob3LT1JiZb8DJtiHkGFS4dW4O3jjRBbsnkJRDobnXJahbMKjleGhtOf7t3Qb0jwTg9AURjHopn3XY8O0/nkaPwyfg9pydb8CHF3sxOxrMXFyaQWMW/Bb8xaUZSNNIaUcoy3JdsgwTQcfgCLbuuYjWgREw0S7Oxh4Xyp7ciwgAlYxBaaYWrQNuKKQSzMw34GJURb1xAqrsY8HhSd51e6bTjmCEIxsi8p2xsaFELGWix/SNNsXJJFx2RK+Sx8U9xGJexLPbuqceN87MSSktT575mumZaB3gOqtJ8JEYS9Igl6mTw+YO4t7lHIfn7K0fIBhhIWNZgZGP7RXhj2/STEc4ajt9IymtG69pY0GQrHRWLpPgWMsQrVfodgbQ7YzPz3sCEeyutdKA4F2vnqC6H81l7oS6qrGcmoQpi6S2yCf5PikGumLz4u3TVqoKngykGIdft7Ck1IT+Ee73ds9o3OOKzYsr0cAXv2ls59F2uPxhASdnIras333SQf8/GOGu62wnV84eW8sCgFdByKIxKmAcCIfHTeY7HhAF+rEgYUZNA8NwRdqxdjkRS5kYtt02A9/9cx0irHiMbCwQYSswoN3Jn3XYkaaW42dfmS26ROBr7rr8Yey70Iftd8wFAFiHPYKJcs20bGy/YzRw+sRNVbSegn99sYFqEng9UN+LnfcspHyiNe02SJSalJqjrukA58u//s02X/kNWFudDZZl4Q9GoJAy8AXDVEBVwgAPri7FiopMvHPGgkh0u1zCgGG5ISQDN9B1CgnKs3S4ebYZe8/3Yk6BAc19I3hobTnWTjfDoJILyqcBbmng8gahkEoQCIfwRk0XpAyglElQYFKjoceFAacfLxxsxosHm3GKFwuQAMjSK+DycqW6EoZr8JIwQqYmjZwBC+A7N5RhRUUmOm0ebJxtxrA7CJ1KhgGnFyKtCRRufwC3zMlDmloOs0GJc13DeGrjdKrORUp/Y69Nq5TQ0myFlDsnrUKKQJil//68gAWQoZHjn68vQUvfCIxaBTRyCXyBMCQMV/6+tCylDCGqzAaUZ2txqs2GLL0Cy8syoFPJRQOIJFhYYFQjTS1Hda4BBpUc9ywrRl66Gh9c7MWgOwC7JygoC+f//qOGXrAsUGTSIF0txw82TMeMvDQsK89Ada4BYLn075LSDHx7TYXgHOYVGvGdGypouwFBbCvDT95vgMMXwuBIAAYV97chtw/nuoZh+fhPA1t/+NR/JrsvqSiSFQL4A4AccO/cTpZldzAMYwLwVwAl4FTJ7mRZ1s4wDANgB4CNADwAvsmy7Jnovu4F8K/RXf+EZdnXxzo2X0WdpKH4LdF50WDciVY7LDafoEU9yOtAIvPxSCCCxl4XVYrK1inQz2Ow4s+6fLdNr5JHlca5v4VZrtWdLM/5Cz7++xWBkCQmwgIPrKqIc13JNbGshEod+oPhlJXNh71hGnwTY4VOlCK+d3kZZZcmVejEO5tgl/2YkEswptFLBbGEvsCoURvyBGGx+VBg1KDO4kB5phZDIwGEWODPNV2iacZEVbz8+0iej5iqeazXFqvD29TnxP6LfTCq5aJBTz65UmOvC4+uqxR4H1OlfUp6qdRyCZUzOHK5nyPJ0aSlpIadyjIkBOBxlmXPMAyjB1DLMMxHAL4J4BDLsj9nGOZJAE8C+AGADeA0TisBLAGnuL4kaly2AlgI7p2qZRjmXZZlE/qwYhyc/NqHBSXp+OOJTjx+YxWe29uYaDcJQdx7MZp0EjG22j04fDk+ws/vbhwPKs1agU4IAIBBnGr2n2s6xrXfv5y6gl3H2qGQMnD5Q1hWloFnvzQbPQ4vJTl55raZyE1TY8fBJgFHaKZWIbpUUsslWFxqFCVqGat8PRGIoYgVUR4PkpVut/SPUCNL1vlAYgX7sdb2xHiI1Z4QAqHNy4riujn5v3164wxUmQ0JYxdkPLt8QehV4qQ8U4Ftt82gNSHHWoaw41AzVldlQq+UIuJzp6QIldRYRAWQe6L/72IYphFAPoDbwWmgAsDrAA6DMxa3A/gDy+kinmQYJp1hmNzodz9iWdYGAFGDczOAPyc6tjcYxu7aLlSatfRG8y1tl82D5jI3ctPUuH9VqUBbMhVwxX4M1kxP7J629I+IskDFGopUZ83fHm3H6mnuuLToo+sqaZMXAIHQEB9qGRPXWAUgLsZwuInTzDjZNhTHTs6fybps3oQxlbsXFeG+laXUsLkDIfQMcwHV3DQ1HvlTLeqso/U8fCNQnskVSomBfIfwjMReDgOO1t8TmyvE6H02quUY8YcQjLCQSzivzaxX4v5Vpdi2pwGD7gA9n+pcPZ7eOEO0iIpPgRdLq0c8BDHPjKTVAdB4kZimCPltl82DJ96qo0aaL3ZlUMspbUKiXqXx9DGJfZfw1s4tMGJuAWc43z/XDZc/DKlanzXmDqMYV4CTYZgSAPMB1AAwEyV1lmV7GIYhrkw+gC7ezyzRbYm2xx7jAQAPAIA6p3xMkSH+rDARcC88i7dPW+NcVBIxXlGRgUfePIM+VyBOq4MPk1aJ6yszsaeuW7AEikVuukowmxA4vUH6kHccakaxSY0rNm9cVF/MUAAcAcyIPyzwLDYtKMCKigzqWRD2az77+J0LC7H/Qi9Otg2isXcEajmDxaUZKMvU4b6VpXEVmsCo0eEbCkDoLSQyFHwk8kxYQNRQAMCKikx4gxEUmtTYXWulXKUAR1HY3OfGmulZ2F1rRQTAklITtt8xF4UmDe7eeYJKApLxxFe4B0CDysmUx8Q6Tseqfo1dbpAgO38MW+2ehNWf//rOBRxpHkRdpx2v//OSBHd09FgkSE6uPdaD+pcvVOFAQ++Y+4lFysaCYRgdgL8BeIxlWScXmhD/qsg2doztwg0suxPATgConDmXzY5hIuKD/3DOWVKLyGdq5ZBKGDh9IdwyJwcnWm1j1v3npqnxoy9yLlyWTi5gtC7PVKNjyAuZhOsMrOsaxvIKk8Btj6XW90Tb6fneETebWQXckHy+zAf/cBp2bwgKKddNyQBQyBgEQyzC4Bi5XvnGItHzLzRp4oqp+AOxy+bBu3UWqjHhDbI40jSII02DaOhxYFm5CVqFFLPyDPCHwqjvdmFZuRFfua4AfzjeLlBzn2qQpY6EAf7pujycaLXhhulZ6Bj04mynDXIJg68vLcQfj3UiBG4wV5q1eP9cN4waOfRKKWrabdj0m2N46avXiUoCJlK9T6a3ItZx+rczndhxqFUwgeWnK/GXB5YLJCKyDUrsPt2FnUdbsbjUCIWUgXXYjT11HD/rBxd78ej6UcX1U+1DVFfl0wT6KgRdNg+cXo57s6bdhptfPIp/3zQblWYtMrVy7D3fjRUVHMVCSYYGthE/utz2lBp6UjIWDMPIwRmKP7Es+9/RzX0Mw+RGvYpcAGRhbwHAD/kWAOiObl8Ts/3wWMdNU8vj2mv5IC7eqfYhfP/tcwm/p1NIUZGti9Ma7bJ5ALYJW/9ej4ffrEN1rg6BYAStgx7qRTT1OiktHZ8dGQD0KgXCrBcqqQQjgTBaBtxoHRAKv8SGNWIbywhJ7JJSEzUQ5JoeefMMXL4QfSFDYc66sgAeXF0hWGOfah8aU681EXbXjhqKWNS02ymPZceQGzKpBMEI54kZNQo4oueVKOYxWRDPI8Ii2ucTwq8/bhUsud440UkD2CFwHKIkVkEYw/ucAZpeThb8rTRrcffOE7hjYT71MFLFrz+O5721Dvtp8JkY6btePQFrNN1OJpb/PtNNWcFGAqMBawCCgrJkK13iCZICLncgjJ/tuxTl5ORa+IlnuP8iN64ZuTolXYGkFZzR7MbvADSyLPtL3p/eBXBv9P/vBbCHt/0ehsNSAI7ocuVDADcyDGNkGMYI4MbotoQgLeqJqtBIldpzextpazbArZf5VZQjAS6zEFvFR1zDoeigauwZoe5zKErDLtYYpZYz2DDLTKskF5WYkKmVoyJLi7mFXLpSAiA/TYli02iNiJQREpQA3Cy26bp8NFgduP4XH6Pkyb2Y+eP92LqnHnUWhyBIFwG3Vi/P1GJFRQZOtQ/hlh1HseNQM55465xot2iySr6xXob8dBUqo/oUMgmDOxbmQx+Vk9z+YRMduJMxFGMNQF4DKc14TM81cGQw0epTfpxIIWWwcbaZVqwuKR3VBqlpt+He1+I1R2NBSuJ/90mHYHkCcDP8Xa+ewKl2cWbury+N7z/JT1fGLUs2zjZDIeE8ptVVGdArpbhp5mhCQikdXZaeah/C4IgfSinnmH9zxdj9JisqMrCk1IT5xWlQySQwKGV4auN0bF5WxHG2ZnGcrZsWFIybeyQVz2IFgM0ALjAMQyKITwP4OYC3GIb5ZwCdAO6I/m0fuLRpC7jU6RYAYFnWxjDMswA+i37v30iwMxH4qVMxV5C4ihtmmXFlyAN7NEvSOugR0K1LGa482qSVUxkAgJtFjGo5rbvPS1NgwBVAMMIF0lZWZmJpmQn/9l4jHax8l5/zTLgg6K+/sQCLS7kg2SNvnkGdxQFPIAKDWgKDSganL4QvXyee1uyyeynPIgC4AxFYh72oztGjJEOD8xYHnYns3iDs3iD2X+hFQ4+TpjoHRvzYsrwEdZ12LPnpQXznBi4tmkgThc/DmAg3zhjVI7U6/Dh8aZDKST5xUxXthJ0MxpopC01clahWIUV+ugrO3hG6cCVBZ37vTSDMYt+FPrrsuuvVE4L9kfEwlpDQA6tKsfNoOx5YVQqbOyR40ZMJEhtUSgCjmhwKKYNjT66P+94fjndS8aCLFidCERaf8JauLEDlHv3BMPVWH11XmZS0hsRgGrod8IUikEsZyq5GOFsJe9vTt1TjWMsQvrfdHq+VKIJUsiGfInHFf9waIZoFeSjBvl4D8FoqJwYkli8kINud0ReID37cIMwCTn8I+y/2ocpsoDf8t0fbBb8LR0ZnqmyDinPTWO480jVyPPdP8RV4J9qGBI1XhSYN5hcZUWdx0Bcb4NKQa6aJB50fv7EKD7x+GsO8zlSnLwS9L4hb5uaix+FFt4Nrb4/E/O6Jt85hYMSPJzdMg80domtb0mC0ZXmJKGcHMbRblpeIZi7UcoZrNsoaZev2hyM0xVto0qA6Nw1b362PaqRObZk3wHl28qgocFvUwzpvcUAapfxWyyXIT1ejZcANCbiGw83LiuiLv6zciFPtNhoYW101KgQVOwnx78emBYVYOz0nbjknFtTkgwSu3ztnxaA7iBm54t49X9piSERcKTdNjZw0NRdvWcBNLhVmXUppVfIdk1ZGxwAxcmSpS671yOV+zC8STymL4Zou904kX8ifGci6HQAONfai0+6DFMCXF+Th8KUBKGUSqBQyVJl1yDEI+SxIaa5RI4fdE8R1xUZ0DHoABnjohnI097nxl1OdGHQHMOgOxInG7K61wDrsQ366SjCA7ltZCpcviLOdw7AOe+ELReANRvDHE52iWZ3FpRmo23YTAFB5Osuwl0buSayEDHq1jIHLH0S/yweTVoHFpSasnc5xajT1OnGmc5h6FmPl98nn1ttm0rQfyS54g1yvAcDSma2xx8X1fjARbJiVh++8cQaD7gA2XZd/VYzFBauTenSZWgWsDj8ydAq0RM9HJZei0qxFy4AbEXAxitqOYTT3uWn5NLlnW5aXYOttM2kAMLauhT/xJPJmY4OafC5Twm9qUI8KLp23ukS9mEhE3J9KU8mQl67GM7fPpLP/ioqMccWi+DGYtdNz6D72X+iFyxfEa5+2Y8PsHNpOUGdxQKo1ptQ487ko9yaltASEto6UrTq9QRxtGkBzP8cdMbcgDc/cPhvfv3k6/vn6Mqydno0+px9fX1osuOGkNPeG6Zlo7hvB7fPykJeuwT3LivHWZ13YU2eF0xuERiHBuupsVJq1uP/3n+Gt012YnquHSi6hpdXuQAh3vXoSrx5uhTlNAblUii0rSxCMROAPRpBjUOHpW6qTMmbnGzW4e3ER5hWmU0bvIZcf/U4/1HIJgmEWoQjQ0ONCbYcNrYMeNPS4YFDJcdOsHNw6Nw/fWlWGeYXGMZnLY8uBSUn44tJ0XLA4sKI8Az//yhzsPNoGd0xRSV2XA502D5qjxWwKqSSOcwEADCoZ/NFcs1ouQbZOLmA4lzBAhlYOTzCC8ky1oPcFEKbQyO8c3iCkDOdh+YIRtPaP0vazAHINSly0OjDNrEe/ywdfMIzFpUb8YEM1Zbb+zdE2rK7Kwk2zRkmL0tRyFBjVONs5jPmF6ZhVYMCjf67Dn2quoMqsE31uZF91XcN0LHL7sKHP6cc3VxShdcAjGKsA8IsEVHf+UARfnJOHfpcfM/IMuGlWDt46bYn7faogzzjfqMHZzmHsOt6Buq5h5KWr8dj6Kpy9wnFlOE/vcafCZ3FNG4tnt//nttPqBQDLNdAQxPY6kGpLUtLa6/SjoduJpWWcIXnkzTPYc647bj/kZv5s3yU09rrQ3DeCjxr70Wnz4HDTIAJhFhFwtG03zszFq0daYfOEYPMEsafOioZuF67YvLjc48Kfa7rgDUbgC0Vw+PIgatptePecFResThSbNHjlGwswIy/1gFK+UYM7ot2xdywsxHfXVaI4Q0P7Fb5QbUZzn4unARJBcYYW979+Gr8/1o6957ux/cPLyE1X0h6RsUDuxdKyTDx0QyW+NJ+7twoZ4qj9JQB+eddc1LQNIRCO4FurSvFp02Bc/MHPK0oJ8choCLh6Cu47GoUMrhQIglgIl2JmgwJgWWRqOUNktXtgGfahx+GF0xemnsW/fGEagMR9MgAEhmTfhV7UWRwYHAmI9nSQfbm8QbAsoFJIUJ1r4Oo5FhfjsfVVWDPNLHo8i92Nhh4XlpSmw+4OIBhhUZ2rwxfn5AEAfnO0TWB8wHL3MhVym9heFf65EokKMmmuqMiEQSXHB7t+2ZaKFMA1bSx+9suXtylm3QhEIninrpuyNRE2KHOaEnvP92JFRQZYlkW/0w9PMAyzXonLfSMwqORo6HFyhgLcjeIbCwK+Lkdlth53LSrE27ysQn66Gh829EGtkMIXHdxhFpBKGHgCYUglTNzsC3AxEIAzXqnMDIkeNEGV2QB3IIwPG/rwhRk5ND4BcGndTpsHdRYHbJ4grMM+uANhAYsTwFHBbdn1GcwGzogkO+bHlwZQ026DViFBsYmjeXtkXTm+sqAIH9T3omPIg+a+Efz0n2ZhaCSApzZOg1mvQkWWFjNyDcjSKWF3B6BRSrFlZTEGnAGEIyzS1TIoZBKUmjTwBsPYvKwIBekayungEymHlYBrGZdJGIRZzlsZ9obwjSXF+CSatQhFQA2HUS2HTiXDj26thlouw2ufttMZW+xa+S/27II0XLQ6kKVX4umNQo/w9eNt2PzbGuSmKyGXSvH+hR7q3fGfcewShNzr+1aUYesXZ+KOhUV4aG0lHltfhW8sLUF5lo56NkR3JE0tx9nOYYEBGQvE6wbLab+S55qmlmPNtGysmZYdp5vzzNYfW7Zt27ZzzB3jGjcWz/zHf27Tzr0ZQ+4AOu1eauGJuMu5rmF81NiPvHQ1tt42E/OKONf9kXUVVIBnRt5o114gFMGB+l6caB0SWGn+LE7cNofXj3NdDswrNODZL81GXroamxbm46LFAW8whAKjGv9rdRm8gQi+sbQQbf1uuHwhsOBSfhqFFDkGFZb/n/bOPT6q8s7/75NkkkzuN0hCEkggCTdBbgqCKwqKApZ29xelW9daWtf+frarttKL+NsV6hbtFn/KttbKWq2rthZxW5VLRVBBuQTlbgjkTi4kgWSSzCQzmUyS8/vjzHNy5syZzEwuEGw+rxevMDPn8pznPM/3+T7fy+ebm8zf5Y7pt+iMSnNW0ex3UGgH9JkLrao9Y3F+Mg/fms8XdW309MrqZLspP4U73SsWwJpXPqOuVck2/M6NE722dM/tOcu9LxVyoKyJBROTmTZO2arljo3lw3OX+JcleeoqPSE5ioPlTdS1dpI3NpZn7p5FfmocUeGhvHPiAj29MidqWrE6u+l09ZJoDqerR9mSVbc43FsLWeHatHfj7O7l9AUrna5eJiSZ6XR1kxYXye3TUlm/ajrpCWZ+UXAty6anUm2xMzU9lrKLHczOSiAirC8OJj81ltXXjefBWyaREBXOkimpqlFPaJzadyH6XytI6lodHKtuZd2KqV5G7XtfKsTm7OFoVQs///sZHit2vNnEkcpmvvfGMf5yvI7tp+vVvtX3tX4M/MsfjvHXokZmZyWwalafbSsY7UKMD2d3L7/dX8GZC1YmJEex9fNanwvChg0b6gMRFiPawDkuwczYhEjumpfBofIW1YgorNIrZqSy83SjGivwyqeVFFZaaG538sqa61X7xLdvzPGI+QcMMwi1eOIrM3jiKzPUz2LAaA2Ugg9gwcRk9v9kiZp7kDsmhn2lTRTMTfXIGPQV36+1xPfHNg6eBqw/fXeR1+9bvjnPI0mtw9nN8s37mZ4ex8O35qsxEnfNy2DDu0U0WB3MyoxX+/DXe8vplpWALC3J79zsBM7UewYppcebWTYtDVuni7c+O+839P7Dcxe9KAaFZ6OquQOHqxcJxdbxozsm9xkqE6PY+lmNRyj0n757A0cqm7F0uMhOMXO8uoWb81Notbu4/6YcShs72HW6gVcOVnG4oplHl+WrRj09PZ6Rd0R4EDbuKGbx5LEe70zLIWGUFfrM7hI1mU14IMDTqKxNRCtt7MDqcKnnaNMAxDuPM5t8Zr7qjxVG/zP1VtXVq893GQhGtLDQRnA+onFXC6u0mKwicEZEWpZd6uBbLxdS3+bwyjEIlcBsCqWk0cr2U3Vq1qoRKYmRGrl5TwllF9tZt3IqSdFhhIfCbz8uIyk6TI0uFCvc7w9W8UZhlWF2qnhpRyqb2XfuIgVzMlSL+qsHK/jZu8Ve0Z+mEHjjnxf0y7EoIvgeXpoHoE6C4nqbqkordVNbPITnPf9VCPTVDA1BCfAR9qBZmfGcqG3z8AgZFSvSYnFeCocrmnG67RR6QREeKtHhNlw63JqQjMLq/fyH5YyNjSAiTGLb0Rps7mhRbe0TEVOg5UhdszCbZ3eXegSziTGy6a5rPQSgQF5qtBKAlqq4iZVYHDNOVzzpCZFs3lvKS59U8IuCGdw5M0MlmRE1OcRnMV7uvWE8TlcPuakxHgljWkEvFi/R9jULs9U+7nB2s/rFQ+piKIKowD/rutZD8+iyfA6UNaveEKvD5RHpe7K2had2niXEHBuQ/3RECwt/0HbgWl24t69Eph5Ziejc9UWjEuykiZHQQ7/iaJOBntldQm2LXU1G2/R+CTfmpaghtMq33tmpGe5EMgGxCkWYQtVBten9Ei9BAUoMiK+2GvUJKMlJRfVWpqfHedxXDKC3j9Vi7ez2Sn7rBbZ+VkNspKK25qbGqCtsIEiJMbFk6hgOVTRjCpHIHRtN+aV2j/7o6pGRZeM8/46ubvaVKu5YbZj92LgID+pBAGunk9cP1TBnQgK2TpeHoADFoyIyl41W1uc/LKeutZPnPyznzpkZyns+Wsf8nCSV61SETQvNUlvER89I9vDSPEMCaC2Ediw0i4K5mdS3pfHM7hKqmhV2rDMX2tQ4CaEZLcpN7teFqmd2E88r0tKFtqEcW0NdayehMckBvdQRbbPYsmXL+gceeMDn71r334TkKP7iZsoCyB0ThcUg4EWgYG4G994wgeb2Lp9FgPWWbC1r1k+WT2HOhAQ+LWlCkuAnyyfzzRtyqLjYzoVWOxKKFhOmY5z63X3zPLwiRoWIoyNC+ORck1eWXSgwLjGSX+w6y4HyJnYXNfDS/gomjolWz9X2SbzZxLLpafzTggksm56mfifsMjdPHssfj1SrOR56OF29bPjqNcRFmvj2ohwvw+CkMTE0tDmoa7F7ZeOumJHGa4eqcbh66ZXB0dXtVb09TPLNxuXo6iE6IhRnt0yiOYyeXuU67Q4XFZfaeb2wmv85Vssd16Tyh8IarM4eunsUj9AZAx7QL2rbWHOjcbW0NwrP09TexZjYCP5pwQQyE82cuWB1BzIlkx4Xga2zm8fvnKp6lqIjQjha1cJ3/i6bQ+UWMhPNJEab+KKujQNll/jFX8/R5nBy82TvEAaxBXl0WT63TFE8JtuO1qrU/NdlJ5IcHcF9Cyfg6Orl3hvG89O3T3OusZ3SRlu/JQb0Xg+taxy312Z+TjL3LJhAXmoMJ2taqfnglfOBMGVJSsDlyMS8efPkzz//fEDnarcQ9W0OntldQqRJYl9JMwVzx7HprtlD3No+iAAnEZ674d0iXjlYpQYGDfR683OSvCj8QdkXD7Tm5vLN+32S7E5Nj2XXw970/0bI/ukOj8+mECWd/GB5E2ZTGDkpUWrtl752JzJtXCy/P1BNTHioR8i7FrERobQ7e7yEp/hNaHTLr0klLc7MKweruDk/RSVsFvD13rWTV4Ts64OtfEH7rgEPu02IBBVPrfQ6Z/WLh9SISrGd3ry3VN2GLJ+eisXuUikFtPa2qemxbLl3nt8ALSPuDv24FJAk6agsy951JXQY0ZrFQ0/8cv1LTZN4fm8pD90anGFGu8IKb8ebR2qoa3UQG2kKquhMsNBrJJPGxBAXaerXIxLI9VZfl0V7p4tWu4vrchKZlh5LTHgYP1k+xW+wly/kp8aw45QxB8dv7pnj97piYkWaJM43KwTCEkoJwAPlFuZlJ1He1EFaXCQNViezMuO5MS+FM/U2lk1L44lV1yDLsF9TADl3TDSpcZFkJkZi63TxozsmkyZMKUgAACAASURBVJ2i5MhoWxkZFqJqdHGRJr67eBKzxycSF2ni+0vygF5O1PRRE56pt/HeiVruW+ipYYjxIZ7VV+CWHttP1fGrvWWkx0ey5sZs5uckY+t0cdHaicPVy7cWjTfULPTapHi/9ywYz7gEM1/UtXG0upVqix1rZzdvH6tjVmY8aXGRfHHBGpALVe/1At8xJoF6Q0a0ZhGRnien3/ccAFVPe0voYKFfQUYKBlLNfaTA12qltSsII5swrAGG8Qc2p4vYCJPHaq5fIVe/eIDCylbm5yR4eYN8aQR6reeRWyf5LPunbbu/97Ho6b2qPSWQJK/+YKQJC81Ca5QNlElL9NsDN+VQ1eTw6Fv9dQLVLK5KYbHhvdO8eqCayakxLJiU4ldVHOnwNeGGEsMlkNTs1U4nrx+qJndsDBFhoV7cIQOFmJAZCZGG3Cba59IapLV9ufat42w7esHjvFAJoiNCsTt7+P7S/oWHL2w/VceT7xUzZ0IC61ZMG1S/Go2BYN6ZrzGkpfoDPLZL4thAhcWI9oZEuKtgSyiDBlmmuaNLNZQVN7ZT3Nju1/c80hFsMeKBwB/z00AhPAz5j+90l41U8kX8eW0ChRG7lRba5zIqZA2w6a7ZZCaa1YJK4M5E7lRsJM/tKR+QsDBiUh8ojMaAcE1bHS6/ti4tl6j+ulaHS9UsjJLnAsWIFhb5qbFEuOn+9SxVAtHhoUGxGY1E+KvmPhQYboH04C0Tef7DclWzGEitTyP4mpB696lYfX1NqkduncKfj11QCzONNGgZ0ta+dZJHl+WrwVn6IC0jbP2shsJKC1s/qwmolMBAxtuIFhagrCw/ePOET+bsjq4edp1uGFE2iJGI4RZIj9w6ZUCr80AxEE3pR3dM9skAH0zh5OGEtrC14JoQsS79QQSraYPWhhojXljcOTODJ98rptGmpEBHhEpqRKDA20drsTldzMiI41d7y0iICmfjP8zwMA5pDWbr3ylSjxnMABHcE/pIvVEMDfpzYQ5EUypt7PD521BtmwYLwbGSmxqjMswH8ozTM+IobrAxPcM4w3gobFYj2sAp4iy0hqQFE5N44l3jgkKxEaEq78GszHgiTKGqf3pCkpm/nzPOY986mPgE6POXw+Ct4aPwhtY4F2z/+qLOW/rMxz6pAIc7/iYQDHRS+zuvPyP6l8Ib4iso69WDFfzHX8+RNzZGCdaRYfaEBFWziI4IwxQaojIq+UKUu5DNrMxYIkwmQ5eq9iWAp8tpVLMYXgQTHCWOF+9n444z7CpqZPn0VF64t28eiIUnPT6CksZ2lU9DYOt3+8+9uVrRnzD5UgsLf9CuSBkJEdS1erM4CSo9LWLdUYRjY0202l2EhkjIskRnd6+hyynYwTyKoYOREBfkxA8vzePNI9U02pykxkZQ+Lg3aa6A3q0qAetXTQ244vpwYbDbhmDO/1K4ToPB9lN1/OtfiujpkZkzIYGCORnERipBKCK7LjMxksLKVqamxfC9Jbmse/s0VmcP0eEhdHT1quHGF21uIdIjqnQoQiIsRMmkFDTtokoUGKe8D+aFawWRtqarv1XvSgd4DbUA1T+PCKzLSjSz7VgdhyuamZYep4bTixT/vNTofl2uApvums2MjHjWv1us1mR54t1ikmMihswtOhAE4zbVQt8/MHSu8i+NsHjyvTOqpvBxSZPH3kwQ6xbMzfKI+hODQXRwQ5uD8xaHzzKF3b1wsLyJfaVNxJlNfH5eqYJmCpF8lqwbaGyDNv37tUNKwZ2NO4r9ZjMOVTzFQIWOvtyhKAGpjeLMS41WhZ+WmPaVTys5Vt3Kv35lqvpu9M8juBmcrng1V2ZaepwqJDyNoFlqXc/+cN/CiUxNj+fuF/vqijy18yzXZiaOmMjaQN+Htn/8caMEiy+NsJgzIZFdXzQSHiqxdKpnKrU/y7ngx9D67l85UMne4kavFPNeGfUlCLp1QYKix1DFNohkTVHUqD8M1T0HKnS0dVxFUebnPyqluxd2FzVQ3GBTSzqufvEw12cnUljVwp+P1VLdosRAPPneGVVY6J/n3hvGU9ti5/6bctTJLIQQ9An+pCgTu4oaqWuxk5EY5XeSXZ+TzIZVU9m48yyRYaE8tmLKsAWyBQK9JyRQTUNbrmCobS9+bRaSJGUB/w2kodAcbJFlebMkSeuBfwYuuQ9dJ8vyTvc5jwHfQRkTD8my/L77+zuAzSjZ1i/Jsvx0f/cOxmYxnOq3Uo/0BMfOt/K9JRMvSzyBUOcb2hxUNdu9Km9fjvsPZgv1rZePePFKGNmJBLRV6GPDQ/ndmuu8Sk3q7RJiAmuzegVfQ2psBI02p5rJORBv1eXe0vV3v8FmLveHITNwuuuYpsuyfEySpFjgKPA14G6gXZblTbrjpwF/BK4HxgF7APGWSoDbUOqefgb8oyzLZ3zdezAp6qO4cli79YS6X9Yj1AeHhTlM8qgQPysz3mPLJYzWaxZmqyuumFD/57XPVc/HymvT1QQqS0d30HU3YPCCcqDn9pcfMpDnCBRDZuB01ymtd//fJklSMdCf5eerwJuyLDuBSkmSylAEB0CZLMsV7ga+6T7Wp7AYiQh0MAi36oU2u2ownZoWw5ZvXnfF97/DjTJ3PREtIsNCyEmJ9ihGFCopxMY2Zw8yEpPGRKkFjdITInn2gxK1n7XbEX3/1btLO9a3dVLa2EFdayc/33GWx1dO4fqc5KDV8cFsPwZzrtEWMtjrDac2FJTNQpKkbGA2UIhSA/X7kiR9E/gceFSW5RYUQaKtQFtLn3Cp0X0/P9gG+4p70Kb1DvVeTcu9mTs2JiArs5a0VaC4oZ2/f/4A6786jdLGDk7UWNhX0szi/GRmZSUZvuDhVoUHmrbvq101Fju5qTFYOrqwdHTR1dNLV49MZ3cvCyYmExEWovbL38/J4ObJY/jJtlN0dPUyOyuBm/L6SjxqJ0l/4errVk71SOl+4eMyunpkD8q7YKCftMH00WBsRkbPGOz1htPOErCwkCQpBngbeESWZaskSS8AT6J4m54EngG+jXFdVBnjgtleCqkkSQ8ADwCMH+9NH6btDOgbUIK5eTjCdrXcm0BAVuZHl+V7aRagVBzXliQEpS7rvhKFI1EvOIbbyOav2K8v+GqX4K/0FYsCOV6u1aNVSrWs2EiTuh+vsdgDDnXWlxV8fOUU1fA8EOgnreijJ94pIs5sHLzn69zBItjrDWfCYEBBWZIkmYDtwPuyLP8/g9+zge2yLF/jNm4iy/JT7t/eB9a7D10vy/Lt7u89jjOCkc3iSmsWA+FpEOcfr24l3BSi1lEVmoURtJPtatMsgm1vMOdcCQIjcc9j5y2qEXZcXDgH1912We4/3BhKA6cEvApYZFl+RPN9utuegSRJPwDmy7L8dUmSpgN/oM/AuRfIQ9E4SlAqr9ehGDi/Ictyka97/y0YOGssdlb/9gAXrF0e33/y41uuCtvG5fYY6PkrtXj1YIWqUQxHBGawjFuDweXs16GM4FwE3AucliRJ5PeuA/5RkqRZKFuJKuC7ALIsF0mStBXFcNkNfE92871LkvR94H0U1+nL/QmKvxVkJUVxcN1t1Fjs/N1/fKR+f7K2hVW/+pRWh0vdq4UADw3jAO0P/gokweWJRdDGEeihpedPjonwIqwdLPTBer/5qGLY3sXl6FfxTgk1hQdyvJEdwQOyLH8qy7Iky/JMWZZnuf/tlGX5XlmWZ7i/XyW0DPc5P5dleZIsy5NlWd6l+X6nLMv57t9+PrBH9I0jlc2sfvEQRyqNVXs9aix2nv2ghBqLcY2RoUIg98lKimJSihmASSlmJSJVIyhACXL5T03W7OWEGLzbNDVgQdkb+7LhbD9Vx6Kn97L9VJ3PPvD3zvS/C/uE0RZk7e35xEaEsvb2fJ587wx1rZ08+d7QONtqLHZum5bqYZBLiTaRt24nz+05O2T3EH1UMDeTNQuz1dQC7TEb3i1iw7tF6vf99WF/Y0+UUwyNTvRmFTaAX2FxJdHV3RvUZNYa6wKBrwkw1Aj0PnvXLqHq6ZXsXbuEOROMw5R7geXP7Wfuz94n+6c7WPT0nmEXduBbKAgDnDYNXLwzYch9audZn32wcUcxhZUW/uUPx9Xn0F7D3zvVHnvfwomc3nAH9y2cqPbflPQ41r51guXP7Wft1hMD7qttR2vZ9UUjDy3No+rplTy8NI8LVqUC+q+GSIBv3qPEWWzeU6KWLHzlYJXaZ0cqmyl44QCvHKzy+P6Jd4pUA6xRu4dqjI/ocO8We1dQqlh/KqoRLgf35UDvs27FNPJT41iUm8xj/3OK8kt9g1wbq1DX6mTpMx+x99HhtXH4s8prizsXN9iwOlwe/JkiR0PfB4LspdHmVGuQalVwf+9UhEGXNFjVWhvX5ySr/SciPkHpt4zEgXkr9O+wYG6m6qKNiggN+npGEPEp4q/+ns/sLqHRpti2ZmXG9/WlUHcM/JD9jT0RUv7DTS2NXj8aYESnqM+cNUde88s3R0Qiz0iA8KoUXbBy0eqg2d5XSezhpXkkRYexcedZ5F4ZV68S9JSXGkNlUwed3TIp0SamZ8TTZncNGfu2gJ5FOtCwZKMs1WCMeyIMWoR36w2fNRY7m/eWUFRnZfo4pTj0UI0lo0I+g4E/T48v/hR/5/nrz79pPouRhuGybC/d9CHlTQ4SIkN576GbWLF5v8oU5g9mk0RiVDiNVqcafj1pTBRP/cPMAYUV91f7Yyig9XQsmZLm5T7XZrIa5ZRcjgXnStMD+IK/MgPjk6P/tvgstPD30oaLuszXsYO1bPtaOfauXeJx3Nrb8/nZe8Xq5JeAKWnemsXnVS10dPXgaPMkBSq/ZGf1i4eRgRPVLbz6ncADbPtj1h4otH3479uLcfXCv28vxtLR7dGfok/F6q7NpxAEuMHyQgwEVzJLtT8IobooN1kNodcHNwaCL6Ww8PfSBvpSgzlPX89C+zdYBBpled/CiSyZkuaXfOZIZTNPvFNEm6PLQ7OAvpDafaVN1FjsA14hhcqcHh9JdESYSkSkNYTqNQJQDJ65Y2N4+FZP24UIhnL1KhXRY8JD+fWHpRwub2LT3bMM3bkiJR6goW1gJQCCWSDE+7V2Osl/fCcP3hJchvJwaCbafv7Bn06okcNaARIamzIukGt9KYWFv8k50MkbzHn6xKfBrDTBGG4DWeGvz0lm1yN9BY+f23PWg8hY4IFXP2PXDxYH1VYxOPedu8iJ2jaP/Bgtm5iY0BnuujDC2yHOERwUoPShdhV8dncp7W6ikcKqFtUwCood4c0j1SyfnkqVxvNR0uid3BYIglkgxHvOW7cDVy88/2FwxYuGQzPR93NGQqTHmHz2gxJCoxPSA7nWl1JY+JucA528wZw3lDkC+tyHoYYsG3vQixvbyf7pDhLMobQ6egiR4CE/pf7E4CyYq2wJtJqFESGRXrMQ7r9Fuck++7BNU3QnI14Z/EJIvXnkPI22Lk7VtfHs6lnc9/IRHK5ewk0DixIYyMIyZ0IihZUtPt3fQ3kvX9AXYcpLjeal/ZUeBEo1FjtWh4seu/ViINccNXCOghqLnZX/uV8t5xcI5uck4nT1enlVBqtKC2NcamwEf5eXonovxGptCoFnvz5LKQ0xPoF1K5Uao+K85dekcqq2TfVQXIlcEtEHSdFh/OajCo92Xi4YGTX134nPtb/5Vn239ZLfrciosBgF4F3F+5E/HueC1V3YKUzC2W08TsJDIDI8DEdXN+MSzPzyrmsHNSm3n6rzqEAnBvbat06w7WgdBXMz2HTXLK82g+8K41cK2irrkWEhXJsZj7O7l/tvyqG0sWNY2+qrbopnBrDSZz9cfs1pubtrpr9rfim3IaMIHlqVX8lX6aPPFwNvy75SHN2e53X1Qlen8uV5i4Mn3inysIcEi9cOVauCYmparCoIHl6aT0ZClIeKrt/jjyQPBIjSmydx9Sp8HoVVCsGztoTi5WyziArdvLdUtR/94LZ8ftjj6vJ/9qiwGEUAEIJEmzK/71wDJ2ptXsdetBkXsA4UggdEG3ikDX4abAnDy4k7Z2YwNjaSe/7rsFetXgnFjjBc8GUs1fZZsIlko9uQUQwK+rRtgKqnVw7pPeb//AMabV2kxoZT+PjVxyFx0398SLVB9fbBls/sD4HYjoK1WYxqFlcRRmqEoB7ZP93BpBSzV9DYQCHKPAjvgjbsPc4cRlWT3aPWiBYjoc9MocaemMJKC5P/7w7S4sxkp0T7DMPX83T4eyatUVcf16Kt3bL95AUSo0zUhYVHBvIcIzrrdBSe2LjjDJv3lrJxx+XjOPaXQh5iRKIIlDcNLAjKCOtWTGN+diJ7zlzkuT1nefnTSrYdq6O4wUZhZQuNNifr3j5t2GaRhi0Me1cCs8cn+PzN2a3YevaVNHGito2NO4q90sq1PB3gP5PUKFNXnPODP51g895Sntp5lvKmDlrsLkIiogLy8Y4Ki6sIx6pbPf5eDogU8o07jCvXP7R00rAPoqykKI5Vt+DqlXn+w3IqLnkHWGkrowdLVTDceHhpPsuvSSXRbGJCkrnfY3NTY7yEgZanA/rnEAHF7jM/J4kVM1JVoVkwN5P5OUnUtXYyPyeJx1ZMISVaMVXI3V0BSfbRbchVhCnpMTTanExJVwJrLkcMgc0dAFXc0Mb2U3UeLr8aix1ZDmHfj2/xYPkaLIxcouOToihvspM7NoazDd6G1ZyUvkkoqpbde8N4rs1M9CD+vRLbkqykKF74p748LSM7T5RJYvV1E1SXJvQZI+9bONGDJtBfwF96vJkFE5N580gtxQ02tezlpruu9Xh2UdHth5vqywJ5jlHN4irCJ6XNHn+HagXtj02pskn5ztmt1P/UrnjDRR6kva74v6VDEVpd3TIrZnoTOxU39FU/E7VDShs7vMh5LhfhUX8w2rnZXTJPrJpuWAoiWDY31RPivpGI2tT3hSp0Rl2nXz7ct3A8rx6s5r6FSomEYMl+fKG/nIT7Fo3n9weqGRsbzmMrpqiaBXi64YwyGMUAD6aqeo3FTl2rnVmZ8SzKTSY9XtEY/nK8lhaHCyR470S913lxkX0ENEbuQaNiRVcK12bGqm7nEBT2s8X5yax+8RCPLsvnQFmzYbmLQGMytEliu043APhMCqyx2ANOJLuqXaf6gWD0WdQLrW/rHDLCF20m3/MfloMEG1ZN97p2oCpvMMdpn6e/SMBArvnqwQqe3nWWpOgIZmbEkxZvZm52giEvhD8YqdaJZhNp8ZFq5mdqbAS/+sZstXK6toTDvTeMp7Sxg5IGK7uKFOImQaBTY7FT8MIBlSVKWxdVYPk1qaTFmam41M7ZBhsP3jKRqiYHx6tb1FqnIjU7WIoBX8cabQMDOS//8Z109ciESLDvRwrDmWAtn5UZz+zxir1RG2Uprmd0z+2n6tTw9zU35nhUrBdMYfNzkth017Ve82Tb0Vp++s2VOOtLfZiq+3BVaBb6FyA6LCvRzLZjdVgdiooqBgag0rMJSjVgyAoQ6TP5fF3bXxahCDaamRGvTpD+Vg/98zz53hkabV2GXA2BZDBuer8Eh0tWn2FXUSO7z/RlgQbTVxGhoOfdaXG4aHG4VBarRpuTx94+TX5qDLuKGilpsFLVbKe4wUZJg5UWR7dqdNM/ixAU4CkoUqLD+cq1ysKo7RvhQQAlnkHP4ZCXGt0vy1Ug/Se4MoRNINDzHrxlIr/5qIIHb5moCpRHl+Xz2NunOddoU4Wb+E17nSfeKaK4wcYT7xSx5ZvzePnTSt787DwOl8yuokYsdheFlRa16NaahdnMz0misNJiSFtYMDeTH3W0eqtqBrgqhIX+BfQRmsTy8NI8D57F+TlJKsnHotxkrA6XuhJr1fXBGLq0GZNCszDaChTMzcTqcFHSYOVrv/7US7MRhLbdPbKXdVu0T5J6+d0nVay9PV+9nnie9IRIdn3R6HVOwdxMQ1VcX1x37e35PL3rLGNjI/nRHZMpbezwygINtJ+eWa0kd8VEhFLutnNEhoVw57Xp3D0vS83+bOt0qd6c/aVNJMcowkE4M67LSaS+tZPcsTHqylowN5PD5U0ccYdLT0mLoatbJtwUwoZV00mPN/Pyp5UUzMngks1J0QUrk8ZGY3V0e1DpLcpN5nCFsuoKboendp41FBaBbFcEf6g2kzOQ8x65dQr/a854th2tVbdqB8qaaXe6cLh6iY0I9X2+hm9Tv3hMSDKTlWRmWno2y2ekqe8a+rQTfRuzkqLosTVd8NlYDa4KYaF/SPGSpmfEqSHIAt++McdDuPjidgiGO0CbgKNX0/vjXtQyNIO39qEltNVfR7QvPFTyqNupfZ4ai5381DiPAWGUKyEi9cRqI37XW9kFtG0RhLh67UXPPym0nEa3o0IC/vs716vC8Zd3zVSPb2538vMdZ+no6mF8RJh6/NS0WNYsyvHa/mQlRbFgUgqFVS2GVdSf/aCEVw5W8fDSPDbdPcuDD3TZ9DRO1rbw9S2HmJkRT2GlhQNlzR59r4VWOBqNC+3vRvkqgVIT6Jmq9BmzvgTzhlXTPexUu4sacPX0MisrgdhIk9oP+oLQ2jYNlD7hqhAWeoiXpKUJE0ICApPuBXMzqWuxs+/cRRbl9l9pWyvBg1XTk6LDiAwLITw0hHtv8Kzdem1mIgVzs1Tma/150eEhpMdHcr7ZjtPVw7deLuTJr83wsGYvyk1m7VsneXRZvs/n1hq8tKuNEVluoNDS/N85M4Mp6XE02prUDNXF+SkefXrnzAxVCD37QQldPTLzc5J4dFm+usq3Ndg4UNZs+C5UFiqHi817S7E6XKrQ0D/3otxkdhc1qP8X1+/ukSmYk8HuogbqWuN4dvUsDpQ1q8a/GoudtW+d9BCoegxV8prRu9Jrb0ZanZbb5NkPSlR70MqZ4yiYm6n2yXC4iP0KC0mSIoH9QIT7+G2yLD8hSVIO8CaQBBwD7pVluUuSpAjgv4G5QDOwWpblKve1HgO+A/QAD8my/H4gjdS/IC3Ljy+rcSCsRjUtDk7Utvmd+EL9B0/NIhBs2V9JZ3cvnd29lDZ2ePym3ztqX+6W/ZV0dPVSdsmuahcflzR5sEKBN+We0XNrVxLtRNQKQVuni4zEKK+tiqCL1wsg/cocZVK8EWlxkVyXnURspIkjlc2qsW3X6QY1ZgMJCuZkEBtp4qKtkziziVBJIjwshLoWu6HlXjzD9lN17D4TSUObQ227/n0fKGtWJ9GBsmYeuCmHTe+X8OAtE9l5upHiBhvFDTZqLA4PwbDtaK1aGtFXXMZQeVP0q7vRexNaXV2L3YM5TLvVFONSz8imnRtDldkaiGbhBJbIstzuLpD8qSRJu4AfAs/KsvymJEm/RRECL7j/tsiynCtJ0teBXwCrJUmaBnwdmI5SA3WPJEn5orRhf/C3YhoxMAWCQF2Peqq6YGjfH1sxhfXvnCEhyqSyFglo99B6gfjYiik8+d4Z5kxIZHJaDP+1v5LrspO8nm8w7tOCuZkq/d3O0/XYXb28eqCS1s5uShqsvHDvPJ8qq1ZTAPi0tAlQQpfPW5SK88Lg/MbhKpo6XB7np8aG02jr4s3PJBwumdTYcKpbHJRd6qDogpUt35xnuMruO3dR1RLWLMz20C7F8dpJtCg3mWd2K8ZOS0e3R1br3fOyPDQt/V4efC9Ug0GguR5CuBZdsLLtWJ36vbY9Wo1aC19zZjAah19hISu+VRFfa3L/k4ElwDfc37+KUin9BeCr9FVN3wb82l1c+avAm7IsO4FKSZLKUIonH/LXBl8vSHyvVaeDwXDT1YEyqUobO9i8t9RLxT5Q1qzuofUvVz8ZfVHZDeYZspKi+NU35nio3q1uborP3MbEQGHr8pb5YrBrBUWi2URoiKR6NxwuxbKZHh9JV7dMi8NFcYONlz+tVAX0kcpm7v/951id3UhAXEQojTYncWaTGpMgtiVCM5qbncDzH5bz3skLNHV0qdpCVlKU6rkAT01LO56EABqOuAxtrsd9Cyf6tJ/FRpoAmJ4Rx7LpaYaLoq9zfc2ZwfB8BmSzkCQpFDgK5ALPA+VAqyzLggqlFhAjOwOoAZBluVuSpDYg2f39Yc1ltecMClp1WksKO1IQiGY0lJydwSArKYpNd13L2rdOcOx8CxOSoyi7ZOe6nOD4I5dfM5ZdX3hSOQqviIApRFICqwCzKQSHxgdadrGDrOQo9XctntldgtWpDDUZxXOi9x4JW4Yw4mrd2hkJkWqMQSAI1C4x0FV67e35qmYBvseHfgtoZF9Kig7DbArhL8dr/dre+rtXIAgqKEuSpATgz8C/Aa/Ispzr/j4L2CnL8gxJkoqA22VZrnX/Vo6iQfwMOCTL8uvu73/nPudt3T0eAB4AGD9+/Nzz58/7bddgDHWj8MRQGMa2n6rjsbdPGxY8CpPga7Mz6Ojq9nD7SkCkKYREs4mFuSkelcOOVDbzTy8VqsliG1ZN9fLiGKVg9xcw1x8C7YO1W0+w7VgdBXMy2HT3rICv7wsifkgEqOk9PkKAhYdKPL5yCvctnMi8J3ermltylImj/7Ys6OcLtCJZUN4QWZZbJUn6GFgAJEiSFObWLjIB4autBbKAWkmSwoB4wKL5XkB7jvYeW4AtoERwBtKu4Shw87eKYDUcowH+8dlL2Jw9JESG0dbZjdkUgt2tRXTLqMa66Igw3j1xga4eGRncmkaX18S7PieZ1++f32/SnJERd6AlBQPtA31t0kDQnyASAVenalvU7Zloh7DDvHaoiq4emQ3vFpMcE+FR86XZ7q2VGcGXO9wf/CaSSZI0xq1RIEmSGbgVKAY+Agrch90HvOP+/7vuz7h//9Bt93gX+LokSRFuT0oecCTglo5iREJ4Y7RJZmLytHZ2IwNd3b3Mz04kLATmZyeqE2XTXbN4/f75JET2rVljY415WIRt5nIxdAeC+2/KISMhkvtvyvF/sBv9JrK5A64cLtnDBSl/UwAAEI5JREFUIwN9C+KNeSmAkk/y/T+coMXuUuO0pqbHMJwIRLNIB1512y1CgK2yLG+XJOkM8KYkSf8OHAd+5z7+d8BrbgOmBcUDgizLRZIkbQXOAN3A9wLxhIxiZEN4Y7SahfBAWDu7KK5vpxdYMTONP/3vhYBS1EiEOz9y6xTVqArwy7uuvUJPEjy02a2Boj+bwYZV09UqbmnxxrwXD96Sy8clTR7fyaDmfgQCX+5wf7iqE8lGMbJRY7Gz9JmP6eqRkVAG9eL8ZA6WNas1QEo3ruT/vP4Zu764yPycBCDkstb4GAyGOvBJXE+kL2hrfoAiZP9zbzm9BlN20pgofr9m/oDaEajNYpTPYhRDDuF6BHh85RRiI0L7aqiWNKtcmnMmJFJjsZMWF8Xi/BSOnW8dcoargfBBBIr6NgeHK5qpH2AdVT0271GMmA1tDtYszMbqcFFjsavP8J97+gTF8ms8OT3KLwUePjDQPrkqwr2Ft8PmdBEbYVITjEY9ICMTYl8uYh9+963ruPvFPq/5prtmeaRIa5Ohwnwk5Q22LdBnLBR5LQ/clIOlo3vAmoE2A3QwtVIEjrsT7Eoa28lPi1Pre4ASiGUKgd5eZYXf9UUj8W4DssAfj1QzNzvBr2F3oLEWI15Y6OP1AbUDtaHKQ+G6ulIYKno8tXJ5QiTR4X2Vy+vbHGzcUUx8lImqpg46nN2s/+r0AXsL/EGfx3G4oo/sd1ZWnIe3QVj5/+dYLW2d3WSnRA/pFsQo+1bUQ9WmsQ8oxkXS/R0Eaix2tYZrW6dLjfbV2hVO1rTwcUkTvaDGkRTMzWDX6QY6unro7O7lh386ybWZif0KP2unkxBJ+RsMRrywEPH6szLjyU2NITaizzDz5pFqGm3OoFxXIxH6/I6BQpQc1FYuBzhTb/X6zldq9lBAGwkpoioFY5OWY1Ic+8Sq6SyfkTYo1i9f9gOjfAmR3anVLAYCfQboYPDyp5U0tXcRHR5CU3sXB8qavXJHtPFE2hT0u+dlsfrFw4rnqUf2yh/SZwi/cbiGXhneOFzDE1+ZEXAbR7yB88+79/s0ImlX5PR4M996uZDyJjuJUWF8bfY4tn6mEL1esnUxPSOOiSkxQW9ZtKt1Wpx5QOc/9vZp2jpdrF81zXCCinvoq3A9+V4xU9Jj6ezq4YsLVn58R75hSrlA7roddOsYpArmZHD3dVn9ahY1Fjsbd57h2PlWHrxl4oBU8+2n6li79SSd3TJT02LY8s3rhtToJ4Kt9IluAtqivwVzMw23qP4Mkr7qgw519qYRIXFdi10J8Jqboaa9B3O/Vw9W8B9/Pcd12UkemcnQV3M1IyGSAz9dqnqj7lmQRVxkRMC1Tke8sAjUG6LPQA2RMLQa6y3M/iDozobifPGyjKCvcK0tqisQGxHK6Q13+LzX2reOs+2oEucmvA+Cmq4/bHi3SN3SmU1KYlcg5wlovR4C87MTmTYuHhicTUn0y9S0WIobbExKiaa8qcOrfXqqODEWgqHTC6TyeKDoT8iI/hbcHJv3lrJmYbb6+1Db4PSahcDfbEWygrmZ/PlYDdUtnYSHSqyalc6u0w1emkWwKqfIUhSaxUDOF5qFnmhF337t38dWTOEn207ToUnQErkEvlBj6cuFEPUh9Gq/PyRFhVPX1reXDSSUftvRWg9BAVBY1aIWAhbp7+LZNu44w7HqVo8qYr4mlzjnz8eVIKZGq/KMhyv6eCgANX1caCHa1O3Ne0rYdqyOuhZ7v7YtbRaw/v7BvvdAjIjHq1u4/6Yc5uckqdsKYdQMVMD5Q43FTmljB28+cINhdTKAH25qa+rvGgKh69evH3BDhhtbtmxZ/8ADDwR0bLzZRKujm8JKCz0y3DYtjY3/MBNnt8zPvnoN9y3M4ebJY4l3G0cDhcLxkEJpYzu2ThdnLliZNCbG8DpiYmUmmok3m6ix2Hm/qJE1N2aTEBXOkimpHucdqWzmh1tPMiE5imnj4rlhUrJ63sufVFJv7SQ8RMLZIxMdHsKqWRnkp8b5vJ85PISTNa3cNS8DWZZ4dFm+urr3h0ljYnj/dD02Zw+drh56ZSiqa+VAWRPvnbzAX4saOVHTyvHzLSzKTVGfQdx/UW4ybxlEJE5Ni6WpvYvw0BDeOXmBMxes1LY4eOtYHR1dPZysaeU7NyrbKlE5LC7SxA2T+iZrvFn5/Prh87TYXUSaQnF299LU3uXVHnGNaoudnyyfwqpZGcSbTfzmozIarE7CQ0P4+vXjvdopsPXzWt4+Vse4BLPaBnH/YMZNjcXO4YpmZmclcM+CCV7nThoTw/HzSvr+hRYHR6tbiYkIo6ldaeM9C8bzflGjYX8EgxqLndW/Pciesxc5U9fGXfOUbIsn31M0m/ZOF8kxEbz/xm8urv+3f33B3/WuCmFRY7Hz3Acl7Dt3iajwUF76pJJ95y55TdrMRDPIqC9JSPfBdDjAc27atjP1Nk7UtPq8nhisyIob7HBFM7/dX0G1xc7bx+q8zrv/1c85UdvGsfMtWDq61En/8qeVvF5YTaerF6d7xXb1yBypaOaBxZO87hcXaSIz0czPtxdTerGDS9YuTl+w8kVdG4vzx/gd6PFmE//vgxJcPbK6deuVoa61kwark9TYCDq6emiwOj2eQdx/XIKZs/VWnDqDyVdmjmNx/hjuWTCeaoudwkoLs7MSSI+LwNbZzeN3TlWFX2aimbhIZUU1au/pulbO1NtYMmUMMeFhNFidXu3JTDRz/HwLR6tbQYabJ48FYOKYaEobbaQnRHK0qsWnsPfXBl/QC+2XP63kt/srWJw/htuvSTPs73MNyli6MS+F26al0dTuZNuxOhqsTsYlKBrsQNqibc/himY+KVM8UfWtnfzDnEwP4Wnr7Gbv2Yu0n/hr77899uNn/F33qtiGaH3xZ+qtqg3AKB1dWN837ynheE0r4xPN/PHIeV7cV861mfEkRYdT0tiukr0G46ablBKtsDm1KoEygsZey4c4PtHM20drsTq7WT49lTULs2mwOtQ6GB5wu9zaOl0eKmvB3Ez+cryW8+7K2+Eh0NULCVGezNdaFXnb0VrV42GxK1wRxQ02Vmzez1P/a4Zfz8eP71DIe0MlifauXsKAuTmJTEuPZ/mMNNWboSVDLmmwkhobQVJ0GFPSYzlS2aIGX83KiPOgshf1MHyp1WIb8fKnlVRcaqfogpXrchJZs0ihtr97XpYH36XYGunJb2aPT+REbRvHq1vUbcr1OcksnjxW7WMxbvRbn4HSBOi3HHmp0WQkRJKXGu3zHG3IdVZSFBveLVL6LTN+wJQFQkgI0qE1C7OZn53IZ1Ut9Lrb+YPb8lm3cqpHiH6g25CrQlhomY+0A1e/jxQvTUtMq0WhjtBl445iFk8e63dfKF6sCMMtbrCRkRCl3kdEHBbryupVNdvJT4tTBZ2e/Ea43rR5FaBMnCVTUtXzlk5LxdLh8nLR6eMVBOvVgolJ1Ld2cq7Rhs3ZE5CbdGp6PDMzE7G6yWfm5iR5uHFFu4VxbndRg/q8v/mo3IOqH+BEnZWspKig6N30AVq7vmjE0qFQ25+obuFYdQtJ0WEsmZKmTjR9BuW3b8xRFxStwXNRbjJ1LXbKLrarQnswRDBa6O0arx2qpq61k9cOVfv0fumrnOuFx0Cg7T9hr9LXCQFvwqQffpkqkulT0H1pA6IzFuUms/WzGj4+d4mmji5MIRKuXpmM+Ejsrh5a7C7MphByU2MCGizauAHtvUTClJjEG3cU02LvUjWC6Rlx5KVGExMeisPVgyR5qun9sVwtn5HG8eoWcsfGMCMzji37K7lo6/Q6TjsQfvWNOR4rpdYK3h9qLHYefP0YTR1djE8yMz8niXtvGO+xYov7aHk0BQRZrxGCMRBqBZ7ZFMLNk8eomsVLn1TQ0aWwS1k6ug0roIHyrkRym56uMCMxim3H6lShPRx8moo2Y8bpilf7UO/uNao34kuTCMZ1q11UtcbooSJWuiqERTCwOlzsOt3Aw7fm8/CtCglrUnQYW/ZX8tiKKYyNjfSIzdBTufcHvdDKSorymOx/+f6NfXEVDhc3Tx7DS/sraXd7NF7cV+GTHk+PA2XNnKhtY/HksWzZX+mzxkV/rE56aj5f2Ha0lqYOZXGJjgjjT9+9QSV2ER4EcZ+pabHMyozn/ptyOFrVV819X4mxsDCiqusvxmHdyqmG9S6SosNUdqmp6fGq1yI93uyVQdkfXaH2/8PBTrbtaC3bjtbx8NI8lU5RX4LBqN6Itg98cYAGEj+idycPJN7EF740wkIfFi72pT+4LZ/VLx5S1cI/ffcGjwnuayCLSX+p3cnY2AgeuS3PZ6lALQ6UNVPepKQsv3aoWh0Y4JurwQjaQa6tntXfcaIftAMhkECjgrmZqoq+buVUwJvYpWBupsegL23sUAdmjcXO8WoLJ2qtarsyEjyf1Z/Kb/S7fgsjAtKe/aBEFQZGtHe+6AqHI5lMixqLHavDxZqF2R7CSV+CwajeCBj3gd4upd2maV2sm/eW8sLHZSqDlrieEd3kQMlvvjTCQhsWPnt8oseL8MeAbfSSntldok56a2e3WidDe4wRtBNPaC/IeEzEQKAd5FlJUT41BP3qqH8WrR1H8FAaMVbr4w+EEUz0WX2bA6erh8V5KVQ2dfBG4XnyUqPVdl2yeeYZzMyM8/jsT+XvTwPQ5nUUzO0j0TVi9jbqE199M9QQk9Oo9KC2vuvJ2ha2Ha0hLzXaMKZET3qjz6MR0GocL3xc5lGMSn/8YLdacBUKC6NQWaGmmk0Sl2yd2Jwu6tscHnH009Lj2PpZDbtON6gBMKJU3703jPcggH31YAWnaltJiAyjo6ub8NAQnK4eIsIk9p1rZN+5i6xbOdVjAIhVfPOeEoouWMlOiWLrZzVqOb3rchLViuCB4Lk9Z/n13nISokw4XD10dPUSKsG/LJ2ELId45VuIQacPLNJqBMIa7k8bUULNz5CdEs3arSdpdbiwuVmvQlGKvgD86K2TvLS/knizySOQC1DJe31V+NIGBomaIrGRJvW9ic/xUaGs2LyfaelxFFa1eNTQ8BVsJd6DENDCxjVUNgowTv7TCzHtO4K+xMf/PlhJi6Obf/3zF36D0rTQboNF3o04/vGVU9RtmrZvjYymX0ryG0mSbMA57XehsSnjQqMT0nvcxVzF/0MiY5Kl0DDVt9jr6mwPMUXG6P+v/Sz3dHdJoWHhva7O9u7mWvU+prETZ0khIaHIsoxSxsALva7Odrmr0ybu32NruiDa5ut5xHGBPHt46qQ5hvd2t0n7TNrravtHvVeoKTw0Kj6lx97WhIHlW3+OaUz2DG1fGkLbN3JvL1KIJzeKDF2NZUcN26O5p693o7+PLMuyJEmS7HJ2SKaI6J6O1vqQcHOcZIqIll3ODldzzVn9tcX1tO92kEgBmgDCkjMnh5giY4yub/RsPXbrReTenh57W5MpOWuqFBISJvf2drsuVpzUnuNnjKj37w/6+wcw7ibIsjzG33VHumZxLhAGn+GEJEmf/6234Urff7QNI+P+o0xZoxjFKALCqLAYxShGERBGurDYcqUbwGgbRsL9YbQNV/z+I9rAOYpRjGLkYKRrFqMYxShGCEaFxShGMYqAMCosRjGKUQSEUWExilGMIiCMCotRjGIUAeH/A5oHf7J39A/DAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.spy(adjacency, markersize=1)\n",
    "plt.title('adjacency matrix')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 1\n",
    "\n",
    "What is the maximum number of links $L_{max}$ in a network with $N$ nodes (where $N$ is the number of nodes in your network)? How many links $L$ are there in your collected network? Comment on the sparsity of your network."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum number of links in our network : 11085570\n",
      "Number of links in our network : 37274\n",
      "The sparsity of our network is : 0.3362%\n"
     ]
    }
   ],
   "source": [
    "# Our graph is directed therefore :\n",
    "Lmax = n_nodes * (n_nodes - 1)\n",
    "print(\"Maximum number of links in our network : {}\".format(Lmax))\n",
    "links = np.count_nonzero(adjacency)\n",
    "print(\"Number of links in our network : {}\".format(links))\n",
    "sparsity = links * 100 / Lmax\n",
    "print(\"The sparsity of our network is : {:.4f}%\".format(sparsity))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The maximum number of links $L_{max}$ in an network with $N$ nodes is $L_{maxUndirected}=\\binom{N}{2}=\\frac{N(N-1)}{2}$ and $L_{maxDirected}=N(N-1)$.\n",
    "\n",
    "We can see that our network is very sparse which makes sense as we work on a flights routes dataset. Thus, small airports have few connections compared to the big ones."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 2\n",
    "\n",
    "Is your graph directed or undirected? If it is directed, convert it to an undirected graph by symmetrizing the adjacency matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our graph is directed. (Some airports don't have the same number of incoming and outgoing flights.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We symmetrize our network by summing our weights (the number of incoming and outgoing flights)\n",
    "# since it allows us to keep the maximum number of information.\n",
    "adjacency_sym = adjacency + adjacency.T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 3\n",
    "\n",
    "In the cell below save the features dataframe and the **symmetrized** adjacency matrix. You can use the Pandas ``to_csv`` to save the ``features`` and Numpy's ``save`` to save the ``adjacency``. We will reuse those in the following milestones."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "features.to_csv('features.csv')\n",
    "np.save('adjacency_sym', adjacency_sym)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 4\n",
    "\n",
    "Are the edges of your graph weighted?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Yes, the weights of the symmetrized adjacency matrix are the total number of outgoing and incoming flights from each node."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 5\n",
    "\n",
    "What is the degree distibution of your network? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The degree of a node is the sum of its corresponding row or column in the adjacency matrix.\n",
    "degree =  [(line > 0).sum() for line in adjacency_sym]\n",
    "\n",
    "assert len(degree) == n_nodes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Execute the cell below to see the histogram of the degree distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAADTZJREFUeJzt3X+o3fddx/Hna4mdsNXpzJ2UJlkyzcAwxJZLV6jMiZ2mHTQKU1MQJ5TlH+sPNoWMSin1n7WighjFiMU5dLVO54LLyGRWFLE1qWu7JiEuy6q9tqzZrFWRrau+/eOc1tPbe3O/Nzn3Hu/7Ph8Q7jnf8+Gczyff3Ge+93vyPUlVIUnq5TWznoAkafqMuyQ1ZNwlqSHjLkkNGXdJasi4S1JDxl2SGjLuktSQcZekhrbO6oW3bdtWu3btmtXLS9KG9Mgjj3y5quZWGjezuO/atYuTJ0/O6uUlaUNK8k9DxnlaRpIaMu6S1JBxl6SGjLskNWTcJakh4y5JDRl3SWrIuEtSQ8Zdkhqa2RWql2PXoU/O7LWf/NC7Z/bakjSUR+6S1JBxl6SGjLskNWTcJakh4y5JDRl3SWrIuEtSQ8Zdkhoy7pLUkHGXpIaMuyQ1ZNwlqSHjLkkNGXdJasi4S1JDg+KeZF+Ss0nOJTm0xOM7kzyY5LNJHk9y8/SnKkkaasW4J9kCHAZuAvYCtybZu2jYLwIPVNU1wAHgN6c9UUnScEOO3K8DzlXV+ap6Abgf2L9oTAHfNL79BuDp6U1RkrRaQ/6bvauBpybuLwBvXzTmLuDTSX4aeB1w41RmJ0m6JEOO3LPEtlp0/1bg96pqO3Az8JEkr3ruJAeTnExy8sKFC6ufrSRpkCFxXwB2TNzfzqtPu9wGPABQVX8HfCOwbfETVdWRqpqvqvm5ublLm7EkaUVD4n4C2JNkd5IrGL1henTRmH8Gvh8gyXcyiruH5pI0IyvGvapeBG4HjgNnGP2rmFNJ7k5yy3jYB4D3JXkM+Cjwk1W1+NSNJGmdDHlDlao6BhxbtO3OidungRumOzVJ0qXyClVJasi4S1JDxl2SGjLuktSQcZekhoy7JDVk3CWpIeMuSQ0Zd0lqyLhLUkPGXZIaMu6S1JBxl6SGjLskNWTcJakh4y5JDRl3SWrIuEtSQ8Zdkhoy7pLUkHGXpIaMuyQ1ZNwlqSHjLkkNGXdJasi4S1JDxl2SGjLuktSQcZekhoy7JDVk3CWpIeMuSQ0Zd0lqyLhLUkPGXZIaMu6S1JBxl6SGjLskNTQo7kn2JTmb5FySQ8uM+dEkp5OcSvKH052mJGk1tq40IMkW4DDwLmABOJHkaFWdnhizB/ggcENVPZfkTWs1YUnSyoYcuV8HnKuq81X1AnA/sH/RmPcBh6vqOYCqena605QkrcaQuF8NPDVxf2G8bdJbgbcm+dskDyXZN60JSpJWb8XTMkCW2FZLPM8e4J3AduBvkrytqv7tFU+UHAQOAuzcuXPVk5UkDTPkyH0B2DFxfzvw9BJjPlFVX6+qLwJnGcX+FarqSFXNV9X83Nzcpc5ZkrSCIXE/AexJsjvJFcAB4OiiMX8GfB9Akm2MTtOcn+ZEJUnDrRj3qnoRuB04DpwBHqiqU0nuTnLLeNhx4CtJTgMPAr9QVV9Zq0lLki5uyDl3quoYcGzRtjsnbhfw/vEvSdKMeYWqJDVk3CWpIeMuSQ0Zd0lqyLhLUkPGXZIaMu6S1JBxl6SGjLskNWTcJakh4y5JDRl3SWrIuEtSQ8Zdkhoy7pLUkHGXpIaMuyQ1ZNwlqSHjLkkNGXdJasi4S1JDxl2SGjLuktSQcZekhoy7JDVk3CWpIeMuSQ0Zd0lqyLhLUkPGXZIaMu6S1JBxl6SGjLskNWTcJakh4y5JDRl3SWrIuEtSQ8ZdkhoaFPck+5KcTXIuyaGLjHtPkkoyP70pSpJWa8W4J9kCHAZuAvYCtybZu8S4K4GfAR6e9iQlSasz5Mj9OuBcVZ2vqheA+4H9S4z7JeBe4KtTnJ8k6RIMifvVwFMT9xfG216W5BpgR1X9+cWeKMnBJCeTnLxw4cKqJytJGmZI3LPEtnr5weQ1wK8BH1jpiarqSFXNV9X83Nzc8FlKklZlSNwXgB0T97cDT0/cvxJ4G/BXSZ4ErgeO+qaqJM3OkLifAPYk2Z3kCuAAcPSlB6vq+araVlW7qmoX8BBwS1WdXJMZS5JWtGLcq+pF4HbgOHAGeKCqTiW5O8ktaz1BSdLqbR0yqKqOAccWbbtzmbHvvPxpSZIuh1eoSlJDxl2SGjLuktSQcZekhoy7JDVk3CWpIeMuSQ0Zd0lqyLhLUkPGXZIaMu6S1JBxl6SGjLskNWTcJakh4y5JDRl3SWrIuEtSQ8Zdkhoy7pLUkHGXpIaMuyQ1ZNwlqSHjLkkNGXdJasi4S1JDxl2SGjLuktSQcZekhoy7JDVk3CWpIeMuSQ0Zd0lqyLhLUkPGXZIaMu6S1JBxl6SGjLskNTQo7kn2JTmb5FySQ0s8/v4kp5M8nuQzSd48/alKkoZaMe5JtgCHgZuAvcCtSfYuGvZZYL6qvgv4GHDvtCcqSRpuyJH7dcC5qjpfVS8A9wP7JwdU1YNV9V/juw8B26c7TUnSagyJ+9XAUxP3F8bblnMb8KnLmZQk6fJsHTAmS2yrJQcmPw7MA9+7zOMHgYMAO3fuHDhFSdJqDTlyXwB2TNzfDjy9eFCSG4E7gFuq6mtLPVFVHamq+aqan5ubu5T5SpIGGBL3E8CeJLuTXAEcAI5ODkhyDfDbjML+7PSnKUlajRXjXlUvArcDx4EzwANVdSrJ3UluGQ/7ZeD1wB8neTTJ0WWeTpK0Doacc6eqjgHHFm27c+L2jVOelyTpMniFqiQ1ZNwlqSHjLkkNGXdJasi4S1JDxl2SGjLuktSQcZekhoy7JDVk3CWpIeMuSQ0Zd0lqyLhLUkPGXZIaMu6S1JBxl6SGjLskNWTcJakh4y5JDRl3SWrIuEtSQ8Zdkhoy7pLUkHGXpIaMuyQ1ZNwlqSHjLkkNGXdJamjrrCew0ew69MmZvO6TH3r3TF5X0sbkkbskNWTcJakh4y5JDRl3SWrIuEtSQ8Zdkhoy7pLUkHGXpIa8iGmDmNXFU+AFVNJGNOjIPcm+JGeTnEtyaInHX5vkj8aPP5xk17QnKkkabsUj9yRbgMPAu4AF4ESSo1V1emLYbcBzVfUdSQ4A9wA/thYT1vrzIxekjWfIaZnrgHNVdR4gyf3AfmAy7vuBu8a3Pwb8RpJUVU1xrtpk/EtFunRD4n418NTE/QXg7cuNqaoXkzwPfCvw5WlMUlpPs3x/Y1b8C62fIXHPEtsWH5EPGUOSg8DB8d3/THJ2wOsvto3N+ZfGZlz3ZlwzzGDduWc9X21Zm3F/X8qa3zxk0JC4LwA7Ju5vB55eZsxCkq3AG4B/XfxEVXUEODJkYstJcrKq5i/nOTaizbjuzbhmcN2znsd6Wss1D/nXMieAPUl2J7kCOAAcXTTmKPDe8e33AH/p+XZJmp0Vj9zH59BvB44DW4D7qupUkruBk1V1FPhd4CNJzjE6Yj+wlpOWJF3coIuYquoYcGzRtjsnbn8V+JHpTm1Zl3VaZwPbjOvejGsG172ZrNma49kTSerHz5aRpIY2VNxX+hiELpI8meRzSR5NcnK87Y1J/iLJ58dfv2XW87xcSe5L8mySJya2LbnOjPz6eN8/nuTa2c388iyz7ruS/Mt4nz+a5OaJxz44XvfZJD84m1lfniQ7kjyY5EySU0l+dry97f6+yJrXZ19X1Yb4xejN3C8AbwGuAB4D9s56Xmu01ieBbYu23QscGt8+BNwz63lOYZ3vAK4FnlhpncDNwKcYXVNxPfDwrOc/5XXfBfz8EmP3jv+svxbYPf4e2DLrNVzCmq8Crh3fvhL4x/Ha2u7vi6x5Xfb1Rjpyf/ljEKrqBeClj0HYLPYDHx7f/jDwQzOcy1RU1V/z6ushllvnfuD3a+Qh4JuTXLU+M52uZda9nP3A/VX1tar6InCO0ffChlJVz1TVP4xv/wdwhtGV7W3390XWvJyp7uuNFPelPgbhYr9RG1kBn07yyPiqXoBvq6pnYPSHBnjTzGa3tpZb52bY/7ePT0HcN3Hard26x58aew3wMJtkfy9aM6zDvt5IcR/0EQdN3FBV1wI3AT+V5B2zntD/A933/28B3w58N/AM8Cvj7a3WneT1wJ8AP1dV/36xoUts25DrXmLN67KvN1Lch3wMQgtV9fT467PAxxn9aPall34sHX99dnYzXFPLrbP1/q+qL1XVf1fV/wC/w//9ON5m3Um+gVHk/qCq/nS8ufX+XmrN67WvN1Lch3wMwoaX5HVJrnzpNvADwBO88iMe3gt8YjYzXHPLrfMo8BPjf0VxPfD8Sz/Od7DofPIPM9rnMFr3gYz+Q5zdwB7g79d7fpcrSRhdyX6mqn514qG2+3u5Na/bvp71O8qrfPf5ZkbvOH8BuGPW81mjNb6F0TvmjwGnXlono49Q/gzw+fHXN856rlNY60cZ/Vj6dUZHLbctt05GP7IeHu/7zwHzs57/lNf9kfG6Hh9/k181Mf6O8brPAjfNev6XuObvYXSK4XHg0fGvmzvv74useV32tVeoSlJDG+m0jCRpIOMuSQ0Zd0lqyLhLUkPGXZIaMu6S1JBxl6SGjLskNfS/CmmyQ7cnxbIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "weights = np.ones_like(degree) / float(n_nodes)\n",
    "plt.hist(degree, weights=weights);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is the average degree?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The average degree is 11.4592\n"
     ]
    }
   ],
   "source": [
    "avg_deg = np.mean(degree)\n",
    "print(\"The average degree is {:.4f}\".format(avg_deg))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 6\n",
    "\n",
    "Comment on the degree distribution of your network."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAECCAYAAADw0Rw8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAFsxJREFUeJzt3X+MHHd5x/HPc5cN3bOqOxvyB1nb1MSpXVM3tTglFEsViQAb6BHLhBIH/gBFsagEUlFlyVGRwh9FMbIq0agp0dFYKSrND9LI2MLIf9RBUSOo7GBHJE3cuJEgd5aIIblria/4bD/9Y28vc3szuzM7uze/3i8pcm52du57HvnZ733n+T6PubsAAOU1lPUAAACDRaAHgJIj0ANAyRHoAaDkCPQAUHIEegAoOQI9AJQcgR4ASq7vgd7M/sDMHjKzJ83sL/p9fQBAMrECvZkdMrPXzeyFtuM7zeysmZ0zs/2S5O4vufsXJf25pPH+DxkAkETcGf0jknYGD5jZsKQHJX1M0hZJe8xsy8Jrn5T075L+rW8jBQD0JFagd/dnJL3RdvhmSefc/VV3vyTpMUm3L5x/xN0/KOmz/RwsACC5a1K8tyHptcDXU5JuMbMPSdot6R2SjkW92cz2StorSatWrXr/5s2bUwwFAKrnueee+5W7X9ftvDSB3kKOubv/SNKPur3Z3SclTUrS+Pi4nzp1KsVQAKB6zOzncc5Lk3UzJWld4Ou1ks6nuB4AYADSBPqTkm40sw1mdq2kOyUdSXIBM5sws8nZ2dkUwwAAdBI3vfJRST+WtMnMpszsbne/LOlLko5LeknSE+7+YpJv7u5H3X3v6Oho0nEDAGKKtUbv7nsijh9ThweuAIDspXkYm7nDp6d18PhZnZ+Z0/Vjde3bsUm7tjWyHhYA5Eqmgd7MJiRNbNy4MfF7D5+e1r1P/Uxz81ckSdMzc7r3qZ9JUubBng8gAHmSaVGzNGv0B4+fXQzyLXPzV3Tw+Nl+Da8nrQ+g6Zk5ud7+ADp8ejrTcQGorsJWrzw/M5fo+ErJ6wcQgOoqbKC/fqye6PhKyesHEIDqyjTQp8mj37djk+q14SXH6rVh7duxqV/D60leP4AAVFdh1+h3bWvo/t1b1RiryyQ1xuq6f/fWzB965vUDCEB1FTq9cte2RuaBvV1rPGTdAMiLQgf6vMrjBxCA6irsGj0AIJ7CrtEDAOIpbHolACAeAj0AlByBHgBKjkAPACVH1g0AlBxZNwBQcizdAEDJsTO2QGhoAqAXBPqCyHNHLQD5xtJNQdDQBECvyLopCBqaAOgVWTcFQUMTAL1i6aYgaGgCoFc8jC0IGpoA6BWBvkBoaAKgFyzdAEDJEegBoOQI9ABQcgR6ACg5NkwBQMmxYQoASo6lGwAoOQI9AJQcG6YqhHr2QDUR6CuCevZAdbF0UxHUsweqi0BfEdSzB6qLQF8R1LMHqotAXxHUsweqi4exFUE9e6C6CPQVQj17oJqodQMAJUetGwAoOR7GAkDJEegBoOR4GItI1MYByoFAj1DUxgHKg6UbhKI2DlAeBHqEojYOUB4EeoSiNg5QHgR6hKI2DlAePIxFKGrjAOVBoEckauMA5cDSDQCUHIEeAEqOpRukxg5aIN8I9EiFHbRA/rF0g1TYQQvk30ACvZntMrNvm9n3zeyjg/geyAd20AL5FzvQm9khM3vdzF5oO77TzM6a2Tkz2y9J7n7Y3e+R9HlJn+nriJEr7KAF8i/JjP4RSTuDB8xsWNKDkj4maYukPWa2JXDKVxdeR0mxgxbIv9iB3t2fkfRG2+GbJZ1z91fd/ZKkxyTdbk3fkPRDd/9p/4aLvNm1raH7d29VY6wuk9QYq+v+3Vt5EAvkSNqsm4ak1wJfT0m6RdKXJX1Y0qiZbXT3h9rfaGZ7Je2VpPXr16ccBrLEDlog39IGegs55u7+gKQHOr3R3SclTUrS+Pi4pxwHACBC2qybKUnrAl+vlXQ+5TUBAH2UdkZ/UtKNZrZB0rSkOyXdFffNZjYhaWLjxo0ph4G8Y/cskJ0k6ZWPSvqxpE1mNmVmd7v7ZUlfknRc0kuSnnD3F+Ne092Puvve0dHRpONGgbR2z07PzMn19u7Zw6enF1/ffuCENuz/gbYfOLF4HEB/xJ7Ru/ueiOPHJB3r24hQOt12z3YqocBvAkB61LrBwHXaPZvmQwBAPJnWujGzCTObnJ2dzXIYGLBOu2fTfAgAiCfTQM8afTV02j3b64cAgPioXomB67R7ttcPAQDxsUaPFRG1e7ZbE/LgGr1EHR2gF5kGevLoIfX+IQAgHnPPvvrA+Pi4nzp1KuthAEChmNlz7j7e7TzW6AGg5Aj0AFBy5NEDQMmRRw8AJcfSDQCUHIEeAEqODVMoJKpaAvER6FE4rfr2VLUE4sl0w1RgZ+w9r7zySmbjQLFsP3BC0yGFzcbqNa16xzXM8lEZhdgwRdYNehFVvXJmbj6yixVQZTyMReHErV5J7XqgiUCPwgkrbRyF2vUAD2NRQGFVLS9euqw3L84vO5fa9QCBHgXVXtq4PRNHonY90EI9epRCVO16qZmlQyYOqox69CitqFl+q40hUHSFSK8EBung8bNLgrxEJg6qiUCP0orKuCETB1VDoEdpRWXckImDqiHQo7TC8u3JxEEVkV6J0grLxLl183U6ePysvvL4GbJwUBkEepRaMN+eqpeoKgI9KqNTFk7SQE89fBQJG6ZQGf3KwuE3AxQNZYpRGf3Kwon6zeAvHz+j7QdOUBoZuUPWDSqjX1k4nX4DoA4+8ohAj8rYta2h+3dvVWOsLpPUGKv3VA6h228A7L5F3vAwFpXSXvWyF/t2bFpWQ6cdu2+RJwR6IKFgfn5Y71qJ3bfIF5ZugB7s2tbQs/tv0zc/88fsvkXuMaMHUoiqg0+aJfKEevRAH7GRCispbj16Aj0qrT0w37r5Oj398oWeAnVYoxOT5Gpm+BD00W9xAz1LN6issB2u//yTXyy+nnTHa9hGqtY0it2zyBKBHpUVFpjbzc1f0V898by+8vgZjdZrMpNmLs6Hzva7pVT2WlcHSCvTrBszmzCzydnZ2SyHgYqKm+t+xV0uaWZuXm9enJcrfAdsnJRK8uuRBWrdoLLS5rq374ANK7HQ7+8J9II8elRWnMDcTXCGHiyxIDUfxAZ1yq8/fHpa2w+c0Ib9P6AwGvqONXpUVlQHqlbWzZCZrnTJSmufobc3OomTaknZYwwa6ZVAhLB0yaB6bbinomjtth84EVpKoTFW17P7b1scC/n5aEd6JZBS+4y/W9ZNr7o1RGHGj7QI9EAH/ah2GRQ2M79+rB46o28tC/WzBSKqiYexwAppzcynZ+aWpGjeuvm6joXR+tUCEdVFoAdWSNTM/OmXL0Q2RDl8elpD1p6/09QtVZNMHrSwdAOskE4z87AlotZvAGGZP91KIbOujyBm9MCAtWbWUfltUTPzqBINw2Zds306reujepjRAwMUJ0UzamYe9RvAVfeus3LW9RHEjB4YoE6F07o1J4+a6ccpoxB1zpAZa/YVRKAHBihqBm2Snt1/W8eZeViJhrhtCqPKO7QKtIUVZUN5EeiBAUozKw/WzmnPxkn63uGQzB3W7KuDNXpggPbt2LRsjT5J8/CoDVtxSiIE37th/w9Cr8+afTUQ6IEBGkTz8F5SJ7vtvkW59T3Qm9l7Jf21pFF3v6Pf1weKpt9lFHopiZD2NwsUW6w1ejM7ZGavm9kLbcd3mtlZMztnZvslyd1fdfe7BzFYANHLLdMzc5E7YdOs96P44s7oH5H095K+0zpgZsOSHpT0EUlTkk6a2RF3/89+DxLA26KWYSQtHg9bzun3bxYojlgzend/RtIbbYdvlnRuYQZ/SdJjkm7v8/gAtNm3Y9Oy7lVhyKpBS5r0yoak1wJfT0lqmNk7zewhSdvM7N6oN5vZXjM7ZWanLly4kGIYQLXs2taILKfQjqwaSOkexoZNKtzdfy3pi93e7O6TkialZoepFOMAKqfRYfkmiKwaSOlm9FOS1gW+XivpfLrhAIgjTmNzsmrQkmZGf1LSjWa2QdK0pDsl3ZXkAmY2IWli48aNKYYBVE+3xuattodfefyMvnbkxVgtEJP0paWHbbHEag5uZo9K+pCkd0n6paT73P1hM/u4pG9KGpZ0yN2/3ssgaA4O9E8vTc3D3hPV/DzJuRisuM3B42bd7HH3d7t7zd3XuvvDC8ePufvvu/sNvQZ5AP3VqWKmFJ6Nk6R+PbXuiyfTomZmNmFmk7Ozs1kOAyiVOJk27eckqV9PrfviybTWjbsflXR0fHz8nizHAWQl6Vp32PnS0rX6sZGa3rw43/H7tmfjJKmFQ92c4qFMMZCR1lr39MxcrBrxYefv+97z2vfk80uO/eb/Lqs2HL2lKiwbJ0nt+zR18pENAj2QkaRr3WHnz191zV/xZcdWXXvNYl2bsXpNq0dqHWvcJKmFQ92c4sl06Yb0SlRZ0rXuJGvgs3PzOnPfRxONJ0ktHOrmFEumM3p3P+rue0dHR7McBpCJpN2nkqyBs16OIBqPABk4fHpab/328rLjrbXusIeuYTXla0MmmZYs36RdLz98elpfO/KiZubefqC7eqSm+ybeJ6m/TVSwMmJtmBo0NkyhSqI2NAWDadSGJGl5oA071mvwPXx6Wvu+97zmry6PC0MmDQ/Zsg8V1uezE3fDFDN6YIVFbWgaufYa7drW0PYDJyIf0j67/7bIB6T9GltYkJekqy5dbXvw262zFfKBDVPACuv2sDXLDUm9fA82SuUfD2OBFdbtYWs/Hsb2qpfvwYPf/COPHlhh3TYcZbkhad+OTc0HvCGGTMs2YrFRqhhYowdWWFiJ4eAD1G6v90NU6YXW9yDrplzIugEqhjLD5dHXMsUAyoMyw9VD1g1QMZQZrh6yboCKyTKrB9lg6QaoGMoMVw8PY4EKaq9nM2TNna9jC03Fg43Epea6/vTMnIbNdMV98c9GW+ZNVDZPPxqs8KB4OUogAOjot5evLv5/q+pBMKWy1dgkWDTtii/9s9UspSWYzdN67dTP39C/Pje97LgUXrqhPSuo2/nojqUboIK6NRBvCWts0q6VsROVzfPof7yWusEKWUHp0HgEqKB+Z9h0ut6ViOXhQTZewVJk3QAV1O8Mm+vH6pHXHLbwkgqDbLyCpVi6ASooLPMmTG3IOjYal97O2InK5tlzy7pEWT5kBfUfD2OBCgrW0wlm06TNummd254tM/6eNbGzaFai1k/VkF4JAAVFrRsAgCSWboBSibNhabTD8kySpZLWNadn5mSSWmsDq0dq+sQfvVtPv3whclmoNYY3L84vee9IbUjvqA13HVsv4+309xP39aJi6QYoiajyw596f2PJhqV2tSFbsimq9b5OZYujGpz3W9jYehlv1JiD7yti+WaWboCKSbJhKShsU1S3DUpxN1ylFTa2XsYrdd+IVeaNWpQpBkoiakNR1IalXq/X7bWsdBtTnpuyDxobpoCSSLphqdfrdXstK93GlOem7IPG0g1QEkk2LAWFbYrqtkEp7oartMLG1st4pXw3ZR80sm6Akui00Si4YakfWTftG66KkHWTh6bsWSHrBgAKiqwbAIAkAj0AlB6BHgBKjkAPACVH1g2QI0WqtdLeYLyVMfPmxfll5YyDmTeXLl/Rxflmv9pWU/J6bUhz81eXXH+k1pyHXmw7LjUze7a8+3f1k1ffXLz+nlvW6W92bV3yd/g7tSH99vLVxZ64kkIbmgd/jtb1W5lDZaiLQ9YNkBNFqrVy+PS09n3vec1fzT5+BG2/YY1++ovZruUZWn+vkmL/HHmsi0PWDVAwRaq1cvD42dwFeUl69r/fiFWDJ9jQPO7PUeS6ODQHB3KiSLVW8jimpHr5GYpaF4daN0BOFKnWSh7HlFSnhuad3hP8M+r1vGHpBsiJItVa2bdjU7MufM5sv2FNrBo8wYbmcX+OItfFIdADObFrW0P3796qxlhdpmZ2SB4fxErNsR789E0aq9cWj43UhrR6pPl1q2Jm+59j9dpiNo3UzLqRmlk37UZqQ0vODVo9UtP2G9Ysuf7nPrBe373nT5b8HdZrQ2qP48G/17Cfo3X9z31gfeS9KNK9ksi6AYDCIusGACCJQA8ApUegB4CSI9ADQMkR6AGg5Aj0AFByBHoAKDkCPQCUHIEeAEqOQA8AJUegB4CSI9ADQMn1vfGIma2S9A+SLkn6kbt/t9/fAwAQX6xAb2aHJP2ZpNfd/Q8Dx3dK+jtJw5L+0d0PSNot6Ul3P2pmj0si0AMV06lxdvC10XpNZtLMxfnF86Rmq77pmbnIpuJBtSGp/bBJ+uANa/Ti+f8NbV4eZtW1w3rr0pXF72mSwmr7Dpl01y3rNf6eNV3H2Wp+3hir69bN1+nply9oemZu2bVXj9R038T7BlbmOFaZYjP7U0m/kfSdVqA3s2FJ/yXpI5KmJJ2UtEfS7ZJ+6O5nzOxf3P2ubtenTDFQHp0aZ0ta9lpQbcgkk+avZF8+vZvhIdOVPvbNrQ2bDt5xU6JgH7dMcawZvbs/Y2a/13b4Zknn3P3VhW/4mJpBfkrSWklnxDMAoHK6Nc7u1Lw7jw3Ho/QzyEvND7eDx88OZFafJhA3JL0W+Hpq4dhTkj5lZt+SdDTqzWa218xOmdmpCxcupBgGgDzp1Dg7r82z82JQfz9pHsaGNVp0d39L0he6vdndJyVNSs2lmxTjAJAj14/VNR0SsFqNs8NeQ9OgmounmdFPSVoX+HqtpPPphgOg6Do1zg57Lag2ZKoN56/peJjhPjdHrw3bwJqLpwn0JyXdaGYbzOxaSXdKOpLkAmY2YWaTs7OzKYYBIE86Nc5uf22sXtPqkdrieQc/fZMO3nGTGgsz26im4kFhh03S9hvWRDYvD7Pq2uEl3zMqjA+Z9LkPrNfffrr7OFufBY2x+mKz8bBrrx6pJX4Qm0TcrJtHJX1I0rsk/VLSfe7+sJl9XNI31UyvPOTuX+9lEGTdAEBy/c662RNx/JikYwnHBgBYQZmmP7J0AwCDl2mgd/ej7r53dHQ0y2EAQKmxoQkASo5ADwAl1/fqlUmY2YSkCUn/Y2avLBwelRS1aB/12rsk/ar/I0yt08+S5XWTvj/u+XHO6+X+dnqNez/Y93Pvk1vJe/+eWO9091z9J2ky6WuSTmU97qQ/S5bXTfr+uOfHOa+X+8u9595z79NdN49LN5H1cbq8lkeDGm/a6yZ9f9zz45zX6/3l3vfnutz7wcvdvY+1YSrvzOyUx9g0gPLh3lcX9z6+PM7oezGZ9QCQGe59dXHvYyrFjB4AEK0sM3oAQAQCPQCUHIEeAEqudIHezFaZ2T+Z2bfN7LNZjwcrx8zea2YPm9mTWY8FK8vMdi38m/++mX006/HkTSECvZkdMrPXzeyFtuM7zeysmZ0zs/0Lh3dLetLd75H0yRUfLPoqyb1391fd/e5sRop+S3jvDy/8m/+8pM9kMNxcK0Sgl/SIpJ3BA2Y2LOlBSR+TtEXSHjPbomZLw1bT8uh28yiKRxT/3qNcHlHye//VhdcRUIhA7+7PSHqj7fDNks4tzOIuSXpM0u1q9rJdu3BOIX4+REt471EiSe69NX1D0g/d/acrPda8K3IgbOjtmbvUDPANSU9J+pSZfUvF2zqNeELvvZm908wekrTNzO7NZmgYsKh/91+W9GFJd5jZF7MYWJ5lWr0ypbDeve7ub0n6wkoPBisq6t7/WhL/yMst6t4/IOmBlR5MURR5Rj8laV3g67WSzmc0Fqws7n11ce97UORAf1LSjWa2wcyulXSnpCMZjwkrg3tfXdz7HhQi0JvZo5J+LGmTmU2Z2d3uflnSlyQdl/SSpCfc/cUsx4n+495XF/e+fyhqBgAlV4gZPQCgdwR6ACg5Aj0AlByBHgBKjkAPACVHoAeAkiPQA0DJEegBoOQI9ABQcv8PRDgwf76LLLYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "test = np.unique(degree,return_counts=True)\n",
    "ax = plt.gca()\n",
    "ax.scatter(test[0],test[1])\n",
    "ax.set_yscale(\"log\")\n",
    "ax.set_xscale(\"log\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our degree distribution follows a power law distribution, hence our network seems to be scale free as we saw in the lecture."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 7\n",
    "\n",
    "Write a function that takes as input the adjacency matrix of a graph and determines whether the graph is connected or not."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def BFS(adjacency, labels, state):\n",
    "    \"\"\"\n",
    "    return a component with an array of updated labels (the visited nodes during the BFS)\n",
    "    for a given adjacency matrix\n",
    "    \n",
    "    :param adjacency: The adjacency matrix where to find the component\n",
    "    :param labels: An array of labels (0 : the node is not yet explored, 1: it is explored)\n",
    "    :param state: The # of the component we are looking for\n",
    "    :return: updated labels array and the component found\n",
    "    \"\"\"\n",
    "    queue = []\n",
    "    # current node is the first one with a label to 0\n",
    "    current_node = np.argwhere(labels == 0).flatten()[0]\n",
    "    labels[current_node] = state\n",
    "    queue.append(current_node)\n",
    "\n",
    "    current_component = []\n",
    "    current_component.append(current_node)\n",
    "    \n",
    "    while len(queue) != 0:\n",
    "        # all the weight of the other nodes for a given one\n",
    "        all_nodes = adjacency_sym[current_node]\n",
    "        # all the nodes reachable from the current one that are not yet labeled\n",
    "        neighbours = np.argwhere((all_nodes > 0) & (labels == 0)).flatten()\n",
    "        # add those nodes to the queue and to the component\n",
    "        queue += list(neighbours)\n",
    "        current_component += list(neighbours)\n",
    "        \n",
    "        for i in neighbours:\n",
    "            # we update the labels array\n",
    "            labels[i] = state\n",
    "        if len(queue) > 1:\n",
    "            # and we update the queue and the current node\n",
    "            queue = queue[1:]\n",
    "            current_node = queue[0]\n",
    "        else :\n",
    "            queue = []\n",
    "    \n",
    "    return np.array(labels), current_component    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def connected_graph(adjacency):\n",
    "    \"\"\"Determines whether a graph is connected.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    adjacency: numpy array\n",
    "        The (weighted) adjacency matrix of a graph.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    bool\n",
    "        True if the graph is connected, False otherwise.\n",
    "    \"\"\"\n",
    "    #Run the BFS, find a component and see if all the nodes are in it. If so, the graph is connected.\n",
    "    first_labels = np.zeros(n_nodes, dtype=int)\n",
    "    labels, component = BFS(adjacency, first_labels, 1)\n",
    "    return labels.sum() == n_nodes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Is your graph connected? Run the ``connected_graph`` function to determine your answer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Is our graph connected ? False\n"
     ]
    }
   ],
   "source": [
    "connected = connected_graph(adjacency)\n",
    "print(\"Is our graph connected ? {}\".format(connected))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 8\n",
    "\n",
    "Write a function that extracts the connected components of a graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_components(adjacency):\n",
    "    \"\"\"Find the connected components of a graph.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    adjacency: numpy array\n",
    "        The (weighted) adjacency matrix of a graph.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    list of numpy arrays\n",
    "        A list of adjacency matrices, one per connected component.\n",
    "    \"\"\"\n",
    "    \n",
    "    #Find the first component\n",
    "    components = []\n",
    "    first_labels = np.zeros(n_nodes, dtype=int)\n",
    "    labels, component = BFS(adjacency, first_labels, 1)\n",
    "    components.append(component)\n",
    "    current_state = 2\n",
    "    \n",
    "    #Redo BFS while we haven't found all the components\n",
    "    while (labels > 0).sum() != n_nodes:\n",
    "        labels, component = BFS(adjacency, labels, current_state)\n",
    "        components.append(component)\n",
    "        current_state += 1\n",
    "        \n",
    "    return np.array(components)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How many connected components is your network composed of? What is the size of the largest connected component? Run the ``find_components`` function to determine your answer. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of connected components : 7\n",
      "Size of largest connected component : 3304\n"
     ]
    }
   ],
   "source": [
    "components = find_components(adjacency_sym)\n",
    "print(\"Number of connected components : {}\".format(len(components)))\n",
    "size_compo = [len(compo) for compo in components]\n",
    "print(\"Size of largest connected component : {}\".format(np.max(size_compo)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 9\n",
    "\n",
    "Write a function that takes as input the adjacency matrix and a node (`source`) and returns the length of the shortest path between that node and all nodes in the graph using Dijkstra's algorithm. **For the purposes of this assignment we are interested in the hop distance between nodes, not in the sum of weights. **\n",
    "\n",
    "Hint: You might want to mask the adjacency matrix in the function ``compute_shortest_path_lengths`` in order to make sure you obtain a binary adjacency matrix. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Find all the neighbours of a given node\n",
    "def neighbours(node, adjacency_sym):\n",
    "    n = adjacency_sym[node]\n",
    "    neighbours = np.argwhere(n != 0).flatten()\n",
    "    return neighbours"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_shortest_path_lengths(adjacency, source):\n",
    "    \"\"\"Compute the shortest path length between a source node and all nodes.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    adjacency: numpy array\n",
    "        The (weighted) adjacency matrix of a graph.\n",
    "    source: int\n",
    "        The source node. A number between 0 and n_nodes-1.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    list of ints\n",
    "        The length of the shortest path from source to all nodes. Returned list should be of length n_nodes.\n",
    "    \"\"\"   \n",
    "    \n",
    "    adjacency[adjacency != 0] = 1\n",
    "    shortest_path_lengths = np.ones(adjacency.shape[0]) * np.inf\n",
    "    shortest_path_lengths[source] = 0\n",
    "    visited = [source]\n",
    "    queue = [source]\n",
    "    while queue:\n",
    "        node = queue[0]\n",
    "        queue = queue[1:]\n",
    "        neighbors = neighbours(node, adjacency)\n",
    "        neighbors = np.setdiff1d(neighbors,visited).tolist()\n",
    "        neighbors = np.setdiff1d(neighbors,queue).tolist()\n",
    "        queue += neighbors\n",
    "        visited += neighbors\n",
    "        shortest_path_lengths[neighbors] = shortest_path_lengths[node] + 1\n",
    "    return shortest_path_lengths"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 10\n",
    "\n",
    "The diameter of the graph is the length of the longest shortest path between any pair of nodes. Use the above developed function to compute the diameter of the graph (or the diameter of the largest connected component of the graph if the graph is not connected). If your graph (or largest connected component) is very large, computing the diameter will take very long. In that case downsample your graph so that it has 1.000 nodes. There are many ways to reduce the size of a graph. For the purposes of this milestone you can chose to randomly select 1.000 nodes. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The diameter of the graph is 5.0\n"
     ]
    }
   ],
   "source": [
    "max_component = components[np.argmax(size_compo)]\n",
    "adjacency_max = adjacency_sym[max_component, :]\n",
    "adjacency_max = adjacency_max[:, max_component]\n",
    "\n",
    "longest = []\n",
    "a = adjacency_max[:1000,:1000]\n",
    "\n",
    "for node in range(len(a)):\n",
    "    short = compute_shortest_path_lengths(a, node)\n",
    "    longest.append(max(short[ short < np.inf ]))\n",
    "    \n",
    "print(\"The diameter of the graph is {}\".format(max(longest)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 11\n",
    "\n",
    "Write a function that takes as input the adjacency matrix, a path length, and two nodes (`source` and `target`), and returns the number of paths of the given length between them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_paths(adjacency, source, target, length):\n",
    "    \"\"\"Compute the number of paths of a given length between a source and target node.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    adjacency: numpy array\n",
    "        The (weighted) adjacency matrix of a graph.\n",
    "    source: int\n",
    "        The source node. A number between 0 and n_nodes-1.\n",
    "    target: int\n",
    "        The target node. A number between 0 and n_nodes-1.\n",
    "    length: int\n",
    "        The path length to be considered.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    int\n",
    "        The number of paths.\n",
    "    \"\"\"\n",
    "    \n",
    "    adjacency[adjacency != 0] = 1\n",
    "    adjacency = np.linalg.matrix_power(adjacency, length)\n",
    "    return adjacency[source, target]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test your function on 5 pairs of nodes, with different lengths."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "0\n",
      "0\n",
      "1\n",
      "0\n"
     ]
    }
   ],
   "source": [
    "print(compute_paths(adjacency_sym, 0, 10, 1))\n",
    "print(compute_paths(adjacency_sym, 0, 10, 2))\n",
    "print(compute_paths(adjacency_sym, 0, 10, 3))\n",
    "print(compute_paths(adjacency_sym, 23, 67, 2))\n",
    "print(compute_paths(adjacency_sym, 15, 93, 4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 12\n",
    "\n",
    "How many paths of length 3 are there in your graph? Hint: calling the `compute_paths` function on every pair of node is not an efficient way to do it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of path of length 3: 158232235\n"
     ]
    }
   ],
   "source": [
    "# We could have used np.linalg.matrix_power(a, 3), but for performance reasons we prefered to make \n",
    "#the multiplication explicitely.\n",
    "a = adjacency_sym.copy()\n",
    "a[a != 0] = 1\n",
    "a = a @ a @ a\n",
    "print(\"Number of path of length 3: {}\".format(a.sum()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 13\n",
    "\n",
    "Write a function that takes as input the adjacency matrix of your graph (or of the largest connected component of your graph) and a node and returns the clustering coefficient of that node. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_clustering_coefficient(adjacency, node):\n",
    "    \"\"\"Compute the clustering coefficient of a node.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    adjacency: numpy array\n",
    "        The (weighted) adjacency matrix of a graph.\n",
    "    node: int\n",
    "        The node whose clustering coefficient will be computed. A number between 0 and n_nodes-1.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    float\n",
    "        The clustering coefficient of the node. A number between 0 and 1.\n",
    "    \"\"\"\n",
    "    \n",
    "    adjacency[adjacency >1]=1\n",
    "    neighbors = adjacency[node]\n",
    "    indices = np.argwhere(neighbors > 0).flatten()\n",
    "    \n",
    "    if len(indices) <2:\n",
    "        return 0\n",
    "    \n",
    "    #Take the neighbors of the neighbors, and find wich one are linked together\n",
    "    neiofnei = adjacency[indices]\n",
    "    test = neiofnei * neighbors\n",
    "    count =  sum(test.flatten())\n",
    "    #Compute the clustering coefficient for the node. Since each link is counted twice, we don't multiply the formula by 2.\n",
    "    clustering_coefficient = count / (len(indices)*(len(indices)-1))    \n",
    "    \n",
    "    return clustering_coefficient"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 14\n",
    "\n",
    "What is the average clustering coefficient of your graph (or of the largest connected component of your graph if your graph is disconnected)? Use the function ``compute_clustering_coefficient`` to determine your answer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adjacency_component = adjacency_sym[components[0]]\n",
    "adjacency_component = adjacency_component.T[components[0]].T\n",
    "N = len(components[0])\n",
    "count = 0\n",
    "for i in range(N):\n",
    "    count+=compute_clustering_coefficient(adjacency_component,i)\n",
    "\n",
    "avg_coeff = count/N\n",
    "\n",
    "print(\"The average coefficient of our largest connected component is : {}\".format(str(avg_coeff)))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}