D measure

"""
schieber.py
-----------

Python implementation of the distance method in 'Quantification of network
structural dissimilarities', by Schieber et al.

"""

import numpy as np
import networkx as nx
from collections import Counter
from scipy import sparse
from scipy.stats import entropy


def jensen_shannon(dists):
    """Jensen-Shannong entropy of a family of distributions.

    dists is a N by M matrix where each row is the ditribution over a set
    of M elements.

    """
    size = dists.shape[0]
    vec = np.log(dists.sum(axis=0)) - np.log(size)
    first_term = (-1/size) * dists.dot(vec).sum()
    second_term = entropy(dists.T).mean()
    return first_term - second_term


def nnd(graph, dists=None):
    """Compute Network Node Dispersion (NND)."""
    if dists is None:
        dists = node_distance(graph)
    diam = dists.shape[1]
    return jensen_shannon(dists) / np.log(diam + 1)


def node_distance(graph):
    """All shortest path distances.

    The nodes must be labeled by integers from 0 to graph.order() - 1.

    """
    size = graph.order()
    if size < 2:
        return 1

    result = np.zeros((size, size)) # does this need to be sparse?
    dists = nx.shortest_path_length(graph)
    dists = np.array([[dists[n1][n2] if dists[n1][n2] < np.inf else size
                       for n2 in dists[n1]]
                      for n1 in dists])

    for idx, row in enumerate(dists):
        counts = Counter(row)
        result[idx] = [counts[l] for l in range(size)]

    diam = (result.sum(axis=0) > 0).sum()
    result = result[:, :diam]

    return result / size


def alpha_centrality(graph, normalize=False):
    """Bonacich centrality."""
    size = graph.order()
    degrees = graph.degree()
    degrees = np.array([degrees[n] for n in graph.nodes()]) / (size - 1)
    alpha = 1 / size
    exogenous = degrees
    mat = sparse.identity(size) - alpha * nx.adjacency_matrix(graph).T
    res = sparse.linalg.inv(mat.asformat('csc')).dot(exogenous)
    return res if not normalize else res / res.sum()


def pad(array, num_cols):
    """Pad with all-zero columns."""
    rows = array.shape[0]
    cols_to_add = num_cols - array.shape[1]
    return np.hstack([array, np.zeros((rows, cols_to_add))])


def schieber(graph1, graph2, w1, w2, w3=None, complement=False):
    """Distance between two graphs. See eqn 2 in the paper."""
    dists1 = node_distance(graph1)
    dists2 = node_distance(graph2)
    if dists1.shape[1] > dists2.shape[1]:
        pad(dists2, dists1.shape[1])
    elif dists2.shape[1] > dists1.shape[1]:
        pad(dists1, dists2.shape[1])

    first_term = np.vstack([dists1.mean(axis=0), dists2.mean(axis=0)])
    first_term = w1 * np.sqrt(jensen_shannon(first_term) / np.log(2))

    second_term = np.sqrt(nnd(graph1, dists1)) - np.sqrt(nnd(graph2, dists2))
    second_term = w2 * np.abs(second_term)

    if w3 is not None:
        alpha1 = alpha_centrality(graph1, normalize=True)
        alpha2 = alpha_centrality(graph2, normalize=True)
        all_alphas = np.vstack([alpha1, alpha2])
        third_term = np.sqrt(jensen_shannon(all_alphas) / np.log(2))
        if complement:
            alpha_comp1 = alpha_centrality(nx.complement(graph1), normalize=True)
            alpha_comp2 = alpha_centrality(nx.complement(graph2), normalize=True)
            all_alphas = np.vstack([alpha_comp1, alpha_comp2])
            third_term += np.sqrt(jensen_shannon(all_alphas) / np.log(2))
        third_term = w3 * third_term / 2
        return first_term + second_term + third_term
    else:
        return first_term + second_term


def main():
    """Compute distance between pre-computed graphs."""
    graph = nx.karate_club_graph()
    print(schieber(graph, graph, 0.5, 0.5))
    print(schieber(graph, graph, 0.45, 0.45, 0.1))
    print(schieber(graph, graph, 0.45, 0.45, 0.1, True))



if __name__ == '__main__':
    main()