Ticket #7288: trac_7288.patch

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

    # HG changeset patch
    # User Nathann Cohen <nathann.cohen@gmail.com>
    # Date 1268382978 -3600
    # Node ID f28f7f29eaf0b2432c15571286384acb35d91208
    # Parent  deb9406f163f6048dd0b93813d158541393cf71b
    [mq]: hu.patch
    
    diff -r deb9406f163f -r f28f7f29eaf0 sage/graphs/generic_graph.py
    a b  
    30223022            v = pv
    30233023        return B, C
    30243024               
    3025     def edge_cut(self,s,t,value_only=True,use_edge_labels=False):
     3025    def edge_cut(self, s, t, value_only=True, use_edge_labels=False, vertices=False):
    30263026        r"""
    30273027        Returns a minimum edge cut between vertices `s` and `t`
     3028        represented by a list of edges.
     3029
     3030        A minimum edge cut between two vertices `s` and `t` of self
     3031        is a set `U` of edges of minimum weight such that the graph
     3032        obtained by removing `U` from self is disconnected.
    30283033        ( cf. http://en.wikipedia.org/wiki/Cut_%28graph_theory%29 )
    3029         represented by a list of edges.
    3030    
    3031        
    3032         INPUT:
    3033 
    3034         - ``s`` -- Source vertex
    3035 
    3036         - ``t`` -- Sink vertex
    3037    
    3038         - ``value_only`` (boolean)
    3039 
    3040             - When set to ``True``, only the value of a maximal
    3041               flow is returned.
    3042 
    3043             - When set to ``False``, a list of edges in the minimum cut is
    3044               also returned.
    3045 
    3046         - ``use_edge_labels`` (boolean)
    3047             - When set to ``True``, computes a weighted minimum cut
    3048               where each edge has a weight defined by its label. ( if
    3049               an edge has no label, `1` is assumed )
    3050 
    3051             - when set to ``False`` (default), each edge has weight `1`.
    3052    
    3053         EXAMPLE:
     3034       
     3035        INPUT:
     3036
     3037        - ``s`` -- source vertex
     3038        - ``t`` -- sink vertex
     3039        - ``value_only`` -- boolean (default: True). When set to
     3040          True, only the weight of a minimum cut is returned.
     3041          Otherwise, a list of edges of a minimum cut is also returned.
     3042        - ``use_edge_labels`` -- boolean (default: False). When set to
     3043          True, computes a weighted minimum cut where each edge has
     3044          a weight defined by its label (if an edge has no label, `1`
     3045          is assumed). Otherwise, each edge has weight `1`.
     3046        - ``vertices`` -- boolean (default: False). When set to True,
     3047          also returns the two sets of vertices that are disconnected by
     3048          the cut. Implies ``value_only=False``.
     3049
     3050        OUTPUT:
     3051
     3052        real number or tuple, depending on the arguments given
     3053        (examples are given below)
     3054
     3055        EXAMPLES:
    30543056       
    3055         A basic application in a ``PappusGraph``::
    3056        
    3057            sage: g=graphs.PappusGraph()
    3058            sage: g.edge_cut(1,2,value_only=True) # optional - requires Glpk or COIN-OR/CBC
     3057        A basic application in the Pappus graph::
     3058       
     3059           sage: g = graphs.PappusGraph()
     3060           sage: g.edge_cut(1, 2, value_only=True) # optional - requires GLPK or COIN-OR/CBC
    30593061           3.0
    30603062
    3061         If the graph is a path with weighted edges, the edge cut between the two ends
    3062         is the edge of minimum weight ::
     3063        If the graph is a path with randomly weighted edges::
    30633064
    30643065           sage: g = graphs.PathGraph(15)
    30653066           sage: for (u,v) in g.edge_iterator(labels=None):
    30663067           ...      g.set_edge_label(u,v,random())
     3068
     3069        The edge cut between the two ends is the edge of minimum weight::
     3070
    30673071           sage: minimum = min([l for u,v,l in g.edge_iterator()])
    3068            sage: minimum == g.edge_cut(0,14,use_edge_labels=True)                                # optional - requires Glpk or COIN-OR/CBC
    3069            True
    3070            sage: [value,[[u,v]]] = g.edge_cut(0,14,use_edge_labels=True, value_only=False)       # optional - requires Glpk or COIN-OR/CBC
    3071            sage: g.edge_label(u,v) == minimum                                                    # optional - requires Glpk or COIN-OR/CBC
     3072           sage: minimum == g.edge_cut(0, 14, use_edge_labels=True) # optional - requires GLPK or COIN-OR/CBC
     3073           True
     3074           sage: [value,[[u,v]]] = g.edge_cut(0, 14, use_edge_labels=True, value_only=False) # optional - requires GLPK or COIN-OR/CBC
     3075           sage: g.edge_label(u, v) == minimum # optional - requires GLPK or COIN-OR/CBC
     3076           True
     3077
     3078        The two sides of the edge cut are obviously shorter paths::
     3079
     3080           sage: value,edges,[set1,set2] = g.edge_cut(0, 14, use_edge_labels=True, vertices=True)  # optional - requires Glpk or COIN-OR/CBC
     3081           sage: g.subgraph(set1).is_isomorphic(graphs.PathGraph(len(set1)))                     # optional - requires Glpk or COIN-OR/CBC
     3082           True
     3083           sage: g.subgraph(set2).is_isomorphic(graphs.PathGraph(len(set2)))                     # optional - requires Glpk or COIN-OR/CBC
     3084           True
     3085           sage: len(set1) + len(set2) == g.order()                                                # optional - requires Glpk or COIN-OR/CBC
    30723086           True
    30733087        """
    30743088        from sage.numerical.mip import MixedIntegerLinearProgram
    3075         g=self
    3076         p=MixedIntegerLinearProgram(maximization=False)
    3077         b=p.new_variable(dim=2)
    3078         v=p.new_variable()
    3079 
     3089        g = self
     3090        p = MixedIntegerLinearProgram(maximization=False)
     3091        b = p.new_variable(dim=2)
     3092        v = p.new_variable()
     3093
     3094        if vertices:
     3095            value_only = False
    30803096        if use_edge_labels:
    3081             weight=lambda x: 1 if x==None else x
    3082         else:
    3083             weight=lambda x: 1
     3097            weight = lambda x: 1 if x == None else x
     3098        else:
     3099            weight = lambda x: 1
    30843100
    30853101        # Some vertices belong to part 1, others to part 0
    3086         p.add_constraint(v[s],min=0,max=0)
    3087         p.add_constraint(v[t],min=1,max=1)       
     3102        p.add_constraint(v[s], min=0, max=0)
     3103        p.add_constraint(v[t], min=1, max=1)       
    30883104
    30893105        if g.is_directed():
     3106
    30903107            # we minimize the number of edges
    3091             p.set_objective(sum([weight(w)*b[x][y] for (x,y,w) in g.edges()]))
    3092            
    3093             # Adjacent vertices can belong to different parts only if the
    3094             # edge that connects them is part of the cut
    3095             [p.add_constraint(v[x]+b[x][y]-v[y],min=0,max=0) for (x,y) in g.edges(labels=None)]
    3096            
    3097            
    3098         else:
    3099             # we minimize the number of edges
    3100             p.set_objective(sum([weight(w)*b[min(x,y)][max(x,y)] for (x,y,w) in g.edges()]))
     3108            p.set_objective(sum([weight(w) * b[x][y] for (x,y,w) in g.edges()]))
    31013109
    31023110            # Adjacent vertices can belong to different parts only if the
    31033111            # edge that connects them is part of the cut
    31043112            for (x,y) in g.edges(labels=None):
    3105                 p.add_constraint(v[x]+b[min(x,y)][max(x,y)]-v[y],min=0)
    3106                 p.add_constraint(v[y]+b[min(x,y)][max(x,y)]-v[x],min=0)
     3113                p.add_constraint(v[x] + b[x][y] - v[y], min=0, max=0)
     3114
     3115        else:
     3116            # we minimize the number of edges
     3117            p.set_objective(sum([weight(w) * b[min(x,y)][max(x,y)] for (x,y,w) in g.edges()]))
     3118            # Adjacent vertices can belong to different parts only if the
     3119            # edge that connects them is part of the cut
     3120            for (x,y) in g.edges(labels=None):
     3121                p.add_constraint(v[x] + b[min(x,y)][max(x,y)] - v[y], min=0)
     3122                p.add_constraint(v[y] + b[min(x,y)][max(x,y)] - v[x], min=0)
    31073123
    31083124        p.set_binary(v)
    31093125        p.set_binary(b)
     
    31113127        if value_only:
    31123128            return p.solve(objective_only=True)
    31133129        else:
    3114             obj=p.solve()
    3115             b=p.get_values(b)
     3130            obj = p.solve()
     3131            b = p.get_values(b)
     3132            answer = [obj]
    31163133            if g.is_directed():
    3117                 return [obj,[(x,y) for (x,y) in g.edges(labels=None) if b[x][y]==1]]
    3118             else:
    3119                 return [obj,[(x,y) for (x,y) in g.edges(labels=None) if b[min(x,y)][max(x,y)]==1]]
    3120        
    3121     def vertex_cut(self,s,t,value_only=True):
    3122         r"""
    3123         Returns a minimum vertex cut between non-adjacent vertices `s` and `t`
     3134                answer.append([(x,y) for (x,y) in g.edges(labels=None) if b[x][y] == 1])
     3135            else:
     3136                answer.append([(x,y) for (x,y) in g.edges(labels=None) if b[min(x,y)][max(x,y)] == 1])
     3137
     3138            if vertices:
     3139                v = p.get_values(v)
     3140                l0 = []
     3141                l1 = []
     3142                for x in g.vertex_iterator():
     3143                    if v.has_key(x) and v[x] == 1:
     3144                        l1.append(x)
     3145                    else:
     3146                        l0.append(x)
     3147                answer.append([l0, l1])
     3148            return tuple(answer)
     3149       
     3150    def vertex_cut(self, s, t, value_only=True, vertices=False):
     3151        r"""
     3152        Returns a minimum vertex cut between non adjacent vertices `s` and `t`
     3153        represented by a list of vertices.
     3154
     3155        A vertex cut between two non adjacent vertices is a set `U`
     3156        of vertices of self such that the graph obtained by removing
     3157        `U` from self is disconnected.
    31243158        ( cf. http://en.wikipedia.org/wiki/Cut_%28graph_theory%29 )
    3125         represented by a list of vertices.
    3126    
    3127        
    3128         INPUT:
    3129    
    3130         - ``value_only`` : When set to ``True``, only the cardinal of the
    3131           minimum cut is returned
    3132    
    3133         EXAMPLE:
    3134 
    3135         A basic application in a ``PappusGraph``::
    3136        
    3137            sage: g=graphs.PappusGraph()
    3138            sage: g.vertex_cut(1,16,value_only=True) # optional - requires Glpk or COIN-OR/CBC
     3159   
     3160       
     3161        INPUT:
     3162   
     3163        - ``value_only`` -- boolean (default: True). When set to
     3164          True, only the size of the minimum cut is returned
     3165        - ``vertices`` -- boolean (default: False). When set to
     3166          True, also returns the two sets of vertices that
     3167          are disconnected by the cut. Implies ``value_only``
     3168          set to False.
     3169
     3170        OUTPUT:
     3171
     3172        real number or tuple, depending on the arguments given
     3173        (examples are given below)
     3174   
     3175        EXAMPLE:
     3176
     3177        A basic application in a Pappus graph::
     3178       
     3179           sage: g = graphs.PappusGraph()
     3180           sage: g.vertex_cut(1, 16, value_only=True) # optional - requires Glpk or COIN-OR/CBC
    31393181           3.0
    31403182
    3141         Or in a bipartite complete graph `K_{2,8}` between the first 2
    3142         vertices, in which case the minimum cut is the other set of 8
    3143         vertices ::
    3144 
    3145            sage: g = graphs.CompleteBipartiteGraph(2,8)
    3146            sage: [value, vertices] = g.vertex_cut(0,1, value_only=False) # optional - requires Glpk or COIN-OR/CBC
    3147            sage: print value                                             # optional - requires Glpk or COIN-OR/CBC
     3183        In the bipartite complete graph `K_{2,8}`, a cut between the two
     3184        vertices in the size `2` part consists of the other `8` vertices::
     3185
     3186           sage: g = graphs.CompleteBipartiteGraph(2, 8)
     3187           sage: [value, vertices] = g.vertex_cut(0, 1, value_only=False) # optional - requires Glpk or COIN-OR/CBC
     3188           sage: print value # optional - requires Glpk or COIN-OR/CBC
    31483189           8.0
    3149            sage: vertices == range(2,10)                                 # optional - requires Glpk or COIN-OR/CBC
     3190           sage: vertices == range(2,10) # optional - requires Glpk or COIN-OR/CBC
     3191           True
     3192
     3193        Clearly, in this case the two sides of the cut are singletons ::
     3194
     3195           sage: [value, vertices, [set1, set2]] = g.vertex_cut(0,1, vertices=True) # optional - requires Glpk or COIN-OR/CBC
     3196           sage: len(set1) == 1 # optional - requires Glpk or COIN-OR/CBC
     3197           True
     3198           sage: len(set2) == 1 # optional - requires Glpk or COIN-OR/CBC
    31503199           True
    31513200        """
    31523201        from sage.numerical.mip import MixedIntegerLinearProgram
    3153         g=self
     3202        g = self
    31543203        if g.has_edge(s,t):
    31553204            raise ValueError, "There can be no vertex cut between adjacent vertices !"
    3156 
    3157         p=MixedIntegerLinearProgram(maximization=False)
    3158         b=p.new_variable()
    3159         v=p.new_variable()
     3205        if vertices:
     3206            value_only = False
     3207
     3208        p = MixedIntegerLinearProgram(maximization=False)
     3209        b = p.new_variable()
     3210        v = p.new_variable()
    31603211
    31613212        # Some vertices belong to part 1, some others to part 0
    3162         p.add_constraint(v[s],min=0,max=0)
    3163         p.add_constraint(v[t],min=1,max=1)
     3213        p.add_constraint(v[s], min=0, max=0)
     3214        p.add_constraint(v[t], min=1, max=1)
    31643215
    31653216        # b indicates whether the vertices belong to the cut
    3166         p.add_constraint(b[s],min=0,max=0)
    3167         p.add_constraint(b[t],min=0,max=0)
     3217        p.add_constraint(b[s], min=0, max=0)
     3218        p.add_constraint(b[t], min=0, max=0)
    31683219
    31693220        if g.is_directed():
     3221
    31703222            p.set_objective(sum([b[x] for x in g.vertices()]))
     3223
    31713224            # adjacent vertices belong to the same part except if one of them
    31723225            # belongs to the cut
    3173             [p.add_constraint(v[x]+b[y]-v[y],min=0) for (x,y) in g.edges(labels=None)]
     3226            for (x,y) in g.edges(labels=None):
     3227                p.add_constraint(v[x] + b[y] - v[y], min=0)
    31743228
    31753229        else:
    31763230            p.set_objective(sum([b[x] for x in g.vertices()]))
    31773231            # adjacent vertices belong to the same part except if one of them
    31783232            # belongs to the cut
    31793233            for (x,y) in g.edges(labels=None):
    3180                 p.add_constraint(v[x]+b[y]-v[y],min=0)
    3181                 p.add_constraint(v[y]+b[x]-v[x],min=0)
     3234                p.add_constraint(v[x] + b[y] - v[y],min=0)
     3235                p.add_constraint(v[y] + b[x] - v[x],min=0)
    31823236
    31833237        p.set_binary(b)
    31843238        p.set_binary(v)
     
    31863240        if value_only:
    31873241            return p.solve(objective_only=True)
    31883242        else:
    3189             obj=p.solve()
    3190             b=p.get_values(b)
    3191             return [obj,[v for v in g if b[v]==1]]
     3243            obj = p.solve()
     3244            b = p.get_values(b)
     3245            answer = [obj,[x for x in g if b[x] == 1]]
     3246            if vertices:
     3247                v = p.get_values(v)
     3248                l0 = []
     3249                l1 = []
     3250                for x in g.vertex_iterator():
     3251                    # if the vertex is not in the cut
     3252                    if not (b.has_key(x) and b[x] == 1):
     3253                        if (v.has_key(x) and v[x] == 1):
     3254                            l1.append(x)
     3255                        else:
     3256                            l0.append(x)
     3257                answer.append([l0, l1])
     3258            return tuple(answer)
     3259
    31923260
    31933261    def vertex_cover(self,algorithm="Cliquer",value_only=False,log=0):
    31943262        r"""
    3195         Returns a minimum vertex cover of the graph
    3196         ( cf. http://en.wikipedia.org/wiki/Vertex_cover )
    3197         represented by the list of its vertices.
    3198 
     3263        Returns a minimum vertex cover of self represented
     3264        by a list of vertices.
     3265       
    31993266        A minimum vertex cover of a graph is a set `S` of
    32003267        its vertices such that each edge is incident to at least
    32013268        one element of `S`, and such that `S` is of minimum
    32023269        cardinality.
    3203 
    3204         A vertex cover is also defined as the complement of an
    3205         independent set.
    3206 
    3207         As an optimization problem, it can be expressed as :
     3270        ( cf. http://en.wikipedia.org/wiki/Vertex_cover )
     3271
     3272        Equivalently, a vertex cover is defined as the
     3273        complement of an independent set.
     3274
     3275        As an optimization problem, it can be expressed as follows:
    32083276
    32093277        .. MATH::
    32103278            \mbox{Minimize : }&\sum_{v\in G} b_v\\
     
    32133281   
    32143282        INPUT:
    32153283       
    3216         - ``algorithm`` -- Select an algorithm
    3217        
    3218             - ``"Cliquer"`` (default) will return a minimum vertex cover
     3284        - ``algorithm`` -- string (default: ``"Cliquer") indicating
     3285          which algorithm is performed. It can be one of those two values.
     3286          - ``"Cliquer"`` will compute a minimum vertex cover
    32193287              using the algorithm Cliquer.
    3220 
    3221             - ``"MILP"`` will return a minimum vertex cover through a Mixed
    3222               Integer Linear Program ( requires packages GLPK or CBC )
    3223 
    3224         - ``value_only`` --
    3225             - If ``value_only = True``, only the cardinality of a
    3226               minimum vertex cover is returned.
    3227             - If ``value_only = False`` ( default ), a minimum vertex cover
    3228               is returned as the list of its vertices
    3229    
    3230         - ``log`` ( integer ) --
    3231           As vertex cover is an `NP`-complete problem, its solving may take some time
    3232           depending on the graph. Use ``log`` to define the level of verbosity you want
    3233           from the linear program solver.
    3234 
    3235           By default ``log=0``, meaning that there will be no message printed by the solver.
    3236 
    3237           Only useful if ``algorithm="MILP"``.
     3288            - ``"MILP"`` will compute a minimum vertex cover through a mixed
     3289              integer linear program (requires packages GLPK or CBC).
     3290        - ``value_only`` -- boolean (default: False). If set to True,
     3291          only the size of a minimum vertex cover is returned. Otherwise,
     3292          a minimum vertex cover is returned as a list of vertices.
     3293        - ``log`` -- non negative integer (default: 0) precising the level
     3294          of verbosity you want from the linear program solver. Since the
     3295          problem of computing a vertex cover is `NP`-complete, its solving
     3296          may take some time depending on the graph. A value of 0 means
     3297          that there will be no message printed by the solver. Only useful
     3298          if ``algorithm="MILP"``.
    32383299   
    32393300        EXAMPLES:
    32403301
    3241         On a ``PappusGraph`` ::
    3242        
    3243            sage: g=graphs.PappusGraph()
     3302        On the Pappus graph ::
     3303       
     3304           sage: g = graphs.PappusGraph()
    32443305           sage: g.vertex_cover(value_only=True)
    32453306           9
    32463307
    3247         Obviously, the two methods return the same results ::
    3248 
    3249            sage: g=graphs.RandomGNP(10,.5)
    3250            sage: vc1 = g.vertex_cover(algorithm="MILP")     # optional requires GLPK or CBC
    3251            sage: vc2 = g.vertex_cover(algorithm="Cliquer")  # optional requires GLPK or CBC
     3308        The two algorithms should return the same result::
     3309
     3310           sage: g = graphs.RandomGNP(10,.5)
     3311           sage: vc1 = g.vertex_cover(algorithm="MILP") # optional requires GLPK or CBC
     3312           sage: vc2 = g.vertex_cover(algorithm="Cliquer") # optional requires GLPK or CBC
    32523313           sage: len(vc1) == len(vc2) # optional requires GLPK or CBC
    32533314           True
    32543315        """
    3255 
    3256 
    3257         if algorithm=="Cliquer":
     3316        if algorithm == "Cliquer":
    32583317            independent = self.independent_set()
    3259            
    32603318            if value_only:
    3261                 return self.order()-len(independent)
     3319                return self.order() - len(independent)
    32623320            else:
    32633321                return set(self.vertices()).difference(set(independent))
    32643322
    3265         elif algorithm=="MILP":
     3323        elif algorithm == "MILP":
    32663324
    32673325            from sage.numerical.mip import MixedIntegerLinearProgram
    3268             g=self
    3269             p=MixedIntegerLinearProgram(maximization=False)
    3270             b=p.new_variable()
     3326            g = self
     3327            p = MixedIntegerLinearProgram(maximization=False)
     3328            b = p.new_variable()
    32713329   
    32723330            # minimizes the number of vertices in the set
    32733331            p.set_objective(sum([b[v] for v in g.vertices()]))
    32743332   
    32753333            # an edge contains at least one vertex of the minimum vertex cover
    3276             [p.add_constraint(b[u]+b[v],min=1) for (u,v) in g.edges(labels=None)]
     3334            for (u,v) in g.edges(labels=None):
     3335                p.add_constraint(b[u] + b[v], min=1)
    32773336           
    32783337            p.set_binary(b)
    32793338   
    32803339            if value_only:
    3281                 return p.solve(objective_only=True,log=log)
     3340                return p.solve(objective_only=True, log=log)
    32823341            else:
    32833342                p.solve()
    3284                 b=p.get_values(b)
    3285                 return set([v for v in g.vertices() if b[v]==1])
     3343                b = p.get_values(b)
     3344                return set([v for v in g.vertices() if b[v] == 1])
    32863345        else:
    32873346            raise ValueError("Only two algorithms are available : Cliquer and MILP.")
    32883347
     
    33723431        N=range(n)
    33733432
    33743433        # First and second constraints
    3375         [p.add_constraint(sum([x[v][i] for i in N]),min=1,max=1) for v in self]
    3376         [p.add_constraint(sum([x[v][i] for v in self]),min=1,max=1) for i in N]
     3434        for v in self:
     3435            p.add_constraint(sum([x[v][i] for i in N]),min=1,max=1)
     3436
     3437        for i in N:
     3438            p.add_constraint(sum([x[v][i] for v in self]),min=1,max=1)
    33773439       
    33783440        # Definition of d_v
    3379         [p.add_constraint(sum([i*x[v][i] for i in N])-d[v],max=0,min=0) for v in self]
     3441        for v in self:
     3442            p.add_constraint(sum([i*x[v][i] for i in N])-d[v],max=0,min=0)
    33803443
    33813444        # The removed vertices cover all the back arcs ( third condition )
    3382         [p.add_constraint(d[u]-d[v]+n*(b[(u,v)]),min=0) for (u,v) in self.edges(labels=None)]
     3445        for (u,v) in self.edges(labels=None):
     3446            p.add_constraint(d[u]-d[v]+n*(b[(u,v)]),min=0)
    33833447
    33843448        p.set_binary(b)
    33853449        p.set_binary(x)
     
    34803544        N=range(n)
    34813545
    34823546        # First and second constraints
    3483         [p.add_constraint(sum([x[v][i] for i in N]),min=1,max=1) for v in self]
    3484         [p.add_constraint(sum([x[v][i] for v in self]),min=1,max=1) for i in N]
     3547        for v in self:
     3548            p.add_constraint(sum([x[v][i] for i in N]),min=1,max=1)
     3549
     3550        for i in N:
     3551            p.add_constraint(sum([x[v][i] for v in self]),min=1,max=1)
    34853552       
    34863553        # Definition of d_v
    3487         [p.add_constraint(sum([i*x[v][i] for i in N])-d[v],max=0,min=0) for v in self]
     3554        for v in self:
     3555            p.add_constraint(sum([i*x[v][i] for i in N])-d[v],max=0,min=0)
    34883556
    34893557        # The removed vertices cover all the back arcs ( third condition )
    3490         [p.add_constraint(d[u]-d[v]+n*(b[u]+b[v]),min=0) for (u,v) in self.edges(labels=None)]
     3558        for (u,v) in self.edges(labels=None):
     3559            p.add_constraint(d[u]-d[v]+n*(b[u]+b[v]),min=0)
    34913560
    34923561        p.set_binary(b)
    34933562        p.set_binary(x)
     
    37083777        a sink `t` to all the vertices of the second set, then computing
    37093778        a maximum `s-t` flow ::
    37103779
    3711             sage: g = graphs.CompleteBipartiteGraph(4,4)
    3712             sage: g.add_edges([('s',i) for i in range(4)])
     3780            sage: g = DiGraph()
     3781            sage: g.add_edges([('s',i) for i in range(4)])
     3782            sage: g.add_edges([(i,4+j) for i in range(4) for j in range(4)])
    37133783            sage: g.add_edges([(4+i,'t') for i in range(4)])
    37143784            sage: [cardinal, flow_graph] = g.flow('s','t',integer=True,value_only=False) # optional - requries GLPK or CBC
    37153785            sage: flow_graph.delete_vertices(['s','t'])                                  # optional - requries GLPK or CBC
    3716             sage: flow_graph.edges(labels=None)                                          # optional - requries GLPK or CBC
    3717             [(0, 5), (1, 4), (2, 7), (3, 6)]
    3718 
     3786            sage: len(flow_graph.edges(labels=None))                                     # optional - requries GLPK or CBC   
     3787            4
    37193788       
    37203789        """
    37213790        from sage.numerical.mip import MixedIntegerLinearProgram
     
    37383807            p.set_objective(flow_sum(x))
    37393808
    37403809            # Elsewhere, the flow is equal to 0-
    3741             [p.add_constraint(flow_sum(v),min=0,max=0) for v in g if v!=x and v!=y]
     3810            for v in g:
     3811                if v!=x and v!=y:
     3812                    p.add_constraint(flow_sum(v),min=0,max=0)
    37423813
    37433814            # Capacity constraints
    3744             [p.add_constraint(flow[u][v],max=capacity(w)) for (u,v,w) in g.edges()]
     3815            for (u,v,w) in g.edges():
     3816                p.add_constraint(flow[u][v],max=capacity(w))
    37453817
    37463818            if vertex_bound:
    3747                 [p.add_constraint(sum([flow[uu][vv] for (uu,vv) in g.outgoing_edges([v],labels=None)]),max=1) for v in g.vertices() if v!=x and v!=y]
     3819                for v in g.vertices():
     3820                    if v!=x and v!=y:
     3821                        p.add_constraint(sum([flow[uu][vv] for (uu,vv) in g.outgoing_edges([v],labels=None)]),max=1)
    37483822
    37493823        else:
    37503824            # This function return the balance of flow at X
     
    37543828            p.set_objective(flow_sum(x))
    37553829
    37563830            # Elsewhere, the flow is equal to 0
    3757             [p.add_constraint(flow_sum(v),min=0,max=0) for v in g if v!=x and v!=y]
     3831            for v in g:
     3832                if v!=x and v!=y:
     3833                    p.add_constraint(flow_sum(v),min=0,max=0)
    37583834
    37593835            # Capacity constraints
    3760             [p.add_constraint(flow[u][v]+flow[v][u],max=capacity(w)) for (u,v,w) in g.edges()]
     3836            for (u,v,w) in g.edges():
     3837                p.add_constraint(flow[u][v]+flow[v][u],max=capacity(w))
    37613838           
    37623839            if vertex_bound:
    3763                 [p.add_constraint([flow[X][v] for X in g[v]],max=1) for v in g if v!=x and v!=y]
     3840                for v in g:
     3841                    if v!=x and v!=y:
     3842                        p.add_constraint([flow[X][v] for X in g[v]],max=1)
    37643843           
    37653844
    37663845        if integer:
     
    39394018        p=MixedIntegerLinearProgram(maximization=True)
    39404019
    39414020        b=p.new_variable(dim=2)
    3942         [p.set_objective(sum([weight(w)*b[min(u,v)][max(u,v)] for (u,v,w) in g.edges()]))]
     4021        p.set_objective(sum([weight(w)*b[min(u,v)][max(u,v)] for (u,v,w) in g.edges()]))
    39434022
    39444023
    39454024        # for any vertex v, there is at most one edge incident to v in the maximum matching
    3946         [p.add_constraint(sum([b[min(u,v)][max(u,v)] for u in g.neighbors(v)]),max=1) for v in g.vertices()]
     4025        for v in g.vertices():
     4026            p.add_constraint(sum([b[min(u,v)][max(u,v)] for u in g.neighbors(v)]),max=1)
    39474027
    39484028        p.set_binary(b)
    39494029   
     
    40194099       
    40204100        # For any vertex v, one of its neighbors or v itself is in
    40214101        # the minimum dominating set
    4022         [p.add_constraint(b[v]+sum([b[u] for u in g.neighbors(v)]),min=1) for v in g.vertices()]
     4102        for v in g.vertices():
     4103            p.add_constraint(b[v]+sum([b[u] for u in g.neighbors(v)]),min=1)
    40234104
    40244105
    40254106        if independent:
    40264107            # no two adjacent vertices are in the set
    4027             [p.add_constraint(b[u]+b[v],max=1) for (u,v) in g.edges(labels=None)]
     4108            for (u,v) in g.edges(labels=None):
     4109                p.add_constraint(b[u]+b[v],max=1)
    40284110
    40294111        # Minimizes the number of vertices used
    40304112        p.set_objective(sum([b[v] for v in g.vertices()]))
  • sage/graphs/graph.py

    diff -r deb9406f163f -r f28f7f29eaf0 sage/graphs/graph.py
    a b  
    27342734        else:
    27352735            raise NotImplementedError, "Minimum Spanning Tree algorithm '%s' is not implemented."%algorithm
    27362736
     2737
     2738    def _gomory_hu_tree(self, vertices=None):
     2739        r"""
     2740        Returns the Gomory-Hu tree associated to self (private function)
     2741
     2742        This function is the privatre counterpart of ``gomory_hu_tree``,
     2743        with the difference that it accepts an optional argument
     2744        needed for recursive computations, which the user is not
     2745        interested in defining by himself.
     2746
     2747        See the documentation of ``gomory_hu_tree`` for more information.
     2748
     2749        INPUT:
     2750
     2751        - ``vertices`` -- a set of ''real'' vertices, as opposed to the
     2752          fakes one introduced during the computations. This variable is
     2753          useful for the algorithm, and I see no reason for the user to
     2754          change it.
     2755
     2756        EXAMPLE:
     2757
     2758        This function is actually tested in ``gomory_hu_tree``, this
     2759        example is only present to have a doctest coverage of 100%
     2760
     2761            sage: g = graphs.PetersenGraph()
     2762            sage: t = g._gomory_hu_tree()      # optional - requires GLPK or CBC
     2763
     2764        """
     2765        from sage.sets.set import Set
     2766       
     2767        # The default capacity of an arc is 1
     2768        capacity = lambda label: label if label is not None else 1
     2769
     2770        # Keeping the graph's embedding
     2771        pos = False
     2772
     2773        # Small case, not really a problem ;-)
     2774        if self.order() == 1:
     2775            return copy(g)
     2776
     2777        # This is a sign that this is the first call
     2778        # to this recursive function
     2779        if vertices is None:
     2780            # Now is the time to care about positions
     2781            pos = self.get_pos()
     2782
     2783            # if the graph is not connected, returns the union
     2784            # of the Gomory-Hu tree of each component
     2785            if not self.is_connected():
     2786                g = Graph()
     2787                for cc in self.connected_components_subgraphs():
     2788                    g = g.union(cc._gomory_hu_tree())
     2789                g.set_pos(self.get_pos())
     2790                return g
     2791            # All the vertices is this graph are the "real ones"
     2792            vertices = Set(self.vertices())
     2793
     2794        # There may be many vertices, though only one which is "real"
     2795        if len(vertices) == 1:
     2796            g = Graph()
     2797            g.add_vertex(vertices[0])
     2798            return g
     2799
     2800        # Take any two vertices
     2801        u,v = vertices[0:2]
     2802
     2803        # flow = connectivity between u and v
     2804        # edges = min cut
     2805        # sets1, sets2 = the two sides of the edge cut
     2806        flow,edges,[set1,set2] = self.edge_cut(u, v, use_edge_labels=True, vertices=True)
     2807
     2808        # One graph for each part of the previous one
     2809        g1,g2 = self.subgraph(set1), self.subgraph(set2)
     2810
     2811        # Adding the fake vertex to each part
     2812        g1_v = Set(set2)
     2813        g2_v = Set(set1)
     2814        g1.add_vertex(g1_v)
     2815        g1.add_vertex(g2_v)
     2816
     2817        # Each part of the graph had many edges going to the other part
     2818        # Now that we have a new fake vertex in each part
     2819        # we just say that the edges which were in the cut and going
     2820        # to the other side are now going to this fake vertex
     2821
     2822        # We must preserve the labels. They sum.
     2823
     2824        for x,y in edges:
     2825            # Assumes x is in g1
     2826            if x in g2:
     2827                x,y = y,x
     2828            # If the edge x-g1_v exists, adds to its label the capacity of arc xy
     2829            if g1.has_edge(x, g1_v):
     2830                g1.set_edge_label(x, g1_v, g1.edge_label(x, g1_v) + capacity(self.edge_label(x, y)))
     2831            else:
     2832                # Otherwise, creates it with the good label
     2833                g1.add_edge(x, g1_v, capacity(self.edge_label(x, y)))
     2834            # Same thing for g2
     2835            if g2.has_edge(y, g2_v):
     2836                g2.set_edge_label(y, g2_v, g2.edge_label(y, g2_v) + capacity(self.edge_label(x, y)))
     2837            else:
     2838                g2.add_edge(y, g2_v, capacity(self.edge_label(x, y)))
     2839
     2840        # Recursion for the two new graphs... The new "real" vertices are the intersection with
     2841        # with the previous set of "real" vertices
     2842        g1_tree = g1._gomory_hu_tree(vertices=(vertices & Set(g1.vertices())))
     2843        g2_tree = g2._gomory_hu_tree(vertices=(vertices & Set(g2.vertices())))
     2844
     2845        # Union of the two partial trees ( it is disjoint, but
     2846        # disjoint_union does not preserve the name of the vertices )
     2847        g = g1_tree.union(g2_tree)
     2848
     2849        # An edge to connect them, with the appropriate label
     2850        g.add_edge(g1_tree.vertex_iterator().next(), g2_tree.vertex_iterator().next(), flow)
     2851
     2852        if pos:
     2853            g.set_pos(pos)
     2854
     2855        return g
     2856
     2857    def gomory_hu_tree(self):
     2858        r"""
     2859        Returns the Gomory-Hu tree associated to self.
     2860
     2861        Given a tree `T` with labelled edges representing capacities, it is very
     2862        easy to determine the maximal flow between any pair of vertices :
     2863        it is the minimal label on the edges of the unique path between them.
     2864
     2865        Given a graph `G`, the Gomory-Hu tree `T` of `G`, is a tree
     2866        with the same set of vertices, and such that the maximal flow
     2867        between any two vertices is the same in `G` and in `T`.
     2868        (see http://en.wikipedia.org/wiki/Gomory%E2%80%93Hu_tree)
     2869
     2870        OUTPUT:
     2871
     2872        graph with labeled edges
     2873
     2874        EXAMPLE:
     2875
     2876        Taking the Petersen graph::
     2877
     2878            sage: g = graphs.PetersenGraph()
     2879            sage: t = g.gomory_hu_tree() # optional - requires GLPK or CBC
     2880
     2881        Obviously, this graph is a tree::
     2882
     2883            sage: t.is_tree()  # optional - requires GLPK or CBC
     2884            True
     2885
     2886        Note that if the original graph is not connected, then the
     2887        Gomory-Hu tree is in fact a forest::
     2888
     2889            sage: (2*g).gomory_hu_tree().is_forest() # optional - requires GLPK or CBC
     2890            True
     2891            sage: (2*g).gomory_hu_tree().is_connected() # optional - requires GLPK or CBC
     2892            False
     2893
     2894        On the other hand, such a tree has lost nothing of the initial
     2895        graph connectedness::
     2896
     2897            sage: all([ t.flow(u,v) == g.flow(u,v) for u,v in Subsets( g.vertices(), 2 ) ]) # optional - requires GLPK or CBC
     2898            True
     2899
     2900        Just to make sure, let's check the same is true for two vertices
     2901        in a random graph::
     2902
     2903            sage: g = graphs.RandomGNP(20,.3)
     2904            sage: t = g.gomory_hu_tree() # optional - requires GLPK or CBC
     2905            sage: g.flow(0,1) == t.flow(0,1) # optional - requires GLPK or CBC
     2906            True
     2907
     2908        And also the min cut::
     2909
     2910            sage: g.edge_connectivity() == min(t.edge_labels()) # optional - requires GLPK or CBC
     2911            True
     2912        """
     2913        return self._gomory_hu_tree()
     2914       
     2915
    27372916    def two_factor_petersen(self):
    27382917        r"""
    27392918        Returns a decomposition of the graph into 2-factors.
  • sage/graphs/graph_generators.py

    diff -r deb9406f163f -r f28f7f29eaf0 sage/graphs/graph_generators.py
    a b  
    44084408            sage: g = graphs.DegreeSequenceBipartite([2,2,2,2,1],[5,5]) # optional - requires GLPK or CBC
    44094409            Traceback (most recent call last):
    44104410            ...
    4411             ValueError: The two partitions must sum to the same value.
    4412         """
    4413        
    4414         from sage.combinat.partition import Partition
     4411            ValueError: There exists no bipartite graph corresponding to the given degree sequences
     4412        """
     4413       
     4414        from sage.combinat.integer_vector import gale_ryser_theorem
    44154415        from sage.graphs.bipartite_graph import BipartiteGraph
    44164416
    4417         s1.sort(reverse=True)
    4418         s2.sort(reverse=True)
    4419        
    4420         p1 = Partition(s1)
    4421         p2 = Partition(s2)
    4422 
    4423         return BipartiteGraph(p1.gale_ryser_theorem(p2))
     4417        s1 = sorted(s1, reverse = True)
     4418        s2 = sorted(s2, reverse = True)
     4419
     4420        m = gale_ryser_theorem(s1,s2)
     4421       
     4422        if m is False:
     4423            raise ValueError("There exists no bipartite graph corresponding to the given degree sequences")
     4424        else:
     4425            return BipartiteGraph(m)
    44244426
    44254427    def DegreeSequenceConfigurationModel(self, deg_sequence, seed=None):
    44264428        """