%
% Problem Epage 202
% Example Epage 190
% 
% 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 bank;   % <- two matrices forge/genuine

% put the data together (stacked so that they appear as ONE data set): 
X = [genuine;forge];
% plot the raw data: 
%figure; imagesc(X); colorbar

% 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')';

% 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; 
  plot( Zp1(1:100,1), Zp1(1:100,2), 'k.' ); hold on; 
  plot( Zp1(101:end,1), Zp1(101:end,2), 'rx' );
  title( sprintf('PPEDA bank projection=%d',si) );
  xlabel( '\alpha^*' ); ylabel( '\beta^*' ); 
  saveas( gcf, sprintf('ppeda_bank_proj_%d',si), 'epsc' ); 
end