% plot spontaneous rate in two different states: awake and anesthetized
% plot firing rate in response to maskers in two different states: awake and anesthetized
load('sp_context.mat')
dur_echo=1.3311;
dur_distress=1.9266;

c_=sp_context(:,2)==0;   % unit classification
aw_=sp_context(:,3)==0;  % awake data 
an_=sp_context(:,3)==1;  % anesthetized data 
spont_=sp_context(:,7)==0;    % spontaneous firing rate
cont_=sp_context(:,7)==1;     % context firing rate intact
cont_chi=sp_context(:,7)==2;  % context firing rate chimeras
cont_pt=sp_context(:,7)==3;   % context firing rate pure tones 

% submatrices 
aw_spont=sp_context(logical(c_.*aw_.*spont_),:); 
an_spont=sp_context(logical(c_.*an_.*spont_),:);

aw_cont=sp_context(logical(c_.*aw_.*cont_),:);
an_cont=sp_context(logical(c_.*an_.*cont_),:);

an_chi=sp_context(logical(c_.*an_.*cont_chi),:);


%% calculate spontaneous firing rate

% awake data - spontaneous & context 
units=unique(aw_spont(:,1));
aw_spontR=nan(numel(units),1);
aw_contR=nan(numel(units),2);
for u=1:numel(units)
    u_=aw_spont(:,1)==units(u);
    u_2=aw_cont(:,1)==units(u);
    rate=nanmean(aw_spont(u_,8));
    rate2=nanmean(aw_cont(u_2,5));   % nav
    rate3=nanmean(aw_cont(u_2,6));   % com 
    aw_spontR(u,1)=rate;  
    aw_contR(u,1)=rate2;  
    aw_contR(u,2)=rate3;  
end
aw_spontR=aw_spontR./((3.5-0.2));
aw_contR(:,1)=aw_contR(:,1);
aw_contR(:,2)=aw_contR(:,2);

% anesthetized data - spontaneous 
units=unique(an_spont(:,1));
an_spontR=nan(numel(units),1);
for u=1:numel(units)
    u_=an_spont(:,1)==units(u);
    rate=nanmean(an_spont(u_,8));
    an_spontR(u,1)=rate;   
end
an_spontR=an_spontR./((3.5-0.2));

% anesthetized data - context 
units=unique(an_cont(:,1));
an_contR=nan(numel(units),2);
for u=1:numel(units)
    u_2=an_cont(:,1)==units(u);
    rate2=nanmean(an_cont(u_2,5));   % nav
    rate3=nanmean(an_cont(u_2,6));   % com 
    an_contR(u,1)=rate2;  
    an_contR(u,2)=rate3;  
end
an_contR(:,1)=an_contR(:,1);
an_contR(:,2)=an_contR(:,2);

% anesthetized data - context CHIMERAS 
units=unique(an_chi(:,1));
an_chiR=nan(numel(units),2);
for u=1:numel(units)
    u_2=an_chi(:,1)==units(u);
    rate2=nanmean(an_chi(u_2,5));   % nav
    rate3=nanmean(an_chi(u_2,6));   % com 
    an_chiR(u,1)=rate2;  
    an_chiR(u,2)=rate3;  
end
an_chiR(:,1)=an_chiR(:,1); % flipped 
an_chiR(:,2)=an_chiR(:,2);

% firing rate in response to probes
int=load('all_data.mat'); int=int.all_di; 
c_=int(:,2)==0;
an_=int(:,3)==1;
int=int(logical(c_.*an_),18:19);
chi=load('chi_data.mat'); chi=chi.all_di; 
c_=chi(:,2)==0;
an_=chi(:,3)==0;
chi=chi(logical(c_.*an_),18:19);
pt=load('pt_data.mat'); pt=pt.all_di; 
c_=pt(:,2)==0;
an_=pt(:,3)==1;
pt=pt(logical(c_.*an_),11:12);
probe_an=[int; chi].*20 ; %; chi; pt].*20;             % x20: spikecounts to firing rate
aw=load('all_data.mat'); aw=aw.all_di; 
c_=aw(:,2)==0;
aw_=aw(:,3)==0;
probe_aw=aw(logical(c_.*aw_),18:19).*20; % x20: spikecounts to firing rate 


% compare non-paired samples 
pv1=ranksum(aw_spontR,an_spontR);
pv2=ranksum(aw_contR(:,1),an_contR(:,1));
pv3=ranksum(aw_contR(:,2),an_contR(:,2));
pv4=signrank(aw_contR(:,1),aw_contR(:,2));
pv5=signrank(an_contR(:,1),an_contR(:,2));
pv6=ranksum(an_contR(:,1),an_chiR(:,1));
pv7=ranksum(an_contR(:,2),an_chiR(:,2));
pv8=ranksum(probe_aw(:,1),probe_an(:,1));
pv9=ranksum(probe_aw(:,2),probe_an(:,2));
pv10=signrank(probe_aw(:,1),probe_aw(:,2));
pv11=signrank(probe_an(:,1),probe_an(:,2));

%% plotting
close all;
red =[0.8500,0.3250, 0.0980];
blue=[0,0.4470,0.7410];
red_light=[1,0.5,0.5].*red;
blue_light=[1,0.5,0.5].*blue;
colors=[[0.3 0.3 0.3];[0.5 0.5 0.5];blue_light;red_light;blue;red];

y1=nan(58,6);
y1(1:numel(aw_spontR),1)=aw_spontR; 
y1(1:numel(an_spontR),2)=an_spontR;
y1(1:numel(aw_contR(:,1)),3)=aw_contR(:,1);
y1(1:numel(an_contR(:,1)),4)=an_contR(:,1);
y1(1:numel(aw_contR(:,2)),5)=aw_contR(:,2); 
y1(1:numel(an_contR(:,2)),6)=an_contR(:,2);

x = 1:6;

figure(1); set(gcf,'Position',[400 200 250 170])
ax = axes();
hold(ax);
for i=1:6
    h=boxchart(x(i)*ones(size(y1(:,i))), y1(:,i), 'BoxFaceColor', colors(i,:),'Notch','on');
    h.MarkerStyle='.';
    h.MarkerColor=colors(i,:);
end
hold on
plot(1:2,nanmedian(y1(:,1:2)),'o','LineWidth',1,'Color',[0.3 0.3 0.3],'MarkerSize',3,'MarkerFaceColor',[0.3 0.3 0.3])
plot(3,nanmedian(y1(:,3)),'o','LineWidth',1,'Color',blue_light,'MarkerSize',3,'MarkerFaceColor',blue_light)
plot(4,nanmedian(y1(:,4)),'o','LineWidth',1,'Color',red_light,'MarkerSize',3,'MarkerFaceColor',red_light)
plot(5,nanmedian(y1(:,5)),'o','LineWidth',1,'Color',blue,'MarkerSize',3,'MarkerFaceColor',blue)
plot(6,nanmedian(y1(:,6)),'o','LineWidth',1,'Color',red,'MarkerSize',3,'MarkerFaceColor',red)
hold off
xlim([0.5 6.5])
ylabel('# spikes')
xticks(1:6)
xticklabels({'awake','anesthetized','awake','anesthetized','awake','anesthetized'});

set(gca, 'box', 'off')
set(gca, 'Color','none')
set(gca,'linewidth',1);set(gca,'fontsize',8);

xt = get(gca, 'XTick');
line1=165;
line2=190;
line3=215;

sl=0.75; 
% plot lines for comparison 
hold on 
plot(xt(1:2), [1 1]*line1, '-k','LineWidth',sl) 
plot([1;1]*xt(1),[line1,line1-3],'-k', 'LineWidth',sl);
plot([1;1]*xt(2),[line1,line1-3],'-k', 'LineWidth',sl);

plot(xt(3:4), [1 1]*line1, '-k','LineWidth',sl) 
plot([1;1]*xt(3),[line1,line1-3],'-k', 'LineWidth',sl);
plot([1;1]*xt(4),[line1,line1-3],'-k', 'LineWidth',sl);

plot(xt(5:6), [1 1]*line1, '-k','LineWidth',sl) 
plot([1;1]*xt(5),[line1,line1-3],'-k', 'LineWidth',sl);
plot([1;1]*xt(6),[line1,line1-3],'-k', 'LineWidth',sl);

plot([xt(3) xt(5)], [1 1]*line2, '-k','LineWidth',sl) 
plot([1;1]*xt(3),[line2,line2-3],'-k', 'LineWidth',sl);
plot([1;1]*xt(5),[line2,line2-3],'-k', 'LineWidth',sl);

plot([xt(4) xt(6)], [1 1]*line3, '-k','LineWidth',sl) 
plot([1;1]*xt(4),[line3,line3-3],'-k', 'LineWidth',sl);
plot([1;1]*xt(6),[line3,line3-3],'-k', 'LineWidth',sl);

ylim([-5 170])

% plot level of significance
if pv1 <= 0.05 && pv1 > 0.01
    text(mean(xt(1:2))-0.05, line1+5, '*','HorizontalAlignment','center')
elseif pv1 <=0.01 && pv1 > 0.001
    text(mean(xt(1:2))-0.05, line1+5, '**','HorizontalAlignment','center')
elseif pv1 <= 0.001
    text(mean(xt(1:2))-0.05, line1+5, '***','HorizontalAlignment','center')
else 
    text(mean(xt(1:2))-0.05, line1+12, 'ns','FontSize',8,'HorizontalAlignment','center')
end

if pv2 <= 0.05 && pv2 > 0.01
    text(mean(xt(3:4))-0.05, line1+5, '*','HorizontalAlignment','center')
elseif pv2 <=0.01 && pv2 > 0.001
    text(mean(xt(3:4))-0.05, line1+5, '**','HorizontalAlignment','center')
elseif pv2 <= 0.001
    text(mean(xt(3:4))-0.05, line1+5, '***','HorizontalAlignment','center')
else 
    text(mean(xt(3:4))-0.05, line1+12, 'ns','FontSize',8,'HorizontalAlignment','center')
end


if pv3 <= 0.05 && pv3 > 0.01
    text(mean(xt(5:6))-0.05, line1+5, '*','HorizontalAlignment','center')
elseif pv3 <=0.01 && pv3 > 0.001
    text(mean(xt(5:6))-0.05, line1+5, '**','HorizontalAlignment','center')
elseif pv3 <= 0.001
    text(mean(xt(5:6))-0.05, line1+5, '***','HorizontalAlignment','center')
else 
    text(mean(xt(5:6))-0.05, line1+15, 'ns','FontSize',8,'HorizontalAlignment','center')
end


% if pv4 <= 0.05 && pv4 > 0.01
%     text(mean(xt(3:5))-0.05, line2+5, '*','HorizontalAlignment','center')
% elseif pv4 <=0.01 && pv4 > 0.001
%     text(mean(xt(3:5))-0.05, line2+5, '**','HorizontalAlignment','center')
% elseif pv4 <= 0.001
%     text(mean(xt(3:5))-0.05, line2+5, '***','HorizontalAlignment','center')
% else 
%     text(mean(xt(3:5))-0.05, line2+12, 'ns','FontSize',8,'HorizontalAlignment','center')
% end
% 
% if pv5 <= 0.05 && pv5 > 0.01
%     text(mean(xt(4:6))-0.05, line3+5, '*','HorizontalAlignment','center')
% elseif pv5 <=0.01 && pv5 > 0.001
%     text(mean(xt(4:6))-0.05, line3+5, '**','HorizontalAlignment','center')
% elseif pv5 <= 0.001
%     text(mean(xt(4:6))-0.05, line3+5, '***','HorizontalAlignment','center')
% else 
%     text(mean(xt(4:6))-0.05, line3+12, 'ns','FontSize',8,'HorizontalAlignment','center')
% end

% exportgraphics(gcf,'E:\Users\User\Desktop\delay paper\anesthetized\figures_paper\preparation_rate_cont.pdf',...
%     'Resolution',300','ContentType','vector','BackgroundColor','none')

%% plotting probes 
colors=[blue_light;blue;red_light;red];

y2=nan(58,4);
y2(1:numel(probe_aw(:,1)),1)=probe_aw(:,1); 
y2(1:numel(probe_an(:,1)),2)=probe_an(:,1);
y2(1:numel(probe_aw(:,2)),3)=probe_aw(:,2); 
y2(1:numel(probe_an(:,2)),4)=probe_an(:,2);

x = 1:4;

figure(2); set(gcf,'Position',[700 200 170 170])
ax = axes();
hold(ax);
for i=1:4
    h=boxchart(x(i)*ones(size(y2(:,i))), y2(:,i), 'BoxFaceColor', colors(i,:),'Notch','on');
    h.MarkerStyle='.';
    h.MarkerColor=colors(i,:);
end
hold on
plot(1:2,nanmedian(y2(:,1:2)),'-o','LineWidth',1,'Color',blue_light,'MarkerSize',3,'MarkerFaceColor',blue_light)
plot(3:4,nanmedian(y2(:,3:4)),'-o','LineWidth',1,'Color',red_light,'MarkerSize',3,'MarkerFaceColor',red_light)
hold off
xlim([0.5 4.5])
ylabel('firing rate ')
xticks(1:4)
xticklabels({'awake','anesthetized','awake','anesthetized'});

set(gca, 'box', 'off')
set(gca, 'Color','none')
set(gca,'linewidth',1);set(gca,'fontsize',8);

xt = get(gca, 'XTick');
line1=290;
line2=340;
line3=390;

sl=0.75; 
% plot lines for comparison 
hold on 
plot(xt(1:2), [1 1]*line1, '-k','LineWidth',sl) 
plot([1;1]*xt(1),[line1,line1-10],'-k', 'LineWidth',sl);
plot([1;1]*xt(2),[line1,line1-10],'-k', 'LineWidth',sl);

plot(xt(3:4), [1 1]*line1, '-k','LineWidth',sl) 
plot([1;1]*xt(3),[line1,line1-10],'-k', 'LineWidth',sl);
plot([1;1]*xt(4),[line1,line1-10],'-k', 'LineWidth',sl);

% plot([xt(1) xt(3)], [1 1]*line2, '-k','LineWidth',sl) 
% plot([1;1]*xt(1),[line2,line2-10],'-k', 'LineWidth',sl);
% plot([1;1]*xt(3),[line2,line2-10],'-k', 'LineWidth',sl);
% 
% plot([xt(2) xt(4)], [1 1]*line3, '-k','LineWidth',sl) 
% plot([1;1]*xt(2),[line3,line3-10],'-k', 'LineWidth',sl);
% plot([1;1]*xt(4),[line3,line3-10],'-k', 'LineWidth',sl);

ylim([-5 320])

% plot level of significance
if pv8 <= 0.05 && pv8 > 0.01
    text(mean(xt(1:2))-0.05, line1+20, '*','HorizontalAlignment','center')
elseif pv8 <=0.01 && pv8 > 0.001
    text(mean(xt(1:2))-0.05, line1+20, '**','HorizontalAlignment','center')
elseif pv8 <= 0.001
    text(mean(xt(1:2))-0.05, line1+20, '***','HorizontalAlignment','center')
else 
    text(mean(xt(1:2))-0.05, line1+28, 'ns','FontSize',8,'HorizontalAlignment','center')
end

if pv9 <= 0.05 && pv9 > 0.01
    text(mean(xt(3:4))-0.05, line1+20, '*','HorizontalAlignment','center')
elseif pv9 <=0.01 && pv9 > 0.001
    text(mean(xt(3:4))-0.05, line1+20, '**','HorizontalAlignment','center')
elseif pv9 <= 0.001
    text(mean(xt(3:4))-0.05, line1+20, '***','HorizontalAlignment','center')
else 
    text(mean(xt(3:4))-0.05, line1+28, 'ns','FontSize',8,'HorizontalAlignment','center')
end


% if pv10 <= 0.05 && pv10 > 0.01
%     text(mean(xt(1:3))-0.05, line2+20, '*','HorizontalAlignment','center')
% elseif pv10 <=0.01 && pv10 > 0.001
%     text(mean(xt(1:3))-0.05, line2+20, '**','HorizontalAlignment','center')
% elseif pv10 <= 0.001
%     text(mean(xt(1:3))-0.05, line2+20, '***','HorizontalAlignment','center')
% else 
%     text(mean(xt(1:3))-0.05, line2+28, 'ns','FontSize',8,'HorizontalAlignment','center')
% end
% 
% if pv11 <= 0.05 && pv11 > 0.01
%     text(mean(xt(2:4))-0.05, line3+20, '*','HorizontalAlignment','center')
% elseif pv11 <=0.01 && pv11 > 0.001
%     text(mean(xt(2:4))-0.05, line3+20, '**','HorizontalAlignment','center')
% elseif pv11 <= 0.001
%     text(mean(xt(2:4))-0.05, line3+20, '***','HorizontalAlignment','center')
% else 
%     text(mean(xt(2:4))-0.05, line3+28, 'ns','FontSize',8,'HorizontalAlignment','center')
% end

% exportgraphics(gcf,'E:\Users\User\Desktop\delay paper\anesthetized\figures_paper\preparation_rate_probe.pdf',...
%     'Resolution',300','ContentType','vector','BackgroundColor','none')

%% plotting chimeras 
green_light=[0.4863,0.3059,0.7098];   % green: [0.3922,0.8314,0.0745];
green=[0.2980,0.1294,0.5098];   % gree: [0.4667,0.6745,0.1882];
colors=[blue_light;green;red_light;green_light];

y3=nan(46,4);
y3(1:numel(an_contR(:,1)),1)=an_contR(:,1); 
y3(1:numel(an_chiR(:,1)),2)=an_chiR(:,1);
y3(1:numel(an_contR(:,2)),3)=an_contR(:,2); 
y3(1:numel(an_chiR(:,2)),4)=an_chiR(:,2);

x = 1:4;

figure(3); set(gcf,'Position',[900 200 170 170])
ax = axes();
hold(ax);
for i=1:4
    h=boxchart(x(i)*ones(size(y3(:,i))), y3(:,i), 'BoxFaceColor', colors(i,:),'Notch','on');
    h.MarkerStyle='.';
    h.MarkerColor=colors(i,:);
end
hold on
plot(1:2,nanmedian(y3(:,1:2)),'-o','LineWidth',1,'Color',blue_light,'MarkerSize',3,'MarkerFaceColor',blue_light)
plot(3:4,nanmedian(y3(:,3:4)),'-o','LineWidth',1,'Color',red_light,'MarkerSize',3,'MarkerFaceColor',red_light)
hold off
xlim([0.5 4.5])
ylabel('firing rate ')
xticks(1:4)
xticklabels({'nav','fast nav','com','slow com'});

set(gca, 'box', 'off')
set(gca, 'Color','none')
set(gca,'linewidth',1);set(gca,'fontsize',8);

xt = get(gca, 'XTick');
line1=170;


sl=0.75; 
% plot lines for comparison 
hold on 
plot(xt(1:2), [1 1]*line1, '-k','LineWidth',sl) 
plot([1;1]*xt(1),[line1,line1-3],'-k', 'LineWidth',sl);
plot([1;1]*xt(2),[line1,line1-3],'-k', 'LineWidth',sl);

plot(xt(3:4), [1 1]*line1, '-k','LineWidth',sl) 
plot([1;1]*xt(3),[line1,line1-3],'-k', 'LineWidth',sl);
plot([1;1]*xt(4),[line1,line1-3],'-k', 'LineWidth',sl);

ylim([-5 185])

% plot level of significance
if pv6 <= 0.05 && pv6 > 0.01
    text(mean(xt(1:2))-0.05, line1+4, '*','HorizontalAlignment','center')
elseif pv6 <=0.01 && pv6 > 0.001
    text(mean(xt(1:2))-0.05, line1+4, '**','HorizontalAlignment','center')
elseif pv6 <= 0.001
    text(mean(xt(1:2))-0.05, line1+4, '***','HorizontalAlignment','center')
else 
    text(mean(xt(1:2))-0.05, line1+10, 'ns','FontSize',8,'HorizontalAlignment','center')
end

if pv7 <= 0.05 && pv7 > 0.01
    text(mean(xt(3:4))-0.05, line1+4, '*','HorizontalAlignment','center')
elseif pv7 <=0.01 && pv7 > 0.001
    text(mean(xt(3:4))-0.05, line1+4, '**','HorizontalAlignment','center')
elseif pv7 <= 0.001
    text(mean(xt(3:4))-0.05, line1+4, '***','HorizontalAlignment','center')
else 
    text(mean(xt(3:4))-0.05, line1+10, 'ns','FontSize',8,'HorizontalAlignment','center')
end