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 cluster phase method across subjects for each frequency band.
In [5]:
def disco_cpm(indata,intype,outdir,subjects):
''' Cluster phase method 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_cpm'
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)
levels = indata.drop(['id','t'],axis=1).shape[1] # Number of frequency bands
CPM = {} # Group rho timeseries 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)
# Create table of subjects x samples for current level
SUBJ = pd.DataFrame({})
for S in subjects:
temp = indata[indata.id==S]['L'+str(L).zfill(2)]
SUBJ = pd.concat([SUBJ,temp],axis=1,ignore_index=True)
SUBJ = pd.DataFrame.transpose(SUBJ)
# ADAPTED FROM CLUSTERPHASE_DO.M
# Compute phase for each timeseries using Hilbert transform
hilbert = sig.hilbert(SUBJ)
arctan2 = np.arctan2(np.real(hilbert),np.imag(hilbert))
TSphase = np.unwrap(arctan2)
# Compute mean running or cluster phase
znumber = np.exp(1j*TSphase)
ztotal = np.mean(znumber,axis=0)
zangle = np.angle(ztotal)
clusterphase = np.unwrap(zangle)
# Compute relative phases between timeseries phase and cluster phase
znumber = np.exp(1j*(TSphase-clusterphase))
ztotal = np.mean(znumber,axis=1)
TSRPM = np.angle(ztotal)
TSRPM = np.tile(TSRPM,(len(clusterphase),1)).transpose() # So valid dimensions
# Compute cluster amplitute in rotation frame
znumber = np.exp(1j*(TSphase-clusterphase-TSRPM))
ztotal = np.mean(znumber,axis=0)
TSrhoGRP = np.absolute(ztotal)
CPM['L'+str(L).zfill(2)] = TSrhoGRP
CPM = pd.DataFrame(CPM)
T = pd.DataFrame({'t':indata[indata.id==subjects[0]].t})
timeseries = pd.concat([T,CPM],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 [6]:
cpm_cf1_avg_lin_zal = disco_cpm(cf1_avg_lin_zal,'coif1_average_linear_zalign','',['1b','2b','3b','4b','5b','6b'])
In [7]:
cpm_cf1_avg_lin_zal.head()
Out[7]:
Heat map of cluster phase method at each frequency band.
In [2]:
# Load synchrony timeseries
df = sio.loadmat('matlabfiles/disco_cpm_coif1_mean_linear_zalign.mat',variable_names='CPM')
df = pd.DataFrame.transpose(pd.DataFrame(df['CPM']))
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 = 'Cluster Phase Synchronization (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)