An implementation of the FKT-algorithm

This post is dedicated to give you an implementation (in SAGE) of the FKT-algorithm. I coded this while spending some time on the ideas i explained in this blog post. I searched several hours to find an existing implementation without success. Since it took me a while to code this on my own, i think it is worth to present it, therewith other can use it in their own work. It works well, but if you find an error, please leave a note.

The FKT-algorithm can be used to count the number of perfect matchings in a planar graph. And a graph is planar if and only if it does not have K$_5$ or K$_{3,3}$ as a minor.

A paper that will be presented at the CCC conference this year is called:

Counting the Number of Perfect Matchings in $K_5$-free Graphs [in polynomial time] Simon Straub (Ulm University), Thomas Thierauf (Aalen University), Fabian Wagner (Ulm University)

I did not manage to read it yet, but i am curious to see, how they circumvent the K$_{3,3}$ minor problem.

########################################################################
# AUTHOR: Dr. Christian Schridde
# E-MAIL: christianschridde [at] googlemail [dot] com
#
# DESCRIPTION: Implementation of the FKT-algorithm
#
# INPUT:  Adjacency matrix A of a undirected loop-free planar graph G
# OUTPUT: Skew matrix M, such that PerfMatch(G) = Sqrt(Determinat(M))
########################################################################

def FKT(A):
 # make some local copies
 B_graph = copy(A);
 B_digraph = copy(A);

 # FIRST: REMOVE SELF-LOOPS
 for i in range(0,B_graph.nrows()):
  B_graph[i,i] = 0;
  B_digraph[i,i] = 0;

 # SECOND: TRANSFORM B_graph TO AN UNDIRECTED GRAPH
 G = DiGraph(B_graph);
 G = G.to_undirected();
 B_graph = G.adjacency_matrix();

 # THIRD: COMPUTE A PLANAR EMBEDDING
 if (G.is_planar(set_embedding=true, set_pos=true) == false):
  print "G is not planar";
  return null;
 embd = G.get_embedding();
 # each face is listed in clockwise order, except the first one which is the outer face
 faces = G.trace_faces(embd);
     
 # FOURTH: GET THE SPANNING TREE
 T1 = G.min_spanning_tree(algorithm='Prim_fringe');
 T1 = Graph(T1);

 # FIFTH: ASSIGN AN ARBITRARY DIRECTION TO THE EDGES
 DG = DiGraph(B_digraph);
 DG = DG.to_undirected();
 B_digraph = DG.adjacency_matrix();
 for i in range(0,B_digraph.nrows()):
  for j in range(0,B_digraph.nrows()):
   if (B_digraph[i,j] == 1):
    if (B_digraph[j,i] == 1):
     r = randint(0,1)
     if (r==0):
      B_digraph[j,i] = 0;
     else:
      B_digraph[i,j] = 0;
 G = DiGraph(B_digraph);

 # SIXTH: FIND FACE THAT HAS ONLY ONE EDGE NOT IN T1
 edgesT1 = T1.edges(labels=false);
 adj_T1 = T1.adjacency_matrix();
     
 # give the edges in T1 the orientation from B_digraph
 for edge in edgesT1:
  if (B_digraph[edge[0],edge[1]] == 0):
   adj_T1[edge[0],edge[1]] = 0;
  else:
   adj_T1[edge[1],edge[0]] = 0;
             
 T1 = DiGraph(adj_T1);
 edgesT1 = T1.edges(labels=false);      
     
 # remove the first face, which is the outer face
 # QUESTION: Is it always correct to remove this face?
 faces.pop(0);
     
 while (len(faces) > 0):
  index = -1;
  for face in faces:
   countMissingEdges = 0; 
   missingEdge = 0;
   index += 1;
   for edge in face:
    try:
     idx1 = edgesT1.index(edge);
    except ValueError:
     try:
      idx2 = edgesT1.index(reverseEdge(edge));
     except ValueError:
      countMissingEdges += 1;
      missingEdge = edge;
     else:
      doNothing();
    else:
     doNothing();

   if (countMissingEdges == 1):
    # in this face, only one edge is missing.        
    # Place the missing edge such that the total number
    # of clockwise edges of this face is odd
    # add this edge to the spanning tree
    if (is_odd(numberOfClockwiseEdges(face,edgesT1))):
     # insert counterclockwise in adj_T1;
     if (isClockwise(missingEdge,face) == false):
      adj_T1[missingEdge[0],missingEdge[1]] = 1;
     else:
      adj_T1[missingEdge[1],missingEdge[0]] = 1;
    else:
     # insert clockwise in adj_T1
     if (isClockwise(missingEdge,face) == true):
      adj_T1[missingEdge[0],missingEdge[1]] = 1;
     else:
      adj_T1[missingEdge[1],missingEdge[0]] = 1;
                             
    # rebuild the graph
    T1 = DiGraph(adj_T1);
    edgesT1 = T1.edges(labels=false);
                             
    # remove the face that was found
    faceFound = faces.pop(index);
    break;

 return toSkewSymmetricMatrix(adj_T1);



Also, you need the following helper functions:

############################
# Returns true if the given edge is
# clockwise oriented regarding the given face
def isClockwise(e,face):
  try:
   face.index(e);
  except ValueError:
   return false;
  else:
   return true;


############################
# This is a placeHolder function
def doNothing():
 return 0;


############################ 
# Reverses a given edge
def reverseEdge(edge):
 return (edge[1],edge[0]);


############################
# Inputs are the face and the orientedEdges from the spanning-tree T1
# Note, that all edges in face are clockwise (since it is the result
# from the obtained embedding).
# It returns how many of the edges from T1 that are part of the
# given face are actually clockwise.
def numberOfClockwiseEdges(face, edgesT1):
 clockwise = 0;
 for edge in face:
   try:
     edgesT1.index(edge);
   except ValueError:
     doNothing();
   else:
     clockwise += 1;
 return clockwise;


############################
# Transforms a given matrix A to a 
# skewSymmetric Matrix
def toSkewSymmetricMatrix(A):
  for i in range(0,A.nrows()):
    for j in range(0,A.nrows()):
      if (A[i,j] == 1):
        A[j,i] = -1;
  return A;