koln network params

compute_network_params.py

#/usr/bin/env python

import numpy as np
import pyunicorn
import matplotlib.pyplot as pl
pl.ion()
data = np.load('all_fmri_data.npy')

def compute_xcorr(data, thresh=95, binary=False):
    '''
    Compute the cross correlation for the given data.
    '''
    xcorr = np.corrcoef(data)
    # remove diagonal elements
    xcorr = xcorr - np.eye(len(xcorr))
    # compute the threshold
    thresh = np.percentile(xcorr, thresh)
    if binary:
        print 'computing binary matrix..'
        xcorr = xcorr >= thresh
    else:
        print 'computing weighted matrix'
        xcorr[xcorr <= thresh] = 0

    return xcorr

#channel_names = range(0, xcorr.shape[0])
#channel_names = [str(i) for i in channel_names]

def compute_network_params(adj_matrix, link_weights=None, channel_names=None):
    '''
    Compute the various network parameters for a given
    adjacency matrix.
    '''
    # initializing the network
    net = pyunicorn.Network(adj_matrix)
    surr = net.ErdosRenyi(n_nodes=net.N, n_links=net.n_links)

    if link_weights:
        link_weights = np.abs(adj_matrix)
        net.set_link_attribute(attribute_name="weight", values=link_weights)

    if channel_names:
        net.set_node_attribute(attribute_name="label", values=channel_names)

    # network parameters for surrogate data (random)
    surr_deg = surr.degree()
    surr_closeness = surr.closeness()
    surr_global_clustering = surr.global_clustering()
    surr_transitivity = surr.transitivity()
    surr_avg_path_length = surr.average_path_length()

    # network params for real data
    deg = net.degree()
    closeness = net.closeness()
    global_clustering = net.global_clustering()
    transitivity = net.transitivity()
    avg_path_length = net.average_path_length()

    #path_lengths = net.path_lengths()
    #local_clustering = net.local_clustering()
    #weighted_local_clustering = net.weighted_local_clustering()
    #between = net.betweenness()
    return deg, closeness, global_clustering, transitivity, avg_path_length, \
        surr_deg, surr_closeness, surr_global_clustering, surr_transitivity, surr_avg_path_length

# initialise some arrays, quick fix hard code expected shapes
xcorr_all = np.zeros((10, 90, 90))
degree = np.zeros((10, 90))
closeness = np.zeros((10, 90))
global_clustering = np.zeros((10))
transitivity = np.zeros((10))
avg_path_length = np.zeros((10))

surr_degree = np.zeros((10, 90))
surr_closeness = np.zeros((10, 90))
surr_global_clustering = np.zeros((10))
surr_transitivity = np.zeros((10))
surr_avg_path_length = np.zeros((10))


# compute netweork parameters for all subjects together
for subj in range(10):
    print 'Running for subject %d' % (subj)
    sub = data[subj]
    xcorr_all[subj] = compute_xcorr(sub, binary=True)
    #compute_network_params(xcorr_all[subj])
    degree[subj], closeness[subj], global_clustering[subj], transitivity[subj], avg_path_length[subj], surr_degree[subj], surr_closeness[subj], surr_global_clustering[subj], surr_transitivity[subj], surr_avg_path_length[subj] = \
        compute_network_params(xcorr_all[subj], link_weights=None, channel_names=None)

labels = np.load('labels.npy')

for i in range(10):
    pl.figure('degree')
    pl.plot(degree[i], '-')
    #pl.plot(surr_degree[i], '--')
    pl.xlabel('Vertices')
    pl.ylabel('Degree of cross correlation network')
    pl.title('Degree')

    pl.figure('closeness')
    pl.plot(closeness[i], '-')
    #pl.plot(surr_closeness[i], '--')
    pl.xlabel('Vertices')
    pl.ylabel('Closeness of cross correlation network')
    pl.title('Closeness')

fig, (ax1, ax2, ax3) = pl.subplots(3)
ax1.set_title('Global clustering')
ax1.plot(global_clustering, 'k')
ax1.plot(surr_global_clustering, '--', color='k')
ax1.set_xlabel('subjects')
ax1.set_ylabel('global clustering')

ax2.set_title('Transitivity')
ax2.plot(transitivity, 'b')
ax2.plot(surr_transitivity, '--', color='b')
ax2.set_xlabel('subjects')
ax2.set_ylabel('transitivity')

ax3.set_title('Average path length')
ax3.plot(avg_path_length, 'y')
ax3.plot(surr_avg_path_length, '--', color='y')
ax3.set_xlabel('subjects')
ax3.set_ylabel('Average path length')

pl.tight_layout()
pl.show()