clear all;
load ('R_2R.mat');

%get parameters
BinBase=R_2R.Param.BinBase;
BinDura=R_2R.Param.BinDura;
bins=R_2R.Param.bins;
binint=R_2R.Param.binint;
binstart=R_2R.Param.binstart;
PStatBins=0.01; %using more stringent cutoff to reduce pre-delivery noise

%divide neurons up by region
NAneurons=strcmp(R_2R.Ninfo(:,3),'NA');
VPneurons=strcmp(R_2R.Ninfo(:,3),'VP');

%which bins bound the area examined for reward-selectivity? (in seconds)
Time1=0.4; %seconds
Time2=3; %seconds
Bin1=(((Time1-BinDura(2)/2)-binstart)/binint); %convert to bin name
Bin2=(((Time2-BinDura(2)/2)-binstart)/binint); %convert to bin name

%sorting bin -- which bin the neurons' activity is sorted on for heatmap(in seconds)
SortBinTime=1.1; %seconds
SortBin=round((((SortBinTime-BinDura(2)/2)-binstart)/binint)); %convert to bin name

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];
exc=[0 113/255 188/255];
inh=[240/255 0 50/255];
NAc=[0.5  0.1  0.8];
VP=[0.3  0.7  0.1];

%% bins analysis

%first and last bin aren't included because you can't compare to the previous/subsequent bin
%this axis plots the bins on their centers
xaxis=linspace(binstart+binint+BinDura(2)/2,binstart+(bins-2)*binint+BinDura(2)/2,bins-2);

for i=2:(bins-1) %no including first or last bin because can't compare to previous/subsequent bin
    
    %finds out whether firing is stronger (high excitation or lower inhibition) for 1 or 2
    for k=1:length(R_2R.Ninfo) %runs through neurons
        if R_2R.BinStatRwrd{i,1}(k).IntSig < PStatBins %if neuron is significant for this bin
            if R_2R.BinStatRwrd{i-1,1}(k).IntSig < PStatBins || R_2R.BinStatRwrd{i+1,1}(k).IntSig < PStatBins %either previous or subsequent bin must be significant
                if R_2R.BinStatRwrd{i,1}(k).R1Mean > R_2R.BinStatRwrd{i,1}(k).R2Mean %if firing is greater on sucrose trials than maltodextrin
                    R_2R.BinRewPref{i,1}(k,1)=1;
                else
                    R_2R.BinRewPref{i,1}(k,1)=2; %otherwise maltodextrin must be greater than sucrose
                end
            else
                R_2R.BinRewPref{i,1}(k,1)=0; %if not significant in 2 consecutive bins
            end
        else
            R_2R.BinRewPref{i,1}(k,1)=0; %if no sig reward modulation
        end
        
    end
    %find how many NAc neurons have significant reward modulation in each bin
    NN1perBin(i,1)=sum(R_2R.BinRewPref{i,1}(NAneurons,1)==1); %sucrose pref
    NN2perBin(i,1)=sum(R_2R.BinRewPref{i,1}(NAneurons,1)==2); %malto pref
    NNperBin(i,1)=sum(R_2R.BinRewPref{i,1}(NAneurons,1)>0); %either
    %normalize to number of neurons in population
    NN1norm=NN1perBin./sum(NAneurons);
    NN2norm=NN2perBin./sum(NAneurons);
    NNnorm=NNperBin./sum(NAneurons);
    
    %find how many VP neurons have significant reward modulation in each bin
    NV1perBin(i,1)=sum(R_2R.BinRewPref{i,1}(VPneurons,1)==1); %sucrose pref
    NV2perBin(i,1)=sum(R_2R.BinRewPref{i,1}(VPneurons,1)==2); %malto pref
    NVperBin(i,1)=sum(R_2R.BinRewPref{i,1}(VPneurons,1)>0); %either
    %normalize to number of neurons in population
    NV1norm=NV1perBin./sum(VPneurons);
    NV2norm=NV2perBin./sum(VPneurons);
    NVnorm=NVperBin./sum(VPneurons);
   

end

%Cumulative reward selectivity
cumsel=zeros(length(R_2R.Ninfo),bins); %set up result matrix
cumsel1=zeros(length(R_2R.Ninfo),bins); %set up result matrix for sucrose
cumsel2=zeros(length(R_2R.Ninfo),bins); %set up result matrix for maltodextrin
%note, this has to be corrected to +1 because of the way the bin matrix is
%layed out (bin 0 is the first column, not the 0th column)
for i=1:length(R_2R.Ninfo) %go through each neuron
    for j=2:bins-1 %look in each bin we did analysis on
        if R_2R.BinRewPref{j,1}(i,1)>0 || cumsel(i,j-1)==1
            cumsel(i,j) = 1;
        end
        if R_2R.BinRewPref{j,1}(i,1)==1 || cumsel1(i,j-1)==1
            cumsel1(i,j) = 1;
        end
        if R_2R.BinRewPref{j,1}(i,1)==2 || cumsel2(i,j-1)==1
            cumsel2(i,j) = 1;
        end        
    end
end

%plotting number of significantly modulated neurons across time
%NAc
subplot(2,3,1);
hold on;
plot(xaxis,NNnorm(2:bins-1),'Color', total,'linewidth',1.5);
plot(xaxis,NN1norm(2:bins-1),'Color',sucrose,'linewidth',1.5);
plot(xaxis,NN2norm(2:bins-1),'Color',maltodextrin,'linewidth',1.5);
plot([Time1 Time1],[-1 1],':','color','k','linewidth',0.75);
plot([Time2 Time2],[-1 1],':','color','k','linewidth',0.75);
axis([xaxis(1) xaxis(end) 0 0.4]);
legend('Total','Suc > mal','Mal > suc','Location','northeast')
ylabel('Fraction of population');
title('NAc reward-selective neurons');

%VP
subplot(2,3,2);
hold on;
plot(xaxis,NVnorm(2:bins-1),'Color', total,'linewidth',1.5);
plot(xaxis,NV1norm(2:bins-1),'Color',sucrose,'linewidth',1.5);
plot(xaxis,NV2norm(2:bins-1),'Color',maltodextrin,'linewidth',1.5);
plot([Time1 Time1],[-1 1],':','color','k','linewidth',0.75);
plot([Time2 Time2],[-1 1],':','color','k','linewidth',0.75);
axis([xaxis(1) xaxis(end) 0 0.4]);
legend('Total','Suc > mal','Mal > suc','Location','northeast')
ylabel('Fraction of population');
title('VP reward-selective neurons');

%plotting cumulative selectivity in each region
%NAc
subplot(2,3,4);
hold on;
plot(xaxis,sum(cumsel(NAneurons,2:bins-1),1)/sum(NAneurons),'Color', total,'linewidth',1.5);
plot(xaxis,sum(cumsel1(NAneurons,2:bins-1),1)/sum(NAneurons),'Color', sucrose,'linewidth',1.5);
plot(xaxis,sum(cumsel2(NAneurons,2:bins-1),1)/sum(NAneurons),'Color', maltodextrin,'linewidth',1.5);
plot([Time1 Time1],[0 1],':','color','k','linewidth',0.75);
plot([Time2 Time2],[0 1],':','color','k','linewidth',0.75);
axis([xaxis(1) xaxis(end) 0 0.75]);
xlabel('Seconds from reward delivery');
ylabel('Cumulative fraction');

%VP
subplot(2,3,5);
hold on;
plot(xaxis,sum(cumsel(VPneurons,2:bins-1),1)/sum(VPneurons),'Color', total,'linewidth',1.5);
plot(xaxis,sum(cumsel1(VPneurons,2:bins-1),1)/sum(VPneurons),'Color', sucrose,'linewidth',1.5);
plot(xaxis,sum(cumsel2(VPneurons,2:bins-1),1)/sum(VPneurons),'Color', maltodextrin,'linewidth',1.5);
plot([Time1 Time1],[0 1],':','color','k','linewidth',0.75);
plot([Time2 Time2],[0 1],':','color','k','linewidth',0.75);
axis([xaxis(1) xaxis(end) 0 0.75]);
xlabel('Seconds from RD');
ylabel('Cumulative fraction');

%plot difference between the two regions for each measure
%reward selective neurons per bin
subplot(2,3,3);
hold on;
plot(xaxis,NNnorm(2:bins-1)-NVnorm(2:bins-1),'Color', total,'linewidth',1.5);
plot(xaxis,NN1norm(2:bins-1)-NV1norm(2:bins-1),'Color',sucrose,'linewidth',1.5);
plot(xaxis,NN2norm(2:bins-1)-NV2norm(2:bins-1),'Color',maltodextrin,'linewidth',1.5);
axis([xaxis(1) xaxis(end) -0.3 0.3]);
plot([xaxis(1) xaxis(end)],[0 0],'color','k','linewidth',0.75);
plot([Time1 Time1],[-1 1],':','color','k','linewidth',0.75);
plot([Time2 Time2],[-1 1],':','color','k','linewidth',0.75);
legend('Total','Suc > mal','Mal > suc','Location','northeast')
title('Difference (NAc - VP)');
ylabel('Fraction of population');

%cumulative selectivity
subplot(2,3,6);
hold on;
plot(xaxis,sum(cumsel(NAneurons,2:bins-1),1)/sum(cumsel(NAneurons,bins-1)),'Color', NAc,'linewidth',1.5);
plot(xaxis,sum(cumsel(VPneurons,2:bins-1),1)/sum(cumsel(VPneurons,bins-1)),'Color', VP,'linewidth',1.5);
plot([Time1 Time1],[-1 1],':','color','k','linewidth',0.75);
plot([Time2 Time2],[-1 1],':','color','k','linewidth',0.75);
axis([xaxis(1) xaxis(end) 0 1]);
legend('NAc','VP','location','southeast');
xlabel('Seconds from RD');
ylabel('Cumulative distribution');
title('Reward-selective onsets');

%statistical test comparing onset of reward-selective responses
for i=1:length(cumsel)
    if sum(cumsel(i,:))>0
        onsetcolumn=min(find(cumsel(i,:)));
        onset(i,1)=round((onsetcolumn-1)*binint+binstart+BinDura(2)/2,1);
    else
        onset(i,1)=NaN;
    end
end        
[~,R_2R.OnsetSig{1,1},R_2R.OnsetSig{1,2}]=anovan(onset,{NAneurons,R_2R.Ninfo(:,4)},'varnames',{'Region','Subject'},'random',[2],'nested',[0 0;1 0],'display','off','model','full');
if cell2mat(R_2R.OnsetSig{1,1}(2,7))<0.05
    plot((xaxis(end)-xaxis(1))/2+xaxis(1),0.95,'*','color','k','linewidth',1);
end

%% individual rat data
figure;

%split data up into individual rats
for i=1:length(R_2R.Ninfo)
    A=char(R_2R.Ninfo(i,4));
    Rats(i,:)=cellstr(A);
end

C=unique(Rats);

for i=1:length(C)
    Sel=strcmp(C(i),Rats);
    for j=2:bins-1
        NRat1(j,1)=sum(R_2R.BinRewPref{j,1}(Sel,1)==1)/sum(Sel);
        NRat2(j,1)=sum(R_2R.BinRewPref{j,1}(Sel,1)==2)/sum(Sel);
        NRat(j,1)=sum(R_2R.BinRewPref{j,1}(Sel,1)>0)/sum(Sel);
    end
    
    %graphing
    subplot(4,3,i);
    hold on;
    plot(xaxis,NRat(2:bins-1),'Color', total,'linewidth',1.5);
    plot(xaxis,NRat1(2:65),'Color',sucrose,'linewidth',1.5);
    plot(xaxis,NRat2(2:65),'Color',maltodextrin,'linewidth',1.5);
    xlabel('Seconds from reward delivery');
    title([cell2mat(C(i)) ' (n=' num2str(sum(Sel)) ')'])
    axis([xaxis(1) xaxis(end) 0 0.5]);
    legend('Total neurons','Suc > mal Hz','Mal > suc Hz','Location','northeast')
    ylabel('Fraction of total population');

end

%% plotting sucrose-selective neurons

figure;

%color map
[magma,inferno,plasma,viridis]=colormaps;
colormap(plasma);
c=[-100 2000];ClimE=sign(c).*abs(c).^(1/4);%colormap

%events we're looking at
ev1=strcmp('RD1', R_2R.Erefnames);
ev2=strcmp('RD2', R_2R.Erefnames);

%setting up parameters
Xaxis=[-2 5];
inttime=find(R_2R.Param.Tm>=Xaxis(1) & R_2R.Param.Tm<=Xaxis(2));
paramtime=R_2R.Param.Tm(inttime);

%find all neurons with greater firing for sucrose
for i = 1:(Bin2-Bin1+1)
    %the added +1 is necessary because bin 20 is the 21st entry in the matrix
    Pref1(:,i)=R_2R.BinRewPref{Bin1+i}==1; %get neurons that have greater firing for sucrose in any of the bins bounded above
    Resp11(:,i)=Pref1(:,i)&cat(1,R_2R.BinStatRwrd{Bin1+i,1}.SucRespDir)==1; %get neurons with excitation to sucrose
    Resp12(:,i)=Pref1(:,i)&cat(1,R_2R.BinStatRwrd{Bin1+i,1}.MalRespDir)==-1;%get neurons with inhibition to maltodextrin
end

Sel=sum(Pref1,2)>0;  %all neurons selective for sucrose in any bin
Sel2=sum(Resp11,2)>0; %all selective neurons sucrose excited in any bin
Sel3=sum(Resp12,2)>0; %all selective neurons malto inhibited in any bin

%-----------------NAc--------------------------------------------
subplot(4,40,[15:23]); %heatmap of sucrose preferring neurons' response to sucrose
Neurons=R_2R.Ev(ev1).PSTHz(Sel&NAneurons,inttime); %get the firing rates of neurons of interest
SucResp=cat(1,R_2R.BinStatRwrd{SortBin+1,1}.R1Mean); %sucrose responses
TMP=SucResp(Sel&NAneurons); %find the magnitude of the excitations for sucrose for this bin
TMP(isnan(TMP))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
[~,SORTimg]=sort(TMP);
Neurons=Neurons(SORTimg,:); %sort the neurons by magnitude
imagesc(paramtime,[1,sum(Sel&NAneurons,1)],Neurons,ClimE);
title(['Sucrose trials']);
ylabel('Individual neurons');
hold on;
plot([0 0],[0 sum(Sel&NAneurons)],':','color','k','linewidth',0.75);

subplot(4,40,[27:35]); %heatmap of sucrose preferring neurons' response to maltodextrin
Neurons=R_2R.Ev(ev2).PSTHz(Sel&NAneurons,inttime); %get the firing rates of neurons of interest
Neurons=Neurons(SORTimg,:); %sort the neurons same as before
imagesc(paramtime,[1,sum(Sel&NAneurons,1)],Neurons,ClimE);
title(['Maltodextrin trials']);
hold on;
plot([0 0],[0 sum(Sel&NAneurons)],':','color','k','linewidth',0.75);

%plot sucrose preferring neurons to sucrose
psth1=nanmean(R_2R.Ev(ev1).PSTHz(Sel&NAneurons,inttime),1); 
sem1=nanste(R_2R.Ev(ev1).PSTHz(Sel&NAneurons,inttime),1); %calculate standard error of the mean
up1=psth1+sem1;
down1=psth1-sem1;

%plot sucrose preferring neurons to malt
psth2=nanmean(R_2R.Ev(ev2).PSTHz(Sel&NAneurons,inttime),1); 
sem2=nanste(R_2R.Ev(ev2).PSTHz(Sel&NAneurons,inttime),1); %calculate standard error of the mean
up2=psth2+sem2;
down2=psth2-sem2;

%plotting
subplot(4,3,1);
title(['Suc>mal (n=' num2str(sum(Sel&NAneurons)) ' of ' num2str(sum(NAneurons)) ')'])
hold on;
plot(paramtime,psth1,'Color',sucrose,'linewidth',1);
plot(paramtime,psth2,'Color',maltodextrin,'linewidth',1);
patch([paramtime,paramtime(end:-1:1)],[up1,down1(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
patch([paramtime,paramtime(end:-1:1)],[up2,down2(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
plot([-2 5],[0 0],':','color','k','linewidth',0.75);
plot([0 0],[-1 10],':','color','k','linewidth',0.75);
axis([-2 5 -1 4]);
ylabel('Mean firing (z-score)');
legend('Suc trials','Mal trials');

%display inhibitions and excitations
both = sum(Sel2&Sel3&NAneurons);
excite = sum(Sel2&NAneurons)-both;
inhib = sum(Sel3&NAneurons)-both;

subplot(4,40,40);
b = bar([inhib both excite; 0 0 0],'stacked');
b(1,1).FaceColor=maltodextrin;
b(1,2).FaceColor=total;
b(1,3).FaceColor=sucrose;
axis([1 1.2 0 both+excite+inhib]);
ylabel('Inh / Both / Excited');
set(gca,'xtick',[])


%----------------VP----------------
subplot(4,40,[80+15:80+23]); %heatmap of sucrose preferring neurons' response to sucrose
Neurons=R_2R.Ev(ev1).PSTHz(Sel&VPneurons,inttime); %get the firing rates of neurons of interest
SucResp=cat(1,R_2R.BinStatRwrd{SortBin+1,1}.R1Mean); %sucrose responses
TMP=SucResp(Sel&VPneurons); %find the magnitude of the excitations for sucrose for this bin
TMP(isnan(TMP))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
[~,SORTimg]=sort(TMP);
Neurons=Neurons(SORTimg,:); %sort the neurons by magnitude
imagesc(paramtime,[1,sum(Sel&VPneurons,1)],Neurons,ClimE);
title(['Sucrose trials']);
ylabel('Individual neurons');
hold on;
plot([0 0],[0 sum(Sel&VPneurons)],':','color','k','linewidth',0.75);

subplot(4,40,[80+27:80+35]); %heatmap of sucrose preferring neurons' response to maltodextrin
Neurons=R_2R.Ev(ev2).PSTHz(Sel&VPneurons,inttime); %get the firing rates of neurons of interest
Neurons=Neurons(SORTimg,:); %sort the neurons same as before
imagesc(paramtime,[1,sum(Sel&VPneurons,1)],Neurons,ClimE);
title(['Maltodextrin trials']);
hold on;
plot([0 0],[0 sum(Sel&VPneurons)],':','color','k','linewidth',0.75);

%plot sucrose preferring neurons to sucrose
psth1=nanmean(R_2R.Ev(ev1).PSTHz(Sel&VPneurons,inttime),1); 
sem1=nanste(R_2R.Ev(ev1).PSTHz(Sel&VPneurons,inttime),1); %calculate standard error of the mean
up1=psth1+sem1;
down1=psth1-sem1;

%plot sucrose preferring neurons to malt
psth2=nanmean(R_2R.Ev(ev2).PSTHz(Sel&VPneurons,inttime),1); 
sem2=nanste(R_2R.Ev(ev2).PSTHz(Sel&VPneurons,inttime),1); %calculate standard error of the mean
up2=psth2+sem2;
down2=psth2-sem2;

%plotting
subplot(4,3,7);
title(['Suc>mal (n=' num2str(sum(Sel&VPneurons)) ' of ' num2str(sum(VPneurons)) ')'])
hold on;
plot(paramtime,psth1,'Color',sucrose,'linewidth',1);
plot(paramtime,psth2,'Color',maltodextrin,'linewidth',1);
patch([paramtime,paramtime(end:-1:1)],[up1,down1(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
patch([paramtime,paramtime(end:-1:1)],[up2,down2(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
plot([-2 5],[0 0],':','color','k','linewidth',0.75);
plot([0 0],[-1 10],':','color','k','linewidth',0.75);
axis([-2 5 -0.6 4]);
ylabel('Mean firing (z-score)');
legend('Suc trials','Mal trials');

%display inhibitions and excitations
both = sum(Sel2&Sel3&VPneurons);
excite = sum(Sel2&VPneurons)-both;
inhib = sum(Sel3&VPneurons)-both;

subplot(4,40,120);
b = bar([inhib both excite; 0 0 0],'stacked');
b(1,1).FaceColor=maltodextrin;
b(1,2).FaceColor=total;
b(1,3).FaceColor=sucrose;
axis([1 1.2 0 both+excite+inhib]);
ylabel('Inh / Both / Excited');
set(gca,'xtick',[])

%% plotting maltodextrin-selective neurons

%find all neurons with greater firing for maltodextrin
for i = 1:(Bin2-Bin1+1)
    %the added +1 is necessary because bin 20 is the 21st entry in the matrix
    Pref2(:,i)=R_2R.BinRewPref{Bin1+i}==2; %get neurons that have greater firing for maltodextrin in any of the bins bounded above
    Resp21(:,i)=Pref2(:,i)&cat(1,R_2R.BinStatRwrd{Bin1+i,1}.SucRespDir)==-1; %get neurons with inhibition to sucrose
    Resp22(:,i)=Pref2(:,i)&cat(1,R_2R.BinStatRwrd{Bin1+i,1}.MalRespDir)==1;%get neurons with excitation to maltodextrin
end

Mel=sum(Pref2,2)>0;  %all neurons selective for malotdextrin in any bin
Mel2=sum(Resp21,2)>0; %all selective neurons sucrose inhibited in any bin
Mel3=sum(Resp22,2)>0; %all selective neurons malto excited in any bin

%-----------------NAc--------------------------------------------
subplot(4,40,[40+15:40+23]); %heatmap of maltodextrin preferring neurons' response to sucrose
Neurons=R_2R.Ev(ev1).PSTHz(Mel&NAneurons,inttime); %get the firing rates of neurons of interest
MalResp=cat(1,R_2R.BinStatRwrd{SortBin+1,1}.R2Mean); %sucrose responses
TMP=MalResp(Mel&NAneurons); %find the magnitude of the excitations for sucrose for this bin
TMP(isnan(TMP))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
[~,SORTimg]=sort(TMP);
Neurons=Neurons(SORTimg,:); %sort the neurons by magnitude
imagesc(paramtime,[1,sum(Mel&NAneurons,1)],Neurons,ClimE);
title(['Sucrose trials']);
xlabel('Seconds from RD');
ylabel('Individual neurons');
hold on;
plot([0 0],[0 sum(Mel&NAneurons)],':','color','k','linewidth',0.75);

subplot(4,40,[40+27:40+35]); %heatmap of maltodextrin preferring neurons' response to maltodextrin
Neurons=R_2R.Ev(ev2).PSTHz(Mel&NAneurons,inttime); %get the firing rates of neurons of interest
Neurons=Neurons(SORTimg,:); %sort the neurons same as before
imagesc(paramtime,[1,sum(Mel&NAneurons,1)],Neurons,ClimE);
title(['Maltodextrin trials']);
xlabel('Seconds from RD');
hold on;
plot([0 0],[0 sum(Mel&NAneurons)],':','color','k','linewidth',0.75);

%plot sucrose preferring neurons to sucrose
psth1=nanmean(R_2R.Ev(ev1).PSTHz(Mel&NAneurons,inttime),1); 
sem1=nanste(R_2R.Ev(ev1).PSTHz(Mel&NAneurons,inttime),1); %calculate standard error of the mean
up1=psth1+sem1;
down1=psth1-sem1;

%plot sucrose preferring neurons to malt
psth2=nanmean(R_2R.Ev(ev2).PSTHz(Mel&NAneurons,inttime),1); 
sem2=nanste(R_2R.Ev(ev2).PSTHz(Mel&NAneurons,inttime),1); %calculate standard error of the mean
up2=psth2+sem2;
down2=psth2-sem2;

%plotting
subplot(4,3,4);
title(['Mal>suc (n=' num2str(sum(Mel&NAneurons)) ' of ' num2str(sum(NAneurons)) ')'])
hold on;
plot(paramtime,psth1,'Color',sucrose,'linewidth',1);
plot(paramtime,psth2,'Color',maltodextrin,'linewidth',1);
patch([paramtime,paramtime(end:-1:1)],[up1,down1(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
patch([paramtime,paramtime(end:-1:1)],[up2,down2(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
plot([-2 5],[0 0],':','color','k','linewidth',0.75);
plot([0 0],[-2 7],':','color','k','linewidth',0.75);
axis([-2 5 -1.2 3]);
ylabel('Mean firing (z-score)');
xlabel('Seconds from RD');
legend('Suc trials','Mal trials','Location','northeast');

%display inhibitions and excitations
both = sum(Mel2&Mel3&NAneurons);
excite = sum(Mel3&NAneurons)-both;
inhib = sum(Mel2&NAneurons)-both;

subplot(4,40,80);
b = bar([inhib both excite; 0 0 0],'stacked');
b(1,1).FaceColor=sucrose;
b(1,2).FaceColor=total;
b(1,3).FaceColor=maltodextrin;
axis([1 1.2 0 both+excite+inhib]);
ylabel('Inh / Both / Excited');
set(gca,'xtick',[]);


%----------------VP----------------
subplot(4,40,[120+15:120+23]); %heatmap of maltodextrin preferring neurons' response to sucrose
Neurons=R_2R.Ev(ev1).PSTHz(Mel&VPneurons,inttime); %get the firing rates of neurons of interest
MalResp=cat(1,R_2R.BinStatRwrd{SortBin+1,1}.R2Mean); %sucrose responses
TMP=MalResp(Mel&VPneurons); %find the magnitude of the excitations for sucrose for this bin
TMP(isnan(TMP))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
[~,SORTimg]=sort(TMP);
Neurons=Neurons(SORTimg,:); %sort the neurons by magnitude
imagesc(paramtime,[1,sum(Mel&VPneurons,1)],Neurons,ClimE);
title(['Sucrose trials']);
xlabel('Seconds from RD');
ylabel('Individual neurons');
hold on;
plot([0 0],[0 sum(Mel&VPneurons)],':','color','k','linewidth',0.75);

subplot(4,40,[120+27:120+35]); %heatmap of maltodextrin preferring neurons' response to maltodextrin
Neurons=R_2R.Ev(ev2).PSTHz(Mel&VPneurons,inttime); %get the firing rates of neurons of interest
Neurons=Neurons(SORTimg,:); %sort the neurons same as before
imagesc(paramtime,[1,sum(Mel&VPneurons,1)],Neurons,ClimE);
title(['Maltodextrin trials']);
xlabel('Seconds from RD');
hold on;
plot([0 0],[0 sum(Mel&VPneurons)],':','color','k','linewidth',0.75);

%plot maltodextrin preferring neurons to sucrose
psth1=nanmean(R_2R.Ev(ev1).PSTHz(Mel&VPneurons,inttime),1); 
sem1=nanste(R_2R.Ev(ev1).PSTHz(Mel&VPneurons,inttime),1); %calculate standard error of the mean
up1=psth1+sem1;
down1=psth1-sem1;

%plot malt preferring neurons to malt
psth2=nanmean(R_2R.Ev(ev2).PSTHz(Mel&VPneurons,inttime),1); 
sem2=nanste(R_2R.Ev(ev2).PSTHz(Mel&VPneurons,inttime),1); %calculate standard error of the mean
up2=psth2+sem2;
down2=psth2-sem2;

%plotting
subplot(4,3,10);
title(['Mal>suc (n=' num2str(sum(Mel&VPneurons)) ' of ' num2str(sum(VPneurons)) ')'])
hold on;
plot(paramtime,psth1,'Color',sucrose,'linewidth',1);
plot(paramtime,psth2,'Color',maltodextrin,'linewidth',1);
patch([paramtime,paramtime(end:-1:1)],[up1,down1(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
patch([paramtime,paramtime(end:-1:1)],[up2,down2(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
plot([-2 5],[0 0],':','color','k','linewidth',0.75);
plot([0 0],[-2 7],':','color','k','linewidth',0.75);
axis([-2 5 -1.5 2.5]);
ylabel('Mean firing (z-score)');
xlabel('Seconds from RD');
legend('Suc trials','Mal trials','Location','northeast');

%display inhibitions and excitations
both = sum(Mel2&Mel3&VPneurons);
excite = sum(Mel3&VPneurons)-both;
inhib = sum(Mel2&VPneurons)-both;

subplot(4,40,160);
b = bar([inhib both excite; 0 0 0],'stacked');
b(1,1).FaceColor=sucrose;
b(1,2).FaceColor=total;
b(1,3).FaceColor=maltodextrin;
axis([1 1.2 0 both+excite+inhib]);
ylabel('Inh / Both / Excited');
set(gca,'xtick',[]);

%% creating and saving new variables
R_2R.RSinfo=R_2R.Ninfo(Sel|Mel,1:2); %reward-selective neurons
R_2R.RNSinfo=R_2R.Ninfo(Sel==0 & Mel==0,1:2); %reward-nonselective neurons
R_2R.SucN=Sel;
R_2R.MalN=Mel;
R_2R.Param.SelectiveBins=[Bin1 Bin2];


% proportion summary and chi squared test for reward selective neurons
R_2R.RegSelX2(1,1)=sum((Sel | Mel) & NAneurons)/sum(NAneurons); %how many selective NAc neurons
R_2R.RegSelX2(1,2)=sum((Sel | Mel) & VPneurons)/sum(VPneurons); %how many selective VP neurons
[~,R_2R.RegSelX2(2,1),R_2R.RegSelX2(2,2)]=crosstab(Sel | Mel,NAneurons);

% proportion summary and chi squared test for max neurons in any bin
R_2R.RegMaxSelX2(1,1)=max(NNnorm); %how many selective NAc neurons
R_2R.RegMaxSelX2(1,2)=max(NVnorm); %how many selective VP neurons
[~,R_2R.RegMaxSelX2(2,1),R_2R.RegMaxSelX2(2,2)]=crosstab(cat(1,ones(max(NNperBin),1),zeros(sum(NAneurons)-max(NNperBin),1),ones(max(NVperBin),1),zeros(sum(VPneurons)-max(NVperBin),1)),NAneurons);


save('R_2R.mat','R_2R');

%divide neurons up by region
NAneurons=strcmp(R_2R.Ninfo(:,3),'NA');
VPneurons=strcmp(R_2R.Ninfo(:,3),'VP');

Xaxis=[-1 2];
Ishow=find(R_2R.Param.Tm>=Xaxis(1) & R_2R.Param.Tm<=Xaxis(2));
time1=R_2R.Param.Tm(Ishow);
Xaxis2=[-2 2];
Ushow=find(R_2R.Param.Tm>=Xaxis2(1) & R_2R.Param.Tm<=Xaxis2(2));
time2=R_2R.Param.Tm(Ushow);

inh=[0.1 0.021154 0.6];
exc=[0.9 0.75 0.205816];

Cue=strcmp('Cue', R_2R.Erefnames);
PE=strcmp('PECue', R_2R.Erefnames);
RD=strcmp('RD', R_2R.Erefnames); %should this be LickR or RD?
%% Plotting heatmaps of all neurons to each event
figure;

%color map
[magma,inferno,plasma,viridis]=colormaps;
colormap(plasma);
c=[-100 2000];ClimE=sign(c).*abs(c).^(1/4);%colormap

for i=1:2
    
    if i==1 Reg=NAneurons; end %1 is NAc
    if i==2 Reg=VPneurons; end %2 is VP
    
    Sel=Reg&R_2R.Ev(Cue).RespDir<2; %gets all neurons

    %sort cue heatmap by magnitude of response
    Neurons=R_2R.Ev(Cue).PSTHz(Sel,Ishow); %get the firing rates of neurons of interest
    TMP=R_2R.Ev(Cue).Meanz(Sel); %find the magnitude of the inhibitions for this event
    TMP(isnan(TMP))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
    [~,SORTimg]=sort(TMP);
    Neurons=Neurons(SORTimg,:); %sort the neurons by magnitude

    %heatmap
    subplot(4,6,[1 7]+(i-1)*12);
    imagesc(time1,[1,sum(Sel,1)],Neurons,ClimE); title('Cue responses');
    xlabel('Seconds post cue');
    ylabel('Individual neurons sorted by response strength');
    hold on;
    plot([0 0],[0 sum(Sel)],':','color','k','linewidth',0.75);

    %sort PE heatmap by magnitude of response
    Neurons=R_2R.Ev(PE).PSTHz(Sel,Ushow); %get the firing rates of neurons of interest
    TMP=R_2R.Ev(PE).Meanz(Sel); %find the magnitude of the inhibitions for this event
    TMP(isnan(TMP))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
    [~,SORTimg]=sort(TMP);
    Neurons=Neurons(SORTimg,:); %sort the neurons by magnitude

    %heatmap
    subplot(4,6,[3 9]+(i-1)*12);
    imagesc(time2,[1,sum(Sel,1)],Neurons,ClimE); title('PE responses');
    xlabel('Seconds post PE');
    hold on;
    plot([0 0],[0 sum(Sel)],':','color','k','linewidth',0.75);

    %sort RD heatmap by magnitude of response
    Neurons=R_2R.Ev(RD).PSTHz(Sel,Ushow); %get the firing rates of neurons of interest
    TMP=R_2R.Ev(RD).Meanz(Sel); %find the magnitude of the inhibitions for this event
    TMP(isnan(TMP))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
    [~,SORTimg]=sort(TMP);
    Neurons=Neurons(SORTimg,:); %sort the neurons by magnitude

    %heatmap
    subplot(4,6,[5 11]+(i-1)*12);
    imagesc(time2,[1,sum(Sel,1)],Neurons,ClimE); title('Reward responses');
    xlabel('Seconds post RD');
    hold on;
    plot([0 0],[0 sum(Sel)],':','color','k','linewidth',0.75);

    %% Plotting neurons that respond to the cue

    %inhibitions
    Sel=Reg&R_2R.Ev(Cue).RespDir<0; %Find neurons that have an inhibitory response to this event

    %average firing rate
    subplot(4,6,2+(i-1)*12);
    psthI=nanmean(R_2R.Ev(Cue).PSTHz(Sel,Ishow),1); %find the average firing rate of the neurons at the time of the event
    semI=nanste(R_2R.Ev(Cue).PSTHz(Sel,Ishow),1); %calculate standard error of the mean
    upI=psthI+semI;
    downI=psthI-semI;
    hold on;
    patch([time1,time1(end:-1:1)],[upI,downI(end:-1:1)],inh,'EdgeColor','none');alpha(0.5);
    plot(time1,psthI,'Color',inh,'linewidth',1); title('Mean firing (z-score)'); %create plot of avg firing rate
    plot([-1 2],[0 0],':','color','k','linewidth',0.75);
    plot([0 0],[-3 9],':','color','k','linewidth',0.75);
    axis([-1 2 min(downI)-0.2 max(upI)+0.2]);
    title(['Cue inhibitions (' num2str(round(100*sum(Sel)/sum(Reg))) '%)'])
    ylabel('Z-score');

    %excitation
    Sel=Reg&R_2R.Ev(Cue).RespDir>0; %Find neurons that have an excitatory response to this event

    %average firing rate
    subplot(4,6,8+(i-1)*12);
    psthE=nanmean(R_2R.Ev(Cue).PSTHz(Sel,Ishow),1);
    semE=nanste(R_2R.Ev(Cue).PSTHz(Sel,Ishow),1); %calculate standard error of the mean
    upE=psthE+semE;
    downE=psthE-semE;
    hold on;
    patch([time1,time1(end:-1:1)],[upE,downE(end:-1:1)],exc,'EdgeColor','none');alpha(0.5);
    plot(time1,psthE,'Color',exc,'linewidth',1);
    plot([-1 2],[0 0],':','color','k','linewidth',0.75);
    plot([0 0],[-3 9],':','color','k','linewidth',0.75);
    axis([-1 2 min(downE)-0.2 max(upE)+0.3]);
    title(['Cue excitations (' num2str(round(100*sum(Sel)/sum(Reg))) '%)'])
    ylabel('Z-score');

    %% Plotting neurons that respond to PE

    %inhibitions
    Sel=Reg&R_2R.Ev(PE).RespDir<0; %Find neurons that have an inhibitory response to this event

    %average firing rate
    subplot(4,6,4+(i-1)*12);
    psthI=nanmean(R_2R.Ev(PE).PSTHz(Sel,Ushow),1); %find the average firing rate of the neurons at the time of the event
    semI=nanste(R_2R.Ev(PE).PSTHz(Sel,Ushow),1); %calculate standard error of the mean
    upI=psthI+semI;
    downI=psthI-semI;
    hold on;
    patch([time2,time2(end:-1:1)],[upI,downI(end:-1:1)],inh,'EdgeColor','none');alpha(0.5);
    plot(time2,psthI,'Color',inh,'linewidth',1); title('Mean firing (z-score)'); %create plot of avg firing rate
    plot([-2 2],[0 0],':','color','k','linewidth',0.75);
    plot([0 0],[-3 9],':','color','k','linewidth',0.75);
    axis([-2 2 min(downI)-0.2 max(upI)+0.2]);
    title(['PE inhibitions (' num2str(round(100*sum(Sel)/sum(Reg))) '%)'])
    ylabel('Z-score');

    %excitation
    Sel=Reg&R_2R.Ev(PE).RespDir>0; %Find neurons that have an excitatory response to this event

    %average firing rate
    subplot(4,6,10+(i-1)*12);
    psthE=nanmean(R_2R.Ev(PE).PSTHz(Sel,Ushow),1);
    semE=nanste(R_2R.Ev(PE).PSTHz(Sel,Ushow),1); %calculate standard error of the mean
    upE=psthE+semE;
    downE=psthE-semE;
    hold on;
    patch([time2,time2(end:-1:1)],[upE,downE(end:-1:1)],exc,'EdgeColor','none');alpha(0.5);
    plot(time2,psthE,'Color',exc,'linewidth',1);
    plot([-2 2],[0 0],':','color','k','linewidth',0.75);
    plot([0 0],[-3 9],':','color','k','linewidth',0.75);
    axis([-2 2 min(downE)-0.2 max(upE)+0.3]);
    title(['PE excitations (' num2str(round(100*sum(Sel)/sum(Reg))) '%)'])
    ylabel('Z-score');

    %% Plotting neurons that respond to reward

    %inhibitions
    Sel=Reg&R_2R.Ev(RD).RespDir<0; %Find neurons that have an inhibitory response to this event

    %average firing rate
    subplot(4,6,6+(i-1)*12);
    psthI=nanmean(R_2R.Ev(RD).PSTHz(Sel,Ushow),1); %find the average firing rate of the neurons at the time of the event
    semI=nanste(R_2R.Ev(RD).PSTHz(Sel,Ushow),1); %calculate standard error of the mean
    upI=psthI+semI;
    downI=psthI-semI;
    hold on;
    patch([time2,time2(end:-1:1)],[upI,downI(end:-1:1)],inh,'EdgeColor','none');alpha(0.5);
    plot(time2,psthI,'Color',inh,'linewidth',1); title('Mean firing (z-score)'); %create plot of avg firing rate
    plot([-1 2],[0 0],':','color','k','linewidth',0.75);
    plot([0 0],[-3 9],':','color','k','linewidth',0.75);
    axis([-1 2 min(downI)-0.2 max(upI)+0.2]);
    title(['RD inhibitions (' num2str(round(100*sum(Sel)/sum(Reg))) '%)'])
    ylabel('Z-score');

    %excitation
    Sel=Reg&R_2R.Ev(RD).RespDir>0; %Find neurons that have an excitatory response to this event

    %average firing rate
    subplot(4,6,12+(i-1)*12);
    psthE=nanmean(R_2R.Ev(RD).PSTHz(Sel,Ushow),1);
    semE=nanste(R_2R.Ev(RD).PSTHz(Sel,Ushow),1); %calculate standard error of the mean
    upE=psthE+semE;
    downE=psthE-semE;
    hold on;
    patch([time2,time2(end:-1:1)],[upE,downE(end:-1:1)],exc,'EdgeColor','none');alpha(0.5);
    plot(time2,psthE,'Color',exc,'linewidth',1);
    plot([-1 2],[0 0],':','color','k','linewidth',0.75);
    plot([0 0],[-3 9],':','color','k','linewidth',0.75);
    axis([-1 2 min(downE)-0.2 max(upE)+0.3]);
    title(['RD excitations (' num2str(round(100*sum(Sel)/sum(Reg))) '%)'])
    ylabel('Z-score');
end

%% plotting firing rate on sucrose and maltodextrin trials
figure;

%color map
[magma,inferno,plasma,viridis]=colormaps;
colormap(plasma);
c=[-100 2000];ClimE=sign(c).*abs(c).^(1/4);%colormap

Xaxis=[-1 4];
Ishow=find(R_2R.Param.Tm>=Xaxis(1) & R_2R.Param.Tm<=Xaxis(2));
time1=R_2R.Param.Tm(Ishow);
Xaxis2=[-1 4];
Ushow=find(R_2R.Param.Tm>=Xaxis2(1) & R_2R.Param.Tm<=Xaxis2(2));
time2=R_2R.Param.Tm(Ushow);

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];
exc=[0 113/255 188/255];
inh=[240/255 0 50/255];

RD=strcmp('RD', R_2R.Erefnames);
RD1=strcmp('RD1', R_2R.Erefnames);
RD2=strcmp('RD2', R_2R.Erefnames);

for i=1:2
    
    if i==1 Reg=NAneurons; end %1 is NAc
    if i==2 Reg=VPneurons; end %2 is VP

    %% Plotting heatmaps of sucrose and maltodextrin response
    Sel=Reg&R_2R.Ev(RD).RespDir<2; %gets all neurons

    %sort RD1 heatmap by magnitude of response
    Neurons=R_2R.Ev(RD1).PSTHz(Sel,Ushow); %get the firing rates of neurons of interest
    TMP=R_2R.Ev(RD).Meanz(Sel); %find the magnitude of the inhibitions for this event
    TMP(isnan(TMP))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
    [~,SORTimg]=sort(TMP);
    Neurons=Neurons(SORTimg,:); %sort the neurons by magnitude

    %heatmap
    subplot(2,40,[1:9]+(i-1)*40);
    imagesc(time2,[1,sum(Sel,1)],Neurons,ClimE); title('Sucrose trials');
    ylabel('Individual neurons sorted by RD response');
    xlabel('Seconds from RD');
    hold on;
    plot([0 0],[0 sum(Sel)],':','color','k','linewidth',0.75);

    %sort RD2 heatmap the same way
    Neurons=R_2R.Ev(RD2).PSTHz(Sel,Ushow); %get the firing rates of neurons of interest
    [~,SORTimg]=sort(TMP);
    Neurons=Neurons(SORTimg,:); %sort the neurons by magnitude
    xlabel('Seconds from RD');


    %heatmap
    subplot(2,40,[13:21]+(i-1)*40);
    imagesc(time2,[1,sum(Sel,1)],Neurons,ClimE); title('Maltodextrin trials');
    xlabel('Seconds from RD');
    hold on;
    plot([0 0],[0 sum(Sel)],':','color','k','linewidth',0.75);

    %% Plotting average firing rate of reward excitations
    Sel=Reg&(R_2R.Ev(RD).RespDir>0); %| R_2R.Ev(LickR1).RespDir>0 | R_2R.Ev(LickR2).RespDir>0); %Find neurons that have an inhibitory response to this event

    %average firing rate to sucrose
    psth1=nanmean(R_2R.Ev(RD1).PSTHz(Sel,Ishow),1);
    sem1=nanste(R_2R.Ev(RD1).PSTHz(Sel,Ishow),1); %calculate standard error of the mean
    up1=psth1+sem1;
    down1=psth1-sem1;

    %average firing rate to maltodextrin
    psth2=nanmean(R_2R.Ev(RD2).PSTHz(Sel,Ishow),1);
    sem2=nanste(R_2R.Ev(RD2).PSTHz(Sel,Ishow),1); %calculate standard error of the mean
    up2=psth2+sem2;
    down2=psth2-sem2;

    %plotting
    subplot(4,5,[9 10]+(i-1)*10);
    hold on;
    plot(time1,psth1,'Color',sucrose,'linewidth',1);
    plot(time1,psth2,'Color',maltodextrin,'linewidth',1);
    patch([time1,time1(end:-1:1)],[up1,down1(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
    patch([time1,time1(end:-1:1)],[up2,down2(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
    plot([-1 4],[0 0],':','color','k','linewidth',0.75);
    plot([0 0],[-3 9],':','color','k','linewidth',0.75);
    axis([-1 4 min(down1)-0.2 max(up1)+0.3]);
    title(['Reward-excited (n=' num2str(sum(Sel)) ')'])
    ylabel('Mean firing (z-score)');
    xlabel('Seconds from RD');
    legend('Suc trials','Mal trials');

    % %onset
    % subplot(2,4,4);
    % cumul(R_2R.Ev(38).Onsets(Sel,1),[R_2R.Param.RespWin(1,1):0.04:R_2R.Param.RespWin(38,2)],2,maltodextrin);


    %% Plotting average firing rate of reward inhibitions
    Sel=Reg&(R_2R.Ev(RD).RespDir<0); %| R_2R.Ev(LickR1).RespDir<0 | R_2R.Ev(LickR2).RespDir<0); %Find neurons that have an inhibitory response to this event

    %avg firing rate to sucrose
    psth1=nanmean(R_2R.Ev(RD1).PSTHz(Sel,Ushow),1);
    sem1=nanste(R_2R.Ev(RD1).PSTHz(Sel,Ushow),1); %calculate standard error of the mean
    up1=psth1+sem1;
    down1=psth1-sem1;

    %avg firing rate to maltodextrin
    psth2=nanmean(R_2R.Ev(RD2).PSTHz(Sel,Ushow),1);
    sem2=nanste(R_2R.Ev(RD2).PSTHz(Sel,Ushow),1); %calculate standard error of the mean
    up2=psth2+sem2;
    down2=psth2-sem2;

    %plotting
    subplot(4,5,[4 5]+(i-1)*10);
    hold on;
    plot(time2,psth1,'Color',sucrose,'linewidth',1);
    plot(time2,psth2,'Color',maltodextrin,'linewidth',1);
    patch([time2,time2(end:-1:1)],[up2,down2(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
    patch([time2,time2(end:-1:1)],[up1,down1(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
    plot([-1 4],[0 0],':','color','k','linewidth',0.75);
    plot([0 0],[-3 9],':','color','k','linewidth',0.75);
    axis([-1 4 min(down1)-0.2 max(up1)+0.2]);
    title(['Reward-inhibited (n=' num2str(sum(Sel)) ')'])
    ylabel('Mean firing (z-score)');
    xlabel('Seconds from RD');
    legend('Suc trials','Mal trials','location','southeast');
end