Ticket #9923: trac_9923.patch

File trac_9923.patch, 10.8 KB (added by ncohen, 9 years ago)
  • sage/graphs/digraph.py

    # HG changeset patch
    # User Nathann Cohen <nathann.cohen@gmail.com>
    # Date 1285887698 -7200
    # Node ID 5c9ed0fbc0efa3956f509f6cfccc6039eeb8ac1a
    # Parent  f30bb70155c938e77cdea1da6c71e688fc3502b5
    trac 9923 -- Minimum Feedback Arc/Edge set through constraint generation
    
    diff -r f30bb70155c9 -r 5c9ed0fbc0ef sage/graphs/digraph.py
    a b  
    10911091        return sorted(self.out_degree_iterator(), reverse=True)
    10921092
    10931093
    1094     def feedback_edge_set(self, value_only=False, solver=None, verbose=0):
     1094    def feedback_edge_set(self, constraint_generation= True, value_only=False, solver=None, verbose=0):
    10951095        r"""
    10961096        Computes the minimum feedback edge set of a digraph
    10971097        (also called feedback arc set).
     
    11141114          - When set to ``False``, the ``Set`` of edges of a minimal edge set
    11151115            is returned.
    11161116
     1117        - ``constraint_generation`` (boolean) -- whether to use
     1118          constraint generation when solving the Mixed Integer Linear
     1119          Program (default: ``True``).
     1120
    11171121        - ``solver`` -- (default: ``None``) Specify a Linear Program (LP)
    11181122          solver to be used. If set to ``None``, the default one is used. For
    11191123          more information on LP solvers and which default solver is used, see
     
    11251129        - ``verbose`` -- integer (default: ``0``). Sets the level of
    11261130          verbosity. Set to 0 by default, which means quiet.
    11271131
    1128         This problem is solved using Linear Programming, which certainly
    1129         is not the best way and will have to be updated. The program to be
    1130         solved is the following:
     1132        ALGORITHM:
     1133
     1134        This problem is solved using Linear Programming, in two
     1135        different ways. The first one is to solve the following
     1136        formulation:
    11311137
    11321138        .. MATH::
    11331139
     
    11511157        The number of edges removed is then minimized, which is
    11521158        the objective.
    11531159
     1160        (Constraint Generation)
     1161 
     1162        If the parameter ``constraint_generation`` is enabled, a more
     1163        efficient formulation is used :
     1164 
     1165        .. MATH::
     1166 
     1167            \mbox{Minimize : }&\sum_{(u,v)\in G} b_{(u,v)}\\ 
     1168            \mbox{Such that : }&\\ 
     1169            &\forall C\text{ circuits }\subseteq\in G, \sum_{uv\in C}b_{(u,v)}\geq 1\\
     1170 
     1171        As the number of circuits contained in a graph is exponential,
     1172        this LP is solved through constraint generation. This means
     1173        that the solver is sequentially asked to solved the problem,
     1174        knowing only a portion of the circuits contained in `G`, each
     1175        time adding to the list of its constraints the circuit which
     1176        its last answer had left intact.
     1177
    11541178        EXAMPLES:
    11551179
    11561180        If the digraph is created from a graph, and hence is symmetric
     
    11741198           sage: (u,v,l) = g.edge_iterator().next()
    11751199           sage: (u,v) in feedback or (v,u) in feedback
    11761200           True
     1201
     1202        TESTS:
     1203 
     1204        Comparing with/without constraint generation::
     1205 
     1206            sage: g = digraphs.RandomDirectedGNP(10,.3)
     1207            sage: x = g.feedback_edge_set(value_only = True)
     1208            sage: y = g.feedback_edge_set(value_only = True, constraint_generation = False)
     1209            sage: x == y
     1210            True
    11771211        """
     1212
     1213        if constraint_generation:
     1214            from sage.graphs.generic_graph_pyx import minimum_feedback_arc_set
     1215            obj,mfas = minimum_feedback_arc_set(self, constraint_log = verbose, solver_log = verbose)
     1216            if value_only:
     1217                return obj
     1218            else:
     1219                return mfas
    11781220       
    11791221        from sage.numerical.mip import MixedIntegerLinearProgram, Sum
    11801222       
  • sage/graphs/generic_graph_pyx.pyx

    diff -r f30bb70155c9 -r 5c9ed0fbc0ef sage/graphs/generic_graph_pyx.pyx
    a b  
    949949        sage_free(u)
    950950        sage_free(v)
    951951
     952##############################
     953# Minimum feedback arc set   #
     954##############################
     955from sage.numerical.backends.generic_backend cimport GenericBackend
     956
     957cpdef 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
     1062cdef 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
    9521157
    9531158cpdef tuple find_hamiltonian( G, long max_iter=100000, long reset_bound=30000, long backtrack_bound=1000, find_path=False ):
    9541159    r"""
     
    10571262    answers are still satisfiable path (the algorithm would raise an
    10581263    exception otherwise)::
    10591264
    1060         sage: for i in range(200):
    1061         ...      g = graphs.RandomGNP(20,.1)
    1062         ...      _ = fh(G,find_path=True)
     1265        sage: for i in range(50):              # long
     1266        ...      g = graphs.RandomGNP(20,.4)   # long
     1267        ...      _ = fh(G,find_path=True)      # long
    10631268    """
    10641269
    10651270    from sage.misc.prandom import randint