function varargout=dopDataDropCalc(LRdata,filters,varargin)
% dopOSCCI Mark II: dopDataDropCalc ~ 1-06-2011 (last edit)
%
% calculate raw data drop out
%

% -------------------------------------------------------------------------
%     '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 Raw data dropout:');
try
    inputs.varargin=varargin;
    inputs.defaults=struct(...);
        'screen',[],...);% epoch screen column
        'sampleRate',25,...); % sample rate in hertz to figure out timing
        'filterZero',0,...); % event marker point/position within epoch
        'evType',2,... % event type: 0 = base, 1 = stim only, 2 = both
        'baseLRdata',[]); % for evType stimOnly see dopTaskDef - base data pre-filtered for baseline
    %         'GoR',2); % Goodness of Recording criterion
    ddc=setGetInputsStruct(inputs);
    
    ddc.translation='data index calculation';
    ddc.structAbbreviations='ddc';
    ddc.LRdata=LRdata; clear LRdata
    ddc.filter=filters; clear filters
    
    ddc.data=cat(3,ddc.LRdata,diff(flipdim(ddc.LRdata,3),1,3),mean(ddc.LRdata,3)); % L:R:difference data
    
    % epoch x sample x 4 array
    ddc.nEpochs=size(ddc.data,1);
    ddc.out.epoch.n=ddc.nEpochs;
    %% epoch screening?
    ddc.scrn=ones(ddc.nEpochs,1)==1; % use all by default
    if ~isempty(ddc.screen)
        if length(ddc.screen)==ddc.nEpochs
            ddc.scrn=logical(ddc.screen);
        else
            fprintf('\n%s\n\t%s\n\t%s\n',...
                'Mismatch in Screen filter length:',...
                [num2str(ddc.nEpochs),' Epochs, screen size = ',...
                num2str(size(ddc.screen))],...
                ' > using all epochs (no screening).');
        end
    end
    ddc.data=ddc.data(ddc.scrn,:,:); % update data, removing screened epochs
    ddc.nEpochs=size(ddc.data,1);
    ddc.out.epoch.n=ddc.nEpochs;
    %% for each filter
    ddc.fns=fieldnames(ddc.filter); % filter names
    for i=1:length(ddc.fns)
        ddc.fn=ddc.fns{i};
        switch ddc.fn
            case {'zeroPt','zeroValue'} % don't do anything for this
                
            case  {'left','right','avg','dif'} % index filters
                ddc.tmp.fnsj=fieldnames(ddc.filter.(ddc.fn));
                for j=1:numel(ddc.tmp.fnsj);ddc.tmp.fnj=ddc.tmp.fnsj{j};
                    ddc.tmp.fnsjj=fieldnames(ddc.filter.(ddc.fn).(ddc.tmp.fnj));
                    for jj=1:numel(ddc.tmp.fnsjj);ddc.tmp.fnjj=ddc.tmp.fnsjj{jj};
                        ddc.tmp.vn=[ddc.fn,ddc.tmp.fnj]; % variable name
                        [ddc.out.(ddc.tmp.vn) ddc.out.all.(ddc.tmp.vn)]=...
                            dropCalc(ddc.data(:,ddc.filter.(ddc.fn).(ddc.tmp.fnj).(ddc.tmp.fnjj).filter,:));
                        ddc.out.(ddc.tmp.vn).time=ddc.filter.(ddc.fn).(ddc.tmp.fnj).(ddc.tmp.fnjj).sec;
                        ddc.out.all.(ddc.tmp.vn).time=ddc.out.(ddc.tmp.vn).time;
                    end
                end
            otherwise % calculate inddces
                %% > mean, std, min, max
                [ddc.out.(ddc.fn) ddc.out.all.(ddc.fn)]=dropCalc(ddc.data(:,ddc.filter.(ddc.fn),:));
                %                 ddc.out.(ddc.fn).mean=mean(ddc.data(:,ddc.filter.(ddc.fn),:),2);
                %                 ddc.out.(ddc.fn).std=std(ddc.data(:,ddc.filter.(ddc.fn),:),1,2);
                %                 ddc.out.(ddc.fn).min=min(ddc.data(:,ddc.filter.(ddc.fn),:),[],2);
                %                 ddc.out.(ddc.fn).max=max(ddc.data(:,ddc.filter.(ddc.fn),:),[],2);
                %                 ddc.out.(ddc.fn).rms=sqrt(mean(ddc.data(:,ddc.filter.(ddc.fn),:).^2,2)...
                %                     /size(ddc.data(:,ddc.filter.(ddc.fn),:),2));
                
        end
    end
    %% > Stim Only Event type:
    % for stim only event type, input baseline data
    if ddc.evType==1 && ~isempty(ddc.baseLRdata)
        ddc.baseData=cat(3,ddc.baseLRdata,ddc.baseLRdata(:,:,1)-ddc.baseLRdata(:,:,2));%diff(ddc.baseLRdata,1,3));
        ddc.baseData=ddc.baseData(logical(ddc.scrn),:,:); % update data, removing screened epochs
        ddc.fn='base';
        [ddc.out.(ddc.fn) ddc.out.all.(ddc.fn)]=dropCalc(ddc.baseData(:,:,:));
        %         ddc.out.(ddc.fn).mean=mean(ddc.baseData(:,:,:),2);
        %         ddc.out.(ddc.fn).std=std(ddc.baseData(:,:,:),1,2);
        %         ddc.out.(ddc.fn).min=min(ddc.baseData(:,:,:),[],2);
        %         ddc.out.(ddc.fn).max=max(ddc.baseData(:,:,:),[],2);
        %         ddc.out.(ddc.fn).rms=sqrt(mean(ddc.baseData(:,:,:).^2,2)...
        %             /size(ddc.baseData(:,:,:),2));
    end
    
    %% index description
    % include a description of the inddces
    ddc.out.note=sprintf(['%s\n\t%s\n\t\t%s\n\t\t%s\n\t\t%s\n'...% divisions
        '\t%s\n\t\t%s\n\t\t%s\n\t\t%s\n\t\t%s\n\t\t%s\n\t\t%s\n\t\t%s\n'],...
        'Baseline corrected data indices for left,right, and difference data',...
        'includes divisions: ','epoch (all data)', 'base (baseline), and ',...
        'poi (period of interest)',...
        'summaries:','mean', 'std (standard deviation)', 'min', 'max', ...
        'rms (root mean square)',...
        'liv (laterality index variability from Vadikolias et al., 2007)',...
        'gor (goodness of recording from Knecht et al., 2001)');
    
    %% set ouputs
    varargout{1}=ddc.out;
    fprintf('\t%s\n','Done');
catch err
    %% catch dopOSCCI error
    save(dopOSCCIdebug(mfilename));rethrow(err);
end

clear ddc inputs
end
%%-------------------
%% embedded function
%% > inddces
% date in, inddces calculations out:
function [ep all]=dropCalc(data)
ep=[]; % epoch - number for each epoch (left, right, diff)
ep.n=size(data,1);
ep.count = sum(data==0,2);
ep.pct = sum(data==0,2)/size(data(:,:,:),2)*100; % percentage dropout

% ep.mean=mean(data(:,:,:),2);
% ep.std=std(data(:,:,:),1,2);
% ep.sem=std(data(:,:,:),1,2)/sqrt(ep.n);
% ep.CI95=ep.sem*1.96;
% ep.min=min(data(:,:,:),[],2);
% ep.max=max(data(:,:,:),[],2);
% ep.rms=sqrt(mean(data(:,:,:).^2,2)...
%     /size(data(:,:,:),2));
% ep.med=median(data(:,:,:),2);
% ep.iqr=iqr(data(:,:,:),2);
%
% % gor = goodness of recording (see Knecht et al., 2001, Brain)
% ep.gorMean = ep.rms(:,:,:);
% ep.gorMax = ep.rms(:,:,:);


all=[]; % all epochs (left, right, diff, average)
all.n=ep.n;
if all.n==1;
    all.count = ep.count;
    all.pct = ep.pct;
    %     all.mean=ep.mean;
    %     all.std=ep.std;
    %     all.sem=ep.sem;
    %     all.CI95=all.sem*1.96;
    %     all.min=ep.min;
    %     all.max=ep.max;
    %     all.rms=ep.rms;
    %     all.med=ep.med;
    %     all.iqr=ep.iqr;
    %     all.gorMean = ep.rms;
    %     all.gorMax = ep.rms;
else
    all.count = sum(ep.count); % total number of samples
    all.pct = sum(ep.pct); % total of baseline affected - NOT average
    %     all.mean=mean(ep.mean);
    %     all.std=mean(ep.std);
    %     all.sem=mean(ep.sem);
    %     all.CI95=all.sem*1.96;
    %     all.min=min(ep.min);
    %     all.max=max(ep.max);
    %     all.rms=mean(ep.rms);
    %     all.med=mean(ep.med);
    %     all.iqr=mean(ep.iqr);
    %     all.gorMean = mean(ep.rms);
    %     all.gorMax = max(ep.rms);
end

clear data;
end
