%% initial analysis of event-evoked activity - Early training dataset with session by session hierarchy

% need RAW datafiles and Coord datafile with electrode placement
% need to save excel sheets in the specified folder
% the excel sheets say which events to make PSTHs for and what windows to
% use for statistical testing


%clear all; clc;

clear all; clc;
global Dura Baseline Tm Tbase BSIZE Tbin DuraITI TmITI
tic
path='C:\Users\yvandae1\Documents\MATLAB\DT Nex Files\RESULTdt.xls'; %excel file with windows etc


%Main settings
SAVE_FLAG=1;
BSIZE=0.01; %Do not change
Dura=[-25 20]; Tm=Dura(1):BSIZE:Dura(2);
DuraITI=[-50 20]; TmITI=DuraITI(1):BSIZE:DuraITI(2);
%Baseline=[-22 0]; Tbase=Baseline(1):BSIZE:Baseline(2); %now defined line 98
Tbin=-0.5:0.005:0.5; %window used to determine the optimal binsize
PStat=0.05; %for comparing pre versus post windows, or event A versus event B
MinNumTrials=10;
R=[];R.Ninfo={};NN=0;
Nneurons=[];
Neur=0;
BinSize=0.05;

load('RAWearlytraining.mat');
RAW=RAW;

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


% List of events to analyze and analysis windows EXTRACTED from excel file
[~,Erefnames]=xlsread(path,'Windows','a3:a13'); % cell that contains the event names
prewin =  xlsread(path,'Windows','b3:c13'); %for DT5
EventWin =  xlsread(path,'Windows','l3:m13'); %to compute MeanZ in relevant window
PreEventWin = xlsread(path,'Windows','d3:e13');
PostEventWin = xlsread(path,'Windows','f3:g13');
BLWin=    xlsread(path,'Windows','h3:i13');
RespWin=  xlsread(path,'Windows','j3:k13');
%[~,Session] = xlsread(path,'Windows','a16:a33');
%[~,Session] = xlsread(path,'Windows','c18:c33');
[~,Session] = xlsread(path,'Windows','c18:c30');

%Settings for Response detection w/ signrank
WINin=[0.03,0.1]; %Onset and offset requirements in sec.
%WINin=[-2 -4]; %negative values define onset requirement by the number of consecutive bins and not by duration
CI=[.95,.95];

%%

R.Erefnames= Erefnames;
R.Session= Session;
SessionName = {RAW.SessionName};
Ncum=0;

%Finds the total number of neurons accross all sessions
for i=1:length(Session) %loops through sessions
    Ses=strcmp(Session(i),SessionName);
    SesIndex=find(Ses==1);
    NN=0; Neur=0;
    R.Ses(i).Ninfo=[];
    
    for y=1:length(SesIndex)
        R.Ses(i).Ninfo=cat(1,R.Ses(i).Ninfo,RAW(SesIndex(y)).Ninfo);
        Neur=Neur+size(RAW(SesIndex(y)).Nrast,1);
        Nneurons(i)=Neur;
    end
    
    load('Coord_earlytraining.mat');
    R.Ses(i).Coord=part(i).Coord;
    %R.Ses(i).Coord=ones(Nneurons(i),4); %comment this line and uncomment the 2 lines above when you have real coordinates.
    R.Ses(i).Structure=part(i).Coord(:,4);
    R.Ses(i).Ninfo=cat(2, R.Ses(i).Ninfo, mat2cell([1:Nneurons(i)]',ones(1,Nneurons(i)),1));
    
    for k=1:length(Erefnames)
        R.Ses(i).Ev(k).PSTHraw(1:Nneurons(i),1:length(Tm))=NaN(Nneurons(i),length(Tm));
        R.Ses(i).Ev(k).PSTHrawBL(1:Nneurons(i),1:501)=NaN(Nneurons(i),501);
        R.Ses(i).Ev(k).PSTHz(1:Nneurons(i),1:length(Tm))=NaN(Nneurons(i),length(Tm));
        R.Ses(i).Ev(k).rawMeanz(1:Nneurons(i),1)=NaN;
        R.Ses(i).Ev(k).rawMeanzPre(1:Nneurons(i),1)=NaN;
        R.Ses(i).Ev(k).Meanraw(1:Nneurons(i),1)=NaN;
        R.Ses(i).Ev(k).Meanz(1:Nneurons(i),1)=NaN;
        R.Ses(i).Ev(k).MeanzPRE(1:Nneurons(i),1)=NaN;
        R.Ses(i).Ev(k).BW(1:Nneurons(i),1)=NaN;
        R.Ses(i).Ev(k).ttestPreEvent(1:Nneurons(i),1)=NaN;
        R.Ses(i).Ev(k).ttestPostEvent(1:Nneurons(i),1)=NaN;
        R.Ses(i).Ev(k).RespDirPre(1:Nneurons(i),1)=NaN;
        R.Ses(i).Ev(k).RespDirPost(1:Nneurons(i),1)=NaN;
        R.Ses(i).Ev(k).NumberTrials(1:Nneurons(i),1)=NaN;
        R.Ses(i).Ev(k).MaxVal(1:Nneurons(i),1)=NaN;
        R.Ses(i).Ev(k).MaxTime(1:Nneurons(i),1)=NaN;
    end
end


% preallocating
R.Param.Tm=Tm;
R.Param.Tbin=Tbin;
R.Param.Dura=Dura;
R.Param.Baseline=Baseline;
R.Param.PStat=PStat;
R.Param.MinNumTrials=MinNumTrials;
R.Param.path=path;
R.Param.prewin=prewin;
R.Param.postwinPreEvent=PreEventWin;
R.Param.postwinPostEvent=PostEventWin;
R.Param.BLwin=BLWin;
R.Param.RespWin=RespWin;
R.Param.ResponseReq=WINin;
R.Param.SmoothTYPE=SmoothTYPE;
R.Param.SmoothSPAN=SmoothSPAN;

NS=0;
TNN=0;

%% runs the main routine

for i=1:length(Session) %loops through sessions
    Ses=strcmp(Session(i),SessionName);
    SesIndex=find(Ses==1);
    NN=0;
    NS=NS+1;
    
    for y=1:length(SesIndex)
        for j= 1:size(RAW(SesIndex(y)).Nrast,1) %Number of neurons per session
            NN=NN+1; %neuron counter
            TNN=TNN+1; %neuron counter across all session
            
            if R.Ses(i).Structure(NN)~=0 %to avoid analyzing excluded neurons
                
                ev2=NaN(size(RAW(SesIndex(y)).Eint{1,3},1),2);
                [PSR1iti,N1iti]=MakePSR0int(RAW(SesIndex(y)).Nrast(j),RAW(SesIndex(y)).Eint{1,1},RAW(SesIndex(y)).Eint{1,2},{2},{1});% makes trial by trial rasters for baseline based on ITI
                [PSR1trial,N1trial]=MakePSR0int(RAW(SesIndex(y)).Nrast(j),RAW(SesIndex(y)).Eint{1,3},RAW(SesIndex(y)).Eint{1,4},{2},{2});% makes trial by trial rasters for each trial based on time interval
                
                for l=1:length(RAW(SesIndex(y)).Eint{1,3}) %loops through trials
                    ev2(l,1)=sum(PSR1iti{l}<RAW(SesIndex(y)).Eint{1,2}(l) & PSR1iti{l}>RAW(SesIndex(y)).Eint{1,1}(l))/(RAW(SesIndex(y)).Eint{1,2}(l)-RAW(SesIndex(y)).Eint{1,1}(l)-20);
                    ev2(l,2)=sum(PSR1trial{l}<RAW(SesIndex(y)).Eint{1,4}(l) & PSR1trial{l}>RAW(SesIndex(y)).Eint{1,3}(l))/(RAW(SesIndex(y)).Eint{1,4}(l)-RAW(SesIndex(y)).Eint{1,3}(l));
                end
                
                R.WF(TNN,:)=RAW(SesIndex(y)).waveforms(j,:);
                R.WFanalysis(TNN,1)=mean(ev2(:,1)); % mean FR during Baseline in the middle of each ITI (first and last 10s excluded)
                proportionISI=0;
                for k=3:1:size(RAW(SesIndex(y)).Nrast{j,1},1) % loop thru timestamps waveform
                    ISIn=RAW(SesIndex(y)).Nrast{j,1}(k)-RAW(SesIndex(y)).Nrast{j,1}(k-1); 
                    ISIo=RAW(SesIndex(y)).Nrast{j,1}(k-1)-RAW(SesIndex(y)).Nrast{j,1}(k-2);
                    calcul(k-2)=2*abs(ISIn-ISIo)/(ISIn+ISIo);
                    if ISIo>0.5
                        proportionISI=proportionISI+1;
                    end
                end
                R.CV2(TNN,1)=mean(calcul);%CV2
                R.CV2(TNN,2)=100*proportionISI/(size(RAW(SesIndex(y)).Nrast{j,1},1)-1);% of ISI>0.5s
                R.SESSION(TNN,1)=NS;
                
                for k=1:length(Erefnames) %loops thorough the events
                    EvInd=strcmp(Erefnames(k),RAW(SesIndex(y)).Einfo(:,2)); %find the event id number from RAW
                    if sum(EvInd)==0
                        fprintf('HOWDY, CANT FIND EVENTS FOR ''%s''\n',Erefnames{k});
                    end
                    
                    R.Ses(i).Ev(k).NumberTrials(NN,1)=length(RAW(SesIndex(y)).Erast{EvInd});
                    if  ~isempty(EvInd) && R.Ses(i).Ev(k).NumberTrials(NN,1)>= MinNumTrials
                        
                        Tbase=prewin(k,1):BSIZE:prewin(k,2);
                        [PSR0,N0]=MakePSR04(RAW(SesIndex(y)).Nrast(j),RAW(SesIndex(y)).Erast{EvInd},prewin(k,:),{1});% makes collapsed rasters for baseline.
                        [PSR1,N1]=MakePSR04(RAW(SesIndex(y)).Nrast(j),RAW(SesIndex(y)).Erast{EvInd},Dura,{1});% makes collpased rasters. PSR1 is a cell(neurons)
                        %                     if strcmp(Erefnames(k),'BaselineITI')
                        %                         [PSR2,N2]=MakePSR04(RAW(SesIndex(y)).Nrast(j),RAW(SesIndex(y)).Erast{EvInd},DuraITI,{2});% makes trial by trail rasters for BL based on ITI. PSR1 is a cell(neurons, trials)
                        %                     else
                        [PSR2,N2]=MakePSR04(RAW(SesIndex(y)).Nrast(j),RAW(SesIndex(y)).Erast{EvInd},Dura,{2});% makes trial by trail rasters. PSR1 is a cell(neurons, trials)
                        
                        if ~isempty(PSR0{1}) || ~isempty(PSR1{1}) %to avoid errors, added on 12/28 2011
                            
                            %Fixed bin size
                             [PTH1,BW1,~]=MakePTH07(PSR1,repmat(N1, size(RAW(SesIndex(y)).Nrast{j},1),1),{2,0,BinSize});%-----DP used here
                             [PTH0,~,~]=MakePTH07(PSR0,repmat(N0, size(RAW(SesIndex(y)).Nrast{j},1),1),{1,0,BinSize});%BW1 reinjected here to make sure PTH0 & PTH1 have the same BW
                            % calculate MeanZ from PSTH before smoothing
                            if sum(PTH0,2)~=0
                                unsmoothPSTHz(1,1:length(Tm))=normalize(PTH1,PTH0,0);
                                R.Ses(i).Ev(k).rawMeanz(NN,1)=nanmean(unsmoothPSTHz(1,Tm>EventWin(k,1) & Tm<EventWin(k,2)),2);
                                R.Ses(i).Ev(k).rawMeanzPre(NN,1)=nanmean(unsmoothPSTHz(1,Tm>PreEventWin(k,1) & Tm<PreEventWin(k,2)),2);
                            else
                                R.Ses(i).Ev(k).rawMeanz(NN,1)=NaN;
                                R.Ses(i).Ev(k).rawMeanzPre(NN,1)=NaN;
                            end
                            
                            PTH1=smooth(PTH1,SmoothSPAN,SmoothTYPE)';
                            PTH0=smooth(PTH0,SmoothSPAN,SmoothTYPE)';
                            
                            %------------- Fills the R.Ses(i).Ev(k) fields --------------
                            R.Ses(i).Ev(k).BW(NN,1)=BW1;
                            R.Ses(i).Ev(k).PSTHraw(NN,1:length(Tm))=PTH1;
                            R.Ses(i).Ev(k).PSTHrawBL(NN,1:length(PTH0))=PTH0;
                            R.Ses(i).Ev(k).Meanraw(NN,1)=nanmean(R.Ses(i).Ev(k).PSTHraw(NN,Tm>PostEventWin(k,1) & Tm<PostEventWin(k,2)),2);
                            if sum(PTH0,2)~=0
                                R.Ses(i).Ev(k).PSTHz(NN,1:length(Tm))=normalize(PTH1,PTH0,0);
                                R.Ses(i).Ev(k).Meanz(NN,1)=nanmean(R.Ses(i).Ev(k).PSTHz(NN,Tm>EventWin(k,1) & Tm<EventWin(k,2)),2);
                                R.Ses(i).Ev(k).MeanzPRE(NN,1)=nanmean(R.Ses(i).Ev(k).PSTHz(NN,Tm>PreEventWin(k,1) & Tm<PreEventWin(k,2)),2);
                            else
                                R.Ses(i).Ev(k).PSTHz(NN,1:length(Tm))=NaN(1,length(Tm));
                                R.Ses(i).Ev(k).Meanz(NN,1)=NaN;
                                R.Ses(i).Ev(k).MeanzPRE(NN,1)=NaN;
                            end
                            
                            %------------------ firing (in Hz) per trial in pre/post windows ------------------
                            %used to make the between events comparisons and Response detection in a single window----
                            ev(k).pre=NaN(size(RAW(SesIndex(y)).Erast{EvInd},1),1);
                            ev(k).PreEvent=NaN(size(RAW(SesIndex(y)).Erast{EvInd},1),1);
                            ev(k).PostEvent=NaN(size(RAW(SesIndex(y)).Erast{EvInd},1),1);
                            for m=1:size(RAW(SesIndex(y)).Erast{EvInd},1) %loops through trials
                                ev(k).pre(m)=sum(PSR2{m}<prewin(k,2) & PSR2{m}>prewin(k,1))/(prewin(k,2)-prewin(k,1));
                                ev(k).PreEvent(m)=sum(PSR2{m}<PreEventWin(k,2) & PSR2{m}>PreEventWin(k,1))/(PreEventWin(k,2)-PreEventWin(k,1));
                                ev(k).PostEvent(m)=sum(PSR2{m}<PostEventWin(k,2) & PSR2{m}>PostEventWin(k,1))/(PostEventWin(k,2)-PostEventWin(k,1));
                            end
                            
                            %---------------------------Response detection w/ SignRank on pre/post windows---------------------------
                            
                            if  ~isempty(EvInd) && R.Ses(i).Ev(k).NumberTrials(NN,1)>= MinNumTrials && nanmean(ev(k).pre(:,1),1)>0.5 %avoid analyzing sessions where that do not have enough trials
                                
                                [~,R.Ses(i).Ev(k).ttestPreEvent(NN,1)]=ttest(ev(k).pre, ev(k).PreEvent); %Signrank used here because it is a dependant sample test
                                if R.Ses(i).Ev(k).ttestPreEvent(NN,1)<PStat
                                    R.Ses(i).Ev(k).RespDirPre(NN,1)=sign(mean(ev(k).PreEvent)-mean(ev(k).pre));
                                else R.Ses(i).Ev(k).RespDirPre(NN,1)=0;
                                end
                                
                                [~,R.Ses(i).Ev(k).ttestPostEvent(NN,1)]=ttest(ev(k).pre, ev(k).PostEvent); %Signrank used here because it is a dependant sample test
                                if R.Ses(i).Ev(k).ttestPostEvent(NN,1)<PStat
                                    R.Ses(i).Ev(k).RespDirPost(NN,1)=sign(mean(ev(k).PostEvent)-mean(ev(k).pre));
                                else R.Ses(i).Ev(k).RespDirPost(NN,1)=0;
                                end
                                
                                
                                if R.Ses(i).Ev(k).RespDirPre(NN,1)>0 || R.Ses(i).Ev(k).RespDirPost(NN,1) >0
                                    Search=find(Tm>=PreEventWin(k,1) & Tm<=PostEventWin(k,2));
                                    [R.Ses(i).Ev(k).MaxVal(NN,1),MaxInd]=max(R.Ses(i).Ev(k).PSTHraw(NN,Search));
                                    R.Ses(i).Ev(k).MaxTime(NN,1)=Tm(Search(1)+MaxInd-1);
                                else
                                    Search=find(Tm>=PreEventWin(k,1) & Tm<PostEventWin(k,2));
                                    [R.Ses(i).Ev(k).MaxVal(NN,1),MaxInd]=min(R.Ses(i).Ev(k).PSTHraw(NN,Search));
                                    R.Ses(i).Ev(k).MaxTime(NN,1)=Tm(Search(1)+MaxInd-1);
                                end
                            else R.Ses(i).Ev(k).RespDirEvent(NN,1)=0;
                            end
  
                        end %if ~isempty(PSR0{1}) || ~isempty(PSR1{1})
                    end %if EvInd=0 OR n(trials) < MinNumTrials fills with NaN
                end %Events
                
                fprintf('Neuron ID # %d\n',NN);
            elseif R.Ses(i).Structure(NN)~=0
            end %exclusion: IF R.Structure(NN)~=0 to avoid analyzing excluded neurons
        end %neurons: FOR j= 1:size(RAW(SesIndex(y)).Nrast,1)
    end %SesIndex
end %sessions: FOR i=1:length(RAW)

save('RearlyTraining.mat','R');


toc