% 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 4: The fundamental theorem of linear algebra.

% 1. Consider the matrix U:

U = [2 0 2;3 1 2;1 1 0];

% Find the basis for the left null space of U. How would you represent this
% subspace of R3?

LN = null(U')
LN'*U
C = orth(U)
LN'*C


% 2. Calculate the A x B product using a columns x rows multiplication:

A = [1 1 2;0 2 1;2 1 2];
B = [0 1 1;2 1 1;0 1 2];

% In particular, show the individual column x row products that are added
% up to form the final matrix.

A*B
A(:,1)*B(1,:)+A(:,2)*B(2,:)+A(:,3)*B(3,:)

% 3. Consider the linear transformation represented by the matrix:
%
A = [-6 2 0 4;-15 5 0 10;-7 3 1 3;...
     -1 1 1 -1;1 1 2 -4;6 0 3 -9] 
rank(A)

% Let's identify the independent columns. We can accomplish this by doing
% Gaussian elimination on A until we achieve the 'row reduced echelon form'
% of the matrix.

[R,jb] = rref(A)

% a. Can you identify one possible choice of independent rows? Use the following
% two commands, and explain the result obtained. 
At = A'
[Rt,jbt] = rref(At)

% b. find a basis for the four subspaces of A.
C = A(:,jb)
R = A(jbt,:)'
N = null(A,'r')
LN = null(A','r')
% N_orth = null(A)
% LN_orth = null(A')
R'*N
C'*LN
LN'*C

% 4. Vector u1 = [3 7 5] is the 1st basis vector of an orthogonal basis U
% in R3. Using what you have learnt so far about the fundamental spaces of
% a matrix and about rotations:

% 1. Find two other possible vectors of an orthogonal basis U. 
% 2. Display the standard basis V and the new basis U.
 
u1 = [3 7 5]
u1 = u1/norm(u1)
rank(u1)

[u23] = null(u1)
U = [u1' u23]
show_frame('new',eye(3))
show_frame('add',U,'--r','--g','--b','u1','u2','u3')

% c. Find the transformation matrix Q that would allow the representation
% of any vector in the rotated frame U.

Q = inv(U)


% 5. Consider the matrix equation A*x = b for:

A = [-6 2 0 4;-15 5 0 10;-7 3 1 3;...
     -1 1 1 -1;1 1 2 -4;6 0 3 -9]
b = [14 35 14 0 -7 -21]'

r = rank(A)

% a. Certainly one possible solution for x is [1 2 3 4]. Verify that this
% is true.
x = [1 2 3 4]'
A*x

% b. Using mldivide
x = A\b


% c. Using Gauss-Jordan
[R,jb] = rref([A b])
x1 = [-3.5 -3.5 0 0]'
 
% d. Using the null space of A can you find two other exact solutions? 

N_a = null(A,'r')
A*(x1+N_a(:,1))
A*(x1+N_a(:,2))

% e. How many exact solutions do you think exist for this system of linear
% equations?

% f. Is it possible to find a solution that would have the shortest possible
% length? Notice that using the normal equation is not an option here
% because A'A is not invertible

x = inv(A'*A)*A'*b
x = (A'*A)\A'*b

% However, we can find the independent columns of A
[R,jb] = rref([A])
A_ind = A(:,jb)
[R,jb] = rref([A_ind b])
x = inv(A_ind'*A_ind)*A_ind'*b
x = A_ind\b

% Using the pseudoinverse instead
x = pinv(A)*b


% 6. Consider the following transformation matrix:
clear,clc
% p = [1 1 1]';
% theta = pi
% [Q] = Q_from_p_and_theta_alt(p,theta)
% [Q_alt] = Q_from_p_and_theta(p,theta)

Q = [0 1 0;
     0 0 1;
     1 0 0];
 
% Any transformation matrix can be represented as a rotation around a
% single axis p by an angle ?.  Use the function 'p_and_theta_from_Q.m' to
% identify this axis and angle, and then the function 'Q_from_p_and_theta'
% to calculate back the original transformation matrix. You can also try
% the original algorithm developed by Euler to carry out these operations:
% they are provided as the functions 'p_and_theta_from_Q_alt.m' and
% 'Q_from_p_and_theta_alt'.
 
[p,theta] = p_and_theta_from_Q(Q)
[Q_bk] = Q_from_p_and_theta(p,theta)

[p_alt,theta_alt] = p_and_theta_from_Q_alt(Q)
[Q_bk_alt] = Q_from_p_and_theta_alt(p_alt,theta_alt)

% a. What type of rotation is this?
% Q = [0 1 0;
%      0 0 1;
%      1 0 0];
% The columns of Q are the coordinates of V in the U basis, and the rows of
% Q are the coordinates of U in V basis, therefore:

theta_deg = 180*theta/pi


% b. Be ready to discuss the functions you have used.
% c. Is there another name for the transformation Q?

% Q is a 'permutation' matrix.

