function [tp, fp, tn, fn, split] = roc(label, confidence, valid_splits)
%ROC ROC curves for binary classification
%
% [tp, fp, tn, fn, split] = cfmat(label, conf, valid_splits)
%
% Input:
% - label         Labels from {0, 1}, any size.
% - confidence    Confidence for positive class, corresponding to label.
% - valid_splits  Consider only valid splits?
%
% Output:
% - tp     Number of true positives corresponding to the split.
% - fp     Number of false positives.
% - tn     Number of true negatives.
% - fn     Number of false negatives.
% - split  Corresponding confidence split thresholds.
%

if nargin < 3 || isempty(valid_splits)
    valid_splits = true;
end
assert(isscalar(valid_splits));
valid_splits = logical(valid_splits);

assert(numel(label) == numel(confidence));

% Convert labels to binary 0/1 labels.
label = logical(label(:));
confidence = confidence(:);

n = numel(label);
% Sort confidence for positive class.
[sorted_conf, ix] = sort(confidence);
% Reorder labels according to sorted confidence.
label = label(ix);
% Compute FN as the number of examples of given label up to given confidence,
% i.e., the decision boundary is just above current confidence.
% Extra first row correspond to all decisions being positive.
fn = [0; cumsum(label)];
tn = [0; cumsum(~label)];
neg = tn(end);
pos = fn(end);
fp = neg - tn;
tp = pos - fn;
assert(all(n == fn + tn + tp + fp));

if ~valid_splits
    split = [sorted_conf(1) - eps(sorted_conf(1)); ...
             (sorted_conf(1:end-1) + sorted_conf(2:end)) / 2; ...
             sorted_conf(end) + eps(sorted_conf(end))];
    return;
end

% Split only along unique confidence values.
[unique_conf, ia] = unique(sorted_conf);
% Create decision boundary before each unique confidence (first occurence).
% One is not subtracted from ia due to the extra row added to fn, tn, fp,
% and tp but not to sorted_conf.
% Add an extra row for all decisions being negative.
tp = tp([ia; end]);
fp = fp([ia; end]);
tn = tn([ia; end]);
fn = fn([ia; end]);
assert(all(n == fn + tn + tp + fp));
split = [unique_conf(1) - eps(unique_conf(1)); ...
         (unique_conf(1:end-1) + unique_conf(2:end)) / 2; ...
         unique_conf(end) + eps(unique_conf(end))];
end
