How to modify Johnson's elementary cycles algorithm to cap maximum cycle length?

Viewed 936

I would like to modify the networkx implementation of Johnson's algorithm for finding all elementary cycles in a graph (also copied below) so that is does not search for cycles larger than some maximum length.

def simple_cycles(G):
    def _unblock(thisnode,blocked,B):
        stack=set([thisnode])
        while stack:
            node=stack.pop()
            if node in blocked:
                blocked.remove(node)
                stack.update(B[node])
                B[node].clear()

    # Johnson's algorithm requires some ordering of the nodes.
    # We assign the arbitrary ordering given by the strongly connected comps
    # There is no need to track the ordering as each node removed as processed.
    subG = type(G)(G.edges_iter()) # save the actual graph so we can mutate it here
                              # We only take the edges because we do not want to
                              # copy edge and node attributes here.
    sccs = list(nx.strongly_connected_components(subG))
    while sccs:
        scc=sccs.pop()
        # order of scc determines ordering of nodes
        startnode = scc.pop()
        # Processing node runs "circuit" routine from recursive version
        path=[startnode]
        blocked = set() # vertex: blocked from search?
        closed = set() # nodes involved in a cycle
        blocked.add(startnode)
        B=defaultdict(set) # graph portions that yield no elementary circuit
        stack=[ (startnode,list(subG[startnode])) ]  # subG gives component nbrs
        while stack:
            thisnode,nbrs = stack[-1]
            if nbrs:
                nextnode = nbrs.pop()
#                    print thisnode,nbrs,":",nextnode,blocked,B,path,stack,startnode
#                    f=raw_input("pause")

                if nextnode == startnode:
                    yield path[:]
                    closed.update(path)
#                        print "Found a cycle",path,closed
                elif nextnode not in blocked:
                    path.append(nextnode)
                    stack.append( (nextnode,list(subG[nextnode])) )
                    closed.discard(nextnode)
                    blocked.add(nextnode)
                    continue
            # done with nextnode... look for more neighbors
            if not nbrs:  # no more nbrs
                if thisnode in closed:
                    _unblock(thisnode,blocked,B)
                else:
                    for nbr in subG[thisnode]:
                        if thisnode not in B[nbr]:
                            B[nbr].add(thisnode)
                stack.pop()
                assert path[-1]==thisnode
                path.pop()

        # done processing this node
        subG.remove_node(startnode)
        H=subG.subgraph(scc)  # make smaller to avoid work in SCC routine
        sccs.extend(list(nx.strongly_connected_components(H)))

Of course, I'd also accept a suggestion that differs from the implementation above but runs in similar time. Also, my project uses networkx, so feel free to use any other function from that library, such as shortest_path.

(Note: not homework!)

Edit

Dorijan Cirkveni suggested (if I understood correctly):

        if len(blocked) >= limit + 1:
            continue

        elif nextnode == startnode:
            yield path[:]

However, that doesn't work. Here's a counterexample:

G = nx.DiGraph()
G.add_edge(1, 2)
G.add_edge(2, 3)
G.add_edge(3, 1)
G.add_edge(3, 2)
G.add_edge(3, 4)

my_cycles = list(simple_cycles(G, limit = 3)) # Modification
nx_cycles = list(nx.simple_cycles(G)) # Original networkx code
print("MY:", my_cycles)
print("NX:", nx_cycles)

Will output

MY: [[2, 3]]
NX: [[1, 2, 3], [2, 3]]

Also, if we substitute blocked by stack or path, the result will be correct for this example, but will give the wrong answer for other graphs.

2 Answers

This is a highly modified version of this code, but at least it is working.

def simple_cycles(G, limit):
    subG = type(G)(G.edges())
    sccs = list(nx.strongly_connected_components(subG))
    while sccs:
        scc = sccs.pop()
        startnode = scc.pop()
        path = [startnode]
        blocked = set()
        blocked.add(startnode)
        stack = [(startnode, list(subG[startnode]))]

        while stack:
            thisnode, nbrs = stack[-1]

            if nbrs and len(path) < limit:
                nextnode = nbrs.pop()
                if nextnode == startnode:
                    yield path[:]
                elif nextnode not in blocked:
                    path.append(nextnode)
                    stack.append((nextnode, list(subG[nextnode])))
                    blocked.add(nextnode)
                    continue
            if not nbrs or len(path) >= limit:
                blocked.remove(thisnode)
                stack.pop()
                path.pop()
        subG.remove_node(startnode)
        H = subG.subgraph(scc)
        sccs.extend(list(nx.strongly_connected_components(H)))

You only need to change two things:

  1. The definition line (obviously)

    def simple_cycles(G,limit):

  2. Add an overriding condition somewhere in the next node processor (example below:)

            ...
            if blocked.size>=limit+1:
                pass
            elif if nextnode == startnode:
                yield path[:] ...
    

    Bonus: Using == instead of >= will result in the function running as there is no limit when a negative value is used, as opposed to not returning any nodes.

Related