"""
ipsen_mikhailov.py
--------------------------
Graph distance based on paper:
Evolutionary reconstruction of network
Available here:
https://journals.aps.org/pre/abstract/10.1103/PhysRevE.66.046109
author: Guillaume St-Onge
email: guillaume.st-onge.4@ulaval.ca
Submitted as part of the 2019 NetSI Collabathon.
"""
import numpy as np
import networkx as nx
from .base import BaseDistance
from scipy.sparse.csgraph import laplacian
from scipy.linalg import eigh
from scipy.integrate import quad
from ..utilities.graph import unweighted
[docs]class IpsenMikhailov(BaseDistance):
"""Compares the spectrum of the Laplacian matrices."""
[docs] @unweighted
def dist(self, G1, G2, hwhm=0.08):
"""Compare the spectrum ot the associated Laplacian matrices.
The results dictionary also stores a 2-tuple of the underlying
adjacency matrices in the key `'adjacency_matrices'`.
Parameters
----------
G1, G2 (nx.Graph)
two networkx graphs to be compared.
hwhm (float)
half with at half maximum of the lorentzian kernel.
Returns
-------
dist (float)
the distance between G1 and G2.
Notes
-----
Requires undirected networks.
References
----------
.. [1] https://journals.aps.org/pre/abstract/10.1103/PhysRevE.66.046109
"""
# get the adjacency matrices
adj1 = nx.to_numpy_array(G1)
adj2 = nx.to_numpy_array(G2)
self.results['adjacency_matrices'] = adj1, adj2
# get the IM distance
dist = _im_distance(adj1, adj2, hwhm)
self.results['dist'] = dist
return dist
def _im_distance(adj1, adj2, hwhm):
"""Computes the Ipsen-Mikhailov distance for two symmetric adjacency
matrices
Base on this paper :
https://journals.aps.org/pre/abstract/10.1103/PhysRevE.66.046109
Note : this is also used by the file hamming_ipsen_mikhailov.py
Parameters
----------
adj1, adj2 (array): adjacency matrices.
hwhm (float) : hwhm of the lorentzian distribution.
Returns
-------
dist (float) : Ipsen-Mikhailov distance.
"""
N = len(adj1)
# get laplacian matrix
L1 = laplacian(adj1, normed=False)
L2 = laplacian(adj2, normed=False)
# get the modes for the positive-semidefinite laplacian
w1 = np.sqrt(np.abs(eigh(L1)[0][1:]))
w2 = np.sqrt(np.abs(eigh(L2)[0][1:]))
# we calculate the norm for both spectrum
norm1 = (N - 1) * np.pi / 2 - np.sum(np.arctan(-w1 / hwhm))
norm2 = (N - 1) * np.pi / 2 - np.sum(np.arctan(-w2 / hwhm))
# define both spectral densities
density1 = lambda w: np.sum(hwhm / ((w - w1) ** 2 + hwhm**2)) / norm1
density2 = lambda w: np.sum(hwhm / ((w - w2) ** 2 + hwhm**2)) / norm2
func = lambda w: (density1(w) - density2(w)) ** 2
return np.sqrt(quad(func, 0, np.inf, limit=100)[0])