% Copyright (c) 2017, Domenico L. Gatti
% All rights reserved.
% 
% Redistribution and use in source and binary forms, with or without 
% modification, are permitted provided that the following conditions are 
% met:
% 
%     * Redistributions of source code must retain the above copyright 
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright 
%       notice, this list of conditions and the following disclaimer in 
%       the documentation and/or other materials provided with the 
%       distribution
%       
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
% IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
% CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
% EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
% PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
% PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
% NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%
%% General dependencies
% We always start from the CODE directory and we add to the path
% subdirectories containing various tools described in the book chapters.
addpath(genpath('../GENERAL_SCRIPTS_FUNCTIONS'));
addpath(genpath('../DATABASE'));
% addpath(genpath('../TOOLBOXES'));

%% CHAPTER 10: Principal component analysis

%% Covariance matrix
close all
clear

nvars = 10;
nobs = 20;
X = rand(nvars,1)
X_mat = zeros(nvars,nobs)
for i = 1:nobs
    X_mat(:,i) = normrnd(X,0.05)
end
X_mean = mean(X_mat,2)
cX_mat = X_mat - X_mean(:,ones(1,nobs))
% or simply
cX_mat = X_mat - X_mean

X1_var = cX_mat(1,:)*cX_mat(1,:)'/(nobs-1)
X1_var = var(cX_mat(1,:))
X12_var = cX_mat(1,:)*cX_mat(2,:)'/(nobs-1)

cov_X = cX_mat*cX_mat'/(nobs-1)
cov_X = cov(cX_mat')
cov_X = cov(cX_mat',1) % using nobs instead of nobs-1
imagesc(cov_X)

a = [5 1]'
b = [3 5]'

a'*a
sqrt(a'*a)
norm(a)
norm(b)

cos_theta = (a'*b)/(norm(a)*norm(b))
theta = (180/pi)*acos(cos_theta)
theta = acosd(cos_theta)

ca = a - mean(a);
cb = b - mean(b);
cos_theta = (ca'*cb)/(norm(ca)*norm(cb))
corr_ab = corr(a,b)
theta = acosd(cos_theta)

% Notice:
% angle = 0   --> full correlation
% angle = 180 --> full anticorrelation
% angle = 90  --> no correlation

cov([a a],1) % covariance
var(a,1)     % variance
(ca'*ca)/2   % variance = sum((a-mean(a))*(a-mean(a)))/nobs

% bisectant of 1st quadrant
v = [1 1]'
% left null space of v
ov = null(v','r')
% Projection matrix
P = ov*inv(ov'*ov)*ov'
ca = P*a
cb = P*b

% Dimensionality reduction upon centering
A = [2 4 6 8;2 1 1 4; 2 3 2 5;2 2 4 4 ;2 4 2 4]
meanA = mean(A,2);
cA = A - repmat(meanA,1,4)
% cA = A - meanA
rank(A)
rank(cA)

rank(X_mat)
rank(cX_mat)
rank(cov_X)
R = chol(cov_X)
R'*R

% but
nvars = 20;
nobs = 20;
X = rand(nvars,1)
X_mat = zeros(nvars,nobs)
for i = 1:nobs
    X_mat(:,i) = normrnd(X,0.05)
end
X_mean = mean(X_mat,2)
cX_mat = X_mat - X_mean
cov_X = cov(cX_mat')

rank(X_mat)
rank(cX_mat)
rank(cov_X)
R = chol(cov_X)
R'*R
rcond(cov_X)

%%
close all
clear

x_norm = normrnd(zeros(1,1000),2.5);
y_norm = normrnd(zeros(1,1000),5.0);
data = [x_norm;y_norm];
 
PCA_1 = figure;
set(PCA_1,'Unit','Normalized','Position', [0 0.2 0.6 0.8]);
plot(data(1,:),data(2,:),'.','Linewidth',0.5,'MarkerEdgeColor','b',...
    'MarkerSize',15,'MarkerFaceColor','g')
xlabel('X');ylabel('Y');
grid on
box on
set(gca,'Ylim',[-15 15],'XLim',[-5 +5])
axis equal
hold on

[nvar,nobs] = size(data);

% Here we really don't need to center the data to calculate the covariance
% matrix because we generated all the points around 0, but for sake of
% mathematical rigor we do it anyway.
mean_data = mean(data,2);
c_data = data - mean_data(:,ones(1,nobs)); 
% c_data = data - mean_data; 
C = (c_data*c_data')/(nobs-1)
cov(data')

s_data_1 = sort(data(1,:),2);
s_data_2 = sort(data(2,:),2);
s_data = [s_data_1;s_data_2];
cov(s_data')

% Eigen decomposition of the covariance matrix
[S,D] = eig(C);

k = size(S,2);
[~,i] = sort(diag(D),'descend');
W = S(:,i);
R = D(i,i);
 
 for n=1:k
     W(:,n)=sign(mean(W(:,n)))*W(:,n);
 end

[S,D] = pcacov(C)

line(10*[0 S(1,1)],10*[0 S(1,2)],'LineWidth',4,'Color','red');
line(10*[0 S(2,1)],10*[0 S(2,2)],'LineWidth',4,'Color','green');

scores = S'*data;
cov(scores')

%%
close all

v1 = [1;-1];
v2 = [1;1];
v1 = v1/norm(v1);
v2 = v2/norm(v2);
cos_theta = v1'*v2;
acos(cos_theta)*180/pi;

% Each point in the 2D space is going to be a linear combination of the
% new axes v1 and v2: 
data = [v1 v2]*[x_norm;y_norm];

PCA_2 = figure;
set(PCA_2,'Unit','Normalized','Position', [0 0.2 0.6 0.8]);
plot(data(1,:),data(2,:),'.','Linewidth',0.5,'MarkerEdgeColor','b',...
    'MarkerSize',15,'MarkerFaceColor','g')
xlabel('X');ylabel('Y');
grid on
box on
set(gca,'Ylim',[-15 15],'XLim',[-5 +5])
axis equal
hold on

C = cov(data')

[S,D] = pcacov(C);
line(10*[0 S(1,1)],10*[0 S(1,2)],'LineStyle','-','LineWidth',4,'Color','red');
line(10*[0 S(2,1)],10*[0 S(2,2)],'LineStyle','-','LineWidth',4,'Color','green');
hold off

scores = S'*data;
cov(scores')
D

mean_data = mean(data,2);
c_data = data - mean_data(:,ones(1,nobs)); 
C = (c_data*c_data')/(nobs-1)
scores = S'*c_data;
cov(scores')

% figure;plot(scores(1,:),scores(2,:),'.','Linewidth',0.5,'MarkerEdgeColor','b',...
%     'MarkerSize',15,'MarkerFaceColor','g')
% xlabel('S1');ylabel('S2');
% grid on
% box on
% axis equal
% hold off

% Dimensionality reduction
one_dim_data = S(:,1)*S(:,1)'*c_data + mean_data(:,ones(1,nobs));
PCA_3 = figure;
set(PCA_3,'Unit','Normalized','Position', [0 0.2 0.6 0.8]);
plot(one_dim_data(1,:),one_dim_data(2,:),'.','Linewidth',0.5,'MarkerEdgeColor','b',...
    'MarkerSize',15,'MarkerFaceColor','g')
xlabel('X');ylabel('Y');
grid on
box on
set(gca,'Ylim',[-15 15],'XLim',[-5 +5])
axis equal
hold on

%% A 3-D example
% clear
% close all
% nobs = 100;
% x_norm = normrnd(zeros(1,nobs),5);
% y_norm = normrnd(zeros(1,nobs),1.0);
% z_norm = normrnd(zeros(1,nobs),0.1);
% a = (0.05:0.05:5) ; b = (0.1:0.1:10) ; c = (0.1:0.1:10);
% A = [a + x_norm; b + y_norm; c + z_norm];


% dlmwrite('../DATABASE/PCA_3D_data_1.txt',A,'\t');

%%
A = dlmread('../DATABASE/PCA_3D_data.txt');
[nvar,nobs] = size(A);

plot3(A(1,:),A(2,:),A(3,:),'o','Linewidth',0.5,'MarkerEdgeColor','b',...
    'MarkerFaceColor','g')
xlabel('X');ylabel('Y');zlabel('Z');
grid on
box on
hold off

C = cov(A')

%%
close all
 
[S,D] = pcacov(C);

% Next we calculate the 'scores' (the representation of the centered data 
% in the eigenvector basis):

mean_A = mean(A,2);
cA = A - mean_A(:,ones(1,nobs));

% Centered data using a projection matrix
es = ones(1,nobs);
plane = null(es);
P = plane/(plane'*plane)*plane';
cA_p = P*A';

cA_s = S'*cA;

%%
% For comparison we plot the original ...

figure
scatter3(A(1,:),A(2,:),A(3,:));
hold on
x = -2:1:12;
y = zeros(1,15);
z = zeros(1,15);
plot3(x,y,z,'--r');
x = zeros(1,15);
y = -2:1:12;
z = zeros(1,15);
plot3(x,y,z,'--r');
x = zeros(1,15);
y = zeros(1,15);
z = -2:1:12;
plot3(x,y,z,'--r');
grid on
box on
axis equal
set(gca,'XLim',[-12 20],'YLim',[-2 13],'ZLim',[-2 13]);

xlabel('x data','FontSize',12,'FontWeight','n');
ylabel('y data','FontSize',12,'FontWeight','n');
zlabel('z data','FontSize',12,'FontWeight','n');

%%
% and the centered data
figure;
scatter3(cA(1,:),cA(2,:),cA(3,:));
box on
hold on
x = -10:1:10;
y = zeros(1,21);
z = zeros(1,21);
plot3(x,y,z,'--r');
x = zeros(1,13);
y = -6:1:6;
z = zeros(1,13);
plot3(x,y,z,'--r');
x = zeros(1,9);
y = zeros(1,9);
z = -4:1:4;
plot3(x,y,z,'--r');

%%
% Here we also plot the three eigenvectors. In order to do this we resort 
% to the parametric vector equation of a line:
% L(t) = P1 + t*V  ; where V = P2 - P1
% In many cases P1 = [x1 y1 z1] = [0 0 0] and V = P2 = [x2 y2 z2]:
% Substituting we obtain the parametric form of a line intersecting two 
% points P1 [x1 y1 z1] and P2 [x2 y2 z2]:
% L(t) = P1 + t*V = P1 + t*(P2-P1)
% L(t) = (1-t)*P1 + t*P2 ; clearly P1 can be the origin [0 0 0] or any 
% other point (for example the centroid of the points cloud). 

origin = zeros(size(cA,1),1);

% Since we have plotted centered data, in our case we can simply use:
% L(t) = P1 + t*V where P1 is the origin at [0 0 0] and V is the 
% eigenvector.
% We can make the eigenvectors appear longer or shorter at will by either
% decreasing the sampling of t of by decreasing t itself:

for i = 1:15
    t = i-8;
    t1 = t;
    t2 = t*0.8;
    t3 = t*0.6;
pc1(i,:) = origin + t1*S(:,1);
pc2(i,:) = origin + t2*S(:,2);
pc3(i,:) = origin + t3*S(:,3);
end
plot3(pc1(:,1),pc1(:,2),pc1(:,3),'-b');
plot3(pc2(:,1),pc2(:,2),pc2(:,3),'-m');
plot3(pc3(:,1),pc3(:,2),pc3(:,3),'-g');
grid on
axis equal
set(gca,'XLim',[-13 16],'YLim',[-10 10],'ZLim',[-6 6]);

xlabel('centered x data','FontSize',12,'FontWeight','n');
ylabel('centered y data','FontSize',12,'FontWeight','n');
zlabel('centered z data','FontSize',12,'FontWeight','n');

%%
% Finally we plot the centered data in the eigenvector basis
figure;
scatter3(cA_s(1,:),cA_s(2,:),cA_s(3,:));
box on
grid on
hold on
%%

rcA_s = [cA_s(1:2,:) ; zeros(1,nobs)];

% Notice that we have lost one dimension (the z coordinate, which is zero
% everywhere now). 

scatter3(rcA_s(1,:),rcA_s(2,:),rcA_s(3,:),'or');figure(gcf)
axis equal
% set(gca,'XLim',[-16 18],'YLim',[-10 10],'ZLim',[-2 +2]);

xlabel('pc 1 data','FontSize',12,'FontWeight','n');
ylabel('pc 2 data','FontSize',12,'FontWeight','n');
zlabel('pc 3 data','FontSize',12,'FontWeight','n');
hold off

%%
% Next we carry out a change of basis from the eigenvector basis back to
% the standard basis and we add back the mean of each column:

rcA = S*rcA_s 
rA = rcA + mean_A(:,ones(1,nobs))

figure;
scatter3(A(1,:),A(2,:),A(3,:));figure(gcf)
hold on

% Plot the axes
x = -2:1:12;
y = zeros(1,15);
z = zeros(1,15);
plot3(x,y,z,'--r');
x = zeros(1,15);
y = -2:1:12;
z = zeros(1,15);
plot3(x,y,z,'--r');
x = zeros(1,15);
y = zeros(1,15);
z = -2:1:12;
plot3(x,y,z,'--r');

%%
% Plot the three eigenvectors at the center of the data cloud. 
for i = 1:size(A,1)
origin(i,1) = mean(A(i,:));
end
% Since in our case origin shift does not change P2-P1, we can simply use:
% L(t) = P1 + t*V where P1 is the new origin (the centroid of the points) 
% and V is the eigenvector.
for i = 1:15
    t = i-8;
    t1 = t;
    t2 = t*0.8;
    t3 = t*0.6;
pc1(i,:) = origin + t1*S(:,1);
pc2(i,:) = origin + t2*S(:,2);
pc3(i,:) = origin + t3*S(:,3);
end
plot3(pc1(:,1),pc1(:,2),pc1(:,3),'-b');
plot3(pc2(:,1),pc2(:,2),pc2(:,3),'-m');
plot3(pc3(:,1),pc3(:,2),pc3(:,3),'-g');

%%
% Plot the reduced data
scatter3(rA(1,:),rA(2,:),rA(3,:),'or');figure(gcf)
grid on
box on
axis equal
set(gca,'XLim',[-12 20],'YLim',[-2 13],'ZLim',[-2 13]);
xlabel('x data','FontSize',12,'FontWeight','n');
ylabel('y data','FontSize',12,'FontWeight','n');
zlabel('z data','FontSize',12,'FontWeight','n');
hold off

%% PCA as a projection operation

% Projection matrix:
P = S(:,1:2)*((S(:,1:2)'*S(:,1:2))\S(:,1:2)');
P = S(:,1:2)*S(:,1:2)'
rA_p = P*A;

% Plot the original data
figure
scatter3(A(1,:),A(2,:),A(3,:),'*g');
grid on
box on
axis equal
set(gca,'XLim',[-12 20],'YLim',[-2 13],'ZLim',[-2 13]);
xlabel('x data','FontSize',12,'FontWeight','n');
ylabel('y data','FontSize',12,'FontWeight','n');
zlabel('z data','FontSize',12,'FontWeight','n');
hold on

% Plot the dimension reduced data 
scatter3(rA(1,:),rA(2,:),rA(3,:),15,'*',...
    'Linewidth',0.5,'MarkerEdgeColor','r',...
    'MarkerFaceColor','w');figure(gcf);figure(gcf)
grid on
box on
axis equal
set(gca,'XLim',[-12 20],'YLim',[-2 13],'ZLim',[-2 13]);
xlabel('x data','FontSize',12,'FontWeight','n');
ylabel('y data','FontSize',12,'FontWeight','n');
zlabel('z data','FontSize',12,'FontWeight','n');

% Plot the projected data
scatter3(rA_p(1,:),rA_p(2,:),rA_p(3,:),'*k');figure(gcf)
axis equal
set(gca,'XLim',[-12 20],'YLim',[-2 13],'ZLim',[-2 13]);
xlabel('x data','FontSize',12,'FontWeight','n');
ylabel('y data','FontSize',12,'FontWeight','n');
zlabel('z data','FontSize',12,'FontWeight','n');

%%
origin = mean_A;
shift = origin-P*origin;
rA_ps = rA_p + (shift(:,ones(1,nobs)));
scatter3(rA_ps(1,:),rA_ps(2,:),rA_ps(3,:),80,'s',...
    'Linewidth',0.5,'MarkerEdgeColor','m',...
    'MarkerFaceColor','n');figure(gcf)
hold off

%% Fluorescence changes produced by calcium binding to Atp11p
close all
clear, clc
importfile('../DATABASE/Fluorescence_data.txt');
A = Fluorescence_data;
% openvar A
ca_conc=[0 2.6 7.8 17.8 38.8 80.8 163.8 331 664 ...
1331 2664 3997 5330 6663 7996]; wl =[300:400];
[XI,YI] = meshgrid(ca_conc,wl); 
Fluor_A = figure;
axes0 = axes('Parent',Fluor_A,'YDir','reverse');
view(axes0,[-77.5 12]);
box(axes0,'on');
grid(axes0,'on');
hold(axes0,'all');
surf(XI,YI,A,'Parent',axes0);
[nwls,nobs] = size(A);
xlabel('[Ca], micromolar'), ylabel('wavelength, nm');box on

mean_A = mean(A,2);
cA = A - mean_A(:,ones(1,nobs));
% cA = A - mean_A;
rank(A), rank(cA)

% Centering of the data by projection
equisect = ones(1,15);
H = null(equisect)
size(H)
P = H/(H'*H)*H'
p_cA = (P*A')'

Centered_A = figure;
axes1 = axes('Parent',Centered_A,'YDir','reverse');
view(axes1,[-77.5 12]);
box(axes1,'on');
grid(axes1,'on');
hold(axes1,'all');
surf(XI,YI,cA,'Parent',axes1);
xlabel('[Ca], micromolar'), ylabel('wavelength, nm');box on


Mean_A = figure;
axes2 = axes('Parent',Mean_A,'YDir','reverse');
view(axes2,[-77.5 12]);
box(axes2,'on');
grid(axes2,'on');
hold(axes2,'all');
surf(XI,YI,mean_A(:,ones(1,nobs)),'Parent',axes2);
xlabel('[Ca], micromolar'), ylabel('wavelength, nm');box on

% PCA
C = (cA*cA')/(nobs-1);
rank(cA)
openvar C
cov_A = cov(A');
openvar cov_A
[S,D,percent] = pcacov(cov_A);
openvar D

Scree = figure;
% semilogy(D(1:14)/max(D),'-ro',...
%     'Linewidth',1.0,'MarkerEdgeColor','m',...
%     'MarkerFaceColor','c'); hold on
semilogy(percent(1:14),'-ro',...
    'Linewidth',1.0,'MarkerEdgeColor','m',...
    'MarkerFaceColor','c'); 
ylabel('log(% variance)  ')
xlabel('eigenvalue index  ')
xlim([0 15]),ylim([1e-4,200])
grid on

% The covariance of the data is completely wiped out when the original data 
% or the centered data are represented in eigenvector basis
% A_s = S' * A;
% cov_A_s = cov(A_s');
% openvar cov_A_s
% cov_A_s(cov_A_s<=1e-9) = 0;
% 
% A_data_s = figure;
% axes1 = axes('Parent',A_data_s,'YDir','reverse');
% view(axes1,[-77.5 12]);
% box(axes1,'on');
% grid(axes1,'on');
% hold(axes1,'all');
% nS = nwls;
% [XS,YS] = meshgrid(ca_conc,[1:nS]); 
% surf(XS,YS,A_s,'Parent',axes1);
% xlabel('[Ca], micromolar'), ylabel('S index');box on
% ylim([1,nS])

cA_s = S' * cA;
cov_cA_s = cov(cA_s');
f_cov_cA_s = fix(cov_cA_s);
openvar f_cov_cA_s

Normal_modes = figure;
plot(wl,S(:,1),'b');
hold on
plot(wl,S(:,2),'r');
plot(wl,S(:,3),'g');
xlabel('Wavelength (nm)  '),ylabel('S units'  )

cA_data_s = figure;
axes1 = axes('Parent',cA_data_s,'YDir','reverse');
view(axes1,[-77.5 12]);
box(axes1,'on');
grid(axes1,'on');
hold(axes1,'all');
nS = nwls;
[XS,YS] = meshgrid(ca_conc,[1:nS]); 
surf(XS,YS,cA_s,'Parent',axes1);
xlabel('[Ca], micromolar'), ylabel('S index');box on
ylim([1,nS])

Scores_cA = figure;
plot(ca_conc, cA_s(1,:),'-ob', 'Linewidth',1.0,...
'MarkerEdgeColor','b','MarkerFaceColor','c');
hold on
plot(ca_conc, cA_s(2,:),'-or', 'Linewidth',1.0,...
'MarkerEdgeColor','r','MarkerFaceColor','y');
xlabel(['[Ca] (' texlabel('mu') 'M)  '])
ylabel('SCORES'  )
hold off

% Hyperbola fit of the scores
xvec = ca_conc';
yvec = (cA_s(1,:)' - min(cA_s(1,:)))/1200;
Equil_bind = figure;
plot(xvec,yvec,'ob'); hold on

f = fittype('b*(x/(a + x))');
[Hyperb_1,GOF_1] = fit(xvec,yvec,f,...
'StartPoint',[40 1]);
u1 = coeffvalues(Hyperb_1)
lsq_1 = u1(2)*xvec./(u1(1)+xvec);
plot(xvec,lsq_1,'--m');

f = fittype('b*(x/(a + x)) + d*(x/(c + x))');
[Hyperb_2,GOF_2] = fit(xvec,yvec,f,...
'StartPoint',[20 0.5 2000 0.5]);
u2 = coeffvalues(Hyperb_2)
lsq_21 = u2(2)*xvec./(u2(1)+xvec);
lsq_22 = u2(4)*xvec./(u2(3)+xvec);
plot(xvec,lsq_21,'--b',xvec,lsq_22,'--r')
lsq_2 = lsq_21 + lsq_22;
plot(xvec,yvec,'ob',xvec,lsq_2,'-c')
xlabel(['[Ca] (' texlabel('mu') 'M)  '])
ylabel('Fractional Saturation'  )

% Bring back the data in the original space removing the noise
fA = S(:,1)*cA_s(1,:) + mean_A(:,ones(1,nobs));
Filtered_A = figure;
axes3 = axes('Parent',Filtered_A,'YDir','reverse');
view(axes3,[-77.5 12]);
box(axes3,'on');
grid(axes3,'on');
hold(axes3,'all');
surf(XI,YI,fA,'Parent',axes3);
xlabel('[Ca], micromolar'), ylabel('wavelength, nm');box on

% PCA filtering as a least squares operation
S1 = S(:,1)
P = S1/(S1'*S1)*S1';
P = S1*inv(S1'*S1)*S1';
P = S1*S1';
p_fA = P*cA; % 'part' of S1 that looks most like cA

Projection_fA = figure;
axes4 = axes('Parent',Projection_fA,'YDir','reverse');
view(axes4,[-77.5 12]);
box(axes4,'on');
grid(axes4,'on');
hold(axes4,'all');
surf(XI,YI,p_fA,'Parent',axes4);
xlabel('[Ca], micromolar'), ylabel('wavelength, nm');box on

% Least squares solution
p_fA = S1*inv(S1'*S1)*S1'*cA;
x = inv(S1'*S1)*S1'*cA % solution for S1*x = cA
x = S1'*cA 

LS_Scores = figure;
plot(ca_conc, x,'-ob', 'Linewidth',1.0,...
'MarkerEdgeColor','b','MarkerFaceColor','c');
xlabel(['[Ca] (' texlabel('mu') 'M)  '])
xlabel('[Ca] (\muM)')
ylabel('SCORES'  )

%% Essential dynamics of Atp12p
% To run this section we need to download the MDTOOLBOX developed by by
% Yasuhiro Matsunaga (https://github.com/ymatsunaga/mdtoolbox), which
% contains a very sophisticated set of tools for the analysis of MD
% simulations. After downloading the toolbox and installing it in our
% TOOLBOXES directory, we can add it to our path:
addpath(genpath('TOOLBOXES/MDTOOLBOX'))

close all
clear, clc
trj_mat = dlmread('../DATABASE/md_trj_matrix.txt');
savefile = 'MD_trj_pca';
[nframes,ncoords] = size(trj_mat);
natoms = ncoords/3;


[Aref,A] = meanstructure(trj_mat);
Aref = Aref'; A = A';
[nvar,nobs] = size(A);

a1 = A(:,1);
a1_x = a1(1:3:end);
a1_y = a1(2:3:end);
a1_z = a1(3:3:end);
 
A1_Frame = figure;
plot3(a1_x,a1_y,a1_z,'-bo','Linewidth',1.0,...
'MarkerEdgeColor','r','MarkerFaceColor','y')
box('on')
grid('on')
xlabel('X coord')
ylabel('Y coord')
zlabel('Z coord')

mean_A = mean(A,2);
cA = A - mean_A(:,ones(1,nobs));
C = (cA*cA')/(nobs-1); 
C2 = (cA'*cA)/(nvar-1); 
MSSF = diag(C2)

[S,D,percent] = pcacov(C);
sum(D)
sum(MSSF)
cA_s = S' * cA;

x_S1 = S(1:3:end,1);
y_S1 = S(2:3:end,1);
z_S1 = S(3:3:end,1);
Vibration = figure; 
plot(x_S1,'r')
hold on
plot(y_S1,'g')
plot(z_S1,'b')
legend('x\_s1','y\_s1','z\_s1','Location','Best')
xlabel('Ca number')
ylabel('Coefficient')

energy = cumsum(D);
Energy = figure; plot(energy/energy(end))
xlabel('Eigenvalue index   ')
ylabel('E_j/E_t_o_t_a_l   ')

figure; plot3(cA_s(1,:),cA_s(2,:),cA_s(3,:),'ro')
xlabel('Eigenmode 1 scores   ')
ylabel('Eigenmode 2 scores   ')
zlabel('Eigenmode 3 scores   ')
grid('on');box('on')
hold on

CL = clusterdata(cA_s(1:3,:)','linkage','ward','savememory','on',...
    'maxclust',5);
scatter3(cA_s(1,:),cA_s(2,:),cA_s(3,:),20,CL,'fill')

CL1 = CL==1; CL2 = CL==2; CL3 = CL==3;CL4 = CL==4; CL5 = CL==5;
CL1_size = sum(CL1); CL2_size = sum(CL2); CL3_size = sum(CL3); CL4_size = sum(CL4);
CL5_size = sum(CL5);

rcA = S(:,1:15)*cA_s(1:15,:);
rA = rcA + mean_A(:,ones(1,nobs));

rcA_ed1 = S(:,1)*cA_s(1,:)
rA_ed1 = rcA_ed1 + mean_A(:,ones(1,nobs))

ED1 = figure;
for i = 1:29:nobs
    F = rA_ed1(:,i);
    x = F (1:3:end);
    y = F (2:3:end);
    z = F (3:3:end);
plot3(x,y,z,'-bo','Linewidth',1.0,...
'MarkerEdgeColor','r','MarkerFaceColor','y');
hold on
box('on')
grid('on')
drawnow expose; pause(.1)
end

%%
first = 1;
MSF = zeros(natoms,1);
RMSF = zeros(natoms,1);
dRMSF = zeros(natoms,1);
C_diag = diag(C);
 
first = 1;
RMSF = zeros(natoms,1);
dRMSF = zeros(natoms,1);
 
for i = 1:natoms
 second = first+1;
 last = first+2;
 % Direct determination from the RMSF formula:
 dRMSF(i) = sqrt(((cA(first,:)*cA(first,:)' + ...
                  cA(second,:)*cA(second,:)' + ...
                  cA(last,:)*cA(last,:)')/nobs));
 % or from the covariance matrix:
 MSF(i) = sum(C_diag(first:last));
 RMSF(i) = sqrt(sum(C_diag(first:last)));
 first = first+3;
end

total_RMSF = sqrt(sum(MSF)) 
total_RMSF = sqrt(sum(diag(C)))
total_RMSF = sqrt(sum(D))

C_ed1 = (rcA_ed1*rcA_ed1')/(nobs-1);
total_RMSF_ed1 = sqrt(sum(diag(C_ed1)))
total_RMSF_ed1 = sqrt(D(1))

% save(savefile);

%% Special topic: Rotational ambiguity of two-way PCA
clear, clc, close all
load('../DATABASE/Fluorescence_data.txt');
X = Fluorescence_data;
ca_conc = [0 2.6 7.8 17.8 38.8 80.8 163.8 331 664 1331 2664 3997 5330 6663 7996];
wl =[300:400];
[XI,YI] = meshgrid(ca_conc,wl);
figure;surf(XI,YI,X);

% Center the data along each row and calculate the covariance matrix
nobs = size(X,2)
mX = mean(X,2)
cX = X - mX
covX = cX*cX'/(nobs-1)
figure;surf(XI,YI,cX);

% Loading and scores (given the variables along the columns of X)
[S,D] = pcacov(covX)
scores = S'*cX;

% The eigenvalues are the variance of the data represented in eigenvector
% basis.
var_scores = var(scores')
figure;plot(wl,S(:,1))

% Noise filtered data set
fcX = S(:,1:2)*scores(1:2,:)
figure;surf(XI,YI,fcX);

% This is the same as:
S_red = S;
scores_red = scores;
S_red(:,3:4) = 0;
scores_red(3:4,:)= 0;
fcX_red = S_red*scores_red;
figure;surf(XI,YI,fcX_red);

% Here we introduce an orthogonal transformation matrix that will rotate
% the relevant vectors of the basis of the eigenvector space, for example
% by 30 degrees counterclockwise
theta = pi/6;
Q = [cos(theta) sin(theta);-sin(theta) cos(theta)];
Srot1 = S(:,1:2)*Q';

% The new basis is orthogonal
Srot1'*Srot1

% Here we recalculate the scores based on the new basis
scoresrot1 = Srot1'*cX;

% Here we calculate the percentage of the variance explained by the new
% basis
var_scores(1:2)
var_scoresrot1 = var(scoresrot1')
sum_var_scores = sum(var_scores,2)
sum_var_scoresrot1 = sum(var_scoresrot1,2)

% Noise filtered data set
fcXrot1 = Srot1*scoresrot1
figure;surf(XI,YI,fcXrot1);

% This result is not limited to an orthogonal transformation matrix, but
% applies to any invertible matrix
T = rand(2);
Srot2 = S(:,1:2)*inv(T);

% The new basis is not orthogonal we can't simply take the transpose,
% but we need to calculate a pseudoinverse
scoresrot2 = pinv(Srot2)*cX;
% Variance explained
var_scoresrot2 = var(scoresrot2')
sum_var_scoresrot2 = sum(var_scoresrot2,2)

% Noise filtered data
fcXrot2 = Srot2*scoresrot2
figure;surf(XI,YI,fcXrot2);

% Instead of rotating just the selected loading columns, here we use a
% large matrix to rotate the entire eigenvector basis
T_full = rand(101);
Srot2_full = S*inv(T_full);

% Vectors in the rotated basis are no longer orthogonal
% sqrt(diag(Srot2_full'*Srot2_full))
Srot2_full'*Srot2_full

% Here we represent the data in the new rotated oblique basis:
scoresrot2_full = Srot2_full\cX;
% Variance explained
var_scoresrot2_full = var(scoresrot2_full')
sum_var_scoresrot2_full = sum(var_scoresrot2_full,2)
var_scoresrot2_full = var(scoresrot2_full(1:2,:)')
sum_var_scoresrot2_full = sum(var_scoresrot2_full,2)

% Here we use the entire basis to reconstruct the data
fcXrot2_full = Srot2_full*scoresrot2_full;
figure;surf(XI,YI,fcXrot2_full);

% But if we use only the first 2 loadings and scores the result is junk;
% thus only the full basis provides the correct result.
scoresrot2_full_red = scoresrot2_full;
scoresrot2_full_red(3:end,:) = 0;
Srot2_full_red = Srot2_full;
Srot2_full_red(:,3:end) = 0;
fcXrot2_full_red = Srot2_full_red*scoresrot2_full_red;
figure;surf(XI,YI,fcXrot2_full_red);
fcXrot2_full_red = Srot2_full(:,1:2)*scoresrot2_full(1:2,:);
figure;surf(XI,YI,fcXrot2_full_red);
fcXrot2_full_red = Srot2_full*scoresrot2_full;
figure;surf(XI,YI,fcXrot2_full_red);

% Here we plot the two spectral components for the three cases (unrotated,
% rotated orthogonal, rotated oblique:
figure;
subplot(2,1,1);plot(wl,[S(:,1) Srot1(:,1) Srot2(:,1)]), 
legend('original basis','rotated orthog. basis ','rotated non-orthog. basis')  
ylabel('Fluorescence')
title('1st spectral component')
subplot(2,1,2);plot(wl,[S(:,2) Srot1(:,2) Srot2(:,2)]), legend show
legend('original basis','rotated orthog. basis ','rotated non-orthog. basis','Location','Best')  
xlabel('wl')
ylabel('Fluorescence')
title('2nd spectral component')

%% Loading and scores given the variables along the rows of X.
clear, clc, close all
load('../DATABASE/Fluorescence_data.txt');
X = Fluorescence_data;
ca_conc = [0 2.6 7.8 17.8 38.8 80.8 163.8 331 664 1331 2664 3997 5330 6663 7996];
wl =[300:400];
[XI,YI] = meshgrid(ca_conc,wl);

Xt = X';
nobs = size(X,2)
mXt = mean(Xt);
cXt = Xt - mXt;
cX = cXt';
covXt = cXt'*cXt/nobs;
[S2,D2] = pcacov(covXt)
loadings2 = S2(:,1:2);
scores2 = cXt*loadings2;
cfXt = scores2*loadings2';
figure;surf(XI,YI,cfXt');
corr_loadings2 = corr(loadings2)
corr_scores2 = corr(scores2)

% Varimax objective function
% D = size(S2_rot1,1)
% GAMMA = 1
% OF = sum(D*sum(S2rot1.^4,1) - GAMMA*sum(S2rot1.^2,1).^2)

% Here we rotate the first 2 columns of the loading matrix
[S2rot1, T1] = rotatefactors(loadings2,'Method','varimax');
S2rot1 = loadings2*T1;

% Here we make the table of the rotated loadings.
cnames = {'1st spectr. comp.','2nd spectr. comp.','1st spectr. comp.','2nd spectr. comp.','1st spectr. comp.','2nd spectr. comp.',};
% rnames = {'x1','x2','x3','x4','x5','x6','x7','x8','x9','x10','x11','x12','x13','x14','x15','x16','SUM Y'};

Rotated_loadings = figure('Position',[300 300 700 700]);
t = uitable(Rotated_loadings,'Data',[S2rot1(1:34,:) S2rot1(35:68,:) [S2rot1(69:101,:);0 0]],'ColumnName',cnames,'FontWeight','bold','FontSize',9,'RowStriping','off','ColumnEditable',true); 
t.Position(3) = t.Extent(3);
t.Position(4) = t.Extent(4);

% Here we recalculate the scores based on the rotated loadings:
scores2rot1 = cXt*S2rot1;

% Here we generate the filterted data
cfXtrot1 = scores2rot1*S2rot1';
figure;surf(XI,YI,cfXtrot1');
corr_scores2rot1 = corr(scores2rot1*T1')
corr_scores2rot1 = inv(T1'*T1)

% Here we plot the rotated loadings:
figure;
plot(wl,[S2rot1(:,1) S2rot1(:,2)])
legend('1st spectral component','2nd spectral component')
title('Varimax Rotation')
xlabel('wl')
ylabel('Fluorescence')

% Other rotations.
% Method Orthomax
[S2rot2, T2] = rotatefactors(loadings2,'Method','orthomax','Coeff',0);
scores2rot2 = cXt*S2rot2;
cfXtrot2 = scores2rot2*S2rot2';
figure;surf(XI,YI,cfXtrot2');
corr_scores2rot2 = corr(scores2rot2)
corr_loadings2 = corr(loadings2)
corr_loadings2rot2 = corr(S2rot2)
T2'*T2

% Method Promax
[S2rot3, T3] = rotatefactors(S2(:,1:2),'Method','promax');
scores2rot3 = cXt*S2rot3;
cfXtrot3 = scores2rot3*S2rot3';
figure;surf(XI,YI,cfXtrot3');

% Method Parsimax
[S2rot4, T4] = rotatefactors(S2(:,1:2),'Method','parsimax');
scores2rot4 = cXt*S2rot4;
cfXtrot4 = scores2rot4*S2rot4';
figure;surf(XI,YI,cfXtrot4');

% Here we plot the two spectral components for the three cases:
figure;
subplot(2,1,1);plot(wl,[S2rot1(:,1) S2rot2(:,1) S2rot3(:,1) S2rot4(:,1)]), 
legend('varimax','orthomax','promax','parsimax')  
ylabel('Fluorescence')
title('1st spectral component')
subplot(2,1,2);plot(wl,[S2rot1(:,2) S2rot2(:,2) S2rot3(:,2) S2rot4(:,2)]), legend show
legend('varimax','orthomax','promax','parsimax')  
xlabel('wl')
ylabel('Fluorescence')
title('2nd spectral component')

%% Special topic: Norms and condition number
clear, clc, close all

a = [2 3 5 7]'
b = [5 1 0 3]'
norm(a+b)
norm(a) + norm(b)

A = randi(10,4)
B = randi(10,4)
norm(A+B,2)
norm(A,2) + norm(B,2)

norm(A*B,2)
norm(A,2)*norm(B,2)

% Symmetric matrices
C = [-53   45    53    0;
      45  100    16   -4;
      53   16  -108   44;
       0   -4    44  -46];
[Q,D] = eig(C)
x = [1 2 3 4]'
norm(x)
norm(Q'*x)
norm(D*Q'*x)
norm(Q*D*Q'*x)
Dmax = eye(4)*max(abs(D(:)))
norm(Dmax*Q'*x)
norm(Dmax*Q'*x)/norm(x)
norm(C,2)

% Non-symmetric matrices
A = randi(10,4)
sqrt(max(eig(A'*A)))
norm(A,2)

% All matrices
max(svd(A))

% Inverses and condition number
A = randi(10,4);
rank(A)
inv(A)
A1 = A; A1(:,4) = A(:,3)+ ones(4,1)*1e-6
rank(A1)
inv(A1)
A2 = A; A2(:,4) = A(:,3)+ ones(4,1)*1e-10
rank(A2)
inv(A2)
A3 = A; A3(:,4) = A(:,3)+ ones(4,1)*1e-14
rank(A3)
inv(A3)

norm(inv(A))
sqrt(1/min(eig(A'*A)))

c_1 = cond(A,1)
c_1 = norm(A,1)*norm(inv(A),1)
c_2 = cond(A,2)
c_fro = cond(A,'fro')
c_inf = cond(A,inf)

rc_1 = 1/cond(A,1)
rc_1 = rcond(A)

%% Special topic: Independent component analysis

% CDF and PDF
clear, clc, close all

% Normal distribution with 0.05 sigma
nobs = 1000000;
mu = 10;
sigma = 0.05'
pd = makedist('Normal','mu',mu,'sigma',sigma);
y = random(pd,[nobs,1]);
y_max = max(y);
y_min = min(y);
y_step = 0.01;
y_range = floor(y_min):y_step:ceil(y_max);
y_length = length(y_range); 
cdf = zeros(y_length,1);
for i = 1:y_length
    cdf(i) = sum(y<=y_range(i));
end
cdf = cdf/nobs;
cdf_1 = figure
plot(y_range,cdf,'-r','LineWidth',2)
xlabel('Y range')
ylabel('Cumulative probability (CDF)')
hold on

D = toeplitz([0 1 zeros(1,y_length-2)]);
D = triu(D);
D = D-D';
D = D/(2*y_step);
% openvar D
% pdf = D*cdf';
pdf = gradient(cdf,y_step)
plot(y_range,[0 ; pdf(2:end-1); 0],'-b','LineWidth',2);
pdf_area = trapz(y_range,pdf)
legend({'CDF','PDF as derivative'},'Location','Best')

% Here we add the histogram in which the 'area' of each rectangle is the
% count of y's per y_step. By similarity with physical density
% (mass/volume), pdf is the probability density (change in counts/change in
% range).
hist_step = sigma/2;
edges = [y_min y_min+hist_step:hist_step:y_max-hist_step y_max];
histogram(y,edges,'Normalization','pdf','FaceColor','green','FaceAlpha',0.1)
legend({'CDF','PDF as derivative','PDF as histogram'},'Location','Best')

% Normal distribution with 0.5 sigma
nobs = 1000000;
mu = 10;
sigma = 0.5'
pd = makedist('Normal','mu',mu,'sigma',sigma);
y = random(pd,[nobs,1]);
y_max = max(y);
y_min = min(y);
y_step = 0.01;
y_range = floor(y_min):y_step:ceil(y_max);
y_length = length(y_range); 
cdf = zeros(y_length,1);
for i = 1:y_length
    cdf(i) = sum(y<=y_range(i));
end
cdf = cdf/nobs;
cdf_2 = figure
plot(y_range,cdf,'-r','LineWidth',2)
xlabel('Y range')
ylabel('Cumulative probability (CDF)')
hold on

D = toeplitz([0 1 zeros(1,y_length-2)]);
D = triu(D);
D = D-D';
D = D/(2*y_step);
% openvar D
% pdf = D*cdf';
pdf = gradient(cdf,y_step)
plot(y_range,[0 ; pdf(2:end-1); 0],'-b','LineWidth',2);
pdf_area = trapz(y_range,pdf)
legend({'CDF','PDF as derivative'},'Location','Best')

% Here we add the histogram in which the 'area' of each rectangle is the
% count of y's per y_step. By similarity with physical density
% (mass/volume), pdf is the probability density (change in counts/change in
% range).
hist_step = sigma/2;
edges = [y_min y_min+hist_step:hist_step:y_max-hist_step y_max];
histogram(y,edges,'Normalization','pdf','FaceColor','green','FaceAlpha',0.1)
legend({'CDF','PDF as derivative','PDF as histogram'},'Location','Best')

%% Independence versus uncorrelatedness
close all
y1 = random('normal',0, 10,[1,10000]); 
y2 = random('normal',15, 30,[1,10000]);
y1_bins = internal.stats.histbins(y1);
y2_bins = internal.stats.histbins(y2);
y1_nbins = length(y1_bins);
y2_nbins = length(y2_bins);
scatterhist(y1,y2,'nbins',[y1_nbins,y2_nbins])
y12 = [y1;y2];
y1_pdf = hist(y1,y1_nbins)
y2_pdf = hist(y2,y2_nbins)
y12_pdf_joint = hist3(y12',[y1_nbins,y2_nbins]);
y12_pdf_prod = y1_pdf'*y2_pdf; % outer product

% Normalize the densities so that the integral is 1
y12_pdf_joint = y12_pdf_joint/trapz(trapz(y12_pdf_joint),2);
y12_pdf_prod = y12_pdf_prod/trapz(trapz(y12_pdf_prod),2);
Joint_pdf_1 = figure; 
imagesc(y12_pdf_joint)
set(gca,'Ydir','normal')
Prod_pdf_1 = figure;
imagesc(y12_pdf_prod)
set(gca,'Ydir','normal')

% Theoretical pdf of two independent gaussian variables 
mu = [0 15];
Sigma = cov(y12');
[Y1,Y2] = meshgrid(y1_bins,y2_bins);
F = mvnpdf([Y1(:) Y2(:)],mu,Sigma);
F = reshape(F,y2_nbins,y1_nbins);
Prod_pdf_surf = figure;
h = surf(y1_bins,y2_bins,F);
xlabel('x1'); ylabel('x2'); zlabel('Probability Density');
Prod_pdf_calc = figure;
imagesc(get(h,'ZData')')
set(gca,'Ydir','normal')

% Covariance of two independent gaussian variables
nobs = size(y12,2)
Ey1y2 = y12(1,:)*y12(2,:)'/nobs
Ey1Ey2 = sum(y12(1,:),2)/nobs*sum(y12(2,:),2)/nobs
emp_cov = Ey1y2 - Ey1Ey2
cov(y12',1)

%% uncorrelated but dependent
% Binomial distribution number generation test
% for i = 1:10000
%     r(i) = binornd(10,0.5);
% end
% hist(r)
nobs = 10000;
close all
y1 = random('normal',0, 1,[1,nobs]); 
chsign = binornd(1,0.5,1,nobs);
chsign(chsign == 0) = -1;
sum(chsign)
y2 = chsign.*y1;
y12 = [y1;y2];

Dep_distr = figure;
scatterhist(y1,y2)
y1_nbins = length(internal.stats.histbins(y1));
y2_nbins = length(internal.stats.histbins(y2));

% Empirical joint density (middle)
y12_pdf_joint = hist3(y12',[y1_nbins,y2_nbins])
% Normalize the density so that the P integral is 1
y12_pdf_joint = y12_pdf_joint/trapz(trapz(y12_pdf_joint),2);
Joint_pdf_1 = figure; imagesc(y12_pdf_joint)
set(gca,'Ydir','normal')

% Calculated joint density from product of empirical marginal densities (bottom)
y1_pdf = hist(y1,y1_nbins)
y2_pdf = hist(y2,y2_nbins)
y12_pdf_prod = y1_pdf'*y2_pdf
y12_pdf_prod = y12_pdf_prod/trapz(trapz(y12_pdf_prod),2);
Prod_pdf_1 = figure; imagesc(y12_pdf_prod)
set(gca,'Ydir','normal')

% Covariance
emp_cov = cov(y12',1)
Ey1y2 = y12(1,:)*y12(2,:)'/nobs
Ey1Ey2 = sum(y12(1,:),2)/nobs*sum(y12(2,:),2)/nobs
emp_cov = Ey1y2 - Ey1Ey2

%%
close all
clear, clc
a1 = [1;0.5]; a2 = [0;1.0]; 
a1 = a1/norm(a1); a2 = a2/norm(a2);
cos_theta = a1'*a2;acos(cos_theta)*180/pi
s1 = random('uniform',-8 , 8,[1,1000]); 
s2 = random('uniform',-8 , 8,[1,1000]); 

% Here we produce linear combinations of the original unit vectors:
X = [a1 a2]*[s1;s2];
S = [s1;s2];
ICA_1 = figure;
set(ICA_1,'Unit','Normalized','Position',[0 0.2 0.6 0.8])
plot(X(1,:),X(2,:),'.','Linewidth',0.5,'MarkerEdgeColor','b',...
'MarkerSize',15,'MarkerFaceColor','g')
line(1*[0 3*a1(1)],1*[0 3*a1(2)],'LineWidth',4,'Color','red');
line(1*[0 3*a2(1)],1*[0 3*a2(2)],'LineWidth',4,'Color','green');
xlabel('X1');ylabel('X2');grid on; 
box on;axis equal;hold on

% PCA of the ensemble
% covX = cov(X');
% [Ec,Dc] = eig(covX)
% or simply:
[E,~,D] = pca(X'); 
line(2*[0 E(1,1)],2*[0 E(2,1)],'LineWidth',3,'Color','yellow');
line(2*[0 E(1,2)],2*[0 E(2,2)],'LineWidth',3,'Color','magenta');

% Prior to running the following lines of code install the Toolboxes:
% FastICA (http://research.ics.aalto.fi/ica/fastica/) and RADICAL
% (http://people.cs.umass.edu/~elm/ICA/) in the TOOLBOXES directory, and
% add them to your path using 'addpath', i.e.:
% addpath(genpath('../TOOLBOXES/FastICA'))
% addpath(genpath('../TOOLBOXES/RADICAL'))

% FASTICA
[S1, A1, W1] = fastica(X)
% RADICAL:
[S3,W3] = RADICAL(X)
A3 = inv(W3)

ICA_2 = figure;
set(ICA_2,'Unit','Normalized','Position',[0 0.2 0.6 0.8])
plot(X(1,:),X(2,:),'.','Linewidth',0.5,'MarkerEdgeColor','b',...
'MarkerSize',15,'MarkerFaceColor','g')
xlabel('X1');ylabel('X2');grid on; 
box on;axis equal;hold on

line([0 10*E(1,1)],[0 10*E(2,1)],'LineStyle','--','LineWidth',4,'Color',[1 .5 0]);
line([0 10*E(1,2)],[0 10*E(2,2)],'LineStyle','--','LineWidth',4,'Color',[1 .5 0]);
line([0 2*A1(1,1)],[0 2*A1(2,1)],'LineStyle','-','LineWidth',4,'Color','green');
line([0 2*A1(1,2)],[0 2*A1(2,2)],'LineStyle','-','LineWidth',4,'Color','green');
line([0 2*A3(1,1)],[0 2*A3(2,1)],'LineStyle','-','LineWidth',4,'Color','cyan');
line([0 2*A3(1,2)],[0 2*A3(2,2)],'LineStyle','-','LineWidth',4,'Color','cyan');
line([0 2*a1(1)],[0 2*a1(2)],'LineWidth',4,'Color','red');
line([0 2*a2(1)],[0 2*a2(2)],'LineWidth',4,'Color','red');

ICA_2B = figure;
set(ICA_2B,'Unit','Normalized','Position',[0 0.2 0.6 0.8])
plot(S3(1,:),S3(2,:),'.','Linewidth',0.5,'MarkerEdgeColor','b',...
'MarkerSize',15,'MarkerFaceColor','g')
xlabel('X1');ylabel('X2');grid on; 
box on;axis equal;hold on
cov(S3')

%% Normal Distributions

close all
clear, clc
a1 = [1;0.5]; a2 = [0;1.0]; 
a1 = a1/norm(a1); a2 = a2/norm(a2);
cos_theta = a1'*a2;acos(cos_theta)*180/pi
s1 = random('normal',zeros(1,1000),1.0);
s2 = random('normal',zeros(1,1000),1.0);
cov(s1',s2'),corr(s1',s2')

% Here we produce linear combinations of the original unit vectors:
X = [a1 a2]*[s1;s2];
ICA_3 = figure;
set(ICA_3,'Unit','Normalized','Position',[0 0.2 0.6 0.8])
plot(X(1,:),X(2,:),'.','Linewidth',0.5,'MarkerEdgeColor','b',...
'MarkerSize',15,'MarkerFaceColor','g')
line(1*[0 3*a1(1)],1*[0 3*a1(2)],'LineWidth',4,'Color','red');
line(1*[0 3*a2(1)],1*[0 3*a2(2)],'LineWidth',4,'Color','green');
xlabel('X1');ylabel('X2');grid on; 
box on;axis equal;hold on

% PCA of the ensemble
% covX = cov(X');
% [Ec,Dc] = eig(covX)
% or simply:
[E,~,D] = pca(X'); 
line(3*[0 E(1,1)],3*[0 E(2,1)],'lineStyle','--','LineWidth',4,'Color',[1 .5 0]);
line(3*[0 E(1,2)],3*[0 E(2,2)],'lineStyle','--','LineWidth',4,'Color',[1 .5 0]);

% FASTICA:
[S1, A1, W1] = fastica(X);
% RADICAL:
[S3,W3] = RADICAL(X);
A3 = inv(W3);

ICA_4 = figure;
set(ICA_4,'Unit','Normalized','Position',[0 0.2 0.6 0.8])
plot(X(1,:),X(2,:),'.','Linewidth',0.5,'MarkerEdgeColor','b',...
'MarkerSize',15,'MarkerFaceColor','g')
xlabel('X1');ylabel('X2');grid on; 
box on;axis equal;hold on

line([0 3*E(1,1)],[0 3*E(2,1)],'LineStyle','--','LineWidth',4,'Color',[1 .5 0]);
line([0 3*E(1,2)],[0 3*E(2,2)],'LineStyle','--','LineWidth',4,'Color',[1 .5 0]);
line([0 2*A1(1,1)],[0 2*A1(2,1)],'LineStyle','-','LineWidth',4,'Color','green');
line([0 2*A1(1,2)],[0 2*A1(2,2)],'LineStyle','-','LineWidth',4,'Color','green');
line([0 2*A3(1,1)],[0 2*A3(2,1)],'LineStyle','-','LineWidth',4,'Color','cyan');
line([0 2*A3(1,2)],[0 2*A3(2,2)],'LineStyle','-','LineWidth',4,'Color','cyan');
line([0 2*a1(1)],[0 2*a1(2)],'LineWidth',4,'Color','red');
line([0 2*a2(1)],[0 2*a2(2)],'LineWidth',4,'Color','red');

%% An example of ICA

addpath(genpath('GENERAL_SCRIPTS_FUNCTIONS'));
addpath(genpath('TOOLBOXES/fastMI'));
addpath(genpath('TOOLBOXES/KernelMI'));

close all
pd1 = makedist('uniform','lower',-3,'upper',3);
s1 = random(pd1,1,1000);
pd2 = makedist('uniform','lower',-2.5,'upper',2.5);
s2 = random(pd2,1,1000);
S = [s1;s2];
A = [2 3;2 1];
X = A*S;

ICA_5 = figure;
set(ICA_5,'Unit','Normalized','Position',[0 0.2 0.6 0.8])
plot(X(1,:),X(2,:),'.','Linewidth',0.5,'MarkerEdgeColor','b',...
    'MarkerSize',15,'MarkerFaceColor','g')
xlabel('X1');ylabel('X2');
xlim([-15 15]);ylim([-12 12]);
ax = gca;
ax.XAxisLocation = 'origin';
ax.YAxisLocation = 'origin';
box on, grid on
axis equal
hold on
line([0 4*A(1,1)],[0 4*A(2,1)],'LineStyle','-','LineWidth',4,'Color',[1 .0 0]);
line([0 4*A(1,2)],[0 4*A(2,2)],'LineStyle','-','LineWidth',4,'Color',[1 .0 0]);

ICA_5b = figure;
set(ICA_5b,'Unit','Normalized','Position',[0 0.2 0.6 0.8])
plot(S(1,:),S(2,:),'.','Linewidth',0.5,'MarkerEdgeColor','b',...
    'MarkerSize',15,'MarkerFaceColor','g')
xlabel('X1');ylabel('X2');
ax = gca;
ax.XAxisLocation = 'origin';
ax.YAxisLocation = 'origin';
box on, grid on
axis equal

% Centering
mX = mean(X,2);
cX = X - mX(:,ones(1,1000));

% Whitening: we also set the whitened data as global variable so we can
% pass it to the minimizer.
[E,~,D] = pca(cX');
white = (E*diag(D.^-0.5)*E')
global X_white 
X_white = white*cX;
cov(X_white')

ICA_6 = figure;
set(ICA_6,'Unit','Normalized','Position',[0 0.2 0.6 0.8])
plot(X_white(1,:),X_white(2,:),'.','Linewidth',0.5,'MarkerEdgeColor','b',...
    'MarkerSize',15,'MarkerFaceColor','g')
xlabel('X1');ylabel('X2');
xlim([-2.5 2.5]);ylim([-2.5 2.5]);
ax = gca;
ax.XAxisLocation = 'origin';
ax.YAxisLocation = 'origin';
box on, grid on
axis equal

% Here we look at the distributions
histogram(X_white(1,:))
histogram(X_white(2,:))

% We check the MI before any rotation is applied (angle = 0 degrees) using
% the function mi_from_angle. Notice that if the method used is either
% 'ksdensity' or 'ctransform', the MI value reported is actually a value of
% rho correlation (dependency) for the gaussian copulafit. The recommended
% method is 'ctransform'.

global method
method = 'ctransform';

MI = mi_from_angle(0,X_white)
 
% Here we test the effects on MI of a rotation by a small angle.
phi = -8;
Q = [cosd(phi) sind(phi);-sind(phi) cosd(phi)];
X_q = (Q*X_white);

ICA_7 = figure;
set(ICA_7,'Unit','Normalized','Position',[0 0.2 0.6 0.8])
plot(X_q(1,:),X_q(2,:),'.','Linewidth',0.5,'MarkerEdgeColor','b',...
    'MarkerSize',15,'MarkerFaceColor','g')
xlabel('X1');ylabel('X2');
xlim([-2.5 2.5]);ylim([-2.5 2.5]);
ax = gca;
ax.XAxisLocation = 'origin';
ax.YAxisLocation = 'origin';
box on, grid on
axis equal
 
MI = mi_from_angle(phi,X_white)
 
% Here we set up a multistart minimization to avoid remaining trapped in
% local minima: for this purpose we use MATLAB fminsearch function, which
% calls mi_from_angle at each step of the minimizer.
 
start_vec = [-3:-3:-45];
nstart = length(start_vec);
phi = zeros(1,nstart);
fval = zeros(1,nstart);
options = optimset('TolFun',1e-8, 'TolX',1e-8);
 
for i = 1:nstart
    start_point = start_vec(i);
    [phi(i),fval(i),~,~] = fminsearch(@mi_from_angle,start_point,options);
end
 
% Here we plot the multistart search and identify the rotation angle
% associated with the smallest value of MI.
ICA_8 = figure;
set(ICA_8,'Unit','Normalized','Position',[0 0.2 0.6 0.8]);
scatter(start_vec,fval,100,'r','LineWidth', 2.5,'MarkerFaceColor','y');
box on, grid on
xlabel('Start angle');ylabel('Function value')
[~,min_ind] = min(fval);
best_phi = phi(min_ind);
% best_phi = mod(best_phi,-45);

% Best rotation 
Q = [cosd(best_phi) sind(best_phi);-sind(best_phi) cosd(best_phi)];

% Here we combine whitening and rotation to obtain the separating matrix W.
W = Q*white

% Independent components and empirical mixing matrix from ICA.
Z = W*cX;
A_ica = inv(W)
S_ica = Z + W*mX(:,ones(1,1000));

% FASTICA:
[S_fastica, A_fastica, W_fastica] = fastica(X);

ICA_9 = figure;
set(ICA_9,'Unit','Normalized','Position',[0 0.2 0.6 0.8])
plot(S_ica(1,:),S_ica(2,:),'.','Linewidth',0.5,'MarkerEdgeColor','b',...
    'MarkerSize',15,'MarkerFaceColor','g')
xlabel('S1');ylabel('S2');
xlim([-2. 2.]);ylim([-2. 2.]);
ax = gca;
ax.XAxisLocation = 'origin';
ax.YAxisLocation = 'origin';
box on, grid on
axis equal

ICA_10 = figure;
set(ICA_10,'Unit','Normalized','Position',[0 0.2 0.6 0.8])
plot(S_fastica(1,:),S_fastica(2,:),'.','Linewidth',0.5,'MarkerEdgeColor','b',...
    'MarkerSize',15,'MarkerFaceColor','g')
xlabel('S1');ylabel('S2');
xlim([-2. 2.]);ylim([-2. 2.]);
ax = gca;
ax.XAxisLocation = 'origin';
ax.YAxisLocation = 'origin';
box on, grid on
axis equal

% Here we check the MI between independent components
MI_ica = kernelmi(S_ica(1,:),S_ica(2,:))
MI_fastica = kernelmi(S_fastica(1,:),S_fastica(2,:))

rho_ica = mi_from_angle(0,S_ica)
rho_fastica = mi_from_angle(0,S_fastica)

