Ticket #12716: trac_12716_MILP.patch

sage/graphs/graph_decompositions/vertex_separation.pyx

@@ -1 +1 @@
 r"""
 Vertex separation
 
-This module implements several algorithms to compute the vertex
-separation of a digraph and the corresponding ordering of the
-vertices. Given an ordering `v_1, ..., v_n` of the vertices of `V(G)`,
-its *cost* is defined as:
+This module implements several algorithms to compute the vertex separation of a
+digraph and the corresponding ordering of the vertices. It also implements tests
+functions for evaluation the width of a linear orderings.
+
+Given an ordering
+`v_1,\cdots, v_n` of the vertices of `V(G)`, its *cost* is defined as:
 
 .. MATH::
 
@@ -16 +18 @@
 
     c'(S) = |N^+_G(S)\backslash S|
 
-The *vertex separation* of a digraph `G` is equal to the minimum cost
-of an ordering of its vertices.
+The *vertex separation* of a digraph `G` is equal to the minimum cost of an
+ordering of its vertices.
 
 **Vertex separation and pathwidth**
 
@@ -25 +27 @@
 `G` a digraph `D` with the same vertex set, and in which each edge `uv` of `G`
 is replaced by two edges `uv` and `vu` in `D`. The vertex separation of `D` is
 equal to the pathwidth of `G`, and the corresponding ordering of the vertices of
-`D` encodes an optimal path-decomposition of `G`.
-
+`D`, also called a *layout*, encodes an optimal path-decomposition of `G`.
 This is a result of Kinnersley [Kin92]_ and Bodlaender [Bod98]_.
 
-**References**
 
-.. [Kin92] *The vertex separation number of a graph equals its path-width*,
-  Nancy G. Kinnersley,
-  Information Processing Letters,
-  Volume 42, Issue 6, Pages 345-350,
-  24 July 1992
+**This module contains the following methods**
 
-.. [Bod98] *A partial k-arboretum of graphs with bounded treewidth*,
-  Hans L. Bodlaender,
-  Theoretical Computer Science,
-  Volume 209, Issues 1-2, Pages 1-45,
-  6 December 1998
+===================================  ===========================================
 
-**Authors**
+:meth:`path_decomposition`           Returns the pathwidth of the given graph
+                                     and the ordering of the vertices resulting
+                                     in a corresponding path decomposition
 
-    - Nathann Cohen
+:meth:`vertex_separation`            Returns an optimal ordering of the vertices
+                                     and its cost for vertex-separation
+
+:meth:`vertex_separation_MILP`       Computes the vertex separation of `G` and
+                                     the optimal ordering of its vertices using
+                                     an MILP formulation
+
+:meth:`lower_bound`                  Returns a lower bound on the vertex
+                                     separation of `G`
+
+:meth:`is_valid_ordering`            Test if the linear vertex ordering `L` is
+                                     valid for (di)graph `G`
+
+:meth:`width_of_path_decomposition`  Returns the width of the path decomposition
+                                     induced by the linear ordering `L` of the
+                                     vertices of `G`
+===================================  ===========================================
+
 
 Exponential algorithm for vertex separation
 -------------------------------------------
@@ -114 +125 @@
 
     \max c'(\{v_1\}),c'(\{v_1,v_2\}),...,c'(\{v_1,...,v_n\})\geq\max c'_1,...,c'_n
 
-where `c_i` is the minimum cost of a set `S` on `i`
-vertices. Evaluating the `c_i` can take time (and in particular more
-than the previous exact algorithm), but it does not need much memory
-to run.
+where `c_i` is the minimum cost of a set `S` on `i` vertices. Evaluating the
+`c_i` can take time (and in particular more than the previous exact algorithm),
+but it does not need much memory to run.
 
-Methods
+
+MILP formulation for the vertex separation
+------------------------------------------
+
+We describe bellow a mixed integer linear program (MILP) for determining an
+optimal layout for the vertex separation of `G`. This MILP is an improved
+version of the formulation proposed in [SP10]_.
+
+VARIABLES:
+
+- ``x_v^t`` -- Variable set to 1 if either `v` has an in-neighbor `u` such that
+  `u\in S_{t'}=\{v_1, v_2,\cdots, v_{t'}\}` and `v\in N_G^+(S_{t'})\backslash
+  S_{t'}` with `t'\leq t`, or `v\in S_t`. It is set to 0 otherwise.
+
+- ``y_v^t`` -- Variable set to 1 if `v\in S_t`, and 0 otherwise. The order of
+  `v` in the layout is the smallest `t` such that `y_v^t==1`.
+
+- ``u_v^t`` -- Variable set to 1 if `v\in N_G^+(S_t)\backslash S_t` has an
+  in-neighbor in `S_t`, and set to 0 otherwise.
+
+- ``z`` -- Objective value to minimize. It is equal to the maximum over all step
+  `t` of the number of vertices such that `u_v^t==1`.
+
+MILP formulation:
+
+.. MATH::
+    :nowrap:
+
+    \begin{alignat}{2}
+    \intertext{Minimize:}
+    &z&\\
+    \intertext{Such that:}
+    x_v^t &\leq x_v^{t+1}& \forall v\in V,\ 0\leq t\leq n-2\\
+    y_v^t &\leq y_v^{t+1}& \forall v\in V,\ 0\leq t\leq n-2\\
+    y_v^t &\leq x_w^t& \forall v\in V,\ \forall w\in N^+(v),\ 0\leq t\leq n-1\\
+    \sum_{v \in V} y_v^0 &\leq 1&\\
+    \sum_{v \in V} y_v^{t+1} &\leq 1+\sum_{v \in V} y_v^t& 0\leq t\leq n-2\\
+    \sum_{v \in V} y_v^{N-1} &= n&\\
+    x_v^t-y_v^t&\leq u_v^t & \forall v \in V,\ 0\leq t\leq n-1\\
+    \sum_{v \in V} u_v^t &\leq z& 0\leq t\leq n-1\\
+    0 \leq x_v^t &\leq 1& \forall v\in V,\ 0\leq t\leq n-1\\
+    0 \leq u_v^t &\leq 1& \forall v\in V,\ 0\leq t\leq n-1\\
+    y_v^t &\in \{0,1\}& \forall v\in V,\ 0\leq t\leq n-1\\
+    0 \leq z &\leq n&
+    \end{alignat}
+
+The vertex separation of `G` is given by the value of `z`, and the order of
+vertex `v` in the optimal layout is given by the smallest `t` for which
+`y_v^t==1`.
+
+
+REFERENCES
+----------
+
+.. [Bod98] *A partial k-arboretum of graphs with bounded treewidth*, Hans
+  L. Bodlaender, Theoretical Computer Science 209(1-2):1-45, 1998.
+
+.. [Kin92] *The vertex separation number of a graph equals its path-width*,
+  Nancy G. Kinnersley, Information Processing Letters 42(6):345-350, 1992.
+
+.. [SP10] *Lightpath Reconfiguration in WDM networks*, Fernando Solano and
+  Michal Pioro, IEEE/OSA Journal of Optical Communication and Networking
+  2(12):1010-1021, 2010.
+
+
+AUTHORS
+-------
+
+- Nathann Cohen (2011-10): Initial version and exact exponential algorithm
+
+- David Coudert (2012-04): MILP formulation and tests functions
+
+
+
+METHODS
 -------
 """
 
@@ -154 +238 @@
     documentation).
 
     .. NOTE::
-        
+
         This method runs in exponential time but has no memory constraint.
 
 
@@ -217 +301 @@
 
     return min
 
-####################
-# Exact algorithms #
-####################
+################################
+# Exact exponential algorithms #
+################################
 
-def path_decomposition(G, verbose = False):
+def path_decomposition(G, algorithm = "exponential", verbose = False):
     r"""
     Returns the pathwidth of the given graph and the ordering of the vertices
     resulting in a corresponding path decomposition.
@@ -230 +314 @@
 
     - ``G`` -- a digraph
 
+    - ``algorithm`` -- (default: ``"exponential"``) Specify the algorithm to use
+      among
+
+      - ``exponential`` -- Use an exponential time and space algorithm. This
+        algorithm only works of graphs on less than 32 vertices.
+
+      - ``MILP`` -- Use a mixed integer linear programming formulation. This
+        algorithm has no size restriction but could take a very long time.
+
     - ``verbose`` (boolean) -- whether to display information on the
       computations.
 
@@ -240 +333 @@
 
     .. NOTE::
 
-        Because of its current implementation, this algorithm only works on
-        graphs on less than 32 vertices. This can be changed to 54 if necessary,
-        but 32 vertices already require 4GB of memory.
+        Because of its current implementation, this exponential algorithm only
+        works on graphs on less than 32 vertices. This can be changed to 54 if
+        necessary, but 32 vertices already require 4GB of memory.
 
     EXAMPLE:
 
-    The vertex separation of a circuit is equal to 2::
+    The vertex separation of a cycle is equal to 2::
 
         sage: from sage.graphs.graph_decompositions.vertex_separation import path_decomposition
         sage: g = graphs.CycleGraph(6)
-        sage: path_decomposition(g)
-        (2, [0, 1, 2, 3, 4, 5])
+        sage: pw, L = path_decomposition(g); pw
+        2
+        sage: pwm, Lm = path_decomposition(g, algorithm = "MILP"); pwm
+        2
 
     TEST:
 
@@ -267 +362 @@
     from sage.graphs.graph import Graph
     if not isinstance(G, Graph):
         raise ValueError("The parameter must be a Graph.")
-    
+
     from sage.graphs.digraph import DiGraph
-    return vertex_separation(DiGraph(G), verbose = verbose)
+    if algorithm == "exponential":
+        return vertex_separation(DiGraph(G), verbose = verbose)
+    else:
+        return vertex_separation_MILP(DiGraph(G), verbosity = (1 if verbose else 0))
+
 
 def vertex_separation(G, verbose = False):
     r"""
@@ -296 +395 @@
 
     EXAMPLE:
 
-    The vertex separation of a circuit is equal to 2::
+    The vertex separation of a circuit is equal to 1::
 
         sage: from sage.graphs.graph_decompositions.vertex_separation import vertex_separation
         sage: g = digraphs.Circuit(6)
         sage: vertex_separation(g)
-        (1, [0, 1, 2, 3, 4, 5])        
+        (1, [0, 1, 2, 3, 4, 5])
 
     TEST:
 
@@ -316 +415 @@
 
     Graphs with non-integer vertices::
 
+        sage: from sage.graphs.graph_decompositions.vertex_separation import vertex_separation
         sage: D=digraphs.DeBruijn(2,3)
         sage: vertex_separation(D)
         (2, ['000', '001', '100', '010', '101', '011', '110', '111'])
@@ -454 +554 @@
         return a
     else:
         return b
+
+
+#################################################################
+# Function for testing the validity of a linear vertex ordering #
+#################################################################
+
+def is_valid_ordering(G, L):
+    r"""
+    Test if the linear vertex ordering `L` is valid for (di)graph `G`.
+
+    A linear ordering `L` of the vertices of a (di)graph `G` is valid if all
+    vertices of `G` are in `L`, and if `L` contains no other vertex and no
+    duplicated vertices.
+
+    INPUT:
+
+    - ``G`` -- a Graph or a DiGraph.
+
+    - ``L`` -- an ordered list of the vertices of ``G``.
+
+
+    OUTPUT:
+
+    Returns `True` if `L` is a valid vertex ordering for `G`, and False
+    oterwise.
+
+
+    EXAMPLE:
+
+    Path decomposition of a cycle::
+
+        sage: from sage.graphs.graph_decompositions import vertex_separation
+        sage: G = graphs.CycleGraph(6)
+        sage: L = [u for u in G.vertices()]
+        sage: vertex_separation.is_valid_ordering(G, L)
+        True
+        sage: vertex_separation.is_valid_ordering(G, [1,2])
+        False
+
+    TEST:
+
+    Giving anything else than a Graph or a DiGraph::
+
+        sage: from sage.graphs.graph_decompositions import vertex_separation
+        sage: vertex_separation.is_valid_ordering(2, [])
+        Traceback (most recent call last):
+        ...
+        ValueError: The input parameter must be a Graph or a DiGraph.
+
+    Giving anything else than a list::
+
+        sage: from sage.graphs.graph_decompositions import vertex_separation
+        sage: G = graphs.CycleGraph(6)
+        sage: vertex_separation.is_valid_ordering(G, {})
+        Traceback (most recent call last):
+        ...
+        ValueError: The second parameter must be of type 'list'.
+    """
+    from sage.graphs.graph import Graph
+    from sage.graphs.digraph import DiGraph
+    if not isinstance(G, Graph) and not isinstance(G, DiGraph):
+        raise ValueError("The input parameter must be a Graph or a DiGraph.")
+    if not isinstance(L, list):
+        raise ValueError("The second parameter must be of type 'list'.")
+
+    return len( set( G.vertices() ).symmetric_difference( set(L) ) ) == 0
+
+
+####################################################################
+# Measurement functions of the widths of some graph decompositions #
+####################################################################
+
+def width_of_path_decomposition(G, L):
+    r"""
+    Returns the width of the path decomposition induced by the linear ordering
+    `L` of the vertices of `G`.
+
+    If `G` is an instance of Graph, this function returns the width `pw_L(G)` of
+    the path decomposition induced by the linear ordering `L` of the vertices of
+    `G`. If `G` is a DiGraph, it returns instead the width `vs_L(G)` of the
+    directed path decomposition induced by the linear ordering `L` of the
+    vertices of `G`, where
+
+    .. MATH::
+
+        vs_L(G) & =  \max_{0\leq i< |V|-1} | N^+(L[:i])\setminus L[:i] |\\
+        pw_L(G) & =  \max_{0\leq i< |V|-1} | N(L[:i])\setminus L[:i] |\\
+
+    INPUT:
+
+    - ``G`` -- a Graph or a DiGraph
+
+    - ``L`` -- a linear ordering of the vertices of ``G``
+
+    EXAMPLES:
+
+    Path decomposition of a cycle::
+
+        sage: from sage.graphs.graph_decompositions import vertex_separation
+        sage: G = graphs.CycleGraph(6)
+        sage: L = [u for u in G.vertices()]
+        sage: vertex_separation.width_of_path_decomposition(G, L)
+        2
+
+    Directed path decomposition of a circuit::
+
+        sage: from sage.graphs.graph_decompositions import vertex_separation
+        sage: G = digraphs.Circuit(6)
+        sage: L = [u for u in G.vertices()]
+        sage: vertex_separation.width_of_path_decomposition(G, L)
+        1
+
+    TESTS:
+
+    Path decomposition of a BalancedTree::
+
+        sage: from sage.graphs.graph_decompositions import vertex_separation
+        sage: G = graphs.BalancedTree(3,2)
+        sage: pw, L = vertex_separation.path_decomposition(G)
+        sage: pw == vertex_separation.width_of_path_decomposition(G, L)
+        True
+        sage: L.reverse()
+        sage: pw == vertex_separation.width_of_path_decomposition(G, L)
+        False
+
+    Directed path decomposition of a circuit::
+
+        sage: from sage.graphs.graph_decompositions import vertex_separation
+        sage: G = digraphs.Circuit(8)
+        sage: vs, L = vertex_separation.vertex_separation(G)
+        sage: vs == vertex_separation.width_of_path_decomposition(G, L)
+        True
+        sage: L = [0,4,6,3,1,5,2,7]
+        sage: vs == vertex_separation.width_of_path_decomposition(G, L)
+        False
+
+    Giving a wrong linear ordering::
+
+        sage: from sage.graphs.graph_decompositions import vertex_separation
+        sage: G = Graph()
+        sage: vertex_separation.width_of_path_decomposition(G, ['a','b'])
+        Traceback (most recent call last):
+        ...
+        ValueError: The input linear vertex ordering L is not valid for G.
+    """
+    if not is_valid_ordering(G, L):
+        raise ValueError("The input linear vertex ordering L is not valid for G.")
+
+    vsL = 0
+    S = []
+    neighbors_of_S_in_V_minus_S = []
+
+    for u in L:
+
+        # We remove u from the list of neighbors of S
+        if u in neighbors_of_S_in_V_minus_S:
+            neighbors_of_S_in_V_minus_S.remove(u)
+
+        # We add vertex u to the set S
+        S.append(u)
+
+        if G._directed:
+            Nu = G.neighbors_out(u)
+        else:
+            Nu = G.neighbors(u)
+
+        # We add the (out-)neighbors of u to the list of neighbors of S
+        for v in Nu:
+            if (not v in S) and (not v in neighbors_of_S_in_V_minus_S):
+                neighbors_of_S_in_V_minus_S.append(v)
+
+        # We update the cost of the vertex separation
+        vsL = max( vsL, len(neighbors_of_S_in_V_minus_S) )
+
+    return vsL
+
+
+##########################################
+# MILP formulation for vertex separation #
+##########################################
+
+def vertex_separation_MILP(G, integrality = False, solver = None, verbosity = 0):
+    r"""
+    Computes the vertex separation of `G` and the optimal ordering of its
+    vertices using an MILP formulation.
+
+    This function uses a mixed integer linear program (MILP) for determining an
+    optimal layout for the vertex separation of `G`. This MILP is an improved
+    version of the formulation proposed in [SP10]_. See the module's
+    documentation for more details on this MILP formulation.
+
+    INPUTS:
+
+    - ``G`` -- a DiGraph
+
+    - ``integrality`` -- (default: ``False``) Specify if all variables must be
+      integral of if some variables can be relaxed.
+
+    - ``solver`` -- (default: ``None``) Specify a Linear Program (LP) solver to
+      be used. If set to ``None``, the default one is used. For more information
+      on LP solvers and which default solver is used, see the method
+      :meth:`solve<sage.numerical.mip.MixedIntegerLinearProgram.solve>` of the
+      class
+      :class:`MixedIntegerLinearProgram<sage.numerical.mip.MixedIntegerLinearProgram>`.
+
+    - ``verbose`` -- integer (default: ``0``). Sets the level of verbosity. Set
+      to 0 by default, which means quiet.
+
+    OUTPUT:
+
+    A pair ``(cost, ordering)`` representing the optimal ordering of the
+    vertices and its cost.
+
+    EXAMPLE:
+
+    Vertex separation of a De Bruijn digraph::
+
+        sage: from sage.graphs.graph_decompositions import vertex_separation
+        sage: G = digraphs.DeBruijn(2,3)
+        sage: vs, L = vertex_separation.vertex_separation_MILP(G); vs
+        2
+        sage: vs == vertex_separation.width_of_path_decomposition(G, L)
+        True
+        sage: vse, Le = vertex_separation.vertex_separation(G); vse
+        2
+
+    The vertex separation of a circuit is 1::
+
+        sage: from sage.graphs.graph_decompositions import vertex_separation
+        sage: G = digraphs.Circuit(6)
+        sage: vs, L = vertex_separation.vertex_separation_MILP(G); vs
+        1
+
+    TESTS:
+
+    Comparison with exponential algorithm::
+
+        sage: from sage.graphs.graph_decompositions import vertex_separation
+        sage: for i in range(20):
+        ...       G = digraphs.RandomDirectedGNP(10,0.2)
+        ...       ve,le = vertex_separation.vertex_separation(G)
+        ...       vm,lm = vertex_separation.vertex_separation_MILP(G)
+        ...       if ve != vm:
+        ...          print "The solution is not optimal!"
+
+    Giving anything else than a DiGraph::
+
+        sage: from sage.graphs.graph_decompositions import vertex_separation
+        sage: vertex_separation.vertex_separation_MILP([])
+        Traceback (most recent call last):
+        ...
+        ValueError: The first input parameter must be a DiGraph.
+    """
+    from sage.graphs.digraph import DiGraph
+    if not isinstance(G, DiGraph):
+        raise ValueError("The first input parameter must be a DiGraph.")
+
+    from sage.numerical.mip import MixedIntegerLinearProgram, Sum, MIPSolverException
+    p = MixedIntegerLinearProgram( maximization = False, solver = solver )
+
+    # Declaration of variables.
+    x = p.new_variable( integer = integrality, dim = 2 )
+    u = p.new_variable( integer = integrality, dim = 2 )
+    y = p.new_variable( integer = True, dim = 2 )
+    z = p.new_variable( integer = True, dim = 1 )
+
+    N = G.num_verts()
+    V = G.vertices()
+
+    #  (2) x[v][t] <= x[v][t+1]   for all v in V, and for t:=0..N-2
+    #  (3) y[v][t] <= y[v][t+1]   for all v in V, and for t:=0..N-2
+    for v in V:
+        for t in xrange(N-1):
+            p.add_constraint( x[v][t] - x[v][t+1], max = 0 )
+            p.add_constraint( y[v][t] - y[v][t+1], max = 0 )
+
+    #  (4) y[v][t] <= x[w][t]  for all v in V, for all w in N^+(v), and for all t:=0..N-1
+    for v in V:
+        for w in G.neighbors_out(v):
+            for t in xrange(N):
+                p.add_constraint( y[v][t] - x[w][t], max = 0 )
+
+    #  (5) sum_{v in V} y[v][0] <= 1
+    p.add_constraint( Sum([ y[v][0] for v in V ]), max = 1 )
+
+    #  (6) sum_{v in V} y[v][t+1] <= sum_{v in V} y[v][t] + 1  for t:=0..N-2
+    for t in xrange(N-1):
+        p.add_constraint( Sum([ y[v][t+1] - y[v][t] for v in V ]), max = 1 )
+
+    #  (7) sum_{v in V} y[v][N-1] = N
+    p.add_constraint( Sum([ y[v][N-1] for v in V ]), min = N )
+    p.add_constraint( Sum([ y[v][N-1] for v in V ]), max = N )
+
+    #  (8) u[v][t] >= x[v][t]-y[v][t]    for all v in V, and for all t:=0..N-1
+    for v in V:
+        for t in xrange(N):
+            p.add_constraint( x[v][t] - y[v][t] - u[v][t], max = 0 )
+
+    #  (9) z >= sum_{v in V} u[v][t]   for all t:=0..N-1
+    for t in xrange(N):
+        p.add_constraint( Sum([ u[v][t] for v in V ]) - z['z'], max = 0 )
+
+    # (10) 0 <= x[v][t] and u[v][t] <= 1
+    # (11) y[v][t] in {0,1}
+    for v in V:
+        for t in xrange(N):
+            p.add_constraint( x[v][t], min = 0 )
+            p.add_constraint( x[v][t], max = 1 )
+            p.add_constraint( u[v][t], min = 0 )
+            p.add_constraint( u[v][t], max = 1 )
+            p.set_binary( y[v][t] )
+
+    # (12) 0 <= z <= |V|
+    p.add_constraint( z['z'], min = 0 )
+    p.add_constraint( z['z'], max = N )
+
+    #  (1) Minimize z
+    p.set_objective( z['z'] )
+
+    try:
+        obj = p.solve( log=verbosity )
+
+        taby = p.get_values( y )
+        tabz = p.get_values( z )
+        # since exactly one vertex is processed per step, we can reconstruct the sequence
+        seq = []
+        for t in xrange(N):
+            for v in V:
+                if (taby[v][t] > 0) and (not v in seq):
+                    seq.append(v)
+                    break
+        vs = int(round( tabz['z'] ))
+
+    except MIPSolverException:
+        if integrality:
+            raise ValueError("Unbounded or unexpected error")
+        else:
+            raise ValueError("Unbounded or unexpected error. Try with 'integrality = True'.")
+
+    del p
+    return vs, seq
+