Ticket #8870: trac_8870.patch

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

    # HG changeset patch
    # User Nathann Cohen <nathann.cohen@gmail.com>
    # Date 1277154381 -7200
    # Node ID 0282855fb58cbfb584276aa9b10e2b6236e479f1
    # Parent  9336b249c84d56a12b08766fd90101bf63f3e0e4
    trac 8870 - Multicommodity flow
    
    diff -r 9336b249c84d -r 0282855fb58c sage/graphs/generic_graph.py
    a b  
    40664066        from sage.numerical.mip import MixedIntegerLinearProgram
    40674067        g=self
    40684068        p=MixedIntegerLinearProgram(maximization=True)
    4069         flow=p.new_variable(dim=2)
     4069        flow=p.new_variable(dim=1)
    40704070
    40714071        if use_edge_labels:
    40724072            from sage.rings.real_mpfr import RR
    40734073            capacity=lambda x: x if x in RR else 1
    40744074        else:
    40754075            capacity=lambda x: 1
    4076 
    4077         # maximizes z, which is the flow leaving from x
    40784076           
    40794077        if g.is_directed():
    40804078            # This function return the balance of flow at X
    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:
    40874111            for v in g:
    40884112                if v!=x and v!=y:
    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)
    41214114
    41224115        if integer:
    41234116            p.set_integer(flow)
    41244117
    41254118
    41264119        if value_only:
    4127             return p.solve(objective_only=True, solver=solver, log=verbose)
    4128 
    4129         obj=p.solve(solver=solver, log=verbose)
     4120            return p.solve(objective_only=True)
     4121
     4122        obj=p.solve()
     4123
    41304124        flow=p.get_values(flow)
    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   
    41344233        if g.is_directed():
    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)
    41464261                else:
    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
    41624284                else:
    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
    41664385
    41674386    def edge_disjoint_paths(self, s, t):
    41684387        r"""