Ticket #10497: trac_10497.patch

File trac_10497.patch, 16.0 KB (added by ncohen, 10 years ago)
  • sage/graphs/generic_graph.py

    # HG changeset patch
    # User Nathann Cohen <nathann.cohen@gmail.com>
    # Date 1292781822 -3600
    # Node ID bb2c0761c613beac78ffc61f6642c2fd014c9459
    # Parent  7697ba9f3ea0dd2881997620a3834a186d3c3157
    trac 10497 - Computation of Hamiltonian Cycles through constraint generation
    
    diff -r 7697ba9f3ea0 -r bb2c0761c613 sage/graphs/generic_graph.py
    a b  
    49064906        else:
    49074907            return g
    49084908
    4909     def traveling_salesman_problem(self, weighted = True, solver = None, verbose = 0):
     4909
     4910    def traveling_salesman_problem(self, weighted = True, solver = None, constraint_generation = True, verbose = 0, verbose_constraints = False):
    49104911        r"""
    49114912        Solves the traveling salesman problem (TSP)
    49124913
    4913         Given a graph (resp. a digraph) `G` with weighted edges,
    4914         the traveling salesman problem consists in finding a
    4915         Hamiltonian cycle (resp. circuit) of the graph of
    4916         minimum cost.
    4917 
    4918         This TSP is one of the most famous NP-Complete problems,
    4919         this function can thus be expected to take some time
    4920         before returning its result.
    4921 
    4922         INPUT:
    4923 
    4924         - ``weighted`` (boolean) -- whether to consider the
    4925           weights of the edges.
    4926 
    4927               - If set to ``False`` (default), all edges are
    4928                 assumed to weight `1`
    4929 
    4930               - If set to ``True``, the weights are taken into
    4931                 account, and the edges whose weight is ``None``
    4932                 are assumed to be set to `1`
     4914        Given a graph (resp. a digraph) `G` with weighted edges, the traveling
     4915        salesman problem consists in finding a hamiltonian cycle (resp. circuit)
     4916        of the graph of minimum cost.
     4917
     4918        This TSP is one of the most famous NP-Complete problems, this function
     4919        can thus be expected to take some time before returning its result.
     4920
     4921        INPUT:
     4922
     4923        - ``weighted`` (boolean) -- whether to consider the weights of the
     4924          edges.
     4925
     4926              - If set to ``False`` (default), all edges are assumed to weight
     4927                `1`
     4928
     4929              - If set to ``True``, the weights are taken into account, and the
     4930                edges whose weight is ``None`` are assumed to be set to `1`
    49334931
    49344932        - ``solver`` -- (default: ``None``) Specify a Linear Program (LP)
    49354933          solver to be used. If set to ``None``, the default one is used. For
     
    49394937          of the class
    49404938          :class:`MixedIntegerLinearProgram <sage.numerical.mip.MixedIntegerLinearProgram>`.
    49414939
     4940        - ``constraint_generation`` (boolean) -- whether to use constraint
     4941          generation when solving the Mixed Integer Linear Program (default:
     4942          ``True``).
     4943
    49424944        - ``verbose`` -- integer (default: ``0``). Sets the level of
    49434945          verbosity. Set to 0 by default, which means quiet.
    49444946
     4947        - ``verbose_constraints`` -- whether to display which constraints are
     4948          being generated.
    49454949
    49464950        OUTPUT:
    49474951
    4948         A solution to the TSP, as a ``Graph`` object whose vertex
    4949         set is `V(G)`, and whose edges are only those of the
    4950         solution.
     4952        A solution to the TSP, as a ``Graph`` object whose vertex set is `V(G)`,
     4953        and whose edges are only those of the solution.
    49514954
    49524955        ALGORITHM:
    49534956
    4954         This optimization problem is solved through the use
    4955         of Linear Programming.
     4957        This optimization problem is solved through the use of Linear
     4958        Programming.
    49564959       
    49574960        NOTE:
    49584961
    4959         - This function is correctly defined for both graph and digraphs.
    4960           In the second case, the returned cycle is a circuit of optimal
    4961           cost.
     4962        - This function is correctly defined for both graph and digraphs.  In
     4963          the second case, the returned cycle is a circuit of optimal cost.
    49624964
    49634965        EXAMPLES:
    49644966
     
    49934995            ...
    49944996            ValueError: The given graph is not Hamiltonian
    49954997
    4996         One easy way to change is is obviously to add to this
    4997         graph the edges corresponding to a Hamiltonian cycle.
     4998        One easy way to change is is obviously to add to this graph the edges
     4999        corresponding to a Hamiltonian cycle.
    49985000     
    4999         If we do this by setting the cost of these new edges
    5000         to `2`, while the others are set to `1`, we notice
    5001         that not all the edges we added are used in the
    5002         optimal solution ::
     5001        If we do this by setting the cost of these new edges to `2`, while the
     5002        others are set to `1`, we notice that not all the edges we added are
     5003        used in the optimal solution ::
    50035004
    50045005            sage: for u, v in g.edges(labels = None):
    50055006            ...      g.set_edge_label(u,v,1)
     
    50145015            sage: sum( tsp.edge_labels() ) < 2*10
    50155016            True
    50165017
    5017         If we pick `1/2` instead of `2` as a cost for these new edges,
    5018         they clearly become the optimal solution
     5018        If we pick `1/2` instead of `2` as a cost for these new edges, they
     5019        clearly become the optimal solution
    50195020
    50205021            sage: for u,v in cycle.edges(labels = None):
    50215022            ...      g.set_edge_label(u,v,1/2)
     
    50245025            sage: sum( tsp.edge_labels() ) == (1/2)*10
    50255026            True
    50265027
    5027         """
    5028        
    5029         from sage.numerical.mip import MixedIntegerLinearProgram, Sum
    5030 
    5031         p = MixedIntegerLinearProgram(maximization = False, solver = solver)
    5032 
    5033         f = p.new_variable()
    5034         r = p.new_variable()
    5035 
    5036         # If the graph has multiple edges
     5028        TESTS:
     5029
     5030        Comparing the results returned according to the value of
     5031        ``constraint_generation``. First, for graphs::
     5032
     5033            sage: from operator import itemgetter
     5034            sage: n = 20
     5035            sage: for i in range(20):
     5036            ...       g = Graph()
     5037            ...       g.allow_multiple_edges(False)
     5038            ...       for u,v in graphs.RandomGNP(n,.2).edges(labels = False):
     5039            ...            g.add_edge(u,v,round(random(),5))
     5040            ...       for u,v in graphs.CycleGraph(n).edges(labels = False):
     5041            ...            if not g.has_edge(u,v):
     5042            ...                g.add_edge(u,v,round(random(),5))
     5043            ...       v1 = g.traveling_salesman_problem(constraint_generation = False)
     5044            ...       v2 = g.traveling_salesman_problem()
     5045            ...       c1 = sum(map(itemgetter(2), v1.edges()))
     5046            ...       c2 = sum(map(itemgetter(2), v2.edges()))
     5047            ...       if c1 != c2:
     5048            ...           print "Erreur !",c1,c2
     5049            ...           break
     5050
     5051        Then for digraphs::
     5052
     5053            sage: from operator import itemgetter
     5054            sage: n = 20
     5055            sage: for i in range(20):
     5056            ...       g = DiGraph()
     5057            ...       g.allow_multiple_edges(False)
     5058            ...       for u,v in digraphs.RandomDirectedGNP(n,.2).edges(labels = False):
     5059            ...            g.add_edge(u,v,round(random(),5))
     5060            ...       for u,v in digraphs.Circuit(n).edges(labels = False):
     5061            ...            if not g.has_edge(u,v):
     5062            ...                g.add_edge(u,v,round(random(),5))
     5063            ...       v2 = g.traveling_salesman_problem()
     5064            ...       v1 = g.traveling_salesman_problem(constraint_generation = False)
     5065            ...       c1 = sum(map(itemgetter(2), v1.edges()))
     5066            ...       c2 = sum(map(itemgetter(2), v2.edges()))
     5067            ...       if c1 != c2:
     5068            ...           print "Erreur !",c1,c2
     5069            ...           print "Avec generation de contrainte :",c2
     5070            ...           print "Sans generation de contrainte :",c1
     5071            ...           break
     5072
     5073        """
     5074        ###############################
     5075        # Quick cheks of connectivity #
     5076        ###############################
     5077
     5078        # TODO : Improve it by checking vertex-connectivity instead of
     5079        # edge-connectivity.... But calling the vertex_connectivity (which
     5080        # builds a LP) is way too slow. These tests only run traversals written
     5081        # in Cython --> hence FAST
     5082
     5083        if self.is_directed():
     5084            if not self.is_strongly_connected():
     5085                raise ValueError("The given graph is not hamiltonian")
     5086
     5087        else:
     5088            # Checks whether the graph is 2-connected
     5089            if not self.strong_orientation().is_strongly_connected():
     5090                raise ValueError("The given graph is not hamiltonian")           
     5091           
     5092       
     5093
     5094        ############################
     5095        # Deal with multiple edges #
     5096        ############################
     5097
    50375098        if self.has_multiple_edges():
    50385099            g = self.copy()
    50395100            multi = self.multiple_edges()
     
    50605121        else:
    50615122            g = self
    50625123
     5124        from sage.numerical.mip import MixedIntegerLinearProgram, Sum
     5125        from sage.numerical.mip import MIPSolverException
     5126
     5127        weight = lambda l : l if (l is not None and l) else 1
     5128
     5129        ####################################################
     5130        # Constraint-generation formulation of the problem #
     5131        ####################################################
     5132
     5133        if constraint_generation:
     5134           
     5135            p = MixedIntegerLinearProgram(maximization = False,
     5136                                          solver = solver,
     5137                                          constraint_generation = True)
     5138
     5139           
     5140            # Directed Case #
     5141            #################
     5142            if g.is_directed():
     5143
     5144                from sage.graphs.all import DiGraph
     5145                b = p.new_variable(binary = True, dim = 2)
     5146
     5147                # Objective function
     5148                if weighted:
     5149                    p.set_objective(Sum([ weight(l)*b[u][v]
     5150                                          for u,v,l in g.edges()]))
     5151
     5152                # All the vertices have in-degree 1 and out-degree 1
     5153                for v in g:
     5154                    p.add_constraint(Sum([b[u][v] for u in g.neighbors_in(v)]),
     5155                                     min = 1,
     5156                                     max = 1)
     5157                    p.add_constraint(Sum([b[v][u] for u in g.neighbors_out(v)]),
     5158                                     min = 1,
     5159                                     max = 1)
     5160
     5161                # Initial Solve
     5162                try:
     5163                    p.solve(log = verbose)
     5164                except MIPSolverException:
     5165                    raise ValueError("The given graph is not hamiltonian")
     5166                   
     5167
     5168                while True:
     5169                    # We build the DiGraph representing the current solution
     5170                    h = DiGraph()
     5171                    for u,v,l in g.edges():
     5172                        if p.get_values(b[u][v]) > .5:
     5173                            h.add_edge(u,v,l)
     5174
     5175                    # If there is only one circuit, we are done !
     5176                    cc = h.connected_components()
     5177                    if len(cc) == 1:
     5178                        break
     5179
     5180                    # Adding the corresponding constraint                   
     5181                    if verbose_constraints:
     5182                        print "Adding a constraint on set",cc[0]
     5183
     5184
     5185                    p.add_constraint(Sum(b[u][v] for u,v in
     5186                                         g.edge_boundary(cc[0], labels = False)),
     5187                                     min = 1)
     5188
     5189                    try:
     5190                        p.solve(log = verbose)
     5191                    except MIPSolverException:
     5192                        raise ValueError("The given graph is not hamiltonian")
     5193
     5194            # Undirected Case #
     5195            ###################
     5196            else:
     5197
     5198                from sage.graphs.all import Graph
     5199                b = p.new_variable(binary = True)
     5200                B = lambda u,v : b[(u,v)] if u<v else b[(v,u)]
     5201
     5202                # Objective function
     5203                if weighted:
     5204                    p.set_objective(Sum([ weight(l)*B(u,v)
     5205                                          for u,v,l in g.edges()]) )
     5206
     5207                # All the vertices have degree 2
     5208                for v in g:
     5209                    p.add_constraint(Sum([ B(u,v) for u in g.neighbors(v)]),
     5210                                     min = 2,
     5211                                     max = 2)
     5212
     5213                # Initial Solve
     5214                try:
     5215                    p.solve(log = verbose)
     5216                except MIPSolverException:
     5217                    raise ValueError("The given graph is not hamiltonian")
     5218
     5219                while True:
     5220                    # We build the DiGraph representing the current solution
     5221                    h = Graph()
     5222                    for u,v,l in g.edges():
     5223                        if p.get_values(B(u,v)) > .5:
     5224                            h.add_edge(u,v,l)
     5225
     5226                    # If there is only one circuit, we are done !
     5227                    cc = h.connected_components()
     5228                    if len(cc) == 1:
     5229                        break
     5230
     5231                    # Adding the corresponding constraint
     5232                    if verbose_constraints:
     5233                        print "Adding a constraint on set",cc[0]
     5234
     5235                    p.add_constraint(Sum(B(u,v) for u,v in
     5236                                         g.edge_boundary(cc[0], labels = False)),
     5237                                     min = 2)
     5238
     5239                    try:
     5240                        p.solve(log = verbose)
     5241                    except MIPSolverException:
     5242                        raise ValueError("The given graph is not hamiltonian")
     5243
     5244
     5245
     5246            # We can now return the TSP !
     5247            answer = self.subgraph(edges = h.edges())
     5248            answer.set_pos(self.get_pos())
     5249            answer.name("TSP from "+g.name())
     5250            return answer
     5251
     5252        #################################################
     5253        # ILP formulation without constraint generation #
     5254        #################################################
     5255
     5256
     5257
     5258        p = MixedIntegerLinearProgram(maximization = False, solver = solver)
     5259
     5260        f = p.new_variable()
     5261        r = p.new_variable()
     5262
    50635263        eps = 1/(2*Integer(g.order()))
    50645264        x = g.vertex_iterator().next()
    50655265
     
    50695269            # returns the variable corresponding to arc u,v
    50705270            E = lambda u,v : f[(u,v)]
    50715271
    5072             # All the vertices have in-degree 1
    5073             [p.add_constraint(Sum([ f[(u,v)] for u in g.neighbors_in(v)]), min = 1, max = 1) for v in g]
    5074 
    5075             # All the vertices have out-degree 1
    5076             [p.add_constraint(Sum([ f[(v,u)] for u in g.neighbors_out(v)]), min = 1, max = 1) for v in g]
     5272            # All the vertices have in-degree 1 and out-degree 1
     5273            for v in g:
     5274                p.add_constraint(Sum([ f[(u,v)] for u in g.neighbors_in(v)]),
     5275                                 min = 1,
     5276                                 max = 1)
     5277
     5278                p.add_constraint(Sum([ f[(v,u)] for u in g.neighbors_out(v)]),
     5279                                 min = 1,
     5280                                 max = 1)
    50775281
    50785282
    50795283            # r is greater than f
     
    51015305            E = lambda u,v : f[R(u,v)]
    51025306
    51035307            # All the vertices have degree 2
    5104             [p.add_constraint(Sum([ f[R(u,v)] for u in g.neighbors(v)]), min = 2, max = 2) for v in g]
     5308            for v in g:
     5309                p.add_constraint(Sum([ f[R(u,v)] for u in g.neighbors(v)]),
     5310                                 min = 2,
     5311                                 max = 2)
    51055312
    51065313            # r is greater than f
    51075314            for u,v in g.edges(labels = None):
     
    51195326                p.add_constraint(Sum([ r[(u,v)] for u in g.neighbors(v)]),max = 1-eps)
    51205327
    51215328
    5122         weight = lambda u,v : g.edge_label(u,v) if (g.edge_label(u,v) is not None and g.edge_label(u,v) != {}) else 1
     5329
    51235330
    51245331        if weighted:
    5125             p.set_objective(Sum([ weight(u,v)*E(u,v) for u,v in g.edges(labels=None)]) )
     5332            p.set_objective(Sum([ weight(l)*E(u,v) for u,v,l in g.edges()]) )
    51265333        else:
    51275334            p.set_objective(None)
    51285335
    51295336        p.set_binary(f)
    51305337
    5131         from sage.numerical.mip import MIPSolverException
     5338
    51325339
    51335340        try:
    51345341            obj = p.solve(log = verbose)