close all;
clear;
% plot spiekcount and spontaneous for paired data set 
load('sp_context_pair.mat') 
dur_echo=1.3311;
dur_distress=1.9266;

c_=sp_context(:,2)==0;   % unit classification
aw_=sp_context(:,3)==0;  % awake data
particular1=logical(aw_.*c_);
an_units=sp_context(particular1,1)+1;  % anesthetized units 
particular2=ismember(sp_context(:,1),unique(an_units));

spont_=sp_context(:,7)==0;    % spontaneous firing rate
cont_=sp_context(:,7)==1;     % context firing rate intact

% submatrices 
aw_spont=sp_context(logical(particular1.*spont_),:); 
an_spont=sp_context(logical(particular2.*spont_),:);

aw_cont=sp_context(logical(particular1.*cont_),:);
an_cont=sp_context(logical(particular2.*cont_),:);


%% calculate average firing rate from spike counts across trials

% 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=mean(aw_spont(u_,8));
    rate2=mean(aw_cont(u_2,5));   % nav
    rate3=mean(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 & context
units=unique(an_spont(:,1));
an_spontR=nan(numel(units),1);
an_contR=nan(numel(units),2);
for u=1:numel(units)
    u_=an_spont(:,1)==units(u);
    u_2=an_cont(:,1)==units(u);
    rate=mean(an_spont(u_,8));
    rate2=mean(an_cont(u_2,5));   % nav
    rate3=mean(an_cont(u_2,6));   % com 
    an_spontR(u,1)=rate;   
    an_contR(u,1)=rate2;  
    an_contR(u,2)=rate3; 
end
an_spontR=an_spontR./((3.5-0.2));
an_contR(:,1)=an_contR(:,1);
an_contR(:,2)=an_contR(:,2);

y=[aw_spontR,an_spontR,aw_contR(:,1),an_contR(:,1),aw_contR(:,2),an_contR(:,2)];

% firing rate in response to probes
load('pair_data.mat'); 
c_=all_di(:,2)==0;
aw_=all_di(:,3)==0;
part=logical(c_.*aw_);
an_=[0; part];
an_(end)=[];
probe_aw=all_di(part,18:19).*20;
probe_an=all_di(logical(an_),18:19).*20;

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


%% plot spontaneous and firing during context in same plot  
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;
lines_color=[0.7 0.7 0.7];  % light gray

bar_colors=[[0.3 0.3 0.3];[0.5 0.5 0.5];blue;red;blue_light;red_light];

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(y(:,i))), y(:,i), 'BoxFaceColor', bar_colors(i,:),'Notch','off');
    h.MarkerStyle='.';
    h.MarkerColor=bar_colors(i,:);
end
% add lines per unit
hold on
for i=1:numel(units)
    l=plot([1,2], [aw_spontR(i), an_spontR(i)],'color',lines_color);
    uistack(l,'bottom')
    l=plot([3,4], [aw_contR(i,1), an_contR(i,1)],'color',lines_color);
    uistack(l,'bottom')
    l=plot([5,6], [aw_contR(i,2), an_contR(i,2)],'color',lines_color);
    uistack(l,'bottom')
end
hold off

xlim([0.5 6.5])
ylabel('# spikes')
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=180;
line2=200;
line3=220;

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 180])

% 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',7,'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',7,'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+12, 'ns','FontSize',7,'HorizontalAlignment','center')
end


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

exportgraphics(gcf,'E:\Users\User\Desktop\delay paper\anesthetized\figures_paper\sameunit_inj_rate.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','off');
    h.MarkerStyle='.';
    h.MarkerColor=colors(i,:);
end
% add lines per unit
hold on
for i=1:numel(units)
    l=plot([1,2], [probe_aw(i,1), probe_an(i,1)],'color',lines_color);
    uistack(l,'bottom')
    l=plot([3,4], [probe_aw(i,2), probe_an(i,2)],'color',lines_color);
    uistack(l,'bottom')
end
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=220;
% 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 240])

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

if pv7 <= 0.05 && pv7 > 0.01
    text(mean(xt(3:4))-0.05, line1+10, '*','HorizontalAlignment','center')
elseif pv7 <=0.01 && pv7 > 0.001
    text(mean(xt(3:4))-0.05, line1+10, '**','HorizontalAlignment','center')
elseif pv7 <= 0.001
    text(mean(xt(3:4))-0.05, line1+10, '***','HorizontalAlignment','center')
else 
    text(mean(xt(3:4))-0.05, line1+18, '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\sameunit_rate_probe.pdf',...
    'Resolution',300','ContentType','vector','BackgroundColor','none')