clearvars -except ppp
SC_group_level_stat_input.subnum=0;
for sub = {'IR29','IR35','IR51','IR54', 'IR57' ,'IR75','IR84'}    
    sub=char(sub)
    clearvars -except sub SC_group_level_stat_input ppp

    cd(['\\mycloudpr4100\ivan\HaoxinAnalyses\' sub])
    load([sub '_log_EMOPSData.mat'])
    load([sub '_HaoxinROI.mat'])
    try
        load([sub '_evnt_clean.mat'])
    catch
        ok_epochs_test=1:290;
        ok_epochs_study=1:148;
    end

    try
        MTL_roi.EC_Stan;
    catch
        MTL_roi.EC_Stan=[];
    end
    try
        MTL_roi.PRC_Stan;
    catch
        MTL_roi.PRC_Stan=[];
    end
    try
        MTL_roi.PHC_Stan;
    catch
        MTL_roi.PHC_Stan=[];
    end
    try
        MTL_roi.Sub_Stan;
    catch
        MTL_roi.Sub_Stan=[];
    end
    try
        MTL_roi.AMY_Stan;
    catch
        MTL_roi.AMY_Stan=[];
    end
    %     if sub=='IR75'
    %         disp('here we go')
    %         figure
    %     end
    chan_ROI = unique([MTL_roi.CA3_Stan MTL_roi.CA1_Stan]); % MTL_roi.EC_Stan MTL_roi.PRC_Stan MTL_roi.PHC_Stan  MTL_roi.CA3_Stan MTL_roi.CA1_Stan  MTL_roi.AMY_Stan

    if ~isempty(chan_ROI)
        SC_group_level_stat_input.subnum=SC_group_level_stat_input.subnum+1;
        SC_group_level_stat_input.channum(SC_group_level_stat_input.subnum) = 0;
        %
        %         cd(['/Users/jacklin/Desktop/JLab/EMOP/Data/' sub '/Similarity'])
        %         Sim_sub_level_Output_R = load(['Sim_sub_level_LFA_LCR_w_shuffle1000_Output_' sub '.mat'], 'Sim_sub_level_Output', 'trial_condition', 'time_middle_ecd_rtv');
        %         time_middle_ecd_rtv = Sim_sub_level_Output_R.time_middle_ecd_rtv;
        %         Sim_sub_level_Output_R = Sim_sub_level_Output_R.Sim_sub_level_Output;
        %          R_map_avg=[];
        %         for ch = chan_ROI
        %             R_map_avg = cat(4,R_map_avg,Sim_sub_level_Output_R.true_sameitem_r{find(ch==Sim_sub_level_Output_R.chan)});
        %         end
        %         tROI = find(time_middle_ecd_rtv>=0.5 & time_middle_ecd_rtv<=1.5); % AMY
        %         [m,which_ch] = max(mean(R_map_avg(tROI,:,:,:),[1 2 3],'omitnan'));

        for ch = chan_ROI
            ch
            tic
            SC_group_level_stat_input.channum(SC_group_level_stat_input.subnum) = SC_group_level_stat_input.channum(SC_group_level_stat_input.subnum) +1;
            cd(['\\mycloudpr4100\ivan\HaoxinAnalyses\' sub ])
            load(['TFR_EcdRtv_ForRippletrigerRSA_ch' num2str(ch)],  'Ecd_trial_list', 'Rtv_trial_list', 'same_pair_TFR_ecd_bined', 'same_pair_zTFR_ecd_bined', 'time_middle_ecd', 'trial_condition', 'wavefVec')
            not_good_trial = find(~ismember(Ecd_trial_list,ok_epochs_study) | ~ismember(Rtv_trial_list,ok_epochs_test));
            t_RSA_ecd = time_middle_ecd(find(time_middle_ecd >= 0 & time_middle_ecd<=2 )) ;
            t_RSA_ecd_rt = time_middle_ecd(find(time_middle_ecd >= 0 & time_middle_ecd<=5 )) ;  % (time_middle_ecd >= 2 & time_middle_ecd<=5 )
            TFR1 = same_pair_zTFR_ecd_bined(:,find(time_middle_ecd >= 0 & time_middle_ecd<=2 ),: ) ;
            TFR2 = same_pair_zTFR_ecd_bined(:, find(time_middle_ecd >= 0 & time_middle_ecd<=5 ) ,: ); % (time_middle_ecd >= 2 & time_middle_ecd<=5 )
            
            Similarity_same_pair_ecdwresponse = nan(length(t_RSA_ecd),length(t_RSA_ecd_rt),147) ;
            Similarity_next_pair_ecdwresponse = nan(length(t_RSA_ecd),length(t_RSA_ecd_rt),147) ;
            for trial_num = 1:147
                Similarity_same_pair_ecdwresponse(:,:,trial_num) = atanh(corr(TFR1(:,:,trial_num),TFR2(:,:,trial_num),'Type','Spearman')) ;
%                 if trial_num<147
                    for ntriallist = setdiff(1:147,trial_num) 
                    a(:,:,ntriallist) =  atanh(corr(TFR1(:,:,trial_num),TFR2(:,:,ntriallist),'Type','Spearman')) ;
                    end
                     a(:,:,trial_num)  = nan; 
                     
                Similarity_next_pair_ecdwresponse(:,:,trial_num) = mean(a,3,'omitnan') ;
%                 else
%                     Similarity_next_pair_ecdwresponse(:,:,trial_num) = nan;
%                 end
            end
            

            SC_group_level_stat_input.Similarity_same_pair_ecdwresponse{SC_group_level_stat_input.subnum}(:,:,:,SC_group_level_stat_input.channum(SC_group_level_stat_input.subnum)) = Similarity_same_pair_ecdwresponse;
            SC_group_level_stat_input.Similarity_next_pair_ecdwresponse{SC_group_level_stat_input.subnum}(:,:,:,SC_group_level_stat_input.channum(SC_group_level_stat_input.subnum)) = Similarity_next_pair_ecdwresponse;
            SC_group_level_stat_input.trial_condition{SC_group_level_stat_input.subnum}(:,SC_group_level_stat_input.channum(SC_group_level_stat_input.subnum))=trial_condition;
            SC_group_level_stat_input.t_RSA_ecd = t_RSA_ecd ;
            SC_group_level_stat_input.t_RSA_ecd_rt = t_RSA_ecd_rt;
            toc
            clearvars -except SC_group_level_stat_input ch ok_epochs_study ok_epochs_test sub chan_ROI tt
        end
    end
end

clearvars -except sub SC_group_level_stat_input

%% Include the ripple event
cd('\\mycloudpr4100\ivan\HaoxinAnalyses')
load('R01_ripple_stats_output','output')
output(14,8) = 1;
cd('C:\Users\elvis\Dropbox\EMOP\BehaviourResults')
StimandArousal = importStim_and_Arousal('Stim_and_Arousal.csv', 1, 192);
load('Stim_Valence.mat')


group_level_stat_input.subnum=0;
group_level_stat_input.bhv=[];
group_level_stat_input.Visitation=[];
group_level_stat_input.Visitation_bhv=[];
group_level_stat_input.ripple_summary = [];
group_level_stat_input.ripple_presence = [];
group_level_stat_input.ripple_rate = [];
for sub = {'IR29','IR35','IR51','IR54' ,'IR57', 'IR75','IR84'}    
    sub=char(sub)
    clearvars -except SC_group_level_stat_input ch ok_epochs_study ok_epochs_test sub chan_ROI tt output group_level_stat_input StimandArousal Stim_Valence
    
    cd(['\\mycloudpr4100\ivan\HaoxinAnalyses\' sub])
    load([sub '_log_EMOPSData.mat'])
    load([sub '_evnt_clean.mat'])
    load([sub 'Ecd_trial_num.mat'])
    group_level_stat_input.subnum = group_level_stat_input.subnum+1;
    EMOPSData.testVar(:,end+1) = group_level_stat_input.subnum ;
    EMOPSData.testVar(:,end+1) = 0;
    EMOPSData.testVar(ok_epochs_test,end) = 1;
    EMOPSData.testVar(:,end+1) = Ecd_trial_num;
    bhv = EMOPSData.testVar;
    
    img = EMOPSData.testImage;
    aro = nan(290,1);
    val = nan(290,1);
    sti_list = table2array(StimandArousal(:,1));
    for nn=1:290
        sti_curt = str2num(img{nn}(1:5)) ;
        if ismember(sti_curt,sti_list)
            aro(nn) = table2array(StimandArousal(find(sti_curt==sti_list),2)) ;
        end
        if ismember(sti_curt,Stim_Valence(:,1))
            val(nn) = Stim_Valence(sti_curt==Stim_Valence(:,1),2) ;
        end
    end
    EMOPSData.testVar(:,end+1) = aro ;
    EMOPSData.testVar(:,end+1) = val ;
    
    
    load([sub '_HaoxinROI.mat'])
    %   MTL_roi.CA3_Stan MTL_roi.CA1_Stan   MTL_roi.AMY_Stan
    chan_ROI = unique([MTL_roi.CA3_Stan MTL_roi.CA1_Stan]);
    
    load(['PCA_tROI_Overlay_output_d6_Nov6_ch' num2str(chan_ROI(1))],'Ecd_trial_list','Rtv_trial_list','trial_condition')
    
    
    Fs=2000;
    nsub = group_level_stat_input.subnum;
    
    load(['R00_Output_RippleDetection_' sub],'ripple_encoding_ch','ripple_retrieval_ch','chan_ROI','MTL_roi','trials_time_axis_encoding','trials_time_axis_retrieval')
    %        good_ripple_chan = (sum(output(output(:,1)==nsub,2:3),2) >= 150 ) & output(output(:,1)==nsub,8)==1; good_ripple_chan = good_ripple_chan';
    [~,good_ripple_chan1] = max(sum(output(output(:,1)==nsub,2:3),2))  ; good_ripple_chan =zeros(1,sum(output(:,1)==nsub))'; good_ripple_chan(good_ripple_chan1)=1;
    good_ripple_chan = good_ripple_chan & (output(output(:,1)==nsub,8)==1) & sum(output(output(:,1)==nsub,2:3),2) >= 150 ;
    if 0~=sum(good_ripple_chan)
        ripple_ev_rtv = zeros(length(Ecd_trial_list),2*Fs);
        ripple_ev_ecd = zeros(length(Ecd_trial_list),4*Fs);
        ripple_ev_rtv_locked2rt = zeros(length(Ecd_trial_list),[4+2]*Fs);    % The 4*Fs+1 samples corresponds to the response time
        ripple_ev_ecd_locked2rt = zeros(length(Ecd_trial_list),[4+2]*Fs);
        %% extract retrieval
        ripple_trial=[];ripple_trial_time=[];
        for nch = find(good_ripple_chan==1),
            ripple_trial = cat(2,ripple_trial, ripple_retrieval_ch(nch).trial);
            ripple_trial_time = cat(2,ripple_trial_time,ripple_retrieval_ch(nch).trial_time);
        end
        rt_rtv = nan(1,147);
          for rpl_ind = 1:length(ripple_trial)
            if (ripple_trial_time(rpl_ind)>0 & ismember(ripple_trial(rpl_ind),Rtv_trial_list) & ~isnan(ripple_trial_time(rpl_ind)))
                trialind = find(ripple_trial(rpl_ind)==Rtv_trial_list);
                rt_rtv(trialind) = EMOPSData.testVar(ripple_trial(rpl_ind),8); if rt_rtv(trialind)<0.01, rt_rtv(trialind)=4;end;
                if rt_rtv(trialind)<2, rt_rtv(trialind)=2;end
                iti_duration = EMOPSData.testVar(ripple_trial(rpl_ind),9);
                if (iti_duration + rt_rtv(trialind) ) >= (ripple_trial_time(rpl_ind)/Fs)
                    if (round((4-rt_rtv(trialind))*Fs)+ripple_trial_time(rpl_ind)) > 12000
                        disp('opppps')
                    end
                    ripple_ev_rtv_locked2rt(trialind,round((4-rt_rtv(trialind))*Fs)+ripple_trial_time(rpl_ind)  ) = 1;
                end
            end
        end
        bad_rpl_event = find(1==isnan(ripple_trial_time) | ripple_trial_time<0 | ripple_trial_time>2*Fs | ...
            ~ismember(ripple_trial,Rtv_trial_list));
        ripple_trial(bad_rpl_event) = []; ripple_trial_time(bad_rpl_event) = [];
        for rpl_ind = 1:length(ripple_trial)
            trialind = find(ripple_trial(rpl_ind)==Rtv_trial_list);
            ripple_ev_rtv(trialind,ripple_trial_time(rpl_ind) ) = 1;
        end
        
        %% extract encoding
        ripple_trial=[];ripple_trial_time=[];
        for nch = find(good_ripple_chan==1),
            ripple_trial = cat(2,ripple_trial, ripple_encoding_ch(nch).trial);
            ripple_trial_time = cat(2,ripple_trial_time,ripple_encoding_ch(nch).trial_time);
        end
        rt_ecd = nan(1,147);
        for rpl_ind = 1:length(ripple_trial)
            if (ripple_trial_time(rpl_ind)>0 & ismember(ripple_trial(rpl_ind),Ecd_trial_list) & ~isnan(ripple_trial_time(rpl_ind)))
                trialind = find(ripple_trial(rpl_ind)==Ecd_trial_list);
                rt_ecd(trialind) = EMOPSData.studyVar(ripple_trial(rpl_ind),5);if rt_ecd(trialind)<0.01, rt_ecd(trialind)=4;end
                if rt_ecd(trialind)<2, rt_ecd(trialind)=2;end
                iti_duration = EMOPSData.studyVar(ripple_trial(rpl_ind),6);
                if (iti_duration + rt_ecd(trialind) ) >= (ripple_trial_time(rpl_ind)/Fs)
                    if (round((4-rt_ecd(trialind))*Fs)+ripple_trial_time(rpl_ind)) > 12000
                        disp('opppps')
                    end
                    ripple_ev_ecd_locked2rt(trialind,round((4-rt_ecd(trialind))*Fs)+ripple_trial_time(rpl_ind)  ) = 1;
                    ripple_ev_ecd_locked2rt(trialind,1:round((6-rt_ecd(trialind))*Fs)  ) = 0;
                end
                
            end
        end
        bad_rpl_event = find(1==isnan(ripple_trial_time) | ripple_trial_time<0 | ripple_trial_time>4*Fs | ...
            ~ismember(ripple_trial,Ecd_trial_list));
        ripple_trial(bad_rpl_event) = []; ripple_trial_time(bad_rpl_event) = [];
        for rpl_ind = 1:length(ripple_trial)
            trialind = find(ripple_trial(rpl_ind)==Ecd_trial_list);
            ripple_ev_ecd(trialind,ripple_trial_time(rpl_ind) ) = 1;
        end
        %% extract other propeties
            
        nanind = find(1==isnan(trial_condition)) ;
        arolevel = EMOPSData.testVar( Rtv_trial_list,14) ;
        thre_aro = 4.13; LFA_trial = find(trial_condition==1); RPT_trial = find(trial_condition==0);
        trial_condition_emo = nan(1,147); trial_condition_emo(arolevel>thre_aro)=2; trial_condition_emo(arolevel<=thre_aro)=1;
        trial_condition_emo(nanind) = nan;
        trial_condition_emo(RPT_trial) = 0;
        trial_condition_arolevel = arolevel;
        
        
        group_level_stat_input.ripple_raster_rtv{nsub} = ripple_ev_rtv_locked2rt;
        group_level_stat_input.ripple_raster_ecd{nsub} = ripple_ev_ecd_locked2rt;
        group_level_stat_input.trial_condition{nsub} = trial_condition ;
        group_level_stat_input.trial_condition_emo{nsub} = trial_condition_emo ;
        group_level_stat_input.trial_condition_arolevel{nsub} = trial_condition_arolevel;
        group_level_stat_input.ripple_raster_rtv0{nsub} = ripple_ev_rtv;
        group_level_stat_input.ripple_raster_ecd0{nsub} = ripple_ev_ecd;
        group_level_stat_input.rt_rtv{nsub}=EMOPSData.testVar(Rtv_trial_list,8);
        group_level_stat_input.rt_ecd{nsub}=EMOPSData.studyVar(Ecd_trial_list,5);
        %
    else
        ripple_summary(1:8,1:3) = nan;
        tmp_re_rtv(1:4,1:2*2000) = nan;
        tmp_re_ecd(1:4,1:4*2000) = nan;
        tmp_re_ecd_rt(1:4,1:6*2000) = nan;
        tmp_re_rtv_rt(1:4,1:6*2000) = nan;
        ripple_presence=[];
        ripple_rate=[];
    end
    
    
end
clearvars -except SC_group_level_stat_input ch ok_epochs_study ok_epochs_test sub chan_ROI tt output group_level_stat_input



sublist=[2:5 6 7];
trial_cond_outcome=[];trial_cond_aro=[];trial_cond_arovalue=[];tmp_ecd = [];sub = []; trial_RSA = [];trial_rt = [];
trial_RSA_next=[];
for nsub = sublist
    trial_cond_outcome = cat(2,trial_cond_outcome,group_level_stat_input.trial_condition{nsub});
    trial_cond_aro = cat(2,trial_cond_aro,group_level_stat_input.trial_condition_emo{nsub});
    trial_cond_arovalue = cat(2,trial_cond_arovalue,group_level_stat_input.trial_condition_arolevel{nsub});
    tmp_ecd = cat(1,tmp_ecd,group_level_stat_input.ripple_raster_ecd0{nsub}) ;
    trial_rt = cat(1,trial_rt,group_level_stat_input.rt_ecd{nsub}) ;

    trial_RSA = cat(3,trial_RSA,mean(SC_group_level_stat_input.Similarity_same_pair_ecdwresponse{nsub},4,'omitnan'));
    trial_RSA_next = cat(3,trial_RSA_next,mean(SC_group_level_stat_input.Similarity_next_pair_ecdwresponse{nsub},4,'omitnan'));
    sub = cat(1,sub,ones(147,1)*nsub);
end
% ind_include = find(trial_cond_outcome==1 | trial_cond_outcome==2);
ind_include = find(~isnan(trial_cond_outcome));


Fs=2000;
taxis_ripple = (2*Fs+1):4*Fs ;  % (2*Fs+1):4*Fs ;
tmp_ripple = tmp_ecd(ind_include,taxis_ripple) ;
trial_RSA = trial_RSA(:,:,ind_include) ;
trial_RSA_next = trial_RSA_next(:,:,ind_include) ;
sub = sub(ind_include) ;
trial_cond_outcome = trial_cond_outcome(ind_include) ;
trial_cond_aro = trial_cond_aro(ind_include);
trial_cond_arovalue = trial_cond_arovalue(ind_include);

% As Null distribution (Regardless of Ripple existance )
for nsub = 2:7
    Null_RSA_sub(:,:,nsub-1) = mean(trial_RSA(:,:,sub==nsub),3,'omitnan') ;
    Null_RSA_next_sub(:,:,nsub-1) = mean(trial_RSA_next(:,:,sub==nsub),3,'omitnan') ;
end


 trial_rt = trial_rt(ind_include);

tt = SC_group_level_stat_input.t_RSA_ecd_rt;
tt_ecd = SC_group_level_stat_input.t_RSA_ecd ; 
Ntp = length(tt_ecd);
for k = 1:size(trial_RSA ,3 ) 
    for n = 1:(Ntp+5)
        trial_RSA( max(n-5,1) : min(n+5,Ntp) , max(n-5,1) : min(n+5,Ntp+5) , k ) = nan; 
        trial_RSA_next( max(n-5,1) : min(n+5,Ntp) , max(n-5,1) : min(n+5,Ntp+5) , k ) = nan; 
    end
end


% RSA_rw = cat(2,RSA_rw(:,[1 1 1 1 1],:),RSA_rw,RSA_rw(:,[end end end end end],:)) ;
% tt = cat(2, tt(1)-(0.05:-0.01:0.01), tt , tt(end) + (0.01:0.01:0.05) ) ;

nrr = 0 ;
rw_rsa = []; rw_sub = []; rw_tt = []; rw_cond = []; rw_cond_aro = [];rw_cond_arovalue = [];
rw_rsa_DifferentTrial = [];
T = 0.25 ; Fs_rsa = 100;
for n1 = 1:size(tmp_ripple,1)-1
    rpp_ts = find(tmp_ripple(n1,:)>=1) ;
%    rpp_ts = rpp_ts(rpp_ts>0*Fs & rpp_ts<=2*Fs);
    rpp_ts = rpp_ts(rpp_ts>0*Fs & rpp_ts<=(trial_rt(n1)-2+1.25)*Fs);
    if ~isempty(rpp_ts)
        for ntp = rpp_ts
            nrr = nrr+1 ;
            ntt =  2 + (ntp/Fs);       % I changed here  ntt =  (ntp/Fs);   
%             ind = find(tt>=(ntt-T) & tt<=(ntt+T));
%             ind = ind(1:49);
            ind = round(((ntt-T)*Fs_rsa+1):(ntt+T)*Fs_rsa);
            ind = ind(1:49);
            nan_ind = find(ind <=0 | ind>=length(tt)) ;
            ind(nan_ind) = 1;
            rw_rsa(:,:,nrr) = trial_RSA(:,ind,n1 ) ;
            rw_rsa(:,nan_ind,nrr) = nan;
            rw_rsa_DifferentTrial(:,:,nrr) = trial_RSA_next(:,ind,n1) ;
            rw_rsa_DifferentTrial(:,nan_ind,nrr) = nan;
            
            rw_sub(nrr) = sub(n1); rw_tt(nrr) = ntt;
            rw_cond(nrr) = trial_cond_outcome(n1);
            rw_cond_aro(nrr) = trial_cond_aro(n1);
             rw_cond_arovalue(nrr) = trial_cond_arovalue(n1);
        end
    end
end

figure, imagesc(median( rw_rsa -  rw_rsa_DifferentTrial   ,3 , 'omitnan'))
colormap jet
colorbar

%% Grand average on current trial - different trial difference 
for nsub = 1:length(sublist)
    rw_rsa_difference(:,:,nsub) = mean( rw_rsa(:,:,rw_sub==sublist(nsub)) -  rw_rsa_DifferentTrial(:,:,rw_sub==sublist(nsub)) ,3,'omitnan'  ); 
    figure(50), subplot(2,3,nsub),
    imagesc(rw_rsa_difference(:,:,nsub))
    colormap jet
colorbar
end

%% Stats on grand avg level 
Null_rsa_difference = Null_RSA_sub - Null_RSA_next_sub;
Null_rsa_responseAvg = Null_rsa_difference(:,226:(225+49),:);
Null_rsa_responseAvg(:,1:49,:) =repmat( mean( Null_rsa_difference(:,400:450,:) ,2,'omitnan' ) ,[1 49 1]);
% Null_rsa_responseAvg(:,1:49,:) =repmat( mean( Null_rsa_difference(:,225:(400-25),:) ,2,'omitnan' ) ,[1 49 1]);
[h,p,~,stats] = ttest(rw_rsa_difference,Null_rsa_responseAvg,'Dim', 3,'tail','right');
figure, imagesc(linspace(-0.25,0.25,49) , linspace(0,2,200),p), caxis([0 0.1])

[labeledA, numRegions] = bwlabel(h , 8);
true_t_val = stats.tstat;
for n_cluster = 1:numRegions , true_t_val(labeledA==n_cluster)=sum(true_t_val(labeledA==n_cluster));end
figure, imagesc(true_t_val)

Nperm = 1000 ;
for nperm = 1:Nperm
    rw_rsa_difference_shuffle = [];
    perm_ind = randi(2,[size(rw_rsa,3) 1]) -1;
    perm_ind=zeros(size(perm_ind));
    perm_ind(randperm(size(rw_rsa,3) , round(0.05*size(rw_rsa,3)))) = 1; 
    rw_rsa_shuffle = rw_rsa;
    rw_rsa_shuffle(:,:,perm_ind==1) = rw_rsa_DifferentTrial(:,:,perm_ind==1);
    rw_rsa_DifferentTrial_shuffle = rw_rsa_DifferentTrial;
    rw_rsa_DifferentTrial_shuffle(:,:,perm_ind==1)  = rw_rsa(:,:,perm_ind==1) ;
    
    for nsub = 1:length(sublist)
        rw_rsa_difference_shuffle(:,:,nsub) = mean( rw_rsa_shuffle(:,:,rw_sub==sublist(nsub)) -  rw_rsa_DifferentTrial_shuffle(:,:,rw_sub==sublist(nsub)) ,3,'omitnan'  );
    end
    [h,p,~,stats] = ttest(rw_rsa_difference_shuffle,Null_rsa_responseAvg,'Dim', 3,'tail','right');  
    
    [labeledA, numRegions] = bwlabel(h , 8);
    shuffle_t_val = stats.tstat;
    for n_cluster = 1:numRegions , shuffle_t_val(labeledA==n_cluster)=sum(shuffle_t_val(labeledA==n_cluster));end
    t_cluster_max(nperm) = max(shuffle_t_val,[],'all');
    
end

figure, imagesc(linspace(-0.25,0.25,49) , linspace(0,2,200) , mean(   true_t_val<reshape(t_cluster_max,1, 1 ,Nperm)  ,3) )

caxis([0 0.1])

p_cluster =  mean(   true_t_val<reshape(t_cluster_max,1, 1 ,Nperm)  ,3)  ; 



% Plot the figure 
 figure, imagesc(linspace(-0.25,0.25,49) , linspace(0,2,200),mean(rw_rsa_difference,3))
colormap jet
colorbar
 dim1 = linspace(-0.25,0.25,49) ; dim2 =  linspace(0,2,200); 
 BW = p_cluster<0.05; BW = imfilter(BW,[1 1;1 1]); 
 [B,~] = bwboundaries(BW, 'noholes') ; boundary = B{1} ; 
 hold on, plot(dim1(boundary(:,2)) , dim2(boundary(:,1)) , 'k', 'LineWidth',2) 






