function [W, H, nmfV] = NMF(V, parameter)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Name: NMF
% Date: May 2018
% Programmer: Christian Dittmar
%
% Given a non-negative matrix V, find non-negative templates W and activations
% H that approximate V.
%
% References:
% [1] Lee, DD & Seung, HS. "Algorithms for Non-negative Matrix Factorization"
%
% [2] Andrzej Cichocki, Rafal Zdunek, Anh Huy Phan, and Shun-ichi Amari
% "Nonnegative Matrix and Tensor Factorizations: Applications to
% Exploratory Multi-Way Data Analysis and Blind Source Separation"
% John Wiley and Sons, 2009.
%
% Input:  V               K x M non-negative matrix to be facorized%
%         parameter.
%           costFunc      Cost function used for the optimization, currently
%                         supported are:
%                         'EucDdist' for Euclidean Distance
%                         'KLDiv' for Kullback Leibler Divergence
%                         'ISDiv' for Itakura Saito Divergence
%           numIter       Number of iterations the algorithm will run.
%           numComp       the rank of the approximation
%
% Output: W               K x R non-negative templates
%         H               R x M non-negative activations
%         nmfV            cell-array with approximated component matrices
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% If you use the 'NMF toolbox' please refer to:
% [3]  Patricio López-Serrano, Christian Dittmar, Yiğitcan Özer, and Meinard
%     Müller
%     NMF Toolbox: Music Processing Applications of Nonnegative Matrix
%     Factorization
%     In Proceedings of the International Conference on Digital Audio Effects
%     (DAFx), 2019.
%
% License:
% This file is part of 'NMF toolbox'.
% https://www.audiolabs-erlangen.de/resources/MIR/NMFtoolbox/
% 'NMF toolbox' is free software: you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the
% the Free Software Foundation, either version 3 of the License, or (at
% your option) any later version.
%
% 'NMF toolbox' is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
% Public License for more details.
%
% You should have received a copy of the GNU General Public License along
% with 'NMF toolbox'. If not, see http://www.gnu.org/licenses/.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% check parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~isfield(parameter,'costFunc')
  parameter.costFunc = 'KLDiv';
end
if ~isfield(parameter,'numIter')
  parameter.numIter = 30;
end
if ~isfield(parameter, 'fixW')
    parameter.fixW = 0;
end

% get important params
[K, M] = size(V);
R = parameter.numComp;
L = parameter.numIter;

% initialization of W and H
if iscell(parameter.initW)
    W = cell2mat(parameter.initW);
elseif isnumeric(parameter.initW)
    W = parameter.initW;
end;

H = parameter.initH;

% create helper matrix of all ones
onesMatrix = ones(K, M);

% normalize to unit sum
V = V./(eps+sum(V(:)));

% maiin iterations
for iter = 1:L

  % compute approximation
  Lambda = eps+(W*H);

  % switch between pre-defined update rules
  switch parameter.costFunc

    case 'EucDist' % euclidean update rules
        if ~parameter.fixW
            W = W .* (V * H' ./ ((Lambda * H') + eps));
        end
      H = H .* (W' * V ./ ((W' * Lambda) + eps));

    case 'KLDiv' % Kullback Leibler divergence update rules
        if ~parameter.fixW
            W = W .* ((V ./ Lambda) * H') ./ (onesMatrix*H'+eps);
        end
      H = H .* (W' * (V ./ Lambda)) ./ (W'*onesMatrix+eps);

    case 'ISDiv' % Itakura Saito divergence update rules
        if ~parameter.fixW
            W = W .* ((Lambda.^-2 .* V) * H') ./ (Lambda.^-1 * H' + eps);
        end
      H = H .* (W' * (Lambda.^-2 .* V)) ./ (W' * Lambda.^-1 + eps);

    otherwise
      error('Unknown cost function');

  end

  % normalize templates to unit sum
  if ~parameter.fixW
    normVec = sum(W,1);
    W = bsxfun(@times,W,1./(eps+normVec));
  end

end

% compute final output approximation
for r = 1:R
  nmfV{r} = W(:,r)*H(r,:);
end

end