% Percentage signs signify comments to Matlab. Everything after is for you,
% Matlab ignores it. But it is useful for comments to yourself.

%Such as the fact that Matlab is an interpreted language and that lines are
%read in sequential order, according to the numbers on the far left. 

%% Double percentages indicate code "segments". It's also autobolded 
%It allows you to execute a whole paragraph of code in one go. We will use
%this extensively, e.g. one paragraph on vectors as dot products, one on
%vectors as projections, and so on. Suggestion: Organize *your* code in
%logical segments like this. 

%% Every script file should start with a header with the following information: 
%1) What is this (What does it do, what does it assume and what does it
%output)
%2) Who did it? 
%3) How can I reach them? 
%4) When was this done (Version history)

%% Pretend header: Lab 1 - from 0 to 0.01
%This script file illustrates MATLAB's capabilities as a graphing linear
%algebra calculator 
%You - and your email 
%09/08/2017

%% 0 Initialization step - skipping it today because we didn't introduce functions yet

%% 1 Arithmetic operations - everything works like you would expect

5 + 2 %Plus
5 - 2 %Minus
5 * 2 %Multiplication
5 / 2 %Division
5 ^ 2 %Exponentiation
5 ^ 0.5 %Square root

%% Order of operations and blank space: Also everything as expected
(5 + 4/2)^2 %Let's do some contrast
5 + 4 / 2^2 %This works

5 + 2
5 +             2
5 + - 2 %It simply ignores the blankspace
%Ctrl-Enter: Run the highlighted line
%Command-Enter: Run the highlighted segment
%Green arrow: Run whole program
%Caveats: On Macs. Check your shortcuts. You can set all of that in the
%preferences. 

%% 2 Variables - in MATLAB, we store valus (usually numbers) in these

%ans:
%Every time you do a calculation without explicit assignment to a variable,
%Matlab creates a temporary variable "ans" for you. It is updated every
%time you run something. So it makes for poor permanent storage. But it is
%very useful for iterative or recursive processes. 
5 + 5 %ans will become 10 
ans + 5 %ans will become 15
ans + ans %ans will become 30

%More permanent storage: 
A = 5
B = 2.5
%Now these variables have been created in the workspace (=RAM) and you can
%use them for other computations after.
C = 2*B
D = 5*C + A + B

%The equal sign in MATLAB means an assignment operation. It's radical. It means take whatever is on the
%right side of it, and replace the left side with it. So if you had
%important information on the left, it is lost. 

%It is *not* as the looks would suggest a test of equality. It
%suspiciously looks like it, but it is not. One of the major sources of
%logical errors if you code at 3 am, racing against a deadline. 

%Test for equality in Matlab: Double == 
A == A %True
A == B %False
A == C %True

%Returns Booleans (0 or 1)

%Now you know how to make variables. We should talk what you can name them.
%In principle, almost anything you want. But most choices are unwise. 
%Importantly: You can't have blank spaces in them, you should not overwrite
%existing functions. I recommend shortened camelCase. Descriptive, but easy
%to parse

%Bad: size = 500; %Unclear what this is the size of and overwriting the
%function "size"
%Better: 
numParticipants = 500; %Avoids underscores. More power (n = 500)
num_participants = 500; %Weak? Do whatever you want. It's kind of like tabs and spaces

%% 3 Functions and operators: If variables are your matter, functions are your forces

%Recall that we did square root with ^0.5. There is also a function for
%that: 
9^0.5 %This works
sqrt(9) %This also works

%There are over 7000 built-in functions in Matlab, there is basically a
%function for everything. 

%Everything in Matlab is case-sensitive. So perhaps name your functions and
%variables in a way that is easy to spell - it's the most common mistake.
5 + 2 %With operator
plus(5,2) %With function

%We will use this one a lot:
abs(2)
abs(-2)

%Bring up the help. It is very helpful
help abs %Brings up quick and dirty inline help
%doc abs %Opens up an encyclopedia. Very extensive

%If you don't even know what you are looking for, use the function browser.
%You can even drag and drop from that

%Functions can be concatenated, chained, nested and so on. We will do a lot
%of that
sqrt(-3) %This throws us into the complex plane, where we will spend a lot of time. 
sqrt(abs(-3)) %Let's keep it real 
%The innermost function of a nest of functions is executed first

rand(1) %Makes a random number, and a different one each time. Drawing from a uniform distribution from 0 to 1 with replacement
randn(1) %Draws from a standard normal distribution with replacement
%The argument is the dimensionality of the matrix of numbers coming out
%If you get the same numbers as each other, we're using the same seed - I
%can't get into this right now. Later

%Some things that can go wrong:
abs(-2)
abs = abs(-4) %Say I want to store this
abs(-2) %Now abs is defined as a variable, and we are trying to index the -2th element
%This throws a s yntax error - albeit a something cryptic one. But that's
%fine. Worse would be a logical error: 
abs(1) %Yields 4, not 1. Because it is now interpreted as the 1st element of the variabl abs, which is 4.
%This is a logical error. 

%How do you fix this, so you can use it again? 
clear abs
abs(-2) %Works again
%If you used to overwrite clear, you were out of luck. Since Matlab 2016b,
%they fixed that somehow. With a backdoor hack. That might not always work.
%In that case, go to the workspace and delete clear manually. 

%clear = 5 %Don't assign clear to anything
%P = 7

%% 4 Vectors! Vectors are stacks of numbers

%Row vector with 5 elements
A = [1, 2, 3, 4, 2] %1) We just introduced a 2nd set of parantheses: Square brackets. It denotes a set of numbers
%, means horizonatl concatenation. We could have omitted it. It's silent.
%It's implied. 

%Column vector with 3 elements
B = [3; 5; 4] %Semicolon *here* means vertical concatenation 

%Dimensionality: To be a vector, all but one of the dimensions has to be 1
%How to assess dimensionality
size(A)
size(B)

%Never call anything size

%Dimensionality: 1st: Number of rows, 2nd: Number of columns, 3rd: Sheets,
%4th: Workbooks, 5th to nth: Unnamed, but you can have as many as you want.
%A vector has 2 and one of them is 1. It's just a stack of numbers. 

%Another use for the semicolon: Let's create two vectors with 1000 numbers
A = [1:1000] %Colon means: "from, to"
B = [1:1000]; %Semicolon at the end of a line means "echo off". It is still created, but not shown. Much faster, much cleaner

A = 27; %You can easily reassign variables. Everything on the left is overwritten/discarded/thrown away

%For the sake of completeness: The 3rd use of the semicolon is to
%concatenate commands in the same line: 
C = [1; 2; 3]; D = sum(C); %All the semicolons used in this line have a different meaning


%% 5 Vector operations

%Converting a row vector to a column vector and vice versa
B = 1:10; size(B)
C = B' %Transposing with an operator
size(C)
D = transpose(B) %Transposing with a function 
size(D)

%Most functions work *natively* with vectors in Matlab
mean(B) %Mean function with argument B
sum(B)

%% Vector length
ranVec = [1 5]; %Pretty random 
magVec = sqrt((ranVec(1).^2 + ranVec(2).^2))
%Regular parantheses serve 3 functions: Indexing, like here and giving
%function arguments, like above and it strucures order of operations
%Charlie: Here carrot will work because there is only one element. But if
%you want element-wise carrot (you usually do), get into the habit of
%putting a dot in front of it.

%Trust, but verify - always check whether your code does what is supposed
%to.
knownVec = [3 4]; %Pretty random 
magVec2 = sqrt((knownVec(1).^2 + knownVec(2).^2))

%If you reuse code logic more than once, it is probably to write a
%function, *NOT* copy and paste. We'll do that next time. We will also
%generalize this to vectors with more numbers. 

%% Creating a unit vector
newVec = [7 4]; %Random new vector
magVec3 = sqrt((newVec(1).^2 + newVec(2).^2))
uniVec = newVec./magVec3 %Dot divide is element-wise division
magVec4 = sqrt((uniVec(1).^2 + uniVec(2).^2)) %Checking that this worked

%format long displays more digits, internally it works with whatever
%precision your variable has, e.g. single or double or integer. 
%format short makes the output display of the numbers shorter

%% Scalar multiplication of a vector
%Let's say we have our unit vector, but we now want to make it longer
scal = 5; 
scaledVec = scal .* uniVec
magVec4 = sqrt((scaledVec(1).^2 + scaledVec(2).^2)) %Checking that this worked

%% Creating a null vector
nullify = 0; 
nullVec = nullify .* uniVec
magVec5 = sqrt((nullVec(1).^2 + nullVec(2).^2)) %Checking that this worked

%% Vector addition
vec1 = [1 4]; 
vec2 = [5 7]; 
vec3 = vec1 + vec2 %Works as expected

%% 6 Vector plotting 
%3 ways of doing this, they all have subtle differences: Line, plot,
%quiver. %Next time, we will make this a function. We will call it eero. 
%For obvious reasons

%% a) As lines
figure %Open a figure
line([0 vec1(1)],[0 vec1(2)],'color',[89/255 0 179/255]) %Vector as a purple line. The last argument is the color as RGB values
pause %Pause pauses the program flow until you press a key. In this case this is useful to build up the figure. So you see them sequentially
line([0 vec2(1)],[0 vec2(2)],'color','b'); %2nd vector in blue
pause %Pause again
line([0 vec3(1)],[0 vec3(2)],'color','r'); %2nd vector in blue

%This works. Totally legit, but some people like arrows or dots at the end
%of the vector. So let's do that. 

%% b) By using "plot"
figure %Open a new figure
plot([0 vec1(1)],[0 vec1(2)],'color','g') %Vector as a green plot line
hold on
plot(vec1(1),vec1(2),'.','markersize',20,'color','k')
pause(1)
plot([0 vec2(1)],[0 vec2(2)],'color',[255/255 140/255 0]) %Vector as an orange plot line
plot(vec2(1),vec2(2),'.','markersize',20,'color','k')
pause(1)
plot([0 vec3(1)],[0 vec3(2)],'color','m') %Vector as a magenta plot line
plot(vec3(1),vec3(2),'.','markersize',20,'color','k')
%In contrast to line, plot by default overwrites what is already in the
%figure. It only doesn't do that if you literally hold on.
%Hold on is needed for plot (but not for line), otherwise you will only
%see the last thing you plotted. Pause is only needed to build up the
%figure

%% c) Several people urged to use arrows. For good reason. We'll use quiver for that
figure
ceruleanBlue = [42, 82, 190]./255; %Cerulean blue
hotPink = [255,105,180]./255; %Hot pink
quiver(0,0,vec1(1),vec1(2),0,'color',ceruleanBlue) %Vector as arrow in cerulean blue
hold on %quiver works like plot, not like line, so we need to hold on
pause
quiver(0,0,vec2(1),vec2(2),0,'color','g') %Vector as arrow in green
pause
quiver(0,0,vec3(1),vec3(2),0,'color',hotPink) %Vector as arrow in green

%% 7 Vectors as operators
%The dot product is the sum of the element-wise vector multiplication (if
%you don't do anything else)

%Use case 1: Summer
%Let's say we record from all 302 neurons of C Elegans and pretend that
%they fired action potentials. We want to know the expected spike count of
%all action potentials at a given time in the C elegans brain if it had a
%brain. Let's say there are 7 cell classes
neuralFiringRates = [0.5 1 2 4 8 0.25 3]; %Firing rates of cell classes
numberOfNeurons = [45; 15; 35; 20; 50; 36; 101]; %Number of neurons per class

%Check that we count all neurons
totalNumberOfNeurons = sum(numberOfNeurons)
expectedSpikeCount = neuralFiringRates * numberOfNeurons %With an operator
%with a function: 
expectedSpikeCount = dot(neuralFiringRates,numberOfNeurons) %With a function

%How about the average firing rate? 
averageSpikeCount = expectedSpikeCount./totalNumberOfNeurons %This works, using a previous result
averageSpikeCount = dot(neuralFiringRates,(numberOfNeurons./totalNumberOfNeurons)) %Same output. Averager

%% 8 Dot products and the angle between vectors
%dot product = m1 * m2 * cos(angle1 - angle2) %Result of trig identities
%Solve this for angle between
vec1 = [0 0.5]; %v
vec2 = [1 1]; %Will be direction of u
magVec1 = sqrt((vec1(1).^2 + vec1(2).^2)) %Magnitude of vector 1
magVec2 = sqrt((vec2(1).^2 + vec2(2).^2)) %Magnitude of vector 2
dotProduct = dot(vec1,vec2); %Using a function 
angleBetween = rad2deg(acos(dotProduct./(magVec1.*magVec2)))

%% 9 Dot product as projection
%Gives the length of the piece of v in the direction of the unit vector
% (v * u) * u
% Length of projection times direction of unit vector
% v * u = p = magnitude of v times cos of the angle between
% p = magnitude of v * cos angle? 
% Then: p * u

uVec = vec2./magVec2; %Making it an unit vector

p = magVec1 * cos(deg2rad(angleBetween))
projVec = p * uVec; %That's the actual projected vector

figure
quiver(0,0,vec1(1),vec1(2),0,'color','r')
hold on
quiver(0,0,uVec(1),uVec(2),0,'color','b')
quiver(0,0,projVec(1),projVec(2),0,'color','k','linewidth',2)
axis equal
axis square

%% 10 Orthogonal vectors
% (- sin theta) or (sin theta)
% (cos theta)      (-cos theta)
orthVec = [sin(deg2rad(angleBetween)) -1 * cos(deg2rad(angleBetween))] %orthogonal vector
%Arithmetic check:
dot(uVec,orthVec)
quiver(0,0,orthVec(1),orthVec(2),0,'color','k','linewidth',1)

%% Owed
%Distance from line
%Spike plot
%Functions (your own)
%All of this with matrices, maybe SVD