% Boston House Price Prediction with Linear Model
%
% The Boston Housing Dataset comes from: http://lib.stat.cmu.edu/datasets/boston
%
% Two methods are implemented in this demo:
%   1. Least Square Method: 
%       pros: closed-form solution -- high accuracy
%       cons: expensive to compute the inverse -- impractical for many large scale problems
%   2. Gradient Descent:
%       pros: each iteration is inexpensive -- the only possible ways for many large scale problems
%       cons: a set of new hyperparameters on learning rate is involved
%
% Dec 1, 2018 by
% Johnny Chen <johnnychen94@hotmail.com>

verbose = true; % set true to show training/iteration progress
predict = @(X,W,b) W' * X + b;

%%%%%%%%%% BEGIN of demo %%%%%%%%%%%
[X,Y] = readdata();
[X,Y] = sortdata(X,Y);
X = normalize(X,2); % necessary for gradient descent

% groundtruth
plot(Y);

% least square method
[W,b] = solver_LS(X,Y);
hold on; plot(predict(X,W,b)); hold off;
fprintf(sprintf('Loss of Least Squared:%f\n', norm(predict(X,W,b) - Y)));

% gradient descent method
init_lr = 1e-2;
lr_drop_period = 200;
[W,b] = solver_GD(X,Y, init_lr, lr_drop_period, verbose);
hold on; plot(predict(X,W,b)); hold off;
fprintf(sprintf('Loss of Gradient Descent:%f\n', norm(predict(X,W,b) - Y)));

xlabel('n');ylabel('House Price');
legend('GroundTruth','LeastSquare','GradientDescent');
%%%%%%%%%% END of demo %%%%%%%%%%%%



function [sorted_X,sorted_Y] = sortdata(X,Y)
% sortdata sorts [X,Y] in order of Y
    [sorted_Y, I] = sort(Y);
    sorted_X = X(:,I);
end


function [W,b] = solver_LS(X,Y)
% solver_LS uses Least Square Method to analytically solve the problem
    n = size(X,2);
    
    LHS = [X * X', sum(X,2); sum(X,2)', n];
    RHS = [X * Y'; sum(Y)];
    rst = LHS\RHS; % inv(LHS)*RHS
    
    W = rst(1:end-1);
    b = rst(end);
end


function [W,b] = solver_GD(X,Y, lr, lr_drop_period, verbose)
% solver_GD uses Gradient Descent Method to iteratively solve the problem
    % params
    fps          = 1e4; % plot frequency
    rel_tol      = 1e-4; % relative tolerance
    threshold    = 1; % gradient clipping threshold
    lr_drop_rate = 1e-1; % learning rate drop rate
    
    [n_in,num_sample] = size(X);
    n_out = size(Y,1);
    % helpers
    DW    = @(W,b) X * X' * W + b * sum(X,2) - X * Y'; % gradient of W
    Db    = @(W,b) sum(X,2)' * W + num_sample * b - sum(Y,2)'; % gradient of b
    clip  = @(grad,thres) min(max(grad, -thres),thres); % gradient clipping
    mse   = @(W,b) norm(W'*X + b -Y)/num_sample;
    
    % init
    W = 1e-3 .* randn(n_in,n_out);
    b = 1e-3 .* randn(n_out);
    cur_eps = mse(W,b);
    cur_rel_eps = cur_eps;
    
    iter = 0;
    if verbose
        hfig = figure;
        h = animatedline();
        xlabel('iter');ylabel('Mean Squared Error');
        title('Gradient Descent Progress')
        addpoints(h,iter,cur_eps);
        a = tic;
    end
    
    while cur_rel_eps > rel_tol
        iter = iter + 1;
        
        % use gradient clipping to avoid gradient explosion
        W = W - clip(lr * DW(W,b), threshold);
        b = b - clip(lr * Db(W,b), threshold);
        
        pre_eps = cur_eps;
        cur_eps = mse(W,b);
        cur_rel_eps = abs(pre_eps - cur_eps) / pre_eps;
        
        if mod(iter, lr_drop_period) == 0
            lr = lr * lr_drop_rate;
        end
        
        if verbose
            addpoints(h, iter, cur_eps);
            if toc(a) > 1/fps % control refresh rate
                drawnow; a = tic;
            end
        end
    end
    
    if verbose
        close(hfig)
    end
end


function [X,Y] = readdata()
% readdata reads boston housing dataset
    dataurl    = 'http://lib.stat.cmu.edu/datasets/boston';
    filename   = 'BostonHousing.txt';
    
    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)';
    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