Generating Linear Extensions of a Partial Order

Question

Given a set $S$ and a partial order $\prec$ over $S$, I'm looking for a way to "efficiently" generate a list of linear extensions of $\prec$. Suppose the partial order is given by a List of pairs $\{x,y\}$ such that $x,y\in S$. For example, if $S = \{a,b,c\}$, then $\{\{a,c\},\{b,c\}\}$ defines a partial order where $a \prec c$ and $b \prec c$.

Essentially, given a list $S$ and a list of pairs $P$, I want to generate a list of permutations which respect the orders of the pairs in $P$.

linearExtensions[{a, b, c, d, e}, {{a, b}, {c, b}, {c, d}, {e, d}}]
(* {{a, c, b, e, d}, {a, c, e, b, d}, {a, c, e, d, b}, {a, e, c, b, d},
    {a, e, c, d, b}, {c, a, b, e, d}, {c, a, e, b, d}, {c, a, e, d, b},
    {c, e, a, b, d}, {c, e, a, d, b}, {c, e, d, a, b}, {e, a, c, b, d},
    {e, a, c, d, b}, {e, c, a, b, d}, {e, c, a, d, b}, {e, c, d, a, b}} *)

Since $\{a,b\}$ is in $P$, the permutation $\{b, c, a, d, e\}$ is not a linear extension because $b$ comes before $a$.

I've written two different functions which do the job, both using pattern matching, but I've been using Mathematica for less than a week and am still getting into the functional mindset. I'm interested to see how seasoned Mathematica users would tackle this problem.

My first approach was to use NestWhile (edit: this is ugly! I did not know about the Fold[] function when I wrote this):

linearExtensions[set_, po_] := 
    Module[{patterns},
        patterns = {___, #[[1]], ___, #[[2]], ___}& /@ po;
        First@NestWhile[{Cases[#[[1]], First[#[[2]]]], Rest[#[[2]]]}&,
            {Permutations[set], patterns}, Length[#[[2]]] > 0&]
    ]

My second approach, which turned out to be significantly slower, was to expand a set of rules from the partial order and use Select on the list of permutations.

linearExtensions[set_, po_] :=
    Module[{poQ},
        poQ[rule_] := And @@ (MatchQ[rule, {___, #[[1]], ___, #[[2]], ___}]& /@ po);
        Select[Permutations[set], poQ[#]&]
    ]

Note that on a list of length $n$ there can be $\Omega(n!)$ linear extensions, so by "efficient," I don't mean polynomial-time.

Answer 1

Could use integer linear programming. In Mathematica this can be done with Reduce[]. One way (probably not the best) to set this up is shown below. It uses an array of 0-1 variables, where a 1 in position (j,k) will indicate that the kth element of the input variables goes in position j of a particular ordering.

Caveat: I may have mixed up rows and columns.

consistentOrders[elems_, pairorders_] := Module[
  {n = Length[elems], vars, x, fvars, c1, c2, c3, c4, constraints, 
   ineqs, solns},
  ineqs[{a1_, a2_}, n_, v_] := 
   Table[Total[Take[v[[All, a1]], j]] >= 
     Total[Take[v[[All, a2]], j]], {j, 1, n - 1}];
  vars = Array[x, {n, n}];
  fvars = Flatten[vars];
  c1 = Map[0 <= # <= 1 &, fvars];
  c2 = Thread[Total[vars] == 1];
  c3 = Thread[Total[Transpose@vars] == 1];
  c4 = Flatten[
    Map[ineqs[#, n, vars] &, pairorders /. Thread[elems -> Range[n]]]];
  constraints = Join[c1, c2, c3, c4];
  solns = Reduce[constraints, fvars, Integers];
  solns = 
   solns /. {(_ == 0) :> Sequence[], aa_ == 1 :> aa, And -> List, 
     Or -> List};
  Sort[solns /. x[i_, j_] :> elems[[j]]]
  ]

--- edit ---

The first three constraint sets are fairly standard for this type of 0-1 programming. The fourth constraint subset arises as follows. The idea is that if the jth list element must precede the kth, then the 1 in column j must occur in an earlier row than the one in column k. So for every 1<=m<=n-1 (n=dimension) the sum of the first m entries in col j >= corresponding sum in col k.

--- end edit --- Example:

consistentOrders[{a, b, c, d, 
  e}, {{a, b}, {c, b}, {c, d}, {e, d}}]

(* Out[83]= {{a, c, b, e, d}, {a, c, e, b, d}, {a, c, e, d, b}, {a, e, c,
   b, d}, {a, e, c, d, b}, {c, a, b, e, d}, {c, a, e, b, d}, {c, a, e,
   d, b}, {c, e, a, b, d}, {c, e, a, d, b}, {c, e, d, a, b}, {e, a, c,
   b, d}, {e, a, c, d, b}, {e, c, a, b, d}, {e, c, a, d, b}, {e, c, d,
   a, b}} *)

Bigger example:

vars = {a, b, c, d, e, f, g, h, i, j, k, l};
porderlist = {{a, c}, {b, c}, {f, g}, {g, e}, {d, a}, {h, i}, {i, 
    d}, {g, h}, {g, i}, {h, j}, {h, k}, {k, j}, {k, c}};

Timing[ss = consistentOrders[vars, porderlist];]

(* Out[81]= {60.03, Null}

In[82]:= Length[ss]

Out[82]= 12840 *)

Answer 2

A more concise version of the first approach, using Fold[] instead of NestWhile[]:

linearExtensions[set_, po_] :=
    Fold[Cases, Permutations[set], {___, #1, ___, #2, ___} & @@@ po]

Answer 3

My first modest attempt:

linearExtensions[set_List, po_?MatrixQ] :=
    Select[Permutations[set], Complement[po, Subsets[#, {Last[Dimensions[po]]}]] === {} &]

My second modest attempt:

linearExtensions[set_List, po_?MatrixQ] :=
  Select[Permutations[set], And @@ Map[Function[p, LongestCommonSequence[#, p] === p], po] &]

Answer 4

I like @Daniel's approach using 0/1 integer programming. Instead of Reduce, however, one can achieve more optimal performance using SatisfiabilityInstances. In the code below, I've written the main function (AllTotalOrderings) for determining all total orderings of a partially ordered graph. It is then easy to encode the desired linearExtensions function as a simple call to AllTotalOrderings, for a suitably constructed graph. To further increase performance I've made use of CanonicalGraph (so as to introduce the constraints to the SAT solver in the most efficient order) and TransitiveReductionGraph (to minimize the number of constraints given to the SAT solver.) Note that I'm using the undocumented Method->"BDD" to call SatisfiabilityInstances, and BooleanConvert to preprocess the constraints into a natively-optimized Boolean Decision Diagram. The ordering is extracted from the output of the SAT solver using ArrayReshape and the AdjacencyLists property of SparseArray.

AllTotalOrderings[g_Graph]:=With[{canonicalgraph=CanonicalGraph[g],v=VertexCount[g]},
With[{orderedvlist=Values[KeySort[First[FindGraphIsomorphism[canonicalgraph,g]]]]},
Map[orderedvlist[[Flatten[ArrayReshape[SparseArray[Boole[#]],{v,v}]["AdjacencyLists"]]]]&,
Block[{Pos},SatisfiabilityInstances[
BooleanConvert[And[
And@@Table[BooleanCountingFunction[{1},v]@@Table[Indexed[Pos,{i,j}],j,v}],{i,v}],
And@@Table[BooleanCountingFunction[{1},v]@@Table[Indexed[Pos,{j,i}],{j,v}],{i,v}],
And@@Function[{i,j},And@@Table[And[Implies[
Indexed[Pos,{r,i}],Nor@@Table[Indexed[Pos,{k,j}],{k,1,r}]]],{r,1,v}]
]@@@EdgeList[TransitiveReductionGraph[canonicalgraph]]],"BFF"]//Head,
All,Method->"BDD"]]]]];

and

linearExtensions[set_List, po_?MatrixQ]:=AllTotalOrderings[Graph[set,DirectedEdge@@@po]];

This give dramatically improved performance. For the example in question,

With[{vars = {a, b, c, d, e, f, g, h, i, j, k, l},
porderlist = {{a, c}, {b, c}, {f, g}, {g, e}, {d, a}, {h, i}, {i, d}, {g, h}, {g, i}, {h, j}, {h, k}, {k, j}, {k, c}}},
Timing[Length[linearExtensions[vars, porderlist]]]]

(*Out[] = {1.53125, 12840}*)