13-markov.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "IN_COLAB = 'google.colab' in sys.modules\n",
    "if IN_COLAB:\n",
    "    !apt-get install -y graphviz-dev\n",
    "    !pip install pygraphviz\n",
    "    !git clone https://github.com/cs357/demos-cs357.git\n",
    "    !mv demos-cs357/figures/ .\n",
    "    !mv demos-cs357/additional_files/ ."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.linalg as la\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "import sys\n",
    "sys.path.append('./additional_files')\n",
    "sys.path.append('./figures')\n",
    "from graph import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Review of graphs and adjacency matrices\n",
    "\n",
    "How can we write the graph below as an adjacency matrix? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPEAAAFbCAYAAAAEO2g1AAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOydd1yT5/rGr0z2RplOEGSpuEHBgbgoKu5RRati1bautupprbWeHls91raOVmtbEbVurQtRREEEHBUVEAcoQ/Y0zIQkz+8PD/xqVYTkTd4Qnu/nwx9N8t7PFcrl/cz74RBCCCgUSouFy7YACoWiHNTEFEoLh5qYQmnh8NkWQGn5lJaWIi8vD+Xl5aipqYFYLEZ1dTX4fD6MjIzA4/FgamoKKysrWFlZgcfjsS1Zq6AmpjQJQgju37+P27dvIzk5GcnJyXj48CFycnJQW1vb5Dg8Hg9WVlbo3Lkz3N3d4e7ujm7duqFPnz7Q1dVV4TfQXjh0dpryJjIzM3H69GlERUXh6tWrKC4uho6ODlxcXODm5gYXFxe0a9cONjY2sLW1hbm5OXR1daGjowN9fX3U1dWhsrIScrkcpaWlKCgoQG5uLnJzc/Hw4UPcv38fycnJKC0thY6ODvr27YvBgwdj9OjR6NevHzgcDtu/ghYBNTHlJXJychAaGorjx4/jr7/+gomJCQYNGoTBgwfD19cX3bt3B5/PbAcuMzMT0dHRiImJweXLl/HkyRPY2dlh3LhxmDFjBry8vBhtT9ugJqaAEIKLFy/ip59+wpkzZ2Bqaorx48cjKCgIQ4cOhVAoVKuepKQknDhxAseOHcO9e/fQvXt3LFiwADNnzoShoaFatbQICKVVc/HiRdKnTx8CgPTq1Yvs3LmTVFdXsy2rgVu3bpGQkBBiYGBALCwsyNq1a0l5eTnbsjQKauJWSkJCAunduzfhcDhk3LhxJDExkW1JjVJUVERWr15NjIyMiKWlJfn555+JTCZjW5ZGQE3cyigrKyMLFiwgXC6XDBkyROPN+0+Ki4vJihUriEAgIH379iV//fUX25JYh5q4FREbG0vat29PrK2tyb59+9iWoxRJSUnEx8eHCAQCsnHjRiKXy9mWxBrUxK2EDRs2ED6fT9555x1SVFTEthxGkMvl5NtvvyUCgYCMHDmSlJWVsS2JFaiJtRypVEpCQkIIj8cjW7Zs0cqMlZCQQOzt7Ym7uzvJzs5mW47aoUtMWoxUKsXkyZMRERGBP/74A2PGjGFbksrIzs7GqFGjIBKJcPnyZTg4OLAtSW1QE2sphBDMmTMHR48eRUREBAYMGMC2JJVTVlaG4cOHo6ysDNeuXYOVlRXbktQCPcWkpaxZswZ//PEHjh071ioMDABmZmY4d+4ceDweRo8e3aw93S0ZamIt5NKlS9iwYQO2b9+OESNGsC1HrbRp0wbh4eF48uQJPv74Y7blqAXandYyysrK4O7ujoEDB+LQoUNsy2GNw4cPY+rUqTh9+jQCAgLYlqNSqIm1jI8//hhhYWF4+PAhTE1N2ZbDKtOnT8etW7eQkpICgUDAthyVQbvTWsTTp0+xbds2rF27lhUDi8ViXLhwARs3bkRcXBxkMpnaNfydDRs2IDs7Gzt37mRVh6qhJtYitm7dCltbW8yfP1/tbRcWFsLFxQVZWVl47733cPLkSYwdO5ZVI3fo0AEhISH47rvvIJfLWdOhclhboaYwikQiIW3atCHr169Xe9symYwMHDiQjBkzpuE1qVRKOnToQFauXKl2PX8nNTWVACCXLl1iVYcqoZlYS7h06RKKi4sRHBys9rZjYmIQGxv7Ug+Ax+MhODgY27ZtQ1VVldo11dO1a1f069cPBw4cYE2DqqEm1hKuXr0KZ2dntGvXTu1tHz9+HADg4eHx0uvu7u6oqqrCuXPn1K7p7wwbNgyxsbGsalAl1MRaQkJCAmtlbNLS0gAANjY2L73etm1bAMCjR4/UrunveHt749GjRygpKWFVh6qgJtYSMjIy4OzszErbBQUF4PF4r5Tx0dfXBwDk5eWxIasBJycnEEKQlZXFqg5VQU2sJZSUlMDCwoKVtt9U96p+Ztra2lqdcl6h/vdSXFzMqg5VQU2sJVRXV0NPT4+Vttu1aweZTAaxWPzS6xUVFQAAV1dXNmQ1YGBgAACsTrCpEmpiLcHMzAxlZWWstO3i4gLgxXHAv1Of+dg2cWlpKQDA3NycVR2qgppYS7C0tERRURErbc+dOxc6Ojq4du3aS6//9ddf6NGjB5ycnFjRVU/978XS0pJVHaqCmlhLcHNzQ2JiIittW1tb44MPPsCmTZtA/rcVv7a2FqdPn8avv/4KLpfdP7Pbt29DR0dHawsF0LuYtARvb298/fXXIISwcv3Jpk2bwOfzMWbMGAwfPhx5eXn4/PPP0bNnT7Vr+Sfx8fHo3bs3dHR02JaiEugpJi0hKSkJ3bp1w9WrVzFw4EDWdMhkMhQXF2tMVQ2pVIoOHTpg3rx5WLduHdtyVALtTmsJHh4e8PT0xG+//caqjvpbDzWF8PBw5OXlYdasWWxLURk0E2sRO3bswCeffILHjx/D1taWbTkagZ+fHzgcDiIjI9mWojKoibWI2tpadO3aFf7+/vjll1/YlsM64eHhGD16NOtDDFVDTaxl7Nu3D7Nnz0ZcXBz69u3LthzWqKmpQe/evdGlSxecPHmSbTkqhZpYyyCEIDAwEA8ePMDt27dhbGzMtiRWWLRoEf744w8kJiaiY8eObMtRKXRiS8vgcDjYvXs3KioqMG/ePO2uaPEGDhw4gJ9//hm7du3SegMD1MRaibW1NQ4dOoRTp05h6dKlbMtRKxcuXMCcOXOwYsUKTJo0iW056oGNciIU9XD06FHC4/HIp59+qpV3MP2TixcvEkNDQzJz5sxW8X3roSbWcsLCwohAICAzZ84kEomEbTkqY9++fUQoFJLp06dr9fd8HdTErYCIiAhiZGREBg4cSDIzM9mWwygSiYR88sknhMPhkI8//rhVZeB66Ji4FTB8+HDExcWhtLQUnp6eOHbsGNuSGCEtLQ0+Pj7YsWMHfvvtN2zatImVfeNsQ03cSnB3d8fNmzcxYcIETJw4EWPHjkVGRgbbshSitrYW69atg4eHB2pra3Hr1i3Mnj2bbVnswXZXgKJ+oqKiiIuLC9HX1yerV68mRUVFbEtqElKplISGhhIHBwdiZGRE/vvf/7a68e/roCZupUgkErJ48WJiYmJCDA0NyapVq0h2djbbsl5LdXU1+fXXX0mXLl0In88nc+bMIc+ePWNblsZATdxKOXToENHV1SVTpkwhmzZtItbW1oTH45ExY8aQM2fOaESGS0pKIkuXLiVmZmZEKBSS2bNnk7S0NLZlaRzUxK2Q77//nnC5XPLRRx8RmUxGCCFELBaTQ4cOkSFDhhAOh0PMzc1JcHAw+fPPP0lFRYVadEmlUnLz5k3yr3/9i3Tt2pUAIJ07dyYbNmwg+fn5atHQEqF7p1sRMpkMS5cuxfbt2/HFF1/gyy+/fO3nnj59iuPHj+P48eNISEgAl8tFz5494ePjg379+sHd3R1dunQBn9+8wjB//fUXevXq1fDfZ86cQVFREfLz8xEbG4vY2FiIRCJ07NgR48ePx/jx4+Hl5cV6eR9Nh5q4lVBbW4vg4GD8+eef2LNnD6ZOndqk5woKChAdHY2rV6/iypUrSE1NhUwmg0AggKOjIzp16gRra2vY29vD2NgYRkZG4PP5MDIyglgsRnV1NcRiMcrLy/H9999jwIABKCoqwsOHDxuqc9rY2MDX17fhx93dXZW/Cq2DmrgVUFJSgrFjxyI1NRUnT56Ej4+PwrFqa2uRmpqK1atXIzk5GcOGDUN+fj5ycnIgEolQUVGBuro6VFZWQigUwsDAALq6uuDz+cjOzoa3tzf69u2LLl26QCqVYsmSJRg6dCgiIiKandkp/4PNvjxF9aSnpxMnJyfSqVMnkpqayljcd955h8yYMaPJn1+yZAkBQIYNG9bwWlxcHAFAOBwO+fDDDxnT1tqggw0t5vr16/Dy8oKJiQni4+PRtWtXxmKnpaU1uQQsIQSHDx8GAERFReHZs2cAAIFA0PD+1q1b8fPPPzOmrzVBTaylnDx5EkOHDoWnpycuXbrEaPE6uVyOp0+fNtnEt27darhUjcfjYe/evQD+38T1fPDBB4iKimJMZ2uBmlgL+fHHHzFhwgRMmzYNZ86cgZGREaPxs7OzIRaL4ejo2KTPnzhxouHGxLq6OuzatQuEkFdMTAhBUFBQw1WplKZBTaxFEELw5ZdfYunSpVizZg12796tksmi9PR0AGiyiQ8dOgSJRNLw35mZmbh27dorJpbL5aiursbIkSNRXl7OnGAth5pYSxCLxZg+fTo2bNiAsLCwN64BM0FaWhoMDQ0bLhFvjOTkZDx58uSl1wQCAX777bdXTAy8KPaelZWFSZMmNVyNSmkcamItoKysDMOHD8f58+cRERGBGTNmqLS99PR0dOnSpUmfPX78+Ctmraurw8GDB1/Kzv98//Lly/jkk0+U1toaoCZu4Tx9+hTe3t5IT0/HlStXMHjwYJW3mZaW1uSu9OHDh1FXV/fK62KxGBEREW98TiaTYcuWLdi1a5fCOlsL1MQtmJs3b8LLywsCgQAJCQno3r27WtpNT09vkokzMjKQkpLyxvf37dv31hiLFi3C5cuXm6WvtUFN3EK5cOEC/Pz84OHhgdjYWNjb26ut7SdPnjRpeeno0aNvnFiTy+W4efPmW2PI5XJMnDgRT58+bbbO1gI1cQvkt99+Q0BAACZMmIBz586ptUB8fn4+KioqmpSJDx8+3OjkVGMz5/XjaFtbW3z44Yd0S2ZjsLpfjNIs5HI5Wbt2LeFwOGTt2rWsaLh69SoB8NYCAjk5OYTD4RAATf7hcrmEy+USHo9HJk6cSC5evNgqC981F2riFoJYLCYzZswgfD6f7Nq1izUdYWFhREdH563m2r59e5PNy+FwCIfDIb6+vuTTTz8lAMjDhw/V9I1aPrSP0gIoLy9HUFAQbt26hVOnTmHUqFGsacnKyoK9vf1bq0q6uLhg165dMDU1BQDo6uriX//6F5ycnLB69WpwOBzo6Ohg0qRJmDBhAkJCQtC+fXvIZDKEhobi0KFDWLNmjTq+UouHHkXUcHJychAQEIDCwkKcPXsWnp6erOpZvHgxUlJScOXKlWY/6+XlBS8vL3z33XeNfu7DDz9EZGQkUlNTFVTZuqATWxpMUlIS+vfvD5lMhoSEBNYNDADPnj1Taia8KXWhp02bhgcPHuDu3bsKt9OaoCbWUCIjIzFw4EA4OTkhNjYW7du3Z1sSgBeHH9q1a6fQs03t9Hl5eaFjx444ePCgQu20NqiJNZA9e/Zg9OjRGDduHMLDw2FiYsK2pAaePXsGOzs7hZ9vSibmcDiYMmUKDhw40GTjt2aoiTWMb7/9Fu+99x4WLlyIPXv2NBzh0wTEYjGKi4sV7k43x5BTp05FVlYWEhISFGqrNUFNrCFIpVIsWLAAn332GXbs2IEffvhB4+4VevbsGQghSpm4qd+pR48ecHV1pV3qJkBNrAFUVlZizJgx2L9/P/7880+8//77bEt6LfVldVQ9sVXP5MmTcfToUcjlcoXbaw1QE7NMXl4efH19kZiYiOjoaAQEBLAt6Y1kZ2dDKBQ26Rzx62ju+HbChAnIzc2lXeq3QE3MIikpKejfvz/EYjESEhJeKqyuidRPailTzL05mdjd3R3Ozs44fvy4wu21BqiJWSIqKgoDBw6EnZ0doqOj0aFDB7YlvZWcnBylutKKzDQHBQXh2LFjdJa6EaiJWeDIkSMICAiAn58fLl26BEtLS7YlNYns7Gyljzw2d7Ju/PjxyMjIwJ07d5RqV5uhJlYzP/zwA6ZOnYqQkBAcPnwYenp6bEtqMsru1lIkm/bu3RsdOnSgXepGoCZWEzKZDIsWLcKKFSvw448/4ocffmhxF4Upu9GjOUtM9XA4HAQFBeHIkSMKt6vttKy/ohZKVVUVxo0bh9DQUBw/fhyLFy9mW1KzqaurQ1FRkVImBprfnQZedKkfPnxID0S8AWpiFZOfn49BgwYhPj4eFy5cwJgxY9iWpBCFhYWQy+WwsbFROIaik1MDBgyAjY0N7VK/AWpiFZKWlgYfHx+Ul5cjLi4OAwYMYFuSwhQUFACAUtfBSCQShbaRcrlcBAYG4tSpUwq3rc1QE6uI+Ph4eHl5wcLCAvHx8XBycmJbklKwaWIACAgIwK1btxp0UP4famIVcOzYMfj5+cHHxweXL19GmzZt2JakNAUFBdDT01PqXieJRPLaWx+agp+fH4RCIc6fP69w+9oKNTHD/PDDD5g8eTLmz5+Po0ePtqglpMYoKChQ+mbFuro6hTOxgYEBBg0ahLNnzyqlQRuhJmYImUyGDz/8EMuWLcOaNWta5BJSYzBhYmW60wAwevRoREREvPZGidaM9vyVsUhtbS2mT5+OX375BQcOHFDpZWZsoQkmfueddyASiXDt2jWldGgb1MRKUlJSgmHDhiEyMhIXL17E1KlT2ZakEjTBxJ07d4azszPtUv8DamIlePLkCby9vZGbm4tr167Bx8eHbUkqg+0xcT0BAQHUxP+AmlhBrl+/Di8vL5iYmCA+Ph5du3ZlW5JKKSwsVMrEUqkUcrlc4dnpegICApCamtpw0TmFmviNJCYmvvG9kydPYujQofD09MSlS5eUzlCajkwmQ0lJidJrxACUzsQ+Pj4wMjLChQsXlIqjTVATv4ajR482VNv4Jz/++CMmTJiAadOm4cyZM0qtm7YUiouLIZPJNMLEAoGgYf2d8gJq4n8gFouxfPlyVFVVYcSIEcjMzATwYt/vl19+iaVLl2LNmjXYvXt3q7mpj6ndWoDyJgaAIUOGICoqitbe+h/UxP/ghx9+QG5uLgghKC8vh7+/P/Lz8zF9+nRs2LAB+/bt08olpMbQNBMPHToUJSUlSEpKUjqWNtA6UkkTKSoqwldffdVwp25dXR0yMjIwePBgFBYW4sKFCxg0aBDLKtVPQUEBhEKhUkXsmTRxjx49YGlpiaioKHTv3l3peC0dmon/xhdffNHwx1ZPXV0d0tPT4ePjo9VLSI1RUFCAtm3bKlUHu36XlbKz08CLU02+vr50XPw/qIn/R2pqKnbt2vXaLX1SqRRnzpzB559/zoIy9ikuLlb6EAeTmRh4MS6+cuUKpFIpI/FaMtTE/2PJkiXg8XhvfF8ul2PDhg3Ytm2bGlVpBiUlJUoX8xOLxQAAHR0dJiRh6NChqKiowF9//cVIvJaMRo2JZTIZCgoKUFBQgPLycshkMlRUVEAqlUJfXx86OjrQ09ODqakpbGxsYG5uzki7ERERuHjx4ls/x+FwsGTJEjg4OLB60be6KS0tVfp3XVVVBeDFaSQmcHV1hZ2dHaKiotCvXz9GYrZUWDFxTU0Nbt68iXv37iE5ORkpKSlIT09HYWFhw6RSU9DV1YW9vT2cnZ3h7u4ONzc39OzZE66urk0ev8lksoYs/Ka2BQIB6urq0KtXL3zwwQcYPHhwkzVqAyUlJXB1dVUqRnV1NQBAX1+fCUkAgEGDBiE6OhqrV69mLGZLRC0mlsvluH79Os6dO4crV67g5s2bEIvFMDc3h7u7Ozw8PDBx4kTY2NjA1tYWVlZWMDc3B5fLhZGREfh8PqqrqyEWi1FbW4vS0lLk5uYiLy8P2dnZuH//Ps6fP48tW7ZAIpGgTZs28PHxwZAhQzBmzJhG7/bdtWsXHj9+/MqaI4fDAZfLhVAoxLvvvotFixahR48eqv5VaSRMZGJVmHjAgAFYvXo1ZDJZo0MhrYeokNjYWLJw4UJia2tLABAHBwcyd+5csnfvXpKZmcl4e3V1deTmzZtk8+bNZMyYMcTExIRwOBzSu3dvsmHDBpKTk/PS50UiEbGwsCAcDocAIACIQCAgAIiHhwfZuXMnqaioYFxnS8POzo589913SsUIDQ0lurq6DCl6wZ07dwgAcvfuXUbjtjQYN3FFRQXZvn078fDwIABIt27dyLp168i9e/eYbuqtiMViEh4eTkJCQoilpSXh8/kkKCiIXLhwgRBCyMqVKwkAwuFwCI/HI7q6umTu3Lnk1q1bateqyejp6ZE9e/YoFeOnn34i5ubmDCl6gUwmIyYmJmTHjh2Mxm1pMGbiiooK8v333xNra2uiq6tLJk2aRC5evMhUeKURi8Xk8OHDZNiwYYTD4ZCuXbsSPp9PABBXV1eyY8cO8vz5c7Zlahw1NTUEADl9+rRScTZv3kzs7e0ZUvX/DB8+nMyYMYPxuC0JpcfEcrkcO3fuxJo1a1BXV4ePPvoIy5YtY2zmmCmEQiEmTZqESZMm4fbt2xg3bhykUilcXV2xZ88e9OnT543PyuVy5Ofno6SkBCUlJSguLm6ooaynpwdLS0tYWFjAzs6O0TGfJlBSUgIAjMxOMzUz/XcGDBiA33//nfG4LQmlTHz79m28//77uHPnDpYsWYLVq1drnHlfR48ePZCUlIQnT55g+fLl6N+/P0JCQrBhwwYQQpCQkIDr16/j3r17ePToEdLS0hrWORuDw+GgQ4cOcHV1RZ8+feDn54f+/fszskuJLZgycU1NjUr+gRswYADWrl2r9D1RLRmFNnsQQrBx40b0798fenp6SExMxKZNm1qEgYEX2/ZMTEzg6emJqKgohIaG4uTJk3B1dYWFhQVGjx6N/fv3QyAQYNy4cdi5cydiYmKQkpKCgoICSKVSkBdDEVRVVSErKwu3b9/GiRMnEBISAjMzM+zZswe+vr4wNzfH7NmzERMT0yKv5ywtLQUAWFhYKBWnurpaJZm4f//+4PP5iI+PZzx2i6G5/e/S0lIyfPhwIhAIyMaNG4lcLme2g88SRUVFZPTo0YTH45G1a9cyEvPx48fkxx9/JL169WqY8T5+/HiL+p0dO3aMcDgcUldXp1ScefPmkeHDhzOk6mV69uxJlixZopLYLYFmmTgrK4u4ubmRdu3akRs3bqhKE2vI5XKyefNmwuPxyMKFC4lUKmUs9p07d8jkyZMJh8Mh/fr1I/fv32cstirZtWsXMTExUTrO9OnTybhx4xhQ9Coffvgh6d27t0pitwSa3J1OS0uDt7c3uFwu4uPjG50IaqlwOBwsX74cR48exZ49ezB16lTGNth3794dhw4dwu3btwEAvXr1wtatWxmJrUpKS0uV7koDqutOAy/GxXfu3EFlZaVK4ms6TTJxXl4eRowYAVtbW8TExCh9vaWmM27cOERERODcuXNYsGABo2PZHj16IDY2FqtXr8bSpUuxbNkyjR4rM7FbC3gxO62qmfuBAwdCKpXi5s2bKomv6bzVxNXV1Rg1ahSEQiHOnj0LU1NTdehiHR8fHxw5cgRhYWFYt24do7H5fD7WrFmDP/74Az/99BMWLlzIaHwmYcrEqszEdnZ2aN++fastKv/WJaYVK1YgKysLiYmJSh9Ha2mMHj0aW7duxaJFizB48GDGDz5MnjwZenp6GDduHNzd3fHBBx8wGp8JSkpKGOtOq3INfcCAAa3WxI1m4tOnT2Pnzp345Zdf0KFDB3Vp0igWLFiA8ePHY+bMmSgvL2c8fmBgINavX49ly5Y1WiaXLcrKymBmZqZ0HHWYOC4urlmn4LSFN5pYIpFg2bJlmD59OiZMmKBOTRrHrl27IBaL8Z///Ecl8VevXg0vLy988MEHGjc+Li8vZ2QIpQ4Ti0Qi3L9/X2VtaCpvNPGOHTuQm5uLDRs2qFPPK+Tn5+PKlSusajAzM8OaNWuwdevWhhK2TMLhcPDdd98hISEBp06dYjy+MohEIhgbGzMSR5U1uj08PGBsbNwqu9SvNbFcLseWLVuwcOFCtGvXTt2aALyoPPnxxx+jc+fOOHHiBCsa/s6CBQvQtm1bbN++XSXxe/fujREjRmDXrl0qia8oTJiYEILKykpG/jF4EzweD3379qUmricyMhJZWVmYP3++uvU0kJGRgVmzZqGmpoY1DX9HKBRizpw52Lt3r8rux507dy4iIiKQk5OjkviKIBKJlCpVCwCVlZWQyWQqNTHQeie3XmvigwcPwsvLi9VLwvr06aNxl5TNmTMHhYWFKiuVGhgYCKFQiMjISJXEby41NTWQSCRKm08kEgGA0v8YvI0BAwbg6dOnGvWPoDp4rYmvXr0Kf39/dWvReDp06ABHR0fExsaqJL5QKETPnj1x/fp1lcRvLkyZ7/nz5wCg8kzs5eXVKg9DvGLi4uJipKenw8vLiw09Go+3t7dK/0j69OmjMWVYmTJf/T8GqjaxoaEh3N3dW12X+hUTZ2ZmghACZ2dnNvRoPE5OTsjIyFBZfBsbGxQWFqosfnNgynzqMjHwIhtrSk9GXbw2EwPKnx/VViwsLBoOyqsCc3NzlJWVqSx+c6jPxMp2p0UiETgcjlquge3Tpw8SExNfuY5Hm3nFxPWzwXp6emoX0xIwNDRsKISuCvT19VUavznUZ1BlzScSiWBoaAguV/UXjvTt2xe1tbWt6sbEV36r9VvsNCUbaBolJSUqrWBSWVkJQ0NDlcVvDiKRCPr6+kqXF3r+/LnKZ6brcXFxgbGxMW7cuKGW9jSBV0xc340uKipSu5iWQFFRkUqHGiUlJYzsVWYCpszH1K6vpsDlctGrV69WdSzxFRN36dIFurq6GrEZv743UFtby7KS/+f27dvw8PBQWfxHjx6hS5cuKovfHJgyX0VFhdpMDLzoUrfqTKyjowNPT0/ExcWxoaeB8PBwLFmyBABw8uRJ7N69G/n5+axqIv+rhKnK5bfbt2/D3d1dZfGbAxO7tYAXGV2dJu7Tpw9SU1MbxvTazmtnGvz8/HD69GlWj3WNGjUKBw8eBCEEBQUFmDdvHqytrVnTAwAxMTEoLS3F0KFDVRK/oKAASUlJKovfXJgyH1P/GDSVvn37Qi6XN5RC0nZea+I5c+YgJyenSdd9tiZ+/fVX9OnTR2WZ8uDBgzA0NNSYWxeZMp86x8QA0K5dO9ja2raaLvVrTdy5c2f4+vpiy5Yt6tajsTx79gxHjx7FvHnzVBJfJpNh586dmDp1qvZaRgsAACAASURBVMrK2DQXJo8hqtPEwIsudWuZ3Hrjwt1XX32FCxcu4MKFC+rUo7GsWbMGVlZWCA4OVkn8ffv24fHjx1i+fLlK4isCU2eA1T0mBl6YuFVnYgDw9fVFYGAgli9frlGzw2yQkJCAsLAwfP3119DR0WE8flFREVauXIl58+Zp1MmtqqoqRtas1T07DbwYF2dlZSEvL0+t7bJBo1totm/fjtzcXHz66afq0qNxVFZWYtasWRg2bBimTZvGeHyZTIbg4GDo6uri22+/ZTy+MlRXVzOyc0/VG2ReR9++fcHlcltFl7pRE7dr1w47duzAtm3bcOjQIXVp0hjkcjmCg4NRUVGB0NBQcDgcxuMvXLgQV65cwaFDh9Serd4GEzcZisViVFVVqd3EJiYmcHR0bBUmfmvJ2qlTp+L69esIDg6GpaUl/Pz81KFLI/jggw9w7tw5XLhwAVZWVm/8XFFRETgcTrNK+orFYsybNw+HDx/GkSNH0K9fPyYkMwoTxe2YulVREVrNpo+m3PUik8nIjBkziJGREYmKilLZnTKaglwuJytWrCA8Ho+cOHHilfezs7PJvn37yIIFC4iDgwPhcDjkwoULTY7/8OFD0qtXL2JsbEwiIiKYlM4ourq6ZO/evUrFSEpKIgBYuXvqhx9+IKampi3qAjtFaPKFahKJhEydOpXo6OiQP/74Q5WaWEUsFpPp06cToVBI9u/fTwgh5NGjR+TXX38lwcHBpF27dgQA4XK5REdHhwAgAEhycvJbYz9//pysWbOG6OjokF69epHHjx+r+usojFQqJQDI0aNHlYoTHR1NAJD8/HyGlDWd+Ph4AoA8fPhQ7W2rk2bdiiiTycjy5csJh8MhK1euJBKJRFW6WCEjI4N4e3sTQ0NDsnTpUvLuu+8SGxsbAoDweDzC5/MbTPvPn9LS0jfGffLkCVm9ejUxMzMjpqamZMuWLUpfFapqRCIRAUDCw8OVinP8+HECgIjFYoaUNZ3a2lqio6NDwsLC1N62Omn2/cSEEPLrr78SAwMD0r9/f5KWlsa0JrXz4MEDsmDBAiIQCAiXyyUACJ/PJzwe742m/fuPUCh8qcsmFovJlStXyGeffUb69u1LuFwusbGxIV999RUpKytj8Zs2nfz8fAKAREdHKxVn9+7dxMjIiCFVzadXr15k6dKlrLWvDt46sfU63nvvPXh5eWHatGnw8PDA6tWr8emnn6pkDVXVPH36FBMmTEBKSgo6duyI58+fo7y8HHK5HHK5vEkxjIyMsGzZMjx9+hRPnjxBeno6ampq4OjoCH9/f3z++ecYNWoU+HyFft2sUF1dDQCMTGyxWSWme/fuuHv3LmvtqwOFSy24uLjg5s2bWLduHTZu3Ah3d3fs27evxdyFU1RUhFWrVsHNzQ1yuRw2NjbIyMhATU0NBg8eDFtbW3A4HPB4vLfG0tHRwY0bN2BkZIRx48bhp59+wtOnT/H48WPs2LEDgYGBLcrAABqqiyhr4rKyMlZmputpDSZWqDv9T7Kzs0lwcDDh8/nE2dmZ/P7776SmpoaJ0IyTlZVFPvnkE2JoaEjatm1LtmzZQiQSCTl58mRD97h+7Nu5c2fSvn37l1775w+HwyETJ05k+2sxzvXr1wkAkpGRoVSc+fPnE39/f4ZUNZ/6ibXs7GzWNKgaRkxcz+PHj8ns2bOJUCgk5ubmZOnSpU2atVU1EomEnD59mgQGBhIej0esra3Jf//7X1JVVfXS54YMGUIEAkGDQblcLuHxeERXV5e0a9eOcDicl97H/8bDS5YsYembqY7Lly8TAKSwsFCpOBMmTCBTpkxhSFXzKSsrIxwOh5w+fZo1DaqG0cpljo6O+P3335GVlYWPP/4Yf/75J9zd3eHi4oLPP/8ct27dUlt3u6KiAqdOncLs2bNhbW2NMWPGoLq6Gn/88QcyMzOxYsWKV7qKW7ZseUmfXC6HTCZDbW0tCgoKQAiBvr4+OBzOS3Wn2D7nrAqY6k6zseXy75iamqJDhw5a3aVWyUDNysoKq1evxsqVK3Ht2jUcP34c+/btw9dffw0TExMMHDgQAwcORM+ePeHu7g5bW1ul2pNKpXj8+DGSk5ORkJCAq1evIjExEXK5HJ06dcKYMWOwdu1adOzYsdE43bt3R3BwMPbt2/fKfUv1JVArKioAvLjJkMvlQiKRwMbGRin9mkh1dTU4HI7Se6dLS0tZL3+s7eNilc62cLlc+Pj4wMfHB1u2bEFycjKio6MRExODrVu3Ijc3F8CLLXlOTk6wtrZGu3bt0LZtW5iYmEBHRwf6+vrQ0dFBRUUFpFIpKioqIBKJkJ2djYKCAmRlZeHRo0eQSCTg8/lwcXHBoEGD8Mknn8DX1xfr16/HL7/8gqlTp77VxADw9ddf4+DBg2+8NK1+xloikYDD4WDMmDFwdHRk7HemKVRVVUFXV1fpMrOlpaWsF/7r3r07Dh48yKoGVaLWKVN3d3e4u7tj8eLFAF50tZKSkpCSkoK0tDTk5+fjzp07KCgogEgkatg8L5FIYGhoCIFAACMjIxgbG8POzg7W1tbw9PRE165d4ebmBldX11eWubZu3QqxWIygoCCEh4dj0KBBjWq0sbHB6tWrsW7duka7/hwOB99++y0++eQT5X8xGkh1dTUjxQnYXmICXpj43//+t8ovOmcNtgfl6kAqlZIpU6YQY2NjcuPGjbd+vrq6mtjY2DRs/PjnD4/HI/Pnz1eDcvbYtGkTad++vVIxampqCABy6tQphlQpRmpqKgFAbt++zaoOVaH6kvwaAI/HQ1hYGHx8fDB8+PC3luPV09PDxo0bQQh55T2BQABvb2/s2LFDVXI1gtraWujq6ioVo/4EE9uZ2NHRETo6Orh//z6rOlRFqzAx8MJ8R44cQY8ePTBy5Eg8ePCg0c/PmDEDnp6eL23SEAgE6Ny5M+Li4rB9+3ZVS2YViUQCoVCoVAw2jyH+HT6fD0dHR6SmprKqQ1W0GhMDLzLsmTNn0KVLFwwfPrzR2w05HA62bdvWMC7m8/mwsrJCdHQ0Nm/ejGXLlmHVqlVqUq5+mDBxfZ1wTViCc3FxoSbWFgwMDHDmzBlYWlrC39+/0RpMXl5eCAoKanju0qVLsLKywpIlSxAaGorNmzdj4cKFTd5j3ZJgwsSFhYUQCoVqrTn9JlxcXGh3WpswNTVFREQEBAIBRowY0ehVpZs2bYKJiQnOnTsHJyenhtdnzpyJ48ePIzQ0FBMmTNC6YoJMmLigoABt27ZlvKyRIri4uCAtLQ1isZhtKYzTKk0MAG3atMGFCxdQWVmJYcOGoby8/LWf69y5M5KSkuDt7f3Ke4GBgQgPD8fly5cREBDQsBFEG2AqEzdW1kiduLq6QiqVIj09nW0pjNNqTQwA9vb2uHjxIgoLCxEQEPDGe4HbtWv3xhiDBg1CVFQUkpOT4efn13BJe0uHqUysKSZ2dnYGl8vVynFxqzYxADg4OCAiIgIPHz7ExIkTFeoW9+zZE/Hx8SgtLYWvry+ys7NVoFS9SCQSpc+Ha5KJdXV1YWdnRzOxtuLu7o6IiAhcv34d48ePV2jc1LlzZ1y9ehUCgQA+Pj54+PChCpSqDybHxJqCg4MDnjx5wrYMxqEm/h+9evXCpUuXcP36dYwdO1ahjGxjY4MrV67Azs4O3t7euH79ugqUqgdtGxMDL/6hpSbWcjw9PREZGYmbN29i3LhxChnZzMwMFy5cQN++fTFs2LAWe7OksiYmhKCoqIhmYjVATfwPmDCygYEB/vzzT7zzzjsIDAzE0aNHVaBUtShr4rKyMkgkEo3LxJmZmW88odZSoSZ+DZ6enrh48SJu3ryJoKAghYwsFAqxf/9+vPfee5g6dSr27t2rAqWqQ1kTFxQUAIBGmdjBwQFSqVQrJh7/DjXxG+jZsycuXryIGzduKGxkLpeLHTt2YOXKlZgzZw527typAqWqgSkTa1J3unPnzgCgdV3qllWCUc3UG9nf3x9BQUE4ceKEQid7vv76axgbG2PhwoWoqqrSqDuI34SyJi4sLASPx2vW/VSqxsLCAsbGxo3umW+J0Ez8FuqNfP36dYUzMgCsXLkSW7duxccff9wiDk5IJJKX6og1l4KCAlhYWDSp5K86sbOzw7Nnz9iWwSg0EzcBpjLy4sWLIRAIsHDhQgDAN998w7RUxqirq1M6E2vSeLgee3t75OTksC2DUWgmbiK9evVqyMjjx49XOCOHhIQgLCwMmzdvxuLFi19beEATkMlkStXXKiwsRJs2bRhUxAx2dnbUxK2ZeiMnJCQovLMLAKZPn479+/fjl19+wYIFCzTyKCMhRKnTR7m5uUpXMVUF2tidpiZuJr169cLZs2dx7do1TJo0SWEjT548GSdOnEBYWBhmzpwJqVTKsFJ2yc7ObvTgCFvQTEwB8KJYwPnz5xETE4OxY8c2XD7WXAICAnDixAmcOHECM2bM0KhNCMpm4mfPnsHe3p5BRcxgZ2eH0tJShf+faSLUxAri5eWFqKgo/PXXXxg5ciREIpFCcUaOHInz588jPDxcqdlvplHGxNXV1SgpKdFYEwNoqHmuDVATK0HPnj0RExOD9PR0+Pn5NVohpDF8fX0RFRWF+Ph4jBo1CpWVlQwrVS/1Y05NNHH95pOioiKWlTAHNbGSuLi4ICoqCvn5+fD391f4j6N3796IjIxEcnIyRo8erXBmZwplMnG9iTVxTFxfPre0tJRlJcxBTcwAzs7OiI2NhUgkwqBBgxTuqnl6eiIqKgqPHj3CiBEj8Pz5c4aVqodnz55BR0dHo3Zr1aOvrw89PT2Fe02aCDUxQ3To0AFXr14Fl8vFkCFDFN5k7+HhgejoaGRlZSk11lYWZTOxvb29RhTIex3m5ubUxJTXY2Njg0uXLkFXVxc+Pj5IS0tTKI6zszOuXLmCrKwsDB06FGVlZQwrVS05OTkaOR6ux8LCgpqY8masrKxw5coVWFlZwcfHBykpKQrF6dKlCy5fvoy8vDwEBASoPSMrk4mzs7OpidUINbEKqK/u0blzZwwdOhT37t1TKI6TkxMuX76MzMxMtXetle1Oa+KkVj3UxJQmYWJiggsXLsDDwwODBw9WuN4Wm0ZWlGfPnjWsx2oi5ubmLW6I0hjUxCrEwMAAZ8+ehY+PD0aMGIG4uDiF4vzdyKNGjVKLkRXNxLW1tSguLtbo7rSBgcEba4y3RKiJVYyOjg4OHz6MwYMHY8SIEbh8+bJCcZycnBAVFYWMjAyMGjVKLbdNyGQylJWVoaysDNnZ2Xjy5MlbN6Lk5OSAEKLR3WlqYkqz0dHRwZEjR/DOO+9g9OjROHnypEJxnJ2dG4w8cuRIxoy8aNEitGnTBkZGRjAyMgKXy0VeXh5WrFgBc3NzmJubo3379nBwcHjrCSBN3q1Vj4GBQYvfFfd3qInVhEAgwIEDB/D+++9j4sSJ+OWXXxSKU2/kp0+fMmbk+utnKisrUVlZ+cYzzo6OjujatWujsbKysqCjo6ORZ4nroZmYojAcDgdbtmzB119/jQULFuDLL79UKI6zszMuX76Mp0+fMtK1DgwMhKmpaaOfEQqFmDZt2ltjpaeno1OnTkoVFFA1hoaG1MQU5aivt7V+/XosWbJEoeoe9UZ+8uSJ0kYWCoWYNWtWo+V4JBIJxo8f/9ZY6enpcHR0VFiLOqjPxJpaVaW5UBOzxOLFixEWFoaffvoJwcHBChUFcHZ2xsWLF5GWlobAwEClzsjOnj0bEonkje/b2dmhR48eb42TlpYGBwcHhXWoAwMDAxBCUFNTw7YURqAmZpHp06cjPDwcJ06cwPjx4xX6o3Jzc8OlS5eQkpKC8ePHN2rExvD09ISbm9trl5Wa2pUGWoaJ9fT0AICamMIMfn5+uHTpEuLi4jBq1CiFTi65ubkhMjISN27cwJQpUxQu9TN//vzXjmWb2pUWiUQoLi7W+O50fRldmUzGshJmoCbWAPr27dtQXGDgwIFvPMookUjeeKiie/fuOHv2LCIjIzFnzhyFiu+9++67rzVx27Zt0a9fv7c+X69N0zNxvYm1pa4ZNbGG4OrqiqtXr0IikWDgwIGvmJUQgrlz52LmzJlvnJDx8vLCyZMncfToUcyfP7/ZEzcWFhYICAh4qWi8UCjE5MmTmzTbnJaWBh6Ph44dOzarXXXD578ot05NTGGcjh07IiYmBqampvDx8cHdu3cb3lu7di0OHDiAhISERm9Z9PPzw6FDh7B3716FrouZO3fuSwX7JBIJgoKCmvRseno62rVrp/S9xqpG27rTIBSNo6KiggwbNoyYmZmR2NhYsnv3bgKAACAcDofY29uTmpqaRmMcOXKE8Hg8sn79+ma1XVdXRywtLRvaMzExIXV1dU16du7cuWTYsGHNao8Nbt++TQCQR48esS2FEWgm1kAMDQ1x6tQp+Pj4YMGCBQgJCWl4jxCCvLw8/Pjjj43GqN8V9sUXX+C///1vk9vm8/l47733IBAIIBAIEBQU1ND9fBstYWYa+P/uNM3EFJWTmJhI9PX1CZfLbciM9T96enokLy/vrTF++OEHwuFwyM6dO5vc7oMHDxraOXv2bJOfs7e3Jxs3bmzy59kiJSWFACBJSUlsS2EEeqGahpKTk4ORI0dCIpG8dqZZKpXiyy+/xM8//9xonI8++gilpaVYuHAhjIyMmrTe6+zsDE9PTzx48ADGxsaIjIyEXC5vWP4yMTEBl8uFsbExLCwsYGtrCw6Hg9zc3BaRiesn6TTx+hxF4BCiJXvPtAiRSIT+/fsjLS2t0VshuFwu7ty5Aw8Pj7fG/PTTT/H999/j2LFjCAwMbHi9oKAACQkJSElJQVJSElJTU5GVldXsQ/OmpqaQyWTw9fXFgAED4Obmhn79+mnkzYgpKSlwd3dHSkoKXF1d2ZajNDQTaxgSiQTjxo1DamrqWz/L4/GwbNkyREZGvvWz3377LcrLyzF58mRs2rQJycnJiI6OxoMHD8DhcNCxY0e4ublh1KhR6NChA8zMzJCfn4+goCAYGRmBw+E0HJIoLy8HIQQikQhFRUXIy8tDTk4OMjMzkZycjF9++QUZGRkghMDFxQW+vr7w9/fHyJEjYWBgoPTvSFnqx8KadneywrDbm6e8jmvXrpGZM2cSoVBI+Hw+4XA4r4yJ//5z7ty5RuPV1NSQsLAwEhgYSIRCIeHxeMTb25usWrWKhIeHE5FIxPh3eP78OTl37hxZtWoV8fLyIjwej+jp6ZFx48aR/fv3k9raWsbbbCraNjtNTazBlJeXk507d5KuXbsSAEQoFL5iYB6PRxwcHIhEInnl+cePH5MVK1YQCwsLIhAIyOjRo8nu3btJYWGh2r9LYWEh2b17Nxk9ejQRCATE0tKSfPzxxyQtLU3tWm7dukUAkPT0dLW3rQqoiVsIt27dIvPnz2/Izn83MpfLJdu3b2/4bHp6OgkJCSF8Pp/Y2dmRlStXkqysLBbVv0xeXh755ptvSMeOHQmXyyWTJk0iDx48UFv7169fJwBIRkaG2tpUJdTELYyioiKyadMm0rFjRwKACASChk0ZaWlpZO7cuYTH4xEXFxeyf/9+IpPJ2Jb8RqRSKQkLCyPOzs6Ex+OR+fPnk+LiYpW3GxcXRwCQ7OxslbelDqiJWyhyuZxERkaSiRMnNmRmPT09YmdnR/bt26fR5v0nUqmU7N27l9ja2hILCwuye/duIpfLVdbe1atXCQCSm5ursjbUCTVxC6eoqIj4+/sTLpdLzM3Nyd27d9mWpDDPnz8nS5cuJXw+nwQGBqosK1+5coUAIPn5+SqJr27otssWTFxcHHr06IFHjx4hJiYGRUVFsLGxYVuWwhgbG2PLli24fPky7t69ix49eiAhIYHxduovcq8vDtDSoSZuoZw8eRLDhg1Dz549kZiYiAEDBoDL5Wp0lcmmMnDgQCQmJqJ79+7w8/PDqVOnGI1fX9FDV1eX0bhsQU3cAtmzZw8mTpyI2bNn48SJEzAzM2NbEuOYm5vjzz//xLvvvovx48dj7969jMWuqakBj8fT+COTTYXu2GphnDx5EvPmzcOqVavw73//m205KoXH42Hnzp0wNzfH3LlzYWZm9tKWUUWpqanRmiwMUBO3KK5fv45p06YhJCRE6w38dzZs2IDS0lJMmTIF0dHR6NOnj1LxampqtGY8DNDudItBJBJh+vTpGDJkCLZt28a2HLWzY8cO+Pj4YPr06UoXy6+traUmpqifDz/8ENXV1dizZ49G366gKng8HkJDQyESifDRRx8pFYtmYoraSUhIQFhYGH7++We0bduWbTmsYW1tjZ9//hmhoaG4ceOGwnG0zcT0PHELYODAgRAIBApfi6os5eXl+PXXX5GVlYWAgAD4+fmxeozP19cXHA4H0dHRCj2/fPlyxMfHIz4+nmFl7EAzsYYTHx+Pa9euYcOGDay0X1pait69e+Pu3btITk7GqFGj4O3tzYqWev7zn/8gJiYGN2/eVOj5iooKGBkZMayKPaiJNZzffvsNHh4e6N+/PyvtHz58GDdu3MDevXtx6dIlfPnll7hx4wauXbvGih7gRc+kW7du+O233xR6npqYojZkMhmOHj2KOXPmsNK+RCLBiBEjYG5u3vDarFmzALzYIskms2bNwqFDhxSqk1VRUcG6fiahJtZgkpKSUF5eDn9/f1baFwqF6NSp00uv3bt3D++8806T6nqpEj8/P5SVlSElJaXZz4pEIpqJKeohPj4eJiYmGlHMjRCCw4cPY9WqVfjpp5/YlgMPDw8YGRkhLi6u2c/STExRG5mZmXB0dGR9XbiqqgoLFizAnDlzcP/+fXh4eCg8qcQUPB4PDg4OyMzMbPazdExMURslJSWwsLBgWwYMDAywa9cuVFRUYMuWLaioqMDChQvZlgVLS0uUlJQ0+znanaaojerqao3alMDlcrF06VKMHz8eiYmJEIvFrOoxMDBAZWVls5+j3WmK2jAzM2t2EXd14O/vD3Nzc+jo6LCqo7S0tNk9FYlEArFYTDMxRT1YWlqiqKiIbRmvkJyczMiRQGUpKipqtolFIhEA9pfImISaWINxc3PD48ePFeoyMkFNTQ2+/vprJCcnN7xWUlKCxMREbNmyhRVN9VRUVODx48dwd3dv1nP1Y2hNmGtgCmpiDcbb2xtSqZS1mWC5XI5jx46hW7du6Nu3L7744gvs378f586dg4mJCSua6rl+/TpkMlmzt4Bqo4lpUQANxs7ODs7Ozjh58iSGDBmi9vYNDAxw+/ZtlJeXQygUQl9fX+0a3sTJkyfh4uLS7MKApaWlAPDSLrSWDs3EGk5wcDD279/P6kywqampRhm4trYWBw4cUGg7aklJCfT19TVq1l9ZqIk1nODgYIhEIoSFhbEtRWPYs2cPqqqqMHPmzGY/W1JSolVZGKAm1nhsbW0REhKCtWvXoqqqim05rFNZWYl169Zh0aJFsLa2bvbziixLaTrUxC2AL774ApWVla2qON6bWL9+PWpqavDZZ58p9DzNxBRWaNu2LTZt2oSNGzciKiqKbTmsER0djc2bN2Pz5s2wtLRUKAbNxBTWCAkJQVBQEGbNmoVnz56xLUftZGVlYdq0aZgwYQLmzp2rcBxN2Y/OJNTELYjdu3fDzMwMI0aMaFgqaQ0UFxdjxIgRaNOmDXbt2qVULNqdprCKqakpzp8/j6qqKowcOVIjt2QyTWFhIUaOHAmxWIzz588rvcmEdqcprGNnZ4fIyEiUlJRgwIABePLkCduSVEZaWhq8vb1RXl6OixcvMnLjI83EFI3A0dERcXFxMDY2Rr9+/XD27Fm2JTHO6dOn0b9/f5ibmyMuLg4ODg5KxxSLxaiqqqKZmKIZWFlZ4cqVKwgICEBgYCCWL1/ecGVnS6ampgbLli3D2LFjERgYiMuXLzNWML9+37S2ZWKwfMk5hQH27t1LjIyMiIODAzl79izbchTm9OnTpHPnzsTY2Jjs27eP8fhJSUkEALl//z7jsdmEZmItYObMmUhNTUXv3r0bMvPt27fZltVkbt261aC7X79+SE1NxYwZMxhvRxtPMAG0O6012NnZ4eDBg4iMjERBQQF69+6NsWPHslrk/W1cvXoVgYGB6NOnD0pKShAVFYUDBw7A1tZWJe3Vm1jbLmWnJtYy/Pz8cOPGDZw5cwaFhYUYOHAgPDw8sG3bNoWKyjFNSUkJtm7dCnd3d/j6+qK4uBjnzp1DQkKCyo9blpaWwtjYGAKBQKXtqB22+/MU1XLz5k0yb948YmBgQPh8Phk2bBj56aefSFZWlto0ZGZmkh07dpBhw4YRPp9PDA0Nyfz588mtW7fUpoEQQr755hvSqVMntbapDuitiK2EiooKnD17FsePH0d4eDgqKyvRqVMn+Pr6YsCAAejWrRtcXV2VLiBXUVGBlJQU3Lt3D9euXUNMTAwyMjJgaGiI0aNHY8KECRg9ejQMDQ0Z+mZNZ+XKlYiKimK9ZjbTUBO3QmpraxEfH4+YmBhER0fjxo0bqKqqAofDQceOHdG+fXvY29vD2toabdq0gampKTgcDkxNTQG8uOqUEILy8nIUFRUhPz8fz549Q2ZmJjIzM0EIgaGhIfr27QtfX18MGjQI/fv3h66uLqvfe968ecjOzkZERASrOpiGmpgCQgiePn2K5ORkpKSkIDs7Gzk5OcjLy0NJSQmeP38OuVyO8vJyAC+2f3K5XJiYmMDS0hLW1tawt7eHvb09EhMTkZ6ejlu3boHD4bD8zV5m3Lhx0NfXx4EDB9iWwijUxBRGOX/+PEaNGoXU1FR07dqVbTkv0adPHwwePBibNm1iWwqj0NlpCqMMHz4c7dq1w759+9iW8gq5ubmws7NjWwbjUBNTGIXL5WLGjBnYs2cPZDIZ23IakMlkKCgoUNkaNJtQE1MYZ86cOcjNzdWoKiT5+fmQyWQ0E1MoTcHJyQn9+vVD2aKQQAAAEXlJREFUaGgo21IayM3NBQCaiSmUphIcHIzjx483zGizTU5ODjgcDiNnkjUNamKKSpg2bRoA4MiRIywreUFOTg4sLCxYX6tWBdTEFJVgYmKCsWPHakyXWltnpgFqYooKCQ4OxrVr1/DgwQO2pSA3N1crx8MANTFFhWjSmnFOTg7NxBRKc9GkNWOaiSkUBdGUNWNqYgpFQTRhzbi8vBxlZWXo1KkTaxpUCTUxReWwvWaclpYGAIyUvdVEqIkpKoftNeP09HTw+Xy0b9+elfZVDTUxReWwvWacnp6O9u3ba19trf9BTUxRC2yuGaenp2ttVxqgJqaoifo147CwMLW3TU1MoTAAl8vFu+++i9DQULWvGVMTUygMMXv2bOTm5uLSpUtqa7O2tha5ubnUxBQKEzg5OaF///5qneB68uQJ5HI5NTGFwhTqXjNOT08HAK3d6AFQE1PUzNSpU8HlcnH48GG1tJeeng4rKyuli+JrMtTEFLWi7jVjbZ/UAqiJKSwQHByMuLg4tawZP378GI6Ojipvh02oiSlqx9/fX21rxvfv34erq6vK22ETamKK2lHXmnFlZSWePXtGTUyhqAJ1rBnfv38fhBC4uLiorA1NgJqYwgrqWDO+f/8+dHV1tXp5CaAmprCIqteM6y914/F4KomvKVATU1hD1WvG9+/f1/quNEBNTGERVa8Zp6amUhNTKKqmsTVjkUikcNza2lpkZGRo/cw0QE1MYZl/rhkXFhbi+++/h7u7O9atW6dw3AcPHkAmk7UKE/PZFkBp3XC5XEybNg07d+7EvXv3EB4eDg6HA7lcjiFDhigc9/79+xAIBFq/WwugJqawSEpKCsLCwrBz506IRCKcP3++YfOHQCBAXV2dwrFTU1PRpUsXra2r9XeoiSlqZ9euXfj++++RmpoKoVAIiUQCAJBKpQ2fIYQ0vK4IKSkpraIrDdAxMYUFLCwsGiay3mRUuVyulInv3LmD7t27K/x8S4KamKJ2JkyYgK+++gpc7pv//AghCnenKyoqkJGRgW7duikqsUVBTUxhhc8++wxTpkwBn//6EZ0y3el79+6BEEIzMYWiSjgcDn777Tf06NEDQqHwtZ+pra1VKPbdu3dhamqqtTc+/BNqYgpr6Orq4vTp0zA3N39tRlYmE3fr1g0cDkdZiS0CamIKq1hbWyM8PBx8Pv8V04nFYoVi3r17t9WMhwFqYooG0KNHDxw4cOCV1xUxsVwuR3JyMjUxhaJugoKCsH79+peysSKz0+np6aisrGw1k1oANTFFg/jXv/6FqVOnNoyPFcnE9+7dA5fLhZubG9PyNBa6Y4uiMdTPWD9+/Bi3bt16KRPX1NQgNzcXJSUlEIlEkMvleP78OYAXRxq5XC6MjY1x5coVODg4wMDAgK2voXY4hBDCtggKpZ6CggKEh4fjo48+glwuh6OjI7KyslBWVtasOGZmZmjfvj1cXV3h4eEBNzc39OvXD1ZWVipSzh7UxBRWqaqqQnh4OCIjIxEdHY0HDx6Aw+HA1tYWIpEIixcvRocOHWBnZwdbW1u0adMGRkZG4HA4MDU1BQCUl5eDEAKRSISioiLk5eUhJycHmZmZSE5ORkpKCjIyMhqK5vn6+sLf3x8jR47UioxNTUxRO7W1tTh69CiOHDmCixcvQiKRoF+/fvD19cWgQYMwYMAAGBkZ4erVq/Dx8WGkTZFIhGvXriEmJgbR0dG4ceMGhEIhRowYgUmTJmHChAnQ0dFhpC11Q01MURtpaWn4+eefsWfPHohEIvj7+2P8+PEYM2YM2rRpo1YtRUVFOHXqFI4fP46LFy/CxMQEs2fPxvvvv9/yrn0hFIqKSU9PJyEhIYTP5xM7OzuycuVKkpWVxbasBvLy8sg333xDOnbsSLhcLpk0aRJ58OAB27KaDM3EFJVRXFyMVatWYc+ePXBycsLnn3/eUOFSE5HJZPjjjz/w73//G2lpaXjvvfewYcMGWFhYvPGZ2tpaFBQUoKioCMXFxaioqGh4z8TEBG3btoWFhQVsbGzeeNhDWaiJKYxDCMGvv/6KVatWQVdXF99++y2mTZumseb9JzKZDAcOHMCqVasgFovx7bff4r333sOTJ08QHx+PGzduIDU1FY8fP0ZWVhaaYiEdHR107doVbm5uGDhwIPz8/ODk5MSIXmpiCqMUFxdjzpw5OH/+PD788EOsW7euxd4NLBKJsHbtWmzbtg2dO3fGo0ePIBQK4enpCXd3d3Tp0gVdunSBnZ0dLCws0LZtWxgbGzc8X1ZWhuLiYpSUlCAtLQ3Jycm4d+8eYmNjUVFRgc6dO2PWrFmYPXs2OnTooLhQ9nryFG3j2rVrxM7OjnTo0IHExsayLYcxrl69SmxtbUnbtm1JdHS00vHq6upIbGwsWbFiBbGysiJcLpdMmTKF3L9/X6F41MQURjhx4gTR09MjgYGBpLS0lG05jFNSUkICAgKIvr4++fPPPxmLK5FIyOHDh4m7uzvh8Xhk0aJFpLKyslkxqIkpSvP7778THo9HFi5cSKRSKdtyVIZUKiUhISGEx+OR0NBQRmPLZDISGhpKzM3NiaOjI7l58//au/+YqOs/DuDPz915nXHJoZeDvDkJGSpIEUWHCeqU4m7Mabr5s24ROp1uzdaarVZ/teXWyrTVdG7qTKV9URhiUZLKKR5wiKOk0kPxQAQPKOTuFJPj+f3D5Te+HHXAffiovR/b/fP5fO71fjL2uvvc59fbGfJ7RRMLI1JYWEi1Ws13331X6SijZvPmzdRoNCwuLg577WvXrjErK4sREREsLS0N6T2iiYVhq6yspE6n4/r165WOMurWrl3LsWPHsrq6Ouy179y5Q5vNRq1Wy2+++eYftxdHp4Vh6e7uRkpKChISElBSUvLAnD4Kl0AgAKvVisuXL6O2tjbsR+BJIjc3F4WFhaiqqkJCQsKg24omFobFZrPh+++/R11dHSZOnKh0HEW0tbXhqaeegtVqxe7du8Ne//bt25g7dy58Ph/OnTs3+MUiYd8XEB56DoeDkiSxqKhI6SiKO3z4MCVJYlVVlSz1L126RJ1Ox88++2zQbcQ3sTBks2fPxpgxY3DixAlFxvd6vThw4AAaGxsxdepUrFy5Eo8++qgiWQAgMzMTkiShvLxclvrvvPMOdu7ciaampuC3Tsry8SE8tM6cOUMAdDgcioz/66+/Mjo6mvHx8dRqtQTAuLg4tra2KpKHvHsxCABZDnKRZEdHB3U6HXfv3h10vWhiYUjy8vI4c+ZMxca3WCysq6sjSXo8Hubl5REAc3NzFctEksnJyVy3bp1s9ZcvX87Zs2cHXSeaWAhZb28vDQYDP/nkE0XGr6mp4VdffdVv2bVr16hSqTht2jRFMv3p448/ZlRUFAOBgCz1Dx06RLVaze7u7gHr/l3nBYQR+emnn9DV1YWsrCxFxp8yZQpWrlzZb1lMTAxSU1MRFRWlSKY/zZ8/H7///jvq6+tlqZ+eno5AIACn0zlgnWhiIWQOhwORkZGKzfs7YcKEoFOzNDc3w2KxKJDof2bOnInHHnsMZ86ckaV+TEwMJk2ahNra2gHrRBMLIXO73Zg6dep9dWGH3W6HRqPBpk2bFM2hVqsRFxcHt9st2xhPPPEEPB7PgOX3z39DuO91dnb+7VMuRlsgEMD777+P4uJi6PV6pePAaDSis7NTtvpRUVFBH90rmlgI2c2bNzF27FilY9zz1ltv4c0330RKSorSUQAAERER8Pl8stb3+/0DlosmFkI22DeBEnbu3ImUlBQsXLhQ6Sj3/Pbbb7Luqfh8vqB7HKKJhZAZjUa0t7crHQOFhYUgiVdffbXfcrmumApVe3u7rE3c2dkZ9Ci8mItJCFliYiJcLteg3wijoaysDFu2bMHq1avx+eefA7j72/jnn39GUlIS5syZo0gur9cLl8uFpKQkWeqThMvlwvr16wesE00shGzWrFno7e2F0+nEvHnzRn382tpaLFq0CH6/H1VVVf3W6XQ6tLS0jHqmP1VVVSEQCGDWrFmy1He5XPB6vUE/JMTutBCySZMmISEhAUVFRYqM/8wzz8Dn84F3rzTs97p16xbGjx+vSC4AKCoqwvTp0xETEyNL/WPHjmHcuHF49tlnB6wTTSwMic1mw/79+4c1d/DDqqenBwcOHMBrr70m2xgHDx6ExWIJek+xaGJhSGw2G7q7u7Fv3z6lo9w39uzZA7/fj1deeUWW+nV1daioqMDrr78edL24n1gYso0bN6KwsBAXL158KKYGHQmfz4f4+HgsX74cn376qSxj5OTkoLW1FTU1NUEvOxV3MQlDdv36dY4bN46bN29WOori3n77bUZGRrK9vV2W+gUFBZQkiWVlZYNuI5pYGJYdO3ZQpVLxhx9+UDqKYk6ePEm1Ws1du3bJUr+hoYFGo5Fr1qz52+3E7rQwbEuXLkVlZSUqKythMpmUjjOqmpqaYDabkZGRga+//jrs9dva2jB37lzo9XqUl5f/7c8W0cTCsHV1dSEjIwN9fX04deqUoqd4RlNHRwcyMjKg1Wpht9sRGRkZdDuS+OWXX4Z862ZDQwMsFgvUajVOnDjxj6etxNFpYdgMBgNKS0vh9/uRnZ19X1ySKTePx4Ps7Gzcvn0bpaWl/Ro4EAigtrYWW7duxaJFixAVFYX09PQh1c/Pz0dqaioMBgPsdnto551l2ZkX/lVcLheffPJJxsfH89KlS0rHkY3L5WJcXBzj4uLY0NDAP/74gxUVFfzoo4/40ksvMSIiggA4ZswYqtVqAmBsbGxItc+fP0+LxUJJkrhhwwbeunUr5FyiiYWwaGtrY2pqKo1GI0tKSpSOE3bFxcWcMGECZ8yYwffee4/Z2dn3mlar1VKlUhHAgFdmZuagNfv6+lhWVsaXX36ZKpWKycnJtNvtQ84mmlgIG6/XS5vNRkmSuGnTJt68eVPpSCPW1NTEtLQ0ArjXqFqtlpIkBW3av75UKhVXrVrVr57H4+HBgweZm5tLk8lEAExPT2dBQcGwH7InDmwJYbdv3z5s2LABEydOxLZt22C1WpWONCwlJSV444034PF48PTTT8PhcECSJPT29ob0/kceeQRpaWmIj49HY2MjLl++jObmZqjVapjNZmRlZWHx4sUjv/Np5J9VgjDQ1atXuWzZMgJgTk4Oz549q3SkkDmdTlqtVgLgihUr2NLSQpJsbW3lBx98wIiICGo0mn/8JtZqtZw8eTJffPFFrlu3jlu2bOGRI0fo9XrDmlc0sSCrsrIyPvfcc5QkiQsXLuTp06eVjjQou93OnJwcAuDzzz/P48ePB93uxo0b3Lp1K41GI9Vq9aC71pIkMT8/X/bcoomFUXH06FGazWYCYFJSErdv386Ojg6lY7Gjo4Pbtm1jYmIiAdBsNoc0JzBJ9vT0cO/evYyNjaUkSUEPbpWXl8v8F4gmFkaZ0+lkXl7evV3SBQsW8Msvv2RTU9OoZXC73fziiy+4YMECajQa6vV6rlmzhjU1NcOqd+fOHe7fv//eB8Ffd7UvXLgQ5vQDiQNbgiK8Xi+OHj2Kw4cP49tvv4XP50NsbCwyMzPxwgsvIDk5GTNmzBjx5N1erxf19fX48ccfUVFRAbvdjitXrkCv18NqtWLJkiWwWq1hedwQSXz33Xf48MMPcfr0aQB3J2MP9wTk/080saC4np4eOBwO2O12lJeXo7q6Gn6/H5IkYcqUKZg8eTJMJhOio6Px+OOPw2AwQJIkGAwGAHcv/ySJrq4utLe3o62tDVevXoXb7Ybb7QZJ6PV6pKWlITMzE3PmzIHZbIZOp5Ptb6qursb27dtH5b5r0cTCfYckGhsbcf78edTX16O5uRktLS1obW1FZ2cnbty4gb6+PnR1dQG4e/mnSqVCZGQkjEYjoqOjYTKZYDKZkJiYiKSkJMTGxga/F/chIJpYEB5w4gYIQXjAiSYWhAecaGJBeMBpAPxH6RCCIAzffwH6TLwxwQu4RgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "make_graph_adj_random(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember that:\n",
    "\n",
    "$A_{ij} = 1$ if the node $j$ has an outgoing edge (arrow) going to node $i$ \n",
    "\n",
    "$A_{ij} = 0$ otherwise\n",
    "\n",
    "Hence the rows indicate incoming edges to a specific node and columns indicate outgoing edges from a specific node.\n",
    "\n",
    "Write the adjacency matrix A of the above graph:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "A = np.array([[0,0,0,0],\n",
    "              [0,1,1,0],\n",
    "              [0,1,0,1],\n",
    "              [1,1,0,1]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can call the helper function:\n",
    "```python    \n",
    "graph_matrix(mat, mat_label=None, show_weights=True, round_digits=3)\n",
    "# mat: 2d numpy array of shape (n,n) with the adjacency matrix\n",
    "# mat_label: 1d numpy array of shape (n,) with optional labels for the nodes\n",
    "# show_weights: boolean - option to display the weights of the edges\n",
    "# round_digits: integer - number of digits to display when showing the edge weights\n",
    "```\n",
    "to check if you get the same graph. Since in this example all the edge weights are zero, you should use `show_weights=False`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAI8AAAFbCAYAAADsu90eAAAABmJLR0QA/wD/AP+gvaeTAAAgAElE\nQVR4nO2deVyU19n+r2dWRllkQBgURJFNNiOKEY1oNGpc4E2MaawacUlU2sbgVtpqYiWNqSQWq8a3\nJrFxicWgvxA1alCrgriFCqkiVlHDIsg2YBhZZmDm/v1h4XVBhZl55swM8/18/MNnhvtcD16e/dyH\nIyKCDRudZ6+AtQIblovNPDb0xmYeG3ojYi3AnNFqtaioqEBFRQXu3r0LrVYLlUqFlpYWdOvWDVKp\nFDKZDD169ICHhwfkcjlrySbFZh4AjY2NyM7OxqVLl5CXl4crV67g5s2bqKyshFar7XAcOzs7eHp6\nIiAgACEhIQgODkZ4eDiCgoLAcRyPb8AGriuOtnQ6HS5cuIDDhw/j1KlTyM7Ohlqthlwub/tHDwgI\ngIeHB3r16gV3d3fI5XIIBAI4ODhAJBKhoaEBarUaTU1NqKmpQVlZGe7cuYOSkhLk5+fjypUruHr1\nKjQaDXr27ImRI0fixRdfRExMDPr06cP6V2AM9nYp85w5cwa7d+/G/v37UVZWhv79+2P06NEYNWoU\nRo0aZfR/1JaWFvz444/IzMxERkYGMjIyUFdXh8GDB+O1117D7Nmz0atXL6OWaUL2gqwclUpFn376\nKYWGhhIACgsLozVr1tClS5dMrkWtVtORI0dowYIF5OrqSiKRiF599VU6evSoybUYgVSrNY9KpaIN\nGzaQQqEgOzs7ev311+nYsWOsZbWhVqspNTWVXnrpJeI4jsLCwig1NZV0Oh1raR3F+syj1Wppy5Yt\n5OLiQo6OjrRq1SpSKpWsZT2VixcvUkxMDHEcR0OHDqUffviBtaSOYF3muXjxIkVERJBYLKbly5eb\nvWkeJScnh0aPHk0CgYAWLVpEtbW1rCU9Deswj06no3Xr1pFYLKaoqCjKy8tjLUlvdDod7dq1ixQK\nBXl7e9OZM2dYS3oSlm+empoaGj9+PInFYkpKSrKkPsNTqaqqosmTJ5NIJKKkpCTWctrDss1TXFxM\nwcHB5OXlZSn9hE6h0+lo/fr1JBQKKS4ujlpaWlhLehDLNU9BQQF5enpSaGgo3b59m7UcXklLSyOZ\nTEbTpk2j5uZm1nJasUzzlJWVkY+PDw0dOtTcO5VGIzMzk7p160bz5s0zl6bZ8sxTX19PAwcOpMDA\nQKqqqmItx6QcOnSIxGIxrV69mrUUIks0z6JFi8jZ2ZkKCwtZS2HC3/72NxIIBHTy5EnWUizLPAcO\nHCCO42jfvn2spTBl2rRp5OnpybrJTrWYhVGNRoOgoCAMGzYMX331FWs5TKmtrUVAQADmzJmDpKQk\nVjIsZ2E0OTmZZDIZFRcXMym/qamJ0tPTad26dXTmzBnmw+aNGzeSnZ0dy+Y71SK2oep0OiQnJyMu\nLg5eXl4mL7+yshIDBgxAcXEx5s2bh2+//Rb/8z//06mNYsZm4cKFcHNzw6effspMg0XUPOnp6QSA\nrl69avKytVotvfDCCxQTE9P2rKWlhby9vSkhIcHkeh5k9erV5O7uThqNhkXxltFhnjt3LkVGRjIp\n++TJkwSADh48+NDz999/n7p370737t1joouIqLCwkDiOo/T0dBbFW0azdfr0aYwbN45J2d988w0A\nIDQ09KHnISEhqK+vx+HDh1nIAgB4e3vD19cXWVlZTMo3e/NUV1fj5s2biIyMZFL+jRs3AAAeHh4P\nPXdzcwMAXL9+3eSaHmT48OE4d+4ck7LN3jxFRUUgIgQEBDApv6KiAkKhEBKJ5KHn3bp1AwDcuXOH\nhaw2/P39UVhYyKRsszdPdXU1AMDFxYVJ+fb29u0+bx1pKRQKU8p5DBcXFyiVSiZlm715GhsbAQAy\nmYxJ+V5eXtBqtVCr1Q89V6lUAICgoCAWstqwt7dHfX09k7LN3jzOzs4A7s+qsmDAgAEAgJKSkoee\nt9aIrM2jVCqZnVQ1e/O0NldVVVVMyp8/fz6kUinOnDnz0POLFy/iueeeg7+/PxNdrVRVVTFr0s3e\nPH5+frCzs0Nubi6T8hUKBX7zm9/g448/Bv13GbCpqQkHDx7Etm3bIBCw/RXm5OQ8No1gKsz+rLpU\nKsWgQYNw9uxZzJo1i4mGjz/+GCKRCDExMRg/fjzu3LmDVatWITw8nImeVogI58+fx+rVq5kJMHtW\nrVpFnp6ezBcjW1paqLy8nKmGBzl16hQBoMuXL7Mo3jJmmOfOnYvS0lIcO3aMqQ6hUAh3d3emGh5k\n27ZtiIiIQEhICJPyLcI8Pj4+iIqKQnJyMmspZsPt27exb98+vPXWW+xEsKjv9CEjI4MAsFoENDvm\nzJlDffv2paamJlYSLGNVvZXo6GgKDg6mxsZG1lKYcu7cORIKhbR7926WMixnGypwf6Ju4MCBmDVr\nFjZu3MhaDhPu3buH8PBw+Pj44MiRIywzjlnONtRWUlJSiOM42rNnD2spJker1dLUqVNJoVCYw6jP\nspqtVuLj40kqldLx48dZSzEpcXFxZGdnR5mZmaylEFmqebRaLc2cOZMcHBzoxIkTrOXwjk6no2XL\nlpFQKKS0tDTWclqxTPMQEWk0Gpo+fTpJpVJKSUlhLYc31Go1zZgxgyQSCesO8qNYrnmI7tdAS5cu\nJY7jKCEhgdVGcN4oLCyk4cOHk6Ojo1mlxPsvlm2eVrZt20bdu3enYcOG0Y0bN1jLMQp79+4lZ2dn\nCg4ONtdkVdZhHiKi/Px8GjhwIMlkMkpMTGQ5eWYQt27doilTphAAWrBgATU0NLCW9CSsxzxE9/tB\nSUlJZG9vT76+vrRr1y7mi6kdpbKykhISEkgmk1FQUJA5JDJ4FtZlnlZKSkooNjaWRCIRBQQE0Jdf\nfmm2s9LFxcW0YsUKsre3Jzc3N0pOTraUvpt1mqeVgoICmjNnDkkkEpLL5RQfH28W/QeNRkMHDx6k\n6OhoEgqFpFAo6JNPPqH6+nrW0jqDdZunlfLyclq7di3169ePAFBgYCCtXLmSsrOzTdas1dXV0f79\n+yk2NpbkcjlxHEdjx46l1NRUUqvVJtFgZCxrbctQdDodzpw5g2+++QZpaWkoKiqCk5MTXnjhBbzw\nwgsIDw9HSEiIwfdBtLS0oKCgAHl5eTh//jxOnz6N3Nxc6HQ6REZGYurUqZg6dSr69u1rnBdjQ9e6\nuORR8vLykJGRgczMTGRlZaGsrAwAIJfL4e/vD4VCAS8vL7i5ucHJyQlSqbTtnq3We7dUKhXq6upQ\nUlKCiooKFBcX4/r169BoNBCJRBgwYABGjRqFqKgoREVFmdVmMgPp2uZ5FKVSicuXL+PKlSu4ceMG\nysvLUVpaioqKCtTV1UGtVuPevXtobm6Gvb09xGIxHBwc4OjoiN69e0OhUMDT0xOBgYEIDg5GUFAQ\npFIp69fiC5t5OktqaireeOMN2H5tttuNbRiAzTw29MZmHht6YzOPDb2xmceG3tjMY0NvbOaxoTc2\n89jQG5t5bOiNzTw29MZmHht6YzOPDb2xmceG3tjMY0NvbOaxoTc289jQG5t5bOiNzTw29MZmHht6\nYzOPDb2xmceG3tjMY0NvbOaxoTc289jQG5t5bOiNzTw29MZmHht6YzOPDb2xmceG3tjMY0NvbOax\noTdmf0EtS+7du4dr16499OzWrVsA7l+N/SBisRhhYWEm02YO2JI7PYW6ujq4u7ujqanpmd995ZVX\nkJaWZgJVZoMtudPTcHR0xKRJkyASPbuC/uUvf2kCReaFzTzPYObMmdBqtU/9jkwmw5QpU0ykyHyw\nmecZTJ48Gd27d3/i52KxGNOmTUO3bt1MqMo8sJnnGUilUkybNg0SiaTdz5ubmzFjxgwTqzIPbObp\nADNmzIBGo2n3sx49emDs2LEmVmQe2MzTAcaMGQMXF5fHnovFYsyYMQNisZiBKvbYzNMBhEIhZs6c\n+VjT1dzc3CVHWa3YzNNBfvnLXz7WdHl4eGDEiBGMFLHHZp4OMmzYMPTp06ft72KxGLNnzwbHcQxV\nscVmnk4wa9astv5NV2+yAJt5OsUvf/lLNDc3AwB8fX0xcOBAxorYYjNPJwgJCUFgYCAAIDY2lrEa\n9thW1Z+CVqtFRUUFKioqcPfuXWi1WgwePBjXrl2DQqHA8ePHIZPJ0KNHD3h4eEAul7OWbFJsq+oA\nGhsbkZ2djUuXLiEvLw9XrlzBzZs3UVlZ+cx1rQexs7ODp6cnAgICEBISguDgYISHhyMoKMgaO9Zd\n874tnU6HCxcu4PDhwzh16hSys7OhVqshl8vb/tEDAgLg4eGBXr16wd3dHXK5HAKBAA4ODsjOzsbA\ngQOhVqvR1NSEmpoalJWV4c6dOygpKUF+fj6uXLmCq1evQqPRoGfPnhg5ciRefPFFxMTEPDRqs2C6\nlnnOnDmD3bt3Y//+/SgrK0P//v0xevRojBo1CqNGjTL6P2pLSwt+/PFHZGZmIiMjAxkZGairq8Pg\nwYPx2muvYfbs2QbfZ8qQvVZ/u7FKpaJPP/2UQkNDCQCFhYXRmjVr6NKlSybXolar6ciRI7RgwQJy\ndXUlkUhEr776Kh09etTkWoyA9V6NrVKpaMOGDaRQKMjOzo5ef/11OnbsGGtZbajVakpNTaWXXnqJ\nOI6jsLAwSk1NJZ1Ox1paR7G+q7F1Oh22bt2K9957D83NzVi8eDGWLFli1iOhnJwcrFmzBgcPHkRE\nRAQ2b96MiIiIJ35fp9OhvLwcSqUSSqUS1dXVbXeeymQyuLq6wsXFBb179+Zzn5F1NVsXL16kiIgI\nEovFtHz5clIqlawldYqcnBwaPXo0CQQCWrRoEdXW1lJNTQ0dPnyYVq9eTa+++ioFBweTVColAM/8\nw3Ec9e3blyZNmkSrV6+mzMxM0mg0xpJrHTUPEeHjjz/GqlWrEBkZiS1btiA4OJi1LL0gIuzevRsr\nVqwAx3EoLy8HEcHX1xfh4eHw8/ODn58ffHx84OLi0lbLCIVCAEBDQ0NbbVRcXNw28svKykJRURHs\n7e3x2muvYd68eRg5cqQhUwiWX/PU1NTQ+PHjSSwWU1JSkiX1GZ5KVVUVTZo0iYRCIa1evdooMQsK\nCmjjxo00ePBgAkChoaH0zTff6Ps7s+wOc3FxMQUHB5OXlxf98MMPrOUYHZ1OR+vXryehUEhxcXHU\n0tJitNg//vgj/eIXvyCO4+j555+n/Pz8zoawXPMUFBSQp6cnhYaG0u3bt1nL4ZW0tDSSyWQ0bdo0\nam5uNmrs3Nxcev7550kmk9HGjRs786OWaZ6ysjLy8fGhoUOHUm1tLWs5JiEzM5O6detG8+bNM3rT\n3NzcTImJiSQQCCg+Pr6j8S3PPPX19TRw4EAKDAykqqoq1nJMyqFDh0gsFhutD/QoX3/9NUmlUlq4\ncGFHvm555lm0aBE5OztTYWEhaylM+Nvf/kYCgYBOnjzJS/wDBw6QQCCgTZs2PeurlmWeAwcOEMdx\ntG/fPtZSmDJt2jTy9PTkrcn+8MMPSSQSUU5OztO+ZjnmUavV1L9/f5o5cyZrKcypqamhnj170ooV\nK3iJr9PpaOTIkTR8+PCn9X8sxzzJyckkk8mouLiYSfm1tbX0ySef0OLFiyk9Pd2ow2Z92LhxI9nZ\n2fHWfGdnZ5NAIKBvv/32SV+xDPNotVrq06cPLV26lEn5SqWS+vfvT2+++SaNGTOGBAIBDR06lImW\nVtRqNfXp04e32oeIaOLEiTRp0qQnfWwZ5klPTycAdPXqVSbl/+///u9D62SJiYkEgLKyspjoaWX1\n6tXk7u5uzPWqh9i3bx8JhcInzaNZhnnmzp1LkZGRTMpWq9V069ath54VFhYSACZ7gh7VwXEcpaen\n8xJfrVaTTCaj7du3t/dxqkWcnjh9+jTGjRvHpGyJRIJ+/fo99OzSpUuYMmUKQkNDmWhqxdvbG76+\nvsjKyuIlvkQiQXh4OC5cuNDu52Z/eqK6uho3b95EZGQkaykgIuzduxdr1qxBeno6azkAgOHDh+Pc\nuXO8xY+IiMDZs2fb/czsa56ioiIQEQICApjqqK+vx8KFCzF37lzk5+cjNDQU2dnZTDUBgL+/PwoL\nC3mL7+HhgcrKynY/M3vzVFdXA0C7KU5MSffu3fHZZ59BpVIhOTkZKpUKcXFxTDUB938vSqWSt/hy\nuRy1tbXtfmb25mlsbARwf3ulOSAQCBAfH4+pU6ciNzcXarWaqR57e3vU19fzFr9bt25PjG/25nF2\ndgaAJ7qfFePGjYNcLodUKmWqQ6lU8ro/+969e7C3t2/3M7M3T2tzVVVVxVjJw+Tl5SE6Opq1DFRV\nVfHapCuVyrb/wI9i9ubx8/ODnZ0dcnNzmZTf2NiIDz/8EHl5eW3PlEolcnNzkZyczETTg+Tk5PA6\nZXD9+nX4+fm1+5nZm0cqlWLQoEFPHC7yjU6nw//7f/8PYWFhGDp0KN5//33s3r0bhw8fhpOTExNN\nrRARzp8/z+s0Rk5ODkJCQp4owOxZtWoVeXp6Ml2MrK2tpfr6emblt8epU6cIAF2+fJmX+OXl5cRx\nHH333XftfWwZM8xz585FaWkpjh07xkxDjx49zC5R97Zt2xAREfHkmsFA9uzZA3t7e4wePbrdzy3C\nPD4+PoiKijKLPoa5cPv2bezbtw9vvfUWL/G1Wi22bt2K6dOnPzkDPi/1HQ9kZGQQAN4WAS2NOXPm\nUN++fampqYmX+Nu3byeRSPS0nQyWsareSnR0NAUHB1NjYyNrKUw5d+4cCYVC2r17Ny/xKysryd3d\nnRYtWvS0r1mWeYqLi8nZ2Zneeecd1lKYoVKpyM/PjyZMmMDL6diWlhaaOHEieXt7088///y0r1qW\neYiIUlJSiOM42rNnD2spJker1dLUqVNJoVBQeXk5L/HffvttkslkdP78+Wd93fLMQ0QUHx9PUqmU\njh8/zlqKSYmLiyM7OzvKzMx86vcqKys7faatqamJZs2aRRKJhPbv39+RH7FM82i1Wpo5cyY5ODjQ\niRMnWMvhHZ1OR8uWLSOhUEhpaWmPfV5SUkJfffUVLVy4kPr3708cx3Uq29i1a9do8ODB5Ojo2JkB\niWWah4hIo9HQ9OnTSSqVUkpKCms5vKFWq2nGjBkkkUjaOsjXr1+nbdu2UWxsLHl5eREAEggED+Xt\nycvLe2bsn3/+md577z2SSqU0ePBgKigo6Iw0yzUP0f0aaOnSpcRxHCUkJPC2EZwVhYWFNHz4cLK3\nt6f4+HiaNWsWeXh4EAASCoUkEomemNippqbmiXFv3bpFv//978nZ2Zl69OhBycnJ+iRQsGzztLJt\n2zbq3r07DRs2jG7cuMFajsG0tLTQwoULSSwWk0AgIAAkEolIKBR2KCOYRCJ5aCSmVqvp1KlTtHLl\nSho6dCgJBALy8PCgxMREQ06dWod5iIjy8/Np4MCBJJPJKDExkbfJM765desWTZkyhQBQREQEyeVy\n4jiuzUQd+ePi4kLvvvsuxcTEUEhICMlkMgJAvr6+FBcXRwcOHDBGqhbrMQ/R/X5QUlIS2dvbk6+v\nL+3atYv5yc6OUllZSQkJCSSTySgoKKgtkYFaraYdO3aQj48PcRzXodqnV69eFBkZSTNnzqRVq1bR\n9u3b6aeffjK2ZOsyTyslJSUUGxtLIpGIAgIC6MsvvzTbWeni4mJasWIF2dvbk5ubGyUnJ7fbd9Nq\ntXTgwAEaNGhQWzPWnnE4jqNp06aZQrp1mqeVgoICmjNnDkkkEpLL5RQfH9+hUQjfaDQaOnjwIEVH\nR5NQKCSFQkGffPJJh7d8nD59ml5++WXiOI7EYvFj/Z13332X5zcgIms3Tyvl5eW0du1a6tevHwGg\nwMBAWrlyJWVnZxu9WXvSkkFdXR3t37+fYmNj2/oxY8eOpdTUVFKr1XqVlZ2dTVOnTn3IRBKJhD76\n6CNDXqGjdA3ztKLVaikzM5Pi4+PJ29ubAJCTkxNNnjyZPvroI0pPT6fS0lKDykhOTqa7d+9Sfn4+\npaam0tKlSykiIoJEIhEJBAIaMWIErV+/3qh9kOvXr9Pbb7/dZqAnHA82NtaRh1lf8vLykJGRgczM\nTGRlZaGsrAzA/bNK/v7+UCgU8PLygpubG5ycnCCVStGtWzdIpVKoVCq0tLRApVKhrq4OJSUlqKio\nwJkzZ3Dv3j1otVqIRCIMGDAAo0aNQlRUFKKiouDu7s7b+5SXl2PDhg2Ijo42xcW5XevWm2fxpz/9\nCV988QWWL1+Omzdvory8HKWlpaioqEBdXR3UajXq6+uh0Whgb28PsVgMBwcHODo6onfv3lAoFPjn\nP/+J0tJS7NixA7/4xS+YH83hEZt5HuS5557D0KFD8dlnn+kdw8fHBz/99BP8/Pxw+fJlqzaPRWxD\nNQU5OTn497//bfDdoa3ny27duoUPPvjAGNLMFlvN81/effddHDlyBNeuXdP7PoampqaHjkULBAJc\nuHABQ4YMMZZMc8JW8wCARqNBSkoKYmNjDboLtKKi4qG/CwQCzJ49GxqNxlCJZonNPAC+++47KJVK\nvPnmmwbFedQ8LS0tuH79Ov785z8bFNdcsZkHwI4dOzB27FiD7xgtLy9/7JlWq0ViYiKz49J80uXN\nU1lZiSNHjhjcUQbu1zwi0ePJ1jiOQ2xsLJqbmw0uw5zo8ub56quvIJPJ8Oqrrxocq7y8vO3StAdp\naWlBfn4+1q9fb3AZ5kSXN8/OnTvxxhtvGOUocUVFBXQ6XbufabVavP/++8jPzze4HHOhS5vHWHM7\nrZSXl6OlpeWJnxMRZs+eDa1Wa5TyWNOlzbNjxw74+flh+PDhRol3+/ZtPG3ajIhw8eJFbNq0ySjl\nscbsU+nyRXNzM1JSUvDuu+8aNLfzII+Otlr7P1qtFnK5HOPGjcOYMWMwZswYo5THmi5rHmPN7TxI\na+ZWAHB1dUVYWBhOnz6N8+fPW+Usc5c1z44dOzBmzBiD53ZaaW5uRnR0NEaPHo1Ro0YhMDAQ5eXl\n8PDwQF1dnVHKMDe6pHkqKytx+PBh/P3vfzdaTLFYjJSUlIeeKRQK9O7dG9nZ2VbTVD1Il+ww7969\nGzKZDFOnTuW9rIiICLPIFM8HXdI8rRu1TJEmbsiQIfjXv/7Fezks6HLmyc3NNerczrOIiIhAUVHR\nE+9vsGS6nHl27NgBX19fU+zxBXC/5uE4ziprny5lnubmZvzjH/8weN9OZ5DL5fDx8bHKfk+XMk/r\n3M7s2bNNWu6QIUNs5rF0jD2301GGDBmCnJwck5ZpCrqMeVrndkzVUX6QkJAQ3Llzh9d7sVjQZcxj\nyrmdRwkODgYAXL161eRl80mXMY8p53YexdPTE05OTrhy5YrJy+aTLmEeU8/tPArHcQgMDLSZxxIx\n9dxOewQHB1vVLkKgC5iHxdxOewQHB9tqHkuD1dzOowQHB6O8vPyhPT+WjtWbh9XczqMEBQUBsK4R\nl1Wbh+XczqNY44jLqs3Dcm7nUTiOw4ABA/Cf//yHtRSjYRXmaWxsbPc0Jsu5nfbo27cvCgsLWcsw\nGlZhnh9//BG9evXCkiVL8O9//xsA+7md9vD29kZRURFrGUbDKvYwNzc3o7q6Gps3b8aGDRsQHBwM\nV1dX9OvXj+nczqNYm3msouZRq9UA0HZaMz8/H1lZWSgsLMTYsWOxd+9es8iR4+3tjdraWvz888+s\npRgFqzDPo/0dIoJWqwUR4fTp03jjjTfg7u6OJUuW4O7du4xU3u/zAEBxcTEzDcbEKs3zIC0tLSAi\n3L17F2VlZXBycjKhsofx9vYGAKtpuqyiz9PabD0JsViMoKAgfPnll0yXKLp37w4XFxerMY/V1DwC\nQfuvIhKJ0KNHDxw6dMgshux9+/a1mcec0Gg07dYoHMdBIBDg0KFD6N27NwNlj+Pt7W01cz1WbR7g\nfvKmiIgIEyt6Mraax8xobm5+zDwCgQAffPAB3njjDUaq2sfDwwN37txhLcMoWEWH+dGaRywWIyYm\nBn/4wx8YqmofNze3tizxlo5V1DwPmqd1ZLVz506mI6sn0bNnTzQ1NUGlUrGWYjBWYZ7WuRyhUAgX\nFxd8//33ZjGyag83NzcAsIqz62bVbGm1WlRUVKCiogJ3796FVqttu9eq9Z4rmUyGHj16wMPDA3K5\nHMD9mkej0UAikeDAgQNQKBSM3+TJtJqnqqoK/fv3Z6zGMJiYp7GxEdnZ2bh06RLy8vJw5coV3Lx5\nE5WVlZ3KFGpnZwdPT08QETiOQ1xcHLp169b2d3PExcUFAKxiO6pJbr3R6XS4cOECDh8+jFOnTiE7\nOxtqtRpyuRwhISEIDg5GQEAAPDw80KtXL7i7u0Mul0MgEMDBwQEikQgNDQ1Qq9VoampCTU0NysrK\ncOfOHZSUlOCrr77Czz//DKVSCY1Gg549e2LkyJF48cUXERMTw3wL6qNIJBJs374dM2bMYC3FEPby\nesdoVlYWxcXFUa9evQgA9e/fn+bPn087d+6koqIio5WTmZlJOp2OmpubKTs7m9avX08xMTHk5ORE\nHMfRkCFD6KOPPjL4/lBj4eLiQlu2bGEtw1CMf0GtSqWiTz/9lEJDQwkAhYWF0Zo1a+jSpUvGLuqZ\nqNVqOnLkCC1YsIBcXV1JJBLRq6++SkePHjW5lgfp168f/fnPf2aqwQgYzzwqlYo2bNhACoWC7Ozs\n6PXXX6djx44ZK7zBqNVqSk1NpZdeeok4jqOwsDBKTU194lXWfPLcc8/R73//e5OXa2QMN49Wq6Ut\nW7aQi4sLOTo60qpVq0ipVBpDHG9cvHiRYmJiiOM4Gjp0KP3www8mLX/UqFH0q1/9yqRl8oBh5rl4\n8SJFRESQWCym5cuXm71pHiUnJ4dGjx5NAoGAFi1aRLW1tSYpd/z48TR//nyTlMUjqXpNEhIRkpKS\nMGzYMMhkMuTm5uLjjz9um3exFAYNGoQTJ05gx44d+Pbbb/Hcc8/h7NmzvJcrlUqfuQfJIuis3Wpq\namj8+PEkFospKSmJSZ+BD6qqqmjy5MkkEokoKSmJ17KmTZtGv/jFL3gtw/wRNHkAABsZSURBVASk\ndmqSsKSkBBMnTkRdXR3OnDljVlsdDMXV1RUHDx5EcnIyfvvb3+Knn37Cpk2b2r18zVAkEgnq6+uN\nHtfUdNg8N27cwIsvvghnZ2ecO3fObDZXGROO47B06VL4+PhgxowZqKqqQkpKSrtXPxqCVCpFTU2N\nUWOyoEN9njt37mDChAno1asXMjMzrdI4D/LKK68gPT0dhw8fxsKFC596h5Y+CASCJ94IaEk80zwN\nDQ2YOHEiJBIJDh06hB49ephCF3NGjhyJvXv3YteuXVizZo1RY5MZr711hmeaZ9myZSguLsb3338P\nV1dXU2gyGyZNmoRNmzbhgw8+wKlTp4wWt0uY5+DBg9i6dSs+//zztjNHXY2FCxdi6tSpePPNN412\nYJCInnjaw5J44htoNBosWbIEM2bMwGuvvWZKTWbHZ599BrVajbVr1xolnk6ns+6aZ8uWLSgrK8NH\nH31kSj1tqFQqbN26Fb/73e/wxRdfoKGhgYkOAHB2dsZ7772HTZs2GeXkg06ns4qap91JQq1WS336\n9KGlS5eactKpjf/85z+kUCjIz8+PJBJJ23aOO3fuMNFDdH9htU+fPrRixQqDY02dOpXeeOMNI6hi\nSvvLE8ePH0dxcTHefvttE1v5PkuWLEF6ejquX7+O27dv46233sLNmzexcuVKJnqA+xN7c+fOxc6d\nO596Nr4jNDQ0oHv37kZSxo52zbNnzx5ERkYiMDDQ1Hpw8eJFzJw5E2FhYQDunzZITEyEQCAwybrT\n05g7dy4qKytx8uRJg+I0NDSY7Qb9ztCueU6fPo1x48aZWguA+ycqH92e6eHhgcGDB8PZ2ZmJpla8\nvb3h6+uLrKwsg+JYbc1TXV2NmzdvIjIykoUeuLi4tDsSaV1XY83w4cNx7tw5g2I0NDRAJpMZSRE7\nHjNPUVERiAgBAQEs9LRLZmYmRCIRlixZwloK/P39DU5UUFNTw7wWNQbt1jzA/x0RYY1Wq8X777+P\nAwcOwN7enrUcuLi4GHRvFhFBqVRaxWz9Y+ZpbGwEALOpVpcvX46lS5di0KBBrKUAAOzt7Q3aTvHz\nzz+jubnZOs3TWp3W1taaXMyjfPbZZxg0aBBiYmJYS2lDqVQatGOytWa3SvO0NlesMzmkpaWBiB67\ncCQjI4ORovtUVVUZ1KRbk3ke2+Xk5+cHOzs75Obmtl1vaGqOHz+OdevWYdasWdi8eTOA+32f/Px8\nhISEYNSoUUx0AUBOTg5CQ0P1/vnKykpwHGed5pFKpRg0aBDOnj2LWbNmmVxQTk4OXnnlFdTX1+PC\nhQsPfWZnZ4fS0lKTa2qFiHD+/HmsXr1a7xhFRUVwdXW13knCsWPH4uDBg51KOmAswsPDce/ePRDR\nY38aGxuZntDIzMxETU0NxowZo3eMoqIiq9ne0q555s6di9LSUhw7dszUesyabdu2ISIiAiEhIXrH\nsHrz+Pj4ICoqCsnJyabWY7bcvn0b+/btw1tvvWVQHKs3DwAkJibi6NGjOHr0qCn1mC3vvfce3N3d\nDb5Fx5rM89RDf9HR0RQcHEyNjY08bw0xb86dO0dCoZB2795tUJzq6moCQOnp6UZSxpSnn1UvLi4m\nZ2dneuedd0wlyOxQqVTk5+dHEyZMMPh0bEZGBgGgkpISI6ljytPPqnt5eWHLli3YvHkzvv76a9NU\nhWaETqdDbGwsVCoVduzYYfC+4/z8fDg6OlrNubdnHoWcPn06Lly4gNjYWLi6umLs2LGm0GUW/OY3\nv8Hhw4dx9OhRuLu7GxwvPz8fQUFBVrH5HUDHEh1otVqaOXMmOTg40IkTJ/iuDpmj0+lo2bJlJBQK\nKS0tzWhxx4wZQ/PmzTNaPMZ0PD+PRqOh6dOnk1QqpZSUFD5FMUWtVtOMGTNIIpEY3EF+EJ1OR87O\nzrRx40ajxWRM55I7abVaWrp0KXEcRwkJCaTRaPgSxoTCwkIaPnw4OTo6Gj0l3tWrVwmAybOQ8Yh+\nmcG2bdtG3bt3p2HDhtGNGzeMLYoJe/fuJWdnZwoODqa8vDyjx//yyy/Jzs6O1Gq10WMzQv+0cvn5\n+TRw4ECSyWSUmJhITU1NxhRmMm7dukVTpkwhALRgwQJqaGjgpZxFixbR8OHDeYnNCMNyEmo0GkpK\nSiJ7e3vy9fWlXbt2UUtLi7HE8UplZSUlJCSQTCajoKAgOnnyJK/lDRo0iJYsWcJrGSbGOKl0S0pK\nKDY2lkQiEQUEBNCXX35ptrPSxcXFtGLFCrK3tyc3NzdKTk7mve9WW1tLIpGIUlNTeS3HxBg3iXdB\nQQHNmTOHJBIJyeVyio+P56X/0Fk0Gg0dPHiQoqOjSSgUkkKhoE8++YTq6+tNUv4333xDQqGQqqur\nTVKeiTB+BngiovLyclq7di3169ePAFBgYCCtXLmSsrOzTdas1dXV0f79+yk2NpbkcjlxHEdjx46l\n1NRUk3da4+LiaOjQoSYt0wSk8npxiU6nw5kzZ/DNN98gLS0NRUVFcHJywgsvvIAXXngB4eHhCAkJ\nQa9evfQuo66uDufPn4eXlxfy8vJw/vx5nD59Grm5udDpdIiMjMTUqVMxdepU9O3b13gv1wl8fX0x\nffp0/OlPf2JSPk/sNcmtN63k5eUhIyMDmZmZyMrKQllZGQBALpfD398fCoUCXl5ecHNzg5OTE6RS\nads9W633bqlUKtTV1aGkpAQVFRXIy8vD7du3odPpIBKJMGDAAIwaNQpRUVGIiooyyrKCIRQWFqJf\nv344deoU073XPLDXpPdthYSEICQkBL/+9a8B3D/GsnPnTqxatQqDBg1CbW0tfvzxR1RUVKCurg5q\ntRr19fXQaDSwt7eHWCyGg4ND2+KiQqGAq6srSkpKsH//fkyYMAFSqdSUr/RMvvvuOzg6OjI7vs0n\nTG/6c3FxQa9evaBWq7F582a9Eh55eXmBiHDy5EmzOt/Vyr59+zBlyhRIJBLWUowO8/RU1dXVbRez\ndZaCggLcvn0bAPD555/j3r17xpZnEBUVFcjKyrLatHxmYR59zzClp6e3JdhuamrCjh07jCnNYNLS\n0mBnZ4eXX36ZtRResGjzfP/9920JtnU6HT755BOjJ9w2hK+++grR0dFWcUarPczCPD179uz0zzU3\nN+PkyZNtZ8uICIWFhTh+/LixJerFjRs3cPbsWYM3zJszZmEefWqes2fPPpYhVSQSYcOGDcaSZhC7\ndu2Cm5sbXnrpJdZSeMNizZOenv7YCKalpQVHjhzB9evXjSVPL1paWrB9+3bMnj3b6JeemBPMzVNb\nW6vXEeJDhw5Bo9E89lwkEmHr1q3GkKY3+/fvx+3bt7Fo0SKmOviGuXlUKlWnM35VV1fj8uXL7X7W\n3NyMzz77jOmwfdOmTZgyZQp8fHyYaTAFzM2jT1rZZ51ibWxsxK5duwyRpTd5eXnIzMzEO++8w6R8\nU8LUPFqtFk1NTZ1OK/vg/E576HQ6/OUvf2EybN+wYQMCAwO7xBElpuZpHS11puYhIhw+fPipWdiJ\nCDdu3MA///lPgzV2htraWqSkpOCdd96xnrNZT4HpUKA1MWRnap7Lly+3pWZ7lO7du8PJyQnOzs5w\ndXXFtWvXTDpU3rp1KyQSCd58802TlckSizOPnZ0dUlJS4OzsDLlcDmdnZ2zZsgUnTpzAjz/+yJfU\nZ9LU1ITNmzdj7ty5ZpHy1xQwNY8+zZa/vz/8/f0feubo6GjwZSKG8vnnn0OpVGLZsmVMdZgSpn0e\nfWqe9pBIJO3O+ZiKpqYmJCUlYeHChVaTxKAjMDVPa8JwOzs7g+KIxWKm5vnss89QXV2NFStWMNPA\nAuZDdQAGX3wvkUiYNVuttU5cXFyXqnUAxuZpvVvc0CsTWTZbW7duhVKpxPLly5mUzxKrMA+rZqux\nsRHr1q3Dr371K4NOgFgqtmbLADZs2ACVSoXf/va3Ji/bHLCamsfU5qmoqMCf//xnJCQkMD/ewwqr\nMA+Lq6b/8Ic/wNHREUuXLjVpueYE00lCSzXPv//9b+zYsQO7du2y2v3JHcEqah6tVmtwv6kzLF++\nHBEREZg+fbrJyjRHmNY8rVsmDF2BNmXNk5aWhn/+8584ffp0l1g5fxpMax6xWAwABnd2TWUetVqN\n3/72t5g+fTpGjBjBe3nmDlPztG5gN3SOxlTN1tq1a1FeXo5169bxXpYlwLTZMpZ5TFHzXLt2DevW\nrcO6devg5eXFa1mWgq3m6QBEhLi4OISGhuI3v/kNb+VYGkxrntZ0KGq12qA4fNc8n3/+OTIzM/HD\nDz+YdFRn7lhFzaPPCYyOUlFRgd/97neIj49HeHg4L2VYKlZhnvr6eoM3lD2Jd999F46OjlizZg0v\n8S0Zps1W6zzJtWvXUFdXh5qaGtTW1qKmpga/+tWv2obyz4Iv86SlpeHrr7/GoUOHeDOnRWOq1Jkt\nLS0UExNDgYGB5O7uTjKZjAC0+ycgIKBTsRcsWEBjx441qt7Kykpyd3en+fPnGzWuFfH0y9qMiVAo\nRL9+/XDt2jVUVFS0bUF9FIlE0ulkSA0NDUavGeLi4iCTyfCXv/zFqHGtCZP2eTpyBFej0WDChAmd\nimvsZmvHjh1IS0vD3//+dzg6OhotrrVhUvP0798fEyZMeGpfRiQSISoqqlNxjWme0tJSLFmyBIsX\nL8aLL75olJjWislHW/Hx8U9cy+I4DiNGjOi0EYxlHiLCW2+9BXd3d6xdu9bgeNaOyUdb48ePh6+v\nL27evPlYIgKRSIRJkyZ1Oua9e/eMckpzw4YNOH78OM6cOQOZTGZwPGvH5DUPx3GIj49vd0a4ubkZ\n48eP73TM6upquLi4GKTr4sWL+N3vfofExEQMHTrUoFhdBhZjPJVKRd27d39siC6Xy/W6u1wqldKu\nXbsM0hMQEECjR4+2mPvCzADTDdUfxN7eHgsWLHio4ywSiTBx4sROb7BSqVRQq9V6p+MF7g/La2tr\n8Y9//MO2dtUJmC1PLF68uO3oDXB/cbOzQ3QAqKqqAgC9zfP3v/8du3fvxrZt2+Dh4aFXjK4KM/P0\n7dsXL7/8clvtQ0R65dJpzdWjj3muXr2KxYsXY8WKFZgyZUqnf76rw3Rh9MFhe0BAgF7/81trns4m\nAr937x5ee+01hIWFWds9WCaD6cLouHHjEBAQgGvXriE6OhparRYVFRWoqKjA3bt3odVq2+7Zar13\nSyaToUePHvDw8IBcLkd1dTXs7Ow6Pc8zf/58KJVKHD16tMMLsDYegUU3vaGhgTIyMmjTpk00cuTI\ntpGWUCh84mJpe3/s7OzI3d2dnJycKCEhgXbu3El5eXnPHLGtX7+eBAIBpaenm+iNrRJ+r4lsRafT\n4cKFCzh8+DBOnTqF7OxsqNVqyOVyBAUFITc3F3/84x/Rp08f9OrVC+7u7m3XKDk4OEAkEqGhoQFq\ntRpNTU2oqalBWVkZ7ty5g5KSEuTn5+PKlSu4evUqNBoNevbsiZEjR+LFF19ETEwM+vTp06bl3Llz\nGD16NBITE5GQkMD3q1sze3mtebKysiguLo569epFAKh///40f/582rlzJxUVFbV97/vvvzdKec3N\nzZSdnU3r16+nmJgYcnJyIo7jaMiQIfTRRx9Rbm4ueXh40CuvvKLXfJKNhzD+7cYqlYo+/fRTCg0N\nJQAUFhZGa9asoUuXLhm7qGeiVqvpyJEjtGDBAnJ1dSWRSER+fn6UlpZmci1WiPHMo1KpaMOGDaRQ\nKMjOzo5ef/11OnbsmLHCG4xarabU1FR66aWXiOM4CgsLo9TUVFsNpD+G93l0Oh22bt2K9957D83N\nzVi8eDGWLFmi12UkpiInJwdr1qzBwYMHERERgc2bNyMiIuKJ39fpdCgvL4dSqYRSqUR1dXXboq5M\nJoOrqytcXFzQu3fvrpT4wLA+z8WLFykiIoLEYjEtX76clEqlUSxtKnJycmj06NEkEAho0aJFVFtb\nS0qlkg4fPkzvv/8+vfLKKxQUFERSqbRDoz+O46hfv340efJk+uMf/0inT58mjUbD+jX5Qr+ah4jw\n8ccfY9WqVYiMjMSWLVsQHBxsRFObDiLC7t27sWLFCnAch/LychAR/Pz8MGjQIPj5+cHf3x/9+vWD\ni4tLWy3TugbW0NDQVhsVFxcjPz8feXl5yMrKQnFxMRwcHPDaa69h3rx5eOGFF6wpOULna56amhoa\nP348icViSkpKspo+Q1VVFU2aNImEQiGtXr3aKDGvX79Of/3rXyk8PJwAUGhoKKWlpVnL76xzHebi\n4mIKDg4mLy8v+uGHH/gSxQydTkfr168noVBIcXFxRt2ekZubS6+//jpxHEfPP/885efnGy02Izpu\nnoKCAvL09KTQ0FC6ffs2n6KYk5aWRjKZjKZNm0bNzc1GjZ2bm0vPP/88yWQy2rhxo1Fjm5iOmaes\nrIx8fHxo6NChVFtby7cosyAzM5O6detG8+bNM3oz09zcTImJiSQQCCg+Pt5Sm7Fnm6e+vp4GDhxI\ngYGBVFVVZQpRZsOhQ4dILBYbrQ/0KF9//TVJpVJauHAhL/F55tnmWbRoETk7O1NhYaEpBJkdf/vb\n30ggENDJkyd5iX/gwAESCAS0adMmXuLzyNPNc+DAAeI4jvbt22cqQWbJtGnTyNPTk7cm+8MPPySR\nSEQ5OTm8xOeJJ5tHrVZT//79aebMmaYUZJbU1NRQz549acWKFbzE1+l0NHLkSBo+fLgl9X+ebJ7k\n5GSSyWRUXFxsSkGPcefOHd6ajM6wceNGsrOz4635zs7OJoFAQN9++y0v8XmgffNotVrq06cPLV26\n1NSC2qisrKRly5aRTCajxYsXM9PRilqtpj59+vBW+xARTZw4kSZNmsRbfCPT/tGb48ePo7i4GG+/\n/baJZ7z/j8LCQsyePfuJ2TRMjUQiwdy5c7Fz507e7rmYP38+0tPTUVpaykt8Y9Ouefbs2YPIyEgE\nBgaaWk8bERERTMtvj7lz56KyshInT57kJX50dDQkEgmOHz/OS3xj0655Tp8+jXHjxplai9nj7e0N\nX19fZGVl8RJfIpEgPDwcFy5c4CW+sXnMPNXV1bh58yYiIyNZ6DF7hg8fjnPnzvEWPyIiAhcvXuQt\nvjF5zDxFRUUgIgQEBLDQY/b4+/ujsLCQt/geHh6orKzkLb4xabfmAWBw1glrxcXFBUqlkrf4crkc\ntbW1vMU3Jo+Zp3V0Y8tP0z729vZt98HzQbdu3XiNb0weM4+zszMAWIz7TY1SqeR1f7axElWZgsfM\n09pctZ4Bt/EwVVVVvDbpSqWy7T+wufOYefz8/GBnZ4fc3FwWeh6itfZrampirOT/yMnJQWhoKG/x\nr1+/Dj8/P97iG5PHzCOVSjFo0CCcPXuWhZ42jhw5gnfffRcA8O233+KLL75AeXk5U01EhPPnz/M6\njZGTk4OQkBDe4huV9hYtVq1aRZ6enrYUa49w6tQpAkCXL1/mJX55eTlxHEffffcdL/GNTPtrW3Pn\nzkVpaSmOHTtmYiubN9u2bUNERARvNcOePXtgb2+P0aNH8xLf2LRrHh8fH0RFRSE5OdnUesyW27dv\nY9++fXjrrbd4ia/VarF161ZMnz7dci5JeVKdlJGRQQBsOWz+y5w5c6hv377U1NTES/zt27eTSCSi\nq1ev8hKfB56+DTU6OpqCg4OpsbHRVILMknPnzpFQKKTdu3fzEr/1hp1FixbxEp8nnm6e4uJicnZ2\npnfeecdUgswOlUpFfn5+NGHCBF62iLa0tNDEiRPJ29ubfv75Z6PH55Fnn55ISUkhjuNoz549phBk\nVmi1Wpo6dSopFAoqLy/nJf7bb79N3bp1o/Pnzxs9Ps907NBffHw8SaVSOn78ON+CzIq4uDiys7Oj\nzMzMp36vvLycKisrOxW7qamJZs2aRVKplA4cOGCITFZ0zDxarZZmzpxJDg4OdOLECb5FMUen09Gy\nZctIKBS2m0WsuLiYdu3aRW+//Tb5+PgQADp69GiH4//nP/+h8PBwcnR07NTPmRkdP6uu0Who+vTp\nJJVKKSUlhU9RTFGr1TRjxgySSCRtHeRr167RF198QbNnzyZPT08CQAKB4KG8PXl5ec+MfffuXVq5\nciVJpVIaMmQIFRQU8P06fNK5LBlarZaWLl1KHMdRQkKC1SUuKiwspOHDh5O9vT3Fx8fTrFmzSKFQ\nEAASCoUkEomemNippqbmiXFv3rxJCQkJ1KNHD+rRowclJycbPYECA/TLSbht2zbq3r07DRs2jG7c\nuGFsUSZHq9XSwoULSSwWk0AgIAAkEok6nBdaIpE8NBJramqiEydO0O9//3saMmQICQQC6t27N33w\nwQd09+5dhm9qVPRPaJmfn08DBw4kmUxGiYmJvE2e8c2tW7doypQpBIAiIiJILpcTx3HEcVyHk4m7\nuLjQ4sWLacqUKQ+lofP396df//rX9N1331lDTfMohmVD1Wg0lJSURPb29uTr60u7du2ymMXUyspK\nSkhIIJlMRkFBQW2nUpuammjHjh3k4+NDHMd1qPbx8PCgESNG0Jtvvknvv/8+7dq166E801aKcVLp\nlpSUUGxsLIlEIgoICKAvv/zSbGeli4uLacWKFWRvb09ubm6UnJzcbt9Nq9XSgQMH2lLCPam/w3Ec\nTZs2jcGbMMe4SbwLCgpozpw5JJFISC6XU3x8fIdGIXyj0Wjo4MGDFB0dTUKhkBQKBX3yySdUX1/f\noZ8/ffo0vfzyy8RxHInF4sf6O+ZwHJoBxs8AT3R/0mzt2rXUr18/AkCBgYG0cuVKys7ONlmzVldX\nR/v376fY2Ni2fszYsWMpNTWV1Gq1XjH/9a9/0bRp00ggELSZSCKR0Nq1a42s3iLgxzytaLVayszM\npPj4ePL29iYA5OTkRJMnT6aPPvqI0tPTqbS01OBympubKT8/n1JTU2np0qUUERFBIpGIBAIBjRgx\ngtavX08//fST4S/0XwoKCtpGZwBo+/btRottQZjm1ptW8vLykJGRgczMTGRlZaGsrAzA/bNK/v7+\nUCgU8PLygpubG5ycnCCVStvu2Wq9d0ulUqGurg4lJSWoqKhAcXExrl+/Do1GA5FIhAEDBmDUqFGI\niopCVFQU3N3deXufiooK/PWvf8XkyZMxYsQI3soxU/aa1DyPolQqcfnyZVy5cgU3btxAeXk5SktL\nUVFRgbq6OqjVatTX10Oj0cDe3h5isRgODg5wdHRE7969oVAo4OnpicDAQAQHByMoKAhSqZTV63Q1\n2JrHhkWzl+kdozYsG5t5bOiNzTw29EYEYC9rETYskvP/HyIL2lLTHX08AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "graph_matrix(A,show_weights=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Transition matrices\n",
    "\n",
    "### Using graphs to represent the transition from one state to the other"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After collecting data about the weather for many years, you observed that the chance of a rainy day occurring after one rainy day is 50% and that the chance of a rainy day after one sunny day is 10%. \n",
    "\n",
    "The graph that represents the transition from the weather on day 1 to the weather on day 2 can be expressed as an adjacency matrix, where the edge weights are the probabilities of weather conditions. We call that the transition matrix.\n",
    "\n",
    "Write the transition matrix ${\\bf A}$ for the weather observation described above, and use the helper function to plot the graph. Your graph will look better if you assign the labels for your nodes, for example, use the label `['sunny','rainy']`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Properties of a transition (Markov) matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- $A_{ij}$ entry of a transition matrix has the probability of transitioning from state $j$ to state $i$\n",
    "\n",
    "- Since the entries are probabilities, they are always non-negative real numbers, and the columns should add up to one.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The weather today is sunny. What is the probability of sunny day tomorrow?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The answer to this question is quite trivial, and in this example, we can directly get that information from the graph.\n",
    "But how could we obtain the same as a matrix-vector multiplication?\n",
    "\n",
    "Write the numpy array ${\\bf x0}$ representing the current state, which is the weather today. Your vector should be consistent with the transition matrix. If the first column of ${\\bf A}$ corresponded to transitioning **from** a sunny day, the first entry of the state vector should be sunny."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can now obtain the probabilities for tomorrow weather by performing a matrix-vector multiplication:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The weather today (Thursday) is sunny. What is the probability of rain on Sunday? Should I keep my plans for a barbecue?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What is the probability of sunny days in the long run? \n",
    "\n",
    "Let's run power iteration for the equivalent of 20 days. We will store all the state vectors ${\\bf x}$ as columns of a matrix, so that we can plot the results later.\n",
    "\n",
    "We will store all iterations in `allx` in order to view the running probabilities later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "its = 20\n",
    "allx = np.zeros((2,its))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will start with the initial state of rainy day:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x = np.array([0,1])\n",
    "allx[:,0] = x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What are the probabilities after 20 days?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print('Probabilities of initial state:', allx[:,0])\n",
    "print('Probabilities after 20 days:', allx[:,-1])\n",
    "plt.plot(allx.T)\n",
    "plt.xlabel('')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What if we were to start with a random current state? Would the weather probabilities change in the long run? Write a similar code snippet, but now starting at a random initial state. Remember that the sum of the probabilities has to be equal to 1 (or that columns and state vectors should sum to 1). Think about normalizing. Which norm would satisfy this property?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print('Probabilities of initial state:', allx[:,0])\n",
    "print('Probabilities after 20 days:', allx[:,-1])\n",
    "plt.plot(allx.T)\n",
    "plt.xlabel('')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "You can run your code snippet above a few times. What do you notice?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Summarizing:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- You started with an initial weather state ${\\bf x_0}$.\n",
    "\n",
    "- Using the transition matrix, you can get the weather probability for the next day: ${\\bf A  x_0} = {\\bf x_1}$.\n",
    "\n",
    "- In general, you can write ${\\bf A  x_n} = {\\bf x_{n+1}}$.\n",
    "\n",
    "- Predictions for the weather on more distant days are increasingly inaccurate.\n",
    "\n",
    "- Power iteration will converge to a **steady-state** vector, that gives the weather probabilities in the long-run.\n",
    "\n",
    "$${\\bf A  x^*} = {\\bf x^*}$$\n",
    "\n",
    "- ${\\bf x^*}$ is the long-run equilibrium state, or the eigenvector corresponding to the eigenvalue $\\lambda = 1$.\n",
    "\n",
    "- This steady-state condition ${\\bf x^*}$ does not depend on the initial state."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example 2: What are students doing during CS 357 lectures?\n",
    "\n",
    "Consider the following graph of states. Assume this is a model of the behavior of a student after 10 minutes of a lecture. For example, if a student at the beginning of a lecture is participating in the class activity, there is a probability of 20% that they will be surfing the web after 10 minutes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/StudentLecture.png\" alt=\"Sparse\" style=\"width: 500px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write the transition matrix ${\\bf A}$ corresponding to the graph above. Build your matrix so that the columns are given in the following order:\n",
    "```python\n",
    "activity_names = ['lecture', 'web', 'hw', 'text']\n",
    "```\n",
    "Plot the graph using the `graph_matrix` helper function, to make sure you are indeed building the correct matrix. You can use the list above to label the nodes. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Build the initial state vector ${\\bf x}$. Recall that this array should follow the same order as the matrix, here defined in `activity_names`.  In the initial state (suppose that we are at the 0 minute mark of the lecture, i.e. the beginning of the lecture), students have an 80% probability of participating in lecture, 10% probability of texting friends and 10% probability of surfing the web. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is the probability that students will be working on their HW after 10 minutes?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is the probability that students will be sending text messages after 30 minutes?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Would your answer above change if you were to start from a different initial state? Try for example x0 = [0, 0.2, 0.1, 0.7]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Steady-state will not depend on the initial state, but before achieving steady-state, the probabilities will be different!\n",
    "\n",
    "Get the steady-state for this problem:\n",
    "- Start from a random initial state. \n",
    "- Make sure the probability of the initial state sums to 1. \n",
    "- Assume 20 iterations of power method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "its = 20\n",
    "allx = np.zeros((4,its))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.plot(allx.T)\n",
    "plt.legend(activity_names)\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Do you think you achieved steady-state at iteration 40?\n",
    "\n",
    "Note that iteration 40 corresponds to the 400th minute of a lecture, which in this case does not exist (luckily? :-))\n",
    "\n",
    "Does the class achieve steady-state at the end of lecture (about 70 minutes or 7 steps)?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use `numpy.linalg.eig` to see how the steady-state (solution of power iteration) is the eigenvector of ${\\bf A}$ corresponding the the eigenvalue 1. The eigenvectors are given in the columns of ${\\bf U}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "lambdas, U = la.eig(A)\n",
    "print(lambdas)\n",
    "print(U[:,0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not what you expected? Remember that if ${\\bf u}$ is an eigenvector of ${\\bf A}$, then $\\alpha {\\bf u}$ is also en eigenvector. Can you show that the eigenvector from the power iteration solution is the same as the one obtained using `numpy.linalg.eig`?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example 3: Page Rank"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose we have 6 different websites that we want to rank. They appear in the dictionary `name_mapping`, that indicates the website string name for each node in the graph (or position in the transition matrix)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "name_mapping={0:'Google', 1:'Facebook', 2:'Amazon', 3:'Wikipedia', 4:'Twitter', 5:'YouTube'}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The list `edges` contains the outgoing link information. Each list entry `[x,y]` indicates that a node x (or website x) has a link to node y (or website y). The code snippet below construct the adjacency matrix, and plots the corresponding graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "edges = [\n",
    "    [0,1], [0,5],\n",
    "    [1,2], [1,3],\n",
    "    [2,3], [2,4], [2,5],\n",
    "    [3,0],\n",
    "    [4,5],\n",
    "    [5,0]]\n",
    "\n",
    "A = np.zeros((6,6))\n",
    "for i, j in edges:\n",
    "    A[j, i] = 1\n",
    "\n",
    "graph_matrix(A, list(name_mapping.values()), show_weights=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Definition of the **PageRank model**:\n",
    "\n",
    " “PageRank can be thought of as a model of user behavior. We assume there is a random surfer who is given a web page at random and keeps clicking on links, never hitting “back”…”\n",
    " \n",
    " http://infolab.stanford.edu/~backrub/google.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    " Let's put numbers to the statement above and see how it applies to the example above. Suppose a random surfer is initially at website 0. His initial state would be:\n",
    "\n",
    "[1,0,0,0,0,0]\n",
    "\n",
    "The web surfer will click on links on current websites, and he will keep doing this until achieving steady-state (he will never hit \"back button\"). So for the example above, if he is on the Google website, he will either be going to Facebook or YouTube next. Since there are two links, we will assume that each website has the same probability of being the next one, so there is a probability of 50% he will go to Facebook and 50% he will go to YouTube.\n",
    "\n",
    "We now need to modify the adjacency matrix above in order to get a matrix that satisfy the properties of the Markov matrix, specifically the one that the column sum is equal to 1 (here representing the probabilities of the outgoing links).\n",
    "\n",
    "Here is what your adjacency matrix looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How should you change the matrix `A`, to model the behavior described above? Make the Markov matrix M."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "graph_matrix(M, list(name_mapping.values()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the steady-state for this problem:\n",
    "- Start from a random initial state. \n",
    "- Make sure the probability of the initial state sums to 1. \n",
    "- Assume 20 iterations of power method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.plot(allx.T)\n",
    "plt.legend(list(name_mapping.values()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At steady-state, what is the most likely page the user will end up at, when starting from a random page?\n",
    "\n",
    "Creating the list of strings `ranking` with the rank of the websites (the string names), starting from the mostly likely the user will end up at, to the least likely. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What happens when the web surfer ends up in a website without any outgoing link? \n",
    "\n",
    "Let's repeat the example above, with a simple change to the outgoing links. Page 5 (YouTube) no longer has an outgoing link."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "edges = [\n",
    "    [0,1], [0,5],\n",
    "    [1,2], [1,3],\n",
    "    [2,3], [2,4], [2,5],\n",
    "    [3,0],\n",
    "    [4,5]]\n",
    "\n",
    "A = np.zeros((6,6))\n",
    "for i, j in edges:\n",
    "    A[j, i] = 1\n",
    "\n",
    "graph_matrix(A, list(name_mapping.values()), show_weights=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that you can no longer use the method above to obtain the Markov matrix, since the sum of the fifth column is now zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Once the web surfer reaches a page without outgoing links, we can assume that he will not stay on that webpage forever!\n",
    "\n",
    "We then assume that the web surfer will move to any of the other webpages with equal probability. Construct the Markov matrix using this assumption."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the steady-state for this problem:\n",
    "- Start from a random initial state. \n",
    "- Make sure the probability of the initial state sums to 1. \n",
    "- Assume 20 iterations of power method\n",
    "\n",
    "Did your ranking changed?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### One remaining issue!\n",
    "\n",
    "The Markov matrix does not guarantee a unique solution! Let's take a look at this example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "B = np.array([[0,0,1,0,0],[1,0,0,0,0],[0,1,0,0,0],[0,0,0,0,1],[0,0,0,1,0]])\n",
    "graph_matrix(B,  show_weights=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The matrix `B` has two eigenvectors corresponding to the same eigenvalue of 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "l,v = la.eig(B)\n",
    "print('The eigenvalues are:')\n",
    "print(l)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some of the eigenvalues are complex, but we are interested in the two eigenvectors corresponding to the eigenvalue 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "v[:,2]/la.norm(v[:,2],1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "v[:,3]/la.norm(v[:,3],1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hence both\n",
    "$$\n",
    "{\\bf x^*} = \\begin{bmatrix}\n",
    "    0 \\\\\n",
    "    0 \\\\\n",
    "    0 \\\\\n",
    "    0.5 \\\\\n",
    "    0.5\n",
    "\\end{bmatrix}\n",
    "\\qquad \\textrm{and} \\qquad \n",
    "{\\bf x^*} = \\begin{bmatrix}\n",
    "    0.33 \\\\\n",
    "    0.33 \\\\\n",
    "    0.33 \\\\\n",
    "    0 \\\\\n",
    "    0\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "are possible solutions!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Perron-Frobenius theorem (CIRCA 1910):\n",
    "If ${\\bf M}$ is a Markov matrix with ***all positive entries***, then ${\\bf M}$  has unique steady-state vector ${\\bf x}^*$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Definition of the **PageRank model** (complete):\n",
    "\n",
    "“PageRank can be thought of as a model of user behavior. We assume there is a random surfer who is given a web page at random and keeps clicking on links, never hitting 'back', but **eventually gets bored and starts on another random page**.\"\n",
    " \n",
    "To model the behavior of a web surfer getting bored, the proposed **Google matrix** is defined as:\n",
    " \n",
    " \n",
    " $$ G = \\alpha \\, {\\bf M} + (1-\\alpha) \\,\\frac{1}{n} \\mathbf{1} $$\n",
    " \n",
    " \n",
    "where $\\mathbf{1}$ is a matrix of ones with same dimension $n$ as the Markov matrix ${\\bf M}$. We divide $\\mathbf{1}$ by $n$ to enforce that the columns sum to 1. In this model, a  web surfer clicks on a link on the current page with probability $\\alpha$ and opens a random page with probability $(1-\\alpha)$, where $0< \\alpha < 1$. A typical value for $\\alpha$ is $0.85$. \n",
    "Note that the Google matrix ${\\bf G}$ has all entries greater than zero, and guarantees a unique solution.\n",
    "\n",
    "\n",
    "\n",
    "Construct the Google matrix ${\\bf G}$ for the Markov matrix above, using $\\alpha = 0.85$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "graph_matrix(G, list(name_mapping.values()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the steady-state for this problem:\n",
    "- Start from a random initial state. \n",
    "- Make sure the probability of the initial state sums to 1. \n",
    "- Assume 20 iterations of power method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Overall complexity of PageRank problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In general, for the dense Google matrix with shape $(n,n)$, the cost of each matrix-vector multiplication during the power iteration method would be $O(n^2)$. \n",
    "\n",
    "$$ {\\bf x}_{i+1} = {\\bf G \\, x}_{i} $$\n",
    "\n",
    "\n",
    "In a PageRank problem, $n$ is very large! But let's look at each part of the Google matrix:\n",
    "\n",
    " $$ {\\bf G} = \\alpha \\, {\\bf M} + (1-\\alpha) \\,\\frac{1}{n} \\mathbf{1} $$\n",
    " \n",
    "The Markov Matrix for the PageRank problem will likely be very sparse, since webpages will only have a few number of outgoing links. Hence, the calculation \n",
    "\n",
    "$$ {\\bf x}_{i+1} = {\\bf M \\, x}_{i} $$\n",
    "\n",
    "will take advantage of sparse operations, which can be in general approximated with $O(n)$ complexity.\n",
    "\n",
    "We need to look at the complexity of performing\n",
    "\n",
    "$$ {\\bf x}_{i+1} = \\mathbf{1} \\, {\\bf x}_{i} $$\n",
    "\n",
    "Since $\\mathbf{1}$ is a dense matrix, some would quickly jump to the conclusion that the complexity is $O(n^2)$. But it is not! Can you think of how this operation can be done efficiently?\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}