function varargout = process_weighted_merging( varargin )
% PROCESS_WEIGHTED_MERGING: Merges left and right ABR channels.
% Merging process is using weighted averaging process.
% Method taken from "Acoustic change complex for assessing speech
% discrimination in normal-hearing and hearing-impaired infants"
% By Ching et al. (2023) Clin Neurophysiol, 149, 121-132.
% Author: by MinChul Park (June 2022)
% University of Canterbury | Te Whare Wānanga o Waitaha
% Christchurch | Ōtautahi
% New Zealand | Aotearoa
% Contributors: François Tadel and Raymundo Cassani

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Weighted merging';
    sProcess.FileTag     = 'WMerging';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Custom Processes';
    sProcess.Index       = 1000;
    sProcess.Description = 'https://www.sciencedirect.com/science/article/pii/S1388245723002195?via%3Dihub';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data'};
    sProcess.OutputTypes = {'data'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;

    % Description
    sProcess.options.info.Comment = ['Merges left and right ABR channels through weighted averaging process.<BR>' ... 
                                     'Method from "Acoustic change complex for assessing speech discrimination<BR>'...
                                     'in normal-hearing and hearing-impaired infants".<BR>' ...
                                     'Ching et al. (2023) Clin Neurophysiol, 149, 121-132."<BR><BR>'...
                                     'Merges the left and right ABR channels and produces one ABR waveform.<BR>'];
    sProcess.options.info.Type    = 'label';
    sProcess.options.info.Value   = [];
end


%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
     Comment = sProcess.Comment;
end


%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputs) %#ok<DEFNU>
    % Initialize returned list of files
    OutputFiles = {};

    % ===== LOAD THE DATA =====
    % Read the first file in the list, to initialize the loop
    DataMat = in_bst(sInputs(1).FileName, [], 0);
    epochSize = size(DataMat.F);
    Time = DataMat.Time;
    % Initialize the load matrix: [Nchannels x Ntime x Nepochs]
    AllMat = zeros(epochSize(1), epochSize(2), length(sInputs));
    % Reading all the input files in a big matrix
    for i = 1:length(sInputs)
        % Read the file #i
        DataMat = in_bst(sInputs(i).FileName, [], 0);
        % Check the dimensions of the recordings matrix in this file
        if ~isequal(size(DataMat.F), epochSize)
            % Add an error message to the report
            bst_report('Error', sProcess, sInputs, 'One file has a different number of channels or a different number of time samples.');
            % Stop the process
            return;
        end
        % Add the current file in the big load matrix
        AllMat(:,:,i) = DataMat.F;
    end
    
    % ===== PROCESS =====
    % This is where the actual process of data manipulation and calculation takes place.
    Nchannels = epochSize(1); % Number of channels. 
    Nepochs   = length(sInputs); % Number of epochs.  
    Ntime     = epochSize(2); % Number of time samples. In the ABR it is 493
    EpochW    = ones(Nchannels,1,Nepochs)./var((AllMat),0,2); % Calculation of epoch specific weighting.
    WEpoch    = EpochW.*(AllMat); % Application of epoch specific weighting to all epochs. 
    WAvg      = sum(WEpoch,3)./sum(EpochW,3); % Calculation of weighted averaging through the weighted epochs. 
    WAvgL     = WAvg(4,:); % Selects the left channel; this will be swapped if the channel rows are swapped
    WAvgR     = WAvg(5,:); % Selects the right channel; this will be swapped if the channel rows are swapped
 
    AvgWAvg   = mean(WAvg,2); % Mean of the weighted average ABR. 
    Squared   = (WAvg - AvgWAvg).^2; % (EEG value at time sample N - AvgWAvg)^2
    VarABR    = (1/(Ntime+1))*sum(Squared,2); % (1/Number of time samples + 1)*the sum of Squared
    VarL      = VarABR(4,:); % Selects the left channel
    VarR      = VarABR(5,:); % Selects the right channel
    WFactorL  = VarR/(VarR+VarL); % Weighting factor - left
    WFactorR  = VarL/(VarR+VarL); % Weighting factor - right
    WMerging = (WFactorL*WAvgL) + (WFactorR*WAvgR); % Final weighted merging calculation

    WAvg(4,:) = WMerging;
    WAvg(5,:) = WMerging;

    % ===== SAVE THE RESULTS =====
    % Get the output study (pick the one from the first file)
    iStudy = sInputs(1).iStudy;
    % Create a new data file structure
    DataMat             = db_template('datamat');
    DataMat.F           = WAvg;
    DataMat.Comment     = sprintf('WMerging (%d)', length(sInputs)); % Names the output file as 'WMerging' with the number of epochs used to generate the file.
    DataMat.ChannelFlag = ones(epochSize(1), 1);   % List of good/bad channels (1=good, -1=bad)
    DataMat.Time        = Time;
    DataMat.DataType    = 'recordings';
    DataMat.nAvg        = length(sInputs);         % Number of epochs that were averaged to get this file
    % Create a default output filename 
    OutputFiles{1} = bst_process('GetNewFilename', fileparts(sInputs(1).FileName), 'data_WMerging');
    % Save on disk
    save(OutputFiles{1}, '-struct', 'DataMat');
    % Register in database
    db_add_data(iStudy, OutputFiles{1}, DataMat);
end
