load('dualstate_freqPowerPart.mat')

n_theta = 5:6; 
Fs=2000;
ripple_tmp_rate_ecd = squeeze(group_level_stat_input.tmp_re_ecd(1,:,:)) ;
ripple_tmp_rate_ecd=ripple_tmp_rate_ecd(1:2*Fs,:);
avg_theta_power = squeeze(mean(TFR_theta(n_theta,:,1:7),[1],'omitnan') );
avg_theta_power = avg_theta_power(Fs+1:(3*Fs),:) ;
% figure
state_theta_SWR=[];state_theta_SWR_grand=[];
output=[]; output_sub=[];
z_SWR_grand=[];z_nonSWR_grand=[];
thre1 = 0 ;
thre2 = -0 ;
for nsub=2:7
    rippletrial = sum( group_level_stat_input.ripple_raster_ecd0{nsub}(:, 1:3.5*Fs) ,2,'omitnan');
    rippletrial = (rippletrial>0.1);

    ripple_rate_ripple_trial = group_level_stat_input.ripple_raster_ecd0{nsub}(find(rippletrial==1), 1:2.5*Fs);
    
    theta_band_encoding=[];
    theta_band_encoding_rt=[];
    theta_band_encoding = SC_group_level_stat_input_SpectralPower.theta_band_encoding{nsub};
    
    theta_band_encoding = (theta_band_encoding - mean(theta_band_encoding(:,(Fs/2):6*Fs,:,:),2,'omitnan'))./std(theta_band_encoding(:,(Fs/2):6*Fs,:,:),[],2,'omitnan') ;
    theta_band_encoding = theta_band_encoding - mean(theta_band_encoding(:,(Fs/2):Fs,:,:),2,'omitnan') ;
    theta_band_encoding = mean(theta_band_encoding,4,'omitnan');
    
    

    theta_power_ripple_trial = squeeze(mean(theta_band_encoding(n_theta,(Fs+1):3.5*Fs,:),1,'omitnan'));
    theta_power_ripple_trial = (theta_power_ripple_trial(:,find(rippletrial>=1)))';
    

    
    theta_state_trial1 = theta_power_ripple_trial>=thre1 ;
    theta_state_trial2 = -(theta_power_ripple_trial<=thre2) ;
    theta_state_trial = theta_state_trial1 + theta_state_trial2 ; 
    d = theta_state_trial(find(ripple_rate_ripple_trial==1)) ;

    
    
%     state_theta_SWR = cat(1,state_theta_SWR,d);
    state_theta_SWR_grand(nsub,:) = [sum(d==-1) sum(d==1)];
%     output = cat(1,output,[reshape(theta_state_trial,[],1) reshape(ripple_rate_ripple_trial,[],1)]); 
%     output_sub = cat(1,output_sub,ones(length(reshape(theta_state_trial,[],1)),1) * nsub) ; 
    %     [r_coef(nsub),r_pval(nsub)] = corr(theta_power_ripple_trial' ,ripple_rate_ripple_trial' ,'type','Spearman')  ;
    %     subplot(2,4,nsub),plot(theta_power_ripple_trial),yyaxis right,  plot(ripple_rate_ripple_trial)
    %     title(num2str(r_coef(nsub)))

end

[p,h,stat]=signrank(state_theta_SWR_grand(2:7,1) , state_theta_SWR_grand(2:7,2),'method','approximate','tail','right')
fVec = 3:15; 

figure, bar([median(state_theta_SWR_grand(2:7,1)) median(state_theta_SWR_grand(2:7,2))],'FaceColor',[0.7 0.7 0.7])
hold on, errorbar([median(state_theta_SWR_grand(2:7,1)) median(state_theta_SWR_grand(2:7,2))] , [mysem(state_theta_SWR_grand(2:7,1)) mysem(state_theta_SWR_grand(2:7,2)) ],'LineStyle','none','Color','k','LineWidth',2)
hold on, plot(state_theta_SWR_grand(2:7,:)','LineWidth',2,'Marker','o') , xlim([0 3]) 
set(gca,'LineWidth',2)
title(num2str(p))

z_state_theta_SWR_grand = state_theta_SWR_grand./(sum(state_theta_SWR_grand,2)); 
figure, bar([median(z_state_theta_SWR_grand(2:7,1)) median(z_state_theta_SWR_grand(2:7,2))],'FaceColor',[0.7 0.7 0.7])
hold on, errorbar([median(z_state_theta_SWR_grand(2:7,1)) median(z_state_theta_SWR_grand(2:7,2))] , [mysem(z_state_theta_SWR_grand(2:7,1)) mysem(z_state_theta_SWR_grand(2:7,2)) ],'LineStyle','none','Color','k','LineWidth',2)
hold on, plot(z_state_theta_SWR_grand(2:7,:)','LineWidth',2,'Marker','o') , xlim([0 3]) 
set(gca,'LineWidth',2)
title(num2str(p))
ylim([0.3 0.7])

[p,h,stat]=signrank(z_state_theta_SWR_grand(2:7,1) , z_state_theta_SWR_grand(2:7,2),'method','approximate','tail','right')



%% Gamma state and aSWR state 
state_theta_SWR_grand=[];state_gamma_SWR_grand=[];
thre1=0; thre2=0;
for nsub=2:7
    rippletrial = sum( group_level_stat_input.ripple_raster_ecd0{nsub}(:, 1:3.5*Fs) ,2,'omitnan');
    rippletrial = (rippletrial>0.1);
    %rippletrial = (rippletrial'>=0.1) & ismember(group_level_stat_input.trial_condition{nsub},[1 2 ]);
    ripple_rate_ripple_trial = group_level_stat_input.ripple_raster_ecd0{nsub}(find(rippletrial==1), 1:2.5*Fs);
    
    theta_band_encoding=[];
    theta_band_encoding_rt=[];
    theta_band_encoding = SC_group_level_stat_input_SpectralPower.theta_band_encoding{nsub};
    
    theta_band_encoding = (theta_band_encoding - mean(theta_band_encoding(:,(Fs/2):6*Fs,:,:),2,'omitnan'))./std(theta_band_encoding(:,(Fs/2):6*Fs,:,:),[],2,'omitnan') ;
    theta_band_encoding = theta_band_encoding - mean(theta_band_encoding(:,(Fs/2):Fs,:,:),2,'omitnan') ;
    theta_band_encoding = mean(theta_band_encoding,4,'omitnan');
    
    
    
    gamma_band_encoding=[];
    gamma_band_encoding = SC_group_level_stat_input_SpectralPower.gamma_band_encoding{nsub};
    gamma_band_encoding = (gamma_band_encoding - mean(gamma_band_encoding(:,(Fs/2):6*Fs,:,:),2,'omitnan'))./std(gamma_band_encoding(:,(Fs/2):6*Fs,:,:),[],2,'omitnan') ;
    gamma_band_encoding = gamma_band_encoding - mean(gamma_band_encoding(:,(Fs/2):Fs,:,:),2,'omitnan') ;
    gamma_band_encoding = mean(gamma_band_encoding,4,'omitnan');
    
    
    
%     xx= mean(theta_band_encoding(:,(Fs+1):3*Fs,:),[2 3],'omitnan');
%      [pk,loc] = findpeaks(xx);
%      n_theta = (max(1,loc-4):min(13,loc+4) ); 

    theta_power_ripple_trial = squeeze(mean(theta_band_encoding(n_theta,(Fs+1):3.5*Fs,:),1,'omitnan'));
    theta_power_ripple_trial = (theta_power_ripple_trial(:,find(rippletrial>=1)))';
    gamma_power_ripple_trial = squeeze(mean(gamma_band_encoding(5:end,(Fs+1):3.5*Fs,:),1,'omitnan'));
    gamma_power_ripple_trial = (gamma_power_ripple_trial(:,find(rippletrial>=1)))';
    
   
    
    
%     thre1 = quantile(theta_power_ripple_trial,0.75,'all') ;
%     thre2 = quantile(theta_power_ripple_trial,0.25,'all') ;
    
    theta_state_trial1 = theta_power_ripple_trial>=thre1 ;
    theta_state_trial2 = -(theta_power_ripple_trial<=thre2) ;
    theta_state_trial = theta_state_trial1 + theta_state_trial2 ; 
    d = theta_state_trial(find(ripple_rate_ripple_trial==1)) ;
    state_theta_SWR_grand(nsub,:) = [sum(d==-1) sum(d==1)];
    
    %thre1 = quantile(theta_power_ripple_trial,0.5,'all'); thre2=thre1; 
    gamma_state_trial1 = gamma_power_ripple_trial>=thre1 ;
    gamma_state_trial2 = -(gamma_power_ripple_trial<=thre2) ;
    gamma_state_trial = gamma_state_trial1 + gamma_state_trial2 ; 
    d = gamma_state_trial(find(ripple_rate_ripple_trial==1)) ;
    state_gamma_SWR_grand(nsub,:) = [sum(d==-1) sum(d==1)];
    
    
    

end

[p,h,stat]=signrank(state_gamma_SWR_grand(2:7,1) , state_gamma_SWR_grand(2:7,2),'method','approximate','tail','right')

z_state_gamma_SWR_grand = state_gamma_SWR_grand./(sum(state_gamma_SWR_grand,2)); 
figure, bar([median(z_state_gamma_SWR_grand(2:7,1)) median(z_state_gamma_SWR_grand(2:7,2))],'FaceColor',[0.7 0.7 0.7])
hold on, errorbar([median(z_state_gamma_SWR_grand(2:7,1)) median(z_state_gamma_SWR_grand(2:7,2))] , [mysem(z_state_gamma_SWR_grand(2:7,1)) mysem(z_state_gamma_SWR_grand(2:7,2)) ],'LineStyle','none','Color','k','LineWidth',2)
hold on, plot(z_state_gamma_SWR_grand(2:7,:)','LineWidth',2,'Marker','o') , xlim([0 3]) 
set(gca,'LineWidth',2)
title(num2str(p))
ylim([0.3 0.7])


figure, bar([median(state_gamma_SWR_grand(2:7,1)) median(state_gamma_SWR_grand(2:7,2))],'FaceColor',[0.7 0.7 0.7])
hold on, errorbar([median(state_gamma_SWR_grand(2:7,1)) median(state_gamma_SWR_grand(2:7,2))] , [mysem(state_gamma_SWR_grand(2:7,1)) mysem(state_gamma_SWR_grand(2:7,2)) ],'LineStyle','none','Color','k','LineWidth',2)
hold on, plot(state_gamma_SWR_grand(2:7,:)','LineWidth',2,'Marker','o') , xlim([0 3]) 
set(gca,'LineWidth',2)
title(num2str(p))
