% This script shows examples from 
% the topics covered in the Matlab lesson 1
% 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

% 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

%%%%%%%%%%%%%%%%%%%%%%
% 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))

% Function
[mean,stdev]=stat(y)