Ticket #8403: trac_8403-rebased.patch

File trac_8403-rebased.patch, 6.2 KB (added by rlm, 11 years ago)

apply before part 2

  • sage/graphs/generic_graph.py

    # HG changeset patch
    # User Nathann Cohen <nathann.cohen@gmail.com>
    # Date 1267379885 -3600
    # Node ID a3c4858b2b0d21f6511dcf499387967ad2f21881
    # Parent  5190a7e83c64077af77be7a41425523c622592f8
    trac 8403 : Graph.steiner_tree
    
    diff -r 5190a7e83c64 -r a3c4858b2b0d sage/graphs/generic_graph.py
    a b  
    32193219            v = pv
    32203220        return B, C
    32213221
     3222
     3223    def steiner_tree(self,vertices, weighted = False):
     3224        r"""
     3225        Returns a tree of minimum weight connecting the given
     3226        set of vertices.
     3227
     3228        Definition :
     3229
     3230        Computing a minimum spanning tree in a graph can be done in `n
     3231        \log(n)` time (and in linear time if all weights are
     3232        equal). On the other hand, if one is given a large (possibly
     3233        weighted) graph and a subset of its vertices, it is NP-Hard to
     3234        find a tree of minimum weight connecting the given set of
     3235        vertices, which is then called a Steiner Tree.
     3236       
     3237        `Wikipedia article on Steiner Trees
     3238        <http://en.wikipedia.org/wiki/Steiner_tree_problem>`_.
     3239
     3240        INPUT:
     3241
     3242        - ``vertices`` -- the vertices to be connected by the Steiner
     3243          Tree.
     3244         
     3245        - ``weighted`` (boolean) -- Whether to consider the graph as
     3246          weighted, and use each edge's label as a weight, considering
     3247          ``None`` as a weight of `1`. If ``weighted=False`` (default)
     3248          all edges are considered to have a weight of `1`.
     3249
     3250        .. NOTE::
     3251
     3252            * This problem being defined on undirected graphs, the
     3253              orientation is not considered if the current graph is
     3254              actually a digraph.
     3255
     3256            * The graph is assumed not to have multiple edges.
     3257
     3258        ALGORITHM:
     3259
     3260        Solved through Linear Programming.
     3261
     3262        COMPLEXITY:
     3263
     3264        NP-Hard.
     3265
     3266        Note that this algorithm first checks whether the given
     3267        set of vertices induces a connected graph, returning one of its
     3268        spanning trees if ``weighted`` is set to ``False``, and thus
     3269        answering very quickly in some cases
     3270       
     3271        EXAMPLES:
     3272
     3273        The Steiner Tree of the first 5 vertices in a random graph is,
     3274        of course, always a tree ::
     3275
     3276            sage: g = graphs.RandomGNP(30,.5)
     3277            sage: st = g.steiner_tree(g.vertices()[:5])              # optional - requires GLPK, CBC or CPLEX
     3278            sage: st.is_tree()                                       # optional - requires GLPK, CBC or CPLEX
     3279            True
     3280
     3281        And all the 5 vertices are contained in this tree ::
     3282
     3283            sage: all([v in st for v in g.vertices()[:5] ])          # optional - requires GLPK, CBC or CPLEX
     3284            True
     3285
     3286        An exception is raised when the problem is impossible, i.e.
     3287        if the given vertices are not all included in the same
     3288        connected component ::
     3289
     3290            sage: g = 2 * graphs.PetersenGraph()
     3291            sage: st = g.steiner_tree([5,15])
     3292            Traceback (most recent call last):
     3293            ...
     3294            ValueError: The given vertices do not all belong to the same connected component. This problem has no solution !
     3295        """
     3296
     3297        if self.is_directed():
     3298            g = Graph(self)
     3299        else:
     3300            g = self
     3301
     3302        if g.has_multiple_edges():
     3303            raise ValueError("The graph is expected not to have multiple edges.")
     3304
     3305        # Can the problem be solved ? Are all the vertices in the same
     3306        # connected component ?
     3307        cc = g.connected_component_containing_vertex(vertices[0])
     3308        if not all([v in cc for v in vertices]):
     3309            raise ValueError("The given vertices do not all belong to the same connected component. This problem has no solution !")
     3310
     3311        # Can it be solved using the min spanning tree algorithm ?
     3312        if not weighted:
     3313            gg = g.subgraph(vertices)
     3314            if gg.is_connected():
     3315                st = g.subgraph(edges = gg.min_spanning_tree())
     3316                st.delete_vertices([v for v in g if st.degree(v) == 0])
     3317                return st
     3318
     3319        # Then, LP formulation
     3320        from sage.numerical.mip import MixedIntegerLinearProgram
     3321        p = MixedIntegerLinearProgram(maximization = False)
     3322
     3323        # Reorder an edge
     3324        R = lambda (x,y) : (x,y) if x<y else (y,x)
     3325       
     3326        # edges used in the Steiner Tree
     3327        edges = p.new_variable()
     3328
     3329        # relaxed edges to test for acyclicity
     3330        r_edges = p.new_variable()
     3331
     3332        # Whether a vertex is in the Steiner Tree
     3333        vertex = p.new_variable()
     3334        for v in g:
     3335            for e in g.edges_incident(v, labels=False):
     3336                p.add_constraint(vertex[v] - edges[R(e)], min = 0)
     3337
     3338        # We must have the given vertices in our tree
     3339        for v in vertices:
     3340            p.add_constraint(sum([edges[R(e)] for e in g.edges_incident(v,labels=False)]), min=1)
     3341
     3342        # The number of edges is equal to the number of vertices in our tree minus 1
     3343        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)
     3344
     3345        # There are no cycles in our graph
     3346
     3347        for u,v in g.edges(labels = False):
     3348            p.add_constraint( r_edges[(u,v)]+ r_edges[(v,u)] - edges[R((u,v))] , min = 0 )
     3349
     3350        eps = 1/(5*Integer(g.order()))
     3351        for v in g:
     3352            p.add_constraint(sum([r_edges[(u,v)] for u in g.neighbors(v)]), max = 1-eps)
     3353               
     3354
     3355        # Objective
     3356        if weighted:
     3357            w = lambda (x,y) : g.edge_label(x,y) if g.edge_label(x,y) is not None else 1
     3358        else:
     3359            w = lambda (x,y) : 1
     3360
     3361        p.set_objective(sum([w(e)*edges[R(e)] for e in g.edges(labels = False)]))
     3362
     3363        p.set_binary(edges)           
     3364        p.solve()
     3365
     3366        edges = p.get_values(edges)
     3367
     3368        st =  g.subgraph(edges=[e for e in g.edges(labels = False) if edges[R(e)] == 1])
     3369        st.delete_vertices([v for v in g if st.degree(v) == 0])
     3370        return st
     3371
    32223372    def edge_disjoint_spanning_trees(self,k, root=None, **kwds):
    32233373        r"""
    32243374        Returns the desired number of edge-disjoint spanning