function varargout=dopDataIndexCalc(LRdata,filters,varargin)
% dopOSCCI Mark II: dopDataIndexCalc ~ 10-11-2010 (last edit)
%
% calculate data indices/descriptives
%
% expects difData to be in form X (epoch) by Y (samples) and filter to be
% single row of Y (samples) size logical array.
%
% mean
% std
% min
% max
% rms = root mean square
%
% LI - variability:
%   Vadikolias et al. 2007, J Cereb Blood Flow Metab 27: 1870-1877.
%
% Goodness of recording
%   Knecht et al., (2001). Behaviorual relevance of atypical language
%       lateralization in healthy subjects. Brain, 124, 1657-1665.
%
% "To control for goodness of the Doppler sonographic recording, the root
% mean square of the averaged CBFV fluctuations in the control condition
% (-17 to -3 s before the cueing tone) was calculated.  All recordings
% with root mean square values greater than 2% of the averaged CBFV were
% discarded."
%
% Base & Stim Only event types: (see dopTaskDef) 10-11-10
% can't calculate the LI-variabilty for base or stim only (currently)
% can't calculate the Goodness of Recording for stim only
%
fprintf('\n\t%s\n\t\t%s',['Running ',mfilename],...
    'Calculating Event Data Indices/Descriptives:');
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
    dic=setGetInputsStruct(inputs);
    
    dic.translation='data index calculation';
    dic.structAbbreviations='dic';
    dic.LRdata=LRdata; clear LRdata
    dic.filter=filters; clear filters
    
    dic.data=cat(3,dic.LRdata,diff(flipdim(dic.LRdata,3),1,3),mean(dic.LRdata,3)); % L:R:difference data

    % epoch x sample x 4 array
    dic.nEpochs=size(dic.data,1);
    dic.out.epoch.n=dic.nEpochs;
    %% epoch screening?
    dic.scrn=ones(dic.nEpochs,1)==1; % use all by default
    if ~isempty(dic.screen)
        if length(dic.screen)==dic.nEpochs
            dic.scrn=logical(dic.screen);
        else
            fprintf('\n%s\n\t%s\n\t%s\n',...
                'Mismatch in Screen filter length:',...
                [num2str(dic.nEpochs),' Epochs, screen size = ',...
                num2str(size(dic.screen))],...
                ' > using all epochs (no screening).');
        end
    end
    dic.data=dic.data(dic.scrn,:,:); % update data, removing screened epochs
    dic.nEpochs=size(dic.data,1);
    dic.out.epoch.n=dic.nEpochs;
    %% for each filter
    dic.fns=fieldnames(dic.filter); % filter names
    for i=1:length(dic.fns)
        dic.fn=dic.fns{i};
        switch dic.fn
            case 'zeroPt' % don't do anything for this
                
            case  {'left','right','avg','dif'} % index filters   
                dic.tmp.fnsj=fieldnames(dic.filter.(dic.fn));
                for j=1:numel(dic.tmp.fnsj);dic.tmp.fnj=dic.tmp.fnsj{j};
                    dic.tmp.fnsjj=fieldnames(dic.filter.(dic.fn).(dic.tmp.fnj));
                    for jj=1:numel(dic.tmp.fnsjj);dic.tmp.fnjj=dic.tmp.fnsjj{jj};
                            dic.tmp.vn=[dic.fn,dic.tmp.fnj]; % variable name
                            [dic.out.(dic.tmp.vn) dic.out.all.(dic.tmp.vn)]=...
                                indCalc(dic.data(:,dic.filter.(dic.fn).(dic.tmp.fnj).(dic.tmp.fnjj).filter,:));
                           dic.out.(dic.tmp.vn).time=dic.filter.(dic.fn).(dic.tmp.fnj).(dic.tmp.fnjj).sec;
                           dic.out.all.(dic.tmp.vn).time=dic.out.(dic.tmp.vn).time; 
                    end
                end
            otherwise % calculate indices
                %% > mean, std, min, max
                [dic.out.(dic.fn) dic.out.all.(dic.fn)]=indCalc(dic.data(:,dic.filter.(dic.fn),:));
                %                 dic.out.(dic.fn).mean=mean(dic.data(:,dic.filter.(dic.fn),:),2);
                %                 dic.out.(dic.fn).std=std(dic.data(:,dic.filter.(dic.fn),:),1,2);
                %                 dic.out.(dic.fn).min=min(dic.data(:,dic.filter.(dic.fn),:),[],2);
                %                 dic.out.(dic.fn).max=max(dic.data(:,dic.filter.(dic.fn),:),[],2);
                %                 dic.out.(dic.fn).rms=sqrt(mean(dic.data(:,dic.filter.(dic.fn),:).^2,2)...
                %                     /size(dic.data(:,dic.filter.(dic.fn),:),2));
                
        end
    end
    %% > Stim Only Event type:
    % for stim only event type, input baseline data
    if dic.evType==1 && ~isempty(dic.baseLRdata)
        dic.baseData=cat(3,dic.baseLRdata,dic.baseLRdata(:,:,1)-dic.baseLRdata(:,:,2));%diff(dic.baseLRdata,1,3));
        dic.baseData=dic.baseData(logical(dic.scrn),:,:); % update data, removing screened epochs
        dic.fn='base';
        [dic.out.(dic.fn) dic.out.all.(dic.fn)]=indCalc(dic.baseData(:,:,:));
        %         dic.out.(dic.fn).mean=mean(dic.baseData(:,:,:),2);
        %         dic.out.(dic.fn).std=std(dic.baseData(:,:,:),1,2);
        %         dic.out.(dic.fn).min=min(dic.baseData(:,:,:),[],2);
        %         dic.out.(dic.fn).max=max(dic.baseData(:,:,:),[],2);
        %         dic.out.(dic.fn).rms=sqrt(mean(dic.baseData(:,:,:).^2,2)...
        %             /size(dic.baseData(:,:,:),2));
    end
    %% Laterality Index Variability
    % not sure about this one yet - they use a different procedure in the
    % paper - single trial tasks...
    dic.out.livTrans='liv = Laterality Index Variability (see Vadikolias et al., 2007,J Cereb Blood Flow Metab)';
    if isfield(dic.out,'epoch') && isfield(dic.out,'base') && ~isempty(dic.out.base)
        dic.out.liv=cat(3,(dic.out.epoch.max(:,1,1)-dic.out.epoch.min(:,1,1)... % left max-min
            ./dic.out.base.mean(:,1,1)),... % left mean base
            (dic.out.epoch.max(:,1,2)-dic.out.epoch.min(:,1,2)... % right max-min
            ./dic.out.base.mean(:,1,2))); % right mean base
        dic.out.liv(:,1,3)=dic.out.liv(:,1,1)-dic.out.liv(:,1,2);
%         dic.out.all.liv=mean(dic.out.liv);
    end
    %% Goodness of Recording
    % "To control for goodness of the Doppler sonographic recording, the root
    % mean square of the averaged CBFV fluctuations in the control condition
    % (-17 to -3 s before the cueing tone) was calculated.  All recordings
    % with root mean square values greater than 2% of the averaged CBFV were
    % discarded."
    %
    % via e-mails and deductions...
    % gor = rms(baseline)
    % if gor > 2 then this is >2% of average CBFV (set to 100 during
    % normalisation).
    %
%     dic.out.gor.trans=fprintf('gor = goodness of recording (see Knecht et al., 2001, Brain)\n',...
%         'Left and Right channels\n');
%     %     dic.out.gorStd=sqrt(mean(dic.out.base.std.^2,2)/dic.nEpochs);
%     if isfield(dic.out,'base')
%         %         dic.out.gor=mean(sqrt(mean(dic.out.base.std.^2,2)/dic.nEpochs));
%         dic.out.gor.mean=mean(dic.out.base.rms(:,:,1:2));
%         dic.out.gor.max=max(dic.out.base.rms(:,:,1:2));
%     end
    %     dic.out.gor=dic.out.gorMax(1:2)<2; % keep gor < 2%
    
    %% index description
    % include a description of the indices
    dic.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}=dic.out;
    fprintf('\t%s\n','Done');
catch err
    %% catch dopOSCCI error
    % 
save(dopOSCCIdebug(mfilename));rethrow(err); 
%     dopOSCCIcrash(err,dic,mfilename);
end

clear dic inputs
end
%%-------------------
%% embedded function
%% > indices
% date in, indices calculations out:
function [ep all]=indCalc(data)
ep=[]; % epoch - number for each epoch (left, right, diff)
ep.n=size(data,1);
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));
ep.med=median(data(:,:,:),2);
ep.iqr=iqr(data(:,:,:),2);

% gor = goodness of recording (see Knecht et al., 2001, Brain)
% avgData = mean(data);
% sqrt(mean(mean(data).^2);
% these don'r really make any sense for each trial/epoch
ep.gorMean = ep.rms(:,:,:);
ep.gorMax = ep.rms(:,:,:); 


all=[]; % all epochs (left, right, diff, average)
all.n=ep.n;
if all.n==1;
    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.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);
     %% Goodness of Recording
    % "To control for goodness of the Doppler sonographic recording, the root
    % mean square of the averaged CBFV fluctuations in the control condition
    % (-17 to -3 s before the cueing tone) was calculated.  All recordings
    % with root mean square values greater than 2% of the averaged CBFV were
    % discarded."
    %
    % via e-mails and deductions...
    % gor = rms(baseline)
    % if gor > 2 then this is >2% of average CBFV (set to 100 during
    % normalisation).
    %
    all.gorMean = sqrt(mean(mean(data).^2));%mean(ep.rms); % same as rms...
    all.gorMax = sqrt(mean(max(data).^2));
end

clear data;
end
