{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# PageRank - A numpy / Jupyter / matplotlib example\n",
"\n",
"In this Jupyter notebook we consider Google's PageRank computation. Nodes correspond to webpages and edges to links from one webpage to another webpage. We consider a very simple graph with only six _nodes_ and where every node has one or two outgoing _edges_. The original description of the PageRank computation can be found in the paper below that contains an overview of the original infrastructure of the Google search engine.\n",
"\n",
"* Sergey Brin and Lawrence Page.\n",
" _The Anatomy of a Large-Scale Hypertextual Web Search Engine_.\n",
" Seventh International World-Wide Web Conference (WWW 1998), April 14-18, 1998, Brisbane, Australia.\n",
" [http://ilpubs.stanford.edu:8090/361/]\n",
" \n",
"\n",
"\n",
"## Random surfer (simplified)\n",
"The PageRank of a node (web page) is the ratio of the time one visits a node by performing an _infinite random traversal_ of a graph where one starts in node 1, and in each step performs:\n",
"\n",
" * with probability 1/6 jumps to a random page (probability 1/6 for each node)\n",
" * with probability 5/6 follows an outgoing edge to an adjacent node (selected uniformly)\n",
" \n",
"The above can be simulated by using a dice: Roll a dice. If it shows 6, jump to a random page by rolling the dice again to figure out which node to jump to. If the dice shows 1-5, follow an outgoing edge - if two outgoing edges roll the dice again and go to the lower number neighbor if it is odd."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Graph represenation - _adjacency matrix_"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[0 1 0 0 0 0]\n",
" [0 0 0 1 0 0]\n",
" [1 1 0 0 0 0]\n",
" [0 1 0 0 1 0]\n",
" [0 1 0 0 0 1]\n",
" [0 1 0 0 0 0]]\n"
]
}
],
"source": [
"import numpy as np\n",
"\n",
"# Adjacency matrix of the above directed graph \n",
"# (note that the rows/colums are 0-indexed, whereas in the figure the nodes are 1-indexed)\n",
"\n",
"G = np.array([[0, 1, 0, 0, 0, 0],\n",
" [0, 0, 0, 1, 0, 0],\n",
" [1, 1, 0, 0, 0, 0],\n",
" [0, 1, 0, 0, 1, 0],\n",
" [0, 1, 0, 0, 0, 1],\n",
" [0, 1, 0, 0, 0, 0]])\n",
"\n",
"n = G.shape[0] # number of rows in G\n",
"\n",
"degree = np.sum(G, axis=1, keepdims=True) # column vector with row sums = out-degrees\n",
"\n",
"# The below code handles sinks, i.e. nodes with outdegree zero\n",
"# this has no effect on the graph above, since no sinks\n",
"\n",
"G = G + (degree == 0) # add edges from sinks to all nodes \n",
"degree = np.sum(G, axis=1, keepdims=True)\n",
"\n",
"print(G)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Simulate random walk (random surfer model)"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[1], [3], [0, 1], [1, 4], [1, 5], [1]]\n",
"[0.039108 0.353342 0.027621 0.322408 0.162235 0.095286]\n"
]
}
],
"source": [
"from random import randint\n",
"\n",
"STEPS = 1000000\n",
"\n",
"# adjacency_list[i] is a list of all j where (i, j) is an edge of the graph.\n",
"adjacency_list = [[j for j, e in enumerate(row) if e] for row in G]\n",
"\n",
"count = np.zeros(n) # histogram over number of node visits\n",
"state = 0 # start at node with index 0\n",
"for _ in range(STEPS):\n",
" count[state] += 1 # increment count for state\n",
" if randint(1, 6) == 6: # original paper uses 15% instead of 1/6\n",
" state = randint(0, 5)\n",
" else:\n",
" state = adjacency_list[state][randint(0, degree[state] - 1)]\n",
"\n",
"print(adjacency_list, count / STEPS, sep='\\n')"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEWCAYAAABbgYH9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAfT0lEQVR4nO3df7RdZX3n8ffHhGIqEAJEB5PYoERHoBokBjpYi2KTVFqDU9DYjqQ2bSzFX8uu1YH+AmHSgamCCy3YWCIJY4GIpaQqxQgCpcbADQ2E8GO4I1FiIokmhqAlTuJn/tjP1Z3LuTcnyd335J58Xmuddfb57v0857txeb959vOcvWWbiIiIofaiTicQERHdKQUmIiIakQITERGNSIGJiIhGpMBEREQjUmAiIqIRKTARw0jSGZLWdzqPvSXpekn/o2yPyHOI4ZcCEwc9Sesk/Yek5yR9r/wxPazTee0PSX8n6Zra50Mk/WiA2GmdyTK6XQpMROW3bB8GTAVOBi7qcD77617g12qfpwHfAd7cLwawariSioNLCkxEje3vAXdQFRoAJJ0l6d8lPSvpaUmX1PZNlmRJcyV9R9L3Jf15bf+YMiLaKulR4I3175P0Wkl3S/qhpLWS3lHbd72kayTdXkZX/ybpP0n6ZOnvcUknD3Aq9wCvlXRM+fyrwE3AS/rFVtj+f+X7vlBGcNsk3SvpxHb+m0n6kKRHJU1s5/g4eKTARNSUP5K/AfTWwj8CzgOOBM4Czpd0dr+mbwJeA5wJ/JWk15b4xcCrymsmMLf2XYcA/wx8FXgp8EHg85JeU+v3XcBfAMcAO4AVwIPl8y3Ala3Ow/Z64NtURQSqkcu/At/oF7u31ux2YErJ5UHg8636rpP0l8DvAb9WvjPiZ1JgIir/JGk78DSwiaowAGD7bttrbP/U9sPAjex++QngY7b/w/ZDwEPA60v8XcAC21tsPw1cXWtzGnAYcLntn9i+C/gS8J7aMbfaXmX7eeBW4HnbS2zvAm6mupw3kHuAN0t6ETAd+CZVkemLnV6O6TvPRba3294BXAK8XtLYAfqWpCupiuZbbG8eJI84SKXARFTOtn04cAbwn6lGCABIOlXS1yVtlrQN+KP6/uJ7te0fUxUOgJdTFa0+365tvxx42vZP++2fUPv8TG37P1p8Hmwxwr1Uo5RfBr5l+8fAfbXYGGBlOcdRki6X9H8lPQusK330P88+RwLzgf9pe9sgOcRBLAUmosb2PcD1wMdr4X8AlgGTbI8FPgOozS43ApNqn19R294ATCqjifr+7+5l2gO5l2okdRbVyAVgbcnnLOCBMjIC+B1gNvA2YCwwucQHOs+twG8Cn5N0+hDlG10mBSbihT4J/Lqkvon+w4Ettp+XNJ3qj3G7lgIXSRpX5nc+WNu3kmp+50/LkuEzgN+imozfb7Z7qUY8H6YUGFfP51hZYvX5l8Op5nh+APwi8Ndt9H838LvArZJOHYqco7ukwET0U+YTlgB/WUJ/DFxa5mj+iqpotOtjVJe9nqKazL+h9j0/Ad5Btajg+8A1wHm2H9/fc6i5FxgP/Fst9q9UE/n1ArOk5Pld4FGq+Zo9sr0ceB+wTNIpQ5FwdA/lgWMREdGEjGAiIqIRKTAREdGIFJiIiGhECkxERDRidKcTOFAcc8wxnjx5cqfTiIgYUVatWvV92+Nb7UuBKSZPnkxPT0+n04iIGFEkfXugfblEFhERjUiBiYiIRqTAREREI1JgIiKiESkwERHRiBSYiIhoRApMREQ0orECI+nFku6X9JCktZI+VuKXSPqupNXl9fZam4sk9Up6QtLMWvwUSWvKvqslqcQPlXRzia+UNLnWZq6kJ8trLhERMaya/KHlDuCttp+TdAhwn6Tby76rbNefGIikE4A5wIlUj5L9mqRXl2ePX0v1eNZvAl8BZgG3A/OArbaPlzQHuAJ4t6SjqJ6pPg0wsErSMttbGzzfiIioaazAlCfnPVc+HlJegz18ZjZwk+0dwFOSeoHpktYBR9heASBpCXA2VYGZDVxS2t8CfLqMbmYCy21vKW2WUxWlG4fsBLvc5Au/3OkU2rLu8rM6nUJEDKDRORhJoyStBjZR/cFfWXZ9QNLDkhZJGldiE4Cna83Xl9iEst0/vlsb2zuBbcDRg/QVERHDpNECY3uX7anARKrRyElUl7teBUwFNgKfKIerVReDxPe1zc9Imi+pR1LP5s2bBz2XiIjYO8Oyisz2D4G7gVm2nymF56fAZ4Hp5bD1wKRas4nAhhKf2CK+WxtJo4GxwJZB+uqf10Lb02xPGz++5c1AIyJiHzW5imy8pCPL9hjgbcDjko6tHfZO4JGyvQyYU1aGHQdMAe63vRHYLum0Mr9yHnBbrU3fCrFzgLvK3M8dwAxJ48oluBklFhERw6TJVWTHAosljaIqZEttf0nSDZKmUl2yWge8H8D2WklLgUeBncAFZQUZwPnA9cAYqsn9vtVo1wE3lAUBW6hWoWF7i6TLgAfKcZf2TfhHRMTwaHIV2cPAyS3i7x2kzQJgQYt4D3BSi/jzwLkD9LUIWLQXKUdExBDKL/kjIqIRKTAREdGIFJiIiGhECkxERDQiBSYiIhqRAhMREY1IgYmIiEakwERERCNSYCIiohEpMBER0YgUmIiIaESTN7uMiAblqaNxoMsIJiIiGpECExERjUiBiYiIRqTAREREI1JgIiKiESkwERHRiBSYiIhoRApMREQ0IgUmIiIa0ViBkfRiSfdLekjSWkkfK/GjJC2X9GR5H1drc5GkXklPSJpZi58iaU3Zd7Uklfihkm4u8ZWSJtfazC3f8aSkuU2dZ0REtNbkCGYH8FbbrwemArMknQZcCNxpewpwZ/mMpBOAOcCJwCzgGkmjSl/XAvOBKeU1q8TnAVttHw9cBVxR+joKuBg4FZgOXFwvZBER0bzGCowrz5WPh5SXgdnA4hJfDJxdtmcDN9neYfspoBeYLulY4AjbK2wbWNKvTV9ftwBnltHNTGC57S22twLL+XlRioiIYdDoHIykUZJWA5uo/uCvBF5meyNAeX9pOXwC8HSt+foSm1C2+8d3a2N7J7ANOHqQvvrnN19Sj6SezZs378+pRkREP40WGNu7bE8FJlKNRk4a5HC16mKQ+L62qee30PY029PGjx8/SGoREbG3hmUVme0fAndTXaZ6plz2orxvKoetBybVmk0ENpT4xBbx3dpIGg2MBbYM0ldERAyTJleRjZd0ZNkeA7wNeBxYBvSt6poL3Fa2lwFzysqw46gm8+8vl9G2SzqtzK+c169NX1/nAHeVeZo7gBmSxpXJ/RklFhERw6TJB44dCywuK8FeBCy1/SVJK4ClkuYB3wHOBbC9VtJS4FFgJ3CB7V2lr/OB64ExwO3lBXAdcIOkXqqRy5zS1xZJlwEPlOMutb2lwXONiIh+Giswth8GTm4R/wFw5gBtFgALWsR7gBfM39h+nlKgWuxbBCzau6wjImKo5Jf8ERHRiBSYiIhoRApMREQ0IgUmIiIakQITERGNSIGJiIhGpMBEREQjUmAiIqIRKTAREdGIFJiIiGhECkxERDQiBSYiIhqRAhMREY1IgYmIiEakwERERCNSYCIiohEpMBER0YgUmIiIaEQKTERENCIFJiIiGtFYgZE0SdLXJT0maa2kD5f4JZK+K2l1eb291uYiSb2SnpA0sxY/RdKasu9qSSrxQyXdXOIrJU2utZkr6cnymtvUeUZERGujG+x7J/Anth+UdDiwStLysu8q2x+vHyzpBGAOcCLwcuBrkl5texdwLTAf+CbwFWAWcDswD9hq+3hJc4ArgHdLOgq4GJgGuHz3MttbGzzfiIio2asRjKRxkl7XzrG2N9p+sGxvBx4DJgzSZDZwk+0dtp8CeoHpko4FjrC9wraBJcDZtTaLy/YtwJlldDMTWG57Sykqy6mKUkREDJM9FhhJd0s6oowKHgI+J+nKvfmScunqZGBlCX1A0sOSFkkaV2ITgKdrzdaX2ISy3T++WxvbO4FtwNGD9NU/r/mSeiT1bN68eW9OKSIi9qCdEcxY288C/xX4nO1TgLe1+wWSDgO+CHyk9HMt8CpgKrAR+ETfoS2ae5D4vrb5ecBeaHua7Wnjx48f9DwiImLvtFNgRpfLVO8CvrQ3nUs6hKq4fN72PwLYfsb2Lts/BT4LTC+Hrwcm1ZpPBDaU+MQW8d3aSBoNjAW2DNJXREQMk3YKzMeAO4Be2w9IeiXw5J4albmQ64DHbF9Zix9bO+ydwCNlexkwp6wMOw6YAtxveyOwXdJppc/zgNtqbfpWiJ0D3FXmae4AZpQ5o3HAjBKLiIhh0s4qso22fzaxb/tbbc7BnA68F1gjaXWJ/RnwHklTqS5ZrQPeX/pdK2kp8CjVCrQLygoygPOB64ExVKvHbi/x64AbJPVSjVzmlL62SLoMeKAcd6ntLW3kHBERQ6SdAvMp4A1txHZj+z5az4V8ZZA2C4AFLeI9wEkt4s8D5w7Q1yJg0WA5RkREcwYsMJJ+BfgvwHhJH63tOgIY1XRiERExsg02gvkF4LByzOG1+LNU8x0REREDGrDA2L4HuEfS9ba/PYw5RUREFxjsEtknbX8E+LSkVr8heUejmUVExIg22CWyG8r7xwc5JiIioqXBLpGtKu/39MXKb0om2X54GHKLiIgRbFjuRRYREQefxu9FFhERB6dG70UWEREHr3YKzKXsw73IIiLi4LbHW8XY/gLwhdrnbwG/3WRSEREx8g32O5g/tf2/JH2K1s9S+VCjmUVExIg22AjmsfLeMxyJREREdxnsdzD/XDYftv3vw5RPRER0iXYm+a+U9LikyySd2HhGERHRFfZYYGy/BTgD2AwslLRG0l80nVhERIxs7YxgsP0921cDfwSsBv6q0awiImLEa+dWMa+VdImkR4BPA98AJjaeWUREjGjtPDL5c8CNwAzbGxrOJyIiukQ7P7Q8bTgSiYiI7tLWHMy+kDRJ0tclPSZpraQPl/hRkpZLerK8j6u1uUhSr6QnJM2sxU8piwt6JV0tSSV+qKSbS3ylpMm1NnPLdzwpaW5T5xkREa01VmCAncCf2H4tcBpwgaQTgAuBO21PAe4snyn75gAnArOAaySNKn1dC8wHppTXrBKfB2y1fTxwFXBF6eso4GLgVGA6cHG9kEVERPMGLDCSbijvH96Xjm1vtP1g2d5OdWeACcBsYHE5bDFwdtmeDdxke4ftp4BeYHq5k/MRtlfYNrCkX5u+vm4Bziyjm5nActtbbG8FlvPzohQREcNgsBHMKZJ+Cfh9SePKpa2fvfbmS8qlq5OBlcDLbG+EqggBLy2HTQCerjVbX2ITynb/+G5tbO8EtgFHD9JXREQMk8Em+T8D/AvwSmAVoNo+l/geSToM+CLwEdvPlumTloe2iHmQ+L62qec2n+rSG694xSsGyisiIvbBgCMY21eX+ZNFtl9p+7jaq93icghVcfm87X8s4WfKZS/K+6YSXw9MqjWfCGwo8Ykt4ru1kTQaGAtsGaSv/ue40PY029PGjx/fzilFRESb2rlVzPmSXi/pA+X1unY6LnMh1wGP2b6ytmsZ0Leqay5wWy0+p6wMO45qMv/+chltu6TTSp/n9WvT19c5wF1lnuYOYEa5tDcOmFFiERExTPb4OxhJH6K6jNQ3Avm8pIW2P7WHpqcD7wXWSFpdYn8GXA4slTQP+A5wLoDttZKWAo9SrUC7wPau0u584HpgDHB7eUFVwG6Q1Es1cplT+toi6TLggXLcpba37OlcIyJi6LTzS/4/AE61/SMASVcAK4BBC4zt+2g9FwJw5gBtFgALWsR7gJNaxJ+nFKgW+xYBiwbLMSIimtPO72AE7Kp93sXAhSMiIgJo/15kKyXdWj6fTXVpKiIiYkDt3IvsSkl3A2+iGrm8L0+4jIiIPWlnBEP5Rf6DDecSERFdpMl7kUVExEEsBSYiIhoxaIGRNErS14YrmYiI6B6DFpjyQ8cfSxo7TPlERESXaGeS/3mqX+MvB37UF7T9ocayioiIEa+dAvPl8oqIiGhbO7+DWSxpDPAK208MQ04REdEF9riKTNJvAaupng2DpKmSljWdWEREjGztLFO+hOq59j8EsL0aOK7BnCIiogu0U2B22t7WL/aCp0NGRETUtTPJ/4ik3wFGSZoCfAj4RrNpRUTESNfOCOaDwInADuBG4FngI00mFRERI187q8h+DPx5edCYbW9vPq2IiBjp2nlk8hupngx5ePm8Dfh926sazi0iDjKTLxwZP7lbd/lZnU5hRGhnDuY64I9t/yuApDdRPYTsdU0mFhERI1s7czDb+4oLgO37gFwmi4iIQQ04gpH0hrJ5v6S/o5rgN/Bu4O7mU4uIiJFssBHMJ8prKvBq4GKqH12+FviVPXUsaZGkTZIeqcUukfRdSavL6+21fRdJ6pX0hKSZtfgpktaUfVdLUokfKunmEl8paXKtzVxJT5bX3Db/W0RExBAacARj+y372ff1wKeBJf3iV9n+eD0g6QRgDtVy6JcDX5P06vK4gGuB+cA3ga8As4DbgXnAVtvHS5oDXAG8W9JRVMVwGtWIa5WkZba37uf5RETEXmhnFdmRwHnA5Prxe7pdv+1766OKPZgN3GR7B/CUpF5guqR1wBG2V5RclgBnUxWY2VQjKoBbgE+X0c1MYLntLaXNcqqidGObuURExBBoZ5L/K1TFZQ2wqvbaVx+Q9HC5hDauxCYAT9eOWV9iE8p2//hubWzvBLYBRw/S1wtImi+pR1LP5s2b9+OUIiKiv3aWKb/Y9keH6PuuBS6junR1GdUcz+8DanGsB4mzj212D9oLgYUA06ZNy/3VIiKGUDsjmBsk/aGkYyUd1ffaly+z/YztXbZ/CnyW6i7NUI0yJtUOnQhsKPGJLeK7tZE0GhgLbBmkr4iIGEbtFJifAH8DrODnl8d69uXLJB1b+/hOoG+F2TJgTlkZdhwwBbjf9kZgu6TTyvzKecBttTZ9K8TOAe6ybeAOYIakceUS3IwSi4iIYdTOJbKPAsfb/v7edCzpRuAM4BhJ66lWdp0haSrVJat1wPsBbK+VtBR4FNgJXFBWkAGcT7UibQzV5P7tJX4d1eiql2rkMqf0tUXSZcAD5bhL+yb8IyJi+LRTYNYCP97bjm2/p0X4ukGOXwAsaBHvAU5qEX8eOHeAvhZR3T8tIiI6pJ0CswtYLenrVLfsB/a8TDkiIg5u7RSYfyqviIiItrXzPJjFw5FIRER0l3Z+yf8ULX5HYvuVjWQUERFdoZ1LZNNq2y+mmljfp9/BRETEwWOPv4Ox/YPa67u2Pwm8dRhyi4iIEaydS2RvqH18EdWI5vDGMoqIiK7QziWyT9S2d1L9QPJdjWQTERFdo51VZPv7XJiIiDgItXOJ7FDgt3nh82AubS6tiIgY6dq5RHYb1bNWVlH7JX9ERMRg2ikwE23PajyTiIjoKu3crv8bkn658UwiIqKrtDOCeRPwe+UX/Tuonhhp269rNLOIiBjR2ikwv9F4FhER0XXaWab87eFIJCIiuks7czARERF7LQUmIiIakQITERGNSIGJiIhGNFZgJC2StEnSI7XYUZKWS3qyvI+r7btIUq+kJyTNrMVPkbSm7Ltakkr8UEk3l/hKSZNrbeaW73hS0tymzjEiIgbW5AjmeqD/HQAuBO60PQW4s3xG0gnAHODE0uYaSaNKm2uB+cCU8urrcx6w1fbxwFXAFaWvo4CLgVOB6cDF9UIWERHDo7ECY/teYEu/8GxgcdleDJxdi99ke4ftp4BeYLqkY4EjbK+wbWBJvzZ9fd0CnFlGNzOB5ba32N4KLOeFhS4iIho23HMwL7O9EaC8v7TEJwBP145bX2ITynb/+G5tbO+kuiHn0YP09QKS5kvqkdSzefPm/TitiIjo70CZ5FeLmAeJ72ub3YP2QtvTbE8bP358W4lGRER7hrvAPFMue1HeN5X4emBS7biJwIYSn9givlsbSaOBsVSX5AbqKyIihlE79yIbSsuAucDl5f22WvwfJF0JvJxqMv9+27skbZd0GrASOA/4VL++VgDnAHfZtqQ7gL+uTezPAC5q/tQiInY3+cIvdzqFtqy7/KxG+m2swEi6ETgDOEbSeqqVXZcDSyXNA74DnAtge62kpcCjwE7gAtu7SlfnU61IGwPcXl4A1wE3SOqlGrnMKX1tkXQZ8EA57lLb/RcbREREwxorMLbfM8CuMwc4fgGwoEW8BzipRfx5SoFqsW8RsKjtZCMiYsgdKJP8ERHRZVJgIiKiESkwERHRiBSYiIhoRApMREQ0IgUmIiIakQITERGNSIGJiIhGpMBEREQjUmAiIqIRKTAREdGIFJiIiGhECkxERDQiBSYiIhqRAhMREY1IgYmIiEakwERERCNSYCIiohEpMBER0YgUmIiIaERHCoykdZLWSFotqafEjpK0XNKT5X1c7fiLJPVKekLSzFr8lNJPr6SrJanED5V0c4mvlDR5uM8xIuJg18kRzFtsT7U9rXy+ELjT9hTgzvIZSScAc4ATgVnANZJGlTbXAvOBKeU1q8TnAVttHw9cBVwxDOcTERE1B9IlstnA4rK9GDi7Fr/J9g7bTwG9wHRJxwJH2F5h28CSfm36+roFOLNvdBMREcOjUwXGwFclrZI0v8ReZnsjQHl/aYlPAJ6utV1fYhPKdv/4bm1s7wS2AUf3T0LSfEk9kno2b948JCcWERGV0R363tNtb5D0UmC5pMcHObbVyMODxAdrs3vAXggsBJg2bdoL9kdExL7ryAjG9obyvgm4FZgOPFMue1HeN5XD1wOTas0nAhtKfGKL+G5tJI0GxgJbmjiXiIhobdgLjKSXSDq8bxuYATwCLAPmlsPmAreV7WXAnLIy7Diqyfz7y2W07ZJOK/Mr5/Vr09fXOcBdZZ4mIiKGSScukb0MuLXMuY8G/sH2v0h6AFgqaR7wHeBcANtrJS0FHgV2AhfY3lX6Oh+4HhgD3F5eANcBN0jqpRq5zBmOE4uIiJ8b9gJj+1vA61vEfwCcOUCbBcCCFvEe4KQW8ecpBSoiIjrjQFqmHBERXaRTq8i6zuQLv9zpFNqy7vKzOp1CRBwkMoKJiIhGpMBEREQjUmAiIqIRKTAREdGIFJiIiGhECkxERDQiy5TjoJGl5BHDKyOYiIhoRApMREQ0IgUmIiIakQITERGNSIGJiIhGpMBEREQjUmAiIqIRKTAREdGIFJiIiGhECkxERDQiBSYiIhqRAhMREY3o6gIjaZakJyT1Srqw0/lERBxMurbASBoF/C3wG8AJwHskndDZrCIiDh5dW2CA6UCv7W/Z/glwEzC7wzlFRBw0ZLvTOTRC0jnALNt/UD6/FzjV9gdqx8wH5pePrwGeGPZEB3cM8P1OJzGEuu18oPvOqdvOB7rvnA608/kl2+Nb7ejmB46pRWy3amp7IbBweNLZe5J6bE/rdB5DpdvOB7rvnLrtfKD7zmkknU83XyJbD0yqfZ4IbOhQLhERB51uLjAPAFMkHSfpF4A5wLIO5xQRcdDo2ktktndK+gBwBzAKWGR7bYfT2lsH7OW7fdRt5wPdd07ddj7Qfec0Ys6nayf5IyKis7r5EllERHRQCkxERDQiBeYA1G23uJG0SNImSY90OpehIGmSpK9LekzSWkkf7nRO+0vSiyXdL+mhck4f63ROQ0HSKEn/LulLnc5lKEhaJ2mNpNWSejqdz55kDuYAU25x83+AX6daav0A8B7bj3Y0sf0g6c3Ac8AS2yd1Op/9JelY4FjbD0o6HFgFnD3C/zcS8BLbz0k6BLgP+LDtb3Y4tf0i6aPANOAI27/Z6Xz2l6R1wDTbB9IPLQeUEcyBp+tucWP7XmBLp/MYKrY32n6wbG8HHgMmdDar/ePKc+XjIeU1ov/1KWkicBbw953O5WCVAnPgmQA8Xfu8nhH+x6ubSZoMnAys7Gwm+69cTloNbAKW2x7p5/RJ4E+Bn3Y6kSFk4KuSVpVbXR3QUmAOPHu8xU0cGCQdBnwR+IjtZzudz/6yvcv2VKq7XkyXNGIvZ0r6TWCT7VWdzmWInW77DVR3ib+gXH4+YKXAHHhyi5sRoMxTfBH4vO1/7HQ+Q8n2D4G7gVkdTmV/nA68o8xZ3AS8VdL/7mxK+8/2hvK+CbiV6pL6ASsF5sCTW9wc4MqE+HXAY7av7HQ+Q0HSeElHlu0xwNuAxzub1b6zfZHtibYnU/1/6C7b/63Dae0XSS8pi0qQ9BJgBnBAr8xMgTnA2N4J9N3i5jFg6Qi8xc1uJN0IrABeI2m9pHmdzmk/nQ68l+pfxavL6+2dTmo/HQt8XdLDVP/IWW67K5b2dpGXAfdJegi4H/iy7X/pcE6DyjLliIhoREYwERHRiBSYiIhoRApMREQ0IgUmIiIakQITERGNSIGJGGEkTe6WO1NHd0uBiYiIRqTARHRAGYU8Jumz5fkrX5U0RtJUSd+U9LCkWyWNK8efUp7VsgK4oNbPKEl/I+mB0ub9HTupiH5SYCI6Zwrwt7ZPBH4I/DawBPjvtl8HrAEuLsd+DviQ7V/p18c8YJvtNwJvBP5Q0nHDkn3EHqTARHTOU7ZXl+1VwKuAI23fU2KLgTdLGtsvfkOtjxnAeeU2+yuBo6kKV0THje50AhEHsR217V3AkQMcJwZ+ZIOAD9q+YygTixgKGcFEHDi2AVsl/Wr5/F7gnnL7/G2S3lTiv1trcwdwfnl8AJJeXe60G9FxGcFEHFjmAp+R9IvAt4D3lfj7gEWSfkxVVPr8PTAZeLA8RmAzcPbwpRsxsNxNOSIiGpFLZBER0YgUmIiIaEQKTERENCIFJiIiGpECExERjUiBiYiIRqTAREREI/4/9mP6K/As51cAAAAASUVORK5CYII=\n",
"text/plain": [
"