In [1]:
import time
import os
import sqlite3
import pandas as pd
import numpy as np
import scipy as sp
import scipy.io as sio
import scipy.signal as sig
import plotly
plotly.offline.init_notebook_mode()
Load decomposed timeseries.
In [2]:
cf1_avg_lin_zal = pd.read_pickle('disco_decompose_coif1_average_linear_zalign.pkl')
cf1_avg_lin_zal.head()
Out[2]:
Compute intersubject phase synchronization across subjects for each frequency band.
In [77]:
def disco_ips(indata,intype,outdir,subjects):
''' Intersubject phase synchronization for group synchrony measure
indata = Input data table (id,t,v0,v1...vn)
intype = ('coif','db','sym') x ('ave','dec') x ('lin','nea','cub') x ('vec','xal','yal','zal')
outdir = Output data directory
subjects = List of subjects to include in analysis
Returns: timeseries = Group synchrony values (id,t,v) '''
function = 'disco_ips'
print(time.strftime("%m/%d/%Y"),time.strftime("%H:%M:%S"),'Running',function)
print(time.strftime("%m/%d/%Y"),time.strftime("%H:%M:%S"),'Subjects =',subjects)
pairs = len(subjects) * (len(subjects)-1) / 2 # Number of subject pairs
levels = indata.drop(['id','t'],axis=1).shape[1] # Number of frequency bands
IPS = {} # Intersubject phase synchronization for each frequency band (samples x levels)
for L in range(levels):
print(time.strftime("%m/%d/%Y"),time.strftime("%H:%M:%S"),'Level =',L)
# ADAPTED FROM CALCULATEPHASESYNCH.M (ISC TOOLBOX)
AD = [] # Angular distance at each time point for each subject pair
for s1 in range(len(subjects)):
# Hilbert transform to obtain corresponding analytic signal
t1 = sig.hilbert(indata[indata.id==subjects[s1]]['L'+str(L).zfill(2)])
for s2 in range(len(subjects)):
if s2 > s1:
# print(time.strftime("%m/%d/%Y"),time.strftime("%H:%M:%S"),subjects[s1],subjects[s2])
t2 = sig.hilbert(indata[indata.id==subjects[s2]]['L'+str(L).zfill(2)])
AD.append(np.angle(t1)-np.angle(t2))
if pairs == 1:
absAD = np.absolute(AD) # Absolute angular distance (Dissimilarity measure)
else:
absAD = np.sum(np.absolute(AD),0) / pairs # Average of all subject-pairwise absolute angular distances
IPS['L'+str(L).zfill(2)] = 1 - absAD / np.pi # Normalized version of (average) absolute angular distance(s)
IPS = pd.DataFrame(IPS)
T = pd.DataFrame({'t':indata[indata.id==subjects[s1]].t})
timeseries = pd.concat([T,IPS],axis=1)
outfile = function+'_'+intype+'.pkl'
print(time.strftime("%m/%d/%Y"),time.strftime("%H:%M:%S"),'Saving',outfile)
timeseries.to_pickle(os.path.join('.',outdir,outfile))
print(time.strftime("%m/%d/%Y"),time.strftime("%H:%M:%S"),'Done')
return timeseries
In [78]:
ips_cf1_avg_lin_zal = disco_ips(cf1_avg_lin_zal,'coif1_average_linear_zalign','',['1b','2b','3b','4b','5b','6b'])
In [80]:
ips_cf1_avg_lin_zal.head()
Out[80]:
Heat map of intersubject phase synchronization at each frequency band.
In [2]:
# Load synchrony timeseries
df = sio.loadmat('matlabfiles/disco_ips_coif1_mean_linear_zalign.mat',variable_names='IPS')
df = pd.DataFrame.transpose(pd.DataFrame(df['IPS']))
In [3]:
# Load condition times and names
timedb = sqlite3.connect('disco_millisecond_conditions.db')
query = "SELECT * FROM conditions;"
timedf = pd.read_sql_query(query,timedb)
timedb.close()
In [ ]:
samplingrate = 62.50 # Sampling rate in HZ (1000ms / 16ms)
minutes = np.arange(0,df.shape[0])/samplingrate/60 # Time in minutes (X)
frequency = [str('%.3f' %i) for i in samplingrate/2**np.arange(0,df.shape[1])] # Frequencies of levels (Y)
# Condition lines
conditions = []
for t in range(0+1,len(timedf)-1): # Not including first and last time points
X = (timedf['msec'][t]-timedf['msec'][0])/1000/60 # Convert to minutes from first time point
L = dict(x0=X,x1=X,y0=0-.5,y1=len(frequency)-.5,type='line',line=dict(width=1.5,dash='dot'))
conditions.append(L)
# Synchrony values (Z)
values = []
for f in range(len(frequency)):
values.append(df[f].tolist())
data = [plotly.graph_objs.Heatmap(
z = values,
x = minutes,
y = frequency,
colorscale = 'Viridis')]
layout = plotly.graph_objs.Layout(
title = 'Intersubject Phase Sychrony (Pipeline: C=Z-align, I=Linear, D=Average, W=Coiflet/Short)',
xaxis = dict(title='Time (minutes)',hoverformat='.3f',ticks='',zeroline=False),
yaxis = dict(title='Frequency (Hz)',hoverformat='.3f',ticks='',autorange='reversed',type='category'),
shapes = conditions)
print(''); print(timedf['name'][0:len(timedf)-1].tolist()) # Condition names
fig = plotly.graph_objs.Figure(data=data,layout=layout)
plotly.offline.iplot(fig)