% Copyright (c) 2017, Domenico L. Gatti
% All rights reserved.
% 
% Redistribution and use in source and binary forms, with or without 
% modification, are permitted provided that the following conditions are 
% met:
% 
%     * Redistributions of source code must retain the above copyright 
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright 
%       notice, this list of conditions and the following disclaimer in 
%       the documentation and/or other materials provided with the 
%       distribution
%       
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
% IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
% CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
% EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
% PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
% PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
% NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%
%% General dependencies
% We always start from the CODE directory and we add to the path
% subdirectories containing various tools described in the book chapters.
addpath(genpath('../GENERAL_SCRIPTS_FUNCTIONS'));
% addpath(genpath('../DATABASE'));
% addpath(genpath('../TOOLBOXES'));

%% CHAPTER 9: Eigenvalues and eigenvectors 

% Problem 2. find the characteristic polynomial of the matrix A:

A = [3 1;1 3]
syms x
charpoly(sym(A),x)
ch_poly = charpoly(sym(A))
D = roots(double(ch_poly))
S1 = null(A-D(1)*eye(2)) , S2 = null(A-D(2)*eye(2))

A = [1 -1 0;-1 2 -1;0 -1 1]
syms x
charpoly(sym(A),x)
ch_poly = charpoly(sym(A))
D = roots(double(ch_poly))
S1 = null(A-D(1)*eye(3)) , S2 = null(A-D(2)*eye(3)) , S3 = null(A-D(3)*eye(3))

A = [1 -1 0 2;-1 2 -1 1;0 -1 1 -2;1 -2 -1 0]
rank(A)
syms x
charpoly(sym(A),x)
ch_poly = charpoly(sym(A))
D = roots(double(ch_poly))
S1 = null(A-D(1)*eye(4)) , S2 = null(A-D(2)*eye(4)) ,
S3 = null(A-D(3)*eye(4)) , S4 = null(A-D(4)*eye(4)) ,

% Find the nullspace of A using Gaussian elimination
A = [1 1 2 3 1 1 1;2 2 8 10 3 2 1;3 3 10 13 1 2 3]
% A = [1 1 2 3 1 1 1;2 2 8 10 3 2 1;3 3 6 9 3 3 3;2 2 4 6 2 2 2]
[m,n] = size(A);
r = rank(A);
x_ind = [1:n];
[~,pivot_ind] = rref(A);
free_ind = setdiff(x_ind,pivot_ind);
% npivots = length(pivot_ind);
nfree = length(free_ind);
Ap = -A(:,pivot_ind);
Af = A(:,free_ind);
X = Ap\Af;
I = eye(nfree);
NA = zeros(n,n-r);
% for i = 1:nfree
%     NA(pivot_ind,i) = X(:,i);
%     NA(free_ind,i) = I(:,i);
% end
NA(pivot_ind,:) = X;
NA(free_ind,:) = I;
NA

% NA = null(sym(A))
% NA = null(A)

% Gram-Schmidt
s1 = NA(:,1);s2 = NA(:,2);s3 = NA(:,3);s4 = NA(:,4);
s2 = s2-s1*inv(s1'*s1)*s1'*s2;
s3 = s3-[s1 s2]*inv([s1 s2]'*[s1 s2])*[s1 s2]'*s3;
s4 = s4-[s1 s2 s3]*inv([s1 s2 s3]'*[s1 s2 s3])*[s1 s2 s3]'*s4;

s1 = s1/norm(s1);
s2 = s2/norm(s2);
s3 = s3/norm(s3);
s4 = s4/norm(s4);

oNA = [s1 s2 s3 s4]
oNA'*oNA
oNA*oNA'
A*oNA

% Problem 3. Find the eigenvalues, the eigenvectors, and the inverse of A:

A1 = eye(4)*5
A2 = -ones(4)
Aprod = A1*A2
Asum = A1 + A2
A = sym(Asum)
[S,D] = eig(A)
inv(A)

% Problem 4. True or false: if the columns of S (eigenvectors of A) are
% linearly independent, then:
% a) A is invertible
% b) A is diagonalizable
% c) S is invertible
% d) S is diagonalizable

A = randi(20,[5 5])
A = A'*A
[S,D] = eig(A)
rank(S)
[S2,D2] = eig(S)
rank(S2)

% True = b,c,d

% Problem 6. Which of the four subspace of a projection matrix contains
% eigenvectors with lambda = 1? Which subspace contains eigenvectors vith
% lambda = 0? Can you say if a projection matrix can be diagonalized?

A = [3 6 7;2 4 3;1 2 1]
rank(A)
[Sl,D] = eig(A)
[Sr,D] = eig(A')

if rank(A'*A)<size(A'*A,1)
    disp('The projection matrix does not exist')
else
Pa = A/(A'*A)*A'
[Sp,Dp] = eig(Pa)
end

% but:
A = randi(10,[3 2])
rank(A)

if rank(A'*A)<size(A'*A,1)
    disp('The projection matrix does not exist')
else
Pa = A/(A'*A)*A'
[Sp,Dp] = eig(Pa)
end

Sp*Dp/Sp

