util_robustness.py

# code adapted from
# https://github.com/tsantanam/networkx-robustness/blob/v0.0.7/src/networkx_robustness/networkx_robustness.py
# with some shortcuts

import networkx as nx
import random

def simulate_random_attack(G=None, attack_fraction=0.1, stepSize=10, weight=None):
    """
    Simulate random attack on a network
    :param G: networkx graph
    :param attack_fraction: fraction of nodes to be attacked (default: 0.1)
    :param stepSize: number of nodes to be removed in each step (default: 20)
    :param weight: weight of edges (default: None)
    :return: initial (float), frac (list), apl (list)
    """
    # copy the graph to avoid changing the original graph
    G = G.copy()
    G = G.to_undirected() # make sure the graph is undirected for the robustness analysis
    # get the  number of nodes
    G_nodes = G.number_of_nodes()

    # get the initial fraction of nodes in the largest connected component
    lcc_size = len(max(nx.connected_components(G), key=len))
    initial = lcc_size / G_nodes
    try:
        aveShortPath = nx.average_shortest_path_length(G, weight=weight)
    except nx.NetworkXError:
        aveShortPath = -1

    # initialize lists
    frac = [initial]
    apl = [aveShortPath]
    numNodes = [G_nodes]

    # simulate random attack in batches
    for i in range(0, int(G_nodes * attack_fraction), stepSize):
        G.remove_nodes_from(random.sample(list(G.nodes()), stepSize))
        # get the largest component
        Gc = G.subgraph(max(nx.connected_components(G), key=len))
        # calculate metrics
        numNodes.append(G.number_of_nodes())
        frac.append(Gc.number_of_nodes() / G.number_of_nodes())
        if Gc.number_of_nodes() > 1:
            apl.append(nx.average_shortest_path_length(Gc, weight=weight))
        else:
            apl.append(0)
    return frac, apl, numNodes

def simulate_degree_attack(G=None, attack_fraction=0.1, stepSize=10, weight=None):
    """
    Simulate random attack on a network
    :param G: networkx graph
    :param attack_fraction: fraction of nodes to be attacked (default: 0.1)
    :param stepSize: number of nodes to be removed in each step (default: 20)
    :param weight: weight of edges (default: None)
    :return: initial (float), frac (list), apl (list)
    """
    # copy the graph to avoid changing the original graph
    G = G.copy()
    G = G.to_undirected() # make sure the graph is undirected for the robustness analysis
    # get the  number of nodes
    G_nodes = G.number_of_nodes()

    # get the initial fraction of nodes in the largest connected component
    lcc_size = len(max(nx.connected_components(G), key=len))
    initial = lcc_size / G_nodes
    try:
        aveShortPath = nx.average_shortest_path_length(G, weight=weight)
    except nx.NetworkXError:
        aveShortPath = -1

    # initialize lists
    frac = [initial]
    apl = [aveShortPath]
    numNodes = [G_nodes]

    #get the degree of each node
    degree = nx.degree_centrality(G)
    # sort the nodes by degree
    degree = sorted(degree, key=degree.get, reverse=True)

    # simulate degree attack
    for i in range(0, int(G_nodes * attack_fraction)):
        # remove the node with the highest degree
        G.remove_nodes_from(degree[i:i+stepSize])
        # get the largest connected component of G
        Gc = G.subgraph(max(nx.connected_components(G), key=len))
        # calculate metrics
        numNodes.append(G.number_of_nodes())
        frac.append(Gc.number_of_nodes() / G.number_of_nodes())
        if Gc.number_of_nodes() > 1:
            apl.append(nx.average_shortest_path_length(Gc, weight=weight))
        else:
            apl.append(0)
    return frac, apl, numNodes

def molloy_reed(G=None):
    """
    Compute the Molloy-Reed criterion for a network
    :param G: networkx graph
    :return: Molloy-Reed criterion
    """
    # get the average squared degree
    avg_sq_degree = sum([d ** 2 for n, d in G.degree()]) / G.number_of_nodes()
    # get the average degree
    avg_degree = sum([d for n, d in G.degree()]) / G.number_of_nodes()
    # compute the Molloy-Reed criterion
    molloy_reed = avg_sq_degree/avg_degree

    return molloy_reed

def critical_threshold(G=None):
    """
    Compute the critical threshold for a network
    :param G: networkx graph
    :return: critical threshold
    """
    # get the average squared degree
    avg_sq_degree = sum([d ** 2 for n, d in G.degree()]) / G.number_of_nodes()
    # get the average degree
    avg_degree = sum([d for n, d in G.degree()]) / G.number_of_nodes()
    # compute the Molloy-Reed criterion
    molloy_reed = avg_sq_degree/avg_degree
    # compute the critical threshold
    critical_threshold = 1 - (1/(molloy_reed-1))

    return critical_threshold