Ticket #11279: trac_11279.patch

File trac_11279.patch, 21.8 KB (added by ncohen, 8 years ago)
  • doc/en/reference/graphs.rst

    # HG changeset patch
    # User Nathann Cohen <nathann.cohen@gmail.com>
    # Date 1304282108 -7200
    # Node ID 21fa710c545fdd76a157e1677a77d1f6e78e617d
    # Parent  e79daf00e92ab5227088289de977a30ae6beef7e
    trac 11279 -- Convexity in graphs
    
    diff --git a/doc/en/reference/graphs.rst b/doc/en/reference/graphs.rst
    a b  
    4848   sage/graphs/pq_trees
    4949   sage/graphs/matchpoly
    5050   sage/graphs/graph_decompositions/vertex_separation
     51   sage/graphs/convexity_properties
    5152   sage/graphs/graph_latex
    5253   sage/graphs/graph_list
    5354
  • module_list.py

    diff --git a/module_list.py b/module_list.py
    a b  
    322322    Extension('sage.graphs.graph_decompositions.vertex_separation',
    323323              sources = ['sage/graphs/graph_decompositions/vertex_separation.pyx']),
    324324
     325    Extension('sage.graphs.convexity_properties',
     326              sources = ['sage/graphs/convexity_properties.pyx']),
     327
    325328    Extension('sage.graphs.generic_graph_pyx',
    326329              sources = ['sage/graphs/generic_graph_pyx.pyx'],
    327330              libraries = ['gmp']),
  • new file sage/graphs/convexity_properties.pxd

    diff --git a/sage/graphs/convexity_properties.pxd b/sage/graphs/convexity_properties.pxd
    new file mode 100644
    - +  
     1include "../misc/bitset_pxd.pxi"
     2
     3cdef class ConvexityProperties:
     4    cdef int _n
     5    cdef list _list_integers_to_vertices
     6    cdef dict _dict_vertices_to_integers
     7    cdef bitset_t * _cache_hull_pairs
     8
     9    cdef list _vertices_to_integers(self, vertices)
     10    cdef list _integers_to_vertices(self, integers)
     11    cdef _bitset_convex_hull(self, bitset_t hull)
     12    cpdef hull(self, list vertices)
     13    cdef _greedy_increase(self, bitset_t bs)
     14    cpdef hull_number(self, value_only = *, verbose = *)
  • new file sage/graphs/convexity_properties.pyx

    diff --git a/sage/graphs/convexity_properties.pyx b/sage/graphs/convexity_properties.pyx
    new file mode 100644
    - +  
     1r"""
     2Convexity properties of graphs
     3
     4This class gathers the algorithms related to convexity in a graph.
     5
     6AUTHOR: Nathann Cohen
     7"""
     8
     9include "../misc/bitset.pxi"
     10from sage.numerical.backends.generic_backend cimport GenericBackend
     11from sage.numerical.backends.generic_backend import get_solver
     12
     13cdef class ConvexityProperties:
     14    r"""
     15    This class gathers the algorithms related to convexity in a graph.
     16
     17    **Definitions**
     18
     19    A set `S \subseteq V(G)` of vertices is said to be convex if for all `u,v\in
     20    S` the set `S` contains all the vertices located on a shortest path between
     21    `u` and `v`. Alternatively, a set `S` is said to be convex if the distances
     22    satisfy `\forall u,v\in S, \forall w\in V\backslash S : d_{G}(u,w) +
     23    d_{G}(w,v) > d_{G}(u,v)`.
     24
     25    The convex hull `h(S)` of a set `S` of vertices is defined as the smallest
     26    convex set containing `S`.
     27
     28    It is a closure operator, as trivially `S\subseteq h(S)` and `h(h(S)) =
     29    h(S)`.
     30
     31    **What this class contains**
     32
     33    As operations on convex sets generally involve the computation of distances
     34    between vertices, this class' purpose is to cache that information so that
     35    computing the convex hulls of several different sets of vertices does not
     36    imply recomputing several times the distances between the vertices.
     37
     38    In order to compute the convex hull of a set `S` it is possible to write the
     39    following algorithm.
     40
     41    *For any pair `u,v` of elements in the set `S`, and for any vertex `w`*
     42    *outside of it, add `w` to `S` if `d_{G}(u,w) + d_{G}(w,v) = d_{G}(u,v)`.*
     43    *When no vertex can be added anymore, the set `S` is convex*
     44
     45    The distances are not actually that relevant. The same algorithm can be
     46    implemented by remembering for each pair `u, v` of vertices the list of
     47    elements `w` satisfying the condition, and this is precisely what this class
     48    remembers, encoded as bitsets to make storage and union operations more
     49    efficient.
     50
     51    .. NOTE::
     52
     53        * This class is useful if you compute the convex hulls of many sets in
     54          the same graph, or if you want to compute the hull number itself as it
     55          involves many calls to :meth:`hull`
     56       
     57        * Using this class on non-conected graphs is a waste of space and
     58          efficiency ! If your graph is disconnected, the best for you is to
     59          deal independently with each connected component, whatever you are
     60          doing.
     61
     62    **Possible improvements**
     63
     64    When computing a convex set, all the pairs of elements belonging to the set
     65    `S` are enumerated several times.
     66
     67    * There should be a smart way to avoid enumerating pairs of vertices which
     68      have already been tested. The cost of each of them is not very high, so
     69      keeping track of those which have been tested already may be too expensive
     70      to gain any efficiency.
     71
     72    * The ordering in which they are visited is currently purely lexicographic,
     73      while there is a Poset structure to exploit. In particular, when two
     74      vertices `u, v` are far apart and generate a set `h(\{u,v\})` of vertices,
     75      all the pairs of vertices `u', v'\in h(\{u,v\})` satisfy `h(\{u',v'\})
     76      \subseteq h(\{u,v\})`, and so it is useless to test the pair `u', v'` when
     77      both `u` and `v` where present.
     78
     79    * The information cached is for any pair `u,v` of vertices the list of
     80      elements `z` with `d_{G}(u,w) + d_{G}(w,v) = d_{G}(u,v)`. This is not in
     81      general equal to `h(\{u,v\})` !
     82
     83    Nothing says these recommandations will actually lead to any actual
     84    improvements. There are just some ideas remembered while writing this
     85    code. Trying to optimize may well lead to lost in efficiency on many
     86    instances.
     87
     88    EXAMPLE::
     89
     90        sage: from sage.graphs.convexity_properties import ConvexityProperties
     91        sage: g = graphs.PetersenGraph()
     92        sage: CP = ConvexityProperties(g)
     93        sage: CP.hull([1,3])
     94        [1, 2, 3]
     95        sage: CP.hull_number()
     96        3
     97
     98    TESTS::
     99
     100        sage: ConvexityProperties(digraphs.Circuit(5))
     101        Traceback (most recent call last):
     102        ...
     103        ValueError: This is currenly implemented for Graphs only.Only minor updates are needed if you want to makeit support DiGraphs too.
     104    """
     105   
     106    def __init__(self, G):
     107        r"""
     108        Constructor
     109
     110        EXAMPLE::
     111
     112            sage: from sage.graphs.convexity_properties import ConvexityProperties
     113            sage: g = graphs.PetersenGraph()
     114            sage: ConvexityProperties(g)
     115            <sage.graphs.convexity_properties.ConvexityProperties object at ...>
     116        """
     117        from sage.graphs.digraph import DiGraph
     118        if isinstance(G, DiGraph):
     119            raise ValueError("This is currenly implemented for Graphs only."+
     120                             "Only minor updates are needed if you want to make"+
     121                             "it support DiGraphs too.")
     122       
     123        # Cached number of vertices
     124        cdef int n = G.order()
     125        self._n = n
     126
     127        cdef int i = 0
     128        cdef int j,k
     129
     130        # Temporary variables
     131        cdef dict d_i
     132        cdef dict d_j
     133        cdef int d_ij
     134        self._dict_vertices_to_integers = {}
     135        self._list_integers_to_vertices = []
     136
     137        # Remembering integers instead of the labels, and building dictionaries
     138        # in both directions.
     139        for v in G:
     140            self._dict_vertices_to_integers[v] = i
     141            self._list_integers_to_vertices.append(v)
     142            i = i + 1
     143
     144
     145        # Computation of distances between all pairs. Costly.
     146        cdef dict distances = G.distance_all_pairs()
     147
     148        # _cache_hull_pairs[u*n + v] is a bitset whose 1 bits are the vertices located on a shortest path from vertex u to v
     149        #
     150        # Note that  u < v
     151        self._cache_hull_pairs = <bitset_t *> sage_malloc(((n*(n-1))>>1)*sizeof(bitset_t))
     152        cdef bitset_t * p_bitset = self._cache_hull_pairs
     153
     154        # Filling the cache
     155        #
     156        # The p_bitset variable iterates over the successive elements of the cache
     157        #
     158        # For any pair i,j of vertices (i<j), we built the bitset of all the
     159        # elements k which are on a shortest path from i to j
     160
     161        for 0<= i < n-1:
     162            # Caching the distances from i to the other vertices
     163            d_i = distances[self._list_integers_to_vertices[i]]
     164
     165            for i < j < n:
     166                # Caching the distances from j to the other vertices
     167                d_j = distances[self._list_integers_to_vertices[j]]
     168
     169                # Caching the distance between i and j
     170                d_ij = d_i[self._list_integers_to_vertices[j]]
     171
     172                # Initializing the new bitset
     173                bitset_init(p_bitset[0], n)
     174                bitset_set_first_n(p_bitset[0], 0)
     175
     176                # Filling it
     177                for 0<= k < n:
     178                    if ((d_i[self._list_integers_to_vertices[k]]
     179                         + d_j[self._list_integers_to_vertices[k]])
     180                        == d_ij):
     181                        bitset_add(p_bitset[0], k)
     182
     183                # Next bitset !
     184                p_bitset = p_bitset + 1
     185
     186   
     187    def __destruct__(self):
     188        r"""
     189        Destructor
     190
     191        EXAMPLE::
     192
     193            sage: from sage.graphs.convexity_properties import ConvexityProperties
     194            sage: g = graphs.PetersenGraph()
     195            sage: ConvexityProperties(g)
     196            <sage.graphs.convexity_properties.ConvexityProperties object at ...>
     197
     198        """
     199        cdef bitset_t * p_bitset = self._cache_hull_pairs
     200        cdef int i
     201
     202        for 0 <= i < ((self._n*(self._n-1))>>1):
     203            bitset_free(p_bitset[0])
     204            p_bitset = p_bitset + 1
     205
     206        sage_free(self._cache_hull_pairs)
     207           
     208    cdef list _vertices_to_integers(self, vertices):
     209        r"""
     210        Converts a list of vertices to a list of integers with the cached data.
     211        """
     212        cdef list answer = []
     213        for v in v:
     214            answer.append(self._dict_vertices_to_integers[v])
     215        return answer
     216
     217    cdef list _integers_to_vertices(self, integers):
     218        r"""
     219        Converts a list of integers to a list of vertices with the cached data.
     220        """
     221
     222        cdef list answer = []
     223        for v in integers:
     224            answer.append(self._list_integers_to_vertices[v])
     225        return answer
     226
     227    cdef _bitset_convex_hull(self, bitset_t hull):
     228        r"""
     229        Computes the convex hull of a list of vertices given as a bitset.
     230
     231        (this method returns nothing and modifies the input)
     232        """
     233        cdef int count
     234        cdef int tmp_count
     235        cdef int i,j
     236
     237        cdef bitset_t * p_bitset
     238
     239        # Current size of the set
     240        count = bitset_len(hull)
     241
     242        while True:
     243           
     244            # Iterating over all the elements in the cache
     245            p_bitset = self._cache_hull_pairs
     246
     247            # For any vertex i
     248            for 0 <= i < self._n-1:
     249
     250                # If i is not in the current set, we skip it !
     251                if not bitset_in(hull, i):
     252                    p_bitset = p_bitset + (self._n-1-i)
     253                    continue
     254
     255                # If it is, we iterate over all the elements j
     256                for i < j < self._n:
     257
     258                    # If both i and j are inside, we add all the (cached)
     259                    # vertices on a shortest ij-path
     260
     261                    if bitset_in(hull, j):
     262                        bitset_union(hull, hull, p_bitset[0])
     263
     264                    # Next bitset !
     265                    p_bitset = p_bitset + 1
     266
     267           
     268            tmp_count = bitset_len(hull)
     269
     270            # If we added nothing new during the previous loop, our set is
     271            # convex !
     272            if tmp_count == count:
     273                return
     274   
     275            # Otherwise, update and back to the loop
     276            count = tmp_count
     277
     278    cpdef hull(self, list vertices):
     279        r"""
     280        Returns the convex hull of a set of vertices.
     281
     282        INPUT:
     283
     284        * ``vertices`` -- A list of vertices.
     285
     286        EXAMPLE::
     287
     288            sage: from sage.graphs.convexity_properties import ConvexityProperties
     289            sage: g = graphs.PetersenGraph()
     290            sage: CP = ConvexityProperties(g)
     291            sage: CP.hull([1,3])
     292            [1, 2, 3]
     293        """
     294        cdef bitset_t bs
     295        bitset_init(bs, self._n)
     296        bitset_set_first_n(bs, 0)
     297
     298        for v in vertices:
     299            bitset_add(bs, self._dict_vertices_to_integers[v])
     300
     301        self._bitset_convex_hull(bs)
     302
     303        #cdef list answer = bitset_list(bs)
     304        cdef list answer = self._integers_to_vertices(bitset_list(bs))
     305
     306        bitset_free(bs)
     307
     308        return answer
     309
     310    cdef _greedy_increase(self, bitset_t bs):
     311        r"""
     312        Given a bitset whose hull is not the whole set, greedily add vertices
     313        and stop before its hull is the whole set.
     314
     315        NOTE:
     316
     317        * Counting the bits at each turn is not the best way...
     318        """
     319        cdef bitset_t tmp
     320        bitset_init(tmp, self._n)
     321
     322
     323        for 0<= i < self._n:
     324            if not bitset_in(bs, i):
     325                bitset_copy(tmp, bs)
     326                bitset_add(tmp, i)
     327                self._bitset_convex_hull(tmp)
     328                if bitset_len(tmp) < self._n:
     329                    bitset_add(bs, i)
     330   
     331
     332    cpdef hull_number(self, value_only = True, verbose = False):
     333        r"""
     334        Computes the hull number and a correspondin generating set.
     335
     336        The hull number `hn(G)` of a graph `G` is the cardinality of a smallest
     337        set of vertices `S` such that `h(S)=V(G)`.
     338
     339        INPUT:
     340
     341        * ``value_only`` (boolean) -- whether to return only the hull number (default) or
     342          a minimum set whose convex hull is the whole graph.
     343
     344        * ``verbose`` (boolean) -- whether to display information on the LP.
     345
     346        **COMPLEXITY:**
     347
     348        This problem is NP-Hard [CHZ02]_, but seems to be of the "nice"
     349        kind. Update this comment if you fall on hard instances `:-)`
     350
     351        **ALGORITHM:**
     352
     353        This is solved by linear programming.
     354
     355        As the function `h(S)` associating to each set `S` its convex hull is a
     356        closure operator, it is clear that any set `S_G` of vertices such that
     357        `h(S_G)=V(G)` must satisfy `S_G \not \subseteq C` for any *proper*
     358        convex set `C \subsetneq V(G)`. The following formulation is hence
     359        correct
     360
     361        .. MATH::
     362
     363            \text{Minimize :}& \sum_{v\in G}b_v\\
     364            \text{Such that :}&\\
     365            &\forall C\subsetneq V(G)\text{ a proper convex set }\\
     366            &\sum_{v\in V(G)\backslash C} b_v \geq 1
     367
     368        Of course, the number of convex sets -- and so the number of constraints
     369        -- can be huge, and hard to enumerate, so at first an incomplete
     370        formulation is solved (it is missing some constraints). If the answer
     371        returned by the LP solver is a set `S` generating the whole graph, then
     372        it is optimal and so is returned. Otherwise, the constraint
     373        corresponding to the set `h(S)` can be added to the LP, which makes the
     374        answer `S` infeasible, and another solution computed.
     375
     376        This being said, simply adding the constraint corresponding to `h(S)` is
     377        a bit slow, as these sets can be large (and the corresponding constrait
     378        a bit weak). To improve it a bit, before being added, the set `h(S)` is
     379        "greedily enriched" to a set `S'` with vertices for as long as
     380        `h(S')\neq V(G)`. This way, we obtain a set `S'` with `h(S)\subseteq
     381        h(S')\subsetneq V(G)`, and the constraint corresponding to `h(S')` --
     382        which is stronger than the one corresponding to `h(S)` -- is added.
     383
     384        This can actually be seen as a hitting set problem on the complement of
     385        convex sets.
     386
     387        EXAMPLE:
     388
     389        The Hull number of Petersen's graph::
     390
     391            sage: from sage.graphs.convexity_properties import ConvexityProperties
     392            sage: g = graphs.PetersenGraph()
     393            sage: CP = ConvexityProperties(g)
     394            sage: CP.hull_number()
     395            3
     396            sage: generating_set = CP.hull_number(value_only = False)
     397            sage: CP.hull(generating_set)
     398            [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
     399
     400        REFERENCE:
     401
     402        .. [CHZ02] F. Harary, E. Loukakis, C. Tsouros
     403          The geodetic number of a graph
     404          Mathematical and computer modelling
     405          vol. 17 n11 pp.89--95, 1993
     406        """
     407   
     408        cdef int i
     409        cdef list constraint # temporary variable to add constraints to the LP
     410
     411        if self._n <= 2:
     412            if value_only:
     413                return self._n
     414            else:
     415                return self._list_integers_to_vertices
     416           
     417        cdef GenericBackend p = <GenericBackend> get_solver(constraint_generation = True)
     418
     419        # Minimization
     420        p.set_sense(False)
     421
     422        # We have exactly n binary variables, all of them with a coefficient of
     423        # 1 in the objective function
     424        p.add_variables(self._n, 0, None, True, False, False, 1, None)
     425
     426        # We know that at least 2 vertices are required to cover the whole graph
     427        p.add_linear_constraint(zip(range(self._n), [1]*self._n), 2, None)
     428
     429        # The set of vertices generated by the current LP solution
     430        cdef bitset_t current_hull
     431        bitset_init(current_hull, self._n)
     432
     433        # Which is at first empty
     434        bitset_set_first_n(current_hull,1)
     435
     436        while True:
     437           
     438            # Greedily increase it to obtain a better constraint
     439            self._greedy_increase(current_hull)
     440
     441            if verbose:
     442                print "Adding a constraint corresponding to convex set ",
     443                print bitset_list(current_hull)
     444
     445            # Building the corresponding constraint
     446            constraint = []
     447            for 0 <= i < self._n:
     448                if not bitset_in(current_hull, i):
     449                    constraint.append((i,1))
     450
     451            p.add_linear_constraint(constraint, 1, None)
     452
     453            p.solve()
     454
     455            # Computing the current solution's convex hull
     456            bitset_set_first_n(current_hull,0)
     457
     458            for 0 <= i < self._n:
     459                if p.get_variable_value(i) > .5:
     460                    bitset_add(current_hull, i)
     461
     462            self._bitset_convex_hull(current_hull)
     463
     464            # Are we done ?
     465            if bitset_len(current_hull) == self._n:
     466                break
     467
     468        bitset_free(current_hull)
     469
     470        if value_only:
     471            return <int> p.get_objective_value()
     472
     473        constraint = []
     474        for 0 <= i < self._n:
     475            if p.get_variable_value(i) > .5:
     476                constraint.append(i)
     477
     478        return self._integers_to_vertices(constraint)
  • sage/graphs/graph.py

    diff --git a/sage/graphs/graph.py b/sage/graphs/graph.py
    a b  
    30313031        from matchpoly import matching_polynomial
    30323032        return matching_polynomial(self, complement=complement, name=name)
    30333033
     3034    ### Convexity
     3035
     3036    def convexity_properties(self):
     3037        r"""
     3038        Returns a ``ConvexityProperties`` objet corresponding to ``self``.
     3039
     3040        This object contains the methods related to convexity in graphs (convex
     3041        hull, hull number) and caches useful information so that it becomes
     3042        comparatively cheaper to compute the convex hull of many different sets
     3043        of the same graph.
     3044
     3045        .. SEEALSO::
     3046
     3047            In order to know what can be done through this object, please refer
     3048            to module :mod:`sage.graphs.convexity_properties`.
     3049
     3050        .. NOTE::
     3051
     3052            If you want to compute many convex hulls, keep this object in memory
     3053            ! When it is created, it builds a table of useful information to
     3054            compute convex hulls. As a result ::
     3055           
     3056                sage: g = graphs.PetersenGraph()
     3057                sage: g.convexity_properties().hull([1, 3])
     3058                [1, 2, 3]
     3059                sage: g.convexity_properties().hull([3, 7])
     3060                [2, 3, 7]
     3061
     3062            Is a terrible waste of computations, while ::
     3063
     3064                sage: g = graphs.PetersenGraph()
     3065                sage: CP = g.convexity_properties()
     3066                sage: CP.hull([1, 3])
     3067                [1, 2, 3]
     3068                sage: CP.hull([3, 7])
     3069                [2, 3, 7]
     3070
     3071            Makes perfect sense.
     3072        """
     3073        from sage.graphs.convexity_properties import ConvexityProperties
     3074        return ConvexityProperties(self)
     3075
    30343076    ### Centrality
    30353077
    30363078    def centrality_betweenness(self, normalized=True):
  • sage/numerical/backends/generic_backend.pxd

    diff --git a/sage/numerical/backends/generic_backend.pxd b/sage/numerical/backends/generic_backend.pxd
    a b  
    2828    cpdef bint is_variable_binary(self, int)
    2929    cpdef bint is_variable_integer(self, int)
    3030    cpdef bint is_variable_continuous(self, int)
    31 
    3231    cpdef problem_name(self, char * name = *)
    3332    cpdef row_bounds(self, int index)
    3433    cpdef col_bounds(self, int index)
     
    3635    cpdef col_name(self, int index)
    3736    cpdef variable_upper_bound(self, int index, value = *)
    3837    cpdef variable_lower_bound(self, int index, value = *)
     38
     39cpdef GenericBackend get_solver(constraint_generation = ?, solver = ?)
  • sage/numerical/backends/glpk_backend.pyx

    diff --git a/sage/numerical/backends/glpk_backend.pyx b/sage/numerical/backends/glpk_backend.pyx
    a b  
    179179                glp_set_col_kind(self.lp, n_var - i, GLP_IV)
    180180
    181181            if obj:
    182                 self.bjective_coefficient(n_var - i - 1, obj)
     182                self.objective_coefficient(n_var - i - 1, obj)
    183183
    184184            if names is not None:
    185185                glp_set_col_name(self.lp, n_var, names[number - i - 1])