%Plotting decoding analysis
clear all;
load('UnitMdlDSSucrose.mat')
load('UnitMdlDSShuffSucrose.mat')
load('UnitMdlPavSucrose.mat')
load('UnitMdlPavShuffSucrose.mat')
if ~exist('RPPSall'),load('RPPSall.mat'); end
if ~exist('RPSall'),load('RPSall.mat'); end
load('UnitMdlDSAlcohol.mat')
load('UnitMdlDSShuffAlcohol.mat')
load('UnitMdlPavAlcohol.mat')
load('UnitMdlPavShuffAlcohol.mat')
if ~exist('RvppasALL'),load('RvppasALL.mat'); end
if ~exist('RvpppasALL'),load('RvpppasALL.mat'); end

sucrose=[1  0.6  0.1];
maltodextrin=[.9  0.3  .9];
water=[0.00  0.75  0.75];
total=[0.3  0.1  0.8];
DS=[0    0.4470    0.7410];
Pav=[0.8500    0.3250    0.0980];
DSAlcohol=[51/255 160/255 44/255];%[0.4940    0.1840    0.5560];
PavAlcohol=[180/255 41/255 230/255];%[0.9290    0.6940    0.1250];

DSShuff=[0 0 0];
PavShuff=[0 0 0];

xaxis=linspace(-0.75,3.15,13);

selection1=RPPSall.Structure==10 & RPPSall.CSPlusRatio>=.5 & (RPPSall.CSPlusRatio./(RPPSall.CSPlusRatio+RPPSall.CSMinusRatio))>=.7;% & RPPSall.Ev(1).RespDir==1;
selection2=RPSall.Structure==10 & RPSall.CSPlusRatio>=.5 & (RPSall.CSPlusRatio./(RPSall.CSPlusRatio+RPSall.CSMinusRatio))>=.7;% & RPSall.Ev(1).RespDir==1;
selection3=RvpppasALL.Structure==10 & RvpppasALL.CSPlusRatio>=.5 & (RvpppasALL.CSPlusRatio./(RvpppasALL.CSPlusRatio+RvpppasALL.CSMinusRatio))>=.7;% & RvpppasALL.Ev(1).RespDir==1;
selection4=RvppasALL.Structure==10 & RvppasALL.CSPlusRatio>=.5 & (RvppasALL.CSPlusRatio./(RvppasALL.CSPlusRatio+RvppasALL.CSMinusRatio))>=.7;% & RvppasALL.Ev(1).RespDir==1;

%Figure code

figure;

%get average accuracy for each bin
for i = 1:13
    AvgAccDSSucrose(1,i)=nanmean(UnitMdlAccDSSucrose(selection1,i)); %average accuracy DS
    SEMAccDSSucrose(1,i)=nanste(UnitMdlAccDSSucrose(selection1,i),1); %SEM
    AvgAccDSShuffSucrose(1,i)=nanmean(UnitMdlAccDSShuffSucrose(selection1,i)); %average accuracy DSShuff
    SEMAccDSShuffSucrose(1,i)=nanste(UnitMdlAccDSShuffSucrose(selection1,i),1); %SEM
    AvgAccPavSucrose(1,i)=nanmean(UnitMdlAccPavSucrose(selection2,i)); %average accuracy DSShuff
    SEMAccPavSucrose(1,i)=nanste(UnitMdlAccPavSucrose(selection2,i),1); %SEM
    AvgAccPavShuffSucrose(1,i)=nanmean(UnitMdlAccPavShuffSucrose(selection2,i)); %average accuracy DSShuff
    SEMAccPavShuffSucrose(1,i)=nanste(UnitMdlAccPavShuffSucrose(selection2,i),1); %SEM
    
    DSSigSucrose(1,i) = ttest2(UnitMdlAccDSSucrose(:,i),UnitMdlAccDSShuffSucrose(selection1,i),'Vartype','unequal'); %welch's t-test
    PavSigSucrose(1,i) = ttest2(UnitMdlAccPavSucrose(:,i),UnitMdlAccPavShuffSucrose(selection2,i),'Vartype','unequal'); %welch's t-test
end

%get normalized accuracy increase over shuffled for each neuron, then take
%average and find SEM
for i = 1:13
    for j = 1:length(UnitMdlAccDSSucrose)
        NormAccDSSucrose(j,i) = UnitMdlAccDSSucrose(j,i)-nanmean(nanmean(UnitMdlAccDSShuffSucrose(:,i)));
    end
    for k = 1:length(UnitMdlAccPavSucrose)
        NormAccPavSucrose(k,i) = UnitMdlAccPavSucrose(k,i)-nanmean(nanmean(UnitMdlAccPavShuffSucrose(:,i)));
    end
    
    AvgNormDSSucrose(1,i)=nanmean(NormAccDSSucrose(selection1,i)); %average norm accuracy DS
    SEMNormDSSucrose(1,i)=nanste(NormAccDSSucrose(selection1,i),1); %SEM
    AvgNormPavSucrose(1,i)=nanmean(NormAccPavSucrose(selection2,i)); %average norm accuracy DS
    SEMNormPavSucrose(1,i)=nanste(NormAccPavSucrose(selection2,i),1); %SEM
    
    %NormSig(1,i) = ttest2(NormAccDS(:,i),NormAccPav(:,i),'Vartype','unequal'); %welch's t-test
    [~,NormSig(1,i)] = ranksum(NormAccDSSucrose(selection1,i),NormAccPavSucrose(selection2,i)); %wilcoxon's ranksum test
end

%prepare shading
upSEMDS1=AvgAccDSSucrose+SEMAccDSSucrose;
downSEMDS1=AvgAccDSSucrose-SEMAccDSSucrose;
upSEMPav1=AvgAccPavSucrose+SEMAccPavSucrose;
downSEMPav1=AvgAccPavSucrose-SEMAccPavSucrose;
upSEMDSShuff1=AvgAccDSShuffSucrose+SEMAccDSShuffSucrose;
downSEMDSShuff1=AvgAccDSShuffSucrose-SEMAccDSShuffSucrose;
upSEMPavShuff1=AvgAccPavShuffSucrose+SEMAccPavShuffSucrose;
downSEMPavShuff1=AvgAccPavShuffSucrose-SEMAccPavShuffSucrose;
upSEMNormDS1=AvgNormDSSucrose+SEMNormDSSucrose;
downSEMNormDS1=AvgNormDSSucrose-SEMNormDSSucrose;
upSEMNormPav1=AvgNormPavSucrose+SEMNormPavSucrose;
downSEMNormPav1=AvgNormPavSucrose-SEMNormPavSucrose;

% %plotting decoder accuracy over time
% subplot(3,4,1); %accumbens
% hold on;
% plot(xaxis,AvgAccDSSucrose(1:13),'Color', DS,'linewidth',3);
% plot(xaxis,AvgAccDSShuffSucrose(1:13),'Color', DSShuff,'linewidth',3);
% %plot(xaxis,DSSigSucrose-0.51,'*','Color', 'k');
% patch([xaxis,xaxis(end:-1:1)],[upSEMDS,downSEMDS(end:-1:1)],DS,'EdgeColor','none');alpha(0.5);
% patch([xaxis,xaxis(end:-1:1)],[upSEMDSShuff,downSEMDSShuff(end:-1:1)],DSShuff,'EdgeColor','none');alpha(0.5);
% xlabel('Seconds post cue');
% legend('Sucrose DS','Shuffled','Location','best');
% %title('Sucrose DS single unit decoding accuracy');
% ylabel('Single Unit Decoding Accuracy');
% axis([-.75 3.15 0.45 0.7]);
% 
% subplot(3,4,2); %Pav
% hold on;
% plot(xaxis,AvgAccPavSucrose(1:13),'Color', Pav,'linewidth',3);
% plot(xaxis,AvgAccPavShuffSucrose(1:13),'Color', PavShuff,'linewidth',3);
% %plot(xaxis,PavSig-0.51,'*','Color', 'k');
% patch([xaxis,xaxis(end:-1:1)],[upSEMPav,downSEMPav(end:-1:1)],Pav,'EdgeColor','none');alpha(0.5);
% patch([xaxis,xaxis(end:-1:1)],[upSEMPavShuff,downSEMPavShuff(end:-1:1)],PavShuff,'EdgeColor','none');alpha(0.5);
% xlabel('Seconds post cue');
% %title('Sucrose CS+ single unit decoding accuracy');
% legend('Sucrose CS+','Shuffled','Location','northeast');
% axis([-.75 3.15 0.45 0.7]);

%get average accuracy for each bin
for i = 1:13
    AvgAccDSAlcohol(1,i)=nanmean(UnitMdlAccDSAlcohol(selection3,i)); %average accuracy DS
    SEMAccDSAlcohol(1,i)=nanste(UnitMdlAccDSAlcohol(selection3,i),1); %SEM
    AvgAccDSShuffAlcohol(1,i)=nanmean(UnitMdlAccDSShuffAlcohol(selection3,i)); %average accuracy DSShuff
    SEMAccDSShuffAlcohol(1,i)=nanste(UnitMdlAccDSShuffAlcohol(selection3,i),1); %SEM
    AvgAccPavAlcohol(1,i)=nanmean(UnitMdlAccPavAlcohol(selection4,i)); %average accuracy DSShuff
    SEMAccPavAlcohol(1,i)=nanste(UnitMdlAccPavAlcohol(selection4,i),1); %SEM
    AvgAccPavShuffAlcohol(1,i)=nanmean(UnitMdlAccPavShuffAlcohol(selection4,i)); %average accuracy DSShuff
    SEMAccPavShuffAlcohol(1,i)=nanste(UnitMdlAccPavShuffAlcohol(selection4,i),1); %SEM
    
    DSSigAlcohol(1,i) = ttest2(UnitMdlAccDSAlcohol(:,i),UnitMdlAccDSShuffAlcohol(selection3,i),'Vartype','unequal'); %welch's t-test
    PavSigAlcohol(1,i) = ttest2(UnitMdlAccPavAlcohol(:,i),UnitMdlAccPavShuffAlcohol(selection4,i),'Vartype','unequal'); %welch's t-test
end

%get normalized accuracy increase over shuffled for each neuron, then take
%average and find SEM
for i = 1:13
    for j = 1:length(UnitMdlAccDSAlcohol)
        NormAccDSAlcohol(j,i) = UnitMdlAccDSAlcohol(j,i)-nanmean(nanmean(UnitMdlAccDSShuffAlcohol(:,i)));
    end
    for k = 1:length(UnitMdlAccPavAlcohol)
        NormAccPavAlcohol(k,i) = UnitMdlAccPavAlcohol(k,i)-nanmean(nanmean(UnitMdlAccPavShuffAlcohol(:,i)));
    end
    
    AvgNormDSAlcohol(1,i)=nanmean(NormAccDSAlcohol(selection3,i)); %average norm accuracy DS
    SEMNormDSAlcohol(1,i)=nanste(NormAccDSAlcohol(selection3,i),1); %SEM
    AvgNormPavAlcohol(1,i)=nanmean(NormAccPavAlcohol(selection4,i)); %average norm accuracy DS
    SEMNormPavAlcohol(1,i)=nanste(NormAccPavAlcohol(selection4,i),1); %SEM
    
    %NormSig(1,i) = ttest2(NormAccDS(:,i),NormAccPav(:,i),'Vartype','unequal'); %welch's t-test
    [~,NormSig(1,i)] = ranksum(NormAccDSAlcohol(selection3,i),NormAccPavAlcohol(selection4,i)); %wilcoxon's ranksum test
end

%prepare shading
upSEMDS=AvgAccDSAlcohol+SEMAccDSAlcohol;
downSEMDS=AvgAccDSAlcohol-SEMAccDSAlcohol;
upSEMPav=AvgAccPavAlcohol+SEMAccPavAlcohol;
downSEMPav=AvgAccPavAlcohol-SEMAccPavAlcohol;
upSEMDSShuff=AvgAccDSShuffAlcohol+SEMAccDSShuffAlcohol;
downSEMDSShuff=AvgAccDSShuffAlcohol-SEMAccDSShuffAlcohol;
upSEMPavShuff=AvgAccPavShuffAlcohol+SEMAccPavShuffAlcohol;
downSEMPavShuff=AvgAccPavShuffAlcohol-SEMAccPavShuffAlcohol;
upSEMNormDS=AvgNormDSAlcohol+SEMNormDSAlcohol;
downSEMNormDS=AvgNormDSAlcohol-SEMNormDSAlcohol;
upSEMNormPav=AvgNormPavAlcohol+SEMNormPavAlcohol;
downSEMNormPav=AvgNormPavAlcohol-SEMNormPavAlcohol;

% %plotting decoder accuracy over time
% subplot(3,4,5); %accumbens
% hold on;
% plot(xaxis,AvgAccDSAlcohol(1:13),'Color', DSAlcohol,'linewidth',3);
% plot(xaxis,AvgAccDSShuffAlcohol(1:13),'Color', DSShuff,'linewidth',3);
% %plot(xaxis,DSSigAlcohol-0.51,'*','Color', 'k');
% patch([xaxis,xaxis(end:-1:1)],[upSEMDS,downSEMDS(end:-1:1)],DSAlcohol,'EdgeColor','none');alpha(0.5);
% patch([xaxis,xaxis(end:-1:1)],[upSEMDSShuff,downSEMDSShuff(end:-1:1)],DSShuff,'EdgeColor','none');alpha(0.5);
% xlabel('Seconds post cue');
% legend('Suc(alc) DS','Shuffled','Location','northeast');
% %title('Alcohol DS single unit decoding accuracy');
% axis([-.75 3.15 0.45 0.7]);
% ylabel('Single Unit Decoding Accuracy');
% 
% subplot(3,4,6); %Pav
% hold on;
% plot(xaxis,AvgAccPavAlcohol(1:13),'Color', PavAlcohol,'linewidth',3);
% plot(xaxis,AvgAccPavShuffAlcohol(1:13),'Color', PavShuff,'linewidth',3);
% %plot(xaxis,PavSig-0.51,'*','Color', 'k');
% patch([xaxis,xaxis(end:-1:1)],[upSEMPav,downSEMPav(end:-1:1)],PavAlcohol,'EdgeColor','none');alpha(0.5);
% patch([xaxis,xaxis(end:-1:1)],[upSEMPavShuff,downSEMPavShuff(end:-1:1)],PavShuff,'EdgeColor','none');alpha(0.5);
% xlabel('Seconds post cue');
% %title('Alcohol CS+ single unit decoding accuracy');
% legend('Suc(alc) CS+','Shuffled','Location','northeast');
% axis([-.75 3.15 0.45 0.7]);

% combined single unit decoding over time
subplot(2,2,1);
hold on;
plot(xaxis,AvgAccDSAlcohol(1:13),'Marker','o','Color', DSAlcohol,'linewidth',1);
plot(xaxis,AvgAccPavAlcohol(1:13),'Marker','o','Color', PavAlcohol,'linewidth',1);
plot(xaxis,AvgAccDSShuffAlcohol(1:13),'Marker','o','Color', DSShuff,'linewidth',1);
plot(xaxis,AvgAccPavShuffAlcohol(1:13),'Marker','o','Color', PavShuff,'linewidth',1);
%plot(xaxis,DSSigAlcohol-0.51,'*','Color', 'k');
patch([0 0.3 0.3 0],[0 0 1 1],[0.5 0.5 0.5],'EdgeColor','none');alpha(0.5);
patch([xaxis,xaxis(end:-1:1)],[upSEMDS,downSEMDS(end:-1:1)],DSAlcohol,'EdgeColor','none');alpha(0.5);
patch([xaxis,xaxis(end:-1:1)],[upSEMDSShuff,downSEMDSShuff(end:-1:1)],DSShuff,'EdgeColor','none');alpha(0.5);
patch([xaxis,xaxis(end:-1:1)],[upSEMPav,downSEMPav(end:-1:1)],PavAlcohol,'EdgeColor','none');alpha(0.5);
patch([xaxis,xaxis(end:-1:1)],[upSEMPavShuff,downSEMPavShuff(end:-1:1)],PavShuff,'EdgeColor','none');alpha(0.5);


%this adds all 4 plots
% plot(xaxis,AvgAccDSSucrose(1:13),'Color', DS,'linewidth',1);
% plot(xaxis,AvgAccDSShuffSucrose(1:13),'Color', DSShuff,'linewidth',1);
% patch([xaxis,xaxis(end:-1:1)],[upSEMDS1,downSEMDS1(end:-1:1)],DS,'EdgeColor','none');alpha(0.5);
% patch([xaxis,xaxis(end:-1:1)],[upSEMDSShuff1,downSEMDSShuff1(end:-1:1)],DSShuff,'EdgeColor','none');alpha(0.5);
% plot(xaxis,AvgAccPavSucrose(1:13),'Color', Pav,'linewidth',1);
% plot(xaxis,AvgAccPavShuffSucrose(1:13),'Color', PavShuff,'linewidth',1);
% patch([xaxis,xaxis(end:-1:1)],[upSEMPav1,downSEMPav1(end:-1:1)],Pav,'EdgeColor','none');alpha(0.5);
% patch([xaxis,xaxis(end:-1:1)],[upSEMPavShuff1,downSEMPavShuff1(end:-1:1)],PavShuff,'EdgeColor','none');alpha(0.5);

xlabel('Seconds post cue');
legend('Suc(alc) DS','Suc(alc) CS+','Shuffled','Location','northeast');
%title('Alcohol DS single unit decoding accuracy');
axis([-.75 3.15 0.45 0.7]);
ylabel('Single Unit Decoding Accuracy');

%plotting cdf plots for DS at time of max predictive power
subplot(2,2,2);
DSbin=4;
[cdfDS,xDS] = ecdf(UnitMdlAccDSAlcohol(selection3,DSbin));
[cdfPav,xPav] = ecdf(UnitMdlAccPavAlcohol(selection4,DSbin));
[cdfDSsuc,xDSsuc] = ecdf(UnitMdlAccDSSucrose(selection1,DSbin));
[cdfPavsuc,xPavsuc] = ecdf(UnitMdlAccPavSucrose(selection2,DSbin));
%subplot(4,3,[3 6]);
hold on;
plot(xDSsuc,cdfDSsuc,'Color',DS,'linewidth',1);
plot(xPavsuc,cdfPavsuc,'Color',Pav,'linewidth',1);
plot(xDS,cdfDS,'Color',DSAlcohol,'linewidth',1);
plot(xPav,cdfPav,'Color',PavAlcohol,'linewidth',1);
plot([.5 .5], [0 1],'LineStyle',':','Color','k');
plot([AvgAccDSSucrose(1,4) AvgAccDSSucrose(1,4)], [0 1],'Color',DS,'linewidth',1);
plot([AvgAccDSAlcohol(1,4) AvgAccDSAlcohol(1,4)], [0 1],'Color',DSAlcohol,'linewidth',1);
plot([AvgAccPavSucrose(1,4) AvgAccPavSucrose(1,4)], [0 1],'Color',Pav,'linewidth',1);
plot([AvgAccPavAlcohol(1,4) AvgAccPavAlcohol(1,4)], [0 1],'Color',PavAlcohol,'linewidth',1);
axis([0.25 1 0 1]);
title(['Cumulative Distribution 0 to .3s']);
xlabel('Decoding Accuracy')
ylabel('Proportion of Population')
legend('Sucrose DS','Sucrose CS+','Suc(alc) DS','Suc(alc) CS+','Location','best');


%Plotting ensemble decoding

clear all;
DS=[0    0.4470    0.7410];
Pav=[0.8500    0.3250    0.0980];
DSAlcohol=[51/255 160/255 44/255];%[0.4940    0.1840    0.5560];
PavAlcohol=[180/255 41/255 230/255];%[0.9290    0.6940    0.1250];
NumNeurons=[1 5 10 25 50 100];

load('PoolMdlAccDSSucrose.mat','PoolMdlAccDSSucrose')
load('PoolMdlAccDSSucroseShuff.mat','PoolMdlAccDSSucroseShuff')
load('PoolMdlAccPavSucrose.mat','PoolMdlAccPavSucrose')
load('PoolMdlAccPavSucroseShuff.mat','PoolMdlAccPavSucroseShuff')
load('PoolMdlAccDSAlcohol.mat','PoolMdlAccDSAlcohol')
load('PoolMdlAccDSAlcoholShuff.mat','PoolMdlAccDSAlcoholShuff')
load('PoolMdlAccPavAlcohol.mat','PoolMdlAccPavAlcohol')
load('PoolMdlAccPavAlcoholShuff.mat','PoolMdlAccPavAlcoholShuff')


for i=1:6    
AvgAccDSSucrose(i)=nanmean(PoolMdlAccDSSucrose{i,1}); %average accuracy
SEMAccDSSucrose(i)=nanste(PoolMdlAccDSSucrose{i,1},1); %SEM

AvgAccPavSucrose(i)=nanmean(PoolMdlAccPavSucrose{i,1}); %average accuracy
SEMAccPavSucrose(i)=nanste(PoolMdlAccPavSucrose{i,1},1); %SEM

AvgAccDSAlcohol(i)=nanmean(PoolMdlAccDSAlcohol{i,1}); %average accuracy
SEMAccDSAlcohol(i)=nanste(PoolMdlAccDSAlcohol{i,1},1); %SEM

AvgAccPavAlcohol(i)=nanmean(PoolMdlAccPavAlcohol{i,1}); %average accuracy
SEMAccPavAlcohol(i)=nanste(PoolMdlAccPavAlcohol{i,1},1); %SEM

AvgAccDSSucroseShuff(i)=nanmean(PoolMdlAccDSSucroseShuff{i,1}); %average accuracy
SEMAccDSSucroseShuff(i)=nanste(PoolMdlAccDSSucroseShuff{i,1},1); %SEM

AvgAccPavSucroseShuff(i)=nanmean(PoolMdlAccPavSucroseShuff{i,1}); %average accuracy
SEMAccPavSucroseShuff(i)=nanste(PoolMdlAccPavSucroseShuff{i,1},1); %SEM

AvgAccDSAlcoholShuff(i)=nanmean(PoolMdlAccDSAlcoholShuff{i,1}); %average accuracy
SEMAccDSAlcoholShuff(i)=nanste(PoolMdlAccDSAlcoholShuff{i,1},1); %SEM

AvgAccPavAlcoholShuff(i)=nanmean(PoolMdlAccPavAlcoholShuff{i,1}); %average accuracy
SEMAccPavAlcoholShuff(i)=nanste(PoolMdlAccPavAlcoholShuff{i,1},1); %SEM
end



subplot(2,2,3); %ensemble decoding results
errorbar(NumNeurons, AvgAccDSSucrose,SEMAccDSSucrose,'Color',DS,'linewidth',1);
hold on;
errorbar(NumNeurons, AvgAccPavSucrose,SEMAccPavSucrose, 'Color',Pav,'linewidth',1);
errorbar(NumNeurons, AvgAccDSAlcohol,SEMAccDSAlcohol, 'Color',DSAlcohol,'linewidth',1);
errorbar(NumNeurons, AvgAccPavAlcohol,SEMAccPavAlcohol, 'Color',PavAlcohol,'linewidth',1);
% errorbar(NumNeurons, AvgAccDSSucroseShuff,SEMAccDSSucroseShuff, 'Color','k');
% errorbar(NumNeurons, AvgAccPavSucroseShuff,SEMAccPavSucroseShuff, 'Color','k');
% errorbar(NumNeurons, AvgAccDSAlcoholShuff,SEMAccDSAlcoholShuff, 'Color','k');
% errorbar(NumNeurons, AvgAccPavAlcoholShuff,SEMAccPavAlcoholShuff, 'Color','k');
legend('Sucrose DS','Sucrose CS+','Suc(alc) DS','Suc(alc) CS+','Location','best');
plot([-5 105], [.5 .5],'LineStyle',':','Color','k');
xlabel('Neurons used');
ylabel('Mean Accuracy (+/- SEM)');
xticks([1 5 10 25 50 100])
%xticklabels({'-3\pi','-2\pi','-\pi','0','\pi','2\pi','3\pi'})
axis([-5 105 0.4 1]);

%make matrices with bins at all levels
A=[];
B=[];
C=[];
D=[];
A1=[];
B1=[];
C1=[];
D1=[];

for v=1:length(PoolMdlAccDSSucrose)
A=cat(1,A,PoolMdlAccDSSucrose{v,1});
B=cat(1,B,PoolMdlAccPavSucrose{v,1});
C=cat(1,C,PoolMdlAccDSAlcohol{v,1});
D=cat(1,D,PoolMdlAccPavAlcohol{v,1});
A1=cat(1,A1,PoolMdlAccDSSucroseShuff{v,1});
B1=cat(1,B1,PoolMdlAccPavSucroseShuff{v,1});
C1=cat(1,C1,PoolMdlAccDSAlcoholShuff{v,1});
D1=cat(1,D1,PoolMdlAccPavAlcoholShuff{v,1});
end

subplot(2,2,4); %summary of accuracy
errorbar(1,nanmean(A),nanste(A,1),'o','Color',DS,'linewidth',1);
hold on;
errorbar(2,nanmean(C),nanste(C,1),'o', 'Color',DSAlcohol,'linewidth',1);
errorbar(3,nanmean(B),nanste(B,1),'o', 'Color',Pav,'linewidth',1);
errorbar(4,nanmean(D),nanste(D,1),'o', 'Color',PavAlcohol,'linewidth',1);
plot([-5 105], [.5 .5],'LineStyle',':','Color','k');
plot([1 2], [0.93 0.93],'Color','k','linewidth',1);
plot([3 4], [0.86 0.86],'Color','k','linewidth',1);
plot(1.5,.95,'*','Color','k');
plot(3.5,.88,'*','Color','k');
axis([0.5 4.5 0.4 1]);
set(gca,'xtick',[])
ylabel('Mean Accuracy (+/- SEM)');
%xlabel('DS/Suc DS/Suc(a) CS+/Suc CS+/Suc(a)');
legend('Sucrose DS','Suc(alc) DS','Sucrose CS+','Suc(alc) CS+','Location','best');

%% STATS

%Plotting single All decoding
bins = 1; %number of bins
binint = 0.3; %spacing of bins
binstart = 0; %start time of first bin relative to event
Dura=[0 0.3]; %size of bin used for decoding:
NumNeurons=[1 5 10 25 50 100]; %matrix of how many neurons to use on each iteration
repetitions=20; %how many times to run the analysis

%ANOVA for effect of ensemble size, task and reward on accuracy 0-.3 sec
size=[];
for i=1:length(NumNeurons)
    size=cat(1,size,i*ones(repetitions,1));
end
size2=cat(1,size,size,size,size);
%shuffd=cat(1,zeros(length(A),1),zeros(length(B),1),ones(length(C),1),zeros(length(D),1));
reward=cat(1,ones(length(A)+length(B),1),zeros(length(C)+length(D),1));
task=cat(1,ones(length(A),1),zeros(length(B),1),ones(length(C),1),zeros(length(D),1));
[~,AccAnova,AccStats]=anovan(cat(1,A,B,C,D),{reward,task,size2},'varnames',{'reward','task','ens size'},'model','full');
%multcompare(AccStats,'Dimension',[1:3])

rewardtask=cat(1,zeros(length(A),1),ones(length(B),1),2*ones(length(C),1),3*ones(length(D),1));
[~,AccAnova,AccStats]=anovan(cat(1,A,B,C,D),{rewardtask,size2},'varnames',{'reward and task','ens size'},'model','full');


size=[];
for i=1:length(NumNeurons)
    size=cat(1,size,i*ones(repetitions,1));
end
size1=cat(1,size,size,size,size,size,size,size,size);
shuffd=cat(1,zeros(length(A)+length(B)+length(C)+length(D),1),ones(length(A)+length(B)+length(C)+length(D),1));
reward=cat(1,ones(length(A)+length(B),1),zeros(length(C)+length(D),1),ones(length(A)+length(B),1),zeros(length(C)+length(D),1));
task=cat(1,ones(length(A),1),zeros(length(B),1),ones(length(C),1),zeros(length(D),1),ones(length(A),1),zeros(length(B),1),ones(length(C),1),zeros(length(D),1));
[~,AccAnova,AccStats]=anovan(cat(1,A,B,C,D,A1,B1,C1,D1),{reward,task,size1,shuffd},'varnames',{'reward','task','ens size','shuffled vs real'},'model','full');


rewardtask=cat(1,zeros(length(A),1),ones(length(B),1),2*ones(length(C),1),3*ones(length(D),1));
[~,AccAnova,AccStats]=anovan(cat(1,A,B,C,D),{rewardtask,size2},'varnames',{'reward and task','ens size'},'model','full');


%without shuff vs true
% size3=cat(1,size,size);
% region3=cat(1,ones(length(A),1),zeros(length(C),1));
% [~,AccAnovaNoShuff,~]=anovan(cat(1,A,C),{region3,size3},'varnames',{'region','ens size'},'display','off','model','full');

%% Sinlge Unit Stats

if ~exist('RPPSall'),load('RPPSall.mat'); end
if ~exist('RPSall'),load('RPSall.mat'); end
load('UnitMdlDSAlcohol.mat')
load('UnitMdlDSShuffAlcohol.mat')
load('UnitMdlPavAlcohol.mat')
load('UnitMdlPavShuffAlcohol.mat')
if ~exist('RvppasALL'),load('RvppasALL.mat'); end
if ~exist('RvpppasALL'),load('RvpppasALL.mat'); end

sucrose=[1  0.6  0.1];
maltodextrin=[.9  0.3  .9];
water=[0.00  0.75  0.75];
total=[0.3  0.1  0.8];
DS=[0    0.4470    0.7410];
Pav=[0.8500    0.3250    0.0980];
DSAlcohol=[0.4940    0.1840    0.5560];
PavAlcohol=[0.9290    0.6940    0.1250];
DSShuff=[0 0 0];
PavShuff=[0 0 0];

xaxis=linspace(-0.75,3.15,13);

selection1=RPPSall.Structure==10 & RPPSall.CSPlusRatio>=.5 & (RPPSall.CSPlusRatio./(RPPSall.CSPlusRatio+RPPSall.CSMinusRatio))>=.7;
selection2=RPSall.Structure==10 & RPSall.CSPlusRatio>=.5 & (RPSall.CSPlusRatio./(RPSall.CSPlusRatio+RPSall.CSMinusRatio))>=.7;
selection3=RvpppasALL.Structure==10 & RvpppasALL.CSPlusRatio>=.5 & (RvpppasALL.CSPlusRatio./(RvpppasALL.CSPlusRatio+RvpppasALL.CSMinusRatio))>=.7;
selection4=RvppasALL.Structure==10 & RvppasALL.CSPlusRatio>=.5 & (RvppasALL.CSPlusRatio./(RvppasALL.CSPlusRatio+RvppasALL.CSMinusRatio))>=.7;

%LME for the 2 alc exposed
task=cat(1,ones(length(UnitMdlAccDSAlcohol(:,4)),1),ones(length(UnitMdlAccDSAlcohol(:,4)),1),zeros(length(UnitMdlAccPavAlcohol(:,4)),1),zeros(length(UnitMdlAccPavAlcohol(:,4)),1));
shuffd=cat(1,ones(length(UnitMdlAccDSAlcohol(:,4)),1),zeros(length(UnitMdlAccDSAlcohol(:,4)),1),ones(length(UnitMdlAccPavAlcohol(:,4)),1),zeros(length(UnitMdlAccPavAlcohol(:,4)),1));
rat=cat(1,RvpppasALL.Rat,RvpppasALL.Rat,RvppasALL.Rat,RvppasALL.Rat);
tbl=table(cat(1,UnitMdlAccDSAlcohol(:,4),UnitMdlAccDSShuffAlcohol(:,4),UnitMdlAccPavAlcohol(:,4),UnitMdlAccPavShuffAlcohol(:,4)),task,shuffd,rat,'variablenames',{'Accuracy','Task','Shuffled','Subject'});
lme=fitlme(tbl,'Accuracy~Task+Shuffled+(Task:Shuffled)+(1|Subject)');%+(1|Subject:Task:Shuffled)');

%all 4
load('UnitMdlDSSucrose.mat')
load('UnitMdlDSShuffSucrose.mat')
load('UnitMdlPavSucrose.mat')
load('UnitMdlPavShuffSucrose.mat')
if ~exist('RPPSall'),load('RPPSall.mat'); end
if ~exist('RPSall'),load('RPSall.mat'); end
task=cat(1,ones(length(UnitMdlAccDSAlcohol(:,4)),1),ones(length(UnitMdlAccDSAlcohol(:,4)),1),zeros(length(UnitMdlAccPavAlcohol(:,4)),1),zeros(length(UnitMdlAccPavAlcohol(:,4)),1),ones(length(UnitMdlAccDSSucrose(:,4)),1),ones(length(UnitMdlAccDSSucrose(:,4)),1),zeros(length(UnitMdlAccPavSucrose(:,4)),1),zeros(length(UnitMdlAccPavSucrose(:,4)),1));
shuffd=cat(1,ones(length(UnitMdlAccDSAlcohol(:,4)),1),zeros(length(UnitMdlAccDSAlcohol(:,4)),1),ones(length(UnitMdlAccPavAlcohol(:,4)),1),zeros(length(UnitMdlAccPavAlcohol(:,4)),1),ones(length(UnitMdlAccDSSucrose(:,4)),1),zeros(length(UnitMdlAccDSSucrose(:,4)),1),ones(length(UnitMdlAccPavSucrose(:,4)),1),zeros(length(UnitMdlAccPavSucrose(:,4)),1));
rat=cat(1,RvpppasALL.Rat,RvpppasALL.Rat,RvppasALL.Rat,RvppasALL.Rat,RPPSall.Rat,RPPSall.Rat,RPSall.Rat,RPSall.Rat);
exposure=cat(1,ones(length(UnitMdlAccDSAlcohol(:,4)),1),ones(length(UnitMdlAccDSAlcohol(:,4)),1),ones(length(UnitMdlAccPavAlcohol(:,4)),1),ones(length(UnitMdlAccPavAlcohol(:,4)),1),zeros(length(UnitMdlAccDSSucrose(:,4)),1),zeros(length(UnitMdlAccDSSucrose(:,4)),1),zeros(length(UnitMdlAccPavSucrose(:,4)),1),zeros(length(UnitMdlAccPavSucrose(:,4)),1));
tbl=table(cat(1,UnitMdlAccDSAlcohol(:,4),UnitMdlAccDSShuffAlcohol(:,4),UnitMdlAccPavAlcohol(:,4),UnitMdlAccPavShuffAlcohol(:,4),UnitMdlAccDSSucrose(:,4),UnitMdlAccDSShuffSucrose(:,4),UnitMdlAccPavSucrose(:,4),UnitMdlAccPavShuffSucrose(:,4)),task,shuffd,rat,exposure,'variablenames',{'Accuracy','Task','Shuffled','Subject','Exposure'});
lme=fitlme(tbl,'Accuracy~Task+Exposure+Task:Exposure+(1|Subject)');%+(1|Subject:Task:Exposure)');

% [~,AccAnova,AccStats]=anovan(cat(1,UnitMdlAccDSAlcohol(:,4),UnitMdlAccDSShuffAlcohol(:,4),UnitMdlAccPavAlcohol(:,4),UnitMdlAccPavShuffAlcohol(:,4)),{task,shuffd,rat},'varnames',{'task','shuff vs true','subject'},'random',[3],'nested',[0 0 0;0 0 0;1 0 0],'model','full');
% %no subject
% [~,AccAnova,AccStats]=anovan(cat(1,UnitMdlAccDSAlcohol(:,4),UnitMdlAccDSShuffAlcohol(:,4),UnitMdlAccPavAlcohol(:,4),UnitMdlAccPavShuffAlcohol(:,4)),{task,shuffd},'varnames',{'task','shuff vs true',},'model','full');
% %no shuff vs true
% fitlme
% task=cat(1,ones(length(UnitMdlAccDSAlcohol(:,4)),1),zeros(length(UnitMdlAccPavAlcohol(:,4)),1));
% rat=cat(1,RvpppasALL.Rat,RvppasALL.Rat);
% [~,AccAnova,AccStats]=anovan(cat(1,UnitMdlAccDSAlcohol(:,4),UnitMdlAccPavAlcohol(:,4)),{task,rat},'varnames',{'task','subject'},'random',[2],'nested',[0 0;1 0],'model','full');