Economic_Networks_PS.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "56128938",
   "metadata": {},
   "source": [
    "Based on the sheets of the OII Course: Statistical Analysis of Networks \\\n",
    "**Professor**: Prof. Renaud Lambiotte\n",
    "\n",
    "**Author**: John Pougué-Biyong \\\n",
    "Contact: john.pougue-biyong@maths.ox.ac.uk\n",
    "\n",
    "**Tutor**: Nicola Pedreschi\n",
    "\n",
    "*25 June 2024*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1122bcf2",
   "metadata": {},
   "source": [
    "# 1. Introduction to networks and basic metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3832e92e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#!pip install scipy==1.8.1\n",
    "import scipy\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ef972af",
   "metadata": {},
   "source": [
    "A common library for network analysis in Python is [NetworkX](https://networkx.org/). We first [install the library ](https://networkx.org/documentation/stable/install.html) and import it into our notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b89ed192",
   "metadata": {},
   "outputs": [],
   "source": [
    "#!pip install networkx --upgrade\n",
    "import networkx as nx"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "260eeefc",
   "metadata": {},
   "source": [
    "## Basics of NetworkX\n",
    "NetworkX enables to easily create a network, add/remove nodes & edges, and subsequently analyse and draw the network. The library documentation can be found [here](https://networkx.org/documentation/stable/tutorial.html). Let's do it! We first create an undirected network."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69181daf",
   "metadata": {},
   "outputs": [],
   "source": [
    "G = nx.Graph()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5bd498e",
   "metadata": {},
   "source": [
    "We add 5 edges"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b3286aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "list_of_edges = [(1, 2), (1, 3), (1, 4), (2, 3), (4, 2)]\n",
    "G.add_edges_from(list_of_edges)\n",
    "print(G.edges) # display the edges in G"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "854be952",
   "metadata": {},
   "source": [
    "We can also add a node as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8b2754b",
   "metadata": {},
   "outputs": [],
   "source": [
    "G.add_node(5)\n",
    "print(G.nodes) # display the nodes in G"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa8f5c94",
   "metadata": {},
   "source": [
    "The method `print` provides a short summary of the network..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26e17201",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(G)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de11213e",
   "metadata": {},
   "outputs": [],
   "source": [
    "pip install decorator==5.0.9\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd9d90a7",
   "metadata": {},
   "source": [
    "... And the `draw` method displays the network."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1966e208",
   "metadata": {},
   "outputs": [],
   "source": [
    "nx.draw(G, with_labels=True, font_weight='bold', node_color='grey')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d28e20f7",
   "metadata": {},
   "source": [
    "We will see later how to [enrich the network with node/edge attributes](https://networkx.org/documentation/stable/tutorial.html#adding-attributes-to-graphs-nodes-and-edges), and [display the network with all this information](https://networkx.org/documentation/stable/tutorial.html#drawing-graphs)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33b59e83",
   "metadata": {},
   "source": [
    "## Network analysis with NetworkX\n",
    "Once a network is built, we can analysis it by calling the plethora of [algorithms available in the NetworkX library](https://networkx.org/documentation/stable/reference/algorithms/index.html). As a matter of example, let's analyse the [Zachary's karate club network](https://en.wikipedia.org/wiki/Zachary%27s_karate_club)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e775c3d",
   "metadata": {},
   "outputs": [],
   "source": [
    "K = nx.karate_club_graph()\n",
    "print(K)\n",
    "nx.draw(K, with_labels=True, font_weight='bold', node_color='grey')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5860acf8",
   "metadata": {},
   "source": [
    "### Degree distribution\n",
    "Let's display the degree distribution of K. We can access the degree of a node `n` simply with `K.degree[n]`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3d678f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "K.degree[3]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb93c4d8",
   "metadata": {},
   "source": [
    "Let's compute the degrees then! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a21c57cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "degrees = [K.degree[n] for n in K.nodes]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "787379ad",
   "metadata": {},
   "source": [
    "... And print the degree distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65133a67",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(degrees, bins=20)\n",
    "plt.xlabel('degree')\n",
    "plt.ylabel('#')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "04748e1d",
   "metadata": {},
   "source": [
    "## *Bonus*: number of walks between two nodes\n",
    "The Dijsktra's algorithm uses a heuristic to find the shortest path between 2 nodes in a connected graph. Another important property to quantify the indirect relation between two nodes is the number of walks between them. The number of walks of length k between nodes $i$ and $j$ is the $(i,j)$-entry of the $k$-th power of the adjacency matrix $A$, i.e  $A^k(i, j)$.\n",
    "\n",
    "Let's verify this property for any graph of our choice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "625418b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = nx.adjacency_matrix(G)\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d8535ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(A.dot(A).dot(A))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "516f36d1",
   "metadata": {},
   "source": [
    "# 2. Random Graph models & Network motifs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a878c45",
   "metadata": {},
   "source": [
    "## Degree-preserving graph randomisation\n",
    "Degree-preserving graph randomisation is essential to assess whether some patterns in a given graph are inherent to the graph structure, or if they are just the result of *chance*. \n",
    "If the patterns persist even after edge randomisation, they are less likely to be *accidental*. Let's randomise the [Zachary's karate club network](https://en.wikipedia.org/wiki/Zachary%27s_karate_club)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "369cb00f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#fill in the details\n",
    "\n",
    "G = nx.karate_club_graph()\n",
    "print(G)\n",
    "\n",
    "# DRAW THE GRAPH USING NX.DRAW function check the inputs in the documentation \n",
    "# https://networkx.org/documentation/stable/reference/drawing.html\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33abf405",
   "metadata": {},
   "outputs": [],
   "source": [
    "degree_sequence_before_randomisation = sorted((d for n, d in G.degree()), reverse=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5dd90525",
   "metadata": {},
   "source": [
    "The idea is to perform a randomisation via operations that preserve the degrees, as described in https://en.wikipedia.org/wiki/Degree-preserving_randomization. This approach is an alternative to generating a graph via the configuration model for instance.\n",
    "To preserve the degree of each node, the trick is to repeatedly take 2 random edges and swap one of their nodes respectively. We also need to make sure that the 2 selected edges do not have a node in common, as this could introduce a self-loop on the node they have in common (you can convince yourself of this with a simple sketch)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80a821e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "#fill in the details\n",
    "\n",
    "import random as rd \n",
    "def randomise_graph(graph, nb_trials=50):\n",
    "    for _ in range(nb_trials):\n",
    "        selected = rd.sample(list(graph.edges()), 2)\n",
    "        edgeA, edgeB = selected[0], selected[1]\n",
    "        intersection = list(set(edgeA) & set(edgeB))\n",
    "        if len(intersection) == 0:\n",
    "            new_edgeA = (edgeA[0], edgeB[1])\n",
    "            \n",
    "            # add a line for the second new edge after the swap, call it new_edgeB\n",
    "            \n",
    "            \n",
    "            if new_edgeA not in graph.edges() and new_edgeB not in graph.edges():\n",
    "                graph.remove_edges_from([edgeA, edgeB])\n",
    "                \n",
    "                #the previous line removed the old eges. Now add a line that adds the new edges\n",
    "                \n",
    "                \n",
    "                \n",
    "    return graph \n",
    "\n",
    "\n",
    "#apply randomise_graph function to graph G, set number of trails to 50\n",
    "\n",
    "\n",
    "degree_sequence_after_randomisation = sorted((d for n, d in G.degree()), reverse=True) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10e492a6",
   "metadata": {},
   "source": [
    "Let's make sure that the degree distribution is well conserved..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5ad6f57",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(10, 4))\n",
    "\n",
    "ax0.hist(degree_sequence_before_randomisation, histtype='stepfilled', color='grey', alpha=.4)\n",
    "ax0.set_title(\"Degree distribution before randomisation\")\n",
    "ax0.set_xlabel(\"degree\")\n",
    "ax0.set_ylabel(\"# nodes\")\n",
    "\n",
    "ax1.hist(degree_sequence_after_randomisation, histtype='stepfilled', color='orange', alpha=.4)\n",
    "ax1.set_title(\"Degree distribution after randomisation\")\n",
    "ax1.set_xlabel(\"degree\")\n",
    "\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab66e63f",
   "metadata": {},
   "source": [
    "### Number of triangles\n",
    "We are going to check if the number of triangles is conserved after randomisation (see [nx.triangles](https://networkx.org/documentation/stable/reference/algorithms/generated/networkx.algorithms.cluster.triangles.html) for documentation). Intuitively, redirecting edges should break triangles. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3bfb0fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_number_triangles(graph):\n",
    "    triangles_per_node = nx.triangles(graph) \n",
    "    # note that each triangle is counted 3 times (for each of its nodes)\n",
    "    number_of_triangles = sum(triangles_per_node.values()) / 3 \n",
    "    return number_of_triangles\n",
    "actual_number_of_triangles = compute_number_triangles(nx.karate_club_graph())\n",
    "print('actual number of triangles:', actual_number_of_triangles)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "898cb1ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "#fill in the details\n",
    "\n",
    "iterations = 100\n",
    "expected_number_of_triangles = 0\n",
    "\n",
    "#use the for loop randomises the graph (using the function from before) and computes the expected number\n",
    "#of triangles over all iterations (using the function above).\n",
    "\n",
    "for _ in range(iterations):\n",
    "    \n",
    "\n",
    "    \n",
    "\n",
    "print('expected number of triangles:', expected_number_of_triangles / iterations) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "537b4030",
   "metadata": {},
   "source": [
    "In expectation over several randomisations, the number of triangles has decreased. This means that the selected social network displays an intrinsic **social mechanism** that is quite common: **the friend of my friend is likelier to be my friend, compared to an individual selected at random in the population**. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "610dbabc",
   "metadata": {},
   "source": [
    "# 3. Community structure\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b35478ab",
   "metadata": {},
   "source": [
    "## Accuracy of the Louvain algorithm\n",
    "On SBMs with two communities, we explore the accuracy of Louvain to find back the two communities as a function of the parameters $p_{in}$ and $p_{out}$ of the affinity matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54ef4181",
   "metadata": {},
   "outputs": [],
   "source": [
    "#pip install python-louvain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75fbaf23",
   "metadata": {},
   "outputs": [],
   "source": [
    "import networkx as nx\n",
    "from community import community_louvain\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.metrics.cluster import normalized_mutual_info_score as nmi_score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3a18b04",
   "metadata": {},
   "outputs": [],
   "source": [
    "#fill in details!\n",
    "\n",
    "def detect_communities(G):\n",
    "    \"\"\" \n",
    "    This function applies the Louvain algorithm to a given nx.Graph G. \n",
    "    \"\"\"\n",
    "    #use the community_louvain.best_partition function to find partition and then store communities in a\n",
    "    #separate list that is returned by the function.\n",
    "\n",
    "    \n",
    "    \n",
    "    return communities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6bdabc59",
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_accuracy_louvain(left_size, right_size, it=3):\n",
    "    \"\"\" \n",
    "    This function applies the Louvain algorithm to many network simulations, \n",
    "    for different values of p_in and p_out.\n",
    "    left and right_size are the sizes of the two communities.\n",
    "    \"\"\"\n",
    "    ground_truth = [0] * left_size + [1] * right_size\n",
    "    community_sizes = [left_size, right_size]\n",
    "    p_in_range = np.arange(.2, 1., .1)\n",
    "    p_out_range = np.arange(.1, .9, .1)\n",
    "    accuracy = {}\n",
    "    for p_in in p_in_range:\n",
    "        accuracy[p_in] = {}\n",
    "        for p_out in p_out_range:\n",
    "            accuracy[p_in][p_out] = 0\n",
    "            affinity_matrix = [[p_in, p_out], \n",
    "                               [p_out, p_in]]\n",
    "            for _ in range(it):\n",
    "                G = nx.stochastic_block_model(community_sizes, affinity_matrix)\n",
    "                communities = detect_communities(G)\n",
    "                accuracy[p_in][p_out] += nmi_score(communities, ground_truth) / it\n",
    "    return accuracy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61b083a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "#fill in the details\n",
    "\n",
    "left_size = 100\n",
    "right_size = 100\n",
    "\n",
    "#calculate accuracy using the function above, set iterations to 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1247a5b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw_results(accuracy):\n",
    "    heatmap = np.zeros((len(accuracy), len(accuracy)))\n",
    "    point_cloud = {'x': [], 'y': []}\n",
    "    for idx, x in enumerate(accuracy):\n",
    "        for idy, y in enumerate(accuracy[x]):\n",
    "            acc = accuracy[x][y]\n",
    "            heatmap[idx][idy] = acc\n",
    "            point_cloud['x'] += [x / y]\n",
    "            point_cloud['y'] += [acc]\n",
    "\n",
    "    ax0 = sns.heatmap(heatmap, linewidth=0.5)\n",
    "    ax0.set_xticklabels(range(20, 91, 10))\n",
    "    ax0.set_yticklabels(range(10, 81, 10))\n",
    "    ax0.set_ylabel('P_in (%)')\n",
    "    ax0.set_xlabel('P_out (%)')\n",
    "    ax0.set_title('NMI as a function of P_in & P_out')\n",
    "    plt.show()\n",
    "\n",
    "    plt.scatter(point_cloud['x'], point_cloud['y'])\n",
    "    plt.xlabel('P_in/P_out')\n",
    "    plt.ylabel('nmi')\n",
    "    plt.title('NMI as a function of P_in / P_out')\n",
    "    plt.show() \n",
    "\n",
    "draw_results(accuracy)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30de6dbf",
   "metadata": {},
   "source": [
    "We observe that **the Louvain algorithm is a good community detection method for assortative networks**, i.e. networks where nodes connect more to nodes in their community, than with nodes in other communities ($p_{in} > p_{out}$)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4754f6e",
   "metadata": {},
   "source": [
    "# 4. Linear dynamics and time-scale separation\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74cd62e7",
   "metadata": {},
   "source": [
    "The master equation for continuous time random walk is $x(t) = e^{t L} x_0$, where $L$ is the normalised laplacian.\n",
    "\n",
    "In this section, we try reproduce Figure B below: we let the master equation evolve, when the walks starts all on the President (top row) or the Instructor (bottom row) at $t=0$. We plot the distribution of location of the walks at time $t$ for 3 different times. We use colors to distinguish between blue and red nodes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62993c2e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Image\n",
    "Image(\"img.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1210318",
   "metadata": {},
   "outputs": [],
   "source": [
    "#fill in the details \n",
    "\n",
    "from scipy.sparse.csgraph import laplacian\n",
    "from scipy.linalg import expm\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "G = nx.karate_club_graph()\n",
    "color_mapping = []\n",
    "for idx in range(34):\n",
    "    side = G.nodes[idx][\"club\"]\n",
    "    if side == \"Mr. Hi\":\n",
    "        color_mapping += [\"blue\"]\n",
    "    elif side == \"Officer\":\n",
    "        color_mapping += [\"red\"] \n",
    "\n",
    "        \n",
    "#get the adjacency matrix of G (you want it dense) and call it A\n",
    "\n",
    "\n",
    "#get the laplacian matrix of A (you want to set the normed parameter to true) and call it L\n",
    "\n",
    "\n",
    "\n",
    "x0_president = np.array([0] * 33 + [1]).T\n",
    "x0_instructor = np.array([1] + [0] * 33).T\n",
    "\n",
    "fig, ax = plt.subplots(ncols=3, nrows=2, figsize=(12, 8))  \n",
    "node_ids = list(range(1, 35))\n",
    "for _, t in enumerate([.38, 5.42, 105.56]):\n",
    "    exptL = expm(-t*L)\n",
    "    xt_president = exptL.dot(x0_president)\n",
    "    xt_instructor = exptL.dot(x0_instructor)\n",
    "    ax[0][_].bar(node_ids, xt_president, color=color_mapping)\n",
    "    ax[1][_].bar(node_ids, xt_instructor, color=color_mapping)\n",
    "    ax[1][_].set_xlabel('Node ID')\n",
    "ax[0][0].set_ylabel('Probability')\n",
    "ax[1][0].set_ylabel('Probability')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}