Source code for networkmeasures

"""
Defines network measures for quantum mutual information matrices.

**Authors**

* David L. Vargas (original version)
* Logan Hillberry
* D. Jaschke (adaption for OSMPS module)
"""
import numpy as np
#import scipy.stats as stats

def networkmeasures(resdict):
    """
    Set the network measures to a result dictionary.

    **Arguments**

    resdict : dictionary
        Contains the results of an MPS simulation. The mutual information
        matrix measurement is required. This measurement is stored as key
        ``MIM``. (Could extend function such that all single and two site
        reduced density matrices are sufficient as well.)

    **Details**

    The following flags will be set (or overwritten if existent):

    * ``nwm_clustering``: Clustering coefficient
      (:py:func:`networkmeasures.clustering`)
    * ``nwm_density``: Density of the mutual information matrix
      (:py:func:`networkmeasures.density`)
    * ``nwm_disparity``: Disparity of the mutual information matrix
      (:py:func:`networkmeasures.disparity`)
    * ``nwm_pearson``: Pearson coefficient and 2p taisl of the mutual
      information matrix (:py:func:`networkmeasures.disparity`)
    """
    if('MIM' not in resdict):
        raise KeyError("Key ``MIM`` is required for setting network measures.")

    resdict['nwm_clustering'] = clustering(resdict['MIM'])
    resdict['nwm_density'] = density(resdict['MIM'])
    resdict['nwm_disparity'] = disparity(resdict['MIM'])
    resdict['nwm_pearson'] = pearson(resdict['MIM'])

    return


def clustering(matrix):
    """
    Calculates the clustering coefficient as it is defined in
    equation (7.39) of Mark Newman's book on networks (page 199).

    **Arguments**

    matrix : 2d numpy array
        Contains the mutual information matrix.
    """
    matrix2 = np.linalg.matrix_power(matrix, 2)
    matrix3 = np.linalg.matrix_power(matrix, 3)

    # Zero out diagonal entries. So we do not count edges as
    # connected triples.
    matrix2 -= np.diag(np.diag(matrix2))
    denominator = np.sum(matrix2)

    numerator = np.trace(matrix3)

    # if there are no closed paths of length three the clustering
    # is automatically set to zero.
    if(numerator == 0.):
        # Covers case where everything is zero
        return 0.
    else:
        return numerator / denominator


def density(matrix):
    """
    Calculates density, also termed connectance in some literature.
    Defined on page 134 of Mark Newman's book on networks.

    **Arguments**

    matrix : 2d numpy array
        Contains the mutual information matrix.
    """
    ll = matrix.shape[0]
    lsq = ll * (ll - 1)
    return np.sum(matrix) / lsq


def disparity(matrix):
    """
    Disparity defined on page 199 of doi:10.1016/j.physrep.2005.10.009
    Equation (2.39), Here I take the average of this quantity over the
    entire network

    **Arguments**

    matrix : 2d numpy array
        Contains the mutual information matrix.
    """
    ll = matrix.shape[0]

    numerator = np.sum(matrix**2, axis=0) / ll
    denominator = np.sum(matrix, axis=0)

    # Logical Check
    logos = (denominator > 0)
    numerator = numerator[logos]
    denominator = denominator[logos]
    denominator = denominator**2

    # Check for zero denominator.
    if(np.sum(denominator) == 0.):
        # sum(denominator) == 0 ==> logos.shape[0] == 0
        return np.nan
    else:
        return np.sum(numerator / denominator)


def pearson(matrix):
    """
    Calculates the Pearsons correlation coefficient and the 2-tailed p-value.
    For definitions see scipy.stats.pearsonr. Function returns a tuple
    with matrix containing Pearson correlation coefficients and a second
    matrix with the 2-tailed p-values.

    **Arguments**

    matrix : 2d numpy array
        Contains the mutual information matrix.
    """
    ll = matrix.shape[0]
    pearsonR = np.zeros((ll, ll))
    pvalues = np.zeros((ll, ll))

    # Pearson coefficients should be symmetric
    for ii in range(ll):
        for jj in range(ii + 1, ll):
            r, p = stats.pearsonr(matrix[ii, :], matrix[jj, :])
            pearsonR[ii][jj] = r
            pearsonR[jj][ii] = r
            pvalues[ii][jj] = p
            pvalues[jj][ii] = p

    return (pearsonR, pvalues)