Finding subgraphs with constraint programming and Google OR-tools

Here's a graph:

Suppose you want to find all subgraphs that look like this:

(That is, two adjacent edges.)

One way to implement such a search is to use Constraint Programming, and Google has a library that includes a constraint programming solver: or-tools. OR-tools has a Python SWIG wrapper which makes it (sort of) easy to use in Python. Unfortunately most of the online documentation is for the C++ interface, but all the function names are the same.

Here's an example of how to use or-tools for the graph problem above, step by step. I made heavy use of the quickstart guide.

Import the python library:

from ortools.sat.python import cp_model

Represent the graph by the set of (bidirectional) edges it contains.

edges = [ (1,2), (1,3), (1,4), (3,4), (5,6) ]
numVertices = 6

Create a model with a variable for each vertex in the subgraph we're trying to find:

model = cp_model.CpModel()

# Find two adjacent edges
a = model.NewIntVar( 1, numVertices, "a" )
b = model.NewIntVar( 1, numVertices, "b" )
c = model.NewIntVar( 1, numVertices, "c" )

This syntax means that a can take on the integer values 1 through 6. We'll interpret that as a vertex in our graph.

Constrain the model so that the only allowed assignments are of vertexes connected by an edge. We need to allow the edge to be used in either direction, though. (We could just as easily solve the problem for directed graphs by not including the reversed edges.)

reversedEdges = [ (x,y) for (y,x) in edges ]
model.AddAllowedAssignments( [a,b], edges + reversedEdges )
model.AddAllowedAssignments( [b,c], edges + reversedEdges )
# If we use inequality then we get both a-b-c and c-b-a as solutions.
model.Add( a < c )

As the comment says, we want to prohibit x-y-x as a solution but if we just add a constraint that a and c are are not the same, then we'll get the same subgraph twice.

The solver will give us solutions one at a time (if we want) so we need a callback to save or print them:

class SolutionPrinter( cp_model.CpSolverSolutionCallback ):
    def __init__( self ):
        cp_model.CpSolverSolutionCallback.__init__(self)
        self.count = 0

    def OnSolutionCallback( self ):
        self.count += 1
        print( "{} -- {} -- {}".format( self.Value( a ),
                                        self.Value( b ),
                                        self.Value( c ) ) )

Under the hood SWIG is doing something really cool: this is a Python class that's "derived" from a C++ class! The C++ code gets a C++ object, and that C++ object's methods call the corresponding Python code.

The quickstart guide has a documentation error as I write this; it calls the method NewSolution but pydoc says it is OnSolutionCallback. You'll get an obscure error message from SWIG if you get it wrong.

Finally, instantiate our printer and a solver, and tell it we want all the solutions:

sp = SolutionPrinter()
solver = cp_model.CpSolver()
status = solver.SearchForAllSolutions( model, sp )
print( "Status:", status )
print( "Total number of subgraphs:", sp.count )

Output:

2 -- 1 -- 3
1 -- 3 -- 4
2 -- 1 -- 4
1 -- 4 -- 3
3 -- 1 -- 4
Status: 2
Total number of subgraphs: 5

Internally, the constraint solver is translating the problem into a boolean satisfiability (SAT) instance. We can get the number of variables used

>>> print( "Number of SAT variables:", solver.NumBooleans() )
47

but I can't figure out how to dump the model itself. There's a debug variable that could be set in the C++ code: https://github.com/google/or-tools/blob/master/ortools/sat/cp_model_solver.cc, but I don't know if I can access that variable from Python. I'm unclear how to set any of the non-debug settings for the SAT solver, of which there are a ton.