In [ ]:
%matplotlib inline
Test script for trying plotting routines involving degree + Connectivity circle plots.
In [ ]:
import numpy as np
import os.path as op
import mne
from jumeg import get_jumeg_path
from jumeg.connectivity import plot_degree_circle, plot_lines_and_blobs
import matplotlib.pyplot as plt
orig_labels_fname = get_jumeg_path() + '/data/desikan_label_names.yaml'
yaml_fname = get_jumeg_path() + '/data/desikan_aparc_cortex_based_grouping.yaml'
con_fname = get_jumeg_path() + '/data/sample,aparc-con.npy'
# real connectivity
con = np.load(con_fname)
con = con[0, :, :, 2] + con[0, :, :, 2].T
degrees = mne.connectivity.degree(con, threshold=0.2)
# test known connections
# con = np.zeros((68, 68))
# con[55, 47] = 0.9 # rostralmiddlefrontal-rh - posteriorcingulate-rh
# con[46, 22] = 0.6 # lateraloccipital-lh - posteriorcingulate-lh
# con = con + con.T
# degrees = mne.connectivity.degree(con, threshold=0.2)
fig, ax = plot_lines_and_blobs(con, degrees, yaml_fname, orig_labels_fname,
figsize=(8, 8), node_labels=True)
In [ ]:
# n_nodes = len(degrees)
# node_order_size = 68
# fig = None
# subplot = 111
# radsize = 5.
# degsize = 10.
# alpha = 0.5
# n_lines = 50
# vmin, vmax = None, None
# node_width = None
# arrow = False
# colormap = 'Blues'
# linewidth = 1.5
# def plot_lines_and_blobs(con, degrees, yaml_fname, orig_labels_fname,
# node_order_size=68, fig=None, subplot=111,
# color='b', cmap='Blues', tight_layout=False,
# alpha=0.5, fontsize_groups=6, textcolor_groups='k',
# radsize=1., degsize=1, show_group_labels=True,
# linewidth=1.5, n_lines=50, node_width=None,
# arrow=False, out_fname='lines_and_blobs.png',
# vmin=None, vmax=None, figsize=None, node_labels=False,
# show=True):
# '''
# Plot connectivity circle plot with a centrality index per node shown as
# blobs along the circulference of the circle, hence the lines and the blobs.
# '''
#
# import yaml
# with open(orig_labels_fname, 'r') as f:
# orig_labels = yaml.safe_load(f)['label_names']
#
# n_nodes = len(degrees)
#
# assert n_nodes == len(orig_labels), 'Mismatch in node names and number.'
# assert n_nodes == len(con), 'Mismatch in n_nodes for con and degrees.'
#
# # read the yaml file with grouping of the various nodes
# if op.isfile(yaml_fname):
# with open(yaml_fname, 'r') as f:
# labels = yaml.safe_load(f)
# else:
# print('%s - File not found.' % yaml_fname)
# sys.exit()
#
# # make list of label_names (without individual cortex locations)
# label_names = [list(lab.values())[0] for lab in labels]
# label_names = [la for l in label_names for la in l]
#
# lh_labels = [name + '-lh' for name in label_names]
# rh_labels = [name + '-rh' for name in label_names]
#
# # save the plot order
# node_order = list()
# node_order.extend(lh_labels[::-1]) # reverse the order
# node_order.extend(rh_labels)
# assert len(node_order) == node_order_size, 'Node order length is correct.'
#
# # the respective no. of regions in each cortex
# # yaml fix order change
# group_bound = [len(list(key.values())[0]) for key in labels]
# group_bound = [0] + group_bound[::-1] + group_bound
# group_boundaries = [sum(group_bound[:i+1]) for i in range(len(group_bound))]
#
# # remove the first element of group_bound
# # make label colours such that each cortex is of one colour
# group_bound.pop(0)
#
# # remove the last total sum of the list
# group_boundaries.pop()
#
# from mne.viz.circle import circular_layout
# node_angles = circular_layout(orig_labels, node_order, start_pos=90,
# group_boundaries=group_boundaries)
#
# if node_width is None:
# # widths correspond to the minimum angle between two nodes
# dist_mat = node_angles[None, :] - node_angles[:, None]
# dist_mat[np.diag_indices(n_nodes)] = 1e9
# node_width = np.min(np.abs(dist_mat))
# else:
# node_width = node_width * np.pi / 180
#
# # prepare group label positions
# group_labels = [list(lab.keys())[0] for lab in labels]
# grp_lh_labels = [name + '-lh' for name in group_labels]
# grp_rh_labels = [name + '-rh' for name in group_labels]
# all_group_labels = grp_lh_labels + grp_rh_labels
#
# # save the group order
# group_node_order = list()
# group_node_order.extend(grp_lh_labels[::-1])
# group_node_order.extend(grp_rh_labels)
# n_groups = len(group_node_order)
#
# group_node_angles = circular_layout(group_node_order, group_node_order,
# start_pos=90.)
#
# import matplotlib.pyplot as plt
# import matplotlib.path as m_path
# import matplotlib.patches as m_patches
# if fig is None:
# fig = plt.figure(figsize=figsize)
#
# # use a polar axes
# if not isinstance(subplot, tuple):
# subplot = (subplot,)
#
# ax = plt.subplot(*subplot, polar=True, facecolor='white')
#
# cortex_colors = ['m', 'b', 'y', 'c', 'r', 'g',
# 'g', 'r', 'c', 'y', 'b', 'm']
#
# label_colors = []
# for ind, rep in enumerate(group_bound):
# label_colors += [cortex_colors[ind]] * rep
# assert len(label_colors) == len(node_order), 'Number of colours do not match'
#
# # the order of the node_colors must match that of orig_labels
# # therefore below reordering is necessary
# reordered_colors = [label_colors[node_order.index(orig)]
# for orig in orig_labels]
#
# # first plot the circle showing the degree
# theta = np.deg2rad(node_angles)
# radii = np.ones(len(node_angles)) * radsize
# c = ax.scatter(theta, radii, c=reordered_colors, s=degrees * degsize,
# cmap=cmap, alpha=alpha)
#
# ax.grid(False)
# ax.set_xticks([])
# ax.set_yticks([])
# ax.spines['polar'].set_visible(False)
#
# # handle 1D and 2D connectivity information
# if con.ndim == 1:
# if indices is None:
# raise ValueError('indices has to be provided if con.ndim == 1')
# elif con.ndim == 2:
# if con.shape[0] != n_nodes or con.shape[1] != n_nodes:
# raise ValueError('con has to be 1D or a square matrix')
# # we use the lower-triangular part
# indices = np.tril_indices(n_nodes, -1)
# con = con[indices]
# else:
# raise ValueError('con has to be 1D or a square matrix')
#
# # Draw lines between connected nodes, only draw the strongest connections
# if n_lines is not None and len(con) > n_lines:
# con_thresh = np.sort(np.abs(con).ravel())[-n_lines]
# else:
# con_thresh = 0.
#
# # get the connections which we are drawing and sort by connection strength
# # this will allow us to draw the strongest connections first
# con_abs = np.abs(con)
# con_draw_idx = np.where(con_abs >= con_thresh)[0]
#
# con = con[con_draw_idx]
# con_abs = con_abs[con_draw_idx]
# indices = [ind[con_draw_idx] for ind in indices]
#
# # now sort them
# sort_idx = np.argsort(con_abs)
# con_abs = con_abs[sort_idx]
# con = con[sort_idx]
# indices = [ind[sort_idx] for ind in indices]
#
# # get the colormap
# if isinstance(cmap, str):
# colormap = plt.get_cmap(cmap)
#
# # Get vmin vmax for color scaling
# if vmin is None:
# vmin = np.min(con[np.abs(con) >= con_thresh])
# if vmax is None:
# vmax = np.max(con)
# vrange = vmax - vmin
#
# # We want to add some "noise" to the start and end position of the
# # edges: We modulate the noise with the number of connections of the
# # node and the connection strength, such that the strongest connections
# # are closer to the node center
# nodes_n_con = np.zeros((n_nodes), dtype=np.int)
# for i, j in zip(indices[0], indices[1]):
# nodes_n_con[i] += 1
# nodes_n_con[j] += 1
#
# # initialize random number generator so plot is reproducible
# rng = np.random.mtrand.RandomState(seed=0)
#
# n_con = len(indices[0])
# noise_max = 0.25 * node_width
# start_noise = rng.uniform(-noise_max, noise_max, n_con)
# end_noise = rng.uniform(-noise_max, noise_max, n_con)
#
# nodes_n_con_seen = np.zeros_like(nodes_n_con)
# for i, (start, end) in enumerate(zip(indices[0], indices[1])):
# nodes_n_con_seen[start] += 1
# nodes_n_con_seen[end] += 1
#
# start_noise[i] *= ((nodes_n_con[start] - nodes_n_con_seen[start]) /
# float(nodes_n_con[start]))
# end_noise[i] *= ((nodes_n_con[end] - nodes_n_con_seen[end]) /
# float(nodes_n_con[end]))
#
# # scale connectivity for colormap (vmin<=>0, vmax<=>1)
# con_val_scaled = (con - vmin) / vrange
#
# # convert to radians for below code
# node_angles = node_angles * np.pi / 180
# # Finally, we draw the connections
# for pos, (i, j) in enumerate(zip(indices[0], indices[1])):
# # Start point
# t0, r0 = node_angles[i], radsize
#
# # End point
# if arrow:
# # make shorter to accomodate arrowhead
# t1, r1 = node_angles[j], radsize - 1.
# else:
# t1, r1 = node_angles[j], radsize
#
# # Some noise in start and end point
# t0 += start_noise[pos]
# t1 += end_noise[pos]
#
# verts = [(t0, r0), (t0, radsize/2.), (t1, radsize/2.), (t1, r1)]
# codes = [m_path.Path.MOVETO, m_path.Path.CURVE4, m_path.Path.CURVE4,
# m_path.Path.LINETO]
# path = m_path.Path(verts, codes)
#
# color = colormap(con_val_scaled[pos])
#
# if arrow:
# # add an arrow to the patch
# patch = m_patches.FancyArrowPatch(path=path,
# arrowstyle=arrowstyle,
# fill=False, edgecolor=color,
# mutation_scale=10,
# linewidth=linewidth, alpha=1.)
# else:
# patch = m_patches.PathPatch(path, fill=False, edgecolor=color,
# linewidth=linewidth, alpha=1.)
#
# ax.add_patch(patch)
#
# # draw node labels
# if node_labels:
# angles_deg = 180 * node_angles / np.pi
# for name, angle_rad, angle_deg in zip(orig_labels, node_angles,
# angles_deg):
# if angle_deg >= 270:
# ha = 'left'
# else:
# # Flip the label, so text is always upright
# angle_deg += 180
# ha = 'right'
#
# ax.text(angle_rad, radsize + 0.2, name, size=8,
# rotation=angle_deg, rotation_mode='anchor',
# horizontalalignment=ha, verticalalignment='center',
# color='k')
#
# if show:
# plt.show()
#
# if tight_layout:
# fig.tight_layout()
#
# if out_fname:
# fig.savefig(out_fname, dpi=300.)
#
# return fig, ax
#