1
|
1 '''
|
|
2 This file contains functions used to calculate feedback vertex set for directed graph
|
|
3
|
|
4 The major function is the FVS_local_search(), which calculate maximum sub topological ordering of a graph
|
|
5 The algorithm is given by "Applying local search to feedback vertex set problem"
|
|
6 Author of the paper Philippe Galinier, Eunice Lemamou, Mohamed Wassim Bouzidi
|
|
7 The code mimic the pseudocode given in Page 805 in that paper
|
|
8 The code is written by Gang Yang, Penn State University
|
|
9 '''
|
|
10
|
|
11
|
|
12 import networkx as nx
|
|
13 import random
|
|
14 import math
|
|
15 import array
|
|
16
|
|
17
|
|
18
|
|
19 #the following two function calculate the position for the candidate given the existing topological ordering S
|
|
20 def get_position_minus(candidate_incoming_neighbour, S):
|
|
21 '''
|
|
22 get_position_minus return position as just after its numbered in-coming neighbours
|
|
23 As in the paper, the function return i_minus(v)
|
|
24 '''
|
|
25 position = 1
|
|
26 for x in range(len(S)-1,-1,-1): #we loop through index from the end as we try to find the largest index of numbered incoming neighbour
|
|
27 if S[x] in candidate_incoming_neighbour:
|
|
28 position = x+2 #2 comes from we need to put candidate after the incoming neighbour and also the position count from 1 instead of 0
|
|
29 return position
|
|
30 return position #put the candidate in the first position if there is no numbered incoming neighbour
|
|
31
|
|
32
|
|
33
|
|
34 def get_position_plus(candidate_outgoing_neighbour, S):
|
|
35 '''
|
|
36 get_position_plus return position as just before its numbered out-going neighbours
|
|
37 As in the paper, the function return i_plus(v)
|
|
38 '''
|
|
39 position = 1+len(S)
|
|
40 for x in range(len(S)): #we loop through index from the beginning as we try to find the smallest index of numbered outgoing neighbour
|
|
41 if S[x] in candidate_outgoing_neighbour:
|
|
42 position = x+1 #1 comes from the fact position count from 1 instead of 0
|
|
43 return position
|
|
44 return position #put the candidate in the first position if there is no numbered outgoing neighbour
|
|
45
|
|
46
|
|
47
|
|
48 def FVS_local_search(G_input, T_0, alpha, maxMvt, maxFail, randomseed=None):
|
|
49 '''
|
|
50 Returns an maximum sub topological ordering of a DiGraph G.
|
|
51 FVS is G_input \ the topological ordering
|
|
52 A topological ordering of a graph is an ordering of vertices such that
|
|
53 the starting point of every arc occurs earlier in the ordering than
|
|
54 the endpoint of the art.
|
|
55
|
|
56 This algorithm is a fast approximation based on simulating annealing(SA) of a noval local search strategy [1]_.
|
|
57 The code follows the pseudocode given in Page 805 in that paper.
|
|
58
|
|
59
|
|
60 Parameters
|
|
61 ----------
|
|
62 G : NetworkX Graph/DiGraph, result for MultiGraph is not tested
|
|
63 T_0 : the initial temperature in SA
|
|
64 alpha : the cooling multiplier for temperatue in the geometric cooling regime
|
|
65 maxMvt_factor : maxMvt_factor times network size is the number of iterations for the inner loop given a fixed temperatue
|
|
66 maxFail : FVS_local_search stops when maxFail number of outloops (temperatue) fail to improve the result
|
|
67 randomseed: random seed for generating random numbers
|
|
68
|
|
69 Returns
|
|
70 -------
|
|
71 An approximation of the maximum ordering of the given graph as a list.
|
|
72
|
|
73
|
|
74 Notes
|
|
75 -----
|
|
76 The code is written by Gang Yang, Department of Physics, Penn State University
|
|
77
|
|
78
|
|
79 References
|
|
80 ----------
|
|
81 ..[1] Galinier, P., Lemamou, E. & Bouzidi, M.W. J Heuristics (2013) 19: 797. doi:10.1007/s10732-013-9224-z
|
|
82
|
|
83 '''
|
|
84 #set a random seed
|
|
85 random.seed(randomseed)
|
|
86 G=G_input.copy()
|
|
87 N= len(G.nodes()) #number of nodes in the graph
|
|
88 edges=G.edges()
|
|
89
|
|
90
|
|
91 #Initialization
|
|
92 T = T_0 #set initial temperatue
|
|
93 nbFail = 0 #Outer loop counter to record failure times
|
|
94 S = [] #list to record the ascending ordering
|
|
95 S_optimal = [] #list to record the optimal ordering
|
|
96
|
|
97
|
|
98 #calculate parent and child node for each node
|
|
99 parent = [{} for i in range(N)]
|
|
100 child = [{} for i in range(N)]
|
|
101
|
|
102 for edge in edges:
|
|
103 child[int(edge[0])][int(edge[1])]=None
|
|
104 parent[int(edge[1])][int(edge[0])]=None
|
|
105
|
|
106 #all the nodes that are self_loops
|
|
107 self_loops = [edge[0] for edge in edges if edge[0]==edge[1]]
|
|
108 print(self_loops)
|
|
109 #all the nodes that is not in S
|
|
110 unnumbered=[x for x in range(N) if x not in self_loops]
|
|
111 N_unnumbered = len(unnumbered)
|
|
112 while nbFail< maxFail:
|
|
113 nbMvt = 0 #Inner loop counter to record movement times
|
|
114 failure = True #if cardinal of S increase after one inner loop, failure will be set to false
|
|
115 while nbMvt < maxMvt:
|
|
116 candidate_index = random.randint(0, N_unnumbered-1)
|
|
117 candidate = unnumbered[candidate_index] #random pick a node from all unnumbered node
|
|
118 position_type = random.randint(0,1) #random choose a position type
|
|
119 #calculate incoming/outgoing neighbour for the candidate node, store as keys in the dict
|
|
120 candidate_incoming_neighbour = parent[candidate]
|
|
121 candidate_outgoing_neighbour = child[candidate]
|
|
122 #see how to find the position on Page 803 of the paper
|
|
123 #position_type=1 means just after incoming neighbours
|
|
124 if position_type==1:
|
|
125 position = get_position_minus(candidate_incoming_neighbour,S)
|
|
126 #position_type=0 means just before outgoint neighbours
|
|
127 elif position_type==0:
|
|
128 position = get_position_plus(candidate_outgoing_neighbour,S)
|
|
129
|
|
130 #first, insert the candidate to the given position
|
|
131 S_trail = S[:] #copy the list
|
|
132 S_trail.insert(position-1,candidate)
|
|
133
|
|
134 #second remove all the conflict
|
|
135 #break the sequence into two parts: before and after the candidate node and
|
|
136 S_trail_head= S_trail[:position-1]
|
|
137 S_trail_tail= S_trail[position:]
|
|
138 #determine conflict node See page 801
|
|
139 if position_type==1:
|
|
140 CV_pos=[] #conflict before the newly inserted node in the topological ordering
|
|
141 for x in range(len(S_trail_head)):
|
|
142 nodetemp=S_trail_head[x]
|
|
143 #print nodetemp,candidate_outgoing_neighbour,nodetemp in candidate_outgoing_neighbour
|
|
144 if nodetemp in candidate_outgoing_neighbour:
|
|
145 CV_pos.append(nodetemp)
|
|
146 conflict=CV_pos #there won't be conflict after the inserted node as the node inserted after its incoming neighbour
|
|
147 elif position_type==0:
|
|
148 CV_neg=[] #conflict after the newly inserted node in the topological ordering
|
|
149 for x in range(len(S_trail_tail)):
|
|
150 nodetemp=S_trail_tail[x]
|
|
151 #print nodetemp,candidate_incoming_neighbour,nodetemp in candidate_incoming_neighbour
|
|
152 if nodetemp in candidate_incoming_neighbour:
|
|
153 CV_neg.append(nodetemp)
|
|
154 conflict=CV_neg #there won't be conflict before the inserted node as the node inserted before its outgoing neighbour
|
|
155 #finally remove the conflict node
|
|
156 N_conflict=len(conflict)
|
|
157 if N_conflict>0:
|
|
158 for i in range(N_conflict):
|
|
159 S_trail.remove(conflict[i])
|
|
160
|
|
161 #third, evaluate the move
|
|
162 #delta_move=-len(S_trail)+len(S)
|
|
163 delta_move=N_conflict-1
|
|
164 #accept all the move that is not harmful, otherwise use metrospolis algorithm
|
|
165 if delta_move<=0 or math.exp(-delta_move/float(T))>random.random():
|
|
166 S = S_trail[:]
|
|
167 #update unnumbered nodes
|
|
168 unnumbered.remove(candidate) #remove the node just inserted
|
|
169 if N_conflict>0: #add all the conflict node just removed
|
|
170 for i in range(N_conflict):
|
|
171 unnumbered.append(conflict[i])
|
|
172 N_unnumbered+=delta_move
|
|
173 nbMvt = nbMvt+1
|
|
174 #update S_optimal only when there is increase in cardinal of the sequence
|
|
175 if len(S)>len(S_optimal):
|
|
176 S_optimal = S[:]
|
|
177 failure = False
|
|
178 if N_unnumbered==0:
|
|
179 return S_optimal
|
|
180
|
|
181 #Increment the failure times if no progress in size of the sequence
|
|
182 if failure==True:
|
|
183 nbFail+=1
|
|
184 else: #otherwise reset the num of failure times
|
|
185 nbFail=0
|
|
186 #shrink the temperatue by factor alpha
|
|
187 T=T*alpha
|
|
188 #print T
|
|
189 #print nbFail
|
|
190 return S_optimal
|
|
191
|
|
192
|
|
193
|