| 952 | ############################## |
| 953 | # Minimum feedback arc set # |
| 954 | ############################## |
| 955 | from sage.numerical.backends.generic_backend cimport GenericBackend |
| 956 | |
| 957 | cpdef minimum_feedback_arc_set(g, constraint_log = 0,solver_log = 0): |
| 958 | r""" |
| 959 | Returns a minimum feedback arc set of the given digraph. |
| 960 | |
| 961 | INPUT: |
| 962 | |
| 963 | - ``g`` -- a digraph |
| 964 | |
| 965 | - ``constraint_log`` (boolean) -- whether to prin informations |
| 966 | about the constraints generation. ``0`` by default. |
| 967 | |
| 968 | - ``solver_log`` (integer) -- log lever for GLPK. ``0`` by |
| 969 | default. |
| 970 | |
| 971 | ALGORITHM: |
| 972 | |
| 973 | Mixed Integer Linear Program with constraint generation. |
| 974 | |
| 975 | Each edge has a corresponding binary variable, set to 1 if the |
| 976 | edge is included in the MFAS, and set to 0 otherwise. We want to |
| 977 | minimize their number. |
| 978 | |
| 979 | We first solve the problem without any constraint, which means all |
| 980 | the variables have a value of 0, and that no edges are to be |
| 981 | removed. We are likely to find a circuit in this graph, so we look |
| 982 | for it using the method ``find_cycle``, then add to our LP that |
| 983 | the sum of the variables must be at least 1 over this cycle (it |
| 984 | has to be broken !) |
| 985 | |
| 986 | EXAMPLES: |
| 987 | |
| 988 | Disjoint union of cycles (this method is tested through |
| 989 | ``DiGraph.feedback_edge_set`` too):: |
| 990 | |
| 991 | sage: from sage.graphs.generic_graph_pyx import minimum_feedback_arc_set |
| 992 | sage: g = digraphs.Circuit(5)*10 |
| 993 | sage: minimum_feedback_arc_set(g)[0] |
| 994 | 10.0 |
| 995 | """ |
| 996 | |
| 997 | from sage.numerical.backends.generic_backend import get_solver |
| 998 | |
| 999 | cdef int n = g.order() |
| 1000 | cdef int m = g.size() |
| 1001 | cdef list indices |
| 1002 | cdef list values |
| 1003 | |
| 1004 | # Associating its variable number to each edge of the graph |
| 1005 | cdef dict edge_id = dict( [(x,y) for (y,x) |
| 1006 | in enumerate(g.edges(labels = False))] ) |
| 1007 | |
| 1008 | cdef GenericBackend p |
| 1009 | p = get_solver(constraint_generation = True) |
| 1010 | p.add_variables(m) |
| 1011 | p.set_sense(-1) |
| 1012 | |
| 1013 | # Variables are binary, and their coefficient in the objective is |
| 1014 | # 1 |
| 1015 | for 0 <= i < m: |
| 1016 | p.set_variable_type(i,0) |
| 1017 | |
| 1018 | p.set_objective([1]*m) |
| 1019 | |
| 1020 | p.set_verbosity(solver_log) |
| 1021 | |
| 1022 | p.solve() |
| 1023 | |
| 1024 | # For as long as we do not break because the digraph is |
| 1025 | # acyclic.... |
| 1026 | while (1): |
| 1027 | |
| 1028 | # Find a cycle |
| 1029 | c = find_cycle(p, g, edge_id) |
| 1030 | |
| 1031 | # If there is none, we are done ! |
| 1032 | if not c: |
| 1033 | break |
| 1034 | |
| 1035 | if constraint_log: |
| 1036 | print "Adding a constraint on cycle : ",c |
| 1037 | |
| 1038 | indices = [] |
| 1039 | for 0<= i < len(c)-1: |
| 1040 | indices.append(edge_id[(c[i],c[i+1])]) |
| 1041 | |
| 1042 | indices.append(edge_id[(c[len(c)-1],c[0])]) |
| 1043 | |
| 1044 | values = [1] * len(c) |
| 1045 | |
| 1046 | p.add_linear_constraint(indices, values, -1, 1) |
| 1047 | |
| 1048 | p.solve() |
| 1049 | |
| 1050 | # Getting the objective back |
| 1051 | obj = p.get_objective_value() |
| 1052 | |
| 1053 | # listing the edges contained in the MFAS |
| 1054 | edges = [(u,v) for u,v in g.edges(labels = False) |
| 1055 | if p.get_variable_value(edge_id[(u,v)]) > .5] |
| 1056 | |
| 1057 | # Accomplished ! |
| 1058 | return obj,edges |
| 1059 | |
| 1060 | |
| 1061 | # Find a cycle in the digraph given by the previous solution returned |
| 1062 | cdef list find_cycle(GenericBackend lp, g, dict edges_id): |
| 1063 | r""" |
| 1064 | Returns a cycle if there exists one in the digraph where the edges |
| 1065 | removed by the current LP solution have been removed. |
| 1066 | |
| 1067 | If there is none, an empty list is returned. |
| 1068 | """ |
| 1069 | |
| 1070 | # Vertices that have already been seen |
| 1071 | cdef dict seen = dict( [(v,False) for v in g] ) |
| 1072 | |
| 1073 | # Activated vertices |
| 1074 | cdef dict activated = dict( [(v,True) for v in g] ) |
| 1075 | |
| 1076 | # Vertices whose neighbors have already been added to the stack |
| 1077 | cdef dict tried = dict( [(v,False) for v in g] ) |
| 1078 | |
| 1079 | # Parent of a vertex in the discovery tree |
| 1080 | cdef dict parent = {} |
| 1081 | |
| 1082 | # The vertices left to be visited |
| 1083 | cdef list stack = [] |
| 1084 | |
| 1085 | # We try any vertex as the source of the exploration tree |
| 1086 | for v in g: |
| 1087 | |
| 1088 | # We are not interested in trying de-activated vertices |
| 1089 | if not activated[v]: |
| 1090 | continue |
| 1091 | |
| 1092 | stack = [v] |
| 1093 | seen[v] = True |
| 1094 | |
| 1095 | # For as long as some vertices are to be visited |
| 1096 | while stack: |
| 1097 | |
| 1098 | # We take the last one (depth-first search) |
| 1099 | u = stack.pop(-1) |
| 1100 | |
| 1101 | # If we never tried it, we put it at the end of the stack |
| 1102 | # again, but remember we tried it once. |
| 1103 | if not tried[u]: |
| 1104 | tried[u] = True |
| 1105 | seen[u] = True |
| 1106 | stack.append(u) |
| 1107 | |
| 1108 | # If we tried all the path starting from this vertex |
| 1109 | # already, it means it leads to no circuit. We can remove |
| 1110 | # it, and go to the next one |
| 1111 | else: |
| 1112 | seen[u] = False |
| 1113 | activated[u] = False |
| 1114 | continue |
| 1115 | |
| 1116 | # For each out-neighbor of the current vertex |
| 1117 | for uu in g.neighbors_out(u): |
| 1118 | |
| 1119 | # We are not interested in using the edges removed by |
| 1120 | # the previous solution |
| 1121 | if lp.get_variable_value(edges_id[(u,uu)]) > .5: |
| 1122 | continue |
| 1123 | |
| 1124 | # If we have found a new vertex, we put it at the end |
| 1125 | # of the stack. We ignored de-activated vertices. |
| 1126 | if not seen[uu]: |
| 1127 | if activated[uu]: |
| 1128 | parent[uu] = u |
| 1129 | stack.append(uu) |
| 1130 | |
| 1131 | # If we have already met this vertex, it means we have |
| 1132 | # found a circuit ! |
| 1133 | else: |
| 1134 | |
| 1135 | # We build it, then return it |
| 1136 | answer = [u] |
| 1137 | tmp = u |
| 1138 | while u != uu: |
| 1139 | u = parent.get(u,uu) |
| 1140 | answer.append(u) |
| 1141 | |
| 1142 | answer.reverse() |
| 1143 | return answer |
| 1144 | |
| 1145 | # Each time we have explored all the descendants of a vertex v |
| 1146 | # and met no circuit, we de-activate all these vertices, as we |
| 1147 | # know we do not need them anyway. |
| 1148 | for v in g: |
| 1149 | if seen[v]: |
| 1150 | seen[v] = False |
| 1151 | activated[v] = False |
| 1152 | tried[v] = True |
| 1153 | |
| 1154 | # No Cycle... Good news ! Let's return it. |
| 1155 | return [] |
| 1156 | |