%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Matlab Demos
% (adapted from http://www.stanford.edu/class/cs223b/matlabIntro.html)
%
% Stefan Roth <roth (AT) cs DOT brown DOT edu>, 09/10/2003
% Alex Balan  <alb  (AT) cs DOT brown DOT edu>, 09/14/2004
% Leonid Sigal<ls   (AT) cs DOT brown DOT edu>, 09/14/2005
% Greg Shakhnarovich <gregory (AT) cs DOT the same place> 09/04/2006
% Deqing Sun <dqsun (AT) cs DOT brown DOT edu>  09/14/2009
% Jason Pacheco <pachecoj (AT) cs DOT brown DOT edu>  02/02/2011
% Soumya Ghosh <sghosh (AT) cs DOT brown DOT edu>  02/02/2011

% Last modified: 09/14/2009
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% The symbol "%" is used to indicate a comment (for the remainder of
% the line).

% Within Matlab interpreter window, you can evaluate a portion of code by
% selecting it and pressing F9.

% help <topic> is useful; even if you don't immediately see the information
% you are looking for, the "SEE ALSO" pointers often lead to it.
help help
help general

% this is

A = [1, 2; % note: once the matrix input has started, Matlab allows you
           % to break the line, and correctly inteprets the rest of the input
     3, 4];

% A semicolon at the end of a statement means that Matlab will not
% display the result of the evaluated statement. If the ";" is omitted
% then Matlab will display the result.  This is also useful for
% printing the value of variables, e.g.

A

% Matlab's command line is a little like a standard shell:
% - Use the up arrow to recall commands without retyping them (and
%   down arrow to go forward in the command history).
% - In linux, C-a moves to beginning of line (C-e for end), C-f moves
%   forward a character and C-b moves back (equivalent to the left and
%   right arrow keys), C-d deletes a character, C-k deletes the rest of
%   the line to the right of the cursor, C-p goes back through the
%   command history and C-n goes forward (equivalent to up and down
%   arrows), Tab tries to complete a command.
% You can execute a linux command in matlab by prepending "!" to it
!ls -l
!touch testing
!rm testing

% Simple debugging:
% If the command "dbstop if error" is issued before running a script
% or a function that causes a run-time error, the execution will stop
% at the point where the error occurred. Very useful for tracking down
% errors.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The basic types in Matlab are scalars (usually double-precision
% floating point), vectors, and matrices:

N = 5                        % A scalar

A = [1 2; 3 4]               % Creates a 2x2 matrix
B = [1 2; 3 4]               % The simplest way to create a matrix is
                             % to list its entries in square brackets.
                             % The ";" symbol separates rows;
                             % the (optional) "," separates columns.

v = [1 0 0]                  % A row vector
v = [1; 2; 3]                % A column vector
v = v'                       % Transpose a vector (row to column or
                             %   column to row)
v = 1:.5:3                   % A vector filled in a specified range:
v = [-4:4]                   %   [start:stepsize:end]
v = [-4:0.1:4]               %   default stepsize is 1
v = linspace(-4,4,100)       %   (brackets are optional)

v = []                       % Empty vector

isempty(v)
isempty([1 2])

b = zeros(3, 2)             % Matrix of zeros
c = ones(2, 3)              % Matrix of ones
d = eye(4)                  % Identity matrix
d2 = diag([1:4]/4)          % Diagonal matrix
e = rand(64, 128);          % Random matrix (uniform) between [0,1]
f = randn(64, 128);         % Random matrix (normal) N(0,1)
g = zeros(128);             % Matrix of zeros (128x128)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Visualization: plotting  etc.

x = [0 2 1 3 4];             % Basic plotting
plot(x);                     % Plot x versus its index values starting from 1


figure;                      % Open new figure.

% by default, Matlab will plot using current figure.
% figures are identified by handles -- integers starting from 1.
% New figure is created with the lowest available handle.


x = pi*[-1:1/24:1]           % Vector of arguments for functions to be ploted
plot(x, sin(x));
xlabel('radians');           % Assign label for x-axis
ylabel('sin value');         % Assign label for y-axis
title('dummy title');        % Assign plot title

figure(1);                   % Use the figure with handle 1 for
                             % subsequent plotting commands
subplot(1, 2, 1);            % Multiple functions in separate graphs
plot(x, sin(x));             %   (see "help subplot")
axis square;                 % Make visible area square
subplot(1, 2, 2);
plot(x, 2*cos(x));
axis square;

figure(2);                   % Use the figure with handle 2
plot(x, sin(x));
hold on;                     % Multiple functions in single graph
plot(x, 2*cos(x), '--');     % '--' chooses different line pattern
legend('sin', 'cos');        % Assigns names to each plot
hold off;                    % Stop putting multiple figures in current
                             %   graph
pause
close all


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% (C) Indexing vectors and matrices.
% Warning: Indices always start at 1 and *NOT* at 0!

v = [1 2 3]
v(3)                         % Access a vector element

v(5)   % this causes an error (out of range)
% however, it is OK to add an element:
v(5)=5    % note how the 4th element got the default value of zero


m = [1 2 3 4; 5 7 8 8; 9 10 11 12; 13 14 15 16]
m(3, 2)                      % Access a matrix element
                             %       matrix(ROW #, COLUMN #)

% a range of indices can be specified. Note the row transposition effect
% due to the negative step -1 in the range expression:
m(3:-1:1, [1 3])

% a matrix can be visualized using image or imagesc commands
figure;
imagesc(m)
% colors correspond to values according to colormap (help colormap)
% the color "legend" can be displayed explicitly:
colorbar

% more matix stuff:

m = rand(128,256);
size(m)                      % Returns the size of a matrix
size(m, 1)                   % Number of rows
size(m, 2)                   % Number of columns
length(m)                    % max(size(m))
m1 = zeros(size(m));         % Create a new matrix with the size of m,
                             % fill it with zeros

% ----  Simple operations on vectors and matrices

% ---- Reshaping and assembling matrices:

a = [1 2; 3 4; 5 6]          % A 3x2 matrix
b = a(:)                     % Make 6x1 column vector by stacking
                             %   up columns of a
sum(a(:))                    % Useful:  sum of all elements

c = reshape(a, [2,3])    % Make 2 x 3 matrix out of a 3 x 2 matrix
                             %   elements are arranged column-major
			     % note: reshape will throw error if sizes
                             % don't add up
% Matrix concatenation
d = [c -c]

e = repmat(a, 2, 2)         % Create a 4x4 matrix

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% (B) Vector Operations
% Built-in Matlab functions that operate on vectors

a = [1 4 6 3]                % A row vector
sum(a)                       % Sum of vector elements
mean(a)                      % Mean of vector elements
var(a)                       % Variance of elements
std(a)                       % Standard deviation
max(a)                       % Maximum
min(a)                       % Minimum

% If a matrix is given, then these functions will operate on each column
%   of the matrix and return a row vector as result
a = [1 2 3; 4 5 6]           % A matrix
mean(a)                      % Mean of each column
max(a)                       % Max of each column
max(max(a))                  % Obtaining the max of a matrix
mean(a, 2)                   % Mean of each row (second argument specifies
                             %   dimension along which operation is taken)
sum(a,2)                     % sum of each row
sum(a)                       % the default dimension is always 1 (i.e.,
                             %   sum of each column in this case).
			     % Note: the result of such an operation is
                             % normally a vector!


% ' performs matrix transposition

[1 2 3] * [4 5 6]'           % 1x3 row vector times a 3x1 column vector
                             %   results in a scalar.  Known as dot product
                             %   or inner product.  Note the absence of "."

[1 2 3]' * [4 5 6]           % 3x1 column vector times a 1x3 row vector
                             %   results in a 3x3 matrix.  Known as outer
                             %   product.  Note the absence of "."

% Other elementary functions:
% >> abs(x)       % absolute value of x
% >> exp(x)       % e to the x-th power
% >> fix(x)       % rounds x to integer towards 0
% >> log10(x)     % common logarithm of x to the base 10
% >> rem(x,y)     % remainder of x/y
% >> sqrt(x)      % square root of x
% >> sin(x)       % sine of x; x in radians
% >> acoth(x)     % inversion hyperbolic cotangent of x

% Other element-wise arithmetic operations include e.g. :
%   floor, ceil, ...

help elfun   % get a list of all available elementary functions

% ---- element-wise versus matrix operations:
a = [1 2 3; 4 5 6]
b = [10 20 30; 40 50 60]

% "x .<op> b" means perform operation <op> on each pair of matching (by
% indices) elements of a and b. This will cause error if the sizes are
% not equal.

% compare this command
a.*b
% to this (illegal) command
a*b

b./5

% if one of the arguments is a scalar, it is used for each element:
a.^2
2.^a


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% "Cell Arrays"
% Cell arrays in Matlab are multidimensional arrays whose elements can
% be anything. Cell arrays are created by enclosing a
% miscellaneous collection of things in curly braces, fg:

S = 'Hello World!'
R = rand(2,3,2)

C = {R -1 S}

% To retrieve the contents of one of the cells, use subscripts in
% curly braces.
C{3}

% Like numeric arrays also cell arrays can be mulitdimensional
CM{2,2} = S;
CM{1,1} = 'hallo';
CM

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% (4) Control statements & vectorization

% Syntax of control flow statements:
%
% for VARIABLE = EXPR
%     STATEMENT
%      ...
%     STATEMENT
% end
%
%   EXPR is a row vector here, e.g. 1:10 or -1:0.5:1 or [1 4 7]
%     - it should never be a column vector
%
%
% while EXPRESSION
%     STATEMENTS
% end
%
% if EXPRESSION
%     STATEMENTS
% elseif EXPRESSION
%     STATEMENTS
% else
%     STATEMENTS
% end
%
%   (elseif and else clauses are optional, the "end" is required)
%
%   EXPRESSIONs are usually made of relational clauses, e.g. a < b
%   The operators are <, >, <=, >=, ==, ~=  (almost like in C(++))
%
%   Element-wise Logical operators :
%     and: &
%     or: |
%     not: ~
%     exclusive or: xor
%   For short-circuit: &&, ||

% Warning:
%   Loops run very slowly in Matlab, because of interpretation overhead.
%   This has gotten somewhat better in recent versions, but you should
%   nevertheless try to avoid them by "vectorizing" the computation,
%   i.e. by rewriting the code in form of matrix operations.  This is
%   illustrated in some examples below.

R = [4 5] > [3 7]
any(R)                  % True if any element of a vector is nonzero
all(R)                  % True if all elements of a vector are nonzero


% Examples:
for i=1:2:7                  % Loop from 1 to 7 in steps of 2
  i                          % Print i
end

for i=[5 13 -1]              % Loop over given vector
  if (i > 10)                % Sample if statement
    disp('Larger than 10')   % Print given string
  elseif i < 0               % Parentheses are optional
    disp('Negative value')
  else
    disp('Something else')
  end
end


% Here is another example: given an mxn matrix A and a 1xn
% vector v, we want to subtract v from every row of A.

m = 50000; n = 10; A = ones(m, n); v = 2 * rand(1, n);
%
% Implementation using loops:
tic                          % start stopwatch
for i=1:m
  A(i,:) = A(i,:) - v;
end
toc                          % stop stopwatch

% We can compute the same thing using only matrix operations
tic                          % start stopwatch
A = ones(m, n) - repmat(v, m, 1);   % This version of the code runs faster
toc                          % stop stopwatch

%%%% New matlab version does code optimiation about for loop, yet avoding
%%%%    for loops can still save some time, esp. if you have to use an
%%%%    older version


% We can vectorize the computation even when loops contain
%   conditional statements.
%
% Example: given an mxn matrix A, create a matrix B of the same size
%   containing all zeros, and then copy into B the elements of A that
%   are greater than zero.

% Implementation using loops:
tic                          % start stopwatch
B = zeros(m,n);
for i=1:m
  for j=1:n
    if A(i,j)>0
      B(i,j) = A(i,j);
    end
  end
end
toc                          % stop stopwatch

% All this can be computed w/o any loop!
tic                          % start stopwatch
B = zeros(m,n);
ind = find(A > 0);           % Find indices of positive elements of A
                             %   (see "help find" for more info)
B(ind) = A(ind);             % Copies into B only the elements of A
                             %   that are > 0
toc                          % stop stopwatch

% The indeces "find" returns can be converted back to row-column index
% using ind2sub.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% fprintf : Can be used to print output to the screen or to a file.

% Print to the screen
x=pi;
fprintf('\nThis is a test, x=%g \n',x)
% This command prints the number x in
% either exponential or fixed point formant
% (whichever is shorter) and then starts a new
% line.

disp(x)
disp(['X = ' num2str(x)])       % Display by string concatenation

% Print to a file called 'output.txt':
fid=fopen('output.txt') %fid=fopen('output.txt'. 'w')
fprintf(fid,'hello %g', x)

% Get input from the keyboard, e.g.,
a = input('Left Endpoint a = ?')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%(6) Creating scripts or functions using m-files:
%
% Matlab scripts are files with ".m" extension containing Matlab
% commands.  Variables in a script file are global and will change the
% value of variables of the same name in the environment of the current
% Matlab session.  A script with name "script1.m" can be invoked by
% typing "script1" in the command window.

% Functions are also m-files. The first line in a function file must be
% of this form:
% function [outarg_1, ..., outarg_m] = myfunction(inarg_1, ..., inarg_n)
%
% The function name should be the same as that of the file
% (i.e. function "myfunction" should be saved in file "myfunction.m").
% Have a look at myfunction.m and myotherfunction.m for examples.
%
% Functions are executed using local workspaces: there is no risk of
% conflicts with the variables in the main workspace. At the end of a
% function execution only the output arguments will be visible in the
% main workspace.

a = [1 2 3 4];               % Global variable a
b = myfunction(a .^ 2)       % Call myfunction which has local
                             %   variable a
a                            % Global variable a is unchanged

% NOTE:
% When writing a long Matlab statement that becomes too long for a
% single line use "..." at the end of the line to continue on the next
% line.  E.g.

[c, d] = ...
  myotherfunction(a, b)      % Call myotherfunction with two return
                             % values

%    %%%In "myfunction.m" %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    function y = myfunction(x)
%    % Function of one argument with one return value
%
%    a = [1 1 1 1];              % Have a global variable of the same name
%    y = a + x;
%    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%    %%%In "myotherfunction.m" %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    function [y, z] = myotherfunction(a, b)
%    % Function of two arguments with two return values
%
%    y = a + b;
%    z = a - b;
%    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%(5) Saving your work

who                          % List variables in workspace
whos                         % List variables w/ info about size, type, etc.

save myfile                  % Saves all workspace variables into
                             %   file myfile.mat
save myfile a b              % Saves only variables a and b

clear a b                    % Removes variables a and b from the
                             %   workspace
clear                        % Clears the entire workspace

load myfile a                % Loads variable a

load myfile                  % Loads variable(s) from myfile.mat

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%(6) Numerical Tricks

% Here is an example of UNDERFLOW. Multiplying a large number of even
% moderately small numbers can cause an underflow. 
p1 = repmat(0.1,1000,1)
log(prod(p1))
%Instead convert to the log domain
sum(log(p1))

% suppose we want to solve Ax = b for x
A = [ [ 1, 2 ]; [ 3, 5 ] ];
b = [ 1; 5 ];

% we COULD simply invert A
x = inv( A ) * b

% If A is large then this is SLOW and numerically unstable
% Alternatively, we could use the left-divide operator \
x = A\b

% Working in Log-Domain is more numerically stable.  For example, suppose
% we want to multiply a bunch of numbers and take the log
v = repmat( 1/eps, [1,20] )
log( prod( v(:) ) )

% This is an example of a numerical OVERFLOW.  Remember we can transform the
% log of a product to the sum of logs
sum( log( v(:) ) )    % no overflow




