{
 "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": 23,
   "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": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1], [3], [0, 1], [1, 4], [1, 5], [1]]\n",
      "[0.039309, 0.353161, 0.027576, 0.322223, 0.162165, 0.095566]\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": 3,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAEWCAYAAACwtjr+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHjlJREFUeJzt3X24V2Wd7/H3J9Ry8gFM9BDQYEWND1OoO6Rj05gWolYwkzbanGDMGcow7WquU3hmJp9yDp0p7bIHG0oUHEcyy5EpjBgfc1Jko4gietghxQ6SrSBiJh7oe/5Y9x6X29/+7d9+uPfa/Pi8rmtdv7W+6173+i678uu91r3XUkRgZmaW02uqTsDMzJqfi42ZmWXnYmNmZtm52JiZWXYuNmZmlp2LjZmZZediYzaIJJ0gqb3qPHpL0nWSvpTWd8trsGq52NgeT9J6Sb+T9Lyk36R/se5XdV79IemfJX2rtL23pN92E5tUTZa2J3GxMSt8KCL2AyYARwMXVpxPf90D/GlpuwX4FfDeLjGAFYOVlO25XGzMSiLiN8ASiqIDgKTTJD0k6TlJGyRdXNo3TlJImiHpV5KelvR3pf37ppHSVkmPAe8qn0/S4ZLukvSspNWSPlzad52kb0m6LY26/lPSf5P0tdTf45KO7uZS7gYOl3Rw2v4TYCHw+i6x+yLi/6XzfT+N7LZJukfSkY38M5N0vqTHJI1ppL3tmVxszErSvzBPAdpK4d8C04HhwGnAuZKmdTn0PcDbgZOAL0o6PMUvAt6SlpOBGaVz7Q38O/BT4BDgM8ANkt5e6vejwN8DBwM7gPuAB9P2zcAVta4jItqBX1IUFChGND8Dft4ldk/psNuA8SmXB4EbavVdJukfgL8C/jSd06wmFxuzwr9J2g5sADZTFAkAIuKuiHgkIn4fEauAG3nlLSqASyLidxHxMPAw8M4U/yhweURsiYgNwFWlYyYB+wFzIuKliLgD+BFwVqnNLRGxIiJeBG4BXoyIBRGxC/gexS2/7twNvFfSa4CJwP0UBaczdnxq03md8yJie0TsAC4G3inpwG76lqQrKAro+yKio04eZi42Zsm0iNgfOAH4I4qRAwCSjpN0p6QOSduAT5X3J78prb9AUUQA3khRwDr9srT+RmBDRPy+y/7Rpe2nSuu/q7FdbyLDPRSjlz8G1kXEC8C9pdi+wLJ0jcMkzZH0C0nPAetTH12vs9NwYCbwvyNiW50czAAXG7NXiIi7geuAr5TC/wosAsZGxIHAtwE12OUmYGxp+02l9Y3A2DTKKO//dS/T7s49FCOs0yhGNACrUz6nAcvTiAngY8BU4P3AgcC4FO/uOrcCHwSulXT8AOVrTczFxuzVvgZ8QFLnJIH9gS0R8aKkiRT/Ym7UTcCFkkak50GfKe1bRvE86PNpGvIJwIcoHuT3W0S0UYyELiAVmyi+KbIsxcrPa/aneCb0DPAHwD820P9dwF8Ct0g6biBytublYmPWRXr+sAD4hxT6NHBpeqbzRYoC0qhLKG6NPUkxEeD60nleAj5MMSHhaeBbwPSIeLy/11ByDzAS+M9S7GcUkwDKxWZByvPXwGMUz3d6FBFLgbOBRZKOHYiErTnJH08zM7PcPLIxM7PsXGzMzCw7FxszM8vOxcbMzLLbq+oEhoqDDz44xo0bV3UaZma7lRUrVjwdESN7audik4wbN47W1taq0zAz261I+mXPrXwbzczMBoGLjZmZZediY2Zm2bnYmJlZdi42ZmaWnYuNmZll52JjZmbZudiYmVl22YqNpNdJekDSw5JWS7okxa+T9KSklWmZkOKSdJWkNkmrJB1T6muGpLVpmVGKHyvpkXTMVZKU4gdJWpraL5U0Itd1mplZz3K+QWAHcGJEPC9pb+BeSbelff8zIm7u0v4UYHxajgOuBo6TdBBwEdACBLBC0qKI2JrazKT40NNiYApwGzAbuD0i5kianba/kPFam8642T+uOoWGrJ9zWtUpmFkDso1sovB82tw7LfW+1DYVWJCOux8YLmkUcDKwNCK2pAKzFJiS9h0QEfelT90uAKaV+pqf1ueX4mZmVoGsz2wkDZO0EthMUTCWpV2Xp1tlV0p6bYqNBjaUDm9PsXrx9hpxgEMjYhNA+j2km/xmSmqV1NrR0dHn6zQzs/qyFpuI2BURE4AxwERJRwEXAn8EvAs4iJdvb6lWF32I9ya/uRHREhEtI0f2+NJSMzPro0GZjRYRzwJ3AVMiYlO6VbYDuBaYmJq1A2NLh40BNvYQH1MjDvBUus1G+t08oBdkZma9knM22khJw9P6vsD7gcdLRUAUz1IeTYcsAqanWWmTgG3pFtgSYLKkEWlW2WRgSdq3XdKk1Nd04NZSX52z1maU4mZmVoGcs9FGAfMlDaMoajdFxI8k3SFpJMVtsJXAp1L7xcCpQBvwAnA2QERskXQZsDy1uzQitqT1c4HrgH0pZqF1znabA9wk6RzgV8AZ2a7SzMx6lK3YRMQq4Oga8RO7aR/ArG72zQPm1Yi3AkfViD8DnNTLlM3MLBO/QcDMzLJzsTEzs+xcbMzMLDsXGzMzy87FxszMsnOxMTOz7FxszMwsOxcbMzPLzsXGzMyyc7ExM7Pscr4bzcwy8tdUbXfikY2ZmWXnYmNmZtm52JiZWXYuNmZmlp2LjZmZZediY2Zm2bnYmJlZdi42ZmaWnYuNmZlll63YSHqdpAckPSxptaRLUvwwScskrZX0PUn7pPhr03Zb2j+u1NeFKf6EpJNL8Skp1iZpdile8xxmZlaNnCObHcCJEfFOYAIwRdIk4MvAlRExHtgKnJPanwNsjYi3Alemdkg6AjgTOBKYAnxL0jBJw4BvAqcARwBnpbbUOYeZmVUgW7GJwvNpc++0BHAicHOKzwempfWpaZu0/yRJSvGFEbEjIp4E2oCJaWmLiHUR8RKwEJiajunuHGZmVoGsz2zSCGQlsBlYCvwCeDYidqYm7cDotD4a2ACQ9m8D3lCOdzmmu/gb6pyja34zJbVKau3o6OjPpZqZWR1Zi01E7IqICcAYipHI4bWapV91s2+g4rXymxsRLRHRMnLkyFpNzMxsAAzKbLSIeBa4C5gEDJfU+WmDMcDGtN4OjAVI+w8EtpTjXY7pLv50nXOYmVkFcs5GGylpeFrfF3g/sAa4Ezg9NZsB3JrWF6Vt0v47IiJS/Mw0W+0wYDzwALAcGJ9mnu1DMYlgUTqmu3OYmVkFcn48bRQwP80aew1wU0T8SNJjwEJJXwIeAq5J7a8BrpfURjGiORMgIlZLugl4DNgJzIqIXQCSzgOWAMOAeRGxOvX1hW7OYWZmFchWbCJiFXB0jfg6iuc3XeMvAmd009flwOU14ouBxY2ew8zMquE3CJiZWXYuNmZmlp2LjZmZZediY2Zm2bnYmJlZdi42ZmaWnYuNmZll52JjZmbZudiYmVl2LjZmZpadi42ZmWXnYmNmZtm52JiZWXYuNmZmlp2LjZmZZediY2Zm2bnYmJlZdi42ZmaWnYuNmZll52JjZmbZZSs2ksZKulPSGkmrJV2Q4hdL+rWklWk5tXTMhZLaJD0h6eRSfEqKtUmaXYofJmmZpLWSvidpnxR/bdpuS/vH5bpOMzPrWc6RzU7gbyPicGASMEvSEWnflRExIS2LAdK+M4EjgSnAtyQNkzQM+CZwCnAEcFapny+nvsYDW4FzUvwcYGtEvBW4MrUzM7OK9KrYSBoh6R2NtI2ITRHxYFrfDqwBRtc5ZCqwMCJ2RMSTQBswMS1tEbEuIl4CFgJTJQk4Ebg5HT8fmFbqa35avxk4KbU3M7MK9FhsJN0l6QBJBwEPA9dKuqI3J0m3sY4GlqXQeZJWSZonaUSKjQY2lA5rT7Hu4m8Ano2InV3ir+gr7d+W2nfNa6akVkmtHR0dvbkkMzPrhUZGNgdGxHPAnwPXRsSxwPsbPYGk/YAfAJ9N/VwNvAWYAGwCvtrZtMbh0Yd4vb5eGYiYGxEtEdEycuTIutdhZmZ910ix2UvSKOCjwI9607mkvSkKzQ0R8UOAiHgqInZFxO+B71DcJoNiZDK2dPgYYGOd+NPAcEl7dYm/oq+0/0BgS29yNzOzgdNIsbkEWELx3GS5pDcDa3s6KD0juQZYExFXlOKjSs3+DHg0rS8CzkwzyQ4DxgMPAMuB8Wnm2T4UkwgWRUQAdwKnp+NnALeW+pqR1k8H7kjtzcysAnv13IRNEfFfkwIiYl2Dz2yOBz4OPCJpZYr9L4rZZBMobmutBz6Z+l0t6SbgMYqZbLMiYheApPMoCt4wYF5ErE79fQFYKOlLwEMUxY30e72kNooRzZkN5GtmZpk0Umy+DhzTQOwVIuJeaj87WVznmMuBy2vEF9c6LiLW8fJtuHL8ReCMevmZmdng6bbYSHo38N+BkZI+V9p1AMUIw8zMrCH1Rjb7APulNvuX4s/x8nMSMzOzHnVbbCLibuBuSddFxC8HMSczM2sy9W6jfS0iPgt8Q1Ktv1H5cNbMzMysadS7jXZ9+v3KYCRiZmbNq95ttBXp9+7OWHq1zNiIWDUIuZmZWZMYlHejmZnZni37u9HMzMyyvhvNzMwMGis2l9KHd6OZmZl16vF1NRHxfeD7pe11wEdyJmVmZs2l3t/ZfD4i/o+kr1P7WzDnZ83MzMyaRr2RzZr02zoYiZiZWfOq93c2/55WV0XEQ4OUj5mZNaFGJghcIelxSZdJOjJ7RmZm1nR6LDYR8T7gBKADmCvpEUl/nzsxMzNrHo2MbIiI30TEVcCngJXAF7NmZWZmTaWR19UcLuliSY8C3wB+DozJnpmZmTWNRj4LfS1wIzA5IjZmzsfMzJpQI3/UOWkwEjEzs+bV0DObvpA0VtKdktZIWi3pghQ/SNJSSWvT74gUl6SrJLVJWiXpmFJfM1L7tZJmlOLHpgkLbelY1TuHmZlVI1uxAXYCfxsRhwOTgFmSjgBmA7dHxHjg9rQNcAowPi0zgauhKBzARcBxwETgolLxuDq17TxuSop3dw4zM6tAt8VG0vXp94K+dBwRmyLiwbS+neKNBKOBqcD81Gw+MC2tTwUWROF+YHh62/TJwNKI2BIRW4GlwJS074CIuC8iAljQpa9a5zAzswrUG9kcK+kPgU9IGpFuTf3X0puTSBoHHA0sAw6NiE1QFCTgkNRsNLChdFh7itWLt9eIU+ccZmZWgXoTBL4N/AR4M7ACUGlfpHiPJO0H/AD4bEQ8lx6r1GxaIxZ9iDdM0kyK23C86U1v6s2hZmbWC92ObCLiqvS8ZV5EvDkiDistjRaavSkKzQ0R8cMUfirdAiP9bk7xdmBs6fAxwMYe4mNqxOudo+s1zo2IlohoGTlyZCOXZGZmfdDI62rOlfROSeel5R2NdJxmhl0DrImIK0q7FgGdM8pmALeW4tPTrLRJwLZ0C2wJMDndyhsBTAaWpH3bJU1K55repa9a5zAzswo08gaB84EbKJ57HALcIOkzDfR9PPBx4ERJK9NyKjAH+ICktcAH0jbAYmAd0AZ8B/g0QERsAS4Dlqfl0hQDOBf4bjrmF8BtKd7dOczMrAKNvEHgr4HjIuK3AJK+DNwHfL3eQRFxL7WfqwCcVKN9ALO66WseMK9GvBU4qkb8mVrnMDOzajTydzYCdpW2d9F9ETEzM3uVRt+NtkzSLWl7GsWzGDMzs4Y08m60KyTdBbyHYkRztr/caWZmvdHIyIb0JoAHM+diZmZNKue70czMzAAXGzMzGwR1i42kYZL+Y7CSMTOz5lS32ETELuAFSQcOUj5mZtaEGpkg8CLwiKSlwG87gxFxfraszMysqTRSbH6cFjMzsz5p5O9s5kvaF3hTRDwxCDmZmVmTaeRFnB8CVlJ82wZJEyQtyp2YmZk1j0amPl8MTASeBYiIlcBhGXMyM7Mm00ix2RkR27rEevVFTDMz27M1MkHgUUkfA4ZJGg+cD/w8b1pmZtZMGhnZfAY4EtgB3Ag8B3w2Z1JmZtZcGpmN9gLwd+mjaRER2/OnZWZmzaTHYiPpXRRfydw/bW8DPhERKzLnZmZ7kHGzd48/51s/57SqU9gtNfLM5hrg0xHxMwBJ76H4oNo7ciZmZmbNo5FnNts7Cw1ARNwL+FaamZk1rNuRjaRj0uoDkv6ZYnJAAH8B3JU/NTMzaxb1RjZfTcsE4G3ARRR/4Hk48O6eOpY0T9JmSY+WYhdL+rWklWk5tbTvQkltkp6QdHIpPiXF2iTNLsUPk7RM0lpJ35O0T4q/Nm23pf3jGvxnYWZmmXQ7somI9/Wz7+uAbwALusSvjIivlAOSjgDOpJhi/UbgPyS9Le3+JvABoB1YLmlRRDwGfDn1tVDSt4FzgKvT79aIeKukM1O7v+jntZiZWT80MhttODAdGFdu39MnBiLinl6MKqYCCyNiB/CkpDaKV+QAtEXEupTLQmCqpDXAicDHUpv5FKOuq1NfF6f4zcA3JCki/NYDM7OKNDJBYDFFoXkEWFFa+uo8SavSbbYRKTYa2FBq055i3cXfADwbETu7xF/RV9q/LbV/FUkzJbVKau3o6OjHJZmZWT2NTH1+XUR8boDOdzVwGcVEg8songl9AlCNtkHtYhh12tPDvlcGI+YCcwFaWlo88jEzy6SRkc31kv5G0ihJB3UufTlZRDwVEbsi4vfAd3j5Vlk7MLbUdAywsU78aWC4pL26xF/RV9p/ILClL/mamdnAaKTYvAT8E3AfL99Ca+3LySSNKm3+GdA5U20RcGaaSXYYMB54AFgOjE8zz/ahmESwKD1/uRM4PR0/A7i11NeMtH46cIef15iZVauR22ifA94aEU/3pmNJNwInAAdLaqeYOn2CpAkUt7XWA58EiIjVkm4CHgN2ArMiYlfq5zxgCTAMmBcRq9MpvgAslPQl4CGKNx2Qfq9Pkwy2UBQoMzOrUCPFZjXwQm87joizaoSvqRHrbH85cHmN+GKKSQpd4+t4+TZcOf4icEavkjUzs6waKTa7gJWS7qT4zADQ89RnMzOzTo0Um39Li5mZWZ808j2b+YORiJmZNa9G3iDwJDX+TiUi3pwlIzMzazqN3EZrKa2/juLhe5/+zsbMzPZMPf6dTUQ8U1p+HRFfo3gvmZmZWUMauY12TGnzNRQjnf2zZWRmZk2nkdtoXy2t76T4Y8yPZsnGzMyaUiOz0fr7XRszM9vDNXIb7bXAR3j192wuzZeWmZk1k0Zuo91K8U2YFZTeIGBmZtaoRorNmIiYkj0TMzNrWo18YuDnkv44eyZmZta0GhnZvAf4q/QmgR0UX8KMiHhH1szMzKxpNFJsTsmehZmZNbVGpj7/cjASMTOz5tXIMxszM7N+cbExM7PsXGzMzCw7FxszM8suW7GRNE/SZkmPlmIHSVoqaW36HZHiknSVpDZJq8pvmpY0I7VfK2lGKX6spEfSMVdJUr1zmJlZdXKObK4Dur55YDZwe0SMB25P21BMrx6flpnA1VAUDuAi4DhgInBRqXhcndp2Hjelh3OYmVlFshWbiLgH2NIlPBWYn9bnA9NK8QVRuB8YLmkUcDKwNCK2RMRWYCkwJe07ICLui4gAFnTpq9Y5zMysIoP9zObQiNgEkH4PSfHRwIZSu/YUqxdvrxGvd45XkTRTUquk1o6Ojj5flJmZ1TdUJgioRiz6EO+ViJgbES0R0TJy5MjeHm5mZg0a7GLzVLoFRvrdnOLtwNhSuzHAxh7iY2rE653DzMwq0si70QbSImAGMCf93lqKnydpIcVkgG0RsUnSEuAfS5MCJgMXRsQWSdslTQKWAdOBr/dwDjOzQTVu9o+rTqEh6+eclv0c2YqNpBuBE4CDJbVTzCqbA9wk6RzgV8AZqfli4FSgDXgBOBsgFZXLgOWp3aUR0Tnp4FyKGW/7ArelhTrnMDOzimQrNhFxVje7TqrRNoBZ3fQzD5hXI94KHFUj/kytc5iZWXWGygQBMzNrYi42ZmaWnYuNmZll52JjZmbZudiYmVl2LjZmZpadi42ZmWXnYmNmZtm52JiZWXYuNmZmlp2LjZmZZediY2Zm2bnYmJlZdi42ZmaWnYuNmZll52JjZmbZudiYmVl2LjZmZpadi42ZmWXnYmNmZtlVUmwkrZf0iKSVklpT7CBJSyWtTb8jUlySrpLUJmmVpGNK/cxI7ddKmlGKH5v6b0vHavCv0szMOlU5snlfREyIiJa0PRu4PSLGA7enbYBTgPFpmQlcDUVxAi4CjgMmAhd1FqjUZmbpuCn5L8fMzLozlG6jTQXmp/X5wLRSfEEU7geGSxoFnAwsjYgtEbEVWApMSfsOiIj7IiKABaW+zMysAlUVmwB+KmmFpJkpdmhEbAJIv4ek+GhgQ+nY9hSrF2+vEX8VSTMltUpq7ejo6OclmZlZd/aq6LzHR8RGSYcASyU9Xqdtrect0Yf4q4MRc4G5AC0tLTXbmJlZ/1UysomIjel3M3ALxTOXp9ItMNLv5tS8HRhbOnwMsLGH+JgacTMzq8igFxtJr5e0f+c6MBl4FFgEdM4omwHcmtYXAdPTrLRJwLZ0m20JMFnSiDQxYDKwJO3bLmlSmoU2vdSXmZlVoIrbaIcCt6TZyHsB/xoRP5G0HLhJ0jnAr4AzUvvFwKlAG/ACcDZARGyRdBmwPLW7NCK2pPVzgeuAfYHb0mJmZhUZ9GITEeuAd9aIPwOcVCMewKxu+poHzKsRbwWO6neyZmY2IIbS1GczM2tSVc1GayrjZv+46hQasn7OaVWnYGZ7KI9szMwsOxcbMzPLzsXGzMyyc7ExM7PsXGzMzCw7FxszM8vOU59tj+Dp6WbV8sjGzMyyc7ExM7PsXGzMzCw7FxszM8vOxcbMzLJzsTEzs+xcbMzMLDsXGzMzy87FxszMsnOxMTOz7FxszMwsOxcbMzPLrmmLjaQpkp6Q1CZpdtX5mJntyZqy2EgaBnwTOAU4AjhL0hHVZmVmtudqymIDTATaImJdRLwELASmVpyTmdkeSxFRdQ4DTtLpwJSI+Ou0/XHguIg4r0u7mcDMtPl24IlBTbS+g4Gnq05igDXbNTXb9UDzXVOzXQ8MvWv6w4gY2VOjZv14mmrEXlVVI2IuMDd/Or0nqTUiWqrOYyA12zU12/VA811Ts10P7L7X1Ky30dqBsaXtMcDGinIxM9vjNWuxWQ6Ml3SYpH2AM4FFFedkZrbHasrbaBGxU9J5wBJgGDAvIlZXnFZvDcnbe/3UbNfUbNcDzXdNzXY9sJteU1NOEDAzs6GlWW+jmZnZEOJiY2Zm2bnYDEHN9qodSfMkbZb0aNW5DARJYyXdKWmNpNWSLqg6p/6Q9DpJD0h6OF3PJVXnNBAkDZP0kKQfVZ3LQJC0XtIjklZKaq06n97yM5shJr1q5/8CH6CYwr0cOCsiHqs0sX6Q9F7geWBBRBxVdT79JWkUMCoiHpS0P7ACmLa7/m8kScDrI+J5SXsD9wIXRMT9FafWL5I+B7QAB0TEB6vOp78krQdaImIo/UFnwzyyGXqa7lU7EXEPsKXqPAZKRGyKiAfT+nZgDTC62qz6LgrPp82907Jb/1eopDHAacB3q87FCi42Q89oYENpu53d+F9kzU7SOOBoYFm1mfRPuuW0EtgMLI2I3fp6gK8Bnwd+X3UiAyiAn0pakV61tVtxsRl6GnrVjlVP0n7AD4DPRsRzVefTHxGxKyImULxtY6Kk3fZ2p6QPApsjYkXVuQyw4yPiGIq32c9Kt6d3Gy42Q49ftbMbSM82fgDcEBE/rDqfgRIRzwJ3AVMqTqU/jgc+nJ5xLAROlPQv1abUfxGxMf1uBm6huOW+23CxGXr8qp0hLj1QvwZYExFXVJ1Pf0kaKWl4Wt8XeD/weLVZ9V1EXBgRYyJiHMX/f+6IiP9RcVr9Iun1aTIKkl4PTAZ2q9mdLjZDTETsBDpftbMGuGk3fNXOK0i6EbgPeLukdknnVJ1TPx0PfJziv5hXpuXUqpPqh1HAnZJWUfzHztKIaIrpwk3kUOBeSQ8DDwA/joifVJxTr3jqs5mZZeeRjZmZZediY2Zm2bnYmJlZdi42ZmaWnYuNmZll52JjtpuRNK5Z3qBtew4XGzMzy87FxqwCaXSyRtJ30jdkfippX0kTJN0vaZWkWySNSO2PTd+buQ+YVepnmKR/krQ8HfPJyi7KrA4XG7PqjAe+GRFHAs8CHwEWAF+IiHcAjwAXpbbXAudHxLu79HEOsC0i3gW8C/gbSYcNSvZmveBiY1adJyNiZVpfAbwFGB4Rd6fYfOC9kg7sEr++1MdkYHr6PMAy4A0URcxsSNmr6gTM9mA7Suu7gOHdtBPdf2ZCwGciYslAJmY20DyyMRs6tgFbJf1J2v44cHd67f82Se9J8b8sHbMEODd98gBJb0tvBTYbUjyyMRtaZgDflvQHwDrg7BQ/G5gn6QWKAtPpu8A44MH06YMOYNrgpWvWGL/12czMsvNtNDMzy87FxszMsnOxMTOz7FxszMwsOxcbMzPLzsXGzMyyc7ExM7Ps/j/cB9MQuw+94AAAAABJRU5ErkJggg==\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": 24,
   "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",
      "[[1]\n",
      " [1]\n",
      " [2]\n",
      " [2]\n",
      " [2]\n",
      " [1]]\n",
      "[[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": 25,
   "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": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzsnXd8XWX5wL/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": [
       "<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": 49,
   "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": 26,
   "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.38 µs ± 444 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": [
      "2.85 µs ± 140 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.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}