function [ A, B, cost ] = betaNMF( X, beta, K, it, curve, A, B )
%
% Coded by D. Kitamura (d-kitamura@ieee.org) on 31 Jul, 2015 (ver1.0).
% Modified by D. Kitamura (d-kitamura@ieee.org) on 18 Nov, 2018 (ver1.1).
%
% This function decomposes input matrix by nonnegative matrix
% factorization (NMF) based on beta-divergence criterion and
% multiplicative update rules.
% The input nonnegative data matrix X will be decomposed as
% X = A*B approximately, where matrices A and B are often called
% "basis matrix" and "activation matrix," respectively.
% All the entries in X, A, and B must be nonnegative values.
% All the update rules implemented in this function are based on
% majorization-minimization algorithm.
%
% see also
% C. Fevotte, J. Idier, "Algorithms for nonnegative matrix
% factorization with the beta-divergence,"
%
% M. Nakano, H. Kameoka, J. Le Roux, Y. Kitano, N. Ono, and S. Sagayama,
% "Convergence-guaranteed multiplicative algorithms for nonnegative matrix
% factorization with beta-divergence," Proc. MLSP, pp. 283-288, 2010.
%
% See also
% http://d-kitamura.net
%
% [syntax]
%   [ A, B, cost ] = betaNMF( X )
%   [ A, B, cost ] = betaNMF( X, beta )
%   [ A, B, cost ] = betaNMF( X, beta, K )
%   [ A, B, cost ] = betaNMF( X, beta, K, it )
%   [ A, B, cost ] = betaNMF( X, beta, K, it, curve )
%   [ A, B, cost ] = betaNMF( X, beta, K, it, curve, A )
%   [ A, B, cost ] = betaNMF( X, beta, K, it, curve, A, B )
%
% [inputs]
%        X: observed nonnegative matrix (I x J)
%     beta: beta value of beta-divergence (default: 2)
%           (0: Itakura-Saito divergence, 1: Generalized KL divergence, 2: Squared Euclidian distance)
%        K: number of bases (default: ceil(J/10))
%       it: number of iterations (default: 100)
%    curve: draw convergence curve or not (true (1) or false (0), default: false)
%        A: initial value of basis matrix A
%        B: initial value of activation matrix B
%
% [outputs]
%        A: Decomposed basis matrix (I x K)
%        B: Decomposed activation matrix (K x J)
%     cost: Convergence of cost function (it+1 x 1)
%

% Set default values and check errors
[I,J,M] = size(X);
if M ~= 1
    error('This function supports only 2D matrix.\n');
end
if (nargin < 2)
    beta = 2;
end
if (nargin < 3)
    K = ceil(J/10);
end
if (nargin < 4)
    it = 100;
end
if (nargin < 5)
    curve = false;
end
if (nargin < 6)
    A = rand( I, K ); % initialize basis matrix with nonnegative random values
end
if (nargin < 7)
    B = rand( K, J ); % initialize activation matrix with nonnegative random values
end
if ( isreal(X) == 0 || min(min(X)) < 0 )
    error('All the entries of input matrix X must be nonnegative.\n');
end
if ( size(A,1) ~= I || size(A,2) ~= K || size(B,1) ~= K || size(B,2) ~= J )
    error('The size of input initial matrix A or B is wrong.\n');
end
X = max( X, eps ); % machine epsilon flooring for avoiding zero-division in update rules
E = ones( I, J ); % ones matrix whose size is the same as X
if ( beta < 1 ) % exponent value that is used in the update rules
    phi = 1/( 2 - beta ) ;
elseif ( beta < 2 )
    phi = 1 ;
else
    phi = 1/( beta - 1 ) ;
end

% Iterate update rules
fprintf('Iteration:     ');
if ( curve == true )
    cost = zeros( it+1, 1 );
    cost(1) = calcCost_local(X, A*B, beta);
    for i = 1 : it
        [ A, B ] = updateRules_local( X, A, B, E, beta, phi );
        cost(i+1) = calcCost_local(X, A*B, beta);
        fprintf('\b\b\b\b%4d', i);
    end
    figure;
    semilogy( (1:it+1), cost );
    set(gca,'FontName','Times','FontSize',16);
    xlabel('Number of iterations','FontName','Arial','FontSize',16);
    ylabel('Value of cost function','FontName','Arial','FontSize',16);
else
    cost = 0;
    for i = 1 : it
        [ A, B ] = updateRules_local( X, A, B, E, beta, phi );
        fprintf('\b\b\b\b%4d', i);
    end
end
fprintf('\n');
end

%%% Local functions %%%%
function [ A, B ] = updateRules_local( X, A, B, E, beta, phi )
if beta == 2 % update rules based on squared Euclidean distance
    A = A .* (( X*B.' )./( A*(B*B.') ));
    A = max(A,eps);
    B = B .* (( A.'*X )./( A.'*A*B ));
    B = max(B,eps);
elseif beta == 1 % update rules based on generalized KL divergence
    A = A .* (( (X./(A*B))*B.' )./( E*B.' ));
    A = max(A,eps);
    B = B .* (( A.'*(X./(A*B)) )./( A.'*E ));
    B = max(B,eps);
elseif beta == 0 % update rules based on IS divergence
    AB = A*B;
    A = A .* sqrt(( (X.*((AB).^(-2)))*B.' )./( ((AB).^(-1))*B.' ));
    A = max(A,eps);
    AB = A*B;
    B = B .* sqrt(( A.'*(X.*((AB).^(-2))) )./( A.'*((AB).^(-1)) ));
    B = max(B,eps);
else % update rules based on beta-divergence for general value of beta
    AB = A*B;
    A = A .* ((( (X.*((AB).^(beta-2)))*B.' )./( ((AB).^(beta-1))*B.' )).^(phi));
    A = max(A,eps);
    AB = A*B;
    B = B .* ((( A.'*(X.*((AB).^(beta-2))) )./( A.'*((AB).^(beta-1)) )).^(phi));
    B = max(B,eps);
end
end

function [ cost ] = calcCost_local( X, AB, beta )
if beta == 0 % IS divergence
    cost = sum( sum( X./AB - log(X./AB) - 1 ) );
elseif beta == 1 % generalized KL divergence
    cost =  sum( sum( X .*log(X./AB) - X + AB ) );
elseif beta == 2 % squared Euclidean distance
    cost = sum( sum( ( X - AB ).^2 ) );
else % beta-divergence for general value of beta
    cost = sum(sum( abs( ( ( X.^(beta) ) ./...
        ( beta * ( beta-1 ) ) ) + ( ((AB).^(beta)) ./ beta ) - ...
        (X .* (AB).^(beta-1))./(beta-1) ) ) ) ;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% EOF %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%