%% Applications of Eigen Values and Vectors: Part 1
%
% Calendar: Refer this link for due dates                                                                     NOTE: This part deals with Ex 1,2,3.
%
% Objective: In this part of the lab we will see how eigen values and vectors can be used to find the principal components of a dataset. To understand PCA we need to understand the following statistical terms: standard deviation, variance and covariance.
% Standard Deviation(SD):  The average distance from the mean of the data set to a point. The way to calculate it is to compute the squares of the distance from each data point to the mean of the set, add them all up, divide by n-1, and take the positive square root.
% \mathit{\mathbf{s}}=\sqrt{\;\frac{\sum_{\mathit{\mathbf{i}}=1}^{\mathit{\mathbf{n}}} {\left({\mathit{\mathbf{X}}}_{\mathit{\mathbf{i}}\;} -\bar{\mathit{\mathbf{X}}} \right)}^2 }{\mathit{\mathbf{n}}-1}}
% Variance (V): Variance is another measure of the spread of data in a data set. In fact it is almost identical to the standard deviation. The formula is this:
%  {\mathit{\mathbf{s}}}^2 = \frac{\sum_{\mathit{\mathbf{i}}=1}^{\mathit{\mathbf{n}}} {\left({\mathit{\mathbf{X}}}_{\mathit{\mathbf{i}}\;} -\bar{\mathit{\mathbf{X}}} \right)}^2 }{\mathit{\mathbf{n}}-1}
% Covariance(C): Standard deviation and variance only operate on 1 dimension, so that you could only calculate the standard deviation for each dimension of the data set independently of the other dimensions. However, it is useful to have a similar measure to find out how much the dimensions vary from the mean with respect to each other. Covariance is such a measure. Covariance is always measured between 2 dimensions. If you calculate the covariance between one dimension and itself, you get the variance. So, if you had a 3-dimensional data set (x,y,z), then you could measure the covariance between the x and y dimensions, the x and z dimensions, and the y and z dimensions. Measuring the covariance between x and x, or y and y, or z and z would give you the variance of the x, y and z dimensions respectively.
% \textrm{Cov}\left(\mathit{\mathbf{x}},\mathit{\mathbf{y}}\right)= \frac{\sum_{\mathit{\mathbf{i}}=1}^{\mathit{\mathbf{n}}} {\left({\mathit{\mathbf{X}}}_{\mathit{\mathbf{i}}\;} -\bar{\mathit{\mathbf{X}}} \right)\left({\mathit{\mathbf{Y}}}_{\mathit{\mathbf{i}}} -\bar{\mathit{\mathbf{Y}}} \right)} }{\mathit{\mathbf{n}}-1}
% The covariance matrix shows the relationship between different features(i.e. columns) in the dataset. For a three dimentional dataset the covariance matrix would look something like this:
% \textrm{Cov}\left(\mathit{\mathbf{x}},\mathit{\mathbf{y}},\mathit{\mathbf{z}}\right)\left\lbrack \begin{array}{c}
% \textrm{cov}\left(\mathit{\mathbf{x}},\mathit{\mathbf{x}}\right)\;\;\;\;\;\textrm{cov}\left(\mathit{\mathbf{x}},\mathit{\mathbf{y}}\right)\;\;\;\;\;\textrm{cov}\left(\mathit{\mathbf{x}},\mathit{\mathbf{z}}\right)\\
% \textrm{cov}\left(\mathit{\mathbf{y}},\mathit{\mathbf{x}}\right)\;\;\;\;\;\textrm{cov}\left(\mathit{\mathbf{y}},\mathit{\mathbf{y}}\right)\;\;\;\;\;\;\textrm{cov}\left(\mathit{\mathbf{y}},\mathit{\mathbf{z}}\right)\\
% \textrm{cov}\left(\mathit{\mathbf{z}},\mathit{\mathbf{x}}\right)\;\;\;\;\;\;\textrm{cov}\left(\mathit{\mathbf{z}},\mathit{\mathbf{y}}\right)\;\;\;\;\;\;\;\textrm{cov}\left(\mathit{\mathbf{z}},\mathit{\mathbf{z}}\right)
% \end{array}\right\rbrack
% If the value is positive, then that indicates that both dimensions increase together. If the value is negative, then as one dimension increases, the other decreases. In the last case, if the covariance is zero, it indicates that the two dimensions are independent of each other.
%% Exercise 1: Eigen Values and Vectors
%
% The eigenvalues and vectors describe the essense of how a matrix transforms a vector. To see the effect we can generate a set of vectors of unit length representing all directions:

t = 0:0.1:2*pi+0.1;
x = [cos(t); sin(t)];
plot(x(1,:),x(2,:),'LineWidth',1.5)
axis equal

% Let's Transform the vectors using some matrix A and see how it effects the vectors:

A = [1 .2
     .2 1];
y = A*x
plot(y(1,:),y(2,:),'LineWidth',1.5)
axis equal

% Using the eigenvectors and values we can compactly capture the effect of the matrix, and draw the transform. The transform here means the direction in which the vectors were affeced the most. You can plot the major and minor axis of this ellipse to get the visualize the magnitude of change.

[X,D] = eig(A)
plot(y(1,:),y(2,:),'LineWidth',1.5)
axis equal
hold on
plot(D(1,1)*[0 X(1,1)],D(1,1)*[0 X(2,1)],'g','LineWidth',1.5)
plot(D(2,2)*[0 X(1,2)],D(2,2)*[0 X(2,2)],'r','LineWidth',1.5)
hold off

% Try a few of these transforms on your own (no need to turn in), to get an idea of the effect of different matrices. Remember that (symmetric) positive definite matrices will give you real ON eigenvectors. Unsymmetric matrices may give complex results. To think about: In the above example we only rotated and scaled. More generally we can also add translation \mathit{\mathbf{y}}=\mathit{\mathbf{A}}\ldotp \mathit{\mathbf{x}}+\mathit{\mathbf{t}}. You might try on your own to make a routine which can plot any ellipse given its parametric form (as in Heath Computer problem 3.5). A more common problem in practice is that we have data that we want to somehow characterize. The data may have significant variation along some axis and less along others. A common model is that of a Gaussian distribution. If there is little variation along a particular axis, that dimension could be omitted altogether.
%

%You can create your own Gaussian Distribution using the following code:
% Gaussian distr, Box-Mueller method
n = 400
U = rand(2,n);
v = 2*pi*U(2,:);
G = ones(2,1)*sqrt(-2.*log(U(1,:))).*[cos(v);sin(v)];

plot(G(1,:),G(2,:),'.')
axis equal

% Transform
A = [1 .5
     .5 0.5];
t = [8; 1];
Y = A*G+t*ones(1,n);
plot(Y(1,:),Y(2,:),'.')
axis equal

%

[X,D] = eig(A)
D=sqrt(D);
plot(Y(1,:),Y(2,:),'.')
axis equal
hold on
plot(2*D(1,1)*[0 X(1,1)]+8,2*D(1,1)*[0 X(2,1)]+1,'r','LineWidth',1.5)
plot(2*D(2,2)*[0 X(1,2)]+8,2*D(2,2)*[0 X(2,2)]+1,'g','LineWidth',1.5)
hold off

%
%% Principal Componet Analysis (PCA)
%
% Principal components are the orthogonal vectors which correspond to highest variance of the data. Eigen values are the variances in the direction of the corresponding eigen vector.
% When using PCA we need to standardized the data by transforming it to having mean = 0. This is important because in PCA we will be looking at how different features contribute to the spread of the data. Therefore, different features must be comparable in terms of scale. Here we show how to standardize the data:
% \textrm{For}\;\textrm{any}\;\textrm{matrix}\;{\mathit{\mathbf{A}}}_{\mathit{\mathbf{m}}*\mathit{\mathbf{n}}}
% {\mathit{\mathbf{A}}}_{\textrm{std}} =\frac{\left(\mathit{\mathbf{A}}-\bar{\mathit{\mathbf{A}}} \right)}{s}
% \textrm{Cov}\left(\mathit{\mathbf{x}},\mathit{\mathbf{y}}\right)={\mathit{\mathbf{A}}}^{\mathit{\mathbf{T}}} \mathit{\mathbf{A}}
% Here, the eigen vectors of this covariance matrix are the principal components. Can you explain why (based on the definition of Covariance stated above)?
%
% In this exercise you will use PCA as described in the lab sessions, seminars and below. The purpose is to recover the main axes of the distribution and plot the axes with length (2*sqrt of eigenvalues). You can plot the axes as in the Gaussian example above, or plot an ellipse corresponding to the axes.
% Three data sets Y1, Y2 and Y3 are provided in EllipsePoints.mat. Show how to find the major and minor axes of those.

%%%%            Enter your code here        %%%%
%Add relevant comments and answers to the questions asked


%%%%            End             %%%%

%
%%
%% Exercise 2: A moving wave image
%
% The effect of image motion can be obtained by modulating a sin/cos basis. Considering the following formula:
% \;\sin \left(\textrm{ax}+\textrm{bt}\right)=\sin \left(\textrm{ax}\right)\cos \left(\textrm{bt}\right)+\cos \left(\textrm{ax}\right)\sin \left(\textrm{bt}\right)
% where a and b are constants, x is the spatial variable and t is time. We can interpret the above formula as the modulation of a sin and cos spatial basis with time varying coefficients - cos(bt), sin(bt). The example can be extended in 2D. Download the sincos_basis.mat data. B contains the sin/cos image basis (as columns) and Y is an example set of temporal modulation coefficients (each set of coefficients is a column, hence Y(:,1) are the coefficients for time 1, Y(:,2) for time 2 etc.). Show the modulated images on the screen in quick succession resulting from the different coefficients in Y. You can use imshow, then drawnow to force immediate rendering to the screen.
% You can use the function renderim.m Ii = renderim(Y(:,i),B,imsize); imshow(Ii,[]); drawnow; pause(0.1);
% Why do you see the pattern moving even though it only has 1 frame?


load sincos_bs.mat;
figure('visible','on')
for i=1:101
Ii=renderim(Y(:,i),B,imsize);
imshow(Ii,[]);
drawnow;
pause(0.01);
end

% Construct a similar basis that has twice the frequency of the given one and show the effect of motion through modulation.
% What is the effect of changing b ?

%%%%            Enter your code here            %%%%
%Add relevant comments and answers to the questions asked


%%%%            End             %%%%

%
%%
%% Exercise 3: Small image motion
% Consider an image \mathit{\mathbf{I}}=\left(\mathit{\mathbf{u}},\mathit{\mathbf{v}}\right). We study the effect of a small shift (horizontal and vertical) using Taylor expansion:
% \mathit{\mathbf{I}}\left(\mathit{\mathbf{u}}+\textrm{du},\mathit{\mathbf{v}}+\textrm{dv}\right)=\mathit{\mathbf{I}}\left(\mathit{\mathbf{u}},\mathit{\mathbf{v}}\right)+\left\lbrack \textrm{Ix}\left(\mathit{\mathbf{u}},\mathit{\mathbf{v}}\right)\textrm{Iy}\left(\mathit{\mathbf{u}},\mathit{\mathbf{v}}\right)\right\rbrack \left\lbrack \textrm{du}\;\textrm{dv}\right\rbrack \prime +\textrm{hot}
% \textrm{GI}=\left\lbrack \textrm{Ix}\;\textrm{Iy}\right\rbrack represents the image gradient. Again we can interpret the above formula using image modulation. In this case the basis B is formed by \left\lbrack \mathit{\mathbf{I}},\;\textrm{Ix},\;\textrm{Iy}\right\rbrack and the coefficients are of the form \left\lbrack 1\;\textrm{du}\;\textrm{dv}\right\rbrack. The effect of the modulation is a small image shift. In the below we compute the derivative by taking differences in spatial coordinates. It is also possible to use a time series of images, akin this example of a bus.
% Download the immotion_basis.mat data. The data is in the same format as before. Visualize the motions produced by the given coefficients the same way as for the example above. (You can slow down the rendering speed to see the effect better by using a more finely spaced set of t-values, e.g. -3:0.2:3, appropriatetely inserted into the Y matrix of coefficients.
% Repeat the procedure for another image of your choice. To calculate the image gradient (derivatives) use the Matlab function 'gradient'. Study the effect of changing the magnitudes of the coefficients in Y. What happens if you pick larger Y values (e.g. 5:10)? Why? What type of image does this procedure work best for? Explain appealing to the Taylor expansion what happens when it doesn't work.
%

%%%%            Enter your code here            %%%%
%Add relevant comments and answers to the questions asked


%%%%            End             %%%%

%
%%
%% Application of PCA: Image denoising
%
% PCA can help seperate the signal from the noise using the pricipal components. The idea is this: any components with variance much larger than the effect of the noise should be relatively unaffected by the noise. So if you reconstruct the data using just the largest subset of principal components, you should be preferentially keeping the signal and throwing out the noise. It does not remove the noise completely unless some other filters and methods like Local Grouping Pixels are used. For more detailed and better filtering algorithm you can check "Two-stage image denoising by PCA with local pixel grouping". The MATLAB code for the following paper is also available.
% So let's begin with some handwritten digits dataset available on MATLAB.

 % These are datasets for training machine learning models but let us use
 % them for our PCA. Here only need the XTrain data as we are not going to
 % performing a classification of the data.
 
[XTrain,YTrain,anglesTrain] = digitTrain4DArrayData;

%Now let's save our Xtrain data and forget about the other data
save('digits.mat','XTrain');

% The code snippet above allows you to save the data in the your present working directory. Once you run that section you can comment the codes as you will always have the 'digits.mat' file in your workspace. Now let's start having some fun with our data.

%Let's load the data we just saved
load digits.mat;

%Let's see how the digits look like before adding noise to them.
%Since we have 5000 digits let's only check the first 56 digits
montage(XTrain(:,:,1,1:56))
%Allright, Let's make things ugly here by adding some noise to it.
noisy=imnoise(XTrain,'gaussian',0.3,0.05);

%Let's see how the data looks now
montage(noisy(:,:,1,1:56))
% Now that we have added noise to the data. Let's use the noisy data to
% perform our PCA here.

% To take principal components we need to first flatten the images into
% columns. The current size of the XTrain and noisy is (28,28,1,5000) that
% contains 5000 images of 28*28 dimentions.
% Let's make the dimentions (784,5000)
noisy_flat=reshape(noisy,784,5000);

% Now I will be using the built-in PCA function for ease. Please do not use
% it for your lab exercises
[coeff, score, latent, tsquared, explained, mu] = pca(noisy_flat,'Algorithm','SVD','centered',true,'Rows','all','NumComponents',100);

%%The above line takes only 100 components from 5000 images as it is the
%%componet with most variance. You can check the score values to check it
%%for yourself.

denoised = score * coeff' + repmat(mu, 784,1);
%Here we project the data we have on these 100 pricipal components

%%Let's check how our denoising worked.
montage(reshape(denoised(:,1:56),[28,28,1,56]))

% You can see you were able to remove most of the noise from the data. You can change the number of components and list your observations. What happens when we increase the number of components to 400? What happens when we reduce it to 20?
%

%%%%                Type your answer here               %%%%%%
%% Answer:


%%%%%%                 End                      %%%%%

%
% PCA in Collaborative Filtering: You can check this blogpost for application of PCA in Netflix and Google type recommendation systems.
%
%
%