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

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

%figure; imagesc( X ); colorbar; 

% maybe preprocessing this matrix so that every vector has zero mean and unit variance: 
% none:             
% range_trans:      seems to give rather "boxy" nonlinear projections
% range_trans1:     
% zscore:           
% whitening:        
pp='none'; 
%addpath( '../Chapter1' ); X = range_trans(X); pp='range_trans'; 
%addpath( '../Chapter1' ); X = range_trans(X,1); pp='range_trans1'  % <- 
%X = zscore(X);                                    % <- roughly spherical clusters 
%addpath( '../../../../Duda_Hart_Stork/BookSupplements/ClassificationToolbox/Src' ); X = Whitening_transform(X.').'; 

[n,p] = size(X);
% For m = 5 random starts, find the N = 5 best projection planes.
%
N = 2;	  % <- 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*')
    fnp=['spam','_',pp];
    if( 1 & (i <= 2) ) saveas( gcf, ['../../WriteUp/Graphics/Chapter4/prob_4_7_',fnp,'_ppeda_si_',num2str(i)], 'epsc' ); end 
end