function [G,H,Z,zz,PRE,prow,pcol] = NMF(X,D,K,L)
% NMF_nodiag - non-negative matrix factorization 
%              K-means used to cluster the basis vectors
% 
% Input:
%   X  - the matrix to factorize
%   D  - number of basis vectors to generate
%   K  - number of clusters for rows
%   L  - number of clusters for columns
%
% Results:
%   G    - a set of r basis vectors
%   H    - represenations of the columns of X in 
%          the basis given by G
%   Z    - the sum-of-squares loss function value
%   zz   - a vector containing the loss function value for each restart
%   PRE  - percentage reduction of error
%   prow - partition of the row objects
%   pcol - partition of the column objects
%

%---------------------------------------
% Make sure program is usable
%---------------------------------------

if ndims(X) ~= 2
    error('X must be 2 dimensionsal');
end
[N, C] = size(X);
if prod(size(D)) ~= 1 | D < 1
    error('D must be a positive scalar');
end

%---------------------------------------
% Initialization
%---------------------------------------
state = 1;                 % fix state for testing
rand('state', state);
eps = .00000001;
MAXIT = 5000;
restarts = 20;
tic;
zbest = 99999999999;
zz = zeros(restarts,1);
for ijkl = 1:restarts
    
    G = rand(N,D);   % Initialize G
    H = rand(D,C);   % Initialize H
    err = []; delta = 99999999;
    iterations = 1;
    GH = G*H; GH = GH;
    err = [err sum(sum((X-GH).^2))];
%---------------------------------------
% Descent Algorithm of Lee & Seung (2001)
%---------------------------------------

    while delta > eps && iterations < MAXIT
        iterations = iterations + 1;
        %% Update G
        num = X*H';
        denom = (G*H)*H'+eps;
        G = G.*(num./denom);
        %% Update H
        num = G'*X;
        denom = G'*(G*H)+eps;
        H = H.*(num./denom);
        %% Standardize       
        H = H .* repmat(sqrt(sum(G.^2,1)+eps)', [1 C]); 
        G = G ./ repmat(sqrt(sum(G.^2,1)+eps), [N 1]);  
        GH = G*H; GH = GH;
        err = [err sum(sum((X-GH).^2))];
        delta = err(iterations-1)-err(iterations);
    end
    Z = err(iterations); zz(ijkl) = Z;
    if Z < zbest
        zbest = Z; Hbest = H; Gbest = G;
    end
end
Z = zbest; H = Hbest; G = Gbest;
% Check
GH = G*H; GH = GH;
Z = sum(sum((X-GH).^2));
Xbar = (1./(N.*C)).*(sum(sum(X)));
Xmat = repmat(Xbar,N,C);
Xss  = ((X-Xmat).^2);
TSS = sum(sum(Xss));
PRE = 1 - (Z./TSS);
% ###################################################################
% Use K-means to Cluster G and H
% ###################################################################
X = [];
X = G;
[n,dd] = size(X);
c = K;
sbest = 9999999999999999; itermax = 500; irep = 0; toggle = 1;
while irep < itermax
    irep = irep + 1;
    q = randperm(n);
    q = q(1:c);
    centroids = X(q,:);
    pp = zeros(1,n);
    for i = 1:n
        distbest = 999999999999;
        for k = 1:c
            diff = X(i,:)-centroids(k,:);
            dist = sum(diff.^2);
            if dist < distbest
                distbest = dist;
                ksel = k;
            end
        end
        pp(i) = ksel;
    end
    [p, totwcss] = hkmeans(X, pp, toggle);
    if totwcss < sbest
        sbest = totwcss; pbest = p;
    end
end 
prow = pbest;
X = H';
[m,dd] = size(X);
c = L;
sbest = 9999999999999999; itermax = 500; t = 0; irep = 0; toggle = 1;
while irep < itermax
    irep = irep + 1;
    q = randperm(m);
    q = q(1:c);
    centroids = X(q,:);
    pp = zeros(1,m);
    for i = 1:m
        distbest = 999999999999;
        for k = 1:c
            diff = X(i,:)-centroids(k,:);
            dist = sum(diff.^2);
            if dist < distbest
                distbest = dist;
                ksel = k;
            end
        end
        pp(i) = ksel;
    end
    [p, totwcss] = hkmeans(X, pp, toggle);
    if totwcss < sbest
        sbest = totwcss; pbest = p;
    end
end 
pcol = pbest;

toc
