%% Lecture 3, ACM 11
%% Review of Matrices
% We'll need to discuss eigenvalues, so first a quick review of matrices.
% This is on the blackboard, so for a reference, see MathWorld or
% Wikipedia.
% Topics: diagonalizable, normal, jordan canonical form
%   orthogonal/unitary, symmetric/hermitian matrices.
% Are eigendecompositions unique (or, when are they unique, and in 
% what sense?)?
%   Schur form, SVD

%% Structures
% Structures are a datatype, and are very useful for *organizing* several
% variables, especially in conjunction with arrays.  We'll start with the
% easiest use of a structure, where we have a single variable (no arrays).  
% The distinguishing feature of a structure is that it can have "fields",
% and as many fields as you wish.  To make a variable a structure, you can use
% the "struct" command, but this is usually not necessary. If a variable
% doesn't yet exist (or is set to something like an empty matrix), then
% just referring to a field will tell Matlab to make the structure for you.
%  Access fields is done with the "." mark.  For example,
clc; clear;
myStructure.month = 'October';  % this will create the variable "myStructure"
myStructure.day = 13;           % we create a new field
% The fields can represent have different datatypes
% If we display myStructure, this is what we see:
myStructure
% Or we can see just a particular field:
fprintf('The current month is %s\n', myStructure.month )
% There are some useful commands associated with structures.  To see them,
% see the "see also" section after running the "help struct" command.  In
% particular, we can use "fieldnames" to see a list of fields in the
% struct.  For much more information, in the Help browser (the graphical
% one), in the "Contents" tab, select "Matlab >> Programming Fundaments >>
% Built-In Classes (datatypes) >> Structures".
%% Dynamic field names of structures
% Suppose we have many fields, and want to access them all in an automated
% fashion.  In the example below, we'll make a "structure of structures"
% (an alternative would be to make a "cell" of structures, or to use 
% structures in combination with array).
clc; clear;
% Make a smaller level structure
experiment.date = 'Tuesday, October 13';
experiment.hypothesis = 'Undergraduates need 2 hours of sleep per night';
experiment.conclusion = true;

% Make a higher level structure
sociologyClass.experiment1 = experiment;

% Now, suppose we wish to add many experiments
for n= 2:8
    fieldName = sprintf('experiment%d',n);
    sociologyClass.(fieldName) = [];       % we don't yet have data
end

% We used the format   sociologyClass.(fieldName), and NOT the format
%                      sociologyClass.fieldName
% The parenthesis tell Matlab that we want to access the value of the
% variable inside "fieldName", and that the field is not actually called
% "fieldName".
sociologyClass
%% Arrays of structures, and structures of arrays
% The usage of "structures of structures" is very versatile, and the
% low-level structures need not look identical (i.e. they can have different
% fields).  If all the low-level structures are identical, then we can use
% an array.  There are two options: an array of structs, and a struct of
% arrays.
clc; clear;
% === A struct of arrays ===
structArray.dayOfMonth(1) = 15;
structArray.dayOfMonth(2) = 7;
structArray.dayOfMonth(3:5) = [6, 5, 20];
structArray.digitsOfPi(1:7) = [3 1 4 1 5 9 6];
structArray.name(1,:) = 'Bob';
structArray.name(2,:) = 'Sam';
% structArray.name(2,:) = 'Nancy';  % would not work, because
                                    % Nancy is a different length than Bob
structArray

% === An array of structs ===
arrayStruct(1).dayOfMonth = 15;
arrayStruct(2).dayOfMonth = 7;
% arrayStruct(3:5).dayOfMonth = [6, 5, 20];   this does not work!
arrayStruct(5).dayOfMonth = 20;

arrayStruct(1).name = 'Bob';
arrayStruct(2).name = 'Nancy';
arrayStruct
%%
% Note that we actually have 5 spots allocated for the "name" field in
% arrayStruct, since we referenced the 5th element in referring to
% "dayOfMonth".
arrayStruct.name

%%
% So, which is better, an array of structs, or a struct of arrays?  Both
% can be useful.  For small amounts of data, use either.  For large
% amounts of data, a struct of arrays is preferred, because this makes one
% structure and many arrays, and an array is a more "efficient" datatype
% than structure (however, this overhead can be irrelevant if the value 
% of the field takes up a lot of memory).  An array of structures makes 
% many structures, and is more inefficient.  For example, 
clc; clear;
for j=1:40000
    structArray.field(j) = pi;
    arrayStruct(j).field = pi;
end
%%
% To see how big these variables are, we again use stuctures!
w = whos('structArray','arrayStruct') 
% the "whos" command returns an array of structs, with fields like "name",
% "size", and "bytes".
for j = 1:length(w)
    fprintf('Variable %s takes up %.3f MB of RAM\n', w(j).name, w(j).bytes/(1024^2) );
end

%% Cells
% Cells are another Matlab datatype.  A cell is like an array, except that
% every element of an array must be the same datatype (e.g. a character, or
% a double).  Cells are not restricted like this.  In syntax, the
% difference is that we use { and }; these curly braces are used for both
% creating and access elements of a cell (whereas for arrays, we have
% separate notation for both these).  For example, 
clc; clear
myCell = {'First entry is a string', -.364, [1 2; 3 4], ones(5) }
% Like structures, you can explicitly create a cell (use the command
% "cell"), but this is rarely necessary.

% This is a useful way to visualize the contents of a cell"
clf; cellplot(myCell); title('from cellplot');
%%
% There are numerous functions for converting between datatpes (e.g.
% "cell2mat"), especially between strings and cells.  Cells are a useful
% way to store strings (which are 1-D arrays of characters), because
% storing a 2-D array of characters requires you to pad short strings with
% zeros to make everything line up.   Note that both cells and structs can be
% nested and intermixed.
%%
% For much more information on cells, in the Help browser (the graphical
% one), in the "Contents" tab, select "Matlab >> Programming Fundaments >>
% Built-In Classes (datatypes) >> Cell Arrays".
%%
% A word of caution: it is possible to use ( ) to access cell indices, but
% this doesn't do what you might think: instead of getting the element
% inside the cell, you actually get a single cell containg the element.
myCell = {'This is a string','Ceci n''est pas une corde' };

% To demonstrate, let's call the "deblankl" funuction, which is just a
% random Matlab function that operates on strings (it takes away whitespace
% and makes everything lowercase):
 disp( deblankl( myCell{1} ) );       % this is fine
%disp( deblankl( myCell(1) ) );       % this will give an error message

disp(myCell{1})     % with the curly brackets, we get a string
disp(myCell(1))     % with parenthesis, we get a new cell containing a string

%% Functions
%{
    Functions are different than scripts.  They take a specified input
(usually), and give an ouput (usually), or perhaps print to the screen.  In 
general, most functions don't have "side-effects".  This means that, if 
you disregard their output, then calling the function won't affect the rest
of the code.  But, this is not always true, and not always desired; we'll
see this effect when we use "global" variables.  As an example, the function
"tic" has global effect.
    A function is defined in a separate file (unlike the anonymous
functions that we saw last lecture), and this file should have the same
name as the function (although not always required).  The syntaxt is:

function output = functionName( input )
[code...]

or, for several inputs/outputs:

function [output1, output2] = functionName( inputA, inputB)
[code...]

In the rest of the function body, all the variables that weren't passed
as inputs are "local".  This means they are completely different variables
than variables of the same name elsewhere.  

There is no need to "end" the function, but you can if you wish; to do this,
just type "end" at the very end.  You do not need to have inputs or outputs,
so a line like:

function functionName

is valid.

Functions can have "sub-functions".  These are available ONLY to the main
function, and not to any other part of MATLAB.   To make a sub-function,
below the rest of the function code just use the same function syntax.  You
can make as many subfunctions as you want.  Again, there is no need to "end"
a subfunction, because Matlab can determine this automatically; but, if you do,
use the "end" keyword.  If you do this once in the function, you must do it
for all subfunctions.

One major advantage to Matlab is its flexibility in dealing with different
inputs and outputs.  You do not always need to pass all the arguments into
a Matlab function, and you do not always need to read out all possible output 
values.  Inside the function, you can see how many parameters were passed
in with the "nargin" command (it returns an integer), and the number of
output variables is returned by the "nargout" command.  These are useful
for setting default values of parameters.

Another way to handle default arguments is by using structures.  The
disadvantage of "nargin" is that if we define a function like

function out = functionName( inA, inB, inC, inD)

and want to pass in values for "inA" and "inD", while using defaults
for "inB" and "inC", then we still need to give values to "inB" and "inC" 
(otherwise "inD" will be interpreted as "inB").  Usually, we assign a 
value like [ ] (the empty matrix), and check for this in the body of
the code using an "isempty" test.  But this is inelegant.  Structures can
be used to avoid this.  For example,

function output = functionName( essentialInput, options )
if ~isfield(options,'tolerance'), options.tolerance = 1e-6; end
...

If a certain option isn't set, then we gracefully set it to a default
value.

%}
%%
% Examples: to see subfunctions (and use of nargin and nargout), type "edit
% fzero".  To see many more examples, just edit your favorite Matlab
% function.  Many of the numerical functions show their source code
% (whereas functions that do low-level operations, such as loading a file,
% will not generally show the source code).
%% Documentation
% In addition to adding comments to your code, there is a special section
% that is very useful.  After the "function output = functionName(input)"
% line, the following commented lines are the "help" for that function.
% If you type "help functionName", it will show every comment that is
% written below the first line, until the first non-comment line (which
% includes any blank line).  This is very useful!  Use this, especially if
% someone else will need to use your code.
%% Path
% When you type the name of a function, Matlab needs to know where to
% search for it.  It can't search the entire hard drive, so it restricts
% itself to directories in the "path".  To see which directories are in your
% path, just type "path", although you're likely to get more output than
% you want.  There are many functions relating to the path.  The most
% useful are "addpath", which adds a certain path to the directory,
% and "pathtool", which is a graphical tool that lets you manipulate 
% the path.  The current directory is always in the path.  For any directory
% in the path, if it has a a subdirectory called "private", then the files 
% in this subdirectory are also available, but only to files in the parent
% directory.  This allows you to have many functions with the same name,
% but in different directories, without getting mixed up.
%
% To see where a function resides, use the "which" command.
% 
% For example, if I want to use a Preconditioned Conjugate Gradient solver
% (this is a solver for huge linear systems), I can use Matlab's "pcg"
% function.  This resides in C:\Program
% Files\MATLAB\R2008a\toolbox\matlab\sparfun
which pcg
if ispc
    winopen('C:\Program Files\MATLAB\R2008a\toolbox\matlab\sparfun')
end
% We see that there is a subdirectory called "private", which contains
% functions like "interapp.m".  This file, "interapp", is only available to
% functions in the "sparfun" directory (such as pcg).  I cannot see it!
which interapp
% "which" is also useful if you have several files with the same name.
% Matlab will call the first one it finds -- in other words, the *order* of
% directores in the path is important!  Use the "-all" flag for "which" to
% see all functions with the same name (e.g. "which -all times" will show
% all the different functions called "times"; we say that this function is
% "overloaded", meaning that there a different function is called depending
% on the type of input).
%% Global variables
% Global variables are variables that can be accessed anywhere.  However,
% before a function can access a global variable, it must declare it.
% And to change the value of a global variable in the main workspace,
% you also must declare the variable.  To declare a variable, use the
% form:

global myGlobalVariable1 myGlobalVariable2
% (do NOT separate these with commas!)

% for examples, type "doc global".  When you "clear" variables,
% global variables are unaffected unless you type "clear global"
% or "clear all".  (Note: a useful command, similar to "clear", is "close",
% and is used to close figure windows; "close all" closes all figure
% windows).
%%
% On a related note, there are several "workspaces".  This is mainly for
% advanced use of Matlab, so isn't very important.  To use this feature,
% the functions "evalin" and "assignin" are useful.  If, for example, you
% want to evaluate the "whos" command to get a list of variables in the
% main workspace, but you call "whos" from a function, then you can't see
% the main workspace variables unless you tell Matlab to use the correct
% workspace (via the "evalin" command).
%% switch statements
% Used to simplify "if" and "elseif" statements.  See the help
% documentation:
doc switch
%% funfun
% Can a function take another function as an input?  Yes!  We did it last
% week using anonymous functions.  If we want to use a normal function, we
% can pass the function as an argument using a "function handle".  Recall
% our "ezplot" example from lecture 1.  This can accept a function handle,
% and use the function to plot.  To get the "handle" to a funuction, just
% type "@" before the function name.  Simple!
clf; ezplot( @cos )
% To see the most common Matlab functions that operate on other functions,
% type "help funfun".  These include functions that will find the minimum
% of other functions, as well as ODE solvers.
%
% To evaluate a function handle, usually you can treat it like a function.
% For example,
f = @cos;
f(pi)
% Sometimes you want to use "feval" (this is more robust; it accepts not
% only function handles, but strings with the name of a fuction).  The
% syntax is:
feval(f,pi)
%% Minimization/Maximization, Optimization
% We are not going to cover the optimization toolbox.  The conventional
% wisdom is that Matlab is very good at numerical linear algebra and
% related operations, but is relatively poor at optimization.  But there
% are many excellent 3rd party Matlab applications, some of them free, that
% do optimization very well.  A particularly simple package that I
% recommend is "cvx" by Michael Grant and Stephen Boyd.
%
% For minimizing nonlinear functions, we can classify the algorithms
% according to how much information we have about the function.  Do we know
% the derivative?  Do we know the second derivative?  
%
% The basic functions for this in Matlab are "fminbnd", "fminsearch" and
% "fzero".  If you have the optimization toolbox, then "fsolve" is another
% useful function.  These do not require the derivative.
%
% Most of these functions use structures to pass in additional arguments,
% as mentioned above.

% Example:
% What's the minimum of cosine, in the range 0 to pi?
clear;clc;
opts.TolX = 1e-2;
xmin = fminbnd( @cos, 0, pi, opts );
fprintf('Using a tolerance of %.2e, error is %.3e\n', opts.TolX, xmin - pi);
opts.TolX = 1e-9;
xmin = fminbnd( @cos, 0, pi, opts );
fprintf('Using a tolerance of %.2e, error is %.3e\n', opts.TolX, xmin - pi);
%%
% If you have the derivative information, and it's not computational
% expensive to compute the derivative at a point, then a good starting
% point is Newton's Method.  This basically forms a linear model of your
% function, and finds the exact minimum of the linear model every step.
%
% Be careful: Newton's Method sometimes refers to minimization, and
% sometimes to root-finding.  The are related, since minimization of f is
% the same as finding a root of the derivative of f (if f has continuous
% derivatives).  For root-finding, Newton's method needs a derivative; for
% minimization, it actually needs a second-derivative.
%% ODEs
%{
    ODE stands for Ordinary Differential Equation.  These differential
equations can involve several variables (often, we think of a single variable,
but allow it to be a vector).  They also involve derivatives of this vector 
variable (say, x) with respect to a new variable (say, t).  Example:

x'' + 3*x' + t*x = -t + 9

where x' means d/dt (x).

    Types of ODEs: autonomous, linear, homogenous, constant-coefficient.
Note that all ODEs with higher derivatives can be turned into systems
of more variables, but only with first derivatives (for example, if we have
x'' = 3, we turn this into:
y' = 3;
x' = y;

    Basic ODEs can be solved with the matrix exponential.  To solve
the matrix exponential, you may have learned to diagonalze the matrix.  
This is works well if, for example, if the matrix is symmetric, and we want
to find the value at many time points.  But if we only need a single time
point, there are much better ways to do it in practice for an arbitrary matrix.
For a general matrix, use Matlab's "expm" function.

    For tough ODEs that are not so simple, we can approximate the solution
with numerical methods.  The standard Matlab solver is "ode45".  Type 
"doc ode45" for lots of good info.  To learn more about the algorithms,
take ACM106.
%}

%% a stiff ODE
% This is perhaps the simplest possible example of a stiff ODE.  In this
% example, the two variables decouple, so there is no reason to solve for
% them together (we could add a week coupling, although the solution
% wouldn't be as simple).  But, suppose we solve for both variables at
% the same time.  Because the two variables (y_1 and y_2) 
% live on such different time scales, the standard solver ODE45 is
% confused, and solves for both variables on the small time scale.  This
% requires very many time steps (since we have to run it for a long time --
% due to the variable on the long time scale).
%
% $$\frac{d}{dt} y_1 = -\alpha y_1 $$
%
% $$\frac{d}{dt} y_2 = -\beta y_2$$
%
% $$ \alpha << \beta $$
alpha = 1;
beta = 100000;
f = @(t,y) [-alpha*y(1); -beta*y(2)];
tic; [T,Y] = ode45(f,[0,1],[1;1] ); toc
tic; [Ts,Ys] = ode15s(f,[0,1],[1;1] ); toc
fprintf('Number of time samples used in ODE45:  %d\n', length(T) );
fprintf('Number of time samples used in ODE15s: %d\n', length(Ts) );
fprintf('l_infinity error via ODE45:  %.2e %.2e\n',...
    norm( Y(:,1) - exp(-alpha*T),'inf'), norm( Y(:,2) - exp(-beta*T),'inf') )
fprintf('l_infinity error via ODE15s: %.2e %.2e\n',...
    norm( Ys(:,1) - exp(-alpha*Ts),'inf'), norm( Ys(:,2) - exp(-beta*Ts),'inf') )
% We find that the "ode15s" solver is much better at handling the different
% scales than the "ode45" solver is.  "doc ode45" explains when to use the
% different solvers.

%% Plotting, part II
% A common scenario for plotting is wanting to combine two variables on one
% plot, but have different scales for both of them.  Matlab provides the
% "plotyy" command for this, and it will also give us a chance to see
% "handle graphics" for the first time.  Handles are a tool that Matlab
% uses to identify objects.  We've already seen handles when we used
% "fopen" and when we used function handles. 
clear;clc; clf
% Example of two variables with different scales
t = 1949:2008;
stockMarket = 500*(randn(60,1)' + .05*(t-1949));
population = exp( .1*(t-1949) );
% plot(t,stockMarket)
% plot(t,population)
plotyy(t,stockMarket,t,population)
ylabel('Dow Jones Industrial Average');
xlabel('year');
%%
% Now, let's be a bit more advanced.  We'll put the stock market data
% on a normal plot, but plot the world population on a log plot, since we
% know the growth is exponential.
[ax,h1,h2] = plotyy(t,stockMarket,t,population,@plot,@semilogy);
ylabel('Dow Jones Industrial Average');
xlabel('year');
% How can we set the ylabel for the population data?  Use handles!
ylabel( ax(2), 'World Population')
% What if we want to adjust the line thickness for one of the graphs?  (in
% the help browser, search for "lineseries" to see possible options that
% can be used in this way):
set(h2,'linewidth',6)
set(h1,'marker','+')
legend([h1,h2],'Stock Market','Population','location','best')
title('A fancier graph!')
% Warning: do NOT try to edit the yticks and ylabels on a "plotyy" graph
% because it will foul everything up (the "plotyy" function goes to great
% lengths to make the left and right tick marks align).
%%
% If handle graphics scare you, then you can still get nice plots.  In the
% figure window, click the "View" menu and then "Property Editor".  Now,
% click in various items in the figure to see menus relating to their
% properties.  This is very powerful, but the disadvantage is that it is
% not automated, so if you close the plot, you must start over again.
%
% There are many other graphic tools.  For example, the linkdata button
% (also available on the command line as "linkdata on" and "linkdata off")
% will refresh the plot with new values whenever the variables that created
% it are updated (by default, this does not happen!  plots are "frozen"
% with the current value of a variable).

%% FFT (Fast Fourier Transform)
% You should refresh yourself on the FFT -- see a good 
% textbook or the internet.  It is often
% convenient to think of the FFT as a matrix; when the FFT operates on a
% vector of length N, we can instead think of a N x N matrix F that
% multiplies the vector.  Of course, we don't usually do it this way
% because a matrix-vector multiplication takes O(N^2) operations, while the
% FFT exploits redundancies and takes O(N*logN) operations.  Hence the
% name "Fast".
clear; clc;
N = 4;
F = fft( eye(N))
% Observations: F is a complex matrix.  And, a lot of the entries look
% similar.  Also, the matrix is (almost) unitary:
disp('F*F'' is:')
F*F'
disp('F''*F is:')
F'*F
disp('Is F symmetric? F - F.'' is:')
F - F.'
% It's not exactly unitary in this form because we don't get the identity,
% we get a constant (8) times the identity.  But this is of little
% importance, and in fact the IFFT (inverse FFT) will account for this
% constant.  Because F is complex, be very careful to distinguish adjoint
% (') from transpose(.')
%%
disp('F*I is:')
I = ifft(eye(N));
F*I
%%
% Now, let's take the FFT of something interesting.
N = 100;
t = linspace(0,2.7, N);
w = 3;              % linear frequency
y = sin(w*2*pi*t);
f = fft(y);
%%
% What if we plot f?
clf; plot(f)
% Huh?  Confusing.  This is because Matlab is plotting the Real part vs.
% the Imaginary part, which is not what we want. 
%%
% Instead, try this:
figure(1); plot( abs(f),'o:' ); title('fft')
% We're taking the FFT of a pure sinusoid.  Shouldn't we expect to have
% exactly one large component of the FFT, corresponding to the frequency of
% the sinusoid?  Yes and no.  First, notice the TWO peaks.  This is because
% our frequency w shows up, but so does -w.  For a real input signal, the
% FFT gives redundant information, and the first half of the signal is the
% complex conjugate of the second half, if we flip around the midpoint.
%%
% Matlab has a useful tool called "fftshift" that will do this flipping.
% To undo the effects of fftshift, use "ifftshift" (fftshift and ifftshift
% are either the same, or off by one, depending on if N is even or odd).
fftshift( 1:8 )         % to give an idea of what fftshift does
figure(2);
plot( [-N/2:0, 1:N/2-1], abs(fftshift(f)), 'ro:' )
title('using FFTshift')
%% 
% A useful way to plot:
semilogy( abs(f(1:N/2+1)), ':o' ); title('Periodogram')
%%
% Or have Matlab do it for you:
periodogram(y)  % requires signal-processing toolbox
% If you tell it the sampling rate, it will use the correct units.
%% Making a plot with error bars
% This is easy: use "errorbar" (see "help specgraph")
clear;clf;clc
x = 1:10;
y = 3*x.^2 + 20*randn(1,10);
err = 30*randn(1,10);
errorbar(x,y,err,':')
