function [ W, H ] = NNDSVD( A, K, type )
%
% Corded by D. Kitamura (d-kitamura@nii.ac.jp) on 29 Dec. 2015.
%
% This function calculates initial basis and coefficient matrices
% of NMF based on nonnegative double singular value decomposition
% proposed in the following literature:
% C. Boutsidis, E. Gallopoulos, "SVD based initialization: a head
% start for nonnegative matrix factorization," Pattern Recognition,
% vol.41, pp.1350-1362, 2008.
% The process is completely the same as Table 1 in the above paper.
%
% see also
% http://d-kitamura.sakura.ne.jp/index.html
%
% [syntax]
%   [ W, H ] = NNDSVD( A )
%   [ W, H ] = NNDSVD( A, K )
%   [ W, H ] = NNDSVD( A, K, type )
%
% [inputs]
%      A: input nonnegative matrix (M x N)
%      K: number of bases (nonnegative scalar, default: N/10)
%   type: type of algorithm ('o', 'a', or 'r', default: 'o')
%         'o': original NNDSVD
%         'a': NNDSVDa (NNDSVD + replace zeros with an average of A)
%         'r': NNDSVDar (NNDSVD + replace zeros with uniform random values in the range of [0, mean(mean(A))])
%
% [outputs]
%      W: initialized basis matrix (M x K)
%      H: initialized coefficient matrix (K x N)
%

% Check errors and set default values
[M,N] = size(A);
if min(min(A)) < 0
    error('All the entries of input matrix must be nonnegative.');
end
if (nargin < 2)
    K = N/10;
end
if (nargin < 3)
    type = 'o';
end

W = zeros(M, K);
H = zeros(K, N);
[U, S, V] = svds(A, K); % SVD and extracting factors corresponding top k singular values
if min(U(:,1)) < 0 % The first singular vector sometimes has opposite direction (all negative vectors)
    W(:,1) = sqrt(S(1,1))*(-U(:,1));
else
    W(:,1) = sqrt(S(1,1))*U(:,1);
end
if min(V(:,1)) < 0 % The first singular vector sometimes has opposite direction (all negative vectors)
    H(1,:) = sqrt(S(1,1))*(-V(:,1)).';
else
    H(1,:) = sqrt(S(1,1))*V(:,1)';
end
for j = 2:K
    x = U(:,j);
    y = V(:,j);
    [xp, xn] = local_posneg(x);
    [yp, yn] = local_posneg(y);
    xpnorm = norm(xp);
    ypnorm = norm(yp);
    mp = xpnorm*ypnorm;
    xnnorm = norm(xn);
    ynnorm = norm(yn);
    mn = xnnorm*ynnorm;
    if mp > mn
        u = xp/xpnorm;
        v = yp/ypnorm;
        sigma = mp;
    else
        u = xn/xnnorm;
        v = yn/ynnorm;
        sigma = mn;
    end
    W(:,j) = sqrt(S(j,j)*sigma)*u;
    H(j,:) = sqrt(S(j,j)*sigma)*v.';
end
a = mean(mean(A));
if strcmp(type, 'o') %NNDSVD
    W = max(W,eps);
    H = max(H,eps);
elseif strcmp(type, 'a') %NNDSVDa
    W = W + (W==0).*a; % Replace zero elements with an average of A
    H = H + (H==0).*a; % Replace zero elements with an average of A
elseif strcmp(type, 'r') %NNDSVDar
    W = W + (W==0).*rand(M,K).*a/100; % Replace zeros with uniform random values in the range of [0, a]
    H = H + (H==0).*rand(K,N).*a/100; % Replace zeros with uniform random values in the range of [0, a]
end
end

%%% subfunction for getting positive and negative parts %%%%
function [xp,xn] = local_posneg(x)
xp = (x >= 0).*x;
xn = (x < 0).*(-x);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% EOF %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%