figure_decoding.py

#%matplotlib qt
import numpy as np
from config_AVTest import MEG_data_path,group_name,Ids,task,day
import os.path as op
import matplotlib.pyplot as plt
import seaborn as sns
import mne
from EZ_stats import calculate_significance, plot_generalization_matrix


# load all scores from all subjects
all_GE_congruency_scores = []
all_FI_congruency_scores = []
all_VV_scores = []

all_evks_reactivation_map = {'FI': [], 'GE': [], 'VV': []}
all_scores={'GE_GE': [], 'FI_FI': [], 'GE_FI': [], 'FI_GE': [], 'VV_VV': [], 'score_diff_FI_FI_GE_GE': [], 'score_diff_FI_GE_GE_FI':[]}


for subject_id in Ids:
    subject = group_name+"%d" % subject_id
    fname=op.join(MEG_data_path,subject,task+'_%d'%(day+subject_id)+'_tsss_mc.fif')
    GE_congruency_score = np.load(fname.replace('tsss_mc.fif','TD_GE.npy'))
    FI_congruency_score = np.load(fname.replace('tsss_mc.fif','TD_FI.npy'))
    VV_score = np.load(fname.replace('tsss_mc.fif','TD_VV.npy'))

    all_GE_congruency_scores.append(GE_congruency_score)
    all_FI_congruency_scores.append(FI_congruency_score)
    all_VV_scores.append(VV_score)

    # reactivation maps
    FI_evoked_reactivation_map = mne.read_evokeds(fname.replace('tsss_mc.fif','FI_reactivation_map-ave.fif'))[0]
    GE_evoked_reactivation_map = mne.read_evokeds(fname.replace('tsss_mc.fif','GE_reactivation_map-ave.fif'))[0]
    VV_evoked_reactivation_map = mne.read_evokeds(fname.replace('tsss_mc.fif','VV_reactivation_map-ave.fif'))[0]
    all_evks_reactivation_map['FI'].append(FI_evoked_reactivation_map)
    all_evks_reactivation_map['GE'].append(GE_evoked_reactivation_map)
    all_evks_reactivation_map['VV'].append(VV_evoked_reactivation_map)


    # TG
    scores_GE_GE=np.load(fname.replace('tsss_mc.fif','TG_GE_GE.npy'))
    scores_FI_FI=np.load(fname.replace('tsss_mc.fif','TG_FI_FI.npy'))
    scores_GE_FI=np.load(fname.replace('tsss_mc.fif','TG_GE_FI.npy'))
    scores_FI_GE=np.load(fname.replace('tsss_mc.fif','TG_FI_GE.npy'))
    score_VV_VV=np.load(fname.replace('tsss_mc.fif','TG_VV_VV.npy'))
    score_diff_FI_FI_GE_GE = scores_FI_FI - scores_GE_GE
    score_diff_FI_GE_GE_FI = scores_FI_GE - scores_GE_FI

    all_scores['GE_GE'].append(scores_GE_GE)
    all_scores['FI_FI'].append(scores_FI_FI)
    all_scores['GE_FI'].append(scores_GE_FI)
    all_scores['FI_GE'].append(scores_FI_GE)
    all_scores['VV_VV'].append(score_VV_VV)
    all_scores['score_diff_FI_FI_GE_GE'].append(score_diff_FI_FI_GE_GE)
    all_scores['score_diff_FI_GE_GE_FI'].append(score_diff_FI_GE_GE_FI)

all_GE_congruency_scores=np.array(all_GE_congruency_scores)
all_FI_congruency_scores=np.array(all_FI_congruency_scores)
all_VV_scores=np.array(all_VV_scores)

mean_GE_congruency_scores = np.mean(all_GE_congruency_scores, axis=0)
mean_FI_congruency_scores = np.mean(all_FI_congruency_scores, axis=0)
mean_VV_scores = np.mean(all_VV_scores, axis=0)

#standard error of the mean
se_GE_congruency_scores = np.std(all_GE_congruency_scores, axis=0) / np.sqrt(len(Ids))
se_FI_congruency_scores = np.std(all_FI_congruency_scores, axis=0) / np.sqrt(len(Ids))
se_VV_scores = np.std(all_VV_scores, axis=0) / np.sqrt(len(Ids))

times = np.linspace(-0.2, 1.0, GE_congruency_score.shape[0])

corr_method='none'
significant_GE = calculate_significance(all_GE_congruency_scores, times, [0.2,1.0], corr_method)
significant_FI = calculate_significance(all_FI_congruency_scores, times, [0.2,1.0], corr_method)
significant_VV = calculate_significance(all_VV_scores, times, [0.0,1.0], corr_method)
significant_diff = calculate_significance(all_FI_congruency_scores - all_GE_congruency_scores, times, [0.3,1.0], corr_method, 0.0,0)

corr_method='tfce'
significant_GE_corr = calculate_significance(all_GE_congruency_scores, times, [0.2,1.0], corr_method)
significant_FI_corr = calculate_significance(all_FI_congruency_scores, times, [0.2,1.0], corr_method)
significant_VV_corr = calculate_significance(all_VV_scores, times, [0.0,1.0], corr_method)
significant_diff_corr = calculate_significance(all_FI_congruency_scores - all_GE_congruency_scores, times, [0.2,1.0], corr_method, 0.0,0)

colors = sns.color_palette()# "colorblind"
# plot the results for GE Congruency
fig, ax = plt.subplots(figsize=(8, 4))

# plot the results for FI Congruency
ax.plot(times, mean_FI_congruency_scores, color=colors[0], label='FI Congruency', linewidth=2)
ax.fill_between(times, (mean_FI_congruency_scores - se_FI_congruency_scores), (mean_FI_congruency_scores + se_FI_congruency_scores), color=colors[0], alpha=.2)
ax.axhline(.5, color='k', linestyle='--', label='chance')
ax.plot(times[significant_FI], np.ones_like(times[significant_FI])*0.47, 'o', color=colors[0], label='Significant FI', markersize=1)
ax.plot(times[significant_FI_corr], np.ones_like(times[significant_FI_corr])*0.47, 'o', color=colors[0], label='Significant GE (TFCE)', markersize=4)
ax.set_xlabel('Times')
ax.set_ylabel('AUC')  # Area Under the Curve
ax.legend()
ax.axvline(.0, color='k', linestyle='-')
ax.set_title('Sensor space decoding - FI Congruency')

# plot the results for GE Congruency
ax.plot(times, mean_GE_congruency_scores, color=colors[1], label='GE Congruency', linewidth=2)
ax.fill_between(times, (mean_GE_congruency_scores - se_GE_congruency_scores), (mean_GE_congruency_scores + se_GE_congruency_scores), color=colors[1], alpha=.2)
ax.axhline(.5, color='k', linestyle='--', label='chance')
ax.plot(times[significant_GE], np.ones_like(times[significant_GE])*0.465, 'o', color=colors[1], label='Significant GE', markersize=1)
ax.plot(times[significant_GE_corr], np.ones_like(times[significant_GE_corr])*0.465, 'o', color=colors[1], label='Significant GE (TFCE)', markersize=4)
ax.set_xlabel('Times')
ax.set_ylabel('AUC')  # Area Under the Curve
ax.legend()
ax.axvline(.0, color='k', linestyle='-')
ax.set_title('Sensor space decoding - GE Congruency')

# plot the results for Difference

#fig, ax = plt.subplots(figsize=(9, 3))
# ax.plot(times, mean_GE_congruency_scores - mean_FI_congruency_scores, color=colors[2], label='Difference')
# ax.fill_between(times, (mean_GE_congruency_scores - mean_FI_congruency_scores - se_GE_congruency_scores + se_FI_congruency_scores), (mean_GE_congruency_scores - mean_FI_congruency_scores + se_GE_congruency_scores - se_FI_congruency_scores), color=colors[2], alpha=.1)
# ax.axhline(.0, color='k', linestyle='--', label='no difference')
ax.plot(times[significant_diff], np.ones_like(times[significant_diff])*0.46, 'o', color=colors[2], label='Significant Difference', markersize=1)
ax.plot(times[significant_diff_corr], np.ones_like(times[significant_diff_corr])*0.46, 'o', color=colors[2], label='Significant Difference (TFCE)', markersize=4)
ax.set_xlabel('Times')
ax.set_ylabel('Difference in AUC')  # Area Under the Curve
ax.legend()
ax.axvline(.0, color='k', linestyle='-')
ax.set_title('Sensor space decoding - Difference')

# plot the results for VV
fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(times, mean_VV_scores, color=colors[3], label='VV', linewidth=2)
ax.fill_between(times, (mean_VV_scores - se_VV_scores), (mean_VV_scores + se_VV_scores), color=colors[3], alpha=.2)
ax.axhline(.5, color='k', linestyle='--', label='chance')
ax.plot(times[significant_VV], np.ones_like(times[significant_VV])*0.47, 'o', color=colors[3], label='Significant VV', markersize=1)
ax.plot(times[significant_VV_corr], np.ones_like(times[significant_VV_corr])*0.47, 'o', color=colors[3], label='Significant VV (TFCE)', markersize=4)
ax.set_xlabel('Times')
ax.set_ylabel('AUC')  # Area Under the Curve
ax.legend()
ax.axvline(.0, color='k', linestyle='-')
ax.set_title('Sensor space decoding - VV')

# Print the TFCE corrected significant time windows
print("TFCE Significant Time Windows:")
print("GE:", times[significant_GE_corr][[0, -1]] if any(significant_GE_corr) else "No significant time points")
print("FI:", times[significant_FI_corr][[0, -1]] if any(significant_FI_corr) else "No significant time points")
print("VV:", times[significant_VV_corr][[0, -1]] if any(significant_VV_corr) else "No significant time points")
print("Difference:", times[significant_diff_corr][[0, -1]] if any(significant_diff_corr) else "No significant time points")



#grand average sensor space
FI_evoked_reactivation_map = mne.grand_average(all_evks_reactivation_map['FI'])
GE_evoked_reactivation_map = mne.grand_average(all_evks_reactivation_map['GE'])
VV_evoked_reactivation_map = mne.grand_average(all_evks_reactivation_map['VV'])

# plot the reactivation maps
FI_evoked_reactivation_map.plot_topomap(times=[0.1,0.3,0.5,0.7,0.9],vlim=(0,30),average=0.201,ch_type='grad')
#plt.savefig('FI_activation_topomap.pdf')
GE_evoked_reactivation_map.plot_topomap(times=[0.1,0.3,0.5,0.7,0.9],vlim=(0,30),average=0.201,ch_type='grad')
#plt.savefig('GE_activation_topomap.pdf')
VV_evoked_reactivation_map.plot_topomap(times=[0.1,0.3,0.5,0.7,0.9],vlim=(0,25),average=0.201,ch_type='grad')
#plt.savefig('VV_activation_topomap.pdf')



times = np.linspace(-0.2, 1.0, all_scores['GE_GE'][0].shape[0])
# plot the TG results
plot_generalization_matrix(times,[0.0,1.0], all_scores['VV_VV'], 'VV', 'VV',True,'tfce',0,chance=0.5,vmin=0.4,vmax=0.6)

p_val_corrected_GE_GE=plot_generalization_matrix(times,[0.2,1.0], all_scores['GE_GE'], 'GE', 'GE',True,'tfce',0,chance=0.5,vmin=0.44,vmax=0.56)
p_val_corrected_FI_FI=plot_generalization_matrix(times,[0.2,1.0], all_scores['FI_FI'], 'FI', 'FI',True,'tfce',0,chance=0.5,vmin=0.44,vmax=0.56)

p_val_corrected_GE_GE[np.isnan(p_val_corrected_GE_GE)] = 1
p_val_corrected_FI_FI[np.isnan(p_val_corrected_FI_FI)] = 1
exclud_mask =(p_val_corrected_GE_GE > 0.05) & (p_val_corrected_FI_FI > 0.05)
plot_generalization_matrix(times,[0.2,1.0], all_scores['score_diff_FI_FI_GE_GE'], 'FI_FI', 'GE_GE',True,'tfce',0,chance=0.0, tfce_stat_mask=exclud_mask)

p_val_corrected_FI_GE=plot_generalization_matrix(times,[0.2,1.0], all_scores['FI_GE'], 'FI', 'GE',True,'tfce',0,chance=0.5,vmin=0.44,vmax=0.56)
p_val_corrected_GE_FI=plot_generalization_matrix(times,[0.2,1.0], all_scores['GE_FI'], 'GE', 'FI',True,'tfce',0,chance=0.5,vmin=0.44,vmax=0.56)

p_val_corrected_FI_GE[np.isnan(p_val_corrected_FI_GE)] = 1
p_val_corrected_GE_FI[np.isnan(p_val_corrected_GE_FI)] = 1
exclud_mask=(p_val_corrected_FI_GE > 0.05) & (p_val_corrected_GE_FI > 0.05)
plot_generalization_matrix(times,[0.2,1.0], all_scores['score_diff_FI_GE_GE_FI'], 'FI_GE', 'GE_FI',True,'tfce',0,chance=0.0, tfce_stat_mask=exclud_mask)