%{
*********************************************************** 
(c) Ian Davidson and Davidson Lab. Please cite 

The Cluster Description Problem-Complexity Results, Formulations and Approximations, I Davidson, A Gourru and SS. Ravi
Advances in Neural Information Processing Systems, 2018, 6191-6201
***********************************************************
%}


function x = basic_apart_together_func(T, pi, w, R, P, tagNames)

k = max(pi(:));

[n, t] = size(T);

idx_deleted = [];

%% Together Constraints

if ~isempty(R)
    if ~isequal(R, R')
        error('R must be symmetric');
    end
    [r, c] = find(R);
    G = digraph(r, c);
    H = transclosure(G);
    edges = table2array(H.Edges);
    edges(edges(:, 1) >= edges(:, 2), :) = [];
    if ~isempty(P)
        if ~isequal(P, P')
            error('P must be symmetric');
        end
    
        [crn, ~, ~] = chromnum(P, 0);
        if crn > k
            error('Apart constraints cannot be satisfied');
        end
        
        t_check = zeros(t, t);
        t_check(sub2ind(size(t_check), edges(:, 1), edges(:, 2))) = 1;
        t_check = t_check + t_check';
        if sum(sum(t_check .* P)) ~= 0
            error('Apart-Together constraints conflict');
        end
    end
    M = containers.Map('KeyType','double','ValueType','any');
    for i = 1 : size(edges, 1)
        v1 = edges(i, 1);
        v2 = edges(i, 2);
        if ~isKey(M, v1)
            if ~ismember(v1, idx_deleted)
                M(v1) = [v1, v2];
            end
        else
            M(v1) = [M(v1), v2];
        end
        if ~ismember(v2, idx_deleted) && ~isKey(M, v2)
            idx_deleted = [idx_deleted, v2];
        end
        if ~isempty(P)
            P(v1, P(v2, :) ~= 0) = 1;
        end
        T(T(:, v2) == 1, v1) = 1;
    end
    idx_remained = setxor((1 : t), idx_deleted);
    T(:, idx_deleted) = [];
    w(idx_deleted) = [];
    t = t - length(idx_deleted);
end

%% Apart Constraints
if ~isempty(P)
    
    P(:, idx_deleted) = [];
    P(idx_deleted, :) = [];
   
    [r, c] = find(P);
    rc = [r, c];
    rc(rc(:, 1) >= rc(:, 2), :) = [];
    [p, ~] = size(rc);
end

%% construct A \in R^{n \times kt}
A = zeros(n, k * t);
for i = 1 : n
    start_idx = t * (pi(i) - 1) + 1;
    A(i, start_idx : start_idx + t - 1) = T(i, :);
end
A_b = ones(n, 1);

%% construct B \in R^{t \times kt}
B = zeros(t, k * t);
for i = 1 : t
    b = zeros(1, t);
    b(1, i) = 1;
    B(i, :) = repmat(b, 1, k);
end

%% construct C \in R^{t \times kt}
C = []; C_b = [];
if ~isempty(P)
    C = zeros(k * p, k * t);
    for i = 1 : p
        offset_C = (i - 1) * p;
        c = zeros(1, t);
        c(1, [rc(i, 1), rc(i, 2)]) = 1;
        c_cell = repmat({c}, 1, k);
        c = blkdiag(c_cell{:});
        C(offset_C + 1 : offset_C + k, :) = c;
    end
    C_b = ones(k * p, 1);
end

A = [-A ; B ; C];

b = [-A_b ; w ; C_b];

f = ones(k * t, 1);

intcon = (1 : k * t);

x = intlinprog(f, intcon, A, b, [], [], zeros(k * t, 1), ones(k * t, 1));

x = uint8(x);
%% display tags
if ~isempty(R)
    if ~isempty(x)
        for i = 1 : k
            offset = (i - 1) * t + 1;
            xi = x(offset : offset + t - 1);
            d = ['tags for cluster ', num2str(i), ' :'];
            disp(d)
            for j = 1 : t
                if xi(j) == 1
                    if ~isKey(M, idx_remained(j))
                        disp(tagNames(idx_remained(j)))
                    else
                        idx_tags = M(idx_remained(j));
                        for l = 1 : length(idx_tags)
                            disp(tagNames(idx_tags(l)))
                        end
                    end
                end
            end
        end
    end
else
    for i = 1 : k
        offset = (i - 1) * t + 1;
        xi = x(offset : offset + t - 1);
        d = ['tags for cluster ', num2str(i), ' :'];
        disp(d)
        for j = 1 : t
            if xi(j) == 1
                disp(tagNames(j))
            end
        end
    end
end
