%% generate Figure 10 and Figure 10 supplement 1 - Correlation analysis - Correlation analysis - RR vs FR

%firing analysis during LP and correl with response rate
%Analysis on all neurons without DLS/DMS dissociation
%histogram for correlated neurons in each regions, in each class

clear all;clc;
tic
global Dura Baseline Tm Tbase BSIZE Tbin
SAVE_FLAG=1;
BSIZE=0.05;
Dura=[-25 60]; Tm=-1:BSIZE:60;
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) -20];
postwin=[0 .5];
EventWin=[-0.5 0.5];
postwin2=[Dura(1):0.05:Dura(2)]; %bounds should match Dura
Slopebounds=[-1.5:0.25:1.5];
R2Bounds=[0:0.05:0.5];
PvalBounds=[0:0.01:0.5];
Struct=[10 20];
XI=1;
RunRegress=0;
SAVEFIG=1;
%R=[];R.Ninfo={}; 
NN=0;
pre=[]; post=[];
if ~exist('RAW'), load ('RAWextendedtraining.mat'); RAW=RAW; end
if ~exist('Ev'), load('Rextendedtraining_light.mat'); end
load('C.mat');

%Smoothing
Smoothing=1; %0 for raw and 1 for smoothing
SmoothTYPE='lowess';
SmoothSPAN=5;
if Smoothing~=1, SmoothTYPE='NoSmoothing';SmoothSPAN=NaN; end


%path='C:\FRED\GALLO\Exps\2009-present GoNoGo\FINAL\PLX Files GoNoGo\GoNoGo NEX files\RESULTgn14.xls';
%%

if RunRegress
  % PREALLOCATING
    for i=1:length(RAW)
        EvInd=strmatch('First_LP',RAW(i).Einfo(:,2),'exact'); %finds the index of the event
        Rind=strmatch('ReliableLP_DT',RAW(i).Einfo(:,2),'exact');  %finds the index of the event
        
        if  sum(EvInd)~=0 && sum(Rind)~=0
            DS(i)=length(RAW(i).Erast{EvInd});
            Neur(i)=length(RAW(i).Nrast);
        end
    end
    MaxTrials=max(DS); MaxNeur=sum(Neur); MaxWin=length(postwin);
%     C.Lat=NaN(MaxNeur, MaxTrials); % (neurons, trials)
    C.preRR=NaN(MaxNeur,MaxTrials);  % (neurons, trials)
    C.SEQ=NaN(MaxNeur,MaxTrials); % (neurons, trials)
    C.postwinRR=NaN(MaxNeur,MaxTrials,MaxWin); % (neurons, trials, windows)
    
    
    %
    for i=1:length(RAW) %loops through sessions
        EvInd=strmatch('First_LP',RAW(i).Einfo(:,2),'exact'); %finds the index of the event
        Rind=strmatch('ReliableLP_DT',RAW(i).Einfo(:,2),'exact');  %finds the index of the event
        
        if  sum(EvInd)~=0 && sum(Rind)~=0            
            Dcell2=MakePSR05(RAW(i).Erast(Rind),RAW(i).Erast{EvInd},[-1 60],{2}); % makes trial by trail 
                    
            for j= 1:size(RAW(i).Nrast,1) %Number of neurons within sessions
                NN=NN+1;

                [PSRNeur,Nneur]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},[-1 60],{1});% collapsed raster
                [PSR2,N2]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},Dura,{2});% makes trial by trail rasters. PSR1 is a cell(neurons, trials)
                
                [PTHneur,BW1,~]=MakePTH07(PSRNeur,repmat(N2, size(RAW(i).Nrast{j},1),1),{2,0,BSIZE});
                PTHneur=smooth(PTHneur,SmoothSPAN,SmoothTYPE)';               
                C.PSTHraw(NN,1:length(Tm))=PTHneur;
                Search=find(Tm>=EventWin(1) & Tm<=EventWin(2));
                [C.MaxVal(NN,1),MaxInd]=max(C.PSTHraw(NN,Search));
                C.MaxTime(NN,1)=Tm(Search(1)+MaxInd-1);
                
                for m=1:size(RAW(i).Erast{EvInd},1) %loops through trials
                    SequenceDur=RAW(i).Eint{1,6}(m)-RAW(i).Eint{1,5}(m);
                    C.RR(NN,m)=length(Dcell2{1,m})/SequenceDur;
                    C.preRR(NN,m)=sum(PSR2{1,m}<prewin(2) & PSR2{1,m}>prewin(1))/(abs(prewin(2)-prewin(1))); %neurons, trials
                    C.SEQ(NN,m)=sum(PSR2{1,m}<SequenceDur & PSR2{1,m}>0)/SequenceDur;
                end %m loop
                fprintf('Neuron ID # %d\n',NN);
            end %j loop
        end %if
    end %i loop
    
 

    
    %% Runs the regression analysis
    for NN=1:MaxNeur
        RegX=[C.RR(NN,:)', ones(size(C.RR,2),1)];
        [C.SLOPEpreRR,C.STATSpreRR]=corr(C.preRR(NN,:)',C.RR(NN,:)','rows','pairwise','type','Spearman');
        [C.SLOPEseq(NN,:),C.STATSseq(NN,:)]= corr(C.SEQ(NN,:)',C.RR(NN,:)','rows','pairwise','type','Spearman');
        [C.SLOPEseqpre(NN,:),C.STATSseqpre(NN,:)]= corr((C.SEQ(NN,:)-C.preRR(NN,:))',C.RR(NN,:)','rows','pairwise','type','Spearman');
        fprintf('CORREL Neuron ID # %d\n',NN);
    end
    
    C.SelectionRR(:,1)=Coord(:,4)==10 & C.STATSseq(:,1)<PStat;
    C.SelectionRR(:,2)=Coord(:,4)==20 & C.STATSseq(:,1)<PStat;
    % Saving data in C structure and excel file
    save('C.mat', 'C');

end


%% --------------Plots distribution of correlation variables SINGLE WINDOW -------------
% POST WINDOW
figure;
Struct=[10 20];
clear A;

for i=1:2 % loop through structure
    clear sel1;clear sel2;    
    sel1=Coord(:,4)==i*10; %(R.Class(:,1)==-4); %& R.DSNSStat<PStat & (R.Ev(1).Meanz>R.Ev(2).Meanz);
    sel2= sel1 & C.STATSseq(:,1)<PStat;
    
    
    A(i,:)=cat(2, sum(Coord(:,4)==i*10), sum(sel1), sum(sel2));
        
    h=subplot(2,2,i);
    xmin=floor(min(C.SLOPEseq(sel1,1)));
    xmax=ceil(max(C.SLOPEseq(sel1,1)));
    step=0.05;
    xcenters=xmin:step:xmax;
    Nsel1=hist(C.SLOPEseq(sel1,1),xcenters);
    Nsel2=hist(C.SLOPEseq(sel2,1),xcenters);
    hold on;
    bar(xcenters, Nsel1, 'w', 'EdgeColor', 'k', 'BarWidth', 1);alpha(0.3);
    bar(xcenters, Nsel2, 'k', 'EdgeColor', 'k', 'BarWidth', 1);
    plot([0 0], [0 30],'LineStyle',':','Color','k');
    pval1=signrank(C.SLOPEseq(sel1,1));
    pval2=signrank(C.SLOPEseq(sel2,1));
    MEANsel1=mean(C.SLOPEseq(sel2,1));
    plot([MEANsel1 MEANsel1], [0 10],'Color','r');

        
    subplot(2,2,i+2)
    [N,xcenters2]=hist(C.STATSseq(sel2,1),20);
    bar(xcenters2, N, 'k', 'EdgeColor', 'k', 'BarWidth', 1); alpha(0.3);
    axis ([0 ceil(max(xcenters2)*10)/10 0 ceil(max(N))+1]);
end

A(:,4)=A(:,3)./A(:,2)*100;


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