% Homework 3 for Deeplearning 2018 Fall
%
% Based on the current code, achieve any of the following todo List:
%
% Todo List
%   * Network layer
%       - [] linear function: use `W' * X + b` instead of `W' * X`
%       - [] activation function: tanh
%       - [] activation function: relu
%       - [] loss function: mae(mean-absolute-error)
%            compared with mse(mean-squared-error), it takes 1-norm instead of 2-norm
%       - [] (BONUS) accelaration: rewrite the for-loop using matrix operation in @backwardLoss
%            This requires some knowledge on vectorization, see
%            <https://www.mathworks.com/help/matlab/matlab_prog/vectorization.html>
%   * optimization
%       - [] gradient clipping
%            this solves gradient explosion problem.
%       - [] (BONUS) minibatch SGD
%            this makes it practical to training a large dataset
%            Requirement: 
%               * the batchsize should be controllable: there should be an parameter "BatchSize" in trainingOptions
%       - [] SGD with momentum
%            This is the default optimizer used in DeepLearning
%            see <http://cs231n.github.io/neural-networks-3/#sgd>
%       - [] Adam
%            Most researchers nowadays use Adam by default, instead of SGD with momentum
%            see <http://cs231n.github.io/neural-networks-3/#ada>
%       - [] Adagrad
%            see <http://cs231n.github.io/neural-networks-3/#ada>
%       - [] RMSprop
%            see <http://cs231n.github.io/neural-networks-3/#ada>
%   * Dataset
%       - [] (BONUS) iris dataset: write the dataset reader
%            see <https://archive.ics.uci.edu/ml/datasets/iris>
%            There are two ways of representing a class:
%               * Gray Code:
%                   use 1, 2, 3 to represent class 1, class 2, and class 3, respectively
%               * One Hot:
%                   use (1,0,0), (0,1,0),(0,0,1) to represent class 1, class 2, and class 3, respectively
%            For simplicity, we use the Gray Code here.
%       - [] (BONUS) iris dataset: train
%
% Dec 8, 2018 by
% Johnny Chen <johnnychen94@hotmail.com>


%%%%% BEGIN of demo (template) %%%%%
%----- BEGIN of Description -----%
% This demo tests:
% - Activation function: sigmoid
%----- END of Description -----%

[X,Y] = readData('toy');
W = initializeNetwork(size(X,1),size(Y,1));
trainingOptions = struct(...
    'InitialLearnRate', 1, ...
    'LearnRateDropPeriod', 7000, ...
    'MaxEpoch',10000, ...
    'LearnRateDropFactor', 0.1, ...
    'Verbose', true, ...
    'VerboseFrequency', 50);
W = trainNetwork(X, Y, W, trainingOptions);
Y_hat = predict(X,W);

% show results
Y_hat
n = size(X,2);
plot(1:n, Y,'ro',...
    1:n, Y_hat,'b+');
xlabel('x');ylabel('y');title('prediction');
%%%%% END of demo (template) %%%%%



%%% Network Definition

function W = trainNetwork(X, Y, W, option)
   InitialLearnRate    = option.InitialLearnRate;
   LearnRateDropPeriod = option.LearnRateDropPeriod;
   LearnRateDropFactor = option.LearnRateDropFactor;
   MaxEpoch            = option.MaxEpoch;
   Verbose             = option.Verbose;
   VerboseFrequency    = option.VerboseFrequency;
   
   if Verbose     
       hfig = figure;

       subplot(1,2,1);
       hfig_loss = animatedline();
       xlabel('iter');ylabel('Mean Squared Error');
       title('Gradient Descent Progress')
   end
    
   LearnRates = LearnRateScheduler(InitialLearnRate, MaxEpoch, LearnRateDropPeriod, LearnRateDropFactor);
   NumSample = size(X,2);
   for epoch = 1:MaxEpoch
       % get training data of current iteration
       cur_lr = LearnRates(epoch);
       cur_X  = X;
       cur_Y  = Y;
       
       % training 
       [loss, cur_Y_hat] = forwardLoss(cur_X,cur_Y,W);
       dLdW = backwardLoss(cur_X,cur_Y,cur_Y_hat,W);
       W = update(W, dLdW, cur_lr);
       
       % show training progress
       if Verbose && mod(epoch, VerboseFrequency) == 0
           addpoints(hfig_loss, epoch, loss);
           
           figure(hfig);subplot(1,2,2);
           plot(1:NumSample,cur_Y,'ro',...
               1:NumSample,cur_Y_hat,'b+');
           xlabel('x');ylabel('y');title('prediction');
           
           drawnow;
           
           fprintf(sprintf("Epoch: %d\tLoss: %f\tLearning Rate: %f\n", epoch, loss, cur_lr));
       end
   end
   
   if Verbose
       close(hfig);
   end
end

%%% Weight Initialization

function W = initializeNetwork(inLength,outLength)
% initNet returns initial weights for the so-called "network"
    W = 1e-3 .* rand(inLength,outLength);
end

%%% LearnRateScheduler(s)

function lr = LearnRateScheduler(InitialLearnRate, MaxEpoch, LearnRateDropPeriod, LearnRateDropFactor)
% LearnRateScheduler generates a sequence of learning rate using step decay method
    lr = repmat(InitialLearnRate,MaxEpoch,1);
    n_drop = ceil(MaxEpoch/LearnRateDropPeriod);
    for i = 2:n_drop
        pos_l = (i-1)*LearnRateDropPeriod + 1;
        pos_r = min(i*LearnRateDropPeriod, MaxEpoch);
        lr(pos_l:pos_r) = lr(pos_l-1) * LearnRateDropFactor;
    end
end

%%% forward propogation + backward propogation

function Y = forward(X, W)
% forward outputs the network result with given input X
    V = W' * X;
    Y = logistic(V);
end

function Y = predict(X,W)
% predict predicts the label of X with given network weight W 
    Y = forward(X,W);
    Y = round(Y); % post-processing
end

function [loss, Y_hat] = forwardLoss(X, Y, W)
% forwardloss uses Mean Squared Error to calculate training loss
    Y_hat = forward(X,W);
    loss = 1/2 * mean((Y_hat - Y).^2); % mse
end

function dLdW = backwardLoss(X,Y,Y_hat,W)
% backwardloss calculates the gradient
    dphidx = @(x) exp(x)./(1 + exp(x))^2; % i.e., Sigmoid(x)(1-Sigmoid(x))
    dWTxdW = @(x) x;
    
    n = size(X,2);
    dLdW = zeros(size(W));
    for i = 1:n % this for-loop can be acceralated
        x     = X(:,i);
        y     = Y(:,i);
        y_hat = Y_hat(:,i);
        
        dLdW_i = dphidx(W' * x) * dWTxdW(x) * (y_hat - y);
        dLdW = dLdW + dLdW_i;
    end
    dLdW = dLdW/n;
end

%%% optimizer(s)

function W = update(W, dLdW, lr)
% update updates network weights W
    W = update_sgd(W, dLdW, lr);
end

function W = update_sgd(W, dLdW, lr)
% update_sgd updates network weights W using gradient descent
    W = W - lr * dLdW;
end

function W = update_sgdm(W, dLdW, lr, mom)
% update_sgdm updates network weights W using gradient descent with momentum
    error('Not implemented');
end

function [W, cache] = update_adam(W, dLdW, beta1, beta2, cache)
% update_adam updates network weights W using ADAM method
    error('Not implemented');
end

%%% activation function(s) %%%

function y = logistic(x)
    y = 1 ./ (1 + exp(-x));
end

function y = tanh(x)
    y = (exp(2*x)-1)./(exp(2*x)+1);
end

function y = relu(x)
    y = max(0,x);
end

%%% dataset reader %%%

function [X,Y] = readData(dataset)
% readData returns training data
% 
% Input:
%   dataset : (Required) String
%       Current supported datasets are:
%           'toy'    : toy dataset for hw3
%           'boston' : The Boston Housing Dataset
%           'iris'   : Iris Data Set %TODO: Not suppored yet
% Return:
%    X : Array
%       input data
%    Y : Array
%       label/groundtruth of input data
% Note:
%   * By convention, the last dimension indicates the index of sample, 
%     i.e., X(:,1) stands for the first sample.
    dataReader=struct(...
        'toy',@() readData_toy(),...
        'boston',@() readData_boston(),...
        'iris',@() readData_iris()); %#ok<NASGU>
    
    readData_ = eval(['dataReader.',lower(dataset)]);
    [X,Y] = readData_();
end

function [X,Y] = readData_toy()
% readData_toy returns a toy dataset for homework 3
    X = [ 0 0 1 1; 0 1 0 1; 1 1 1 1];
    Y = [ 0 0 1 1];
end

function [X,Y] = readData_boston()
% readData_Boston reads boston housing dataset
    dataurl    = 'http://lib.stat.cmu.edu/datasets/boston';
    datadir    = 'data';
    filename   = fullfile(datadir,'BostonHousing.txt');
    
    if ~isfolder(datadir)
        mkdir(datadir);
    end
    if ~isfile(filename)
        websave(filename, dataurl);
    end

    assert(isfile(filename), [filename,' not found']);
    try 
        fileID = fopen(filename);
        data = parsedata(fileID);
        X = data(:,1:end-1)';
        Y = data(:,end)';
        fclose(fileID);
    catch err
        fclose(fileID);
        rethrow(err)
    end
    return
    
    function data = parsedata(fileID)
        row_offset  = 23; % records start from line 23
        num_records = 506;
        num_columns = 14;

        % pass description
        for i = 1: row_offset - 1
            fgetl(fileID);
        end
        
        data = zeros(num_records,num_columns);
        i = 1;
        while i <= 506
            % each record contains two lines of file
            part1 = fgetl(fileID);
            part2 = fgetl(fileID);

            record = textscan([part1,part2],'%f'); % cell
            data(i,:) = record{:};
            i = i + 1;
        end
    end
end

function [X,Y] = readData_iris()
    error('Not implemented');
end