%% MATLAB Recipes for Earth Sciences - Chapter 9

%% Section 9.2
clear

%%
data = load('sediments_1.txt');

%%
subplot(2,1,1)
plot(data(:,1:3)), grid
legend('MinA1','MinA2','MinA3')
xlabel('Sample ID'), ylabel('Quantity')
subplot(2,1,2)
plot(data(:,4:6)), grid
legend('MinB1','MinB2','MinB3')
xlabel('Sample ID'), ylabel('Quantity')

%%
for i = 1 : 30
 sample(i,:) = [sprintf('%02.0f',i)];
end
minerals = ['MinA1';'MinA2';'MinA3';'MinB1';'MinB2';'MinB3'];

%%
corrmatrix = corrcoef(data);
corrmatrix = flipud(corrmatrix);

%%
imagesc(corrmatrix), colormap(hot), caxis([-1 1])
title('Correlation Matrix')
axis square, colorbar, hold
set(gca,'XTickLabel',minerals,'YTickLabel',flipud(minerals))

%%
[pcs,newdata,variances] = princomp(data);

%%
pcs(:,1:5)

%%
subplot(2,2,1), plot(1:6,pcs(:,1),'o'), axis([1 6 -1 1])
text((1:6)+0.2,pcs(:,1),minerals,'FontSize',8), hold
plot(1:6,zeros(6,1),'r'), title('PC 1')
subplot(2,2,2), plot(1:6,pcs(:,2),'o'), axis([1 6 -1 1])
text((1:6)+0.2,pcs(:,2),minerals,'FontSize',8), hold
plot(1:6,zeros(6,1),'r'), title('PC 2')
subplot(2,2,3), plot(1:6,pcs(:,3),'o'), axis([1 6 -1 1])
text((1:6)+0.2,pcs(:,3),minerals,'FontSize',8), hold
plot(1:6,zeros(6,1),'r'), title('PC 3')
subplot(2,2,4), plot(1:6,pcs(:,4),'o'), axis([1 6 -1 1])
text((1:6)+0.2,pcs(:,4),minerals,'FontSize',8), hold
plot(1:6,zeros(6,1),'r'), title('PC 4')

%%
plot(pcs(:,1),pcs(:,2),'o'), hold on
text(pcs(:,1)+0.02,pcs(:,2),minerals,'FontSize',14)
plot([-0.8 0.8],[0 0],'r')
plot([0 0],[-0.8 0.8],'r')
xlabel('First Principal Component Loads')
ylabel('Second Principal Component Loads')
hold off

%%
plot(newdata(:,1),newdata(:,2),'+'), hold on
text(newdata(:,1)+0.01,newdata(:,2),sample)
plot([-80 100],[0 0],'r')
plot([0 0],[-60 80],'r')
xlabel('First Principal Component Scores')
ylabel('Second Principal Component Scores')
hold off

%%
percent_explained = 100*variances/sum(variances)

%%
corrcoef(newdata(:,1),newdata(:,2))

%%
plot(1:30,newdata(:,1),1:30,newdata(:,2)), grid, legend('PC1','PC2')
xlabel('Sample ID'), ylabel('Value')

%% Section 9.3
clear

%%
i = (1:0.01:10 * pi)';
[dummy index] = sort(sin(i));

%%
s1(index,1) = i/31; s1 = s1 - mean(s1);
s2 = abs(cos(1.89*i)); s2 = s2 - mean(s2);
s3 = sin(3.43*i);

%%
subplot(3,2,1), plot(s1), ylabel('s_1'), title('Raw signals')
subplot(3,2,3), plot(s2), ylabel('s_2')
subplot(3,2,5), plot(s3), ylabel('s_3')

%%
randn('state',1);

%%
x = [.1*s1 + .8*s2 + .01*randn(length(i),1), ...
 .4*s1 + .3*s2 + .01*randn(length(i),1), ...
 .1*s1 + s3 + .02*randn(length(i),1)];

%%
subplot(3,2,2), plot(x(:,1)), ylabel('x_1'), title('Mixed signals')
subplot(3,2,4), plot(x(:,2)), ylabel('x_2')
subplot(3,2,6), plot(x(:,3)), ylabel('x_3')

%%
[E sPCA D] = princomp(x);
sPCA = sPCA./repmat(std(sPCA),length(sPCA),1);

%%
subplot(3,2,1), plot(sPCA(:,1))
ylabel('s_{PCA1}'), title('Separated signals - PCA')
subplot(3,2,3), plot(sPCA(:,2)), ylabel('s_{PCA2}')
subplot(3,2,5), plot(sPCA(:,3)), ylabel('s_{PCA3}')

%%
A_PCA = E * sqrt(diag(D));
W_PCA = inv(sqrt(diag(D))) * E';

%%
rand('state',1);

%%
div = 0;
B = orth(rand(3, 3) - .5);
BOld = zeros(size(B));

%%
while (1 - div) > eps
 B = B * real(inv(B' * B)^(1/2));
 div = min(abs(diag(B' * BOld))); 
 BOld = B;
 B = (sPCA' * (sPCA * B) .^ 3) / length(sPCA) - 3 * B;
 sICA = sPCA * B;
end

%%
subplot(3,2,2), plot(sICA(:,1)), ylabel('s_{ICA1}'),
 title('Separated signals - ICA')
subplot(3,2,4), plot(sICA(:,2)), ylabel('s_{ICA2}')
subplot(3,2,6), plot(sICA(:,3)), ylabel('s_{ICA3}')

%%
A_ICA = A_PCA * B;
W_ICA = B' * W_PCA;

%% Section 9.4
clear

%%
data = load('sediments_2.txt');
Y = pdist(data);

%%
imagesc(squareform(Y)), colormap(hot)
title('Euclidean distance between pairs of samples')
xlabel('First Sample No.')
ylabel('Second Sample No.')
colorbar

%%
Z = linkage(Y)

%%
dendrogram(Z);
xlabel('Sample No.')
ylabel('Distance')
box on

%%
cophenet(Z,Y)