4081 | | flow_sum=lambda X: sum([flow[X][v] for (u,v) in g.outgoing_edges([X],labels=None)])-sum([flow[u][X] for (u,v) in g.incoming_edges([X],labels=None)]) |
4082 | | |
4083 | | # Maximizes the flow leaving x |
4084 | | p.set_objective(flow_sum(x)) |
4085 | | |
4086 | | # Elsewhere, the flow is equal to 0- |
| 4079 | flow_sum=lambda X: sum([flow[(X,v)] for (u,v) in g.outgoing_edges([X],labels=None)])-sum([flow[(u,X)] for (u,v) in g.incoming_edges([X],labels=None)]) |
| 4080 | |
| 4081 | # The flow leaving x |
| 4082 | flow_leaving = lambda X : sum([flow[(uu,vv)] for (uu,vv) in g.outgoing_edges([X],labels=None)]) |
| 4083 | |
| 4084 | # The flow to be considered when defining the capacity contraints |
| 4085 | capacity_sum = lambda u,v : flow[(u,v)] |
| 4086 | |
| 4087 | else: |
| 4088 | # This function return the balance of flow at X |
| 4089 | flow_sum=lambda X:sum([flow[(X,v)]-flow[(v,X)] for v in g[X]]) |
| 4090 | |
| 4091 | # The flow leaving x |
| 4092 | flow_leaving = lambda X : sum([flow[(X,vv)] for vv in g[X]]) |
| 4093 | |
| 4094 | # The flow to be considered when defining the capacity contraints |
| 4095 | capacity_sum = lambda u,v : flow[(u,v)] + flow[(v,u)] |
| 4096 | |
| 4097 | # Maximizes the flow leaving x |
| 4098 | p.set_objective(flow_sum(x)) |
| 4099 | |
| 4100 | # Elsewhere, the flow is equal to 0 |
| 4101 | for v in g: |
| 4102 | if v!=x and v!=y: |
| 4103 | p.add_constraint(flow_sum(v),min=0,max=0) |
| 4104 | |
| 4105 | # Capacity constraints |
| 4106 | for (u,v,w) in g.edges(): |
| 4107 | p.add_constraint(capacity_sum(u,v),max=capacity(w)) |
| 4108 | |
| 4109 | # No vertex except the sources can send more than 1 |
| 4110 | if vertex_bound: |
4089 | | p.add_constraint(flow_sum(v),min=0,max=0) |
4090 | | |
4091 | | # Capacity constraints |
4092 | | for (u,v,w) in g.edges(): |
4093 | | p.add_constraint(flow[u][v],max=capacity(w)) |
4094 | | |
4095 | | if vertex_bound: |
4096 | | for v in g.vertices(): |
4097 | | if v!=x and v!=y: |
4098 | | p.add_constraint(sum([flow[uu][vv] for (uu,vv) in g.outgoing_edges([v],labels=None)]),max=1) |
4099 | | |
4100 | | else: |
4101 | | # This function return the balance of flow at X |
4102 | | flow_sum=lambda X:sum([flow[X][v]-flow[v][X] for v in g[X]]) |
4103 | | |
4104 | | # Maximizes the flow leaving x |
4105 | | p.set_objective(flow_sum(x)) |
4106 | | |
4107 | | # Elsewhere, the flow is equal to 0 |
4108 | | for v in g: |
4109 | | if v!=x and v!=y: |
4110 | | p.add_constraint(flow_sum(v),min=0,max=0) |
4111 | | |
4112 | | # Capacity constraints |
4113 | | for (u,v,w) in g.edges(): |
4114 | | p.add_constraint(flow[u][v]+flow[v][u],max=capacity(w)) |
4115 | | |
4116 | | if vertex_bound: |
4117 | | for v in g: |
4118 | | if v!=x and v!=y: |
4119 | | p.add_constraint([flow[X][v] for X in g[v]],max=1) |
4120 | | |
| 4113 | p.add_constraint(flow_leaving(v),max=1) |
4131 | | |
4132 | | flow_graph = g.copy() |
4133 | | |
| 4125 | # Builds a clean flow Draph |
| 4126 | flow_graph = g._build_flow_graph(flow, integer=integer) |
| 4127 | |
| 4128 | # Which could be a Graph |
| 4129 | if not self.is_directed(): |
| 4130 | from sage.graphs.graph import Graph |
| 4131 | flow_graph = Graph(flow_graph) |
| 4132 | |
| 4133 | return [obj,flow_graph] |
| 4134 | |
| 4135 | def multicommodity_flow(self, terminals, integer=True, use_edge_labels=False,vertex_bound=False, solver=None, verbose=0): |
| 4136 | r""" |
| 4137 | Solves a multicommodity flow problem. |
| 4138 | |
| 4139 | In the multicommodity flow problem, we are given a set of pairs |
| 4140 | `(s_i, t_i)`, called terminals meaning that `s_i` is willing |
| 4141 | some flow to `t_i`. |
| 4142 | |
| 4143 | Even though it is a natural generalisation of the flow problem |
| 4144 | this version of it is NP-Complete to solve when the flows |
| 4145 | are required to be integer. |
| 4146 | |
| 4147 | For more information, see the |
| 4148 | `Wikipedia page on multicommodity flows |
| 4149 | <http://en.wikipedia.org/wiki/Multi-commodity_flow_problem>`. |
| 4150 | |
| 4151 | INPUT: |
| 4152 | |
| 4153 | - ``terminals`` -- a list of pairs `(s_i, t_i)` or triples |
| 4154 | `(s_i, t_i, w_i)` representing a flow from `s_i` to `t_i` |
| 4155 | of intensity `w_i`. When the pairs are of size `2`, a intensity |
| 4156 | of `1` is assumed. |
| 4157 | |
| 4158 | - ``integer`` (boolean) -- whether to require an integer multicommodity |
| 4159 | flow |
| 4160 | |
| 4161 | - ``use_edge_labels`` (boolean) -- whether to consider the label of edges |
| 4162 | as numerical values representing a capacity. If set to ``False``, a capacity |
| 4163 | of `1` is assumed |
| 4164 | |
| 4165 | - ``vertex_bound`` (boolean) -- whether to require that a vertex can stand at most |
| 4166 | `1` commodity of flow through it of intensity `1`. Terminals can obviously |
| 4167 | still send or receive several units of flow even though vertex_bound is set |
| 4168 | to ``True``, as this parameter is meant to represent topological properties. |
| 4169 | |
| 4170 | - ``solver`` -- Specify a Linear Program solver to be used. |
| 4171 | If set to ``None``, the default one is used. |
| 4172 | function of ``MixedIntegerLinearProgram``. See the documentation of ``MixedIntegerLinearProgram.solve`` |
| 4173 | for more informations. |
| 4174 | |
| 4175 | - ``verbose`` (integer) -- sets the level of verbosity. Set to 0 |
| 4176 | by default (quiet). |
| 4177 | |
| 4178 | ALGORITHM: |
| 4179 | |
| 4180 | (Mixed Integer) Linear Program, depending on the value of ``integer``. |
| 4181 | |
| 4182 | EXAMPLE: |
| 4183 | |
| 4184 | An easy way to obtain a satisfiable multiflow is to compute |
| 4185 | a matching in a graph, and to consider the paired vertices |
| 4186 | as terminals :: |
| 4187 | |
| 4188 | sage: g = graphs.PetersenGraph() |
| 4189 | sage: matching = [(u,v) for u,v,_ in g.matching()] # optional - GLPK, CBC |
| 4190 | sage: h = g.multicommodity_flow(matching) |
| 4191 | sage: len(h) |
| 4192 | 5 |
| 4193 | |
| 4194 | We could also have considered ``g`` as symmetric and computed |
| 4195 | the multiflow in this version instead. In this case, however |
| 4196 | edges can be used in both directions at the same time:: |
| 4197 | |
| 4198 | sage: h = DiGraph(g).multicommodity_flow(matching) # optional - GLPK, CBC |
| 4199 | sage: len(h) |
| 4200 | 5 |
| 4201 | |
| 4202 | An exception is raised when the problem has no solution :: |
| 4203 | |
| 4204 | sage: h = g.multicommodity_flow([(u,v,3) for u,v in matching]) # optional - GLPK, CBC |
| 4205 | Traceback (most recent call last): |
| 4206 | ... |
| 4207 | ValueError: The multiflow problem has no solution |
| 4208 | """ |
| 4209 | |
| 4210 | from sage.numerical.mip import MixedIntegerLinearProgram |
| 4211 | g=self |
| 4212 | p=MixedIntegerLinearProgram(maximization=True) |
| 4213 | |
| 4214 | # Adding the intensity if not present |
| 4215 | terminals = [(x if len(x) == 3 else (x[0],x[1],1)) for x in terminals] |
| 4216 | |
| 4217 | # defining the set of terminals |
| 4218 | set_terminals = set([]) |
| 4219 | for s,t,_ in terminals: |
| 4220 | set_terminals.add(s) |
| 4221 | set_terminals.add(t) |
| 4222 | |
| 4223 | # flow[i][(u,v)] is the flow of commodity i going from u to v |
| 4224 | flow=p.new_variable(dim=2) |
| 4225 | |
| 4226 | # Whether to use edge labels |
| 4227 | if use_edge_labels: |
| 4228 | from sage.rings.real_mpfr import RR |
| 4229 | capacity=lambda x: x if x in RR else 1 |
| 4230 | else: |
| 4231 | capacity=lambda x: 1 |
| 4232 | |
4135 | | for (u,v) in g.edges(labels=None): |
4136 | | # We do not want to see both edges (u,v) and (v,u) |
4137 | | # with a positive flow |
4138 | | if g.has_edge(v,u): |
4139 | | m=min(flow[u][v],flow[v][u]) |
4140 | | flow[u][v]-=m |
4141 | | flow[v][u]-=m |
4142 | | |
4143 | | for (u,v) in g.edge_iterator(labels=None): |
4144 | | if flow[u][v]>0: |
4145 | | flow_graph.set_edge_label(u,v,flow[u][v]) |
| 4234 | # This function return the balance of flow at X |
| 4235 | flow_sum=lambda i,X: sum([flow[i][(X,v)] for (u,v) in g.outgoing_edges([X],labels=None)])-sum([flow[i][(u,X)] for (u,v) in g.incoming_edges([X],labels=None)]) |
| 4236 | |
| 4237 | # The flow leaving x |
| 4238 | flow_leaving = lambda i,X : sum([flow[i][(uu,vv)] for (uu,vv) in g.outgoing_edges([X],labels=None)]) |
| 4239 | |
| 4240 | # the flow to consider when defining the capacity contraints |
| 4241 | capacity_sum = lambda i,u,v : flow[i][(u,v)] |
| 4242 | |
| 4243 | else: |
| 4244 | # This function return the balance of flow at X |
| 4245 | flow_sum=lambda i,X:sum([flow[i][(X,v)]-flow[i][(v,X)] for v in g[X]]) |
| 4246 | |
| 4247 | # The flow leaving x |
| 4248 | flow_leaving = lambda i, X : sum([flow[i][(X,vv)] for vv in g[X]]) |
| 4249 | |
| 4250 | # the flow to consider when defining the capacity contraints |
| 4251 | capacity_sum = lambda i,u,v : flow[i][(u,v)] + flow[i][(v,u)] |
| 4252 | |
| 4253 | |
| 4254 | # Flow constraints |
| 4255 | for i,(s,t,l) in enumerate(terminals): |
| 4256 | for v in g: |
| 4257 | if v == s: |
| 4258 | p.add_constraint(flow_sum(i,v),min=l,max=l) |
| 4259 | elif v == t: |
| 4260 | p.add_constraint(flow_sum(i,v),min=-l,max=-l) |
4147 | | flow_graph.delete_edge(u,v) |
4148 | | |
4149 | | else: |
4150 | | for (u,v) in g.edges(labels=None): |
4151 | | m=min(flow[u][v],flow[v][u]) |
4152 | | flow[u][v]-=m |
4153 | | flow[v][u]-=m |
4154 | | |
4155 | | # We do not want to see both edges (u,v) and (v,u) |
4156 | | # with a positive flow |
4157 | | for (u,v) in g.edges(labels=None): |
4158 | | if flow[u][v]>0: |
4159 | | flow_graph.set_edge_label(u,v,flow[u][v]) |
4160 | | elif flow[v][u]>0: |
4161 | | flow_graph.set_edge_label(v,u,flow[v][u]) |
| 4262 | p.add_constraint(flow_sum(i,v),min=0,max=0) |
| 4263 | |
| 4264 | |
| 4265 | # Capacity constraints |
| 4266 | for (u,v,w) in g.edges(): |
| 4267 | p.add_constraint(sum([capacity_sum(i,u,v) for i in range(len(terminals))]),max=capacity(w)) |
| 4268 | |
| 4269 | |
| 4270 | if vertex_bound: |
| 4271 | |
| 4272 | # Any vertex |
| 4273 | for v in g.vertices(): |
| 4274 | |
| 4275 | # which is an endpoint |
| 4276 | if v in set_terminals: |
| 4277 | for i,(s,t,_) in enumerate(terminals): |
| 4278 | |
| 4279 | # only tolerates the commodities of which it is an endpoint |
| 4280 | if not (v==s or v==t): |
| 4281 | p.add_constraint(flow_leaving(i,v), max = 0) |
| 4282 | |
| 4283 | # which is not an endpoint |
4163 | | flow_graph.delete_edge(v,u) |
4164 | | |
4165 | | return [obj,flow_graph] |
| 4285 | # can stand at most 1 unit of flow through itself |
| 4286 | p.add_constraint(sum([flow_leaving(i,v) for i in range(len(terminals))]), max = 1) |
| 4287 | |
| 4288 | p.set_objective(None) |
| 4289 | |
| 4290 | if integer: |
| 4291 | p.set_integer(flow) |
| 4292 | |
| 4293 | from sage.numerical.mip import MIPSolverException |
| 4294 | |
| 4295 | try: |
| 4296 | obj=p.solve() |
| 4297 | except MIPSolverException: |
| 4298 | raise ValueError("The multiflow problem has no solution") |
| 4299 | |
| 4300 | flow=p.get_values(flow) |
| 4301 | |
| 4302 | # building clean flow digraphs |
| 4303 | flow_graphs = [g._build_flow_graph(flow[i], integer=integer) for i in range(len(terminals))] |
| 4304 | |
| 4305 | # which could be .. graphs ! |
| 4306 | if not self.is_directed(): |
| 4307 | from sage.graphs.graph import Graph |
| 4308 | flow_graphs = map(Graph, flow_graphs) |
| 4309 | |
| 4310 | return flow_graphs |
| 4311 | |
| 4312 | def _build_flow_graph(self, flow, integer): |
| 4313 | r""" |
| 4314 | Builds a "clean" flow graph |
| 4315 | |
| 4316 | It build it, then looks for circuits and removes them |
| 4317 | |
| 4318 | INPUT: |
| 4319 | |
| 4320 | - ``flow`` -- a dictionary associating positive numerical values |
| 4321 | to edges |
| 4322 | |
| 4323 | - ``integer`` (boolean) -- whether the values from ``flow`` are the solution |
| 4324 | of an integer flow. In this case, a value of less than .5 is assumed to be 0 |
| 4325 | |
| 4326 | |
| 4327 | EXAMPLE: |
| 4328 | |
| 4329 | This method is tested in ``flow`` and ``multicommodity_flow``:: |
| 4330 | |
| 4331 | sage: g = Graph() |
| 4332 | sage: g.add_edge(0,1) |
| 4333 | sage: f = g._build_flow_graph({(0,1):1}, True) |
| 4334 | """ |
| 4335 | |
| 4336 | from sage.graphs.digraph import DiGraph |
| 4337 | g = DiGraph() |
| 4338 | |
| 4339 | # add significant edges |
| 4340 | for (u,v),l in flow.iteritems(): |
| 4341 | if l > 0 and not (integer and l<.5): |
| 4342 | g.add_edge(u,v,l) |
| 4343 | |
| 4344 | # stupid way to find Cycles. Will be fixed by #8932 |
| 4345 | # for any vertex, for any of its in-neighbors, tried to find a cycle |
| 4346 | for v in g: |
| 4347 | for u in g.neighbor_in_iterator(v): |
| 4348 | |
| 4349 | # the edge from u to v could have been removed in a previous iteration |
| 4350 | if not g.has_edge(u,v): |
| 4351 | break |
| 4352 | sp = g.shortest_path(v,u) |
| 4353 | if sp != []: |
| 4354 | |
| 4355 | #find the minimm value along the cycle. |
| 4356 | m = g.edge_label(u,v) |
| 4357 | for i in range(len(sp)-1): |
| 4358 | m = min(m,g.edge_label(sp[i],sp[i+1])) |
| 4359 | |
| 4360 | # removes it from all the edges of the cycle |
| 4361 | sp.append(v) |
| 4362 | for i in range(len(sp)-1): |
| 4363 | l = g.edge_label(sp[i],sp[i+1]) - m |
| 4364 | |
| 4365 | # an edge with flow 0 is removed |
| 4366 | if l == 0: |
| 4367 | g.delete_edge(sp[i],sp[i+1]) |
| 4368 | else: |
| 4369 | g.set_edge_label(l) |
| 4370 | |
| 4371 | # if integer is set, round values and deletes zeroes |
| 4372 | if integer: |
| 4373 | for (u,v,l) in g.edges(): |
| 4374 | if l<.5: |
| 4375 | g.delete_edge(u,v) |
| 4376 | else: |
| 4377 | g.set_edge_label(u,v, round(l)) |
| 4378 | |
| 4379 | # returning a graph with the same embedding, the corersponding name, etc ... |
| 4380 | h = self.subgraph(edges=[]) |
| 4381 | h.delete_vertices([v for v in self if (v not in g) or g.degree(v)==0]) |
| 4382 | h.add_edges(g.edges()) |
| 4383 | |
| 4384 | return h |