{
 "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",
    "<img src=\"pagerank.png\" height=\"300\" width=\"300\">\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": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# To make matploblib figures interactive replace 'inline' by 'notebook'\n",
    "%matplotlib inline\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.bar(range(6), count)\n",
    "\n",
    "plt.title(\"Random Walk\")\n",
    "plt.xlabel(\"node\")\n",
    "plt.ylabel(\"number of visits\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Transition matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.  1.  0.  0.  0.  0. ]\n",
      " [0.  0.  0.  1.  0.  0. ]\n",
      " [0.5 0.5 0.  0.  0.  0. ]\n",
      " [0.  0.5 0.  0.  0.5 0. ]\n",
      " [0.  0.5 0.  0.  0.  0.5]\n",
      " [0.  1.  0.  0.  0.  0. ]]\n"
     ]
    }
   ],
   "source": [
    "A = G / degree  # Normalize row-sums to one. Note that 'degree' \n",
    "                # is a n x 1 matrix, whereas G is an n x n matrix.\n",
    "                # The elementwise division is repeated for each column of G.\n",
    "       \n",
    "print(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computation of transition matrix\n",
    "\n",
    "We now want to compute the probability $p^{(0)}_j$ to be in vertex $j$ after $i$ steps. Initially we have $p^{(0)}_0=1$ and $p^{(0)}_j=0$ for $j\\neq 0$. Let $p^{(i)}=(p^{(i)}_0,\\ldots,p^{(i)}_{n-1})$, i.e. $p^{(0)} = (1,0,\\ldots,0)$. We compute a matrix $M$, such that $p^{(i)} = M^i p$ (assuming p is a column vector). If we let $\\textbf{1}_n$ denote the $n\\times n$ matrix with one in each entry, then $M$ can be computed as:\n",
    "\n",
    "$$p^{(i+1)}_j = \\frac{1}{6} \\cdot \\frac{1}{n} + \\frac{5}{6} \\sum_k p^{(i)}_{k} \\cdot A_{k, j}$$\n",
    "\n",
    "$$p^{(i+1)} = M \\cdot p^{(i)}$$\n",
    "\n",
    "$$M = \\frac{1}{6}\\cdot\\frac{1}{n}\\textbf{1}_n + \\frac{5}{6} A^\\mathrm{T}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.03935185]\n",
      " [0.35326184]\n",
      " [0.02777778]\n",
      " [0.32230071]\n",
      " [0.16198059]\n",
      " [0.09532722]]\n"
     ]
    }
   ],
   "source": [
    "ITERATIONS = 20\n",
    "\n",
    "p_0 = np.zeros((n, 1))\n",
    "p_0[0, 0] = 1.0\n",
    "\n",
    "M = 1 / (6 * n) + 5 / 6 * A.T \n",
    "\n",
    "p = p_0\n",
    "prob = p  # 'prob' will contain each computed 'p' as a new column\n",
    "\n",
    "for _ in range(ITERATIONS):\n",
    "    p = M @ p\n",
    "    prob = np.append(prob, p, axis=1)\n",
    "  \n",
    "print(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOydZ3hc1bWw3zVNZdRlyVaxLRfAHWNMS8AGkpgSOrlcTMcQcAJJ+O4NgTQSSEKAlEsSOgFTQiAFQjWGEIoppthgMAYMxpZkWVaxitU1bX8/9pnRSBpNkTSW8ez3efScus/Z54y911l7NVFKYTAYDIbUxTbWHTAYDAbD2GIEgcFgMKQ4RhAYDAZDimMEgcFgMKQ4RhAYDAZDimMEgcFgMKQ4RhAYdhsicqSI1Ix1P4aDiOwnIu+JSLuIfHes+zMQEVEiMn2YbStF5KtDHDtCRDZFOldEfiQif45y3bNF5Pnh9MmwezGCIMWx/mN3i0iHiNSJyH0ikjXW/RopInKyiKwXkTYR2Ski/xGRihFc8gfAy0qpbKXUH0ehfxeIiN96721WX08Y6XVHG6XUq0qp/YY4dr1S6mIAEamwhJEj7PhDSqklu6uvhuFjBIEB4ESlVBYwHzgA+OEY92dEWF/GDwD/C+QCU4DbgMAwrhUc2CYDG4fZH8cQh9ZY7z0PuAf4u4gUJNDeYBgVjCAwhFBK1QHPoQUCACLydWtKpE1EtonIz8OOBb8CzxeRauvL+8dhxzMsDaNFRD4CDgq/n4jMFJGXRaRVRDaKyElhx+4TkdtE5Fnrq/l1EZkgIjdb1/tERA4Y4lHmA1uVUv9Rmnal1KNKqeqwa/8y7F79pqwsLekqEfkA6BSRF4GjgFusvuwrImki8lvruetF5A4RyQi/nnWNOmBFjPceAO4FMoCpQ7UXkW+KyGYRaRaRJ0WkdMCljheRLdbv8BsRsVntponIiyLSZB17SETyBrQ9SEQ+st7tChFJj/RuBvx+PxeRv1ibq61lq/WODrO0ntfCzp8hIv+2+r9JRM4IO3a8df92EdkuIt+P9s4Mo4sRBIYQIlIOHAdsDtvdCZyH/mr9OvAtETllQNPDgf2ArwDXiMhMa//PgGnW3zHA+WH3cgJPAc8DxcB3gIdEJHwa4gzgJ8A4oBdYA7xrbf8T+P0Qj/IuMENE/k9EjhrmVNdS63nzlFJHA68ClyulspRSnwI3Avuihc50oAy4Jqz9BKAArUlcEu1G1hf/xUAH8Fmk9iJyNPBr9DspAaqARwZc6lRgIbAAOBlYFryF1bYUmAlMBH4+oO3Z6N9omvVcP4nW5wgsspZ51jtaM+AZ3cC/gb+if++lwG0iMts65R7gUqVUNjAHeDHB+xtGgBEEBoDHRaQd2AY0oAdwAJRSLyulNiilAkqpD4CHgcUD2l+rlOpWSr0PvA/sb+0/A/iVUqpZKbUNCJ9bPxTIAm5QSnmUUi8CT6MHiCD/UkqtU0r1AP8CepRSDyil/MDf0NNYg1BKbQGORA/Ofwd2DsP28Uel1DalVPfAAyIiwDeB/2c9WztwPXBm2GkB4GdKqd5I17A4VERagTrruU9VSu0aov3ZwL1KqXeVUr3o6bvDBtg9brT6Uw3cbF0TpdRmpdS/rWs1ogXowN/wFut5m4Ff0f93GA1OACqVUiuUUj6l1LvAo8A3rONeYJaI5CilWqzjht2EEQQGgFOsL7EjgRnoL24AROQQEXlJRBpFZBewPPy4RV3Yehd6gAf9Bbot7FhV2HopsM2aFgk/Xha2XR+23h1he8iBXSn1plLqDKVUEXAE+ov1x0OdH4FtUY4VAZnAOmtaqxVYZe0P0mgJsGi8qZTKU0qNU0odqpR6IUr7UsLen1KqA2ii//sa+K5LAUSkWEQesaZc2oC/MPg3jNh2FJkMHBJ8X9Y7Oxut+QCcDhwPVInIKyJy2Cjf3xAFIwgMIZRSrwD3Ab8N2/1X4ElgolIqF7gDPdUQDzvQ0xBBJoWt1wITg/PYYce3J9jtmCil3gEeQ085gJ7uygw7ZcKgRhAtLe9OtCCabQ3keUqpXMvwG0/7eBjYvhY9mAKhqZZC+r+vge+61lr/tXW9eUqpHOAcBv+GQ7Udbn8Hsg14Jex9BaeQvgX6N1JKnYyeNnocrckZdhNGEBgGcjPwNREJGoyzgWalVI+IHAyclcC1/g78UETyLfvDd8KOvYUekH8gIk4RORI4kcHz3gkjIodbhtVia3sGcBLwpnXKerRhtUBEJgBXJHJ9S4u5G/i/sHuUicgxI+17FP4KXCgi80UkDT0V9ZZSqjLsnCutdz0R+B56+gz0b9iBNuSWAVdGuP5lIlIu2mvpR2Ft46URPZ01dYjjTwP7isi51u/tFJGDRDsMuETHHOQqpbxAG+BP8P6GEWAEgaEf1hzyA8BPrV3fBq6zbAjXkNiX2rXoaYataKPwg2H38aAH5+PQX9i3AecppT4Z6TMArda1N4hIB3ra5l/ATdbxB9G2jEqrX4kOegBXoY3qb1rTLS+gDeZJQSn1H/Rv8iha05pGf5sEwBPAOrSgewZtgAX9OywAdln7H4twi7+i38UW6++XEc6J1r8utG3hdWvq59ABx9uBJVafa9HTiTcCadYp5wKV1rtcjtZaDLsJMYVpDAaDIbUxGoHBYDCkOEYQGAwGQ4pjBIHBYDCkOEYQGAwGQ4rzhUtmNW7cOFVRUTHW3TAYDIYvFOvWrdtpBVgO4gsnCCoqKli7du1Yd8NgMBi+UIhI1VDHzNSQwWAwpDhGEBgMBkOKYwSBwWAwpDhfOBuBwWAwBPF6vdTU1NDTEyvRa+qQnp5OeXk5Tqcz7jZGEBgMhi8sNTU1ZGdnU1FRgS4TkdoopWhqaqKmpoYpU6bE3S5pU0Micq+INIjIh0McFxH5o1V67wMRWZCsvhgMhr2Tnp4eCgsLjRCwEBEKCwsT1pCSaSO4Dzg2yvHjgH2sv0uA25PYF4PBsJdihEB/hvM+kiYIlFKrgeYop5wMPGAVF38TyBORkmT155O6Nn773CZaOj3JuoXBYDB8IRlLr6Ey+pfHq6F/2b0QInKJiKwVkbWNjY3Dulnlzi5ueWkztbuGKh9rMBgMu5fKykrmzJkT+8QwVq1axX777cf06dO54YYbRqUfYykIIukvEYsjKKXuUkotVEotLCqKGCEdkwK3C4CWTu+w2hsMBsNY4/f7ueyyy3j22Wf56KOPePjhh/noo49GfN2xFAQ19K+TWk7idVLjpsCtXamaOnuTdQuDwZBiVFZWMnPmTL75zW8ye/ZslixZQne3nnVYv349hx56KPPmzePUU0+lpaUFgHXr1rH//vtz2GGHceutt4au5ff7ufLKKznooIOYN28ed95556D7vf3220yfPp2pU6ficrk488wzeeKJJ0b8HGPpPvokcLmIPAIcAuxSSu1I1s3yM4MagbERGAx7I9c+tZGPattG9ZqzSnP42Ymzo57z2Wef8fDDD3P33Xdzxhln8Oijj3LOOedw3nnn8ac//YnFixdzzTXXcO2113LzzTdz4YUXhvZfeWVf+eh77rmH3Nxc3nnnHXp7e/nyl7/MkiVL+rmBbt++nYkT+76fy8vLeeutt0b8nMl0H30YWAPsJyI1InKRiCwXkeXWKSvRtVE3owuBfztZfQHIy3QhAs1dZmrIYDCMHlOmTGH+/PkAHHjggVRWVrJr1y5aW1tZvHgxAOeffz6rV68etP/cc88NXef555/ngQceYP78+RxyyCE0NTXx2Wef9btXpNLCo+E1lTSNQCm1NMZxBVyWrPsPxG4T8jKcNJupIYNhryTWl3uySEtLC63b7fbQ1FAklFJDDtxKKf70pz9xzDHHDNm+vLycbdv6fGxqamooLS0dRq/7k1K5hvLdLmMsNhgMSSc3N5f8/HxeffVVAB588EEWL15MXl4eubm5vPbaawA89NBDoTbHHHMMt99+O16vHqM+/fRTOjs7+133oIMO4rPPPmPr1q14PB4eeeQRTjrppBH3N6VSTBRkumg2NgKDwbAbuP/++1m+fDldXV1MnTqVFStWALBixQqWLVtGZmZmv6//iy++mMrKShYsWIBSiqKiIh5//PF+13Q4HNxyyy0cc8wx+P1+li1bxuzZI9eEJNKc057MwoUL1XAL01zywFqqm7tYdcWiUe6VwWAYCz7++GNmzpw51t3Y44j0XkRknVJqYaTzU2pqqMDtosloBAaDwdCPlBIE2kbgiWh5NxgMhlQlpQRBoduFL6Bo7/WNdVcMBoNhjyGlBIEJKjMYDIbBpJQgCOYbMnYCg8Fg6COlBEG+22gEBoPBMJCUEgSFliAwsQQGg2FPYDhpqJctW0ZxcXHC7aKRUoIg3wgCg8HwBeeCCy5g1apVo3rNlBIEbpcdl91Gc5cRBAaDYeTs7jTUAIsWLaKgoGBUnyOlUkyICAVWLIHBYNjLePZqqNswutecMBeOi14FbHemoU4WKaURgJ4eajaJ5wwGwyixO9NQJ4uU0ghAVyozqagNhr2QGF/uyWJ3pqFOFqmnEWS6aDHFaQwGQxJJVhrqZJFyGkGh26SiNhgMyScZaagBli5dyssvv8zOnTspLy/n2muv5aKLLhpRX1MqDTXAzS98ys0vfMbmXx2Hw55yCpHBsFdh0lBHxqShjkEwzYSZHjIYDAZNygmCUOI5E0tgMBgMQAoKApNmwmAwGPqTcoJgWGkmenbBb/aBLS8np1MGg8EwhqScICgYjiBoqYTOBti+LjmdMhgMhjEk5QTBsIrTdDTqZXtdEnpkMBgMY0vKCQKXw0Z2miOxxHOdDXrZviM5nTIYDClJommot23bxlFHHcXMmTOZPXs2f/jDH0alHykXUAbBfEOJCAKjERgMhrHH4XDwu9/9jgULFtDe3s6BBx7I1772NWbNmjWi66acRgDDEAQdQY3ACAKDwdDH7k5DXVJSwoIFCwDIzs5m5syZbN++fcTPkZIaQaHbRUN7T/wNQhrBDggEwJaS8tNg2KO58e0b+aT5k1G95oyCGVx18FVRzxmrNNSVlZW89957HHLIISN+zpQc0fIzXbQkkoo6qBEEfNDVlJxOGQyGLyRjkYa6o6OD008/nZtvvpmcnJwRP0NKagQFbidNiaSi7mwEsYPya60gqyh5nTMYDMMi1pd7stjdaai9Xi+nn346Z599NqeddtrwOj2A1NQI3C56vAG6Pf74GnQ0QNEMvW7sBAaDIQbJSkOtlOKiiy5i5syZ/M///M+o9TepgkBEjhWRTSKyWUSujnA8V0SeEpH3RWSjiFyYzP4ECaWZiMeFNOCHrp26ZB0YF1KDwRAX999/P1deeSXz5s1j/fr1XHPNNYBOQ33ZZZdx2GGHkZGRETr/4osvZtasWSxYsIA5c+Zw6aWX4vP5+l3z9ddf58EHH+TFF19k/vz5zJ8/n5UrV464r0lLQy0iduBT4GtADfAOsFQp9VHYOT8CcpVSV4lIEbAJmKCUGnKEHmkaaoDnN9ZxyYPreOryw5lbnhv95I5G+O10WPIreP7HcOQP4chBMs1gMIwBJg11ZPakNNQHA5uVUlusgf0R4OQB5yggW/SkWRbQDPhIMgWJaATBYLKcUnAXGY3AYDDsdSRTEJQB28K2a6x94dwCzARqgQ3A95RSgYEXEpFLRGStiKxtbGwcccdCNQniiSUIuo5mFUP2BGMjMBgMex3JFASRTOMD56GOAdYDpcB84BYRGeQLpZS6Sym1UCm1sKho5B47CSWeC+YZchdDdonRCAwGw15HMgVBDTAxbLsc/eUfzoXAY0qzGdgKzEhinwDISXdikzgFQXBqKKvIaAQGg2GvJJmC4B1gHxGZIiIu4EzgyQHnVANfARCR8cB+wJYk9gkAm03Iz3TFZyPoaAC7C9LztEbQ0QB+U+bSYDDsPSQtoEwp5RORy4HnADtwr1Jqo4gst47fAfwCuE9ENqCnkq5SSu1MVp/CKXC74rcRuItARAsClBYGuQPNHQaDwfDFJKmRxUqplcDKAfvuCFuvBZYksw9DEXfiuY4GLQjAEgTo6SEjCAwGwwiprKzkhBNO4MMPP4zr/J6eHhYtWkRvby8+n49vfOMbXHvttSPuR0qmmAAoyHTxeWNH7BM7GyBrvF7PnqCXxmBsMBjGgLS0NF588UWysrLwer0cfvjhHHfccRx66KEjum5KppgAKMhy0RKXjaAxgkZgBIHBYNj9aahFhKysLEDnHPJ6vUPmLkqElNYIWrq8BAIKm22IF6lUn40AwD1OJ58znkMGwx5H3fXX0/vx6KahTps5gwk/+lHUc3Z3Gmq/38+BBx7I5s2bueyyy0wa6pGQ73bhDyjaeqJ4APW0QsCrg8kAbHbLhTQxjSDQ2YnyJFAIx2AwfGHY3Wmo7XY769evp6amhrfffjtu+0I0UlcjcDsBHUuQZxW0H0R4MFmQYQiCyjOXkrV4EcXf//5wumowGOIg1pd7stjdaaiD5OXlceSRR7Jq1aqE6h5HImU1ggK3/vGi2gnCg8mCZJckNDWkfD56P/+cnlFWWQ0Gw55LstJQNzY20traCkB3dzcvvPACM2aMPAY3dTWCzGCaiShTQ8HKZAM1gqrX476Pb+dOCATw1g4MqjYYDHsz999/P8uXL6erq4upU6eyYsUKQKehXrZsGZmZmf2+/i+++GIqKytZsGABSimKiop4/PHH+11zx44dnH/++fj9fgKBAGeccQYnnHDCiPuasoIgPzQ1FKVSWXjCuSDZE6C7Bbw94EyPeR9ffT0A3traqGqhwWD44lFRUdFvjv77YdO/8+fP58033xzU5sADD+T9998Pbf/85z8HwGazcf3113P99dcPeb958+bx3nvvjULP+5PCU0NxagRih4yCvn3ZpXoZp53AW6cFgertxd9k6h0bDIY9j5QVBJkuB+lOW2wbgXsc2MJeUyioLD47QVAjAMz0kMFg2CNJWUEA2k4QNc1EeDBZkASDyrz1fQLDu317ol00GAyGpJPSgiBmvqHOhgiCIFGNoAF7fj5gNAKDwbBnktKCoCCmIGjsbygGyMgHe1rcGoGvrg7XtKnYcnONRmAwGPZIUl4QDGkjUCry1JAI5MRfqczb0ICzeDzO0lK8241GYDAY9jxSWhDkZ7po7hhCEHg6wNc9WCOAuIPKlFL46upwTJiAs6wUb63RCAwGQx+VlZXDigr2+/0ccMABoxJDACkuCArcLtp7fXh8gcEHIwWTBYkzzYS/tRXl8eAcXxzSCJQaWLbZYDAYEuMPf/gDM2fOHLXrpbwgAGiNND0UCiYrGnwsTo3A16CFiWP8eFxlZQS6uvBb4eEGg+GLz+5OQw1QU1PDM888w8UXXzxqz5GykcUQFlTW5aE4Z0CUcCyNwNMBve2Qlj3k9X11Wlg4xo8Hux3QnkMOy4vIYDCMHq/+/VN2bouj2FQCjJuYxRFn7Bv1nN2dhvqKK67gpptuor29fdSeM6U1gvxgvqFIdoJQwrlIgsCKLm6LPj3ktYLJnJZGAMaF1GDY29idaaiffvppiouLOfDAA0f1GVJaIyjM6tMIBhFMQZ1ZOPhYeMnKoqG/Fnx19SCCo6gIW2YmYILKDIZkEevLPVnszjTUr7/+Ok8++SQrV66kp6eHtrY2zjnnHP7yl78M/wEwGgEALZFiCTobdI4hu3PwsfAi9lHwNtRjH1eIOJ3YcnOxZWYajcBgSAGSlYb617/+NTU1NVRWVvLII49w9NFHj1gIQIprBHmZwQykERLPRQomC5JtFbOP4Tnkq6vHOV5rDyKCs6zMxBIYDClCMtJQJ4uUFgROu42cdEfkVNSRgsmCpGWDKzumRuCrr8c5aVLf/UpLjUZgMOxF7O401OEceeSRHHnkkcPr+ABSemoIoDArjeauSBpBw9AaAVjRxdEHdW99Pc7x40PbWiMwNgKDwbBnkfKCID/TGdlG0NEY2XU0SPaEqBpBoKuLQFubdh21cJaVEmhrwz+Kbl8Gg8EwUlJeEBS4XTQNFATebvC0Rw4mC5IdPd9QyHV0Qn+NAIwLqcFg2LNIeUGQn+karBFECyYLEtQIhkgZ4avviyoO4izV8QfGYGwwGPYkUl4QFGS5aO7y9M8BFKlW8UCyS8Dv0fWLI+CrD4sqtggJAqMRGAyGPQgjCDJdeHwBujz+vp0hjWDc0A2DsQRtkQd1r6URhBuL7YWFSFqaMRgbDIY9irgEgYg8KiJfF5G9TnDkh4rYh00PdcYzNRQ9qMxXX48tJycUUQxWLIFxITUYDBbDSUNdUVHB3LlzmT9/PgsXLhyVfsQ7sN8OnAV8JiI3iMiMeBqJyLEisklENovI1UOcc6SIrBeRjSLySpz9GTUKIwoCa2poqDgC6J9mIgLe+jqc4wcLEuNCajAYRspLL73E+vXrWbt27ahcLy5BoJR6QSl1NrAAqAT+LSJviMiFIhIhBwOIiB24FTgOmAUsFZFZA87JA24DTlJKzQb+a9hPMkzy3RHyDXU0QlouONOHaEXM2sW++gYcxeMH7TcagcGw9zAWaaiTQdyRxSJSCJwDnAu8BzwEHA6cDxwZocnBwGal1Bar/SPAycBHYeecBTymlKoGUEo1JP4II6MgUr6hzoborqMAjjSdi2gIjcBXV0favvsM2u8sK8Pf3Eygq6vftJHBYBgZL913Fw1VW0b1msWTp3LUBZdEPWd3p6EWEZYsWYKIcOmll3LJJdH7Fw/x2ggeA14FMoETlVInKaX+ppT6DpA1RLMyYFvYdo21L5x9gXwReVlE1onIeUPc/xIRWSsiaxsbG+PpctxEtBHECiYLklMaUSNQXi++nTv7GYqDhDyHdsRX89hgMOzZ7M401KAzkL777rs8++yz3HrrraxevXrEzxCvRvBnpdTK8B0ikqaU6lVKDWWtiJRrdaDTvQM4EPgKkAGsEZE3lVKf9muk1F3AXQALFy4c1VqPOekOHDYZbCwujqMMXPaEiGkmfDt3glI4rIRz4TjL+lxI06ZNG3a/DQZDf2J9uSeL3ZmGGqDU+pgsLi7m1FNP5e2332bRokXD6Hkf8RqLfxlh35oYbWqAiWHb5cDAUbMGWKWU6lRK7QRWA/vH2adRQUTId7to6WcjaIhPIxgizYTPiip2DGEsBlOXwGDYm0lWGurOzs5QZbLOzk6ef/75hL2OIhFVIxCRCejpnAwROYC+r/wc9DRRNN4B9hGRKcB24Ey0TSCcJ4BbRMQBuIBDgP9L6AlGgYJMF03BKmU+D/S0Rg8mC5JdAh31EPCDzR7a7a0LppcYrBE4iorA6TTRxQbDXk4y0lDX19dz6qmnAuDz+TjrrLM49thjR9zXWFNDxwAXoL/mfx+2vx34UbSGSimfiFwOPAfYgXuVUhtFZLl1/A6l1Mcisgr4AAigp6A+HPqqyaEgXCMIuY5GCSYLkj0BVEC3ye4b9H0NQY1gsI1AbDacJSVGIzAY9gJ2dxrqqVOn9ms7WkQVBEqp+4H7ReR0pdSjiV7csiusHLDvjgHbvwF+k+i1R5MCt4tP6tr0RjzBZEGCtYvbd/QTBN66esTlwp6XF7GZcSE1GAx7ErGmhs5RSv0FqBCR/xl4XCn1+wjNvnDku520BGsSdO7Uy7imhqzBv20HlB4Q2u2rr8cxfvyQRiFnWSmdq18dSZcNBoNh1Ig1NeS2lkO5iO4VFGTqqSF/QGEP5RmKEUcAYWkm+ruC6qjiwdNCQZylpfgaGwl4PNhcruF222AwEN0TJxVRQ2REjkasqaE7reW1w+zTF4ICtwulYFe3l4Lg1FA8GoG7CMQ2yHPIV99Axrx5QzYLeg75amtxVVQMt9sGQ8qTnp5OU1MThYWFRhighUBTUxPp6VGyIkQg1tTQH2Pc9LsJ3W0PJTyorKCjEZxucLljtALsDm1LCNMIlFJ6amhCdI0AdCyBEQQGw/ApLy+npqaG0Q40/SKTnp5OeXl5Qm1iTQ2tG353vjgUWIKgpcsTX3qJcHJK+mkE/tZWlMcTY2rIVCozGEYDp9M5KAWDIXHi8Rra68m38g01dXjiDyYLkl0CrdWhzb5gssExBEGcE8aD3Y7HuJAaDIY9gFhTQzcrpa4QkacYnB4CpdRJSevZbqQwK1wjaISCqfE3zp4A294KbXrrtHYQKQV1EHE4cIwvxmc0AoPBsAcQa2roQWv522R3ZCwJagTNnZZGMPHg+Btnl0BXE/h6wZEWsVZxJFylZUYjMBgMewSxpobWWctXRMQFzEBrBpuUUp5obb9IpDvtZLrstHR060E9oakhawqoox7yJulaxTYbjnHRI5OdZaV0vvPOCHptMBgMo0O8aai/DnwO/BG4BdgsIscls2O7mwK3C09bI6Dicx0NEoou1lNC3vp6HIWFiDNivZ4QzrIyfHX1KCvBlMFgMIwV8aah/h1wlFJqM4CITAOeAZ5NVsd2NwVuF4GOGr0RTzBZkFB0sZ7v99XV44iQbG4gztJSCATw1jfgKh9YpsFgMBh2H/GmoW4ICgGLLcBuryaWTPIzXdiCCecS0gj6F7H3NdRHTD89kFAsgbETGAyGMSaW19Bp1upGEVkJ/B1tI/gvdJrpvYYCtwvHjmDm0QQEQWYB2JyhoDJvXT2ZBx8Ss1moLoHxHDIYDGNMrKmhE8PW64HF1nojkJ+UHo0RBW4Xab1NWkdKJKBMRGsF7XUEOjsJtLfH9BgCcJRoTcJbazQCg8EwtsTyGrpwd3VkrClwu7AHWlHONCQtJ7HGOSXQvgOv5TrqjJJeIojN5cJRXGwK1BgMhjEnLmOxiKQDFwGzgVA2I6XUsiT1a7eTn+kiTXbhzxyHI9HkVdkToOHjvoI0xbEFAZi6BAaDYc8gXmPxg8AEdMWyV9AVy9qT1amxoMDtZBy78KYXJt7YmhoKRRXHoRGAthMYY7HBYBhr4hUE05VSPwU6rfxDXwfmJq9bu58CdxqF0ka3cziCYAL0tuGr3QbEjioO4iwtxVtXh/L7E7+nwWAwjBLxCoJg1FOriMwBcoGKpPRojChwOxknu2h3DsMGbrmQ+mqqsOXmYsvIiKuZs6wUvF58JoWuwWAYQ+IVBHeJSD7wU+BJ4CPgxqT1agzIz3BQSBu7JHKd4bHsVMkAACAASURBVKhYQWXeHdtxFsfvempcSA0Gw55AXMZipdSfrdVXgARSc35xyJMu7OKnaViCQAeH+RoacZTG/3r6gspqYcGCxO9rMBgMo0C8uYYKReRPIvKuiKwTkZtFZBiT6Xsu9i49PdMYyE28cVAj2NkSV1RxEBNdbDAY9gTinRp6BJ1S4nTgG8BO4G/J6tSYYNUqrvVnJ942LRtld+Nv68IZp+sogC0jA3tBgZkaMhgMY0q8SecKlFK/CNv+pYickowOjRkdWhDU9GYl3lYEn60YVG/UWsWRMC6kBoNhrIlXI3hJRM4UEZv1dwY6++jeg5Vwrqo3jqL1EfAqPVMWrVZxJExQmcFgGGuiCgIRaReRNuBS4K+Ax/p7BPh/ye/ebqSjAT92qrtcw2ru8+kppXhSUIfjLCvDW1uLUoMqgRoMBsNuIVauoWFMmH9B6Wyky5lPc5cPpRSSYJoJX28aAI6iBDKXojUC1duLv6kpclWz3g7we3SWU4PBYEgC8U4NISInichvrb8TktmpMaGzkd60Qrx+RXuvL+Hm3i4bYlfY0xP7so/pOfTkd+CuxeDpSrhPBoPBEA/xuo/eAHwPHUj2EfA9a9/eQ0cDvgz9Rd7SmXg5Zl+HH0eGH+moT6hd1KCyQAC2vASt1fDGHxPuk8FgMMRDvBrB8cDXlFL3KqXuBY619u09dDairII0zcMQBN7WHpwZ/lCBmnhxllkaQSRBsHMTdLdAZiG8djO0bku4XwaDwRCLuKeGgPCQ22FEXe3BKAUdDdiyhy8IfC3tODL9oZKV8WLPysKWmxt5aqjqDb38xgpAwb9/mnC/DAaDIRbxCoJfA++JyH0icj+wDrg+ViMROVZENonIZhG5Osp5B4mIX0S+EWd/RpfeNvD3kparPX4SFQRKKXyNTVojaEvcFdRZWhq5QE31GsgaD1MWwZevgI3/gsrXE77+sOlq3n33MhgMY0ZMQSDafeY14FDgMevvMKXUIzHa2YFbgeOAWcBSEZk1xHk3As8l3PvRokPHEKTna0HQ0pWYIPC3tKC8Xhw5aQlrBKCnhyKWrKxaA5MO0+Uwv/w9yCmHZ6+CwG5IW/3+I3DTVHjz9uTfy2AwjCkxBYHSDu6PK6V2KKWeVEo9oZSKZ7Q7GNislNqilArGHpwc4bzvAI+iU1iMDVZ6ibTcCbjsNpo7vTEa9MdXb1UmK8xL2EYAfRpBv1iC1mpoq4HJX9LbrkxYch3Ub4B3H0j4HgnR9Dk887/gSINVP4RNzyb3fgaDYUyJd2roTRE5KMFrlwHh1s0aa18IESkDTgXuiHYhEblERNaKyNrGZOTut9JLSFYR+W4nzZ29CTUPVSYrHjc8jaC0lEBXF/7W1r6dVWv0ctJhfftmnwaTvgQv/kIbkZOB3wuPXgw2Oyx/DUrnwz8vgh3vJ+d+BoNhzIlXEByFFgafi8gHIrJBRD6I0SZSRNZAJ/ubgauUUlHnOpRSdymlFiqlFhYVFcXZ5QSw0kvgLqbAnTYMjUALEkdp+TCnhiK4kFa/AWk5MH523z4ROO4GPXf/cuxyEL6mJtpWrUosavml66H2XTjxjzBuH1j6CGTkw1//G3aZnEgGw95IvILgOHQdgqOBE4ETrGU0aoCJYdvlwECL6ELgERGpRGc1vW1Mktl1NgICmYUUuJ0J2wi89XVgs+EomQgdddr/PwFCQWXhgqBqDUw8RH+Zh1OyPxx4Prx9FzR8EvW6TXf/me1X/D/ann46vo5sfRVe+z9YcB7Mtn6G7Alw1t90hPPD/62XBoNhryJWrqF0EbkCuBIdO7BdKVUV/Itx7XeAfURkioi4gDPR1c1CKKWmKKUqlFIVwD+BbyulHh/uwwybjgbtq293kJ/pSthryFffgGPcOCSvDAI+6NqZUHtXUCMIupB2NukYgsmHRW5w9E/BlQXP/VC7vg5Bx+rVANRd94vYie26muGxS6BwGhw7IFZwwhz4rxVQvxEevWj3GKsNBsNuI5ZGcD/6q30DWiv4XbwXVkr5gMvR3kAfA39XSm0UkeUisnyY/U0OnY2QpWMICt3DEAR1dTrZnFWgJlGDsS03F1tmZt9gXR20D3wpcgP3ODjyavj8Rfh0VcRTPDU1eLZsIf/cc8Hvp/bqH6KG0lSUgqe+q9/D6feAK0IG1n2+BsfdpO/33I8Ser6E8PvgzTvg1d9Be2JR2gaDYXjEqkcwSyk1F0BE7gHeTuTiSqmVwMoB+yIahpVSFyRy7VGlowHc2vaQ73axq9uLzx/AYY9v5szbUE/alCmhIva01+kpnDgREasuQZggsKdBWZTylQd/E9at0F49047WHj7hj2RpA/lnLSV9xn7s+PFPaL7/AQovvGDwtd69Hz5+Cr72C20cjnbP5i3w5m1QMA0OuWTIU3u3bsVZWootLW3Icwaxa7s2VFdbgXQvXQ8zT4SDLobJX9Y2EoPBMOrEGulCVlPrC3/vpLMhpBEUuHUa6tbu+A3Gvrp6HOMnhAmCYbqQBjWCqjeg7MBBg3s/7E449tfQslUPzAPoXP0qzokTcVVUkHvaaWR99Ss0/v739Gz6tP+JjZ/Cs1fD1CPhsMtjd3TJL2G/42HVVfBp5NCP9hdfZMvxX6f6/Avwt7fHviboa91xuPZOOvUu+M67cMhy+PwluO/rcNuh8NZd0LMrvusNh84mqH4zFFdiMKQKsTSC/a16BKC9gDKsbUGHGOQktXe7i45GsPIM5WdqQdDc6WFcVuyvWX9HJ4GODl2rOKsYEGgbhiAoK6Xrvfe0MXbH+3D4FbEbTf8q7HscrP4t7L80NDUV6O2l8623yDv11FA67ZLrrmPLiSdR+4MfUPGPv2NzucDXC48u0zEKp94Jtjg0IJsdTrsbVhwH/1wGy1bBhLmhw90fbmT7/34f1+TJdH/4IdUXLmPSn+/GnpcX+Xo+D/znWlhzi77ON+6DcdP1sWN+BUf/BD58DN75Mzx7Jbzwc5j3X7DwIiiZF7u/0ehu1UK38lXYuhrqP+w75i7SHlvFs2H8LL1eNAOcGSO7ZzgBKyXJrhrYtU2/24wCnXI8uBzN+xkMQxCrHoE92vG9Ak8neDshS08NFbr7BEE8+Br0PLZzwgT9le4uGp5GUFZGoK0N/6bV2JV/aPvAQI75Fdx6CLxwLZyqo4C71q5FdXfjXnRE6DRHQQElv/wFNd/6Njv/+EeKv/99+M91ULdBu4hmJ1BQJy1LexLd/RXtVnrxfyCnBG9tLdu+tRxHfj6TH3yA7g0fsv1736PqwmVMuvceHPn5/a/TUqmFyfZ1cNA3tbbhTB/wYjLggLP13/Z3Ye098P7fYN19UH4wHHQRzDplcLtI9HboL/7K1Xrg3/E+qAA40mHiwVrojJ+jp7/qN+q/tfeAr0e3FxsUTB0sIPIqIgvR3g5rkLcG+l3b+rZbt0F7rXYuiIYjI0ww5A8WFMGlzaH/LQf/PXvC/zoGbA845veCw6Xfg91aOtLC/sL3p/ed60jT+wM+fY3Q0qttPQHvENth54uA2LUQDC7D1yPtszn0bwHavqUCMf4inIMCxLq/LfJ68DcXsfaFraMsRw3V14/wfeFOHJH2QdhUpwzYjnJ8xgkw74zo/2aGQbw1i/derGCycBsBxJ+KOhRVHCxanz1h2EFlAN73X8IuNj0wxUPhNDjs2/D6H/RcevmBdK5+FXG5cB9ySL9Ts486irwzzqDpnntxT8vC/f4tegDe77iE+0tOqRYG9x4LD/83/m/8g22XLkf19DLx3ntxFBWRffRRlN92KzWXf4fq885n0op7+4rvbHwcnvyuXj/jAZgVKeh8AGUL9N+SX8L6h/Ug/a9LtZ3kgHNg4TIomNJ3vrcbtr3d98W/fZ0efGxOKD8IFl2p8ziVLRxakAT80LxVawsNH2nhsOMD+OhJQoOA0w3FM2DcvtDTBruq9WA/MOhP7Pq95U6ESYdCbjnkTdTbOWV6gOpu1u26mvV614Dt+o1956g43JSdbm38d7m1p5nLDem5uh/B/XaXLn7k69FaYuivR+/vbR+839cL/l59XOz6I8jmBLvDWjr1gB1xv1NroTZr+An49HtWAX3dgA+Uv29fwB9hn4++wdnWN4gP+TfgONB/4A70H7D7rQcGr4eEhfQfsPvtG3jMWgbvHbboL1CibJfHOS4kiBEEnZarp7u/jaA5zlgCb11QIwgKgpJhawQA3k/Wkj5xDqQnMOu26EqdG+jZH8BF/6Zj9WoyDz4YW8bgaYXxV/2AzjWvU/vLm5l69gzsS36RcF9DlMyD/1qBeuhMtp97Ir3V3Uy6607S9tkndErWEUcw8Y7b2fatb1N13vlMuvt2nO/drAfxsoXwjXsgvyKx+2bka+F36Ldg6yvwzj2w5lZ4408w/SvavlL1hhYC/l49UJUeAF/6Lkw5QsdnRPKMioTNrqeqxk3vi60A/TXd8Ak0bOzTHra8rPuWW67/w+aW60E+b6Jez5qgB8TRIBCA3l19gsLv1ZqaM7NvwHdmxjfdNxKUMkb8vQAjCKw8Q8GpoZCNoCNBjSBYtD6nBGrfS7gbIY2gegscfmZo/+Z1DXS19TLvqIlDNYW0bPjKz+CJb+N54Q48W7eSv/TMiKfaMjMp+2o6lfdBffUCSkc4B632WUJd3Vfp3PIhJWcegPtLg6e03IcdxqS772LbJZdQddpxTF5Ui/Nr34WvXKO/EIeLiDZyTz0S2mrxv3Y3O1c8TG/De2TNKiTn6HNx7L9Ep+lIRLDGg8sN5Qfqv7HAZtNCJyM/9rnJxAiBvYIkfy58AQhNDWmNwOWwkZ3miFsj8DXUY8/NxZZuTS1kl2h/fH9iaSrshYWIy4m3LdAvv9A7z2xlzWOf4/PECOLafymUHUjHI3/Qj3PEEZHPe/tuMrpfZdxJh7Jr1cu0Pfd8Qv0cSNOf/0zrKx9S+NXp5PGMNupGINO5mYmLduLv8lP19mw8s5aPTAiEoZSi7Y0NbLnhRZo/tOFNm079Sx189rNVVP36EVqeeq5/HieDwdAPIwhCeYb6chjlu11x2wi8dfV92gBYRlcFCZasFBGcBW68XfZQxtGuNg/NtZ34vAFqNsVIMmezwbE30lnpwTnOjauiYvA59Rvh+Z/APscw7hd3kT5nDnXXXIO3YXiJX9uefZbG3/2enOOPp+jmx2CfY2DlD+CzF/pO6u2Afy2Hx5eTOX8ek+6+DX+3l6pzz8VTXT2s+4bj2baNbZdcyvYrrsA+rpCKv/+Nac8/z9Snn2Lc8uX4andQ99Nr+PTwI6i+9FJ2PfEE/o7RSZMR8Hjo/uADmh/8C9t/8AO2nHwK1cuWUXfddTQ/8AAdq1fjqapC+ZLnea2Uwt/Whq+pCeVJvKCSwQBmakhrBOl52hPCIt/toikBY7FjQrggCAsqyy1PqCvOTB/e9sxQTMP2T/sG/6oPm6iYOy5q+0DxPDobM8mbshNp2ao9XIJ4u3WwVkYenHIb4nJRetNNbD3tNHb8+CdMvOvOkKtpPHS99x61V11NxoIFlPz6esTh1PP99x4H/7gALnpOzx//4wJo2gyLr4bFPyDDZmfyfSuoXnYRVeecy6T77iNt6pRYtxv8rB4Pzffey87b70AcDsb/6Efkn7UUceh/0mnTp1P03e8w7juX0/vxx7StXMmulSupvepqxOUia/Ficr5+PFmLF0e0pQxEKYW3qoruDz6g+/0P6N6wgd6PP0Z5teZnLxpH+oyZ+Ftb2fXU0wTC4yccDlxWTIdr8mS9rKjANaUCR3FxxPeulMLf2oqvoRFfo/XX0NC3Hratevuy5Up6OvbsbGw5Obr6XU6OtZ2NPTsntLTnZGPLzsaenY2kp6N6egj09BDo7tbr3T2oXr0M9HSjuvVx1dPdb5/q7QWnA3E6EZcLcTq1a7K1FKez37F+S4dDJ0T0+1E+P8rn7Vv3+8DnR/l8fev+4H4fymdpyDZBbHaw2/TSZkPsNsu7KPyY3id2y1hss2k7iwroPgQUBAKooNE40GdAVtYxlLKOB182SMgQHfzTH3UhY7HNpu3DItbvbHkcWb+xNkCHfvRBnkhqwPHMQw4m+8gjY/57TRQjCMKCyYIUul00tPfE1dxbX0/6rJl9O4aZZoJAAKejhZ6uzNCu7ZtacKXbKdknj6oNTagzVdTBuuudtShvgKzyADz3E1j6176D/75Ge72c85hOUQGkTZ1C8Q+upP66X9Dy8MMUnHVWXF31VFdT8+3LcJRMoPzWW/qih9OytSfRn78CD5yig78y8uH8J7V3jkX6rFlMuv9+qpcto+q885i84t5+BuZYdL71NnXXXotnyxayjz2W8T+8Gme4VhaGiJA+axbps2ZR9L//S/f69bStfJa2Vc/S/u9/I5mZZB99NDnHH4/78C/rQQzwNTfT/cEH9HzwAd0fbKB7wwYCu3Qwm2RmkjF7NgXnn0f63Hlk7D8Px/jxod9GKYW/pQVPZSWerZV4qqr0emUlnW+80X/gzsjQwmHyZJTfFxrk/Y07Q0ImHFtWFo6iIhzFxWTMn6/Xi4oQl4tARzv+tnYC7W19y9ZWvNXV+NvbdXBfhGvGgzidSEYGtvR0JCMdW7q17nKhej0E2jtQXi/K4+m/DK4PV1txOBCHA7Hb9brdHlpHAH8AFfCDP6AH8kBAC5OggAno/fjjzI8l1uBts2n/Hms9NJBb66FBOzRYD/5TYAmVQNggH3af4DIkOPS2RDku6WlGECSFsGCyIPmZLjbVxY6IVR4P/qYmHVUcJFsbfRN2IW38BGdaN/5OJ4GuLmyZmdRsaqF033wq5hZStaGJ5tpOCsuyhrxE56urEZeLzNMugNeu17mIph0Nm1bpbKWHXa69asKfdelSOl58iYabfoP70MNifp37W1vZdulyCASYdOedg2MDcst0XMJ9J2gPnVPuCBniw0nfb18mP3A/1RdcqL2JVtxL+owZUe/ta2qi4abfsOuJJ3CWlzPxrjvJWrQoaptwRITMAw4g84ADGH/1VXS9s5a2lStpf+452p5+GltODhkHzMfz+Ra8NTW6kc1G2r77krNkCRn7zyN97jzSpk/Tg1GU+zgKCnAUFJC5oH+aEBUI4KuvDwkGT6UWEr2bNiFOB46iYtIqpuAoLgoN8o7i4tB6PNrLUCilUL29+NvaCLS362VHB4HubmwZGUhaGrbQYG8t0zOwpaeFNK2R3Bufb5CQ0F/wYYO8wxFax2ZLSEuN2YdAmLAIBLSWED7wp7Dh2wiCzoZ+kbEABW5nXAFlvsZGUEpHFQfJLNT+0YlqBNVv4HTrrxbvjh148svY1dDN3MXlTJ4zDthE1YdNUQVBx+pXtdvo4u/Bxr9q//pzHoUnvq2f8SvXDGojIpT86ldsPekkaq+6ioq/PoQ4IxtxAx4PNZd/B29NDZNW3BvZDgE6X9GVn+mgoyj/udKmTWPygw9QdcGFVJ1/AZPuuYeMObMHnacCAVr/+U8afvd7Al1dFC6/lHGXXjqiQVHsdtyHHoL70EOY8NOf0LlmDW3PrKR7wwbS58whf+lSPfDPmoUtMzP2BeO9r82Gs6QEZ0kJ7sOGyC6bJEQESU/Xjg3FxbEbjPK9saaJGMX3mVAfwr/0Df0wxuKOxn6GYtA2gm6vn+4YnjpeqyBNv2kJm037iyeaZqJqDc5xufq627eH7ANl++WTlZ/GuIlZVG4YOr21p6YGz9atZC06QgdHHXM9NH4Cdx0Jni44/d4hcxc5xxcz4brr6NmwgZ23Ry4Wp5Si7qc/pWvtWkquv57MhQujP48zIy7XQldFBZMffAC72031hRfS/X7/Smg9mzZRddbZ1F3zM9L33Zepj/+L4iuuGJEQGIg4nWQtWkTpjTcwbeUzlN/8fxRetIzMhQtHVQgYDHsqqS0IfL06KMc92EYAsYPKfPV6+qff1BBY0cUJCAKloHoNzpnaJ91bW8v2T1pIz3JSWKoDnyrmjqPu8130DFE9LZhtNOQ2ut/x2r++s1FXNSvaN2oXco5ZQu7JJ7PzzjvpXr9+0PGdt9zKrieepOh73yX3xBPif7Y4cE2cqIVBXh7Vyy6ia906Ap2d1N94E1tPOx1PVRUlN/yaSQ/cT9q0aaN6b4PBkOqCIOg6OmAOOxhUFsuF1GsFkznHD1CzE00z0VoNbdtxzFoMTieemlpqNrVQtm8+YtNf1ZPnFGp5sbEp8qO8shrnpEl90zUiOovn6ffAgvPj6sb4n/wY5/jxbL/qKgJdXX3de/xxdt56K7mnnkrh8uSUknCWlTH5Lw/iKCqi+puX8PnXT6B5xQryTjuNac+uJO+UU1J6DtdgSCapLQgGBJMFKYgz8Zyvrl7Puebm9j+QU5qYRmAVopEpX8JZUkJrTQsdLb2Uz+gzxBZX5JCR7aRyw2BBEMw2mnXEEf0Hy+zxMPcbcUd/2rOzKbnh13irt1F/402A9tDZ8dNryDz0UEqu/XlSB2Pn+PFMfvABXBMnYs/NZfJf/0rJL64bOnOpwWAYFVLbWBzSCIYpCBrqcYyP4AeePQF6WrXvfjwpHKregLRcKJ6Fs7SU2mYHZEL5fn2CwGYTJs0upPKDnQT8AWxhRXO63lmL6unR9oER4j74YAqWXUjzPffimlLBzttuxzVpEuV//APicsVsP1IcRUVM+ddjYX7XBoMh2RiNAAYZi+MVBN66epwD7QOQeIGa6jUwSReqd5aW0ujNx53rIre4vxCpmDuO3i4fdVvb+u0PuY0ePDqZCYu+9z3S9tuPhhtuRJxOJt55B/ac3Vd6QlLclc9g2N2ktiAIJZzrrxHkpDuxCbTENBYPSC8RJBRUFoedoHMn7Pw0lF/IUVpKU8YkyvbJHTQYTpxVgM0mVA2YHgq5jY6SJ43N5aL0NzeRuXAhE2+/DVd5YhHSBoPhi0VqC4KORnBlD5q+sdmE/MzoaSZUIIC3oaEv/XQ4iWgEwUL1wfxCuRPxunKYEMHNOy3DQck+uf3cSD3btlluo/EHVsVD+r77MvkvD5Ixb4RVwAwGwx5PaguCzoaIUa+gp4eieQ35W1rA6+0rSBNOeL6hWFRZhepLDwBgp78AgOL0toinT54zjubaTtqauoE+t9HRsA8YDIbUJLUFQUfDII+hIPluV1QbQagOQSSNID1XlxiMSyN4A8oXhoK96lpcpHc3kjZEQFrF3EKA0PRQ5+pX+7uNGgwGQ4KktiDo3BlKwDaQgkxXVBtBqDJZJBuBiLYTxIou7u3QZQ8t+0AgoKir7iG/9TM827dHbJI3PpOcogyqPmzq5zZqMBgMwyXFBcHgzKNBYmoEDcHKZEMUfc8uiT01VPO2rsM6WQuCndva6e32MY56fLW1EZuICBVzCqnZ1ELbmndGzW3UYDCkLqkrCPw+Xe91iKmhQreLli4vgYCKeNxbVwd2O45xhZGvH0+aiao1Oje6VZA6WHymOKdnSI0AtBup3xtg638+QNLSRs1t1GAwpCapKwi6dgJqSGNxvtuFP6Bo74lcXcpX36DzwA+VjjinVGsEKrIgAbTH0IS5oXq62ze1kF/iJqe0AO8QGgFA6T55ONLsbNvSM6puowaDITVJXUEwRHqJIAVunYp5qMRzvvq6/umnB5I9Abyd0DtEXQOfB2regUnabdTvC1C7eRfl++XjKC3FV1cfsTAJgN1po7winQbnZNyHm2khg8EwMlJXEAwRTBakwK29eJo7eyMe99Y34IzkOhokVizBjvXg6wnZBxoq2/D1+infLx9XWRkEAqE015EoDmynNz0fzwwzLWQwGEZG6gqCjsFF68MpyAymmYj8Ve6rq8MxYQhDMcQuWVn1hl5aHkM1m1pAoHTfPJylusqZN4qdIO9THT+wvTlyjQGDwWCIl9QVBDE0gnxraihSUJm/o4NAZ+fg9NPhxAoqq14DhdP7CtVvamFceRbpbifOsjKAIe0Egd5e/G+9Qp6zY1C6CYPBYEiUpAoCETlWRDaJyGYRuTrC8bNF5APr7w0R2T+Z/elHR4MO+nJFLv0YTDwXKc1EKJhsKNdRiK4RBAJQ/WZIG/B5/NRtaQtlG3WUaCHirY2sEXS9rd1GJ++XTd3WXXR3DLMwuMFgMJBEQSAiduBW4DhgFrBURGYNOG0rsFgpNQ/4BXBXsvoziM6delpoiCyXmS4H6U5bxKCyPkEQRSNwuXVq6UgaQePHOk21lV+obssu/L4AZZYgsLlcOIqK8G6PrBF0vLoaSUtj+pI5oKB6Y3O0JzUYDIaoJFMjOBjYrJTaopTyAI8AJ4efoJR6QynVYm2+Cey+NJdR8gwFKciMHFQWiiqOZiMAK7o4wmAewT4gNqF0n74CLM6ysiGnhjqtbKPjp48jI8cVtZaxwWAwxCKZgqAM2Ba2XWPtG4qLgGcjHRCRS0RkrYisbWxsHJ3edTQO6ToaJH+IxHOhqOLi6O2HLFlZvUbbEPIrAG0fGF+RjSu9r06Qs7Q0orHYU12Np7JSVyOzCZPnFLLto2b8/kD0vhgMBsMQJFMQRJpziRhdJSJHoQXBVZGOK6XuUkotVEotLCqK/hUfN/FoBO7Iqai99fXY8/KwpadHv0ekNBNK6YjiSYeBCJ4eH/WV7aFpoSDOsjK8dXUov7/f/o7VrwKQtVinna6YU6iL1Xy+K3pfEkQpRXNtJ2qIyGqDwbD3kExBUANMDNsuBwbNdYjIPODPwMlKqd3jAhMIWDaC6F/0Be7Iied8dUMUpBlITok2FodHF7dWQXttyD5Q+1krKqD6laUEcJaVgteLb4AG1PHqapyTJ+GaPBmAiTMLsNkHF6sZKRtfreXh697inzeto6Eqckpsg8Gwd5BMQfAOsI+ITBERF3Am8GT4CSIyCXgMOFcp9WkS+9Kf7mad7G0I19Eg+UPYCHz19ZHTTw8kuwQCXp3TKEiVVYjGsg9s39SC3WFjwtTcfk0juZAGenvpeuttso7oK0LjynBQuk/e2oTsPwAAIABJREFUqNoJWhu6eP2fnzFuYhbtzT3844a1vPzQJ/R0RI6pGE38/gCfvlPH5nUN+L1mustg2B0krXi9UsonIpcDzwF24F6l1EYRWW4dvwO4BigEbrPKMvqUUguT1acQQ9QqHkiB20V7jw+vP4AzrFi8t76e9NmzY98n5EJaC24rOV31G7peQbF2oKrZ1MKEaTk4XP1zFvULKluwAOhzGx2YbbRi7jhe+8dn7GrsJrdoZHmHAgHFf+77GJvdxte/PQ9nuoN3ntrKBy/XsPndBg47ZRozv1yKzTa6NYV9Xj8fv76Dd5+voqNZR3NnZDuZ+aUSZh1eNuLnMhgMQ5M0QQCglFoJrByw746w9YuBi5PZh4jECCYLEowlaOn0UJyj7QHK48Hf1BS/RgDaTjBhrl6vWgMTDwWbjZ4OLztrOjjkxCmDmjqDsQRhLqRBt9GB2UYnzynktX98RtWHO5l3VN9snFKKns4OMrKyY/fVYv2/q6nbsouvXjiLrHz9zIefsQ8zv1zC6kc+5eWHNvHRa7UsOnM/xk8ZeUF7b6+fja9u573nq+lq8zBhai6Ll+6H3W7jw1e3896/t/Huc9VMmlXA7EVlVMwtxGZPjiLb2+2jtb6LrPw0MnNcg2pGGwx7K0kVBHssofQS8QmC5q4+QeBt0G0jFqQZyMCgso5GaPoMDjgbgO2ftYCCsn3zBzW1ZWZiL+ifhbTzldVkHnLwICN13vhMcoszqNrQ1E8QrH5oBe+teoql1/2G8VOnx+zuzpoO3npyC9MOKGLfg/s/X2FZFqf8zwF89k49rz+6mX/etJZZXyrh0FOnkZHliv0uBtDb7WPDSzW8/59t9HR6KZ+Rz5KLZlO6b15oAJ44q4DO1l4+er2Wj16r5dk7NuDOdTHz8FJmfbmU7IIYxvoYeHv97Pi8le2bWqjZ1EpjVVvInONMt5NblEFuUSZ5xRnkFmeQW5xJXnEmGdnOEQsJr8dPd5uHrjYP3e0eutu9iA3SMpy4Mh2kZThwZThIy9TL0dbAlFL4vQG8Hj8Bn8LmEOwOG3aHDZtdjBBMMVJTEHQGBUHk6mRB8kP5hvrsBDEL0oSTFRQEludQsFC9lXF0+yctONLsFFdE/rIOdyH1VFfjqaoi/5xzIp5bMXccG16pwdPjw5XuoH7r56x7+nGUCrDyT7/lnBtuxpk29MDp9wZ4YcVHpLmdLD5rv4gDgYiw78ETqJg7jref2coHL9bw+XuNHHrKNGYdHt90UXeHhw9erOGDl2rwdPuYPLeQhcdVDLKRBHHnpXHQ16dw4LGTqdrYzMbV21m7spJ1KyuZPHcccxaVMXFWQVz39nsD1G3dRc2mFrZvaqF+axsBv8JmE8ZPzeHA4ysoKs+mo7WXXQ1dtDZ0s7Omna3rG/vVpRgsJLQgzinMwNvro7vdQ1eb11p66Gr30G0N+Hrbi6/XH6Wng3Gm2/sJh9B6hgNXpgOxCT5PAF+vH5/Hj9fjx+cJ4A1tB/B5+q9H9uHT2ByC3W4JhjAhYbfWbXYbdqdgs9tQAUXAr/TSWg/49bo+FtDr/rDj1jYCYhP9+wnYbFoIiU0Q61hw3Wbrvx8sPwylUApUQIUEuV63thWh9dA+LLdGCS71dfW23hn8LyAS3BZru7//h1Iq9C5VqFPBY4Qd6//CJdyxMnSvAT9E2A4BZh1RyoIlk4f+4YZJigqCBrA5IWPwl3g4IY0gXBDU6UE9alRxEIcLMsf1aQRVb4AjPVSovmZTC6XTc7E7Ik91OMvK6P1U29BDbqNDVCObPLeQ9/+zjZpPWqiYV8ALd99CRk4OR1+4nKdvvoFX/rKCr170rSG7+vYzW2na3sHx3/7/7Z15eBzVma/fr1epta+2LFm2bEte8L7FO2B2J2NIDCGM50IgDAM3hG1IhiRzc5OZe58Lk2Vu5iYDwxZDwm6DgYTdwYAB73iRd+NFkmWttiRr7+XcP6q61ZLVUpVsYUyf93nqqVNV56tzurv6/E6d5TuTSUzpu4bvSXSx4Npixs/L46Pn9/PBs2Zz0Q0lDC3qvUBvaexg27tllH5USaAzyOhpOcy4ciQ5hdaarRxOB0WTsymanE1TXRu711Wy++NKjuyoIyUrgQsWDmP8vGH4UrvyHgqGqDl6KlLwH/+8kaA/hAjkFKYw9dLh5JdkkDcmHbc3xroSGB3Yp+rbaaxto7GmlcaaNkMkyk9xaFtt30NsBRKT3SSmeEhM8TCkKA1fiofEVOOcL9XYElM8qJCioy1AZ2vA2LcF6AiHWwN0tIev+Wlu6KDzeEvkmlLgcjtweZ24PA7cHicujxO310lCsocUT/iaE7fHEbnm8hiFeigYIuhXBIMhgoEQoYAy98ZxMGBe84cIBZV5LkSgMxApoJ1uB26HIE6jYHfECjsdpo1RRoYL8G77cEEeMvqtUIpQCHMfVZhLdEFthI0CXfq8hhARCFRUAR4JGwW3Co9XUKrrvCIiGtAlJJGfvNs16RIbus5Ha4LqEVDdVCWSPAApGWf2FhyL+BSC5to+3UuEie4jCBN2DW2paQi6zyUo+wTyZ4LLQ0tjByerWhk3Ly+mqXvYMJrXrkUpddqw0Z4MG5OOO8HJ0dJ6mmo2UfX5AZb84H7Gzl3A8QPXsOUvqxk1bSajps86zfb454189vZRxs/Po2hy329J0WQNS+bqe6dxcHMNH688wKqHtjB+fh5zrxkdEZNTJ9r57O2j7P74OKGQomTWEKZfOYLMvCTL6fQkNTuROdeMZtY3iji8vY7SD4+xfvUhNr52mKKpOeSOSKHyYAOVBxrwtxs176z8JC5YOIyCsRkMK07H63NbTs/pdJBuNgtxQfcV6aJF4lR9O26v0yzoPSSmuElMdg9an0aYcGEmZ7n5SBM/xKcQWJhMBpBuFhbRrqgDVVVIYiKOVIsdpWE3E+1NULUTFv4jAMf2G541es4fiMadn4/q6MB/rJLWDRtJv/bamHGdLgeF4zM5/NkRdp14isKJUxg3/0IAFtxwE2U7t/H2I7/lpl/+Dl9alysLf0eQNSt2k5yRwIJri619pihEhOJZQxgxKYvNfznC9jXlHPqslplLRnKisoV966tAYNzcPKZfUUhajs92Gn195jEzchkzI5eTVS3sWlfJ3k+O8/nWGtKH+CiZPZSCsRnkl6T3+5Yz4DxEi8Q5Ilzr1WgGSnwKQXNNvx3FAG6ng9QEV7dJZf6aaty5udY701KGGgJQsRFUqGv+wN6TeH0usofHbhoJDyFtfHW1MWz0wkUx4wKMmJTNno9WQKiDS7733yN5dLndLLnrh/zpx/fw9iO/5Zof/Sxy7ZOXD9JY18Y1907Dkzjwx8GT4GLesjGMm5fHRy/s5+OVB3G6HVxwYT7TLis8447d/sgYmsSCa4uZc/UoOtuC3ZqINBpN38SnELTUwhAL8wCArGRvNzcTgarqvhek6UnqMOMN5Mg6Y6H64V0L1Q8rTu+zkzM8qazhpZXGsNFZpzfrRONwlBPy76Nw8lVkDuvu1il7+AgWLb+Z91c8yvZ332Tq5Uso211P6QfHmGK2lZ8NMvOSWHr3VKqPNJGalfiFF8gutxOXO3Z7v0ajOZ34W5hGKUMI+plMFibD5+7WRxCorrbWURwmZajxJrD7VRg6GbwpNNW10VTXTsG4vgtfd77xRhCoqup12Gg0gc5O1j33GC5vJjim9xpn2pV/w8gp0/ngj09QeeAQf316LxlDfcy5epT1z2MBEWFoUZqulWs05wnxJwTtDRDs7HcyWZjMpC43EyoUwl9Tg9vK0NEw4UllJw5F/AuF+wd6OprriTM5OdIXEe1Wojc2vvoSDVXHGb/wBmrK2mhtOt01hohwxR334PZ6eeWhB2ltbOPSmyfoGrRGE+fEnxBYnEwWJsPX5XgueOIEBALWHM6FSYkSjaj1BxJT3JZGzoSbh2INGwU4UVnBxtUvMW7+hUy9fL65WE3vTuiSMzKZdOlNtJ+qJCuvlNwRZz47WKPRnN/EnxBYnEwWJjPZcEWtlMJfHV6Qxo4QRA0PLZyLUopje09SMDbDUoezd8wYvMVjYg4bVUqx5on/xOXxctGNt5IzPAVfmocjMbyRtjR2sH+Tj6Ss6VTsXkP5rh3WP4tGo/lKEodCYM3PUJhMn4fOQIjWzmDUEpU2hCApB8QJWcWQnENDdSstjZ39NguFGfqz/0HhihUxr+9dt5ay0h0suOEmktIzoharqT9tsRqlFGv/tBd/Z5Bv/uhuMobm8cbvf0N7c7P1z6PRaL5yxJ8Q2G0aippdPCAhcDghZxwUXw4Ybqeh//6BMM6UFFxZWb1ea29uZu0fn2DomBImX3pF5PzISdl0tgc5frD7YjV7PjnOkZ31zL1mNENGZrHkzvtpbTjJe4//3piUdJb4fMsGnv6nu3jhFw/wyUvPUL5rB4HO0/ssBoPWpkZOVlUSCtlz4aDRxDPxN3y0pcYYxunLtBQ9K0oIvFXV4HTGLJhjcuu74DTuU7HvJMkZ3rPiVnnd80/R1tTEsp/8Cw5HV4dvwbgMHC7hyM66yIS1pro21r14gPyx6Uy+2FgaeuiYEuZdt5x1zz/NqOmzmLBo8Rnlp7Gmir+ueJRDWzaSOawAv8PB+lUv8OnK53C6XOQVj6NgwiSGT5hEXslY3B7vGaXX1nyK6kMHqf78ANWHD1J96CBNtcYbn8vtISO/gOyCQrKGjyB7eCFZBSNIy8lFHGe3/hPo7KT5RD0tjQ243G68viS8SUl4fUk4nLojXvPlJ/6EoLnG8P/jsPYHzYjyQJpTXY0rNxex++f2GJ3CKqQ4tq+BkZOyzti7Y+X+vWx/7y1mLFlK7sjuwz89CS7ySzI4urOeBdcWo0KKNU/tAYHFN47v5opg1tXLOLxtC2uefJj8cRNIy7UxIsok4Pez+fWX2fDyC4jDwaK/u4XpVy3F6XLR3tLMsb27Kd+9k4rdpWx4+QXWrzKEYeiYsQy/ICwM4/oUhvbm5khhb2wHaDQdAAKkDRnK0DFjmXLZEhJTUqk/Vk59+VHK95SyZ93aSDyX10tWfqEhDMNHmEJRSEpWzmm/iVKKjtYWmuvraD5Rz6kT9ebeOA6faz8VewU3l9drCEOiLyIOHl8SXp/POB8VVkoR6Owg0NmJv8PYG8fR54zN39lJwIzj7+yAUAin243D5cblcnWF3S4cLjdOlxuneb572IXD6SQUDBIKBgwHcd3CAeM4FDTP99hCQUQcOBwO01mcA3GYx+FNehxHXTfcY4QdxClUKNTjXMj07xOKXAcie8OnkPQIm47rwn6GupwNRTmYk0g60MNZnFKGv5/IoeE5rreXZunhRKjbMxTtb6h75Eg6kWBPD4A9ndqZFM+ee8YVtt6IPyFoqbXcPwBGHwEY/ob81VW4+1uwvg/qK1tob/GT38/8gf4IBYO899jvSM7IZN51y3uNM2JiFutePEBDTStHdtRReaCBxTeOJzWr+5uIw+FkyZ3/yFM/vJM3fvcbrv+f/8dWLfbojm2sefJhTh4/RsnX5nPRTX9PSlZXR3xCUjKjZ8xm9AxjIl1Ha0uUMOxkw8svsn7V86YwlDB8wiQKJkxCxEH1oQORgr+h+njknmm5QxhSNIbJl17FkKIx5I4a3eeaCx2tLdSVl1FfcZT68jLqKso4sn0ruz5YE4njSfSRVTCc1OxcWpsaIwV+oKPjtPv50tJJzsgiJSubYSXjSM7IIjkrm6T0DIJ+Px2tLVFbKx0tLXS2ttDR1kp7SzONNdV0tLbQ2dpKwB+7yUzEgcvjweX14vJ4cHu8uDxm2JtAYkpq5Jw4hGAgQNDvN/YBP6GAn4A/QKitrdv50+IFgzidLsTpNBzCOQ1xcDicOFzm3hm1RZ13ulyRAjsUChIK+CGkCIVCZqFu7GMdizhMp3AO07NouAB3RMKRcw4HgpgVmS7vbV2O48yC3XQQFxGNqHNECnUVEY1wAR3tPE567KNFJEK0iEQdR9KE048jHuvCt41dIewuKka47dSpmPHPhPgTguYay5PJoGcfQQ3eYvv+eMJE+gfOcBbv1jdfo7bsCEvv+wmexN593IycZAjBZ++Wse/TKoqmZDNubu+1/dScXC793h288btfs3H1S8xZ9p1+83DqRB1rn36C/Z9+RPrQPJb9+BeMnDqjXzuvL4lR02dFnN91tLZybN8uynftpGJPKRtWv8T6l1/olrcho8Yw8eLLGDK6mCFFo0lMsTfk1etLIn/sePLHju92vq35FPXlR7uJRPWhg/jSM8gZOYpR02eSnJlNcmYWKeY+KSMTl9u6w7r+CAb8plg0GwV/uND3enE4XXpdAM0XQvwJQUsNZI22HD01wYXLIYYQVFWRtGD+gJOu2HeStNzEM/K701RXyycvPsOo6bMYM3tuzHhpOT4yhvrY/VElCcluLlo+rs9CZfzCizn02WY+WfksI6ZMI2/M2F7jBQMBPnvrdT556VlUMMi8by9n1t8sw+UZ2Cxir8/HqGmzGDWtSxgq9+9BgNxRY/Cl9u7W+myQmJxCwfiJFIyfOGhp9IfT5caXmjaon1Oj6Y/4EgKlulxQW0REyEjy0HyigVBrq3X30z0IBUNU7j9J8ayB2Yd5f8V/oZRi8c2391tbHDEpm5NVZVy8fJwldw+XfO8Oju3bzRv/71f8t4f+A09C92akir27WPPEw9SVHaFo2kwW33w76XZmWVvA6/NRZOHNQqPRnD3iSwg6WyDQZksIwOgn8FfZWJmsF2rLmulsD1oeNtobBzdv4OCm9Sz82++Sltu/oMy4cgQFYzMYMdHaKKeEpGSu+v59vPgvP+H9FY9xxe13AdDa2MCHz6xg1wfvkZKdw9L7f8qYmXN0s4VG8xUhvoTA5mSyMBlJbtTh8FrFA+ssrth3Ahh4/4C/vZ2//uERsgoKmfH1ayzZJCS5LYtAmOETJjF76TI2vrqSoqnTaW1qYt3zT+Fv72D21dcy51vfwd2H8zuNRnP+EV9CYHMyWZisJC/OenN8uh0X1FEc23eSzGFJA/bI+emq5zhVV8v1P38Qp2twf7Z5317O0Z3beP3fHwSgcOJkFt9yB1n5wwc1XY1Gc26ILyGIvBHYaxrKSHLjaTBq9K4BDB8N+kMcP9jIhAXDbNsC1JYdYctfVjPx4su+kI5Np8vNkh/8kDVPPszEiy9j3LxFuhlIo/kKE19C0GwKgc03gkyfB9VUjzMjA4fX/mzY6iONBPyhAfUPqFCI9x77PR5fEouW32zbfqBkDsvnun/+X19YehqN5twRX0Jg0/NomMwkD7Q1IjkD7B/YexIRyC9JP+1aoLOTjtYW2lua6WgxJyG1NJvnWjhRUUbl/j1cccc9tsfPazQajRXiSwiaayAxE5z2JgRlJHlwtDcRHFVkKX6gM0DpzsPs33GA2iNldJRXolQrv37gdcTfjsPfhnS24/C3I6FAn/dSDicdBZN4ujYbVm63lW+NRvPV4qKxuSyZlNd/RJvElxC01NgeMQTGG4GrrYHO9O4jcMrKati5eS/HDhylueoYgcZapKMeR+AkEO390k27O5mWU178Ti9+ZyadyeGwsXU6E3ocG/uQw/yJDva+voBGo4kfRmT1v5jVQIgzIaizPYcAwNt6gurs4Ww80cpfbrsPZ1sTDn8DotqjYgkOZxpBTyqBjDwcmWkkFqSTXZJFZmEqTudAPV6G+o+i0WjighGpg1MexJcQNNdA3hRLUf0dHbz5yH9RtrmUjs46yPFD+yHcnYkE3D5aUlM5leSjPqWV6vQGqtJqCbp6+ZGOmptGo9GcIbdMvIWSGSVn/b7xJQQWPI+uW7WSnX9+n7a2OpRqAVx43DlkjxvNyGWXkDVU+4TRaDTnhswEa+uo2CV+hMDfDh1NvTYN7d+4gQ//8DynTtYSUg2A4HJkkzWymMvv/AdyC/REKo1G89VlUIVARK4Efgs4gceVUg/2uC7m9SVAK/BdpdTWQclMD/cSNRXlvPUfj3CyvJJAqA5QOCWDtOwSFn33ekpmfW1QsqHRaDRfNgZNCETECfweuAyoADaJyGtKqd1R0a4Cis3ta8DD5v7s01xLW8DBm6/s4thv/p5Ofy0QQCSJpKRCJn39YuYvu3ZQktZoNJovM4P5RjAbOKiUOgQgIs8DVwPRQnA18LQylu9ZLyLpIpKnlDp++u3OjJf+8znKD12KUrsADwmeXEbMnsQVt92GewCzhTUajearwmAKQT5QHnVcwem1/d7i5APdhEBEbgNuAygsLBxQZoaOGsbxw/XkjB7LlffeRUb2wJec1Gg0mq8SgykEvXkp67n8s5U4KKUeBR4FmDlzZi9LSPfPwlvvY+GtA7HUaDSarzYDneVkhQogerhNAVA5gDgajUajGUQGUwg2AcUiUiQiHuA7wGs94rwG3CgGc4DGwegf0Gg0Gk1sBq1pSCkVEJE7gbcxho8+qZTaJSK3m9cfAd7AGDp6EGP46BfnZ1mj0Wg0wCDPI1BKvYFR2EefeyQqrIDvD2YeNBqNRtM3g9k0pNFoNJrzAC0EGo1GE+doIdBoNJo4RwuBRqPRxDli9NeeP4hILQP38J8N1J1B8mdifz7ansu0z0fbc5n2+Wh7LtOOx888QinV+8pcSqm42YDN58r+fLQ9X/Otv6/zw/Z8zff5+pn72nTTkEaj0cQ5Wgg0Go0mzok3IXj0HNqfj7bnMu3z0fZcpn0+2p7LtOPxM8fkvOss1mg0Gs3ZJd7eCDQajUbTAy0EGo1GE+fEjRCIyJUisk9EDorIAzZtnxSRGhEptWk3XETeF5E9IrJLRO62aZ8gIhtFZLtp/ws79uY9nCLymYj82abdERHZKSLbRGSzTdt0EVkpInvNzz7Xhu1YM83w1iQi99iwv9f8rkpF5DkRSbBhe7dpt6u/NHt7JkQkU0TeFZED5j7Dpv11ZtohEZlp0/aX5ve9Q0ReEZF0G7b/atptE5F3RGSYVduoa/eLiBKRbBvp/lxEjkX91kvsfGbz/A/M//UuEfk3G2m/EJXuERHZZsN2qoisD/83RGS2DdspIvKp+d96XURSY9j2WnbYecZsMRhjUr9sG4Yb7M+BUYAH2A5MsGG/CJgOlNpMNw+YboZTgP020xUg2Qy7gQ3AHJt5uA94FvizTbsjQPYAv++ngFvNsAdIP4PfrQpjIoyV+PnAYSDRPH4R+K5F24lAKeDD8Mr7HlBs55kA/g14wAw/ADxk0348MBZYC8y0aXs54DLDD8VKO4ZtalT4LuARq7bm+eEY7uaPxnpmYqT7c+B+i79Pb/YXm7+T1zzOtZPvqOu/Bn5mI913gKvM8BJgrQ3bTcCFZvgW4F9j2PZadth5xuxs8fJGMBs4qJQ6pJTqBJ4HrrZqrJT6EDhhN1Gl1HGl1FYzfArYg1FYWbVXSqlm89BtbpZ790WkAPg68LjlTJ8hZg1nEfAEgFKqUynVMMDbXQJ8rpSyM5PcBSSKiAujULe64t14YL1SqlUpFQA+AL4ZK3KMZ+JqDBHE3F9jx14ptUcpta+/jMawfcfMN8B6jNX+rNo2RR0mEeMZ6+N/8O/Aj2LZ9WNriRj2dwAPKqU6zDg1dtMWEQG+DTxnw1YB4Zp8GjGesRi2Y4EPzfC7wLIYtrHKDsvPmB3iRQjygfKo4wpsFMhnAxEZCUzDqNXbsXOar601wLtKKTv2/xfjDxqyk6aJAt4RkS0icpsNu1FALfAHs0nqcRFJGkD6YKxq1+sftDeUUseAXwFlwHGMFe/esWheCiwSkSwR8WHU9Ib3Y9OTIcpcYc/c59q0P1vcArxpx0BE/reIlAPLgZ/ZsFsKHFNKbbeXxQh3ms1STw6gmaMEWCgiG0TkAxGZNYD0FwLVSqkDNmzuAX5pfl+/An5sw7YUWGqGr8PCM9aj7BiUZyxehEB6OfeFjZsVkWRgFXBPj9pXvyilgkqpqRg1vNkiMtFimt8AapRSW2xn2GC+Umo6cBXwfRFZZNHOhfE6/LBSahrQgvEKawsxljddCrxkwyYDo8ZUBAwDkkTk76zYKqX2YDSpvAu8hdF8GOjT6EuIiPwUI9/P2LFTSv1UKTXctLvTYlo+4KfYEI4ePAyMBqZiCPevbdq7gAxgDvBD4EWzhm+HG7BR2TC5A7jX/L7uxXz7tcgtGP+nLRhNPp19RT6TssMO8SIEFXRX3gKsNxmcESLixvghn1FKvTzQ+5jNK2uBKy2azAeWisgRjKawxSLyJxvpVZr7GuAVjOY1K1QAFVFvLisxhMEuVwFblVLVNmwuBQ4rpWqVUn7gZWCeVWOl1BNKqelKqUUYr/R2aokA1SKSB2Due22qGCxE5CbgG8ByZTYiD4BnidFc0QujMUR3u/mcFQBbRWSoFWOlVLVZ0QkBj2H9GQtTAbxsNqFuxHjz7bWzujfM5sNvAS/YTPcmjGcLjIqK5XwrpfYqpS5XSs3AEKDP+8hfb2XHoDxj8SIEm4BiESkya5rfAV4b7ETN2skTwB6l1G8GYJ8THv0hIokYBd1eK7ZKqR8rpQqUUiMxPu9flVKWascikiQiKeEwRkekpRFTSqkqoFxExpqnLgF2W7HtwUBqamXAHBHxmd/9JRhtq5YQkVxzX4hRQNhN/zWMQgJz/6pN+wEjIlcC/wQsVUq12rQtjjpcivVnbKdSKlcpNdJ8ziowOjirLKabF3X4TSw+Y1GsBhab9yrBGJhgxzPnpcBepVSFzXQrgQvN8GJsVBiinjEH8M/AIzHixSo7BucZOxs9zufDhtHmux9DgX9q0/Y5jFdXP8bD/j2LdgswmqB2ANvMbYmNdCcDn5n2pcQY2WDhPhdhY9QQRjv/dnPbNYDvayqw2cz3aiDDpr0PqAfSBvBZf4FRkJUCf8QcUWLR9iMM0doOXGL3mQCygDUYBcMaINOm/TfNcAdQDbxtw/YgRj9Y+DmLNfKnN9tV5ve1A3gdyB/I/4Cmrs7vAAACx0lEQVQ+RprFSPePwE4z3deAPJvflwf4k5n3rcBiO/kGVgC3D+B3XgBsMZ+TDcAMG7Z3Y5RD+4EHMb079GLba9lh5xmzs2kXExqNRhPnxEvTkEaj0WhioIVAo9Fo4hwtBBqNRhPnaCHQaDSaOEcLgUaj0cQ5Wgg0cYeINJv7kSLyt2f53j/pcfzJ2by/RjMYaCHQxDMjAVtCICLOfqJ0EwKllOWZzRrNuUILgSaeeRDDadk2MdYxcIrh13+T6QjtHwBE5CLTN/yzGBOgEJHVpkO+XWGnfCLyIIbn020i8ox5Lvz2Iea9S01f9NdH3XutdK3f8EzYX46IPCgiu828/OoL/3Y0cYPrXGdAozmHPIDhD/8bAGaB3qiUmiUiXuBjEQl7L50NTFRKHTaPb1FKnTBdf2wSkVVKqQdE5E5lOAnsybcwZlxPwfCHs0lEwu6IpwEXYLgu+BiYLyK7MWYaj1NKKYmx0IxGczbQbwQaTReXAzeabr83YEznD/vh2RglAgB3ich2DN//w6PixWIB8JwynKxVY6x3EHabvFEpVaEM52vbMJqsmoB24HER+RZgy3+QRmMHLQQaTRcC/EApNdXcilTXegYtkUgiF2E4LJurlJqC4Q+qvyUx+3KP3BEVDmKsNBbAeAtZhbH4yFu2PolGYwMtBJp45hSGT/gwbwN3mO5/EZGSGIvqpAEnlVKtIjIOwx9+GH/YvgcfAteb/RA5GKu4bYyVMdMPfZpS6g2MhVB6a27SaM4Kuo9AE8/sAAJmE88K4LcYzTJbzQ7bWnpfCvAt4HYR2QHsw2geCvMosENEtiqllkedfwWYi+GxUgE/UkpVmULSGynAqyKSgPE2ce/APqJG0z/a+6hGo9HEObppSKPRaOIcLQQajUYT52gh0Gg0mjhHC4FGo9HEOVoINBqNJs7RQqDRaDRxjhYCjUajiXP+P+TTOhkgLnKeAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "x = range(ITERATIONS + 1)\n",
    "for node in range(n):\n",
    "    plt.plot(x, prob[node], label=\"node %s\" % node)\n",
    "\n",
    "plt.xticks(x)\n",
    "plt.title(\"Random Surfer Probabilities\")\n",
    "plt.xlabel(\"Iterations\")\n",
    "plt.ylabel(\"Probability\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Repeated squaring\n",
    "Exploiting that for $k$ being a power of two, \n",
    "$$M \\cdot (M \\cdots (M \\cdot (M \\cdot p_\\mathrm{init}))) \n",
    "= M^k \\cdot p_\\mathrm{init}\n",
    "= M^{2^{\\log k}} \\cdot p_\\mathrm{init}\n",
    "= (\\cdots(((M^2)^2)^2)\\cdots)^2 \\cdot p_\\mathrm{init}$$ \n",
    "we can compute the PageRank vector with a reduced number of matrix multiplications."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.03935185]\n",
      " [0.35332637]\n",
      " [0.02777778]\n",
      " [0.32221711]\n",
      " [0.16203446]\n",
      " [0.09529243]]\n"
     ]
    }
   ],
   "source": [
    "from math import ceil, log\n",
    "\n",
    "MP = M\n",
    "for _ in range(1 + int(log(ITERATIONS, 2))):\n",
    "    MP = MP @ MP\n",
    "    \n",
    "p = MP @ p_0 \n",
    "\n",
    "print(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computing eigen vector\n",
    "\n",
    "We want to find the vector $p$ where $M p = p$, i.e. we want to find an _eigen vector_ $p$ for the _eigen value_ $\\lambda=1$, where $|p|=1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.03935185 0.3533267  0.02777778 0.32221669 0.16203473 0.09529225]\n"
     ]
    }
   ],
   "source": [
    "eigenvalues, eigenvectors = np.linalg.eig(M)\n",
    "\n",
    "#print(\"eigenvalues =\", eigenvalues)\n",
    "#print(\"eigenvectors =\", eigenvectors)\n",
    "idx = eigenvalues.argmax()  # find the largest eigen value (= 1)\n",
    "p =  np.real(eigenvectors[:, idx])  # .real returns the real part of complex numbers\n",
    "p /= p.sum()  # normalize p to have sum 1\n",
    "\n",
    "print(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Note on practicality\n",
    "In practice an explicit matrix for billions of nodes is infeassable, since the number of entries would be order of $10^{18}$. Instead one has to work with sparse matrices (in Python modul <tt>scipy.sparse</tt>) and stay with repeated multiplication."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time considerations\n",
    "\n",
    "You can arrive at the same results in many ways, but some are more efficient than others. In the example below a matrix is created with all values equal to 1/n in two different ways:\n",
    "\n",
    "  * First 1 / n is computed, and then filled into all entries of a new matrix.\n",
    "  * First a matrix is construct containing 1 in each entry, and then all entries are divided by the value n. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.59 µs ± 154 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "np.full((n, n), 1 / n)  # n x n matrix with 1/n in all entries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.86 µs ± 187 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "np.ones((n, n)) / n  # same as above, just slower"
   ]
  }
 ],
 "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.8.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}