% Pre-defining parameters.
lambda = 0.000001; % 0.000001 for Covertype dataset; 0.01 for MNIST.
ita = 3; % 3 for polynomial-decay averaging (in the paper); 0 for simple averaging (baseline 1).
num_rounds = 2 ^ 15;

% Reading Covertype data.
raw_data = csvread('covtype.data');
[n, m] = size(raw_data); % Should be [581012, 55], second dimension being 54 attributes + 1 label.
                         % Note that in our source code the usage of m and n
                         % is different from the paper; i.e. n is the number
                         % of samples while m is the dimension of vectors.
X = raw_data(:, 1:m-1);

% Normalizing features (optional).
% X = (X - min(X)) ./ (max(X) - min(X));

y_raw = raw_data(:, m); % Labels come in [1, 2, ..., 7], which will be
                        % converted to [1, -1, -1, ..., -1] for our purpose.
twos = ones(n, 1) * 2;
y_interm = min(twos, y_raw); % [1, 2, ..., 7] => [1, 2, 2, ..., 2].
y = 3 - y_interm * 2; % [1, 2] => [1, -1].

w = zeros(m-1, 1);

% Reading MNIST data.
%{
X = csvread('MNIST_X_train.csv'); % This dataset is already normalized to [0, 1].
[n, m] = size(X); % Should be [60000, 784].
y_raw = csvread('MNIST_y_train.csv'); % Labels come in [0, 1, ..., 9], which will be
                                      % converted to [1, -1, -1, ..., -1] for our purpose.
y_interm = min(ones(n, 1), y_raw); % [0, 1, ..., 9] => [0, 1, 1, ..., 1].
y = 1 - y_interm * 2; % [0, 1] => [1, -1].
w = zeros(m, 1);
%}

w_avg = w;
t = 1;
next_t_to_print = t; % Less-printing mode for MUCH better efficiency.

out_file = fopen('out.csv', 'w'); % Preparing output file, for keeping (log2(T), log2(F)) values.

while t <= num_rounds
    %{
    % Always-printing mode. VERY SLOW.
    F = sum(max(zeros(n, 1), ones(n, 1) - y .* (X * w_avg))) / n + (w_avg' * w_avg) * lambda / 2;
    fprintf(out_file, '%f, %f\n', log(t) / log(2), log(F) / log(2)); % Outputting log2(t), log2(F).
    fprintf('%d/%d\n', t, num_rounds); % Printing progress to console.
    %}

    if t >= next_t_to_print % Less-printing mode.
        next_t_to_print = next_t_to_print * 2; % Only printing at exponents of 2.
        F = sum(max(zeros(n, 1), ones(n, 1) - y .* (X * w_avg))) / n + (w_avg' * w_avg) * lambda / 2;
        fprintf(out_file, '%f, %f\n', log(t) / log(2), log(F) / log(2)); % Outputting log2(t), log2(F).
        fprintf('%d/%d\n', t, num_rounds); % Printing progress to console.
    end

    % Starting the SGD algorithm in the paper.
    i = randi([1 n]);
    g = lambda * w - (y(i) * X(i, :) * (X(i, :) * w * y(i) <= 1))';
    w = w - (1 / (lambda * t)) * g;
    t = t + 1;

    % The paper: ita = 3; baseline 1: ita = 0;
    w_avg = (1 - (ita + 1) / (ita + t)) * w_avg + ((ita + 1) / (ita + t)) * w;

    % Baseline 2: taking last element.
    % w_avg = w;
end

fclose(out_file);