function  challengeExampleContinuationIdea(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.)
tic;
[~,config] = wfdbloadlib;
verbose = false;
% 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);
description = squeeze(struct2cell(siginfo));
description = description(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 = get_index(description,'ECG');
bp_ind = get_index(description,'BP');

% 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 = [];
cleanThre = 0.1;
gqrs(recordName,N,N0,ecg_ind(1),threshold,'pecg');
% For more information, run 'help gqrs'.
iwantblood = true; %% dummy variable to pretend whe dont have blood information
% If there is a blood pressure signal, there is more processing to do.
[rr_ecg,~] = ann2rr(recordName,'pecg');
mlMu = median(rr_ecg);
mlStd = std(rr_ecg);
selEcg = noOutlier(rr_ecg,mlMu,mlStd,cleanThre);
cleanrr_ecg = rr_ecg(selEcg);
heartMu = median(cleanrr_ecg);
ecg_sample = rdann(recordName,'pecg');
if (~isempty(bp_ind)&& iwantblood)  
	fprintf('blood there\n');
	% Using the Toolbox's wabp function, annotate the first ABP, ART or BP signal
	% with wabp to generate a second annotation file (recordName.wabp). (wabp
	% finds this signal itself -- we don't need to tell it which signal to use.)
	wabp(recordName);
	% For more information, run 'help wabp'.

	% Using the Toolbox's ann2rr function, read the BP annotations from wabp, and
	% fill an array with the inter-beat (RR) intervals.
	[rr_abp,~] = ann2rr(recordName,'wabp');
	% For more information, run 'help ann2rr'.

	% Fill another array with the RR intervals from the ECG annotations.

	% Compare the variances of the BP- and ECG-based RR intervals, and assume
	% that the RR series with the lower variance is best.
	[~,best] = min([var(rr_abp) var(rr_ecg)]);
	
	% If the BP-based RR intervals are best, overwrite recordName.qrs
	if(best == 1)
		% This example assumes, naively, that the delay from the QRS annotations
		% to the BP annotations is always 200 ms.  We need to determine the
		% number of sample intervals corresponding to 200 ms (which is the
		% sampling frequency divided by 5).
		


		fprintf('blood is best\n');
		trust_ecg_sample = ecg_sample(selEcg);

		%%now bloodpresure
		[rr_abp,~] = ann2rr(recordName,'wabp');
		blood_sample = rdann(recordName,'wabp');
		if verbose
		fprintf('pre wrann\n');
		end
		mwrann(recordName,annName,blood_sample);
		if verbose
		fprintf('post wrann\n');
		end
		trust_blood_sample = blood_sample;
		delta = nan*zeros(numel(trust_blood_sample),1);
		counter = 1;
		pointerEcg = numel(trust_ecg_sample);
		for i=numel(trust_blood_sample):-1:1
			%%get ecg sample just before it
			while trust_ecg_sample(pointerEcg) >= trust_blood_sample(i) && pointerEcg >=2
				pointerEcg = pointerEcg-1;
			end
			if trust_ecg_sample(pointerEcg) < trust_blood_sample(i) %%recheck for pointerEcg ==1
				delta(i) = trust_blood_sample(i) - trust_ecg_sample(pointerEcg);
			end
		end
		%%outlier stuff
		%%if distance bigger than heartfrequency raus damit
		delta = delta(delta<heartMu);
		muDelta = mean(delta);

		% Using the Toolbox's rdann function, read the BP annotations into a
		% set of vectors.
		trust_blood_sample = trust_blood_sample-round(muDelta);
		% For more information, run 'help rdann'.

		% Adjust the annotation times (sample numbers) by -200 ms.

		% Using the Toolbox's wrann function, write the time-shifted BP annotations
		% into recordName.qrs (replacing the original ECG annotations).
		if verbose
		fprintf('pre wrann\n');
		end
		mwrann(recordName,annName,trust_blood_sample);
		if verbose
		fprintf('post wrann\n');
		end
		% For more information, run 'help wrann'.
	end% If the ECG-based RR intervals are best, recordName.qrs does not need to
	% be changed.
end % blood not available
if isempty(bp_ind) || not(iwantblood) || best==2
	%%get the most out of the ecg signal as possible
	curHeart = heartMu;
	fprintf('ecg is used\n');
	%%initialize our outlieer detectors
	
	%based on the heartBeat
	
	%based on noise in the signal
	[tm,signal] = rdsamp(recordName);
	%get first derivative
	derivative = nan*zeros(numel(tm)-1,1);
	for i=1:(numel(tm)-1)
		derivative(i) = signal(i+1)-signal(i);
	end
	windowLength = heartMu;
	windowHalf = round(heartMu/2);
	stdWindow = nan*zeros(numel(ecg_sample),1); 
	for i=1:numel(ecg_sample)
		lb = max((ecg_sample(i)-windowHalf),1);
		ub = min((ecg_sample(i)+windowHalf),numel(derivative));
		stdWindow(i) = std(derivative(lb:ub));
	end
	%%based on the last heartbeat and the heartMu search for the closest ecg annotation
	%%always accept first sample
	sel = false(numel(ecg_sample),1);
	sel(1) = true;
	wantSamples = nan*zeros(numel(ecg_sample)+100,1);
	wantSamplesCounter = 1;
	addSamples = nan*zeros(numel(ecg_sample),1);
	addSamplesCounter = 0;
	curSampleAdd = false;
	i = 1;
	while i<=numel(ecg_sample)-1
		%%we already believe that this sample is valid
		if curSampleAdd
			curSample = addSamples(addSamplesCounter); 
			proceed = true;
		else
			proceed = sel(i);
			curSample = ecg_sample(i);
		end
		if proceed
			%%based on the variance at the next sample do we trust the next sample?
			if nonParametricOneSidedlikelihood(stdWindow(i+1),stdWindow) > 0.3
				%yes we trust the sample
				j=i+1;
			else
				%%based on the current heart rate where should the next sample lie
				wantSample = curSample+curHeart;
				wantSamples(wantSamplesCounter) = wantSample;
				wantSamplesCounter = wantSamplesCounter +1;
				%%get the closest sample to nextSample
				j=i+1;
				curDis = abs(ecg_sample(j)-wantSample);
				while j<numel(ecg_sample) && abs(ecg_sample(j+1)-wantSample)<curDis
					%fprintf('I was here for %d marker',i);
					j=j+1;
					curDis = abs(ecg_sample(j)-wantSample);
				end
			end 
			%%check if the selected sample is plausible at all, if not replace 
			%with our best guess
			if ecg_sample(j) > curSample+2*curHeart
				addSamplesCounter = addSamplesCounter+1;
				addSamples(addSamplesCounter) =  curSample+curHeart;
				curSampleAdd = true;
			else
				curSampleAdd = false;
				sel(j) = true;
				curHeart = ecg_sample(j)-ecg_sample(i);
				i = i +1;
			end
			%if ecg_sample(i) > 23300
			%	notchosen = ecg_sample(setdiff((i+1):(j+1),j));
			%	wrann(recordName,'dbg',[ecg_sample([i,j]);notchosen;wantSample],[repmat('N',2,1);repmat('B',numel(notchosen),1);repmat('L',1,1)]);
			%	error('debug here');
			%end
			
			
		else
			i = i+1;
		end
	end
	cleanEcg = ecg_sample(sel);
	wantSamples = wantSamples(not(isnan(wantSamples)));
	if verbose
		fprintf('pre wrann\n');
	end
	toc
	mwrann(recordName,annName,cleanEcg);
	if verbose
		fprintf('post wrann\n');
	end
	%wrann(recordName,'hrt',1:curHeart:5000);
	%wrann(recordName,'wnt',wantSamples(not(isnan(wantSamples))));
end
 
end % of function
%%%%%%%%%%%% Helper Function get_index %%%%%%%%%%%%%%%%%%%%%
% Inputs:
%   description  strings containing signal descriptions from wfdbdesc
%   pattern      a string to be found in the descriptions
% Output:
%   ind          indices of the descriptions that contain the pattern

function ind = get_index(description,pattern)
	M = length(description);
	tmp_ind = strfind(description,pattern);
	ind = [];
	for m = 1:M
		if(~isempty(tmp_ind{m}))
			ind(end+1) = m;
		end
	end
end

function sel = noOutlier(x,mu,std1,threShold)
	keepFun = @(x)normcdf(x,mu,std1) > threShold & 1-normcdf(x,mu,std1) > threShold; %just the likelihood under the model
	sel = keepFun(x);
end

function sel = isOutlier(x,mu,std1,threShold)
	outl = @(x)1-normcdf(x,mu,std1) < threShold; %just the likelihood under the model
	sel = outl(x);
	sel(isnan(sel)) = false;
end

%function mwrann(recordName,annName,samples)
%	N = numel(samples);
%	wrann(recordName,annName,samples,repmat('N',N,1),repmat(0,N,1),repmat(0,N,1),repmat(0,N,1))
%end

function lh = oneSidelikelihood(x,mu,std1)
	lh = 1-normcdf(x,mu,std1);
end

function lh = nonParametricOneSidedlikelihood(x,data)
	lh = 1-sum(x>data)/numel(data);
end

function mwrann(recordName,annName,cleanEcg);
	wrann(recordName,annName,cleanEcg);
end
