% 
% Written by:
% -- 
% John L. Weatherwax                2005-08-14
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

close all; drawnow; 
clc; 
clear; addpath( '../../Code/eda_data' ); addpath( '../../Code/eda_toolbox' ); 

load yeast; 

% group everthing togther into one matrix: 
X = data; 

[n,p] = size(X);
% For m = 5 random starts, find the N = 5 best projection planes.
%
N = 5;	  % <- specify MORE projection planes ... 

m = 5;	
% Set values for other constants.
c = tan(80*pi/180);
half = 30;
% These will store the results for the N structures.
astar = zeros(p,N);   
bstar = zeros(p,N);
ppmax = zeros(1,N);

% Sphere the data.
muhat = mean(X);
[V,D] = eig(cov(X));
Xc = X-ones(n,1)*muhat;
Z = ((D)^(-1/2)*V'*Xc')';

% Now do the PPEDA: Find a structure, remove it, and look for another one.
Zt = Z;
for i=1:N
    % Find one structure using the chi-squared index:
    [astar(:,i),bstar(:,i),ppmax(i)] = ppeda(Zt,c,half,m,'chi');
    % Find a second structure using the "moment" index: 
    %[astar(:,i),bstar(:,i),ppmax(i)] = ppeda(Zt,c,half,m,'mom');
    % Now remove the structure.
    Zt = csppstrtrem(Zt,astar(:,i),bstar(:,i));
end

% Now project and see the structure: 
% 
for i=1:N
    proj1 = [astar(:,i), bstar(:,i)];
    Zp1 = Z*proj1;
    figure
    plot(Zp1(:,1),Zp1(:,2),'k.'),title(['Structure ',num2str(i)])
    xlabel('\alpha*'),ylabel('\beta*')
    if( (i <= 2) ) saveas( gcf, ['../../WriteUp/Graphics/Chapter4/prob_4_7_yeast_ppeda_si_',num2str(i)], 'epsc' ); end 
end