function [ x_ret, y_ret ] = hits( A , mode )
%HITS   Calculates auth.- and hub-scores with HITS
%
% [ x_ret, y_ret ] = HITS( A , mode ) where A is the adjacency 
%     matrix of the network.  mode  results in the output to be
%   - 's' : directly the scores (unsorted)
%   - 'so': the scores sorted in decr. order
%   - 'id': Docids (sorted in decr. order of importance) (default)
%     where x stands for auth.- and y for hub-scores.
%
%  Last change:  10/2/2005 - Florian Knorn

if nargin < 1 || nargin > 2
    error('Please input adjacency matrix (and mode) !');
else
	n = size(A,1);
	if n ~= size(A,2)
		error('Please input SQUARE matrix !');
	end
	if nargin == 1
		mode = 'id';
	end
end

eps = 10^(-floor(log10(n)+1)-3); % error tolerance
err = 42;				% initial error
x(:,1) = ones(n,1)/n;	% initial auth. scores

A = double(A);			% make sure A isn't a logical
ATA = A'*A;				% create authority matrix

% power method iterations
k = 1;
while err > eps
    k = k + 1;			% increase counter
    x(:,k) = ATA * x(:,k-1); % update x
    x(:,k) = x(:,k) / sum( x(:,k) ); % normalize x
    err = abs(  sum( x(:,k)-x(:,k-1) )  ); % error
end

% calculate y algebraically
x = x(:,k); y = A*x; y = y / sum(y);

% return results
if strcmp(mode,'s') && ~strcmp(mode,'os') && ~strcmp(mode,'so')
	x_ret = x; y_ret = y;
elseif strcmp(mode,'so') || strcmp(mode,'os')
	x_ret = sort(x, 'descend');
	y_ret = sort(y, 'descend');
elseif strcmp(mode,'id')
	[xtemp,xi] = sort(x, 'descend');
	[ytemp,yi] = sort(y, 'descend');
	x_ret = xi; y_ret = yi;
else
	error('unknown mode !');
end

% % sort the results
% [ x_sorted , x_docids ] = sort(x, 'descend');
% [ y_sorted , y_docids ] = sort(y, 'descend');
% 
% % display results
% disp('place      authscore        docid        |         hubscore        docid');
% for k = 1:n
%     disp(sprintf('%3i.        %0.6f         %3i         |         %0.6f         %3i', ...
%         k , x_sorted(k) , x_docids(k), y_sorted(k), y_docids(k)));
% end