VentralPallidumAlcohol2019/AnalysisAndFigures/SupportingPrograms/ResDetect08.m

function  [DATout]=ResDetect08(Rin,Rbasein,INTin,NORMTYPin,WINin,CI) 
% This function detects responses in peri-event histograms. Confidence intervals are extracted from the
% distribution of the baseline period.
% Arpil 2012, Fred: This version was adapted from RespDetect07 (now handles single-neuron PTH)
% This version contains the helper functions "sigtestin", "normalizein", "findonset03in" 
% Inputs:
%       Rin       is a matrix e.g. PTH1(neurons, bins)
%       Rbasein   is the baseline version of Rin (used to normalize)
%       INTin     is the window in which responses are searched
%       NORMTYPin for the function normalizein - should be set to zero
%       WINin     are onset and offset requirements; by default it is [0.1 0.3]
%       CI        is the confidence interval (inhibitions<lower bound, excitation>Upper bound) 
% Outputs:
%       DATout.R= Rin
%             .LUb(neurons, 1:2): Lower and Upper Bounds = Threshold ofr excitation and inhibitions
%             .Rnorm(neurons, bins)= z-score transformed Rin
%             .Rnormmask(neurons, bins) = z-score with all bins that are within the confidence intervals are zeroed
%             .RFLAG(neurons, 1): +1 for excitation, -1 for inhibition, 0 for no resp.
%             .CFLAG(neurons, 1)= 0 or 1. 1 means that an inhibition follows an excation or vice-versa 
%             .Onsets(neurons,1:2): First and Second columns stores the onset of excitations and inhibitions respectively
%             .Offsets(neurons,1:2): same as onset
%             .Paramnames={'Dura','Baseline','Bsize','INT','WINs'};
%             .Params={Dura,Baseline,BSIZE,INTin,WINin};


global Dura Tm Tbase BSIZE Erefnames 

%Allows to handle single neurones cases where R(bins, 1) instead of R(1,neurons)
if size(Rin,2)==1
    R=Rin'; Rbase=Rbasein';
else
    R=Rin;Rbase=Rbasein;
end

if isempty(WINin)
    WINON=ceil(0.1/BSIZE);WINOFF=ceil(0.3/BSIZE);
    WINFLAG=0;
elseif size(WINin,1)==1
    WINON=ceil(WINin(1)/BSIZE);WINOFF=ceil(WINin(2)/BSIZE); %same detection for all neurons
    WINFLAG=0;
else
    WINONmat=ceil(WINin(:,1)/BSIZE);WINOFFmat=ceil(WINin(:,2)/BSIZE); % different for each neuron
    WINFLAG=1;
end

if size(INTin,1)==1
    Iref=find(Tm>=INTin(1) & Tm<=INTin(2));
    INTFLAG=0;
else
    INTmat=INTin;
    INTFLAG=1;
end

% ------------ Determines the thresholds and makes a mask ------------
%ALPHA=0.05/(size(R,2)-length(Tmind)); %Bonferroni correction
if length(CI)==1, CI0(1)=CI;CI0(2)=CI; else CI0=CI; end %size(CI0)=(1x2)
if CI0(1)>0 
    %LUb= Lower and Upper Bound = threshold determined by the confidence intervals
    LUb=(prctile(Rbase',[CI0(1)/2 100-CI0(2)/2]))'; % Outside LUb range consider significant
    if size(LUb,2)==1, LUb=LUb'; end %To handle single neuron cases
else
    %Ali's stuff - I think it uses it to go beyond 0 & 100% confidence intervals
    %it can be reached by using negative CI values
    LUb(:,1)=min(Rbase,[],2)-(max(Rbase,[],2)-min(Rbase,[],2))*abs(CI0(1))/200;
    LUb(:,2)=max(Rbase,[],2)+(max(Rbase,[],2)-min(Rbase,[],2))*abs(CI0(2))/200;
end
[lmask,umask]=sigtestin(R,LUb);
Rnorm=normalizein(R,Rbase,NORMTYPin);%used to put (-) to bins below lower bound
Rnormmask=Rnorm.*(umask+lmask);

%Outupts R, LUb, Normmask and Rnorm
DATout.R=R;DATout.LUb=LUb;DATout.Rnormmask=Rnormmask;DATout.Rnorm=Rnorm;

for j=1:size(R,1) %loops through neurons
    if mod(j,11)==0
        fprintf('.');
    elseif mod(j,100)==0
        fprintf('%d',j)
    end
    if WINFLAG
        WINON=WINONmat(j);WINOFF=WINOFFmat(j);
    end
    if INTFLAG
        Iref=find(Tm>INTmat(j,1) & Tm<INTmat(j,2));
    end
    %PInd1 and NInd1 are the Tm indexes of the first excitation and inhibition respectively
    PInd1=findonset03in(double(Rnormmask(j,Iref)>0),1,WINON);
    NInd1=findonset03in(double(Rnormmask(j,Iref)<0),1,WINON);
    
    
    %% Flag settings
    %------------ Finds the onsets/offset of the FIRST response detected ------------
    Ponset=NaN;Nonset=NaN;Poffset=NaN;Noffset=NaN;TFLAG=[0 0];%termination flag
    if PInd1<NInd1
        % RFLAG= Response Flag: +1:excitation, -1:inhibition, 0:No response
        RFLAG=1;
        % Ponset is the excitation onset
        Ponset=Tm(Iref(PInd1));
    elseif PInd1>NInd1
        RFLAG=-1;
        % Nonset is the inhibition onset
        Nonset=Tm(Iref(NInd1));
    elseif isinf(PInd1) && isinf(NInd1)
        RFLAG=0;
    end
    
    % CFLAG= conflict flag=1 if both excitation and inhibition detected
    if ~isinf(PInd1) && ~isinf(NInd1)
        CFLAG=1;
        Ponset=Tm(Iref(PInd1));Nonset=Tm(Iref(NInd1));
    else
        CFLAG=0;
    end
    
    %------------ Finds the onsets/offset of the SECOND response detected ------------
    if ~isinf(PInd1)
        PInd2=findonset03in(double(Rnormmask(j,:)<=0),Iref(PInd1)+1,WINOFF);
        if ~isinf(PInd2)
            Poffset=Tm(PInd2);
            if Poffset<-0.03 % 0.025 is the interp1 offset with 0.05 bin
                TFLAG(1)=1;
            end
        else
            Poffset=Tm(end)+sqrt(-1); % if offset not detected it is beyond Tm(end)
            PInd2=length(Tm);
        end
    end
    if ~isinf(NInd1)
        NInd2=findonset03in(double(Rnormmask(j,:)>=0),Iref(NInd1)+1,WINOFF);
        if ~isinf(NInd2)
            Noffset=Tm(NInd2);
            if Noffset<-0.03 % 0.025 is the interp1 offset with 0.05 bin
                TFLAG(2)=1;
            end
        else
            Noffset=Tm(end)+sqrt(-1); % if offset not detected it is beyond Tm(end)
            NInd2=length(Tm);
        end
    end
    
    % setting outputs (new change to work for LP,Ent and Ext while keeping DS,NS results intact
    if (CFLAG==0 && sum(TFLAG)==1) || (CFLAG==1 && sum(TFLAG)==2)  % early termination response not valid (set as imaginary)
        DATout.RFLAG(j,1)=RFLAG*sqrt(-1);
    else
        DATout.RFLAG(j,1)=RFLAG;
    end
    DATout.RFLAG(j,1)=RFLAG;
    DATout.CFLAG(j,1)=CFLAG;
    % offset and onsets should be swaped and negated
    DATout.Onsets(j,:)=[Ponset,Nonset];DATout.Offsets(j,:)=[Poffset,Noffset];
end

% setting outputs
DATout.Paramnames={'Dura','Bsize','INT','WINs'};
DATout.Params={Dura,BSIZE,INTin,WINin};

%fprintf('\n response detection completed \n')


%----- HELPER FUNCTIONS ----

function [LMask,UMask]=sigtestin(Matin,Bin)
UMask=zeros(size(Matin));LMask=zeros(size(Matin));
for k=1:size(Matin,1)
    TMP=Matin(k,:);B=Bin(k,:); 
    LMask(k,find(TMP<(B(1)-eps)))=1;
    UMask(k,find(TMP>(B(2)+eps)))=1;
end

function Mnormal=normalizein(M,Mbase,MD)
%MD is the type of normalization. 
%   0 = Z-score  >>(X-mean)/SD
%   1 = normalized to response max
%   2 = Baseline substracted
%   3 = normalized by SD only >> X/SD

%Added April 2012 to handle vectors
if size(M,2)==1, M=M'; Mbase=Mbase'; end
Mnormal=zeros(size(M));
for k=1:size(M,1)
    if MD==0
        Mnormal(k,:)=(M(k,:)-nanmean(Mbase(k,:)))/nanstd(Mbase(k,:)); % Zscore
    elseif MD==1
        Tmp=M(k,:)-nanmean(Mbase(k,:));
        Mnormal(k,:)=Tmp/max(abs(Tmp));% mean corrected normalized to RESPONSE MAX
    elseif MD==2
       Mnormal(k,:)=(M(k,:)-nanmean(Mbase(k,:)));
    elseif MD==3
            Mnormal(k,:)=M(k,:)/nanstd(Mbase(k,:)); % Zscore no demean
    end
end


function Indout=findonset03in(Rin,Indin,Dur)
% Finds the index of the first element when 'Dur' consecutive elements are equal to 1
% Rin must contain double and not logical values. Therefore, the result of a find function
% must be converted as follow: double(find(MAT>X))
% Indin: looks from index Indin onward
% Dur: Indetfies onset only if response is Dur long
Indout=Inf;
if ~isinf(Indin)
    R=Rin(Indin:end);
    TMP=find(conv(ones(1,Dur),R)==Dur)-(Dur-1);
    if ~isempty(TMP)
        Indout=TMP(1)+Indin-1;
    end
end