%Lab 3: SVD and projections 

%1 Verify the things about the SVD that were in the intro lecture

%Start from a blank slate
clear all
close all
clc

%Putting everything we just wrote in a loop to show that it *really* works
%Loops are - perhaps the easiest and easiest to understand - way to tell
%computers to do something repeatedly 
numTries = 1e5;
equalityTracker = 0; %Initalize as 0

for ii = 1:numTries %This opens a for loop and iterates the counter ii with the right side


M = randn(2,2); %Creating a random 2x2 matrix M by drawing from a standard normal distribution with replacement
[U, S, V] = svd(M); %This will do the SVD and give you U, S and V, with som conventions
%S will always contain positive values in order of magnitude, U and V will be adjusted accordingly, you don't have to worry about it
x = randn(2,2); %Our random 2x2 input matrix

%Proof by Matlab (numerical)
C = M*x; %Original matrix M matrix-multiplying x
D = U * S * V' * x; %This is the SVD, basically. Output of C and D should match, if Matlab did it right, and SVD is valid as a method
%We also know that the SVD should yield a sum of outer products, that
%should *also* match C (and D)
E = (S(1,1).*(U(:,1)*V(:,1)') * x ) ... %... = Ellipses allow you to continue a statement in the next line
  + (S(2,2).*(U(:,2)*V(:,2)') * x ); 

%Dealing with excessive numerical precision (and getting rid of round
%errors)
sigDig = 4; 
C = round(C,sigDig); %Rounding C to the desired number of significant digits
D = round(D,sigDig); %Same
E = round(E,sigDig); %Again

%Test for equality
C == D;
C == E;
D == E; 

% == returns whether each element matched. We want a single number
%to denote whether the whole matrix is the same. Logic: If the number of
%elements in the matrix matches the sum of the booleans. 

equalityTracker = equalityTracker + (sum(sum(C == D)) == numel(C));
equalityTracker = equalityTracker + (sum(sum(C == E)) == numel(C));
equalityTracker = equalityTracker + (sum(sum(D == E)) == numel(D));


end %Matlab executes everything in between

equalityTracker  %Absolute value of matches
equalityTracker == 3*numTries %Does it match what we expect?

%% New segment:
%Bsxfun applies a function to both inputs
%This is the first time we pass a function to a function 
%You do this via function handles
%They are denoted by @
bsxfun(@eq,C,D) %Testing for equality with bsxfun and function handles

%Another example to generalize. Say we want to scale it (normalize)
A = pascal(7)
%Idea: Divide each element by the max of A
%Let's do this with bsxfun. We could do without. 
B = bsxfun(@rdivide,A,max(max(A))) %Useful because fast without loop

%% Visualize what each part of the SVD does
close all
figure
%Create a random rectangle and plot it
M = randn(1,4) %Corner X, corner y, width, height
xTrace = [M(1,1), M(1,1)+M(1,3), M(1,1) + M(1,3), M(1,1), M(1,1)]; %x trace
yTrace = [M(1,2), M(1,2), M(1,2) + M(1,4), M(1,2) + M(1,4), M(1,2)]; %y trace
plot (xTrace,yTrace,'color','b','linewidth',3)
axis equal
xlim([-5 5])
ylim([-5 5])
hold on %So we can see both
plot(0,0,'.','color','k','markersize',20)
pause
sFx = 2; %scale factor x
sFy = 3; %singular value y
S = [sFx 0; 0 sFy]; %Diagonal matrix S, stretch one dimension by scale factor x, the other by y
rot = 10; %Rotation in degrees
theta = rot/180*pi; %Convert degrees to radians (which Matlab thinks in)
%Put the traces into a matrix so we can rotate it
R(1,:) = xTrace;
R(2,:) = yTrace; 
%We can now use an orthogonal matrix to rotate. Let's call it VTranspose
vTranspose = [cos(theta) sin(theta); -sin(theta) cos(theta)]; %This will help us rotate
rotatedXY = vTranspose*R ; %Matrix multiply to get the rotation
plot(rotatedXY(1,:),rotatedXY(2,:),'color',[89/255 0 179/255],'linewidth',3)
rotatedAndStretchedXY = S*rotatedXY;
pause
plot(rotatedAndStretchedXY(1,:),rotatedAndStretchedXY(2,:),'color','r','linewidth',3)
%xlim([min(min(rotatedAndStretchedXY) 3])
%ylim([-3 3])
%Fix this after class to make it flexible 
rot2 = -45;
theta2 = rot2/180*pi; %Convert degrees to radians (which Matlab thinks in)
U = [cos(theta2) sin(theta2); -sin(theta2) cos(theta2)]; %This will help us rotate
finalXY = U *rotatedAndStretchedXY; %Boom
pause
plot(finalXY(1,:),finalXY(2,:),'color','g','linewidth',3)

%Exercise for the reader: Put this into a function that takes the colors
%and the rotations (in degrees) and the singular values as input argument

%% 2 Saving and loading of variables
whos %Tells you what is in your working memory (the workspace)
save MyFirstVariables 
clear all
whos
load MyFirstVariables
whos

%% Error checking (important for HW, 2nd program control flow function we introduce)
clc
A = randn(2,3) %Make a random 2x3 matrix A
B = randn(2,3) %Same for B

try %Try this
C = A*B %This will not work out of the box because the dimensions where they touch don't match
catch %If you encounter an error, do this. Catcher in the MATLAB
    disp('You''re matrix dimensions didn''t match = we fixed that for you by transposing')  %Alert the user to this
    C = A*B'
end

%Catch is a conditional switch. The code does everything from "try" to
%"end". IF (and only if) an error is encountered, the statements between
%catch and end are executed. 

%% 3 Creating numbered vectors without for loops
%For instance, a matrix of odd numbers from 1 to 100, square roots and
%squares. Without a for loop.
%Dimensionality: 50x3
format shortg
n = 1:2:100; %This gets us all odd numbers from to 100
temp1 = sqrt(n); %2nd operation: Square root
temp2 = n.^2; %3rd operation: Square
M = [n; temp1; temp2]' %Assemble the final matrix by vertical concatenation of the rows

%% 4 Projections - the movie
%Use for loops again to make a rotation movie and understand projections
%vec = [0.75; 0]; %An arbitrary vector
vec = [0.75; -0.75]; %An arbitrary vector

%base1 = [0 1]; %x base
%base2 = [1 0]; %y base
base1 = [1 1] %New basis 1
base2 = [1 -1] %New basis 2

%We are on our road to projections. We need a lot of stuff along the way.
%We need unit vectors, length of the vectors, dot products, etc.
%Luckily, we already made all the functions last time
magVec1 = vecLen(vec); %Magnitude of vector
%magVec2 = vecLen(base1); %2nd one
%magVec3 = vecLen(base2); %3rd one

uVec = makeUnitVec(base1); %Unit vector from base 1
uVec2 = makeUnitVec(base2); %Unit vector from base 2

figure
set(gcf,'color','w') %Make it white

startPoint = 0; %We can change this later
endPoint = pi/2; %Rotate by 90 degrees, we will change this later
steps = 100; %In 100 steps
%Getting the increment theta (per iteration)
temp = linspace(startPoint,endPoint,steps); %Get step number of points linearly spaced between startPoint and endPoint
temp2 = diff(temp); 
for ii = 1:2*steps-1 %Do the rotation steps number of times

    if ii > 99
    theta = -temp2(1); %If we are tapped out, rotate back
    else
    theta = temp2(1); %Otherwise, leave as is

    end

    eero2(base1,'b',2)
    hold on
    eero2(base2,'b',2)
    eero2(vec,'r',3)
    axis equal
    axis square
    xlabel('x component')
    ylabel('y component')
    box off %Tufte would rotate himself if we left the box on
    
    dotProduct = dot(vec,base1); %Use the inbuilt function dot
    angleBetween = angBetVec(vec,base1); %Right? 
    p = magVec1 * cos(deg2rad(angleBetween)); %Length of projection
    projVec = p * uVec; %This is actual projected vector, onto base 1
    eero(projVec,'k',4)
    
    %We need to do all of this again for base 2
    %We should put this into a function, so we don't make a mistake
    dotProduct2 = dot(vec,base2); %Use the inbuilt function dot
    angleBetween2 = angBetVec(vec,base2); %Right? 
    p2 = magVec1 * cos(deg2rad(angleBetween2)); %Length of projection
    projVec2 = p2 * uVec2; %This is actual projected vector, onto base 2
    eero(projVec2,'k',4)

    xlim([-0.1 1.1])
    ylim([-1.1 1.1])
    set(gca,'xtick',[-1 -0.5 0 0.5 1])
    set(gca,'ytick',[-1 -0.5 0 0.5 1])
    hold off
    
    R = [cos(theta) -sin(theta); sin(theta) cos(theta)]; %Our rotation matrix R
    vec = R*vec;
    pause(0.05)
end %Close the loop