diff --git a/m10-epidemics/m10-networkepidemiology.ipynb b/m10-epidemics/m10-networkepidemiology.ipynb deleted file mode 100755 index ada370f..0000000 --- a/m10-epidemics/m10-networkepidemiology.ipynb +++ /dev/null @@ -1 +0,0 @@ -{"cells":[{"attachments":{},"cell_type":"markdown","metadata":{"id":"_oaQGMHSBltD"},"source":["\n"," \n"," \n","
\n"," Run in Google Colab\n"," \n"," View on Github\n","
"]},{"attachments":{},"cell_type":"markdown","metadata":{"id":"et3I3l-iBltF"},"source":["# Network epidemiology\n","\n","The goal of this assignment is to implement a working SIS model on a network and reproduce a figure from this week's readings. We will start with an SI model as a template which you can then expand upon to build the SIS model.\n","\n","## SI model\n","\n","We will put together a function that will be able to take as input networkx graphs. In the SI model, a fraction of nodes begin as infected and then new infections spread throughout the network across the links with some probability $\\beta$. If the dice rolls out of the nodes favor then it becomes infected and can then infect other nodes that it neighbors.\n","\n","In a connected graph, the SI model should eventually infect everyone as $t \\rightarrow \\infty$ because there will always be a non-zero probability of transmission. Alternatively, the SIS model will reach an equilibrium point where there is a balance between infections and reversions to susceptibility.\n","\n","Let's build the SI model!"]},{"cell_type":"code","execution_count":2,"metadata":{"executionInfo":{"elapsed":408,"status":"ok","timestamp":1648464750012,"user":{"displayName":"Yong Yeol Ahn","photoUrl":"https://lh3.googleusercontent.com/a/default-user=s64","userId":"00405984523192108505"},"user_tz":240},"id":"W2gd8rYsBltF"},"outputs":[],"source":["# We will be using numpy and networkx for our function\n","import networkx as nx\n","import numpy as np"]},{"attachments":{},"cell_type":"markdown","metadata":{"id":"lc9CvcP-CHtu"},"source":["The simulation first requires a network. Once a network is given, we can then randomly choose a small fraction of initially infected nodes. We also need to keep track of every node's state (infected or susceptible). In each time step, we go through every node and figure out what would be the next state of the node based on its current state and neighbors. For instance, if the node is already infected, then it will remain infected regardless of the states of the neighbors (SI model). \n","\n","We can play with each step first. Let's first create a network. "]},{"cell_type":"code","execution_count":4,"metadata":{"executionInfo":{"elapsed":416,"status":"ok","timestamp":1648465079814,"user":{"displayName":"Yong Yeol Ahn","photoUrl":"https://lh3.googleusercontent.com/a/default-user=s64","userId":"00405984523192108505"},"user_tz":240},"id":"7WjEkBsDJ5Ic"},"outputs":[],"source":["G = nx.barabasi_albert_graph(10**4, 5)"]},{"attachments":{},"cell_type":"markdown","metadata":{"id":"7U8uddN0LPf3"},"source":["We can create an attribute to store the nodes' states. "]},{"cell_type":"code","execution_count":7,"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"elapsed":263,"status":"ok","timestamp":1648465170308,"user":{"displayName":"Yong Yeol Ahn","photoUrl":"https://lh3.googleusercontent.com/a/default-user=s64","userId":"00405984523192108505"},"user_tz":240},"id":"xzPd867cLPRQ","outputId":"9adbd7a0-343e-4475-c631-87a52d095c55"},"outputs":[{"data":{"text/plain":["[0]"]},"execution_count":7,"metadata":{},"output_type":"execute_result"}],"source":["nx.set_node_attributes(G, {node: [0] for node in G.nodes()}, 'inf')\n","G.nodes[0]['inf']"]},{"attachments":{},"cell_type":"markdown","metadata":{"id":"YPA1eSbCLmGg"},"source":["We are using a list to store the state so that we can keep a time-series of every node's state. This also makes the simulation easier. \n","\n","Now we can randomly choose a certain fraction of nodes and mark them as initially infected using `choice()` function in numpy. "]},{"cell_type":"code","execution_count":10,"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"elapsed":8,"status":"ok","timestamp":1648465374102,"user":{"displayName":"Yong Yeol Ahn","photoUrl":"https://lh3.googleusercontent.com/a/default-user=s64","userId":"00405984523192108505"},"user_tz":240},"id":"aIMIziOZL-8w","outputId":"dc2b6bf8-30b8-48dc-8c1f-e7b8dbc32fd3"},"outputs":[{"data":{"text/plain":["{1455, 3868, 3878, 5390, 5753, 7700, 7901, 8819, 9000, 9647}"]},"execution_count":10,"metadata":{},"output_type":"execute_result"}],"source":["init_infected_fraction = 0.001\n","init_infected = set(np.random.choice(G.nodes(), \n"," size=int(len(G) * init_infected_fraction), \n"," replace=False))\n","init_infected"]},{"attachments":{},"cell_type":"markdown","metadata":{"id":"Da9t1uoCMgEq"},"source":["We should set these nodes' states to be `1` and then we can begin the simulation. "]},{"cell_type":"code","execution_count":11,"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"elapsed":510,"status":"ok","timestamp":1648465499752,"user":{"displayName":"Yong Yeol Ahn","photoUrl":"https://lh3.googleusercontent.com/a/default-user=s64","userId":"00405984523192108505"},"user_tz":240},"id":"RgO-hI6AMpUq","outputId":"afc43957-b2c4-4fc8-d03b-7444e085c9ed"},"outputs":[{"data":{"text/plain":["[1]"]},"execution_count":11,"metadata":{},"output_type":"execute_result"}],"source":["for node in init_infected:\n"," G.nodes[node]['inf'][0] = 1\n","G.nodes[1455]['inf']"]},{"attachments":{},"cell_type":"markdown","metadata":{"id":"Cuh5w7hPNULr"},"source":["But how should we update the states and run the simulation? \n","\n","There is an important issue here. If we immediately change the state of a particular node, that may affect the updating of other nodes! Imagine an extreme case where we are simulating a highly infectious disease on a chain-like network where node 0 is connected to node 1, node 1 is connected to node 0 and 2, and so on. In the simulation we are going through each node one by one from node 0 to node N and update them immediately. But it happens to be that node 0 was infected. Then, it can propagate to many nodes within a single time step because it is possible that node 1 gets updated to be infected by node 0, and then node 2 gets infected by node 1, and so on. This is clearly unrealistic. \n","\n","If we change a node's state, then the next node will be updating with respect to a network that is now in a different state!\n","\n","When modelling discrete-time dynamical systems there are generally two different update strategies: synchronous and asynchronous updating. In the asynchronous updating, a random node is picked and its state is updated according to the current network state. In the synchronous updating, there is a global time clock that all nodes are synced to, so nodes only update according to the state of the network at the _current time-step_ and all nodes move to the next time step simultaneously. \n","\n","Choosing the updating scheme can have a huge impact on dynamics. We will be using the synchronous updating scheme, which means we need to store the _next state_ of the system before updating everyone all at once. There are many ways to accomplish this, but for the sake of simplicity, we will just keep the whole history. "]},{"cell_type":"code","execution_count":12,"metadata":{"executionInfo":{"elapsed":208,"status":"ok","timestamp":1648467910491,"user":{"displayName":"Yong Yeol Ahn","photoUrl":"https://lh3.googleusercontent.com/a/default-user=s64","userId":"00405984523192108505"},"user_tz":240},"id":"BunVNgrZBltG"},"outputs":[],"source":["def run_SI(graph, tmax: int, beta: float, initial_inf: float):\n"," \"\"\"Runs the SI model on the given graph with given parameters. \n","\n"," Parameters\n"," ----------\n"," graph : networkx graph object\n"," The network (graph) on which the simulation will be run\n"," tmax : int\n"," The maximum time step that we will run the model \n"," beta : float\n"," The transmission probability\n"," initial_inf : float\n"," The initial fraction of infected nodes\n","\n"," Returns\n"," -------\n"," list[float]\n"," the time-series of the fraction of infected nodes\n"," \"\"\"\n"," # First lets generate a set of initially infected nodes.\n"," init_infected = set(\n"," np.random.choice(graph.nodes(), size=int(len(graph) * initial_inf), replace=False)\n"," )\n"," \n"," # The code below uses a dictionary comprehension to generate a dictionary\n"," # with keys=nodes and values=a list of 0's and 1's. The 1 is for infected\n"," # and 0 is for susceptible. We then give that dictionary to networkx's\n"," # attribute function which then gives all the nodes the 'inf' attribute.\n"," nx.set_node_attributes(\n"," graph, \n"," {node: ([1] if node in init_infected else [0]) for node in graph.nodes()},\n"," 'inf'\n"," )\n"," \n"," # Now we need to loop through for `tmax` time step. One time step equals to \n"," # updating the whole network once. \n"," for t in range(tmax):\n"," for node in graph.nodes():\n"," \n"," # Now we check if the node if susceptible to infection\n"," # If it is, we need to determine the probability of it switching\n"," # and then switch it for the next time-step\n"," if graph.nodes[node]['inf'][t] == 0:\n"," \n"," # First determine how many infected neighbors the node has at time t:\n"," num_inf_neighbors = np.sum(\n"," [ graph.nodes[neighbor]['inf'][t] for neighbor in graph.neighbors(node)]\n"," )\n"," \n"," # Instead of drawing a bunch of random numbers for each neighbor\n"," # we can just calculate the cumulative probability of getting\n"," # infected since these events are independent and then just\n"," # draw 1 random number to check against:\n"," if np.random.random() < (1 - (1 - beta)**num_inf_neighbors):\n"," # If infection occurs we add a 1 to the state list of the node.\n"," # Note that by doing this we don't change how the other \n"," # nodes update, because they will be using time index t not t+1\n"," graph.nodes[node]['inf'].append(1)\n"," \n"," else:\n"," # If no infection occurs, then just append the current state (0)\n"," graph.nodes[node]['inf'].append(0)\n"," \n"," # Similarly, if the node is already infected it can't change back\n"," # So we append the current state if it wasn't susceptible\n"," else:\n"," graph.nodes[node]['inf'].append(1)\n"," # the function returns a time series of the fraction of infected in each time step. \n"," return [ np.mean([ graph.nodes[node]['inf'][t] for node in graph.nodes() ]) for t in range(tmax)]"]},{"attachments":{},"cell_type":"markdown","metadata":{"id":"Z-DB-9N4BltH"},"source":["And there we have our SI model. The function is mostly comments, there are only a dozen lines of code involved in the whole process. Lets give it a run:"]},{"cell_type":"code","execution_count":13,"metadata":{"executionInfo":{"elapsed":195,"status":"ok","timestamp":1648467998110,"user":{"displayName":"Yong Yeol Ahn","photoUrl":"https://lh3.googleusercontent.com/a/default-user=s64","userId":"00405984523192108505"},"user_tz":240},"id":"7vJ4gcT0BltH"},"outputs":[],"source":["# Lets generate a random graph for testing\n","rnd_graph = nx.erdos_renyi_graph(100, 0.1)\n","\n","# We want to make sure that the graph is connected, so we will only take the largest\n","# connected component, as disconnected parts can't be infected or transmit infection:\n","largest_component = max(nx.connected_components(rnd_graph), key=len)\n","# above returns a set of nodes, so we use it to creat a subgraph\n","largest_component = rnd_graph.subgraph(largest_component)"]},{"cell_type":"code","execution_count":15,"metadata":{"colab":{"base_uri":"https://localhost:8080/","height":57},"executionInfo":{"elapsed":9,"status":"ok","timestamp":1648468005007,"user":{"displayName":"Yong Yeol Ahn","photoUrl":"https://lh3.googleusercontent.com/a/default-user=s64","userId":"00405984523192108505"},"user_tz":240},"id":"RqFS7albBltH","outputId":"94922e5c-40e7-4d67-de4e-1b4644baad14"},"outputs":[{"data":{"application/vnd.google.colaboratory.intrinsic+json":{"type":"string"},"text/plain":["'Graph with 100 nodes and 486 edges'"]},"execution_count":15,"metadata":{},"output_type":"execute_result"}],"source":["nx.info(largest_component)"]},{"attachments":{},"cell_type":"markdown","metadata":{"id":"5o2NLrC1Wgot"},"source":["Now we can run it and plot it. "]},{"cell_type":"code","execution_count":16,"metadata":{"executionInfo":{"elapsed":310,"status":"ok","timestamp":1648468021624,"user":{"displayName":"Yong Yeol Ahn","photoUrl":"https://lh3.googleusercontent.com/a/default-user=s64","userId":"00405984523192108505"},"user_tz":240},"id":"kbaM5FSxBltH"},"outputs":[],"source":["import matplotlib.pyplot as plt"]},{"cell_type":"code","execution_count":18,"metadata":{"colab":{"base_uri":"https://localhost:8080/","height":303},"executionInfo":{"elapsed":1081,"status":"ok","timestamp":1648468028065,"user":{"displayName":"Yong Yeol Ahn","photoUrl":"https://lh3.googleusercontent.com/a/default-user=s64","userId":"00405984523192108505"},"user_tz":240},"id":"3qjDbTvVBltH","outputId":"3bb26a59-fd89-4339-a190-3b62923222ce"},"outputs":[{"data":{"text/plain":["[]"]},"execution_count":18,"metadata":{},"output_type":"execute_result"},{"data":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAf5UlEQVR4nO3deXxU9b3G8c+XsIYdEvYlIIR9D4uIYqsooqJV24Ki4gLaW711q2Jt1VrbKlavtXptUVGKyuKOiKLXuqFllT0QCBAIkSUhkJCdJL/7RwYbQ0IGMsmZmTzv1ysvZub8knk4nHk4OfM7Z8w5h4iIhL46XgcQEZHAUKGLiIQJFbqISJhQoYuIhAkVuohImKjr1RNHRUW5mJgYr55eRCQkrVmzJs05F13eMs8KPSYmhtWrV3v19CIiIcnMdle0TIdcRETCRKWFbmazzeygmW2qYLmZ2TNmlmhmG8xsaOBjiohIZfzZQ38FGH+S5RcBPX1f04Hnqx5LREROVaWF7pz7Ekg/yZDLgH+6EsuBFmbWPlABRUTEP4E4ht4RSC51f6/vsROY2XQzW21mq1NTUwPw1CIiclyNvinqnJvlnItzzsVFR5c760ZERE5TIAo9Behc6n4n32MiIlKDAjEPfRFwm5nNB0YCGc65fQH4uSIiVeac43DOMZIOZbP7UDbJ6bkUFhV7mum8Pm0Z1LlFwH9upYVuZvOAc4EoM9sLPATUA3DO/R1YAkwAEoEc4IaApxQROQnnHKlH80k6lPN9cScdymH3oWx2H8rhaF7hD8abeRTUp02zht4UunNuciXLHfDLgCUSESlHcbFjX2Yeu9N8ZZ2eze604wWeQ+6xou/HRtQxOrdsRNfWjRnWpSVdWzcmJiqSLq0a07lVIxrUjfDwb1J9PDv1X0SkrMKiYlKO5JJ0KIc9pfaykw7lsCc9h4LC/xwqqR9Rh86tGhHTujGjz4giJiqypLhbR9KhRSPqRdS+E+FV6CK1TFGxw8uPniwsdqQcyS0p6rScHxwe2Xs4l8Li/2RrWK8OMa0b0z2qMT/u3YaurSPp2qpkb7t980ZE1PH42EmQUaGLhLniYkf8vky+2JbKl9tSWbP78A9K02tNGtSla+tI+nVozoQB7Ylp3ZiurSOJiWpMm6YNMK8PeIcQFbpIGErPLuCr7am+Ek8jLSsfgH4dmnHjmG40beDdS79OHaN984bfHx5p1bi+SjtAVOgiYaCwqJj1e4/wRUJJiW9IycA5aBlZj3NioxkbG83ZPaOJbtrA66hSjVToIiFqX0YuX24rKfBl29PIzCukjsGQLi258/xYxsZG079jcx1nrkVU6CIhJD27gH98uYPPt6aScOAoAO2aNeSi/u0Z2yuas86IonlkPY9TildU6CIh4ottqdzzxnoOZxcwsnsrrhzWm7GxbYht20THoAVQoYsEvbxjRfx5yRbm/Hs3sW2b8MoNw+nXobnXsSQIqdBFgtimlAzuWLCOxINZ3HBWDPeN703DeuF5lqNUnQpdJAgVFTv+8eUO/ueTbbRqXJ+5N43g7J665LScnApdJMgkp+dw98L1rExKZ8KAdvzpJwNoEVnf61gSAlToIkHCOcc7a1N48L3NADz500FcMbSj3vAUv6nQRYLAkZwCHnhnEx9s3MfwmJY89bPBdG4V6XUsCTEqdBGPLduext1vrCM9u4B7x/filnPO0MlAclpU6CIeyTtWxMyPEpj99S56tGnCS9cPp39HTUeU06dCF/FA/HeZ3LFgLdsOZDF1dAwzLtJ0RKk6FbpIDcrMO8brK/bw1MfbaB5Zjzk3jmBsrKYjSmCo0EWqUXGxY/N3mXy5PZUvElJZs+cwRcWO8f3a8acrBtCqsaYjSuCo0EUCLC0rn2Xb077/QIlD2QUA9O/YjFvHdufcXm2I69pS0xEl4FToIlVUWFTM2uT/XIt8Y0oGAK0a1+ecnlGM7RXNmB66FrlUPxW6yGlIOeK7FnlCKl8npnE0v5CIOsbQLi24e1wsY3tF079Dc+po+qHUIBW6yCnYlZbNL15dw9b9Jdci79C8IRcPbM/Y2GhG94iieSNdi1y8o0IX8VNxsePeN9fz3ZFcfntxH8bGRtOjja5FLsFDhS7ip9dW7GZV0mGeuGogP43r7HUckRPU8TqASChIOZLLYx9u5eyeUVw1rJPXcUTKpUIXqYRzjgfe2YgD/vSTATrEIkFLhS5SiXfXpfB5Qiq/vrCXroAoQU2FLnISaVn5/P79eIZ2acF1Z8Z4HUfkpFToIifx8KLN5OQX8fiVA3VJWwl6KnSRCny8eT+LN+zj9h/3oGfbpl7HEamUCl2kHBm5x/jde5vo3a4pt4w9w+s4In7RPHSRcjz24RZSj+bzwnVx1K+r/R4JDdpSRcr4JjGNeSuTmXZ2dwZ2auF1HBG/qdBFSsktKGLG2xuJaR3JHefHeh1H5JTokItIKU99ksCe9BzmTRtFo/r6SDgJLX7toZvZeDNLMLNEM5tRzvIuZvaZma01sw1mNiHwUUWq17rkI7y0bBdXj+zCmWe09jqOyCmrtNDNLAJ4DrgI6AtMNrO+ZYb9FljonBsCTAL+N9BBRapTQWEx9725gTZNGzLjot5exxE5Lf7soY8AEp1zO51zBcB84LIyYxzQzHe7OfBd4CKKVL/nP99BwoGjPHp5f5o11DXNJTT5U+gdgeRS9/f6HivtYWCKme0FlgC3l/eDzGy6ma02s9WpqamnEVck8LYdOMqzn21n4qAOnN+3rddxRE5boGa5TAZecc51AiYAc83shJ/tnJvlnItzzsVFR0cH6KlFTl9RseO+tzbQpEFdHrq07JFEkdDiT6GnAKWv5t/J91hpNwELAZxz/wYaAlGBCChSneZ8k8TaPUd4eGI/WjfRhzhLaPOn0FcBPc2sm5nVp+RNz0VlxuwBzgMwsz6UFLqOqUhQS07P4YmlCfy4dxsmDurgdRyRKqu00J1zhcBtwFJgCyWzWTab2SNmNtE37G5gmpmtB+YBU51zrrpCi1SVc477395IRB3j0cv760MrJCz4dWKRc24JJW92ln7swVK344GzAhtNpPq8sWYvyxLTePTy/nRo0cjrOCIBoVP/pdY5mJnHo4vjGdGtFVeP6OJ1HJGAUaFLrfPge5vJLyzmsSsGUEcfWiFhRIUutcprK3bz0eb93Dkulu7RTbyOIxJQKnSpNT6JP8Dv3t3Ej3pFc/OYbl7HEQk4FbrUCmt2H+b2ed8yoGNznrtmKHUjtOlL+NFWLWFvR2oWN89ZRbtmDZk9dTiR9XXVaAlPKnQJaweP5nH97JVE1DHm3DhCZ4NKWNOuioSto3nHuOHlVaRnFzB/+ii6tm7sdSSRaqVCl7BUUFjML179lq37j/LS9XH6bFCpFXTIRcKOcyVXUFyWmMZjVwzg3F5tvI4kUiNU6BJ2Zi5N4J21KdxzQSw/jetc+TeIhAkVuoSVOd8k8fznO5gyqgu//FEPr+OI1CgVuoSNDzfu4+H3N3NB37b8fqKuoCi1jwpdwsLKXen8asE6hnZpyTOThxCha7RILaRCl5C3/cBRbp6zik4tG/HidXE0rBfhdSQRT6jQJaTty8jl+tkraVAvgjk3jKBl4/peRxLxjApdQlZG7jGmzl5FZl4hr9wwnM6tIr2OJOIpFbqEpPzCIm6Zu5qdaVn849ph9OvQ3OtIIp7TmaIScoqLHXcvXM/ynek8/fPBnNUjyutIIkFBe+gScv64ZAuLN+zj/ot6c/mQjl7HEQkaKnQJKXOX7+alZbuYOjqG6ed09zqOSFBRoUvIWLM7nUfe38yPekXzu0v66sQhkTJU6BISDmbmceur39KhRSOe/rlOHBIpj94UlaBXUFjML177lqy8Ql69aSTNI+t5HUkkKKnQJej9YXE8a3Yf5tmrh9CrXVOv44gELR1ykaD2xupk5i7fzfRzunPJwA5exxEJaip0CVob92bwwLubOKtHa+69sJfXcUSCngpdgtKhrHxufXUN0U0a8MykIdSN0KYqUhkdQ5egU1hUzO3z1pKalc9bt46mdZMGXkcSCQna7ZGgM3NpAt/sOMQfL+/PgE66RouIv1ToElTeX/8ds77cybWjuurzQEVOkQpdgkbC/qPc++YG4rq25HeX9PU6jkjIUaFLUMjIPcYtc1fTpGFd/veaodSvq01T5FTpVSOeKy523DF/LXsP5/L8NUNp06yh15FEQpIKXTz39Kfb+SwhlYcu7UtcTCuv44iELL8K3czGm1mCmSWa2YwKxvzMzOLNbLOZvR7YmBKuPok/wDOfbueqYZ2YMqqr13FEQlql89DNLAJ4DhgH7AVWmdki51x8qTE9gfuBs5xzh82sTXUFlvCxMzWLuxasY0DH5jx6eX9dDlekivzZQx8BJDrndjrnCoD5wGVlxkwDnnPOHQZwzh0MbEwJN1n5hdwydw316tbh79cOo2G9CK8jiYQ8fwq9I5Bc6v5e32OlxQKxZva1mS03s/GBCijhxznHvW+uZ0dqFs9OHkLHFo28jiQSFgJ16n9doCdwLtAJ+NLMBjjnjpQeZGbTgekAXbp0CdBTS6j5+xc7WbJxP7+Z0JvR+oBnkYDxZw89BSh9yl4n32Ol7QUWOeeOOed2AdsoKfgfcM7Ncs7FOefioqOjTzezhLCPNu3niaVbuWRge6adrc8EFQkkfwp9FdDTzLqZWX1gErCozJh3Kdk7x8yiKDkEszOAOSUMfLBhH798/VsGdW7B41cO1JugIgFWaaE75wqB24ClwBZgoXNus5k9YmYTfcOWAofMLB74DPi1c+5QdYWW0LNo/Xf89/y1DOncgn/eOILGDXShT5FAM+ecJ08cFxfnVq9e7clzS816Z+1e7l64nuExrZg9dbjKXKQKzGyNcy6uvGU6U1Sq1cLVydy1cD2jurfm5RtU5iLVSa8uqTbzVu7h/rc3cnbPKF64Lk5zzUWqmfbQpVrMXb6b+9/eyLm9olXmIjVEe+gScK98vYuH34/n/D5teO6aoTSoqzIXqQkqdAmoF7/ayaMfbOHCfm3522Rd11ykJqnQJWCe/3wHj3+0lYsHtOfpSYOpF6EyF6lJKnQJiGf/tZ2/fLyNiYM68NTPBlFXZS5S41ToUiXOOf766Xae/r/tXDGkI0/8dBARdXQGqIgXVOhy2pxzPPnxNp79LJGrhnXi8SsHqsxFPKRCl9PinOOxj7byjy92MnlEZ/54+QDqqMxFPKVCl1PmnOPRD7bw0rJdTBnVhUcm9leZiwQBFbqcEuccv38/nle+SWLq6BgeurSvrpooEiRU6HJK/v7FTl75Jombx3TjgYv7qMxFgojmlonfNuw9wpMfJzBhQDuVuUgQUqGLX3IKCvnV/HVEN23An34yQGUuEoR0yEX88sj78SQdyub1m0fRIrK+13FEpBzaQ5dKfbRpH/NXJXPr2DM484zWXscRkQqo0OWk9mfkMePtjQzo2Jw7z4/1Oo6InIQKXSpUXOy4a+E68o8V89dJg3XlRJEgp1eoVOiFr3byzY5DPHRpX7pHN/E6johUQoUu5dqUksFfPk5gfL92/Hx4Z6/jiIgfVOhygpyCQv57/lpaN27An6/QFEWRUKFpi3KCPyzewq60bF67aSQtG2uKokio0B66/MDSzfuZt3IP08/pzugeUV7HEZFToEKX7x3IzGPGWxvo37EZd4/r5XUcETlFKnQBSqYo3r1wPbnHivjrpCGaoigSgvSqFQBeWraLZYlpPHhJP87QFEWRkKRCFzalZDBz6VYu6NuWySM0RVEkVKnQa7ncgiJ+NX8tLSPr89iVAzVFUSSEadpiLffoB/HsSM3m1ZtG0kpTFEVCmvbQa7FP4g/w2oo9TDu7G2N6aoqiSKhToddSBzPzuO+tDfRt34x7LtQURZFwoEKvhYqLHXe/sZ6cgkKemTyYBnUjvI4kIgGgQq+FZn+9i6+2p/Hbi/vSo01Tr+OISICo0GuZL7elMvOjBM7v05ZrRnbxOo6IBJBmudQSeceKePyjrbz8dRI92jTh8St1FUWRcOPXHrqZjTezBDNLNLMZJxl3pZk5M4sLXESpqvjvMpn47DJe/jqJqaNjWHz7GFo3aeB1LBEJsEr30M0sAngOGAfsBVaZ2SLnXHyZcU2BXwErqiOonLqiYseLX+3kLx8n0CKyPnNuHMHY2GivY4lINfHnkMsIINE5txPAzOYDlwHxZcb9AXgc+HVAE8ppSTmSy10L1rFiVzrj+7XjT1cM0IlDImHOn0LvCCSXur8XGFl6gJkNBTo75z4wMxW6x95bl8Jv391EcbFj5lUD+emwTjpeLlILVPlNUTOrAzwFTPVj7HRgOkCXLpphEWgZOcf47XubeH/9dwzr2pL/+dlgurSO9DqWiNQQfwo9BSh9Cb5OvseOawr0Bz737QW2AxaZ2UTn3OrSP8g5NwuYBRAXF+eqkFvK+CYxjbvfWE/q0XzuuSCWW8eeQd0IzUoVqU38KfRVQE8z60ZJkU8Crj6+0DmXAXx/IRAz+xy4p2yZS/XIO1bEX5Ym8OKyXXSPaszb/zWagZ1aeB1LRDxQaaE75wrN7DZgKRABzHbObTazR4DVzrlF1R1Syrd1fyZ3zF/H1v1HmTKqC7+Z0IfI+jq1QKS28uvV75xbAiwp89iDFYw9t+qx5GSKix2zv97FzI8SaNaoLrOnxvHj3m29jiUiHtPuXIhJzy7gtte/5ZsdhxjXty2PXTFAJwmJCKBCDzmPvL+Z1UmHeeyKAfx8eGdNRxSR72kaRAhZsfMQ7677jlvHdmfSiC4qcxH5ARV6iCgsKuahRZvp2KIRvzi3h9dxRCQIqdBDxKvLd7N1/1F+d0kfGtXXB1KIyIlU6CEgLSufJz/Zxtk9o7iwXzuv44hIkFKhh4CZH20lt6CIhy7tp+PmIlIhFXqQW7vnMAtX7+WmMd3o0aaJ13FEJIip0INYUbHjwfc207ZZA24/r6fXcUQkyKnQg9iCVclsTMngNxP60KSBThkQkZNToQepw9kFzFy6lRHdWjFxUAev44hICFChB6knP0ngaF4hv5+oN0JFxD8q9CC0KSWD11bs4dpRXenTvpnXcUQkRKjQg0xxsePB9zbRKrI+d46L9TqOiIQQFXqQeWdtCt/uOcJ9F/WmeaN6XscRkRCiQg8imXnH+POHWxncuQVXDe3kdRwRCTGaCxdEnv5kO4ey83l56nDq1NEboSJyarSHHiQS9h9lzr+TmDyiCwM6Nfc6joiEIBV6EHDO8dCiTTRtWJdfX9DL6zgiEqJU6EFg8YZ9LN+Zzj0X9KJl4/pexxGREKVC91h2fiF//GAL/To0Y/KILl7HEZEQpjdFPfbsZ4nsz8zjuWuGEqE3QkWkCrSH7qEdqVm8+NVOrhrWiWFdW3odR0RCnArdI845Hl60mYZ1I7hvfG+v44hIGFChe+Tj+AN8tT2NO8fFEt20gddxRCQMqNA9kHesiEfej6dX26Zcd2ZXr+OISJjQm6IeeP7zHaQcyWXetFHUjdD/qSISGGqTGrbtwFGe/2IHlw7qwJlntPY6joiEEe2h1xDnHPNWJvOHxfE0rh/BbybojVARCSwVeg1Iy8rnvjc38OnWg5zdM4onrhpEu+YNvY4lImFGhV7NPt1ygPve2kBmXiEPXtKXqaNjdCVFEakWKvRqklNQyB8Wb2Heyj30ad+M16cNJrZtU69jiUgYU6FXg3XJR7hzwTqSDmVzy9ju3DUulgZ1I7yOJSJhToUeQIVFxTz32Q6e+dd22jVryLxpoxjVXTNZRKRmqNADJCktmzsXrmPtniNcPrgDv7+svz4TVERqlAq9ipxzLFiVzCOL46lbx3hm8hAmDurgdSwRqYX8OrHIzMabWYKZJZrZjHKW32Vm8Wa2wcw+NbNacT77oax8ps9dw4y3NzK4cws+uuMclbmIeKbSPXQziwCeA8YBe4FVZrbIORdfathaIM45l2NmvwBmAj+vjsDB4rOtB/n1mxvIzD3Gby/uw41nddN0RBHxlD+HXEYAic65nQBmNh+4DPi+0J1zn5UavxyYEsiQwSTvWBF//GALc5fvpne7prx68wh6t2vmdSwREb8KvSOQXOr+XmDkScbfBHxY3gIzmw5MB+jSJfQ+bi07v5AbX1nFil3pTDu7G3df0IuG9TQdUUSCQ0DfFDWzKUAcMLa85c65WcAsgLi4OBfI565uWfmF3PDyStbsPsxfJw3mssEdvY4kIvID/hR6CtC51P1Ovsd+wMzOBx4Axjrn8gMTLzhk5h1j6uyVrN+bwd8mD+Xige29jiQicgJ/ZrmsAnqaWTczqw9MAhaVHmBmQ4B/ABOdcwcDH9M7GbnHuPallWzYm8FzVw9RmYtI0Kp0D905V2hmtwFLgQhgtnNus5k9Aqx2zi0CngCaAG+YGcAe59zEasxdI47kFDDlpRUk7D/K81OGMa5vW68jiYhUyK9j6M65JcCSMo89WOr2+QHO5bn07AKmvLiCxNQsZl0bx496t/E6kojISelM0XKkZeUz5cUV7ErL5oXr4hgbG+11JBGRSqnQyzh4NI9rXlhB8uEcZk8dzlk9oryOJCLiFxV6KQcy85j8wnL2Z+Tx8tQR+sxPEQkpKnSffRm5XP3CCg5m5jHnxhEMj2nldSQRkVOiQgdSjuQyedZyDmcX8M+bRjKsa0uvI4mInLJaX+jJ6TlMfmE5mbnHmHvzSAZ3buF1JBGR01KrC333oWyufmEFWfmFvD5tFP07Nvc6kojIaau1hb4rLZvJs5aTX1jE69NG0q+DylxEQlutLPTEg1lc/cJyiood86aP0uVvRSQs1KpCLyp2fLb1IDPe3gjA/Omj6Nm2qcepREQCo1YU+pGcAhauTmbu8t0kp+fSqWUjXrlhBD3aNPE6mohIwIR1oW9KyeCf/07ivXXfkV9YzMhurbj/oj6M69uWehF+fZyqiEjICLtCLygs5sNN+5jzTRLf7jlCo3oRXDWsE9edGUOvdjq8IiLhK2wKfV9GLq+v2MO8lcmkZeXTLaoxD17SlyuHdaJ5o3pexxMRqXYhXejOOZbvTGfu8iSWbj5AsXOc17sN150Zw5geUdSpY15HFBGpMSFZ6Nn5hby9NoW5/05i24EsWkTW4+Yx3ZgyqiudW0V6HU9ExBMhV+gLVu3h0cVbOJpfSP+OzZh51UAmDupAw3oRXkcTEfFUyBV6xxaRnNenDdeNjmFI5xb4PvJORKTWC7lCH9MzijE99aETIiJlaTK2iEiYUKGLiIQJFbqISJhQoYuIhAkVuohImFChi4iECRW6iEiYUKGLiIQJc85588RmqcDu0/z2KCAtgHECTfmqRvmqLtgzKt/p6+qciy5vgWeFXhVmtto5F+d1joooX9UoX9UFe0blqx465CIiEiZU6CIiYSJUC32W1wEqoXxVo3xVF+wZla8ahOQxdBEROVGo7qGLiEgZKnQRkTAR1IVuZuPNLMHMEs1sRjnLG5jZAt/yFWYWU4PZOpvZZ2YWb2abzexX5Yw518wyzGyd7+vBmsrne/4kM9voe+7V5Sw3M3vGt/42mNnQGszWq9R6WWdmmWZ2R5kxNb7+zGy2mR00s02lHmtlZp+Y2Xbfny0r+N7rfWO2m9n1NZTtCTPb6vv3e8fMWlTwvSfdFqo548NmllLq33FCBd970td7NeZbUCpbkpmtq+B7a2QdVolzLii/gAhgB9AdqA+sB/qWGfNfwN99tycBC2owX3tgqO92U2BbOfnOBRZ7uA6TgKiTLJ8AfAgYMApY4eG/9X5KTpjwdP0B5wBDgU2lHpsJzPDdngE8Xs73tQJ2+v5s6bvdsgayXQDU9d1+vLxs/mwL1ZzxYeAeP7aBk77eqytfmeVPAg96uQ6r8hXMe+gjgETn3E7nXAEwH7iszJjLgDm+228C51kNfcioc26fc+5b3+2jwBagY008dwBdBvzTlVgOtDCz9h7kOA/Y4Zw73TOHA8Y59yWQXubh0tvZHODycr71QuAT51y6c+4w8AkwvrqzOec+ds4V+u4uBzoF8jlPVQXrzx/+vN6r7GT5fN3xM2BeoJ+3pgRzoXcEkkvd38uJhfn9GN9GnQG0rpF0pfgO9QwBVpSz+EwzW29mH5pZvxoNBg742MzWmNn0cpb7s45rwiQqfhF5uf6Oa+uc2+e7vR9oW86YYFiXN1LyG1d5KtsWqtttvsNCsys4ZBUM6+9s4IBzbnsFy71eh5UK5kIPCWbWBHgLuMM5l1lm8beUHEYYBPwNeLeG441xzg0FLgJ+aWbn1PDzV8rM6gMTgTfKWez1+juBK/ndO+jm+prZA0Ah8FoFQ7zcFp4HzgAGA/soOawRjCZz8r3zoH89BXOhpwCdS93v5Hus3DFmVhdoDhyqkXQlz1mPkjJ/zTn3dtnlzrlM51yW7/YSoJ6ZRdVUPudciu/Pg8A7lPxaW5o/67i6XQR865w7UHaB1+uvlAPHD0X5/jxYzhjP1qWZTQUuAa7x/YdzAj+2hWrjnDvgnCtyzhUDL1Tw3J5ui77+uAJYUNEYL9ehv4K50FcBPc2sm28vbhKwqMyYRcDx2QRXAf+qaIMONN/xtpeALc65pyoY0+74MX0zG0HJ+q6R/3DMrLGZNT1+m5I3zzaVGbYIuM4322UUkFHq0EJNqXCvyMv1V0bp7ex64L1yxiwFLjCzlr5DChf4HqtWZjYeuBeY6JzLqWCMP9tCdWYs/b7MTyp4bn9e79XpfGCrc25veQu9Xod+8/pd2ZN9UTILYxsl734/4HvsEUo2XoCGlPyqngisBLrXYLYxlPzqvQFY5/uaANwK3OobcxuwmZJ37JcDo2swX3ff8673ZTi+/krnM+A53/rdCMTV8L9vY0oKunmpxzxdf5T857IPOEbJcdybKHlf5lNgO/B/QCvf2DjgxVLfe6NvW0wEbqihbImUHHs+vg0en/XVAVhysm2hBtffXN/2tYGSkm5fNqPv/gmv95rI53v8lePbXamxnqzDqnzp1H8RkTARzIdcRETkFKjQRUTChApdRCRMqNBFRMKECl1EJEyo0EVEwoQKXUQkTPw/wLLGKnbVnjgAAAAASUVORK5CYII=","text/plain":["
"]},"metadata":{"needs_background":"light"},"output_type":"display_data"}],"source":["plt.plot(run_SI(largest_component, tmax=20, beta=0.05, initial_inf=0.1))"]},{"cell_type":"code","execution_count":19,"metadata":{"colab":{"base_uri":"https://localhost:8080/","height":303},"executionInfo":{"elapsed":1697,"status":"ok","timestamp":1648468087684,"user":{"displayName":"Yong Yeol Ahn","photoUrl":"https://lh3.googleusercontent.com/a/default-user=s64","userId":"00405984523192108505"},"user_tz":240},"id":"mf5smaShBltI","outputId":"a9f6dc8a-f848-40d9-9db3-42042255d08a"},"outputs":[{"data":{"text/plain":["[]"]},"execution_count":19,"metadata":{},"output_type":"execute_result"},{"data":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxU9b3/8deHkABZgISwJ+wIRGWNyNJWW0VxpXVpUXGjgljxaq3ttbWl1m4/621rucUqLhVEFsVWsdKiuLS1CpKwhzWyJWxhCWEJ2b+/P2bwxpiQgUzmzEzez8cjj5w55yTzfhxm3px8zzLmnENERCJfM68DiIhIcKjQRUSihApdRCRKqNBFRKKECl1EJEo09+qJU1NTXY8ePbx6ehGRiJSdnX3QOde+tmWeFXqPHj3Iysry6ulFRCKSme2sa5mGXEREooQKXUQkSqjQRUSiRL2FbmYvmFmBma2vY7mZ2XQzyzWztWY2NPgxRUSkPoHsob8IjD3N8iuAvv6vycCfGh5LRETOVL2F7pz7F3D4NKuMA2Y7n2VAWzPrHKyAIiISmGCMoXcF8qo9zvfP+wIzm2xmWWaWdeDAgSA8tYiInBLS89CdczOBmQCZmZm6b6+IRLyqKsfJ8krfV1klxWW+6eKyCkrKfY+Lyyo/N31J/w4MSm8b9CzBKPTdQHq1x2n+eSIiYaWqylF0spzjpRWcKKvgRKmveD/7XlZJcWkFJ0r903Wsc9K/rLisktKKqjPO0SGpRdgW+iJgqpnNBy4Eipxze4Pwe0VETquyynGkuIzDJ8o4eNz3/fCJUg6dKOOQ//GhE6WfTRcWl1EV4NhAQlwM8S2a+77HNSehRQxt4+Po0tb3OD4uhlZxMbSKjakx3ZxWcc1oFVv7Oi2bx9CsmTXK9qi30M1sHnAxkGpm+cBPgVgA59zTwGLgSiAXKAbubJSkIhL1ThV0XYX8+ekyjpymoNvGx5KSEEe7hDh6t0/kgp6+6ZSEOJJaxtZa2Ke+N2bpNqZ6C905d1M9yx1wb9ASiUjUOVlWSX5hMfmFJ8k/cpKDx0prLejC4jLq+lTM2go61V/QKYktfNOJvsfJ8XHExjS96yY9uzmXiESP0opK9hwpIb+wmLzDJ8nzl3feYd/3g8dLv/AzbeNjaZcQR7uEFvRun8jwanvQ7RJb+KYTfcuT42Np3gQL+kyp0EUkIFVVjm0Hj7M6r4hdh074Cttf3PuOlnxuz7p5M6NL21akp7Tikv4dSE9pRXpKPGnJrUhLjqddQpwKuhGo0EWkVifLKlmdd4SVuwrJ3lnIyl2FHCkuB8AMOrduSVpKPKN6p5KW/H+FnZ4ST6fWLYmJwDHoSKdCFxEA9hadJHtnIVk7fOW9Yc9RKvxHHPt0SOTyjE4M657MkG5t6d4ugbjm2sMONyp0kSaovLKKTXuPkb3zMFk7C1m5s5A9RSUAtIxtxuD0ttx9US9fgacnk5wQ53FiCYQKXaSJ2HHwBO9s2M97mwpYnXeEk+WVAHRu05Jh3ZOZ1D2ZYd2TGdC5dZM8QyQaqNBFopRzjrX5RbyzYT9vb9jHlv3HAejfKYlvXZDOMH+Bd2nbyuOkEiwqdJEoUlZRxbJth3h7wz6Wbihg39ESYpoZF/RIZtrVGYzJ6Eh6SrzXMaWRqNBFItzRknI+2HyAdzbs54NNBRwrraBVbAwXndOeMRkd+Vr/DhoDbyJU6CIRaG/RSZZu2M/bG/azbNshyisdqYlxXHl+Zy47tyOj+6TSMjbG65gSYip0kQhRXFbBG6v3MH9FHmvyjgDQMzWBiaN7MiajI0O6Jevc7yZOhS4S5rYfPMFLH+/k1ew8jpVU0L9TEt+/vB+Xn9uR3u0TMVOJi48KXSQMVVY53ttUwOyPd/DvrQeJjTGuOK8zt43szrDuySpxqZUKXSSMHD5RxoIVecxZtpPdR07SqXVLHhxzDuOHp9MhqaXX8STMqdBFwsDqvCPM/ngHf1u7l7KKKkb2asePrxrAmIyOuomVBEyFLuKRkvJK/rZ2Ly99vIM1+UUkxMXwrcx0bh3ZnXM6JnkdTyKQCl0kxPIOFzNn+U5eWZFHYXE5fTok8ti4c/nGkK4ktYz1Op5EMBW6SIjsKyrhd+9sZmF2PmbGZRkduXVkd0b2aqeDnBIUKnSRRna0pJxn/vkpz3+4naoquHN0T+76ck86t9E9VCS4VOgijaSsooq5y3cy/b1cDp8oY9zgLjx0WT/dS0UajQpdJMiccyxet4/fLNnEzkPFjOrdjh9eMYDz09p4HU2inApdJIiWbzvEr/6+iTV5R+jfKYkX77yAi85przFyCQkVukgQbN1/jMf/sYmlGwvo1LolT9wwkOuGpuneKhJSKnSRBth/tIQnl25hwYo8EuKa84Ox/Zg4uqfudCieUKGLnIXjpRXM/OenPPvv7VRUVXHHqJ5M/VofUnTfcfGQCl3kDFRWOeYu38mTS7dy6EQZ1wzqwvcv60e3djpzRbynQhcJ0PHSCu6ft4p3NxUwolcKL1wxgEHpbb2OJfIZFbpIAPILi7lrVhZbC47z2LhzuXVEd525ImFHhS5Sj5W7Cpk8O4vSiipevPMCvty3vdeRRGqlQhc5jTdW7+b7C9fSuU1L5k/OpE8H3QVRwpcKXaQWVVWOJ5duYfp7uQzvmcIzE4aRrDNYJMyp0EVqOFlWyUML1/DW2r18MzONX3z9fOKa60MmJPyp0EWqKThawqTZWazdXcQPr+jP5K/00sFPiRgqdBG/9buLuGtWFkdLypl5ayZjMjp6HUnkjKjQRYAlOft4YP5qkuNjWThlFBldWnsdSeSMBTQwaGZjzWyzmeWa2cO1LO9mZu+b2SozW2tmVwY/qkjwOef40wefMmVONud0SuL1qaNV5hKx6t1DN7MYYAYwBsgHVpjZIufchmqr/Rh4xTn3JzPLABYDPRohr0jQlFZU8qO/rOe1lflcM6gLT9wwUDfVkogWyJDLcCDXObcNwMzmA+OA6oXugFO7NW2APcEMKRJsh46XMmVONit2FPLApX25/5K+OvgpES+QQu8K5FV7nA9cWGOdR4G3zew+IAG4tLZfZGaTgckA3bp1O9OsIkGxdf8xJs5aQcHRUv73piFcM6iL15FEgiJYJ9feBLzonEsDrgReMrMv/G7n3EznXKZzLrN9e10+LaG3bNshrnvqI06WVTF/8giVuUSVQPbQdwPp1R6n+edV921gLIBz7mMzawmkAgXBCCkSDB9uPchds1eQlhzPrInD6dq2ldeRRIIqkD30FUBfM+tpZnHAeGBRjXV2AZcAmNkAoCVwIJhBRRri/c0FTJy1gh7tEpg/eYTKXKJSvYXunKsApgJLgI34zmbJMbPHzOxa/2rfAyaZ2RpgHnCHc841VmiRM/HOhv3cPTubvh0SmTdpBKmJLbyOJNIoArqwyDm3GN+piNXnTas2vQEYHdxoIg23eN1e/mveKs7t2obZdw6nTXys15FEGo2uFJWo9cbq3Tz4yhoGp7flxTsvIKmlylyim24hJ1FpYXY+DyxYTWb3ZGZPHK4ylyZBe+gSdeZ9sosf/XUdo3un8uxtmbSK09Wf0jRoD12iyuyPd/DDv6zjonPa89ztKnNpWrSHLlHjuX9v4xdvbWRMRkf+ePMQWjRXmUvTokKXqDDj/VyeWLKZK8/vxB/GDyE2Rn98StOjQpeI5pzjD+9u5cmlWxk3uAu/vXEQzVXm0kSp0CViOed4YslmnvrgU24Ylsbj1w8kppnumChNlwpdIpJzjl++tZHnPtzOTcO78cuvn0czlbk0cSp0iThVVY6fvZnDrI93cvvI7jx67bm6l7kIKnSJMFVVjkdeX8+8T3Yx6cs9+dGVA1TmIn4qdIkov/77RuZ9sot7v9qbhy7rpzIXqUanA0jEmP3xDp7993ZuH9ldZS5SCxW6RIR3N+7n0UU5XDqgA9Ou0Zi5SG1U6BL21uUXMXXuKs7t0obpNw3RqYkidVChS1jbfeQkE2etICUhjufvyCQ+Tod9ROqid4eEraMl5dz5508oKa/k5bsupENSS68jiYQ17aFLWCqvrOI7c1ay7cAJnp4wjHM6JnkdSSTsaQ9dwo5zjkf+uo4Pcw/yxA0DGd0n1etIIhFBe+gSdma8n8srWfn81yV9uTEz3es4IhFDhS5h5fVVu/mft7dw3ZCufPfSvl7HEYkoKnQJG8u2HeIHC9cyolcK/+/6gTrXXOQMqdAlLOQWHOful7JJT2nFMxMyiWuul6bImdK7Rjx38Hgpd774CbExxot3DqdNfKzXkUQiks5yEU+dLKvkrllZHDhWyvzJI0lPifc6kkjEUqGLZyqrHA8sWMWa/CM8PWEYg9Pbeh1JJKJpyEU88+vFG1mSs58fX5XB5ed28jqOSMRToYsnZn20g+c+3M4do3owcXQPr+OIRAUVuoTc0g37+dmbvlvh/uTqDJ2eKBIkKnQJqZw9Rdw3T7fCFWkMKnQJmaLicqbMyaZNq1iev123whUJNr2jJCSq/Ge07CsqYcHdI+nQWrfCFQk27aFLSPzx/Vze33yAaVdnMLRbstdxRKKSCl0a3QebC/j90i18Y0hXJozo7nUckailQpdGlXe4mAcWrKZfxyR+9Y3zdUaLSCNSoUujKSmv5Dsvr6SyyvH0hGG0iovxOpJIVAuo0M1srJltNrNcM3u4jnW+aWYbzCzHzOYGN6ZEokcX5bBudxG/++ZgeqQmeB1HJOrVe5aLmcUAM4AxQD6wwswWOec2VFunL/BDYLRzrtDMOjRWYIkMC1bsYv6KPO79am/GZHT0Oo5IkxDIHvpwINc5t805VwbMB8bVWGcSMMM5VwjgnCsIbkyJJOvyi/jJGzl8qU8qD47p53UckSYjkELvCuRVe5zvn1fdOcA5ZvYfM1tmZmNr+0VmNtnMssws68CBA2eXWMJa4Yky7nk5m9SEOP4wfrCuBBUJoWAdFG0O9AUuBm4CnjWzL9wL1Tk30zmX6ZzLbN++fZCeWsKF73a4qyk4WspTE4bRLrGF15FEmpRACn03UP2j19P886rLBxY558qdc9uBLfgKXpqQ6e9u5Z9bDjDtmgzd21zEA4EU+gqgr5n1NLM4YDywqMY6r+PbO8fMUvENwWwLYk4Jc+9vKmD6e1u5fmgat1zYzes4Ik1SvYXunKsApgJLgI3AK865HDN7zMyu9a+2BDhkZhuA94HvO+cONVZoCS+nLh7q36k1v/j6ebp4SMQjAd2cyzm3GFhcY960atMOeND/JU1ISXklU+ZkU+UcT08YqouHRDykuy3KWXPO8ZPX15Oz5yjP355J93a6eEjES7r0X87a/BV5vJqdz31f68MlA3TxkIjXVOhyVtbkHeGnb+Tw5b6pPHDpOV7HERFU6HIWDp8o4zsvr6R9Ugumj9fHyImEC42hyxmprHLcP38VB46V8uqUkSQnxHkdSUT8VOhyRqa/u5V/bz3Ir687n0G6eEgkrGjIRQL24daDTH9vK9cN7cr4C9Lr/wERCSkVugRk/9ESHliwij7tE3XxkEiY0pCL1Kuisor75q3iRGkl8yYNJT5OLxuRcKR3ptTr90u38Mn2w/z2xkH07ZjkdRwRqYOGXOS0PthcwIz3P+WbmWlcPyzN6zgichoqdKnT3qKTfHfBavp3SuJn157ndRwRqYcKXWpVXlnFfXNXUVZRxYxbdNMtkUigMXSp1f+8vZmsnYX8YfxgerdP9DqOiARAe+jyBe9u3M8z/9zGzRd2Y9zgmh8fKyLhSoUun7P7yEm+9+oaMjq3ZtrVGV7HEZEzoEKXz5RVVDF17koqKh1P3TKUlrEaNxeJJBpDl888/o9NrNp1hBk3D6VHqj6sQiTSaA9dAFiSs4/nP9zO7SO7c9XAzl7HEZGzoEIX8g4X89CraxiY1oYfXTXA6zgicpZU6E1caUUl985dCcCMm4fSornGzUUilcbQm7hfL97E2vwinp4wjPSUeK/jiEgDaA+9CVu8bi8vfrSDiaN7Mva8Tl7HEZEGUqE3UTsPneC/F65lUHpbHr6iv9dxRCQIVOhNUEl5Jd95eSXNmhkzbh5CXHO9DESigcbQm6Cf/20DOXuO8vztmaQla9xcJFpo16yJeX3Vbl5evou7v9KLSwZ09DqOiASRCr0JWb+7iP9+bS3De6bw0OX9vI4jIkGmQm8iDh4vZfLsLNolxPHULUOJjdE/vUi00Rh6E1BeWcW9L6/k0IkyFk4ZRWpiC68jiUgjUKE3Ab98ayPLtx/m998axPlpbbyOIyKNRH93R7lXsvJ48aMd3PWlnnxjiD7kWSSaqdCj2Kpdhfz4r+sZ3aedLh4SaQJU6FGq4FgJU+Zk06F1C/5401Ca6yCoSNTTGHoUKquo4p45Kzl6soLX7hlFckKc15FEJAQC2m0zs7FmttnMcs3s4dOsd72ZOTPLDF5EOVOPvplD9s5CnrhxIBldWnsdR0RCpN5CN7MYYAZwBZAB3GRmX/j0YDNLAu4Hlgc7pATu5eU7mbt8F/dc3JurB3bxOo6IhFAge+jDgVzn3DbnXBkwHxhXy3o/Bx4HSoKYT85A1o7DPLooh4v7teehy3QlqEhTE0ihdwXyqj3O98/7jJkNBdKdc2+d7heZ2WQzyzKzrAMHDpxxWKnbvqISpsxZSVpyPH8YP4SYZuZ1JBEJsQaf+mBmzYDfAd+rb13n3EznXKZzLrN9+/YNfWrxKymv5O452Zwsq2DmrcNo0yrW60gi4oFACn03kF7tcZp/3ilJwHnAB2a2AxgBLNKB0dBwzvHj19ezJu8Iv/vWYPp2TPI6koh4JJBCXwH0NbOeZhYHjAcWnVronCtyzqU653o453oAy4BrnXNZjZJYPmfWRztYmJ3P/Zf05fJz9TFyIk1ZvYXunKsApgJLgI3AK865HDN7zMyubeyAUrePPz3Ez9/ayKUDOnL/JX29jiMiHgvowiLn3GJgcY150+pY9+KGx5L65BcWc+/clfRoF8/vvzWIZjoIKtLk6XrwCHSyrJK7X8qmvKKKZ2/LJKmlDoKKiC79jzjOOR7+y1o27PV9Jmiv9oleRxKRMKE99AjzzL+28cbqPTx0WT++1l+fCSoi/0eFHkHe2bCfx/+xiasGduY7F/f2Oo6IhBkVeoTYuPco989fxcCubfjtjYMw00FQEfk8FXoEOHCslLtmZdG6ZSwzb8ukZWyM15FEJAzpoGiYKymvZMqcbA6dKOXVu0fRsXVLryOJSJhSoYcx5xw/+ss6sncW8tQtQ/UBzyJyWhpyCWN/+uen/GXVbh4ccw5Xnt/Z6zgiEuZU6GFqSc4+nliymWsGdeG+r/XxOo6IRAAVehjK2VPEdxesZmBaW564YaDOaBGRgKjQw0zBsRImzcqiTatYnr11mM5oEZGA6aBoGCkp992jpbC4nFenjKSDzmgRkTOgQg8Tzjkefm0tq3Yd4ekJQzmvq85oEZEzoyGXMPHUB5/y+uo9PHTZOYw9T2e0iMiZU6GHgX+s38sTSzYzbnAX7v2qzmgRkbOjQvfY+t1FfHfBGgant+Xx63VGi4icPRW6hwqOljBpdhbJ8bHMvE1ntIhIw+igqEdKyiuZ9FI2R4rLWXjPSDok6YwWEWkYFboHnHN8f+Fa1uQd4Zlbh3FuF53RIiINpyEXD/zxvVzeXLOHH4ztx+XndvI6johECRV6iL2+aje/fWcL1w3pyj0X6VOHRCR4VOgh9NbavTz4ympG9ErhV9edrzNaRCSoVOgh8o/1+/iv+asY1j2Z52+/QGe0iEjQqdBD4N2N+7lv3koGprXhz3cOJ6GFjkWLSPCp0BvZB5sLuGfOSgZ0bs2sicNJVJmLSCNRoTeiD7ceZPJL2fTpkMhLEy+kdctYryOJSBRToTeSZdsOcdfsFfRKTWDOXRfSJl5lLiKNS4XeCLJ2HGbiiytIS45nzl0XkpIQ53UkEWkCVOhBtmpXIXf8eQWdWrdk7l0XkprYwutIItJEqNCDaF1+Ebe98AntEuOYO2mEPnFIREJKhR4kOXuKmPD8ctq0imXupBF0aqMyF5HQUqEHweZ9x5jw3HIS4mKYN2kEXdu28jqSiDRBKvQGyi04xi3PLSOueTPmThpBekq815FEpIlSoTfAtgPHuenZ5YAxd9IIeqQmeB1JRJowFfpZ2nnoBDc/u5yqKse8SRfSu32i15FEpIkLqNDNbKyZbTazXDN7uJblD5rZBjNba2bvmln34EcNH3mHi7n52eWUVlTy8qQL6dsxyetIIiL1F7qZxQAzgCuADOAmM8uosdoqINM5NxBYCPwm2EHDRX5hMTc/t4xjJeW89O0L6d+ptdeRRESAwPbQhwO5zrltzrkyYD4wrvoKzrn3nXPF/ofLgLTgxgwPy7Yd4to//ocjxb4yP6+rPjpORMJHIIXeFcir9jjfP68u3wb+XtsCM5tsZllmlnXgwIHAU3rMOcfsj3cw4bnltI2P5fV7RzMova3XsUREPieo93I1swlAJnBRbcudczOBmQCZmZkumM/dWEorKpn2eg4LsvL4Wv8OPDl+sO6aKCJhKZBC3w2kV3uc5p/3OWZ2KfAIcJFzrjQ48bxVcLSEKXOyWbnrCPd+tTcPjulHTDN9bJyIhKdACn0F0NfMeuIr8vHAzdVXMLMhwDPAWOdcQdBTemB13hHufimLoycrmHHzUK4a2NnrSCIip1VvoTvnKsxsKrAEiAFecM7lmNljQJZzbhHwBJAIvOr/4ONdzrlrGzF3o3otO58f/nUdHZJa8No9o8joojNZRCT8BTSG7pxbDCyuMW9atelLg5zLExWVVfxq8SZe+M92RvZqx4xbhupe5iISMfQBl36FJ8qYOm8l/8k9xB2jevDIVQOIjdGFtCISOVTowKZ9R5k0O4v9RaX85oaBfDMzvf4fEhEJM02+0P++bi/fe3UNiS2aM//uEQztlux1JBGRs9JkC72qyvHk0i1Mfy+XwelteebWYXTUJwyJSARrkoV+rKSc7y5Yw9KN+7lxWBo///p5tIyN8TqWiEiDNLlC37TvKFPnrmL7wRM8ek0Gt4/qgf9USxGRiNZkCv1YSTlPLt3Kix/toE2rWF6aOJxRfVK9jiUiEjRRX+jOOd5cu5df/G0DB46XMv6Cbvzg8n4k6/xyEYkyUV3ouQXHmPZGDh99eojzurZm5m2ZDNZdEkUkSkVloReXVTD93Vye/3AbrWJj+PnXz+Pm4d10Yy0RiWpRVejOOZbk7OOxNzewp6iEG4al8fAV/UlNbOF1NBGRRhc1hb794Al+uiiHf205QP9OSUy/aQiZPVK8jiUiEjIRX+gl5ZU89X4uT/9zG3HNmzHt6gxuG9md5roPi4g0MRFd6O9u3M+jb+aQd/gk4wZ34ZErB9BBV3uKSBMVkYWed7iYn72Zw9KNBfTtkMi8SSMY2bud17FERDwVcYX+yoo8fvLGemKaGT+8oj8Tv9RTt7kVESECC71HagKXDOjAT67OoHObVl7HEREJGxFX6MN7pjC8p85eERGpSWMVIiJRQoUuIhIlVOgiIlFChS4iEiVU6CIiUUKFLiISJVToIiJRQoUuIhIlzDnnzRObHQB2nuWPpwIHgxgn2JSvYZSv4cI9o/Kdve7Oufa1LfCs0BvCzLKcc5le56iL8jWM8jVcuGdUvsahIRcRkSihQhcRiRKRWugzvQ5QD+VrGOVruHDPqHyNICLH0EVE5IsidQ9dRERqUKGLiESJsC50MxtrZpvNLNfMHq5leQszW+BfvtzMeoQwW7qZvW9mG8wsx8zur2Wdi82syMxW+7+mhSqf//l3mNk6/3Nn1bLczGy6f/utNbOhIczWr9p2WW1mR83sgRrrhHz7mdkLZlZgZuurzUsxs3fMbKv/e3IdP3u7f52tZnZ7iLI9YWab/P9+fzWztnX87GlfC42c8VEz213t3/HKOn72tO/3Rsy3oFq2HWa2uo6fDck2bBDnXFh+ATHAp0AvIA5YA2TUWOc7wNP+6fHAghDm6wwM9U8nAVtqyXcx8DcPt+EOIPU0y68E/g4YMAJY7uG/9T58F0x4uv2ArwBDgfXV5v0GeNg//TDweC0/lwJs839P9k8nhyDbZUBz//TjtWUL5LXQyBkfBR4K4DVw2vd7Y+Wrsfy3wDQvt2FDvsJ5D304kOuc2+acKwPmA+NqrDMOmOWfXghcYmYWinDOub3OuZX+6WPARqBrKJ47iMYBs53PMqCtmXX2IMclwKfOubO9cjhonHP/Ag7XmF39dTYL+HotP3o58I5z7rBzrhB4Bxjb2Nmcc2875yr8D5cBacF8zjNVx/YLRCDv9wY7XT5/d3wTmBfs5w2VcC70rkBetcf5fLEwP1vH/6IuAtqFJF01/qGeIcDyWhaPNLM1ZvZ3Mzs3pMHAAW+bWbaZTa5leSDbOBTGU/ebyMvtd0pH59xe//Q+oGMt64TDtpyI7y+u2tT3WmhsU/3DQi/UMWQVDtvvy8B+59zWOpZ7vQ3rFc6FHhHMLBF4DXjAOXe0xuKV+IYRBgH/C7we4nhfcs4NBa4A7jWzr4T4+etlZnHAtcCrtSz2evt9gfP97R125/qa2SNABfByHat4+Vr4E9AbGAzsxTesEY5u4vR752H/fgrnQt8NpFd7nOafV+s6ZtYcaAMcCkk633PG4ivzl51zf6m53Dl31Dl33D+9GIg1s9RQ5XPO7fZ/LwD+iu/P2uoC2caN7QpgpXNuf80FXm+/avafGoryfy+oZR3PtqWZ3QFcDdzi/w/nCwJ4LTQa59x+51ylc64KeLaO5/b0tejvj+uABXWt4+U2DFQ4F/oKoK+Z9fTvxY0HFtVYZxFw6myCG4D36npBB5t/vO15YKNz7nd1rNPp1Ji+mQ3Ht71D8h+OmSWYWdKpaXwHz9bXWG0RcJv/bJcRQFG1oYVQqXOvyMvtV0P119ntwBu1rLMEuMzMkv1DCpf55zUqMxsL/AC41jlXXMc6gbwWGjNj9eMy36jjuQN5vzemS4FNzrn82hZ6vQ0D5vVR2dN94TsLYwu+o9+P+Oc9hu/FC9AS35/qucAnQK8QZvsSvj+91wKr/V9XAlOAKf51pgI5+I7YLwNGhTBfL//zrh9gcu4AAACsSURBVPFnOLX9quczYIZ/+64DMkP875uAr6DbVJvn6fbD95/LXqAc3zjut/Edl3kX2AosBVL862YCz1X72Yn+12IucGeIsuXiG3s+9Ro8ddZXF2Dx6V4LIdx+L/lfX2vxlXTnmhn9j7/wfg9FPv/8F0+97qqt68k2bMiXLv0XEYkS4TzkIiIiZ0CFLiISJVToIiJRQoUuIhIlVOgiIlFChS4iEiVU6CIiUeL/A0cqciutBnmbAAAAAElFTkSuQmCC","text/plain":["
"]},"metadata":{"needs_background":"light"},"output_type":"display_data"}],"source":["# Since it is a random process we want to do a couple\n","# sample runs to smooth out the curve\n","\n","plt.plot( np.mean([run_SI(rnd_graph, tmax=20, beta=0.05, initial_inf=0.1) for i in range(50)], axis=0) )"]},{"attachments":{},"cell_type":"markdown","metadata":{"id":"71yaTNYFBltI"},"source":["The axis argument in `numpy.mean` tells you which axis to apply the average over, since we have a two-dimensional array (time on one axis and trials on the other). If I picked `axis=1` instead, it would have run the average over time rather than the number of trials.\n","\n","This curve is much smoother than the previous one. You will find that this sort of averaging over trials is necessary when dealing with noisy or random models.\n","\n","We can see that at 10% initial infected population and an infection rate of 5% we infect the whole 100 node network within 20 time steps. Most of the growth occurs in the middle after the disease ramps up, and then slows as most of the population is already infected."]},{"attachments":{},"cell_type":"markdown","metadata":{"id":"K0X_c5sIBltI"},"source":["## Building the SIS model\n","The example SI model should give you a good starting point from which to create the SIS version. In the SIS model, infected nodes can transform back to susceptible nodes, which means you will have one additional parameter that needs to be provided as an argument to the model. Lets call this probability of reversion `mu`. \n","\n","The SI model implementation is simple but far from optimal, it will be slower to run on larger and more dense graphs. If you want more of a challenge try comming up with an SIS version that can run efficiently on larger graphs. This could be done by relying more heavily on numpy by vectorizing operations and/or by using other faster libraries. \n","\n","Here are your goals:\n","\n","1. Create an SIS version of the function.\n","2. Plot your model's results using a sparse random graph and play with the parameters to get a feel for how `mu`, `beta`, and `initial_inf` change the equilibrium point of the system. The equilibrium point occurs when the system settles on a stable fraction of infected (see Fig 10.7 in Barabasi's textbook). Also take note of how long it takes for the system to reach equilibrium.\n","3. Finally, construct a graph like the figure from Barabasi's book that shows the difference between Erdos-Renyi graphs and Scale-free graphs for the SIS-model. The Y-axis in the figure will be the equilibrium point of the system. This will generally be the last time point of your simulation, assuming you run it long enough to let it reach equilibrium. The X-axis is the parameter `lambda` which is just `beta / mu`. The exact location of the critical point for the SIS model on the ER graph will vary depending upon parameters, but the key take-away is that the Scale-free graph's is lower (and eventually vanishes depending upon scaling exponent). If you are not getting this relationship, you can try running SIS on Erdos-Renyi and Scale-free networks with more nodes. Make sure to adjust the connection probability (p) in the ER graph function so that the two networks have a similar number of edges. Lastly, remember to use averaging to smooth the curves over many trials for each data-point. Note: The BA algorithm only generates exponents of 3. You can generate a directed scale-free graph with varying power-law exponent using networkx's [`scafe_free_graph`](https://networkx.github.io/documentation/stable/reference/generated/networkx.generators.directed.scale_free_graph.html?%20scale_free_graph#networkx.generators.directed.scale_free_graph) function. However, it needs to be converted to an undirected graph. You can make a power-law exponent of ~2.5 with the following parameters `alpha=0.35`, `beta=0.60`, `gamma=0.05`, `delta_in=0.4`, `delta_out=0.4`. \n","\n","4. When you are done submit your notebook to Canvas."]},{"cell_type":"code","execution_count":null,"metadata":{"id":"ahg8aNqxBltI"},"outputs":[],"source":[]},{"cell_type":"code","execution_count":null,"metadata":{"id":"MMgR40-UBltI"},"outputs":[],"source":[]}],"metadata":{"anaconda-cloud":{},"colab":{"collapsed_sections":[],"name":"m10-networkepidemiology.ipynb","provenance":[{"file_id":"https://github.com/yy/netsci-course/blob/master/m10-epidemics/epidemics_assignment.ipynb","timestamp":1648429000871}]},"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.6.7"}},"nbformat":4,"nbformat_minor":0} diff --git a/m12-diffusion/threshold_model_assignment.ipynb b/m12-diffusion/threshold_model_assignment.ipynb deleted file mode 100755 index cfdd738..0000000 --- a/m12-diffusion/threshold_model_assignment.ipynb +++ /dev/null @@ -1,315 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "\n", - " \n", - " \n", - "
\n", - " Run in Google Colab\n", - " \n", - " View on Github\n", - "
" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# The Threshold model\n", - "We will build a basic threshold model and then use it to explore the dynamical space of the model on a real network and then to compare spreading processes between real graphs and randomized versions of those graphs. In the reading we learned that simple contagions and complex contagions differ in the way they spread which can have a huge impact on the resulting dynamics of the system. This is particularly important when considering social systems.\n", - "\n", - "In Granovetter's paper we saw how adjusting the standard deviation of thresholds in a crowd induced a phase-transition where the crowd went from exhibiting no/limited rioting behavior to a full blown riot. However, the Granovetter model assumes perfect mixing in the population. Granovetter was aware of this and alluded to extensions of the model that could include friendship information and social ties. Lets try implementing a network version of the model that takes into account binary social ties between individuals." - ] - }, - { - "cell_type": "code", - "execution_count": 1, - "metadata": {}, - "outputs": [], - "source": [ - "import networkx as nx\n", - "import numpy as np\n", - "import random\n", - "%matplotlib inline\n", - "import matplotlib.pyplot as plt" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "metadata": {}, - "outputs": [], - "source": [ - "# We will take three arguments for our function. The first\n", - "# being a networkx graph, the second being the mean of the\n", - "# threshold distribution, the third being the standard\n", - "# deviation of the threshold distribution\n", - "def threshold_model(graph, mean, std):\n", - " \n", - " # First lets create our distribution of thresholds.\n", - " # Numpy has many different kinds of distributions, and\n", - " # you can easily replace 'normal' with a 'gamma' or 'beta'\n", - " # or 'uniform' or any other of your choosing, though the\n", - " # parameters will change.\n", - " threshold_distribution = np.random.normal(mean, std, size=len(graph))\n", - " # I create a dictionary with keys=nodes and values=thresholds\n", - " # that will store the threshold value for each node for future reference\n", - " node_thresholds = { node: threshold_distribution[i] for i, node in enumerate(graph.nodes()) }\n", - " \n", - " # Now lets select a seed. We will pick a random node from the graph\n", - " # along with its immediate neighbors as the set of initially active nodes\n", - " # (alternatively we could just pick one node at random as the seed)\n", - " seed_node = np.random.choice(graph.nodes())\n", - " seed_neighbors = list(graph[seed_node].keys())\n", - " seed_neighbors.append(seed_node)\n", - " seeds = set(seed_neighbors)\n", - " \n", - " # In our model we will only allow for activation, not deactivation, so we\n", - " # can cut down on processing time a bit by only looping over inactive nodes\n", - " # when updating. The inactive set is just the rest of the nodes less the seeds:\n", - " inactive_nodes = set(graph.nodes()) - seeds\n", - " \n", - " # We want as our output the fraction of nodes that are active in the graph \n", - " # when the process has reached equilibrium. We can check if the process\n", - " # is at equilibrium if between two time steps no new nodes become active.\n", - " # So we will use the 'previous_set' variable to keep track of this\n", - " # and iterate through a 'while loop' that exits when equilibrium is \n", - " # reached (since we can only activate nodes we are guaranteed to reach\n", - " # some equilibrium, even if it takes a long time)\n", - " previous_set = set([])\n", - " while (previous_set != inactive_nodes):\n", - " previous_set = set(inactive_nodes)\n", - " \n", - " # We will be using an asynchronous update scheme that randomly\n", - " # orders the inactive nodes and then updates them in that sequence.\n", - " # First turn the inactive nodes into a list and shuffle it using\n", - " # the shuffle function from python, then we iterate over that list.\n", - " update_sequence = list(inactive_nodes)\n", - " random.shuffle(update_sequence)\n", - " for node in update_sequence:\n", - " # We add up the number of active neighbors that the node has.\n", - " # You can easily extend this to work with weighted graphs by\n", - " # substituting 1 with edge weights. \n", - " num_active_neighbors = np.sum([ 0 if neighbor in inactive_nodes else 1 \n", - " for neighbor in graph.neighbors(node)])\n", - " \n", - " # Now we check if the activity from our neighbors pushes us over\n", - " # the threshold...\n", - " if (num_active_neighbors >= node_thresholds[node]):\n", - " # if it does, then we just remove that node from the inactive list\n", - " inactive_nodes.remove(node)\n", - " \n", - " # Return the fraction of active nodes at equilibrium\n", - " return 1 - len(inactive_nodes) / len(graph)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Our threshold model is complete. It is only one variant of many. For instance, with a few small modifications we could change it to use the fraction of neighbors instead of an absolute threshold.\n", - "\n", - "Lets go ahead and load in a real graph and then apply this model to it and randomized versions of the [Americal college football graph](http://www-personal.umich.edu/~mejn/netdata/) which I choose because it has a distict community structure based on the geographic location of the teams:" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Name: \n", - "Type: Graph\n", - "Number of nodes: 115\n", - "Number of edges: 613\n", - "Average degree: 10.6609\n" - ] - } - ], - "source": [ - "football = nx.read_gml('football.gml')\n", - "football = nx.Graph(football)\n", - "print(nx.info(football))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We can do two different randomizations: one in which we hold only the number of nodes and edges roughly constant, and one where we hold the degree sequence constant. We can then compare the results." - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": {}, - "outputs": [], - "source": [ - "# Lets pick a range of standard deviations and means to loop over\n", - "# You will have to adjust this range and play with it so that it fits\n", - "# your graph. I explored the dynamical space of the system\n", - "# before in order to figure out what range of values to pick,\n", - "# that part will be left to you in the homework.\n", - "stds = np.linspace(0.01, 8.0, 30)\n", - "means = np.linspace(1.0, 6.0, 30)\n", - "# And of course we need the degree sequence for the graph\n", - "deg_seq = list(dict(nx.degree(football)).values())" - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "" - ] - }, - "execution_count": 6, - "metadata": {}, - "output_type": "execute_result" - }, - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8VPW9//HXZ5bsCclkARWQQK2oCCgBr5WKWypaKrhj9V7XUutWr0uLbX+29ddrLXrv9daqLaL2tmLdF6q4L7VQtYALFdD+FFECCiEJZF9m5vP742RCCFkmk5mcyeTzfDy+jznnzJkznwnkPSffc873iKpijDEmtXjcLsAYY0z8WbgbY0wKsnA3xpgUZOFujDEpyMLdGGNSkIW7McakIAt3Y4xJQRbuxhiTgizcjTEmBfnceuOioiIdN26cW29vjDFD0po1a3aoanFf67kW7uPGjWP16tVuvb0xxgxJIvJZNOtZt4wxxqQgC3djjElBFu7GGJOCXOtzN8aYrtra2qioqKC5udntUlyXkZHB6NGj8fv9Mb3ewt0YkzQqKirIzc1l3LhxiIjb5bhGVamqqqKiooLS0tKYttFnt4yI3Cci20Xkgx6eFxH5tYh8LCJrReTwmCoxxgx7zc3NFBYWDutgBxARCgsLB/QXTDR97r8HZvfy/EnAAe1tAXB3zNUYY4a94R7sEQP9OfTZLaOqb4jIuF5WmQv8QZ379b0lIvkiso+qfjGgynrw+PUPUPnYX6nz5VDry6POl0O9P7djusGbDT38UDwaIidYT06wntyOxzpygvVkhpsQVTq/Uuh8C8JobkcYzT9GN9sRoWzRvzHj9COjeL0xxvQtHn3u+wGbO81XtC/bK9xFZAHO3j1jx46N6c3q/vwGl266p8fnQ3iooaCjpdNCPjspoIY86mJ6z8Hw5Lmb0dOe6el7yRgzSCoqKrj88stZv3494XCYOXPmcOutt5KWlrbHelu3buWqq67isccec6nS3sUj3LuLo253c1V1MbAYoKysLKY7c1+w4XfQdDtUV0NNjdM6TXurqymqqaEo8lxaGhQUQH5+7485OeBp76XqnLA9Te/94bpf1t1ruix7uPRk5lW+xOO/q+SMS/u8qtgYkyCqymmnncb3vvc9nn76aUKhEAsWLODHP/4xt956a8d6wWCQfffdN2mDHeIT7hXAmE7zo4Gtcdhu90QgK8tpo0cn7G0G1Xdnkf6LN/jgB7/npH+9nuxstwsyZnh69dVXycjI4MILLwTA6/Xy3//935SWllJaWsprr71Gc3MzDQ0N3HfffcyZM4cPPviAxsZGLrjgAj788EMOOuggNm3axJ133klZWZlrnyUe4b4MuEJEHgKOAHYlqr89VZVdcD5v/uJ2zqq7n9tuvY6f/sz6Zoy5+uqree+99+K6zalTp3L77bf3+Py6deuYNm3aHsvy8vIYO3YswWCQN998k7Vr1xIIBNi0aVPHOnfddRcFBQWsXbuWDz74gKlTp8a17lhEcyrkn4A3gQNFpEJELhaRS0Xk0vZVlgMbgY+Be4DLElZtipowYQJPF6VzMBv4yy//xubNfb/GGBN/qtrtWSqR5eXl5QQCgb2eX7FiBfPnzwdg0qRJTJ48OeG19iWas2XO6eN5BS6PW0XDVOO3vkXd/f/LhcF7uOGGo3jgAbcrMsZdve1hJ8ohhxzC448/vsey2tpaNm/ejNfrJbuHPlPt7piby2xsmSRx9Mkn8yBhzpKHeWbpTt56y+2KjBl+jj/+eBobG/nDH/4AQCgU4tprr+WCCy4gKyurx9fNnDmTRx55BID169fzj3/8Y1Dq7Y2Fe5I49thjWQKkh5q5NO9PXH01hMNuV2XM8CIiPPnkkzz66KMccMABfPWrXyUjI4Obb76519dddtllVFZWMnnyZH71q18xefJkRowYMUhVd0/c+nOirKxM7WYde5peVsaDH35IoOirFH32Dn/8I5x3nttVGTN4NmzYwEEHHeR2Gf0WCoVoa2sjIyODTz75hOOPP55//vOfe50b31/d/TxEZI2q9nkaju25J5Hyb3yDO5qaKPzsXc47aA0LF0JDg9tVGWP60tjYyMyZM5kyZQqnnnoqd99994CDfaAs3JNIeXk5fwiHCaWlsejAJWzZAosWuV2VMaYvubm5rF69mvfff5+1a9dy0kknuV2ShXsy+drXvkZrZiZrxo9nn1eW8m+nN7BoEXz+uduVGWOGGgv3JJKens7RRx/NHU1NUFfHfx7hHH1fuNDlwowxQ46Fe5IpLy/ngc8+o+0rX6HoyXu47jr405/gb39zuzJjzFBi4Z5kysvLAXhv2jR4801uOGUd++4L//7vdmqkMSZ6Fu5JZtKkSZSUlHBfMAh+P1kPLuGXv4S//x2WLnW7OmNS35dffsn8+fOZMGECBx98MCeffDL//Oc/+72dv/71rxxyyCFMnTqVLVu2cMYZZySg2p5ZuCcZj8fDCSecwJMrVqDz5sEf/sB5ZzQzfTp2aqQxCaaqnHrqqRxzzDF88sknrF+/nptvvplt27b1e1tLly7luuuu47333mO//fYb9OGBLdyTUHl5Odu2bePT44+H6mo8y57i9tth61b41a/crs6Y1PXaa6/h9/u59NJLO5ZNnTqVmTNncv311zNp0iQOPfRQHn74YQBef/11jjnmGM444wwmTpzIueeei6qyZMkSHnnkEW666SbOPfdcNm3axKRJkwDnnPizzjqLyZMnc/bZZ3PEEUeQiAs64zHkr4mzE044AYCn6uq4prQU7rmHr70yn3POgZtvhqoq+NnPoHiQ7+vR071HjEmEq6+GOI/4y9Sp0Nt4ZB988MFeQ/4CPPHEE7z33nu8//777Nixg+nTp3P00UcD8O6777Ju3Tr23XdfjjrqKFauXMkll1zCihUrmDNnDmeccYYrwwNbuCeh0aNHM3HiRF565RWuufhi+MlP4JNPuPvuCRQWwt13O/3vP/kJXHklpKfH/l7bt8PKlc7NrCI3tOr6GJmurYVDD4XTT3fawQfH7zMbk8xWrFjBOeecg9frZeTIkcyaNYtVq1aRl5fHjBkzGN1+46CpU6eyadMmZs6c2eu2vv/97wOJHR7Ywj1JlZeXs2TJElruvJP0G2+EJUsY8ctfcscdcNllcP31TrvrLucq1tNPj36vOhiEF16Ae++FP//ZmY/weiEQcO48GAhASQnM2H8bcz/7NTPX/ZZdH+bx+o1HcteNX+OLcUdyyPzJnHqWn6lTo3//hgZYuxbeeQc+/BC+9z37ojB7c2HEXw455JBu+8Z7G4MrvdPeldfrJdj5F6obgzWel/W5J6ny8nKampr422efwTe/Cb//PbS1AXDQQfDMM/Dii5CdDWeeCUcfDatW9b7NTz6BH/8Y9t8f5syBFSvgqqvgzTdh0yZnz7ytzdmb/+gjePOPH7N87KXc9ez+nPjOL8mefTT7fquMs0f9hd9wJY9vKuOHt+Sz6/Bj+G3BDdx32p9Z/fyOPU7ZrKmBV1+F226Dc891as/Lg699Da64An7zG3jwwYT9GI3pl+OOO46WlhbuueeejmWrVq2ioKCAhx9+mFAoRGVlJW+88QYzZsyI6T0Ga3hg23NPUrNmzcLr9fLyyy9z7He+4+xiP/sszJvXsU55Obz7Ltx3H/yf/wMzZjijSP7yl7tvL9vYCI8/7uyl/+Uvzj3AZ8+GO+5wAr7bsY3WrHGO3D7+OPh8cP75cN118NWvAuBVhc2b4c038bz6Joe8+DdmfnYbvieD8CR86v0KH+/zdZ5o+SYPVH6DenIBp6bDD4ezz4bDDnPa1Kmwa1eif5rGRCcy5O/VV1/NLbfcQkZGBuPGjeP222+nvr6eKVOmICIsWrSIUaNG8eGHH/b7PS677DLOP/98Jk+ezGGHHZa44YFV1ZU2bdo0Nb076qijdPr06aptbar77qt68sk9rrtrl+oNN6imp6tmZjrTl16qmpenCqoTJqj+x3+oVlT0sIFwWPXFF1WPP955QV6e6g9/qLp1a3TFNjZq7bNv6Dvzf6Vv7TNPd3oKVEGDXr/uOLxca2/+teqnn+71snHjVP/1X6N7C5P61q9f73YJCRcMBrWpqUlVVT/++GPdf//9taWlpdt1u/t5AKs1ioy1cE9iP/vZz1REtKqqSvXHP1b1eFQ//7zX13z6qer8+c6/bEaG6nnnqb72mmoo1MMLgkHVhx5SPeww50X77KO6aJHzbTEQbW2qf/mL6nXXqU6c6GwbVA85RHXhQtWVK1WDQZ0yRfWUUwb2ViZ1DIdwr62t1WnTpunkyZP10EMP1eXLl/e4roV7ilq5cqUC+uijj6pu3Oj8c/3851G99pNPVGtq+lipvl71pJOc7R54oOqSJarNzQMvvDv/7/+p/td/qR57rKrP57xnUZE+N/Lf9KzpGxPznmbIGQ7h3h8DCXc7oJrEpk+fTm5uLi+//DKUljqd7IsXQ0VFn68dPx7y83tZoaoKTjjBOW3mN7+B9evh4osHdl5lb77yFWeAnFdfhcpKeOghOPFEvrHtj5zw2b2JeU9jhjEL9yTm9/s59thjeemll5wFP/2pc/Tx8MPhtddi3/DmzfD1rztHYx99FC6/3DnSOljy852jqg88QG3mSHIatw/eexszTFi4J7ny8nI2btzIxo0b4aijnPMdi4qcve5bb3V6svtjwwZnO1u2wPPPw2mnJabwKDVklZDXbOFuTLxZuCe5yFAEHXvvEyc6Q0Sefjr84AdwxhnOCerRePttmDkTWlud8yKPOSYxRfdDU14JBcHt/f6OMsb0zsI9yR144IGMHj3a6XePyMmBhx+G//xPePppOOIIZ4+8N88/D8cd53SJrFzpnGCeBFrzSyhhO/X1bldijMPr9TJ16tSOdssttwBwzDHHcOCBBzJlyhSmT5/Oe/Ee+CbO7CKmJCcilJeX89RTTxEKhfB6vZEn4JprYNo0OOssmD4d7r/fuVy1qwcfdC5EmjQJnnsORo0a3A/Ri1DACfedOyE31+1qjIHMzMweg3vp0qWUlZVx//33c/311+/+izoJ2Z77EFBeXk5NTQ3vvPPO3k/OmuUM0jJ5shPy112352Ax//M/znX/M2fC668nVbADUFJCHnXUbmtyuxJjonbkkUeyZcsWt8vole25DwHHH3884PS7T58+fe8V9tvPCe5rr3W6alavdk41/PWvnbEITjvNGUYyI2NwC4+CZ1QJAI2fVULZWJerMUnFjTF/gaampj2G4b3hhhs4++yz91jn+eefZ16noUCSkYX7EFBSUsKUKVN4+eWX+dGPftT9SmlpzoAxRxwBCxbAhAnOwDILFjhDR0a6c3rQ0tLCtm3b9mhffvnlHvM7duzgyCOP5Morr4zbMKW+fZ1wb63YDli4G/f11i1z7rnn0tDQQCgU6v4v6SRi4T5ElJeX8+tf/5rGxkaysrJ6XvG885wumosvdkYGu/HGbsfira2tZdmyZTzyyCOsWLGCmpqabjc3YsQIRo4cyciRIxk3bhxLly5lyZIlzJo1i6uuuopTTjkFny/2/0YZY51wb9tip0OaLtwY87cPS5cuZcqUKSxcuJDLL7+cJ554wu2SemThPkSUl5dz22238cYbbzB79uzeV548udvxfxsaGnjmmWd4+OGHWb58OS0tLYwePZozzzyTMWPGMGrUqI4gj7SMLl051dXV3HvvvfzmN7/h9NNPZ+zYsVx22WVccsklFBYW9vtzZe3v3E4q/KWFuxka/H4/v/jFL5gwYQIbNmzgoIMOcruk7kUzRkEimo0t0z8NDQ2alpam11xzTb9e19jYqI899pieeeaZmpmZqYDus88+etVVV+nKlSs11OOIYr1ra2vTJ554Qo899lgFNCMjQy+55BJ9//33+1ff9jpV0Fdn/yqmOkxqSYaxZTwej06ZMqWj/fCHP1RV1VmzZumqVas61rvtttv0oosuSmgtAxlbxvbch4isrCxmzpzJU089xfjx4xERPB7PXi2yPBgM8tJLL7Fs2TIaGhooKSnhwgsv5KyzzmLmzJm7T6mMkc/n49RTT+XUU0/lH//4B3fccQcPPPBAR5fNRRddxMSJExk/fjyFhYVID7dpyijMppFMPNWVA6rHmHgJhULdLn/99df3mL/22msHoZrYWbgPIaeffjqXX345V1xxRVTrFxYW8u1vf5uzzz6bWbNmDahvvDeHHnooixcv5pZbbuHee+/lzjvv5Pzzz+94Picnh9LSUkpLSxk/fvwe0+PGjaPKU0L6TuuWMSaeRKO47ltEZgP/A3iBJap6S5fnxwL/C+S3r7NQVZf3ts2ysjJdvXp1rHUPWzU1NQSDQcLhMKpKOBzeq0WWjxs3Dr/fP+g1hkIh1q1bx6effrpX27hxI42NjR3righrvIchxSVM3frcoNdqkktS92G7oLufh4isUdWyvl7b566ciHiBO4FyoAJYJSLLVHV9p9V+AjyiqneLyMHAcmBc9B/BRKugoMDtEvrk9XqZPHlyt6dLqiqVlZUdYX/BBRewQ/MZ12B77sahqj124w0n0ex49yaaK1RnAB+r6kZVbQUeAuZ2rQPIa58eAWwdUFUmZYkIJSUlHHHEEcyfP5+ioiKqvCNsZEgDQEZGBlVVVQMOtqFOVamqqtrrbLX+iKYTdj9gc6f5CuCILuv8DHhRRK4EsoETYq7IDCuBQIAdO7MpaN7uDF9se2zD2ujRo6moqKCy0g6wZ2RkMDpyp/sYRBPu3f22df1aPQf4var+p4gcCfxRRCapaniPDYksABYAjB1rVyMaJ9wrP80gTVudoYsTcRd4M2T4/X5KS0vdLiMlRNMtUwGM6TQ/mr27XS4GHgFQ1TeBDKCo64ZUdbGqlqlqWXFxcWwVm5QSCATYHvlvuN26ZoyJl2jCfRVwgIiUikgaMB9Y1mWdz4HjAUTkIJxwt7+rTJ8CgQBfhJw/8EJfWLgbEy99hruqBoErgBeADThnxawTkZtE5JT21a4FviMi7wN/Ai7Q4X5ExEQlEAhQ0doKQNNnFu7GxEtUV7W0n7O+vMuyGztNrweOim9pZjgIBAJsDW0CoPnz7eS4W44xKcNu1mFcFQgE2EELYCNDGhNPFu7GVYFAgDbqqKbARoY0Jo4s3I2rAoEAsJPtlNjZMsbEkYW7cZUT7rvYTgneKgt3Y+LFwt24qnO4+21kSGPixsLduKpzuGfUWrgbEy8W7sZV2dnZ+P1ClaeQzMYqCAbdLsmYlGDhblwlIgQCAaq8I/CgUFXldknGpAQLd+M6J9xznRk7Y8aYuLBwN64LBALs8GQ6MxbuxsSFhbtxnTMyZPtIGBbuxsSFhbtxXSAQYFtkmDkLd2PiwsLduC4QCLC9tZEgXgt3Y+LEwt24LhAI0BaqppJiC3dj4sTC3biu84VMIRs8zJi4sHA3rus8eFjY7sZkTFxYuBvXdd5zV+uWMSYuLNyN6/YYGbLabr1rTDxYuBvXde6W8TbUQVOT2yUZM+RZuBvXdd5zB6DS9t6NGSgLd+O6vLw8ROp2h7v1uxszYBbuxnUej4f8fLFwNyaOLNxNUigs9LOdImfGwt2YAbNwN0mhsLCAahv215i4sXA3SSEQCNDsDdPiy7JwNyYOLNxNUoicDrkzrcTC3Zg4sHA3SSEQCBAK1VDttXA3Jh4s3E1ScMK9ikqxcDcmHizcTVKIdMt8EbZhf42JBwt3kxQiV6lubWsPd9U+X2OM6ZmFu0kKkXCvaB0JbW2wa5fbJRkzpFm4m6QQ6Zb5Ukc6C6xrxpgBsXA3SWGvwcMs3I0ZEAt3kxQs3I2JLwt3kxTy8/OxcDcmfqIKdxGZLSIficjHIrKwh3XOEpH1IrJORB6Mb5km1fl8PrKzg+ywwcOMiQtfXyuIiBe4EygHKoBVIrJMVdd3WucA4AbgKFWtEZGSRBVsUld+vrClwU9LToB0C3djBiSaPfcZwMequlFVW4GHgLld1vkOcKeq1gCoqv1mmn4LBLwANGXbhUzGDFQ04b4fsLnTfEX7ss6+CnxVRFaKyFsiMjteBZrho7jYD0Bdlg1BYMxARRPu0s2yrpcP+oADgGOAc4AlIpK/14ZEFojIahFZXWn3yTRdFBVlAUF22ciQxgxYNOFeAYzpND8a2NrNOk+rapuqfgp8hBP2e1DVxapapqplxcXFsdZsUlRhYQCRWqr9Fu7GDFQ04b4KOEBESkUkDZgPLOuyzlPAsQAiUoTTTbMxnoWa1BcIBFCtoZJiqKqCYNDtkowZsvoMd1UNAlcALwAbgEdUdZ2I3CQip7Sv9gJQJSLrgdeA61W1KlFFm9QUuZDpi1D76ZA7drhajzFDWZ+nQgKo6nJgeZdlN3aaVuCa9mZMTHYPHlboLNi+HUaNcrUmY4Yqu0LVJI3I4GGfN9uFTMYMlIW7SRqRPffPmmwIAmMGysLdJI3d4d7eFWOnyxoTMwt3kzQ6Dqg2FaM+n+25GzMAFu4maRQUFAA7CeNFi2wIAmMGwsLdJI309HTS0poBCAbsQiZjBsLC3SSV3NwwAC0jLNyNGQgLd5NURoxwhi1qyrFwN2YgLNxNUikocP5L1tvIkMYMiIW7SSqFhc5F07vSS6C+HhobXa7ImKHJwt0klZKSNABq/O0XMtm57sbExMLdJJWRIzMAqJT2IaGta8aYmFi4m6RSUjICaOaLUMBZYOFuTEws3E1SiQwe9mlD+428LNyNiYmFu0kqkSEINjWMcBZYuBsTEwt3k1Qi4f5lXSZkZ1u4GxMjC3eTVCLdMjt3AiV2rrsxsbJwN0klsudeV+excDdmACzcTVKJhHtDgw+KbWRIY2Jl4W6SSmZmJl5vPU1NabbnbswAWLibpCIiZGa20taWTqiwPdxV3S7LmCHHwt0knZycEADNeSUQDOIcXTXG9IeFu0k6eXnOnnpDtt0o25hYWbibpJPffnHqrnQLd2NiZeFukk4g4AWgJs3C3ZhYWbibpFNU5AegymPhbkysLNxN0hk5Mh2ALS02vowxsbJwN0lnn32yAPhiRxsEAhbuxsTAwt0knX33zQFg27YWu5DJmBhZuJukU1xcANRRWdlq4W5MjCzcTdKJjC9TVRVywt3uo2pMv1m4m6QTCfedO9X23I2JkYW7STqRcK+tFSfcq6qcYQiMMVGzcDdJJzc3F9hFfb3XCXeAHTtcrcmYocbC3SQdESE9vYnGRv/ucLeuGWP6JapwF5HZIvKRiHwsIgt7We8MEVERKYtfiWY4yshopbk53cLdmBj1Ge4i4gXuBE4CDgbOEZGDu1kvF7gKeDveRZrhJzu7jdbWLAt3Y2IUzZ77DOBjVd2oqq3AQ8Dcbtb7v8AioDmO9ZlhKjc3jGqaM6Y7WLgb00/RhPt+wOZO8xXtyzqIyGHAGFV9Jo61mWFsRPuwMrskH3w+C3dj+imacJdulnXc90xEPMB/A9f2uSGRBSKyWkRWV9qFKaYXBQXOf7udu8TOdTcmBtGEewUwptP8aGBrp/lcYBLwuohsAv4FWNbdQVVVXayqZapaVlxcHHvVJuV1DPtbFbRwNyYG0YT7KuAAESkVkTRgPrAs8qSq7lLVIlUdp6rjgLeAU1R1dUIqNsNCcXEaAFu21Fu4GxODPsNdVYPAFcALwAbgEVVdJyI3icgpiS7QDE+jRmUAsHVrg4W7MTHwRbOSqi4HlndZdmMP6x4z8LLMcLfvvtkAfPllMxQXW7gb0092hapJSvvt54zpvn17+5juDQ1OM8ZExcLdJKXRo/OB8O4DqmBD/xrTDxbuJikVFTkjQ1ZXh+0qVWNiYOFuktKIESOAXezahYW7MTGwcDdJyev14vXWU1fnsXA3JgYW7iZp+f2NNDT4nLNlwMLdmH6wcDdJKyOjmaamdMjOdpqFuzFRs3A3SSszs43WVudiJruQyZj+sXA3SSsnJ0hbm3Mxk4W7Mf1j4W6SVl6eEg7noIqFuzH9FNXwA8a4IT9fAB91dWHySkrgrbdg7Vrw+/tu0t1I1cYMHxbuJmkVFnoB+PzzWiaNGeNcoTplSt8v9HjgRz+Cm26ykDfDloW7SVqRMd03b65l0rXXwrRp0NICbW29t3ffhV/8wtmDv7Hb8e2MSXkW7iZpjRzZadjfnLEwZ050LwyH4eKL4ac/hcxMuP76BFZpTHKycDdJa9SoTAC2bm3s3ws9HliyBJqa4Ac/gIwMuPLKBFRoTPKycDdJa/ToTsP+9pfXC3/8o9ONc9VVTsB/5ztxrtCY5GWnQpqkNWZMHgCVla2xbcDvh4cegtmz4bvfhQceiGN1xiQ323M3SWvs2BEAVFeHYt9Iejo88YTTX3/++c78mWfGqcJOGhthzRp4+22nz//yy50hE4xxiYW7SVp5eX6glZoaHdiGMjNh2TI48UT49redLppvfSv27YXD8M9/Oufdv/2209auhVCnL6G77oK774aTThpY7cbEyLplTNISAY+njtraOJyrnp0Ny5fDYYfBGWfAiy9G/9rt2+HZZ53TKr/xDQgE4KCD4MIL4cEHobAQFi50vkC2bYO//hWysuDkk+Gcc5xlxgwy23M3Sc3na6S+3hufjeXlwfPPw3HHwbx58NxzMGvWnus0NMA778Df/767bdrkPOfxwOTJMH8+HHGE0yZOdJZ3VlLinGu/aJFzvv3zz8Ntt8FFF9lFVWbQiOoA/+SNUVlZma5evdqV9zZDR17eR3i9VdTUfC1+G62sdEL988/h/vuhpsYJ8VWr4IMPnG4XgHHjYMYMp02f7lxE1d9+9I8+ggUL4I034Oij4Xe/c74QjImRiKxR1bK+1rM9d5PUMjNbqatLj+9Gi4vhlVecsD3rLGdZIOCE+Lx5u8M8cgeogTjwQHjtNedL5PrrneETfvQjpxsnPc6fy5hObM/dJLXS0nfZvDmDYPCg+G98xw5YsQIOPRTGj098l8m2bfDv/w5/+pOz9754MXz964l9T5Nyot1ztwOqJqnl5YUIhXJIyE5IUZGzpz5hwuD0hY8c6RyAfe45aG52/nJ49NHEv68ZlizcTVIbMQJgBPX19W6XEj+zZzt9++PHw733ul2NSVEW7iapFRR4gTwqK6vdLiW+srPhtNNCzXN3AAAM/0lEQVTg1Vdh1y63qzEpyMLdJLWiIueY/+bNKRiAc+c6QxQ//7zblZgUZOFuklpxcRoAFRV1LleSAEce6Zy589RTbldiUpCFu0lqu4f9bXC5kgTwep1hEJYvh9YYB0czpgcW7iap7bNPFgDbtjW7XEmCzJsHtbXw+utuV2JSjIW7SWoDGtN9KDjhBGccmqefdrsSk2Is3E1Si9xqb8eOoMuVJEhmpjNa5dNPg0sXFJrUZOFukppznjvU1ITdLSSR5s2DLVuc8eCNiRMLd5PUIuGe0qeCf/ObzsFVO2vGxFFU4S4is0XkIxH5WEQWdvP8NSKyXkTWisgrIrJ//Es1w1FaGng8zdTVpfB+SGGhM8aMhbuJoz5/Y0TEC9wJnAQcDJwjIgd3We1doExVJwOPAYviXagZvvz+Jhoa/G6XkVhz58K6dfDxx25XYlJENLtDM4CPVXWjqrYCDwFzO6+gqq+pamP77FvA6PiWaYazjIxmmppSfHjcue2/UnbWjImTaMJ9P2Bzp/mK9mU9uRh4rrsnRGSBiKwWkdWVlZXRV2mGtaysVlpbM9wuI7FKS527PFm4mziJJty7Gwu123O2ROQ8oAy4tbvnVXWxqpapallxcXH0VZphLScnRDicS1NTk9ulJNa8ebBypXOnKGMGKJpwrwDGdJofDWztupKInAD8GDhFVVP0ihPjhrw8BfKprk6xkSG7mjvXucXfM8+4XYlJAdGE+yrgABEpFZE0YD6wrPMKInIY8DucYN8e/zLNcJafL8CI1A/3ww6DMWPsrBkTF32Gu6oGgSuAF4ANwCOquk5EbhKRU9pXuxXIAR4VkfdEZFkPmzOm3woLvQyLcBdxumZeegkaG/te35heRHWDbFVdDizvsuzGTtMnxLkuYzoUFaUBWWzbVuN2KYk3dy7ccQe8+KIT9MbEKIWvDDGpYuRI5zTIlBz2t6ujj4b8fOuaMQNm4W6SXmTY3y++GAZdFX6/MxzBM89AMEUHSzODwsLdJL2SEmfPffv2FB3Tvat586Cqyjkt0pgYWbibpFdQ4FxqUVk5TPZkTzzRGVTHLmgyA2DhbpJeZGTI6uphEu65uc5NPJ56ysZ4NzGzcDdJb/eY7sMo6ObNg08/hQ8+cLsSM0RZuJukl5/vPNbWDqP/rt/6lnPeu501Y2I0jH5bzFCVm+s81td73S1kMI0aBf/yL9bvbmIW1UVMxrjJ64W0NOeGHTNnzqStrY1gMNjjYzAYJDc3l8LCwj1aUVHRXssyMjKoqamhqqqKqqoqqqur93iMTNfU1KCq+Hw+fD4ffr+/Y7q7ea/X29G6zkea3+/nuOOOY968eaSlpe39wefOhYULYfNmZ1gCY/pB1KUDNmVlZbp69WpX3tsMPSNHNuPxvM4hh9y2R5h29+j1eqmrq+sI5x07dlBVVcXOnTujei+fz0dhYSGBQKDjMRAI4PF49voS6TofmQ6FQj22yPMNDQ3s3LmTkpISLrroIhYsWEBpaenuQj76CCZOdK5YveKKBP1kzVAjImtUtazP9SzczVBw6KFwwAHwxBOxbyMYDFJTU9MR9lVVVTQ3N3eEd2RvPicnB5HuRrqOr3A4zIsvvshvf/tb/vznP6OqnHjiiXz3u99lzpw5+Hw+J9zHjHHGmzEGC3eTYmbOhPR0eOUVtytJjIqKCpYsWcKSJUvYsmUL++23H5dccgnXVFaSt3ixM8Z75MiyGdYs3E1K+eY3Yds2SPX/MsFgkGeffZbf/va3vPDCCxwJrFTl9e98h4+mTaO+vp76+noaGho6pjvPNzQ0kJ6eTk5ODrm5ueTm5nY7nZOTg8/no66ujl27dlFbW0ttbW2303V1dYTDYbxeLx6Pp+OYQWS682N3xyE6H4/ovCwzM5OsrCyys7PJzs7udTrSMjMzyczMxOsdvIPrqkpbW1tHi3S/dW0AGRkZZGZmkpGR0dF8vvge2rRwNynl3HPh7beH1/2jN27cyJLFi7l60SJWq/IrILu9jfD7CaSlUZCWxgifjzyfj1yPh2yPhwYRtqnyZTjMlrY2Klpb2dTYyObWVup6eb+srCxGjBhBXl4eeXl5HdO5ubl4PB7C4XDHcYPIdNfHyPGE7g5ydw3GpqYmGhsbaWnp/7190tPTO8K+c/CnpaX1eqA7skxVO96/t8empiaCAxzjx+fz7RH2mZmZ/PznP+ecc86JaXvRhrudLWOGhBEjoKYGNm50blYUCnXfIs+1tOzdmpv3ng8GndPJPZ6+H32+vZvfv/eytDTIyIDMTKdlZe2ezsx0nvNEcRLy+PHjufmWWwhVV3PyPfdwcucn29qc1tA+UmZGBmRnO29QWwu7dnW7TU1PJxQI0FZQQDAvD292Nr7cXHy5uXiys/csPFJsZqZzypJI9w12T3f9QXm93U/7fFBQAMXFBNPSaGxspLGxkYaGho7WeT4StJH1Oodw5+nIF0hzc/NeB747f8EAe/wlkJWVRV5e3h5fGJHn0tLS8Pv9e/zV0d08QHNz816tqalpr2WDcZtRC3czJBQXQ3U1TJgQv22KOBmj6rRwePCu9k9PdzIzPx8KC50WCHQ/XXz2rYw8fC6akUk4MxvNzCKcmd3RNDOr49si8lloaUF2VOKr3o6nqhJv1XZ81dvxVjuPvppKfLXVUFNFuKWCUGsz4ZYmPC1NSHMTMogjUnqzs8kpKianqBgtKUGKi51/8JHOtBTkO9/ufn/fLfJtDD1/AUVad18+PR1IV3W+TJubd7empj3nQyHnmz09vdum6RmE/emoxxvVl/tAWbibIeGqq2D8eGfa63V+h73e7pvHs/fvVkbG3vM+X/e/y13DPhze/RdBMLi7tbXtOR9Z1trq/N731CK50NgIO3c6A0BWVTl/lVRXO3+h7PklMwL4Zj9/Yuk4tzseHcuPGy9BMmkig2YyacJLCEH3asAe815CHc1HEB/BPaYj82m0UkANxVRS0rCd4oZKij+rpIStFPM+xVSSweDfijmEhxBegvgI4yUkXnzaRiZNeBjYN78AXiCIl79++y5mLV0Ql5p7YuFuhoTCQjj//MF5r849DG4IhXaHfnW181jXW2d5F911K3W3THX3F1LXFgz6aGvLbW97ftn0NB2Z93j67uaKfKlGvji/CMPWTl+k4ZDia6ojq6ESf8POjsK0rQ2JfKu2L5Pg7kcPYef9RXc3D4go3vZ5EZz1NIQnHETCITwaxBsOIhrCEwrufk5DhLxpBH0ZBH2Z7Y+7W5svk5DfmcbjwRduxR9uwRdq6Xj0hVvwh9qn21vJN6bG/h8kShbuxiQZr3d3l8zwJUBeezOxsLFljDEmBVm4G2NMCrJwN8aYFGThbowxKcjC3RhjUpCFuzHGpCALd2OMSUEW7sYYk4JcGxVSRCqBz2J8eRGwI47lDAX2mYcH+8zDw0A+8/6q2ufIY66F+0CIyOpohrxMJfaZhwf7zMPDYHxm65YxxpgUZOFujDEpaKiG+2K3C3CBfebhwT7z8JDwzzwk+9yNMcb0bqjuuRtjjOnFkAp3EblPRLaLyAdu1zJYRGSMiLwmIhtEZJ2IfN/tmhJNRDJE5O8i8n77Z/652zUNBhHxisi7IvKM27UMBhHZJCL/EJH3RGS12/UMBhHJF5HHROTD9t/pIxP2XkOpW0ZEjgbqgT+o6iS36xkMIrIPsI+qviMiucAaYJ6qrne5tIQREQGyVbVeRPzACuD7qvqWy6UllIhcA5QBeao6x+16Ek1ENgFlqjpsznEXkf8F/qqqS0QkDchS1Z2JeK8hteeuqm8A1W7XMZhU9QtVfad9ug7YAOznblWJpY769ll/exs6eyExEJHRODdKXeJ2LSYxRCQPOBq4F0BVWxMV7DDEwn24E5FxwGHA2+5WknjtXRTvAduBl1Q11T/z7cAPgLDbhQwiBV4UkTUikti7RSeH8UAlcH9799sSEclO1JtZuA8RIpIDPA5craq1bteTaKoaUtWpwGhghoikbDeciMwBtqvqGrdrGWRHqerhwEnA5e3drqnMBxwO3K2qhwENwMJEvZmF+xDQ3u/8OLBUVZ9wu57B1P5n6+vAbJdLSaSjgFPa+6AfAo4TkQfcLSnxVHVr++N24ElghrsVJVwFUNHpr9DHcMI+ISzck1z7wcV7gQ2q+l9u1zMYRKRYRPLbpzOBE4AP3a0qcVT1BlUdrarjgPnAq6p6nstlJZSIZLefIEB718Q3gJQ+C05VvwQ2i8iB7YuOBxJ2YoQvURtOBBH5E3AMUCQiFcBPVfVed6tKuKOAfwX+0d4HDfAjVV3uYk2Jtg/wvyLixdkBeURVh8XpgcPISOBJZ98FH/Cgqj7vbkmD4kpgafuZMhuBCxP1RkPqVEhjjDHRsW4ZY4xJQRbuxhiTgizcjTEmBVm4G2NMCrJwN8aYFGThbowxKcjC3RhjUpCFuzHGpKD/D3Q0ynkF7gC6AAAAAElFTkSuQmCC", - "text/plain": [ - "
" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "# I will make three plots. The first will be the original network,\n", - "# the second will be the configuration model, and the last will be\n", - "# an Erdos-Renyi graph that has the same average degree.\n", - "\n", - "# The code is pretty compressed, but it uses an outside list\n", - "# comprehension to iterate through each mean in the means list to\n", - "# build a list of equilibrium points.\n", - "# The inside list comprehension is just used to get the average, since\n", - "# this is a random process it will give different values for each\n", - "# iteration (both in the simulation and the graph generation).\n", - "plt.plot(means, [ np.mean([threshold_model(football, mean, 0.01) \n", - " for i in range(300) ]) \n", - " for mean in means], color='black', label='Orig')\n", - "plt.plot(means, [ np.mean([threshold_model(nx.configuration_model(deg_seq), mean, 0.01) \n", - " for i in range(300) ])\n", - " for mean in means], color='blue', label='Config')\n", - "plt.plot(means, [ np.mean([threshold_model(nx.erdos_renyi_graph(115, 0.09), mean, 0.01) \n", - " for i in range(300) ]) \n", - " for mean in means], color='red', label='ER')\n", - "plt.legend(loc='upper right')" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "For these runs I kept the standard deviation low so that the threshold is effectively the same for everyone. This results in a rather sharp critical point where the system spreads very easily and then doesn't spread at all. However, you might notice that there are several plateaus as well. Think about what causes them? Why the ER graph is the only one with the last plateau?" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Lets hold the means constant and vary the standard deviation and set the `mean=4.0`:" - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "" - ] - }, - "execution_count": 7, - "metadata": {}, - "output_type": "execute_result" - }, - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd4VNXWwOHfTkgMvQSkhRIg9CqhCCoIqDQbilJEsWEv14KoVz+Va792ERVFUBEbqCgBLggaQEU6aYChhVBD6OnJrO+PnUASUibJTCZlvc9znsnMnNlnTRjW7Kyzz95GRFBKKVWxeHk6AKWUUq6nyV0ppSogTe5KKVUBaXJXSqkKSJO7UkpVQJrclVKqAio0uRtjZhpjDhtjwvN53hhj3jXGRBtjthhjLnB9mEoppYrCmZ77LGBoAc8PA4Iyt0nA9JKHpZRSqiQKTe4iEgocLWCXq4HPxfoLqGOMaeyqAJVSShVdFRe00RTYm+1+bOZjBwp6Uf369aVly5YuOLxSSlUe69evPyIiDQrbzxXJ3eTxWJ5zGhhjJmFLNzRv3px169a54PBKKVV5GGP2OLOfK0bLxALNst0PAPbntaOIfCwiwSIS3KBBoV88SimliskVyX0BcHPmqJm+wAkRKbAko5RSyr0KLcsYY+YCA4H6xphY4P8AHwAR+RAIAYYD0UAicKu7glVKKeWcQpO7iIwt5HkB7nNZREoppUpMr1BVSqkKSJO7UkpVQJrclVKqAtLkrlxKBH79FaZNg5gYT0ejVOWlyV25xPHj8O670KEDDBkC998PLVrAJZfAhx9CfLynI1SqctHkrkpk82aYNAmaNoWHHoI6dWD2bIiKgv/8B44cgXvugUaN4Mor4euvITHR01ErVfEZO5Kx9AUHB4tOP1A+paTAvHnwwQewejX4+cG4cXDvvdCzZ859RewXwFdfwdy5EBsL1avDtdfa11x2GVRxxSQYSlUSxpj1IhJc6H6a3JWzYmLgo4/gk0/g8GFo08b2yidOhHr1Cn+9wwErV8KcOfDdd7aU06AB3HEH/Otf9melVME0uSuXSEmBn36CmTPhf/8DY2DkSNtLv+wy8CpmYS8lBRYvhlmzbPt+fra889hjEBDg0regVIXibHLXmrvK05YttobepAnceCNERsK//w07d9pkfMUVxU/sAOedB1dfDT/8YNu+4QZ4/31o1QruvBOio133XpSqjDS5qzOOH7d19OBg6NbNjnIZMgSWLIFdu+CFF+wIGFdr39724KOjbe/9iy+gXTsYOxbCwlx/PKUqA03ulZwILF8O48dD48Zw332QlgbvvAP798M338Dll4N3WjJ8/LHNwuHhkJ7u8lhatrS99927bXnml1+ga1e46ipYs8blh1OqQtOaeyU3dSo8+6wdwjhuHNx+O/ToYWvrgM3+P/4Ijzxis26WqlXtjj172q5+z562C+7t7bLYjh2D996zXzRHj9ovmenTbelGqcpKT6iqQh05AoGBMHiwHaZYtWquHSIibOH911+hc2d4+207oH39eli3zt5u2AAJCXb/atXOJvyLLrJFdV/fEsd5+rQdpfPCC5CRAW++aevyJq81wJSq4JxN7oiIR7aePXuK8qzHHxfx8hKJjMz1xNGjIg8+KOLtLVK3rsj774ukpeXdSHq6SESEyOef29f07y9SrZoIiDRsKPLssyL79rkk3j17RIYMsU1fcYXI3r0uaVapcgVYJ07kWE3uldSBAyJVq4pMmJDtwfR0kQ8/FPH3t1n/nntEjhwpeuPp6SKLF4uMGCFijEiVKiJjxoisXi3icJQo7owMkWnT7PdH7dr2O6WETSpVrjib3PWEaiX18suQmmrr7YC9uig4GO6+25ZgNmywQ2f8/YveuLe3HSv5yy+wfTs88AAsWgT9+9tjzJqFJCWxZs0apkyZwvz5821PwwleXnaM/ebNNsybb4ZRo+DQoaKHqVSF5sw3gDs27bl7zt69Ir6+IrffLiIHD9peNYg0aybyzTfu6QqfOiUyfbqktW0rAhLv7S0vggSAANKvXz/5448/itRkerrI66+LnHeeSP36It995/qwlSpr0J67ys+LL9pBMM9N3G170z/+CP/3f7B1q72ayMVnKtPS0vhx2TKuXrQIv+hoBgFhtWszxRhivL2J7tULv23b6NevH6NHjybaySuYvL3tkMkNG+z4+9Gj7Yifo0ddGr5S5ZMz3wDu2LTn7hm7don4+Ig8NyZKpGlTe8J0zRq3HCssLEweeeQRadCggQDSqFEjmTx5skRFRZ0N5rHHRGrUEAHZFhQkQ/z8xMfHRx588EGJi4tz+lipqSIvvGDL+40biyxc6Ja3pJTHoSdUVV5uu02kt88GSfdvYEezbNni0vbj4+Nl2rRp0qtXLwHEx8dHRo0aJT///LOk5Tfi5uhRkalTbW0F5J+GDWWkMVK7Vi159dVXJSkpyenjb9gg0rmz/WTfcYfIiRMuemNKlRGa3NU5tm8X6e/1hyT61rb19e3bXdJuamqq/PTTTzJq1Cjx8fERQLp06SJvvfWWHD582PmGEhJE3nnHxgays2ZNGQMS2KyZfPHFF5KRkeFUM8nJIk88YQf8tGghsmJFsd6WUmWSJnd1jpcGL5NTVJe0VkF20HgJOBwOWb9+vTz00ENnyi4NGjSQhx9+WDZu3FiyQFNTRWbNEmnfXgRkr6+vTALp2727LF++3OlmVq8WadPGfsofekgkMbFkYSlVFmhyVznsee8nScZX9jfoYkfIFNP+/fvl9ddfl86dOwsgvr6+cv3118vPP/8sqampLoxY7KD2+fPF0auXCMhRLy/5AuT1Hj1ku5Mja06fFrn/fvtJb9fObacXlCo1mtzVWXPmSLrxlr+9esuR7fHFaiIjI0PGjRsnXl5eAkjfvn1l+vTpcvToURcHmweHQ2TZMkkfO1YSqlcXAckA2dGwoZx+4gmRtWvtF0EBli611R4vL5GnnxZJSXF/2Eq5gyZ3ZX30kTiMkeUMlBceP1nsZubNmyeA3HPPPbJ161YXBlhEGRlydPFi+aVXL1mTmeQFxHH++SK33CLy9df2BG0ejh8XmTjRfuq7dRPZvLl0Q1fKFTS5K3uFD8jfDUdIw1qJ+eW8QjkcDunVq5e0bt1a0tPTXRtjCURFRcn4yy+X8SA/VqsmyZlDKsXbW2TgQJH33hOJjT3ndT/9ZAcK+fiIvPyyvRhKqfJCk3tl99xzIiDxQ24QH1Jk6tTiN7Vs2TIB5KOPPnJdfC60fPly6dGjh3iB3Nq+vey56SaRDh3sxxtELrxQ5L//Fdm588xr4uJErr/ePn3ppS6b20wpt9PkXplt3mz/aSdMkOFXpIu/v8jJ4ldkZPDgwdK4cWNJTk52XYwulpGRIbNnz5amTZsKIBMnTpTTa9fa8fM9epxN9BdcIPKf/4hERYnDIfLpp3YSsvr1RX75xdPvQqnCaXKvzMaPF6lRQ/5adFRA5NVXi9/U33//LYC8/vrrrovPjRISEuSpp54SLy8vad26tazJGh6zY4ctU/XtezbRd+woMnWqbFt7Qrp1sw89/LAdJ69UWaXJvbLaudPWnB99VAYNEjn/fDscsLiuvfZaqVu3rpwsSdffA37//Xdp1qyZVKlSRV5++eWc5wr27hV5912RAQPsf4H69SX1v+/Iv+5NFrAd/W3bPBa6UgXS5F5Z3XuviK+vrPp2n4DIW28Vv6nIyEgB5JlnnnFdfKXo6NGjMnr0aAFk4MCBsjev1T3WrhUZNMj+VwgMlHWPzBH/uhlSvbrI7NmlH7NShdHkXhkdPCji5ydyxx1y8cUiTZqIFGFalnPccsstUq1atSJN4FXWOBwOmTlzplSvXl3q1q0r8+bNy2snu7hIZm0mpcsF8kjXpQIiN91UsvMVSrmas8ldp/ytSN55B1JSOH7n46xcaRe18PMrXlMxMTHMmTOHO++8k/r167s2zlJkjOHWW29l48aNtG7dmuuuu4677rqLhKx1X+1OdnGRDRvgiy/wPRnPG1suI7r15UTO2UiPHnbJWKXKE03uFcWJEzBtGlx/PWEpbQG7VnVxvfHGGwA8+uijrojO44KCgli9ejVPPPEEM2bMoGfPnmzcuDHnTl5ecNNNsG0bvPkmrY6tZ71cwGv7xjK27w7+8580MjI8E79SRWVsL7/0BQcHyzrtDrnOa6/BE0/AunV8uLYn99wDu3fbRSyKKi4ujhYtWnDjjTfy2WefuTxUT1u+fDkTJkwgLi6Ou+++G4AjR45w5MgR4uPjz9x6JyQwGfgXUAVvFnIl4S0uZ8Lca2h5YWOPvgdVeRlj1otIcKH7aXKvAJKTITAQunSB//2PBx6AWbPg5MniLar0zDPP8OKLLxIZGUn79u1dHm5ZEB8fz6RJk/jhhx+oXbs29evXx9/fH39//zM/Z90GGEOr7+dRc9lqAsSWcw63CKbBxJGYK0faP5G89I9gVTpcmtyNMUOBdwBv4BMReSXX882B2UCdzH2miEhIQW1qcnehjz6yC1svXw6XXsqgQZCYCH/9VfSmTp48SYsWLRg0aBDz5s1zfaxljMPhwMvJxLx71y4euWwC7Xe0ZSRb6ctfeCHQuDEMHw4jR8KQIVCjhpujVpWZs8m90E+1McYbmAYMAzoCY40xHXPt9m/gWxHpAYwBPih6yKpY0tNtSaZ3bxg4EICICOjUqXjNffTRRxw/fpwnn3zSdTGWYc4mdoCWgYHMCVvKvpuF/lxAY7OHu6vNYm/Li+G77+Daa8HfH66+GjZtcmPUShXOmU92byBaRHaKSCrwNXB1rn0EqJX5c21gv+tCVAX6/nvYuROmTAFjiIuDw4eLl9yTk5N58803GTJkCMHBhXYMKqWqVasya9ZMPvigE0e9hzErozvN//yGW0Yc4dSCFXD//bBqFVxwAUyYAHv2eDpkVUk5k9ybAnuz3Y/NfCy754CbjDGxQAjwgEuiUwUTgVdegfbtbW8R22uH4iX32bNnc/DgwUrTay8uYwz33HMPv//+MXXrjqRKlanM+dabDvcMZOnQN2DHDpg82X7xtm0Ljz0GR496OmxVyTiT3PM6JZe7UD8WmCUiAcBw4AtjzDltG2MmGWPWGWPWxcXFFT1aldPixbB5sx0lk1leKG5yT09P57XXXqN3795ceumlLg60YurXrx8bN67lwguXkpHRm6Skg1x+Odz2SB2WX/4KaRHbYdw4ePNNaN0aXn/dnvxWqhQ4k9xjgWbZ7gdwbtnlduBbABH5E/ADzrnyRUQ+FpFgEQlu0KBB8SJWZ73yCgQE2ASSKSICateGprn/tirE999/z86dO3nyyScxxRliU0k1atSIX3/9lQcf7M/Ro4E0bfotX30lDB4MDS5oxtjkz1j40mbSevWzvfm2beHzz9EB88rtCruEFagC7AQCAV9gM9Ap1z6LgImZP3fAJn9TULs6/UAJrV4teU0ec8klIv36Fa0ph8Mh3bp1kw4dOkhGIcvVqfx9+eWXUrVqValVq7Fceuk7MmTIbjn/fMeZ9UMe7rZcDjQLtv9uXbuKLFpkpz5Qqghw5dwy2FLLdmAH8HTmYy8AV2X+3BFYnZn4NwGXF9amJvcSuvJKkXr1ckz56HCI+PuL3Hln0ZpauHChADJr1iwXB1n5bNmyRW6++WapU6eOAFKtWg0ZNOgpueaaCOnYMV0MGXIDX0uMT6Cdx6Z1e5FnnxUJC9NEr5zibHLXi5jKo/Bwe8HSc8/B//3fmYcPHrRDrt9+Gx56yPnmLr74YmJiYoiOjsbHx8f18VZCaWlp/Pbbb8ybN48ff/yRQ4cOcd5559G//000bHgHh2M70+aPudzomMsAfsdLHNChA4weDTfcUPyxrKrCc3acu84KWR5NmCBSvbrIkSM5Hl62zP4ttnSp802tXLlSAHn33XddHKTKkp6eLitXrpR//etf0rx5cwHE29tbLrroBunTJ1wackAeq/6BxAYNFIeXl/1H7NDB9ujDwz0dvipj0Cl/K6hdu2wB91//Ouepd96x/6IHDjjf3IgRI6R+/fqSkJDguhhVvhwOh6xbt06eeuopad++vQAC3aVq1b8ERHo23SGb73rXLvBtzNlE/+9/2/MsaWmefgvKw5xN7lqWKW8eeMBON7Bzpx0pk81dd9mh1UeOODenTHh4OF26dOGFF17gmWeecVPAqiA7duxg4cKF/PLLQlasqEV6+ktAEA0bbuKJWzZye91D1FqyBEJDweGwQ6GGDLFTFF9xBTRv7um3oEqZThxWER0+bKd5HDsWZs485+mLLrLD3UNDnWvu5ptvZv78+cTExFCvXj0XB6uK6vTp0yxevJw33khizZorEKkJzKB79x+4+co+3NqsGXXWrIElSyA21r6offuziX7AAKhWzaPvQbmfy+aWUWXIp5/ai2AmTz7nKRF7ntXZ83AxMTHMnTuXO++8UxN7GVGjRg2uv/4q/vzzRg4frs3YsfF4ed1JWNh8Hpmawfn3PcuElBT+/v57e0HDm2/aL/uPPrITl9WrB5ddZuf1P3jQ029HeZgm9/JkwQI7QVge0/Du32/X63A2ub/55psAPPLII66MULlI/fqGr746n8hIb4YPrw68iMOxh2++uYw+fe+mz6238mWDBqT89JOd2mDJErjvPti3z85v06SJnUhOE32lpcm9vIiLgzVrYMSIPJ/Omnagc+fCm4qPj2fGjBmMHz+eZs2aFf4C5THt2tnv9C1b4I47fPHxmQBsZMuW6UyYsJBmzVrz7Msvs79zZ3jjDYiMtB+G//s/+5nJnug/+EATfSWiyb28WLzY1l6GD8/z6aLMKfP++++TmJjI5DzKO6ps6tIFPvwQ9u0zvPUWBAT0AOZy8uQWpk71pnnzvowZM4bVq1fjaN/eJveICFure/ZZe77mvvvsvBSXXqqJvhLQE6rlxZgx8Ntvtv6Sxxzkd9xhe3iHDxfcTEJCAi1atKBfv34sWLDAPbEqt3M4bCXm/fdh0SLBGAfe3j+SlvYW1aptokOH9nTs2JEOHTqcuW2VlESVH36Ab7+FqCjbUKdOMGiQ3QYMgLp1PfvGVKF0tExFkp4ODRrYxSDyGCUD0LcvVK0KK1YU3NR7773Hgw8+yKpVq+jfv78bglWlLToapk+HmTOF48cN552XQI0aW0lN/YNTp5YBa4ED+Pr60rZtWzp27Mgl/v4MTk6m3b59mFWr7NJdxth56AcNsr37iy/WVaXKIE3uFUloqO1Vff89XHfdOU+L2OHPN99se3L5SUtLIygoiICAAFatWuXGgJUnJCTAvHl2rZC1ayEs7Ozkk3XqJFC//i68vTdw4sQyDh78BThG8+bNufeOO5jUvTt1N260SzX++SekpkKVKvYEflav/sILoXp1j75HpdMPVCyTJ4v4+IicOJHn03v22AsZp08vuJkvv/xSAPn555/dEKQqaxIS7EWtb78tMm6cSFCQ/ZxkbS1bHpN27d4QqCq+vr4yYcIEWbNmjX3hsmUiTz0l0revvSIaRKpUsfcnTxb55ReRY8c8/RYrJfQK1Qqkc2do2BB+/TXPpxctsudZQ0PtX9J5ERG6deuGw+Fgy5YtRVo7VFUcx47BunW2Z//997BxI9SunU7r1svZuvUBEhO3ExwczH333ceNN95I1apV4dQp+OMP+P13+yH7+29IS7NlnO7d4ZJL7HbxxbZ8qNxKe+4Vxe7dttf0xhv57vL663aX+Pj8m8ma1nf27NluCFKVRw6HyO+/i4waJeLlJeLt7ZCePbdLy5Y3CiD+/v4yefJk2bVrV84XJiaKLF8u8txzIoMGiVStevbPgXbt7MR2774r8uefIklJHnlvFRnac68gPvjADmHbutUOes7DxIl25MSBA/k3M2DAAHbt2sWOHTt0Wl91jl277PmaTz6BkyehQ4cT1KjxKevXP4VIKh06dKB79+706NGDHj160L17d/z9/e2LU1PtnwOhobaH//ffcOiQfa5KFejaFXr1Ort17GgfV8WiJ1QrihEjYNs2+OeffGcD69XLnlBdtizvJv766y8uvPBC3nrrLR5++GE3BqvKu1OnYPZsePdd+5Fr1Cidzp1Dgc/YuvU3YrPmtAGaNWt2JtFn3bZo0cIuuhwba2s/Wdu6dfYSarDz3/ToAcHBdnROjx52LntN+E7R5F4RJCaCvz9MmgTvvJPnLg4H1KwJd95pF+nIy7XXXsvvv/9OTEwMNXRom3KCw2HP5bz9tu00eHnZySivvPIkzZuvZ9u2dWzatImNGzeybds2HA4HAHXq1KFLly506dKFrl270qVLFzp37kytGjXsmM3sCX/TJvsZB/Dzsz38Cy44u3XuDOed58HfQtnkbHLXr8qybMUKO1FYPlMOAOzZY/9/5Hdl6tatW/npp5945plnNLErp3l52Y/diBG2IvjllzBnDjzwQC2qVr2Uq6++lPHjYdYsSEtLJDw8nI0bN7Jp0ybCwsL48ssvOXny5Jn2WrZseSbZd7nwQrredRdBrVpRZdcu2LDh7DZ3rr0UF2xPvlMnuzVsCOefn/emM2HmSXvuZdm998Lnn0N8fL49mF9+gSuvhNWroV+/c5+//fbbmTt3LjExMdSvX9/NAauKTMQOgf/yS3uRa3y8/cPyhhvgppvsMPisyqGIEBMTQ1hYGFu2bDlzu23bNjIyMgAv/Px86dSpE926daNr165069aNbl27Uvf48ZwJf9s2O09OQkLegVWvbpN8w4YwbBjceitU4DmTtCxT3olAy5a2Hvnjj/nu9uqrMGWKHeJWp07O5/bt20dgYCB33XUX7733nnvjVZVKair873+2N//TT5CUBIGB9sLW1FSbh/PaEhOF06eF9HQvmjbdSZ0633Lo0DSOHMlZy8+e8IODg2nVqpVtIC7OzrGR17Z7tz2hCzB0qK1VjhwJFWwAgZZlyruICIiJgX//u8DdwsPtXFC5EzvAW2+9hcPh4NFHH3VTkKqy8vW1eXPkSHsS9ocfbKJfuNBWSapXP7vVq5f9vqFaNdu9//77VkRETKFevSe4884EevfexNGjf7B582Y2b97MokWLMnv50KZNG4YNG8awYcMYOHCgHX+fl1274LPP7DQdo0bZHv3EiXD77dC2LWD/qgAwzixXVp45M17SHZuOcy/EK6/YccOxsQXu1qOHyOWXn/v40aNHpUaNGjJu3Dg3BahUyWRk2Athb7jBXoANIpdcIjJnjh0en5ycLBs2bJB3331Xhg0bJn5+fgKIn5+fDB06VN555x3Zvn173o2np4ssXCipI0aII/MK2+1Nmsh/2reXpvXqSd26dWXixIkSEhIiKSkppfvGSwgd517OXXKJHXC8aVO+u2Rk2Hmd7r3XTuWd3UsvvcTTTz/Npk2b6Natm5uDVapkDh+2He6PP7bLA/v7wy232IFiWZd3JCUlERoayqJFi1i0aBHbt28HoHXr1gwbNowhQ4Zw8uRJwsLCCA8PJzw8nL1799IQuAW40xjaiJDo68umJk1Yu38/MampJFSrRtBFF9HnqqvoPXIkvk2blulhmVpzL8+OHbOXcT/xBLz4Yr67RUdDUJBdfe+2284+7nA4aNKkCT169GDRokWlELBSruFw2LnLPv7YlnrS0+0MB82bQ+PGObeMjFgiIpaxevV8fvttGUlJSQD4+vrSsWNHOnfunGNr3qwZZuVKmDEDQkORQ4cwqal5xpFaqxZVmjbFq1EjaNXKlnTatbNbq1a2LuUhWnMvz5Yssd3yAoZAQv4LdOzevZtDhw4xatQoNwWolHtkjacfMsRe5PrZZ3YZg9277UiduLjsewcAEzFmIg0aCAEBCbRsKVxwQXXat/eiXTubk7MupAXs7JYDBgBgROxfx4cOkRoby5alS4lYvpz9mzZR5+RJApKS6HDkCE3WraPaqVNn2/D2tmePsxJ+1m3TpvYNGJP/bdbPtWrZObrdSJN7WbRwof1E9ulT4G7h4fa2Y8ecj4eFhQHQpUsXd0SnVKlo2NCOBJsy5exjaWk26R84YLf9+7N+NuzfX4PoaHt5SHr62df4+5+bh1u3hoAAQ716tTG1a+Pbti3BgwYRDKSkpLB06VK+/fZbFixYwIlTp6gDBAE9q1enT926dMrIoMXGjfgvW4Z3Pr3/Ak2fDnffXcLfUME0uZc1GRn20sChQ20PoQAREfbP1Zo1cz4enpn1Ozm7WrZS5YSPDwQE2C0/aWm2p79tG2zffvZ2yRJ70VV2fn52idmmTbNv5xEQMJK77x7J8887SEqKZc+eSLZt28rWrVv5LCqKrVu3cvjwYQzQFOji40O3hg1p0rix3Ro1onGjRjRq2BA/X187tNnhOHub3/StLqTJvaz5+297dUghJRmwyT2vBbHDwsJo2bIlNXNnfaUqAR8fey4qKOjc506dsol+507Yty/ntnatrfOnpGR/hRfQHG/v5tSpM5Q6deyw4y5doFq1VESOkZp6mJMJ+/jlRDRH963nwLo/EFkEpAEQEBBA+/btc2zdGzXC/9zwXEqTe1mzcKGtyV1xRYG7pafby8Lz2i0sLExLMkrloWZN6NnTbnkRgaNHcyb9Y8fg+PFztwMHfDl+vCHHjzckMTHn/zdjBH//JGrWjMPLaw9RUZGsWrWe5OSvgJ28996/uf/++9z6XjW5lzULF9p5BOrVK3C3HTvslYC5Ky8pKSls376da665xo1BKlUxGWNr9P7+dh4zZ6Wm2nMBu3bZbedOw65d1di1qwU7d7Zg//5Lcux/5MgxF0d+Lk3uZcm+fXZc+8svF7pr1snU3Ml927ZtpKen0zmveo1Syi18fe10Ns2a2UtUcktOtpP87dxpk/8ll9R1e0ya3MuSkBB762S9Hew02NnpSBmlyh4/v7PD5EuLLqRZloSE2K9+J3rdERH2Worci9GHh4fj4+ND28x5NJRSlZMm97IiJQWWLrW9dicmNIqIyHsO97CwMNq1a4evB6+gU0p5nib3siI01E5p6kRJJjXVjt3NL7lrSUYppcm9rFi40BbmBg0qdNd//rFDIXMn95MnTxITE6PJXSmlyb3MWLjQrnTgxJJh+c0pk3Vlqo6UUUo5ldyNMUONMduMMdGZZy5hAAAdyElEQVTGmCn57HODMSbSGBNhjPnKtWFWcNu32ykenSjJgE3uXl7Qvn3Ox7OSu/bclVKFDoU0xngD04DLgFhgrTFmgYhEZtsnCHgS6C8ix4wx57sr4App4UJ7O3y4U7tHRNiJj3JPKhcWFkaNGjVo3ry5iwNUSpU3zvTcewPRIrJTRFKBr4Grc+1zJzBNRI4BiMhh14ZZwS1caAesBwY6tXtBI2U6d+6Ml5dW25Sq7JzJAk2Bvdnux2Y+ll1boK0xZrUx5i9jzFBXBVjhJSbakTJO9tpTUuwJ1dzJXUQIDw/XkoxSCnDuCtW8Bl3nXr6pCna644HYGfRXGmM6i8jxHA0ZMwmYBGjpIMtff9k5Si+91Kndt22zswLnTu4HDx4kPj5eT6YqpQDneu6xQLNs9wOA/Xns85OIpInILmAbNtnnICIfi0iwiAQ3aNCguDFXLCtX2ouW+vd3aveskTK5c7hOO6CUys6Z5L4WCDLGBBpjfIExwIJc+/wIXApgjKmPLdPsdGWgFVZoKHTrZieJdkJEhF3DI/fsAjoMUimVXaHJXUTSgfuBJUAU8K2IRBhjXjDGXJW52xIg3hgTCawAHheReHcFXWGkptqFIfOaRi4fERF2EYLzzsv5eFhYGA0bNkT/IlJKgZOzQopICBCS67Fns/0swCOZm3LW+vWQlFSkJbciIvKeZ1qnHVBKZadj5jxp5Up762RyT0qy1zrlPpmakZFBZGSkJnel1Bma3D0pNNRO8NywoVO7b91qlwHLXVbfuXMnSUlJWm9XSp2hyd1TMjJg1aoi19vh3J67jpRRSuWmyd1TwsLgxIki1dsjI6FKlXNXdQ8PD8cYQ8eOHV0cpFKqvNLk7ilZ9fYi9Nyjomxi9/HJ+XhYWBitWrWieu5lmZRSlZYmd08JDYXmzaFFC6dfEhV17kyQoCNllFLn0uTuCSI2uReh156WBjt2nLsgdnJyMv/884+eTFVK5aDJ3RP++QcOHy5SvT062q6+lLvnHhUVhcPh0J67UioHTe6eEBpqb4vQc9+61d7m7rnrAh1KqbxocveE0FBo0MCOcXdSVJS9zd1zDwsLw9fXlzZt2rgwQKVUeafJ3ROy6u0mr9mU8xYVBQEBUKNGzsfDwsLo0KEDPrmH0CilKjVN7qUtJgb27ClSvR1sWSZ3SQbQBTqUUnnS5F7aijG+XcQm99wlmWPHjhEbG6sjZZRS59DkXtpCQ6FWrbyndsxHbCycPq0nU5VSztPkXtpCQ+Gii+yKG07SkTJKqaLS5F6aDh+2mbqI9faCRsrUrl2bgIAAFwWolKooNLmXplWr7G0R6u1gvw/q1Dl3ZuCwsDA6d+6MKcKoG6VU5aDJvTSFhoKfHwQHF+llUVG2JJM9h4sI4eHhejJVKZUnTe6lKTQULrwQfH2L9LK8Jgzbt28fx48f13q7UipPmtxLy4kTsHlzkevtx47BoUPnnkzVBTqUUgXR5F5a/vgDHI5i1dvh3J571kgZLcsopfKiyb20hIbaZZT69i3Sy7JGyuTVc2/SpAn16tVzUYBKqYpEk3tpCQ21J1KLuFrS1q22RB8YmPNxXaBDKVUQTe6lISkJ1q4tcr0dbM+9bduc1zylp6cTFRWlJRmlVL40uZeGNWvsUkpFrLdD3hOGRUdHk5KSoj13pVS+NLmXhtBQO0i9f/8ivSw5GXbuzP9kqiZ3pVR+NLmXhtBQO1FY3bpFetk//9gBNnmdTPXy8qJDXnMAK6UUmtzdLy0N/vyzWPX2/CYMCwsLo02bNlStWtUFASqlKiJN7u62YQMkJhar3h4VZas5bdvmfFynHVBKFUaTu7tlLYZdzJ57ixZQrdrZxxITE4mOjtZ6u1KqQJrc3S001Ha9GzUq8kvzmlMmMjISEdHkrpQqkCZ3d3I47DS/xei1OxywbVv+C3RoWUYpVRBN7u4UHg7Hjxer3h4TY699yutkqp+fH23atHFRkEqpikiTuztl1duLeTIV8l59qWPHjngXYZk+pVTlo8ndnVauhGbN7FnRIspvwjAdKaOUcoYmd3cRsT33iy/OuYSSk7Zuhfr17ZYlPj6eAwcO6MlUpVShNLm7S3Q0HDxYrJIM5D1SRhfoUEo5S5O7u/z+u70tZnLPa8KwrOSuZRmlVGGcSu7GmKHGmG3GmGhjzJQC9rveGCPGmKKtAF0RffutrbXn7n474cgRu+V+6YoVKwgICKBJkyYuClIpVVEVmtyNMd7ANGAY0BEYa4zpmMd+NYEHgTWuDrLciY2FZcvg5puLVW/P62RqamoqS5cuZcSIEZhitKmUqlyc6bn3BqJFZKeIpAJfA1fnsd9U4DUg2YXxlU9ffGFPqN5yS7FenteEYatWreL06dMMHz7cBQEqpSo6Z5J7U2BvtvuxmY+dYYzpATQTkV8KasgYM8kYs84Ysy4uLq7IwZYLIjB7Nlx0EbRuXawmoqKgalVo3vzsYyEhIfj6+jJo0CAXBaqUqsicSe551QDkzJPGeAFvAY8W1pCIfCwiwSIS3KBBA+ejLE/WrLHzBkycWOwmoqKgXTvwyvavExISwoABA6hRo0bJY1RKVXjOJPdYoFm2+wHA/mz3awKdgd+MMbuBvsCCSntSdfZs2+0ePbrYTWzdmvNk6q5du4iKitKSjFLKac4k97VAkDEm0BjjC4wBFmQ9KSInRKS+iLQUkZbAX8BVIrLOLRGXZcnJ8PXXMGoU1KpVrCYSE2HPnpz19pCQEABGjBjhiiiVUpVAocldRNKB+4ElQBTwrYhEGGNeMMZc5e4Ay5UFC+xEYcU8kQqwfbst2+dO7m3atCEoKMgFQSqlKoMqzuwkIiFASK7Hns1n34ElD6ucmj0bAgKgBCc9c08YlpSUxPLly5k0aZILAlRKVRZ6haqrHDgAixfDhAlQghkbo6LsidSspfV+++03kpOTtd6ulCoSTe6uMmeOXWGjBCUZsCdTW7WC886z90NCQqhatSoDBgxwQZBKqcpCk7sriMCsWXDhhXYMYwlknzBMRAgJCWHw4MH4+fmVPE6lVKWhyd0VNmyAiIgS99ozMuwJ1ayTqdu3b2fnzp1aklFKFZkmd1eYNcvWUW68sUTN7NoFqalne+5ZQyA1uSulikqTe0mlpMBXX8E110CdOiVqKveEYQsXLqRTp060KMZKTkqpyk2Te0ktXAhHj5ZouoEsWROGtW8Pp06dIjQ0VHvtSqli0eReUrNnQ+PGcNllJW4qKgoaNoS6deHXX38lLS1Nk7tSqlg0uZfE4cMQElLise1ZoqLOlmRCQkKoWbMm/fv3L3G7SqnKR5N7SXz1FaSnl3iUDNjRlFkThmUNgbz88svx8fFxQaBKqcpGk3tJzJoFvXpBx3MWpiqyQ4fstDQdOti1Uvft26clGaVUsWlyL65Nm2DzZpf02iHn6ktZQyCHDRvmkraVUpWPJvfimj0bfH1hzBiXNJd9wrCQkBAuuOACGjdu7JK2lVKVjyb34khLs3PJXHkl+Pu7pMmoKKheHapXP8Yff/yhJRmlVIloci+ORYsgLs4lY9uzZJ1MXbr0f2RkZGhyV0qViCb34pg9G84/H664wmVNZg2DDAkJoV69evTu3dtlbSulKh9N7kUVHw8//ww33QQuGqZ46hTExkL79g4WLVrE0KFD8XbBuHmlVOWlyb2o5s61NXcXjZIB2LbN3vr47CQuLk5LMkqpEtPkXhQi8Nln0KMHdO3qsmazRsrExCzBGMPQoUNd1rZSqnLS5F4UixfbudtdvJ7p1q1QpQqsWTOHvn374u+iEThKqcrLqQWyFXYljSeegNat4bbbXNp0VBS0bJnO+vV/8cILL7i0baXKk7S0NGJjY0lOTvZ0KB7n5+dHQEBAsacg0eTurDlzICwMvv7aXrzkQlFRULv2fkRE6+2qUouNjaVmzZq0bNkSY4ynw/EYESE+Pp7Y2FgCAwOL1YaWZZyRnAzPPAPBwTB6tEubPn4coqMhKWkdjRo1onv37i5tX6nyJDk5GX9//0qd2AGMMfj7+5foLxjtuTtj2jSIibEnU71c+304f76dWDIm5n1Gjx6Gl4vbV6q8qeyJPUtJfw+aSQpz7Bi8+CIMHQqDBrm8+blzoWnTJE6fXsGIESNc3r5SqmhiY2O5+uqrCQoKonXr1jz00EOkpqaes9/+/fu5/vrrPRChczS5F+bVV23t5JVXXN70oUOwfDm0aPEnVapUYciQIS4/hlLKeSLCqFGjuOaaa/jnn3/Yvn07p0+f5umnn86xX3p6Ok2aNOH777/3UKSF07JMQfbuhXfegfHjoVs3lzf/3XfgcMCRI+9z0UUXUbt2bZcfQynlvOXLl+Pn58ett94KgLe3N2+99RaBgYEEBgayYsUKkpOTSUhIYObMmYwcOZLw8HASExOZOHEiW7dupUOHDuzevZtp06YRHBzssfeiyb0gzz1ns+/UqW5pfu5cCAw8xfbtP/DAA++55RhKlVcPP/wwmzZtcmmb3bt35+233873+YiICHr27JnjsVq1atG8eXPS09P5888/2bJlC/Xq1WP37t1n9vnggw+oW7cuW7ZsITw8vEwMjNCyTH4iIuxKS/fdBy1burz5PXvgjz8gMXEmgYGB3HnnnS4/hlKqaEQkzxOZWY9fdtll1KtX75znV61axZjMtR06d+5MVxdewV5c2nPPz5NPQo0akKvW5irffGNvDx16h+++e43zzjvPLcdRqrwqqIftLp06dWLevHk5Hjt58iR79+7F29ub6tWr5/k6ESmN8IpEe+55WbnSzvw4ZYrLFuPIbc6cDHx81tO/fxOuu+46txxDKVU0gwcPJjExkc8//xyAjIwMHn30USZOnEi1atXyfd1FF13Et99+C0BkZCRhYWGlEm9BNLnnJgKTJ0OTJvDQQ245xNatsGWLN2lpn/Pmm2/quF6lyghjDD/88APfffcdQUFBtG3bFj8/P1566aUCX3fvvfcSFxdH165defXVV+natavHB0hoWSa3H36Av/6CGTOggG/qkvjooxNATa69NkMX5VCqjGnWrBk///zzOY9PnDiRidlWX2vZsiXh4eGAnQfmyy+/xM/Pjx07djB48GBatGhRWiHnSZN7dmlpttbeoYNLl9DLTgQ+/TQRL6+NvP32ZLccQylVuhITE7n00ktJS0tDRJg+fTq+Lp6Dqqg0uWc3cyZs3w4//mjn4HWDL74I59SpzgwfvobmzZu75RhKqdJVs2ZN1q1b5+kwctCae5aEBDuuvX9/uOoqtxxCRHjqqc1AGtOmDXbLMZRSCrTnftZbb8HBgzBvHrjpBOd3381j376L6dLlAC1baq9dKeU+TvXcjTFDjTHbjDHRxpgpeTz/iDEm0hizxRjzqzHGs2cSiiouDl57Da65Bvr1c8shUlJSePjhb4HmPP54gFuOoZRSWQpN7sYYb2AaMAzoCIw1xnTMtdtGIFhEugLfA6+5OlC3SU+3KyslJkIhw51K4r333uPAgQH4+mZwzTVaDVNKuZczWaY3EC0iO0UkFfgauDr7DiKyQkQSM+/+BZSPrqkI3Hsv/PILvP++HSXjBnFxcbzwwkv4+o7n6qu9qVnTLYdRSrnAwYMHGTNmDK1bt6Zjx44MHz6c7du3F7mdlStX0qlTJ7p3786+fftKfXpgZ5J7U2BvtvuxmY/l53ZgUV5PGGMmGWPWGWPWxcXFOR+lu0ydasezP/003H232w7z/PPPk5DQh9TUOowd67bDKKVKSES49tprGThwIDt27CAyMpKXXnqJQ4cOFbmtOXPm8Nhjj7Fp0yaaNm1a+tMDi0iBGzAa+CTb/QnAe/nsexO2535eYe327NlTPOrTT0VA5JZbRBwOtx0mMjJSvL29pX37P6RWLZGkJLcdSqlyLzIy0qPH//XXX+Xiiy8+53GHwyGPPfaYdOrUSTp37ixff/21iIisWLFCBgwYINddd520a9dOxo0bJw6HQ2bMmCF169aVli1byrhx42TXrl3SqVMnERFJSEiQ0aNHS5cuXeSGG26Q3r17y9q1a/OMJ6/fB7BOCsmvIuLUaJlYoFm2+wHA/tw7GWOGAE8DA0QkpQTfN+4XEgKTJsEVV9ieuxsv/3/88cepVq0u+/b1YdQo8PNz26GUqlAefhhcPOMv3btDQfORhYeHnzPlL8D8+fPZtGkTmzdv5siRI/Tq1YtLLrkEgI0bNxIREUGTJk3o378/q1ev5o477mDVqlWMHDmS66+/3iPTAztTllkLBBljAo0xvsAYYEH2HYwxPYCPgKtE5LDrw3ShtWvtItfdutnVMnx83HaopUuXsnDhQq699mNOnfLSkoxS5dSqVasYO3Ys3t7eNGzYkAEDBrB27VoAevfuTUBAAF5eXnTv3j1HIs+vrdKYHrjQnruIpBtj7geWAN7ATBGJMMa8gP3zYAHwOlAD+C5zEqwYEXHPlUAlsWMHjBgB558PCxfizjObWbPJBQYGkpBwFfXru2UJVqUqLA/M+EunTp3yrI1LAVP6Zp+u29vbm/T09AKPUVBbruTUmDwRCRGRtiLSWkRezHzs2czEjogMEZGGItI9c3NfYo+NtatcFPUXFBdnF7l2OGDxYmjUyD3xYf/xnn/+ecLCwnj++f8SEuLN6NFu/SNBKeUCgwYNIiUlhRkzZpx5bO3atdStW5dvvvmGjIwM4uLiCA0NLfakf6U1PXD5G3A9Y4adIiAoCJ5/3vbGC5OQACNHwr59dthju3ZuCy89PZ277rqLqVOnMn78eLy8riUpCS3JKFUOZE35u3TpUlq3bk2nTp147rnnGDduHF27dqVbt24MGjSI1157jUbF7CCW2vTAzpx1dcdW7NEyJ06IfPaZyKBBIsbYES/9+olMny4SH3/u/mlpIiNGiHh5ifz0U/GO6aSTJ0/KsGHDBJAnn3xSMjIyZORIkYAAkYwMtx5aqQrB06NlSkN6erokZQ6bi46OlhYtWkhKSkqe+5ZktEz567nXqmWn4/31V4iJgVdegRMn4J57oHFjGDXKzsmeknL2IqWFC+GDD9w2IRjA/v37GTBgAEuWLOHDDz/kpZde4vhxL5YsgRtvBK/y95tWSrlBYmIiF110Ed26dePaa6912/TA5XvisIAAeOIJu3LSpk3wxRfw1Vc2udetC8HBsHQp/PvfcNddbgsjIiKC4cOHEx8fz88//8zw4cMBmD/fThGvJRmlVJbSmh64YvQnjYEePeDNN+0J10WLYNgwWLUKbr8dXnjBbYdesWIF/fv3JzU1ldDQ0DOJHWDuXGjTBi64wG2HV0qpPJXvnnteqlSxo2KGDrWTgnl7u+0ipS+//JLbbruNoKAgQkJCciyrdeAArFhh/2jQJVKVUqWtYvTc81Olilsyq4jw4osvMmHChDNXpGVP7ImJ8NhjtuSfea2CUkqVqorXc3eztLQ07r33Xj755BPGjx/Pp59+muMihrAwm9AjI+GZZ6Bj7smRlVKqFGhyd8KpU6dYu3Yta9asYcGCBfz11188/fTTTJ06lcwrchGBDz+ERx6B2rXhf/+Dyy7zcOBKqSLz9vamS5cuZ+6PGTOGKVOmMHDgQA4cOICfnx++vr7MmDHDbfPCuIIm91wyMjKIiIhgzZo1Z7aIiIgzlwwHBQUxc+ZMbr311jOvOXYM7rjDjo654gqYPRsaNvTUO1BKlUTVqlXZlM+MZXPmzCE4OJjPPvuMxx9/nKVLl5ZydM6r8MldREhMTOTUqVOcPn2a06dP5/nzvn37+Pvvv1m7di0JCQkA1KtXjz59+nD99dfTp08fevfuTb169XK0v3o1jBsH+/fD66/bnruOaVeqYrvwwgt5/fXXPR1Ggcpdcl+3bh0rVqw4k5Tzus3+c0JCglMT9fj4+NC9e3duu+02+vTpQ58+fWjduvWZsktuGRnw8svw3HPQooWd7qZXLxe/WaUqM0/M+QskJSXlKLc8+eST3HjjjTn2Wbx4Mddcc41rY3Oxcpfcf//9dyZPngxAjRo1qFmz5pnbmjVr0rhxY9q2bZvjsRo1apyzb+7HatWqRZUqzv069u2DCRPsUMdx42D6dHvhrFKq/CuoLDN+/HgSEhLIyMhgw4YNpRxZ0ZS75H7vvfdy1113Ua1aNbxKsf4hAocOwfLl8OCDkJQEn30Gt9yi49iVcgtPzPlbiDlz5tCtWzemTJnCfffdx/z58z0dUr7KXXKvWrWq24+RkgJRUbB5M2zZcvY2a9nXbt3g66+hfXu3h6KUKmN8fHz4z3/+Q+vWrYmKiqJDhw6eDilP5S65F0dGhp319/Tps7fZt4QEOHjQjlHfvBm2brWvAbssXufOcOWVNql37Qr9+oEb5vlRSpUBuWvuQ4cO5ZVXXsmxT9WqVXn00Uf573//y6efflraITql3CX3Tz+F//7XrrmRtWVk5H8/Kcluzmje3Cbvq68+m8iDguwMBkqpyiEjq2eXy2+//Zbj/qOPPloK0RRfuUvuDRpAly424Xp5nd1y38/aqlWD6tWhRo2cW+7H6ta1Fx8ppVRFUO6S+1VXuXVadqWUqhD0chullKqANLkrpcoUZy46rAxK+nvQ5K6UKjP8/PyIj4+v9AleRIiPj8fPz6/YbZS7mrtSquIKCAggNjaWuKyLSioxPz8/AgICiv16Te5KqTLDx8eHwMBAT4dRIWhZRimlKiBN7kopVQFpcldKqQrIeOqstDEmDthTzJfXB464MBxX0biKRuMqGo2r6MpqbCWJq4WINChsJ48l95IwxqwTkWBPx5GbxlU0GlfRaFxFV1ZjK424tCyjlFIVkCZ3pZSqgMprcv/Y0wHkQ+MqGo2raDSuoiursbk9rnJZc1dKKVWw8tpzV0opVYByldyNMUONMduMMdHGmCmejieLMWamMeawMSbc07FkZ4xpZoxZYYyJMsZEGGMe8nRMAMYYP2PM38aYzZlxPe/pmLIzxngbYzYaY37xdCxZjDG7jTFhxphNxph1no4nizGmjjHme2PM1szP2YVlIKZ2mb+nrO2kMeZhT8cFYIz5V+ZnPtwYM9cYU/yZwQo7VnkpyxhjvIHtwGVALLAWGCsikR4NDDDGXAKcBj4Xkc6ejieLMaYx0FhENhhjagLrgWs8/TszxhiguoicNsb4AKuAh0TkL0/GlcUY8wgQDNQSkZGejgdscgeCRaRMjdk2xswGVorIJ8YYX6CaiBz3dFxZMvPGPqCPiBT3uhpXxdIU+1nvKCJJxphvgRARmeWO45WnnntvIFpEdopIKvA1cLWHYwJAREKBo56OIzcROSAiGzJ/PgVEAU09GxWIdTrzrk/mViZ6GcaYAGAE8ImnYynrjDG1gEuATwFEJLUsJfZMg4Ednk7s2VQBqhpjqgDVgP3uOlB5Su5Ngb3Z7sdSBhJVeWGMaQn0ANZ4NhIrs/SxCTgMLBWRMhEX8DYwGXB4OpBcBPifMWa9MWaSp4PJ1AqIAz7LLGN9Yoyp7umgchkDzPV0EAAisg/4LxADHABOiMj/3HW88pTcTR6PlYneXllnjKkBzAMeFpGTno4HQEQyRKQ7EAD0NsZ4vJxljBkJHBaR9Z6OJQ/9ReQCYBhwX2Yp0NOqABcA00WkB5AAlKVzYb7AVcB3no4FwBhTF1ttCASaANWNMTe563jlKbnHAs2y3Q/AjX/SVBSZNe15wBwRme/peHLL/DP+N2Coh0MB6A9clVnf/hoYZIz50rMhWSKyP/P2MPADtkzpabFAbLa/ur7HJvuyYhiwQUQOeTqQTEOAXSISJyJpwHygn7sOVp6S+1ogyBgTmPmNPAZY4OGYyrTME5efAlEi8qan48lijGlgjKmT+XNV7Id+q2ejAhF5UkQCRKQl9vO1XETc1rNyljGmeuYJcTLLHpcDHh+ZJSIHgb3GmHaZDw0GPD7AIZuxlJGSTKYYoK8xplrm/83B2PNgblFuVmISkXRjzP3AEsAbmCkiER4OCwBjzFxgIFDfGBML/J+IfOrZqADbE50AhGXWtwGeEpEQD8YE0BiYnTmSwQv4VkTKzLDDMqgh8IPNB1QBvhKRxZ4N6YwHgDmZHa6dwK0ejgcAY0w17Mi6uzwdSxYRWWOM+R7YAKQDG3HjlarlZiikUkop55WnsoxSSiknaXJXSqkKSJO7UkpVQJrclVKqAtLkrpRSFZAmd6WUqoA0uSulVAWkyV0ppSqg/wcquBbGU8FFQQAAAABJRU5ErkJggg==", - "text/plain": [ - "
" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "plt.plot(stds, [ np.mean([threshold_model(football, 4, std) \n", - " for i in range(300) ]) \n", - " for std in stds], color='black', label='Orig')\n", - "plt.plot(stds, [ np.mean([threshold_model(nx.configuration_model(deg_seq), 4, std) \n", - " for i in range(300) ])\n", - " for std in stds], color='blue', label='Config')\n", - "plt.plot(stds, [ np.mean([threshold_model(nx.erdos_renyi_graph(115, 0.09), 4, std) \n", - " for i in range(300) ]) \n", - " for std in stds], color='red', label='ER')\n", - "plt.legend(loc=\"lower right\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Here we can see a similar transition as the one in Granovetter's paper, which is fascinating! It means that at least approximately, some of the results from the well-mixed case appear in our network version of the model. While spreading results in the same outcomes for larger stardard deviations, below 2 we can see that the original graph and configuration model are notably worse than the ER version of the graph, suggesting that some aspect of the degree distribution in the football network is deterimental to spreading of complex contagions. Additionally, we can see that there is a cross-over point around 1 below which the real graph is better for spreading than the configuration model, and above which it is worse, suggesting that there are finer details in the real graph beyond degree distribution that effect the spreading process. How significant these differences are would require a more thorough analysis and perhaps some hypothesis testing." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Assignment\n", - "1. Pick a real-world graph and load it into networkx.\n", - "2. Explore the parameters of the model using just the real-world graph. Find a range of means and standard deviations that captures the full range of dynamics the model has on the network. In general it is best to start with a wide range of values and then hone in on regions of interest. Though you only have to vary one parameter at a time, matplotlib also has a [contour](http://matplotlib.org/api/pyplot_api.html?highlight=contour#matplotlib.pyplot.contour) function that will let you vary both simultaneously while plotting the equilibrium points along a z-axis. You can do this if you are feeling adventurous (though I recommend googling a tutorial first since it works differently than a normal plot function). If making a contour is too difficult, but you want to see how one parameters varies as another does, you can always make overlapping plots, similar to what I have made above but using different parameter values for each plot.\n", - "3. Once you have a good grasp on where the action is, compare your real graph with a configuration model and a comparable ER graph like I do above. Use the intuition you built up in step 2 for picking values of the mean to hold constant at, as different means can place you in different dynamical regimes resulting in very different behavior.\n", - "4. Provide a written summary of what you found in steps 2 and 3 in a markdown section below. Make note of any interesting spreading behavior your graph has and try to explain it if you can. Describe any similarities or differences between the real graph and the random graphs and consider why those similarities or differences exist.\n", - "5. Submit your completed notebook to Canvas." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.6.7" - } - }, - "nbformat": 4, - "nbformat_minor": 1 -}