Trac #17540: Poset.dimension

With this code, Sage at last can compute the dimension of a Poset.

http://en.wikipedia.org/wiki/Order_dimension

Nathann

See https://groups.google.com/d/topic/sage-devel/qEE9YUAPA4g/discussion

URL: http://trac.sagemath.org/17540
Reported by: ncohen
Ticket author(s): Nathann Cohen
Reviewer(s): Dima Pasechnik

src/sage/combinat/posets/posets.py

@@ -40,6 +40,7 @@
     :meth:`~FinitePoset.coxeter_transformation` | Returns the matrix of the Auslander-Reiten translation acting on the Grothendieck group of the derived category of modules.
     :meth:`~FinitePoset.cuts` | Returns the cuts of the given poset.
     :meth:`~FinitePoset.dilworth_decomposition` | Returns a partition of the points into the minimal number of chains.
+    :meth:`~FinitePoset.dimension` | Return the dimension of the poset
     :meth:`~FinitePoset.disjoint_union` | Return the disjoint union of the poset with ``other``.
     :meth:`~FinitePoset.dual` | Returns the dual poset of the given poset.
     :meth:`~FinitePoset.evacuation` | Computes evacuation on the linear extension associated to the poset ``self``.
@@ -2404,6 +2405,185 @@ def is_EL_labelling(self, f, return_raising_chains=False):
         else:
             return True
 
+    def dimension(self, certificate=False):
+        r"""
+        Return the dimension of the Poset.
+
+        The (Dushnik-Miller) dimension of a Poset defined on a set `X` of points
+        is the smallest integer `n` such that there exists `P_1,...,P_n` linear
+        extensions of `P` satisfying the following property:
+
+        .. MATH::
+
+            u\leq_P v\ \text{if and only if }\ \forall i, u\leq_{P_i} v
+
+        For more information, see the :wikipedia:`Order_dimension`.
+
+        INPUT:
+
+        - ``certificate`` (boolean; default:``False``) -- whether to return an
+          integer (the dimension) or a certificate, i.e. a smallest set of
+          linear extensions.
+
+        .. NOTE::
+
+            The speed of this function greatly improves when more efficient MILP
+            solvers (e.g. Gurobi, CPLEX) are installed. See
+            :class:`MixedIntegerLinearProgram` for more information.
+
+        **Algorithm:**
+
+        As explained [FT00]_, the dimension of a poset is equal to the (weak)
+        chromatic number of a hypergraph. More precisely:
+
+            Let `inc(P)` be the set of (ordered) pairs of incomparable elements
+            of `P`, i.e. all `uv` and `vu` such that `u\not \leq_P v` and `v\not
+            \leq_P u`. Any linear extension of `P` is a total order on `X` that
+            can be seen as the union of relations from `P` along with some
+            relations from `inc(P)`. Thus, the dimension of `P` is the smallest
+            number of linear extensions of `P` which *cover* all points of
+            `inc(P)`.
+
+            Consequently, `dim(P)` is equal to the chromatic number of the
+            hypergraph `\mathcal H_{inc}`, where `\mathcal H_{inc}` is the
+            hypergraph defined on `inc(P)` whose sets are all `S\subseteq
+            inc(P)` such that `P\cup S` is not acyclic.
+
+        We solve this problem through a :mod:`Mixed Integer Linear Program
+        <sage.numerical.mip>`.
+
+        EXAMPLES:
+
+        According to Wikipedia, the poset (of height 2) of a graph is `\leq 2`
+        if and only if the graph is a path::
+
+            sage: G = graphs.PathGraph(6)
+            sage: P = Poset(DiGraph({(u,v):[u,v] for u,v,_ in G.edges()}))
+            sage: P.dimension()
+            2
+
+        The actual linear extensions can be obtained with ``certificate=True``::
+
+            sage: P.dimension(certificates=True) # not tested -- architecture-dependent
+            [[(0, 1), 0, (1, 2), 1, (2, 3), 2, (3, 4), 3, (4, 5), 4, 5],
+            [(4, 5), 5, (3, 4), 4, (2, 3), 3, (1, 2), 2, (0, 1), 1, 0]]
+
+        According to Schnyder's theorem, the poset (of height 2) of a graph has
+        dimension `\leq 3` if and only if the graph is planar::
+
+            sage: G = graphs.CompleteGraph(4)
+            sage: P = Poset(DiGraph({(u,v):[u,v] for u,v,_ in G.edges()}))
+            sage: P.dimension()
+            3
+
+            sage: G = graphs.CompleteBipartiteGraph(3,3)
+            sage: P = Poset(DiGraph({(u,v):[u,v] for u,v,_ in G.edges()}))
+            sage: P.dimension() # not tested - around 4s with CPLEX
+            4
+
+        TESTS:
+
+        Empty Poset::
+
+            sage: Poset().dimension()
+            0
+            sage: Poset().dimension(certificate=1)
+            []
+
+        References:
+
+        .. [FT00] Stefan Felsner, William T. Trotter,
+           Dimension, Graph and Hypergraph Coloring,
+           Order,
+           2000, Volume 17, Issue 2, pp 167-177,
+           http://link.springer.com/article/10.1023%2FA%3A1006429830221
+        """
+        if self.cardinality() == 0:
+            return [] if certificate else 0
+
+        from sage.numerical.mip import MixedIntegerLinearProgram, MIPSolverException
+        P = Poset(self._hasse_diagram) # work on an int-labelled poset
+        hasse_diagram = P.hasse_diagram()
+        inc_graph = P.incomparability_graph()
+        inc_P = inc_graph.edges(labels=False)
+
+        # Current bound on the chromatic number of the hypergraph
+        k = 1
+
+        # cycles is the list of all cycles found during the execution of the
+        # algorithm
+
+        cycles = [[(u,v),(v,u)] for u,v in inc_P]
+
+        def init_LP(k,cycles,inc_P):
+            r"""
+            Initializes a LP object with k colors and the constraints from 'cycles'
+
+                sage: init_LP(1,2,3) # not tested
+            """
+            p = MixedIntegerLinearProgram(constraint_generation=True)
+            b = p.new_variable(binary=True)
+            for (u,v) in inc_P: # Each point has a color
+                p.add_constraint(p.sum(b[(u,v),i] for i in range(k))==1)
+                p.add_constraint(p.sum(b[(v,u),i] for i in range(k))==1)
+            for cycle in cycles: # No monochromatic set
+                for i in range(k):
+                    p.add_constraint(p.sum(b[point,i] for point in cycle)<=len(cycle)-1)
+            return p,b
+
+        p,b = init_LP(k,cycles,inc_P)
+
+        while True:
+            # Compute a coloring of the hypergraph. If there is a problem,
+            # increase the number of colors and start again.
+            try:
+                p.solve()
+            except MIPSolverException:
+                k += 1
+                p,b = init_LP(k,cycles,inc_P)
+                continue
+
+            # We create the digraphs of all color classes
+            linear_extensions = [hasse_diagram.copy() for i in range(k)]
+            for ((u,v),i),x in p.get_values(b).iteritems():
+                if x == 1:
+                    linear_extensions[i].add_edge(u,v)
+
+            # We check that all color classes induce an acyclic graph, and add a
+            # constraint otherwise.
+            okay = True
+            for g in linear_extensions:
+                is_acyclic, cycle = g.is_directed_acyclic(certificate=True)
+                if not is_acyclic:
+                    okay = False # one is not acyclic
+                    cycle = [(cycle[i-1],cycle[i]) for i in range(len(cycle))]
+                    cycle = [(u,v) for u,v in cycle if not P.lt(u,v) and not P.lt(v,u)]
+                    cycles.append(cycle)
+                    for i in range(k):
+                        p.add_constraint(p.sum(b[point,i] for point in cycle)<=len(cycle)-1)
+            if okay:
+                break
+
+        linear_extensions = [g.topological_sort() for g in linear_extensions]
+
+        # Check that the linear extensions do generate the poset (just to be
+        # sure)
+        from itertools import combinations
+        n = P.cardinality()
+        d = DiGraph()
+        for l in linear_extensions:
+            d.add_edges(combinations(l,2))
+
+        # The only 2-cycles are the incomparable pair
+        if d.size() != (n*(n-1))/2+inc_graph.size():
+            raise RuntimeError("Something went wrong. Please report this "
+                               "bug to sage-devel@googlegroups.com")
+
+        if certificate:
+            return [[self._list[i] for i in l]
+                    for l in linear_extensions]
+        return k
+
     def rank_function(self):
         r"""
         Return the (normalized) rank function of the poset,