Network_Science.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "590460b7",
   "metadata": {},
   "source": [
    "# **Resilience of Delhi Road Networks to Traffic Disruptions**\n",
    "\n",
    "---\n",
    ">### **Problem Statement**:\n",
    "\n",
    "* Road networks are complex systems that are essential for movement of people and goods.\n",
    "* They are also highly interconnected, making them vulnerable to disruptions caused by accidents, natural disasters, or infrastructure failures.\n",
    "* These disruptions can have a significant impact on the economy and society, and can lead to many problems.\n",
    "\n",
    ">### **Objective of Project**:\n",
    "\n",
    "* The objective of this project is to examine the Resilience of road networks to traffic disruptions from a network science perspective.\n",
    "* This will involve analyzing the network’s structure and investigating strategies for improving the network’s robustness.\n",
    "\n",
    ">### **Background**:\n",
    "\n",
    "* Network science is a field of study that examines the structure and function of complex networks, with Road networks being a prime example.\n",
    "*   This knowledge allows us to assess the resilience of Road networks, referring to their ability to maintain connectivity and function at an acceptable level despite disruptions.\n",
    "\n",
    ">### **Outcome**:\n",
    "\n",
    "* The outcome of this project will be a better understanding of the network science factors that contribute to the resilience of road networks.\n",
    "* This knowledge can be used to develop strategies for improving the resilience of existing networks, and to design new Road networks that are more resilient to disruptions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "579183f3",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import networkx as nx\n",
    "from geopy.distance import geodesic\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "routes = pd.read_csv('./GTFS/routes.txt')\n",
    "stop_times = pd.read_csv('./GTFS/stop_times.txt')\n",
    "stops = pd.read_csv('./GTFS/stops.txt')\n",
    "trips = pd.read_csv('./GTFS/trips.txt')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b6a5c41",
   "metadata": {},
   "source": [
    "## Data Preprocessing\n",
    "\n",
    "* Data preprocessing involves cleaning, organizing, and structuring these datasets to facilitate effective analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83810d19",
   "metadata": {},
   "outputs": [],
   "source": [
    "routes = routes[['agency_id','route_id', 'route_long_name']]\n",
    "\n",
    "print(\"Routes Data:\")\n",
    "routes.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6ec74bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "stop_times = stop_times[['trip_id', 'stop_id', 'stop_sequence']]\n",
    "\n",
    "print(\"Stop Times Data:\")\n",
    "stop_times.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26df908b",
   "metadata": {},
   "outputs": [],
   "source": [
    "stops = stops[['stop_code', 'stop_id', 'stop_lat' , 'stop_lon', 'stop_name']]\n",
    "\n",
    "print(\"Stops Data:\")\n",
    "stops.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "adab9cc4",
   "metadata": {},
   "outputs": [],
   "source": [
    "trips = trips[['route_id', 'service_id', 'trip_id']]\n",
    "\n",
    "print(\"Trips Data:\")\n",
    "trips.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7fbbc7b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a dictionary of route_id to route_long_name\n",
    "route_name = dict(zip(routes['route_id'], routes['route_long_name']))\n",
    "\n",
    "# Create a dictionary of trip_id to route_id\n",
    "trip_to_route = dict(zip(trips['trip_id'], trips['route_id']))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db57569f",
   "metadata": {},
   "source": [
    "## Create Directed Graph\n",
    "\n",
    "* We create a directed graph. In this graph the Nodes represent stops and Edges are the roads connecting these stops.\n",
    "* We assign attributes like Stop ID, Stop Name, Stop Latitude, Stop Longitude, Degree to these Nodes. \n",
    "* The Edges are inserted between stops. We assign attributes like Trips, Routes, and Distance between the two Stops to these Edges."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c168449",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a directed graph\n",
    "G = nx.DiGraph()\n",
    "\n",
    "# Insert Nodes\n",
    "for index, row in stops.iterrows():\n",
    "    stop_id = row['stop_id']\n",
    "    G.add_node(stop_id, type='stop', stop_name=row['stop_name'], stop_lat=float(row['stop_lat']), stop_lon=float(row['stop_lon']), in_degree=0, out_degree=0, degree=0)    \n",
    "\n",
    "# Insert Edges\n",
    "prev_trip_id = None\n",
    "prev_stop_id = None\n",
    "prev_stop_seq = 0\n",
    "\n",
    "for _, row in stop_times.iterrows():    \n",
    "    cur_trip_id = row['trip_id']\n",
    "    cur_stop_id = row['stop_id']\n",
    "    cur_stop_seq = int(row['stop_sequence'])\n",
    "    \n",
    "    if prev_stop_id is not None:   \n",
    "        if prev_trip_id == cur_trip_id:\n",
    "            if (prev_stop_seq + 1) == cur_stop_seq:\n",
    "                if prev_stop_id != cur_stop_id:\n",
    "                    route_id = trip_to_route[cur_trip_id]\n",
    "                    if G.has_edge(prev_stop_id, cur_stop_id):\n",
    "                        if route_id in G[prev_stop_id][cur_stop_id]['trip_dict']:\n",
    "                            G[prev_stop_id][cur_stop_id]['trip_dict'][route_id] += (cur_trip_id,)\n",
    "                        else:\n",
    "                            G[prev_stop_id][cur_stop_id]['trip_dict'][route_id] = (cur_trip_id,)\n",
    "                    else:    \n",
    "                        G.add_edge(prev_stop_id, cur_stop_id, trip_dict={route_id: (cur_trip_id,)}, distance=0)\n",
    "                    G.nodes[prev_stop_id]['out_degree'] += 1\n",
    "                    G.nodes[prev_stop_id]['degree'] = G.nodes[prev_stop_id]['in_degree'] + G.nodes[prev_stop_id]['out_degree']\n",
    "                    G.nodes[cur_stop_id]['in_degree'] += 1\n",
    "                    G.nodes[cur_stop_id]['degree'] = G.nodes[cur_stop_id]['in_degree'] + G.nodes[cur_stop_id]['out_degree']  \n",
    "                    \n",
    "        else:\n",
    "            prev_stop_seq = 0\n",
    "    else:\n",
    "        prev_stop_seq = 0\n",
    "     \n",
    "    prev_stop_id = cur_stop_id\n",
    "    prev_trip_id = cur_trip_id\n",
    "    prev_stop_seq = cur_stop_seq\n",
    "    \n",
    "print(\"Number of nodes:\", G.number_of_nodes())\n",
    "print(\"Number of edges:\", G.number_of_edges())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4718a4d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "print( len(G[488][233]['trip_dict']))\n",
    "print(G[488][233]['trip_dict'][142])\n",
    "print( G.nodes[488]['out_degree'])\n",
    "print( G.nodes[488]['in_degree'])\n",
    "print( G.nodes[488]['degree'])\n",
    "print(G.nodes[488])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19a065a1",
   "metadata": {},
   "source": [
    "## Degree distribution\n",
    "\n",
    "* Degree distribution is a statistical measure that describes how the number of connections (degrees) is distributed among nodes in a network.\n",
    "\n",
    "* The degree of a node is the number of edges connected to it. Degree distribution describes the likelihood that a randomly chosen node in the network has a certain degree.\n",
    "\n",
    "* Understanding degree distribution helps reveal important properties of the network, such as its **connectivity, robustness, and overall structure**.\n",
    "\n",
    "* Networks can have different degree distributions. Two common types are Normal Distribution and Power-law Distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86068cb5",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Extract degrees of all nodes\n",
    "degree_sequence = [G.nodes[node]['degree'] for node in G.nodes]\n",
    "\n",
    "# Create a scatter plot\n",
    "plt.scatter(range(len(degree_sequence)), degree_sequence, color='red', alpha=0.7)\n",
    "plt.title('Degree Scatter Plot', fontsize=15, fontweight='bold')\n",
    "plt.xlabel('Node ID', fontsize=12, fontweight='bold')\n",
    "plt.ylabel('Degree', fontsize=12, fontweight='bold')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ba6788c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot histogram\n",
    "plt.hist(degree_sequence, bins=20, color='red')\n",
    "plt.title('Degree Distribution Histogram', fontsize=15, fontweight='bold')\n",
    "plt.xlabel('Degree', fontsize=12, fontweight='bold')\n",
    "plt.ylabel('Number of Nodes', fontsize=12, fontweight='bold')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08aab375",
   "metadata": {},
   "source": [
    "## Power Law Distribution\n",
    "\n",
    "*   Some real-world networks exhibit a power-law degree distribution, meaning that a small number of nodes have very high degrees, while the majority have relatively low degrees.\n",
    "\n",
    "*   The power-law distribution is given by $(p_k \\sim k^{-\\alpha})$ and the ${\\alpha}$ is its degree exponent.\n",
    "\n",
    "*   The power-law exponent ${\\alpha}$ is a parameter that characterizes the shape of the power-law distribution. A smaller α indicates a steeper decline in the distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8df18693",
   "metadata": {},
   "outputs": [],
   "source": [
    "import powerlaw\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fit = powerlaw.Fit(degree_sequence, discrete=True)\n",
    "fit.plot_ccdf(color='r', linewidth=1)\n",
    "\n",
    "plt.title('Power Law Fit for Degree Distribution', fontsize=15, fontweight='bold')\n",
    "plt.xlabel('Degree (log k)', fontsize=12, fontweight='bold')\n",
    "plt.ylabel('Probability Density (log pk)', fontsize=12, fontweight='bold')\n",
    "plt.show()\n",
    "\n",
    "alpha = fit.power_law.alpha\n",
    "print(f\"Power-law exponent (alpha): {alpha}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93982bce",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fit the power-law distribution to the data\n",
    "fit = powerlaw.Fit(degree_sequence, xmin=1)\n",
    "\n",
    "# Plot the power-law fit\n",
    "fit.plot_pdf(color='r', linewidth=1)\n",
    "\n",
    "# Customize the plot\n",
    "plt.title('Power Law Fit for Degree Distribution', fontsize=15, fontweight='bold')\n",
    "plt.xlabel('Degree', fontsize=12, fontweight='bold')\n",
    "plt.ylabel('Probability Density', fontsize=12, fontweight='bold')\n",
    "plt.legend(['Power-law Fit', 'Degree Distribution'])\n",
    "plt.show()\n",
    "\n",
    "alpha = fit.power_law.alpha\n",
    "print(f\"Power-law exponent (alpha): {alpha}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d49ee40",
   "metadata": {},
   "source": [
    "## Hubs\n",
    "\n",
    "* Node degree centrality is a measure of the importance of a node in a network based on its degree, i.e., the number of edges connected to it.\n",
    "* Hubs in a network refer to nodes that have a significantly higher node degree centrality compared to other nodes. \n",
    "* They play a crucial role in maintaining connectivity and influencing the flow of information.\n",
    "* The concept of hubs can be related to the **power-law distribution** by considering nodes with high degrees.\n",
    "* The centrality of a node i in a network can be represented as: $$C_i = \\frac{\\text{{Degree of Node }} i}{\\text{{Maximum Degree in the Network}}}$$\n",
    "* This provides a normalized measure of how central a node is based on its degree."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f5fb02e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate node degree centrality\n",
    "node_degree_centrality = {node: data['degree'] for node, data in G.nodes(data=True)}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f6a7b20",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the top 50 nodes based on degree centrality\n",
    "dc_top_nodes = sorted(node_degree_centrality, key=node_degree_centrality.get, reverse=True)[:50]\n",
    "\n",
    "for node in dc_top_nodes[:5]:\n",
    "    print(G.nodes[node])\n",
    "\n",
    "# Create a subgraph\n",
    "dc_subgraph = G.subgraph([node for node in dc_top_nodes if G.nodes[node]['type'] == 'stop'])\n",
    "dc_pos = nx.kamada_kawai_layout(dc_subgraph)\n",
    "\n",
    "dc_node_labels = {node: data['stop_name'] for node, data in dc_subgraph.nodes(data=True)}\n",
    "dc_node_colors = 'red'\n",
    "\n",
    "# Plot the subgraph\n",
    "fig, ax = plt.subplots(figsize=(30, 30))\n",
    "nx.draw(dc_subgraph, pos=dc_pos, labels=dc_node_labels, with_labels=True, node_size=800, node_color=dc_node_colors, alpha=1, font_size=20, font_color='black', font_weight='bold', edge_color='black', ax=ax)\n",
    "plt.title('Hubs', fontsize=40, fontweight='bold')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fcb00d09",
   "metadata": {},
   "source": [
    "## Betweenness Centrality\n",
    "\n",
    "* Betweenness centrality is a measure of a node's centrality in a network based on its role in facilitating the flow of information or influence between other nodes.\n",
    "* Nodes with high betweenness centrality have a significant influence on the communication or interaction between other nodes in the network.\n",
    "* Betweenness centrality is calculated by determining the number of shortest paths that pass through a node, relative to the total number of shortest paths in the network.\n",
    "* Nodes with high betweenness centrality often act as bridges or bottlenecks in the network, as they control the flow of information between different parts of the network."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df393036",
   "metadata": {},
   "outputs": [],
   "source": [
    "from geopy.distance import geodesic\n",
    "\n",
    "# Function to calculate distance between two points using Haversine formula\n",
    "def haversine_distance(lat1, lon1, lat2, lon2):\n",
    "    coords_1 = (lat1, lon1)\n",
    "    coords_2 = (lat2, lon2)\n",
    "    return geodesic(coords_1, coords_2).km\n",
    "\n",
    "# Iterate through edges and calculate distance\n",
    "for edge in G.edges(data=True):\n",
    "    source_node, target_node, data = edge\n",
    "    source_coords = (G.nodes[source_node]['stop_lat'], G.nodes[source_node]['stop_lon'])\n",
    "    target_coords = (G.nodes[target_node]['stop_lat'], G.nodes[target_node]['stop_lon'])\n",
    "    \n",
    "    # Calculate distance using Haversine formula\n",
    "    distance = haversine_distance(*source_coords, *target_coords)\n",
    "    \n",
    "    # Assign distance as edge weight\n",
    "    G[source_node][target_node]['distance'] = distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73381b41",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate node betweenness centrality based on the 'distance' attribute\n",
    "node_betweenness_centrality = nx.betweenness_centrality(G, weight='distance')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9a0b062",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the top 50 nodes based on betweenness centrality\n",
    "bc_top_nodes = sorted(node_betweenness_centrality, key=node_betweenness_centrality.get, reverse=True)[:50]\n",
    "\n",
    "for node in bc_top_nodes[:5]:\n",
    "    print(G.nodes[node])\n",
    "\n",
    "# Create a subgraph\n",
    "bc_subgraph = G.subgraph([node for node in bc_top_nodes if G.nodes[node]['type'] == 'stop'])\n",
    "bc_pos = nx.kamada_kawai_layout(bc_subgraph)\n",
    "\n",
    "bc_node_labels = {node: data['stop_name'] for node, data in bc_subgraph.nodes(data=True)}\n",
    "bc_node_colors = 'red'\n",
    "\n",
    "# Plot the subgraph\n",
    "fig, ax = plt.subplots(figsize=(30, 30))\n",
    "nx.draw(bc_subgraph, pos=bc_pos, labels=bc_node_labels, with_labels=True, node_size=800, node_color=bc_node_colors, alpha=1, font_size=20, font_color='black', font_weight='bold', edge_color='black', ax=ax)\n",
    "plt.title('Betweenness Centrality', fontsize=40, fontweight='bold')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3576769b",
   "metadata": {},
   "source": [
    "## Disruption Scenarios\n",
    "\n",
    "* Simulate disruptions like road closures and measure their impact on network connectivity.\n",
    "* Visualize the network before and after the disruption to provide a clear representation of the changes."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63904d53",
   "metadata": {},
   "source": [
    ">### **Road Closure**\n",
    "\n",
    "* Road closure, resulting from construction, maintenance, new projects and infrastructure improvements, is a prolonged disruption that can last for a few days, months, or even years.\n",
    "* In order to simulate this kind of disruption we will remove some Edges that are associated with the Nodes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db54f678",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a copy of the graph for road closure scenario\n",
    "G_rc = G.copy()\n",
    "\n",
    "# Identify the top 10 hubs\n",
    "top_hubs = sorted(node_degree_centrality, key=node_degree_centrality.get, reverse=True)[:10]\n",
    "\n",
    "# Function to remove at most 3 edges associated with the top hubs\n",
    "def road_closure(graph, node):\n",
    "    remove_edges = []\n",
    "    no_of_edges = 1\n",
    "    \n",
    "    for successor in graph.successors(node):\n",
    "        if no_of_edges <= 3:\n",
    "            remove_edges.append((node, successor))\n",
    "            no_of_edges += 1\n",
    "        else:\n",
    "            break\n",
    "            \n",
    "    for edge in remove_edges:\n",
    "        src = edge[0]\n",
    "        destn = edge[1]\n",
    "        edge_data = graph[src][destn]\n",
    "        trip_dict = edge_data.get('trip_dict', {})\n",
    "        no_of_trips = 0\n",
    "        \n",
    "        for route_id in trip_dict:\n",
    "            no_of_trips += len(trip_dict[route_id])\n",
    "\n",
    "        graph.nodes[src]['out_degree'] -= no_of_trips\n",
    "        graph.nodes[src]['degree'] -= no_of_trips\n",
    "        graph.nodes[destn]['in_degree'] -= no_of_trips\n",
    "        graph.nodes[destn]['degree'] -= no_of_trips\n",
    "\n",
    "        # Remove the edge\n",
    "        graph.remove_edge(*edge)\n",
    "\n",
    "# Remove edges for the top hub nodes\n",
    "for hub in top_hubs:\n",
    "    road_closure(G_rc, hub)\n",
    "\n",
    "# Print updated graph information\n",
    "print(\"Number of nodes:\", G_rc.number_of_nodes())\n",
    "print(\"Number of edges:\", G_rc.number_of_edges())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2326ba8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Road Closure Analysis\n",
    "fig, axes = plt.subplots(2, 2, figsize=(20, 20))\n",
    "road_closure_degree_sequence = [G_rc.nodes[node]['degree'] for node in G_rc.nodes]\n",
    "\n",
    "# Plot the original degree distribution scatter plot\n",
    "axes[0, 0].scatter(range(len(degree_sequence)), degree_sequence, color='red', alpha=0.7)\n",
    "axes[0, 0].set_title('Scatter Plot Before Road Closure Disruption', fontsize=15, fontweight='bold')\n",
    "axes[0, 0].set_xlabel('Node ID', fontsize=12, fontweight='bold')\n",
    "axes[0, 0].set_ylabel('Degree', fontsize=12, fontweight='bold')\n",
    "\n",
    "# Plot the disrupted degree distribution scatter plot\n",
    "axes[0, 1].scatter(range(len(road_closure_degree_sequence)), road_closure_degree_sequence, color='blue', alpha=0.7)\n",
    "axes[0, 1].set_title('Scatter Plot After Road Closure Disruption', fontsize=15, fontweight='bold')\n",
    "axes[0, 1].set_xlabel('Node ID', fontsize=12, fontweight='bold')\n",
    "axes[0, 1].set_ylabel('Degree', fontsize=12, fontweight='bold')\n",
    "\n",
    "# Plot the original histogram\n",
    "axes[1, 0].hist(degree_sequence, bins=20, color='red', edgecolor='black')\n",
    "axes[1, 0].set_title('Degree Distribution Histogram Before Road Closure Disruption', fontsize=15, fontweight='bold')\n",
    "axes[1, 0].set_xlabel('Degree', fontsize=12, fontweight='bold')\n",
    "axes[1, 0].set_ylabel('Number of Nodes', fontsize=12, fontweight='bold')\n",
    "\n",
    "# Plot the original histogram\n",
    "axes[1, 1].hist(road_closure_degree_sequence, bins=20, color='blue', edgecolor='black')\n",
    "axes[1, 1].set_title('Degree Distribution Histogram After Road Closure Disruption', fontsize=15, fontweight='bold')\n",
    "axes[1, 1].set_xlabel('Degree', fontsize=12, fontweight='bold')\n",
    "axes[1, 1].set_ylabel('Number of Nodes', fontsize=12, fontweight='bold')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6694ec7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(20, 10))\n",
    "\n",
    "# Plot the original hubs graph\n",
    "dc_top_nodes = sorted(node_degree_centrality, key=node_degree_centrality.get, reverse=True)[:50]\n",
    "dc_subgraph = G.subgraph([node for node in dc_top_nodes if G.nodes[node]['type'] == 'stop'])\n",
    "dc_pos = nx.kamada_kawai_layout(dc_subgraph)\n",
    "dc_node_labels = {node: data['stop_name'] for node, data in dc_subgraph.nodes(data=True)}\n",
    "dc_node_colors = 'red'\n",
    "dc_subgraph = G.subgraph([node for node in dc_top_nodes if G.nodes[node]['type'] == 'stop'])\n",
    "dc_pos = nx.kamada_kawai_layout(dc_subgraph)\n",
    "nx.draw(dc_subgraph, pos=dc_pos, labels=dc_node_labels, with_labels=True, node_size=400, node_color=dc_node_colors, alpha=1, font_size=12, font_color='black', font_weight='bold', edge_color='black', ax=axes[0])\n",
    "axes[0].set_title('Hubs Before Road Closure Disruption', fontsize=15, fontweight='bold')\n",
    "\n",
    "# Plot the disrupted hubs graph\n",
    "road_closure_node_degree_centrality = {node: data['degree'] for node, data in G_rc.nodes(data=True)}\n",
    "road_closure_dc_top_nodes = sorted(road_closure_node_degree_centrality, key=road_closure_node_degree_centrality.get, reverse=True)[:50]\n",
    "road_closure_dc_subgraph = G_rc.subgraph([node for node in dc_top_nodes if G_rc.nodes[node]['type'] == 'stop'])\n",
    "road_closure_dc_pos = nx.kamada_kawai_layout(road_closure_dc_subgraph)\n",
    "road_closure_dc_node_labels = {node: data['stop_name'] for node, data in road_closure_dc_subgraph.nodes(data=True)}\n",
    "road_closure_dc_node_colors = 'blue'\n",
    "nx.draw(road_closure_dc_subgraph, pos=road_closure_dc_pos, labels=road_closure_dc_node_labels, with_labels=True, node_size=400, node_color=road_closure_dc_node_colors, alpha=1, font_size=12, font_color='black', font_weight='bold', edge_color='black', ax=axes[1])\n",
    "axes[1].set_title('Hubs After Road Closure Disruption', fontsize=15, fontweight='bold')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a50923f0",
   "metadata": {},
   "source": [
    ">### **Malignant Accidents**\n",
    "\n",
    "* A Road accident may temporarily close a particular route. This results in temporary cancellation of some trips to stops along the route.\n",
    "* To simulate this kind of disruption, we remove some trips associated with the stops connecting the route."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1d2e13f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a copy of the graph for Malignant accident scenario\n",
    "G_ma = G.copy()\n",
    "\n",
    "# Identify the top 20 hubs\n",
    "top_hubs = sorted(node_degree_centrality, key=node_degree_centrality.get, reverse=True)[:20]\n",
    "\n",
    "# Function to remove at most 5 edges associated with the top hubs\n",
    "def road_accident(graph, node):\n",
    "    remove_edges = []\n",
    "    no_of_edges = 1\n",
    "    \n",
    "    for successor in graph.successors(node):\n",
    "        if no_of_edges <= 5:\n",
    "            remove_edges.append((node, successor))\n",
    "            no_of_edges += 1\n",
    "        else:\n",
    "            break\n",
    "    \n",
    "    for edge in remove_edges:\n",
    "        src = edge[0]\n",
    "        destn = edge[1]\n",
    "        edge_data = graph[src][destn]\n",
    "        trip_dict = edge_data.get('trip_dict', {})\n",
    "        no_of_trips = 0\n",
    "        no_of_routes = 0\n",
    "        \n",
    "        # Create a list of trips to remove\n",
    "        trips_to_remove = []\n",
    "        for route_id in trip_dict:\n",
    "            if no_of_routes < 15:\n",
    "                no_of_trips += len(trip_dict[route_id])\n",
    "                trips_to_remove.append(route_id)\n",
    "                no_of_routes += 1\n",
    "            else:\n",
    "                break\n",
    "        \n",
    "        # Remove keys from the dictionary\n",
    "        for key in trips_to_remove:\n",
    "            trip_dict.pop(key, None)\n",
    "\n",
    "        graph.nodes[src]['out_degree'] -= no_of_trips\n",
    "        graph.nodes[src]['degree'] -= no_of_trips\n",
    "        graph.nodes[destn]['in_degree'] -= no_of_trips\n",
    "        graph.nodes[destn]['degree'] -= no_of_trips\n",
    "        graph[src][destn]['trip_dict'] = trip_dict\n",
    "\n",
    "# Remove edges for the top hub nodes\n",
    "for hub in top_hubs:\n",
    "    road_accident(G_ma, hub)\n",
    "\n",
    "# Print updated graph information\n",
    "print(\"Number of nodes:\", G_ma.number_of_nodes())\n",
    "print(\"Number of edges:\", G_ma.number_of_edges())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f915fe6",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(G.nodes[149])\n",
    "print(G_ma.nodes[149])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f148a13",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Road Accidents Analysis\n",
    "fig, axes = plt.subplots(2, 2, figsize=(20, 20))\n",
    "road_accident_degree_sequence = [G_ma.nodes[node]['degree'] for node in G_ma.nodes]\n",
    "\n",
    "# Plot the original degree distribution scatter plot\n",
    "axes[0, 0].scatter(range(len(degree_sequence)), degree_sequence, color='red', alpha=0.7)\n",
    "axes[0, 0].set_title('Scatter Plot Before Road Accident Disruption', fontsize=15, fontweight='bold')\n",
    "axes[0, 0].set_xlabel('Node ID', fontsize=12, fontweight='bold')\n",
    "axes[0, 0].set_ylabel('Degree', fontsize=12, fontweight='bold')\n",
    "\n",
    "# Plot the disrupted degree distribution scatter plot\n",
    "axes[0, 1].scatter(range(len(road_accident_degree_sequence)), road_accident_degree_sequence, color='blue', alpha=0.7)\n",
    "axes[0, 1].set_title('Scatter Plot After Road Accident Disruption', fontsize=15, fontweight='bold')\n",
    "axes[0, 1].set_xlabel('Node ID', fontsize=12, fontweight='bold')\n",
    "axes[0, 1].set_ylabel('Degree', fontsize=12, fontweight='bold')\n",
    "\n",
    "# Plot the original histogram\n",
    "axes[1, 0].hist(degree_sequence, bins=20, color='red', edgecolor='black')\n",
    "axes[1, 0].set_title('Degree Distribution Histogram Before Road Accident Disruption', fontsize=15, fontweight='bold')\n",
    "axes[1, 0].set_xlabel('Degree', fontsize=12, fontweight='bold')\n",
    "axes[1, 0].set_ylabel('Number of Nodes', fontsize=12, fontweight='bold')\n",
    "\n",
    "# Plot the disrupted histogram\n",
    "axes[1, 1].hist(road_accident_degree_sequence, bins=20, color='blue', edgecolor='black')\n",
    "axes[1, 1].set_title('Degree Distribution Histogram After Road Accident Disruption', fontsize=15, fontweight='bold')\n",
    "axes[1, 1].set_xlabel('Degree', fontsize=12, fontweight='bold')\n",
    "axes[1, 1].set_ylabel('Number of Nodes', fontsize=12, fontweight='bold')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "143f1b97",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(20, 10))\n",
    "\n",
    "# Plot the original hubs graph\n",
    "dc_top_nodes = sorted(node_degree_centrality, key=node_degree_centrality.get, reverse=True)[:50]\n",
    "dc_subgraph = G.subgraph([node for node in dc_top_nodes if G.nodes[node]['type'] == 'stop'])\n",
    "dc_pos = nx.kamada_kawai_layout(dc_subgraph)\n",
    "dc_node_labels = {node: data['stop_name'] for node, data in dc_subgraph.nodes(data=True)}\n",
    "dc_node_colors = 'red'\n",
    "dc_subgraph = G.subgraph([node for node in dc_top_nodes if G.nodes[node]['type'] == 'stop'])\n",
    "dc_pos = nx.kamada_kawai_layout(dc_subgraph)\n",
    "nx.draw(dc_subgraph, pos=dc_pos, labels=dc_node_labels, with_labels=True, node_size=400, node_color=dc_node_colors, alpha=1, font_size=12, font_color='black', font_weight='bold', edge_color='black', ax=axes[0])\n",
    "axes[0].set_title('Hubs Before Road Accident Disruption', fontsize=15, fontweight='bold')\n",
    "\n",
    "# Plot the disrupted hubs graph\n",
    "road_accident_node_degree_centrality = {node: data['degree'] for node, data in G_ma.nodes(data=True)}\n",
    "road_accident_dc_top_nodes = sorted(road_accident_node_degree_centrality, key=road_accident_node_degree_centrality.get, reverse=True)[:50]\n",
    "\n",
    "for node in road_accident_dc_top_nodes[:5]:\n",
    "    print(G.nodes[node])\n",
    "\n",
    "road_accident_dc_subgraph = G_ma.subgraph([node for node in road_accident_dc_top_nodes if G_ma.nodes[node]['type'] == 'stop'])\n",
    "road_accident_dc_pos = nx.kamada_kawai_layout(road_accident_dc_subgraph)\n",
    "road_accident_dc_node_labels = {node: data['stop_name'] for node, data in road_accident_dc_subgraph.nodes(data=True)}\n",
    "road_accident_dc_node_colors = 'blue'\n",
    "nx.draw(road_accident_dc_subgraph, pos=road_accident_dc_pos, labels=road_accident_dc_node_labels, with_labels=True, node_size=400, node_color=road_accident_dc_node_colors, alpha=1, font_size=12, font_color='black', font_weight='bold', edge_color='black', ax=axes[1])\n",
    "axes[1].set_title('Hubs After Road Accident Disruption', fontsize=15, fontweight='bold')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3daf4c06",
   "metadata": {},
   "source": [
    ">### **Natural Disasters**\n",
    "\n",
    "* Natural Disaster can vary significantly in duration and impact on road networks.\n",
    "* Short term disruptions like heavy rains have temporary impact like delay in trips. \n",
    "* But Long term disruptions like floods may last for a few days and may result in road closures, make some stops inaccessible and cancellation of all the trips along that route.\n",
    "* To simulate this kind of disruption, we remove some Nodes and make those places inaccessible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3589391d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a copy of the graph for road closure scenario\n",
    "G_nd = G.copy()\n",
    "\n",
    "# Identify the top 5 hubs\n",
    "top_hubs = sorted(node_degree_centrality, key=node_degree_centrality.get, reverse=True)[:5]\n",
    "\n",
    "# Function to remove the hubs\n",
    "def natural_disaster(graph, node):\n",
    "    remove_edges = []\n",
    "    \n",
    "    for successor in graph.successors(node):\n",
    "        remove_edges.append((node, successor))\n",
    "            \n",
    "    for edge in remove_edges:\n",
    "        src = edge[0]\n",
    "        destn = edge[1]\n",
    "        edge_data = graph[src][destn]\n",
    "        trip_dict = edge_data.get('trip_dict', {})\n",
    "        no_of_trips = 0\n",
    "        \n",
    "        for route_id in trip_dict:\n",
    "            no_of_trips += len(trip_dict[route_id])\n",
    "\n",
    "        graph.nodes[src]['out_degree'] -= no_of_trips\n",
    "        graph.nodes[src]['degree'] -= no_of_trips\n",
    "        graph.nodes[destn]['in_degree'] -= no_of_trips\n",
    "        graph.nodes[destn]['degree'] -= no_of_trips\n",
    "\n",
    "        # Remove the edge\n",
    "        graph.remove_edge(*edge)\n",
    "    # Remove the hub\n",
    "    graph.remove_node(node)\n",
    "\n",
    "# Remove the top hub nodes\n",
    "for hub in top_hubs:\n",
    "    natural_disaster(G_nd, hub)\n",
    "\n",
    "# Print updated graph information\n",
    "print(\"Number of nodes:\", G_nd.number_of_nodes())\n",
    "print(\"Number of edges:\", G_nd.number_of_edges())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f2d1ef5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Natural Disaster Analysis\n",
    "fig, axes = plt.subplots(2, 2, figsize=(20, 20))\n",
    "natural_disaster_degree_sequence = [G_nd.nodes[node]['degree'] for node in G_nd.nodes]\n",
    "\n",
    "# Plot the original degree distribution scatter plot\n",
    "axes[0, 0].scatter(range(len(degree_sequence)), degree_sequence, color='red', alpha=0.7)\n",
    "axes[0, 0].set_title('Scatter Plot Before Natural Disaster Disruption', fontsize=15, fontweight='bold')\n",
    "axes[0, 0].set_xlabel('Node ID', fontsize=12, fontweight='bold')\n",
    "axes[0, 0].set_ylabel('Degree', fontsize=12, fontweight='bold')\n",
    "\n",
    "# Plot the disrupted degree distribution scatter plot\n",
    "axes[0, 1].scatter(range(len(natural_disaster_degree_sequence)), natural_disaster_degree_sequence, color='blue', alpha=0.7)\n",
    "axes[0, 1].set_title('Scatter Plot After Natural Disaster Disruption', fontsize=15, fontweight='bold')\n",
    "axes[0, 1].set_xlabel('Node ID', fontsize=12, fontweight='bold')\n",
    "axes[0, 1].set_ylabel('Degree', fontsize=12, fontweight='bold')\n",
    "\n",
    "# Plot the original histogram\n",
    "axes[1, 0].hist(degree_sequence, bins=20, color='red', edgecolor='black')\n",
    "axes[1, 0].set_title('Degree Distribution Histogram Before Natural Disaster Disruption', fontsize=15, fontweight='bold')\n",
    "axes[1, 0].set_xlabel('Degree', fontsize=12, fontweight='bold')\n",
    "axes[1, 0].set_ylabel('Number of Nodes', fontsize=12, fontweight='bold')\n",
    "\n",
    "# Plot the disrupted histogram\n",
    "axes[1, 1].hist(natural_disaster_degree_sequence, bins=20, color='blue', edgecolor='black')\n",
    "axes[1, 1].set_title('Degree Distribution Histogram After Natural Disaster Disruption', fontsize=15, fontweight='bold')\n",
    "axes[1, 1].set_xlabel('Degree', fontsize=12, fontweight='bold')\n",
    "axes[1, 1].set_ylabel('Number of Nodes', fontsize=12, fontweight='bold')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7fc16941",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(20, 10))\n",
    "\n",
    "# Plot the original hubs graph\n",
    "dc_top_nodes = sorted(node_degree_centrality, key=node_degree_centrality.get, reverse=True)[:50]\n",
    "dc_subgraph = G.subgraph([node for node in dc_top_nodes if G.nodes[node]['type'] == 'stop'])\n",
    "dc_pos = nx.kamada_kawai_layout(dc_subgraph)\n",
    "dc_node_labels = {node: data['stop_name'] for node, data in dc_subgraph.nodes(data=True)}\n",
    "dc_node_colors = 'red'\n",
    "nx.draw(dc_subgraph, pos=dc_pos, labels=dc_node_labels, with_labels=True, node_size=400, node_color=dc_node_colors, alpha=1, font_size=12, font_color='black', font_weight='bold', edge_color='black', ax=axes[0])\n",
    "axes[0].set_title('Hubs Before Natural Disaster Disruption', fontsize=15, fontweight='bold')\n",
    "\n",
    "# Plot the disrupted hubs graph\n",
    "natural_disaster_node_degree_centrality = {node: data['degree'] for node, data in G_nd.nodes(data=True)}\n",
    "natural_disaster_dc_top_nodes = sorted(natural_disaster_node_degree_centrality, key=natural_disaster_node_degree_centrality.get, reverse=True)[:50]\n",
    "natural_disaster_dc_subgraph = G_nd.subgraph([node for node in dc_top_nodes if node in G_nd.nodes and G_nd.nodes[node]['type'] == 'stop'])\n",
    "natural_disaster_dc_pos = nx.kamada_kawai_layout(natural_disaster_dc_subgraph)\n",
    "natural_disaster_dc_node_labels = {node: data['stop_name'] for node, data in natural_disaster_dc_subgraph.nodes(data=True)}\n",
    "natural_disaster_dc_node_colors = 'blue'\n",
    "nx.draw(natural_disaster_dc_subgraph, pos=natural_disaster_dc_pos, labels=natural_disaster_dc_node_labels, with_labels=True, node_size=400, node_color=natural_disaster_dc_node_colors, alpha=1, font_size=12, font_color='black', font_weight='bold', edge_color='black', ax=axes[1])\n",
    "axes[1].set_title('Hubs After Natural Disaster Disruption', fontsize=15, fontweight='bold')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c331fc1",
   "metadata": {},
   "source": [
    "## Robustness Metrics\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1dda0087",
   "metadata": {},
   "outputs": [],
   "source": [
    "import networkx as nx\n",
    "import numpy as np\n",
    "\n",
    "# Function to calculate and print degree-based robustness metrics\n",
    "def calculate_degree_robustness_metrics(graph, name):\n",
    "    print(f\"\\nDegree-Based Robustness Metrics for {name}:\\n\")\n",
    "    \n",
    "    # Degree sequence of the graph\n",
    "    degree_sequence = [graph.degree(node) for node in graph.nodes]\n",
    "    \n",
    "    # Average degree\n",
    "    avg_degree = np.mean(degree_sequence)\n",
    "    print(f\"Average Degree: {avg_degree:.5f}\")\n",
    "\n",
    "    # Standard deviation of degree\n",
    "    degree_std_dev = np.std(degree_sequence)\n",
    "    print(f\"Standard Deviation of Degree: {degree_std_dev:.5f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f2b2615",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Calculate and print degree-based robustness metrics for the original graph (G)\n",
    "calculate_degree_robustness_metrics(G, \"Original Graph (G)\")\n",
    "\n",
    "# Calculate and print degree-based robustness metrics for Road closures (G_rc)\n",
    "calculate_degree_robustness_metrics(G_rc, \"Modified Graph after Road Closures (G_rc)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f00cf7d",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Calculate and print degree-based robustness metrics for the original graph (G)\n",
    "calculate_degree_robustness_metrics(G, \"Original Graph (G)\")\n",
    "\n",
    "# Calculate and print degree-based robustness metrics for Malignant Accidents (G_ma)\n",
    "calculate_degree_robustness_metrics(G_ma, \"Modified Graph after Malignant Accidents (G_ma)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "466aefea",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate and print degree-based robustness metrics for the original graph (G)\n",
    "calculate_degree_robustness_metrics(G, \"Original Graph (G)\")\n",
    "\n",
    "# Calculate and print degree-based robustness metrics for Natural Disasters (G_nd)\n",
    "calculate_degree_robustness_metrics(G_nd, \"Modified Graph after Natural Disasters (G_nd)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rerouting Algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "71aa522d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Function to find the shortest path between two nodes in a graph\n",
    "def find_shortest_path(graph, source, target):\n",
    "    try:\n",
    "        shortest_path = nx.shortest_path(graph, source=source, target=target, weight='distance')\n",
    "        return shortest_path\n",
    "    except nx.NetworkXNoPath:\n",
    "        return None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "938616ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find the top 20 nodes in terms of degree in the original graph (G)\n",
    "top_nodes_original = sorted(G.nodes, key=lambda x: G.degree(x), reverse=True)[:20]\n",
    "\n",
    "# Find the top 20 nodes in terms of degree in the disrupted graph (G_rc)\n",
    "top_nodes_disrupted = sorted(G_rc.nodes, key=lambda x: G_rc.degree(x), reverse=True)[:20]\n",
    "\n",
    "# Define positions for the nodes using Kamada-Kawai layout\n",
    "pos_original = nx.kamada_kawai_layout(G)\n",
    "pos_disrupted = nx.kamada_kawai_layout(G_rc)\n",
    "\n",
    "# Find shortest path between two nodes in the original graph (G)\n",
    "original_shortest_path = find_shortest_path(G, source=top_nodes_original[0], target=top_nodes_original[1])\n",
    "\n",
    "# Find shortest path between two nodes in the disrupted graph (G_rc)\n",
    "disrupted_shortest_path = find_shortest_path(G_rc, source=top_nodes_disrupted[0], target=top_nodes_disrupted[1])\n",
    "\n",
    "# Plot the original graph with only the top 20 nodes\n",
    "plt.figure(figsize=(10, 10))\n",
    "nx.draw(G, pos_original, nodelist=top_nodes_original, with_labels=True, font_weight='bold', node_size=500, node_color='red', font_size=8)\n",
    "nx.draw_networkx_edges(G, pos_original, edgelist=[(u, v) for u, v in original_shortest_path], edge_color='black', width=2)\n",
    "plt.title('Original Graph with Shortest Path (Top 20 Nodes)', fontsize=15)\n",
    "plt.show()\n",
    "\n",
    "# Plot the disrupted graph with only the top 20 nodes\n",
    "plt.figure(figsize=(10, 10))\n",
    "nx.draw(G_rc, pos_disrupted, nodelist=top_nodes_disrupted, with_labels=True, font_weight='bold', node_size=500, node_color='blue', font_size=8)\n",
    "nx.draw_networkx_edges(G_rc, pos_disrupted, edgelist=[(u, v) for u, v in disrupted_shortest_path], edge_color='black', width=2)\n",
    "plt.title('Disrupted Graph with Shortest Path (Top 20 Nodes)', fontsize=15)\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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}