nperm = 1000

T_obs_bin,clusters_bin,clusters_pb_bin,H0_bin = mne.stats.spatio_temporal_cluster_test(X_bin,threshold=None,n_permutations=nperm,out_type='mask')

T_obs_ste,clusters_ste,clusters_pb_ste,H0_ste = mne.stats.spatio_temporal_cluster_test(X_ste,threshold=None,n_permutations=nperm,out_type='mask')

"""
Explanation: Analyse de clusters
End of explanation
"""


def extract_electrodes_times(clusters,clusters_pb,tmin_ind=500,tmax_ind=640,alpha=0.005,evoked = ev_bin_dev):

    ch_list_temp = []
    time_list_temp = []

    for clust,pval in zip(clusters,clusters_pb):
        if pval < alpha:

            for j,curline in enumerate(clust[tmin_ind:tmax_ind]):

                for k,el in enumerate(curline):
                    if el:                    
                        ch_list_temp.append(evoked.ch_names[k])
                        time_list_temp.append(evoked.times[j+tmin_ind])

    return np.unique(ch_list_temp),np.unique(time_list_temp)


channels_deviance_ste,times_deviance_ste=extract_electrodes_times(clusters_ste,clusters_pb_ste)

channels_deviance_bin,times_deviance_bin=extract_electrodes_times(clusters_bin,clusters_pb_bin)

print(channels_deviance_bin),print(times_deviance_bin)

print(channels_deviance_ste),print(times_deviance_ste)

times_union = np.union1d(times_deviance_bin,times_deviance_ste)

ch_union = np.unique(np.hstack([channels_deviance_bin,channels_deviance_ste]))



print(ch_union)

#Selecting channels 
epochs_bin_dev_ch = epochs_bin_dev.pick_channels(ch_union)
epochs_bin_std_ch = epochs_bin_std.pick_channels(ch_union)
epochs_ste_dev_ch = epochs_ste_dev.pick_channels(ch_union)
epochs_ste_std_ch = epochs_ste_std.pick_channels(ch_union)

X_diff = [epochs_bin_dev_ch.get_data().transpose(0, 2, 1) - epochs_bin_std_ch.get_data().transpose(0, 2, 1),
     epochs_ste_dev_ch.get_data().transpose(0, 2, 1) - epochs_ste_std_ch.get_data().transpose(0, 2, 1)]

X_diff_ste_bin = X_diff[1]-X_diff[0]


epochs_bin_dev_ch.plot_sensors(show_names=True)
plt.show()

roi = ['E117','E116','E108','E109','E151','E139','E141','E152','E110','E131','E143','E154','E142','E153','E140','E127','E118']

roi_frontal = ['E224','E223','E2','E4','E5','E6','E13','E14','E15','E20','E21','E27','E28','E30','E36','E40','E41']

len(roi_frontal),len(roi)

"""
Explanation: On récupère les channels trouvés grace a l'analyse de clusters
End of explanation
"""


from scipy.stats import ttest_1samp
from mne.stats import bonferroni_correction,fdr_correction

def ttest_amplitude(X,times_ind,ch_names,times):

    # Selecting time points  and averaging over time 
    amps = X[:,times_ind,:].mean(axis=1)
    
    T, pval = ttest_1samp(amps, 0)
    alpha = 0.05

    n_samples, n_tests= amps.shape
    threshold_uncorrected = stats.t.ppf(1.0 - alpha, n_samples - 1)

    reject_bonferroni, pval_bonferroni = bonferroni_correction(pval, alpha=alpha)
    threshold_bonferroni = stats.t.ppf(1.0 - alpha / n_tests, n_samples - 1)

    reject_fdr, pval_fdr = fdr_correction(pval, alpha=alpha, method='indep')
    
    mask_fdr = pval_fdr < 0.05
    mask_bonf = pval_bonferroni < 0.05

    print('FDR from %02f to %02f' % ((times[times_ind[0]]),times[times_ind[-1]]))
    for i,curi in enumerate(mask_fdr):
        if curi:
            print("Channel %s, T = %0.2f, p = %0.3f " % (ch_names[i], T[i],pval_fdr[i]))
    
    
    print('Bonferonni from %02f to %02f' % ((times[times_ind[0]]),times[times_ind[-1]]))
    for i,curi in enumerate(mask_bonf):
        if curi:
            print("Channel %s, T = %0.2f, p = %0.3f " % (ch_names[i], T[i],pval_bonferroni[i]))
    
    
    
    return T,pval,pval_fdr,pval_bonferronia

def ttest_amplitude_roi(X,times_ind,ch_names_roi,times):

    print(X.shape)
    # Selecting time points  and averaging over time 
    amps = X[:,times_ind,:].mean(axis=1)
    
    # averaging over channels
    amps = amps.mean(axis=1)
    
    T, pval = ttest_1samp(amps, 0)
    alpha = 0.05

    n_samples, _, n_tests= X.shape
    
    print('Uncorrected from %02f to %02f' % ((times[times_ind[0]]),times[times_ind[-1]]))
    print("T = %0.2f, p = %0.3f " % (T,pval))
    
    
    return T,pval,pval_fdr,pval_bonferroni

"""
Explanation: One sample ttest FDR corrected  (per electrode)
End of explanation
"""


toi = np.arange(0.28,0.44,0.001)
toi_index = ev_bin_dev.time_as_index(toi)

wsize = 20
wstep = 10

toi

"""
Explanation: Tests de 280 a 440, par fenetres de 20 ms avec chevauchement de 10 ms
End of explanation
"""


all_toi_indexes = []

for i in range(14):
    print(toi[10*i],toi[10*i + 20])
    cur_toi_ind = range(10*i+1,(10*i+21))
    all_toi_indexes.append(ev_bin_dev.time_as_index(toi[cur_toi_ind]))
    
print(toi[10*14],toi[10*14 + 19])
cur_toi_ind = range(10*14+1,(10*14+19))
all_toi_indexes.append(ev_bin_dev.time_as_index(toi[cur_toi_ind]))

"""
Explanation: Printing and preparing all time windows
End of explanation
"""


for cur_timewindow in all_toi_indexes:
    T,pval,pval_fdr,pval_bonferroni = ttest_amplitude(X_diff_ste_bin,cur_timewindow,epochs_bin_dev_ch.ch_names,times=epochs_bin_dev_ch.times)

"""
Explanation: Tests on each time window
End of explanation
"""


#Selecting channels 

epochs_bin_dev = _matstruc2mne_epochs(mat_bin_dev).crop(tmax=tcrop)
epochs_bin_std = _matstruc2mne_epochs(mat_bin_std).crop(tmax=tcrop)

epochs_ste_dev = _matstruc2mne_epochs(mat_ste_dev).crop(tmax=tcrop)
epochs_ste_std = _matstruc2mne_epochs(mat_ste_std).crop(tmax=tcrop)

mne.equalize_channels([epochs_bin_dev,epochs_bin_std,epochs_ste_dev,epochs_ste_std])

epochs_bin_dev_ch = epochs_bin_dev.pick_channels(roi)
epochs_bin_std_ch = epochs_bin_std.pick_channels(roi)
epochs_ste_dev_ch = epochs_ste_dev.pick_channels(roi)
epochs_ste_std_ch = epochs_ste_std.pick_channels(roi)

X_diff_roi = [epochs_bin_dev_ch.get_data().transpose(0, 2, 1) - epochs_bin_std_ch.get_data().transpose(0, 2, 1),
     epochs_ste_dev_ch.get_data().transpose(0, 2, 1) - epochs_ste_std_ch.get_data().transpose(0, 2, 1)]

X_diff_ste_bin_roi = X_diff_roi[1]-X_diff_roi[0]

for cur_timewindow in all_toi_indexes:
    T,pval,pval_fdr,pval_bonferroni = ttest_amplitude_roi(X_diff_ste_bin_roi,cur_timewindow,roi,times=epochs_bin_dev_ch.times)

grav_bin_dev = epochs_bin_dev_ch.average()
grav_bin_std = epochs_bin_std_ch.average()
grav_ste_dev = epochs_ste_dev_ch.average()
grav_ste_std = epochs_ste_std_ch.average()

evoked_bin = mne.combine_evoked([grav_bin_dev, -grav_bin_std],
                            weights='equal')

evoked_ste = mne.combine_evoked([grav_ste_dev, -grav_ste_std],
                            weights='equal')


mne.viz.plot_compare_evokeds([grav_bin_std,grav_bin_dev,grav_ste_std,grav_ste_dev],picks=[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16])
plt.show()

mne.viz.plot_compare_evokeds([evoked_bin,evoked_ste],picks=[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16])
plt.show()

"""
Explanation: On a channel subset (ROI) - average over channels
Parietal roi
End of explanation
"""


#Selecting channels 

epochs_bin_dev = _matstruc2mne_epochs(mat_bin_dev).crop(tmax=tcrop)
epochs_bin_std = _matstruc2mne_epochs(mat_bin_std).crop(tmax=tcrop)

epochs_ste_dev = _matstruc2mne_epochs(mat_ste_dev).crop(tmax=tcrop)
epochs_ste_std = _matstruc2mne_epochs(mat_ste_std).crop(tmax=tcrop)

mne.equalize_channels([epochs_bin_dev,epochs_bin_std,epochs_ste_dev,epochs_ste_std])

epochs_bin_dev_ch = epochs_bin_dev.pick_channels(roi_frontal)
epochs_bin_std_ch = epochs_bin_std.pick_channels(roi_frontal)
epochs_ste_dev_ch = epochs_ste_dev.pick_channels(roi_frontal)
epochs_ste_std_ch = epochs_ste_std.pick_channels(roi_frontal)

X_diff_roi = [epochs_bin_dev_ch.get_data().transpose(0, 2, 1) - epochs_bin_std_ch.get_data().transpose(0, 2, 1),
     epochs_ste_dev_ch.get_data().transpose(0, 2, 1) - epochs_ste_std_ch.get_data().transpose(0, 2, 1)]

X_diff_ste_bin_roi = X_diff_roi[1]-X_diff_roi[0]

for cur_timewindow in all_toi_indexes:
    T,pval,pval_fdr,pval_bonferroni = ttest_amplitude_roi(X_diff_ste_bin_roi,cur_timewindow,roi,times=epochs_bin_dev_ch.times)



grav_bin_dev = epochs_bin_dev_ch.average()
grav_bin_std = epochs_bin_std_ch.average()
grav_ste_dev = epochs_ste_dev_ch.average()
grav_ste_std = epochs_ste_std_ch.average()

evoked_bin = mne.combine_evoked([grav_bin_dev, -grav_bin_std],
                            weights='equal')

evoked_ste = mne.combine_evoked([grav_ste_dev, -grav_ste_std],
                            weights='equal')


mne.viz.plot_compare_evokeds([grav_bin_std,grav_bin_dev,grav_ste_std,grav_ste_dev],picks=[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16])
plt.show()

mne.viz.plot_compare_evokeds([evoked_bin,evoked_ste],picks=[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16])
plt.show()

mne.viz.plot_compare_evokeds?

from scipy import stats
from mne.stats import bonferroni_correction,fdr_correction
T, pval = ttest_1samp(X_diff_ste_bin, 0)
alpha = 0.05

n_samples, n_tests,_ = X_diff_ste_bin.shape
threshold_uncorrected = stats.t.ppf(1.0 - alpha, n_samples - 1)

reject_bonferroni, pval_bonferroni = bonferroni_correction(pval, alpha=alpha)
threshold_bonferroni = stats.t.ppf(1.0 - alpha / n_tests, n_samples - 1)

reject_fdr, pval_fdr = fdr_correction(pval, alpha=alpha, method='indep')
#threshold_fdr = np.min(np.abs(T)[reject_fdr])

masking_mat = pval<0.05

Tbis = np.zeros_like(T)
Tbis[masking_mat] = T[masking_mat]


plt.matshow(Tbis.T,cmap=plt.cm.RdBu_r)
plt.colorbar()
plt.show()

plt.matshow(-np.log10(pval).T)
plt.colorbar()

"""
Explanation: Frontal roi
End of explanation
"""


