\n",
"The following code may take quite some time to execute\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": 77,
"metadata": {},
"outputs": [],
"source": [
"n_runs = 100\n",
"p_list = np.linspace(0, 0.1, 10)\n",
"nb_infected = [repeat_spreading(G, p, n_runs) for p in p_list]\n",
"nb_infected_random = [repeat_spreading(G_random, p, n_runs) for p in p_list]\n",
"nb_infected_random_config = [repeat_spreading(G_random_config, p, n_runs) for p in p_list]\n",
"G_scale_free = ig.Graph.Barabasi(n=G.vcount(), m=int(np.mean(G.degree())))\n",
"nb_infected_scale_free = [repeat_spreading(G_scale_free, p, n_runs) for p in p_list]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us plot the results"
]
},
{
"cell_type": "code",
"execution_count": 78,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 78,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEGCAYAAABLgMOSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd3xT1f/H8ddJ0qZ7l1LKKJQlS0bZW3CgAgoo4gJRERQV9KcoCCiCAxAURWQpQwX9ylJERZbIpgiykRYKFOigeydNzu+PpFilQoGmN23P8/HoI8lNcvsGSj6999zzOUJKiaIoiqJcL53WARRFUZTySRUQRVEU5YaoAqIoiqLcEFVAFEVRlBuiCoiiKIpyQwxaB7gZQUFBMjw8XOsYiqIo5cq+ffsuSSmDb3Y/5bqAhIeHExUVpXUMRVGUckUIcaY09qNOYSmKoig3RBUQRVEU5YaoAqIoiqLcEFVAFEVRlBuiCoiiKIpyQ1QBURRFUW6IKiCKoijKDSnX80AURVGcjSwowHTuHKZTp8g/dQqsEp2XJ3ovL3ReXug8PdF5ev29zdMT4e6OEELr6NdNFRBFUZQbYM3LwxQbS35MDKaYU7bbUzGYYs8gzebr25lOd7m46L0KC4y92BQpNJe3e3naX1v42P5eT0+Eq6tj/sDFUAVEURTlKiyZmZhiYsiPOUX+qb+LhTkuDgoX5NPpcKlRHWOdCLy6dcO1TgTGiDq41qmDcHHBmp2NNSsLa1YWlqws++NsrNn2bYWP7a+xZmdjyczAfPHi5cfW7Oy/v99VCFdXe+GxFRa9p+c/HuvcXUrt70YVEEVRKj0pJZbkZFuRiIm2FQl7sShITLz8OuHigmvt2rg1aYxvnz62IhFRF9fwWuiMxv/cv85ohICAm8totWLNybUVnSsKUo690BRfkAoSE7CcOIw1Mx1r3nUeHV2FKiCKolQa0mrFfOEiplMx/yoWp7Cmp19+nc7DA9eICDw7dMA1og7GiAiMdergUr06wqDNx6bQ6dB72U5xlUhBPsRsgsMr4PhuMGeDZxVofB/cc6JUMqkCoihKhSPNZkxnz9rHJU79XSxOxyJzcy+/Th8QgLFOHXzuust+NBGBMSICQ0hIuRzUxlIAp3+Dwyvh+A+Qlw7u/tB0ADTpD+GdQKcHppfKt1MFRFGUckdKiSU1lYL4eMwJCRQkJGC+GH/5yifTmTNQUHD59YZqoRjrRODZuvXf4xMRERj8/TX8U5QSqxXO7rAVjaNrIOcSGH2g4T22olGnG+hLb9yjKFVAFEVxKrKggIKkJFtRiE+gICEec0KirVgkJlAQbysYV1zppNfjWqMGrhERePfoYR/EjsBYpzY6zxKe9ikvpITz+2ynp46sgsyLYHCHBr1sRaNuT3Bxc3gMVUAURSkz1ry8vwtDYgLm+HhbQUgsLBYJFFy6ZPutugjh6oqhalVcQkJwb94cl6ohGKqEYKgagkvVqhhCQjAEBmo2PlEmpIT4Q/aisRLSzoLeFerdAY3vh/p3gdGrTCNV4L9tRVHKipQSa2bmv44abEcL5oR4CuxHEJYiA9WFdN7eGEKq4BJSFWP9eriEhGAIqWorEiG2L72fX/kckygNSSdsRePwCkiOBqGHiO7Q7XXbaSo33yveYrVK0nLNXMrKJykz//Lt+YxUdqd+VWrRVAFRFKXEpJTk/fknmVu2UHCxsEjEY05MRObkXPF6fWAgLiEhuISF4d6yBS4htqMFW3GoiktIlYp3eqk0pJyyjWkcWQUJhwGBDO9EbqsRxIfdTkKBF0lZ+VyKSiEp6yKXMvNtj+2FIjnLRIG16JwRicH7EG5Vf0Dos0otpiogiqJcU8GlS6Sv+Z60lSsxxcSAXm87aqgSgrFhQ7y6di1SGGzFwVAlGF0Zzoouj6SUZOUXkJSZT3rCGYwn1lDlzFqCMo4AcNLYmM2ez/CDuTUnTnphOm4FDv9jHwadIMjLSLC3kWAvI41CfS4/DvIygksK38XO4o+knTQMaMjE9hNpStNSya8KiKIoxZJmM1lbt5K2chVZW7aAxYJ78+ZUfXsSPr16ofcq2/Pt/8lqgb0L4ffpoDOAVxXwqmq79a4KXiG2r8v3q4Dhvyf9laboxCyiEzNJyjL941RS4X1rZgI95C5663fSRmebm3HQWpv51kfY4dYZi0d1gr2N1Pcy0tHbSJCX6+VCUVggfN1d0OmuPL1ntppZfGQxc/+YixCCVyJf4eFbHsagK72PfYcVECGEG7AVMNq/z3dSyolCiNrAciAA+AN4TEppEkIYgSVAKyAZGCiljHVUPkVRipcfHU3aylWkr1mDJTkZfVAQgU8MwbdfP4x16mgd758uHIC1o+DCfgjvDH41ITMe0s/B+SjIvgQU0/7D3f9fhcVedIre96piG1+4gbGXcyk5TPvlBN//eeHyNiEg0NOVcA8T9+r30lm3lXqG/eiwkuFdl9jao7E0up9q1Rvyqocr+mKKQkntT9zPpJ2TiE6L5rYat/F629ep6ln1hvf3Xxx5BJIP3CalzBJCuADbhBA/AS8BM6WUy4UQnwFPAnPst6lSyrpCiIeA94GBDsynKIqdJSuLjHXrSF+xktw//wSDAa9uXfHr1x+vzp0QLo6ZR3DD8jNh0xTYMxc8gqD/Qtvlq//+sLeYITsJshIgM8F2m5VgKzKF98/utD1nyb/y+xjcihSaEHthKbxfpAB5BIHeQGq2iY83RbN0Vyx6nWBk97r0alqVKq4mAs5tQH90FcRsBGsBBNSBVi9D4374hDTCpxT+WtLz05m5byYrTq6gqmdVZnWfRfea3Uthz8UTsgTNuW76mwjhAWwDRgA/AlWllAVCiPbAm1LKO4UQv9jv7xRCGIB4IFheJWBkZKSMiopyeH5FqYik1UrO3ijSV64g45f1yLw8XOtG4NevP759emMICtI64pWkhGM/wE9jbHMfIp+AHhPB3e/m95uXBlmJ9uKSCFnx/7xf+Fxe2pVvFzpyDf6cMXuTYPXFKzCMhnXr4uUfAud2wV/rbQXKpzo06Wf7Cm1+Q0c3xceXrD21lulR00nPT+fRWx7l2ebP4uHiUezrhRD7pJSRN/t9HToGIoTQA/uAusBsIAZIk1IWThGNA8Ls98OAcwD24pIOBAKX/rXPYcAwgJo1azoyvqJUSOaLF0lfvZq0laswnzuHzssL37598evfD7emTZ33ctm0s7DuFfjrZwhpAg8ugRqtS2ffQthOa7n7Q3CDq7/WnAfZiZCZgCUzngNHj/PnsRO45V6ikXcu7bxyccs7APvW2440vEKg1RDbEVL11qAr3XX8YtNjmbx7Mrsv7qZpUFPm3j6XhgENS/V7/BeHFhAppQVoLoTwA1YBtxT3MvttcT+1Vxx9SCnnAfPAdgRSSlEVpUKz5ueTtXEjaStXkb19O0iJR7t2BD8/Eu/bb0fn7q51xP9mMcOuT2HLe7bHt78N7UY4rD3HNbm4IX1r8FuCG++tz+V4fCNurd6e1+++heZ1Av9+ndUKuam2oyOdvtRjmCwmFh5eyIKDC3DVuzKu7TgeqP8Aegd8r/9SJldhSSnThBBbgHaAnxDCYD8KqQ4UjjLFATWAOPspLF8gpSzyKUpFlXf0KGkrVpK+di3W9HQMoaEEjRiBb7/7ca1eXet413ZuD/wwChKPQIO7oddU8KuhaaTD59N596djbI9OpmaAB5883IJ7moZeeeSm04FnYPE7uUl74/cyaeckYjNiuSv8Ll5t/SrBHsEO+V5X48irsIIBs714uAM9sQ2MbwYGYLsSazCwxv6W7+2Pd9qf33S18Q9FUYpXkJpKxtofSVu5kvxjxxCurnj37Ilv/354tmuH0Jfdb6g3LDcVNrwJ+xaBTxgM/ApuuVfTSOdScpi+/gRrDlzA38OFN3s34uG2tXA1lO4pqatJzUtletR0vo/5njCvMOb0nEOnsE5l9v3/zZFHIKHAYvs4iA74Vkq5VghxFFguhJgM7AcW2l+/EFgqhIjGduTxkAOzKUqFIi0WsnfsIG3lSrI2bESazbg1bkzIhPH43nMPet8r2104JSnh0P/gl7GQkwztnoPur4PRW7NIqdkmZm+OZsnOMwgBz3WP4JmuEfi4ld0pNCklq6NX88G+D8g2ZfNU06cY1mwY7gZtTz06rIBIKQ8CLYrZfgpoU8z2POABR+VRlIrIdOYMaatWkb56DQXx8ej9/PAb9BB+/frh1rBsBlJLTXIM/PgSnNoCYa3g0RUQeqtmcfLMFhbtiGX25miy8wsY0Ko6o2+vT6hv2X5ox6TFMGnnJP5I/IOWVVoyvt146vrXLdMM/0XNRFeUcsaak0PGL+tJX7mSnL17QafDs3MnQl5/Ha/u3cpf+5CCfNj2Ifz+gW2G+N3TIXKoQwaeS8Jilazaf54Z609wIT2P2xpWYcxdDWlQtWyPgvIK8ph3cB5fHPkCD4MHb3V4i/vq3odOlN0ps2tRBURRygEpJbkHDpC+ciUZ637Cmp2NS62aBI8eje99fXEJCdE64o05vRXWjrZ1mW3cD+561zYxTyO//ZXEu+uOcTw+k2bVfZn+4K10iCj7+TA7zu9g8u7JnMs8R+86vXk58mUC3R0zIH8zVAFRFCdWkJpK+sqVpK1YienUKYSHBz533YVf/364t2zpvHM2riX7EvwyDg4uB/9w2+mquj01i3P4fDrv/XScbdGXqBHgzqxBLbi3aWixPaYc6VLuJabumcpPsT8R7hPOgjsW0Da0bZlmuB6qgCiKE7JmZ5OyZAnJCz/HmpWFe8uWhE6ZjPedd6H3Ksftz61W2L8Ufp0Apmzo/H/Q5f/ARZvB4LjUHD5Y/xer9p/Hz8OF8fc24tF2NTEayvb0mVVa+e6v7/hw34fkWfJ49tZnGdp0KEZ92TR9vFGqgCiKE5EmE6n/+x+X5nyG5dIlvHr2IPiFF3CrX1/raDcv8ZhtTse5XVCzA9w7E6poM9CflmO7smrxDtuVVSO6RTC8awS+7mU/OfFEygkm7ZrEwaSDtKnahjfavUFt39plnuNGqAKiKE5AWq1krPuJpI8+wnzuHB6tW1Plk49xb95c62g3z5QDW6fCjo/B6AN9Z0PzR0qtD9T1yDNbWLIzlk82RZOZX8CAlrYrq6r5lf0RUI45hzl/zmHp0aX4uPrwTqd3uLfOveXqtKQqIIqiISkl2du2kThjJvnHjmFs2JAa8+bi2blzufog+U9/rYd1L9v6WDV/xNaGxEGzs6/GapWsPnCeD9b/xfm0XLo1CGbMXQ25JbQ0euBev9/O/caU3VO4mH2R/vX6M7rVaHyN5WSuThGqgCiKRnIPHCBxxkxy9uzBpXp1qk2bhs89dyNKudmeJjIuws9j4OgaCKoPQ36EcMfMmM4x57Djwg52XdyFi84Ffzd//N38CTAG4O/mT0y8ZMFviZy4UECTMD+mDmhGx7radBqOz47n/T3vs+HsBiJ8I1h01yJahbTSJEtpUAVEUcpYfkwMSR9+SOavG9AHBhLyxhv4P/gAorzN3yiO1QJ7F8DGt8FqhtvegA4vlPoKgGl5aWyJ28LGsxvZeWEn+ZZ8PF08kVKSU3Dl2uz4go+vjgyjL9MOB7Ag2lZk/I3+fxcctwD8jH4EuAVcfs6llBo2WqwWlh1fxsf7P8YiLbzY8kUGNxpcavvXiiogilJGzBcvkjR7NukrV6FzdyfohecJHDwYnWc5vqqqqKKrA0bcBvd8YFs0qZRczLrIpnOb2HR2E/sS9mGRFkI8Quhfrz+31byNViGtMOgMnLqUyvQN+/nleDSe7vnc0cyLJjX0ZJjSSM1PJS0vjZS8FKLToknLSyMtPw1Z3KqFgJeL1+UC849iYwzAz81ebIz+l+97GDyuOPV4JPkIb+14i2Mpx+hYrSPj2o2jhre2DSFLiyogiuJglrQ0Ls2bT+qXX4KUBDz2GIHPDMMQEKB1tNJR0tUBr5OUklPpp9h4diMbz27kaPJRAOr41mFok6H0qNmDRoGNLn9gp+eYmb3lGIt2xALwVMfOPNu1Lr4eV/8t32K1kG5KJzUv1faVb7tNyUshLd9WbFLzUonPjudYyjFS81IxW83F7stV5/qPgmPUG9l6fisBbgFM6zKNO8PvrBhjW3ZlsiKho6gVCRVnZs3JIWXJUpIXLsSalYVv374EPz8Sl7Cwa7+5PJASjn0PP71mXx1wKPSYcFOrA5otFvZcOMCGMxvZfnELF3POARDm3oBw97aEGCIR5ipk5JnJyC0gI89MZl4BGblmLmXlY7JY6deiOi/dUZ8wB11ZJaUk25x9udAULTb/vp9uSqd9aHueb/k8Pq7aDNgXp1ysSKgolZE0m0lbsYKk2bOxJF3C67bbCB71YsWYy1Eo9YxtdcCTv0BIUxi4FKpHUmCxkpVjuvzhnpFr/seHve1xwT+2p+flkmo9Tq7hAFaPI+hcMpBShyU7goLM+yjIuoXjBb4cByAHb7c4fNxc8HYz4OPuQpifO7eEeuPv4Ur/ltVpVM2xH9RCCLxcvfBy9aowp6JulCogilJKpNVK5s8/k/jRR5jPnMW9VSuqfPQRHi1bah2tVCX+/jl+m17HCnzj/TTLc+8hbWk6Gbk/k22yXPW9QoC30YC3hxUXr5MUuB0ky/0gFnLRY6SWW3Ma+nSgWUA7QrwC8HFzwcfdYL91wctoQF/G7UWU/6YKiKLcJCkl2dt3kDRjBnlHj2KsX5/qn83Bq2vXCnW+G0sBSStfocqRz9ktGzPH72XyPcOo6W6gif0DvvAD39vNBR/7EULhNqsui6iE7WyO+5mdF3aSYcnH1+jLvdXvoEfNHrSv1h43g5vWf0rlOqgCoig3IffQIRI/mEHOrl24hIVRber7+NxzT/lY9e965KSQtuQRguN38K3+Xto88ymLqlx74pvtyqmf2XhoI/sS9mGVVqp6VmVA/QH0qNmDFlVaYNCpj6HySv3LKcoNyD91mqSPPiLzl1/QBwQQMnYsfg8NLH9rcZREwhFyljyIR1Y80z1e5LERYwnxKf5IQUpJTFrM5SunjqUcA6CuX12ebPIkPWr1oFFAo4p1ZFaJqQKiKNfBnJDApU9mk7ZyJTqjkaCRIwkYMqR8d8i9mqPfY14xjMwCIzMDpvLa04/h5/HPImmVVg5dOsTGsxvZdHYTZzLOANAsuBmjW42mR80e1PKppUV6xcFUAVGUErCkp5M8fz4pS79EWq34P/IwQc88gyHQ+Rb5KRVWK/z2Pvz2HoetdVlU/W3eHXIHHq62jwyzxcze+L1sPLuRzec2k5SbhEEYaBPahscbPU63Gt2o4lFF4z+E4miqgCjKVVhzc0lZ+iXJCxZgzczEt09vgp5/AdfqFWQuR3HyM5GrnkEc/5HvLF3Y0XAc0wa2wdWgo8BawFfHvmLuwblkmjJxN7jTKawTPWr2oHP1zk4110FxPFVAFKUY0mwmbeUqLs2eTUFiIl7duhE8ejRuDSrQXI7ipJxCLhuENekkk82PYW39DNP7NEGnE/yV+hcTt0/kcPJhulTvwgP1H6BdaDt15VQl5rACIoSoASwBqgJWYJ6U8iMhxJvA00CS/aVjpZTr7O95HXgSsAAvSCl/cVQ+RSmOlJLMX34h6cOPMMXG4t6iBWEzZ+DRqvx2TC2xmE3I/z1BtsnCMNMY2tx2Py/2qIfZambe/nksPLQQH6MP07pO485aFaslh3JjHHkEUgC8LKX8QwjhDewTQvxqf26mlHJ60RcLIRoBDwGNgWrABiFEfSnl1WcmKUopkFYr2b//TtLHn5B3+DDGenWp/umneHXvVvE/KKWEXZ8i17/BOX1NHsl9kad638bgDuEcSDzAxB0TOZV+it51evNq61fxc7vxViVKxeKwAiKlvAhctN/PFEIcA6524rgvsFxKmQ+cFkJEA22AnY7KqCgFqamkr1xF6vLlmM+dw1AtlNB338W3T++KN5ejOOY8WwfdP5exw7UDI7Ke5u2Bbbm9sT/v73mfr459RVXPqszpOYdOYY5Zz0Mpv8pkDEQIEQ60AHYDHYGRQojHgShsRymp2IrLriJvi6OYgiOEGAYMA6hZs6ZDcysVV+6hQ6R+9TUZ69YhTSbcI1sRPOpFfG6/vWKsy1ESGRdg+SNw4Q8+d32Yabn38ung1hi9o+n3/ZOczzrPoIaDeLHli3i6VNDLlJWb4vACIoTwAlYAo6SUGUKIOcDbgLTffgAMBYo7T3BFq2Ap5TxgHti68Toqt1LxWPPyyPhxHanLlpF3+DA6Dw98+/fD/6FBFX9w/N/O7YFvHsWan8VrhjH8bGrFZ4Mb8Gv8LNbsWkO4TziL71pMy5CK1cdLKV0OLSBCCBdsxeMrKeVKACllQpHn5wNr7Q/jgKKtLasDFxyZT6kcTGfOkLr8G9JWrsSano5r3QhCxr+Bb9++6L28tI5X9v5YCj++RJ5HKA+bxxDnUouX7stm4r4hpOWn8XTTp3nm1mcw6kt3FUGl4nHkVVgCWAgck1LOKLI91D4+AnA/cNh+/3vgayHEDGyD6PWAPY7Kp1Rs0mIh67etpH79NdnbtoHBgHfPnvgPGoRHm9YVf2C8OBYz/DIO9swlpWpH7j4/FBc/Ay2arOGDP7dwS8AtfHb7ZzQMaKh1UqWccOQRSEfgMeCQEOKAfdtYYJAQojm201OxwDMAUsojQohvgaPYruB6Tl2BpVyvguRk0r5bQdo332C+cAFDlSoEPT8SvwEP4BJSiWdGZyfD/wZD7O/E1B3CPcd6UKX6ccy+q9mXZGJUy1EMbjxYNTZUrotakVAp96SU5O4/QOqyZWT+/DPSbMajXTv8Bw3C+7buCJerL2la4cUfhuWDIDOB7Y3G89ifflQJ/4Fs/TFahbTizfZvEu4brnVKpQypFQmVSs+ak0P62rWkLltO/rFj6Ly88Bs4EP9BD2GMiNA6nnM4shpWj0C6+bK8yWdM+Gsv3hGLwODC+MjxDKg/AJ3QaZ1SKadUAVHKnfxTp0ldvoz0VauxZmZibNCAqm+9he+996DzVJebArZmiFvega3TkNXbMMH3SVbEL8St6lk6hHVmYvsJVPWsqnVKpZxTBUQpF2RBAZmbNpG6bBk5O3eBiws+d96J/8ODcG/RonIOiv+XvAxY9QycWEf+rY8wMDuA6Oz3cPf04K2O73JPnXvU35dSKlQBUZyaOTGRtO++I+2bbylISMAQGkrwqFH4DeiPIShI63jOJzkGlg2C5Gj+6PISw8/uJFf8Tj3PLiy4dxKB7hW0/byiCVVAFKcjpSQ3KorUZcvIWP8rFBTg2bEjVSdOwKtLF4RB/dgWK3oDfDeUXKFnZrtHWHZ2BdYCbx6sPZGJPQZonU6pgNT/RMVpWLKyyfjhe1K/Xkb+yZPofHwIePRR/B8aiGt4uNbxnJeUsPMT+HUCu6vWZ4KfNxfit1CQ3o7JXcfQv3ldrRMqFZQqIIrm8k+eJHXZMtJXr8Gak4Nbo0aETpmMz913o3N31zqeczPnwg8vknH4W2ZENGeFJRldlieWhOHMe+BBOtcL1jqhUoGpAqJoQprNZG7YQOrXy8jZuxfh6opPr174PzwIt2bN1CBvSaSfh+UPsyn9BJPr1CfZmoY+4zZk2h18PaQjzWuotuuKY6kCopQ5S1YWsQ8OxHTqFC7Vq1Pllf/Dt18/DP7+WkcrP87u4tK3j/Gep45fQoKp4RaG9a978RLhLH2mDXWreGudUKkEVAFRylzC5CmYYmMJmzkD7zvuqBzrbpQiGfUFa3+bwPuBfuQYDNxT7XFWb6lHNT8vlj7ZljA/ddpPKRuqgChlKuOX9aSvXk3QsyPw6dVL6zjli8XMhR9HMen8L2wP8qN5YBM6B77A+z+k0ijUh0VPtCbQS3XQVcqOKiBKmTEnJhI/cSJuTZoQNGKE1nHKFWtWIsv/N4APSQZPL16PfIW8lHZMWnOc9nUCmT84Ei+j+u+slC3VBEcpE1JKLr7xBta8PKpNnaoaHF6HU9E/M2R5D97VpdLStx6r7v+R+LhIJq09zp2NQ/jiidaqeCiaUD91SplIW76c7K2/EzL+DYx1amsdp9xYvf0dJp38Gg89vNN4OL1aPMubPxzhy11nGRhZgyn3N8GgV78HKtpQBURxuPzTp0l4fyqenTrh//DDWscpNzYf/46JJ7+mtdXAe32W4+NXn1HfHGDtwYsM7xrBmLsaqMudFU2pAqI4lDSbufDqGHRGI6FTpqgPvBI6EL+PV3ZNopHZwqw+34BPXZ5cvJffT15i7N0NGdZFtatXtKcKiOJQl+bOI+/QIcI+/LByrwh4HU6ln2Lkr88QUmBidvOXyfeqxxPzd3MwLo2pA5rxYGQNrSMqCqAKiOJAuX/+yaU5c/Dt2wefu+7UOk65kJiTyPCfh2Iw5/CZdwtM9R7lsbk7OZOSw5xHW3FnY7WGh+I8VAFRHMKak8OFV8dgCKlCyBtvaB2nXMg0ZTLi12dIz03miywdfv1mMeCLvVxMz2PRE63pEKHa1yvORRUQxSESpk3DdPYsNRctQu+t2mpci8liYtTmUZxKi2F2QhINHlzBU6tOE52UpYqH4rTU9X9KqcvaupW0ZcsJGDIEz7ZttI7j9KzSyrht49gTv4dJiUl0aP08Uw77s/lEEm/1aaw66ipOSxUQpVQVpKZyYdw4jPXrEzx6lNZxnJ6Ukml7p/Fz7M+8lJ5Lb/8mfOX2EJ9vP80THcN5tF0trSMqyn9yWAERQtQQQmwWQhwTQhwRQrxo3x4ghPhVCHHSfutv3y6EELOEENFCiINCiJaOyqY4hpSS+AkTsKalU23aVHSurlpHcnqLjyzmy2Nf8qjVkyFZeextNZUJa0/QvUEwb9zTSOt4inJVjjwCKQBellLeArQDnhNCNAJeAzZKKesBG+2PAXoB9exfw4A5DsymOED66jVk/rqB4FEv4taggdZxnN7aU2v5YN8H3OlRi1fOHCOhyzsMXZNI3WAvZg1qgV6n5swozs1hBURKeVFK+Yf9fiZwDAgD+gKL7S9bDNxnv98XWCJtdgF+QohQR+VTSpcp7jwJkyfjERlJwJAhWsdxejsv7GT89vG09mvAO0e3Y270AAN31sBo0LFgcFLlFUEAACAASURBVCTebqpXmOL8ymQMRAgRDrQAdgMhUsqLYCsyQOHssjDgXJG3xdm3/Xtfw4QQUUKIqKSkJEfGVkpIWixceG0MANXef0+t73ENx5KPMWrzKGp71+Kj0ydw8a3BsEsPcTE9j7mPRVIjwEPriIpSIg4vIEIIL2AFMEpKmXG1lxazTV6xQcp5UspIKWVkcLC6OsUZpHzxBblR+wgZ/wYuYVfUfKWIuMw4RmwYga/Rlzkmb7zSzzPL7zV+O5vPtAHNaFVLrcqolB8OLSBCCBdsxeMrKeVK++aEwlNT9ttE+/Y4oGiPhurABUfmU25e3vHjJH40C+877sC3b1+t4zi1lLwUhm8Yjtlq5rMafQg5+j17aj3DzOO+vNijHn2bq+KrlC+OvApLAAuBY1LKGUWe+h4YbL8/GFhTZPvj9qux2gHphae6FOdkzc/nwiuvovfzpepbb6pGiVeRY85h5MaRxGfHM7v1WOpsfI+UoEgGHe9A71urMapnPa0jKsp1c+RM9I7AY8AhIcQB+7axwHvAt0KIJ4GzwAP259YBdwPRQA7whAOzKaUgaeaH5J88SY358zD4q1Mv/6XAWsArW1/hSPIRZnaZTvNf38OCoH/CEJrVCGDagGaq+CrlksMKiJRyG8WPawD0KOb1EnjOUXmU0pW9axcpixbh//AgvDp31jqO05JSMmnnJLbGbWV8u/HcFrMLzkcxXv8SJs8w5j3eCjcXddGBUj5dtYAIITIpZiC7kJTSp9QTKU7PkpHBhdfH4hoeTpX/+z+t4zi12Qdmsyp6FcNvHc6DxmrI34exwdiTNblt+e6pSKp4u2kdUVFu2FULiJTSG0AIMQmIB5ZiO6p4BFAd8iqp+MmTKUhMJHzZ1+g81CWn/+XbE98y9+Bc+tfrz7P1ByE/60SiSzVGZzzMx4NbcEuo+v1LKd9KOoh+p5TyUyllppQyQ0o5B+jvyGCKc8r46Scyvv+BoBEjcG/WTOs4Tmvj2Y1M2T2FrtW78kbbcYi1o7BmxvNU1ghG39OS2xqGaB1RUW5aSQuIRQjxiBBCL4TQCSEeASyODKY4H3NCAhfffAu3Zs0IGv6M1nGc1v7E/YzZOoYmgU2Y2mUqhj+Xw9E1TDM9QNM23RjaMVzriIpSKkpaQB4GHgQS7F8P2LcplYS0Wrn4+likyWSbbW5QS8kUJyYthpEbRxLqGconPT7BI/08lnWvsMPamMPhg3mrT2N1xZVSYZToU0BKGYutV5VSSaV+vYzsHTuo+uZEjLVrax3HKSVkJzB8w3Bc9a7M6TkHf4Mn+d/0IbdAz4feLzP/kda46NUKCkrFUaKfZiFEfSHERiHEYfvjZkIItU5pJZEfE0PitGl4du2C38CBWsdxShmmDIZvGE6mKZM5PedQ3bs6+esnYUw6yFtiBFOfuAtfD9UgUalYSvrr0HzgdcAMIKU8CDzkqFCK85AmExdeeRWduzvVJk9Wp1+KkW/J58VNLxKbEcuH3T+kYUBDCk5uwrjnY5ZZejDwsRGEB3lqHVNRSl1JT2R7SCn3/OvDo8ABeRQnkzRnDnlHjxI26yMMqnnlFazSytjfxxKVEMX7nd+nXWg7ZPYlcr55ikRrNVzvfY92dQK1jqkoDlHSI5BLQogI7JMKhRADANWnqoLL2b+f5Lnz8L3/fnzuuEPrOE5HSsnUvVNZf2Y9/xf5f9xd526QkrOLnsRoTuf3Zu/Rv219rWMqisOU9AjkOWAe0FAIcR44jW0yoVJBWbOzuTDmNVxCQwkZN1brOE7piyNf8NWxr3i80eMMbmzrD3r8x49omLSFbwNHMLhfH40TKopjlbSASCllTyGEJ6CTUmYKIdSlOBVYwnvvYz53jlpLl6D38tI6jtP5IeYHZu6bSa/avXg58mUATh2NInzvZPa5tODeZyahU0vSKhVcSU9hrQCQUmbbl6cF+M4xkRStZW7aTNr//kfgU0/iERmpdRyns/38diZsn0Dbqm2Z3HEyOqEjKTUdy/+eJEe4U2PoEjyMrlrHVBSHu1YzxYZAY8BXCNGvyFM+gOoCVwEVJCdzcfx4jA0bEvT881rHcTpHko8westoIvwi+LD7h7jqXckzW9g57wX6yFhi7/iC8Go1tY6pKGXiWqewGgD3An5A7yLbM4GnHRVK0YaUkovjJ2DNyKDaF5+jc1W/RRd1LuMcz254Fn+jP3N6zsHL1QspJZ8vns+zuas5E/EI4R36XXtHilJBXKsb7xpgjRCivZRyZxllUjSSvmIFWZs2UWXMGNzqq6uHikrOTWb4huFYpZXPbv+MYA/bJc3zf9rFA+emkOwVQa2HPtA4paKUrZKOgQwXQvgVPhBC+AshPndQJkUDprNniX/nXTzatiVg8ONax3EqhcvRJuYk8kmPT6jta7t+ZM3+OOrtHIOfLo+Ax5eCi7vGSRWlbJW0gDSTUqYVPpBSpgItHBNJKWvSYuHCmNcQej3V3n0HoVP9mgqZrWZe+u0ljqYcZVrXadwafCsAf5xN5eDKqXTX/wl3vI0IaaxxUkUpeyW9jFcnhPC3Fw6EEAHX8V7FySXPX0Du/v1UmzYVl2rVtI7jNKSUvLnjTbaf386b7d+kW41uAMSl5jB18Xcs1n+Nqc7tuLZTre2VyqmkReADYIcQ4jtss9EfBKY4LJVSZnKPHCHpk0/w7nUXPvfeq3Ucp/Lx/o/5PuZ7nm3+LP3r29ZPy8wzM+KL7cy0fIjewx9D/89A9QdTKqmStnNfIoSIAm7DtqRtPynlUYcmUxzOmpfHhVfHYAgIIHTiRNUosYjlx5cz/9B8BtQfwPBmwwGwWCUvLNvPg6lzqauPg34rwTNI46SKop3rOdkdAGRLKT8Gkq41E10I8bkQIrGwBbx925tCiPNCiAP2r7uLPPe6ECJaCHFCCHHndf9JlOuWOGMGppgYQt99B72f37XfUElsOLOBd3a/Q7ca3RjXdtzlwjrlx2PoTv7MY/pfof1IqNtD46SKoq0SHYEIISYCkdjmhXwBuABfAh2v8rZFwCfAkn9tnymlnP6v/TfC1h6+MVAN2CCEqC+lVMvmOkjW9u2kLlmK/6OP4tXxav+Mlcu+hH2M2TqGZsHNbMvR6mz/Rb7cdYa12/9gi+cCCG4KPSZonFRRtFfSI5D7gT5ANoCU8gLgfbU3SCm3Aikl3H9fYLmUMl9KeRqIBtqU8L3KdbKkpXHx9bG41qlDlZdf0jqO04hOjeb5Tc8T5h3GJ7d9grvBdlnutpOXePP7Q3zhtxB3YYL+C8Fg1DitomivpAXEJKWU/N3O/WZWxxkphDhoP8Xlb98WBpwr8po4+zbFAeInvU1BSgrVpk5F567mLgCYLCZGbxmNm96Nz3p+hp+b7ZRedGIWI77ax6s+G2mc9wfirnchuIHGaRXFOZS0gHwrhJgL+AkhngY2YFul8HrNASKA5tjWEymculvc6K0sbgdCiGFCiCghRFRSUtINRKjc0tf+SMa6dQSPfA73JmruQqGFhxcSmxHL5I6TqeZlu5Q5JdvEk4v30kwfy9OmpdDwXmg1RNugiuJErtVM0Wg/rTRdCHE7kIFtHGSClPLX6/1mUsqEIvueD6y1P4wDahR5aXXgwn/sYx62tUmIjIwstsgoxTNfvEj8W2/h3rw5gU89pXUcp3Em4wwLDi6gV3gvOoR1ACC/wMLwpftIS0/jl4DPEDII+nysLtlVlCKuNYi+E2gphFgqpXwMuO6iUZQQIlRKWbiS4f1A4RVa3wNfCyFmYBtErwfsuZnvpfyTtFq58PpYpMVCtfffQxjUPFCwTRacsmsKRr2RV9u8ennb2JWH2RObwpYGP+J25jQ8vgY8AjROqyjO5VqfIq5CiMFAh3+1cwdASrnyv94ohFgGdAOChBBxwESgmxCiObbTU7HAM/b9HBFCfAscxbbW+nPqCqzSlbp0KTm7dlF10lu41qqldRyn8dPpn9h5cSfj2o4jyN02p+Oz306x4o84Pml+jvDj/4OOo6BOV42TKorzEbax8f94UohO2JaufRDbUUJRUko51IHZrikyMlJGRUVpGaFcyD95ktP9B+DZsSPVP52tJgzaZZgy6LOqD6GeoXx595fodXp+PnyR4V/+weONDLx1YRgioDYMXQ8G1dpeqTiEEPuklDe9Wty12rlvA7YJIaKklAtv9pspZU+azZx/dQw6Ly9C356kikcRs/6YRWp+Kp/2/BS9Ts+huHRGfXOAVjW8mWh5H2Ex2y/ZVcVDUYpT0lYmC4UQHYDwou+RUv57kqDiZNLX/kj+sWOEffghhiDVdqPQoaRDfHviWx655REaBTYiz2xh5LI/CPQ0srjBLvTbtkHf2RAYoXVURXFaJZ2JvhTb5bcHgMKxCcmVs8wVJyKtVpIXLMBYvz7ed96hdRynUWAtYNKuSQR7BDOyxUgAPt0czZnkHL6/zw2v9e9D4/uh+SMaJ1UU51bSS3EigUbyagMmitPJ2rIFU0wM1aZNU6euilh2fBnHU44zo9sMPF08iUnKYs5vMQxs6kez3c+BV1W4d6a6ZFdRrqGkEwkPA1UdGUQpXVJKkufOwyUsDJ9ed2kdx2nEZ8fzyf5P6BzWmZ41eyKlZPzqw7i56HnT+CWknYH+88Hd/9o7U5RKrqRHIEHAUSHEHiC/cKOUso9DUik3LTcqitw//yRk/BtqzkcR7+95H6u0MrbtWIQQrN5/nh0xyczrnIv73q9tl+zW6qB1TEUpF0r6yfKmI0Mope/S/PnoAwLw63fF9J1K67dzv7Hh7AZebPki1b2rk55jZvKPR2lV3ZPbT00Ev1rQdYzWMRWl3CjpVVi/OTqIUnryTpwge+vvBI96UTVLtMstyOWd3e8Q4RvB4EaDAZi2/jgp2SZ+vHUXYt9JeOQ7cPXQOKmilB/X6oWVSfFNDQW2iYQ+Dkml3JTk+QvQeXjgP2iQ1lGcxtw/53Ih+wKL7lqEi96FA+fS+Gr3WUa3NBBy4GNodB/Uu13rmIpSrlxrIuFV1/xQnI8pLo6MdesIGDwYva+v1nGcwsnUkyw+spj7695Pq5BWFFisjFt1iCperjybMxv0rnDXe1rHVJRy53qWtFXKgZTPvwC9noAhg7WO4hSs0srkXZPxcvVidKvRACzZeYYjFzKY0/wMhtObocd48AnVOKmilD+qgFQgBcnJpK1YgW/fPriEhGgdxymsiV7DH4l/8FKrl/B38yc+PY8Zv/5Fr7rutDj6PoQ2h9aqtb2i3AhVQCqQlKVLkSYTgUOf1DqKU0jJS+GDfR/QskpL7qt7HwBvrz2K2WJlqv9qRM4l6P0h6PQaJ1WU8kkVkArCkpVN6tfL8O7ZE2Od2lrHcQozomaQbcpmQvsJCCHYciKRHw9d5O3IPLwPLYE2w6BaC61jKkq5pQpIBZH27bdYMzIIfFqdjgHYG7+XNTFrGNJkCBF+EeSZLUxYc4S6QW4MuPgBeFeF7uO0jqko5ZqaolwBWE0mUhYtwqNdO9ybNdM6jubMFjOTd00mzCuMYc2GATB7czRnU3LY3PEIun2H4IHF4KauQleUm6GOQCqAjO+/pyAxUa1zbrfoyCJOpZ9ibNuxuBvciU7M4rPfYniiiQu1D34I9e6ARn21jqko5Z4qIOWctFhIXvg5xka34NlR9XA6l3mOuQfncnut2+lSvcvlZonuLnpe43OQVrh7muq0qyilQBWQci5z40ZMp08T9PTTlb5lu5SSKbunoBd6xrS29bRac+ACO08lM6vlRYzRP0HXV8E/XNugilJBqAJSjkkpSZ6/AJeaNfG+Qy0Ytf7Meraf387zLZ4nxDPkcrPENtXd6BozDYJvgQ7Pax1TUSoMVUDKsZzde8g7dIjAoUMR+so9lyHLlMX7e97nloBbeKjhQwBM/cXWLHF2tV8Q6XG2RaL0LhonVZSKQ12FVY4lz5+PPigI3/vv0zqK5j458AmXci8x67ZZGHQG9p9N5es9Z3m1eQHBhxZAi8egVnutYypKheKwIxAhxOdCiEQhxOEi2wKEEL8KIU7ab/3t24UQYpYQIloIcVAI0dJRuSqK3CNHyN6+nYDBj6MzGrWOo6kjyUdYdnwZAxsMpElQE3uzxMOEeLkwLH0WuPvB7ZO0jqkoFY4jT2EtAv69luprwEYpZT1go/0xQC+gnv1rGDDHgbkqhJSFC9F5eeH/0ENaR9GUxWrh7Z1vE+AWwAstXwBszRKPXsxgfpOj6C9EwR1TwCNA46SKUvE4rIBIKbcCKf/a3BdYbL+/GLivyPYl0mYX4CeEUO1R/4PpzBkyfv4F/0EPofeu3B33vznxDUeSjzCm9Ri8Xb2JT8/jg/Un6B1hoMmxGRDeGW6t3EVWURylrAfRQ6SUFwHst1Xs28OAc0VeF2ffdgUhxDAhRJQQIiopKcmhYZ1V8udfIAwG/B97TOsomkrMSWTW/ll0qNaBO8PvBGzNEguskne9liNMOXDPDDXnQ1EcxFmuwiruf3hxKyEipZwnpYyUUkYGBwc7OJbzKUhKIn3VKnzvuw+XKlWu/YYKbNreaZgtZsa1HYcQgs32ZonvNU/B68RK6DQagutrHVNRKqyyLiAJhaem7LeJ9u1xQI0ir6sOXCjjbOVCypKlyIICAp8cqnUUTW0/v52fY3/m6WZPU9OnJnlmCxPXHKFBkAt9L3wA/rWh80tax1SUCq2sL+P9HhgMvGe/XVNk+0ghxHKgLZBeeKpL+ZslM5PUZcvwvvMOXGvV0jqOZvIK8pi8azLhPuEMbWIrpIXNEn9vvQvdoRh4bBW4uGuc9MaYzWbi4uLIy8vTOopSzrm5uVG9enVcXBwz/8lhBUQIsQzoBgQJIeKAidgKx7dCiCeBs8AD9pevA+4GooEc4AlH5SrPUpcvx5qVVembJs4/NJ+4rDgW3rEQV73r5WaJwxpZqXH0M2gyACJu0zrmDYuLi8Pb25vw8PBK355GuXFSSpKTk4mLi6N2bcesEeSwAiKlHPQfT/Uo5rUSeM5RWSoCa34+KUuW4NmhA+6NG2sdRzOn0k7x+eHP6V2nN21C2xRplqjj/wo+BYM73PmO1jFvSl5enioeyk0TQhAYGIgjLzZylkF05RrSV6/BknSJwGFPax1FM1JK3t71Nh4GD16OfBmA1QfOs/NUMnOaxeB69nfoOQG8y/968Kp4KKXB0T9HqoCUA9JiIfnzhbg1bYpH27Zax9HMD6d+ICohitGtRhPoHkh6jpkpPx6jY5ieDtEzICwSWlXuiwsUpSypAlIOZP76K+YzZwl86qlK+5tpen460/dO59bgW+lXrx/wd7PEj4PXIHJTbc0SdepHujTo9XqaN29++eu9994rlf3efffdpKWlXfdzV7No0SJGjhx5s9GUG6CaKTo5KSXJ8+bjGh6Od88rho8qjZn7ZpJhymB8u/HohO5ys8TxzTIJOP41tB8JoWo539Li7u7OgQMHSn2/69atu2KblBIpZbHPKc5NFRAnl71jB3lHjxI6+e1K27J9f+J+VpxcwZDGQ2gQ0OBys8RQLz2DUz4CnzDo9rrWMR3irR+OcPRCRqnus1E1Hyb2vrELMcLDw3n44YfZvHkzZrOZefPm8frrrxMdHc0rr7zC8OHD2bJlCxMmTCAwMJATJ07QpUsXPv30U3Q6HeHh4URFRZGVlUWvXr3o3r07O3fuZPXq1XTt2pWoqCiCgoJYsmQJ06dPRwhBs2bNWLp0KT/88AOTJ0/GZDIRGBjIV199RUhI+R/vKs/U8b6TS56/AEOVKvj06aN1FE2YrWYm7ZxEqGcoI24dAcBie7PEhQ2i0CcdhV5TweilcdKKJTc39x+nsL755pvLz9WoUYOdO3fSuXNnhgwZwnfffceuXbuYMGHC5dfs2bOHDz74gEOHDhETE8PKlSuv+B4nTpzg8ccfZ//+/dQqMq/pyJEjTJkyhU2bNvHnn3/y0UcfAdCpUyd27drF/v37eeihh5g6daoD/waUklBHIE4s99Ahcnbtosorr6BzddU6jia+PPol0WnRzOo+Cw8XD+LT85ix/gT96lhpeGI2NLgbbrlX65gOc6NHCjfraqew+th/mWnatClZWVl4e3vj7e2Nm5vb5TGMNm3aUKdOHQAGDRrEtm3bGDBgwD/2U6tWLdq1a3fF/jdt2sSAAQMICgoCICDA1kk5Li6OgQMHcvHiRUwmk8PmNiglp45AnFjy/AXofHzwG/ig1lE0cSHrAnP+nEP3Gt3pXrM7AJPWHqHAamWycbGtgVov9VtoWTPa15/R6XSX7xc+LigoAK68fLS4iz88PT2L3b+UstjXP//884wcOZJDhw4xd+5cNVPfCagC4qTyT50m89df8X94EHqvynd6RkrJO7ttEwJfb2Mb39h8IpF1h+KZ2SwOj9O/2sY9/GpcbTeKRvbs2cPp06exWq188803dOrUqcTv7dGjB99++y3JyckApKTYVoVIT08nLMzWpHvx4sX/+X6l7KgC4qRSvvgc4epKQCVt2b7p3CZ+i/uN55o/R6hX6OVmiY2DBL3OzYSQJtBuhNYxK6x/j4G89tpr135TEe3bt+e1116jSZMm1K5dm/vvv7/E723cuDHjxo2ja9eu3Hrrrbz0kq0p5ptvvskDDzxA586dL5/eUrQlbF1EyqfIyEgZFRWldYxSZ05IILrn7fg/MICqRQYmK4tsczZ9V/fF1+jL8nuX46Jz4YP1J/h4UzQ7Wmyg2rEv4MlfoUZrraM6xLFjx7jlllu0jnHDtmzZwvTp01m7dq3WURSK/3kSQuyTUkbe7L7VILoTSlm8BKxWAoZWzlnVnx74lIScBKZ3nY6LzuVys8SRt2RT7fgiaDWkwhYPRSlPVAFxMpb0dNKWL8enVy9cq1fXOk6ZO55ynK+OfcWA+gNoXqU5UkreWH0ITxfBi7mfgkcg9JyodUzlKrp160a3bt20jqGUAVVAnEzqsuVYc3IIfOpJraOUOau08vbOt/E1+jKq5SjA1ixx16kUvm15GJej+6HfAnD31zipoiigBtGdijUvz9ayvUtn3Bo21DpOmfvur+84eOkg/xf5f/gafUnPMTN57TG6hVlpHfMJ1OkGTQdcazeKopQRVUCcSNrKlVhSUgh6uvK1bL+Ue4kP931Im6ptuLeObWLg1F+Ok5pj4iO/bxEF+XDPDKikzSQVxRmpAuIkZEEBKZ9/gfutt+IeedMXR5Q706Omk2fJ4412byCEuNwscVKTRHxjvofOL0NghNYxFUUpQhUQJ5Hx8y+Y4+IIHPZ0pWvZvvPCTn489SNPNn2S2r61LzdLrOkleDjpIwisC51GaR2zUils596kSRN69+59Q23WixMbG0uTJk1KZV+K9lQBcQJSSpIXLMA1IgKv7t21jlOm8i35TNk9hRreNXiqqW2t98JmiQsitqJLO21b58NgvMaelNJU2Avr8OHDBAQEMHv2bK0jKU5IXYXlBLK3bSP/+HFC330XUckWRPr80OecyTjD3J5zMeqNXEzPZcb6Ewyqk0PdvxZAs4egdhetY2rnp9cg/lDp7rNqU+hV8gWi2rdvz8GDBwHIysqib9++pKamYjabmTx5Mn379iU2NpZevXrRqVMnduzYQVhYGGvWrMHd3Z19+/YxdOhQPDw8/tHSJC8vjxEjRhAVFYXBYGDGjBl0796dRYsWsXr1aiwWC4cPH+bll1/GZDKxdOlSjEYj69atu9xgUdFW5fq0clLJ8+ZjqFoV33vu1jpKmTqTcYb5h+bTK7wXHcI6APD22qMUWK1MFAsRrp5wx2SNU1ZuFouFjRs3Xu7A6+bmxqpVq/jjjz/YvHkzL7/8MoXdLE6ePMlzzz3HkSNH8PPzY8WKFQA88cQTzJo1i507d/5j34VHNYcOHWLZsmUMHjz4coPEw4cP8/XXX7Nnzx7GjRuHh4cH+/fvp3379ixZsqSs/vjKNagjEI3l7N9Pzt69hLz+GqIStWyXUvL2rrdx07vxaptXgb+bJS5o9hduf+2E3h+BV7DGSTV2HUcKpamwF1ZsbCytWrXi9ttvB2z/bmPHjmXr1q3odDrOnz9PQkICALVr16Z58+YAtGrVitjYWNLT00lLS6Nr164APPbYY/z0008AbNu2jeeffx6Ahg0bUqtWLf766y8AunfvfrlNvK+vL7179wZsLeQLj4YU7WlyBCKEiBVCHBJCHBBCRNm3BQghfhVCnLTfVorZYskLFqL39cVvQOWa37Du9Dp2X9zNCy1fIMg9iDyzhQlrDtMiyEKPuI+hRlto8bjWMSutwjGQM2fOYDKZLh8tfPXVVyQlJbFv3z4OHDhASEjI5aOGoq3d9Xo9BQUF/9maHeBqffj+3Sa+aAv5wpbxiva0PIXVXUrZvEhDr9eAjVLKesBG++MKLT8mhqyNG/F/5BF0/7E2QkX027nfmLJ7Ck0Cm/BA/QcA+GRTNOdScvms6g+IvHTbwHklGw9yRr6+vsyaNYvp06djNptJT0+nSpUquLi4sHnzZs6cOXPV9/v5+eHr68u2bdsAWwEq1KVLl8uP//rrL86ePUuDBg0c94dRSp0z/Q/tCxQ2+f//9u47PKoyffj4956SRjAkhBYIVYrUhCoi0gMoVvCHiJRVRBZQl11EWdQVfyKiKLwguksREPCVorC8gIgoVQQCAaNINUBCh0BCejIzz/vHHEKAAGlTAs/nunJx5tT7PMzMPc8p95kPPOHBWNwicfYcxM+P4AHPeToUt8ix5/BR9EeM/GkkYWXC+LDDh5hNZo6cS+U/m//kH/UvUOnIEmg7Aip55kl82o0iIyNp1qwZX3/9Nf3792fXrl20bNmSRYsW0aAAFRPmzp3LiBEjaNu2Lf7+/rnjhw8fjt1up0mTJvTt25d58+Zd0/PQvJ9HyrmLyFHgEqCA/yilZopIklKqXJ55LimlbjiMJSJDgaEA1atXb3G7X0DeKuf0aY50iyK4Xz8qKL3B0wAAGpFJREFUj/unp8NxuYSUBMZsGsPvib/zTP1nGN1qNL5mX5RS9Ju1nUOnLrKz/HgstgwYsR187p4e2fVKezl3zbvcieXc2ymlTolIReAHETlQ0AWVUjOBmeB8HoirAnS1i/Ocna3ygwd5OBLXW3tsLeO3jUdEmNJxCl1rdM2dtnyPs1ji8qY7sBw6CP0W39XJQ9NKE48kEKXUKePfcyKyHGgNnBWRKkqp0yJSBTjnidjcwXbpEpeWLiXokYexGo/ovBNl2jKZFD2JZYeW0bRCUz586EOqBl7d3+T0HCas3k9UWCYRcTPhvkehfg8PRqxpWmG4/RyIiJQRkbJXhoEo4HdgJXDl5/gg4L/ujs1dLn31FSo9nZAX7tyS7X8m/Um/1f1YdmgZzzd+nnk95l2TPAAmfX+AS+lZfFxmAWIyQ49JHopW07Si8EQPpBKw3Li0zwJ8pZRaKyLRwBIReQGIB572QGwu50hP59KChQR26oRfvXqeDqfEKaVYcWQFE3dOxN/iz+ddP+fBqg/eMF9M/CX+7854PrzvGGXjNkD3iRB05/bGNO1O5PYEopSKA5rlMz4R6OLueNwt6ZtvsSclUf7FIZ4OpcSl5aTxv9v/l9Vxq2lduTUT20+kYkDFG+bbdewiY5bFUivQQe9z06FyU2g91AMRa5pWHPpOdDdSOTkkzv0C/xYtCGje3NPhlKj9ift5bfNrJKQkMCJiBC82eRGzyXzNPIfPpjBp7UHW7z9LxbK+rKj9A6aDZ6HfV2DWb0VNK2286T6QO97lNWuwnTp9R/U+lFIs2r+I/mv6k5GTwZyoOQxrNuya5HE6OYPXl8XSfepmtsclMjqqHpueK0fYoQXQaghUbeHBPdDyUxrLuX/55Zc0btyYRo0a0bBhQyZPnlzkdfXr14+mTZsyZcoU3n77bdavX1+CkeZv6tSppKen575++OGHS6zdXUX/7HMT5XCQOHs2vnXrEmjUBSrtkrOSefvnt/kp4SceqvYQ77V7j2C/q7fuJGfk8O9Nf/LF1qM4lGLwA7UY2fleQhJjYOUrUKYidHnLg3ug3cyVUiYAgwYNYsaMGYwbN87DUd3cd999x9SpU1m3bh1hYWFkZmayYMGCIq3rzJkzbNu27bZ32ReWUgqlFKabVFiYOnUqzz33HAEBAQCsWbOmRLfvCjqBuEnqpk1kHT5C2IeT7ogHRu09t5cxm8dwPuM8o1uOZmDDgbn7lZljZ8Evx/l0wxGSM3J4IiKMf0TVJzzzICx/Fo784EweT/0H/II8vCfebdLOSRy4WODbpAqkQUgDXm/9eoHnLw3l3CdOnMjkyZMJCwsDnFWDXzQeDb13716GDRtGeno6derU4YsvviA4OJiOHTvSpk0bNmzYQFJSEnPmzKF9+/ZERUVx7tw5IiIimD59OnPmzKFXr1706dOHNWvW8Pe//53Q0FCaN29OXFwcq1at4p133iEwMJDRo0cD0LhxY1atWgVAz5496dSpE7/88gsrVqzggw8+IDo6moyMDPr06cP48eOZNm0ap06dolOnToSGhrJhwwZq1qzJrl27CA0N5ZNPPuGLL74AYMiQIfztb3+7ZZu7iz6E5SaJs2ZjDQvjnp49PR1KsTiUg9m/zWbw2sGYxMSCngsY1GgQIoLdofhm9wm6fLyJCWv20yy8HKtefpCpnXwJX/cizOwIJ3dBt3fh1V+hTmdP7452G6WlnPvvv/9Oixb5HwodOHAgkyZNIjY2liZNmjB+/PjcaTabjZ07dzJ16tTc8StXrqROnTrs3buX9u3b586bmZnJSy+9xHfffcfWrVs5f/58gdrw4MGDDBw4kD179lCjRg0mTJjArl27iI2NZdOmTcTGxvLKK68QFhbGhg0b2LBhwzXL7969m7lz57Jjxw62b9/OrFmz2LNnzy3b3F10D8QN0nfvJiMmhkpvvolYrZ4Op8guZFxg3NZxbDu1je41u/Ovtv+irE9ZlFJsPHSeSd8d4MCZFBpXvYcP+zSlXXAybPwH/LYMfMtCx3/C/X8Fv3s8vSulRmF6CiXpTinnfv32Bw0axNNPX71D4Kmnnrom3ls5cOAAtWvXplatWoDzPMnMmTNvG0ONGjW4//77c18vWbKEmTNnYrPZOH36NH/88QdNmza96fJbt27lySefpIxRcPWpp55iy5YtPPbYY/m2uTvpBOIGibNmYw4OplzvpzwdSpFtP72dsVvGkpKdwlv3v8XT9Z5GRPg1IYmJ3+1ne9xFqocEMK1fJL3CczBteQf2fuV8FO2Do+CBlyFAP0WutLhyDiQ5OZlevXoxY8YMXnnllWvKuVutVmrWrHnTcu4ZGRluK+feqFEjdu/eTefOhevVXlnvlfLzt3KreC0WCw6HI/f1lTYBcr/4AY4ePcrkyZOJjo4mODiYwYMHXzNvYbebX5u7kz6E5WKZBw+RunEjwQOew+TGY5MlxeawMX3PdIauG0pZn7J89chX/E/9/+FYYjojFsXw+IyfOXw2lfGPNWL9i/V47MQnmD5tAbFLoM1LzkNVXf+lk0cpVVrKuY8dO5YxY8Zw5swZALKyspg2bRpBQUEEBwezZcsWABYsWJDbGymsBg0aEBcXl/srf/HixbnTatasSUxMDAAxMTEcPXo033VcvnyZMmXKEBQUxNmzZ3N7YwBly5YlJSXlhmUeeughVqxYQXp6OmlpaSxfvvyaQ2uepHsgLpY4ZzYSEEDIs896OpRCO5N2htc3v07MuRieuPcJxrYeS1qmmTdX/MbXOxPwsZh4pUtdhra8h8DoT2HGLHDYIHIAPPSavrP8DnF9OfdHH32Uli1bEhERUeBy7ldOonfv3j13/PDhwxk2bBhNmjTBYrEUq5z7ww8/zNmzZ+natWtur+f5558HYP78+bkn0WvXrs3cuXOLtA1/f38+++wzevToQWhoKK1bt86d1rt3b7788ksiIiJo1aoV9W5SZaJZs2ZERkbSqFEjateuTbt27XKnDR06lJ49e1KlSpVrzoM0b96cwYMH525vyJAhREZGuv1wVX48Us69pLRs2VLt2rXL02HcVM7JkxyJ6k7IgAFUesMzx7KLamPCRt78+U1y7Dm8ef+bdKrWk5mb45i9JY4sm4N+rcN5tV1FKvw2C7Z/Djnp0PQZ6DAGQmp5OvxSTZdz916pqakEBgailGLEiBHUrVuXUaNGeTqsW7oTy7nfFRLnzgOTiZBSVLI9x57DJ7s/YeH+hTQIacD77Sax7YCJDl9tIDEtm0eaVGF0x6rU+nMBzJkGmcnQ6EnnCfIKd15tL03La9asWcyfP5/s7GwiIyN56aWXPB2SR+kE4iK2ixdJWraMoEcfxVq5sqfDKZCEywm8tvk19iXu45n6/Wjg+ywvzD5K/MV02tQK4YuoWjQ7vQwWTYH0C1CvJ3QeB5WbeDp0TXOLUaNGeX2Pw510AnGRSwsXobKyKP/C854OpUDWHl3LO7+8g0lMvNTgXb7fWYFZJ/fRoHJZ5g1sRofUtci3f4GU01C7E3R+E6oVuwesaVopphOICzjS0ri4aBGBXTrjW6eOp8O5pQxbBpN2TuKbw99QN6gxPhcHMHm5omq5bD7p3YgnzFswrRsBSfFQvS30ng01byzPrmna3UcnEBe4tHQpjuRkQod4d9HEP5P+ZPSm0RxJOkINcy9itrclyN/CuJ61GRQUg8+W1yDxCIRFQq8pUKcL3AFlWDRNKxk6gZQwlZ3NxXnzCWjdGn/jDlFvc+WhTxN2TACHL9knXuBwRn3+2qEmI6seoszWgXBuH1RsCH0XQYNHdOLQNO0G+kbCEpa8ajW2M2cobxRy8zZpOWm8tul13t72Nlmp4Vw6NJIn6nVi29OK1xP+SplvB4I9C3rPgWE/w329dPK4C02YMIFGjRrRtGlTIiIi2LFjR6HXUZTS7VfKyF/584Z7HbSb0z2QEpRbsv2++yjzYLvbL+Bmv53bx4j1f+dS9mmyznejfcVneKd3MuF73oDl2yCoOjw+w3k/h37A013rl19+YdWqVcTExODr68uFCxfIzs52y7bzlpHPj81mw2LR701vof8nisCemkZOQjzZx+PJTognJz6e7PgEso8fx3bmDGEfT/aqku0Oh4O3N/6H/8bPxGELoLr9H3zcsRoN978P/90AZavAIx9D5ECw+Hg6XC2PM++/T9b+ki3n7ntfAyr/8583nX769GlCQ0Nz7woPDQ3NnRYdHc2rr75KWloavr6+/PjjjyQmJjJgwADS0tIA+PTTT3nggQeuWafdbueNN95g48aNZGVlMWLEiALfQzFv3jxWr15NZmYmaWlp/PTTT3z00UcsWbKErKwsnnzyydxKugsXLmTatGlkZ2fTpk0bPvvsM8xm8222oBWVTiD5UEphT0oyEoPzLzdJxMdjT0y8Zn5z+fL4hIdTpk1r/Jo05Z4ePdwar8OhSM7I4WTyZeIuniHh8jlOXj7H2fTzJGYkcjJjP5nWP/DJaczE+v3pdmIBsnotBJSHqAnQ6gWwlr46XZprREVF8e6771KvXj26du1K37596dChA9nZ2fTt25fFixfTqlUrLl++jL+/PxUrVuSHH37Az8+Pw4cP069fP66vEDFnzhyCgoKIjo4mKyuLdu3aERUVlVvZ9oorVYDBWd13+fLlgLNXFBsbS0hICOvWrePw4cPs3LkTpRSPPfYYmzdvpkKFCixevJiff/4Zq9XK8OHDWbRoEQMHDnRPw92F7toEohwObOfPk338ODkJCVd7E8fjyU5IwHFdUTNLlSr4hIdTtnMnrOHV8aleHZ8a1bGGh2MODCz5+JTicqaNU0ZSiE8+y8mUc5xNvUBiZiLJ2RdJtSWR6UjCJpcRcypizr+qp8lShkfKPcmEzOOYf+oPvkHO+zjaDHOWWde81q16Cq4SGBjI7t272bJlCxs2bKBv37588MEHtGjRgipVqtCqVSsA7rnHWZY/LS2NkSNHsnfvXsxmc25J9rzWrVtHbGwsy5YtA5xl1g8fPnxDArnZIaxu3brlPkRq3bp1rFu3jsjISMBZXuTw4cPExsaye/fu3PgyMjKoWLFiCbWKlh+vSyAi0gP4P4AZmK2U+qCo61I2GzmnT5N9PD7PIacEcuKPkx2fgMrKujqzxYK1ahg+4dUJiogwkoORJKpVw1TEIm/XS8uycTLpMnGXznA86QwnU85zNs3ZU0i6khTsyeRIMmJOQcxZ+a7HonzxN5ehsiWAcuYQypurUtHiR2WLL1WsvlSz+lDFJJRX4Jt8CvZOB4s/tB8ND4wE/+B816tp4DyZ3bFjRzp27EiTJk2YP38+zZs3z/fQ7JQpU6hUqRK//vorDocDPz+/G+ZRSjF9+vRriikWRt6S6Eopxo4de8MhsOnTpzNo0CAmTpxYpG1ohedVCUREzMAMoBtwAogWkZVKqT9utowjK8vZg4hPIDv+ODnGYabshHhyTp6CPDX+xdcXa3g41urV8XugHeZq4UhYZQirgj0kGLso7HYbqfZsbLYc7A4bNlsStrgL2Ow27HbnOLvDRo7dht2Wg8N4bXfYsTlyyLJlcT7tPBcyErmUncTlnMukqFTSyCBdMskwZ5NjduS7L2XsEKygPIryNgcVsmxUsNsIzcki1G6jvN1BqN1OiN1Ogc9UWPzBNxDuH+58LkeZ0Nsvo93VDh48iMlkom7duoDzkbA1atSgQYMGnDp1iujoaFq1akVKSgr+/v4kJydTrVo1TCYT8+fPx26337DO7t278/nnn9O5c2esViuHDh2iatWq1ySGgurevTtvvfUW/fv3JzAwkJMnT2K1WunSpQuPP/44o0aNomLFily8eJGUlBRq1KhR7DbR8udVCQRoDRxRSsUBiMjXwONAvgkkdf8+9jeLuOZa5AwfOB8M58rBuZZwtpxwNhjOBMOlwGzschSHHMUBOBziTFMnXLQ3AveY7YTaHYTb7ZS32wmxOwhRJkKUmRCTlVCzHxUt/pS3BOBjDXB+4Vv9wBoAFj/nuQmrvzEc4Jxm8b863up/62W86GS+Vjqkpqby8ssvk5SUhMVi4d5772XmzJn4+PiwePFiXn75ZTIyMvD392f9+vUMHz6c3r17s3TpUjp16pRvUhgyZAjHjh2jefPmKKWoUKECK1asKFJ8UVFR7N+/n7Zt2wLOQ24LFy6kYcOGvPfee0RFReFwOLBarcyYMUMnEBfyqnLuItIH6KGUGmK8HgC0UUqNzG/+OsEBauwT9UkuZ+JyOQvJ5Sxk+ZkwmUyYMCEIIsawiPGvCbMIghmTCCYxO8cZ00y5f8Z4MTtfm8zOYZNzmslkxiIWxGTCLBZMJhNmswUfkw9h91SgerkwKpWtjNU3MM8XewCYrfpLXbslXc5dK0l3Uzn3/L5Zr8lwIjIUGApQvXp1hszd4464NE3TtOt4253oJ4DwPK+rAafyzqCUmqmUaqmUalmhQgW3BqdpmqZd5W0JJBqoKyK1RMQHeAZY6eGYNM3tvOnQslZ6ufp95FUJRCllA0YC3wP7gSVKqX2ejUrT3MvPz4/ExESdRLRiUUqRmJiY72XVJcXbzoGglFoDrPF0HJrmKdWqVePEiROcP3/e06FopZyfnx/VqlVz2fq9LoFo2t3OarXecIe2pnkjrzqEpWmappUeOoFomqZpRaITiKZpmlYkXnUnemGJSApw0NNxeIlQ4IKng/ASui2u0m1xlW6Lq+orpYpdiru0n0Q/WBK3498JRGSXbgsn3RZX6ba4SrfFVSKy6/Zz3Z4+hKVpmqYViU4gmqZpWpGU9gQy09MBeBHdFlfptrhKt8VVui2uKpG2KNUn0TVN0zTPKe09EE3TNM1DdALRNE3TisRrE4iI9BCRgyJyRETeyGe6r4gsNqbvEJGaeaaNNcYfFJHu7ozbFYraFiLSTUR2i8hvxr+d3R17SSvO+8KYXl1EUkVktLtidpVifkaaisgvIrLPeH+4rmSrGxTjM2IVkflGG+wXkbHujr2kFaAtHhKRGBGxGU+BzTttkIgcNv4G3XZjSimv+wPMwJ9AbcAH+BVoeN08w4F/G8PPAIuN4YbG/L5ALWM9Zk/vk4faIhIIM4YbAyc9vT+eaos8078BlgKjPb0/HnxfWIBYoJnxuvxd/Bl5FvjaGA4AjgE1Pb1PLm6LmkBT4EugT57xIUCc8W+wMRx8q+15aw+kNXBEKRWnlMoGvgYev26ex4H5xvAyoIuIiDH+a6VUllLqKHDEWF9pVeS2UErtUUpdeaLjPsBPRHzdErVrFOd9gYg8gfNDcSc8Y6Y4bREFxCqlfgVQSiUqpexuitsVitMWCigjIhbAH8gGLrsnbJe4bVsopY4ppWIBx3XLdgd+UEpdVEpdAn4AetxqY96aQKoCCXlenzDG5TuPcj6IKhnnL6mCLFuaFKct8uoN7FFKZbkoTncocluISBngdWC8G+J0h+K8L+oBSkS+Nw5ljHFDvK5UnLZYBqQBp4F4YLJS6qKrA3ah4nz/FXpZby1lIvmMu/5645vNU5BlS5PitIVzokgjYBLOX56lWXHaYjwwRSmVanRISrvitIUFeBBoBaQDP4rIbqXUjyUbotsUpy1aA3YgDOdhmy0isl4pFVeyIbpNcb7/Cr2st/ZATgDheV5XA07dbB6j+xkEXCzgsqVJcdoCEakGLAcGKqX+dHm0rlWctmgDfCgix4C/Af8UkZGuDtiFivsZ2aSUuqCUSsf5BNDmLo/YdYrTFs8Ca5VSOUqpc8DPQGmul1Wc779CL+utCSQaqCsitUTEB+dJr5XXzbMSuHKVQB/gJ+U8E7QSeMa46qIWUBfY6aa4XaHIbSEi5YDVwFil1M9ui9h1itwWSqn2SqmaSqmawFTgfaXUp+4K3AWK8xn5HmgqIgHGl2kH4A83xe0KxWmLeKCzOJUB7gcOuCluVyhIW9zM90CUiASLSDDOIxbf33IJT181cIurCR4GDuG8omCcMe5d4DFj2A/n1TRHcCaI2nmWHWcsdxDo6el98VRbAG/iPL67N89fRU/vj6feF3nW8Q6l/Cqs4rYF8BzOiwl+Bz709L54qi2AQGP8PpxJ9DVP74sb2qIVzt5GGpAI7Muz7PNGGx0B/nK7belSJpqmaVqReOshLE3TNM3L6QSiaZqmFYlOIJqmaVqR6ASiaZqmFYlOIJqmaVqR6ASi3XVExC4ie0XkdxFZKiIBhVw+tZDzz7u+6qkxvqWITDOGB4vIp8bwMBEZmGd8WGG2p2nuohOIdjfKUEpFKKUa4yyeNyzvROOmMpd/NpRSu5RSr+Qz/t9KqS+Nl4NxltnQNK+jE4h2t9sC3CsiNY3nQXwGxADhItLPeE7E7yIyKe9CIvKxUYjwRxGpYIx7UUSiReRXEfnmup5NVxHZIiKHRKSXMX9HEVl1fUAi8o6IjDZ6LS2BRUaP6RERWZ5nvm4i8m3JN4mmFYxOINpdyyjj0RP4zRhVH/hSKRUJ5OAsQNkZiABaGeXgAcoAMUqp5sAm4F/G+G+VUq2UUs2A/cALeTZXE2fJkEeAf0sBHuCklFoG7AL6K6UicNasuu9KwgL+Aswt9I5rWgnRCUS7G/mLyF6cX87xwBxj/HGl1HZjuBWwUSl1XjnLfy8CHjKmOYDFxvBCnJVtARobvYzfgP5AozzbXKKUciilDuN8JkmDwgatnGUjFgDPGXXO2gLfFXY9mlZSvLWcu6a5Uobxiz6XUeI9Le+oQqzvSj2gecATSqlfRWQw0DGfeW72uqDmAv8PyASWGslN0zxC90A0LX87gA4iEioiZqAfzsNV4PzcXLmq6llgqzFcFjgtIlacPZC8nhYRk4jUwfm40YMFjCPFWC8AyvmEyVM4C2XOK9QeaVoJ0z0QTcuHUuq0iIwFNuDsjaxRSv3XmJwGNBKR3TifbNfXGP8WzsRzHOd5lbJ5VnkQZwKqBAxTSmUW8MFW83CeM8kA2iqlMnAeTquglCrNJdi1O4CuxqtppYxxv8gepdSc286saS6kE4imlSJGrycN6KZK9/PttTuATiCapmlakeiT6JqmaVqR6ASiaZqmFYlOIJqmaVqR6ASiaZqmFYlOIJqmaVqR/H8Dv43OBBCr3QAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.plot(p_list, nb_infected, label='Empirical')\n",
"plt.plot(p_list, nb_infected_random, label='Random')\n",
"plt.plot(p_list, nb_infected_random_config, label='Random Configuration')\n",
"plt.plot(p_list, nb_infected_scale_free, label='Scale Free')\n",
"plt.xlim(0, 0.1)\n",
"plt.xlabel('Probability')\n",
"plt.ylabel('Infected')\n",
"plt.legend(loc='best')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The most clearest difference is that the scale-free network is much more \"infectious\", just because of its structure. This is a well-known result: ordinary random graphs have some threshold (here around 0.03) below which almost nobody is infected, and above which a sizeable part of the network becomes infected. In the scale-free graph this is radically different, and even for the smallest probabilities, infections always spread to a sizeable part of the network. The difference around 0.03 is substantial: more than 200 infections in a scale free graph, while not even 50 infections in the other graphs.\n",
"\n",
"Nonetheless, for the highschool social network, the random graph is a very reasonable approximation, and the configuration model is an even better approximation. This is important because for further analysis we could use the random approximation, which makes it significantly easier and simpler to get further insights.\n",
"\n",
"Remarkably, the presence of a clear group structure in terms of classes does not seem to have any effect. The group structure is completely absent in the configuration model, yet they show almost identical results. This is quite different compared to the opinion dynamics, which depend clearly on the group structure.\n",
"\n",
"There are many other questions that are of interest here: how does the spreading depend on which node is infected first? How is it best contained? How many nodes need to be vaccinated so that spreading is unlikely? We won't go into these question here, but as you can imagine, they are of great importance."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Community detection"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A common task in network analysis is *community detection*: finding groups of nodes that are relatively densely connected within communities, but sparsely between communities. It provides a sort of birds-eye view of the network, which helps in the interpretation of the network. Communities often seem to share some characteristic (or in for example biological networks, some function).\n",
"\n",
"The method that has been most prominent in the literature is called modularity. It is based on comparing the real network to a random network using the configuration null-model. The idea is that we want to find as many edges within communities compared to the expected number of such edges. Let us illustrate this on the highschool network. We use the recently introduced Leiden algorithm to detect communities."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Modularity"
]
},
{
"cell_type": "code",
"execution_count": 79,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 79,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"partition = leidenalg.find_partition(G, leidenalg.ModularityVertexPartition, weights='weight')\n",
"ig.plot(partition)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As you can see, modularity finds a number of groups that seems to uncover the classes to some extent. The actual overlap between the two partitions (one partition based on the classes and the other partition based on modularity) can be calculated using something called the *normalised mutual information* (NMI). This value reaches 1 if the two correspond exactly with each other, and reaches 0 if there is no correspondence."
]
},
{
"cell_type": "code",
"execution_count": 80,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.9350187443979149"
]
},
"execution_count": 80,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"G.vs['modularity'] = partition.membership\n",
"class_partition = ig.VertexClustering.FromAttribute(G, 'class')\n",
"partition.compare_to(class_partition, 'nmi')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### CPM"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Indeed modularity thus seems to uncover the classes quite well. Nonetheless, some classes seem to be split in multiple communities, probably because some social groups within the classes exist. Let us examine one particular class from closer by."
]
},
{
"cell_type": "code",
"execution_count": 81,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 81,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"H_class = class_partition.subgraph(4)\n",
"H_class['layout'] = H_class.layout_fruchterman_reingold(weights='weight')\n",
"H_class.es['width'] = 0.01*np.array(H_class.es['weight'])\n",
"partition_of_class = ig.VertexClustering.FromAttribute(H_class, 'modularity')\n",
"ig.plot(partition_of_class)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see that only a few people are not in the same group as the rest of the class. They are also clearly not so much connected to the rest of the class. However, there does seem to be some more particular structure within this group. Let us try modularity again."
]
},
{
"cell_type": "code",
"execution_count": 82,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 82,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"partition_of_class2 = leidenalg.find_partition(H_class, leidenalg.ModularityVertexPartition, weights='weight')\n",
"ig.plot(partition_of_class2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The results are interesting: it suggest there are indeed some smaller groups within the classes. But the definition of 'community' is now no longer clear. What would be called a group at one level (i.e. the class at the overall level of the school) would no longer be considered a group at a lower level, as modularity tends to split it up. This is related to the problem of the *resolution limit* in community detection. Modularity may 'hide' smaller communities, hidden away in the communities it detects.\n",
"\n",
"Both the groups within a class, as well as taking classes as groups themselves would be valid social groups. Clearly they have a relevance in some way. Indeed, what resolution is \"better\" is not something that has a single unique answer. However, the downside of modularity is that the resolution in a sense depends on the size of the network we are looking at.\n",
"\n",
"We will now investigate another method which does not suffer from this problem, which is called the *Constant Potts Model* (CPM). Although it does require us to choose a resolution parameter, it has a definition of communities that is independent of the size of the network."
]
},
{
"cell_type": "code",
"execution_count": 83,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 83,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"resolution_parameter = 1.15\n",
"partition = leidenalg.find_partition(G, leidenalg.CPMVertexPartition, weights='weight', resolution_parameter=resolution_parameter)\n",
"ig.plot(partition)"
]
},
{
"cell_type": "code",
"execution_count": 84,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 84,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"G.vs['CPM'] = partition.membership\n",
"H_class = ig.VertexClustering.FromAttribute(G, 'class').subgraph(4)\n",
"H_class['layout'] = H_class.layout_fruchterman_reingold(weights='weight')\n",
"H_class.es['width'] = 0.01*np.array(H_class.es['weight'])\n",
"partition_of_class = ig.VertexClustering.FromAttribute(H_class, 'CPM')\n",
"ig.plot(partition_of_class)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"CPM finds that this class is completely part of a single community. If we rerun CPM on this particular network, we should find a highly similar community structure."
]
},
{
"cell_type": "code",
"execution_count": 85,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 85,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"partition_of_class2 = leidenalg.find_partition(H_class, leidenalg.CPMVertexPartition, weights='weight', \n",
" resolution_parameter=resolution_parameter)\n",
"ig.plot(partition_of_class2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we could change the resolution parameter if we want to find smaller subgroups within the classes."
]
},
{
"cell_type": "code",
"execution_count": 86,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 86,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"resolution_parameter = 25\n",
"partition_of_class2 = leidenalg.find_partition(H_class, leidenalg.CPMVertexPartition, weights='weight', \n",
" resolution_parameter=resolution_parameter)\n",
"ig.plot(partition_of_class2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If we apply this new resolution parameter also at the higher level, we should again get similar results."
]
},
{
"cell_type": "code",
"execution_count": 87,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 87,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"partition = leidenalg.find_partition(G, leidenalg.CPMVertexPartition, weights='weight', resolution_parameter=resolution_parameter)\n",
"ig.plot(partition)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Indeed, now each class seems to be partitioned in multiple smaller groups. In other words, the definition of what is a 'community' is independ of whether we look at the level of a single class, or at the level of an entire school.\n",
"\n",
"However, it can be tricky to try to find resolutions that are of interest. It seems impossible to do this automatically without creating problems, so it really depends on your own interpretation of the network and of the results. Nonetheless, there are some tools to help you. One is to construct a so-called *resolution profile*. We simply detect for each resolution within a certain range the partition. This can actually be done reasonably efficiently using bisectioning, which explains also the name of the procedure:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
\n",
"Constructing a resolution profile may take some time.\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": 88,
"metadata": {},
"outputs": [],
"source": [
"opt = leidenalg.Optimiser()\n",
"resolution_profile = opt.resolution_profile(G, leidenalg.CPMVertexPartition, resolution_range=(0, 25), weights='weight')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us plot the results, and see how they correspond to to the classes as they are present in the data."
]
},
{
"cell_type": "code",
"execution_count": 89,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAD8CAYAAAAxDXBhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3de3hU1b3/8feXBEUDooBUKv6EHtAWFSNG1OOFqFXRHsEiKmB/olLDqbeq9Xh56k8R9PFSj1KeekuFghcERE/hVFoOXkDleEs1KpciKcUaoVIuWoQiBL6/P/YemExmMslkkj3JfF7PMw8za69Ze2c2e76zLnstc3dERERaWruoD0BERPKTApCIiERCAUhERCKhACQiIpFQABIRkUgoAImISCQK02UwsynAvwHr3P3IMO1o4HGgI7AauMTd/xFuuw0YA+wErnP3+WH6YOCXQAHwpLvfF6b3BmYAXYD3gf/r7tvNbG/gKeBYYANwsbuvTne87dq183322aehf7+IiABbt251d2/RSomluw/IzE4FvgaeigtA7wE3ufsiM7sC6O3u/8/M+gHPAQOBbwMvA4eFRX0CnAlUA+8BI919mZnNAl509xlm9jjwobs/ZmZXAf3d/d/NbATwQ3e/ON0fVFRU5Fu2bGn0ByEiks/MbKu7F7XkPtNGO3d/HdiYkHw48Hr4fAFwQfh8KDDD3b9x978AVQTBaCBQ5e6r3H07QY1nqJkZcDowO3z/NOD8uLKmhc9nA2eE+UVEpA3ItLq1BBgSPr8QOCR8fjDwWVy+6jAtVXpX4Et3r0lIr1VWuP2rML+IiLQBmQagK4CrzeyPQCdge5ierIbiGaTXV1YdZlZmZhVmVlFTU5Msi4iI5Ji0gxCScfc/AWcBmNlhwA/CTdXsqQ0B9ATWhM+Tpa8H9jezwrCWE58/Vla1mRUCnanbFBg7nnKgHII+oMTtO3bsoLq6mm3btjXyL5WodejQgZ49e9K+ffuoD0UyoGsv9+TSNZVRADKz7u6+zszaAbcTjIgDmAtMN7OHCAYh9AXeJajN9A1HvH0OjABGubub2WvAcIJ+odHAnLiyRgNvhdtf9QxnTq2urqZTp0706tULdSO1Hu7Ohg0bqK6upnfv3lEfjmRA115uybVrKm0TnJk9RxAEDjezajMbA4w0s0+APxHUWH4D4O5LgVnAMuAPwNXuvjOs3VwDzAeWA7PCvAC3ADeaWRVBH8/kMH0y0DVMvxG4NdM/ctu2bXTt2lUXQCtjZnTt2lW/nlsxXXu5JdeuqbQ1IHcfmWLTL1Pkvwe4J0n6PGBekvRVBKPkEtO3EQxwaJRt26C0tHbanXdCu3bpL4AuXeDAAxu7R2lO+uJqBuXlMH167bRRo6CsrFl2p3OYW3LpfGgmhNDmzfDpp7BiRcMef/9748rv2LFj2jwTJ05k69atGf4FDbd69WqOPPLIrJZ57rnn8uWXX9abp7S0lIqKijrplZWVzJtX57eJNJfp06Gycs/rysq6AUmkBbS5ANShAyxcWPtx0EFw+OH1Pw49FDp1atg+tm6FjUmHQzRNJgFo586d2T+QDMybN4/9998/o/cqAEWguHjPBVJcHASh0tL6H+XlUR1tk5gZP/vZz3a/fvDBBxk3bhwA48aNw8yoqqravf3hhx/GzHb/WOrVqxfr16+vU667c/rpp/OPf/wj5b7HjRvHgw8+mKW/pHFuuukmXn311Uj23VBtLgBl6sAD0wep2GPffTPfz8KFCyktLWX48OF897vf5ZJLLsHdmTRpEmvWrOG0007jtNNOA+B//ud/OPHEExkwYAAXXnghX3/9NRBcEOPHj+fkk0/m+eefp7S0lFtuuYWBAwdy2GGH8cYbbwBBTeeUU05hwIABDBgwgP/93/+t99iuuuoq5s6dC8APf/hDrrjiCgAmT57M7bffDsAzzzzDwIEDKS4uZuzYsbsDYPxFOmHCBL773e9y5plnMnLkyFoX4PPPP1/rOLdv384dd9zBzJkzKS4uZubMmZl/uJKZUaOCIFSfRYtg7NhWGZD23ntvXnzxxaRBBOCoo45ixowZu1/Pnj2bfv36pS133rx5HH300ey3335ZO9Zsuvbaa7nvvvuiPox6ZTQKrjW7/vrarQ+ZiFVSYoGouBgmTmz4+z/44AOWLl3Kt7/9bU466SQWL17Mddddx0MPPcRrr71Gt27dWL9+PXfffTcvv/wyRUVF3H///Tz00EPccccdQDCU8s033wTg8ccfp6amhnfffZd58+Zx11138fLLL9O9e3cWLFhAhw4dWLlyJSNHjkzaBBZz6qmn8sYbbzBkyBA+//xz1q5dC8Cbb77JiBEjWL58OTNnzmTx4sW0b9+eq666imeffZZLL710dxkVFRW88MILfPDBB9TU1DBgwACOPfbY3duTHef48eOpqKjgV7/6VcM/RMmesrL0/T+J/UaLFgWP+LR0/UjZuPgSNeDiKywspKysjIcffph77qnTPc3555/PnDlzuP3221m1ahWdO3du0BDlZ599lrK4v/epp57iwQcfxMzo378/Tz/9dK38v/71rykvL2f79u306dOHp59+mn333Zfnn3+eu+66i4KCAjp37szrr7/O0qVLufzyy9m+fTu7du3ihRdeoG/fvjzzzDNMmjSJ7du3c/zxx/Poo48CMGbMGCoqKjAzrrjiCm644QYOPfRQNmzYwN/+9jcOOuighnyaLS7vAlC27Nq1JxBt2hT0C9XHPcjz17/CkUcOZMuWnqxcCb16FfPWW6s58MCT2bEDqqpgwwZ47bW3+fjjZZSUnATAjh3bKS4+kRUrYMcOKCm5ePc+t26FkpJhrFgBnTsfy8qVq1mxAjZv3sGECdewfHklBQUFrF79CStWQHU1JGu5O+WUU5g4cSLLli2jX79+bNq0ibVr1/LWW28xadIkpk2bxh//+EeOO+44AP75z3/SvXv3WmW8+eabDB06lNiEsOedd16t7cOGDQPg2GOPZfXq1Q3/wCW5ZAMKUmnKQIPEIJUsIMXy5aCrr76a/v37c/PNN9fZtt9++3HIIYewZMkS5syZw8UXX8xvfvObtGUuXryYJ554AoClS5dyzz33sHjxYrp168bGJG30w4YN48orrwTg9ttvZ/LkyVx77bWMHz+e+fPnc/DBB+/uR3388cf56U9/yiWXXML27dvZuXNnyh+ARxxxBJ9//jlLliwBqNUXO2DAABYvXswFF1xQ53hyQd4FoMbUVFL5+9+b1gfUvv3eu5+3a1dAstkb3J1//dczeeih55KWse++tecMjJVZULCnvKlTH6Zr128xZ86H7Nq1i6OP7gAEIwWTBaCDDz6YTZs28Yc//IFTTz2VjRs3MmvWLDp27EinTp1wd0aPHs29996b8m9Ld6vW3nvXPU5pgtiAgoY0ocVqLA3Jn05iQCot3dOPFO/OO6Fbt6CNOxsXX4b2228/Lr30UiZNmkSy2fJHjBjBjBkzmD9/Pq+88kqDAtDGjRvpFHYcv/rqqwwfPpxu3boB0KVLlzr5lyxZwu23386XX37J119/zdlnnw3ASSedxGWXXcZFF120+wfaiSeeyD333EN1dTXDhg2jb9++vPLKK0l/AJ533nmsWrWKa6+9lh/84AecddZZu/fZvXt31qxZU+dYckXeBaBsOPDAxg/XNgv6j9auhY4dg+cABxwAPXoEr7t06cS3vrWZ3r270aXLCdx779UUFFTRp08ftm7dSnV1NYcddhjt20OfPsF1DUFTYK9eQRnr10P79sHz9u2/4ogjevK977XjN7+Zxs6dOzn88KAGlMqJJ57IxIkTefXVV9mwYQPDhw9n+PDhAJxxxhkMHTqUG264ge7du7Nx40Y2b97MoYceuvv9J598MmPHjuW2226jpqaGl156afevvlQ6derE5s2bG/eByh6xAQX1ia+xFBcHtaFsSlXeN98Ew0tT/WJrwXsfrr/+egYMGMDll19eZ9t5553Hf/zHf1BSUtLgPp3CwkJ27dpFu3btcPe0w5svu+wyfvvb33L00UczdepUFobn7PHHH+edd97hpZdeori4mMrKSkaNGsXxxx/PSy+9xNlnn82TTz5Z7w/ADz/8kPnz5/PII48wa9YspkyZAgT3YeXy8jQahJBDysrKOOecczjttNM48MADmTp1KiNHjqR///6ccMIJ/OlPf2pUeVdddRXTpk3jhBNO4JNPPqGoKP1M66eccgo1NTX06dOHAQMGsHHjRk455RQA+vXrx913381ZZ51F//79OfPMM3f3E8Ucd9xxDBkyhKOPPpphw4ZRUlJC586d693naaedxrJlyzQIoTmVldUeGprtprLE8mOPLl1SDy9truGkKXTp0oWLLrqIyZMn19m2zz77cP/99/Pzn/+8weUdfvjhrFq1Cgh+nM2aNYsNGzYAJG2C27x5Mz169GDHjh08++yzu9P//Oc/c/zxxzN+/Hi6devGZ599xqpVq/jOd77Dddddx5AhQ/joo48444wzmD17NuvWrdu9j08//ZT169eza9cuLrjgAiZMmMD777+/u+xPPvkk67dcZFPa9YBam2TrAS1fvpzvfe97ER1R7lmxIrj2mzKarz5btnxNUVFH/vnPrfzoR6cyfnw5RxwxIOPy1qxZzl131T1/zXjvZOsRa/JKVwOKSL3XXqwTM9Yc0Ew6duy4ewTpF198Qe/evbn55psZN24c48aNo2PHjtx000213lNaWsqDDz5ISUkJvXr1oqKiYnfzWsyECRPo0aMHP/7xjwGYNm0av/jFLygoKOCYY45h6tSptcp/7LHHeOCBBzj00EM56qij2Lx5M1OnTmXYsGGsXLkSd+eMM85g4sSJ3HfffTzzzDO0b9+egw46iOnTp9OlSxdmzpzJvffey65du2jfvj2PPPII++yzD5dffjm7du0C4N577+Wcc85hx44d9O/fn48//pjCwtqNXcnOSxTrASkA5aGm9mGl87OfjeLPf17GN99s4/zzRzN27G1NKi9ZAFq0CAYNytnv3eYXa1KL9efk6AeRNgAl+yXUSqYkWbt2LZdeeikLFiyI+lCS+q//+i/ef/99JkyYUGdbrgQg9QHloUz6sBrjd7/L7l31u3bV/X5N7Otu81KNOhs0KPv9OS0lSUf97qGlrSAA9ejRgyuvvJJ//OMfOXkvUE1NTa0bcHORApBIa5A4ei0WeFpzG2SyX0Lp7mfIMRdddFHUh5DShRc2eirNFpc3Aagho1Qk99TXRJxs1G825Oz3eg43tdVH115uyaVul7wYBdehQwc2bNiQUx+8pBdbu6RDhw51tjVk9phMaF7O7Mro2tu6NbMZfyWt+q6pKOTFIAStyth6tfTqjTk7qCxnD6x+jb72Nm+GLVuC+4eg/mHckpFU15QGITST9u3b58TqfyL5JuNrr7w8mPwUgv6umJxtH5VM5EUAEmmof1tTzvfXTYfSLBWoL8zMxD6zVjTfnDReXjTBiTRU5f6l9Pm6ko4nZ6GDKX6odKLGBqZW2gSXVbH55mJTCSkQZZWa4ERyQFXHYoqz8UWfaqbqZEsZpBO78zafxe53ymQpCMlJqgGJxKncvxSA4i8XNt9OGrOEQjx9yQbS3ZSrzygjmoonCxSApClaJABJdsUCUmKTp4JRo0QRgPLiPiARacNiM3E/8cSe4KMbuloF9QGJJPj66+zNsKAf4S0ofpG8vJsssHVSABKJ0/1b2SurvrEGCkwiCkAitXy7R/DItUFwoKDVaIsWJa8J6YPMGQpAIs0kvkUoXiaD4HQPZiOlWqIiMforGEVKo+BE4uXoDZ85elitT3z0T3WjcJ4GJY2CExFpTrERc4mj5mIWLQrmoCsvj+LocpqZDTazFWZWZWa3Jtn+f8zsNTP7wMw+MrNz05WpACQi+Sk+GMUHJdAQ7gRmVgA8ApwD9ANGmlm/hGy3A7Pc/RhgBPBounIVgEREYsrKglpRbLXD2EM1ooFAlbuvcvftwAxgaEIeB2Jrk3cG1qQrNG0AMrMpZrbOzJbEpRWb2dtmVmlmFWY2MEw3M5sUVtE+MrMBce8ZbWYrw8fouPRjzezj8D2TLFw60cy6mNmCMP8CMzsg3bGKiDRZ4mqH+XNTa2H4fR57xHeEHQx8Fve6OkyLNw74kZlVA/OAa9PusAEHNRX4FfBUXNoDwF3u/vuwne8BggnszwH6ho/jgceA482sC3AnUEIQJf9oZnPdfVOYpwx4OzzowcDvgVuBV9z9vrC98VbglgYcr0iblGoJ8jztM28+icMXS0uTD+luex98jbuXpNiWbE31xBFsI4Gp7v6fZnYi8LSZHenuu1LtMG0NyN1fBzYm2XGyqtZQ4CkPvA3sb2Y9gLOBBe6+MQw6C4DB4bb93P0tD4bjPQWcH1fWtPD5tLh0kbyTagny/PlxHqFRo+oOVsi/D74aOCTudU/qNrGNAWYBuPtbQAegW32FZnof0PXAfDN7kCCI/WuYnqqaVl96dZJ0gG+5+1oAd19rZt1THUxYVSwD2GuvvTL8k0RyV6p7ijTjTAtI9uHn3wf/HtDXzHoDnxMMMki82eqvwBnAVDP7HkEA+nt9hWY6COEnwA3ufghwAzA5TE9VTWtseqO4e7m7l7h7SWGh7q0VkRZSXp4XAxXcvQa4BpgPLCcY7bbUzMab2ZAw28+AK83sQ+A54DJPc6Nppt/Wo4Gfhs+fB54Mn6eqplVTe5HjnsDCML1nkvwAX5hZj7D20wNYl+GxiohkX2XlnptZO3cO/m1bfUK1uPs8gn76+LQ74p4vA05qTJmZ1oDWALFG0dOBleHzucCl4Wi4E4Cvwma0+cBZZnZAOJrtLGB+uG2zmZ0Qjn67FJgTV1ZstNzouHQRkWjFOuUGDQruHUrWQSdppa0BmdlzBLWXbuHwujuBK4FfmlkhsI2w/4UgOp4LVAFbgcsB3H2jmU0gaEcEGO/usYENPyEYabcPwei334fp9wGzzGwMQdvihRn/lSIi2ZTYL5RfAxKyJm0AcveRKTYdmySvA1enKGcKMCVJegVwZJL0DQQdWiJSj1TDs+O1vRHD0haox16kFUs16XO8ZMs/KCBJLlAAEmnFUg3Pjpe4/ENl5Z73ikRJAUhap8RvVf2kTynZjf3SDBLbQvV/Mi0FIGmdpk8PLvji4syXGE0mVqZIYyS2haqa2SAKQNJ6FRcHU+hnssRofWU2pGNFJJ6qmRlRAJLWryEdISKScxSARPKQuiskFygAieSZxBbGTLvQFLTSSHeDlj5ALM1cca1OUVGRb9myJerDkOYWu7AXLozyKNqETLrQYlOgxa9SoO/TOOk+1GQfYDIt+KGa2VZ3L2qRncX2qQAkrZICUKQSv1/r+z5VYEqiIVE/NiKzhf6PKwBlgQJQnlAAyimpvk/T/dBXcKpHC/8fjyIAqQ9IRJos1UDE+n7o61YZUQASkWZT3wh53SrTAG18uKICkIhILsqD2RUUgEQkMslGKrexH/mZy4PZFRSARCQSyWY8aoM/8rMr3b1FxcUwcWKLHU5TKQCJSCSS9Q+1wR/52dMG5yhUABKRnBL7ka+muARtcM5DBSARyRmxH/mx6YGgzX3nShzdiCqtk25EbdPKy2Hs2OC5pvtpGboRVVpGNtfPiYoWjmvTYkEm2XQ/CkBth2pA+ai0tG18gevncF5Rpbd5qQYkLacFJzkUEUmmXdQHICIi+UkBSERajdgQ7fLyqI9EskFNcCLSKiQO0Y4foKDuwNZJNSARaRXKyoJuyyeeqD00u7Ky9Q/qzFcKQCLSqsQCUexRXBzUiNQs1/ooAIlIqxZrmhs7Nugfij0UkHJf2gBkZlPMbJ2ZLYlLm2lmleFjtZlVxm27zcyqzGyFmZ0dlz44TKsys1vj0nub2TtmtjIsd68wfe/wdVW4vVe2/mgRaTvKytQs11o1pAY0FRgcn+DuF7t7sbsXAy8ALwKYWT9gBHBE+J5HzazAzAqAR4BzgH7AyDAvwP3Aw+7eF9gEjAnTxwCb3L0P8HCYT0SkjmTNcpJdqSoRCXkuMrNlZrbUzNL+BEgbgNz9dWBjip0ZcBHwXJg0FJjh7t+4+1+AKmBg+Khy91Xuvh2YAQwN3386MDt8/zTg/LiypoXPZwNnhPlFRKQFpalExPL0BW4DTnL3I4Dr05Xb1D6gU4Av3H1l+Ppg4LO47dVhWqr0rsCX7l6TkF6rrHD7V2H+OsyszMwqzKyipqYmWRYREclc0kpEQp4rgUfcfROAu69LV2hTA9BI9tR+AJLVUDyD9PrKqpvoXu7uJe5eUlioW5tERDJQGPshHz7i76xKVYmIdxhwmJktNrO3zWwwaWT8bW1mhcAw4NiEgzok7nVPYE34PFn6emB/MysMaznx+WNlVYf76kyKpkAREWmyGncvSbGtIRWCQqAvUErwXf6GmR3p7l+m2mFTakDfB/7k7tVxaXOBEeEItt7hwbwLvAf0DUe87UUwUGGuB1NxvwYMD98/GpgTV9bo8Plw4FVva1N3i4i0DvVVLuLzzHH3HeEYgBUEMSClhgzDfg54CzjczKrNLDZKbQS1m99w96XALGAZ8AfganffGdZurgHmA8uBWWFegFuAG82siqCPZ3KYPhnoGqbfCCQddSEiIs0uaSUiIc9vgdMAzKwbQZPcqvoK1XpA+UgLq0gbp//ijZduPSAzOxeYCBQAU9z9HjMbD1S4+9xwlPJ/EtyCsxO4x91n1LdP9diLiEha7j4PmJeQdkfccydorbqxoWVqKh4REYmEApCItEmaoDT3qQlORNqcUaOCADR2bN054bR2UO5QABKRNicWYBKDT+JidgpG0dIouHykIUKSp8rL9wSfRYuCfwcNUiCC9KPgmmWfCkB5SAFIZHcwig9EkL/BKIoApEEIIpKXki3xHes30uCFlqEAJCJ5LX4toSeeCNK0mF3LUAASEQmVldVeWVWalwKQiIhEQgFIREQioQAkIiKR0I2oIiIJKiv33K0A+Ts0u7kpAImIxBk1qvbrysrgXwWg7FMAEhGJU1ZWO9jE14Qku9QHJCIikVAAEhGRSCgAiYhIJBSAREQkEgpAIiJpaHXV5qEAJCJSj9iw7LFjgxFxCkTZowAkIlKPsrI9SzbElmtQIMoOBSARkTQS1w6qrNSSDdmgACQi0kCxQFRcHPWRtA2aCSET8QvLt0aVlbqCRCRyqgFlYvr0PRNEtUbFxXUnvBIRaWGqAWWquDioi4uISEZUAxIRkUikDUBmNsXM1pnZkoT0a81shZktNbMH4tJvM7OqcNvZcemDw7QqM7s1Lr23mb1jZivNbKaZ7RWm7x2+rgq398rGHywiIrmhITWgqcDg+AQzOw0YCvR39yOAB8P0fsAI4IjwPY+aWYGZFQCPAOcA/YCRYV6A+4GH3b0vsAkYE6aPATa5ex/g4TCfiIi0EWkDkLu/DmxMSP4JcJ+7fxPmWRemDwVmuPs37v4XoAoYGD6q3H2Vu28HZgBDzcyA04HZ4funAefHlTUtfD4bOCPMLyIibUCmfUCHAaeETWOLzOy4MP1g4LO4fNVhWqr0rsCX7l6TkF6rrHD7V2H+OsyszMwqzKyipqYmWRYRkayKLdutGREyl+kouELgAOAE4Dhglpl9B0hWQ3GSBzqvJz9pttVOdC8HygGKioqS5hERyZbYXQxarrtpMq0BVQMveuBdYBfQLUw/JC5fT2BNPenrgf3NrDAhnfj3hNs7U7cpUESkxeXjjAipBpIlyTfczNzMStKVmWkA+i1B3w1mdhiwF0EwmQuMCEew9Qb6Au8C7wF9wxFvexEMVJjr7g68BgwPyx0NzAmfzw1fE25/NcwvIiItKM1Asvh8nYDrgHcaUm5DhmE/B7wFHG5m1WY2BpgCfCccmj0DGB3WhpYCs4BlwB+Aq919Z9iHcw0wH1gOzArzAtwC3GhmVQR9PJPD9MlA1zD9RiBlxBURiUqe9AUlHUiWJN8E4AFgW0MKTdsH5O4jU2z6UYr89wD3JEmfB8xLkr6K4I9LTN8GXJju+EREotLG+oIKzawi7nV52L8OyQeSHR//ZjM7BjjE3X9nZjc1aIdNOVoRkXxWVhY8SkujPpKsqHH3VP029Q4KM7N2BPdrXtaYHWoqHhERSSfVQLKYTsCRwEIzW00wQnpuuoEICkAiIpJO0oFksY3u/pW7d3P3Xu7eC3gbGOLuFcmLCygAiYhIvVINJDOz8WY2JNNy1QckIiJpJRtI5u53pMhb2pAyVQMSEZFIKACJiEgkFIBERCQSCkAiIhIJBSAREYmEApCIiERCAUhERCKhACQiIpFQABIRkUgoAImISCQUgEREJBIKQCIiEglNRioikgWxpbljRo1q9SukNjsFIBGRJootzR3TRpbobnYKQCIiTRRbmjumjSzR3ezUByQiIpFQABIRkUgoAImISCQUgEREJBIKQCIiEgkFIBERiYQCkIiIREIBSEREIqEAJCIikUgbgMxsipmtM7MlcWnjzOxzM6sMH+fGbbvNzKrMbIWZnR2XPjhMqzKzW+PSe5vZO2a20sxmmtleYfre4euqcHuvbP3RIiISvYbUgKYCg5OkP+zuxeFjHoCZ9QNGAEeE73nUzArMrAB4BDgH6AeMDPMC3B+W1RfYBIwJ08cAm9y9D/BwmE9ERNqItAHI3V8HNjawvKHADHf/xt3/AlQBA8NHlbuvcvftwAxgqJkZcDowO3z/NOD8uLKmhc9nA2eE+UVEpA1oSh/QNWb2UdhEd0CYdjDwWVye6jAtVXpX4Et3r0lIr1VWuP2rML+IiLQBmQagx4B/AYqBtcB/hunJaiieQXp9ZdVhZmVmVmFmFTU1NcmyiIhIjskoALn7F+6+0913Ab8maGKDoAZzSFzWnsCaetLXA/ubWWFCeq2ywu2dSdEU6O7l7l7i7iWFhVphQkSkNcgoAJlZj7iXPwRiI+TmAiPCEWy9gb7Au8B7QN9wxNteBAMV5rq7A68Bw8P3jwbmxJU1Onw+HHg1zC8iIm1A2uqCmT0HlALdzKwauBMoNbNigiax1cBYAHdfamazgGVADXC1u+8My7kGmA8UAFPcfWm4i1uAGWZ2N/ABMDlMnww8bWZVBNhtw9gAAAgfSURBVDWfEU3+a0VEJGdYW6tUFBUV+ZYtW5p3J7HlDhcubN79iEir1Bq/Isxsq7sXteQ+1WEiItIMKivrX5p71Kjay3jnOjMbDPySoBXrSXe/L2H7jcCPCVq//g5c4e6f1lempuIREcmyUaOguDj19spKmD695Y6nqdJMJhDzAVDi7v0J7t18IF25qgGJiGRZWVn9tZv6akY5avdkAgBmNoNgsoBlsQzu/lpc/reBH6UrVDUgEREBKIzdTxk+4kNoqskEUhkD/D7tDjM7ThERaYrEPqIc6BOqcfeSFNsaMzHAj4ASYFC6HSoAiYi0sFGjar+urAz+zeFBCakmE6jFzL4P/BwY5O7fpCtUAUhEpIUl9hG1gj6h3ZMJAJ8T3JdZK4ya2THAE8Bgd1/XkELVByQiIvUKJ4SOTSawHJgVTjww3syGhNl+AXQEng/XiZubrlzVgEREckCy+4ZyoF9ot3Ddt3kJaXfEPf9+Y8tUABIRiVhinxC0in6hJlMAEhGJWLL7hlpBv1CTqQ9IREQioQAkIiKRUAASEZFIKACJiEgkFIBERCQSCkAiIhIJBSAREYmEApCIiERCAUhERCKhACQiIpFQABIRkUgoAImISCQUgEREJBIKQCIiEgktxyAikqMSF6nLpQXqskEBSEQkByUuUtcWF6hTABIRyUGJi9S1xQXq1AckIiKRUAASEZFIpA1AZjbFzNaZ2ZIk224yMzezbuFrM7NJZlZlZh+Z2YC4vKPNbGX4GB2XfqyZfRy+Z5KZWZjexcwWhPkXmNkB2fmTRUQkFzSkBjQVGJyYaGaHAGcCf41LPgfoGz7KgMfCvF2AO4HjgYHAnXEB5bEwb+x9sX3dCrzi7n2BV8LXIiLSRqQNQO7+OrAxyaaHgZsBj0sbCjzlgbeB/c2sB3A2sMDdN7r7JmABMDjctp+7v+XuDjwFnB9X1rTw+bS4dBERaQMy6gMysyHA5+7+YcKmg4HP4l5Xh2n1pVcnSQf4lruvBQj/7V7P8ZSZWYWZVdTU1GTwF4mISEtr9DBsM9sX+DlwVrLNSdI8g/RGcfdyoBygqKio0e8XEZGWl0kN6F+A3sCHZrYa6Am8b2YHEdRgDonL2xNYkya9Z5J0gC/CJjrCf9dlcKwiIpKjGh2A3P1jd+/u7r3cvRdBEBng7n8D5gKXhqPhTgC+CpvP5gNnmdkB4eCDs4D54bbNZnZCOPrtUmBOuKu5QGy03Oi4dBERaQMaMgz7OeAt4HAzqzazMfVknwesAqqAXwNXAbj7RmAC8F74GB+mAfwEeDJ8z5+B34fp9wFnmtlKgtF29zXuTxMRkVxmweCztqOoqMi3bNnSvDuJzYmxcGHz7kdEJNTcXztmttXdi5qn9OQ0E4KIiERCAUhERCKhACQiIpHQcgwx5eUwfXrD8lZWQnFx8x6PiEgOMbPBwC+BAuBJd78vYfveBLPZHAtsAC5299X1ldnmAtCh27ZltnDGokXBv4MGpc9bXFx3tSgRkTbKzAqARwhGJFcD75nZXHdfFpdtDLDJ3fuY2QjgfuDi+sptcwEoY4MGtb31bkWkTUlcojtRcTFMnNgsux4IVLn7KgAzm0EwX2d8ABoKjAufzwZ+ZWbm9Qy1bnMB6NMOHTQ8WkTanBZodCk0s4q41+XhNGeQfD7P4xPevzuPu9eY2VdAV2B9yh02+ZBFRKTZJS7R3Qxq3L0kxbaGzNvZ6Lk9NQpORETSSTWfZ9I8ZlYIdCb5Uj67KQCJiEg67wF9zay3me0FjCCYrzNe/Pydw4FX6+v/ATXBiYhIGmGfzjUEE0sXAFPcfamZjQcq3H0uMBl42syqCGo+I9KVq7ngREREc8GJiEj+UAASEZFIKACJiEgkFIBERCQSbW4QgpntAv7ZiLcUADuzsOvGltPQ/Ony1bc91bZk6YlphUBNA46vOWTjnGRSRnOfk6acD4junOgaSZ7W1q6Rfdy9ZSsl7p7XD4LpJlq8nIbmT5evvu2ptiVLT0wjGFrZas9JJmU09zlpyvmI8pzoGkmelo/XSLYfaoKD/46onIbmT5evvu2ptiVLz9bnkA3ZOJZMymjuc5LP5yOTcnSNpBbVNZJVba4JTrLDzCo89bxQEgGdk9yi89F0qgFJKuXps0gL0znJLTofTaQakIiIREI1IBERiYQCkIiIREIBSEREIqEAJI1mZt8xs8lmNjvqY8lXZlZkZtPM7NdmdknUxyO6LjKhAJRnzGyKma0zsyUJ6YPNbIWZVZnZrfWV4e6r3H1M8x5p/mnkuRkGzHb3K4EhLX6weaIx50TXReMpAOWfqcDg+AQzKwAeAc4B+gEjzayfmR1lZr9LeHRv+UPOG1Np4LkhWBL5szBbNqbJkeSm0vBzIo2kFVHzjLu/bma9EpIHAlXuvgrAzGYAQ939XuDfWvYI81djzg1QTRCEKtEPyWbTyHOyrGWPrvXTf1wBOJg9v6Yh+HI7OFVmM+tqZo8Dx5jZbc19cHku1bl5EbjAzB4jB6ZUyTNJz4mui8ZTDUgALElayjuU3X0D8O/NdzgSJ+m5cfctwOUtfTACpD4nui4aSTUggeAX3CFxr3sCayI6FqlN5yb36JxkiQKQALwH9DWz3ma2FzACmBvxMUlA5yb36JxkiQJQnjGz54C3gMPNrNrMxrh7DXANMB9YDsxy96VRHmc+0rnJPTonzUuTkYqISCRUAxIRkUgoAImISCQUgEREJBIKQCIiEgkFIBERiYQCkIiIREIBSEREIqEAJCIikVAAEhGRSPx/Y93Rp0OsHYUAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig = plt.figure()\n",
"ax1 = fig.add_subplot(111)\n",
"ax2 = ax1.twinx()\n",
"ax1.step([part.resolution_parameter for part in resolution_profile], \n",
" [part.total_weight_in_all_comms() for part in resolution_profile], 'blue', label='Internal weight')\n",
"ax2.step([part.resolution_parameter for part in resolution_profile], \n",
" [part.compare_to(class_partition, 'nmi') for part in resolution_profile], 'red', label='NMI (classes)')\n",
"ax1.legend(loc='upper left')\n",
"ax2.legend(loc='upper right')\n",
"plt.xscale('log')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This reveals that resolutions somewhere between 0.1 and 0.3 seem to give identical partitions. Hence, a partition in this range is not very sensitive to small changes in the resolution parameter. This might be an indication that there is something of interest at that range."
]
},
{
"cell_type": "code",
"execution_count": 90,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 90,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"resolution_parameter = 0.2\n",
"partition = leidenalg.find_partition(G, leidenalg.CPMVertexPartition, weights='weight', resolution_parameter=resolution_parameter)\n",
"ig.plot(partition)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If we take a closer look, we see that the classes can actually be divided in three groups: PC/PSI (Physics & Chemistry/engineering), BIO (biology) and MP (Mathematics & Physics). This suggest that the pupils tend to connect more to other people with a similar educational background or interest."
]
},
{
"cell_type": "code",
"execution_count": 91,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[Counter({'PSI*': 32, 'PC': 44, 'PC*': 39}),\n",
" Counter({'2BIO1': 35, '2BIO2': 34, '2BIO3': 40}),\n",
" Counter({'MP*1': 29, 'MP*2': 38, 'PSI*': 1, 'MP': 33}),\n",
" Counter({'PSI*': 1}),\n",
" Counter({'2BIO1': 1})]"
]
},
"execution_count": 91,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from collections import Counter\n",
"[Counter(H.vs['class']) for H in partition.subgraphs()]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Negative links"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Modularity does not work correct if there are negative weights. However, this can be corrected for. Let us again work with the international relations network."
]
},
{
"cell_type": "code",
"execution_count": 92,
"metadata": {},
"outputs": [],
"source": [
"G = ig.Graph.Read('data/international_relations_2000_2015.gml')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In order to properly detect communities with negative links, we need to split the network into two parts: one with all the positive weights, and one with all the negative weights. Then we want to maximize the number of internal links in the positive network, and minimize the number of internal links in the negative network. This is represented below by the respectively positive and negative `layer_weight`."
]
},
{
"cell_type": "code",
"execution_count": 93,
"metadata": {},
"outputs": [],
"source": [
"G_positive = G.subgraph_edges(G.es.select(weight_gt=0), delete_vertices=False)\n",
"G_negative = G.subgraph_edges(G.es.select(weight_lt=0), delete_vertices=False)\n",
"G_negative.es['weight'] = -np.array(G_negative.es['weight'])\n",
"\n",
"partition_positive = leidenalg.ModularityVertexPartition(G_positive, weights='weight')\n",
"partition_negative = leidenalg.ModularityVertexPartition(G_negative, weights='weight')\n",
"opt = leidenalg.Optimiser()\n",
"diff = opt.optimise_partition_multiplex([partition_positive, partition_negative], layer_weights=[1.0, -1.0])\n",
"mod_partition = ig.VertexClustering(G, membership=partition_positive.membership)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us examine the results."
]
},
{
"cell_type": "code",
"execution_count": 94,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 94,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"G['layout'] = G.layout_fruchterman_reingold()\n",
"ig.plot(mod_partition, vertex_size=8, \n",
" edge_color=['blue' if e['weight'] > 0 else 'red' for e in G.es])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Modularity still has some issues with degree. For example, the long paths connected through negative links are sometimes put in the same community. But this is a bit strange, since they should be in separate communities.\n",
"\n",
"Again, CPM seems to work somewhat better in this regard. We can simply work with CPM directly, as it has no problems dealing with negative weights. By setting the resolution parameter to 0, we find communities that are positive connected within, but negatively connected in between. In fact, this corresponds to finding communities that are socially balanced."
]
},
{
"cell_type": "code",
"execution_count": 95,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 95,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"CPM_partition = leidenalg.find_partition(G, leidenalg.CPMVertexPartition, weights='weight', resolution_parameter=0)\n",
"\n",
"G['layout'] = G.layout_fruchterman_reingold()\n",
"ig.plot(CPM_partition, vertex_size=8, \n",
" edge_color=['blue' if e['weight'] > 0 else 'red' for e in G.es])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Conclusion"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This concludes the introduction in the network aspects of computational social science. You can keep this notebook to look up things during the assignment. Also, don't forget you can always can help within the notebook using Shift-Tab and Tab. If you need more help on `igraph`, you can also go to either the [reference manual](http://igraph.org/python/doc/tutorial/tutorial.html) or the [introduction](http://igraph.org/python/doc/tutorial/tutorial.html). Finally, you can always ask one of the lecturers to help you out."
]
}
],
"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.3"
}
},
"nbformat": 4,
"nbformat_minor": 1
}