% 
% Written by:
% -- 
% John L. Weatherwax                2005-08-04
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
% Epage 64 - 65
%
%-----

close all; clc; clear; 

rand('seed',0);
randn('seed',0);

X = [ 1, 1 ; 2, 2; 2, 1; 1, -1; 1, -2; 2, -2; -1, 1; -1, 2; -2, 1 ];
labels = [ 1, 1, 1, 2, 2, 2, 3, 3, 3 ]; 

inds = find( labels==1 );
plot( X(inds,1), X(inds,2), '.g' ); axis( [-3,+3,-3,+3] ); hold on; 
inds = find( labels==2 );
plot( X(inds,1), X(inds,2), '.k' );
inds = find( labels==3 );
plot( X(inds,1), X(inds,2), '.b' );

[N,p] = size(X);
M     = length(unique(labels));

% prepend a column of ones to the data: 
%
X      = [ X, ones(size(X,1),1) ]; 
X_orig = X;

X_Kesler = kesler_construction(X,labels);
X        = X_Kesler; 

% Train a perceptron so that w^T x > 0 for all samples:
%
w_i = rand(M*(p+1),1);

rho = 0.5; 
maxIters = 1000;

Niters = 0; 
while( 1 ) % we assume the data is linearly seperable

  % Find the set J (the missclassified samples) with this weight vector:
  predicted_class = X * w_i;
  Y = find( predicted_class < 0 ); % find the indices of misclassified vectors 
  if( isempty(Y) ) % we are done ... everything is classified correctly!
    break;
  end
  
  delta_w = sum( -X( Y, : ), 1 ).';

  fprintf('Niter= %10d; updating w_i; size= %10.6f\n',Niters,norm(delta_w));
  
  w_i = w_i - rho * delta_w;
  
  Niters = Niters + 1; 
  if( Niters > maxIters ) 
    fprintf('max number of iterations= %10d exceeded\n',maxIters);
    break
  end
end

% draw the decision boundary w^T x = 0 for each sample: 
%
x1_grid = linspace( min(X(:,1)), max(X(:,1)), 50 );

w1      = w_i( 1:3, 1 );
w1_db   = ( -w1(3) - w1(1) * x1_grid ) / w1(2);
plot( x1_grid, w1_db, '-g' ); hold on; 

w2      = w_i( 4:6, 1 );
w2_db   = ( -w2(3) - w2(1) * x1_grid ) / w2(2);
plot( x1_grid, w2_db, '-k' ); 

w3      = w_i( 7:9, 1 );
w3_db   = ( -w3(3) - w3(1) * x1_grid ) / w3(2);
plot( x1_grid, w3_db, '-b' ); 

title('Kesler construction computed decision line');

fn = ['chap_3_example_2_class_regions.eps'];
saveas(gcf,['../../WriteUp/Graphics/Chapter3/',fn],'epsc');

% Lets verify that we have a valid multidensional solution 
% i.e. one that our discriminant functions correctly classify each point 
%
for ii=1:M
  inds = find( labels==ii );
  X_ii = X_orig( inds, : );
  N_ii = length(inds);
  for si=1:N_ii
    fprintf('class= %3d; w1^T x= %8.4f; w2^T x= %8.4f; w3^T x= %8.4f\n', ii, dot( w_i(1:3,1), X_ii(si,:) ), dot( w_i(4:6,1), X_ii(si,:) ), dot( w_i(7:9,1), X_ii(si,:) ) );
  end
end