function challenge(recordname)
    [~,config] = wfdbloadlib;
    description = squeeze(struct2cell(wfdbdesc(recordname)));
    ind_ecg = get_index(description(5,:),'ECG');
    ind_bp = get_index(description(5,:),'BP');

    [~, ecg, fs_ecg] = rdsamp(recordname, ind_ecg);
    tol_samples = floor(0.15 * fs_ecg);

    ecg_n = normalize(filter_ecg(ecg));
    
    gqrs(recordname, [], [], ind_ecg);
    [ann_ecg,type,subtype,chan,num] = rdann(recordname, 'qrs');
    ann = find_rhythm(ann_ecg, ecg_n, fs_ecg, 30*fs_ecg);

    if ~isempty(ind_bp)
        %disp('a');
 
        wabp(recordname, [], [], [], ind_bp);
        ann_bp = rdann(recordname, 'wabp');
        [~, bp, fs_bp] = rdsamp(recordname, ind_bp);
        ann_bp(bp(ann_bp)>mean(bp)) = [];
 
        %disp('b');
 
        ann_bp = round(ann_bp*fs_ecg/fs_bp);
 
        max_match = 0;
        best_dist = 0;
        best_delay = 0;
        best_matching = [];
        k = 0;
        for delay = fs_ecg:-25:round(0.1*fs_ecg)
            matching = find_matching(ann, ann_bp - delay, tol_samples);
            if isempty(matching), continue; end
            mean_dist = mean(matching(:,3));
            if size(matching, 1) > max_match || (size(matching, 1) == max_match && mean_dist < best_dist)
                max_match = size(matching, 1);
                best_dist = mean_dist;
                best_delay = delay;
                best_matching = matching;
            end
            k = k + 1;
        end
 
        %disp('c');
 
        if isempty(best_matching), return; end
 
        ann_bp = round(ann_bp - best_delay);
 
        %disp('d');
 
        median_b2b = median(diff(ann_bp));
        has_bp = false(1, length(ecg));
        for j=1:length(ann_bp)
            a = max(1, round(ann_bp(j) - 0.6*median_b2b));
            b = min(length(ecg), round(ann_bp(j) + 0.6*median_b2b));
            has_bp(a:b) = true;
        end
 
        %disp('e');
 
        ann = round(ann_union(ann(~has_bp(ann)), ann_bp, tol_samples));
    end
    
    %disp('f');
 
    ann = round(ann);
    ann(ann<0|ann>length(ecg)) = [];
    ann = reshape(ann, length(ann), 1);
 
    % wrann crashes randomly for some reason on the test machine
    % (works just fine here on my linux of course)
    % saving to ascii then use the wrann shell utility
    %wrann(recordname, 'qrs', ann);
    fd = fopen([recordname '.ascii'], 'w');
    for j=1:length(ann)
        time = ann(j)/fs_ecg;
        m = floor(time/60);
        s = floor(time - m*60);
        ms = round((time - floor(time)) * 1000);
        fprintf(fd, '%5d:%02d.%03d %8d     N    0    0    0\n', m, s, ms, ann(j));
    end
    fclose(fd);
end

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 w = normalize(v, p)
    if nargin < 2
        p = inf;
    end
    w = v - mean(v);
    w = w / norm(w, p);
end

% Filters from Pan-Tompkins algorithm
function ecg = filter_ecg(ecg)
    % low-pass filtering
    % LPF (1-z^-6)^2/(1-z^-1)^2
    b = [1 0 0 0 0 0 -2 0 0 0 0 0 1];
    a = [1 -2 1];
    h_LP = filter(b, a, [1 zeros(1, 12)]);
    ecg = conv(ecg, h_LP, 'same');
    
    % high-pass filtering
    % HPF = Allpass-(Lowpass) = z^-16-[(1-z^-32)/(1-z^-1)]
    b = [-1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 32 -32 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1];
    a = [1 -1];
    h_HP = filter(b, a, [1 zeros(1, 32)]);
    ecg = conv(ecg, h_HP, 'same');
end

function matching = find_matching(a, b, tol)
    assert(isvector(a) && isvector(b));
    ia = 1;
    ib = 1;
    n = 1;
    dist = tol + 1;
    matching = zeros(min(length(a), length(b)), 3);
    while ia <= length(a) && ib <= length(b)
        aa = a(ia);
        bb = b(ib);
        if abs(aa - bb) < dist
            dist = abs(aa - bb);
            matching(n,:) = [ia, ib, dist];
        elseif matching(n, 1) ~= 0
            if matching(n, 1) == ia, ia = ia + 1; end
            if matching(n, 2) == ib, ib = ib + 1; end
            dist = tol + 1;
            n = n + 1;
            continue;
        end
        if aa <= bb, ia = ia + 1; end
        if aa >= bb, ib = ib + 1; end
    end
    if n <= size(matching, 1) && matching(n, 1) ~= 0, n = n + 1; end
    matching(n:end,:) = [];
end

function rhythm = find_rhythm(ann, normalized_ecg, fs, window_size)
    max_skips = 2;
    min_ratio = 0.6;
    max_ratio = 1.66;
    min_rr = 0.2*fs;
    max_rr = 2.5*fs;
    
    ecg_filt = abs(maxfilt(normalized_ecg, floor(0.07*fs)));
    
    n = 1;
    rhythm = zeros(1, length(ann));
    
    rr_window = diff(ann(ann<=window_size));
    
    j = 1;
    while j < length(ann)
        ref_rr = median(rr_window(rr_window>=min_rr&rr_window<=max_rr));
        
        k = j + 1;
        while (k <= length(ann)) && (ann(k) - ann(j) <= (max_skips+1.5)*ref_rr), k = k + 1; end
        k = k - 1;
        
        spans = ann(j+1:k) - ann(j);
        if isempty(spans), j = j + 1; continue; end
        costs = abs(smod(spans, ref_rr));
        skips = round(spans / ref_rr) - 1;
        wrongskips = skips<0|skips>max_skips;
        rrs = zeros(1, length(spans));
        rrs(~wrongskips) = spans(~wrongskips) ./ (skips(~wrongskips) + 1);
        for l=1:length(costs)
            if wrongskips(l), continue; end
            if skips(l)==0, continue; end
            power = min(ecg_filt(ann(j) + round(rrs(l)*(1:skips(l)))));
            wrongskips(l) = power < 0.3;
        end
        costs(wrongskips) = inf;
        if all(isinf(costs)), j = j + 1; continue; end

        [~, ix] = min(costs);

        if rrs(ix) >= min_ratio*ref_rr && rrs(ix) <= max_ratio*ref_rr
            if n == 1
                rhythm(n) = ann(j);
                n = n + 1;
            end
            rhythm(n:n+skips(ix)) = ann(j) + round(rrs(ix)*(1:skips(ix)+1));
            n = n + 1 + skips(ix);
            
            rr_window(1) = [];
            rr_window(end+1) = rrs(ix);
        end
        k = j + ix;
        j = k;
    end
    
    rhythm(n:end) = [];
end

function z = smod(x, y)
    z = mod(x, y);
    big = z>y/2;
    z(big) = z(big)-y;
end

function y = maxfilt(x, n)
    assert(isvector(x));
    m = length(x);
    nh = floor(n/2);
    y = zeros(size(x));
    for j=1:m
        a = max(1, j-nh);
        b = min(m, j+n-nh-1);
        y(j) = max(x(a:b));
    end
end

function ann = ann_union(ann_1, ann_2, tol)
    ann = zeros(1, length(ann_1) + length(ann_2));
    n = 1;
    
    j = 1;
    k = 1;
    while j <= length(ann_1) && k <= length(ann_2)
        a1 = ann_1(j);
        a2 = ann_2(k);
        if a1 < a2 - tol
            ann(n) = a1;
            n = n + 1;
            j = j + 1;
        elseif a2 < a1 - tol
            ann(n) = a2;
            n = n + 1;
            k = k + 1;
        else
            ann(n) = 0.5*(a1 + a2);
            n = n + 1;
            j = j + 1;
            k = k + 1;
        end
    end
    
    ann(n:n+length(ann_1)-j) = ann_1(j:end);
    n = n + length(ann_1) - j + 1;
    ann(n:n+length(ann_2)-k) = ann_2(k:end);
    n = n + length(ann_2) - k + 1;
    ann(n:end) = [];
end
