%% 

clear
clc

%% Define a matrix 

% by entering numbers directly
A = [1 2 3; 4 5 6; 7 8 9]

A = [1 2 3;
     4 5 6;
     7 8 9]

% or by combining row vectors
a = 1:5
b = 2:6
A2 = [a; b]

% or by combining column vectors
a = a'
b = b'
c = (3:7)'

B = [a b c]

% or by reshaping a vector
C = reshape(1:12, 3, 4)

% or by replicating matrices
C3 = [C C C]
C13 = repmat(C, 1, 3)

C31 = repmat(C, 3, 1)

% or randomly
R = rand(3,5)

% zero matrix
O = zeros(3,4)

% identity matrix
I = eye(4)

% create a diagonal matrix from a vector
D = diag(1:2:10)

whos

save mat_file A B C* R
% load mat_file

%% Indexing and assigning values in a matrix
clc

A = reshape(1:9, 3, 3)

A(1:2,:) % first two rows
A(1,:) = 0 % set first row = 0

A(1:2,1:2) %top 2x2 block
A(:,end) % last column

A(:) % convert to vector in column fashion

A(1:4:end) % extract all diagonal entries
A(1:4:end) = -1 % modify all diagonal entries to zero
% for nxn matrix, digonal is A(1:n+1:end)

diag(A) % also extract the digonal of A, 
        % but cannot modify values
%diag(A) = -2 % error!!!

A(A<0) = 0 % set all negative values to 0

%% matrix properties

A = reshape(1:15,3,[])

[m,n] = size(A) % vector
size(A,1)
size(A,2)

numel(A)

A = A(:,1:end-1) % drop last column
A(:,end) = [] % alternative way

% A is now a square matrix
issymmetric(A) % logical value

%% single matrix operations

A + 3 % add 3 to all entries

A.^2 % entrywise square

A^2 % square of A (matrix power)
A*A

transpose(A)
A'

sum(A)
sum(A,1) % sum along 1st dimension (column)
ones(1, size(A,1))*A % not efficient!

sum(A,2) % sum along 2nd dimension (row)
A*ones(size(A,2),1) % not efficient!

sum(A,'all') % add all entries 
sum(sum(A)) 
sum(A(:)) % A(:) is a vector!
ones(1,size(A,1))* A * ones(size(A,2),1) % not efficient!

mean(A,1)

min(A,[],1)
max(A,[],1)

trace(A)
sum(diag(A))

det(A)

D = eig(A) % eigenvalue decompostion
[V,D] = eig(A) % both eigenvalues and eigenvectors

% The following functions should NOT
% be used on large matrices

rank(A)

inv(A) % inverse

%% matrix-vector and marix-matrix operations

B = magic(3)

A + B

A - B

A.*B % Hadamard product

A*B % matrix product

a = [1 0 -1]
a*A

A * a'

%% techniques for multiplying large matrices

n = 5000;

A = rand(n,n); 
B = rand(n,n);
x = rand(n,1);

% wall clock time
tic; y1 = (A*B)*x; toc  % slow way

tic; y2 = A*(B*x); toc % much faster!

norm(y1-y2, Inf) % largest difference in absolute value

% another way to measure time
t0 = cputime;
(A*B)*x;
t1 = cputime-t0

t0 = cputime;
A*(B*x); % much faster!
t2 = cputime-t0

% matrix mutliplication involving a diagonal matrix
tic; diag(x)*A; toc  % slow method

tic; 
repmat(x,1,n).*A; % much faster
% x.*A            % equivalent
toc
%% sparse matrices

A = zeros(5,5);
A(1,1) = 2; A(1,2) = -1; A(2,1) = 1;

A

A_s = sparse(A) 
% A = full(A_s) 

% sparse matrices are very efficient, and
% all the previous functions apply to sparse matrices

A_s.^2

A_s^2

sum(A_s,2)

B = reshape(1:10, 5, 2)

A_s*B
A*B

%% 3D arrays

A = zeros(3,3,2)
A(:,:,1) = 1
A(:,:,2) = rand(3)

size(A)
ndims(A)
numel(A)

B = reshape(1:18,[3 3 2])
B(1,3,2)

sum(B,3) % add the slices together

% use a cell array to store matrices
C = {B(:,:,1), B(:,:,2)}
