%Script lab 2 - Matrices and matrix operations

%We already did functions, which we will use now

%% 0 Vectors: Stacks of numbers
a = [1 2 3] %Row vector
b = [1;2;3] %Column vector

%Vectors: lower case

%Indexing of vectors: Simple
a(3) %Indexing of vectors is simple, because only one dimension is not 1
b(2)

%% 1 Matrices: Stacks of vectors

%a) Matrix creation - denoted by upper case

A = [1 2 3; 4 5 6; 7 8 9] %Stacking row vectors

size(A) %Dimensionality
numel(A) %Total number of elements

%b) Accessing elements and reassigning (indexing)
%Because multiple dimensions are now in play, it is not as simple as with
%vectors

%Entire 2nd row
A(2,:) %Colon (here!) means "all"

%Last 2 elements of the 3rd row
A(3,2:3)
A(3,2:end) %There are often many ways to do the same thing

%First column
A(:,1)

%3rd element of 3d row, and reassigning it to 55
A(3,3) = 55

%Deleting the entire 3rd row
A(3,:) = []; %If you delete from a matrix, you have to delete entire rows and columns, you can't delete individual elements

%Vectorize all surviving elements of A (linearize or flatten)
B = A(:) %Without any other inputs, : flattens a matrix

%By default, Matlab operates column-wise. 
%We can make it do it row-wise, but we need to tell it to do so, by using
%the transpose
B = A'; C = B(:)

%If you want the result to be row-vector, you need to transpose the output
C = B(:)'
%Order of operations matters

%c) Matrix operations
%Seamless. This is the #1 strength of MATLAB. Everything is optimized to
%do fast matrix operations. Why? Every matrix operation is usually a
%one-liner (Matlab reads line by line), no matter how big the matrix is. 
mean(A) %Without arguments, this is column-wise
mean(A,2) %But you specify the direction in which you want to take the mean

%Matrix addition
A(3,:) = [7 8 9] %I want our square matrix back
B = ones(3) %3x3 square matrix of ones
C = A+B %As long as A and B have the same dimensionality, matrix addition works fine

%d) Holes in a matrix
%Just like nature abhors a vacuum, MATLAB abhors holes in a matrix
rtP1 = [ 1000 1500 500 1200] %Reaction time of participant 1 in 4 trials
rtP2 = [ 2000 1700 3000 700] %Reaction time of participant 2 in 4 trials
RTs = [rtP1;rtP2] %Matrix that represents all participants

TM = mean(RTs,1) %Per trial mean
PM = mean(RTs,2) %Per participant mean

rtP2 = [2000 1700 3000]
%RTs = [rtP1;rtP2] %Matrix that represents all participants

%This throws an error because MATLAB doesn't know HOW to stack the row
%vectors. We don't know where the hole is.
%To fix this, we need to represent the hole explicitly. 
rtP2 = [2000 1700 3000 nan]
RTs = [rtP1;rtP2]

%Nans allow to represent missing values. But at a terrible price.
%Whatever dimension you take the mean over, as before, the entire result
%becomes a nan. Nan makes a nan whatever it touches. 
TM = mean(RTs,1) %Per trial mean
PM = mean(RTs,2) %Per participant mean

%Nans allow you represent holes, but you still can't compute with them.
%You have to address them separately.

%How to handle nans? 
%1) Matlab in recent years recognized the pain imposed by nans, so there
%are new functions that compute as if the nans were not there
nanmean(RTs,1)
nanmean(RTs,2)
%Ignores the nans

%Find the nans and destroy them
%Question: Can we locate the nans in our matrix? 
lI = find(isnan(RTs)==1) %If you ask find to return 1 output argument, it yields the linear index of the nans, as if the matrix was flattened
[rr,cc] = find(isnan(RTs)==1) %If you ask for 2 or more output arguments, it yields the rows of columns of the matrix that are nan
%Package this as a function "nanfinder", on your own. 
%If you don't remember this, you can always use sub2ind to make a linear
%index out of subscripts (subscripts is what you use to access matrix
%elements), and ind2sub does the opposite. 

%2) More traditional option: Throw out all the missing data - the entire
%trial or the entire participant? 
%Entire trial: 
RTs(:,4) = []; %Get rid of the 4th column
%or
RTs(2,:) = []; %Get rid of the 2nd participant
%Both is not necessary

%% 2 Plotting of matrices as vectors (last time we plotted vectors)
%What if you have more than one? You can represent that as a matrix
%a) with plotvec
close all
figure
M = [1; 1]; %A matrix of column vectors
plotv(M) %plotv is an ancient function that plots each column as a vector
M = [1 -1 -1 1; 1 1 -1 -1]
plotv(M)
%plotv is so old that we can't get a handle on it. So how do we take
%control? By finding the objects that correpond to the vectors and getting
%a handle on those. Matlab creates an object for every graphical element
%you make, handle or not. So you can always take control, even if you can't
%explicitly declare a handle. 
findobj
%findobj gives you a list of all existing graphics objects
%Find all the lines:
lineHandle = findobj('type','line'); %Finds and gets a handle on all the lines
set(lineHandle,'color','k','linewidth',7) %Black and thick

%% b) As a spike plot - you could easily package this as a function 
numVec = 15;
Y2 = randn(1,numVec); %Draw a specified number of random numbers from a standard normal distribution without replacement
%This will represent the end points of our vectors
Y1 = zeros(1,numVec); %Draw a specific number of zeros, which will represent our vector starting positions
Y = [Y1; Y2];
X = 1:numVec; %Our x base
X = [X;X]; %Because we will use vertical lines to represent vectors, this will do - each element of a vertical line has the same x coordinate
figure
q = line(X,Y,'color','k','linewidth',2); 
hold on
plot(X,Y2,'.','color','m','markersize',16); %Adding a little red dot on top
title(['Spike plot of vectors of a matrix with dimensionality: ', num2str(numVec)])
xlim([0 numVec+1])
%Adding a horizon line - a horizontal line each vector starts from, at 0,
%from beginning to plot to end
line([min(xlim) max(xlim)],[0 0], 'color','r','linestyle',':') %Adding a thin red line to represent the horizon
%Importantly, this should not share properties of how we represent vectors,
%so as not to confuse ourselves (linestyle, color, thickness)

%% 3 Revisiting inner products and matrix multiplication

%Inner products
rowVector = [1 2 3 4];
columnVector = [5;6;7;8];
innerProduct = rowVector*columnVector
dot(rowVector,columnVector) %Same
sum(rowVector.*columnVector') %Same

%Matrix multiplication 
A = [1 2; 3 4];
B = [5 6; 7 8];
C = A*B
D = B*A
%Matrix multiplication does not commute. Neither do dot products, as a
%special case of matrix multiplication. Only element-wise multiplication
%commutes
E1 = A.*B
E2 = B.*A
E1 == E2

%% 4 Doing something useful with that
%Tuesday!
%Also introduce bsxfun, also reshaping and function handles,
%per request


