# HG changeset patch
# User Nathann Cohen <nathann.cohen@gmail.com>
# Date 1334072918 7200
# Node ID de6bc52cf4a69cc72a8dd2aa4087d0ffcb07c1bb
# Parent cb71d94752d432671881ce3bfbcf78f9d8248379
MILP formulation and test functions for vertex separation  reviewer's patch
diff git a/sage/graphs/graph_decompositions/vertex_separation.pyx b/sage/graphs/graph_decompositions/vertex_separation.pyx
a

b


3  3  
4  4  This module implements several algorithms to compute the vertex separation of a 
5  5  digraph and the corresponding ordering of the vertices. It also implements tests 
6   functions for evaluation the width of a linear orderings. 
 6  functions for evaluation the width of a linear ordering. 
7  7  
8  8  Given an ordering 
9  9  `v_1,\cdots, v_n` of the vertices of `V(G)`, its *cost* is defined as: 
… 
… 

12  12  
13  13  c(v_1, ..., v_n) = \max_{1\leq i \leq n} c'(\{v_1, ..., v_i\}) 
14  14  
15   Where 
 15  Where 
16  16  
17  17  .. MATH:: 
18  18  
… 
… 

112  112  
113  113  Because of its current implementation, this algorithm only works on graphs 
114  114  on less than 32 vertices. This can be changed to 64 if necessary, but 32 
115   vertices already require 4GB of memory. 
 115  vertices already require 4GB of memory. Running it on 64 bits is not 
 116  expected to be doable by the computers of the next decade `:D` 
116  117  
117  118  **Lower bound on the vertex separation** 
118  119  
… 
… 

133  134  MILP formulation for the vertex separation 
134  135   
135  136  
136   We describe bellow a mixed integer linear program (MILP) for determining an 
137   optimal layout for the vertex separation of `G`. This MILP is an improved 
138   version of the formulation proposed in [SP10]_. 
 137  We describe below a mixed integer linear program (MILP) for determining an 
 138  optimal layout for the vertex separation of `G`, which is an improved version of 
 139  the formulation proposed in [SP10]_. It aims at building a sequence `S_t` of 
 140  sets such that an ordering `v_1, ..., v_n` of the vertices correspond to 
 141  `S_0=\{v_1\}, S_2=\{v_1,v_2\}, ..., S_{n1}=\{v_1,...,v_n\}`. 
139  142  
140   VARIABLES: 
 143  **Variables:** 
141  144  
142    ``x_v^t``  Variable set to 1 if either `v` has an inneighbor `u` such that 
143   `u\in S_{t'}=\{v_1, v_2,\cdots, v_{t'}\}` and `v\in N_G^+(S_{t'})\backslash 
144   S_{t'}` with `t'\leq t`, or `v\in S_t`. It is set to 0 otherwise. 
145  145  
146    ``y_v^t``  Variable set to 1 if `v\in S_t`, and 0 otherwise. The order of 
 146   `y_v^t`  Variable set to 1 if `v\in S_t`, and 0 otherwise. The order of 
147  147  `v` in the layout is the smallest `t` such that `y_v^t==1`. 
148  148  
149    ``u_v^t``  Variable set to 1 if `v\in N_G^+(S_t)\backslash S_t` has an 
150   inneighbor in `S_t`, and set to 0 otherwise. 
 149   `u_v^t`  Variable set to 1 if `v\not \in S_t` and `v` has an inneighbor in 
 150  `S_t`. It is set to 0 otherwise. 
151  151  
152    ``z``  Objective value to minimize. It is equal to the maximum over all step 
 152   `x_v^t`  is morally equal to "`y_v^t+u_v^t`". It is set to 1 if either 
 153  `v\in S_t` or if `v` has an inneighbor in `S_t`. It is set to 0 otherwise. 
 154  
 155   `z`  Objective value to minimize. It is equal to the maximum over all step 
153  156  `t` of the number of vertices such that `u_v^t==1`. 
154  157  
155   MILP formulation: 
 158  **MILP formulation:** 
156  159  
157  160  .. MATH:: 
158  161  :nowrap: 
… 
… 

179  182  vertex `v` in the optimal layout is given by the smallest `t` for which 
180  183  `y_v^t==1`. 
181  184  
182   
183  185  REFERENCES 
184  186   
185  187  
… 
… 

244  246  
245  247  EXAMPLE: 
246  248  
247   On a cycle:: 
 249  On a circuit:: 
248  250  
249  251  sage: from sage.graphs.graph_decompositions.vertex_separation import lower_bound 
250  252  sage: g = digraphs.Circuit(6) 
… 
… 

464  466  
465  467  sage_free(neighborhoods) 
466  468  
467   _sig_off 
 469  _sig_off 
468  470  
469  471  return k, order 
470  472  
… 
… 

577  579  
578  580  OUTPUT: 
579  581  
580   Returns `True` if `L` is a valid vertex ordering for `G`, and False 
 582  Returns ``True`` if `L` is a valid vertex ordering for `G`, and ``False`` 
581  583  oterwise. 
582  584  
583  585  
… 
… 

619  621  if not isinstance(L, list): 
620  622  raise ValueError("The second parameter must be of type 'list'.") 
621  623  
622   return len( set( G.vertices() ).symmetric_difference( set(L) ) ) == 0 
 624  return set(L) == set(G.vertices()) 
623  625  
624  626  
625  627  #################################################################### 
… 
… 

631  633  Returns the width of the path decomposition induced by the linear ordering 
632  634  `L` of the vertices of `G`. 
633  635  
634   If `G` is an instance of Graph, this function returns the width `pw_L(G)` of 
635   the path decomposition induced by the linear ordering `L` of the vertices of 
636   `G`. If `G` is a DiGraph, it returns instead the width `vs_L(G)` of the 
 636  If `G` is an instance of :mod:`Graph <sage.graphs.graph>`, this function 
 637  returns the width `pw_L(G)` of the path decomposition induced by the linear 
 638  ordering `L` of the vertices of `G`. If `G` is a :mod:`DiGraph 
 639  <sage.graphs.digraph>`, it returns instead the width `vs_L(G)` of the 
637  640  directed path decomposition induced by the linear ordering `L` of the 
638  641  vertices of `G`, where 
639  642  
… 
… 

703  706  raise ValueError("The input linear vertex ordering L is not valid for G.") 
704  707  
705  708  vsL = 0 
706   S = [] 
707   neighbors_of_S_in_V_minus_S = [] 
 709  S = set() 
 710  neighbors_of_S_in_V_minus_S = set() 
708  711  
709  712  for u in L: 
710  713  
711   # We remove u from the list of neighbors of S 
712   if u in neighbors_of_S_in_V_minus_S: 
713   neighbors_of_S_in_V_minus_S.remove(u) 
 714  # We remove u from the neighbors of S 
 715  neighbors_of_S_in_V_minus_S.discard(u) 
714  716  
715  717  # We add vertex u to the set S 
716   S.append(u) 
 718  S.add(u) 
717  719  
718  720  if G._directed: 
719  721  Nu = G.neighbors_out(u) 
720  722  else: 
721  723  Nu = G.neighbors(u) 
722  724  
723   # We add the (out)neighbors of u to the list of neighbors of S 
 725  # We add the (out)neighbors of u to the neighbors of S 
724  726  for v in Nu: 
725   if (not v in S) and (not v in neighbors_of_S_in_V_minus_S): 
726   neighbors_of_S_in_V_minus_S.append(v) 
 727  if (not v in S): 
 728  neighbors_of_S_in_V_minus_S.add(v) 
727  729  
728  730  # We update the cost of the vertex separation 
729  731  vsL = max( vsL, len(neighbors_of_S_in_V_minus_S) ) 
… 
… 

742  744  
743  745  This function uses a mixed integer linear program (MILP) for determining an 
744  746  optimal layout for the vertex separation of `G`. This MILP is an improved 
745   version of the formulation proposed in [SP10]_. See the module's 
746   documentation for more details on this MILP formulation. 
 747  version of the formulation proposed in [SP10]_. See the :mod:`module's 
 748  documentation <sage.graphs.graph_decompositions.vertex_separation>` for more 
 749  details on this MILP formulation. 
747  750  
748  751  INPUTS: 
749  752  
… 
… 

815  818  p = MixedIntegerLinearProgram( maximization = False, solver = solver ) 
816  819  
817  820  # Declaration of variables. 
818   x = p.new_variable( integer = integrality, dim = 2 ) 
819   u = p.new_variable( integer = integrality, dim = 2 ) 
820   y = p.new_variable( integer = True, dim = 2 ) 
 821  x = p.new_variable( binary = integrality, dim = 2 ) 
 822  u = p.new_variable( binary = integrality, dim = 2 ) 
 823  y = p.new_variable( binary = True, dim = 2 ) 
821  824  z = p.new_variable( integer = True, dim = 1 ) 
822  825  
823  826  N = G.num_verts() 
… 
… 

858  861  
859  862  # (10) 0 <= x[v][t] and u[v][t] <= 1 
860  863  # (11) y[v][t] in {0,1} 
861   for v in V: 
862   for t in xrange(N): 
863   p.add_constraint( x[v][t], min = 0 ) 
864   p.add_constraint( x[v][t], max = 1 ) 
865   p.add_constraint( u[v][t], min = 0 ) 
866   p.add_constraint( u[v][t], max = 1 ) 
867   p.set_binary( y[v][t] ) 
 864  if not integrality: 
 865  for v in V: 
 866  for t in xrange(N): 
 867  p.add_constraint( x[v][t], min = 0 ) 
 868  p.add_constraint( x[v][t], max = 1 ) 
 869  p.add_constraint( u[v][t], min = 0 ) 
 870  p.add_constraint( u[v][t], max = 1 ) 
 871  
 872  # (11) y[v][t] in {0,1} 
 873  p.set_binary( y ) 
868  874  
869  875  # (12) 0 <= z <= V 
870   p.add_constraint( z['z'], min = 0 ) 
871  876  p.add_constraint( z['z'], max = N ) 
872  877  
873  878  # (1) Minimize z 
… 
… 

875  880  
876  881  try: 
877  882  obj = p.solve( log=verbosity ) 
878   
879   taby = p.get_values( y ) 
880   tabz = p.get_values( z ) 
881   # since exactly one vertex is processed per step, we can reconstruct the sequence 
882   seq = [] 
883   for t in xrange(N): 
884   for v in V: 
885   if (taby[v][t] > 0) and (not v in seq): 
886   seq.append(v) 
887   break 
888   vs = int(round( tabz['z'] )) 
889   
890  883  except MIPSolverException: 
891  884  if integrality: 
892  885  raise ValueError("Unbounded or unexpected error") 
893  886  else: 
894  887  raise ValueError("Unbounded or unexpected error. Try with 'integrality = True'.") 
895  888  
896   del p 
 889  taby = p.get_values( y ) 
 890  tabz = p.get_values( z ) 
 891  # since exactly one vertex is processed per step, we can reconstruct the sequence 
 892  seq = [] 
 893  for t in xrange(N): 
 894  for v in V: 
 895  if (taby[v][t] > 0) and (not v in seq): 
 896  seq.append(v) 
 897  break 
 898  vs = int(round( tabz['z'] )) 
 899  
897  900  return vs, seq 
898   