% generate forest compatibility matrix
function M = ForestMatrix(n)
    P = SetPartition(n);
    bell = size(P,1);
    M = zeros(bell, bell);
    for k = 1:bell
        for l = 1:bell
            M(k,l) = ForestJoin(P, k, l);
        end
    end
end

% Checks if two partitions are compatible
function output = ForestJoin(P, i,j)
    a = size(P{i},2); % size of partition i
    b = size(P{j},2); % size of partition j
    M = zeros(a,b); % intersect. size between sets in the partitions
    % computing entries of M
    for k = 1:a
        for l = 1:b
            M(k,l) = length(intersect(P{i}{k}, P{j}{l}));
            if M(k,l) > 1 % two sets k, l induce a cycle
                output = 0;
                return
            end
        end
    end
    % check for cycles that use more than 2 sets
    G = graph([zeros(a,a) M; M' zeros(b,b)]);
    output = 1 - hascycles(G);
end