function challenge(recordName)
%recordName='data/153';


%% Init

fprintf('\tInit...');pause(0);drawnow;
% Using the WFDB Toolbox's wfdbloadlib function, initialize the Toolbox
% configuration to the default values.  (This is required for Octave.)
[~,config] = wfdbloadlib;

annName = 'qrs'; % This annotator name is required by the Challenge rules.
fprintf('done\n');pause(0);drawnow;


%% Load Signals
% 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:
fprintf('\tLoading Signals...');pause(0);drawnow;
[tm,signal,samplerate] = rdsamp(recordName);

% 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,:);


%samplerate = str2num(strtok(siginfo(1).SamplingFrequency));
%assert(~isnull(samplerate));

ecg_idxs = get_index(description,'ECG');
if isempty(ecg_idxs)
	fprintf('No ECG Found\n');
	ecg_idxs = 1;
end

fprintf('done\n');pause(0);drawnow;


%For now choose only the first ecg...
ecg_idxs = ecg_idxs(1);



for ii=1:length(ecg_idxs)
	ecg_idx = ecg_idxs(ii);
	fprintf('\tWorking on ECG in Signal %i (%i/%i)...\n', ecg_idx, ii, length(ecg_idxs));pause(0);drawnow;

	raw_ecg = signal(:,ecg_idx);


%% ECG Handling (Noise etc...
%tic
	fprintf('\tRemoving NaNs from the ECG...');pause(0);drawnow;
	nonan_ecg = ecg_interpolate_nan(raw_ecg);
	fprintf('done\n');pause(0);drawnow;
%toc

%tic
	fprintf('\tRemoving Powerline Interference from the ECG...');pause(0);drawnow;
	hd1 = create_filter_cheby(8,[57.5 62.5] ./ (samplerate / 2), 'stop');
	if (samplerate > 240)
		hd2 = create_filter_cheby(8,[118 122] ./ (samplerate / 2), 'stop');
		nopower_ecg = filtfilthd(hd2, filtfilthd(hd1, nonan_ecg));
	else
		nopower_ecg = filtfilthd(hd1, nonan_ecg);
	end
	%clear('hd1', 'hd2');
	fprintf('done\n');pause(0);drawnow;
%toc

%tic
	fprintf('\tCalculate ECG Baseline...');pause(0);drawnow;
	ecg_baseline = ecg_calculate_baseline(nopower_ecg, samplerate);
	nobaseline_ecg = nopower_ecg - ecg_baseline;
	fprintf('done\n');pause(0);drawnow;
%toc

%tic
	fprintf('\tDetect R Peaks...');pause(0);drawnow;
	list_r_peaks = find_r_peaks(nobaseline_ecg, samplerate);
	fprintf('done\n');pause(0);drawnow;
%toc


	if true 
	    fprintf('\tEstimate Noise...');pause(0);drawnow;
%tic
	    noise_estimate = ecg_estimate_noise(nobaseline_ecg, list_r_peaks, samplerate) + ecg_baseline;
%toc
	    stride_length = samplerate;
	    window_radius = 5;
%tic
	    [ecg_signal_power, ecg_noise_power, ecg_snr, ecg_noise_color, ecg_corrected_snr] = calculate_signal_properties (nobaseline_ecg, noise_estimate, samplerate, stride_length, window_radius);
	    %peak_dynamics_vector = calculate_peak_dynamics(list_r_peaks, s_nobaseline_ecg(list_r_peaks), samplerate, stride_length, window_radius, length(s_raw_ecg));    
%toc
	    fprintf('done\n');pause(0);drawnow;
	else
	    ecg_corrected_snr = [];
	    stride_length = 0;
	end

%tic
	fprintf('\tSelect Probably Correct Beats...');pause(0);drawnow;
	[list_clean_r_peaks, no_r_peak_areas_start, no_r_peak_areas_end, removed_beats] = clean_r_peaks(list_r_peaks, nobaseline_ecg, ecg_corrected_snr, stride_length, samplerate);
	fprintf('done\n');pause(0);drawnow;
%toc
end


if (length(no_r_peak_areas_start) > 0)
    fprintf('\tFind Cardiac Artifacts in other Signals...');pause(0);drawnow;
    other_sig_idxs = setdiff(1:length(description), [ecg_idx]);

    results = struct('signal_id', {}, 'peak_idxs', {}, 'cardiac_ratio', {});
    cur_result_id = 0;

    for ii=other_sig_idxs
%tic
	try
	    raw_signal = signal(:,ii);
	    [cardiac_shape, peak_offset, cardiac_ratio] = find_cardiac_artifact_shape(raw_signal, list_clean_r_peaks, samplerate);
	    if cardiac_ratio < 0.3
%toc
		continue
	    end
	    signal_xcov = xcov(raw_signal, cardiac_shape);
	    signal_xcov = signal_xcov((ceil(length(signal_xcov) / 2)-peak_offset):(end-peak_offset));
	    mean_peak_height = mean(signal_xcov(list_r_peaks));
	    [signal_peaks, signal_peak_vals] = find_peaks(signal_xcov, mean_peak_height / 3);
	    cur_result_id = cur_result_id +1;
	    results(cur_result_id).signal_id = description{ii};
	    results(cur_result_id).peak_idxs = signal_peaks;
	    results(cur_result_id).cardiac_ratio = cardiac_ratio;
	catch
		fprintf('CATCH BLOCK SOMETHING IS FISHY!!!!!\n');
	    	cur_result_id = cur_result_id +1;
		results(cur_result_id).signal_id = 'Broken';
		results(cur_result_id).peak_idxs = [];
		results(cur_result_id).cardiac_ratio = 0;
	end_try_catch
%toc
    end
    fprintf('done\n');pause(0);drawnow;



    fprintf('\tEstimate Peak Positions...');pause(0);drawnow;
    combined_peaks = removed_beats;
    for ii=1:length(results)
        combined_peaks = [combined_peaks, reshape(results(ii).peak_idxs, 1, [])];
    end
    possible_peaks=zeros(0,1);
    for ii=1:length(no_r_peak_areas_start)
        selected_peaks = (combined_peaks > no_r_peak_areas_start(ii)) &...
                        (combined_peaks < no_r_peak_areas_end(ii));
        possible_peaks = [possible_peaks, reshape(combined_peaks(selected_peaks),1,[])];
    end

    if ~isempty(possible_peaks)
%tic
            density_steps = floor(samplerate/100);
%            reduced_rate = samplerate/density_steps;
            density_domain = 1:density_steps:length(nobaseline_ecg);
            density_filter = false(length(nobaseline_ecg),1);
            for ii=1:length(no_r_peak_areas_start)
                density_filter(no_r_peak_areas_start(ii):no_r_peak_areas_end(ii)) = true;
            end
            assert(max(density_domain) <= length(density_filter));
            density_domain=density_domain(density_filter(density_domain));

            density_estimation = estimate_density_gauss_kernel(possible_peaks, ...
							    density_domain, ...
							    10)  ...
				* length(possible_peaks) * sqrt(2*pi)*10;
%toc
	    [density_peak_idxs,density_peak_heights] = find_peaks(density_estimation, 1.5);
            [estimated_peak_idx, estimated_peak_height] = achieve_min_peak_distance(density_domain(density_peak_idxs),density_peak_heights, round(0.5 * samplerate));
%	    estimated_peak_idx = round(estimated_peak_idx * density_steps);
    else
	    estimated_peak_idx = [];
    end
    fprintf('done\n');pause(0);drawnow;
else
    estimated_peak_idx = [];
end

fprintf('\tCombine Peak Positions...');pause(0);drawnow;
r_peak_positions = unique([reshape(estimated_peak_idx, 1, []), reshape(list_clean_r_peaks, 1, [])]);
fprintf('done\n');pause(0);drawnow;


if false
    close all;
    figure;
    plot(noise_estimate);
    figure;
    plot(nobaseline_ecg)
    hold on
    plot(estimated_peak_idx, nobaseline_ecg(estimated_peak_idx), 'ro');
    %plot(combined_peaks, nobaseline_ecg(combined_peaks), 'rx');
    plot(r_peak_positions, nobaseline_ecg(r_peak_positions), 'gx');
    %plot(removed_beats, nobaseline_ecg(removed_beats), 'go');
    figure;
    plot(ecg_corrected_snr);
    figure;
    ecg_corrected_snr_interp1 = interp1(1:stride_length:max((length(ecg_corrected_snr)-1)*stride_length, length(nobaseline_ecg)), ecg_corrected_snr, 1:length(nobaseline_ecg));
    plot(ecg_corrected_snr_interp1, 'r');
    hold on;
    plot(nobaseline_ecg / max(nobaseline_ecg) * max(ecg_corrected_snr_interp1));
end
%figure
%for ii = 1:length(results)
%  plot(results(ii).peak_idxs, 0 , 'rx');
%end
%plot(list_r_peaks, 0, 'go');


wrann(recordName, annName, r_peak_positions);




