Ticket #8403: trac_8403.patch

File trac_8403.patch, 6.2 KB (added by ncohen, 11 years ago)
  • sage/graphs/generic_graph.py

    # HG changeset patch
    # User Nathann Cohen <nathann.cohen@gmail.com>
    # Date 1267379885 -3600
    # Node ID 1a774059334288a72aaf5e897fb146eb56e41016
    # Parent  556bb66e4c6dbb92a4ee37c1750d82a5c6298eeb
    trac 8403 : Graph.steiner_tree
    
    diff -r 556bb66e4c6d -r 1a7740593342 sage/graphs/generic_graph.py
    a b  
    30863086
    30873087            v = pv
    30883088        return B, C
     3089
     3090
     3091    def steiner_tree(self,vertices, weighted = False):
     3092        r"""
     3093        Returns a tree of minimum weight connecting the given
     3094        set of vertices.
     3095
     3096        Definition :
     3097
     3098        Computing a minimum spanning tree in a graph can be done in `n
     3099        \log(n)` time (and in linear time if all weights are
     3100        equal). On the other hand, if one is given a large (possibly
     3101        weighted) graph and a subset of its vertices, it is NP-Hard to
     3102        find a tree of minimum weight connecting the given set of
     3103        vertices, which is then called a Steiner Tree.
     3104       
     3105        `Wikipedia article on Steiner Trees
     3106        <http://en.wikipedia.org/wiki/Steiner_tree_problem>`_.
     3107
     3108        INPUT:
     3109
     3110        - ``vertices`` -- the vertices to be connected by the Steiner
     3111          Tree.
     3112         
     3113        - ``weighted`` (boolean) -- Whether to consider the graph as
     3114          weighted, and use each edge's label as a weight, considering
     3115          ``None`` as a weight of `1`. If ``weighted=False`` (default)
     3116          all edges are considered to have a weight of `1`.
     3117
     3118        .. NOTE::
     3119
     3120            * This problem being defined on undirected graphs, the
     3121              orientation is not considered if the current graph is
     3122              actually a digraph.
     3123
     3124            * The graph is assumed not to have multiple edges.
     3125
     3126        ALGORITHM:
     3127
     3128        Solved through Linear Programming.
     3129
     3130        COMPLEXITY:
     3131
     3132        NP-Hard.
     3133
     3134        Note that this algorithm first checks whether the given
     3135        set of vertices induces a connected graph, returning one of its
     3136        spanning trees if ``weighted`` is set to ``False``, and thus
     3137        answering very quickly in some cases
     3138       
     3139        EXAMPLES:
     3140
     3141        The Steiner Tree of the first 5 vertices in a random graph is,
     3142        of course, always a tree ::
     3143
     3144            sage: g = graphs.RandomGNP(30,.5)
     3145            sage: st = g.steiner_tree(g.vertices()[:5])              # optional - requires GLPK, CBC or CPLEX
     3146            sage: st.is_tree()                                       # optional - requires GLPK, CBC or CPLEX
     3147            True
     3148
     3149        And all the 5 vertices are contained in this tree ::
     3150
     3151            sage: all([v in st for v in g.vertices()[:5] ])          # optional - requires GLPK, CBC or CPLEX
     3152            True
     3153
     3154        An exception is raised when the problem is impossible, i.e.
     3155        if the given vertices are not all included in the same
     3156        connected component ::
     3157
     3158            sage: g = 2 * graphs.PetersenGraph()
     3159            sage: st = g.steiner_tree([5,15])
     3160            Traceback (most recent call last):
     3161            ...
     3162            ValueError: The given vertices do not all belong to the same connected component. This problem has no solution !
     3163        """
     3164
     3165        if self.is_directed():
     3166            g = Graph(self)
     3167        else:
     3168            g = self
     3169
     3170        if g.has_multiple_edges():
     3171            raise ValueError("The graph is expected not to have multiple edges.")
     3172
     3173        # Can the problem be solved ? Are all the vertices in the same
     3174        # connected component ?
     3175        cc = g.connected_component_containing_vertex(vertices[0])
     3176        if not all([v in cc for v in vertices]):
     3177            raise ValueError("The given vertices do not all belong to the same connected component. This problem has no solution !")
     3178
     3179        # Can it be solved using the min spanning tree algorithm ?
     3180        if not weighted:
     3181            gg = g.subgraph(vertices)
     3182            if gg.is_connected():
     3183                st = g.subgraph(edges = gg.min_spanning_tree())
     3184                st.delete_vertices([v for v in g if st.degree(v) == 0])
     3185                return st
     3186
     3187        # Then, LP formulation
     3188        from sage.numerical.mip import MixedIntegerLinearProgram
     3189        p = MixedIntegerLinearProgram(maximization = False)
     3190
     3191        # Reorder an edge
     3192        R = lambda (x,y) : (x,y) if x<y else (y,x)
     3193       
     3194        # edges used in the Steiner Tree
     3195        edges = p.new_variable()
     3196
     3197        # relaxed edges to test for acyclicity
     3198        r_edges = p.new_variable()
     3199
     3200        # Whether a vertex is in the Steiner Tree
     3201        vertex = p.new_variable()
     3202        for v in g:
     3203            for e in g.edges_incident(v, labels=False):
     3204                p.add_constraint(vertex[v] - edges[R(e)], min = 0)
     3205
     3206        # We must have the given vertices in our tree
     3207        for v in vertices:
     3208            p.add_constraint(sum([edges[R(e)] for e in g.edges_incident(v,labels=False)]), min=1)
     3209
     3210        # The number of edges is equal to the number of vertices in our tree minus 1
     3211        p.add_constraint(sum([vertex[v] for v in g]) - sum([edges[R(e)] for e in g.edges(labels=None)]), max = 1, min = 1)
     3212
     3213        # There are no cycles in our graph
     3214
     3215        for u,v in g.edges(labels = False):
     3216            p.add_constraint( r_edges[(u,v)]+ r_edges[(v,u)] - edges[R((u,v))] , min = 0 )
     3217
     3218        eps = 1/(5*Integer(g.order()))
     3219        for v in g:
     3220            p.add_constraint(sum([r_edges[(u,v)] for u in g.neighbors(v)]), max = 1-eps)
     3221               
     3222
     3223        # Objective
     3224        if weighted:
     3225            w = lambda (x,y) : g.edge_label(x,y) if g.edge_label(x,y) is not None else 1
     3226        else:
     3227            w = lambda (x,y) : 1
     3228
     3229        p.set_objective(sum([w(e)*edges[R(e)] for e in g.edges(labels = False)]))
     3230
     3231        p.set_binary(edges)           
     3232        p.solve()
     3233
     3234        edges = p.get_values(edges)
     3235
     3236        st =  g.subgraph(edges=[e for e in g.edges(labels = False) if edges[R(e)] == 1])
     3237        st.delete_vertices([v for v in g if st.degree(v) == 0])
     3238        return st
    30893239               
    30903240    def edge_cut(self, s, t, value_only=True, use_edge_labels=False, vertices=False, solver=None, verbose=0):
    30913241        r"""