constraint programming mesh network

Question

I have a mesh network as shown in figure. Now, I am allocating values to all edges in this sat network. I want to propose in my program that, there are no closed loops in my allocation. For example the constraint for top-left most square can be written as -

E0 = 0 or E3 = 0 or E4 = 0 or E7 = 0, so either of the link has to be inactive in order not to form a loop. However, in this kind of network, there are many possible loops.

For example loop formed by edges - E0, E3, E7, E11, E15, E12, E5, E1.

Now my problem is that I have to describe each possible combination of loop which can occur in this network. I tried to write constraints in one possible formula, however I was not able to succeed.

Can anyone throw any pointers if there is a possible way to encode this situation? Just for information, I am using Z3 Sat Solver.

Answer 1

The following encoding can be used with any graph with N nodes and M edges. It is using (N+1)*M variables and 2*M*M 3-SAT clauses. This ipython notebook demonstrates the encoding by comparing the SAT solver results (UNSAT when there is a loop, SAT otherwise) with the results of a straight-forward loop finding algorithm.

Disclaimer: This encoding is my ad-hoc solution to the problem. I'm pretty sure that it is correct but I don't know how it compares performance-wise to other encodings for this problem. As my solution works with any graph it is to be expected that a better solution exists that uses some of the properties of the class of graphs the OP is interested in.

Variables:

I have one variable for each edge. The edge is "active" or "used" if its corresponding variable is set. In my reference implementation the edges have indices 0..(M-1) and this variables have indices 1..M:

def edge_state_var(edge_idx):
    assert 0 <= edge_idx < M
    return 1 + edge_idx

Then I have an M bits wide state variable for each edge, or a total of N*M state bits (nodes and bits are also using zero-based indexing):

def node_state_var(node_idx, bit_idx):
    assert 0 <= node_idx < N
    assert 0 <=  bit_idx < M
    return 1 + M + node_idx*M + bit_idx

Clauses:

When an edge is active, it links the state variables of the two nodes it connects together. The state bits with the same index as the node must be different on both sides and the other state bits must be equal to their corresponding partner on the other node. In python code:

# which edge connects which nodes
connectivity = [
    ( 0,  1), # edge E0
    ( 1,  2), # edge E1
    ( 2,  3), # edge E2
    ( 0,  4), # edge E3
    ...
]

cnf = list()

for i in range(M):
    eb = edge_state_var(i)
    p, q = connectivity[i]
    for k in range(M):
        pb = node_state_var(p, k)
        qb = node_state_var(q, k)
        if k == i:
            # eb -> (pb != qb)
            cnf.append([-eb, -pb, -qb])
            cnf.append([-eb, +pb, +qb])
        else:
            # eb -> (pb == qb)
            cnf.append([-eb, -pb, +qb])
            cnf.append([-eb, +pb, -qb])

So basically each edge tries to segment the graph it is part of into a half that is on one side of the edge and has all the state bits corresponding to the edge set to 1 and a half that is on the other side of the edge and has the state bits corresponding to the edge set to 0. This is not possible for a loop where all nodes in the loop can be reached from both sides of each edge in the loop.