% centroid_plot.m
%
% Jason Blevins <jrb11@duke.edu>
%
% Modification of centroid.m to produce plots of each iteration.  The
% data matrix must contain two dimensional normalized data since a
% compassplot is used.  In other words, A must have two columns with
% rows of unit length.
%
% $Id: centroid_plot.m,v 1.2 2004/07/28 00:00:52 jrblevin Exp $

function [B,V,Z]=centroid_plot(A);

num_iteration = 100;

[n,m] = size(A);

Aold = A;
B = [];
V = [];
Cv = [];
Z = [];

p = rank(A);

   hData=compass(A(:,1),A(:,2));
   set(hData,'LineWidth',2);
   set(hData,'Color','k');
   hold on;

   % Plot the standard basis vectors
   %I = eye(2);
   %hEye = compass(I(:,1),I(:,2));
   %set(hEye,'LineWidth',3);
   %set(hEye,'Color','k');

for i = 1:min(p,num_iteration)
   
   %disp(['***** ITERATION ', num2str(i), '*****']);
   [A_next,z,b,v,npmu] = new_decom(A);
   B = [B,b];
   Z = [Z,z];
   V = [V,v];
   Cv = [Cv,npmu];
   A = A_next;
   
   hFactor = compass(v(1,:),v(2,:));
   set(hFactor, 'LineWidth', 3);
   set(hFactor, 'Color', 'r');
end

hold off;

%disp('***** END *****');

% Note that Rold = R_next + Z*inv(diag(S))*Z';
%
% Note also that S/n is the approximate singular value

function [A_next,z,b,v,npmu] = new_decom(A);
%
% This code is for Algorithm 8.1
%
% Input:
%
%		A = indexing matrix
%
% Output:
%
%		npmu = n*(centroid values)
%		v = centroid factors
%		b = loadings
%

% Uncomment for debug output
A;

[n,m] = size(A);
tol = eps*max([n,m]);

for i = 1:n
   d(i,1) = norm(A(i,:))^2;		% diagonal of AA'
end

% Uncomment for debug output
d;

% this is the start value 

z = ones(n,1);
k = 1;
w = A'*z;

while ~isempty(k)
   %disp('***** STEP *****');
   e = d - sign(z).*(A*w);			%Step 1   
   [y,k] = max(e);					%Step 2
   
   if y > 0 & abs(y) > tol    
      z(k) = -z(k);					%Step 3
      if z(k) == 1					%Step 4
         w = w + 2*A(k,:)';
      else
         w = w - 2*A(k,:)';
      end
   else
      k = [];
   end
end

%disp('***** DATA *****');

npmu = w'*w;
v = w/sqrt(npmu);					%Centroid factors
b = A*v;							%Loadings
A_next = A - b*v';