%% Vector basics

%% a vector is an ordered sequence of numbers
x = [2 4 6 8]

y = [2:2:8]
z = [8:-2:2]

% getting an element out of a vector
y(3)

% basic operations
x/2
2*x
whos

% the number of elements of a vector
length(y)
% note: this is the number of elements, not the "Euclidean length"

% loop over a vector
for i=1:length(y)
    2*y(i)
end

% plotting vectors
plot(y)
hold on   % don't overwrite this
plot(z, 'r')  % plot additional line in red colour

whos  % look at your workspace
clear all  % clear your workspace
whos  % look at your (now empty) workspace
close all  % close your figures

%% angle between vectors
x = [0 1];
y = [1 1];
plot([0 x(1)], [0 x(2)], 'r', [0 y(1)], [0 y(2)], 'k');
% what is the angle between these two vectors?

% cosine of angle between vectors
cosxy = x*y' / (norm(x)*norm(y));

% get the angle from cosine
acosxy = acos( cosxy )

% Why does is this not what we expected?
% tranform radians to degrees:
acosxy * (180/pi)


%% plotting functions of vectors

x = [-5:0.1:5]

% square every element of x
y = x.^2
plot(y)

% look at the x-axis of the plot - what's wrong?
% this is better
plot(x, y);
plot(x, y, 'x');

% a few more things
hold on
plot(x, (x-1).^2, 'r')
hold on
plot(x, x.^2+5, 'k')
plot(x, 2*(x-1).^2, 'g')

close all

% plot sine and cosine
sinx = sin(x);
cosx = cos(x);

plot(x, sinx, 'k', x, cosx, 'r')

% correlation
corr(sinx, cosx)
% produces lots of NaNs because dimension not appropriate, better:
corr(sinx', cosx')
% sine and cosine are orthogonal - why is this correlation not exactly 0?

%% Vector multiplication
% we can start afresh from here

% Define two parallel vectors:
x = [1 0 1 0]
y = [2 0 2 0]

% Multiply vectors
x*y  % doesn't work...

% we have to "transpose" the vector on the right
x*y'

% transposing turns row (horizontal) vectors into column (vertical) vectors and vice versa
% this is because vectors are special types of matrices - we'll come to this later
x'

% now try
x'*y

% it really matters what you are doing! a times b is not always b times a!

%% Functions for two vectors

% correlation between two vectors
% try
corr(x, y)
% and you'll get a strange result
% the input should be column vectors, so do it like this

corr(x', y')

% ttest
ttest(x', y')
% produces a suspiciously significant result...
% check how the ttest works:
help ttest

% this produces the p-value
[h, p] = ttest(x', y')

% if you want the t-value:
[h, p, ci, stats] = ttest(x', y')
stats.tstat  % "stats" is a structure, with multiple elements

% for simple measures
sum(x)
mean(x)
std(x)


%% Matrices

% Define a matrix
x = [1 1 1; 2 2 2]

% Getting elements out of matrices
x(1,2)
x(:,1)
x(2,:)

% Getting matrix dimensions
size(x) % returns number of rows and columns, respectively
[m,n] = size(x)

%% Matrices c'd
% Matrix transpose
x
y = x'

% visualise a matrix
imagesc(x)
imagesc(y)

% adding/multiplying values etc.
x + 2
2 * x
x.^2

% Define another matrix
z = [1 2 3; 4 5 6]

% Subtract z from x
x - z

% Add z to x
x + z

% !!! Elementwise multiplication (NOT "matrix multiplication": "*")
x .* z

% Try... (matrix multiplication)
x * z
% doesn't work because dimensions don't match (we'll get to this later)

% Elementwise division (NOT "matrix division/inversion": "/")
x ./ z

% Try... (matrix inversion)
x / z

%% Multiplying matrix and vector
% Define a vector
v = [1 1 1]'

% Mutliply matrix x with vector v (and obtain sum within each row of x)
x
v
x * v

% matrix and vector dimensions must agree
v * x % does not work

%% interpretation as sum of column vectors
% x * v results in the sum of column vectors in x
x * v
% produces the same result as
1*x(:,1) + 1*x(:,2) + 1*x(:,3)

%% Matrix multiplication
% Multiply two matrices
x * y

% Multiply two matrices
y * x

% !!! Note: x*y is not necessarily y*x!!!

% Define new matrix
z = [1 1; 2 2]

% Works:
y*z
z*x

% Doesn't work:
z*y
x*z

%% Make your own matrices

% Identity matrix
id = eye(3)

% Example:
x*id

id*x'

% diagonal matrix
di = eye(3)
di(1,1) = 1
di(2,2) = 2
di(3,3) = 3

% Multiplies columns by a factor:
x*di

% Multiplies rows by a factor
di*y

%% Some useful functions for matrices

% sum/average across columns of a matrix
% note: most Matlab function operate column-by-column
sum( x )
mean( x )

% sum across rows of a matrix
sum( x' )
mean( x' )

% get the diagonal of a (square) matrix
mydiag = diag(di)

% dimensions of a matrix, size()
[m,n] = size(x)

% dimension of a vector
diaglen = length(mydiag)

% create a matrix/vector with ones
ones(10)   % creates 10x10 matrix
ones(10,1) % creates 10x1 column vector

% create a matrix/vector with zeros
zeros(10)   % creates 10x10 matrix
zeros(1,10) % creates 1x10 row vector

% create matrix/vector with standard uniform random numbers in interval 0-1
rand(3)

% create matrix/vector with standard normal random numbers in interval 0-1
randn(3)

% correlation among columns of a matrix
R = rand(5,3)
% correlation among columns
corr(R)
% correlation among rows
corr(R')
% visualise a matrix
imagesc( corr(R') )
colorbar


% "reshaping" a matrix or a vector
x = 1:10
reshape(x,5,2)
reshape(x,2,5)

% Kronecker product ("concatenates" matrices/vectors)
kron(id,y)

%% EMEG example
load EMEGdata.mat
whos
data
data.evoked
data.evoked.epochs

% get only EEG data
mydat = data.evoked.epochs(307:376,:)
plot( mydat' )

% find bad channel
max(max(mydat))
max(max(mydat'))
[a,b] =max(max(mydat'))

% remove bad channel
mydat(63,:) = []
whos

plot( mydat' )

[m,n] = size( mydat );

% create "average reference operator"
avg_op = (eye(m)-ones(m)/m);

% this subtracts the mean across electrodes at each sample (column)
mydatref = (eye(69)-ones(69)/69)*mydat;

% this should look different:
plot( mydatref' )

% a different way of plotting 2D data
imagesc(mydatref)
colorbar

% plot correlation matrix among electrodes
imagesc(corr(mydatref'))
colorbar