%% This program is the code of principal component accumulation (PCAcc) method 

clear;

%% loading data
load Leukemial
%load Lung

%% dividing group
dd_c = data_c(:,2:end);
y_c = data_c(:,1);
n_c = find(diff(y_c)~=0); 

dd_v = data_v(:,2:end);
y_v = data_v(:,1);
n_v = find(diff(y_v)~=0); 

sp.A = dd_c(1:n_c,:);
sp.B = dd_c(n_c+1:end,:);

[grps, X] = analgroup(sp);

%% Classifcation using classical PCA
maxrank = 18; % The parameter of PCA; Leukemial - 18 and Lung - 13

[R_c,R_v] = svdpca(X,dd_v,maxrank);
disprotpcsallpoints(R_c,grps,3,4);
title('Distribution of the training set in PC plot')
xlabel('PC 3'); ylabel('PC 4')

[maxR,pcind,RT,T,P,T_v,maxR_v] = rotation(R_c,grps,maxrank,R_v);

% training
[Rs1 meanR111 meanR211 stdR111 stdR211] = dispsumR(maxR,grps);
title('Distribution of the training set on the projected PC axis')
xlabel('Projected PC axis'); 
[Error1,error11,error21,index1] = disclass(maxR,y_c,meanR111,meanR211);

% prediction
[Error1_v,error11_v,error21_v,index1_v] = disclass(maxR_v,y_v,meanR111,meanR211);
figure;plot(maxR_v(1:n_v), ones(size(maxR_v(1:n_v))), 'o');
hold;plot(maxR_v(n_v+1:end), ones(size(maxR_v(n_v+1:end))),'go');
title('Distribution of the test set on the projected PC axis')
xlabel('Projected PC axis'); 


%% Classification using AccPc method
% selecting PC subspace
[maxRs PP W W_v] = stats(RT,T,T_v,P,grps);

% accumulating the accepted PC subspaces
S = []; sumR = zeros(size(W(:,1))); sumR_v = zeros(size(W_v(:,1)));%

for i=1:size(W,2)  
   sumR = sumR+W(:,i);
   sR = (sumR-min(sumR))./(max(sumR)-min(sumR));
   [rt2 meanR11 meanR21 stdR11 stdR21] = dispsumRs(sR,grps);
   S = [S;rt2]; 
   sumR_v = sumR_v+W_v(:,i);
   sR_v = (sumR_v-min(sumR_v))./(max(sumR_v)-min(sumR_v));
end
figure;plot(S,'o-');
title('Variation of Rs with the number of the accumulated PC subspaces')
xlabel('Number of PC subspaces'); ylabel('Accumulated Rs');

% training
[Rs2 meanR1 meanR2 stdR1 stdR2] = dispsumR(sR,grps);
title('Distribution of the training set on the accumulated PC axis')
xlabel('Accumulated PC axis');
[Error3,error13,error23,index3] = disclass(sR,y_c,meanR1,meanR2);

% prediction
[Error3_v,error13_v,error23_v,index3_v] = disclass(sR_v,y_v,meanR1,meanR2);
figure;plot(sR_v(1:n_v), ones(size(sR_v(1:n_v))), 'o');
hold;plot(sR_v(n_v+1:end), ones(size(sR_v(n_v+1:end))),'go');
title('Distribution of the test set on the accumulated PC axis')
xlabel('Accumulated PC axis');
