% This script shows examples from 
% the topics covered in the Matlab lessons
% Ay/By 199 - Methods for Computational Sciences
% Apr. 2009, Ciro Donalek (donalek@astro.caltech.edu)

% Run the commands one by one in the Matlab Command Window

;

%%%%%%%%%%%%%%%%%%
% Matrices
%%%%%%%%%%%%%%%%%%

% Durer's Matrix (Magic Matrix)
% Why is it Magic?

% enter Durer's matrix "by hand"
% put a semicolon at the end to not display the results
% (good habit, expecially when working with big matrices)
A = [16 3 2 13; 5 10 11 8; 9 6 7 12; 4 15 14 1]

% sum, diag and fliplr are matrix and array operations
sum(A) % column sum

% to sum over the rows we need to do the transpose
sum(A') % row sum

sum(diag(A)) % diagonal sum

% fliplr can be used to extract the antidiagonal
help fliplr
sum(diag(fliplr(A))) % antidiagonal sum

% all these sum are the same, that's why it is
% called magic matrix

% array operations applied to matrices
% extract the maximum value from the matrix A
max1=max(A)
max2 = max(max(A)) % max(max1)

% column and subscriptions
a=[2:2:10]; % create an array
A(1:3,3) % access part of a matrix
A([2 4],:) % access 2nd and 4th row, all the columns
A(:,2:end) % all rows, column 2 to the last

% exceeding the indexes, this give an error
t = A(4,5)

% but this add a column!
X = A;
X(4,5) = 17

% Element by element operation
eyeD=eye(4);
prodMat=A*eyeD; % usual matrix product
prodElem=A.*eyeD; % element by element product
prodMat
prodElem

% building matrices by blocks
MB1 = zeros (2,5) % build a 2x5 matrix with all 0
MB2 = [magic(3), ones(3,2); MB1] % 5x5 matrix

% load from an external file
%transients=load('artifacts.dat'); % double precision array: see WORKSPACE


%%%%%%%%%
% Types 
%%%%%%%%%

% FORMAT does not affect how MATLAB computations are done
% it is just for display
help format

a=[0.5 0 -3];
format rat % approx with integers ratio
a
format % return to the default display

% to see the list of keywords in Matlab
iskeyword
clear all

% Structures
my.age=35;
my.country='Italy';
my.numbers=[7 9 11 13 15];

my % double click also in the workspace

% access structure elements
my.numbers(1)
my.country(1:2)

clear my % delete my from the workspace/memory

% cell array
my={35, 'Italy', [7 9 11 13 15]}
my{1}
my{2}
my{2}(1)

ischar(my{2}) % example of logical subscripting

clear my

%%%%%%%%%%%%%%%%%%%%%%%
% Structures of Arrays
%%%%%%%%%%%%%%%%%%%%%%%

a.x = 1:10;  
a.y = sin(a.x); % sin of a.x in radians

a.x(1)

%%%%%%%%%%%%%%%%%%%%%%%%
% Arrays of Structures  
%%%%%%%%%%%%%%%%%%%%%%%%

people(1).name ='Ciro';  
people(1).age = 35;
people(1).country='Italy';
people(2).name = 'Matusalem';
people(2).alt_name='Methuseleh'; % value missing for Ciro 
people(2).age = 969;
people(2).country='Unknown';
people

% people(2) ?
% people(1) ?

% isempty: true for empty array 
isempty(people(1).alt_name)
isempty(people(2).alt_name)

% what if I need the field names?
% fieldnames: returns a cell array of strings containing 
% the structure field names associated with the structure

people_fields=fieldnames(people)

% How can I get all and only the names?
% solution 1
people(1).name
people(2).name

% solution 2
people(:).name

% solution 3
a=people(:).name

% solution 4 -> OK
a={people(:).name}

% what if we want the ages as numbers?
ages={people(:).age} % still a cell array
mean(ages) % ?
ages2=cell2mat({people(:).age}) % array of doubles!
mean(ages2) % ?

% other examples with a cell array
soccer{1}='Italy';
soccer{2}='Championship';
soccer{3}=[1934 1938 1982 2006];
cellplot(soccer)

% nesting cell arrays
soccer{4}={'Italy','France','Spain','Germany'};
soccer{4}{1}(1:3)

%%%%%%%%%%%
% Function
%%%%%%%%%%%

[mean,stdev]=stat(y)

%%%%%%%%%%%%%%%%%%%
% function handles
%%%%%%%%%%%%%%%%%%%

% reference to an anonymous function
my_sqrt = @(x) x.^0.5;
my_sqrt(25)

% function handle array 
trigFun = {@sin, @cos, @tan};
plot(trigFun{2}(-pi:0.01:pi))

figure
for i=1:size(trigFun,2)
    subplot(size(trigFun,2),1,i)
    plot(trigFun{i}(-pi:0.01:pi))
end

% other tips
% Deleting Rows and Columns
% You can delete rows and columns from a matrix using just 
% a pair of square brackets. 

X = magic(3)
% to delete the second column of X
X(:,2) = []

%%%%%%%%%%%%%%%%%%%%%%
% Loops and Relations
%%%%%%%%%%%%%%%%%%%%%%

a
comp1=[a>2] % look at the results, it's an array

% for loop, examples
% loop 1
x1 = []; n=5;
for i = 1:n 
  x1=[x1,i^2];
end

% loop 2
x2 = []; n=5;
for i = 1:n 
  x2(i)=i^2;
end
x2

% loop 3
x3 = []; n=5;
for i = n:-1:1 
  x3(i)=i^2;
end
x3

%loop 4
x4 = []; n=5;
for i = n:-1:1 
  x4(n-i+1)=i^2;
end
x4

% Matrix comparison

A = [1 2 3; 4 5 6; 7 8 9];
B = [1 10 11; 12 13 14; 15 16 17];

A==B % check WHERE A and B are equals

if A~=B
   i=0;
else
   i=1;
end
% i = ? 

if any(any(A~=B))
   i=0;
else
   i=1;
end
%i = ?

help any

isequal(A,B) % check if the two matrices are actually equal

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Loops vs Built-in Functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Built-in functions are optimized and much faster
% than loops

% random: 100000 random values in [0,1]
random=rand(1,100000);

% create a new array with elements in random > 0.5
% solution 1: for loop
tic
real1=[];
j=1;
n=size(random,2)
for i=1:n
    if (random(i)>0.5)
       real1(j)=random(i);
       j=j+1;
    end
end
toc

% solution 2: built-in function (find)
tic
real2=[];
real2=random(find(random(:)>0.5));
toc

help find

% Logical Subscription
x = [2.1 1.7 1.6 1.5 NaN 1.9];
y=x(isfinite(x))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Vectorizing and Preallocation 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;

% creates a matrix by computation
A = rand(1000); 
B = rand(1000); 
tic
for j = 1:1000 
   for k = 1:1000; 
     X1(j,k) = sqrt(A(j,k)) * B(j,k); 
   end 
end
sum1=sum(sum(X1))
toc

%using vectorization
tic
X2 = sqrt(A).* B;
sum2=sum(sum(X2))
toc

%using preallocation
tic
X3=zeros(1000,1000);
for j = 1:1000 
   for k = 1:1000; 
     X3(j,k) = sqrt(A(j,k)) * B(j,k); 
   end 
end
sum3=sum(sum(X3))
toc

clear all

%%%%%%%%
% Plot
%%%%%%%%

% family of sine curves
x=[0:.2:20];
y=sin(x)./sqrt(x+1);
y(2,:)=sin(x/2)./sqrt(x+1);
y(3,:)=sin(x/3)./sqrt(x+1);
plot(x,y)
% legend
legend('signal 1','signal 2','signal 3')

% interactive plotting
% restores the configuration of tools the last time you were using them
plottools % eg: area 

% random: 100000 random values in [0,1]
random=rand(1,100000);
% does the random generator work well?
hist(random)

close all
clear all

% subplots
x = 0:pi/100:2*pi;
y1 = sin(x);
y2 = sin(x+.25); % sum a constant to each element of the array
y3 = sin(x+.5);
subplot(2,1,1); % first subplot
plot(x,y1,x,y2,x,y3);
axis tight; % sets the axis limits to the range of the data
w1 = cos(x);
w2 = cos(x+.25);
w3 = cos(x+.5);
subplot(2,1,2);
plot(x,w1,x,w2,x,w3);
axis tight; 

%%%%%%%%%%%%
% 3D Plots
%%%%%%%%%%%%

% PEAKS is a function of two variables, obtained by translating and
% scaling Gaussian distributions
%Z = peaks(20);
%Z = rand(20);
Z = magic(20);
figure(1); % Select window 
h = surf(Z); % 3D colored surface
colormap hot % set colormap look up table
shading interp % controls the color shading
xlabel('X Axis') % annotate the graph
ylabel('Y Axis')
zlabel('Function Value')
title('My First 3D Plot')

% Using Meshgrid
% what the function meshgrid does is
a=[1 2 3];
[X Y]=meshgrid(a)
clear X Y

% consider the sin(r)/r or sinc function. 
[X,Y] = meshgrid(-8:.5:8);
%The matrix R contains the distance from the center of the matrix
R = sqrt(X.^2 + Y.^2) + eps; % Adding eps prevents the divide by zero
Z = sin(R)./R; % values of the sinc function
% plots the colored parametric mesh defined by four matrix arguments
mesh(X,Y,Z) % C=Z, so color is proportional to mesh height

% lights and view adjustment
figure;
surf(X,Y,Z,'FaceColor','interp',...
	'EdgeColor','none',...
	'FaceLighting','phong')
daspect([5 5 1]) % data aspect ratio
axis tight
view(-50,30)
camlight left % create or set position of a light

% Sphere
sphere % generates a sphere consisting of 20-by-20 faces.
n=5
sphere(n) % draws a surf plot of an n-by-n sphere in the current figure.

%%%%%%%%%%%%%%%%%
% CREATING A GUI
%%%%%%%%%%%%%%%%%

guide
%guide astroneural

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Example: Image Processing Toolbox
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

coin_count

%%%%%%%%%%%%%
% Profiling
%%%%%%%%%%%%%

%help profile
profile on
% plot(magic(35))
coin_count
profile viewer
profsave(profile('info'),'profile_results') % save a static version of the HTML profile report

% run desktop -> profiler

% help m_lint
% mlint matlab_class_1
mlint crea_mat

clear all

%%%%%%%%%%%
% Examples
%%%%%%%%%%%

% a scalar is subtracted from a matrix by subtracting it from each element
A=[1 2 3; 4 5 6];
A-2

% use find to build a new array
load artifacts.dat
% true_obj: data with 1 in the 5th column
true_obj=artifacts(find(artifacts(:,5)==1),:);

% how to easily count NaNs
my_data=[1 0.5 NaN 90 NaN]; 
my_NaNCount = sum(isnan(my_data))

% Simple data analysis
% count.dat is 24-by-3 array containing hourly traffic counts (rows) 
% at three intersections (columns) for a single day.
load count.dat
% finding outliers at intersection 3
i3=count(:,3);
bin_counts = hist(i3); % Histogram bin counts
N = max(bin_counts); % Maximum bin count
mu3 = mean(i3); % Data mean
sigma3 = std(i3); % Data standard deviation
hist(i3) % Plot histogram
hold on
plot([mu3 mu3],[0 N],'r','LineWidth',2) % Mean
X = repmat(mu3+(1:2)*sigma3,2,1);
Y = repmat([0;N],1,2);
plot(X,Y,'g','LineWidth',2) % Standard deviations
legend('Data','Mean','Stds') 
hold off
% some of the data are more than two standard deviations above the mean
% If you identify these data as errors (not features), 
% replace them with NaN values as follows:
outliers = (i3 - mu3) > 2*sigma3;
i3m = i3; % Copy c3 to c3m
i3m(outliers) = NaN; % Add NaN values

% shape of a distribution
figure
hist(count)
legend('Intersection 1','Intersection 2','Intersection 3')

% polynomial curve fitting
clear all
close all
x = [1 2 3 4 5]; y = [5.5 43.1 128 290.7 498.4];
% A third degree polynomial that approximately fits the data is
p = polyfit(x,y,3)
% Compute the values of the polyfit estimate over a finer range, 
% and plot the estimate over the real data values for comparison:
x2 = 1:.1:5;
y2 = polyval(p,x2);
plot(x,y,'o',x2,y2)
grid on