python - Chu-Liu Edmond's algorithm for Minimum Spanning Tree on Directed Graphs -

i find minimum spanning tree (mst) on weighted directed graph. have been trying use chu-liu/edmond's algorithm, have implemented in python (code below). simple, clear description of algorithm can found here. have 2 questions.

1. is edmond's algorithm guaranteed converge on solution?

   i concerned removing cycle add cycle. if happens, algorithm continue trying remove cycles forever.

   i seem have found example happens. input graph shown below (in code). algorithm never finishes because switches between cycles [1,2] , [1,3], , [5,4] , [5,6]. edge added graph resolve cycle [5,4] creates cycle [5,6] , vice versa, , [1,2] , [1,3].

   i should note not implementation correct.

2. to resolve issue, introduced ad hoc patch. when edge removed remove cycle, permanently remove edge underlying graph g on searching mst. consequently, edge cannot added again , should prevent algorithm getting stuck. change, guaranteed find mst?

   i suspect 1 can find pathological case step lead result not mst, have not been able think of one. seems work on simple test cases have tried.

code:

import sys  # --------------------------------------------------------------------------------- #  def _reverse(graph):     r = {}     src in graph:         (dst,c) in graph[src].items():             if dst in r:                 r[dst][src] = c             else:                 r[dst] = { src : c }     return r  # finds cycles in graph using tarjan's algorithm def strongly_connected_components(graph):     """     tarjan's algorithm (named discoverer, robert tarjan) graph theory algorithm     finding connected components of graph.      based on: http://en.wikipedia.org/wiki/tarjan%27s_strongly_connected_components_algorithm     """      index_counter = [0]     stack = []     lowlinks = {}     index = {}     result = []      def strongconnect(node):         # set depth index node smallest unused index         index[node] = index_counter[0]         lowlinks[node] = index_counter[0]         index_counter[0] += 1         stack.append(node)          # consider successors of `node`         try:             successors = graph[node]         except:             successors = []         successor in successors:             if successor not in lowlinks:                 # successor has not yet been visited; recurse on                 strongconnect(successor)                 lowlinks[node] = min(lowlinks[node],lowlinks[successor])             elif successor in stack:                 # successor in stack , hence in current connected component (scc)                 lowlinks[node] = min(lowlinks[node],index[successor])          # if `node` root node, pop stack , generate scc         if lowlinks[node] == index[node]:             connected_component = []              while true:                 successor = stack.pop()                 connected_component.append(successor)                 if successor == node: break             component = tuple(connected_component)             # storing result             result.append(component)      node in graph:         if node not in lowlinks:             strongconnect(node)      return result  def _mergecycles(cycle,g,rg,g,rg):     allinedges = [] # edges entering cycle outside cycle     mininternal = none     mininternalweight = sys.maxint      # find minimal internal edge weight     n in cycle:         e in rg[n]:             if e in cycle:                 if mininternal none or rg[n][e] < mininternalweight:                     mininternal = (n,e)                     mininternalweight = rg[n][e]                     continue             else:                 allinedges.append((n,e)) # edge enters cycle      # find incoming edge minimum modified cost     # modified cost c(i,k) = c(i,j) - (c(x_j, j) - min{j}(c(x_j, j)))     minexternal = none     minmodifiedweight = 0     j,i in allinedges: # j vertex in cycle, candidate vertex outside cycle         xj, weight_xj_j = rg[j].popitem() # xj vertex in cycle goes j         rg[j][xj] = weight_xj_j # put item in dictionary         w = rg[j][i] - (weight_xj_j - mininternalweight) # c(i,k) = c(i,j) - (c(x_j, j) - min{j}(c(x_j, j)))         if minexternal none or w <= minmodifiedweight:             minexternal = (j,i)             minmodifiedweight = w      w = rg[minexternal[0]][minexternal[1]] # weight of edge entering cycle     xj,_ = rg[minexternal[0]].popitem() # xj vertex in cycle goes j     rem = (minexternal[0], xj) # edge remove     rg[minexternal[0]].clear() # popitem() should delete 1 edge j, ensure      # remove offending edge rg     # rg[minexternal[0]].pop(xj, none) #highly experimental. throw away offending edge, never again      if rem[1] in g:         if rem[0] in g[rem[1]]:             del g[rem[1]][rem[0]]     if minexternal[1] in g:         g[minexternal[1]][minexternal[0]] = w     else:         g[minexternal[1]] = { minexternal[0] : w }      rg = _reverse(g)  # --------------------------------------------------------------------------------- #  def mst(root,g):     """ chu-liu/edmond's algorithm      arguments:      root - root of mst     g - graph in mst lies      returns: graph representation of mst      graph representation same 1 found at:     http://code.activestate.com/recipes/119466/      explanation copied verbatim here:      input graph g assumed have following     representation: vertex can object can     used index dictionary.  g     dictionary, indexed vertices.  vertex v,     g[v] dictionary, indexed neighbors     of v.  edge v->w, g[v][w] length of     edge.     """      rg = _reverse(g)      g = {}     n in rg:         if len(rg[n]) == 0:             continue         minimum = sys.maxint         s,d = none,none          e in rg[n]:             if rg[n][e] < minimum:                 minimum = rg[n][e]                 s,d = n,e          if d in g:             g[d][s] = rg[s][d]         else:             g[d] = { s : rg[s][d] }      cycles = [list(c) c in strongly_connected_components(g)]      cycles_exist = true     while cycles_exist:          cycles_exist = false         cycles = [list(c) c in strongly_connected_components(g)]         rg = _reverse(g)          cycle in cycles:              if root in cycle:                 continue              if len(cycle) == 1:                 continue              _mergecycles(cycle, g, rg, g, rg)             cycles_exist = true      return g  # --------------------------------------------------------------------------------- #  if __name__ == "__main__":      # example of input works     root = 0     g = {0: {1: 23, 2: 22, 3: 22}, 1: {2: 1, 3: 1}, 3: {1: 1, 2: 0}}      # example of input causes infinite cycle     root = 0     g = {0: {1: 17, 2: 16, 3: 19, 4: 16, 5: 16, 6: 18}, 1: {2: 3, 3: 3, 4: 11, 5: 10, 6: 12}, 2: {1: 3, 3: 4, 4: 8, 5: 8, 6: 11}, 3: {1: 3, 2: 4, 4: 12, 5: 11, 6: 14}, 4: {1: 11, 2: 8, 3: 12, 5: 6, 6: 10}, 5: {1: 10, 2: 8, 3: 11, 4: 6, 6: 4}, 6: {1: 12, 2: 11, 3: 14, 4: 10, 5: 4}}      h = mst(int(root),g)      print h      s in h:         t in h[s]:             print "%d-%d" % (s,t) 

don't ad hoc patches. concede implementing contraction/uncontraction logic not intuitive, , recursion undesirable in contexts, here's proper python implementation made production quality. rather perform uncontraction step @ each recursive level, defer end , use depth-first search, thereby avoiding recursion. (the correctness of modification follows complementary slackness, part of theory of linear programming.)

the naming convention below _rep signifies supernode (i.e., block of 1 or more contracted nodes).

#!/usr/bin/env python3 collections import defaultdict, namedtuple   arc = namedtuple('arc', ('tail', 'weight', 'head'))   def min_spanning_arborescence(arcs, sink):     good_arcs = []     quotient_map = {arc.tail: arc.tail arc in arcs}     quotient_map[sink] = sink     while true:         min_arc_by_tail_rep = {}         successor_rep = {}         arc in arcs:             if arc.tail == sink:                 continue             tail_rep = quotient_map[arc.tail]             head_rep = quotient_map[arc.head]             if tail_rep == head_rep:                 continue             if tail_rep not in min_arc_by_tail_rep or min_arc_by_tail_rep[tail_rep].weight > arc.weight:                 min_arc_by_tail_rep[tail_rep] = arc                 successor_rep[tail_rep] = head_rep         cycle_reps = find_cycle(successor_rep, sink)         if cycle_reps none:             good_arcs.extend(min_arc_by_tail_rep.values())             return spanning_arborescence(good_arcs, sink)         good_arcs.extend(min_arc_by_tail_rep[cycle_rep] cycle_rep in cycle_reps)         cycle_rep_set = set(cycle_reps)         cycle_rep = cycle_rep_set.pop()         quotient_map = {node: cycle_rep if node_rep in cycle_rep_set else node_rep node, node_rep in quotient_map.items()}   def find_cycle(successor, sink):     visited = {sink}     node in successor:         cycle = []         while node not in visited:             visited.add(node)             cycle.append(node)             node = successor[node]         if node in cycle:             return cycle[cycle.index(node):]     return none   def spanning_arborescence(arcs, sink):     arcs_by_head = defaultdict(list)     arc in arcs:         if arc.tail == sink:             continue         arcs_by_head[arc.head].append(arc)     solution_arc_by_tail = {}     stack = arcs_by_head[sink]     while stack:         arc = stack.pop()         if arc.tail in solution_arc_by_tail:             continue         solution_arc_by_tail[arc.tail] = arc         stack.extend(arcs_by_head[arc.tail])     return solution_arc_by_tail   print(min_spanning_arborescence([arc(1, 17, 0), arc(2, 16, 0), arc(3, 19, 0), arc(4, 16, 0), arc(5, 16, 0), arc(6, 18, 0), arc(2, 3, 1), arc(3, 3, 1), arc(4, 11, 1), arc(5, 10, 1), arc(6, 12, 1), arc(1, 3, 2), arc(3, 4, 2), arc(4, 8, 2), arc(5, 8, 2), arc(6, 11, 2), arc(1, 3, 3), arc(2, 4, 3), arc(4, 12, 3), arc(5, 11, 3), arc(6, 14, 3), arc(1, 11, 4), arc(2, 8, 4), arc(3, 12, 4), arc(5, 6, 4), arc(6, 10, 4), arc(1, 10, 5), arc(2, 8, 5), arc(3, 11, 5), arc(4, 6, 5), arc(6, 4, 5), arc(1, 12, 6), arc(2, 11, 6), arc(3, 14, 6), arc(4, 10, 6), arc(5, 4, 6)], 0))