%
% Problem Epage 202
% Example Epage 189
% 
% Written by:
% -- 
% John L. Weatherwax                2008-02-20
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

clear all; 
close all; 
clc; 

addpath('../../Code/CSTool');
load ppdata;   % <- one matrix X

% plot the RAW data: 
figure; plotmatrix(X); title( 'the original data' ); drawnow; 

% for m random starts, find the best projection plane using "N" structure removals.
N = 5; 

d = size(X,2); 
% number of random starts: 
m = 4; 
c = tan(80*pi/180); 
% number of steps with no increase: 
half = 30; 

% store the N structures: 
astar = zeros(d,N); 
bstar = zeros(d,N); 
ppmax = zeros(1,N); 

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

% plot the RAW sphered data: 
figure; plotmatrix(Z); title( 'the original data' ); drawnow; 

% now do the PPEDA.  Find a structure, remove it, and look for another one
Zt = Z; 
for i=1:N
  [astar(:,i),bstar(:,i),ppmax(i)] = csppeda(Zt,c,half,m); 
  % now remove the structure
  Zt = csppstrtrem(Zt,astar(:,i),bstar(:,i)); 
end

% project and see the structure: 
for si=1:N,
  proj1 = [astar(:,si), bstar(:,si)];
  Zp1 = Z*proj1; 
  figure; 
  plot( Zp1(:,1), Zp1(:,2), 'k.' ); hold on; 
  title( sprintf('PPEDA ppdata projection=%d',si) );
  xlabel( '\alpha^*' ); ylabel( '\beta^*' ); 
  saveas( gcf, sprintf('ppeda_ppdata_proj_%d',si), 'epsc' ); 
end