% linear algebra basics

switch opt
    
    case 'a'
    
    % properties of matrix multiplication
    m = 5
    p = 3;
    n = 4;
    a = rand([m p])
    b = rand([p n])
    
    % must be conformable to do matrix multiply
    a*b
    size(a*b) 
    
    % transpose of a product: reverse order and transpose each
    (a*b)' - b'*a'
    
  case 'b'
    
    % square matrices are always conformable
    n = 5;
    a = rand([n n])
    b = rand(size(a));
    
    % matrix multiplication is not cummutative
    a*b
    
    % ... is not equal to ...
    b*a
    
  case 'c'
    
    % mean and variance of a data matrix easily done
    % in terms of matrix multiplication
    
    % data matrix; each column is e.g. a time series
    m = 20;
    n = 5;
    d = rand([m n]);
    
    % row vector of means of each column of the data matrix
    d_mean = mean(d,1);
    
    % d'*d give covariance matrix
    % diagonal is sum of squares of each time series (each column of d)
    % divide by size(d,1)=N (length of time series)
    d_variance = diag(d'*d)'/size(d,1) - d_mean.^2
    
    % compare to:
    tmp = var(d,1,1)
    
  case 'd'

    % projection on to orthogonal basis set
    %
    % make a random square matrix of dimension n
    n = 5;
    a = rand([n n]);
    c  = a*a'; 
    % (didn't explain in class, but the above step makes the matrix 
    %  positive definite so that the eigenvectors will be real)
    
    % this arbitrary matrix has columns that are not orthogonal
    disp('multiply two colums of a matrix (not equal to zero)')
    c(:,1)'*c(:,2)
    
    % but its eigenvectors will be orthogonal
    % get the eigenvectors (more on this later in class)
    [e,l] = eig(c)
    
    disp(['Multiplying e(1)'' by each other e(k)']) 
    for k=1:size(e,1)
      disp([ 'k= ' int2str(k) '    ' num2str(e(:,1)'*e(:,k))])
      pause(0.5)
    end
    % the eigenvectors are a spanning or basis set so that     
    % any vector can be projected (expressed as a linear weighted
    % sum) on to the set of basis vectors
    
    % random vector, for example
    f = rand([n 1])
    
    % can therefore be represented as a weighted sum of the basis vectors
    % f = sum {k=1,n} a(k)*e(k) = e*a
    
    % a(k) is obtained from
    k=3;
    e(:,k)'*f
    
    % or the entire vector of weights from
    a = e'*f;
    
    % because e' = inv(e)
    e'-inv(e)
    
    % verify that this is correct
    disp('misfit in projection of f onto basis vector set e:')
    misfit = e*a-f    % should be all zeros
    
    disp(['|max misfit| is ' num2str(max(abs(misfit)))])

end