%% Code for final figures for pavlovian versus instrumental comparison
if ~exist('RvppaALL'),load('RvppaALL.mat'); end
if ~exist('RvpppaALL'),load('RvpppaALL.mat'); end
DarkPurple=[123/255 50/255 148/255];
LightPurple=[194/255 165/255 207/255];
LightGreen=[166/255 219/255 160/255];
DarkGreen=[77/255 172/255 38/255];
DarkPink=[208/255 28/255 139/255];
MyBlue=[0.4940    0.1840    0.5560];
MyOrange=[0.9290    0.6940    0.1250];

BSIZE=0.01;
Dura=[-0.5 1]; Tm=Dura(1):BSIZE:Dura(2);
Baseline=[-2 0]; Tbase=Baseline(1):BSIZE:Baseline(2); %used to calculate Z-scores
Tbin=-1:0.005:1; %window used to determine the optimal binsize
PStat=0.05;
prewin=[Dura(1) 0];
postwin=[.1 .3];
postwin2=[Dura(1):0.05:Dura(2)]; %bounds should match Dura
Slopebounds=[-.5:.05:.5];
R2Bounds=[0:0.05:0.5];
PvalBounds=[0:0.01:0.5];
Struct=10;


%% Figure 3. VP neurons encode the value of discriminative stimuli for alcohol availability but fail to encode the value of Pavlovian alcohol cues. 

BrainRegion=1; %10 for Core, 100 for dms and 1000 for shell
LeftEv=1;          RightEv=2;
ExcInh=1;        % 1 for excitations, 2 for inhibitions -- Used to find Exc or Inh onsets in R.Ev(.).Onsets

Pstat=0.05;
stepsize=0.3;      % Used for the scatter plot projections

%For COLOR-CODED PSTH only (not Average PSTH)
Norm=0;          % 0 = no norm - max accross all the selected neurons
% 1 each neuron is normalized wh its own max
% 2 for baseline substraction

if LeftEv==1, Xaxis=[-0.25 .5]; elseif LeftEv==3, Xaxis=[-0.25 .5]; elseif LeftEv==5; Xaxis=[-3 3]; elseif LeftEv==7; Xaxis=[-1 2]; elseif LeftEv==9; Xaxis=[-0.25 1]; end
if ExcInh==1, Yaxis=[-5 10]; elseif ExcInh==2, Yaxis=[-5 5]; end

%----------- COLORS ---------
myblue=[0 113/255 188/255];
mypurple=[200/255 0 255/255];
mydarkblue=[0 0/255 255/255];
TickSize=[0.015 0.02];

DarkPurple=[123/255 50/255 148/255];
LightPurple=[194/255 165/255 207/255];
LightGreen=[166/255 219/255 160/255];
DarkGreen=[77/255 172/255 38/255];
DarkPink=[208/255 28/255 139/255];



Ishow=find(RvppaALL.Param.Tm>=Xaxis(1) & RvppaALL.Param.Tm<=Xaxis(2));
time=RvppaALL.Param.Tm(Ishow);
if ExcInh==1, A=1; elseif ExcInh==2; A=-1;end
figure;
for i=1:2 %plots 3 different conditions (eg LP+No0 / LP+No+ / LP0No+)
    
    selectionInst=RvpppaALL.Structure==10 & ~isnan(RvpppaALL.Ev(LeftEv).Meanz) & ~isnan(RvpppaALL.Ev(RightEv).Meanz) & RvpppaALL.CSPlusRatio>=.5 & (RvpppaALL.CSPlusRatio./(RvpppaALL.CSPlusRatio+RvpppaALL.CSMinusRatio))>=.7;%) & RvpppaALL.Ev(1).RespDir==0;% & CvpppaALL.STATSpost>=0.05;
    selectionPav=RvppaALL.Structure==10 & ~isnan(RvppaALL.Ev(LeftEv).Meanz) & ~isnan(RvppaALL.Ev(RightEv).Meanz) & RvppaALL.CSPlusRatio>=.5 & (RvppaALL.CSPlusRatio./(RvppaALL.CSPlusRatio+RvppaALL.CSMinusRatio))>=.7;% & RvppaALL.Ev(1).RespDir==0;% & CvppaALL.STATSpost>=0.05;
    
    %------------- SORTING -------------
    SortEv=LeftEv; NonSortEv=RightEv;
    if i==1;
        selection=selectionInst;
        TMP=-RvpppaALL.Ev(SortEv).Meanz(selection,ExcInh);
        %TMP=RvpppaALL.Ev(SortEv).RespDir(selection,ExcInh);
    else;
        selection=selectionPav;
        TMP=-RvppaALL.Ev(SortEv).Meanz(selection,ExcInh);
        %TMP=RvppaALL.Ev(SortEv).RespDir(selection,ExcInh);
    end
    
    TMP(isnan(TMP))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
    [~,SORTimg]=sort(TMP);
    
    %NNids are the indexes for the 3 different selections
    if i==1; NNid1=cell2mat(RvpppaALL.Ninfo(selection,3)); NNid1=NNid1(SORTimg); end
    if i==2; NNid2=cell2mat(RvppaALL.Ninfo(selection,3)); NNid2=NNid2(SORTimg); end
    
    if i==1,
        selection=selectionInst; SortEv=LeftEv; NonSortEv=RightEv;
        % ------------- NORMALIZATION -------------
        if Norm==0
            imgLeft=RvpppaALL.Ev(LeftEv).PSTHz(selection,Ishow);
            imgRight=RvpppaALL.Ev(RightEv).PSTHz(selection,Ishow);
            E=[-200 400];Clim=sign(E).*abs(E).^(1/4);%ColorMap
        elseif Norm~=0
            imgLeft=normalize(RvpppaALL.Ev(LeftEv).PSTHz(selection,Ishow),RvpppaALL.Ev(LeftEv).PSTHz(selection,(RvpppaALL.Param.Tm<0)),Norm); %normalize to max response
            imgRight=normalize(RvpppaALL.Ev(RightEv).PSTHz(selection,Ishow),RvpppaALL.Ev(LeftEv).PSTHz(selection,(RvpppaALL.Param.Tm<0)),Norm); %normalize to max response of the left event
            if Norm==1, Clim=[0 1]; elseif Norm==2, Clim=[-5 10]; end
        end
        
        imgLeft=imgLeft(SORTimg,:);
        imgRight=imgRight(SORTimg,:);
        
    else if i==2,
            selection=selectionPav; SortEv=LeftEv; NonSortEv=RightEv;
            if Norm==0
                imgLeft=RvppaALL.Ev(LeftEv).PSTHz(selection,Ishow);
                imgRight=RvppaALL.Ev(RightEv).PSTHz(selection,Ishow);
                E=[-200 400];Clim=sign(E).*abs(E).^(1/4);%ColorMap
            elseif Norm~=0
                imgLeft=normalize(RvppaALL.Ev(LeftEv).PSTHz(selection,Ishow),RvppaALL.Ev(LeftEv).PSTHz(selection,(RvppaALL.Param.Tm<0)),Norm); %normalize to max response
                imgRight=normalize(RvppaALL.Ev(RightEv).PSTHz(selection,Ishow),RvppaALL.Ev(LeftEv).PSTHz(selection,(RvppaALL.Param.Tm<0)),Norm); %normalize to max response of the left event
                if Norm==1, Clim=[0 1]; elseif Norm==2, Clim=[-5 10]; end
            end
            
            imgLeft=imgLeft(SORTimg,:);
            imgRight=imgRight(SORTimg,:);
            
        end
    end
    
    %---------------- HEAT MAPS ----------------------
    subplot(3,4,i+((i-1)*3));
    imagesc(time,[1 size(imgLeft,1)],imgLeft,Clim); %colorbar;axis tight;
    colormap(parula);
    ylabel('Neuron #');
    xlabel('Time from Cue (sec)');
    if i==1;
        title('DS');
    elseif i==2;
        title('CS+');
    end
    
    subplot(3,4,i+1+((i-1)*3));
    imagesc(time,[1 size(imgRight,1)],imgRight,Clim) ;%colorbar;axis tight;
    colormap(parula);
    %title(Erefnames(RightEv));
    ylabel('Neuron #');
    xlabel('Time from Cue (sec)');
    if i==1;
        title('NS');
    elseif i==2;
        title('CS-');
    end
    
    subplot(3,2,i*2);
    hold on;
    if i==1 %PPS
        if     LeftEv==1, Stat=RvpppaALL.CSStat;
        elseif LeftEv==3, Stat=RvpppaALL.CSPlusResponseStat;
        end
        selection1=selectionInst & Stat>=0.05;
        selection2=selectionInst & Stat<0.05 & RvpppaALL.Ev(LeftEv).Meanz>RvpppaALL.Ev(RightEv).Meanz;
        selection3=selectionInst & Stat<0.05 & RvpppaALL.Ev(LeftEv).Meanz<RvpppaALL.Ev(RightEv).Meanz;
        scatter(RvpppaALL.Ev(LeftEv).Meanz(selection1),RvpppaALL.Ev(RightEv).Meanz(selection1),30, 'k');%,'filled');
        scatter(RvpppaALL.Ev(LeftEv).Meanz(selection2),RvpppaALL.Ev(RightEv).Meanz(selection2),30, DarkPink);
        scatter(RvpppaALL.Ev(LeftEv).Meanz(selection3),RvpppaALL.Ev(RightEv).Meanz(selection3),30, DarkGreen);
        %selection4= R.Ev(RightEv).RespDir==1 & R.Ev(6).RespDir==1;
        %scatter(R.Ev(LeftEv).Meanz(selection4),R.Ev(RightEv).Meanz(selection4),30,'g','filled');
        xlabel('DS Response Magnitude (z-score)');
        ylabel('NS Response magnitude (z-score)');
        legend('DS=NS','DS>NS','DS<NS','Location','northwest','AutoUpdate','off');
    else %PS
        if     LeftEv==1, Stat=RvppaALL.CSStat;
        elseif LeftEv==3, Stat=RvppaALL.CSResponseStat;
        end
        selection1=selectionPav & Stat>=0.05;
        selection2=selectionPav & Stat<0.05 & RvppaALL.Ev(LeftEv).Meanz>RvppaALL.Ev(RightEv).Meanz;
        selection3=selectionPav & Stat<0.05 & RvppaALL.Ev(LeftEv).Meanz<RvppaALL.Ev(RightEv).Meanz;
        scatter(RvppaALL.Ev(LeftEv).Meanz(selection1),RvppaALL.Ev(RightEv).Meanz(selection1),30, 'k');%,'filled');
        scatter(RvppaALL.Ev(LeftEv).Meanz(selection2),RvppaALL.Ev(RightEv).Meanz(selection2),30, DarkPink);
        scatter(RvppaALL.Ev(LeftEv).Meanz(selection3),RvppaALL.Ev(RightEv).Meanz(selection3),30, DarkGreen);
        Min=floor(min([RvppaALL.Ev(LeftEv).Meanz(selection2);RvppaALL.Ev(RightEv).Meanz(selection2)]));
        Max=ceil(max([RvppaALL.Ev(LeftEv).Meanz(selection2);RvppaALL.Ev(RightEv).Meanz(selection2)]));
        xlabel('CS+ Response Magnitude (z-score)');
        ylabel('CS- Response magnitude (z-score)');
        legend('CS+=CS-','CS+>CS-','CS+<CS-','Location','northwest','AutoUpdate','off');
    end
    Min=floor(min([RvpppaALL.Ev(LeftEv).Meanz(RvpppaALL.Structure==10);RvpppaALL.Ev(RightEv).Meanz(RvpppaALL.Structure==10)]));
    Max=ceil(max([RvpppaALL.Ev(LeftEv).Meanz(RvpppaALL.Structure==10);RvpppaALL.Ev(RightEv).Meanz(RvpppaALL.Structure==10)]));
    line([Min Max],[Min Max],'Color', 'k');
    line([Min Max], [0 0],'Color', 'k');
    line([0 0], [Min Max],'Color', 'k');
    
    
    
    subplot(10,6,58); %Histograms for scatterplots
    if     LeftEv==1, Stat=RvpppaALL.CSStat;
    elseif LeftEv==3, Stat=RvpppaALL.CSPlusResponseStat;
    end
    axis([-10 10 0 .25]);
    edges = [-10:1:9];
    selection1=selectionInst & Stat>=0.05;
    selection2=selectionInst & Stat<0.05 & RvpppaALL.Ev(LeftEv).Meanz>RvpppaALL.Ev(RightEv).Meanz;
    selection3=selectionInst & Stat<0.05 & RvpppaALL.Ev(LeftEv).Meanz<RvpppaALL.Ev(RightEv).Meanz;
    Delta1= RvpppaALL.Ev(LeftEv).Meanz(selection1)-RvpppaALL.Ev(RightEv).Meanz(selection1);
    Delta2= RvpppaALL.Ev(LeftEv).Meanz(selection2)-RvpppaALL.Ev(RightEv).Meanz(selection2);
    Delta3= RvpppaALL.Ev(LeftEv).Meanz(selection3)-RvpppaALL.Ev(RightEv).Meanz(selection3);
    Delta1Hist=hist(Delta1,edges)/314;
    Delta2Hist=hist(Delta2,edges)/314;
    Delta3Hist=hist(Delta3,edges)/314;
    edges = [-10:1:10];
    histogram('FaceColor','k', 'BinEdges',edges,'BinCounts',Delta1Hist);
    hold on;
    histogram('FaceColor',DarkPink, 'BinEdges',edges,'BinCounts',Delta2Hist);
    histogram('FaceColor',DarkGreen, 'BinEdges',edges,'BinCounts',Delta3Hist);
    
    
    subplot(10,6,59); %Histograms for scatterplots
    if     LeftEv==1, Stat=RvppaALL.CSStat;
    elseif LeftEv==3, Stat=RvppaALL.CSrespStat;
    end
    axis([-10 10 0 .25]);
    edges = [-10:1:9];
    selection1=selectionPav & Stat>=0.05;
    selection2=selectionPav & Stat<0.05 & RvppaALL.Ev(LeftEv).Meanz>RvppaALL.Ev(RightEv).Meanz;
    selection3=selectionPav & Stat<0.05 & RvppaALL.Ev(LeftEv).Meanz<RvppaALL.Ev(RightEv).Meanz;
    Delta1= RvppaALL.Ev(LeftEv).Meanz(selection1)-RvppaALL.Ev(RightEv).Meanz(selection1);
    Delta2= RvppaALL.Ev(LeftEv).Meanz(selection2)-RvppaALL.Ev(RightEv).Meanz(selection2);
    Delta3= RvppaALL.Ev(LeftEv).Meanz(selection3)-RvppaALL.Ev(RightEv).Meanz(selection3);
    Delta1Hist=hist(Delta1,edges)/392;
    Delta2Hist=hist(Delta2,edges)/392;
    Delta3Hist=hist(Delta3,edges)/392;
    edges = [-10:1:10];
    histogram('FaceColor','k', 'BinEdges',edges,'BinCounts',Delta1Hist);
    hold on;
    histogram('FaceColor',DarkPink, 'BinEdges',edges,'BinCounts',Delta2Hist);
    histogram('FaceColor',DarkGreen, 'BinEdges',edges,'BinCounts',Delta3Hist);
    
    selection1=selectionInst & ~isnan(RvpppaALL.Ev(LeftEv).Meanz) & ~isnan(RvpppaALL.Ev(RightEv).Meanz)& RvpppaALL.Ev(1).RespDir==1;% & CvpppaALL.STATSpost>=0.05;
    selection2=selectionPav & ~isnan(RvppaALL.Ev(LeftEv).Meanz) & ~isnan(RvppaALL.Ev(RightEv).Meanz)& RvppaALL.Ev(1).RespDir==1;% & CvppaALL.STATSpost>=0.05;
    
    %subplot of DS vs NS excited neurons
    if i==1
        selection=selection1;
        subplot(6,4,i+16);
        psthLeft=nanmean(RvpppaALL.Ev(LeftEv).PSTHz(selection,Ishow),1);
        semLeft=nanste(RvpppaALL.Ev(LeftEv).PSTHz(selection,Ishow),1);
        upLeft=psthLeft+semLeft;
        downLeft=psthLeft-semLeft;
        hold on;
        
        plot(time,psthLeft,'Color',MyBlue,'linewidth',1.5);
        
        psthRight=nanmean(RvpppaALL.Ev(RightEv).PSTHz(selection,Ishow),1);
        semRight=nanste(RvpppaALL.Ev(RightEv).PSTHz(selection,Ishow),1);
        upRight=psthRight+semRight;
        downRight=psthRight-semRight;
        
        plot(time,psthRight,'k','linewidth',1.5);
        
        patch([time,time(end:-1:1)],[upLeft,downLeft(end:-1:1)],MyBlue,'EdgeColor','none');alpha(0.5);
        patch([time,time(end:-1:1)],[upRight,downRight(end:-1:1)],'k','EdgeColor','none');alpha(0.5);
        plot([0 0], [-10 20],'LineStyle',':','Color','k');
        plot([-10 20], [0 0],'LineStyle',':','Color','k');
        if LeftEv==3
            plot([1 1], [-10 20],'LineStyle',':','Color','k');
            plot([2.5 2.5], [-10 20],'LineStyle',':','Color','k');
        elseif LeftEv==5
            plot([-2 -2], [-10 20],'LineStyle',':','Color','k');
            plot([.5 .5], [-10 20],'LineStyle',':','Color','k');
        elseif LeftEv==7
            plot([1 1], [-10 20],'LineStyle',':','Color','k');
            plot([5 5], [-10 20],'LineStyle',':','Color','k');
        end
        axis([Xaxis,Yaxis]);
        %xlabel('Time from Cue (sec)');
        ylabel('Average Response (z-score)');
        legend('DS','NS','Location','northeast');
        
    else if i==2
            selection=selection2;
            subplot(6,4,i+19);
            psthLeft=nanmean(RvppaALL.Ev(LeftEv).PSTHz(selection,Ishow),1);
            semLeft=nanste(RvppaALL.Ev(LeftEv).PSTHz(selection,Ishow),1);
            upLeft=psthLeft+semLeft;
            downLeft=psthLeft-semLeft;
            hold on;
            
            plot(time,psthLeft,'Color',MyOrange,'linewidth',1.5);
            
            psthRight=nanmean(RvppaALL.Ev(RightEv).PSTHz(selection,Ishow),1);
            semRight=nanste(RvppaALL.Ev(RightEv).PSTHz(selection,Ishow),1);
            upRight=psthRight+semRight;
            downRight=psthRight-semRight;
            
            plot(time,psthRight,'k','linewidth',1.5);
            
            patch([time,time(end:-1:1)],[upLeft,downLeft(end:-1:1)],MyOrange,'EdgeColor','none');alpha(0.5);
            patch([time,time(end:-1:1)],[upRight,downRight(end:-1:1)],'k','EdgeColor','none');alpha(0.5);
            plot([0 0], [-10 20],'LineStyle',':','Color','k');
            plot([-10 20], [0 0],'LineStyle',':','Color','k');
            if LeftEv==3
                plot([1 1], [-10 20],'LineStyle',':','Color','k');
                plot([2.5 2.5], [-10 20],'LineStyle',':','Color','k');
            elseif LeftEv==5
                plot([-2 -2], [-10 20],'LineStyle',':','Color','k');
                plot([.5 .5], [-10 20],'LineStyle',':','Color','k');
            elseif LeftEv==7
                plot([1 1], [-10 20],'LineStyle',':','Color','k');
                plot([5 5], [-10 20],'LineStyle',':','Color','k');
            end
            axis([Xaxis,Yaxis]);
            
            xlabel('Time from Cue (sec)');
            ylabel('Average Response (z-score)');
            legend('CS+','CS-','Location','northeast');
        end
    end
        
        
        selection1=selectionInst & ~isnan(RvpppaALL.Ev(LeftEv).Meanz) & ~isnan(RvpppaALL.Ev(RightEv).Meanz)& RvpppaALL.Ev(1).RespDir==-1;
        selection2=selectionPav & ~isnan(RvppaALL.Ev(LeftEv).Meanz) & ~isnan(RvppaALL.Ev(RightEv).Meanz)& RvppaALL.Ev(1).RespDir==-1;
        
        if i==1
            selection=selection1;
            subplot(6,4,i+17);
            psthLeft=nanmean(RvpppaALL.Ev(LeftEv).PSTHz(selection,Ishow),1);
            semLeft=nanste(RvpppaALL.Ev(LeftEv).PSTHz(selection,Ishow),1);
            upLeft=psthLeft+semLeft;
            downLeft=psthLeft-semLeft;
            hold on;
            
            plot(time,psthLeft,'Color',MyBlue,'linewidth',1.5);
            
            psthRight=nanmean(RvpppaALL.Ev(RightEv).PSTHz(selection,Ishow),1);
            semRight=nanste(RvpppaALL.Ev(RightEv).PSTHz(selection,Ishow),1);
            upRight=psthRight+semRight;
            downRight=psthRight-semRight;
            
            plot(time,psthRight,'k','linewidth',1.5);
            
            patch([time,time(end:-1:1)],[upLeft,downLeft(end:-1:1)],MyBlue,'EdgeColor','none');alpha(0.5);
            patch([time,time(end:-1:1)],[upRight,downRight(end:-1:1)],'k','EdgeColor','none');alpha(0.5);
            plot([0 0], [-10 20],'LineStyle',':','Color','k');
            plot([-10 20], [0 0],'LineStyle',':','Color','k');
            if LeftEv==3
                plot([1 1], [-10 20],'LineStyle',':','Color','k');
                plot([2.5 2.5], [-10 20],'LineStyle',':','Color','k');
            elseif LeftEv==5
                plot([-2 -2], [-10 20],'LineStyle',':','Color','k');
                plot([.5 .5], [-10 20],'LineStyle',':','Color','k');
            elseif LeftEv==7
                plot([1 1], [-10 20],'LineStyle',':','Color','k');
                plot([5 5], [-10 20],'LineStyle',':','Color','k');
            end
            axis([Xaxis,Yaxis]);
            legend('DS','NS','Location','northwest');
            
            %xlabel('Time from Cue (sec)');
            %ylabel('Average Response (z-score)');
            
        else if i==2
                selection=selection2;
                subplot(6,4,i+20);
                psthLeft=nanmean(RvppaALL.Ev(LeftEv).PSTHz(selection,Ishow),1);
                semLeft=nanste(RvppaALL.Ev(LeftEv).PSTHz(selection,Ishow),1);
                upLeft=psthLeft+semLeft;
                downLeft=psthLeft-semLeft;
                hold on;
                
                plot(time,psthLeft,'Color',MyOrange,'linewidth',1.5);
                
                psthRight=nanmean(RvppaALL.Ev(RightEv).PSTHz(selection,Ishow),1);
                semRight=nanste(RvppaALL.Ev(RightEv).PSTHz(selection,Ishow),1);
                upRight=psthRight+semRight;
                downRight=psthRight-semRight;
                
                plot(time,psthRight,'k','linewidth',1.5);
                
                patch([time,time(end:-1:1)],[upLeft,downLeft(end:-1:1)],MyOrange,'EdgeColor','none');alpha(0.5);
                patch([time,time(end:-1:1)],[upRight,downRight(end:-1:1)],'k','EdgeColor','none');alpha(0.5);
                plot([0 0], [-10 20],'LineStyle',':','Color','k');
                plot([-10 20], [0 0],'LineStyle',':','Color','k');
                if LeftEv==3
                    plot([1 1], [-10 20],'LineStyle',':','Color','k');
                    plot([2.5 2.5], [-10 20],'LineStyle',':','Color','k');
                elseif LeftEv==5
                    plot([-2 -2], [-10 20],'LineStyle',':','Color','k');
                    plot([.5 .5], [-10 20],'LineStyle',':','Color','k');
                elseif LeftEv==7
                    plot([1 1], [-10 20],'LineStyle',':','Color','k');
                    plot([5 5], [-10 20],'LineStyle',':','Color','k');
                end
                axis([Xaxis,Yaxis]);
                
                xlabel('Time from Cue (sec)');
                %ylabel('Average Response (z-score)');
                legend('CS+','CS-','Location','northwest');
            end
        end
    
end

subplot(9,2,[14 18]);
%Plot of auROC descriptions for CS versus DS on response versus no response
%(percent of neurons!)
sel1=RvpppaALL.Structure~=0;
sel2=RvppaALL.Structure~=0;
axis([0 1 0 .25]);
%histogram(-RvpppaALL.CSauROC(RvpppaALL.Structure==10),[0:0.05:1],'Normalization','probability');
cdfplot(-RvpppaALL.CSauROC(RvpppaALL.Structure==10 & RvpppaALL.CSPlusRatio>=.5 & (RvpppaALL.CSPlusRatio./(RvpppaALL.CSPlusRatio+RvpppaALL.CSMinusRatio))>=.7));
hold on;
cdfplot(-RvppaALL.CSauROC(RvppaALL.Structure==10 & RvppaALL.CSPlusRatio>=.5 & (RvppaALL.CSPlusRatio./(RvppaALL.CSPlusRatio+RvppaALL.CSMinusRatio))>=.7));
%histogram(-RvppaALL.CSauROC(RvppaALL.Structure==10),[0:0.05:1],'Normalization','probability');
MEANsel1=nanmedian(-RvpppaALL.CSauROC(RvpppaALL.Structure==10 & RvpppaALL.CSPlusRatio>=.5 & (RvpppaALL.CSPlusRatio./(RvpppaALL.CSPlusRatio+RvpppaALL.CSMinusRatio))>=.7));
MEANsel2=nanmedian(-RvppaALL.CSauROC(RvppaALL.Structure==10 & RvppaALL.CSPlusRatio>=.5 & (RvppaALL.CSPlusRatio./(RvppaALL.CSPlusRatio+RvppaALL.CSMinusRatio))>=.7));
plot([MEANsel1 MEANsel1], [0 1],'Color',MyBlue);
plot([MEANsel2 MEANsel2], [0 1],'Color',MyOrange);
plot([.5 .5], [0 1],'Color','k');
legend('DS','CS+','Location','northwest');
xlabel('Cue Identity auROC');
ylabel('Proportion of Neurons');



%% Figure 2. VP neurons respond to cues and reward seeking in both tasks.

figure;
%DS
subplot(4,5,11);
exc=sum((RvpppaALL.Ev(1).RFLAG(selectionInst)==1 | RvpppaALL.Ev(1).RespDir(selectionInst)==1) & RvpppaALL.Ev(1).RFLAG(selectionInst)~=-1 & RvpppaALL.Ev(1).RespDir(selectionInst)~=-1);
inh=sum(RvpppaALL.Ev(1).RespDir(selectionInst)==-1 | RvpppaALL.Ev(1).RFLAG(selectionInst)==-1);
no=sum(RvpppaALL.Ev(1).RFLAG(selectionInst)==0 & RvpppaALL.Ev(1).RespDir(selectionInst)==0);
response=[exc inh no];
labels = {'Excited','Inhibited','No Response'};
pie(response);
colormap([252/255 206/255 46/255;      %// excitation yellow
    50/255 67/255 186/255;      %// inhibition blue
    46/255 183/255 164/255]);      %// no response teal
title('DS');

%CS+
subplot(4,5,12);
exc=sum((RvppaALL.Ev(1).RFLAG(selectionPav)==1 | RvppaALL.Ev(1).RespDir(selectionPav)==1) & RvppaALL.Ev(1).RFLAG(selectionPav)~=-1 & RvppaALL.Ev(1).RespDir(selectionPav)~=-1);
inh=sum(RvppaALL.Ev(1).RespDir(selectionPav)==-1 | RvppaALL.Ev(1).RFLAG(selectionPav)==-1);
no=sum(RvppaALL.Ev(1).RFLAG(selectionPav)==0 & RvppaALL.Ev(1).RespDir(selectionPav)==0);
response=[exc inh no];
pie(response);
title('CS+');

%DS port entry
subplot(4,5,13);
exc=sum((RvpppaALL.Ev(8).RFLAG(selectionInst)==1 | RvpppaALL.Ev(8).RespDir(selectionInst)==1) & RvpppaALL.Ev(5).RFLAG(selectionInst)~=-1 & RvpppaALL.Ev(5).RespDir(selectionInst)~=-1);
inh=sum(RvpppaALL.Ev(8).RespDir(selectionInst)==-1 | RvpppaALL.Ev(8).RFLAG(selectionInst)==-1);
no=sum(RvpppaALL.Ev(8).RFLAG(selectionInst)==0 & RvpppaALL.Ev(8).RespDir(selectionInst)==0);
response=[exc inh no];
labels = {'Excited','Inhibited','No Response'};
pie(response);
title('post-DS Port Entry');

%CS+ port entry
subplot(4,5,14);
exc=sum((RvppaALL.Ev(8).RFLAG(selectionPav)==1 | RvppaALL.Ev(8).RespDir(selectionPav)==1) & RvppaALL.Ev(5).RFLAG(selectionPav)~=-1 & RvppaALL.Ev(5).RespDir(selectionPav)~=-1);
inh=sum(RvppaALL.Ev(8).RespDir(selectionPav)==-1 | RvppaALL.Ev(8).RFLAG(selectionPav)==-1);
no=sum(RvppaALL.Ev(8).RFLAG(selectionPav)==0 & RvppaALL.Ev(8).RespDir(selectionPav)==0);
response=[exc inh no];
h=pie(response);
title('post-CS+ Port Entry');

%CS+ reward
subplot(4,5,15);
exc=sum((RvppaALL.Ev(9).RFLAG(selectionPav)==1 | RvppaALL.Ev(9).RespDir(selectionPav)==1) & RvppaALL.Ev(7).RFLAG(selectionPav)~=-1 & RvppaALL.Ev(7).RespDir(selectionPav)~=-1);
inh=sum(RvppaALL.Ev(9).RespDir(selectionPav)==-1 | RvppaALL.Ev(9).RFLAG(selectionPav)==-1);
no=sum(RvppaALL.Ev(9).RFLAG(selectionPav)==0 & RvppaALL.Ev(9).RespDir(selectionPav)==0);
response=[exc inh no];
h=pie(response);
title('post-CS+ Reward');