%% MCHI, cv4_2015, file 1 (cv4_2015_1) - Least squares
%
% Least squares - regression
%
% Author:       Martin Isoz
% Organisation: UCT Prague
% Date:         06. 10. 2015
%
% License: This code is published under MIT License, please do not abuse
% it.
%

%% Clear the environment
close all
clear all
clc

%% Insert data
tic
x     = -3:1:3;
y = [20.3 11.5 6.2 1.9 0.1 -0.05 2.1];

m = 3;                                                                      %number of parameters
n = numel(y);                                                               %number of points

%% Prepare right hand side
xMat = repmat(x,m,1);                                                       %create matrices with input data
yMat = repmat(y,m,1);

pVec = (0:m-1)';                                                            %vector with powers
pMat = repmat(pVec,1,n);                                                    %create a matrix of it

b    = sum(yMat.*xMat.^pMat,2);                                             %prepare rhs

%% Prepare the problem matrix
xMat = reshape(repmat(xMat,m,1),m,m,n);                                     %!!3D MATRIX!!

pMat = repmat((meshgrid(pVec) + meshgrid(pVec)'),1,1,n);                    %!!3D MATRIX!!

A    = sum(xMat.^pMat,3);                                                   %prepare problem matrix

%% Solve the system (get the coefficients)
aVec = A\b;
aVec = flipud(aVec);                                                        %because of MATLAB

toc
%% Show the solution
figure('Units','Pixels','Position',[100 50 1000 700],'Visible','on');
plot(x,y,'ro','MarkerSize',10)
hold on
plot(linspace(x(1),x(end),100),polyval(aVec,linspace(x(1),x(end),100)))
title(...
    {...
        'Polynomial regression ',...
        sprintf('$y_{fit} = \\sum_{i = 0}^{%d}a_{i}x^i, \\vec{a} = [a_0,a_1,\\dots,a_%d]^{T} = %s^{T}$',m-1,m-1,strcat(mat2str(aVec,3)))...
    },...
    'FontSize',20,'FontWeight','Bold','interpreter','latex')
xlabel('$x$',...
    'FontSize',20,'FontWeight','Bold','interpreter','latex')
ylabel('$y$',...
    'FontSize',20,'FontWeight','Bold','interpreter','latex')

%% Do the same thing again for more data
%  NOTE: the degree of the polynomial is also a thing to play with
nPts = 50000;                                                               %number of points
nMag = 2;                                                                   %noise magnitude
xVec = linspace(-3,3,nPts);
% -- var A
coefs= [1;-3;2];                                                            %vector of polynomial coefs
% -- var B
% coefs= [1;-3;2;4;3;1;-2];                                                 %vector of polynomial coefs
% -- var C
% polD = 8;
% coefs= round(10*rand(1,polD+1)-10);
yVec = polyval(coefs,xVec);
yVec = yVec + nMag*rand(size(yVec))-1;                                      %add noise

m = numel(coefs);                                                           %number of parameters (degree of the polynomial)
n = numel(yVec);                                                            %number of points

tic;
xMat = repmat(xVec,m,1);                                                    %create matrices with input data
yMat = repmat(yVec,m,1);

pVec = (0:m-1)';                                                            %vector with powers
pMat = repmat(pVec,1,n);                                                    %create a matrix of it

b    = sum(yMat.*xMat.^pMat,2);                                             %prepare rhs

xMat = reshape(repmat(xMat,m,1),m,m,n);                                     %!!3D MATRIX!!

pMat = repmat((meshgrid(pVec) + meshgrid(pVec)'),1,1,n);                    %!!3D MATRIX!!

A    = sum(xMat.^pMat,3);                                                   %prepare problem matrix

aVec = A\b;
aVec = flipud(aVec);                                                        %because of MATLAB

toc

figure('Units','Pixels','Position',[100 50 1000 700],'Visible','on');
plot(xVec,yVec,'r+','MarkerSize',10)
hold on
plot(xVec,polyval(aVec,xVec));
title(...
    {...
        'Polynomial regression ',...
        sprintf('$y_{fit} = \\sum_{i = 0}^{%d}a_{i}x^i, \\vec{a} = [a_0,a_1,\\dots,a_%d]^{T} = %s^{T}$',m-1,m-1,strcat(mat2str(aVec,3))),...
        sprintf('$y_{org} = \\sum_{i = 0}^{%d}k_{i}x^i, \\vec{k} = [k_0,k_1,\\dots,k_%d]^{T} = %s^{T}$',m-1,m-1,strcat(mat2str(coefs,3)))...
    },...
    'FontSize',20,'FontWeight','Bold','interpreter','latex')
xlabel('$x$',...
    'FontSize',20,'FontWeight','Bold','interpreter','latex')
ylabel('$y$',...
    'FontSize',20,'FontWeight','Bold','interpreter','latex')

%% Do the same thing again, using the normal equations
%  NOTE: the degree of the polynomial is also a thing to play with
tic;
pVec = (0:m-1)';                                                            %vector with powers

A    = repmat(xVec',[1,m]).^repmat(pVec',[n,1]);                             %prepare matrix for normal equations

aVec = (A'*A)\(A'*yVec');                                                    %get the coefficients from normal equations
aVec = flipud(aVec);                                                        %because of MATLAB

toc

figure('Units','Pixels','Position',[100 50 1000 700],'Visible','on');
plot(xVec,yVec,'r+','MarkerSize',10)
hold on
plot(xVec,polyval(aVec,xVec));
title(...
    {...
        'Polynomial regression ',...
        sprintf('$y_{fit} = \\sum_{i = 0}^{%d}a_{i}x^i, \\vec{a} = [a_0,a_1,\\dots,a_%d]^{T} = %s^{T}$',m-1,m-1,strcat(mat2str(aVec,3))),...
        sprintf('$y_{org} = \\sum_{i = 0}^{%d}k_{i}x^i, \\vec{k} = [k_0,k_1,\\dots,k_%d]^{T} = %s^{T}$',m-1,m-1,strcat(mat2str(coefs,3)))...
    },...
    'FontSize',20,'FontWeight','Bold','interpreter','latex')
xlabel('$x$',...
    'FontSize',20,'FontWeight','Bold','interpreter','latex')
ylabel('$y$',...
    'FontSize',20,'FontWeight','Bold','interpreter','latex')