-
Notifications
You must be signed in to change notification settings - Fork 0
/
FKT.py
177 lines (138 loc) · 4.96 KB
/
FKT.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
"""
Created on Wed May 15 11:54:57 2019
@author: daryl
#######################################################################
# BASED On An Original Sage Version by:
# 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: The number of perfect matchings of G
########################################################################
"""
import networkx as nx #Requires at least netwrokx 2.3+
import matplotlib.pyplot as plt
import random
import math
import numpy as np
import time
#Helper Functions
def doNothing():
return 0;
def find_faces(embd):
#Returns a list of faces of the planar embedding by
#the edges that bound the face
ed_list = list(embd.edges())
faces=[]
for ed in embd.edges():
if ed in ed_list:
faces.append(embd.traverse_face(ed[0],ed[1]))
for i in range(len(faces[-1])):
ed_list.remove((faces[-1][i],faces[-1][(i+1)%len(faces[-1])]))
face_edges=[]
for face in faces:
face_edges.append([])
for i in range(len(face)):
face_edges[-1].append((face[i],face[(i+1)%len(face)]))
return face_edges
def toSkewSymmetricMatrix(A):
#Skew--symmetrize a matrix
A[(A==1).T] = -1
return A
def numberOfClockwiseEdges(face, edgesT1):
#Iterate over edges of a face to determine
#the number of positive orientations
clockwise = 0
for edge in face:
try:
edgesT1.index(edge)
except ValueError:
doNothing()
else:
clockwise += 1
return clockwise
def isClockwise(e,face):
#Checks orientation of an edge on a face
try:
face.index(e);
except ValueError:
return False
else:
return True
#Main Function
def FKT(A):
n = len(A)
B_graph = A[:]
G = nx.Graph(A)
tf, embd = nx.check_planarity(G)
if embd is None:
return 0
faces = find_faces(embd)
T1 = nx.minimum_spanning_tree(G)
T1 = nx.Graph(T1)
mask = np.random.randint(2, size=(n, n))
mask = ((mask + mask.T) == 1)
B_digraph = A * mask
G = nx.DiGraph(B_digraph)
edgesT1 = T1.edges();
adj_T1 = (nx.adjacency_matrix(T1)).todense();
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 = nx.DiGraph(adj_T1)
edgesT1 = list(T1.edges())
if embd is not None:
faces.sort(key=len)
faces.reverse()
faces.pop(0)
if embd is not None:
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((edge[1], edge[0]));
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 ((numberOfClockwiseEdges(face, edgesT1)) % 2 == 1):
# 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 = nx.DiGraph(adj_T1);
edgesT1 = list(T1.edges());
# remove the face that was found
faceFound = faces.pop(index);
break;
try:
return math.sqrt(np.linalg.det(toSkewSymmetricMatrix(adj_T1)));
except ValueError:
pass