{
 "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\">\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": 10,
   "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=1)  # creates a 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=1)\n",
    "\n",
    "print(G)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulate random walk (random surfer model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1], [3], [0, 1], [1, 4], [1, 5], [1]]\n",
      "[0.039417, 0.352732, 0.027955, 0.321849, 0.162372, 0.095675]\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 = [[col for col in range(n) if G[row, col]] for row in range(n)]\n",
    "\n",
    "count = [0] * 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\n",
    "    if randint(1, 6) == 6:  # original paper uses 15% instead of 1/6\n",
    "        state = randint(0, 5)\n",
    "    else:\n",
    "        d = len(adjacency_list[state])\n",
    "        state = adjacency_list[state][randint(0, d - 1)]\n",
    "\n",
    "print(adjacency_list, [c / STEPS for c in count], sep=\"\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAEWCAYAAACwtjr+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAHjlJREFUeJzt3X2YH2V97/H3xwBK5SFBFk5MYoMaLQ+VAGuIB2sRNARQk1awYI9JkTaKQfCy19Fw2sqT9OCpghc+xEYJJJQS8YGSajCmPEqFkA2EQAicrCGaNZEsJIQgEk7i9/wx95Zh+e1vf5vde2fzy+d1XXP9Zr5zzz3fwUu+3DP3zigiMDMzy+k1VSdgZmbNz8XGzMyyc7ExM7PsXGzMzCw7FxszM8vOxcbMzLJzsTEbRJJOlNRRdR59Jel6SV9M67vlNVi1XGxsjydpnaTfSXpe0m/Sv1j3qzqv/pD0z5K+WdreW9Jve4hNrCZL25O42JgVPhgR+wHjgWOAiyrOp7/uAf60tN0K/Ap4T7cYwPLBSsr2XC42ZiUR8RtgMUXRAUDS6ZIekvScpPWSLintGyspJE2X9CtJT0v6u9L+fdNIaYukx4B3ls8n6XBJd0l6VtIqSR8q7bte0jcl3ZZGXf8p6b9J+mrq73FJx/RwKXcDh0s6OG3/CbAAeH232H0R8f/S+b6XRnZbJd0j6chG/plJukDSY5JGN9Le9kwuNmYl6V+YpwLtpfBvgWnAcOB04DxJU7sd+m7g7cDJwBckHZ7iFwNvScspwPTSufYG/h34KXAI8GngRklvL/X7EeDvgYOB7cB9wINp+/vAVbWuIyI6gF9SFBQoRjQ/A37eLXZP6bDbgHEplweBG2v1XSbpH4C/Av40ndOsJhcbs8K/SdoGrAc2URQJACLiroh4JCJ+HxErgZt45S0qgEsj4ncR8TDwMHB0in8EuCIiNkfEeuCa0jETgf2AKyPipYi4A/gRcHapzS0RsTwiXgRuAV6MiPkRsRP4LsUtv57cDbxH0muACcD9FAWnK3ZCatN1nXMjYltEbAcuAY6WdGAPfUvSVRQF9L0R0VknDzMXG7NkakTsD5wI/BHFyAEAScdLulNSp6StwCfL+5PflNZfoCgiAG+kKGBdfllafyOwPiJ+323/qNL2U6X139XYrjeR4R6K0csfA2sj4gXg3lJsX2BpusZhkq6U9AtJzwHrUh/dr7PLcGAG8L8jYmudHMwAFxuzV4iIu4HrgS+Xwv8KLATGRMSBwLcANdjlRmBMaftNpfUNwJg0yijv/3Uf0+7JPRQjrNMpRjQAq1I+pwPL0ogJ4KPAFOB9wIHA2BTv6Tq3AB8ArpN0wgDla03Mxcbs1b4KvF9S1ySB/YHNEfGipAkU/2Ju1M3ARZJGpOdBny7tW0rxPOhzaRryicAHKR7k91tEtFOMhC4kFZsovimyNMXKz2v2p3gm9AzwB8A/NtD/XcBfArdIOn4gcrbm5WJj1k16/jAf+IcU+hRwWXqm8wWKAtKoSylujT1JMRHghtJ5XgI+RDEh4Wngm8C0iHi8v9dQcg/QAvxnKfYzikkA5WIzP+X5a+Axiuc7vYqIJcA5wEJJxw1Ewtac5I+nmZlZbh7ZmJlZdi42ZmaWnYuNmZll52JjZmbZ7VV1AkPFwQcfHGPHjq06DTOz3cry5cufjoiW3tq52CRjx46lra2t6jTMzHYrkn7ZeyvfRjMzs0HgYmNmZtm52JiZWXYuNmZmlp2LjZmZZediY2Zm2bnYmJlZdtmKjaTXSXpA0sOSVkm6NMWvl/SkpBVpGZ/iknSNpHZJKyUdW+pruqQ1aSl/w/04SY+kY66RpBQ/SNKS1H6JpBG5rtPMzHqXc2SzHTgpIo4GxgOTJU1M+/5nRIxPy4oUOxUYl5YZwGwoCgfF9+CPp/iO+sWl4jE7te06bnKKzwJuj4hxwO1p28zMKpLtDQLpi4DPp82901Lv4zlTgPnpuPslDZc0kuKb8EsiYjOApCUUhesu4ICIuC/F5wNTgdtSXyemfucBdwGfH6hr2xOMnfXjqlNoyLorT686BTNrQNZnNpKGSVoBbKIoGEvTrivSrbKrJb02xUYB60uHd6RYvXhHjTjAoRGxESD9HjKAl2VmZn2UtdhExM6IGA+MBiZIOgq4CPgj4J3AQbw84lCtLnYh3jBJMyS1SWrr7Ozsy6FmZtYHgzIbLSKepbiVNTkiNkZhO3AdxXMYKEYmY0qHjQY29BIfXSMO8FS6BUf63dRDXnMiojUiWltaen1pqZmZ7aKcs9FaJA1P6/sC7wMeLxUBUTxjeTQdshCYlmalTQS2pltgi4FJkkakiQGTgMVp3zZJE1Nf04BbS311zVqbXoqbmVkFcn5iYCQwT9IwiqJ2c0T8SNIdklooboOtAD6Z2i8CTgPagReAcwAiYrOky4Flqd1lXZMFgPOA64F9KSYG3JbiVwI3SzoX+BVwZrarNDOzXuWcjbYSOKZG/KQe2gcws4d9c4G5NeJtwFE14s8AJ/cxZTMzy8RvEDAzs+xcbMzMLDsXGzMzy87FxszMsnOxMTOz7FxszMwsOxcbMzPLzsXGzMyyc7ExM7PsXGzMzCw7FxszM8su54s4zSwjf03Vdice2ZiZWXYuNmZmlp2LjZmZZediY2Zm2bnYmJlZdi42ZmaWnYuNmZll52JjZmbZudiYmVl22YqNpNdJekDSw5JWSbo0xQ+TtFTSGknflbRPir82bben/WNLfV2U4k9IOqUUn5xi7ZJmleI1z2FmZtXIObLZDpwUEUcD44HJkiYCXwKujohxwBbg3NT+XGBLRLwVuDq1Q9IRwFnAkcBk4JuShkkaBnwDOBU4Ajg7taXOOczMrALZik0Unk+be6clgJOA76f4PGBqWp+Stkn7T5akFF8QEdsj4kmgHZiQlvaIWBsRLwELgCnpmJ7OYWZmFcj6zCaNQFYAm4AlwC+AZyNiR2rSAYxK66OA9QBp/1bgDeV4t2N6ir+hzjm65zdDUpukts7Ozv5cqpmZ1ZG12ETEzogYD4ymGIkcXqtZ+lUP+wYqXiu/ORHRGhGtLS0ttZqYmdkAGJTZaBHxLHAXMBEYLqnr0wajgQ1pvQMYA5D2HwhsLse7HdNT/Ok65zAzswrknI3WIml4Wt8XeB+wGrgTOCM1mw7cmtYXpm3S/jsiIlL8rDRb7TBgHPAAsAwYl2ae7UMxiWBhOqanc5iZWQVyfjxtJDAvzRp7DXBzRPxI0mPAAklfBB4Crk3trwVukNROMaI5CyAiVkm6GXgM2AHMjIidAJLOBxYDw4C5EbEq9fX5Hs5hZmYVyFZsImIlcEyN+FqK5zfd4y8CZ/bQ1xXAFTXii4BFjZ7DzMyq4TcImJlZdi42ZmaWnYuNmZll52JjZmbZudiYmVl2LjZmZpadi42ZmWXnYmNmZtm52JiZWXYuNmZmlp2LjZmZZediY2Zm2bnYmJlZdi42ZmaWnYuNmZll52JjZmbZudiYmVl2LjZmZpadi42ZmWXnYmNmZtllKzaSxki6U9JqSaskXZjil0j6taQVaTmtdMxFktolPSHplFJ8coq1S5pVih8maamkNZK+K2mfFH9t2m5P+8fmuk4zM+tdzpHNDuBvI+JwYCIwU9IRad/VETE+LYsA0r6zgCOBycA3JQ2TNAz4BnAqcARwdqmfL6W+xgFbgHNT/FxgS0S8Fbg6tTMzs4r0qdhIGiHpHY20jYiNEfFgWt8GrAZG1TlkCrAgIrZHxJNAOzAhLe0RsTYiXgIWAFMkCTgJ+H46fh4wtdTXvLT+feDk1N7MzCrQa7GRdJekAyQdBDwMXCfpqr6cJN3GOgZYmkLnS1opaa6kESk2ClhfOqwjxXqKvwF4NiJ2dIu/oq+0f2tq3z2vGZLaJLV1dnb25ZLMzKwPGhnZHBgRzwF/DlwXEccB72v0BJL2A34AfCb1Mxt4CzAe2Ah8patpjcNjF+L1+nplIGJORLRGRGtLS0vd6zAzs13XSLHZS9JI4CPAj/rSuaS9KQrNjRHxQ4CIeCoidkbE74FvU9wmg2JkMqZ0+GhgQ53408BwSXt1i7+ir7T/QGBzX3I3M7OB00ixuRRYTPHcZJmkNwNrejsoPSO5FlgdEVeV4iNLzf4MeDStLwTOSjPJDgPGAQ8Ay4BxaebZPhSTCBZGRAB3Amek46cDt5b6mp7WzwDuSO3NzKwCe/XehI0R8V+TAiJibYPPbE4APgY8ImlFiv0vitlk4ylua60DPpH6XSXpZuAxiplsMyNiJ4Ck8ykK3jBgbkSsSv19Hlgg6YvAQxTFjfR7g6R2ihHNWQ3ka2ZmmTRSbL4GHNtA7BUi4l5qPztZVOeYK4ArasQX1TouItby8m24cvxF4Mx6+ZmZ2eDpsdhIehfw34EWSZ8t7TqAYoRhZmbWkHojm32A/VKb/Uvx53j5OYmZmVmveiw2EXE3cLek6yPil4OYk5mZNZl6t9G+GhGfAb4uqdbfqHwoa2ZmZtY06t1GuyH9fnkwEjEzs+ZV7zba8vR7d1csvVpmTESsHITczMysSQzKu9HMzGzPlv3daGZmZlnfjWZmZgaNFZvL2IV3o5mZmXXp9XU1EfE94Hul7bXAh3MmZWZmzaXe39l8LiL+j6SvUftbMBdkzczMzJpGvZHN6vTbNhiJmJlZ86r3dzb/nlZXRsRDg5SPmZk1oUYmCFwl6XFJl0s6MntGZmbWdHotNhHxXuBEoBOYI+kRSX+fOzEzM2sejYxsiIjfRMQ1wCeBFcAXsmZlZmZNpZHX1Rwu6RJJjwJfB34OjM6emZmZNY1GPgt9HXATMCkiNmTOx8zMmlAjf9Q5cTASMTOz5tXQM5tdIWmMpDslrZa0StKFKX6QpCWS1qTfESkuSddIape0UtKxpb6mp/ZrJE0vxY9LExba07Gqdw4zM6tGtmID7AD+NiIOByYCMyUdAcwCbo+IccDtaRvgVGBcWmYAs6EoHMDFwPHABODiUvGYndp2HTc5xXs6h5mZVaDHYiPphvR74a50HBEbI+LBtL6N4o0Eo4ApwLzUbB4wNa1PAeZH4X5geHrb9CnAkojYHBFbgCXA5LTvgIi4LyICmN+tr1rnMDOzCtQb2Rwn6Q+Bj0sakW5N/dfSl5NIGgscAywFDo2IjVAUJOCQ1GwUsL50WEeK1Yt31IhT5xxmZlaBehMEvgX8BHgzsBxQaV+keK8k7Qf8APhMRDyXHqvUbFojFrsQb5ikGRS34XjTm97Ul0PNzKwPehzZRMQ16XnL3Ih4c0QcVloaLTR7UxSaGyPihyn8VLoFRvrdlOIdwJjS4aOBDb3ER9eI1ztH92ucExGtEdHa0tLSyCWZmdkuaOR1NedJOlrS+Wl5RyMdp5lh1wKrI+Kq0q6FQNeMsunAraX4tDQrbSKwNd0CWwxMSrfyRgCTgMVp3zZJE9O5pnXrq9Y5zMysAo28QeAC4EaK5x6HADdK+nQDfZ8AfAw4SdKKtJwGXAm8X9Ia4P1pG2ARsBZoB74NfAogIjYDlwPL0nJZigGcB3wnHfML4LYU7+kcZmZWgUbeIPDXwPER8VsASV8C7gO+Vu+giLiX2s9VAE6u0T6AmT30NReYWyPeBhxVI/5MrXOYmVk1Gvk7GwE7S9s76bmImJmZvUqj70ZbKumWtD2V4lmMmZlZQxp5N9pVku4C3k0xojnHX+40M7O+aGRkQ3oTwIOZczEzsyaV891oZmZmgIuNmZkNgrrFRtIwSf8xWMmYmVlzqltsImIn8IKkAwcpHzMza0KNTBB4EXhE0hLgt13BiLggW1ZmZtZUGik2P06LmZnZLmnk72zmSdoXeFNEPDEIOZmZWZNp5EWcHwRWUHzbBknjJS3MnZiZmTWPRqY+XwJMAJ4FiIgVwGEZczIzsybTSLHZERFbu8X69EVMMzPbszUyQeBRSR8FhkkaB1wA/DxvWmZm1kwaGdl8GjgS2A7cBDwHfCZnUmZm1lwamY32AvB36aNpERHb8qdlZmbNpNdiI+mdFF/J3D9tbwU+HhHLM+dmZnuQsbN2jz/nW3fl6VWnsFtq5JnNtcCnIuJnAJLeTfFBtXfkTMzMzJpHI89stnUVGoCIuBfwrTQzM2tYjyMbScem1Qck/TPF5IAA/gK4K39qZmbWLOqNbL6SlvHA24CLKf7A83DgXb11LGmupE2SHi3FLpH0a0kr0nJaad9FktolPSHplFJ8coq1S5pVih8maamkNZK+K2mfFH9t2m5P+8c2+M/CzMwy6XFkExHv7Wff1wNfB+Z3i18dEV8uByQdAZxFMcX6jcB/SHpb2v0N4P1AB7BM0sKIeAz4UuprgaRvAecCs9Pvloh4q6SzUru/6Oe1mJlZPzQyG204MA0YW27f2ycGIuKePowqpgALImI78KSkdopX5AC0R8TalMsCYIqk1cBJwEdTm3kUo67Zqa9LUvz7wNclKSL81gMzs4o0MkFgEUWheQRYXlp21fmSVqbbbCNSbBSwvtSmI8V6ir8BeDYidnSLv6KvtH9rav8qkmZIapPU1tnZ2Y9LMjOzehqZ+vy6iPjsAJ1vNnA5xUSDyymeCX0cUI22Qe1iGHXa08u+VwYj5gBzAFpbWz3yMTPLpJGRzQ2S/kbSSEkHdS27crKIeCoidkbE74Fv8/Ktsg5gTKnpaGBDnfjTwHBJe3WLv6KvtP9AYPOu5GtmZgOjkWLzEvBPwH28fAutbVdOJmlkafPPgK6ZaguBs9JMssOAccADwDJgXJp5tg/FJIKF6fnLncAZ6fjpwK2lvqan9TOAO/y8xsysWo3cRvss8NaIeLovHUu6CTgROFhSB8XU6RMljae4rbUO+ARARKySdDPwGLADmBkRO1M/5wOLgWHA3IhYlU7xeWCBpC8CD1G86YD0e0OaZLCZokCZmVmFGik2q4AX+tpxRJxdI3xtjVhX+yuAK2rEF1FMUugeX8vLt+HK8ReBM/uUrJmZZdVIsdkJrJB0J8VnBoDepz6bmZl1aaTY/FtazMzMdkkj37OZNxiJmJlZ82rkDQJPUuPvVCLizVkyMjOzptPIbbTW0vrrKB6+79Lf2ZiZ2Z6p17+ziYhnSsuvI+KrFO8lMzMza0gjt9GOLW2+hmKks3+2jMzMrOk0chvtK6X1HRR/jPmRLNmYmVlTamQ2Wn+/a2NmZnu4Rm6jvRb4MK/+ns1l+dIyM7Nm0shttFspvgmznNIbBMzMzBrVSLEZHRGTs2diZmZNq5FPDPxc0h9nz8TMzJpWIyObdwN/ld4ksJ3iS5gREe/ImpmZmTWNRorNqdmzMDOzptbI1OdfDkYiZmbWvBp5ZmNmZtYvLjZmZpadi42ZmWXnYmNmZtllKzaS5kraJOnRUuwgSUskrUm/I1Jckq6R1C5pZflN05Kmp/ZrJE0vxY+T9Eg65hpJqncOMzOrTs6RzfVA9zcPzAJuj4hxwO1pG4rp1ePSMgOYDUXhAC4GjgcmABeXisfs1LbruMm9nMPMzCqSrdhExD3A5m7hKcC8tD4PmFqKz4/C/cBwSSOBU4AlEbE5IrYAS4DJad8BEXFfRAQwv1tftc5hZmYVGexnNodGxEaA9HtIio8C1pfadaRYvXhHjXi9c7yKpBmS2iS1dXZ27vJFmZlZfUNlgoBqxGIX4n0SEXMiojUiWltaWvp6uJmZNWiwi81T6RYY6XdTincAY0rtRgMbeomPrhGvdw4zM6tII+9GG0gLgenAlen31lL8fEkLKCYDbI2IjZIWA/9YmhQwCbgoIjZL2iZpIrAUmAZ8rZdzmJkNqrGzflx1Cg1Zd+Xp2c+RrdhIugk4EThYUgfFrLIrgZslnQv8CjgzNV8EnAa0Ay8A5wCkonI5sCy1uywiuiYdnEcx421f4La0UOccZmZWkWzFJiLO7mHXyTXaBjCzh37mAnNrxNuAo2rEn6l1DjMzq85QmSBgZmZNzMXGzMyyc7ExM7PsXGzMzCw7FxszM8vOxcbMzLJzsTEzs+xcbMzMLDsXGzMzy87FxszMsnOxMTOz7FxszMwsOxcbMzPLzsXGzMyyc7ExM7PsXGzMzCw7FxszM8vOxcbMzLJzsTEzs+xcbMzMLLtKio2kdZIekbRCUluKHSRpiaQ16XdEikvSNZLaJa2UdGypn+mp/RpJ00vx41L/7elYDf5VmplZlypHNu+NiPER0Zq2ZwG3R8Q44Pa0DXAqMC4tM4DZUBQn4GLgeGACcHFXgUptZpSOm5z/cszMrCdD6TbaFGBeWp8HTC3F50fhfmC4pJHAKcCSiNgcEVuAJcDktO+AiLgvIgKYX+rLzMwqUFWxCeCnkpZLmpFih0bERoD0e0iKjwLWl47tSLF68Y4a8VeRNENSm6S2zs7Ofl6SmZn1ZK+KzntCRGyQdAiwRNLjddrWet4SuxB/dTBiDjAHoLW1tWYbMzPrv0pGNhGxIf1uAm6heObyVLoFRvrdlJp3AGNKh48GNvQSH10jbmZmFRn0YiPp9ZL271oHJgGPAguBrhll04Fb0/pCYFqalTYR2Jpusy0GJkkakSYGTAIWp33bJE1Ms9CmlfoyM7MKVHEb7VDgljQbeS/gXyPiJ5KWATdLOhf4FXBmar8IOA1oB14AzgGIiM2SLgeWpXaXRcTmtH4ecD2wL3BbWszMrCKDXmwiYi1wdI34M8DJNeIBzOyhr7nA3BrxNuCofidrZmYDYihNfTYzsyZV1Wy0pjJ21o+rTqEh6648veoUzGwP5ZGNmZll52JjZmbZudiYmVl2LjZmZpadi42ZmWXn2Wi2R/CMQbNqeWRjZmbZudiYmVl2LjZmZpadi42ZmWXnYmNmZtm52JiZWXYuNmZmlp2LjZmZZediY2Zm2bnYmJlZdi42ZmaWnYuNmZll52JjZmbZNW2xkTRZ0hOS2iXNqjofM7M9WVMWG0nDgG8ApwJHAGdLOqLarMzM9lxNWWyACUB7RKyNiJeABcCUinMyM9tjKSKqzmHASToDmBwRf522PwYcHxHnd2s3A5iRNt8OPDGoidZ3MPB01UkMsGa7pma7Hmi+a2q264Ghd01/GBEtvTVq1i91qkbsVVU1IuYAc/Kn03eS2iKiteo8BlKzXVOzXQ803zU12/XA7ntNzXobrQMYU9oeDWyoKBczsz1esxabZcA4SYdJ2gc4C1hYcU5mZnuspryNFhE7JJ0PLAaGAXMjYlXFafXVkLy910/Ndk3Ndj3QfNfUbNcDu+k1NeUEATMzG1qa9TaamZkNIS42ZmaWnYvNENRsr9qRNFfSJkmPVp3LQJA0RtKdklZLWiXpwqpz6g9Jr5P0gKSH0/VcWnVOA0HSMEkPSfpR1bkMBEnrJD0iaYWktqrz6Ss/sxli0qt2/i/wfoop3MuAsyPisUoT6wdJ7wGeB+ZHxFFV59NfkkYCIyPiQUn7A8uBqbvr/0aSBLw+Ip6XtDdwL3BhRNxfcWr9IumzQCtwQER8oOp8+kvSOqA1IobSH3Q2zCOboafpXrUTEfcAm6vOY6BExMaIeDCtbwNWA6OqzWrXReH5tLl3Wnbr/wqVNBo4HfhO1blYwcVm6BkFrC9td7Ab/4us2UkaCxwDLK02k/5Jt5xWAJuAJRGxW18P8FXgc8Dvq05kAAXwU0nL06u2disuNkNPQ6/asepJ2g/4AfCZiHiu6nz6IyJ2RsR4irdtTJC0297ulPQBYFNELK86lwF2QkQcS/E2+5np9vRuw8Vm6PGrdnYD6dnGD4AbI+KHVeczUCLiWeAuYHLFqfTHCcCH0jOOBcBJkv6l2pT6LyI2pN9NwC0Ut9x3Gy42Q49ftTPEpQfq1wKrI+KqqvPpL0ktkoan9X2B9wGPV5vVrouIiyJidESMpfj/zx0R8T8qTqtfJL0+TUZB0uuBScBuNbvTxWaIiYgdQNerdlYDN++Gr9p5BUk3AfcBb5fUIencqnPqpxOAj1H8F/OKtJxWdVL9MBK4U9JKiv/YWRIRTTFduIkcCtwr6WHgAeDHEfGTinPqE099NjOz7DyyMTOz7FxszMwsOxcbMzPLzsXGzMyyc7ExM7PsXGzMdjOSxjbLG7Rtz+FiY2Zm2bnYmFUgjU5WS/p2+obMTyXtK2m8pPslrZR0i6QRqf1x6Xsz9wEzS/0Mk/RPkpalYz5R2UWZ1eFiY1adccA3IuJI4Fngw8B84PMR8Q7gEeDi1PY64IKIeFe3Ps4FtkbEO4F3An8j6bBByd6sD1xszKrzZESsSOvLgbcAwyPi7hSbB7xH0oHd4jeU+pgETEufB1gKvIGiiJkNKXtVnYDZHmx7aX0nMLyHdqLnz0wI+HRELB7IxMwGmkc2ZkPHVmCLpD9J2x8D7k6v/d8q6d0p/pelYxYD56VPHiDpbemtwGZDikc2ZkPLdOBbkv4AWAuck+LnAHMlvUBRYLp8BxgLPJg+fdAJTB28dM0a47c+m5lZdr6NZmZm2bnYmJlZdi42ZmaWnYuNmZll52JjZmbZudiYmVl2LjZmZpbd/wdbAsPOjR+rpAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1859cb0>"
      ]
     },
     "metadata": {},
     "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": 15,
   "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",
    "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": 31,
   "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 = 5 / 6 * A.T + 1 / (6 * n) * np.ones((n, n))\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",
    "    \n",
    "    prob = np.append(prob, p, axis=1)\n",
    "  \n",
    "print(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzsnXd8XWX5wL/PXUludtKkWZ3Q0m3pYIiUglCGDBFEQJBhGYoTRQUVQRAZ/hRFQJGNCqIgIJaCyGiFMlpoCxRKS0nbpM1us5O73t8f77k3N8mdSW4Smvf7+dzP2ee899zkPOfZopTCYDAYDAYA20gPwGAwGAyjByMUDAaDwRDCCAWDwWAwhDBCwWAwGAwhjFAwGAwGQwgjFAwGg8EQwggFw7AhIktFpGqkxzEQROQAEXlbRFpF5FsjPZ6+iIgSkf0HeGyliBwdZdvhIrI50r4icpWI3B3jvF8WkecGMibDyGGEwhjH+ifvFJE2EakRkftFJGukxzVYROQUEVkvIi0i0iAi/xWRyYM45Q+Al5RS2Uqp3w3B+M4XEb9131ussZ442PMONUqp1UqpA6Jsu0EptRxARCZbgskRtv0vSqllwzVWw9BghIIB4CSlVBYwHzgQuHKExzMorDfmB4HvAbnAFOAOIDCAcwUfcpOA9wY4HkeUTWus+54H3AM8KiIFSRxvMAw5RigYQiilaoBn0cIBABH5nGU2aRGRnSJyTdi24NvheSKyw3oj/3HY9gxL89gjIpuAxeHXE5GZIvKSiOwVkfdE5OSwbfeLyB0i8oz1Nv2KiJSIyK3W+T4QkQOjfJX5wMdKqf8qTatS6jGl1I6wc18fdq1eZi1Le/qhiGwE2kXkBeBI4PfWWKaLSJqI/Mr63rUi8gcRyQg/n3WOGuC+OPc9ANwLZABTox0vIheJyFYRaRKRp0SkrM+pThCRbdbvcIuI2Kzj9hORF0Sk0dr2FxHJ63PsYhHZZN3b+0QkPdK96fP7XSMif7YWV1nTvdY9OtTShv4Xtv8MEfmPNf7NInJG2LYTrOu3iki1iHw/1j0zpA4jFAwhRKQCOB7YGra6HfgK+m32c8DXROTzfQ79DHAA8FngahGZaa3/GbCf9TkWOC/sWk7gX8BzQDHwTeAvIhJuqjgD+AkwDugG1gBvWcv/AH4d5au8BcwQkd+IyJEDNIedZX3fPKXUUcBq4BtKqSyl1IfATcB0tADaHygHrg47vgQoQGsYF8e6kKUJLAfagC2RjheRo4Bfou9JKbAdeKTPqU4FFgELgFOAC4OXsI4tA2YCE4Br+hz7ZfRvtJ/1vX4Sa8wRWGJN86x7tKbPd8wE/gP8Ff17nwXcISKzrV3uAS5RSmUDc4AXkry+YYgwQsEA8ISItAI7gTr0wxwApdRLSql3lFIBpdRG4GHgiD7HX6uU6lRKbQA2AJ+y1p8B/EIp1aSU2gmE2+IPAbKAG5VSHqXUC8DT6IdFkH8qpdYppbqAfwJdSqkHlVJ+4G9oU1c/lFLbgKXoB/WjQMMAfCW/U0rtVEp19t0gIgJcBHzX+m6twA3AmWG7BYCfKaW6I53D4hAR2QvUWN/7VKVUc5Tjvwzcq5R6SynVjTbxHdrHT3KTNZ4dwK3WOVFKbVVK/cc6Vz1amPb9DX9vfd8m4Bf0/h2GghOBSqXUfUopn1LqLeAx4HRruxeYJSI5Sqk91nbDCGCEggHg89Yb2lJgBvpNHAAROVhEXhSRehFpBi4N325REzbfgX7Yg34z3Rm2bXvYfBmw0zKdhG8vD1uuDZvvjLAc9SGvlHpNKXWGUqoIOBz9JvvjaPtHYGeMbUWAG1hnmb72Aiut9UHqLWEWi9eUUnlKqXFKqUOUUs/HOL6MsPunlGoDGul9v/re6zIAESkWkUcss0wL8Gf6/4YRjx1CJgEHB++Xdc++jNaIAE4DTgC2i8jLInLoEF/fkCBGKBhCKKVeBu4HfhW2+q/AU8AEpVQu8Ae0OSIRdqNNFUEmhs3vAiYE7d5h26uTHHZclFJvAo+jzRKgTWLusF1K+h0EscoHN6CF0mzroZ6nlMq1nMaJHJ8IfY/fhX6wAiFzTCG971ffe73Lmv+ldb55Sqkc4Bz6/4bRjh3oePuyE3g57H4FzUxfA/0bKaVOQZuWnkBreIYRwAgFQ19uBY4RkaCzORtoUkp1ichBwNlJnOtR4EoRybf8Fd8M2/Y6+uH8AxFxishS4CT628mTRkQ+Yzlli63lGcDJwGvWLuvRTtkCESkBvpPM+S3t5k/Ab8KuUS4ixw527DH4K3CBiMwXkTS0uep1pVRl2D5XWPd6AvBttIkN9G/YhnYClwNXRDj/ZSJSITr66aqwYxOlHm3ymhpl+9PAdBE51/q9nSKyWHSwgUt0TkOuUsoLtAD+JK9vGCKMUDD0wrI5Pwj81Fr1deDnls/hapJ7g7sWbYr4GO1QfijsOh70g/p49Jv3HcBXlFIfDPY7AHutc78jIm1o084/gZut7Q+hfR+V1riSfQAC/BDtkH/NMsk8j3a2pwSl1H/Rv8ljaA1sP3r7MACeBNahhd6/0c5b0L/DAqDZWv94hEv8FX0vtlmf6yPsE2t8HWhfxCuWeeiQPttbgWXWmHehTY43AWnWLucClda9vBStzRhGADFNdgwGg8EQxGgKBoPBYAhhhILBYDAYQhihYDAYDIYQRigYDAaDIcQnrtDWuHHj1OTJk0d6GAaDwfCJYt26dQ1WMmdMPnFCYfLkyaxdu3akh2EwGAyfKERke/y9jPnIYDAYDGEYoWAwGAyGEEYoGAwGgyHEJ86nYDAYDEG8Xi9VVVV0dcUrSDt2SE9Pp6KiAqfTOaDjjVAwGAyfWKqqqsjOzmby5MnoNhdjG6UUjY2NVFVVMWXKlAGdI2XmIxG5V0TqROTdKNtFRH5ntRfcKCILUjUWg8Gwb9LV1UVhYaERCBYiQmFh4aA0p1T6FO4Hjoux/XhgmvW5GLgzhWMxGAz7KEYg9Gaw9yNlQkEptQpoirHLKcCDVmP114A8ESlN1XjerGzippUfYKrCGgwGQ3RGMvqonN4tAKvo3VowhIhcLCJrRWRtfX39gC62Yede7nzpI1q6fAM63mAwGIaayspK5syZE3/HMFauXMkBBxzA/vvvz4033jjkYxpJoRBJx4n4Gq+UuksptUgptaioKG6WdkQKMl0A7Gn3DOh4g8FgGGn8fj+XXXYZzzzzDJs2beLhhx9m06ZNQ3qNkRQKVfTuC1tB8n1hEybfEgqNRigYDIYhorKykpkzZ3LRRRcxe/Zsli1bRmdnJwDr16/nkEMOYd68eZx66qns2bMHgHXr1vGpT32KQw89lNtvvz10Lr/fzxVXXMHixYuZN28ef/zjH/td74033mD//fdn6tSpuFwuzjzzTJ588skh/U4jGZL6FPANEXkEOBhoVkrtTtXFCtxGUzAY9mWu/dd7bNrVMqTnnFWWw89Omh1zny1btvDwww/zpz/9iTPOOIPHHnuMc845h6985SvcdtttHHHEEVx99dVce+213HrrrVxwwQWh9Vdc0dMu+5577iE3N5c333yT7u5uDjvsMJYtW9YrtLS6upoJE3repSsqKnj99deH9DunMiT1YWANcICIVInIV0XkUhG51NplBboX7FZ0E/Svp2os0GM+auowQsFgMAwdU6ZMYf78+QAsXLiQyspKmpub2bt3L0cccQQA5513HqtWreq3/txzzw2d57nnnuPBBx9k/vz5HHzwwTQ2NrJly5Ze14oUKDPU0Vcp0xSUUmfF2a6Ay1J1/b6EhILRFAyGfZJ4b/SpIi0tLTRvt9tD5qNIKKWiPsSVUtx2220ce+yxUY+vqKhg586e+JyqqirKysoGMOrojJnaR26XHZfDZsxHBoMh5eTm5pKfn8/q1asBeOihhzjiiCPIy8sjNzeX//3vfwD85S9/CR1z7LHHcuedd+L1egH48MMPaW9v73XexYsXs2XLFj7++GM8Hg+PPPIIJ5988pCOfcyUuRARCtwuoykYDIZh4YEHHuDSSy+lo6ODqVOnct999wFw3333ceGFF+J2u3tpBcuXL6eyspIFCxaglKKoqIgnnnii1zkdDge///3vOfbYY/H7/Vx44YXMnj20GpJ80pK5Fi1apAbaZOeE366mLC+du89bPMSjMhgMI8H777/PzJkzR3oYo45I90VE1imlFsU7dsyYj0D7FUxIqsFgMERnTAmF/EyX8SkYDAZDDMaUUCjMND4Fg8FgiMWYEgr5bhctXT68/sBID8VgMBhGJWNKKBRk6k5Ee0wCm8FgMERkTAmF/FBRPO8Ij8RgMBhGJ2NKKJisZoPBMJoYSOnsCy+8kOLi4qSPSxQjFAwGg+ETxPnnn8/KlStTdv6xJRTcpiiewWAYOoa7dDbAkiVLKCgoSNl3GjNlLiDcp2CEgsGwz/HMj6DmnaE9Z8lcOD52d7PhLJ09HIwpTcFpt5Gd7jDmI4PBMGQMZ+ns4WBMaQqg/QpGKBgM+yBx3uhTxXCWzh4OxpSmADqBzeQpGAyGVJKq0tnDwZjTFAozXdS0dI30MAwGwz5OKkpnA5x11lm89NJLNDQ0UFFRwbXXXstXv/rVIRv3mCqdDfD9v2/g1a0NvHrlZ4dwVAaDYSQwpbMjY0pnJ0GwfPYnTRgaDAbDcDDmhEK+20W3L0Cn1z/SQzEYDIZRx5gTCoUmq9lgMBiiMuaEQv5AhELTNrhxItR/mKJRGQwGw+hgzAmFYPnspIRC7XvQ1Qw1G1M0KoPBYBgdjEGhoBNNkspVaKvT09bdKRiRwWAwjB7GnlAIFsVLpqdCe72ettakYEQGg2Gskmzp7J07d3LkkUcyc+ZMZs+ezW9/+9shH9OYS17LTndgtwlN7d2JHxQSCkZTMBgMI4fD4eD//u//WLBgAa2trSxcuJBjjjmGWbNmDdk1xpymYLMJ+W5ncppCyHxkNAWDwdDDcJfOLi0tZcGCBQBkZ2czc+ZMqqurh/Q7jTlNAXQCW1Lls4OaQsuu1AzIYDAMmpveuIkPmj4Y0nPOKJjBDw/6Ycx9Rqp0dmVlJW+//TYHH3zwkH7nMacpgE5gS6rRTrimYDKhDQZDGCNROrutrY3TTjuNW2+9lZycnCH9PmNWU9hS15b4Ae31IHbwderQ1Iy81A3OYDAMiHhv9KliuEtne71eTjvtNL785S/zhS98YWCDjsHY1BSSMR95u6C7BYpm6GXjVzAYDHFIVelspRRf/epXmTlzJpdffnlKxp5SoSAix4nIZhHZKiI/irB9ooi8KCJvi8hGETkhleMJUpipeyoEAgmYgtot01HJXD01EUgGgyEBHnjgAa644grmzZvH+vXrufrqqwFdOvuyyy7j0EMPJSMjI7T/8uXLmTVrFgsWLGDOnDlccskl+Hy+Xud85ZVXeOihh3jhhReYP38+8+fPZ8WKFUM67pSVzhYRO/AhcAxQBbwJnKWU2hS2z13A20qpO0VkFrBCKTU51nkHWzob4N7/fczPn97E2z89JlT2IipV6+Duo+DYX8KzV8Ln74T5Zw/q+gaDYWgwpbMjM1pLZx8EbFVKbVNKeYBHgFP67KOAoJckFxiW8J6CYP2jRJzNQU2hdJ6eGk3BYDDsw6RSKJQDO8OWq6x14VwDnCMiVcAK4JuRTiQiF4vIWhFZW19fP+iBBYVCQn6FYDhq3kRIzzU+BYPBsE+TSqEQycXe11Z1FnC/UqoCOAF4SET6jUkpdZdSapFSalFRUdGgB1aQTKXUYDhqZjFklxpNwWAw7NOkUihUARPClivobx76KvAogFJqDZAOjEvhmIAky2e310NaDjjTIbvEaAoGg2GfJpVC4U1gmohMEREXcCbwVJ99dgCfBRCRmWihMHj7UBxCRfES8Sm01UGmpZ1kl0KL0RQMBsO+S8qEglLKB3wDeBZ4H3hUKfWeiPxcRE62dvsecJGIbAAeBs5Xw9A8OcNlJ8NpT9ynkFWs57NLoa0GAoHUDtBgMBhGiJRmNCulVqAdyOHrrg6b3wQclsoxRKMg05VYUby2Oig6QM9nl0LABx2NkDV434bBYBjbVFZWcuKJJ/Luu+8mtH9XVxdLliyhu7sbn8/H6aefzrXXXjukYxqTGc0A+ZnOxMpnt9eFaQolemqczQaDYQRIS0vjhRdeYMOGDaxfv56VK1fy2muvDek1xqxQKMhMo6kjjqbg90Lnnt4+BTDOZoPBAAx/6WwRISsrC9A1kLxeb9RaSgNlTBbEAyhwO6lsaI+9UzBHISQUjKZgMIxWam64ge73h7Z0dtrMGZRcdVXMfYa7dLbf72fhwoVs3bqVyy67zJTOHiryM13xQ1KDQmGQ5iN/ayvK709yhAaD4ZPAcJfOttvtrF+/nqqqKt54442E/RGJMoY1BRdt3T66fX7SHPbIO7UFNQVLKNidWmtIQigoj4ePjj6Gou98m/yzzhrkqA0GQzTivdGniuEunR0kLy+PpUuXsnLlyqT6PMdjzGoKBVk6V2FvLL9CsO5ReKRRkgls3rp6/M3NdA2xWmswGEYvqSqdXV9fz969ewHo7Ozk+eefZ8aMGUM69jGtKYDOah6fkx55p/ASF0GSLHXhq9UCxLvLtPI0GMYSDzzwAJdeeikdHR1MnTqV++67D9Clsy+88ELcbncvrWD58uVUVlayYMEClFIUFRXxxBNP9Drn7t27Oe+88/D7/QQCAc444wxOPPHEIR33mBUKCZW6aK8HpxvSsnrWZZfA7g0JX8dXWwuAd4ibaxsMhpFn8uTJvWz63//+90Pz8+fPjxguunDhQjZs6HmGXHPNNQDYbDZuuOEGbrjhhqjXmzdvHm+//fYQjDw6Y9d8lIhQCC9xESS7TK/3J5D4BnhrLKGwezfDkKxtMBgMg2LMC4U9seofhSeuBckuAVSPaSkOQU1BdXXhb2oayFANBoNh2BizQiEvwwnE0xTqI2gKySWweS2hAMaEZDAYRj9jVig47DZyM5xxfAqRzEfJ5Sr4amux5+cDxtlsMBhGP2NWKAAUxkpgC/itwnd9zUdBTSExoeCtrSFjwQI9bzQFg8EwyhnTQiE/0xXdp9DRBCrQOxwVtOYg9oSEggoE8NXVkzZ1KracHLzVRlMwGAyjm7EtFNwuGtuiCIVIiWsANlvCCWz+xkbw+XCUjMdZXm40BYPB0IvKysoBZSP7/X4OPPDAIc9RgDEuFAoyndE1hUiJa0GySxLSFLy1+hzO8eNxlpUZn4LBYBgSfvvb3zJz5syUnHuMC4U09rR7I+cP9C2GF052aUKagq9ORx45xo/HWV6Gt7ra5CoYDPsQw106G6Cqqop///vfLF++PCXfacxmNIPWFDz+AO0eP1lpfW5FSFOI0GEtuwS2vxL3/N4aLTgclqYQ6Ogg0NyMPS9vsEM3GAx9WP3ohzTsbBvSc46bkMXhZ0yPuc9wl87+zne+w80330xra+uQftcgY1pTyA/WP4rkV2ivA7sL0nP7b8su1c13vNGrIQL4auvAbsdRWIizvBwwYakGw77GcJbOfvrppykuLmbhwoUp+z5jWlMotCqlNnV4mFjo7r0xmLgWqcxteAJbwZT+2y18NTU4iosRux1nWRkAnupq0mfNGpLxGwyGHuK90aeK4Syd/corr/DUU0+xYsUKurq6aGlp4ZxzzuHPf/7zwL9AH4ymAOyJlKsQKXEtSCiBLbZfwVtXi7NY+ySCQsFnNAWDYZ8nVaWzf/nLX1JVVUVlZSWPPPIIRx111JAKBBjjmkLMonjt9ZA1PvKBCSaw+WpqSZuu317seXnY3G48JizVYBgTpKJ09nAwpoVCzPLZbfUwfm7kAxPQFJRSeGtryTz8M4BuuO0sN2GpBsO+xHCXzg5n6dKlLF26dGADj8GYNh9lpzlw2oWmvrkKSlmaQhTzUUY+ONKhNfoDPtDWhurowDm+JLTOWVZuspoNBsOoZkwLBREh3+3q71Po3AMBb+TENX1g3KzmYMlsx/geE5TRFAwGw2hnTAsF0H6Fxr5CIVbiWpA4CWzB5jrOknChUE6guRl/29DGUhsMBsNQMeaFQkRNIVbiWpA4pS4iagpWBJIxIRkMhtHKmBcKBVmu/j6FUDG8QWgKtVY2c3HPOUJCYZeJQDIYDKMTIxQiagqW+SimplAKnjboaom42Vdbhz0/H1tYYksoq9loCgaDYZSSkFAQkcdE5HMiss8JkfxMF3s7vfgDYYXq2ut0z4SMgugHxmnL6autxVFS0mudvbAQSUszzmaDwQAMrHT25MmTmTt3LvPnz2fRokVDPqZEH/J3AmcDW0TkRhGZkchBInKciGwWka0i8qMo+5whIptE5D0R+WuC4xkyCjNdKAV7w01I7fWQOU73TohGnLac3tqebOYgIqJLaJsENoPBMAhefPFF1q9fz9q1a4f83AkJBaXU80qpLwMLgErgPyLyqohcICLOSMeIiB24HTgemAWcJSKz+uwzDbgSOEwpNRv4zoC/yQAJJrD16qvQVh89HDVIIprC+P4Z0aavgsGw7zASpbNTTcIZzSJSCJwDnAu8DfwF+AxwHrA0wiEHAVuVUtus4x8BTgE2he1zEXC7UmoPgFKqLvmvMDgKgpVS2709K9vroieuBYmhKQQ8HvxNTThKIgiF8nK63n9/wOM1GAyRefH+u6jbvm1Iz1k8aSpHnn9xzH2Gu3S2iLBs2TJEhEsuuYSLL449vmRJ1KfwOLAacAMnKaVOVkr9TSn1TSArymHlwM6w5SprXTjTgeki8oqIvCYix0W5/sUislZE1tbX1ycy5ITJz9SKTlN7d8/KRDSFtCxIy4moKfjqejqu9cVZVoa/qYlAjEqKBoPhk8Nwls4GXSn1rbfe4plnnuH2229n1apVQ/p9EtUU7lZKrQhfISJpSqlupVQ0T0ek+rB92445gGloTaMCWC0ic5RSe3sdpNRdwF0AixYtGtLWZYWZOjoopCkolZimAFauQn9TkC/UXKek3zZnuRWWuns3aVOnDnDUBoOhL/He6FPFcJbOBiizQtuLi4s59dRTeeONN1iyZMkARh6ZRB3N10dYtybOMVXAhLDlCqDvE7QKeFIp5VVKfQxsRguJYSPPrTWFkE+huxV8XfE1BYha6sJrJa45x/c/R09YqnE2Gwz7Kqkqnd3e3h7quNbe3s5zzz2XdPRSPGJqCiJSgjb5ZIjIgfS8/eegTUmxeBOYJiJTgGrgTHQEUzhPAGcB94vIOLQ5aWiNgnFId9rJdNlpDHZfS6TERZDsUtjRXzb6rBIXfUNSwWQ1GwxjhVSUzq6treXUU08FwOfzcfbZZ3PccRGt7gMmnvnoWOB89Fv+r8PWtwJXxTpQKeUTkW8AzwJ24F6l1Hsi8nNgrVLqKWvbMhHZBPiBK5RSjQP6JoOgIMvVoymESlyMi39gUFNQqleHNl9dLeJ2Y8vq725xFBWB02k0BYNhH2C4S2dPnTq117GpIKZQUEo9ADwgIqcppR5L9uSWH2JFn3VXh80r4HLrM2IUuF09PRWCJS4SMh+Vgd+jq6q6exLdvDW1OMePj2g7FLsdZ0mJCUs1GAyjknjmo3OUUn8GJotIvwe3UurXEQ77xJGfGS4UkjEfWeahll29hEK0HIUgzvJyoykYDIZRSTxHc6Y1zQKyI3z2CQrcrh6fQls9IOBOxHwUOYHNW1sbMRw1iElgMxiGDm1wMAQZ7P2IZz76ozW9dlBXGeUUZIb5FNrr9Fu/PYFo3QgJbMrvx1dfH0dTKMNXV0fA48Hmcg1m6AbDmCY9PZ3GxkYKCwujhnqOJZRSNDY2kp6ePuBzxDMf/S7OAL414CuPIvIzXXR4/HR5/aS31SXmT4CIvZp9jY3g80XMZg7iLNNhqb7du3FNmjTgcRsMY52KigqqqqoY6qTWTzLp6elUVFQM+Ph4r8PrBnzmTxAFYfWPSmP1Zu6LIw3chb00BV9t9GzmID19FXYZoWAwDAKn09mvDIRhcCQSfbTPk2/VP2ps81DaVgflCxM/OLu0t1CoC3Zc65+jEMQksBkMhtFKPPPRrUqp74jIv+hfogKl1MkpG9kwUpgVVim1vT6xyKMgfdpyeq0SF5GymYM4xxeDzWaczQaDYdQRz3z0kDX9VaoHMpIENYXm5mbdTS2RxLUg2SVQ+15o0VdbBw4H9sLCqIeI04mjZLzRFAwGw6gjnvlonTV9WURcwAy0xrBZKeWJdewniaBPoXOv5TBO1NEM2nzUVgsBP9js+GprcBQXIbEa9GCFpZpSFwaDYZSRaOnszwEfAb8Dfg9sFZHjUzmw4SQ3w4lNwNus/QHJmY9KQQVCSW/e2jqcxdGdzEFc5eV4dhlNwWAwjC4SrZL6f8CRSqmlSqkjgCOB36RuWMOL3SbkuV0EQnWPEow+gp4Ethb91u+rqYlYCK8vjrIyfLV1KJ8v2eEaDAZDykhUKNQppbaGLW8Dhr1LWirJdzuRZEpcBAnLVVBK4a2ri+lkDuIsKwO/H69VUdVgMBhGA/Gij75gzb4nIiuAR9E+hS+iS2PvMxRkunAEhcJANIXW3QRaW1EdHTHDUYO4gmGpu6pxVfRtSGcwGAwjQ7zoo5PC5muBI6z5eiA/JSMaIQoyXaQ3NUF6rk5KS5TMIhAbtNbgqw3mKCSoKYAJSzUYDKOKeNFHFwzXQEaagkwXbm8jFCRhOgJdIylrPLTuDpmCnAn6FMAksBkMhtFFQj2aRSQd+CowGwhVWlJKXZiicQ07+W4Xuf49qKyiiM2lY2IlsPl8QU0hfvSRzeXCUVRkNAWDwTCqSNTR/BBQgu7E9jK6E1trqgY1EhRkuiikGV96EolrQbJLobUmlM3sKE5M29B9FYxQMBgMo4dEhcL+SqmfAu1WPaTPAXNTN6zhpyDTRaG00OkqiL9zX4KaQm0d9oKChMthm74KBoNhtJGoUPBa070iMgfIBSanZEQjREE65Ek7bc6BCIVS6GjEV7M7ZsnsvjjLy/Du3o0KBJK/psFgMKSARIXCXSKSD/wUeArYBNyUslGNAMU2bQ1rlrzkD7ZyFby7qxPKZg7iLC8HrxefqQVvMBhGCQk5mpVSd1uzLwNTUzeckaOAvQA0Dkgo6EgiX109GQsXJ3ye+XSYAAAgAElEQVRYKCy1elfM/gsGg8EwXCRa+6hQRG4TkbdEZJ2I3Coi0cuAfgLJ82uhUB/ISf7g7BICfvA3tyb1cDd9FQwGw2gjUfPRI+iyFqcBpwMNwN9SNaiRIK27AYDd/oEIhVJ8nXYAHMmYj0p1NrRxNhsMhtFCQuYjoEApdV3Y8vUi8vlUDGikCNY9qvJkJX+wuwBfl86CTsbRbHO7sRcUGE3BYDCMGhLVFF4UkTNFxGZ9zgD+ncqBDTvt9XSQQV1XorckDBG8SkctJesbMGGpBoNhNBGvIF4rugCeAJcDf7Y22YA24GcpHd1w0lZHqz2PpvaB9Q7y+bKB5oTKZofjLC+n+8MPB3RNg8FgGGpivhYrpbKVUjnW1KaUclgfm1JqAMb3UUx7Pe3OggELBW93OjYn2LOSMz8FNQWl+rXA1nTuge62AY3JYDAYkiVRnwIicjKwxFp8SSn1dGqGNEK019OdNo6m5gFqCh02HO7kk9CcZWWo7m78jY04xkUosfHn00EELnwO4rT4NBgMhsGSaEjqjcC30Ulrm4BvW+v2Hdrq8GWMo6XLh9ef/MPd1+bHke5N+q0+FJYaya/QuReq10HVm7DxkaTHZDAYDMmS6KvnCcAxSql7lVL3AsdZ6/YN/D7oaCRgNdfZ2+GNc0B/vM1dODP80JZcJzVneYy+CjtfBxS4C+H5a6B7n6pBaDAYRiHJ2CPCU31zh3ogI0pHI6CwWW04k/UrKL8f3542HG4/tO5O6lhnrL4K218FmxO++IAWNqtuSercBoPBkCyJCoVfAm+LyP0i8gCwDrgh3kEicpyIbBaRrSLyoxj7nS4iSkQWJTieoaVdt5t25urIoWSFgq+hEQIBrSm0JCcU7NnZ2HJyIpfQ3rEGyubDlMPhU2fDmjug8aOkzj8oOpqG71oGg2FUEFcoiIgA/wMOAR63PocqpWIauUXEDtwOHA/MAs4SkVkR9ssGvgW8nvToh4o2LRQy8rVQ2NORpFCos5rrDEBTgGBfhT6agrcTqt+CiYfq5aN/ptuEPvvjpM8/IF6+BW7ZD97/1/Bcz2AwjAriCgWlYyWfUErtVko9pZR6UilVk8C5DwK2KqW2KaU86FIZp0TY7zrgZqArmYEPKVY2c1aBNuUkrSkEezNnu6A1kVvTm4gJbNXrIOCFSZ/Wy9klsOT78OEzsPX5pK+RFDteg5duALsLHrtIj8VgMIwJEjUfvSYiiZf/1JQDO8OWq6x1IUTkQGBCvPBWEblYRNaKyNr6VJSZtjSF7EJdiyhZoRDqzVxcNDBNoawMb3V171yF7Wv0dMLBPesO+TrkT4GVV4E/eWd4QnQ1a0GQOwG+9ipkFcHDZ8HenfGPNRgMn3gSFQpHogXDRyKyUUTeEZGNcY6J1Oo49NQTERvwG+B78S6ulLpLKbVIKbWoqKgowSEnQXs92NNwZeaRne4YmKbgdGIfXzYwTaG8jEBHB4Hm5p6VO16F4lngDmv640iDY2+Ahs3wxp/intezYwdtq1YlPhCl4OnLoaUaTrsHCveDs/+uTVl/PQO6WpL4VgaD4ZNIokLheHQfhaOAk4ATrWksqoAJYcsVQLiNJBuYA7wkIpVon8VTI+Jsbq+HrGIQoSDTlbRPwVtbg7OoCMkphdbk6xiFIpCCJiS/D3a+0eNPCOeA42G/o+ClG6G9IeZ56275FTsvuZSON99MbCAb/wbv/gOOvBImWIph8Qw44wGo3wz/uECPzWAw7LPEFAoiki4i3wGuQOcmVCultgc/cc79JjBNRKaIiAs4E921DQClVLNSapxSarJSajLwGnCyUmrtYL7QgGirAytHId/tGoCmUIdj/HirV3ONfuNOgmACmyfobK59BzxtPf6EcETguBv19heu67/dQnk8tL/6KijFrh/+CH9rnByHpm3w7+/BpMPgM5f33rbfUXDir7Uv45kfJP39DAbDJ4d4msIDwCLgHbS28H+Jnlgp5QO+ATwLvA88qpR6T0R+bpXMGD2012lNASjMHIBQqKnRJbOzS8HXBV17kzo+qCn4gppC0J8QSVMAKDoADroY1j0AuyNb8TreXk+gvZ3Ci5bjramh9hcxIoj9XnhsOdjs8IW79LQvC8+HT38L1t4Dr92R4DcbAJ4OnY/x2h90RrfBYBhW4tU+mqWUmgsgIvcAbyRzcqXUCmBFn3VXR9l3aTLnHlLa6qF0PgD5mS7e35247VwphbeujqylS0O9mmmtgYz8hM9hz8vD5nb3aAo7XoW8iZBbHv2gpT+Edx6FZ34IF6zQGkT4V1r1MjidFF5yKeJ00nDHnWQdeSQ5xy7rf66XfqkjjL74AORWRL/m0dfCno91WGz+ZJjxuai7dm/ZgmvKFMSRcHktqPsA/n4+1L+vl5+/BuaeDouX63wNg8GQcuJpCqEQF+vNf98jEOjxKQAFmS6akvApBFpaUJ2dumS21as52QgkEcFZboWlKqU1hYkRTEfhZOTDUT/RAuS9x/ttbl+1GvfChdizMhn3ta+RPmcONT/7Gd66ut47frwaVv8aDjwXZsfpm2Szwal3QdmBWrPY9XbE3Zr+8he2nXQyVd/5DgFPAvdSKXjrIbhrKXQ0wDmPwyWrYN4Z8O5jcNcR8KejYP1ftdM7VbTs0r6czj2pu4bBMMqJ9xr3KREJvjYLkGEtCzqF4ZNfPrtzDyg/ZGqhkO920eUN0OHx4XbFf8v1WjkKzvHFPZpCklnNAI6yMp3V3LhVPxgnRTEdhbPgPFh7Lzx3NUw/HlxuPabdu+nesoXiz+uHvDidlN18Mx9/4Qvs/vFPmHDXHxERnbH8+MU6yui4BOsbutxw1iNw92fhr2fCRf/tpV20vvQStb+4gbRp+9P2/H+p+uY3qfjd77ClpUU+X3erjnh651GYsgS+8Kee+3jy72DZdbDhEXjzbnjia/DsVTD/y7DoQj3uwdBWD5Wr4eNVetq4tWdbTrmO/ho/W3+KZ8G46eBwDe6a4fi9WhA1V+mpwwUZBTriLDh1RLlvBkOKiPnUU0pFMC7vY1glLsjSjubCTP1P39TuSUgohBLXSkrCzEfJCwVXeTmd6zfoekcQX1MAbfs/7ia4/wR45bc6aghoW70agKwlh4d2TZs6heIrvk/tddez95FHyD/zTPjXt7SWdNZ/IC2JPhDZ4+HsR+GeZfDXL8GFKyEtm65Nm6i+/Hukz5jBpIcepPlfT1NzzTVUfe3rVNz+e2wZGb3Ps3ujNhft+RiO/Akcfnl/f0Z6Lhx8ifahVP5PC4fX/wBrfq8d4IuXw7RjwZ6AmapzD1S+0iME6jbp9a5sLYQXng8FU6HhQ6jdpLdve0knEQLYHFow9BUWuRX9zHeAzvlortI5Hs079XxoWqX/TlScirzOTEtI5PcWFn2nIuBpD/u0WdOOsPn23vNea5vfp4WPI10LJke6XranxVgfnHdq4Rbw6vMEvNayL2x932Vrv4APxAZi1/fWZtfLNnvkdTaHtd5aB1rLVIE4nwj7BHuHiVjnijQvkdfrd2Ir4EKFBV5EWNdvWxihvxnpsxxj+4LzYP/Pxv6bGSRJGHz3UazEtVD0kSUU9rR7qUjALRASCsXjwZkB6XkDzmoONDfj/3A1dvc4GDctsQMnHwazT4VXboUDz4G8CbSvXo2jtBTX/vv32jX/7LNpe/Elam+6GXfWbtLe/xccc502ByXL+Flwxv3wlzPgHxfiXfobdl76Ney5uVT84U5smZnkn/klxOlk909+ws5Lv8aEO+/A5nbrf44379Zv/e5COO9p/T1iIaJrQE05XN/ftx6EtffBI2frt/qFF8CCr2iBFaS7VZviKldpQbB7I6DAkQETD9H+iilHaH9SL6ES5ivxebQGUbcJat/VwmLn6zp0N0harr4fBVN1mHDw4d/dxzdlc2o/Ue4Efd3cCv3Jm6BNj34PdDZpDa6zCTr29Flugr07tHDr3EtY2k+0mwauLK3duTKtT5YWInkTrG2Z+mHr67Y+XeAPm/d1a+Hm6+6z3qOnAa/18HZqAWFzWFOnvqcR1zu1QLE5rIe0HwLWR3msafi68KlPm3yVX38/sVmf8PlInz7bgd4P8UDvB3rE9WHz9BEavR7s0bbRMx/87UI/YQQhEmk5ySCWgWCEglXiImg+Ksh0AiTsV+iVzQw6AmmA9Y8AvJtexz77kMhvntE45jrY/Az856eoU+6i/ZVXyTnpJG0iCkNEKP3FL/j4pBPZdeOdTF5+BHLoN5Iea4j9j4YTbsH/xPfYefcXCLQHmPTXv+IsLg7tknfaFxCng10/upIdF13MhFtvxP7fH8IHT+s3/M/fCZmFyV03uwSO+IEOnf1wpRYwL14PL98IM0/SWd+Vq3XtKOXX5ToqDoKlV2qhUr4wcbOMw6Uf+ONnaSESpHMv1L0Pde9pQVH7Hmz9r9Y48ydpIZc7wXrwT9AP4czioWuUFPDrh3VQYCjV+8HvytQvKcn8HQ0EpVJ/DcOwYoRCUCiEHM36YdHU3p3Q4b7aWuyFhYjLsjXnlA5YUwAtZNJP6DEdvbuqmrQMB9MWj492qH7gHPYdePlGOlyfIdDR0ct01Os6hbmUHA7V/3bQ0PRpigb5kFIHnkf1DQ/RXVPDhO+dSvoB0/vtk3vyyYjTSfX3vs+O045h4pJ67Cf8QpftGMz17Q6YeaL+NH6E74U7qHvwKXztq8iZV0L2Zy/DPvtomHCQfkAOJRl52uSUiO8nFdjs+o0/PON9JDACYZ/D9Hdsq9NqbLpuF1HgDvoUEqst5K2rxTk+7IE9WE2h3R7KTwgEFGv++RGvPflR9B7OQQ77NuRU0Pbo7eBw4D74kMj7/ffn5GRvJvfIRTTc+2c6N2xIeqxBlFLUXH897e/XUvK5iWRV/wE2r+y/YyBATvZmyg9roqsBdmw8CP+sc4bsrVkFAuz5z5t8dPNqWioz8DCF3Sv3sOWqJ9n5q0dpfvYFAu3tQ3Itg2Ffx2gK7VY2s/WAyk53YLcJexJMYPPV1Ibe8vUJrKzmQCCph569sBBx2PB2ZUDJPAAadrbi6fTh6fSxZ3cHBWWZ0U/gcsOy62j/2w9wT98Pe1aEfbf+VztoF1/E+CU/o/2UU6j+wQ+Y+s9/alt/kjTdex97H/kbhRctJ/+bl8J9J8A/LoQLn4HST+md2hvgn5fA1ufJOeYU5PQzqP7+VWw//wIm3nsPjoLBvel2vf8+u6+5hq4NG3EfegglV1+Na/Jkut59j5YVK2h55hnaXngBSU8na+lSck44nqwlS7Clpw/qugCBjg66Nm2ic8NGOt95B8+2bThKxuOaPBnXpEm4Jk8mbfJkHKWlSIr6ayul8O/VdmZ7VhbidKbkOoaxgxEKbfUhJzOAzSbku500JioUamvJWBDmqM0u1XbsjoaQSSoRRARnFngDhSGnZ9UHPfHyle82xBYKgLfgELqbnRRP+0g7I8MT6NobdEhn0UxYdh12ZwZlN97IjvPOp/bmmym95pqExwrQ8uxz1N1yC9nHHUfRd7+rBWAoVPVLcNELuiHQY8v1WD73a1h0IdkiVGTmU3XZZew47zwm3nsvjgEUOfS3tdNw2200PfQQ9vx8ym65mZwTTwz5UTLmziFj7hyKr/g+nW+/Tcu/V9Dy7LO0rlyJLTOT7KM/S84JJ5D56U8n9CBVfj/dH31E18aNdG58h86NG+nesgX8fkBreq7998NXX0/H2nWojo7QsZKWhmviRFyTtaAIfSZN0i8DEUwwKhDA39SEr75ef+rqeubr6/Fay/76BpS3R6sVtxt7djb2nGxsWdnYcrKxZ+fo5dA0G3tODrYsvZ84nQS6ulBdXQQ6u1BdnQQ6uwh0dfZe19Wt13V26f07O1EeD+JygtOJzeVCnE79Cc27EFefqdNpzTtRfj/4/SifH+X36Xmvr2c+uN7n672PT9937DZEbHpqs4PdjtgEbHbEbjmVg9tsNr3OZtdmr4CORFKBgPbjBgKooDM5oPptQwW0xm4p7WILj04KfvT/co9jW69DxPqdrcglsM7Vc75eDm5rWfXZnr3sGNwHDiAwJAmMUAgrcRGkINOVkKYQ6OrCv3dvH/NRWFhqEkKBzj040zvxdvW8OVdv3kNBWSYiwvZ3GlmwbFLMU7St/h8AWUV74aWb4Hgr90ApePIy7Rw9958h+3rmQQdRcOEFNN1zL1lLl5K9dGliQ92wgV0/+AEZ8+dTduMve96Cc0p1qOq9x8Ldx+jigAX7wTn/gJK5oeOzPnMYE/74B3Z+7ets/8p5TLz/vt73MAZKKVqf+w+1N9yAr66OvC+dQfF3v4s9N3KHWLHZcC9ciHvhQsZfdSUdb7xB84oVtD73H5qffAp7bi7Zy5aR87kTcC9ejNh1SKy3pobOjRtDQqDr3XcJWA96W04OGXPnknXkUjLmzSNj7lwc48b1GqOvrh5PZSWe7ZV4Krfjqayk+6NttL70MoQ9xG1ZWVpATJxAoLOr58Hf0BASOOHYcnNxFhfhKCrSWkiRnkds+NtaCbS04m9tsaat+Bsa8XxcSaClRde/inDORBCXC8nIwJaeji09PTQvTqcWGs0t+LxelNeL8ngiTgeE04nY7frjcIDDoX8j63ciEEAF/OAPaGGhlJ4GAtY2vZ5AnNDf0BcV/YJjs/U83PvMI9LzAA8E9DM79IDv+ShrfL3Wh18nOA2Fv+plibHdNXWKEQopp61evz2Hke9OLKvZZ2UHO8aX9KwMZTXX9JhQEmHH6zjdProatYPb7wuwa+teZh5WhivNzlvP7aCr3Ut6ZvS32rbVq3CUleI68hh44y4dd188Q0fnfLgSjr9Fx9aHUfTtb9O++n/s/slPyXjqybjmHE9VFTu/9nUcxcVU3HF7fzNMyRz44v26B8O8L8EJv4qYA5F5yCFM/NNd7Lz4Eraf+xUm3X9fbzNclGvXXHcd7S+vIm3GDCp+eysZ8xMvfyEOB5mf/jSZn/406uqraXvlFVpWPEPLv//N3r//Hfu4caTPnkX3+x+EflucTtJnziT31FPJmDeX9HnzcE2aFNMcJCI4xxfjHF9M5sEH9dqmfD68u3drgWEJC09lJZ3vvofN7dYP++nTQw97hyUAHEXFOIrGRU8CTAClFKqzE39ra0hIBFpbCXR3Y8twY8tIR9IzsGVYD/70nmlQWA7m2gSFhterM929Xv1m73CEpmK36we/wzHkJreggFCWgJDwh3zwwW8Y40JBKUtT6G2+KMh0sbWuLe7h3hodZeQcH6YRDDSBbcerOLPBv62VQGcntdXd+DwBKg7IJyPbxbqV29n5fhPTFkV+o1YeDx2vrtGhqJ/9Brz3BKz8ke6/8OyPdfjnQRf1O87mclF2yy1Unn46u6++morbbov6z+FvbmbnxZeg/H4m/PEP0QXItGPgyp1xI37cixYx8d572LH8Iraf+xUmPnA/ror+tZeUx0PjfffTcOediM3G+Ct/RP6Xv5xcXaU+iMtF9pFHkn3kkQQ6O2l7eRUtK1bQve0j3AcfrDWAeXNJmzkTm2vospjF4cA1YQKuCRPg8MgRYqlCRBC3W/uPEtTMhvLauFyhKL2RyIqVoAYwAtf+JDG2hUJXs04YyuwtFPITrJTqqw1qCmH/YFnFgCRf6mL7GpzlE2BDE95du6jebAOBsml5uDIcpGc6qXynIapQ6Hjr7Z5Q1MxxOrt55Y/gwZN1VvApt0cNH0w/YDpF3/0udTffTPPj/yTvtC/020d5PFR969t4du5k4j13kzZ1auzvk2AIaMb8+Uy87z52LF/O9nPOZdID9+Oa1GMma3/jDWqu/Tmejz4ie9kyxl91Jc6SkhhnTB5bRgY5xx1LznHHDul5DYZPImM7JDXYpCazt+2/0Gq0EwjEDgP11WpNoZf5yO7UQiYZTcHbCbvexnmANjd5d+2i6oM9FE3IJj3Tic0mTJxTwI53m6KOqW31KnA6e0JRFy+HcQfoPIxT7+ynDfWl4PzzcB90ELW/+AWeqqpe25RS7L76Z3S8/jpl119H5kEHRTnLwMiYO4dJ99+H6u5m+znn0r1tG76mJp3w9pXzUF1dVPzhTip+99shFwgGg6E3Y1wo9K57FCTf7SKgoKUrtnPMW1uHLTOzf/hnMCw1UarWQsCLc84SADp37KLm42YqDuiJHpo8Zxxd7V5qP45c1rt91Srcixb2jMXuhLMe1hFB+x8ddwhis1H2yxvAZmPXD3+ko0IsGu68k+YnnmDcZZeRe8opiX+vJEifOZOJD9yPUort55zLtuNPoPnppym8+GKmPv2vhJ3gBoNhcIxtoRCqe9Q/+giIa0LSzXUivLnmlCWnKexYAwiOTx0NTic125oJ+BTlM3qEwoRZBYhN2P5O/xacuirqVrIOX9J7Q+F+un1ngjjLyym5+qd0rltH4z33AtD8r3/R8LvbyD3lZMZ947LEv9MASJ8+nUkPPoCkpZE2bRpT//k4xZd/t38hPYPBkDLGtk+hT4mLIOFCYWoMq4vOZo4QdppdouvuJMr2V6F4FpI1DmdJCVV1Ol+idL+eMMv0TCel++VS+U4jh3y+d8notlX9q6IOlJyTTqL1hRepv+02bG43dTfdhHvxYkquu25YojPSpk5l//8+n7JkL4PBEJux/Z/XVqcTTNy9C7IlrinU9vYnBMku1QLHn0Bstt8HVW+Gaug4y8qo68ymeHIOrvTeMnvS3EIaq9toberq/TWCoaj7DbK/ADpKpORnV+PIy6P2+utxVlRQcdvvhjQCJ+4YjEAwGEaMsf3f116nBUKfGv6h8tkxchWUz4evoQFHNE0BBW218cdQs1HXtA/2Yy6bSLO9iIoZ/et2T56jE6S2v9vYMw4rFDXr8CVD9ibvyM+n7Fe/wr14MRP++AfseXlDcl6DwTD6GdtCoa2+nz8BeorixSp14WtsBL8/cjRMdqmeJuJs3rFGTyfpyqh7c/ZDiY2yqf0TvvJL3eSMS+/lV+h46y0dinrEkn77D4bMgw9i0kMP4po4cUjPazAYRjdjWyhESFwDyHDZyXDaY5a66NVcpy8hoZCAs3n7q5A3STungQaKsQW8FGb0r+opIkyaM46qD/bg8+jooLZVqxGnk8yDD45/LYPBYIjD2BYKbXURNQXQfoVY5bNDvZlLYgmFOJqCUrDjtZCWAFDXmkFu8zZUXeRjJ80txOcNULVZF8trX72KjEULsWXGLpZnMBgMiTC2hUJ7Q79s5iD5mc6YPgWf1XHNEalcgLtQ92ho2RX7+g1bdDVVy5/Q1ealqSlA/p7NeKurIx5SPj0Ph8vG9ncbo4eiGgwGwwAZuyGpnnbwtkfN9M13u2L7FOpqEacTe36ERs42G2QlkMC241U9tTSF6g/1239+8xa8uyKXkXA47VTMKKDynQbmim6QMxShqAaDwQBjWVOIkrgWpDBO+WxvTS2O4uLo4ZPZJfF9CtvXaE2lcH8AqjbvwZlmJ9/dHVVTAJg8t5C2pm52r1qPs6xsSEJRDQaDAcayUIiSuBYkP45Q8NXWRs5mDpJIr+Ydr8LEQ0KF6qo376FsWh5p5SV4q6ObniZZoalV2z1kLjnclPw1GAxDxtgVCiFNIbL5qMDtorXbh8cXuTmHt7YmcjZzkHi9mpurYe8OmKhNR+17u9lT00H5Afk4y8rw7IquKWTlp1FQINRnTydrifEnGAyGoWPsCoVQMbwo0UdZ0RPYlFL4ausih6MGyS6Brr26AmokQvkJ2skcjCaqOCAfZ3k5vto6lM8X9fTFgV00507FPmdB9DEYDAZDkoxdodBmmY9iaAoQudRFoLkZ1dWFI1I4apB4uQrbXwVXFozXbSqrN+8hze1gXEWW7kDm9+OtiZ4Rnf/hSyA2qiq7ou5jMBgMyTJ2hUJ7nW5sb4/c3jJU6iKCUPBazXVi9hUOdWCL4lfYsQYmHAR2HQBWtXkP5dPzEZvgKi/X14liQvLu2kX6pldIc/ipfKcx4j4Gg8EwEFIqFETkOBHZLCJbReRHEbZfLiKbRGSjiPxXRGJ3ph9KYiSuQU9RvEhhqRGb6/Ql1Ks5gqbQ0QR1m0L+hJaGTlobuyi3+icEexV7d0V2NretWo2gmHhADjveayTgT7ApucFgMMQhZUJBROzA7cDxwCzgLBGZ1We3t4FFSql5wD+Am1M1nn7ESFyDHqEQyacQymaO6WiOoSnsfF1PI/gTAByl2vQULSy1bfVqnGVlTDl0Mt0dPmqiNN4xGAyGZEmlpnAQsFUptU0p5QEeAXq17VJKvaiU6rAWXwP6d21PFVHqHgXJy9BmpUg+BV9NLYjgKIrRbCE9FxwZkbOat78KNieULwS0PyEjx0V+qRsAW1oajqKiiJqC8njoWLOGzCWHM3F2IbYojXcMBoNhIKRSKJQDO8OWq6x10fgq8EykDSJysYisFZG19fX1QzO6KBVSgzjsNnIznBF9Cr66WuyFhUisHgMi0dty7lgD5QvAmYFSiqrNe6g4IL9XvoGzrCxirkKoKuqSJaRlOCidlmv8CgaDYchIpVCIlFEVseu8iJwDLAJuibRdKXWXUmqRUmpRUay380TxdkF3c9xm9oWZkUtdeGtrYzuZg2RHSGDzdMCut0P1jvbWdtDR7OnVjxl0a8xImkLby6t6VUWdNGccTbvaaWmMEvo6QFRA0bSrHaUi/mQGg2EfJZVCoQqYELZcAfR7yonI0cCPgZOVUt0pHE8PwWzmGJoCWFnNEXwKuuNaAkIhJ0ICW/VaCPhC9Y6qPtD+hPJ+QqEM7+7dqEBvJ3Lb6lW4Fy8KVUWdPFd3jds+xNrCa09+xMM/f52nb9vA3tqO+AcYDIZ9glQKhTeBaSIyRURcwJnAU+E7iMiBwB/RAqEuhWPpTZzEtSD57sjls321tZFLZvclqCmEv21vXwMITNBv+tWb95BdkE7OuPRehzrLy8HrxRdmLvPu2oVn60dkhlVFzRaLarYAACAASURBVBvvJqcoY0hNSLu37uWt53ZQul8uNduaefi613ntiY/wdvuH7BrR8Hb72fTKLirfaSAQMFqKwTDcpKxKqlLKJyLfAJ4F7MC9Sqn3ROTnwFql1FNoc1EW8HfLnr5DKXVyqsYUoi0xTaEg08m71c291gW6uvA3N8fOZg6SXaIrsXa3aMcz6HpH42dDRh4qoKj6cA9TPlXUr35RKCy1ujpkqmpbtRroXRVVRJg8t5D3Vu3C2+3Hmda7tWiyeLp8PP/A+2QXpHPiNz+Ft9vPmsc/Yt3K7Wx+vYbPfHEaUw/sP97B0t3p492Xq1j//E662rQgzipIY/Znypl5WCmZuWlDej2DwRCZlJbOVkqtAFb0WXd12PzRqbx+VEKaQmyfQkFmGk3tHpRSoYdgqONaopoCaG0hPRf8Ptj5Jsw/G4CG6ja62339/AkQLhR2wQJdyiIYiuqa2rus9uQ549j4QhVVm/cwZd640PpAwI+3q4s0d+INeF59/CNaGjo59fIDcaU7cKU7OPqCWcw6vIxVD3/IyrveZcLMfA7/0nTySwbf2KerzcuGF3ay8cUqPJ0+Js0pZMGxk+hs8/Deqmpef2obbz79MVPmj2P24eXaIW9LTQHAzlYPrU1d5BRmkJ4VOanRYNjXGZv9FOIUwwtSkOnE4w/Q7vGTlaZvVbD0RGKO5mCuwm4oOgBqNmjNwcpPqLbyE8qnxxAKlrM54PHQvmYNuaec3O8tvWxaHo40O9vfaeglFP59683s2PQO5974W3LGxXfQb3+vkfdWVTP/6AmUTes9prL98zjjqkW8u6qa15/cxiPXvcH8oyew8PjJuNKT/zNqb+5mw/M7eWdVNb5uP1MPLGLR8ZMpmpgd2me/A4vZW9fBptW7eP/V3Xz0Vj25RRnMPrycGZ8uISMrRvRXAnR3eNm1ZS9Vm/dQvXkPjdU9LVDT3A5yizLILXaTW5xBXnBa5B60wFBK4e3209HiobPFQ2erl842D3aHDVeGgzS3/rjSe6ZDLQhVQOHzBvB5/AQCCrvdhs0h2B02bHYxlXfHMGNTKLQ3gCsbnBkxd8t395S6CAoFX12w41qMbOYgfdtybreK4FmZzFWb95A33k1Wfn/TiM3txp6fH0pg61y3DtXREbHLmt1pY+LMAra/2xjSaraufZ0PX38FgJV3/IYv/uT66L0fgK52Ly8++D75pZkcfErkBj82u415R05g/4XjWfP4Vt56dgcfvlHLYadPY78FiZmUWpu6ePu5HWx6ZRcBX4Bpi8ez4LhJFJZlRdw/r9jNp0/bn4NOnsK2t+t5d1U1rz6+ldee+oj9FxQze0k5pfvlJnRtb7ef3Vt7hED9jlaU0vevdL9cDvn8ePLHZ9LS2ElzXSfN9R3UbGtmy9raXnFz4QIjr7hHcGTnp9Pd4aOz1UNHqyfsoe+ho9Xba9nnTSILXdBaW4adtAynFhQZDtIyHLjcegrg8/jxevSD3tfdM+/t9lvb/Pg8AXzd/rjXtzt6hITdLtidNmx2m14OCg9LgKiAIuC3PgHVsxyaD0Te5ldgE0TAZtOCSGxY097zNms/CdsPQAUAFEppQaeUFrpY097LPeus2woSnOrz62W9Ui9Lz37W35hIbzdh8Px6JPTaqBRh23r7yCQYoBn2p9vvzzhshQAHnTyF6YsTePYMgjEqFGInrgUJL3UxoUAnlnlr9AM+ZjZzkHBNAXTSWv5kyCkl4A+wa8teph8U/QcOD0ttW7Vah6IecnDEfSfNLWTb+noaq9vIHefkhfv+QGHFRA487kSev/sO1j79TxaffFrUa6165EM6W7187v/bO+/wuKoz4f/eqdKMumzLVnWXkbtxxQWwHWrWlEDwLnmAJQkbEiCwS7JkU75ks4Uk5MuXB7JkCTWE0GsSuunVxsbGvchVtiXZ6nU05Xx/3DujkTwj3Tu2ADPn9zz3uefeOe89Z+aeOe+p7/ud6bjcA89L+HI8LLuqiqrFJbz1yHZe+sMmSifls2Rl8iGlliOdrHtxH9s+MH6/yvkjmXV2BXkjfAOmFcXldjJx7kgmzh1Jw6F2Nr99iO3vH2bH6joKiv1MXlxC5bwivL7eVnwoGKZud2tMCdTtbSUSVjicQtGYHGafN5rSSfkUjc7F6U6uMMPBCK0NXTTXd9FS30lLfRfN9YkVRn/EIWRmucnM9uDLcZNXlIsv22Nee8jM8ZjXbsIhRU9XiEBnkJ6uMIGuIIHOEIGuED2dIeOzrhCBzhBtjd00dPXeA3B5nLg9DlweZ5+wL8eDy+PAbd53eZ2xa7fXiTiEcChCJKQIhyKEwxEioQjhoCIcjvT9LBQhbIYj4QihHmVU2A7B5XEYYadRiTvMc+zaIYjT0Rt2GD+diihUxKywI3EVvBmORJT5Wd84YNSZEl+hx8Iy4GcAKLOaVnGVeSzcV4H0iauIUyD0ZiQ+KH0r/d540ef1OfW7d+yH0Xwcb+/YCumpFAaxexSlIIFRvFBdPY6srNiS0AHx+MGba/QUIhFj09rEcwCo399GsDuccD4hiru4mMDOnUaWo0tRfYkr0YopxtLUvRsb6Gx8nbajR7jsZ7+gpLKKfZ+s551HHqR86gyKxhzrpW3nR3XsXFPHvBVj+gzfDMaocblc+oM5sbH/R36+munLyph9Xu+QUuOhDta+uJeda+pwOB1MXlzCzLPKyS7IGOTpySkszmLJZRNZcOE4dn5Ux+a3DvL2ozt4/6ldTJhTRM6wTA7uaOJwdQvhYAQRGF6Rw4zlZZRU5jNqXJ6tCXmn20H+SH9ChRevMNobu/H6XX0q/Qy/e8jmQKKoiOrTktVojof0VAodR2DYhEGjRZVCYx+lUGttkjlK9kjD1MXRHdDVmGA+IS+pqLukhPY336Sn5iA9u6rJ+8olSeP6c70ML89m5+rN1G5/hilnnkXppMkAfOma6zj8vet4/vbb+Np//wa3t7dC7mgJ8ObD2xkxOodZZ9u3R+hwCFPPKGXcrBF88Ew1H79sDCnNPreCmm1NVK8/gsvjZMbycqYvLzuhq4jcXidVC4upWlhM/b5WNr99iB2rawn1RCgsyWLK4hJKJuUb3uwyh6aoD6QwPi2GWulo0ov0VArt9VCxcNBo+QmM4gXr6nFbWY4aJWrqYv97xnV576a1wpIsMrOTdwfdxcWoQICWZ58BIOv0gb2sVUzJ591H/4A308eSy6+K3c/MzuHsb9/Ek//5Y9566D6WXX0tYHSPX39wG6GeCMuvOgWHM/VtK74cD0uvOIWqRcW89cgO3nx4B55MF7PPHc30pWVDvppnREUOIypyWHjJeCIhpVcPaTQpkn5KIRw0WuyDbFwDyPa6cDulj6mLUG0t3kWLrKeXUwx73zUmmf3DoXAc4WCEw9UtTFk8kCkocwMb0Pz4E7hLSvCMGTNg/GD3ZlToEJWnXUVmdk6fz0ZPm8mp51/A2r89y5iZsxk7cw5b3jnEvk0NLL5swglr6Y4cm8slt8ymtrqFwtKsIWuhJyOVlVAajaaX9HOy02FaFB1kOSoYY7T5Pk9sTkGFQoSOHsVlZZI5SvZIY6J5//uGvSMRavcYY90lk5LPJ4Bh6gIMReRfsnjAMePO1hbWv/gwLm8ZoUhlwjiLVl7JsPLRvHTnb6ndXcu7T+yipDKfqaefWOO0DocM6ZCNRqMZOtJQKVgzcRGlwO+JzSmEjh6FSAS3leWoUbJHQSQILQd67R1tb0LE2F8wENG9CkDCpajxvPWnewl2dzN+3mUc2NJEOIHjHZfHw/nX30ygs4Onbv0loFh25Sl6TFqj0cRIP6Vg0cRFlHxfr1G82G5muz2FKOW9k8zDK3IGbUk7s7Nx5OQMuBQV4MDmT9j85irmrLiYSQuq6OkKUVvdkjDusPLRjJ11AV0tO6ioOnxcq4A0Gs0Xj/RTCjELqcMGjmdSkNVrPjvmcW2kzZ4CGJvlRk4lGDDWzg+0FDUe7/jx+E87LelS1FAwyCt3/w+5I4qYd9FXKTulAIdTklpNbTjYTk11Of78CWx753EaDh5IGE+j0aQnaagUbA4fxc0phGqjPQU7q49MpVA2FxxODu9qJhJRlFQOPHQUpfSO2yn+VXIvpR/95SmaDtWw9Opv4fZm4Ml0UTwhj70JvLGFQxFevX8LGT43X/m3f8XlzeD5228jHDrWEqxGo0lP0k8ptNcbbjI9ic0q9Cff76G5K0g4ogjV1yFuN858a618ALKKIDMfJpwFGPMJDqcwapw1peAqKMCZk5Pws+baw3z41KNMnLeQsTPnxO6PnjqMptpOWo70dbzz0fN7OXqgnTMun8Tw8pGc9U/XU7+nmncfe8j69xkEFYmw6fVXuO+fr+XJ//oJHz7zOId2bCMcCp2wNAaivamR1iP1x/ih0Gg01ki/5SEdRwwTFxZ3fxb6PSgFzZ09BE3nOrZ2jro8cONGcBtLPg9ub6JoTM5xm7hWSrHqvt/jcDk546pv9vmsYkoh7zy+k32bjjLtTMPPUe2eFta+uI9JC0Yydoax8mrCnAVMXXY2a557kjHTZ1E2edpx5enIvj28es+dHNq+hRFjxtHWcJR3Hn4AALc3g+LKUyirmkrZ5KkUjZ2A03V8xa+juYm63buord5J3e6d1O2ppqOpMZZeYWkZhWUVDCstN85lFWQVFJ7wnb/B7m7aGhvoam3BnZGB1+fD4/Pj9flwOI7vPWs0nzbppxQsmriIEr+BzV1XZ283cxSvYToi0BnkyP42Zp832v4z+rHjg3fYu34tZ151DdkFfedH8op85BX52LexgWlnlhHsCbPq/q348zws+urEPnHPvOKb1GzZyAu/+w1X/PJ2MrKs9aDiCXR28v4TD7Huhb+Q4c/i7GtvZPKSpYjDQWdLMzVbN3Fgy0ZqtmzinUf+CIDL66WksoqyqqmUVk1l5LjxOF3JN5x1tjRTu3sndbt3xY72RnPeRISC4lLKp0xn5NjxON0eGmr201Czjz0ff8TmN16NPcfr81NQWsawqLIorWBYeQW+3LxjlIVSiq62VtobG2JHW+PR3nDDUdqbGgh0dJAMd0YmXp8Pr89vHr0Ko/eeH4/Ph4pECAYChHoChAIBQsGevtc9PQR7jHMoet+MA+B0u3E6XcbZFT33hh0uNy6XC4fL3SeOw+EgHAqhImHCoTAqEiYS7j3C4TAqHCYSCRMJhYhEIuY5TCQSQcSBw+FAoodI3LWz37VxRK8NY3IKFYmYdpAivfeUQqmIYcaDqP0jowcY6wlKnEVXMyym8SERiTNwJ3HG7YzPoulAvO2j6HWckTsVTf/Y9yv9jBr15iUuDtI/cvTBvcH+BrT6G9wzmb78HEbPOPXYjJxA0k8pdByBvHLL0Qt8UVMXQQrq6sicMiXlpA/tbEYpKB1kf8JgBDo7eP2BPzBizDhmnHV+wjgVUwrZ+GYNPd0hPnh2N811nVxw44xjVjy5MzI477qbefgn3+PVe/6H82/4nuWWtFKK7e+9xRsP3kNHcxPTlp3Nor+/ksysXvtJvtw8Js5fxMT5xoa/ztYWQ0ls3kjN1r5Konii0ZMoPWUywUDArPx3Ure7mrYGc4GACPmjSiirmkrR2PEUjR3PiNFj8WQmN6zX2dpiKIkD+zl6YB8NNfvZufp9Nq56KRYnIyubwtJy/PkFdDY3xSr/cLDffIsI/rx8sgsKyR9VTNnkqWQVDCO7oBBfTi7BngCBzk4CHR0EOjvo6eroc93Z2kpz3WG6Ozro6ewYcFhNHA7cXi8ujxeXxxM7u71ePD4fvry82D2ASChEOBgkHDbPoRChnh4CnZ2EQ8Z1OBgkEgoSCoWIhILGdSSC0+nC4XLicDgRpxOn0zg7nE4cThcOh6M37OwNO10Sq6wj4RCRYMSo2CNmRW8qjmgco9KPGPdMhWIsiZZeC6iGOdRYOHbP4UDMePGW5WJWUHtNoZoKRcXKaa9Bu14FE1Ue0co6XmH0vu44hWOe4/4AZg76XSfQHrF7MWt6/Z6f6P3Hf2aGA11D7xo3/ZRCez2UWNe0+X6j9drYHiCnrg7X8tT9AtVsb8LldlA0OjflZwC8++if6Ghu4sKbf4TDmXh4omJqIRtWHeCDZ3ez8fUapi0tpXRSQcK4I8dPZMEl/8C7jz7I2FlzqFp85qB5aDxUw6p77mT/pg0UjR3PhTf/iJHjJw4q58vJZeK8hUycZ5gZ6Wxt4eDWzWZPYiPvPvpgn/j5o0oomVQVpwDG4U2yEmugNH1VUymrmhq7p5Sis6U5piQaDuznaM1+6vfsIiu/kFHjK8kqKCS7oJCswmFk5ReSVVCIPy//uIe94jEqbUNxOF3OPgrgRKaj0VglvUpdJAydRy2vPAIo9BsG3FrqG1CBgDWT2Uk4uL2JUeMHNtM8GLXVO1n/0t+YcdZ5A1bCxePzcGc42fh6DXlFPhZceKx11HjmXngJezesZdU9d1JSWUVuEvtOwUA3Hz79GGueewq318uyq69l2pfOSXns3JeTy4R5pzFhnrGxr6utlYPbt+LNzGTEmHG2vMbZQcwWvz8vn4qpM4YkDSsYCsCDP+/4eo8azYkivVYfdTYaXjlszCnkmfb5uw4ZPhFsLUeNT7q1h4aDHZRY3J+QiEgkzKt3/w5fbi6LVl4xYFyny0F5VQHiEJZfVYXLM3Cl7XA4Ofc7/wLA83f8mkg4fEycXR99yP3/8m0+fPoxJi1cwj/+5vfMOPv8EzqZmpmdw/jZ8yibPG3IFIJGo0lOevUUbG5cA8hwO/F7nDHnOqkqhYM7DFPZpZWJh3CssP6l56nbvYvzv/t9SxXmwksmMO3MUorGJF7S2p/cEUUs+/q1vHDHr1n9zOPM/8pKAFrqa3ntvv9l97o1FJaWc9n/uZXSqtTnVjQazeeXNFMK9jauRcn3e4jUG7KWfDMn4OD2JjwZToaX21/dA9De2MC7j/6RimkzqVyw2JJMdkGGbTMWpyw6g93r1vDeE3+mtGoKNVs28eHTjyFOJ6d/7WpmnrtCj3VrNF9g0uvfbdPuUZRCvwfH9iMggmv44NZVE1GzvYniCXkp+yx4/Y93Ew6FWPb1a4fUw5aIsPwb3+bQjq08+tNbAJg4fxFnXPENsgut97A0Gs3JSXophVhPwV7Fnu/34GluwDmsEHHbd97S1thNS30XU5YM7D8hGXvWr2XH+29z2lcvJ39k8eACx0mGP4vzb/g+7z3+EHP+7mJGT5815GlqNJrPB+mlFNrrwemBDGsmJqIU+Dz4WhrsmcyOIzafkML+hGBPgFX33kl+cSlzViR3x3miKak8hUt/9B+fWnoajebzQXophY4jhnMdm8MvBX4P2e3NuCZVpZTswW1NZPjdFBb3nU9QShHqCcQ2NnV3dBDobDeuzXuHd22npa6WS3/8X7hS6KVoNBqNHdJLKbTXW/K41p98v4eCzmZkuLW5iO6uAJ98XE315l0c3XeAUM0hQq4A/33jkzh6upBgN45gN9LTjahjl37GoxwuOicu4o5tCrZtsJ13jUbzxeGimaUsGFc4pGmkl1LoqDesltqk0BUhO9hFML/3ZYTDYXZXH2bLuh0c3r2PjtpDhNuO4Ag0IuFmIM5Kp3jpdvvpCnvpcXoJurIIer0EneZ17Mjod+0h4jBf0c5jTWFrNJr0Yv7YoVUIkHZK4SgUTR08Xj98dXvYXjqZ3RurCXzzRpxdrThDzaB64mI5EVcOoYxcIlmlOIfl4i/NZ/ikYeSO9ONMcdVRH+Wi0WjSmkkFCazynWDSRykoZc4pWFtW2dHazAu338nhLdX0hI5CYQjaduB2+Al6Mmjz59HmD3E0u526/CbqcxuIOBK8sJ3modFoNMfJj+f/mLLKsiFNI32UQnczhHsG3LgWCoV47YH72fXWR3QHjqJUN+Alw1vE8JmTGbtiMTkF2vSCRqP5bBiRmbrtNaukj1IYYOPaupdfZM3jf6OzrYGIagWcuJ3DKKqczvk3XE9Wvr0lrBqNRnOyMqRKQUTOAX4LOIG7lVK39vvcC/wROBVoAC5TSu0dksz027i2b8smXvvfB2ipryMcMbx1OR2FFI6azPJr/5HSiZOGJBsajUbzeWbIlIKIOIHfAV8CaoA1IvKcUmpLXLSvA01KqfEishL4BXDZkGSovZ62Hjd/vesljux/kGD4CBDBITlk545j7sq/Y8bS1H0laDQazReBoewpzAV2KaV2A4jII8AFQLxSuAD4qRl+ArhDREQlcl10nNz/u1doaFwCbEQkk8yMEiYtnceSy7+GSxt402g0GmBolUIJcCDuugaYlyyOUiokIi1AIdBnUb6IXANcA1Bebt2VZjw5w7NpawlQMm0S5173bTJT8EWs0Wg0X3SGUikksiXRvwdgJQ5KqbuAuwBmz56dUi/i4n//ZSpiGo1Gk1YMpee1GiB+QW0pcChZHBFxAblA4xDmSaPRaDQDMJRKYQ0wQUTGiIgHWAk81y/Oc8CVZvgS4LWhmE/QaDQajTWGbPjInCO4DngJY0nqvUqpzSLy78BHSqnngHuAB0VkF0YPYeVQ5Uej0Wg0gzOky26UUs8Dz/e795O4cDdw6VDmQaPRaDTWGcrhI41Go9GcZGiloNFoNJoYWiloNBqNJoZWChqNRqOJISfbClAROQLsS1F8GP12S3+K8ukm+1mmfTLKfpZp6+98csger3yFUmpwf8RKqbQ5MJbCfiby6SZ7suZb/176O39eZU+EvJVDDx9pNBqNJoZWChqNRqOJkW5K4a7PUD7dZD/LtE9G2c8ybf2dTw7ZEyE/KCfdRLNGo9Foho506yloNBqNZgC0UtBoNBpNjLRRCiJyjohsF5FdInKLDbl7RaReRDalkGaZiLwuIltFZLOIfNemfIaIrBaRDab8z2zKO0XkYxH5q72cg4jsFZGNIrJeRD6yKZsnIk+IyDbzuy+wKFdpphc9WkXkRptp32T+VptE5GERybAh+11TbvNg6SYqFyJSICKviMhO85xvQ/ZSM92IiMxOIe1fmb/3JyLytIjk2ZD9uSm3XkReFpFiq7Jxn90sIkpEhtlI96cicjDufZ9n5zub9683/9ebRSShJ60kaT8al+5eEVlvQ3aGiHwQ/W+IyFwbstNF5H3zv/UXEclJIpuw7rBaxo6LoV7z+nk4MEx3VwNjAQ+wAaiyKLsEmAVsSiHdUcAsM5wN7LCarikjQJYZdgMfAvNtyP8z8GfgrynkfS8wLMXf+wHgG2bYA+Sl+M5qMTbcWJUpAfYAmeb1Y8BVFmWnAJsAH4b14FeBCXbKBfBL4BYzfAvwCxuypwCVwBvAbLtlEjgLcJnhX9hMOycufAPwe6uy5v0yDBP5+5KVmSTp/hS42eL7SSR/pvmevOb1CDv5jvv818BPbKT7MnCuGT4PeMOG7BrgdDN8NfDzJLIJ6w6rZex4jnTpKcwFdimldiuleoBHgAusCCql3iJFb3BKqcNKqXVmuA3YilFxWZVXSql289JtHpZWBohIKXA+cLetTB8nZstnCYavDJRSPUqp5hQetQyoVkrZ3b3uAjLF8OTn41hvf8k4BfhAKdWplAoBbwIXJYucpFxcgKEQMc8XWpVVSm1VSm23ktEk8i+b+Qb4AMPToVXZ1rhLP0nK2AD/hd8A308mN4isJZLIXwvcqpQKmHHq7aYtIgJ8FXjYhqwCoi38XJKUsSSylcBbZvgV4CtJZJPVHZbK2PGQLkqhBDgQd12Djcr5RCAio4GZGK19O3JOs2tbD7yilLIq//8w/qgRO+nFoYCXRWStiFxjQ24scAS4zxy6ultE/Cmkv5Ikf9RkKKUOArcB+4HDQItS6mWL4puAJSJSKCI+jBZg2SAy/SlSSh0283IYGGFT/kRxNfCCHQER+U8ROQBcDvxksPhxciuAg0qpDfayGOM6c+jq3hSGQiYCi0XkQxF5U0TmpJD+YqBOKbXThsyNwK/M3+s24Ac2ZDcBK8zwpVgoY/3qjiEvY+miFCTBvU9tLa6IZAFPAjf2a5UNilIqrJSagdHymysiUyyk92WgXim1NqUMGyxUSs0CzgW+IyJLLMq5MLrMdyqlZgIdGN1cy4jhvnUF8LhNuXyMltQYoBjwi8jXrMgqpbZiDLu8AryIMcQYGlDoc4iI/BAj3w/ZkVNK/VApVWbKXWcxLR/wQ2wokX7cCYwDZmAo8V/blHcB+cB84HvAY2bL3w5/j83GB0YP5Sbz97oJs1dskasx/k9rMYaFegaKfDx1R6qki1Kooa9GLsX6sMJxISJujJf6kFLqqVSfYw7BvAGcYyH6QmCFiOzFGCpbKiJ/spneIfNcDzyNMQRnhRqgJq5H8wSGkrDDucA6pVSdTbnlwB6l1BGlVBB4CjjNqrBS6h6l1Cyl1BKMbr+d1iNAnYiMAjDPCYczhgoRuRL4MnC5MgedU+DPJBnSSMA4DAW8wSxrpcA6ERlpRVgpVWc2eiLAH7BexqLUAE+Zw6yrMXrFCSe6E2EOMV4MPGoz3SsxyhYYDRfL+VZKbVNKnaWUOhVDGVUPkL9EdceQl7F0UQprgAkiMsZsha4EnhvqRM1Wyz3AVqXU/01Bfnh0FYmIZGJUetsGk1NK/UApVaqUGo3xXV9TSllqMZtp+UUkOxrGmMS0tPpKKVULHBCRSvPWMmCL1bRNUmm9gTFsNF9EfOZvvwxjLNYSIjLCPJdjVBZ28/AcRoWBeX7WpnzKiMg5wL8CK5RSnTZlJ8RdrsBCGQNQSm1USo1QSo02y1oNxuRorcV0R8VdXoTFMhbHM8BS81kTMRY12LEguhzYppSqsZnuIeB0M7wUG42HuDLmAH4E/D5JvGR1x9CXsRM9c/15PTDGiHdgaOYf2pB7GKNrG8Qo9F+3IbsIY5jqE2C9eZxnQ34a8LEpv4kkKyQGecYZ2Fx9hDEvsME8Ntv5vUz5GcBHZr6fAfJtyPqABiA3xff8M4xKbRPw2HnYkwAAAwlJREFUIObKFIuyb2MosA3AMrvlAigEVmFUEquAAhuyF5nhAFAHvGQz7V0Y82bRcpZsBVEi2SfN3+sT4C9ASSr/BQZYsZYk3QeBjWa6zwGjbH5nD/AnM+/rgKV28g3cD3wrhfe8CFhrlpMPgVNtyH4Xox7aAdyKaVUigWzCusNqGTueQ5u50Gg0Gk2MdBk+0mg0Go0FtFLQaDQaTQytFDQajUYTQysFjUaj0cTQSkGj0Wg0MbRS0KQdItJunkeLyD+c4Gf/W7/r907k8zWaoUYrBU06MxqwpRRExDlIlD5KQSlleUe1RvN5QCsFTTpzK4ZBtfVi+GFwiuGXYI1ppO2fAETkDNO2/Z8xNlshIs+YxgI3Rw0GisitGBZa14vIQ+a9aK9EzGdvMm3pXxb37Dek1//EQ1H7PSJyq4hsMfNy26f+62jSEtdnnQGN5jPkFgx7/l8GMCv3FqXUHBHxAu+KSNTK6lxgilJqj3l9tVKq0TQ/skZEnlRK3SIi1ynDgGF/LsbY6T0dwz7PGhGJmlCeCUzGMJ/wLrBQRLZg7HKepJRSksRpjkZzotE9BY2ml7OAK0xT5R9imBSI2gVaHacQAG4QkQ0YvgvK4uIlYxHwsDIMwNVh+GuImnperZSqUYZhuPUYw1qtQDdwt4hcDNiyZ6TRpIpWChpNLwJcr5SaYR5jVK8/ho5YJJEzMIypLVBKTcewTzWY28+BTDoH4sJhDA9qIYzeyZMYjlRetPVNNJoU0UpBk860Ydi0j/IScK1pshgRmZjEQVAu0KSU6hSRSRj2/KMEo/L9eAu4zJy3GI7hnW51soyZdvRzlVLPYzh1STQkpdGccPScgiad+QQImcNA9wO/xRi6WWdO9h4hsbvDF4FvicgnwHaMIaQodwGfiMg6pdTlcfefBhZgWNZUwPeVUrWmUklENvCsiGRg9DJuSu0rajT20FZSNRqNRhNDDx9pNBqNJoZWChqNRqOJoZWCRqPRaGJopaDRaDSaGFopaDQajSaGVgoajUajiaGVgkaj0Whi/H8j9UM3xfCSSAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x321e370>"
      ]
     },
     "metadata": {},
     "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": 29,
   "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(int(ceil(log(ITERATIONS + 1, 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": 27,
   "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",
    "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": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.18 µs ± 86.6 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": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.37 µs ± 184 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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}