function [pr, pc, Z, vaf, attract] = tmklmp(X, K, L)
% tmklmp - two-mode KL-means partitioning
% 
% Input:
%   X  - the two-mode input matrix 
%   K  - number of clusters for rows
%   L  - number of clusters for columns
%
% Results:
%   Z    - the sum-of-squares loss function value
%   vaf  - varianc-accounted for
%   pr - partition of the row objects
%   pc - partition of the column objects
%   attract - percentage of restarts tielding best VAF
%

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

if ndims(X) ~= 2
    error('X must be 2 dimensionsal');
end
[n,m] = size(X);

% Outer Loop for Multiple Restarts
tic;
state = 2;                  % fix state for testing
rand('state', state);
ZZ = 999999999999;

restarts = 500;
for ijkl = 1:restarts

% Make the initial row cluster assignments
q = randperm(n);
q = q(1:K);
rcent = X(q,:);
r = zeros(n,1);
for k = 1:K
    r(q(k)) = k;
end
for i = 1:n
    if r(i) == 0
        mindist = 999999999999;
        for k = 1:K
            dist = sum((X(i,:)-rcent(k,:)).^2);
            if dist < mindist
                mindist = dist;
                ksel = k;
            end
        end
        r(i) = ksel;
    end
end

% Make the initial column cluster assignments
q = randperm(m);
q = q(1:L);
ccent = X(:,q);
c = zeros(m,1);
for k = 1:L
    c(q(k)) = k;
end
for i = 1:m
    if c(i) == 0
        mindist = 99999999999;
        for k = 1:L
            dist = sum((X(:,i)-ccent(:,k)).^2);
            if dist < mindist
                mindist = dist;
                ksel = k;
            end
        end
        c(i) = ksel;
    end
end

% Compute the cluster centroids
v = zeros(K,L);
for k = 1:K
    for l = 1:L
        v(k,l) = mean(mean(X(r==k,c==l)));
    end
end
% Compute the Objective Function
Zbest = 0;
for i = 1:n
    for j = 1:m
        Zbest = Zbest + (X(i,j)-v(r(i),c(j))).^2;
    end
end
trig = 0;

while trig == 0
    trig = 1;
    % Reassign row objects
    mindist = 999999999999.*ones(n,1);
    for i = 1:n
        for k = 1:K
            dist = 0;
            for j = 1:m
              dist = dist + (X(i,j)-v(k,c(j))).^2;
            end
            if dist < mindist(i)
                mindist(i) = dist;
                ksel = k;
            end
        end
        r(i) = ksel;
     end

    % Ensure that there are no empty row clusters.
    checksum = 0;
    while checksum < K
        checksum = 0;
        for k = 1:K
            nk = sum(r==k);
            if nk >= 1
                checksum = checksum + 1;
                continue
            else
                [dsel,isel] = max(mindist);
                r(isel) = k;
                mindist(isel) = 0;
                break
            end
        end
    end

    % Compute the cluster centroids
    v = zeros(K,L);
    for k = 1:K
        for l = 1:L
            v(k,l) = mean(mean(X(r==k,c==l)));
        end
    end

    % Reassign column objects
    mindist = 999999999999.*ones(m,1);
    for j = 1:m
        for k = 1:L
            dist = 0;
            for i = 1:n
                dist = dist + (X(i,j)-v(r(i),k)).^2;
            end
            if dist < mindist(j)
                mindist(j) = dist;
                ksel = k;
            end
        end
        c(j) = ksel;
    end

    % Ensure that there are no empty column clusters.
    checksum = 0;
    while checksum < L
        checksum = 0;
        for k = 1:L
            nk = sum(c==k);
            if nk >= 1
                checksum = checksum + 1;
                continue
            else
                [dsel,jsel] = max(mindist);
                c(jsel) = k;
                mindist(jsel) = 0;
                break
            end
        end
    end
    
    % Compute the cluster centroids
    v = zeros(K,L);
    for k = 1:K
        for l = 1:L
            v(k,l) = mean(mean(X(r==k,c==l)));
        end
    end

% Compute the Objective Function
    Z = 0;
    for i = 1:n
        for j = 1:m
            Z = Z + (X(i,j)-v(r(i),c(j))).^2;
        end
    end

    if Z < Zbest - .0000001;
        Zbest = Z;
        trig = 0;
        pr = r;
        pc = c;
    end
end

if Zbest < ZZ
    ZZ = Zbest;
    prb = pr;
    pcb = pc;
    attract = 1;
elseif Zbest < ZZ+.0001
    attract = attract + 1;
end

end
    

% Compute final statistics
Z = ZZ;
pr = prb;
pc = pcb;
LL = zeros(n.*m,1);
idx = 0;
for i = 1:n
    for j = 1:m
        idx = idx + 1;
        LL(idx) = X(i,j);
    end
end
grandvar = var(LL).*(n.*m-1);
vaf = (grandvar-Z)./grandvar;
attract = attract./restarts;
toc

