Ticket #8894: trac_8894.patch

File trac_8894.patch, 8.6 KB (added by ncohen, 10 years ago)
  • sage/graphs/graph.py

    # HG changeset patch
    # User Nathann Cohen <nathann.cohen@gmail.com>
    # Date 1284856788 -7200
    # Node ID 44bebdf79947d9472ea151aaba549bae33bbff61
    # Parent  4a9f7b41ae22e247c2a6675c7854e419f7340d09
    trac 8894 -- Graph.topological_minor
    
    diff -r 4a9f7b41ae22 -r 44bebdf79947 sage/graphs/graph.py
    a b  
    29582958        f.write( print_graph_eps(self.vertices(), self.edge_iterator(), pos) )
    29592959        f.close()
    29602960
     2961    def topological_minor(self, H, vertices = False, paths = False, solver=None, verbose=0):
     2962        r"""
     2963        Returns a topological `H`-minor from ``self`` if one exists.
     2964
     2965        We say that a graph `G` has a topological `H`-minor (or that
     2966        it has a graph isomorphic to `H` as a topological minor), if
     2967        `G` contains a subdivision of a graph isomorphic to `H` (=
     2968        obtained from `H` through arbitrary subdivision of its edges)
     2969        as a subgraph.
     2970
     2971        For more information, see the `Wikipedia article on graph
     2972        minor
     2973        <http://en.wikipedia.org/wiki/Minor_%28graph_theory%29>`_.
     2974
     2975        INPUT:
     2976
     2977        - ``H`` -- The topological minor to find in the current graph.
     2978
     2979        - ``solver`` -- (default: ``None``) Specify a Linear Program (LP)
     2980          solver to be used. If set to ``None``, the default one is used. For
     2981          more information on LP solvers and which default solver is used, see
     2982          the method
     2983          :meth:`solve <sage.numerical.mip.MixedIntegerLinearProgram.solve>`
     2984          of the class
     2985          :class:`MixedIntegerLinearProgram <sage.numerical.mip.MixedIntegerLinearProgram>`.
     2986
     2987        - ``verbose`` -- integer (default: ``0``). Sets the level of
     2988          verbosity. Set to 0 by default, which means quiet.
     2989
     2990        OUTPUT:
     2991
     2992        The topological `H`-minor found is returned as a subgraph `M`
     2993        of ``self``, such that the vertex `v` of `M` that represents a
     2994        vertex `h\in H` has ``h`` as a label (see
     2995        :meth:`get_vertex <sage.graphs.generic_graph.get_vertex>`
     2996        and
     2997        :meth:`set_vertex <sage.graphs.generic_graph.set_vertex>`),
     2998        and such that every edge of `M` has as a label the edge of `H`
     2999        it (partially) represents.
     3000
     3001        If no topological minor is found, this method returns
     3002        ``False``.
     3003
     3004        ALGORITHM:
     3005
     3006        Mixed Integer Linear Programming.
     3007
     3008        COMPLEXITY:
     3009
     3010        Theoretically, when `H` is fixed, testing for the existence of
     3011        a topological `H`-minor is polynomial. The known algorithms
     3012        are highly exponential in `H`, though.
     3013
     3014        .. NOTE::
     3015       
     3016            * This function can be expected to be *very* slow, especially
     3017              where the topological minor does not exist.
     3018
     3019            (CPLEX seems to be *much* more efficient than GLPK on this
     3020            kind of problem)
     3021
     3022        EXAMPLES:
     3023
     3024        Petersen's graph has a topological `K_4`-minor::
     3025
     3026            sage: g = graphs.PetersenGraph()
     3027            sage: g.topological_minor(graphs.CompleteGraph(4))
     3028            Subgraph of (Petersen graph): Graph on ...
     3029
     3030        And a topological `K_{3,3}`-minor::
     3031
     3032            sage: g.topological_minor(graphs.CompleteBipartiteGraph(3,3))
     3033            Subgraph of (Petersen graph): Graph on ...
     3034
     3035        And of course, a tree has no topological `C_3`-minor::
     3036
     3037            sage: g = graphs.RandomGNP(15,.3)
     3038            sage: g = g.subgraph(edges = g.min_spanning_tree())
     3039            sage: g.topological_minor(graphs.CycleGraph(3))
     3040            False
     3041        """
     3042
     3043        # Useful alias ...
     3044        G = self
     3045
     3046        from sage.numerical.mip import MixedIntegerLinearProgram, Sum, MIPSolverException
     3047        p = MixedIntegerLinearProgram()
     3048
     3049        # This is an existence problem
     3050        p.set_objective(None)
     3051
     3052        #######################
     3053        # Vertex representant #
     3054        #######################
     3055        #
     3056        # v_repr[h][g] = 1 if vertex h from H is represented by vertex
     3057        # g from G, 0 otherwise
     3058
     3059        v_repr = p.new_variable(binary = True, dim = 2)
     3060
     3061        # Exactly one representant per vertex of H
     3062        for h in H:
     3063            p.add_constraint( Sum( v_repr[h][g] for g in G), min = 1, max = 1)
     3064
     3065        # A vertex of G can only represent one vertex of H
     3066        for g in G:
     3067            p.add_constraint( Sum( v_repr[h][g] for h in H), max = 1)
     3068
     3069
     3070        ###################
     3071        # Is representent #
     3072        ###################
     3073        #
     3074        # is_repr[v] = 1 if v represents some vertex of H
     3075       
     3076        is_repr = p.new_variable(binary = True)
     3077
     3078        for g in G:
     3079            for h in H:
     3080                p.add_constraint( v_repr[h][g] - is_repr[g], max = 0)
     3081
     3082
     3083        ###################################
     3084        # paths between the representents #
     3085        ###################################
     3086        #
     3087        # For any edge (h1,h2) in H, we have a corresponding path in G
     3088        # between the representants of h1 and h2. Which means there is
     3089        # a flow of intensity 1 from one to the other.
     3090        # We are then writing a flow problem for each edge of H.
     3091        #
     3092        # The variable flow[(h1,h2)][(g1,g2)] indicates the amount of
     3093        # flow on the edge (g1,g2) representing the edge (h1,h2).
     3094
     3095        flow = p.new_variable(binary = True, dim = 2)
     3096
     3097        # This lambda function returns the balance of flow
     3098        # corresponding to commodity C at vertex v v
     3099
     3100        flow_in = lambda C, v : Sum( flow[C][(v,u)] for u in G.neighbors(v) )
     3101        flow_out = lambda C, v : Sum( flow[C][(u,v)] for u in G.neighbors(v) )
     3102
     3103        flow_balance = lambda C, v : flow_in(C,v) - flow_out(C,v)
     3104
     3105        for h1,h2 in H.edges(labels = False):
     3106           
     3107            for v in G:
     3108               
     3109                # The flow balance depends on whether the vertex v is
     3110                # a representant of h1 or h2 in G, or a reprensentant
     3111                # of none
     3112
     3113                p.add_constraint( flow_balance((h1,h2),v) == v_repr[h1][v] - v_repr[h2][v] )
     3114
     3115        #############################
     3116        # Internal vertex of a path #
     3117        #############################
     3118        #
     3119        # is_internal[C][g] = 1 if a vertex v from G is located on the
     3120        # path representing the edge (=commodity) C
     3121
     3122        is_internal = p.new_variable(dim = 2, binary = True)
     3123
     3124        # When is a vertex internal for a commodity ?
     3125        for C in H.edges(labels = False):
     3126            for g in G:
     3127                p.add_constraint( flow_in(C,g) + flow_out(C,g) - is_internal[C][g], max = 1)
     3128
     3129
     3130        ############################
     3131        # Two paths do not cross ! #
     3132        ############################
     3133
     3134        # A vertex can only be internal for one commodity, and zero if
     3135        # the vertex is a representent
     3136
     3137        for g in G:
     3138            p.add_constraint( Sum( is_internal[C][g] for C in H.edges(labels = False))
     3139                              + is_repr[g], max = 1 )
     3140
     3141        # (The following inequalities are not necessary, but they seem
     3142        # to be of help (the solvers find the answer quicker when they
     3143        # are added)
     3144
     3145        # The flow on one edge can go in only one direction. Besides,
     3146        # it can belong to at most one commodity and has a maximum
     3147        # intensity of 1.
     3148
     3149        for g1,g2 in G.edges(labels = None):
     3150           
     3151            p.add_constraint(   Sum( flow[C][(g1,g2)] for C in H.edges(labels = False) )
     3152                              + Sum( flow[C][(g2,g1)] for C in H.edges(labels = False) ),
     3153                                max = 1)
     3154
     3155
     3156        # Now we can solve the problem itself !
     3157
     3158        try:
     3159            p.solve(solver = solver, log = verbose)
     3160
     3161        except MIPSolverException:
     3162            return False
     3163
     3164
     3165        minor = G.subgraph()
     3166
     3167        is_repr = p.get_values(is_repr)
     3168        v_repr = p.get_values(v_repr)
     3169        flow = p.get_values(flow)
     3170
     3171       
     3172        for u,v in minor.edges(labels = False):
     3173            used = False
     3174            for C in H.edges(labels = False):
     3175
     3176                if flow[C][(u,v)] + flow[C][(v,u)] > .5:
     3177                    used = True
     3178                    minor.set_edge_label(u,v,C)
     3179                    break
     3180            if not used:
     3181                minor.delete_edge(u,v)
     3182
     3183        minor.delete_vertices( [v for v in minor
     3184                                if minor.degree(v) == 0 ] )
     3185
     3186        for g in minor:
     3187            if is_repr[g] > .5:
     3188                for h in H:
     3189                    if v_repr[h][v] > .5:
     3190                        minor.set_vertex(g,h)
     3191                        break
     3192
     3193        return minor
     3194
    29613195
    29623196    ### Cliques
    29633197