function varargout=dopDataLIcalc(LRdata,poiFilter,varargin)
% dopOSCCI Mark II: dopDataLIcalc ~ 06-Dec-2011 (last edit)
%
% calculate Laterality Index (etc.) for each event channel within period of
% interest
%
% expects difData to be in form X (epoch) by Y (samples) and filter to be
% single row of Y (samples) size logical array.
%
% want to do a couple of things:
% - calculate overall LI by finding peak difference in POI and taking
% average dif around 2 seconds (actWin) of thins point).
% - split this up for overall and by trial
%

% -------------------------------------------------------------------------
%     'dopOSCCI' summarises functional transcranial Doppler ultrasonography
%     (fTCD) data.  The primary function of the software is to assess the
%     hemispheric lateralization of cognitive function.
%
%     Copyright (C) 2011 the Chancellor, Masters and Scholars of the
%     University of Oxford.
% 
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
% 
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/> 
%     or write to the Free Software Foundation, Inc., 51 Franklin Street,
%     Fifth Floor, Boston, MA 02110-1301, USA.
%
%     Authors: Nic Badcock, Georgina Holt, Anneka Holden, & Dorothy Bishop.
%       contact: nicholas.badcock@gmail.com
%
%     Type 'dopOSCCIlicense' for license details or see dopOSCCIlicense.txt
%

fprintf('\n\t%s\n\t\t%s',['Running ',mfilename],...
    'Calculating Event Laterality Index:');
lic=[];
try
    inputs.varargin=varargin;
    inputs.defaults=struct('actWin',2,...);%,... % activation window duration in seconds
        'screen',[],...);% epoch screen column
        'sampleRate',25,...); % sample rate in hertz to figure out timing
        'filterZero',0); % event marker point/position within epoch
    lic=setGetInputsStruct(inputs);
    
    lic.translation='laterality index calculation';
    lic.structAbbreviations='lic';
    lic.LRdata=LRdata; clear LRdata
    lic.filter=poiFilter; clear poiFilter
    
    lic.data=diff(flipdim(lic.LRdata,3),1,3); % difference data
    
    lic.actWinSamples=lic.actWin/(lic.sampleRate^-1); % seconds to samples
    
    lic.nEpochs=size(lic.data,1);
    %% epoch screening?
    lic.scrn=ones(lic.nEpochs,1)==1; % use all by default
    if ~isempty(lic.screen)
        if length(lic.screen)==lic.nEpochs
            lic.scrn=logical(lic.screen);
        else
            fprintf('\n%s\n\t%s\n\t%s\n',...
                'Mismatch in Screen filter length:',...
                [num2str(lic.nEpochs),' Epochs, screen size = ',...
                num2str(size(lic.screen))],...
                ' > using all epochs (no screening).');
        end
    end

    
    %% dummy data
    % in case something is missing
    lic.tmp.allList={'LImean','LIstd','LIsem','LI95pctCI','LImin','LImax',...
            'LImedian','LIiqr','LIepochs','LIpeakSample','LIpeakSec',...
            'LIt1p','t1t','t1df','t1ci','t1stats','LIt1cohen','LIt1sig',...
            'LIt2p','t2t','t2df','t2ci','t2stats','LIt2cohen','LIt2sig',...
            'LIrp','LIrsig','rstats','LIrzval','LIrZr','LIepochs'};
        for i=1:length(lic.tmp.allList)
            if strfind(lic.tmp.allList{i},'ci') %&& ~strfind(lic.tmp.allList{i},'95pct')
                lic.out.all.(lic.tmp.allList{i})=[999 999];
            elseif strfind(lic.tmp.allList{i},'sig')
                lic.out.all.(lic.tmp.allList{i})=0;
            elseif strfind(lic.tmp.allList{i},'stats')
                if strcmp(lic.tmp.allList{i},'rstats')
                    lic.out.all.(lic.tmp.allList{i}).zval=999;
                    lic.out.all.(lic.tmp.allList{i}).ranksum=999;
                else
                    lic.out.all.(lic.tmp.allList{i}).tstat=999;
                    lic.out.all.(lic.tmp.allList{i}).df=999;
                    lic.out.all.(lic.tmp.allList{i}).sd=999;
                end
            else
                lic.out.all.(lic.tmp.allList{i})=999;
            end
        end
        
    lic.data=lic.data(lic.scrn,:); % update data, removing screened epochs
    if sum(isnan(lic.data))
        fprintf('\t%s\n','NaN found in data - removing row');
        lic.tmp.scrn=sum(isnan(lic.data(:,:)),2);
        lic.data(lic.tmp.scrn~=0,:)=[];
    end
    lic.nEpochs=size(lic.data,1);
    if lic.nEpochs %sum(lic.scrn)
        %% Overall: 'all'
        %% >find peak in poi
        lic.meanData=mean(lic.data);
        if size(lic.data,1)==1
            lic.meanData=lic.data;
        end
        lic.poiPeakSample=find(abs(lic.meanData)==max(abs(lic.meanData(:,lic.filter))),1,'first');
        
        % make sure it's within limits - sometimes multiple values are the same
        % and are outside period of interest
        if lic.poiPeakSample<find(lic.filter,1,'first');
            lic.tmp.pps=find(abs(lic.meanData)==max(abs(lic.meanData(:,lic.filter))));
            j=2;
            lic.poiPeakSample=lic.tmp.pps(j);
            while and(lic.poiPeakSample<find(lic.filter,1,'first'),j<length(lic.tmp.pps))
                j=j+1; lic.poiPeakSample=lic.tmp.pps(j);
            end
        elseif lic.poiPeakSample>find(lic.filter,1,'last');
            lic.tmp.pps=find(abs(lic.data)==max(abs(lic.meanData(:,lic.filter))));
            j=length(lic.tmp.pps)-1;
            lic.poiPeakSample=lic.tmp.pps(j);
            while and(lic.poiPeakSample<find(lic.filter,1,'first'),j>1)
                j=j-1; lic.poiPeakSample=lic.tmp.pps(j);
            end
        end
        lic.poiFilter=zeros(1,size(lic.data,2));
        lic.actWinValues=[lic.poiPeakSample-lic.actWinSamples*.5...
            lic.poiPeakSample+lic.actWinSamples*.5]; % lower and upper activation window values
        
        % if value is outside of available data range, move it back
        if lic.actWinValues(2)>size(lic.data,2)
            lic.actWinValues=lic.actWinValues-(lic.actWinValues(2)-size(lic.data,2));
        elseif lic.actWinValues(1)<1
            lic.actWinValues=lic.actWinValues+abs(lic.actWinValues(1))+1;
        end
        lic.poiFilter(lic.actWinValues(1):lic.actWinValues(2))=1; % filter
        lic.poiFilter=logical(lic.poiFilter);
        % include a description of the indices
        lic.out.all.note=sprintf(['%s\n\t%s\n\t%s\n\t%s\n\t%s\n\t%s\n'...
            '\t%s\n\t\t%s\n\t\t%s\n\t\t%s\n'... % LIt1*
            '\t%s\n\t\t%s\n\t\t%s\n\t\t%s\n'... % LIt2*
            '\t%s\n\t\t%s\n\t\t%s\n'],... % LIr*
            'all = all included epochs',...
            'LImean/median = laterality index (average of peak differences across epochs)',...
            'LIsem = standard error of the mean for Laterality Index',...
            'LIstd/iqr = standard deviation/inter-quartile range of LIs',...
            'LIpeakSample = sample number for peak difference within epoch',...
            'LIpeakSec = sample time for peak difference within epoch',...
            'LIt1* = one-sample t-test for LI dif to zero ',...
            '(sig = y/n, p = p-value, ci = 95% confidence interval, ',...
            'stats = t-value, degrees of freedom, & standard deviation, ',...
            'cohen = Cohen''s d effect size)',...
            'LIt2* = paired-sample t-test for LI dif between Left and Right ',...
            '(sig = y/n, p = p-value, ci = 95% confidence interval, ',...
            'stats = t-value, degrees of freedom, & standard deviation, ',...
            'cohen = Cohen''s d effect size)',...
            'LIr* = Wilcoxon ranksum test of dif between Left and Right ',...
            ['(p = p-value, sig = y/n, ',...
            'stats = ranksum value and z-value, '],...
            'Zr = Z to r conversion, effect size)');
        
        % mean of differences for each epoch
        lic.out.all.LIs=mean(lic.data(:,lic.poiFilter),2);
        
        % > left and right for stats tests
        lic.out.all.LIsLeft=mean(lic.LRdata(lic.scrn,lic.poiFilter,1),2);
        lic.out.all.LIsRight=mean(lic.LRdata(lic.scrn,lic.poiFilter,2),2);
        
        % parametric descriptives
        lic.out.all.LImean=mean(lic.out.all.LIs);
        lic.out.all.LIstd=std(lic.out.all.LIs);
        lic.out.all.LIsem=lic.out.all.LIstd/sqrt(lic.nEpochs);
        lic.out.all.LI95pctCI=lic.out.all.LIsem*1.96;
        
        lic.out.all.LImin=min(lic.out.all.LIs);
        lic.out.all.LImax=max(lic.out.all.LIs);
        
        % non-parametric descriptives
        lic.out.all.LImedian=median(lic.out.all.LIs);
        try
            lic.out.all.LIiqr=iqr(lic.out.all.LIs); % inter-quartile range
            % t-test:
        % > one-sample
        [lic.out.all.LIt1sig lic.out.all.LIt1p ...
            lic.out.all.t1ci lic.out.all.t1stats]=ttest(lic.out.all.LIs,0);
        lic.out.all.t1ci = abs(diff([mean(lic.out.all.t1ci) lic.out.all.t1ci(1)]));
        
        lic.out.all.t1t=lic.out.all.t1stats.tstat;
        lic.out.all.t1df=lic.out.all.t1stats.df;
        lic.out.all.t1sd=lic.out.all.t1stats.sd;
               % >> One-sample Kolmogorov-Smirnov test - normality
        % rejecting even though looks good - some of the other tests don't
        %         [lic.out.all.ks1_norm lic.out.all.ks1_p lic.out.all.ks1_stat ...
        %             lic.out.all.ks1_cv]=kstest(lic.out.all.LIs);
        %         [lic.out.all.lt1_norm lic.out.all.lt1_p lic.out.all.lt1_stat ...
        %             lic.out.all.lt1_crit]=lillietest(lic.out.all.LIs);
        
        % > two-sample: dependent
        [lic.out.all.LIt2sig lic.out.all.LIt2p ...
            lic.out.all.t2ci lic.out.all.t2stats]=ttest(lic.out.all.LIsLeft,lic.out.all.LIsRight);
        lic.out.all.t2ci = abs(diff([mean(lic.out.all.t2ci) lic.out.all.t2ci(1)]));
        
        lic.out.all.t2t=lic.out.all.t2stats.tstat;
        lic.out.all.t2df=lic.out.all.t2stats.df;
        lic.out.all.t2sd=lic.out.all.t2stats.sd;
    if numel(lic.out.all.LIs)>3 
            [lic.out.all.jb1_norm lic.out.all.jb1_p lic.out.all.jb1_stat ...
                lic.out.all.jb1_crit]=jbtest(lic.out.all.LIs);
    end
    if numel(lic.out.all.LIs)>3 
    [lic.out.all.jb2_Lnorm lic.out.all.jb2_Lp lic.out.all.jb2_Lstat ...
        lic.out.all.jb2_Lcrit]=jbtest(lic.out.all.LIsLeft);
    [lic.out.all.jb2_Rnorm lic.out.all.jb2_Rp lic.out.all.jb2_Rstat ...
        lic.out.all.jb2_Rcrit]=jbtest(lic.out.all.LIsRight);
    end
        % Mann-Whitney ranksum test
        [lic.out.all.LIrp,lic.out.all.LIrsig,lic.out.all.rstats] = ...
            ranksum(lic.out.all.LIsLeft,lic.out.all.LIsRight);
        lic.out.all.LIranksum=lic.out.all.rstats.ranksum;
        

        % > Z to effect size = r value
        lic.out.all.LIrZr=999; % dummy value
        lic.out.all.rstats.zval=999;
        if isfield(lic.out.all.rstats,'zval')
            lic.out.all.LIrzval=lic.out.all.rstats.zval;
            lic.out.all.LIrZr=lic.out.all.rstats.zval/sqrt(lic.out.all.LIepochs*2);
        end
        catch err
            disp(err)
%             lic.out.all.LIiqr=999;
        end
        lic.out.all.LIepochs=lic.nEpochs;
        lic.out.all.LIpeakSample=(lic.poiPeakSample-lic.filterZero);
        lic.out.all.LIpeakSec=(lic.poiPeakSample-lic.filterZero)*(lic.sampleRate^-1);
        if numel(lic.out.all.LIs)>1        
        
        
     
        
        % > cohen's d for one-sampled t-test
        lic.out.all.LIt1cohen=abs(lic.out.all.LImean)/std(lic.out.all.LIs);
        
        
        
        % > cohen's d for one-sampled t-test
        lic.out.all.LIt2cohen=...
            (mean(lic.out.all.LIsLeft)-mean(lic.out.all.LIsRight))/...
            sqrt((((lic.out.all.LIepochs-1)*std(lic.out.all.LIsLeft).^2)+...
            ((lic.out.all.LIepochs-1)*std(lic.out.all.LIsRight).^2))/...
            lic.out.all.LIepochs*2);
        
        % > normality L & R
%         [lic.out.all.ks2_Lnorm lic.out.all.ks2_Lp lic.out.all.ks2_Lstat ...
%             lic.out.all.ks2_Lcv]=kstest(lic.out.all.LIsLeft);
%         [lic.out.all.ks2_Rnorm lic.out.all.ks2_Rp lic.out.all.ks2_Rstat ...
%             lic.out.all.ks2_Rcv]=kstest(lic.out.all.LIsRight);
    
        end
        %     %% saving categories
        %     % define categories to save data with
        %     lic.out.all.sav.saveDescrip={'LIepochs','LIpeakSec','LImean','LIstd','LIsem',...
        %         'LImin','LImax','LImedian','LIiqr'};
        %     lic.out.all.sav.save1sample={'LIt1sig','LIt1cohen'};
        %     lic.out.all.sav.save2sample={'LIt2sig','LIt2cohen'};
        %     lic.out.all.sav.saveRank={'LIrsig','LIrZr'};
    else
        %% no acceptable trials...
        fprintf('\t%s\n','No suitable trials available');
        lic.out.all.note='no acceptable epochs';
    end
    %% by Epoch: 'epoch'
    lic.out.epoch.note='epoch = single epoch single (e.g., by trial)';
    lic.ep.trans='summaries for each Epoch';
    
    [lic.ep.mxs lic.ep.poiPeakSample]=max(abs(lic.data(:,lic.filter)),[],2); % max values for each epoch
    lic.ep.poiPeakSample=lic.ep.poiPeakSample+find(lic.filter,1,'first'); % adjust for filter
    %     for i=1:lic.nEpochs;
    %         lic.ep.poiPeakSample(i)=find(abs(lic.data(i,:))==lic.ep.mxs(i),1,'first');
    
    %         lic.ep.poiFilter=zeros(1,size(lic.data,2));
    
    lic.ep.actWinValues=[lic.ep.poiPeakSample-lic.actWinSamples*.5...
        lic.ep.poiPeakSample+lic.actWinSamples*.5]; % lower and upper activation window values
    
    % if value is over available data range, move down
    lic.ep.grFilt=lic.ep.actWinValues(:,2)>size(lic.data,2);
    lic.ep.actWinValues(lic.ep.grFilt,:)=...
        bsxfun(@minus,lic.ep.actWinValues(lic.ep.grFilt,:),...
        lic.ep.actWinValues(lic.ep.grFilt,2)-size(lic.data,2));
    % if value is under available data range, move it up
    lic.ep.lsFilt=lic.ep.actWinValues(:,1)<1; % less filter
    lic.ep.actWinValues(lic.ep.lsFilt,:)=...
        bsxfun(@plus,lic.ep.actWinValues(lic.ep.lsFilt,:),...
        abs(lic.ep.actWinValues(lic.ep.lsFilt,1))+1);
    
    %         lic.ep.poiFilter(lic.ep.actWinValues(1):lic.ep.actWinValues(2))=1; % filter
    
    % create filter full rows of the same numbers 1:number of columns
    lic.ep.poiFilter=bsxfun(@times,ones(size(lic.data)),1:size(lic.data,2));
    
    lic.ep.poiFilter(bsxfun(@lt,lic.ep.poiFilter,lic.ep.actWinValues(:,1)))=0;
    lic.ep.poiFilter(bsxfun(@gt,lic.ep.poiFilter,lic.ep.actWinValues(:,2)))=0;
    lic.ep.poiFilter(lic.ep.poiFilter>0)=1;
    
    lic.out.epoch.n=lic.nEpochs;
    for i=1:lic.nEpochs;
        lic.out.epoch.LImean(i)=mean(lic.data(i,logical(lic.ep.poiFilter(i,:))));
        lic.out.epoch.LIstd(i)=std(lic.data(i,logical(lic.ep.poiFilter(i,:))));
        lic.out.epoch.LIsem(i)=lic.out.epoch.LIstd(i);
        lic.out.epoch.LI95pctCI(i)=lic.out.epoch.LIsem(i)*1.96;
        lic.out.epoch.LImin(i)=min(lic.data(i,logical(lic.ep.poiFilter(i,:))));
        lic.out.epoch.LImax(i)=max(lic.data(i,logical(lic.ep.poiFilter(i,:))));
        lic.out.epoch.LImed(i)=median(lic.data(i,logical(lic.ep.poiFilter(i,:))));
        try
        lic.out.epoch.LIiqr(i)=iqr(lic.data(i,logical(lic.ep.poiFilter(i,:))));
        catch
            lic.out.epoch.LIiqr(i)=999;
        end
    end
    %         lic.out.epoch.LIepochs(i)=1;
    lic.out.epoch.LIpeakSample=(lic.ep.poiPeakSample-lic.filterZero);
    lic.out.epoch.LIpeakSec=(lic.ep.poiPeakSample-lic.filterZero)*(lic.sampleRate^-1);
    %         lic.out.epoch.LIsem(i)=std(mean(lic.data(i,logical(lic.ep.poiFilter))))/sqrt(1);
    
    %% report
    %     fprintf('X is %4.2f meters or %8.3f mm\n', 9.9, 9900, B)
    lic.tmp.sig={'not ',''};
    
    if ~sum(lic.out.all.LIt1sig==[1 0]); lic.out.all.LIt1sig=0; end
    if ~sum(lic.out.all.LIt2sig==[1 0]); lic.out.all.LIt2sig=0; end
    if ~sum(lic.out.all.LIrsig==[1 0]); lic.out.all.LIrsig=0; end
    
    lic.out.report=sprintf('%s\n%s\n\t%s\n\t%s\n%s\n\t%s\n',...
        sprintf('Number of epochs = %u',lic.out.all.LIepochs),...
        sprintf('LI mean = %3.2f (sd %3.2f), sem=%3.2f, min=%3.2f, max=%3.2f.',...
        lic.out.all.LImean,lic.out.all.LIstd,lic.out.all.LIsem,lic.out.all.LImin,lic.out.all.LImax),...
        sprintf([lic.tmp.sig{lic.out.all.LIt1sig+1},'dif to Zero: t',...
        '(%u)=%3.2f, p=%3.2f, Cohen''s d=%3.2f'],...
        lic.out.all.t1stats.df,lic.out.all.t1stats.tstat,lic.out.all.LIt1p,lic.out.all.LIt1cohen),...
        sprintf(['Left ',lic.tmp.sig{lic.out.all.LIt2sig+1},'dif Right: t',...
        '(%u)=%3.2f, p=%3.2f, Cohen''s d=%3.2f'],...
        lic.out.all.t2stats.df,lic.out.all.t2stats.tstat,lic.out.all.LIt2p,lic.out.all.LIt2cohen),...
        sprintf('LI median = %3.2f (inter-quartile range %3.2f).',...
        lic.out.all.LImedian,lic.out.all.LIiqr),...
        sprintf(['Left ',lic.tmp.sig{lic.out.all.LIrsig+1},'dif Right: U',...
        '=%3.2f, p=%3.2f, r=%3.2f.'],...
        lic.out.all.rstats.ranksum,lic.out.all.LIrp,lic.out.all.LIrZr));
    
    %% set ouputs
    varargout{1}=lic.out;
    fprintf('\t%s\n','Done');
catch err
    save(dopOSCCIdebug(mfilename));rethrow(err);
end

clear lic inputs


