% centroidsign.m
%
% $Id: centroidsign.m,v 1.5 2004/07/25 16:29:25 jrblevin Exp $
%
% Centroid method modified to compute cost and time and begin with a
% given set of sign vectors.
%
% Based on code by Moody T. Chu <chu@math.ncsu.edu>
% Modified by Jason R. Blevins <jrb11@duke.edu>
%
% Usage:
%
%   [B,V,Z,count,time,Count,Time] = centroidsign(A,Z,k)
% 
% Parameters:
%
%   A      - Matrix m-by-l                             (Required)
%   Z0     - Initial sign vectors (default=ones(m,l))  (Optional)
%   k      - for a rank-k truncation (default=100)     (Optional)
%
% Returns:
%
%   B,V,Z  - 
%   count  - 
%   time   - 
%   Count  - 
%   Time   - 

function [B,V,Z,count,time,Count,Time] = centroidsign(A,varargin);

% Initialize variables
Aold = A;
B = [];
V = [];
Cv = [];
Z = [];
% Count (Uppercase) is a list of the number of sign vector searches
% performed at each successive step in the rank reduction of A.
Count = [];
Time = [];

[n,m] = size(A);
p = rank(A);

% Parse optional parameter list varargin
args = length(varargin);

if (args > 0)
   Z0 = varargin{1};
else
   Z0 = ones(size(A));
end;

if (args > 1)
  num_iteration = varargin{2};
else
  num_iteration = p;
end;

% This is the outermost loop, performed p = rank(A) times.  The i-th
% sign vector is passed to the new_decom() function as a starting
% point for the sign vector search.
for i = 1:min(p,num_iteration)
  [A_next,z,b,v,npmu,count,time] = new_decom(A,Z0(:,i));
  B = [B,b];
  Z = [Z,z];
  V = [V,v];
  Cv = [Cv,npmu];
  A = A_next;
  Count = [Count, count];
  Time = [Time, time];
end

% Sum the individual iteration counts to get an overall total
[p,q]=size(Count);
count = Count * ones(q,1);
time = Time * ones(q,1);

% 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,count,time] = new_decom(A,z0);

% Finds the most significant centroid factor in A according to
% Algorithm 8.1.
%
% Parameters:
%
%  A       - indexing matrix (n by l)
%  z0      - initial sign vector
%
% Returns:
%
%  z       - resulting sign vector
%  b       - loading vector
%  v       - centroid factor
%  npmu    - n*(centroid value)
%  count   - number of sign vector searches performed
%  time    - time elapsed during computation

% Initialize the sign vector search counter to 0
count = 0;
time = clock;

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

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

% This is the requested starting sign vector
z = z0;
k = 1;
w = A'*z;

while ~isempty(k)

  % Increment the sign vector search counter
  count = count + 1;

  % Step 1
  e = d - sign(z).*(A*w);

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

npmu = w'*w;

% Centroid factors
v = w/sqrt(npmu);

% Loadings
b = A*v;

% Rank Reduction
A_next = A - b*v';

time = etime(clock,time);