class12_pagerank_python3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CS446/546 - Class Session 12 - Pagerank/Katz/Eigenvector centrality\n",
    "\n",
    "In this class session we are going to compute the outgoing-edge PageRank centrality of each gene (vertex) in a human gene regulatory network (a directed graph) from a landmark paper on human gene regulation (Neph et al., Cell, volume 150, pages 1274-1286, 2012; see PDF on Canvas)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using Pandas `read_csv`, read in the ifle `shared/neph_gene_network.txt`, which has two columns of text (first column is the regulator gene, second column is the target gene), into a data frame. The file has no header and is tab-delimited.  Assign the column names of the dataframe to be `regulator` and `target`, respectively."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's load the Python packages that we will need for this exercise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas\n",
    "import igraph\n",
    "import numpy\n",
    "import matplotlib.pyplot\n",
    "import random"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using `pandas.read_csv`, read the file `shared/neph_gene_network.txt`; name the two columns of the resulting data frame, `regulator` and `target`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "edge_list_neph = pandas.read_csv(\"shared/neph_gene_network.txt\",\n",
    "                                 sep=\"\\t\", \n",
    "                                 names=[\"regulator\",\"target\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the edge-list data into a directed `igraph.Graph` object `neph_graph`, using `igraph.Graph.TupleList`. *Make sure to reverse the columns of the data frame when you input the data frame into `Graph.TupleList` since we want the outgoing pagerank centrality not incoming pagerank centrality.* Print out a summary of the graph, using the `igraph.Graph.summary` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'IGRAPH DN-- 538 47945 -- \\n+ attr: name (v)'"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "neph_graph = igraph.Graph.TupleList(edge_list_neph[[\"target\",\"regulator\"]].values.tolist(), directed=True)\n",
    "neph_graph.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the pagerank centrality measures of all vertices, using `igraph.Graph.pagerank`.  Use the resulting object to initialize a `numpy.array`, `pageranks`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "pageranks = numpy.array(neph_graph.pagerank())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which vertex has highest pagerank centrality in the gene regulatory network, and what is its pagerank centrality value? (think `numpy.max` and `numpy.argmax`). Get a `VertexSet` sequence using the `igraph.Graph.vs` property, and then index into that sequence using Python indexing:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.00749112688039\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "igraph.Vertex(<igraph.Graph object at 0x7f6dbbb9be58>, 122, {'name': 'FOS'})"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(numpy.max(pageranks))\n",
    "neph_graph.vs[numpy.argmax(pageranks)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the in-degree of all vertices in the graph, and scatter plot `log(degree)` vs. `log(pagerank)`.  (Do you see why we used `in` here?  Note the column swapping we did earlier). Note-- you will have to eliminate one vertex that has zero in-degree."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAikAAAGBCAYAAACjNCEAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3X18XGWd///XlXB/l4ZUq98FVqBJvuANbVqWFnoDMSFt\nWCvuT8EpzfpdFAR6g3WpKN5gFVelFRcKIrLiD8kyi/pYv9u1adPGQkuXFkrSCkK/TNrityyrWJq2\ny51IJ9f3j+uczJnJTDKZzGROJu/n4zGPzJw558w1c5rMp9f1uT6XsdYiIiIiEjZlxW6AiIiISDoK\nUkRERCSUFKSIiIhIKClIERERkVBSkCIiIiKhpCBFREREQklBioiIiISSghQREREJJQUpIiIiEkoK\nUkRERCSUFKSIiIhIKJVkkGKMOc0Ys9kY86wx5jljzFXFbpOIiIgMjSnFBQaNMROAU621u4wx7wJ2\nAv/TWvtakZsmIiIiWTqq2A0oBGvtK8Ar3v39xpjfA+MBBSkiIiKjREkO9wQZY6YCx1hrXyx2W0RE\nRCR7oQhSjDEzjTGrjTEvG2N6jTHz0uxzgzFmrzHmLWPMdmPMjCzOWwU8CFxTiHaLiIhI4YQiSAFO\nxOWN3AD0S5IxxlwJfA+4GagF1gNrjTGnBfZZbIzZYYzpMsYcZYw5BvhX4FvW2idH4k2IiIhI/oQu\ncdYY0wtcbq1dHdi2Ddhirb0psG0H0Gat/XKG80SBXdbab2TxmlVAE/A74E/DewciIiJjynHA+4B2\na+2BfJ449ImzxpijgSlAarDxa+DCDMdcBHwCeMYY8zFc70yLtfa5DC/TBPxzflosIiIyJl0FPJzP\nE4Y+SMHNyikH/pCy/Q/Ae9IdYK39D4b23n4H0NrayjnnnJNDE/Nr6dKlfP/73w/F+YZybDb7DrbP\nQM9nei7d9nx/hsNRiLbkes6hHlfIa5qv7cUQlt9RXc/8CMv1HOqxYfmbu2vXLhYsWADed2k+jYYg\nJZN8jlP9CeCcc86hrq4uj6fNTUVFRV7bMZzzDeXYbPYdbJ+Bns/0XLrt+f4Mh6MQbcn1nEM9rpDX\nNF/biyEsv6O6nvkRlus51GPD9DfXk/d0ibAkzg7kVSBO/16T99K/d6UkRCKR0JxvKMdms+9g+wz0\nfKbn8v155Vsh2pfrOYd6XCGv6Wi9nhCe31Fdz/wIy/Uc6rFj4W/uaE6c7QTWZUqcHeJr1gGdnZ2d\noYnsZXjmzZvH6tWrB99RRg1d09Ki61k6urq6mDJlCsAUa21XPs8diuEeY8yJwETAeJvOMsacB/RY\na18C7gB+Yox5EngKV/ekGri8GO0VERGRwgtFkAJMBR7F5ZlYXE0UcIXYrrbW/swYcypwO27Y57fA\nXC+AEelnNHQxy9DompYWXU/JRuiGe4rBH+6ZNWsWFRUVRCIR/QKJiIgMIBqNEo1GOXz4MJs3b4YC\nDPcoSEE5KSIiIrkqZE7KaJjdIyIiImOQghQREREJJQUpIiIiEkoKUkRERCSUFKSIiIhIKClIERER\nkVAKSzG3UFi6dKnqpIiIjHGxWIw9e/YwceJEqquri92c0ArWSSkU1UlBdVJERMKkWEFCT08P8+e3\n0N7e1retqamZaLSVysrKEWvHaKM6KSIiUvJ6enqYM+cyamtraW5upqamhjlzLuPgwYMj8vrz57fQ\n0bENaAX2Aa10dGwjElkwIq8v/SlIERGRUChmkBCLxWhvbyMevwu4CjgduIp4/E7a29vo7u4ueBuk\nPwUpIiJSdMUOEvbs2ePdm5XyzGwAdu/eXdDXl/QUpIiISNEVO0g4++yzvXubU57ZBMDEiRML+vqS\nnoIUEREpumIHCTU1NTQ1NVNevgQ33PQS0Ep5+Y00NTVrlk+RKEgREZGiC0OQEI220tAwDWgBzgBa\naGiYRjTaWvDXlvRUJ0VEREIhGm0lEllAe3tL37aGhuYRCxIqKytZt24N3d3d7N69W3VSQkB1UkjU\nSZk1a5aKuYmIFJmChNEhWMxt8+bNUIA6KQpSUDE3ERGRXKmYm4iIiIw5ClJEREQklBSkiIiISCgp\nSBEREZFQUpAiIiIioaQgRUREREJJQYqIiIiEkoIUERERCSWVxQ9YunSpKs6KiIhkIVhxtlBUcRZV\nnBURKbZYLMaePXtUCn8UUsVZEREpST09PcyZcxm1tbU0NzdTU1PDnDmXcfDgwWI3TUJAQYqIiBTN\n/PktdHRsA1qBfUArHR3biEQWJO0Xi8VYu3Yt3d3dxWimFIlyUkREpChisRjt7W24AOUqb+tVxOOW\n9vYWuru7qaqqYv78Fm8/p6mpmWi0lcrKyiG/noaURhf1pIiISFHs2bPHuzcr5ZnZAOzevTvrnpaB\naEhp9FKQIiIiRXH22Wd79zanPLMJgPLyctrb24jH78L1tJyO62m5k/b2tkGHfvwhossv/5thBzpS\nHBruERGRoqipqaGpqZmOjiXE4xbXg7KJ8vIbaWhoJh6Pe3tm7mlJN2zT09PTb4hooCGl4Dk0JBQu\n6kkREZGiiUZbaWiYBrQAZwAtNDRMIxptHbSnZeLEiWnPmTxE9KC3NXOgAxoSCisFKSIiUjSVlZWs\nW7eGWCxGW1sbsViMdevWUFlZ2dfTUl6+BBdwvAS0Ul5+I01NzWl7Ovxk3MQQ0TTvmYEDnXzkvkj+\nabhHRESKrrq6Om3QEY22EoksoL29pW9bQ4Ob3ZNO/2TcGqAZWAL0H1Kqrq7OapaRhn6KQ0GKiIiE\nlt/T0t3dze7duwfNFUkeIvIDjlagHjek5AQDnU2bNnlbh5b7IoWnIEVEREIvU09LqszJuPuYPn02\nt9xyc1+g4+ehJBJsg4ENDJb7IoWnICVACwyKiIx+6YaIpk+fzerVv0wqAJech/IAsJhMQ0LSnxYY\nHCFaYFBEpLT09PRw+eV/w+OPb+rbFqxUG4vFqK2tJZGHchBYAAy/su1YowUGRUREAgZby2f+/Bae\neOJZMs3W6Z9gWwmswR/iAVi16h8VoBSZghQRESm6bBcQzKaeSf9pyP0r1WauwbKv755fQ0WKR0GK\niIgUzVCLqA1Uz8QPdDZv9gOPzLN1ampqmDFjNrCQYA0WuBGYBChhNgyUOCsiIkWTHHTMAjbT0bGE\nSGQB69atSdp3sHomtbVtJPsZ8PeBx8mzdVav/iXV1edw4EBLYJ9JlJX9jsZGJcyGgXpSRESkKLIZ\nlgkabNVkWISrh+K7CZgMuNyU1Eq1lZWVdHfv8npUfDtpbLwwY7E4GVnqSRERkaIYLOhILaKWvlAb\nJJJdn8UPSPxeGTecMwnoTVuptrKykscffyzrYnEyshSkiIhIUZSV+Z352RVRy1SozfWgTPPuJw8F\nubonLaxfv57GxsaMbcm2WJyMLA33iIjIiPKTZefMmYP7GkpOXh1oAcF0qybDa8D53h7pe2WOHDlS\niLciBaYgRURERlRysuxO4EyCQUdDwzSi0da005L379/PjTcuYv369bS1tfHUU09RVfUu4CfeHgOv\ndiyji4Z7RERkxKSfobMDWAksY/369UyZMsUra5+YrXPJJQ0YY9i4cUPftqamZt555x0OHnwbmAg8\ng+uVUWn7UqEgRURERkzmZNkrgWUcOXIk7bTkRx/9NMYcT+pU5Xj8EG4Wz4nAhcC/kWm1Yxl9FKSI\niMiIGWyGTnl5eZqelvOBt7H2xyTXR3kN13Nye+A8fhaD4R/+4Vt86Utfyv+bkBGjnBQRERkx/gyd\n8vIlpEuWjcfj3p7BnpZMvS8/B04mWH3WPS4DLF/5ytd5+umnC/VWZAQoSAlYunQp8+bNIxqNFrsp\nIiIlK90MnfPOO4vbblueYU2ddNtiwEbgHoKF4OBuoBdYSW/v8Zx//gUDltmX3EWjUebNm8fSpUsL\n9hrGWluwk48Wxpg6oLOzs5O6urpiN0dEpKhisRh79uwpeGGz7du3c911C+nq2t63zU+G3bSpk3j8\nThK1UD6DMcdj7Spv2124ZNt9uADF9xIu8GkDDgAtlJVV0Nh4Ub8y+5IfXV1dTJkyBWCKtbYrn+dW\nT4qIiABDX+xvuL761a/zm9/sIXWxwAMHDjB9+gcJ9rTU18/ikkvOD2xb6Z0l/ZRjN9vH1Ujp7f14\n2jL7En5KnBUREWBoi/0N10CLBf7mN252zsyZs1m06AYmT57c16Pjl69/+eWXueaa60idcuxWMW4G\nqr1zA/wYKGPHjh2aijzKqCdFRESGvNjfcA2+WOAynnjiWR544MGkwKK6upqzzz6bvXv3AnFSC8G5\nn9/BBSiLgfNwvS4nc/fdP8jre5DCU5AiIiJZLfaXT+kTZCExXHNNvwApOBz17W9/29tvNS6J9hHc\nVOWdwIdwAct/A7/B1VE5k8cf36Qhn1FGQYqIiAwaNOS7rHymqcjJwzUuQIpGo3R3d6cMR00LtLca\nuAJ4ClgBGOAU4KckpibvA8ryHmxJYSlIERGRQeuXFCKXI/1igdO81+8B5gFw6623UlNT4w1H3YLr\nMdkGTAJSg5zluByVH5A8NflOoJejjlIq5miiIEVERID0QYO/2F8hVFZWsm7dGmKxGHV1UykrqwAi\nwOvAh4EXgR8C9YGjbgL+2rt/DZA8Cwje8J4beDXkdIsXSvgoSBERESA5aGhrayMWi7Fu3RoqKysL\n+rrV1dV0dKynsfEiEgHHTlyhttXe/WCw4gcWi3HDUbNxOSkrcL0oAKsC+4E/bDV+/PgRnWYtw6Mg\nRUREklRXVzN37twRna7rB0hPPfUUZ5/tv+7LuKJsdwEP43JO/AJurcA43JDPs8C3gW9SUVHlHbsC\nqMH1yNyHGxYq4/OfXxbIa0nUZolEFozE25Qh0uCciIiExs0338KePfv8R97PJbgcFXDDPRtxQYbF\n9bysAJZRUXEqhw8f8R6/G9gPfBN4DLgY2MiWLZtIV5ulvb2F7u5u1VEJGfWkiIhIKMRiMR59dCNw\nAomejkm4eijBRQS3AQtI1FR5NwCHD/fg6qYsAz6FP/XYreUzN/BKIzPNWoZPQYqIyBgWpgTSTZs2\n4QKKVbiejrdI5KakztRpwy0mCPBH72cZyVOOE1OPjbmVGTP8oOYRYC2JnJXCTLOW4VOQIiISYoUK\nIgq5Ts/w2+z3dAxWlfZ2YBxlZd/i7LNrcAFOcsVcf+rx+973Xh588AGqqibgelqacTkrkykrW1yw\nadYyPApSRERCqNCL/SUXRstPAmk2bR4ogDnjjDO8e35BucGq0q4ELMcfb7nyyo9729IHNCeccDw3\n3LCYQ4feIbmn5UUqK48t2DRrGSZr7Zi/AXWA7ezstCIiYdDU1GzLy0+10Gphn4VWW15+qm1qah72\nuV944QULeOe2gdtDFrCxWCzvbT5w4IBtamr2XtfdmpqabU9PT9/xbW1tFsosnOq1ZZ+FSRYqAo8f\n8p5vTmpz4rbCQpuFWIbn8/uexdrOzk7/862z+f5+zvcJR+NNQYqIhEmhggifCwbwvvSD599nAdvW\n1pb3Ns+YMXvQoCtxjkkpgcVRKY+bLfQktRlOsFCest8kC+Ms1Ae25e89i1PIIEXDPSIiIVPoxf4K\nsU7PYG3esmXToCssJ0rz78NNI34QOAs4AnzGO9/NwBrALzDnD/38JXASqUM5cCxuXR/fyKxNJPmh\nOikBS5cupaKigkgkQiQSKXZzRGSMSg4irgo8k58vVD8Y6OhYQjxucYHEJsrLb6ShIbcE0sHa7KQP\nYB577DF2797NxIkTiUZbiUQW0N6+LLBfGfBP3v3vAu24xQN/AyzCBSy7SK1/kqijchNVVROYOnVK\nXt/zWBeNRolGoxw+fLhwL5LvrpnReEPDPSISMon8jkQuRr5yUqy1tqenZ9AckXy1edy48RmGgu71\nclASbbjkkgZbX98Y2Fbm5aQkhoncY/+4cwYdyoH3WMBu37497+9ZCjvco54UEZEQSvQotPRta2ho\nztssFL8MfXd3d18vxnB7E9K1edy4CRw8+DaJFYsTvRjwOeBkXB2UWcBmHn3008BxuF6R03CVYv06\nKZDcQ/JB3NBPcFZQul6c7wFXsX///ry/ZyksBSkiIiFUiCAinerq6rydN7XN5eXlNDU14QKOZlyV\n2JaUo75FIrA4H3gb+LG3ba23PVOdlK+QCERm0z8IWgRMwNVPSQyT5fM9S2EpcVZEJMSKsdjfcPlt\njsfj3pZZuLyRNbhFAs8J7H0TcBlwkP7F2wark7IPuBEXAP0SmEZiFeUW4BjgC8AiZs6cPao+Q3EU\npIiISEGkn0X0deD39F+L5wrgO4H9e4CluK+phd5+L3k/F3nbl+ECEn+xwVQHgL8HXmPRohvy9K5k\nJClIERGRgkhMKV6CCyQ24dbcuQb4K5JL128EniWRu/JhXPByLzCF5B6SCvwhHDfNeA3wCWArycHP\nKbhem15WrfpB3qr1yshRkCIiIkMylLV5otFWGhr8YZh6b+sK3Lo5FwM/AwyJhQU3Ah/ALSx4F3At\n8GsghutRMbggxPeWd+6NJBYm9IOfVbipyeewdeuzwyr5L8WhIEVERLKSy3pClZWV3HXX93n/+z+E\nMRUk93TsBCK4hFpI5K58MfDYV40LOE5JOcfxJGb3ZEqw/Xq/wnEyOihIERGRrAx1UcJgUPPcc89g\nbbCn43zcsE8viWqyd+Nm9JR7j4O5LDFcb8k9JPeW3E2iZyVTgu1k8lWtV0aWpiCLiMigYrEY7e1t\npFZ1jcct7e0tdHd3J82eicViRCJXsXNnNy7BdQWup6MHNzzTFjj7A97P23ErG/fivp4WkphSfL+3\nT6beksn0n4Lsz/yp9tqt8vejjXpSRERkUNmuJxTsPenqepre3sO4acfgejpacAmxwSGbCly+Sisw\nDjf1+AjwGomE2ZWBcwT5vSXH038K8gTcjKFWystvpKlJ5e9HGwUpIiIyqGwXJUw3JORm7YwHrsP1\noCQvNOgeb8TN+LmTRL2U4FDQJlyviD9TKDgd+QzgCeAS3CKDZ3rH7AI+BLTQ0DAtb9V6ZeRouEdE\nRAaVzaKEmYaEEmXsfZmGbHYH7oP7f/TPvPt+UJJatbYM8Be4u9l7rWXArdx2223U1dWp/P0opp4U\nERHJSvJ0YjekEuyhGGxIKCHTkM3EwP1jcD0pU0kUdFsD/BCX33ICcDRuNWS/x+Zk3OyffwTKuOKK\nK0ZdtV5Jpp4UERHJymDrCSUPCaVb6A/6J8Ruwg3h1ANPkqgm+2dv/6uB7wI30L8H5V7S99gcTX39\nhxWclAAFKSIiMiQDL9BXBiymfxBSxkknnczrr1tczkgw4CjH5aRsxH0t3QvU4gKUa3EBykpcAu6t\nuAJuvcDclNd2PTann/4X/OIXjwzvTUooKEgREZG8cMM9vbjpwMEgZBqwjddfP0wiX6Ubl4PyHC6H\n5H7gkHf/xyRmBAHcgQtSynABjb9OT/oem89+9jNUVlbm8Z1JsSgnRURE8iIx3HMFyXko20h83fj5\nKtW4npDzA/s14YKQF0ieHTQOt6bPyUAcOM073yKSZ/q4Hpsrrrgin29LikhBioiI5MX48eOpqpqA\nK6L2DMmBxsm4rxw/afYpXIBysff4GuACXBCSWlX2TlwJ/a/gemqacL0p/ro9fl2UN0m/GrKMVgpS\nREQkL+bPb6Gn503gbfov9nc3LsC4HjccdAHwtHfku7yfb3k/M80Oerf388e4YGRyyn4XAVal70uI\nghQREcmZvyLy+vXraW9vw9rrvGcyBRpvAC+S6GWZhJvJ0wo85u2zCpez4vNnB/3R++l/dW3zzvsI\nbm2fvwNU+r6UKHFWRESGrKenh49+9GNs2RKseWJwCwRC5mnIvSSGc2K4YZxWXH6Kn2y7wrvV4/Jb\nbsEFM7cBxwLH4eqofB83k+g+XICyiLq68zX1uIQoSBERKQGxWIw9e/aMSHXVnp4eamrO5cCBVwJb\nJwHP45JYJ9F/sb9FwFnAXhK9LH7gcjrJa/rMwgU5C0n0rvQA5+BK3d+IW4zw3biS+i246ctl3Hff\nD/L7ZqWoFKSIiIxiPT09zJ/f4pWjd5qamolGW4c9DTcY+Fhr++5fffU1HDjwJ5IDiiW4RQHvwa2x\nk658/UPAPNz6PatJrITsDwVlKqd/Ai7J9iFckuwRb5+t+EM8xpzApZdezNSpU4f1niVcFKSIiIxi\nyQv6uYCho2MJkcgC1q1bk9M5+wc+ZbhhmqCB1ufZhitN/4+4QOJT3vEveOdaggs8/DYvxyXDZspj\nORf4Ja4sPiS+uqbj98ZcdNH5WkCwBClIEREZpTIt6BePW9rbW1i/fj3xeHzIQ0DJgc8DwA5cMuss\n7+cK+gcU5+ECkAWBbc24nhNwBd2+igtm/uyd189D8YOhTHks38MFKItwNVPuw80IeoWysiVcdNFs\nNm9+LOv3J6OHsVZzyo0xdUDnrFmzqKioIBKJEIlEit0sEZEBrV27lubmZtwsmdMDzzyLywtJ9H5k\nOwQUi8Wora3FBSjn48rTB4OgWJpt4KYDv4gb7vGHgBbjapf8GReclOGSXt/02nwdrtflLhLB0F0k\n57G85p2/F1foLR74mb+hLRm6aDRKNBrl8OHDbN68GWCKtbYrry9irR3zN6AOsJ2dnVZEZLR44YUX\nLGCh1YIN3CZZqPC277PQasvLT7VNTc2DnrOtrc075z4L/v3HvPsx7/z13vkf8va7PUM7HvK2B28r\nUn76x/RYaE7Zt97Cfd5rneRtm2cBe//999tYLDYCn7IMprOz079mdTbP38+qkyIiMkrV1NTQ1NRM\nefkSEuXhV+Cm9SZXbY3H76S9vY3u7u7MJyR1JeMqXO/HxbihmxrgMu/2Golqr1/wjsmUU7IMeNC7\nf6V3ruUpx1TihnT8IZ77gV/jFhi8G3jd2+5qpdx7732MHz9+wPcio5+CFBGRUcYvoNbd3U002kpD\nwzSyDRgGq8bqBz5lZYtw5edPJrm8/VbgFurq6rwj7iGx/s7mlLP5Acc1uJwUf59WXAd2umP2JbU3\n+f5xuOEh2Lmzm0gkmP8ipUiJsyIio8RA041fffVVdu/eTXl5OU1NTWRKQs2mGms02spZZ1Vz6NAB\nMs3iueWWL/Lxj38c1+PhV45NVxulHreYILgeFH+fVlwQtDDNMZMCxyTaDn/yfk6lt/dG2ttb6O7u\nVvG2EqYgRURklBhsurH/Zd3U1ExHxxLi8cSXf3n5jTQ0NGf1hb5//34vQIFMPTInnHACM2bM8irO\ntpK5NsqVgcetuKAluM+xKY/LcQXfWkkELjd6x2309rkXmAC4niEFKaVLwz0iIqOAP904Hr+LwXJN\n+g8BtdDQMC3rOiJ79uwJPEo/hHPUUUexePFCb9ssEjklMRL5JwBfIpEvswbX62KAv/Gefzuwb5l3\ne53k1Y2n4crjgwtcpjKUniEZvdSTIiIyCiQCh8y5Jn6PQmVlJevWraG7u5vdu3cPuU5KInn2g7jh\nl9ThmDIuvfRSZszwc0WCQ0vVwJPe/V4SgYZvEq7E/QXAvwa2H42rJDsR2I0r9nYd8BFcYOPXSFkF\ntA6pZ0hGLwUpIiKjQPKsm+xyTaqrq3P6EveTZ9vbO3C9HqlBxk+BZ9i6dQlVVRM4dCh5aMkPZFxw\nsxs34+jdwO+Af/DOc7P3069me8Q79w5vexxY6d1g3LjxHDr0KvAhABoamlVhdgzQcI+IyCiQfrqx\n61Foahpej0JwtpDvm9/8Oq4I24+Bdm/rSlwQ8UHgfOLxqzlw4BVOPtkPZPxek7dwgYfFDd3cjCuN\nvxw4nuTZQifjvoqOwRWD85/7KcacQl3dVGKxGAcP7icWi9HW1kYsFmPdujUq4DYGqCdFRGSUiEZb\niUQW0N6e6NkYTo/CYLOFnFnAb737V+CGaoKl7OHQoYPAF4Ff4VZCPh5XQTa4mvE7uEqzd5N+zZ+j\ngfclPWetpasr8V5z7RmS0Us9KSIio4Sfa5KvHoXk2UKuZ6OjYxuRyIKU4aXg/RZcrZLU3pB/xAUz\nvaQWknOByZveOTIVfLsOV4RuQ7/nBqvtIqVLQYqIyChTXV3N3Llzhz3EM9BsIWNMYHjpKdwU4Btw\nPSjJx7gg5E/Ap72zZwpEIHPBt494Py/FVbQ9iGbwiIIUEZESly7nJJvZQslTmTeSWOwvUxBygfcz\nUyBicEm1iZwaVwPFXyQR4BxcT019XvJtZHRTkCIiUqJ6enqYM+cyamtraW5upqamhjlzLuPgwYMp\nwzlBid6LyspK7rrr+3z2s5/lXe96Ny5/JPMxiTV+kpN7E8XYLDCe5CTbQySmGJ8D7MLVVtnJ9Okf\n1AyeMU6JsyIiJWqwCrUDVaatqqpi1qxLePzxzbg8E6es7BhgMb29mUrZt9K/8mwj8Hvv/l7gJOB6\n4K9xgcxC3CygXbj/O58IwC233KwZPGOcelJEREpQNhVqB6pMO39+C48/vhWoIJgk29t7HL29/01y\nb8jrJErZvw5EcIHIWbgg5iXv+S/gAp4fArfjAic/p8Wvi3IyrmCbclFEPSkiIiUp2wq1qZVprbX8\n4he/CExL/jHppwwbEsM/ce9nsPdkAi4wuYRET8ztuP8bfyhtm+BcXMDSwrnnflC5KJJbT4oxJuP6\n2MaYFbk3R0RE8iGbnBNfdXU1F1xwAYsXf47a2lquvfbawP6ZkmRTi7L5gciDuPV7NgMfwA3dpE5X\n/tu0bXL/b3bnv/baTyOS63DP3caYy1I3GmO+jxuMFBGRIhpqhdr+NVO+4D2TKUn2GyRPQ/6at/0t\n4HNALa5uymvAw7jhH39oZyduaCeYWFuGK43vzt/c3Dy8D0BKQq5ByieBfzbGzPQ3GGNW4coRXpKP\nhomIyPBkuxpycv7KXFxhNX9oZiGJIGcF8FncUM9UYC3gT2v+pLf/jcBWkntPtpH4/6vfE7OMRE7L\nGbiemOeARdTXN2qoR4Acc1KsteuMMdcB/2aMuRS4GvgocIm1NpbPBoqISG78CrXr169n27ZtTJ8+\nncbGxn77JeevBCvKfgg3NNNCYiFAX33gcTMwz3v8NpnzWLpJrJB8NvB5XCXa27zzLwPK+O53/UUI\nZazLOXF1oH9nAAAgAElEQVTWWvsvxphKYAuwH5htrVXtYhGRIonFYuzZs4eJEydSXV094No8wam9\nifyVR3AVZVtJBBk7cAmt/4Urd++vx7MYmIz7P+pi4Ne4tXd+R+Y8lh8B9+MCkj24Xhq88/9f3MKF\n29i/f3/uH4KUlKyDFGPMHRme+iPuX/ENxhgArLWfH37TREQkVWogApkXCnznnXfYtKmTTHVSkpXh\nVimG5CAjhqtfEgxcgr0jP8SVyW/BBSjggpirAufw81hWAifgelxW4oKTiSTqq7jZQZp6LL6h9KRM\nzrB9D3BK4HmbYT8REcnRQL0i6Yq2bdiwkN7ew6QGF/G4pb29he7u7r7elkjkKlzg8AHcUE8wyBh4\nKjPsJjnP5N+Ba3HF267EBSg3AucBv8Gty/NzXArj6f3OV1c3Vfko0ifrIMVaq4RYEZEiyVQ9dt68\nj7FlyyZSg5HeXn8GTfrg4pOfnE9Hx3rmz29h504/+fV/4YKOhbj/b87GzdABV2DtGlyvByR6R8px\nPSng1vf5P979ZcDNuOCnEfj/cAm5c3FBSvrelvvuuzfrz0RKn7FWHR/GmDqgs7Ozk7q6umI3R0Qk\nSSwWo7a2luRABIJDJG4mTbBnYhNuLZ30x5SVVXDhheexZctmb9vDuLySY4A3SCTFpibMngf8DXCH\nt+/+wH4nk5y3sgiXGFvrtc/gVje+DNdjcyfJ5finpRmGkrDr6upiypQpAFOstV35PHeuxdxONsas\nMMZ0GmP2GGP2Bm/5bKCIyFg3WPVYJ7WeyUtAGcbcgOvV8IMRt+pwb+/dXoDin/fruJk51+GCkk3A\n+cA4kqcT/w6Xu/I68Gdv22PeMfeQXDtllbfPs97+vd7+3yF5+nH6qdEiuc7ueQg3Sf5BXOKsiIgU\nSHL12P5DJDNmzGbr1uSFAsvKFjNu3Hh6ev6IG/ZZift/6YdJrLFD4Lynevcvw9VDeQrYTuaEWUgE\nJWu9xwMFUXFc6mKidP7UqeezbNlNTJ48WXkoklauQcqluCnH2/PZGBERcYKzePzqsZlWLI5GW4lE\nFtDenggAKisncOjQOwRzWGAJcDRuiOZ+AGbOnM0TTywhHv+id+R/evt/1Xs8UODhPzdwEJVwNfDX\nlJUt58ILp/D4449l/4HImJRrkPIKbmBRRETyKNMsnnvvvZvrr1+UFIj4AYpftM1fKLC8vJympiYy\n94Ikzv3OO0eorf0fPP/8zbielk/jhn18AwUe/nM1uIJui0kk3G7CBUXHAu/ghno+BUBjY7OGdiQr\nuQYpXwWWG2P+zlr753w2SERkLMs0i+f66xf1W7E4dYikurqa6upq1q4dbPhlGS6geIRt227Frbdj\ncYmtx+Eqxs7CVZENzvTxpxM344q7BZ+bh0u8Da6EXIb7mjkKOMLy5bcSiUQ0tCNZyzVI+Xtc/95+\nY8yLJKd+Y63VFBkRkSHy19AZrLbJYF/yg+WwwCdwCbJtwaNwNVHuCRyzEVf+Phh41OMCkptxs4CC\nz52C+zp4x3vciwt8jgX+WwGKDFmuQcr/zmsrRERk0Fk8u3fv7vcln64C7fjx4zn66ON5553UXpBF\nuNk6y3AzboL5KotwPR/B164EVuNm4EzAjfRv9G7+1ORv4FZA/iauk70V+CtcvZWJuLV6WpgxY7YC\nFBmyXBcYXD74XiIiMhSD9YAEy8UPVIH2ox/9GO+88xauPknq8Is/vThTvsrPcJ3lya8NrwHn4Erk\nGxId6F8jUSMFXJBzOomib8cBsHjxDYN/ACIpcqqTIiIi+efP4ikvX4ILIl4CWikvv5Gmpuaknojk\n3BVXw2TDhq3MnHlxoP7JatzaO1OBCuCnuMoRkDlf5WtJr52YEfQnb9skXA9LsHbKycBPvONT67W4\nIGfy5Ewrq4hkllNPijHmKODLwMeB9+DqIvex1p6a7jgRERlYuunE/iweX6bcld5ey3PPBXtONuMK\nsj0d2DcWeC5dvspbJPe+TCCRY/I6sLPf6yZ6YY4GbiA4xORPk9ZQj+Qi15yU7+JWh/oe8G3coORp\n3rZv5qdpIiJjT+p04nSzeLKrQDsJ1wtydcq+/nThJfTPV/E71//du/9dXFCyDPhr4Fe4Qm+n4Qq4\n+SsY+697NMa8hbWZAyyRocg1SJkPfMZau8YY803gEWvtXmPMM7h/rXcNfLiIiAxkoFk8g8/emY1b\ncfh9uEqzqfu2kn7WzpXAF4CvAPcBj+NyT1Z4t8m44OXiwHHNuNk+AMuxdhk/+clPmDBhQtoAS2Qo\ncs1JqcL9BoBbmMHPmFqLq6ksIiIFVFd3PmVli0nOH/FrmPwSuBDXCwLuT/3CwL5rcLNvwPWSxHA1\nTq4F7vaO+zhwIq5Hxvcb3J/7YD7K1sDrXgnAvn37mDt3rgIUGbZcg5R9wLu9+924EBxgOik1U4rB\nGHOSMeZpY0yXMWaXMebzxW6TiMhw9fT0MGfOZdTW1tLVtZ3e3sMEF+lzCa3LvZ9rcL0f8OCDP2Hc\nuKNT9vU70heTmIkDiaGbl3C1U/xg5DHSLyJ4F65C7XL8npzp06fn943LmJVrkPK/SQQm9wArjDG7\nSE4dL6Y3gIu8onKTgUXGmPcWuU0iUuJisRhr166lu7u7IOf/xCc+yfr1T5DoyfgpcBKJP+V7cImy\nHwbuo6zsNs499wP8+c9/ZsWKb3v7fB5Yj6tfAplm4zg7cUHIVcCb3rZMeTCPAIuoqppAY2Njrm9R\nJEmudVJuCtx/yBjzEq4XZQ/w8zy1LWfWWkti8YkTcENSr2c+QkQkdwPVLKmsrMzLa8RiMTZu3ED/\nmTX/hcsjuRk3pLPZ+/kYvb3w/PO/5ZprrvH2L8PVTmnEDfFU0r/svZ9A63eKZ7uI4EqqqiawffvW\nvLxfEcihJ8UYc5Qx5ifGmPf526y1j1lrv22t/ZkXIBSdMabCGLMT99+Nu6y1rxW7TSJSmtLVLOno\n2EYksiBvr7Fpkx8M+EFDDy4F8Ave4+/iSt034/JKwNVGSa1nshjXwVyLWyf2TJKHgV4jUc4eEj0t\nwUUEg3kwizjzzImsX7+eV1/9A2eeeWbe3rPIkIMUa+0R4GPAkXw1whgz0xiz2hjzsjGm1xgzL80+\nNxhj9hpj3jLGbDfGzBiknYettZNwv4ELjTFnD7S/iEgu/Jol8bg/LOJyNeLxO2lvbyvA0I8fNLQA\nyYGRe7zAa0MvsIrk/JG7cTVPYt5+kL7g2z7gh7jO9mDC7TzcsE8iqJk5cxKdnU9piEcKIteclF+S\nmHOWDyfiBj/9KkBJjDFX4mqy3IwL/9cDa40xpwX2WWyM2eEly/YNY1lr9+Myvqbmsb0iIkB26+1k\nY6B8lp6eHn7601bcn+zFuITYNhL5In4Q8kVvu9+TEmxTDy6F0OICDf98a71tT5OcFOvP9HmNRFBy\nHXARLkn2BKZNu5DNmx/L25CWSKpcg5TdwG3GmIeNMTcZY5YEb0M9mbV2nbX2a9bafyPRxxi0FLjH\nWvtza+0+a+2XvTZcHzjHKmvtZC9ZttIYcwKAMeYU3F+L53N4nyIiA0quWRLUf72ddIIzdpqbm6mp\nqWHOnMs4ePBg3z7z57ewdetvgXNxAYY/xJNp6OcXJAIQXwvwf+g//LMEN1QELjgJasb1yMwDzvO2\nbQRupanpYtrafjXgexMZrlyLuX0GOIxLlk2da2bJYzE3Y8zRwBRcVdugX+MKAaRzBvBjYwy4Os13\nW2ufHey1li5dSkVFRdK2SCRCJBIZarNFZIzw19vp6FhCPJ5cDn769Fl9PSmZaoYk57O4FYk7OpYQ\niSxg3bo1KSXwm3HDNH6Crp/EGhz68Vc1XogLQE7ABR/9y+gnytk/4G2b7b1GKy6p1s+DWd3X3rq6\nqdx3371MnarO6bEoGo0SjUaTth0+fLhwL2itDdUNL2wPPH6vt60uZb+bgF15es06wHZ2dloRkaHq\n6emxTU3NFvetbwFbVTUh6XFTU7Pt6elJOu6FF17wnm+1YAO3hyxgt2/fbtva2rx99lk4YMF/nTIL\nFRZuH/Ac7mYC5wjus8/bvsy732qh0kK9d3yF9zrYurqpdvv27UX6hCXMOjs7/X9rSd/T+bgNaxVk\nbwbNucEckBEUillEIiL+ejuxWIy2tjZmzpzNoUPvMNhsn8HyWT772etThpOCPSY7cZ3GqUM/yec4\n9tjjSXSaZ6qJcg3Jxdk24ifFPvJIlFgsRmfndvWeyIjLdRXkk4B/Aj7hbaoG9hpjvgcctNbelqf2\nAbwKxHGrLQe9F/hDHl9HRGRYqqursdby+OObSB1aicct7e0tdHd39w39DLYGT1fX0xhjaGpqZsOG\nhV6F2eB5Twf24upXpj/H22+/BZyCWwxwMf1rotSTruLstddey0033aTS9lJUufakfB/4S1xOyFuB\n7R24lZDzxlr7DtBJosKt7xLgiXy+lojIcA1ltk9NTQ3nnvtBkqf5+mvw1PftH422MmmSHyz4543h\n8ky+B4xPOccPcamDvv/GJc2mltF/DX+9nQQX3PzoRz9Km8QrMpJyDVI+Biyy1j5J8rDLcySH5Fkx\nxpxojDnPGOOvZHWW99hPNb8DuN4Y8wljzF8aY27zXueHObZfRKQgsp3t48/qef75Z0me5tsCTMP/\n/97EiROprKwkGv3nlPP6wdBXcUW1g0XZFgLHkTyTZxxu9eKVuGRafyLlF0gtzuYWFSxMUTqRocg1\nSDkB+GOG7X/O4XxTgR24HhOL+69BF24yPtbanwF/D9yO++9AEzDXWvtSDq8lIjJsmeqa+LN9ysuX\nEPzyLy+/kaam5r7hk+RZPb/DBQ8n4QKUCGVlN1NXN5UXX3yRtWvX9g37JM7rBxr7gR/g/oTGgB+R\neSHAjbjpxPfh/tR+HfgTyQHSmd5+hS5KJ5KFXLJtcf8luM67/xpwpnd/JbAu39m9hb7hze6ZNWuW\n/chHPmIffvjhIWU2i8jYceDAgX4zeVJn7qSb7RPcJ3lWzwsW2izEAjNyypKO9R/X1zfa+vrGlOdS\nZ+20pdkWnMnTlnLff80lAx7X1tZWrI9cQurhhx+2H/nIR+ysWbMKNrsn1y/1C4FDwJ24nJTv4ibS\nvw5MyXcjC31DU5BFJEtNTc22vPxUL8Bw03bLy0+1TU3N/faNxWK2ra3Ntre327a2NhuLxay1NjCt\nuD4l2Ki3YGxZWWXS+eFUC5P6XicWi9mbb745MLU4OP144GnNycFQLBCw/GjA4/y2i6QK3RRka+0T\nwEzcIOce4CO4zKzp1trOXM4pIhJ2Q12np6qqijvvvJumpqakarLjx4/HjbbvIDlvZAdg6O29BTgV\nNxRzFe7/gzuJx7/Ut9Lypk1bcLN2ziY5afYp4Fhcbkkw12QJLhn3SVxibjMute9fvNZmOm4RM2fO\n1iwfKYqc65tYV8H1U3lsi4hIqGUzcyf4ZZ6pmuwbb9xE8gKAkFwBdlng3M3Ad7z77wbgscceY9u2\n//DOW4PLY2kJHNOIq9AQ3FaGyzXZiEuM/RJuNeSd3vP/hPtK+B9Jx1VVTeDf/u2XGT8TkULKqSfF\nGHNKhtvJxphj8t1IEZEwGMo6PQP1umzZ4h+fPthxQUpwZWO/ssNzAPzxj/68hVnA+cAc3ArGN3lt\n+VvveDeDZ8aM2axadad3zDm4wGQ28CLJPTkn4eYm0Hdcd/cuLSAoRZPr7J5DwME0t0PAW8aY/2uM\nWW6MKc9PM0VEii/bmTsweK+LM1gF2Lm4WTd+4HA7UM6bb76ZcnwrbnXild75/dk6vwFa2br1WX71\nq7Ve218BziL9DKBVQC/3338/sViMxx/XCsdSXLkGKZ8CXsbNX2v2bl8H/hP327UKV9own5VnRUSK\nLhptpaHBH15x03YbGqYRjbYm7TdYr8vMmbP7BTsuH6QMl/L3ovcafi9HorejtTXq7bfY2/afwH+l\nvM7pwGkEc2Zuu205kyadjatSC5kCqL/4i79QDoqEQq45KS3AQmvtvwe2tRtjdgKfs9Z+2BizD1eE\n7UvDbaSISFj46/R0d3eze/duysvLicfjvPrqq2l6HfxAIliKfglQxh13rOQrX7mV9vZg3kg9rgLs\nF4C/wq0K0n/l4n37/GP8InBlwMkkr4K8BLdi8hqCawF1dT0deL30pfSDw1YiRZXLlCDgTeCMNNvP\nAN707v8F8Od8T0cqxA3VSRGRIRqsXsrA04xd3ZF169Z525Z504FTpwsPvHKxMeMC+w0+5bisrMLC\nCm/7JG9q80PeOd2qxzNnzi7uByujRpjrpMSA5Wm2Lwde8O5PAv6Q7wYX4obqpIjIEA1WLyW5YFvM\nphZsi8Vidvny5YMEItnWOxnoHDd57SyzieJxUy2c7AUqiQCqqmpCUlE6kWyErk4KLoX8i8aY7caY\nH3i3p4CbceXrwaWcP5Lj+UVEQiubeinJSbZPAh8AnkxKsr3gggu8M6bPW/nABz5IWZmfdxJcfNCv\ncZJNEu5KTjzR4hJlHwBqgadxQ0XP9O2tmTwSRjnlpFhrVxtjJgOfBs71Nj8O/C9r7fPePvfnp4ki\nIuEy2Mydxx57jN27d3PbbcuB5LyThobmviTbpqYmxo0bz6FDC0nOW1nEUUcdR2VlFb29z5Fc76QZ\nF6xAIhCZRvrcl/OAFl5//ZskF4/z81YWA4dZv34djY2Nw/hERApjOMXcnifRayIiMmYkz9wJJp6u\nAcq49tpr+7Y0NTWzfft29u/fz8SJE/tmzfT09DB/fguHDr0KlJMciIzjyJE4W7bsBH6K66n5O9y6\nrhHcCiRrcIGIH7QsIH0wU0lvby8uGTd98bj3ve99w/g0RAon5yAFwBhTBbwH9xvWx1r7TPojRERG\nP38op6NjCfF4sPficxhzCtbeTbDCLNzKunVrks6RqEa7Ale87WbcysbTgb8Ear3z+EHF07jZP6mz\ngVwg4oKWzV5bvosLSnwTvJ/ZVcoVCYtcK86+1xjzKC6sfwbXhxi8iYiUtHT1UuDtQGCReV2f5JyW\n93tbFwJfw5W0TzecVIlbxxUWL17sbbva2+7b5/1MraP5ivdz8Eq5ImGSa+LsPd6xU3GrIDcCn8SV\nRZyXn6aJiISXXy8lFovR1tbG/ff7aXiZeyt8yTkt6Yq+DVwIbvHixRkr31ZVTaC8/B9Stn/H2z54\npVyRMMk1SPkwsMxauwOXMr7PWvtz4POoeJuIjFKxWIy1a9f2W814INXV1cydO5dZs/zgJHNvhX/+\n8vLywL41uPyRYADhViQ2JnlF4mBQkany7fbtW4e0PbVSrkiYGOvqhAztIGPeAiZaa182xvwBuNxa\nu80Y85fALmvtCfluaCEZY+qAzlmzZlFRUUEkEiESiRS7WSIyQvwk1vb2tr5tM2bMYvHihUyePDnr\nnoaZMy/miSeeobf3Lvw8lfLyG5k1q45jjjkm6fxVVRM4ePBtentX4Wbh/C2JFYmhvt7Nttm4cUOg\nTbNZvfqXSdOE/cq3waTcXLaLDFU0GiUajXL48GE2b94MMMVa25XXF8mluArQBdR799fi1viegMvW\neinfxVwKfUPF3ETGtOTCbM/0K3IWrCSbTnL12bJ+x9bXN6Yt/FZVNSFp3xkzZttHHnnExmKxvvPO\nnDl7SG0RGWmFLOaW6+yefwT+h3f/VqAdl8EVxy0wKCIyKvhJrIk1ci4jsZhfYoZOJLKg3wwdX2Km\nzg+BnwEb+557443X2bJlM6lr8MTjlgMHWli/fj1HjhxJ27Mxf34LTzzx7JDaIlJKci3m9tPAwx24\n+XLvB35nrf19PhomIjIUsViMPXv2DHkYIzmJNQYEAxbwA4r29ha6u7v7nTs5yPn/gf9Ien7LlicC\n5w9yCbVHjhxh7ty5ad9PcvA0eFtESk2uibMYYxYbY14A/gQcBB4ErshXw0REMgkmuPb09DBnzmXU\n1tbS3NxMTU0Nc+ZcxsGDB7M6V3JhtoEryQZn6Pg2bdoUeLQRV+uklURvjJ+iN7Tpv4NVtU3XFpFS\nk2udlG8C3wL+Bdc3ehkQBb5hjLktf80TEUlIF5DU1JzrDbUkAoOOjm1EIguyOmfyGju/9bYOHlD4\nbUlUl/Vf7zsE66T4FRtS1+AZbPpvcvA0cFtESlYuiSzAq8BH02y/HHg134kzhb6hxFmRUaH/ysMr\nBlwp2E9AHUxPT09K4muFd459Fh5KWt04c1tavePq065G/P73f2jICbCJ1xi4LSLFFMbE2XKCc+US\ndjDMUvsiIumkz9Hwq7UOr9y7X5itu7ubHTt2cPfdP+Dxx9MvCpi5LYm1cKAbt0ox+D0fv/zlL/ra\nlG3eTDTaSiSyIOMChSKlLteAohW4jv6F264B/nlYLRIRSSN9jkamhf5yGxKprq6murqaK664YsB6\nIoPli8CPcMXZNmHMYi65pLHvHNbavnySwQKVYPCk2iYyFg2n1+MGY8wcYBvuvw/TcH8xfmqMucPf\nyVr7+eE1UUQkNUfjfFyS60RgEm7dm8RCf+XlN9LQMLxy737AMnhb+gdHsNK7waWXup6PdAXjmprc\nc8HibENti0gpy7Xi7KNZ7mqttfVDfoERpoqzIqNDfX0jjz76OPB2YOuxVFaewsGD+/u2ZPvlP5DB\npjR/+MOX8uij27F2FX5w5HpNzueHP7ynX8/HnDmX0dGxzVtU0NU8KS9fQkPDNNU8kVEptBVnS+2G\nEmdFRoX6+kZrzLikZFVjxtn6+kYbi8VsW1tb1smymSRXj82c5Fpf32jh2KT94FhbX9/Y75wvvPBC\nXhJ8RcKokImzOddJEREZSbFYjI0bN2Dt3QSn+Fq7qm99m7lz5w5rWCQWi9HYOMeb0rwCV/5pZb8p\nzX5b4MckCsDFgH9i48YN/RYoVM0TkdwoSBGRUaGQX/TB+itdXduJx88AlgGfAm4iHj+D9va2vuAj\nuS3VwFzvZ/q2qOaJSG4UpIjIqFDIL/rE2jvLcH8W/WqxwZ9lfcHHUNuSXDAu+4JuImOdapqIjBG5\nrm0TFv4XfUfHEuLx/M3kSa55Mh43zHMX6eqfHHXUUX1tmTJlKl1dn8HaZ4BFA7YlFotx9dWf4o03\n3mDLFtU8EcmWghSREjecqa9hU4jiZslDN78N3A9KLAa4Z88eLrjgIg4ceMV77nbge0C8X1vSffYz\nZ85m0aIbmDx58qgMFkVGkoZ7REpcYigjt7VtwsQvbhaLxWhrayMWi7Fu3ZphBVtlZf6fwc0kF4cL\nSgzjuADlTyQPB53EuHHj+7Ul3Wf/xBPP8sADDypAEcmCelJESlim8u3xuKW9vYXu7u5R+WWZj+Jm\nyb0cZbiCcHcD9cBigsXhYBFVVRPYu3ev14PSvxz+oUMtbNiwgcbGRqB0P3uRkaSeFJESpqmvmSX3\ncuwEzsStu7MROOzdP8P7eSYHDrxCW5s/bJP+89y6dWvfFn32IsOnIEWkhGnqa3p+L4er/no+8J/A\nz3BJswAPkVz/ZDUA73rXu7zn03+e06dP79uiz15k+BSkiJSw0Tb1tb29nW984xts2LChoK+T6OV4\nAKgFmoEaYK23fSfJ9U9cYHHllVdSVTUBNzSU+Dz94SB/qAdG32cvEkr5LmE7Gm+oLL6UsJ6enqzK\nvBfT7t27bVXVhKQ2VlVNsHv37i3I67ky9WUWKpNK7LvHZbasrMIrWb/PwkO2vPxU29TUbK21du/e\nvVm3dTR89iLDVciy+DktMFhqtMCgjAXd3d39Fr0Li/Hj3+PNmLkHf/E9WEhV1XG8+uof8v56sViM\n2tpakpNa8R63MGPGbLZs2dS3Nd2U7Q0bNrB161amT5+e1IOSTpg/e5FcjcQCgwpSSAQpnZ2d1NXV\nFbs5ImNKe3s7c+bMIVPAsH79+kGDgKFau3Ytzc3NuGnBpweeeQk4g7a2NiZOnKjAQiQLXV1dTJky\nBQoQpCgnRUSK6sknn/TuDT5jJl+ySWqtrq4e9oKFIjI8ClJEpKguuOAC797gM2byRUmtIqODghQR\nKZhYLMbatWv7Vg9Op6mpKesZM/kUjbbS0DCNYD2UhoZpWktHJERUcVZKTpgW0hurbRnqekHbt2/l\n/POnc+BAYk2eqqoJbN+e/6Een19iX0mtIiGW7+lCo/GGpiCXhAMHDoRmuudYb0tTU7MtLz81aXpv\ncBpvJuvXr7fLly+369evL1jbRCS/CjkFuegBQhhuClJKQ65fjGpLfrkaJHivZwO3hyxgY7FYQV5X\nRIqjkEGKclKkJCSXOb8KN630KuLxO2lvbxswJ0JtyS+tWSMi+aIgRUpCmL4Yx3pbtGaNiOSLghQp\nCWH6YhzrbampqaG+vhFjFhGcrWPMYurrG5WcKiJZ0+weKQl+3YuOjiXE4xbXU7CJ8vIbaWgY2boX\nyW35L2AC8EfKy79d5LaM7Odi7Vu46b3+42ML9loiUprUkyIlI0x1L37wg1WMG3c08AXgU8Ayxo07\nmnvvvXvE2zLSn0ssFmPjxg3Aj4EY0Ob9/Cc2btwwojk5IjK6qSdFSkaY6l7ccMNiDh16Bzfc4RbM\nO3RoCddfv4h169aMaFtG+nNJzoM5HfBf6zjA5cFoyEdEsqEgRUpOdXV1Ub8E/Rk1yQvmXUU8bmlv\nb6G7u7so7RupzyU5Dya4YKASZ0VkaBSkBCxdupSKigoikQiRSKTYzZEcFbvKazYzakq5JyFM+UEi\nUjjRaJRoNMrhw4cL9hrGumJmY5oxpg7o7OzspK6urtjNkRwNtRR7ocRiMWpra0nuScF73EIsFiv5\nL+qDBw8SiSwo+rUQkcLr6upiypQpAFOstV35PLcSZ6VkzJ/fQkfHNlwwsA9opaNjG5HIghFth1bY\nTeTBxGIx2traiMVirFu3RgGKiAyJhnukJIQtDyQabfV6EhJTcBsamsfcCrvFzg8SkdFNQYqUhLDl\ngYRpppGIyGilIEVKQlhnlKgnQUQkd8pJkZKgPBARkdKjIEVKRpgqzoqIyPBpuEdKhvJARERKi4IU\nKW0wELQAAAx/SURBVDnKAxERKQ0a7hEREZFQUpAiIiIioaQgRUREREJJQYqIiIiEkoIUERERCSUF\nKSIiIhJKClJEREQklBSkiIiISCgpSBEREZFQUsXZgKVLl1JRUUEkEiESiRS7OSIiIqEVjUaJRqMc\nPny4YK9hrLUFO/loYYypAzo7Ozupq6srdnNERERGja6uLqZMmQIwxVrblc9za7hHREREQklBioiI\niISSghQREREJJQUpIiIiEkoKUkRERCSUFKSIiIhIKClIERERkVBSkCIiIiKhpCBFREREQklBioiI\niISSghQREREJJQUpIiIiEkoKUkRERCSUFKSIiIhIKClIERERkVBSkCIiIiKhpCBFREREQklBioiI\niISSghQREREJJQUpIiIiEkoKUkRERCSUFKSIiIhIKClIERERkVBSkCIiIiKhpCBFREREQumoYjcg\nTJYuXUpFRQWRSIRIJFLs5oiIiIRWNBolGo1y+PDhgr2GsdYW7OSjhTGmDujs7Oykrq6u2M0REREZ\nNbq6upgyZQrAFGttVz7PreEeERERCSUFKSIiIhJKClJEREQklBSkiIiISCgpSBEREZFQUpAiIiIi\noaQgRUREREJJQYqIiIiEkoIUERERCSUFKSIiIhJKClJEREQklBSkiIiISCgpSBEREZFQUpAiIiIi\noaQgRUREREJJQYqIiIiEkoIUERERCSUFKSIiIhJKClJEREQklBSkiIiISCgpSBEREZFQUpAiIiIi\noaQgRUREREJJQYqIiIiEkoIUERERCSUFKSIiIhJKClJEREQklBSkiIiISCgpSBEREZFQUpAiIiIi\noaQgRUREREJJQYqIiIiEkoIUERERCSUFKSIiIhJKClJEREQklBSkiIiISCgpSBEREZFQUpAiIiIi\noVTSQYox5nhjzO+MMbcXuy0iIiIyNCUdpABfBrYWuxEy8qLRaLGbIHmma1padD0lGyUbpBhjJgK1\nwNpit0VGnv4Alh5d09Ki6ynZKNkgBVgJfAkwxW6IiIiIDF0oghRjzExjzGpjzMvGmF5jzLw0+9xg\njNlrjHnLGLPdGDNjgPPNA16w1u72NxWq7YWQ7/9hDOd8Qzk2m30H22eg5zM9F/b/kRWifbmec6jH\nFfKajtbrCeH5HdX1zI+wXM+hHjsW/uaGIkgBTgR2AjcANvVJY8yVwPeAm3FDOOuBtcaY0wL7LDbG\n7DDGdAEXAJ80xuzF9ah8xhjzlcK/jfzQL8zQngv7H0EFKfpSK9T5dD3zIyzXc6jHjoW/uUeN2CsN\nwFq7DlgHYIxJ1+uxFLjHWvtz7/GXjTHNwPW45FistauAVYFjvuyd71PA+621tw3QhOMAdu3aNZy3\nkTeHDx+mq6srFOcbyrHZ7DvYPgM9n+m5dNvz/RkORyHakus5h3pcIa9pvrYXQ1h+R3U98yMs13Oo\nx4blb27gu/O4QRs9RMbafh0XRWWM6QUut9au9h4fDbwJfNRa2xbYbyUwxVp7ySDn84OULwywz3zg\nn/PRfhERkTHqKmvtw/k8YSh6UgYxHigH/pCy/Q/AewY72Fr7YBav0Q5cBfwO+NMQ2yciIjKWHQe8\nD/ddmlejIUjJJG9dQNbaA0Beoz8REZEx5IlCnDQsibMDeRWI07/X5L30710RERGREhH6IMVa+w7Q\nCdSnPHUJBYrcREREpPhCMdxjjDkRmEiinslZxpjzgB5r7UvAHcBPjDFPAk8B1wDVwOXFaK+IiIgU\nXihm9xhjZgOP0j/P5EFr7dXePtfh6qS8B/gt8Dlr7X+MaENFRERkxIQiSBERERFJFfqclLAwxvyr\nMabHGPOzYrdFhscY89fGmN8aY54zxny62O2R4dPvZ+kwxpxmjNlsjHnW+x29qthtktwZY04yxjxt\njOkyxuwyxnx+SMerJyU7xphZwMnAp6y1VxS7PZIbY8yxuOHCabgigTuAi621mik2iun3s3QYYyYA\np1prdxlj3oVbMuV/WmtfK3LTJAdeFfljrLVvG2OOA54HLrLW/j6b49WTkiVr7Wbg9WK3Q4btr4Au\na+0Ba+1bwL8DDUVukwyTfj9Lh7X2FWvtLu/+fuD3uKKeMgpZ523v4QnAnxnC76qCFBlrTgNeDjx+\n2dsmIiFjjJmK+1/4i8Vui+TOGFNhjNkJ7APuGkqvWEkGKcaYmcaY1caYl40xvcaYeWn2ucEYs9cY\n85YxZrsxZkYx2irZ03UtPbqmpSWf19MYUwU8iCs5IUWQr+tprT1srZ0EnAksNMacnW0bSjJIAU7E\njWPeQJry+caYK4Hv4aY01wLrgbXGmNMC+yw2xuzwkn1CUU9Ghn9d6d9zchrwn4VqsAwqH9dUwiMv\n19MYcwzwr8C3rLVPFrrRklFefz+94bvHgKlZt8BaW9I3oBeYl7JtG7AyZdsO3C/EQOe6GPh5sd+T\nbrlfV+BYYDdujPsEYBfw3mK/H92G/7uq389w3YZzPYEo8LVivwfdhn89gXcBJ3j3T8FNXPhgtq9b\nqj0pGRljjgamABtTnvo1cOEAx60DHgHmGmP2GWOmFK6VMlTZXlfrErg+h4vmt+N+wbLKMpeRNZTf\nVf1+hl+219MYcxHwCeDyQG/2+0eupZKNIfx+ngE84eWkbAXusdY+m+3rjMVhjPFAOf0XJ/wD/Rcx\n7GOtnVPIRsmwZX1drbW/An41Qu2S3A3lmur3M/yyup7WVRIfi99No02217MTmJTri4y5npQBqGBM\nadJ1LT26pqVF17O05PV6jsUg5VUgTv9ek/fSPyKU0UPXtfTompYWXc/SMiLXc8wFKdbad4BOoD7l\nqUuAJ0a+RZIPuq6lR9e0tOh6lpaRup4lOe5njDkRmAgYb9NZxpjzgB5r7UvAHcBPjDFPAk/h5uFX\nA5cXo72SHV3X0qNrWlp0PUtLKK5nsac1FWiq1GzcdKl4yu2BwD7XAS8Cb+FmeVxU7Hbrpus61m66\npqV10/UsrVsYrqcWGBQREZFQGnM5KSIiIjI6KEgRERGRUFKQIiIiIqGkIEVERERCSUGKiIiIhJKC\nFBEREQklBSkiIiISSgpSREREJJQUpIiIiEgoKUgRERGRUFKQIiIFZYx51BhzR7HbISKjj4IUEZH/\n1769hFpVxXEc//4qQySQmlwhnYT0sNSwB5UDEyoteiglEUEmEVSDiBoIPTAHFRR3lN6C6EGDglKI\nIAkKqkmhkBKRSqVBJXeQEOIdaIH/Bntf2By8oXj07uj7GZ2z/mutve6BCz/WXktSLxlSJP2nJTln\nutcg6fQwpEgamiSzkryb5HCSA0meHKjPSPJykt+THEmyK8mtA30eTvJrkokkHyR5IsmfnfqGdty6\nJPuAI53a+iT7khxNsjfJAwNzL0jySbu+Q0k+TDJymn4OSafIkCJpmMaApcAK4EbgBmBJp/4OcDWw\nCpgPvAVsTXIFQJJbgNeAl4DLgG3AM0ANPGc+cCdwO7CwHfsCcB+wDrgI2AiMJVnZ1ucAXwLfAIuA\n64GZwEdD+cslDV2qBv/3JenkJZkN/AHcU1Uft23nAePAG8Am4EdgTlUd7IzbBuypqqeSbAWOVdWa\nTv1NYHVVXdB+3wCsB0aq6nDbNgs4CCytql2dsWNtv7uTbASurKq7OvWRdn2Lq+r74f8qkk6F73Il\nDcslwNnA9smGqppI8kP7dQnN7u3+JOmMOxc43H5eALw9MO8OYPVA2/7JgNIZNxP4amDuGcB37eer\ngJuTdMdBs0szDzCkSD1jSJE0LJPhYKrt2bOAo8DiTt9JEyf5rL+OMzfATTQ7Kl1HO322AM8e5/nj\nJ/l8SWeAIUXSsOwFjgHXAd3XPZcDXwO7aHZNzq+qnVPMsRu4ZqDt2hN49m6aA7Rzq2rHFH12AndU\n1S8nMJ+kHvDgrKShqKpDwHvAaJKlSS6mORh7rK3/1NbfT3JbkguTLGpv76xop3kdWJXk0STzkqyl\nOSD7r4fnqmoCGAU2J7k3ydwkl7Y3he5vu20GRtrbR4vaPsuTjA28IpLUE4YUScP0GM2uyac0N2m2\nA9926g/SBJVXgf3A58By4ABAVX0GPAI8DeyhCSiv0LlmPJWqeg54EXge+Lldxxrgt7Y+TnPzaBbw\nBbCPJhT9Xd4gkHrJ2z2Sei3JJmBhVS2b7rVIOrM8kyKpV5I8TrPDMgGsBB6i2aGR9D/jToqkXkmy\nBVgGzKZ5VTNWVaPTuypJ08GQIkmSesmDs5IkqZcMKZIkqZcMKZIkqZcMKZIkqZcMKZIkqZcMKZIk\nqZcMKZIkqZcMKZIkqZf+AfVdcKoMhveOAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f6dbb4921d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot\n",
    "ax = matplotlib.pyplot.gca()\n",
    "ax.set_xscale(\"log\")\n",
    "ax.set_yscale(\"log\")\n",
    "degrees = numpy.array(neph_graph.indegree())\n",
    "inds_keep = numpy.where(degrees > 0)\n",
    "ax.scatter(degrees[inds_keep], \n",
    "           pageranks[inds_keep])\n",
    "#ax.scatter(neph_graph.indegree(), pageranks)\n",
    "matplotlib.pyplot.xlabel(\"degree\")\n",
    "matplotlib.pyplot.ylabel(\"pagerank\")\n",
    "matplotlib.pyplot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See if you can calculate the pagerank centrality yourself, using the matrix inversion (Eq. 7.19 from Newman). Test your function on a small directed graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pagerank(g):\n",
    "    # N is the number of vertices\n",
    "    N = len(g.vs)\n",
    "    \n",
    "    # alpha is the damping parameter\n",
    "    alpha = 0.85\n",
    "    \n",
    "    # beta = (1-alpha)/N\n",
    "    beta = (1-alpha)/N\n",
    "    \n",
    "    # compute the out-degree of each vertex\n",
    "    degree_values = g.degree(mode=\"out\")\n",
    "    \n",
    "    # get a floating-point adjacency matrix M in the Newman format (take transpose from igraph format)\n",
    "    M = numpy.matrix(g.get_adjacency().data).transpose().astype(float)\n",
    "    \n",
    "    # or each column in 0,N-1:\n",
    "    for j in range(0,N):\n",
    "        \n",
    "        # get the out degree of the vertex as a float\n",
    "        degree_value = float(degree_values[j])\n",
    "        \n",
    "        # if degree is nonzero, normalize the column of M\n",
    "        if degree_value > 0:\n",
    "            M[:,j] /= degree_value\n",
    "        else:\n",
    "            # set the column to zero\n",
    "            M[:,j] = 0\n",
    "            \n",
    "    # use Equation on the board\n",
    "    pr = numpy.linalg.inv(numpy.diag([1.]*N) - alpha * M) * numpy.matrix([beta]*N).transpose()\n",
    "    pr = pr/numpy.sum(pr)\n",
    "    \n",
    "    return(pr.transpose())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test out your function on a small 5-vertex graph:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.2778626  0.1480916  0.2778626  0.1480916  0.1480916]]\n",
      "[0.2778625954198473, 0.14809160305343513, 0.2778625954198473, 0.14809160305343513, 0.1480916030534351]\n"
     ]
    }
   ],
   "source": [
    "g = igraph.Graph.Barabasi(n=5, m=2) \n",
    "print(pagerank(g))\n",
    "print(g.pagerank())"
   ]
  }
 ],
 "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}