% 
% 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 posse; 

%X = croix;    % <- a cross
%X = struct2;  % <- a L shape
%X = boite;    % <- a donut 
%X = groupe;   % <- four clusters  
%X = curve;    % <- a curve
X = spiral;   % <- a sprial

[n,p] = size(X);
% For m = 5 random starts, find the N = 5 best projection planes.
%
N = 2;	  % <- specify the number of 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.'),
    %fnp='croix';        % <- very spherical clusters ... 
    %fnp='struct2';
    %fnp='boite';         % <- initial small cluster ... other projections are much more spherical
    %fnp='groupe';        % <- initial two clusters in the first and second projection
    %fnp='curve';         % <- initial cluster contains a "hole"
    fnp='spiral';        % <- very spherical clusters ... sometimes with and initial spiral
    title([fnp,': Structure ',num2str(i)]); xlabel('\alpha*'),ylabel('\beta*')
    if( 1 & (i <= 2) ) saveas( gcf, ['../../WriteUp/Graphics/Chapter4/prob_4_7_posse_',fnp,'_ppeda_si_',num2str(i)], 'epsc' ); end 
end