%% Code for final figures for pavlovian versus instrumental comparison
clear all;
if ~exist('RvppasALL'),load('RvppasALL.mat'); end
if ~exist('RvpppasALL'),load('RvpppasALL.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];
DS=[51/255 160/255 44/255];
Pav=[180/255 41/255 230/255];

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 1 VP neurons encode cue value in both tasks

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(RvppasALL.Param.Tm>=Xaxis(1) & RvppasALL.Param.Tm<=Xaxis(2));
time=RvppasALL.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=RvpppasALL.Structure==10 & ~isnan(RvpppasALL.Ev(LeftEv).Meanz) & ~isnan(RvpppasALL.Ev(RightEv).Meanz) & RvpppasALL.CSPlusRatio>=.5 & (RvpppasALL.CSPlusRatio./(RvpppasALL.CSPlusRatio+RvpppasALL.CSMinusRatio))>=.7;% & strcmp(RvpppasALL.Ninfo(:,1),'PA168Sess13_sorted-01.nex');%'PA167Sess07_sorted.nex''PA168Sess13_sorted-01.nex' & RvpppasALL.Ev(1).RespDir==0;% & CvpppasALL.STATSpost>=0.05;
    selectionPav=RvppasALL.Structure==10 & ~isnan(RvppasALL.Ev(LeftEv).Meanz) & ~isnan(RvppasALL.Ev(RightEv).Meanz) & RvppasALL.CSPlusRatio>=.5 & (RvppasALL.CSPlusRatio./(RvppasALL.CSPlusRatio+RvppasALL.CSMinusRatio))>=.7;% & RvppasALL.Ev(1).RespDir==0;% & CvppasALL.STATSpost>=0.05;
    
    %------------- SORTING -------------
    SortEv=LeftEv; NonSortEv=RightEv;
    if i==1;
        selection=selectionInst;
        TMP=-RvpppasALL.Ev(SortEv).Meanz(selection,ExcInh);
        %TMP=RvpppasALL.Ev(SortEv).RespDir(selection,ExcInh);
    else;
        selection=selectionPav;
        TMP=-RvppasALL.Ev(SortEv).Meanz(selection,ExcInh);
        %TMP=RvppasALL.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(RvpppasALL.Ninfo(selection,3)); NNid1=NNid1(SORTimg); end
    if i==2; NNid2=cell2mat(RvppasALL.Ninfo(selection,3)); NNid2=NNid2(SORTimg); end
    
    if i==1,
        selection=selectionInst; SortEv=LeftEv; NonSortEv=RightEv;
        % ------------- NORMALIZATION -------------
        if Norm==0
            imgLeft=RvpppasALL.Ev(LeftEv).PSTHz(selection,Ishow);
            imgRight=RvpppasALL.Ev(RightEv).PSTHz(selection,Ishow);
            E=[-200 400];Clim=sign(E).*abs(E).^(1/4);%ColorMap
        elseif Norm~=0
            imgLeft=normalize(RvpppasALL.Ev(LeftEv).PSTHz(selection,Ishow),RvpppasALL.Ev(LeftEv).PSTHz(selection,(RvpppasALL.Param.Tm<0)),Norm); %normalize to max response
            imgRight=normalize(RvpppasALL.Ev(RightEv).PSTHz(selection,Ishow),RvpppasALL.Ev(LeftEv).PSTHz(selection,(RvpppasALL.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=RvppasALL.Ev(LeftEv).PSTHz(selection,Ishow);
                imgRight=RvppasALL.Ev(RightEv).PSTHz(selection,Ishow);
                E=[-200 400];Clim=sign(E).*abs(E).^(1/4);%ColorMap
            elseif Norm~=0
                imgLeft=normalize(RvppasALL.Ev(LeftEv).PSTHz(selection,Ishow),RvppasALL.Ev(LeftEv).PSTHz(selection,(RvppasALL.Param.Tm<0)),Norm); %normalize to max response
                imgRight=normalize(RvppasALL.Ev(RightEv).PSTHz(selection,Ishow),RvppasALL.Ev(LeftEv).PSTHz(selection,(RvppasALL.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=RvpppasALL.CSStat;
        elseif LeftEv==3, Stat=RvpppasALL.CSPlusResponseStat;
        end
        selection1=selectionInst & Stat>=0.05;
        selection2=selectionInst & Stat<0.05 & RvpppasALL.Ev(LeftEv).Meanz>RvpppasALL.Ev(RightEv).Meanz;
        selection3=selectionInst & Stat<0.05 & RvpppasALL.Ev(LeftEv).Meanz<RvpppasALL.Ev(RightEv).Meanz;
        scatter(RvpppasALL.Ev(LeftEv).Meanz(selection1),RvpppasALL.Ev(RightEv).Meanz(selection1),30, 'k');%,'filled');
        scatter(RvpppasALL.Ev(LeftEv).Meanz(selection2),RvpppasALL.Ev(RightEv).Meanz(selection2),30, DarkPink);
        scatter(RvpppasALL.Ev(LeftEv).Meanz(selection3),RvpppasALL.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');
    else %PS
        if     LeftEv==1, Stat=RvppasALL.CSStat;
        elseif LeftEv==3, Stat=RvppasALL.CSResponseStat;
        end
        selection1=selectionPav & Stat>=0.05;
        selection2=selectionPav & Stat<0.05 & RvppasALL.Ev(LeftEv).Meanz>RvppasALL.Ev(RightEv).Meanz;
        selection3=selectionPav & Stat<0.05 & RvppasALL.Ev(LeftEv).Meanz<RvppasALL.Ev(RightEv).Meanz;
        scatter(RvppasALL.Ev(LeftEv).Meanz(selection1),RvppasALL.Ev(RightEv).Meanz(selection1),30, 'k');%,'filled');
        scatter(RvppasALL.Ev(LeftEv).Meanz(selection2),RvppasALL.Ev(RightEv).Meanz(selection2),30, DarkPink);
        scatter(RvppasALL.Ev(LeftEv).Meanz(selection3),RvppasALL.Ev(RightEv).Meanz(selection3),30, DarkGreen);
        Min=floor(min([RvppasALL.Ev(LeftEv).Meanz(selection2);RvppasALL.Ev(RightEv).Meanz(selection2)]));
        Max=ceil(max([RvppasALL.Ev(LeftEv).Meanz(selection2);RvppasALL.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');
    end
    Min=floor(min([RvpppasALL.Ev(LeftEv).Meanz(RvpppasALL.Structure==10);RvpppasALL.Ev(RightEv).Meanz(RvpppasALL.Structure==10)]));
    Max=ceil(max([RvpppasALL.Ev(LeftEv).Meanz(RvpppasALL.Structure==10);RvpppasALL.Ev(RightEv).Meanz(RvpppasALL.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=RvpppasALL.CSStat;
    elseif LeftEv==3, Stat=RvpppasALL.CSPlusResponseStat;
    end
    axis([-10 10 0 .25]);
    edges = [-10:1:9];
    selection1=selectionInst & Stat>=0.05;
    selection2=selectionInst & Stat<0.05 & RvpppasALL.Ev(LeftEv).Meanz>RvpppasALL.Ev(RightEv).Meanz;
    selection3=selectionInst & Stat<0.05 & RvpppasALL.Ev(LeftEv).Meanz<RvpppasALL.Ev(RightEv).Meanz;
    Delta1= RvpppasALL.Ev(LeftEv).Meanz(selection1)-RvpppasALL.Ev(RightEv).Meanz(selection1);
    Delta2= RvpppasALL.Ev(LeftEv).Meanz(selection2)-RvpppasALL.Ev(RightEv).Meanz(selection2);
    Delta3= RvpppasALL.Ev(LeftEv).Meanz(selection3)-RvpppasALL.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=RvppasALL.CSStat;
    elseif LeftEv==3, Stat=RvppasALL.CSrespStat;
    end
    axis([-10 10 0 .25]);
    edges = [-10:1:9];
    selection1=selectionPav & Stat>=0.05;
    selection2=selectionPav & Stat<0.05 & RvppasALL.Ev(LeftEv).Meanz>RvppasALL.Ev(RightEv).Meanz;
    selection3=selectionPav & Stat<0.05 & RvppasALL.Ev(LeftEv).Meanz<RvppasALL.Ev(RightEv).Meanz;
    Delta1= RvppasALL.Ev(LeftEv).Meanz(selection1)-RvppasALL.Ev(RightEv).Meanz(selection1);
    Delta2= RvppasALL.Ev(LeftEv).Meanz(selection2)-RvppasALL.Ev(RightEv).Meanz(selection2);
    Delta3= RvppasALL.Ev(LeftEv).Meanz(selection3)-RvppasALL.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(RvpppasALL.Ev(LeftEv).Meanz) & ~isnan(RvpppasALL.Ev(RightEv).Meanz)& RvpppasALL.Ev(1).RespDir==1;% & CvpppasALL.STATSpost>=0.05;
    selection2=selectionPav & ~isnan(RvppasALL.Ev(LeftEv).Meanz) & ~isnan(RvppasALL.Ev(RightEv).Meanz)& RvppasALL.Ev(1).RespDir==1;% & CvppasALL.STATSpost>=0.05;
    
    %subplot of DS vs NS excited neurons
    if i==1
        selection=selection1;
        subplot(6,4,i+16);
        psthLeft=nanmean(RvpppasALL.Ev(LeftEv).PSTHz(selection,Ishow),1);
        semLeft=nanste(RvpppasALL.Ev(LeftEv).PSTHz(selection,Ishow),1);
        upLeft=psthLeft+semLeft;
        downLeft=psthLeft-semLeft;
        hold on;
        
        plot(time,psthLeft,'Color',DS,'linewidth',1.5);
        
        psthRight=nanmean(RvpppasALL.Ev(RightEv).PSTHz(selection,Ishow),1);
        semRight=nanste(RvpppasALL.Ev(RightEv).PSTHz(selection,Ishow),1);
        upRight=psthRight+semRight;
        downRight=psthRight-semRight;
        
        plot(time,psthRight,'k','linewidth',1.5);
        %legend('DS','NS','Location','northwest');
        patch([time,time(end:-1:1)],[upLeft,downLeft(end:-1:1)],DS,'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)');
        
    else if i==2
            selection=selection2;
            subplot(6,4,i+19);
            psthLeft=nanmean(RvppasALL.Ev(LeftEv).PSTHz(selection,Ishow),1);
            semLeft=nanste(RvppasALL.Ev(LeftEv).PSTHz(selection,Ishow),1);
            upLeft=psthLeft+semLeft;
            downLeft=psthLeft-semLeft;
            hold on;
            
            plot(time,psthLeft,'Color',Pav,'linewidth',1.5);
            
            psthRight=nanmean(RvppasALL.Ev(RightEv).PSTHz(selection,Ishow),1);
            semRight=nanste(RvppasALL.Ev(RightEv).PSTHz(selection,Ishow),1);
            upRight=psthRight+semRight;
            downRight=psthRight-semRight;
            
            plot(time,psthRight,'k','linewidth',1.5);
            %legend('CS+','CS-','Location','northwest');
            patch([time,time(end:-1:1)],[upLeft,downLeft(end:-1:1)],Pav,'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)');
        end
    end
        
        
        selection1=selectionInst & ~isnan(RvpppasALL.Ev(LeftEv).Meanz) & ~isnan(RvpppasALL.Ev(RightEv).Meanz)& RvpppasALL.Ev(1).RespDir==-1;% & CvpppasALL.STATSpost>=0.05;
        selection2=selectionPav & ~isnan(RvppasALL.Ev(LeftEv).Meanz) & ~isnan(RvppasALL.Ev(RightEv).Meanz)& RvppasALL.Ev(1).RespDir==-1;% & CvppasALL.STATSpost>=0.05;
        
        if i==1
            selection=selection1;
            subplot(6,4,i+17);
            psthLeft=nanmean(RvpppasALL.Ev(LeftEv).PSTHz(selection,Ishow),1);
            semLeft=nanste(RvpppasALL.Ev(LeftEv).PSTHz(selection,Ishow),1);
            upLeft=psthLeft+semLeft;
            downLeft=psthLeft-semLeft;
            hold on;
            
            plot(time,psthLeft,'Color',DS,'linewidth',1.5);
            
            psthRight=nanmean(RvpppasALL.Ev(RightEv).PSTHz(selection,Ishow),1);
            semRight=nanste(RvpppasALL.Ev(RightEv).PSTHz(selection,Ishow),1);
            upRight=psthRight+semRight;
            downRight=psthRight-semRight;
            
            plot(time,psthRight,'k','linewidth',1.5);
            legend('DS','NS','Location','northwest');
            patch([time,time(end:-1:1)],[upLeft,downLeft(end:-1:1)],DS,'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)');
            
        else if i==2
                selection=selection2;
                subplot(6,4,i+20);
                psthLeft=nanmean(RvppasALL.Ev(LeftEv).PSTHz(selection,Ishow),1);
                semLeft=nanste(RvppasALL.Ev(LeftEv).PSTHz(selection,Ishow),1);
                upLeft=psthLeft+semLeft;
                downLeft=psthLeft-semLeft;
                hold on;
                
                plot(time,psthLeft,'Color',Pav,'linewidth',1.5);
                
                psthRight=nanmean(RvppasALL.Ev(RightEv).PSTHz(selection,Ishow),1);
                semRight=nanste(RvppasALL.Ev(RightEv).PSTHz(selection,Ishow),1);
                upRight=psthRight+semRight;
                downRight=psthRight-semRight;
                
                plot(time,psthRight,'k','linewidth',1.5);
                legend('CS+','CS-','Location','northwest');
                patch([time,time(end:-1:1)],[upLeft,downLeft(end:-1:1)],Pav,'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)');
            end
        end
    
end

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

task=cat(1,ones(length(RvpppasALL.CSauROC(selectionInst)),1),zeros(length(RvppasALL.CSauROC(selectionPav)),1));
rat=cat(1,RvpppasALL.Rat(selectionInst),RvppasALL.Rat(selectionPav));
tbl=table(cat(1,-RvpppasALL.CSauROC(selectionInst),-RvppasALL.CSauROC(selectionPav)),task,rat,'variablenames',{'auROC','Task','Subject'});
%lme=fitlme(tbl,'auROC~Task+(1|Subject)');%+(1|Subject:Exposure)');


%stats on firing of excited and inhibited neurons
excitedNsP=selectionPav & (RvppasALL.Ev(1).RFLAG==1 | RvppasALL.Ev(1).RespDir==1) & RvppasALL.Ev(1).RFLAG~=-1 & RvppasALL.Ev(1).RespDir~=-1;
excitedNsI=selectionInst & (RvpppasALL.Ev(1).RFLAG==1 | RvpppasALL.Ev(1).RespDir==1) & RvpppasALL.Ev(1).RFLAG~=-1 & RvpppasALL.Ev(1).RespDir~=-1;
inhibNsI=selectionInst & (RvpppasALL.Ev(1).RFLAG==-1 | RvpppasALL.Ev(1).RespDir==-1) & RvpppasALL.Ev(1).RFLAG~=1 & RvpppasALL.Ev(1).RespDir~=1;
inhibNsP=selectionPav & (RvppasALL.Ev(1).RFLAG==-1 | RvppasALL.Ev(1).RespDir==-1) & RvppasALL.Ev(1).RFLAG~=1 & RvppasALL.Ev(1).RespDir~=1;

% [~,p,~,stats]=ttest(RvppasALL.Ev(1).Meanraw(excitedNsP),RvppasALL.Ev(2).Meanraw(excitedNsP)); %CS excited
% [~,p,~,stats]=ttest(RvppasALL.Ev(1).Meanraw(inhibNsP),RvppasALL.Ev(2).Meanraw(inhibNsP)); %CS inhibited
% [~,p,~,stats]=ttest(RvpppasALL.Ev(1).Meanraw(excitedNsI),RvpppasALL.Ev(2).Meanraw(excitedNsI)); %DS excited
% [~,p,~,stats]=ttest(RvpppasALL.Ev(1).Meanraw(inhibNsI),RvpppasALL.Ev(2).Meanraw(inhibNsI)); %DS inhibited
%---------------- COLOR SCALE ----------------------
subplot(10,6,60);
imagesc(time,[1 size(imgLeft,1)],imgLeft,Clim);colorbar;
%ylabel('z-score');

%
% if SAVEFIG==1
%     Figpath='C:\Users\jricha95\Documents\Documents (3)\Experiments\New Experiments\PS PPS';
%     exportfig(figure(1), sprintf('%s\\Figure1cmyk.eps',Figpath), 'FontMode', 'fixed','FontSize', 12,'LineWidth', 0.5, 'Width', 10, 'format','eps2','color','cmyk','renderer', 'painters','resolution',300);
% %     %end





%% Figure 5 pie charts

figure;


%Pie charts of response patterns
%DS
subplot(4,5,11);
exc=sum((RvpppasALL.Ev(1).RFLAG(selectionInst)==1 | RvpppasALL.Ev(1).RespDir(selectionInst)==1) & RvpppasALL.Ev(1).RFLAG(selectionInst)~=-1 & RvpppasALL.Ev(1).RespDir(selectionInst)~=-1);
inh=sum(RvpppasALL.Ev(1).RespDir(selectionInst)==-1 | RvpppasALL.Ev(1).RFLAG(selectionInst)==-1);
no=sum(RvpppasALL.Ev(1).RFLAG(selectionInst)==0 & RvpppasALL.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((RvppasALL.Ev(1).RFLAG(selectionPav)==1 | RvppasALL.Ev(1).RespDir(selectionPav)==1) & RvppasALL.Ev(1).RFLAG(selectionPav)~=-1 & RvppasALL.Ev(1).RespDir(selectionPav)~=-1);
inh=sum(RvppasALL.Ev(1).RespDir(selectionPav)==-1 | RvppasALL.Ev(1).RFLAG(selectionPav)==-1);
no=sum(RvppasALL.Ev(1).RFLAG(selectionPav)==0 & RvppasALL.Ev(1).RespDir(selectionPav)==0);
response=[exc inh no];
pie(response);
title('CS+');

%DS port entry
subplot(4,5,13);
exc=sum((RvpppasALL.Ev(8).RFLAG(selectionInst)==1 | RvpppasALL.Ev(8).RespDir(selectionInst)==1) & RvpppasALL.Ev(5).RFLAG(selectionInst)~=-1 & RvpppasALL.Ev(5).RespDir(selectionInst)~=-1);
inh=sum(RvpppasALL.Ev(8).RespDir(selectionInst)==-1 | RvpppasALL.Ev(8).RFLAG(selectionInst)==-1);
no=sum(RvpppasALL.Ev(8).RFLAG(selectionInst)==0 & RvpppasALL.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((RvppasALL.Ev(8).RFLAG(selectionPav)==1 | RvppasALL.Ev(8).RespDir(selectionPav)==1) & RvppasALL.Ev(5).RFLAG(selectionPav)~=-1 & RvppasALL.Ev(5).RespDir(selectionPav)~=-1);
inh=sum(RvppasALL.Ev(8).RespDir(selectionPav)==-1 | RvppasALL.Ev(8).RFLAG(selectionPav)==-1);
no=sum(RvppasALL.Ev(8).RFLAG(selectionPav)==0 & RvppasALL.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((RvppasALL.Ev(9).RFLAG(selectionPav)==1 | RvppasALL.Ev(9).RespDir(selectionPav)==1) & RvppasALL.Ev(7).RFLAG(selectionPav)~=-1 & RvppasALL.Ev(7).RespDir(selectionPav)~=-1);
inh=sum(RvppasALL.Ev(9).RespDir(selectionPav)==-1 | RvppasALL.Ev(9).RFLAG(selectionPav)==-1);
no=sum(RvppasALL.Ev(9).RFLAG(selectionPav)==0 & RvppasALL.Ev(9).RespDir(selectionPav)==0);
response=[exc inh no];
h=pie(response);
title('post-CS+ Reward');

%% Figure 5 behavior
clear all;

load('RAWvppsALL.mat')
load('RAWvpppsALL.mat')
load('RAWvppasALL.mat')
load('RAWvpppasALL.mat')
load('RPPSall.mat')
load('RPSall.mat')
load('RvppasALL.mat')
load('RvpppasALL.mat')


SessRat=[];
for task=1:4
    if task==1 %inst suc
        RAW=RAWvpppsALL;
        R=RPPSall;
    elseif task==2 %inst suc(alc)
        RAW=RAWvpppasALL;
        R=RvpppasALL;
    elseif task==3 %pav suc
        RAW=RAWvppsALL;
        R=RPSall;
    elseif task==4 %pav suc(alc)
        RAW=RAWvppasALL;
        R=RvppasALL;
    end
    
    %which neurons are included
    selection=R.Structure==10 & R.CSPlusRatio>=.5 & (R.CSPlusRatio./(R.CSPlusRatio+R.CSMinusRatio))>=.7;
    
    %which sessions are included
    allsessions=unique(R.Ninfo(:,1)); %all sessions from this task
    [includedsessions,index]=unique(R.Ninfo(selection,1)); %all included sessions
    included=ismember(allsessions,includedsessions); %which of all sessions are included
    
    %get rat for each session
    Rats=R.Rat(index)+task*100; %make sure rats with same name from different experiment aren't accidentally combined
    SessRat=cat(1,SessRat,Rats,Rats); %twice to represent the CSPlus data and then the CSMinus data

    
    
    RAW=RAW(included);
    
    for i=1:length(RAW)
        if task==1 | task==3
            Ev1=strcmp('CSPlus1', RAW(i).Einfo(:,2));
            Ev2=strcmp('CSMinus1', RAW(i).Einfo(:,2));
        else
            Ev1=strcmp('CSPlus', RAW(i).Einfo(:,2));
            Ev2=strcmp('CSMinus', RAW(i).Einfo(:,2));
        end
        PE=strcmp('PortEntry', RAW(i).Einfo(:,2));
        
        CSPluses=RAW(i).Erast{Ev1};
        CSMinuses=RAW(i).Erast{Ev2};
        PortEntries=RAW(i).Erast{PE};
        
        CSPLatency=[];
        for j=1:length(CSPluses)
            if sum(PortEntries>CSPluses(j))>0
                CSPLatency(j,1)=min(PortEntries(PortEntries>CSPluses(j)))-CSPluses(j);
            else
                CSPLatency(j,1)=100;
            end
        end
        CSPLatency(CSPLatency>10)=[];
        
        CSMLatency=[];
        for j=1:length(CSMinuses)
            if sum(PortEntries>CSMinuses(j))>0
                CSMLatency(j,1)=min(PortEntries(PortEntries>CSMinuses(j)))-CSMinuses(j);
            else
                CSMLatency(j,1)=100; %will be removed
            end
        end
        CSMLatency(CSMLatency>10)=[];        
        
        CSPMean{task,1}(i,1)=mean(CSPLatency,'omitnan');
        CSMMean{task,1}(i,1)=mean(CSMLatency,'omitnan');
        CSPlusProb{task,1}(i,1)=length(CSPLatency)/length(CSPluses);
        CSMinusProb{task,1}(i,1)=length(CSMLatency)/length(CSMinuses);
        
    end
    
end

%% plotting

figure;

DS=[0    0.4470    0.7410];
Pav=[0.8500    0.3250    0.0980];
DSAlcohol=[51/255 160/255 44/255];
PavAlcohol=[180/255 41/255 230/255];
black=[0 0 0];
Colors=[DS;black;DSAlcohol;black;Pav;black;PavAlcohol;black];


%probability
data=[CSPlusProb{1,1};CSMinusProb{1,1};CSPlusProb{2,1};CSMinusProb{2,1};CSPlusProb{3,1};CSMinusProb{3,1};CSPlusProb{4,1};CSMinusProb{4,1}];
group=[ones(size(CSPlusProb{1,1}));2*ones(size(CSMinusProb{1,1}));3*ones(size(CSPlusProb{2,1}));4*ones(size(CSMinusProb{2,1}));5*ones(size(CSPlusProb{3,1}));6*ones(size(CSMinusProb{3,1}));7*ones(size(CSPlusProb{4,1}));8*ones(size(CSMinusProb{4,1}))];
task=[ones(size(CSPlusProb{1,1}));ones(size(CSMinusProb{1,1}));ones(size(CSPlusProb{2,1}));ones(size(CSMinusProb{2,1}));2*ones(size(CSPlusProb{3,1}));2*ones(size(CSMinusProb{3,1}));2*ones(size(CSPlusProb{4,1}));2*ones(size(CSMinusProb{4,1}))];
exposure=[ones(size(CSPlusProb{1,1}));ones(size(CSMinusProb{1,1}));2*ones(size(CSPlusProb{2,1}));2*ones(size(CSMinusProb{2,1}));ones(size(CSPlusProb{3,1}));ones(size(CSMinusProb{3,1}));2*ones(size(CSPlusProb{4,1}));2*ones(size(CSMinusProb{4,1}))];
cue=[ones(size(CSPlusProb{1,1}));2*ones(size(CSMinusProb{1,1}));ones(size(CSPlusProb{2,1}));2*ones(size(CSMinusProb{2,1}));ones(size(CSPlusProb{3,1}));2*ones(size(CSMinusProb{3,1}));ones(size(CSPlusProb{4,1}));2*ones(size(CSMinusProb{4,1}))];

subplot(1,2,1);
hold on;
boxplot(data,group,'colors',Colors,'symbol','');

for i=1:8
    scatter(i+((rand(sum(group==i),1)-0.5)/5),data(group==i),[],Colors(i,:),'filled')
end

plot([4.5 4.5],[0 0.9],'--','color','k');
title('Port Entry Probability');
xticks(1:8);
xticklabels({'DS','NS','DS','NS','CS+','CS-','CS+','CS-'});
text(0.8,0.6,'Sucrose','color',DS);
text(2.8,0.8,'Suc(alc)','color',DSAlcohol);
text(4.8,0.5,'Sucrose','color',Pav);
text(6.8,0.6,'Suc(alc)','color',PavAlcohol);
text(4.3,0.97,'ns','color','k');
text(2,0.9,'<- ns ->','color','k','horizontalalignment','center');
text(6,0.85,'<- ns ->','color','k','horizontalalignment','center');
axis([0 9 0 1]);

%stats
%[~,ProbAnova,ProbStats]=anovan(data,{task,cue,SessRat},'varnames',{'Task','Cue','Subject'},'random',3,'nested',[0 0 0;0 0 0;1 0 0],'model','full');
[~,ProbAnova,ProbStats]=anovan(data,{task,exposure,cue},'varnames',{'Task','Exposure','Cue'},'model','full','display','off');
%[~,ProbAnova,ProbStats]=anovan(data(cue==1),{task(cue==1),exposure(cue==1)},'varnames',{'Task','Exposure'},'model','full','display','off');
%multcompare(ProbStats,'dimension',[1 2]);

%latency
data=[CSPMean{1,1};CSMMean{1,1};CSPMean{2,1};CSMMean{2,1};CSPMean{3,1};CSMMean{3,1};CSPMean{4,1};CSMMean{4,1}];
group=[ones(size(CSPMean{1,1}));2*ones(size(CSMMean{1,1}));3*ones(size(CSPMean{2,1}));4*ones(size(CSMMean{2,1}));5*ones(size(CSPMean{3,1}));6*ones(size(CSMMean{3,1}));7*ones(size(CSPMean{4,1}));8*ones(size(CSMMean{4,1}))];

subplot(1,2,2);
hold on;
boxplot(data,group,'colors',Colors,'symbol','');

for i=1:8
    scatter(i+((rand(sum(group==i),1)-0.5)/5),data(group==i),[],Colors(i,:),'filled')
end

plot([4.5 4.5],[0 8],'--','color','k');
title('Port Entry Latency (sec)');
xticks(1:8);
xticklabels({'DS','NS','DS','NS','CS+','CS-','CS+','CS-'});
xticklabels({'DS','NS','DS','NS','CS+','CS-','CS+','CS-'});
text(0.8,7,'Sucrose','color',DS);
text(2.8,6,'Suc(alc)','color',DSAlcohol);
text(4.8,7.5,'Sucrose','color',Pav);
text(6.8,7,'Suc(alc)','color',PavAlcohol);
text(2,1.5,'<- ns ->','color','k','horizontalalignment','center');
text(6,1.5,'<- ns ->','color','k','horizontalalignment','center');
text(4.3,8.5,'*','color','k');

axis([0 9 0 9]);

%stats
%[~,LatAnova,LatStats]=anovan(data,{task,cue,SessRat},'varnames',{'Task','Cue','Subject'},'random',3,'nested',[0 0 0;0 0 0;1 0 0],'model','full');
[~,LatAnova,LatStats]=anovan(data,{task,exposure,cue},'varnames',{'Task','Exposure','Cue'},'model','full','display','off');
%[~,LatAnova,LatStats]=anovan(data(cue==1),{task(cue==1),exposure(cue==1)},'varnames',{'Task','Exposure'},'model','full','display','off');

%multcompare(LatStats,'dimension',[1 2]);