| 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 |