function newbinary(train_data, test_data, options)

    if (nargin ~= 3)
        if (nargin == 2)
            options = '';
        else
            fprintf('Incorrect number of input parameters\n');
            return;
        end
    end

    tmp = strfind(options, '-R');
    if isempty(tmp) == 0
        fprintf(1, '-R is not supported\n');
        exit;
    end
    
    bias = 0;
    options_split = strsplit(options, ' ');
    tmp = strfind(options_split, '-B');
    for i = 1:size(tmp, 2),
        if tmp{i} == 1,
            bias = str2num(options_split{i+1});
            options = [strjoin(options_split(1,1:i-1)) ' ' strjoin(options_split(1,i+2:end))];
            break;
        end
    end
    
    if isequal(class(train_data), 'char')
        [y, x, map] = libsvmread_ml(train_data);
        [ty, tx, tmap] = libsvmread_ml(test_data);
    else
        x = train_data.x;
        y = train_data.y;
        map = train_data.map;
        tx = test_data.x;
        ty = test_data.y;
        tmap = test_data.map;
    end

    x = x';
    tx = tx';

    % eliminate instances that have duplicated label
    y = full(y);
    ty = full(ty);

    t0all = clock;

    rand('seed', 1);

    % fill x
    probk = size(x, 1);
    testk = size(tx, 1);
    k = max(probk, testk);
    x = [x; zeros(k-probk, size(y, 1))];
    tx = [tx; zeros(k-testk, size(ty, 1))];

    if bias > 0,
        x=[x; ones(1, size(x,2))*bias];
        tx=[tx; ones(1, size(tx,2))*bias];
    end
    
    d = size(y, 2);

    w = zeros(size(x, 1), d);
    b = zeros(d, 1);

    % obtain (w,b) for each label
    for i = 1:d
        fprintf(1, 'INFO: training label %d (i %d)...\n', map(i), i);
        t0 = clock;
        [w(:, i), b(i)] = train_one_label(2*y(:, i)-1, x, options);
        fprintf(1, 'INFO: b = %g\n', b(i));
        fprintf(1, 'INFO: Training Time: %.6f seconds\n', etime(clock, t0));
    end

    fprintf(1, 'INFO: predicting...\n');

    fprintf(1, 'INFO: calculate measurement\n');

    % we calculate measurement for 1+train(101) (the first row in the paper)

    F = 0;
    tp_sum = 0;
    fp_sum = 0;
    fn_sum = 0;
    tn_sum = 0;

    for i = 1:d
        % deal with this label
        label = map(i);

        test_idx = find(tmap == label);

        if isempty(test_idx)
            orig = zeros(size(ty, 1), 1);
        else
            orig = ty(:, test_idx);
        end

        pre = ((w(:,i)'*tx)' + b(i)) > 0;

        tp = full(sum(orig == +1 & pre == +1));
        fn = full(sum(orig == +1 & pre == 0));

        tn = full(sum(orig == 0 & pre == 0));
        fp = full(sum(orig == 0 & pre == +1));

        this_F = 0;
        if tp ~= 0 || fp ~= 0 || fn ~= 0
            this_F = (2*tp) / (2*tp + fp + fn);
        end
        F = F + this_F;

        tp_sum = tp_sum + tp;
        fp_sum = fp_sum + fp;
        fn_sum = fn_sum + fn;
        tn_sum = tn_sum + tn;

        fprintf(1, 'INFO: label %d:\tF %.6f\ttp %d\tfp %d\ttn %d\tfn %d\n', label, this_F, tp, fp, tn, fn);
    end

    micro = (2*tp_sum) / (2*tp_sum + fp_sum + fn_sum);
    macro = F / d;

    fprintf(1, 'INFO: tp_sum %d fp_sum %d fn_sum %d tn_sum %d\n', tp_sum, fp_sum, fn_sum, tn_sum);
    fprintf(1, 'INFO: microaverage: %.6f\n', micro);
    fprintf(1, 'INFO: macroaverage: %.6f\n', macro);

    fprintf(1, 'INFO: Total Time: %.6f seconds\n', etime(clock, t0all));

end

function [w, b] = train_one_label(proby, x, options)

    fbr_list = [0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8];

    nr_fold = 3;

    l = size(proby, 1);

    perm = randperm(l)';

    f_list = zeros(size(fbr_list));

    % cv
    for fold = 1:nr_fold
        train_id = [1:floor((fold-1)*l/nr_fold) floor(fold*l/nr_fold)+1:l]';
        valid_id = [floor((fold-1)*l/nr_fold)+1:floor(fold*l/nr_fold)]';

        y = proby(perm(train_id));
        validy = proby(perm(valid_id));

        % scutfbr
        [scutfbr_w, scutfbr_b_list] = scutfbr(y, x(:, perm(train_id)), fbr_list, options);
        wTx = (scutfbr_w'*x(:, perm(valid_id)))'; % +b

        for i = 1:size(fbr_list, 2)
            F = fmeasure(validy, 2*(wTx > -scutfbr_b_list(i))-1);
            f_list(i) = f_list(i) + F;
        end
    end

    best_fbr = fbr_list(find(f_list == max(f_list), 1, 'last'));
    if max(f_list) == 0
        best_fbr = min(fbr_list);
        fprintf(1, 'INFO: train_one_label: F all 0\n');
    end

    % final model
    [w, b] = scutfbr(proby, x, best_fbr, options);
    fprintf(1, 'INFO: train_one_label: best_fbr %.1f\n', best_fbr);

end

function [w, b_list] = scutfbr(proby, probx, fbr_list, options)
    b_list = zeros(size(fbr_list, 2));

    nr_fold = 3;

    l = size(proby, 1);

    perm = randperm(l)';

    % cv
    for fold = 1:nr_fold
        train_id = [1:floor((fold-1)*l/nr_fold) floor(fold*l/nr_fold)+1:l]';
        valid_id = [floor((fold-1)*l/nr_fold)+1:floor(fold*l/nr_fold)]';

        y = proby(perm(train_id));
        validy = proby(perm(valid_id));

        % scut
        scut_w = do_train(y, probx(:, perm(train_id)), options);

        wTx = (scut_w'*probx(:, perm(valid_id)))';
        scut_b = 0;

        start_F = fmeasure(validy, 2*(wTx > -scut_b)-1);

        [sorted_wTx, sorted_wTx_index] = sort(wTx, 1, 'ascend'); 

        tp = sum(validy == 1); fp = size(validy,1)-tp; fn = 0;
        cut = -1;
        best_F = (2*tp) / (2*tp + fp + fn);
        for i = 1:size(validy, 1)
            if validy(sorted_wTx_index(i)) == -1,
                fp = fp -1;
            else
                tp = tp -1; fn = fn + 1;
            end
            F = (2*tp) / (2*tp + fp + fn);
            
            if F >= best_F
                best_F = F;
                cut = i;
            end
        end

        % modify b
        if best_F > start_F
            if cut == -1 % i.e., all +1
                val = sorted_wTx(1);
                scut_b = - (val - eps(val)); % predict all +1
            elseif cut == size(validy, 1)
                val = sorted_wTx(end);
                scut_b = - (val + eps(val));
            else
                scut_b = - (sorted_wTx(cut) + sorted_wTx(cut+1))/2;
            end
        end

        F = fmeasure(validy, 2*(wTx > -scut_b)-1);

        % compare to fbr_list
        for i = 1:size(fbr_list, 2)
            if F > fbr_list(i)
                b_list(i) = b_list(i) + scut_b;
            else
                b_list(i) = b_list(i) - max(wTx);
            end
        end
    end

    % final model
    b_list = b_list / nr_fold;
    w = do_train(proby, probx, options);

end

function w = do_train(y, x, options);

    model = train(y, x, options, 'col');

    w=model.w'; % n by nr_class

    % Liblinear flips +1/-1 labels so +1 is always the first label,
    % but not if all labels are -1.
    % For our usage, we need +1 to always be the first label,
    % so the check is necessary.
    if model.Label(1) == -1
        fprintf(1, 'warning: Label(1) = -1');
        w = -w;
    end

end

function F = fmeasure(original, predict)

    F = 0;

    tp = full(sum(original == +1 & predict == +1));
    fn = full(sum(original == +1 & predict == -1));
    fp = full(sum(original == -1 & predict == +1));

    if tp ~= 0 || fp ~= 0 || fn ~= 0
        F = (2*tp) / (2*tp + fp + fn);
    end

end
