%% generate Figure 10 and Figure 10 supplement 1 - Correlation analysis

%firing analysis post lever insertion and correl with response latency
%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.01;
Dura=[-2 4]; 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=[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;
XI=1;
RunRegress=1;
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');

%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('LeverInsertion',RAW(i).Einfo(:,2),'exact'); %finds the index of the event
        Rind=strmatch('First_LP',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.pre=NaN(MaxNeur,MaxTrials);  % (neurons, trials)
    C.post=NaN(MaxNeur,MaxTrials); % (neurons, trials)
    C.postwin=NaN(MaxNeur,MaxTrials,MaxWin); % (neurons, trials, windows)
    
    
    %
    for i=1:length(RAW) %loops through sessions
        EvInd=strmatch('LeverInsertion',RAW(i).Einfo(:,2),'exact'); %finds the index of the event
        Rind=strmatch('First_LP',RAW(i).Einfo(:,2),'exact');  %finds the index of the event
        
        if  sum(EvInd)~=0 && sum(Rind)~=0
            Dcell=MakePSR04(RAW(i).Erast(Rind),RAW(i).Erast{EvInd},[-1 60],{2,'first'});
            D=cell2mat(Dcell); %inds=find(~isnan(D));
            D(isnan(D))=60; %trials with no response are set to 10.25 to allow to use them in the correlation
            for j= 1:size(RAW(i).Nrast,1) %Number of neurons within sessions
                NN=NN+1;
                C.Lat(NN,1:length(RAW(i).Erast{EvInd}))=D;
                [PSR2,N2]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},Dura,{2});% makes trial by trail rasters. PSR1 is a cell(neurons, trials)
                for m=1:size(RAW(i).Erast{Rind},1) %loops through trials
                    Resplat=RAW(i).Eint{1,8}(m)-RAW(i).Eint{1,7}(m);
                    C.pre(NN,m)=sum(PSR2{1,m}<prewin(2) & PSR2{1,m}>prewin(1))/(abs(prewin(2)-prewin(1))); %neurons, trials
                    C.postLAT1(NN,m)=sum(PSR2{1,m}<postwin(2) & PSR2{1,m}>postwin(1))/(abs(postwin(2)-postwin(1)));
                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.Lat(NN,:)', ones(size(C.Lat,2),1)];
        [C.SLOPEpre,C.STATSpre]=corr(C.pre(NN,:)',C.Lat(NN,:)','rows','pairwise','type','Spearman');
        [C.SLOPEpostLAT1(NN,:),C.STATSpostLAT1(NN,:)]= corr(C.postLAT1(NN,:)',C.Lat(NN,:)','rows','pairwise','type','Spearman');
        [C.SLOPEpostpre(NN,:),C.STATSpostpre(NN,:)]= corr((C.post(NN,:)-C.pre(NN,:))',C.Lat(NN,:)','rows','pairwise','type','Spearman');

        fprintf('CORREL Neuron ID # %d\n',NN);
    end
    
    % Saving data in C structure and excel file
    save('C.mat', 'C');
    Cmat=[];
%    Cmat=cat(2,C.STATSpre, C.SLOPEpre, C.STATSpost, C.SLOPEpost);
    %xlswrite(path,Cmat,'Correl', 'F3');
end

%% --------------Plots distribution of correlation variables SINGLE WINDOW -------------
% POST WINDOW

XI=1; clear A;

for i=1:2
    figure(1)
    sel1=Coord(:,4)==i*10;
    sel2= sel1 & C.STATSpostLAT1(:,1)<PStat;
    A(1,:)=cat(2, sum(sel1), sum(sel2));     
    h=subplot(2,2,i);
    xmin=floor(min(C.SLOPEpostLAT1(sel1,1)));
    xmax=ceil(max(C.SLOPEpostLAT1(sel1,1)));
    step=0.05;
    xcenters=xmin:step:xmax;
    Nsel1=hist(C.SLOPEpostLAT1(sel1,1),xcenters);
    Nsel2=hist(C.SLOPEpostLAT1(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');
    %axis([-10 10 0 20]) %note: their is an outlier @-14 for shell
    pval1=signrank(C.SLOPEpostLAT1(sel1,1));
    pval2=signrank(C.SLOPEpostLAT1(sel2,1));
    MEANsel1=mean(C.SLOPEpostLAT1(sel2,1));
    plot([MEANsel1 MEANsel1], [0 10],'Color','r');
    %text(-5.8,9.8,sprintf('p=%d', pval1), 'Parent',h);
    %text(-5.8,8.8,sprintf('p=%d', pval2), 'Parent',h);    
    subplot(2,2,i+2)
    [N,xcenters2]=hist(C.STATSpostLAT1(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
save('C.mat', 'C');