PLV and BI calculation from epochs data in a label.

gistfile1.py

##################################################
##DESCRIPTION##
#The script stat_BI_distribRecording.py loads .mat files saved as in the script BifurcInd_reponse_LG.py. BI, PLV and power distribution used in resampling statistics (Busch_2009) are registered to be used in the second scrit stat_BI_resamplingTest_fromMat.py
##################################################

import os

from copy import deepcopy
import mne
from mne.minimum_norm import apply_inverse, make_inverse_operator,apply_inverse_epochs

from mne.event import merge_events
from scipy.signal import hilbert
from mne.viz import plot_topo
from mne.time_frequency.tfr import cwt,morlet,induced_power
from mne.minimum_norm import source_induced_power
from mne.stats import bonferroni_correction, fdr_correction
import cmath
from math import sqrt, atan2, pi, floor,exp

import numpy as np
import matplotlib.pyplot as plt
import scipy
import scipy.io
from scipy.io import loadmat
import pylab as pl
import random
from scipy import linalg

################################### Function definition ########################################
 
# Average data using a ponderation filter found by svd (singulate value decomposition) and return the averaged evoked data and the induced averaged induced data. Use the spatial filter created on the total data on each condition.
# data(trials,vertices,times)
# d(trials, times)
# d_induced(trials, times)
def average_svd(data_tot,data_cond1,data_cond2):

	# Compute the ponderation filter from total data
	avg = data_tot.mean(0) # average across trials
	U,s,V = linalg.svd(avg) # apply svd
	spatial_filter = U[:,0] # ponderation filter

	# Apply the filter to each condition
	d_tot = np.array([np.dot(spatial_filter,e) for e in data_tot])
	d_cond1 = np.array([np.dot(spatial_filter,e) for e in data_cond1])
	d_cond2 = np.array([np.dot(spatial_filter,e) for e in data_cond2])

	# Compute induced data
	d_tot_mean = np.mean(d_tot, axis=0)[None,:]
	d_tot_ind = d_tot - d_tot_mean # substract evoked response
	d_cond1_mean = np.mean(d_cond1, axis=0)[None,:]
	d_cond1_ind = d_cond1 - d_cond1_mean # substract evoked response
	d_cond2_mean = np.mean(d_cond2, axis=0)[None,:]
	d_cond2_ind = d_cond2 - d_cond2_mean # substract evoked response

	return (d_tot, d_cond1, d_cond2, d_tot_ind, d_cond1_ind, d_cond2_ind )
################################### Variable definition #########################################

sujet_ordre = ['pf120155','pe110338','cj100142','jm100042', 'jm100109','sb120316', 'tk130502', 'sl130503', 'rl130571', 'bd120417','rb130313', 'mp140019']
sujetinit_ordre =['pf','pe', 'cj','jm42','jm09','sb','tk', 'sl', 'rl','bd','rb','mp14']
PSS =  [-0.039, -0.023,  -0.027, -0.049, -0.083, -0.003, -0.043, -0.012, -0.012, +0.043, +0.023, +0.033] 
JND1 = [-0.118, -0.101,  -0.115, -0.077, -0.112,-0.075, -0.089, -0.097, -0.116, -0.121, -0.067, -0.132]
JND2 = [+0.040, +0.055,  +0.060, -0.021, -0.054,+0.070, +0.003, +0.073, +0.093, +0.208, +0.112, +0.198]

frequency_band = ['delta','theta','alpha','lbeta','hbeta','gamma']
frequency_list = [2., 5.5, 10., 16., 25., 40.]
n_cycles_list = [3., 3.4, 3.8, 5.3, 4.6, 5]

decim = 4
Fs = 1000

tmin, tmax =-3.3, 2
tmin_stat, tmax_stat =-1, 0 # To change: temporal window on the plot
dur = floor((tmax-tmin)*1000/decim)+1
dur_stat = floor((tmax_stat-tmin_stat)*1000/decim)+1
times = np.arange(int(tmin*1000), int(tmax*1000+1),decim)


wdir = '/neurospin/meg/meg_tmp/PhaseTime_Anne_2013/'
datapath_stc = wdir + 'stc/'
datapath_mat = wdir + 'data_mat/Distributions_resamplingTest/'

cond = 'pss'
hemi ='rh'
label_name = 'Brodmann.42-' + hemi

######################################## MAIN ################################################

#Initialize array for BI distribution per subject
Nloop = 500
BI = np.empty((len(sujet_ordre),Nloop,len(frequency_list),dur))
BI_toTest = np.empty((len(sujet_ordre),len(frequency_list),dur))

for s, subj in enumerate(sujet_ordre):
	print subj
	data_path = wdir + 'maxfilter/'+subj+'/'
	subj_init = sujetinit_ordre[s]

	### Data of each trials in the label of interest ###

	# Read data from epochs created by stc_BIpowPLV_fband.py (equalized trials)
	epochs_A = mne.read_epochs(wdir + 'epochs/' +  subj + '_equalizedEpochs_' + cond + '_firstStimLocked_Afirst')
	epochs_V = mne.read_epochs(wdir + 'epochs/' +  subj + '_equalizedEpochs_' + cond + '_firstStimLocked_Vfirst')
	epochs_tot = mne.read_epochs(wdir + 'epochs/' +  subj + '_equalizedEpochs_' + cond + '_firstStimLocked_tot')

	# compute noise covariance from empty room data
	emptyroom_raw = mne.fiff.Raw('/neurospin/meg/meg_tmp/PhaseTime_Anne_2013/maxfilter/'  +subj+'/'+ subj_init +'_empty_sss.fif')  
	noise_cov = mne.compute_raw_data_covariance(emptyroom_raw)

	# compute dSPM solution
	fname_fwd = data_path + subj + '_phase1_trans_sss_filt140_raw-ico5-fwd.fif'
	forward = mne.read_forward_solution(fname_fwd, surf_ori=True)

	# create inverse operator
	inverse_operator = make_inverse_operator(epochs_tot.info, forward, noise_cov, loose=0.4, depth=0.8)

	# Read label of interest
	datapath_label = wdir + 'freesurfer_segmentation/' + subj + '/label/'
	label = mne.read_label(datapath_label + label_name + '.label')

	# Compute inverse solution
	snr = 3.0
	lambda2 = 1.0 / snr ** 2
	dSPM = True
	stcs_A = apply_inverse_epochs(epochs_A, inverse_operator, lambda2, method='dSPM', label =label ,nave=1, pick_ori="normal")
	stcs_V = apply_inverse_epochs(epochs_V, inverse_operator, lambda2, method='dSPM', label =label ,nave=1, pick_ori="normal")
	stcs_tot = apply_inverse_epochs(epochs_tot, inverse_operator, lambda2, method='dSPM', label =label ,nave=1, pick_ori="normal")

	# Store data
	tmp_tot = np.empty((len(stcs_tot),stcs_tot[0].shape[0],stcs_tot[0].shape[1]))
	for i in range(len(stcs_tot)):
		tmp_tot[i,:,:] = stcs_tot[i].data

	tmp_A = np.empty((len(stcs_A),stcs_A[0].shape[0],stcs_A[0].shape[1]))
	tmp_V = np.empty((len(stcs_V),stcs_V[0].shape[0],stcs_V[0].shape[1]))
	for i in range(len(stcs_A)):
		tmp_A[i,:,:] = stcs_A[i].data
		tmp_V[i,:,:] = stcs_V[i].data

	# Average via svd on the vertices (selected label)	
	data_tot, data_A, data_V, data_tot_ind, data_A_ind, data_V_ind = average_svd(tmp_tot, tmp_A, tmp_V)

	### Compute BI to test ###
	# Rque: the BI in stc_BIpowPLV_fband.py is calculated on each source then averaged so should be slightly different from this BI calculated from svd-averaged trials on one label.
	pow_A, PLV_A = induced_power(data_A[:,None,:], Fs=Fs/decim, frequencies=frequency_list, n_cycles=n_cycles_list, n_jobs=8, use_fft=False, decim=1, zero_mean=True)
	pow_V, PLV_V = induced_power(data_V[:,None,:], Fs=Fs/decim, frequencies=frequency_list, n_cycles=n_cycles_list, n_jobs=8, use_fft=False, decim=1, zero_mean=True)
	pow_tot, PLV_tot = induced_power(data_tot[:,None,:], Fs=Fs/decim, frequencies=frequency_list, n_cycles=n_cycles_list, n_jobs=8, use_fft=False, decim=1, zero_mean=True)

    # Compute bifurcation index
	BI_toTest[s] = (PLV_A[0] - PLV_tot[0])*(PLV_V[0] - PLV_tot[0])