Ticket #12716: trac_12716-combined.patch

File trac_12716-combined.patch, 22.6 KB (added by dcoudert, 8 years ago)
  • sage/graphs/graph_decompositions/vertex_separation.pyx

    # HG changeset patch
    # User dcoudert <david.coudert@inria.fr>
    # Date 1334084333 -7200
    # Node ID 7326fcc2254e51398b0d27375cb929ca4ad742a3
    # Parent  7d8ac0ee679de0e5eeba4966fe280f5f65295472
    Trac #12716 -- MILP for vetex separation
    * * *
    Trac #12716 -- MILP formulation and test functions for vertex separation
    * * *
    MILP formulation and test functions for vertex separation -- reviewer's patch
    * * *
    Trac #12716 -- modifications
    
    diff --git a/sage/graphs/graph_decompositions/vertex_separation.pyx b/sage/graphs/graph_decompositions/vertex_separation.pyx
    a b  
    11r"""
    22Vertex separation
    33
    4 This module implements several algorithms to compute the vertex
    5 separation of a digraph and the corresponding ordering of the
    6 vertices. Given an ordering `v_1, ..., v_n` of the vertices of `V(G)`,
    7 its *cost* is defined as:
     4This module implements several algorithms to compute the vertex separation of a
     5digraph and the corresponding ordering of the vertices. It also implements tests
     6functions for evaluation the width of a linear ordering.
     7
     8Given an ordering
     9`v_1,\cdots, v_n` of the vertices of `V(G)`, its *cost* is defined as:
    810
    911.. MATH::
    1012
    1113    c(v_1, ..., v_n) = \max_{1\leq i \leq n} c'(\{v_1, ..., v_i\})
    1214
    13 Where 
     15Where
    1416
    1517.. MATH::
    1618
    1719    c'(S) = |N^+_G(S)\backslash S|
    1820
    19 The *vertex separation* of a digraph `G` is equal to the minimum cost
    20 of an ordering of its vertices.
     21The *vertex separation* of a digraph `G` is equal to the minimum cost of an
     22ordering of its vertices.
    2123
    2224**Vertex separation and pathwidth**
    2325
     
    2527`G` a digraph `D` with the same vertex set, and in which each edge `uv` of `G`
    2628is replaced by two edges `uv` and `vu` in `D`. The vertex separation of `D` is
    2729equal to the pathwidth of `G`, and the corresponding ordering of the vertices of
    28 `D` encodes an optimal path-decomposition of `G`.
    29 
     30`D`, also called a *layout*, encodes an optimal path-decomposition of `G`.
    3031This is a result of Kinnersley [Kin92]_ and Bodlaender [Bod98]_.
    3132
    32 **References**
    3333
    34 .. [Kin92] *The vertex separation number of a graph equals its path-width*,
    35   Nancy G. Kinnersley,
    36   Information Processing Letters,
    37   Volume 42, Issue 6, Pages 345-350,
    38   24 July 1992
     34**This module contains the following methods**
    3935
    40 .. [Bod98] *A partial k-arboretum of graphs with bounded treewidth*,
    41   Hans L. Bodlaender,
    42   Theoretical Computer Science,
    43   Volume 209, Issues 1-2, Pages 1-45,
    44   6 December 1998
     36.. csv-table::
     37    :class: contentstable
     38    :widths: 30, 70
     39    :delim: |
    4540
    46 **Authors**
     41    :meth:`path_decomposition` | Returns the pathwidth of the given graph and the ordering of the vertices resulting in a corresponding path decomposition
     42    :meth:`vertex_separation` | Returns an optimal ordering of the vertices and its cost for vertex-separation
     43    :meth:`vertex_separation_MILP` | Computes the vertex separation of `G` and the optimal ordering of its vertices using an MILP formulation
     44    :meth:`lower_bound` | Returns a lower bound on the vertex separation of `G`
     45    :meth:`is_valid_ordering` | Test if the linear vertex ordering `L` is valid for (di)graph `G`
     46    :meth:`width_of_path_decomposition` | Returns the width of the path decomposition induced by the linear ordering `L` of the vertices of `G`
    4747
    48     - Nathann Cohen
    4948
    5049Exponential algorithm for vertex separation
    5150-------------------------------------------
     
    101100
    102101    Because of its current implementation, this algorithm only works on graphs
    103102    on less than 32 vertices. This can be changed to 64 if necessary, but 32
    104     vertices already require 4GB of memory.
     103    vertices already require 4GB of memory. Running it on 64 bits is not
     104    expected to be doable by the computers of the next decade `:-D`
    105105
    106106**Lower bound on the vertex separation**
    107107
     
    114114
    115115    \max c'(\{v_1\}),c'(\{v_1,v_2\}),...,c'(\{v_1,...,v_n\})\geq\max c'_1,...,c'_n
    116116
    117 where `c_i` is the minimum cost of a set `S` on `i`
    118 vertices. Evaluating the `c_i` can take time (and in particular more
    119 than the previous exact algorithm), but it does not need much memory
    120 to run.
     117where `c_i` is the minimum cost of a set `S` on `i` vertices. Evaluating the
     118`c_i` can take time (and in particular more than the previous exact algorithm),
     119but it does not need much memory to run.
    121120
    122 Methods
     121
     122MILP formulation for the vertex separation
     123------------------------------------------
     124
     125We describe below a mixed integer linear program (MILP) for determining an
     126optimal layout for the vertex separation of `G`, which is an improved version of
     127the formulation proposed in [SP10]_. It aims at building a sequence `S_t` of
     128sets such that an ordering `v_1, ..., v_n` of the vertices correspond to
     129`S_0=\{v_1\}, S_2=\{v_1,v_2\}, ..., S_{n-1}=\{v_1,...,v_n\}`.
     130
     131**Variables:**
     132
     133
     134- `y_v^t` -- Variable set to 1 if `v\in S_t`, and 0 otherwise. The order of
     135  `v` in the layout is the smallest `t` such that `y_v^t==1`.
     136
     137- `u_v^t` -- Variable set to 1 if `v\not \in S_t` and `v` has an in-neighbor in
     138  `S_t`. It is set to 0 otherwise.
     139
     140- `x_v^t` -- Variable set to 1 if either `v\in S_t` or if `v` has an in-neighbor
     141  in `S_t`. It is set to 0 otherwise.
     142
     143- `z` -- Objective value to minimize. It is equal to the maximum over all step
     144  `t` of the number of vertices such that `u_v^t==1`.
     145
     146**MILP formulation:**
     147
     148.. MATH::
     149    :nowrap:
     150
     151    \begin{alignat}{2}
     152    \intertext{Minimize:}
     153    &z&\\
     154    \intertext{Such that:}
     155    x_v^t &\leq x_v^{t+1}& \forall v\in V,\ 0\leq t\leq n-2\\
     156    y_v^t &\leq y_v^{t+1}& \forall v\in V,\ 0\leq t\leq n-2\\
     157    y_v^t &\leq x_w^t& \forall v\in V,\ \forall w\in N^+(v),\ 0\leq t\leq n-1\\
     158    \sum_{v \in V} y_v^{t} &= t+1& 0\leq t\leq n-1\\
     159    x_v^t-y_v^t&\leq u_v^t & \forall v \in V,\ 0\leq t\leq n-1\\
     160    \sum_{v \in V} u_v^t &\leq z& 0\leq t\leq n-1\\
     161    0 \leq x_v^t &\leq 1& \forall v\in V,\ 0\leq t\leq n-1\\
     162    0 \leq u_v^t &\leq 1& \forall v\in V,\ 0\leq t\leq n-1\\
     163    y_v^t &\in \{0,1\}& \forall v\in V,\ 0\leq t\leq n-1\\
     164    0 \leq z &\leq n&
     165    \end{alignat}
     166
     167The vertex separation of `G` is given by the value of `z`, and the order of
     168vertex `v` in the optimal layout is given by the smallest `t` for which
     169`y_v^t==1`.
     170
     171REFERENCES
     172----------
     173
     174.. [Bod98] *A partial k-arboretum of graphs with bounded treewidth*, Hans
     175  L. Bodlaender, Theoretical Computer Science 209(1-2):1-45, 1998.
     176
     177.. [Kin92] *The vertex separation number of a graph equals its path-width*,
     178  Nancy G. Kinnersley, Information Processing Letters 42(6):345-350, 1992.
     179
     180.. [SP10] *Lightpath Reconfiguration in WDM networks*, Fernando Solano and
     181  Michal Pioro, IEEE/OSA Journal of Optical Communication and Networking
     182  2(12):1010-1021, 2010.
     183
     184
     185AUTHORS
     186-------
     187
     188- Nathann Cohen (2011-10): Initial version and exact exponential algorithm
     189
     190- David Coudert (2012-04): MILP formulation and tests functions
     191
     192
     193
     194METHODS
    123195-------
    124196"""
    125197
     
    154226    documentation).
    155227
    156228    .. NOTE::
    157        
     229
    158230        This method runs in exponential time but has no memory constraint.
    159231
    160232
    161233    EXAMPLE:
    162234
    163     On a cycle::
     235    On a circuit::
    164236
    165237        sage: from sage.graphs.graph_decompositions.vertex_separation import lower_bound
    166238        sage: g = digraphs.Circuit(6)
     
    217289
    218290    return min
    219291
    220 ####################
    221 # Exact algorithms #
    222 ####################
     292################################
     293# Exact exponential algorithms #
     294################################
    223295
    224 def path_decomposition(G, verbose = False):
     296def path_decomposition(G, algorithm = "exponential", verbose = False):
    225297    r"""
    226298    Returns the pathwidth of the given graph and the ordering of the vertices
    227299    resulting in a corresponding path decomposition.
     
    230302
    231303    - ``G`` -- a digraph
    232304
     305    - ``algorithm`` -- (default: ``"exponential"``) Specify the algorithm to use
     306      among
     307
     308      - ``exponential`` -- Use an exponential time and space algorithm. This
     309        algorithm only works of graphs on less than 32 vertices.
     310
     311      - ``MILP`` -- Use a mixed integer linear programming formulation. This
     312        algorithm has no size restriction but could take a very long time.
     313
    233314    - ``verbose`` (boolean) -- whether to display information on the
    234315      computations.
    235316
     
    240321
    241322    .. NOTE::
    242323
    243         Because of its current implementation, this algorithm only works on
    244         graphs on less than 32 vertices. This can be changed to 54 if necessary,
    245         but 32 vertices already require 4GB of memory.
     324        Because of its current implementation, this exponential algorithm only
     325        works on graphs on less than 32 vertices. This can be changed to 54 if
     326        necessary, but 32 vertices already require 4GB of memory.
    246327
    247328    EXAMPLE:
    248329
    249     The vertex separation of a circuit is equal to 2::
     330    The vertex separation of a cycle is equal to 2::
    250331
    251332        sage: from sage.graphs.graph_decompositions.vertex_separation import path_decomposition
    252333        sage: g = graphs.CycleGraph(6)
    253         sage: path_decomposition(g)
    254         (2, [0, 1, 2, 3, 4, 5])
     334        sage: pw, L = path_decomposition(g); pw
     335        2
     336        sage: pwm, Lm = path_decomposition(g, algorithm = "MILP"); pwm
     337        2
    255338
    256339    TEST:
    257340
     
    267350    from sage.graphs.graph import Graph
    268351    if not isinstance(G, Graph):
    269352        raise ValueError("The parameter must be a Graph.")
    270    
     353
    271354    from sage.graphs.digraph import DiGraph
    272     return vertex_separation(DiGraph(G), verbose = verbose)
     355    if algorithm == "exponential":
     356        return vertex_separation(DiGraph(G), verbose = verbose)
     357    else:
     358        return vertex_separation_MILP(DiGraph(G), verbosity = (1 if verbose else 0))
     359
    273360
    274361def vertex_separation(G, verbose = False):
    275362    r"""
     
    296383
    297384    EXAMPLE:
    298385
    299     The vertex separation of a circuit is equal to 2::
     386    The vertex separation of a circuit is equal to 1::
    300387
    301388        sage: from sage.graphs.graph_decompositions.vertex_separation import vertex_separation
    302389        sage: g = digraphs.Circuit(6)
    303390        sage: vertex_separation(g)
    304         (1, [0, 1, 2, 3, 4, 5])       
     391        (1, [0, 1, 2, 3, 4, 5])
    305392
    306393    TEST:
    307394
     
    316403
    317404    Graphs with non-integer vertices::
    318405
     406        sage: from sage.graphs.graph_decompositions.vertex_separation import vertex_separation
    319407        sage: D=digraphs.DeBruijn(2,3)
    320408        sage: vertex_separation(D)
    321409        (2, ['000', '001', '100', '010', '101', '011', '110', '111'])
     
    364452
    365453    sage_free(neighborhoods)
    366454
    367     _sig_off 
     455    _sig_off
    368456
    369457    return k, order
    370458
     
    454542        return a
    455543    else:
    456544        return b
     545
     546
     547#################################################################
     548# Function for testing the validity of a linear vertex ordering #
     549#################################################################
     550
     551def is_valid_ordering(G, L):
     552    r"""
     553    Test if the linear vertex ordering `L` is valid for (di)graph `G`.
     554
     555    A linear ordering `L` of the vertices of a (di)graph `G` is valid if all
     556    vertices of `G` are in `L`, and if `L` contains no other vertex and no
     557    duplicated vertices.
     558
     559    INPUT:
     560
     561    - ``G`` -- a Graph or a DiGraph.
     562
     563    - ``L`` -- an ordered list of the vertices of ``G``.
     564
     565
     566    OUTPUT:
     567
     568    Returns ``True`` if `L` is a valid vertex ordering for `G`, and ``False``
     569    oterwise.
     570
     571
     572    EXAMPLE:
     573
     574    Path decomposition of a cycle::
     575
     576        sage: from sage.graphs.graph_decompositions import vertex_separation
     577        sage: G = graphs.CycleGraph(6)
     578        sage: L = [u for u in G.vertices()]
     579        sage: vertex_separation.is_valid_ordering(G, L)
     580        True
     581        sage: vertex_separation.is_valid_ordering(G, [1,2])
     582        False
     583
     584    TEST:
     585
     586    Giving anything else than a Graph or a DiGraph::
     587
     588        sage: from sage.graphs.graph_decompositions import vertex_separation
     589        sage: vertex_separation.is_valid_ordering(2, [])
     590        Traceback (most recent call last):
     591        ...
     592        ValueError: The input parameter must be a Graph or a DiGraph.
     593
     594    Giving anything else than a list::
     595
     596        sage: from sage.graphs.graph_decompositions import vertex_separation
     597        sage: G = graphs.CycleGraph(6)
     598        sage: vertex_separation.is_valid_ordering(G, {})
     599        Traceback (most recent call last):
     600        ...
     601        ValueError: The second parameter must be of type 'list'.
     602    """
     603    from sage.graphs.graph import Graph
     604    from sage.graphs.digraph import DiGraph
     605    if not isinstance(G, Graph) and not isinstance(G, DiGraph):
     606        raise ValueError("The input parameter must be a Graph or a DiGraph.")
     607    if not isinstance(L, list):
     608        raise ValueError("The second parameter must be of type 'list'.")
     609
     610    return set(L) == set(G.vertices())
     611
     612
     613####################################################################
     614# Measurement functions of the widths of some graph decompositions #
     615####################################################################
     616
     617def width_of_path_decomposition(G, L):
     618    r"""
     619    Returns the width of the path decomposition induced by the linear ordering
     620    `L` of the vertices of `G`.
     621
     622    If `G` is an instance of :mod:`Graph <sage.graphs.graph>`, this function
     623    returns the width `pw_L(G)` of the path decomposition induced by the linear
     624    ordering `L` of the vertices of `G`. If `G` is a :mod:`DiGraph
     625    <sage.graphs.digraph>`, it returns instead the width `vs_L(G)` of the
     626    directed path decomposition induced by the linear ordering `L` of the
     627    vertices of `G`, where
     628
     629    .. MATH::
     630
     631        vs_L(G) & =  \max_{0\leq i< |V|-1} | N^+(L[:i])\setminus L[:i] |\\
     632        pw_L(G) & =  \max_{0\leq i< |V|-1} | N(L[:i])\setminus L[:i] |\\
     633
     634    INPUT:
     635
     636    - ``G`` -- a Graph or a DiGraph
     637
     638    - ``L`` -- a linear ordering of the vertices of ``G``
     639
     640    EXAMPLES:
     641
     642    Path decomposition of a cycle::
     643
     644        sage: from sage.graphs.graph_decompositions import vertex_separation
     645        sage: G = graphs.CycleGraph(6)
     646        sage: L = [u for u in G.vertices()]
     647        sage: vertex_separation.width_of_path_decomposition(G, L)
     648        2
     649
     650    Directed path decomposition of a circuit::
     651
     652        sage: from sage.graphs.graph_decompositions import vertex_separation
     653        sage: G = digraphs.Circuit(6)
     654        sage: L = [u for u in G.vertices()]
     655        sage: vertex_separation.width_of_path_decomposition(G, L)
     656        1
     657
     658    TESTS:
     659
     660    Path decomposition of a BalancedTree::
     661
     662        sage: from sage.graphs.graph_decompositions import vertex_separation
     663        sage: G = graphs.BalancedTree(3,2)
     664        sage: pw, L = vertex_separation.path_decomposition(G)
     665        sage: pw == vertex_separation.width_of_path_decomposition(G, L)
     666        True
     667        sage: L.reverse()
     668        sage: pw == vertex_separation.width_of_path_decomposition(G, L)
     669        False
     670
     671    Directed path decomposition of a circuit::
     672
     673        sage: from sage.graphs.graph_decompositions import vertex_separation
     674        sage: G = digraphs.Circuit(8)
     675        sage: vs, L = vertex_separation.vertex_separation(G)
     676        sage: vs == vertex_separation.width_of_path_decomposition(G, L)
     677        True
     678        sage: L = [0,4,6,3,1,5,2,7]
     679        sage: vs == vertex_separation.width_of_path_decomposition(G, L)
     680        False
     681
     682    Giving a wrong linear ordering::
     683
     684        sage: from sage.graphs.graph_decompositions import vertex_separation
     685        sage: G = Graph()
     686        sage: vertex_separation.width_of_path_decomposition(G, ['a','b'])
     687        Traceback (most recent call last):
     688        ...
     689        ValueError: The input linear vertex ordering L is not valid for G.
     690    """
     691    if not is_valid_ordering(G, L):
     692        raise ValueError("The input linear vertex ordering L is not valid for G.")
     693
     694    vsL = 0
     695    S = set()
     696    neighbors_of_S_in_V_minus_S = set()
     697
     698    for u in L:
     699
     700        # We remove u from the neighbors of S
     701        neighbors_of_S_in_V_minus_S.discard(u)
     702
     703        # We add vertex u to the set S
     704        S.add(u)
     705
     706        if G._directed:
     707            Nu = G.neighbors_out(u)
     708        else:
     709            Nu = G.neighbors(u)
     710
     711        # We add the (out-)neighbors of u to the neighbors of S
     712        for v in Nu:
     713            if (not v in S):
     714                neighbors_of_S_in_V_minus_S.add(v)
     715
     716        # We update the cost of the vertex separation
     717        vsL = max( vsL, len(neighbors_of_S_in_V_minus_S) )
     718
     719    return vsL
     720
     721
     722##########################################
     723# MILP formulation for vertex separation #
     724##########################################
     725
     726def vertex_separation_MILP(G, integrality = False, solver = None, verbosity = 0):
     727    r"""
     728    Computes the vertex separation of `G` and the optimal ordering of its
     729    vertices using an MILP formulation.
     730
     731    This function uses a mixed integer linear program (MILP) for determining an
     732    optimal layout for the vertex separation of `G`. This MILP is an improved
     733    version of the formulation proposed in [SP10]_. See the :mod:`module's
     734    documentation <sage.graphs.graph_decompositions.vertex_separation>` for more
     735    details on this MILP formulation.
     736
     737    INPUTS:
     738
     739    - ``G`` -- a DiGraph
     740
     741    - ``integrality`` -- (default: ``False``) Specify if variables `x_v^t` and
     742      `u_v^t` must be integral or if they can be relaxed. This has no impact on
     743      the validity of the solution, but it is sometimes faster to solve the
     744      problem using binary variables only.
     745
     746    - ``solver`` -- (default: ``None``) Specify a Linear Program (LP) solver to
     747      be used. If set to ``None``, the default one is used. For more information
     748      on LP solvers and which default solver is used, see the method
     749      :meth:`solve<sage.numerical.mip.MixedIntegerLinearProgram.solve>` of the
     750      class
     751      :class:`MixedIntegerLinearProgram<sage.numerical.mip.MixedIntegerLinearProgram>`.
     752
     753    - ``verbose`` -- integer (default: ``0``). Sets the level of verbosity. Set
     754      to 0 by default, which means quiet.
     755
     756    OUTPUT:
     757
     758    A pair ``(cost, ordering)`` representing the optimal ordering of the
     759    vertices and its cost.
     760
     761    EXAMPLE:
     762
     763    Vertex separation of a De Bruijn digraph::
     764
     765        sage: from sage.graphs.graph_decompositions import vertex_separation
     766        sage: G = digraphs.DeBruijn(2,3)
     767        sage: vs, L = vertex_separation.vertex_separation_MILP(G); vs
     768        2
     769        sage: vs == vertex_separation.width_of_path_decomposition(G, L)
     770        True
     771        sage: vse, Le = vertex_separation.vertex_separation(G); vse
     772        2
     773
     774    The vertex separation of a circuit is 1::
     775
     776        sage: from sage.graphs.graph_decompositions import vertex_separation
     777        sage: G = digraphs.Circuit(6)
     778        sage: vs, L = vertex_separation.vertex_separation_MILP(G); vs
     779        1
     780
     781    TESTS:
     782
     783    Comparison with exponential algorithm::
     784
     785        sage: from sage.graphs.graph_decompositions import vertex_separation
     786        sage: for i in range(10):
     787        ...       G = digraphs.RandomDirectedGNP(10, 0.2)
     788        ...       ve, le = vertex_separation.vertex_separation(G)
     789        ...       vm, lm = vertex_separation.vertex_separation_MILP(G)
     790        ...       if ve != vm:
     791        ...          print "The solution is not optimal!"
     792
     793    Comparison with Different values of the integrality parameter::
     794
     795        sage: from sage.graphs.graph_decompositions import vertex_separation
     796        sage: for i in range(10):
     797        ...       G = digraphs.RandomDirectedGNP(10, 0.2)
     798        ...       va, la = vertex_separation.vertex_separation_MILP(G, integrality = False)
     799        ...       vb, lb = vertex_separation.vertex_separation_MILP(G, integrality = True)
     800        ...       if va != vb:
     801        ...          print "The integrality parameter change the result!"
     802
     803    Giving anything else than a DiGraph::
     804
     805        sage: from sage.graphs.graph_decompositions import vertex_separation
     806        sage: vertex_separation.vertex_separation_MILP([])
     807        Traceback (most recent call last):
     808        ...
     809        ValueError: The first input parameter must be a DiGraph.
     810    """
     811    from sage.graphs.digraph import DiGraph
     812    if not isinstance(G, DiGraph):
     813        raise ValueError("The first input parameter must be a DiGraph.")
     814
     815    from sage.numerical.mip import MixedIntegerLinearProgram, Sum, MIPSolverException
     816    p = MixedIntegerLinearProgram( maximization = False, solver = solver )
     817
     818    # Declaration of variables.
     819    x = p.new_variable( binary = integrality, dim = 2 )
     820    u = p.new_variable( binary = integrality, dim = 2 )
     821    y = p.new_variable( binary = True, dim = 2 )
     822    z = p.new_variable( integer = True, dim = 1 )
     823
     824    N = G.num_verts()
     825    V = G.vertices()
     826
     827    # (2) x[v][t] <= x[v][t+1]   for all v in V, and for t:=0..N-2
     828    # (3) y[v][t] <= y[v][t+1]   for all v in V, and for t:=0..N-2
     829    for v in V:
     830        for t in xrange(N-1):
     831            p.add_constraint( x[v][t] - x[v][t+1] <= 0 )
     832            p.add_constraint( y[v][t] - y[v][t+1] <= 0 )
     833
     834    # (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
     835    for v in V:
     836        for w in G.neighbors_out(v):
     837            for t in xrange(N):
     838                p.add_constraint( y[v][t] - x[w][t] <= 0 )
     839
     840    # (5) sum_{v in V} y[v][t] == t+1 for t:=0..N-1
     841    for t in xrange(N):
     842        p.add_constraint( Sum([ y[v][t] for v in V ]) == t+1 )
     843
     844    # (6) u[v][t] >= x[v][t]-y[v][t]    for all v in V, and for all t:=0..N-1
     845    for v in V:
     846        for t in xrange(N):
     847            p.add_constraint( x[v][t] - y[v][t] - u[v][t] <= 0 )
     848
     849    # (7) z >= sum_{v in V} u[v][t]   for all t:=0..N-1
     850    for t in xrange(N):
     851        p.add_constraint( Sum([ u[v][t] for v in V ]) - z['z'] <= 0 )
     852
     853    # (8)(9) 0 <= x[v][t] and u[v][t] <= 1
     854    if not integrality:
     855        for v in V:
     856            for t in xrange(N):
     857                p.add_constraint( 0 <= x[v][t] <= 1 )
     858                p.add_constraint( 0 <= u[v][t] <= 1 )
     859
     860    # (10) y[v][t] in {0,1}
     861    p.set_binary( y )
     862
     863    # (11) 0 <= z <= |V|
     864    p.add_constraint( z['z'] <= N )
     865
     866    #  (1) Minimize z
     867    p.set_objective( z['z'] )
     868
     869    try:
     870        obj = p.solve( log=verbosity )
     871    except MIPSolverException:
     872        if integrality:
     873            raise ValueError("Unbounded or unexpected error")
     874        else:
     875            raise ValueError("Unbounded or unexpected error. Try with 'integrality = True'.")
     876
     877    taby = p.get_values( y )
     878    tabz = p.get_values( z )
     879    # since exactly one vertex is processed per step, we can reconstruct the sequence
     880    seq = []
     881    for t in xrange(N):
     882        for v in V:
     883            if (taby[v][t] > 0) and (not v in seq):
     884                seq.append(v)
     885                break
     886    vs = int(round( tabz['z'] ))
     887
     888    return vs, seq