function challenge(recordName)
%
% challenge(recordName)
%
% Octave-compatible code for the 2014 PhysioNet/CinC Challenge
%
% This function takes only one input argument:
%
% recordName
%       String specifying the name of the record to process.  Do not include
%       the '.dat' or '.hea' suffix in recordName.
%
% This function has no output arguments, but it writes an annotation file
% named "recordName.qrs" in the current directory, containing a annotation
% for each detected beat in the record.
%
% Dependencies:
%
%       1) This function requires the WFDB Toolbox for MATLAB and Octave.
%          For information on how to install the toolbox, please see:
%              http://www.physionet.org/physiotools/matlab/wfdb-app-matlab/
%
%       2) The Toolbox is supported only on 64-bit MATLAB 2013a and 2013b
%          and on GNU Octave 3.6.4 and later, on Linux, Mac OS X, and Windows.
%
% This file was written as part of a sample entry for the 2014 Challenge;  for
% the complete sample entry including all other required components, see
%    http://physionet.org/challenge/2014/octave/
% For a description of the Challenge and how to participate in it, see
%    http://physionet.org/challenge/2014/
%
% Version 0.91
%
% Written by Ikaro Silva, December 10, 2013.
% Last Modified: February 14, 2014

% Using the WFDB Toolbox's wfdbloadlib function, initialize the Toolbox
% configuration to the default values.  (This is required for Octave.)
[~,config] = wfdbloadlib;
% For more information, run 'help wfdbloadlib'.

annName = 'qrs'; % This annotator name is required by the Challenge rules.

% Although the code in this example does not directly access the raw samples
% of the record, they can be loaded into an array for processing using the
% WFDB Toolbox's rdsamp function, like this:
%   [tm,signal] = rdsamp(recordName);
% For more information, run 'help rdsamp'.

% The Toolbox's wfdbdesc function reads metadata for the record, including
% the names of the signals.
siginfo = wfdbdesc(recordName);
descriptions = squeeze(struct2cell(siginfo));
descriptions = descriptions(5,:);
% For more information, run 'help wfdbdesc'.

% Using the helper function defined below, find the ECG and blood pressure
% signals available in the record.
[ecg_ind, bp_ind] = get_ECG_BP_index(descriptions);


% Using the Toolbox's gqrs function, annotate the first ECG signal to generate
% an annotation file (recordName.qrs), using default values for gqrs's N, N0,
% and threshold parameters.
N = [];
N0 = [];
threshold = [];

Delay_min = 0.1;
Delay_max = 0.4;

if ~isempty(ecg_ind)
    gqrs(recordName,N,N0,ecg_ind(1),threshold,annName);
    [ANN_GQRS,type,subtype,chan,num] = rdann(recordName,'qrs');
    
    [TM, signal, Fs] = rdsamp(recordName,[]);
    SamPer=1/Fs;
    
    ECG=signal(:,ecg_ind(1));
    ANN_GQRS=adjustGQRSByPeak(ANN_GQRS, ECG, Fs);

    % If there is a blood pressure signal, there is more processing to do.
    if (~isempty(bp_ind)) 
%         ECG=signal(:,ecg_ind(1));
        BP=signal(:,bp_ind);
        
        All_ANN_WABP=cell(1,length(bp_ind));
        AllBeatQ=cell(1,length(bp_ind));
        NumGoodBPPeaks = zeros(1,length(bp_ind));
        for i=1:length(bp_ind)
            [All_ANN_WABP{i}, AllBeatQ{i}] = getWABPjSQI(BP(:,i), Fs);
            NumGoodBPPeaks(i) = length(find(AllBeatQ{i}==0));
        end
        [NumMax, BestBP_ind] = max(NumGoodBPPeaks);
        ANN_WABP = All_ANN_WABP{BestBP_ind};
        BeatQ = AllBeatQ{BestBP_ind};
        Good_ANN_WABP = ANN_WABP(find(BeatQ==0));
        
        % Convert samples to time
%         ANN_GQRS_TM=ANN_GQRS*SamPer;
%         warning('Marker 1')
%         warning(['Good_ANN_WABP length: ', num2str(length(Good_ANN_WABP))])
%         warning(['ANN_GQRS length: ', num2str(length(ANN_GQRS))])
        if length(Good_ANN_WABP)>length(ANN_GQRS)/5
            
%             warning('Marker 2')
            [Val_QRS, CorresWABP] = sandwich(ANN_GQRS, Good_ANN_WABP);

            %QRS validated by corresponding ABP onset
%             Val_QRS = ANN_GQRS(SandwichedQRSIdx);
%             CorresWABP = Good_ANN_WABP(CorresWABPIdx);
            
            %ABP onset which does not correspond to a validated QRS
%             NotCorresWABPIdx = 1:length(Good_ANN_WABP);
%             NotCorresWABPIdx(CorresWABPIdx)=[];
%             NotCorresWABP=Good_ANN_WABP(NotCorresWABPIdx);

%             warning('Marker 3')
            NotCorresWABP = setdiff(Good_ANN_WABP, CorresWABP);
            
%             warning('Marker 4')
            if length(Val_QRS)==length(CorresWABP)
                Delay = CorresWABP - Val_QRS;
            else
                Delay = CorresWABP - Val_QRS(1:length(CorresWABP));
            end
            
            if isempty(Delay) || median(Delay)*SamPer < Delay_min || median(Delay)*SamPer > Delay_max
                    delta_t = round(Fs/5);
                    ANN_QRS = Good_ANN_WABP - delta_t;
            else
                [Delay_refined, Delay_Normal_Idx, Delay_Large_Idx, Delay_Small_Idx, Delay_Left_max, Delay_Right_max] = hampel(Delay);
                Val_QRS(Delay_Large_Idx) = CorresWABP(Delay_Large_Idx) - Delay_refined(Delay_Large_Idx);
                Val_QRS(Delay_Small_Idx) = CorresWABP(Delay_Small_Idx) - Delay_refined(Delay_Small_Idx);

            
%             warning('Marker 5')
                if ~isempty(Delay_refined)
                    ABP_Derived_QRS = NotCorresWABP - round(median(Delay_refined));
                else
                    ABP_Derived_QRS = NotCorresWABP - round(Fs/5);
                end
                ABP_Derived_QRS = ABP_Derived_QRS(ABP_Derived_QRS>0);
            
%             warning('Marker 6')
                NotVal_QRS = setdiff(ANN_GQRS, Val_QRS);
                NotValGood_QRS=[];
                for i=1:length(NotVal_QRS)
                    if min(abs(Good_ANN_WABP-NotVal_QRS(i))) > round(2*Fs)
                        NotValGood_QRS = [NotValGood_QRS; NotVal_QRS(i)];
                    end
                end
            
%             warning('Marker 7')
                ANN_QRS = [Val_QRS; ABP_Derived_QRS; NotValGood_QRS];
                ANN_QRS = sort(ANN_QRS);
            end

            %Add Hampel filter
            leftAddQRS=[];
            rightAddQRS=[];
            middleAddQRS=[];

            RR = ANN_QRS(2:end)-ANN_QRS(1:end-1);
            [RR_refined, RR_Normal_Idx, RR_Large_Idx, RR_Small_Idx, RR_Left_max, RR_Right_max] = hampel(RR);
            if ANN_QRS(1)>RR_Left_max
                tmp_RR=RR_refined(1);
                N_RR=floor(ANN_QRS(1)/tmp_RR); %how many QRS peaks to add to the left
                leftAddQRS = ANN_QRS(1)-(N_RR:-1:1)'*tmp_RR;
            end

            if length(TM)-ANN_QRS(end)>RR_Right_max
                tmp_RR=RR_refined(end);
                N_RR=floor((length(TM)-ANN_QRS(end))/tmp_RR); %how many QRS peaks to add to the right
                rightAddQRS = ANN_QRS(end)+(1:N_RR)'*tmp_RR;
            end

            middleAddQRS=[];
            for i=1:length(RR_Large_Idx)
                tmp_RR = RR_refined(RR_Large_Idx(i));
                N_RR=round(RR(RR_Large_Idx(i))/tmp_RR);
                if N_RR>1
                    middleAddQRS=[middleAddQRS; ANN_QRS(RR_Large_Idx(i))+(1:N_RR-1)'*tmp_RR];
                end
            end

            ANN_QRS=sort([leftAddQRS; ANN_QRS; middleAddQRS; rightAddQRS]);
            ANN_QRS=ANN_QRS(ANN_QRS>0 & ANN_QRS<length(TM));
%             warning('Marker 8')
            if ~isempty(ANN_QRS)
                wrann(recordName,annName,ANN_QRS,[],[],[],[]);
            end
%         else
%             warning('Marker 8')
%             wrann(recordName,annName,ANN_GQRS,[],[],[],[]);
        end
        
%     else
%         wrann(recordName,annName,ANN_GQRS,[],[],[],[]);
    end

elseif (~isempty(bp_ind)) %There is no recrod without ECG but with BP.
    wabp(recordName, [], [], [], bp_ind(1));
    [ANN_WABP,type,subtype,chan,num] = rdann(recordName,'wabp');
    [TM, signal, Fs] = rdsamp(recordName,[]);
    delta_t = round(Fs/5);
    ANN_QRS = ANN_WABP - delta_t;
    wrann(recordName,annName,ANN_QRS,[],[],[],[]);
else
    ANN_QRS=[];
    for i=1:size(descriptions,1)
        gqrs(recordName,N,N0,i,threshold,annName);
        [ANN_GQRS,type,subtype,chan,num] = rdann(recordName,'qrs');
        if length(ANN_GQRS)>length(ANN_QRS)
            ANN_QRS=ANN_GQRS;
        end
    end
    wrann(recordName,annName,ANN_QRS,[],[],[],[]);
end
    

