%
% 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 posse;   % <- 6 matrices matrices boite/croix/curve/groupe/spiral/struct2

% put the data together (stacked so that they appear as ONE data set): 
X = [boite;croix;curve;groupe;spiral;struct2];

%--
% for plotting the "truth" we define some variables ...
%
% the number of known classes:
nClass = 6; 
% the number of samples per class: 
nSamps = size(boite,1); 
% the dimension of the feature space: 
d = size(X,2); 

% plot the raw data: 
if( 0 )
  figure; imagesc(X); colorbar
  % plot boxplots to determine how sepearable this data is:
  for fi=1:d,
    figure;
    boxplot( [boite(:,fi),croix(:,fi),curve(:,fi),groupe(:,fi),spiral(:,fi),struct2(:,fi)] );
    set(gca,'XTickLabel',{'boite','croix','curve','groupe','spiral','struct2'});
    title( sprintf('feature %d',fi) ); 
  end
  % plot Andrews curves 
  figure; hold on; col=['rbgmkyc'];
  for ci=1:nClass,
    inds=(nSamps*(ci-1)+1):(nSamps*ci);
    rinds=randperm(length(inds));       % <- randomly downsample a few instances 
    inds=inds(rinds(1:50)); 
    csandrews( X(inds,:), '-', col(ci) );
  end
  % plot the classes in parallel coordinates: 
  figure; hold on; col=['rbgmkyc'];
  for ci=1:nClass,
    inds=(nSamps*(ci-1)+1):(nSamps*ci);
    rinds=randperm(length(inds));       % <- randomly downsample a few instances 
    inds=inds(rinds(1:50)); 
    csparallel( X(inds,:), '-', col(ci) );
  end
end


% the number of structure removals: 
N = 5; 
% the number of random starts: 
%m = 10; 
%m = 20; 
m = 50; 

c = tan(80*pi/180); 
% number of steps with no increase: 
%half = 30; 
half = 50; 

% 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 found structure
  Zt = csppstrtrem(Zt,astar(:,i),bstar(:,i)); 
end

% project and see the structure: 
col=['rbgmkyc'];
for si=1:N,
  proj1 = [astar(:,si), bstar(:,si)];
  Zp1 = Z*proj1; 
  % map back in physical dimensions: 
  
  figure; hold on;       
  for ci=1:nClass,
    inds=(nSamps*(ci-1)+1):(nSamps*ci);
    plot( Zp1(inds,1), Zp1(inds,2), [col(ci),'.'] );
  end
  title( sprintf('PPEDA posse projection=%d',si) );
  xlabel( '\alpha^*' ); ylabel( '\beta^*' ); 
  saveas( gcf, sprintf('ppeda_posse_proj_%d',si), 'epsc' ); 
end