% cd_z_path.m
%
% An implementation of the general centroid method with sign vector
% search path plotting. Based on centroid.m b
%
% Jason R. Blevins <jrb11@duke.edu>
%
%
% Input:         A           Data matrix
%
% Output:        B           Loading matrix
%                V           Factor matrix
%                Z           Matrix of sign vectors
%                Zpathlist   List of sign vector search paths
%
% $Id: cd_z_path.m,v 1.2 2004/07/27 23:52:02 jrblevin Exp $

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

num_iteration = 100;

[n,m] = size(A);

Aold = A;
B = [];
V = [];
Cv = [];
Z = [];
Zpathlist = list();

p = rank(A);

for i = 1:min(p,num_iteration)
%   disp(['***** ITERATION ', num2str(i), '*****']);
   [A_next,z,b,v,npmu,Zpath] = new_decom(A);
   B = [B,b];
   Z = [Z,z];
   V = [V,v];
   Cv = [Cv,npmu];
   A = A_next;
   Zpathlist = append(Zpathlist,Zpath);
end

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

% Function:      new_decom
% Description:   Perform one iteration of the general centroid
%                algorithm
%
% Input:         A        data matrix
%
% Output:        A_next   the rank-reduced data matrix
%                v        centroid factors
%                b        factor loadings
%                npmu     n*(centroid values)
%                Zpath    Matrix containing all choices of z

function [A_next,z,b,v,npmu,Zpath] = new_decom(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

% Choose the starting sign vector
z = ones(n,1);

% Initialize a matrix to store the sign vector path and store
% the starting sign vector
Zpath = z;

k = 1;
w = A'*z;

j = 1;

while ~isempty(k)
   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

      % Store the second choice of the sign vector:
      Zpath=[Zpath,z];

      if z(k) == 1					%Step 4
         w = w + 2*A(k,:)';
      else
         w = w - 2*A(k,:)';
      end
      
      j = j + 1;
      
   else
      k = [];
   end
end

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

% Display the Zpath for this iteration
%Zpath