edX ML SVM with Cross-Validation

% edX ML HW 8 problem 2, version 2, with home-grown K-fold binning

% Return values:
%     minCs:        A 100x1 vector of C value with lowest Ecv for
%                   each repetition
%
%     C_Ecv_table:  N_rep x |Cvals| matrix of Ecv for each C value in
%                   each repetition.

function [minCs, C_Ecv_table] = p7hw8v2()

    Q = 2;       % polynomial degree
    K = 10;      % K-fold cross validation
    N_rep = 100; % Repeat experiment N_rep times

    Dtrain = importdata('features.train');

    % Comparing digit 1 vs digit 5 classification

    % Pull out subset of training data where digit is '1' or '5'.
    Dtrain_1vs5  = Dtrain((Dtrain(:,1) == 1) | (Dtrain(:,1) == 5),:);
    y_train      = label_digits(1, Dtrain_1vs5(:,1));

    % The C values we'll compare, from smallest to largest
    Cvals = [0.0001 0.001 0.01 0.1 1];

    % Maintain a N_rep x |Cvals| matrix of Ecv for each C value in each rep.

    C_Ecv_table = zeros(N_rep, numel(Cvals));
    C_Ecv = zeros(1, numel(Cvals));     % one row in table, pre-allocated

    % pre-allocate return value
    minCs = zeros(N_rep,1);

    for rep=1:N_rep
        fprintf('================== RUN %d of %d ===================\n', rep, N_rep);
        c_idx = 1;
        for C = Cvals
            fprintf('========\nC: %f\n', C);
            svmopts = sprintf('-s 0 -t 1 -d %d -g 1 -r 1 -c %f -h 0', Q, C);
            fprintf('*** svmopts: "%s"\n', svmopts);
            Ntrain = size(y_train,1);

            % Generate set of indices for K-fold CV
            I = randperm(Ntrain);

            Eval_sum = 0;   % Running sum of K values for computing Ecv
            for k=1:K
                % The interval [a,b] runs through the K bins used
                % for the validation set.
                a = round(Ntrain*(k-1)/K)+1;
                b = round(Ntrain*k/K);

                % Validation set, bin k of K:
                ycv_val = y_train(I(a:b));
                Dcv_val = Dtrain_1vs5(I(a:b),:);

                % Training set, eveything not in bin k:
                ycv_train = y_train;     % Copy entire training set labels.
                ycv_train(I(a:b)) = [];  % Now delete the validation bin.

                Dcv_train = Dtrain_1vs5;  % Repeat for features.
                Dcv_train(I(a:b),:) = [];

                dbg_idxs = I(a:b);
                fprintf('\n**** CV %d of %d; bin window %d to %d, indexes [%d %d %d %d ... ****\n', ...
                    k, K, a, b, dbg_idxs(1:4));

                fprintf('**** Training on %d points ****\n', numel(ycv_train));
                model = svmtrain(ycv_train, Dcv_train(:,2:3), svmopts);

                fprintf('**** Testing on %d points ****\n', numel(ycv_val));
                [pred_val]  = svmpredict( ycv_val, Dcv_val(:,2:3), model);

                % Calculate the error ourselves, as a sanity check.
                N_mismatch_val = numel(find(pred_val ~= ycv_val));
                Eval = N_mismatch_val/size(ycv_val, 1);

                Eval_sum = Eval_sum + Eval;
                fprintf('** Eval = %f **\n', Eval);
            end

            Ecv = Eval_sum/K;
            fprintf(' ** (%d) Ecv = %f for C = %f **\n', c_idx, Ecv, C);

            C_Ecv(c_idx) = Ecv;
            c_idx = c_idx + 1;
        end     % for C = Cvals
        C_Ecv_table(rep,:) = C_Ecv;

        [~, I] = min(C_Ecv);
        minCs(rep) = Cvals(I);
        fprintf('*** minCs(%d) = %f ***\n', rep, minCs(rep));
    end      % for rep=1:N_rep

    % Display output here, for convenience.
    fprintf('\n');
    tabulate(minCs)

    fprintf('\n');
    disp(sum(C_Ecv_table)/N_rep)  % Average Ecv for each C-value column.
end


% Given an Mx1 matrix of digits, return an
% Mx1 matrix y where y(i) = +1 or -1 as
% digits(i) = digit_one
function y = label_digits( digit_one, digits)

    y = repmat(-1, size(digits,1), 1);
    y( logical(digits == digit_one)) = +1;
end