%% Lecture 4, ACM 11
%% Making a plot with error bars
% This is easy: use "errorbar" (see "help specgraph")
clear;clf;clc
x = 1:10;
randn('state',5);
y = 3*x.^2 + 20*randn(1,10);
err = 30*randn(1,10);
errorbar(x,y,err,':');  % these are just as an example; they have no real meaning
text(1,250,'Data generated via 3*x^2 + noise');
%% Basic regression
% Using the data from the error-bar plot, let's fit it to a polynomial
p = polyfit(x,y,2);     % we want a degree-2 polynomial

estimate = polyval(p, 1:.1:10 );
hold on
plot(1:.1:10,estimate,'r')
legend('data','fit')
text(4,10, ['via polyfit: ' poly2str(p,'x')])
title('Error bar plot')

% You can also do this graphically; in the figure window, go to "Tools" and
% pick "Basic Fitting".
%% Better regression by doing it ourselves
% First, we recreate the results of polyfit "by hand"
% For a general quadratic fit, we do the following:
x = x(:); y = y(:);   % make sure x and y are column vectors
A = [x.^0, x, x.^2 ] 
coefficients = A\y    % the backslash solver minimizes the l_2 norm of the residual
% This is the same as what polyfit gave us!

%%
% Now, suppose we know that the polynomial that best explains our data has no
% constant term.  We can use this info to improve our estimate
% Now, we get rid of the 1's column:
coefficients = A(:,2:3)\y
coefficients = [flipud(coefficients);0];
text(4,-10,['forcing a zero intercept: ' poly2str(coefficients,'x')])

% And if we know it's exactly a quadratic, we would do this:
coefficients = A(:,3:3)\y
coefficients = [flipud(coefficients);0;0];
text(4,-30,['forcing a quadratic: ' poly2str(coefficients,'x')])

%% 3D plotting
% Suppose we want a 3D plot of sin(x)*cos^2*(y).  The simplest way to set
% this up is using "meshgrid".
clc;clear; N = 5;
f = @(x,y) sin(x).*(cos(y).^2);
x = linspace(0, 2*pi, N );
y = linspace(0, 2*pi, N );
[xGrid,yGrid] = meshgrid(x,y)
%%
% Here are four types of plots that are useful for 2D functions
N = 100;
x = linspace(0, 2*pi, N );
y = linspace(0, 2*pi, N );
[xGrid,yGrid] = meshgrid(x,y);
F = f(xGrid,yGrid);
% For 3D plotting tools, type "help graph3d"
figure(1); clf; subplot(2,2,1)
plot3(xGrid,yGrid,F,'o'); title('plot3')

subplot(2,2,2)
mesh(xGrid,yGrid,F); title('mesh') % better!

subplot(2,2,3)
surf(xGrid,yGrid,F); title('surf')
colorbar('location','westoutside')  % this adds a legend describing the
                                    % significance of each color
subplot(2,2,4)
contour(xGrid,yGrid,F,20);  title('contour');
% "contour" is listed in "help specgraph" along with many other useful
% functions

% In the figure window, click the "Rotate 3D" icon (let the mouse hover
% over the icons to see their names), and you can rotate these objects.
%%
% Here are some cool 3D plots that give an idea of what Matlab is capable
% of.
clf; xpklein
%%
teapotdemo
%%
% transpdemo  % this crashed my computer!
%%
% volvec  % this behaves funny too, but gives an idea of how to display 3D stream lines.  type "edit volvec" for details.
%% vector field plots
% If you have the values of the vectors, then something like "quiver" will
% do what you want.  Take our function "f" from earlier, and find the
% gradient field.
dfdx = @(x,y) cos(x).*(cos(y).^2);
dfdy = @(x,y) -2*cos(x).*cos(y).*sin(y);
% let's make a smaller grid so that it's not too crowded
close;clf;[xGrid,yGrid] = meshgrid(x(1:5:end),y(1:5:end));
quiver(xGrid,yGrid,dfdx(xGrid,yGrid), dfdy(xGrid,yGrid) )
title('Gradient field, using "quiver" command')
xlim([0,2*pi]); ylim([0,2*pi]);

%% displaying images
% Matlab can display most common image file formats, including jpeg.  To
% see all the availabe functions, type "help specgraph" and see the "Images
% display and file I/O" section.  The basic command to load an image is
% "imread()".  This will open the file and convert it to either a 2D array,
% or a 3D array if it has color.  The 3D array will basically have 3 copies
% of the 2D array, with each copy representing a different color (these do 
% not have to be colors -- they could be chrominance, luminance, etc.;  often,
% these are referred to as "channels").  The data
% in the array will probably be an unsigned integer, and not a floating
% point number; this is done to keep the RAM requirements as small as
% possible.  It's always possible to convert back to to a floating point
% number using "double()", but if you just want to display the image, this
% is unnecessary, because the display functions expect unsigned integer
% datatypes.

% This file should be on everyone's computer
% (type "which street1.jpg" to see where it's located)
clear; X = imread('street1','jpeg');
% X = imread('street1.jpg');    % this command work work also
whos X
% We see that X is a 3D array, with type "unsigned integer, 8 bits"
%%
% To actually view the image, we can use "image" or "imagesc".  "imagesc"
% is identical to "image", except that the values will be scaled to take up
% the full range of colors available.  If you need to convert something
% back to an integer form, use "uint8()".  To see all possible conversions,
% type "help datatypes" and "help isa"
imagesc(X); title(' file ''street1.jpg'', found in the ''demos'' directory');
%%
% We can view only one channel of the image, if we wish
imagesc( X(:,:,1) ); title('displaying only one channel');
%% Importing data the very easy way
% Use Matlab's graphical importer:
uiimport
%%
% Also, here's a useful function to see what variables are inside a .mat
% file (in this example, there's a file called "myFile.mat")
whos -file myFile
% and a nice way to see (and modify) the values if a variable is with the
% "openvar" command.
s = randn(9,1);
openvar s       % the form   openvar('s')   also works
%% Sparse matrices
% In a few fields, the term "sparse" has become a buzzword.  What does it
% mean?  There's no quantitative definition, but the qualitative definition
% is easy: a vector or matrix is "sparse" when it has many zero entries.
% In Matlab, there are two ways to store a matrix.  The first way, which
% we've been using so far, is to store all entries.  If we know the
% dimensions of the matrix, and have a convention of how the entries are
% stored, then there is no need to store the indices.  The second method,
% which Matlab calls the "sparse" datatype, is to only store the nonzero
% entries.  But, with this scheme, we need to store the indices also, so
% there's a penalty in memory.  It is very important to realize that
% Matlab's definition of "sparse" only refers to the way a matrix is
% stored.  The same matrix can be sparse, or non-sparse (aka "full").
% Here's an example:
clc; clear;
A = randn(5)        % Method 1 of storing a matrix, aka "full" storage
B = sparse(A)       % Method 2 of storing a matrix, aka "sparse" storage
% The "sparse" function converts a matrix to the sparse storage format.  We
% can go the other direction using the "full" command, e.g.
% A = full(B);
%%
% If we look at the memory of the two variables, we see that the sparse
% datatype takes more memory, because we also have to store the
% coordinates.  
whos
%%
% Now, suppose we have a matrix that really is sparse.  THEN the sparse
% matrix format has the advantage:
clear; A = zeros(100);
A(45,60) = 1;  % method 1
B = sparse(A)  % method 2
whos
%%
% We can create a sparse matrix on its own, i.e. we don't have to convert
% from a full matrix first.  The simplest method is the following:
clear;
B = sparse(100,100); % makes an empty 100 x 100 sparse matrix
% Note: B = sparse(100) is NOT what you want.  This would make a 1 x 1 sparse
% matrix, with the (1,1) entry being 100.

B(45,60) = 4
%%
% Just like with the default matrix storage, we can "expand" a sparse
% matrix after we define it, and this also penalizes us in terms of memory
% (since we have to re-allocate space in RAM every time we extend the
% matrix).  We can allocate memory and define a sparse matrix all at once,
% again using the "sparse" command.  Here's an example:
clear;
A = rand(45,50);   % uniform random numbers in [0,1]
[i,j,s] = find( A .* (A < .01)  );
B = sparse(i,j,s,45,50)
%%
% Some useful functions for sparse matrices (type "help sparfun" for a full
% list).

% Visualize the sparsity pattern. 
spy(B); title('Sparsity pattern of our matrix B');
% How many non-zero entries:
fprintf('matrix B has %d nonzero entries\n', nnz(B) );
% Is the matrix stored in the sparse format?  (note: this does NOT imply
% that the matrix really is sparse, only that it is stored in the sparse
% format).
fprintf('matrix A is sparse?  %d\n', issparse(A) );
fprintf('matrix B is sparse?  %d\n', issparse(B) );

%% Eigenvalues
% Matlab has a good eigenvalue solver called "eig".  This is different than
% "eigs" (which we will cover in a minute; the "s" stands for "sparse", and
% not to pluralize "eig").  "eig" does different things, depending on the
% number of outputs.  If we only have one output, it just returns the
% eigenvalues in a vector.  If we have two outputs, it gives the
% eigenvalues still (this time, as the diagonal of a matrix), as well as
% the eigenvectors.
clear; A = randn(4);
disp('The eigenvalues of A are:')
eig(A)
% (recall that eigenvalues of a real matrix are allowed to be complex; if
% the matrix is symmetric/Hermitian, then we're guaranteed to get real
% eignevalues).
disp('The eigenvalues of A are on the diagonal of matrix D');
disp('  and the eigenvectors are the columns of matrix V');
[V,D] = eig(A)

% If you don't need the eigenvectors, don't ask for them!  It will slow
% down the computation.
%%
% The command "eigs" finds only the few largest (or smallest) eigenvalues
% of a matrix.  It uses an entirely different kind of algorithm; this
% algorithm is very well-suited to sparse matrices.  Here, we'll
% demonstrate with a very large sparse matrix.
B = sprandsym(1000,.01);    % this is a 1000 x 1000 sparse symmetric matrix
                            % 1% of its entries are randomly selected to
                            % be nonzero; on those entries, the value is
                            % chosen from a uniform distribution.
spy(B); title('Sparsity pattern of the random symmetric matrix');
% Now, ask for the top and bottom 8 eigenvalues:
% (by "top", we mean those with the greatest magnitued; remember, the
% singular values may be complex)
opts = []; opts.disp = 0;  % tell the algorithm to be quiet!
opts.maxiter = 1500;       % tell the algorithm to be patient
disp('largest eigenvalues:')
D = eigs(B,8,'LM',opts)
disp('smallest eigenvalues:')
D = eigs(B,8,'SM',opts)

% Note: you can use this command to also get the eigenvectors, not just the
% eigenvalues
%%
% To see the biggest/smallest eigenvalues, in terms of the number and not
% the absolute value (this only works of eigenvalues are real, i.e. the
% matrix is symmetric), we use "LA" and "SA" instead of "LM" and "SM"
disp('largest eigenvalues:')
D = eigs(B,8,'LA',opts)
disp('smallest eigenvalues:')
D = eigs(B,8,'SA',opts)
%% Singular Value Decomposition (SVD)
% The SVD is vaguely like a a generalization of the eigenvalue
% decomposition; it gives slightly different information, and it is also
% more general (i.e. it applies to non-square matrices).  We factor a
% matrix A as the product of three matrices: A = U * S * V'.  S need not be
% square, but it is zero except on the diagonal.  There are important variations,
% called the "thin" or "compact" svd.  The diagonal entries of S are called
% the singular values, and they are always real and positive.  By convention, they
% are always listed in decreasing order.  In Matlab, the function to
% generate the SVD is called "svd()", and there is also a sparse version
% ("svds()").  Like "eig", we can ask for the singular values only, or the
% singular values and the singular vectors (the vectors refer to U and V).
clear; A = randn(5);
disp('The singular values of A are:')
s = svd(A)
disp('The singular values and singular vectors of A are:');
[U,S,V] = svd(A)
fprintf('Norm of A - U*S*V'' is %.3e\n', norm( A - U*S*V') )
%%
% The SVD is a robust method to determine the rank of a matrix.  The rank
% is simply the number of nonzero singular values.  In practice, there is a
% little difficulty because we may find tiny singular valies (1e-15) that
% are not quite zero; do these count or not?  In general, Matlab will apply
% some thresholding (probably after normalizing by the largest singular
% value).
clear; T = randn(5,2);  A = T*T';  % A is a 5 x 5 matrix, with rank 2
disp('The singular values of A are:')
s = svd(A)
disp('The rank of A is:')
r = rank(A)
%%
% The sparse version, "svds", works very similarly to "eigs".  It is not
% one of Matlab's finest algorithms, from a numerical analysis
% point-of-view.  It simply calls "eigs" on a modified version of the
% matrix.
%%
% If you are interested in how the SVD or eigenvalue decomposition is
% performed, for the non-sparse case, Cleve Moler has a fantastic visual
% demonstration (for ACM 106, this is vey useful):
% http://www.mathworks.com/moler/ncm/eigsvdgui.m
%%
% The SVD is a key ingredient in many techniques, such as Principal
% Component Analysis (PCA) and Latent Semantic Analysis (LSA).  Cleve
% Moler's textbook includes a very nice routine that demonstrates the use
% of the SVD in approximating matrices (in this case, the matrices are
% pictures).  This could lead to a good project...  
% The following file can be found at:
% http://www.mathworks.com/moler/ncm/imagesvd.m
imagesvd
%% Numerical Integration, aka quadrature
% Numerical integration is the process of approximating the value of an
% integral.  We'll only cover definite integrals with finite endpoints, but
% it can be adapted for infinite endpoints.  The general technique, called
% quadrature (not to be confused with other meanings of "quadrature"), is
% to sample the function f at a few points, interpolate these points
% piece-wise with a polynomial of degree n, integrate this polynomial
% exactly, and return this as our approximation.  The interpolation and
% exact integration can be combined into one step, so no actual
% interpolation really happens.  It is also not necessary to interpolate
% with a plain polynomial, but we will not cover this.
%
% The simplest case is using polynomials of degree n = 0, i.e. use a
% piece-wise constant approximation.  This is basically a Riemann sum.  We
% can use linear polynomials to get the midpoint rule.
%
% Matlab can do all this for us, and it will chose gridpoints adaptively.  The
% workhorse routine is "quad()", but there are other variants (type "help
% funfun" and see the quadrature section).  Simpler methods are "trapz" and
% using "sum".  Here's an example:

f = @(x) exp(-x);   % this is the integrand
a = 0;              % The integral is from a to b
b = 3;
integralExact = 1 - exp(-3);
fprintf('Exact integral is %.4f\n', integralExact );
% Use a very basic method, with N points
for N = [50, 100, 200, 1000, 10000]
    RiemannSum = sum( f( linspace(a,b,N) ) ) *(b-a) / N;
    fprintf('Riemann sum, with N = %5d, gives %.4f\n',N,RiemannSum);
end
% Use another very basic method, with N points
for N = [50, 100, 200, 1000, 10000]
    trapRule = trapz( f( linspace(a,b,N) ) ) *(b-a) / N;
    fprintf('Trapezoidal Rule, with N = %5d, gives %.4f\n',N,trapRule);
end
% and use a smarter method:
[adaptiveSimpson,N] = quad(f,a,b);
fprintf('Using matlab''s "quad", with N = %d, gives %.4f\n',N,adaptiveSimpson);
%% dialog boxes, etc.
% Matlab can popup dialog boxes (see "dialog" or "msgbox").  Here's an
% example of a waitbar -- it graphically displays a sliding counter.
h = waitbar(0,'Computing...');
for i = 1:300
    pause(.01);         % wait .01 seconds
    % (usually, a lengthy computation goes here)
    waitbar(i/300,h);
end
close(h);
%%
% For other length computations that update a graph at every step, use the
% command "drawnow" afer a "plot" command to force the graph to
% refresh, otherwise Matlab might wait to draw the graph until all the
% computation is done.  For the most efficient way to update a graph,
% change the XData using handle graphics, or link the variables (see last
% lecture for how to do this).
%% "while" statement; "break"
% The "while" statement is very similar to a "for" statement, and usually
% you can interchange the two (with appropriate modifications), but
% sometimes one or the other makes more sense.  The idea is that the code
% is executed UNTIL some condition stops.  The format is:
%
% while [test]
%   [body of code]
% end
%
% and it will execute until the test if false.  The simplest example is
% like a "for" loop:
clc;
j = 1;
while j < 10
    j = j + 1;
end
fprintf('j is now %d\n', j)
%%
% The example above doesn't take full advantage of a "while" statement.
% In general, the body of the code would something very complicated.
% Here's an example
seed = 0;
mn = 0;
while (mn < .6) && (seed < 10000)
    seed = seed + 1;
    rand('state',seed);
    mn = mean( rand(15,1) );
end
fprintf('Took %d tries to find a random vector with mean >= .6\n',seed);
%%
% Here's how we can do the same using a "for" loop, with a "break"
% statement.  Whenever a "break" statement is called, the most local "for"
% loop will quite.
for seed = 1:10000
    rand('state',seed);
    mn = mean( rand(15,1) );
    if mn >= .6
        break;
    end
	seed = seed + 1;
end
fprintf('Took %d tries to find a random vector with mean >= .6\n',seed);

%% Dealing with errors
% There are two things we'd like to do with errors: (1) when we write our
% own code, we might want to "throw" an error, and (2) when Matlab
% generates an error, we might want to do something special.  First, let's
% cover (1).  When we throw an error, execution will stop.
% The simplest way to throw an error is just using "error('message')"

% error('This is an error message')
%%
% We can also give an ID to the message (this is useful if some program is
% going to do something about the error):

% error('lectureNotes:ErrMsg1','Another error message')
%%
% For objective (2), here's how to see some information about an error.  We
% use the function "lasterror" as follows:
errInfo = lasterror
% The "stack" field is itself a structure,and contains information about
% where the error occurred (e.g. which file, which function, and which
% line).  This doesn't work when the code was executed in "cell mode"!
%%
% Also as a part of objectve (2), we can handle Matlab's errors using
% "try" and "catch".  Here's how:
try
    size( thisVariableDoesNotExist );   % because the variable doesn't exist,
                                        % Matlab will generate an error.
catch
    % If an error is thrown, then the code here will be executed:
    errInfo = lasterror;
    fprintf('Message identifier is: %s\n',errInfo.identifier);
    fprintf('Message string is: %s\n',errInfo.message);
    % and if we want to le the error continue, we can "rethrow" it:
    rethrow( lasterror );
end