function [P, Q, Z, newfp, attract] = tmbp(A, K, L)
% ########################################################################
% RELOCATION ALGORITHM FOR 2-MODE CLUSTERING OF BINARY NETWORK MATRIX
%   - ALGORITHM IS GREEDY - FIND BEST MOVE BEFORE ACCEPTING
%   - ALGORITHM DESIGNED FOR INDUCTIVE BLOCKMODELING
%   - ALGORITHM WILL COLLECT EQUALLY-WELL FITTING PARTITIONS
%   - ALGORITHM MINIMIZES THE SUM OF 1'S IN NULL BLOCKS + 0'S IN COMPLETE
% ALGORITHM INPUTS:
%   A = AN N-BY-M BINARY NETWORK MATRIX
%   K = THE NUMBER OF CLUSTERS FOR ROW OBJECTS
%   L = THE NUMBER OF CLUSTERS FOR COLUMN OBJECTS
% ALGORITHM OUTPUTS
%   Z = THE BEST-FOUND VALUE OF THE CRITERION FUNCTION
%   NEWFP = THE NUMBER OF EQUALLY-WELL FITTING PARTITIONS FOUND
%   P = A MATRIX CONTAINING EQUALLY WELL-FITTING PARTITIONS OF ROW OBJECTS
%   Q = A MATRIX CONTAINING EQUALLY WELL-FITTING PARTITIONS OF COL OBJECTS
%   attract = the percentage of repetitions yielding minimum inconsistencies
% ########################################################################
%---------------------------------------
% Make sure program is usable
%---------------------------------------
if ndims(A) ~= 2
    error('X must be 2 dimensionsal');
end
[N, M] = size(A); nk = ceil(N/K); nl = ceil(M/L); tic;

state = 2;                  % fix state for testing
rand('state', state);
Zbest = 9999999;
repetitions = 5000;
for ijkl = 1:repetitions
    p = (1:K); p = repmat(p,1,nk); v = randperm(nk.*K); p = p(v); p = p(1:N);
    q = (1:L); q = repmat(q,1,nl); v = randperm(nl.*L); q = q(v); q = q(1:M);
    Z=0; nr=zeros(1,K); nc = zeros(1,L); z1 = zeros(K,L); z0 = zeros(K,L);
    if length(unique(p)) < K || length(unique(q)) < L
        continue
    end
    for k = 1:K
        nr(k) = sum(p==k); 
        for l = 1:L
            nc(l) = sum(q==l); 
            B = A(p==k,q==l); z1(k,l) = sum(sum(B)); [u,v] = size(B);
            z0(k,l) = 0;
            for i = 1:N
                for j = 1:M
                    if p(i) == k && q(j) == l && A(i,j) == 0
                        z0(k,l) = z0(k,l) + 1;
                    end
                end
            end
            Z = Z + min(z0(k,l),z1(k,l));
        end
    end
    dbest = -99999;
    while dbest < 0
        dbest = 999999;
        for i = 1:N
            k1 = p(i);
            if nr(k1) == 1
                continue
            end
            for k2 = 1:K
                if k2 == k1
                    continue
                end
                delta = 0;
                for l = 1:L
                    delta = delta - min(z1(k1,l),z0(k1,l))-min(z1(k2,l),z0(k2,l));
                    z1trial(k1,l) = z1(k1,l); z0trial(k1,l) = z0(k1,l);
                    z1trial(k2,l) = z1(k2,l); z0trial(k2,l) = z0(k2,l);
                end
                for j = 1:M
%                    if i == j 
%                        continue
%                    end
                    l = q(j);
                    if A(i,j) == 1
                        z1trial(k1,l) = z1trial(k1,l)-1; z1trial(k2,l) = z1trial(k2,l)+1;
                    else
                        z0trial(k1,l) = z0trial(k1,l)-1; z0trial(k2,l) = z0trial(k2,l)+1;
                    end
                end
                for l = 1:L
                    delta = delta + min(z1trial(k1,l),z0trial(k1,l))+min(z1trial(k2,l),z0trial(k2,l));
                end
                if delta < dbest
                    dbest = delta; imove = 1; isel1 = i; ksel1 = k1; ksel2 = k2;
                end
            end
        end
        for j = 1:M
            l1 = q(j);
            if nc(l1) == 1
                continue
            end
            for l2 = 1:L
                if l2 == l1
                    continue
                end
                delta = 0;
                for k = 1:K
                    delta = delta - min(z1(k,l1),z0(k,l1))-min(z1(k,l2),z0(k,l2));
                    z1trial(k,l1) = z1(k,l1); z0trial(k,l1) = z0(k,l1);
                    z1trial(k,l2) = z1(k,l2); z0trial(k,l2) = z0(k,l2);
                end
                for i = 1:N
 %                   if i == j 
 %                       continue
 %                   end
                    k = p(i);
                    if A(i,j) == 1
                        z1trial(k,l1) = z1trial(k,l1)-1; z1trial(k,l2) = z1trial(k,l2)+1;
                    else
                        z0trial(k,l1) = z0trial(k,l1)-1; z0trial(k,l2) = z0trial(k,l2)+1;
                    end
                end
                for k = 1:K
                    delta = delta + min(z1trial(k,l1),z0trial(k,l1))+min(z1trial(k,l2),z0trial(k,l2));
                end
                if delta < dbest
                    dbest = delta; imove = 2; jsel1 = j; lsel1 = l1; lsel2 = l2;                   
                end
            end
        end
        if dbest < 0
            Z = Z + dbest;
            if imove == 1
                nr(ksel1) = nr(ksel1) - 1; nr(ksel2) = nr(ksel2) + 1;
                p(isel1) = ksel2;
                for j = 1:M
                    l = q(j);
 %                   if j == isel1
 %                       continue
 %                   end
                    if A(isel1,j) == 1
                        z1(ksel1,l) = z1(ksel1,l)-1; z1(ksel2,l) = z1(ksel2,l)+1;
                    else
                        z0(ksel1,l) = z0(ksel1,l)-1; z0(ksel2,l) = z0(ksel2,l)+1;
                    end
                end
            else
                nc(lsel1) = nc(lsel1) - 1; nc(lsel2) = nc(lsel2) + 1;
                q(jsel1) = lsel2;
                for i = 1:N
                    k = p(i);
 %                   if i == jsel1
 %                       continue
 %                   end
                    if A(i,jsel1) == 1
                        z1(k,lsel1) = z1(k,lsel1)-1; z1(k,lsel2) = z1(k,lsel2)+1;
                    else
                        z0(k,lsel1) = z0(k,lsel1)-1; z0(k,lsel2) = z0(k,lsel2)+1;
                    end
                end
            end
        end
    end
    
    if Z < Zbest
        Zbest = Z; P = p; Q = q; newfp = 1; attract = 1;
    elseif Z == Zbest
        attract = attract + 1;
        for iii = 1:newfp
            inew = 0;
            for i = 1:N-1
                for j = i+1:N
                    if p(i) == p(j) && P(iii,i) ~= P(iii,j)
                        inew = 1;
                        continue
                    end
                end
                if inew == 1
                    break
                end
            end
            if inew == 1
                break
            end
            for i = 1:M-1
                for j = i+1:M
                    if q(i) == q(j) && Q(iii,i) ~= Q(iii,j)
                        inew = 1;
                        continue
                    end
                end
                if inew == 1
                    break
                end
            end
            if inew == 0
                break
            end
        end
        if inew == 0
            continue  %break
        end
        if inew == 1
            newfp = newfp + 1;
            P = [P; p]; Q = [Q;q];
        end
    end
end
Z = Zbest;
attract = attract./repetitions;
toc
                        
        
            
            
