{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Community detection algorithms" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On this seminar your are asked to implement simple community detection algorightm. It is called [Markov Cluster Algorithm](http://micans.org/mcl/) (MCL).\n", "\n", "It consists of the following steps:\n", "\n", "**Input:** Transition matrix $T = D^{-1}A$\n", "\n", "**Output:** Adjacency matrix $M^*$\n", "\n", "\n", "1. Set $M = T$\n", "2. **repeat:**\n", " 3. *Expansion Step:* $M = M^p$ (usually $p=2$)\n", " 4. *Inflation Step:* Raise every entry of $M$ to the power $\\alpha$ (usualy $\\alpha=2$)\n", " 5. *Renormalize:* Normalize each row by its sum\n", " 6. *Prunning:* Remove entries that are close to $0$\n", "7. **until** $M$ converges\n", "8. $M^* = M$\n", "\n", "\n", "\n", "As a result you should get a cluster matrix s.t. elements of the cluster correspont to nonzero elements of the rows of the matrix. \n", "\n", "* What advantages and disadvantages does this method have?\n", "* Run this method for network [1](https://www.dropbox.com/s/so0ly2ozh2pzxp6/network1.mat?dl=0), [2](https://www.dropbox.com/s/xcswyhoeehq95v2/network2.mat?dl=0) and [3](https://www.dropbox.com/s/cwshsfr2d8fn470/network3.mat?dl=0).\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "plt.xkcd()\n", "import networkx as nx\n", "import scipy\n", "import scipy.io\n", "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQEAAAD8CAYAAAB3lxGOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX2cHEWZ+L+9O/syb7uz79mEFwOEEALcQUggMRAQEFGD\nKBEQCIi8yPGmoj8O0JOggIAeyCkHUe/gVEAQCASUQ1FA3hLC+1vIhZcg5mWzr7M7Lzu7s/P8/ujp\nTm9vT2/P7kzv7Ka/n898Zqaruqq6up6qp6qeqgIPDw8PDw8PDw8PDw8PDw8PDw8Pj2IgIjUicpWI\nvCAiL4nIzSIy3cX4fSKyXETuFZE/iMgNIjJTRCpE5AIReUpEXhWRO0RkrgvpCYjID0XkRRH5o4gs\nyKbnNhF5WUSeFZHLRCTgQloOEZHVIvKMiJyevbariPw8m5bnReQKEQkVOR2zRORaEfm9iNwqIgdl\nr+8jIv+dfT9Pi8hFIlJVhPgVETldRJ4UkT+JyLHZfPhZtsw+LyJXikhIRBqzZejF7Oc6EWkoYFrK\nReQbIvI3EXlUROaLyAwR+UG2/P5eRM4VkSoR2T2bXy+LyHMicrmIBAuVloIgIvuKyDax5nAX4p8m\nImtyxJ+LM4uYnjIRecxhOjZKEStLETnGIs5bcqRlk4jsVqR0fCVHnGtzXH9FROoKGH+lqBXNeOgV\nkX8qQFpqRG0YxsMHIrJrIfJm3Ihao70tIvL888/LMcccI4sWLZKHH35YS2ybiFQXMX5FVO1DNm3a\nJOedd558/vOfl9/97nd6bm3cuFG+/OUvy7x58+TWW2+VTCajOe1epDRdLCLS1tYmS5YskRdeeEFP\ny29/+1tZsGCBfP7zn5fXX39du/xQkdJxoBbBbbfdJj//+c+HlaJ7771XDj30UDnuuOPk5Zdf1i7/\nbxHSsYeISCaTkTvuuENOPvlkueGGG/T3kMlk5LbbbpN58+bJiSeeKP/3f/+npeVXBUzDv4uIxONx\nOeuss+SJJ57Q8+G+++6TQw89VD7zmc/ISy+9pF//05/+JIcffrgcccQR8uSTT2qX3xKRsnGm5QER\nkfb2djnhhBNk/fr1epyrVq2S448/Xs466yx577339Ov33HOPHHLIIfLZz35WXn31Ve3yHwqVP+NC\nRE4QEXnvvfdk2rRp8sgjj8gLL7wge+yxh7zyyitaYk8vYvzTtQzdfffd5ec//7k89thjMmPGDHnp\npZckkUjI7rvvLitXrpTXX39dFi9eLL/97W+1dF1fhPTM1gr2F77wBWlsbJQHHnhAL2z777+/rF27\nVu6//37ZbbfdJB6Pa2nZowhpeUZE5I477pADDjhAtm/frheq1atXyz777CMvvPCCPPTQQ7LLLrtI\nb2+v5jynwOn4nohaAR533HHywgsvyPLly+Xmm28WEZH//M//lEWLFsnLL78sd955p+y9996STqe1\ntBREBReR90VEFi1aJOFwWB566CEREXn00Udl9uzZ8vzzz8vDDz8su+yyi0SjUVm3bp3ssssu8pe/\n/EX++te/yi677CKbNm3S0nTMONJRJiLpVCole+65pwSDQXnttdfUGu9Xv5KDDz5YHn/8cbnmmmvk\nwAPVOvzBBx+UfffdV9asWSMPPvig7LrrrtLX16elZe9C5M+4EJHbRUQuu+wy+fGPf6wXsptvvlm+\n9a1vaX//s4jxV4lIYmhoSDo6OkQTwAMOOECeeuop+fWvfy1nnnmmnq4nnnhCjj76aO3vMwVOi0+y\n3ZI777xTjjnmGFm+fLleCSxatGiYVrBs2TK5//77tb+nFTgtC0VEurq6ZM6cOdLe3i5GPvWpTxlb\nNzn99NONlePZBU7LrSIiF154oaxcuVJERFauXCnnnHOOWor33ls2btyop+WII46Qp556Svt7bIHS\n8JKIyMcffyxf+tKXZNWqVSIicvTRR8tf/vIXPe4zzjhDLzO//vWv9euXX365XH/99drfq8aRjnJR\ntWP5+OOP5aCDDtJb9mg0KqlUSkRE1q5dK3vttZeIiCxZskT+9re/6Wk55ZRTjJquo27tuFQXB3wC\n4K233mLhwoX6xT333JOPPvpI+9tarMgVRUkBl5WVldHQ0EB/fz/nn38+gUCAT37yk7z11lsceuih\nuv899tjDmK5pBU7OdcAhf//737nyyiu5/fbbURRFd1y/fj3z5893Ky1nA9x2221kMhkWLVrE4Ycf\nztq1awEs8+Xvf/97sdLSD3D66afz3e9+ly9+8YtcffXVXHjhhaRSKbq7u9lrr72KnZaLAHbZZRe6\nu7upr68HcudDscqzoihDwHlWaampqaGyspK3336b5cuX853vfAeAt99+m0MOOWRcaSl2JZABUBSF\ndDqtX0ylUtTU1Gh/e4uchn0BNm7cyIIFCwD485//jM/nG5GugYEBamtrC54uEakBvplOpznllFPY\nf//9efbZZ/nggw946aWXiEajZDIZMplM0dOSZTbAn/70J5YvX87jjz/O2Wefzamnngq4/r6aADZt\n2sS0adNYunQpu+++O3fccQeKojA4ODjMczHSoijKGmAjQGdnJ9OmTdOuW+ZDMfNHUZSHtd+dnZ20\ntLTobr/61a847rjjuPbaa/n617+uXx8aGhpXWopdCWwFmDlzJuvXr9cvvv7668ycOXOYn2IgIpXA\n+UNDQxx//PFcccUVrFy5klBIne1yMV3VQEV/fz/HHHMMc+fOZd26dbS1tbFmzRref/99Zs6cyYYN\nG9xIC0AlqAXm+OOPZ+bMmZx55pl0d3czMDAwIl/eeOONYqblIICf/OQn3H333Xzta19j9erV3H77\n7fh8PsLhMFu37oiyGGkRkX2BWQDJZJJgUJ1hmzlzJu++++6IuItZbkRkifZ7aGiIqip1NvRvf/sb\n119/PWvXrmXZsmW6f3Ma3ZItx4jIGSIiTz/9tOyzzz6yadMmeeWVV2TatGmyYcMGrd/ymSLGv7eI\nyPvvvy+77767PP3007J69Wr54x//KP39/bJlyxZpbW2VdevWyccffywHHHCA3h8Ukf9X4LTonWpj\nH1MbE7juuutk2bJl0tXVJffee6/ssssuMjg4qHkt2Bx0Ni33iYicf/75el/2lVdekTlz5oiIyE03\n3SRLly6Vzs5OWbVqlbS2tur9UREpaPdN1Pl/OfTQQ+Wxxx7T09LS0iKZTEYuvfRSOeeccyQajcrt\nt98uc+fO1WYO+kXEX6A0PCOiztjMnTtX/vCHP0gmk5Gf/vSn8rnPfU46Ozvl4YcflmnTpkl/f7/c\nd999csghh8i2bdvk6aeflsbGRuO4yoHjSEeZiHwkIvKPf/xDamtr5e233xYRkSuuuEIuvvhi+dOf\n/iQPPvigPmB44403ygknnCBdXV1y//33y/Tp02VgYEBLS3Mh8mdciEitiCRERH75y1/KbrvtJrNm\nzZK77rpLS+RHIlJexPg/IaLODhx11FHy+c9/Xr7yla/Il770JX0w5X//939l3333lRkzZsh1111n\nnCIsaN83+4K/JKpBjIiIXHvttXo6+vv75dvf/rY0NTXJwoULZc0a3bThnkKmI5uWU0VEPvzwQznw\nwAPlwAMPlNmzZ8tzzz0nIiKpVEouv/xyaW5ulgULFsizzz6rpWVVEdJyk4g6hbzffvvJrFmzZL/9\n9pPHH39cRER6e3vlvPPOk8bGRjnyyCPljTfe0NLyswKmoVtE5Mtf/rIsXrxYlixZIu3t7TIwMCBX\nXHGFNDc3y/z58+WZZ54REXVw+cYbb5TW1lbZb7/95JFHHtHStG6c6ajVysLRRx8tixcvlk9/+tOS\nTqflnnvukSOOOEJOPPFEOeOMM+TEE0+UTCYj/f39ctlll0lzc7Mccsgh8vzzz2tp+X2h8iefB9hN\nVIuq20XkX7UWQ0SOF2uGRORTLqTrbxZxb7K4ZuTcIqfpc6PEb0znLkWI3ydZ+4mhoSFpb283Vn65\n+IeIfKIIaWnOhi0i4iQdIqrtSX0B05DLQCofeiVr5TjOtDw8WkRZnrZx+7vkYdiljO7FUcIvBv7D\nwukMRVF+I+rc8g8BQe17vQZcqyjKxkLEP0raGoErgIOBGBAF3gf6gDeALwILgPeAKuC67EBRMdOk\noI7QX5BNx/eBAPCvwAbgU8D9wI2KonQXKQ0h4AfAt4DXgX8C7ga+A8wFrkQdLDsaeAC4QVGUziKl\npQW1fJwMvAn8M3AfsAL13VwKvAUclk3jzYqixAoYvx817w8HKoAXUcdxnkN9T++h5sODwA2os14r\ngC3AMcD/opbnfxQgLQ2oeT8PKAfWoD53N2oZ6QBSwMPAAcA64BLUd3UUO95V13jTkk+iTxBRa/B7\n7rlH/uVf/kUefPBBY41+sGuJ8cgbUbUCvxSxW5YP2QrSY7IgqlnuGyIi3/jGN+SII46QX/ziF7Jo\n0SK56aabtErg/olOp4eHR27GVeuKyD7A+m3btjF//nw2bNhAIBDg448/ZtGiRXz88ccACSCkKIoU\nIsFuIqodeAh1Lns6EAEas58A4AfC2f81qN2Jyuy39juY/fhRVc0yRuZ7GhjKflKoXYQ+II6af1HU\nrkwf6txvHOhBnQLqyLr1AG1ATFGUDJMAEfGh5l8AqANqTd8hoAE1bwPZT3X2em32O5i9XgX4sM5f\njTQwmP2kUQ2VYqh5rOVvj+E7hZrfbezI9wTQCfxDUZR4IfKh2IiIYid/vnGGPwvg1Vdf5bDDDiMQ\nUFe+7rrrriQSCaLRKLW1tQFU4ekWkdEqgnPZ8UK0Tz9q5vdmf2svUAwfUF++9inPflehCl8laoEJ\nZ/9rv8OoAl6X/dQb/regFjQAMpkMiUSCzs5O2tvbicfjJJNJYrEYnZ2dRKNR+vv7GRgYIJVK0d/f\nz+DgIIlEgr6+PpLJJOl0eoRRkKIo+Hw+X3l5ua+8vJyqqqpAOByuq6mpwe/3EwqFqK2tJRQKEQ6H\nCYfDBINB6uvrmTZtGvX19ZSVDTP3yIhIG2p/tQu1kmjP/u5ErVC6st9R1IKezOb1QDZv02QNvVAF\nqgy1rJQbvisYXuGF2CGM/uw7bwFmAM3ZvGzIXg+zQ4CHEY/H6enpIRqN0tvbS1dXF9FolEQiQSKR\noL+/n1gsRjQaJR6PE4vFSCQSpFIpPX9FBHNR0/K5oqLCV1FR4ff5fPj9fgKBAMFgkFAoRE1NDbW1\ntfp3VVUVkUiE5uZmampqCAQC+Hw7REZE+rN52J399GQ/8Wwea//bs/562FGmtXzXyvRQNq/L2dF4\naBWflr/12W/tU5PNywhwljkvTeRs8MdbCVSD+uIqKyuHOZSVlRmvDVx11VXKiSeeSEtLC9XV1bS0\ntNDa2kp9fT2BQIDa2lpqa2t/qb2QYDBoLtwFIZ1OE4/H6e3tJRaL0d7eTldXFz09PXR1ddHe3k53\ndzfbt2+ns7NTF/ItW7aMsF6zQlEUqqqqqK6upqKiQi9gfr8fn89HWVmZ/hERhoaGSKVSDA0NkU6n\nSaVS9PX10dfXpxfu0fD7/TQ0NNDU1EQkEilrzaJda2xspL6+noaGBmpra6mrqyMQCAwzWy4EmUyG\ngYEBkskk0WiUtrY2tmzZQltbm56/3d3dxGIxYrEYvb29ujBrgu8kjwHC4TChUIhQKEQgEKCyslLP\nX0VR9G9Az2etYh4cHCSdTpNMJkkkEnpl4oTKykoaGxuZMWMGjY2N1bW1tdX19fUtkUiESCSiV9h1\ndXVEIhHq6upobGyktraWioqKMeetiOhlNh6P62Wkr6+PaDRKNBqlr6+PeDxOIpHQ8zaRSJBMJgkE\nAtfE4/GfKIrSYw57vJXAdoDp06cb7bnp7u5GRPD7/aDWdJufe+652v7+fh544AFqa2vp7R3dotHv\n91NdXa23gppgaSa/2osWEb2FzWQyDA0NkclkSKVSJJNJBgYG9IwbrZD5fD7q6+tpbGyksbGR1tZW\n5s6dy/Tp02loaKCuro6WlhaCwaDeUtfX1xOJRKiqqhrWUhSCdDqtv3ztpcfjcTo6Oti2bZsuVB0d\nHXR2dtLV1cXzzz9PW1sbyWTS9jm1PNWeQxMkTZgAPW/T6bReUQ0NDTE4OKhrPalUSteM7KisrKSh\noUEX3pqaGnbffXcCgQBGIaqvr6empoaamhoikQj19fV6q11dXU0gEKC8vLDjmCKiaxiaAGlrF9rb\n2/VKOR6Ps337djZv3kx7ezvvvfee3ogYzXet0BqEqqoqqqqq9PJdUVFBeXk5IkI6ndYrUk2gtfwd\nLXxQGyG/369rM+vWreOLX/wiBx100HeB74pIhaIoaeM94y2xbwDMnz+f9evXs379eubMmcPPf/5z\nli5dqvlZCxxcVVVFZ6c6wxSNRonFYmzbto2uri6969Db26vXylqNp6ncsViMZDKp1+JGlS+r6ukt\nbHl5OWVlZXpGV1ZWEgwGCYfDuvpXU1NDOBymoaGBhoYGIpEIDQ0NhMPhYRpIPB6nq6uLXXdV92i4\n66672LhxI7FYjO9///uAuvilu7ubP/wh/yXc5eXl+sdKiLRKSVtIkg/JZJLt27fT1dVFR0cH3d3d\nRKNRuru76e7upq+vj1QqRSKRIBaLMTAwQDqd1vPXmLea4Pl8PsrLy6moqKCyspLq6mo9f7WKxO/3\nU1tbS1NTEzNmzKClpYW6ujr8fr9j7WPlypV89NFH9Pb2cu211+b97PmiKAqnnXYa4XCYu+66a5jb\n3XffTWtrK4FAgKOOOsryfq0S0SpsrYvY1dVFZ2enrunEYjH6+/v1BiqVSjEwMKDnd2VlJZFIRC+n\nfr+fqqoqvbHRuoNa1zAcDmtatK5xmvO4urqari59xtCPOuax49nHm3ki8gjw+dWrV3PhhRdSU1OD\nz+dj1apV7LHHHgBfA2446aSTmt544w3effddjjvuOHw+H62trdTW1hIMBolEIoTDYc4+e8dK1XXr\n1lFdXa0/bEPD2KxnN23apGsC2qqwu+66i76+Pr2mj0aj/Nd//deIe40qpfG/1bXRhzxGYhXeRJFM\nJvH5fONSWwvFRORLrvdYSu/o7rvv1tX9Sy+91NE9J510EprsAa2KomwraKJE3QOuX0Skr69PNm7c\nKENDQ9r04FpR9/F7/5xzzpFp06aJiMjs2bNtjJ1EH/BzgpVf8zWr3wwfWMz5MYdp5TYeCh2el5ad\nOy1WbkbZE3VGbxjjHnlTFOVdYD/gsVAoxF577fWPrDr9U+BoRVEGgZ6amhp9HEBbxTeGuPSPhpi6\nBbnUTTs3Yxjm67n8WrmVIoUe/PMoPIV8R1Zl0yh7qLMKwyjIKJaiKO8BnxXV5FEB4oqiGDu4iUAg\nQDKZRET0pZoeHh7Fxyh7iqKMEL6CDmXb2JbHg8EgIkIymbTVBLTRfu03DG+RzW5GzDVgvn11q/tz\n3esk/rGkodBkX/yEpsHDnmK/F6PsBQKBEcJX7E1FNGLhcFj9EYt5moCHh4sYZQ8LA63CTmrnplfb\nKisajbJ8+fIRHoytVa4Rebt++2j9KrP7VVep+0FeffXV424p7e4bLUwtHVoai9EqOA3TmJaJplTT\nUqx3NJa0OMUoe83NzQU7ryEvROSm1atXCyDr1q0bdVQTrEfm7UZBze65wswVdq777eLLdd9YKWRY\nHsVhMr4jk+yN2DHLLU2gT1NJ+vpUOwVzyyyGFl3y6Ieb73d6zW6+30n4udKYD2YNpBBhehSXUnpH\nTrUSk+yFze6udQeMCbnpppvGNIhndNP8W3Uj7MhV+Rh/O+0ejKUiyxW3h0e+OC0/pkpg4gcG+/r6\njFsie3h4FBlTJTBC+NzSBGLajEA8HrecIrRTyQuh1ufSPKxabasWPZ/pSq+FnxzsLO/NKHuoS5OH\n4ZomoAl+X1+f8VANDw+PImOUPQx7ZGi4pQl01dTUUF5eTldXF5/73Od0B6vaOJ8pu3wH9ezCdhJv\nvmMQbuMZBjlnZ8kjo+yhbkwyDNc0AUVRCAaDjjdv8PDwKAwm2RvRHXBLE0gBVFVV5dwpJ99aOR8N\nIt9xA6fGSU60ELcpZuvmaRmjU6p5ZJC9KrObW5VAAnYsZDBiJ0hWGZrLmtDptFw+3QFjmFaDh6Uk\n/G5QagW71Jhoa0I7DLI34ug2t7oDCVBHKbMjlB4eHi5ikL0JWzsQB3XPQG3XWG3aws4Sz8ogyOxm\n16LnCjfX/U7SUqo1/ViZ6Oea6PgLSSk/gyZ7TOAU4QCo+6cNDAyM2JnYw8OjuGiyR/ZYeiNuaQJD\noG6aOTQ0hM/ns10fkKsf7nQwzukUn5NwzO5TqeWCwj9Hvv3iqZKPpY4me6jnGgzDLU1AoLRG0T08\ndibsZM8tTaAcYGhoaMR+8fka74ylJc53t1in5selOEVYDMZivOVROAqhfRpkb0QgbmkCFQCDg4P6\neIDdxp9mjH5FZESl4GRzUe0+42alZtXVHLadW67NTp08y2jupVap5PN8HoWnEPlvkL0Bs5tblUA1\nQH9/P1VVI2wVPDw8ioxB9vrNbm51B/ygHm6RPZrMkVXfWLEa/LMbEHQyfTjeNQ654srXvdSYagOl\nUxWD7I045irvSkBEqoFTgIWop9c+DTyIao74dWAJar/jT8B/KYqSyrqRSqU8TcDDYwIwyN6I7kBe\nlYCIHAA8BMw0XL6ArEUgww0RvgB8U0SOIdsdSKVSVFdXs2HDBvbZZ/hBKE4Hnez6y3bTeXa7AJmv\n2RkQ5ZvmqY6XB8WlUKbImuxh0R1wPCYgIo3Ak8DMN998k4svvphvf/vbbNq0CbLnqL/00kucfPLJ\nnH766bz99tsAs4Bfk93SKBaLEQqFPNNhDw+X0WSPrPWukXwGBpcC9U899RTHH388+++/PzU1NRx1\n1FGICG+88QZLly7lhBNO4Nhjj+Uzn/mMdgrx4cAhg4ODxONxamtriUajOUc8jSP3Gnajo1ZuViP5\n5tkBo99c5FrUVIoj+B6TDyflqBBagFH2gKjZPZ/uQAfAfvvtx9q1a2lubqajo4Prr7+eTCbDzTff\nzHXXXcdXvvIVAF588UXuueceLrroIoBDolE17kgkwpFHHunI4s9Kvc/HzcnqQ6tKxKobYbeuoJh4\nA2/u49ZqQLfeqVH2gG6zez6awGPAe42NjTQ3N9PV1cWpp57K2WefTXl5Oa+99hpHHHGE7nn27Nm8\n//772t8ZHR0dAPrx4prweXPQHh7FxSR7XWb3fCoBIdufeOaZZ5g3bx6LFi3illtuASCdTg/zPDQ0\nZDxuTOnuViug+vp6brjhhhFGO4qisGLFCkuV3WjYY/xt/m/GzlgolwaQyy1X2MXuFniV5NgZ67uZ\nTPltlB0rmbriiiswyh4w4rzQfLoDnwT+6aOPPmLZsmU88sgjLFiwQHfcbbfdeP/995k5U504eOed\ndzjooIM056Genh4AfUzAw8Oj+EQiEYyyh8WYQD6awAEAjz32GJ/97GcJhUKsXbuWDz/8EIATTjiB\nm266iVQqxVtvvcWqVas44YQTtHv7tNqooaFB2/BwBMZzAcFZK2hlDmzuaowWZq4WPdfaAeNUopXm\nsjNSinkwmVr0sTKajNTX12OUPcbZHVAAWlpaeP755zn77LP53ve+x8knn0wikeCss85i1qxZNDU1\nceyxx/Lv//7vtLS0ALwBJLUNRkOhELfffrs3JuDhUUC0rrT5c+6552KUPaBvzJGIyEwRGZSRrDX+\nGRwclEwmY7y0WET+eOONNwogvb29snTpUoHRDwG1wu4+88fKbbzx2YU1WZhs6Z9MaS1FjLInIseb\nZdvxmICiKB+KyHHAlcCuqAMMCeBe4APgWeAcn8+3FTgO1ZT43xRFeUdEIj09PZSXlxMKhVi9enXO\nrcPEZlputGlADfM1Kz9OLQ/N8RmvT9bpu8mS3smav6WGUfawmCLMy2xYUZQngCcsnFZmv28VkbKs\n34zBva6np4dIJFJy/UYPj6mOSfZ6zO4FX0VoEn6Nuq6uLm2KIp+wAGdagtW1fDclsVt9mCuO0cL0\nGDtTIT+LWTacGjWZZG9cA4PjoTYajXpnEHp4TAAm2Su+JmBGREJAdUdHh6UmkKvPncst15iAlZbg\nRAOwqk3tzIatzI29Lo41noZU3Gd3GrZB9gbYseJXxw1NYAbA5s2bmTFjhgvReXh4GDHI3lZFUSZk\nj8EGUGuj5uZmYPiqPg3tv/Fj5Wa+pmFlZmxnNmx0y2UubGdCbGWc5DGcicqXUjRcmggWLVoEDJO9\nNit/bmwvFhgcHGRgYGDYqUNOVhFqOO0q5ONmhZNBQCdukwW3Vsu5jfH9T8Xnc8rzzz+PSfYsjwR3\nQxMIaoeQBgIjTkDy8PAoIibZGzEeAO5oAjW9vb0AhMNhYHznAOQyBBqt1s+leVhpJU61jMneykz2\n9I9GIZ4vn8HNUiwTJtmzNBl2QxNo2b59OwBNTU0uROfh4aFhkr0Ry4jBJU1AW8pYV1eX05MT7cDO\nWMipQZC5lbeLd7QpwnzGEEqRUmy53MLps+eTP6WYlybZs6wE3NAEqlKplPrD227cw8NVTLI3Yqdh\ncEcTiGibiNTW1nLhhRc6rjGtWmEn/XU7w55c91iFOZomUYo1fz6MNoYy2Z/PjtGebaoYOhlljxxj\nAm5UAmFtcKKmpoZbb73VtoCNtuLPyeo/K4vBXPEY02LVVSjlQlDMglrKz+0GpfT843nPRtnDwmQY\n3OkOREy1kYeHh0uUlCagKIq2ntnRCsBc1wrNWOMdT5cmVxh23R27gc98p1zN4VjhdKC1UOSrdU0V\ndd0J43lGk+xNmCYQisViBINBz5TTw8NlTLI3YcZCgXg8btx+3DGTubY3V3hOdyuy+m3GbrrSSZ7Z\nTXMa/TgZWykU+cYxGcuEEwpd5k2yN2Fmw6FEIuGZDHt4TAAm2RtxLDm4owlUG85Gz4vJWts73QtR\nw0n/26nJtJVGkG8atPvHagzl1hhCseOYCAr9PCbZm7DugC+dTuPz2Uc1FV6q1VSmneruRNCtws9n\nleRofuy6JmN9J6XYfdhZMcle2sqPG92BisHBQSoqKlyIysPDw4hJ9gat/LihCTiqBKZCzW6nntu1\n3sb785mC/v5sAAAgAElEQVT+G82wyhymVcueT/qKzVTQBieKXHnnpBIYkyYgIoeKyOMi8rqI3CIi\nYRE5Q0T+KiLvishqETk6691Rd8DDw6PwOOkO5C2ZInII8ILh0gHAJSZvs4GlInJj9h7Kytza2Hji\nWhQnuxwVwsgo34E+u5WTo6XXTSbqfU0FzcPODN8ge1bHAeSnCYhIALgH4Be/+AVLly7V3bZu3cqZ\nZ57Jvvvuy+WXX05/fz/AZZq7Zyjk4TExGCp/y5oi3+b5emDmG2+8wTXXXMObb74JqLXNF77wBWbP\nns3DDz/Mli1b+OEPfziedI8Lu1WDbsWvYbfppXGzUqu+fC7jIc3NfL9ZWxgtzFzXxvIMk5GJLiel\nguNKINvHv3hgYIAzzjiDH//4x7rb2rVrqays5Morr2TWrFn86Ec/4n/+53/M9ztO1GQtWFbCqV0z\nFjg7Nw0r4c4l8OZwrAb+rNJpF99o6bMLy6O0MFh+WgqVo0pARFqAuwGuvvpqDjzwQI455hjd/e23\n32bBggX6/+nTp9PR0UEmo3ZBFEXRf3t4eLiHSfYs5d2pJnA20PTEE09w4403Mnv2bH77298Si8V4\n+umnKS8vZ3Bwx+xDJpPB7/frAxLGKTOfz2fZ6iiKwooVK4Cp0bqM5RmsWvexhm2leVi19Lm0Co/C\nU0wNd8WKFSPebyaTMRugjV0TAIKgblF0+eWX09PTw/r160mlUvz5z39mjz324J133tE9v/7663zi\nE5/YEUlZmV4buTlL4OGxM5PJZIbJHjnk3ekU4X8B/3LYYYfVHXbYYQB0dnby2GOPcc0115BOp/ng\ngw944IEHWLx4Mf/6r//Kaaedpt075PP5ytNpdYrS5/MN0xrM5LuuvJRwYoxj/p3LTz7xODUocjKF\n6WSF4WR9PxON2/mm2Qhosocq7wNmf46aZUVRPgDmABcBNwEEg0EOPfRQNWSfj9WrV/Ozn/2Mgw8+\nmAMPPJBLLtFNB5JVVVUMDKhxJxIJy4EsEdG7Ax4eHvmxYsWKEfJUXV2NUfaASqt7HRsLKYrSBtwK\nICL7VVdXf/p3v/ud7r7//vvz1FNPmW+7BjgnEAiEEgnLBUwjKERtOVFGIMU0x3W6UtAcf670OXHz\nWvzJi1ZeAoEABtkLYrG70Fg76KcCvwJeB24GWoAVwDuoFcXjwKcVRfk3IB4KhYjFLPczKAoTNbhl\nHPixmuc3TNXkdMuF1SCeMV7z9KQ57Fz329kZWD2X2c2jNNHKk0n2LHf2GZNBv6IoncC5pstXZz9m\nYsFg0NVKwMPDQ8UkeyErP26s6onW1tbS19fH0NAQ5eXlLkSZm4leV2AVrxNb/tGumd3s1hrYDSQ6\nSXOx824q2fRPNCbZs9zu2435upi2y7DTcQEPD4/CYJI9S03AjUqgv7q6Wv3Rb3kKkqsUc7zAzkBH\nw6m5r52dfy7jn7EYGzlZc+A2O5vBUjHz2yR7lucAulEJJLSNDj1NwMPDXUyyV7iBwTzRK4Fk0nKz\n0ymJuf/tdCwg15Si3diAXUtiNbVoFaaT+ArJVO/35/N8xcwDk+xZbvntRiXQqx0/1tXV5UJ0E4dT\n68BcA225lv8aw7YKczQ7gXwqIreEcqzrKiZLpVGMdI4lD0yyV2/lx43uQGd9vRq3dla6h4eHO5hk\nr8bKjxuaQCx7Iqp+QupUw2oKLld3wK7VtrPXN4aRK36r+IzX7NJpFWYptbqllJaJYix5YJK9iJUf\nNzSBbk0lmaqVgIdHqWKSvQnTBKJTfUzAqtV3Mu1jN16Qyzw4V3xmP7n+j3afUWPxyJ9SG7swyV6d\nlR9XNIFQKERFRYV+VrqHh4c7mGTP0mLQFU0A1BppZxkYHE9/2klL4mTcwElrP17z4VKjFNJZqLgL\n+SwG2bPUBNyoBLoBIpHIuCuBUhusMmM11ZfPUuJ81XPjgJ+TQUMn6xBKQZDGykSluRh5VsiwDLI3\nYWsHVFOlYNCzGPTwmAAMsjdhFoN9oKok4x0TmCytk93AoNEtn4E9Oy3B6j6n6v1EGQsVEifGUMV8\nrlLPM4PshUWkXFGUIaN70TUBRVHSQLympsabIvTwmABMsjdiJaFbp4TGQqFQ0GpjkcncB7XDboDO\n7ObE3Njoz+q603UEue4bLY5Sxsn+Cjszpt2FQmQH6zXc2v+7z+0txjw8PFRMshc2u7ulCfTU1dXR\n2dk5oh/qZLXdZMRJ3z7f1juXBuEUJ33mfGYzPCYHJtkbsYjILU2gs76+noGBgZzLiafaRhK5Ngtx\nsuzXuHGIMQzzNfO9Vm75dA+M4RifwWNyY5K9EesH3KoEopGIGnd3d7dLUXp4eIBqJwC67I0wGHKr\nEtje0tICQFtbm0tROqfQW2qNtt2Xk/jMrb+dH/Nvp3FYxeNU8/AoLMXc1s0ke9PM7nlXAiKiiMh3\nROQjEflYRM7NXvuSiKwVkYSIbBaRfxcRrf+xefr06eqPzZvH/DAeHh75Y5K9GWb3vAYGRT3f/MfA\ntw2XfwGcASw2XPMDlwJfFpFDgDatNmpvb88nSldw0so5WYuf7+BmPv5Hi8/JjkR2hkFTYe3AZKaY\neWuSvUaze76awIXAt1OpFCeccAKrV6/Wri9OJBJccMEF1NTUMG/ePF588UWAXYEbgB5tSaO3ktDD\nw11MsjdiTwHHmoCIhIGfiQhnnnkm06dPZ+nSpbr797//fZLJJFu3buXFF1/klFNOYcOGDVRUVCwH\nVvv9fqDwm4261UqZw3e6ys/Jir18V/zZPbPVlOtYTINLfbGW20xmbcgke36zez7dgfMAnnnmGRKJ\nBL/5zW+GZcxvfvMbNmzYQDAY5Mgjj2TGjBmsW7eORYsWAcyorFQPRDWckFoQSvWlOLF/cGLBZwzP\nicWgXVrM8Y5mN1DodQWTuWKZrOkGMMneiLMH8ukOnAZwww038MorrxAOh5k1axYvvvgiHR0dhEIh\nfSoC1MEIw0xAdVlZGX6/37Ma9PBwGZPsjVg7kE8lsBfA2rVrWb16Nclkkn/7t3/jO9/5DlVVVSNO\nF4rFYjQ0NGh/20Bd0hiPx1m/fv2waTTts2LFijE84uTA6VRfrmu5phu1/+Z4jNN8udys0md1X6GY\nzK3peCm2BeaKFSssZWru3LnADtnDYjlxPpVARvsxZ84cFEXhmGOOYdu2bYTDYRRFYdu2bQCk02le\nffVVZs2apd3yMUBVVRUDAwO6euLh4VFcUqkUsEP2gBHCl08lsBFg0aJF3Hnnnfo4wOLFi1EUheXL\nl3PBBRfw8ssvc/HFF3PggQfS2toK8AHwf1pC+vv79fPRJjtODILM5sN2prpWLbqdGa+VdmHXyts9\ng53RkGcsNHas3q2bGCuBrLY+duETkUtFRDZu3ChHHnmkzJw5U5YtWybt7e0iIpJKpeS6666TRYsW\nycUXX6xfF5Fvi0iLiMg+++wjy5Ytk6kCIGrWWF/P5Wblz+zX7GYXj5M0jhZPPmF6OKdU8tMgeyMM\ndRx3VETEDzwNzM+j7ngIOAl1+WLnAQccwB577MFDDz2URxCTE6cmu+Z7rK7l8p9vOsSgpdhdG2tc\n+aSnWOFPBiYiDwyyF1UUZdgiIsdThIqiJEXkCFRLwOXAE8A84H+AlaiVw4XAe8A/Ab8GViuKIiIy\nBFBeXs7Q0JBF6IVhoguY3XJhK3/5huXEfz5hTsRy4bFWZFOp0piIZzHIXrnZLS+zYUVREsA12Y+Z\ntdmPFWkAn89HOp3OJ0oPD48CYJC9ETLv1qYiBdMEnFjLTRSjaQBO3PLxMx7/5nvyvT+f7oPTNQ52\nXSGn6XCaFnPco9032bHTBNxaSuxpAh4eE0gpaAIZUC2XMpnMaH5tmcq1tfHZpkLrNJ5W2ug311To\naGFbtfpTLY+dYpC9ES/FLU1AwNu3zsNjorCTPbc0AY88KcXWycl0pbFvb+5z22kGoz3veKcyx7oq\nsxQpdHrdqgTKADKZDD6fe/VOqb7cUk1XIRhN3c5nmtJO5XcStpWb3X2ThbGUGzvZc6s74AN1TYGb\nlYCHh4eKQfZGjMy7JZETUgk4qTEnwhBlsmoAxchPO43AiZbgND4nXRO3cVMjtKsE3NIEKgFvBaGH\nxwRhkL2U2c2tZjkE6h4DodCIPQ0Kwlhb9MnaKk80ufr2VgZBuVZBWoWTK66xDuw5CX+iyoCb8Rpk\nL252c0sT8IO6x1kgEHApSg8PDw2D7CXMbm5pAlUA/f39VFWN2OKsIJRaiz6VZwBg5HNZ9dGdzAQ4\nXWCVT5hOwrHyN9bpxsmAQfYmrDtQKyL09fURDo84FLVkGc+g4WQqIPliJehmd8hv6bIVTqYIc6VP\n8+ttk6am1yB7fWZ3t7oDtYlEgnQ6TV3diKPQPDw8iohJ9nrM7m5pAvXaQaSTqRKYbDV+sXEyGGfX\nHbC7zxx+Ljcra0S7rkmu+6faHgUaVvlrkr0RlYBrmoB28pB2GoqHh4c7mGRvwjSBhs7OTmByaQI7\nE/lOpeVqtUfrhzsJ34mJ71j3DCgFI6FiYpUvJtnrNru7pQnUd3V1AdDYOOI8RA8PjyJikr0Os7tb\nmkBTR4cad319/ShePSYCJ1NnThbfWO1WlO9CIqfmybn85rtz01QdG9CeyyR7IzQBtyqBRu1I8ubm\nZpeinHq4PT893m3Q8lHn7boRow0YOrFDsBvAnIoYn80ke21mv651Bzo7O/H7/foJqR4eHu5gkr1O\ns3tBNQERKQMagbiiKEYb5elbt26lpaWlkNEBk9N6a6xM9DM6UfVHa2HzeYaxGhQ5CcOoeUx0vhYb\nk+xtN7vnpQmISJmIfF1EXhaRlIh8KCLfEZFKEbkY+AequhETkYdFZO9sxdCydetWpk+fPt7n8fDw\nyBOT7G0xuzuuBEREAX4B3A4chLo8+BPAj1Htkf8DaO3s7NTOPzseWA8cDJT19PQUZVDQamXaRGG1\nIeZUQnu+XKbCTvrZdnlkfpfG/3bv2RymlV/tWq6VjGN5b5PlfRtkL64oyoi1A/loAp8Czo7FYpx8\n8snMnDmTY489lt7eXj2i448/njlz5rDXXnuxatUqLfybQB2cMBxV7uHh4RIG2euycs+nEvg8wE9+\n8hNmzZrFhx9+yBFHHMF1110HwPe+9z1mzZpFW1sbTz31FBdddBGxWAzgk6AOTkx1G4F8WqvJyFi1\nLqtnN2oVudyM8Ro1DeMnX8MkY5x22oETSkkLtcMgeyMOI4X8BgarQG3xw+EwQ0NDtLW1UVlZiYhw\n33338d5776EoCnvuuSf//M//zHPPPcexxx5LJpMhmUwWbUORycBYC8tkHLyym+rLZ0rRyQpFqwE+\nq01MJmM+FgKT7MWs/ORTCVQAfOMb32Du3LncfPPNVFZW8tprr9Hb20tlZSU1NTW658bGRn3hgmax\n5JkMe3i4i0n2RqwbgPy6AzMAVq5cyUknncSGDRs444wzOOusswgEAsTjw3ct6unpoampCYBt27YB\nMG3aNG6//fYRqqD2WbFiRX5PWGBKUWUvJZXTavDPLn1O1G3z1mG5VP987jP6cfI8pfje82XFihUj\nnueSSy4ZJnvAVqt789EE9gG4//77eeaZZ2htbeXqq6+mqakJn8/H9OnTWb9+PXPmzGFgYIB169ax\n7777AsNtlzdu3Dj2J/Xw8HBMJBIZdd0A5KcJtAPss88+rFy5ko0bN/KrX/2KefPmoSgK559/Pl/9\n6le5//77+cpXvsInP/lJWltbAbo0LSEYDGqDhSWFueXxGI5VC5urFXU6UJdvq51r+tBqatBKg3Ci\nsUw1gsEgRtkjx5iAY0Tk2yIi27Ztk0suuUSOO+44ueSSS2Tz5s0iIpLJZOT3v/+9nHbaaXLjjTdK\nPB6XLH+79957BZC33npLShFA1Ef0sELLn3w/+YZpdrO7z/h/vM81lTHJ3kVWsp1Pd+BnwJEtLS2f\nu+WWW0Y4KorCsmXLWLZsmfHybUBNX5+6rVmp7i8ok6wVmOiFRHaj/FZuTs1/raYKzfGbwzbe5yQ/\n8t2jYLLPKphkr9fKj+NKQFGUARE5HjgR1WZgHao14B+AB4HpwNmoZolzgLsURXlJRB7SugA78xRh\nISmFAulkNaDmpn2PdTWg0/UIToTZ6r5c6RnNTyEpVnwm2Rv3FCGKomSA32c/Zj4GVlhc17cWM04h\nenh4FB+T7I1PExgHgVgsRlVV1ZQ/jNRYm092NdIOO1XaqkU3agT5rAx0mof5dCM0N/OAYq40uEWx\n4jPJnqUm4MZ+AkHv5CEPj4nBJHsjTh8CdzSBcG9vb8kOChYSuwGsqYBV62nVj8/VMlv5t9MqrLQL\nK3INRNrd53Rg0AmlrPWZZG/CNIHwZDt5yMNjqmCSvRGnD4E7mkBVKpUq2hmEHtZM1DSi01F4u2nE\nfBcJ5brP6WzEePOokHlc6Pdmkr1+Kz9uVALlQ0NDU35QcKLJpUoXAyvBs3Kzw4llod1qQKuwnAwo\nGv2MNe3FpNBpMMle2sqPG92B8qGhIcrLy12IysPDw4hJ9oas/LjRPPvS6fSkrQRKedDHiBvpy3cq\nzW760C6cfLQEuynJ0cLaGTDJnmUl4MqW4yJCWZlbu5t7eHhomGQvY+XHtY56qdXITlv4fA1XdgZG\nM7TJZ4DOeE8+ZsZWfXsn78bptONUwqA1WT6wa81zJmNZCXl4eBQZTfZE3f5/BG5oAlJWVqak05YD\nkxNGvq3AZGs1JkpzsTMasjPjNf7OZ4ZjtMVCdubCO4N2V1ZWhkH2yrDoErhRCWTKysrKPU3AXYo9\nRWgWKjvBdWKv73TaMV/BzeV/rN3TQnQj3Kx8ysrKjFq4pSbgRndg0OfzUWqagIfHzoBJ9iwbfTc0\ngaTf769OJpMuROVRTOyMd5y08qO1vvnsB+DEMtFu1eJYW+HxaiDjiXss+P1+DLLnx2IRkRuaQDwc\nDpfk3oIeHlMdk+xZLuBxQxOIW21JPtnZmaaYrMjXTl9jrFOuufr0VuMTVvFO9JZsE4VJ9izX87uh\nCQxUVFQwODjoQlQeHh5GTLJXYeXHDU0gVV1dTX+/5QImR5TiVM541p6X0nPkQz7GP1b3jRaWkxF8\np2Fp10vNSM1tTLJXbeXHlUqgqqpKO658p2eyVgDgTAUfTfDGMmhoFeZYNxzZmXj22WcxyZ7len5X\nugPaoaXeNKGHh3sMDAxgkr1KK3+jagIiogB7A/2KonyUvRZE3VY8CazP7kKMiFQBc1FXK72lKMoQ\n2e4AqBsc+Hy+vFXinbk2LyWsLP5G8zeWsPOJz2ma8sGuGzKaWz6bnljF52Tg1OmGtqlUCpPsWWoC\ntpWAiOwB/A6Yn/1/B/B34N/YoUW8LyLfAGYB1wDB7PU2EbkU6K2trQXUQ0qzxyF5eHgUmd7eXkyy\nZ3kseM7ugIjMA94A5m/evFk7yeQs4CoRKXvttdf4+OOPAfYEHgVuBoJvv/0277//PkALcBfQ3NDQ\nAOw4mHQiW3arHWU8PGBH2dCMoYwfo1uh4tEYLUwr/07u6ejowCR79VZ+7cYEVgLBRx55hDlz5rBm\nzRo9sMMPP5yvfe1rLFmyhB/+8IeAurXx8ccfz0knncRxxx3HRRddpCX0MO3kIc9gyMPDPfr6+jDJ\nnqUabtcd2FtEuOOOO9h11131mujaa69l4cKF3HjjjcRiMQ444ABOPfVUHnroIcLhMG+99RaDg4Ms\nXLiQF154gUWLFpX7/X4ASsF0uJBaSClOXU51JirP89lBKZf/fA2rtP9OxyOsrj/99NOALnt+qzjt\nNIHHFUXhwQcfpKmpSd+x9IEHHuD//b//B6jnmx199NE888wzPPjgg1x22WUoikJlZSVLly7lqaee\nAvRjkaeE1aB5EwuvAnAHs6o+kemw6i6YpzHNKrz2386PFeaKxfw/VxdGwyR7loeB2lUCZ2s/Ojs7\naW5u1n83NTXpnmpra+nr62PLli3ssssu+vVIJEJvr3r0mea/ra3N9oE9PDwKi0n2Wqz82FUC+kEF\n3d3demB+v59EYsdCpPb2dlpaWohEIvrAn3Z92rRpAPrgRHd3N3/4wx+G1YrGz4oVK8bynK4y0S1R\noZksA6Vu5LtVi25FrjyzarWt/JqfxWn3wE4bWrFixQh5Wrp06TDZAyJWcdhVApebIweYP38+f/3r\nXwF1J9OnnnqKgw8+mPnz5/Pkk0/qD/nnP/+ZBQsWAOoihoqKihFahIeHR/Ho6ekZJnuApfBZDgxm\nDYQuExGuv/56EokEl112Gf/xH//BN7/5Tc4991zefvttnnjiCQ466CD22GMPLrjgAj772c/S3t7O\nq6++Snl5OQsXLgR1Z6Gy6dOns2XLFlpbW20TPplt6ycjO0teOxlQtDPw0bAzmTb7M35blWs7QyKn\nBkHaNSstesuWLZSVlaHJHmApfDmfQkTaRKR51apVJBIJfD4fX/7ylykvL+fNN9/k3nvvZe+99+bU\nU0/VTzj58MMPufPOO5k+fTpnnnmmZq10M/CtefPm0drayqOPPporSv3hd5aC6eEexagE7LCyCswn\nbqeVgF04NTU19Pb2YpC9hKIoI6YJ7SqBk1CtBTU/W8lRkwD3o1oM/pPp+i+Bi4H+JUuWADumLDym\nLpN16tTpFJwTYc4VjtM05GumnOt+AJPs+bLm/Do57QQURblPRNYBh2T9fQx0AB8AX0TdpmhP4CFF\nUd4XER9wPOp5Z/sCjyqK8lY2QdGGhobaDRs2jJoJHpOfsQh/qWuA+c7xO/FvvM/cfTCHYwzLSkvI\nFW93dzcNDQ0YZK8RGDZNZ7t2QFGUD4EPLZzutvCbBh7M/l1tcm5rbm6ufe655+yi8/DwKDBtbW00\nNzdjkL0G8qkECkh3XV0dXV1dZDIZ70iyIjJZVfGxpreQGoTTaTu7+53se5DvCsNc2oFTQyOT7DWY\n/bklje1NTU2k02ltIZKHh4dLmGRvwiqBrvp6dQFTdr5y3EyUkUupG9dMNWOm0Sj0WhCzQY7d+7Yz\nBDIb7liFZWXya+WWK2y7Zze6mWRvwiqB7khENVbq6elxKUoPDw9QTfhBl71as7trYwKF0gQmus9b\niq1sqY+sTxbsjHGMWK30M/s199tHGxNwQq54c/nVwjfJ3girQbcqgXbNXHi8lcDOXNhzFZ6dOU/y\nxU4A7ab4jFN4+QzimcOzCjtXesxh2sVnNX2oYZK9EbsLudUd6NM2N/AGBj083MUkeyNOIXJLE+ir\nqakB0JcXlzIT3eXIhZvpmapdDKeDaeZrTqf1ck0RGrFS4+3us9MAnWglJtmrMcftliYQm0obi3h4\nTCZMsjfiKDK3KoGkz+ejoqJi2F4EpcrONs1mxMn0086E1bSecUov15SdVRmy24cg1xRhrv0DzPea\npyS1MD75yU9ikr0RC4hcGxMA9YRUb0zAw8M9NHkzyN6EjQmoW50Gg153wIJSGoMohTRMFE5G60cb\nN3Ayyu9kHwK7RUJGv3ZLkI0YZG+EJuBWJTAAUFlZycDAgEtRTh5KVfBKqXIaL06exShITgb4zGFb\nxWPnZrxutz4gl4DbLTM2Y5C9EUeRuTYmAOr+hKWw7biHx86GQfZGbDvuliaQAswnpLrKVGrV3CKf\nvCr1/M0nXaOt4c/l36nKn2/a8rUUtArTIHsjziN0SxPIAJSVlZHJZFyK0sPDQ8MgeyNk3i1NYAig\nvLycoaFhOxu51oKUagvllFKftivltDnFyRoAK/9mv1Y4Mf7JFY+T8HONF1RWVpJKpYyyV26+19ME\nPDymMJq8lYImkLMSmAotiBsUet18ocN0CzfS7nR1npNR+nxG8J36cTKVqVGSlcBkLHj54qnuxcNO\nhR7vc9mtATD6sVpZOBqj2RBYpcF8X76b2RjDspM9Vzf7K+UdeTw8pjJ2sjeqJiAifuBwoB94VlGU\nIREJA/ujGgG9qShKKuu3FViIOhD4jKIo2uGEStZ9p6gIJqKlncwq/ngp1DM7KZtOW3QNKwMku2lG\nO/LRguxWHJqxrQRE5BDg98Cu2Ut/EZGXgcsM3vpE5HvAHsA3TPd/F7iebCWQyWT004o8PDzcw072\nckqkiHwaeBzg1VdfZcaMGTQ3Nx8FHJXJZHj55ZcJBALMnTs3DNwC6gGljz/+OH6/nyOPPBJFUa7N\nxnFdNswJ0QRKvY9eCKb687mB3TSelS2/k70GrMJ3YuxjvMduKjJXmDaawIjpObsxgZ8C3Hbbbcyf\nP5833ngDUA85XLhwIZdeeilf/OIXOf/88wH1pJMFCxZw66238oMf/IClS5dqo5FXA82gVhKeJuDh\n4T4G2Rsyu9lVAjNEhHfffZcDDzxQr1nWr1/Pt771LZ555hnWrFnDHXfcgYjw05/+lCVLlvDHP/6R\nJ598klgsxhNPPKGF9TmAwcFBKioqCvt0BuzWaxeTXPF6TE6s1vAb1+2bsVrvb7XXQD77Axjvt9tr\nIFeY5nQaZG/QHKddJXC3oijccsst+P1+AgF1Q5KjjjqKU045BYA777yTBQsWoCgKjz32GF//+tf1\nBBx99NGsWbNGC2s2QH9/P37/iPULBcMuc4uJ2/F6lU5xMQqshtXGIbk28jBP+VlN++WaFrR6t2ZB\nN4aZq7LS7tM2GTXIXr/5ee1080uA80HdpbSlpUV3iMfjfPOb32TNmjWsXq0eO9je3k5jY6PuJxgM\n0tHRof8FdYMDbdNDDw+P4mGueAyyN2JXH0d2Aj09PTQ0qAeXJJNJFi9eTCQSYd26dcycOROAhoYG\n2tvb9Xu2bt1qrDh6AWKxGKFQiDfffNNSpVmxYkWej7pzMlEaj5GpqI04aYUhd+vrND+caAdGv3bh\nGGXHSqYOOOAAYIfskWclcKP2w+fz6eaH999/P3vttRc//vGPqa6u1j0vXryYRx99FFD7H48++ihH\nHHGE5rwVdqgk/f0jNBIPD48ioMmoXXfAshIQEQU4T0T4+te/Tnd3N6eddhrxeJzNmzfz6quvsmDB\nAl6KljQAAAmVSURBVObMmcPRRx/NwMAAF198Mbfeeivnnnsuhx12GLNnz9ZqoV7gnUwmQzKZJBgM\njrrZ6ES3MhMd/2SgFLSRQmM3wGfErB3YjRdYhTXW8mUXZi5CoRBG2QNGCJ/lmICiKCIibYqizDzv\nvPNYvnw5Pp+PQCDAhRdeyJIlS6ivrycUCtHW1harqKgI7bnnnrz22musWrWKE088kWOPPVYL7gdA\ntba3YCgUYsmSJVOuAHl4TCRXX301MHKQMRQKYZQ9svt9GrEbGDwHeGjevHna7qQbgVnhcJiFCxfq\nnmbMmPFr4AJQDz4866yzjGF8F7gZOCUWixkTYosbFcR4TTg9ph5WI/Z2La/V+EGu+0cbM7DSNszX\n7YyT7Mrx1q1bgTFUAoqi/FVEZgBzs/46gTbUE0xmAB2oewcmgO8BPcCBqNsX7Q78RVGU9mxiAlpt\npB2EMNG4Leg7s23/ZMHOYjBfrKwKx7rmINe6ACs3c6Vz9tlnc/nllwO67OWlCaAoSh+wxnS5C9iU\n45ZXst8vmK5HotEooB6JtGbNGg499FC7qIuOcSZiImclSiUd5vi9tOygVNPihEgkglH2sJgdcAUR\n+cFf/vIXAeTJJ5+Uq666SgBRnSYGLf6JTEMppaOYaRlLmBOdL8a4jWlxmi6zH6v78gnHKi1O4rv6\n6qvFKHsissIsn27tJxDQZgQCgYB3AImHh0sEAgGMsgeMED63VvM0dnWpWwvU19cTjUa9vnEOpuLY\ngRSwr+0Gdv13q2cxuuVaYWi+N5cfs5sxTDvMYWn3NTY2YpQ9IGq+1y1NoFUboWxtbdVHKz08PIqL\nUd5aW1sha7hnxC1NoKWtrY1gMEgwGKStrc2laMdPTU0NTU1NvP/++0UJX0QYHBzUt2Iv5RayEDh9\nvgsvvJBQKKSbqxeTV199lVgsRl9fHz09Pdx+++3D3E8//XR6e3uJx+MsWrSI559/fpi78Znyae2t\nsMqfXHlm1Fispis1t8cff1yXPdQZvuF+HaVsnIjIR1/96ld3e/LJJ/noo49455132Hfffd2IOifL\nly8nHA7T1NSkG1rYkclkSCQSdHZ20t7eTjweJ5lMEovF6OzsJBqN0t/fz8DAAKlUiv7+fgYHB0kk\nEvT19ZFMJkmn02QymWE7LiuKgs/no7y8nPLycqqqqgiHw9TU1OD3+wmFQtTW1hIKhQiHw4TDYYLB\nIPX19UybNo36+nrKygqj0AWDQdLpNOl0esT5EIVg3bp1zJ8/39ZPPB6np6eHaDRKb28vXV1dRKNR\nEokEiUSC/v5+YrEY0WiUeDxOLBYjkUiQSqX0/JUcVn7aEd0VFRX4fD59dWwwGCQUClFTU0Ntba3+\nXVVVRSQSobm5mZqaGgKBwLD9MJLJJNFolO7ubrq7u1m0aJHudsMNN9DT08OPfvQj/dr8+fOJxWKs\nX79+1Lxqbm4mEAiwadOmYddPO+00QqGQnt6rrroKUM35Q6EQwWCQQCBAJBIhHA7T2NhIMplk3333\n5aOPPgLYU1GUD4blzaipKQAi8uqXvvSlf960aROvvPIKV155JbW1tbS2tlJfX08gEKC2tpba2lr9\nhQSDwYIVbiPpdJp4PE5vby+xWIz29na6urro6emhq6uL9vZ2uru72b59O52dnbqQb9myhcHBEUux\nR6AoClVVVVRXV1NRUaE/j9/vx+fzUVZWpn9EhEwmw9DQEENDQ6TTaVKpFH19ffT19emFezT8fj8N\nDQ00NTURiURobW2ltbVVv9bY2Eh9fT0NDQ3U1tZSV1dHIBAouGl0JpNhYGBAF462tja2bNlCW1ub\nnr/d3d3EYjFisRi9vb26MGuC7ySPQT1qWxOGQCBAZWWlnr+KoujfgJ7P6XSawcFBBgcHSafTJJNJ\nEomEXpk4obKyksbGRmbMmEFjYyO1tbXU19cTiUSIRCJ6hV1XV0ckEqGurk73N569NEREL7PxeFwv\nI319fUSjUaLRKH19fcTjcRKJhJ63WsXZ1dXFK6+8AhBRFGXYuIBrBvK//OUvG//7v/+7vbKykpde\nesnRDIHf76e6ulpvBTXB8vl8w1609pK1z9DQEJlMhlQqRTKZZGBgQM+40QqZz+ejvr6exsZGGhsb\nCYfD1NfXM336dBoaGqirq6OlpYVgMKi31FohqKqqKvjOSel0Wn/52kuPx+N0dHSwbds2Xag6Ojro\n7Oykq6tLFzy7w199Pp+ep9pzaIKkCRMMFyCtohoaGmJwcFDXelKplK4Z2VFZWUlDQ8Owlqy2tlZv\nuTQhqq+vp6amhpqaGiKRCPX19XqrXV1dTSAQoLx8xEE640JEdA1DE6BUKkV3dzft7e16pRyPx9m+\nfTubN2/WtRStERlNe9IahKqqKqqqqvTyXVFRQXl5OSJCOp3WK1JNoLX8daKdKYqC3+8flrfBYJAz\nzzxzxjnnnBMCNiqKMkxNcq0SEJHZwLva/1gsxrZt2+jq6iKRSOiZr9XKWo2nqdyxWIxkMqnX4maV\nT6sUysrKKC8vp6ysTM/oyspKgsEg4XBYz5SamhrC4TANDQ00NDQQiURoaGggHA6bNZAEqoHUVqAd\n1TKyDdXyKok65dIFdKMevDqQ/e5H3cUlgWqgkQTSqHu8ZRRFyYi6UKsM9WioctQxmmogDISAQPa7\nxvQdAhqBViCS/d8E1GM4fz6ZTLJ9+3a6urro6Oigu7t7mPra19dHKpUikUgQi8UYGBjQuwPGvC0r\nK9O7LNp3RUUFlZWVVFdX6/mrVSR+v5/a2lqampqYMWMGLS0t1NXV4ff7rbSPvmweRrN5qH33ZvO6\nN5u/iay/VDbvo9nvePYzYMhfyeaxZMuels8+oCL78WXz12/I60g277XvaqAh+wkb3kcLMD0bDtk4\ndO3G2EXs6uqis7NT13RisRj9/f16A5VKpRgYGNDzu7y8nMrKSr2c+v1+qqqq9MZG6w5qXcNwOKxr\n0ZrGmUPDG6EBaLhZCSiohbgBtQC3oBbiWtSCq2V8yPAJsEMotJdSgSowVn2FTPYzlP1OoQrfAGpB\n0QpcArVw9aEWuHbUQtVp+nQoilL4znEREZEq1ALajJrXLUAdap7Xoeaz9gmjmnkHUN9BFapwmNUZ\nQRWwIcP3IGq+9rMjf2PsyPNeVNPyzcA21HzuQRXovqzf2GTLX41sedYqjvrsdx07ymo9atluyP6u\nY0dZrwYqUSsg7bcmi0Oo+alVeknUPNUaoyg7ym8MNS97UfNXewfad48WhqIo6WLkg4eHh4eHh4eH\nh4eHh4eHx+Tl/wM8Vj/PTskoOgAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "data = scipy.io.loadmat('network1.mat')\n", "A = data['A'].astype('float')\n", "plt.spy(A)\n", "comm = data['Comm']\n", "# Continue here\n", "#\n" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def MCL(A, tol, p, alpha):\n", " step = 1\n", " col_sums = A.sum(axis = 0)\n", " T = A / col_sums[np.newaxis, :]\n", " M = T\n", " while(1):\n", " print 'step ', step\n", " step += 1\n", " # Expancion step:\n", " M1 = np.linalg.matrix_power(M, p)\n", " # Inflation step:\n", " M1 = np.power(M1, alpha)\n", " col_sums = M1.sum(axis = 0)\n", " M1 = M1 / col_sums[np.newaxis, :]\n", " M1[M1<=tol] = 0\n", " if np.linalg.norm(M - M1) == 0:\n", " return M1\n", " else:\n", " M = M1.copy()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.9" } }, "nbformat": 4, "nbformat_minor": 0 }